pmc.ncbi.nlm.nih.gov

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

  • ️Mon Oct 04 2021

Abstract

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.

Subject terms: Conservation biology, Genomics

Introduction

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 extinction1. 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 km2 and directly or indirectly killed millions of native animals2. The impacts of fire on genetic diversity are nontrivial3, and the extent of the recent fires means that substantial portions of many species’ ranges were impacted2. 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 change4–6.

The greater glider (Petauroides volans, (Kerr7) 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 shelter8. Greater gliders have a widespread distribution primarily associated with eucalypt forests along the Great Dividing Range from northern Queensland to southern Victoria9.

Figure 1.

Figure 1

Greater gliders from the northern (top left), central (bottom left) and southern (right) groups identified through DArTseq showing morphological differences that are typical of our dataset. Greater gliders of the type shown on the right have several pelage colour morphs including white and light grey. Photos by Denise McGregor (top left) and Jasmine Vink (bottom left and right).

Data from long-term monitoring sites have revealed alarming declines and localized population extinctions of greater gliders over the past twenty years9–12. As a result, this once common and abundant species is now listed as vulnerable under the National Environmental Protection and Biodiversity Conservation Act13 and vulnerable globally under the IUCN’s Red List of threatened species14. Owing to the recent widespread fires and the history of unresolved population declines, the conservation status of the currently recognised single species, P. volans, is being reconsidered. The division of the greater glider into multiple species would have substantial conservation implications. In particular, population declines and recent fire impacts, which are mostly associated with populations south of Queensland, may extend across a much greater proportion of the distribution of P. volans than previously thought, and have a more substantial combined effect.

Petauroides volans is regarded as the sole member of the genus Petauroides and the only gliding member of the family Pseudocheiridae13,14. There are currently two subspecies of P. volans recognised by the Australian Government for the purpose of conservation legislation, P. v. minor and P. v. volans13. They have been described on the basis of morphological differences in body size, pelage colouration, and geographic distribution15,16. Petauroides volans volans is found south of the Tropic of Capricorn, weighs between 1000–1700 g and occurs in several colour morphs ranging from completely black, to black with a white belly and chest, to lighter grey dorsal pelage or white all over15,17. Petauroides volans minor has a relatively small body mass of 650–1100 g, occurs from the Tropic of Capricorn to just north of Cairns, and also has light ventral pelage, but the dorsal hair has dark brown tips and lighter subterminal bands16,18.

Recently, Jackson19 and Jackson and Groves20 recognised P. volans as comprising three distinct species, geographically delineated across the range of Petauroides. This raises the subspecies P. volans minor, to species status (P. minor) and resurrects the name P. armillatus, to apply to populations in mid-Queensland from the Eungella Range to just north of Townsville19. Petauroides armillatus was once considered a subspecies (P. volans armillatus) by Thomas in 192321 but subsequently synonymized within minor22. The multi-species taxonomy of Jackson19 and Jackson and Groves20 retains Petauroides volans as the most southerly distributed species from Bundaberg in Queensland (Qld) to Victoria (Vic), and includes two subspecies (P. volans incus and P. volans volans).

The claim that the three taxa P. volans, P. minor and P. armillatus are morphologically and genetically distinct at the level of species has been controversial because it lacks supporting published evidence and was referenced by Jackson and Groves20 from a personal communication with the late Dr Kenneth Aplin and a brief conference abstract by Arbogast et al.23. Neither of these sources have published data to support the distinctiveness of the proposed P. minor, P. armillatus, and P. volans. For the first time, we test Jackson and Groves’s20 three species designation using molecular sequence data, DArTseq24, and report the genetic population structure of greater gliders at multiple points across their natural range.

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 equivocal18,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 (F1,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 groups (Supplementary Table S4 and Fig. S1). A canonical variates analysis (CVA) of the measured traits from the three groups representing samples collected from the northern sites (Taravale and Blackbraes), central site (Redcliffe Vale), and southern sites (Bendoc and Wombat) demonstrated a high degree of separation (Fig. 2, CV1 = 88.4% and CV2 = 11.6%).

Figure 2.

Figure 2

Canonical variate analyses (CVA) of body measurements of 56 greater gliders from three regions (Northern, Central, and Southern). The shaded circles are the group mean.

We further investigated differences between greater glider morphology across locations using Tukey’s post-hoc 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 (t54 = 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 (t54 = − 0.23, p = 0.971), and Central and Southern greater gliders differed significantly in every attribute except ear length (t54 = − 5.92, p = 0.502) and ear width (t54 = 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).

Table 1.

Means and (standard deviations) of the measured morphological traits from sampled greater gliders across three regions. Mass unit is in kilograms, body length and tail length are measured in centimetres and all other measured traits are in millimetres.

Region N Mass Head length Head width Body length Tail length Knee to heal Ear length Ear width
Northern 27 0.68 (0.06) 57.21 (2.99) 34.19 (1.61) 27.89 (2.47) 40.78 (3.23) 104.19 (4.91) 35.9 (2.06) 26.57 (3.49)
Central 15 0.87 (0.09) 64.17 (3.21) 37.26 (1.51) 28.27 (3.31) 45.02 (2.13) 112.59 (5.02) 40.37 (4.02) 29.88 (2.87)
Southern 14 1.36 (0.16) 67.05 (2.34) 39.14 (2.79) 36.11 (2.95) 51.75 (2.82) 121.04 (3.82) 41.62 (3.28) 31 (2.46)

Table 2.

Tukey’s Post-hoc test between regions for each measured morphological trait.

Comparison Northern-Central Central-Southern Northern–Southern
Diff. in mean (SE) t value p value Diff. in mean (SE) t value p value Diff. in mean (SE) t value p value
Body mass − 0.19 (0.03) − 5.66 **** 0.50 (0.04) 13.13 **** 0.68 (0.03) 20.35 ****
Head length − 6.95 (0.94) − 7.43 *** 2.88 (1.08) 2.67 * 9.83 (0.96) 10.27 ***
Head width − 3.08 (0.63) − 4.91 **** 1.88 (0.72) 2.59 * 4.95 (0.64) 7.72 ****
Body length − 0.21 (0.91) − 0.23 0.97 7.84 (1.05) 7.45 **** 8.05 (0.93) 8.63 ****
Tail length 5.65 (0.93) − 6.10 **** 6.73 (1.07) 6.31 **** 12.38 (0.95) 13.09 ****
Knee to heel − 10.77 (1.51) − 7.11 **** 8.45 (1.75) 4.84 **** 19.21 (1.55) 12.41 ****
Ear length − 5.71 (0.97) − 5.92 **** 1.25 (1.11) 1.12 0.50 6.96 (0.99) 7.06 ****
Ear width − 4.31 (1.00) − 4.32 *** 1.12 (1.15) 0.97 0.60 5.44 (1.02) 5.32 ****

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 Northern cluster and the Central cluster, and are likely instances of contemporary hybridization, a proposition that we test in a separate analysis.

Figure 3.

Figure 3

Principal component analysis (PCA) of genetic data coloured by collapsed sample regions. PC1 explains 59.5% of the variation in the dataset and PC2 explains 20.8% of the variation. The potential hybrids are the samples from the Northern group that are plotted between that group and the Central group.

Excluding the putative hybridizations, fixed differences analysis26,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.

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).

Sites Bendoc Blackbraes Redcliffe Taravale Wombat
Bendoc 0 24 26.7 25 14.9
Blackbraes 41 0 32.7 31 20.9
Redcliffe 17 26 0 33.7 23.7
Taravale 39 0 23 0 21.9
Wombat 2 33 10 31 0

Table 4.

Private alleles between final OTUs, taken pairwise, for OTUs vertically listed (lower matrix)) and horizontally listed (upper matrix). This is based on 11,317 loci, that is, post filtering. We did not filter on minor allele frequency. The count of private alleles includes the count of fixed differences.

Northern Central Southern
Northern 0 4307 5149
Central 5717 0 4035
Southern 6527 4003 0

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 Hybrids28 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. 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).

