Asexual reproduction in bad times? The case of Cladocora caespitosa in the eastern Mediterranean Sea.

We analysed the patterns of genetic variability of eastern Mediterranean populations of the scleractinian coral Cladocora caespitosa, from the Aegean and Levantine seas, using 19 polymorphic microsatellite loci, 11 of which were newly characterized. The observed genetic pattern reflects a scenario of isolation by environment: F ST comparisons showed a higher degree of genetic differentiation between the two Cypriot populations that are separated by only 11 km than between these two Levantine populations and the Aegean population in Greece, which are separated by 1300 km. We hypothesize that local-scale oceanographic factors influenced the dispersal of planulae between the geographically close populations, playing a crucial role in the genetic structure of this coastal coral. Yet, despite being characterized as a species with limited dispersal and high self-recruitment, large-scale migration does eventually occur as first-generation migrants were identified between the most distant populations. In line with previous findings of reproductive plasticity in C. caespitosa, we also found localized differences in reproduction mode (sexual vs. asexual) within a geographically limited context. Several individuals were identified as clones, indicating the predominance of asexual reproduction in one of the Cypriot populations. We interpret this predominance either as a direct response to or as an indirect consequence of perturbations suffered by this C. caespitosa population. These perturbations are caused by unfavourable environmental conditions that threatened local survival, in particular water temperature changes and windstorm swells. Asexual reproduction may be a mechanism used by C. caespitosa to counteract mortality events and recolonize devastated areas, and likely accounts for the occasional high levels of clonality and low levels of genetic diversity. Local adaptations such as these should therefore be considered in conservation and management strategies to maintain and preserve the gene pool of this endangered species. SUPPLEMENTARY INFORMATION: The online version contains supplementary material available at (10.1007/s00338-020-02040-3).

Introduction

The colonial scleractinian coral Cladocora caespitosa is distributed across the entire Mediterranean basin, with scant occurrences in the adjacent eastern Atlantic Ocean (Zibrowius 1980, 1983). Cladocora caespitosa shows a patchy distribution, commonly as isolated colonies, although at some localities, it develops monospecific reef-like structures (e.g. Peirano et al. 1998). This versatile coral hosts, in its tissues, endosymbiotic dinoflagellates of the genus Symbiodinium (Casado-Amezúa et al. 2011) and occurs in a wide range of habitats, from shallow photophilic communities to circalittoral assemblages (Bellan-Santini et al. 2002). Although not resistant to strong hydrodynamics and waves, intense currents and a complex and irregular coastline and bottom topography favour the settlement of this coral (Chefaoui et al. 2017). Cladocora caespitosa is the only recent zooxanthellate taxon reminiscent of the true coral reef ecosystems that inhabited the Mediterranean basin up to the late Miocene (Vertino et al. 2014; Kersting and Linares 2019). It has been present in the Mediterranean since the early Miocene (Vertino et al. 2014), and putative morphological species of C. caespitosa have been found in Pliocene and early Pleistocene deposits of the Mediterranean basin (D’Alessandro and Bromley 1995; Peirano et al. 1998; Dornbos and Wilson 1999; Spadini 2015; Borghi 2019). Most occurrences, however, date to the late Pleistocene as C. caespitosa was especially common in marine terraces formed during the last interglacial (MIS5e, ca 125 kyr BP) (Cuerda 1975; Copat Marconi et al. 1982; Peirano et al. 1998; Amorosi et al. 2014).

Sexual reproduction in C. caespitosa is generally synchronous and seasonal; however, notable differences in some reproductive traits have been found between populations in the Adriatic Sea (central Mediterranean) and those in the western and eastern Mediterranean. For instance, in the Adriatic, hermaphroditic colonies reproduce at the beginning of the summer, when temperatures begin to rise, typically coinciding with a full moon (Kružić et al. 2008). By contrast, in the western and eastern Mediterranean, gonochoric colonies release gametes at the end of summer, when temperatures begin to fall (Kersting et al. 2013b; Hadjioannou 2019).

During sexual reproduction, polyps release mucus-coated eggs and sperm bundles and fertilization takes place in the surrounding water (Kružić et al. 2008). Larvae do not appear to typically disperse great distances, given the clumped distribution of the species (Kersting and Linares 2012). Asexual reproduction could increase the species’ tendency for local retention, as occurs in other corals due to polyp budding or fragmentation (Kružić et al. 2008). The negative buoyancy of the eggs may also favour retention mechanisms (Kersting et al. 2014). Overall, the reproductive traits of this coral appear to limit dispersal and promote self-recruitment, such that long-distance dispersion is sporadic and recruitment rates are low (Casado-Amezúa et al. 2014; Kersting et al. 2014).

The low recruitment rate, recovery potential and dispersal capability of this long-lived species, combined with the long-term impacts of climate change (Kersting et al. 2013a), have led to its decline in distribution and overall abundance. Human impacts in coastal areas and the effects of heat waves and severe storms have also contributed to the species’ decline. Moreover, given that C. caespitosa may already be living near its upper thermal limit (Rodolfo-Metalpa et al. 2006), changes such as those anticipated by the present climate scenario, e.g., more heat wave episodes in the Mediterranean (Russo et al. 2014; Garrabou et al. 2019), may increase the species’ susceptibility to temperature increases, even small ones. For these reasons, C. caespitosa has been listed as an endangered species on the IUCN Red List (Casado-Amezúa et al. 2015).

