Introduction

Black corals (Cnidaria: Anthozoa: Hexacorallia: Antipatharia) in the Mediterranean Sea are found within mesophotic assemblages together with gorgonians (Bo et al. 2009a, 2019; Grinyó et al. 2016), forming complex multi-species three-dimensional forests (Gori et al. 2017). These highly diverse ecosystems are usually found between 60 and 200 m depth (Bo et al. 2015, 2019; Grinyó et al. 2016; Ingrassia et al. 2016; Gori et al. 2017), with the greatest diversity and abundance observed along shelf edge environments (100–200 m) (Gori et al. 2017). This is most likely a result of a mid-domain effect (Colwell and Lees 2000; Grinyó et al. 2016) where shallow and deep species with overlapping distributions merge into mixed assemblages. Importantly, however, the abundance of Mediterranean black corals increases with depth, and they become more prevalent compared to other coral species in deeper hard-bottom environments within cold-water coral assemblages, including white coral ecosystems and offshore seamount summits (Bo et al. 2011, 2020; Grinyó et al. 2018; Chimienti et al. 2019).

Due to their arborescent morphology and complex three-dimensional structure, antipatharians serve as important ecosystem engineers (sensu Jones et al. 1994). As such, they enhance the spatial heterogeneity of the seabed and attract highly diverse associated fauna (Bo et al. 2012, 2015), including numerous species of commercial interest, such as the lobsters Palinurus elephas and Homarus gammarus and the fish Conger conger (Bo et al. 2009a; Chimienti et al. 2020). Furthermore, black corals play an important ecological role in the benthic-pelagic coupling by capturing plankton and particulate organic matter suspended in the water (Gili and Coma 1998; Rossi et al. 2017; Coppari et al. 2020a).

Seven species of black corals, belonging to five families and five genera, are currently known to inhabit the Mediterranean Sea: (1) Antipathes dichotoma, Antipathidae, (2) Leiopathes glaberrima, Leiopathidae (3) Antipathella subpinnata and A. wollastoni, Myriopathidae, (4) Parantipathes larix and Parantipathes sp., Schizopathidae and (5) Phanopathes sp., Aphanipathidae (Bo et al. 2018, 2020; Bo and Bavestrello 2019; Corbera et al. 2019).

Among Mediterranean black corals, A. subpinnata (Ellis and Solander 1786) is the most frequently observed (Bo et al. 2018; Bo and Bavestrello 2019). It typically inhabits the mesophotic zone with a peak of abundance at depths between 60 and 150 m (Bo et al. 2009a), although it can be found down to 600 m depth on hard-bottom substrate (Deidun et al. 2015).

Black corals are accidentally targeted by artisanal and recreational fisheries because the rocky elevations on which they thrive often represent fishing grounds (Bo et al. 2009a, 2019; Chimienti et al. 2020). Due to their arborescent morphology and erect habitus, black corals are among the most frequently found coral species in the fishing bycatch (Mytilineou et al. 2014; Deidun et al. 2015). Fishing impact can cause the direct removal or partial damage to the colonies, which consequently may lead to overgrowth of the skeleton by various fast-growing organisms (Mortensen et al. 2005; Bo et al. 2014). These impacts can have far-reaching and long-lasting effects on population dynamics of Mediterranean black corals, especially when considering their low growth rates (Bo et al. 2015).

For these reasons, many black coral species are categorized as ‘threatened’ by the International Union for Conservation of Nature (IUCN) Red List of Mediterranean Anthozoa (Otero et al. 2017). A. subpinnata in particular is considered as ‘Near Threatened’ by IUCN. This antipatharian is also included in the Annex II of a Barcelona Convention protocol for the Specially Protected Area of Mediterranean Interest (SPAMI) and listed as a representative species that contributes to the formation of hard-bottom coral gardens by the International Council for the Exploration of the Sea (ICES). These coral gardens are considered as Vulnerable Marine Ecosystems (VMEs) by ICES and, due to their low resilience to mechanical impacts, should receive priority of protection (FAO 2009).

Genetic connectivity provides a valuable insight into the processes of population maintenance and replenishment following environmental disturbance (Hastings and Botsford 2006; Bernhardt and Leslie 2013) and is often used as a proxy for population resilience. There are numerous biological factors that affect population connectivity patterns in benthic marine species (e.g., cold-water corals), such as (1) sexual (broadcast spawning vs. brooding) and asexual reproductive strategies, (2) the pelagic larval duration (PLD) of the offspring, (3) larval buoyancy, (4) motility behavior, and (5) life-stage mortality (Cowen et al. 2006; Treml et al. 2012, 2015; Costantini et al. 2018). For instance, populations of broadcast spawning coral species have been documented to exhibit higher admixture rates compared to coral species with internal fertilization (Thomas et al. 2019).

