Reassessment of the subspecific status of the Australian Wet Tropics yellow-bellied glider, Petaurus australis

Keywords: aspidonym, Endangered species, Evolutionarily Significant Unit, marsupial, phylogenetics, subspecies, systematics, taxonomy.

Petaurid gliders are widespread in their distribution across native forest regions of eastern and northern Australia (and northward to New Guinea and adjacent islands), but they have suffered declines in population size due to habitat clearance for agriculture and forestry, resulting in habitat fragmentation that is reducing migration and gene flow across the landscape (Goldingay et al. 2013; Malekian et al. 2015). Their distribution was also impacted historically through Pleistocene climatic changes, which led to substructuring of populations following the replacement of wet forest by open sclerophyll woodlands (e.g. Petaurus norfolcensis (Kerr, 1792); Pavlova et al. 2010).

The yellow-bellied glider, Petaurus australis Shaw, 1791, occurs in only a single isolated population in South Australia and is patchily distributed from Victoria through to northern Queensland. It is listed as Endangered or Vulnerable in three of the four states in which it occurs and hence there is a need for conservation management programs to be implemented. Such programs, which may require translocation of animals among populations to boost genetic diversity, should be based on a sound taxonomy of the species that takes into account the genetic distinctiveness of populations and their history of connectivity (Weeks et al. 2011).

The Wet Tropics (WT) population of the yellow-bellied glider in northern Queensland, referred to as Petaurus australis unnamed subspecies, is the most northerly population of the species in Australia, and is isolated from more southerly populations on the Clarke Ranges by the dry Burdekin Gap, spanning a distance of approximately 400 km (Fig. 1; Brown et al. 2006; Van Dyck and Strahan 2008; Jackson and Groves 2015). The WT population is further subdivided into three subpopulations – on the Mt Windsor Tableland, Mt Carbine Tableland and the Greater Cardwell Range (GCR; Fig. 2) – and is restricted to a narrow band of tall eucalypt forest on the western margins of the rainforest (Winter et al. 2004; Williams 2006). The three WT subpopulations are separated by natural habitat breaks of dry eucalypt woodland not used by the gliders. The Mt Windsor Tableland and Mt Carbine Tableland populations are separated by a gap of 15 km and a saddle of 750 m in elevation of extremely dry woodland. Seventy five km separates the Mt Carbine Tableland and the GCR subpopulations, with only a narrow band of suitable habitat on the western slopes of the Lamb Range that is not known to be occupied by the gliders. The GCR sub-population is further subdivided by natural breaks in the narrow band of the glider’s natural habitat, exacerbated by recent clearing for agriculture, particularly at the northern end in the Wondecla region.

Fig. 1.  Distribution of Petaurus australis and the Wet Tropics (WT) population, modified from Goldingay (2008) with permission from the author. Also shown is the location of the Burdekin gap.

Fig. 2.  Distribution of Petaurus australis within the Wet Tropics bioregion showing the three sub-populations, the Mt Spurgeon location of the type specimen and the Wondecla location of P. australis samples used in the current study. Imagery from Google Earth.

The WT population was first listed as a subspecies by Russell (1983) and allocated the subspecific name reginae Thomas, 1923. McKay (1988) also recognised P. a. reginae as a distinct subspecies comprising the WT population, and Flannery (1994) explicitly restricted the name reginae to the WT portion of the range, recognising it as a distinct WT endemic subspecies, P. a. reginae. However, this name was first applied (based largely on fur colour differences) to populations, not from the Wet Tropics, but from southern Queensland, with the type specimen coming from Gin Gin near Bundaberg (Thomas 1923), south of the Burdekin Gap (Fig. 1). Fur colour was later shown to vary with the age of individuals (Goldingay and Kavanagh 1990; Goldingay et al. 2001), darkening as the animal matures, initially casting doubt on the subspecific status of reginae and the need for the use of trinomial designations in P. australis. Nevertheless, it was proposed by Maxwell et al. (1996) that the WT population be called an ‘unnamed subspecies’ and this proposal was recognised as such by Groves (2005) and Clayton et al. (2006). Brown et al. (2006) carried out molecular (mitochondrial DNA, mtDNA) and morphological analyses of P. australis and concluded that the subspecies reginae was no longer valid, as the southern Queensland population was genetically more closely related to P. australis australis from NSW, Vic and SA compared to individuals from WT. Individuals from the latter population formed a reciprocally monophyletic group relative to a second group comprising individuals from all other populations, from southern Queensland to SA (Brown et al. 2006). Sequence data from one mtDNA locus, NADH Dehydrogenase subunit 4 (ND4) gene) from four animals, including just one sample from the GCR region, representing the largest sub-population in WT, was considered insufficient to delineate the WT population as a distinct subspecies, though the authors recommended it be treated as a distinct Evolutionarily Significant Unit (ESU) for the purposes of conservation management. The WT population is currently listed as Endangered under the Queensland Nature Conservation Act 1992 and Endangered under the Commonwealth Environment Protection and Biodiversity Conservation Act 1999, both of which recognise its subspecific status, albeit unnamed.

In the current study, we obtained additional samples from the GCR subpopulation of North Queensland to supplement those obtained by Brown et al. (2006) and carried out phylogenetic analyses of mtDNA and nuclear gene markers to further assess the subspecific status of the WT population. We also examined patterns of morphological variation in museum specimens of Petaurus australis collected throughout the range of the species. We employed criteria for defining subspecies using the concept first proposed by Avise and Ball (1990) and O’Brien and Mayr (1991) who defined subspecies as ‘groupings of populations within species that share a unique geographic range or habitat and are distinguishable from other subdivisions of the species by multiple, independent genetically based traits’ (as stated in Frankham et al. 2002, p. 617).

