# Symbiotic status does not preclude hybridisation in Mediterranean octocorals

Didier Aurelle, Anne Haguenauer, Marc Bally, Frédéric Zuberer, Jean-Baptiste Ledoux, Stéphane Sartoretto, Cédric Cabau, Lamya Chaoui, Hichem Kara, Sarah Samadi, et al.

► To cite this version:

Didier Aurelle, Anne Haguenauer, Marc Bally, Frédéric Zuberer, Jean-Baptiste Ledoux, et al. Symbiotic status does not preclude hybridisation in Mediterranean octocorals. 2024. hal-04568927

HAL Id: hal-04568927

https://hal.science/hal-04568927

Preprint submitted on 6 May 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
Symbiotic status does not
preclude hybridisation in
Mediterranean octocorals

Didier Aurelle<fnref n="1" />,<fnref n="2" />*, Anne Haguenauer<fnref n="3" />, Marc Bally<fnref n="1" />,
Frédéric Zuberer<fnref n="4" />, Jean-Baptiste Ledoux<fnref n="5" />, Stéphane
Sartoretto<fnref n="6" />, Cédric Cabau<fnref n="7" />, Lamya Chaoui<fnref n="8" />, Hichem
Kara<fnref n="8" />, Sarah Samadi<fnref n="2" />, Pierre Pontarotti<fnref n="9" />,<fnref n="10" />,<fnref n="11" />

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

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

<footnote n="1">¹ Aix Marseille Univ, Université de Toulon, CNRS, IRD, MIO, Marseille, France</footnote>

<footnote n="2">² Institut Systématique Evolution Biodiversité (ISYEB), Muséum national d'Histoire naturelle, CNRS, Sorbonne Université, EPHE, Université des Antilles, CP 26, 75005 Paris, France.</footnote>

<footnote n="3">³ CNRS - Délégation Provence et corse, Marseille, France</footnote>

<footnote n="4">⁴ Aix Marseille Univ, CNRS, IRD, INRAE, OSU Inst. PYTHEAS, Marseille, France</footnote>

<footnote n="5">⁵ CIIMAR/CIMAR, Centro Interdisciplinar de Investigação Marinha e Ambiental, Universidade do Porto, Porto, Portugal.</footnote>

<footnote n="6">⁶ Ifremer, LITTORAL, 83500 La Seyne-sur-Mer, France</footnote>

<footnote n="7">⁷ Sigenae, GenPhySE, Université de Toulouse, INRAE, ENVT, 31326, Castanet Tolosan, France.</footnote>

<footnote n="8">⁸ Laboratoire Bioressources marines. Université d'Annaba Badji Mokhtar, Annaba - Algérie.</footnote>

<footnote n="9">⁹ Aix Marseille Univ, MEPHI, Marseille, France.</footnote>

<footnote n="10">¹⁰ IHU Méditerranée Infection, Marseille, France.</footnote>

<footnote n="11">¹¹ CNRS SNC5039</footnote>

<footnote>*Corresponding author</footnote>

<footnote>Correspondence: [didier.aurelle@univ-amu.fr](mailto:didier.aurelle@univ-amu.fr)</footnote>
Symbiotic status does not
preclude hybridisation in
Mediterranean octocorals

Didier Aurelle<fnref n="1" />,<fnref n="2" />*, Anne Haguenauer<fnref n="3" />, Marc Bally<fnref n="1" />,
Frédéric Zuberer<fnref n="4" />, Jean-Baptiste Ledoux<fnref n="5" />, Stéphane
Sartoretto<fnref n="6" />, Cédric Cabau<fnref n="7" />, Lamya Chaoui<fnref n="8" />, Hichem
Kara<fnref n="8" />, Sarah Samadi<fnref n="2" />, Pierre Pontarotti<fnref n="9" />,<fnref n="10" />,<fnref n="11" />

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

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

<footnote n="1">¹ Aix Marseille Univ, Université de Toulon, CNRS, IRD, MIO, Marseille, France</footnote>

<footnote n="2">² Institut Systématique Evolution Biodiversité (ISYEB), Muséum national d'Histoire naturelle, CNRS, Sorbonne Université, EPHE, Université des Antilles, CP 26, 75005 Paris, France.</footnote>

<footnote n="3">³ CNRS - Délégation Provence et corse, Marseille, France</footnote>

<footnote n="4">⁴ Aix Marseille Univ, CNRS, IRD, INRAE, OSU Inst. PYTHEAS, Marseille, France</footnote>

<footnote n="5">⁵ CIIMAR/CIMAR, Centro Interdisciplinar de Investigação Marinha e Ambiental, Universidade do Porto, Porto, Portugal.</footnote>

<footnote n="6">⁶ Ifremer, LITTORAL, 83500 La Seyne-sur-Mer, France</footnote>

<footnote n="7">⁷ Sigenae, GenPhySE, Université de Toulouse, INRAE, ENVT, 31326, Castanet Tolosan, France.</footnote>

<footnote n="8">⁸ Laboratoire Bioressources marines. Université d'Annaba Badji Mokhtar, Annaba - Algérie.</footnote>

<footnote n="9">⁹ Aix Marseille Univ, MEPHI, Marseille, France.</footnote>

<footnote n="10">¹⁰ IHU Méditerranée Infection, Marseille, France.</footnote>

<footnote n="11">¹¹ CNRS SNC5039</footnote>

<footnote>*Corresponding author</footnote>

<footnote>Correspondence: didier.aurelle@univ-amu.fr</footnote>
**ABSTRACT**

Understanding how species can form and remain isolated in the marine
environment is still an active research area. Correctly delimiting species if
also of interest for biodiversity conservation and for a wide range of
biological studies. Here we study the differentiation and the possibility of
hybridization among three temperate octocorals : *Eunicella cavolini*,
*E. singularis* and *E. verrucosa*. These species can occur in sympatry and
morphologically intermediate individuals have been observed. Among
these three species, *E. singularis* is the only one known to show a
mutualistic symbiosis with photosynthetic Symbiodiniaceae. As symbiotic
relationships can be involved in species barriers, 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 current 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. Our data suggest an intermediate abundance of Symbiodiniaceae
in the hybrids of the two species. We discuss the evolution of the
Symbiodiniaceae / cnidarian symbiosis in the light of our results.

**Keywords:** speciation, hybridization, symbiosis, transcriptome, octocoral
**ABSTRACT**

Understanding how species can form and remain isolated in the marine
environment is still an active research area. Correctly delimiting species if
also of interest for biodiversity conservation and for a wide range of
biological studies. Here we study the differentiation and the possibility of
hybridization among three temperate octocorals : *Eunicella cavolini*,
*E. singularis* and *E. verrucosa*. These species can occur in sympatry and
morphologically intermediate individuals have been observed. Among
these three species, *E. singularis* is the only one known to show a
mutualistic symbiosis with photosynthetic Symbiodiniaceae. As symbiotic
relationships can be involved in species barriers, 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 current 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. Our data suggest an intermediate abundance of Symbiodiniaceae
in the hybrids of the two species. We discuss the evolution of the
Symbiodiniaceae / cnidarian symbiosis in the light of our results.

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

2

3 As corner stones of evolutionary biology, species and speciation still raise a
4 wealth of questions fueled by the technological and conceptual
5 advancements in genomics. Species can be defined as a part of a
6 genealogical network, and such definition should be clearly distinguished
7 from recognition criteria (Samadi & Barberousse, 2006). Genomic data, either
8 from complete or partial representation of genomes, allow testing species
9 hypotheses, understanding speciation scenarios, and go deeper in the
10 analysis of interactions between biodiversity and environments. Species are
11 indeed hypotheses to be tested, if possible through the integration of
12 independent criteria, and independent molecular markers (Pante et al.,
13 2015b). Sound species delimitation and identification is useful, among others,
14 to better estimatespecies range and biodiversity patterns (Coelho et al.,
15 2023; Muir et al., 2022), to avoid biases in connectivity studies (Pante et al.,
16 2015b), and adaptive abilities (Brener-Raffalli et al., 2022). Species
17 delimitation can be problematic in the context of the grey zone of speciation,
18 where different delimitation criteria may bring contradictory conclusions (De
19 Queiroz, 2007). This is the case for cryptic species, i.e. entities assigned to
20 the same nominal species though reproductively isolated (Cahill et al., 2023).
21 Conversely, gene flow can be observed among morphologically differentiated
22 entities, and this can be an important factor in adaptation and evolution
23 (Leroy et al., 2020). The analysis of a high number of independent markers
24 now provides the required analytical power for species delimitation in such
25 difficult cases, for example when mitochondrial barcoding is not informative
26 (McFadden et al., 2011). This wealth of information allows demographic
27 inferences and testing speciation scenarios (De Jode et al., 2023; Roux et al.,
28 2016). Incomplete speciation might still be inferred, and it's an interesting
29 situation to study the speciation process itself (Roux et al., 2016).

30 In the marine realm, the question of speciation was considered as particularly
31 puzzling. From a fundamental point of view, understanding the drivers of
32 speciation is not easy for species with high effective size, and important gene
33 flow in a supposedly open environment (Faria et al., 2021; Mayr, 2001;
34 Palumbi, 1992). Difficulties in sampling and rearing organisms also limits
35 experiments to test reproductive isolation (Faria et al., 2021). Important
36 progress has nevertheless been made in various marine organisms to better
37 understand spatial patterns of genetic structure. One can use models of
38 oceanographic connectivity to understand the observed genetic
39 differentiation, and highlight potential barriers to gene flow (Reynes et al.,
40 2021). Barriers to gene flow can also be inferred from clines in allele
41 frequencies (Gagnaire et al., 2015). The observation of hybrid zones in
42 marine species unveils heterogeneity in this environment, and shows the
43 importance of the interactions between exogenous (i.e. dependent on
44 different environments) and intrinsic (genetic incompatibilities) barriers to
45 gene flow (Bierne et al., 2011). Some marine species also provided important
46 models to study the role of gametes recognition in reproductive isolation
47 (Palumbi, 1999).

48 In this context, a poorly investigated topics remains the role of symbiotic
49 interaction in speciation. There are various examples of the involvement of
50 microbial species in reproductive isolation, especially in insects (Brucker &
51 Bordenstein, 2012). For marine species, the role of microbial communities has
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

38

39

40

41

42

43

44

45

46

47

48

49

50

51

As corner stones of evolutionary biology, species and speciation still raise a wealth of questions fueled by the technological and conceptual advancements in genomics. Species can be defined as a part of a genealogical network, and such definition should be clearly distinguished from recognition criteria (Samadi & Barberousse, 2006). Genomic data, either from complete or partial representation of genomes, allow testing species hypotheses, understanding speciation scenarios, and go deeper in the analysis of interactions between biodiversity and environments. Species are indeed hypotheses to be tested, if possible through the integration of independent criteria, and independent molecular markers (Pante et al., 2015b). Sound species delimitation and identification is useful, among others, to better estimates species range and biodiversity patterns (Coelho et al., 2023; Muir et al., 2022), to avoid biases in connectivity studies (Pante et al., 2015b), and adaptive abilities (Brener-Raffalli et al., 2022). Species delimitation can be problematic in the context of the grey zone of speciation, where different delimitation criteria may bring contradictory conclusions (De Queiroz, 2007). This is the case for cryptic species, i.e. entities assigned to the same nominal species though reproductively isolated (Cahill et al., 2023).

Conversely, gene flow can be observed among morphologically differentiated entities, and this can be an important factor in adaptation and evolution (Leroy et al., 2020). The analysis of a high number of independent markers now provides the required analytical power for species delimitation in such difficult cases, for example when mitochondrial barcoding is not informative (McFadden et al., 2011). This wealth of information allows demographic inferences and testing speciation scenarios (De Jode et al., 2023; Roux et al., 2016). Incomplete speciation might still be inferred, and it's an interesting situation to study the speciation process itself (Roux et al., 2016).

In the marine realm, the question of speciation was considered as particularly puzzling. From a fundamental point of view, understanding the drivers of speciation is not easy for species with high effective size, and important gene flow in a supposedly open environment (Faria et al., 2021; Mayr, 2001; Palumbi, 1992). Difficulties in sampling and rearing organisms also limits experiments to test reproductive isolation (Faria et al., 2021). Important progress has nevertheless been made in various marine organisms to better understand spatial patterns of genetic structure. One can use models of oceanographic connectivity to understand the observed genetic differentiation, and highlight potential barriers to gene flow (Reynes et al., 2021). Barriers to gene flow can also be inferred from clines in allele frequencies (Gagnaire et al., 2015). The observation of hybrid zones in marine species unveils heterogeneity in this environment, and shows the importance of the interactions between exogenous (i.e. dependent on different environments) and intrinsic (genetic incompatibilities) barriers to gene flow (Bierne et al., 2011). Some marine species also provided important models to study the role of gametes recognition in reproductive isolation (Palumbi, 1999).

In this context, a poorly investigated topics remains the role of symbiotic interaction in speciation. There are various examples of the involvement of microbial species in reproductive isolation, especially in insects (Brucker & Bordenstein, 2012). For marine species, the role of microbial communities has
been mainly explored at the light of adaptative evolution, with the concept of
hologenome which comprises the genomes of the host and associated
microbiome (Rosenberg & Zilber-Rosenberg, 2018). Shallow water
scleractinian corals (hexacorals) are usually associated with photosynthetic
zooxanthellae (Cairns, 2007). These zooxanthellae, of the Symbiodiniaceae
family, correspond to different genera and species (LaJeunesse et al., 2018).
Changes in associated Symbiodiniaceae can impact the thermotolerance of
the coral holobiont, and the possibility of adaptation facing climate change
(Berkelmans & van Oppen, 2006; van Oppen & Medina, 2020). Inferences
from the phylogeny of scleractinian corals have shown multiple acquisitions
of the symbiotic state, but there is still an important diversity of non
symbiotic corals (Cairns, 2007; Campoy et al., 2020). Despite the diversity of
Symbiodiniaceae, the diversity in symbiotic state (presence / absence,
species diversity), and the possibility of shifting at different evolutionary
scales, there is no information on the potential involvement of these
symbionts in reproductive isolation, or even on the possibility of gene flow
among species with different symbiotic state. The symbiotic interactions
between Anthozoans (hexacorals and octocorals) and Symbiodiniaceae
presents important mutualistic benefits especially from a nutritional point of
view (Furla et al., 2005). This interaction requires specific adaptations for the
animal host, as for example protection against oxygen produced by
photosynthesis (Furla et al., 2005). The association with Symbiodiniaceae can
also range from mutualism to parasitism (Lesser et al., 2013; Sachs & Wilcox,
2006). Therefore, one can envision that hybrids could be impaired by a
modification in host - symbiont interactions. The presence or type of
Symbiodiniaceae could also be involved in genetic incompatibilities with the
host genome, as previously observed with bacterial species (Bordenstein,
2003; Brucker & Bordenstein, 2012).
Here we will study species limits between species with different symbiotic
states in *Eunicella* gorgonians (octocorals). In shallow conditions (above 50 m
depth), three *Eunicella* species are mainly present in the Mediterranean sea:
*Eunicella cavolini*, *E. singularis*, and *E. verrucosa*. These three species have
partially overlapping ranges, and they can be observed in sympatry, as is the
case in the area of Marseille (France). *Eunicella singularis* hosts
Symbiodiniaceae corresponding to temperate clade A (Forcioli et al., 2011;
LaJeunesse et al., 2018; Porro, 2019), whereas the two other species are
devoided of these symbionts (Carpine & Grasshoff, 1975). The
Symbiodiniaceae contribute to the carbon metabolism of *E. singularis*, but a
non-symbiotic *aphyta* morph of this species has already been observed (Gori
et al., 2012). While the lack of variability in mitochondrial DNA does not allow
to distinguish these three species (Calderon et al., 2006), an initial study
using two nuclear introns, and considering morphologically intermediate
individuals, suggested the possibility of hybridisation between *E. singularis*
and *E. verrucosa* (Aurelle et al., 2017). Demographic inferences based on the
transcriptome sequences of *E. cavolini* and *E. verrucosa* indicated the
possibility of current gene flow between these two species, but *E. singularis*
was not analysed there (Roux et al., 2016). Here, we will go further on these
questions with the following objectives: i) estimate the genomic
differentiation among these three species, ii) test the possibility of
hybridisation according to symbiotic state and genetic similarities, and iii)
infer scenarios of speciation.
been mainly explored at the light of adaptative evolution, with the concept of
hologenome which comprises the genomes of the host and associated
microbiome (Rosenberg & Zilber-Rosenberg, 2018). Shallow water
scleractinian corals (hexacorals) are usually associated with photosynthetic
zooxanthellae (Cairns, 2007). These zooxanthellae, of the Symbiodiniaceae
family, correspond to different genera and species (LaJeunesse et al., 2018).
Changes in associated Symbiodiniaceae can impact the thermotolerance of
the coral holobiont, and the possibility of adaptation facing climate change
(Berkelmans & van Oppen, 2006; van Oppen & Medina, 2020). Inferences
from the phylogeny of scleractinian corals have shown multiple acquisitions
of the symbiotic state, but there is still an important diversity of non
symbiotic corals (Cairns, 2007; Campoy et al., 2020). Despite the diversity of
Symbiodiniaceae, the diversity in symbiotic state (presence / absence,
species diversity), and the possibility of shifting at different evolutionary
scales, there is no information on the potential involvement of these
symbionts in reproductive isolation, or even on the possibility of gene flow
among species with different symbiotic state. The symbiotic interactions
between Anthozoans (hexacorals and octocorals) and Symbiodiniaceae
presents important mutualistic benefits especially from a nutritional point of
view (Furla et al., 2005). This interaction requires specific adaptations for the
animal host, as for example protection against oxygen produced by
photosynthesis (Furla et al., 2005). The association with Symbiodiniaceae can
also range from mutualism to parasitism (Lesser et al., 2013; Sachs & Wilcox,
2006). Therefore, one can envision that hybrids could be impaired by a
modification in host - symbiont interactions. The presence or type of
Symbiodiniaceae could also be involved in genetic incompatibilities with the
host genome, as previously observed with bacterial species (Bordenstein,
2003; Brucker & Bordenstein, 2012).
Here we will study species limits between species with different symbiotic
states in *Eunicella* gorgonians (octocorals). In shallow conditions (above 50 m
depth), three *Eunicella* species are mainly present in the Mediterranean sea:
*Eunicella cavolini*, *E. singularis*, and *E. verrucosa*. These three species have
partially overlapping ranges, and they can be observed in sympatry, as is the
case in the area of Marseille (France). *Eunicella singularis* hosts
Symbiodiniaceae corresponding to temperate clade A (Forcioli et al., 2011;
LaJeunesse et al., 2018; Porro, 2019), whereas the two other species are
devoided of these symbionts (Carpine & Grasshoff, 1975). The
Symbiodiniaceae contribute to the carbon metabolism of *E. singularis*, but a
non-symbiotic *aphyta* morph of this species has already been observed (Gori
et al., 2012). While the lack of variability in mitochondrial DNA does not allow
to distinguish these three species (Calderon et al., 2006), an initial study
using two nuclear introns, and considering morphologically intermediate
individuals, suggested the possibility of hybridisation between *E. singularis*
and *E. verrucosa* (Aurelle et al., 2017). Demographic inferences based on the
transcriptome sequences of *E. cavolini* and *E. verrucosa* indicated the
possibility of current gene flow between these two species, but *E. singularis*
was not analysed there (Roux et al., 2016). Here, we will go further on these
questions with the following objectives: i) estimate the genomic
differentiation among these three species, ii) test the possibility of
hybridisation according to symbiotic state and genetic similarities, and iii)
infer scenarios of speciation.
104 **Material and methods**

105 **Species distribution**

106 *E. singularis* and *E. cavolini* are only present in the Mediterranean Sea,
107 whereas *E. verrucosa* is present both in the Eastern Atlantic Ocean and the
108 Mediterranean Sea (Carpine & Grasshoff, 1975). In the Atlantic, *E. verrucosa*
109 can be found from Ireland, West coasts of Britain to the north-west Africa
110 (Readman & Hiscock, 2017). *Eunicella verrucosa* is also present in the North
111 Western Mediterranean Sea with a patchy distribution, in Sardinia (Canessa et
112 al., 2022), and possibly in the Adriatic and Aegean Seas (Chimienti, 2020). In
113 the Mediterranean Sea, it can be observed from shallow conditions (20-40 m)
114 up to 200 m depth: this presence in relatively deep conditions can lead to
115 possible underestimation of this species range (Chimienti, 2020; Fourt &
116 Goujard, 2012; Sartoretto & Francour, 2011). *Eunicella cavolini* is present in
117 the Western Mediterranean, Adriatic and Aegean Seas, from 5 to 200 m depth
118 (Carugati et al., 2022; Sini et al., 2015). As previously mentioned,
119 *E. singularis* is the only Mediterranean octocoral known to habour
120 Symbiodiniaceae (but see Bonacolta et al., 2024): these Symbiodiniaceae
121 correspond to the temperate clade A (Casado-Amezúa et al., 2016; Forcioli
122 et al., 2011), now corresponding to the *Philozoon* genus, which is sister to
123 the *Symbiodinium* genus (LaJeunesse et al., 2018, 2022). *Eunicella singularis*
124 can be found in the Western Mediterranean and Adriatic Seas, and less
125 frequently in the Eastern Mediterranean (Gori et al., 2012). It is usually
126 observed up to top 40 m depth, but deeper occurrences (up to 70 m) have been
127 mentioned, which correspond to the azooxanthellate aphyta morph (without
128 Symbiodiniaceae; Gori et al., 2012). In the area of Marseille, these three
129 species can be observed in sympatry, sometimes at the same depth range,
130 up to 20 m depth (Sartoretto & Francour, 2011).

131

132 **Sampling**

133 Samples of *E. cavolini*, *E. singularis*, and *E. verrucosa* have been collected by
134 scuba diving in the Mediterranean (for the three species), and in the Atlantic
135 (*E. verrucosa* only; Figure 1) in 2016. In the area of Marseille, the three
136 species have been sampled in sympatry. Four morphologically intermediate
137 individuals (i.e. intermediate colors and branching patterns between
138 *E. cavolini* and *E. singularis*) were collected in the area of Marseille as well:
139 the objective was to test their hybrid status (Aurelle et al., 2017). Figure S1
140 presents some examples of morphology for the different species and
141 morphologically intermediate individuals in the area of Marseille. Sampling
142 was non destructive, with only a few centimeters of a branch collected at
143 each time. Sampling was performed with authorizations from the local
144 authorities.

145

146 **Phylogenetic relationships**

147 To put the three *Eunicella* species studied here in a phylogenetic context,
148 and to test their proximity, we performed a phylogenetic analysis by using
149 the mitochondrial MutS locus, which is widely used in octocorals
150 (McFadden et al., 2011), on the basis of sequences available in GenBank.
151 The methods of the phylogenetic reconstruction are detailed in the legend
152 of supplementary Figure S2, and the list of sequences in Table S1.