Peripheral populations usually show less (reduced) genetic variability (Lesica and Allendor 1995); therefore, they are of special concern and should be surveyed in order to determine their conservation needs and to maintain the genetic pool of certain species. Although the genetic parameters of western Mediterranean and Adriatic populations of C. caespitosa have been analysed (Casado-Amezúa et al. 2014; López-Márquez et al. 2019, respectively), similar studies have not been carried out on populations in the eastern basin. In addition to their marginal condition, eastern populations of C. caespitosa have been recently impacted by warming events and windstorms (Hadjioannou et al. 2016; Jiménez et al. 2016; Hadjioannou 2019). Cyprus is a key area of historical significance for C. caespitosa due to the species’ prolonged presence in its coastal waters during the Plio-Pleistocene (Dornbos and Wilson 1999) and likely continuous presence since the last interglacial (Galil et al. 2016). Given the phylogeographic and biogeographic implications of its historical presence in Cyprus, we chose to perform genetic assessments on current Cypriot populations of C. caespitosa. We analyse the genetic variability of two populations from Cyprus and one from Greece using both newly developed and previously described hypervariable microsatellite markers (Casado-Amezúa et al. 2011). The previous markers were used to identify population structure patterns and to assess connectivity among populations in the western Mediterranean (Casado-Amezúa et al. 2014) and the Adriatic (López-Márquez et al. 2019). Although particular reproductive features may have produced a pronounced population genetic structure, these markers showed a low level of genetic differentiation among locations, indicating a low power of resolution. We, therefore, designed new microsatellites to improve the analyses and to provide greater information on the genetic parameters characterizing endangered C. caespitosa populations in the eastern Mediterranean.

Materials and methods

Study site, sample collection, DNA extraction, PCR amplification and microsatellite isolation

The study area was comprised of three localities in the eastern Mediterranean basin: one in Greece along the northern Aegean coast (Nea Peramos, Greece = NEA) and two along the southeastern coast of Cyprus (Liopetri = LIO and Kryo Nero = KRY) (Fig. 1; Table 1). All necessary permits were obtained for field studies from the relevant authorities. The two coral communities sampled in Cyprus, which are separated by ~ 11 km and are located at the same depth (~ 3–4 m deep), have been monitored and studied since 2012 (Jimenez et al. 2016; Hadjioannou et al. 2016, 2019; Hadjioannou 2019). Kryo Nero (KRY) is situated at a small cove protected by cliffs in the southeast area of Cyprus. It hosts approximately 100 C. caespitosa colonies, and the average colony size is 42.5 ± 29.8 cm (Hadjioannou 2019). The population in Liopetri (LIO), which is west of KRY, is comprised of more than 200 small colonies (with an average size of 30.25 ± 13.6 cm) that are settled in marine terraces close to the coastline, which are exposed to incoming swells. Surface water circulation along these two coastal localities is similar. During most of the year, a near-shore current flows eastward from Cape Pyla (west of LIO) to KRY and is deflected offshore (in a SE direction) where it meets the southward coastal current coming from the easternmost point of the island. Both localities are exposed to large north and northwest swells, primarily during the winter months; however, turbulence conditions are more pronounced in KRY than in LIO (Menna et al. 2012, Poulain et al. 2012, Politikos et al. 2017 and personal communication from D. Hayes). This is due to differences in the geomorphology of the coastline: KRY is comprised of a series of steep cliffs, small coves and sea caves, whereas in LIO, the shoreline is more uniform or rectilinear. The environment inhabited by the two populations also has contrasting features. KRY is a naturally low-nutrient area, whereas LIO is a high-nutrient one due to the impact of anthropogenic activities, such as inflow from a nearby fish farm hatchery and agriculture-related impacts (Jimenez et al. 2016; Hadjioannou 2019).

Fig. 1

Sampling locations of Cladocora caespitosa in the Aegean and Levantine seas. The general surface circulation (grey arrows) and predominant local currents (dark blue arrows) of the studied areas are shown. Main eddies are indicated by the dotted lines.

Adapted from Menna et al. 2012, Poulain et al. 2012, Politikos et al. 2017 and personal communication by D. Hayes

Table 1

Location and number of Cladocora caespitosa samples (N) collected from the eastern Mediterranean

Location	Label	GPS coordinates	N
Nea Peramos, Greece	NEA	40°49′31.9″N 24°20′01.9″E	31
Liopetri, Cyprus	LIO	34°57′30.2″N 33°54′05.7″E	31
Kryo Nero, Cyprus	KRY	34°58′57.0″N 34°01′00.8″E	30

The Aegean population (NEA) is ~ 1300 km from the Cypriot ones and has been monitored since 2015. It consists of > 50 small colonies attached to rocks (~ 3 m deep) that are adjacent to the shore. In addition, there are small free-living coral nodules or coralliths (sensu Glynn 1974) of C. caespitosa, mixed with rhodoliths, scattered among rubble in a large Posidonia oceanica meadow that extends to deeper waters (3–10 m). The study site is located on the southern flank of Eleftere Bay. Prevailing north and north-eastern winds favour a southward coastal circulation inside the bay, and turbulent conditions in the shallower areas of the Posidonia meadow.

Individual polyps from a minimum of 30 randomly selected C. caespitosa colonies were collected by SCUBA diving, at a depth of 3 m, during October and November of 2017. Samples were collected from colonies that were at least 1 m apart in order to avoid sampling clones in each of the three localities. Samples were stored at 4–10 °C in vials of absolute ethanol for subsequent genetic analyses.

