A total of 541 individuals from 20 Mediterranean localities were collected by SCUBA diving. Due to the close proximity of Bonassola and Framura (both in the Ligurian Sea), and the relatively low number of specimens sampled there, the two localities were analysed as a single population (BON). Four localities were sampled in the eastern Mediterranean: one each in Crete and Greece and two in Cyprus. Seven localities were sampled in the central Mediterranean in the Ionian Sea or along both coasts of the Adriatic Sea in Italy, Croatia and Montenegro. Nine localities were sampled in the western Mediterranean, in Spain and western Italy (Table 1; Figure 1).

Some of these localities were previously analysed (Casado-Amezúa et al., 2014; López-Márquez et al., 2019, 2021), however, with two different sets of primers (the first set consisted of eight primer pairs and was used by Casado-Amezúa et al., 2011 and López-Márquez et al., 2019; the second set, used by López-Márquez et al., 2021, consisted of 11 new pairs). López-Márquez et al. (2021) previously analysed three of the four eastern Mediterranean populations (NEA, KRY and LIO) using both primer sets. The results of these previous studies served as a reference and were included here in order to provide a complete picture of the species distribution across the entire Mediterranean. To facilitate comparisons, several specimens representing the analysed populations were also genotyped to ensure the uniform allele assignation with both primer sets. As the Adriatic (OTR, TOG, TRE, KOR and BOK) and Ionian (POC) populations were already analysed using the first primer set (López-Márquez et al., 2019), here we only analysed them using the second set. The newly sampled populations from the western Mediterranean (GAV, ESP, MEN, BON and PAL), the Adriatic (POR) and the eastern Mediterranean (CRE), together with the populations previously studied by Casado-Amezúa et al. (2014) (CDP, CON and PUN), were genotyped using both primer sets.

In order to avoid sampling clonal individuals, colonies collected by SCUBA diving were separated by at least 1 m. Some polyps from each colony were carefully excised and stored in vials with absolute ethanol and preserved at 4°C before processing. All necessary permits were obtained for the field studies.

The QIAGEN Biosprint 15 DNA Blood Kit (Qiagen) was used for DNA extraction and purification, following the manufacturer’s protocol. Each DNA sample was diluted to a final concentration of 0.3 ng/μl.

The DNA samples were genotyped with the 19 polymorphic microsatellite loci previously isolated for C. caespitosa [the eight by Casado-Amezúa et al. (2011) and the 11 by López-Márquez et al. (2021)]. The 19 primers pairs were combined in five multiplex reactions at a concentration of between 0.2 and 0.4 µM, and mixed with 1x Qiagen Multiplex PCR Master Mix, 0.3 ng of DNA and water to a total volume of 7 µl. To facilitate the genotyping, the forward primer from each pair was end-labelled with either NED, VIC, PET or 6-FAM, and the reverse primers were pig-tailed with 5′-GTTTCTT-3’ (López-Márquez et al., 2021). The PCR cycling profile included an initial denaturation step at 95°C for 5 min, followed by 35 cycles at 94°C for 30 s, 56°C for 90 s and 72°C for 30 s, and a final extension step at 72°C for 10 min.

GENEMAPPER software v4.0 (Applied Biosystems) was used to analyse the electropherograms. The presence of null alleles was checked with MICRO-CHECKER v.2.2.3 (Van Oosterhout et al., 2004).

By performing the option “multilocus matches” in GenAlEx 6.5 (Peakall and Smouse, 2006), we calculated the number of both unique (N _g ) and non-unique (identical) multilocus genotypes per site, which allowed us to determine the number of clones present in our data set. To analyse clonal structure parameters, and to calculate genotypic richness (standardised), genotypic evenness and genotypic diversity, we followed the methodology of Aranceta-Garza et al. (2012). Also, we classified populations according to the ratio of asexual to sexual reproduction they displayed, which was based on both genotypic evenness and genotypic diversity as in Baums et al. (2006) who delimited four groups: sexual, mostly sexual, mostly asexual and asexual. We also calculated the D index (Baums et al., 2006) as: D = 1 − (     ∑ n i ( n i − 1 )   N ( N − 1 ) ) Where n _i is the number of individuals of genotype i, and N the total number of individuals in the population. D values range between zero and one. If D equals zero, the entire population belong to the same unit of clonal growth.

