Genetic evidence supports three previously described species of greater glider, Petauroides volans, P. minor, and P. armillatus

The identification and classification of species are essential for effective conservation management. This year, Australia experienced a bushfire season of unprecedented severity, resulting in widespread habitat loss and mortality. As a result, there has been an increased focus on understanding genetic diversity and structure across the range of individual species to protect resilience in the face of climate change. The greater glider (Petauroides volans) is a large, gliding eucalypt folivore. This nocturnal arboreal marsupial has a wide distribution across eastern Australia and is considered the sole extant member of the genus Petauroides. Differences in morphology have led to suggestions that the one accepted species is actually three. This would have substantial impacts on conservation management, particularly given a recent history of declining populations, coupled with extensive wildfires. Until now, genetic evidence to support multiple species has been lacking. For the first time, we used DArT sequencing on greater glider tissue samples from multiple regions and found evidence of three operational taxonomic units (OTUs) representing northern, central and southern groups. The three OTUs were also supported by our morphological data. These findings have important implications for greater glider management and highlight the role of genetics in helping to assess conservation status.

Effective conservation management relies on accurate taxonomic classification and a robust understanding of the genetic structure of populations. A lack of knowledge about the genetic structure of species across their range can result in an inability to properly manage and protect species from extinction 1 . This is especially true in the wake of a natural disaster, when wildlife management decisions need to be made quickly and under challenging circumstances.
The catastrophic 2019-2020 bushfire season in Australia burnt over 97,000 km 2 and directly or indirectly killed millions of native animals 2 . The impacts of fire on genetic diversity are nontrivial 3 , and the extent of the recent fires means that substantial portions of many species' ranges were impacted 2 . As a result, the conservation status of a number of species are being revisited and the importance of protecting genetic diversity within species is becoming more widely recognised for its central role in conserving species' resilience to anthropogenic climate change [4][5][6] .
The greater glider (Petauroides volans, (Kerr 7 ) is a large, nocturnal gliding marsupial endemic to Australia (Fig. 1). Strictly arboreal, they are vulnerable to habitat loss and disturbance due to their specialised diet of eucalypt leaves and obligate dependence on mature trees with large hollows for shelter 8 . Greater gliders have a widespread distribution primarily associated with eucalypt forests along the Great Dividing Range from northern Queensland to southern Victoria 9 .
Data from long-term monitoring sites have revealed alarming declines and localized population extinctions of greater gliders over the past twenty years [9][10][11][12] . As a result, this once common and abundant species is now listed as vulnerable under the National Environmental Protection and Biodiversity Conservation Act 13 and vulnerable globally under the IUCN's Red List of threatened species 14 . Owing to the recent widespread fires and the history

Results
Summary of morphological data. We collected eight morphological measurements from 63 greater gliders sampled from five locations (Supplementary materials Table S1 and S2 available online). Morphological measurements were not available for the 12 museum samples (Supplementary Table S3). Only adult animals were used in the morphological analyses, resulting in 56 individuals. Data on sex based differences from other studies are equivocal 18,25 . Using linear models for unbalanced data with backward elimination method, we found that measurements for males and females in our dataset were not significantly different (p > 0.05) in seven out of eight attributes, with the exception that females had longer tails (F 1,52 = 7.14, p = 0.01). We combined male and female data within locations for the remaining analyses. A principal components analysis (PCA) of morphological measurements broadly clustered into three groups with two samples showing some overlap (PC1 = 70%, PC2 = 8.3%, Supplementary Fig. S1 and Table S4). Body mass was the strongest contributor to differences among  We further investigated differences between greater glider morphology across locations using Tukey's posthoc multiple comparison test. We found no significant difference in greater glider body measurements between the two southern sites (Wombat and Bendoc) for any of the eight measured attributes (Supplementary Table S5). Similarly, greater gliders from the two northern sites (Taravale and Blackbraes) differed only in that Taravale gliders had significantly longer tails (t 54 = 4.99, p < 0.0001, Supplementary Table S5). We then collapsed the two southern sites into a single region called Southern and the two northern sites into a single region called Northern, and explored morphological differences in greater gliders among the three regions (Northern, Central, and Southern; summary statistics provided in Table 1). The Northern and Southern greater gliders differed significantly in every measured morphological attribute (Table 2). Northern and Central greater gliders differed significantly in every measured attribute except body length (t 54 = − 0.23, p = 0.971), and Central and Southern greater gliders differed significantly in every attribute except ear length (t 54 = − 5.92, p = 0.502) and ear width (t 54 = 0.97, p = 0.595; Table 2). Northern greater gliders were the smallest of the three groups in almost every measured attribute and Southern greater gliders were the largest, with the Central intermediate, but still significantly different in most measured attributes from the Northern and Southern groups. Descriptions of pelage colouration, which was distinctly different between regions, can be found in the supplementary materials (Supplementary Table S2 and Fig. S2).
Summary of geographic population structuring and metadata. Genetic population structure analysis revealed three diagnosable aggregations that we refer to as operational taxonomic units (OTUs). The raw data consisted of 78 individuals scored for SNP polymorphisms at 45,245 loci. Filtering, as outlined in Methods, resulted in only high-quality SNP loci being retained with high coverage across individuals and loci, ensuring completeness of the final dataset and minimising the amount of missing data. This yielded SNP calls for 11,317 loci for 75 individuals. Following this filtering, these data are considered to be highly reliable.
Principal components analysis of this dataset revealed three distinct genetic groups corresponding to geographic locations of Taravale/Blackbraes (Northern), Bendoc and Wombat (Southern) and Redcliffe Vale/ Museum (Central) (Fig. 3). Two individuals from Taravale (Northern) lay intermediate in position between the  Excluding the putative hybridizations, fixed differences analysis 26,27 identified four putative operational taxonomic units (OTUs) diagnosable on the basis of one or more fixed allelic differences-the Northern (Taravale and Blackbraes combined), and the two Southern, Bendoc and Wombat. Wombat (n = 6) and Bendoc (n = 9) were distinguished by only 2% (233) fixed differences in comparison with the other populations that differed between 11 and 30% (1373-4767). Test of significance demonstrated that the difference between Wombat and Bendoc was not significant (i.e. did not exceed the false positive rate for fixed differences, given the respective sample sizes-see Georges et al. 26 : R scripts are provided in the online data repository). We therefore combined Wombat and Bendoc to yield a final three putative OTUs, representing North, Central, and Southern populations (Fig. 2). Percent fixed differences between OTUs are presented in Table 3 and counts of private alleles between the final OTUs taken pairwise are presented in Table 4. Southern greater gliders had the highest expected heterozygosity (0.11), a measure of genetic diversity, substantially higher than Northern gliders (0.04, p < 0.0001) and marginally higher than Central gliders (0.10, p < 0.02). Northern gliders had the lowest heterozygosity (0.04), significantly lower also than Central gliders (0.10, p < 0.0001).
Hybrids exist between the two Queensland genetic groups. The two individuals collected at Taravale (Northern group) had intermediate positions in the PCA falling between the Northern group and the Central group. Comparing these individuals with likelihood bins associated with parental genotypes, F1 hybrids, backcrosses of the F1 to the parentals, and F2 hybrids using New Hybrids 28 showed that one individual had a genotype consistent with an F1 hybrid between the Northern group and the Central group (T1) and the other a genotype consistent with a backcross of an F1 hybrid to the Northern group (T5) ( Table 5). We suspected based on the proximity of capture, ages, and evidence of an F1 and backcross genotypes, that individuals T1 and T5 might be parent-offspring. Analysis of pedigree incompatibility estimates for this pair and pairs taken at random from our pool of individuals, indicated that a parent-offspring relationship between T1 and T5 was highly likely. Table 2. Tukey's Post-hoc test between regions for each measured morphological trait. ****p < 0.0001, ***p < 0.001, **p < 0.01, *p < 0.05.  www.nature.com/scientificreports/ The parental assignments arose naturally from the analysis, and were not defined a priori. The analysis identified one F1 (T1) and one F1 × P0 (T5).

Discussion
Greater gliders are currently classified as a single species, P. volans 13,14 . We used DArTseq on samples from multiple locations across the range of P. volans, which revealed three operational taxonomic units (OTUs) representing Northern, Central and Southern groups. These broadly correspond with the species designations summarised by Jackson and Groves 20 and provide the first genetic support for multiple greater glider species. The morphological measurement from the sampled greater gliders from the five study areas also align with these three OTUs. The Northern and Southern populations differed significantly in all morphological measurements, and this is in keeping with previously reported differences in the ear, body size, and pelage colouration of the currently recognised subspecies (P. volans minor and P. volans volans) 16 . There was some overlap between animals from the Central population in Redcliffe Vale and the other two groups in morphology, but the Central population still differed significantly in a number of measurements from the Northern and the Southern animals. We also found evidence of hybridisation between the Northern and Central groups and one of those hybrids had morphological measurements that were intermediate between those groups as well. As discussed by de Queiroz 29 , delimiting species is complicated by the numerous species concepts. The Genetic Species Concept recognises that species are genetically isolated but not reproductively isolated, allowing for hybridisation, which is likely to occur during the speciation process 30 . Petauroides volans has a long and complicated taxonomic history, which has been extensively reviewed by Maloney 31 and Maloney and Harris 32 . The greater glider was first identified to western science by Kerr 7 , with further variability in morphology among collected individuals described by Collett 18 who suggested a "variety" called Petauroides volans minor based on differences in skull size, body length and pelage colour. The greater Table 3. Percent fixed differences between populations taken pairwise (lower matrix and diagonal, loci = 11,317) and average number of individuals scored (upper matrix) to estimate the fixed differences. Note that Taravale and Blackbraes have no fixed allelic differences and so are amalgamated into a single OTU (Northern); the 2% fixed difference between Bendoc and Wombat were not statistically significant and so were amalgamated into a single OTU (Southern).    T2, B9, T3, B2, B4, T17, B12, T7, B14, T18, B13, B3, T15, B5, T9, B1, T6, T11, B7, T14, T10, T4, Bjuv, T13, T12,  T16, B15, B10, T8 20,22 . Based on personal communications with the late Dr Kennith Aplin and unpublished data in a conference presentation by Arbogast et al. 23 , Jackson and Groves 19,20 recognised three species within Petauroides and two subspecies. Until now, there has been no published genetic evidence to support either multiple or single species designations, and P. volans has been legislated for the purpose of management and conservation as a single species 13,14 . Our results support a population structure consistent with the multi-species designation of Jackson and Groves 20 , with a Northern, Central and Southern species of greater glider. We are not able to resolve subspecies structure from our data. While the three OTUs roughly align with the species designations described by Jackson 19 , there are some notable differences. In particular, the Central OTU that corresponds with the proposed "central" species, P. armillatus, extends much farther south than anticipated based on the distribution described in Jackson 19 . Although we had a relatively small, non-continuous sample, it is also worth noting that we found no evidence of the Southern OTU that corresponds with P. volans in Queensland, despite the range proposed by Jackson 19 extending to Bundaberg. There is no information provided by Jackson 19 to explain how the distributions in that text were determined. Notably, Jackson's distribution for P. volans did not agree with the range proposed by his source, Arbogast 23 , who does not extend the distribution of P. volans into Queensland, and that agrees with our data.
An important caveat is that the findings from this study arose in the process of conducting unrelated research to look at relationships between variations in physiology and climate at the northern and southern distribution of what was thought to be one species at the time (McGregor et al., in prep). Given the distinct morphological differences between the northern and southern samples, we decided to use genetics to investigate the population structure. We did not comprehensively sample greater gliders throughout their range and no samples were obtained in New South Wales. Hypothesised species can represent arbitrary slices across geographic clines 33 , and we have been cautious in interpreting our data in the context of existing work that defines the three species and draws support from the SNP data 19,20 . In particular, we have not attempted to describe the new taxa de novo using genetic data alone, but rather to examine the genetic data to see if it is consistent or inconsistent with the taxonomic concepts put forward by Jackson and Groves 19,20 of three species. Our data are certainly inconsistent with the samples representing panmictic or close to panmictic populations, and in the view of the magnitude of the observed fixed differences and the very substantial counts of private alleles held in each of our three OTUs, the genetic divergences we observed are very unlikely to be explained by punctuated sampling along a cline. Additional surveys should be undertaken to obtain a more complete understanding of genetic structure across the range of greater gliders. While our data does not enable an exhaustive assessment of Petauroides systematics, it does provide the first genetic evidence to support multiple extant species in Petauroides, and this finding has important conservation implications.
Over the past decade, several studies have reported alarming declines in greater glider populations in the Blue Mountains, NSW and Central Highlands, Vic and localized extinctions in other areas, such as Booderee NP, NSW 9-12 . Population declines have been attributed to the cumulative effects of land clearing and climate change, including logging, lower than average rainfall and increased occurrence of intense fires 11,[34][35][36] . This year, the most destructive wildfire season on record burned through substantial portions of greater glider habitat in south-eastern Australia. The extent of greater glider habitat impacted by fire and repercussions to populations are still under investigation. However, as a result of the recent widespread habitat destruction and the history of unresolved population declines, the conservation status of the currently recognised single species, P. volans, is being reassessed.
The division of the greater glider into multiple species reduces the presumed, widespread distribution of P. volans, further increasing conservation concern for this species and highlighting the lack of information about the other Petauroides spp. Ecological studies have focused primarily on populations in NSW and Vic with only a few published studies of P. minor 16 and no prior studies of P. armillatus. The knowledge that there is now genetic support for multiple species, with distributions that are probably much smaller than the previously designated range of P. volans, should be a consideration in future conservation status decisions and management legislation.

Methods
Sample collection. We sampled tissue from wild greater gliders from four locations that broadly represent the northern and southern distribution of the Petauroides latitudinal and longitudinal geographic range as part of a separate study to investigate relationships between animal physiology and climate (Fig. 4). The sites included Mount Zero-Taravale Australian Wildlife Sanctuary (19° 07′ 18″ S, 146° 04′ 42″ E, n = 18) and Blackbraes National Park (19° 34′ 39″ S, 144° 05′ 05″ E, n = 15) in North Queensland, and Bendoc State Forest (37° 10′ 35″ S, 148° 56′ 52″ E, n = 9) and Wombat State Forest (37° 29′ 50″ S, 144° 09′ 23″ E, n = 6) in Victoria. We then conducted additional field sampling in Redcliffe Vale, Queensland (21° 06′ 57″ S, 148° 06′ 58″ E, n = 18) (Fig. 4) as this was the suggested location of the proposed P. armillatus 19,20 . In addition, 12 museum specimen tissue samples were obtained from the Queensland Museum from greater gliders collected from northern, central and southern Queensland to investigate genetic structure in that area more broadly (Fig. 4). Additional information about the climate, geography and vegetation at each sample location can be found in Supplementary Material, (Table S1 online).
Wild greater gliders were located through spotlight searches using high-powered, handheld torches (Ledlenser P7, Zweibrüder Optoelectronics GmbH and Co., Solingen, Germany) to detect greater glider eye shine. Individuals were then captured using a gas-powered, tranquilizer dart-gun (Montech Black Wolf; Tranquil Arms

Morphology data investigation. Principal components analysis of the eight measured morphological
traits for greater gliders from the five sites (Taravale, Blackbraes, Redcliffe Vale, Bendoc, and Wombat) was performed in R, using the "prcomp" function. The plot was generated using the ggplot2 package 37 . We then used a canonical variate analysis (CVA) in R package MASS 38 to analyse regional group structure (Northern, Central, and Southern) in the multivariate data. To explore whether there were significant differences in measured traits between sexes and account for unequal sample sizes, we used linear models with backward elimination variable selection method and Tukey's post-hoc multiple comparison test. Linear models included each morphological trait as the response variable and region (Northern, Central, Southern), sex (male, female) and the interaction between region and sex. Insignificant variables were eliminated until only significant variables remained. Based on these results, sex was pooled and Tukey's post-hoc multiple comparison test was used to compare measured morphological traits between regions. We also explored differences between sites with Tukey's post-hoc multiple comparison test.   www.nature.com/scientificreports/ Ltd, Canberra, Australia), which uses a combination of complexity reduction using restriction enzymes, implicit fragment size selection and next generation sequencing 39 , as described in detail by Kilian et al. 40 and Georges et al. 26 . Essentially, DArTseq is an implementation of sequencing complexity-reduced representations 41 and more recent applications on next generation sequencing platforms 42,43 . To achieve the most appropriate complexity reduction (the fraction of the genome represented, controlling average read depth and number of polymorphic loci), four combinations of restriction enzymes (Pstl enzyme combined with either Hpall, Sphl, Nspl or Msel) were evaluated and restriction enzyme combination of Pstl (recognition sequence 5′-CTGCAIG-3′) and Sphl (5′-GCATGIC-3′) was selected. Amplification using polymerase chain reaction (PCR) 26,44 and the conditions applied are as described in Georges et al. 26 . After PCR, equimolar amounts of amplification products from each sample were pooled and applied to cBot (Illumina) bridge PCR for sequencing on the Illumina Hiseq 2500. The sequencing (single end) was run for 77 cycles to yield sequence tags of 20-69 bp after removing adaptors. SNP genotyping. Sequences generated from each lane were processed using proprietary DArT Pty Ltd analytical pipelines as described by Georges et al. 26 . In particular, one third of samples were processed twice from DNA, using independent adaptors, to allelic calls as technical replicates, and scoring consistency (repeatability) was used as the main selection criterion for high quality/low error rate markers. The DArT analysis pipelines have been tested against hundreds of controlled crosses to verify Mendelian behaviour of the resultant SNPs as part of their commercial operations. The resultant data set contained the SNP genotypes and various associated metadata of which CloneiD (unique identity of the sequence tag for a locus), repAvg (proportion of technical replicate assay pairs for which the marker score is identical), CallRate (proportion of individuals scored at a particular locus) and SnpPosition (position in the sequence tag at which the defined SNP variant base occurs) are of particular relevance to our analyses.
Additional SNP filtering. The SNP data and associated metadata were read into a genlight object ({ade-genet} 45 ) to facilitate processing with package dartR (version 1.8) 46 . We first removed all but one SNP from each sequence tag (12,782 SNPs removed) and retained only those loci supported by a read depth between 5 × and 100 × (8895 loci removed). Three individuals were removed from the dataset owing to an exceptionally poor call rate of less than 0.5 (MS9, MS12, MS8) and resultant monomorphic loci removed from the dataset. Loci with a repeatability less than 0.99 were removed (2380 loci) and finally, loci with a call rate of less than 0.95 were removed (9870 loci). We regard the data remaining after this additional filtering (11,317 SNP markers for 75 individuals) as highly reliable.
Visualization. Genetic similarity among individuals, populations and colour morphs was visualized using ordination (principal coordinates analysis or PCoA 47 ) as implemented in the gl.pcoa and gl.pcoa.plot functions of dartR. A scree plot of eigenvalues guided the number of informative axes to examine 48 , taken in the context of the average percentage variation explained by the original variables (using the gl.pcoa.scree function in dartR).
Diagnosable units. Diagnosable units are populations or aggregations of populations that can be diagnosed by one or more fixed allelic differences. Individuals that were intermediate between two aggregations in the PCoA (T1 and T5) were removed as suspected hybrids or introgressed individuals and examined separately. A fixed difference analysis as implemented in dartR was undertaken, and pairs of populations with no fixed allelic differences were progressively amalgamated to yield a putative set of diagnosable operational taxonomic units (OTUs). Putative OTUs that were distinctive because of false positives (owing to sampling error) were identified using gl.fixed.diff() with test = TRUE in dartR (see Georges et al. 26 , R scripts are provided in the online data repository), and pairs of populations for which the count of fixed differences did not exceed the estimated false positive rate were also amalgamated. Counts of private alleles were obtained from gl.report.pa in dartR.
Genetic diversity. Expected heterozygosity was used as a measure of relative genetic diversity. Heterozygosity was obtained for each population from allele frequencies using the gl. report.heterozygosity function of dartR and pairwise comparisons of heterozygosity between populations were tested for significance using the gl.test.heterozygosity function in dartR (significance evaluated by re-randomizaton with 10,000 replicates).
Hybridisation. The genotypes of suspected hybrids/introgressed individuals (T1 and T5) were examined using New Hybrids 28 without specifying parental source populations. Briefly, New Hybrids uses simulation to characterize likelihood bins for each of the parental populations, F1 hybrids generated by crossing the parentals, F2 hybrids, and backcrosses between the F1 hybrids and the parental populations. These bins are used to estimate the likelihood of an individual belonging to each bin, and these likelihoods are rendered and scaled to produce posterior probabilities of bin membership. Parameters were set as ThetaPrior = Jeffreys, PiPrior = Jeffreys, burnin = 10,000, sweeps = 10,000 and the default genotype frequency classes (P0, P1, F1, F2, F1 × P0, F1 × P1). New Hybrids gives a posterior probability of individual membership in each of the genotype frequency classes, allowing effective assignment of first generation hybrids (F1), second generation hybrids (F2) and backcrosses of F1 hybrids to the parental populations.

Relatedness. Two individuals with genotypes intermediate between Taravale/Blackbraes and Redcliffe Vale
(T1 and T5) were found in the same cluster of trees in Taravale. Given their status and proximity, we wished to examine the possibility that they were in a parent-offspring relationship. We did this by assuming they were Scientific Reports | (2020) 10:19284 | https://doi.org/10.1038/s41598-020-76364-z www.nature.com/scientificreports/ parent-offspring and examining the frequency of pedigree inconsistencies (e.g. one individual homozygous reference, the other homozygous alternate). Theoretically, the count of pedigree inconsistencies should be zero for true parent-offspring pairs, but errors in calling the SNPs generates false positives. To overcome this, we generated counts of pedigree inconsistencies for all pairs of individuals to generate a null expectation and identified putative pairs of parent-offspring as outliers using software to be included in the next release of dartR (gl.report. parent.offspring).

Data availability
The datasets generated and analysed during the current study have been uploaded as related manuscript files that are available to the editors and reviewers. The raw data and all scripts we used for analysing it are publicly available through The Tropical Data Hub Research Data Repository (https:// resea rch. jcu. edu. au/ resea rchda ta) at https:// doi. org/ 10. 25903/ xx67-2m46.