Living upside down: patterns of red coral settlement in a cave

Study area and experimental design

The study was conducted inside the Colombara (or Marcante) cave (Lat 44°18′35″N, Long 9°10′37″E) located on the east coast of the Gulf of San Fruttuoso (Italy). The cave, 10 m long, stands from 34 to 39 m depth on a rocky cliff southward oriented. The cave walls host a rich assemblage of sessile invertebrates typical of the sublittoral caves of the northwest Mediterranean Sea (Morri et al., 1986) with many sponges, corals, bryozoans, polychaetes, and tunicates.

Field experiments were approved by the Marine Protected Area of Portofino and by the Università Politecnica delle Marche (Authorization n. 3/2011 (n. prot. 449/2-1-5.) and authorization no 4/2012 (Protoc. No 409/2-1-1)).

In June 2010, about three weeks before the start of red coral spawning (Santangelo et al., 2003), sixteen 20 × 20 cm white PVC tiles, drilled in the centre, were fixed inside the cave by steel screws. PVC was selected owing to the success in previous experiment on larval recruitment (C Cerrano, pers. comm., 2010; but see also Kennedy et al., 2017). Moreover, having positive buoyancy, the risk of detachment from the ceiling was avoided.

In order to test if recruitment is affected by orientation, one plot of four tiles was located on the right vertical wall, another plot on the left vertical wall, and two plots on the ceiling of the cave (Fig. 1). Each tile was attached to the rock, 1 m from the entrance and at a minimum distance of 30 cm from any another tiles, to avoid possible mutual interference. The red coral population distribution into the cave is patchy, showing an average density of 349 ± 215 col/m² (Cattaneo-Vietti, Bavestrello & Senes, 1993). Tiles were attached in low-density areas of the red coral population trying to limit as much as possible the breakage of surrounding colonies.

In February 2012, after two reproductive events (summer 2010 and summer 2011), the tiles were removed (n = 14 as two located on one vertical wall were lost, Fig. 1). Recovered tiles were fixed in 80% ethanol and transferred to the laboratory. A picture of each tile was taken with a NIKON camera on a stereomicroscope NIKON SMZ 1500 and analysed with IMAGE J software version 1.24o (Schneider, Rasband & Eliceiri, 2012) to study the spatial distribution of the red coral individuals. The position of each individual was marked and size (diameter) was measured by averaging the minimum and maximum width. All the individuals were then removed from the tiles and for each of them; the number of polyps was counted under the dissecting microscope. Polyps were removed from each individual and stored in plastic tubes with 80% ethanol at 4 °C for the subsequent DNA extraction.

Based on the literature on red coral early life stages (Bramanti et al., 2005) the age of each individual was estimated on the basis of its height: individuals with an encrusting button shape and colonies with height equal to zero were assigned to the cohort 2011 (hereafter recruits); individuals that developed in height forming a small branch were assigned to the cohort 2010 (hereafter juveniles) (Fig. 2).

To check if a sort of cave-effect affected the pattern of recruitment, in the same period an additional series of plots with four PVC tiles were deployed on a vertical cliff out of the cave (Punta del Faro) where the red coral population has the same range of density as that of the cave (425 ± 100 colonies × m²) (Bavestrello et al., 2015). Plots were fixed at 35, 45, 55 and 70 m depth.

Measuring growth performances

The size structure of red coral individuals was obtained by analysing the frequency distributions of the basal diameter (for both recruits and juveniles) and of colony height (only for juveniles). The correlations between the parameters (diameter vs. height; diameter vs. number of polyps; height vs. number of polyps) were analysed by Pearson’s coefficient. Moreover, we tested the temporal and spatial variations of the abundance of individuals by two separate one-way ANOVAs based on 9,999 permutations with Cohort (two levels; fixed; recruits and juveniles) and Position (two levels; fixed; ceiling and walls) as factors using PRIMER v6 software program (Clarke & Gorley, 2006).

Measuring genetic variability and structuring

Due to the small number of C. rubrum individuals found on the tiles deployed on the walls, the molecular analysis was carried out only on the individuals occurring on the tiles deployed on the ceiling of the cave.

