Genome-wide SNPs detect fine-scale genetic structure in threatened populations of squirrel glider Petaurus norfolcensis

Australian arboreal mammals are experiencing significant population declines, particularly due to land clearing and resulting habitat fragmentation. The squirrel glider, Petaurus norfolcensis, is a threatened species in New South Wales, with a stronghold population in the Lake Macquarie Local Government Area (LGA) where fragmentation due to urbanization is an ongoing problem for the species conservation. Here we report on the use of squirrel glider mitochondrial (385 bp cytochrome b gene, 70 individuals) and nuclear DNA (6,834 SNPs, 87 individuals) markers to assess their population genetic structure and connectivity across 14 locations sampled in the Lake Macquarie LGA. The mitochondrial DNA sequences detected evidence of a historical genetic bottleneck, while the genome-wide SNPs detected significant population structure in the Lake Macquarie squirrel glider populations at scales as fine as one kilometer. There was no evidence of inbreeding within patches, however there were clear effects of habitat fragmentation and biogeographical barriers on gene flow. A least cost path analysis identified thin linear corridors that have high priority for conservation. These areas should be protected to avoid further isolation of squirrel glider populations and the loss of genetic diversity through genetic drift.

2016), and researchers have yet to take advantage of the power of next-generation sequencing technology and single nucleotide polymorphisms (SNPs).
SNPs are the most common form of sequence variation, characterised by single base nucleotide changes in the genome that differ among individuals (Brumfield et al. 2003). Advances in next-generation sequencing means that tens of thousands of bi-allelic, co-dominant SNP markers can be detected in the genome (Kumar et al. 2012). SNPs have wider genome coverage than microsatellites such that population samples of four individuals are enough to identify population structure (Liu et al. 2005;Shi et al. 2010;Lah et al. 2016;Dussex et al. 2018). Genome-wide SNP discovery with next-generation sequencing is gaining importance in conservation genetic research on threatened species given that small sample sizes are common. Additionally, genome-wide SNPs hold an advantage over microsatellites when it comes to inferring demographic history and genetic diversity estimates for conservation management purposes (Zimmerman et al. 2020). Genome-wide SNPs have evaluated the success of threatened Australian marsupial translocation programs in recent years, including that of the greater bilby, Macrotis lagotis and western barred bandicoot, Perameles bougainville (White et al. 2018). Furthermore, SNPs have effectively informed conservation management of threatened Australian marsupials including the koala, Phascolarctos cinereus (Kjeldsen et al. 2016) and the Tasmanian devil, Sarcophilus harrisii (Wright et al. 2019). SNPs should therefore be extremely valuable in assessing the effects of habitat fragmentation on squirrel glider populations at a fine-scale.
Lake Macquarie Local Government Area (LGA) is a nationally significant area within the distribution of squirrel gliders since it is believed to contain the highest density of the species across their range in the eastern seaboard of Australia (estimated 5000 individuals in the Lake Macquarie -Wyong area, concluded from the Wildlife Atlas records in New South Wales) (Smith 2002;Fallding 2015). Because of this, it is essential to monitor squirrel gliders in the area to ensure habitat fragmentation and biogeographical barriers do not negatively impact upon the conservation status further. Here we use mitochondrial DNA and genome-wide SNP markers to examine the demographic history and effect of habitat fragmentation on the population structure of threatened squirrel gliders in the Lake Macquarie LGA. We hypothesised that: (1) squirrel gliders in isolated patches would have low levels of genetic diversity and limited gene flow; (2) squirrel glider populations isolated by biogeographical barriers (Lake Macquarie and the Pacific Ocean) would be highly differentiated from other populations; and (3) landscape heterogeneity (vegetation and land use) would impact the genetic structure and differentiation of squirrel canopy that exceed their maximum glide distance of 70 m or by short trees that do not support a glide angle of 28.5 degrees (Van der Ree, 2002;Van Der Ree et al., 2004;Goldingay and Taylor, 2009). Additionally, squirrel gliders generally avoid remnant trees in the urban matrix as they are sensitive to light pollution, road density and associated noise (Francis et al., 2015). Limited access to food, as a result of land clearing and fragmentation, reduces the effective population size of squirrel gliders (Sharpe 2004). Evidence suggests that small isolated patches must be linked to a larger metapopulation to prevent inbreeding depression and localised extinction (Goldingay and Sharpe, 2004;Van Der Ree et al., 2004). Many studies have relied on radio-tracking to infer connectivity of patches and dispersal of squirrel gliders (van der Ree and Bennett, 2003;Sharpe and Goldingay, 2007;Goldingay et al., 2010;Brearley et al., 2011;Soanes et al., 2015), however few have examined genetic structure and gene flow to make this inference.
The emergence of landscape genetics (Manel et al. 2003) contributed substantially to the field of ecology, evolution, and species conservation in recent decades. Among others, landscape genetics combines population genetics, landscape ecology and spatial analytical techniques to quantify the effects of landscape composition, configuration, and matrix quality on gene flow using neutral genetic data (Balkenhol et al. 2015). The field has rapidly moved from the more traditional isolation-by-distance analyses to explain differences in the genetic structure of populations that drill down further into movement and dispersion dynamics. Estimates of gene flow can be tested against measures of effective geographic distance to find estimates of landscape resistance that best fit with the genetic data. Landscape genetics has proven useful to quantify the effects of habitat fragmentation on genetic structure, and identify landscape features that are important for maintaining genetic connectivity between populations within heterogeneous landscapes (e.g. Sharma et al. 2013), as well as pinpoint barriers to gene flow (Cushman and Lewis 2010).
Large-scale studies have compared the genetic diversity of squirrel gliders across the entire east coast of Australia (Pavlova et al. 2010;Taylor et al. 2011), but key squirrel glider hotspots remain unsampled. Additionally and to our knowledge, only one published study has examined the effect of fragmentation on a fine-scale, in Mackay and Brisbane (Queensland) (Goldingay et al. 2013). Assessing the fine-scale genetic structure of squirrel glider populations would assist with the development and use of local conservation management strategies in areas where they are threatened. Thus far, squirrel glider research has only used mitochondrial (cytochrome b) and microsatellite DNA markers to examine genetic structure (Pavlova et al. 2010;Taylor et al. 2011;Goldingay et al. 2013;Dudaniec et al. of squirrel gliders in the Lake Macquarie area associated with an understory of Xanthorrhoea and Banksia species, and scribbly gum (Eucalyptus haemastoma or racemosa), smooth-barked apple (Angophora costata) and red bloodwood trees (Corymbia gummifera). Once potential sites were located, we referred to BioNet records (NSW Office of Environment and Heritage 2021), and used spotlighting surveys and/or infrared cameras (Gracanin et al. 2019) to confirm glider presence at a site prior to live trapping.

Sample collection and DNA extraction
We conducted live trapping at 32 locations over a fouryear period, with squirrel gliders captured at 14 locations ( Fig. 1). Each site was subject to at least one week of live trapping using 12 traps per site (two rows of six with each trap spaced 50 m apart), and further focus was given to sites that captured squirrel gliders during this period to increase sample size. Once squirrel gliders were detected at a site, we re-trapped the area every six months to a year across the survey period. We used a combination of Mawbey traps, Elliot B traps, cage traps and Winning and King pipe traps to live trap gliders (Mawbey 1989;Quin 1995;Winning and King 2008). All traps were secured to trees so that the trap entrance was at least two meters from ground height. Each trap contained a bait ball made from a mixture of peanut butter, honey, and oats (Winning and King 2008), leaves for insulation, and we sprayed a 1:4 ratio of honey:water around each trap entrance and up the tree to a height of 6 m as a further attractant (Sharpe and Goldingay 2007).
Traps were checked each morning at sunrise. We weighed, sexed, measured (tail length, right hind foot length, head width and head length) and marked each squirrel glider with a metal ear tag or an ear punch combination. Before releasing the squirrel gliders, we collected DNA samples from 94 squirrel gliders in the form of ear tissue (or in the early stages of the project, a buccal swab, n = 30), using a 2 mm metal ear punch (Able Scientific, Australia) that was previously sprayed with 70% ethanol and flamed for sterilization. Once cooled, a 2 mm tissue punch was taken from the outside edge of the ear and stored in sterilized vials containing 80-95% ethanol. DNA was kept at -20 °C prior to DNA extraction. The processing of each squirrel glider was limited to a maximum time of five-minutes, and individuals were released on the same tree in which they were caught.
We extracted genomic DNA from the buccal swab and ear tissue samples using the One-4-ALL Genomic Miniprep Kit (BioBasic Inc. Ontario, Canada) and DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany) as per the manufacturer's instructions. glider populations. This is the first study to report on the genetic structure of the Lake Macquarie LGA squirrel glider metapopulation, despite the LGA being a stronghold population for the species (Smith 2002;Fallding 2015).

