Chasing genetic structure in coralligenous reef invertebrates: patterns, criticalities and conservation issues.

Conservation of coastal habitats is a global issue, yet biogenic reefs in temperate regions have received very little attention. They have a broad geographic distribution and are a key habitat in marine ecosystems impacted by human activities. In the Mediterranean Sea coralligenous reefs are biodiversity hot spots and are classified as sensitive habitats deserving conservation. Genetic diversity and structure influence demographic, ecological and evolutionary processes in populations and play a crucial role in conservation strategies. Nevertheless, a comprehensive view of population genetic structure of coralligenous species is lacking. Here, we reviewed the literature on the genetic structure of sessile and sedentary invertebrates of the Mediterranean coralligenous reefs. Linear regression models and meta-analytic approaches are used to assess the contributions of genetic markers, phylum, pelagic larval duration (PLD) and geographical distance to the population genetic structure. Our quantitative approach highlight that 1) most species show a significant genetic structure, 2) structuring differs between phyla, and 3) PLD does not appear to be a major driver of the structuring. We discuss the implication of these finding for the management and conservation, suggesting research areas that deserve attention, and providing recommendations for broad assessment and monitoring of genetic diversity in biogenic reefs species.

Introduction

Biogenic temperate reefs include coral forests, maerl beds, coralligenous reefs and share several features with tropical reefs[1,2]. They support biodiversity by providing habitats, feeding grounds, recruitment and nursery sites for a myriad of species, and they host some of the most productive and diverse assemblages. However, they are also among the most threatened habitats globally[3-6] owing primarily to human disturbances. They are exposed to long-term degradation, which affects species composition and abundance, and can cause dramatic declines in diversity and functioning[3,7,8]. Among temperate biogenic reefs, coralligenous formations, occurring in lower subtidal or twilight Mediterranean benthic environments (between 20 and 200 m), are considered one of the most valuable and diverse habitats[9-12] (Fig. 1). In addition to being threatened by a variety of human activities[5,12,13], episodes of mass mortality of benthic invertebrates over hundreds of kilometres on the Northwestern Mediterranean shores have been recorded[14]. Several international bodies (e.g., UNEP RAC-SPA, IUCN, EU), as well as local authorities, set the conservation of coralligenous reefs as a priority. The EU Habitats Directive, the Bern Convention, and the Barcelona Convention Action Plan have included coralligenous habitats in conservation plans, and the scientific literature has stressed the need to implement effective conservation policies[3,5,12,15-17]. To quantify the impacts of a variety of environmental and anthropogenic stressors, several studies applied a meta-analytical approach to data on the spatial and temporal distribution of coralligenous species. Giakoumi et al.[15] have also stressed the importance of connectivity (e.g., dispersal of individuals among discrete populations) in planning the conservation of coralligenous reefs, as already suggested for other marine habitats[18-20].

Figure 1

Mediterranean coralligenous outcrops (Photo by Simonepietro Canese).

Connectivity is essential to population resilience following disturbance and underlies the dynamics of populations and metapopulations[18,21]. Its importance has been recently recognised at the policy level, as highlighted, for example, by the European Union “Marine Strategy Framework Directive” and by the Strategic Plan for Biodiversity 2011–2020 by The Convention on Biological Diversity (as a tool for the preservation of biodiversity and the development of marine spatial planning).

Different methods can be used to quantify connectivity (e.g., direct observation using tagging, biophysical and ecological modelling and population genetic analyses). Genetic methods are frequently used to assess population connectivity when it is difficult to measure dispersal directly, as in the case of sessile species in the subtidal marine realm[21,22]. Multispecies, multiscale and multimarker studies on genetic connectivity have been proven to be effective in supporting biodiversity conservation planning[18,23,24]. However, most of these studies focused on tropical habitats and included comparative phylogeography, descriptive[25,26] and meta-analytical reviews[27,28]. Only a few studies have investigated patterns of genetic connectivity in temperate biogenic reefs[20,25,29-31] and in Mediterranean coralligenous reefs[31,32].

