Population structure and connectivity in the Mediterranean sponge Ircinia fasciculata are affected by mass mortalities and hybridization.

Recent episodes of mass mortalities in the Mediterranean Sea have been reported for the closely related marine sponges Ircinia fasciculata and Ircinia variabilis that live in sympatry. In this context, the assessment of the genetic diversity, bottlenecks and connectivity of these sponges has become urgent in order to evaluate the potential effects of mass mortalities on their latitudinal range. Our study aims to establish (1) the genetic structure, connectivity and signs of bottlenecks across the populations of I. fasciculata and (2) the hybridization levels between I. fasciculata and I. variabilis. To accomplish the first objective, 194 individuals of I. fasciculata from 12 locations across the Mediterranean were genotyped at 14 microsatellite loci. For the second objective, mitochondrial cytochrome c oxidase subunit I sequences of 16 individuals from both species were analyzed along with genotypes at 12 microsatellite loci of 40 individuals coexisting in 3 Mediterranean populations. We detected strong genetic structure along the Mediterranean for I. fasciculata, with high levels of inbreeding in all locations and bottleneck signs in most locations. Oceanographic barriers like the Almeria-Oran front, North-Balearic front and the Ligurian-Thyrrenian barrier seem to be impeding gene flow for I. fasciculata, adding population divergence to the pattern of isolation by distance derived from the low dispersal abilities of sponge larvae. Hybridization between both species occurred in some populations that might be increasing genetic diversity and somewhat palliating the genetic loss caused by population decimation in I. fasciculata.

Introduction

In the past decades, massive mortalities of marine invertebrates have been reported in all oceans and seas, caused by a variety of factors from contaminants to thermal-related diseases, and seriously threaten the biodiversity of the areas where they occurred (Pérez ; Laboy-Nieves ; Sala and Knowlton, 2006; Webster, 2007; Garrabou ). From a genetic point of view, massive mortality events usually produce bottlenecks (dramatic reductions of the effective population size), loss of rare alleles and thus an overall reduction in genetic diversity. Over time, the strong effect of genetic shift on reduced populations can ultimately drive populations to inbreeding depression (Frankham, 2005; Hauser and Carvalho, 2008). Eventually, if genetic diversity remains low and inbreeding increases, detrimental effects will be observed in the affected population and the overall species fitness (Hauser and Carvalho, 2008; Frankham, 2010). In addition, the reduction of genetic diversity is predicted to decrease species resilience to environmental stressors or perturbations and reduce their adaptive evolutionary potential (Spielman ; Frankham, 2005). In this sense, the assessment of genetic diversity and structure, connectivity and inbreeding of the affected species is key to estimating the effects of environmental perturbations (for example, thermal stress), disease outbreaks and mass mortalities on population demographics, and hence their potential for recovery.

Sponges are among the most simple and primitive metazoans, yet they have various biological and ecological properties that make them an influential component of marine ecosystems. Marine sponges provide refuge for many small invertebrates and are critical to benthic–pelagic coupling across a wide range of habitats (see, for example, Ilan ). In some ecosystems, sponges show an abundance and diversity of species that can comprise most of the living mass (see, for example, Diaz and Rützler, 2001; McClintock ). In addition, they have essential roles in the biogeochemical cycle of silica (Maldonado ) and dissolved organic matter (de Goeij ). Unfortunately, like many other marine invertebrates, sponges are currently threatened by numerous environmental perturbations caused by global climate change, especially in near-shore ecosystems (Hughes ; Garrabou ; Webster ; Pérez and Vacelet, 2014). The vulnerability of these ecosystems is increased by higher harvesting rates of economically important species, human recreational activities and coastal development, terrestrial runoff and the presence of heavy metals and pollutants in harbors and marinas (Commendatore ). Sessile invertebrates, like sponges, are also particularly susceptible to temperature anomalies and continuous exposure to waterborne pathogens and pollutants because of their immobility (see, for example, Cebrian ; Garrabou ).

Reports of sponge disease have increased dramatically in recent years, with sponge populations decimated at some sites in the Mediterranean and Caribbean (see Webster, 2007 for a review). In particular, Mediterranean and Adriatic dictyoceratids (especially Ircinia spp.) have been consistently reported as affected by episodes of mass mortalities since the 1990s (see, for example, Garrabou ; Maldonado ; Cebrian ; Stabili ). Three species of the genus Ircinia are widely distributed in the Mediterranean Sea, with populations occurring in shallow rocky coasts: Ircinia fasciculata, Ircinia variabilis and Ircinia oros. Ircinia species are aspiculate sponges, chemically diverse and harbor complex host-specific communities of bacterial and cyanobacterial symbionts (Erwin ,2012b; Pita ). Both I. fasciculata and I. variabilis are affected by a disease-like condition (Stabili ) reported by several authors in 2008 in the North and South Western Mediterranean Sea (Maldonado ; Cebrian ). Whereas a monitored population of I. fasciculata in the southern coast of Spain (Granada) experienced 30% of mortality (Maldonado ), two monitored populations of I. fasciculata in Cabrera and Scandola Marine Reserves were decimated with ∼90% mortality (Cebrian ). This disease-like condition, which was characterized by small round pustules that developed into large areas of necrotic tissue and eventually resulted in the death of the sponges, was linked to abnormally high seawater temperatures (>24 °C; Cebrian ). Moreover, I. fasciculata disease probably resulted from a proliferation of several bacterial pathogens (Maldonado ). Given the large impact of the disease-like condition on some populations, Cebrian predicted local extinctions of the species in the studied areas. In situations under environmental stress, extinction rates could also be increased substantially by high inbreeding levels that are predicted to decrease species resilience to environmental stressors and therefore reduce the adaptive evolutionary potential of the species (Spielman ; Frankham, 2005). Thus, determining the genetic diversity and level of connectivity in sponge populations of this species is crucial when designing management strategies to prevent its extinction. Sponge populations are usually highly structured in ‘stable conditions' (that is, healthy populations), with relatively low gene flow and high levels of inbreeding (see, for example, Duran ; Dailianis ; Chaves-Fonnegra ; Giles ; Pérez-Portela ). Therefore, even in ‘stable conditions' the recovery of sponge populations after massive reduction of the effective population size can be limited by their own biology, dispersal potential and oceanographic currents, and one could even predict lower recovery rates in stressed populations.