After excluding the individuals identified as clones from the data set, we used GenAlEx 6.5 and GENEPOP v4.0 (Raymond and Rousset, 1995) to calculate allelic diversity (N _a ), observed (H _o ) and expected (H _e ) heterozygosity and the index F _IS , which is commonly used as an inbreeding coefficient, and test for Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium (LD). We corrected the significance of p values with the sequential Bonferroni method (Rice, 1989). Prior to performing further analyses with the data set, we look for the presence of loci under selection using two methodologies: a neutrality test performed in ARLEQUIN v3.5 (Excoffier et al., 2005) and a Bayesian approach using BAYESCAN (Foll, 2012). Finally, we estimated effective population size (N _e ) with NE ESTIMATOR 2 (Do et al., 2014), following the linkage disequilibrium method (Waples and Do, 2008).

To assess population structure and differentiation in C. caespitosa, we used various methods, including Wright’s fixation index (F _ST ), principal coordinates analysis (PCoA), isolation by distance (IBD), BARRIER analyses, the analysis of molecular variance (AMOVA), STRUCTURE and discriminant analysis of principal components (DAPC).

Population differentiation was estimated with F _ST between pairwise sampling sites through Weir and Cockerham’s estimators in GENETIX v.4.05 (Belkhir et al., 2004). The negative values obtained because of mathematical artefacts (F _ST cannot be negative) were set to zero. Standardised F _ST values (F’ _ST ) were calculated in GenAlEx, assuming that each population has different alleles for each locus, which allowed us to estimate the maximum possible value of differentiation (Hedrick, 2005). To visualise a possible pattern of genetic structure, a PCoA was performed in GenAlEx with the obtained F _ST values.

To determine whether the genetic structure was driven by population distribution and geographic distance, we quantified IBD as the correlation between linearised F _ST (F _ST /1- F _ST ) and log-transformed geographic distance in kilometres using a Mantel test (Mantel, 1967) with 10,000 permutations in GenAlEx. Geographic distance was calculated as the shortest path between sampling locations across the sea. These analyses were also performed independently for the populations in the three Mediterranean sub-basins (western and eastern Mediterranean and the Adriatic Sea). Furthermore, to identify potential barriers to gene flow for C. caespitosa, such as oceanographic fronts or currents, pairwise F _ST values and population geographic coordinates were included in the analyses performed in BARRIER v2.2 (Manni et al., 2004). To calculate barrier robustness, we generated 100 resampled bootstrap matrices in R (R Development Core Team, 2017), using an R script provided by Eric Petit (UMR ECOBIO CNRS, Paimpont, France).

Population genetic structure was assessed by a Bayesian clustering approach performed in STRUCTURE 2.3.4 (Pritchard et al., 2000). This program calculates population allele frequencies and, on the basis of HWE estimates, the probability of an individual belonging to one of the obtained genetic clusters. As we assumed that individuals can have ancestors from different locations, an admixture model was implemented with correlated allele frequencies and location as a prior. The analysis was run with a burn-in of 10,000 iterations and 100,000 Markov chain Monte Carlo, with a putative K of up to 25 (six clusters more than the number of sampling sites considered in the analysis) and 20 replicate runs. To infer the number of genetic clusters (best value of K), we applied three methods. We first used STRUCTURE HARVESTER (Earl and vonHoldt, 2012) and CLUMPAK (Kopelman et al., 2015) to evaluate the optimal K value following the method proposed by Evanno et al. (2005) to calculate ΔK. CLUMPAK was also employed to identify the K for which Pr (K = k) is the highest using ln (Pr (X/K)). Although the Evanno ΔK method may efficiently capture the uppermost level of structure, it may underestimate the number of groups (K) (Puechmaille, 2016). Therefore, we also estimated the number of clusters with the method of Puechmaille (2016) using StructureSelector (Li and Liu, 2018). We used the CLUMPAK web server to find the best alignment of the STRUCTURE results across the 20 replicates for each K.

