Population Genomics and Lagrangian Modeling Shed Light on Dispersal Events in the Mediterranean Endemic Ericaria zosteroides (=Cystoseira zosteroides) (Fucales)

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.

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 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(Slatkin, , 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., 2005Thibaut et al., , 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?

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 northwestern 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 subseapipelines 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 m 2 ). 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 posthoc pairwise Wilcoxon test. The distribution of the length was A significant deficit or excess in heterozygotes (P < 0.001) is indicated as " ‡" or " †,"respectively.
Frontiers in Marine Science | www.frontiersin.org 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.
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.

DNA Extraction and SNP Discovery
Total genomic DNA was isolated with the Nucleospin R 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 (r 2 > 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 F ST 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). F ST 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 F ST 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  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 (H o ), expected heterozygosity (H e ), and F IS (Weir and Cockerham, 1984). We tested for differences in H e between populations with a pairwise Wilcoxon test adjusted for multiple comparisons (Bonferroni adjustment). Positive and negative F IS values were observed, so we tested departure from panmixia across all loci and populations using the Global Hardy-Weinberg test [Score (U) test] , 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  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 F ST 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 holdout 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 spatiotemporal 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 highperformance 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 P ij = N i−j /N i where N i is the total number of released particles from site i and N i−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 P ij . 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).

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 (R 2 adj ) 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, R 2 adj and the p-value for the F-tests were used to evaluate the validity of the models.

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  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 .

Density (ind.m −2 )
Axis height (cm) 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).

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.  (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.

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 (H e ) 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 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 H e and H o 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 F IS ranged from slightly negative values at ARM20 (F IS = -0.02) and ARM40 (F IS = -0.02) to positive values, notably at RP (F IS = 0.04) and ALT (F IS = 0.04).

Genetic Structure and Assignment Tests
The overall estimated F ST was 0.09, and it was slightly lower on a local scale within the Bay of Marseille (F ST = 0.08). Populations were significantly genetically differentiated from each other (exact G test for genic differentiation, P < 0.001), and pairwise F ST ranged from 0.05 to 0.14. The highest F ST values were observed for the comparison of the RP population with other ones (Figure 4). The lowest genetic differentiation (F ST = 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).
The principal component analysis (PCA), and sNMF clustering gave similar results and were congruent with the F ST 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.
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.
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.
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, FIGURE 4 | Heatmap of pairwise F ST 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).
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 (R 2 adj ) 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 R 2 adj = 0.11 against R 2 adj = 0.03 for both geographical distance and density predictors, and R 2 adj = 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.
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 , 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 intraspecific 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 R 2 adj = 11%) than geographical distance (adjusted R 2 adj = 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  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.
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.  Environmental predictors of the neutral genetic structure selected by the stepwise forward selection (ordiR2step) as highlighted in bold are included in the db-RDA framework. All predictors retained in the models are significant at P < 0.001 using ANOVA. The adjusted coefficient of determination (R 2 adj ) and the P-value of the model are reported. Seawater temperature was not selected as a significant predictor in the global db-RDA. 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 F IS 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 longliving 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 (H e ) 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(Capdevila et al., , 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.