Study area
Our study was undertaken in the Lake Macquarie LGA, ~ 130 km north of Sydney, Australia (Fig. 1). Lake Macquarie LGA has experienced a high level of land clearing resulting in fragmentation of forests that were recognised as critical squirrel glider habitat (Smith 2002;Fallding 2015). In addition to habitat fragmentation, squirrel gliders also face added pressure from biogeographical barriers in the landscape, namely a large natural saltwater lake (Lake Macquarie, 120 km 2 ) that presents a barrier to gene flow in the LGA as well as the Pacific Ocean on the east coast ( Fig. 1).
We selected squirrel glider trapping locations in the Lake Macquarie LGA through a combination of spatial modelling and visual identification of suitable vegetation. Sites required dens (large hollow bearing trees) and eucalypts (Corymbia maculata, Eucalyptus robusta, Eucalyptus tereticornis), acacias (Acacia irrorata) and banksias (Banksia spinulosa; Banksia integrifolia) since squirrel gliders are known to feed on their sap, gum and nectar (Ball et al. 2009). Additionally, Smith and Murray (2003) found higher densities The 30-meter-wide powerline easement is located at WR sequence of individuals. Haplotype diversity and nucleotide diversity were calculated for the Lake Macquarie and Forster regions. Haplotype diversity (h) gave the probability that two randomly sampled alleles would be different (probability 0-1), while nucleotide diversity (π) determined whether the differences between haplotypes were minor (approaching zero) or major (approaching 1). Tajima's D and Fu's Fs tests were used to test the neutrality of the mitochondrial DNA sequences within populations and regions (Tajima 1989;Fu 1997), and a mismatch distribution analysis was conducted to determine whether squirrel gliders had undergone a recent range expansion (Schneider and Excoffier 1999). All of these metrics were calculated in DnaSP v 5.10 (Librado and Rozas 2009). A median-joining network was constructed in PopART to visualise relationships among haplotypes (Bandelt et al. 1999).
The mutation rate of the cytochrome b gene was calculated to gain insight into the demographic history of squirrel gliders. To do this we trimmed 70 cytochrome b sequences and aligned them to closely related species. This included a 374 bp cytochrome b sequence from a sugar glider (Petaurus breviceps, GenBank Accession Number FJ657661.1), a mahogany glider (Petaurus gracilis, Gen-Bank Accession Number HQ284045.1) and a leadbeaters possum (Gymnobelideus leadbeateri, GenBank Accession Number NC_053898.1). The alignment was imported to BEAUti v2.6.6 to create a file for BEAST analysis (Bouckaert et al. 2019). We selected a General Time Reversible