The two target species of the present study are I. fasciculata (Pallas, 1766) and I. variabilis (Schmidt, 1862). These species are phylogenetically closely related (see Erwin ) with an external appearance so similar and plastic that it can lead to uncertainties in identification, even though they are generally considered different species (Pronzato ). From a morphological point of view, I. variabilis contains inorganic inclusions in its skeletal fibers, whereas I. fasciculata lacks such inclusions, prompting Pronzato to transfer I. fasciculata to the genus Sarcotragus and rename the species, S. fasciculatus (Pallas, 1766). From a molecular point of view, a recent analysis using the nuclear ribosomal gene internal transcribed spacer 2 (ITS-2) indicated that I. variabilis and I. fasciculata were indeed different but closely related species with a divergence of 2.97%±0.38 (Erwin ), even though analyses of mitochondrial sequences lacked enough resolution to differentiate among them. In light of these results, we decided to maintain the name I. fasciculata in this study until a thorough analysis is performed to solve its true generic status. Each species (I. variabilis and I. fasciculata) harbors several specific microbes within their symbiotic content (Erwin ), and interestingly, even though these species live in sympatry, they present slightly different habitat preferences and large differences in growth rates (Turon ). In marine sperm broadcasters, even when species show large phylogenetic divergences, hybridization is relatively common in areas of secondary contact (see, for example, Bierne ; Willis ; Harper and Hart, 2007; Nydam and Harrison, 2011). In particular, hybridization/introgression among sponges has long been suspected and previously suggested in closely related species with sympatric distributions (Dailianis ; Escobar ). In this context, assessing the effect of potential hybridization on overall genetic diversity of both Ircinia species is fundamental to understanding their genetic structure and connectivity.

The main objectives of this study are twofold: (1) to describe the genetic structure of I. fasciculata across the Mediterranean and Adriatic seas in order to reveal connectivity patterns among populations and detect the potential effects of mass mortalities on the effective population size and (2) to identify whether hybridization occurs between I. fasciculata and I. variabilis and, if it does, to which extent.

Materials and methods

Sampling and DNA extraction

We collected ∼2 cm3 of tissue from 194 specimens of I. fasciculata of similar size (>500 cm3) between 2010 and 2013 in 12 locations across the Mediterranean and Adriatic seas (see details in Table 1). Similar sizes were selected in an effort to avoid sampling sponges from different cohorts or sponges that had recently been diseased (small fragments sometimes remained attached to the rock if the disease did not kill the whole sponge and could be confused with juveniles). Sampling was performed by SCUBA diving within the first 10 m of each location. Only six healthy individuals were sampled from Scandola Marine Reserve (Table 1) given the extreme decimation this population experienced from 2008 to 2010 (Cebrian ). Variation in geographical distance ranged from 21 to 3006 km for the studied area. Sponge fragments were preserved in absolute ethanol that was replaced with fresh ethanol at least three times within 48 h and stored at −20 °C until further processed. DNA was extracted with the DNeasy Blood & Tissue kit (Qiagen, Hilden, Germany) following the manufacturer's instructions with a minor modification concerning overall cell lysis time (that is, incubation was conducted overnight) and the final DNA elution step (performed twice using 50 μl of buffer EB each time).

Table 1

Details on collection sites and number of individuals per location for Ircinia fasciculata

Label	Location	Coordinates	N	Group
CAI	Es Caials, Cap de Creus, Spain	42° 17' 06.99" N, 3° 17' 48.57" E	14	Western Mediterranean
ESC	L'Illa Mateua, L'Escala, Spain	42° 6' 43.45" N, 3° 10' 23.68" E	16	Western Mediterranean
TOS	Mar Menuda, Tossa de Mar, Spain	41° 43' 17.81" N, 2° 56' 23.09" E	17	Western Mediterranean
BLA	S'Agulla, Blanes, Spain	41° 40' 54.16" N, 2° 48' 57.42" E	24	Western Mediterranean
CAL	Calafat, Tarragona, Spain	40º 55' 16.93" N, 0º 50' 30.61" E	12	Western Mediterranean
ALI	Benidorm Island, Vilajoiosa, Alicante, Spain	38º 30' 00.37" N, 0º 08' 12.38" E	10	Western Mediterranean
CAB	Illot L'Imperial, Cabrera National Park, Spain	39º 07' 33.14" N, 2º 57' 34.32" E	18	Pre-Balear
COR	Reserve Naturelle de Scandola, Corsica, France	42° 21′ 31.34" N, 8° 33′ 16.42" E	6	Pre-Balear
TAR	Las Palomas, Tarifa, Spain	36° 0′ 42.74′′ N, 5° 35′ 48.17′′ W	20	Alboran Sea
GRA	Piedra del Hombre, Almuñécar, Granada, Spain	36º 43' 16.0'' N, 3º 44' 13'' W	20	Alboran Sea
NAP	Bacoli, Napoli, Italy	40° 47′ 39.84′′ N, 14° 5′ 16.16′′ E	17	Thyrrenian Sea
CRO	Hvar Island, Croatia	43°78′ 22.92′′ N, 16° 43′ 38.42′′ E	20	Adriatic Sea

Microsatellite amplification and analysis

We genotyped all individuals of I. fasciculata from all locations at 14 microsatellite loci (3IFAS, 7IFAS, 8IFAS, 10IFAS, 22IFAS, 25IFAS, 36IFAS, 40IFAS, 41IFAS, 42IFAS, 44IFAS, 47IFAS, 48IFAS and 49IFAS) previously described elsewhere (Riesgo ), using the PCR conditions described therein. The sizes of the fluorescently labeled PCR products were estimated using an internal size marker (GeneScan 500 LIZ, Applied Biosystems, Foster City, CA, USA) on an ABI Prism 7700 Sequencer (Applied Biosystems) and analyzed with PeakScanner v1.0.

Genetic diversity in I. fasciculata populations

Estimations for the observed (Ho) and expected (He) heterozygosity, and the fixation index (FIS), commonly used as an inbreeding coefficient, were performed using GenAlEx 6.5 (Peakall and Smouse, 2006). We used Genepop web version 4.0.10 (Raymond and Rousset, 1995) to obtain values for departure of Hardy–Weinberg equilibrium (HWE) by locus and population (sampling site) using a probability test with level of significance determined by the following Markov chain parameters: 5000 dememorization steps, 1000 batches and 5000 iterations per batch. The total number of alleles per locus and population, allele richness, number of private alleles and genetic (gene) diversity (expected frequency of heterozygotes within subpopulations, including a correction for sampling bias stemming from sampling a limited number of individuals per population) were calculated with GENODIVE (Meirmans and Van Tienderen, 2004) and FSTAT 2.9.3.2 (Goudet, 1995). Tests for linkage disequilibrium for each locus were previously performed elsewhere (Riesgo ).

