# Symbiosis, hybridization and speciation in Mediterranean octocorals (Octocorallia, Eunicellidae)

Didier Aurelle, Anne Haguenauer, Marc Bally, Frédéric Zuberer, Dorian Guillemain, Jean-Baptiste Ledoux, Stéphane Sartoretto, Cédric Cabau, Rachel Lapeyre, Lamya Chaoui, et al.

► To cite this version:

Didier Aurelle, Anne Haguenauer, Marc Bally, Frédéric Zuberer, Dorian Guillemain, et al. Symbiosis, hybridization and speciation in Mediterranean octocorals (Octocorallia, Eunicellidae). Biological Journal of the Linnean Society, 2024, 143 (4), 10.1093/biolinnean/blae116. hal-04859569

HAL Id: hal-04859569

https://hal.science/hal-04859569v1

Submitted on 30 Dec 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Symbiosis, hybridization and speciation in Mediterranean octocorals
(Octocorallia, Eunicellidae)

Didier Aurelle¹,²*, Anne Haguenauer³, Marc Bally¹, Frédéric Zuberer⁴, Dorian Guillemain⁴, Jean-Baptiste Ledoux⁵, Stéphane Sartoretto⁶, Cédric Cabau⁷, Rachel Lapeyre⁸, Lamya Chaoui⁹, Hichem Kara⁹, Sarah Samadi², Pierre Pontarotti¹⁰,¹¹,¹²

¹ Aix Marseille Univ, Université de Toulon, CNRS, IRD, MIO, Marseille, France

² Institut Systématique Evolution Biodiversité (ISYEB), Muséum national d'Histoire naturelle, CNRS, Sorbonne Université, EPHE, Université des Antilles, CP 26, 75005 Paris, France.

³ CNRS - Délégation Provence et corse, Marseille, France

⁴ Aix Marseille Univ, CNRS, IRD, INRAE, OSU Inst. PYTHEAS, Marseille, France

⁵ CIIMAR/CIMAR, Centro Interdisciplinar de Investigação Marinha e Ambiental, Universidade do Porto, Porto, Portugal.

⁶ Ifremer, LITTORAL, 83500 La Seyne-sur-Mer, France

⁷ Sigenae, GenPhySE, Université de Toulouse, INRAE, ENVT, 31326, Castanet Tolosan, France.

⁸ MGX-Montpellier GenomiX, Univ. Montpellier, CNRS, INSERM, Montpellier France

⁹ Laboratoire Bioressources marines. Université d'Annaba Badji Mokhtar, Annaba - Algérie.

¹⁰ Aix Marseille Univ, MEPHI, Marseille, France.

¹¹ IHU Méditerranée Infection, Marseille, France.

¹² CNRS SNC5039

Correspondence: didier.aurelle@univ-amu.fr

<figure><img src="image_2.png" /></figure>

CC-BY 4.0 https://creativecommons.org/licenses/by/4.0/

<footnote>*Corresponding author</footnote>
Understanding how species can form and remain isolated in the marine environment still stimulates active researches. Here we study the differentiation and the possibility of hybridization among three temperate octocorals: *Eunicella cavolini*, *E. singularis* and *E. verrucosa*. Morphologically intermediate individuals have been observed between them. Among these three species, *E. singularis* is the only one described in mutualistic symbiosis with photosynthetic Symbiodiniaceae. The symbiosis between Symbiodiniaceae and scleractinian corals is well studied, especially in the context of the response to anthropogenic climate change. Nevertheless, the potential role of symbiotic interactions in speciation processes remains unknown in cnidaria. We tested here the possibility of hybridization between symbiotic and non-symbiotic *Eunicella* species. Through multivariate analyses and hybrid detection, we prove the existence of on-going gene flow between *E. singularis* and *E. cavolini*, with the observation of F1 and F2 hybrids, and backcrosses. Demographic inferences indicate a scenario of secondary contact between these two species. Despite current gene flow, these two species appear genetically well differentiated. Our data also suggest an intermediate abundance of Symbiodiniaceae in the hybrids of the two parental populations. We discuss the evolution of the Symbiodiniaceae / cnidarian symbiosis in the light of our results.

**Keywords:** speciation, hybridization, symbiosis, transcriptome, RAD sequencing, octocoral
1 Introduction

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

18

19

20

21

22

23

24

25

26

27

28

29

30

31

32

33

34

35

36

37

As corner stones of evolutionary biology, species and speciation still raise a wealth of questions fuelled by the technological and conceptual advancements in genomics. Genomic data allow testing hypotheses about species boundaries and origins. Named species are indeed hypotheses, built on available data, that can be rejected or validated through the integration of additional data and / or the use of additional criteria based on evolutionary concepts (Pante et al., 2015b). Sound species delimitations are useful, among others, to better estimate species range and biodiversity patterns (Muir et al., 2022; Coelho et al., 2023), to avoid biases in studies of connectivity (Pante et al., 2015b), and of adaptive abilities (Brener-Raffalli et al., 2022). However, proposing sound species delimitation can be problematic because different delimitation criteria may bring contradictory conclusions about species boundaries (the Grey Zone of de Queiroz, 2007). This grey zone corresponds to puzzling cases such as the absence of gene flow among morphologically undifferentiated sets of organisms (i.e. cryptic species, Cahill et al., 2024), or conversely, the detection of gene flow among sets of organisms recognized, based on morphological distinctiveness, as distinct species (Leroy et al., 2020). Evolutionary inferences, based on genomic data, allow testing scenarios of speciation and current gene flow: this provides a better understanding on the origin and persistence of species at the light of genomic divergence (Roux et al., 2016; De Jode et al., 2023).

In the marine realm, the question of speciation is considered as particularly confusing. Notably, how new species can originate from populations with large effective size associated to high level of gene flow is still abundantly debated in the literature (e.g. Palumbi, 1992; Mayr, 2001; Faria et al., 2021). Difficulties in sampling and rearing organisms also hamper experiments to test reproductive isolation (Faria et al., 2021). Important progresses in methodologies now allow to better understand spatial patterns of genetic structure in marine organisms, for example through the study of oceanographic connectivity (Reynes et al., 2021), clines in allele frequencies (Gagnaire et al., 2015), and hybrid zones (Bierne et al., 2003).

In this context, the role of symbiotic interactions in reproductive isolation remains poorly investigated. There are various examples of the involvement of microbial species in reproductive isolation, especially in insects (Brucker and Bordenstein, 2012). For marine species, microbial communities have been mainly explored in light of adaptative evolution (Rosenberg and Zilber-Rosenberg, 2018). Shallow water scleractinian corals (hexacorals) are usually associated with various species of photosynthetic zooxanthellae, in the family Symbiodiniacae (Cairns, 2007; LaJeunesse et al., 2018). Changes in associated Symbiodiniaceae can impact the thermotolerance of the coral holobiont, and the possibility
38 of adaptation facing climate change (Berkelmans and van Oppen, 2006; van Oppen & 39 Medina, 2020). Inferences from the phylogeny of Anthozoans (hexacorals and octocorals) 40 have shown multiple acquisitions of the symbiotic state throughout evolution (Cairns, 2007; 41 Campoy et al., 2020, Mc Fadden et al., 2021). The symbiotic interactions between 42 Anthozoans and Symbiodiniaceae provide important mutualistic benefits especially from a 43 nutritional point of view (Furla et al., 2005). These interactions require specific adaptations 44 for the animal host, as for example protection against oxygen produced by photosynthesis 45 (Furla et al., 2005). Therefore, one can hypothesize that in hybrids such adaptations could 46 be modified and a breakdown of symbiosis could occur, leading to reduced fitness. The 47 association with Symbiodiniaceae can range from mutualism to parasitism (Sachs and 48 Wilcox, 2006; Lesser et al., 2013; see also Matz, 2024), and a change in the genomic 49 background in hybrid hosts could modify the nature of symbiosis as well. The presence of 50 Symbiodiniaceae could then be involved in genetic incompatibilities with the host genome, 51 as previously observed with bacterial species (Bordenstein, 2003; Brucker and Bordenstein, 52 2012). All these observations raise the question of the potential role of Symbiodiniaceae in 53 speciation and reproductive isolation in Anthozoans. This topic has been poorly explored up 54 to now. In *Plexaura* octocorals, two incompletely isolated species have been shown to 55 present different populations of Symbiodiniaceae, questioning their role in species 56 boundaries (Pelosi et al., 2020).

57 Here we explore the robustness of species limits between named species of the gorgonian 58 genus *Eunicella* (Octocorallia, Eunicellidae) documented as displaying different symbiotic 59 relationships. In shallow conditions (above 50 m depth), three *Eunicella* species are mainly 60 present in the Mediterranean Sea: *Eunicella cavolini* (Koch, 1887), *E. singularis* (Esper, 61 1971), and *E. verrucosa* (Pallas, 1766). These three species have partially overlapping 62 ranges, and they can be observed in sympatry in the area of Marseille (France). *Eunicella* 63 singularis hosts Symbiodiniaceae corresponding to the *Philozoon* genus (Forcioli et al., 64 2011; LaJeunesse et al., 2018, 2022; Porro, 2019), whereas the two other gorgonian species 65 are devoided of these symbionts (Carpine and Grasshoff, 1975). The Symbiodiniaceae 66 contribute to the carbon metabolism of *E. singularis*, but a non-symbiotic *aphyta* morph has 67 already been observed (Gori et al., 2012). The lack of variability in mitochondrial DNA does 68 not allow to distinguish these three species (Calderón et al., 2006), and a study using two 69 nuclear introns suggested the possibility of hybridization between *E. cavolini* and 70 *E. singularis* (Aurelle et al., 2017). Moreover, demographic inferences based on a large 71 number of nuclear loci in *E. cavolini* and *E. verrucosa* indicated the possibility of current 72 gene flow between these two species (Roux et al., 2016). However, these data are 73 incomplete because neither individuals identified as *E. singularis*, nor individuals that are 74 morphologically difficult to attribute to a named species (potential hybrids) have been
75

