Discovered and disappearing? Conservation genetics of a recently named Australian carnivorous marsupial

Abstract Five new species within the Australian carnivorous marsupial genus Antechinus have recently been named, at least two of which are threatened. Important facets of the habitat use and extinction risk of one of these new species, the buff‐footed antechinus, A. mysticus, are not well understood. Previous research has suggested that the species utilizes a broad range of inter‐connected forest habitats in southeast Queensland (Qld), Australia. Based on this potentially connected habitat, we predicted that A. mysticus should have low population genetic structure, particularly in relation to its congener, the spatially restricted, high altitude, closed‐forest A. subtropicus. We genotyped nine microsatellite loci for six populations of A. mysticus, sampled throughout their known range in eastern Australia, and compared them with four proximate populations of A. subtropicus. Surprisingly, genetic structuring among southeast Qld populations of A. mysticus was moderate to high and similar to that between A. subtropicus populations. We postulate that all A. mysticus populations have declined recently (<100 generations), particularly the northernmost southeast Qld population, which may be at risk of extinction. Our results suggest that A. mysticus is limited to a more scattered and fragmented distribution than previously thought and may be in decline. The identification of population decline in this study and recently in other Antechinus suggests the extinction risk of many Australian mammals should be reassessed.

have been widely studied as a model in breeding biology (Naylor, Richardson, & Mcallan, 2008). Such interest has also encouraged a number of genetic studies on Antechinus, with a principal focus on deeper level systematics (Mitchell et al., 2014;Westerman et al., 2015), taxonomy (Baker, Mutton, Mason, & Gray, 2015;Van Dyck, 2002) and a range of single population studies mainly focusing on habitat use (Banks et al., 2005;Kraaijeveld-Smit, Lindenmayer, Taylor, Macgregor, & Wertheim, 2007;Lada, Thomson, Mac Nally, & Taylor, 2008). However, genetic studies of the genus have generally not been undertaken across a species' total geographic range. Such a comprehensive evaluation may provide invaluable information on the life history, habitat use, and extinction risk of a species. Recently, five new Antechinus species have been named based on combined phylogenetic, morphology and breeding biology data, and at least two of these are at risk of extinction (Baker, Mutton, & Hines, 2013;Baker, Mutton, Hines, & Van Dyck, 2014;Baker, Mutton, & Van Dyck, 2012;Baker et al., 2015). Currently, little is known about the ecology or population genetics of these five species. Consequently, there is a vital need for fine-scale, yet geographically detailed molecular studies on the suite of newly discovered species. Understanding species' ranges, population structure, demography, and habitat preference is crucial for effective conservation and is particularly important in Australia where mammals have suffered an extremely high rate of extinction in the last 200 years (Woinarski, Burbidge, & Harrison, 2015).
Based on the fragmented nature of A. subtropicus' habitat, and the potential connectivity of A. mysticus' habitats, we expected that the relative connectivity of the two species to be reflected in differing population genetic signatures. The lower altitude, drier forest habitats of A. mysticus would presumably have been more connected in recent geological time than the comparative "islands" of high altitude vine forest and rainforest which A. subtropicus favor (Byrne et al., 2011). Over time, such geographic partitioning could lead to a relatively higher degree of genetic structuring among A. subtropicus populations in southeast Qld. Mitochondrial cytochrome b (Cytb) results support this contention, with A. mysticus showing very little diversity in southeast Qld, comprising just two haplotypes, only one base pair different (0.2%), across its range (Baker et al., 2012;Mutton, 2017). In contrast, A. subtropicus showed four times as much genetic divergence (0.8%) over only half the geographic distance (75 km; Mutton, 2017) and greater divergence in Cytb has been recorded within antechinus populations and across similar geographic distances in a number of species (Beckman, Banks, Sunnucks, Lill, & Taylor, 2007;Mutton, 2017). The two species, which are not closely related within the genus, show a deep interspecies divergence (13.3%-14.3% at Cytb -Baker et al. (2012)). The differing habitat preference of these species also appears to hold in known areas of sympatry. The two species can co-occur at mid-altitude (~250 m) and transitional environments/ecotones where open forest and closed, vine forest, and rainforest communities merge, but each occur in the absence of congeners in higher altitude rainforest (A. subtropicus) and lower altitude open forest (A. mysticus) of southeast Qld (Baker et al., 2012;Mutton et al., 2017).
To test the hypothesis of greater population connectivity in A. mysticus than A. subtropicus, we sampled the majority of known populations of A. mysticus and representative populations throughout the range of A. subtropicus. We compared the species' population structures using nine microsatellite loci and examined levels of genetic diversity within and between populations, which is essential for effective conservation management. Microsatellites evolve more rapidly than mitochondrial markers, allowing more recent divergence and finer scale genetic structuring to be revealed. Mitochondrial genes can also be biased by sex-specific dispersal patterns, as they are maternally inherited. In contrast, microsatellites are inherited from both parents and population size and stochastic factors can have a significant influence on genetic estimates (Selkoe & Toonen, 2006). These attributes make microsatellites an ideal marker to test the ideas posed in this study.