The most commonly used metric to estimate connectivity based on genetic data is Wright’s fixation index (F), summarising genetic variation among geographically separated populations. The use of F to assess differentiation has been criticized[33] since it is affected by within- population heterozygosity and tends to underestimate differentiation between populations as heterozygosity increases. Moreover, the assumption that pattern of spatial genetic differentiation (F) primarily reflect variation in gene flow is questioned, since spatial genetic structure is also shaped by the demographic history of the populations[34]. Nevertheless, up to now is still the most suitable used metric[28,35-37]. Based on this metric, several theoretical patterns of genetic structure have been defined (the island model[38], isolation by distance[39], the metapopulation model[40] and the chaotic genetic pattern[41]). Empirical evidence, summarised in a recent review[42], indicates that these theoretical patterns could be observed in the same species, depending on the spatial and/or temporal scale of analysis. Understanding the drivers and the processes underlying these patterns of genetic structure is challenging. Abiotic (e.g., hydrodynamism and substrate geomorphology) and biotic factors, related to the life history traits (e.g., reproductive strategies and larval ecology), drive patterns of genetic structuring. However, specific information on causal relationships between habitat features and the observed structure is still lacking. Only recently has the seascape ecology approach been applied in genetic studies[43]. Furthermore, for most marine invertebrates, information on larval ecology (type of development, duration, and behaviour) is limited, and most information refers to the pelagic larval duration (PLD). PLD is commonly considered the best proxy for effective larval dispersal; therefore, its influence on connectivity has been extensively studied[28,35-37,44]. However, contrasting conclusions on the relationship between PLD and genetic structure have been reached[28,31,36,37]. A weak correlation between PLD and genetic differentiation has been highlighted in recent reviews on fishes and invertebrates[28,35], suggesting that other life history traits[45] and/or differences in effective population sizes[36] could affect estimates of genetic differentiation. Nevertheless, these studies agree that direct-developers and benthic species with low mobility showed a relatively lower genetic structure to mobile species with planktonic larvae[29,46].

The role of the biological and ecological features of a species in determining the genetic structure may be tangled by the difference in the resolutions of the molecular markers used[36,47]. Meta-analyses on the effect of dispersal traits on population differentiation in a wide range of species (from angiosperms to crustaceans) in temperate habitats (from benthic to pelagic environments) have recently been published[29,31]. These studies show that genetic structuring was stronger in species with duration of dispersal phases shorter then 15 days compared to species with higher motility and dispersal, regardless of the presence of oceanographic discontinuities[31].

In the present study, we have summarised the literature on the genetic structure of sessile and sedentary benthic invertebrates with negligible adult mobility that are commonly found in Mediterranean coralligenous habitats. Our aims were 1) to identify features and/or patterns in genetic diversity and structure and 2) to understand processes generating spatial genetic structuring among populations. The quantitative review has been carried out using two different statistical approaches (linear regression modelling and meta-analysis) to investigate the roles of genetic markers, phylum, PLD and geographical distance in determining the observed patterns of genetic structure. To eliminate the confounding effect of the presence of geographical discontinuities (i.e., recognized biogeographic barriers[10,31]), we screened all the literature and analysed separately data referring to populations sampled within single biogeographic areas. Finally, we discussed the implications of these findings for management and conservation of Mediterranean coralligenous reefs.

Results

The full data set comprised 78 records from 50 papers published between 2000 and 2016 (Supplementary Table S1). The number of papers using microsatellite loci ranged from 0 papers in 2001 to 6 papers in 2016, while studies using mitochondrial DNA markers were evenly distributed in time (Fig. 2). The only study with nuclear markers developed by transcriptome sequences[48] was published in 2016. The investigated species occur across the whole Mediterranean Sea, although most papers were focused on the Northwestern Mediterranean coast (from Gibraltar to the Ionian Sea). Only a few studies included samples collected in the Southern and Eastern Mediterranean Sea (e.g., Tunisian, Greek, Turkish and Cyprus coasts; Fig. 3).

Figure 2

Type of molecular markers used across time (beginning 2000-end 2016) based on the 50 papers considered in the manuscript. Some authors analysed more molecular markers in the same publication. Nuclear: nuclear sequences; micro: microsatellite loci; allozyme; mtDNA: mitochondrial DNA; trans: transcriptome-based nuclear markers.

Figure 3