153

154 **Transcriptome sequencing and assembly**
104 **Material and methods**

105 **Species distribution**

106 *E. singularis* and *E. cavolini* are only present in the Mediterranean Sea,
107 whereas *E. verrucosa* is present both in the Eastern Atlantic Ocean and the
108 Mediterranean Sea (Carpine & Grasshoff, 1975). In the Atlantic, *E. verrucosa*
109 can be found from Ireland, West coasts of Britain to the north-west Africa
110 (Readman & Hiscock, 2017). *Eunicella verrucosa* is also present in the North
111 Western Mediterranean Sea with a patchy distribution, in Sardinia (Canessa et
112 al., 2022), and possibly in the Adriatic and Aegean Seas (Chimienti, 2020). In
113 the Mediterranean Sea, it can be observed from shallow conditions (20-40 m)
114 up to 200 m depth: this presence in relatively deep conditions can lead to
115 possible underestimation of this species range (Chimienti, 2020; Fourt &
116 Goujard, 2012; Sartoretto & Francour, 2011). *Eunicella cavolini* is present in
117 the Western Mediterranean, Adriatic and Aegean Seas, from 5 to 200 m depth
118 (Carugati et al., 2022; Sini et al., 2015). As previously mentioned,
119 *E. singularis* is the only Mediterranean octocoral known to habour
120 Symbiodiniaceae (but see Bonacolta et al., 2024): these Symbiodiniaceae
121 correspond to the temperate clade A (Casado-Amezúa et al., 2016; Forcioli
122 et al., 2011), now corresponding to the *Philozoon* genus, which is sister to
123 the *Symbiodinium* genus (LaJeunesse et al., 2018, 2022). *Eunicella singularis*
124 can be found in the Western Mediterranean and Adriatic Seas, and less
125 frequently in the Eastern Mediterranean (Gori et al., 2012). It is usually
126 observed up to top 40 m depth, but deeper occurrences (up to 70 m) have been
127 mentioned, which correspond to the azooxanthellate aphyta morph (without
128 Symbiodiniaceae; Gori et al., 2012). In the area of Marseille, these three
129 species can be observed in sympatry, sometimes at the same depth range,
130 up to 20 m depth (Sartoretto & Francour, 2011).

131

132 **Sampling**

133 Samples of *E. cavolini*, *E. singularis*, and *E. verrucosa* have been collected by
134 scuba diving in the Mediterranean (for the three species), and in the Atlantic
135 (*E. verrucosa* only; Figure 1) in 2016. In the area of Marseille, the three
136 species have been sampled in sympatry. Four morphologically intermediate
137 individuals (i.e. intermediate colors and branching patterns between
138 *E. cavolini* and *E. singularis*) were collected in the area of Marseille as well:
139 the objective was to test their hybrid status (Aurelle et al., 2017). Figure S1
140 presents some examples of morphology for the different species and
141 morphologically intermediate individuals in the area of Marseille. Sampling
142 was non destructive, with only a few centimeters of a branch collected at
143 each time. Sampling was performed with authorizations from the local
144 authorities.

145

146 **Phylogenetic relationships**

147 To put the three *Eunicella* species studied here in a phylogenetic context,
148 and to test their proximity, we performed a phylogenetic analysis by using
149 the mitochondrial MutS locus, which is widely used in octocorals
150 (McFadden et al., 2011), on the basis of sequences available in GenBank.
151 The methods of the phylogenetic reconstruction are detailed in the legend
152 of supplementary Figure S2, and the list of sequences in Table S1.

153

154 **Transcriptome sequencing and assembly**
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 supplementary material, Table S2.

We used the BLAT software (Kent, 2002) to remove potential Symbiodiniaceae sequences in the obtained transcriptomes. We used the transcriptome of the type A1 *Symbiodinium* (Genbank accession number GAKY01000000) (Baumgarten et al., 2013) to search for Symbiodiniaceae sequences separately in individual *Eunicella* transcriptomes. The output of BLAT was analysed with the *blat_parser.pl* script.

## Analysis of the presence of Symbiodiniaceae

We analysed the presence of Symbiodiniaceae in *Eunicella* gorgonians in two ways. First, we counted the number of reads corresponding to the aforementioned Symbiodiniaceae transcriptome with Salmon (Patro et al., 2017). Second, we used the percentage of assembled sequences (i.e. 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) on individual transcriptomes to try to identify the Symbiodiniaceae genera present in the different samples.

As our results pointed to the potential and unexpected presence of Symbiodiniaceae in *E. cavolini* (see Results), we wanted to explore this topic further by using the 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*. This experiment and its results are presented in Supplementary File S2.

## SNPs calling and filtering

We produced bam files by mapping the reads with bwa option mem (Li & Durbin, 2009). We used as reference the meta transcriptome filtered for Symbiodiniaceae sequences. 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 (Gayral et al., 2013; Tsagkogeorga et al., 2012). 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 (see below), and will be referred as “all-CS” for the *cavolini / singularis* comparison, “all-CV” for the *cavolini / verrucosa* comparison, and “all-SV” for the *singularis / verrucosa* comparison.

We filtered the “all sites” vcf file obtained with reads2snp with vcftools (Danecek et al., 2011) for an analysis of genetic diversity and differentiation.

We only retained biallelic sites which have been genotyped in all individuals
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 supplementary material, Table S2.

We used the BLAT software (Kent, 2002) to remove potential Symbiodiniaceae sequences in the obtained transcriptomes. We used the transcriptome of the type A1 *Symbiodinium* (Genbank accession number GAKY01000000) (Baumgarten et al., 2013) to search for Symbiodiniaceae sequences separately in individual *Eunicella* transcriptomes. The output of BLAT was analysed with the *blat_parser.pl* script.

## Analysis of the presence of Symbiodiniaceae

We analysed the presence of Symbiodiniaceae in *Eunicella* gorgonians in two ways. First, we counted the number of reads corresponding to the aforementioned Symbiodiniaceae transcriptome with Salmon (Patro et al., 2017). Second, we used the percentage of assembled sequences (i.e. 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) on individual transcriptomes to try to identify the Symbiodiniaceae genera present in the different samples.

As our results pointed to the potential and unexpected presence of Symbiodiniaceae in *E. cavolini* (see Results), we wanted to explore this topic further by using the 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*. This experiment and its results are presented in Supplementary File S2.

## SNPs calling and filtering

We produced bam files by mapping the reads with bwa option mem (Li & Durbin, 2009). We used as reference the meta transcriptome filtered for Symbiodiniaceae sequences. 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 (Gayral et al., 2013; Tsagkogeorga et al., 2012). 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 (see below), and will be referred as “all-CS” for the *cavolini / singularis* comparison, “all-CV” for the *cavolini / verrucosa* comparison, and “all-SV” for the *singularis / verrucosa* comparison.

We filtered the “all sites” vcf file obtained with reads2snp with vcftools (Danecek et al., 2011) for an analysis of genetic diversity and differentiation.

We only retained biallelic sites which have been genotyped in all individuals
207 (no missing data). To minimize correlation among loci, we retained SNPs
208 separated by at least 1 kb: we will thereafter refer to this dataset as the
209 "polymorphic sites" dataset. From this "polymorphic sites" dataset, we built a
210 dataset focused on the differentiation between *E. cavolini* and *E. singularis*:
211 we excluded *E. verrucosa* samples and we retained the first percent of the
212 loci with the highest F<sub>ST</sub> between *E. cavolini* and *E. singularis*. This last
213 dataset will be referred as "1% SNPs" dataset. The characteristics of the
214 different datasets are summarised in Table S3.
215

216

Genetic differentiation and analysis of hybrids

217 We analysed the genetic structure and differentiation among species with the
218 "polymorphic sites" dataset. We used the LEA R package to estimate
219 ancestry coefficients (Frichot et al., 2014; Frichot & François, 2015). We
220 tested K values from 1 to 10, with 10 replicates for each K. To analyse the
221 genetic differences among individuals, we performed a Principal Component
222 Analysis (PCA) with the R package adegenet. The pairwise F<sub>ST</sub> among species
223 were computed with the R package Genepop (Rousset, 2008; Rousset et al.,
224 2020), after conversion of the vcf file with PGDSpider (Lischer & Excoffier,
225 2012). We present here the F<sub>ST</sub> estimates of Weir & Cockerham (1984). The
226 distribution of F<sub>ST</sub> among loci was obtained with vcftools.

227 The hybrid status (e.g. first generation hybrids) of morphologically
228 intermediate individuals was analysed with the newhybrids software
229 (Anderson & Thompson, 2002). We used the genepopedit R package to
230 prepare the input file from genepop format (Stanley et al., 2017). Following
231 the results of the LEA and PCA analyses, we focused here on the comparison
232 between *E. cavolini*, *E. singularis* and potential hybrids. The newhybrids
233 analysis had difficulties to converge with such a high number of loci
234 compared to the number of markers (see discussion here:
235 https://github.com/erigande/newhybrids/issues/5). We therefore used the "1%
236 SNP" dataset (with the 1% highest F<sub>ST</sub> between *E. cavolini* and *E. singularis*
237 for the newhybrids analysis. As a prior, we used individuals with the lowest
238 levels of admixture in LEA as potential parental individuals: this corresponded
239 to three over five individuals for *E. cavolini*, and six over eight individuals for
240 *E. singularis* (see results). The newhybrids analysis was repeated five times
241 with different seeds to test the stability of the results.
242

243

Scenarios of speciation

244 We tested scenarios of speciation with the Demographic Inferences with Linked
245 Selection (DILS) pipeline (Csilléry et al., 2012; Fraïsse et al., 2021; Pudlo et al.,
246 2016). The DILS pipeline allows the analysis of two species scenarios only: we
247 therefore performed separate analyses for the three two-species
248 comparisons, with the "all-CS", "all-CV", and "all-SV" pairwise datasets. We
249 first tested the priors proposed by Monnet et al. (2023) according to observed
250 statistics of polymorphism and divergence: we then increased the range of
251 priors to get better results on the goodness-of-fit tests. We used the same
252 priors for all analyses, with different numbers of sequences per gene and per
253 sample according to the dataset (Table S4). For all pairwise comparisons, we
254 performed two DILS analyses: one with constant population sizes, and one
255 with variable population sizes.
256

256

Results

257
207 (no missing data). To minimize correlation among loci, we retained SNPs
208 separated by at least 1 kb: we will thereafter refer to this dataset as the
209 "polymorphic sites" dataset. From this "polymorphic sites" dataset, we built a
210 dataset focused on the differentiation between *E. cavolini* and *E. singularis*:
211 we excluded *E. verrucosa* samples and we retained the first percent of the
212 loci with the highest F<sub>ST</sub> between *E. cavolini* and *E. singularis*. This last
213 dataset will be referred as "1% SNPs" dataset. The characteristics of the
214 different datasets are summarised in Table S3.
215

216

Genetic differentiation and analysis of hybrids

217 We analysed the genetic structure and differentiation among species with the
218 "polymorphic sites" dataset. We used the LEA R package to estimate
219 ancestry coefficients (Frichot et al., 2014; Frichot & François, 2015). We
220 tested K values from 1 to 10, with 10 replicates for each K. To analyse the
221 genetic differences among individuals, we performed a Principal Component
222 Analysis (PCA) with the R package adegenet. The pairwise F<sub>ST</sub> among species
223 were computed with the R package Genepop (Rousset, 2008; Rousset et al.,
224 2020), after conversion of the vcf file with PGDSpider (Lischer & Excoffier,
225 2012). We present here the F<sub>ST</sub> estimates of Weir & Cockerham (1984). The
226 distribution of F<sub>ST</sub> among loci was obtained with vcftools.

227 The hybrid status (e.g. first generation hybrids) of morphologically
228 intermediate individuals was analysed with the newhybrids software
229 (Anderson & Thompson, 2002). We used the genepopedit R package to
230 prepare the input file from genepop format (Stanley et al., 2017). Following
231 the results of the LEA and PCA analyses, we focused here on the comparison
232 between *E. cavolini*, *E. singularis* and potential hybrids. The newhybrids
233 analysis had difficulties to converge with such a high number of loci
234 compared to the number of markers (see discussion here:
235 https://github.com/erigande/newhybrids/issues/5). We therefore used the "1%
236 SNP" dataset (with the 1% highest F<sub>ST</sub> between *E. cavolini* and *E. singularis*)
237 for the newhybrids analysis. As a prior, we used individuals with the lowest
238 levels of admixture in LEA as potential parental individuals: this corresponded
239 to three over five individuals for *E. cavolini*, and six over eight individuals for
240 *E. singularis* (see results). The newhybrids analysis was repeated five times
241 with different seeds to test the stability of the results.
242

243

Scenarios of speciation

244 We tested scenarios of speciation with the Demographic Inferences with Linked
245 Selection (DILS) pipeline (Csilléry et al., 2012; Fraïsse et al., 2021; Pudlo et al.,
246 2016). The DILS pipeline allows the analysis of two species scenarios only: we
247 therefore performed separate analyses for the three two-species
248 comparisons, with the "all-CS", "all-CV", and "all-SV" pairwise datasets. We
249 first tested the priors proposed by Monnet et al. (2023) according to observed
250 statistics of polymorphism and divergence: we then increased the range of
251 priors to get better results on the goodness-of-fit tests. We used the same
252 priors for all analyses, with different numbers of sequences per gene and per
253 sample according to the dataset (Table S4). For all pairwise comparisons, we
254 performed two DILS analyses: one with constant population sizes, and one
255 with variable population sizes.
256

256

Results
258 **Phylogenetic relationships**

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

266

267 **Presence of Symbiodiniaceae**

268 The analysed samples showed low numbers of reads counts for the
269 Symbiodiniaceae transcriptome (between 1868 and 58406 reads; Table S5).
270 The proportion of assembled sequences corresponding to Symbiodiniaceae
271 with BLAT was also very low (between 0.00276 and 0.03686; Table S5).
272 Significant differences were observed among species in both cases (Kruskal-
273 Wallis test, p = 0.047 for reads counts, and p = 0.002 for the proportions of
274 assembled sequences). The pairwise Wilcoxon-Test showed significant
275 differences only for the comparisons of proportions of assembled sequences
276 involving *E. singularis* (Table S6; Figure S3). The frequency of assembled
277 sequences corresponding to Symbiodiniaceae was higher in *E. singularis*
278 compared to other species: in *E. singularis*, it varied between 0.004 and
279 0.037 (mean = 0.022), whereas it reached 0.004 in one sample of *E. cavolini*
280 and it was lower than 0.004 in all other non-singularis samples. The mean
281 values of reads counts and assembled sequences in the hybrids were lower
282 than in *E. cavolini* but higher than in *E. verrucosa*; the corresponding pairwise
283 tests were not significant.

284 The blast analysis with the LSU, ITS and psbA sequences of *Philozoon* only
285 retrieved some sequences in the transcriptomes of *E. singularis*: no hit was
286 observed in the other samples. Regarding the pilot study of 18S rDNA
287 metabarcoding, a diversity of 92 Operational Taxonomic Units (OTUs)
288 corresponding to Symbiodiniaceae was observed in *E. singularis*, with a single
289 OTU largely dominant in abundance. The same OTU was also observed in E.
290 cavolini with a low abundance of reads, but still representing 99% of all 12 to
291 13 Symbiodiniaceae OTUs detected in the two analysed colonies. A blast
292 analysis of this shared OTU on Genbank indicated that it was phylogenetically
293 related to the clade A of Symbiodiniaceae (Supplementary file S2).
294

294

295 **Genetic differentiation and analysis of hybrids**

296 We obtained 31369 SNPs for the "polymorphic sites" dataset. With this
297 dataset, the highest F<sub>ST</sub> values were observed for the comparisons between
298 *E. verrucosa* and all other samples (F<sub>ST</sub> > 0.43; table S7). The F<sub>ST</sub> between *E.
299 cavolini* and *E. singularis* was much lower (0.21), and the lowest F<sub>ST</sub> values
300 were observed for hybrids compared to these two species (F<sub>ST</sub> around 0.07 in
301 both case). These differences corresponded to different distributions of F<sub>ST</sub>
302 over SNPs for the three inter-specific comparisons (Figure S4). Both
303 comparisons with *E. verrucosa* involved a more heterogeneous distribution of
304 F<sub>ST</sub>, and more SNPs reaching the maximum estimate of 1 (1364 SNPs for
305 *E. cavolini* vs *E. verrucosa*, and 1641 SNPs for *E. singularis* vs *E. verrucosa*),
306 than for the comparison between *E. cavolini* and *E. singularis* (F<sub>ST</sub> = 1 for 41
307 SNPs). If we consider the 1% SNPs with the highest F<sub>ST</sub> estimates, 52 were
308 shared by both comparisons involving *E. cavolini* (i.e. *E. cavolini* vs
309 *E. singularis* and *E. cavolini* vs *E. verrucosa*), amounting to 0.17% of all SNPs,
258 **Phylogenetic relationships**

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

266

267 **Presence of Symbiodiniaceae**

268 The analysed samples showed low numbers of reads counts for the
269 Symbiodiniaceae transcriptome (between 1868 and 58406 reads; Table S5).
270 The proportion of assembled sequences corresponding to Symbiodiniaceae
271 with BLAT was also very low (between 0.00276 and 0.03686; Table S5).
272 Significant differences were observed among species in both cases (Kruskal-
273 Wallis test, p = 0.047 for reads counts, and p = 0.002 for the proportions of
274 assembled sequences). The pairwise Wilcoxon-Test showed significant
275 differences only for the comparisons of proportions of assembled sequences
276 involving *E. singularis* (Table S6; Figure S3). The frequency of assembled
277 sequences corresponding to Symbiodiniaceae was higher in *E. singularis*
278 compared to other species: in *E. singularis*, it varied between 0.004 and
279 0.037 (mean = 0.022), whereas it reached 0.004 in one sample of *E. cavolini*
280 and it was lower than 0.004 in all other non-singularis samples. The mean
281 values of reads counts and assembled sequences in the hybrids were lower
282 than in *E. cavolini* but higher than in *E. verrucosa*; the corresponding pairwise
283 tests were not significant.

284 The blast analysis with the LSU, ITS and psbA sequences of *Philozoon* only
285 retrieved some sequences in the transcriptomes of *E. singularis*: no hit was
286 observed in the other samples. Regarding the pilot study of 18S rDNA
287 metabarcoding, a diversity of 92 Operational Taxonomic Units (OTUs)
288 corresponding to Symbiodiniaceae was observed in *E. singularis*, with a single
289 OTU largely dominant in abundance. The same OTU was also observed in E.
290 cavolini with a low abundance of reads, but still representing 99% of all 12 to
291 13 Symbiodiniaceae OTUs detected in the two analysed colonies. A blast
292 analysis of this shared OTU on Genbank indicated that it was phylogenetically
293 related to the clade A of Symbiodiniaceae (Supplementary file S2).
294

294

295 **Genetic differentiation and analysis of hybrids**

296 We obtained 31369 SNPs for the "polymorphic sites" dataset. With this
297 dataset, the highest F<sub>ST</sub> values were observed for the comparisons between
298 *E. verrucosa* and all other samples (F<sub>ST</sub> > 0.43; table S7). The F<sub>ST</sub> between *E.
299 cavolini* and *E. singularis* was much lower (0.21), and the lowest F<sub>ST</sub> values
300 were observed for hybrids compared to these two species (F<sub>ST</sub> around 0.07 in
301 both case). These differences corresponded to different distributions of F<sub>ST</sub>
302 over SNPs for the three inter-specific comparisons (Figure S4). Both
303 comparisons with *E. verrucosa* involved a more heterogeneous distribution of
304 F<sub>ST</sub>, and more SNPs reaching the maximum estimate of 1 (1364 SNPs for
305 *E. cavolini* vs *E. verrucosa*, and 1641 SNPs for *E. singularis* vs *E. verrucosa*),
306 than for the comparison between *E. cavolini* and *E. singularis* (F<sub>ST</sub> = 1 for 41
307 SNPs). If we consider the 1% SNPs with the highest F<sub>ST</sub> estimates, 52 were
308 shared by both comparisons involving *E. cavolini* (i.e. *E. cavolini* vs
309 *E. singularis* and *E. cavolini* vs *E. verrucosa*), amounting to 0.17% of all SNPs,
116 top 1% SNPs were shared by both comparisons involving *E. singularis* (0.37% of all SNPs), and 1042 top 1% SNPs were shared by both comparisons involving *E. verrucosa* (3.32% of all SNPs). All the results point to the higher differentiation of *E. verrucosa* compared to the two other species, with more outlier SNPs for the divergence of this species.

The cross-entropy analysis with LEA did not give a clear signal, with the best clustering solution corresponding to K = 2 or K = 3 clusters (Figure S5). The barplots of the coancestry coefficients for K = 2 to 4 are presented in Figure 2. At K = 2, the first distinction was observed between *E. verrucosa* and all other samples. The K = 3 analysis further separated *E. cavolini* and *E. singularis*, and the morphologically intermediate individuals appeared well 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 (site 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 their conspecifics from the northern part of the Mediterranean, without admixture. The low signal of introgression in two *E. singularis* individuals from Banyuls appeared to correspond to *E. cavolini* from the northern and the southern part, or only the southern part of the Mediterranean depending on the individuals. A very low signal of *E. verrucosa* ancestry was observed in most *E. singularis* individuals, but not in *E. cavolini*.

The PCA first separated *E. verrucosa* from other samples on the first axis (33.2% of variance; Figure S6). The second axis (13% of variance) separated *E. cavolini* and *E. singularis*, with the potential hybrids in intermediate position between them.

The newhybrids analysis indicated that the morphologically intermediate individuals, which displayed intermediate ancestry with LEA, were indeed hybrids with a probability of one in all five iterations of the analysis. One individual was a first-generation hybrid, another one was a second-generation hybrid, and the two other ones corresponded to backcrossing with *E. singularis* (Figure 2). In the same analysis, the *E. cavolini* and *E. singularis* individuals not included as priors for parental species (see Figure 2 for the individuals used as priors), were indeed inferred as parental with a probability of one.

## Scenarios of speciation

The Supplementary File S3 gives the complete results of the DILS analysis including estimated parameters distribution, comparisons between observed and expected joint Site Frequency Spectrum (jSFS), and the PCA comparing the observed dataset with the prior and posterior distributions. 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 S7). 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 1). 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 was inferred with strong
116 top 1% SNPs were shared by both comparisons involving *E. singularis* (0.37% of all SNPs), and 1042 top 1% SNPs were shared by both comparisons involving *E. verrucosa* (3.32% of all SNPs). All the results point to the higher differentiation of *E. verrucosa* compared to the two other species, with more outlier SNPs for the divergence of this species.

The cross-entropy analysis with LEA did not give a clear signal, with the best clustering solution corresponding to K = 2 or K = 3 clusters (Figure S5). The barplots of the coancestry coefficients for K = 2 to 4 are presented in Figure 2. At K = 2, the first distinction was observed between *E. verrucosa* and all other samples. The K = 3 analysis further separated *E. cavolini* and *E. singularis*, and the morphologically intermediate individuals appeared well 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 (site 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 their conspecifics from the northern part of the Mediterranean, without admixture. The low signal of introgression in two *E. singularis* individuals from Banyuls appeared to correspond to *E. cavolini* from the northern and the southern part, or only the southern part of the Mediterranean depending on the individuals. A very low signal of *E. verrucosa* ancestry was observed in most *E. singularis* individuals, but not in *E. cavolini*.

The PCA first separated *E. verrucosa* from other samples on the first axis (33.2% of variance; Figure S6). The second axis (13% of variance) separated *E. cavolini* and *E. singularis*, with the potential hybrids in intermediate position between them.

The newhybrids analysis indicated that the morphologically intermediate individuals, which displayed intermediate ancestry with LEA, were indeed hybrids with a probability of one in all five iterations of the analysis. One individual was a first-generation hybrid, another one was a second-generation hybrid, and the two other ones corresponded to backcrossing with *E. singularis* (Figure 2). In the same analysis, the *E. cavolini* and *E. singularis* individuals not included as priors for parental species (see Figure 2 for the individuals used as priors), were indeed inferred as parental with a probability of one.

## Scenarios of speciation

The Supplementary File S3 gives the complete results of the DILS analysis including estimated parameters distribution, comparisons between observed and expected joint Site Frequency Spectrum (jSFS), and the PCA comparing the observed dataset with the prior and posterior distributions. 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 S7). 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 1). 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 was inferred with strong
support (p ≥ 0.99) for all analyses. In the case of current migration (between *E. cavolini* and *E. singularis*), a genomic heterogeneity in migration rates was inferred (p ≥ 0.82). The inferred parameters for the different scenarios are presented in Supplementary Table S8. 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 at 62 039 generations (median estimate), which translates in around 85% of time 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. Regarding effective size, the comparisons all indicated lower values for *E. singularis* compared to both other species (around 200 000 to 300 000 for *E. singularis*, around 600 000 to 750 000 for both other species). The estimated effective sizes were of similar order for *E. cavolini* and *E. verrucosa*. This aligns well with the observed estimates of nucleotide diversity around 0.007-0.009 for *E. cavolini*, 0.005-0.006 for *E. singularis*, and 0.007-0.009 for *E. verrucosa* (variations depend on the comparison which may change the retained dataset; Supplementary File S3). 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

