Population Differentiation and Species Formation in the Deep Sea: The Potential Role of Environmental Gradients and Depth

Ecological speciation probably plays a more prominent role in diversification than previously thought, particularly in marine ecosystems where dispersal potential is great and where few obvious barriers to gene flow exist. This may be especially true in the deep sea where allopatric speciation seems insufficient to account for the rich and largely endemic fauna. Ecologically driven population differentiation and speciation are likely to be most prevalent along environmental gradients, such as those attending changes in depth. We quantified patterns of genetic variation along a depth gradient (1600-3800m) in the western North Atlantic for a protobranch bivalve ( Nuculaatacellana ) to test for population divergence. Multilocus analyses indicated a sharp discontinuity across a narrow depth range, with extremely low gene flow inferred between shallow and deep populations for thousands of generations. Phylogeographical discordance occurred between nuclear and mitochondrial loci as might be expected during the early stages of species formation. Because the geographic distance between divergent populations is small and no obvious dispersal barriers exist in this region, we suggest the divergence might reflect ecologically driven selection mediated by environmental correlates of the depth gradient. As inferred for numerous shallow-water species, environmental gradients that parallel changes in depth may play a key role in the genesis and adaptive radiation of the deep-water fauna.


Introduction
How species form is one of the most fundamental questions in evolutionary biology. Over the past two decades considerable progress has been made in identifying the scales, mechanisms, and driving forces of species formation in terrestrial and shallow-water ecosystems (e.g. [1][2][3][4][5][6][7]). However, little is known about these processes in the deep ocean, arguably the largest evolutionary realm on Earth with few obvious barriers to gene flow.
A growing body of evidence suggests that ecological speciation, defined as "the process by which barriers to gene flow evolve between populations as a result of ecologically based divergent selection between environments" [6], may be one of the key mechanisms of species formation in marine ecosystems (e.g. [7,[34][35][36][37]). Population divergence is well known to occur along environmental gradients and may lead to the formation of new species [6,38]. Even weak selective gradients, as might occur with environmental gradients, can promote strong population divergence despite gene flow among continuously distributed populations [39]. Adaptation to local selective regimes can result in environment-phenotype mismatches such that larvae dispersing from their natal environment to a contrasting one would not survive to reproduce, effectively isolating populations [31,40]. Numerous theoretical and empirical studies suggest selection along environmental gradients (e.g. temperature, moisture, altitude, salinity) promotes adaptation to different suites of abiotic and biotic conditions and ultimately may impede gene flow, leading to speciation (reviewed in [6,41,42]).
While considerable evidence exists for each of these mechanisms influencing population structure of shallow-water organisms, evidence of them operating in the deep sea is limited, apart from hydrothermal vents and other chemosynthetic ecosystems that have been more intensively studied (reviewed in [43,44]). Several interesting patterns have begun to emerge from the few phylogeographic studies of deep-sea organisms in non-chemosynthetic environments. The most distinctive is that isolation by depth appears to be much greater than isolation by distance [45][46][47][48][49][50]. For example, population divergence based on mitochondrial markers was much greater for protobranch bivalves separated by 3km of depth than 10,000 km of geographic distance [48,51]. Another interesting pattern is that population divergence appears to decrease with depth, suggesting that continental margins might be the primary site of adaptive radiation for deep-sea organisms [47,[52][53][54]. Probably the most surprising result to emerge is that population divergence can occur on extremely small scales despite the lack of obvious oceanographic or topographic features that might impede gene flow [55,56]. The small-scale divergence is often associated with depth differences and likely reflects the strong environmental gradients that attend changes in depth. In some cases the divergence is sufficient to be suggestive of cryptic species [57][58][59][60].
These emerging phylogeographic patterns suggest that the environmental gradients paralleling changes in depth likely play an important role in the formation of new species in deep-water ecosystems [47]. Increasing depth is associated with changes in a wide variety of environmental variables including temperature, hydrostatic pressure, oxygen, hydrodynamics, habitat heterogeneity, and the nature and amount of food [61]. Singly or in combination, these environmental changes are thought to influence the bathymetric distribution of organisms and shape many of the major macroecological patterns involving alpha and beta diversity [62][63][64]. While their potential ecological roles have long been appreciated, their influence on evolutionary processes and dynamics remains poorly understood.
Here we document patterns of connectivity and quantify the scale and geography of population divergence in a common protobranch bivalve Nucula atacellana Schenck 1939 (formerly Deminucula atacellana) distributed across a depth gradient (1600-3800 m) in the western North Atlantic. Previous work using mtDNA (16S) identified strong genetic divergence among populations above and below 3000m [47,55]. These results were surprising because there were no obvious topographic or oceanographic features that might isolate populations from different depth regimes. Moreover, the distance separating these regions is less than 100 km, very likely within the dispersal window of N. atacellana's demersal pelagic larvae. Several explanations might account for the divergence including idiosyncrasies of mtDNA (e.g. smaller effective population size, gender-biased dispersal), selection due to environmental changes associated with depth, or the presence of bathymetrically separated cryptic species. To better evaluate these alternative explanations, we quantify phylogeographic patterns using five loci, including both mitochondrial and nuclear markers. Recent work has stressed the importance of using multiple loci because mutational and coalescent stochasticity can lead to incongruent patterns among independent loci (e.g. [65][66][67][68][69]), and phylogeographic patterns often differ between nuclear and mitochondrial loci (reviewed in [24,70]). Our results indicate that Nucula atacellana has diverged across the depth gradient with very limited gene flow among bathyal and abyssal populations for more than 0.5 MY, possibly indicative of incipient speciation.

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.

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 Differentiation & Speciation in a Deep-Sea Bivalve PLOS ONE | www.plosone.org 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 χ 2 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 locusby-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 8 ) 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 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 7 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 χ 2 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.
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).

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~1 21,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).
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.

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   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][99][100]), and such evidence has been found recently in several marine taxa [52,54,[58][59][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 Differentiation & Speciation in a Deep-Sea Bivalve PLOS ONE | www.plosone.org 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][114][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][62][63]. Any of these gradients, singly or in combination, might lead to divergence, and they are frequently invoked as mediating adaptation (e.g. [118][119][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  [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][125][126][127][128][129][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 deepsea species, numerous studies have documented strong bathymetric divergence suggestive of cryptic species [55,[57][58][59][60]. In addition, divergence is consistently much greater among populations separated vertically than those separated horizontally [45][46][47][48][49][50][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 depthrelated 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][137][138]) and animals (reviewed in [139], also [140][141][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][144][145][146][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][148][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][152][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. 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.  Supporting Information Figure S1. Haplotype network for COI. Circle size indicates number of individuals possessing that haplotype. Small circles represent unsampled haplotypes required to connect the network. Squares indicate the most likely ancestral haplotype.
Haplotypes are shaded shallow and deep as in Figure 1.
(TIF) Figure S2. Haplotype network for CAL. Haplotype shape, size, and coloring are as in Figure S1. (TIF) Figure S3. Haplotype network for MAC. Haplotype shape, size, and coloring are as in Figure S1. (TIF) Figure S4. Haplotype network for DAC3. Haplotype shape, size, and coloring are as in Figure S1. (TIF) Figure S5. Haplotype network for DAC6. Haplotype shape, size, and coloring are as in Figure S1. (TIF) Figure S6. Population demographic history and migration estimates from IM for nuclear 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. Coloring of shallow and deep is as in Figure 1.
(TIF )   Table S1. PCR reaction mixtures and thermocycler conditions for amplified loci. PCRs were performed in 50µL reactions consisting of 1X GoTaq Flexi buffer with loading dye (Promega, Madison, WI), 2.5mM MgCl 2 , 2pmol dNTPs, 1.2pmol of each primer, 2µL genomic DNA, and 1 U of Taq polymerase (Promega). Conditions specific to each locus are given below; all protocols had an initial denaturation of 94°C for 3 min., 35 cycles of (denaturation at 94°C for 30 sec., annealing at the indicated temperature for 45 sec., extension at 72°C for 1 min), final extension at 72°C for 3 min., and a final hold at 4°C. (DOC)