Sixty-nine samples were collected from multiple sources and locations throughout the North Atlantic and Mediterranean Sea at depths of 30 to 2,312 m ( Figure 1 and Tables 1 , S1 ). Most specimens used for sequencing were preserved onboard the ship in 95% EtOH and stored either at room temperature or frozen at -20°C ( Table S1 ). However, the samples collected off Canada were frozen on board the ship at -20°C and then transferred to 100% EtOH and stored at room temperature ( Table S1 ).

Small sections of the axis (1-2 polyps in length) were treated with an 8% sodium hypochlorite solution to dissolve organic matter, followed by one rinse with 3% hydrogen peroxide, one rinse of deionized water, and three rinses of 95% ethanol. Sclerites were then sorted using an ultrafine bristle paintbrush and assessed using a JEOL 5600LV Scanning Electron Microscope and a compound microscope ( Supplemental File 1 ). Sclerite sizes were measured with JEOL visualization software. Specimens have been given tentative species names or designations based on morphological characters (Deichmann, 1936; Grasshoff, 1977; Table 2 ) combined with RAD-Seq results; however, we acknowledge that all await further taxonomic confirmation as a comparison to type material is essential. Therefore, we retained the use of cf. (to confer) or aff. (affinity to) for some of the morphospecies designations ( Table 2 ).

DNA was extracted using a Qiagen DNeasy kit following the manufacturer’s protocol. DNA concentration was assessed using a Qubit v 2.0 fluorometer with a Broad Range Assay Kit. Quality was assessed by running 100 ng of DNA for each sample on a 1% agarose gel and checked with a NanoDrop spectrophotometer. Normalized (20 ng per ul), high-quality (230/260 and 260/280 ratios >1.8) DNA was sent to Floragenex Inc (Eugene, OR) for RAD-Seq library preparation. Libraries were constructed for each of the 69 samples using the 6-cutter PstI enzyme. Due to some degradation of DNA, two libraries were constructed for a subset of the samples (n=47). Each library was sequenced (100 bp SE reads) on one and a quarter lanes of an Illumina HiSeq2500 (University of Oregon’s Genomics and Cell Characterization Core Facility lab).

In addition to RAD-Seq, two gene regions (mtMutS, 28S rDNA) used widely for barcoding in octocorals were PCR-amplified and Sanger-sequenced using published primers and protocols (McFadden et al., 2011; McFadden et al., 2014). mtMutS and 28S rDNA sequences of Paramuricea were aligned separately using the E-INS-i method in MAFFT (100PAM/k=2, gap open penalty 1.53, Katoh and Toh, 2008) and then concatenated. Additional, outgroup sequences from the octocoral families Acanthogorgiidae (Acanthogorgia), Gorgoniidae (Pacifigorgia), and other plexaurid genera (Thesea, Villogorgia, Muriceides, Bebryce, Plexaura) were downloaded from GenBank and included in the analyses to determine the relationship of sample LII-501 (an unidentified species in the family Plexauridae) to the other Paramuricea samples. Plexaura kuna was used as the outgroup, as this genus has been previously shown to fall outside of the other included taxa in larger datasets (Wirshing et al., 2005; McFadden et al., 2006).

MOTHUR v1.44 (Schloss et al., 2009) was used to delimit molecular operational taxonomic units (MOTUs) based on genetic distance thresholds applied to each gene separately and to the concatenated dataset of mtMutS+28S (e.g., McFadden et al., 2014; Quattrini et al., 2019). The distseq method was used to calculate uncorrected pairwise distances, counting a string of gaps as one (calc=onegap) without penalizing gaps at the ends (countends=F). Phylogenetic trees were constructed separately for mtMutS, 28S rDNA, and the concatenated dataset using maximum likelihood (IQTree; Nguyen et al., 2015) ( Supplemental Figure S1 ). We used the built-in program ModelFinder (-m MFP+Merge, Kalyaanamoorthy et al., 2017) to identify the best partitioning scheme and models of evolution (Bayesian Information Criterion, BIC) to use in these analyses (mtMutS: K2P+R2; 28S rDNA: TPM2u+F+G4). Both ultrafast bootstrapping (-bb 1000, Hoang et al., 2018) and the Sh-like approximate likelihood ratio test (-alrt 1000, Guindon et al., 2010) were performed.

