Speciation is generally viewed as an irreversible process, although habitat alterations can erase reproductive barriers if divergence between ecologically differentiated species is recent. Reversed speciation might also occur if geographical contact is established between species that have evolved the same reproductive isolating barrier in parallel. Here, we demonstrate a loss of intrinsic reproductive isolation in a clade of scincid lizards as a result of parallel body size evolution, which has allowed for gene flow where large-bodied lineages are in secondary contact. An mtDNA phylogeny confirms the monophyly of the Plestiodon skiltonianus species complex, but rejects that of two size-differentiated ecomorphs. Mate compatibility experiments show that the high degree of body size divergence imposes a strong reproductive barrier between the two morphs; however, the strength of the barrier is greatly diminished between parallel-evolved forms. Since two large-bodied lineages are in geographical contact in the Sierra Nevada Mountains of California, we were also able to test for postzygotic isolation under natural conditions. Analyses of amplified fragment length polymorphisms show that extensive gene exchange is occurring across the contact zone, resulting in an overall pattern consistent with isolation by distance. These results provide evidence of reversed speciation between clades that diverged from a common ancestor more than 12 Myr ago.
A general pattern in animals is that reproductive compatibility decreases as the time since divergence increases, leading to the formation of new species (Coyne & Orr 2004). Exceptions to this pattern can occur if traits that determine reproductive compatibility evolve in phylogenetically independent lineages, causing descendant lineages with similar derived morphologies to be more reproductively compatible with each other than they are with their morphologically divergent sister lineages, a process called parallel speciation (Schluter & Nagel 1995). Since parallel speciation increases reproductive compatibility between some lineages, it is a natural process that has the potential to reverse speciation if the parallel-evolved forms come into geographical contact. Documented examples of parallel speciation involve lineages that originated relatively recently (tens to hundreds of thousands of years ago) and have been important in demonstrating a role for natural selection in speciation (Funk 1998; McPeek & Wellburn 1998; Rundle et al. 2000; Nosil et al. 2002; McKinnon et al. 2004; Rolán-Alvarez et al. 2004). However, spatial separation of the independently evolved forms precludes investigation of the evolutionary consequences of parallel speciation in most systems (but see Holloway et al. 2006).
Our study examined the basis of premating isolation in scincid lizards of the Plestiodon skiltonianus species complex, and tested for gene exchange between populations in secondary geographical contact that evolved large body size independently. The species complex is broadly distributed in western North America, extending from British Columbia to Baja California, and east into the Great Basin (Taylor 1935; Rodgers & Fitch 1947). Two ecomorphs are distinguishable, and these are classified into three species: a widespread, small-bodied form (P. skiltonianus and the geographically disjunct Plestiodon lagunensis) that occurs in relatively mesic habitats along the Pacific coast and at higher elevations, and a large-bodied form (Plestiodon gilberti) that occurs in more xeric habitats found inland and at lower elevations (Rodgers & Fitch 1947; Tanner 1957; Morrison et al. 1999; figure 1). Phylogenetic analyses of mtDNA sequences confirm the monophyly of the species complex, but reject the monophyly of each ecomorph (Richmond & Reeder 2002). Instead, P. gilberti consists of three clades that are nested among lineages of the small-bodied P. skiltonianus. The two forms co-occur in parts of southern California, where they behave as separate species. Sympatric populations maintain their distinctiveness in colour pattern and body size: juvenile skiltonianus have iridescent blue tails whereas gilberti have reddish-pink tails, the striped dorsal colour pattern is retained in adult skiltonianus but lost in gilberti, and average adult body size in skiltonianus reaches approximately 60% of that in gilberti (Rodgers & Fitch 1947; Stebbins 2003; J. Q. Richmond 2004, unpublished data). In addition to a lack of morphological intermediates, no mitochondrial introgression has been detected between the two forms (Richmond & Reeder 2002; this study).
Phylogenetic reconstruction of snout–vent length (SVL), a common proxy for overall size, shows that body size has increased independently in each of the three gilberti lineages (Richmond & Reeder 2002; Richmond 2006). Two of these lineages (the Sierran and southwestern clades; figure 1) approach each other geographically and mtDNA haplotypes belonging to the two lineages have been recovered within a distance of 11.2 km of one another near the juncture of the Sierra Nevada and Tehachapi Mountains in southern California, where the skinks are continuously distributed. This contact zone also coincides with a shift in juvenile tail colour from pink to bright blue (figure 1), and putative intergrades have been collected in the same area.
Parallelism, occurrence of possible intergrades and a common tendency for organisms to mate assortatively by size led us to ask whether separate lineages with the same ecomorphology are reproductively compatible, despite their independent origins. To this end, we used no-choice tests of mate compatibility to estimate the effects of three factors on reproductive isolation: body size divergence; degree of phylogenetic relatedness; and the geographical distance separating parent populations. All of these are important predictors of reproductive isolation in other organisms (Coyne & Orr 1989, 2004; Tilley et al. 1990; Endler & Houde 1995; Strecker & Kodric-Brown 1999; Schluter 2000, 2001; de Kort & ten Cate 2001; Nosil et al. 2002; Wong et al. 2004). Support for parallel speciation is gained if mate compatibility depends on body size divergence. This result would implicate natural selection as a driving force in the evolution of reproductive isolation in this group, as selection is the best explanation of parallelism in ecologically similar populations (Endler 1986; Schluter & Nagel 1995; Rundle et al. 2000).
Secondary contact between the Sierran and southwestern clades of P. gilberti provides a natural experiment for determining whether clades of large-bodied skinks are isolated by postzygotic factors, as might be expected given the length of time they have been diverging (Coyne & Orr 2004). To test this hypothesis, we used amplified fragment length polymorphisms (AFLPs) to assess the spatial genetic structuring of multiple independent markers in populations that span the mtDNA contact zone. If the two clades represent true biological species, gene flow should be low or non-existent across the contact zone (i.e. high Fst values), reflecting their ancestral separation. Alternatively, the parallel speciation model predicts that the loss of prezygotic reproductive isolation has led to the reestablishment of gene flow and subsequent merger of the lineages.
Our findings show that the probability of successful copulation decreases as body size divergence increases, indicating that ecomorphological differences cause reproductive isolation while simultaneously restoring mate compatibility between lineages that have evolved the same ecomorphology in parallel. Spatial structuring of AFLP data shows that extensive gene flow occurs between two gilberti lineages where their ranges adjoin in the Sierra Nevada Mountains of California. As a consequence, non-sister lineages in secondary contact now behave as a single biological species in the wild, despite their separation from a common ancestor more than 12 Myr ago.
2. Material and methods
(a) Mate compatibility analyses
No-choice mate compatibility trials were conducted using a total of 126 individuals (n=51 males, n=75 females; n=88 gilberti, n=38 skiltonianus). In each trial, a single male and female were placed in a testing chamber at approximately 28°C and observed for 45 min or until copulation was complete. Individuals were divided into 20 trial series, and each individual was used in only one trial series. Each trial series used a maximum of five individuals of each sex, chosen from five different mtDNA clades. Mate compatibility was assessed for all pairwise combinations of opposite-sex individuals within a trial series (figure 1 in electronic supplementary material), resulting in a maximum of five trials per individual and 25 trials per trial series, all conducted within a 5 day period. The order of trials within trial series was randomized subject to the constraint that each individual be used only once per day.
To induce simultaneous sexual receptivity and stimulate reproductive behaviour, we used hormone injections in many individuals. Such treatments have been shown to cause sexual behaviour that is consistent with normal breeding behaviour in the scincid lizard Plestiodon laticeps (Cooper & Vitt 1993) and the phrynosomatid lizard Gambelia wislizenii (Fallahpour 2004). Further details of these procedures are outlined in the electronic supplementary material. All trials were conducted at the University of Connecticut using skinks that were hand captured in California in 2003–2004 (California State Permit ID# SC-002294), and trials were conducted from March to August of both years. All work was approved by the Office of Research Compliance at the University of Connecticut (IACUC protocol no. A04-072).
The index of pair sexual isolation IPSI, which ranges from −1 to 1, was used to estimate the degree of sexual isolation between ecomorphs (Rolán-Alvarez & Caballero 2000). Negative values indicate disassortative mating, 0 indicates random mating and 1 indicates complete isolation. Significance was assessed using 10 000 bootstrap replicates in the software package JMating (Carvajal-Rodriguez & Rolán-Alvarez 2006).
To further explore the biological basis for sexual isolation between ecomorphs, we used a Bayesian approach to logistic regression and model comparison as implemented in WinBUGS v. 1.4.1 (Spiegelhalter et al. 2003). The effects of three factors on the probability of copulation were evaluated: body size difference; geographical distance separating parent populations; and evolutionary relatedness. We also tested whether unmeasured differences between individuals (random effects of individuals) contributed to the variation in copulation success and whether patterns were the same across ecomorph pairings. To control for the experimental design, we evaluated the effects of hormone use, trial series and day within trial series (to control for fatigue).
The contribution to copulation success of both biological and experimental factors was evaluated by comparing the fit of a series of models including and excluding each parameter of interest. Model fit was assessed with the deviance information criterion (DIC, a penalized fit measure) and comparisons followed a standard convention: if the ΔDIC is less than 2.0, models have equal support; from 4 to 7, the model with the higher score has less support than the competing model; if greater than 10, the higher-scoring model has no support (Burnham & Anderson 1998). Those variables whose inclusion improved model fit were considered to be significant predictors of copulation success. We further tested the relative effect of each explanatory variable by examining whether 0 was included in the 95% credibility interval (CI) of the slope coefficient estimated from the posterior distribution. Exclusion of 0 from the 95% CI indicates that there is a less than 5% chance that the true value of the coefficient is 0, and thus supports the conclusion that the variable is a predictor of copulation success. The Markov Chain Monte Carlo (MCMC) procedures used to approximate posterior distributions are described in the electronic supplementary material.
Data for the three predictor variables were calculated as follows. For body size, we used the absolute value of the difference in body size between paired individuals (male SVL−female SVL) minus the quantity Ω. Ω is the size difference that maximizes the probability of copulation and was estimated in the regression analysis. The geographical distances were estimated as the great circle distance separating collecting sites based on coordinates determined using a global positioning system. Phylogenetic distances were calculated from the posterior distribution of branch lengths of ultrametric trees of mtDNA haplotypes (n=196 unique haplotypes) inferred using BEAST v. 1.4.1 (Drummond & Rambaut 2003; Drummond et al. 2006). Details are provided in the electronic supplementary material.
We also determined the degree of prezygotic reproductive isolation expected in the wild by measuring the SVL of 935 reproductively mature museum specimens from throughout the range of both ecomorphs (figure 2b) and estimating the posterior distribution of δ, the difference between body size divergence and the size difference that maximized copulation success (Ω), within and between clades. Size divergences between males and females of each pair of focal taxa (n male skiltonianus=188, n male gilberti=257; n female skiltonianus=211, n female gilberti=279) were estimated in WinBUGS v. 1.4.1. The posterior distribution of δ was calculated as the difference between 10 000 samples from the posterior distributions of the difference between male and female SVL and 10 000 samples of the posterior distribution of Ω.
(b) Divergence date estimates
As we were interested in the evolutionary timeframe over which mate compatibility can be restored, a Bayesian relaxed clock method was used to infer the time to common ancestry of the three P. gilberti lineages. Neither fossils nor geologic events could be used to calibrate the rate of molecular evolution, so divergence dates were estimated using two substitution rates derived from cytochrome b sequences of lizard species with body sizes similar to those of the skinks (Tarentola delalandii: 8.25×10−3 substitutions per site per million years, Gubitz et al. 2000; Anolis occulatus: 7.15×10−3 substitutions per site per million years, Malhotra & Thorpe 2000). Nucleotide substitution models, coalescence models, priors and MCMC conditions for all BEAST analyses are presented in the electronic supplementary material.
(c) Estimates of gene flow
AFLPs were generated using standard protocols (Vos et al. 1995) with slight modifications (see the electronic supplementary material). Markers were size separated on an ABI3100 genetic analyser and Genotyper v. 3.6 NT (Applied Biosystems) software was used for locus identification and scoring. An initial subset of individuals from geographically disparate localities was used to determine a set of target peaks (i.e. categories; size tolerance ±0.5 bp) that were distinct and reproducible between runs. Peak presence or absence in each category was manually verified by eye, and ambiguous peaks (small peaks in the same-sized bin where other samples had large peaks) were coded as (?).
We divided individuals into eight geographical clusters and used a Bayesian approach to estimate θB (an analogue of Wright's FST, Wright 1951) between each neighbouring pair of clusters, and between pairs separated by increasing distance across the mtDNA contact zone. The posterior distribution for θB was estimated using vague normal priors for all parameters and the ‘f free’ model, which is appropriate for dominant markers, in the software package Hickory v. 1.0.3 (Holsinger et al. 2002; Holsinger & Lewis 2004; Holsinger & Wallace 2004). To further assess the geographical structuring of AFLPs, Mantel tests were performed using the pairwise matrices of geographical distances and Jaccard distances of AFLP data in the software package Isolation by Distance Web Service (Jensen et al. 2005). Jaccard distances were used because this measure does not assume that band absence indicates homology.
3. Results and discussion
(a) Mate compatibility analyses
A total of 227 mating trials were scored for copulation success (figure 2). Overall, 54.2% of trials between gilberti individuals (n=119), 55.0% of trials between skiltonianus individuals (n=21) and only 5.7% of trials between individuals of different ecomorphs (n=87) resulted in successful copulation. As expected if parallel speciation has occurred, mating is strongly assortative by ecomorph (IPSI =0.87±0.06, p<0.0001).
Only one of the three main factors examined, body size difference relative to the optimal size difference between male and female, was a strong predictor of copulation success. All models including a body size term (table 1, models A–E) provided a substantially better fit to the data than models lacking it (table 1, models F–H). Furthermore, when divergence in body size was included in the model, inclusion of genetic and geographical distance did not improve model fit (table 1). The regression coefficient for body size was negative, while for genetic and geographical divergence it was non-distinguishable from zero (i.e. the 95% CI included zero). However, the coefficient for phylogenetic distance was strongly biased towards negative values (table 1), suggesting that more closely related individuals may have a higher probability of mating after differences in body size have been accounted for. Thus, we conclude that mate compatibility in the P. skiltonianus complex is highly dependent on body size divergence, is not affected by the geographical separation of populations and may increase slightly with increased relatedness.
Inclusion of two other factors, intercept differences across trial series and individuals, provided large improvements in model fit (table 1, compare models A, C and D). This shows that there were differences across trial series and that individuals varied in their propensity to mate, even after differences in body size were accounted for. Two experimental factors, trial order and hormone use, did not improve model fit (table 2 in electronic supplementary material), indicating that results were consistent across the 5 days of each trial series and across trials with and without hormones.
We further tested whether body size divergence influences the probability of copulation in all trial types by asking whether gilberti–skiltonianus pairs showed results consistent with those from gilberti–gilberti and skiltonianus–skiltonianus pairs. Allowing different slope coefficients for each trial type (table 1, model B) provides a small improvement in model fit (compare model C), suggesting that there are slight differences in the effects of body size divergence in different trial types (see also figure 2). However, this model confirms the universality of the importance of size difference, because estimates of all three slope coefficients are negative (table 3 in electronic supplementary material). Thus, all evidence points to a decline in copulation success as body size divergence increases, regardless of what kind of pair of individuals is being considered.
The probability of successful copulation was maximized at 73% when males were 5.6 mm larger than females from snout to vent (the Ω parameter; ±0.9 mm, 95% CI=3.8–7.4 mm; table 1). Copulation success diminished rapidly as the difference in body size increased, falling to less than 5% once the size difference between paired individuals deviated from this optimum by 15.0 mm (figure 2a). As expected given these results, heterotypic matings resulted from pairing the largest P. skiltonianus with the smallest P. gilberti; all five trials involving heterotypic pairs that differed from the optimum size difference by less than 10.0 mm resulted in successful copulation. Examination of museum specimens (figure 2b) shows that individuals involved in successful heterotypic matings (mean skiltonianus size=74.1±5.2 mm; mean gilberti size=77.7±5.7 mm) were at the extremes of the size range for their ecomorph.
Analysis of museum specimens also confirmed that mating between the geographically adjoining Sierran (Si) and southwestern (Sw) clades of P. gilberti is not constrained by ecomorph divergence; body size differences fall within a range in which successful copulation is common (δSiSw=2.38±4.31 mm; δSwSi=−2.32±4.70 mm; subscripts indicate male and female clade identity, respectively; overall, 66.7% of values are within 5.0 mm of the optimal divergence and 95.3% are within 10.0 mm). By contrast, 99.5% of the opposite-ecomorph comparisons surpass the 15.0 mm size threshold outside of which copulation success is very limited (δsg=28.67±3.18 mm; δgs=−22.03±3.11 mm). Thus, while the degree of body size divergence makes mating between P. gilberti and P. skiltonianus highly unlikely, size-dependent reproductive barriers do not prevent mating among P. gilberti clades in the wild.
(b) Divergence date estimates
Since other lineages in which parallel speciation has been documented are thought to have diverged relatively recently (Nagel & Schluter 1998; West-Eberhard 2003), we tested whether this pattern also holds true for P. gilberti. Using the 8.25×10−3 substitutions per site per my rate calibration, the Bayesian relaxed clock method places the most recent common ancestor of the three gilberti lineages at 12.7±0.1 Myr ago (95% CI =8.8–17.0 Myr ago), and the divergence separating the Inyo and southwestern gilberti lineages at 5.5±0.1 Myr ago (95% CI=4.1–7.1 Myr ago; figure 1). These dates provide maximum age estimates for the initial loss of reproductive compatibility between gilberti lineages. The southwestern clade is inferred to have diverged from its closest skiltonianus relative 4.6±0.1 Myr ago (95% CI=3.2–6.4 Myr ago), while divergence between Inyo gilberti and its sister lineage is more recent (1.7±0.02 Myr ago; 95% CI=1.1–2.5 Myr ago). These dates provide minimum age estimates for the restoration of mate compatibility between gilberti lineages. Full results, including those for the 7.15×10−3 calibration rate, are reported in electronic supplementary material, table 1.
(c) Gene flow estimates in natural populations
While the mate compatibility experiments suggest that gene flow is possible, the age of divergence makes it likely that postzygotic barriers to gene flow have evolved (Coyne & Orr 2004). Therefore, we used AFLP markers to examine genetic population structure across the mtDNA contact zone separating the Sierran and southwestern gilberti lineages. One hundred fifty-seven markers between 91 and 465 bp were unambiguously scored in the initial screening, 105 of which were polymorphic in P. gilberti (n=126 individuals). The total pool of polymorphic loci shows a positive association between levels of genetic subdivision and geographical separation across the contact zone (figure 3b). Additionally, values of θB for nearest-neighbour populations spanning the contact zone were non-distinguishable from those of nearest-neighbour populations outside the contact zone (figure 3a). Neither of these results is expected under a model in which lineages are reproductively isolated; instead, both support a model of gene flow across the contact zone. Analyses comparing pairwise AFLP distances and geographical distances separating individuals also showed a positive linear relationship between genetic structuring and geography (Jaccard similarity index; r=0.3048, P=0.001, Mantel test). Taken together, these results are consistent with a stepping-stone model of genetic isolation by geographical distance (Slatkin 1993), and suggest that formerly independent gilberti lineages in the southern Sierra Nevada are collapsing into a single biological species as a result of geographical contact following parallel speciation.
One question that arises is why the mtDNA boundary persists even when there are no obvious barriers to gene exchange. The most probable explanation is that skinks have low dispersal capabilities and strong female philopatry (Rutherford & Gregory 2003). Both theoretical and organismal examples have demonstrated that short-distance dispersal can lead to a pattern of isolation by distance when data are averaged over a large number of markers, while maintaining deep phylogeographic breaks in any individual marker (such as mtDNA; Avise 2000; Wake & Jockusch 2000; Irwin 2002; Irwin et al. 2005). Thus, patterns of contemporary gene flow may cause a merger with respect to nuclear genes while simultaneously preserving the deeper history that persists in the mtDNA phylogeny. In theory, the discordance between genetic markers could instead be explained by the repeated, unidirectional introgression of P. skiltonianus mtDNA into P. gilberti. This would obfuscate the true species boundaries in the mtDNA phylogeny and cast doubt on our conclusion that the size-dependent isolation has evolved in parallel. This alternative is implausible for two reasons: it would require that size divergence did not impose the same constraints on premating isolation in historical populations, and it fails to explain the tail colour transition that corresponds with the mtDNA clade boundaries (figure 1). Thus, parallel speciation followed by gene flow is the best explanation for the patterns of morphological and genetic divergence in P. gilberti.
Our findings show that premating isolation in the skiltonianus group depends on the degree of body size divergence between populations, and that reproductive compatibility can be restored between lineages as a result of homoplastic evolution of similar body sizes in ecologically similar populations. In Gilbert skinks, this has led to the reestablishment of phyletic connectivity between lineages that arose through independent speciation events separated by millions of years, and thus to a loss of the species boundary that once separated them. Reversed speciation has also been documented between recently evolved forms as a result of human-influenced environmental modification (Seehausen et al. 1997; Seehausen 2006; Taylor et al. 2006), and in cases of independent polyploidization events from within diploid ancestors (Holloway et al. 2006 and references therein). The example of Gilbert skinks greatly extends the timeframe over which non-polyploid reversals of speciation have been documented, and highlights the possibility that even lineages that are at present fully isolated by prezygotic barriers may regain their ability to mate as a result of parallel phenotypic evolution in similar environments. These results may have broad evolutionary implications given that parallelism is common and closely related species show a tendency to diverge in size (Schluter & Nagel 1995). If mate compatibility can be restored over long evolutionary time-scales, reversal of speciation may be a significant, but only transiently detectable, source of reticulation in the tree of life.
This research was supported by NSF-DDI Grant DEB0308969, the Theodore Roosevelt Memorial Fellowship of the American Museum of Natural History and Sigma Xi. We thank B. Alexander, B. Branciforte, C. Feldman, R. Hansen, J. LaBonte, T. Papenfuss, M. Polihronakis and D. Wood for help with field collecting. B. Goldman provided assistance with the mate compatibility pilot study. E. Adams, K. Holsinger and P. Lewis introduced us to Bayesian regression modelling and provided helpful critiques of our analyses, and A. Drummond and E. Rólan-Alvarez advised us on the use of software. R. Colwell, H. Greene, K. Holsinger, M. Polihronakis, D. Wake and two anonymous reviewers provided valuable comments on the manuscript.