Sampling effort per phylum in the Mediterranean Sea. All study sites sampled in the 50 Mediterranean studies considered in this review have been grouped into 20 clusters (areas) based on their Euclidian distance. In each panel, circles are centred on the centroid of each area and sized proportionally to the absolute maximum number of species sampled (circle sizes are shown in the top right corner), while the colour gradient represents the number of study sites sampled (reference scale shown in the bottom right corner). The maximum number of species sampled per phylum is reported in parentheses in the panel title. Numbers on the axes are latitude and longitude. The map has been created with R software[76] (R Development Core Team; www.r-project.org) using the package ‘mapdata’ v2.2–6.

The 50 analysed papers refer to 24 species and three species complexes (Ophioderma longicauda, Amphipholis squamata and Ophiothrix fragilis) associated with coralligenous assemblages (Supplementary Table S1, but see[11]). The dataset includes four major taxonomic groups: Porifera, Cnidaria, Echinodermata and Tunicata (Supplementary methods). Cnidaria and Porifera are the most represented taxa, with eight species each, followed by Tunicata, with six species. Seventeen of the 27 species (24 species + 3 species complex) were analysed at different spatial scales and using different molecular markers.

Most of the species are brooding (79.3%), with one (Paramuricea clavata) that is an external brooder and six that are broadcaster (Supplementary Table S1; Supplementary Table S2). The studied taxa included a wide range of pelagic larval durations: brooding species that lack a dispersive larval stage (e.g., the brittle star Amphipholis squamata) and species with long pelagic larval durations (e.g., the sea star Marthasterias glacialis; Supplementary Table S1, Supplementary Table S2). Regarding larval feeding types, 72.4% of the species are lecithotrophic, three are planktotrophic and two have direct development. The Ophioderma longicauda species complex showed both lecithotrophic and direct development, depending on the lineages/genetic clusters considered (Supplementary Table S1). For three cnidarians, the larval feeding type was not known.

Genetic patterns within the Mediterranean Sea

Three patterns of genetic structure (i.e., no genetic structure, chaotic, and IBD; Supplementary Table S1) have been found. Two cnidarians (Leptopsammia pruvoti and Eunicella singularis) and three echinoderms (the brittle stars Ophioderma longicauda C3, Ophiothrix sp. II and the starfish Marthasterias glacialis) showed no genetic structure at the analysed spatial scales (but see Discussion). Ten species (37.03%) showed a chaotic genetic pattern, and six species (22.2%) showed an isolation-by-distance pattern (Supplementary Table S1). Eight species (29.62%) showed a mixed pattern of genetic structure (IBD/chaotic), depending on the molecular markers and on the spatial scale analysed. Moreover, within the three species complexes (Ophioderma longicauda, Ophiothrix fragilis, and Amphipholis squamata), contrasting patterns of genetic structure were observed depending on the lineage[48,49]. The most recent paper[50], which reanalyses previous data, showed a clear chaotic genetic pattern within the Amphipholis squamata species complex.

Most studies were not designed to test genetic structuring across Mediterranean oceanographic discontinuities, but populations were haphazardly sampled across the species distribution range. Nevertheless, the majority of taxa showed strong genetic breaks across their geographic range of distribution, with some common geographical patterns of genetic differentiation. The role of oceanographic fronts in setting connectivity patterns was recently investigated[31]. However, we found new areas of discontinuities that affect several species. For example, a barrier to gene flow in the Ligurian-Tyrrhenian transition[51,52] was observed in two species, but not in others (e.g., Chondrosia reniformis[20]). Moreover, most of the observed species show a significant genetic structure across the Western and Eastern Mediterranean transition, even if the location of the shift varies, depending on the species.

Estimated dispersal distance at regional/local scales

We tested for the presence of a significant isolation-by-distance pattern by comparing pairwise genetic and geographic distances in populations sampled within a defined biogeographic areas (i.e. Mediterranean sub-basins; see Materials and column U in the Supplementary Table S1). We observed significant isolation-by-distance patterns in five species[27,35,53], for which we estimated the dispersal distance. For C. rubrum and C. crambe we included in the analysis the estimated dispersal distance calculated by Ledoux et al.[54] and Calderon et al.[55], respectively. They found an estimated dispersal distance ranging between 20 and 30 cm in C. rubrum and of approximately 35 cm in C. crambe. For all the seven species considered, we estimated a dispersal distance ranging from 0.2 to 300 m (Fig. 4, Supplementary Table S3). No correlation between the estimated dispersal and the taxonomic groups was found.