Table 5.

Results of the New Hybrids analysis. IDs, the identifiers used for each individual glider; P0, parental population Northern; F1, F1 hybrid; F2, F2 hybrid; F1 × P0, backcross of an F1 with the Northern parental population; F1 × P1, backcross of an F1 with the Central parental population. Count is the number of individuals tallied for each class. The body of the table contains the posterior probabilities of group membership, calculated from the likelihood distributions generated by simulation by New Hybrids28.

IDs Pop Count P0 P1 F1 F2 F1  × P0 F1 × P1
T5 Northern 1 0 0 0 0 1 0
T1 Northern 1 0 0 1 0 0 0
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, B6, B8, MS2 Northern 32 1 0 0 0 0 0
ER1. E9. ER17, ER2, ER10, ER10A, ER3, ER11, ER4, ER12, ER5, ER13, ER6, ER14, ER7, ER15, ER8, ER16, MS3, MS7, MS1, MS10, MS11, MS4, MS5, MS6 Central 26 0 1 0 0 0 0

Discussion

Greater gliders are currently classified as a single species, P. volans13,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 Groves20 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 Queiroz29, 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 process30.

Petauroides volans has a long and complicated taxonomic history, which has been extensively reviewed by Maloney31 and Maloney and Harris32. The greater glider was first identified to western science by Kerr7, with further variability in morphology among collected individuals described by Collett18 who suggested a “variety” called Petauroides volans minor based on differences in skull size, body length and pelage colour. The greater glider was more formally described as one species by Oldfield Thomas21 with four subspecies (Petauroides volans typicus, P. v. minor, P. v. incanus and P. v. armillatus)31. Petauroides v. armillatus was then synonymised with P. v. minor in 1934 by Iredale and Troughton who also changed the genus to Schoinobates resulting in three renamed subspecies S. v. minor, S. v. incanus, and S. v. volans (previously typicus)19,20. Petauroides was restored by McKay in 1982 who also synonymised P. volans incanus with P. v. volans in 198820,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 Groves19,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 species13,14. Our results support a population structure consistent with the multi-species designation of Jackson and Groves20, 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 Jackson19, 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 Jackson19. 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 Jackson19 extending to Bundaberg. There is no information provided by Jackson19 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, Arbogast23, 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 clines33, 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 data19,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 Groves19,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, NSW9–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 fires11,34–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. minor16 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. armillatus19,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).

