Historic hybridization and persistence of a novel mito-nuclear combination in red-backed voles (genus Myodes)

Background The role of hybridization in generating diversity in animals is an active area of discovery and debate. We assess hybridization across a contact zone of northern (Myodes rutilus) and southern (M. gapperi) red-backed voles using variation in skeletal features and both mitochondrial and nuclear loci. This transect extends approximately 550 km along the North Pacific Coast of North America and encompasses 26 populations (n = 485). We establish the history, geographic extent and directionality of hybridization, determine whether hybridization is ongoing, and assess the evolutionary stability of novel genomic combinations. Results Identification of M. rutilus and M. gapperi based on the degree of closure of the post-palatal bridge was concordant with the distribution of diagnostic nuclear MYH6 alleles; however, an 80 km zone of introgressed populations was identified. The introgressant form is characterized by having mitochondrial haplotypes closely related to the northern M. rutilus on a nuclear background and morphological characteristics of southern M. gapperi. Conclusion Introgression appears to have been historic as pure populations of M. rutilus are now isolated to the north from introgressants or pure M. gapperi by the LeConte Glacier. As we do not find pure M. rutilus or M. gapperi individuals throughout the distribution of the introgressant form, it appears that the introgressants are a self-sustaining entity not requiring continued hybridization between pure parental forms to generate this novel combination of characters.


Background
The evolutionary significance of hybridization has been widely recognized in some taxa such as plants, but our understanding of how this process contributes to animal diversity, especially in vertebrates is relatively limited [1][2][3]. Ecological and demographic settings known to contribute to hybridization between otherwise well-defined species include newly established contact (such that behavioral, pre-zygotic filters to breeding may not exist) and low density of one or both parental species such that conspecific mating opportunities are limited [4]. Such biogeographic and demographic conditions are known to characterize many areas of interspecific, post-glacial contact where the leading edges of the expanding ranges of post-glacial colonizers meet [5][6][7][8]. Such regions provide a unique opportunity to examine the origin and maintenance of hybrid forms within an increasingly well-understood biogeographic and temporal framework.
Hybrid zones between species are thought to be maintained by two primary classes of models that predict 1) differential fitness between pure parental and hybrid individuals and 2) differential degrees of spatial overlap between pure parental species. The 'tension zone' model [4] posits a dynamic balance between selection against hybrids and dispersal of hybrids into the zone while the 'bounded superiority' model [9] holds that hybrids are superior to pure parental types in a limited set of environments. As such, the tension zone model necessitates continued generation of hybrid offspring from pure parentals so such systems would be characterized by sympatry of pure parental and hybrid forms on a spatial scale that encompasses the dispersal distances of the focal species. Alternatively, if hybrids have a selective advantage, even in a limited set of environmental conditions, we would predict little to no spatial overlap between pure parental and hybrid individuals. If the distribution of hybrid-appropriate habitat is at a scale that greatly exceeds the dispersal distance of pure parentals, maintaining occupancy of these areas would require that the introgressed form become independently sustaining. It is through this latter scenario that introgressive hybridization could establish stable, evolutionarily independent populations.
A particularly common hybrid form is a pattern of mitochondrial introgression leading to novel cytonuclear combinations. Because mitochondrial DNA is maternally inherited, this pattern can exist in animals where females are the homogametic sex (i.e. mammals). Due to reproductive inferiority of the heterogametic sex (Haldane's rule), the perpetuation of such novel cytonuclear combinations is prevented in animals with heterogametic females [10]. A plausible model for a pattern of mitochondrial introgression would be one wherein interspecific hybrids initially backcross with the most available pure parental species and then become a self-perpetuating entity, perhaps even displacing pure parentals. Under such a model, introgressant hybrids would be characterized by a mitochondrial genome of one parental taxon on a genomic background composed predominantly of that of a second taxon. Phylogenetic relationships and overall levels of diversity of mitochondrial types in the new introgressant hybrid form relative to the parental taxa could provide a great deal of insight into the origin of the introgressant form as well as subsequent factors contributing to the maintenance of these novel combinations.
Where the distributions of the northern red-backed vole (Myodes rutilus) and southern red-backed vole (Myodes gapperi) meet, we have an opportunity to examine interspecific interactions at the leading edge of two expanding ranges. As the ice sheets retreated, M. rutilus expanded its distribution south from Beringia, while M. gapperi colonized northward into Canada approximately 13,500 years before present (Ka) [11,12]. Hybridization between the two species has been hypothesized based on discordance between morphological and mitochondrial traits in a limited number of field-collected specimens [13], convergence in allozyme variation in areas of contact [14] as well as successful interspecific crosses in the laboratory [15]. Despite their capacity to interbreed, M. rutilus and M. gapperi are not sister lineages and are separated by approximately 9% sequence divergence in the mitochondrial cytochrome b gene [16]. Our goals were to 1) establish the geographic extent and directionality of any introgression between these taxa in southeast Alaska, 2) assess whether hybridization between M. rutilus and M. gapperi is ongoing and/or historic, and 3) assess the evolutionary stability of any novel genomic combinations.

