Speciation‐by‐depth on coral reefs: Sympatric divergence with gene flow or cryptic transient isolation?

Abstract The distributions of many sister species in the sea overlap geographically but are partitioned along depth gradients. The genetic changes leading to depth segregation may evolve in geographic isolation as a prerequisite to coexistence or may emerge during primary divergence leading to new species. These alternatives can now be distinguished via the power endowed by the thousands of scorable loci provided by second‐generation sequence data. Here, we revisit the case of two depth‐segregated, genetically isolated ecotypes of the nominal Caribbean candelabrum coral Eunicea flexuosa. Previous analyses based on a handful of markers could not distinguish between models of genetic exchange after a period of isolation (consistent with secondary contact) and divergence with gene flow (consistent with primary divergence). Analyses of the history of isolation, genetic exchange and population size based on 15,640 new SNP markers derived from RNAseq data best support models where divergence began 800K BP and include epochs of divergence with gene flow, but with an intermediate period of transient isolation. Results also supported the previous conclusion that recent exchange between the ecotypes occurs asymmetrically from the Shallow lineage to the Deep. Parallel analyses of data from two other corals with depth‐segregated populations (Agaricia fragilis and Pocillopora damicornis) suggest divergence leading to depth‐segregated populations may begin with a period of symmetric exchange, but that an epoch of population isolation precedes more complete isolation marked by asymmetric introgression. Thus, while divergence‐with‐gene flow may account for much of the differentiation that separates closely related, depth‐segregated species, it remains to be seen whether any critical steps in the speciation process only occur when populations are isolated.


| INTRODUC TI ON
Most multi-cellular marine organisms found in shallow temperate or tropical waters possess a pelagic larval stage in their life cycle (Thorson, 1950). Such planktonic larvae may endow broad dispersal that should impede the divergence of populations and the formation of new species. Mayr (1954) proposed a strict allopatric solution to this problem. Surveying the geographic distributions of congeneric sea urchins, Mayr saw patterns that agreed "completely with terrestrial animals in which geographic speciation has been demonstrated", that is, fully allopatric distributions of most congeneric species and sympatry only for taxa deemed long-separated.
The geographic distributions of close relatives, however, do not all fit a model of strict allopatric speciation. Examining the geographic ranges of over 1,000 species of New World seaweeds, Pielou (1978) found the overlap between congeners to be greater than that of more distantly related taxa and noted an excess of congeners where the range of one species was nested within that of another. Subsequent molecular phylogenetic analyses found similar patterns (e.g. Hellberg, 1998;Hodge & Bellwood, 2016;Krug, 2011;Tavera et al., 2012;Taylor & Hellberg, 2005), again suggesting that the processes leading to species formation occurred within a relatively small region. Even some of Mayr's examples have eroded over time, with molecular analyses revealing some of the widespread "species" he examined to be composed of more geographically restricted cryptic taxa (e.g. Echinometra, Landry et al., 2003).
How then to explain the formation of new marine species with pelagic larval development without deep geographic divides? One way is to recognize that the geographic ranges of species are dynamic, such that the degree of overlap between sister species that we see today reveals little about their ranges during steps leading towards genetic isolation. Distant colonization events, geological change or the alternation of sea level or prevailing currents may produce transiently isolated populations. Divergence between these populations and their ancestral source could initiate reproductive isolation that was either continued or completed upon secondary contact. The genetic imprint of such diverge and reconnection may be well preserved in coastal species with limited dispersal (e.g. Marko, 1998). For species with pelagic larvae, such transient allopatry has been inferred by overlaying phylogenetic relationships with geographical history (Hellberg, 1998) or by polarizing phylogeographic patterns (Landry et al., 2003), but such suggestions are indirect at best.
More recently, interest has grown for an alternative scenario requiring no geographic isolation: ecological speciation, where divergent selection for performance in different habitats or on different resources produces reproductive isolation as a pleiotropic by-product (Potkamp & Fransen, 2019). Hosts (Faucci et al., 2007;Hurt et al., 2013), prey (Moura et al., 2015) and exposure to waves and predators (Rolán-Alvarez, 2007) have all been indicated as ecological drivers of such divergence, but the most common factor has been depth (Carlon & Budd, 2002;Gaither et al., 2018;Ingram, 2011;Prada & Hellberg, 2013). Indeed, depth has been noted as the most common ecological factor differentiating cryptic marine species (Knowlton, 1993).
Whether the course of speciation is completed in allopatry or sympatry, it must surmount a common barrier to differentiation: gene flow. Divergent selection can overwhelm the homogenizing effects of gene flow while allopatry escapes them altogether. Analyses that reveal the history of gene flow between populations (Sousa & Hey, 2013), then, should be able to distinguish between these alternatives. Simplistically, continuous isolation would favour an allopatric model while exchange throughout the divergence history would point to sympatric speciation. Actual results, however, do not map so easily to classic alternative models (see Bird et al., 2012). Still, emerging genomic analyses can expose the degree of isolation and exchange that have led to the emergence of genetically differentiated ecotypes (Foote et al., 2019). These analyses can also test for whether bouts of genetic exchange were symmetrical or directionally biased (Bertola et al., 2020) and whether exchange is uniform across the genome or there are "genomic islands" that experience less exchange (Rougemont et al., 2017). Finally, they can also be used to infer timing of demographic expansions (Prada et al., 2016;Takeuchi et al., 2020), that may be telling with regard to biotic and environmental changes that produce major changes in population size.
Here, we revisit a previous demographic analysis of two depth-segregated, genetically isolated ecotypes of the Caribbean candelabrum coral Eunicea flexuosa (Lamouroux, 1821). The two ecotypes, Shallow and Deep, differ in several aspects of colony and spindle morphology (Prada et al., 2008). Both also show greater rates of survival in their native depth than when reciprocally transplanted (Prada & Hellberg, 2013). The ecotypes co-occur at intermediate depths along with a small percentage of colonies of mixed heritage . Despite the small scale (<100 m) over which this habitat-associated replacement occurs, each ecotype on its own appears panmictic at a Caribbean-wide spatial scale (Prada & Hellberg, 2013). Analyses based upon one mitochondrial and three nuclear loci (Prada & Hellberg, 2013) weighed strongly against strict isolation and suggested recent exchange was primarily from the Shallow ecotype to the Deep, but could not distinguish among more subtle demographic models with varying degrees of isolation, exchange and population growth. Here, we employ >3,000 times more markers, along with more powerful demographic analyses that make use of the site frequency spectrum, to test over 100 alternative scenarios for divergence between the two ecotypes. We apply these same analyses to two other recent data sets reporting genetic differentiation between shallow and deep coral populations to test whether any trends in patterns of isolation, the timing and direction of genetic exchange, and population growth appear general. Together, our results suggest that even when gene flow connects deep and shallow populations over most of the course of their divergence, bouts of transient isolation still occur.