The R package adegenet v2.1.1 (Jombart, 2008) was used to perform the DAPC. This method divides genetic variability of variance between and within groups and then optimises the variance between groups and minimises it within groups (Jombart, 2012), without considering whether population are in HWE.

To quantify the molecular variance, we ran the AMOVA in ARLEQUIN v3.11 (Excoffier et al., 2005) with 1000 permutations for all the populations and with the groups inferred by the STRUCTURE analyses. Putative first generation migrants were identified with a Bayesian assignment method (Rannala and Mountain, 1997) in GENECLASS2 (Piry et al., 2004).

Finally, we tested for recent reductions in effective population size using the allele frequency data in BOTTLENECK v1.2.02 (Cornuet and Luikart, 1996). Simulations were run under three mutation models: the infinite alleles model (IAM), the stepwise mutation model (SMM) and the two-phase mutation model (TPM) (Di Rienzo et al., 1994), with 10,000 iterations and a descriptor of allele frequency distribution (‘mode-shift’). Two statistical tests, the Sign test (Cornuet and Luikart, 1996) and the Wilcoxon sign-rank test (Luikart et al., 1998), were used to test each model.

Of the 541 individuals genotyped, 482 (89%) had unique multilocus genotypes (N _g ) (Table 2). A total of 59 ramets belonging to 25 genets was detected. Clonality percentages of the populations varied from 0% (CDP, COL, PAL, OTR, TOG, TRE, POR and NEA) to 68.75% in CRE. Genotypic evenness (G _o /N _g ) ranged from 0.21 for LIO to one for CDP, COL, PAL, OTR, TOG, TRE, POR and NEA. The lowest values of genotypic richness (N _g /N) were detected for CRE and LIO. Even though these two populations, which have more clonal individuals than the others, had similar values of genetic richness, the observed genotypic diversity (G _o ) of LIO was half that of CRE, reflecting a difference in genet size (LIO had fewer but bigger genets). Genotypic diversity (G _o /G _e ) varied from 0.08 for LIO, indicating “mostly asexual” reproduction (also for CRE and MEN), to one for CDP, COL, PAL, OTR, TOG, TRE, POR and NEA, indicating “sexual reproduction”. The others populations presented “mostly sexual” reproduction by the same criterion. The population with the lowest D index value was LIO, agreeing with the results obtained for genotypic diversity and genotypic richness. Due to the low number of unique multilocus genotypes detected for CRE (N _g = 5), we excluded this population from further analysis.

Linkage disequilibrium (LD) among loci was detected in some populations (in KRY for all loci, and in LIO, KOR, POC, COL, POR, TOG and OTR for some of them). Nevertheless, LD was not observed in any of the loci for all the populations; therefore, all loci were considered independent.

Standardised allelic richness (N _a ) across loci for each population ranged from 3.167 for LIO to 4.859 for NEA (Table 3). The mean value across all loci and populations was 4.412. Observed heterozygosity (H _o ) ranged from 0.477 in COL to 0.587 in LIO. Only LIO presented a H _o value greater than expected (H _e ) (0.587 and 0.496, respectively) and, therefore, a negative value for F _IS (-0.195), indicating a heterozygote excess. The rest of the populations presented lower values of H _o than H _e and positive F _IS values ranging from 0.042 in GAV to 0.16 in BON, indicating a heterozygosity deficit (or a homozygosity excess) and a high degree of inbreeding. All populations, except GAV and MEN, deviated from Hardy-Weinberg equilibrium (HWE). Effective population size estimates were relatively low, especially for MEN (N _e = 4.9), LIO (N _e = 0.4) and KRY (N _e = 4.9). In contrast, BON, TOG, KOR and NEA showed a theoretical infinite value.