Sequence reads were de-multiplexed and trimmed for low quality using the process_radtags program in Stacks v 1.35 (Catchen et al., 2013). Illumina adapters and barcodes were excluded from reads, and lengths were truncated (-t) to 91bp. The rescue barcodes and RAD-tags options were enabled (-r). Reads with ambiguous bases (-c) and low-quality scores (-q, <10) were discarded. The default sliding window (-w 0.15) and average score (-s 10) for filtering were used. ipyrad v0.9.6 (Eaton and Overcast, 2020) was then used for additional filtering and clustering of sequences to identify homologous loci. ipyrad enables the identification of homologous loci across highly divergent samples by including indel variation in the clustering process (Eaton, 2014). Parameters were left at default values, based on recommendations in the ipyrad manual (https://ipyrad.readthedocs.io/). They included the maximum number of low-quality bases (phred quality score <20) per read of five, removal of low-quality reads (with more than five of bases with a phred quality score <20), short (< 35 bp) reads, and reads with more than five Ns, in addition to a strict filter for adapters. To detect potential paralogs, we changed the maximum number of heterozygous sites retained in a locus to 25% across samples to reduce the chances of clustering paralogs with fixed differences. We also used a clustering threshold of 85% based on previous RAD-Seq studies of octocorals (Pante et al., 2015; Herrera and Shank, 2016; Quattrini et al., 2019). Different datasets ( Table 3 ) were generated in step seven of the ipyrad pipeline for analyses outlined below.

A phylogeny of all taxa (n=69 individuals) was constructed using IQTree with a concatenated data matrix of no more than 25% missing taxa per locus ( Table 3 ). The best fit model (-m MFP, Kalyaanamoorthy et al., 2017) was chosen (TVM+F+R5) using BIC followed by maximum likelihood analyses with ultrafast bootstrapping (-bb 1000, Hoang et al., 2018) and the Sh-like approximate likelihood ratio test (-alrt 1000, Guindon et al., 2010). RAxML v 8 (Stamatakis, 2006) was also used to construct a phylogeny using a GTR+G model and rapid bootstrapping; a total of 20 ML searches and 200 bootstrap replicates were performed. Trees were visualized with FigTree and rooted to the unidentified Plexauridae (sample LII-501) as this was an appropriate outgroup supported by our mtMutS and 28S results ( Supplemental Figure S1 ).

Structure v2.3 (Pritchard et al., 2000) was used to define putative species of Paramuricea to help guide the coalescent-based SNAPP analyses (below). We used an unlinked SNP dataset (u.str file from ipyrad) with a maximum of 25% missing taxa per locus ( Table 3 ) for the deep-sea (> 200 m depth) clade. Structure was run in parallel with StrAutoParallel v 1.0 (Chhatre and Emerson, 2017) using an admixture model with correlated allele frequencies. Burnin was set to 25,000 followed by 1,000,000 MCMC generations. The inferred number of populations (K) was set from 1 to 12, and 10 runs of each K were conducted. Because the commonly used ΔK method often finds only the uppermost genetic structure (Evanno et al., 2005; Janes et al., 2017), which is a likely scenario when examining several metapopulations, we chose three K values to plot based on the three increases apparent in ΔK values ( Supplemental File 2 ) from Structure Harvester (Earl and Vonholdt 2012). The ten runs for each of those K values (K=4, 6, and 10) were aligned with CLUMPP v 1.2 (Jakobsson and Rosenberg, 2007), and the resulting indivq file was input into Distruct v. 1.1 (Rosenberg, 2004) for the graphical display of individual population assignments.

Coalescent-based SNAPP v 1.3 (Bryant et al., 2012) analysis was used to infer the species tree using unlinked SNPs ( Table 3 ). Prior to SNAPP analysis, the full dataset was subsampled in ipyrad (step 7) to ≤ 3 individuals per population per putative species with no more than 10% missing taxa per locus (SNAPP_subset1 dataset). Using the unlinked_SNPs_to_nexus.py script (https://github.com/brunoasm/usnps_to_nexus written by B. deMedeiros), unlinked SNPs (usnps file from ipyrad) were converted into a binary nexus file. BEAUti v2.6 (Bouckaert et al., 2019) was used to create an xml file. Mutation rates (u=3.333, v=0.588) were set based on Quattrini et al. (2019), and the coalescent rate parameter was estimated during the MCMC. The Yule tree prior was set to a gamma distribution with default alpha (2.0) and beta (200) parameters, and the ancestral population sizes were set at the default gamma distribution priors. One run (MCMC=2M generations, samplefreq=1000) was performed in BEAST2 v2.6 (Bouckaert et al., 2019) using CIPRES (Miller et al., 2010). Convergence (effective sample size, ESS, values > 200) and mixing were checked in Tracer v. 1.7 (Rambaut et al., 2018). Tree topologies were drawn in DensiTree v2.2 (Bouckaert, 2010). A maximum clade credibility tree and posterior probabilities were generated with Tree Annotator v 2.3 (burnin=25%, mean heights, Bouckaert et al., 2014). We also repeated the SNAPP analysis with another subsampling of the dataset (SNAPP_subset2) to examine potential bias in subsampling.

Divergence dating estimation was conducted in SNAPP v 1.3 following the tutorial (https://github.com/ForBioPhylogenomics/tutorials/tree/main/divergence_time_estimation_with_snp_data) and based on the methods of Stange et al. (2018) to determine the timing of diversification of Paramuricea into deep water. The ruby script snapp_prep.rb was used to generate an xml file based on unlinked SNPs. The SNAPP maximum clade credibility tree was used as a topology constraint. Two relative dates were used based on 95% highest posterior density (HPD) node ages from Poliseno et al. (2017). We approximated these 95% HPD intervals for the ingroup and the clade containing P. clavata and P. grayi using lognormal distributions in real space. Lognormal distributions were centered around 7.7 and 4.2 MYA with standard deviations of 0.28 and 0.4 MYA for the ingroup (4.2-12.8 MYA 95% CI) and the clavata/grayi (1.8-8.5 MYA 95% CI) clades, respectively. A strict clock was also used, mutation rates were set as noted above, and the effective population sizes were unlinked. One run (MCMC=1M generations, samplefreq=1000) was performed in BEAST2 v2.6 (Bouckaert et al., 2019) using CIPRES (Miller et al., 2010). Convergence and ESS values > 200 were confirmed in Tracer v. 1.7. A maximum clade credibility tree with 95% HPD ages was generated with Tree Annotator v 2.3 (burnin=25%, mean heights).

To examine the ancestral depth ranges of Paramuricea, we used the Bayesian dispersal- extirpation-cladogenesis (DEC) model (Ree and Smith, 2008) implemented in RevBayes (Höhna et al., 2016). Ancestral ranges of shallow-mesophotic (10-200 m) and deep-sea (> 200 m) depths were estimated following the guidelines in the online tutorial on simple analysis of historical biogeography (https://revbayes.github.io/tutorials/biogeo/biogeo_simple.html). The ultrametric SNAPP tree constructed without the clock model, with the topology that matched the concatenated analyses (undated-SNAPP_subset1), was input into RevBayes. To estimate depth ranges for each putative species/population, we used the depths from sample collections supplemented by depth records in the literature (Deichmann, 1936; Grasshoff, 1977) in cases where we were confident about species identities. We also estimated the ancestral ranges of Paramuricea in particular water masses to determine whether water masses facilitated the dispersal of deep-sea species across the North Atlantic Ocean. Water masses that are most often found at particular depths in each of the regions (following Rivas et al., 2005; Townsend et al., 2015; Wilson et al., 2015; Georgian et al., 2016; Aldama-Campino and Döös, 2020; Liu and Tanhua, 2021) were used in analyses: Antarctic Intermediate Water (AAIW), Labrador Sea Water (LSW), North Atlantic Central Water (NACW), Mediterranean Overflow Water (MOW), and North Atlantic Deep Water (NADW) (see Figure 1 ). In RevBayes, the number of generations (MCMC) was set to 5000. The R package RevGadgets (Tribble et al., 2022) was then used to generate phylogenies with the 1st through 3rd most probable ancestral ranges shown (as pie charts) before and after cladogenic events.

We also modeled continuous trait evolution of temperature and depth ranges to further elucidate the ancestral states of the environment in the deep-sea clade. We used the contMap function in the R package phytools (Revell, 2012). This mapping approach uses a fast estimation of the most likely ancestral states at each node. For this analysis, we used the IQTree tree generated with concatenated data for all 69 individuals. Depth of collection was used as the depth data. Temperature data were taken from modal temperatures of water masses (following Rivas et al., 2005; Townsend et al., 2015; Wilson et al., 2015; Aldama-Campino and Döös, 2020; Liu and Tanhua, 2021) where those water masses occur in each region. We chose to use these water mass temperatures instead of temperatures taken at the time of collection because we expect that long-term temperatures, as reflected in water masses, are more important in the evolution of this genus than a “snapshot” that is captured during a single collection time.

Diagnostic sclerites were apparent for each shallow/mesophotic species and for the monophyletic groups/genetic clusters of individuals from deep waters as identified in RAD-Seq analysis (see below). Hereafter, we refer to these species as: P. clavata, P. grayi, P. macrospina, P. sp. 1- 4, P. cf. placomus, P. cf. biscaya, and P. aff. biscaya. All putative species exhibited differences in the morphology of the calicular thornscales and/or sclerites in the coenenchyme of the axis or calyx ( Table 2 , Supplemental File 1 ). P. aff. biscaya was similar to P. cf. biscaya, but P. aff. biscaya differed by the presence of lacinated plates with rounded projections in the coenenchyme. P. sp. 1 and 2 were also morphologically similar but differed in that P. sp. 2 had triangular projections in thornscales. P. sp. 3 contained unique branched spindles in the coenenchyme, which were not observed in the other clades. Two morphological variants for P. cf. placomus were observed, with differences in the size of the thornscale projection. However, genetics did not distinguish these as separate lineages (see below).

We obtained mtMutS and 28S rDNA sequences for 59 and 56 Paramuricea individuals, respectively. The mtMutS alignment was 735 bp, 28S was 652 bp, and the concatenated alignment was 1,387 bp. Although a monophyletic deep-sea (> 200 m) clade was recovered in all phylogenies, the topologies were largely incongruent with respect to the placement of the mesophotic species P. macrospina and several putative species within the deep-sea clade ( Supplemental Figure 1 ). In addition, most nodes were poorly supported (< 90% b.s.) in each phylogeny, except the node for the entire Paramuricea clade.

Based on a 0.3% genetic distance threshold, the mtMutS barcode identified seven molecular operational taxonomic units (MOTUs): one MOTU each for the Plexauridae outgroup, P. clavata. P. grayi, P. macrospina; and three MOTUs in the deep-sea clade, but none was congruent with the RAD-Seq phylogeny or Structure analyses ( Figure 2 ). For 28S, seven MOTUs were designated at a 0.3% genetic distance threshold: one MOTU each for Plexauridae, P. macrospina, P. clavata/grayi; and four MOTUs in the deep-sea clade. The MOTUs based on 28S for the deep-sea clade corresponded largely with monophyletic clades of the RAD-Seq phylogeny but did not correspond with groupings of the Structure analysis ( Figure 2 ). Based on a 0.3% genetic threshold, the concatenated mtMutS+28S barcode designated nine MOTUs: one MOTU each for Plexauridae, P. macrospina, P. clavata, and P. grayi; and five MOTUs in the deep-sea clade, none of which were congruent with the RAD-Seq phylogeny or Structure analyses ( Figure 2 ).

After trimming and filtering in both stacks and ipyrad, 307M reads were retained for 69 Paramuricea individuals. On average, 13% of reads were removed due to unmatched RAD tags, low quality and unambiguous bases, and short reads. The mean number of reads retained per sample was 4,453,204 ± 3,589,920 SD. The dataset, allowing for 25% missing data per locus, included 2,266 ± 600 SD loci per sample. In total, 29,202 SNPs were obtained, and the alignment matrix consisted of 249,438 bp.

The phylogeny based on a concatenated RAD-Seq data alignment (249,438 bp) was well supported at most nodes, based on rapid and ultrafast bootstrapping and the Sh-like approximate likelihood ratio test implemented in IQTree ( Figure 2 ). The RAxML and IQTree topologies were largely congruent. The only exceptions were the relationships among populations of P. cf. placomus. RAxML recovered a sister relationship between the Ireland and GoM samples and IQtree recovered a sister relationship between GoM and NE United States samples.

The undated-SNAPP_subset1 species tree ( Figures 3A, B ) inferred from 958 unlinked SNPs ( Table 3 ) was strongly supported (> 95%) at only 50% of the nodes. Still, the recovered relationships largely matched those in the concatenated tree analyses. The only exception was that P. placomus clade was sister to P. sp. 1 and 2 in the SNAPP analysis with 100% support ( Figures 3A, B ). Relationships recovered in the undated-SNAPP_subset2 species tree, inferred from 802 SNPs and a different subset of individuals ( Supplemental Figure S3 ), were incongruent with other topologies. Inferences from this dataset placed the unidentified plexaurid (LII-501) as sister to P. clavata/grayi (90% pp support) and P. aff. biscaya sister to P. sp. 3 (70% pp support, Supplemental Figure S3 ). The dated-SNAPP_subset1 and SNAPP_subset2 species trees also produced somewhat dissimilar topologies ( Figure 3C and Supplemental Figure 3 ) from the undated analyses. In the dated-SNAPP_subset 1 tree, P. sp. 4 was sister to P. sp. 3 (68% pp support, Figure 3C ) whereas the dated-SNAPP_subset2 tree matched the undated-SNAPP_subset2 tree ( Supplemental Figure 3 ).

A monophyletic deep-sea (> 200 m) clade was recovered with strong support in all tree topologies inferred from RAD-Seq data ( Figures 2 , 3 and Supplemental Figure S1 ). The deep-sea samples formed seven well-supported monophyletic clades within the deep-sea clade, four of which contained individuals from multiple locations (i.e., Ireland and the GoM). However, these monophyletic clades were not recovered in the mtMutS+28S phylogeny ( Supplemental Figure S1 ).

Structure analyses were conducted on the deep-sea clade (n=63 individuals) using 3,821 unlinked SNPs ( Table 3 and Figure 2 ). The ΔK method of Evanno indicated that the most likely number of K groups was 4 ( Supplemental Figure S2 ), corresponding to: 1) P. sp. 1 and 2, 2) P. cf. placomus, 3) P. sp. 3, and 4) a cluster of P. cf. biscaya, P. aff. biscaya, and P. sp. 4. However, because ΔK often finds only the uppermost structure in a hierarchical system, we plotted the membership probability for K=6 and 10 as indicated by the ΔK plot ( Supplemental Figure S2 ). Both analyses further separated P. biscaya, P. aff. biscaya, and P. sp. 4 into three distinct clusters. K=6 indicated that P. sp. 1 and sp. 2 are grouped, whereas the K=10 and K=4 analysis indicated some differentiation between these two morphospecies. Thus, we chose the cluster designations as indicated in the K=10 Structure plot for subsequent analyses. These cluster designations also matched morphospecies designations.