Figure 4

Estimated dispersal of the seven species showing an isolation-by-distance pattern (IBD < 0.05) using the method of Kinlan and Gaines (2003). For Corallium rubrum and Crambe crambe we used the already published values of mean dispersal distance calculated by Ledoux et al.[55] and by Calderon et al.[54] using an “individual based” sampling.

Quantitative assessment of the genetic structure

The meta-analysis revealed that the regression slope between pairwise F and geographic distance was positive and significant (slope value ± standard error = 0.0076 ± 0.0023, F(d.f 1,30) = 10.56, P = 0.0028; Fig. 5) when all the species were included in the analyses (i.e., not clustered based on zero, short or long PLD). We detected significant differences in F /geographic distance relationships among phyla, with Cnidaria appearing more structured than Echinodermata and Porifera (Fig. 5), while we did not detect any pattern related to the PLD or to the genetic marker (Table 1). We found similar results when we analysed only the species with a long PLD (i.e., ≥1 day), although differences among phyla were less evident (Table 1).

Figure 5

Forest plot summarising the results of the meta-analysis. Grey dots represent the effect size (ES) of single observations, their size represents the weight assigned to the observation while the horizontal lines represent the 95% confidence interval (CI.) of the ES. Red dots and lines represent the ES and the 95% CI of the levels of a factor (e.g. factor marker with levels Microsatellites and mtDNA). Grey and red dots refer to the black and red axes respectively. On the left, each observation is uniquely identified by its row number in the META-set followed by the literature reference. On the right, the ES [ES ± 95% CI] and the species of each observation are reported. A confidence interval including the zero (dashed line) indicates a non-significant ES.

Table 1

Meta-analytical models for the population structure at the Mediterranean scale. AICc value are presented for models fitted using the maximum likelihood estimator on datasets with (All PLD) and without species with no or short PLD (Long PLD). Significance of single factors was tested as a log likelihood ratio (LRT) against the null model. Degree of freedom are reported as “d.f”. For each model it is reported which phyla were included in the analysis and the sample size (N).

All PLD (no, short and long PLD)	AICca	ΔAICcb	d.f.	LRT	p-value	Phylac	N
Null model	−137.3853	0.014
Marker	−137.3996	0	1	2.475	0.116	C,E,P,T	31
Maximum PLD	−136.4934	0.906	1	1.568	0.210	C,E,P,T	31
Marker + Maximum PLD	−136.0300	1.370
Phylum			2	7.121	0.028	C,E,P	30
Long PLD only
Null model	−113.6858	0
Marker	−112.6997	0.986	1	1.759	0.185	C,E,P	21
Maximum PLD	−112.1282	1.557	1	1.188	0.276	C,E,P	21
Marker + Maximum PLD	−111.3647	2.321
Phylum			2	5.399	0.067	C,E,P	21

aLowest AICc values are in bold.

bDifference in AIC(c) between a model and the best model.

cC = Cnidaria, E = Echinodermata, P = Porifera, T = Tunicata.

Using the ANCOVA approach, we investigated the global F as a function of marker, maximum PLD and maximum geographical distance. In the analysis of the dataset including all species (with zero, short or long PLD), the best model (Marker × Max. PLD + Max. Distance, Table 2) showed a significant interaction between “Marker” and “Maximum PLD” (LRT = 9.557, d.f. = 1, P = 0.002), while the geographical distance was not significant (LRT = 3.64, d.f. = 1, P = 0.056). As a rule of thumb, models having a ΔAIC < 2 are considered to be alternative models with evidence of support[56]. The only other model with ΔAIC < 2 included the “Marker × Max. PLD” interaction (Table 2). The top model indicates a small but positive effect of the geographical distance on F values. PLD has a negative effect on F values obtained using mtDNA markers, while it has a positive effect when using microsatellites (Fig. 6). Restricting the analysis only to the species with a long PLD, the best model is the one with the additive effect between “Marker” (LRT = 5.446, d.f. = 1, P = 0.020) and “Maximum Distance” (LRT = 3.341, d.f. = 1, P = 0.068) (Table 2). The second and third models with ΔAICc < 2 included the factor “Marker”. The third model retains the “Marker × Max. PLD” interaction, consistently with the analysis performed on all the observations.

