Ethics Statement

No specific permissions were required for collection of specimens, because they were collected in international waters below the continental shelf. Collection of specimens did not involve endangered or protected species.

Cruise and specimen collection

On a research cruise in 2008 on board the R/V Endeavor, specimens were collected along a transect closely following the Gay Head—Bermuda transect sampled by Hessler and Sanders [71]. At most stations two epibenthic sled tows were conducted; sediments were sieved and sorted live on board, in a chilled room (2°C) using chilled seawater to minimize stress to organisms. Following sorting, the remaining bulk samples were preserved in 95% ethanol and kept at –20°C. Nucula atacellana specimens sorted on board were either preserved by flash-freezing or by placing in 95% ethanol, and stored at –80°C. Additional specimens were sorted from the bulk samples after the cruise, also using chilled ethanol to slow DNA degradation. N. atacellana was collected at 9 of 20 stations (Figure 1, Table 1), across a depth range of 1600–3800m.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 1. Map of sampled stations.

The red box in the inset shows the location of the depth transect along the slope, rise, and abyssal plain of the Northwest Atlantic, with station names and depths indicated. Stations are color-coded according to a genetic separation between shallow (red) and deep (blue; see Results).

https://doi.org/10.1371/journal.pone.0077594.g001

DNA extraction and locus amplification

Genomic DNA was extracted using the QIAamp Mini DNA Extraction Kit (Qiagen, Valencia, CA), using the standard animal tissue protocol with 2 sequential elutions of 100µL. PCR amplifications of mitochondrial COI and four noncoding nuclear loci (an actin intron (MAC), a calmodulin intron (CAL), and two noncoding anonymous fragments (DAC3 and DAC6)) were performed separately. Anonymous DNA markers were obtained by digestion of genomic DNA with the restriction enzyme AluI (NEB, Ipswich, MA), agarose gel selection of 1.0-1.5 kb fragments, cloning and sequencing; candidate markers were screened for polymorphism by sequencing a subset of the above specimens. Nuclear introns were selected based on a previous survey of introns in protobranch bivalves [72]. Standard PCR reaction mixtures were employed, and thermocycler conditions optimized for each locus (Table S1). In the few cases of poor amplification under these conditions, reamplification was performed with both a new negative control and reamplification of the original negative control, using nested primers where possible.

Sequencing, Heterozygote Detection, and Alignment

All successful amplifications produced single PCR bands as seen by gel electrophoresis, except for some individuals heterozygous at CAL for a 68bp indel that allowed separation of the two alleles in the gel. These alleles were gel purified and sequenced separately; all other single PCR bands were sequenced regardless of heterozygous status. Bi-directional sequencing was performed by Agencourt (a Beckman-Coulter company, Beverly, MA). The two reads for each individual were trimmed, aligned, and manually edited using Sequencher 5.0.1 (Gene Codes Corp., Ann Arbor, MI).

Individual base pairs were considered heterozygous if a clear double peak of near-equal height existed in both chromatograms, in the context of otherwise low or nonexistent background. Heterozygotes possessing alleles of different lengths (polymorphic indels) were ascertained by initially clear, single-peaked chromatograms that became almost-totally double-peaked (except for runs of a single nucleotide) while maintaining well-shaped peaks and regular spacing. The two juxtaposed sequences were deconvolved with the online program Indelligent [73]; each strand was deconvolved separately and the estimated alleles realigned to each other for editing and quality control. Remaining heterozygous positions were phased using PHASE 2.1.1 [74,75], employing the Parent-Independent Mutation (PIM) model for sites containing indels or more than two bases. Any uncertain phases were estimated with a second run using haplotypes phased with certainty 1.00 as knowns (–k option). Sequences were aligned using the CLUSTAL algorithm [76] in BioEdit with default alignment parameters. For nuclear loci, both alleles of all individuals were included in the alignment. Final alignments were trimmed and checked manually.

Basic and Within-Locus Analyses

To ensure the noncoding status of nuclear loci, they were checked for the potential presence of coding sequences by BLAST searches against the GenBank nucleotide database, GenBank’s ORF finder, the Gene Ontology database BLAST2GO [77], and AUGUSTUS [78]. Potential RNA secondary structure formation was assessed with Mfold [79].