Divergence dating of the SNAPP_subset1 and SNAPP_subset 2 species trees indicated that Paramuricea arose approximately 8 MYA (5-13 MYA 95% CI). Around 4.5 MYA (1.9-7.0 MYA 95% CI), P. macrospina and the deep-sea clade diverged from one another ( Figure 3C ). At about 2 MYA (1.0-3.3 MYA 95% CI), both Paramuricea clavata and P. grayi diverged from one another. The deep-sea clade diversified into the deep North Atlantic around 1 to 2 MYA (0.5-3.0 MYA 95% CI).

Several putative species within the deep-sea clade, supported by tree inferences from RAD-Seq data and the SNAPP species trees ( Figures 2 , 3 and Supplemental Figure 3 ), comprised specimens from different geographic locations across the North Atlantic. These results indicate broad distributions (from Ireland to the Gulf of Mexico) for at least three species (P. cf. biscaya, P. cf. placomus, P. sp 3). Furthermore, it appears that species diverged first across depth and/or temperature gradients ( Figure 4 ), and then spread into different geographical regions multiple times. A strong pattern of diversification associated with depth and temperature was apparent across the phylogeny. This pattern suggests a general diversification trend from onshore to offshore and warmer to colder temperatures ( Figure 4 ). Many species were restricted to relatively narrow modal temperature and/or depth ranges. The exceptions were P. sp. 3 and P. cf. placomus, found across wide depth ranges (800-1600 m and 230-1300 m, respectively) yet narrow temperature ranges (4-6°C and 8-10°C) ( Figure 4 ).

