Evolution of mitochondrial and nuclear genomes in Pennatulacea

We examine the phylogeny of sea pens using sequences of whole mitochondrial genomes and the nuclear ri- bosomal cluster generated through low coverage Illumina sequencing. Taxon sampling includes 30 species in 19 genera representing 13 families. Ancestral state reconstruction shows that most sea pen mitochondrial genomes have the ancestral gene order, and that Pennatulacea with diverse gene orders are found in a single clade. The monophyly of Pennatulidae and Protoptilidae are rejected by both the mitochondrial and nuclear dataset, while the mitochondrial dataset further rejects monophyly of Virgulariidae, and the nuclear dataset rejects monophyly of Kophobelemnidae. We show discordance between nuclear ribosomal gene cluster phylogenies and whole mitochondrial genome phylogenies and highlight key Pennatulacea taxa that could be included in cnidarian genome-wide studies to better resolve the sea pen tree of life. We further illustrate how well frequently sequenced markers capture the overall diversity of the mitochondrial genome and the nuclear ribosomal genes in sea pens.


Introduction
The subclass Octocorallia represents cnidarians that have polyps with eight tentacles and eight mesenteries. They can inhabit different regions, depths and oceans (Daly et al., 2007). They comprise the orders Alcyonacea (soft corals and sea fans), Pennatulacea (sea pens) and Helioporacea (blue corals) with approximately 3500 extant species.
Pennatulaceans are intriguing and fascinating animals. They were first reported in the early part of the seventeenth century, and explorers puzzled over whether they might be plants with worms in their roots, that could transform into stones, and could not determine what type of creature they were (Kerr, 1824, in Williams, 2011. How sea pens evolved and how sea pen species are related are still far from being understood. The order Pennatulacea currently has 235 extant accepted species, placed in 15 families and 40 genera (Pérez et al., 2016, García-Cárdenas et al., 2019, López-González and Drewery, 2022, Worms Editorial Board (2022. Through evolutionary history, Pennatulacea have evolved a range of morphologies (with and without axis, with and without sclerites, clubshaped, flattened leaf-shape) and have developed from three to five differentiated and specialized polyp types (Williams et al., 2012). The primary polyp (oozooid) is modified into the rachis (the distal region where the secondary polyps rise) and peduncle (the proximal region that anchors the colony in mostly soft substrates [Williams and Alderslade, 2011]). Their differentiated polyps are the main morphological feature that characterize the Order Pennatulacea (synapomorphy) and distinguish them from all the other octocorals (Williams et al., 2012).
Pennatulaceans have successfully adapted to mainly soft substrate habitats (Williams and Alderslade, 2011), and to different bathymetric ranges globally (from shallow coastal waters down to 6000 m) (Williams, 2011). The first attempts to investigate the phylogeny of sea pens were made by Kölliker (1870) who suggested that species that have sessile polyps were the most primitive (e.g. Umbellula and Protoptilum), and those that have polyps leaves and ridges were derived. Kükenthal (1915) proposed two suborders: Subselliflorae (polyps arranged in polyp leaves) and Sessiliflorae (polyps raised directly from the rachis). He also suggested that the earliest-diverging taxa were those with radial symmetry, e.g., family Veretillidae, believing symmetry evolved from radial to bilateral. Williams (1997) wrote an assessment of sea pen phylogeny, with a detailed review, and pointed out the inadequacy of the suborder classification, and highlighted the need for a systematic review.
One of the earliest molecular studies to incorporate multiple sea pens, found Pennatulacea to be polyphyletic (Berntson et al., 2001), with Umbellula falling distant from the main grouping of sea pens. This study, based on 18S rRNA, contrasted with that of McFadden et al. (2006), based on mitochondrial ND2 and mutS, which recovered Pennatulacea as monophyletic. Both phylogenies were part of higher-level systematic studies and included only a limited number of sea pen representatives. The first molecular phylogenetic study to focus on sea pens combined mitochondrial MutS and NADH dehydrogenase subunit 2 (ND2) sequences (Dolan et al., 2013). The study included 11 families and 14 genera of mostly deep-sea species. This study rejected the monophyly of both suborders and five families (Umbellulidae, Pennatulidae, Kophobelemnidae, Pteroeididae, Protoptilidae). Subsequent studies, including shallower species from the northwestern Pacific Ocean found the families Virgulariidae and Scleroptilidae also nonmonophyletic (Kushida and Reimer, 2019) and that Veretillidae and Echinoptilidae were not ancestral (Kushida and Reimer, 2019). The morphology of polyp leaves, which formed the basis of the suborder classification, was determined to result from convergent evolution (Dolan et al., 2013, Kushida andReimer, 2019). García-Cárdenas et al. (2020), combining the two mitochondrial genes with nuclear 28S rRNA, found similar problems with the existing classification, as did López-González and Drewery (2022) when combining three mitochondrial genes with nuclear 28S rRNA. These studies revealed some consistent clades, but none provided robust support for deeper nodes. Furthermore, the inclusion of 28S rRNA into concatenated datasets appeared to generate instability in the tree, with ML and BI topologies differing in the placement of certain taxa.
Complete mitochondrial genomes can potentially provide a phylogenetic solution since they can easily be assembled from next generation sequencing at low cost (Lavrov, 2007;Kayal et al., 2012), and their utility has been proven at different taxonomic levels (Kayal and Lavrov, 2008;Brockman and McFadden, 2012;Figueroa and Baco, 2013;Kayal et al., 2013;Lavrov and Pett, 2016), although phylogenomics has shown some conclusions inferred from mitogenomics to be misleading (e.g., Medina et al., 2006;Kayal et al., 2013). In 2019, 18 pennatulacean mitochondrial genomes (Hogan et al., 2019) were described from shallow Illumina sequencing. This increased the number of sea pen mitochondrial genomes available for phylogenetic analysis from two to twenty, and several others have become available since then (e.g., Muthye et al., 2022). Herein, the objective is to use those mitochondrial genome sequences in tandem with ribosomal nuclear data assembled from the original Illumina output and from Sequence Read Archives available on GenBank in an attempt to improve understanding of pennatulacean evolution, including their relationships to other octocorals. The specific aims of this study were to use phylogenetic methods to: 1) re-assess the monophyly and position of Pennatulacea within Octocorallia, and identify the most suitable outgroup for further analyses; 2) investigate patterns of gene-order distribution and diversity in the octocorals including sea pens; 3) use whole mitochondrial genomes and nuclear ribosomal cluster sequences to investigate the phylogenetic relationships within Pennatulacea; 4) to investigate which genes are most suitable for progressing molecular systematics in sea pens.