## Species relationships and differentiation

As already observed with COI (Aurelle et al., 2017), the mitochondrial MutS marker did not allow to discriminate the three species. This is the consequence of the usually slow evolution of mitochondrial DNA in octocorals (Muthye et al., 2022). The use of transcriptome sequences first confirmed that *E. cavolini* and *E. singularis* are sister species, with *E. verrucosa* being more distantly related. The $F_{ST}$ estimate between *E. verrucosa* and the two other species was more than twice higher than the estimate between *E. cavolini* and *E. singularis*. This had been previously suggested with two intron sequences, but incomplete lineage sorting for these markers did not allow a formal conclusion (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 shown in other genera by the use of RAD sequencing or ultra conserved elements and exons (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. One can note that counter examples exist which show that the presence of this locus is not the only factor explaining the slow evolution of mitochondrial DNA in octocorals (Muthye et al., 2022).
support (p ≥ 0.99) for all analyses. In the case of current migration (between *E. cavolini* and *E. singularis*), a genomic heterogeneity in migration rates was inferred (p ≥ 0.82). The inferred parameters for the different scenarios are presented in Supplementary Table S8. 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 at 62 039 generations (median estimate), which translates in around 85% of time 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. Regarding effective size, the comparisons all indicated lower values for *E. singularis* compared to both other species (around 200 000 to 300 000 for *E. singularis*, around 600 000 to 750 000 for both other species). The estimated effective sizes were of similar order for *E. cavolini* and *E. verrucosa*. This aligns well with the observed estimates of nucleotide diversity around 0.007-0.009 for *E. cavolini*, 0.005-0.006 for *E. singularis*, and 0.007-0.009 for *E. verrucosa* (variations depend on the comparison which may change the retained dataset; Supplementary File S3). 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

## Species relationships and differentiation

As already observed with COI (Aurelle et al., 2017), the mitochondrial MutS marker did not allow to discriminate the three species. This is the consequence of the usually slow evolution of mitochondrial DNA in octocorals (Muthye et al., 2022). The use of transcriptome sequences first confirmed that *E. cavolini* and *E. singularis* are sister species, with *E. verrucosa* being more distantly related. The $F_{ST}$ estimate between *E. verrucosa* and the two other species was more than twice higher than the estimate between *E. cavolini* and *E. singularis*. This had been previously suggested with two intron sequences, but incomplete lineage sorting for these markers did not allow a formal conclusion (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 shown in other genera by the use of RAD sequencing or ultra conserved elements and exons (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. One can note that counter examples exist which show that the presence of this locus is not the only factor explaining the slow evolution of mitochondrial DNA in octocorals (Muthye et al., 2022).
The divergence between *E. singularis* and *E. cavolini* would then have occurred after the divergence from *E. verrucosa*. One can note the highest divergence of *E. singularis* from *E. verrucosa* compared to *E. cavolini* (pairwise FST 0.529 and 0.432 respectively): this could point to a more complex evolutionary history than a single divergence, for example involving different demographic histories, or gene flow from unsampled taxa (Tricou et al., 2022).

## Speciation scenarios

The speciation scenarios inferred from DILS are in agreement with the observed pattern of differentiation. For the results obtained here with pairwise analyses, a much higher divergence time was observed for both comparisons involving *E. verrucosa* than for those involving other species. A scenario of constant isolation was strongly supported in both cases, with the possibility of ancestral migration. The current isolation of *E. verrucosa* from *E. cavolini* is 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 interesting to explore the reasons for the discrepancy between this study and the present one, which are both based on transcriptome datasets but obtained from different samples and sequencing platforms.

*Eunicella verrucosa* is currently more widely distributed in the North Eastern Atlantic Ocean, and less frequent in the Mediterranean Sea, whereas both other species analyzed here 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). One could then envision 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 studied here, 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 suppose a lower hypothesis of 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). As a comparison, this is of the same order as estimates of divergence time for sister scleractinian species which can go from around 1 Ma to 4-5 Ma (Aurelle et al., 2024; Johnston et al., 2017). Note that higher generation time would point to divergence time older than 5 Ma, where the Messinian crisis (between 5.5 and 6 Ma; Rouchy & Caruso, 2006) could have played a role in initiating divergence. Obviously, more data on generation time and mutation rate are needed to tentatively link speciation times with fluctuations in paleo-environments. With a median estimate around 400 000 generations, the divergence time between *E. cavolini* and *E. singularis* would be 2.5 times more recent, around 800 000 years. The median time of secondary contact between these two species would be 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 (such as near Marseille), their range do not completely overlap. For example *E. cavolini* is
The divergence between *E. singularis* and *E. cavolini* would then have occurred after the divergence from *E. verrucosa*. One can note the highest divergence of *E. singularis* from *E. verrucosa* compared to *E. cavolini* (pairwise FST 0.529 and 0.432 respectively): this could point to a more complex evolutionary history than a single divergence, for example involving different demographic histories, or gene flow from unsampled taxa (Tricou et al., 2022).

420

## Speciation scenarios

The speciation scenarios inferred from DILS are in agreement with the observed pattern of differentiation. For the results obtained here with pairwise analyses, a much higher divergence time was observed for both comparisons involving *E. verrucosa* than for those involving other species. A scenario of constant isolation was strongly supported in both cases, with the possibility of ancestral migration. The current isolation of *E. verrucosa* from *E. cavolini* is 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 interesting to explore the reasons for the discrepancy between this study and the present one, which are both based on transcriptome datasets but obtained from different samples and sequencing platforms.

*Eunicella verrucosa* is currently more widely distributed in the North Eastern Atlantic Ocean, and less frequent in the Mediterranean Sea, whereas both other species analyzed here 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). One could then envision 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 studied here, 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 suppose a lower hypothesis of 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). As a comparison, this is of the same order as estimates of divergence time for sister scleractinian species which can go from around 1 Ma to 4-5 Ma (Aurelle et al., 2024; Johnston et al., 2017). Note that higher generation time would point to divergence time older than 5 Ma, where the Messinian crisis (between 5.5 and 6 Ma; Rouchy & Caruso, 2006) could have played a role in initiating divergence. Obviously, more data on generation time and mutation rate are needed to tentatively link speciation times with fluctuations in paleo-environments. With a median estimate around 400 000 generations, the divergence time between *E. cavolini* and *E. singularis* would be 2.5 times more recent, around 800 000 years. The median time of secondary contact between these two species would be 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 (such as near Marseille), their range do not completely overlap. For example *E. cavolini* is
nearly absent at the West of the Rhone on French coasts, 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* (Carugati et al., 2022; Gori et al., 2012). 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.

### 473 Hybridisation and species barriers

The analysis of genetic ancestry with LEA and the hybrid inferences both confirmed that morphologically intermediate individuals are indeed hybrids between *E. singularis* and *E. cavolini*, with the identification of one F1, one F2 and one backcross with *E. singularis*: the last two categories indicate that first generation hybrids can be fertile and can participate in reproduction. The possibility of current gene flow is confirmed with the DILS analysis, with a probability of 0.87 and 0.88 for models with constant and variable population sizes respectively. Importantly the DILS analysis did not include hybrid individuals, which indicates that gene flow between these two species indeed goes further than the aforementioned hybrid levels. Accordingly, the LEA analysis also indicated a small level of *E. cavolini* ancestry in some *E. singularis* individuals, but the reverse was not observed: this would point to asymmetric gene flow, more frequent towards *E. singularis*, but this does not agree with the results of DILS. Quantifying hybrid frequency in situ may be difficult, as we don't know the range of morphological variations in hybrids. Nevertheless, the ease to find such hybrids in the area studied here, as well as similar observations in other sites near Marseille (S. Sartoretto, pers. com.) indicates that hybridization is not rare at an evolutionary scale.

The alternation of parental populations with mixed populations would point to a mosaic hybrid zone (Bierne et al., 2003), where hybrids could form in different areas and from different genetic compositions of parental species. As hybridization between *E. cavolini* and *E. singularis* had not been reported before, the presence of hybrids has probably been overlooked up to now. Communicating on this subject towards scientists and diving associations might help sampling other potential hybrid zones: this would be useful to analyse the spatial and ecological distribution of hybrids. One interesting question in this context is whether changes in selection regimes induced by human activities can change the outcome of hybridization (Ålund et al., 2023). For example, Mediterranean octocorals are impacted by mortality events linked with climate change (Estaque et al., 2023; Sini et al., 2015), and it would be interesting to compare the thermotolerance of hybrids and parental individuals. In scleractinian corals, interspecific hybridisation has been reported to enhance the survival under elevated temperature conditions (Chan et al., 2018).

Regarding the comparison between *E. cavolini* and *E. singularis*, the DILS analysis pointed to an homogeneity of gene flow. The net divergence between *E. cavolini* and *E. singularis* (0.0018) puts these two species below the grey zone of speciation in the analyses of Roux et al. (2016). They may correspond to an early stage of speciation where genomic islands of differentiation did not expand in the genome of these species (see Peñalba et al., 2024, and references therein). One can note that in the meta-analysis evolutionary inferences of De Jode et al. (2023), an heterogeneity of gene flow was more frequently inferred than homogeneity for scenarios of
nearly absent at the West of the Rhone on French coasts, 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* (Carugati et al., 2022; Gori et al., 2012). 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.

472

### 473 Hybridisation and species barriers

The analysis of genetic ancestry with LEA and the hybrid inferences both confirmed that morphologically intermediate individuals are indeed hybrids between *E. singularis* and *E. cavolini*, with the identification of one F1, one F2 and one backcross with *E. singularis*: the last two categories indicate that first generation hybrids can be fertile and can participate in reproduction. The possibility of current gene flow is confirmed with the DILS analysis, with a probability of 0.87 and 0.88 for models with constant and variable population sizes respectively. Importantly the DILS analysis did not include hybrid individuals, which indicates that gene flow between these two species indeed goes further than the aforementioned hybrid levels. Accordingly, the LEA analysis also indicated a small level of *E. cavolini* ancestry in some *E. singularis* individuals, but the reverse was not observed: this would point to asymmetric gene flow, more frequent towards *E. singularis*, but this does not agree with the results of DILS. Quantifying hybrid frequency in situ may be difficult, as we don't know the range of morphological variations in hybrids. Nevertheless, the ease to find such hybrids in the area studied here, as well as similar observations in other sites near Marseille (S. Sartoretto, pers. com.) indicates that hybridization is not rare at an evolutionary scale.

The alternation of parental populations with mixed populations would point to a mosaic hybrid zone (Bierne et al., 2003), where hybrids could form in different areas and from different genetic compositions of parental species. As hybridization between *E. cavolini* and *E. singularis* had not been reported before, the presence of hybrids has probably been overlooked up to now. Communicating on this subject towards scientists and diving associations might help sampling other potential hybrid zones: this would be useful to analyse the spatial and ecological distribution of hybrids. One interesting question in this context is whether changes in selection regimes induced by human activities can change the outcome of hybridization (Ålund et al., 2023). For example, Mediterranean octocorals are impacted by mortality events linked with climate change (Estaque et al., 2023; Sini et al., 2015), and it would be interesting to compare the thermotolerance of hybrids and parental individuals. In scleractinian corals, interspecific hybridisation has been reported to enhance the survival under elevated temperature conditions (Chan et al., 2018).

Regarding the comparison between *E. cavolini* and *E. singularis*, the DILS analysis pointed to an homogeneity of gene flow. The net divergence between *E. cavolini* and *E. singularis* (0.0018) puts these two species below the grey zone of speciation in the analyses of Roux et al. (2016). They may correspond to an early stage of speciation where genomic islands of differentiation did not expand in the genome of these species (see Peñalba et al., 2024, and references therein). One can note that in the meta-analysis evolutionary inferences of De Jode et al. (2023), an heterogeneity of gene flow was more frequently inferred than homogeneity for scenarios of
secondary contact, which points to a diversity of stages of divergence in the
tested species.

The distribution of F<sub>ST</sub> between *E. cavolini* and *E. singularis*, indicated a slight heterogeneity, with few loci with F<sub>ST</sub> higher than 0.5 for example. Conversely, both comparisons with *E. verrucosa* showed much more loci with high F<sub>ST</sub>, as expected with higher divergence time. Among the list of the most highly differentiated loci, more overlap was also observed for the two comparisons involving *E. verrucosa* than for the other pairwise comparisons: this may indicate that few areas of potential incompatibilities with *E. verrucosa* are involved in the divergence between *E. cavolini* and *E. singularis*. Such incompatibilities could explain the persistence of interspecific differentiation despite a quite long time since secondary contact. These conclusions should be investigated more thoroughly with genome wide analysis of differentiation, and by taking into account the alternative hypotheses explaining genomic islands of divergence such as background selection or local adaptation (Quilodrán et al., 2020).

533

534 Evolution of symbiosis

As previously discussed, we clearly demonstrated here the possibility of gene
flow between a symbiotic (i.e. hosting Symbiodiniaceae) and a non-symbiotic
octocoral species: the hybrid status then does not lead to complete isolation,
even if Symbiodiniaceae could be involved in genetic incompatibilities. The
methods used here did not aim at a precise quantification of
Symbiodiniaceae, and one can note the general low levels of sequences
corresponding to these symbionts, even in *E. singularis*: this may be due to
difficulties in extracting the RNA of the symbionts with this protocol or to loss
of symbionts before extractions. This is not a general result: for example, a
whole transcriptome sequencing of the octocorals *Heliopora coerulea* and
*Briareum asbestinum* recovered 29% and 17.2% of Symbiodiniaceae
sequences respectively (Guzman et al., 2018; Rivera-García et al., 2019).
Despite these limits we observed a clear signal of 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 tested with more samples. These results indicate a breakdown of
symbiosis following hybridization with potential consequences on the fitness
of hybrids. In *E. singularis*, the transmission of Symbiodiniaceae seems to
occur both vertically, through ovules, and horizontally, from the environment
(Forcioli et al., 2011). Vertical transmission may change the fate of hybrids
depending on the species of the mother, as it can change the initial load in
Symbiodiniaceae: this hypothesis is nevertheless difficult to test as these
species can not be reproduced in aquarium. More generally, the question of
the link between symbiosis and the fitness of hybrids would require a
dedicated study involving more precise estimates of Symbiodiniaceae
abundance (e.g. with quantitative PCR), and of physiological parameters such
as photosynthetic and respiration rates in controlled conditions (Ezzat et al.,
2013).

Our results also question the evolution and significance of octocoral /
Symbiodiniaceae symbiosis. A phylogenetic study of scleractinians has shown
that in these hexacorals, the transition between symbiotic and non-symbiotic
states happened repeatedly, but mostly in the direction of the acquisition of
symbiosis with Symbiodiniaceae, with very low rates of transition in the
517 secondary contact, which points to a diversity of stages of divergence in the
tested species.
518

519 The distribution of F<sub>ST</sub> between *E. cavolini* and *E. singularis*, indicated a slight
520 heterogeneity, with few loci with F<sub>ST</sub> higher than 0.5 for example. Conversely,
521 both comparisons with *E. verrucosa* showed much more loci with high F<sub>ST</sub>, as
522 expected with higher divergence time. Among the list of the most highly
differentiated loci, more overlap was also observed for the two comparisons
523 involving *E. verrucosa* than for the other pairwise comparisons: this may
524 indicate that few areas of potential incompatibilities with *E. verrucosa* are
525 involved in the divergence between *E. cavolini* and *E. singularis*. Such
526 incompatibilities could explain the persistence of interspecific differentiation
527 despite a quite long time since secondary contact. These conclusions should
528 be investigated more thoroughly with genome wide analysis of
529 differentiation, and by taking into account the alternative hypotheses
530 explaining genomic islands of divergence such as background selection or
531 local adaptation (Quilodrán et al., 2020).
532

533

534 Evolution of symbiosis

535 As previously discussed, we clearly demonstrated here the possibility of gene
536 flow between a symbiotic (i.e. hosting Symbiodiniaceae) and a non-symbiotic
537 octocoral species: the hybrid status then does not lead to complete isolation,
538 even if Symbiodiniaceae could be involved in genetic incompatibilities. The
539 methods used here did not aim at a precise quantification of
540 Symbiodiniaceae, and one can note the general low levels of sequences
541 corresponding to these symbionts, even in *E. singularis*: this may be due to
542 difficulties in extracting the RNA of the symbionts with this protocol or to loss
543 of symbionts before extractions. This is not a general result: for example, a
544 whole transcriptome sequencing of the octocorals *Heliopora coerulea* and
545 *Briareum asbestinum* recovered 29% and 17.2% of Symbiodiniaceae
546 sequences respectively (Guzman et al., 2018; Rivera-García et al., 2019).
547 Despite these limits we observed a clear signal of higher Symbiodiniaceae
548 concentration in *E. singularis* than in *E. cavolini* and *E. verrucosa*.
549 Interestingly, the hybrids showed a lower frequency of Symbiodiniaceae than
550 *E. singularis*, and possibly than *E. cavolini*, though this last result remains to
551 be tested with more samples. These results indicate a breakdown of
552 symbiosis following hybridization with potential consequences on the fitness
553 of hybrids. In *E. singularis*, the transmission of Symbiodiniaceae seems to
554 occur both vertically, through ovules, and horizontally, from the environment
555 (Forcioli et al., 2011). Vertical transmission may change the fate of hybrids
depending on the species of the mother, as it can change the initial load in
556 Symbiodiniaceae: this hypothesis is nevertheless difficult to test as these
557 species can not be reproduced in aquarium. More generally, the question of
558 the link between symbiosis and the fitness of hybrids would require a
559 dedicated study involving more precise estimates of Symbiodiniaceae
560 abundance (e.g. with quantitative PCR), and of physiological parameters such
561 as photosynthetic and respiration rates in controlled conditions (Ezzat et al.,
562 2013).
563