The DEC model indicated that a likely vicariance event approx. 4.5 MYA drove the divergence between the mesophotic Mediterranean lineage of P. macrospina and the deep-sea clade ( Figure 5A and Supplemental Figure 4 ). In addition, this divergence was likely from a common ancestor with a broad depth range and an estimated ancestral, optimal temperature state of about 10°C. In contrast, the divergence between P. clavata and P. grayi was likely in sympatry, with P. clavata diversifying in mesophotic depths of the Mediterranean Sea from an ancestor with a broad depth range. P. grayi retained its ancestral depth range ( Figure 5A ).

We further examined the evolution and distribution of the deep-sea clade with respect to contemporary water mass distributions using a DEC model and the SNAPP_subset1 trees. Two clades, which were reciprocally monophyletic in the SNAPP_subset 1 trees, occurred in either present-day upper (NACW, MOW) and intermediate water masses (LSW) or intermediate (AAIW, LSW) and lower water masses (NADW) ( Figure 5B and Supplemental Figure S5 ). Five out of seven putative deep-sea species were found to occupy only one water mass, whereas P. cf. placomus and P. sp. 3 were found in two and three water masses, respectively. The DEC model suggested that P. cf. biscaya and P. aff. biscaya diverged in sympatry in NADW. The DEC model also indicated that P. sp. 4 and P. cf. biscaya/aff. biscaya diverged from a common ancestor with an ancestral range of LSW and NADW. This divergence was likely due to a vicariance event around 750K to 1 MYA. Also, at about 1 MYA, P. sp. 1 and 2 diverged from a common ancestor with an ancestral range consisting of LSW and NACW.

