Márcio

G Coelho

Gareth A Pearson
email: gpearson@ualg.pt

Joana R H Boavida

Diogo Paulo

Didier Aurelle

Sophie Arnaud-Haond

Daniel Gómez-Gras

Nathaniel Bensoussan

Paula López-Sendino

Carlo Cerrano

Márcio A G Coelho
email: macoelho@ualg.pt

| Diogo

| Didier Aurelle

| Daniel Gómez-Gras

Silvija Kipson

Tatjana Bakran-Petricioli

Eliana Ferretti

Cristina Linares

Joaquim Garrabou

Ester A Serrão

CIIMAR Jean-Baptiste Ledoux

/ Cimar

Not out of the Mediterranean: Atlantic populations of the gorgonian Paramuricea clavata are a separate sister species under further lineage diversification

Keywords: cryptic diversity, incomplete lineage sorting, Octocorallia, phylotranscriptomics, species delimitation

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

| INTRODUC TI ON

The species is the fundamental unit of analysis in ecology and evolution, from field and behavioral ecology to population genetics, phylogeography, and systematics [START_REF] De Queiroz | Ernst Mayr and the modern concept of species[END_REF]. Thus, species delimitation sets the tone for research questions in a myriad of fields, with critical implications for biodiversity assessment and conservation planning [START_REF] Hey | Understanding and confronting species uncertainty in biology and conservation[END_REF][START_REF] Pante | Species are hypotheses: Avoid connectivity assessments based on pillars of sand[END_REF][START_REF] Sukumaran | Multispecies coalescent delimits structure, not species[END_REF]. However, identifying species boundaries is challenging in many taxa for several reasons, including phenotypic plasticity [START_REF] Pica | Distribution and phenotypic variability of the Mediterranean gorgonian Paramuricea macrospina (cnidaria: Octocorallia)[END_REF][START_REF] Pino-Bodas | Phenotypical plasticity and homoplasy complicate species delimitation in the Cladonia gracilis group (Cladoniaceae, Ascomycota)[END_REF][START_REF] Sánchez | Phenotypic plasticity and morphological integration in a marine modular invertebrate[END_REF], lack of diagnostic morphological characters [START_REF] Erpenbeck | Speculation with spiculation? -Three independent gene fragments and biochemical characters versus morphology in demosponge higher classification[END_REF][START_REF] Mcfadden | Species boundaries in the absence of morphological, ecological or geographical differentiation in the Red Sea octocoral genus Ovabunda (Alcyonacea: Xeniidae)[END_REF], homoplasy [START_REF] Baker | Phylogeny, character evolution, and a new classification of the Calamoid palms[END_REF][START_REF] Morrow | Molecular phylogenies support homoplasy of multiple morphological characters used in the taxonomy of Heteroscleromorpha (porifera: Demospongiae)[END_REF][START_REF] Waeschenbach | First molecular estimate of cyclostome bryozoan phylogeny confirms extensive homoplasy among skeletal characters used in traditional taxonomy[END_REF], low-resolution of DNA barcoding markers [START_REF] Calderón | Evaluation of the utility of COI and ITS markers as tools for population genetic studies of temperate gorgonians[END_REF][START_REF] Huang | Slow mitochondrial COI sequence evolution at the base of the metazoan tree and its implications for DNA barcoding[END_REF][START_REF] Shearer | Slow mitochondrial DNA sequence evolution in the Anthozoa (cnidaria)[END_REF], incomplete lineage sorting (ILS) in recently diverged taxa [START_REF] Aurelle | Fuzzy species limits in Mediterranean gorgonians (cnidaria, Octocorallia): Inferences on speciation processes[END_REF][START_REF] Meyer | Disentangling incomplete lineage sorting and introgression to refine species-tree estimates for Lake Tanganyika cichlid fishes[END_REF][START_REF] Pollard | Widespread discordance of gene trees with species tree in Drosophila: Evidence for incomplete lineage sorting[END_REF], and introgressive hybridization [START_REF] Combosch | Trans-Pacific RAD-seq population genomics confirms introgressive hybridization in eastern Pacific Pocillopora corals[END_REF][START_REF] Padilla-García | The challenge of species delimitation in the diploid-polyploid complex Veronica subsection Pentasepalae[END_REF]reviewed in Harrison & Larson, 2014). The latter phenomenon is generally the result of complex biogeographic histories in which populations that diverged in allopatry come into secondary contact. Identifying barriers to gene flow and understanding how genetic differences accumulate over time is an inherently challenging task, particularly for marine systems in which these barriers are less obvious than for their terrestrial counterparts.

The Atlantic-Mediterranean transition provides an interesting geographic setting to study population differentiation and speciation in the marine realm [START_REF] Borsa | Infraspecific zoogeography of the Mediterranean: Population genetic analysis on sixteen Atlanto-Mediterranean species (fishes and invertebrates)[END_REF][START_REF] Patarnello | Pillars of Hercules: Is the Atlantic-Mediterranean transition a phylogeographical break[END_REF]. This narrow region has a complex biogeographical history, involving dramatic geological events that affected the entire Mediterranean Sea, including the Messinian Salinity Crisis and subsequent Zanclean flood (5.96 and 5.33 Ma ago, respectively;[START_REF] Garcia-Castellanos | Catastrophic flood of the Mediterranean after the Messinian salinity crisis[END_REF][START_REF] Garcia-Castellanos | Messinian salinity crisis regulated by competing tectonics and erosion at the Gibraltar arc[END_REF][START_REF] Krijgsman | Chronology, causes and progression of the Messinian salinity crisis[END_REF], as well as sea level fluctuations during glacial cycles [START_REF] Vacchi | Multiproxy assessment of Holocene relative sea-level changes in the western Mediterranean: Sea-level variability and improvements in the definition of the isostatic signal[END_REF]. These, in turn, have contributed to multiple episodes of species extinction-recolonization or allopatry-secondary contact between the Atlantic and Mediterranean basins, which have shaped contemporary patterns of genetic differentiation across this range [START_REF] Alberto | Genetic differentiation and secondary contact zone in the seagrass Cymodocea nodosa across the Mediterranean-Atlantic transition region[END_REF][START_REF] Duranton | The origin and remolding of genomic islands of differentiation in the European sea bass[END_REF]reviewed in Patarnello et al., 2007). Within the Mediterranean, several oceanographic discontinuities exist, among which is the Almeria-Oran Front (AOF), a large-scale density front at the eastern edge of the East Alboran gyre that extends southwards from Almeria (SE Spain) to Oran (Algeria) [START_REF] Millot | Circulation in the Western Mediterranean Sea[END_REF][START_REF] Tintore | A study of an intense density front in the eastern Alboran Sea: The Almeria-Oran front[END_REF]. The AOF constitutes a barrier to gene flow in many species, often resulting in complex phylogeographical patterns [START_REF] Patarnello | Pillars of Hercules: Is the Atlantic-Mediterranean transition a phylogeographical break[END_REF], Marinha e Ambiental, Universidade do Porto,Portugal. Email: jbaptiste.ledoux@gmail.com ability to understand species diversification. We investigated population differentiation and species limits in Atlantic (Iberian Peninsula) and Mediterranean lineages of the octocoral genus Paramuricea previously identified as P. clavata. We used a diverse set of molecular markers (microsatellites, RNA-seq derived single-copy orthologues [SCO] and mt-mutS [mitochondrial barcode]) at 49 locations. Clear segregation of Atlantic and Mediterranean lineages was found with all markers. Species-tree estimations based on SCO strongly supported these two clades as distinct, recently diverged sister species with incomplete lineage sorting, P. cf. grayi and P. clavata, respectively. Furthermore, a second putative (or ongoing) speciation event was detected in the Atlantic between two P. cf. grayi color morphotypes (yellow and purple) using SCO and supported by microsatellites. While segregating P. cf. grayi lineages showed considerable geographic structure, dominating circalittoral communities in southern (yellow) and western (purple) Portugal, their occurrence in sympatry at some localities suggests a degree of reproductive isolation. Overall, our results show that previous molecular and morphological studies have underestimated species diversity in Paramuricea occurring in the Iberian Peninsula, which has important implications for conservation planning. Finally, our findings validate the usefulness of phylotranscriptomics for resolving evolutionary relationships in octocorals.

| 3 of 21

COELHO et al.
particularly in taxa with restricted dispersal such as seagrasses [START_REF] Alberto | Genetic differentiation and secondary contact zone in the seagrass Cymodocea nodosa across the Mediterranean-Atlantic transition region[END_REF], sponges [START_REF] Dailianis | Genetic diversity of the imperilled bath sponge Spongia officinalis Linnaeus, 1759 across the Mediterranean Sea: Patterns of population differentiation and implications for taxonomy and conservation[END_REF][START_REF] Riesgo | Genetic diversity, connectivity and gene flow along the distribution of the emblematic Atlanto-Mediterranean sponge Petrosia ficiformis (Haplosclerida, Demospongiae)[END_REF], and corals [START_REF] Casado-Amezúa | Genetic assessment of population structure and connectivity in the threatened Mediterranean coral Astroides calycularis (Scleractinia, Dendrophylliidae) at different spatial scales[END_REF][START_REF] Mokhtar-Jamaï | From global to local genetic 20457758[END_REF]. Indeed, genetic discontinuities associated with the AOF are even reported for species with high dispersal potential such as the European sea bass (Dicentrarchus labrax; [START_REF] Duranton | The origin and remolding of genomic islands of differentiation in the European sea bass[END_REF][START_REF] Lemaire | Maintenance of genetic differentiation across a transition zone in the sea: Discordance between nuclear and cytoplasmic markers[END_REF][START_REF] Naciri | Genetic study of the Atlantic/Mediterranean transition in sea bass (Dicentrarchus labrax)[END_REF][START_REF] Tine | European sea bass genome and its variation provide insights into adaptation to euryhalinity and speciation[END_REF].

Branching octocorals of the genus Paramuricea (Octocorallia: Paramuriceidae) are important constituents of many shallow and deep-sea communities of benthic suspension feeders (e.g., corals, sponges, and bivalves) in the Atlantic Ocean and Mediterranean Sea, often forming extensive, dense assemblages supporting high levels of biodiversity [START_REF] Carpine | Les Gorgonaires de la Méditerranée[END_REF][START_REF] Cerrano | Gold coral (Savalia savaglia) and gorgonian forests enhance benthic biodiversity and ecosystem functioning in the mesophotic zone[END_REF][START_REF] Doughty | Insights into the population dynamics of the deep-sea coral genus Paramuricea in the Gulf of Mexico[END_REF][START_REF] Grasshoff | Die Gorgonarien des östlichen Nordatlantik und des Mittelmeeres III. Die Familie Paramuriceidae (Cnidaria, Anthozoa)[END_REF][START_REF] Grinyó | Diversity, distribution and population size structure of deep Mediterranean gorgonian assemblages (Menorca Channel, Western Mediterranean Sea)[END_REF][START_REF] Nestorowicz | Identifying habitats of conservation priority in the São Vicente submarine canyon in southwestern Portugal[END_REF][START_REF] Ponti | The understorey of gorgonian forests in mesophotic temperate reefs[END_REF][START_REF] Walting | Biology of deep-water octocorals[END_REF]. Currently, seven valid species are known to occur along the NE Atlantic and Mediterranean [START_REF] Carpine | Les Gorgonaires de la Méditerranée[END_REF][START_REF] Grasshoff | Die Gorgonarien des östlichen Nordatlantik und des Mittelmeeres III. Die Familie Paramuriceidae (Cnidaria, Anthozoa)[END_REF]WoRMS, 2022), with a further two putatively new species recently discovered after re-examination of natural history collections [START_REF] Sampaio | Natural history collections as a basis for sound biodiversity assessments: Plexauridae (Octocorallia, Holaxonia) of the Naturalis CANCAP and tyro Mauritania II expeditions[END_REF]. Two of these, Paramuricea clavata and Paramuricea macrospina, were historically regarded as endemic to the Mediterranean Sea [START_REF] Carpine | Les Gorgonaires de la Méditerranée[END_REF][START_REF] Grasshoff | Die Gorgonarien des östlichen Nordatlantik und des Mittelmeeres III. Die Familie Paramuriceidae (Cnidaria, Anthozoa)[END_REF][START_REF] Pica | Distribution and phenotypic variability of the Mediterranean gorgonian Paramuricea macrospina (cnidaria: Octocorallia)[END_REF] despite several recent studies extending the distribution ranges of P. clavata to adjacent Atlantic coastlines [START_REF] Boavida | Comparison of small remotely operated vehicles and diveroperated video of circalittoral benthos[END_REF]Cúrdia et al., 2012;[START_REF] Pilczynska | Genetic diversity increases with depth in red gorgonian populations of the Mediterranean Sea and the Atlantic Ocean[END_REF]Pilczynska, Cocito, et al., 2017) and of P. macrospina potentially to the Cape Verde archipelago [START_REF] Sampaio | Natural history collections as a basis for sound biodiversity assessments: Plexauridae (Octocorallia, Holaxonia) of the Naturalis CANCAP and tyro Mauritania II expeditions[END_REF].