Post-palatal bridge morphology
There was a distinct break in character state of the postpalatal bridge near the Stikine River (Table 1). Individuals from Jap Creek (locality 7, Figure 1) and north had incomplete post-palatal bridges, which is the same character state observed in the M. rutilus from the reference sample (interior Alaska), and from throughout its range (Runck in prep). Individuals from Mallard Slough (locality 8, Figure 1) and south had complete bridges, characteristic of the M. gapperi reference sample (Minnesota) and from specimens of M. gapperi throughout its range ( Figure  1; Runck in prep). There were, however, two individuals (localities 3 & 4), north of Mallard Slough that possessed complete post-palatal bridges, and two individuals (localities 14 & 23) south of Mallard Slough that had incomplete post-palatal bridges. In interior Alaska, 5 of 46 M. rutilus individuals had complete post-palatal bridges suggesting a low incidence of natural variation within this character in the northern red-backed vole.

Nuclear gene MYH6 genetic variation
Phylogenetic analyses of sequences of MYH6 revealed two clades with an average uncorrected pairwise divergence of 1.94% ( Figure 2). Maximum likelihood, neighbor-joining, and Bayesian analyses produced similar topologies. Individuals from Jap Creek (locality 7) and northward had alleles that formed a clade with M. rutilus from western Alaska and Russia ( Figure 2). Individuals from Mallard Slough (locality 8) and southward had alleles that formed a clade with M. gapperi from Minnesota and British Columbia. Uncorrected sequence divergence between the outgroup taxa and the M. gapperi and M. rutilus clades were 2.98% and 1.80%, respectively. None of the sequenced individuals was heterozygous at the diagnostic nucleotide positions.
Further sampling of the species-specific SNP using RFLP analyses resulted in the digestion pattern diagnostic of M. rutilus in all individuals from Jap Creek (7) northward, while all individuals from Mallard Slough (8) southward had MYH6 alleles diagnostic of M. gapperi. We found no individuals heterozygous for species-specific alleles.

