Genetic connectivity supports recovery of gorgonian populations affected by climate change

Climate‐induced threats are increasingly affecting marine populations worldwide. In the last few decades, several gorgonian species have been affected by mass mortality events in the north‐west Mediterranean, putatively linked to local sea temperature increases during heatwaves. For many benthic sessile species, recovery after disturbances depends upon larval supply shaping the connections among populations. In the Ligurian Sea, genetic analyses showed that some Paramuricea clavata populations recovered after mass mortality events; however, the patterns of connectivity, and the potential role of migration in supporting the recovery of the populations of P. clavata, across its distribution range within the Ligurian Sea are still unknown. In this study, the population genetic structure and migration patterns of P. clavata populations have been analysed across seven sites in the Ligurian Sea, some of which have been affected by mass mortality events. Evidence of a population bottleneck was found in most of the populations studied. Significant genetic differentiation was found among P. clavata populations, reflecting habitat fragmentation at a regional scale, except for two populations found 20 km apart. Continuing gene flow between distant populations was also revealed. Empirical data suggest that gene flow among populations may have contributed to support their recovery from mass mortality events. The study identified populations in the central part of the Ligurian Sea that can be strategic for the regional persistence of the species. These findings highlight that the preservation of key populations could maintain connectivity and gene flow in the metapopulation, and increase the resilience of the species.

. 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 (Doney et al., 2012). 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 (Lacy, 1987). Furthermore, within this context, genetic drift may amplify the loss of genetic variation in isolated populations, and potentially lead to inbreeding (Nei, Maruyama, & Chakraborty, 1975). Nonetheless, the loss of genetic diversity caused by a population collapse can be offset or reversed by migration processes (Crow & Kimura, 1970;Keller et al., 2001). 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 (Almany et al., 2009). 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 (Lipcius et al., 2008;Shanks, Grantham, & Carr, 2003), 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 (Bramanti, Iannelli, Fan, & Edmunds, 2015;Caplins et al., 2014;Gilmour, Smith, Heyward, Baird, & Pratchett, 2013;Magris, Pressey, Weeks, & Ban, 2014;Underwood, Smith, van Oppen, & Gilmour, 2007;Vollmer & Palumbi, 2007). 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 (Hughes, Inouye, Johnson, Underwood, & Vellend, 2008;Kovach et al., 2015;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 (Waser & Strobeck, 1998).
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 (Cupido, Cocito, Sgorbini, Bordone, & Santangelo, 2008; but see Rivetti, Fraschetti, Lionello, Zambianchi, & Boero, 2014). Among the species affected, gorgonians (Cnidaria, Octocorallia) were the most damaged, with a local loss of more than 50% in density and biomass (Bramanti, Magagnini, De Maio, & Santangelo, 2005;Cerrano et al., 2000;Coma et al., 2009;Cupido et al., 2008;Garrabou et al., 2009;Santangelo et al., 2015). Gorgonians are characterized by long life span and slow growth, which along with the low fecundity of small colonies delays their recovery after disturbance (Linares, Coma, Garrabou, Díaz, & Zabala, 2008), even when negative density-dependent control of recruitment is removed (Santangelo et al., 2015). 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 (Ponti et al., 2014).
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 (Rivetti et al., 2014). The impact of reported mass mortality events on the demographic structure and dynamics of several populations has been assessed Santangelo et al., 2015;Santangelo, Bramanti, & Iannelli, 2007). The recovery of these populations has been attributed to local replenishment (Santangelo et al., 2015), overlooking the role of larval immigration from other source populations, based on the assumption of low larval dispersal in P.
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 et al., 2015;Pérez-Portela et al., 2016;Pilczynska, Cocito, Boavida, Serrão, & Queiroga, 2016). 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 (Pilczynska et al., 2016). 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 (Huete-Stauffer et al., 2011).
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.   (Table   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% (Cupido et al., 2008;Santangelo et al., 2015). A reduction of two-thirds of the population density was reported in Portofino (Cerrano et al., 2005), whereas in Punta Mesco, P. clavata was also heavily affected by mortality down to depths of 25 m (Peirano, Sgorbini, Cupido, Lombardi, & Cocito, 2009). Although to our knowledge no literature has been pub-    Linares et al., 2008), according to Pinardi and Masetti (2000) and Rossi et al. (2014)

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 Weir and Cockerham's (1984) fixation index in GENEPOP. The presence of null alleles was estimated using the expectation maximization (EM) algorithm (Dempster, Laird, & Rubin, 1977) in FREENA (Chapuis & Estoup, 2007).
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 (Cupido et al., 2012) were tested using BOTTLENECK 1.2.02, following the two-phase model (TPM; Piry, Luikart, & Cornuet, 1999), with 95% single step mutation and a variance among multiple steps of approximately 12. A Wilcoxon's test was used, as suggested by Piry et al. (1999), for a data set of less than 20 loci. P values were corrected for multiple comparisons using false discovery rate (FDR) corrections (Benjamini & Hochberg, 1995).

| Population genetic structure
Population genetic structure was estimated using both the global and the pairwise F ST estimates of Weir and Cockerham (1984) 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 (Benjamini & Hochberg, 1995) 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 (Mantel, 1967), 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 (Hardy & Vekemans, 2002). 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' (Sokal & Wartenberg, 1983). 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 (Falush, Stephens, & Pritchard, 2003Pritchard, Stephens, & Donnelly, 2000). 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 Falush et al. (2007). 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 \/) (Earl & von Holdt, 2011 Q17 ). For the chosen K, the results were averaged using CLUMPP 1.1.2 (Jakobsson & Rosenberg, 2007), and the graphical representation was performed in DISTRUCT 1.1 (Rosenberg, 2004). A hierarchical analysis of molecular variance (AMOVA; 1000 permutations) was performed in ARLEQUIN 3.5 (Excoffier, 2005) 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 Rannala and Mountain (1997) (n = 10000; α = 0.05) in GENECLASS 2.0 (Piry et al., 2004). 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 (Beerli & Felsenstein, 2001). 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 Beerli (2009) 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   F2 2; for K = 2 to K = 5, see Figure S2 Table S1).  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.    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.   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 (Cupido et al., 2012). 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 (Aguilar, Jessup, Estes, & Garza, 2008;Spong & Hellborg, 2002). 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 (Cerrano et al., 2005). 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, Parravicini et al. (2010) 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 (Mokhtar-Jamaï et al., 2013), and for the Mediterranean red coral, Corallium rubrum (Costantini, Fauvelot, & Abbiati, 2007).

| Patterns of migration
The analyses of first-generation migrants and long-term migration These findings, combined with the weak correlation found between genetic and geographic distance, demonstrate, as previously proposed by Pérez-Portela et al. (2016), 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 (Lortie et al., 2004).
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 (Pinardi & Masetti, 2000). This flow separation results in two hydrodynamic provinces (Rossi, Ser-Giacomi, López, & Hernandez-García, 2014). 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 (Pilczynska et al., 2016). 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 ;Meehl et al., 2000), 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 (Bramanti et al., 2005), 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 (Olds et al., 2016); 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 (Padrón & Guizien, 2016). 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