Mariana Padrón
email: mariana.padron@gmail.com

Federica Costantini

Lorenzo Bramanti

Katell Guizien

Marco Abbiati

| Federica Costantini

| Katell Guizien

| Marco Abbiati

Genetic connectivity supports recovery of gorgonian populations affected by climate change

Keywords: benthos, climate change, coastal, coral, dispersal, disturbance, genetics, invertebrates

come

 [START_REF] Hoegh-Guldberg | The impact of climate change on the world's marine ecosystems[END_REF]. Trends of increasing seawater temperature, decreasing ocean pH, and sea-level rise, along with the higher frequency and magnitude of extreme climatic events (e.g. storms and heatwaves), are threatening marine species and habitats [START_REF] Doney | Climate change impacts on marine ecosystems[END_REF]. The threshold of populations and habitat resilience to climatic disturbance is still an open question in contemporary ecology and conservation biology, however (Wernberg, Russell, Thomsen, & Connell, 2014).

Severe and sudden reductions in population size (i.e. demographic bottleneck) caused by mass mortality events are expected to lead to a loss of genetic diversity, primarily allelic diversity, and eventually cause a genetic bottleneck [START_REF] Lacy | Loss of genetic diversity from managed populations: Interacting effects of drift, mutation, immigration, selection and population subdivision[END_REF]. Furthermore, within this context, genetic drift may amplify the loss of genetic variation in isolated populations, and potentially lead to inbreeding [START_REF] Nei | The bottleneck effect and genetic variability in populations[END_REF]. Nonetheless, the loss of genetic diversity caused by a population collapse can be offset or reversed by migration processes [START_REF] Crow | An introduction to population genetics theory[END_REF][START_REF] Keller | Immigration and the ephemerality of a natural population bottleneck: Evidence from molecular markers[END_REF]. Therefore, connectivity among populations is an asset for the recovery of the predisturbance demographic structure and genetic diversity of the population.

In marine sessile benthic species, migration between distant populations ensured by larval dispersal is expected to be a key aspect of resilience [START_REF] Almany | Connectivity, biodiversity conservation and the design of marine reserve networks for coral reefs[END_REF]. The role of migration patterns in the recovery after disturbance is controversial, however: although it is possible to forecast wide larval dispersal ranges with biophysical modelling [START_REF] Lipcius | Importance of metapopulation connectivity to restocking and restoration of marine species[END_REF][START_REF] Shanks | Propagule dispersal distance and the size and spacing of marine reserves[END_REF], gene flow, as assessed with molecular tools, appears to be more restricted, as a result of filtering by settlement and recruitment processes (Pineda, Hare, & Sponaugle, 2007 Q6 ). Despite several works assuming the importance of connectivity to the persistence of sessile benthic species populations, empirical evidence on how patterns of larval dispersal confer resilience to populations remains rare [START_REF] Bramanti | Using demographic models to project the effects of climate change on scleractinian corals: Pocillopora damicornis as a case study[END_REF][START_REF] Caplins | Landscape structure and the genetic effects of a population collapse[END_REF][START_REF] Gilmour | Recovery of an isolated coral reef system following severe disturbance[END_REF][START_REF] Magris | Integrating connectivity and climate change into marine conservation planning[END_REF][START_REF] Underwood | Multiple scales of genetic connectivity in a brooding coral on isolated reefs following catastrophic bleaching[END_REF][START_REF] Vollmer | Restricted gene flow in the Caribbean staghorn coral Acropora cervicornis: Implications for the recovery of endangered reefs[END_REF]. Investigating the genetic diversity and the structure of marine populations can provide estimations of their ability to withstand perturbations, and of their rate of recovery [START_REF] Hughes | Ecological consequences of genetic diversity[END_REF][START_REF] Kovach | Genetic diversity is related to climatic variation and vulnerability in threatened bull trout[END_REF]Van Oppen & Gates, 2006 Q7 ). In particular, the assignment of individuals to the population of origin can provide more detailed inference on larval migration patterns among populations [START_REF] Waser | Genetic assignment tests as indices of interpopulation dispersal[END_REF].

During the last two decades, several mass mortality events coinciding with anomalous temperature increases (heatwaves) have affected benthic sessile species along the rocky shores of the northwestern Mediterranean Sea [START_REF] Cupido | Response of a gorgonian (Paramuricea clavata) population to mortality events: Recovery or loss?[END_REF]; but see [START_REF] Rivetti | Global warming and mass mortality of benthic invertebrates in the Mediterranean Sea[END_REF]. Among the species affected, gorgonians (Cnidaria, Octocorallia) were the most damaged, with a local loss of more than 50% in density and biomass [START_REF] Bramanti | Recruitment, early survival and growth of the Mediterranean red coral Corallium rubrum (L 1758), a 4-year study[END_REF][START_REF] Cerrano | A catastrophic mass-mortality episode of gorgonians and other organisms in the Ligurian Sea (Northwestern Mediterranean), summer 1999[END_REF][START_REF] Coma | Global warming-enhanced stratification and mass mortality events in the Mediterranean[END_REF][START_REF] Cupido | Response of a gorgonian (Paramuricea clavata) population to mortality events: Recovery or loss?[END_REF][START_REF] Garrabou | Mass mortality in Northwestern Mediterranean rocky benthic communities: Effects of the 2003 heat wave[END_REF][START_REF] Santangelo | Effects of increased mortality on gorgonian corals (Cnidaria, Octocorallia): Different demographic features may lead affected populations to unexpected recovery and new equilibrium points[END_REF]. Gorgonians are characterized by long life span and slow growth, which along with the low fecundity of small colonies delays their recovery after disturbance [START_REF] Linares | Size distribution, density and disturbance in two Mediterranean gorgonians: Paramuricea clavata and Eunicella singularis[END_REF], even when negative density-dependent control of recruitment is removed [START_REF] Santangelo | Effects of increased mortality on gorgonian corals (Cnidaria, Octocorallia): Different demographic features may lead affected populations to unexpected recovery and new equilibrium points[END_REF]. Similarly to other habitat-forming species, gorgonians have an important role in promoting and stabilizing associated assemblages, such as biogenic reefs. Therefore, the conservation of these assemblages also depends on the capability of gorgonian species to respond to climate change and increasing human disturbance [START_REF] Ponti | Ecological shifts in Mediterranean coralligenous assemblages related to gorgonian forest loss[END_REF].