Mitochondrial DNA: cytochrome b gene
We used genomic DNA samples extracted from buccal swabs and ear tissue for mitochondrial DNA analysis. We amplified a 385 base pair fragment of the mtDNA cytochrome b gene using primer pair L14724 (5'-CGA AGC TTG ATA TGA AAA ACC ATC GTT G-3') and H15149 (5'-AAA CTG CAG CCC CTC AGA ATG ATA TTT GTC CTC A-3') as described by Kocher et al. (1989) and Irwin et al. (1991). PCR products were treated in the post-PCR lab via the ExoSAP-IT™ enzymatic method (Affymetrix Inc. USB Corp., Cleveland, Ohio) prior to sequencing with an ABI PRISM BigDyeTM Terminator Cycle Sequencing Kit (Applied Biosystems, USA). We cleaned the sequencing products with a precipitation cleanup method and sequences were generated with a 3130xl Genetic Analyzer (Applied Biosystems, Australia), at the University of Wollongong, School of Biological Sciences. We visually inspected and manually edited the sequences using the program Chro-masPro v 1.33 (Technelysium Pty Ltd, Brisbane) and aligned them using BioEdit v 7.2.6.1 (Hall 1999). Cytochrome b sequences were obtained from 55 squirrel glider samples in the Lake Macquarie LGA, and from 15 buccal swab samples that were donated to the project from the Forster region, approximately 130 km northeast of Lake Macquarie. Therefore, cytochrome b sequences were generated for a total of 70 individual squirrel gliders (Table 1).
All mtDNA haplotypes were identified as a unique arrangement of nucleotides in the cytochrome b DNA Table 1 Cytochrome b data from DNA extracted from buccal swabs and ear tissue (385 bp): squirrel glider sample size (n), haplotypes, haplotype diversity (h), nucleotide diversity (π) and Tajima's D and Fu's F values for locations in Lake Macquarie and Forster*. Refer to Fig. 1 Fig. 1. Only buccal swabs were collected from these locations, so they were not included in the nuclear DNA analyses. General locations are as follows: 151.529,152.528,152.542,152.480 sequencing to detect a large number of SNPs at a relatively low cost (Kilian et al. 2012). Additionally, DArTseq optimises the technology for each species by selecting the most appropriate complexity reduction method. For this squirrel glider study, the PstI-SphI method was used prior to PCR amplification and sequencing of fragments on an Illumina Hiseq2500. This PstI-SphI restriction enzyme combination has proved extremely efficient in detecting SNPs in Australian marsupials (Schultz et al. 2018;Kjeldsen et al. 2019;Wright et al. 2019), including sugar gliders (Petaurus breviceps) (Knipler et al. 2021).
DArTseq aligned the squirrel glider sequences to the Leadbeaters possum (Gymnobelideus leadbeateri) reference genome and retained SNPs with a minimum sequence identity of 70%. Calling quality was assured with a high average read depth per locus (over 20 reads per locus average across all markers). Following this, we filtered SNPs and individuals using dartR 1.1.11 (Gruber et al. 2018) in R 4.0.2 (R Core Team 2015). Duplicate samples of lower quality were removed from the dataset so that only unique individuals remained. Loci were removed if the call rate was < 0.95, if they had reproducibility < 0.99, and if the hamming distance was < 0.2. Loci were also removed if they deviated from Hardy-Weinberg Equilibrium (p < 0.05) and the SNP data set was analysed to ensure no monomorphic loci remained. The resulting dataset was used for inferring demographic history as it contained linked loci (8,580 SNPs retained). We used TASSEL v5.2.78 (Bradbury et al. 2007) to test for neutrality in the form of Tajima's D (± SD) using a step size of 100 and a window size of 500 for populations with ≥ 4 samples. Current effective population sizes (Ne) were calculated for each population using the larger partially filtered dataset of SNPs and NeESTIMATOR v2.1 (Do et al. 2014). A critical threshold value of 0.05 was selected and random mating was assumed. After this, we further filtered the dataset to remove linked loci using the "gl.filter.secondaries" function in the dartR package. We conducted an OutFLANK outlier test in R to infer the neutrality of SNPs (Whitlock and Lotterhos 2015). No loci were identified as being under selection, so all were retained, and this filtered dataset was used for further population genomic analyses (6,834 SNPs).
We undertook all subsequent analyses in R 4.0.2 (R Core Team 2015) unless otherwise stated. We calculated observed and expected heterozygosity for each location to gain insight into the genetic diversity of squirrel glider populations. Observed heterozygosity (H obs ) was calculated with the R package dartR, while expected heterozygosity (H ex ) was calculated with the package adegenet 2.1.3 (Jombart 2008). We also calculated the inbreeding coefficient for each squirrel glider location (F IS = 1-H obs /H ex ), and the overall average F IS for the total region. The nucleotide (maximum-likelihood GTR) site model, with a strict molecular clock and calibrated yule model for priors.
In priors, sugar gliders and squirrel gliders were combined as most recent common ancestors, and monophyly was enforced. To calibrate the tree, we selected a Normal distribution, and the mean was centered at 1.73 million years with a standard deviation of 0.5 million years (30% of the divergence estimate). This mean was taken from Meredith et al. (2010) after examining genetic research and fossil information provided by Malekian et al. (2010) and Meredith et al. (2010). Finally, all taxa were combined as most recent common ancestors and monophyly was enforced. The mean was centered at 15.97 million years before present as that's when the Leadbeaters possum is estimated to have diverged from gliders (Meredith et al. 2010). A standard deviation of 5 million years was selected (30% of the divergence estimate). The parameters for the molecular clock rate (clockRate) and the speciation rate (birthRateY) of the Yule tree prior were entered as Alpha = 0.001 and Beta = 1000. Markov chain Monte Carlo runs were set to 10,000,000 generations. The output file was generated and run through BEAST v2.6.6 (Bouckaert et al. 2019). The resulting tree was visualized in FigTree v1.4.4 (Rambaut 2019) and the log file was analyzed in Tracer v1.7.2 to estimate the mutation rate of the cytochrome b gene in squirrel gliders (Rambaut et al. 2018). We used this mutation rate to estimate the average number of generations taken for a single mutation to be acquired, assuming a generation time of 4.5 years as per Pavlova et al. (2010).
Finally, we used the mutation rate from the previous analysis to estimate the female effective population size of squirrel gliders in the Lake Macquarie LGA. We displayed the female effective population size with a Bayesian skyline plot in BEAST as per Drummond et al. (2005). Changes to the female effective population size were examined through time.