564 Our results also question the evolution and significance of octocoral /
565 Symbiodiniaceae symbiosis. A phylogenetic study of scleractinians has shown
566 that in these hexacorals, the transition between symbiotic and non-symbiotic
567 states happened repeatedly, but mostly in the direction of the acquisition of
568 symbiosis with Symbiodiniaceae, with very low rates of transition in the
reverse direction (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 Mediterranean Sea, all octocoral species are non-symbiotic, except for *E. singularis*. The most parsimonious scenario here would be an acquisition of symbiosis in *E. singularis* during or following its divergence from *E. cavolini*. The symbiotic status of *E. singularis* nevertheless could be facultative as non-symbiotic colonies of *E. singularis* have been observed between 40 and 60 m depths (Gori et al., 2012). Additionally, experimental physiological studies have demonstrated the nutritional plasticity of *E. singularis* which is able to use either heterotrophy or autotrophy for its metabolism (Ezzat et al., 2013). Nevertheless, in natural conditions, autotrophy seems to provide an important contribution to the metabolism of *E. singularis*, and the collapse of photosynthetic capacities in too warm conditions could contribute to mortality events in this species (Coma et al., 2015). The question of symbiosis could be reversed as well: why don't we observe any symbiosis with Symbiodiniaceae in *E. cavolini*? This species can be observed in shallow conditions (less than 10 m depth) where there is enough light for photosynthesis, and in syntopy with *E. singularis*. The availability of preys or particulate organic matter may provide enough energy to *E. cavolini* in its habitat, but this species may have never engaged in mutualistic interaction with Symbiodiniaceae. Interestingly we observed a low rate of sequences related to Symbiodiniaceae in the transcriptomes of *E. cavolini* (and even lower, but not null in *E. verrucosa*). This could either correspond to a signal from free living Symbiodiniaceae, or to background, transient, associations with the cnidarian. In addition, a Symbiodiniaceae OTU that is common to *E. singularis* and *E. cavolini* was identified among the microeukaryotes associated with the two species, which probably corresponds to the Symbiodiniaceae species symbiotic with *E. singularis*, and is related to strains observed in symbiosis with other cnidarians. Rare Symbiodiniaceae strains can also be observed in symbiotic hexacorals, probably with low impact on the hosts physiology (Lee et al., 2016). Molecular markers also allowed to evidence the presence of Symbiodiniaceae in species previously supposed to be asymbiotic, as in the Mediterranean octocoral *Paramuricea clavata*, and in several Hawaiian antipatharian species (Bonacolta et al., 2024; Wagner et al., 2011). These results, and our observations in *Eunicella* species, obviously underline the dynamic nature of interactions between Symbiodiniaceae and cnidarians: the establishment of these symbiotic interactions may be preceded by various types of interactions, more or less stable, and more or less mutualistic. The development of effective symbiosis, with stable relationships, and higher abundance of symbiont, would require specific adaptation from both partners. We can see here that even if on a macro-evolutionary scale, the acquisition of symbiosis is much more frequent than its loss, on a micro-evolutionary scale the gene flow between the *Eunicella* species analysed here did not lead to the full development of symbiosis in *E. cavolini*.

Conclusions and perspectives
reverse direction (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 Mediterranean Sea, all octocoral species are non-symbiotic, except for *E. singularis*. The most parsimonious scenario here would be an acquisition of symbiosis in *E. singularis* during or following its divergence from *E. cavolini*. The symbiotic status of *E. singularis* nevertheless could be facultative as non-symbiotic colonies of *E. singularis* have been observed between 40 and 60 m depths (Gori et al., 2012). Additionally, experimental physiological studies have demonstrated the nutritional plasticity of *E. singularis* which is able to use either heterotrophy or autotrophy for its metabolism (Ezzat et al., 2013). Nevertheless, in natural conditions, autotrophy seems to provide an important contribution to the metabolism of *E. singularis*, and the collapse of photosynthetic capacities in too warm conditions could contribute to mortality events in this species (Coma et al., 2015). The question of symbiosis could be reversed as well: why don't we observe any symbiosis with Symbiodiniaceae in *E. cavolini*? This species can be observed in shallow conditions (less than 10 m depth) where there is enough light for photosynthesis, and in syntopy with *E. singularis*. The availability of preys or particulate organic matter may provide enough energy to *E. cavolini* in its habitat, but this species may have never engaged in mutualistic interaction with Symbiodiniaceae. Interestingly we observed a low rate of sequences related to Symbiodiniaceae in the transcriptomes of *E. cavolini* (and even lower, but not null in *E. verrucosa*). This could either correspond to a signal from free living Symbiodiniaceae, or to background, transient, associations with the cnidarian. In addition, a Symbiodiniaceae OTU that is common to *E. singularis* and *E. cavolini* was identified among the microeukaryotes associated with the two species, which probably corresponds to the Symbiodiniaceae species symbiotic with *E. singularis*, and is related to strains observed in symbiosis with other cnidarians. Rare Symbiodiniaceae strains can also be observed in symbiotic hexacorals, probably with low impact on the hosts physiology (Lee et al., 2016). Molecular markers also allowed to evidence the presence of Symbiodiniaceae in species previously supposed to be asymbiotic, as in the Mediterranean octocoral *Paramuricea clavata*, and in several Hawaiian antipatharian species (Bonacolta et al., 2024; Wagner et al., 2011). These results, and our observations in *Eunicella* species, obviously underline the dynamic nature of interactions between Symbiodiniaceae and cnidarians: the establishment of these symbiotic interactions may be preceded by various types of interactions, more or less stable, and more or less mutualistic. The development of effective symbiosis, with stable relationships, and higher abundance of symbiont, would require specific adaptation from both partners. We can see here that even if on a macro-evolutionary scale, the acquisition of symbiosis is much more frequent than its loss, on a micro-evolutionary scale the gene flow between the *Eunicella* species analysed here did not lead to the full development of symbiosis in *E. cavolini*.

Conclusions and perspectives
We here demonstrated the semi-permeable nature of species barriers between two octocoral species with and without mutualistic interaction with Symbiodiniaceae. This opens the way to further researches. On these model species it will be useful to estimate the frequency and spatial extent of hybrid zone: 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 search for islands of divergence, and to look for the impact on 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 between species of symbiotic scleractinian corals (e.g. Combosch & Vollmer, 2015): it would then be interesting to study the dynamics of symbiosis in these cases when different Symbiodiniaceae species are involved.

## Acknowledgements:

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. Camille Roux, Jonathan Romiguier and Christelle Fraïsse were of a great help for the analysis of speciation scenarios. We thank the SCUBA diver team of INSU/OSU Pytheas (Dorian Guillemain) for their help during the fieldwork. 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 for his 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). Part of this work has been performed during a CNRS detachment position of D. Aurelle at the ISYEB laboratory.

## Data availability
We here demonstrated the semi-permeable nature of species barriers between two octocoral species with and without mutualistic interaction with Symbiodiniaceae. This opens the way to further researches. On these model species it will be useful to estimate the frequency and spatial extent of hybrid zone: 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 search for islands of divergence, and to look for the impact on 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 between species of symbiotic scleractinian corals (e.g. Combosch & Vollmer, 2015): it would then be interesting to study the dynamics of symbiosis in these cases when different Symbiodiniaceae species are involved.

## Acknowledgements:

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. Camille Roux, Jonathan Romiguier and Christelle Fraïsse were of a great help for the analysis of speciation scenarios. We thank the SCUBA diver team of INSU/OSU Pytheas (Dorian Guillemain) for their help during the fieldwork. 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 for his 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). Part of this work has been performed during a CNRS detachment position of D. Aurelle at the ISYEB laboratory.

## Data availability
670 The raw sequences are available in Genbank under BioProject ID
671 PRJNA1037721 and SRA accession number SRA PRJNA1037721.

672 The scripts used in this study are available at
https://doi.org/10.5281/zenodo.10966625

674

675 Conflict of interest disclosure

676 The authors declare that they comply with the PCI rule of having no financial
677 conflicts of interest in relation to the content of the article

678

679

680 References

681

682

Ålund, M., Cenzer, M., Bierne, N., Boughman, J. W., Cerca, J., Comerford, M. S., Culicchi, A., Langerhans, B., McFarlane, S. E., & Möst, M. H. (2023). Anthropogenic Change and the Process of Speciation. *Cold Spring Harbor Perspectives in Biology*, *15*(12), a041455.

Anderson, E., & Thompson, E. A. (2002). A model-based method for identifying species hybrids using multilocus genetic data. *Genetics*, **160**(3), 1217-1229.

Arnold, M. L., Bulger, M. R., Burke, J. M., Hempel, A. L., & Williams, J. H. (1999). Natural hybridization: How low can you go and still be important? *Ecology*, **80**(2), 371-381.

Aurelle, D., Haguenauer, A., Blaise, C., Reynes, L., Arnaud-Haond, S., Boavida, J., Cabau, C., Klopp, C., Lundalv, T., Noûs, C., Sartoretto, S., Wienberg, C., Jiménez, C. E., & Orejas, C. (2024). On the specific status of eastern Mediterranean Dendrophyllia corals (Cnidaria, Anthozoa): Genetic characterization and speciation scenarios. *Zoologica Scripta*, **53**(2), 235-247. https://doi.org/10.1111/zsc.12643

Aurelle, D., Pivotto, I. D., Malfant, M., Topcu, N. E., Masmoudi, M. B., Chaoui, L., Kara, H. M., Coelho, M. A., Castilho, R., & Haguenauer, A. (2017). Fuzzy species limits in Mediterranean gorgonians (Cnidaria, Octocorallia): Inferences on speciation processes. *Zoologica Scripta*, **46**(6), 767-778.

Baumgarten, S., Bayer, T., Aranda, M., Liew, Y. J., Carr, A., Micklem, G., & Voolstra, C. R. (2013). Integrating microRNA and mRNA expression profiling in Symbiodinium microadriaticum, a dinoflagellate symbiont of reef-building corals. *BMC Genomics*, **14**(1), 704. https://doi.org/10.1186/1471-2164-14-704

Berkelmans, R., & van Oppen, M. J. H. (2006). 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*, **273**(1599), 2305-2312.

Bierne, N., Borsa, P., Daguin, C., Jollivet, D., Viard, F., Bonhomme, F., & David, P. (2003). Introgression patterns in the mosaic hybrid zone between Mytilus edulis and M. galloprovincialis. *Molecular Ecology*, **12**, 447-462.

Bierne, N., Welch, J., Loire, E., Bonhomme, F., & David, P. (2011). The coupling hypothesis: Why genome scans may fail to map local adaptation genes. *Molecular Ecology*, **20**(10), 2044-2072.

Bonacolta, A. M., Miravall, J., Gómez-Gras, D., Ledoux, J., López-Sendino, P., Garrabou, J., Massana, R., & Del Campo, J. (2024). Differential apicomplexan presence predicts thermal stress mortality in the Mediterranean coral Paramuricea clavata. *Environmental Microbiology*, **26**(1), e16548.
670 The raw sequences are available in Genbank under BioProject ID
671 PRJNA1037721 and SRA accession number SRA PRJNA1037721.

672 The scripts used in this study are available at
https://doi.org/10.5281/zenodo.10966625

674

675 Conflict of interest disclosure

676 The authors declare that they comply with the PCI rule of having no financial
677 conflicts of interest in relation to the content of the article

678

679

680 References

681

682

Ålund, M., Cenzer, M., Bierne, N., Boughman, J. W., Cerca, J., Comerford, M. S., Culicchi, A., Langerhans, B., McFarlane, S. E., & Möst, M. H. (2023). Anthropogenic Change and the Process of Speciation. *Cold Spring Harbor Perspectives in Biology*, *15*(12), a041455.

Anderson, E., & Thompson, E. A. (2002). A model-based method for identifying species hybrids using multilocus genetic data. *Genetics*, **160**(3), 1217-1229.

Arnold, M. L., Bulger, M. R., Burke, J. M., Hempel, A. L., & Williams, J. H. (1999). Natural hybridization: How low can you go and still be important? *Ecology*, **80**(2), 371-381.

Aurelle, D., Haguenauer, A., Blaise, C., Reynes, L., Arnaud-Haond, S., Boavida, J., Cabau, C., Klopp, C., Lundalv, T., Noûs, C., Sartoretto, S., Wienberg, C., Jiménez, C. E., & Orejas, C. (2024). On the specific status of eastern Mediterranean Dendrophyllia corals (Cnidaria, Anthozoa): Genetic characterization and speciation scenarios. *Zoologica Scripta*, **53**(2), 235-247. https://doi.org/10.1111/zsc.12643

Aurelle, D., Pivotto, I. D., Malfant, M., Topcu, N. E., Masmoudi, M. B., Chaoui, L., Kara, H. M., Coelho, M. A., Castilho, R., & Haguenauer, A. (2017). Fuzzy species limits in Mediterranean gorgonians (Cnidaria, Octocorallia): Inferences on speciation processes. *Zoologica Scripta*, **46**(6), 767-778.

Baumgarten, S., Bayer, T., Aranda, M., Liew, Y. J., Carr, A., Micklem, G., & Voolstra, C. R. (2013). Integrating microRNA and mRNA expression profiling in Symbiodinium microadriaticum, a dinoflagellate symbiont of reef-building corals. *BMC Genomics*, **14**(1), 704. https://doi.org/10.1186/1471-2164-14-704

Berkelmans, R., & van Oppen, M. J. H. (2006). 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*, **273**(1599), 2305-2312.

Bierne, N., Borsa, P., Daguin, C., Jollivet, D., Viard, F., Bonhomme, F., & David, P. (2003). Introgression patterns in the mosaic hybrid zone between Mytilus edulis and M. galloprovincialis. *Molecular Ecology*, **12**, 447-462.

Bierne, N., Welch, J., Loire, E., Bonhomme, F., & David, P. (2011). The coupling hypothesis: Why genome scans may fail to map local adaptation genes. *Molecular Ecology*, **20**(10), 2044-2072.

Bonacolta, A. M., Miravall, J., Gómez-Gras, D., Ledoux, J., López-Sendino, P., Garrabou, J., Massana, R., & Del Campo, J. (2024). Differential apicomplexan presence predicts thermal stress mortality in the Mediterranean coral Paramuricea clavata. *Environmental Microbiology*, **26**(1), e16548.
Bordenstein, S. (2003). Symbiosis And The Origin Of Species. In *Insect Symbiosis* (pp. 283–303). CRC Press.

Brener-Raffalli, K., Vidal-Dupiol, J., Adjeroud, M., Rey, O., Romans, P., Bonhomme, F., Pratlong, M., Haguenauer, A., Pillot, R., Feuillassier, L., Claereboudt, M., Magalon, H., Gélin, P., Pontarotti, P., Aurelle, D., Mitta, G., & Toulza, E. (2022). Gene expression plasticity and frontloading promote thermotolerance in *Pocillopora* corals. *Peer Community Journal, 2*. https://doi.org/10.24072/pcjournal.79

Brucker, R. M., & Bordenstein, S. R. (2012). Review: Speciation by symbiosis. *Trends in Ecology & Evolution, 27*(8), 443. ScienceDirect. https://doi.org/10.1016/j.tree.2012.03.011

Cabau, C., Escudié, F., Djari, A., Guiguen, Y., Bobe, J., & Klopp, C. (2017). Compacting and correcting Trinity and Oases RNA-Seq de novo assemblies. *PeerJ*, 5, e2988.

Cahill, A. E., Meglécz, E., & Chenuil, A. (2023). Scientific history, biogeography, and biological traits predict presence of cryptic or overlooked species. *Biological Reviews, 99*, 546-561. https://doi.org/10.1111/brv.13034

Cairns, S. D. (2007). Deep-water corals: An overview with special reference to diversity and distribution of deep-water scleractinian corals. *Bulletin of Marine Science, 81*(3), 311-322.

Calderon, I., Garrabou, J., & Aurelle, D. (2006). 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, 336*, 184-197.

Campoy, A. N., Addamo, A. M., Machordom, A., Meade, A., Rivadeneira, M. M., Hernández, C. E., & Venditti, C. (2020). The Origin and Correlated Evolution of Symbiosis and Coloniality in Scleractinian Corals. *Frontiers in Marine Science, 7*.

Canessa, M., Bavestrello, G., Bo, M., Enrichetti, F., & Trainito, E. (2022). Filling a gap: A population of Eunicella verrucosa (pallas, 1766)(Anthozoa, Alcyonacea) in the Tavolara-punta Coda Cavallo marine protected area (Ne Sardinia, Italy). *Diversity, 14*(5), 405.

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

Carugati, L., Moccia, D., Bramanti, L., Cannas, R., Follesa, M. C., Salvadori, S., & Cau, A. (2022). Deep-Dwelling Populations of Mediterranean Corallium rubrum and Eunicella cavolini: Distribution, Demography, and Co-Occurrence. *Biology, 11*(2). https://doi.org/10.3390/biology11020333

Casado-Amezúa, P., Terrón-Sigler, A., Pinzón, J. H., Furla, P., Forcioli, D., Allemand, D., Ribes, M., & Coma, R. (2016). General ecological aspects of Anthozoan-Symbiodinium interactions in the Mediterranean Sea. *The Cnidaria, Past, Present and Future: The World of Medusa and Her Sisters, 375-386*.

Chan, W. Y., Peplow, L. M., Menéndez, P., Hoffmann, A. A., & Van Oppen, M. J. (2018). Interspecific hybridization may provide novel opportunities for coral reef restoration. *Frontiers in Marine Science, 5*, 160.

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

Coelho, M., Pearson, G., Boavida, J., Paulo, D., Aurelle, D., Arnaud-haond, S., Gómez-gras, D., Bensoussan, N., López-sendino, P., Cerrano, C., Kipson, S., Bakran-petricioli, T., Ferretti, E., Linares, C., Garrabou, J., Serrão, E., & Ledoux, J. (2023). Not out of the Mediterranean: Atlantic populations of the gorgonian Paramuricea clavata are a separate sister species under further
Bordenstein, S. (2003). Symbiosis And The Origin Of Species. In *Insect Symbiosis*
(pp. 283–303). CRC Press.

Brener-Raffalli, K., Vidal-Dupiol, J., Adjeroud, M., Rey, O., Romans, P., Bonhomme, F., Pratlong, M., Haguenauer, A., Pillot, R., Feuillassier, L., Claereboudt, M., Magalon, H., Gélin, P., Pontarotti, P., Aurelle, D., Mitta, G., & Toulza, E. (2022). Gene expression plasticity and frontloading promote thermotolerance in *Pocillopora* corals. *Peer Community Journal, 2*. https://doi.org/10.24072/pcjournal.79

Brucker, R. M., & Bordenstein, S. R. (2012). Review: Speciation by symbiosis. *Trends in Ecology & Evolution, 27*(8), 443. ScienceDirect. https://doi.org/10.1016/j.tree.2012.03.011

Cabau, C., Escudié, F., Djari, A., Guiguen, Y., Bobe, J., & Kropp, C. (2017). Compacting and correcting Trinity and Oases RNA-Seq de novo assemblies. *PeerJ*, 5, e2988.

Cahill, A. E., Meglécz, E., & Chenuil, A. (2023). Scientific history, biogeography, and biological traits predict presence of cryptic or overlooked species. *Biological Reviews, 99*, 546-561. https://doi.org/10.1111/brv.13034

Cairns, S. D. (2007). Deep-water corals: An overview with special reference to diversity and distribution of deep-water scleractinian corals. *Bulletin of Marine Science, 81*(3), 311-322.

Calderon, I., Garrabou, J., & Aurelle, D. (2006). 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, 336*, 184-197.

Campoy, A. N., Addamo, A. M., Machordom, A., Meade, A., Rivadeneira, M. M., Hernández, C. E., & Venditti, C. (2020). The Origin and Correlated Evolution of Symbiosis and Coloniality in Scleractinian Corals. *Frontiers in Marine Science, 7*.

Canessa, M., Bavestrello, G., Bo, M., Enrichetti, F., & Trainito, E. (2022). Filling a gap: A population of Eunicella verrucosa (pallas, 1766)(Anthozoa, Alcyonacea) in the Tavolara-punta Coda Cavallo marine protected area (Ne Sardinia, Italy). *Diversity, 14*(5), 405.

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

Carugati, L., Moccia, D., Bramanti, L., Cannas, R., Follesa, M. C., Salvadori, S., & Cau, A. (2022). Deep-Dwelling Populations of Mediterranean Corallium rubrum and Eunicella cavolini: Distribution, Demography, and Co-Occurrence. *Biology, 11*(2). https://doi.org/10.3390/biology11020333

Casado-Amezúa, P., Terrón-Sigler, A., Pinzón, J. H., Furla, P., Forcioli, D., Allemand, D., Ribes, M., & Coma, R. (2016). General ecological aspects of Anthozoan-Symbiodinium interactions in the Mediterranean Sea. *The Cnidaria, Past, Present and Future: The World of Medusa and Her Sisters, 375-386*.

Chan, W. Y., Peplow, L. M., Menéndez, P., Hoffmann, A. A., & Van Oppen, M. J. (2018). Interspecific hybridization may provide novel opportunities for coral reef restoration. *Frontiers in Marine Science, 5*, 160.

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

Coelho, M., Pearson, G., Boavida, J., Paulo, D., Aurelle, D., Arnaud-haond, S., Gómez-gras, D., Bensoussan, N., López-sendino, P., Cerrano, C., Kipson, S., Bakran-petricioli, T., Ferretti, E., Linares, C., Garrabou, J., Serrão, E., & Ledoux, J. (2023). Not out of the Mediterranean: Atlantic populations of the gorgonian Paramuricea clavata are a separate sister species under further
lineage diversification. *Ecology and Evolution*, 13, e9740.
https://onlinelibrary.wiley.com/doi/full/10.1002/ece3.9740#

Coma, R., Llorente-Llurba, E., Serrano, E., Gili, J.-M., & Ribes, M. (2015). Natural heterotrophic feeding by a temperate octocoral with symbiotic zooxanthellae: A contribution to understanding the mechanisms of die-off events. *Coral Reefs*, 34, 549-560.

Combosch, D. J., & Vollmer, S. V. (2015). Trans-Pacific RAD-Seq population genomics confirms introgressive hybridization in Eastern Pacific Pocillopora corals. *Molecular Phylogenetics and Evolution*, 88, 154-162.

Csilléry, K., François, O., & Blum, M. G. B. (2012). abc: An R package for approximate Bayesian computation (ABC). *Methods in Ecology and Evolution*, 3(3), 475-479. https://doi.org/10.1111/j.2041-210X.2011.00179.x

Danecek, P., Auton, A., Abecasis, G., Albers, C. A., Banks, E., DePristo, M. A., Handsaker, R. E., Lunter, G., Marth, G. T., & Sherry, S. T. (2011). The variant call format and VCFtools. *Bioinformatics*, 27(15), 2156-2158.

De Jode, A., Le Moan, A., Johannesson, K., Faria, R., Stankowski, S., Westram, A. M., Butlin, R. K., Rafajlović, M., & Fraïsse, C. (2023). Ten years of demographic modelling of divergence and speciation in the sea. *Evolutionary Applications*, 16(2), 542-559.

De Queiroz, K. (2007). Species Concepts and Species Delimitation. *Systematic Biology*, 56(6), 879-886. https://doi.org/10.1080/10635150701701083

Erickson, K. L., Pentico, A., Quattrini, A. M., & McFadden, C. S. (2021). New approaches to species delimitation and population structure of anthozoans: Two case studies of octocorals using ultraconserved elements and exons. *Molecular Ecology Resources*, 21(1), 78-92.

Estaque, T., Richaume, J., Bianchimani, O., Schull, Q., Mérigot, B., Bensoussan, N., Bonhomme, P., Vouriot, P., Sartoretto, S., & Monfort, T. (2023). Marine heatwaves on the rise: One of the strongest ever observed mass mortality event in temperate gorgonians. *Global Change Biology*, 29, 6159-6162.
https://doi.org/10.1111/gcb.16931

Ezzat, L., Merle, P.-L., Furla, P., Buttler, A., & Ferrier-Pagès, C. (2013). The Response of the Mediterranean Gorgonian Eunicella singularis to Thermal Stress Is Independent of Its Nutritional Regime. *PLoS ONE*, 8(5), e64370.
https://doi.org/10.1371/journal.pone.0064370

Faria, R., Johannesson, K., & Stankowski, S. (2021). Speciation in marine environments: Diving under the surface. *Journal of Evolutionary Biology*, 34(1), 4-15.

Forcioli, D., Merle, P. L., Caligara, C., Ciosi, M., Muti, C., Francour, P., Cerrano, C., & Allemand, D. (2011). Symbiont diversity is not involved in depth acclimation in the Mediterranean sea whip Eunicella singularis. *Marine Ecology Progress Series*, 439, 57-71.

Fourt, M., & Goujard, A. (2012). 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.

Fraïsse, C., Popovic, I., Mazoyer, C., Spataro, B., Delmotte, S., Romiguier, J., Loire, E., Simon, A., Galtier, N., & Duret, L. (2021). DILS: Demographic inferences with linked selection by using ABC. *Molecular Ecology Resources*, 21(8), 2629-2644.