Among the gorgonian species affected by mass mortalities, the red gorgonian Paramuricea clavata (Risso, 1826) suffered the most drastic reductions in density and biomass over a broad spatial scale, reaching up Q8 to 80% of colonies affected in the Gulf of Genoa in 1993 [START_REF] Rivetti | Global warming and mass mortality of benthic invertebrates in the Mediterranean Sea[END_REF]. The impact of reported mass mortality events on the demographic structure and dynamics of several populations has been assessed [START_REF] Cupido | Unexpected long-term population dynamics in a canopyforming gorgonian coral following mass mortality[END_REF][START_REF] Santangelo | Effects of increased mortality on gorgonian corals (Cnidaria, Octocorallia): Different demographic features may lead affected populations to unexpected recovery and new equilibrium points[END_REF][START_REF] Santangelo | Population dynamics and conservation biology of the over-exploited Mediterranean red coral[END_REF]. The recovery of these populations has been attributed to local replenishment [START_REF] Santangelo | Effects of increased mortality on gorgonian corals (Cnidaria, Octocorallia): Different demographic features may lead affected populations to unexpected recovery and new equilibrium points[END_REF], overlooking the role of larval immigration from other source populations, based on the assumption of low larval dispersal in P. et al., 2015;[START_REF] Pérez-Portela | Lonely populations in the deep: Genetic structure of red gorgonians at the heads of submarine canyons in the north-western Mediterranean Sea[END_REF][START_REF] Pilczynska | Genetic diversity and local connectivity in the Mediterranean red gorgonian coral after mass mortality events[END_REF]. In the Ligurian Sea, genetic analyses conducted on a local scale have suggested that, after mass mortality events, the recovery of P. clavata populations dwelling below 20 m in depth was the result of larval supply from deeper populations (>60 m in depth), presumably not affected by mass mortality [START_REF] Pilczynska | Genetic diversity and local connectivity in the Mediterranean red gorgonian coral after mass mortality events[END_REF]. Nonetheless, although decreases in mortality rates with depth have been observed for this species (Linares et al., 2005), it has been shown that below 20 m of depth P. clavata populations can be affected by mass mortalities [START_REF] Huete-Stauffer | Paramuricea clavata (Anthozoa, Octocorallia) loss in the Marine Protected Area of Tavolara (Sardinia, Italy) due to a mass mortality event[END_REF].

In the present work we extended the geographical scale (from local to regional) and the bathymetrical range (down to 40 m) of sampling in order to test: (i) whether P. clavata populations dwelling between 25 and 40 m of depth in the Ligurian Sea show signatures of a genetic bottleneck as a consequence of mass mortality events; and (ii) the occurrence of gene flow among those populations.

To answer these questions, the genetic diversity, genetic structure, and regional migration patterns of the coastal P. clavata populations in the Ligurian Sea were analysed. T1 1). All samples were preserved in 95% ethanol and stored at 4°C.

The sampling design of P. clavata populations includes sites in which populations of the species were affected by mass mortality events: Portofino, Punta Mesco, and Portovenere. In the Ligurian Sea, mass mortality events as a result of heatwaves were reported in 1999 and 2003. In Portovenere, the pre-disturbance density of 35 individuals/m 2 was reduced by approximately 70% [START_REF] Cupido | Response of a gorgonian (Paramuricea clavata) population to mortality events: Recovery or loss?[END_REF][START_REF] Santangelo | Effects of increased mortality on gorgonian corals (Cnidaria, Octocorallia): Different demographic features may lead affected populations to unexpected recovery and new equilibrium points[END_REF]. A reduction of two-thirds of the population density was reported in Portofino [START_REF] Cerrano | Gorgonian population recovery after a mass mortality event[END_REF], whereas in Punta Mesco, P. clavata was also heavily affected by mortality down to depths of 25 m [START_REF] Peirano | Decadal monitoring of corralligenous and bio-construction organisms in the eastern Ligurian sea (NW Mediterranean)[END_REF]. Although to our knowledge no literature has been published on mass mortality events of P. clavata in other locations within the Ligurian Sea, recreational divers reported observing dead gorgonian colonies in Calafuria in 1999 (L. De Nigris, pers. comm.). Moreover, other sessile invertebrates were affected by the increases in water temperature in the area of Bergeggi [START_REF] Parravicini | Consequences of seawater temperature anomalies on a Mediterranean submarine cave ecosystem[END_REF].

| Microsatellite genotyping

| Genetic diversity

Repeated multilocus genotypes were searched using GENODIVE 2.0b23 [START_REF] Meirmans | GENOTYPE and GENODIVE: two programs for the analysis of genetic diversity of asexual organisms[END_REF]. Observed and [START_REF] Nei | Molecular evolutionary genetics[END_REF] unbiased expected heterozygosity (H O and H E , respectively) were estimated using GENETIX 4.05 (Belkhir, Borsa, Chikhi, Raufaste, & Bonhomme, 2004). Allelic richness Q12 (A r ) was estimated using FSTAT 2.9.3.2 [START_REF] Goudet | FSTAT, a program to estimate and test gene diversities and fixation indices (version 2.9.3[END_REF], whereas private allelic Q13 richness (A p ) was calculated with a rarefaction procedure using HP-RARE software [START_REF] Kalinowski | HP-RARE 1.0: a computer program for performing rarefaction on measures of allelic richness[END_REF]. Linkage disequilibrium was tested for all loci using GENEPOP 4.2 [START_REF] Rousset | GENEPOP'007: a complete re-implementation of the GENEPOP software for Windows and Linux[END_REF], as implemented for online use (http://genepop.curtin.

edu.au Q14

). The test for Hardy-Weinberg equilibrium was computed for each locus using the null hypothesis of no heterozygote deficiency, and the level of significance was calculated with a Markov chain algorithm with default parameters (100 batches and 1000 iterations per batch). Single and multilocus F IS values were estimated using [START_REF] Weir | Estimating F-statistics for the analysis of population structure[END_REF] fixation index in GENEPOP. The presence of null alleles was estimated using the expectation maximization (EM) algorithm [START_REF] Dempster | Maximum likelihood from incomplete data via EM algorithm[END_REF] in FREENA [START_REF] Chapuis | Microsatellite null alleles and the estimation of population differentiation[END_REF].

