Disentangling the Genetic Relationships of Three Closely Related Bandicoot Species across Southern and Western Australia

: The taxonomy of Australian Isoodon bandicoots has changed continuously over the last 20 years, with recent genetic studies indicating discordance of phylogeographic units with current taxonomic boundaries. Uncertainty over species relationships within southern and western Isoodon , encompassing I. obesulus , I. auratus , and I. fusciventer , has been ongoing and hampered by limited sampling in studies to date. Identiﬁcation of taxonomic units remains a high priority, as all are threatened to varying extents by ongoing habitat loss and feral predation. To aid diagnosis of conservation units, we increased representative sampling of I. auratus and I. fusciventer from Western Australia (WA) and investigated genetic relationships of these with I. obesulus from South Australia (SA) and Victoria (Vic) using microsatellite markers and mitochondrial DNA. mtDNA analysis identiﬁed three major clades concordant with I. obesulus (Vic), I. auratus , and I. fusciventer ; however, I. obesulus from SA was polyphyletic to WA taxa, complicating taxonomic inference. Microsatellite data aided identiﬁcation of evolutionarily signiﬁcant units consistent with existing taxonomy, with the exception of SA I. obesulus . Further, analyses indicated SA and Vic I. obesulus have low diversity, and these populations may require more conservation efforts than others to reduce further loss of genetic diversity.


Introduction
Threatened species often require active management to ensure their persistence. An important first step in threatened species conservation is identifying and defining appropriate units of management to direct conservation effort to where it is needed the most [1,2]. Confusion about taxonomic boundaries of subspecies and other conservation units has the potential to lead to inappropriate management decisions that may have detrimental consequences. For example, mixing populations that are genetically and evolutionarily distinct could result in the loss of local adaptation or outbreeding depression [3,4]. Conversely, managing populations as separate units, particularly when population sizes are small, may lead to increased levels of inbreeding and loss of adaptive potential due to random genetic drift [5]. Developing an appropriate taxonomic classification that reflects the actual nature of species or subspecies is therefore essential to guide management decisions to improve conservation outcomes.
The concept of "evolutionarily significant units" (ESU) has been used to aid management of populations in recent decades. The concept was introduced by Ryder [6] to resolve "subspecies" into manageable units where an ESU is defined as a group of organisms with high genetic and ecological distinctiveness. With increasing use of molecular tools, Moritz [7] proposed a genetic criterion for defining an ESU as a group of organisms that are reciprocally monophyletic for mtDNA haplotypes and have diverged allele frequencies at nuclear loci, a definition that places emphasis on historical isolation of populations. A sub-category under an ESU is "management unit" (MU), which recognises populations that do not show reciprocal monophyly for mtDNA haplotypes but have diverged haplotype and allele frequencies at mitochondrial or nuclear loci to accommodate more contemporary population structuring [7]. The use of reciprocal monophyly as a criterion for defining an ESU, whilst popular, has been criticised as being too stringent [8,9]. Fraser and Bernatchez's [8] concept of "adaptive evolutionary conservation" offers a more flexible approach, which defines an ESU as "a lineage demonstrating highly restricted gene flow from other such lineages within the higher organization level (or lineage) of the species". These lineages are sufficiently isolated that each lineage has limited or no impact on the evolution, genetic variance, and demography of other lineages [8]. This approach integrates multiple ESU concepts, including Moritz's, to fulfil conservation goals, which under differing circumstances will warrant the use of some concepts over others.
The short-nosed bandicoots, genus Isoodon (family Peramelidae), is a small-to-mediumsized Australian marsupial that can be found across Australia. Bandicoots are terrestrial, omnivorous marsupials with a bounding gait and forelimbs adapted for digging for food [10]. As digging mammals, they play an important role as ecosystem engineers, contributing to altered soil nutrient and moisture dynamics that enhance plant germination and growth [11]. Their "rat-like" appearance and cryptic nature make them less attractive for conservation efforts. However, they are known for their complex biology of adaptive traits and evolutionary diversity [11,12]. These traits enable some Isoodon species to adapt to urban environments [13,14]. However, facing pressure from introduced predators and human disturbances, many Isoodon species have dramatically reduced in distribution since European settlement and persist in remnants of their formal ranges.
The taxonomy of the genus Isoodon has changed multiple times over the last 70 years and is currently still in flux. Up to eleven different species of Isoodon have been recognised at various times [15]; however, only three are formally recognised currently: I. macrourus, I. obesulus, and I. auratus [16,17], with additional taxa I. peninsulae [18] and I. fusciventer [19] recently proposed, both raised from subspecies (of I. obesulus) to full species rank. Most recently, comprehensive geographic and taxonomic sampling in a phylogeographic study of I. obesulus and I. auratus has indicated further complexity in the relationships of species and subspecies within these taxa [12]. One major finding was that the geographic distribution of I. o. obesulus is significantly more restricted than previously recognised, with specimens of I. o. obesulus from South Australia (SA; encompassing Mount Lofty Ranges, Fleurieu Peninsula, and Kangaroo Island) genetically differentiated from I. o. obesulus from south-eastern Australia and more closely allied with I. fusciventer from Western Australia (WA; [12,20]). Further, consistent with early molecular studies [15,21], Cooper et al. [12] failed to resolve a clear pattern of reciprocal monophyly between I. fusciventer (southern WA) and I. auratus (northern WA) in mitochondrial DNA, evidence that has previously been used to suggest the species be synonymised [15,21]. This suggestion based on genetic data contrasts with morphological data indicating that I. obesulus and I. auratus can be readily distinguished by a range of characters, including size, skull and teeth characters, and fur colour [15,19,22,23]. While monophyly could not be resolved, Cooper et al. [12] identified a lack of mitochondrial DNA haplotype sharing between the taxa as well as fixed allozyme differences and private nuclear gene haplotypes, suggesting a putative pattern of paraphyly of these taxa.
Here, we focus on identifying conservation units within the "southern and western" group of Isoodon bandicoots, WA taxa Isoodon auratus, I. a. barrowensis and I. fusciventer, and southern Australian taxa I. obesulus (Vic and SA), all of which are of conservation concern. Isoodon species have dramatically declined in their abundance and distribution in the last 220 years since European settlement of Australia [24,25], primarily as a result of introduced predators such as feral cats, foxes and black rats, altered fire regimes, disease, and habitat loss from clearing of native vegetation and habitat modification [26][27][28]. Isoodon obesulus, I. fusciventer, and I. auratus are listed as Near Threatened, Least Concern, and Near Threatened, respectively (International Union for Conservation of Nature red list/Environment Protection and Biodiversity Conservation (EPBC) Act, Australia). The subspecies I. o. obesulus in eastern and south-eastern Australia is listed as Endangered, although re-circumscription as suggested by analyses in Cooper et al. [12] means this subspecies is now more geographically restricted requiring re-evaluation of its conservation status; whilst other subspecies, I. a. auratus and I. a. barrowensis, are listed as Vulnerable under the EPBC act. These species/subspecies are currently subjected to different levels of conservation action to mitigate threats [26], with I. a. auratus and I. o. obesulus particularly targeted for translocations and reintroductions to island and mainland feral predator-free areas to increase their population size [29]. Clarification of taxa and conservation units is required to support conservation assessment and listing across the group and to ensure future conservation management activities are appropriately applied.
In this study, we expand sampling of representative I. auratus from northern WA and I. fusciventer from southern WA, as well as introduce additional molecular data from hyper-variable nuclear loci (microsatellites), to provide further resolution to the relationships amongst southern and western Isoodon bandicoots as raised in Cooper et al. [12]. Incorporating existing molecular data from Cooper et al. [12] and Li et al. [20] with our expanded sampling, we use a combination of phylogenetic and hierarchical population genetic approaches to identify putative ESUs across the southern and western Isoodon bandicoots, I. auratus, I. fusciventer, I. obesulus SA, and I. o. obesulus Vic. Within each regional population, we investigate geographical sub-structuring and assess patterns of mitochondrial and microsatellite genetic diversity to inform priorities for conservation genetic management.

