Genetic divergence between isolated populations of the North Island New Zealand Rifleman (Acanthisitta chloris granti) implicates ancient biogeographic impacts rather than recent habitat fragmentation

Abstract This research investigates the extent and causal mechanisms of genetic population divergence in a poorly flighted passerine, the North Island Rifleman or Titipounamu (Acanthisitta chloris granti). While this species has a historically widespread distribution, anthropogenic forest clearance has resulted in a highly fragmented current distribution. We conducted analyses of mitochondrial DNA (COI and Control Region) and 12 nuclear DNA microsatellites to test for population divergence and estimate times of divergence. diyabc and biogeobears were then used to assess likely past dispersal scenarios based on both mtDNA and nDNA. The results reveal several significantly divergent lineages across the North Island of New Zealand and indicate that some populations have been isolated for extensive periods of time (0.7–4.9 mya). Modeling indicated a dynamic history of population connectivity, with a drastic restriction in gene flow between three geographic regions, followed by a more recent re‐establishment of connectivity. Our analyses indicate the dynamic influence of key geological and climatological events on the distribution of genetic diversity in this species, including support for the genetic impact of old biogeographic boundaries such as the Taupo Line and Cockayne's Line, rather than recent anthropogenic habitat fragmentation. These findings present a rare example of an avian species with a genetic history more like that of flightless taxa and so provide new general insights into vicariant processes affecting populations of passerines with limited dispersal.