RAD-Seq was highly effective for reconstructing evolutionary relationships among Paramuricea. Most nodes were well supported, and reciprocally monophyletic clades that generally corresponded with morphospecies were recovered in the phylogenetic analyses. Furthermore, structure analysis indicated several clusters with little to no admixture between them. In general, each SNAPP species tree showed slight incongruence among relationships, all with lower node support, compared to the concatenated phylogenies. SNAPP can be sensitive to the number of samples and SNPs and the proportion of missing data. We set the maximum proportion of missing data in our datasets to 10%. This resulted in datasets with relatively few SNPs (< 1000), which likely explains the observed incongruences (Schmidt-Lebuhn et al., 2017; O’Connell and Smith, 2018). Regardless, most of our analyses converged on the same evolutionary relationships ( Figures 2 and 3A, B , Supplemental Figure S3 ). Of note, many of the samples used in the analysis contained some amount of high molecular weight DNA, but most of them also had a high amount of degraded DNA. Yet, we still recovered 100s to 1000s of unlinked SNPs per sample. Thus, our results and others (Herrera and Shank, 2016; Quattrini et al., 2019) indicate the utility of RAD-Seq in octocorals that have some level of DNA degradation.

The monophyletic clades in the phylogenetic analysis and cluster designations from Structure analyses helped guide morphospecies designations and vice versa; and thus, conducting these analyses in concert resulted in higher confidence of putative species designations. In all but one case, morphological evaluation agreed with genetic clade designations. Two morphological variants for P. cf. placomus were found, but structure analysis across different values of K indicated that these are likely one species. In another case, P. sp. 1 and 2 were quite similar morphologically ( Supplemental File 1 ), and there was some degree of admixture between them. Still, they formed well-supported monophyletic clades that diverged into different environments around 1-2 MYA. We remain conservative in assigning names to many of these putative species because comparison to and incorporation of type material in a systematic framework is essential. Our analyses included at least nine out of 19 valid species of Paramuricea, but we acknowledge that our current study is only a first step towards a robust systematic framework, which is essential to illuminate the systematics of this genus.

