Study area
Sugar gliders require tree hollows for sleeping and nesting and primarily rely on the sap, gum and nectar of eucalypt, acacia, and banksia species (Smith 1982; Lindenmayer 2002). Together the Lake Macquarie LGA and neighbouring Newcastle LGA contains 47, 100 ha of native forest (55% of total area), the majority of which is comprised of medium open eucalypt forest that is suitable for glider species (Department of Agriculture Water and Environment 2021). Due to the eucalypt, acacia and banksia habitat spread across the landscape, Lake Macquarie LGA holds the most abundant squirrel glider population in New South Wales (NSW) (Smith 2002) and similarly holds a large population of sugar gliders (Smith and Murray 2003). This location is recognized as being the most genetically diverse for squirrel gliders in Australia, and is located 130 km north of Sydney, NSW (Pavlova et al. 2010) (Fig. 1). Additionally, Lake Macquarie LGA contains the largest coastal saltwater lake in Australia, (120 km2) a potential biogeographical barrier to glider gene flow.
Live trapping and DNA collection
Live trapping was conducted at 32 sites from 2017 to 2020, with 113 sugar glider individuals caught across 12 of the sites (Fig. 1). A combination of Mawbey traps, Elliot B traps, cage traps and Winning and King pipe traps were used to live trap gliders (Mawbey 1989; Quin 1995; Winning and King 2008). The Winning and King pipe traps were secured one metre from the base of the tree and the bottom was filled with leaves for bedding and a bait ball made from a mixture of peanut butter, honey and oats (Winning and King 2008). The other traps were instead secured to wooden planks that were drilled into tree trunks two metres above the ground. Each of these trap types contained a mix of leaves for bedding and a bait ball. A 1:4 ratio of honey water was then sprayed up and down the tree to a height of six metres as well as around the entrance to each trap as an olfactory attractant (Sharpe and Goldingay 2007). 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 meters apart), and further focus was given to locations that required a larger sample size for intended genetic analyses. Trapping ceased when no new individuals were detected after one week or when a previously marked individual was recaught three times within a week.
Traps were checked each morning at sunrise. When a sugar glider was caught, body measurements were recorded. These measurements included weight, tail length, right hind foot length, head width, head length and sex to monitor the body condition of any recaptured animals. Individuals were given a unique identifying number in the form of a metal ear tag or a unique ear marking combination. Before being released, their DNA was collected in the form of an ear biopsy. A 2 mm metal ear punch (Able Scientific, Australia) was sprayed with 70% ethanol and flamed for sterilization. Once cooled, a small clipping was taken from the outside edge of the ear and stored in sterilized vials containing 95% ethanol. These were kept at – 20°C prior to DNA extraction. The processing of each sugar glider was limited to five-minutes to avoid unnecessary disturbance and following approved animal ethics protocol (AE15/11 and AE19/02). When the processing of an individual was complete, the glider was safely released onto the tree where it was caught.
SNP genotyping and filtering
Genomic DNA was extracted from ear tissue using DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany) as per the manufacturer’s instructions. A set of 96 samples were selected for analysis after DNA concentration and purity were assessed using agarose gel electrophoresis and a NanoDrop 2000 Spectrophotometer (Thermo Scientific). These samples were plated and transported to Diversity Arrays Technology (DArTseq), University of Canberra, Australia, for concentration and sequencing (Kilian et al. 2012). DArTseq technology has effectively assisted with the research of Australian marsupials in recent years, including koalas (Phascolarctos cinereus) (Schultz et al. 2018; Kjeldsen et al. 2019), bare-nosed wombats (Vombatus ursinus) (Martin et al. 2019) and Tasmanian devils (Sarcophilus harrisii) (Wright et al. 2019). DArTseq utilises next generation sequencing and genome complexity reduction methods to call bi-allelic, codominant SNP markers (Kilian et al. 2012). For this species, a double digest was conducted with PstI and SphI restriction enzymes. The PstI overhang was compatible with the barcode adapter while the SphI overhang was compatible with the reverse adaptor and the flowcell attachment region (Elshire et al. 2011; Kilian et al. 2012). The resulting DNA fragments were then amplified with PCR and sequenced on Illumina Hiseq2500. Reads were cleaned, barcodes were removed, and sequences were aligned to the Leadbeater’s possum reference genome (Gymnobelideus leadbeateri “, GCA_011680675.1 LBP_v1”) using DArTseq analytical pipelines. Once the process was complete, genome-wide SNPs were retained with a minimum sequence identity of 70%. For additional information regarding the DArTseq process, refer to Kilian et al. 2012.
Final filtering was conducted in R 4.0.2 (R Core Team 2015) using the dartR package (Gruber et al. 2018). SNPs were filtered on call rate (0.95 threshold) and hamming distance (0.2 threshold). Monomorphic loci and linked loci were removed, as well as those that significantly deviated from Hardy-Weinberg Equilibrium (p < 0.05).
Population structure and genetic diversity analyses
All analyses were undertaken in R 4.0.2 (R Core Team 2020) unless otherwise stated. For sampling locations with two or more individuals, average observed (Hobs) and expected (Hex) heterozygosity was calculated at each SNP using the package hierfstat (Goudet 2005). These values were used to assess the level of inbreeding (FIS = (Hex-Hobs)/Hex) at different geographic locations. Additionally, F-statistics in the form of average FIS and FST were generated for the overall region.
A Principal Component Analysis (PCOA) ordination enabled visualisation of genetic differentiation between individuals and between sampling locations (Gower 1966). This was generated with the gl.pcoa.plot function in dartR (Gruber et al. 2018). Similarly, Pairwise FST and an Analysis of Molecular Variance (AMOVA) were used to examine population structure and population differentiation. Pairwise FST and the corresponding p values were calculated between the 12 sampling locations using the gl.fst.pop function in the R package dartR 1.1.11. 100 bootstraps were performed across loci to generate the p-values. The AMOVA was conducted through the R package poppr 2.9.0 using 9,999 permutations (Zhian Kamvar 2021).
A STRUCTURE analysis examined individual-based structure and admixture. The program STRUCTURE 2.3.4 was employed to determine the appropriate number of clusters (K) and to visualise the ancestry proportions per individual (Pritchard et al. 2000; Pritchard et al. 2003). The program tested genetic clusters K = 1 to 15 with 8 repeats of K. Each run had a 10,000-length burn-in period followed by 10,000 Monte Carlo Markov Chain replications. The results were then uploaded to Structure Harvester (Web v0.6.94) and the most probable K was chosen based on the Evanno method and the largest Delta K value (Evanno et al. 2005; Earl and VonHoldt 2012). Individuals with q values > 0.8 were considered pure to a cluster while individuals with q values < 0.8 were considered admixed.
Isolation by distance and least cost path analyses
An isolation by distance (IBD) Mantel test with 9, 999 permutations assessed whether there was a relationship between geographic distance and genetic distance of populations. This analysis took the log of straight-line distances between populations using the gl.ibd function in the dartR package and compared them to genetic distances in the form of pairwise FST (that is, FST/(1-FST)).
Additionally, three least cost path analyses were performed using the gl.genleastcost function in dartR to account for habitat suitability and barriers. All rasters (friction matrixes) were created in ARCMAP 10.7.1. The first friction matrix accounted for biogeographical barriers by assigning “NoData” values to the ocean and lake from a land use spatial layer (NSW Government 2007), indicating a complete barrier to dispersal. Everything else was assigned a cell value of 1 to allow easy dispersal. This matrix was used in a Mantel test and will be hereon referred to as Mantel test B (biogeographical barriers).
The second friction matrix built on Mantel test B by assigning a cell value of 200 to four-lane highways obtained from the NSW Department of Industry, indicating unsuitable habitat and high effort for dispersal. This was used in a Mantel test and given the name Mantel test BH (biogeographical barriers and highways). The third and final friction matrix was created with cells ranging from optimal habitat (1) to unsuitable habitat (200) and given the name BHH (biographical barriers, highways and habitat). First, Plant Community Type (PCT) spatial layers were taken from Bell et al. (2016) and Eco Logical Australia Pty Ltd (2003) and vegetation types were divided into five categories with unique cell cost values based on the habitat and dietary preferences of sugar gliders (Smith 1982; Smith and Murray 2003). These included highly suitable vegetation (e.g. Eucalyptus haemastoma or racemose, Angophora costata, understorey of Banksia spp and Xanthorrhoea spp, cell cost value 1), suitable vegetation (10), moderately suitable vegetation (20), moderately unsuitable vegetation (75) and unsuitable vegetation (e.g. saltmarsh, sedgeland, cell cost value 100). Next, land use spatial layers were taken from the NSW Government (2007) and divided into the following classes and cell cost values: Urban residential (80), industrial (100), roads (100), railways (80), 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) and the Lake and Ocean (‘NoData’, complete biogeographical barrier to dispersal). Highly suitable vegetation was given a 25-meter buffer so that two patches either side of a barrier (e.g., road/river) would intersect with a glide distance threshold of 50 meters, indicating crossing potential. Major highways were overlayed with a value of 200 with the assumption that sugar gliders would avoid busy four-lane highways (Pacific Motorway, Pacific Highway).
New geographic distances were calculated using the three friction matrices (cell factor = 15, function = mean) and the least cost paths between spatial coordinates of individuals. The new pairwise geographic distances were then tested for correlation with genetic distances and number of neighbours = 8. The results from the four spatial analyses (IBD and three least cost paths) were compared to find the best fit, as recommended by Milanesi et al. (2016).