Sample Collection, Microsatellite Genotyping, and Mitochondrial DNA Sequencing
We collated microsatellite and mitochondrial DNA (mtDNA) profiles from previously published sources with new sampling for a total of 731 and 218 additional samples, respectively. Isoodon obesulus, I. fusciventer, and I. auratus samples were distributed across multiple locations in Western Australia, South Australia, and Victoria and across multiple years from 2002 to 2018 ( Figure 1, Table S1). Microsatellite profiles for WA specimens of I. auratus and I. fusciventer (n = 172) were predominantly obtained from previously published datasets in Ottewell et al. [30] and Ottewell et al. [13], respectively, with some additional new sampling from southwest WA (I. fusciventer) and the Kimberley region (I. a. auratus) (Table S1). Ear biopsy samples were opportunistically obtained by environmental consultants during urban fauna relocations and during population monitoring by the Australian Wildlife Conservancy (Kimberley samples). A small number of samples were obtained from roadkill animals by citizen scientistsunder the Department of Biodiversity, Conservation and Attractions license to use animals for scientific purposes (permit no. U10/2018). DNA extraction and microsatellite genotyping were conducted as described in Ottewell et al. [30] using 10 microsatellite primers (B3-2, B20-5, B34-2, Ioo2, Ioo4, Ioo6, Ioo7, Ioo8, Ioo10, and Ioo16). Briefly, amplification reactions were carried out using the Qiagen Multiplex Kit, following the manufacturer's instruction (Qiagen, Hilden, Germany). PCR products were separated on an Applied Biosystems 3130xl capillary sequencer (Foster City, CA, USA). Fragment sizes were determined with GeneScan-500 LIZ size standard (Applied Biosystems) and scored in Genemapper v5.0 software (Applied Biosystems). Microsatellite profiles from published datasets were re-scored with new samples to ensure consistency in allele scoring. Microsatellite genotypes from SA and Vic were obtained from Li et al.'s [20] study (n = 559). Due to the addition of M13 tags to primers used in this study, which added 30 bp to tandem repeat sizes, we re-amplified a subset of 15 samples (n = 5 from each of Mount Burr, Mount Lofty Ranges, Grampians populations) in our WA laboratory to calibrate genotype scores between SA/Vic and WA datasets.
primers used in this study, which added 30 bp to tandem repeat sizes, we re-amplified a subset of 15 samples (n = 5 from each of Mount Burr, Mount Lofty Ranges, Grampians populations) in our WA laboratory to calibrate genotype scores between SA/Vic and WA datasets. Sequences from the highly variable non-coding mitochondrial control region (Dloop) were obtained from Li et al. [20] (n = 49) and Cooper et al. [12] (n = 23, Table S1) with additional sequences obtained for a subset of the WA samples used in microsatellite genotyping (n = 146). New specimens were sequenced using primers L15999M and H16498M [31] with amplifications carried out in 25 μL reactions, which included ~10-20 ng of template DNA, 0.5 μM of each primer, 1.25 mM MgCl2, 0.2 μM of each dNTPs, 0.05 μL of 10% BSA, 1× Reaction Buffer, and one unit of Invitrogen Taq polymerase. Reactions were run under the following conditions: 94 °C for 3 min followed by 35 cycles of denaturation 94 °C for 30 s, annealing at 55 °C for 45 s, and extension at 72 °C for 60 s, with a final extension of 72 °C for 7 min. Products were visualised on 1% agarose gels and Sanger-sequenced at a commercial service (Australian Genome Research Facility). Sequences were edited and aligned using SEQUENCHER v5.2 (Gene Codes Corporation, Ann Arbor, MI, USA) with consensus sequences from each of the datasets then aligned using CLUSTAL W with default parameters [32] and sequences trimmed to the same length in MEGA v10.0.5 [33].