The global value of F _ST revealed a low but significant level of genetic differentiation (F _ST global = 0.043, p = 0). Pairwise F _ST values ranged from zero for CDP vs. GAV, CDP vs. COL and CDP vs. PUN (that is, among the western populations) to 0.134 for COL vs. LIO (i.e., between one western and one eastern population) (Figure 2; Table 4). The standardised F _ST value for COL vs. LIO was 0.279. Overall, three ranges of F _ST values were detected: from 0 to 0.02 (usually below 0.01) for within region comparisons, 0.02 to 0.06 for comparisons between populations from different regions and 0.06 to 0.13 for comparisons between any of the populations and the Cypriot ones (LIO and KRY), including between them. Notably, we observed relatively low F _ST values between the Greek population (NEA) and the Ionian (POC, F _ST = 0.039) or the Sardinian (PAL, F _ST = 0.037) ones, and high values between KOR and the other Adriatic populations.

The first two axes of the PCoA explained 55.36% of the variation in F _ST (Figure 3). In this analysis, the populations grouped according to the three main divisions of the Mediterranean basin (western, central and eastern). The western populations from the Spanish peninsular coasts (CDP, PUN and COL) grouped with those from the Balearic Formentera Islands (ESP and GAV). Populations located on the eastern side of the western sub-basin (PAL, BON and MEN) were also differentiated mainly by the first axis. The Adriatic populations, representing a distinctive sub-basin in the middle of the Mediterranean Sea, formed the second group, though with a clear separation between the northern populations of KOR and POR and the rest of the Adriatic populations. With respect to the eastern populations, LIO and KRY in Cyprus were separated from NEA in Greece by the second axis, though all three were clearly separated from one another.

To test whether an IBD pattern was present, we performed a Mantel test with the entire data set and with each of the three groups detected by the STRUCTURE analysis (see below). A pattern of IBD was only observed when all of the populations were analysed together (R ² = 0.28; p = 0.001). No significant association between genetic differentiation (F _ST ) and geographic distance was observed when the populations were analysed separately according to sub-basin (Supplementary Figure S1).

