- 1Aix Marseille Univ, Université de Toulon, Centre National de la Recherche Scientifique, Institut de Recherche pour le Développement, Institut Méditerranéen d’Océanologie, Marseille, France
- 2Institut de Systématique, Evolution, Biodiversité (ISYEB), Muséum National d’Histoire Naturelle, CNRS, Sorbonne Université, École Pratique des Hautes Études, Paris, France
- 3IRL 3614, Evolutionary Biology and Ecology of Algae, CNRS, Sorbonne Université, Université Catholique, Université Australe du Chili, Roscoff, France
- 4Institut Français de Recherche pour l’Exploitation de la Mer, Zone Portuaire de Brégaillon, La Seyne-sur-mer, France
Dispersal is a central process that affects population growth, gene flow, and ultimately species persistence. Here we investigate the extent to which gene flow occurs between fragmented populations of the deep-water brown algae Ericaria zosteroides (Turner) Greville (Sargassaceae, Fucales). These investigations were performed at different spatial scales from the bay of Marseille (western Provence) to Corsica. As dispersal of zygotes is shown to be limited over distances beyond a few meters, we used a multidisciplinary approach, based on Lagrangian modeling and population genomics to test the hypothesis that drifting of fertile parts of thallus (eggs on fertile branches), mediated by ocean currents, enable occasional gene flow between populations. Therefore we assessed the respective contribution of oceanographic connectivity, geographical isolation, and seawater temperatures to the genetic structure of this species. The genetic structure was assessed using 10,755 neutral SNPs and 12 outlier SNPs genotyped by dd-RAD sequencing in 261 individuals of E. zosteroides. We find that oceanographic connectivity is the best predictor of genetic structure, while differentiation in outlier SNPs can be explained by the depth of populations, as emphasized by the minimum seawater temperature predictor. However, further investigations will be necessary for clarifying how depth drives adaptive genetic differentiation in E. zosteroides. Our analyses revealed that local hydrodynamic conditions are correlated with the very high divergence of one population in the Bay of Marseille. Overall, the levels of gene flow mediated by drifting were certainly not sufficient to counteract differentiation by local genetic drift, but enough to allow colonization several kilometers away. This study stresses the need to consider secondary dispersal mechanisms of presumed low dispersal marine species to improve inference of population connectivity.
Introduction
Genetic diversity is a prerequisite for the adaptive evolution of species. Accordingly, surveys of genetic diversity of rare and vulnerable species should be on the agenda of management and conservation frameworks (Hoban et al., 2020). This involves, among other requirements, the need to estimate the levels of genetic connectivity of species and their effect on the genetic structure of natural populations (Slatkin, 1985, 1987). Here, we will use the following definition of genetic connectivity: the degree to which the dispersal of individuals or propagules affect gene flow, with repercussions on evolutionary processes within populations, while demographic connectivity refers to the contribution of dispersal to the rate of population growth (Lowe and Allendorf, 2010). Genetic connectivity can counteract the potential loss of genetic diversity due to local inbreeding (i.e., inbreeding connectivity, Lowe and Allendorf, 2010), genetic drift, and the fixation of deleterious mutation (Keller and Waller, 2002). Demographic connectivity can maintain the stability of populations structured across space by increasing the range of distribution and enhancing the resilience of populations following disturbances (Gotelli, 1991; Lowe and Allendorf, 2010). Genetic connectivity can also counteract the evolution of local adaptation (Lenormand, 2002). The interplay between genetic connectivity, demography and local selection is then at the heart of the evolution of populations. In the marine environment, many organisms disperse through the release of many small propagules, which consequently require the use of indirect methods such as genetic inference and/or Lagrangian models for assessing dispersal. Combining these indirect methods in a multidisciplinary approach was seen as an effective way to create the most complete picture of population connectivity, notably in marine species with low dispersal potential for which connectivity may be driven by unsuspected occasional long-distance dispersal events (for review see Kinlan and Gaines, 2003).
Next−generation sequencing (NGS) has turned the study of genetic connectivity from the use of tens to thousands of markers, going up to whole individual genomes (Gagnaire, 2020). The use of this high number of markers allows more precise analysis of genetic differentiation, and increased power to detect genetic differences among populations (Waples, 1998; Gagnaire et al., 2015). Information about the spatial distribution of adaptive genetic diversity can also be considered in management and conservation strategies (Xuereb et al., 2020). Combining population-genomic inferences and Lagrangian dispersal models has helped to elucidate, with unprecedented clarity, how environmental factors and geography drove the genetic structure of marine species (e.g., Benestan et al., 2016; Dalongeville et al., 2017; Paterno et al., 2017; Xuereb et al., 2018a,b; Bracco et al., 2019; Coscia et al., 2020). To our knowledge, such integrative frameworks at genome-scale have rarely been applied to brown algal species. The scale of the study is a critical point in such investigations of connectivity. According to the characteristics of the studied species, we can define a fine scale (from a few centimeters and a dozen meters within a site) that can be used to study population structure and delimitation (demes) (e.g., Ledoux et al., 2010), a local scale (from hundreds of meters to a few kilometers) where the genetic structure can be influenced by local hydrodynamic currents (e.g., Palumbi, 2003) and by the stochasticity of recruitment, and a regional or continental scale where permanent oceanic currents are the main drivers of connectivity (e.g., White et al., 2010). At a local scale, the high variation in hydrodynamic conditions can, in some cases, be caused by different wind direction and intensity while permanent oceanic currents over large regions are mainly driven by salinity, temperature and coriolis effects.
In the Mediterranean Sea, Fucales (Cystoseira C. Agardh, Ericaria Stackhouse, Fucus Linnaeus, Gongolaria Boehmer, and Sargassum C. Agardh) (see Molinari Novoa and Guiry, 2020 for an update of the taxonomy of the Cystoseira sensu lato) form marine forests and structure rocky reef assemblages. Most of the time, populations are patchy, from a few individuals to larger populations dwelling on rocks or boulders surrounded by detrital soft bottoms, coralligenous and Posidonia oceanica meadows (e.g., Ballesteros, 2006; Hereu et al., 2008; Boudouresque et al., 2017). Very few studies have focused on the connectivity of Mediterranean Fucales at different scales, and all concerned shallow water species (i.e., Susini et al., 2007; Robvieux, 2013; Thibaut et al., 2016; Buonomo et al., 2017; Medrano et al., 2020). Previous studies have reported that populations of Ericaria amentacea (C. Agardh) Molinari and Guiry (synonym: Cystoseira amentacea) were genetically differentiated at local and regional spatial scales, from 100 m to tens of kilometers (Susini et al., 2007; Thibaut et al., 2016; Buonomo et al., 2017). Similar levels of genetic differentiation were also observed within the Ericaria selaginoides (Linnaeus) Molinari and Guiry complex (formerly called Cystoseira tamariscifolia complex) (Bermejo et al., 2018). Such high levels of genetic differentiation could have been expected given the low dispersal ability of these species, which is directly related to the short interval between external fertilization (i.e., mating process between motile male flagellate gametes and large female, non-motile gametes) and fixation of fertilized zygotes to the substrate. Large and non-buoyant zygotes of several species of Fucales (Cystoseira spp., Ericaria spp., and Gongolaria spp.) sink and adhere to the substrate a few hours after fertilization (Guern, 1962; Falace et al., 2018). This mating pattern is expected to only allow restricted dispersal events (less than 40 cm–10 m in Ericaria species; Mangialajo et al., 2012; Capdevila et al., 2018). Importantly, the local genetic structure of Ericaria amentacea in the Bay of Marseille was better explained by oceanographic distance rather than spatial distance (Thibaut et al., 2016). Furthermore, and non-exclusively, the high population genetic differentiation in Ericaria and Gongolaria species is probably driven by the balance between restricted gene flow and strong genetic drift, as reported in some Mediterranean animal species (e.g., Corallium rubrum, Ledoux et al., 2010). Up to now, connectivity studies in Ericaria and Gongolaria have involved relatively small numbers of genetic markers (e.g., about ten microsatellites).
The connectivity of the deep-water species Ericaria zosteroides (C. Agardh) Molinari and Guiry (basionym: Cystoseira zosteroides C. Agardh) dwelling down in the lower zone of the circalittoral (down to 80 m depth) has never been investigated. Ericaria zosteroides is an endemic perennial and habitat-forming Mediterranean species and one of the most representative of the brown seaweeds of deep-water assemblages (Giaccone, 1973; Ballesteros, 1990; Ballesteros et al., 2009; Boudouresque et al., 2017). In addition to the uprooting of entire individuals by storms, the season’s shedding of branches constitutes a vector for long-distance dispersal of eggs. During the last decades, anthropogenic disturbances (e.g., uprooting by fishing gear, increasing water turbidity, Thibaut et al., 2005, 2015), episodic climatic events and storms (Capdevila et al., 2015) have drastically fragmented the populations of this species, with possible repercussions on the connectivity among populations. Furthermore, rising seawater temperatures due to climate change could have dramatic repercussions on the population recovery of this vulnerable species (Capdevila et al., 2019). As reported for other Ericaria species from the upper sublittoral zone, E. zosteroides is also able to colonize artificial habitats (Bonhomme et al., 2014, Thibaut pers. obs.) at distances greater than the dispersal scale of its zygotes (Capdevila et al., 2018). This can indicate that occasional long-distance dispersal events have a great influence on the connectivity of E. zosteroides.
In the present study, we sampled ten sites situated less than 1 km to several hundreds of km apart to investigate the pattern of spatial structure at different scales. At the local scale (a few kilometers), we aimed at assessing the potential for distance dispersal connecting the different sites, perhaps driven by ocean currents. We therefore conducted inference of genetic structure with RAD sequencing (dd-RADseq; Peterson et al., 2012) in a complementary way with Lagrangian modeling to unravel how gene flow is mediated by drifting of fertile parts of the thallus. This multidisciplinary approach enabled us to answer the following questions: (i) Is E. zosteroides spatially structured at different spatial scales into genetically distinct populations? (ii) What is the respective influence of demography, oceanographic currents, and geographical distance on the genetic structure? (iii) Can we identify the sources of recruitment?
Materials and Methods
Sampling and Demographic Structure of Populations
Two hundred and sixty-one individuals of Ericaria zosteroides were sampled at 10 sites in 2017 and 2018 by scuba diving, and a small fragment of the thallus (2–3 cm) was preserved in silica gel until DNA extraction. Populations were sampled in the north-western Mediterranean Sea from the Bay of Marseille (western Provence) to Corsica (Table 1 and Figure 1). We included in the study a population settled on artificial substrates, the subsea-pipelines of the company Alteo (ALT) (Bonhomme et al., 2014). Furthermore, we focused on a site, L’Armoire, which is an isolated rocky bank from 22 to 40 m depth surrounded by sand where a continuous population of E. zosteroides is thriving over the whole of all the depth range of the bank (ARM). This rare spatial configuration offered a great opportunity to study vertical genetic connectivity. For this purpose, we separately sampled individuals at 22–24 m depth (ARM20) and 38–40 m depth (ARM40). The density of individuals of E. zosteroides was measured using random quadrats (1 m2). A total of 214 quadrats were performed during the study. The length of the major axis of each individual is measured in situ with a plastic ruler and an accuracy of ±5 mm. Individuals of less than 1 cm are considered as recruits of the year. From this length distribution, we defined cohorts according to discrete classes of 1 cm (see Hereu et al., 2008 for details); in total 289 individuals were measured. Because of the difficulties of working at these depths, we could not measure the length of the axis at ARM40 and LC, and we could not measure the demographic parameters at SC. Densities were compared with a non-parametric Kruskal-Wallis test followed by a post-hoc pairwise Wilcoxon test. The distribution of the length was compared pairwise with a Kolmogorov-Smirnov test. The mean size of the main axis was compared using ANOVA followed by a pairwise Tuckey’s test. In this study we define three scales: the fine scale for individuals within the same site placed up to 20 m apart, the local scale from Marseille to La Ciotat up to 29 km, and the regional scale from western Provence to Corsica, with a minimum distance by sea of 216 km between these two areas.
Table 1. Genetic diversity within populations of Ericaria zosteroides based on 10,755 neutral SNPs (n, number of individuals successfully genotyped; Ho, observed heterozygosity; He, expected heterozygosity and SE among loci; FIS, average inbreeding coefficient).
Figure 1. Sampling map of Ericaria zosteroides and genetic ancestry estimates derived from sNMF at K = 10 clusters using the set of 10,755 putatively neutral SNPs. The pie plots indicate ancestry averaged over individuals per site, with different colors indicating admixture. The dotted box indicates that ARM20 and ARM40 populations have the same geographical position but differ by depth. Genetic ancestry is also reported at the individual level. Populations are labeled according to the nomenclature defined in Table 1.
DNA Extraction and SNP Discovery
Total genomic DNA was isolated with the Nucleospin® 96 plant kit II (Macherey-Nagel, Düren, Germany), according to the manufacturer’s protocol. Then, SNP discovery was performed using double digest RAD sequencing libraries (dd-RADseq; Peterson et al., 2012). The protocol used to build ddRAD-seq libraries for sequencing is identical to that previously described in Reynes et al. (2020). Two hundred and seventy-one individuals of Ericaria zosteroides were randomly distributed across two sequencing libraries, including 10 replicates (one per sampling site), and 100 ng of genomic DNA per individual was double digested using the PstI and HhaI enzymes (NEB). Finally, dd-RADseq libraries were sequenced with the paired-end method (2 × 150 bp) on an Illumina Hiseq 4000 platform (Génome Québec Innovation Centre, McGill Univ., Montreal, Canada).
Paired-end raw reads were checked for control quality using FASTQC v.0.11.7 (Andrews, 2010) and trimmed to 137 bp after the removal of adapters sequencing by Trimmomatic (Bolger et al., 2014). High-quality reads were demultiplexed using Stacks’s process_radtags (Catchen et al., 2011). Additionally, one individual genome of E. zosteroides was paired-end sequenced in short reads (2 x 150 bp) on an Illumina NovaSeq (Genewiz, Inc., South Plainfield, NJ). Then, de novo genome assembly was performed using SOAPdenovo v.2.41, with a minimum contig length of 1000 bp. The resulting draft assembly of the E. zosteroides nuclear genome was used as a reference genome to map paired-end reads of the RAD-sequencing experiment using the BWA-mem algorithm in BWA v.0.7.17 (Li and Durbin, 2009). Aligned reads were then filtered using Sambamba v.0.6.6 (Tarasov et al., 2015) to discard all possible multi-mapped reads to the genome. RAD loci were built from sets of alignments using the reference mode of the Stacks v2.4 pipeline (Rochette et al., 2019), then SNP discovery was done with the Bayesian genotype caller model (BGC; Maruki and Lynch, 2017) implemented in the Stacks v2.4. pipeline.
SNP Filtering
Initially, we ran the population module of Stacks v2.4 requiring a locus to be genotyped at a frequency greater than 80% in at least 8 out of 10 sampling sites, while deleting SNPs with more than 50% of heterozygosity. Then, stringent quality filters were applied to the dataset (see Supplementary Table 1) using PlinkQC, the R version of PLINK v1.9 (Chang et al., 2015). This procedure consisted of i) discarding individuals outlying missing genotype frequency >25% and those outlying three SD away from the mean rate of heterozygosity ii) dropping SNPs with missingness rate >5% and SNPs with minor allele frequency (MAF) <1%. Following this step, SNP error rate between replicates ranged from 0.009 (OURS) to 0.03 (SC) with an average rate of 0.01. Also, individual relatedness was checked using the KING-robust kinship estimator (Manichaikul et al., 2010), implemented in the Bigsnpr v1.3.0 R package (Privé et al., 2018). We applied a KING-relatedness cutoff of 0.354 to discard duplicate samples and up to a third degree of relationship. Finally, SNP clumping was performed using plink to remove associated SNPs (r2 > 0.3) using the MAF as a statistic of importance. At the end of this filtering process, the dataset comprised 11,067 SNPs among 229 individuals, with a genotyping rate of 98%.
Species Delimitation
Ericaria zosteroides is present in sympatry with other species of Cystoseira sensu lato, mainly Gongolaria montagnei var. compressa and E. funkii (e.g., Hereu et al., 2008; Ballesteros et al., 2009). As outlined by Pante et al. (2015), evolutionary units (e.g., species and populations) must be properly delimited to conduct a robust inference of population connectivity (see also Boero, 2010; Boudouresque et al., 2020). Hence, we tested for the existence of potential cryptic species in E. zosteroides, and subsequently validated our morphological inspection by performing an analysis of species delimitation (see Supplementary Material 1).
Outlier Loci
We searched for loci potentially influenced by selection, to build a neutral dataset used to infer population structure and connectivity. The putatively non-neutral (selected) loci were analyzed separately to test for adaptation to local conditions (see Supplementary Material 2). We applied three tests to minimize the number of false positives loci assumed to be under selection. First, we used the FST method described in Martins et al. (2016) and implemented in the LEA R package (Frichot et al., 2014; Frichot and François, 2015). This method starts by defining the optimum number of populations with ancestry estimates (Frichot et al., 2014). FST values are then converted to z-scores and significant outliers are detected after controlling for a False Discovery Rate (FDR) of 0.05. Secondly, we used the Pcadapt v.4.3.1 R package (Luu et al., 2017). Pcadapt analyses the population structure by PCA, then it estimates to what extent each SNP is related to the first K principal components to identify outlier loci. Here, outlier tests were performed from the first six principal components, according to the percentage of explained variance for each PC. P-values were then computed by making a Gaussian approximation for each PC and by estimating the standard deviation of the null distribution. Finally, P-values were corrected at 0.05 FDR. The third approach used here, implemented in Bayescan v. 2.1, estimates directly the probability that each locus is subject to selection using a Bayesian method (Foll and Gaggiotti, 2008). The test was run with default settings (5,000 iterations, 10 thinning intervals, 20 pilot runs—5,000 iterations each, 10 prior odds). SNPs were identified as significant outliers after FDR correction of 0.10.
For the following analyses of genetic variation, we selected the set of putatively neutral SNPs by excluding those identified as an outlier by at least two out of the three methods. The set of putatively outlier SNPs used for the following analyses of local adaptation (see Supplementary Material 2) was built by keeping all SNPs identified as an outlier by the three methods.
Population Genetic Variation
Spatial population structure was assessed with putatively neutral SNPs at fine (one site, two depths), local (six sites) as well as regional scales (ten sites). Two complementary approaches were used. First, a principal component analysis (PCA) was performed by singular value decomposition (SVD) of the centered genotype matrix, as implemented in Bigsnpr v. 1.3.0 and Bigstatsr v. 1.2.3 R packages (Privé et al., 2018). The missing genotypes were imputed according to the mean allele frequencies using the function snp_fastImputeSimple before performing the PCA. Second, the function sNMF of the LEA R package (Frichot et al., 2014; Frichot and François, 2015) was used to estimate individual ancestry proportions. SNMF was run with 5,000 iterations and 10 repetitions for K value from 1 to 15. The best K value was defined using the cross-entropy criterion (see Supplementary Figure 1). Genetic differentiation was measured by calculating pairwise FST values (Weir and Cockerham, 1984) between populations using the Genepop v. 4.7.5 R package (Rousset, 2008) and performing an exact test of genic differentiation (Raymond and Rousset, 1995) using the Markov chain method with default parameters, then a combination of all tests across loci (Fisher’s method) was performed for each population pair.
Genetic diversity within populations was estimated using Genepop by estimating the average observed heterozygosity (Ho), expected heterozygosity (He), and FIS (Weir and Cockerham, 1984). We tested for differences in He between populations with a pairwise Wilcoxon test adjusted for multiple comparisons (Bonferroni adjustment). Positive and negative FIS values were observed, so we tested departure from panmixia across all loci and populations using the Global Hardy-Weinberg test [Score (U) test] (Rousset and Raymond, 1995), implemented in Genepop. The test was performed twice: once by considering the alternative hypothesis of heterozygote excess and a second time with heterozygote deficiency. The P-value of each test was approximated by a Markov chain (MCMC) algorithm with the following settings: dememorization: 10,000, batch: 100, iterations per batch: 5,000. Then, a multiple score test (Rousset and Raymond, 1995) was performed to obtain a global P-value per population.
Individual Assignment Analysis
Individual assignment tests were used to test for contemporary gene flow and admixture between shallow (ARM20) and deep (ARM40) populations of Armoire, as well as to investigate the potential source population(s) for the colonization of ALT. Assignment tests were conducted using the Assigner v. 0.5.7 R package (Gosselin et al., 2016) and the implemented gsi_sim function (Anderson et al., 2008). We used the Training, Holdout and Leave-one-out (THL) method to avoid strong grading bias related to the fact that the selection of informative loci and population assignment is performed from the same set of individuals (Anderson, 2010). Specifically, individuals were distributed into two sets: a training set (including 50% of individuals per site), a holdout-set (including the other half of individuals per site). These analyses were standardized by considering the minimal sample size (i.e., 16 individuals). The selection of informative genetic markers was performed from the training dataset, thus identifying 14 sets of neutral SNPs ranked based on decreasing overall FST values (50, 100, 200, 500, 1,000, 2,000, 3,000, 4,000, 5,000, 6,000, 7,000, 8,000, 9,000, and all SNPs). Finally, the individual assignment was performed using separately the 14 sets of SNPs for the hold-out set of individuals, and confidence intervals were obtained by bootstrapping individuals (i.e., 16 individuals per population were randomly chosen 10 times).
Hydrodynamic Modeling
Oceanographic connectivity was assessed at local scale, from the Bay of Marseille to La Ciotat (western Provence), which is an area where the hydrodynamic processes showed high spatio-temporal heterogeneity (Pradal and Millet, 2006). The MARS3D model (3D hydrodynamic Model for Applications at Regional Scale, IFREMER, Lazure and Dumas, 2008) in its RHOMA configuration (Pairaud et al., 2011; Fraysse et al., 2013) was used to generate three-dimensional oceanic circulation patterns. For our study, the model used 120 x 252 horizontal grid points at 400 m resolution and 30 sigma layers, with a corresponding time step of 30 s, took into account the forcing by the North-Western Mediterranean general circulation, using a nesting strategy with the large scale MARS3D-MENOR configuration (1.2 km resolution) (Nicolle et al., 2009). We expected that oceanographic connectivity was highly variable according to the year of simulation, as such Lagrangian modeling was performed over seven years from the ocean circulation patterns of the years 2008 to 2017. Hydrodynamic model were forced hourly by different atmospheric models (MM5, Arome, and WRF models). To ensure the robustness of our results, we focused on this long period considering a wide variety of typical atmospheric forcing and oceanic features of the region of Marseille in the calculation of the mean pattern of oceanographic connectivity. Notably, strong NW and SE winds (Pradal and Millet, 2006), the intrusion of the Northern Mediterranean Current (Barrier et al., 2016; Ross et al., 2016; Millet et al., 2018), and diluted waters from the Rhone River plume (Fraysse et al., 2014).
Lagrangian Modeling
We first tested if dispersal of fragments of fertile thalli allows genetic connectivity between spatially distant populations, which requires that fertile fragments be transported to another population where they produce viable gametes which contribute to fertilization and reproduction in the target population. Detached secondary branches with receptacles are commonly observed around the E. zosteroides populations. Lagrangian simulations were performed with Ichthyop v. 3 targeting the breeding season of E. zosteroides, known to take place from May to September. Simulations were performed on the high-performance computing (HPC) cluster of the OSU Pythéas Institute (Marseille) and were repeated every day for the first 28 days of May, June, July, and August, amounting to a total of 784 release events per site over seven years. ICHTHYOP was parameterized to constraint dispersal of non-buoyant particles around the seafloor. Specifically, 1,000 particles were released each day from each of six locations within a radius of 400 m at the depth (±2 m) of the study sites. The computational time step (dt) was set to 400 s and particle position was recorded in the NetCDF output files every 1.33 h. Pairwise oceanographic connectivity probabilities were defined as Pij = Ni–j/Ni where Ni is the total number of released particles from site i and Ni–j the total number of released particles from site i arriving in site j after 30 days of transport time; these probabilities were computed in MATLAB. To take into account the life span of the species, particles recorded at more than 100 m depth were not considered in the computation of Pij. Finally, the probabilities of oceanographic connectivity between all pairs of sites were averaged over 7 years of simulation to obtain the mean connectivity matrix and the mean time of transport.
Environmental Predictors of Genetic Structure
To test the effects of environmental parameters on the genetic structure of E. zosteroides, we selected four sets of explanatory variables: oceanographic connectivity, geographical distance, population density and seawater temperature. Geographical distances were computed as distance-based Moran’s eigenvector maps or dbMEM (Dray et al., 2006), which accounts for the effect of spatial variation among sites. To do this, pairwise Euclidean distances between sites were calculated using the codep v. 0.9.1 R package, then the distance matrix was used as input data for computing dbMEM eigenvectors with the adespatial v.0.3.8 R package. We kept the default settings for truncations. Oceanographic connectivity was computed as Asymmetric Eigenvector Maps or AEM (Blanchet et al., 2008), which turned out to be extremely useful for modeling asymmetric processes, such as ocean currents, and testing their effects on genetic structure (e.g., Benestan et al., 2016; Xuereb et al., 2018a; Coscia et al., 2020). For this purpose, we transformed the mean oceanographic connectivity matrix (Figure 7) from Lagrangian simulations into a nodes-to-edges matrix, which records the presence/absence of connectivity link between nodes. Specifically, fifteen connectivity links (i.e., edges) were reported among the six study sites considered (i.e., nodes). Edges were weighted by the probabilities of dispersal, and this weighted matrix was then used as input data to compute AEM predictors with the adespatial v.0.3.8 R package. Seawater temperature predictors were computed by extracting the temperature variables of the MARS3D model. Specifically, we considered minimum (Min_T) and maximum (Max_T) seawater temperatures, the annual mean (Mean_T), and variance (Var_T). These parameters were extracted from the site’s coordinates accounting for the depth of populations and were averaged over seven years of simulations. For the density predictor, we used the mean value of populations reported in our study (Table 2).
Table 2. Demographic structure of Ericaria zosteroides with the density of populations (n, number of quadrats per site; average, average density, and SE), the size of individuals (n, number of individuals; average, average with SE; and Max, maximum size of the main axis), and age: estimate range of the age of individuals according to a growth rate estimation for the axis of 0.5 to 0.9 cm.year–1.
Distance Based RDA: Linking Environment to Genetic Structure
At a local scale (Marseille-La Ciotat), partial and global distance based RDA (Legendre and Anderson, 1999) were used to investigate, respectively, the separate and joint effects of oceanographic connectivity (AEM1; AEM-2; AEM-3; AEM-4; AEM-5), geographical distance (dbMEM-1), population demography (density), and seawater temperature (Min-T, Max-T, Mean-T, Var-T) on the response variable: individual coordinates of principal coordinate analysis (PCoA) based on the pairwise Euclidean genetic distance matrix. These analyses were carried out using the Vegan v. 2.5.6 R package (Oksanen et al., 2017) and repeated twice with the set of 10,755 neutral and 12 outlier SNPs. We considered all independent predictors, with Variance Inflation Factor (VIF) below 10 to avoid multicollinearity effects (pairwise correlations of the initial dataset of 11 environmental predictors were reported in Supplementary Table 2). Following the VIF threshold, Min-T and Var-T predictors were excluded in partial and global distance based RDA while the density predictor was only excluded in the whole model. Then, we applied stepwise forward selection to select the best predictors using the ordiR2step R function. The ordiR2step performs model choice solely on the adjusted coefficient of determination (R2adj) and selects the model with the highest one. The statistical significance of partial and global Distance Based RDA and their respective predictors were determined by analysis of variance (ANOVA) with 1000 permutations. Finally, R2adj and the p-value for the F-tests were used to evaluate the validity of the models.
Results
Demographic Structure
Populations of Ericaria zosteroides are spatially isolated in the studied zone, without a continuum between them. The density of populations was highly variable and significantly different between sites [KW test, H = 125.25, df (7), p < 0.001] (Figure 2). The mean population density ranged from 1.23 ± 0.25 ind.m–2 recorded in the Bay of Marseille at RP, to 11 ± 0.62 ind.m–2 at ARM20 (Table 2). The size distribution of individuals was statistically only similar in 8 out of the 21 pairwise comparisons (Kolmogorov-Smirnov pairwise test, p > 0.05) (see Supplementary Table 3). Only one cohort was observed in five populations (TIB, VEY, OUR, ALT, ARM20), two cohorts in the RP population, and three cohorts in CAV (see Supplementary Figure 2). The mean size of the main axis ranged from 3.5 ± 0.2 cm to 7.4 ± 0.5 cm in ALT and OUR, respectively. The mean size differed among sites [F(6, 282) = 11.91, p < 0.001]. Individuals at ALT were significantly smaller than all the other populations studied, except for the comparison with TIB (Tukey test) (Figure 2). According to the estimated growth rate, which ranged from 0.5 to 0.9 cm.year–1 (Ballesteros et al., 2009; Navarro et al., 2011), the maximum size of the main axis seems to indicate that the oldest individuals reach 8–14 years on the pipeline (ALT), against 11 to 31 years in other populations (Table 2).
Figure 2. Demographic structure of populations of Ericaria zosteroides. (A) spatial distribution of the population density (ind.m– 2), all pairwise comparisons are significant (pairwise Wilcoxon test, p < 0.001), except for those indicated by p > 0.05. (B) Distribution of the height of the main axis of individuals, the significant difference between populations (Tukey test, *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001) are reported.
Bioinformatic and Genotyping
The draft genome of Ericaria zosteroides was sequenced (51.98 Gb of raw data), and de novo assembled into 123,092 contigs of a mean length of 2,872 bp. The ddRAD-sequencing experiment enabled us to generate and filter 974,562,924 paired-end high-quality reads among the 261 individuals of E. zosteroides, consisting of an average of 3,596,173 (±1,356,820 SE) high-quality reads per individual. On average, 75.50% (min: 27.98%; max: 91.34% per individual) of ddRADseq reads were successfully mapped to the E. zosteroides genome. Following the SNP calling process by Stacks v.2.4, mapped reads were assembled into 314,538 loci in which 721,374 SNPs were identified. The SNP filtering procedure (see Supplementary Table 1) reduced the initial set of 721,374 SNPs to 21,606 SNPs, while also reducing the number of individuals from 271 to 229. Finally, clumping on MAF (r2 > 0.3) allowed retention of 11,067 SNPs.
From the set of 11,067 SNPs, BayeScan detected 39 outlier SNPs, while pcadapt and sNMF identified 895 and 482 outliers, respectively. Twelve outliers were shared between the three outlier tests (Figure 3). Furthermore, 312 SNPs detected by at least two of the three methods were excluded from the set of 11,067 SNPs and those remaining (10,755 SNPs) constitute the neutral dataset.
Figure 3. Venn diagram denoting overlap between SNPs identified as an outlier with sNMF, Pcadapt, and BayeScan. The set of putative selected SNPs was built by considering the 12 SNPs shared among these three methods.
Species Delimitation
The results of the analysis of species delimitation are presented in Supplementary Material 1. All samples identified as E. zosteroides included in this analysis were grouped in a clade. The divergence among samples from the different populations analyzed here was much lower than the divergence among samples from different species.
Genetic Variation Within Populations
The genetic diversity (He) of most of the populations differed among sites (Pairwise Wilcoxon Test, Bonferroni adjusted: p < 0.001), except for 8 of the 45 pairwise comparisons (see Supplementary Table 4). The lowest levels of observed (Ho) and expected (He) heterozygosity were recorded in the population of RP (Ho = 0.07, He = 0.08) in the Bay of Marseille, while the highest ones were reported in Corsica at SC (Ho = 0.13, He = 0.12). The level of genetic diversity (He) on the pipelines (ALT) was similar to those reported in neighboring populations, such as CAV (pairwise Wilcoxon test, Bonferroni adjusted: p = 0.86). On the isolated Armoire reef bank, the shallowest (ARM20) and deepest (ARM40) sites did not show any significant difference in genetic diversity (see Supplementary Table 4), but the average values of He and Ho were low in this locality compared to other sites (Table 1). We observed a significant excess in heterozygotes at ARM20 and ARM40, while 5 of the 10 populations were characterized by significant deficits in heterozygotes (Table 1). The mean FIS ranged from slightly negative values at ARM20 (FIS = -0.02) and ARM40 (FIS = -0.02) to positive values, notably at RP (FIS = 0.04) and ALT (FIS = 0.04).
Genetic Structure and Assignment Tests
The overall estimated FST was 0.09, and it was slightly lower on a local scale within the Bay of Marseille (FST = 0.08). Populations were significantly genetically differentiated from each other (exact G test for genic differentiation, P < 0.001), and pairwise FST ranged from 0.05 to 0.14. The highest FST values were observed for the comparison of the RP population with other ones (Figure 4). The lowest genetic differentiation (FST = 0.05) was observed between CAV and TIB populations, as well as between shallow (ARM20) and deep (ARM40) populations at the Armoire site (Figure 4).
Figure 4. Heatmap of pairwise FST between populations of Ericaria zosteroides using the set of 10,755 SNPs putatively neutral SNPs. Pairwise genetic differentiation was reported as significant among populations (exact G test for genic differentiation, p < 0.001).
The principal component analysis (PCA), and sNMF clustering gave similar results and were congruent with the FST values. The first axis of PCA (PC1, 19% of the variance) allowed a distinction between ARM20/ARM40 individuals and other ones, while the second axis (PC2, 13% of the variance) discriminated the population of RP from other ones (Figure 5A). Additional PC axes showed the separation of the other populations (see Supplementary Figure 3). By performing the PCA only on individuals from the Bay of Marseille, ALT, and LC (Figure 5B), the first axis (PC1, 18% of the variance) differentiated individuals from RP from other sites, while the second axis (PC2, 17% of the variance) showed differentiation between three groups: one with individuals from ALT and LC, the second with individuals of CAV and TIB, the third with individuals of VEY. However, one individual of VEY was grouped with those of CAV and TIB.
Figure 5. Principal component analysis (Axes 1 and 2) computed using the set of 10,755 putative neutral SNPs (A) with all individuals (n = 229) and (B) restricted to individuals (n = 143) sampled from the Bay of Marseille, Cassis to La Ciotat.
According to the minimum cross-entropy criterion, K = 10 seemed to be the most informative number of clusters in sNMF (see Supplementary Figure 1). All individuals were clustered according to their populations of origin, except at ARM40 for which 17% of admixture was detected with ARM20. The admixture rate of ALT with LC was relatively low (6%) but higher than that observed with other populations (admixture range: 2–3%) (Figure 1).
The proportion of individuals assigned to their populations of origin reached 100% using the set of 10,755 neutral SNPs, except for ARM40 for which one individual was assigned to ARM20 (see Supplementary Figure 4). In contrast, no individual of ARM20 was assigned to ARM40. The success of the assignment increased with the number of SNPs considered (see Supplementary Figure 5). The overall assignment success ranged from 73% with the top-ranked 50 SNPs to 99% with the top-ranked 500 SNPs, and reached a plateau at 100% by using more than 500 SNPs.
Oceanographic Connectivity
According to the results of Lagrangian simulations, dispersal events among populations of E. zosteroides were scarce. The mean probability of oceanographic connectivity (P) was below 0.01 for most of the pairwise comparisons (Figure 6A). Across the Bay of Marseille, the simulations indicated that the population of RP lacked oceanographic connectivity with other populations. The highest probabilities of connection, still very low, were observed among the CAV, TIB, and VEY populations (P ranging from 0.04 to 0.06). Regarding the connectivity with the ALT population on pipelines, there was only a unidirectional connection from LC to ALT (Figure 7), with P = 0.02, and the corresponding mean time of transport (around 10 days; Figure 6B) was the lowest for the comparisons tested here.
Figure 6. Pairwise oceanographic connectivity from the source area to sink area. (A) The average probability of dispersal between populations with (B) the associated mean time of particle transport (in day). NA value of the diagonal matrices denotes that self-recruitment was not considered; another NA indicates the absence of dispersal between populations to compute transport time.
Figure 7. An example of a probability density map of the distribution of particles recorded at the local scale (Marseille–La Ciotat) over Lagrangian simulations of the year 2014. Probabilities are color-coded in each mesh cell (resolution of 200 meters), indicating no particle record (white) from low (blue) to high probability (yellow). For each panel (A–F), the black arrow indicates the source area considered among the six sites, materialized as a red dot.
Factors Explaining the Genetic Structure at a Local Scale (Marseille, ALT, and LC)
For the neutral genetic structure of E. zosteroides, oceanographic connectivity relating to AEM-1, AEM-2, and AEM-4 predictors, as well as geographical distance relating to the db-MEM-1 predictor, were selected by the stepwise forward selection (Table 3), and included in the global distance based RDA. The global distance based RDA was significant (P < 0.001) with an adjusted coefficient of determination (R2adj) of 0.11. When partitioning the respective effects of oceanographic connectivity, geographical distance, density and temperature in partial distance based RDA, all of them were significant predictors of neutral genetic structure, but oceanographic connectivity explained the greatest amount of variance with R2adj = 0.11 against R2adj = 0.03 for both geographical distance and density predictors, and R2adj = 0.05 for temperature predictors (Max-T and Mean-T). However, the temperature was not reported as a significant predictor of neutral genetic structure in the global Distance Based RDA (Table 3). The first axis of the global distance based RDA explained 30.3% of the total variance and it was mainly explained by oceanographic connectivity with the AEM-1 predictor (Figure 8). This axis separated the individuals of RP from other ones. The second significant axis (accounting for 27.5% of the total variance) was predominantly defined by geographical distance with db-MEM1, secondly by oceanographic connectivity with AEM-2 and AEM-4 predictors. The geographical distance enabled the distinction between ALT-LC populations and the other ones of the Bay of Marseille, while oceanographic connectivity drove the genetic differentiation of the CAV, TIB, and VEY populations.
Figure 8. Distance-based RDA (db-RDA) using the set of 10,755 putatively neutral SNPs at the Marseille-La Ciotat site. The significant predictors are indicated by the directional vectors in the biplot (RDA axes 1 and 2) and individuals are represented as dots and are colored according to their population reported on the sampling map (bottom left corner).
For the adaptive genetic structure (i.e., 12 potentially selected loci), seawater temperature (Min-T) was selected by the stepwise forward selection and included in the global Distance Based RDA. The Min-T predictor contributed to the differentiation between the deepest population (LC, depth 45 m) and the other ones (see Supplementary Material 2).
Discussion
Populations of Ericaria zosteroides are genetically differentiated from each other, including at short distances, as reported in the Mediterranean related species Ericaria amentacea and Gongolaria elegans (synonym: Cystoseira elegans) (Susini et al., 2007; Thibaut et al., 2016; Buonomo et al., 2017; Medrano et al., 2020). Such local genetic differentiation has already been observed in other brown algae, such as in the Mediterranean Kelp Laminaria rodriguezii (Reynes et al., 2020), the Atlantic Kelp Laminaria digitata (Billot et al., 2003; Robuchon et al., 2014), Saccharina latissima (Mao et al., 2020), Undaria pinnatifida (Grulois et al., 2011; Guzinski et al., 2018; Salamon et al., 2020), as well as in the Fucoids Hormosira banksii (Coleman et al., 2018), and Fucus ceranoides (Neiva et al., 2012). These studies, from low to high-throughput genotyping approaches, point to the important genetic structuring generally observed in these species. Unexpectedly, the highest level of differentiation observed here did not correspond to the sites separated by the largest geographical distance [i.e., the Corsica sample from all other sites (∼216 km)], but to the RP site from all others of the Bay of Marseille, sometimes separated by a spatial distance of only a few kilometers (∼3 km). Such unexpected patterns of genetic structure can be the result of incorrect species identification (Pante et al., 2015). For example, the differentiation of populations of the Lithophyllum calcareous red algae in the same area is shaped by a combination of cryptic species and intra-specific differentiation (De Jode et al., 2019). Within the Ericaria genus, the use of microsatellite loci challenged the traditional species delimitation in the E. selaginoides species complex (Bermejo et al., 2018, as E. tamariscifolia). However, morphological characters of E. zosteroides are clear-cut, and cryptic species within our sampling are unlikely as indicated by the molecular comparison with other congeneric species. All our analysis of connectivity then corresponds to intra-specific comparisons.
We clearly showed that bottom currents can play a pivotal role in the connectivity of the populations via the dispersal of fragments. Thus, the genetic structure of E. zosteroides was better explained by oceanographic connectivity (adjusted R2adj = 11%) than geographical distance (adjusted R2adj = 3%) using a variation partitioning method, a commonly used method in this kind of study (e.g., Benestan et al., 2016; Dalongeville et al., 2017; Xuereb et al., 2018a,b; Coscia et al., 2020). Fucoids and kelps are usually known to disperse via sea-surface rafting of fragments of the thallus (e.g., Coleman et al., 2011; Hawes et al., 2017; Durrant et al., 2018). Unlike the very shallow species E. amentacea, for which dispersal at the sea-surface has been highlighted directly with zygotes on a limited time dispersal (12 h) or with fragments entangled in floating debris or via seagull (Thibaut et al., 2016; Buonomo et al., 2017), the dispersal of E. zosteroides is most probably done via fertile fragments of thallus rolling on the sea-floor. Dispersal of zygotes of E. zosteroides seems to be limited to 10 m (Capdevila et al., 2018), and vegetative reproduction does not occur in Ericaria spp. This mode of dispersal leads to episodic gene flow among spatially distant populations. Dispersal events above a few kilometers distance appear to be rare and require on average more than 10 days of transport time. In addition to the randomness of the occurrence of currents, topography and unfavorable habitats to the development of E. zosteroides (Posidonia oceanica beds, sand, coralligenous) can become traps for these fragments and therefore considerably limit gene flow between populations, below what is suggested by our modeling approach.
The highly divergent population of RP in the Bay of Marseille is isolated by oceanographic connectivity and characterized by a low level of genetic variation and significant heterozygote deficiency. This population has the lowest individual density. All of these elements support the fact that the genetic diversity in RP is driven by strong local genetic drift and probably high self-recruitment. A meta-analysis of the genetic structure of populations in this area has shown for several species a barrier to gene flow near the Frioul Archipelago, where RP is located (Cahill et al., 2017), but differences in sampling locations limit a more precise comparison with these results. The evolutionary history of the populations should be considered as well. The RP population could belong to a different genetic cluster that was previously differentiated, and the contact zone of this cluster with other populations could be trapped in this area of reduced oceanographic connectivity (Bierne et al., 2011).
More generally, the strong genetic structure of populations of E. zosteroides could be explained by another factor: the priority effect (De Meester et al., 2002). This effect seems to be relevant to explain high levels of genetic differentiation in Fucoids (e.g., Neiva et al., 2012; Thibaut et al., 2016) but further investigations, including additional sampling will be needed to test this specific hypothesis for E. zosteroides populations. The priority effect states that gene flow can be restricted not by dispersal but by the density and the size of individuals, which can limit new recruitments. Within E. zosteroides populations the recruitment is highly regulated by density-dependent effects (e.g., Ballesteros et al., 2009; Capdevila et al., 2015). If the mean density of our populations is within the range of other E. zosteroides populations measured in other North-Western Mediterranean populations (Thibaut et al., 2005; Hereu et al., 2008; Ballesteros et al., 2009; Capdevila et al., 2015), the size range of the main axis is much smaller than those recorded in the MPA of Scandola in Corsica (Ballesteros et al., 2009) and Port-Cros (Hereu et al., 2008), and similar to those observed at the Balearic and Catalonia sites (Capdevila et al., 2015). The estimated age of the individuals studied here is lower than those recorded in the MPA of Scandola (Corsica) (Ballesteros et al., 2009). When individuals of E. zosteroides are extirpated during episodic disturbance events, the loss of individuals and free available space can promote recruitment (e.g., Navarro et al., 2011; Capdevila et al., 2015; Capdevila Lanzaco, 2017). The population studied on the isolated Armoire bank (ARM) had the highest density of individuals recorded in the study. Additionally, no recruits were observed here. The high genetic differentiation observed for this location (the second highest) may be driven by the combined effects of spatial isolation and priority effect, the latter reducing the success of zygote installation. In the case of RP, despite low density, oceanographic connectivity probably remains too low to allow recruitment. To test the effect of demography on the genetic structure of this species, a dedicated study should be performed with the analysis of more populations with different densities, if possible not correlated with the oceanographic connectivity.
The ARM site was also characterized by an excess of heterozygotes compared to panmixia, which is usually not observed in Fucoids (Coleman and Brawley, 2005; Thibaut et al., 2016; Buonomo et al., 2017; Medrano et al., 2020). The mating system of Fucoids should rather create homozygous zygotes through selfing. The excess of heterozygotes observed here, but at this site only, could be explained by several non-mutually exclusive hypotheses. If few breeders contribute to the next generation, allelic frequencies may differ between individuals contributing to the male and female gamete pools because of binomial sampling error: this can create heterozygosity excess (Pudovkin et al., 1996). This effect can be more important in the case of self-incompatibility (Stoeckel et al., 2006). The selective advantage of heterozygotes compared to homozygotes (i.e., heterosis) could explain the negative FIS as has been suggested for the red algae Gracilaria chilensis (Guillemin et al., 2008). To test this last hypothesis, it would be interesting to study genotypic change between recruits and the oldest individuals. Heterozygote excess can also be observed under partial clonal reproduction (e.g., Stoeckel and Masson, 2014), but this has never been observed in E. zosteroides.
The observation of a significant genetic differentiation at ARM among quite close depths (20 vs. 40 m) indicates a structuring along an environmental gradient which corresponds to differences in light intensity and quality, and in thermal conditions (Pratlong et al., 2018). The observation of a vertical genetic structure at that scale has already been reported in this area in three octocoral species based on microsatellite data (Ledoux et al., 2010; Mokhtar-Jamaï et al., 2011; Cánovas-Molina et al., 2018), and more recently with RAD sequencing in Corallium rubrum (Pratlong et al., 2018). The level of genetic differentiation reported between shallow and deep individuals (ARM) are of the same order of magnitude as populations at the same depth (between CAV and TIB) but 5.8 km apart. This differentiation according to depth can be driven by a combination of limited dispersal of sinking non-buoyant zygotes and fragments (Guern, 1962; Capdevila et al., 2018) and genetic drift, by an effect of local adaptation for some loci, or by intrinsic genetic incompatibilities coupled with environmental barriers along this depth gradient (Bierne et al., 2011). Interestingly, our results suggest a possible effect of temperature on the genetic structure at outlier loci: nevertheless, this subject would require the analysis of more comparisons among different depths and thermal conditions.
Several studies had previously reported the colonization of new habitats by marine species previously thought to be low dispersal, such as Ericaria species, which can colonize habitats several kilometers away (Thibaut et al., 2014). The individuals of ALT (estimation from 8 to 14 years) were younger than those from other populations studied here. For this long-living species, the rate of natural mortality is extremely low in absence of disturbance (lower than 2%, Ballesteros et al., 2009). These results seem to indicate that ALT is one of the youngest populations of the study, which corroborates well with the first report of the species in 2014 (Bonhomme et al., 2014) since the installation of the pipelines in 1967. The high density of this younger population does not preclude the possibility of founder effects during colonization, with the reduction in genetic variation. However, the genetic variation (He) was not reduced in the founder population in comparison to the levels reported in other ones in the Bay of Marseille. This finding was similar to those in the recently founded populations of the Fucoid Gongolaria elegans, which showed no signs of reduction in genetic diversity (Medrano et al., 2020). According to pairwise oceanographic connectivity and genetic structure, the population of La Ciotat (LC) 9 km away could be the donor population.
The application of genomics in conservation strategy will be essential for monitoring species and taking decisive actions by prioritizing the evolutionary potential of populations (von der Heyden, 2017; Xuereb et al., 2020). The use of genomic information coupled with in situ measurements of the demography and modeling of dispersal is seen to be highly effective to unravel occasional connectivity events in supposed low dispersal marine species. Based on our findings, we propose that such an integrative framework must be the proper way to unravel the relative influence of environmental factors on the genetic structure of marine forests. For E. zosteroides in some cases, the species can disperse over long distances and this new information implies another strategy for it is conservation. Thus, apart from including species in the list of protection (all the Mediterranean Fucales have been listed under the Barcelona convention as strictly protected), none of these lists of species have been transcribed into national law and the difficulty of identifying species makes this measure ineffective (Verlaque et al., 2019). Currently, the marine forests are surveyed in MPAs but there is no conservation plan to effectively protect these populations (lack of management of the impact of herbivores and fishing activities on marine forests). Ericaria zosteroides is one of the most long-living erect seaweeds in the Mediterranean, and although rarely impacted by severe natural disturbances (Capdevila et al., 2015, 2019) most of the time, the species is impacted by cumulative direct anthropogenic stressors such as the decline in water transparency (Thibaut et al., 2005) and fishing activities (Hereu et al., 2008). Our results show that some populations of E. zosteroides are spatially isolated and genetically differentiated but some of them remain connected through occasional dispersal. Conservation should be focused at the population level and with a priority of protection for the most isolated populations. The most effective measure is to ban fishing activities and then to set up surveillance of the increase in turbidity mainly due to untreated sewage outfall. None of the studied populations are within a No Take Zone, except VEY. Even though the other populations near Marseille—TIB, VEY, ALT, LC—are in the core area of the Parc National des Calanques, and RP, CAV are within its “adjacent area,” all the populations are subject to anthropic disturbances. Our results stress that conservation strategies should be focused on the implementation of the proposed measures, allowing the preservation in priority of the evolutionary potential and the ecological function of these isolated populations. Monitoring the temporal connectivity patterns of these populations will be necessary to determine whether isolated and vulnerable populations need assisted gene flow from a restoration perspective.
Data Availability Statement
The datasets presented in this study can be found in online repositories. This data can be found here: Genbank, with BioProject ID PRJNA722107 (http://www.ncbi.nlm.nih.gov/bioproject/722107).
Author Contributions
TT and DA devised the project and were responsible for the main conceptual ideas. AB, SR, and SS performed the sampling and enabled the acquisition of the demographic parameters. CC, CP, and LR performed the analyses using Lagrangian models. SM and LR performed the RAD-sequencing experiment and analyzed the sequencing reads. LR performed the research, analyzed the data, and wrote the manuscript in close collaboration with TT, DA, C-FB, MVe, and MVa. All authors provided critical feedback and helped shape the research, analysis and manuscript.
Funding
This project was funded by the MARFOR Biodiversa/0004/2015 project (http://marfor.eu/), and by the Agence Nationale de la Recherche, Grant/Award No.: ANR-18-CE32-0001 (Clonix2D). LR also received support from the GDR Génomique Environnementale (GDR3692). The project leading to this publication has received funding from the European FEDER Fund under project 1166-39417. This work was also funded by a Ph.D. Fellowship from Région PACA and the Calanques National Park.
Conflict of Interest
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
Acknowledgments
We would like to thank the diving unit of OSU Institut Pythéas, the divers of the National Institute for Ocean Science (IFREMER), and the divers of GIS Posidonie for sampling Ericaria zosteroides populations, the Biogenouest genomics core facility (Genomer Plateforme Génomique at the Station Biologique de Roscoff) for their technical support, the staff of the “Cluster de calcul intensif HPC” Platform of OSU Institut Pythéas (Aix-Marseille University, INSU-CNRS) for providing the computing facilities, and finally Laura Benestan for her invaluable support in producing Asymmetric Eigenvector Maps (AEM). Finally, we would like to thank Michael Paul for English language proofreading.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmars.2021.683528/full#supplementary-material
References
Anderson, E. C. (2010). Assessing the power of informative subsets of loci for population assignment: standard methods are upwardly biased. Mol. Ecol. Resour. 10, 701–710. doi: 10.1111/j.1755-0998.2010.02846.x
Anderson, E. C., Waples, R. S., and Kalinowski, S. T. (2008). An improved method for predicting the accuracy of genetic stock identification. Can. J. Fish. Aquat. Sci. 65, 1475–1486. doi: 10.1139/F08-049
Andrews, S. (2010). FastQC: A Quality Control Tool for High Throughput Sequence Data. (Version 0.11.7). Available online at: http://www.bioinformatics.babraham.ac.uk/projects/fastqc (accessed November 24, 2020).
Ballesteros, E. (1990). Structure and dynamics of the community of Cystoseira zosteroides (Turner) C. Agardh (Fucales, Phaeophyceae) in the Northwestern Mediterranean. Sci. Mar. 54, 217–229.
Ballesteros, E. (2006). Mediterranean coralligenous assemblages: a synthesis of present knowledge. Oceanogr. Mar. Biol. 44, 123–195.
Ballesteros, E., Garrabou, J., Hereu, B., Zabala, M., Cebrian, E., and Sala, E. (2009). Deep-water stands of Cystoseira zosteroides C. Agardh (Fucales, Ochrophyta) in the Northwestern Mediterranean: insights into assemblage structure and population dynamics. Estuar. Coast. Shelf Sci. 82, 477–484.
Barrier, N., Petrenko, A. A., and Ourmières, Y. (2016). Strong intrusions of the Northern mediterranean current on the eastern Gulf of Lion: insights from in-situ observations and high resolution numerical modelling. Ocean Dyn. 66, 313–327. doi: 10.1007/s10236-016-0921-7
Benestan, L., Quinn, B. K., Maaroufi, H., Laporte, M., Clark, F. K., Greenwood, S. J., et al. (2016). Seascape genomics provides evidence for thermal adaptation and current-mediated population structure in American lobster (Homarus americanus). Mol. Ecol. 25, 5073–5092. doi: 10.1111/mec.13811
Bermejo, R., Chefaoui, R. M., Engelen, A. H., Buonomo, R., Neiva, J., Ferreira-Costa, J., et al. (2018). Marine forests of the Mediterranean-Atlantic Cystoseira tamariscifolia complex show a southern Iberian genetic hotspot and no reproductive isolation in parapatry. Sci. Rep. 8:10427. doi: 10.1038/s41598-018-28811-1
Bierne, N., Welch, J., Loire, E., Bonhomme, F., and David, P. (2011). The coupling hypothesis: why genome scans may fail to map local adaptation genes. Mol. Ecol. 20, 2044–2072. doi: 10.1111/j.1365-294X.2011.05080.x
Billot, C., Engel, C. R., Rousvoal, S., Kloareg, B., and Valero, M. (2003). Current patterns, habitat discontinuities and population genetic structure: the case of the kelp Laminaria digitata in the English Channel. Mar. Ecol. Prog. Ser. 253, 111–121.
Blanchet, F. G., Legendre, P., and Borcard, D. (2008). Modelling directional spatial processes in ecological data. Ecol. Modell. 215, 325–336. doi: 10.1016/j.ecolmodel.2008.04.001
Boero, F. (2010). The study of species in the era of biodiversity: a tale of stupidity. Diversity 2, 115–126.
Bolger, A. M., Lohse, M., and Usadel, B. (2014). Trimmomatic: a flexible trimmer for illumina sequence data. Bioinformatics 30, 2114–2120. doi: 10.1093/bioinformatics/btu170
Bonhomme, P., Goujard, A., Javel, A., Grondin, J., and Boudouresque, C. F. (2014). “Unexpected artificial-reef-like effect due to a Mediterranean pipeline and the conservation of two circalittoral emblematic species: Centrostephanus longispinus and Cystoseira zosteroides,” in Proceedings of the 2nd Mediterranean Symposium on the Conservation of Coralligenous and Other Calcareous Bio-Concretions, eds C. Bouafif, H. Langar, and A. Ouerghi (Tunis: RAC/SPA Publication), 43–48.
Boudouresque, C. F., Blanfuné, A., Harmelin-Vivien, M., Personnic, S., Ruitton, S., Thibaut, T., et al. (2017). “Where seaweed forests meet animal forests: the examples of macroalgae in coral reefs and the Mediterranean coralligenous ecosystem,” in Marine Animal Forests, eds S. Rossi, L. Bramanti, A. Gori, and C. Orejas (Switzerland: Springer International Publishing), 369–396.
Boudouresque, C. F., Blanfuné, A., Ruitton, S., and Thibaut, T. (2020). “Macroalgae as a tool for coastal management in the Mediterranean Sea,” in Handbook of Algal Science, Technology and Medicine, ed. O. Konur (London: Elsevier and Academic Press).
Bracco, A., Liu, G., Galaska, M. P., Quattrini, A. M., and Herrera, S. (2019). Integrating physical circulation models and genetic approaches to investigate population connectivity in deep-sea corals. J. Mar. Syst. 198:103189. doi: 10.1016/j.jmarsys.2019.103189
Buonomo, R., Assis, J., Fernandes, F., Engelen, A. H., Airoldi, L., and Serrão, E. A. (2017). Habitat continuity and stepping-stone oceanographic distances explain population genetic connectivity of the brown alga Cystoseira amentacea. Mol. Ecol. 26, 766–780. doi: 10.1111/mec.13960
Cahill, A. E., Jode, A. D., Dubois, S., Bouzaza, Z., Aurelle, D., Boissin, E., et al. (2017). A multispecies approach reveals hot spots and cold spots of diversity and connectivity in invertebrate species with contrasting dispersal modes. Mol. Ecol. 26, 6563–6577. doi: 10.1111/mec.14389
Cánovas-Molina, A., Montefalcone, M., Bavestrello, G., Masmoudi, M. B., Haguenauer, A., Hammami, P., et al. (2018). From depth to regional spatial genetic differentiation of Eunicella cavolini in the NW Mediterranean. C. R. Biol. 341, 421–432. doi: 10.1016/j.crvi.2018.09.002
Capdevila, P., Hereu, B., Salguero−Gómez, R., Rovira, G., Medrano, A., Cebrian, E., et al. (2019). Warming impacts on early life stages increase the vulnerability and delay the population recovery of a long-lived habitat-forming macroalga. J. Ecol. 107, 1129–1140. doi: 10.1111/1365-2745.13090
Capdevila, P., Linares, C., Aspillaga, E., Navarro, L., Kersting, D. K., and Hereu, B. (2015). Recruitment patterns in the Mediterranean deepwater alga Cystoseira zosteroides. Mar. Biol. 164, 1665–1674.
Capdevila, P., Linares, C., Aspillaga, E., Riera, J. L., and Hereu, B. (2018). Effective dispersal and density-dependence in mesophotic macroalgal forests: insights from the Mediterranean species Cystoseira zosteroides. PLoS One 13:e0191346. doi: 10.1371/journal.pone.0191346
Capdevila Lanzaco, P. (2017). Life History, Population Dynamics and Conservation of Underwater Mediterranean Forests: Insights From the Long-Lived Alga Cystoseira zosteroides = Història de Vida, Ecologia de Poblacions i Conservació dels Boscos Submergits del Mediterrani: el cas de l’alga longeva Cystoseira zosteroides. Doctoral dissertation. Barcelona: Universitat de Barcelona, 274.
Catchen, J. M., Amores, A., Hohenlohe, P., Cresko, W., and Postlethwait, J. H. (2011). Stacks: building and genotyping Loci De Novo from short-read sequences. G3 1, 171–182. doi: 10.1534/g3.111.000240
Chang, C. C., Chow, C. C., Tellier, L. C., Vattikuti, S., Purcell, S. M., and Lee, J. J. (2015). Second-generation PLINK: rising to the challenge of larger and richer datasets. GigaScience 4:7. doi: 10.1186/s13742-015-0047-8
Coleman, M. A., and Brawley, S. H. (2005). Are life history characteristics good predictors of genetic diversity and structure? A case study of the Intertidal alga Fucus spiralis (heterokontophyta; Phaeophyceae). J. Phycol. 41, 753–762. doi: 10.1111/j.0022-3646.2005.04136.x
Coleman, M. A., Clark, J. S., Doblin, M. A., Bishop, M. J., and Kelaher, B. P. (2018). Genetic differentiation between estuarine and open coast ecotypes of a dominant ecosystem engineer. Mar. Freshw. Res. 70, 977–985.
Coleman, M. A., Roughan, M., Macdonald, H. S., Connell, S. D., Gillanders, B. M., Kelaher, B. P., et al. (2011). Variation in the strength of continental boundary currents determines continent-wide connectivity in kelp. J. Ecol. 99, 1026–1032. doi: 10.1111/j.1365-2745.2011.01822.x
Coscia, I., Wilmes, S. B., Ironside, J. E., Goward−Brown, A., O’Dea, E., Malham, S. K., et al. (2020). Fine-scale seascape genomics of an exploited marine species, the common cockle Cerastoderma edule, using a multimodelling approach. Evol. Appl. 13, 1854–1867. doi: 10.1111/eva.12932
Dalongeville, A., Andrello, M., Mouillot, D., Lobreaux, S., Fortin, M.-J., Lasram, F., et al. (2017). Geographic isolation and larval dispersal shape seascape genetic patterns differently according to spatial scale. Evol. Appl. 11, 1437–1447.
De Jode, A., David, R., Haguenauer, A., Cahill, A. E., Erga, Z., Guillemain, D., et al. (2019). From seascape ecology to population genomics and back. Spatial and ecological differentiation among cryptic species of the red algae Lithophyllum stictiforme/L. cabiochiae, main bioconstructors of coralligenous habitats. Mol. Phylogenet. Evol. 137, 104–113.
De Meester, L., Gómez, A., Okamura, B., and Schwenk, K. (2002). The monopolization hypothesis and the dispersal–gene flow paradox in aquatic organisms. Acta Oecol. 23, 121–135. doi: 10.1016/S1146-609X(02)01145-1
Dray, S., Legendre, P., and Peres-Neto, P. R. (2006). Spatial modelling: a comprehensive framework for principal coordinate analysis of neighbour matrices (PCNM). Ecol. Model. 196, 483–493. doi: 10.1016/j.ecolmodel.2006.02.015
Durrant, H. M. S., Barrett, N. S., Edgar, G. J., Coleman, M. A., and Burridge, C. P. (2018). Seascape habitat patchiness and hydrodynamics explain genetic structuring of kelp populations. Mar. Ecol. Progr. Ser. 587, 81–92. doi: 10.3354/meps12447
Falace, A., Kaleb, S., De La Fuente, G., Asnaghi, V., and Chiantore, M. (2018). Ex situ cultivation protocol for Cystoseira amentacea var. stricta (Fucales, Phaeophyceae) from a restoration perspective. PLoS One 13:e0193011. doi: 10.1371/journal.pone.0193011
Foll, M., and Gaggiotti, O. (2008). A genome-scan method to identify selected loci appropriate for both dominant and codominant markers: a Bayesian perspective. Genetics 180, 977–993. doi: 10.1534/genetics.108.092221
Fraysse, M., Pairaud, I., Ross, O. N., Faure, V., and Pinazo, C. (2014). Intrusion of Rhone River diluted water into the Bay of Marseille: generation processes and impacts on ecosystem functioning. J. Geophys. Res. 119, 6535–6556. doi: 10.1002/2014JC010022
Fraysse, M., Pinazo, C., Faure, V. M., Fuchs, R., Lazzari, P., Raimbault, P., et al. (2013). Development of a 3D coupled physical-biogeochemical model for the Marseille coastal area (NW Mediterranean Sea): what complexity is required in the coastal zone? PLoS One 8:e80012. doi: 10.1371/journal.pone.0080012
Frichot, E., and François, O. (2015). LEA: an R package for landscape and ecological association studies. Methods Ecol. Evol. 6, 925–929. doi: 10.1111/2041-210X.12382
Frichot, E., Mathieu, F., Trouillon, T., Bouchard, G., and François, O. (2014). Fast and efficient estimation of individual ancestry coefficients. Genetics 196, 973–983. doi: 10.1534/genetics.113.160572
Gagnaire, P., Broquet, T., Aurelle, D., Viard, F., Souissi, A., Bonhomme, F., et al. (2015). Using neutral, selected, and hitchhiker loci to assess connectivity of marine populations in the genomic era. Evol. Appl. 8, 769–786. doi: 10.1111/eva.12288
Gagnaire, P. A. (2020). Comparative genomics approach to evolutionary process connectivity. Evol. Appl. 13, 1320–1334.
Giaccone, G. (1973). Écologie et chorologie des Cystoseira de Méditerranée. Rapp. Comm. Intl. Mer Médit. 22, 49–50.
Gosselin, T., Anderson, E. C., and Bradbury, L. (2016). Assigner: Assignment Analysis with GBS/RAD Data Using R. Rpackage Version 0.4.1.
Gotelli, N. J. (1991). Metapopulation models: the rescue effect, the propagule rain, and the core-satellite hypothesis. Am. Nat. 138, 768–776. doi: 10.1086/285249
Grulois, D., Lévêque, L., Viard, F., Frangoudes, K., and Valero, M. (2011). Mosaic genetic structure and sustainable establishment of the invasive kelp Undaria pinnatifida within a bay (Bay of St-Malo, Brittany). Cah. Biol. Mar. 52:485.
Guern, M. (1962). Embryologie de quelques espèces du genre Cystoseira Agardh 1821 (Fucales). Vie et Milieu. 13, 649–679.
Guillemin, M. L., Faugeron, S., Destombe, C., Viard, F., Correa, J. A., and Valero, M. (2008). Genetic variation in wild and cultivated populations of the haploid–diploid red alga Gracilaria chilensis: how farming practices favor asexual reproduction and heterozygosity. Evolution 62, 1500–1519.
Guzinski, J., Ballenghien, M., Daguin−Thiébaut, C., Lévêque, L., and Viard, F. (2018). Population genomics of the introduced and cultivated Pacific kelp Undaria pinnatifida: marinas—not farms—drive regional connectivity and establishment in natural rocky reefs. Evol. Appl. 11, 1582–1597.
Hawes, N. A., Taylor, D. I., and Schiel, D. R. (2017). Transport of drifting fucoid algae: nearshore transport and potential for long distance dispersal. J. Exp. Mar. Biol. Ecol. 490, 34–41. doi: 10.1016/j.jembe.2017.02.001
Hereu, B., Mangialajo, L., Ballesteros, E., and Thibaut, T. (2008). On the occurrence, structure and distribution of deep-water Cystoseira (Phaeophyceae) populations in the Port-Cros National Park (north-western Mediterranean). Eur. J. Phycol. 43, 263–273. doi: 10.1080/09670260801930330
Hoban, S., Bruford, M., D’Urban Jackson, J., Lopes-Fernandes, M., Heuertz, M., Hohenlohe, P. A., et al. (2020). Genetic diversity targets and indicators in the CBD post-2020 global biodiversity framework must be improved. Biol. Conserv. 248:108654. doi: 10.1016/j.biocon.2020.108654
Keller, L., and Waller, D. M. (2002). Inbreeding effects in wild populations. Trends Ecol. Evol. 17, 230–241. doi: 10.1016/S0169-5347(02)02489-8
Kinlan, B. P., and Gaines, S. D. (2003). Propagule dispersal in marine and terrestrial environments: a community perspective. Ecology 84, 2007–2020. doi: 10.1890/01-0622
Lazure, P., and Dumas, F. (2008). An external–internal mode coupling for a 3D hydrodynamical model for applications at regional scale (MARS). Adv. Water Resour. 31, 233–250. doi: 10.1016/j.advwatres.2007.06.010
Ledoux, J. B., Garrabou, J., Bianchimani, O., Drap, P., Féral, J. P., and Aurelle, D. (2010). Fine−scale genetic structure and inferences on population biology in the threatened Mediterranean red coral, Corallium rubrum. Mol. Ecol. 19, 4204–4216. doi: 10.1111/j.1365-294X.2010.04814.x
Legendre, P., and Anderson, M. J. (1999). Distance-based redundancy analysis: testing multispecies responses in multifactorial ecological experiments. Ecol. Monogr. 69:24.
Lenormand, T. (2002). Gene flow and the limits to natural selection. Trends Ecol. Evol. 17, 183–189.
Li, H., and Durbin, R. (2009). Fast and accurate short read alignment with burrows-wheeler transform. Bioinformatics 25, 1754–1760. doi: 10.1093/bioinformatics/btp324
Lowe, W. H., and Allendorf, F. W. (2010). What can genetics tell us about population connectivity? Mol. Ecol. 19, 3038–3051. doi: 10.1111/j.1365-294X.2010.04688.x
Luu, K., Bazin, E., and Blum, M. G. B. (2017). pcadapt: an R package to perform genome scans for selection based on principal component analysis. Mol. Ecol. Resour. 17, 67–77. doi: 10.1111/1755-0998.12592
Mangialajo, L., Chiantore, M., Susini, M. L., Meinesz, A., Cattaneo-Vietti, R., and Thibaut, T. (2012). Zonation patterns and interspecific relationships of fucoids in microtidal environments. J. Exp. Mar. Biol. Ecol. 412, 72–80. doi: 10.1016/j.jembe.2011.10.031
Manichaikul, A., Mychaleckyj, J. C., Rich, S. S., Daly, K., Sale, M., and Chen, W.-M. (2010). Robust relationship inference in genome-wide association studies. Bioinformatics 26, 2867–2873. doi: 10.1093/bioinformatics/btq559
Mao, X., Augyte, S., Huang, M., Hare, M. P., Bailey, D., Umanzor, S., et al. (2020). Population genetics of sugar kelp throughout the Northeastern United States using genome-wide markers. Front. Mar. Sci. 7:694. doi: 10.3389/fmars.2020.00694
Martins, H., Caye, K., Luu, K., Blum, M. G. B., and François, O. (2016). Identifying outlier loci in admixed and in continuous populations using ancestral population differentiation statistics. Mol. Ecol. 25, 5029–5042. doi: 10.1111/mec.13822
Maruki, T., and Lynch, M. (2017). Genotype calling from population-genomic sequencing data. G3 7, 1393–1404. doi: 10.1534/g3.117.039008
Medrano, A., Hereu, B., Mariani, S., Neiva, J., Pagès-Escolà, M., Paulino, C., et al. (2020). Ecological traits, genetic diversity and regional distribution of the macroalga Treptacantha elegans along the Catalan coast (NW Mediterranean Sea). Sci. Rep. 10:19219. doi: 10.1038/s41598-020-76066-6
Millet, B., Pinazo, C., Banaru, D., Pagès, R., Guiart, P., and Pairaud, I. (2018). Unexpected spatial impact of treatment plant discharges induced by episodic hydrodynamic events: modelling Lagrangian transport of fine particles by Northern current intrusions in the bays of Marseille (France). PLoS One 13:e0195257. doi: 10.1371/journal.pone.0195257
Mokhtar-Jamaï, K., Pascual, M., Ledoux, J. B., Coma, R., Féral, J. P., Garrabou, J., et al. (2011). From global to local genetic structuring in the red gorgonian Paramuricea clavata: the interplay between oceanographic conditions and limited larval dispersal. Mol. Ecol. 20, 3291–3305. doi: 10.1111/j.1365-294X.2011.05176.x
Molinari Novoa, E. A., and Guiry, M. D. (2020). Reinstatement of the genera gongolaria boehmer and ericaria stackhouse (Sargassaceae, Phaeophyceae). Notul. Algar. 171, 1–10.
Navarro, L., Ballesteros, E., Linares, C., and Hereu, B. (2011). Spatial and temporal variability of deep-water algal assemblages in the Northwestern Mediterranean: the effects of an exceptional storm. Estuar. Coast. Shelf Sci. 95, 52–58. doi: 10.1016/j.ecss.2011.08.002
Neiva, J., Pearson, G. A., Valero, M., and Serrão, E. A. (2012). Fine-scale genetic breaks driven by historical range dynamics and ongoing density-barrier effects in the estuarine seaweed Fucus ceranoides. BMC Evol. Biol. 12:78. doi: 10.1186/1471-2148-12-78
Nicolle, A., Garreau, P., and Liorzou, B. (2009). Modelling for anchovy recruitment studies in the Gulf of Lions (Western Mediterranean Sea). Ocean Dyn. 59, 953–968. doi: 10.1007/s10236-009-0221-6
Oksanen, J., Blanchet, F. G., Friendly, M., Kindt, R., Legendre, P., Mcglinn, D., et al. (2017). Community ecology package ‘vegan’. Version 2.4-4. R Package 2, 1–192.
Pairaud, I. L., Gatti, J., Bensoussan, N., Verney, R., and Garreau, P. (2011). Hydrology and circulation in a coastal area off Marseille: validation of a nested 3D model with observations. J. Mar. Syst. 88, 20–33. doi: 10.1016/j.jmarsys.2011.02.010
Palumbi, S. R. (2003). Population genetics, demographic connectivity, and the design of marine reserves. Ecol. Appl. 13(Suppl.1), 146–158. doi: 10.1890/1051-07612003013[0146:PGDCAT]2.0.CO;2
Pante, E., Puillandre, N., Viricel, A., Arnaud-Haond, S., Aurelle, D., Castelin, M., et al. (2015). Species are hypotheses: avoid connectivity assessments based on pillars of sand. Mol. Ecol. 24, 525–544. doi: 10.1111/mec.13048
Paterno, M., Schiavina, M., Aglieri, G., Ben Souissi, J., Boscari, E., Casagrandi, R., et al. (2017). Population genomics meet Lagrangian simulations: oceanographic patterns and long larval duration ensure connectivity among Paracentrotus lividus populations in the adriatic and ionian seas. Ecol. Evol. 7, 2463–2479. doi: 10.1002/ece3.2844
Peterson, B. K., Weber, J. N., Kay, E. H., Fisher, H. S., and Hoekstra, H. E. (2012). Double Digest RADseq: an inexpensive method for de novo SNP discovery and genotyping in model and non-model species. PLoS One 7:e37135. doi: 10.1371/journal.pone.0037135
Pradal, M. A., and Millet, B. (2006). Hétérogénéité spatiale du forçage du vent: fonctionnement des récifs artificiels et circulation des eaux dans la baie sud de Marseille. C. R. Biol. 329, 541–550. doi: 10.1016/j.crvi.2006.04.002
Pratlong, M., Haguenauer, A., Brener, K., Mitta, G., Toulza, E., Garrabou, J., et al. (2018). Separate the wheat from the chaff: genomic scan for local adaptation in the red coral Corallium rubrum. BioRxiv doi: 10.1101/306456
Privé, F., Aschard, H., Ziyatdinov, A., and Blum, M. G. B. (2018). Efficient analysis of large-scale genome-wide data with two R packages: bigstatsr and bigsnpr. Bioinformatics 34, 2781–2787. doi: 10.1093/bioinformatics/bty185
Pudovkin, A. I., Zaykin, D. V., and Hedgecock, D. (1996). On the potential for estimating the effective number of breeders from heterozygote-excess in progeny. Genetics 144, 383–387.
Raymond, M., and Rousset, F. (1995). An exact test for population differentiation. Evolution 49, 1280–1283. doi: 10.1111/j.1558-5646.1995.tb04456.x
Reynes, L., Thibaut, T., Mauger, S., Blanfuné, A., Holon, F., Cruaud, C., et al. (2020). Genomic signatures of clonality in the deep water kelp Laminaria rodriguezii. Authorea Preprint 30, 1806–1822. doi: 10.22541/au.159654446.63647641
Robuchon, M., Le Gall, L., Mauger, S., and Valero, M. (2014). Contrasting genetic diversity patterns in two sister kelp species co−distributed along the coast of Brittany, France. Mol. Ecol. 23, 2669–2685.
Robvieux, P. (2013). Conservation des Populations de Cystoseira en Régions Provence-Alpes-Côte-d’Azur et Corse. Doctoral Dissertation. Nice: University of Nice-Sophia Antipolis.
Rochette, N. C., Rivera−Colón, A. G., and Catchen, J. M. (2019). Stacks 2: analytical methods for paired-end sequencing improve RADseq-based population genomics. Mol. Ecol. 28, 4737–4754. doi: 10.1111/mec.15253
Ross, O. N., Fraysse, M., Pinazo, C., and Pairaud, I. (2016). Impact of an intrusion by the Northern current on the biogeochemistry in the eastern Gulf of Lion, NW Mediterranean. Estuar. Coast. Shelf Sci. 170, 1–9. doi: 10.1016/j.ecss.2015.12.022
Rousset, F. (2008). genepop’007: a complete re-implementation of the genepop software for Windows and Linux. Mol. Ecol. Resour. 8, 103–106. doi: 10.1111/j.1471-8286.2007.01931.x
Rousset, F., and Raymond, M. (1995). Testing heterozygote excess and deficiency. Genetics 140, 1413–1419.
Salamon, M., Levêque, L., Ballenghien, M., and Viard, F. (2020). Spill-back events followed by self-sustainment explain the fast colonization of a newly built marina by a notorious invasive seaweed. Biol. Invas. 22, 1411–1429.
Slatkin, M. (1987). Gene flow and the geographic structure of natural populations. Science 236, 787–792. doi: 10.1126/science.3576198
Stoeckel, S., Grange, J., Fernández−Manjarres, J. F., Bilger, I., Frascaria-Lacoste, N., and Mariette, S. (2006). Heterozygote excess in a self−incompatible and partially clonal forest tree species—Prunus avium L. Mol. Ecol. 15, 2109–2118.
Stoeckel, S., and Masson, J. P. (2014). The exact distributions of FIS under partial asexuality in small finite populations with mutation. PLoS One 9:e85228. doi: 10.1371/journal.pone.0085228
Susini, M. L., Thibaut, T., Meinesz, A., and Forcioli, D. (2007). A preliminary study of genetic diversity in Cystoseira amentacea (C. Agardh) Bory var. stricta Montagne (Fucales, Phaeophyceae) using random amplified polymorphic DNA. Phycologia 46, 605–611. doi: 10.2216/06-100.1
Tarasov, A., Vilella, A. J., Cuppen, E., Nijman, I. J., and Prins, P. (2015). Sambamba: fast processing of NGS alignment formats. Bioinformatics 31, 2032–2034. doi: 10.1093/bioinformatics/btv098
Thibaut, T., Blanfuné, A., Boudouresque, C. F., and Verlaque, M. (2015). Decline and local extinction of fucales in the French Riviera: the harbinger of future extinctions? Medit. Mar. Sci. 16, 206–224.
Thibaut, T., Blanfuné, A., Markovic, L., Verlaque, M., Boudouresque, C. F., Perret-Boudouresque, M., et al. (2014). Unexpected abundance and long-term relative stability of the brown alga Cystoseira amentacea, hitherto regarded as a threatened species, in the north-western Mediterranean Sea. Mar. Pollut. Bull. 89, 305–323. doi: 10.1016/j.marpolbul.2014.09.043
Thibaut, T., Bottin, L., Aurelle, D., Boudouresque, C. F., Blanfuné, A., Verlaque, M., et al. (2016). Connectivity of populations of the seaweed Cystoseira amentacea within the Bay of Marseille (Mediterranean Sea): genetic structure and hydrodynamic connections. Cryptogam. Algol. 37, 233–255. doi: 10.7872/crya/v37.iss4.2016.233
Thibaut, T., Pinedo, S., Torras, X., and Ballesteros, E. (2005). Long-term decline of the populations of fucales (Cystoseira spp. and Sargassum spp.) in the Albères coast (France, North-western Mediterranean). Mar. Pollut. Bull. 50, 1472–1489. doi: 10.1016/j.marpolbul.2005.06.014
Verlaque, M., Boudouresque, C. F., and Perret-Boudouresque, M. (2019). Mediterranean seaweeds listed as threatened under the Barcelona convention: a critical analysis. Sci. Rep. Port Cros Natl. Park 33, 179–214.
von der Heyden, S. (2017). Making evolutionary history count: biodiversity planning for coral reef fishes and the conservation of evolutionary processes. Coral Reefs 36, 183–194. doi: 10.1007/s00338-016-1512-2
Waples, R. S. (1998). Separating the wheat from the chaff: patterns of genetic differentiation in high gene flow species. J. Hered. 89, 438–450.
Weir, B. S., and Cockerham, C. C. (1984). Estimating F-statistics for the analysis of population structure. Evolution 38, 1358–1370. doi: 10.1111/j.1558-5646.1984.tb05657.x
White, C., Selkoe, K. A., Watson, J., Siegel, D. A., Zacherl, D. C., and Toonen, R. J. (2010). Ocean currents help explain population genetic structure. Proc. R. Soc. B 277, 1685–1694. doi: 10.1098/rspb.2009.2214
Xuereb, A., Benestan, L., Normandeau, É, Daigle, R. M., Curtis, J. M. R., Bernatchez, L., et al. (2018a). Asymmetric oceanographic processes mediate connectivity and population genetic structure, as revealed by RADseq, in a highly dispersive marine invertebrate (Parastichopus californicus). Mol. Ecol. 27, 2347–2364. doi: 10.1111/mec.14589
Xuereb, A., D’Aloia, C. C., Andrello, M., Bernatchez, L., and Fortin, M. J. (2020). Incorporating putatively neutral and adaptive genomic data into marine conservation planning. Conserv. Biol. doi: 10.1111/cobi.13609
Xuereb, A., Kimber, C. M., Curtis, J. M. R., Bernatchez, L., and Fortin, M. J. (2018b). Putatively adaptive genetic variation in the giant California sea cucumber (Parastichopus californicus) as revealed by environmental association analysis of restriction-site associated DNA sequencing data. Mol. Ecol. 27, 5035–5048. doi: 10.1111/mec.14942
Keywords: connectivity, Fucales, Mediterranean, marine forest, population genomics, Lagrangian modelling, dispersal, gene flow
Citation: Reynes L, Aurelle D, Chevalier C, Pinazo C, Valero M, Mauger S, Sartoretto S, Blanfuné A, Ruitton S, Boudouresque C-F, Verlaque M and Thibaut T (2021) Population Genomics and Lagrangian Modeling Shed Light on Dispersal Events in the Mediterranean Endemic Ericaria zosteroides (=Cystoseira zosteroides) (Fucales). Front. Mar. Sci. 8:683528. doi: 10.3389/fmars.2021.683528
Received: 21 March 2021; Accepted: 27 April 2021;
Published: 02 June 2021.
Edited by:
Manuel Aranda, King Abdullah University of Science and Technology, Saudi ArabiaReviewed by:
Ludwig Triest, Vrije University Brussel, BelgiumAlicia Dalongeville, Université de Montpellier, France
Copyright © 2021 Reynes, Aurelle, Chevalier, Pinazo, Valero, Mauger, Sartoretto, Blanfuné, Ruitton, Boudouresque, Verlaque and Thibaut. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Lauric Reynes, lauric.reynes@mio.osupytheas.fr