| INTRODUC TI ON
Intraspecific phylogeographic analyses provide vital information on the distribution of genetic diversity and levels of genetic connectivity across a species' distribution (Bermingham & Moritz, 1998;Frankham et al., 2002). Molecular indices of diversity may reveal hidden patterns of diversification among populations that are not necessarily evident using phenotypic measures (e.g., Burbridge et al., 2003;Daugherty et al., 1990) and may inform on patterns of historical gene flow. Such studies provide data regarding patterns of individual movement, population-level gene flow, and potential processes of speciation. The same data also provide essential information to conservation managers. Species that appear widespread may in fact be represented by a series of isolated sub-populations with varying demographic and diversity characteristics, and their genetic and phenotypic diversity may therefore be at higher risk of extinction than they first appear (Frankham et al., 2002;Koumoundouros et al., 2009). Small isolated populations with low genetic diversity are at greater risk of the effects of inbreeding, genetic drift, and potentially extinction, each of which are critical contributors to specieslevel evolutionary and management processes (Bakker et al., 2010;Frankham, 2005;Frankham et al., 1999Frankham et al., , 2002. The analysis of population genetic variation is particularly important for less dispersive species. In the presence of geographic barriers to gene flow, species with high dispersal potential may be able to surmount these barriers, whereas sedentary species or species with limited dispersal potential are at greater risk of population isolation, and therefore genetic fragmentation (e.g., Claramunt et al., 2012).
New Zealand's geological and human history makes it a particularly interesting place for phylogeographic studies of terrestrial species, especially those with low dispersal potential (Cooper & Cooper, 1995;Cooper & Millener, 1993;Goldberg et al., 2008;Trewick et al., 2017). This dynamic geological history includes extensive periods of land submersion, montane uplift, and volcanism, as well as repeated glacial cycles, all of which have significantly modified land connectivity and distributions of terrestrial species (Alloway et al., 2007;Fleming, 1962;Goldberg et al., 2008;McGlone, 2005;McGlone et al., 2001). Geological events dated back millions of years have had some of the most significant effects on both flora and fauna in New Zealand. For example, the lower North Island experienced a drastic reduction in land area during a period of extensive land submersion in the Pliocene, significantly restricting distributions of land species up until 2-3 mya (Bunce et al., 2009;Ellis et al., 2015;McGlone, 1985). Following this period of drastic range restriction, uplift along the tectonic plate boundary created an extensive geological barrier to dispersal for nonflighted land species (Cockayne, 1911;Ellis et al., 2015). A summary of the region's most important paleogeographic events is presented in Figure 1 and addressed in more detail in the Discussion. More recently, anthropogenic impacts have included extensive deforestation and the introduction of invasive mammalian predators, which have had significant detrimental impacts on both the abundance and distribution of native species (Gibbs, 2006;Holdaway, 1999;Veitch et al., 2011).
Multiple studies within New Zealand have demonstrated the significant impacts of ancient paleogeographic events on both biogeographic (among-species) patterns and on phylogeographic (within-species) patterns of nonavian terrestrial species (Ellis et al., 2015;Trewick et al., 2017). Molecular analyses have generally indicated that phylogeographic patterns are dominated by genetic division between North and South Islands, with relatively little substructure within the North Island itself (summarized in Trewick et al., 2017). Most population divergence within each of the two major islands appears to have occurred relatively recently for birds, due to either anthropogenic factors (e.g., Tracy & Jamieson, 2011) F I G U R E 1 A summary of the significant geological events impacting the distribution of land species in the North Island of New Zealand over three time periods, including the marine submersion of the lower North Island creating the Taupo Line (from McGlone, 1985;Wardle, 1963 (upper), and Rogers & McGlone, 1989 (lower)), and uplift of the axial ranges creating Cockayne's Line (Cockayne, 1911). Sampling locations for this study are also shown or other geologically recent, postglacial factors (e.g., Dussex et al., 2014;Miller & Lambert, 2006;Weston & Robertson, 2015).
The New Zealand Rifleman (Acanthisitta chloris) presents a unique focal species for research into the causal mechanisms of genetic divergence in passerines with lower dispersal potential. The Rifleman is a member of the endemic New Zealand wren family (Acanthisittidae), now considered as the sister group to all other passerines (Barker et al., 2004;Ericson et al., 2002). The family has contained some of the least flighted passerines in the world, including the extinct flightless Lyall's Wren (Traversia lyalli; Millener, 1989). The two remaining extant species within the acanthisittid wrens are the Rifleman (Acanthistta chloris) and the Rock Wren (Xenicus gilviventris). Both species have reduced dispersal ability due to their reduced tail and wing morphology and are therefore likely to have limited long-range dispersal potential (Higgins et al., 2001). The Rifleman is currently categorized into two insular sub-species, the North Island Rifleman (Acanthisitta chloris granti) and the South Island Rifleman (Acanthisitta chloris chloris). The North Island sub-species was once found throughout the North Island of New Zealand (Gill, 1996;Higgins et al., 2001), but is now limited to a highly fragmented distribution, being restricted predominantly to discontinuous highaltitude mountain ranges and wooded offshore islands (Robertson et al., 2007). Given the unique phylogenetic position of this species, studies on patterns of past dispersal are particularly informative to both biogeographic and phylogeographic investigations across the world, providing insight into species movements in a dynamic environment.
This investigation aimed to analyze the molecular diversity among geographically separated populations of the North Island sub-species of Rifleman. We addressed the following questions: (a) Is there evidence for significant genetic divergence among currently fragmented and isolated North Island populations? and (b) Is genetic isolation of populations likely to be related to past geological and climatological factors, or due to more recent anthropogenic forces of change? Given their low dispersal capability, and their previously widespread North Island distribution, we predicted that North Island Rifleman populations would show evidence of genetic population divergence following an isolation-by-distance pattern. We used a combination of mtDNA and nDNA to analyze population genetic divergence and to test several different models of past dispersal to determine which model best explains the current genetic distribution for Rifleman.  Figure 1). These areas represented Insular, Western, Central, Eastern Ranges, Eastern Coastal, and Southern locations, respectively. Rifleman were caught within their territories using 24 mm gauge mist-nets and conspecific lure calls between 2009 and 2012. Blood samples were collected using brachial venipuncture with a 29 gauge sterilized needle and a 20 µl capillary tube. Blood was stored in Queen's lysis buffer (Seutin et al., 1991)

| DNA extraction
Whole genomic DNA was extracted from blood samples using either standard phenol-chloroform extraction procedures or using a Qiagen DNEasy Blood and Tissue Kit. The Animal Tissue Modification protocol of the Qiagen Manual was used with the following modifications due to low DNA yields: (a) 160 µl of SET buffer replaced the recommended 200 µl of PBS buffer, (b) 40 µl of Proteinase K replaced the recommended 20 µl, (c) 30 µl of blood in Seutin buffer replaced 10 µl of whole blood, (d) incubation was carried out overnight as opposed to 10 min, and (e) final elution stages were done with two rounds of 100 µl of buffer AE as opposed to 200 µl.