| Sampling, RNA extraction, library preparation and Illumina sequencing
To quantify genetic variation across depths, we sampled 59 colonies from both the Shallow (32) and the Deep ecotypes (27)  we sequenced the 5′ end of the diagnostic mitochondrial msh gene using a previous protocol (Prada & Hellberg, 2013). In addition, we sampled one individual from the closely related species Plexaura homomalla (kukenthali form) (Bayer, 1961;Sánchez & Wirshing, 2005).
We preserved all collected individuals in RNA later in at least 1:4 volume. We kept samples at 6°C in a cooler with ice for 12 hr, then changed the RNAlater and stored the samples at −80°C. Samples were shipped on dry ice to LSU and later to UC Davis, where RNA extractions were made.
We chose to use RNAseq to generate SNPs for our demographic analysis for two reasons. First, corals host many microbial symbionts, including single-celled algae that can number >30,000 per cell. DNA from these symbionts must be filtered out or they can confound demographic analyses of their host. Algal symbionts can be screened if bleached or naturally symbiont-free colonies are available to create a reference (e.g. Posbic Leydet et al., 2018), but that was not the case here. Genomic resources could be used to ensure SNPs are host specific, but these are lacking for taxa phylogenetically close to E. flexuosa. Transcriptomic data provided the best option. Second, these RNAseq data come from a larger reciprocal transplant experiment aimed at testing changes in gene expression.
To isolate high-quality mRNA, we dissected ~25 mg of tissue and used the NEB ® magnetic isolation kit following the manufacturer's instructions. After a few cleaning steps, we produced enriched poly(A) + RNA by hybridizing it with oligo (d) 25 magnetic beads. We then coupled the extraction kit with the NEBNext ® RNA Library Preparation Kit following the manufacturer's recommendations (except we used ½ reactions) with a final PCR enrichment step.
To generate an unfolded site frequency spectrum, we built a transcriptome for P. homomalla. We combined all cleaned unmatched reads from P. homomalla to assemble the transcriptome using Trinity 2.1.1 (Grabherr et al., 2011) with normalization and Trimmomatic flags. Mapping efficiency was almost identical (>85% after cleaning) for either the E. flexuosa Shallow or Deep types as sequence divergence between lineages is <1%. To avoid mapping reads against multi-copy genes, we identified single copy sequences by running OrthoFinder (Emms & Kelly, 2015) and kept those for subsequent analysis.
After assembling our mapping transcriptome, we aligned all cleaned reads to it using the BWA mem algorithm version 0.7.12 (Li & Durbin, 2009) and then generated sorted bam files using SAMtools 1.9. To avoid biases from misalignments and mapping, we conser- . The sorted BAM files generated from BWA were used to feed GATK. In GATK, we initially Add read groups, sort and create indexes using the feature AddOrReplaceReadGroups, then marked PCR duplicates with picard.jar. We then used a new GATK tool called SplitNCigarReads developed specially for RNAseq, which splits reads into exon segments (getting rid of Ns but maintaining grouping information) and hard-clipping any sequences overhanging into the intronic regions.
To call SNPs, we used Freebayes v 1.2.0 (Garrison & Marth, 2012). VCF files were then filtered using vcftools v 0.1.17 (Danecek et al., 2011) and vcflib (Garrison & Marth, 2012). We eliminated under-sequenced individuals with <50% of the SNPs and filtered SNPs within 30 bp of indels (using bcftools v 1.9-174-g4caf1fd). We retained only SNPs that were biallelic, had a minDP of 5, had a minQ of 30 and a max-missing <0.9 per population. To avoid false positives, we retained SNPs in which each allele is present in at least 20% of the reads at that site (i.e. a 20%-80% balance between variants for each SNP). We controlled for Hardy-Weinberg disequilibrium within each population using a p-value exclusion threshold of.001 and retained only one SNP per 1,000 bp (approximately one per scaffold), avoiding linked SNPs.