Table 2

ANCOVA model for the population structure at the Mediterranean scale. AIC value for analyses with (All PLD) and without (Long PLD) species with no or short PLD. Sample size is reported as N within brackets.

	All PLD (N = 45) (no, short and long PLD)		Long PLD only (N = 34)
	AICca	ΔAICcb	AICca	ΔAICcb
Full model
Marker × Max. PLD + Max. Distance	122.08	0	Not testedc
Interaction only
Marker × Max. PLD	122.90	0.82	91.18	1.09
Additive effect
Marker + Max. PLD + Max. Distance	128.82	6.74	92.38	2.28
Marker + Max. PLD	129.18	7.10	93.38	3.28
Marker + Max. Distance	126.36	4.28	90.10	0
Single factor
Marker	127.98	5.90	90.67	0.58
Maximum PLD	130.07	7.99	95.62	5.52
Maximum Distance	127.94	5.86	92.78	2.68

aLowest AICc values are in bold.

bDifference in AIC(c) between a model and the best model.

cDue to reduced sample size.

Figure 6

Global FST as a function of Maximum PLD and Maximum distance. The two surfaces depict the predictions of the regression model Log(FST) ~ Marker × Maximum PLD + Maximum distance. The two predictor variables are the x and y axes at the base of the cubes while the z-axis represents the response variable in the predictor scale (i.e. exponentiated). The red lines highlight how Global FST varies with Maximum PLD while the Maximum distance is held constant (here the red dashed lines correspond to the mean value of distance of studies for each marker). The model equation is Log(FST) = −3.0686 + 0.0080 × Maximum PLD + 0.0007 × Maximum distance for Microsatellites; Log(FST) = −2.9392–0.0897 × Maximum PLD + 0.0007 × Maximum distance for mtDNA.

Discussion

Our study showed that among the roughly 300 species of benthic sessile and low mobility invertebrates inhabiting the Mediterranean coralligenous habitats[11], genetic diversity and structure have been investigated for less than 30 species. Most studies address sponges, cnidarians and echinoderms, while there are no studies on bryozoans and polychaetes, despite their importance as bioconstructors[11].

The quantitative approach revealed that 1) most species show a significant genetic structure within the biogeographic areas considered, 2) the structuring differs between phyla, 3) PLD does not appear to be a major driver of the structuring, and 4) the genetic marker used did not introduce a bias in the analyses.

Technical bias, such as uninformative molecular markers and/or low sampling efforts (in terms of populations and individuals), may affect the ability to detect genetic differentiation. This could be the case for Leptopsammia pruvoti[57] and Marthasterias glacialis[58], where the sample size is small. Regarding other species, the absence of genetic structure could be attributed to the broadcast reproductive mode (Ophioderma longicauda LI (C3)[48]), to the large effective population size (Ophiothrix sp. II[59]) or to the small spatial scales investigated (Eunicella singularis[60]).

Our qualitative synthesis highlights that, at the larger scale, a constraint to dispersal acting across multiple species was located in the Eastern-Western Mediterranean transitions. This area of discontinuity is well known for several marine species[30,31]. Genetic shifts have been found in different locations according to the species: across the Siculo-Tunisian strait[61], in the mid-Ionian Sea[20], or in the eastern Ionian[49,62]. Other critical areas of discontinuities were observed[31] and were species specific or shared by few species (e.g., the Balearic front[63,64] and Almeria Oran front[65,66]).

Sessile and sedentary species inhabiting the Mediterranean coralligenous are highly structured, with an estimated dispersal distance within the range of hundreds of metres. These results are in accordance with the patchy distribution of the coralligenous outcrops in the Mediterranean Sea. Coralligenous outcrops are widely spread across the Mediterranean Sea but with a high spatial heterogeneity in distribution and species composition at small scales[16,67]. Both these features contribute to explain the patchy distribution of the organisms, and can influence processes acting in the early life stage of the organisms (recruitment pattern, pre and post settlement processes). Moreover, stochastic processes have a major influence on sessile and sedentary invertebrates with low dispersal capability (see below). Genetic differentiation was observed at smaller distances in Cnidaria, where the majority of populations were separated by less than 152 km, than in Echinodermata and Porifera, where the majority of populations were separated by more than 382 km with a median distance of 175 km (Supplementary Figure S1).