Phylogenetic Analysis
Phylogenetic analyses of mtDNA sequences were carried out using maximum likelihood in RAXML v7.0.3 [34], neighbour-joining in MEGA v10.0.5 [33] and Bayesian inference in BEAST v1.10.4 [35]. In RAXML, we applied a single model of evolution, Genera Time Reversible (GTR) model [36] with unequal variation at sites modelled using a Gamma (G) distribution [37] to the sequence data. Rapid bootstrap analysis was applied  Sequences from the highly variable non-coding mitochondrial control region (Dloop) were obtained from Li et al. [20] (n = 49) and Cooper et al. [12] (n = 23, Table S1) with additional sequences obtained for a subset of the WA samples used in microsatellite genotyping (n = 146). New specimens were sequenced using primers L15999M and H16498M [31] with amplifications carried out in 25 μL reactions, which included ~10-20 ng of template DNA, 0.5 μM of each primer, 1.25 mM MgCl2, 0.2 μM of each dNTPs, 0.05 μL of 10% BSA, 1× Reaction Buffer, and one unit of Invitrogen Taq polymerase. Reactions were run under the following conditions: 94 °C for 3 min followed by 35 cycles of denaturation 94 °C for 30 s, annealing at 55 °C for 45 s, and extension at 72 °C for 60 s, with a final extension of 72 °C for 7 min. Products were visualised on 1% agarose gels and Sanger-sequenced at a commercial service (Australian Genome Research Facility). Sequences were edited and aligned using SEQUENCHER v5.2 (Gene Codes Corporation, Ann Arbor, MI, USA) with consensus sequences from each of the datasets then aligned using CLUSTAL W with default parameters [32] and sequences trimmed to the same length in MEGA v10.0.5 [33].

Phylogenetic Analysis
Phylogenetic analyses of mtDNA sequences were carried out using maximum likelihood in RAXML v7.0.3 [34], neighbour-joining in MEGA v10.0.5 [33] and Bayesian inference in BEAST v1.10.4 [35]. In RAXML, we applied a single model of evolution, General Time Reversible (GTR) model [36] with unequal variation at sites modelled using a Gamma (G) distribution [37] to the sequence data. Rapid bootstrap analysis was applied indicates sampling locations. Each pie chart represents haplotype frequencies within each sampling site, and size of pie chart represents numbers of haplotypes. Colour codes indicate the sites where haplotypes were mostly found. Pattern codes represent individual haplotypes found within each site.
Sequences from the highly variable non-coding mitochondrial control region (D-loop) were obtained from Li et al. [20] (n = 49) and Cooper et al. [12] (n = 23, Table S1) with additional sequences obtained for a subset of the WA samples used in microsatellite genotyping (n = 146). New specimens were sequenced using primers L15999M and H16498M [31] with amplifications carried out in 25 µL reactions, which included~10-20 ng of template DNA, 0.5 µM of each primer, 1.25 mM MgCl 2 , 0.2 µM of each dNTPs, 0.05 µL of 10% BSA, 1× Reaction Buffer, and one unit of Invitrogen Taq polymerase. Reactions were run under the following conditions: 94 • C for 3 min followed by 35 cycles of denaturation 94 • C for 30 s, annealing at 55 • C for 45 s, and extension at 72 • C for 60 s, with a final extension of 72 • C for 7 min. Products were visualised on 1% agarose gels and Sanger-sequenced at a commercial service (Australian Genome Research Facility). Sequences were edited and aligned using SEQUENCHER v5.2 (Gene Codes Corporation, Ann Arbor, MI, USA) with consensus sequences from each of the datasets then aligned using CLUSTAL W with default parameters [32] and sequences trimmed to the same length in MEGA v10.0.5 [33].