| Demographic inference
To quantify variation in demographic history, we used ordinary differential equations to model the evolution of allele frequencies as implemented in moments (Jouganous et al., 2017). While many loci may be subject to subtle selective divergence between close taxa (Westram et al., 2018), models in moments can include parameters for genomic islands of divergence and estimate the proportion of SNPs within them, which should both be a less disruptive way of caging potential effects of selection than removing outlier loci before estimating the AFS and provide a picture of how the proportion of loci residing in genomic islands changes through time (Roux et al., 2013).
Briefly, the models for two diverging populations can include up to three historical epochs, during which effective population size can vary from other epochs. Parameters for genetic exchange can also vary among epochs and may be symmetrical or asymmetrical. SNPs can be partitioned into two classes: those belonging to genomic islands that experience enhanced differentiation or limited exchange and the rest of the genome. Critically, multiple models not only allow for divergence with gene flow (parameters for genetic exchange in all epochs) and various forms of secondary contact (genetic exchange subsequent to an epoch of isolation) but can account for changes in population size in ancestral and post-split populations. Simulations (Momigliano et al., 2020) suggest that failing to account for such past changes in population size can bias model choice towards secondary contact.
To rank the different models according to the log-likelihoods and number of parameters, we estimated the evidence ratio from the best model following (Anderson, 2008). Specifically, for each model, we ran moments at least 10 times for each of the five bootstrap replicates and recorded the inferred model parameters with the lowest likelihood score across the 107 models. We then inferred the Akaike information criterion (AIC) scores, differences between each model against the best model, the relative likelihood of each model given the data, the model probability and the evidence ratio in favour of each model. We used the evidence ratio to rank each model and followed Anderson (2008) for model selection. Demographic curves are scaled by a per-generation mutation rate of 1.02 × 10 -9 (Prada, DeBiasse, et al., 2014;Prada et al., 2016;Voolstra et al., 2011). Assuming a conservative annual growth rate of 4 cm (Beiring & Lasker, 2000;Prada et al., 2008;Yoshioka & Yoshioka, 1991), the generation time for E. flexuosa is 5 years.

| Patterns of shared ancestry and admixture
To understand the structure of our genomic data, we first transformed our VCF into a genlight object and assessed missingness and distribution of SNPs using the R package Adeneget (Jombart et al., 2010). We then inferred population genetic structure using principal components analysis of genotypes (PCA) and inferred genetic clusters (k = 2) in adeneget v. 2. 1.1 for a dataset of 15,640 SNPs. We plotted samples using PC1 and PC2. To better understand the structure of the populations across depth gradients, we performed a discriminant analysis of principal components (DAPC) in Adegenet v. 2.1.1, partitioning the variance into that between-and within-groups to maximize discrimination. The number of clusters inferred was estimated by 1,000 iterations of the K-means clustering between K = 1 and 3 after retaining all PCs, and selecting the optimal number of PCs by 1,000 replicates of the a-score. We also inferred neighbour-joining trees using the R package phangorn v2.5.5 (Schliep, 2010).