Arlequin 3.5 [80] was used to compute basic indices and statistics for each locus separately: the number of haplotypes (Nhap), haplotypic diversity (H), and nucleotide diversity (π). Tajima’s D (tested at α=0.05), and Fu’s Fs (tested at α=0.02) were computed as basic tests of neutrality and demographic stability. Note that Arlequin excludes gapped positions (e.g. indels) when determining haplotypes; therefore, haplotype counts and diversity differ from other estimates. Based on initial indications of strong genetic separation at COI between a shallow group (stations 4–9a, 1600–2700m) and a deep group (stations 10–18/18a, 2800–3800m), indices and statistics were computed for each individual station with n≥3, for stations pooled among the shallow and deep group, and for all pooled individuals. For nuclear loci, deviations from Hardy-Weinberg equilibrium (HWE) were determined in Arlequin using default settings, computed among whole haplotypes. Within each locus, estimated recombination rates were determined between successive base pairs in PHASE [81,82], assuming a threshold of >5x background. Linkage disequilibrium (LD) among the nuclear loci was tested in Arlequin, using 1000 dememorization steps and 1,000,000 steps in the Markov chain.

As an additional test of selection, the McDonald-Kreitman test was performed by hand on COI, between the shallow and deep groups. Statistical significance was determined by computing the χ² statistic and p-value for the 2-by-2 contingency table of differences: Fixed vs. Variable and Synonymous vs. Nonsynonymous [83].

Population Clustering Analyses

Structure v2.3.4 [67] was used to determine the most likely number of populations (K) and to assign individuals to putative populations. The admixture model was employed, estimating separate α’ s for each population, and setting λ=1 (the Dirichlet parameter for allele frequencies) for all populations. Allele frequencies among putative populations were modeled as uncorrelated (discussed in 84), and the chain was run with a burn-in of 100,000 steps followed by 500,000 steps. Twenty replicate runs per K were conducted for K=1 to K=10, and Structure Harvester [85] was used to choose K using the delta K criterion [86]. Because K=1 cannot be evaluated using delta K, the method of Pritchard et al. [67] for choosing K was also calculated. Output for the chosen K was analyzed in CLUMPP [87] using greedy heuristic searches with 5000 random permutations, and resulting admixture proportions were plotted using distruct [88].

Haplotype networks were constructed for individual loci in TCS [89], treating gaps as a 5^th base (except for COI, where they represented missing data) and increasing the connection limit until all haplotypes were incorporated into a single network. For CAL, a large 68-bp indel (see Results) required two networks, one for the four “deletion haplotypes” and one for the remaining “insertion haplotypes”.

The population clustering determined by Structure (see Results) was tested in Arlequin by AMOVA on all five loci, nesting individuals within stations, and stations within the shallow and deep groups. AMOVAs on multilocus and locus-by-locus pairwise differences were calculated with significance assessed from a null distribution of 1000 randomizations; multilocus pairwise Φ_ST and Φ’_ST values (Φ_ST standardized by its maximum attainable value [90]) were also calculated and tested for significance within this AMOVA framework. To test for isolation-by-distance (IBD) within each population, we regressed Slatkin’s linearized Φ_ST against pairwise measures of (1) the log of geographic great-circle distance between stations and (2) the log of depth difference between stations, separately within the shallow and deep groups identified by Structure. Regression was performed via partial Mantel tests in Arlequin [91], removing the effect of depth on distance, and of distance on depth. Significance was assessed from a null distribution of 1000 random permutations.

Demographic and Population Genetic Analyses