Mitochondrial genetic variation and character concordance
We found 46 unique mitochondrial cytochrome b haplotypes in our entire sample. Patterns of variation across haplotypes were characteristic of functional mitochondrial genes, with average base frequencies of A (29.3) G (13.5) C (29.6) T (27.6), a 3.64 transition/transversion ratio, and a gamma distribution of 0.75 of changes across classes of codon sites. In addition, the distribution of the variable amino acid residues was consistent with the model of variable and conserved regions in cytochrome b [17].
Phylogenetic reconstruction using maximum likelihood, Bayesian, and neighbor-joining methods produced simi-lar topologies consisting of two highly supported clades, A and B (bootstrap support of 100; Figure 3). Average pairwise sequence divergence between southeast Alaska individuals in clades A and B is 7.2% (uncorrected p distance). Clade A haplotypes extended from Siberia, through western Alaska, and along the southeast Alaskan coast north of Hut Point (15; Figure 3). Clade B haplotypes were found in Minnesota and British Columbia as well as the southeast Alaskan coast from Hut Point (15) southward to Gwent Cove (19) on the mainland, Tyee (11) to Bond Bay (22), on the Cleveland Peninsula, and Wrangell (23), Etolin (24), and Revillagigedo (25,26) islands. Through additional RFLP screening of cytochrome b, we determined that all individuals north of the Stikine River (localities 1-8) and individuals up to 80 km south of the Stikine River (localities 9 -15) had Clade A mitochondrial haplotypes. Clade B haplotypes were Individuals examined for morphological and molecular data. Locality numbers correspond to Figure 1. Numbers of individuals examined for cytochrome b and MYH6 include direct sequencing and RFLP analyses. Contact of M. gapperi and introgressants is found at localities 11, 12, and 15.
Individuals with Clade A haplotypes that were also characterized by the post-palatal bridge morphology of M. rutilus and MYH6 alleles concordant with this morphology were only found north of Jap Creek (7). Hereafter, we refer to these individuals as M. rutilus. Individuals with Clade B haplotypes and the post-palatal bridge morphology of M. gapperi and MYH6 alleles concordant with this morphology were found on Etolin (23), Wrangell (24), Revillagigedo (25,26) islands, the Cleveland Peninsula (20 -22), and south of Hut Point (15), (hereafter M. gapperi). All individuals from Mallard Slough (8) south to Berg Bay (10) and Unuk and Chickamin rivers (13 & 14) had the post palatal bridge morphology and MYH6 sequence characteristic of M. gapperi but a set of cytochrome b haplotypes more closely related to pure M. rutilus (hereafter introgressants; Figure 1). Introgressants were found in sympatry with M. gapperi only at localities 11, 12, and 15 ( Figure 1; Table 1). Though haplotypes found within the introgressant individuals as a group were most closely related to M. rutilus, none was identical to haplotypes found in populations of M. rutilus. A synapomorphic mutation (position 513) was shared by all but one individual (Reflection Lake; locality 12) of the 45 sequenced introgressants.
To further resolve haplotype relationships within the two major clades, we constructed statistical parsimony networks using cytochrome b sequence variation without violating the parsimony criterion, as haplotypes within each clade were ≤ 10 mutational steps away from each other. The Clade A network had two major subclades, one consisting of haplotypes only found in the introgressant form ( Figure 4; indicated in grey) and the other consisting of all pure M. rutilus individuals and one introgressant ( Figure  4). The maximum distance between pure M. rutilus individuals was five mutations (0.008%) and between introgressant individuals six mutations (0.01%). Two ancestral haplotypes were inferred in the Clade B network (internal nodes) and a maximum distance of 5 mutations (0.008%) was found among these M. gapperi haplotypes ( Figure 5).

Demographic history and molecular evolution
Values obtained through Fu's test of selective neutrality were largely negative and significantly different from zero (Table 2), which is expected for populations undergoing recent growth. However, negative values can also be a result of selection. The mismatch distributions were unimodal for M. rutilus, M. gapperi, and the introgressants, which is expected for populations undergoing sudden expansion or under certain selective regimes ( Figure 6).
Estimates of the timing of any expansion events (τ ) were very similar in M. gapperi and the introgressants, but more recent in M. rutilus (

Discussion
After the retreat of Pleistocene glaciers, distributions of many high latitude organisms shifted resulting in new species assemblages and opportunities for genetic and ecological interactions [5][6][7][8]18]. Historical genetic signatures of heterospecific mitochondrial genes may be preserved in hybrid zone populations that are no longer undergoing genetic exchange. Examination of these populations using population and coalescent methods should provide insight into the timing and dynamics of geographical shifts in species' ranges in response to climate change [8,19].

Colonization and hybridization dynamics
Estimates of expansion times into southeast Alaska obtained from the mismatch distribution (t = τ /2u) indicate that these species arrived post-glacially. Working under the assumption that rates of evolution are consistent among these three groups, estimates of time since expansion are similar in M. gapperi and the introgressants, dating to 9.46 Ka and 8.57 Ka, respectively. As introgressants largely reflect the genetic signature of the hybridizing M. rutilus, expansion of both species into southeast Alaska date back to the early Holocene [21,22]. These estimates of expansion into southeast Alaska are notably earlier than the expansion time of those pure populations of M. rutilus that now exist north of the LeConte Glacier (4.31 Ka).
Given our refined view of the distribution of M. rutilus and M. gapperi and potential barriers, there appears to be no contemporary contact between these species in this transect. Consistent with that view, is the lack of MYH6 heterozygotes with alleles diagnostic of pure M. rutilus and M. gapperi that might suggest ongoing hybridization with M. rutilus [23]. Likewise, inspection of the post-palatal bridge, presumably controlled by multiple nuclear loci, revealed no intermediate morphs.
In addition to not finding evidence for contemporary gene exchange between these two species, we did not find the parental species in sympatry. Furthermore, there was not extensive overlap of the introgressants with the pure parentals. M. gapperi and the introgressants occur in sympatry only at the southern edges of the zone of introgressants (localities 11,12,15), and the current tidewater position of the LeConte Glacier prevents contemporary contact of pure M. rutilus with voles in areas farther south that are now occupied by introgressants or, even farther south, by pure M. gapperi. As such, it would appear that the introgressants are self-sustaining populations and not hybrids that are continuously generated from pure parental crosses. An active contact zone does exist, however, at the southern, leading edge of the introgressant distribution and the northern edge of pure M. gapperi (at localities 11,12,15). The direction and degree of genetic exchange at the latter contact zone is the subject of an ongoing study (Runck et al., in prep). Future ecological studies should explore whether abiotic or biotic shifts are associated with the transition between these groups and/or whether direct competitive interactions limit their coexistence.