The Echinodermata analysed in this study present a wide range of pelagic larval duration (see below), different reproductive modes and different nutritional modes, while Cnidaria do not. Previous studies found differences between major taxonomic groups (e.g., vertebrates vs. invertebrates[44,46,68]). Chust et al.[29] and Pascual et al.[31] showed that species with similar life history traits (e.g., benthic sessile or sedentary vs. benthic vagile or pelagic) present a similar genetic pattern, but they did not test this effect across Phyla.

We found only a weak correlation between PLD and genetic structure. Therefore, despite PLD being the better-known larval feature, caution must be used when employing PLD as a proxy of effective dispersal. Several empirical studies[20,37], as well as recent meta-analyses[27,35], have found minimal-to-modest support for the correlation between PLD and genetic structure. Constraints in the estimation of PLD could explain its limited resolution power[28]. In fact, in marine invertebrates, PLD is estimated based on larval features or by laboratory experiments, which may not be representative of field conditions. Our qualitative and quantitative synthesis shows that the ability to predict patterns of connectivity based on PLD is limited and challenged. Other biological traits (e.g., larval swimming ability and reproduction timing[69]), demographic parameters (e.g., fecundity and population size[45]), and recruitment success[70] may shape population genetic structure. These biological traits, together with the oceanographic condition, influence the reproductive success, leading to an increasing temporal and spatial variance of the effective population size (Ne). Moreover, since F is affected by both the migration rate (m) and effective population size (Ne), a variation in Ne could also shape the correlation between F and PLD[36]. Taking in account variability of the effective population size is a key task since it can mitigate the effect of genetic drift (see below).

ANCOVA analysis on mitochondrial DNA markers showed a significant negative correlation between PLD and population structure, suggesting that PLD might be a better predictor of genetic connectivity patterns when markers with low level of polymorphism and mutation rate were used. In contrast, the non-significant correlation observed using microsatellites, which have higher mutation rate, suggested that PLD is not the only factor that drives the genetic structure at a time scale of few generations[71].

In our study, many species (10 of 27) present a chaotic genetic patchiness (sensu Johnson & Black[72]). In the literature, this pattern is mainly related to selection, sweepstakes reproductive success[73], collective dispersal[74], and temporal shifts in local population dynamics[75]. However, understanding the drivers determining this genetic pattern is challenging[75]. For example, local hydrodynamic conditions affecting both larval dispersal and the stochastic delivery of larvae from different sources could determine a high variability in reproductive success or mortality related to pre- and post-settlement processes[72]. The stochasticity of the early life history dynamics[76] is related to the intra- and inter-specific interaction, the reproductive mode and time, the timing of larval release and larval behaviour. Laboratory experiments on life cycles and field manipulative experiments have to be implemented to better understand the primary processes acting at different temporal and spatial scales. These studies are extremely difficult to carry out due to the difficulties in manipulating small larvae and individuals in the field. However, this review, as well as other recent papers[74-76], emphasize the need to evaluate the effective population size and estimate genetic drift in the field. For example, studies that compare the genetic makeup of dispersing larvae with the post-dispersal stages can help to clarify the role of genetic drift (e.g., Riquet et al.[76] for Crepidula fornicata and Costantini et al. submitted for Corallium rubrum).

Our synthesis shows that the sessile and low-mobility species inhabiting the coralligenous habitats have high genetic structuring. Occurrence of structured breeding units and genetically differentiated populations at the scale of kilometres or less suggest that conservation strategies have to consider smaller scales than previously thought[77,78]. A genetically differentiated population is indeed more vulnerable to local extinction when isolated from other populations[22], and isolation seriously undermines its resilience potential. Recovery of species threatened mainly by human activities (e.g., fishing, tourism, marine urban sprawl) that act at a local scale, may be supported by nearby healthy populations or by local protection programmes. However, these actions are helpless for impacts due to climatic changes (e.g., increasing seawater temperature, ocean acidification) that act at global scales and affecting all populations.