| A. mysticus
Samples were collected from across the known geographic range of A. mysticus, including five localities in southeast Qld and one from Eungella, mid-east Qld (Figure 1) (Baker et al., 2012;Mutton et al., 2017). Despite extensive trapping effort (cumulatively more than 10,000 trap nights across eight locations), A. mysticus has not been found in the approximately 700 km straight line distance (SLD) separating southeast Qld from Eungella or from other localities in mideast Qld (Baker et al., 2012Mutton et al., 2017). In southeast Qld, our sampled A. mysticus sites ranged from 15 to 155 km SLD F I G U R E 1 Map of southeast Queensland showing remnant vegetation and the Antechinus sites sampled in the present study. Circles represent A. mysticus populations, triangles represent A. subtropicus populations. The square represents a site at which both species were caught apart and 130 to 230 m above sea level (ASL) (Appendix S1 and S2).
Eungella is located 709-810 km SLD from the other sites and at a much higher altitude (750 m ASL) than other A. mysticus populations (Appendix S1 and S2). Study sites were situated in areas of protected remnant vegetation (Figure 1). The southeast Qld A. mysticus sites were predominantly riparian open woodland, while the mid-east Qld site was dominated by rainforest (Pearce, 2016) (See Appendix S3 for vegetation community information). D'Aguilar was the only site surveyed (and known) at which both species occur (Figure 1). This site is typified by an abrupt merging of riparian open woodland and rainforest communities .

| A. subtropicus
Antechinus subtropicus samples were collected from four sites covering the entire known range of the species, except the southern latitudinal maximum at Border Ranges ( Figure 1). Despite numerous surveys over a 3-year period, including more than 5,000 cumulative trap nights at several locations, we were unable to trap A. subtropicus at this location (Mutton, 2017). Three of our A. subtropicus sites were approximately equidistant (~50 km SLD) (Wrattens, Woondum and Conondale), while the fourth site, D'Aguilar, was located 75-124 km SLD south of these sites ( Figure 1). The sites were located within patches of remnant vegetation ( Figure 1). Excluding D'Aguilar, where the two species co-occur, all A. subtropicus sites were high altitude, ranging from 460-to 708 m ASL (Appendix S1 and S2). The sites were dominated by Eucalyptus species and notophyll vine forest (Appendix S3).

| Mapping
Southeast Qld study site locations were mapped in relation to remnant vegetation cover (DEHP, Queensland Government 2009) and elevation (Geoscience Australia 2011) using QGIS software V 2.12.0 (QGIS Development Team 2017). Remnant vegetation cover described native vegetation that had not been cleared or had been cleared but retained a dominant canopy >70% of the height and >50% of the cover relative to the undisturbed height and cover of that layer and was dominated by species characteristic of the vegetation's undisturbed canopy (DEHP, Queensland Government 2017).
Elevation data were generated using a 1-s SRTM Derived Digital Elevation Model. Elevation was categorized and color-coded into 175 m increments between 0 m and ≥700 m.