Isolated coral populations are particularly vulnerable to increased disturbance since the negligible supply of coral propagules from adjacent sites is likely to decrease their genetic diversity and consequently reduce their resilience and capacity to adapt to environmental perturbations (Costantini et al. 2017; Thomas et al. 2019). Therefore, quantifying genetic variability between coastal and offshore populations and their levels of connectivity will provide valuable insight into which populations can serve as a potential reservoir of genetic diversity to others and could prove critical for effective policy-making and resource management (Turner et al. 2019). Currently, the population genetics of Mediterranean black corals has not been investigated, due to field sampling constraints to collect colonies at mesophotic depths and the lack of species-specific genetic markers with high resolutive power (Costantini et al. 2017). Restriction site-associated DNA sequencing (RAD-Seq) methods (Baird et al. 2008) are particularly valuable to study these non-model species because the large number of single nucleotide polymorphisms (SNPs) obtained with RAD-Seq allows for more robust inferences of phylogenetic and population genomic structuring compared to traditional markers, even in biological data obtained from small sample sizes (Schopen et al. 2008; Reitzel et al. 2013; Everett et al. 2016; Hodel et al. 2017).

Here, a 2bRAD genotyping approach (Wang et al. 2012) was used to genotype A. subpinnata individuals from lower mesophotic (offshore) and upper mesophotic (coastal) populations to detect current genetic diversity and to identify geographic and bathymetric barriers to gene flow. RAD-Seq data were integrated with the environmental features of the sites and reproductive behaviour of the species to understand their influence on the observed genetic pattern. The main objective of our study was to evaluate the genetic connectivity within and among Mediterranean populations of A. subpinnata. This paper highlights the value and vulnerability of Mediterranean A. subpinnata underwater forests, as well as the importance of using genomic techniques to underpin informed management and conservation measures, particularly for remote deep-sea ecosystems where the collection of population connectivity data with different approaches (i.e., population demography, marking and recapture of larvae, tracking of larvae using Lagrangian modeling) is more challenging (Miller and Gunasekera 2017).

Material and methods

Sample collection

All the sampling activities were carried out between April 2017 and July 2018 (Fig. 1; Table 1). Six A. subpinnata populations, located along the Italian coastline at different depths (ranging from 55 to 210 m), were targeted for this study. The sampling locations included one Ligurian seamount (Santa Lucia; 164–210 m depth), one Tyrrhenian canyon (Posada; 152–156 m depth), and four coastal populations (55–75 m depth): three from Liguria (from East to West, Portofino, Savona, Bordighera, distanced apart 34 NM and 50 NM, respectively) and one from Sicily (Favignana). Colonies from the two deep populations were collected using a ROV MultiPluto grabber (GayMarine), whilst coastal populations were sampled using technical divers, for a total of 55 samples of A. subpinnata (Table 1). A detailed description of the sites can be found in Supplementary Material (file ESM1.docx of the Supplementary Material).

Fig. 1
figure 1

Field sampling. In situ pictures of the sampling operations by means of the ROV MultiPluto grabber (a, b) and Antipathella subpinnata forests (c, Bordighera, d, Santa Lucia Seamount)

Full size image
Table 1 Sampling localities for Antipathella subpinnata including code, date of sampling, geographic coordinates (latitude and longitude), depth: shown as categories (coastal/offshore) and expressed in meters, average number of reads per population (± SD), number of colonies sampled in each sampling locality (N. col), number of 2bRAD genotyped individuals (Ind) and unique multi-locus genotypes (MLG), observed and expected heterozygosity values (HO and HE), and inbreeding coefficient values (FIS)
Full size table

Laboratory work

DNA extraction

Out of the 55 sampled colonies, intact polyps were subsampled successfully from 49 A. subpinnata individuals and immediately preserved in 95% ethanol (EtOH) until genetic analysis. Total genomic DNA was extracted using the MinElute PCR Purification Kit (QIAGEN) following the manufacturer’s protocol. This kit is optimized for the extraction of degraded DNA (i.e., ancient DNA, see Angelici et al. 2019) and has yielded higher quantities of DNA compared to the other two extraction methods tested: (1) EUROGOLD Tissue-DNA Mini Kit (EuroClone) and (2) CTAB (cetyl trimethylammonium bromide)-based DNA extraction. Quality of extracted high molecular weight gDNA was visualized on a 1% agarose gel stained with GelRed (Invitrogen).

2bRAD library preparation