Major efforts should be dedicated to understanding local-scale ecological processes, to developing seascape genetic approaches at the local scale, to enhancing restoration plans and to implementing the creation of Sites of Community Importance (e.g., Kelly & Palumbi[37], Micheli et al.[79], van der Heyden et al.[26]). Characterisation of the genetic structure of species inhabiting fragile Mediterranean coralligenous reefs is also needed to fill the gaps in genetic knowledge on some of the most important bio-constructors (e.g., Bryozoa, Polychaeta[11]).

We believe that comparative phylogeographic and population genetic studies, along with studies on demography and reproductive features, are fundamental to understanding the functioning of coralligenous habitats. To successfully implement these approaches, the scientific community should prioritise sampling design protocols[45]. Moreover, international initiatives promoting collaboration could substantially contribute to this aim (e.g., Regional Activity Centre for Specially Protected Areas UNEP—RAC/SPA—meeting on coralligenous and bio-concretions, Union for the Mediterranean initiatives) by establishing a network of scientists and extending the studies to the eastern and southern Mediterranean Sea.

Promoting public awareness and implementing effective communication based on coralligenous flag species with high aesthetic and naturalistic values could help stimulate interest for the conservation of these complex, yet unique, Mediterranean habitats.

Methods

Literature search and data extraction

We searched the ISI Web of Knowledge database (from 2000 to 2016, cutoff date 31 December 2016) for articles addressing the genetic variation and population structure of marine benthic invertebrates with low (negligible) adult mobility living mainly (but not exclusively) in the Mediterranean coralligenous reefs. The combination of keywords used in our queries is provided in the Supplementary methods. We examined the title and abstract of each resulting article and retained only those articles deemed eligible for the scopes of our research (exclusion criteria: primer and technical notes, no population genetic studies and papers dealing with deep populations (below a 50 m depth); Fig. 7). For each investigated species, we extracted qualitative and quantitative information regarding life history traits, molecular markers, and fixation indexes. Then, we merged the data in a unique dataset (hereafter “the full dataset”; see Supplementary methods for dataset and method description). Each dataset row represents an “observation” uniquely defined by the combination of study, species, molecular marker, sub-basin in which data were extrapolated, and, when indicated, genetic lineage (but see Supplementary Table S1). Whenever possible, we extracted values of both global and pairwise F as response variables to assess if population genetic structure is affected by genetic markers, pelagic larval duration, phylum and geographical distance. F is a measure of population differentiation ranging between one, indicating fixation of different alleles in each population (i.e., no gene flow), and zero, indicating the absence of differentiation (i.e., high gene flow[80]). Despite controversy over the use of Wright’s fixation index (F[81]) as a measure of relative differentiation[33,34], F is commonly used in population genetics and is therefore the most viable parameter to include in a multi-study analysis (e.g.,[28,35,44]). We observed a strong relationship between global and mean pairwise values in the studies where these values were both present (R2 = 0.96, P = 0.0001). Therefore, in those studies where pairwise F were the only available values, we estimated the global F for the species using the following linear fit: global F = 0.9733 × (mean pairwise F) + 0.0052. For each study we extracted pairwise F values to compare populations occurring within geographic areas delimited by recognized oceanographic or biogeographic barriers (see[31]). Maximum geographical distances included in the analysis range from less than 1 km in study carried out at a very fine scale, to about 1000 km in the Northwestern Mediterranean basin (Column U in Supplementary Table S1).

Figure 7

Flow diagram summarising the process for the inclusion of studies in the descriptive analysis and in the meta-analyses.

In addition to F we calculated the maximum distance between localities from the coordinates provided in each study. We extracted data from publication plots, maps, and tables or using the software Plot Digitizer 2.6.4 when information was graphical (e.g., isolation-by-distance figures). For studies providing geographic coordinates of sampling locations, but not geographic distances, we calculated the least-cost distance between sampling points using Google Earth.

All F values derived from mitochondrial data were converted using the equation provided by Baco et al.[46] to accurately compare haploid mitochondrial values with data from diploid nuclear markers[82].

Thus, using a linear regression between pairwise F (transformed using Rousset[39]) and the natural logarithm of the geographic distance between locations in each individual observation, we calculated the IBD slope and its associated standard error (columns N and O in Supplementary Table S1). Linear regressions were performed using the R programming language[83].