The need for robust species delimitation analyses is exemplified by a recent population genomic study of P. biscaya in the Gulf of Mexico, which utilized some of the same samples herein (Galaska et al., 2021). Based on >12K SNPs, Galaska et al. (2021) found two lineages of P. biscaya, lineage 1 and lineage 2, represented in our study as P. cf. biscaya and P. aff. biscaya, respectively. Lineage 1 was generally found deeper (mean depth 1,896 m) than lineage 2 (mean depth 1,454 m); however, these lineages co-occur at some intermediate-depth sites and 15% of assessed specimens had an admixed ancestry (see Figure 1 in Galaska et al., 2021). Interestingly, the morphology of the specimens examined in this study differs between these two lineages, with P. aff. biscaya containing plates with projecting, rounded spines in the calyx and P. cf. biscaya containing plates without the protruding spines. These findings raise several questions. First, are these lineages presently diverging in sympatry in the Gulf of Mexico, as the DEC model suggests? Or, did they diverge in allopatry and are now undergoing hybridization in secondary contact? On the other hand, it is curious that we found little to no admixture between these lineages in the present study, begging whether this is a result of SNP calling bias. In our study, SNPs called in the deep-sea clade were from a range of species that diverged over 1 to 2 MY. Therefore, we likely included SNPs associated with conserved regions. In contrast, Galaska et al. (2021) had three times as many SNPs from five times more specimens in just two sister lineages. Further investigation is needed to distinguish the degree of admixture from incomplete lineage sorting and elucidate the divergence process between these two groups.