| Additional datasets
We also explored two other systems on coral reefs for which (a) shallow and deep populations of what has been a single nominal species are genetically differentiated and (b) genomic data (>1,000 SNP loci) were available.
The first is the Agaricia fragilis system in Bermuda . Their final VCF file had been filtered for minor allele frequencies, which could distort the SFS. We used their VCF file "afra_2f.vcf", then applied filters as described for E. flexuosa. The

| RE SULTS
A total of 1.79 × 10 9 100-bp PE raw reads were generated from the 59 individuals of E. flexuosa (Table S1). After trimming and cleaning, 3.30 × 10 8 reads were aligned to the reference transcriptome of P. homomalla, resulting in the mapping of 1.65 × 10 8 reads. Transcriptomewide nucleotide diversity was high: 0.115 for the Shallow ecotype and 0.113 for the Deep. After filtering out paralogues and bad quality SNPs, then retaining only one SNP per 1 Kb, our final data set consisted of 15,640 SNPs. PCA suggests that the two ecotypes of E. flexuosa are genetically distinct (Figure 1), with Deep type individuals being more tightly clustered than Shallow ones.
The best-fit demographic model for the recent history of the E. flexuosa pair (sc3ielsm1, Figure 2) was far better supported (evidence ratio > 10 40 ) than any alternative (Table S2). The 2D allele frequency spectrum (Figure 3) showed a high density of SNPs with similar allele frequencies between Shallow and Deep populations (around the diagonal), but the presence of markers at the edges of the spectrum suggest that some regions of the genome display strong differences in allele frequencies between the ecotypes.
The narrow range of the residuals between the model and the data ( Figure S1) suggests strong power for the best model to predict the observed data.
This best-fit model ( and Figure S2). Both populations experienced population expansions beginning 221 Kya. Gene flow during this more recent epoch has been symmetrical and at levels lower than before. The best-fit model also included a term allowing differential levels of exchange between genomic islands of differentiation and the rest of the genome. This best model was 2.5 times more likely than the next-best alternative (sc12il) and over five times better than the third best (sc12imlsm2) (Table S3). These next-best alternatives differed from the best in that they included a predivergence episode of population growth.
The probability weight for these three best-fit models accounted for 85% (Table S3). The genetic distinctiveness of the two A. fragilis populations was not as great as for E. flexuosa ( Figure S1) and may contribute to the inability to identify a single, clear-cut best-fit model. Several alternative models for the demographics of reef slope and reef flat populations of Pocillopora damicornis had similar levels of support (Table S4). The best-fit model (IMisc) had a 40% probability. It suggested that divergence began in isolation 163 Kya ( Figure 5), with population sizes remaining stable since then. Genetic exchange has occurred only recently, being higher in the deep-to-shallow direction. The two next-best models, Sc2il and Sc2ilsm, suggest isolation with recent gene flow, but model Sc2i is different in that initial divergence occurred with gene flow. All models (except Sc2ilsm) favour a minor exchange in the direction from slope to flat. All four of these models are within a 10-fold evidence ratio, so distinguishing between those with initial full isolation and isolation with gene flow is equivocal. As with A. fragilis, the slope and flat populations of P. damicornis were less genetically distinctive than in E. flexuosa ( Figure S3) and this data set also contained the fewest SNPs. As for the other two species, all of the best-fit models also included a term allowing differential levels of exchange between genomic islands of differentiation and the rest of the genome.