Figure 4.

Figure 4

Location of the five study areas in eastern Australia in blue triangles (Blackbraes National Park (NP), Mount Zero-Taravale Sanctuary, Redcliffe Vale, Bendock State Forest (SF) and Wombat SF), and location of museum samples in orange squares (MS). The grey shading represents the current distribution of the greater glider from the Australian Species of National Environmental Significance Database49. This map was generated in R using the Australian coastline data from Geoscience Australia50 and multiple R packages (ggplot239, ggsn51, sp52,53, rgdal54, raster55).

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 Company, VIC, Australia) and darts specifically designed for mammals between 400 and 2000 g (0.5 ml; Minidarts, Tranquil Arms Company) containing 30–60 mg Zoletil 100 (Zolazepan and Tiletamine 50:50; Virbac), depending on the estimated body mass of the animal. While still under sedation, captured individuals were weighed and sexed, and reproductive status was assessed. External measurements were taken with vernier callipers using external jaws to measure head length (tip of snout to occipital bone protuberance), head width (widest part of one zygomatic arch to the analogous location on the other side of the head), ear width (widest part of the ear when flattened), ear length (from tragus to the outermost edge of the ear, excluding fur), and knee to heal-hind limb (top of knee to base of the heel with limb flexed to ninety degrees). Body length (occipital bone protuberance to the base of the tail, following the spine, with head in-line with the plane of the body) and tail length (cloaca to the tip of the tail, excluding fur) were measured with a flexible tape measure. Each individual was marked with a PIT tag (AVID Microchip Company, CA, USA) implanted subcutaneously and a tissue sample was clipped from the margin of the ear for DNA analysis. All work involving live animals complied with animal ethics and relevant guidelines and regulations. The animal capture and tissue collection was approved by James Cook University (Animal Ethics Permits A2137, A2545).

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 package37. We then used a canonical variate analysis (CVA) in R package MASS38 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.