Similar to other phylogenetic studies of octocorals (Pante et al., 2015; Herrera and Shank, 2016; Quattrini et al., 2019), our results suggest that the conventionally-used, single-locus barcodes are insufficient to reconstruct evolutionary relationships or delimit species within speciose genera and rapid and/or recent radiations of octocorals. Except for the shallow/mesophotic species, neither mtMutS or 28S rDNA alone nor in combination was sufficient to recover the well-supported, reciprocally-monophyletic clades inferred from RAD-Seq data. There was also no congruence between clusters denoted by Structure and MOTUs in the deep-sea group. In addition, morphospecies contained multiple MOTUs, yet the same MOTU could be found in several putative species. Radice et al. (2016) noted this pattern in mtMutS haplotypes in Paramuricea and concluded there were potentially several species off Canada due to several different mtMutS haplotypes. However, the data presented here support high mtMutS haplotype diversity in one species (P. sp. 4) off Canada. We also documented high haplotype diversity within species collected off Ireland. In contrast, some of those identical haplotypes found in P. sp. 4 or species off Ireland (e.g., type B, type A) were each found in a single species in the Gulf of Mexico. Similar levels of incongruence between mitochondrial COI and RAD-Seq results were also found in the two P. biscaya lineages (Galaska et al., 2021). This incongruence between RAD-Seq and single-locus barcodes is likely caused by incomplete lineage sorting of mitochondrial and nuclear rDNA genes. In conclusion, we suggest that neither mtMutS nor 28S should be used to guide species designations for deep-sea species of Paramuricea as it could lead to inaccurate and potentially confused species identities (also see Thoma, 2013), particularly without accompanying morphological data.

Our results provide a first survey of the phylogeography of a genus of octocoral common to the North Atlantic Ocean and Mediterranean Sea. Diversification of Paramuricea occurred throughout the Pliocene and Pleistocene geological epochs. Notably, the deep-sea clade of Paramuricea and the lineage of the mesophotic Mediterranean species P. macrospina diverged around 4.5 MYA from an ancestor common in both deep and mesophotic depths, as indicated by the DEC model. The DEC model indicated that this divergence occurred via a vicariance event, similar to the findings based on ancestral range models and mitochondrial genomes of Poliseno et al. (2017). In addition, the timing of this event is within the range inferred (5.9-3.3 MYA) by Poliseno et al. (2017), although the phylogenetic relationships differed between studies. A monophyletic deepwater clade was not inferred by Poliseno et al. (2017). Rather, it was split into two clades, with P. macrospina sister to one. Poliseno et al. (2017) postulated that the Messinian salinity crisis, which caused sea-level change and desiccation through evaporation in the Mediterranean Sea, was a likely driver of this divergence event. At the end of the Messinian salinity crisis, approximately 5.3 MYA, the Mediterranean Sea re-filled as the connection to the North Atlantic reopened. Shortly after that, around 4.5 MYA, the restriction of the Central American Seaway reached a critical threshold (Karas et al., 2017). Both events caused drastic changes in ocean circulation and climate (Roveri et al., 2014; Karas et al., 2017), including changes in water mass transport, and ventilation and formation of NADW in the North Atlantic (Steph et al., 2010). Together, these changes in ocean circulation might have caused population isolation and divergence of Paramuricea species across depth.

Two to three million years later, the deep-sea clade of Paramuricea diversified throughout the deep North Atlantic Ocean. Repeated glaciations marked this time in the Pleistocene (2.8 MYA to 11 KYA). During the Mid-Pleistocene Transition, around 1 MYA, a significant shift in the behavior of glacial cycles occurred. This transition period was marked by changes in deep water formation and ocean circulation (Pena and Goldstein, 2014; Poirier and Billups, 2014). In addition, there was increased warming at depth and outflow of Mediterranean Outflow Water from the Mediterranean Sea (Catunda et al., 2021). Changes in deep ocean circulation over the past 1 million years could have sparked the spread and diversification of Paramuricea throughout the deep North Atlantic. These results are consistent with the inferred timing of colonization and spread of Paragorgia arborea in the North Atlantic (Herrera et al., 2012), and thus may be indicative of a critical time for the establishment of modern deep-sea coral diversity.