Sampling, sequencing and data assembly of Pennatulacea
Eighteen Pennatulacea morphospecies were previously collected from the Northeast Atlantic, two from shallow water and sixteen from deep water ( Fig. 1; Hogan et al., 2019). Specimens were preserved in 96 % ethanol or frozen either at − 20 • C or − 80 • C. Tissue samples for genetics were taken from polyps and/or peduncle tissue and preserved in 96 % ethanol and stored at − 20 • C. Vouchers are deposited at National Museums Scotland, Edinburgh as NMS.Z.2019.25.1-NMS.Z.2019.25.18, and mitochondrial genomes were previously assembled following shallow Illumina sequencing and accessioned to Genbank as MK919655-MK919673 (Hogan et al., 2019). Herein, nuclear rDNA was assembled using the de novo method in MITObim (Hahn et al., 2013) which relies on the MIRA assembler, using partial 28S or 5.8S gene reference sequences from Pennatulacea species as initial seeds (Supplementary  Table S1). Genes were annotated in Geneious using octocoral sequences from GenBank as references (Supplementary Table S1). Reads were mapped to the assembled sequences using the Geneious Read Mapper (Kearse et al., 2012) to check coverage.
The dataset nuclear ribosomal cluster (nrDNA) comprised complete 18S, 5.8S and 28S genes, and the internal transcribed spacers ITS1 and ITS2 and was accessioned to GenBank under OL741691-OL741709).
For consistency, all octocoral mitochondrial genome sequences from GenBank were reannotated following the same method. Protein coding genes were first identified through Mitos2 (Bernt et al., 2013), with NCBI RefSeq 81 Metazoa database reference and genetic code 4, which includes Coelenterata. Annotations were edited in Geneious prime (Kearse et al., 2012) with respect to the widest open reading frame (ORF) using cnidarian start (ATG) and stop (TAG, TAA) codons.
Protein-coding genes were aligned separately using MAFFT G-INS-I and translation align settings for protein-coding genes in Geneious Prime (Kearse et al., 2012). All alignments were checked manually. Regions that could not be confidently aligned were removed using GBlocks software (Talavera and Castresana, 2007) applying the less stringent selection, which allows smaller final blocks, gap positions within the final blocks, and less strict flanking positions, and with protein coding-genes treated as codons.
The genes were concatenated into a matrix using Geneious Prime (Kearse et al., 2012). For alignment lengths see Supplementary Table S3. DAMBE (Xia and Xie, 2001) was used to test nucleotide base saturation.
Maximum Likelihood (ML) reconstruction was performed using IQ-TREE 2.0 software (Nguyen et al., 2015) with 1000 bootstraps. Model-Finder (Kalyaanamoorthy et al., 2017) implemented in the IQ-TREE software, computed the best partition and best-fit model with Free-Rate model (+R) and Bayesian information criterion (BIC). ModelFinder starts its partition scheme tests with the full partition model (that can be pre-determined and specified into genes regions, codons, type of genes and DNA, or other appropriate subsets), and subsequently merges these subsets to find the optimal partition scheme. Each codon position of each protein coding gene was considered as a partition under the starting full partition model. Bayesian inference (BI) was conducted in MrBayes v3.2.6 (Ronquist et al., 2012), with a default general time reversible model with settings and model parameters as follows: lset nst = 6 rates = invgamma, with two Markov Chain Monte Carlo (MCMC) runs for 10,000,000 generations, tree sampling every 1000 generations. Output was checked for stationarity in Tracer v.1.7.1 and a burn-in of 25 % applied. Rooting Octocorallia trees is problematic because the closest outgroup taxa lack the MutS gene which contains considerable phylogenetic signal. Following Figuero and Baco (2014), we rooted on Briareum asbestinum as this taxon is considered the earliest diverging within Octocorallia (McFadden et al., 2006, Brockman and McFadden, 2012, Park et al., 2012. Phylogenetic trees were visualized and edited using FigTree v. 1.4.4 (https://tree.bio.ed.ac.uk/software/figtree/). The ancestral gene order was reconstructed using ace in the R package ape (Paradis and Schliep, 2019), and using a randomization test (1000 reps) to test for phylogenetic conservancy. Gene orders were mapped onto phylogenetic trees to illustrate their distribution.

Pennatulacea
Phylogenetic analyses were conducted on 30 putative morphospecies of sea pens (Table 1) Analyses used two separate datasets. The first comprised the mitochondrial genome, with 14 genes as described above for octocorals plus two ribosomal rRNA genes (16S rRNA and 12S rRNA). 12S and 16S rRNA were aligned on the MAFFT server (https://mafft.cbrc.jp/align ment/server/) using Q-INS-i which considers the secondary structure of RNA. This dataset did not include Calibelemnon or Malacobelemnon as we could not recover the mitogenomes from the available SRAs. The second dataset comprised the nuclear ribosomal cluster (rnDNA) of three genes. The nuclear dataset lacked Stylatula since the mitochondrial genome for this species was from GenBank and there was no SRA from which to generate nuclear data. Phylogenetic analyses were performed as described above for octocorals, except we used 1,000,000 generations for Bayesian Inference, as analysis in Tracer v.1.7.1 indicated that this was sufficient.

Mitochondrial dataset
We explored two versions of the data: (1) protein coding genes + mt rDNA; and 2) protein coding genes protein coding genes translated to amino acids. Full details of alignment lengths, final partitions and models are provided in Supplementary Table S3. We used the newly generated octocoral phylogenetic tree, based on Table 1 GenBank Accession numbers for Pennatulacea nuclear ribosomal genes newly assembled from listed Sequence Read Archives.

SRA Accession No
R.I. Hogan et al. whole mitochondrial genomes, to determine the most suitable out-group for the mitochondrial dataset (see results).

Nuclear datasets
We again explored different versions of the nuclear dataset: (1) The nuclear ribosomal cluster comprising 5.8S, 28S, 18S (2) the nuclear gene 28S, and (3) the nuclear gene 28S rRNA trimmed to the region limited by the primers 28S-Far (5 ′ CACGAGACCGATAGCGAA CAAGTA 3 ′ ) and 28S-Rar (5 ′ TCATTTCGACCC TAAGACCTC 3 ′ ) (McFadden and van Ofwegen, 2012) using an extended dataset (Table 2) of 73 sequences comprising 23 genera in 15 families, including all available sea pen sequences from GenBank for this primer region. Full details of alignment lengths, final partitions and models are provided in Supplementary  Table S3.

Sliding window analysis
To investigate the variability of these large nuclear and mitochondrial alignments, we conducted a sliding window analysis in the R package Spider: Species Identity and Evolution in R (Brown et al., 2012), implementing SlideAnalyses with window width 99, and interval width 1 for the nuclear ribosomal data and interval width 3 (i.e., codon) for the mitochondrial genomes, outputting mean K2P distance per window.
The monophyly of Pennatulacea has maximum support in the ML and BI analyses ( Fig. 2) with Junceella fragilis (Family Ellisellidae) as a sister taxon, and the calcaxonians from the families Keratoisididae (Acanella eburnea, Keratoisidinae sp.) and Primnoidae (Narella hawaiiensis) as a sister group (Fig. 3), together with the blue coral Heliopora coerulea (Order Helioporacea). Subsequent mitochondrial analyses are rooted on Junceella fragilis.

Gene order
Ancestral State Reconstruction across the octocoral phylogenetic tree determined gene order A as the ancestral gene order (p = 0.003) and that this character shows significant phylogenetic conservancy (p≪0.001). We mapped the octocoral mitochondrial gene orders and ancestral state reconstruction onto our octocoral phylogenetic tree. Variations in gene order from the ancestral gene order A only occur in the Pennatulacea/Calcaxonia and Scleraxonia-Alcyoniina clades.
The Pennatulacea-Calcaxonia clade has at least four different types of gene arrangements, which include the octocoral gene orders type A, B, G and a bipartite genome from Umbellula sp. 3. Two of these gene orders are found in both pennatulaceans and calcaxonians (type A and B). The Scleraxonia-Alcyoniina clade includes gene orders A, C, D and E.

Mitochondrial phylogeny Pennatulacea
Whole mitochondrial genome (including protein coding and mt rRNA genes) ML and BI phylogenies for Pennatulacea recovered identical topologies of four major clades, each of which received maximum support (Fig. 3). Analysing translated protein coding genes resulted in a  Fig. S2), and the same four highly supported major clades.

Fig. 2.
Octocoral phylogeny based on 73 available (newly sequenced from this study and GenBank samples) pennatulacean mitochondrial genomes based on a Bayesian inference (BI) and maximum likelihood (ML) analysis. All nodes marked * have bootstrap values / Bayesian posterior probabilities of 100/1. Support only shown for deep nodes due to space constraints. Ancestral state reconstructions of gene orders mapped to nodes. Three major clades indicated with pale grey boxes; monophyletic Pennatulacea indicated with darker grey box.
Relationships among these three subclades are less well supported than the subclades themselves, but the internal topology is the same in the amino-acid tree ( Supplementary Fig. S2). In Clade 2, Distichoptilum gracile (Family Protoptilidae) is sister taxon to the two Ptilella species (Family Pennatulidae) (BS 98, PP 1). Renilla muelleri (Family Renillidae) is sister to two other Pennatulidae (two Pennatula species) (BS 100, PP 1), and Actinoptilum molle (family Echinoptilidae) is sister to these (BS 100, PP1). These subclades (Distichoptilum + Ptilella and Actinoptilum,(Pennatula + Renilla)) are sister with maximum support (BS 100, PP 1). These relationships are mirrored in the amino-acid tree ( Figure S2) but with slightly less support. Protoptilum carpenteri and Stylatula elongata are the other members of Clade 4, but their position relative to one another and the clade containing the other six species is not clear, does not receive strong support, and varies between trees ( Fig. 3; Supplementary Fig. S2). Regardless, we recover neither Protoptilidae nor Pennatulidae as monophyletic (note also that Pteroeides is in Clade 1), and the recovery of Stylatula in Clade 2 and Virgularia in Clade 1 renders the family Virgulariidae polyphyletic (Fig. 3).
Clade 3 recovered Funiculina quadrangularis as sister to four morphospecies of Kophobelemnidae (Fig. 3). All nodes within the clade received maximum support (BS 100, PP 1). The internal topology of Clade 3 was identical in other trees, but the node uniting Kophobelemnidae received slightly lower support in analyses of the amino acid dataset (Supplementary Fig. S2).

Nuclear rDNA phylogeny Pennatulacea
Trees built with three nuclear ribosomal genes (28S, 18S, 5.8S) vary in topology from our mitochondrial trees (Fig. 4). The 3-gene nuclear trees do not include Stylatula but have Calibelemnon and Malacobelemnon as additional taxa, but they also do not recover the same four major clades. There is some congruence between nuclear trees and the mitochondrial tree as follows: (i) they recover mitochondrial Clade 4 as a monophyletic grouping, with Balticina cf. finmarchica sister to Pseudumbellula with maximum support (BS 100, PP 1), (ii) they recover mitochondrial Clade 2 as a monophyletic grouping with high to maximum support (BS 97, PP 1), and (iii) within Clade 2, there are two subclades, one comprising Distichoptilum gracile sister to two Ptilella species (BS 99, PP 1), the other encompassing Actinoptilum, Pennatula and Renilla (BS 99, PP 0.99).
The 3-gene nuclear trees recovered two other clades. The first, which is well supported (BS 74, PP 0.99), encompasses all the mitochondrial Clade 1 taxa, but also includes Funiculina quadrangularis, which was sister to Kophopbelemnon in Clade 3 in the mitochondrial trees, and Malacobelemnon daytoni (family Kophobelemnidae, and not included in the mitochondrial trees). The second, comprises only Kophobelemnon species. Deeper relationships were supported with the Kophobelemon clade sister to Clade 4 (BS 100, PP 1), and the two other clades sister to each other (BS 70, PP 1).
Trees built on 28S rRNA only were congruent with trees built with three nuclear ribosomal genes, but less well resolved. Clades 2 and 4 were supported (Clade 2 BS 85, PP 0.99; Clade 4 BS 100, PP 1); Kophobelemnon species formed their own clade (BS 100, PP 1) which was sister to Clade 4 (BS 100, PP 1); but the taxa from mitochondrial clades 1 and 3 that had united into a clade in the nuclear 3-gene tree formed multiple unresolved lineages ( Supplementary Fig. S3).
Trees built on a short region of 28S rRNA defined by the primers 28S-Rar and 28S-Far that had better taxon sampling were better resolved than the complete 28S tree but less well resolved than the nuclear 3-gene tree (Fig. 5). Here, mitochondrial Clade 2 was recovered as monophyletic (BS 91, PP 0.98) and included the genera Alloptilella, Gilibelemnon, and Scytalium, sequences of which were not available for our other analyses. Kophobelemnon was recovered in a supported sistertaxon relationship with Gyrophyllum (BS 99, PP 1) and this clade was sister (BS 100, PP 1) to a maximally supported Clade 4 (BS 100, PP1). A clade comprising most of the mitochondrial Clade 1 taxa plus Malacobelemnon daytoni was also recovered (BS 70, PP 0.91) but the position of some taxa (species of Pteroeides, Funiculina and Virgularia mirabilis) recovered in this clade in the nuclear 3-gene analyses was not resolved by this short region of 28S rRNA.

Sliding window analysis
Sliding window analyses revealed the most variable regions of the mitochondrial genome to be the 3 ′ end of the cytochrome b gene and the mid section of MutS, whereas the most variable region of the three nuclear ribosomal genes is the 5 ′ portion of 28S (Figue 6).

Discussion
We generated the first Pennatulacea phylogenetic trees with two extended datasets for sixteen mitochondrial genes (30 taxa; 17,891 characters) and three nuclear ribosomal regions (30 taxa; max 5,363 characters), allowing us to test the congruence of sea pen nuclear and mitochondrial genomes, and examine the utility of widely-used partial gene regions in sea pen systematics. We also generated an octocoral phylogeny using an additional 53 mitochondrial genomes available in GenBank (83 taxa; 14,712 characters), and mapped the evolution of mitochondrial genome rearrangement onto it.

Mitochondrial gene orders
Pennatulacea was recovered as a monophyletic group in the octocoral phylogeny closely associated with Calcaxonia. A Pennatulacea-Calcaxonia clade was first recovered using ND2 and mtMutS (McFadden et al., 2006), and subsequently confirmed in a mitochondrial genome study (Figueroa and Baco, 2014), and most recently using Ultra Conserved Elements (UCEs) (Quattrini et al., 2020, McFadden et al., 2021. These studies found Ellisellidae to be the sister taxon of Pennatulacea (as did we), which agrees with early opinions of Kölliker (1870), as well as Bayer (1955), who noted the similarity of calcification pattern in these groups. Williams (2019) proposed that Pennatulacea could be incorporated as a family of Calcaxonia, or that a new taxon Actinaxonia could combine Ellisellidae and Pennatulacea as sister taxa. We recovered Helioporacea within the clade populated otherwise by Pennatulacea and Calcaxonia, as found in recent studies using UCEs (Quattrini et al., 2020, McFadden et al., 2021 suggesting that mitochondrial genome evolution mirrors the nuclear genome at this taxonomic level. In fact, our mitochondrial genome phylogeny was extremely congruent with that from UCEs also placing plexaurids, gorgoniids, nephtheids, and the alcyoniid taxon Sinularia as a separate clade, and grouping paragorgiids and the alcyoniids Anthomastus and Paraminabea with the Pennatulacea-Calcaxonia-Helioporacea clade. Given that UCEs sample Where accession numbers not shown on tree, see Table 1. broadly across the nuclear genome, that phylogeny is assumed to be correct. Hogan et al. (2019) established that gene order B of Brugler and France (2008) evolved independently in sea pens and bamboo corals and uncovered new mitochondrial gene orders among sea pens. The ancestral state reconstruction herein confirms the two independent evolution events for gene order B and demonstrates the prevalence of novel rearrangement among deep-sea taxa, particularly in Clade 1 of the sea pen tree. All the variation in sea pen gene arrangements known to date occurs in Clade 1 (Fig. 3), which contains at least four different types (gene order A, B, G and bipartite genome). Mitogenomes from further species and genera are required to see whether this feature is really with bootstrap support (BS) and posterior probabilities (PP) from separate Bayesian inference analysis on nodes. * indicates maximum support of 100/1. The position of taxa in relation to the mitochondrial tree is indicated through highlighting blocks of taxa and indicating the mitochondrial clade to which they belong in large text (right). Where accession numbers not shown on tree, see Table 1. unique to Clade 1, or whether it can be found in other parts of the sea pen tree. In the octocoral tree, mitochondrial genome rearrangement occurs among the deep-sea bamboo corals (family Keratoisididae), and among mostly deep-water taxa in the Scleraxonia-Alcyoniina clade (Fig. 2). It is not clear whether genome rearrangement is more prevalent in deep-water taxa, or simply that more mitogenomes are available for deep-water species. It could simply be that gene order is more conserved in some clades than others, as has been found in other invertebrates, e.g., cephalopods where a single gene order persists throughout Octopodiformes (e.g. Taite et al., in review) in contrast to extensive rearrangement in Decapodiformes (e.g., Fernández-Álvarez et al., 2022), or shrimps, where rearrangements tend to be simple in carideans, but more complex in axiideans and gebiideans, and where it has been speculated that rearrangement might be associated with specialized adaptations (Tan et al., 2017). As more mitochondrial genomes are extracted from SRAs associated with the extensive UCE dataset generated for Cnidaria (Quattrini et al., 2020), as demonstrated by Muthye et al. (2022), it may become possible to test specific hypotheses on mitochondrial gene rearrangement across Cnidaria.

Mitochondrial and nuclear phylogenies of Pennatulacea
Our phylogeny based on mitogenomic data (Fig. 3) recovered four well-supported clades (BS 100, PP 1) which, allowing for variation in taxon sampling, appear to be recovered consistently in other Pennatulacea phylogenies based on mitochondrial mtMutS and ND2 (Dolan et al., 2013, Kushida andReimer, 2019;Fig. 7), suggesting that 5 ′ MutS and 3 ′ ND2, as captured by widely-used primer pairs (Fig. 6) provide an adequate reflection of the mitochondrial genome. Nonetheless, 5 ′ MutS and 3 ′ ND2 are not the most diverse regions of the mitochondrial genomes in sea pens: the 3 ′ end of CytB and the middle regions of MutS provide at least 50 % more genetic diversity. These regions might be the best available for genetic barcoding if suitable primers are designed, and could be used to better discriminate closely related species.
Our phylogeny based on three nuclear genes (Fig. 4) also recovered four well-supported clades, but only two of these clades were the same as those in the mitochondrial tree. Alternative tree topologies have previously been discerned where mitochondrial genes have been concatenated with 28S rRNA (García-Cárdenas et al., 2020;López-González and Drewery, 2022) in a single analysis. In those studies, maximum likelihood has yielded a different tree topology than Bayesian Inference (Fig. 7), and this topological instability could be indicative of incongruence of the mitochondrial and nuclear datasets. Our conflicting, but in both cases well-resolved, clades recovered in each case (mitochondrial and nuclear) from extensive alignments is further evidence that the phylogenetic signal from the nuclear ribosomal cluster is not congruent with that from the mitochondrial genome in sea pens. Supermatrix approaches are often used to combine data, but the mitochondrial and nuclear genomes can have different evolutionary histories. Maternal inheritance and lack of recombination in the mitochondrion may cause the evolutionary signal in the mitochondrial genome to diverge from the nuclear genome (Edwards and Bensch, 2009). Such discordance has been attributed to Introgression, hybridization, and incomplete lineage sorting (e.g., Toews and Brelsford, 2012) and is not uncommon. For example, among octocorals, Frometa et al. (2021) found that nuclear 28S rRNA and mitochondrial MutS sequences were not phylogenetically congruent in work with the genus Swiftia. As the mitochondrial protein coding genes are effectively a single locus, high concordance among mitochondrial genes can lead to high bootstrap values even when the tree topology reflects the evolutionary history of the mitochondrion rather than the organism. Thus, even though our mitochondrial tree has high support we prefer to highlight where the topology conflicts with the nuclear tree rather than provide a combined analysis of mitochondrial and nuclear data. In ribosomal clusters, Fig. 6. Sliding window analysis of relative overall genetic diversity (based on Kimura-2-parameter distance) of A) Pennatulacea mitochondrial genome dataset, with regions of MutS and ND2 commonly amplified through Sanger sequencing highlighted with black boxes; and B) Pennatulacea 3-gene nuclear ribosomal dataset, with region amplified by primers 28S-Far and 28S-Rar highlighted.
despite the concept of concerted evolution, polymorphic copies of ribosomal nuclear genes can exist at significant levels in some organisms, confounding phylogenies of closely related taxa using these regions (Weitemier et al., 2015). We did not find evidence of this herein, for example where closely related or duplicate species were included in our taxon sampling, but it is a drawback of using ribosomal nuclear genes, especially if they are the only nuclear genes used.
The outstanding question, which is not answered here, is whether the nuclear ribosomal cluster is representative of the entire nuclear genome. If it is, then we have shown the utility of a short region of 28S in sampling the nuclear genome. Sequencing just 28S ( Supplementary Fig. 3) did not provide the same resolution as including all three nuclear ribosomal genes (Fig. 4). But, while the added resolution provided by higher taxon sampling (based on the numerous short sequences available on GenBank) of the partial 28S region was also not sufficient to provide complete resolution into the four expected clades, it was close, and additional sequences of this region could provide further resolution in the future. Certainly, this region, bounded by primers 28S-Far and 28S-Rar, captures the greatest overall diversity provided by the nuclear ribosomal cluster genes in sea pens (Fig. 6).
Ours are the only sea pen phylogenies to date to recover congruent support at deep nodes in both ML and BI analyses, probably reflecting the large size of our datasets in terms of alignment lengths. However, our datasets contain fewer taxa (30 versus e.g., 61 in López-González and Drewery, 2022), and increasing the number of taxa within them could increase deep node resolution further. Our mitochondrial tree supported a relationship between clades 1, 3 & 4 to the exclusion of clade 2 (BS 73, PP 1) (Figs. 3 and 7). This latter relationship is well supported in the Ultra Conserved Elements (UCE) cnidarian tree of Quattrini et al (2020). Given the amount of nuclear data in the UCE tree, it is likely that the UCE topology reflects evolutionary history, however, this was part of a much broader study and taxon sampling within Pennatulacea was limited to seven taxa (Fig. 7). The very limited taxon sampling in the UCE study means the UCE tree is also congruent with the different but also well-resolved topology of our 3-gene nuclear ribosomal tree (Fig. 7).
In fact, what is clear from Fig. 7, is that there are some key taxa whose position needs resolving. This could be achieved by sequencing other nuclear genes, for example house-keeping genes, as attempted in Porifera (Hill et al., 2013), or other high copy number genes such as H3 (Histone) as used in cephalopods (Lindgren et al., 2012) which might be retrieved from SRAs, or through the inclusion of key taxa in genome wide studies. Adding Funiculina, Gyrophyllum, Scleroptilum and Malacobelemnon to the existing UCE dataset could resolve their position, and also indicate which, if either, of the mitochondrial and nuclear ribosomal datasets is congruent with the evolution of the nuclear genome. Sanger sequencing of MutS, ND2 and 28S rRNA for sea pen genera not yet included in any trees to date (15 genera are yet to be included) would quickly reveal whether these taxa have conflicting positions in the nuclear ribosomal and mitochondrial phylogenies and should be considered as candidates for a UCE or other genome-wide study.

Convergent evolution in Pennatulacea
In the four main clades of both our mitochondrial and nuclear ribosomal trees, we observed findings similar to those already discussed in previous studies (Williams, 1995, McFadden et al., 2006, Dolan et al., 2013, Kushida and Reimer, 2019, García-Cárdenas et al., 2020 which include that: a) Pennatulacea needs revising, integrating molecular and morphological data at all systematic levels; b) traditional suborders are not monophyletic; c) there is a mix of different colony morphological forms in each clade.
The status of traditional Pennatulacea families in light of molecular data have been discussed in detail by four previous studies (Dolan et al., 2013, Kushida and Reimer, 2019, García-Cárdenas et al., 2020, López-González and Drewery, 2022. To date, at least half of the sea pen families have been recovered as para-or polyphyletic taxa and the polyphyly of several of these are evident in our phylogenetic trees. Phenotypic convergence is clearly evident in molecular phylogenetic work on Pennatulacea (Dolan et al., 2013, Kushida and Reimer, 2019, García-Cárdenas et al., 2020, López-González and Drewery, 2022, with multiple morphological forms apparent in each clade. This suggests that pennatulacean morphology has adapted strongly to environmental pressures, perhaps to cope with the pressure of mass extinctions (Quattrini et al., 2020) and possibly predation, and other environmental parameters that could drive adaptation such as depth, hydrodynamics, and food supply. Sea pens have two characteristics unusual for Anthozoa that may have allowed them to expand into new niches: their adaptation to soft substrates, and their ability to move. Sea pens have become very well adapted to soft substrates which allows them to dwell in habitats not used by most other corals. Yesson et al. (2012) conducted global habitat suitability analyses for cold-water octocorals and predicted that pennatulaceans (suborder Sessiliflorae) had the widest potential habitat range. Another advantage is that sea pens are not a completely sessile organism; they are able to burrow and also to release themselves in the current and move (Chimienti et al., 2018). These two adaptations (soft substrate and locomotion) potentially provide a great advantage in their expanding niche and avoiding unfavourable environmental conditions, dwindling food resources, and possibly even predators. Similar characteristics contributed to crinoids becoming hugely successful in the Jurassic following the evolution of Comatulidae which could move, and, through the development of their cirri, exploit novel substrates (Meyer and Macurda, 1977).

Conclusions
Significantly from this study we find that neither mitochondrial genomes, nor nuclear ribosomal phylogenies, nor morphological taxonomy of the extant pennatulacean fauna, tell the same evolutionary story. We find congruence in some parts of the trees generated from mitochondrial and nuclear ribosomal datasets, but the placement of certain taxa in these trees is starkly different. The large amounts of sequence data generated herein produce similar topologies to analyses limited to one to four gene regions, but provide more support at deeper nodes even though taxon sampling is lower than in many studies relying on Sanger sequencing. We show that in the mitochondrial genome the 3 ′ end of ND2 and the 5 ′ end of MutS capture considerable genetic diversity in sea pens, but that greater genetic diversity might be captured by primers placed in the middle of MutS or at the 3 ′ end of CytB, and that in the nuclear ribosomal cluster, existing primers capture the highest genetic diversity. We make suggestions for future molecular studies that could help resolve sea pen systematics.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.