One of the most endangered vertebrates, the Devils Hole pupfish Cyprinodon diabolis, survives in a nearly impossible environment: a narrow subterranean fissure in the hottest desert on earth, Death Valley. This species became a conservation icon after a landmark 1976 US Supreme Court case affirming federal groundwater rights to its unique habitat. However, one outstanding question about this species remains unresolved: how long has diabolis persisted in this hellish environment? We used next-generation sequencing of over 13 000 loci to infer the demographic history of pupfishes in Death Valley. Instead of relicts isolated 2–3 Myr ago throughout repeated flooding of the entire region by inland seas as currently believed, we present evidence for frequent gene flow among Death Valley pupfish species and divergence after the most recent flooding 13 kyr ago. We estimate that Devils Hole was colonized by pupfish between 105 and 830 years ago, followed by genetic assimilation of pelvic fin loss and recent gene flow into neighbouring spring systems. Our results provide a new perspective on an iconic endangered species using the latest population genomic methods and support an emerging consensus that timescales for speciation are overestimated in many groups of rapidly evolving species.
Simple environments, isolated habitats with minimal biotic interactions, serve as natural laboratories for examining the extremes of ecological and evolutionary processes [1,2]. From adaptive radiation on remote depauperate islands [3–5] to laboratory microcosms [6–8], simple ecosystems may be isolated from complex community-level processes and often provide unique opportunities for testing theory. For example, what are the simplest environments in which speciation can occur without physical barriers [9–12]? How does effective population size affect the risk of extinction or the potential for adaptation [13–16]? Small inbred populations in simple natural environments may also provide a window into the future of managed endangered populations [17–19].
Devils Hole is the smallest natural range of any vertebrate species, comprising a 3.5 × 22 m opening into a flooded fault cavern greater than 130 m deep (figure 1d; [21,22]). The census population size of Cyprinodon diabolis within Devils Hole has varied from only 35–548 individuals since semi-annual monitoring began in 1972 , reaching its record low in 2013 and second lowest point (38) in 2006–2007. In addition to limited spawning habitat (a single 2.6 × 6.1 m shelf : figure 1d), the population is limited by the productivity of Devils Hole, which receives no direct sunlight for two months each year and produces 200 times less algal biomass than neighbouring springs . The population size typically decreases in response to reduced winter food availability and rebounds in the spring. Year-round water temperatures in Devils Hole remain near 32°C ; however, broad thermal tolerances in Cyprinodon from −1.9°C to 45.5°C [25–27], such as Cyprinodon pachycephalus which completes its entire life cycle at 44°C , suggest that diabolis should tolerate higher water temperatures (contra ).
The improbability of a vertebrate species surviving at such low population sizes in a suboptimal habitat, even for a short time, would seem to defy basic principles of conservation biology. For example, the probability of diabolis extinction within 50 years exceeds 80% based on historical census population sizes . Initial time-calibrated mitochondrial phylogenies suggested that Death Valley pupfishes were relicts isolated for 2–3 Myr [30–32] and estimated the age of the diabolis mitochondrial lineage to be 500 000 years [30,31]. However, early studies did not detect any unique mitochondrial haplotypes in diabolis relative to other Death Valley pupfishes , consistent with more recent colonization. Geological evidence strongly indicates that Devils Hole opened to the surface only 60 000 years ago based on the abrupt cessation of mammillary calcite deposition , which places an upper bound on the isolation of diabolis within this extreme habitat. Previous studies attributed overestimation of species divergence times to incomplete lineage sorting [30,31]; however, this cannot account for such a large discrepancy [35,36]. Instead, this conflict appears to be caused by 11 Myr priors placed on the root age of Cyprinodon + Megupsilon + Cualac [30,31] based on the assumption that pupfish populations were isolated by the formation of ancient land barriers . Thus, the discrepancy between geological and mtDNA evidence for the age of diabolis is not due to any underlying conflict in the data (e.g. time-calibration is extremely sensitive to priors ), but rather the assumption that Cyprinodon cannot disperse over land. Devils Hole was never connected to neighbouring springs by surface flow, thus the current consensus is that diabolis colonized approximately 10–20 000 years ago when water levels were higher [22,31].
A recent study estimated the age of diabolis from the coalescence of 12 microsatellite loci using approximate Bayesian computation without assumptions about overland dispersal and placed the age between 90 and 14 400 years (95% credible intervals in table 1: ). However, additional uncertainty in microsatellite mutation rates (10−3–10−6) was not accounted for by use of a point estimate from mammals (10−3), nor was the known failure of a molecular clock for microsatellite mutation rates . Thus, despite substantial attention, diabolis age estimates still range from 90 to 20 000 years [38,40]. Furthermore, additional outstanding questions about the survival of this species within its habitat remain: was Devils Hole colonized only once? Did phenotypic plasticity facilitate adaptation to this harsh environment ? How did so many small Death Valley pupfish populations survive inbreeding depression if they were isolated for millions of years ?
We examined population structure and gene flow among Death Valley populations, estimated the pupfish mutation rate from a new time-calibrated phylogeny containing all major Cyprinodon lineages, and then used this rate to calibrate a demographic model for the divergence between two closely related pupfish populations in Death Valley and approximate an age for diabolis. We sampled all extant Death Valley and Ash Meadows species and subspecies as well as outgroup species from across the Cyprinodon clade (electronic supplementary material, appendix S1). We conducted all analyses at the population level (i.e. each desert spring) and for clarity we refer to these populations by only their subspecies designation and desert spring name. In some cases, the same subspecies was sampled from multiple desert springs and we treated these as different populations. We used double-digest restriction-site-associated DNA sequencing (ddRADseq ) to genotype over 13 000 loci in 56 individuals. This provided 1.3 million sites for demographic inference, 100 000 times more data than previous studies. Our analyses reshape our understanding of population isolation and speciation in this highly endangered group of extremophile fishes and suggest a surprisingly rapid timescale for speciation, genetic assimilation and the evolution of intrinsic reproductive incompatibilities in this group.
2. Material and Methods
(a) Sample collection and bioinformatics
Cyprinodon diabolis is one of the most endangered fish on earth and thus collecting tissue from live animals was impossible. Highly degraded samples putrefied by the 32°C water were obtained from dead specimens collected by National Park staff. Other Death Valley samples used archived specimens from a previous study . DNA was extracted with Qiagen Blood & Tissue kits and ddRADseq libraries  were prepared as described previously . One hundred and ninety-eight million Illumina HiSeq 2000 single-end 100 bp reads were aligned to the Cyprinodon variegatus reference genome (v. 1.0) and called using the Stacks pipeline . Loci genotyped in at least six of the 56 high-coverage individuals and sequenced to a minimum depth of eight reads were used for analyses. Overall mean coverage depth per locus per individual for aligned reads was 138×. Further details are provided in the electronic supplementary material, Supplemental Methods.
(b) Population genetic structure and analyses of gene flow
Population genetic structure among populations was evaluated by principal components analysis of genetic variance and Bayesian unsupervised clustering using Structure . Gene flow among populations was estimated using four complementary methods: formal tests of introgression using ABBA/BABA ratios (D-statistics) , ancestry proportions across independent Structure runs, non-zero migration parameters in our dadi model, and by using Treemix  to estimate secondary gene flow among branches in a maximum-likelihood tree based on the residual covariance in allele frequencies among populations. Although we sampled small numbers of individuals per population, this will not necessarily bias our estimates of genetic differentiation owing to the large number of genetic markers sampled per individual .
(c) Phylogenetic analyses and time-calibration
We used BEAST (v. 1.8.1; ) for phylogenetic inference from 16 567 concatenated 100 bp loci based on a separate Stacks pipeline sampling the most divergent lineages across Cyprinodon. Estimation of divergence times and mutation rates is highly sensitive to the choice of priors  and requires internal calibrations [36,52,53], thus we chose the only well-defined internal calibration event known for Cyprinodon (options discussed in ): the 8000 ± 200 year age of Laguna Chichancanab [55,56], an endorheic basin which contains an endemic species flock of Cyprinodon pupfishes (figure 2d; ). It is highly unlikely that the Chichancanab species flock diverged before the basin formed based on microsatellite evidence [58,59] and because these species cannot coexist with fish predators found in all neighbouring surface waters (the flock is now extinct due to invasive fishes [54,60]). We conservatively placed a lower bound on the stem age of this clade using a normal prior for the divergence between Cyprinodon artifrons (the most closely related species) and the Chichancanab lineage with a mean of 8000 years and standard deviation of 100 years. This age and associated error (95% confidence interval (CI): ± 200 years) were estimated in geological studies of this lake using multiple core samples and multiple lines of evidence, including stable isotope data and shifts from terrestrial to aquatic invertebrate communities [55,56]. Our analysis is still highly dependent on this choice of prior and would support much older ages for diabolis if we used ancient land barriers or other mutation rates for calibration as in previous studies of this group [30,32]. For example, calibration with the human mutation rate would unrealistically place the age of the Laguna Chichancanab species flock at 4.9 Myr (see the electronic supplementary material, Supplemental Methods).
(d) Estimation of the mutation rate in pupfishes
We estimated a median global substitution rate of 5.37 × 10−7 mutations site−1 yr−1 using our original time-calibrated phylogeny for Cyprinodontidae (electronic supplementary material, table S2). This is much greater than the human genome-wide mutation rate estimated from pedigrees (0.5 × 10−9 mutations site−1 yr−1: ). Differences between human generation time (20–30 years) and pupfish generation time (less than 4 months in the wild) only partially account for this difference. It is unknown whether this reflects a bias in our RADseq data or a real biological difference in pupfishes, such as increased metabolism in extreme thermal or hypersaline environments . Alternative phylogenetic models, more stringent filtering, and taxon subsets did not substantially change our estimate of the mutation rate (electronic supplementary material, table S2). This estimate was also robust to completely rerunning the pipeline from raw reads trimmed to 53 bp to remove later positions with decreased read qualities and estimation of a new time-calibrated phylogeny (electronic supplementary material, figure S4 and table S2).
(e) Demographic modelling
Our time-calibrated phylogeny probably overestimates the age of young taxa such as diabolis owing to incomplete lineage sorting of only a single concatenated haplotype . Therefore, we used demographic modelling in dadi  to estimate divergence time from allele frequencies. Rather than using an approximate simulation of the likelihood (e.g. ), dadi directly calculates the likelihood of the joint-site frequency spectrum under various demographic scenarios using a diffusion approximation .
Owing to the extreme degradation of diabolis samples, we were not able to obtain a good estimate of diabolis allele frequencies. Instead, we bounded the age of diabolis from the divergence time between the two main population clusters in Ash Meadows: mionectes and amargosae/shoshone/nevadensis (figure 2a–c) and the older divergence of the salinus/milleri outgroup. Our phylogenetic analyses place the diabolis split in between these two branches (figure 2) or within the mionectes group in the trimmed dataset (electronic supplementary material, figure S4). In both cases, the diabolis divergence time is very close to the mionectes and amargosae/shoshone/nevadensis divergence time. We then approximated the diabolis age by correcting for the relative difference in branch lengths. Nonetheless, this approximation is a weakness of our study and should be interpreted with caution.
(a) Population structure and genetic diversity
We found evidence of genetic structure among five major groups in Death Valley: diabolis, salinus/milleri, mionectes, pectoralis and amargosae/shoshone/nevadensis (figure 2a–c; electronic supplementary material, figure S5). This is consistent with the expectation that genetic drift in small populations will accelerate lineage sorting and accumulation of genetic structure . Indeed, Death Valley populations exhibited up to 40-fold lower genetic diversity than coastal Cyprinodon populations (table 1). The lowest diversity observed in pectoralis was equal to the lowest genetic diversity ever documented in a wild population (: Lynx lynx) and 10-fold lower than a captive pupfish population extinct in the wild for over 15 years (table 1: Cyprinodon alvarezi ). Despite low diversity within each spring, if we pool all Death Valley populations into one large interconnected population, the genetic diversity represented in the entire valley is similar to high diversity coastal populations (table 1).
(b) Time-calibrated phylogeny
Our phylogeny places the stem divergence of Death Valley pupfishes at approximately 10 kyr ago (95% CI: 5–20 kyr ago; figure 2d; electronic supplementary material, figure S4). This interval is centered directly on the most recent flooding of the Death Valley basin from 10 to 35 kyr ago observed in core samples . Death Valley pupfishes shared a common ancestor with Northeastern Mexican species approximately 17 kyr ago (95% CI: 9–34 kyr ago) and Atlantic coastal populations 25 kyr ago (95% CI: 12–41 kyr ago; figure 2d). These ages are consistent with increased population mixing expected from the formation of large pluvial lakes throughout North America (e.g. Palaeolake Manly: figure 1a) during the most recent glacial period 12–25 kyr ago . It is not apparent how low-lying desert populations could have remained isolated within large inland seas because nearly all pupfish species readily interbreed on secondary contact and produce viable offspring (see Discussion).
(c) Demographic analysis of divergence time
Our maximum-likelihood estimate of the divergence between amargosae/shoshone/nevadensis and mionectes was 210 years (95% CI: 60–363; figure 3; electronic supplementary material, table S2). Encouragingly, this interval contains the Great Flood of 1862, the largest recorded rainfall for California and Nevada , which probably flooded Ash Meadows and mixed mionectes with amargosae/shoshone/nevadensis populations downstream (figure 1a). We next compared the median divergence time of diabolis relative to the median divergence between mionectes and amargosae/shoshone/nevadensis in our phylogeny (figure 2d) to rescale an approximate age for diabolis based on our demographic model and our mutation rate estimate for pupfishes (electronic supplementary material, table S2: 5.37 × 10−7 site−1 yr−1). We estimate the age of diabolis to be 255 years (electronic supplementary material, table S2: 95% CI: 105–408), suggesting that the current population colonized only 159 years before the first historical record of Devils Hole pupfish in 1893 . Alternatively, using the slowest and fastest mutation rates across the range of models and datasets we explored, diabolis ranged in age from 830 years to as young as 105 years (younger than the first historical record), respectively (electronic supplementary material, table S2: max and min 95% confidence boundaries).
(d) Recent gene flow
Unexpectedly, we also found evidence of recent gene flow between diabolis and both amargosae and pectoralis supported by formal tests for introgression (electronic supplementary material, table S4: amargosae D-statistic = 0.17, p = 0.0007; pectoralis D-statistic = 0.15, p = 0.002). Cyprinodon diabolis introgression was also consistent with Structure analyses (figure 2c) and Treemix (electronic supplementary material, figure S1). One amargosae individual contained a median of 1.6% diabolis ancestry overall (figure 2c); however, 95% credible intervals of diabolis ancestry proportions for this individual contained zero in seven out of eight independent Structure runs. Introgression from diabolis into amargosae was also supported by the residual covariance in allele frequencies between these populations estimated using Treemix (electronic supplementary material, figure S1).
The small amount of diabolis ancestry within amargosae (figure 2c), lack of significant introgression with other Death Valley and Ash Meadows populations besides pectoralis (electronic supplementary material, table S4), and directionality of gene flow inferred by Treemix (electronic supplementary material, figure S1) is consistent with gene flow from Devils Hole into Amargosa River sometime after the divergence of amargosae/shoshone/nevadensis and mionectes 210 years ago. Biologists transferred 75 diabolis to a spring in southern Ash Meadows in 1946  which might explain the pectoralis signal, but does not appear to explain the amargosae signal since no introgression was observed in the intervening shoshone population.
Overall, genetic covariance among Death Valley populations indicates additional introgression from salinus and pectoralis into nevadensis (electronic supplementary material, figure S1) and suggests overland dispersal on 100-year timescales. Similarly, our demographic model supported significant low levels of migration between amargosae/shoshone/nevadensis and mionectes after divergence (electronic supplementary material, table S3). We also confirmed a previous study  which found that a refuge population of diabolis within a concrete pond hybridized with introduced mionectes (electronic supplementary material, figure S2: 85.3% mionectes ancestry), despite its isolation from neighbouring spring systems.
Our population genomic analyses paint a dramatically different picture of survival in Death Valley pupfishes than previous assumptions. Rather than vicariance-driven speciation of isolated relicts surviving millions of years throughout flooding of the entire region lasting thousands of years (figure 1a), we found evidence for frequent overland dispersal and mixing of Death Valley pupfish populations on a scale of hundreds–thousands of years. Although overland dispersal across kilometres of desert by an aquatic organism is difficult to fathom, this view is much more consistent with established theories of population persistence and known hydrological connectivity over historical and geological timescales. For example, our time-calibrated phylogeny reasonably indicates that Death Valley pupfishes were isolated from other desert basins after the drying of the most recent pluvial lake in the region 10 kyr ago . Our demographic analysis of the most recent mixing between two desert spring populations is coincident with the largest recorded rainfall in the area in 1862  which would have flooded the dry riverbed connecting these two habitats (figure 1a). We explored a range of phylogenetic models which suggest that diabolis was isolated after salinus/milleri and only slightly earlier than mionectes and amargosae/shoshone/nevadensis, and bound its age between 105 and 830 years (min/max 95% confidence boundaries in the electronic supplementary material, table S2). This very recent origin is consistent with early reports of no unique mitochondrial haplotypes in diabolis relative to other Death Valley pupfishes in the region . However, unbiased sampling of genomic diversity by genomic resequencing and divergence time estimation from the distribution of haplotype block lengths [70,71] are needed to conclusively estimate divergence time. Finally, expectations of rapidly increasing genetic load at the low population sizes frequently observed in desert springs , particularly Devils Hole, suggest that frequent dispersal among these populations may be necessary for their long-term persistence. For example, genetic diversity in Death Valley pupfishes was up to 40-fold lower than coastal populations and included the lowest genetic diversity ever observed in a wild population (table 1: Cyprinodon nevadensis pectoralis was equivalent to L. lynx ). Nonetheless, our results suggest that the stem lineage of pupfishes in Death Valley survived at least 5000 years before the appearance of extant populations (figure 2d), possibly by maintaining a stable metapopulation despite frequent local extirpations. This perspective may inform conservation strategies for these fishes.
The young age of diabolis is surprising given its large number of unique traits: reduced aggression, sexual dimorphism and size; larger eyes; darker coloration; and the absence of pelvic fins [21,22]. For example, although Death Valley pupfishes do not show exceptional rates of morphological diversification as in the San Salvador and Chichancanab sympatric radiations (which diversified up to 130 times faster than background rates for certain jaw traits ), morphological diversification of the Death Valley clade containing diabolis was five times faster than background rates in other pupfishes , comparable to rates observed for similar traits in coral reef fishes  and Malawi cichlids  over longer timescales. Some diabolis-like traits were observed in amargosae when reared in the laboratory at 32°C on a restricted diet to simulate conditions in Devils Hole, including pelvic fin loss in up to 86% of fish . However, laboratory-rearing experiments indicate that pelvic fin loss in diabolis has a genetic basis (electronic supplementary material, figure S3). Combined, these observations suggest that colonists of Devils Hole initially lost pelvic fin expression as a plastic response to higher temperatures and reduced resources, only later evolving a broader norm of reaction for pelvic fin loss at lower temperatures. This process, known as genetic assimilation, has rarely been demonstrated in nature [74,75]. Pelvic fin loss may be adaptive in a benthic habitat, as suggested for stickleback . However, the haphazard pelvic fin loss observed in several other lineages of Cyprinodontiform fishes with no known environmental commonalities (including Empetrichthyids, Megupsilon, Orestias, two Rivulus, one Simpsonichthys and one Aphanius species) is consistent with genetic drift.
Our study adds to an emerging consensus that species ages may frequently be overestimated by time-calibrated mtDNA phylogenies [53,70,76]. There appear to be three main reasons. First, similar demographic modelling studies have consistently found that species divergence times were overestimated in mtDNA phylogenies by millions of years, e.g. in flycatchers , due to a systematic bias when sampling a single haplotype per species. Similarly, molecular estimates for the age of polar bears spanned 600 kyr–5 Myr ago using conventional demographic modelling approaches, but were revised to 400 kyr ago, consistent with the fossil evidence, by modelling with dadi . Second, many studies calibrate with distantly related outgroup fossils or ancient vicariance events, such as the breakup of Gondwanaland, assuming no subsequent dispersal was possible [30,32,77]. For example, the age of the modern bowhead whale was revised from 1.1 Myr ago to 150 kyr ago (1 order of magnitude, similar to this study) when calibrating with radiocarbon dates for ancient DNA samples rather than outgroup fossils . Third, substitution rates appear to be at least an order of magnitude slower in older lineages in both coding and non-coding regions, possibly owing to purifying selection or saturation of mutational hotspots, resulting in the overestimation of species ages when these rates are used to calibrate more recently diverged taxa [78,79]. Estimates based on only the single mitochondrial haplotype also contain substantially more uncertainty than multi-locus nDNA estimates.
Overestimation of species divergence times challenges existing assumptions about the timescale for speciation. For example, current models assume that the accumulation of postzygotic intrinsic isolation takes millions of years due to substantial work in many groups [80,81], including pupfishes [82,83]. Nearly all studies of rapid speciation in fishes report no intrinsic incompatibilities in hybrid crosses [84–87], most notably in East African cichlid clades diversifying over millions of years . These assumptions have inspired recent models of diversification, such as ephemeral speciation  and the million-year-waiting-time for macroevolutionary bursts . However, these long timescales conflict with observations of reproductive incompatibilities segregating in natural populations  which predict much shorter times for the evolution of intrinsic isolation [92,93]. Indeed, our results suggest that intrinsic genetic incompatibilities documented in two different pupfish species [94,95] evolved within less than 30 000 years. Similarly, recent studies indicate that non-lethal intrinsic genetic incompabitilities abound after secondary contact in both swordtail fishes  and hominids , suggesting that traditional progeny screens are insufficient for detecting subtle incompatibilities. Overall, our study points to the evolution of intrinsic incompatibilities, phenotypic change, and genetic assimilation on extremely rapid timescales, challenging existing models of speciation [89,90].
Our results also provide context for the management of highly inbred wild and captive populations [17,97] and the recent decline of diabolis. High egg survival rates in captivity (: 55%), rebounding population sizes , and lower genetic diversity in neighbouring stable populations (table 1) all suggest that the decline is not owing to an exceptional genetic load causing mutational meltdown as previously speculated . However, theory predicts that extinction due to increasing genetic load is highly likely for this small population over timescales of 100 years or less [13,42]. Combined with the 60 000 year age of Devils Hole, this suggests that pupfish extinction and gene flow may have happened many times before in this unique habitat and that diabolis may not be the first pupfish population to have survived there. Protecting the connectivity of this region will be essential for this cycle of rebirth to continue.
All datasets (SNP matrices) used for this study are deposited in Dryad: http://dx.doi.org/10.5061/dryad.s30t4.
C.H.M. designed the study, prepared samples for sequencing, conducted phylogenetic and population genetic analyses, and wrote the manuscript. J.E.C. ran demographic analyses with dadi and calculated D-statistics. L.H.S. provided logistical support, funding, samples, and motivation for the study. B.J.T. collected most of the samples. All authors commented on the manuscript and gave final approval for publication.
We have no competing interests.
This study was funded by a research stipend from the Miller Institute for Basic Research in the Sciences to C.H.M. and a small informal expense fund from the National Park Service.
The findings and conclusions in this article are those of the author(s) and do not necessarily represent the views of the US Fish and Wildlife Service.
We thank Anthony Echelle for supplying tissue samples from Indian Spring and the refuge diabolis population; Kevin Wilson, Bailey Gaines, and Jeffrey Goldstein for providing access to preserved diabolis specimens; Patrik Nosil, Anthony Echelle, and two anonymous reviewers for insightful comments on the manuscript. We thank the National Park Service and the US Fish and Wildlife Service for permission to conduct this research.
- Received September 29, 2015.
- Accepted January 5, 2016.
- © 2016 The Author(s)
Published by the Royal Society. All rights reserved.