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

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° 56′ 58″ E, n = 18) (Fig. 4) as this was the suggested location of the proposed P. armillatus^19,20. In addition, 12 museum specimen tissue samples were obtained from the Queensland Museum from greater gliders collected from northern, central and southern Queensland to investigate genetic structure in that area more broadly (Fig. 4). Additional information about the climate, geography and vegetation at each sample location can be found in Supplementary Material, (Table S1 online).

Wild greater gliders were located through spotlight searches using high-powered, handheld torches (Ledlenser P7, Zweibrüder Optoelectronics GmbH and Co., Solingen, Germany) to detect greater glider eye shine. Individuals were then captured using a gas-powered, tranquilizer dart-gun (Montech Black Wolf; Tranquil Arms 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 package³⁷. We then used a canonical variate analysis (CVA) in R package MASS³⁸ 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 sequencing³⁹, as described in detail by Kilian et al.⁴⁰ and Georges et al.²⁶. Essentially, DArTseq is an implementation of sequencing complexity-reduced representations⁴¹ and more recent applications on next generation sequencing platforms^42,43. To achieve the most appropriate complexity reduction (the fraction of the genome represented, controlling average read depth and number of polymorphic loci), four combinations of restriction enzymes (Pstl enzyme combined with either Hpall, Sphl, Nspl or Msel) were evaluated and restriction enzyme combination of Pstl (recognition sequence 5′-CTGCAIG-3′) and Sphl (5′-GCATGIC-3′) was selected.

Amplification using polymerase chain reaction (PCR)^26,44 and the conditions applied are as described in Georges et al.²⁶. 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.²⁶. 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}⁴⁵) to facilitate processing with package dartR (version 1.8)⁴⁶. We first removed all but one SNP from each sequence tag (12,782 SNPs removed) and retained only those loci supported by a read depth between 5 × and 100 × (8895 loci removed). Three individuals were removed from the dataset owing to an exceptionally poor call rate of less than 0.5 (MS9, MS12, MS8) and resultant monomorphic loci removed from the dataset. Loci with a repeatability less than 0.99 were removed (2380 loci) and finally, loci with a call rate of less than 0.95 were removed (9870 loci). We regard the data remaining after this additional filtering (11,317 SNP markers for 75 individuals) as highly reliable.

Visualization

Genetic similarity among individuals, populations and colour morphs was visualized using ordination (principal coordinates analysis or PCoA⁴⁷) as implemented in the gl.pcoa and gl.pcoa.plot functions of dartR. A scree plot of eigenvalues guided the number of informative axes to examine⁴⁸, 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.²⁶, R scripts are provided in the online data repository), and pairs of populations for which the count of fixed differences did not exceed the estimated false positive rate were also amalgamated. Counts of private alleles were obtained from gl.report.pa in dartR.

Genetic diversity

Expected heterozygosity was used as a measure of relative genetic diversity. Heterozygosity was obtained for each population from allele frequencies using the gl. report.heterozygosity function of dartR and pairwise comparisons of heterozygosity between populations were tested for significance using the gl.test.heterozygosity function in dartR (significance evaluated by re-randomizaton with 10,000 replicates).

Hybridisation

The genotypes of suspected hybrids/introgressed individuals (T1 and T5) were examined using New Hybrids²⁸ 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).

Source: Ecology - nature.com