DNA extraction and sequencing

DNA was extracted by Diversity Arrays Technologies (DArT Pty Ltd, Canberra, Australia) using a NucleoMag 96 Tissue Kit (MachereyNagel, Duren, Germany) coupled with NucleoMag SEP (Ref. 744900) to allow automated separation of high-quality DNA on a Freedom Evo robotic liquid handler (TECAN, Miinnedorf, Switzerland). Tissue was first incubated overnight with proteinase K, adjusted in concentration depending on the tissue. Sequencing for SNP genotyping was done using DArTseq (DArT Pty Ltd, Canberra, Australia), which uses a combination of complexity reduction using restriction enzymes, implicit fragment size selection and next generation sequencing39, as described in detail by Kilian et al.40 and Georges et al.26. Essentially, DArTseq is an implementation of sequencing complexity-reduced representations41 and more recent applications on next generation sequencing platforms42,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 ({adegenet}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 PCoA47) 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 examine48, 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 Hybrids28 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 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).

Supplementary information

Acknowledgements

We acknowledge the Traditional Custodians of the land on which the fieldwork was conducted and pay our respects to elders past and present. We thank A. Hemphill, B. Baker, M. Furumi, J. Scott, A. Feck and other volunteers for assistance with fieldwork. We also thank S. Banks, B. Moore, C. Meakin, the Queensland Museum for providing greater glider tissue samples, the Australian Wildlife Conservancy and Regional Operations manager T. White, Queensland Parks and Wildlife Service and especially G. Anchen, and E. Adams, Victoria Environment and Natural Resources and T. Mitchell and A. Kilian for their help with the project. This study was funded by the Winifred Violet Scott Charitable Trust, Skyrail Rainforest Foundation, Wet Tropics Management Authority and the Holsworth Wildlife Research Endowment. Research was conducted under permit by Victoria Department of Environment, Land, Water and Planning (10007842), and Queensland Scientific Purposes Permits (WITK16408715, WISP16408815, WIF416492015, WITK18792718, WISP18803718, WIF418792618).

Author contributions

D.M. collected field data from wild greater gliders, obtained museum sample specimens, collated data, and secured funding support. D.M., A.P., A.G., H.Y., A.K., and K.Y. contributed to manuscript text. A.G. undertook DArTseq analyses and A.P. and A.G. interpreted results from those analyses. H.Y. and K.Y. undertook morphological analyses. K.Y. prepared Fig. 1, H.Y. prepared Fig. 2, A.G. prepared Figs. 3, and A.P. prepared Fig. 4. A.K. assisted with fieldwork preparation and contributed to funding acquisition. A.K. and K.Y. assisted with project development, research support and supervision. All authors reviewed the full manuscript.

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://research.jcu.edu.au/researchdata) at https://doi.org/10.25903/xx67-2m46.

Competing interests

The authors declare no competing interests.

Footnotes

The original online version of this Article was revised: The original version of this Article and accompanying Supplementary Information file contained a repeated error, where the co-ordinates given for the Redcliffe Vale, Queensland study site “(21° 06′ 57″ S, 148° 06′ 58″ E, n = 18)” was incorrectly given as “(21° 06′ 57″ S, 148° 56′ 58″ E, n = 18).”

Publisher's note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Change history

10/4/2021

A Correction to this paper has been published: 10.1038/s41598-021-99388-5

Supplementary information

is available for this paper at 10.1038/s41598-020-76364-z.