analysed. Here, we will go further on these topics with the following objectives: i) estimate
the genomic differentiation among these three species and test for species limits, ii) test
whether symbionts are present or absent in the hybrids, to look for a possible breakdown in
symbiosis, and iii) infer scenarios of speciation. Studying the history of speciation is useful to
infer how divergence happened, and to test the possibility of current and past gene flow.

76

77

78

79

80

Analysing the hybrid status on morphologically intermediate individuals allows to further test
if hybridization is still on-going. We used restriction sites associated DNA sequencing (RAD-
sequencing; Baird et al., 2008) to test species limits and hybridization. We complementary
used transcriptome data for demographic inferences, for the analysis of putative hybrids, and
to test for the presence of Symbiodiniaceae. The results will be useful to better understand
the evolution of these species in different environments and particularly the possible impact
of hybridization in adaptation to changing environment.

81

82

83

84

85

86

87

88

## Material and methods

89

90

## Species distribution

91

92

*Eunicella verrucosa* is present both in the Eastern Atlantic Ocean and the Mediterranean Sea (Carpine and Grasshoff, 1975). In the Atlantic, *E. verrucosa* can be found from Ireland and West coasts of Britain, to Angola (Grasshoff, 1992; Readman and Hiscock, 2017).

93

94

95

*Eunicella verrucosa* has been observed in the North Western Mediterranean Sea, in Sardinia (Canessa et al., 2022), and in the Adriatic and Aegean Seas (Chimienti, 2020). In the Mediterranean Sea, *E. verrucosa* can be observed from shallow conditions (20-40 m) up to 200 m depth (Sartoretto and Francour, 2011; Fourt and Goujard, 2012; Chimienti, 2020).

96

97

98

99

*Eunicella singularis* and *E. cavolini* are only present in the Mediterranean Sea. *Eunicella cavolini* can be observed in the Western Mediterranean, Adriatic and Aegean Seas, from 5 to 200 m depth (Sini *et al.*, 2015; Carugati *et al.*, 2022). *Eunicella singularis* can be found in the Western Mediterranean and Adriatic Seas, and, less frequently, in the Eastern Mediterranean (Gori *et al.*, 2012). It is usually observed up to 40 m depth. *Eunicella singularis* is the only Mediterranean octocoral known to harbour Symbiodiniaceae (but see Bonacolta *et al.*, 2024). These Symbiodiniaceae belong to the temperate clade A (Forcioli *et al.*, 2011; Casado-Amezúa *et al.*, 2016), now corresponding to the *Philozoon* genus (LaJeunesse *et al.*, 2018, 2022). Deep occurrences (up to 70 m) of *E. singularis* have been mentioned, and assigned to the *aphyta* morph, without Symbiodiniaceae (Gori *et al.*, 2012). In the area of Marseille, these three species can be observed in sympatry and at the same depth (Sartoretto and Francour, 2011).

100

101

102

103

104

105

106

107

108

109

110

111
112 **Sampling**

113

114 The sampling for RAD sequencing included 25 specimens identified as *E. cavolini*, 23
115 *E. singularis*, seven *E. verrucosa*, and 12 morphologically intermediate individuals (potential
116 hybrids). These latter individuals displayed intermediate colors and branching patterns
117 between *E. cavolini* and *E. singularis* (Figure S1), and they were analysed to test their hybrid
118 status (Aurelle *et al.*, 2017). The specimens have been sampled by scuba diving at different
119 periods in the area of Marseille, where the three species are present in sympatry (Figure S2;
120 Table S1).

121 For transcriptome sequencing, specimens attributed to *E. cavolini*, *E. singularis*, and
122 *E. verrucosa* have been collected in the Mediterranean (for the three species), and in the
123 Atlantic (*E. verrucosa* only; Table S2; Figure 1) in 2016. The final sampling for
124 transcriptomics included five *E. cavolini*, eight *E. singularis*, three *E. verrucosa*, and four
125 potential hybrids.

126

127 Sampling was non-destructive, with authorizations from the local authorities, including
128 Marine Protected Areas.

129

130 **Mitochondrial MutS**

131 To test the genetic proximity of three *Eunicella* species studied here, we built a tree with
132 mitochondrial MutS sequences (McFadden *et al.*, 2011), available in GenBank. The methods
133 and sequences are detailed in supplementary Figure S3, and Table S3.

134

135 RAD sequencing

136

137 DNA has been extracted with the Macherey-Nagel NucleoSpin DNA RapidLyse kit. RAD
138 library preparation (with the PstI restriction enzyme) and sequencing (Illumina NovaSeq600
139 with 150 nucleotides paired-end sequencing) have been performed at the MGX platform
140 (CNRS). The MGX platform performed control quality, demultiplexing and removal of PCR
141 duplicates with unique molecule identifiers. Potential contaminants have been removed with
142 kraken2 (Wood et al., 2019; Lu et al., 2022). RAD loci have been assembled with ipyrad
143 (Eaton and Overcast, 2020). We tested four assembly strategies to test the robustness of
144 the results: a de novo assembly, with a clustering threshold of 0.85, and assembly on a
145 reference genome, with each of the three available genomes: for *E. cavolini*, *E. singularis*,
146 and *E. verrucosa* (Ledoux et al., in prep).

147 From these datasets, we built four datasets focused on the differentiation between
148 *E. cavolini* and *E. singularis*: we excluded *E. verrucosa* samples and we retained the first
percent of the loci with the highest FST between *E. cavolini* and *E. singularis*. These last datasets will be labelled as "1%" (see characteristics of the different datasets are summarised in Table S4).

Transcriptome sequencing and SNPs calling

Total RNA has been extracted as in Haguenauer et al. (2013). RNAs were sent to the LIGAN genomic platform for sequencing (Lille, France) on four flow cells of Illumina NextSeq 500 (2 x 75 bp). The transcriptomes have been assembled with the *de novo* RNA-Seq Assembly Pipeline (DRAP; Cabau et al., 2017) with Oases (Schulz et al., 2012) and default parameters. We performed an individual assembly, and a meta-assembly to be used as reference. The statistics describing the assembled transcriptomes are given in Table S2. We used the BLAT software (Kent, 2002) and the *blat_parser.pl* script to remove potential Symbiodiniaceae sequences in the obtained transcriptomes, with the transcriptome of the type A1 (Baumgarten et al., 2013) as a reference.

We mapped the reads on the meta transcriptome filtered for Symbiodiniaceae sequences with bwa option mem (Li and Durbin, 2009). The obtained sam files were converted in bam format with samtools 1.9 (Li et al., 2009), and sorted with Picard tools ('Picard Toolkit', 2019). The SNPs calling has been performed with reads2snp 2.0 with default parameters (Tsagkogeorga et al., 2012; Gayral et al., 2013). The obtained dataset, including variable and non variable sites, will thereafter be referred as the "all sites" dataset. We performed separate SNP calls with reads2snp for pairwise comparisons among species and without the potential hybrid samples. These three datasets have been used for demographic inferences, and will be referred as "all-CS" for the *E. cavolini* / *E. singularis* comparison, "all-CV" for the *E. cavolini* / *E. verrucosa* comparison, and "all-SV" for the *E. singularis* / *E. verrucosa* comparison.

For an analysis of genetic differentiation, we filtered the "all sites" vcf file with vcftools (Danecek et al., 2011). We retained biallelic sites, without missing data, and separated by at least 1 kb: this is the "polymorphic sites" dataset. As for RAD sequencing, we built a dataset focused on the differentiation between *E. cavolini* and *E. singularis*, retaining the 1% loci with the highest differentiation between *E. cavolini* and *E. singularis* (Table S4).

Presence of Symbiodiniaceae

We analysed the presence of Symbiodiniaceae in *Eunicella gorgonians* with transcriptome data. First, we counted the number of reads corresponding to the Symbiodiniaceae transcriptome type A1 with Salmon (Patro et al., 2017). Second, we used the percentage of
assembled sequences (contigs) in the *Eunicella* transcriptomes corresponding to
Symbiodiniaceae following the BLAT analysis. We used a Kruskal-Wallis test in R to test for
differences among the four groups of samples (the three *Eunicella* species and the potential
hybrids) for each metric. Additionally, we performed a BLAST analysis with the LSU, ITS and
psbA sequences of *Philozoon* (LaJeunesse et al., 2022) to try to identify the
Symbiodiniaceae genera present in the different samples.

As our data pointed to the unexpected presence of Symbiodiniaceae in *E. cavolini* (see
Results), we further explored this topic with preliminary data from another experiment
dedicated to studying the microbiome of *E. cavolini* and *E. singularis*. This pilot study
involved an analysis of microeukaryotic communities through 18S rDNA metabarcoding on
two colonies of *E. cavolini*, and one *E. singularis* (Supplementary File S2).

197

## Genetic differentiation and analysis of hybrids

199

With RAD sequencing data, we performed the analysis of genetic diversity with the four
datasets including all loci. With transcriptomes, we performed the same analyses with the
"polymorphic sites" dataset. We used the LEA R package to estimate ancestry coefficients
(Frichot et al., 2014; Frichot and François, 2015). We tested K values from 1 to 10, with 10
replicates for each K. To analyse the genetic differences among individuals, we performed a
Principal Component Analysis (PCA) with the R package adegenet (Jombart, 2008). The
pairwise FST (Weir and Cockerham, 1984) estimated among species were computed with the
R package Genepop (Rousset, 2008; Rousset et al., 2020), after conversion of the vcf file
with PGDSpider (Lischer & Excoffier, 2012). The distribution of FST among loci was obtained
with vcftools.

