Multi-Species Phylogeography of Arid-Zone Sminthopsinae (Marsupialia: Dasyuridae) Reveals Evidence of Refugia and Population Expansion in Response to Quaternary Change

Historical population contraction and expansion events associated with Pleistocene climate change are important drivers of intraspecific population structure in Australian arid-zone species. We compared phylogeographic patterns among arid-adapted Dasyuridae (Sminthopsis and Planigale) with close phylogenetic relationships and similar ecological roles to investigate the drivers of phylogeographic structuring and the importance of historical refugia. We generated haplotype networks for two mitochondrial (control region and cytochrome b) and one nuclear (omega-globin) gene from samples distributed across each species range. We used ΦST to test for a genetic population structure associated with the four Pilbara subregions, and we used expansion statistics and Bayesian coalescent skyline analysis to test for signals of historical population expansion and the timing of such events. Significant population structure associated with the Pilbara and subregions was detected in the mitochondrial data for most species, but not with the nuclear data. Evidence of population expansion was detected for all species, and it likely began during the mid-late Pleistocene. The timing of population expansion suggests that these species responded favorably to the increased availability of arid habitats during the mid-late Pleistocene, which is when previously patchy habitats became more widespread. We interpret our results to indicate that the Pilbara region could have acted as a refugium for small dasyurids.


Introduction
In the northern hemisphere, environmental change during the Pleistocene elicited range contraction into southern refugia during repeated glacial maxima (GM) [1], or, expansion in the case of cold-tolerant species [2]. In non-glaciated regions, the absence of widespread ice-sheets has led to idiosyncratic rather than leading-edge phylogeographic patterns, resulting in multiple, species-specific refugia distributed in climatically altered landscapes [3]. In arid regions, habitat shifts were associated with the activation of mobile dune systems as vegetation patterns accommodated climate change and resulted in a complex array of species responses to change during the Quaternary [4], including in Australian deserts [5,6]. comparative phylogeography as they have extensive distributions across a variety of arid habitats, allowing for comparisons across multiple bioregions and habitat types.
In this study, we explore the processes that generate phylogeographic structuring in arid occurring dasyurids, with a particular focus on the role of the Pilbara as an evolutionary refugium. We selected six co-distributed species of Sminthopsis and Planigale that occur in the Pilbara that have shared evolutionary histories and similar life history strategies, but different dispersal potential and habitat preferences. Using DNA sequence data to test for restricted gene flow across the boundary of the Pilbara, and between the IBRA subregions within the Pilbara, we expected higher levels of genetic diversity to occur within the Pilbara if that region had acted as an historic refugium for small dasyurids. Additionally, we test for historical demographic responses to past environmental change within this confamilial group to identify responses driven by GM habitat contractions between species with different dispersal abilities and divergent habitat specificities, particularly since the last GM.

Study Species
Six small dasyurid species (subfamily Sminthopsinae), Sminthopsis longicaudata [33], Sminthopsis macroura [34], Sminthopsis ooldea [35], Sminthopsis youngsoni [36], Planigale sp. 1, and Planigale sp. 2, were selected for this study. The two unnamed Planigale species have previously been determined from molecular data [37] and are morphologically distinct (K.P. Aplin pers. comm.). As it is currently recognized, S. macroura represents a species complex with potentially five named taxa that were all synonymized by Archer [38]. Molecular data have shown that at least three distinct species (S. macroura, S. froggatti, and S. stalkeri) should be recognized, as well three lineages within the 'true' S. macroura, one of which occurs in the northwest of Western Australia (WA) [39] and is somewhat geographically isolated (Figure 1b). The study species each have different habitat preferences and exhibit a range of potential dispersal abilities (Table 1). Both Planigale species are endemic to the greater Pilbara region [40], while the four Sminthopsis occur in the Pilbara and also in the adjacent arid regions (Figure 1b). Table 1. Summary of attributes of the six mammal species chosen for this study (data from [29,32]).

Species
Habitat/Specializations Reproduction Dispersal Potential