The colonies sampled in the present study have been the subject of long-term monitoring by some of the authors; thus, we were fairly confident in our initial species identification. However, to avoid any doubt, we used molecular means to confirm that all collected specimens belong to the same species, namely C. caespitosa. We selected 8–9 individuals from each population to sequence two genes (one mitochondrial and one nuclear). Genomic DNA was extracted and purified using the QIAGEN BioSprint 15 DNABlood Kit (Qiagen), according to the manufacturer’s protocol. The mitochondrial cytochrome oxidase subunit I (COI, 658 base pairs) and the nuclear ITSs (3′ end of the 18S RNA + ITS-1 + 5.8S + ITS-2 + 5′ end of the 28S RNA, 654 base pairs after alignment) were amplified with the same primers and conditions as described by Merino-Serrais et al. (2012). Due to the presence of multiple copies of the ITSs fragment, the PCR products were cloned. Amplicons were gel-extracted, ligated into the pGEM(R)-T Easy Vector System and transformed in JM109 competent cells (Promega). Recombinant colonies were identified by white-blue selection on ampicillin-coated agar plates. Subsequently, up to eight colonies per specimen were further validated through amplification with ITSs primers. Finally, up to six of the positive amplicons were sequenced using M13 universal primers. Haplotype networks were visualized using HapView software (Salzburger et al. 2011). Sequences were deposited in GenBank (accession numbers: MW032500 to MW032523 and MW033407 to MW033523).

For the microsatellite isolation, genomic DNA was extracted from a naturally bleached colony from Espardell Island (Formentera, Spain) in order to avoid the risk of zooxanthella DNA contamination. Library preparation of the microsatellite-enriched fragments was performed at the Sequencing and Genotyping Facility at Cornell Life Sciences Core Laboratory Center, as previously detailed by López-Márquez et al. (2016). A total of 250 potential contigs were identified as containing microsatellites. Using QDD software, sequences with sufficient flanking areas were selected for primer design and then filtered by selecting for hexa-, penta-, tetra-, tri- or dinucleotide repeat motifs. Primers were designed using PRIMER 3 v.0.4.0 (Untergrasser et al. 2012) with the same criteria used by López-Márquez et al. (2016). A total of 23 primer pairs were selected and combined in five tetraplex and one triplex reactions in an initial screen of potential microsatellites. Genomic DNA was then extracted and purified from the rest of the 92 sampled individuals. After quantification, each DNA sample was diluted to a final concentration of 0.3 ng/µl. Forward primers were fluorescently end-labelled with 6-FAM, NED, VIC or PET (Table 2). Amplifications were carried out in a total volume of 7 µl containing 1X Qiagen Multiplex PCR Master Mix, 0.2–0.3 µM of each primer, 0.30 ng of DNA template and water. The PCR cycling profile was 95 °C for 15 min, 35 cycles at 94 °C for 30 s, 56 °C for 90 s and 72 °C for 30 s, and a final extension at 72 °C for 10 min. Fluorescently labelled PCR products were run on an ABI PRISM 3730 DNA Sequencer (Applied Biosystems), and alleles sizes were determined according to the ABI GS-500 LIZ standard. A total of 19 polymorphic microsatellite loci were utilized to characterize genetic connectivity: 8 specifically developed for C. caespitosa by Casado-Amezúa et al. (2011) and 11 novel ones developed by Illumina MiSeq next-generation sequencing for this study.

Table 2

Cladocora caespitosa microsatellite characterization

Locus name	Primer sequences	Repeat motif
*Cc-L2	F: 6-FAM-CGTGTAAATGCCACCAAACA R: TTACTTTGGGTGCCCCATTA	(TG)3CGAG(TG)2GG(TG)5AG(TG)7
*Cc-L4	F: PET-CAACAGGAAGCTGAAGCTGA R: CTTGCGCTTGTTC	(AAG)4AAT(AAG)16GAG(AAG)9(AAC)10CAC(AAC)3
*Cc-L16	F: NED-TTGCCATTTTAAACAGTTAC R: TCTGAGACTAGAGTGAGTGC	(GGGGGT)5
*Cc-L19	F: 6-FAM-TTTGACGATTATTGTATGCT R: GTATTGCTGTGTTTTTGC	(AAC)7ACCAGCTAC(AAC)2
*Cc-L21	F: NED-AACGTGAAGAAATACAAGTG R: AGAAAAGGACTTTTATAACGA	(AAC)16AGCAACAGC(AAC)4
*Cc-L27	F: NED-CTTAAACTAATTGCAAAAGG R: AAAGTTCAGATAGCAAAGGT	(AAC)8
*Cc-L29	F: VIC-TGTGGGAGCTGGGCCGAGTA R: AAATTCGACGGCAGCACATG	(AAC)12
*Cc-L37	F: VIC-GATTGATGTCAGTCACTCTG R: GCATGAGAAACGAAGATAG	(AAC)3(AAT)4(AAC)4(AAT)4(AAC)4(AAT)3(AAC)7
CcV11	F: NED-TTCAAAATGATGCCAGGACAA R: TTTTCGTCCGGAAAACTACGA	(TTC) 11
CcV18	F: VIC-CGGTTCACATCGCTGACAT R: AACAGGAGACACCCTCAACC	(CAG)6
CcV22	F: NED-TTATGTGGATCGTTAGCCTGT R: GAGTAGGCCACGAAATCCAA	(TCAA)11
CcV23	F: 6-FAM-CCTTATTGCTAAGCGTATGACAC R: TTGACTTTTATTGCTCCTTGCC	(ATTGT)12
CcV26	F: 6-FAM-TTGTTGTGCTTCAGCCTCAG R: TTTGCCAATGTGACCAGGTA	(GATG)7
CcV34	F: PET-GGAAGCCCTGGAACATATTG R: CTCTAGTTACCAGGCATGAACAA	(TCAA)11
CcV35	F: VIC-TGGCCTCGCCTATAAAACTC R: CAAGGGCTAGTTCTGATGCTC	(ATG)10
CcV38	F: 6-FAM-TGATGCAGGTTTATCCGTGT R: CGTGCCGACCACGTTAATA	(AATC)6
CcV42	F: 6-FAM-AAGTTCATTAGTTTTGTTCATTACTGC R: CCACGGAAACCCTAAAAGAA	(CT)3AAATGAAGT(TTGA)8
CcV46	F: PET-GAGACATGTAGAGAATCTGGTGCAG R: AAATAGCAGAATTCATGGATCACG	(TTA)9
CcV47	F: PET-CACGTGTCTTTCTGGCTCTG R: TTCTCATAGCTGCCATGCTG	(ATC)12