Population differentiation in I. fasciculata

We performed six different methods to assess population structure and differentiation in I. fasciculata. Four methods were based on distances: FST estimations, isolation by distance (IBD), BARRIER analyses and the analysis of the molecular variance (AMOVA), and two used a clustering approach: STRUCTURE and the discriminant analysis of principal components (DAPC). Given that the distance methods for measuring population dissimilarity are largely affected by the violation of the HWE (Waples, 2015), we removed five loci from the final data matrix that displayed clear HWE deviation in more than half of the analyzed locations (3IFAS, 10IFAS, 8IFAS, 47IFAS and 49IFAS; see Riesgo ). Thus, all distance-based analyses were computed on a subset of data with the remaining nine microsatellites that did not show deviation from the HWE.

We first assessed the occurrence of population subdivision in large clusters using a Bayesian clustering approach in STRUCTURE (Pritchard ), a software that assigns individuals to clusters probabilistically, under an admixture model and with allele frequencies correlated between samples. The program was run with a burn-in time of 50 000 repetitions and 150 000 iterations (Markov chain Monte Carlo), setting the putative K (predicted number of populations) from 1 to 15 and 10 replicates for each run. The estimation of log probabilities of data Pr(X|K) for each value of K was evaluated by calculating ΔK that accounts for the rate of change in the log probability of data between successive K values (Evanno ). ΔK is currently considered a more reliable predictor of the true number of populations (Evanno ). Calculations and evaluation of ΔK were performed with Structure Harvester (Pritchard ) and graphs were visualized in STRUCTURE (Pritchard ). The CLUMPAK web server (Kopelman ) was then used to find the best alignment of the results from STRUCTURE across the range of K values by averaging the probabilities of each K cluster, enabling the production of a figure with the entire range of K values.

Population differentiation was estimated with FST statistic between pairwise sampling sites using an infinite allele model and the software Arlequin (Excoffier ). We only calculated the FST statistic and not RST because the former performs better when sample size is <50 individuals per site or <20 loci are analyzed (Gaggiotti ). Significance of FST values was evaluated by performing 20 000 permutations in Arlequin and corrected based on the false discovery rate approach described in Narum (2006). The frequency of null alleles was estimated using the expectation maximization algorithm (Dempster ) in FreeNA (Chapuis and Estoup, 2007). As the presence of null alleles in well-differentiated populations is known to yield an overestimation of population differentiation (Chapuis and Estoup, 2007), we repeated our analysis Excluding Null Alleles (Chapuis and Estoup, 2007). In all cases, the corrections only affected the fourth or fifth decimal place in the FST value and consequently the effect of presence of null alleles was disregarded.

To determine the spatial patterns driving differentiation among populations, linearized pairwise FST estimates (FST/1−FST) obtained for nine microsatellite loci (see above) were correlated against log-transformed geographical distances between samples and IBD using a Mantel test in Arlequin (Excoffier ). Geographical distances were estimated as the minimum linear distance between pairs of locations by sea. We tested IBD for the whole data set, and also separately for a data set containing the locations on the Western Mediterranean (CAI, ESC, TOS, BLA, CAL and ALI) and the Pre-Balear island locations (CAB and COR), between which potential dispersal barriers exist (that is, oceanic currents). The effect of barriers in determining the genetic structure of I. fasciculata populations was further evaluated using pairwise FST values and visualized with the software BARRIER v2.2 (Manni ). This method links a matrix of geographical coordinates with their corresponding distance matrix (FST), and applies the Monmonier's maximum distance algorithm to identify ‘barriers' to gene flow among sites, namely the zones where differences between pairs of sites are the largest.