Sminthopsis longicaudata
Specialist. Exposed rock/stony soils, flat-topped hills, lateritic plateaus, sandstone ranges, and breakaways. Striated foot pads for climbing rocky surfaces. Potentially poor, perhaps only on clay substrates  [7]), location and direction of major dune systems [41], approximate location of Last Glacial Maxima (LGM) shoreline [5], and the Pilbara bioregion outlined in black, inset map shows digital elevation model and detail of Interim Biogeographic Regionalization for Australia (IBRA) within the Pilbara.  [7]), location and direction of major dune systems [41], approximate location of Last Glacial Maxima (LGM) shoreline [5], and the Pilbara bioregion outlined in black, inset map shows digital elevation model and detail of Interim Biogeographic Regionalization for Australia (IBRA) within the Pilbara. (b) Distribution maps of the six species (gray points indicate vouchered specimens in Australian museums), and tissues sequenced in this study (black points). Non-WA records were accessed from Atlas of Living Australia.

DNA Extraction and Amplification and Sequence Alignment
Total genomic DNA was obtained from 833 samples held at the WA Museum, the WA Department of Biodiversity, Conservation and Attractions, and the Australian Biological Tissue Collection at the South Australian Museum (see Table S1 Supplementary Materials). Biological material was preserved in 70-100% ethanol or cryo-frozen at −75 • C, and it included liver, muscle, heart, ear, tail, and toe. For DNA extractions, approximately 5-25 mg of biological material was used in Qiagen DNeasy tissue and blood kits according to the manufacturers' instructions, and DNA was eluted into 60-120 µL tris buffer. Several mitochondrial and nuclear loci were screened for suitable degrees of intraspecific variation across the study species and control region (CR), cytochrome b (cytb), and omega globin (ω-globin) were chosen, GenBank accession numbers, including those for non-suitable loci sequences are in Table S1 (MT366215-MT454794). Primers used for CR, cytb, and ω-globin are listed in Table S2.
All regions were amplified via polymerase chain reaction (PCR) as described in [42], following the PCR cycling conditions of [43]. For the amplification of CR for S. youngsoni, the annealing temperature was increased from 50 to 57 • C due to the presence of large repeat regions at the 5 end (see [44]). DNA purification and bidirectional sequencing was carried out at the Australian Genome Research Facility, Perth WA. Forward-reverse assemblies, quality control, and alignment of sequences was performed in GENEIOUS v.9.1.7 [45], as described in [42]. Final sequence length for CR ranged from 560 to 800 bp, and ω-globin ranged from 770 to 790 bp ( Table 2). Only approximately the first 1000 bp of cytb could be amplified for Planigale sp. 1, Planigale sp. 2, and WA S. macroura samples. Table 2. Molecular diversity indices for the concatenated mtDNA (CR + cytb), and the phased nuclear (ω-globin) gene alignments of the six species in this study. Indices are also calculated for species within and outside of the Pilbara where appropriate. Number of samples (n), sequence length in base pairs excluding alignment gaps (L), number of haplotypes or alleles (Hn), haplotype diversity (H), nucleotide diversity (π), haplotypic richness under rarefaction (Hr).

Species/Lineage Boundaries in S. macroura
To determine the geographic extent of the WA lineage of S. macroura, previously identified from three samples by Blacket et al. [39], a maximum likelihood tree was built on the concatenated DNA sequences of a subset of S. macroura samples using the GENEIOUS RAXML plug-in [46]. Two additional samples, Planigale sp. 1 (WAM M57732) and S. longicaudata (WAM M52415), were used as outgroups. The model selecting program PARTITIONFINDER v. 1.1.1 [47] was used and selected GTR + G, with the alignment partitioned by gene and codon position for cytb, and run for 1000 bootstrap replicates. Based on the tree analysis ( Figure 2), only the WA lineage was considered for the subsequent phylogeographic and population-level analysis. additional samples, Planigale sp. 1 (WAM M57732) and S. longicaudata (WAM M52415), were used as outgroups. The model selecting program PARTITIONFINDER v. 1.1.1 [47] was used and selected GTR + G, with the alignment partitioned by gene and codon position for cytb, and run for 1000 bootstrap replicates. Based on the tree analysis (Figure 2), only the WA lineage was considered for the subsequent phylogeographic and population-level analysis.

