Limited genetic connectivity between gorgonian morphotypes along a depth gradient

RESEARCH ARTICLE Limited Genetic Connectivity between Gorgonian Morphotypes along a Depth Gradient Federica Costantini 1 *, Andrea Gori 2,6 , Pablo Lopez-González 3 , Lorenzo Bramanti 4 , Sergio Rossi 5 , Josep-Maria Gili 6 , Marco Abbiati 1,7 1Dipartimento di Scienze Biologiche, Geologiche ed Ambientali (BiGeA) & Centro Interdipartimentale di Ricerca per le Scienze Ambientali (CIRSA), University of Bologna, CoNISMa, Via S. Alberto 163, I-48123, Ravenna, Italy, 2Departament d’Ecología, Facultat de Biologia, Universitat de Barcelona, Av. Diagonal 643, 08028, Barcelona, Spain, 3Biodiversidad y Ecología de Invertebrados Marinos, Departamento de Zoología, Facultad de Biología, Universidad de Sevilla, Av. Reina Mercedes 6, 41012, Sevilla, Spain, 4Sorbonne Universités, UPMC Univ Paris 06, CNRS, Laboratoire d'Ecogéochimie des Environnements Benthiques (LECOB), Observatoire Océanologique, 66650, Banyuls sur Mer, France, 5Institut de Ciència i Tecnologia Ambientals, Universitat Auntònoma de Barcelona, Cerdanyola del Vallés, Spain, 6Institut de Ciències del Mar–CSIC, Pg. Maritim de la Barceloneta 37–49, 08003, Barcelona, Spain, 7Consiglio Nazionale delle Ricerche, Istituto di Scienze Marine, ISMAR, Bologna, Italy *[email protected] Abstract Gorgonian species show a high morphological variability in relation to the environment in which they live. In coastal areas, parameters such as temperature, light, currents, and food availability vary significantly with depth, potentially affecting morphology of the colonies and the structure of the populations, as well as their connectivity patterns. In tropical seas, the existence of connectivity between shallow and deep populations supported the hypothesis that the deep coral reefs could potentially act as (reproductive) refugia fostering re-colonization of shallow areas after mortality events. Moreover, this hypothesis is not so clear accepted in temperate seas. Eunicella singularis is one of the most common gorgonian species in Northwestern Mediterranean Sea, playing an important role as ecosystem engineer by providing biomass and complexity to the coralligenous habitats. It has a wide bathymetric distribution ranging from about 10 m to 100 m. Two depth-related morphotypes have been identified, differing in colony morphology, sclerite size and shape, and occurrence of symbiotic algae, but not in mitochondrial DNA haplotypes. In the present study the genetic structure of E.singularis populations along a horizontal and bathymetric gradient was assessed using microsatellites and ITS1 sequences. Restricted gene flow was found at 30–40 m depth between the two Eunicella morphotypes. Conversely, no genetic structuring has been found among shallow water populations within a spatial scale of ten kilometers. The break in gene flow between shallow and deep populations contributes to explain the morphological variability observed at different depths. Moreover, the limited vertical connectivity hinted that the refugia hypothesis does not apply to E.singularis. Re-colonization of shallow water populations, occasionally affected by mass mortality events, should then be mainly fueled by larvae from other shallow water populations. PLOS ONE | DOI:10.1371/journal.pone.0160678 August 4, 2016 1/20 a11111 OPEN ACCESS Citation: Costantini F, Gori A, Lopez-González P, Bramanti L, Rossi S, Gili J-M, et al. (2016) Limited Genetic Connectivity between Gorgonian Morphotypes along a Depth Gradient. PLoS ONE 11 (8): e0160678. doi:10.1371/journal.pone.0160678 Editor: Lorenzo Zane, University of Padova, ITALY Received: February 2, 2016 Accepted: July 24, 2016 Published: August 4, 2016 Copyright: © 2016 Costantini et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. Data Availability Statement: All relevant data are within the paper and its Supporting Information files. We have obtained the GenBank Accession numbers for all the sequences used in the work (GenBank from KX051564 to KX051577 and from KX002245 to KX002255). Funding: AG was funded by a I3P contract of the Consejo Superior de Investigaciones Cientificas (Ref. I3P-BPD2005). SR wants to thank the support of the Generalitat de Catalunya to MERS (2014 SGR - 1356) and the Marie Curie International Outgoing Fellowship (ANIMAL FOREST HEALTH, project number 327845). The study was funded by a 2003- 2004 INTERREG project between Catalonia- Introduction Marine modular organisms exhibit a large morphological variability, with phenotypic plasticity a likely source of this variability (see [1] for a review). In corals and gorgonians variability has led to considerable confusion regarding species boundaries and taxonomy. Genetic investigations can contribute to the recognition of cryptic species boundaries and population identification. Comprehensive studies combining morphological, ecological and genetic data can clarify whether genetic differentiation underlies morphological variation. Moreover, the integrative approach can help to resolve and ⁄or revise taxonomic affinities among closely related species [2]). Many environmental parameters (e.g. hydrodynamics, salinity, irradiance, trophic resources) can influence organisms’morphology and drive genetic differentiation [3]. Each parameter exerts its influence independently, but they may act synergistically and their effect is more evident along environmental gradients (e.g. latitude, depth). Indeed, it is widely accepted that the depth gradient, integrating several environmental parameters (e.g. temperature, light, hydrodynamics), may affect population structure, colony morphology, and connectivity patterns of corals and gorgonians [4–7]. Occurrence of high connectivity between shallow and deep populations gave rise to the hypothesis that deep coral reefs can potentially act as (reproductive) refugia [8], fostering the re-colonization of shallow areas after mortality events [9]. Recent studies conflict with this hypothesis, by showing that depth related patterns of genetic structuring vary according to species and sites [10,11]. Prada et al.[7] showed that colonies of the tropical gorgonian Eunicea flexuosa exposed to high water motion in shallow habitats compared to those living in deeper habitats are taller, have bigger calices, thicker branches and developed on a single plane. Shallow and deep colonies of E.flexuosa belong to two genetic lineages [7], while within lineages high connectivity was observed among populations across geographic distances of thousands of kilometers. Depth adaptive divergence generates a pattern where neutral genetic divergence is primarily partitioned by habitats, with little geographic structure [12]. Adaptation to local environmental conditions have been also suggested as a driver of genetic structuring across habitats in the tropical corals Seriatopora hystrix [10] and Montastraea cavernosa [13]. Gorgonian corals of the genus Eunicella are among the most representative ecosystem engineers providing biomass and structural complexity in sublittoral communities. Among the 18 recognized Eunicella species, 6 species are present in the Mediterranean Sea: E.singularis,E. cavolinii.E.verrucosa,E.filiformis,E.gazella and E.labiata. The distribution of the last three species is limited to the Strait of Gibraltar and to the Alboran Sea, and they are considered mainly Atlantic. E.cavolinii and E.singularis occur both in the western and eastern Mediterranean Sea, being the later also present in few locations along the Atlantic coast near the Strait of Gibraltar. E.singularis is the only species showing a high morphological variability along a depth gradient from 10 to 70 m in the western Mediterranean Sea [14,15](Fig 1;[16]). A shallow and a deep morphotype of E.singularis have been identified (see [17], [16], [18] for a detailed description of the geographical and bathymetric distribution, and the variability of the two morphotypes). The two morphotypes significantly differ in colony/sclerite shape, as well as in trophic ecology [19] and for the presence/absence of symbiotic algae [16]. Théodor [17] associated the deep morphotype to the aposymbiotic (facultative symbiotic) form of E.singularis, raising a question about its taxonomic status. Just as for the other gorgonians [20–22], mitochondrial DNA in Eunicella species has a very low mutation rate and is not a useful marker to discriminate between the two E.singularis morphotypes as well as among the different Mediterranean species of Eunicella [16,21]. To date, a detailed phylogenetic analysis within the Eunicella genus is missing, and the taxonomic status of the two E.singularis morphotypes has not been clarified yet [16]. The objective of the present study is to test whether along the Depth-Related Genetic Patterns in Eunicella singularis PLOS ONE | DOI:10.1371/journal.pone.0160678 August 4, 2016 2/20 Languedoc Rouissillon ''Pirineus Mediterrànis: La muntanya que uneix'' of the CSIC and CNRS, of the 2009-2013 European project LIFE+ INDEMARES ''Inventario y designación de la red natura 2000 en áreas marinas del estado español'' (LIFE07/NAT/E/ 000732) of the European Union, and by a PRIN 2010-2011 project (prot. 2010Z8HJ5M) on 'Coastal bioconstructions: structure, function and management' Italian Ministry of Education, Univeristy and Research. The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript. Competing Interests: The authors have declared that no competing interests exist. depth gradient morphological differences observed in E.singularis match with genetic differentiation. To this aim a population genetic approach using microsatellites and ITS1 sequence polymorphism was carried out. Moreover, divergence pattern within the genus were investigated by addressing the phylogenetic relationships among the Eunicella species. To assess the level of vertical connectivity, and the variation of the genetic parameters with depth samples of E.singularis were collected every 10 m from 10 to 60 m depth in the only site where it occurs on a vertical cliff (Cap de Creus, Northwestern Mediterranean Sea). Eight further samples were collected in the same area along about 15 km of coast at 15–20 m depth to assess the horizontal pattern of genetic connectivity, and to be able to compare the 2 patterns. Finally, the phylogenetic relationships between E.singularis, the Mediterranean species E.cavolinii and E.verrucosa, and the Atlantic E.racemosa, were analyzed to scale their genetic distances. Fig 1. Colonies of Eunicella singularis from 10 to 60 meter depth sampled in Cap de Creus; Eunicella cavolinii from Elba Island at 18–20 m depth, E. verrucosa from Tarragona at 13–17 m depth, and E.racemosa along the Morocco coast. Reprinted from [16], under a CC BY license, with permission from SPRINGER, original copyright 2012 doi:10.1371/journal.pone.0160678.g001 Depth-Related Genetic Patterns in Eunicella singularis PLOS ONE | DOI:10.1371/journal.pone.0160678 August 4, 2016 3/20 Materials and Methods Ethic statement Permissions to collect the gorgonian samples were requested to the local authorities at each site: the Parc Natural de Cap de Creus for the Cap de Creus MPA; the Societat d'Exploracions Submarines de Tarragona for Tarragona; the Parco Nazionale dell’Arcipelago Toscano for Elba Island, and the Université Ibn Zohr Agadir-Maroc for Taghazout. The field studies did not involve endangered or protected species. Study area and sampling design The study was planned at Cap de Creus (42°18’44’N; 003°19’05”E), where Eunicella singularis is the most common and abundant gorgonian species occurring at high densities in sublittoral rocky bottoms from 10 to 70 m depth [18]. Cap de Creus is characterised by vertical cliffs. However, after a careful search, only in one accessible site (Els Forcats; CCR3) E.singularis occurred along a vertical transect from 10 to 60 m depth. In this site 16 branch fragments were collected from different colonies every 10 m depth (Fig 2,Table 1). Furthermore, in eight locations along 15 km of Cap de Creus coast, branch fragments were collected by SCUBA diving from about 30 colonies of E.singularis at 18–20 m depth (Fig 2;Table 1). Finally, to analyse the genetic diversity and phylogenetic relationship in the Eunicella genus, a small branch fragment was sampled from 5 different colonies of: E.cavolinii at Elba Island (42°49’18”N; 010°09’52”E, Central Mediterranean) at 18–20 m depth, E.verrucosa at Tarragona (41°06’07”N; 001°15’12”E, Northwestern Mediterranean) at 13–17 m depth, and E.racemosa along the Morocco coast (30°32’51”N; 9°43’59”W, Taghazout, East Atlantic Ocean) at 15 m depth. All the collected material was preserved in 80–100% ethanol and stored at 4°C pending analysis. Microsatellite genotyping and ITS1 amplification Total genomic DNA was extracted from three to five polyps per colony fragments following the cetyl trimethyl ammonium bromide (CTAB) procedure [23]. After DNA extraction, all the individuals were genotyped using 6 microsatellite loci developed by [24]: C21, C30, C40 and [25]: EVER1, EVER3, EVER9 following their protocols. Genotyping of individuals was carried out on an ABI 310 Genetic Analyser (Applied Biosystems), using forward primers labelled with FAM, HEX, TAMRA, ROX (Sigma) and LIZ HD500 (Applied Biosystems) as internal size standard through MACROGEN INC. Service. Allele sizing was determined using PEAK SCANNER v.1.0 software from Applied Biosystems, Inc. For ITS1 sequence analysis, only individuals belonging CCR1, CCR2, CCR3 (10–60) and CCR4 were amplified. Polymerase chain reaction (PCR) amplifications of the ITS1 region were carried out using the primers ITS1-new-F and ITS1-new-R [26]. Each 25 μL PCR reaction contained approximately 20 ng DNA, 1× PCR buffer (Promega), 2 mm MgCl 2 , 0.5 μm of each primer, 0.4 mm dNTPs and 1 U of Taq polymerase (Promega). Amplifications were performed on a GeneAMP PCR System 2700 (Applied Biosystems) as follows: an initial denaturation at 95°C for 3 min, 30 cycles including 95°C for 30 s, 57°C for 30 s and 72°C for 60 s. A final extension at 72°C for 7 min was added. PCR products were sent to Macrogen (South Korea) for purification and sequencing since some sequences contained several heterozygous sites, ITS1 sequence types were estimated using PHASE 2.1 [27] on DNASP v. 5 [28], which implements a coalescent-based Bayesian method to infer them. Both alleles of all individuals were included in the alignment (286 reference individuals). Alignment was made manually using the biological sequences alignment editor BIOEDIT v. 7.2.5. Depth-Related Genetic Patterns in Eunicella singularis PLOS ONE | DOI:10.1371/journal.pone.0160678 August 4, 2016 4/20 Genetic variability: microsatellite and ITS Linkage disequilibrium analysis for pairs of loci was based on the likelihood ratio test with the EM algorithm [29] through 10,000 permutation procedures (number of initial conditions for EM, 100) using ARLEQUIN v. 3.5 [30]. Individuals sharing the same multilocus genotype (MLG) were checked using GENALEX v. 6.1 [31]. Identical MLGs can be the result of two different genotypes originated by two distinct sexual reproduction events but sharing the same alleles for all genotyped loci. The unbiased probability of identity (P ID )[32] was computed to test this possibility that two sampled individuals share identical MLG by chance through sexual Fig 2. Map of the sampling areas where the Eunicella species where collected, and position of the 8 locations where Eunicella singularis colonies were sampled along the Cap de Creus coast. In Els Forcats (CCR3) colonies were collected at 6 depths along a sloping rocky bottom from 10 to 60 m depth (named CCR3_10 to CCR3_60). Reprinted from [16], under a CC BY license, with permission from SPRINGER, original copyright 2012. doi:10.1371/journal.pone.0160678.g002 Depth-Related Genetic Patterns in Eunicella singularis PLOS ONE | DOI:10.1371/journal.pone.0160678 August 4, 2016 5/20 reproduction and not to clonal reproduction or sampling errors due to the collections of the same colony. Genetic diversity within site for each locus and over all loci was estimated as observed (H O ) and expected (H S ) heterozygosity using the GENETIX software package v. 4.05 [33]. Allelic richness (Ar) and private allelic richness (Ap) were calculated with a rarefaction procedure using the HP-Rare software [34]. Single- and multilocus F IS values were estimated using Weir & Cockerham’sF[35], and significant departures from the Hardy–Weinberg equilibrium were tested using “Fisher’s exact test”in GENEPOP v. 3.4 [36] as implemented for online use (http://genepop.curtin.edu.au/), with the level of significance determinate by a Markov-chain randomization. The presence of null alleles was estimated using the expectation maximization (EM) algorithm of [29]inFREENA[37]. BOTTLENECK v. 1.2.02 [38] was used to test for recent demographic changes. Significant differences between He and Heq–heterozygosity under mutation-drift equilibrium calculated using two-phase model (TPM) with 95% of the stepwise mutation model (SMM) and variance among multiple steps equal to 12 [38] were tested using Wilcoxon’s signed rank-test. A one-way full-factorial permutation multivariate analysis of variance (PERMANOVA, [39]) was used to test the null hypothesis of no difference in genetic diversity indexes (observed heterozygosity, allelic richness and private allelic richness) between depths using the 6 loci as replicates. Data were normalized (after transformation of private allelic richness as logarithm plus one), and the Euclidean distance calculated between each pair of samples. A total of 9,999 permutations were used under a reduced model. The analysis was followed by pairwise comparisons to test for differences between depths. Sequence genetic diversity within samples was estimated as number of haplotypes (h), haplotype and nucleotide diversity (Hd and π, respectively). All these parameters were calculated through DNASP. To display evolutionary relationships between sequence types (ST), ITS1 sequences were represented in a haplotype network calculated by Median Joining with the software Network v. 4.6.1.1 [40]. Population genetic structure For both molecular markers, the genetic divergence among populations was estimated using Weir & Cockerham [35]F ST estimator in ARLEQUIN. Genotypic differentiation among populations was tested with an exact test implemented in GENEPOP (Markov chain parameters: 1,000 dememorizations, followed by 1,000 batches of 1,000 iterations per batch). Since the presence of null alleles, genetic divergence among samples at microsatellite loci was estimated Table 1. Geographical coordinates of the sampling locations and depth range at which E.singularis colonies were collected. N: numbers of analysed colonies. Location Code Lat. Long. Depth (m) N Farallons CCR1 423.247 33.111 18–21 30 Es Camarallis CCR2 423.230 33.186 18 30 Els Forcats CCR3_10 423.144 33.186 10 16 Els Forcats CCR3_20 ““ 20 16 Els Forcats CCR3_30 ““ 30 16 Els Forcats CCR3_40 ““ 40 16 Els Forcats CCR3_50 ““ 50 16 Els Forcats CC3_60 ““ 60 16 Ses Ielles CCR4 423.119 32.977 17–20 30 Es Caials CCR5 422.841 32.977 15–17 30 Es Cucurucú CCR6 422.733 32.844 15–17 30 La Reparada CCR7 422.511 32.808 17 25 Cap Norfeu CCR8 422.391 32.675 20 30 doi:10.1371/journal.pone.0160678.t001 Depth-Related Genetic Patterns in Eunicella singularis PLOS ONE | DOI:10.1371/journal.pone.0160678 August 4, 2016 6/20 in FREENA using the F ST estimates of Weir [41] and following the so-called ENA method described in [37] which provides unbiased F ST estimates, computed excluding null alleles. To evaluate the number of clusters (K) in the dataset without prior information regarding the geographic distribution of the samples, a Bayesian method implemented in STRUCTURE v. 2.3.4 [42,43] was used under the admixture model and choosing the assumption of correlated allele frequencies and the option of recessive alleles to cope with null alleles, as suggested by Falush [43]. Mean and variance of log likelihoods of the number of clusters for K = 1 to K = 10 were inferred from multilocus genotypes. Ten runs were performed for each K with 500,000 iterations and a burn-in period of 50,000. In order to identify the number of clusters that best fit the data, the resulting output files were then analysed using the Evanno method in STRUCTURE HARVESTER v. 0.6.94, as implemented for online use (http://taylor0.biology.ucla.edu/ structureHarvester/)[44]. Moreover, the STRUCTURE results were summarized using CLUMPAK [45] to obtain the probability of each individual to belong to each cluster. A discriminating analysis of principal components (DAPC) as implemented in the ADEGENET software v. 1.3 [46] was performed. This technique extracts information from genetic datasets (multivariate in nature) by first performing a principal component analysis (PCA) on pre-defined groups or populations, and then using the PCA factors as variables for a discriminating analysis (DA), which seeks to maximize the intergroup component of variation. The optimal number of clusters (populations) was predicted using the k-means clustering algorithm, find clusters, retaining all principal components. In all analyses, 20 principal components of PCA were retained as input to DA, which accounted for approximately 96% of the total genetic variability. Significance levels for multiple comparisons of loci across samples were adjusted using a false discovery rate (FDR) correction for multiple tests [47]. Hierarchical analysis of molecular variance (AMOVA) was conducted in ARLEQUIN in order to quantify genetic variation between two groups depending on the depth: shallow populations and deep populations (CCR3_40, CCR3_50 and CCR_60). Significance of F-statistics was achieved using 99,999 permutations. Individual assignment tests were performed using the program GENECLASS v. 2 [48]. This Bayesian procedure computes the likelihood of a genotype in a given population assuming an equal prior probability density to the allelic frequencies of each locus in each population. Bayesian method [49] and the re-sampling algorithm of Cornuet et al.[50] set at 1,000 individuals were used. An individual was excluded from a given candidate population if its probability of belonging to a particular population was lower than 5%. To evaluate the extent of contemporary (ecological timescale) dispersal of E.singularis among sites, the numbers of first-generation migrant were also calculated through the program GENECLASS v. 2. Phylogenetic analyses To perform a preliminary phylogeny of the genus Eunicella, a dataset with all the individual ITS1 sequences from the different species here examined was created. Distance matrixes of sequence divergence among the species were calculated as p distance (Dp) in MEGA v. 5.05. The best-fit substitution model for the dataset was calculated using JMODELTEST v. 1.1 [51] Mac software, considering 88 substitution models by hLRT calculator with 4-gamma category. Phylogenetic relationship among species were carried out using a Bayesian approached implemented in MrBayes v. 3.1.2 software using the Hasegawa, Kishino and Yano substitution model. The analysis was carried out for 2,000,000 generations, sampling every 1,000 generations. The first 500 of sampled generations were discarded as the burn in. Phylogenetic tree were visualized as posterior probability using FIGTREE v. 1.4.0. ITS1 sequences of Corallium rubrum (GenBank accession number: FJ87608282 and FJ87608283) were used as outgroup. Depth-Related Genetic Patterns in Eunicella singularis PLOS ONE | DOI:10.1371/journal.pone.0160678 August 4, 2016 7/20 Results Genetic variability: microsatellite and ITS 301 individuals were genotyped for the 6 loci. Tests of linkage disequilibrium between loci within populations indicated no significant association of alleles (all p>0.05), confirming that all loci are independent markers. Six shared multilocus genotypes (MLGs) between colonies were found. Five MLGs were encountered twice in CCR3_40, CCR3_50 and CCR3_60. Within CCR3_50 one MLG was encountered five times. In CCR3_50 the probability that each of these genotypes was produced by chance alone was high (P ID = 4.6 e-02). The allelic richness (Ar) based on a minimum sample size of 22 genes, ranged from 2.24 (in CCR3_50) to 4.25 (in CCR7), while the private allelic richness (Ap) ranged from 0 (in CCR3_20 and CCR3_30) to 0.46 (in CCR7). Expected heterozygosity (H E ) and observed heterozygosity (H O ) were lower in CCR3_50 (H E = 0.25, H O = 0.28) and higher in CCR4 (H E = 0.58, H O = 0.57). Three (CCR2, CC3_10 and CCR3_30) out of 13 samples showed significant deviations from HW equilibrium with values of F IS positive and higher than 0.173 (Table 2). Estimated null allele frequencies ranged between loci from 0.022 ± 0.035 in EVER9 to 0.0748 ± 0.087 in C21. The Wilcoxon’s signed rank-test, performed in BOTTLENECK showed no demographic population changes in the recent past for all samples (data not shown). The one-way PERMANOVA showed that the total genetic variability differed by depth (MS = 6.68; F = 2.65; df = 5; p= 0.005). The pairwise test located a threshold of genetic variability at 40 m depth (10–50 t = 2.42, p= 0.0097; 20–50 t = 3.86, p= 0.0001; 30–50 t = 2.88, p= 0.0054). Samples collected in 3 shallow sites (CCR1, CCR2, CCR4) and along the vertical gradient in CCR3 (CCR3_10 to CCR3_60) were analyzed using ITS1 sequence polymorphism. In total, 286 ITS1 sequences corresponding to 143 individuals were analyzed. The length of the amplified ITS1 fragment was 200 bp with 9 variable sites at the positions 14, 21, 114, 119, 158, 159 and 173, corresponding to 9 sequence types (ST1 to ST9; Table 3). The more frequent sequence type was ST1 occurring in all populations, while the second more frequent sequence type was ST7 occurring in CCR3_40 and CCR3_60 with 9 and 5 individuals, respectively (Table 3;Fig 3). Haplotype diversity ranged from 0 in CCR3_50 to 0.653 in CCR3_60 (Table 3). No significant differences between CCR3_10, CCR3_20, CCR3_30 vs. CCR3_40, CCR3_50, CCR3_60 were observed in the haplotype diversity (Student’s t-test, Table 2. Genetic diversity of Eunicella singularis at six microsatellite loci; H O : observed heterozygosity, H E : unbiased expected heterozygosity, Ar: allelic richness, Ap: private allelic richness, F IS : Weir and Cockerham’s[35] estimate fixation index with significant values in bold (0.05 threshold after FDR correction). Site H O H E Ar Ap F IS CCR1 0.50 0.57 3.68 0.18 0.115 CCR2 0.47 0.57 3.85 0.02 0.173 CCR3_10 0.44 0.55 3.43 0.06 0.205 CCR3_20 0.52 0.52 3.21 0.00 0.07 CCR3_30 0.47 0.59 3.42 0.00 0.203 CCR3_40 0.54 0.53 3.39 0.15 -0.028 CCR3_50 0.28 0.25 2.24 0.05 -0.098 CCR3_60 0.49 0.51 3.46 0.13 0.045 CCR4 0.57 0.58 3.70 0.07 0.011 CCR5 0.55 0.58 4.29 0.21 0.053 CCR6 0.52 0.54 3.62 0.05 0.043 CCR7 0.52 0.57 4.25 0.46 0.096 CCR8 0.53 0.58 3.98 0.15 0.083 doi:10.1371/journal.pone.0160678.t002 Depth-Related Genetic Patterns in Eunicella singularis PLOS ONE | DOI:10.1371/journal.pone.0160678 August 4, 2016 8/20 p= 0.07; data not shown). The sequence type network showed a star-shape pattern with a high frequency dominant haplotype with at the center. Derivatives with lower frequencies were connected to the dominant haplotypes by only one mutation step (Fig 3). Population genetic structuring Genotypic differentiation between samples and F ST estimates based on the ENA method gave similar results to those obtained when presence of null alleles was not taken into account. Pairwise F ST values ranged from -0.002 (CCR2 vs. CCR5) to 0.482 (CCR3_50 vs. CCR3_60). After FDR correction, all the pairwise comparison with CCR3_40, CCR3_50 and CCR3_60 were significant (Table 4). During the first round of STRUCTURE, two genetic clustering were identified as more plausible, K = 2 (ΔK = 74.71) (Fig 4A). When K = 2 is considered, the first cluster includes CCR3_40, CCR3_50 and CCR3_60, while the other cluster includes all the remaining samples. When K = 4 is considered, STRUCTURE outline a new genetic pool characterized corresponding to CCR3_50 (Fig 4B). The DAPC confirmed the genetic differentiation observed using Structure. In the data transformation step for PCA analysis, 20 principal components (PCs) were retained, accounting for approximately the 95% of the total genetic variability. The eigenvalues of the DAPC indicated that the first two components explained most of the variation. Along the first axis CCR3_40, CCR3_50 and CCR3_60 were separated from the others samples; the second axis contrasts CCR3_50 with CCR3_40 and CCR3_60 (Fig 5). In fact, 93% of individuals were reassigned to their original clusters CCR3_50 and CCR3_40, while 75% of the individuals were reassigned to CCR_60. The values of the proportions of successful reassignment were below 30% in all the shallow populations suggesting admixture among them. The AMOVA supported the clustering structure showing a significant difference between shallow and deep groups: 15.39%, p<0.05 (Table 5). Of the 301 individuals included in the assignment analysis, 38.2% were assigned to their actual sampling location. In particular, all individuals from CCR3_40, CCR3_50 and CCR3_60 were assigned to their population of origin. The analysis of first generation migrants, based on the Bayesian computation of Rannala et al. [49] implemented in GENECLASS v. 2, classified 6 colonies as originating from a site different to the one they were collected at. None of the 6 colonies were identified as immigrants coming from deeper populations (data not shown). For the ITS1 sequence data set, F ST ranged from 0 (CCR3_10 vs. CCR3_20; CCR3_20 vs. CCR3_30 and CCR3_50 vs. CCR3_10, CCR3_20, CCR3_30) to 0.548 (CCR3_40 vs. Table 3. Sequence types frequencies, including heterozygous individuals in each population of E.singularis.H: number of haplotypes (in the sense of sequence type); Hd: haplotype diversity; πd: nucleotide diversity. CCR1 CCR2 CCR£_10 CCR3_20 CCR3_30 CCR3_40 CCR3_50 CCR3_60 CCR4 ST1 26 22 31 31 31 14 32 16 26 ST2 4 6 ST3 7 ST4 3 ST5 1 1 ST6 1 ST7 18 10 ST8 4 ST9 2 H232222142 Hd 0.024±0.092 0.486±0.085 0.063±0.058 0.063±0.058 0.006±0.058 0.508±0.031 0 0.653±0.057 0.315±0.087 πd 0.00120 0.00398 0.00124 0.00063 0.00031 0.00254 0 0.00415 0.00157 doi:10.1371/journal.pone.0160678.t003 Depth-Related Genetic Patterns in Eunicella singularis PLOS ONE | DOI:10.1371/journal.pone.0160678 August 4, 2016 9/20 extreme thermal stress in shallow waters has recently resulted in mass mortalities of benthic invertebrate species, including gorgonians (e.g. [69]). However, re-colonization and population’s recovery after a mass mortality event cannot be supported by larvae from deeper colonies, but can rely only on larvae produced by nearby shallow water populations. Concerning the taxonomy of the genus Eunicella in the Mediterranean Sea, low levels of mitochondrial genome evolution [20], and low polymorphism of nuclear markers, do not allow to validate the species status neither of the two morphotypes, nor of the three Eunicella species. The large morphological variability of Eunicella genus, as in other Octocorallia genera, leads uncertainty regarding species boundaries and classification. Methods based on next-generation sequencing (e.g. RAD sequencing, [70]) may provide further insights into the evolutionary relationships in recalcitrant taxa such as Eunicella. Supporting Information S1 Dataset. Microsatellite data set. The allele sizes in 6 loci studied in 13 E.singularis populations. (XLS) S1 Materials and Methods. Materials and Methods and Results. DNA extraction of the Eunicella species and PCR amplification of the three mitochondrial markers and their genetic variability. (DOCX) S1 Table. Pairwise genetic p distance among Eunicella species. Pairwise genetic p distance (Dp) among Eunicella species including the sampled populations of Eunicella singularis. (DOCX) Acknowledgments Authors are grateful to L. Terrazzoni for her collaboration in lab work, and to N Viladrich, I Fiorillo and T Garcia for their help during sampling. A Gori was funded by a I3P contract of the Consejo Superior de Investigaciones Cientificas (Ref. I3P-BPD2005). S Rossi wants to thank the support of the Generalitat de Catalunya to MERS (2014 SGR—1356) and the Marie Curie International Outgoing Fellowship (ANIMAL FOREST HEALTH, project number 327845). This work is part of the 2003–2004 INTERREG project between Catalonia-Lan- guedoc Rouissillon ‘‘Pirineus Mediterrànis: La muntanya que uneix”of the CSIC and CNRS, of the 2009–2013 European project LIFE+ INDEMARES ‘‘Inventario y designación de la red natura 2000 en áreas marinas del estado español”(LIFE07/NAT/E/000732) of the European Union, and of the PRIN 2010–2011 project (prot. 2010Z8HJ5M) on ‘Coastal bioconstructions: structure, function and management’Italian Ministry of Education, Univeristy and Research. Author Contributions Conceived and designed the experiments: FC AG JMG PLG SR LB. Performed the experiments: FC AG LB SR. Analyzed the data: FC MA. Contributed reagents/materials/analysis tools: MA JMG. Wrote the paper: FC AG PLG LB SR JMG AM. Depth-Related Genetic Patterns in Eunicella singularis PLOS ONE | DOI:10.1371/journal.pone.0160678 August 4, 2016 16 / 20 References 1. Todd PA. Morphological plasticity in scleractinian corals. Biol Rev. 2008; 83: 315–337. PMID: 18979594 2. Pinzón JH, LaJeunesse TC. Species delimitation of common reef corals in the genus Pocillopora using nucleotide sequence phylogenies, population genetics and symbiosis ecology. Mol Ecol. 2011; 20: 311–325. doi: 10.1111/j.1365-294X.2010.04939.x PMID: 21091563 3. Sanford E, Kelly MW. Local adaptation in marine invertebrates. Ann Rev Mar Sci. 2011; 3: 509–535. doi: 10.1146/annurev-marine-120709-142756 PMID: 21329215 4. Garrabou J, Harmelin J. A 20- year study on lifehistory traits of a harvested longlived temperate coral in the NW Mediterranean: insights into conservation and management needs. J Anim Ecol. 2002; 71: 966–978. doi: 10.1046/j.1365-2656.2002.00661.x 5. Carlon DB, Budd AF. Incipient speciation across a depth gradient in a scleractinian coral? Evolution. 2002; 56: 2227–42. doi: 10.1111/j.0014-3820.2002.tb00147.x PMID: 12487353 6. Costantini F, Rossi S, Pintus E, Cerrano C, Gili J-M, Abbiati M. Low connectivity and declining genetic variability along a depth gradient in Corallium rubrum populations. Coral Reefs. 2011; 30: 991–1003. doi: 10.1007/s00338-011-0771-1 7. Prada C, Schizas N V, Yoshioka PM. Phenotypic plasticity or speciation? A case from a clonal marine organism. BMC Evol Biol. 2008; 8: 47. doi: 10.1186/1471-2148-8-47 PMID: 18271961 8. Glynn PW. Coral reef bleaching: facts, hypotheses and implications. Glob Chang Biol. 1996; 2: 495– 509. 9. Bongaerts P, Ridgway T, Sampayo EM, Hoegh-Guldberg O. Assessing the ‘deep reef refugia’hypothesis: focus on Caribbean reefs. Coral Reefs. 2010; 29: 309–27. doi: 10.1007/s00338-009- 0581-x 10. Bongaerts P, Riginos C, Ridgway T, Sampayo EM, van Oppen MJH, Englebert N, et al. Genetic divergence across habitats in the widespread coral Seriatopora hystrix and its associated Symbiodinium. PLoS One 2010; 5: e10871. doi: 10.1371/journal.pone.0010871 PMID: 20523735 11. van Oppen MJH, Bongaerts P, Underwood JN, Peplow LM, Cooper TF. The role of deep reefs in shallow reef recovery: an assessment of vertical connectivity in a brooding coral from west and east Australia. Mol Ecol. 2011; 20: 1647–60. doi: 10.1111/j.1365-294X.2011.05050.x PMID: 21410573 12. Prada C, Hellberg ME. Long prereproductive selection and divergence by depth in a Caribbean candelabrum coral. Proc Natl Acad Sci U S A. 2013; 110: 3961–3966. doi: 10.1073/pnas.1208931110 PMID: 23359716 13. Serrano X, Baums IB, O’Reilly K, Smith TB, Jones RJ, Shearer TL, et al. Geographic differences in vertical connectivity in the Caribbean coral Montastraea cavernosa despite high levels of horizontal connectivity at shallow depths. Mol Ecol. 2014; 23: 4226–4240. doi: 10.1111/mec.12861 PMID: 25039722 14. Grasshoff M. Flachwasser—Gorgonarien von Europa und Westafrika (Cnidaria, Anthozoa). Frankfurt am Main: Senckenbergische Naturforschende Gesellschaft. Cour Forschunginstitut Senckenb. 1992; 149:135. 15. Carpine C, Grasshoff M. Les Gorgonaires de la Méditerranée. Bull Inst Océanography Fish. 1975; 71:1–140. 16. Gori A, Bramanti L, López-González P, Thoma JN, Gili J-M, Grinyó J, et al. Characterization of the zooxanthellate and azooxanthellate morphotypes of the Mediterranean gorgonian Eunicella singularis. Mar Biol. 2012; 159: 1485–1496. doi: 10.1007/s00227-012-1928-3 17. Théodor J. Contribution a l’étude des gorgones (VIII): Eunicella stricta aphyta sous-espèce nouvelle sans zooxanthelles proche d’une espèce normalement infestée par ces algues. Vie Milieu. 1969; 20: 635e638. 18. Gori A, Rossi S, Berganzo E, Pretus JL, Dale MRT, Gili J-M. Spatial distribution patterns of the gorgonians Eunicella singularis,Paramuricea clavata, and Leptogorgia sarmentosa (Cape of Creus, Northwestern Mediterranean Sea). Mar Biol. 2011; 158: 143–58. doi: 10.1007/s00227- 010-1548-8 19. Gori A, Viladrich N, Gili J-M, Kotta M, Cucio C, Magni L, et al. Reproductive cycle and trophic ecology in deep versus shallow populations of the Mediterranean gorgonian Eunicella singularis (Cap de Creus, northwestern Mediterranean Sea). Coral Reefs 2012; 31: 823–837 doi: 10.1007/s00338-012-0904-1 20. Shearer TL, van Oppen MJH, Romano SL, Wörheide G. Slow mitochondrial DNA sequence evolution in Anthozoa (Cnidaria). Mol Ecol. 2002; 11: 2475–87. doi: 10.1046/j.1365- 294X.2002.01652.x PMID: 12453233 21. Calderòn I, Garrabou J, Aurelle D. Evaluation of the utility of COI and ITS markers as tools for population genetic studies of temperate gorgonians. J Exp Mar Biol Ecol. 2006; 336: 184–97. doi: 10.1016/j. jembe.2006.05.006 Depth-Related Genetic Patterns in Eunicella singularis PLOS ONE | DOI:10.1371/journal.pone.0160678 August 4, 2016 17 / 20 22. McFadden CS, Benayahu Y, Pante E, Thoma JN, Nevarez PA, France SC. Limitations of mitochondrial gene barcoding in Octocorallia. Mol Ecol Resour. 2011; 11:19–31. doi: 10.1111/j.1755-0998.2010. 02875.x PMID: 21429097 23. Winnepenninckx B, Backeljau TRDW, De Wachter R. 1994 Extraction of high molecular weight DNA from molluscs. Trends Genet. 1994; 9: 407. doi: 10.1016/0168-9525(93)90102-N 24. Abdoullaye D, Acevedo I, Adebayo A, Behrmann-Godel J. Benjamin RC, Bock DG.et al. Permanent genetic resources added to molecular ecology resources database 1 August 2009–30 September 2009. Mol Ecol Resour. 2010; 10. https://dx.doi.org/10.1111/j.17550998.2009.02796.x 25. Holland LP, Dawson DA, Horsburgh GJ, Krupa AP, Stevens JR. Isolation and characterization of fourteen microsatellite loci from the endangered octocoral Eunicella verrucosa (Pallas 1766). Conserv Genet Resour. 2013; 5: 825–9. doi: 10.1007/s12686-013-9919-3 26. Costantini F, Fauvelot C, Abbiati M. Genetic structuring of the temperate gorgonian coral (Corallium rubrum) across the western Mediterranean Sea revealed by microsatellites and nuclear sequences. Mol. Ecol. 2007; 16: 5168–82. doi: 10.1111/j.1365-294X.2007.03579.x PMID: 17986195 27. Stephens M, Smith N, Donnelly P. A new statistical method for haplotype reconstruction from population data. Am J Hum Genet. 2001; 68: 978–89. doi: 10.1086/319501 PMID: 11254454 28. Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009; 25: 1451–2. doi: 10.1093/bioinformatics/btp187 PMID: 19346325 29. Dempster A, Laird N, Rubin D. Maximum likelihood from incomplete data via the EM algorithm. J R Stat Soc. 1977; 39: 1–38. doi: 10.1.1.133.4884 30. Excoffier L, Lischer H. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Resour. 2010; 10: 564–7. doi: 10.1111/j.1755-0998. 2010.02847.x PMID: 21565059 31. Peakall R, Smouse PE. Genalex 6: genetic analysis in Excel. Population genetic software for teaching and research. Mol Ecol Notes. 2006: 6: 288–95. doi:j.1471-8286.2005.01155.x 32. Kendall M, Stewart A. 1977 The advanced theory of statistics. Vol. 1. New York: Macmillan 33. Belkhir K, Borsa P, Chikhi L, Raufaste N, Bonhomme F. GENETIX 4.05, Logiciel sous Windows pour la Génétique des Populations. Montpellier, France: Laboratoire Génome, Populations, Interactions, CNRS UMR 5171, Université de Montpellier II. 2004 34. Kalinowski ST. Hprare 1.0: a computer program for performing rarefaction on measures of allelic richness. Mol Ecol Notes. 2005; 5: 187–9. doi: 10.1111/j.1471-8286.2004.00845.x 35. Weir B, Cockerham C. Estimating F-statistics for the analysis of population structure. Evolution. 1984; 38: 1358–70. doi: 10.2307/2408641 36. Raymond M, Rousset F. An exact test for population differentiation. Evolution. 1995; 49: 1280–3. doi: 10.2307/2410454 37. Chapuis M, Estoup A. Microsatellite null alleles and estimation of population differentiation. Mol Biol Evol. 2007; 24: 621–31. doi: 10.1093/molbev/msl191 PMID: 17150975 38. Piry S, Luikart G, Cornuet J. BOTTLENECK: a computer program for detecting recent reductions in the effective size using allele frequency data. J Hered. 1999; 90: 502–3. doi: 10.1093/jhered/90.4.502 39. Anderson M. A new method for non-parametric multivariate analysis of variance. Aust Ecol. 2001; 26: 32–46. 40. Bandelt H, Forster P, Rohl A. Median-joining networks for inferring intraspecific phylogenies. Mol Biol Evol. 1999; 16:37–48. PMID: 10331250 41. Weir B. Genetic Data Analysis II. Methods for discrete population genetic data. Sinauer Associates Sunderland MA. 1996 doi: 10.1111/j.1365-2052.2009.01970.x 42. Pritchard J, Matthew S, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000; 155: 945–59. PMID: 10835412 43. Falush D, Stephens M, Pritchard J. Inference of population structure using multilocus genotype data: dominant markers and null alleles. Mol Ecol Notes. 2007; 7: 574–8. doi: 10.1111/j.1471- 8286.2007. 01758.x PMID: 18784791 44. Earl D, von Holdt B. STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method. Cons Gen Res. 2012; 4: 359–61. doi: 10.1007/s12686- 011-9548-7 45. Kopelman NM, Mayzel J, Jakobsson M, Rosenberg NA, Mayrose I. CLUMPAK: a program for identifying clustering modes and packaging population structure inferences across K. Mol Ecol Resour. 2015; 15: 1179–91. doi: 10.1111/1755-0998.12387 PMID: 25684545 Depth-Related Genetic Patterns in Eunicella singularis PLOS ONE | DOI:10.1371/journal.pone.0160678 August 4, 2016 18 / 20 46. Jombart T, Devillard S, Balloux F. Discriminant analysis of principal components: a new method for the analysis of genetically structured populations. BMC Genet. 2010; 11: 94. doi: 10.1186/1471-2156-11- 94 PMID: 20950446 47. Benjamini Y Hochberg Y. Controlling the False Discovery Rate: a practical and powerful approach to multiple testing. J R Stat Soc Series B. 2015; 57: 289–300. doi: 10.2307/2346101 48. Piry S, Alapetite A, Cornuet J-M, Paetkau D, Baudouin L, Estoup A. GENECLASS2: a software for genetic assignment and first-generation migrant detection. J Hered. 2004; 95: 536–9. doi: 10.1093/ jhered/esh074 PMID: 15475402 49. Rannala B, Mountain JL. Detecting immigration by using multilocus genotypes. Proc Natl Acad Sci. 1997; 94: 9197–201. doi: 10.1073/pnas.94.17.9197 PMID: 9256459 50. Cornuet JM, Piry S, Luikart G, Estoup A, Solignac M. New methods employing multilocus genotypes to select or exclude populations as origins of individuals. Genetics. 1999; 153: 1989–2000. PMID: 10581301 51. Posada D. jModelTest: phylogenetic model averaging. Mol Biol Evol. 2008; 25: 1253–6. doi: 10.1093/ molbev/msn083 PMID: 18397919 52. Wörheide G, Hooper J, Degnan B. Phylogeography of western Pacific Leucetta chagosensis (Porifera: Calcarea) from ribosomal DNA sequences: implications for population history and conservation of the Great Barrier Reef World Heritage Area (Australia). Mol Ecol. 2002; 11:1753–68. PMID: 12207725 53. Gori A, Rossi S, Linares C, Berganzo E, Orejas C, Dale MR, et al. Size and spatial structure in deep versus shallow populations of the Mediterranean gorgonian Eunicella singularis (Cap de Creus, northwestern Mediterranean Sea). Mar Biol. 2011; 158: 1721–32. doi: 10.1007/s00227- 011-1686-7 54. Pey A, Catanéo J, Forcioli D, Merle P-L, Furla P. Thermal threshold and sensitivity of the only symbiotic Mediterranean gorgonian Eunicella singularis by morphometric and genotypic analyses. C R Biol. 2013; 336: 331–41. doi: 10.1016/j.crvi.2013.06.008 PMID: 23932253 55. Pivotto ID, Nerini D, Masmoudi M, Kara H, Chaoui L, Aurelle D. Highly contrasted responses of Mediterranean octocorals to climate change along a depth gradient R. Soc. open sci. 2015; 2: 140493. doi: 10. 1098/rsos.140493 56. Rossi S, Gili J-M. Reproductive features and gonadal development cycle of the soft bottomgravel gorgonian Leptogorgia sarmentosa (Esper 1791) in the NW Mediterranean Sea. Invertebr Reprod Dev. 2009; 53: 175–90. doi: 10.1080/07924259.2009.9652304 57. Rubio A, Barnier B, Jordá G, Espino M, Marsaleix P. Origin and dynamics of mesoscale eddies in the Catalan Sea (NW Mediterranean): insight from a numerical model study. J Geophys Res Ocean. 2009; 114(C6). 58. Théodor J. Contribution a l’étude des Gorgones (VII): Ecologie et comportement de la planula. Vie Milieu. 1967; 18: 291–301. 59. Weinberg S, Weinberg F. The life cycle of a gorgonian: Eunicella singularis (Esper, 1974). Bijdr. Dierk. 1979; 48: 127–40. 60. Padrón M. Evaluation of conservation efficiency for gorgonian species at a regional scale based on an existing Marine Protected Area network and modeling scenarios accounting for hydrodynamics. Université Pierre et Marie Curie, Universitàdi Bologna. 2015 61. Forcioli D, Merle P, Caligara C, Ciosi M, Muti C, Francour P, et al. Symbiont diversity is not involved in depth acclimation in the Mediterranean sea whip Eunicella singularis. Mar Ecol Prog Ser. 2011; 439: 57–71. doi: 10.3354/meps09314 62. Baums IB. A restoration genetics guide for coral reef conservation. Mol Ecol. 2008; 17: 2796–2811 doi: 10.1111/j.1365-294X.2008.03787.x PMID: 18482262 63. Holland LP, Stevens JR. Assessing the genetic connectivity of two octocoral species in the Northeast Atlantic. Natural England, Commissioned Report. 2014 64. Barbará T, Palma-Silva C, Paggi GM, Bered F, Fay MF, Lexer C. Cross-species transfer of nuclear microsatellite markers: Potential and limitations. Mol Ecol. 2007; 16: 3759–67. doi: 10.1111/j.1365- 294X.2007.03439.x PMID: 17850543 65. Quattrini AM, Baums IB, Shank TM, Morrison C, Cordes EE. Testing the depth-differentiation hypothesis in a deepwater octocoral. Proc R Soc B Biol Sci. 2015; 282: 20150008. doi: 10.1098/rspb.2015. 0008 66. Rouchy JM, Caruso A. The Messinian salinity crisis in the Mediterranean basin: a reassessment of the data and an integrated scenario. Sediment Geol. 2006; 188: 35–67. doi: 10.1016/j.sedgeo.2006.02. 005 67. Hinderstein LM, Marr JCA, Martinez FA, Dowgiallo MJ, Puglise KA, Pyle RL, et al. Theme section on ‘Mesophotic Coral Ecosystems: Characterization, Ecology, and Management’. Coral Reefs. 2010; 29: 247–51. doi: 10.1007/s00338-010-0614-5 Depth-Related Genetic Patterns in Eunicella singularis PLOS ONE | DOI:10.1371/journal.pone.0160678 August 4, 2016 19 / 20 68. Holstein DM, Paris CB, Vaz AC, Smith TB. Modeling vertical coral connectivity and mesophotic refugia. Coral Reefs 2015 doi: 10.1007/s00338-015-1339-2 69. Garrabou J, Coma R, Bensoussan N, Bally M, Chevaldonné P, Cigliano M. et al. 2009 Mass mortality in Northwestern Mediterranean rocky benthic communities: effects of the 2003 heat wave. Glob Change Biol. 2009; 15: 1090–1103. doi: 10.1111/j.1365-2486.2008.01823.x 70. Baird NA, Etter PD, Atwood TS, Currey MC, Shiver AL, Lewis ZA, et al. Rapid SNP discovery and genetic mapping using sequenced RAD markers. PLoS One 2008; 3: e3376. doi: 10.1371/journal. pone.0003376 PMID: 18852878 Depth-Related Genetic Patterns in Eunicella singularis PLOS ONE | DOI:10.1371/journal.pone.0160678 August 4, 2016 20 / 20