To estimate the dispersal distance at the regional scale for each species, we re-calculated the isolation-by-distance (IBD) ant tested the pattern (option ‘P-value’ for the Mantel test) between the new matrices of pairwise F and geographical distance using the online version of GENEPOP 4.2 (http://ge-nepop.curtin.edu.au). In species showing a significant IBD, we estimated the dispersal distance applying the method of Kinlan and Gaines[82]. An alternative for the analyses was the Rousset[39] method, however, both methods provide similar results (a strong log-linear relationship[46]), and the Kinlan and Gaines[82] method was already used in other metanalysis[29,46]. We used a power function model (dispersal distance = 0.0016 (IBD slope)-1.0001) following Chust et al.[29]. When a species was analysed with more than one method and/or using both mitochondrial markers and microsatellites, we considered the study carried out at the smallest spatial scale and using microsatellite loci. For Corallium rubrum and Crambe crambe we used the already-published values of dispersal distance calculated using an “individual based” sampling by Ledoux et al.[54] and Calderon et al.[55], respectively.

Descriptive assessment of the genetic structure and dispersal estimates

We undertook a qualitative review of all the studies listed in the full dataset (Supplementary Table S1, including those not considered in the two quantitative approaches) to describe the common phylogeographic and population structure patterns across species. We assigned each observation in the full dataset to three categories of spatial genetic structuring based on the significance of global F and IBD, as found directly in each paper, as follows[68]:

Chaotic structure, defined as genetic differentiation of samples when global F was significant (P < 0.05) but IBD was non-significant (Mantel correlation, P > 0.05);

Isolation by distance (IBD) when the Mantel correlation was significant (P < 0.05);

No genetic structure: when both global F and the Mantel correlation were not significant (P > 0.05).

Quantitative assessment of the genetic structure within the Mediterranean Sea

We investigated the population genetic structure observed for several species of the coralligenous reef as a function of different factors by using both meta-analytical and linear regression approaches. Statistical analyses were restricted to mtDNA sequences and microsatellite observations since the information was insufficient for other marker types. We first performed meta-analyses to investigate whether the genetic structure of a population was affected by the genetic marker used, by the PLD (i.e., the maximum PLD, as it seems better correlated with the population genetic structure compared to the mean larval duration[35]), or by the phylum. In a meta-analysis, the results from different studies are aggregated and compared by identifying a suitable ‘Effect Size’ (ES, i.e., the response variable), and importantly, the relative relevance of individual observations is weighted using the error associated with the ES[84]. Our ES was defined as the IBD slope, describing the relationship between the F and the geographic distance, while we used the standard error of the slope to calculate weights. We performed meta-analyses only on a subset of the full dataset (hereafter “META-set”), for which we were able to obtain the IBD slope and its associated error from at least four pairwise comparisons. The META-set consisted of 31 observations (19 microsatellites and 12 mtDNA). Models were fitted using random-effect models with a maximum likelihood estimator (Table 1 and Supplementary methods) and compared using AICc[56] (Akaike information criterion corrected for small samples). All the analyses were then repeated only on META-set observations with a long PLD (i.e., ≥1 day), excluding observation with no or short (i.e., hours) PLD. The significance of factors was tested by log-likelihood ratio testing between the model of interest and the reduced model (i.e., without that factor, e.g., single factor vs. null model). Models were then refitted using a DerSimonian-Laird estimator[84] to produce a forest plot.

ANCOVA analyses were conducted on a subset of the full dataset that included only observations for which global F was available (“ANCO-set”, N = 45, of which 25 were microsatellites and 20 mtDNA).

We used ANCOVA to investigate the relationship between global F, “Marker”, “Maximum PLD”, and maximum geographical distance (“Maximum Distance”). To meet the ANCOVA assumptions, the response variable “global F” was log transformed, and models were fitted using GLS (Generalised Least Squares, Supplementary methods) with a constant variance function structure using “Marker” as a grouping factor. Assumptions were checked by graphically inspecting the residuals.

Analyses were repeated for long PLD observation only in ANCO-set (N = 34, 24 microsatellites and 11 mtDNA), limiting the number of variables in the models to a maximum of three Statistical analyses were conducted in the R programming language[83] (see Supplementary methods for the R packages used).

Supplementary data

Supplementary Table S1