The hybrid status (e.g. first generation hybrids) of morphologically intermediate individuals
was analysed with the NewHybrids software (Anderson and Thompson, 2002). We used the
genepopedit R package to prepare the input file from genepop format (Stanley et al., 2017).
Following the results of the LEA and PCA analyses, we compared *E. cavolini*, *E. singularis*
and potential hybrids. The NewHybrids analysis has difficulties to converge with a high
number of loci compared to the number of individuals
([https://github.com/erigande/newhybrids/issues/5](https://github.com/erigande/newhybrids/issues/5)). We therefore used the different "1%
SNP" datasets of RAD sequencing and transcriptome datasets (i.e. the most differentiated
loci) for the NewHybrids analysis. As a prior, we used individuals with the lowest levels of
admixture in LEA as potential parental individuals. For the RAD datasets, this corresponded
to ten individuals of each species as priors. For transcriptome sequencing this corresponded
to three individuals for *E. cavolini*, and six individuals for *E. singularis*. Each NewHybrids
analysis was repeated ten times to test the robustness of the results.
## Scenarios of speciation

We tested scenarios of speciation with the Demographic Inferences with Linked Selection (DILS) pipeline (Csilléry, et al., 2012; Pudlo et al., 2016; Fraïsse et al., 2021) on transcriptome data only. Note that with the high number of loci recovered with transcriptomes, the numbers of specimens used here are adequate for robust inferences (Roux et al., 2016). The DILS pipeline allows the analysis of two species scenarios only: we therefore performed separate analyses for the three two-species comparisons, with the “all-CS”, “all-CV”, and “all-SV” pairwise datasets. We did not include the potential hybrids in the analysis, which would have required the consideration of a separate population. The tested scenarios are presented in Figure S4 (see Fraïsse et al., 2021 for details). Briefly, DILS allows testing scenario with current migration (i.e. gene flow), such as isolation / migration or secondary contact, versus scenarios of current isolation (no gene flow), such as complete or ancestral migration (gene flow among ancestral populations).

We used the same priors for all analyses, with different numbers of sequences per gene and per sample according to the dataset (Table S5). For all pairwise comparisons, we performed two DILS analyses: one with constant population sizes, and one with variable population sizes.

## Results

### Mitochondrial MutS

The mitochondrial MutS sequences available in GenBank confirmed the proximity of the three *Eunicella* species analysed here: all sequences were identical for these three species, as well as for three other sequences deposited in GenBank as unidentified *Eunicella* (Figure S3). The closest species to this group was *Eunicella racemosa*. All other *Eunicella* MutS sequences (*E. tricoronata* and *E. albicans*) grouped separately with *Complexum monodi*, but with low bootstrap support.

### Presence of Symbiodiniaceae

The transcriptomes showed low numbers of reads counts aligning on the Symbiodiniaceae transcriptome (1868 to 58406 reads; Table S6). The proportion of contigs corresponding to Symbiodiniaceae with BLAT was also very low (between 0.00276 and 0.03686; Table S6). Significant differences were observed among species in both cases (Kruskal-Wallis test, p =
0.047 for reads counts, and p = 0.002 for the proportions of contigs). The pairwise Wilcoxon-Test showed significant differences only for the comparisons of proportions of contigs involving *E. singularis*, which was higher than in other species (Table S7; Figure 2). The mean values of reads counts and contigs for the Symbiodiniaceae in the hybrids were lower than in *E. singularis* and *E. cavolini* but higher than in *E. verrucosa*, although pairwise tests were not significant.

The BLAST analysis with the LSU, ITS and psbA sequences of *Philozoon* only retrieved corresponding sequences in the transcriptomes of *E. singularis*. Regarding the pilot study of 18S rDNA metabarcoding, a diversity of 92 Operational Taxonomic Units (OTUs) corresponding to Symbiodiniaceae in the Silva database was observed in *E. singularis*, with a single OTU largely dominant in abundance (Supplementary file S2). The same OTU was also observed in *E. cavolini* with a low abundance of reads, but still representing 99% of all 12 to 13 Symbiodiniaceae OTUs detected in the two analysed colonies. A BLAST search in GenBank identified a subset of Symbiodiniaceae sequences related to this OTU. Phylogenetic inference based on these data indicated that this OTU was related to clade A of the Symbiodiniaceae.

## Genetic differentiation and analysis of hybrids

With RAD sequencing we obtained between 12 952 and 29 061 SNPs for the assembly on *E. cavolini* and *E. verrucosa* genomes respectively (Table S4). The $F_{ST}$ estimates from RAD sequencing were highest for the comparisons between *E. verrucosa* and all other samples ($F_{ST}$ between 0.51 and 0.66 depending on dataset; Table S8). The $F_{ST}$ between *E. cavolini* and *E. singularis* was lower ($F_{ST}$ between 0.29 and 0.38), and the lowest $F_{ST}$ values were observed for hybrids compared to these two species ($F_{ST}$ between 0.09 and 0.13). The cross-entropy analysis using LEA with RAD sequencing showed a minimum at $K = 3$ for the four datasets (results not shown). The barplots of coancestry coefficients were very similar for the four datasets, with a separation of the three species, and an admixture between *E. cavolini* and *E. singularis* for the morphologically intermediate individuals (Figure S5). The PCAs on RAD sequencing were very similar for all datasets, with a separation between *E. verrucosa* and all other samples on the first axis (Figure S6). The second axis separated *E. cavolini* and *E. singularis*, with the potential hybrids in intermediate position between them. Projections on axes 3 and 4 resulted mainly in the separations of *E. verrucosa* samples from each other.

With transcriptomes, we obtained 31 369 SNPs for the "polymorphic sites" dataset. With this dataset, the highest $F_{ST}$ values were observed for the comparisons between *E. verrucosa* and all other samples ($F_{ST} > 0.43$; Table S9). The $F_{ST}$ between *E. cavolini* and *E. singularis* was much lower (0.21), and the lowest $F_{ST}$ values were observed for hybrids compared to
these two species ($F_{ST}$ around 0.07 in both cases). These differences corresponded to different distributions of $F_{ST}$ over SNPs (Figure S7). For the 1% SNPs with the highest $F_{ST}$ estimates, 52 SNPs were shared by both comparisons involving *E. cavolini* (i.e. *E. cavolini* vs *E. singularis* and *E. cavolini* vs *E. verrucosa*), 116 top 1% SNPs were shared by both comparisons involving *E. singularis*, and 1042 top 1% SNPs were shared by both comparisons involving *E. verrucosa*.

The cross-entropy analysis using LEA with transcriptomes indicated a best clustering solution corresponding to K = 2 or K = 3 clusters (Figure S8). At K = 2, the first distinction was observed between *E. verrucosa* and all other samples (Figure 3). The K = 3 analysis further separated *E. cavolini* and *E. singularis*, with morphologically intermediate individuals admixed between these two species. Conversely the individuals representative of *E. cavolini* and *E. singularis* presented low levels of admixture, apart from the *E. cavolini* of the site in Algeria (code ANB), and, at a small level, two *E. singularis* individuals from Banyuls (BAN). At K = 4, the two *E. cavolini* individuals from Algeria separated from other *E. cavolini* from the northern part of the Mediterranean.

As with RAD sequencing, the PCA on transcriptome SNPs separated *E. verrucosa* from other samples on the first axis (Figure 4). The second axis separated *E. cavolini* and *E. singularis*, with the potential hybrids in intermediate position between them. The third axis separated the *E. cavolini* samples from Algeria (ANB site) from all other samples (Figure S6).

The NewHybrids analysis on RAD sequencing indicated that all morphologically intermediate individuals, except one, appeared as hybrids: first generation (F1), second generation (F2), or backcrosses with *E. singularis* or *E. cavolini* (Table 1). These samples also appeared admixed on the basis of LEA (Figure S5). One individual identified as a potential hybrid in *situ*, was inferred as a parental *E. singularis*. For four individuals, the hybrid status varied according to the dataset: F2 or backcross with *E. cavolini* in two cases, F1 or F2 in two cases. Potential parental individuals not included in the priors were well inferred as parental with NewHybrids. The NewHybrids analysis with transcriptomes indicated that the morphologically intermediate individuals were hybrids with a probability of one in all ten iterations of the analysis. One individual was a F1 hybrid, another one was F2 hybrid, and the two other ones corresponded to backcrosses with *E. singularis* (Figure 3; Table 1). In the same analysis, the *E. cavolini* and *E. singularis* individuals not included as priors for parental species (see Figure 3 for the individuals used as priors), were indeed inferred as parental with a probability of one, including the *E. cavolini* individual from Algeria (ANB).

## Scenarios of speciation
The average pairwise net divergence estimated from DILS was 0.0018 between *E. cavolini* and *E. singularis*, and around 0.007 for the two comparisons with *E. verrucosa* (Table S9, https://zenodo.org/records/12532817/files/results_DILS_suppl_file.ods?download=1). The DILS analysis indicated the existence of current gene flow between *E. cavolini* and *E. singularis* with high probability, both with constant and variable population sizes (p = 0.87 and 0.88 respectively; Table 2). This possibility of gene flow corresponded to a scenario of secondary contact. Conversely, a model of current isolation was inferred for the comparisons between *E. verrucosa* and each of the two other species, with a probability p ≥ 0.87: in these two cases, the inferred scenario included a period of ancestral migration, though with moderate support (p between 0.61 and 0.69). A genomic heterogeneity in effective size (i.e. variations among loci) was inferred with strong support (p ≥ 0.99) for all analyses. In the case of current gene flow (between *E. cavolini* and *E. singularis*), a genomic heterogeneity in migration rates was inferred (p ≥ 0.82). We repeated the DILS analysis without including the two divergent samples of *E. cavolini* from Algeria: this led to similar results, with inference of secondary contact for the comparison with *E. singularis*, and ancestral migration for the comparison with *E. verrucosa* (results not shown). For parameters inferences, we used the complete datasets, with all *E. cavolini* samples. The inferred parameters for the different scenarios are presented in Supplementary Table S9. We will first present the results obtained for the constant population sizes models. The divergence time between *E. cavolini* and *E. singularis* (median 403 273 generations) was much lower than between *E. cavolini* and *E. verrucosa* (median 1 054 488 generations), and between *E. singularis* and *E. verrucosa* (median 899 098 generations). For the comparison between *E. cavolini* and *E. singularis*, the time of secondary contact was estimated after around 85% of time spent in isolation since divergence. Following secondary contact, the gene flow was similar in both directions for these two species. The duration of ancestral migration roughly corresponded to 6% and 8% of the total time since divergence for the comparison between *E. cavolini* and *E. verrucosa*, and for the comparison between *E. singularis* and *E. verrucosa*, respectively. For these last two cases, the gene flow (forward in time) during ancestral migration was higher towards *E. verrucosa* than in the opposite direction. The estimated effective sizes were of similar order for *E. cavolini* and *E. verrucosa*. Similar results were obtained for the models including variations in effective size, except for the estimate of current gene flow between *E. cavolini* and *E. singularis*: with variable population size, gene flow from *E. singularis* to *E. cavolini* was higher than in the opposite direction.

## Discussion
The three named *Eunicella* species studied here have been previously described with differences in colony morphology, sclerites shape, and in the presence of photosynthetic Symbiodiniaceae (Carpine and Grasshoff, 1975). Our results demonstrate a continuum between *E. cavolini* and *E. singularis*, with morphologically intermediate individuals, on-going gene flow, and hybrids characterised by a reduced frequency of Symbiodiniaceae compared to *E. singularis*. On the other hand, *E. verrucosa* appears genetically isolated from these two species. We will discuss here the differences observed among markers, the outcome of hybridization, the speciation scenarios, and what can be learnt on the evolution of symbiosis.

## Discordances between molecular markers

As previously observed (Aurelle et al., 2017), mitochondrial DNA did not allow to discriminate the three species due to the usually slow evolution of mitochondrial DNA in octocorals (McFadden et al., 2011; Muthye et al., 2022). The use of transcriptome sequences first confirmed the closer proximity between *E. cavolini* and *E. singularis* than with *E. verrucosa*. This had been previously suggested with two intron sequences, but with incomplete lineage sorting (Aurelle et al., 2017). The Mediterranean *Eunicella* then add a new example of the lack of power of mitochondrial DNA to discriminate genetically differentiated octocoral species, as demonstrated in other genera (Erickson et al., 2021; Pante et al., 2015a). The slow rate of evolution of mitochondrial DNA in octocorals has been linked to the presence of the mitochondrial locus MutS, an homolog of a bacterial gene involved in DNA repair. However, there are contradictory examples showing that the presence of this locus is not the only factor explaining the slow evolution of mitochondrial DNA in octocorals (Muthye et al., 2022). More generally, as hybridization can lead to the sharing of mitochondrial DNA among species, the use of multiple independent nuclear loci is required for species discrimination in such cases.

## Incomplete reproductive isolation among two named species

Inferences of genetic ancestry and hybrid status confirmed that morphologically intermediate individuals are indeed hybrids between *E. singularis* and *E. cavolini*, with the identification of F1, F2 and backcrosses with both parental lineages: first generation hybrids can then be fertile. The fact that gene flow indeed goes further than the hybrid levels is confirmed by the DILS analysis, which did not include hybrid individuals. Reproductive isolation is therefore at least partial between these lineages. The ease to find hybrids in the area studied here, as well as similar observations in other sites (S. Sartoretto, pers. com.) indicate that hybridization is not rare on an evolutionary scale. Similarly, transcriptome sequencing has
408 led to infer hybridization among *Plexaura* species on the basis of a small number of samples
409 (Pelosi *et al.*, 2022).

410 The alternation of populations with and without hybrids would point to a mosaic hybrid zone
411 (Bierne *et al.*, 2003), where hybrids could form in different areas and from different parental
412 populations. As, or because, hybridization between *E. cavolini* and *E. singularis* had not
413 been reported before, the presence of hybrids has probably been overlooked up to now. This
414 may be the consequence of previously focusing on colonies with "typical" morphologies. The
415 frequency of hybridization therefore remains to be studied.

416 Our results allow discussing the evolution of genomic divergence among these species. The
417 persistence of genomic differentiation between these lineages in sympatry, despite current
418 gene flow, indicates that intrinsic (i.e. genomic incompatibilities) or extrinsic (e.g. ecology)
419 factors can maintain partial isolation. Difference or overlap in the timing of reproduction
420 should also be considered in contributing to pre-mating isolation (Pelosi et al., 2022). A
421 better characterization of the ecological range of parental and hybrid populations would be
422 useful to test if local adaptation is involved in their distribution. Intrinsic factors such as
423 genetic incompatibilities, potentially coupled with differences in adaptation to local
424 environments, can be present as well (Bierne et al., 2011). A genome wide analysis of
425 differentiation is required to investigate whether divergence between *E. cavolini* and
426 *E. singularis*, is homogeneous along the genome (as suggested by the DILS analysis which
427 inferred a homogeneity of gene flow), or whether genomic islands of differentiation exist
428 (Peñalba et al., 2024). We could then better understand to what stage of divergence the
429 *E. cavolini* / *E. singularis* split corresponds: from intra-specific polymorphism to species
430 separated by semipermeable barriers to gene flow.

431 One interesting question in this context is whether changes in selection regimes induced by
432 human activities can change the outcome of hybridization (Ålund et al., 2023). For example,
433 Mediterranean octocorals are impacted by mortality events linked with climate change (Sini
434 et al., 2015; Estaque et al., 2023), and the impact of these events could be different for
435 hybrids and parental individuals. In scleractinian corals, interspecific hybridization has been
436 reported to enhance the survival under elevated temperature conditions (Chan et al., 2018).
437 Regarding *E. verrucosa*, the more ancient divergence corresponded to much more loci with
438 high FST. Among the list of the most highly differentiated loci, more overlap was also
439 observed for the two comparisons involving *E. verrucosa* than for the other pairwise
440 comparisons: this may indicate that few genomic areas of potential incompatibilities with
441 *E. verrucosa* are involved in the divergence between *E. cavolini* and *E. singularis*.

442

443 Scenarios of speciations

444
The scenarios of speciations inferred with DILS supported the current isolation (no gene flow) of *E. verrucosa* with the two other species with high posterior probability. Conversely current gene flow was strongly supported versus isolation between *E. cavolini* and *E. singularis*. The posterior probabilities for ancestral migration (for *E. verrucosa* versus the two other species), and secondary contact (*E. cavolini* and *E. singularis*), were lower than for inferences on current gene flow. These scenarios were indeed the best ones among those tested here but they might not provide the best possible representation of the evolutionary history. Other models of evolution could be tested for better inferences, for example by including the three species and hybrids, or gene flow from unsampled taxa (Tricou *et al.*, 2022). The current isolation of *E. verrucosa* from *E. cavolini* is also at odds with previous results which showed the possibility of current gene flow between these two species despite an important divergence (Roux *et al.*, 2016). It will be useful to explore the reasons for the discrepancy between this last study and the present one, which are both based on transcriptome datasets but obtained from different samples and sequencing platforms.

*Eunicella verrucosa* is currently widely distributed in the North Eastern Atlantic Ocean, and less frequent in the Mediterranean Sea, whereas both other species are only present in the Mediterranean Sea. The Atlantic / Mediterranean Sea transition does not seem to act as a phylogeographic barrier for *E. verrucosa* (Macleod *et al.*, 2024). We can propose a scenario where the split between *E. verrucosa* and both other species occurred in allopatry between the Atlantic Ocean and the Mediterranean Sea, followed by the colonization of the Mediterranean Sea by *E. verrucosa*. The generation time remains unknown for the *Eunicella* species, and previous studies have shown important variation in the age at first reproduction in gorgonians, from 2 to 13 years (see references in Munro, 2004). If we use a generation time of two years for *Eunicella* species, with a median estimate of divergence time around 900 000 generations for *E. verrucosa* / *E. singularis* and 1 000 000 for *E. verrucosa* / *E. cavolini*, and based on a mutation rate set at 3.10⁻⁹, this would indicate a divergence at least around 2 000 000 years (2 Ma). The divergence time between *E. cavolini* and *E. singularis* would be 2.5 times more recent, around 800 000 years, with a median time of secondary contact around 60 000 generations, corresponding to 15% of the time spent since divergence. It is difficult to infer past distributions of *E. singularis* and *E. cavolini*, but one can note that even if they are currently found in sympatry in different areas, their range do not completely overlap. For example *E. cavolini* is nearly absent at the West of the Rhone estuary on the French coast, whereas *E. singularis* is present there. The ecological range of *E. singularis* and *E. cavolini* is also not completely overlapping, as *E. cavolini* can be observed deeper than *E. singularis* (Gori *et al.*, 2012; Carugati *et al.*, 2022). Therefore one can envision a historical separation of these two species either geographically or ecologically, followed by a secondary contact where gene flow took place. In any case,
482

additional information on generation time, mutation rate and past demographic fluctuations are required to be more precise on the history of these species.

483

484

485

## Evolution of symbiosis

486

487

As previously discussed, we clearly demonstrated here the possibility of gene flow between symbiotic (hosting Symbiodiniaceae) and non-symbiotic octocorals. Symbiodiniaceae could nevertheless be involved in genetic incompatibilities with the genome of some cnidarian hosts, but this would require additional analysis of symbiotic status in hybrids. The methods used here did not aim at a precise quantification of Symbiodiniaceae, and one can note the low levels of sequences corresponding to these symbionts, even in *E. singularis*, which may be due to difficulties in extracting the RNA of the symbionts (but see Guzman *et al.*, 2018; Rivera-García *et al.*, 2019). Despite these limits we observed, as expected, a higher Symbiodiniaceae concentration in *E. singularis* than in *E. cavolini* and *E. verrucosa*. Interestingly, the hybrids showed a lower frequency of Symbiodiniaceae than *E. singularis*, and possibly than *E. cavolini*, though this last result remains to be confirmed. In *E. singularis*, the transmission of Symbiodiniaceae seems to occur both vertically, through ovules, and horizontally, from the environment (Forcioli *et al.*, 2011). Both transmission modes did not restore the levels of Symbiodiniaceae in the hybrids to those of *E. singularis*.

488

489

490

491

492

493

494

495

496

497

498

499

500

501

502

503

504

505

506

507

508

509

510

This suggests a breakdown of or a failure to establish symbiosis for hybrid genotypes, which may impact the fitness of hybrids and consequently the possibility of introgression. The *aphyta* type of *E. singularis* observed in deep conditions indicates a plasticity of symbiotic status apart from hybridization. Nevertheless, here the hybrids were sampled in shallow conditions (10-20 m depth) which underlines the role of hybridization in reducing the extent of symbiosis. More precise estimates of Symbiodiniaceae abundance, and of physiological parameters such as photosynthetic and respiration rates (Ezzat *et al.*, 2013). would help understanding the role of symbionts in hybrids fitness. It would also be interesting to study if the Symbiodiniaceae of the different samples belong to the same population (Pelosi *et al.*, 2022).

511

512

513

514

515

516

517

518

Our results also question the evolution and significance of octocoral / Symbiodiniaceae symbiosis. In scleractinians, the transition between symbiotic and non-symbiotic states happened repeatedly, but mostly in the direction of the acquisition of symbiosis, with very low rates of reversal (Campoy *et al.*, 2020). This could indicate that investing in such mutualistic interactions for the cnidarian would lead to increasingly relying on autotrophy for energetic supply, making reversal to heterotrophy difficult. In octocorals, an evolutionary versatility in symbiotic state seems possible, as in various families and genera, both symbiotic and non-symbiotic species are present (Van Oppen *et al.*, 2005). In the
519 Mediterranean Sea, all octocoral species are non-symbiotic, except for *E. singularis* (but see
520 Bonacolta et al. 2024). The most parsimonious scenario here would be an acquisition of
521 symbiosis in *E. singularis* during or following its divergence from *E. cavolini*. The symbiotic
522 status of *E. singularis* nevertheless could be facultative as previously mentioned for the
523 *aphyta* type (Gori *et al.*, 2012). Additionally, experimental physiological studies have
524 demonstrated the nutritional plasticity of *E. singularis* which is able to use either heterotrophy
525 or autotrophy for its metabolism (Ezzat *et al.*, 2013). Nevertheless, in natural conditions,
526 autotrophy seems to provide an important contribution to the metabolism of *E. singularis*,
527 and the collapse of photosynthetic capacities in too warm conditions can contribute to
528 mortality events in this species (Coma *et al.*, 2015).

529 The question of symbiosis could be reversed as well: why are Symbiodiniaceae not more
530 abundant in *E. cavolini*? This species can be observed in shallow conditions (less than 10 m
531 depth) where there is enough light for photosynthesis, and in syntopy with *E. singularis*. The
532 availability of preys or particulate organic matter may provide enough energy to *E. cavolini* in
533 its habitat, but this species may have never engaged in mutualistic interaction with
534 Symbiodiniaceae. Interestingly we observed a low rate of sequences related to
535 Symbiodiniaceae in the transcriptomes of *E. cavolini* (and even lower, but not null in
536 *E. verrucosa*). This could either correspond to a signal from free living Symbiodiniaceae, or
537 to rare, transient, associations with the cnidarian. In addition, a Symbiodiniaceae OTU that is
538 common to *E. singularis* and *E. cavolini* was identified among the microeukaryotes
539 associated with the two species: this OTU is related to strains observed in symbiosis with
540 *E. singularis* and other Mediterranean cnidarians. Molecular markers also allowed to
541 evidence the presence of Symbiodiniaceae in species previously supposed to be asymbiotic,
542 as in the Mediterranean octocoral *Paramuricea clavata*, and in several Hawaiian
543 antipatharian species (Wagner *et al.*, 2011; Bonacolta *et al.*, 2024). These results, and our
544 observations in *Eunicella* species, obviously underline the dynamic nature of interactions
545 between Symbiodiniaceae and cnidarians: the establishment of symbiosis may be preceded
546 by more or less stable, and more or less mutualistic interactions. The development of
547 effective symbiosis, with stable relationships, and higher abundance of symbiont, would
548 require specific adaptation from both partners. We can see here that even if on a macro-
549 evolutionary scale the acquisition of symbiosis is much more frequent than its loss, on a
550 micro-evolutionary scale the gene flow between the *Eunicella* species considered here has
551 not led to the full development of symbiosis in *E. cavolini*.

552

Conclusions and perspectives
553

554
We demonstrated the lack of genetic isolation between octocorals with contrasted levels of mutualistic interaction with Symbiodiniaceae. Understanding the evolution and adaptation of these species in heterogeneous environments should then consider the possible impact of introgression. We also show that symbiosis is more flexible that previously envisioned in octocorals. For these species it will be useful to estimate the frequency and spatial extent of hybrid zones: does it correlate with particular environments with a coupling between endogenous and exogenous barriers to gene flow (Bierne et al., 2011)? Characterizing the genomic landscape of introgression would help to look for the effects of introgression on adaptation or symbiosis for example. Indeed, even low levels of interspecific gene flow can have important consequences on the evolution of species (Arnold et al., 1999). Finally, various cases of hybridization have been demonstrated in symbiotic anthozoans (e.g. Combosch and Vollmer, 2015; Pelosi et al., 2022): it would then be interesting to study the dynamics of symbiosis in these cases, especially when different Symbiodiniaceae strains are involved.

569

570

## Acknowledgements:

572

We thank the ECCOREV Research Federation (FR 3098) for the financial support of part of this study (https://www.eccorev.fr/). The project leading to this publication has received funding from European FEDER Fund under project 1166-39417. The project leading to this publication has received funding from Excellence Initiative of Aix-Marseille University - A*MIDEX, a French "Investissements d'Avenir" programme. The authors thank the UMR 8199 LIGAN-PM Genomics platform (Lille, France, especially Véronique Dhennin) which belongs to the 'Federation de Recherche' 3508 Labex EGID (European Genomics Institute for Diabetes; ANR-10-LABX-46) and was supported by the ANR Equipex 2010 session (ANR-10-EQPX-07-01; 'LIGAN-PM'). The LIGAN-PM Genomics platform (Lille, France) is also supported by the FEDER and the Region Nord-Pas-de-Calais-Picardie. JBL was supported by the strategic funding UIDB/04423/2020, UIDP/04423/2020 and 2021.00855.CEECIND through national funds provided by FCT -Fundação para a Ciência e a Tecnologia. Reference genomes were obtained with the support from EASI-genomics funded from the European Union's Horizon 2020 research and innovation programme under grant agreement No 824110 to J.-B.L. Camille Roux, Jonathan Romiguier and Christelle Fraïsse were of a great help for the analysis scenarios of speciation. We thank the diving service of INSU/OSU Pytheas for fieldwork, and the Calanques National Park for sampling authorisations. We acknowledge the staff of the "Cluster de calcul intensif HPC" Platform of the OSU Institut Pythéas (Aix-Marseille Université, INSU-CNRS) for providing the computing
facilities. We are grateful to the Genotoul bioinformatics platform Toulouse Occitanie (Bioinfo Genotoul, https://doi.org/10.15454/1.5572369328961167E12) for providing help, computing and storage resources. We thank Christophe Klopp and Marie-Stéphane Trotard for their help. We acknowledge the use of the computing cluster of MNHN (Plateforme de Calcul Intensif et Algorithmique PCIA, Muséum National d'Histoire Naturelle, Centre national de la recherche scientifique, UAR 2700 2AD, CP 26, 57 rue Cuvier, F-75231 Paris Cedex 05, France). Part of the bioinformatics analyses have been performed on the Core Cluster of the Institut Français de Bioinformatique (IFB) (ANR-11-INBS-0013). MGX acknowledges financial support from France Génomique National infrastructure, funded as part of "Investissement d'Avenir" program managed by Agence Nationale pour la Recherche (contract ANR-10-INBS-09). Part of this work has been performed during a CNRS detachment position of D. Aurelle at the ISYEB laboratory.

## Data availability

The transcriptome raw sequences are available in Genbank under BioProject ID PRJNA1037721. The RAD raw sequences are available in Genbank under BioProject ID PRJNA1122331.

The scripts used in this study, the vcf files from RAD sequencing, and the detailed results of the DILS analysis are available at https://doi.org/10.5281/zenodo.14007931

## Conflict of interest disclosure

The authors declare that they have no conflict of interest in relation to the content of the article
Ålund M, Cenzer M, Bierne N, et al. Anthropogenic Change and the Process of Speciation.
*Cold Spring Harbor Perspectives in Biology* 2023;15:a041455.

https://cshperspectives.cshlp.org/content/15/12/a041455.short

Anderson E, Thompson EA. A model-based method for identifying species hybrids using multilocus genetic data. *Genetics* 2002;**160**:1217–1229.

https://academic.oup.com/genetics/article-abstract/160/3/1217/6052497

Arnold ML, Bulger MR, Burke JM, et al. 1999. Natural hybridization: how low can you go and still be important? *Ecology* 1999;**80**:371–381.

https://esajournals.onlinelibrary.wiley.com/doi/abs/10.1890/0012-9658(1999)080[0371:NHHLCY]2.0.CO;2

Aurelle D, Pivotto ID, Malfant M, et al. Fuzzy species limits in Mediterranean gorgonians (Cnidaria, Octocorallia): inferences on speciation processes. *Zoologica Scripta* 2017;**46**:767–778. https://onlinelibrary.wiley.com/doi/abs/10.1111/zsc.12245

Baird NA, Etter PD, Atwood TS, et al. Rapid SNP Discovery and Genetic Mapping Using Sequenced RAD Markers. *PLoS ONE* 2008;**3**:e3376.

https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0003376

Baumgarten S, Bayer T, Aranda M, et al. Integrating microRNA and mRNA expression profiling in Symbiodinium microadriaticum, a dinoflagellate symbiont of reef-building corals.

BMC Genomics 2013;14:704. https://link.springer.com/article/10.1186/1471-2164-14-704

Berkelmans R and van Oppen MJH. The role of zooxanthellae in the thermal tolerance of corals: a ‘nugget of hope’ for coral reefs in an era of climate change. *Proceedings of the Royal Society B: Biological Sciences* 2006;**273**:2305–2312.

https://royalsocietypublishing.org/doi/abs/10.1098/rspb.2006.3567

Bierne N, Borsa P, Daguin C, et al. Introgression patterns in the mosaic hybrid zone between Mytilus edulis and M. galloprovincialis. *Molecular Ecology* 2003;**12**:447–462.

https://onlinelibrary.wiley.com/doi/abs/10.1046/j.1365-294X.2003.01730.x

Bierne N, Welch J, Loire E, et al. The coupling hypothesis: why genome scans may fail to map local adaptation genes. *Molecular Ecology* 2011;**20**:2044–2072.

https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1365-294X.2011.05080.x

Bonacolta AM, Miravall J, Gómez-Gras D, et al. Differential apicomplexan presence predicts thermal stress mortality in the Mediterranean coral Paramuricea clavata. *Environmental Microbiology* 2024;**26**:e16548.

https://enviromicro-journals.onlinelibrary.wiley.com/doi/abs/10.1111/1462-2920.16548

Bordenstein S. Symbiosis And The Origin Of Species. In: *Insect Symbiosis*. CRC Press, 2003;283–303.

Brener-Raffalli K, Vidal-Dupiol J, Adjeroud M, et al. Gene expression plasticity and frontloading promote thermotolerance in *Pocillopora* corals. *Peer Community Journal* 2022;**2**. https://peercommunityjournal.org/articles/10.24072/pcjournal.79/

Brucker RM and Bordenstein SR. Speciation by symbiosis. *Trends in Ecology & Evolution* 2012;**27**:443. https://www.cell.com/trends/ecology-evolution/fulltext/S0169-5347(12)00076-6

Cabau C, Escudié F, Djari A, et al. Compacting and correcting Trinity and Oases RNA-Seq de novo assemblies. *PeerJ* 2017;**5**:e2988. https://peerj.com/articles/2988/
Cahill AE, Meglécz E, Chenuil A. Scientific history, biogeography, and biological traits predict presence of cryptic or overlooked species. *Biological Reviews* 2024;99:546-561.

https://onlinelibrary.wiley.com/doi/abs/10.1111/brv.13034

Cairns SD. Deep-water corals: an overview with special reference to diversity and distribution of deep-water scleractinian corals. *Bulletin of marine Science* 2007;81:311–322.

https://www.ingentaconnect.com/content/umrsmas/bullmar/2007/00000081/00000003/art00002

Calderón I, Garrabou J, Aurelle D. Evaluation of the utility of COI and ITS markers as tools for population genetic studies of temperate gorgonians. *Journal of Experimental Marine Biology and Ecology* 2006;336:184–197.

https://www.sciencedirect.com/science/article/pii/S0022098106002498

Campoy AN, Addamo AM, Machordom A, et al. The Origin and Correlated Evolution of Symbiosis and Coloniality in Scleractinian Corals. *Frontiers in Marine Science* 2020;7.

https://www.frontiersin.org/articles/10.3389/fmars.2020.00461/full

Carpine C, Grasshoff M. Les gorgonaires de la Méditerranée. *Bulletin de l'Institut Océanographique de Monaco* 1975;71(140).

Carugati L, Moccia D, Bramanti L, et al. Deep-Dwelling Populations of Mediterranean Corallium rubrum and Eunicella cavolini: Distribution, Demography, and Co-Occurrence.

*Biology* 2022;11. https://www.mdpi.com/2079-7737/11/2/333

Casado-Amezúa P, Terrón-Sigler A, Pinzón JH, et al. General ecological aspects of Anthozoan-Symbiodinium interactions in the Mediterranean Sea. In: Goffredo S, Dubinsky Z, (ed). The cnidaria, past, present and future: the world of medusa and her sisters. Springer, 2016;375–386. https://link.springer.com/chapter/10.1007/978-3-319-31305-4_24

Chan WY, Peplow LM, Menéndez P, Hoffmann AA & Van Oppen MJ. 2018. Interspecific hybridization may provide novel opportunities for coral reef restoration. *Frontiers in Marine Science* 5: 160. https://www.frontiersin.org/articles/10.3389/fmars.2018.00160/full

Chimienti G. 2020. Vulnerable forests of the pink sea fan Eunicella verrucosa in the Mediterranean Sea. *Diversity* 12: 176.

Coelho M, Pearson G, Boavida J, et al. Not out of the Mediterranean: Atlantic populations of the gorgonian Paramuricea clavata are a separate sister species under further lineage diversification. *Ecology and Evolution* 2023;13.

https://onlinelibrary.wiley.com/doi/abs/10.1002/ece3.9740

Coma R, Llorente-Llurba E, Serrano E, et al. Natural heterotrophic feeding by a temperate octocoral with symbiotic zooxanthellae: a contribution to understanding the mechanisms of die-off events. *Coral Reefs* 2015;34:549–560.

https://link.springer.com/article/10.1007/s00338-015-1281-3

Combosch DJ, Vollmer SV. Trans-Pacific RAD-Seq population genomics confirms introgressive hybridization in Eastern Pacific Pocillopora corals. *Molecular phylogenetics and evolution* 2015;88:154–162.

https://www.sciencedirect.com/science/article/pii/S1055790315000858

Csilléry K, François O, Blum MGB. abc: an R package for approximate Bayesian computation (ABC). *Methods in Ecology and Evolution* 2012;3:475–479.

https://besjournals.onlinelibrary.wiley.com/doi/abs/10.1111/j.2041-210X.2011.00179.x

Danecek P, Auton A, Abecasis G, et al. The variant call format and VCFtools. *Bioinformatics* 2011;27:2156–2158. https://academic.oup.com/bioinformatics/article/27/15/2156/402296

De Jode A, Le Moan A, Johannesson K, et al. Ten years of demographic modelling of divergence and speciation in the sea. *Evolutionary Applications* 2023;16:542–559.

https://onlinelibrary.wiley.com/doi/abs/10.1111/eva.13428
De Queiroz K. Species Concepts and Species Delimitation. Systematic Biology
2007;56:879-886. https://academic.oup.com/sysbio/article-abstract/56/6/879/1653163

Eaton DA, Overcast I. ipyrad: Interactive assembly and analysis of RADseq datasets.

*Bioinformatics* 2020; **36**:2592–2594. https://academic.oup.com/bioinformatics/article-abstract/36/8/2592/5697088

Erickson KL, Pentico A, Quattrini AM et al. New approaches to species delimitation and population structure of anthozoans: Two case studies of octocorals using ultraconserved elements and exons. *Molecular Ecology Resources* 2021; **21**:78–92.

https://onlinelibrary.wiley.com/doi/abs/10.1111/1755-0998.13241

Estaque T, Richaume J, Bianchimani O, et al. Marine heatwaves on the rise: One of the strongest ever observed mass mortality event in temperate gorgonians. *Global change biology* 2023; **29**:6159-6162. https://onlinelibrary.wiley.com/doi/10.1111/gcb.16931

Ezzat L, Merle PL, Furla P, et al. The Response of the Mediterranean Gorgonian Eunicella singularis to Thermal Stress Is Independent of Its Nutritional Regime. *PLoS ONE*
2013; **8**:e64370. https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0064370

Faria R, Johannesson K, Stankowski S. Speciation in marine environments: Diving under the surface. *Journal of Evolutionary Biology* 2021; **34**:4–15. https://academic.oup.com/jeb/article-abstract/34/1/4/7326591

Forcioli D, Merle PL, Caligara C, et al. Symbiont diversity is not involved in depth acclimation in the Mediterranean sea whip Eunicella singularis. *Marine Ecology Progress Series* 2011; **439**:57–71. https://www.int-res.com/abstracts/meps/v439/p57-71/

Fourt M, Goujard A. Rapport final de la campagne MEDSEACAN (Têtes des canyons méditerranéens continentaux) novembre 2008–avril 2010. *Partenariat Agence des aires marines protégées–GIS Posidonie*: 1–218. 2012.

http://paleopolis.rediris.es/benthos/TaP/Rapport_Final_MEDSEACAN.pdf

Fraïsse C, Popovic I, Mazoyer C, et al. DILS: Demographic inferences with linked selection by using ABC. *Molecular Ecology Resources* 2021; **21**:2629–2644.

https://onlinelibrary.wiley.com/doi/abs/10.1111/1755-0998.13323

Frichot E, Mathieu F, Trouillon T, et al. Fast and efficient estimation of individual ancestry coefficients. *Genetics* 2014; **196**:973–983.

https://academic.oup.com/genetics/article/196/4/973/5935614

Frichot E, François O. LEA: an R package for landscape and ecological association studies. *Methods in Ecology and Evolution* 2015; **6**:925–929.

https://besjournals.onlinelibrary.wiley.com/doi/abs/10.1111/2041-210X.12382

Furla P, Allemand D, Shick JM, et al. The symbiotic anthozoan: a physiological chimera between alga and animal. *Integrative and Comparative Biology* 2005; **45**:595–604.

https://academic.oup.com/icb/article-abstract/45/4/595/636401

Gagnaire P, Broquet T, Aurelle D, et al. Using neutral, selected, and hitchhiker loci to assess connectivity of marine populations in the genomic era. *Evolutionary Applications* 2015; **8**:769–786. https://onlinelibrary.wiley.com/doi/abs/10.1111/eva.12288

Gayral P, Melo-Ferreira J, Glemin S, et al. Reference-free population genomics from next-generation transcriptome data and the vertebrate–invertebrate gap. *PLoS Genetics* 2013; **9**:e1003457.

https://journals.plos.org/plosgenetics/article?id=10.1371/journal.pgen.1003457

Gori A, Bramanti L, López-González P, et al. Characterization of the zooxanthellate and azooxanthellate morphotypes of the Mediterranean gorgonian Eunicella singularis. *Marine biology* 2012; **159**:1485–1496. https://link.springer.com/article/10.1007/s00227-012-1928-3
Grasshoff, M. Die Flachwasser-Gorgonarien von Europa und Westafrika (Cnidaria, Anthozoa). *Courier Forschungsinstitut Senckenberg* 1992; **149**. Frankfurt a. M.

Guzman C, Shinzato C, Lu TM *et al*. Transcriptome analysis of the reef-building octocoral, *Heliopora coerulea*. *Scientific Reports* 2018; **8**:8397. https://www.nature.com/articles/s41598-018-26718-5

Haguenauer A, Zuberer F, Ledoux JB *et al*. Adaptive abilities of the Mediterranean red coral *Corallium rubrum* in a heterogeneous and changing environment: from population to functional genetics. *Journal of Experimental Marine Biology and Ecology* 2013; **449**:349–357.
Jombart T. adegenet: a R package for the multivariate analysis of genetic markers.
*Bioinformatics* 2008; **24**:1403-1405.

https://www.sciencedirect.com/science/article/pii/S0022098113003493

Kent WJ. BLAT—the BLAST-like alignment tool. *Genome research* 2002; **12**:656–664.

https://genome.cshlp.org/content/12/4/656.short

Krueger-Hadfield S. marmap. https://www.molecularecologist.com/2015/07/03/marmap/. 2015

LaJeunesse TC, Parkinson JE, Gabrielson PW, *et al*. Systematic revision of Symbiodiniaceae highlights the antiquity and diversity of coral endosymbionts. *Current Biology* 2018; **28**:2570–2580. https://www.cell.com/current-biology/fulltext/S0960-9822(18)30907-2

LaJeunesse TC, Wiedenmann J, Casado-Amezúa P, *et al*. Revival of Philozoon Geddes for host-specialized dinoflagellates, ‘zooxanthellae’, in animals from coastal temperate zones of northern and southern hemispheres. *European Journal of Phycology* 2022; **57**:166–180.
https://www.tandfonline.com/doi/abs/10.1080/09670262.2021.1914863

Leroy T, Louvet JM, Lalanne C, *et al*. Adaptive introgression as a driver of local adaptation to climate in European white oaks. *New Phytologist* 2020; **226**:1171–1182.
https://nph.onlinelibrary.wiley.com/doi/abs/10.1111/nph.16095

Lesser MP, Stat M, Gates RD. The endosymbiotic dinoflagellates (Symbiodinium sp.) of corals are parasites and mutualists. *Coral Reefs* 2013; **32**:603–611.
https://link.springer.com/article/10.1007/s00338-013-1051-z

Li H, Handsaker B, Wysoker A, *et al*. The sequence alignment/map format and SAMtools. *Bioinformatics* 2009; **25**:2078–2079. https://academic.oup.com/bioinformatics/article-abstract/25/16/2078/204688

Li H, Durbin R. Fast and accurate short read alignment with Burrows–Wheeler transform.
*Bioinformatics* 2009; **25**:1754–1760.
https://academic.oup.com/bioinformatics/article/25/14/1754/225615

Lischer HE, Excoffier L. PGDSpider: an automated data conversion tool for connecting population genetics and genomics programs. *Bioinformatics* 2012; **28**:298–299.
https://academic.oup.com/bioinformatics/article/28/2/298/198891

Lu J, Rincon N, Wood DE, *et al*. Metagenome analysis using the Kraken software suite. *Nature protocols* 2022; **17**:2815-2839. https://www.nature.com/articles/s41596-022-00738-y
Macleod KL, Jenkins TL, Witt MJ *et al*. Rare, long-distance dispersal underpins genetic connectivity in the pink sea fan, Eunicella verrucosa. *Evolutionary Applications* 2024; **17**:e13649. https://onlinelibrary.wiley.com/doi/abs/10.1111/eva.13649

Matz MV. Not-so-mutually beneficial coral symbiosis. *Current Biology* 2024; **34**(17):R798-R801.

Mayr E. Wu's genic view of speciation. *Journal of Evolutionary Biology* 2001; **14**:866–867.
https://academic.oup.com/jeb/article-abstract/14/6/866/7322934
McFadden CS, Benayahu Y, Pante E, et al. Limitations of mitochondrial gene barcoding in Octocorallia. *Molecular Ecology Resources* 2011; **11**:19–31.

https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1755-0998.2010.02875.x

McFadden CS, Quattrini AM, Brugler MR, et al. Phylogenomics, origin, and diversification of Anthozoans (Phylum Cnidaria). *Systematic Biology* 2021; **70**:635-647.

https://academic.oup.com/sysbio/article-abstract/70/4/635/6122449

Muir PR, Obura DO, Hoeksema BW, et al. Conclusions of low extinction risk for most species of reef-building corals are premature. *Nature Ecology & Evolution* 2022; **6**:357–358.

https://www.nature.com/articles/s41559-022-01659-5

Munro L. Determining the reproductive cycle of Eunicella verrucosa. *Reef Research: ETR* 2004; **11**.

https://www.marine-bio-images.com/RR_Eunicella_PDFS/Report_RR12Jul2004reproductive%20cycle%20pdf.pdf

Muthye V, Mackereth CD, Stewart JB et al. Large dataset of octocoral mitochondrial genomes provides new insights into mt-mutS evolution and function. *DNA repair* 2022; **110**:103273. https://www.sciencedirect.com/science/article/pii/S1568786422000027

van Oppen MJH, Medina M. Coral evolutionary responses to microbial symbioses. *Philosophical Transactions of the Royal Society B: Biological Sciences* 2020; **375**:20190591.

https://royalsocietypublishing.org/doi/abs/10.1098/rstb.2019.0591

Palumbi SR. Marine speciation on a small planet. *Trends in Ecology & Evolution* 1992; **7**:114–118. https://www.sciencedirect.com/science/article/pii/016953479290144Z

Pante E, Puillandre N, Viricel A, et al. Species are hypotheses: avoid connectivity assessments based on pillars of sand. *Molecular Ecology* 2015a; **24**:525–544.

https://onlinelibrary.wiley.com/doi/abs/10.1111/mec.13048

Pante E, Abdelkrim J, Viricel A, et al. Use of RAD sequencing for delimiting species. *Heredity* 2015b; **114**:450–459. https://www.nature.com/articles/hdy2014105

Pante E, Simon-Bouhet B. marmap: a package for importing, plotting and analyzing bathymetric and topographic data in R. *PLoS one* 2013; **8**:e73051.

https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0073051

Patro R, Duggal G, Love MI, et al. Salmon provides fast and bias-aware quantification of transcript expression. *Nature Methods* 2017; **14**:417–419.

https://www.nature.com/articles/nmeth.4197

Pelosi JA, Bernal MA, Krabbenhoft TJ, et al. Fine-scale morphological, genomic, reproductive, and symbiont differences delimit the Caribbean octocorals *Plexaura homomalla* and *P. kükenthali*. *Coral Reefs*. 2022; **41**:635–653.

https://link.springer.com/article/10.1007/s00338-021-02175-x

Peñalba JV, Runemark A, Meier JI, Singh P, Wogan GO, Sánchez-Guillén R, Mallet J, Rometsch SJ, Menon M, Seehausen O. The Role of Hybridization in Species Formation and Persistence. *Cold Spring Harbor perspectives in biology* 2014; **a041445**.

https://cshperspectives.cshlp.org/content/early/2024/03/01/cshperspect.a041445.short
Picard Toolkit. Broad Institute. GitHub Repository. https://broadinstitute.github.io/picard/. 2019

Porro B. Diversités génétiques chez l'holobiote Anemonia viridis: des morphotypes de l'hôte à la différenciation symbiotique. (Doctoral dissertation, COMUE Université Côte d'Azur (2015-2019). 2019. https://theses.hal.science/tel-02736573

Pudlo P, Marin JM, Estoup A, et al. Reliable ABC model choice via random forests. *Bioinformatics* 2016; **32**:859–866.

https://academic.oup.com/bioinformatics/article-abstract/32/6/859/1744513
Readman J, Hiscock K. Eunicella verrucosa. Pink sea fan.

https://www.marlin.ac.uk/species/detail/1121. 2017

Reynes L, Aurelle D, Chevalier C, et al. Population genomics and Lagrangian modeling shed light on dispersal events in the Mediterranean endemic Ericaria zosteroides (= Cystoseira zosteroides) (Fucales). *Frontiers in Marine Science* 2021;8:683528.

https://www.frontiersin.org/articles/10.3389/fmars.2021.683528/full

Rivera-García L, Rivera-Vicéns RE, Veglia AJ et al. De novo transcriptome assembly of the digitate morphotype of Briareum asbestinum (Octocorallia: Alcyonacea) from the southwest shelf of Puerto Rico. *Marine Genomics* 2019;**47**:100676.

https://www.sciencedirect.com/science/article/pii/S1874778718302393

Rosenberg E, Zilber-Rosenberg I. The hologenome concept of evolution after 10 years.

*Microbiome* 2018;6:78. https://link.springer.com/article/10.1186/S40168-018-0457-9

Rousset F. *genepop'007: a complete re-implementation of the genepop software for Windows and Linux*. *Molecular Ecology Resources* 2008;**8**:103–106.

https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1471-8286.2007.01931.x

Rousset F, Lopez J, Belkhir K. Package ‘genepop’. *R package version* 1. 2020. https://cran.r-project.org/web/packages/genepop/index.html

Roux C, Fraïsse C, Romiguier J, et al. Shedding Light on the Grey Zone of Speciation along a Continuum of Genomic Divergence. *PLOS Biology* 2016;**14**:e2000234.

https://journals.plos.org/plosbiology/article?id=10.1371/journal.pbio.2000234

Sachs JL, Wilcox TP. A shift to parasitism in the jellyfish symbiont *Symbiodinium microadriaticum*. *Proceedings of the Royal Society B: Biological Sciences* 2006;**273**:425–429. https://royalsocietypublishing.org/doi/abs/10.1098/rspb.2005.3346

Sartoretto S, Francour P. Bathymetric distribution and growth rates of Eunicella verrucosa (Cnidaria: Gorgoniidae) populations along the Marseilles coast (France). *Scientia Marina* 2011;**76**:349–355. https://archimer.ifremer.fr/doc/00087/19859/

Schulz MH, Zerbino DR, Vingron M et al. Oases: robust de novo RNA-seq assembly across the dynamic range of expression levels. *Bioinformatics* 2012;**28**:1086–1092.

https://academic.oup.com/bioinformatics/article-abstract/28/8/1086/195757

Sini M, Kipson S, Linares C, et al. The Yellow Gorgonian Eunicella cavolini: Demography and Disturbance Levels across the Mediterranean Sea. *PLoS ONE* 2015;**10**:e0126253.

https://journals.plos.org/plosone/article?id=10.1371/journal.pone.0126253

Stanley RRE, Jeffery NW, Wringe BF, et al. GENEPOPEDIT: a simple and flexible tool for manipulating multilocus molecular data in R. *Molecular Ecology Resources* 2017;**17**:12–18.

https://onlinelibrary.wiley.com/doi/abs/10.1111/1755-0998.12569

Tricou T, Tannier E, de Vienne DM. Ghost lineages can invalidate or even reverse findings regarding gene flow. *PLoS Biology* 2022;**20**:e3001776.

https://journals.plos.org/plosbiology/article?id=10.1371/journal.pbio.3001776

Tsagkogeorga G, Cahais V, Galtier N. The population genomics of a fast evolver: high levels of diversity, functional constraint, and molecular adaptation in the tunicate Ciona intestinalis.

*Genome biology and evolution* 2012;**4**:852–861. https://academic.oup.com/gbe/article-abstract/4/8/852/580636

Van Oppen M, Mieog JC, Sanchez C, et al. Diversity of algal endosymbionts (zooxanthellae) in octocorals: the roles of geography and host relationships. *Molecular Ecology* 2005;**14**:2403–2417. https://onlinelibrary.wiley.com/doi/abs/10.1111/j.1365-294X.2005.02545.x
Wagner D, Pochon X, Irwin L, et al. Azooxanthellate? Most Hawaiian black corals contain Symbiodinium. *Proceedings of the Royal Society B: Biological Sciences* 2011;**278**:1323–1328. https://royalsocietypublishing.org/doi/abs/10.1098/rspb.2010.1681

Weir BS, Cockerham CC. Estimating F-statistics for the analysis of population structure. *Evolution* 1984;**38**:1358–1370. https://www.jstor.org/stable/2408641

Wood DE, Lu J, Langmead B. Improved metagenomic analysis with Kraken 2. *Genome Biology* 2019;**20**:257. https://link.springer.com/article/10.1186/s13059-019-1891-0
Figure 1: map of sampling sites for transcriptomes: A) general view, B) zoom on the area of Marseille. The symbols correspond to different samples: EC *E. cavolini*, ES *E. singularis*, EV *E. verrucosa*, HY potential hybrids. The three letters correspond to the codes of the sampling. The maps have been produced with the marmap R package (Pante & Simon-Bouhet, 2013) and following the tutorial of Krueger-Hadfield (2015).

629

A)

<figure><img src="image_3.png" /></figure>

630

631

B)

632

<figure><img src="image_4.png" /></figure>

633

634
Figure 2: distribution of the frequency of Symbiodiniaceae sequences in the individual transcriptomes according to the species based A) on the number of reads estimated with Salmon, and B) on the proportion of assembled sequences (contigs) with the BLAT analyses.