| Mitochondrial DNA analyses
Mitochondrial DNA (mtDNA) polymerase chain reaction (PCR) targeted two regions of the mitochondrial genome for amplification.
COI PCRs were carried out in 25 µl reactions containing 12.83 µl of water, 2.5 µl of 10× reaction buffer, 1.25 mM of MgCl 2 , 2.5 µl of BSA, 0.2 mM of dNTP's, 0.17 µl of Taq Ti polymerase, 1.25 µl of each primer (1.25 µM) and 3 µl of template DNA (10-40 ng/µl). Thermal cycling conditions for amplification involved an initial partial denaturation phase at 94°C for 2 min and 35 cycles of denaturation at 94°C for 30 s, annealing at 57.5°C for 30 s, and extension at 72°C for 30 s, followed by a final extension step at 72°C for 4 min.
Primers targeting an approximately 600 bp region of the 5′ end of the mt Control Region (mtCR/CR) were designed specifically for Rifleman using a published mitochondrial genome from the South Island sub-species of Rifleman (Accession Number AY325307; Mitchell et al., 2016). The forward primer L16733 (5′-ACTTGGCACCTCCCCAAGACCA-3′) was located within the tRNA-Glu region, and the reverse primer H437 was situated within Domain II of the mtCR (5′-GGGTTGCTGATTTCTCGTGAG-3′). PCR reactions contained the same reagent concentrations as for the COI region as detailed above. Thermal cycling conditions were the same as for the COI region, but used an annealing temperature of 57.0°C. PCR products were visualized on a RedSafe-stained 1.6% agarose gel to check for amplification of single fragments of appropriate length before sequencing. Excess primers and nucleotides were removed from PCR products using the SapEx (Shrimp alkaline phosphatase and exonuclease 1) protocol (Werle et al., 1994), before carrying out cycle sequencing using the Big Dye protocol (Applied Biosystems).
Cleanseq (Agencourt) was used to purify products before sequencing on an ABI3130 automated sequencer, using the reverse primer COIBirdR2 for the COI region and the reverse primer H437 for the Control Region.
Sequences were viewed and aligned manually using Geneious Pro (version 5.5.6, Biomatters Inc.). All variable sites were confirmed by visual inspection of chromatograms. Resulting alignments were edited, and low-quality ends were trimmed to create a 377 bp alignment for the Control Region and a 652 bp alignment for COI.
Separate mitochondrial gene alignments were used for subsequent analyses, plus a concatenated dataset was created by combining Control Region and COI sequences. Variable sites were identified using MEGA (version 5.05) to enable manual haplotype designation.
For each separate alignment, the Akaike information criterion (AIC) was used in jModelTest (Guindon & Gascuel, 2003;Posada, 2008) to select the most appropriate nucleotide substitution model for the data (HKY + G). This model was then used to construct phylogenetic trees in Geneious Pro using neighbor-joining (NJ, 1,000 bootstraps), maximum likelihood (PHYML, 1,000 bootstraps, four substitution rate categories), and Bayesian (MrBayes, 10 7 chain length, subsample freq. 200, burn-in 10 6 ) methods for visual representation of haplotype relationships.
To visually represent relationships among haplotypes, a haplotype network was created for the COI gene, using a median-joining algorithm in Network v4.6.1.0 (Flexus Technology Ltd, 2011).