Phylogeographic Structure Analysis
To investigate phylogeographic patterns, TCS neighborhood-joining haplotype networks were built using the software program POPART [48]. POPART excludes gaps and ambiguous bases from the analysis, so alignments with difficult to resolve indel regions (e.g., the 5′ end of CR) were cleaned using the online program GBLOCKS v.0.19b, using default settings [49] (available through the website phylogeny.fr; [50]). Nuclear sequences were phased using the statistical haplotyping program PHASE [51], and input files were created using the online web-tool seqPHASE

Phylogeographic Structure Analysis
To investigate phylogeographic patterns, TCS neighborhood-joining haplotype networks were built using the software program POPART [48]. POPART excludes gaps and ambiguous bases from the analysis, so alignments with difficult to resolve indel regions (e.g., the 5 end of CR) were cleaned using the online program GBLOCKS v.0.19b, using default settings [49] (available through the website phylogeny.fr; [50]). Nuclear sequences were phased using the statistical haplotyping program PHASE [51], and input files were created using the online web-tool seqPHASE (http://seqphase.mpg.de/seqphase/; [52]) where the allele combinations with probability scores above 0.7 were chosen for the analysis (for S. ooldea and S. longicaudata, we allowed scores from 0.5 to be included to ensure all variation was captured). The final length of the cleaned and trimmed alignments used for analysis is given in Table 2 as L. Then, statistical parsimony TCS networks [53] were built for each loci and the concatenated mtDNA alignments of each species (Figure 3).
Genes 2020, 11, x FOR PEER REVIEW 7 of 21 (http://seqphase.mpg.de/seqphase/; [52]) where the allele combinations with probability scores above 0.7 were chosen for the analysis (for S. ooldea and S. longicaudata, we allowed scores from 0.5 to be included to ensure all variation was captured). The final length of the cleaned and trimmed alignments used for analysis is given in Table 2 as L. Then, statistical parsimony TCS networks [53] were built for each loci and the concatenated mtDNA alignments of each species ( Figure 3).

A Priori Groupings for Geographical Analysis
We tested for significant population structure associated with the Pilbara region and major subregions within the Pilbara among the study species. The IBRA subregions contain distinct geology, landscape features, and vegetation that are important for determining species occurrences and shaping genetic structure [22]. To calculate the equilibrium genetic differentiation summary statistics, such as ΦST, we clustered individuals that were continuously distributed across the landscape into subpopulations. To test for population differentiation across the Pilbara boundary, we grouped each species into "Pilbara" (i.e., those samples that fell within the Pilbara bioregion), and "non-Pilbara" populations (i.e., those outside of the Pilbara). This resulted in two populations, except for S. youngsoni, which had three populations, as the distribution of this species is disjunct across the south of the Pilbara (Figure 1b). The two Planigale species were excluded from this level of analysis due to low sample numbers outside of the Pilbara. To test for population differentiation within the Pilbara, samples were grouped into subpopulations based on the four IBRA subregions-Roebourne, Chichester, Fortescue, and Hamersley (Figure 1a), as these subdivisions capture the geographic complexity of the two major range systems and the low-lying Fortescue river valley, which could impact the dispersal of small mammals. S. ooldea and S. longicaudata were not included in this level of analysis, due to low numbers of samples within the Pilbara that were

A Priori Groupings for Geographical Analysis
We tested for significant population structure associated with the Pilbara region and major subregions within the Pilbara among the study species. The IBRA subregions contain distinct geology, landscape features, and vegetation that are important for determining species occurrences and shaping genetic structure [22]. To calculate the equilibrium genetic differentiation summary statistics, such as Φ ST , we clustered individuals that were continuously distributed across the landscape into subpopulations. To test for population differentiation across the Pilbara boundary, we grouped each species into "Pilbara" (i.e., those samples that fell within the Pilbara bioregion), and "non-Pilbara" populations (i.e., those outside of the Pilbara). This resulted in two populations, except for S. youngsoni, which had three populations, as the distribution of this species is disjunct across the south of the Pilbara (Figure 1b). The two Planigale species were excluded from this level of analysis due to low sample numbers outside of the Pilbara. To test for population differentiation within the Pilbara, samples were grouped into subpopulations based on the four IBRA subregions-Roebourne, Chichester, Fortescue, and Hamersley (Figure 1a), as these subdivisions capture the geographic complexity of the two major range systems and the low-lying Fortescue river valley, which could impact the dispersal of small mammals. S. ooldea and S. longicaudata were not included in this level of analysis, due to low numbers of samples within the Pilbara that were predominantly found in just one IBRA subregion. Then, these aforementioned groupings were used to run global Φ ST and pairwise Φ ST tests on the cytb alignments for each species (only cytb was used due to the low coverage of CR for S. longicaudata) in the program ARLEQUIN v.3.5.2.2 [54]. To account for different sample sizes between subpopulations, we also used a rarefaction approach to calculate the Haplotypic Richness and the number of Private Haplotypes as implemented in the program HP-Rare [55].

Tests for Historical Population Size Changes
We tested for evidence of historical population expansion among the study species and subpopulations. Molecular variation indices and expansion statistics were calculated in DNASP v.5 [56] on the concatenated mtDNA and nDNA alignments used in the network analysis for each species, or haplo group where appropriate based on the network results. Significance values for Ramos-Onsins R2 and Fu's Fs were obtained under coalescent simulations with segregating sites. Rates and timing of population expansion were investigated using the extended Bayesian skyline model (EBSP) [57] in BEAST2 [58]. Skyline analysis assumes that sampling is from a panmictic population [59]; so, where species had structured populations, identified from the mtDNA network analysis, these were analyzed separately, and only those that successfully converged are shown ( Figure 4). PARTITIONFINDER was used to estimate substitution models for each of the three genes and codons for cytb; however, most models failed to converge when cytb was partitioned by codon, so it was not partitioned for any analyses except for S. ooldea. Site models from PARTITIONFINDER were selected, with the exception of estimating the portion of invariable sites (+I), as it can overcomplicate the model and is not considered biologically meaningful for intraspecific level data [60]. For Sminthopsis species, a substitution rate of 0.0105 substitutions/site/Myr and standard deviation (sd) of 0.001 using a normal prior, and a rate of 0.0095 substitutions/site/Myr (sd of 0.0005) for Planigale species were used to capture the ranges in Krajewski et al. [61], who reported different cytb mutation rates for each genus. A uniform prior was used for CR and ω-globin with substitution rates estimated. The EBSP analysis was run for 100 million generations to ensure that the effective sample sizes of posterior parameter estimates were ≥ 200. The EBSP log files were plotted in R using the 'plotEBSP.R' script [57].

.Relationships within the Sminthopsis macroura Complex
Tree analysis on S. macroura samples supported the distinctiveness of a WA lineage (bootstrap support 100%) and indicated that the distribution of the WA lineage is restricted to northwest and central WA (Figure 2). Subsequent downstream analysis focused only on the WA lineage.

Gene Diversity and Haplotype Distribution
All six species had multiple haplotypes (30-130) and high haplotype diversity (0.88-0.99) for the concatenated mtDNA dataset ( Table 2). Haplotype richness and the number of private haplotypes under rarefaction were highly correlated (correlation coefficient = 0.99) and were lower for the Pilbara compared to non-Pilbara samples for S. ooldea, S. longicaudata, and S. youngsoni, and higher within the Pilbara for S. macroura ( Table 2). Haplotypes that formed clades in the network analysis were grouped together and hereafter referred to as 'haplo groups' (Figure 3). Allele diversity for the phased nuclear data was lower than the mtDNA loci and did not correlate to geographic regions in the network analysis (see Figure S1, Supplementary Material).
The restriction of haplo groups to geographic regions occurred in the mtDNA sequence data for S. ooldea, S. macroura, and slightly for Planigale sp. 1, but no obvious geographically restricted haplo groups were detected for S. longicaudata, S. youngsoni, or Planigale sp. 2. Of the 50 mtDNA haplotypes

Relationships within the Sminthopsis macroura Complex
Tree analysis on S. macroura samples supported the distinctiveness of a WA lineage (bootstrap support 100%) and indicated that the distribution of the WA lineage is restricted to northwest and central WA (Figure 2). Subsequent downstream analysis focused only on the WA lineage.

Gene Diversity and Haplotype Distribution
All six species had multiple haplotypes (30-130) and high haplotype diversity (0.88-0.99) for the concatenated mtDNA dataset ( Table 2). Haplotype richness and the number of private haplotypes under rarefaction were highly correlated (correlation coefficient = 0.99) and were lower for the Pilbara compared to non-Pilbara samples for S. ooldea, S. longicaudata, and S. youngsoni, and higher within the Pilbara for S. macroura ( Table 2). Haplotypes that formed clades in the network analysis were grouped together and hereafter referred to as 'haplo groups' (Figure 3). Allele diversity for the phased nuclear data was lower than the mtDNA loci and did not correlate to geographic regions in the network analysis (see Figure S1, Supplementary Material).
The restriction of haplo groups to geographic regions occurred in the mtDNA sequence data for S. ooldea, S. macroura, and slightly for Planigale sp. 1, but no obvious geographically restricted haplo groups were detected for S. longicaudata, S. youngsoni, or Planigale sp. 2. Of the 50 mtDNA haplotypes for S. ooldea, nine were restricted to the Pilbara (halpo group 1) and separated by 28 mutations from the other samples (Figure 3a). There were 76 mtDNA haplotypes for S. youngsoni, of which 25 only occurred in the Pilbara, and the species had three haplo groups; two of these were separated by 15 mutations and occurred on either side of the Gibson Desert (halpo group 1 and 2; Figure S2), and a third group that was 30 mutations from haplo groups 1 and 2 was distributed throughout WA (haplo group 3, Figure 3b). The WA lineage of S. macroura had 75 of the 100 mtDNA haplotypes located in the Pilbara only, and four haplo groups were identified, though with fewer mutations between each when compared to S. ooldea and S. youngsoni. Of these, haplo group 1 was mostly restricted to the Hamersley subregion; by comparison, the most diverse haplo group (3) had few individuals occurring within the Hamersley, instead being distributed widely throughout the Chichester and Carnarvon regions (Figure 3c; Figure S2). Planigale sp. 1 exhibited four haplo groups (Figure 3d), the larger two (1 and 2) being split across the center of the Pilbara with haplo group 2 occupying the eastern half of the Chichester and Fortescue subregions, which was exclusive to the other halpo groups ( Figure S2). S. longicaudata had 26 haplotypes, of which nine were only found in the Pilbara, most of which were distributed throughout the network (Figure 3e). Planigale sp. 2 had lower nucleotide diversity for all three genes compared to the other five species (Table 2), and the phylogeographic structure in the network analysis did not correlate to the IBRA subregions (Figure 3f).  (Table 3).

Population Expansion Analysis
Population size changes were tested using summary statistics (Table 4) and extended Bayesian skyline analysis (Figure 4). Population expansion was detected for the concatenated mtDNA sequence data for all six species and subpopulations with significant and high Fu's Fs values; however, test scores for Tajima's D and Ramos-Onsins R2 were not significant except for R2 for Planigale sp. 1 (Table 4). There was almost no evidence of population expansion with the nuclear loci, except for S. youngsoni, which was significant for R2 and Fs. Coalescent skyline analyses showed population size increases for all haplo groups examined except S. ooldea haplo group 1 (Figure 4). The onset of population expansion occurred between 150 and 200 kya for S. youngsoni, S. macroura, and Planigale sp. 1, and 100 kya for S. ooldea. Expansion began earlier for S. longicaudata, approximately 330 kya, and Planigale sp. 2 showed a steady increase in size, but large credibility intervals indicate uncertainty with population size estimates earlier than 200 kya. Table 4. Expansion statistics for the concatenated mtDNA (CR + cytb) and the phased nuclear (ω-globin) gene alignments of the six species in this study. Some species were grouped according to mtDNA haplo groups from the network analysis that also successfully converged in the skyline analysis.

Discussion
Late Quaternary climate changes were a major driver of species phylogeographic structure throughout arid Australia, resulting in intraspecific variation from multiple population expansion and contraction events across an array of species [5]. Here, we show that the demographic responses among related dasyurid species to past climate change exhibit similar patterns of population expansion during the Pleistocene, suggesting that their shared evolutionary history and adaptation to a climatically unstable and relatively young landscapes has determined patterns of genetic structure in these species.

Pilbara Refugia
The Pilbara is recognized as a biodiversity hotspot [18,22] and is home to a number of endemic mammal species [40]. We hypothesized that there would be higher genetic diversity within Pilbara populations if the region had acted as a refugium for small dasyurids during GM. However, mtDNA haplotypic richness was only markedly higher in the Pilbara for S. macroura, among the four congeneric species tested, when compared to populations outside the Pilbara. This indicates that the Pilbara may have acted as an evolutionary refugia for the WA lineage of S. macroura in the past, but that other regions within the distributions of S. ooldea, S. longicaudata, and S. youngsoni have been equally as, if not more, important as the Pilbara in preserving haplotype diversity over time. What is more significant, as shown by our results, is the impact of the Pilbara geological boundary as a barrier to gene flow for small dasyurids. We show significant population differentiation associated with that boundary in the majority of species, indicating that most species do not readily disperse across the Pilbara boundary, which is likely due to the discontinuity of habitats over this boundary [22].
We have identified the southeastern edge of the Pilbara craton as an important phylogeographic break between Pilbara and non-Pilbara samples of S. ooldea, despite this species occurring continuously from the Pilbara into central South Australia (Figure 1b). This differentiation may be attributed to the disruption of the Acacia-dominated habitat between the edge of the Pilbara craton and the sandy substrates of the adjacent Little Sandy Desert, since S. ooldea prefers Acacia woodland habitats [62]. The disjunction between the rocky Pilbara craton and the dune systems of the Great and Little Sandy Deserts has been identified as a limiting factor in the dispersal of non-sand-preferring species [63], such as the gecko Lucasium stenodactylum, where genetically well-defined clades occur in close proximity at the southeastern edge of the Pilbara Craton [64].
The Pilbara has been proposed as an important Pleistocene refugium during the heightened arid conditions of the last 400,000 years [5,22]. During this time, the topographically complex region was isolated by the surrounding lowlands through the formation of vast sand dune deserts [65] that restricted the range and dispersal of species adapted to the older rocky, clay, or stony substrates, and facilitated the range expansion of species preferring sandy habitats [14,23]. The differentiation of subpopulations across the Pilbara boundary of S. youngsoni, S. macroura, and S. longicaudata was less distinct in the haplotype networks than S. ooldea, but it still produced significant Φ ST results, suggesting greater dispersal or habitat connectivity for these three species. While the Φ ST results suggest differentiation between the Pilbara and eastern populations of S. youngsoni, it is probably significantly weighted by the far eastern samples, as sandy habitats are continuous to the north and east of the Pilbara region ( Figure 1a).
Our data show that the southeastern Pilbara area may have contained a restricted population of S. ooldea in the past, as indicated in the mtDNA (but not in the nuclear DNA). The strongest evidence we found was for the Pilbara as an evolutionary refugium for widely distributed dasyurid species within the WA lineage of S. macroura. Haplotype diversity was higher within the Pilbara than outside of it, indicating the possible dispersal of individuals outside of the Pilbara that could be associated with a population expansion event.

Genetic Structure within the Pilbara
We hypothesized that species genetic structure within the Pilbara would be driven by dispersal ability, the distribution of preferred habitat types, and the location of potential biogeographic barriers. We recovered significant genetic structure between IBRA subregions of the Pilbara for S. macroura, S. youngsoni, and both Planigale species, with almost all pairwise Φ ST comparisons being significantly different. However, the mtDNA networks revealed that the underlying patterns among the study species were not congruent. Within the Pilbara, Planigale sp. 1 and S. macroura have comparable distributions and sample sizes but exhibited different phylogeographic patterns in the mtDNA network analysis. This difference probably reflects their different habitat preferences. Sminthopsis macroura prefers clay/loam soils in the Pilbara lowlands, and the rocky Hamersley Range could be an important linear barrier limiting dispersal from populations in the north, thus resulting in restricted gene flow into the Hamersley subregion. A similar pattern has been observed in species of rock-specialized geckos [22], a pebble-mimic dragon T. cephalus [24], and hummock (Triodia species) grasses [23], as genetic clades mirror the discontinuity of substrate types across the Chichester-Fortescue-Hamersley IBRA subregion boundaries. In contrast, the distribution of haplo groups 1 and 2 in Planigale sp. 1 occurred in a broadly east-west pattern through the central Pilbara ( Figure S2), and closely related haplotypes occur in the Chichester and the eastern Hamersley subregions. The clustering of clades in the eastern or western parts of the Chichester subregion was also observed in other gecko species and attributed to strong associations with habitat types and geology [22].
We did not find any clear phylogeographic structure in the mitochondrial network for Planigale sp. 2, but the Φ ST results between the Hamersley and Fortescue IBRA subregions were significant. This could result from low levels of gene flow across a barrier, such as the mountainous Hamersley Range, leading to haplotype frequency differences, but there is also enough gene flow still occurring to interrupt the development of clear phylogeographic structure. Phylogeographic structure was not detected within the Pilbara samples of S. longicaudata (Figure 3e). This surprising finding is despite this species being integrally associated with rocky substrates that have been shown to drive structure in other taxa, and maybe due to low sample size (n = 12).
The sand specialist in this study, S. youngsoni, occurs in dune systems on the coast, west of the Pilbara, and throughout the expansive longitudinal dunes of the Great Sandy Desert through to the Simpson Desert in South Australia. While these habitats are relatively continuous, we recovered mtDNA population differentiation on either side of the Gibson Desert, which is a desert dominated by hard, lateritic stony plains that are unlikely to be a suitable habitat for S. youngsoni and probably limit its dispersal to narrow corridors along waterlines. The narrowing of sandy habitats along the northwest Pilbara coast can be seen clearly in the distribution of records of S. youngsoni (Figure 1), but these samples still formed a haplo group with the eastern Pilbara samples. However, the Φ ST results between the Hamersley and Chichester/Fortescue subregions suggest some differentiation of populations in the east and western parts of the Pilbara region. Two other small dasyurids, Pseudantechinus macdonnellensis and Sminthopsis psammophila, which similar to S. youngsoni are distributed in sandy deserts, also did not exhibit any clear phylogeographic structure throughout their ranges [42,66]. The more recent formation of sand dune deserts has led to lower lineage age, intraspecific structure, and strong signals of population expansion in several gecko species [14,67] and hummock grasses [23] compared to related species occurring in more geographically structured habitats. However, we did not find significantly lower haplotype diversity in S. youngsoni, compared to other species in this study, suggesting that sandy habitats are just as capable of maintaining high levels of intraspecific diversity in dasyurids as other habitat types within and around the Pilbara.
We hypothesized that species occurring in restricted and patchy habitats would have stronger phylogeographic patterns than those inhabiting more extensive and continuous habitats [40]. However, this did not prove to be the case and is likely due to the high dispersal abilities of small arid-dwelling dasyurids [28]. Two species, S. longicaudata and Planigale sp. 2, that are strongly associated with isolated, patchy habitats, lacked strong population structure, despite our expectations. Both species are substrate specialists (Table 1), and they are infrequently captured on other substrate types [40]. The lack of population structure in our data could indicate that both these species have better dispersal abilities than previously documented and that population expansion has masked the signal of population structure. The lack of structure may also be a result of our data being limited in their power, having only three loci in total that includes two linked on the mitochondrial genome.
Our results add support to a growing number of studies that have found low population differentiation and potential panmixia within mammal species occurring in the Pilbara and other arid regions of Australia [66,68,69]. A recent study on three small mammals occurring in the Pilbara, two widespread habitat generalists (a murid, Pseudomys hermansburgensis and dasyurid, Ningaui timealeyi) and a pebble surface specialist mouse (Pseudomys chapmani) found little evidence for population structure using mitochondrial and microsatellite sequence data [68]. Levy et al. [68] proposed this resulted from the high fecundity and dispersal potential of these species, resulting in panmixia despite the complex geological landscape of the Pilbara presenting numerous potential barriers to gene flow. Panmixia has also been suggested for the mainland Pilbara population of the northern quoll [70], the largest Pilbara dasyurid (<1 kg), which is a species that has been documented to have excellent dispersal abilities [71].

Pleistocene Population Expansion
Byrne et al. [5] proposed that an intraspecific genetic structure may be the result of contraction and expansion from GM refugia, and that the timing of expansion should coincide with the last GM. We tested the hypothesis that population expansion among the study species likely would have occurred after the last GM. However, the timing of the onset of these expansion events we detected was much earlier than the last GM (approximately 15-25 kya), beginning instead around 100-200 kya (Figure 4). This expansion coincides with increasing aridity during peak climate oscillations of the mid-late Pleistocene and is comparable to butcherbirds [6], monsoonal/savannah occurring macropods [72], and clades of a widespread arid-adapted lineage of geckos [11], which experienced population expansion that also began during the mid-late Pleistocene. Our results are concordant with a growing number of studies that suggest population expansions of arid-adapted taxa occurred during the high amplitude climate cycles of the mid-late Pleistocene, coinciding with periods of intense aridification [14,15,23,67,73,74].

Conservation Implications for Arid-Occurring Dasyurids
Cryptic diversity is continuing to be a challenge in the conservation of Australia's species, including mammal groups [37]. Our tree analysis on the S. macroura complex supports the differentiated mitochondrial lineages identified by Blacket et al. [39], with an added signal from an independent nuclear locus. In particular, the strength of the S. froggatti + S. stalkeri group appears to provide evidence of separate species within the current S. macroura group, while preliminary examination of specimens in the WA Museum (L.S. Umbrello pers. obs.) has shown that specimens genotyped as S. froggatti can be readily distinguished from the western lineage of S. macroura. However, resolution of the relationships among the S. macroura clades was not well supported and follows a trend from previous studies of low support for the placement of these clades [39,75]. Higher relationships among the stripe-faced dunnarts require further investigation using multiple independent nuclear loci and comparative morphology. Until these relationships are resolved, we suggest that the western lineage of S. macroura is categorized as a distinct Evolutionary Significant Unit [76].

Conclusions
The arid zone comprises a multitude of complex, structured habitats that have facilitated rapid adaptive radiation and niche partitioning within Dasyuridae [77][78][79][80]. The species in this study were selected based on their broadly overlapping distributions and different substrate preferences [40], and we recovered some evidence of species-specific distributions of mtDNA haplo groups. However, we also found congruence in the broader genetic patterns of small dasyurids within the Pilbara, which was likely a result of shared evolutionary history and adaptations to aridity. Small dasyurids continue to persist in Australian arid habitats despite wide-scale extinctions and declines in other arid-adapted mammals [25], and our data add to a growing number of studies finding low levels of genetic structure in widely occurring arid species (reviewed in [81]). This is in contrast to species-specific responses to Pleistocene environmental change that are driven by differences in species ecological requirements [5,16] and demonstrates the adaptive success of small dasyurids to harsh and changing environmental conditions, which allow them to persist in these environments.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2073-4425/11/9/963/s1, Figure S1: Nuclear DNA networks, Figure S2: Concatenated mtDNA networks and maps showing the geographic distribution of major haplo groups, Table S1: Genbank Accession numbers and specimen voucher information for samples sequenced, Table S2: Primers used for PCR amplification and sequencing.