Origin of mitochondrial signature in introgressants
The introgressant form is characterized by a monophyletic group of mitochondrial haplotypes nested within haplotypes otherwise characteristic of M. rutilus. Nonetheless, the introgressed group possesses a distinct mitochondrial genetic signature from M. rutilus populations across this region. Notably, not only do the introgressants not share any haplotypes with M. rutilus, they also have significantly greater mitochondrial variability than the donor species.
Three plausible scenarios could lead to this distinct genetic signature of novel, but closely related, haplotypes and overall higher mtDNA diversity that characterize the introgressants. First, increased mitochondrial diversity in the introgressants may simply reflect differences in population history when compared to M. rutilus in southeast Alaska. Marginal support for the multiple waves of colonization of M. rutilus hypothesis is reflected in the relationship of the introgressant mtDNA relative to pure M. rutilus in the minimum spanning network and likelihood tree. In both analyses, introgressants form a subgroup, and are not intermixed with the remaining M. rutilus mitotypes suggesting that these two mitotypes were not part of a panmictic population. The genetic footprint of an earlier hybridization event supports the hypothesis of multiple colonizations of the northern red-backed vole along the coast.

Evolution of an introgressant contact zone
Although genetic exchange may have been extensive, our data suggest that a very closely related set of haplotypes (or a single haplotype that subsequently mutated) was captured and maintained in the introgressants. The unimodal distribution of pairwise comparisons and position of the introgressant haplotypes in the network and phylogeny support the hypothesis of a single hybridization event instead of multiple temporally discrete events. Initial genetic exchange may also have been bidirectional, but we only have evidence thus far of the mitochondrion of M. rutilus being maintained on the morphological and presumably nuclear background of M. gapperi. Nonetheless, a more complete view of the nuclear composition of the introgressant form relative to the pure parental forms will provide insight to the extent of backcrossing that occurred in this system.
Once established, the novel mito-nuclear combination of the introgressants either diffused neutrally or was selected for and expanded their distribution while displacing the pure parentals. While we cannot determine which scenario is responsible for the 80 km zone of introgressants, it is notable that all individuals in this region are introgressants. Moreover, several of our findings are consistent with predictions of the bounded superiority hypothesis [9]. The introgressants and pure parentals do not overlap extensively and the introgressants are self-sustaining and are not the result of continual hybridization. Also, consistent with the bounded superiority hypothesis is that the introgressants occupy a limited area of southeast Alaska, and have not established populations west of localities 11 and 12 on the Cleveland Peninsula.
The degree of genetic admixture that may have occurred while the introgressant voles were fairly uncommon relative to pure parentals awaits our more complete sampling of the nuclear genome of this group. Similar long-term persistence and spatial expansion of introgressant or hybrid forms has been documented in snails (genus Cerion), resulting from an ancient hybridization event between a now extinct fossil species and an extant species [31]. Hybrids are hypothesized to have persisted due to the novel genetic combinations that enhanced survival during the time that one of the parental species was eliminated [31].
With the advance of the LeConte Glacier to tidewater approximately 5,000 years ago, any potential for gene flow between northern M. rutilus and the introgressants ceased. However at the southern edge, we do find localities (11, 12, and 15) where introgressants overlap spatially with pure M. gapperi suggesting the potential for ongoing gene flow. Ongoing analyses of microsatellite genetic variation should identify not only the degree of overall distinction between the introgressed form and pure M. gapperi but also the amount of gene flow that characterizes their current contact zone. Likewise, ecological studies may help identify the degree to which these groups compete either directly or indirectly with one another and whether the spatial distribution of each is currently stable or actively shifting.