Population divergence based on the concatenated COI and Control
Region dataset was analyzed in Arlequin using AMOVA, along with pairwise calculations of F ST (using haplotype diversities only) and Ф ST (including nucleotide divergences). Statistical adjustment for multiple simultaneous tests was carried out using false discovery rate correction (Benjamini & Yekutieli, 2001). For the COI region, pairwise comparisons between populations were calculated using the Tamura distance method and a gamma level of 0.015, as the closest approximation to the HKY + G model selected by jModel test. Comparisons based on the Control Region divergence data were calculated using the Tamura distance method and a gamma level of 0.306. beAst 2.6.2 was used to estimate divergence dates by a coalescent method, based on the concatenated alignment, using reference sequences (Appendix S1 Table A1; Bouckaert et al., 2014).
PArtitionFinDer determined that HKY + G was the most appropriate substitution model for each partition of this dataset, based on a run of 100 million repeats (Guindon & Gascuel, 2003;Lanfear et al., 2012). The alignment was partitioned into four sections (COI positions 1, 2, and 3, and the D-loop) to allow for differing evolutionary rates. A lognormal relaxed clock, with a gamma site model, was used. Two sets of species divergence dates (Appendix S1 Table   A2) were used to calibrate the phylogeny, both based on Mitchell et al. (2016). The maximum calibration times, preferred by Mitchell et al. (2016), are based on a Cretaceous maximum constraint on the divergence of Passeriformes and Psittaciformes, a constraint on the age of the Acanthisitta lineage, and not using third codons in their analysis. The minimum calibration times were based on a Paleocene maximum constraint on the root of their passerine tree, no constraint on the age of the Acanthisitta lineage, and including third codons in their analysis (Mitchell et al., 2016). MRCA priors were set as normal distributions based on the divergence time means and confidence intervals determined by Mitchell et al. (2016). Each analysis was run with a chain length of 10 8 , and a burn-in of 10%.