The red gorgonian P. clavata is a major structuring species of Mediterranean coralligenous assemblages, generally restricted to low-light [START_REF] Ballesteros | Mediterranean coralligenous assemblages: A synthesis of present knowledge[END_REF], reaching densities well above 20 colonies m -2 [START_REF] Gori | Spatial distribution patterns of the gorgonians Eunicella singularis, Paramuricea clavata, and Leptogorgia sarmentosa (cape of Creus, northwestern Mediterranean Sea)[END_REF][START_REF] Linares | Size distribution, density and disturbance in two Mediterranean gorgonians: Paramuricea clavata and Eunicella singularis[END_REF][START_REF] Ponti | The understorey of gorgonian forests in mesophotic temperate reefs[END_REF]. The genetic structure of P. clavata in the Mediterranean shows sharp regional genetic breaks, strong differentiation at 10-100 s of meters and significant isolation by distance (IBD) at both regional and local scales [START_REF] Arizmendi-Mejía | Combining genetic and demographic data for the conservation of a Mediterranean marine habitat-forming species[END_REF][START_REF] Ledoux | Postglacial range expansion shaped the spatial genetic structure in a marine habitatforming species: Implications for conservation plans in the eastern Adriatic Sea[END_REF][START_REF] Mokhtar-Jamaï | From global to local genetic 20457758[END_REF][START_REF] Thoma | Molecular and morphological diversity of sea fans with emphasis on deep-sea octocorals of the order Alcyonacea Lamouroux[END_REF]. Spatial genetic structure in P. clavata is congruent with oceanographic barriers to gene flow (AOF: Mokhtar-Jamaï et al., 2011; seasonal cyclonic gyres along the Adriatic: [START_REF] Ledoux | Postglacial range expansion shaped the spatial genetic structure in a marine habitatforming species: Implications for conservation plans in the eastern Adriatic Sea[END_REF], and with the limited dispersal ability inherent to a reproductive strategy of larval brooding on maternal colonies (as observed in other brooding corals: [START_REF] Ayre | Genotypic diversity and gene flow in brooding and spawning corals along the Great Barrier Reef, Australia[END_REF][START_REF] Ledoux | Genetic survey of shallow populations of the Mediterranean red coral [Corallium rubrum (Linnaeus, 1758)]: New insights into evolutionary processes shaping nuclear diversity and implications for conservation[END_REF]Mokhtar-Jamaï et al., 2013;[START_REF] Smilansky | Fine-scale genetic structure in the surface brooding Caribbean octocoral, Antillogorgia elisabethae[END_REF][START_REF] Underwood | Multiple scales of genetic connectivity in a brooding coral on isolated reefs following catastrophic bleaching[END_REF]. Atlantic populations identified as P. clavata along the Portuguese coast are genetically differentiated from Mediterranean ones and also differentiated between western and southern Portugal [START_REF] Pilczynska | Genetic diversity increases with depth in red gorgonian populations of the Mediterranean Sea and the Atlantic Ocean[END_REF]Pilczynska, Cocito, et al., 2017). Interestingly, this differentiation between western and southern populations in Portugal appears to coincide with color morphotypes that dominate circalittoral coral gardens in each region (Pilczynska, Cocito, et al., 2017), with a yellow morphotype prevalent in the Algarve (southern Portugal) and a purple morphotype dominating communities in western Portugal. Although both morphotypes can occur in sympatry, our field surveys indicate that the purple morphotype is extremely rare in the south [START_REF] Dias | et al. coastal fisheries reveals rich coral habitats in southern Portugal[END_REF]; colonies reported as P. clavata) and vice-versa (J. Boavida and M. Coelho, personal observations). Despite these findings, the existence and exact location of the Atlantic-Mediterranean break, and the level of divergence across these regions (i.e., whether it represents within-species population structure or speciation processes) remains unknown in P. clavata due to incomplete geographic sampling and low marker resolution [START_REF] Mokhtar-Jamaï | From global to local genetic 20457758[END_REF][START_REF] Pilczynska | Genetic diversity increases with depth in red gorgonian populations of the Mediterranean Sea and the Atlantic Ocean[END_REF]Pilczynska, Cocito, et al., 2017). Indeed, standard DNA barcoding approaches for species delimitation are poorly suited to resolve closely related or recently diverged octocorals due to slow rates of mitochondrial gene evolution [START_REF] Calderón | Evaluation of the utility of COI and ITS markers as tools for population genetic studies of temperate gorgonians[END_REF][START_REF] Hellberg | No variation and low synonymous substitution rates in coral mtDNA despite high nuclear variation[END_REF][START_REF] Huang | Slow mitochondrial COI sequence evolution at the base of the metazoan tree and its implications for DNA barcoding[END_REF][START_REF] Mcfadden | Limitations of mitochondrial gene barcoding in Octocorallia[END_REF][START_REF] Shearer | Slow mitochondrial DNA sequence evolution in the Anthozoa (cnidaria)[END_REF]. These issues obscure the relationship between Mediterranean and Atlantic populations. For instance, Mediterranean populations of the threatened precious red coral Corallium rubrum typically exhibit strong genetic differentiation at small spatial scales [START_REF] Costantini | Fine-scale genetic structuring in Corallium rubrum: Evidence of inbreeding and limited effective larval dispersal[END_REF]Ledoux, Garrabou, et al., 2010;[START_REF] Ledoux | Genetic survey of shallow populations of the Mediterranean red coral [Corallium rubrum (Linnaeus, 1758)]: New insights into evolutionary processes shaping nuclear diversity and implications for conservation[END_REF].

However, the Atlantic-Mediterranean genetic break associated with the AOF is apparently lacking, presumably due to the extended influence of Atlantic populations in the Mediterranean gene pool, such as those recently rediscovered in SW Portugal [START_REF] Aurelle | Phylogeography of the red coral (Corallium rubrum): Inferences on the evolutionary history of a temperate gorgonian[END_REF][START_REF] Boavida | A well-kept treasure at depth: Precious red coral rediscovered in Atlantic deep coral gardens (SW Portugal) after 300 years[END_REF].

In addition to limited geographic sampling and marker resolution, the evolutionary relationships between P. clavata populations occurring across the Atlantic-Mediterranean transition have also been obscured by the potential occurrence in sympatry of congeneric species with intergrading morphologies that are difficult to discriminate.

For instance, the southern coast of Portugal is part of the distribution range of Paramuricea grayi, a species of deep littoral habitats historically described from the Lusitanian-Mauritanian East Atlantic to the Gulf of Guinea [START_REF] Grasshoff | Die Gorgonarien des östlichen Nordatlantik und des Mittelmeeres III. Die Familie Paramuriceidae (Cnidaria, Anthozoa)[END_REF][START_REF] Grasshoff | Die Flachwasser-Gorgonarien von Europa und Westafrika (Cnidaria, Anthozoa)[END_REF], but which has subsequently been described from Galicia and the Bay of Biscay (northern Iberian Peninsula) [START_REF] Poliseno | Historical biogeography and mitogenomics of two endemic Mediterranean gorgonians (Holaxonia, Plexauridae)[END_REF], among other locations in the NW Atlantic [START_REF] Thoma | Molecular and morphological diversity of sea fans with emphasis on deep-sea octocorals of the order Alcyonacea Lamouroux[END_REF][START_REF] Walting | Biology of deep-water octocorals[END_REF]. Importantly, recent phylogenetic reconstructions based on a mitochondrial marker and RAD-sequencing (RAD-seq) loci have shown that P. grayi is a sister species of P. clavata [START_REF] Poliseno | Historical biogeography and mitogenomics of two endemic Mediterranean gorgonians (Holaxonia, Plexauridae)[END_REF][START_REF] Quattrini | Phylogeography of Paramuricea: The role of depth and water mass in the evolution and distribution of deep-sea corals[END_REF].