We also performed DAPC with the adegenet package (Jombart, 2008) implemented in R (https://www.r-project.org/). DAPC defines clusters using the clustering algorithm k-means on transformed data with principal component analysis. The algorithm k-means is then run sequentially with increasing values of k, and different clustering solutions are compared using Bayesian information criterion. The optimal cluster solution should correspond then to the lowest value of Bayesian information criterion. We applied the DAPC analysis for both the complete matrix of all sampling sites and a reduced set containing only those from the Western Mediterranean (CAI, ESC, TOS, BLA, CAL and ALI) and the Pre-Balear island (CAB and COR) groups.

AMOVA was used to determine the hierarchical distribution of genetic variation. To run this analysis, we a priori defined five groups according to biogeographical areas previously defined by Spalding and partially isolated by oceanographic fronts as found by Schunter and Villamor : Alboran Sea, Western Mediterranean, Pre-Balear islands, Thyrrenian Sea and Adriatic sea (see Table 1). The significance of the AMOVAs was calculated with 16 000 permutations of the original data.

Population assignment and migration in I. fasciculata

In addition, we performed a population assignment analysis calculating the likelihood ratio thresholds for all 12 populations based on the Monte Carlo test with an α of 0.002 and 1000 replicated data sets using GENODIVE version 2.0b23 (Meirmans and Van Tienderen, 2004). This method assigns or excludes reference populations as possible origins of individuals on the basis of multilocus genotypes. Genetic assignment methods allow inferring where individuals originated, providing estimates of real-time dispersal through the detection of immigrant individuals. These analyses are useful in addressing relationships and structure when genetic differentiation at the population level is low. In addition, the detection of last-generation migrants was performed based on the calculations of the likelihood of an individual belonging to a given population that was then done replacing the zero frequencies by a random 0.005 frequency (estimated to outperform tests) in 4000 permutations in GENODIVE.

Demographic effects of mass mortalities in I. fasciculata

In order to detect recent effective population size reductions (bottlenecks) from allele data frequencies in our localities because of mass mortalities we used the software BOTTLENECK (Cornuet and Luikart, 1996). This program is based on the assumption that populations that have gone through a recent reduction of their effective population size show a reduction of the allelic diversity and heterozygosity, even though the allele frequencies are reduced faster than the heterozygosity. The software computes the distribution of the heterozygosity expected from the number of alleles for each population and locus, and it takes into account the sample size under the assumption of mutation–drift equilibrium using three different models: infinite allele model, two-phase model and the stepwise mutation model. The analysis is performed using three tests: a ‘sign test', a ‘standardized differences test' (Cornuet and Luikart, 1996) and ‘Wilcoxon sign-rank test' (Luikart and Cornuet, 1998).

Phylogenetic relationships between I. variabilis and I. fasciculata

In order to shed light into the phylogenetic relationships between I. variabilis and I. fasciculata, we amplified a fragment of the mitochondrial cytochrome c oxidase I (COI) gene, increasing the taxon sampling from only 3 specimens of I. fasciculata (Erwin ) to 15, and including all sequenced species of the family Irciniidae (see below). A fragment of 1093 bp of COI was amplified for 12 specimens of I. fasciculata from the most distant sites (CAL, n=1; ESC, n=2; TAR, n=1, TOS, n=2; NAP, n=4; CRO, n=2) using the primers dgLCO1490 (Meyer ) and COX1-R1 (Rot ). COI sequences of I. variabilis and I. fasciculata from BLA were obtained from GenBank (see accession numbers in Figure 1a). Reactions were performed using 0.5 μl of total DNA (∼0.5–1 ng of DNA), 5.1 μl of pure distilled water, 3.3 μl of Extract-N-Amp PCR Reaction Mix (Sigma-Aldrich, St Louis, MO, USA) and 0.3 μl of each primer (10 μM). The amplification program was as follows: denaturation at 94 °C for 2 min; 30 cycles of denaturation at 94 °C for 60 s, annealing at 45 °C for 60 s and extension at 72 °C for 60 s; with a final extension step at 72 °C for 7 min. PCR products were purified and sequenced using the BigDye Terminator v3.1 cycle sequencing kit (ABI, Foster City, CA, USA), an annealing temperature of 50 °C and using the same primers as for the amplification step. Samples were analyzed on an Applied Biosystems 3730xl Genetic Analyzer available at the University of Barcelona Scientific and Technological Centers. Sequences were edited using Geneious 8.0 (Kearse ) and aligned with MAFFT (Katoh ). The sequence data obtained in this study have been deposited in the NCBI (National Center for Biotechnology Information) GenBank database (accession numbers in progress). Additional sequences of species from the family Irciniidae were obtained from GenBank: I. oros (JN655186, JN655189), I. strobilina (JX306089), I. felix (JX306085), I. ramodigitata (FN552859), I. irregularis (FN552858), Ircinia sp. (LN828727), Sarcotragus spinulosus (HG816026, HG816027, HE591460), Psammocinia bulbosa (FN552812) and Psammocinia sp. (FN552837, FN552855). We included as outgroups Spongia officinalis (HQ830364) and Strepsichordaia sp. (FN552875) that are also members of the order Dictyoceratida, but not from the family Irciniidae.

Figure 1

Assignment of individual genotypes of I. fasciculata to genetically similar clusters (K) as inferred by STRUCTURE for all studied locations with (a) K=2 and (b) K=5. (c) Pie frequency charts showing percentage of populations assigned to each of the five clusters represented in (b) for all locations. Gray color is applied to the percentage of individuals that were assigned to more than one cluster. Pre-Bal, Pre-Balear islands; Thyrre, Thyrrenian Sea; Adria, Adriatic Sea. A full color version of this figure is available at the Heredity journal online.

The evolutionary model that best fitted our data was obtained with jModeltest2 (Darriba ), and evaluated using the Akaike information criterion. The phylogenetic tree was reconstructed using maximum likelihood in PhyML (Guindon and Gascuel, 2003) with 500 bootstrap replicates, 10 initial trees and the evolutionary model HKY85 (see above). Bayesian inference analysis was performed in MrBayes v3.2 (Ronquist and Huelsenbeck, 2003) with two runs, each with three hot chains and one cold chain with 20 million generations, sampling every 2500th generation, using random starting trees. The remaining trees were combined to find the maximum a posteriori probability estimate of phylogeny. This analysis was conducted twice and convergence was tested using Tracer v1.5 (Rambaut and Drummond, 2007). The relationships among the COI sequence haplotypes of both I. fasciculata and I. variabilis were inferred using the TCS method in the program TCS (Clement ) as implemented in Popart 1.7 (http://popart.otago.ac.nz).

Hybridization between I variabilis and I. fasciculata: microsatellite analyses

We performed a fine-scale analysis of the introgression processes between our two target Ircinia species by genotyping 14 individuals from Tossa de Mar, Corsica and Cabrera of I. variabilis (TOS: n=5, COR: n=6 and CAB: n=3, respectively), and randomly selected 36 individuals from the same sampling sites of I. fasciculata that were already genotyped (TOS: n=10, COR: n=6, and CAB: n=10) using 12 microsatellite loci (all listed above except for 7IFAS and 40IFAS, as we were not able to amplify them in I. variabilis). We subsampled TOS and CAB for I. fasciculata in order to maintain the population sizes as even as possible, given that unbalanced sample sizes could hamper the recovery of population structure (Puechmaille, 2016). These three localities (TOS, COR and CAB) were selected because both species live in sympatry. Hybrid detection was achieved by performing an individual assignment to each cluster using STRUCTURE (Pritchard ), running the program with a burn-in time of 50 000 repetitions and 150 000 iterations (Markov chain Monte Carlo), and setting the putative K (predicted number of genetic units) from 1 to 6, with 10 replicates for each K. The estimation of log probabilities of data and evaluation of ΔK were calculated as described above. We then used CLUMPACK (Kopelman ) to find the best alignment of the results across the range of K values. Genetic differentiation between species and sites was visualized with a DAPC and the adegenet package (Jombart, 2008) implemented in R.

Results

Population structure and connectivity of I. fasciculata

None of the 194 individuals analyzed showed identical genotypes, and therefore clonality can be regarded as negligible in the studied populations. The total number of alleles ranged from 2 to 17 per population, and the average number of alleles per population ranged from 3.857 (COR) to 8.000 (CRO) (Supplementary Table 1). Similarly, the average frequency of private alleles differed from 0.143 (COR) to 1.500 (CRO). Private alleles were present in all populations varying from 1 (ALI, average pA frequency=0.429) to a maximum of 4 in CRO (average pA frequency=1.5; Table 2).

Table 2

Descriptors of genetic diversity for all 12 locations of Ircinia fasciculata using the data set containing 14 loci and the corrected data set with only 9 loci

Pop.	N	A		rA		pA		Ho		He		FIS		HWE
		14 Loci	9 Loci	14 Loci	9 Loci	14 Loci	9 Loci	14 Loci	9 Loci	14 Loci	9 Loci	14 Loci	9 Loci	14 Loci	9 Loci
TAR	17	5.857	7.111	4.683	4.666	1.071	0.778	0.587	0.621	0.643	0.64	0.12	0.016	***	NS
GRA	20	5.857	5.556	4.129	3.918	0.786	0.778	0.482	0.524	0.613	0.595	0.245	0.125	***	NS
CAI	14	5.928	5.667	4.683	4.255	1.071	0.444	0.724	0.738	0.652	0.619	−0.074	−0.222	*	NS
ESC	16	5.714	5.778	4.228	4.183	0.786	0.667	0.599	0.610	0.622	0.595	0.073	−0.058	*	NS
TOS	17	6.357	5.889	4.192	3.944	0.857	0.778	0.686	0.699	0.609	0.565	−0.097	−0.221	***	NS
BLA	24	6.857	6.444	4.169	4.055	0.429	0.556	0.638	0.676	0.625	0.604	0	−0.142	NS	NS
CAL	12	7.285	7.556	5.315	5.424	0.857	0.778	0.659	0.645	0.695	0.682	0.096	0.031	**	NS
ALI	10	5.428	5.444	4.507	4.506	0.429	0.556	0.614	0.644	0.648	0.632	0.106	−0.056	*	NS
CAB	18	6.571	5.778	4.537	4.116	0.643	0.444	0.607	0.627	0.658	0.617	0.107	−0.044	***	NS
COR	6	3.857	3.889	3.857	3.889	0.143	0.222	0.655	0.722	0.611	0.608	0.02	−0.233	NS	NS
NAP	20	7.928	8.333	4.876	5.069	0.786	0.889	0.607	0.648	0.672	0.694	0.126	0.013	***	NS
CRO	20	8.000	8.111	4.737	4.735	1.500	1.778	0.648	0.691	0.655	0.638	0.039	−0.103	NS	NS

Abbreviations: A, number of alleles; FIS, inbreeding coefficient; He, expected heterozygosity; Ho, observed heterozygosity; HWE, deviation from Hardy–Weinberg equilibrium; N, number of individuals; NS, not significant; pA, number of private alleles; Pop., population; rA, number of alleles after rarefaction.

*P<0.05; **P<0.01; ***P<0.001.

Genetic (gene) diversity values were very similar across populations, ranging from 0.626 in TOS to 0.728 in CAL (Supplementary Table 1). Expected heterozygosity values (He) showed the same trend as the genetic diversity values in all populations when using all 14 loci and the subset of 9 loci (Table 2). When the matrix containing all 14 loci was analyzed, fixation indexes (FIS) were positive in all cases except for CAI and TOS. When considering the 14 microsatellite markers, all populations except for BLA, COR and CRO deviated from HWE, whereas with the 9 microsatellites subset all populations remained in HWE (Table 2 and Supplementary Table 1).

The optimal number of clusters for the whole data set by the method of Evanno ) in the Bayesian analysis with STRUCTURE was two genetically homogenous groups (K=2) being followed by five (K=5; Figures 1a and b; Supplementary Figure 1). Results for K=2 revealed the division into two differentiated clusters without a clear pattern of subdivision: one cluster corresponded to the populations TAR and GRA (Alboran Sea), as well as some individuals from CAL, CAB, CAL, CAB, NAP and CRO (Table 1), and another cluster containing most individuals sampled from Western Mediterranean locations (CAI, ESC, TOS, BLA, CAL and ALI) and the Pre-Balear island locations (CAB and COR; Figure 1a). When the number of clusters K was selected by the highest likelihood (K=5; Supplementary Figure 1) the populations revealed inner substructuring (Figures 1b and c): (1) the first group (green) was formed mostly by specimens from the east Mediterranean and half of the individuals of NAP, and CAL populations; (2) the second group (yellow) was formed by several individuals from three different geographical areas, half of the populations of TAR, TOS, NAP and CRO; (3) the third group (pink) was formed by half population of TAR, the entire GRA population and some individuals in NAP; (4) the fourth group (blue) comprised CAI, BLA, COR and half of the populations of ESC, TOS, CAL and CAB; and finally (5) the fifth group (red) comprised half ESC and CAB populations (Figures 1b and c). It is noteworthy that in most populations individual genotypes were assigned to a specific cluster, and ‘mixed' genotypes from different clusters were detected at low frequency.

FST pairwise comparisons using 9 loci were always significant (Supplementary Table 2), except for comparisons between CAI and ESC (FST=0.01905), ESC and ALI (FST=0.06504), TAR and CAB (FST=0.00336) and ESC and CRO (FST=0.01581), indicating significant differentiation in the genetic structure among most of the studied populations.

The Mantel tests (Figure 2) revealed IBD only when the Western Mediterranean and the Pre-Balear island locations were analyzed separately (Figure 2b), but it was not significant when all populations were analyzed together (Figure 2a). Indeed, several barriers, which could potentially be oceanographic, were detected and appeared to contribute to population differentiation more than distance (Figure 3a). In our case, three barriers were set a priori given the three major oceanographic waterfronts known in the Mediterranean (Figure 3b). The Monmonier's maximum distance algorithm detected these three barriers and denoted their rank of importance (Figure 3a): first the Almeria-Oran Front (a) separating the Alboran Sea (TAR and GRA) locations from the rest, then the Ligurian-Thyrrenian barrier (b) separating Thyrrenian and Adriatic locations (NAP and CRO) from the western ones, and finally the North-Balearic Front (c) that also separated BLA and CAL from the rest of the West Mediterranean populations.

Figure 2

IBD analyses in I. fasciculata for (a) all sites and (b) only Western Mediterranean and Pre-Balear sites.

Figure 3

Detection of barriers to gene flow in I. fasciculata across the Mediterranean Sea. (a) Barriers detected by software BARRIER ranked a to c in order of magnitude. (b) Map of the Mediterranean seawater circulation indicating major fronts and oceanographic barriers previously known in the Mediterranean. Note that sampling locations for I. fasciculata are indicated in the map.

The DAPC showed separation between the Alboran Sea (TAR and GRA), the Pre-Balear islands (COR and CAB), the Thyrrenian and Adriatic locations (NAP and CRO) and the remaining locations of the Western Mediterranean (Figure 4a). Indeed, the northern locations of the Western Mediterranean (CAI and ESC) and the Pre-Balear islands (COR and CAB) grouped more closely together, being slightly separated from the rest of populations of the Western Mediterranean and Thyrrenian and Adriatic Sea group (Figure 4a). In addition, when analyzing the spatial structure of Western Mediterranean and Pre-Balear island locations, only CAL was clearly separated from the rest of populations (Figure 4b).

Figure 4

Subdivision of the I. fasciculata sites according to DAPC analysis for (a) all areas and (b) only Western Mediterranean and Pre-Balear sites.

When analyzing population differentiation in the different regions using AMOVA, we found significant genetic differentiation among the groups containing Alboran Sea, Western Mediterranean, Pre-Balear islands and Thyrrenian and Adriatic Sea locations (Tables 1 and 3). In addition, significant genetic differentiation among the populations within these groups was also observed (Table 3).

Table 3

Details on the results of the analysis of the molecular variance (AMOVA) performed using three groups for the collection sites of Ircinia fasciculata: Alboran Sea (TAR and GRA), Western Mediterranean (CAI, ESC, TOS, BLA, CAL and ALI), Pre-Balear islands (CAB and COR), Thyrrenan Sea (NAP) and Adriatic Sea (CRO)

Source of variation	Sum of squares	Variance components	Percentage variation	Fixation indices	P-value
Among groups	18.7940	0.0427	1.9200	0.01923	0.01214
Among populations within groups	38.5160	0.0765	3.4500	0.03516	0.00000
Among individuals within populations	348.9400	−0.1818	−8.2000	−0.08663	1.00000
Within individuals	442.5000	2.2809	102.8300	2.28093	0.99496

See Table 1 for full name of locations. Significant P-values are in bold.

The population assignment analysis showed large genetic exchange within each of the groups (Alboran Sea, Western Mediterranean, Pre-Balear islands, Thyrrenian and Adriatic Seas), with few contributions outside the groups, except for CAI and BLA that seemed to have contributed genotypes to almost all locations in the sampled region (Figure 5a). The analysis only detected 5 migrants from the last generation, 1 individual from GRA comprising 15% of NAP, 1 individual from ESC comprising 16.7% of COR, 1 individual from CAL and another from CRO being 20% of ALI and finally 1 individual from TAR contributing to 5% of GRA (Figure 5b).

Figure 5

Population assignment of individuals and last-generation migrants in I. fasciculata. (a) Number of individuals assigned to each population for all sampling locations. (b) Number of last-generation migrants (indicated in percentage of population). The original location of migrants (origin) is color coded following the same scheme as in DAPC of Figure 4a. A full color version of this figure is available at the Heredity journal online.

Demographic effects of mass mortalities in Ircinia fasciculata

All three BOTTLENECK tests identified a population reduction in COR under the infinite allele model (Table 4), but also detected population reductions in other locations depending on the model used. Using two-phase model, the sign test did not identify any bottleneck in any of the studied locations using the Wilcoxon rank test (Table 4), whereas the standardized differences test identified an effective size reduction in BLA, TOS, TAR, GRA, NAP and CRO with the same model (Table 4). Using the stepwise mutation model, the standardized differences test identified effective size reductions in most sites except for COR (Table 4), even though using the same model, the Wilcoxon rank test identified bottlenecks in CAB, BLA, ESC, TAR, NAP and CRO (Table 4).

Table 4

Results of the Bottleneck analysis for the 12 locations (Loc) of Ircinia fasciculata

Loc	Mean N	Mean k	Mean He	Sign test			Standardized differences test			Wilcoxon rank test
				IAM	TPM	SMM	IAM	TPM	SMM	IAM		TPM		SMM
				P	P	P	P	P	P	P 1t	P 2t	P 1t	P 2t	P 1t	P 2t
CAI	28.00	5.93	0.67649	0.105687	0.459598	0.059423	0.153783	0.292739	0.002558	0.178772	0.357544	0.403870	0.807739	0.975281	0.057983
ALI	19.57	5.43	0.68338	0.358475	0.067512	0.060697	0.344748	0.011027	0.000002	0.500000	1.000.000	0.903137	0.216553	0.987732	0.029541
COR	12.00	3.86	0.66667	0.001685	0.019651	0.424299	0.006473	0.028106	0.106599	0.002686	0.005371	0.014771	0.029541	0.052002	0.104004
CAB	34.14	6.57	0.67806	0.54355	0.06084	0.004096	0.424557	0.018419	0.000000	0.291504	0.583008	0.903137	0.216553	0.999420	0.001526
BLA	47.86	6.86	0.63841	0.249174	0.342223	0.004468	0.47137	0.002866	0.000000	0.380432	0.760864	0.821228	0.390991	0.997864	0.005371
ESC	29.29	5.71	0.64423	0.549849	0.160654	0.003971	0.317263	0.129642	0.000117	0.231567	0.463135	0.852112	0.325806	0.997864	0.005371
CAL	23.86	7.29	0.72507	0.545963	0.302962	0.004117	0.450574	0.064873	0.000037	0.334900	0.669800	0.866211	0.295776	0.991699	0.020264
TOS	33.43	6.36	0.62750	0.453395	0.528966	0.021274	0.11262	0.000098	0.000000	0.665100	0.714844	0.729187	0.583008	0.993286	0.016602
TAR	32.43	7.07	0.66405	0.543737	0.161753	0.004721	0.143488	0.000129	0.000000	0.768433	0.501587	0.947998	0.118896	0.999420	0.001526
GRA	34.86	5.86	0.63206	0.361832	0.176926	0.169848	0.408537	0.002919	0.000000	0.524231	1.000.000	0.932373	0.153076	0.975281	0.057983
NAP	37.00	7.93	0.69106	0.549009	0.067966	0.005446	0.138564	0.000009	0.000000	0.451599	0.903198	0.923462	0.172607	0.997864	0.005371
CRO	38.14	8.00	0.67319	0.322208	0.062244	0.000072	0.185037	0.000266	0.000000	0.548401	0.951538	0.982361	0.041870	0.999908	0.000305

Abbreviations: He, expected heterozygosity; IAM, infinite allele model; k, number of alleles; N, sample size; SMM, stepwise mutation model; TPM, two-phase model.

See Table 1 for full name of locations. Significant P-values are in bold.

Phylogenetic relationship between I. variabilis and I. fasciculata

The 19 sequences of the COI marker analyzed from 8 different locations of I. variabilis and I. fasciculata were 1093 base pairs long and, among them, 28 sites showed polymorphism (5.004%). Within the reduced phylogeny of the genus Ircinia, both the maximum likelihood and Bayesian inference analyses showed a large, well-supported clade containing sequences of both I. variabilis and I. fasciculata (Figure 6a). Within this clade, sequences showed no clear grouping pattern between species or among locations (Figure 6a). The TCS haplotype network showed that the main haplotype for I. fasciculata was recovered in three locations (BLA, CAL and NAP) and this same haplotype was shared with three individuals of I. variabilis from BLA (Figure 6b). In addition, seven more haplotypes were recovered for I. fasciculata (Figure 6b).

Figure 6

Hybridization patterns in I. variabilis and I. fasciculata. (a) Maximum likelihood phylogenetic phylogeny based on the COI sequences from the genus Ircinia with bootstrap values (above nodes) and posterior probabilities (below nodes) mapped on shared clades between maximum likelihood and Bayesian inference analyses. Support values over 70 are located over each node. Sequences obtained in this study appear in bold letters. The abbreviations for sampling sites can be found in Table 1. (b) TCS haplotype network for I. fasciculata. (c) DAPC graphical representation of population differentiation in three locations (TOS, COR and CAB) containing both species (I. variabilis and I. fasciculata) using microsatellites. (d) STRUCTURE output of individual assignment to five clusters (K=5) in three locations (TOS, COR and CAB) containing both species (I. variabilis and I. fasciculata) using microsatellites.

Hybridization patterns in I. variabilis and I. fasciculata

The DAPC analysis using microsatellite data of both species collected in the same three sites showed a clear group containing all locations of I. variabilis (IVTOS, IVCAB and IVCOR) and those individuals of I. fasciculata collected in Tossa de Mar (IFTOS), whereas the rest of locations of I. fasciculata appeared separated from them (IFCAB and IFCOR, Figure 6c). In addition, the individual assignment analysis performed in STRUCTURE showed that individuals from the locations TOS and CAB of I. variabilis belonged to the same group (blue) as some individuals of I. fasciculata in the sites TOS, COR and CAB (Figure 6d). Similarly, two individuals in the location of COR of I. variabilis appeared as a ‘mixture' of I. variabilis clusters (Figure 6d). The crimson cluster mostly belonged to individuals of I. fasciculata at the locations of TOS and CAB (Figure 6d). All our results indicate some degree of hybridization between I. variabilis and I. fasciculata, at least in some locations, and especially in TOS.

Discussion

Genetic diversity in I. fasciculata

Genetic diversity (calculated with the parameters gene diversity in GENODIVE and He in GenAlex) was ∼0.6, a value that is similar to those observed in other Mediterranean sponges, including Crambe crambe (∼0.6, Duran ), the endangered species Spongia lamella (∼0.5, Pérez-Portela ), the bath sponge S. officinalis (∼0.8, Dailianis ) and Scopalina lophyropoda (∼0.4; Blanquer and Uriz, 2010). Inbreeding analysis detected a significant heterozygosity deficit in most localities, indicating that they were not in HWE, a situation widespread among most studied sponges (Duran ; Blanquer and Uriz, 2010; Dailianis ; Pérez-Portela ; Chaves-Fonnegra ; Giles ) and many other marine benthic invertebrates in the same geographical area (see, for example, Addison and Hart, 2005; Pérez-Portela and Turon, 2008; Mokhtar-Jamaï ; Pérez-Portela ). One of the reasons for deviations from HWE in marine invertebrates may be associated with the production of planktonic sperm for fertilization (Addison and Hart, 2005), the reproductive strategy of sponges (Maldonado and Riesgo, 2008). It has been suggested that in free spawners, higher mutation rates caused by larger numbers of cell cycles leading to high sperm production would be consistent with more frequent null alleles, causing greater departures from HWE in comparison with copulating species with lower allocation to sperm production (see Addison and Hart, 2005 for discussion). However, as our FreeNA analysis did not detect any significant effect of null alleles, other nonexclusive factors should be considered to explain heterozygote deficiencies, including high levels of inbreeding, selection against heterozygotes, Wahlund effect or a combination of all or some of these factors (Freeland ). Inbreeding is indeed widespread in most benthic invertebrates with a lecithotrophic larva and very limited dispersal capabilities (see, for example, Goffredo ; Mariani ; Pérez-Portela and Turon, 2008; Ledoux ; Mokhtar-Jamaï ; Pérez-Portela ). For instance, the mean distance between the parents and offspring in the red coral Corallium rubrum was estimated to be 20–35 cm (Ledoux ), a similar value to that reported for the Mediterranean sponge C. crambe (Calderon ). On the other hand, genetic substructuring within I. fasciculata populations because of the presence of different cohorts and family clustering (Wahlund effect) is possible as sponges are long lived and particularly I. fasciculata possesses highly philopatric larvae (Mariani ). In summary, both inbreeding and the existence of substructure might be the factors that most likely are causing the deviation of HWE in our I. fasciculata populations.

Population structure and connectivity in I. fasciculata

The populations of I. fasciculata showed genetic structure at both large and small geographical scales. First, two major clusters were detected across the sampling site covered in this study, one grouping the Alboran Sea locations and the Thyrrenian and Adriatic Sea locations, and the other one clustering the locations along the Western Mediterranean coast of Spain and the Balearic and Corsica Islands (Figure 2a). However, the Mantel test performed along the entire sampling range of our study did not reveal a pattern of stepping stone gene flow (derived from a process of IBD), suggesting that other major processes such as the existence of oceanographic barriers (for example, major currents) may be driving population differentiation patterns in this area. In fact, three barriers to gene flow were detected in our data set: the strongest barrier appeared to separate the Alboran Sea locations from those in the rest of the Mediterranean Sea (presumably the Almeria-Oran Front, Figure 4), and this was also correlated with the limited gene flow detected between these areas and the high population differentiation observed between both regions. For many marine organisms, the Almeria-Oran Front acts as a strong barrier for gene flow between Atlantic and Mediterranean populations (see review in Patarnello ), including sponges (see, for example, Pérez-Portela ). In contrast, stepping stone gene flow was observed when analyzing the spatial distribution of genetic divergence exclusively between the Western and Thyrrenian and Adriatic locations. The other two barriers detected corresponded to the North-Balearic Front and the Ligurian-Thyrrenian barrier (Figure 4). As found for other benthic organisms (Villamor ), the Ligurian-Thyrrenian barrier appeared to be a seasonal barrier, and it did not completely prevent the gene flow between populations across it.

Population differentiation was also detected at smaller scales: significant differences were detected between all pairwise comparisons except for the closely located sites of CAI and ESC in the North-Western coast of Spain, and the distantly located pairs of ESC and ALI, TAR and CAB and ESC and CRO (Figure 2c). Although high gene flow may be expected between CAI and ESC (only located 50 km apart), the extensive gene flow detected among the other 3 site pairs (located between 550 and 2250 km apart) does not have a straightforward explanation. In general, populations of sessile invertebrates with relatively low dispersal capabilities in the Mediterranean are highly structured like those of I. fasciculata, including not only sponges (Duran ; Blanquer and Uriz, 2010; Dailianis ; Pérez-Portela ), but also cnidarians and ascidians (Goffredo ; López-Legentil and Turon, 2006; Pérez-Portela and Turon, 2008; Mokhtar-Jamaï ). Even sponges from other distant regions show high values of population differentiation (Chaves-Fonnegra ; Giles ). The overall FST value of genetic differentiation found for I. fasciculata (FST=0.221) was higher than that reported for the closely related dictyoceratid S. officinalis (FST=0.061, Dailianis ), but similar to other Mediterranean demosponges studied using microsatellites (C. crambe: FST=0.18, Duran ; S. lophyropoda: FST=0.122, Blanquer and Uriz, 2010; S. lamella: FST=0.236, Pérez-Portela ) and even corals, such as Balanophyllia europaea (0.202; Goffredo ), indicating that the differentiation patterns in these Mediterranean benthic invertebrates might be equally driven by the inbreeding caused by the low dispersal abilities of their larvae, and the occurrence of oceanographic barriers to their gene flow.

Recent reductions in the effective population size (or bottlenecks) have been detected in almost all the studied sites of I. fasciculata except for CAI, ALI and CAL depending on the model used. During a bottleneck and under the infinite allele model, allelic diversity is reduced faster than the heterozygosity. However, under the strict stepwise mutation model and probably under the two-phase model, heterozygosity excess may not be observed (Cornuet and Luikart, 1996). In our case, the different results in the detection of bottlenecks under different models may be reflecting differences in the resolution of the analyses rather than inconclusive results. In fact, these reductions in population size can well be related to recent episodes of massive individual losses reported for I. fasciculata (see, for example, Garrabou ; Maldonado ; Cebrian ). Even though a massive episode of mortality was observed in ALI in 2008 (Maldonado ), a bottleneck was not detected with our genetic analysis. The most parsimonious explanations for our observation are that (1) ALI was unaffected by genetic drift after massive dye-offs and/or (2) that this location is well connected (by means of gene flow) with other nearby I. fasciculata populations. Indeed, the largest number of migrants from last generation was detected in ALI, suggesting higher migration rates toward this location.

Phylogenetic analyses and hybridization between I. variabilis and I. fasciculata

Previous analysis using COI sequences did not provide enough resolution to delimit I. variabilis and I. fasciculata species, but analyses of the nuclear ribosomal gene ITS-2 showed that they formed different clades with a divergence of 2.97%±0.38 (Erwin ). Similarly, our results show that COI did not separate I. variabilis and I. fasciculata because they form together a well-supported but unresolved clade, with shared haplotypes between both species. This was not surprising, as nucleotide substitution rates in the mitochondrial genome in sponges and other lower metazoans are up to 100 times slower than other animal lineages (Hellberg, 2006; Shearer ). In this sense, although COI has been previously found to be useful to explain speciation and phylogeographic patterns in some haplosclerid sponges (see, for example, López-Legentil and Pawlik, 2009; DeBiasse and Hellberg, 2015), it seems to be unsuitable for resolving phylogenetic relationships in other sponges (see, for example, Duran ; Dailianis ). In light of these results, two possible explanations arise: the first one is that they are separate but closely related species that have undergone secondary contact favoring introgression, and a second explanation may be that I. variabilis and I. fasciculata are the same species but are currently undergoing incipient speciation without reproductive isolation yet occurring, and therefore a certain level of mixing still exists. Both the reciprocal monophyly of ITS-2 sequence clades and the existence of some morphological differences between these two species lend support to our first hypothesis. In addition, it is important to note here that our molecular results do not support the reassignment of I. fasciculata to the genus Sarcotragus, as previously proposed by Pronzato , as I. fasciculata and the type species for Sarcotragus (S. spinulosus) appear in two separated and clearly divergent clades (Figure 6a).

Instead of using COI, several authors have suggested multilocus (for example, nuclear coding genes, amplified fragment length polymorphisms, microsatellites and single-nucleotide polymorphisms) approaches as potentially suitable for resolving species boundaries (Petit and Excoffier, 2009; Hausdorf and Hennig, 2010). In particular, microsatellite loci are characterized by relatively high mutation rates, an assumed neutral selection rate and expected lower introgression potential than mitochondrial markers for at least some organisms (Petit and Excoffier, 2009). Our results using microsatellite loci suggested that some degree of hybridization and introgression occur in at least some populations where both species of Ircinia appear in sympatry (Figures 6b and c). Introgression between closely related species with overlapping ranges and incomplete reproductive barriers (such as the Ircinia species here) is common in marine invertebrates (see, for example, Willis ) even when the divergence between species is relatively large (Nydam and Harrison, 2011). Interestingly, introgression could be increased when intraspecific gene flow decreases (Petit and Excoffier, 2009). Thus, as the number and densities of Ircinia populations decrease in the Mediterranean Sea, a concomitant decrease in gene flow among affected populations may be accompanied by higher hybridization rates between I. variabilis and I. fasciculata populations, further blurring these two species boundaries. Although the relationship between intra- and interspecific gene flow is mostly documented for human populations (Currat ), the hybridization patterns of I. fasciculata reported here might be a first step toward understanding the impact of introgression in the genetic diversity and structure of marine invertebrates after population reductions.

In conclusion, our results agree with previous studies on the genetic differentiation of sessile invertebrates with lecithotrophic larvae and limited dispersal potential in the Mediterranean Sea (Duran ; Goffredo ; Pérez-Portela and Turon, 2008; Blanquer and Uriz, 2010; Dailianis ; Mokhtar-Jamaï ; Pérez-Portela ), highlighting the important role of oceanographic barriers in shaping the genetic structure across their distribution (Patarnello ; Villamor ). In particular, we found that most locations of I. fasciculata are strongly genetically differentiated, with barriers preventing gene flow between many population pairs that correlate with well-known oceanographic barriers. Importantly, our results indicate that the genetic diversity of I. fasciculata is threatened by massive demographic losses in the Mediterranean Sea, and that its recovery is hindered by limited connectivity between populations. Demographic reductions of I. fasciculata in the Alboran sites might be especially worrisome, as these locations are genetically isolated from the rest. Interestingly, the hybridization potential observed between I. fasciculata and I. variabilis may allow the former to maintain a genetic pool large enough to palliate genetic diversity losses. Notably, increased introgression events among these species may offset reductions in intraspecific gene flow resulting from stress and mass mortality events.

Data archiving

All data sets have been archived in PANGAEA https://doi.pangaea.de/10.1594/PANGAEA.860018.