| D ISCUSS I ON
Previous work on Eunicea flexuosa revealed that it is composed of two genetically distinct populations, each adapted to the depth where it is most common (Prada & Hellberg, 2013;Prada et al., 2008).
Demographic analysis (IMa) based on four genetic markers rejected a model of strict allopatric divergence between them. Based on these results, we suggested that this depth-segregated pair may have F I G U R E 1 Principle Components Analysis for SNP data from Deep and Shallow ecotypes of Eunicea flexuosa diverged despite continuous gene flow via immigrant inviability, enhanced by the extended time habitat-specific selection has to act between larval settlement and adulthood (Prada & Hellberg, 2013. While this scenario is plausible, subsequently emerging data, analyses and insights call it into question. First, rejecting strict allopatry need not imply a model of divergence-with-(continuous)-gene-flow, as more complicated histories of transient genetic isolation may still have occurred (Yang et al., 2017). More generally, the inferred monophyly of ecologically distinct sympatric species does not imply that they evolved in situ, as gene flow following the critical events of speciation could also produce such a pattern (Foote et al., 2019).
Here, we have revisited the history of demography and genetic interchange between the Shallow and Deep ecotypes of E. flexuosa using an analysis (moments) that allows for more complex patterns and with over 3,000 times as many genetic markers as before. This large number of markers provides a robust estimate of the site frequency spectrum, the basis for demographic analyses (like moments) that perform best at recent time scales (Patton et al., 2019). The best-fit model of historical demography and interchange between the Shallow and Deep ecotypes of E. flexuosa (Figure 3) began with a long period of divergence with gene flow, but also included an episode of isolation before more recent asymmetrical gene flow primarily from the Shallow ecotype to the Deep.
Few other marine species segregated by depth have had their histories of genetic exchange examined by genomic markers permitting a similar level of demographic detail. The best models inferred for two co-occurring Caribbean sea anemones (Bartholomea annulata Clades 1 and 2) using moments all supported divergence-withgene-flow (Titus et al., 2019), as did fastsimcoal2, although support for the latter was less than three-times better than a model with secondary contact. These two clades of anemones are not known to be depth-segregated; however, the relative abundance of Clade 1 in the murky waters off Bocas del Toro, Panama, may suggest it generally inhabits deeper waters. The best models for divergence between reef flat and reef slope populations of Pocillopora damicornis ( Figure 5) point to isolation followed by recent gene flow, although without strong enough support to draw firm conclusions.

Discriminating between alternatives was likewise difficult for the
Agaricia fragilis data (Figure 4). Modest divergence and relatively high levels of contemporary genetic exchange may leave the demographic history of divergence shrouded in such cases, especially for neutral loci .
The best-fit model for E. flexuosa also suggests a recent (9 Ky BP) population expansion for both ecotypes. Such a pattern is common to many Caribbean reef-associated species (Prada et al., 2016), from crypto-benthic fish (Eytan & Hellberg, 2010) to sea turtles (Reid et al., 2019), and is correlated with an increase in available shelf area tied to sea level rise that followed the last glacial maximum (Bellwood & Wainwright, 2002). Still, dates inferred from such analyses should be interpreted with caution. Estimating generation times for long-lived clonal organisms like E. flexuosa is complicated by the interaction between age and size: larger, older genets may contribute disproportionately to population-wide reproductive output. Selection can also influence molecular dating. RNAseq-derived markers may be under negative (stabilizing) selection. If so, genetic diversity would be lower than for truly neutral markers (Charlesworth, 1996), biasing estimates of population size downward and those for divergence times to be older (although accelerated rates of lineage sorting in regions of reduced N e would also tend to underestimate divergence times).

F I G U R E 2
Best-fit demographic model from moments analysis for the history of demography and genetic exchange between the Shallow (S) and Deep (D) ecotypes of Eunicea flexuosa. This best model (sc3imilm1) includes three epochs: an initial split, with asymmetric migration in the most recent epoch, no migration in the middle epoch, and low levels of symmetric migration in the earliest epoch. Horizontal lines are separated by 12,650 individuals. Standard deviations for date estimates are below in parentheses F I G U R E 3 2D allele frequency spectra from Eunicea flexuosa for the data (left) and the model (right). The colour scale indicates how many SNPs occur for each combination of Deep (vertical axis) and Shallow (horizontal) ecotype alleles The symmetry of gene flow also varied over the course of divergence. While genetic exchange in the early steps towards speciation was inferred to be symmetric for E. flexuosa (Figures 2 and 3), recent post-isolation exchange was asymmetrical. Similar patterns of recent asymmetric introgression have been reported for sea anemones (Titus et al., 2019), snapping shrimp (Hurt et al., 2013), oysters  and flatfish (Souissi et al., 2018). Migration was also strongly asymmetric for sympatric ecomorphs of the European whitefish in aquatic settings (Rougeux et al., 2019). Such asymmetric genetic exchange may result from differences in population size between the two lineages and unaccounted-for population expansions, but also from differences in the strength of reproductive barriers or from adaptive introgression. If barriers between the incipient species are largely complete, genes will slip through only when boosted by selection (e.g. Yang et al., 2018). Our best-fit model for E. flexuosa shows proportion of loci residing in protected genomic islands increasing greatly between the first epoch of divergence with gene flow to the most recent, a pattern also seen and associated with asymmetric gene flow in European sea bass (Tine et al., 2014).
If asymmetric migration signals adaptive introgression, then the direction of gene flow suggests that recently adaptive genes have arisen in shallow habitats more often than deep: gene flow moves from shallow populations to deep ones in depth-segregated corals, including E. flexuosa (Prada & Hellberg, 2013;Figure 2 above) and Agaricia fragilis (Bongaerts et al., 2017, Figure 4 above). A shallow-to-deep direction may indicate that genes favoured in warmer or more brightly lit surface waters have been spreading to deeper water. Notably, the Deep ecotype of E. flexuosa is more sensitive to warm water bleaching than the Shallow (Prada et al., 2010). As ocean warming continue, introgression from shallow populations may bring genotypes that allow deeper individuals to adapt to warmer oceans, although such introgression may also erode the distinctiveness of deep forms already suffering from warming conditions. Such a dilemma underlines the importance of recognizing cryptic species such as the E. flexuosa ecotypes, as ignoring differences in their susceptibility to climate change stressors can undervalue threats to local populations and ultimately to entire reef ecosystems (Fišer et al., 2018).
Depth provides a steep ecological gradient that can generate the strong differential selection required to drive divergence-with-geneflow. Despite this potential, our analysis of the Deep and Shallow ecotypes of E. flexuosa suggest that their divergence is still punctuated by an episode of transient isolation. Divergence-with-gene-flow flanks both sides of this bout of isolation: before it is symmetrical, afterwards it is asymmetrical. Are there critical changes that cement the divergence of populations that occur only (or more readily) in isolation? As of yet, we do not know the genetic basis of the ecological and reproductive traits that segregate the two ecotypes. Outlier analyses for candidate genes associated with phenotypic divergence may help identify them. Given that information, and a genomic sampling of multiple individuals, new approaches could help to map and date selective sweeps driving divergence (e.g. Tournebize et al., 2019).

F I G U R E 4
Best-fit demographic model (sc2ielsm2) from moments analysis for the history of demography and genetic exchange between shallow (S) and deep (D) populations of Agaricia fragilis based on data from ). This best model included two epochs in each population, with asymmetric migration in the first and symmetric in the second. Both populations experienced expansions beginning 55 Kya. Horizontal lines are separated by 74,894 individuals F I G U R E 5 Best-fit demographic model (IMisc) from moments analysis for the history of demography and genetic exchange between deep reef slope (S) and shallower reef flat (F) populations of Pocillopora damicornis based on data from (van Oppen et al., 2018). This best model included two epochs in each population, with symmetric migration in first epoch and asymmetric in second epoch. Horizontal lines are separated by 42,771 individuals Finally, a full understanding of divergence with depth may require knowledge of the interactions between the coral host and its consortium of symbionts. In E. flexuosa, Shallow and Deep ecotypes harbour discrete populations of Brevolium minutum (Prada, McIlroy, et al., 2014) that may help confer their capacity to occupy different light niches. Studies of similar associations with bacteria are just beginning, but we know that some bacteria co-diverge with their hosts (Pollock et al., 2018), that microbiomes vary across depth within species (Glasl et al., 2017) and that different lineages of the common coral symbiont Endozoicomonas rise and fall in correlation with environmental stressors in E. flexuosa (A. Reigel & M. E. Hellberg, unpub. data). Together, these suggest that multiple components of the coral holobiont play a role in adaptation with depth and potentially to giving rise to depth-segregated reef species.

ACK N OWLED G M ENTS
We thank Jen Roach in the laboratory of Andrew Whitehead for help in generating RNA libraries for sequencing, Diana Beltrán for help collecting, Simon Gravel for discussions on which SNPs to include in moments analyses, and two anonymous reviewers for comments that improved the manuscript.

CO N FLI C T O F I NTE R E S T
All authors declare no conflict of interest.

AUTH O R CO NTR I B UTI O N S
C. P. and M. E. H. designed the study and wrote the manuscript; C. P.
performed the field and laboratory research and analysed the data.

PEER R E V I E W
The peer review history for this article is available at https://publo ns.com/publo n/10.1111/jeb.13731.

DATA AVA I L A B I L I T Y S TAT E M E N T
Reads have been deposited in the NCBI Short Read Archive under accession numbers SAMN16478884-SAMN16478945 (BioProject PRJNA669861).