François Monnet, Zoé Postel, Pascal Touzet, Christelle Fraïsse, Yves Van de Peer, Xavier Vekemans, & Camille Roux. (2023). Rapid establishment of species barriers in plants compared to animals. *bioRxiv*, 2023.10.16.562535. https://doi.org/10.1101/2023.10.16.562535
lineage diversification. *Ecology and Evolution*, 13, e9740.
https://onlinelibrary.wiley.com/doi/full/10.1002/ece3.9740#

Coma, R., Llorente-Llurba, E., Serrano, E., Gili, J.-M., & Ribes, M. (2015). Natural heterotrophic feeding by a temperate octocoral with symbiotic zooxanthellae: A contribution to understanding the mechanisms of die-off events. *Coral Reefs*, 34, 549-560.

Combosch, D. J., & Vollmer, S. V. (2015). Trans-Pacific RAD-Seq population genomics confirms introgressive hybridization in Eastern Pacific Pocillopora corals. *Molecular Phylogenetics and Evolution*, 88, 154-162.

Csilléry, K., François, O., & Blum, M. G. B. (2012). abc: An R package for approximate Bayesian computation (ABC). *Methods in Ecology and Evolution*, 3(3), 475-479. https://doi.org/10.1111/j.2041-210X.2011.00179.x

Danecek, P., Auton, A., Abecasis, G., Albers, C. A., Banks, E., DePristo, M. A., Handsaker, R. E., Lunter, G., Marth, G. T., & Sherry, S. T. (2011). The variant call format and VCFtools. *Bioinformatics*, 27(15), 2156-2158.

De Jode, A., Le Moan, A., Johannesson, K., Faria, R., Stankowski, S., Westram, A. M., Butlin, R. K., Rafajlović, M., & Fraïsse, C. (2023). Ten years of demographic modelling of divergence and speciation in the sea. *Evolutionary Applications*, 16(2), 542-559.

De Queiroz, K. (2007). Species Concepts and Species Delimitation. *Systematic Biology*, 56(6), 879-886. https://doi.org/10.1080/10635150701701083

Erickson, K. L., Pentico, A., Quattrini, A. M., & McFadden, C. S. (2021). New approaches to species delimitation and population structure of anthozoans: Two case studies of octocorals using ultraconserved elements and exons. *Molecular Ecology Resources*, 21(1), 78-92.

Estaque, T., Richaume, J., Bianchimani, O., Schull, Q., Mérigot, B., Bensoussan, N., Bonhomme, P., Vouriot, P., Sartoretto, S., & Monfort, T. (2023). Marine heatwaves on the rise: One of the strongest ever observed mass mortality event in temperate gorgonians. *Global Change Biology*, 29, 6159-6162.
https://doi.org/10.1111/gcb.16931

Ezzat, L., Merle, P.-L., Furla, P., Buttler, A., & Ferrier-Pagès, C. (2013). The Response of the Mediterranean Gorgonian Eunicella singularis to Thermal Stress Is Independent of Its Nutritional Regime. *PLoS ONE*, 8(5), e64370.
https://doi.org/10.1371/journal.pone.0064370

Faria, R., Johannesson, K., & Stankowski, S. (2021). Speciation in marine environments: Diving under the surface. *Journal of Evolutionary Biology*, 34(1), 4-15.

Forcioli, D., Merle, P. L., Caligara, C., Ciosi, M., Muti, C., Francour, P., Cerrano, C., & Allemand, D. (2011). Symbiont diversity is not involved in depth acclimation in the Mediterranean sea whip Eunicella singularis. *Marine Ecology Progress Series*, 439, 57-71.

Fourt, M., & Goujard, A. (2012). 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.

Fraïsse, C., Popovic, I., Mazoyer, C., Spataro, B., Delmotte, S., Romiguier, J., Loire, E., Simon, A., Galtier, N., & Duret, L. (2021). DILS: Demographic inferences with linked selection by using ABC. *Molecular Ecology Resources*, 21(8), 2629-2644.

François Monnet, Zoé Postel, Pascal Touzet, Christelle Fraïsse, Yves Van de Peer, Xavier Vekemans, & Camille Roux. (2023). Rapid establishment of species barriers in plants compared to animals. *bioRxiv*, 2023.10.16.562535. https://doi.org/10.1101/2023.10.16.562535
Frichot, E., & François, O. (2015). LEA: an R package for landscape and ecological association studies. *Methods in Ecology and Evolution, 6*(8), 925-929.

Frichot, E., Mathieu, F., Trouillon, T., Bouchard, G., & François, O. (2014). Fast and efficient estimation of individual ancestry coefficients. *Genetics, 196*(4), 973-983.

Furla, P., Allemand, D., Shick, J. M., Ferrier-Pagès, C., Richier, S., Plantivaux, A., Merle, P.-L., & Tambutté, S. (2005). The symbiotic anthozoan: A physiological chimera between alga and animal. *Integrative and Comparative Biology, 45*(4), 595-604.

Gagnaire, P., Broquet, T., Aurelle, D., Viard, F., Souissi, A., Bonhomme, F., Arnaud-Haond, S., & Bierne, N. (2015). Using neutral, selected, and hitchhiker loci to assess connectivity of marine populations in the genomic era. *Evolutionary Applications, 8*(8), 769-786.

Gayral, P., Melo-Ferreira, J., Glemin, S., Bierne, N., Carneiro, M., Nabholz, B., Lourenco, J. M., Alves, P. C., Ballenghien, M., & Faivre, N. (2013). Reference-free population genomics from next-generation transcriptome data and the vertebrate-invertebrate gap. *PLoS Genetics, 9*(4), e1003457.

Gori, A., Bramanti, L., López-González, P., Thoma, J. N., Gili, J.-M., Grinyó, J., Uceira, V., & Rossi, S. (2012). Characterization of the zooxanthellate and azooxanthellate morphotypes of the Mediterranean gorgonian Eunicella singularis. *Marine Biology, 159*(7), 1485-1496.

Guzman, C., Shinzato, C., Lu, T.-M., & Conaco, C. (2018). Transcriptome analysis of the reef-building octocoral, Heliopora coerulea. *Scientific Reports, 8*(1), 8397. https://doi.org/10.1038/s41598-018-26718-5

Haguenauer, A., Zuberer, F., Ledoux, J.-B., & Aurelle, D. (2013). 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, 449*(0), 349-357.
https://doi.org/10.1016/j.jembe.2013.10.010

Johnston, E. C., Forsman, Z. H., Flot, J.-F., Schmidt-Roach, S., Pinzón, J. H., Knapp, I. S., & Toonen, R. J. (2017). A genomic glance through the fog of plasticity and diversification in Pocillopora. *Scientific Reports, 7*(1), 5991.

Kent, W. J. (2002). BLAT—the BLAST-like alignment tool. *Genome Research, 12*(4), 656-664.

Krueger-Hadfield, S. (2015). Marmap. https://www.molecularecologist.com/2015/07/03/marmap/

LaJeunesse, T. C., Parkinson, J. E., Gabrielson, P. W., Jeong, H. J., Reimer, J. D., Voolstra, C. R., & Santos, S. R. (2018). Systematic revision of Symbiodiniaceae highlights the antiquity and diversity of coral endosymbionts. *Current Biology, 28*(16), 2570-2580.

LaJeunesse, T. C., Wiedenmann, J., Casado-Amezúa, P., D'ambra, I., Turnham, K. E., Nitschke, M. R., Oakley, C. A., Goffredo, S., Spano, C. A., & Cubillos, V. M. (2022). Revival of Philozoon Geddes for host-specialized dinoflagellates, 'zooxanthellae', in animals from coastal temperate zones of northern and southern hemispheres. *European Journal of Phycology, 57*(2), 166-180.

Lee, M. J., Jeong, H. J., Jang, S. H., Lee, S. Y., Kang, N. S., Lee, K. H., Kim, H. S., Wham, D. C., & LaJeunesse, T. C. (2016). Most low-abundance "background" Symbiodinium spp. are transitory and have minimal functional significance for symbiotic corals. *Microbial Ecology, 71*, 771-783.

Leroy, T., Louvet, J.-M., Lalanne, C., Le Provost, G., Labadie, K., Aury, J.-M., Delzon, S., Plomion, C., & Kremer, A. (2020). Adaptive introgression as a
Frichot, E., & François, O. (2015). LEA: an R package for landscape and ecological association studies. *Methods in Ecology and Evolution*, **6**(8), 925-929.

Frichot, E., Mathieu, F., Trouillon, T., Bouchard, G., & François, O. (2014). Fast and efficient estimation of individual ancestry coefficients. *Genetics*, **196**(4), 973-983.

Furla, P., Allemand, D., Shick, J. M., Ferrier-Pagès, C., Richier, S., Plantivaux, A., Merle, P.-L., & Tambutté, S. (2005). The symbiotic anthozoan: A physiological chimera between alga and animal. *Integrative and Comparative Biology*, **45**(4), 595-604.

Gagnaire, P., Broquet, T., Aurelle, D., Viard, F., Souissi, A., Bonhomme, F., Arnaud-Haond, S., & Bierne, N. (2015). Using neutral, selected, and hitchhiker loci to assess connectivity of marine populations in the genomic era. *Evolutionary Applications*, **8**(8), 769-786.

Gayral, P., Melo-Ferreira, J., Glemin, S., Bierne, N., Carneiro, M., Nabholz, B., Lourenco, J. M., Alves, P. C., Ballenghien, M., & Faivre, N. (2013). Reference-free population genomics from next-generation transcriptome data and the vertebrate-invertebrate gap. *PLoS Genetics*, **9**(4), e1003457.

Gori, A., Bramanti, L., López-González, P., Thoma, J. N., Gili, J.-M., Grinyó, J., Uceira, V., & Rossi, S. (2012). Characterization of the zooxanthellate and azooxanthellate morphotypes of the Mediterranean gorgonian Eunicella singularis. *Marine Biology*, **159**(7), 1485-1496.

Guzman, C., Shinzato, C., Lu, T.-M., & Conaco, C. (2018). Transcriptome analysis of the reef-building octocoral, Heliopora coerulea. *Scientific Reports*, **8**(1), 8397. https://doi.org/10.1038/s41598-018-26718-5

Haguenauer, A., Zuberer, F., Ledoux, J.-B., & Aurelle, D. (2013). 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*, **449**(0), 349-357.
https://doi.org/10.1016/j.jembe.2013.10.010

Johnston, E. C., Forsman, Z. H., Flot, J.-F., Schmidt-Roach, S., Pinzón, J. H., Knapp, I. S., & Toonen, R. J. (2017). A genomic glance through the fog of plasticity and diversification in Pocillopora. *Scientific Reports*, **7**(1), 5991.

Kent, W. J. (2002). BLAT—the BLAST-like alignment tool. *Genome Research*, **12**(4), 656-664.

Krueger-Hadfield, S. (2015). Marmap. https://www.molecularecologist.com/2015/07/03/marmap/

LaJeunesse, T. C., Parkinson, J. E., Gabrielson, P. W., Jeong, H. J., Reimer, J. D., Voolstra, C. R., & Santos, S. R. (2018). Systematic revision of Symbiodiniaceae highlights the antiquity and diversity of coral endosymbionts. *Current Biology*, **28**(16), 2570-2580.

LaJeunesse, T. C., Wiedenmann, J., Casado-Amezúa, P., D'ambra, I., Turnham, K. E., Nitschke, M. R., Oakley, C. A., Goffredo, S., Spano, C. A., & Cubillos, V. M. (2022). Revival of Philozoon Geddes for host-specialized dinoflagellates, 'zooxanthellae', in animals from coastal temperate zones of northern and southern hemispheres. *European Journal of Phycology*, **57**(2), 166-180.

Lee, M. J., Jeong, H. J., Jang, S. H., Lee, S. Y., Kang, N. S., Lee, K. H., Kim, H. S., Wham, D. C., & LaJeunesse, T. C. (2016). Most low-abundance "background" Symbiodinium spp. are transitory and have minimal functional significance for symbiotic corals. *Microbial Ecology*, **71**, 771-783.

Leroy, T., Louvet, J.-M., Lalanne, C., Le Provost, G., Labadie, K., Aury, J.-M., Delzon, S., Plomion, C., & Kremer, A. (2020). Adaptive introgression as a
driver of local adaptation to climate in European white oaks. *New Phytologist*, 226(4), 1171-1182. https://doi.org/10.1111/nph.16095

Lesser, M. P., Stat, M., & Gates, R. D. (2013). The endosymbiotic dinoflagellates (Symbiodinium sp.) of corals are parasites and mutualists. *Coral Reefs, 32*(3), 603-611. https://doi.org/10.1007/s00338-013-1051-z

Li, H., & Durbin, R. (2009). Fast and accurate short read alignment with Burrows-Wheeler transform. *Bioinformatics*, 25(14), 1754-1760.

Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., Marth, G., Abecasis, G., & Durbin, R. (2009). The sequence alignment/map format and SAMtools. *Bioinformatics*, 25(16), 2078-2079.

Lischer, H. E., & Excoffier, L. (2012). PGDSpider: An automated data conversion tool for connecting population genetics and genomics programs. *Bioinformatics*, 28(2), 298-299.

Macleod, K. L., Jenkins, T. L., Witt, M. J., & Stevens, J. R. (2024). Rare, long-distance dispersal underpins genetic connectivity in the pink sea fan, *Eunicella verrucosa*. *Evolutionary Applications*, 17(3), e13649.
https://doi.org/10.1111/eva.13649

Mayr, E. (2001). Wu's genetic view of speciation. *Journal of Evolutionary Biology*, 14(6), 866-867. https://doi.org/10.1046/j.1420-9101.2001.00336.x

McFadden, C., Benayahu, Y., Pante, E., Thoma, J. N., Nevarez, P. A., & France, S. C. (2011). Limitations of mitochondrial gene barcoding in Octocorallia. *Molecular Ecology Resources*, 11(1), 19-31. https://doi.org/10.1111/j.1755-0998.2010.02875.x

Muir, P. R., Obura, D. O., Hoeksema, B. W., Sheppard, C., Pichon, M., & Richards, Z. T. (2022). Conclusions of low extinction risk for most species of reef-building corals are premature. *Nature Ecology & Evolution*, 6(4), 357-358.

Munro, L. (2004). Determining the reproductive cycle of Eunicella verrucosa. *Reef Research: ETR*, 11.
https://www.marine-bio-images.com/RR_Eunicella_PDFS/Report_RR12Jul2004reproductive%20cycle%20pdf.pdf

Muthye, V., Mackereth, C. D., Stewart, J. B., & Lavrov, D. V. (2022). Large dataset of octocoral mitochondrial genomes provides new insights into mt-mutS evolution and function. *DNA Repair*, 110, 103273.

Palumbi, S. R. (1992). Marine speciation on a small planet. *Trends in Ecology & Evolution*, 7(4), 114-118.

Palumbi, S. R. (1999). All males are not created equal: Fertility differences depend on gamete recognition polymorphisms in sea urchins. *Proceedings of the National Academy of Sciences*, 96(22), 12632-12637.

Pante, E., Abdelkrim, J., Viricel, A., Gey, D., France, S., Boisselier, M.-C., & Samadi, S. (2015a). Use of RAD sequencing for delimiting species. *Heredity*, 114(5), 450-459.

Pante, E., Puillandre, N., Viricel, A., Arnaud-Haond, S., Aurelle, D., Castelin, M., Chenuil, A., Destombe, C., Forcioli, D., & Valero, M. (2015b). Species are hypotheses: Avoid connectivity assessments based on pillars of sand. *Molecular Ecology*, 24(3), 525-544.

Pante, E., & Simon-Bouhet, B. (2013). marmap: A package for importing, plotting and analyzing bathymetric and topographic data in R. *PLoS One*, 8(9), e73051. https://doi.org/10.1371/journal.pone.0073051

Patro, R., Duggal, G., Love, M. I., Irizarry, R. A., & Kingsford, C. (2017). Salmon provides fast and bias-aware quantification of transcript expression. *Nature Methods*, 14(4), 417-419.

Peñalba, J. V., Runemark, A., Meier, J. I., Singh, P., Wogan, G. O., Sánchez-Guillén, R., Mallet, J., Rometsch, S. J., Menon, M., & Seehausen, O. (2024). The Role
driver of local adaptation to climate in European white oaks. *New Phytologist*, 226(4), 1171-1182. https://doi.org/10.1111/nph.16095

Lesser, M. P., Stat, M., & Gates, R. D. (2013). The endosymbiotic dinoflagellates (Symbiodinium sp.) of corals are parasites and mutualists. *Coral Reefs, 32*(3), 603-611. https://doi.org/10.1007/s00338-013-1051-z

Li, H., & Durbin, R. (2009). Fast and accurate short read alignment with Burrows-Wheeler transform. *Bioinformatics*, 25(14), 1754-1760.

Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., Marth, G., Abecasis, G., & Durbin, R. (2009). The sequence alignment/map format and SAMtools. *Bioinformatics*, 25(16), 2078-2079.

Lischer, H. E., & Excoffier, L. (2012). PGDSpider: An automated data conversion tool for connecting population genetics and genomics programs. *Bioinformatics*, 28(2), 298-299.

Macleod, K. L., Jenkins, T. L., Witt, M. J., & Stevens, J. R. (2024). Rare, long-distance dispersal underpins genetic connectivity in the pink sea fan, *Eunicella verrucosa*. *Evolutionary Applications*, 17(3), e13649.
https://doi.org/10.1111/eva.13649

Mayr, E. (2001). Wu's genetic view of speciation. *Journal of Evolutionary Biology*, 14(6), 866-867. https://doi.org/10.1046/j.1420-9101.2001.00336.x

McFadden, C., Benayahu, Y., Pante, E., Thoma, J. N., Nevarez, P. A., & France, S. C. (2011). Limitations of mitochondrial gene barcoding in Octocorallia. *Molecular Ecology Resources*, 11(1), 19-31. https://doi.org/10.1111/j.1755-0998.2010.02875.x

Muir, P. R., Obura, D. O., Hoeksema, B. W., Sheppard, C., Pichon, M., & Richards, Z. T. (2022). Conclusions of low extinction risk for most species of reef-building corals are premature. *Nature Ecology & Evolution*, 6(4), 357-358.

Munro, L. (2004). Determining the reproductive cycle of Eunicella verrucosa. *Reef Research: ETR*, 11.
https://www.marine-bio-images.com/RR_Eunicella_PDFS/Report_RR12Jul2004reproductive%20cycle%20pdf.pdf

Muthye, V., Mackereth, C. D., Stewart, J. B., & Lavrov, D. V. (2022). Large dataset of octocoral mitochondrial genomes provides new insights into mt-mutS evolution and function. *DNA Repair*, 110, 103273.

Palumbi, S. R. (1992). Marine speciation on a small planet. *Trends in Ecology & Evolution*, 7(4), 114-118.

Palumbi, S. R. (1999). All males are not created equal: Fertility differences depend on gamete recognition polymorphisms in sea urchins. *Proceedings of the National Academy of Sciences*, 96(22), 12632-12637.

Pante, E., Abdelkrim, J., Viricel, A., Gey, D., France, S., Boisselier, M.-C., & Samadi, S. (2015a). Use of RAD sequencing for delimiting species. *Heredity*, 114(5), 450-459.

Pante, E., Puillandre, N., Viricel, A., Arnaud-Haond, S., Aurelle, D., Castelin, M., Chenuil, A., Destombe, C., Forcioli, D., & Valero, M. (2015b). Species are hypotheses: Avoid connectivity assessments based on pillars of sand. *Molecular Ecology*, 24(3), 525-544.

Pante, E., & Simon-Bouhet, B. (2013). marmap: A package for importing, plotting and analyzing bathymetric and topographic data in R. *PLoS One*, 8(9), e73051. https://doi.org/10.1371/journal.pone.0073051

Patro, R., Duggal, G., Love, M. I., Irizarry, R. A., & Kingsford, C. (2017). Salmon provides fast and bias-aware quantification of transcript expression. *Nature Methods*, 14(4), 417-419.

Peñalba, J. V., Runemark, A., Meier, J. I., Singh, P., Wogan, G. O., Sánchez-Guillén, R., Mallet, J., Rometsch, S. J., Menon, M., & Seehausen, O. (2024). The Role
of Hybridization in Species Formation and Persistence. *Cold Spring Harbor Perspectives in Biology*, a041445.

*Picard Toolkit*. (2019). Broad Institute. GitHub Repository. https://broadinstitute.github.io/picard/

Porro, B. (2019). *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)).

Pudlo, P., Marin, J.-M., Estoup, A., Cornuet, J.-M., Gautier, M., & Robert, C. P. (2016). Reliable ABC model choice via random forests. *Bioinformatics*, 32(6), 859-866.

Quilodrán, C. S., Ruegg, K., Sendell-Price, A. T., Anderson, E. C., Coulson, T., & Clegg, S. M. (2020). The multiple population genetic and demographic routes to islands of genomic divergence. *Methods in Ecology and Evolution, 11*(1), 6-21. https://doi.org/10.1111/2041-210X.13324

Readman, J., & Hiscock, K. (2017). *Eunicella verrucosa*. Pink sea fan. https://www.marlin.ac.uk/species/detail/1121

Reynes, L., Aurelle, D., Chevalier, C., Pinazo, C., Valero, M., Mauger, S., Sartoretto, S., Blafuné, A., Ruitton, S., Boudouresque, C.-F., Verlaque, M., & Thibaut, T. (2021). Population genomics and Lagrangian modeling shed light on dispersal events in the Mediterranean endemic Ericaria zosteroides (= Cystoseira zosteroides) (Fucales). *Frontiers in Marine Science*, 8, 683528.

Rivera-García, L., Rivera-Vicéns, R. E., Veglia, A. J., & Schizas, N. V. (2019). De novo transcriptome assembly of the digitate morphotype of Briareum asbestinum (Octocorallia: Alcyonacea) from the southwest shelf of Puerto Rico. *Marine Genomics*, 47, 100676. https://doi.org/10.1016/j.margen.2019.04.001

Rosenberg, E., & Zilber-Rosenberg, I. (2018). The hologenome concept of evolution after 10 years. *Microbiome*, 6(1), 78.

Rouchy, J. M., & Caruso, A. (2006). The Messinian salinity crisis in the Mediterranean basin: A reassessment of the data and an integrated scenario. *Sedimentary Geology*, 188, 35-67.

Rousset, F. (2008). *genepop'007: A complete re-implementation of the genepop software for Windows and Linux*. *Molecular Ecology Resources*, 8(1), 103-106. https://doi.org/10.1111/j.1471-8286.2007.01931.x

Rousset, F., Lopez, J., & Belkhir, K. (2020). Package 'genepop.' *R Package Version*, 1(7).