Phylogenetic Analysis
Phylogenetic analyses of mtDNA sequences were carried out using maximum likelihood in RAXML v7.0.3 [34], neighbour-joining in MEGA v10.0.5 [33] and Bayesian inference in BEAST v1.10.4 [35]. In RAXML, we applied a single model of evolution, General Time Reversible (GTR) model [36] with unequal variation at sites modelled using a Gamma (G) distribution [37] to the sequence data. Rapid bootstrap analysis was applied to search for the best-scoring maximum likelihood (ML) tree. In MEGA, a neighbour-joining tree was constructed using the proportion (p) of nucleotide site changes as a distance estimate and a Gamma (G) distribution as the rate among sites. The robustness of the nodes in the tree was assessed by 1000 bootstrap replicates. For BEAST analysis, we ran JMODEL-TEST v2.1.10 [38] to identify the most likely substitution model for our dataset. In BEAST, we used HKY+I+G model as inferred by the previous step with an uncorrelated lognormal relaxed clock. Four independent runs of 50 million generations and sampling every 5000 generations were carried out. The log files were analysed in TRACER v1.7.1 [39] to ensure stationarity had been reached by visually inspecting the traces and each log file had the effective sample sizes >200 for all parameters. Tree files were combined in LOGCOM-BINER v1.10.4 (BEAST package) with the burn-in period set to 2500 applied to each tree. The combined tree file was annotated using TREEANNOTATOR v1.10.4 (BEAST package). ML, neighbour-joining, and Bayesian trees were visualised using functions from the R packages "APE" [40] and "PHYTOOLS" [41] in R version 3.6.2 [42] using I. macrourus as an outgroup.

Population Level Analysis
To investigate genetic relationships between Isoodon from different clades, principal component analysis was carried out for both nuclear and mitochondrial DNA using the "ADE4" package [43] and the "ADEGENET" package [44] in R based on Euclidean genetic distances using the dudi.pco command. We did not scale allele or haplotype frequencies and left missing genotypes as is. The first two principal components were retained.
Multiple methods were used to infer population structure between and within major geographic regions (WA, SA, Vic) from nuclear DNA. Firstly, we used a Bayesian clustering in STRUCTURE 2.3.4 [45] based on a model that assumed admixture of ancestry and correlated allele frequencies. In each analysis, individuals were assigned a membership coefficient, which is the fraction of the genome with membership to a particular genetic cluster. Ten independent runs were performed using 500,000 iterations, with a burn-in period of 50,000 iterations. We compared the likelihood values for different K values (1)(2)(3)(4)(5)(6)(7)(8)(9)(10)(11)(12)(13)(14)(15)(16)(17)(18)(19)(20) and selected K based on the largest decrease in Delta K value [46,47]. Cluster membership coefficients were permuted over multiple runs using the Greedy algorithm in CLUMPP to obtain a mean across replicates [48]. Secondly, we used spatially explicit Bayesian clustering implemented in TESS 2.3.1 [49,50]. TESS differs from STRUCTURE in that it seeks genetic discontinuities in continuous populations and estimates individual admixture proportions using spatial prior information [49]. We used the admixture model with 50,000 sweeps and 10,000 burn-in for K ranging from 2-20 and 10 replicates per K. The optimal K was chosen where the value of the deviance information criterion (DIC) stabilises. The mean estimated cluster membership coefficient of the chosen K in TESS was obtained using CLUMPP as previously. Thirdly, unlike STRUCTURE and TESS, discriminant analysis of principal components (DAPC) does not assume Hardy-Weinberg equilibrium of clusters. DAPC is a multivariate method that partitions genetic variation into principal components to find groups using K-means that minimize within group variation [51]. We ran the find.cluster command using the ADEGENET package in R. We used the number of components that allowed 90% of cumulative variance to be retained. We selected the optimal K based on a combination of K with the lowest Bayesian information criterion (BIC) and after the largest decrease in BIC [51]. In resolving population structure, we undertook a hierarchical approach, firstly analysing all geographic regions together to identify the broadest genetic structuring, then successively removing the most distinctive cluster and repeating analyses to identify sub-structuring and hence putative conservation units down to the local population scale. Note that for regional analyses, we included the Mt Burr population, located in the far south-eastern South Australia (Figure 1), with other Victorian populations, as this population was shown to cluster with I. o. obesulus Vic [12,20].
Analysis of molecular variance (AMOVA) of mitochondrial DNA was applied in ARLEQUIN version 3.0 [52] to identify the distribution of genetic variation between and within hierarchical clusters identified by STRUCTURE. The locus-by-locus variance component was estimated from 16,000 permutations. Genetic divergences between pairs of population samples were quantified for nuclear DNA using Jost's [53] estimate of differentiation (D est ) using GENALEX version 6.5 [54]. While pairwise φ ST values for mitochondrial DNA were quantified using ARLEQUIN version 3.0 with 16,000 permutations [52].
To test whether the observed haplotype patterns within Isoodon spp. populations indicate recent bottlenecks or, conversely, demographic expansion, we used Tajima's D and Fu's F S neutrality tests and mismatch distribution carried out in ARLEQUIN version 3.0 [52]. In general, Tajima's D and Fu's F S neutrality tests were used to assess whether the data were consistent with the population being at neutral mutation-drift equilibrium [55,56]. Departure of Tajima's D or Fu's F S from zero, either negative or positive, indicates that the population of interest is deviated from the assumption of neutrality and/or equilibrium. If D or F S < 0, the population size may be increasing or an indication of purifying selection. If D or F S > 0, the population size may be decreasing (or bottleneck) or an indication of overdominant selection. We used 1000 simulated samples to calculate the level of significance for both tests. Mismatch distribution employs the distribution of pairwise differences between haplotypes from which spatial population expansion can be estimated. Multi-modal mismatch distribution testifies the population's demographic equilibrium or subdivided populations, while uni-modal suggests recent demographic expansion [57,58], or range expansions with a high level of gene flow [59,60]. We used 1000 bootstrap pseudoreplicates to test the model of spatial expansion by comparing the sum of squared deviations (SSD) between observed and simulated data. The Harpending's raggedness index (r) was also used to test for deviation from unimodality. Smaller values indicate a better fit. The significance of the estimate was also obtained from the corresponding p value.