The demographic history of populations was reconstructed in IM v.12.17.09 [92], using the populations determined by Structure and verified by AMOVA. The HKY mutation model was chosen for all loci, with a mutation rate for COI of 0.45%/(lineage·site·million-years), taken from an analysis of arcid bivalves by Marko [26]; mutation rates for nuclear loci were not specified. An Exponential Population Size Change Model was used because Extended Bayesian Skyline Plots (EBSP; see below and Results) indicated an exponentially growing shallow population. Separate analyses were conducted with and without COI. When COI was excluded, a mutation rate for CAL was used that had been estimated in BEAST calibrated with COI. Initial runs of 50,000 burn-in followed by 100,000 steps were performed to determine proper upper bounds for priors on population size (q), splitting time (t), and migration rates (m). Longer runs (>10⁸) employed 10 Metropolis-coupled chains with a two-step increment model (as per the manual) and a burnin of 100,000, and were continued until all ESSs > 70. Three replicate runs with different starting seeds were performed to assess convergence. Parameters were converted to “demographic units” using a heuristic generation time of 10 years. The inferred demographic history was plotted using IMfig. An Extended Bayesian Skyline Plot (EBSP; [93]) was produced in BEAST 1.7.4 [94] as further analysis of demographic history. Convergence was assessed using Tracer 1.5, and demography plotted using scripts written by J. Heled (https://code.google.com/p/beast-mcmc/downloads/detail?name=EBSP.zip&can=2&q=).

To determine the evolutionary history of populations and individuals, two applications of starBEAST [95] were used. In both applications, only the subset of fully sequenced individuals was used (n=74). Separate partitions were created for COI (single haplotype per individual) and each of the nuclear loci (both alleles included for all individuals), with substitution models, clock models, and locus trees unlinked across loci. The “SRD06” mutation model was used for COI, and each nuclear locus was given a GTR model with estimated equilibrium nucleotide frequencies and four categories of gamma-distributed rate variation. All loci were modeled with uncorrelated lognormally distributed clocks, setting the mean COI rate to 1.0 and estimating the others relative to COI. Starting trees were obtained by UPGMA, and a Yule prior was enforced with a piecewise linear population size and a constant root. Default priors were used for all parameters except for relative mutation rate priors for COI and clock mean rate priors, which were set to normal distributions with means and standard deviations of 1. Operators were tuned automatically, with weights adjusted per the BEAST manual. The MCMC chain was run for 10⁷ steps; burnin was determined with Tracer 1.5 and consensus trees obtained with TreeAnnotator 1.7.4. In the first application, a “population tree” was created by assigning each individual to the population inferred with Structure. In the second application, a genealogy of individuals was created by assigning all nine sequences of an individual to that individual.

Within-Locus Indices and Tests

Four nuclear loci and one mitochodrial locus were successfully sequenced from 95 individuals collected from 9 stations along a depth gradient from 1600-3800 m in the western North Atlantic (Table 2). Heterozygous indels were detected in all four nuclear loci: the 68bp indel in CAL was flanked by two indels of 4bp each, and the MAC intron contained six 1bp indels, one 2bp indel, and one 5bp indel. DAC3 contained a short run of TA repeats, and DAC6 contained 4 short indels (1bp, 1bp, 3bp, and 5bp). All sequences were deposited in GenBank (Accessions KC563091-KC563901, Table 2). No significant BLAST matches were found in the four nuclear loci, nor were ORFs or likely RNA secondary structures detected. Among all loci, DAC3 had the most haplotypes, followed by MAC, COI, CAL, and DAC6. Haplotype diversity was consistently very high; however, the deep group showed noticeably lower haplotype diversity at COI (Two-way ANOVA with locus and depth group as factors, p<0.001; Tukey’s post-hoc comparison of deep COI vs. shallow COI p<0.001). Tests of Hardy-Weinberg equilibrium showed no departures from neutral expectations, and tests of LD showed no disequilibrium among the nuclear loci (Table 3). For recombination, one location (bp 11–12 in MAC) showed evidence of a recombination rate 11x above background, but variance in this estimation within the PHASE run was greater than the mean.

Simple tests of neutrality revealed a few significantly negative values for Tajima’s D at several loci (Table 2), with departures from neutrality more common and negative at COI and DAC3, less so at MAC and CAL, and not detected at DAC6. The McDonald-Kreitman test for COI revealed a ratio of polymorphic nonsynonymous to synonymous sites (Pn/Ps) of 0.0482, and a ratio of fixed nonsynonymous to synonymous sites (Dn/Ds) of 0.0625; Pn/Ps< Dn/Ds implies that potential selection is negative. The Neutrality Index was 0.771 (NI, calculated as (Pn/Ps)/(Dn/Ds)), corresponding to a proportion of selected sites, α, of 1-NI=0.229. The 2-by-2 contingency table χ² statistic was 0.0513 (p=0.821), indicating that COI is not under selection.

Population Clustering Analyses

Structure runs tended to exhibit small variance at the highest and lowest Ks (K≤3 and K≥8), with larger variance at intermediate K (Figure 2A); however, significantly lower likelihood scores at intermediate Ks resulted in the clear choice of K=2 based on the Evanno criterion (Figure 2B); applying the Pritchard criterion resulted in K=3, with support for no separation (K=1) essentially zero. Although admixture proportions were relatively variable within the shallow group, shallow vs. deep individuals were generally ascribed to separate groups (red vs. blue respectively; Figure 2C). Structure analysis of just the nuclear loci produced slightly different admixture proportions, but resulted in a clear choice of K=2 by Evanno and Pritchard criteria (Figure 2 D-F), and again essentially zero support for K=1. Assignment proportions for this K=2 configuration indicated that 76 individuals belonged to the shallow group, and 19 to the deep group.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 2. Structure analysis for all loci and just nuclear loci.

A, Mean likelihood scores and standard deviations from 20 replicate runs at each K from Structure analysis of all 5 loci. B, Plot of Delta K model scores using the Evanno et al. method (2005); “Pr” indicates the probability for the best model (K=3) according to the method of Pritchard et al. (2000) and for the Evanno-selected K=2. C, Admixture proportions for the most likely grouping (K=2). D–F, the same analyses and measures for Structure analysis on nuclear loci only.

https://doi.org/10.1371/journal.pone.0077594.g002

The differential admixture of individuals at mitochondrial vs. nuclear loci was apparent in haplotype networks (Figure S1 vs. S2-S5). While all networks showed high allelic diversity, haplotypes of deep individuals were separated more in COI than in nuclear networks.

The AMOVA confirmed the Structure (K=2) results, indicating significant divergence between shallow and deep populations for each locus independently and when all loci were analyzed together, with little divergence within populations (Table 4). The congruence between nuclear loci and COI indicated that the depth-related divergence occurred across all loci and was not exclusive to the mitochondrion. Across all five loci, Φ_ST’ s were generally higher between shallow and deep pairs than among shallow pairs or among deep pairs (Table 5); the standardized Φ’_ST had the same pattern of significant pairwise values (not shown). Particularly among deep stations, significant Φ_ST’ s likely reflect small sample sizes. The significant separation of shallow and deep lineages was highly supported in the starBEAST genealogy (posterior probability 0.99-1.00, Figure 3), with no nodal support for significant substructure within either group. Although strong divergence was detected between depth regimes, we found little spatial structure within (Table 5). Isolation by depth was not detected within the shallow or deep populations; isolation by distance was statistically significant in the shallow group, but did not remain significant when the effect of depth was removed (Table 6).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 3. The starBEAST genealogy.

Branch lengths are proportional to substitutions per site, combined across loci. Bayesian posterior clade probabilities are shown if >0.75, and population assignments are colored as in Figure 1.

https://doi.org/10.1371/journal.pone.0077594.g003

Demographic History

Coalescent reconstruction of demographic history as estimated in IM using all 5 loci revealed an ancestral population that split approximately 95,000 generations in the past, resulting in largely independent shallow and deep populations along the sampled depth gradient (Table 7A, Figure 4). Although no good estimate exists for protobranch generation times, a conservative value of 10 years (an estimate within the range reported by Zardus [96]) translates to a split of 0.95 million years ago (MYA). The demographic estimates indicated an ancestral effective population size of ~412,000, a smaller deep population (N_e~121,000), and a much larger shallow population (N_e~3,998,000), comparable to the relative population sizes produced by the starBEAST population analysis (Figure 4, inset). Migration rates between populations per generation were extremely low (10^-7 -10^-8) and asymmetric with greater dispersal from the shallow to the deep population. Translation of these estimates into demographic units indicates that, of the 4 million individuals in the shallow population, approximately 5 migrate to the deep population each generation. Results were qualitatively similar for the IM analyses excluding COI (Table 7B). Effective populations sizes were lower, the splitting time was more recent and migration rates were somewhat larger (10^-6 -10^-7), with overall migration still extremely low. All three replicate IM runs for each analysis (with and without COI) produced very similar estimates with 95% HPD overlapping extensively for all parameters (not shown).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 4. Population demographic history and migration estimates from IM for all loci.

The gray box indicates the estimated effective population size (Ne) of the ancestral population. Estimated splitting time is indicated by the horizontal line. Descendant shallow and deep populations are represented above the line by polygons whose starting width is the estimated Ne just after the split and whose upper width is the estimated contemporary Ne. Curved dotted arrows represent estimated migration rates per generation, forward in time from source to destination. Demographic history estimation from starBEAST is shown in the inset, with branch thickness proportional to estimated population size. Coloring of shallow and deep is as in Figure 1.

https://doi.org/10.1371/journal.pone.0077594.g004

The EBSP analysis of shallow population history also showed a likely increase in the shallow population size from its ancestral size to a current N_e of 2 to 3 million, over approximately the last 1MY (Figure 5). Median population size of the deep group EBSP indicated population growth starting about 0.023 MYA (not shown), but the 95% highest posterior distribution (HPD) was quite large, including both zero growth and unrealistically high growth. The smaller sample size of the deep population appears to increase the error around demographic reconstruction.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Figure 5. Extended Bayesian Skyline Plot of population size through time for the shallow population.

Median (line) and 95% HPD (red shading) population size are shown.

https://doi.org/10.1371/journal.pone.0077594.g005

A strong genetic break across a depth gradient

The most striking feature of the phylogeographic analysis of N. atacellana is a sharp genetic break at 2700m across just 100m depth and 40km horizontal distance. Although the divergence is most obvious for the mitochondrial locus with no shared haplotypes between shallow and deep populations (Figure S1), the nuclear loci all show significant population structure (Tables 4 and 5, Figures S2–S5) and all analytical results were qualitatively the same when COI was excluded. The location of the break is quite similar to previous findings for N. atacellana [47,55], but more intense sampling narrowed the depth separation between shallow and deep populations to 100m. Within depth regimes (i.e. 1600-2700m, and 2800-3800m), very little genetic differentiation was detected using Structure, AMOVA, or a Mantel test for IBD. Multilocus coalescent modeling suggests that the depth divergence reflects a historical population split some 0.95 MYA, with extremely low gene flow across the break since its inception. When COI was excluded, Structure, AMOVA, and IM still detected divergence between shallow and deep populations; although the estimated splitting time was younger (0.55 MYA; Figure S6). Taken together, therefore, there is strong evidence from both nuclear and mitochondrial loci for a significant genetic break between populations in close proximity along the depth gradient, and little divergence within shallow and deep depth regimes.

Discordance between COI and nuclear loci

It is not surprising that mitochondrial COI exhibits stronger genetic divergence than the nuclear loci because its smaller effective population size should speed the effects of genetic drift and lineage sorting once gene flow is disrupted (e.g. [24,65,97]). Estimated mutation rates at our nuclear loci are on the same order of magnitude as that for COI, but these loci have not attained reciprocal monophyly. This pattern is expected to arise early in the process of speciation. Speciation, ongoing or recent, is often invoked in explaining discordance between nuclear and mitochondrial loci (e.g. [24,70,98-100]), and such evidence has been found recently in several marine taxa [52,54,58-60,101].

An alternative explanation for the stronger mitochondrial divergence is that either COI or another mitochondrial gene is under selection. Significantly negative neutrality indices (Tajima’s D and Fu’s Fs) were detected for some samples and can indicate purifying selection; however, these tests are highly sensitive to fluctuations in demographic parameters such as population size [102,103]. In particular, exponential population growth can cause negative neutrality indices, and indeed there is evidence for growth in N. atacellana, especially in the shallow population (Figure 5) where most of the significantly negative neutrality indices were detected. The McDonald-Kreitman test detected no selection on COI, but provides little insight into selection on other mitochondrial genes. Although we have no evidence of selection on COI, we cannot rule out the possibility that selection is operating on the mitochondrion and could account for the greater divergence at this locus (see [104,105]).

Significant fixed differences in COI are often ascribed to cryptic species, which are commonly revealed when morphologically identified species are analyzed genetically (e.g. [23,46,57,70,106,107]). In N. atacellana the four nuclear loci analyzed display high allelic diversity and some differentiation by depth; however, full sequences of the nuclear small ribosomal subunit (18S) and a 718 bp fragment of the large subunit (28S) were 100% identical among shallow and deep individuals (data not shown). Although not conclusive evidence, these results do suggest that populations have not been isolated long enough for divergence to accumulate in these more slowly evolving genes, indicating that populations of N. atacellana may be at a very early stage of species formation.

What is disrupting gene flow across mid-rise depths?

The distance between the shallow and deep groups (100m depth, 40km distance) is almost certainly within the dispersal window of N. atacellana, which has demersal pelagic larvae that likely spend days to weeks dispersing [96,108]. The amphi-Atlantic distribution of N. atacellana [109] and the lack of genetic divergence across the North Atlantic [48] suggest dispersal distances are probably quite large, as has been found in other deep-sea taxa [110,111]. If connectivity between depth regimes is not limited by distance, then either hydrographic forces or selection (presumably at unsampled mitochondrial or nuclear loci) might be precluding gene flow.

The Deep Western Boundary Current (DWBC) flows south/southwestward in the immediate vicinity and depth of our sampled region, underneath and counter to the Gulf Stream [112,113], providing a possible isolating force to populations on either side. However, while the mean flow of the DWBC is southwest, highly complex small-scale variation is pervasive and probably more important for understanding actual particle trajectories and dispersal of largely passive invertebrate larvae. Drogues released at depth at three-month intervals over three years revealed significant submesoscale coherent vortexes (SCVs), long-lived eddies propagating from the DWBC and departing from its time-averaged southward trajectory [114]. The DWBC also interacts with the Gulf Stream, creating complex, variable, and non-isobathic water movements that could transport larvae from one side of the DWBC to the other [113-115]. Lagrangian simulations of particle releases in the DWBC show high potential for mixing and transport in the sampled region [115], making it unlikely that the DWBC is an effective barrier to gene flow. It is possible that larvae transported in SCVs from relatively cooler abyssal depths into warmer rise/slope depths or vice-versa face environmental challenges (see below) reducing or eliminating population connectivity through phenotype-environment mismatches (sensu [6,31]). It is also probable that the DWBC has waxed and waned through time in response to shifting climate [116,117] and was much stronger at times over the last 0.98 MY (i.e. through the glacial/interglacial cycles of the Pleistocene), possibly initiating the observed split in N. atacellana.

The lack of a clear isolating barrier and the extremely small scale over which divergence occurs suggest that selection might play an important role. A number of environmental gradients parallel changes in depth including temperature, oxygen, salinity, POC-flux, pressure, sediment characteristics, flow regimes, and topographic complexity as well as a suite of faunal characteristics such as the diversity, composition and trophic complexity of sediment communities [61-63]. Any of these gradients, singly or in combination, might lead to divergence, and they are frequently invoked as mediating adaptation (e.g. [118-120]), delimiting bathymetric distributions [62,121], or fostering population divergence and speciation [47,53,63,122,123]. Even weak environmental gradients can initiate divergence [38], and smoothly varying gradients can create sharply divergent taxa with deep phylogenetic splits [39]. In fact, the greater divergence at mitochondrial genes is exactly what we might expect if depth-related selection on mitochondria limited gene flow between depth regimes.

Identifying the precise environmental forces that shape bathymetric patterns of genetic variation will require considerably more research, but the greater divergence in COI compared to the nuclear loci is consistent with depth-related selection on mitochondrial variants. Metabolic processes might be especially sensitive to various depth-related environmental gradients (e.g. temperature, pressure, oxygen, etc.) leading to selection for different mitochondrial variants along the depth gradient. If this selection was strong enough to impede gene flow by selecting against migrants from contrasting depths (e.g. immigrant inviability [40]) it could account for the discordance between mitochondrial and nuclear loci as well as the greater divergence of COI.

A consensus is emerging for both shallow and deep organisms that strong differences among populations from different depths may be caused by environmental gradients that parallel depth (e.g. [37,47,48,54,59,101,124-130]). For example, depth related divergence between populations of the coral Eunicea flexuosa appears to be related to strong environmental selection against ecophenotypes from contrasting depths that reduces gene flow and may ultimately lead to speciation [37]. Similar inferences were made for another shallow-water Caribbean coral Favia fragum [124] and for several deep-water corals [101,125,129,131]. In the coral Seriatopora hystrix, reciprocal transplants of depth-segregated, genetically distinct ecotypes implicated post-settlement selection against migrants from parts of the reef formation at different depths [123,132]. Even pelagic species exhibit depth-related divergence that likely reflects environmental gradients that parallel depth [54]. A rapidly growing body of evidence suggests selection along environmental gradients can lead to speciation despite continued dispersal (reviewed in [6]).

Consistent with depth and its attendant environmental gradients playing an important role in diversification of deep-sea species, numerous studies have documented strong bathymetric divergence suggestive of cryptic species [55,57-60]. In addition, divergence is consistently much greater among populations separated vertically than those separated horizontally [45-51]. For example, genetic divergence in the amphipod Eurythenes gryllus was much greater across a 3.6 km depth gradient than across 4000km at the same depth [45] or even between the Atlantic and Pacific [46]. Finally, we often find closely related congeners separated bathymetrically (e.g. [52,54,109,130,131,133,134]), precisely the pattern that would emerge if species formation was frequently mediated by depth-related environmental gradients. Depth is the most frequently observed habitat difference between sibling species [23].

The environmental gradients imposed by increasing depth in the oceans make an intriguing parallel to altitudinal gradients in terrestrial systems, with greater depth analogous to greater altitude. Numerous studies have documented gradients of lower genetic diversity with increased altitude in plants (reviewed in [135], also [136-138]) and animals (reviewed in [139], also [140-142]). While correlations in terrestrial systems are not always negative or linear, they are frequently accompanied by significant differentiation of highland and lowland clades [143-147], and often implicate the greater importance of vertical vs. horizontal distance. This was exactly the pattern documented in a widespread passerine bird in the Peruvian Andes, which was attributed to altitudinal shifts in selection on mitochondrial variants [148]. There is some evidence that, at least for animals, increased hypoxia at high altitude drives genetic differentiation and isolation-by-altitude [147-149], although adaptive changes in reproductive characteristics have also been found [150].

Another possible explanation for the depth divergence is that it formed in allopatry and the shallow and deep groups experienced secondary contact in the western North Atlantic, resulting in differential introgression of mitochondrial and nuclear genes. However, N. atacellana is widely distributed in the Atlantic, with virtually no divergence between the eastern and western North Atlantic and only modest divergence between the North and South Atlantic [48]. In addition, a similar depth divergence occurs within the Argentine basin but involves different haplotypes. Unfortunately, only formalin-fixed samples are available for the South Atlantic, restricting genetic analyses to mitochondrial loci. We cannot exclude the possibility that divergence was allopatric, but pan-Atlantic phylogeographic analyses indicate the greatest divergence is between shallow and deep groups in the western North Atlantic.

The deep ocean is a vast semi-continuous ecosystem that supports a highly diverse and largely endemic fauna. The evolutionary processes that gave rise to this distinctive fauna, the spatial and temporal scales over which they operate, and the geography and bathymetry of divergence are poorly understood. Given the limited ecological opportunity and the lack of obvious mechanisms that would allow population differentiation and speciation, it is unclear how new species form, especially at a rate sufficient to explain the high levels of diversity. Unraveling how and where evolution unfolds is critical for explaining biogeographic patterns of diversity [63], predicting how deep-sea ecosystems might respond to climate change [151-153], developing conservation and management strategies to mitigate the intense exploitation of deep-sea resources [64,154] and identifying appropriate locations and scales for MPAs [155,156]. Widespread and consistent divergence across depth gradients suggest depth and its concomitant environmental gradients may provide one of the primary mechanisms mediating population differentiation and speciation, especially below the continental shelves.