Roux, C., Fraïsse, C., Romiguier, J., Anciaux, Y., Galtier, N., & Bierne, N. (2016). Shedding Light on the Grey Zone of Speciation along a Continuum of Genomic Divergence. *PLOS Biology*, 14(12), e2000234. https://doi.org/10.1371/journal.pbio.2000234

Sachs, J. L., & Wilcox, T. P. (2006). A shift to parasitism in the jellyfish symbiont *Symbiodinium microadriaticum*. *Proceedings of the Royal Society B: Biological Sciences*, 273(1585), 425-429.

Samadi, S., & Barberousse, A. (2006). The tree, the network, and the species. *Biological Journal of the Linnean Society*, 89(3), 509-521.

Sartoretto, S., & Francour, P. (2011). Bathymetric distribution and growth rates of Eunicella verrucosa (Cnidaria: Gorgoniidae) populations along the Marseilles coast (France). *Scientia Marina*, 76(2), 349-355.

Schulz, M. H., Zerbino, D. R., Vingron, M., & Birney, E. (2012). Oases: Robust de novo RNA-seq assembly across the dynamic range of expression levels. *Bioinformatics*, 28(8), 1086-1092.

Sini, M., Kipson, S., Linares, C., Koutsoubas, D., & Garrabou, J. (2015). The Yellow Gorgonian Eunicella cavolini: Demography and Disturbance Levels across
of Hybridization in Species Formation and Persistence. *Cold Spring Harbor Perspectives in Biology*, a041445.

*Picard Toolkit*. (2019). Broad Institute. GitHub Repository. https://broadinstitute.github.io/picard/

Porro, B. (2019). *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)).

Pudlo, P., Marin, J.-M., Estoup, A., Cornuet, J.-M., Gautier, M., & Robert, C. P. (2016). Reliable ABC model choice via random forests. *Bioinformatics*, 32(6), 859-866.

Quilodrán, C. S., Ruegg, K., Sendell-Price, A. T., Anderson, E. C., Coulson, T., & Clegg, S. M. (2020). The multiple population genetic and demographic routes to islands of genomic divergence. *Methods in Ecology and Evolution, 11*(1), 6-21. https://doi.org/10.1111/2041-210X.13324

Readman, J., & Hiscock, K. (2017). *Eunicella verrucosa*. Pink sea fan. https://www.marlin.ac.uk/species/detail/1121

Reynes, L., Aurelle, D., Chevalier, C., Pinazo, C., Valero, M., Mauger, S., Sartoretto, S., Blafuné, A., Ruitton, S., Boudouresque, C.-F., Verlaque, M., & Thibaut, T. (2021). Population genomics and Lagrangian modeling shed light on dispersal events in the Mediterranean endemic Ericaria zosteroides (= Cystoseira zosteroides) (Fucales). *Frontiers in Marine Science*, 8, 683528.

Rivera-García, L., Rivera-Vicéns, R. E., Veglia, A. J., & Schizas, N. V. (2019). De novo transcriptome assembly of the digitate morphotype of Briareum asbestinum (Octocorallia: Alcyonacea) from the southwest shelf of Puerto Rico. *Marine Genomics*, 47, 100676. https://doi.org/10.1016/j.margen.2019.04.001

Rosenberg, E., & Zilber-Rosenberg, I. (2018). The hologenome concept of evolution after 10 years. *Microbiome*, 6(1), 78.

Rouchy, J. M., & Caruso, A. (2006). The Messinian salinity crisis in the Mediterranean basin: A reassessment of the data and an integrated scenario. *Sedimentary Geology*, 188, 35-67.

Rousset, F. (2008). *genepop'007: A complete re-implementation of the genepop software for Windows and Linux*. *Molecular Ecology Resources*, 8(1), 103-106. https://doi.org/10.1111/j.1471-8286.2007.01931.x

Rousset, F., Lopez, J., & Belkhir, K. (2020). Package 'genepop.' *R Package Version*, 1(7).

Roux, C., Fraïsse, C., Romiguier, J., Anciaux, Y., Galtier, N., & Bierne, N. (2016). Shedding Light on the Grey Zone of Speciation along a Continuum of Genomic Divergence. *PLOS Biology*, 14(12), e2000234. https://doi.org/10.1371/journal.pbio.2000234

Sachs, J. L., & Wilcox, T. P. (2006). A shift to parasitism in the jellyfish symbiont *Symbiodinium microadriaticum*. *Proceedings of the Royal Society B: Biological Sciences*, 273(1585), 425-429.

Samadi, S., & Barberousse, A. (2006). The tree, the network, and the species. *Biological Journal of the Linnean Society*, 89(3), 509-521.

Sartoretto, S., & Francour, P. (2011). Bathymetric distribution and growth rates of Eunicella verrucosa (Cnidaria: Gorgoniidae) populations along the Marseilles coast (France). *Scientia Marina*, 76(2), 349-355.

Schulz, M. H., Zerbino, D. R., Vingron, M., & Birney, E. (2012). Oases: Robust de novo RNA-seq assembly across the dynamic range of expression levels. *Bioinformatics*, 28(8), 1086-1092.

Sini, M., Kipson, S., Linares, C., Koutsoubas, D., & Garrabou, J. (2015). The Yellow Gorgonian Eunicella cavolini: Demography and Disturbance Levels across
the Mediterranean Sea. *PLoS ONE*, 10(5), e0126253.
https://doi.org/10.1371/journal.pone.0126253

Stanley, R. R. E., Jeffery, N. W., Wringe, B. F., DiBacco, C., & Bradbury, I. R. (2017). *genepopedit: A simple and flexible tool for manipulating multilocus molecular data in R*. *Molecular Ecology Resources*, *17*(1), 12–18.
https://doi.org/10.1111/1755-0998.12569

Tricou, T., Tannier, E., & de Vienne, D. M. (2022). Ghost lineages can invalidate or even reverse findings regarding gene flow. *PLoS Biology*, *20*(9), e3001776.

Tsagkogeorga, G., Cahais, V., & Galtier, N. (2012). 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*, *4*(8), 852–861.

van Oppen, M. J. H., & Medina, M. (2020). Coral evolutionary responses to microbial symbioses. *Philosophical Transactions of the Royal Society B: Biological Sciences*, 375(1808), 20190591.
https://doi.org/10.1098/rstb.2019.0591

Van Oppen, M., Mieog, J. C., Sanchez, C., & Fabricius, K. (2005). Diversity of algal endosymbionts (zooxanthellae) in octocorals: The roles of geography and host relationships. *Molecular Ecology*, *14*(8), 2403–2417.

Wagner, D., Pochon, X., Irwin, L., Toonen, R. J., & Gates, R. D. (2011). Azooxanthellate? Most Hawaiian black corals contain Symbiodinium. *Proceedings of the Royal Society B: Biological Sciences*, 278(1710), 1323–1328.

Weir, B. S., & Cockerham, C. C. (1984). Estimating F-statistics for the analysis of population structure. *Evolution*, *38*, 1358–1370.
the Mediterranean Sea. *PLoS ONE*, 10(5), e0126253.
https://doi.org/10.1371/journal.pone.0126253

Stanley, R. R. E., Jeffery, N. W., Wringe, B. F., DiBacco, C., & Bradbury, I. R. (2017). *genepopedit: A simple and flexible tool for manipulating multilocus molecular data in R*. *Molecular Ecology Resources*, *17*(1), 12–18.
https://doi.org/10.1111/1755-0998.12569

Tricou, T., Tannier, E., & de Vienne, D. M. (2022). Ghost lineages can invalidate or even reverse findings regarding gene flow. *PLoS Biology*, *20*(9), e3001776.

Tsagkogeorga, G., Cahais, V., & Galtier, N. (2012). 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*, *4*(8), 852–861.

van Oppen, M. J. H., & Medina, M. (2020). Coral evolutionary responses to microbial symbioses. *Philosophical Transactions of the Royal Society B: Biological Sciences*, 375(1808), 20190591.
https://doi.org/10.1098/rstb.2019.0591

Van Oppen, M., Mieog, J. C., Sanchez, C., & Fabricius, K. (2005). Diversity of algal endosymbionts (zooxanthellae) in octocorals: The roles of geography and host relationships. *Molecular Ecology*, *14*(8), 2403–2417.

Wagner, D., Pochon, X., Irwin, L., Toonen, R. J., & Gates, R. D. (2011). Azooxanthellate? Most Hawaiian black corals contain Symbiodinium. *Proceedings of the Royal Society B: Biological Sciences*, 278(1710), 1323–1328.

Weir, B. S., & Cockerham, C. C. (1984). Estimating F-statistics for the analysis of population structure. *Evolution*, *38*, 1358–1370.
684

**Figure 1:** map of sampling sites; A) general view, B) zoom on the area of

685

Marseille. The symbols present the different samples: EC *E. cavolini*, ES

686

*E. singularis*, EV *E. verrucosa*, HY hybrids, EC-EV *E. cavolini* and *E. verrucosa* in

687

sympatry, ES-HY *E. singularis* and hybrids in sympatry. The three letters

688

correspond to the codes of the sampling. The maps have been produced with the

689

mmapart R package (Pante & Simon-Bouhet, 2013) and following the tutorial of

690

Krueger-Hadfield (2015).

691

A)

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

**Figure 1:** map of sampling sites; A) general view, B) zoom on the area of

685

Marseille. The symbols present the different samples: EC *E. cavolini*, ES

686

*E. singularis*, EV *E. verrucosa*, HY hybrids, EC-EV *E. cavolini* and *E. verrucosa* in

687

sympatry, ES-HY *E. singularis* and hybrids in sympatry. The three letters

688

correspond to the codes of the sampling. The maps have been produced with the

689

mmapart R package (Pante & Simon-Bouhet, 2013) and following the tutorial of

690

Krueger-Hadfield (2015).

691

A)

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

694

B)

695

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

696
693

694

B)

695

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

696
697 **Figure 2:** barplots of coancestry coefficients inferred with the LEA R package.

698 The analysis is based on the "polymorphic sites" dataset, with 31369 SNPs. The

699 red asterisks indicate the individuals used as prior for parental status in the

700 newhybrids analysis. The results of the newhybrids analysis are indicated below

701 the hybrid individuals: F1, 1st generation; F2, 2nd generation; Sbx, backcross with

702 *E. singularis*. Note that the coancestry analysis is based on 31 369 SNPs,

703 whereas the newhybrids analysis is based on 326 SNPs showing high

704 differentiation between *E. cavolini* and *E. singularis*.
697 **Figure 2:** barplots of coancestry coefficients inferred with the LEA R package.

698 The analysis is based on the "polymorphic sites" dataset, with 31369 SNPs. The

699 red asterisks indicate the individuals used as prior for parental status in the

700 newhybrids analysis. The results of the newhybrids analysis are indicated below

701 the hybrid individuals: F1, 1st generation; F2, 2nd generation; Sbx, backcross with

702 *E. singularis*. Note that the coancestry analysis is based on 31 369 SNPs,

703 whereas the newhybrids analysis is based on 326 SNPs showing high

704 differentiation between *E. cavolini* and *E. singularis*.
<figure><img src="image_4.png" /></figure>

E. cavolini

hybrids

E. singularis

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

E. cavolini

hybrids

E. singularis

E. verrucosa
706

707

708

709

710

711

<table><caption>**Table 1:** results of demographic inferences with DILS. The columns indicate the species comparison, the model choice for population size (constant vs. variable), and the results of inferences: current gene flow (migration vs isolation); if current migration, isolation / migration (IM) vs ancestral migration (AM); if no current migration, strict isolation (SI) vs ancestral migration (AM); homogeneity (N-homo) vs heterogeneity in effective size (N-hetero) among loci; homogeneity (M-homo) vs heterogeneity (M-hetero) in gene flow among loci. The probability of each scenario is given in the same case.</caption><thead><tr><th>Comparison</th><th>Population size</th><th>Current gene flow</th><th>IM / SC</th><th>SI / AM</th><th>Heterogeneity effective size</th><th>Heterogeneity gene flow</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>

712

713

714
706

707

708

709

710

711

<table><caption>**Table 1:** results of demographic inferences with DILS. The columns indicate the species comparison, the model choice for population size (constant vs. variable), and the results of inferences: current gene flow (migration vs isolation); if current migration, isolation / migration (IM) vs ancestral migration (AM); if no current migration, strict isolation (SI) vs ancestral migration (AM); homogeneity (N-homo) vs heterogeneity in effective size (N-hetero) among loci; homogeneity (M-homo) vs heterogeneity (M-hetero) in gene flow among loci. The probability of each scenario is given in the same case.</caption><thead><tr><th>Comparison</th><th>Population size</th><th>Current gene flow</th><th>IM / SC</th><th>SI / AM</th><th>Heterogeneity effective size</th><th>Heterogeneity gene flow</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>

712

713

714
**Symbiotic status does not preclude hybridisation in Mediterranean octocorals**

**Supplementary Material S1**

Didier Aurelle<fnref n="1" />,<fnref n="2" />*, Anne Haguenauer<fnref n="3" />, Marc Bally<fnref n="1" />,
Frédéric Zuberer<fnref n="4" />, Jean-Baptiste Ledoux<fnref n="5" />, Stéphane
Sartoretto<fnref n="6" />, Lamya Chaoui<fnref n="7" />, Hichem Kara<fnref n="7" />, Sarah
Samadi<fnref n="2" />, Pierre Pontarotti<fnref n="8" />,<fnref n="9" />,<fnref n="10" />

<figure><img src="image_5.png" /><figcaption>CC-BY 4.0 https://creativecommons.org/licenses/by/4.0/</figcaption></figure>

<footnote n="1">¹ Aix Marseille Univ, Université de Toulon, CNRS, IRD, MIO, Marseille, France</footnote>

<footnote n="2">² Institut Systématique Evolution Biodiversité (ISYEB), Muséum national d'Histoire naturelle, CNRS, Sorbonne Université, EPHE, Université des Antilles, CP 26, 75005 Paris, France.</footnote>

<footnote n="3">³ CNRS - Délégation Provence et corse, Marseille, France</footnote>

<footnote n="4">⁴ Aix Marseille Univ, CNRS, IRD, INRAE, OSU Inst. PYTHEAS, Marseille, France</footnote>

<footnote n="5">⁵ CIIMAR/CIMAR, Centro Interdisciplinar de Investigação Marinha e Ambiental, Universidade do Porto, Porto, Portugal.</footnote>

<footnote n="6">⁶ Ifremer, LITTORAL, 83500 La Seyne-sur-Mer, France</footnote>

<footnote n="7">⁷ Laboratoire Bioressources marines. Université d'Annaba Badji Mokhtar, Annaba - Algérie.</footnote>

<footnote n="8">⁸ Aix Marseille Univ, MEPHI, Marseille, France.</footnote>

<footnote n="9">⁹ IHU Méditerranée Infection, Marseille, France.</footnote>

<footnote n="10">¹⁰ CNRS SNC5039</footnote>

<footnote>*Corresponding author</footnote>

<footnote>Correspondence: didier.aurelle@univ-amu.fr</footnote>
**Symbiotic status does not preclude hybridisation in Mediterranean octocorals**

**Supplementary Material S1**

Didier Aurelle<fnref n="1" />,<fnref n="2" />*, Anne Haguenauer<fnref n="3" />, Marc Bally<fnref n="1" />,
Frédéric Zuberer<fnref n="4" />, Jean-Baptiste Ledoux<fnref n="5" />, Stéphane
Sartoretto<fnref n="6" />, Lamya Chaoui<fnref n="7" />, Hichem Kara<fnref n="7" />, Sarah
Samadi<fnref n="2" />, Pierre Pontarotti<fnref n="8" />,<fnref n="9" />,<fnref n="10" />

<figure><img src="image_6.png" /><figcaption>CC-BY 4.0 https://creativecommons.org/licenses/by/4.0/</figcaption></figure>

<footnote n="1">¹ Aix Marseille Univ, Université de Toulon, CNRS, IRD, MIO, Marseille, France</footnote>

<footnote n="2">² Institut Systématique Evolution Biodiversité (ISYEB), Muséum national d'Histoire naturelle, CNRS, Sorbonne Université, EPHE, Université des Antilles, CP 26, 75005 Paris, France.</footnote>

<footnote n="3">³ CNRS - Délégation Provence et corse, Marseille, France</footnote>

<footnote n="4">⁴ Aix Marseille Univ, CNRS, IRD, INRAE, OSU Inst. PYTHEAS, Marseille, France</footnote>

<footnote n="5">⁵ CIIMAR/CIMAR, Centro Interdisciplinar de Investigação Marinha e Ambiental, Universidade do Porto, Porto, Portugal.</footnote>

<footnote n="6">⁶ Ifremer, LITTORAL, 83500 La Seyne-sur-Mer, France</footnote>

<footnote n="7">⁷ Laboratoire Bioressources marines. Université d'Annaba Badji Mokhtar, Annaba - Algérie.</footnote>

<footnote n="8">⁸ Aix Marseille Univ, MEPHI, Marseille, France.</footnote>

<footnote n="9">⁹ IHU Méditerranée Infection, Marseille, France.</footnote>

<footnote n="10">¹⁰ CNRS SNC5039</footnote>

<footnote>*Corresponding author</footnote>