Total genomic DNA was extracted from each individual (1 to 4 polyps per colony) following the CTAB protocol and purified by standard chloroform procedure. Seven microsatellite loci COR9, COR46, L7, COR58, MIC20, MIC24, MIC26 (Costantini & Abbiati, 2006; Ledoux et al., 2010b) were amplified either as single locus using the protocol by Costantini & Abbiati (2006) and by Ledoux et al. (2010b) or in multiplex using a QIAGEN Multiplex PCR Kit. Genotyping was performed by MACROGEN INC. Allele sizing was determined with Peak Scanner v1.0 software. Genotypic linkage disequilibrium at each pair of loci for each was tested using FSTAT v.2.9.3.2 (Goudet, 2001). Significance of each pairwise comparison was tested using 3,360 permutations of the data.

Scoring errors due to stuttering, large allele dropout and null alleles were controlled with MICROCHECKER version 2.2.3 (Van Oosterhout et al., 2004). Estimated frequencies of putative null alleles were subsequently calculated for each locus using the expectation maximization (EM) algorithm of Dempster, Laird & Rubin (1977) implemented in FREENA (Chapuis & Estoup, 2007). Red coral settlers sharing the same multilocus genotype (MLG) were identified using GENALEX version 6.1 (Peakall & Smouse, 2006). Identical MLGs can result from two distinct sexual reproduction events. To test this, the unbiased probability of identity (P_ID) that two sampled individuals shared identical MLG by chance (Kendall & Stewart, 1977) was computed. The total number of alleles (N_a), observed (H_o), and unbiased expected (H_e) heterozygosities (Nei, 1973) were calculated for each locus and for each cohort of settlers (either recruits and juveniles) using GENETIX software package version 4.03 (Belkhir et al., 2000). Allelic richness (Ar) was calculated after controlling for differences in sample size using a rarefaction approach implemented in HP rare with a common sample size of 54 settlers (Kalinowski, 2005). The f estimator of F_IS (Weir & Cockerham, 1984) was computed, and significant departures from the Hardy–Weinberg equilibrium were tested using Fisher’s exact test in GENEPOP version 3.4 (http://genepop.curtin.edu.au/; Raymond & Rousset, 1995), with the level of significance determined by a Markov chain randomization. Significant differences in genetic diversity (H_o, H_e and Ar) between recruits and juveniles were tested using a Student’s t-test. For all the analyses, when necessary, a correction for multiple tests was applied with a false discovery rate of 0.05 (Benjamini & Hochberg, 1995).

To determine whether individuals were more related than expected under panmixia, the R_XY pairwise relatedness coefficient (Queller & Goodnight, 1989) was computed separately for recruits and for juveniles. This index varies between 0 and 0.5 with R_XY = 0.5 for full-sib relationships, R_XY = 0.25 for half-sibs and R_XY = 0 for unrelated individuals in an infinitely large panmictic population (Peakall & Smouse, 2006). The observed mean and variance of R_XY were compared with their expected distribution under the null hypothesis of panmixia using 1,000 permutations of alleles as implemented in IDENTIX (Belkhir, Castric & Bonhomme, 2002). The null distribution was obtained by a conventional Monte Carlo resampling procedure, which randomly selected 10,000 genotypes without replacement and then recalculating the statistic. Null hypothesis could be rejected with a significance level of 5%, given that the observed value of the statistic was above the 95% level of the resampled statistics.

An exclusion test, performed in GENECLASS 2.0 (Piry et al., 2004), was used to test whether individuals found in one tile was more similar to the individuals settled on the same tile. First, the likelihood that an individual belonged to a particular tile was computed with a Bayesian criterion of Rannala & Mountain (1997). Then, this likelihood was compared to a distribution of likelihoods of 10,000 genotypes simulated from each candidate tile with a Monte Carlo algorithm. An individual was excluded from its tile when the probability of exclusion was greater than 95% (P or α ≤ 0.05). The second part of this Bayesian analysis utilized the probabilities that the individuals excluded from their sampling tile originated from one of the other tiles. Thus, individuals that were excluded from their sampling tile when P ≤ 0.05, were assigned to another tile when P ≥ 0.1.

The value of effective population size (Ne) for each cohort was inferred using the standard linkage disequilibrium method with Waples (2006) correction. The computations were performed with LDNe under the random-mating model, excluding rare alleles with frequencies of less than 0.02 and using the Jackknife option to estimate confidence intervals (Waples & Do, 2008; Waples & Do, 2010).

To account for unbalanced sample sizes between the recruits and juveniles, genotypic differentiation was tested with an exact test (Markov chain parameters: 1,000 dememorizations, followed by 1,000 batches of 1,000 iterations per batch), and the P value of the log-likelihood (G)-based exact test (Goudet et al., 1996) was estimated in GENEPOP.

A Bayesian approach implemented in the program STRUCTURE version 2.2 (Pritchard, Matthew & Donnelly, 2000; Falush, Stephens & Pritchard, 2007) was used to estimate the number of genetic clusters, K, within the entire sample (i.e., recruits + juveniles). Mean and variance of log likelihoods of the number of clusters for K = 1 to K = 10 were inferred by running structure ten times with 1,000,000 repetitions each (burn-in = 100,000 iterations) under the admixture ancestry model and the assumption of correlated allele frequencies among samples. Due to the presence of null alleles, the clustering analysis was conducted on the original data set, using the option of null alleles coded as recessive alleles described in Falush, Stephens & Pritchard (2007). Mean likelihoods of K from ten runs were plotted using STRUCTURE HARVESTER 0.56.3 (available at http://taylor0.biology.ucla.edu/struct_harvest/). Results of all runs were averaged in CLUMPAK server (Kopelman et al., 2015). Moreover, since the Structure analysis could be inflated by the HW disequilibrium, a discriminant analysis of principal components (DAPC), available in the Adegenet package (Jombart, Devillard & Balloux, 2010) for R (R Development Core Team, 2012), was also performed. This technique is designed for multivariate genetic data (multi-locus genotypes). It maximizes the variation between groups by first performing a principal component analysis (PCA) on pre-defined groups or populations and then uses the PCA factors as variables for a discriminant analysis (DA), thus ensuring their independence.

Spatial autocorrelation analyses among individuals were performed with SPAGEDI (Hardy & Vekemans, 2002). The Loiselle’s kinship coefficient (ϕij; Loiselle et al., 1995) was used since it is not dependent on Hardy–Weinberg (HW) equilibrium conditions (Hardy, 1999) and it has been proven to be very effective in determining spatial genetic structure (Vekemans & Hardy, 2004). Two different analyses were performed. To test the spatial autocorrelation within the cave, distance categories were set based on the distance between tiles and considering the distance between two individuals found within the same tile as zero. Then, spatial autocorrelation analyses within those tiles containing more individuals (T4, T5, T7, T8, T9 and T10; see results) were carried out. Taking into account the distance between individuals, the distance categories were set in such a way that the number of pairwise comparisons within each distance category was approximately constant. Statistical significance was based on permuting individual locations among all individuals 10,000 times and calculating upper and lower 95% confidence interval for each distance class.

Red coral recruitment

In the plots positioned into the cave, 372 individuals were observed on the 14 tiles. Corallium rubrum settled on every tile deployed on the ceiling of the cave, but recruitment failed on three out of eight tiles on the walls (T1, T13, T14). In fact, of the 372 individuals, 350 were found on the ceiling tiles, and only 22 on the tiles deployed on the walls, corresponding to a significantly higher density on the ceiling than on walls (ANOVA: F_1,14 = 10.78; P < 0.01, Fig. 3). Of the 350 individuals on the ceiling, 278 were recruits and 72 were juveniles, corresponding to a density of 34.75 ± 23.86 /400 cm² and of 9 ± 6.82 / 400 cm², respectively (ANOVA: F_1,26 = 5.99; P < 0.05, Fig. 2). The 22 individuals found on the wall were all recruits (Fig. 3). On the tiles positioned along the vertical cliff, at all depths no recruits were found when checked in summer 2011.

Growth performances

The size/recruit structure showed a monotonic and decreasing pattern, in which recruits in the first class (recruits with diameter < 1.5 mm) represented the dominant class (Fig. 4A). Only three recruits (1%) had a diameter >6 mm, and they might have been formed by merging of two or more planulae (Fig. 2) as observed by C Cerrano (pers. comm., 2010). The number of polyps per recruit ranged from 1 to 22, with 70.66% of the recruits presenting 1–2 polyps, 22% between 3–4 polyps and 7.33% more than four polyps (mean  ± SD = 2.3 ± 2.1 polyps/recruit). Pearson’s coefficient showed a low correlation between these two variables (r = 0.4) (Fig. 4C).

The size structure of juveniles showed a more variable trend, with diameter values ranging from 0.15 mm to 8.95 mm with an average value of 2.5 ± 1.5 mm (Fig. 5A). The number of polyps ranged from 0 to 22 (average number of 9.21 ± 4.6) (Fig. 5C). Juveniles with 3–4 polyps were 9.7% of the total number, while 90.3% had more than five polyps. The height of the juveniles ranged from 0.5 to 9.15 mm, with a mean of 3.3 ± 1.8 mm (Fig. 5E). Pearson’s correlation coefficients showed no correlation between diameter and number of polyps (r = 0.09, Fig. 5B), nor between diameter and height (r = 0.03, Fig. 5F). A weak correlation between height and number of polyps was observed (r = 0.3, Fig. 5D).

Genetic variability

Individuals in which more then two loci did not amplify due to technical failures (e.g., low DNA quantity, no amplification of after re-amplification) were not included in the dataset. A total of 290 red coral individuals were genotyped. No genotypic disequilibrium was observed between loci (all P > 0.05 after FDR correction). No evidence of scoring errors due to stuttering or large allele dropout was found, according to MICROCHECKER. An excess of homozygotes was detected at all loci due to the presence of null alleles. Null allele frequencies ranged from 0.09 (Mic24) to 0.25 (Mic26 and Cor58). Number of alleles ranged from 10 (Mic20) to 23 (Cor9 and Mic26). The expected and observed heterozygosity ranged from 0.29 (Cor46) to 0.81(Cor9) and from 0.13 (Mic26) to 0.64 (Mic20), respectively. The estimators of F_IS (f) were positive and ranged from 0.21 (Mic20) to 0.86 (Cor46) (Table 1).

Overall, a low genetic variability (H_e = 0.59 ± 0.16 and H_o = 0.29 ± 0.20; mean ± SD) was found. Out of the 290 individuals analysed, 281 different multilocus genotypes (MLGs) were identified. Four MLGs were found twice; one MLG was encountered three times and one four times. Individuals sharing the same MLG always came from the same tile. Tiles where identical MLGs were found were T7, T8 and T10. The individuals sharing the same MLG were between 0.2 cm and 13.65 cm apart (in T7 and T10, respectively). In T8 and T10, one recruit and one juvenile shared a MLG. The unbiased probability of identity (P_ID) that two sampled individuals share identical MLG by chance was 1.5e⁻⁰⁴.

Juveniles and recruits did not show significant differences in terms of H_e and H_o(P values associated with the permutation procedure: P_Ho = 0.65 and P_He = 0.71); while recruits showed a significantly higher allelic richness than juveniles (P_Ar = 0.001).

Relatedness

A high degree of genetic relatedness among individuals was found. The mean observed relatedness was significantly different from the expected relatedness under panmixia (observed mean rxy = 0.027, resampled mean rxy = 0.014, t =1,173, P = 0.001). Pairwise relatedness based on the Queller & Goodnight (1989) coefficient revealed that 25.98% of the pairwise comparisons were involved in one or more parentage relationships, with 8.57% and 17.41% of individuals involved in a full-sib and half-sib relationship, respectively. These percentages were numerically similar considering the relatedness within the two temporal cohorts separately.

The percentage of individuals correctly assigned at the same tile where they were found by the individual assignments method using GENECLASS was around 50%. The effective population size (Ne) was estimated as 68.7 (95% CI [42.5–116.4]) including all the individuals and as 32.5 for recruits (95% CI [20.6–57.8]) and −1,467.8 for juveniles (95% CI [87.1–∞]). The last negative value is expected when the population is sufficiently large that no notable linkage disequilibrium is induced through genetic drift (Waples & Do, 2010).

Spatial and temporal genetic structure

No genetic structuring was observed between recruits and juveniles (F_ST = 0.008, P = 0.16). The clustering method identified K = 2 gene pools based on Evanno’s delta K statistic (Figs. S1, S2) but with high levels of admixture as many individuals could not be undoubtedly assigned to either cluster. Nevertheless, an eastward genetic gradient of differentiation were observed. Regarding the results from the DAPC analysis, the plot of BIC as a function of the number of clusters K (ranging from one to 16) did not present a real minimum (Fig. S3). Nevertheless, when using tiles as groups, the DAPC was in agreement with the Structure results (Fig. 6) with a higher genetic isolation of the titles T3 and T4 with all the other tiles.

Significantly positive kinship coefficients were detected between all the individuals within the single tile (ϕij = 0.066, p < 0.001) indicating that individuals within tiles had a higher genetic relatedness than random pairs of individuals. Within the cave, the autocorrelogram suggested an estimated patch size of less than 40 cm (Fig. 7). The spatial autocorrelation between individuals within each tile showed a significant positive value within a range from 3 to 10 cm.