Genetic Variation Analysis
We assessed the genotypic error rate of microsatellite scores between laboratories by calculating the allele-and locus-specific genotypic error rates [61]. We then tested for Hardy-Weinberg equilibrium (HWE) among loci for each population with sample size ≥10 using GENALEX version 6.51 [54]. We tested for the presence of null alleles at each locus using MICROCHECKER [62]. Estimates of the allelic richness (an estimate of the number of alleles per locus corrected for sample size), number of alleles, gene diversity, inbreeding coefficient (F IS ) were calculated using FSTAT version 2.9.3.2 [63]. The significant deviation of F IS values from Hardy-Weinberg Equilibrium was determined by randomization tests in FSTAT. Differences in gene diversity and allelic richness among regions were statistically tested using a non-parametric Friedman's test with locus as a blocking factor and post-hoc multiple comparisons made using the Conover test with Bonferroni correction in the R package "PMCMR" [64]. Mitochondrial DNA diversity was quantified by calculating the number of haplotypes, gene diversity, and nucleotide diversity in ARLEQUIN version 3.0 [52].

Genetic Data Quality Assessment
Across 218 mtDNA sequenced samples, 96 polymorphic sites, including sites with gaps, were found in the 497 bp D-loop region, giving a total of 77 haplotypes. Across the microsatellite dataset, the amplification success rate was 0.926 per locus. The allele-specific genotypic error rate from 15 re-genotyped samples was 0.107 (0.030 due to allelic drop-out and 0.077 due to false alleles). All loci deviated from HWE at least in one location, but none consistently deviated from HWE across all locations. Likewise, MICROCHECKER detected putative null alleles in all loci, but they were not consistent in all locations. All loci were retained for further analysis.

Phylogeographic Structure across Southern and Western Australia
Phylogenetic analysis of mtDNA using Bayesian methods resolved three major evolutionary groups within our sample (Figure 2). At the highest level, I. o. obesulus from Victoria was distinguished from the remaining SA and WA entities with very high support (PP 1.0; Figure 2). Additional sampling helped us to further resolve two evolutionary groups within the latter clade of southern and western Isoodon, with northern WA taxa (I. a. auratus, I. a. barrowensis), for the most part, distinct from southern WA (I. fusciventer, PP 0.93). However, mtDNA haplotypes from SA I. obesulus were polyphyletic, occurring across both clades. While I. obesulus on Kangaroo Island and Fleurieu Peninsula clustered with southwest WA I. fusciventer, I. obesulus in Mount Lofty Ranges was paraphyletic in both I. fusciventer and I. auratus clades, although with low internal posterior support (PP < 0.90, Figure 2). Similarly, a small number of I. auratus and I. fusciventer haplotypes were placed within opposing clades, again with low support (PP < 0.90, Figure 2). Maximum likelihood and neighbour-joining trees showed a similar pattern of phylogenetic clustering as the Bayesian tree with three major evolutionary groups formed, however, with differing relationships amongst clades ( Figures S1 and S2). In both ML and NJ trees, SA I. obesulus were polyphyletic across I. fusciventer and I. auratus in a similar pattern to the Bayesian tree.
Nuclear (microsatellite) markers further confirmed support for the distinction of Mount Lofty Range, Fleurieu Peninsula, and Kangaroo Island I. obesulus in South Australia from Mount Burr in south-east SA and Victorian I. obesulus, with two main clusters (Mount Burr/Victoria I. obesulus and SA/WA Isoodon) evident in the PCoA (Figure 3a). This latter cluster indicated a geographic cline and genetic overlap from South Australia (I. obesulus), southern Western Australia (I. fusciventer) to northern Western Australia (I. auratus) (Figure 3b), which was similarly supported by PCoA analysis of mtDNA (Figure 3c,d).
At the highest level, all three clustering analyses (Structure, TESS, and DAPC) resolved K = 3 genetic clusters supporting the clear distinction of Mount Burr and Victorian I. o. obesulus from SA/WA (Figure 4a and Figure S3), but which failed to resolve distinction between the SA and WA taxa as currently recognised (i.e., I. auratus, I. fusciventer, SA I. obesulus). SA and northern WA, the two most geographically distant populations, formed the basis of two additional genetic clusters with a pattern of admixture occurring in southern WA (Figure 4a). This admixture pattern was also retained when Victorian I. o. obesulus was removed from analyses ( Figure 4b and Figure S4).
Structuring within WA taxa based on microsatellite data was only revealed upon hierarchical analysis. With SA samples removed, structure analysis recovered K = 3 populations representing southern WA (I. fusciventer), Barrow Island (I. a. barrowensis), and the Kimberley (I. a. auratus) (Figure 4c). These clades are largely represented in phylogenetic analyses with subspecies I. a. barrowensis forming a distinct clade in maximum likelihood and neighbour-joining trees ( Figures S1 and S2) but nested within a group comprising other Kimberley mainland samples in Bayesian analyses (Figure 2).    (Figures 4a and S3), but which failed to resolve distinction between the SA and WA taxa as currently recognised (i.e., I. auratus, I. fusciventer, SA I. obesulus). SA and northern WA, the two most geographically distant populations, formed the basis of two additional genetic clusters with a pattern of admixture occurring in southern WA (Figure 4a). This admixture pattern was also retained when Victorian I. o. obesulus was removed from analyses (Figures 4b and S4).
Structuring within WA taxa based on microsatellite data was only revealed upon hierarchical analysis. With SA samples removed, structure analysis recovered K = 3 popu-