<footnote>Correspondence: didier.aurelle@univ-amu.fr</footnote>
<table><caption>**Table S1:** list of mitochondrial MutS sequences used for the phylogenetic reconstruction with the corresponding Genbank accession numbers.</caption><thead><tr><td>Accession number</td><td>Genus</td><td>species</td></tr></thead><tbody><tr><td>KP036906</td><td>Complexum</td><td>monodi</td></tr><tr><td>NC_035666</td><td>Eunicella</td><td>albicans</td></tr><tr><td>JQ397290</td><td>Eunicella</td><td>cavolini</td></tr><tr><td>JQ397291</td><td>Eunicella</td><td>cavolini</td></tr><tr><td>JQ397292</td><td>Eunicella</td><td>cavolini</td></tr><tr><td>NC_035667</td><td>Eunicella</td><td>cavolinii</td></tr><tr><td>KX051577</td><td>Eunicella</td><td>racemosa</td></tr><tr><td>JQ397293</td><td>Eunicella</td><td>singularis</td></tr><tr><td>JQ397294</td><td>Eunicella</td><td>singularis</td></tr><tr><td>KX051571</td><td>Eunicella</td><td>singularis</td></tr><tr><td>KX051572</td><td>Eunicella</td><td>singularis</td></tr><tr><td>JQ397307</td><td>Eunicella</td><td>sp.</td></tr><tr><td>JQ397308</td><td>Eunicella</td><td>sp.</td></tr><tr><td>JQ397311</td><td>Eunicella</td><td>sp.</td></tr><tr><td>JX203795</td><td>Eunicella</td><td>tricoronata</td></tr><tr><td>NC_062012</td><td>Eunicella</td><td>tricoronata</td></tr><tr><td>JQ397300</td><td>Eunicella</td><td>verrucosa</td></tr><tr><td>JQ397302</td><td>Eunicella</td><td>verrucosa</td></tr><tr><td>JQ397305</td><td>Eunicella</td><td>verrucosa</td></tr><tr><td>JQ397306</td><td>Eunicella</td><td>verrucosa</td></tr><tr><td>NC_073494</td><td>Eunicella</td><td>verrucosa</td></tr><tr><td>KX904973</td><td>Swiftia</td><td>pacifica</td></tr><tr><td>KX905018</td><td>Swiftia</td><td>simplex</td></tr></tbody></table>
<table><caption>**Table S1:** list of mitochondrial MutS sequences used for the phylogenetic reconstruction with the corresponding Genbank accession numbers.</caption><thead><tr><td>Accession number</td><td>Genus</td><td>species</td></tr></thead><tbody><tr><td>KP036906</td><td>Complexum</td><td>monodi</td></tr><tr><td>NC_035666</td><td>Eunicella</td><td>albicans</td></tr><tr><td>JQ397290</td><td>Eunicella</td><td>cavolini</td></tr><tr><td>JQ397291</td><td>Eunicella</td><td>cavolini</td></tr><tr><td>JQ397292</td><td>Eunicella</td><td>cavolini</td></tr><tr><td>NC_035667</td><td>Eunicella</td><td>cavolinii</td></tr><tr><td>KX051577</td><td>Eunicella</td><td>racemosa</td></tr><tr><td>JQ397293</td><td>Eunicella</td><td>singularis</td></tr><tr><td>JQ397294</td><td>Eunicella</td><td>singularis</td></tr><tr><td>KX051571</td><td>Eunicella</td><td>singularis</td></tr><tr><td>KX051572</td><td>Eunicella</td><td>singularis</td></tr><tr><td>JQ397307</td><td>Eunicella</td><td>sp.</td></tr><tr><td>JQ397308</td><td>Eunicella</td><td>sp.</td></tr><tr><td>JQ397311</td><td>Eunicella</td><td>sp.</td></tr><tr><td>JX203795</td><td>Eunicella</td><td>tricoronata</td></tr><tr><td>NC_062012</td><td>Eunicella</td><td>tricoronata</td></tr><tr><td>JQ397300</td><td>Eunicella</td><td>verrucosa</td></tr><tr><td>JQ397302</td><td>Eunicella</td><td>verrucosa</td></tr><tr><td>JQ397305</td><td>Eunicella</td><td>verrucosa</td></tr><tr><td>JQ397306</td><td>Eunicella</td><td>verrucosa</td></tr><tr><td>NC_073494</td><td>Eunicella</td><td>verrucosa</td></tr><tr><td>KX904973</td><td>Swiftia</td><td>pacifica</td></tr><tr><td>KX905018</td><td>Swiftia</td><td>simplex</td></tr></tbody></table>
<table><caption>**Table S2 :** statistics of assembled transcriptomes for individual samples and meta-transcriptomes. The assembly is based on paired-ends sequencing (2 x 75 bp) and the number of raw sequences corresponds to the number of pairs. Contigs indicates the number of contigs for each assembly, with the corresponding N50 and L50. The Lg columns corresponds to the contigs length in bp, with the sum, minimum, mean, median and maximum of Lg. The last two lanes refer to the meta-transcriptome obtained from all individual transcriptomes with or without potential Symbiodiniaceae sequences. See main text for details.</caption><thead><tr><th>Genus</th><th>Species</th><th>Sample</th><th>raw sequences</th><th>contigs</th><th>N50</th><th>L50</th><th>Lg sum</th><th>Lg min</th><th>Lg mean</th><th>Lg median</th><th>Lg max</th></tr></thead><tbody><tr><td>Eunicella</td><td>cavolini</td><td>e-cavol-anb-a</td><td>21432997</td><td>33627</td><td>1978</td><td>7240</td><td>46041698</td><td>201</td><td>1369.19</td><td>1023</td><td>13533</td></tr><tr><td>Eunicella</td><td>cavolini</td><td>e-cavol-anb-d</td><td>20070761</td><td>33624</td><td>2002</td><td>7025</td><td>46288676</td><td>201</td><td>1376.66</td><td>1004</td><td>24422</td></tr><tr><td>Eunicella</td><td>cavolini</td><td>e-cavol-som-a</td><td>22986734</td><td>43541</td><td>1757</td><td>8987</td><td>52113269</td><td>201</td><td>1196.88</td><td>840</td><td>24228</td></tr><tr><td>Eunicella</td><td>cavolini</td><td>e-cavol-vil-a</td><td>31846763</td><td>36908</td><td>2056</td><td>7687</td><td>51709228</td><td>201</td><td>1401.03</td><td>1018</td><td>18381</td></tr><tr><td>Eunicella</td><td>cavolini</td><td>e-cavol-vil-b</td><td>28751407</td><td>34961</td><td>2044</td><td>7470</td><td>48971519</td><td>201</td><td>1400.75</td><td>1040</td><td>19504</td></tr><tr><td>Eunicella</td><td>hybrid</td><td>e-hybri-jpb-a</td><td>34392918</td><td>39407</td><td>2031</td><td>8098</td><td>54290945</td><td>201</td><td>1377.7</td><td>998</td><td>25557</td></tr><tr><td>Eunicella</td><td>hybrid</td><td>e-hybri-mfn-a</td><td>44256795</td><td>40762</td><td>2081</td><td>8451</td><td>57794280</td><td>201</td><td>1417.85</td><td>1039</td><td>25573</td></tr><tr><td>Eunicella</td><td>hybrid</td><td>e-hybri-mfn-b</td><td>34705411</td><td>39672</td><td>2046</td><td>8040</td><td>54738734</td><td>201</td><td>1379.78</td><td>981.5</td><td>16650</td></tr><tr><td>Eunicella</td><td>hybrid</td><td>e-hybri-mfn-e</td><td>36536647</td><td>39532</td><td>2038</td><td>8090</td><td>54685655</td><td>201</td><td>1383.33</td><td>995.5</td><td>25578</td></tr><tr><td>Eunicella</td><td>singularis</td><td>e-singu-ban-a</td><td>44325669</td><td>45364</td><td>1919</td><td>9379</td><td>58576839</td><td>201</td><td>1291.26</td><td>928</td><td>28882</td></tr><tr><td>Eunicella</td><td>singularis</td><td>e-singu-ban-b</td><td>33184944</td><td>38095</td><td>1930</td><td>8114</td><td>50868966</td><td>201</td><td>1335.32</td><td>987</td><td>20211</td></tr><tr><td>Eunicella</td><td>singularis</td><td>e-singu-ban-c</td><td>46271612</td><td>43821</td><td>2023</td><td>9132</td><td>60512898</td><td>201</td><td>1380.91</td><td>1007</td><td>21714</td></tr><tr><td>Eunicella</td><td>singularis</td><td>e-singu-cav-a</td><td>48947180</td><td>51120</td><td>1967</td><td>10031</td><td>65261049</td><td>201</td><td>1276.62</td><td>868</td><td>22527</td></tr><tr><td>Eunicella</td><td>singularis</td><td>e-singu-mfn-a</td><td>52588076</td><td>70114</td><td>1761</td><td>13336</td><td>79649263</td><td>201</td><td>1136</td><td>739</td><td>16808</td></tr><tr><td>Eunicella</td><td>singularis</td><td>e-singu-mfn-b</td><td>43713977</td><td>55035</td><td>1894</td><td>10583</td><td>67120524</td><td>201</td><td>1219.6</td><td>808</td><td>21143</td></tr><tr><td>Eunicella</td><td>singularis</td><td>e-singu-sdo-a</td><td>37444166</td><td>55928</td><td>1741</td><td>10464</td><td>62326140</td><td>201</td><td>1114.4</td><td>715</td><td>16387</td></tr><tr><td>Eunicella</td><td>singularis</td><td>e-singu-sdo-b</td><td>39266148</td><td>72419</td><td>1652</td><td>13837</td><td>78950323</td><td>201</td><td>1090.19</td><td>715</td><td>24245</td></tr><tr><td>Eunicella</td><td>verrucosa</td><td>e-verru-ros-a</td><td>19398629</td><td>31195</td><td>1936</td><td>6630</td><td>41727111</td><td>201</td><td>1337.62</td><td>981</td><td>16974</td></tr><tr><td>Eunicella</td><td>verrucosa</td><td>e-verru-ros-c</td><td>20495748</td><td>31526</td><td>1968</td><td>6729</td><td>42779660</td><td>201</td><td>1356.96</td><td>1005</td><td>16663</td></tr><tr><td>Eunicella</td><td>verrucosa</td><td>e-verru-som-a</td><td>23332185</td><td>33133</td><td>2005</td><td>6944</td><td>45531674</td><td>201</td><td>1374.21</td><td>1005</td><td>25577</td></tr><tr><th colspan="12">Meta transcriptome</th></tr><tr><td></td><td></td><td></td><td></td><td>number of contigs</td><td>retained contigs</td><td>N50</td><td>L50</td><td>Lg sum</td><td>Lg min</td><td>Lg mean</td><td>Lg median</td><td>Lg max</td></tr><tr><td></td><td></td><td>meta</td><td>891354</td><td>68386</td><td>2144</td><td>14309</td><td>9</td><td>201</td><td>1500,62</td><td>1098</td><td>28882</td></tr><tr><td></td><td></td><td>meta no</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td></td><td></td><td>Symb</td><td>300085</td><td>59697</td><td>1975</td><td>12316</td><td>80903965</td><td>201</td><td>1355,24</td><td>967</td><td>25577</td></tr></tbody></table>
<table><caption>**Table S2 :** statistics of assembled transcriptomes for individual samples and meta-transcriptomes. The assembly is based on paired-ends sequencing (2 x 75 bp) and the number of raw sequences corresponds to the number of pairs. Contigs indicates the number of contigs for each assembly, with the corresponding N50 and L50. The Lg columns corresponds to the contigs length in bp, with the sum, minimum, mean, median and maximum of Lg. The last two lanes refer to the meta-transcriptome obtained from all individual transcriptomes with or without potential Symbiodiniaceae sequences. See main text for details.</caption><thead><tr><th>Genus</th><th>Species</th><th>Sample</th><th>raw sequences</th><th>contigs</th><th>N50</th><th>L50</th><th>Lg sum</th><th>Lg min</th><th>Lg mean</th><th>Lg median</th><th>Lg max</th></tr></thead><tbody><tr><td>Eunicella</td><td>cavolini</td><td>e-cavol-anb-a</td><td>21432997</td><td>33627</td><td>1978</td><td>7240</td><td>46041698</td><td>201</td><td>1369.19</td><td>1023</td><td>13533</td></tr><tr><td>Eunicella</td><td>cavolini</td><td>e-cavol-anb-d</td><td>20070761</td><td>33624</td><td>2002</td><td>7025</td><td>46288676</td><td>201</td><td>1376.66</td><td>1004</td><td>24422</td></tr><tr><td>Eunicella</td><td>cavolini</td><td>e-cavol-som-a</td><td>22986734</td><td>43541</td><td>1757</td><td>8987</td><td>52113269</td><td>201</td><td>1196.88</td><td>840</td><td>24228</td></tr><tr><td>Eunicella</td><td>cavolini</td><td>e-cavol-vil-a</td><td>31846763</td><td>36908</td><td>2056</td><td>7687</td><td>51709228</td><td>201</td><td>1401.03</td><td>1018</td><td>18381</td></tr><tr><td>Eunicella</td><td>cavolini</td><td>e-cavol-vil-b</td><td>28751407</td><td>34961</td><td>2044</td><td>7470</td><td>48971519</td><td>201</td><td>1400.75</td><td>1040</td><td>19504</td></tr><tr><td>Eunicella</td><td>hybrid</td><td>e-hybri-jpb-a</td><td>34392918</td><td>39407</td><td>2031</td><td>8098</td><td>54290945</td><td>201</td><td>1377.7</td><td>998</td><td>25557</td></tr><tr><td>Eunicella</td><td>hybrid</td><td>e-hybri-mfn-a</td><td>44256795</td><td>40762</td><td>2081</td><td>8451</td><td>57794280</td><td>201</td><td>1417.85</td><td>1039</td><td>25573</td></tr><tr><td>Eunicella</td><td>hybrid</td><td>e-hybri-mfn-b</td><td>34705411</td><td>39672</td><td>2046</td><td>8040</td><td>54738734</td><td>201</td><td>1379.78</td><td>981.5</td><td>16650</td></tr><tr><td>Eunicella</td><td>hybrid</td><td>e-hybri-mfn-e</td><td>36536647</td><td>39532</td><td>2038</td><td>8090</td><td>54685655</td><td>201</td><td>1383.33</td><td>995.5</td><td>25578</td></tr><tr><td>Eunicella</td><td>singularis</td><td>e-singu-ban-a</td><td>44325669</td><td>45364</td><td>1919</td><td>9379</td><td>58576839</td><td>201</td><td>1291.26</td><td>928</td><td>28882</td></tr><tr><td>Eunicella</td><td>singularis</td><td>e-singu-ban-b</td><td>33184944</td><td>38095</td><td>1930</td><td>8114</td><td>50868966</td><td>201</td><td>1335.32</td><td>987</td><td>20211</td></tr><tr><td>Eunicella</td><td>singularis</td><td>e-singu-ban-c</td><td>46271612</td><td>43821</td><td>2023</td><td>9132</td><td>60512898</td><td>201</td><td>1380.91</td><td>1007</td><td>21714</td></tr><tr><td>Eunicella</td><td>singularis</td><td>e-singu-cav-a</td><td>48947180</td><td>51120</td><td>1967</td><td>10031</td><td>65261049</td><td>201</td><td>1276.62</td><td>868</td><td>22527</td></tr><tr><td>Eunicella</td><td>singularis</td><td>e-singu-mfn-a</td><td>52588076</td><td>70114</td><td>1761</td><td>13336</td><td>79649263</td><td>201</td><td>1136</td><td>739</td><td>16808</td></tr><tr><td>Eunicella</td><td>singularis</td><td>e-singu-mfn-b</td><td>43713977</td><td>55035</td><td>1894</td><td>10583</td><td>67120524</td><td>201</td><td>1219.6</td><td>808</td><td>21143</td></tr><tr><td>Eunicella</td><td>singularis</td><td>e-singu-sdo-a</td><td>37444166</td><td>55928</td><td>1741</td><td>10464</td><td>62326140</td><td>201</td><td>1114.4</td><td>715</td><td>16387</td></tr><tr><td>Eunicella</td><td>singularis</td><td>e-singu-sdo-b</td><td>39266148</td><td>72419</td><td>1652</td><td>13837</td><td>78950323</td><td>201</td><td>1090.19</td><td>715</td><td>24245</td></tr><tr><td>Eunicella</td><td>verrucosa</td><td>e-verru-ros-a</td><td>19398629</td><td>31195</td><td>1936</td><td>6630</td><td>41727111</td><td>201</td><td>1337.62</td><td>981</td><td>16974</td></tr><tr><td>Eunicella</td><td>verrucosa</td><td>e-verru-ros-c</td><td>20495748</td><td>31526</td><td>1968</td><td>6729</td><td>42779660</td><td>201</td><td>1356.96</td><td>1005</td><td>16663</td></tr><tr><td>Eunicella</td><td>verrucosa</td><td>e-verru-som-a</td><td>23332185</td><td>33133</td><td>2005</td><td>6944</td><td>45531674</td><td>201</td><td>1374.21</td><td>1005</td><td>25577</td></tr><tr><th colspan="12" style="text-align: left;">Meta transcriptome</th></tr><tr><th></th><th></th><th></th><th>number of contigs</th><th>retained contigs</th><th>N50</th><th>L50</th><th>Lg sum</th><th>Lg min</th><th>Lg mean</th><th>Lg median</th><th>Lg max</th></tr><tr><td>Eunicella</td><td>meta</td><td></td><td>891354</td><td>68386</td><td>2144</td><td>14309</td><td>9</td><td>201</td><td>1500,62</td><td>1098</td><td>28882</td></tr><tr><td>Eunicella</td><td>meta no</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Eunicella</td><td>Symb</td><td></td><td>300085</td><td>59697</td><td>1975</td><td>12316</td><td>80903965</td><td>201</td><td>1355,24</td><td>967</td><td>25577</td></tr></tbody></table>
<table><caption>**Table S3:** summary of the different datasets; the first four datasets include variable and non variable sites (all sites), while the "polymorphic sites" and the "1% SNPs" datasets only consider SNPs, i.e. variable sites. For the "all" datasets we indicate the number of contigs and the number of sites retained from reads2snp. See main text for details</caption><thead><tr><th>dataset</th><th>samples</th><th>sites</th><th>number of contigs / SNPs</th><th>analyses</th></tr></thead><tbody><tr><td>all sites</td><td>all</td><td>all from reads2snp</td><td>61500 contigs / 101516577 sites</td><td>build SNPs datasets</td></tr><tr><td>all-CS</td><td>cavolini<br>singularis</td><td>/ all from reads2snp</td><td>61947 contigs / 101515803 sites</td><td>speciation scenarios with DILS</td></tr><tr><td>all-CV</td><td>cavolini<br>verrucosa</td><td>/ all from reads2snp</td><td>59702 contigs / 100704015 sites</td><td>speciation scenarios with DILS</td></tr><tr><td>all-SV</td><td>singularis<br>verrucosa</td><td>/ all from reads2snp</td><td>61373 contigs / 101444729 sites</td><td>speciation scenarios with DILS</td></tr><tr><td>polymorphic sites</td><td>all</td><td>polymorphic sites ; no 31369 SNPs<br>missing data</td><td>F<sub>ST</sub>, LEA, PCA</td><td></td></tr><tr><td>1 % SNPs</td><td>without<br>verrucosa</td><td>polymorphic sites ; no 326 SNPs<br>missing data ; 1 %<br>highest F<sub>ST</sub> cavolini /<br>singularis</td><td>newhybrids</td><td></td></tr></tbody></table>
<table><caption>**Table S3 :** summary of the different datasets; the first four datasets include variable and non variable sites (all sites), while the "polymorphic sites" and the "1% SNPs" datasets only consider SNPs, i.e. variable sites. For the "all" datasets we indicate the number of contigs and the number of sites retained from reads2snp. See main text for details</caption><thead><tr><th>dataset</th><th>samples</th><th>sites</th><th>number of contigs / SNPs</th><th>analyses</th></tr></thead><tbody><tr><td>all sites</td><td>all</td><td>all from reads2snp</td><td>61500 contigs / 101516577 sites</td><td>build SNPs datasets</td></tr><tr><td>all-CS</td><td>cavolini<br>singularis</td><td>/ all from reads2snp</td><td>61947 contigs / 101515803 sites</td><td>speciation scenarios with DILS</td></tr><tr><td>all-CV</td><td>cavolini<br>verrucosa</td><td>/ all from reads2snp</td><td>59702 contigs / 100704015 sites</td><td>speciation scenarios with DILS</td></tr><tr><td>all-SV</td><td>singularis<br>verrucosa</td><td>/ all from reads2snp</td><td>61373 contigs / 101444729 sites</td><td>speciation scenarios with DILS</td></tr><tr><td>polymorphic sites</td><td>all</td><td>polymorphic sites ; no 31369 SNPs<br>missing data</td><td>F<sub>ST</sub>, LEA, PCA</td><td></td></tr><tr><td>1 % SNPs</td><td>without<br>verrucosa</td><td>polymorphic sites ; no 326 SNPs<br>missing data ; 1 %<br>highest F<sub>ST</sub> cavolini /<br>singularis</td><td>newhybrids</td><td></td></tr></tbody></table>
<table><caption>**Table S4:** parameters used in the DILS analyses: Max_NA : maximum proportion of missing data ; Lmin : minimum sequence length per gene ; nMin : minimum number of sequences per gene and per species ; jSFS : use of joint Site Frequency Spectrum as an additional set of summary statistics ; constant / variable : consider constant or variable population size ; minimum and maximum values for the following priors : Tsplit : time of split, Ne : population size, M : migration rate. All other priors were kept at default values. For all analyses we used the option for coding regions, we didn't use any outgroup, we used the bimodal model for barriers, and the "normal" computation mode. The last column indicates the code used to describe the corresponding analysis in the text. The ranges of prior were chosen after preliminary analyses where we analysed the goodness of fit of the data to the models and priors. We used a mutation rate of 3.10⁻⁹.</caption><thead><tr><th>dataset</th><th>max_NA</th><th>Lmin</th><th>nMin</th><th>jSFS</th><th>Tsplit</th><th>Ne</th><th>M</th></tr></thead><tbody><tr><td>all-CS</td><td>0.1</td><td>30</td><td>10</td><td>yes</td><td>100 - 2 000 000</td><td>100 - 2 000 000</td><td>0-30</td></tr><tr><td>all-CV</td><td>0.1</td><td>30</td><td>6</td><td>yes</td><td>100 - 2 000 000</td><td>100 - 2 000 000</td><td>0-30</td></tr><tr><td>all-SV</td><td>0.1</td><td>30</td><td>6</td><td>yes</td><td>100 - 2 000 000</td><td>100 - 2 000 000</td><td>0-30</td></tr></tbody></table>
<table><caption>**Table S4:** parameters used in the DILS analyses: Max_NA : maximum proportion of missing data ; Lmin : minimum sequence length per gene ; nMin : minimum number of sequences per gene and per species ; jSFS : use of joint Site Frequency Spectrum as an additional set of summary statistics ; constant / variable : consider constant or variable population size ; minimum and maximum values for the following priors : Tsplit : time of split, Ne : population size, M : migration rate. All other priors were kept at default values. For all analyses we used the option for coding regions, we didn't use any outgroup, we used the bimodal model for barriers, and the "normal" computation mode. The last column indicates the code used to describe the corresponding analysis in the text. The ranges of prior were chosen after preliminary analyses where we analysed the goodness of fit of the data to the models and priors. We used a mutation rate of 3.10⁻⁹.</caption><thead><tr><th>dataset</th><th>max_NA</th><th>Lmin</th><th>nMin</th><th>jSFS</th><th>Tsplit</th><th>Ne</th><th>M</th></tr></thead><tbody><tr><td>all-CS</td><td>0.1</td><td>30</td><td>10</td><td>yes</td><td>100 - 2 000 000</td><td>100 - 2 000 000</td><td>0-30</td></tr><tr><td>all-CV</td><td>0.1</td><td>30</td><td>6</td><td>yes</td><td>100 - 2 000 000</td><td>100 - 2 000 000</td><td>0-30</td></tr><tr><td>all-SV</td><td>0.1</td><td>30</td><td>6</td><td>yes</td><td>100 - 2 000 000</td><td>100 - 2 000 000</td><td>0-30</td></tr></tbody></table>
<table><caption>**Table S5:** frequency of Symbiodiniaceae sequences in the individual transcriptomes on the basis i) of the proportion of raw reads mapped on the Symbiodiniaceae transcriptome, and ii) on the proportion of sequences in individual transcriptomes following the BLAT analysis. "meta" indicate the meta-transcriptome assembly based on all samples. See Table S2 for the codes of samples.</caption><thead><tr><th>Sample</th><th>Species</th><th>Raw reads</th><th>Transcriptome</th></tr></thead><tbody><tr><td>e-cavol-anb-a</td><td><i>E. cavolini</i></td><td>0.0171</td><td>0.00305</td></tr><tr><td>e-cavol-anb-d</td><td><i>E. cavolini</i></td><td>0.0087</td><td>0.00268</td></tr><tr><td>e-cavol-som-a</td><td><i>E. cavolini</i></td><td>0.0087</td><td>0.00426</td></tr><tr><td>e-cavol-vil-a</td><td><i>E. cavolini</i></td><td>0.0184</td><td>0.00350</td></tr><tr><td>e-cavol-vil-b</td><td><i>E. cavolini</i></td><td>0.0255</td><td>0.00333</td></tr><tr><td>e-hybri-jpb-a</td><td>hybrid</td><td>0.0123</td><td>0.00262</td></tr><tr><td>e-hybri-mfn-a</td><td>hybrid</td><td>0.0076</td><td>0.00321</td></tr><tr><td>e-hybri-mfn-b</td><td>hybrid</td><td>0.0079</td><td>0.00270</td></tr><tr><td>e-hybri-mfn-e</td><td>hybrid</td><td>0.0162</td><td>0.00302</td></tr><tr><td>e-singu-ban-a</td><td><i>E. singularis</i></td><td>0.0192</td><td>0.00675</td></tr><tr><td>e-singu-ban-b</td><td><i>E. singularis</i></td><td>0.0140</td><td>0.00410</td></tr><tr><td>e-singu-ban-c</td><td><i>E. singularis</i></td><td>0.0080</td><td>0.00647</td></tr><tr><td>e-singu-cav-a</td><td><i>E. singularis</i></td><td>0.0261</td><td>0.02263</td></tr><tr><td>e-singu-mfn-a</td><td><i>E. singularis</i></td><td>0.0233</td><td>0.03419</td></tr><tr><td>e-singu-mfn-b</td><td><i>E. singularis</i></td><td>0.0129</td><td>0.02745</td></tr><tr><td>e-singu-sdo-a</td><td><i>E. singularis</i></td><td>0.0158</td><td>0.03644</td></tr><tr><td>e-singu-sdo-b</td><td><i>E. singularis</i></td><td>0.0207</td><td>0.03686</td></tr><tr><td>e-verru-ros-a</td><td><i>E. verrucosa</i></td><td>0.0075</td><td>0.00276</td></tr><tr><td>e-verru-ros-c</td><td><i>E. verrucosa</i></td><td>0.0082</td><td>0.00282</td></tr><tr><td>e-verru-som-a</td><td><i>E. verrucosa</i></td><td>0.0098</td><td>0.00279</td></tr><tr><td>meta</td><td></td><td></td><td>0.01393</td></tr></tbody></table>
<table><caption>**Table S5:** frequency of Symbiodiniaceae sequences in the individual transcriptomes on the basis i) of the proportion of raw reads mapped on the Symbiodiniaceae transcriptome, and ii) on the proportion of sequences in individual transcriptomes following the BLAT analysis. "meta" indicate the meta-transcriptome assembly based on all samples. See Table S2 for the codes of samples.</caption><thead><tr><th>Sample</th><th>Species</th><th>Raw reads</th><th>Transcriptome</th></tr></thead><tbody><tr><td>e-cavol-anb-a</td><td><i>E. cavolini</i></td><td>0.0171</td><td>0.00305</td></tr><tr><td>e-cavol-anb-d</td><td><i>E. cavolini</i></td><td>0.0087</td><td>0.00268</td></tr><tr><td>e-cavol-som-a</td><td><i>E. cavolini</i></td><td>0.0087</td><td>0.00426</td></tr><tr><td>e-cavol-vil-a</td><td><i>E. cavolini</i></td><td>0.0184</td><td>0.00350</td></tr><tr><td>e-cavol-vil-b</td><td><i>E. cavolini</i></td><td>0.0255</td><td>0.00333</td></tr><tr><td>e-hybri-jpb-a</td><td>hybrid</td><td>0.0123</td><td>0.00262</td></tr><tr><td>e-hybri-mfn-a</td><td>hybrid</td><td>0.0076</td><td>0.00321</td></tr><tr><td>e-hybri-mfn-b</td><td>hybrid</td><td>0.0079</td><td>0.00270</td></tr><tr><td>e-hybri-mfn-e</td><td>hybrid</td><td>0.0162</td><td>0.00302</td></tr><tr><td>e-singu-ban-a</td><td><i>E. singularis</i></td><td>0.0192</td><td>0.00675</td></tr><tr><td>e-singu-ban-b</td><td><i>E. singularis</i></td><td>0.0140</td><td>0.00410</td></tr><tr><td>e-singu-ban-c</td><td><i>E. singularis</i></td><td>0.0080</td><td>0.00647</td></tr><tr><td>e-singu-cav-a</td><td><i>E. singularis</i></td><td>0.0261</td><td>0.02263</td></tr><tr><td>e-singu-mfn-a</td><td><i>E. singularis</i></td><td>0.0233</td><td>0.03419</td></tr><tr><td>e-singu-mfn-b</td><td><i>E. singularis</i></td><td>0.0129</td><td>0.02745</td></tr><tr><td>e-singu-sdo-a</td><td><i>E. singularis</i></td><td>0.0158</td><td>0.03644</td></tr><tr><td>e-singu-sdo-b</td><td><i>E. singularis</i></td><td>0.0207</td><td>0.03686</td></tr><tr><td>e-verru-ros-a</td><td><i>E. verrucosa</i></td><td>0.0075</td><td>0.00276</td></tr><tr><td>e-verru-ros-c</td><td><i>E. verrucosa</i></td><td>0.0082</td><td>0.00282</td></tr><tr><td>e-verru-som-a</td><td><i>E. verrucosa</i></td><td>0.0098</td><td>0.00279</td></tr><tr><td>meta</td><td></td><td></td><td>0.01393</td></tr></tbody></table>
<table><caption>**Table S6:** p-values of the Pairwise-Wilcoxon test on the frequency of Symbiodiniaceae. A) on the basis of read counts with Salmon; B) on the proportion of assembled sequences with the BLAT analysis</caption><thead><tr><th></th><th>E. cavolini</th><th>hybrids</th><th>E. singularis</th></tr></thead><tbody><tr><td><strong>A)</strong></td><td></td><td></td><td></td></tr><tr><td>hybrids</td><td>0.69</td><td></td><td></td></tr><tr><td>E. singularis</td><td>0.69</td><td>0.36</td><td></td></tr><tr><td>E. verrucosa</td><td>0.57</td><td>0.69</td><td>0.15</td></tr><tr><td><strong>B)</strong></td><td></td><td></td><td></td></tr><tr><td>hybrids</td><td>0.571</td><td></td><td></td></tr><tr><td>E. singularis</td><td>0.019</td><td>0.020</td><td></td></tr><tr><td>E. verrucosa</td><td>0.571</td><td>1</td><td>0.048</td></tr></tbody></table>
<table><caption>**Table S6:** p-values of the Pairwise-Wilcoxon test on the frequency of Symbiodiniaceae. A) on the basis of read counts with Salmon; B) on the proportion of assembled sequences with the BLAT analysis</caption><thead><tr><th></th><th>E. cavolini</th><th>hybrids</th><th>E. singularis</th></tr></thead><tbody><tr><td><strong>A)</strong></td><td></td><td></td><td></td></tr><tr><td>hybrids</td><td>0.69</td><td></td><td></td></tr><tr><td>E. singularis</td><td>0.69</td><td>0.36</td><td></td></tr><tr><td>E. verrucosa</td><td>0.57</td><td>0.69</td><td>0.15</td></tr><tr><td><strong>B)</strong></td><td></td><td></td><td></td></tr><tr><td>hybrids</td><td>0.571</td><td></td><td></td></tr><tr><td>E. singularis</td><td>0.019</td><td>0.020</td><td></td></tr><tr><td>E. verrucosa</td><td>0.571</td><td>1</td><td>0.048</td></tr></tbody></table>
<table><caption>**Table S7:** above diagonal: average net divergence estimated from DILS for the "all" pairwise datasets (the hybrids were not included in the DILS analysis); below diagonal pairwise F<sub>ST</sub> estimated from variable sites only ("polymorphic SNPs" dataset; see main text and Table S3 for details)</caption><thead><tr><th></th><th><em>E. cavolini</em></th><th>hybrids</th><th><em>E. singularis</em></th><th><em>E. verrucosa</em></th></tr></thead><tbody><tr><td><em>E. cavolini</em></td><td>-</td><td>-</td><td>0.0018</td><td>0.0067</td></tr><tr><td>hybrids</td><td>0.069</td><td>-</td><td>-</td><td>-</td></tr><tr><td><em>E. singularis</em></td><td>0.207</td><td>0.073</td><td>-</td><td>0.0070</td></tr><tr><td><em>E. verrucosa</em></td><td>0.432</td><td>0.456</td><td>0.529</td><td>-</td></tr></tbody></table>
<table><caption>**Table S7:** above diagonal: average net divergence estimated from DILS for the "all" pairwise datasets (the hybrids were not included in the DILS analysis); below diagonal pairwise F<sub>ST</sub> estimated from variable sites only ("polymorphic SNPs" dataset; see main text and Table S3 for details)</caption><thead><tr><th></th><th><em>E. cavolini</em></th><th>hybrids</th><th><em>E. singularis</em></th><th><em>E. verrucosa</em></th></tr></thead><tbody><tr><td><em>E. cavolini</em></td><td>-</td><td>-</td><td>0.0018</td><td>0.0067</td></tr><tr><td>hybrids</td><td>0.069</td><td>-</td><td>-</td><td>-</td></tr><tr><td><em>E. singularis</em></td><td>0.207</td><td>0.073</td><td>-</td><td>0.0070</td></tr><tr><td><em>E. verrucosa</em></td><td>0.432</td><td>0.456</td><td>0.529</td><td>-</td></tr></tbody></table>
**Table S8:** estimated parameters for the different evolutionary scenarios for the three pairwise comparisons. We present here the results of estimations for the optimized posterior with the random forests approach implemented in DILS. For each parameter we present the highest posterior density, with the median, and the lower and higher 2.5 % limits. Models: SC: secondary contact; AM: ancestral migration. Parameters: N: effective size; founders<sub>X</sub>: number of founder individuals in species X; T<sub>split</sub>: time of split at which the ancestral population subdivides in two populations; T<sub>SC</sub>: time of secondary contact; T<sub>AM</sub>: time of the end of gene flow for ancestral migration; T<sub>dem_X</sub>: time of demographic event for species X; M<sub>XY</sub>: introgression rate from Y to X. For all parameters, the subscripts indicate the species: A for ancestral, C for *E. cavolini*, S for *E. singularis*, and V for *E. verrucosa*. Times are given in generations, migration in numbers of migrants per generation.
**Table S8:** estimated parameters for the different evolutionary scenarios for the three pairwise comparisons. We present here the results of estimations for the optimized posterior with the random forests approach implemented in DILS. For each parameter we present the highest posterior density, with the median, and the lower and higher 2.5 % limits. Models: SC: secondary contact; AM: ancestral migration. Parameters: N: effective size; founders<sub>X</sub>: number of founder individuals in species X; T<sub>split</sub>: time of split at which the ancestral population subdivides in two populations; T<sub>SC</sub>: time of secondary contact; T<sub>AM</sub>: time of the end of gene flow for ancestral migration; T<sub>dem_X</sub>: time of demographic event for species X; M<sub>XY</sub>: introgression rate from Y to X. For all parameters, the subscripts indicate the species: A for ancestral, C for *E. cavolini*, S for *E. singularis*, and V for *E. verrucosa*. Times are given in generations, migration in numbers of migrants per generation.
A) comparison *E. cavolini* / *E. singularis*