| Taxon sampling
Species were identified by eye based on pelage color as described in Baker et al. (2012), and their identity confirmed using mtDNA (Cytb) sequencing (Mutton, 2017), prior to microsatellite screening.
Genetic samples were obtained through field collection with type A Elliott traps (Elliott Scientific, Vic, Australia). Ear tissue samples were collected from individuals captured in the field and stored in 80% ethanol. Samples from Eungella and D'Aguilar were collected during parallel capture-mark-recapture studies described in Pearce (2016) and Mutton et al. (2017), respectively. Cooloola samples were collected over a number of months in approximately 9,000 trap nights.
Samples from the other sites were collected by trapping for 3-5 days and between 500 and 1,000 trap nights at each site during the winter prebreeding period.
Nine microsatellite loci (see Table 1) were amplified for both species using a Qiagen multiplex kit (Qiagen, Dusseldorf, Germany).
All primers were originally designed for A. agilis, except one (Aa7Q), which was designed for A. flavipes (see Table 1). As allele sizes of some loci overlapped, two multiplexes were amplified separately (Table 1)

| Genetic variation
Allele size was scored in GeneMapper 4.1 (Applied Biosystems). The total number of alleles (A), private (unique) alleles (uA), rare alleles (rA; frequency ≤5%), observed (H o ), and expected (H e ) heterozygosity per population were estimated using GenAlEx 6.502 (Peakall & Smouse, 2006). Allelic richness was standardized for sample size (10 individuals) (AR) and estimated in Fstat 2.9.3.2 (Goudet, 1995). As one locus (Aa70) had low sample sizes at two A. subtropicus sites, it was excluded from AR analysis for this species. Wilcoxon sign-rank tests were implemented in R to compare differences in AR, rA, and H e between populations of each species.
Tests for linkage disequilibrium and departures from Hardy-Weinberg Equilibrium (HWE) for each locus-population combination were carried out in Genepop 4.2 (Rousset, 2008) using 1,000 permutations to test for statistical significance. To assess whether sample effort was adequate to capture the allelic diversity of each population, we undertook rarefaction analysis using the "jackmsatpop" function of the R package PopGenKit 1.0 (Paquette, 2012).
Rarefaction curves were generated by calculating the sample size adjusted allelic diversity at stepwise increases of one until the sample size of each population was reached. Two hundred jackknife replicates were performed at each interval.

| Population structure and genetic bottlenecks
Genetic structure between populations was analyzed using pairwise F ST (Weir & Cockerham, 1984). As all measures of population differentiation have limitations, it has been suggested that multiple population differentiation statistics be examined (Meirmans & Hedrick, 2011) . Therefore, we also calculated D EST in GenAlEx 6.502 (Jost, 2008). Hierarchical structuring of populations was assessed by an Analysis of Molecular Variance (AMOVA; Excoffier, Smouse, & Quattro, 1992). The relationship between geographic distance and genetic distance (pairwise F ST ) was analyzed using Mantel tests (Mantel, 1967). Pairwise F ST , AMOVA, and Mantel tests were undertaken in Arlequin 3.5.2.2 (Excoffier & Lischer, 2010). The pairwise F ST critical value (α) was corrected for multiple tests using the BY False Discovery Rate method, which accounts for Type I error with less loss of power than the Bonferroni adjustment (Narum, 2006).
We also conducted redundancy analysis (RDA), which allowed us to compare with Mantel tests, and more accurately calculates isolation by distance and the influence of geographic variables on genetic structure (Legendre & Fortin, 2010;Meirmans, 2015). RDA was undertaken using the R package Vegan 2.4-5 (Oksanen et al., 2017;R Core Team, 2015). We used allele frequencies as the dependent variable and four independent variables: latitude, longitude, degree of isolation (IS), and minimum distance to neighboring site (DN). IS was measured as the mean distance to the closest three sites. Stepwise AIC comparisons (using the function "ordistep") of the full and nested RDA models were undertaken to determine the optimal model. Genetic evidence for a recent (2-4 Ne generation) reduction in population size was tested using BOTTLENECK 1.2.02 (Piry, Luikart, & Cornuet, 1999). A two-phased mutation model (TPM) incorporating 95% single-step changes, 5% multi-step changes with a variance among multiple-step changes set to 12% was used (Piry et al., 1999). Following the recommendations of Piry et al. (1999), statistical significance was evaluated from 5,000 simulations of the one-tailed Wilcoxon sign-rank test. For each population, Garza and Williamson's (2001) M-ratio was calculated using Arlequin 3.5.2.2.
M-ratio is a measure of the proportion of unoccupied allelic states given the allele size range. It is sensitive to population bottlenecks (often for over 100 generations) as the ratio reduces when alleles are lost due to random drift after a bottleneck occurs (Garza & Williamson, 2001). The inbreeding coefficient (F IS ) was estimated for all loci in each population using Arlequin 3.5.2.2 (Excoffier & Lischer, 2010), with 10,000 permutations of alleles among individuals within a population to test for significance.
To define the number of distinct population groups (K) within each species, Bayesian clustering of individuals without prior assignment to population was performed using the software package STRUCTURE 2.3.4 (Pritchard, Stephens, & Donnelly, 2000). The program was run for 10 iterations of 100,000 generations, with an initial burn-in of 20,000 generations. The value of K was set from 1 to 10, with ten replicates of each K to verify the convergence of the Markov chain. No assumptions were made about the shared descent of populations, allele frequencies were set to uncorrelated, and separate alpha values were used for each population (Pritchard et al., 2000). The program assigns individuals to K clusters, with the user nominating which value of K is most appropriate for their data. Evanno, Regnaut, and Goudet (2005) recommend that the highest value of ∆K be taken as the "true" K value. However, it has also been suggested that the value of K which captures the majority of structure in the dataset be used (Pritchard et al., 2000), a technique which has been implemented in a number of ecological studies