In the STRUCTURE analyses, three genetically differentiated clusters were detected with Evanno’s K method (Figure 4A). The results based on Pritchard’s method, which measures the probability (ln (Pr (X/K)), showed K = 9 as the best value. Under Puechmaille’s method, the populations of C. caespitosa were divided into seven clusters (Figure 4B). Results for K = 3 showed a clear division between the three Mediterranean sub-basins (western, central and eastern). The Tyrrhenian (PAL) and Ligurian (BON) populations showed a certain degree of admixture with Espardell Island (ESP) in Formentera. Other populations also showed admixture including KOR in the Adriatic Sea and NEA in the eastern Mediterranean. The results obtained with Puechmaille’s method (K = 7), though generally similar to those with Evanno’s method, showed a substructure within the basins. For instance, in the western Mediterranean, MEN, PAL and BON showed a high percentage of allocation to a particular cluster (besides the most frequent one in the western sub-basin). In addition, the singularity of KOR within the Adriatic Sea was evident, and LIO, NEA and KRY could each be distinguished as a separate group. When we considered K = 9, the results remain almost the same as those with a K = 7. Analysis of only the western populations revealed three sub-clusters: one differentiating MEN, PAL and BON; another differentiating COL and PUN and a third consisting of GAV and ESP, which are characterised by different degrees of admixture (Figure 4C). Cabo de Palos (CDP) showed similarities to both the COL and PUN and GAV and ESP sub-clusters. Supporting the general results of these analyses, the DAPC showed that the Adriatic and the western Mediterranean populations formed two separate groups, and that the Cypriot populations (LIO and KRY) and NEA were separate from each other and from the other groups (Figure 5A). Moreover, the DAPC of only the western individuals showed some distinction between the island localities MEN and PAL and the Ligurian one (BON) (Figure 5B).

The AMOVA of the three main groups detected by STRUCTURE (western, central and eastern Mediterranean) revealed that 94.67% of the genetic variation originated within populations, with only 2.84% of the variance being attributed to differences among groups. The lowest percentage of variation observed (2.49%) was among populations within groups. Similar results were obtained when seven groups were considered, with 94.75% of the genetic variation originating within populations (Table 5).

In the migration analysis, high rates of first generation migrants were detected (Table 6). The populations with the highest numbers of migrants were CDP (10 migrants), PUN (6), POC (9) and BOK (6). The most important sink populations were ESP (8) and KOR (6). Migrant flow was most important at an intraregional level, with new individuals coming from populations within the same basin. Given that migrant detection analysis using data comprised of low F _ST values may be imprecise, we used the STRUCTURE results for the assignment analyses (Supplementary Table S1). These results showed that the individuals from the Formentera Islands were, in part, assigned to one of the same genetic cluster as those from Menorca Island and the Italian populations (PAL and BON). The population from the Ionian Sea (POC) showed genetic admixture with all three Mediterranean sub-basins. In the Adriatic Sea, the northern populations TRE and POR were assigned to one of the same genetic clusters as KOR, the most differentiated location. Between the Cypriot populations, LIO showed putative migrants from KRY.

The results of the BARRIER analysis revealed six supported barriers to gene flow (bootstrap value >95) for C. caespitosa in the Mediterranean basin. In the western Mediterranean, two main barriers appear to be present: one separating the north-eastern Spanish populations (COL and PUN) from those in the Formentera Islands (GAV and ESP), and another separating ESP (Formentera) from the Menorca population (MEN) (corresponding to the Balearic Front and Mallorca Channel respectively, see Figure 1). A barrier also appears to separate the northernmost western populations (BON and PUN) from one another. A barrier that runs between Sicily and North Africa (Sicilian Channel) separates all of the western populations from the central and eastern ones (Figure 1). In the Adriatic Sea, a north-south barrier appears to divide the eastern coastal populations from the western ones, except in the north (for POR and TRE). Between the eastern and the central populations, a supported barrier was only observed between OTR and NEA (Aegean Front in Figure 1). Finally, in relation to population bottlenecks, only LIO displayed a mode shift, indicating it has experienced a recent bottleneck (Table 7).

Cladocora caespitosa displays both sexual and asexual reproductive modes, and local populations display differences in the prevalence of these modes. Eight populations across the Mediterranean Sea (three western, four Adriatic and one eastern) appear to undergo sexual reproduction exclusively, as no evidence of clones was found. Sexual reproduction is also dominant in most of the other populations as only a low level of asexual reproduction was observed in them. These results are consistent with those of a previous study showing that sexual reproduction is predominant in C. caespitosa (Kružić et al., 2008). By contrast, asexual reproduction is dominant in CRE and LIO (>60% clones), and to a lesser extent, in MEN (38.1% clones). These values are comparable to the level of clonality shown by coral species that reproduce mostly by fragmentation, such as Pavona clavus or Acropora valida (Aranceta-Garza et al., 2012). The high rate of asexual reproduction in the LIO population was recently related to local natural or anthropogenic disturbances (López-Márquez et al., 2021). This phenomenon has also been observed in other scleractinian corals, such as in populations of Pocillopora verrucosa after tropical storms (Aranceta-Garza et al., 2012). Currently, we do not know the specific circumstances that have led to the high level of asexual reproduction in the CRE and MEN populations, though we hypothesise that they also experienced local disturbances that resulted in high mortality rates. In Crete, for instance, tourism has greatly impacted the coastlines (Tsilimigkas et al., 2020), which may contribute to the high level of asexual reproduction observed in C. caespitosa in this area. In Menorca, colonies are located in a sheltered place with a low swell (J. Templado, personal observation), thus bioerosion may be weakening these colonies, leading to polyp detachment. The small estimated effective size of MEN also suggests this locality may have been affected by environmental disturbances that could have triggered asexual reproduction in these colonies, as has been observed in other coral species (Lirman, 2000). Therefore, it is accepted the hypothesis that the clonal structure varies across its geographic range due to different disturbances affecting the sexual/asexual reproductive ratio.

Despite having a low level of genotypic diversity (G _o /G _e ), the genetic diversity (H _o ) of MEN was similar to nearly all of the other analysed populations of C. caespitosa, and analogous to those reported for other tropical corals (Baums, 2008). This finding suggests that the level of sexual reproduction in this population is enough to maintain its diversity (Bengtsson, 2003). In concordance with other studies (Baums, 2008; Polato et al., 2010; Casado-Amezúa et al., 2012), we observed positive values for the inbreeding index, a heterozygosity deficit in most populations and deviation from HWE in all populations, except GAV and MEN. Deviation from HWE in marine invertebrate populations may be associated with their sexual reproductive mode. For instance, inbreeding is widespread in most marine benthic invertebrates with limited dispersal capabilities (Riesgo et al., 2016). In C. caespitosa, fertilisation takes place in the surrounding water of the broadcasting colonies without sexual selection; therefore, inbreeding is prevalent in this sessile species (Addison and Hart, 2005).

Given the life history traits of C. caespitosa, such as its supposedly low dispersal capability and high self-recruitment rates (Kružić et al., 2008; López-Márquez et al., 2019; this study in which a rate greater than 80% was estimated), we hypothesised that the species would have a strong population structure with a high level of differentiation. In contrast, and rejecting this hypothesis, we found that the populations studied had, in general, relatively low F _ST values among some nearby populations, indicating a low level of differentiation. However, we did observe a regional structure in the western Mediterranean according to the F _ST values: populations from the peninsular coasts of Spain (CDP, PUN and COL) and the Balearic island of Formentera (GAV and ESP) formed a single group, whereas the Menorca island population (MEN) was genetically more similar to the populations in the Tyrrhenian (PAL) and Ligurian (BON) sub-basins (see Figure 2). This substructure is likely mediated by dispersal along the coastline and some resistance to larval movement in the open sea since lower F _ST values were observed between the coastal populations than between these and the island ones. Similar results were found in the Adriatic and adjacent Ionian populations (López-Márquez et al., 2019) and in the eastern Mediterranean populations (López-Márquez et al., 2021), where the gene flow resulted to be mainly driven by the shoreline configuration and sea surface current which likely enhanced the dispersal of planulae among populations. Genetic differentiation of the Croatian locality KOR from its neighbours in the Adriatic Sea was evidenced by the relatively high F _ST values observed among these populations. López-Márquez et al. (2019) showed the same results with a different set of markers and inferred that this was due to the influence of a north-south barrier. Although we also observed this latitudinal separation with our data set, the inclusion of POR, the northernmost Adriatic population, suggests this differentiation is not as strong as previously observed. The F _ST values shown between POR and the other Adriatic populations were not as high as those between KOR and these populations. The north-south differentiation observed in the Adriatic Sea is likely mediated by dispersion along the coastline and some dispersal events through subgyres that promote connectivity between the southern populations. Aside from its northern location, KOR stands out as a differentiated cluster. With respect to the eastern Mediterranean populations (KRY and LIO, in the Levantine Sea, and NEA, in the Aegean Sea), pairwise F _ST comparisons indicate that all are differentiated from one another; however, in the STRUCTURE analysis (for K = 3), they form a single cluster.

The population structure resulting from the various analyses (i.e., STRUCTURE, DAPC and PCoA) coincides with the three main divisions of the Mediterranean Sea (western, central and eastern Mediterranean). According to the Mantel test results, this is likely due to geographic distance (IBD) along the eastern–western distribution axis of the Mediterranean populations (see Supplementary Figure S1).

In addition to geographic distance, barriers to gene flow have likely influenced the current population structure of C. caespitosa. The barriers we identified largely coincide with those found in other population genetic studies (Pascual et al., 2017; Constantini et al., 2018). For this reason, the results suggesting the presence of long-distance first generation migrants in NEA and KRY were unexpected. One possible explanation could be the evolutionary convergence of alleles from different locations, which would lead to low pairwise differentiation and, consequently, indicate the presence of migrants. The fact that the standardised F _ST values between the eastern Mediterranean and Adriatic and Ionian populations are higher than for other pairwise comparisons supports this premise. Otherwise, anthropogenic dispersal should be taken into account (Radziejewska et al., 2006): human-mediated transport, for instance, through lithic anchors or ballast stones in Phoenician times could have also led to this surprising result.

Regarding the intra-regional substructure detected by the STRUCTURE and assignment analyses, gene flow is not mediated by distance, but rather by ocean currents and dispersal along the coastlines. Consistent with this, the Balearic Front (BF) and the Mallorca Channel (MC, Ruiz et al., 2009) appear to act as barriers dividing CDP and the two Formentera populations (GAV and ESP) from PUN and COL. However, evidence of migrants and a low to moderate level of genetic differentiation indicate that these barriers are permeable. The NBF is seasonal, and it depends on the strength of the Northern Current (Font et al., 1995). As such, during autumn and winter, its intensity decreases and the barrier is less strong. This event almost matches the spawning period of C. caespitosa in this region (Kersting et al., 2013b), which would allow for some level of gene flow between populations. Another ocean front of importance is the Northern Tyrrhenian Gyre (Poulain et al., 2012), which facilitates connectivity between the populations from the Tyrrhenian and Ligurian seas (PAL and BON, respectively). On the other hand, connectivity among the eastern populations seem to be influenced mainly by a pattern of isolation by environment, consistent with the results shown by López-Márquez et al. (2021).

Lastly, although we only obtained five unique multilocus genotypes for CRE, leading us to exclude it from further analysis, our results on this population’s genetic structure (Supplementary Figure S2) were potentially interesting. Greater sampling effort of this area, however, is needed to confirm this hypothesis.

As discussed here, various factors affect the genetic differentiation of the different populations of C. caespitosa in the Mediterranean. The population structure of this coral mainly appears to be correlated with the division of the main Mediterranean sub-basins, and the barriers between them, together with favoured dispersal along coastlines. Despite the low dispersal ability, high rate of self-recruitment and estimated low effective population sizes of C. caespitosa, the occasional dispersive events that seem to occur between regions are sufficient to maintain the species’ considerable genetic diversity. Under suitable conditions and habitat continuity without barriers to gene flow, connectivity between distant colonies of C. caespitosa could be high. In fact, the two localities on the peninsular coast of Spain (CDP and PUN), which are separated by 400 km, did not present any genetic differentiation (pairwise F _ST value = 0).

Despite its importance in conservation policies, knowledge of the genetic diversity and connectivity patterns of species is often not considered by policymakers. This type of information is especially important for the conservation of endangered structural species. Our results on the genetic structure and population connectivity of C. caespitosa, whose populations are in alarming decline, highlight key aspects that should be included in the conservation decision making process for this species. These aspects include the relatively low dispersal ability and the high rate of self-recruitment of the species, the importance of sporadic dispersal events to maintain diversity and the level of clonality of some populations.

The regression of populations of C. caespitosa is still in progress (Rodolfo-Metalpa et al., 2005), and the high frequency of mortality events in this coral over the last decades possibly exceeds its recovery potential because of its low recruitment rates (Kersting et al., 2014a). Therefore, understanding the general pattern of connectivity among populations of this coral, particularly in light of its discontinuous distribution (as isolated colonies or, more rarely, as beds or banks) and low dispersal ability, becomes increasingly important. Known beds and banks of the coral should be considered as source populations from which larvae can be occasionally exported to other areas. Conversely, populations comprised of a few scattered colonies for which successful fertilisation may prove highly difficult can serve as potential recipient populations. For this reason, strict protection of the marine areas where beds and banks of C. caespitosa are known to be present should be a top priority. According to the available information from the literature and expert observations describing the abundance and morphology of colonies (compiled by Chafaoui et al., 2017), only 31 localities across the entire Mediterranean basin have reported bioherms of the species. Although in this study we cannot appraise whether the current level of protection in established marine protected areas is sufficient for the conservation of the species, we strongly recommend an increase in the coverage of this emblematic species within these areas and measures to preserve it from major threats, particularly those that are mediated by human activities. For instance, although extreme climatic events, such heat waves, appear to be a major threat to this coral (Kersting et al., 2013a; Jiménez et al., 2016), other factors posing a threat include the spread of mucilage and invasive algae (Kersting et al., 2014b), eutrophication due to sewage discharge or fish farming (Kružić and Požar-Domac, 2007), trawling, anchorage and high sedimentation rates caused by dredging or extreme storms (Casado de Amezua et al., 2015; Chefaoui et al., 2017; López-Marquez et al., 2021). In addition to the conservation of the known coral beds, prospecting studies of other Mediterranean localities should be conducted to identify as yet unknown populations of the coral.

Besides local impacts, global changes (sea warming and acidification) also threaten C. caespitosa. Global threats such heat waves or extreme storms interact with local ones such as pollution or coastal degradation, leading to cumulative impacts that are particularly high in coastal ecosystems (Micheli et al., 2013). Although local mitigation strategies cannot directly protect populations from the impact of global threats, they are often the only feasible way to reduce the synergy between the different types of threats, and to help preserve marine coastal ecosystems (Templado, 2014). In this context, protection against local impacts is of the utmost importance.

Another aspect to consider in the conservation of C. caespitosa is the role of asexual reproduction, and its effect on the genetic diversity of populations. For sessile organisms that undergo external fertilisation, asexual reproduction may be the local means of proliferation for populations with a small population size and/or low gamete density (“Allee effect”; Courchamp et al., 2008). Clonal growth, therefore, may allow the species to persist through periods of low sexual recruitment (Lasker and Coffroth, 1999); however, it may also lead to reduced genotypic diversity and, as a result, higher susceptibility to environmental changes (Reusch et al., 2005). From a conservation point of view, knowledge of the genetic and genotypic diversity patterns of structural species with a high potential for asexual reproduction is critical: populations with a high level of genotypic diversity may be better able to withstand environmental changes or extreme climatic events (Reusch et al., 2005) and, conversely, those with a low level of genetic diversity may be more vulnerable to pathogens and parasites (King and Lively, 2012). For C. caespitosa, we found that the ratio of clonal to sexual recruitment is highly variable between localities, and is likely a consequence of a complex interplay of various impacts. Future studies should focus on determining the factors that affect this ratio and how different impacts interact synergistically on it, as well as on the historical causes for the decline of bioherms of the species. Increased clonality may negatively impact the adaptive potential of the species under the current disturbance regime; therefore, it should be closely monitored.

In summary, enforcement of conservation measures that aim to protect the genetic diversity of endangered structural species such C. caespitosa is essential. The results provided here, combined with those from a recent study on another endangered Mediterranean coral, Astroides calycularis (Ledoux et al., 2021), may help build a highly relevant framework in which to study the evolution of Mediterranean marine diversity as it faces the increasing warming of the basin waters.