<table><thead><tr><th></th><th>HPD 0.025</th><th>HPD median</th><th>HPD 0.0975</th></tr></thead><tbody><tr><td colspan="4"><strong>constant size, SC</strong></td></tr><tr><td>N<sub>C</sub></td><td>545985</td><td>633894</td><td>733842</td></tr><tr><td>N<sub>S</sub></td><td>168290</td><td>192199</td><td>225073</td></tr><tr><td>N<sub>A</sub></td><td>537403</td><td>581831</td><td>632310</td></tr><tr><td>T<sub>split</sub></td><td>336413</td><td>403273</td><td>476196</td></tr><tr><td>T<sub>SC</sub></td><td>51536</td><td>62039</td><td>71760</td></tr><tr><td>M<sub>CS</sub></td><td>12</td><td>15</td><td>17</td></tr><tr><td>M<sub>SC</sub></td><td>12</td><td>15</td><td>18</td></tr><tr><td colspan="4"><strong>variable size, SC</strong></td></tr><tr><td>N<sub>C</sub></td><td>531986</td><td>665965</td><td>875780</td></tr><tr><td>N<sub>S</sub></td><td>185826</td><td>222258</td><td>276889</td></tr><tr><td>N<sub>A</sub></td><td>515018</td><td>578861</td><td>640504</td></tr><tr><td>founders<sub>C</sub></td><td>0</td><td>1</td><td>1</td></tr><tr><td>founders<sub>S</sub></td><td>0</td><td>1</td><td>1</td></tr><tr><td>T<sub>dem_c</sub></td><td>250520</td><td>339056</td><td>418963</td></tr><tr><td>T<sub>dem_s</sub></td><td>245400</td><td>350132</td><td>454320</td></tr><tr><td>T<sub>split</sub></td><td>330907</td><td>434060</td><td>542765</td></tr><tr><td>T<sub>SC</sub></td><td>40560</td><td>57552</td><td>75405</td></tr><tr><td>M<sub>CS</sub></td><td>14</td><td>19</td><td>24</td></tr><tr><td>M<sub>SC</sub></td><td>8</td><td>12</td><td>16</td></tr></tbody></table>
A) comparison E. cavolini / E. singularis

<table><thead><tr><th></th><th>HPD 0.025</th><th>HPD median</th><th>HPD 0.0975</th></tr></thead><tbody><tr><td colspan="4"><strong>constant size, SC</strong></td></tr><tr><td>N<sub>C</sub></td><td>545985</td><td>633894</td><td>733842</td></tr><tr><td>N<sub>S</sub></td><td>168290</td><td>192199</td><td>225073</td></tr><tr><td>N<sub>A</sub></td><td>537403</td><td>581831</td><td>632310</td></tr><tr><td>T<sub>split</sub></td><td>336413</td><td>403273</td><td>476196</td></tr><tr><td>T<sub>SC</sub></td><td>51536</td><td>62039</td><td>71760</td></tr><tr><td>M<sub>CS</sub></td><td>12</td><td>15</td><td>17</td></tr><tr><td>M<sub>SC</sub></td><td>12</td><td>15</td><td>18</td></tr><tr><td colspan="4"><strong>variable size, SC</strong></td></tr><tr><td>N<sub>C</sub></td><td>531986</td><td>665965</td><td>875780</td></tr><tr><td>N<sub>S</sub></td><td>185826</td><td>222258</td><td>276889</td></tr><tr><td>N<sub>A</sub></td><td>515018</td><td>578861</td><td>640504</td></tr><tr><td>founders<sub>C</sub></td><td>0</td><td>1</td><td>1</td></tr><tr><td>founders<sub>S</sub></td><td>0</td><td>1</td><td>1</td></tr><tr><td>T<sub>dem_c</sub></td><td>250520</td><td>339056</td><td>418963</td></tr><tr><td>T<sub>dem_s</sub></td><td>245400</td><td>350132</td><td>454320</td></tr><tr><td>T<sub>split</sub></td><td>330907</td><td>434060</td><td>542765</td></tr><tr><td>T<sub>SC</sub></td><td>40560</td><td>57552</td><td>75405</td></tr><tr><td>M<sub>CS</sub></td><td>14</td><td>19</td><td>24</td></tr><tr><td>M<sub>SC</sub></td><td>8</td><td>12</td><td>16</td></tr></tbody></table>
## B) comparison *E. cavolini* / *E. verrucosa*

<table><thead><tr><th></th><th>HPD 0.025</th><th>HPD median</th><th>HPD 0.0975</th></tr></thead><tbody><tr><td colspan="4"><strong>constant size, AM</strong></td></tr><tr><td>N<sub>C</sub></td><td>630969</td><td>744556</td><td>875220</td></tr><tr><td>N<sub>V</sub></td><td>648850</td><td>755298</td><td>920095</td></tr><tr><td>N<sub>A</sub></td><td>698501</td><td>784512</td><td>879664</td></tr><tr><td>T<sub>split</sub></td><td>909392</td><td>1054488</td><td>1225792</td></tr><tr><td>T<sub>AM</sub></td><td>840920</td><td>991118</td><td>1147073</td></tr><tr><td>M<sub>CV</sub></td><td>4</td><td>6</td><td>7</td></tr><tr><td>M<sub>VC</sub></td><td>9</td><td>12</td><td>14</td></tr><tr><td colspan="4"><strong>variable size, AM</strong></td></tr><tr><td>N<sub>C</sub></td><td>777526</td><td>1099410</td><td>1694348</td></tr><tr><td>N<sub>V</sub></td><td>871210</td><td>1230360</td><td>1803231</td></tr><tr><td>N<sub>A</sub></td><td>692366</td><td>793880</td><td>930000</td></tr><tr><td>founders<sub>C</sub></td><td>0</td><td>1</td><td>1</td></tr><tr><td>founders<sub>V</sub></td><td>0</td><td>0</td><td>1</td></tr><tr><td>T<sub>dem_c</sub></td><td>237960</td><td>369260</td><td>496633</td></tr><tr><td>T<sub>dem_v</sub></td><td>335620</td><td>509096</td><td>679348</td></tr><tr><td>T<sub>split</sub></td><td>819714</td><td>1051517</td><td>1367074</td></tr><tr><td>T<sub>AM</sub></td><td>782210</td><td>930590</td><td>1104120</td></tr><tr><td>M<sub>CV</sub></td><td>7</td><td>12</td><td>16</td></tr><tr><td>M<sub>VC</sub></td><td>11</td><td>22</td><td>31</td></tr></tbody></table>
## B) comparison *E. cavolini* / *E. verrucosa*

<table><thead><tr><th></th><th>HPD 0.025</th><th>HPD median</th><th>HPD 0.0975</th></tr></thead><tbody><tr><td colspan="4"><strong>constant size, AM</strong></td></tr><tr><td>N<sub>C</sub></td><td>630969</td><td>744556</td><td>875220</td></tr><tr><td>N<sub>V</sub></td><td>648850</td><td>755298</td><td>920095</td></tr><tr><td>N<sub>A</sub></td><td>698501</td><td>784512</td><td>879664</td></tr><tr><td>T<sub>split</sub></td><td>909392</td><td>1054488</td><td>1225792</td></tr><tr><td>T<sub>AM</sub></td><td>840920</td><td>991118</td><td>1147073</td></tr><tr><td>M<sub>CV</sub></td><td>4</td><td>6</td><td>7</td></tr><tr><td>M<sub>VC</sub></td><td>9</td><td>12</td><td>14</td></tr><tr><td colspan="4"><strong>variable size, AM</strong></td></tr><tr><td>N<sub>C</sub></td><td>777526</td><td>1099410</td><td>1694348</td></tr><tr><td>N<sub>V</sub></td><td>871210</td><td>1230360</td><td>1803231</td></tr><tr><td>N<sub>A</sub></td><td>692366</td><td>793880</td><td>930000</td></tr><tr><td>founders<sub>C</sub></td><td>0</td><td>1</td><td>1</td></tr><tr><td>founders<sub>V</sub></td><td>0</td><td>0</td><td>1</td></tr><tr><td>T<sub>dem_c</sub></td><td>237960</td><td>369260</td><td>496633</td></tr><tr><td>T<sub>dem_v</sub></td><td>335620</td><td>509096</td><td>679348</td></tr><tr><td>T<sub>split</sub></td><td>819714</td><td>1051517</td><td>1367074</td></tr><tr><td>T<sub>AM</sub></td><td>782210</td><td>930590</td><td>1104120</td></tr><tr><td>M<sub>CV</sub></td><td>7</td><td>12</td><td>16</td></tr><tr><td>M<sub>VC</sub></td><td>11</td><td>22</td><td>31</td></tr></tbody></table>
C) comparison *E. singularis* / *E. verrucosa*

<table><thead><tr><th></th><th>HPD 0.025</th><th>HPD median</th><th>HPD 0.0975</th></tr></thead><tbody><tr><td colspan="4"><strong>constant size, AM</strong></td></tr><tr><td>N<sub>S</sub></td><td>263390</td><td>298162</td><td>336536</td></tr><tr><td>N<sub>V</sub></td><td>490519</td><td>592796</td><td>715930</td></tr><tr><td>N<sub>A</sub></td><td>632004</td><td>708517</td><td>790246</td></tr><tr><td>T<sub>split</sub></td><td>741840</td><td>899098</td><td>1091610</td></tr><tr><td>T<sub>AM</sub></td><td>698891</td><td>811827</td><td>934655</td></tr><tr><td>M<sub>SV</sub></td><td>10</td><td>14</td><td>17</td></tr><tr><td>M<sub>VS</sub></td><td>21</td><td>27</td><td>33</td></tr><tr><td colspan="4"><strong>variable size, AM</strong></td></tr><tr><td>N<sub>S</sub></td><td>281023</td><td>386388</td><td>494606</td></tr><tr><td>N<sub>V</sub></td><td>856542</td><td>1165039</td><td>1566087</td></tr><tr><td>N<sub>A</sub></td><td>592428</td><td>697054</td><td>797828</td></tr><tr><td>founders<sub>C</sub></td><td>0</td><td>0</td><td>0</td></tr><tr><td>founders<sub>V</sub></td><td>0</td><td>0</td><td>1</td></tr><tr><td>T<sub>dem_s</sub></td><td>166517</td><td>273546</td><td>374076</td></tr><tr><td>T<sub>dem_v</sub></td><td>226988</td><td>360174</td><td>493360</td></tr><tr><td>T<sub>split</sub></td><td>713634</td><td>926756</td><td>1207281</td></tr><tr><td>T<sub>AM</sub></td><td>454059</td><td>659458</td><td>858558</td></tr><tr><td>M<sub>SV</sub></td><td>3</td><td>4</td><td>6</td></tr><tr><td>M<sub>VS</sub></td><td>1</td><td>1</td><td>2</td></tr></tbody></table>
### C) comparison *E. singularis* / *E. verrucosa*

<table><thead><tr><th></th><th>HPD 0.025</th><th>HPD median</th><th>HPD 0.0975</th></tr></thead><tbody><tr><td colspan="4"><strong>constant size, AM</strong></td></tr><tr><td>N<sub>S</sub></td><td>263390</td><td>298162</td><td>336536</td></tr><tr><td>N<sub>V</sub></td><td>490519</td><td>592796</td><td>715930</td></tr><tr><td>N<sub>A</sub></td><td>632004</td><td>708517</td><td>790246</td></tr><tr><td>T<sub>split</sub></td><td>741840</td><td>899098</td><td>1091610</td></tr><tr><td>T<sub>AM</sub></td><td>698891</td><td>811827</td><td>934655</td></tr><tr><td>M<sub>SV</sub></td><td>10</td><td>14</td><td>17</td></tr><tr><td>M<sub>VS</sub></td><td>21</td><td>27</td><td>33</td></tr><tr><td colspan="4"><strong>variable size, AM</strong></td></tr><tr><td>N<sub>S</sub></td><td>281023</td><td>386388</td><td>494606</td></tr><tr><td>N<sub>V</sub></td><td>856542</td><td>1165039</td><td>1566087</td></tr><tr><td>N<sub>A</sub></td><td>592428</td><td>697054</td><td>797828</td></tr><tr><td>founders<sub>C</sub></td><td>0</td><td>0</td><td>0</td></tr><tr><td>founders<sub>V</sub></td><td>0</td><td>0</td><td>1</td></tr><tr><td>T<sub>dem_s</sub></td><td>166517</td><td>273546</td><td>374076</td></tr><tr><td>T<sub>dem_v</sub></td><td>226988</td><td>360174</td><td>493360</td></tr><tr><td>T<sub>split</sub></td><td>713634</td><td>926756</td><td>1207281</td></tr><tr><td>T<sub>AM</sub></td><td>454059</td><td>659458</td><td>858558</td></tr><tr><td>M<sub>SV</sub></td><td>3</td><td>4</td><td>6</td></tr><tr><td>M<sub>VS</sub></td><td>1</td><td>1</td><td>2</td></tr></tbody></table>
**Figure S1:** examples of morphological diversity in *Eunicella* species in the area of Marseille.

A) example of typical *E. cavolini* colonies (in the foreground)

<figure><img src="image_6.png" /><figcaption>*E. cavolini*</figcaption></figure>

B) example of a typical *E. singularis*

<figure><img src="image_7.png" /></figure>
**Figure S1:** examples of morphological diversity in *Eunicella* species in the area of Marseille.

A) example of typical *E. cavolini* colonies (in the foreground)

<figure><img src="image_7.png" /><figcaption>*E. cavolini*</figcaption></figure>

B) example of a typical *E. singularis*

<figure><img src="image_8.png" /></figure>
C) example of a typical *E. verrucosa*:

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

D) sampling with the presence of *E. cavolini* (EC), *E. singularis* (ES) and a potential hybrid (EH)

<figure><img src="image_9.png" /></figure>
C) example of a typical *E. verrucosa*:

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

D) sampling with the presence of *E. cavolini* (EC), *E. singularis* (ES) and a potential hybrid (EH)

<figure><img src="image_10.png" /></figure>
<figure><img src="image_10.png" /><figcaption>F) morphology of the colony EH-MFN-B with intermediate branching and color between *E. cavolini* and *E. singularis*.</figcaption></figure>

<figure><img src="image_11.png" /></figure>
<figure><img src="image_11.png" /><figcaption>F) morphology of the colony EH-MFN-B with intermediate branching and color between *E. cavolini* and *E. singularis*.</figcaption></figure>

<figure><img src="image_12.png" /></figure>
<figure><img src="image_12.png" /><figcaption>F) morphology of the colony EH-MFN-E with intermediate branching and color between *E. cavolini* and *E. singularis*</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>F) morphology of the colony EH-MFN-E with intermediate branching and color between *E. cavolini* and *E. singularis*</figcaption></figure>