The importance of the environment in the evolution of Paramuricea is clear from the DEC model and ancestral state reconstructions of depth and temperature. These analyses do not unequivocally support a shallow (< 200 m) origin for Paramuricea, as indicated by the estimated ancestral depth range from DEC analysis ( Figure 5 ). However, the phytools analysis ( Figure 4 ) suggested a general trend of diversification into deeper and colder depths, particularly in the deep-sea clade. This indicates a diversification into the deeper bathyal zone and lends support for the onshore-offshore hypothesis (Jablonski et al., 1983). Together, the high lineage diversity found in the bathyal region and the inferred divergence across depths and water masses followed by diversification into different regions, also constitute support for the depth-differentiation hypothesis (Rex and Etter, 2010).

In some cases, the DEC model indicated that divergence was due to vicariance events, so subsequent adaptations to environmental variables must have occurred following speciation. For example, the diversification of the deep-sea clade of Paramuricea in deep and relatively cold water must have subsequently led to biological adaptations that precluded invasion back into mesophotic depths. Changes in hydrostatic pressure and just a few degrees of temperature can elicit biochemical changes (Somero, 1992). In other instances, the DEC model suggested divergence in sympatry (i.e., P. aff. biscaya and P. biscaya), which could be driven by adaptive divergence to environmental variables that co-vary with depth. The latter is difficult to tease apart from secondary contact following allopatry without additional analyses. Nevertheless, species appear to have adaptations to living within relatively narrow temperature and/or depth ranges.

Although ocean circulation changed during the glacial and interglacial periods of the Pleistocene (e.g., Demenocal et al., 1992; Bell et al., 2015; Kim et al., 2021), the contemporary distributions of Paramuricea spp. in water mass type still show an evolutionary signal. In the SNAPP_subset1 species tree, the two reciprocally monophyletic clades of deep-sea Paramuricea were found in either upper/intermediate or deep/intermediate water masses. The common ancestor of Paramuricea sp. 4 and P. cf. biscaya/aff. biscaya had a likely distribution range in LSW-NADW. The DEC model suggested that a vicariance event led to the divergence of each species in either LSW or NADW masses following the Mid-Pleistocene Transition ( Figures 5B , S4 ). Similarly, the ancestral range of the common ancestor of P. sp. 1 and 2 spanned both LSW and NACW, and diverged in allopatry into either LSW or NACW ( Figures 5B , S4 ). Larval dispersal models of deep-sea corals that assume larvae behave as neutrally buoyant particles have indicated that larvae stay close to the depths at which they were spawned and rarely cross water mass density boundaries (Bracco et al., 2019; Liu et al., 2021). Without significant mixing between water masses, water mass boundaries could promote population isolation via physical barriers marked by density differences. When larvae successfully cross into a different water mass, larvae most likely have adaptations that allow their survival under those environmental conditions. Water masses could also connect distant populations as the species appear to occupy large geographic regions within the same water mass. These results collectively highlight the importance of water masses in shaping the evolution of corals throughout the deep North Atlantic Ocean.

We suggest a few considerations for future work on the systematics and biogeography of Paramuricea and octocorals in general. First, the systematics of this genus needs to be reevaluated in a robust species delimitation framework that includes both genomic and morphological data. In particular, type material needs to be imaged and incorporated into genomic studies to help accurately identify species. In this regard, target-capture of conserved elements (Quattrini et al., 2018; Erickson et al., 2021) is a promising method to obtain genomic-scale data from historical museum samples (Untiedt et al., 2021). Putting accurate binomens on individuals is critical for effective restoration and conservation, particularly for comprehensively documenting species’ distributions and population sources and sinks, and thus understanding the degree of population vulnerability to anthropogenic change.

Our analyses used the modal temperatures of the water masses from sites in which species were collected. We recognize there can be more temperature variability with short-term changes in ocean conditions through, for example, eddy formation or seasonal mixing of water masses at particular depths or locations (e.g., Houpert et al., 2016). However, by using modal temperatures of water masses, we could depict the temperature signal over a longer time scale. We expected to see a pattern across the phylogeny if these longer-term temperatures impacted the evolution of the genus. In addition, we gathered depth data from specimen location instead of from the literature or digital databases. This was partly because we were not confident of the accuracy of species identifications in those sources, and we are not yet confident in our determination of species without access to type material. We expect additional patterns to be revealed as the ability to discern Paramuricea species and our knowledge of their distribution improve.