639

A) Read counts with Salmon; mean values per group: *E. cavolini*: 16508; hybrids: 10238; *E. singularis*: 26023; *E. verrucosa*: 4285. Kruskal-Wallis test of the differences among groups: chi-squared = 7.9467, df = 3, p-value = 0.047.

<figure><img src="image_5.png" /></figure>

643

644

B) Assembled sequences with BLAT; mean values per group: *E. cavolini*: 0.0034; hybrids: 0.0029; *E. singularis*: 0.0219; *E. verrucosa*: 0.0028. Kruskal-Wallis test of the differences among groups: chi-squared = 14.352, df = 3, p-value = 0.002.

<figure><img src="image_6.png" /></figure>

648
649 **Figure 3:** barplots of coancestry coefficients inferred with the LEA R package. The analysis 650 is based on the "polymorphic sites" transcriptome dataset. The red asterisks indicate the 651 individuals used as prior for parental status in the NewHybrids analysis. The results of the 652 NewHybrids analysis are indicated below the hybrid individuals: F1, 1st generation; F2, 2nd 653 generation; Sbx, backcross with *E. singularis*. The coancestry analysis is based on 31 369 654 SNPs, whereas the NewHybrids analysis is based on 326 SNPs showing high differentiation 655 between *E. cavolini* and *E. singularis*. 656 657
<figure><img src="image_7.png" /></figure>