DNA samples (n = 30) for P. australis were sourced from the South Australian Museum (from the study of Brown et al. 2006), with the exception of six new glider samples that were collected from the Wondecla area of the GCR in north Queensland WT (Table 1, Fig. 2). The latter gliders were captured at night by hand by RR when they were feeding on sap at cuts made by them on the trunks of Eucalyptus resinifera. They were transferred to a cloth bag on-site and the tip of each ear removed by JW, with a new sharp disposable scalpel blade. Tissue samples were stored in tubes of 100% ethanol. To reduce stress to the animal it was held for the minimum amount of time, no more than 5 min, before being released back onto the tree. Consequently, no body measurements were taken from the GCR population. The animals appeared little concerned by their handling and were regularly seen the following night on the same tree, showing little interest in the capturer and coming again within capture distance. They were also seen weeks later, recognised by the missing ear tips. These six animals were members of two family groups occupying adjacent territories. The tissue samples were sent to the South Australian Museum for genetic analyses. The ear biopsies were collected under ethics and scientific cover of the Threatened Species Unit of the Queensland Department of Environment and Heritage Protection (SA2013/07/432).

Table 1.  Yellow-bellied glider and outgroup specimen numbers, voucher numbers and location (latitude, longitude) from which specimens originated (if known) and GenBank accession numbers for haplotypes.

DNA was extracted from skin or liver tissue using the Gentra Puregene extraction kit and methods specified by the manufacturer (Gentra Systems Inc.).