The divergence between these two species could be either the result of a vicariance event in an Atlantic-Mediterranean ancestor [START_REF] Poliseno | Historical biogeography and mitogenomics of two endemic Mediterranean gorgonians (Holaxonia, Plexauridae)[END_REF]; or of divergence within the Mediterranean in which P. clavata diversified in mesophotic habitats from an ancestor with a broad depth range [START_REF] Quattrini | Phylogeography of Paramuricea: The role of depth and water mass in the evolution and distribution of deep-sea corals[END_REF]. Thus, it is possible that the Iberian Peninsula and adjacent areas form a contact zone between two recently diverged sister species (P. clavata in the Mediterranean and P. grayi in the NE Atlantic). Alternatively, the Atlantic populations of P. clavata reported from Portugal [START_REF] Boavida | Comparison of small remotely operated vehicles and diveroperated video of circalittoral benthos[END_REF]Cúrdia et al., 2012;[START_REF] Pilczynska | Genetic diversity increases with depth in red gorgonian populations of the Mediterranean Sea and the Atlantic Ocean[END_REF]Pilczynska, Cocito, et al., 2017) could be P. grayi misidentified as P. clavata. This highlights the need for further study of Paramuricea (and octocorals in general) along the Atlantic-Mediterranean transition to better understand the evolutionary forces governing population differentiation and speciation in this genus across these important biogeographic regions. The delimitation of spatially significant and unique genetic units within a species is also necessary for emerging conservation policies aiming to preserve intraspecific genetic diversity and the associated adaptive potential. Paramuricea clavata needs such efforts given its key structuring role in ecosystem functioning across the Mediterranean [START_REF] Rossi | Benthic-pelagic coupling: New perspectives in the animal forests[END_REF][START_REF] Rossi | The importance of food pulses in benthicpelagic coupling processes of passive suspension feeders[END_REF], particularly with the growing recognition of its vulnerability to human disturbances and environmental changes [START_REF] Cerrano | A catastrophic mass-mortality episode of gorgonians and other organisms in the Ligurian Sea (North-Western Mediterranean), summer 1999[END_REF][START_REF] Garrabou | Mass mortality in northwestern Mediterranean rocky benthic communities: Effects of the 2003 heat wave[END_REF][START_REF] Garrabou | Sliding toward the of Mediterranean coastal marine rocky ecosystems[END_REF][START_REF] Linares | Immediate and delayed effects of a mass mortality event on gorgonian population dynamics and benthic community structure in the NW Mediterranean Sea[END_REF]Linares, Coma, & Zabala, 2008). Resolving the nature of the divergence (i.e., population-or species-level) across the Atlantic and Mediterranean is a necessary first step.

The aim of this study was to test genetic boundaries and species affinities within the genus Paramuricea across the Atlantic-Mediterranean transition zone and assess the genetic diversity of the different lineages identified in these biogeographic regions. For this aim, the geographic scope includes the western and southern coasts of Portugal (Northeast Atlantic) and a large part of the distribution of P. clavata in the Mediterranean Sea, including the Adriatic and Aegean seas. We combined evidence from three types of molecular markers: (1) single-copy orthologs (obtained from RNA-seq data) for species delimitation and screening for cryptic species diversity;

(2) microsatellite loci to analyze spatial genetic diversity and structure both across and within the Atlantic and Mediterranean regions considering alternatively lineage and population levels; and (3) the mitochondrial gene mt-mutS (assembled from RNA-seq read data) to extend taxonomic coverage of our phylogenetic analysis for the genus. We present evidence showing that previous assignments of Atlantic populations to P. clavata are incongruent with genomic data, and the taxa present along those coasts should instead be regarded as the morphospecies P. cf. grayi.

| MATERIAL S AND ME THODS

2.1 | Phylogenomic inference and species delimitation using RNA-seq data

| Sample collection

Gorgonian samples used for the phylogenomic analyses were subsampled from colonies used as controls in thermo-tolerance experiments by [START_REF] Gómez-Gras | Exploring the response of a key Mediterranean gorgonian to heat stress across biological and spatial scales[END_REF], focusing on four populations: Balun (Kornati, Croatia), Altare (Portofino, Italy), and La Vaca (Medes, Catalonia) for P. clavata; and Baleeira (Sagres, Portugal) for P. cf. grayi (Figure 1 and Table S1). Tissue from each of 6-8 individuals per population was sampled for RNA-seq analysis (Table S1). The samples were preserved in RNAlater and stored at -80°C until further processing. In addition, four samples of P. cf. grayi were collected from two locations in southern (Tavira) and western (P39, Cape Espichel) Portugal. This included two specimens of the yellow color morph (1 from Tavira and 1 from P39) and two of the purple color morph from P39 (Table S1). These samples were flash frozen with liquid nitrogen upon collection and stored at -80°C until further processing.

| RNA extraction and sequencing

Total RNA was extracted from a 0.5 to 1-cm-long branchlet piece.

We coupled the TRIzol™ Plus RNA Purification Kit (Invitrogen) or the NZYol reagent (Nzytech) for tissue homogenization with the RNeasy Mini Kit (Qiagen) for RNA purification following the manufacturer's instructions. Residual DNA was removed using TurboDNAse (Invitrogen) or DNAse I (Qiagen). The quality, purity, and yield of the RNA was assessed using gel electrophoresis and a Qubit® Fluorometer and/or Nanodrop spectrophotometer. RNA-seq library preparation and sequencing was outsourced to BGI Tech Solutions in Hong Kong (100 bp paired-end DNBseq). A total of 34 individuals were processed for RNA-seq and the sequences deposited in NCBI's Sequence Read Archive (BioProject ID: PRJNA847883; BioSample accessions SAMN28899278-SAMN28899331).

| Sequence quality control and filtering