A North Pacific Coast suture zone
The North Pacific Coast has been documented as post-glacial contact zone for several mammalian lineages from the high latitude refugium called Beringia and from multiple refugia that existed south of the continental ice sheets [13]. Within species, independent colonizations into this recently deglaciated region have occurred by at least two divergent lineages of several species such as dusky shrew (Sorex monticolus) [32], long-tailed vole (Microtus longicaudus) [33], black bear (Ursus americanus) [34], and marten (Martes americana) [35,36]. Ermine (Mustela erminea) are represented by three divergent lineages in southeast Alaska; one is hypothesized to be endemic to the region, perhaps surviving in the North Pacific Coast during the Pleistocene [24].
Introgression along the coast has been documented from divergent lineages of marten [35], black bears [37], and now red-backed voles. Contact zones are likely for divergent clades of dusky shrews and long-tailed voles, therefore the narrow strip of mainland along the southeast Alaska coast may be a suture zone, whereby several formerly isolated species have entered the region by discrete colonization routes, and have subsequently come into contact in the same geographic area [38][39][40]. The North Pacific Coast, and in particular, southeast Alaska, possess characteristics [38] commonly tied to suture zones, such as being located between Pleistocene glacial refugia, and nearby low mountain passes acting as corridors for dispersal [13,41] during warm periods.
A renewed interest has emerged in testing the validity of Remington's thirteen North American suture zones through phylogeographic studies [39,40,42]. The North Pacific Coast was not originally identified as a suture zone by Remington [38], but phylogeographic studies repeatedly demonstrate the existence of multiple lineages within species (e.g., shrews, voles) in this region. Populations representing divergent lineages are now in contact there following postglacial expansion [43]. Most of these studies, however, have limited ability to detect hybridization because only mitochondrial genes were assessed. Future phylogeographic studies of these species should employ multiple independent characters to more rigorously assess the influence of geologic and climatic events on structuring diversity along the North Pacific Coast.

Conclusion
Interspecific hybridization between M. rutilus and M. gapperi resulted in the formation of an introgressant group that spans 80 km. The novel mito-nuclear combination of these introgressants likely expanded either by displacing or by colonizing areas left unoccupied by pure parentals in response to the changing climate of the Holocene. Hybridization between the two species is historical, as a region occupied exclusively by the introgressants now separates these two species. Additionally, physical separation sant populations occurred ca. 5,000 Ka thus establishing reproductive isolation of this pure parental species and the introgressants. These introgressive populations appear to be stable as they are not a result of continuous interspecific matings between the pure parental species.

Sampling
The temperate rain forest of coastal southeast Alaska is naturally fragmented by extensive icefields, fjords, and six major rivers, and is isolated from continental North

Morphological data collection
Closure of the post-palatal bridge is the key diagnostic character used to distinguish M. rutilus from M. gapperi ( Figure 7) [11,20,45,46]. All individuals from the transect were examined, but in some instances the character state could not be determined due to skull damage, resulting in a total 328 individuals analyzed. Character states were compared with M. rutilus skulls from interior Alaska (n = 46) and M. gapperi skulls from Minnesota (n = 18). Previous analyses of these specimens and others show that degree of post-palatal closure is not correlated with sex, age, or latitude, (Runck in prep). Specimens were examined under 25× magnification to score the post palatal bridge either as A) complete, with medial shelf connected to lateral parts of the palate, or B) incomplete, with medial shelf disconnected (Figure 7).

Molecular methods
Genomic DNA was extracted from frozen liver of 95 M. rutilus and 390 M. gapperi following a modified salt extraction method [24,47].

Nuclear locus, MYH6
Twenty-seven nuclear genes were screened from the collection of comparative anchor tagged sequences to obtain a diagnostic nuclear perspective on species identification and to test for interspecific gene flow in the contact zone [48]. Five loci yielded PCR products in the 300-1000 bp range. We used the primer set MYH2F and MYH2R to generate a 257 bp fragment of MYH6 (myosin heavy polypeptide 6). Although this primer set was originally designed to amplify MYH2 (myosin heavy polypeptide 2) [48], a nucleotide BLAST search against the entire nucleotide collection in GenBank indicated that the amplified fragments are MYH6.
The identification of species-specific alleles was conducted using broadly distributed individuals of M. rutilus from Russia, Finland, and western Alaska and individuals of M. gapperi from Washington, Minnesota, and North Carolina. Outgroup taxa M. rufocanus and M. glareolus were also sequenced [16]. From these samples we found three conserved single nucleotide polymorphisms (SNPs) that distinguish the alleles of M. rutilus and M. gapperi. We consider this locus to be diagnostic for the two species because these diagnostic nucleotide differences were found throughout the species' ranges.  [16,50]. Cytochrome b fragments were digested using 3.5 units of ALU I, 1.0 ul Buffer2 (New England BioLabs), and 5.0 ul PCR product. Positive and negative controls were included to confirm enzyme activity in each trial. PCR products were digested for 5 hours at 37°C and then visualized on a 2.0% agarose gel stained with ethidium bromide. Digestion of the cytochrome b fragment resulted in two bands for M. rutilus and one band for M. gapperi. Distributions of the two haplotypes were mapped.