| Microsatellite analyses
Twenty microsatellite primers designed for a kinship study in the South Island Rifleman  were trialed for PCR amplification, and the twelve best performing loci were used in this study.  identified and characterized thirty-seven polymorphic microsatellite loci for Kaikoura Tītipounamu. These microsatellites have been used to study kinship in the Kaikoura population (Preston, Briskie, et al., 2013). Primers for twenty of these loci were ordered; after initial trials, twelve loci were selected. These twelve loci amplified more consistently and could be more easily genotyped than the eight rejected loci. Details of the twelve loci are shown in the Appendix S1 (Table A3).
Polymerase chain reaction (PCR) was set up to a total volume of 10 μl, containing 6 μl of water, 1 μl of 10× PCR buffer, 0.5 μl of 25 mM magnesium chloride (MgCl 2 ), 1.5 μl of DNA template (10- This mix was heat-shocked at 95°C for 5 min and then cooled at 4°C. Genotyping was performed using an ABI3130xL Genetic Analyzer (Applied Biosystems). Individuals were genotyped based on band length using the microsatellite plugin in Geneious Pro Version 9 (Biomatters Inc.). Genotypes were binned to ensure a constant difference in allele length names of either four (for tetramer repeats) or two (for dimer repeats). STRUCTURE (version 2.3) was used to estimate the number of divergent population groups in the microsatellite data (Pritchard et al., 2000). A model that allowed for admixture and did not use sampling location as a prior was used. Each run used a 500,000 burn-in period and 1,000,000 chains for Markov chain Monte Carlo (MCMC). STRUCTURE HARVESTER was used to determine probabilistically what number of populations (k) best explains the data (Earl & Vonholdt, 2012) and GenALex (Peakall & Smouse, 2005). Evidence for the effects of a bottleneck on allelic variation within species was tested using the program bottLeneck Version 1.2.02 (Cornuet & Luikart, 1996).

Population divergence was analyzed using PoPGenrePort and
GenePoP. PoPGenrePort (Adamack & Gruber, 2014) was used to estimate pairwise F ST , Hedrick's F′ ST (Hedrick & Goodnight, 2005), and Jost's D (Jost, 2008), the last two of which take into account the downward bias in measuring population divergence caused by high heterozygosity within populations. structure (version 2.3) was used to estimate the number of divergent population groups in the microsatellite data, without location being used as a prior (Pritchard et al., 2000).

| Tests of models of past biogeographic history
Initial genetic results indicated that the patterns of mtDNA and nDNA (nuclear DNA) variation observed within the species revealed significant historical divergences among some populations, but that these genetic divergences may have begun to erode in more recent times. As such, we undertook a test of different biogeographic models of past dispersal, in order to determine whether a model that incorporated a recent change in the pattern of dispersal was the best at explaining the patterns of genetic diversity within the species. Three potential dispersal histories were modeled, based on the geographic pattern of genetic variation and the region's known paleogeographic history. The first null model (H0) simulated the scenario where there were no differences in dispersal restrictions among any locations, and there had been no changes over time. The first alternative hypothesis (H1) simulated the scenario suggested by the evidence of three distinct phylogeographic mtDNA lineages-that is, after initial colonization, there were strong dispersal restrictions between the three main regions (Insular, Mainland, and southeastern), and these were maintained until recent times. The second alternative hypothesis (H2) simulated the scenario suggested by the additional evidence of recent sharing of mtDNA and nDNA between Mainland and southeastern populations, but no recent sharing between the two southeastern locations-that is, after initial colonization, there were strong dispersal restrictions between the three main regions, and this was followed by recent breakdown of the barrier between Mainland and southeastern regions, and a new barrier to dispersal between the two southeastern locations (Southern and Eastern Coastal).
These three scenarios were modeled, and tested for fit to the genetic data, using two complementary methods. The first method used DiyAbc (Cornuet et al., 2010(Cornuet et al., , 2014 to compare the fit of the models to both the mtDNA and nDNA separately and then together, based on summary statistics. Each scenario was modeled through the specification of population divergence or admixture events (detailed in Figure 8). For each set of analyses, 10 7 datasets were simulated under each scenario, using standard models of microsatellite and mtDNA sequence. The set of summary statistics utilized were Tajima's D, mean sample pairwise differences and F ST for mtDNA, and mean genic diversity, F ST , and (dμ) 2 distance for microsatellites.
The posterior probabilities of each scenario were calculated and compared using the 500 closest simulated datasets in the direct approach, and confidence in the scenario choice was evaluated using 500 datasets from the posterior distribution (Cornuet et al., 2014).
The same three potential models of past dispersal were also compared using bioGeobeArs, which uses a phylogenetic approach to model geographic range distribution (Matzke, 2013). Model weight and parameters are provided in Appendix S1 Table A4. The dated mtDNA phylogeny was used as the basis of comparison, simplified to 14 major haplogroups, and the geographic locations of each haplogroup were recorded. The Rifleman evolutionary period was divided into three periods: 6-3, 3-0.2, and 0.2-0.0 mya, based on hypothesized periods of change in dispersal patterns. The three potential dispersal histories (represented in three time-stratified dispersal matrices, Appendix S1 Table A5) were compared: H0-no differences in dispersal restrictions among any locations (all dispersal parameters = 1), and no changes over time; H1-after initial colonization (all dispersal parameters = 1), strong dispersal restrictions between the three regions (insular, mainland and southeastern; dispersal = 0.01); H2-after initial colonization (all dispersal parameters = 1), strong dispersal restrictions between the three regions (dispersal = 0.01), followed by recent breakdown of the barrier between mainland and southeastern regions (dispersal = 1), and a new barrier to dispersal between the two southeastern locations (Southern/Tararua Ranges and Eastern Coastal/Mohi Bush, dispersal = 0.01). The three biogeographic models were compared across all possible base dispersal models: DEC, DEC + J, DIVALIKE, DIVALIKE + J, BAYAREALIKE, and BAYAREALIKE + J (Matzke, 2014). The most likely biogeographic model was chosen using AIC or AICc criteria, which both gave identical results across all possible dispersal models. bioGeobeArs is more usually employed in a biogeographic rather than a phylogeographic context, but in our application of this method, the OTUs represent monophyletic lineages that are analogous to the monophyletic populations more usually employed. Here, we were simply modeling the geographic range evolution of these lineages and were not assuming or testing specific biogeographic dispersal models at cladogenetic events.

| Mitochondrial DNA analyses
Mitochondrial DNA sequence data obtained from 90 Rifleman representing six populations across the North Island of New Zealand showed evidence of strong genetic structuring among geographic locations. In the 652 bp COI region, a total of 26 variable nucleotide sites (all involving substitutions) were observed that defined 17 different haplotypes among the six geographic regions (4% mean sequence divergence). The mtCR region contained 46 variable sites that were used to characterize 33 haplotypes in total (12% mean sequence divergence).
Haplotype and nucleotide diversity varied significantly among populations ( Table 1). The Insular population contained low haplotype diversity, with all individuals genetically identical for the COI region, but variable for the CR. While the highest haplotype diversity was found in the Western population, all haplotypes found there were closely related, resulting in relatively low nucleotide diversity.
The Southern population had relatively high haplotype diversity and the highest nucleotide diversity. Despite having low levels of haplotype diversity, the Eastern Coastal population was characterized by comparatively high nucleotide diversity.
Few haplotypes were shared among locations for both the COI and Control Region. The two mitochondrial regions displayed slightly different geographic patterns and were thus analyzed both separately and together in a concatenated dataset. Of the 17 haplotypes represented by COI sequences, only a single haplotype was shared between more than one population (Figure 2). The most common haplotype was found across the three central North Island populations and the Southern population and lies at a central node in the network, most likely representing an ancestral type (Crandall & Templeton, 1993). Most North Island mainland haplotypes diverged from this ancestral type by only a single substitution (Figure 2, inset).
A highly differentiated clade of haplotypes was associated with the Eastern Coastal and Southern populations. Although these two populations are dominated by haplotypes that are highly divergent from other mainland haplotypes, these populations also possess minor haplotypes that are closely related to the major ancestral clade. The Insular population was characterized by a single haplotype which was highly divergent, but intermediary, from all other COI haplotypes.
The mtCR exhibited high diversity with only one haplotype shared between any two populations (Central and Eastern Range populations). The phylogeographic patterns seen in the mtCR dataset were somewhat different from those seen in COI, but were highly similar to those displayed in the concatenated dataset, and so only the concatenated data are shown here and are best visualized in a phylogenetic tree (Figure 3)

| Microsatellite analyses
No microsatellite loci were found to depart from neutral expectations for tests of selection or departure from Hardy-Weinberg equilibrium. The population with the lowest observed heterozygosity (0.586) and allelic richness (1.55) was the Insular population, which also has low mtDNA diversity. As with the mitochondrial diversity in-

| Tests of models of past biogeographic history
The combined genetic results above revealed significant historical  Figure 8). This was supported by both the mtDNA and nDNA, but most strongly by the nDNA, and there was a low predicted error in these probabilities.
Comparison of the same three models of dispersal over time, using the mtDNA phylogenetic approach of bioGeobeArs, also showed that model H1 did not confer greater likelihood on the data over the null model H0, and that the model that conferred much higher likelihood across all biogeographic models was H2 (Figure 8).
Model H2 had the highest likelihood for all dispersal models (using all six basic dispersal assumptions) and together garnered 86.6% of the model weight across all variants (Figure 8).
Together, these results mean that the timing and direction of likely dispersal events correspond far better with the hypothesized history of model H2, than with H0 or H1, although a different model could exist that better fits the data. Our estimates of the divergence times among populations are likely to have some error associated with them. Firstly, they are derived from the mtDNA, which is effectively a single locus, and thus subject to considerable stochastic error (Edwards & Beerli, 2000).

| D ISCUSS I ON
Secondly, the estimated divergence times are for the mtDNA lineages, which are expected to be earlier than the divergences of the population themselves (Edwards & Beerli, 2000). Thirdly, the Bayesian method of deriving the mtDNA divergence dates from calibration points in the phylogeny may tend to over-estimate recent dates, particularly when the calibration points are from deeper in time. Thus, our estimates should be considered as maximum divergence times and interpreted in that context. At the same time, the estimates are based on a very confident phylogeny, using a broad range of well-established calibration points from the most closely related species, and that include the divergence period of interest ( Figure 4). As such, we believe our divergence time estimates are the most reliable currently available.
Given those provisos, the genetic separation of the Insular and all mainland sites appears to substantially predate human occupation of the New Zealand archipelago (approx. 1,000 years ago) and is consistent with the long geological isolation of this volcanic island (Hauturu-o-Toi/Little Barrier), which is estimated to have emerged approximately 3 mya with subsequent eruptions at approximately 1.6 mya (Lindsay et al., 1999). Land-bridges to the mainland (approx. 40 km to the south) are likely to have formed at least once during the Pleistocene, when glacial maxima would have caused significant lowering of sea levels Hamilton & Atkinson, 1961;Hare et al., 2008;Lindsay et al., 1999;Turbott, 1961 TA B L E 3 Mean node age estimates (mya) and 95% Highest Posterior Densities inferred using BEAST for maximum and minimum sets of parameters F I G U R E 4 Estimated times of divergence among Rifleman mtDNA lineages and across additional Acanthisittid lineages. Divergence times are estimated using beAst analysis and assuming species divergence times based on preferred Maximum Calibration times from Mitchell et al. (2016). Individuals are collapsed within each clade. Divergence times (in mya) are shown above bars representing their 95% HPD (highest posterior density). Posterior probability of node support is shown in italics on branches. Tree lineages are colored by evolutionary rate, ranging from blue (slowest) to red (fastest). *Divergence times at these nodes were set as priors based on Mitchell et al. (2016). **Extinct species between Insular and nearby mainland sites could not be investi-  North-South divide for some species of plants (McGlone, 1985;Wardle, 1963) and nonavian land animals (Buckley et al., 2010;Chapple et al., 2008;Marske et al., 2011). Various explanations for the existence of this biogeographic line include sea strait flooding (McGlone, 1985;Rogers & McGlone, 1989;Wardle, 1963), tectonic uplift (McGlone, 1985, volcanism Lloyd, 2003), and Pleistocene glacial cycles (Buckley et al., 2010;Lloyd, 2003). An alternative biogeographic line, Cockayne's Line, describes a barrier effectively splitting the southern North Island into western and eastern regions (Cockayne, 1911;Ellis et al., 2015). Patterns of endemicity, variation, and population divergence aligning with Cockayne's Line have been found in a range of plants and animals and appear to be a direct result of subdivision across the axial ranges, emergent during the Pleistocene (Baker et al., 1995;Buckley et al., 2010;Bunce et al., 2009;Holzapfel et al., 2002;Marshall et al., 2012;Nielson et al., 2011, reviewed in Ellis et al., 2015; Figure 1). Our data on Rifleman are unique in that they appear to show a dynamic pattern of population connectivity among populations of this sub-species, and we suggest that the distribution of genetic diversity within and among North Island Rifleman populations contains the genetic signal of both the Taupo Line and Cockayne's Line (Figure 8). F I G U R E 7 structure results for k = 6 groupings F I G U R E 8 Likelihood of potential changes in dispersal pattern in Rifleman over time modeled using BioGeoBears and DIYABC. The three models of dispersal compared were: H0-null; H1-after colonization, only barriers between Insular, Mainland and southeastern regions; H2-including recent breakdown of barrier between Mainland and southeast. Model H2 was most likely. *A, B, C refer to illustrated patterns. Likelihoods derived from bioGeobeArs analysis of mtDNA, posterior probabilities derived from DiyAbc analysis of mtDNA and nDNA One of the predominant explanations for the Taupo Line includes the impact of sea-level rise creating an extensive sea strait (the Manawatu Strait) during the Pliocene. This strait effectively restricted dry land to a series of islands in the upper region of the North Island and an emergent range in the southernmost part of the current North Island from approximately 23 mya up until as recent as 1.5-2 mya (Bunce et al., 2009;Ellis et al., 2015;McGlone, 1985). Support for the resulting biogeographic effect by the Taupo Line has rarely been found in vertebrates as it may require evidence of divergence prior to 5 mya (Ellis et al., 2015) which may be absent or obscured by more recent genetic patterns (Trewick & Bland, 2011). Rifleman appear to provide an unusual case in that our data demonstrate possible support for the influence of the Manawatu Strait on patterns of genetic connectivity in terrestrial vertebrates in the North Island, with the existence of the divergent southeastern mtDNA clade in the Eastern Coastal (Mohi Bush) and Southern (Tararua Ranges) populations ( Figure 8).

As the Manawatu Strait subsided and uplift in the southern
North Island occurred, the axial ranges were exposed, connecting the southern North Island and southeastern coastal "islands" approximately 1 mya (Bunce et al., 2009;Ellis et al., 2015;Marske et al., 2009;Trewick & Bland, 2011). This would have allowed Rifleman populations resident in the southern refuge of the North Island (including the Tararua Ranges) to colonize these eastern coastal regions (including Mohi Bush), while areas to the north-west (including our Mainland populations) may have continued to be isolated due to the emerging ranges, along Cockayne's Line. Our analyses are consistent with this hypothesis, with a much shallower mtDNA divergence time found between the Southern (Tararua Ranges) and Eastern Coastal (Mohi Bush) populations (0.7 to 1.5 mya) than with other mainland populations to the north and west.
In turn, our data also suggest that dispersal patterns have changed substantially through time. Although the three major mtDNA clades indicate very old restrictions to female dispersal between the Insular, Mainland, and southeastern regions, it appears that in more recent times, intermittent dispersal has occurred between central mainland populations and those to the south and southeast. This is evidenced by the small proportion of Mainland clade mtDNA haplotypes that are now found in the Southern (Tararua Ranges) and to a lesser extent in the Eastern Coastal (Mohi Bush) populations and by the microsatellite data that suggest extensive recent dispersal between the Mainland region and the Southern population. Dispersal modeling using both DiyAbc and bioGeobeArs supports the hypothesis that dispersal patterns have changed in the recent past, with increased gene flow from Mainland to both Southern (Tararua Ranges) and Eastern Coastal (Mohi Bush) populations, accompanied by reduced dispersal between the Southern and Eastern Coastal populations. The nDNA data also suggest that nuclear gene flow has now eroded much of the divergence of the Southern (Tararua Ranges) population, but little of the Eastern Coastal (Mohi Bush) divergence. One reason for this difference between mtDNA and nDNA divergence is likely to be the greater dispersal of males, which will more rapidly break down nDNA divergence among populations. This evidence for changing patterns of dispersal over time concurs with known paleogeographic patterns and may also demonstrate the dynamic nature of forest distribution throughout the late Pleistocene in response to glacial cycles and geological impacts McGlone, 1985;Trewick & Bland, 2011).
Populations in the central and north-west of the mainland (the Mainland clade) were much more closely related to each other than to the southeastern populations. The shallow phylogenetic differentiation of sequences and the star-shaped COI haplotype network evident within the Mainland clade are both evidence of more recent population expansion (Dixon, 2011;Rogers & Harpending, 1992 The results of this study differ markedly from previous studies on other New Zealand bird species, where fragmented populations within each main island appear to be more recently isolated, due to either recent anthropogenic factors (e.g., Baillie et al., 2014;Tracy & Jamieson, 2011), or postglacial factors (Baillie et al., 2014;Dussex et al., 2014;Miller & Lambert, 2006;Weston & Robertson, 2015).
The deep genetic divisions shown in this study are even remarkable when compared to extinct, flightless New Zealand species such as North Island moa (Dinornithiformes) populations, which diverged only recently during Pleistocene glacial cycles (Baker et al., 1995).
Worldwide, species divergence patterns in similarly less dispersive groups provide us with insights into how population divergence such as that seen in the rifleman could lead to speciation. In lowland neotropical antpittas (Grallariidae), ongoing species divergence extending back into the Oligocene is likely to have resulted simply from slow and constant diversification over long time scales, with low dispersal acting to maintain boundaries between variants (Carneiro et al., 2018). Similarly, in less dispersive Old World Eurylaimides, ancient vicariance events combined with more recent climatic and geological processes have resulted in species divergence, with historic patterns of gene flow evident due to the lack of dispersal (Moyle et al., 2006). Our data on the Rifleman appear to demonstrate a similar effect in that, due to restrictions in dispersal ability, a lack of extensive gene flow has allowed us to view the long-term patterns of ancient diversification within this species. Indeed, in a New Zealand context, Rifleman appear to exhibit biogeographic patterns more synonymous with several studies on nonavian taxa as opposed to avian species. Support for the Taupo Line has been found in beetles (Marske et al., 2011), while some native skink species contain deep population divergence hypothesized to result from Pliocene volcanic activity . In turn, support for Cockayne's Line, and therefore the impact of Pleistocene mountain uplift, has been found in stick insects (Clitarchus sp.; Buckley et al., 2010), geckos (Diplodactylidae; Nielson et al., 2011), cicadas (Kikihia sp.; Ellis et al., 2015), and in the plant genus Dactylanthus (Holzapfel et al., 2002). Rifleman therefore appear to be a remarkable case among birds in that they show the genetic signal of past dispersal barriers on a scale usually seen only in invertebrates and plants.
This study provides a unique perspective on the impact of ancient vicariant processes on the genetic distribution of avian species with reduced dispersal. These results demonstrate how the Rifleman, a species with high antiquity, no near relatives, and low dispersal ability, is an exceptional candidate to provide a glimpse into ancient historical vicariant processes affecting land birds.

ACK N OWLED G M ENTS
The authors would like to thank the following people for assistance The University of Auckland Ethics Committee (R762).

CO N FLI C T O F I NTE R E S T
There are no conflicts of interest for any author listed in this manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
mtDNA and nDNA data are located in GenBank with accession numbers (MW626949-MW627125).