Population-Level Analysis
At the population level, we detected substantial genetic sub-structuring within each geographic region, with the exception of southern WA (Figure 4d). Three clusters were identified in northern WA separating Kimberley islands, Kimberley mainland and Barrow Island (Figure 4d). We detected multiple admixed clusters in southern WA and the modal K value varied between three and four clusters ( Figure S3); in each, however, inland southwest WA was consistently distinct from the Perth population. In South Australia, Kangaroo Island was differentiated from Mount Lofty Ranges, but Fleurieu Peninsula was either assigned to Kangaroo Island (STRUCTURE), admixed between the two (TESS), or assigned to its own cluster (DAPC). In Victoria, the Grampians was differentiated from Mount Burr in south-east SA, while Lower Glenelg either assigned to the same cluster as Mount Burr (STRUCTURE and TESS) or formed its own cluster (DAPC) ( Figure S3).
AMOVA also supported geographic structuring of mtDNA haplotypes when each region was analysed separately and revealed that a large proportion of genetic variation

Population-Level Analysis
At the population level, we detected substantial genetic sub-structuring within each geographic region, with the exception of southern WA (Figure 4d). Three clusters were identified in northern WA separating Kimberley islands, Kimberley mainland and Barrow Island (Figure 4d). We detected multiple admixed clusters in southern WA and the modal K value varied between three and four clusters ( Figure S3); in each, however, inland southwest WA was consistently distinct from the Perth population. In South Australia, Kangaroo Island was differentiated from Mount Lofty Ranges, but Fleurieu Peninsula was either assigned to Kangaroo Island (STRUCTURE), admixed between the two (TESS), or assigned to its own cluster (DAPC). In Victoria, the Grampians was differentiated from Mount Burr in south-east SA, while Lower Glenelg either assigned to the same cluster as Mount Burr (STRUCTURE and TESS) or formed its own cluster (DAPC) ( Figure S3).
AMOVA also supported geographic structuring of mtDNA haplotypes when each region was analysed separately and revealed that a large proportion of genetic variation occurred between sampling sites, ranging from 27.5-73.3% (p < 0.001 in all cases, Table 1). In South Australia, 73.3% of genetic variation was found among populations, while approximately half of the genetic variation was found among populations in northern WA (48.3%) and Victoria (48.4%). In southern WA, more genetic variation was found within population (72%) indicating greater common shared genetic variation among populations than other regions (Table 1), consistent with patterns identified in cluster analyses. Pairwise population D est values based on microsatellite and φ ST values based on mtDNA data indicated a substantial genetic divergence between Isoodon in different regions ( Table 2). The highest differentiation was recorded between SA and Victoria samples with Substantial genetic differentiation was also observed amongst populations within geographic regions, particularly of Kangaroo Island from remaining SA populations (Table 2).