Forward primers were 5′ end-tailed with 5′-TGACGACCCCATGCTACG-3′, and reverse primers were pig-tailed with 5′-GTTTCTT-3′. * Primers previously described by Casado-Amezúa et al. (2011)

Electropherograms were analysed with GENEMAPPER software 3.0 (Applied Biosystems). MICRO-CHECKER v.2.2.3 (Van Oosterhout et al. 2004) was used to assess the presence of null alleles and scoring errors.

Genotype analyses

The first step was to analyse the genotypic diversity of each population. Using the function and parameters of “multilocus matches”, as implemented in GenAlEx 6.0 (Peakall and Smouse 2006, 2012), we calculated the number of unique multilocus genotypes present in each population (Ng). The presence of repeated genotypes can also be detected in this analysis. Thus, with these data, we can calculate the number of samples that are clones.

Clonal structure parameters

To study clonal population structure, we analysed the genotypic richness standardized to sample size by calculating the total number of unique multilocus genotypes per site (Ng) divided by the total number of genotyped individuals (polyps) per site (N). Values close to zero indicate that the majority of the sampled individuals have the same genotype (i.e. are clones), whereas values close to one indicate that most have a unique genotype.

Genotypic evenness was calculated by dividing the observed genotypic diversity of a site (Go) by the number of unique multilocus genotypes at that site (Ng) where Go is the observed genotypic diversity calculated as Go = 1/Σg 2 and g is the frequency of the ith genotype in the population. Evenness values close to zero indicate a single genotype (a clone) dominates the population, whereas values close to one indicate that genotype abundance is distributed evenly among the samples (Aranceta-Garza et al. 2012). Genotypic diversity was calculated by dividing Go by expected genotypic diversity (Ge), which is also the total number of individuals genotyped per site (Ge = N). Genotypic diversity and evenness are both indicators of the proportion of asexual and sexual reproduction in a population, following the criteria of Baums et al. (2006) who classified reproduction of populations as “sexual”, “mostly sexual”, “mostly asexual” or “asexual”. We also calculated a simple index measure, D (Pielou 1969): values of D range between zero and one. A D of one indicates that the population is comprised of genets with only one ramet (each individual has a different multilocus genotype), and a D of zero indicates that the population has only one genet (all individuals have the same multilocus genotype and, thus, belong to the same clone).

Genetic variability