E. cavolini

hybrids

E. singularis

E. verrucosa
<figure><img src="image_8.png" /><figcaption>**Figure 4:** principal Component Analysis based on the "polymorphic SNPs" transcriptome dataset; axis 1 represents 33.2% of the variance, axis 2 represents 13% of the variance</figcaption></figure>

659
660
661
662

663
664
665
**Table 1** : inference of hybrid status with NewHybrids for transcriptome and RAD sequencing.
For transcriptomes, all probabilities were at 1 for the inferred status and for the ten replicates. For RAD sequencing, the results are given for the four datasets (different assembly strategies). If no probability is mentioned for RAD sequencing, the hybrid status was supported by a probability higher than 0.999 over the ten replicates. In the other cases, the numbers indicate the minimal probability threshold over the ten replicates for this status (and the status was coherent over the ten replicates as well, with slight variations in probability). NA indicates an individual which was removed during the filtering of SNPs because of too many missing data. The lines highlighted in grey indicate the cases where different status were inferred depending on the dataset. Bx-ES and Bx-EC indicate backcrosses with *E. singularis* and *E. cavolini* respectively; ES indicates parental *E. singularis*.

<table><thead><tr><th>Individual - RAD sequencing</th><th>de novo</th><th>ref. E. cavolini</th><th>ref. E. singularis</th><th>ref. E. verrucosa</th></tr></thead><tbody><tr><td>EC-X-MFNB</td><td>F2</td><td>F2</td><td>F2</td><td>F2</td></tr><tr><td>EC-X-MFNC</td><td>F2</td><td>NA</td><td>Bx-EC</td><td>NA</td></tr><tr><td>EC-X-MFND</td><td>Bx-ES</td><td>Bx-ES</td><td>Bx-ES</td><td>Bx-ES</td></tr><tr><td>EC-X-MFNE</td><td>Bx-EC</td><td>Bx-EC</td><td>Bx-EC</td><td>Bx-EC</td></tr><tr><td>EC-X-MFNF</td><td>Bx-EC</td><td>F2 > 0.95</td><td>Bx-EC > 0.92</td><td>F2</td></tr><tr><td>EC-X-MFNG</td><td>F2</td><td>F2</td><td>F2</td><td>F2</td></tr><tr><td>EC-X-MFNH</td><td>Bx-ES</td><td>Bx-ES</td><td>Bx-ES > 0.67</td><td>Bx-ES > 0.98</td></tr><tr><td>EC-X-MFNI</td><td>F1</td><td>F1</td><td>F1</td><td>F2</td></tr><tr><td>EC-X-MFNL</td><td>F1</td><td>F1 > 0.99</td><td>F1</td><td>F2</td></tr><tr><td>ES-X-MFNA</td><td>Bx-ES</td><td>Bx-ES</td><td>Bx-ES</td><td>Bx-ES</td></tr><tr><td>ES-X-MFNJ</td><td>F2</td><td>F2 > 0.96</td><td>F2</td><td>F2</td></tr><tr><td>ES-X-MFNK</td><td>ES</td><td>ES</td><td>ES</td><td>ES</td></tr></tbody></table>