Phylogenetic analysis
We analyzed 257 bp of the MYH6 locus using maximum likelihood and neighbor-joining algorithms in PAUP*b10 [52] and Bayesian statistics in Mr. Bayes v3.4 [53]. Individuals of M. rutilus from western Alaska and Russia and individuals of M. gapperi from Minnesota and British Columbia were also included in the analyses. M. rufocanus and M. glareolus sequences were included as outgroups. The model HKY [54] was determined as the simplest model that best fit these data through Modeltest [55] and was used in the analyses. Neighbor-joining analyses used 1000 bootstrap replicates to assess nodal support. The Bayesian analysis started with a random tree and was run for 1.0 × 10 7 generations sampling every 1,000 generations. A consensus of three runs was computed after stationarity was reached.
The first 600 base pairs of cytochrome b from 172 southeast coastal individuals were used in phylogenetic analyses. M. rutilus from Russia and western Alaska and M.  [57] model of nucleotide evolution with invariable sites (0.4036) and gamma distribution (0.7468) was identified using Modeltest as the model that best fit these data, and was used in subsequent analyses. Phylogeographic relationships were reconstructed using maximum likelihood and neighborjoining algorithms in PAUP*b10 [52]. Nodal strength was assessed using 1000 bootstrap replicates. Bayesian statistics were also used to reconstruct phylogeographic relationships in Mr. Bayes v3.4 [53]. Bayesian reconstruction started from a random tree and was run with four heated chains for 1.0 × 10 7 generations, sampling every 1,000 generations. Three runs were conducted and nodal support (posterior probability) was computed from the three runs after stationarity was reached.
Two cytochrome b statistical parsimony networks were constructed for the two clades of haplotypes [58]. Haplotypes were connected based on the absolute number of mutational differences and coalescent theory was used to identify ancestral (internal) and derived (tip) relationships [59] in the program TCS v1.21 [60]. Ten mutational steps between any two haplotypes is the maximum difference allowable in order to reconstruct the relationships while meeting the parsimony criterion with 95% probability [58].
The M. rutilus haplotype network was constructed using 11 haplotypes representing 47 M. rutilus individuals and also included 17 M. rutilus-like haplotypes found in 45 individuals identified morphologically as M. gapperi (introgressants). The M. gapperi haplotype network was constructed using 18 haplotypes representing 80 M. gapperi individuals.

Demographic history and molecular evolution
Variation of cytochrome b sequences was used to test for signals of population expansion, to estimate the time since expansion, and to estimate mutation rates. Haplotype diversity (h) and nucleotide diversity (π) were estimated [61] separately for M. rutilus, M. gapperi, and introgressants in the program Arlequin ver. 3.11 [62]. Differences in levels of diversity were tested using a t-test. The frequency distribution of pairwise differences (mismatch distribution) [63] and Fu's F S statistic [64,65] were calculated to test for population expansion in Arlequin ver. 3.11. Under the Sudden Expansion Model, pairwise differences will have a unimodal distribution [66]. Values of Fu's F S will be negative and indicative of expansion when there is an excess of singleton mutations or a when a gene is under selection [65].
The peak of the unimodal distribution (τ ) [66,67] in the mismatch distribution was used to calculate time since expansion in Arlequin. Applying the equation τ = 2 μt where μ is the mutation rate and t is time in generations [66], allows an estimate of time since expansion. Confidence intervals of τ were calculated using parametric bootstrapping [68]. Using the mutation rate of 7.5% per million years [69], and values of τ obtained from mismatch analyses, we estimated time since expansion for M. gapperi, M. rutilus and introgressants.
To test for constancy in rates of cytochrome b evolution among lineages, we conducted a relative rates test in rrTree v1.1 [70], which calculates synonymous and non-synonymous rates of evolution [71,72]. We also constrained the ML phylogeny to a molecular clock in PAUP*4.0b10 [52] and performed a log likelihood ratio test of log likelihood scores of constrained and unconstrained trees and compared these values to a χ 2 distribution.