References

  • 1.Gutiérrez EE, Helgen KM. Mammalogy: outdated taxonomy blocks conservation. Nature. 2013;495:314. doi: 10.1038/495314e. [DOI] [PubMed] [Google Scholar]
  • 2.Ward M, et al. Impact of 2019–2020 mega-fires on Australian fauna habitat. Nat. Ecol. Evolut. 2020 doi: 10.1038/s41559-020-1251-1. [DOI] [PubMed] [Google Scholar]
  • 3.Potvin DA, et al. Genetic erosion and escalating extinction risk in frogs with increasing wildfire frequency. J. Appl. Ecol. 2017;54:945–954. doi: 10.1111/1365-2664.12809. [DOI] [Google Scholar]
  • 4.Walters AD, Schwartz MK. Population genomics for the management of wild vertebrate populations. In: Rajora OP, editor. Population Genomics. Berlin: Springer; 2020. [Google Scholar]
  • 5.Luikart G, Hohenlohe PA, Allendorf FW. Genomics and the future of conservation genetics. Nat. Rev. Genet. 2010;11:697–709. doi: 10.1038/nrg2844. [DOI] [PubMed] [Google Scholar]
  • 6.Sgrò CM, Lowe AJ, Hoffmann AA. Building evolutionary resilience for conserving biodiversity under climate change. Evolut. Appl. 2011;4:326–337. doi: 10.1111/j.1752-4571.2010.00157.x. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 7.Linné C., Gmelin J. F., Kerr, R. & Archer, J. The Animal Kingdom, or Zoological System, of the Celebrated Sir Charles Linnæus. Containing a Complete Systematic Description, Arrangement, and Nomenclature, of All the Known Species and Varieties of The Mammalia, or Animals Which Give Suck to Their Young. (Edinburgh: Printed for A. Strahan, and T. Cadell, London, and W. Creech, 1792).
  • 8.Tyndale-Biscoe H. Life of Marsupials. Collingwood: CSIRO Publishing; 2005. [Google Scholar]
  • 9.Lindenmayer DB, et al. How to make a common species rare: a case against conservation complacency. Biol. Conserv. 2011;144:1663–1672. doi: 10.1016/j.biocon.2011.02.022. [DOI] [Google Scholar]
  • 10.Lindenmayer DB, et al. Conservation conundrums and the challenges of managing unexplained declines of multiple species. Biol. Conserv. 2018;221:279–292. doi: 10.1016/j.biocon.2018.03.007. [DOI] [Google Scholar]
  • 11.McLean CM, Kavanagh RP, Penman T, Bradstock R. The threatened status of the hollow dependent arboreal marsupial, the greater glider (Petauroides volans), can be explained by impacts from wildfire and selective logging. For. Ecol. Manag. 2018;415–416:19–25. doi: 10.1016/j.foreco.2018.01.048. [DOI] [Google Scholar]
  • 12.Smith P, Smith J. Decline of the greater glider (Petauroides volans) in the lower Blue Mountains New South Wales. Aust. J. Zool. 2018;66:103–114. doi: 10.1071/ZO18021. [DOI] [Google Scholar]
  • 13.Australian Government Department of Agriculture, Water and the Environment Threatended Species Scientific Committee. Conservation Advice Petauroides volans, greater glider. (Australian Government, 2016).
  • 14.Burbidge, A. A. & Woinarski, J. Petauroides volans (amended version of 2016 assessment). The IUCN Red List of Threatened Species. 10.2305/IUCN.UK.2020-1.RLTS.T40579A166500472.en (2020)
  • 15.McKay GM. Greater glider. In: Van Dyck S, Strahan R, editors. The Mammals of Australia. 3. Sydney: Reed New Holland; 2008. pp. 240–242. [Google Scholar]
  • 16.Comport S, Ward S, Foley WJ. Home ranges, time budgets and food-tree use in a high-density tropical population of greater gliders, Petauroides volans minor (Pseudocheiridae : Marsupialia) Wildl. Res. 1996;23:401–419. doi: 10.1071/WR9960401. [DOI] [Google Scholar]
  • 17.Lindenmayer D. Gliders of Australia: A Natural History. Randwick: UNSW Press; 2002. [Google Scholar]
  • 18.Collett R. On a collection of mammals from central and northern Queensland. Zool. Jahrbücher. 1887;2:829–940. [Google Scholar]
  • 19.Jackson SM. Family Pseudocheiridae (Ringtail Possums and Greater Gliders) In: Mittermeier RA, Wilson DE, editors. Handbook of the Mammals of the World. Volume 5: Monotremes and Marsupials. Barcelona: Springer; 2015. pp. 654–658. [Google Scholar]
  • 20.Jackson SM, Groves CP. Taxonomy of Australian Mammals. Clayton: CSIRO Publishing; 2015. [Google Scholar]
  • 21.Thomas O. On some Queensland Phalangeridae. Ann. Mag. Nat. Hist. 1923;11:246–250. doi: 10.1080/00222932308632848. [DOI] [Google Scholar]
  • 22.Jackson SM, Thorington RW., Jr . Gliding Mammals: Taxonomy of Living and Extinct Species. Washington: Smithsonian Institutions Scholarly Press; 2012. [Google Scholar]
  • 23.Arbogast, B., Aplin, K., Armstrong, K., Ogawa, L. & Taylor, A. Cryptic diversity in an iconic Australian marsupial: the greater glider, Petauroides volans. In Joint Meeting of the 57th Scientific Meeting of the Australian Mammal Society with the 91st Annual Meeting of the American Society of Mammalogists. Portland State University, Oregon. USA. 24–28 June 2011 (2011).
  • 24.Jaccoud D. Diversity Arrays: a solid state technology for sequence information independent genotyping. Nucleic Acids Res. 2001;29:25. doi: 10.1093/nar/29.4.e25. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 25.Flannery TF. Possums of the World: A Monograph of the Phalangeroidea. Berkeley: Geo Productions; 1994. [Google Scholar]
  • 26.Georges A, et al. Genomewide SNP markers breathe new life into phylogeography and species delimitation for the problematic short-necked turtles (Chelidae: Emydura) of eastern Australia. Mol. Ecol. 2018;27:5195–5213. doi: 10.1111/mec.14925. [DOI] [PubMed] [Google Scholar]
  • 27.Georges A, Adams M. Electrophoretic delineation of species boundaries within the short-necked freshwater turtles of Australia (Testudines: Chelidae) Zool. J. Linn. Soc. 1996;118:241–260. doi: 10.1111/j.1096-3642.1996.tb01266.x. [DOI] [Google Scholar]
  • 28.Anderson EC, Thompson EA. A model-based method for identifying species hybrids using multilocus genetic data. Genetics. 2002;160:1217–1229. doi: 10.1093/genetics/160.3.1217. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 29.De Queiroz K. Species concepts and species delimitation. Syst. Biol. 2007;56:879–886. doi: 10.1080/10635150701701083. [DOI] [PubMed] [Google Scholar]
  • 30.Baker RJ, Bradley RD. Speciation in mammals and the genetic species concept. J. Mammal. 2006;87:643–662. doi: 10.1644/06-MAMM-F-038R2.1. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 31.Maloney KS. The Status of the Greater Glider “Petauroides volans” in the Illawarra Region. Sydney: University of Wollongong; 2007. [Google Scholar]
  • 32.Maloney KS, Harris J. Early natural history of the greater glider, Petauroides volans (Kerr, 1792) Linn. Soc. N. S. W. 2008;129:39–55. [Google Scholar]
  • 33.Chambers EA, Hillis DM. The multispecies coalescent over-splits species in the case of geographically widespread taxa. Syst. Biol. 2020;69:184–193. doi: 10.1093/sysbio/syz042. [DOI] [PubMed] [Google Scholar]
  • 34.Kearney MR, Wintle BA, Porter WP. Correlative and mechanistic models of species distribution provide congruent forecasts under climate change. Conserv. Lett. 2010;3:203–213. doi: 10.1111/j.1755-263X.2010.00097.x. [DOI] [Google Scholar]
  • 35.Lindenmayer DB, et al. Fire severity and landscape context effects on arboreal marsupials. Biol. Conserv. 2013;167:137–148. doi: 10.1016/j.biocon.2013.07.028. [DOI] [Google Scholar]
  • 36.Youngentob KN, Wood JT, Lindenmayer DB. The response of arboreal marsupials to landscape context over time: a large-scale fragmentation study revisited. J. Biogeogr. 2013;40:2082–2093. doi: 10.1111/jbi.12158. [DOI] [Google Scholar]
  • 37.Wickham H. ggplot2: Elegant Graphics for Data Analysis. R package. Berlin: Springer; 2016. [Google Scholar]
  • 38.Venables WN, Ripley BD. Modern Applied Statistics with S. Berlin: Springer; 2002. [Google Scholar]
  • 39.Sansaloni C, et al. Diversity Arrays Technology (DArT) and next-generation sequencing combined: genome-wide, high throughput, highly informative genotyping for molecular breeding of Eucalyptus. BMC Proc. 2011;5:P54–P54. doi: 10.1186/1753-6561-5-S7-P54. [DOI] [Google Scholar]
  • 40.Kilian A, et al. Diversity arrays technology: a generic genome profiling tehcnology on open platforms. Methods Mol. Biol. 2012;888:67–89. doi: 10.1007/978-1-61779-870-2_5. [DOI] [PubMed] [Google Scholar]
  • 41.Altshuler D, et al. An SNP map of the human genome generated by reduced representation shotgun sequencing. Nature. 2000;407:513–516. doi: 10.1038/35035083. [DOI] [PubMed] [Google Scholar]
  • 42.Baird NA, et al. Rapid SNP discovery and genetic mapping using sequenced RAD markers. PLoS ONE. 2008;3:e3376. doi: 10.1371/journal.pone.0003376. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 43.Elshire RJ, et al. A robust, simple genotyping-by-sequencing (GBS) approach for high diversity species. PLoS ONE. 2011;6:e19379. doi: 10.1371/journal.pone.0019379. [DOI] [PMC free article] [PubMed] [Google Scholar]
  • 44.Unmack PJ, et al. Phylogeography and species delimitation of Cherax destructor (Decapoda: Parastacidae) using genome-wide SNPs. Mar. Freshw. Res. 2019;70:857–869. doi: 10.1071/MF18347. [DOI] [Google Scholar]
  • 45.Jombart T. Adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics. 2008;24:1403–1405. doi: 10.1093/bioinformatics/btn129. [DOI] [PubMed] [Google Scholar]
  • 46.Gruber B, Unmack PJ, Berry OF, Georges A. dartr: an R package to facilitate analysis of SNP data generated from reduced representation genome sequencing. Mol. Ecol. Resour. 2018;18:691–699. doi: 10.1111/1755-0998.12745. [DOI] [PubMed] [Google Scholar]
  • 47.Gower JC. Some distance properties of latent root and vector methods used in multivariate analysis. Biometrika. 1966;53:325–338. doi: 10.1093/biomet/53.3-4.325. [DOI] [Google Scholar]
  • 48.Cattell RB. The scree test for the number of factors. Multivar. Behav. Res. 1966;1:245–276. doi: 10.1207/s15327906mbr0102_10. [DOI] [PubMed] [Google Scholar]
  • 49.Commonwealth of Australia, Australia - Species of National Environmental Significance (Public Grids). Environment Australia https://www.environment.gov.au/fed/catalog/search/resource/details.page?uuid=%7B337B05B6-254E-47AD-A701-C55D9A0435EA%7D (2020).
  • 50.Commonwealth of Australia (Geoscience Australia). GEODATA COAST 100K 2004 data set. Geoscience Australia https://pid.geoscience.gov.au/dataset/ga/61395, https://researchdata.edu.au/geodata-coast-100k-2004 (2004).
  • 51. Santos Baquero, O. ggsn: North symbols and scale bars for maps created with ‘ggplot2’ or ‘ggmap’. R package (2019).
  • 52.Pebesma, E. J. & Bivand, R. S. Classes and methods for spatial data in R. R News 5, 9–13. R package (2005).
  • 53.Bivand RS, Pebesma E, Gomez-Rubio V. Applied Spatial Data Analysis with R. R package. 2. Berlin: Springer; 2013. [Google Scholar]
  • 54.Bivand, R. et al. rgdal: Bindings for the ‘Geospatial’ Data Abstraction Library. R package (2020).
  • 55.Hijmans, R. J. et al. raster: Geographic Data Analysis and Modelling. R package (2020).

Associated Data

This section collects any data citations, data availability statements, or supplementary materials included in this article.

Supplementary Materials

Data Availability Statement

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://research.jcu.edu.au/researchdata) at https://doi.org/10.25903/xx67-2m46.