Individual - transcriptome

<table><tr><td>EH-JPB-a</td><td>F1</td></tr><tr><td>EH-MFN-a</td><td>Bx-ES</td></tr><tr><td>EH-MFN-b</td><td>F2</td></tr><tr><td>EH-MFN-e</td><td>Bx-ES</td></tr></table>
**Table 2:** results of demographic inferences with DILS with transcriptome data. The columns indicate the species comparison, the model choice for population size (constant vs. variable), and the results of inferences: current (on-going) gene flow (migration vs isolation); if current migration is inferred DILS compares isolation / migration (IM) with secondary contact (SC); if no current migration is inferred, the comparison is between strict isolation (SI) and ancestral migration (AM); the last columns give the results of the tests of homogeneity or heterogeneity among loci for inferences in effective size (N-homo vs N-hetero), and gene flow (M-homo vs M-hetero). The probability of each scenario is given in the same case. Homogeneity and heterogeneity indicate no variation or variation among loci respectively.

<table><thead><tr><th>Comparison</th><th>Population size</th><th>Current gene flow</th><th>IM / SC</th><th>SI / AM</th><th>N-hetero / N-homo</th><th>M-hetero / M-homo</th></tr></thead><tbody><tr><td>cavolini / singularis</td><td>constant</td><td>Migration; 0.87</td><td>SC; 0.79</td><td>-</td><td>N-hetero; 0.99</td><td>M-homo; 0.82</td></tr><tr><td>cavolini / singularis</td><td>variable</td><td>Migration; 0.88</td><td>SC; 0.77</td><td>-</td><td>N-hetero; 1</td><td>M-homo; 0.87</td></tr><tr><td>cavolini / verrucosa</td><td>constant</td><td>Isolation; 0.90</td><td>-</td><td>AM; 0.65</td><td>N-hetero; 1</td><td>-</td></tr><tr><td>cavolini / verrucosa</td><td>variable</td><td>Isolation; 0.89</td><td>-</td><td>AM; 0.69</td><td>N-hetero; 1</td><td>-</td></tr><tr><td>singularis / verrucosa</td><td>constant</td><td>Isolation; 0.87</td><td>-</td><td>AM; 0.61</td><td>N-hetero; 1</td><td>-</td></tr><tr><td>singularis / verrucosa</td><td>variable</td><td>Isolation; 0.87</td><td>-</td><td>AM; 0.61</td><td>N-hetero; 1</td><td>-</td></tr></tbody></table>