Restriction-site associated DNA sequencing libraries were constructed with the 2bRAD genotyping approach (Wang et al. 2012) using protocols outlined in Boscari et al. (2019) to obtain genome-wide SNPs. (More information can be found in ESM2.docx file of the Supplementary Material). Obtained 2bRAD tags were sequenced and demultiplexed by Genomix4Life Srl (Baronissi, Salerno, Italy) on an Illumina NextSeq500 sequencing platform (Illumina, San Diego, CA) using a single-end 75-bp high-output flow cell. To control for potential erroneous variants from library preparation and sequencing, eight specimens were sequenced in duplicate using different barcodes to serve as technical replicates for comparison of libraries. This resulted in a total of 57 fastq files (49 samples + 8 technical replicates) for downstream analysis.

Bioinformatics pipeline and data analysis

Quality filtering and trimming of raw reads

Demultiplexed raw reads were initially checked for various sequence quality metrics using FastQC (v0.11.5) (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) (01_FastQC.sh; script in Supplementary Material). Ligation adapters, CspCI restriction enzyme recognition sites, and low-quality basecalls (Phred score < 20) were trimmed to obtain high-quality 2bRAD fragments of uniform length (32 bp), using the Eli Meyer pipeline (https://github.com/Eli-Meyer/2bRAD_utilities) and as done in Boscari et al. (2019).

SNP calling in stacks v2

SNP genotyping was performed in the Stacks v2 software (version 2.3e) (Rochette et al. 2019). After excluding one technical replicate (by random) from each of the samples sequenced in duplicate, and low coverage samples (i.e., less than 100,000 quality-filtered 2bRAD tags), SNP calling was performed on remaining samples using the denovo_map.pl pipeline, due to a lack of a reference genome for A. subpinnata. Stacks loci assembly parameters were set to 3 (for -m, -M, and -n), and the obtained SNPs were exported in genind format using the following options from the Stacks v2 Populations program: (1) −R = 0.8 (to retain loci shared across at least 80% of individuals), (2) –write_single_snp (to keep a maximum of two SNPs per each variant site), (3) -min-maf = 0.01 (to exclude lowly abundant alleles), and (4) -min-mac = 2 (to exclude non-informative monomorphic loci). For more information, refer to the ESM3.docx file and the 4.tests_denovo.sh and 8.genepop_analyses.sh bash scripts in Supplementary Material.

Statistical analysis

SNP data in genind format were imported into R version 3.6.3 (R Core Team 2018) with the R package ‘adegenet’ (Jombart 2008) and additionally filtered to exclude: (1) loci out of Hardy–Weinberg equilibrium (HWE) in ‘pegas’ (Paradis 2010), (2) loci under selection (outliers) in ‘fsthet’ (Flanagan and Jones 2017), (3) individuals with over 30% of missing loci, and (4) any remaining missing loci, which were replaced with mean allele frequency, in ‘poppr’ (Kamvar et al. 2014, 2015). All downstream analysis was performed on two datasets, using (1) all loci and (2) neutral loci only.

Overall population statistics parameters, such as the inbreeding coefficient and observed and expected heterozygosity values, were computed at population level in ‘hierfstat’ (Goudet 2005) to infer if populations are partially clonal. Additional tests for clonal reproduction were carried out in ‘poppr’ to (1) assess whether certain individuals share the same Multi-Locus Genotype (MLG) based on a dissimilarity threshold of 3% and (2) by computing Agapow and Burt index of association (r̄d) (Agapow and Burt 2001) to investigate whether populations exhibit linkage across loci, as a potential indicator that asexual reproduction prevails. The 3% similarity threshold was chosen by assessing at which value the technical replicates start clustering into the same MLG—at values inferior to 3%, the replicates were still being clustered into separate MLGs.

Prior to ‘LEA’ and STRUCTURE analyses, the R genind object was converted into structure format using the genind2structure R function (Jenkins et al. 2019), and non-informative monomorphic loci were removed. The true number of ancestral populations was detected using (1) Bayesian clustering with the STRUCTURE program (Pritchard et al. 2000) with parameters set as in Carreras et al. (2019) and (2) ‘LEA’ R package (Frichot and François 2015). Admixture coefficients were estimated in ‘LEA’ at individual and population levels with the sparse nonnegative matrix factorization (SNMF) algorithm, as in Jenkins et al. (2019). Population structure was also explored through assignment of each individual to its predetermined locality a priori by conducting a discriminant analysis of principal components (DAPC) in ‘adegenet’.

FST statistics analyses included computation of (1) pairwise Nei's FST values (Nei 1987) for each population pair using 999 permutations and with the Benjamini & Yekutieli False Discovery Rate (B-Y FDR) correction to calculate the respective p values, as done in Carreras et al. (2019) and (2) Goudet’s G statistics test (based on Nei’s (Nei 1987) FST values) with 999 permutations (Monte–Carlo test). Both analyses were performed in ‘hierfstat’. This allowed us to compute FST values for each of the population pairs and to test whether factors ‘Populations’ and ‘Depth’ have a significant effect on overall genetic structuring, respectively.

Additional analyses to explore spatial patterns in genetic subdivision were conducted using a hierarchical design with two factors: (1) ‘Depth’ (two levels: offshore, coastal); and (2) ‘Populations’ (six levels: BORD, FAV, POS, PTF, SAV, SLU) nested within Depth. These analyses were deployed to test the effect of ‘Depth’ and ‘Population’ on genetic structuring, and they included: (1) principal components analysis (PCA) with ‘vegan’ (Dixon 2003; Oksanen et al. 2019), (2) pairwise permutational multivariate analysis of variance (PERMANOVA) using ‘pairwise.adonis’ R wrapper function (Martinez Arbizu 2019), and (3) a hierarchical analysis of molecular variance (AMOVA) in ‘poppr’, with 999 permutations computed in ‘ade4′ (Bougeard and Dray 2018).

A more detailed data analysis protocol can be found in the ESM3.docx file of Supplementary Material.

Results

Locus identification

Sequencing yielded 5,480,732 ± 6,608,522 (mean ± SD) raw reads (75 bp) per sample, ranging from 32,488 to 40,385,714 sequences per individual. This large variation in read depth could have been the culprit of the initial DNA quality differing between samples. After discarding sequences with adapters, CspCI recognition sites, and low-quality bases (Phred score < 20), 4,692,400 ± 6,227,770 high-quality 2bRAD tags (32 bp) remained per sample (Table S1: row ‘AVERAGE’ from column ‘Number of raw reads’). Prior to SNP calling, one technical replicate was randomly excluded for each of the eight replicated samples.

In addition, seven individuals with less than 100,000 good-quality 2bRAD tags were also discarded. Excluded technical replicates and low coverage individuals are indicated with TR and LC labels in (Table S1: column ‘Retained individuals’), respectively. The remaining 42 individuals were processed for SNP calling, with an average of 4,315,612 ± 6,350,346 reads per sample (Table S1: row ‘AVERAGE’ from column ‘Number of filtered 2bRAD tags’).

Six hundred and nineteen SNPs were obtained for downstream analysis. Additional five samples containing > 30% of missing loci were removed from the dataset (marked as ‘MD’ in Table S1: column ‘Retained individuals’), and remaining missing loci were replaced with mean allele frequency, which gave the final number of 37 retained individuals. None of the variant loci showed significant departures from HWE after any of the six correction methods for multiple comparisons. Conversely, each of the four outlier loci detection methods recorded a different number of outliers, equal to (1) 39 (Wright’s FST, Wright 1943), (2) 37 (the β based on the variance in allele frequencies, Cockerham and Weir 1993), (3) 75 (Cockerham and Weir’s sample size corrected β (betahat), Cockerham and Weir 1993), and (4) 68 outliers (Weir’s θ (theta) FST estimator, Weir, 1990) (Fig. S1; Table S3). Neutral SNPs were selected by removing 68 loci that were recorded as outliers in two or more outlier loci detection methods, similarly to what was done in Carreras et al. (2019) (Fig. S1; Table S3). All downstream analysis was performed using two datasets, (1) all loci (619 SNPs) and (2) neutral loci only (551 SNPs), in both cases shared across 37 individuals. Results were highly similar with both datasets, and only those obtained with neutral loci are shown in the main text, while those obtained using the ‘all loci’ dataset can be found in the ESM4.docx file of the Supplementary Material.

Population structure and connectivity

The computed inbreeding coefficient (FIS) values were positive in all localities except SLU, and heterozygote deficiency (HO < HE) was recorded in all populations (Table 1). In addition, four individuals were found to form two unique multi-locus genotypes (MLGs): (1) MLG 18 was composed of individuals BORD_7 and PTF_4, and (2) MLG 20 contained clonal SAV_8r and SAV_9 samples (Fig. S2). The linkage disequilibrium analysis recorded significant support (p = 0.001) for the hypothesis that alleles are linked across loci, with the observed Agapow and Burt index of association (r̄d) being outside of the distribution expected under no linkage in each of the sampling localities (Fig. S3). The null hypothesis of no-linkage among loci was therefore rejected. Both ‘LEA’ and STRUCTURE analyses were conducted using structure files with polymorphic loci only, which contained 594 (‘all loci’ dataset) and 551 (‘neutral loci’ dataset) biallelic SNPs. ‘LEA’ admixture analysis detected two subpopulation clusters (K) (Fig. 2). All individuals sampled at offshore localities were associated with Cluster 1 (blue-gray) with high membership assignment scores, while the majority of samples in coastal populations were more likely to originate from Cluster 2 (orange), with the exception of five samples (FAV_9; PTF_1, 2, 5, 6) that were assigned to the ‘offshore’ Cluster 1 (Fig. 2b). Bayesian clustering in STRUCTURE also inferred K = 2 (for both datasets) as the true number of ancestral populations based on the Evanno method (Fig. S4b, ESM4.docx), while K = 7 (‘neutral loci’ dataset) and K = 5 (‘all loci’ dataset) also seemed like possible scenarios according to mean LnP(k) values (Fig. S4c, ESM4.docx). In StructureSelector, Puechmaille method (Puechmaille 2016) detected that either 3 or 4 is the 'true' number of Ks (for both datasets).

Fig. 2
figure 2

(adapted from Boero et al. (2019), Fig. 7, MAW: Modified Atlantic Water), and the compass arrow (in the top left corner) points North. b Individual admixture proportions were presented as bar plots, with individuals being clustered within their originating sampling locality. Sampling localities were ordered according to geographical proximity. Offshore populations (Santa Lucia—SLU; Posada—POS) are emphasized with the black rectangle both on the map and the bar chart

Admixture analysis. Admixture proportions were computed in LEA R package at population and individual level. Two clusters were detected as the true number of Ks (ancestral populations) and are colored in blue-gray (Cluster 1) and orange (Cluster 2). a Mean population admixture proportions were presented as pie charts, and integrated onto a map using geographical coordinates of the sampling localities. Surface ocean currents are shown as black arrows

Full size image

When using the predefined clusters based on geographical sampling location, the multivariate clustering analysis (DAPC) plot showed clear genetic differentiation between offshore and coastal populations, with the SLU population being particularly isolated from other sampling localities. Furthermore, the offshore populations (SLU and POS) also grouped separately from one another, while DAPC clustering between all coastal localities (BORD, FAV, SAV, PTF) overlapped, with slight differentiation between BORD + FAV and SAV + PTF groups (Fig. 3a). Seven principal components (PCs) were retained to perform the DAPC analysis, as suggested by the calculated alpha (α) and cross-validation scores (Fig. S5).

Fig. 3
figure 3

Population structuring analyses of 37 individuals of Antipathella subpinnata sampled across six localities along the Italian coast at 55–210 m depth range, using 551 SNPs. Sampling localities are as follows: Coastal—(1) BORD: Bordighera; (2) FAV: Favignana; (3) PTF: Portofino; (4) SAV: Savona, and offshore—(5) POS: Posada; (6) SLU: Santa Lucia. a Discriminant analysis of principal components (DAPC) was performed by assigning individual samples to their originating sampling localities a priori. b Stacked bar plots show posterior assignment probabilities of individuals to their predetermined locations and were constructed using previously calculated DAPC population membership assignments. c A heatmap was used to visualize Nei’s FST pairwise distances for each population pair, which are shown in the lower (colored) triangle. Respective p-values (computed with 999 permutations, and B-Y corrected) are shown in the upper triangle. Significant B-Y adjusted p-values (< 0.05) were colored in red, and their corresponding Nei’s FST values from the lower colored triangle were marked in bold to show population pairs that significantly differ. d Goudet’s G statistics with 999 permutations was deployed to test the effect of Populations on genetic structuring

Full size image

Similar to the population structuring detected from the DAPC analysis, the composite bar plot also showed that coastal populations exhibit admixture, based on population membership probabilities computed in DAPC. Furthermore, two individuals from coastal localities (PTF_6 and SAV_2) seemed to originate from Posada, an offshore sampling site. In contrast, all individuals from offshore populations were assigned to their original localities with high membership assignment scores (Fig. 3b).

Pairwise FST values ranged between 0.11–0.19 and 0.10–0.17 for SLU/coastal and POS/coastal comparisons, respectively. Corresponding p-values were significant for all offshore/coastal population comparisons, ranging between 0.012 and 0.033. Coastal populations did not display a high level of genetic differentiation between localities (pairwise FST values = 0.0006–0.06), and significant differentiation was only detected between BORD-PTF populations (FST = 0.06, adjusted p = 0.05). Interestingly, the highest FST value was observed between offshore populations (SLU vs. POS: FST = 0.22, adjusted p = 0.012) (Fig. 3c). Goudet’s G statistics test also showed a significant effect of Population (p = 0.001, Fig. 3d) and Depth (p = 0.002, not shown) on overall genetic structuring.

PCA plot showed similar clustering patterns as the DAPC analysis, with the SLU population diverging the most from other sampling localities (Fig. 4a). Separate PCA clustering was also detected when all offshore and coastal localities were grouped together (Fig. 4b). PERMANOVA computed with factor ‘Depth’ showed significant structuring when all offshore and coastal localities were compared (Holm adjusted p = 0.002), although with low R squared value (R2 = 0.088), which means that, while statistically significant, ‘Depth’ does not account for a large portion of the dissimilarity between offshore and coastal populations (Table 2, Factor 1: Depth). Pairwise PERMANOVA computed for each population pair also showed significant genetic differentiation for offshore/coastal contrasts. In particular, this included two population pairs: (1) BORD vs. SLU and (2) SAV vs. SLU, with Holm (Holm 1979) adjusted p-values being equal to 0.03 for both pairwise comparisons (Table 2, Factor 2: Populations).

Fig. 4
figure 4

PCA plot. Population (a) and Depth (b) were used as PCA grouping factors. Sampling localities are as follows: Coastal—(1) BORD: Bordighera; (2) FAV: Favignana; (3) PTF: Portofino; (4) SAV: Savona, and offshore—(5) POS: Posada; (6) SLU: Santa Lucia

Full size image
Table 2 Permutational multivariate analysis of variance: PERMANOVA
Full size table

Finally, the hierarchical AMOVA showed significant genetic differentiation at all three levels of variation: (1) within Populations (p = 0.001), (2) between Populations, nested within Depth (p = 0.002) and (3) between Depths (p = 0.001) (Table S2; AMOVA: Testing significance with 999 Monte Carlo tests), given that the observed phi (ϕ) did not fall within the distribution expected from 999 permutations in either of those cases (Fig. S6). Importantly however, most variance (~ 71.03%) arose from within Populations, while only ~ 15.49% and ~ 13.48% of variance derived from between Populations (nested within Depths) and between Depths, respectively (Table S2; AMOVA: Components of covariance).

Discussion

Genomic data retrieved from the 37 A. subpinnata individuals revealed differing levels of connectivity among offshore (lower mesophotic) and coastal (upper mesophotic) populations in the Tyrrhenian and Ligurian Sea basins. The four coastal populations exhibited a high degree of admixture, suggesting that gene flow within these coastal regions is sufficient to maintain connectivity. In contrast, the two deeper, offshore populations (Posada and Santa Lucia) exhibited significant genetic differentiation from each other, indicating limited or no ongoing gene flow between these geographic regions (170 NM apart). Interestingly, differences among populations of the Mediterranean black coral were also clearly evident between individuals sampled from upper mesophotic, coastal habitats and the lower mesophotic, offshore localities, despite the close geographic proximity in some cases (e.g., Santa Lucia and Portofino, about 45 NM apart).

The patterns of connectivity observed in shallow populations along the Ligurian coast (Fig. 2a) is not surprising considering the strong current which flows westward along the shelf break (Millot 1999; Pinardi et al. 2015), which can influence gene flow by increasing the dispersal distance of coral propagules (Padron et al. 2018). However, physical characteristics of the environment (i.e., ocean currents and habitat suitability) do not solely explain the observed results, as the ecological traits of the species (i.e., organism life-history) are also important in determining population structuring and as such need to be considered to fully comprehend population connectivity (Werner et al. 2007). Numerous biological traits influence population connectivity in anthozoans: the asexual or sexual (brooders and broadcast spawners) mode of reproduction (Thomas et al. 2019), as well as the larval ecology (Cowen et al. 2006; Treml et al. 2012, 2015).

As with other coral species, black corals are known to reproduce asexually by budding new polyps for the colony growth (Pax et al. 1987), but antipatharians are also known to asexually reproduce by fragmentation and/or through polyps bail-out as a response to stressful conditions (Miller and Grange 1997; Parker et al. 1997; Miller 1998; Bo et al. 2009b; Coppari et al. 2019, 2020b). Nevertheless, the successful attachment of black coral fragments or bail-outs has never been documented in the wild (Coppari et al. 2019, 2020b). For the genus Antipathella in particular, polyp bail-out has been observed under stressful laboratory conditions in A. fiordensis (Parker et al. 1997) and in A. subpinnata (Coppari et al. 2020b). For the latter species, an extensive fragmentation event was also recorded in rearing conditions (Coppari et al. 2019). Fragmentation in particular was suggested to be a successful reproductive strategy in A. subpinnata since newly formed coral fragments were characterized by (1) high survival rates (100%), (2) long viability (observations lasted 7 months with the fragments maintaining a healthy phenotype throughout this entire period, showing no signs of growth impairment), and (3) capability to successfully reattach to different substrates (Coppari et al. 2019). This suggests that asexually produced fragments, as described in Coppari et al. (2019), would be capable of dispersing across hundreds of kilometers, which could in part explain the high connectivity observed between coastal populations of A. subpinnata in our study.

Furthermore, linkage disequilibrium (LD) analysis suggests that the studied populations of A. subpinnata could be partially clonal, which potentially provides additional support that asexual reproduction indeed occurs under natural conditions in this species. LD was detected in each of the sampled localities, which is expected for taxa that asexually reproduce (Brown et al. 1980; Agapow and Burt 2001; Goss et al. 2014). Conversely, loci are not expected to be in LD in populations where sexual reproduction prevails, since genetic recombination would yield new genotypes during meiosis/gametogenesis (Brown et al. 1980; Agapow and Burt 2001; Goss et al. 2014). Clonality is further supported by (1) the observed heterozygote deficiency (HO < HE), (2) positive inbreeding coefficient (FIS) values (Table 1) and the (3) MLG analysis (Fig. S2). Interestingly, genetic clones were already detected within the genus Antipathella (Miller 1998). In Miller (1998), clonality was recorded in A. fiordensis using 7 electrophoretic allozymes, with clones occurring very commonly and across long (> 60 km) geographic distances.

Asexual reproduction in corals is known to prevail in non-optimal conditions (Waller and Tyler 2005). For example, Waller and Tyler (2005) observed that, under trawling stress, none of the Lophelia pertusa specimens were sexually mature, while they were all in a reproductive state in a preserved area. Therefore, as observed in other anthozoans, the occurrence of asexual reproduction in black corals (Parker et al. 1997) may play an important role in the survivorship of black corals in non-optimal environments, as colonization of more suitable habitats may occur through dispersal of asexual coral propagules, as observed in experimental conditions for A. subpinnata (Coppari et al. 2019).

In regard to sexual reproduction, corals can have a brooding or broadcast spawning sexual reproductive strategy, which largely determines their dispersal capacity. Since there is currently no evidence of internal fertilization in antipatharians (Wagner et al. 2011), A. subpinnata is most likely a broadcast spawning coral (Gaino and Scoccia 2010), as are other antipatharians (Miller 1998; Rakka et al. 2017). The disappearance of mature gametes observed in the summer period (mainly late August) in tagged colonies of A. subpinnata support this hypothesis, together with the absence of larvae in the female polyps (Bo et al. unpublished data, 2016–2017). Broadcast spawning in A. subpinnata could thereby be another plausible way to explain the high genetic connectivity observed between upper mesophotic (coastal) populations, although specific studies on the larval ecology of this species are still lacking.

What is currently known about black coral larvae is based primarily on studies of shallow-water antipatharians (Brugler et al. 2013). For instance, sexually produced larvae in the congeneric species A. fiordensis are lecithotrophic, negatively buoyant and negatively phototaxic, with low motility and survival time (approximately 10 days) in the water column (Miller 1998). Similar larval characteristics were observed in L. pertusa, a deep-sea scleractinian coral whose larvae also exhibit negative phototaxic behavior two weeks after their formation (Larsson et al. 2014). These life-history characteristics could reduce pelagic larval distance, as indicated by the philopatric behavior and limited recruitment distance exhibited by A. fiordensis (Miller 1998). Interestingly, the negative buoyancy observed in A. fiordensis and L. pertusa larvae does not coincide with larval characteristics of many scleractinian coral species, as upward swimming behavior is known to prevail in larvae of tropical shallow-reef corals (Bongaerts et al. 2017). This positive buoyancy has already been hypothesized to greatly increase the dispersal potential of coral propagules, particularly from deeper to shallow gradients (Bongaerts et al. 2017). Larval characteristics as described in A. fiordensis potentially indicate that the opposite scenario could also occur in some cases and that mesophotic larvae may be effectively trapped at depth. This could explain why offshore populations of A. subpinnata were more genetically distinct than coastal populations. The strong genetic subdivision observed between the two lower mesophotic (off-shore) localities could provide additional support for this hypothesis, considering that the highest FST value was detected for this grouping (Fig. 3c—SLU vs. POS: Nei’s FST = 0.22).

Gene flow between coastal A. subpinnata populations can additionally be explained when considering physical and environmental parameters, such as strong superficial currents in the Ligurian Sea (Pinardi et al. 2015), and the presence of numerous A. subpinnata forests along the western Italian coastline at similar depths (> 70 forests, Bavestrello, pers. comm.), which may be acting as stepping stones for larval dispersal. Our hypothesis that A. subpinnata, probably similar to other black corals, opts for a ‘stepping-stone’ colonization pathway is further supported when considering the fact that rocky environments become less frequent with increasing distance from the coast and separated by important bathymetric gaps (as in the case of seamounts). The inability to use the ‘stepping-stone’ strategy for larval dispersal in offshore habitats could additionally explain the observed genetic isolation between lower mesophotic populations of A. subpinnata.

Despite a clear upper mesophotic (offshore)/lower mesophotic (coastal) genetic structuring, some coastal specimens in our study seemed to originate from offshore localities, particularly from the Posada canyon (Fig. 3b). Submarine canyons are known to promote upwelling of cold, nutrient-rich waters to the sea surface (Allen and Durrieu de Madron 2009), which can facilitate transport of coral propagules to shallower areas (Canepa et al. 2014; Bongaerts et al. 2017). For instance, upwelling from deep-sea marine canyons carries huge blooms of the scyphozoan jellyfish Pelagia noctiluca to the coasts of the Tyrrhenian Sea (Canepa et al. 2014), which are then distributed by horizontal surface currents (Boero et al. 2016). This means that despite their isolation, deep refuges may aid in shallow reef recovery through the provision of coral propagules (Bongaerts et al. 2017). Although this ‘reseeding’ hypothesis remains largely untested (Bongaerts et al. 2017; but see also Bongaerts et al. 2019), it might explain why some individuals from Portofino and Savona localities showed high membership-assignment scores to the Posada canyon subpopulation (even if not directly located above the canyon).

While strong connectivity patterns were observed among the studied populations, it is likely that finer-scale genetic structuring would be revealed by incorporating a greater breadth of sampling localities and additional replicates per population. For instance, future studies could incorporate a field-sampling design where individuals from intermediate populations would be collected, as this would allow examination of a potential stepping-stone dynamics in population structuring. However, sampling of deep-sea coral individuals for genetic analysis involves the deployment of submersibles (such as ROVs) and is thereby expensive and logistically difficult. This explains the small number of representative samples obtained in our study for the ROV sampling localities.

Genetic connectivity recorded between coastal populations indicates that gene flow occurs at higher rates in coastal areas (upper mesophotic), suggesting that these coral gardens may be less susceptible to human impact (i.e., demersal fishing activities). Conversely, the observed genetic isolation in A. subpinnata lower mesophotic (offshore) forests, particularly the one from the Santa Lucia Seamount, suggests they are potentially less resilient to human impact due to a limited influx of coral propagules from adjacent habitats. The observed offshore / coastal population structuring additionally suggests that lower mesophotic A. subpinnata forests seem unlikely to act as a reservoir of genetic diversity to upper mesophotic populations (as observed in other Mediterranean coral species, Costantini et al. 2011, 2016), all of which highlights the vulnerability of Mediterranean A. subpinnata forests, particularly those located at greater depths. It is furthermore interesting to emphasize that a shift in male-to-female ratio was observed during the field sampling operations between upper mesophotic and lower mesophotic A. subpinnata forests, with male individuals becoming more prevalent with increasing depth (Bo, pers. Comm). While these preliminary field observations could have been biased by a lower sampling effort in the offshore localities and are thereby yet to be validated, the male-to-female sex ratio is an important aspect to consider in future studies as an unbalanced sex ratio has been documented to critically decrease the effective size of populations, thereby lowering population resilience to non-optimal environmental conditions (Dubreuil et al. 2010; Rosche et al. 2018). Since remotely operated vehicle footage (Bo et al. 2014) and accidental bycatch data (Mytilineou et al. 2014) have already provided evidence of large impact by fishing activities on Mediterranean black coral forests along the continental platform, but also within canyons and over seamount summits, this paper aims to raise concern and highlight the importance of enforcing management measures to achieve Good Environmental Status of these valuable A. subpinnata underwater forests under the European Marine Strategy Framework Directive.

Internationally, black coral gardens are listed as (i) threatened and declining by the OSPAR Commission (OSPAR Recommendation, 2010/8), are (ii) identified as special ecological features that require protection under the Convention of Biological Diversity (UNEP 2007), are (iii) defined as Vulnerable Marine Ecosystems (VMEs) potentially impacted by fishing activities (FAO 2009), and (iv) the communities in which they thrive fall under the 1170 Reef Habitat type of the European directive 93/43/CEE (Aguiliar and Marin 2013). All black corals are listed in Appendix II of the Convention on International Trade in Endangered Species of Wild Fauna and Flora (CITES), and recently Mediterranean species have been included in the International Union for Conservation of Nature (IUCN) Red List as Endangered or Near Threatened, mainly due to extreme longevity traits, decline in the populations and fishing impact (Otero et al. 2017). Finally, recognizing their fragility and importance, Mediterranean black corals species have been moved in the Annex II of the Barcelona Convention.

Despite this, there are few management actions established in the Mediterranean basin for the protection of cold-water coral ecosystems and even fewer initiatives have been enacted for the upper mesophotic (coastal) assemblages, due to a greater number of socioeconomic interests over the continental platform. On a national level, one Italian conservation initiative devoted to the protection of an A. subpinnata forest was recently promoted in the Southern Tyrrhenian Sea (Calabria regional deliberation #479, 06/11/2012 for the SCI ‘‘Scogli d’Isca’’) (Bo et al. 2014), while the cold-water coral site of Santa Maria di Leuca, also hosting black corals, is a well-known Fishing Restriction Area (FRA) (Rec. GFCM/30/2006/3). The same protection tool was used for another unique area, the canyon system of the Gulf of Lion (France), hosting, among other, lushing black coral forests (Rec. GFCM/33/2009/1; Fabri et al. 2014). Finally, as far as it is known, no coastal marine-protected area (MPA) includes forests of black corals, and this remains a pending managing action to pursue in the near future.