Locus
Primer sequence Reference Multiplex TA B L E 1 Primer sequence and references of the nine microsatellites loci genotyped in this study. The primers were divided into two multiplex groupings (Bryant & Fuller, 2014;Cardoso, Mooney, Eldridge, Firestone, & Sherwin, 2014;Krosch et al., 2013). Thus, ∆K values were plotted and STRUCTURE results were summarized using Structure Harvester (Earl & Vonholdt, 2012). Cluster membership coefficient matrices for each K value were summarized using CLUMPP 1.1.2 (Jakobsson & Rosenberg, 2007) and visualized as admixture graphs with Distruct 1.1 (Rosenberg, 2004).

| Descriptive statistics
In total, 101 A. mysticus (39% female) and 75 A. subtropicus (31% female) individuals were screened for nine microsatellite loci. One locus, Aa7Q, failed to amplify for all A. mysticus, presumably due to a mutation in the primer binding region (Selkoe & Toonen, 2006 Pairwise comparison between all individual loci revealed low rates of significant linkage disequilibrium for A. mysticus (4.8%) and A. subtropicus (4.0%). As there was no consistent pattern of linkage across sites or loci (data not shown), it was deemed unlikely that these low rates of linkage disequilibrium were due to physical linkage (Selkoe & Toonen, 2006 F IS analysis showed evidence of significant inbreeding ( Table 2).
Evidence of a significant genetic bottleneck was also found in one A. subtropicus site (Wrattens), which was also the only A. subtropicus site with an M-ratio substantially below the critical value Garza and Williamson (2001)

| A. mysticus
Pairwise F ST estimates revealed significant structure among all A. mysticus sites ( Stepwise selection revealed that no geographic variables significantly explained genetic variation when RDA included all populations, or when the divergent Eungella and Cooloola populations were both excluded. However, when only the Eungella population was excluded from the analysis, the optimal model included degree of isolation (IS) as an explanatory variable (p = 0.025) (Figure 2a).
This model explained 38.8% of variance. When the Cooloola population was excluded, distance to nearest neighbor (DN) was revealed as a significant explanatory factor (p = 0.016), and the model containing this factor explained 47.1% of observed variation (Figure 2b).
Bayesian clustering analysis revealed a pattern of genetic structure concordant with the F ST analysis. The Evanno et al. (2005) method identified two groups (K = 2) (Appendix S5a). Graphical representation showed one group to largely correspond to Cooloola and the other group to contain all other sites (Figure 3a). However, exploration of additional groupings revealed further population structuring. The geographically isolated Eungella site was revealed at K = 3 as a separate cluster (Figure 3b). Cooloola represented a second cluster which showed a small amount of admixture with the third grouping, which encompassed the other southeast Qld sites (Figure 3b).
Increasing K to higher values increased admixture and did not reveal further biologically informative groupings (data not shown).

| A. subtropicus
Of the three equidistant (~50 km) sites, Conondale was more divergent from the other two sites (F ST = 0.062-0.096) than Woondum and Wrattens were from each other (F ST = 0.044) ( Bayesian clustering analysis also revealed D'Aguilar as the most Graphical representation of K = 3 also revealed D'Aguilar as a largely discrete grouping, but additionally showed Wrattens to be separate from the Woondum/Conondale sites, which formed a third grouping ( Figure 4b).

| D ISCUSS I ON
Microsatellite genotyping indicated limited gene flow among populations of A. mysticus across its range. This is contrary to our expectations based on ecology, distribution, habitat preference, and mitochondrial gene sequencing of this species (Baker et al., 2012;Mutton, 2017). Indeed, the degree of population genetic structuring appears similar to that observed among populations of scattered high altitude, closed-forest A. subtropicus.

| Genetic differentiation and variation between A. mysticus, A. subtropicus, and other carnivorous marsupials
Two highly differentiated A. mysticus populations (Eungella and Cooloola) were revealed in this study. As heterozygosity and allele frequency distribution strongly influence F ST results, comparisons of F ST values between taxa need to be made with caution (Jakobsson, Edge, & Rosenberg, 2013). However, a general comparison of trends can be revealing. It is notable that the differentiation between each of these and all other A. mysticus populations was high, being comparatively greater than that recorded between some mainland F I G U R E 2 RDAs showing the contribution of spatial factors to genetic structure in A. mysticus, optimal models are shown for when (a) Eungella and (b) Cooloola sites were respectively excluded from the analysis. Each site is shown in RDA space and circles show how explanatory variables fall in the RDA space. IS: degree of isolation and DN: distance to nearest neighbour F I G U R E 3 Graphical representation of membership coefficients of the Bayesian STRUCTURE analysis of 7 microsatellite loci for A. mysticus obtained from 6 sites across the known range of the species. Each plot represents different population assignments for K: (a) K = 2; and (b) K = 3. Solid black lines delineate the 6 different sites; each vertical line represents a single individual. Colors represent cluster assignments Australian and Vulnerable Tasmanian populations of the tiger quoll, Dasyurus maculatus (Firestone, Houlden, Sherwin, & Geffen, 2000), and at a similar level to that found between populations of A. agilis with different teat numbers, which may represent incipient species (Appendix S6) (Beckman et al., 2007;Draper, 2012).
Even if the two highly differentiated A. mysticus populations (Eungella and Cooloola) were excluded, genetic variation was significant and comparable between populations of both Antechinus species in southeast Qld. The microsatellite differentiation among these populations in both species was larger than that recorded between populations of the broad-ranging A. flavipes flavipes in southeast Australia (Lada et al., 2008) and comparable to same teat number populations of fragmented A. agilis (Beckman et al., 2007) (Mills, Moro, & Spencer, 2004).
In the present study, levels of genetic variation (AR, H e ) were low to moderate for all populations of both species except A. mysticus at Cooloola, which had much lower variation. These values were comparatively lower than previously reported for other Antechinus (Appendix S6). While variation in sample size can influence these values, it is notable that genetic variation in both species was similar to averages reported in the last review of marsupial population genetic diversity (Eldridge, 2010), in which the majority of studied taxa were listed as threatened species. For instance, levels of genetic variation found in the current study were similar to that found for well-connected populations of two threatened marsupials, the Vulnerable brush-tailed rock-wallaby, Petrogale penicillata (Hazlitt, Goldizen, & Eldridge, 2006) and mainland populations of the Endangered northern quoll, Dasyurus hallucatus (Cardoso et al., 2009). However, genetic variation was higher here than that reported for two other threatened marsupials which occur in wellconnected habitats, the eastern quoll, Dasyurus viverrinus (Cardoso et al., 2014) and Tasmanian devil, Sarcophilus harrisii (Jones, Paetkau, Geffen, & Moritz, 2004;Appendix S6). In summary, the genetic variation of A. mysticus and A. subtropicus appears to be within the reported range of some threatened carnivorous marsupials. Our genetic structure analysis also suggests that A. mysticus may have experienced a decline and be at high risk of localized extinction at Cooloola (see below).

| Genetic differentiation and variation among A. mysticus populations
Given its geographic (>700 km) and genetic (>2% at Cytb) isolation from known conspecifics (Baker et al., 2012;Mutton, 2017) (Weir & Cockerham, 1984). However, few unique and rare alleles were found in the Cooloola population, and the opposite would be expected if the population had long been isolated (Frankham, Briscoe, & Ballou, 2002). In total, these results suggest that the strong differentiation of the Cooloola A. mysticus from other populations is most likely driven by low abundance and a genetic bottleneck rather than long-term divergence from conspecifics or incipient speciation. Consequently, the species may be at risk of localized extinction at Cooloola.
However, incipient speciation cannot be completely ruled out, as it has been suggested for A. agilis populations which show similarly deep microsatellite, but low Cytb divergence (Beckman et al., 2007;Draper, 2012 (Mutton, 2017) and morphological variation within and between Antechinus species (Baker et al., 2015). Furthermore, strong selection pressure on Cytb is often associated with very high altitude environments (Zhang, Lin, Nevo, Yang, & Su, 2013), which was not a factor in the present study. Most likely, the low Cytb and high microsatellite divergence observed within southeast Qld A. mysticus is primarily driven by the different evolutionary speeds of the two marker types (Sunnucks, 2000). Such an explanation has also been invoked to explain a similar discrepancy between Cytb and microsatellite results in A. agilis (Beckman et al., 2007).

| Comparative population structure between A. mysticus and A. subtropicus
It was hypothesized that A. mysticus utilized a broader and more con- this is yet to be tested in Queensland, A. flavipes in southeastern Australia is less structured than either of the species in the present study (Lada et al., 2008).

| Conservation
M-ratios were generally much lower in A. mysticus than A. subtropicus (Table 2). This suggests A. mysticus was previously more connected but has undergone population contractions throughout its range (Garza & Williamson, 2001;Nei, Maruyama, & Chakraborty, 1975). Southeast Qld has experienced a remarkably high rate of land clearing in the 230 years since European colonization (Figure 1).
The lower altitude habitats that A. mysticus favors are also easier to clear and utilize for humans and thus have experienced higher clearing rates (Bradshaw, 2012). This is apparent in two national parks  Table 2).
As well as land clearing, threats such as introduced predators (feral cats, Felis catus and foxes, Vulpes vulpes) and climate change may have caused population reductions in A. mysticus and will likely continue to negatively affect both species and many other mammals (Woinarski et al., 2015). Indeed, a recent microsatellite study identified that land clearing and related anthropogenic changes are driving a substantial decline in the population size of southwestern Australian A. flavipes (Mijangos, Pacioni, Spencer, Hillyer, & Craig, 2017). However, the evidence of decline throughout the range of A. mysticus, which is possibly severe at Cooloola, suggests that A. mysticus is at a more immediate extinction risk than A. subtropicus. Therefore, we believe A. mysticus could warrant threatened species listing. We advocate regular monitoring of A. mysticus and prioritization of habitat protection for this species.
Australia has experienced an extraordinarily high rate of mammalian extinction in the last ~200 years (Woinarski et al., 2015).
Antechinus have previously been considered less at risk of extinction than many Australian mammals in a similar size range.
However, the present study and recent listing of four species of Antechinus as at threat of extinction suggest this is no longer the case. There is little information on the population structure and abundance of many of Australia's largely endemic marsupial mammals (Woinarski, Burbidge, & Harrison, 2014). In light of the recent declines in Antechinus and the ongoing decline of Australian mammals more broadly, it is apparent that the extinction risk of many of these taxa not currently listed as at risk of extinction should now be reassessed.

ACK N OWLED G M ENTS
We express our deep gratitude to Coral Pearce for providing the Eungella samples and for her assistance on multiple field trips.
Thanks also go to Rebecca Wheatley and her collaborators for providing the Cooloola samples. We thank all volunteers for their high-quality assistance and company in the field: David Benfer, Rob

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N
TYM, SJF, and AMB designed the study. TYM collected and analyzed the data. DT constructed the maps and provided the vegetation data. TYM wrote the manuscript with assistance from all authors.