Samples identified as clones were excluded from the population genetic analyses. Allelic diversity (Na), observed (Ho) and expected (He) heterozygosities, test of Hardy–Weinberg equilibrium (HWE), FIS inbreeding coefficient and linkage disequilibrium (LD) were calculated with GENEPOP (Raymond and Rousset 1995) and GENALEX. When necessary, significance of P values was corrected using the sequential Bonferroni method (Rice 1989). Genetic similarity among samples from the different locations was estimated with Wright’s fixation index (FST) through Weir and Cockerham’s estimators in GENETIX v.4.03 (Belkhir et al. 2004). Standardized FST values were calculated by dividing the original FST values by recoded ones, which assume that each population has different alleles for each locus, while maintaining observed allelic frequencies. F’ST values (Meirmans 2006) were also calculated in GENALEX. To visualize population genetic clustering, factorial correspondence analysis (FCA) using the function “AFC 2D” was performed in GENETIX and principal coordinates analysis (PCoA) on FST was performed in GENALEX. To characterize population genetic structure, a Bayesian clustering analysis was performed with STRUCTURE 2.2.3 (Pritchard et al. 2000). An admixture model was implemented with correlated allele frequencies and location specified as a prior. For each cluster (K), 20 replicates were performed in two analyses, one for a K up to 5 and another for a K up to 10. The mean log probability of the data (lnP(K)) was calculated based on 100,000 Markov chain Monte Carlo (MCMC) iterations, following a 10,000 iteration burn-in. We evaluated the optimal value of K by considering both the highest mean likelihood value (L(K)) and the ΔK, as calculated using STRUCTURE Harvester (Earl and vonHoldt 2012) following the method proposed by Evanno et al. (2005). Clumpak (Kopelman et al. 2015) was also used to estimate the best-fit K using two methodologies: the one used in STRUCTURE Harvester (ΔK by Evanno et al. 2005) and another that uses ln(Pr(X|K) values (median values of Ln (Pr Data)) to identify the K for which Pr (K = k) is the highest. However, Evanno et al. (2005) demonstrated that the ΔK value is likely to detect only the uppermost levels of genetic structure patterns and should not be used exclusively to identify the “true” number of clusters (Janes et al. 2017). To identify the K that best fits the data, Earl and vonHoldt (2012) recommend plotting an average across all iterations of the estimated logarithm of the probability of the data, that is, the likelihood of K for each value. If there is not a clear maximum likelihood, the point where the plot curvature reaches a plateau can be used (Pritchard and Wen 2003), although it is sometimes difficult to determine this point. In these cases, the authors of STRUCTURE acknowledge that users should then select a biologically sensible value. Clumpak was used to assess fitness of the results across the range of K values. The software StructureSelector (Li and Liu 2018) was used with the Puechmaille method (2016) to estimate the number of population clusters.

ARLEQUIN v3.11 (Excoffier et al. 2005) was then used to hierarchically quantify the molecular variance (AMOVA, n = 1000 permutations) in the groups inferred by the STRUCTURE analyses. Additionally, evidence of selection acting on the different loci was inferred using two methodologies, a Bayesian approach in BAYESCAN (Foll 2012) with the following parameters: burn-in = 50,000, thinning interval = 30, number of outputted iterations = 5000, number of pilot runs = 50 and length of pilot runs = 5000 and, neutrality tests in ARLEQUIN. To identify possible recent migration among populations, a Bayesian assignment method (Rannala and Mountain 1997) was used as implemented in GENECLASS2 (Piry et al. 2004). To calculate individual probabilities of assignment to each population, a MCMC resampling method with a simulation algorithm (Paetkau et al. 2004) was performed using 10,000 simulated individuals and a type I error threshold of 0.05.

The existence of potential barriers to gene flow was assessed by Delaunay triangulation from GPS coordinates using Monmonier’s (1973) maximum difference algorithm with the pairwise FST matrix in BARRIER v2.2 (Manni et al. 2004). Barrier robustness was assessed with 100 resampled bootstrap matrices in R (using an R package provided by Eric Petit, UMR ECOBIO CNRS, Paimpont, France).

We quantified isolation by distance (IBD) among all locations as the correlation between linearized FST (FST/(1 − FST)) and the log of geographic distance using the Mantel permutation test (Mantel 1967) with 10,000 permutations, as implemented in GENALEX. Geographic distance was computed as the shortest ocean path between sampling locations.

Lastly, we assessed potential recent bottlenecks using BOTTLENECK v1.2.02 (Cornuet and Luikart 1996). We looked for effective population size reductions using the allele frequency data. Three mutation models, under two different statistical tests and a descriptor of the allele frequency distribution (“mode-shift”) (Luikart et al. 1998), were tested with 10,000 iterations. The mutation models used were the infinite alleles model (IAM), the stepwise mutation model (SMM) and the two-phase mutation model (TPM) (Di Rienzo et al. 1994). The statistical tests performed were the Sign test (Cournet and Luikart 1996) and the Wilcoxon sign-rank test (Luikart et al. 1998).

Results

After an initial screen, 11 new microsatellite loci, together with 8 previously described ones, were selected for further genetic analyses. All of the analysed loci were polymorphic in the populations except for two that were monomorphic in the LIO population.

Linkage disequilibrium (LD) among loci was not observed for the Aegean population. However, the Cypriot populations showed LD in all loci. In the previous microsatellite studies, LD was not observed in any of the loci included in the present study and all were considered statistically independent (Casado-Amezúa et al. 2011, 2014, López-Márquez et al. 2019). These contrasting results indicate that the LD signal in the Cypriot populations is due to factors other than physical linkage.

Genotype, haplotype diversity and clonal structure

Only two mitochondrial COI haplotypes (differing by only one of 658 base pairs sequenced) were observed in the specimens. The main haplotype was found in all specimens except three belonging to two of the three localities. Haplotype analysis of the nuclear ITSs region, selected because of its greater variability, mostly showed a lack of biogeographic structure. The maximum value of divergence observed among haplotypes was 3.6%, found between two belonging to the same individual. Comparison of the obtained haplotype sequences with GenBank data for C. caespitosa and other closely related species indicated a divergence around 4.3% between these haplotypes and the most similar sequences, namely those identified as C. caespitosa. None of the specimens analysed, including the three with the minor COI haplotype, showed a greater genetic difference among groups than expected for a single species, thus confirming the assignment of the samples to the studied species. The network analysis showed that haplotype groups from a given location are more similar to those from the other locations than those from the same location (Fig. 2).

Fig. 2

A Mitochondrial haplotype network (COI, 658 bp). B Nuclear haplotype network (ITSs, 654 bp). Haplotypes are coloured according to sampling region and proportional to sample size. Branch lengths are proportional to sequence substitutions

A total of 92 C. caespitosa individuals were analysed for microsatellite loci variation, of which 68 (73.91%) had a unique multilocus genotype (Table 3). Identical multilocus genotypes were observed in KRY and LIO, which indicated a clonality percentage of 13.33% and 35.48%, respectively. Genetic evenness (Go/Ng) was 0.86 for KRY and 0.21 for LIO. The number of unique genotypes (Ng) ranged from 11 at LIO to 31 at NEA. Genotypic richness (Ng/N) was highest for NEA and lowest for LIO. Genotypic diversity (Go/Ge) was high for both KRY and NEA (0.75 and 1, respectively), indicating the dominance of “sexual” reproduction at these sites; for LIO, it was only 0.07, indicating “mostly asexual” reproduction. The D index values (Pielou 1969) were in line with the genotypic diversity results: D values for KRY and NEA were 0.98 and 1, respectively, and 0.58 for LIO. These values were also consistent with the clonal structure of the KRY and LIO populations, which showed that KRY had three clones (genets) with two, three and two specimens (ramets), respectively, whereas LIO had two genets with twenty and two ramets, respectively. No clonal individuals were detected within the Greek population (NEA).

Table 3

Genotypic diversity of Cladocora caespitosa at three eastern Mediterranean locations based on the analysis of 19 microsatellites

	KRY	LIO	NEA
N	30	31	31
Ng	26	11	31
Ng/N	0.86	0.35	1
Go	22.5	2.33	31
Go/Ng	0.86	0.21	1
Ge	30	31	31
Go/Ge	0.75	0.07	1
D	0.98	0.58	1

N = number of polyps (colonies) sampled; Ng = number of unique multilocus genotypes per site; Ng/N = genotypic richness; Go = observed genotypic diversity; Go/Ng = genotypic evenness; Ge = expected genotypic diversity (number of individuals genotyped per site); Go/Ge = genotypic diversity and the D index

Standardized allelic richness across loci for each population ranged from 3.10 for LIO to 4.65 for NEA (Table 4). The mean value across loci and populations was 3.91. A high level of heterozygosity was detected among all populations. Observed and expected heterozygosities ranged from 0.511 for KRY to 0.594 for LIO and from 0.494 for LIO to 0.578 for NEA, respectively. FIS values for KRY and NEA were positive, indicating a heterozygote deficit in both populations; by contrast, LIO showed a high and negative value, indicating a heterozygote excess in this population.

Table 4

Genetic diversity estimates for the three Cladocora caespitosa populations

	Na	Ho	He	FIS
KRY	3.986	0.511	0.557	0.1001*
LIO	3.105	0.594	0.494	− 0.1577*
NEA	4.655	0.546	0.578	0.0685*
Mean	3.915	0.550	0.543	0.0036

Na, standardized number of alleles; Ho, observed heterozygosity; He, expected heterozygosity; FIS, inbreeding coefficient; *populations that are not in Hardy–Weinberg equilibrium (HWE)

Significant deviation from HWE was observed in some of the loci, mainly for the Cypriot populations, particularly LIO. For the Aegean population, only one locus significantly deviated from HWE after Bonferroni correction. Following a correction for null alleles, only two of the loci showed no significant deviation, and no statistically significant differences were observed in pairwise FST and pairwise corrected FST values. Therefore, in order to avoid the artefact of shared alleles among populations, null alleles were not considered in subsequent analyses. The BAYESCAN analysis showed that no loci were under selection, consistent with the results of the ARLEQUIN tests of neutrality.

Population differentiation

The global value of FST revealed a significant level of genetic differentiation (FST global = 0.079, P = 0; standardized FST global = 0.192). Pairwise FST values ranged from 0.059 for KRY versus NEA to 0.104 for LIO versus NEA (standardized FST values ranged from 0.141 to 0.238; F’ST values ranged from 0.143 for KRY versus NEA to 0.236 for LIO versus NEA, similar to the standardized FST ones) (Table 5). Significant differentiation was observed between the two Cypriot populations, KRY and LIO (0.097), which are only separated by 11 km. A similar level of differentiation was observed between LIO and the Aegean population NEA (0.104), which are separated by around 1300 km. Therefore, no significant association between genetic differentiation (FST) and geographic distance across the sea was found (Mantel test, P = 0.516, R = 0.138).

Table 5

Lower diagonal pairwise FST values among Cladocora caespitosa populations. Upper diagonal shows the F’ST. All values are significant (P < 0.05)

	KRY	LIO	NEA
KRY	0	0.219	0.143
LIO	0.097	0	0.236
NEA	0.059	0.104	0

The STRUCTURE analyses, following the Evanno methodology based on ΔK values implemented in STRUCTURE Harvester and Clumpak, showed two genetically differentiated genetic clusters under the parameters “popinfo + location prior” (Fig. 3). However, the results from Clumpak based on probability (ln (Pr(X|K) values) showed K = 5 as the best value for the data in the analysis of K up to 5. In the analysis of K up to 10, K = 9 was the best value. Results for K = 2 showed that the individual assignment of specimens belonging to the Greek population (NEA) was classified in one of the clusters with a probability of at least 99%. Cypriot individuals from KRY were also assigned with high probabilities (> 99%), mostly to the same group as NEA individuals. The LIO specimens were equally split between the two groups detected. For K = 5, NEA specimens were assigned to one of the clusters with at least a 90% probability. Cypriot individuals from KRY were assigned to two other clusters and those from LIO to the remaining two clusters (except one individual that was assigned to one of the KRY clusters). For K = 9, the Greek specimens were equally divided into two clusters. KRY specimens were spread among five clusters, three of which only represent this population. Lastly, LIO specimens comprised two exclusive clusters, except one individual that coincided with one of the KRY clusters. The results of the StructureSelector analysis (which uses the estimators MedMeanK, MedMedK, MaxMeanK and MaxMedK) indicated a greater likelihood of three subpopulations within the clusters (Fig. 3; Puechmaille plots in Supplementary Appendix S1). In this case, individuals belonging to KRY were divided into the three clusters, those of LIO to the second and third clusters and, lastly, NEA specimens were assigned to the third cluster with a probability around 100%.

Fig. 3

Structure results for the three locations in the Levantine (KRY and LIO) and Aegean (NEA) seas when K = 2 following the Evanno methodology based on ΔK values, K = 3 based on the Puechmaille method, K = 5 and K = 9 based on probability (ln (Pr(X|K) values)

Results of the FCA on the genetic structure showed that the two axes explained 11.78% of the variability (6.48% and 5.30%, respectively). The results of the PCoA analysis showed that 66.06% of the variation could be explained by the first axis, which separated KRY and NEA from LIO, while the second axis (33.9% of the variation) mainly distinguished KRY from NEA (Fig. 4).

Fig. 4

Factorial correspondence and principle coordinates analyses. On the left, the results of the FCA showing that the two axes explain 11.78% of the variability (6.48% and 5.30%, respectively). On the right, results of the PCoA of FST values among populations. The first axis explains 66% of the variation, and the second, 34% of the variation

Finally, we assessed potential gene flow barriers, migration and bottlenecks. The results of the BARRIER analysis revealed two supported barriers (bootstrap value = 100): one separating LIO from KRY and the other separating NEA from LIO. In the GENECLASS2 analysis, three individuals, corresponding to only 4.41% of the total number of individuals studied, were identified as potential first-generation migrants. Two KRY individuals were considered migrants from LIO and NEA, and one LIO individual was considered a migrant from KRY.

Results from the BOTTLENECK analysis indicated a possible recent bottleneck in the LIO population: an excess of heterozygosity with significant values was found for all three mutation models under the Wilcoxon statistical test and only one (IAM) under the Sign test. LIO also showed a shifted mode. The other populations (KRY and NEA) showed significant values for some models but not for others, and both populations had a normal L-shaped allele frequency distribution as expected under mutation–drift equilibrium (Table 6).

Table 6

Heterozygosity excess in the analysed populations of Cladocora caespitosa

	Sign test			Wilcoxon test			Mode shift
	IAM	SMM	TPM	IAM	SMM	TPM
KRY	0.32849	0.20417	0.41081	0.02468	0.11266	0.36905	Normal
LIO	0.00034	0.09676	0.06991	0.00010	0.00750	0.00192	Shifted
NEA	0.10251	0.00400	0.55360	0.04776	0.99734	0.52351	Normal

The two statistical tests, Sign and Wilcoxon sign-rank, were conducted under three different mutation models: IAM, SMM and TPM. Significant values are in bold. The mode shift is also indicated

Discussion

Our microsatellite analysis revealed remarkable results related to microscale differentiation and reproductive mode differences between the two geographically close Cypriot populations, LIO and KRY. The LIO population appeared as a differentiated cluster, with the lowest mean number of alleles per locus and expected heterozygosity, the highest pairwise FST values, a negative inbreeding coefficient, indicating heterozygosity excess, and evidence of a recent bottleneck. Moreover, it is the only population to have an asexual index greater than 35% compared with the other populations analysed here and with those in previous studies of C. caespitosa (Casado-Amezúa et al. 2014; López-Márquez et al. 2019). However, its level of genetic diversity is similar to that of the other two populations, likely due to the contribution of some level of sexual reproduction, which can maintain genotypic variation, even at low levels (Bengtsson 2003).

The Cypriot populations, which are located in highly populated areas, have a fragmented distribution pattern along a shoreline that is characterized by a low number of localities with large reefs (Hadjioannou 2019). Abnormal warming events were reported for the area in 2012, 2014 and 2015 (Jimenez et al. 2016; Garrabou et al. 2019; Hadjioannou 2019), as were extreme storms (Hadjioannou et al. 2016). Damage related to the temperature anomaly was severe: at LIO, the mortality rate of the colonies was 93% during the summer of 2012 (Garrabou et al. 2019). In addition, LIO is more exposed to swells and, thus, to subsequent mechanical damage and is more affected by human-mediated impacts, such as the inflow of high concentrations of inorganic nutrients (Hadjioannou 2019).

We hypothesize that these environmental stressors, temperature and hydrodynamics, account for the high percentage of clones observed in the LIO population. Bleaching caused by heat waves and severe storms can weaken colonies, leading to their physical fragmentation. Depending on the extent of fragmentation and damage, a local decrease in the ratio of sexual to asexual reproduction can occur, leading to the dominance of clones (Aranceta-Garza et al. 2012; Baums et al. 2014). For instance, clonality in populations of the Caribbean reef-forming coral Montastrea annularis has been imputed to hurricane disturbance (Foster et al. 2013). In general, fragment survival, which is influenced mainly by environmental factors (Lirman 2000), determines the level of local clonality. Partial mortality of colony tissue can result in small clonal colonies, as was observed in LIO, where the size of many colonies decreased after bleaching due to partial mortality (Jimenez et al. 2016). Indeed, colony mortality was significantly higher at LIO than at KRY during 2015 (due to summer temperature anomalies; Garrabou et al. 2019). This event affected 100% of the monitored colonies at LIO (resulting in an overall loss of ~ 17% of pigmented tissue) and 45% of the colonies at KRY (~ 11% loss of pigmented tissue). Necrosis of pigmented tissue on coral colonies was also recorded in the summer of 2014, but to a lesser extent than in 2015 (Hadjioannou 2019). In addition, at LIO, corallites were much more brittle, easily detaching from the mother colony, and bioerosion was more evident by the abundance of coral rubble observed scattered on the seabed among the colonies (LH, personal observation). All of these observations support our hypothesis for the high level of clonality and the changes in the ratio of sexual to asexual reproduction found at LIO. Although natural and anthropogenic disturbances, such as thermal stress and increased nutrient availability, are known to compromise sexual reproduction (Harrison 2011), thus playing a crucial role in the population structure of organisms, an increase in asexual reproduction by polyp (or ameiotic planulae) budding cannot be ruled out as a direct response to environmental stress. A similar process has been observed in other hermatypic corals, such as Pocillopora verrucosa, which reproduces asexually in response to environmental factors including storms and increased bioerosion by predators such as echinoids (Aranceta-Garza et al. 2012).

Another interesting result from our study is the indication of LD in all loci for the Cypriot populations. LD was not found in any of the loci analysed for western Mediterranean (Casado-Amezúa et al. 2011) and Adriatic (López-Márquez et al. 2019) populations of C. caespitosa nor for the Greek population analysed in this study. Evolutionary factors such as genetic drift, mutation, gene flow, population subdivision, assortative mating or inbreeding could explain the LD observed for the Cypriot populations. Population size reductions can also increase LD due to the loss of some haplotypes after a bottleneck (Slatkin 2008). In fact, LIO can be classified as a bottlenecked population as evidenced by its shifted mode in our allele frequency distribution analyses, thus discriminating it from stable populations (Luikart 1997). Asexual or clonal reproduction can mimic physical linkage over the genome (de Meeûs and Balloux 2004). The Cypriot populations both present some degree of asexual reproduction given the number of clones identified in our analyses (7 and 22 clones out of the 30 and 31 samples analysed for KRY and LIO, respectively). This result, together with the low level of gene flow found among populations, could also contribute to the LD observed for these populations, matching the effects of inbreeding and/or assortative mating. Therefore, the sum of all impacts may be the cause of the LD. The inbreeding coefficients (FIS), which were negative for LIO, indicating a heterozygote excess, support this hypothesis. Given that none of the loci are under selection, other factors, including a small reproductive population size and asexual reproduction (among others) (Stoeckel et al. 2006), likely explain the negative FIS value. We propose that, under this context, the significant and negative FIS for LIO is due to its small reproductive population size and the effect of clonal reproduction on the number of heterozygotes.

We observed statistically significant differentiation among populations. The level of differentiation found between the two Cypriot populations (KRY and LIO) is comparable to the one found between LIO and the geographically distant Greek population (NEA). The lowest level of differentiation was found between NEA and KRY. Thus, as our IBD results show, geographic distance does not determine the connectivity patterns of C. caespitosa in the study area. Coastline configuration, seabed morphology and local hydrodynamic conditions (e.g. eddies, waves and counter currents) are locality specific factors that could prevent connectivity. Results from the study of putative barriers, which analyses geographic distance and genetic differentiation, support the existence of two main barriers, isolating LIO from both KRY and NEA.

The genetic pattern found in this study reflects one of isolation by environment in which the populations with similar environmental conditions separated by 1300 km (KRY and NEA) are genetically closer to each other (FST = 0.06) than are the adjacent populations (KRY and LIO) (FST = 0.11) separated by only 11 km that have different levels of environmental stress. The differential exposure of KRY (cove sheltered) versus LIO (exposed shoreline) populations to extreme weather events and anthropogenic activities likely contributes to the different genetic patterns observed at the two sampling localities (Jimenez et al. 2016; Hadjioannou 2019). Human impact on the environment may change the extent to which factors such as temperature, salinity, turbidity, nutrient levels or wave exposure build ecological boundaries. These boundaries can play an important role in the genetic structure of coastal populations with contrasting habitats (Tisthammer et al. 2020). They can also contribute to a decline of new recruitment in highly affected nearshore populations (Puritz and Toonen, 2011). According to Casado-Amezúa et al. (2014), self-recruitment predominates in this species. Conversely, genetic connectivity among distant populations, although rare, may be related to sporadic and stochastic dispersal events in which regional surface currents become linked during the spawning period.

Overall, our analyses have consistently shown that the Aegean Sea population comprises a single homogenous group. Also, none of the proposed clusters recognize individuals from each of the Cypriot population as a single group. Genetic structure analysis of KRY indicates that almost half of the population presents usual admixture, whereas the other half is constituted by individuals related to differentiated ancestors. For LIO, there is no evidence of a homogeneous origin for all individuals; rather, at least two groups of multilocus genotypes are found.

The genetic structure results of the analysed populations differed depending on the number of clusters specified or the criterion used. Best values for K varied from 2, 3 or 5 to 9, surpassing the number of locations analysed (except for K = 2), indicating contributions from foreign individuals or breaks in the populations such as those produced by bottlenecks, drift and clonal reproduction. In cases in which the premises of the equilibrium population model are violated, it may be difficult to find a unique or natural answer to what is a “correct” K value. Therefore, it usually makes sense to focus on K values that capture most of the structure from the data and that seem biologically sensible. Consequently, we conclude that we have a strong indication of real population structure, despite the different results of the structure analyses. The low levels of migration found support these results. Only two KRY individuals were identified as putative migrants, one from LIO and the other from NEA, and only one LIO individual was considered as a migrant from KRY. The Aegean population showed a lack of immigration.

Our results, taken together, indicate low connectivity among the three analysed populations. We observed not only significant and high levels of genetic differentiation among populations, even between the two populations separated by only 11 km, but also strong barriers, ecological boundaries from human impacts and low levels of migration. Connectivity is an important aspect of species replenishment (Saenz-Agudelo et al. 2011). Connectivity inferences assessed by genetic studies provide data on population resilience; therefore, they are considered a crucial factor in species conservation. We initially hypothesized that C. caespitosa would display a low level of connectivity due to larval behaviour. However, differentiation was much higher than expected and was not related to geographic distance. Moreover, the marginal distribution of these locations, along with recurrent episodes of heat waves and extreme storms and human impacts on the environment, likely accounts for the unusually high degree of asexual reproduction in the Cypriot populations, a trend that, to date, has not been reported in other Mediterranean areas. The easternmost populations of the endangered C. caespitosa are clearly threatened, and efficient measures, such as pollution control, are needed to preserve them, especially in LIO where contamination seems to worsen habitat conditions. Therefore, local adaptations have to be considered in conservation efforts in order to preserve the gene pool of this endangered species.

Below is the link to the Supplementary Information.

Supplementary material 1 (JPEG 102 kb)