Evidence of recent changes in effective population size as a result of the mass mortality events that affected P. clavata populations in the region during the past decade [START_REF] Cupido | Sexual structure of a highly reproductive, recovering gorgonian population: Quantifying reproductive output[END_REF] were tested using BOTTLENECK 1.2.02, following the two-phase model (TPM; [START_REF] Piry | Bottleneck: A program for detecting recent effective population size reductions from allele data frequencies[END_REF], with 95% single step mutation and a variance among multiple steps of approximately 12. A Wilcoxon's test was used, as suggested by [START_REF] Piry | Bottleneck: A program for detecting recent effective population size reductions from allele data frequencies[END_REF], for a data set of less than 20 loci. P values were corrected for multiple comparisons using false discovery rate (FDR) corrections [START_REF] Benjamini | Controlling the false discovery rate: A practical and powerful approach to multiple testing[END_REF].

| Population genetic structure

Population genetic structure was estimated using both the global and the pairwise F ST estimates of [START_REF] Weir | Estimating F-statistics for the analysis of population structure[END_REF] with GENEPOP, and following the excluding null allele (ENA) method in FREENA. The significance of pairwise genetic differentiation among sites was computed with an exact test using the default parameters in GENEPOP (100 batches and 1000 iterations per batch).

For linkage disequilibrium and F ST estimates, FDR corrections [START_REF] Benjamini | Controlling the false discovery rate: A practical and powerful approach to multiple testing[END_REF] were applied in order to adjust significance levels for multiple comparisons.

Isolation by distance (IBD) was tested at a regional scale through correlation between pairwise F ST /(1 -F ST ) values and the logarithm of the geographical distance between sites, with a Mantel test [START_REF] Mantel | The detection of disease clustering and a generalized regression approach[END_REF], using the ISOLDE program implemented in GENEPOP. Geographical distances were measured in GOOGLE EARTH by using set points (determined as the centre of the area of distribution where the samples were collected), and measuring both Euclidean and along-the-coast distances. All of the individuals sampled from a given site shared the same geographical coordinates.

The significance level of correlation between genetic differentiation and geographical distance among samples was computed using 10 000 permutations.

To estimate the extent of gene diffusion, multilocus genotypes were autocorrelated within spatial distance classes using the kinship coefficient of Loiselle Q15 , as implemented in SPAGEDI 1.4. This revealed the distance at which the random effects of genetic drift, not gene flow, are the primary determinants of genetic composition [START_REF] Hardy | SPAGEDI: a versatile computer program to analyse spatial genetic structure at the individual or population levels[END_REF]. Significant positive coefficients suggest that colonies within a specific distance class are more genetically related than colonies randomly taken from any distance class. The first x-intercept of the autocorrelogram, which corresponds to the average distance at which the genetic similarity between any two quadrats drops below 0.0, is termed the 'genetic patch size' [START_REF] Sokal | A test of spatial autocorrelation analysis using an isolation-by-distance model[END_REF]. Spatial distance intervals for kinship comparisons were constructed in such a way that the number of pairwise comparisons within each distance interval was approximately constant (between 3000 and 8000 pairs among distance classes). A total of 14 distance classes were considered: class 1, 0-14 km; class 2, 14-26 km; class To evaluate the number of clusters (K) in the data set without prior information regarding the geographic location of the colonies, a Bayesian method implemented in STRUCTURE 2.3.4 was used [START_REF] Falush | Inference of population structure using multilocus genotype data: linked loci and correlated allele frequencies[END_REF], 2007;[START_REF] Pritchard | Inference of population structure using multilocus genotype data[END_REF]. STRUCTURE was run under the admixture model, choosing the assumption of correlated allele frequencies and the option of recessive alleles to cope with null alleles, as suggested by [START_REF] Falush | Inference of population structure using multilocus genotype data: dominant markers and null alleles[END_REF]. The mean and variance values of log likelihoods were estimated for the number of clusters, from K = 1 to K = 12. Ten runs were performed for each K with 500 000 iterations and a burn-in period of 50 000. In order to identify the number of clusters that best fit the data, the resulting output files were then analysed based on the rate of change in the log probability of data between successive K values, using the Evanno method in STRUCTURE HARVESTER 0.6.94, as implemented online (http://tay-lor0.biology.ucla.edu/structureHarvester/ Q16 \/) [START_REF] Earl | STRUCTURE HARVESTER: a website and program for visualizing STRUCTURE output and implementing the Evanno method[END_REF] Q17 ). For the chosen K, the results were averaged using CLUMPP 1.1.2 [START_REF] Jakobsson | CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure[END_REF], and the graphical representation was performed in DISTRUCT 1.1 [START_REF] Rosenberg | DISTRUCT: A program for the graphical display of population structure[END_REF]. A hierarchical analysis of molecular variance (AMOVA; 1000 permutations) was performed in ARLEQUIN 3.5 [START_REF] Excoffier | Arlequin (version 3.0): An integrated software package for population genetics data analysis[END_REF] using the clusters (K = 6) defined by STRUCTURE.

| Migration patterns

First-generation migrants and long-term migration rates were estimated as proxies of genetic connectivity among the clusters identified by STRUCTURE. First-generation migrants were detected using the Bayesian computation of [START_REF] Rannala | Detecting immigration by using multilocus genotypes[END_REF] (n = 10000; α = 0.05) in GENECLASS 2.0 [START_REF] Piry | GENECLASS2: A Software for Genetic Assignment and First-Generation Migrant Detection[END_REF]. Considering the possibility of long dispersal via stepping stones, the level and direction of gene flow over the last several generations among clusters was estimated by performing a maximum-likelihood based analysis. MIGRATE-N 3.6 was used to estimate theta (Θ) and the mutation-scaled migration rate M based on microsatellite data [START_REF] Beerli | Maximum likelihood estimation of a migration matrix and effective population sizes in n subpopulations by using a coalescent approach[END_REF]. Theta is where N e is the effective population size and μ is the mutation rate per generation and site. M is defined as m/μ, representing the relative contribution of immigration and mutation to the variability brought into the population, whereas m is the fraction of new immigrants found in the population per generation. This calculation Q19 was run five times following the recommendations of [START_REF] Beerli | How to use migrate or why are markov chain monte carlo programs difficult to use?[END_REF], and the results were averaged.

3 | RESULTS

| Genetic diversity

2).

No significant linkage disequilibrium was found among loci (P > 0.05 after FDR correction). Estimates of null allele frequencies varied between 0 for Parcla09, Parcla10, Parcla12, Parcla14, and Par_d, to 0.26 for Parcla12.

Multilocus estimates of F IS ranged between 0.04 (BER1) and 0.33 (CAL1 and CER). All values were significantly different from zero (P < 0.001, mean = 0.20 ± 0.08), with the exception of BER1 (Table 2).

All populations showed significant levels of heterozygosity deficiency; however, significant deviations from Hardy-Weinberg were not found for all loci in all populations. The high frequency of null alleles found for locus Parcla12 in some populations did not explain the levels of heterozygote deficiency, given that heterozygote deficiency was still significant when Parcla12 was excluded from the analysis.

Despite these populations being considered as unaffected by recent mass mortality events, recent changes in effective population size resulting from an excess of the expected heterozygosity (a signal of a population bottleneck) were revealed at almost all locations (BER2, POF2, SES1, POV1, and CER; P < 0.05), with the exception of Punta Mesco and Calafuria (Table S2).

| Population genetic structure

The global multilocus F ST indicated strong genetic differentiation of P.

clavata populations between sites (F ST = 0.154; P = 0.0001). Pairwise F ST values ranged from 0 to 0.181, and almost all were significantly different after FDR correction (Table T3 3). Pairwise Q20 F ST values for sites located less than 50 m apart (BER1 versus BER2, SES1 versus SES2, POV1 versus POV2, and CAL1 versus CAL2), and two sites separated by approximately 20 km (PME versus POV2), but all within a 10-m depth range, did not significantly differ. No significant differences were found between pairwise F ST values and pairwise F ST values corrected for null alleles (Student's t-test, P = 0.163), which suggests that the effect of null alleles was negligible.

A significant but weak pattern of isolation by distance (IBD) was detected for both measures of geographic distance: straight line (R 2 = 0.23; P < 0.01) and along the coast (R 2 = 0.21; P < 0.01) (Figure S1).

Six genetic clusters with different levels of admixture were identified by the Bayesian clustering analysis (K = 6, ∆K = 12.26; Figure S1). The kinship coefficient of Loiselle was positive and significant for the 0-14 km (P < 0.01, 95% CI = -0.0002 to 0.0002) and 14-26 km (P < 0.01, 95% CI = -0.002 to 0.008) distance classes

(Figure F3
3). The genetic patch size was estimated to be 26 km.

| Patterns of recent and long-term migration

Overall, migration was bidirectional among all locations. Recent migration revealed by the analysis of first-generation migrants showed that the clusters of PMP and POF were the major sources of emigrants (40.62 and 28.12%, respectively). On the other hand, the clusters of BER and CAL were identified as the smallest source of emigrants, with 3.12% in both cases (Figure F4 4). In terms of immigrants, the clusters CAL and PMP received the highest proportion (15.9 and 14.81%, respectively), with the latter receiving immigrants for almost all other clusters (except BER). On the contrary, the cluster BER received the lowest number of immigrants (3.33%), all of them coming from the neighbouring cluster POF.

Patterns of long-term migration resulting from the maximum likelihood-based analysis using MIGRATE-N revealed that values of M (M = m/μ; migration/mutation) were > 1 for all locations, meaning that the effect of migration (m) is larger than the effect of mutation (μ) (Table S3). Migration rates among clusters identified by STRUCTURE decreased with increasing distance between locations. Nonetheless, evidence of migration among the most distant locations (>250 km) was found (Table T4 4). The cluster that received most immigrants over multiple generations was PMP (mean M = 18.02 ± 8.17), whereas the CAL cluster was identified as the one receiving the fewest immigrants (mean M = 7.78 ± 4.76). In terms of source of emigrants to all other populations, the SES cluster had the highest value (mean M = 18.68 ± 9.21), and the CER cluster had the lowest (mean Even though migration over such a large distance (>250 km) may be a rare occurrence, the level of migration (>10%) between clusters more than 60 km apart (PMP towards CAL and CER clusters), compared with the relative exchange of migrants among neighbouring clusters (<8%; e.g. POF and SES, 20 km apart), suggests that dispersal over 60 km may be more common than previously thought;

however, the lowest level of migration (3.33%) was also found in one cluster (BER) located less than 60 km away from its closest neighbour (POF).

| DISCUSSION

Genetic evidence of population bottlenecks was found in most of the populations of P. clavata studied in the Ligurian Sea. Contemporary and long-term patterns of migration showed that all of the studied P.

clavata populations are connected by larval dispersal.

| Genetic diversity and structure

The levels of genetic diversity among different P. clavata populations dwelling between 25 and 40 m of depth in the Ligurian Sea were similar to those found among shallow and deep populations in the same region, and in the north-west Mediterranean [START_REF] Arizmendi-Mejía | Combining genetic and demographic data for the conservation of a Mediterranean marine habitatforming Species[END_REF][START_REF] Mokhtar-Jamaï | From global to local genetic structuring in the red gorgonian Paramuricea clavata: The interplay between oceanographic conditions and limited larval dispersal[END_REF][START_REF] Pilczynska | Genetic diversity and local connectivity in the Mediterranean red gorgonian coral after mass mortality events[END_REF]. Departures from Hardy-Weinberg equilibrium, as well as high F IS values and evidence of heterozygote deficit have been previously reported in P.

clavata populations [START_REF] Mokhtar-Jamaï | From global to local genetic structuring in the red gorgonian Paramuricea clavata: The interplay between oceanographic conditions and limited larval dispersal[END_REF][START_REF] Pérez-Portela | Lonely populations in the deep: Genetic structure of red gorgonians at the heads of submarine canyons in the north-western Mediterranean Sea[END_REF][START_REF] Pilczynska | Genetic diversity and local connectivity in the Mediterranean red gorgonian coral after mass mortality events[END_REF]. This is a common occurrence in gorgonian species, mainly attributed to the low dispersal capabilities of the larvae [START_REF] Gutiérrez-Rodríguez | Microsatellite variation reveals high levels of genetic variability and population structure in the gorgonian coral Pseudopterogorgia elisabethae across the Bahamas[END_REF][START_REF] Ledoux | Fine-scale genetic structure and inferences on population biology in the threatened Mediterranean red coral, Corallium rubrum[END_REF].

The global F ST value (F ST = 0.154) was similar to what was previously reported for the species in the Ligurian Sea [START_REF] Pilczynska | Genetic diversity and local connectivity in the Mediterranean red gorgonian coral after mass mortality events[END_REF] Hence, the absence of a genetic bottleneck in Calafuria and Punta Mesco is striking.

Although very unlikely, it is possible that those populations were indeed not affected by mass mortality events: in Punta Mesco, the affected colonies were only reported down to 25 m of depth [START_REF] Peirano | Decadal monitoring of corralligenous and bio-construction organisms in the eastern Ligurian sea (NW Mediterranean)[END_REF], whereas mass mortalities have not been reported in Calafuria. Nonetheless, it is worth considering that high recruitment, either from reproduction of local survivors or through migration, could have cancelled the signature of a genetic bottleneck. Indeed, immigration rates in Calafuria and Punta Mesco were the highest in the region.

Nonetheless, it is also worth mentioning the potential bias resulting from the methodological limitations of reconstructing the demographic history of populations from genetic data, which could have led to underestimating the number of populations that underwent a drastic decline in population size [START_REF] Girod | Inferring Population Decline and Expansion From Microsatellite Data: A simulation-based evaluation of the Msvar method[END_REF][START_REF] Peery | Reliability of genetic bottleneck tests for detecting recent population declines[END_REF].

A closer look at the Punta Mesco and Portovenere cluster (PMP)

showed that mass mortality events of similar intensity were recorded in the two sites of Portovenere over 10 years ago. Also, both populations had started recovering a few years after that, but at different rates [START_REF] Cupido | Sexual structure of a highly reproductive, recovering gorgonian population: Quantifying reproductive output[END_REF]. Results of the present study showed higher immigration in one of the two sites within Portovenere, with the site with lower immigration showing a signature of genetic bottleneck.

Although demographic bottlenecks may be recurrent in Portovenere, migration among coastal populations may help maintain local genetic diversity and offset population collapse, making the signature of a bottleneck undetectable [START_REF] Aguilar | The distribution of nuclear genetic variation and historical demography of sea otters[END_REF][START_REF] Spong | A near-extinction event in lynx: Do microsatellite data tell the tale?[END_REF]. Similarly, in Portofino, a mass mortality event was documented in 1999, and even though the mass mortality event dramatically reduced the reproductive potential of the native population, after 3 years gorgonian colony densities were greater than before the mortality event [START_REF] Cerrano | Gorgonian population recovery after a mass mortality event[END_REF]. The estimated age of first reproduction of P. clavata is 3 years. Therefore, enhanced recruitment following mortality should be attributed to immigration. The large number of immigrants in Portofino supports this hypothesis.

Despite no reports of P. clavata mass mortality events in Bergeggi, [START_REF] Parravicini | Consequences of seawater temperature anomalies on a Mediterranean submarine cave ecosystem[END_REF] reported heatwaves affecting other sessile invertebrates. The genetic bottleneck signature observed in one of the Bergeggi sites suggests that mass mortality events could have affected P. clavata populations in this area before 1999, and its permanence highlights the vulnerability of populations with low immigration rates.

Another potential explanation for the observed deficiency of heterozygotes among the studied populations could be the mixing of different gene pools, which would cause a Wahlund effect. This hypothesis has been proposed already for P. clavata populations in the Mediterranean Sea [START_REF] Mokhtar-Jamaï | Role of evolutionary and ecological factors in the reproductive success and the spatial genetic structure of the temperate gorgonian Paramuricea clavata[END_REF], and for the Mediterranean red coral, Corallium rubrum [START_REF] Costantini | Fine-scale genetic structuring in Corallium rubrum: evidence of inbreeding and limited effective larval dispersal[END_REF]. These findings, combined with the weak correlation found between genetic and geographic distance, demonstrate, as previously proposed by [START_REF] Pérez-Portela | Lonely populations in the deep: Genetic structure of red gorgonians at the heads of submarine canyons in the north-western Mediterranean Sea[END_REF], that other factors rather than geographical distance alone may influence the pattern of genetic structure among P. clavata populations. Effective migration, in fact, is the result of the filtering of dispersal potential related to larval traits (e.g. pelagic larval duration, larval buoyancy, and behaviour) by environmental pre-and post-settlement drivers, such as hydrodynamics, predation pressure, and local geomorphological features [START_REF] Lortie | Rethinking plant community theory[END_REF].

| Patterns of migration

According to these results, the most important sources of migrants were the sites located in the central part of the Ligurian Sea, which is consistent with the large-scale separation of currents that occurs in front of Portovenere Q23 [START_REF] Pinardi | Variability of the large scale general circulation of the Mediterranean Sea from observations and modelling: A review[END_REF]. This flow separation results in two hydrodynamic provinces [START_REF] Rossi | Hydrodynamic provinces and oceanic connectivity from a transport network help designing marine reserves[END_REF]. Although the Ligurian current bordering the narrow continental shelf induces a rapid directional drift along the coast north-west of Portovenere, the detachment Q24 of meso-scale gyres south-east of Portovenere promotes higher regional retention with diffuse dispersal in the Gulf of La Spezia and north Tyrrhenian Sea. The hydrodynamic province identified north-west of Portovenere extends along the French coast to Marseille, suggesting that migrants proceeding from the northwestern Ligurian Sea (e.g. Bergeggi) should be found in that location. This hydrodynamic pattern also supports the observation that the population of Bergeggi receives five times less immigrants than the population of Calafuria, and that a higher level of gene flow occurred towards the eastern part of the Ligurian Sea. It is also worth noting that although connectivity among coastal P.

clavata populations provided an important contribution to the recovery of populations affected by mass mortality events, it may not be the only explanation for their upturn, and migration from deeper populations should not be excluded [START_REF] Pilczynska | Genetic diversity and local connectivity in the Mediterranean red gorgonian coral after mass mortality events[END_REF]. In any case, the present study shows that an effective number of migrants per generation of less than 5% was not sufficient for some populations (Bergeggi and Sestri Levante) to compensate their genetic divergence resulting from the combined effects of mutation and genetic drift.

| Implications for conservation in the context of global change

According to climate change projections, the rising of ocean temperatures and increasing frequency of extreme heatwave events (de Madron et al., 2011 Q25

; [START_REF] Meehl | Trends in extreme weather and climate events: Issues related to modeling extremes in projections of future climate change[END_REF], could lead to more frequent and more severe demographic breakdowns, especially in small and isolated populations. Paramuricea clavata is an extremely vulnerable species to extreme climatic conditions, which may cause severe depletions. It is generally assumed that heatwaves affect mainly shallow populations; however, during the heatwave of September 1999, unusual temperatures were recorded down to at least 32 m in depth [START_REF] Bramanti | Recruitment, early survival and growth of the Mediterranean red coral Corallium rubrum (L 1758), a 4-year study[END_REF], and genetic data demonstrate that populations of P. clavata between 25 and 40 m in depth have also been affected by drastic reductions in population size, highlighting that extreme climatic events impact coastal populations down to at least 40 m in depth.

Nonetheless, this study provides evidence that P. clavata larvae are able to move among populations up to a distance of tens of kilometres, promoting population resilience and recovery from disturbance through the arrival of propagules from other locations.

Connectivity has become an increasingly important element in spatial conservation planning [START_REF] Olds | Quantifying the conservation value of seascape connectivity: A global synthesis[END_REF]; however, only few studies have actually provided empirical data on the role of connectivity in population recovery following disturbance. Larval connectivity is remarkably relevant in sessile marine invertebrates, which Q26 rely on presettlement processes for dispersal. The protection of populations identified as sources of migrants would help the regional persistence of the species. The identification of key source populations may be of great value for the implementation of conservation measures and regional spatial planning in the context of increasing global climaterelated disturbance. Nonetheless, seascape genetic inference can be strongly influenced by the recent demographic history of populations.

Demographic stability can regulate gene flow, when recruitment in benthic populations is limited by a lack of space [START_REF] Padrón | Modelling the effect of demographic traits and connectivity on the genetic structuration of marine metapopulations of sedentary benthic invertebrates[END_REF]. Therefore, this study highlights the importance of integrating demographic and genetic data (see Arizmendi-Mejía et al., 2015 Q27 ).

In the Ligurian Sea, several Marine Protected Areas have been

  Marine biodiversity on a global scale is under growing pressure from environmental alteration caused by human activities, including climate change

  clavata (Mokhtar-Jamaï et al., 2011 Q9 ). Nevertheless, in the last couple of years genetic studies have revealed a mixed pattern of isolation by distance and clustering on a regional scale, with connectivity occurring among neighbouring, fragmented populations Q10 (Arizmendi-Mejía

2

  | MATERIAL AND METHODS 2.1 | Sampling Colonies of the red gorgonian P. clavata were sampled by scuba diving at seven locations (from 25 to 40 m in depth) in the Ligurian Sea in May 2014 (Figure F1 1). All the samples were considered as coastal populations, regardless of depth. Small (<3 cm) branch fragments from 26 were collected at two different sites per location (separated by distances of between 30 m and 1 km), except in Punta Mesco and Cerboli Island. At those two locations, colonies were sampled only at one site because of restricted sampling permission by the Marine Protected Area authorities (Table

  Total genomic DNA was extracted from 372 colonies using the DNA extraction kit NucleoSpin Tissue by Macherey-Nagel. All individuals were genotyped at six polymorphic microsatellite loci developed for P. clavata, as used in Mokhtar-Jamaï et al. (2011): Parcla09, Parcla10, Parcla12, Parcla14, Parcla17 (Molecular Ecology Resources Primer Development Consortium et al., 2010), and Par_d[START_REF] Agell | Isolation and characterization of eight polymorphic microsatellite loci for the Mediterranean gorgonian Paramuricea clavata[END_REF]. Genotyping was carried out on an ABI 3130 Genetic Analyser using an internal size standard (GeneScan 600 LIZ; Applied Biosystems, Foster City CA, USA Q11 ), and allele scoring was performed using PEAK SCANNER SOFTWARE 1.0 (Applied Biosystems).

FIGURE 1

 1 FIGURE 1 Map of the seven sampling locations of Paramuricea clavata in the Ligurian Sea. Black lines represent the simulated surface currents averaged during summer (the spawning season for P. clavata in the Mediterranean Sea;[START_REF] Linares | Size distribution, density and disturbance in two Mediterranean gorgonians: Paramuricea clavata and Eunicella singularis[END_REF], according to[START_REF] Pinardi | Variability of the large scale general circulation of the Mediterranean Sea from observations and modelling: A review[END_REF] and[START_REF] Rossi | Hydrodynamic provinces and oceanic connectivity from a transport network help designing marine reserves[END_REF]

  e μ Q18 for a diploid system with nuclear microsatellite loci,

  2 to K = 5, see Figure S2): (i) BER1 and BER2 comprised the Bergeggi cluster (BER), with a coefficient of population membership of 90%; (ii) POF1 and POF2 formed the Portofino cluster (POF), with a coefficient of population membership of 66%; (iii) SES1 and SES2 comprised the Sestri Levante cluster, with a coefficient of population membership of 77% (SES); (iv) PME, POV1, and POV2 were grouped in the Portovenere and Punta Mesco cluster (PMP), with a coefficient of population membership of 60%; (v) CAL1 Q21 and CAL2 formed the Calafuria cluster (CAL), with a coefficient of population membership of 70%; and (vi) CER made up the Cerboli cluster (CER), with a coefficient of population membership of 80%. The results of the AMOVA confirmed the clustering structure (P < 0.001; Table

FIGURE 2 FIGURE 3 M

 23 FIGURE 2 Population genetic structure of Paramuricea clavata in the Ligurian Sea, as revealed by STRUCTURE. Each individual is represented by a vertical line, where different colours indicate the individual proportion of membership to each cluster (K = 6). Key: Bergeggi (BER), purple; Portofino (POF), green; Sestri Levante (SES), fuchsia; Punta Mesco and Portovenere (PMP), orange; Calafuria (CAL), blue; and Cerboli (CER), pink

FIGURE 4

 4 FIGURE 4 Pattern of first-generation migration among Paramuricea clavata populations in the Ligurian Sea. The pie charts represent recruitment in each cluster. The white portion indicates the propotion of self-recruits, whereas the colours represent immigrants from other locations. The bar plot indicates the contribution of emigrants (%) of each cluster: Bergeggi (BER), purple; Portofino (POF), green; Sestri Levante (SES), fuchsia; Punta Mesco and Portovenere (PMP), orange; Calafuria (CAL), blue; and Cerboli (CER), pink

F

  ST = 0.035;[START_REF] Arizmendi-Mejía | Combining genetic and demographic data for the conservation of a Mediterranean marine habitatforming Species[END_REF] were much lower than those found in this study. Nonetheless, non-significant F ST values were found among most pairs of sites separated by 50-100 m, and generally at similar depths.4.2 | Evidence of population bottleneckIn contrast with with Pilczynska et al. (2016), who found no evidence of a genetic bottleneck in two populations at 20-31 m of depth, a signal of a recent and drastic reduction of P. clavata effective population size was found in most of the populations studied (at depths of 25-40 m), with the exception of Punta Mesco and Calafuria. The findings corroborate previous studies reporting mass mortality events in populations dwelling at depths deeper than 20 m (in Sardinia, at depths of 20-40 m, Huete-Stauffer et al., 2011; in the north-western Mediterranean, at depths of 50-70 m, Pérez-Portela et al., 2016).

(

  2011), who found that despite strong genetic structure, some gene flow could still occur over short and longer distances (up to 160 km) across the Mediterranean Sea. Unlike the study of Arizmendi-Mejía et al. (2015), connectivity among populations of P. clavata in the Ligurian Sea was found between non-neighbouring clusters, but was not homogeneus throughout the region.

  established based on the natural value of the habitats and/or administrative requirements. This study suggests that the spatial arrangement of the Ligurian Sea network of MPAs is effective in maintaining connectivity among P. clavata populations in the study area. Moreover, the distance between MPAs is in the range of the dispersal distance of the species. Therefore, in case of future disturbances, the network of MPAs in the Ligurian Sea would support the recovery of the populations. Nonetheless, it remains necessary to test the deep refugia hypothesis, given that in the event of a disturbance affecting all shallow populations, deep populations would become of utmost importance. Taking this into consideration, it would be thus fundamental to promote the creation of offshore MPAs that could be incorporated into the existing MPA network, and ensure connectivity among shallow and deep populations. Furthermore, the eastward direction of gene flow in the study area suggests that connectivity should be further investigated in populations located west of Bergeggi, to provide a pattern of the connectivity across national borders and foster the creation of transnational MPA networks.Belkhir, K., Borsa, P., Chikhi, L., Raufaste, N., & Bonhomme, F. (2004).GENETIX 4.05, Logiciel sous Windows™ pour la génétique des populations. Montpellier, France: Laboratoie, Génome, Populations, Interactions, CNRS UMR 5000, Université de Montpellier II. .,Coma, R., Díaz, D., Zabala, M., Hereu, B., & Dantart, L. (2005).Immediate and delayed effects of a mass mortality event on gorgonian population dynamics and benthic community structure in the NW Mediterranean Sea. Marine Ecology ProgressSeries, 305, 127-137. ., Russell, B., Thomsen, M., & Connell, S. (2014). Marine biodiversity and climate change. In B. Freedman (Ed.), Global Environmental Change (pp. 181-187). Dordrecht, the Netherlands: Springer.SUPPORTING INFORMATIONAdditional supporting information may be found online in the Supporting Information section at the end of the article. How to cite this article: Padrón M, Costantini F, Bramanti L, Guizien K, Abbiati M. Genetic connectivity supports recovery of gorgonian populations affected by climate change. Aquatic Conserv: Mar Freshw Ecosyst. 2018;1-12. https://doi.org

TABLE 1

 1 Sampling locations of Paramuricea clavata in the Ligurian Sea

	Location	Marine Protected Area	Site	Site ID	Latitude	Longitude	Depth (m)	Samples
	Bergeggi	Area Marina Protetta Isola di Bergeggi	Bergeggi 1	BER1	44.2416	8.4461	37	32
			Bergeggi 2	BER2	44.2416	8.4461	38	30
	Portofino	Area Marina Protetta di Portofino	Punta del Faro	POF1	44.1753	9.1308	36	31
			Isuela	POF2	44.1903	9.1851	40	31
	Sestri Levante	No	Punta Baffe 1	SES1	44.2482	9.3963	30	32
			Punta Baffe 2	SES2	44.2482	9.3963	30	32
	Punta Mesco	Area Marina Protetta delle Cinque Terre	Punta Mesco	PME	44.1247	9.6322	29	26
	Portovenere	Parco Naturale Regionale di Porto Venere	Secca Tinetto	POV1	44.0122	9.5104	25	30
		-Area di Tutela Marina	Parete Grotta	POV2	44.0126	9.5102	25	32
	Calafuria	No	Gorgonie Rose 1	CAL1	43.2812	10.2002	35	32
			Gorgonie Rose 2	CAL2	43.2812	10.2002	25	29
	Cerboli Island	No	Cerboli Island	CER	42.8555	10.5469	30	37

TABLE 2

 2 Parameters of genetic diversity for 12 sites of Paramuricea clavata at six microsatellite loci

	Location	Site	H O	H E	A r	A p	F IS
	Bergeggi	BER1	0.65	0.67	4.96	0.13	0.04
		BER2	0.63	0.73	5.42	0.2	0.17
	Portofino	POF1	0.61	0.66	5.24	0.55	0.22
		POF2	0.64	0.76	6.29	0.5	0.21
	Sestri Levante	SES1	0.59	0.68	6.01	0.34	0.14
		SES2	0.71	0.70	6.34	0.3	0.14
	Punta Mesco	PME	0.56	0.70	6.41	0.51	0.20
	Portovenere	POV1	0.55	0.69	6.32	0.42	0.21
		POV2	0.58	0.69	6.55	0.28	0.10
	Calafuria	CAL1	0.52	0.74	6.25	0.55	0.33
		CAL2	0.54	0.68	5.54	0.42	0.27
	Cerboli	CER	0.49	0.71	6.19	0.62	0.33
		Mean value	0.59 ± 0.06	0.70 ± 0.03	5.96 ± 0.52	0.40 ± 0.15	0.20 ± 0.08
	Abbreviations: H

O , observed heterozygosity; H E , unbiased expected heterozygosity; A r , allelic richness; A p , private allelic richness; F IS ,

[START_REF] Weir | Estimating F-statistics for the analysis of population structure[END_REF]

estimate of

Wright's (1951)

fixation index, with significant values in bold (0.05 threshold after Bonferroni correction).

TABLE 3

 3 Pairwise F ST values among Paramuricea clavata sites in the Ligurian Sea Significant values, after false discovery rate (FDR) correction, are highlighted in bold (P < 0.05).

		BER1	BER2	POF1	POF2	SES1	SES2	PME	POV1	POV2	CAL1	CAL2	CER
	BER1	0
	BER2	0.004	0
	POF1	0.116	0.116	0
	POF2	0.132	0.144	0.066	0
	SES1	0.127	0.131	0.089	0.125	0
	SES2	0.099	0.096	0.069	0.108	0.011	0
	PME	0.158	0.149	0.097	0.143	0.083	0.066	0
	POV1	0.141	0.139	0.095	0.136	0.088	0.092	0.061	0
	POV2	0.125	0.116	0.094	0.122	0.086	0.074	0.037	0.029	0
	CAL1	0.189	0.181	0.139	0.177	0.111	0.081	0.143	0.177	0.164	0
	CAL2	0.112	0.117	0.089	0.139	0.075	0.052	0.103	0.121	0.112	0.032	0
	CER	0.143	0.145	0.074	0.138	0.091	0.079	0.071	0.062	0.094	0.154	0.102	0

TABLE 4

 4

			Destination
			BER	POF	SES	PMP	CAL	CER	Mean
		BER	0	1.60	16.00	13.18	6.81	7.56	9.03 ± 5.65
	Source	POF	1.81	0	14.15	21.48	4.04	10.99	10.49 ± 7.92
		SES	20.40	26.69	0	26.57	4.63	15.13	18.68 ± 9.21
		PMP	16.62	13.30	17.50	0	15.89	27.70	18.20 ± 5.53
		CAL	5.17	5.47	13.79	22.59	0	3.14	10.03 ± 8.12
		CER	5.49	1.62	9.84	6.29	7.50	0	6.15 ± 3.02
		Mean	9.89 ± 8.10	9.73 ± 10.61	14.25 ± 2.88	18.02 ± 8.17	7.78 ± 4.76	12.90 ± 9.37

Long-term migration rates M (= m/μ; migration/mutation) among genetic clusters of Paramuricea clavata in the Ligurian Sea, as revealed by the maximum likelihood-based analysis using MIGRATE-N (M) PADRÓN ET AL.

ACKNOWLEDGMENTS

This work was (co-)funded through a MARES grant. MARES is a joint doctorate programme selected under Erasmus Mundus coordinated by Ghent University (FPA 2011-0016) (for extra information, see http://www.mares-eu.org Q28 ), and by the PRIN 2011 project 'Coastal bioconstructions: structure, function and management' and the PRIN 2015 (MIUR) project 'Resistance and resilience of Adriatic mesophotic biogenic habitats to human and climate change threats' (2015J922E4). The authors are grateful to the administration of the Portofino and Arcipelago Toscano Marine Protected Areas for the authorization to collect valuable samples of Paramuricea clavata. We are grateful to Silvia Cocito et al. for providing the samples from Punta Mesco. We also thank D. Aurelle for constructive comments. We give particular thanks to L. De Nigris, T. Garcia, L. Castellano, and Gruppo Archeosub Labronico for invaluable help in underwater sampling and logistics.

Q3 Q39 MARES Doctoral Programme, Grant/Award Number: FPA 2011-0016; PRIN 2011 project 'Coastal bioconstructions: structure, function and management'; PRIN 2011 project, Grant/ Award Number: 2015J922E4; MARES is a joint doctorate programme, Grant/Award Number: FPA 2011-0016