Raw reads were preprocessed by BGI to trim adaptor sequences and low-quality reads (<Q20) using filter_fq (https://github.com/ great fireb all/filter_fq). The quality of filtered sequence data was then evaluated with two complementary quality control tools, FastQC v0.11.6 (https://www.bioin forma tics.babra ham.ac.uk/proje cts/fastq c/) and FastqPuri [START_REF] Pérez-Rubio | FastqPuri: Highperformance preprocessing of RNA-seq data[END_REF]. Because these analyses detected the presence of overrepresented sequences corresponding to ribosomal RNAs (rRNAs), contaminant rRNA reads were filtered with FastqPuri using the bloom filter approach implemented in the function trimFilterPE (for details see Pérez-Rúbio et al., 2019 and the FastqPuri documentation page: https://github. com/jenge lmann/ Fastq Puri). The input bloom filter was constructed from nuclear rRNA sequences downloaded from NCBI GenBank for the top octocoral blastn hit of the contaminants, which corresponded to the octocoral Dendronephthya gigantea.

| De novo transcriptome assembly and completeness

Reference transcriptomes were assembled de novo with rnaSPAdes using default parameters [START_REF] Bushmanova | Spatial and depth-associated distribution patterns of shallow gorgonians in the Algarve coast (Portugal, NE Atlantic)[END_REF]. In total, we performed four independent transcriptome assemblies using the sequence data of all samples from each population, one for samples from Baleeira (P. cf. grayi) and three for samples from Balun, Altare and La Vaca (P. clavata). The assemblies were subjected to multiple steps of quality control and curation. First, the quality of each assem- S1. Sites sampled for RNA-seq are underlined, with populations that were also sampled for microsatellite genotyping shown in bold. PHL to GAB -all sites between Pharillons in the south of Marseille and Gabiniere in port-Cros (France) as listed in Table S1; PZU to GGL -All sites in northern Corsica (France). the assembled contigs based on evidence from the reads used to generate the assembly, allowing the filtration of "bad" (i.e., poorly assembled) contigs from the assembly [START_REF] Smith-Unna | TransRate: Reference-free quality assessment of de novo transcriptome assemblies[END_REF].

Second, using the "good" contigs obtained with TransRate, we identified putative open reading frames (ORFs) with FragGeneScan [START_REF] Rho | FragGeneScan: Predicting genes in short and error-prone reads[END_REF] and clustered the ORFs at 97% nucleotide sequence identity using VSEARCH [START_REF] Rognes | VSEARCH: A versatile open source tool for metagenomics[END_REF]. Finally, the contigs of the clustered assemblies were aligned against NCBI's nonredundant database for Anthozoa (Taxonomy ID: 6101) using DIAMOND in BLASTX mode [START_REF] Buchfink | Fast and sensitive protein alignment using DIAMOND[END_REF]. Contigs with a top query hit were retained in the final transcriptome assemblies, with the remainder discarded as potential contaminants (e.g., bacteria).

The completeness of the final, curated transcriptome assemblies was evaluated with BUSCO v2/v3 by screening the recovery of single-copy orthologs (queried as predicted proteins) that are expected to be present across higher taxonomic groups [START_REF] Simão | BUSCO: Assessing genome assembly and annotation completeness with single-copy orthologs[END_REF]. The analysis was performed against the Metazoa reference database (OrthoDB 9) using the online platform gVolante v1.2.1 available at https://gvola nte.riken.jp/ [START_REF] Nishimura | gVolante for standardizing completeness assessment of genome and transcriptome assemblies[END_REF].

Given that the transcriptome assembly for P. cf. grayi had the highest TransRate and BUSCO scores (Section 3), all subsequent analyses were conducted using this assembly as the reference.

| Phylogenetic analysis

The phylogenetic relationships among the samples were explored using 261 complete and unduplicated single-copy orthologues (≥180 bp) identified from BUSCO (N = 106) and from a reciprocal BLASTX search between the transcriptome assembly of P. cf. grayi and that of the octocoral Eunicella verrucosa (Coelho, M. A. G., Serrão, E. A., Pearson G. A., unpublished data; for additional information about the samples used see Table S1) using DIAMOND (N = 155).

For the reciprocal blast, we retained contigs with up to five query hits, subsequently selecting only those with a single sequence hit in both results tables (and thus putative single-copy orthologues).

Transcript .fasta sequences for orthologues were extracted for each sample from sorted .bam files after phasing (samtools phase) and variant calling (bcftools mpileup). Allele sequences were obtained from the resulting .vcf files with vcfutils.pl and seqkit fq2fa, aligned with mafft and trimmed with Gblocks in TranslatorX using a custom python script (https://github.com/cymon). The final set of loci used contained no gaps or missing data. The input .bam files were generated by mapping reads onto the reference transcriptome using the RSEM (v1.2.31) wrapper script and Bowtie2 [START_REF] Langmead | Fast gapped-read alignment with bowtie 2[END_REF][START_REF] Li | RSEM: Accurate transcript quantification from RNA-seq data with or without a reference genome[END_REF]. Locus alignments were assigned to "pseudoalleles" (i.e. alleles derived from in silico phasing of short read data stored in the alignment files) and analyzed by maximum likelihood in IQ-TREE 2 [START_REF] Minh | Role of evolutionary and ecological factors in the reproductive success and the spatial genetic structure of the temperate gorgonian Paramuricea clavata[END_REF]. A partitioned analysis was performed by fitting a separate evolutionary model of sequence evolution for each locus using ModelFinder and based on BIC score. Gene and site concordance factor (CF) analysis, which estimates the percentage of gene trees (gCF) or sites (sCF) which agree with the consensus species tree, was performed following [START_REF] Minh | New methods to calculate concordance factors for phylogenomic datasets[END_REF] and the recommendations at http:// www.iqtree.org/doc/. Species delimitation on the same multi-individual dataset as used for the species tree was performed using Species bOundry Delimitation using Astral method (SODA; [START_REF] Rabiee | SODA: Multi-locus species delimitation using quartet frequencies[END_REF] based on quartet tree analysis with ASTRAL-III (v5.7.3;[START_REF] Rabiee | Multi-allele species reconstruction using ASTRAL[END_REF]Zhang et al., 2018). The input ML gene trees from IQ-TREE 2 were used, and a species tree was estimated with ASTRAL-III. SODA uses the mapping of quartet tree constellations onto an ASTRAL species tree with all individuals to assess if/where coalescence of the tree topology is completely random. Where the coalescence is not random, species boundaries are defined based on the quartet frequencies [START_REF] Rabiee | SODA: Multi-locus species delimitation using quartet frequencies[END_REF].

In addition to SODA species delimitation, we reconstructed the species tree under an ML polymorphism-aware model (PoMo;[START_REF] De Maio | PoMo: An allele frequency-based approach for species tree estimation[END_REF]. PoMo is a statistically consistent method for species tree estimation that builds on nucleotide substitution models to incorporate population-level processes from allele frequency information [START_REF] Borges | Consistency and identifiability of the polymorphism-aware phylogenetic models[END_REF]. PoMo was run in IQ-TREE 2 under the HKY (Hasegawa-Kishino-Yano) substitution model, using individual allele count data from the 261 single-copy orthologue alignments. For the analysis, individual counts were assigned to a total of eight populations corresponding to sampling location and color morph. Branch support values were estimated with 1000 ultrafast bootstrap replicates. The species Paramuricea biscaya was used as an outgroup in all the analyses. We downloaded RNA-seq data for four samples from NCBI's SRA (Accessions: SRX4389727-SRX4389730; [START_REF] Deleo | Gene expression profiling reveals deep-sea coral response to the Deepwater Horizon oil spill[END_REF]. Sequence data were filtered and mapped as described above.

Finally, to extend the taxonomic coverage of the analysis to other Paramuricea spp., we conducted an ML tree reconstruction in IQ-TREE 2 with the mitochondrial gene mt-mutS. Sequence data were retrieved from NCBI Genbank and from .bam files generated here by mapping the RNA-seq reads onto the reference mitogenome of P. clavata using RSEM (Genbank Accession: NC_034749) for a subset of the samples (for additional information see Figure S1). New mt-mutS sequences were deposited in NCBI GenBank (Accession numbers: ON804207-ON804214). [START_REF] Boavida | Comparison of small remotely operated vehicles and diveroperated video of circalittoral benthos[END_REF][START_REF] Boavida | A well-kept treasure at depth: Precious red coral rediscovered in Atlantic deep coral gardens (SW Portugal) after 300 years[END_REF].

| Population differentiation based on microsatellite loci

| Sample collection

We extended the search along the Portuguese west coast between 30 and 70 m depth with the support of volunteer divers (through the citizen participation project Deep Reefs). Nine locations were sampled: Berlengas archipelago (2 sites 20 km apart), Cape Espichel (2 sites 400 m apart), Sines (1 site), Sagres (2 sites 300 m apart), Lagos (1 site), and Tavira (1 site; Table S1 and Figure 1). At each site, samples were collected haphazardly taking well-separated colonies (i.e., individuals) of the same height ( ͌ 30 cm height, adult colonies) to avoid sampling the same colony multiple times and to account for the very low rate of asexual reproduction observed in P. clavata/P. cf. grayi [START_REF] Coma | Reproduction and cycle of gonadal development in the Mediterranean gorgonian Paramuricea clavata[END_REF]Pilczynska, Boavida, et al., 2017). At each site, sample size varied from 12 (BAL) to 84 (CAT) individuals according to the species' density and decompression time constraints in deeper dives (Table S1). Samples consisted of 10 cm apical branches that were preserved in 96% ethanol at room temperature until DNA extraction. Two additional locations in the Mediterranean were sampled (VAC and PZU; Table S1) and combined with 35 existing sampling locations in the Mediterranean Sea described by [START_REF] Mokhtar-Jamaï | From global to local genetic 20457758[END_REF]. Overall, the dataset covered 46 locations, nine from the Atlantic and 37 from the Mediterranean, and included 1597 individuals (448 newly sampled; Figure 1 and Table S1). After quality filtering, the microsatellite dataset kept for downstream analyses included 1223 samples from 46 sites (filters: 23% of samples excluded due to missing data, sites with <5 samples and identical multilocus genotypes).

| Analysis of genetic diversity among lineages

Indices of genetic diversity were compared among the three lineages identified with SODA, P. clavata from the Mediterranean and two segregating lineages of P. cf. grayi coincident with color morphotype (yellow or purple) within the Atlantic (see Results). We compared the observed heterozygosity (H o ), gene diversity (H e ; [START_REF] Nei | Analysis of gene diversity in subdivided populations[END_REF], F IS [START_REF] Weir | Estimating F-statistics for the analysis of population structure[END_REF], allelic richness (Ar g ) and private allelic richness (Ap g ). The H o , H e and F IS were computed for each lineage with GENETIX 4.05 [START_REF] Belkhir | [END_REF], whereas Ar g and Ap g were estimated using the rarefaction method [START_REF] Petit | Identifying populations for conservation on the basis of genetic markers[END_REF] implemented in ADZE [START_REF] Szpiech | ADZE: A rarefaction approach for counting alleles private to combinations of populations[END_REF]. Statistical support for differences in Ar g and Ap g were tested with a Kruskal-Wallis rank sum test followed by pairwise Wilcoxon rank sum tests in R. Further details about these procedures can be found in Appendix S1.

| Genetic structure

A Bayesian clustering analyses was conducted with STRUCTURE 2.3.4 [START_REF] Pritchard | Inference of population structure using multilocus genotype data[END_REF] to evaluate the number of genetic clusters (K) from the individual genotypes without assumptions concerning population boundaries. Because of the unbalanced sampling [START_REF] Wang | The computer program structure for assigning individuals to populations: Easy to use but easier to misuse[END_REF], we considered the admixture model with uncorrelated allele frequencies, a separate ALPHA for each population and an initial ALPHA value = 0.022 (1/46 the number of sampling locations). Null homozygous genotypes were considered as missing data [START_REF] Falush | Inference of population structure using multilocus genotype data: Linked loci and correlated allele frequencies[END_REF][START_REF] Falush | Inference of population structure using multilocus genotype data: Dominant markers and null alleles[END_REF]. Ten independent runs were performed for each K ranging between 1 and 30, with a burn-in period of 1,000,000 followed by 250,000 iterations and using the StrAuto pipeline for parallelization [START_REF] Chhatre | StrAuto: Automation and parallelization of STRUCTURE analysis[END_REF]. The K value corresponding to the "upper most hierarchical level of structure" was determined using the likelihood of observing the data (lnP(D)).

CLUMPP 1.1 [START_REF] Jakobsson | CLUMPP: A cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure[END_REF] and DISTRUCT 1.1 [START_REF] Rosenberg | DISTRUCT: A program for the graphical display of population structure[END_REF] were used for graphical output.

A discriminant analysis of principal components (DAPC; [START_REF] Jombart | Discriminant analysis of principal components: A new method for the analysis of genetically structured populations[END_REF] was conducted in ADEGENET [START_REF] Jombart | Adegenet: A R package for the multivariate analysis of genetic markers[END_REF]. Data were transformed into principal components and discriminant analyses were used to maximize variation among-groups while minimizing within group variation. Population locations were used as a group prior. Based on the a-score method, the number of principal components was set to 80, while we retained two discriminant functions. Hierarchical genetic subdivision (i.e., among lineages, among population and among individuals) was analyzed with an analysis of molecular variance (AMOVA; [START_REF] Excoffier | Analysis of molecular variance inferred from metric distances among DNA haplotypes: Application to human mitochondrial DNA restriction data[END_REF][START_REF] Michalakis | A generic estimation of population subdivision using distances between alleles with special reference for microsatellite loci[END_REF] based on the number of different alleles (infinite allele model) as implemented in GenAlex 6.51 [START_REF] Peakall | GenALEx 6.5: Genetic analysis in excel. Population genetic software for teaching and research-an update[END_REF]) and using 999 permutations. Total and pairwise G ST s and D EST s [START_REF] Jost | GST and its relatives do not measure differentiation[END_REF] among lineages and populations were estimated in GenoDive [START_REF] Meirmans | Genotype and genodive: Two programs for the analysis of genetic diversity of asexual organisms[END_REF]. These two measures provide complementary information about genetic structure. Indeed, the G ST can be considered as an estimate of the nearness to fixation while D EST provides information about the level of allelic differentiation [START_REF] Jost | Differentiation measures for conservation genetics[END_REF]. Significance was tested with a permutation procedure (n = 1000) for the two measures.

The evolutionary relationships among populations were characterized using the chord distance (D c ; [START_REF] Cavalli-Sforza | Phylogenetic analysis. Models and estimation procedures[END_REF].

Allele frequencies were bootstrapped 100 times in SEQBOOT in PHYLIP 3.6 [START_REF] Felsenstein | PHYLIP (phylogeny inference package) version 3.6. Distributed by Author[END_REF]. The consensus tree was obtained using GENDIST, NEIGHBOR and CONSENSUS in PHYLIP and drawn with FigTree v.1.4.4 (https://github.com/ramba ut/figtree). S2). Quality assessment of the four de novo transcriptome assemblies obtained for each population showed that those built for P. cf. grayi (Baleeira, Portugal) and the Italian population of P. clavata (Altare) were superior to those obtained for the other two P. clavata populations (La Vaca and Balun; Table S3). Transrate quality metrics for the two top assemblies were generally similar (average Transrate contig score = 0.36 for P. cf. grayi and 0.33 for P. clavata; Table S3). BUSCO analysis of metazoan single-copy orthologues indicated that the curated assemblies (i.e., following ORF clustering and filtering of non-Anthozoa sequences) from both the Baleeira and Altare populations were very complete, with 98% and 97.4% of core proteins detected, respectively (Table S3).

After curation, a total of 73,871 contigs were assembled in P. cf. grayi (Baleeira population), with an N50 of 1452 bp (Table 1).

| Phylogenetic analysis

Our partitioned ML analysis of 261 single-copy orthologues clearly separated the Atlantic (southern-central Portugal) P. cf. grayi and the Mediterranean P. clavata samples with high (Atlantic clade, 98%) or full (Mediterranean) support (Figure 2). Within the Atlantic samples, the yellow morphotype of P. cf. grayi from three distinct locations (including P39 off Cape Espichel where both color morphs occur in sympatry) formed a monophyletic group to the exclusion of the purple color morph (Figures 1 and2). Within the Mediterranean, individuals from Croatia (Adriatic Sea) also formed a fully-supported clade with full support and distinct from the Spanish and Italian samples (Figure 2).

Despite the high bootstrap support for the Atlantic and

Mediterranean clades, concordance analysis of the branches supporting these clades revealed very low gCFs of only 1.15% (Atlantic) and 2.3% (Mediterranean). Concordance factors for sites (sCF) were also low to moderate (45.3% and 58.3%), which together with the short branch lengths supporting these clades, suggests recent divergence with incomplete lineage sorting (ILS). We note, however, that the branch supporting the yellow morphotype clade within the Atlantic samples was somewhat longer than those supporting the Atlantic and Mediterranean clades, had full bootstrap support, and a fairly robust sCF (62.8%; Figure 2).

| Species delimitation

Species delimitation based on the ASTRAL-III species tree and performed with the SODA software supported the clades formed by samples of P. cf. grayi and P. clavata as separate species (Figure 2). The same split was also evident with the barcode marker mt-mutS for which our analyses recovered two distinct haplotypes, one including all P. clavata samples from the Mediterranean (96% bootstrap support) and another including the samples of P. cf. grayi sequenced here together with P. grayi from the northern Iberian Peninsula (100% bootstrap support; Figure S1). Furthermore, within the Atlantic, the purple and yellow morphotypes of P. cf. grayi were all recovered as two distinct species in the SODA analysis (with the caveat that our dataset only included two purple individuals, although they were sympatric with one individual of the yellow morph -site P39). These results were congruent with the PoMo ML species tree, which supported the Mediterranean (P. clavata) and Atlantic (yellow and purple morphotypes of P. cf. grayi) populations as three distinct entities with full bootstrap support (Figure 3). Given that our sampling included just a few individuals of the purple morphotype, we conservatively assign both color morphs of the Atlantic clade to P. cf. grayi, a species first described (as Acanthogorgia grayi) from Madeira, Portugal [START_REF] Johnson | Description of a second species of Acanthogorgia (J. E. GRAY) from Madeira[END_REF], but which has also been documented to occur in southern Portugal [START_REF] Grasshoff | Die Gorgonarien des östlichen Nordatlantik und des Mittelmeeres III. Die Familie Paramuriceidae (Cnidaria, Anthozoa)[END_REF][START_REF] Grasshoff | Die Flachwasser-Gorgonarien von Europa und Westafrika (Cnidaria, Anthozoa)[END_REF][START_REF] Monteiro-Marques | Report of the hydrographic institute. A plataforma continental do Algarve: definição qualitativa das biocenoses de substrato móvel[END_REF]. This species identification is thus based on molecular data presented here and elsewhere [START_REF] Poliseno | Historical biogeography and mitogenomics of two endemic Mediterranean gorgonians (Holaxonia, Plexauridae)[END_REF], as well as on historical taxonomic records documenting P. grayi (but not other Paramuricea species) at our sampling sites in Portugal.

| Population genetics based on microsatellite loci

| Genetic variability

The observed (H o ) and expected (H e ) heterozygosity ranged between 0.44 (TAV) and 0.83 (PZO; mean ± SD = 0.66 ± 0.08) and between TA B L E 1 Summary statistics for the top transcriptome assembly obtained with rnaSPAdes using samples of Paramuricea cf. grayi from Baleeira, Portugal (NE Atlantic). S1. Please note that there are two terminal nodes for each specimen, one for each of two multi-locus "pseudoalleles" (i.e. alleles derived from in silico phasing of the short read data). The pseudoalleles are distinguished by either a zero or a one appended to the taxon name. Photo credits: A. C. Ferreira (P. cf. grayi) and J. Garrabou (P. clavata).

Metric

0.5 (BAL) and 0.8 (PHL; mean ± SD = 0.71 ± 0.07), respectively.

Multilocus F IS were between 0 (e.g., VAC, LAG, BAL) and 0.3 (TAV; mean ± SD = 0.07 ± 0.08). Significant heterozygote deficiencies, after FDR correction, were found in 15 of 46 populations (28%), located both in the Atlantic (e.g., TAV, SIN) and in the Mediterranean (e.g., CAS, RIO). Regarding rarefied allelic (Ar ( 18) ) and private allelic richness (Ap ( 18) ), we found highly contrasting values among populations. Ar (18) ranged between 3.5 (BAL) and 7.92 (PHL; mean ± SD = 6.51 ± 1.1), while Ap (18) was null for 14 populations (all belonging to the Mediterranean, e.g., CPS, CLB, MES) with a maximum of 1.52 for TAV (mean ± SD = 0.17 ± 0.3). Seven of the highest Ap (18) values were observed for Atlantic populations (LAG, BAL, SEG, SIN, RAB; Table S4).

Significant differences among the three lineages (see below) were observed for all the genetic diversity parameters (H o , H e , Ap (18) , Ar (18) )

with the exception of F IS . Pairwise Wilcoxon rank tests showed that significant differences were mostly driven by variation between Mediterranean and the Atlantic populations, which included the two segregating lineages identified (see Tables S4 andS5 for details). The highest H o , H e and A r values were observed in P. clavata, whereas the yellow lineage of P. cf. grayi harbors the highest A p .

| Population structure

The pattern of population structuring was largely congruent with the phylogenetic analyses. The lnP(D) parameter increased until K = 4, plateaued until K = 9 and decreased for higher K values (Figure S2). The colors correspond to species identified by SODA analysis (plus outgroup, P. biscaya). For the NE Atlantic, the species name is given as P. cf. grayi with color morphotypes (lineages) indicated in yellow and purple. The codes of the sites from which the samples originate are embedded in the taxon names (P39, BAL, TAV, VAC, ALT and BALU) and follow the nomenclature used in Figure 1 and Table S1. Photo credits: A. C. Ferreira (P. cf. grayi) and J. Garrabou (P. clavata). Mediterranean (P. clavata; mean membership coefficient = 0.92; Figure 4a).

The DAPC analysis corroborates the results obtained with STRUCTURE. The first two axes represented relatively high levels of the total variation, 40.4% and 10.7%, respectively. The first axis separated P. clavata populations from the Mediterranean and P. cf. grayi populations from the Atlantic (purple and yellow lineages), while the second axis divided the two lineages of P. cf. grayi (Figure 4b). The topology of the tree based on the Chord distance (D c ) supported the distinction between the Atlantic and the Mediterranean populations with a bootstrap value of 85 (Figure 4c).

The AMOVA revealed highly significant genetic structuring between the three groups (i.e. lineages), among populations, and among individuals (p < .001; Table S5). The genetic variation explained by differences among groups was 14%, higher than that explained by differences among populations (10%) or among individuals (6%). The global G ST and D EST were 0.17 and 0.53, respectively (p < .001). Pairwise G ST among lineages were 0.24, 0.12, and 0.11 for S1. Only bootstrap values higher than 75 are shown.

Atlantic Ocean (0.97) Mediterranean Sea (0.96)

K=2 F A R R A B C O R C A T S IN S E G B A L L A G T A V M IR T Y L M T H C P S IB Z C A B C L B M E S C D V P L L M E D M E T V A C C C S S R E P H L R IO R IS M O R P C L P T C C A S G C L G P R M T M P T V IM P G A B P Z U P L U P Z O P M L G G L A L T S L O F U L A Y V

| DISCUSS ION

This study discovered novel species diversity in the gorgonian genus Paramuricea in the Iberian Peninsula, revealing recently diverged species found predominantly in the Mediterranean (P. clavata) and in the Atlantic (P. cf. grayi), where two color morphotypes (yellow and purple) remain genetically distinct in sympatry, indicating also speciation within P. cf. grayi. These discoveries result from integrating genetic data from three independent types of markers (microsatellites, single-copy orthologues, and a mitochondrial gene) to assess the extent of divergence between Mediterranean populations of the gorgonian P. clavata and adjacent Atlantic populations occurring along southern and western

Portugal, suggested to represent the western limit of the species' distribution range [START_REF] Boavida | Comparison of small remotely operated vehicles and diveroperated video of circalittoral benthos[END_REF]Cúrdia et al., 2012;[START_REF] Pilczynska | Genetic diversity increases with depth in red gorgonian populations of the Mediterranean Sea and the Atlantic Ocean[END_REF]. The three phylotranscriptomic approaches used here (ML phylogenetic inference, SODA species delimitation and PoMo species tree reconstruction) strongly support the Atlantic and Mediterranean clades being two distinct sister species: P. cf. grayi in the Atlantic (shared mt-mutS haplotype with published P. grayi), a species previously described for this region [START_REF] Grasshoff | Die Gorgonarien des östlichen Nordatlantik und des Mittelmeeres III. Die Familie Paramuriceidae (Cnidaria, Anthozoa)[END_REF][START_REF] Grasshoff | Die Flachwasser-Gorgonarien von Europa und Westafrika (Cnidaria, Anthozoa)[END_REF][START_REF] Poliseno | Historical biogeography and mitogenomics of two endemic Mediterranean gorgonians (Holaxonia, Plexauridae)[END_REF]; and P. clavata in the Mediterranean. Our phylogenetic results are also congruent with genetic structure at microsatellite loci (Pilczynska, Cocito, et al., 2017; this study), with a recent phylogenetic reconstruction based on mitochondrial genomes assembled from RNA-seq data [START_REF] Coelho | Complete mitochondrial genome of the branching octocoral Paramuricea grayi (Johnson, 1861), phylogenetic relationships and divergence analysis[END_REF], as well as with a phylogenomic study for P. grayi/P. clavata collected within the Mediterranean [START_REF] Quattrini | Phylogeography of Paramuricea: The role of depth and water mass in the evolution and distribution of deep-sea corals[END_REF]. Although based on a limited number of samples, two statistically consistent species tree methods (SODA, based on the ASTRAL species tree, and PoMo, incorporating withinpopulation genetic diversity) revealed the genetic segregation of color morphotypes within Atlantic P. cf. grayi, and thus an additional putative (or ongoing) speciation event between purple and yellow lineages. Although these lineages display considerable geographic structure, with the yellow morphotype largely dominant in the Algarve (south Portugal) and the purple dominating habitats on the west coast of Portugal, they can also be found in sympatry (e.g., site P39).

| Lineage divergence across the Atlantic-Mediterranean transition

Previous genetic studies of Atlantic P. cf. grayi (heretofore P. clavata)

occurring along Portugal and Mediterranean populations of P. clavata identified high levels of population differentiation between the two regions, indicating highly restricted gene flow across the Atlantic-Mediterranean transition [START_REF] Pilczynska | Genetic diversity increases with depth in red gorgonian populations of the Mediterranean Sea and the Atlantic Ocean[END_REF]Pilczynska, Cocito, et al., 2017). Based on comprehensive sampling, our microsatellite analyses corroborate these early findings, revealing a clear segregation between P. cf. grayi and P. clavata. The uniqueness of the Atlantic populations of P. cf. grayi and the relationship with P. clavata was also evident in the STRUCTURE and DAPC analyses (Figure 4), which showed a two-level hierarchical partitioning of genetic variation that agreed with the results of the phylogenetic analyses based on single-copy ortholog alignments (Figures 2 and3), as well as with the statistically significant differences at most genetic parameters analyzed. Private allelic richness (Ap ( 18) ), in particular, was substantially higher in P. cf. grayi, which included seven of the highest Ap (18) values observed across all populations studied. This pattern is more consistent with the long-term persistence of distinct lineages, rather than Atlantic populations being derived as marginal extensions of a Mediterranean expansion. Indeed, all phylogenetic methods, as well as a distance-based approach, indicate that Mediterranean P. clavata is more genetically divergent from the adjacent Atlantic populations of P. cf. grayi (yellow lineage) than from the more geographically distant populations in western Portugal (purple lineage).

Although our analyses did not consider morphological data, the ML analysis of the mt-mutS barcode showed that the samples of the Atlantic lineage(s) identified here have the same haplotype as the morphospecies P. grayi from northern Iberia (Galicia and Bay of Biscay; [START_REF] Poliseno | Historical biogeography and mitogenomics of two endemic Mediterranean gorgonians (Holaxonia, Plexauridae)[END_REF] and distinct from that of the sister species P. clavata (Figure S1). Thus, while we cannot formally reject the occurrence of P. clavata in adjacent Atlantic coral assemblages (e.g., unsampled populations in northern Morocco), our results confirm the identity of all southern and western Iberian populations studied here as P. cf. grayi. This finding agrees with previous morphological identifications documenting the presence of P. grayi in southern Portugal, particularly in circalittoral rocky habitats along the Algarve and in the Gorringe Bank [START_REF] Grasshoff | Die Gorgonarien des östlichen Nordatlantik und des Mittelmeeres III. Die Familie Paramuriceidae (Cnidaria, Anthozoa)[END_REF][START_REF] Grasshoff | Die Flachwasser-Gorgonarien von Europa und Westafrika (Cnidaria, Anthozoa)[END_REF][START_REF] Monteiro-Marques | Report of the hydrographic institute. A plataforma continental do Algarve: definição qualitativa das biocenoses de substrato móvel[END_REF].

The sister relationship between P. grayi and the Mediterranean endemic P. clavata has been recently described in two studies. Based on mt-mutS sequence data, [START_REF] Poliseno | Historical biogeography and mitogenomics of two endemic Mediterranean gorgonians (Holaxonia, Plexauridae)[END_REF], estimated that the divergence between the two species occurred 2.6 Ma (4-1 Ma) or 4.6 Ma (7-3 Ma) assuming a mutation rate of 0.14% or 0.25% Myr-1, respectively. The authors suggest that speciation occurred through vicariance in association with the Messinian (5.96-5.33 Ma) and/or Gelasian (2.6 Ma) crises [START_REF] Krijgsman | Chronology, causes and progression of the Messinian salinity crisis[END_REF][START_REF] Por | Tethys returns to the Mediterranean: Success and limits of tropical re-colonization[END_REF]. Indeed, the diversification of marine taxa in the Mediterranean is thought to be greatly influenced by these two events, and by the subsequent climatic oscillations during the Quaternary, contributing to the high biodiversity and endemism of present-day Mediterranean fauna and flora [START_REF] Bianchi | Mediterranean Sea biodiversity between the legacy from the past and a future of change[END_REF]Bianchi & Morri, 2000;[START_REF] Lejeusne | Climate change effects on a miniature ocean: The highly diverse, highly impacted Mediterranean Sea[END_REF].

More recently, [START_REF] Quattrini | Phylogeography of Paramuricea: The role of depth and water mass in the evolution and distribution of deep-sea corals[END_REF] 4a). This suggests that the Alboran sea is an area where the two species come into secondary contact, presumably following vicariant speciation. Interestingly, the admixture signal indicates that introgression occurs predominantly from the Atlantic (i.e., P. cf. grayi) into the Mediterranean (i.e., P. clavata) gene pool, and involves alleles found in the purple lineage of P. cf. grayi (Figure 4a). However, this inferred directionality could be a sampling artifact, since Atlantic samples from the adjacent northern Moroccan coastline are currently unavailable. Such a pattern of admixture has been described in other Atlantic-Mediterranean sister taxa/lineages [START_REF] Bermejo | Marine forests of the Mediterranean-Atlantic Cystoseira tamariscifolia complex show a southern Iberian genetic hotspot and no reproductive isolation in parapatry[END_REF][START_REF] Duranton | The spatial scale of dispersal revealed by admixture tracts[END_REF][START_REF] Duranton | The contribution of ancient admixture to reproductive isolation between European sea bass lineages[END_REF][START_REF] Tine | European sea bass genome and its variation provide insights into adaptation to euryhalinity and speciation[END_REF] but see [START_REF] Alberto | Genetic differentiation and secondary contact zone in the seagrass Cymodocea nodosa across the Mediterranean-Atlantic transition region[END_REF], and can be explained by the inflow of Atlantic surface seawater into the Mediterranean resulting from water balance deficits inherent during interglacial periods, which contributes to the continued penetration of Atlantic fauna and flora [START_REF] Bianchi | Mediterranean Sea biodiversity between the legacy from the past and a future of change[END_REF][START_REF] Harmelin | Transfers of bryozoan species between the Atlantic Ocean and the Mediterranean Sea via the strait of Gibraltar[END_REF]. While further research is needed to fully resolve the evolutionary history of P. grayi and P. clavata in the Atlantic-Mediterranean, our findings seem to support suggestions that P. clavata is a "neoendemic" species that diverged from a subtropical Atlantic ancestor that colonized the Mediterranean following the Zanclean flood [START_REF] Bianchi | Mediterranean Sea biodiversity between the legacy from the past and a future of change[END_REF].

The divergence ages estimated recently [START_REF] Poliseno | Historical biogeography and mitogenomics of two endemic Mediterranean gorgonians (Holaxonia, Plexauridae)[END_REF][START_REF] Quattrini | Phylogeography of Paramuricea: The role of depth and water mass in the evolution and distribution of deep-sea corals[END_REF] suggest that the Alboran sea (through the Strait of Gibraltar and AOF) has acted as a buffer reducing admixture between P. cf. grayi and P. clavata throughout the Quaternary, as documented previously for other taxa as a consequence of the coupling between oceanographic barriers (here AOF) and genetic barriers [START_REF] Ayari | The hidden side of a major marine biogeographic boundary: A wide mosaic hybrid zone at the Atlantic-Mediterranean divide reveals the complex interaction between natural and genetic barriers in mussels[END_REF].

The low and moderate values of concordance factors obtained for genes (gCFs) and sites (sCFs) supporting the Atlantic and Mediterranean clades provide further evidence for a recent divergence between P. clavata and P. cf. grayi (Figure 2). Gene-tree discordance is generally attributed to two processes: (1) incomplete lineage sorting (ILS) of ancestral polymorphisms; and (2) gene flow between taxa via introgressive hybridization [START_REF] Maddison | Gene trees in species trees[END_REF][START_REF] Maddison | Inferring phylogeny despite incomplete lineage sorting[END_REF][START_REF] Suh | The dynamics of incomplete lineage sorting across the ancient adaptive radiation of Neoavian birds[END_REF]Yu et al., 2013). Introgression is considered a pivotal process countering lineage diversification of both plants and animals, which is often observed in rapidly (or recently) diverging lineages of closely related species [START_REF] Berner | The genomics of organismal diversification illuminated by adaptive radiations[END_REF][START_REF] Mallet | Hybrid speciation[END_REF][START_REF] Seehausen | Hybridization and adaptive radiation[END_REF], including several coral lineages [START_REF] Combosch | Trans-Pacific RAD-seq population genomics confirms introgressive hybridization in eastern Pacific Pocillopora corals[END_REF][START_REF] Mao | The roles of introgression and climate change in the rise to dominance of Acropora corals[END_REF][START_REF] Poliseno | Hybridization and cryptic speciation in the tropical eastern Pacific octocoral genus Pacifigorgia[END_REF][START_REF] Quattrini | A next generation approach to species delimitation reveals the role of hybridization in a cryptic species complex of corals[END_REF][START_REF] Van Oppen | The evolutionary history of the coral genus Acropora (Scleractinia, cnidaria) based on a mitochondrial and a nuclear marker: Reticulation, incomplete lineage sorting, or morphological convergence?[END_REF][START_REF] Vollmer | Hybridization and the evolution of reef coral diversity[END_REF]. Disentangling ILS and introgressive hybridization is often difficult because they generate similar genomic signatures [START_REF] Meyer | Disentangling incomplete lineage sorting and introgression to refine species-tree estimates for Lake Tanganyika cichlid fishes[END_REF]. However, the congruency between species tree methods (i.e., concatenation and statistically consistent ASTRAL and PoMo models) and between phylogenies based on nuclear and mitochondrial loci analyzed here, as well as evident admixture in the Alboran contact zone suggest that ILS is a more likely explanation for the overall low values of gCFs and sCFs supporting the P. cf. grayi and P. clavata clades [START_REF] Toews | The biogeography of mitochondrial and nuclear discordance in animals[END_REF][START_REF] Twyford | Next-generation hybridization and introgression[END_REF].

| Genetic segregation within the Atlantic lineages

Our results clearly show lineage divergence within Iberian P. cf. grayi.

The microsatellite analysis revealed that populations in southern Portugal (Sagres, Lagos and Tavira) were highly differentiated from those sampled along the western coast (Berlengas, Cape Espichel and Sines), with highly significant pairwise G ST s and D EST s and genetic partitioning (Table S6 and Figure 4). This south-west genetic divide is coincident with the distribution of yellow and purple color morphotypes noted in previous studies examining a subset of the populations [START_REF] Pilczynska | Genetic diversity increases with depth in red gorgonian populations of the Mediterranean Sea and the Atlantic Ocean[END_REF]Pilczynska, Cocito, et al., 2017).

Importantly, our ML and species trees support color morph segregation in P. cf. grayi irrespective of geographic origin, including sympatric individuals (P39 site, Figures 123), which is important preliminary evidence for a diagnostic morphological trait in distinguishing the two Atlantic lineages. It is noteworthy that our STRUCTURE analysis identified a signal of admixture between the two lineages at sites where field surveys confirmed them to occur in sympatry (CAT and SIN; Figure 4a). Collectively, these findings and observations suggest a high degree of habitat segregation (and partial reproductive isolation) between the two morphotypes of P. cf. predominantly philopatric recruitment (e.g., [START_REF] Ayre | Genotypic diversity and gene flow in brooding and spawning corals along the Great Barrier Reef, Australia[END_REF][START_REF] Gazulla | Demo-genetic approach for the conservation and restoration of a habitat-forming octocoral: The case of red coral, Corallium rubrum, in the Réserve naturelle de Scandola[END_REF]Ledoux, Garrabou, et al., 2010;[START_REF] Mokhtar-Jamaï | From global to local genetic 20457758[END_REF][START_REF] Smilansky | Fine-scale genetic structure in the surface brooding Caribbean octocoral, Antillogorgia elisabethae[END_REF][START_REF] Underwood | Multiple scales of genetic connectivity in a brooding coral on isolated reefs following catastrophic bleaching[END_REF] but see discussions in Coelho & Lasker, 2016b;[START_REF] Jones | Larval retention and connectivity among populations of corals and reef fishes: History, advances and challenges[END_REF].

Interestingly, the spatial genetic patterns described for P. clavata in the Mediterranean [START_REF] Arizmendi-Mejía | Combining genetic and demographic data for the conservation of a Mediterranean marine habitat-forming species[END_REF][START_REF] Ledoux | Postglacial range expansion shaped the spatial genetic structure in a marine habitatforming species: Implications for conservation plans in the eastern Adriatic Sea[END_REF][START_REF] Mokhtar-Jamaï | From global to local genetic 20457758[END_REF][START_REF] Thoma | Molecular and morphological diversity of sea fans with emphasis on deep-sea octocorals of the order Alcyonacea Lamouroux[END_REF] coincide with the low effective dispersal inherent to surface brooding inferred from field observations and parental genetic analysis [START_REF] Coma | Reproduction and cycle of gonadal development in the Mediterranean gorgonian Paramuricea clavata[END_REF]Mokhtar-Jamaï et al., 2013; but see [START_REF] Guizien | Survive or swim: Different relationships between migration potential and larval size in three sympatric Mediterranean octocorals[END_REF], while future studies could indicate whether P. cf. grayi is also a brooder. We note, however, that reproductive mode in octocorals is relatively plastic even at the congeneric level (e.g., Antillogorgia: [START_REF] Coelho | Reproductive biology of the Caribbean brooding octocoral Antillogorgia hystrix[END_REF]Coelho & Lasker, 2016a;[START_REF] Gutiérrez-Rodríguez | Reproductive biology, development, and planula behavior in the Caribbean gorgonian Pseudopterogorgia elisabethae[END_REF]and Alcyonium: McFadden et al., 2001;reviewed in Kahng et al., 2011).

| CON CLUDING REMARK S

In this study, we investigated population differentiation and speciation in key structural species of octocorals in marine ecosystems of the NE Atlantic and Mediterranean Sea. The results presented advance the understanding of lineage diversification in these important biogeographic regions. Our analyses revealed that past species identifications based on morphological characters and a limited number of genetic markers have underestimated species diversity in Paramuricea occurring in shallow-water and circalittoral habitats of southern Iberia. This has important implications for conservation planning, as misidentifying putative species as isolated populations compromises the effectiveness of the adopted measures. Specifically, our results provide compelling evidence defining species boundaries between Atlantic populations, here assigned to the morphospecies P. cf. grayi, and the Mediterranean P. clavata;

as well as between two segregating lineages within the P. cf. grayi clade that are coincident with coloration of the colonies, thus casting doubts about the taxonomic status of these two entities. While the conservation status of P. grayi remains unassessed, P. clavata has been listed as vulnerable in the IUCN red list of anthozoans in the Mediterranean since 2016 [START_REF] Otero | Overview of the conservation status of Mediterranean anthozoans[END_REF].

Although the broader species limits between P. cf. grayi and P. clavata could be resolved with a single mitochondrial barcode (mt-mutS), this study emphasizes the need for genome-wide multi-locus datasets to resolve recent lineage splitting events in the face of, for example, ILS and/or introgression. The marked spatial genetic structure and differentiation observed between the P. cf. grayi/P. clavata clades reveal historical footprints of divergence across the Atlantic-Mediterranean transition, the corresponding evolutionary scenarios of which could be tested with further population genomic data. Like recent studies, our results emphasize the urgent need to re-evaluate species diversity and the processes underlying diversification within the genus Paramuricea [START_REF] Coelho | Complete mitochondrial genome of the branching octocoral Paramuricea grayi (Johnson, 1861), phylogenetic relationships and divergence analysis[END_REF][START_REF] Poliseno | Historical biogeography and mitogenomics of two endemic Mediterranean gorgonians (Holaxonia, Plexauridae)[END_REF][START_REF] Quattrini | Phylogeography of Paramuricea: The role of depth and water mass in the evolution and distribution of deep-sea corals[END_REF][START_REF] Thoma | Molecular and morphological diversity of sea fans with emphasis on deep-sea octocorals of the order Alcyonacea Lamouroux[END_REF], which despite its ecological importance, likely harbors high levels of cryptic diversity (e.g., several species described as amphi-Atlantic; [START_REF] Sampaio | Natural history collections as a basis for sound biodiversity assessments: Plexauridae (Octocorallia, Holaxonia) of the Naturalis CANCAP and tyro Mauritania II expeditions[END_REF][START_REF] Walting | Biology of deep-water octocorals[END_REF]. Finally, cryptic diversity that went undetected in previous molecular and morphological studies validates the usefulness of phylotranscriptomics for resolving the evolutionary relationships of octocorals, and more generally of diverse lineages of marine taxa, including at shallow phylogenetic scales as examined here. Although phylotranscriptomics can be a reliable (and in some cases cheaper) alternative to several emerging approaches based on genome sequence data (e.g., [START_REF] Herrera | RAD sequencing enables unprecedented phylogenetic resolution and objective species delimitation in recalcitrant divergent taxa[END_REF][START_REF] Quattrini | Universal target-enrichment baits for anthozoan (Cnidaria) phylogenomics: New approaches to long-standing problems[END_REF]; see discussion in [START_REF] Cheon | Is phylotranscriptomics as reliable as phylogenomics?[END_REF], it remains largely underutilized in anthozoans (reviewed in Quek & Huang, 2022).

AUTH O R CO

  bly was evaluated with TransRate v1.0.1 using both sequence-based and read mapping (onto the assembled contigs) metrics[START_REF] Smith-Unna | TransRate: Reference-free quality assessment of de novo transcriptome assemblies[END_REF]. The read-metrics mode measures the quality of 20457758, 2023, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/ece3.9740 by Cochrane France, Wiley Online Library on [30/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons LicenseF I G U R E 1 Map showing the sampling locations of Paramuricea across the NE Atlantic Ocean and Mediterranean Sea (a). Site codes follow the abbreviations of Table

  For a detailed visualization of the location of Mediterranean sites see figure of Mokhtar-Jamaï et al. (2011). AOF -approximate location of the Almeria-Orán oceanographic front. (b) Detail of sampling sites along southern Iberia (NE Atlantic) and Alboran Sea (Mediterranean) separated by the strait of Gibraltar; insets show sites sampled in cape Espichel (upper) and Sagres (lower) with bathymetric isobaths shown as follows: 50 m, 100-500 m (increments of 100 m) and >500 m (increments of 200 m). (c) Locations sampled in Banyuls-Sur-Mer (southern France) and North-Eastern Spain, including the Medes Islands. (d) Locations sampled along the eastern Adriatic Sea (Croatia). Colors of the markers correspond to species identified by SODA analysis (see results). Note that both the purple and yellow morphotypes of the P. cf. grayi were sampled at P39 (cape Espichel, western Portugal).

  Sampling in the Atlantic took place between 2010 and 2013 in Portugal. The occurrence of P. cf. grayi (heretofore Paramuricea aff. clavata) along the southern coast of Portugal was known from deep divers' records, from shallow hard bottoms down to 90 m depth

  2.6.2 | DNA extraction, amplification and genotypingGenomic DNA of the Atlantic samples was extracted with a cetyltrimethylammonium bromide protocol (CTAB;[START_REF] Winnepenninckx | Extraction of high molecular weight DNA from molluscs[END_REF], with purification by standard chloroform:isoamyl alcohol (24:1) followed by DNA precipitation. All samples (N = 384 colonies) were genotyped at six microsatellite loci: Parcla_09, Parcla_10, Parcla_12, Parcla_14, Parcla_17 (Molecular Ecology Resources Primer Development Consortium, 2010) and Par_d[START_REF] Agell | Isolation and characterization of eight polymorphic microsatellite loci for the Mediterranean gorgonian Paramuricea clavata[END_REF]. The dataset from[START_REF] Mokhtar-Jamaï | From global to local genetic 20457758[END_REF] and the additional Mediterranean sites (N = 1243 samples) were extracted using a salting out procedure[START_REF] Mokhtar-Jamaï | From global to local genetic 20457758[END_REF]. All samples were genotyped as in[START_REF] Mokhtar-Jamaï | From global to local genetic 20457758[END_REF] with minor modifications of the PCR protocol (for more details see Appendix S1).

  20457758, 2023, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/ece3.9740 by Cochrane France, Wiley Online Library on [30/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License 3 | RE SULTS 3.1 | Phylogenomic inference and species delimitation using RNA-seq data 3.1.1 | Transcriptome sequencing and quality of the assemblies RNA-sequencing generated a global dataset of 143 Gb high-quality paired-end reads (28.04 ± 0.91 M reads per sample [mean ± SE]), 38 and 105 Gb for P. cf. grayi and P. clavata, respectively (Table

F

  Phylogenetic relationships in Paramuricea determined by ML partition analysis in IQ-TREE 2. The tree was built using an alignment of 261 single-copy orthologues with a separate evolutionary model per gene (see text for details). The colors correspond to species identified by SODA analysis (plus outgroup, P. biscaya); text boxes to the right indicate geographic region from which samples originate. Red arrows indicate selected key branches with ultrafast bootstrap support (1000 replicates); gCF/sCF = gene/ site concordance factors supporting the branch indicated. For the NE Atlantic, the species name is given as P. cf. grayi with color morphotypes indicated in yellow and purple. The codes of the sites from which the samples originate are embedded in the taxon names (P39, BAL, TAV, VAC, ALT and BALU) and follow the nomenclature used in Figure1and Table

  Considering our focus on the genetic structure among Atlantic and Mediterranean populations, we describe only the results for K = 2 and 3. At K = 4, structuring among Mediterranean populations became apparent, which has been previously characterized (Mokhtar-Jamaï et al., 2011). For K = 2, individuals from the Atlantic (purple and yellow lineages; mean membership coefficient = 0.97) were separated from those from the Mediterranean (mean membership coefficient = 0.96; Figure 4a), with a signal of admixture between the purple lineage of P. cf. grayi and the Mediterranean populations of P. clavata close to the Strait of Gibraltar (MIR, TYL and MTH), and extending to southern Spain (CPS) and the Balearic Islands (CAB). The three genetic clusters (K = 3) were concordant with the areas corresponding to western (P. cf. grayi purple lineage; mean membership coefficient = 0.93) and southern (P. cf. grayi yellow lineage; mean membership coefficient = 0.97) Portugal and the F I G U R E 3 Species tree of Paramuricea samples in this study using the ML polymorphism-aware (PoMo) model and HKY substitution model in IQ-TREE 2. The tree was built using allelic counts from 261 single-copy orthologues. Values on the branches indicate ultrafast bootstrap support (1000 replicates).

F

  Spatial genetic structure in Paramuricea: (a) clustering analysis conducted with STRUCTURE considering K = 2 and 3. Each individual is represented by a vertical line partitioned in K-colored segments, which represent the individual membership fraction in K clusters. Thin black vertical lines delineate the different locations of the samples while dashed and thick lines show the segregating lineages of P. cf. grayi (PPL for 'purple' and YLW for 'yellow'; see results below) and the two oceanographic basins (Atlantic vs. Mediterranean), respectively. Samples names are shown below the assignment plots (for abbreviations see Table S1). The mean membership coefficient for each cluster is also shown. (b) Scatter plot of the discriminant analysis of principal components (DAPC) based on an a-score of 80 and considering axes 1 and 2. These two axes represented 51.1% of the total variation in the data. Each dot corresponds to one individual (n = 1223) from each of the 46 locations. Note that the colors correspond to the three clusters with the red dots for the Mediterranean samples of P. clavata and the yellow and purple dots for the Atlantic samples of P. cf. grayi. Only the names of the Atlantic samples are shown. (c) Dendrogram based on the chord distance (Cavalli-Sforza & Edwards, 1967) showing the relationships among the 46 populations. Colors are the same as for K = 3 while the sample names are shown in Table

  Conceptualization (equal); data curation (equal); formal analysis (equal); investigation (equal); methodology (equal); resources (equal); validation (equal); visualization (equal); writingoriginal draft (lead); writing -review and editing (equal). Gareth A. Pearson: Conceptualization (equal); data curation (equal); formal analysis (equal); investigation (equal); methodology (equal); validation (equal); visualization (equal); writing -original draft (equal); writingreview and editing (equal). Joana R. H. Boavida: Conceptualization (equal); data curation (supporting); formal analysis (supporting); funding acquisition (equal); investigation (equal); methodology (supporting); resources (equal); validation (supporting); writing -review and editing (equal). Diogo Paulo: Resources (supporting); writingreview and editing (supporting). Didier Aurelle: Formal analysis (supporting); resources (supporting); writing -review and editing (equal). Sophie Arnaud-Haond: Writing -review and editing (equal). Daniel Gómez-Gras: Resources (supporting); writing -review and editing (equal). Nathaniel Bensoussan: Resources (supporting); writing -review and editing (supporting). Paula López-Sendino: Investigation (equal); resources (supporting); writing -review and editing (equal). Carlo Cerrano: Resources (supporting); writing -review and editing (equal). Silvija Kipson: Resources (supporting); writing -review and editing (equal). Tatjana Bakran-Petricioli: Resources (supporting); writing -review and editing (equal). Eliana Ferretti: Resources (supporting); writing -review and editing (supporting). Cristina Linares: Resources (supporting); writing -review and editing (equal). Joaquim Garrabou: Funding acquisition (equal); resources (equal); writing -review and editing (equal). Ester A. Serrao: Conceptualization (equal); funding acquisition (equal); writing -review and editing (equal). Jean-Baptiste Ledoux: Conceptualization (equal); data curation (equal); formal analysis (equal); investigation (equal); methodology (equal); resources (equal); validation (equal); visualization (equal); writing -original draft (equal); writing -review and editing (equal). ACK N O WLE D G E M ENTS This work was funded by program BIOMARES, a Pew Marine Fellowship, and FCT -Foundation for Science and Technology UIDB/04326/2020, UIDP/04326/2020, LA/P/0101/2020, PTDC/ BIA-CBI/6515/2020; EU-BiodivERsA BiodivRestore-253 (FCT: DivRestore/0013/2020); EU Horizon 2020 research and innovation program (MERCES: Grant no. 689518); as well as the Inaqua

  20457758, 2023, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/ece3.9740 by Cochrane France, Wiley Online Library on [30/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

  20457758, 2023, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/ece3.9740 by Cochrane France, Wiley Online Library on [30/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

	COELHO et al.	| 7 of 21

  Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/ece3.9740 by Cochrane France, Wiley Online Library on [30/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

	COELHO et al.	P. biscaya	20457758, 2023, 1, | 11 of 21
			P. cf. grayi_P39
			33
		100	P. cf. grayi_BAL
		100	P. cf. grayi_TAV
		P. cf. grayi_P39
			P. clavata_BALU
		100
		100	P. clavata_VAC
		P. clavata_ALT
	0.2

  P. clavata vs. yellow lineage of P. cf. grayi, yellow vs. purple P. cf. grayi lineages and P. clavata vs. purple lineage of P. cf. grayi, respectively.Overall, pairwise G ST ranged between 0.01 (PPL vs. MED) and 0.46 (AYV and BAL), with only 31 nonsignificant comparisons (<3% of total; TableS6). Pairwise D EST ranged between 0.01 (PPL vs. MED) and 0.99 (AYV and SEG). All except 33 pairwise comparisons (<5%) were significant. Focusing on the Atlantic populations, the values of the pairwise G ST and D EST were within the range of the values

	observed among Mediterranean populations (0.04 [SEG vs. BAL] to
	0.31 [FAR vs. BAL] for G ST ; 0.06 [SEG vs. BAL] to 0.67 [RAB vs. BAL]
	for D EST ). For the two estimators, all but one pairwise comparison
	(SEG vs. BAL) were significant.

  20457758, 2023, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/ece3.9740 by Cochrane France, Wiley Online Library on [30/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

	COELHO et al.	| 13 of 21

  used RAD-seq data to estimate the timing of divergence between P. clavata and P. grayi

	at 2 Ma (3.3-1.0 Ma 95% CI), with Bayesian dispersal-extirpation-
	cladogenesis modeling suggesting that divergence likely occurred in
	sympatry (though the analysis only included three individuals col-
	lected within the Mediterranean). Our STRUCTURE analysis based
	on microsatellite loci identified a signal of admixture between P. cf.
	grayi and P. clavata at sites in the Alboran sea (MIR, TYL, MTH, lo-
	cated immediately east of the Strait of Gibraltar). Here and in two
	additional sites east of the AOF (CPS and CAB), the coefficients of
	membership to either of the Atlantic/Mediterranean clusters were
	generally lower (Figure

  , 2023, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/ece3.9740 by Cochrane France, Wiley Online Library on [30/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

	grayi, calling for
	additional analyses using genomic data, broader taxon and habi-
	tat (e.g., depth) sampling, as well as experimental approaches (e.g.,
	crossing experiments) to clarify the taxonomic status of the two
	color morphs.
	Within segregating lineages of P. cf. grayi, we have identified
	significant pairwise genetic differentiation between most popula-
	tions studied, including at local spatial scales, from ~20 km (e.g., BAL
	vs. LAG, southern Portugal) to less than 400 m (e.g., COR vs. CAT,
	western Portugal; Table S6). These results suggest restricted gene
	flow between populations of both Atlantic lineages, with rare long-
	distance dispersal. The extent to which such genetic differentiation
	is driven by seascape (i.e., physical environment), life-history traits

of the species or the interplay between both remains speculative at this point. However, these findings raise several important questions about the biology of P. cf. grayi. For instance, strong genetic differentiation at small spatial scales and a pattern of isolation by distance are common in many brooding corals in which dispersal potential is generally low given the short duration of the pelagic larval stage and 20457758

  20457758, 2023, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/ece3.9740 by Cochrane France, Wiley Online Library on [30/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

20457758, 2023, 1, Downloaded from https://onlinelibrary.wiley.com/doi/10.1002/ece3.9740 by Cochrane France, Wiley Online Library on [30/01/2023]. See the Terms and Conditions (https://onlinelibrary.wiley.com/terms-and-conditions) on Wiley Online Library for rules of use; OA articles are governed by the applicable Creative Commons License

0.01 P. cf. grayi_BAL_6_1 P. cf. grayi_BAL_3_0 P. clavata_ALT_8_0 P. clavata_VAC_6_0 P. biscaya SRX4389730_0 P. clavata_BALU_5_0 P. clavata_ALT_8_1 P. clavata_ALT_7_1 P. clavata_BALU_7_0 P. clavata_VAC_3_0 P. biscaya SRX4389729_1 P. clavata_BALU_4_0 P. clavata_VAC_8_1 P. cf. grayi_BAL_4_1 P. cf. grayi_P39_46_0 P. clavata_BALU_6_1 P. cf. grayi_P39_57_0 P. biscaya SRX4389730_1 P. clavata_VAC_5_1 P. clavata_VAC_6_1 P. clavata_VAC_9_1 P. clavata_VAC_1_0 P. clavata_ALT_3_1 P. clavata_ALT_4_0 P. cf. grayi_P39_57_1 P. clavata_ALT_1_1 P. cf. grayi_BAL_8_1 P. biscaya SRX4389728_0 P. cf. grayi_P39_54_1 P. clavata_VAC_7_1 P. clavata_VAC_3_1 P. clavata_VAC_8_0 P. clavata_ALT_4_1 P. cf. grayi_BAL_4_0 P. cf. grayi_BAL_6_0 P. clavata_VAC_4_0 P. clavata_ALT_3_0 P. cf. grayi_P39_46_1 P. cf. grayi_TAV_97_1 P. biscaya SRX4389728_1 P. clavata_BALU_6_0 P. cf. grayi_TAV_97_0 P. clavata_VAC_4_1 P. cf. grayi_BAL_3_1 P. clavata_BALU_5_1 P. biscaya SRX4389729_0 P. clavata_VAC_1_1 P. clavata_VAC_7_0 P. clavata_BALU_3_1 P. cf. grayi_P39_54_0 P. clavata_BALU_3_0 P. clavata_ALT_6_1 P. clavata_ALT_6_0 P. clavata_VAC_5_0 P. clavata_BALU_7_1 P. clavata_VAC_9_0 P. clavata_BALU_1_1 P. cf. grayi_BAL_1_0 P. clavata_ALT_1_0 P. biscaya SRX4389727_1 P. clavata_ALT_7_0 P. clavata_BALU_9_0 P. cf. grayi_BAL_5_1 P. biscaya SRX4389727_0 P. cf. grayi_BAL_1_1 P. clavata_BALU_1_0 P. cf. grayi_BAL_5_1 P. clavata_BALU_9_1 P. clavata_BALU_4_1 P. cf. grayi_BAL_8_0

DATA AVA I L A B I L I T Y S TAT E M E N T

Raw RNA-seq data and mt-mutS sequences have been uploaded to NCBI's SRA (BioProject ID: PRJNA847883) and GenBank databases (Accessions numbers: ON804207-ON804214), respectively. Final sequence alignments and microsatellite data are available on DRYAD (https://doi.org/10.5061/dryad.r7sqv 9sfm).