Population Genetic Diversity
Overall, mtDNA diversity was high and comparable between southern WA, northern WA, and SA Isoodon populations (Table 3). Across all diversity metrics, Victorian populations had low diversity, with lower haplotype diversity than all other populations (Table 3) and similarly low nuclear diversity to SA (H = 0.44-0.61 and 0.42-0.60, respectively; Friedman's test, p < 0.05) indicating these populations may be in genetic decline. Multi-locus F IS values had significantly positive F IS values indicating a possible Wahlund effect (samples were sourced from multiple populations) or non-random mating patterns within all sites with the exception of Kimberley mainland (randomization tests, p < 0.0042, Table 3). Kimberley mainland had the highest number of unique haplotypes, yet most haplotypes occurred at low frequencies. This pattern is concordant with evidence of recent population expansion according to the unimodal mismatch distribution and correlation with the spatial expansion model (Tajima's D = −1.87, Fu's Fs = −5.90, p < 0.05 in both cases; Figure S5). Similarly, the Kimberley mainland retains high microsatellite diversity relative to other populations.

Discussion
In conservation biology, a taxonomic classification system that closely reflects the biology of a species is essential to guide conservation management actions and increase the chance of success. Incorrect designation of species and subspecies could lead to inappropriate management decisions that can have detrimental effects such as outbreeding depression, maladaptation, or separating populations when genetic augmentation could benefit inbred populations [3,5,67]. Through increased spatial sampling and addition of hyper-variable nuclear markers, our study aimed to provide further resolution to the genetic relationships of Isoodon species in Western Australia and South Australia, where current taxonomy does not reflect reported phylogeographic relationships [12]. At the broadest level, increased sampling enabled us to resolve two well-supported mtDNA clades representing WA entities I. auratus (northern WA) and I. fusciventer (southern WA); however, I. obesulus in South Australia remained polyphyletic with SA haplotypes represented in both WA clades, suggesting that incomplete lineage sorting may be confounding phylogenetic analyses across this group. Further highlighting the close relationship amongst these three taxa, analyses of nuclear microsatellite markers indicated a geographic cline stretching from northern WA to South Australia rather than resolving distinct taxonomic entities, with sub-structuring amongst taxa only becoming evident in hierarchical clustering analyses. In all analyses, we detected highly localised patterns of mitochondrial and nuclear genetic diversity indicating limited contemporary gene flow. Taken together, our analyses indicate a complex phylogeographic history for Isoodon sp. in southern and western Australia. We discuss identification of conservation units below, as well as some of the underlying processes that may have contributed to the observed results.

Distinction between I. fusciventer, I. auratus, and I. obesulus (SA)
The relationship of WA taxa I. auratus and I. fusciventer has long been queried, with multiple mtDNA studies over the past 20 years reporting the taxa as polyphyletic, evidence used to suggest the two species be synonymised [12,15,21]. This observation has been controversial given the striking morphological differences between I. auratus and I. fusciventer, with I. fusciventer more similar in size and pelage colour to I. obesulus than I. auratus (Figure 2; [19]), although several case-studies have indicated morphological variation in WA bandicoots can be strongly environmentally determined [68][69][70][71].
Here, with expanded geographic sampling of both I. auratus (northern WA) and I. fusciventer (southern WA), we largely resolved two monophyletic groups in mtDNA (putative ESUs) with strong Bayesian support (PP 0.93). However, a very small number of haplotypes from I. fusciventer were clustered with I. auratus and vice versa, rendering these groups polyphyletic. Further, to add complexity to this relationship, haplotypes from South Australian I. obesulus were polyphyletic with I. auratus and I. fusciventer, with Kangaroo Island and Fleurieu Peninsula samples clustered within I. fusciventer and Mount Lofty Ranges samples polyphyletic between I. auratus and I. fusciventer. Even though SA mtDNA haplotypes were unique, they were most closely related to haplotypes of either I. auratus and I. fusciventer indicating a relictual connection between these lineages. For species undergoing allopatric speciation, coalescent theory predicts that the process is detected in molecular phylogenies through the progression from a pattern of polyphyly to paraphyly to reciprocal monophyly, the timing of which is a function of effective population size with large populations taking longer to achieve reciprocal monophyly than small ones [72,73]. Isoodon auratus, in particular, was historically abundant and widespread across much of central and western Australia, possibly extending into New South Wales [17] before massive range collapse following European settlement, and I. auratus and I. fusciventer currently retain large population sizes relative to other eastern I. obesulus [28]. The close relationships of various SA I. obesulus, I. auratus, and I. fusciventer mtDNA haplotypes likely reflects historical connection amongst these taxa and suggests that time since geographic isolation has been insufficient to allow complete lineage sorting at this locus. This is perhaps unsurprising given the recent diversification of bandicoots in Australia within the late Pleistocene (3.1 Mya), with more recent speciation between I. auratus and I. fusciventer estimated at approximately 1.5 Mya [19]. An opportunity to genetically sample extinct populations from the previous range of these species (e.g., via museum specimens or sub-fossil deposits [74]) could help to more fully illuminate the historical demography and connectivity across this group.
Analyses of nuclear (microsatellite) allele frequency data showed inconsistency with mtDNA results for the relationship of SA I. obesulus with WA taxa. PCoA of microsatellite allele frequencies indicated a geographic cline ranging from northern WA, through southern WA to South Australia reflecting I. fusciventer as a connecting population according to their geographic relationship. In contrast, PCoA of mtDNA showed SA I. obesulus as closely related to both I. auratus and I. fusciventer with the latter two species at the ends of the cline. Clustering analyses of microsatellite data (TESS, Structure, DAPC) consistently resolved K = 2 clusters separating WA (I. auratus, I. fusciventer) from SA (I. obesulus), albeit with a pattern of admixture between the two clusters occurring in southern WA (Structure, TESS) that is often indicative of a pattern of genetic isolation by distance [75,76]. Such a pattern of IBD in microsatellite data may reflect greater male-biased dispersal, whereas restricted female dispersal has led to greater structuring in mitochondrial DNA that is maternally inherited (e.g., [74]). Hierarchical structuring present within the WA cluster was revealed upon stepwise removal of divergent groups in clustering analyses, which then resolved K = 3 clusters within WA representing I. a. auratus (Kimberley), I. a. barrowensis (Barrow Island), and I. fusciventer (southern WA).
Conflicting signals from mitochondrial and nuclear DNA analyses in our study thus make diagnosis of ESUs difficult under Moritz's genetic criteria of reciprocal monophyly, particularly in relation to SA I. obesulus. Under the less restrictive Adaptive Evolutionary Conservation concept of Fraser and Bernatchez [8], requiring only demonstration of highly restricted gene flow amongst lineages with a fundamental goal to preserve adaptive genetic variance within species, our data suggest I. auratus, I. fusciventer, and SA I. obesulus represent separate ESUs, with additional hierarchical substructure identified within I. auratus separating I. auratus auratus (northern WA, Kimberley) and I. a. barrowensis (northern WA, Barrow Island). The high genetic distinction of I. auratus, I. fusciventer and SA I. obesulus in mitochondrial and nuclear DNA is reflective of highly reduced gene flow, likely as a result of their allopatric distribution, and is in line with previous results from allozyme data [12] and morphological differences [15,19,22,23]. Multiple lines of evidence suggest these are distinct evolutionary lineages that warrant separate taxonomic status. While SA I. obesulus appeared related to both I. fusciventer and I. auratus in mtDNA, cluster analysis of nuclear data indicated that SA I. obesulus is more genetically distinct from the remaining WA taxa, which may warrant recognition of this entity at species level, particularly if I. auratus and I. fusciventer are retained as such. Allozyme data presented in Cooper et al. [12] shows a much closer relationship of I. fusciventer and SA I. obesulus, clearly distinct from I. auratus, which is consistent at least with the superficial morphology of the two entities. Further genetic and morphological investigation is required to clarify taxon boundaries amongst these entities, particularly as the SA population is highly restricted and declining and likely to require urgent conservation attention.

Lack of Clarity on the Taxonomic Relationship of SA I. obesulus from Current Genetic Data
Previous studies of the phylogenetic relationships of Isoodon bandicoots have relied on small numbers of samples and/or limited numbers of phylogenetic markers (mitochondrial or nuclear sequencing markers) [12,15,20,21]. Here, with expanded sampling of problematic taxa (I. auratus, in particular) and with the addition of high mutation rate nuclear markers (microsatellites), we aimed to provide greater resolution to the relationships between purported taxa in the "southern and western" clade of closely related Isoodon bandicoots identified in Cooper et al. [12]. Our analyses recovered ESUs and sub-hierarchical ESUs that are largely reflective of current taxonomic designations, although ambiguity remains over the taxonomic identification of SA I. obesulus. There is growing recognition that multi-locus and multi-species coalescent approaches are required to overcome the limitations of current species delimitation methods, particularly when single gene trees are poorly resolved or in conflict with the actual species tree, as may arise due to ancestral lineage sorting, hybridisation, and a range of other population genetic processes [77,78]. Here, we find a disconnect in the results obtained for mitochondrial and nuclear DNA analyses for SA I. obesulus that may be attributed to a mismatch in the tempo of mutation rates in the markers we have used to the evolutionary processes occurring in the group. For example, mtDNA analysis is impacted by incomplete lineage sorting, and microsatellites may have mutation rates that are too rapid to reflect speciation patterns, as they largely reflect contemporary population genetic processes. Genomic sequencing approaches, for example, targeted exon capture [79] enable concurrent sampling of multiple, genome-wide markers with a range of mutation rates and could provide sufficient resolution to resolve speciation patterns across southern and western Australia. Given the relatively recent pattern of speciation and massive range collapse experienced by these taxa investigation of the demographic history of the group will be particularly fascinating, especially as the current distribution boundaries of Isoodon populations are coincident with biogeographic barriers that have been identified elsewhere (e.g., Nullarbor barrier, Murravian barrier; [80]).

Genetic Diversity and Population Substructure
Bandicoot species, particularly across southern Australia, are highly impacted by habitat fragmentation, habitat loss, and declining population size. Here, we found that both Victorian and South Australian populations of I. obesulus showed consistently low genetic diversity across nuclear markers relative to WA populations, and that Victorian populations had low mitochondrial diversity. In contrast, the Kimberley mainland population of I. auratus, which is less impacted by anthropogenic fragmentation effects, retains high genetic diversity with signals suggesting population expansion in the region.
Consistent with this, we also found evidence of restricted gene flow with very high numbers of unique local mtDNA haplotypes, significant and high pairwise D est and φ ST values, and a significant proportion of genetic variation apportioned between, rather than within, populations. We also detected significant genetic sub-structuring within each region via cluster analyses, although populations in southern WA were more admixed than elsewhere. Such a pattern is consistent with the expected limited dispersal range of bandicoots (~3 km; Li et al. [81,82]), but connectivity is also likely moderated by the amount (and fragmentation) of native vegetation in the landscape [13].

Recommendations for Conservation Management
Our mtDNA analysis, haplotype uniqueness, and significant divergence among populations in different regions suggested that I. auratus, I. fusciventer, and I. obesulus in South Australia should be managed as separate units, and that SA I. obesulus requires urgent consideration to assess its taxonomy and conservation status. Several populations included in this study showed evidence of genetic decline, particularly those in Victoria and South Australia. Ongoing habitat fragmentation and degradation in these regions may be contributing to limited gene flow amongst patches and subsequently genetic drift. Improving landscape-scale connectivity will be of high priority for these populations with evidence that low, shrubby vegetation can provide a structural habitat for bandicoots to avoid predation as well as connectivity, even in degraded environments [13,14,83]. Management of introduced predators (foxes, cats) and habitat quality improvement will also be critical factors [84][85][86]. As one of Australia's important digging mammals, bandicoots are also targeted for translocation to restore ecological function in ecological restoration projects and to secure "safe haven" sites to improve their conservation status [87][88][89][90][91][92][93]. Golden bandicoots (I. auratus) are of particular importance for arid zone translocations with those undertaken to date considered successful [30,70]. There is evidence that bandicoots show morphological adaptation to local conditions [68][69][70], which has been shown to have a genetic, as well as an environmental component [71], suggesting that translocation to similar habitats is important. However, all translocations to date have been in north-western Australia and have involved the Barrow Island population of I. auratus, which we demonstrate here to have lower genetic diversity (being an island population) than mainland Kimberley I. auratus. It may be prudent when planning future translocations to consider whether it is feasible to source animals from Kimberley populations to ensure genetic diversity within the species is maintained across multiple sites to reduce risks to further genetic decline.