Nuclear DNA: SNPs
We examined squirrel glider genomic DNA samples on 1% agarose gel electrophoresis and with a NanoDrop 2000 Spectrophotometer (Thermo Scientific) to identify samples with appropriate concentration and purity for next-generation sequencing. The buccal swab genomic DNA samples did not pass this quality control stage. Therefore, we selected 94 tissue DNA samples (belonging to 89 squirrel glider individuals) for sequencing.
We plated and sent genomic DNA to Diversity Arrays Technology (DArTseq), University of Canberra, Australia, for concentration and sequencing (Kilian et al. 2012). DArTseq is a high throughput method that combines genome complexity reduction methods and next-generation distances between pairwise populations and compared them to the genetic distance matrix of populations in the form of F ST /(1-F ST ), to see if there was a significant relationship between geographic distance and genetic distance of populations. The pairwise F ST values from dartR were used in this case. We used the results of this IBD analysis as a baseline, to run least cost path analyses and determine whether they could explain more of the genetic variation in the landscape.
We ran two separate least cost path (LCP) analyses using the "gl.genleastcost" function in the dartR package. LCP approach estimates the shortest distance between target core areas while accounting for resistance to movement (inferred here from the genetic data) (Adriaensen et al. 2003). LCP requires an input resistance surface, also called a friction matrix, and that is a spatial layer reflecting the degree to which a location in the landscape facilitates or impedes movement of a focal species (Zeller et al. 2012).
We created two friction matrices in ARCMAP 10.7.1 (ESRI, California) for each LCP scenario. The first friction matrix only considered water bodies as biogeographical barriers. Hence, the Lake Macquarie and the Pacific Ocean were manually assigned "NoData" cell values, while all terrestrial areas were given a cell value of 1, assuming constant movement capacity across a homogeneous terrestrial landscape within our study area. This analysis will hereafter be referred to as an isolation by barrier analysis.
The second friction matrix used the Plant Community Type (PCT) spatial layers from Bell et al. (2016) and Eco Logical Australia Pty Ltd (2003), and land use spatial layers from the NSW Government (2007). We assigned values to the cells (pixels in the map) based on our opinion, the literature, and a range of test models we ran. Values scaled from 1 for optimal habitat with low resistance that eased movement, to 200 for unsuitable habitat with high resistance that would render movement (and hence dispersal) challenging. We divided vegetation layers into five categories based on the diet and habitat requirements of squirrel gliders (Smith 2002;Smith and Murray 2003;Ball et al. 2009; Niche Environment and Heritage 2013; Payne 2016): highly suitable vegetation (1), suitable vegetation (10), moderately suitable vegetation (20), moderately unsuitable vegetation (75) and unsuitable vegetation (100). Appendix 1 details how the 62 PCTs were assigned. Land use was divided into the following classes: Urban residential (80), roads (100), railways (80) (because there is less traffic and the width is smaller than roads), rural residential without agriculture (70), rural residential with agriculture (100), native/exotic pasture matrix (100), grazing irrigated modified pastures (100), irrigated turf farming (100), mines/quarries and large cleared areas (200), the Lake and Ocean (NoData, aka complete barrier to dispersal). Highly suitable vegetation was given a 35-meter buffer so that two patches either side of a barrier diversity of SNPs (π ± SD) was calculated for each population in TASSEL as a measure of genetic variation.
We examined structure and population differentiation among squirrel gliders in Lake Macquarie LGA using Analyses of Molecular Variance (AMOVA's), and we plotted a Principal Coordinate Analysis (PCoA) to assist with the visualisation of this genetic structure. The PCoA was conducted in dartR. We ran the AMOVA's with the R package poppr 2.9.0 using 9,999 permutations (Zhian Kamvar 2021). The first AMOVA examined the genetic differentiation of squirrel gliders in two regions: west (locations WYB, WYC, WR, SP, OR) and east (locations RA, CROB, GMP, FERNS, OCN, FERNN, DUD, RHL, GSCA) of Lake Macquarie. The second AMOVA examined the genetic differentiation of squirrel gliders across all 14 populations, and the third served as a control with two haphazardly drawn lines dividing the sampling locations into three groups (group 1: WYB, WYC, SP, group 2: WR, OR, CROB, RA, FERNS, OCN, and group 3: GMP, FERNN, DUD, RHL, GSCA). Next, we examined genetic differentiation in the form of pairwise F ST between the 14 locations of squirrel gliders. Pairwise F ST values were generated with the package dartR and 10,000 bootstraps were performed across loci to generate the corresponding p-values. Similarly, the average F ST was calculated for the overall region.
We further examined squirrel glider genetic structure with STRUCTURE 2.3.4 (Pritchard et al. 2000). STRUC-TURE uses a Bayesian clustering approach to assign individuals to groups ("K") where members display similar patterns of genetic variation. These groups are indicative of ancestral populations that individuals may have originated from (Pritchard et al. 2000;Pritchard et al. 2003). We ran the analysis eight times for each K value (1-15) with a 10,000-length burn-in period and 10,000 Monte Carlo Markov Chain replications. The most likely K values were selected after viewing the results in Structure Harvester (Web v0.6.94) and referring to the Evanno method and the largest Delta K value (Evanno et al. 2005;Earl and Von-Holdt 2012).

Landscape genomic analyses
In addition to identifying genetic differentiation and genetic structure, we wanted to identify landscape features that influenced (diverting, preventing, or enhancing) squirrel glider gene flow in areas of the Lake Macquarie LGA. Genomic data was combined with spatially-explicit data to investigate the effect of geographic distance and habitat heterogeneity on genetic distance of squirrel glider populations. First, we undertook an isolation-by-distance (IBD) analysis in the dartR package using 9,999 permutations. This mantel test took the log of the straight-line Euclidean 57.14% of individuals and most of the locations (Table 1). In comparison, haplotype two (H02) was only observed once and it was unique to Lake Macquarie's MP location. Haplotype seven (H07) was also unique to the MP location at Lake Macquarie. Locations SP and OR had a single haplotype at each site (Table 1). Lake Macquarie Tajima's D and Fu's F values indicated population expansion and an excess of rare mutations in the Lake Macquarie squirrel glider populations; however, these values were not significant (-1.07922 and − 2.254 respectively, p > 0.10) (Table 1). Similarly, the results of the mismatch distribution analysis produced a unimodal curve. The overall average number of pairwise differences was k = 0.951 and the raggedness statistic was r = 0.0588. The mutation rate for the squirrel glider cytochrome b gene was calculated to be 1.63% per site per million years. The estimated time for a single mutation was 164,036 years (lower 95% HD = 312,407, upper 95% HD = 102,053), or 2.74329E-05 mutations per generation ( Fig. 2A), and a star-like haplotype network was produced (Fig. 2B). The Bayesian skyline plot depicted a constant female effective population size over the last 200 years (Appendix 2).

Genetic differentiation and structure (SNPs)
DArTseq returned 22,579 loci from 92 squirrel glider samples and 14 locations in the Lake Macquarie LGA. Two squirrel glider samples did not meet quality control standards and were excluded. During the final filtering stage in R, five duplicate squirrel glider samples were removed. Two datasets were generated after stringent filtering: one for inferring demographic history as it contained linked loci (8,580 loci, 87 squirrel glider individuals, 14 locations in the Lake Macquarie LGA), and one for subsequent population (e.g. road/river) would intersect within maximum gliding distance (70 m (Van der Ree and Bennett 2003)), indicating crossing potential. If cells overlapped between vegetation and land use, then the minimum value was selected for the final friction matrix. This analysis will be heron referred to as a least cost path analysis.
We computed new geographic distances using the two friction matrixes (isolation by barrier raster and the least cost path raster, cell factor = 15, function = mean in Arc-GIS) and the spatial coordinates of individuals. These new pairwise geographic distances were then tested for correlation with genetic distances and number of neighbours = 8 with the "gl.genleastcost" function of the R package dartR. We then compared the statistics from each spatial analysis (IBD, isolation by barrier analysis and least cost path analysis) to find the best fit (Milanesi et al. 2016).

Trapping
We trapped squirrel gliders in 14 out of the 32 sites surveyed. We collected ear tissue samples from 94 squirrel glider individuals, while an additional 30 squirrel gliders had buccal swab samples collected. We caught many male and female squirrel gliders over consecutive years in the original location they were trapped, suggesting site fidelity. For example, we recaught a male and female squirrel glider after two years within 50 m of their original detection at location RA, while a male squirrel glider was recaught after two years within 100 m of its original detection at location GSCA. Additionally, we caught one adult male and one adult female squirrel glider on either side of a 30-meter powerline easement during one week of live trapping (location WR).

Cytochrome b
Of the 70 squirrel gliders sequenced, nine polymorphic sites were identified (2.34%), 376 conserved sites (97.66%) and eight haplotypes. The average total nucleotide diversity was low for both Lake Macquarie LGA and Forster (Lake Macquarie π = 0.002 ± SD 0.000, Forster π = 0.002 ± SD 0.003) and total haplotype diversity was moderately high for Lake Macquarie LGA (h = 0.648 ± SD 0.052) but lower in Forster (h = 0.448, ± SD 0.134). A BLASTn search of GenBank failed to find a match for haplotypes H02 (habitat fragment MP, refer Table 1), H04 (locations RHL, GSCA, OCN, SP, WYB) and H07 (habitat fragment MP, refer Table 1), therefore these haplotypes are believed to be unique to this study. Haplotype three (H03) was the most common, occurring in Fig. 2 (A) Maximum-likelihood tree reflecting divergence times of cytochrome b haplotypes in leadbeaters possum (G. leadbeateri, n = 1), mahogany glider (P. gracilis, n = 1), sugar glider (P. breviceps, n = 1) and squirrel glider (P. norfolcensis, n = 70). Timeline in millions of years before present (MYA). (B) MtDNA median-joining network displaying the relationship among eight cytochrome b haplotypes (H01-H08) of squirrel gliders in Lake Macquarie and Forster, NSW (n = 70). A slash in the branches linking two haplotypes indicate one nucleotide substitution and two slashes indicate two nucleotide substitutions. Refer to Table 1 to see where each haplotype was found were that the southwestern squirrel glider populations (OR, SP, WR, WYB, WYC) separate out from the northeastern populations on axis 1, while FERNN and OCN each appear to separate out on axis 2. Axis 4 depicts a clear separation of population OR.
We compared genetic differentiation of squirrel glider populations through pairwise F ST . Statistically significant pairwise F ST values ranged from 0.014 (DUD vs. RHL) to 0.335 (CROB vs. SP) ( Table 4). F ST values greater than 0.15 are indicative of significant differentiation and genetic structure (Frankham et al. 2002), and in this case many of the pairwise comparisons produced results that indicated considerable differentiation ( Table 4). Squirrel gliders at location OR were greatly differentiated from all other areas (F ST 0.149 to 0.273), and CROB, FERNN, OCN and SP displayed substantial genetic differentiation from most other areas (Table 4). As a result, the average F ST was quite high with a value of 0.148. STRUCTURE analyses provided further insight into the genetic structure of squirrel gliders in the Lake Macquarie LGA. The Structure Harvester results showed that Delta K had a peak at K = 2 (Delta K = 47.57), with a small peak again at K = 6 (Delta K = 5.22) (Appendix 3). Results for both are presented here since there can be bias towards finding K = 2 with the Delta K method (Cullingham et al. 2020). Upon visualising the K = 2 STRUCTURE admixture plot, it was clear that the squirrel glider individuals were divided into two clusters, similar to the results on axis 1 of the PCoA plot (populations northeast of the lake and populations southwest of the lake) (Fig. 5). When K = 6, unique clusters appeared within the FERNN, OCN, OR, and WR populations, suggesting that the Lake Macquarie individuals could have been derived from six squirrel glider ancestral sources (Fig. 5). The separation of these populations (FERNN, OCN, OR and WR) reflects the patterns observed in the PCoA plots (Figs. 4 and 5). Admixture was present in all sampling locations; however, some individuals were pure to clusters within FERNN, GSCA, OCN, OR and WR (q value > 0.8).

Spatial analyses
With the isolation by distance analysis, we found a significant relationship between geographic straight-line distances and genetic distances of squirrel glider populations (R 2 = 0.3227, p < 0.01) (Table 5). Second, an isolation by barrier analysis was performed to see if the biogeographic barriers explained more of the genetic variation (Lake Macquarie and the Pacific Ocean). The R 2 value did improve (R 2 = 0.5390, p < 0.001) ( Table 5), but it did not consider habitat heterogeneity. Thirdly, we ran a least cost path analysis to include landscape features and habitat heterogeneity genetic analyses (6,834 loci, 87 squirrel glider individuals, 14 locations in the Lake Macquarie LGA).
Tajima's D was calculated for populations with ≥ 4 samples using the partially filtered dataset that contained linked SNPs. All values were negative, which is indicative of population expansion after a recent genetic bottleneck (Fig. 3). Care should be given when interpreting these results due to small sample sizes (n = 4-16) (Subramanian 2016). Effective population sizes of squirrel gliders in Lake Macquarie populations were very small (Table 2) and should be viewed with caution due to the small sample sizes. These values should be used as a comparison between populations to identify regions of concern, rather than an accurate estimate of effective population size (Zachos et al. 2016).
Using the strictly filtered SNP dataset of 6,834 loci, we compared observed heterozygosity to expected heterozygosity at each location, to understand the levels of squirrel glider genetic diversity. Observed heterozygosity was higher than expected heterozygosity at every location except for location WR, however the values were low indicating low levels of genetic variability (Table 2). There was no evidence of inbreeding (overall FIS = 0.015). The nucleotide diversity of the SNPs displayed similar patterns to the mtDNA results, with population OR, GMP and OCN having low nucleotide diversity. Three AMOVA's examined the genetic variation of squirrel glider populations. Although west and east Lake Macquarie explained 7.79% of the genetic variation, this value resembled the variation observed between the random control groups (5.51%) ( Table 3). When the 14 sampling locations (i.e. putative populations) were included in an AMOVA, further genetic variation was explained (15.12%) ( Table 3) (Fig. 4). The cumulative percentage of genetic variation explained by the first four axes was 23.2%, with axis 1 accounting for 8.6% of the variation (Fig. 4). The most important observations

Table 2
Location information and number of genotyped squirrel gliders (N) that were sampled in the Lake Macquarie area. Mean effective population sizes (Ne) were calculated for locations containing ≥ 4 samples, using the partially filtered dataset containing linked loci (8,580 loci) and NeESTIMATOR (0.05 lowest allele frequency). Nucleotide diversity π, mean observed heterozy-  from a small effective population size (Grant and Bowen 1998;Avise 2000). This is reflected in the star-like haplotype network, which is consistent with what is displayed after a historical genetic bottleneck with many haplotypes branching from the dominant central H03 haplotype (high haplotype diversity) depicting minimal nucleotide differences (low nucleotide diversity). The unimodal curve from the mismatch distribution and negative Tajima's D values for both the mtDNA sequences and the SNP dataset reinforce this concept (Tajima 1989). It was estimated to take 164,036 years for a mutation to occur in the cytochrome b gene. This information paired with the star-like haplotype network indicates that a range expansion occurred in the mid-late Pleistocene, which supports the results of Pavlova et al. (2010). We also used SNPs to examine whether the genetic diversity of squirrel gliders would be low in habitat patches within Lake Macquarie LGA. Low genetic diversity has been observed in arboreal mammal populations within fragmented landscapes such as the Siberian flying squirrel (Pteromys volans) (Lampila et al. 2009) and ringtail possum (Pseudocheirus peregrinus) (Lancaster et al. 2016). Despite these predictions, the genetic diversity values did not differ greatly between the mtDNA and SNP data, and squirrel gliders displayed higher than expected levels of observed heterozygosity within the Lake Macquarie LGA, with no evidence of inbreeding (FIS = 0.015). Higher than expected heterozygosity could potentially be explained by the genetic bottleneck detected with the mtDNA analyses, as bottlenecks lead to a quick loss in rare alleles while heterozygosity decreases slowly over time (Lampila et al. 2009). Additionally, excess heterozygosity in current squirrel glider (biogeographical barriers, vegetation, and land use). When run, it was clear that the least cost paths between populations (Fig. 6) explained the greatest proportion of genetic variation between populations (R 2 = 0.6467, p < 0.001) ( Table 5). Woinarski et al. (2014) lists habitat fragmentation as the main threat to squirrel gliders, which is particularly concerning in the Central Coast of NSW (Smith 2002;Fallding 2015). In this study, we examined whether the effects of habitat fragmentation were detectable on the populations genetic structure. We expected that squirrel gliders in isolated patches would have low levels of genetic diversity and limited gene flow. We also sought to identify areas that would enhance gene flow between populations in our study area and identify areas that present potential barriers to dispersal.

Discussion
Our first hypothesis was partly rejected as we detected high haplotype diversity in squirrel glider populations, despite previous research recording low haplotype diversity in squirrel glider populations 40 km from Lake Macquarie in Port Stephens (h = 0.18), and 20 km from Lake Macquarie in Warnervale (h = 0.50) (Pavlova et al. 2010). There did not appear to be haplotypes grouped with respect to any geographical barrier or feature. MtDNA has a relatively fast mutation rate and is useful for determining the genetic history of populations (Arif et al. 2011). High haplotype diversity and low nucleotide diversity in the Lake Macquarie squirrel glider population points to a historical genetic bottleneck followed by rapid expansion and recovery  Table 3 Analysis of molecular variance (AMOVA) results testing genetic structure in squirrel gliders within the Lake Macquarie LGA. Results are generated from 6,834 loci and 87 individuals. df = degrees of freedom levels of genetic structuring in the Lake Macquarie squirrel glider populations. The F ST values were unusually high considering the fine scale of this study, and they were considerably higher than sugar glider populations in the same study area (FST values ranging from 0.011 to 0.131 based on 11,292 SNPs) (Knipler et al. unpublished data). Despite the two species sharing similar life history traits and requiring similar habitat (Lindenmayer 2002), the difference was considerable (squirrel glider F ST values from 0.014 to 0.335). For example, we found a high level of genetic differentiation between squirrel gliders at locations OR and WR (F ST = 0.149), but in Knipler et al. (unpublished data) the F ST value for sugar gliders between OR and WR was significantly lower (F ST = 0.042). Squirrel gliders had greater pairwise F ST values for every location where squirrel gliders and sugar gliders were present together (OR, SP, WR, WYB, WYC) (Knipler et al., unpublished data). It is possible that squirrel gliders are more specialised and sensitive to the urban matrix and that is why they are experiencing higher levels of differentiation than sugar gliders. This warrants further investigation as population differentiation is characteristic of most specialist, arboreal mammals in fragmented habitats (Barratt et al. 1999;Larsen et al. 2013;Baden et al. 2014;Sgarlata et al. 2016).
While there was an isolation-by-distance pattern of gene flow in the squirrel glider populations in the Lake Macquarie LGA, it did not explain a large proportion of the genetic variation. Microsatellite analyses by Goldingay et al. (2013) for squirrel glider populations in the cities of Mackay and Brisbane, Queensland, found no isolation-by-distance effect, with strong genetic divergence occurring in squirrel glider populations 3 km apart after 30 years of landscape change. Goldingay et al. (2013) suggests that genetic drift in habitat fragments can have a stronger effect on allele frequencies than geographic distance. In the Lake Macquarie LGA, there was evidence of significant genetic differentiation between squirrel glider locations that were geographically close to each other. For example, we detected high levels of genetic differentiation between two locations that were only 1.8 km apart but divided by urban development in the town Whitebridge (FERNN and OCN F ST = 0.179), and locations 900 m apart but divided by urban development in the town Valentine (CROB and RA F ST = 0.131). This contrasts with squirrel glider samples that were collected 2.2 km and 1.5 km from each other but connected by forest (WYB and WYC F ST = 0.087, DUD and GSCA F ST = 0.027). This is where the landscape genetic analyses proved useful. By combining next-generation sequencing and landscape genetic analyses our study was able to explain a large portion of the genetic variation observed in the squirrel glider populations. To our knowledge, just one other study has used resistance models and least cost path populations could be the result of admixture and population growth following the historic bottleneck. This concept is supported by the mtDNA analyses and the nuclear F ST and STRUCTURE results. Taken together, the data suggest that the impacts of urban fragmentation of the Lake Macquarie forest system has not yet produced a substantial loss of heterozygosity with the associated risks of loss of fitness due to inbreeding effects (Frankham et al. 2002). Nevertheless, restriction of squirrel gliders to smaller and more isolated habitat fragments is likely to produce a drastic effect in coming generations. Therefore, we recommend continued monitoring of heterozygosity in isolated populations.
Genetic differentiation (i.e. structure) between fragments was expected to be quite high due to isolation by both the lake and extensive urbanisation. This expectation was based on the fact that squirrel gliders are sensitive to noise pollution in the urban matrix (Francis et al. 2015) and have a maximum glide distance of 70 m (Van der Ree and Bennett 2003). This theory was supported with evidence of high analyses to infer genetic structure of squirrel gliders in the past (Dudaniec et al. 2016), but the proportion of variation explained by their model was lower than ours (R 2 = 0.027 versus R 2 = 0.647). A possible reason for this improvement is the incorporation of genome-wide SNPs rather than microsatellites, and the use of a 35-meter buffer to account for the maximum glide distance of 70 m (Van der Ree and Bennett 2003).
Following the least cost path analysis, we identified thin vegetation strips as possible corridors to enhance squirrel glider gene flow in the following locations: Fishery Point Road (Bonnells Bay), Kilaben Bay (Rathmines), Fernleigh Track (Whitebridge), vegetation west of Fernleigh Tack (Charlestown), Pacific Highway crossing (Jewells to Bennetts Green), vegetation patches between Glenrock State Conservation Area (Highfields) and Blackbutt Nature Reserve (New Lambton), and vegetation along creek lines (Fig. 7). Locations OR and SP are situated in heavily urbanised areas on isolated peninsulas within the lake, so they were also identified as areas of conservation concern.
Our study not only identified corridors, but it also investigated possible barriers to dispersal. A 30-meter-wide powerline easement and the corresponding fragmentation of habitat did not seem to impact squirrel glider movement in location WR. Empirical movement data from GPS collars also confirmed that squirrel gliders can cross two lane roads in our study area (Meyer et al. 2021). However, Holt Table 4 Pairwise  (Knipler et al., unpublished data), and since squirrel gliders have experienced higher levels of genetic structure in the Lake Macquarie LGA then they may also be experiencing genetic and Brady (2011) have reported the reluctance of female Mahogany gliders (Petaurus gracilis) to cross 30-meter powerline easements, so further research is needed to confirm that the crossing observed by a female squirrel glider in this study is not an anomaly. Additionally, the power poles at WR were wooden and attention should be given to areas that use metal power poles, as the squirrel gliders in this study may have been using the wooden poles as steppingstones across the landscape.  Table 5 Models used to test the relationship between geographic distances and genetic distances between squirrel glider populations in the Lake Macquarie LGA. The distances between populations generated by the least cost path model explained the greatest proportion of the genetic variation between populations compared to the two other models Fig. 6 Friction matrix and least cost paths between 14 squirrel glider populations in Lake Macquarie. Raster was generated in ARCMAP and the least cost paths were created with the R package dartR. Cost raster cell values ranged from optimal habitat (1) to unsuitable habitat (200) based on vegetation communities and land use type. The lake and ocean were given "NoData" values and assumed a complete barrier to dispersal, except where the vegetation buffer allowed crossing using maximum gliding distance of 70 m Finally, considering the difficulty in collecting large sample sizes of threatened species, the use of genome-wide SNPs as we have undertaken, has allowed a thorough investigation into the fine-scale structure of squirrel gliders in the Lake Macquarie LGA. This data is vital and will prove valuable for conservation management decisions in the local area.

Supplementary Information
The online version contains supplementary material available at https://doi.org/10.1007/s10592-022-01435-9. differentiation on either side of this road, especially because there are no mature trees in the median strip. Structures such as glider poles (Ball and Goldingay 2008;Goldingay et al. 2019) should therefore be installed to assist the gene flow of sugar gliders, which would also help the squirrel gliders that inhabit the same area.
With our filtered SNP data, we identified corridors that could potentially enhance gene flow of squirrel gliders. Landscape genomics is a rapidly growing field that presents tremendous potential to estimate functional connectivity (Balkenhol et al. 2017). Genetically derived connectivity estimates reflect past landscape permeability due to the time it takes to detect barriers on species gene flow (Cushman et al. 2013), and thus does not necessarily detect current gene flow in a rapidly evolving landscape such as Lake Macquarie LGA. Furthermore, genetic data do not directly convey how animals move through the landscape. Hence, for conservation purposes, it is often recommended to complement connectivity analyses with other sets of empirical data. Movement data in particular would enable us to examine how animals respond to landscape features as they disperse (Cushman 2013), and would provide a more complete picture of the resources needed by squirrel gliders as they cope with habitat change.
We here used LCP to estimate functional connectivity, which assumes animals have perfect knowledge of the landscape, and therefore move across the landscape in an optimal way. In contrast, circuit-theory connectivity is based on random walk and applies concepts related to flow of charge through an electrical circuit to the movement of individuals through the landscape. The output of such models is the probability of movement across each pixel of the landscape and incorporates all possible pathways across the area studied (McRae et al. 2008). As squirrel gliders are unlikely to have a complete knowledge of their surroundings when dispersing, future studies using circuit-theory modelling may be useful to confirm our findings.
Continued fragmentation is the greatest risk to squirrel gliders in our study area, particularly as the effects of urbanisation are exacerbated by the biogeographical barrier of the Lake Macquarie. Conservation of remaining forest habitat should be a priority for the long-term survival of squirrel gliders in this area. Moreover, thin remnant corridors need to be conserved to promote gene flow and connectivity of the species populations in the urban matrix. If unsuccessful, squirrel glider populations may experience further isolation which could in turn lead to a loss of genetic diversity through genetic drift. The Lake Macquarie LGA is a stronghold for the threatened species in Australia, and as such this location warrants priority conservation management. High levels of genetic structure and population differentiation also warrants proactive conservation management in response.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons. org/licenses/by/4.0/.