Our previous molecular analyses of P. australis (Brown et al. 2006) used a ~900-bp segment of the ND4 gene, therefore we also used this marker in the current study. ND4 was amplified using the primers M245 (5′-TGA CTA CCA AAA GCT CAT GTA GAA GC-3′) and M246 (5′-TTT TAC TTG GAT TTG CAC CA-3′), each based on primers ND4 and Leu designed by Arevalo et al. (1994). We also sequenced a portion of the NADH dehydrogenase subunit 2 gene (ND2) using the primers M635 (5′-GCA CCA TTC CAC TTY TGA GT-3′) and M636 (5′-GAT TTG CGT TCG AAT GTA GCA AG-3′) (Osborne and Christidis 2001). Two nuclear genes were also sequenced: a 945 bp fragment of the vonWillebrand factor gene (vWF exon 28), using primers G807 (5′-GAC TTG GCY TTY CTS YTG GAT GG-3′) and G2526 (5′-TTG ATC TCA TCS GTR GCR GGA TTG C-3′ (Amrine-Madsen et al. 2003), and a 700 bp fragment of the ω-globin gene using primers G314 (5′-GGA ATC ATG GCA AGA AGG TG-3′) and G424 (5′-CCG GAG GTG TTY AGT GGT ATT TTC-3′) (Wheeler et al. 2001). Each of these genes have been used in previous phylogenetic studies of petaurids (Malekian et al. 2010a, 2010b).

PCR amplifications were carried out in 25 μL volumes containing 0.1U AmpliTaq Gold^® polymerase (Applied Biosystems), 1 × AmpliTaq Gold Buffer, 0.20 mM dNTPs, 2.5 mM MgCl₂, 0.5 μM of each primer and approximately 100 ng genomic DNA. Thermocycling conditions were: initial activation at 94°C for 3 min; 35 cycles of denaturation at 94°C for 30 s, annealing at 48–55°C for 45 s, and extension at 72°C for 60 s; and a final extension at 72°C for 3 min. PCR products were purified using Millipore MultiScreen PCR 384 Filter Plates (Millipore) and were capillary sequenced by the Australian Genome Research Facility (AGRF) using the ABI Prism Big Dye Terminator Cycle Sequencing Kit (Applied Biosystems). Sequencing was carried out on an ABI 3700 DNA analyser and edited using SEQED ver. 1.0.3 (Applied Biosystems). Sequences were submitted to GenBank (Table 1).

DNA sequences were edited and aligned using the Geneious alignment option within Geneious (version 9.1.2; www.geneious.com). Pairwise distances among mitochondrial haplotypes were determined using Geneious using the HKY-85 (Hasegawa et al. 1985) model of sequence evolution. Geneious was also used to construct Neighbour Joining (NJ) trees using the same model. The robustness of nodes in the NJ trees was assessed by 1000 bootstrap replicates.

The ND2 data were also phylogenetically analysed using Maximum likelihood (ML) as implemented in the program RAxML (version 7.2.8; Stamatakis et al. 2008) provided as a plugin for the program Geneious (version 9.1.2). A single model of evolution, General Time Reversible (GTR) model (Rodríguez et al. 1990) with unequal variation at sites modelled using a Gamma (G) distribution (Yang 1996) was applied to the sequence data. Robustness of branches on the tree was assessed using 500 bootstrap pseudoreplicates. Trees were visualised and rooted using the petaurids Dactylopsila trivirgata and Gymnobelideus leadbeateri as outgroups for analyses including Petaurus species or using a midpoint rooting approach in FigTree (version 1.4.2; http://tree.bio.ed.ac.uk/), for analyses of P. australis ND4 data only.

The vWF sequence data showed low levels of divergence (i.e. haplotypes differing by 1–2 nucleotide sites) within P. australis and, therefore, a haplotype network approach was used to visualise relationships among haplotypes from P. australis populations. Ambiguities in vWF sequences, representing polymorphic variation, were resolved into distinct haplotypes for each individual and an unrooted NJ network, based on HKY-85 distances (Hasegawa et al. 1985) among haplotypes, was derived using Geneious. This network was exported as a NEWICK tree and then converted into a haplotype network using the program Haploviewer (developed by G. Ewing, http://www.cibiv.at/~greg/haploviewer).

An assessment of how genetic variation was partitioned among groups of P. australis populations was determined using Analyses of Molecular Variation (AMOVA), using the program Arlequin (v. 3.5.2.2; Excoffier and Lischer 2010), with default options based on the frequency of pairwise differences (p-distances) among haplotypes in populations. This approach was used to compare different hypotheses of subspecies boundaries for both ND4 and vWF data.

To estimate the coalescent time of P. australis mtDNA haplotypes, molecular clock analyses were conducted on the combined ND2–ND4 sequence data using BEAST version 2.4.7 (Bouckaert et al. 2014). Exemplar sequences from Australian Petaurus species, the New Guinean P. abidi and the outgroup taxon Petauroides volans, were included in the analyses. First, second and third codon positions of each gene, were analysed as separate partitions using the program PartitionFinder (version 1.1.1; Lanfear et al. 2012) to find an appropriate partitioning scheme and models of nucleotide evolution for each partition in the BEAST analyses. The best partitioning scheme based on the Bayesian Information Criterion suggested combining first, second and third codon positions of each gene, resulting in three partitions for phylogenetic analysis: ND2–ND4 first, ND2–ND4 second, ND2–ND4 third, with the models HKY+G, HKY+I and TrN+I (Hasegawa et al. 1985; Tamura and Nei 1993; Yang 1996), respectively. The program BEAUti was used to define the site substitution models for each partition and assign molecular clock models and tree and parameter priors for different analyses. To calibrate the molecular clock a divergence time of 4.46 ± 0.1 Ma (mega annum) for the appearance of Petaurus in the fossil record was used as a minimum bound. These fossils were a mixed sample of teeth from deposits at Hamilton Victoria that were considered by Turnbull et al. (2003) to be near P. australis and P. norfolcensis. Hence, the 4.46 Ma date was placed on the common ancestral node of these two species = common ancestor of all the known extant Petaurus species (Malekian et al. 2010b). A lognormal distribution with an offset of 4.46 Ma, a mean (m) = 1.0 and a standard deviation (s) = 1.25 was used as the prior distribution for this node calibration. A birth-death tree prior, using estimated birth and death rates and uniform priors for each parameter was used in all the analyses. Molecular clock models and trees were linked across each partition. To estimate the time to most recent common ancestor (tmrca) of the P. australis root a uniform prior ranging between 0.0 and 4.46 Ma was placed on this node. Finally, two different molecular clock models were used for comparison: a strict clock model and an uncorrelated lognormal model; the prior distribution of the molecular clock rate under a strict clock model was given a Gamma distribution with alpha = 0.001 and beta = 1000; the prior distribution of the uncorrelated lognormal model used an exponential prior distribution with m = 10.0 for the mean and a Gamma prior distribution with alpha = 0.5396 and beta = 0.3819 for the standard deviation.

BEAST analyses were run for 15 million generations in multiple independent runs to test for convergence, sampling parameters and trees every 1500 generations, using the CIPRES Science Gateway (Miller et al. 2010). The program Tracer 1.6 (Rambaut et al. 2013) was used to evaluate convergence of parameter estimates following a burnin of 10%. Treeannotator (version 2.4.7) was used to generate a maximum clade credibility tree of the 15 000 trees that were sampled, using a burnin of 10% (1500 trees). Trees were visualised and initially prepared for publication using the program FigTree (version 1.4.2).

Museum specimens discussed here are deposited in the collections of the Australian Museum, Sydney (AMS), the American Museum of Natural History, New York (AMNH), the Natural History Museum, London (NHMUK), the Australian National Wildlife Collection, Canberra (ANWC), the Queensland Museum, Brisbane (QM), the Museum of Comparative Zoology at Harvard University, Cambridge, Massachusetts (MCZ), Museum Victoria, Melbourne (NMV), and the South Australian Museum (SAMA). Craniodental measurements were taken with handheld digital calipers graduated to the nearest 0.01 mm (Supplementary Table S1). Measurements and statistical comparisons are based on adult specimens as judged by direct examinations of dental and sutural maturity in skulls. Measurements taken include the greatest length of skull (GLS), condylobasal length (CBL), zygomatic width (ZYG), length of upper molars (M1-4), width of M1 (WM1), maximum width of palate (MWP), maximum breadth of braincase (BBC), height of braincase from the basioccipital floor to the apex of the braincase (HBC), nasal length (NL), nasal maximum width (NW), length of the lower incisor (length inc) and depth of the lower incisor (depth inc) (Supplementary Table S1). The dental formula of each specimen was also examined, following Archer (1984).

Body mass measurements of the southern (latitudes −35.5° to −37.5°; Victoria and southern New South Wales), central (latitudes −25° to −28°; southern Queensland) and northern populations (latitudes −17.4° to −17.9°; north Queensland’s Wet Tropics) were used to perform a box plot analysis and two-way and one-way ANOVAs, comparing each population and sexes within populations. Measurements taken on skulls were analysed using a principal component analysis (PCA). The analysis was performed both on raw measurements and on transformed measurements, after the effects of allometry were removed, using the software PAST (Hammer et al. 2001). Because several skulls were unsexed, we performed the analyses with and without defining the sex, to detect whether they helped or muddled the results. A MANOVA was also performed on each data set to identify if the means were significantly different from one another.

Specimens from the AMNH and MCZ, representing a larger sample of the northern population, did not have all measurements taken (only two) and couldn’t be used in the PCA. To include them and have a larger sample size, we produced an X–Y graph using those two measurements, using both raw and transformed measurements (allometry effects removed), and including and excluding unsexed skulls.

Six animals were caught, ear biopsies taken (GCR1-6), and released, all from the Wondecla area within the GCR subpopulation in WT (Fig. 2). These six samples were PCR-amplified and sequenced for ω-globin (750 bp), ND4 (884 bp), ND2 (696 bp) and vWF (945 bp). The sequences of ω-globin from the GCR specimens were found to be identical to each other and to ω-globin sequences from P. australis from western Victoria (R9 and R10). Due to a lack of variation this marker was dropped from further analyses.

The ND2 sequences from the six GCR specimens were found to be identical, and most closely related to an ND2 haplotype from a specimen (80840) from Mt Windsor National Park, WT (p-distance = 0.077%) and 2.2–2.5% divergent from P. australis from western Victoria and NSW (76609 and 75336). Similarly, the ND4 sequences were identical for the GCR specimens, 0.6% divergent from 80840 from WT and between 1.6 and 1.7% divergent from P. australis from western Victoria and NSW (76609 and 75366) respectively. Phylogenetic analyses of individual genes, ND2 and ND4, produced identical gene trees, providing evidence that these sequences were derived from the same genetic locus and, therefore, most likely mitochondrial in origin (results not shown). Phylogenetic analyses of concatenated ND2 and ND4 data using ML showed strong support (bootstrap support (BS) = 99%) for a monophyletic group comprising all the WT samples to the exclusion of P. australis from NSW and Victoria (Fig. 3). However, the position of the root into the P. australis clade was not well supported (BS = 65), possibly due to the long branch connecting the outgroup taxa (Petauroides volans (Kerr, 1792), P. abidi, P. norfolcensis and P. breviceps Waterhouse, 1839) to the ingroup. These analyses further showed that P. australis is a sister lineage to the other species in the Petaurus genus, an arrangement that was strongly supported (BS = 100%).

Fig. 3.  Maximum likelihood tree based on combined ND2 and ND4 data from P. australis and other species of Petaurus. The tree is rooted using Petauroides volans as an outgroup. Numbers adjacent to branches represent bootstrap values from 500 pseudoreplicates.

Using the new ND4 data from WT, together with ND4 data (n = 30) from P. australis derived by Brown et al. (2006), we further analysed the phylogeographic structure evident within the species using ML analyses (Fig. 4). These analyses showed the presence of two divergent monophyletic groups of ND4 haplotypes, one comprising samples from WT and a second group comprising samples from SQ, NSW, Victoria and SA. Each group was supported by a BS value of 100%. A NJ tree was also converted into a haplotype network to provide an alternative view of the ND4 haplotype relationships, showing the presence of three distinct haplotypes in the WT population that were at least 13 mutational steps away from the most closely related haplotype from a SQ location (Fig. 5a).

Fig. 4.  Maximum likelihood tree based on ND4 data from P. australis and showing reciprocal monophyly of a WT group relative to a group comprising P. australis from SQLD, NSW, Vic and SA. The tree was rooted using a mid-point root. Numbers adjacent to branches represent bootstrap values from 500 pseudoreplicates.

Fig. 5.  Haplotype networks for ND4 (a) and vWF (b), with circles arbitrarily colour coded by region (WT = northern QLD, SQ = southern QLD).

AMOVA analyses of the ND4 data revealed significant F_ST values (P = 0.000) among sub-populations within a two population structure model and that a maximum of variation among groups of 53.65% was obtained when WT was specified as a distinct population relative to a population comprising samples from SQ, NSW, Vic and SA. All other combinations of sub-populations under a two population structure model revealed lower levels for the partitioning of among group genetic variation (<29%). The optimum population structure model, however, comprised three populations: (WT), (SQ, NSW) and (Vic, SA), which explained 76% of among group genetic variation. This genetic variation among SQ-NSW and Vic-SA can largely be attributed to the presence of a single fixed haplotype in the latter population which is at the extreme western end of the distribution of P. australis.

A 945 bp segment of the vWF gene was sequenced from 25 specimens, with 6 specimens showing heterozygosity for alternative haplotypes. These latter sequences showed variation at a single nucleotide site, and hence could be unambiguously resolved into distinct haplotypes, giving a dataset of 31 sequences that were used to construct a haplotype network (Fig. 5b). The network shows the presence of six haplotypes of which two were found exclusively in WT specimens only, and a third haplotype (n = 8) that was found in seven WT specimens and one specimen from NSW.

AMOVA analyses of the vWF data also revealed significant F_ST values among sub-populations within groups for a two population structure model. For a structure comprising WT versus SQ, NSW, Vic/SA (i.e. the subspecies hypothesis), 29.5% of genetic variation was partitioned among the groups, but the optimum model of population structure comprised three populations: (WT), (SQ, NSW) and (Vic, SA), which explained 52.2% of among group genetic variation.

The ND2–ND4 data set was further analysed using molecular clock Bayesian analyses to estimate the coalescent time of mtDNA haplotypes from the WT population relative to those from NSW and Vic. Initial analyses using the Partition Finder substitution models (HKY+G, HKY+I and TrN+I for 1st, 2nd and 3rd codon positions respectively) failed to reach convergence for the parameter of the invariant sites model despite multiple attempts to tune this parameter. Therefore, a Gamma model of rate variation (Yang 1996) was used for each of the partitions to model site specific rate variation. Estimated sample size (ESS) values for all parameters were >312, suggesting an adequate sample of the posterior distribution for each parameter had been obtained. Under a strict molecular clock model the tmrca estimate for P. australis was 169.3 thousand years (95% posterior distribution (PD) ranged from 102.3 to 248.6 thousand years). Under an uncorrelated lognormal relaxed clock model the estimate for the P. australis MRCA was 137.5 thousand years (95% PD from 76.4 to 207.6 thousand years; Supplementary Fig. S1).

Skulls were examined for morphological differences. Skulls in the northern WT population were overall smaller with a shorter snout, compared to the southern population. The dental formula for the upper dentition (I1-3, C1, P1-3, M1-4) didn’t differ from that identified for P. australis by Archer (1984). The lower dentition was, however, different between the northern and southern populations. The southern population had the dental formula reported by Archer (1984), i1-2, p1-3, m1-4, but the northern population had variable numbers of unicuspid teeth (i2-p3), as high as four, as in the southern population, and as low as two.

We compared weights of wild-caught P. australis from across their distribution (see Fig. 6 caption for data sources) and found that they varied significantly between southern and northern populations (two-way ANOVA, F_1,147 = 67.4, P < 0.0001) and between genders (F_1,147 = 9.1, P = 0.003), with the southern population and males being heavier. However, the interaction between population and gender was also significant (F_1,147 = 8.7, P = 0.004), indicating that the relationship between the sexes differs between northern and southern populations (Fig. 6). In the southern population, males are markedly heavier than females (mean 601 cf. 543 g, n 20 and 16), whereas in the northern population the difference is slight (mean 493.2 cf. 492.5 g, n 58 and 55) and not statistically significant (one-way ANOVA, F_1,112 = 0.1, P = 0.94).

Fig. 6.  Weights of yellow-bellied gliders (Petaurus australis) by population and gender. The southern population is from latitudes −35.5° to −37.5° (Victoria and southern New South Wales), the three central population individuals are from latitudes −25° to −28° (southern Queensland) and the northern population are from latitudes −17.4° to −17.9° (north Queensland’s Wet Tropics). Box plots show mean (heavy line), median (light line), 25–75th percentiles (box), 10–90th percentiles (whiskers) and outliers (dots). Data sources: southern Craig (1985), Henry (1985), Goldingay (1992), Goldingay and Kavanagh (1990); central Queensland Museum; northern Goldingay et al. (2001), Quin pers. comm. 2018, R. Russell unpubl. data.

The results of the Principal Component Analysis (PCA) of cranial measurements, show a significant amount of overlap between the southern and northern populations (Fig. 7). Using all the raw measurement available (Fig. 7a), including from unsexed skulls, specimens from the northern population are within the range of the southern population. One unsexed specimen (SAMA M2237) is large and fits perfectly within the measurements of the southern population. The MANOVA found that the two populations were not significantly different (F = 1.725; P = 0.09868). If unsexed specimens are removed (Fig. 7b), the northern population overlaps less so with the southern population, occupying the top left corner of the PCA with negative component 1 values, while the majority of the southern population had positive values on component 1. The MANOVA found that the northern and southern populations have significantly different means (F = 2.997, P = 0.01148). After transforming the measurements by removing the effects of allometry (option selected in PAST), using all specimens (Fig. 7c), the northern population only overlaps slightly with the southern population, showing that the difference between the two populations is not only driven by size, but by selection in the shape of the skull. The MANOVA found a significant difference between the two populations (F = 2.915; P = 0.005693). When unsexed specimens were removed (Fig. 7d), the same results were found, corroborated by the results of the MANOVA (F = 3.222, P = 0.007707).

Fig. 7.  Principal Component Analysis (PCA) of the cranial measurements, with convex hulls. (a) using raw measurements including unsexed skulls, (b) using raw measurements, excluding unsexed skulls; (c) using transformed measurements (allometry removed), including unsexed skulls; (d) using transformed measurements, excluding unsexed skulls. Blue = northern population; Red = southern population. Triangles = males; circles = females; square = unsexed skulls.

Condylobasal length (CBL) is plotted against nasal length (NL) in Fig. 8. The sample size was large enough to be able to provide 95% ellipses, with any samples outside of the ellipses considered an outlier or as being distinct. Using all specimens’ raw measurements (Fig. 8a), most of the specimens of the northern population fall outside the 95% ellipse of the southern population, suggesting that they are distinct. SAMA M2237, which was found previously to be too large to be from the northern population, is just at the edge of the northern population 95% ellipse. Two specimens from the Southern population are found outside of the 95% ellipse, SAMA M12590 (unsexed) and NMV C2396.14 (female), suggesting that they are outliers or distinct. When unsexed specimens are removed (Fig. 8b), there is much less overlap between the two 95% ellipses. NMV C2396.14 (female) is again found outside the ellipse for the southern population, as well as AMS M357 (male). When the measurements are transformed to remove the effects of allometry (Fig. 8c), the two 95% ellipses overlap each other significantly, with only a couple of specimens from the northern population outside of the southern population ellipse. SAMA M12590 (unsexed) and NMV C2396.14 (female) are found again outside of the southern population ellipse, though SAMA M12590 (unsexed) much more than NMV C2396.14 (female). When unsexed individuals are removed (Fig. 8d) the two ellipses overlap much less, with half of the northern population specimens outside of the southern population ellipse. NMV C2396.14 (female) and AMS M357 (male) remain the most distinct within the southern population, though they are within the 95% ellipse.

Fig. 8.  X–Y graph of condylobasal length (CBL) versus nasal length (NL), with 95% ellipses. (a) using raw measurements including unsexed skulls, (b) using raw measurements, excluding unsexed skulls; (c) using transformed measurements (allometry removed), including unsexed skulls; (d) using transformed measurements, excluding unsexed skulls. Blue = northern population; Red = southern population. Triangles = males; circles = females; square = unsexed skulls.

Our phylogenetic and haplotype network analyses based on mtDNA and nuclear gene loci, when combined with the morphometric analyses that show that the WT specimens had a shortening of the rostrum compared to southern specimens, provide strong support for the distinctiveness of the WT population of P. australis. We further discuss these analyses below, provide additional support from ecological studies for the distinctiveness of the WT population, and propose that it should be recognised as a distinct subspecies that we describe below.

Mitochondrial and nuclear gene data from across the range of P. australis in Australia support the presence of two major genetic lineages, one distributed in WT and a second distributed from southern Queensland to SA (Figs 3, 4). These analyses further support the conclusion from Brown et al. (2006) that the WT population represents a distinct ESU, relative to an ESU comprising all other localities in southern QLD, NSW, Vic and SA, using criteria for defining ESUs proposed by Moritz (1994). Under these criteria each ESU is reciprocally monophyletic for mtDNA loci and shows significant allele frequency variation at nuclear DNA markers.

Based on the coalescent time for the mtDNA haplotypes, the two lineages are likely to have been isolated for >130 thousand years, indicative of genetic isolation over multiple ice age cycles. Although equating population divergence dates to the coalescent time of a single genetic locus is subject to error, they nevertheless suggest that the WT lineage has been evolving on an independent evolutionary trajectory compared to southern populations over a long time period that is significant within the evolutionary history of the species. The divergence among mtDNA haplotypes between WT and southern populations of P. australis (>2.2%) is similar to that found between P. norfolcensis and P. gracilis (De Vis, 1883) (>2.0%) and greater than that found at the intra-specific level within P. norfolcensis (Pavlova et al. 2010), which is distributed over a similar range. Genetic studies of the sugar glider P. breviceps reported high levels of intra-specific variation (10.4–12.2%; Malekian et al. 2010b), which has led to the recognition of two additional species (P. notatus and P. ariel) distinct from P. breviceps (Cremona et al. 2021).

Brown et al. (2006) concluded from their examination of morphological material over the full range of Petaurus australis that ‘morphological analysis provided little evidence for discrimination of populations, although NQ specimens were generally smaller in size than southern forms.’ Our results agree partly with Brown et al. (2006), that the northern population is generally smaller in size, but we have found that skull size was not the only factor contributing to their differentiation. Our morphometric analyses that excluded the effects of allometry showed a significant difference between the northern and southern population, meaning that the shape of the skulls was significantly different. This analysis suggests that facial shortening is being selected in the northern population, by the shortening of the nasal length, with evidence of loss of teeth in the lower jaw, with fewer teeth between the front incisors and the molars. The southern population had the dental formula reported by Archer (1984), i1-2, p1-3, m1-4, with some specimens having an extra premolar. The specimens from the WT population had either the same dental formula or had fewer unicuspid teeth between i1 and m1, with the absence of one or two of those teeth. Overall, our results support the WT population as being distinct morphologically from the southern population.

Some specimens do not follow the trend observed in other specimens from the same region. For example, SAMA M2237 from Herberton (Queensland; north-east of Wondecla, see Fig. 2) is very large and fits better with the morphology seen in the southern population. It is possible that the skin and skulls of the specimens were mixed up, and that the skull belongs to a different specimen, perhaps SAMA M2745 from Lower Fitzroy Valley, which has a much smaller skull than what would be expected from that region. NVM C2396.14 from the Otway Range in Victoria also came up as an outlier in the X–Y plot, mainly because it has the largest nasal length of any specimen. All other measurements are within the range of the population, and therefore this specimen isn’t seen here as being distinct morphologically, but rather an outlier in the dataset. SAMA M12590 is, however, a real mystery. It is unsexed and from an unknown location. The comments associated with the specimen say ‘Old museum stock’, and it was previously identified as Petaurus breviceps. Our results suggest that it is so distinct from Petaurus australis, that this specimen probably represents another taxon. Finally, AMS M357 is a male from Orange (NSW), which is very small in size and has a white belly. Both features would suggest that it is distinct from the rest of the southern population specimens and should warrant genetic assessment in the future.

General observations by two of the authors (R. R. and J. W.) of the gliders in the most northerly population isolates in the GCR, the Mt Carbine Tableland and Mt Windsor Tableland subpopulations, noted that the gliders have paler belly fur than the rich yellow common in more southerly populations (Russell 1984) (Fig. 9a). If this is the case, it could be further evidence of a break in gene flow, providing additional support for the taxonomic distinction of the WT population. However, fur colour is often a difficult character on which to base taxonomic separation of populations as it can vary within a population and if based on museum specimens it is subject to change owing to foxing (Wagstaffe and Williamson 1947). Goldingay et al. (2001) disputed the observation of the white belly fur in WT populations as they found in their Nitchaga study population, in the central section of the GCR, that 74% had yellow or dark yellow belly fur (Fig. 9b), though they did not compare the fur colour intensity of Nitchaga animals with those from southern Australia. In addition, the type specimen we are using from Mt Spurgeon, Mt Carbine Tableland subpopulation, collected in 1937, has strong yellow colouration, even allowing for the colour becoming darker on the specimen owing to foxing (Fig. 9c). It is from the same population which RR and JW have participated in seven biennial population surveys, with missing years, since 1997 without noting strong yellow belly colour of the adult gliders. Also, the first two specimens ever collected from the Wet Tropics region, in the Herberton area in 1912, were described by Finlayson (1934) as ‘larger, much darker and more richly yellow’ than ‘three from Rocky Waterhole, Serpentine Creek, about 28 miles north-east of Rockhampton.’ The Herberton area is the locality of Russell’s intensive behavioural study of the gliders (Russell 1984) and he never noted gliders with richly yellow bellies. This raises the question: is the Wet Tropics population trending towards an increasingly higher proportion of individuals with pale belly fur? Is it a case of genetic drift (Messer et al. 2016) occurring within relatively small population fragments of the Wet Tropics gliders? It is beyond the scope of our study to determine whether this is occurring, but, nevertheless, the white belly fur is a prominent feature of WT glider populations, further supporting the suggestion that they are on an independent evolutionary trajectory, compared to populations in the south.

Fig. 9.  Belly fur colouration of an adult P. australis from the Wet Tropics population: (a) adult from Wondecla, Herberton area in Queensland (photo R. Russell); (b) dark yellow belly of an individual from the Nitchaga population, Queensland (courtesy of Darren Quin, Queensland Electrical Commission 1991. Tully-Millstream hydro-electric scheme - Yellow-bellied glider baseline study - Volume 2. Appendix 2, Photograph 9); (c) dark yellow of the type specimen J6352 from Spurgeon population (courtesy of Queensland Museum, Dr John Sheridan).

The WT population of P. australis also exhibits differences in its socioecology from southern populations. These differences are not necessarily indicative of taxonomic differences, but likely reflect different responses to their environment. WT gliders form groups of between two and six individuals (Russell 1984; Goldingay et al. 2001). The social structure (e.g. genetic relatedness of individuals) of these groups is currently unknown. The structure of the WT groups varied, usually containing an adult male and female, but frequently with two adult females as part of the group (Russell 1984; Goldingay et al. 2001). In contrast, analyses of social structure in P. australis populations from Victoria have provided evidence for social and genetic monogamy, with single males and females sharing their home ranges and dens, and with social partners shown, using microsatellite DNA analyses, to be the true parents of juveniles in the study population (Brown et al. 2007). These differences in social behaviour may have resulted from variation in food resource abundance and productivity at the different localities (Goldingay 1992; Goldingay et al. 2001). For example, the availability of nectar and pollen food resources throughout the year may promote larger group sizes and a possible polygynous mating system (Goldingay 1992; Goldingay et al. 2001; Brown et al. 2007). Further research of social behaviour of the WT population and other populations of P. australis in eastern Australia is required, but taken overall, current observations suggest that the behavioural ecology of the WT gliders is likely to be distinct from that in southern populations.

The WT population of P. australis is unusual in that the gliders tap only one species of tree for its sap, Eucalyptus resinifera (Russell 1984; Quin et al. 1996). The only other population that targets only one species, Eucalyptus viminalis, is at the other extremity of the glider’s range, Rennick Forest in Victoria on the border with South Australia (Carthew et al. 1999). Carthew et al. (1999) point out that the restriction to one tree in Rennick Forest may be attributed to the paucity of eucalypt species present, just two with a third present in some locations. This contrasts with the greater range of species available to the WT population, both within their tall eucalypt forest habitat and adjoining rainforest (Quin et al. 1996). Tree species readily available within the WT glider habitat are E. tereticornis, E. grandis, E. moluccana and Corymbia intermedia, all of which have been recorded as tapped by the gliders in southern populations (Mackowski 1988; Eyre and Goldingay 2005). Although Eyre and Goldingay (2005) state that E. tereticornis is not a favoured tree, it is tapped for its sap by the most westerly population in Queensland in the Carnarvon Ranges together with grey gum and spotted gum (T. Eyre pers. comm. October 2018). In all other localities at least two species of tree are tapped for their sap by the gliders (Henry and Craig 1984; Craig 1985; Goldingay 1992; Goldingay and Kavanagh 1993; Eyre and Goldingay 2003). In addition, two species of tree, E. resinifera and E. moluccana, are tapped in the population immediately to the south of the WT population, at Crediton, west of Mackay (J. Winter pers. obs. 1989). There is one example in the Wet Tropics of an E. grandis tapped by the gliders. It was a young coppice trunk arising from a stump and heavily incised by the gliders in the Mt Carbine Tableland sub-population (R. Russell pers. obs. 1990s). It is unclear why the gliders did not use the readily available E. grandis more often within their range.

A definition of a subspecies was proposed by Lidicker (1962). He states that ‘a subspecies is a relatively homogeneous and genetically distinct portion of a species which represents a separately evolving, or recently evolved lineage, with its own evolutionary tendencies, inhabits a definite geographic area, is usually at least partially isolated, and may intergrade gradually, although over a fairly narrow zone, with adjacent subspecies’ (see Owen 1989). Similar definitions were proposed by Avise and Ball (1990) and O’Brien and Mayr (1991), which further emphasised that subspecies should show divergence for multiple independent genetically based traits.

The use of genetic criteria to define taxonomic units, particularly at the species level, has been proposed by multiple authors (e.g. Halt et al. 2009; Cook et al. 2010; Harvey et al. 2015; Westerman et al. 2016). Subspecies, however, are usually defined based on morphological traits, with genetic data used to define alternative conservation units such as ESUs where clear genetic criteria for their recognition have been proposed (e.g. Moritz 1994). However, we concur with the view of Aplin et al. (2015) that the use of ESUs is largely informal and this concept has often not been recognised in legislation for their conservation protection (e.g. the federal Department of Environment in Australia), in contrast to subspecies. The Moritz (1994) criteria for defining ESUs – reciprocal monophyly of mtDNA, reflecting long term population isolation, and allelic/haplotype frequency variation for nuclear gene markers – in addition to consideration of other criteria (geographic isolation or presence of narrow hybrid zones, ecological differences, morphological variation) provides a sound basis for formally defining subspecies. Overall, our study fulfils each of these criteria, with the WT population representing a distinct ESU, with long-term genetic and geographic isolation. It shows distinct morphological differences (shortening of the rostrum) compared to southern populations, smaller size, a trend towards a white belly colour and some ecological differences associated with its habitat and behaviour. We consider that this degree of divergence at multiple independent genetically based traits warrants its recognition as a distinct subspecies.

Since no reliable scientific name has previously been applied to P. australis from northern Queensland (McKay 1988), we here describe the WT population as a new subspecies.

Following the recommendations of Kaiser et al. (2013) and Jackson et al. (2022) and the position adopted by Taxonomy Australia (2021, downloaded from www.taxonomyaustralia.org.au/codes-of-conduct) we consider previous names proposed for yellow-bellied glider (Wet Tropics) to be unscientific and unavailable.

Aspidonym – Petaurus australis brevirostrum subsp. nova Cooper, Travouillon and Helgen

Northern yellow-bellied glider.

QM J6352, adult female, study skin (figured by Van Dyck 1997) and skull, from Mount Spurgeon, Queensland, collected 8 February 1938 by Gabriele Neuhäuser (see Mather (1986) and Van Dyck (1997) for biographical notes on Neuhäuser, and context of collection). The type locality is located today within the boundaries of Mount Spurgeon National Park.

There are additional specimens at AMNH (AMNH 107263–107268, collected by G. Neuhäuser in November and December 1937) and at MCZ (MCZ 29164–29168, collected by P. J. Darlington in July 1932) that were collected from Mount Spurgeon, the type locality (Tate 1952).

Six individuals sampled and released from the Wondecla area of the GCR, in addition to ABTC80841 from Brown et al. (2006) from the same region, which were included in our genetic comparisons, are referred to this new subspecies. In addition, ABTC80838–80840 from Brown et al. (2006) from Mt Windsor Tableland are referred to this new subspecies. DNA for these samples is available from the SA Museum DNA archive.

Three specimens in the Queensland Museum from Nitchaga Creek on the Atherton Tableland, Queensland, are also referred to this subspecies: QM JM8746, adult male, skull and postcranial skeleton, from Nichaga Creek, Atherton Tableland, collected 16 July 1991 by D. Quinn; QM JM8747, adult male, skull and postcranial skeleton, from Nitchaga Creek, Atherton Tableland, collected 29 April 1991 by D. Quin; QM JM8503, adult female, skull and partial postcranial skeleton, from Nitchaga Creek, Atherton Tableland, collected 10 October 1990 by D. Quin.

Two skins at SAMA are referred to this subspecies: SAMA M2236, adult skin (without any accompanying skull), and SAMA M2237, adult skin (with an accompanying skull, which may be mismatched, see above), both from Herberton, Queensland, collected by S. Grant in 1912. Both skins had their skulls removed from their study skins, presumably (based on information from the old specimen register) after arriving at the museum.

This subspecies is restricted to very tall eucalypt forest on the western edge of rainforest in the Wet Tropics Bioregion of north Queensland at upland altitudes above 600–700 m elevation, on tableland or plateau geomorphology. It occurs from the Herbert River Gorge to the Mount Windsor Tableland and is highly likely only found north of Burdekin Gap (Fig. 2). The subspecies Petaurus australis australis is only found south of Burdekin Gap.

Petaurus australis brevirostrum differs from P. a. australis in averaging smaller in skull length, and in having a proportionally shortened facial region of the skull, as reflected in its especially short nasals (Supplementary Table S1); this shortening of the face is also reflected in some specimens by having the number of lower unicuspid teeth (between i1 and m1) sometimes reduced to 2 or 3 on each side, compared to the 4 unicuspids (i2, p1-3) on each side of the lower jaw characterising P. a. australis. The two subspecies are similar in external appearance, but P. a. brevirostrum averages smaller in overall body size (Fig. 6) and has belly fur that is less distinctly yellow, often appearing white overall, compared to P. a. australis (Fig. 9).

The name brevirostrum means short snout and is derived from Latin ‘brevis’ for short and ‘rostrum’ meaning beak or snout. The name brevirostrum refers to the shortened rostrum of gliders from the Wet Tropics. The subspecific name is used as a noun in apposition.

Supplementary material is available online.

All sequence data are available on GenBank (see Table 1 for accession codes) and morphometrics data are provided as supplementary information.

The authors declare no conflicts of interest.

This research did not receive any specific funding, with genetic analyses funded by SJBC’s discretionary research funds in the SA Museum and morphological examinations in museums funded by KMH’s discretionary funds at the University of Adelaide.