Phylogenetics and Mitogenome Organisation in Black Corals (Anthozoa: Hexacorallia: Antipatharia): An Order-Wide Survey Inferred From Complete Mitochondrial Genomes

Black corals (Anthozoa: Antipatharia) are an ecologically and culturally important group of deep-sea cnidarians. However, as the majority of species inhabit depths >50 m, they are relatively understudied. The inaccessibility of well-preserved tissue for species of interest has limited the scope of molecular analysis, and as a result only a small number of antipatharian mitochondrial genomes have been published. Using next generation sequencing, 18 complete and five partial antipatharian mitochondrial genomes were assembled, increasing the number of complete mitochondrial genomes to 22. This includes species from six antipatharian families, four of which were previously unrepresented, enabling the first family-level, full mitochondrial gene analysis over the whole order. The circular mitogenomes ranged in size from 17,681 to 21,669 bp with the large range in size due to the addition of an intron in COX1 in some species and size variation of intergenic regions. All mitogenomes contained the genes standard to all hexacoral mitogenomes (13 protein coding genes, two rRNAs and two tRNAs). The only difference in gene content is the presence of the COX1 intron in five families. The most variable mitochondrial gene is ND4 which may have implications for future barcoding studies. Phylogenetic analysis confirms that Leiopathidae is sister to all other families. Families Antipathidae, Cladopathidae and Schizopathidae are polyphyletic, supporting previous studies that call for a taxonomic revision.


INTRODUCTION
The order Antipatharia Milne Edwards and Haime, 1857, also known as black or thorny corals, is an ecologically and culturally important group of deep-sea anthozoans in the hexacoral subclass (Cnidaria: Anthozoa: Hexacorallia). The order includes 273 species in 45 genera and seven families (Molodtsova and Opresko, 2020), with many new species expected. Black corals are found in every ocean occupying depths from 2-8,600 m (Wagner et al., 2012), but are primarily concentrated at depths >50 m (Cairns, 2007), with ∼50% of known genera inhabiting depths >500 m (Molodtsova and Opresko, 2017). Anipatharians are exclusively colonial and exhibit a diverse range of polyp, skeletal and spine morphologies. They range in height from 5 to 10 cm to >3 m (Wagner et al., 2012). Unlike scleractinian corals, the skeletons of antipatharians are not calcified but are instead rich in protein and chitin [(Bo et al., 2012b); see Ehrlich (2010) for a review] and often appear black when visible. Although not reef forming, antipatharians can play a crucial ecological role in forming coral garden habitats (Bo et al., 2009b) and individuals can support significant epifaunal populations (Love et al., 2007;Roberts et al., 2009). From China to Hawaii, black coral has been utilized in traditional medicine for a wide variety of therapeutic and spiritual uses, while its common use in jewelery has seen commercial harvesting established in the Caribbean, Hawaii and parts of Asia (Wagner et al., 2012). This has resulted in all species of black coral being listed on CITES appendix II. Their use in traditional medicine has led to interest from pharmacologists, with studies showing potential roles in cancer treatments (Qi et al., 2009;Ghandourah and Alorfi, 2019) and the presence of antibacterial properties (Al-Lihaibi et al., 2010).
The use of complete mitochondrial genomes has become an increasingly popular tool for taxonomic and phylogenetic research due to the advances and availability of next generation sequencing technology (Kayal et al., 2013;Crampton-Platt et al., 2016). In addition, metazoan mitogenomes provide a number of genome-level characters that can support phylogenetic hypotheses. These include gene order arrangements, position of introns, insertions and inferred secondary structures of transfer RNAs (tRNA) and ribosomal RNAs (rRNA) (Boore, 2006). In particular, gene order arrangements have proved useful in anthozoan studies by furthering understanding of evolutionary relationships (Brugler and France, 2007;Brockman and McFadden, 2012;Figueroa and Baco, 2014;Hogan et al., 2019) and provided additional resolution when nucleotide and protein sequence phylogenies have proved contradictory (Lin et al., 2014). This is in part due to a number of advantageous features of the mitogenome, including the presence of a near universal homology of mitochondrial genes found across the animal kingdom (Wolstenholme, 1992;Boore and Fuerstenberg, 2008). In addition, the high number of possible gene arrangements suggests that the observed gene order is most likely to be a derived synapomorphy shared among related taxa, rather than the result of a convergent rearrangement (Clayton, 1992;Falkenberg et al., 2007) [however, ancestral state reversions have been observed, see Figueroa and Baco (2014)].
Hexacoral mitogenomes share some unique properties amongst metazoans. Complete mitogenomes have been sequenced for species from each of the five hexacoral orders: Actiniaria (Foox et al., 2015;Zhang et al., 2017), Antipatharia (Brugler and France, 2007;Kayal et al., 2013;Figueroa et al., 2019), Corallimorpharia (Lin et al., 2014), Scleractinia (Medina et al., 2006), and Zoantharia (Sinniger et al., 2007;Chi and Johansen, 2017). Consistent is a reduced set of the 22 tRNAs necessary for mitochondrial translation, with only tRNA −met (Methionine) and tRNA −trp (Tryptophan) being present. In addition, group 1 introns are found within the NADH dehydrogenase subunit 5 gene (ND5) with varying numbers of embedded energy pathway genes specific to each order. A second group 1 intron within the cytochrome c oxidase subunit I gene (COX1), with an embedded homing endonuclease gene (HEG), has also been identified for some species in each of the hexacoral orders (Celis et al., 2017).
The first sequenced antipatharian mitochondrial genome, Chrysopathes formosa was circular and contained 13 energy pathway genes, two ribosomal RNA genes and a reduced number of tRNA genes consistent with all hexacorals (Medina et al., 2006;Brugler and France, 2007;Lin et al., 2014;Foox et al., 2015;Chi and Johansen, 2017). The gene order was described as most similar to the sea anemone Metridium senile, with the exception of a transposition of three continuous genes COX2-ND4-ND6 and lack of a COX1 intron. However, the subsequent fully sequenced antipatharian mitogenomes, Stichopathes luetkeni (Kayal et al., 2013), Myriopathes japonica (GenBank NC_027667.1), Tanacetipathes thamnea (Figueroa et al., 2019) and the partially sequenced Leiopathes glaberrima (Sinniger and Pawloski, 2009) had a COX1 intron containing a HEG gene. With the exception of the COX1 intron, gene content was the same as in C. formosa. The ND5 intron contains two embedded genes; the only copies of ND1 and ND3, as in actiniarians (Foox et al., 2015;Zhang et al., 2017) and zoantharians (Sinniger et al., 2007;Chi and Johansen, 2017).
Molecular approaches to examining phylogenetic relationships within Antipatharia are a relatively recent addition to systematic studies (Lapian et al., 2007;Wagner et al., 2010;Bo et al., 2012a). Nuclear ribosomal genes (rRNA) and their associated internal transcribed spacers (ITS1 and ITS2) have been effective in resolving phylogenetic relationships within Antipatharia (Lapian et al., 2007;Lapian, 2009;Bo et al., 2012aBo et al., , 2018. Mitochondrial DNA has been used in a number of antipatharian phylogenetic studies in addition to nuclear DNA, focusing on intergenic regions (IGRs) spanning trnW-ND2 and COX3-COX1 (Wagner et al., 2010;Brugler et al., 2013;MacIsaac et al., 2013) and also ND5-ND1 MacIsaac et al., 2013). Most studies have focused on clarifying the relationships within or between the families Antipathidae, Aphanipathidae, and Myriopathidae (Lapian et al., 2007;Bo et al., 2009aBo et al., , 2012aLapian, 2009;Wagner et al., 2010), while MacIsaac et al. (2013) was the first to examine relationships within Schizopathidae. To date, only one molecular study has included species from all seven antipatharian families, providing the first full family phylogeny . However, the regions examined in that study were unable to fully resolve the phylogeny, particularly at basal nodes.
This study aims to expand our knowledge of antipatharians through sequencing and examining complete mitochondrial genomes from species representing six families. Previously sequenced antipatharian mitogenomes downloaded from GenBank are re-analyzed alongside those from this study, expanding the analysis to cover all seven families.

Sampling and DNA Extraction
Antipatharian coral samples were collected from surveys around Whittard Canyon, north of the Porcupine Bank and adjacent to the Belgica Mound Province SAC in the NE Atlantic (Irish Margin) at depths of 669-2,817 m, between 2013 and 2017 using a remotely operated vehicle (ROV) Holland I onboard the R/V Celtic Explorer. Species identifications were made by examining morphological traits based on in situ and ex situ photographs with reference to classic and current taxonomical literature. Whole genomic DNA was extracted from tissue using the NucleoSpinVR TM tissue Kit (Macherey-Nagel, Düren, Germany) following manufacturer's instructions.

Sequencing Preparation
Samples were prepared for sequencing following the TruSeq DNA PCR-Free Library Preparation Kit protocol (Illumina Inc., San Diego, CA), with the exception of Tylopathes n. sp. which was prepared following the QIAseq FX DNA Library protocol (QIAGEN). For the TrueSeq samples DNA was fragmented to a target insert size of 550 bp using a Covaris Focused-ultrasonicator M220 TM . Twenty-three samples were sequenced on a single run of an Illumina MiSeq (MiSeq Reagent Kit v3 600-cycle).

Clean up and Assembly
Read quality was assessed using FastQC (Andrews, 2010). Illumina TruSeq index adapters were trimmed and low quality reads (Q<28) removed using cutadapt 2.3 (Martin, 2011). Mitogenome assembly was carried out de novo in MITObim 1.9.1 (Hahn et al., 2013) utilizing the COX1 seed sequence from Chrysopathes formosa (GenBank accession number: NC_008411.1). Complete circular read coverage was confirmed by mapping the clean reads onto the assembled circular sequence in Geneious Prime 2019.2.1 (Kearse et al., 2012) (see Supplementary Material for additional detail).

Genome Annotation
Protein coding regions, tRNAs and rRNAs were identified in Geneious using the live annotate function based on three published antipatharian mitogenomes (GenBank accession numbers: Chrysopathes formosa NC_008411.1, Stichopathes luetkeni JX023266.1, Myriopathes japonica NC_027667.1). Additionally, confirmation was sought using the online tool Mitos2 (Bernt et al., 2013) with NCBI RefSeq 81 Metazoa and genetic code 4 (mold, protozoan and coelenterate mitochondrial code). Annotations were manually edited using the widest possible open reading frame for protein coding regions based on genetic code 4 and cnidarian start (ATG) and stop (TAG, TAA) codons. Intron boundaries were confirmed manually using homologous regions in non-intron containing taxa. The position of tRNAs were confirmed using tRNAscan-SE (Lowe and Chan, 2016).

Genome Composition
Nucleotide composition and GC content were calculated in Geneious. Relative GC and AT skew was calculated using the skew formula of Perna and Kocher (1995). Origins of replication were predicted using the DNA walker method (Lobry, 1996) within GraphDNA (Thomas et al., 2007) which assesses for abrupt switches in polarity in the asymmetry of strand specific frequencies of nucleotide composition. All genomes were assessed for tandem repeats using Tandem Repeats Finder (Benson, 1999) and inverted repeats using EMBOSS einverted (Rice et al., 2000).

Nuclear DNA
Nuclear ribosomal DNA spanning the region partial 18S rRNA, ITS1, 5.8S rRNA, ITS2 and partial 28S rRNA were assembled from the cleaned reads using a combination of MITObim (see above) and read mapping onto a reference sequence (Chrysopathes formosa, GenBank accession number MG023168.1).

Alignment for Phylogenetic Analysis
The 13 mitochondrial protein coding genes and five nuclear regions were aligned separately in MAFFT v7.407 (Katoh and Standley, 2013) using the L-INS-i method; maxiterate: 1000, as this showed the highest pairwise identity (PI) and percentage of identical sites (IS). Each alignment was processed through Gblocks 0.91b (Castresana, 2000) which improves alignments for phylogenetic analysis by removing poorly aligned and divergent regions (see Table S1 in Supplementary Material for Gblock parameters). Concatenated sequences were created for (a) the 13 mitochondrial genes; (b) five nuclear regions; and (c) combined mitochondrial and nuclear regions, using a seqCAT.pl (Bininda-Emonds, 2005) script. Published genomes used in this study were re-analyzed alongside the newly sequenced genomes.

Phylogenetic Analysis
Three phylogenies were constructed from the concatenated sequences; (a) the 13 mitochondrial genes (30 taxa; 12,227 characters); (b) five nuclear regions (29 taxa; 10,28 characters); (c) combined mitochondrial and nuclear regions (29 taxa; 13,255 characters), using maximum likelihood (ML) implemented in IQ tree (Nguyen et al., 2014) and Bayesian inference (BI) statistics implemented in MrBayes 3.2.6 (Ronquist et al., 2012). Due to the lack of a suitable nuclear sequence, Stichopathes luetkeni and Tanacetipathes thamnea were left out of the nuclear dataset phylogeny. As nuclear regions for Myriopathes japonica were unavailable, the congeneric Myriopathes myriophylla (GenBank accession number AM404328.1) was used to create a chimeric sequence for Myriopathes in the combined dataset phylogeny, and as the single representative of Myriopathes in the nuclear phylogenies. Edwardsia timida replaced Nematostella sp. within the outgroup of the nuclear phylogeny, as suitable nuclear sequences were unavailable for Nematostella sp.. A chimeric sequence consisting of nuclear regions for E. timida and 13 mitochondrial genes of Nematostella sp. represented the family Edwardsiidae as an actiniarian outgroup in the combined phylogeny. Models of sequence evolution were determined for each Gblock-edited region using jModelTest 2.1.10 (Guindon and Gascuel, 2003;Darriba et al., 2012) based on AIC scores, and implemented in Bayesian analysis (see Table S2 for models). Partitions for each region were also input to IQtree which runs an automatic model testing and selection using ModelFinder (Kalyaanamoorthy et al., 2017) (Table S2). Both IQtree and MrBayes analyses were carried out on the CIPRES Science Gateway V. 3.3 (Miller et al., 2010). For ML analysis, 1,000 bootstrap replicates were performed within IQtree using the ultrafast bootstrap (Hoang et al., 2017). Bayesian analysis ran for 10 million generations, sampling every 1000 generations on two separate runs and four chains. All runs were visualized in Tracer v1.7.1 (Rambaut et al., 2018). All runs showed an effective sample size >6,000. A burn-in of 1 million was applied and a single target tree obtained using TreeAnnotator v1.10.4 (Bouckaert et al., 2014). Tests of monophyly were performed by running constrained analyses fixing the monophyly of selected groups and comparing likelihoods with the unconstrained tree using an approximately unbiased (AU) test (Shimodaira, 2002).

Outgroup Selection
We have included a representative from each of the four additional hexacoral orders, including three from Actiniaria as outgroups. Phylogenies are rooted at a midpoint between Antipatharia and the hexacoral outgoups to examine Antipatharia monophyly rather than examine hexacoral relationships. (see Table S3 for GenBank accession numbers). Phylogenies with different outgroup combinations were explored to assess ingroup topology consistency (see Figures S1-S3).

Genome Organisation
Eighteen complete and five partial mitochondrial genomes of 23 deep-sea antipatharian morphospecies from the Northeast Atlantic, were assembled and annotated ( Table 1). These represent 14 genera and six families. The 18 complete genomes are circular and range in size from 17,679 (Tylopathes n. sp.) to 21,669 bp (Leiopathes cf. glaberrima) with the size range primarily due to the addition of the COX1 intron and size variation of IGRs. Partial mitogenomes were recovered for five morphospecies. Two of these were almost complete, with Leiopathes montana missing a 156 bp gap (predicted based on genome similarity to L. cf. glaberrima) and Antipathes cf. dichotoma missing a 9 bp gap (predicted based on genome similarities to the three additional Antipathidae samples) and possessing 16 ambiguous sites. Three samples (Chrysopathes cf. micracantha; Cladopathes cf. plumosa; Trissopathes grasshoffi) yielded very limited coverage, resulting in the recovery of 3.2, 48.2, 1.9%, respectively, of a full-size genome (based on the full genome size of Trissopathes cf. tetracrada and Chrysopathes formosa).
All recovered genomes contain 13 energy pathway genes (COX1-3, ND1-6, ND4L, CYTB, ATP6, and ATP8), two tRNAs (tRNA −met and tRNA −trp ) and two subunit rRNAs (12S and 16S) ( Table 2). Genomes from eight samples are interrupted by a group 1 intron in COX1 at the same insertion site, the '3 end of position 884 [see Emblem et al. (2011) for intron nomenclature]. Additionally, all taxa with COX1 introns contain a HEG gene which varies in length from 657 (Myriopathes japonica, Tanacetipathes thamnea and Tylopathes n. sp.) to 1,062 bp (Stichopathes luetkeni). All genomes contain a group 1 intron in ND5 with the same insertion site at the '3 end of position 717, containing the only genome-wide copies of ND1 and ND3. All genes are separated by IGRs which show large size variation ( Table 2; see Table S4 for full species table).
Across all samples the largest IGR is found between ND4 and ND6 in Phanopathes sp. and extends for 597 bp. However, the range for this region between samples is 531 bp, extending for only 66 bp in C. formosa. The smallest IGR is found between ND6 and ATP8 in M. japonica, T. thamnea and Tylopathes n. sp. and extends for 4 bp. The genomes with the largest total intergenic space are L. cf. glaberrima and L. montana both with a length of 2912 bp, while M. japonica has the smallest at 716 bp.
A number of novel ORFs were identified within the IGRs (see Table S5) The largest ORF is found exclusively in Phanopathes sp. in the IGR ND4-ND6 spanning 477 bp. No single ORF was found consistently across all samples, however within families there are conserved ORFs of the same length and at the same location.
Gene content within the sample genomes shows two variations (Figure 1), with the only variation due to the absence of the HEG gene in the families Schizopathidae and Cladopathidae. All genes are transcribed in the same direction and on the same strand.
The DNA walk method found four abrupt changes in base composition across all genomes (Figure 2). These are found in the intergenic regions ND2-12S, 12S-CYTB, trnM-16S and 16S-COX3. An additional reversal in polarity is present only in samples with the COX1 intron, marking the boundaries of the intron. These junctions suggest the genome is split into four blocks (Figure 1); block 1: COX3-ND2, block 2: 12S, block 3: CYTB-trnM, block 4: 16S, and an additional block for samples with the COX1 intron; block 5: COX1 intron. Three of the IGR which show abrupt changes in polarity are rich in A+T content across all genomes; ND2-12S A+T range of 60.5 to 70.5%, trnM-16S A+T range of 58.4 to 71.4% and 16S-COX3 A+T range of 56.4 to 69.4% (data not shown). The only tandem repeats (>15 bp) were a 37 bp sequence found in the IGR COX2-ND4 in two samples S. abyssicola and A. cf. dichotoma, and a 22 bp repeat found in M. japonica within the COX1 intron. Inverted repeats  (>15 bp) were found in six samples (Table S6), with the only consistent repeats being found in the three Leiopathes samples.

Protein Coding Genes and rRNAs
Protein coding gene length was similar across all genomes (

Phylogenies
Phylogenetic analyses (mitochondrial, nuclear and combined) produced largely congruent trees (Figure 3; see Figures S4-S7 for nuclear regions and combined trees). Various combinations of hexacoral outgroups (see Figures S1-S3) confirm the monophyly of Antipatharia and maintain the same ingroup topology. The family Leiopathidae formed a monophyletic group sister to all other antipatharians. The families Antipathidae (in the mitochondrial and combined phylogenies only), Cladopathidae and Schizopathidae (due to the placement of Sibopathes cf. macrospina nesting within a clade containing three Parantipathes spp) did not form monophyletic groups. The genera Bathypathes, Parantipathes, Stauropathes, Stichopathes and Trissopathes did not form monophyletic groups. Constraining the monophyly of each polyphyletic group produced a significantly worse likelihood value tree compared to the unconstrained tree (P < 0.005; see Figure S8).

DISCUSSION
This study increases the total number of complete antipatharian mitochondrial genomes represented in the literature from four to 22, with five additional partial sequences. This includes species from four families previously unrepresented, enabling for the first time, a family-level full mitochondrial gene analysis over the whole order.

Genome Organization, Gene Order, and Variable Regions
Based on our analysis, antipatharian mitogenomes have the same features typical of all published hexacoral order mitogenomes (Lin et al., 2014;Chi and Johansen, 2017;Zhang et al., 2017). Gene order and content are highly conserved across the antipatharian order, consistent with previous findings (Brugler and France, 2007;Kayal et al., 2013) with the absence of the COX1 intron as the main notable variation. Internal gene rearrangements have been reported within all hexacoral orders (Emblem et al., 2011;Lin et al., 2014;Foox et al., 2015), but it appears that antipatharians are an exception, along with zoantharians (Poliseno et al., 2020). The stable structural arrangement in these two orders may be linked to the slow rate of nucleotide substitution in anthozoans (Hellberg, 2006) . Metazoan mitochondrial DNA, and specifically COX1, was originally viewed as an attractive target for barcoding as it was thought to be almost entirely free of introns and showed a greater range of diversity than other mitochondrial genes (Hebert et al., 2003). In the present analysis, COX1 has more identical sites and a higher PI (92.8% PI; 75.6% IS), than seven other mitochondrial genes (see Table 2). ND2 varies most (91.4% PI; 70% IS), however, the presence of indels, ranging from 6 to 20 bp at multiple locations, drives this variability. The highest variability seen in an  , future researchers may wish to consider the addition of ND4 to further increase the resolving power of taxonomic identification.

Origins of Replication and Intergenic Regions
The origins of replication in mitogenomes are often associated with large intergenic spaces, inverted repeats, A+T rich regions and at the junction of conserved gene blocks (Pearson et al., 1996;Brugler and France, 2008;Hogan et al., 2019). Additionally, mitogenomes have two separate origins of replication (Clayton, 1992), one found on each strand; OriH found on the heavy strand and OriL found on the light strand. Brugler and France (2007) investigated the ND2-12S IGR as a possible location of OriH due to its large size. With only one complete antipatharian genome for their analysis, comparisons were made with other hexacoral species, and they ultimately rejected this location as the OriH. However, with the larger data set provided by the current study, we note abrupt compositional changes around this region, the presence of inverted repeats and relatively high conservation (aligning all ND2-12S IGRs produces a 278 bp region with 67.3% PI), all of which suggest the IGR ND2-12S is the most likely candidate for the OriH in Antipatharia.
Defining the exact location of OriL has proved difficult in invertebrates (Wolstenholme, 1992). Brugler and France (2008) proposed candidate regions between 12S-CYTB and trnM-16S for Chrysopathes formosa based on the abrupt changes in base composition at these locations. In the present study 12S-CYTB has a relatively balanced A+T/G+C base content (average 48.6% A+T) and a relatively small IGR (mean 86 bp). In contrast the IGRs trnM-16S and 16S-COX3 both show an abrupt change in base composition and have a high A+T content (mean 65.4 and 64.5%, respectively), while IGR 16S-COX3 has a relatively large IGR (mean 136 bp), making it the preferred candidate for OriL. Further investigation into the secondary structures of the IGRs, including presence of A+T stable stem loops, would be necessary to provide additional evidence for the location of OriH and OriL.
Open read frames of unknown function have not previously been reported for antipatharian mitogenomes. However, they are common in other hexacoral orders (Flot and Tillier, 2007;Emblem et al., 2014;Foox et al., 2015;Chi and Johansen, 2017). Their functionality in actinarians has been speculated to be that of transposon-like elements capable of influencing nearby genes, as possible enhancers or promotors (Zhang et al., 2017). Emblem et al. (2014) demonstrated that ORFs within FIGURE 2 | DNA Walk [GraphDNA (Thomas et al., 2007)] results from vectorial representation of a linearized genome beginning at COX1 (circular symbol) and ending at IGR COX3-COX1 (star symbol). Each of the four nucleotides in the genome sequence is designated a direction as displayed on the compass. Units on axis are cumulative skew measured as number of nucleotides X 10 3 (Window size is displayed at 400 for ease of viewing). (A) Genomes with the COX1 intron. (B) Genomes with the COX1 intron absent. The dotted line represents the general shape of the walk. Abrupt changes in direction are characteristic of a separate evolutionary origin for that particular region. actinarian mitogenomes are likely functional units. The lack of conserved ORFs across antipatharians suggests that many are likely to be small pseudogenes. However, small proteins, often defined as <100 amino acid residues (Su et al., 2013), have been shown to be responsible for important biological processes in other metazoans (Galindo et al., 2007). The largest ORF found in this study, 477 bp in Phanopathes sp. is the most likely candidate for a functional protein due to its size. Transposonlike elements have been demonstrated to play important roles in adaptation to environmental stress (Casacuberta and González, 2013). Considering the environmental extremes of antipatharian deep-sea habitats, the function of these ORFs may well be similar to transposon-like elements offering beneficial metabolic adaptations to the deep-sea.
The potential for non-homologous cross mapping of nuclear DNA of mitochondrial origin (nuMT) to mitochondrial DNA (mtDNA) assemblies [see Maude et al. (2019) for a discussion] was considered in light of the large number of varied ORFs found throughout the mitogenomes. BLAST analysis of the larger ORF (>70 bp) revealed no matches to hexacoral nuclear regions. Read coverage for ORF loci was consistent with coverage for the overall mitogenome (i.e., if similar copies were present in the nuclear genome, we might expect higher coverage). Additionally, we might expect higher sequence variability in these ORFs compared with other coding regions (e.g., greater disagreement of mapped reads in these ORFs), however this was not the case.

Introns
A cyclic evolutionary model of intron gain and loss has previously been proposed for the presence of the COX1 intron in cnidarians (Goddard et al., 2006). Celis et al. (2017) suggest the COX1 intron was inserted early in cnidarian history and rapidly invaded intron-less alleles via horizontal transfer due to the HEG gene. After becoming fixed in the population, it was subject to mutational degradation, losing HEG functionality, but continued to be inherited vertically as a HEG remnant. The effects of purifying selection resulted in the complete deletion of the intron and HEG from COX1. Support for this hypothesis in Antipatharia is evident in the variation of HEG sizes observed across families, with the presence of shorter HEG genes in Myriopathidae and Stylopathidae (657 bp) indicating they are in a later stage of the degradation period as compared to Antipathidae (1,062 bp), Aphanipathidae (981 bp) and Leiopathidae (996 bp). The absence of the COX1 intron within the Cladopathidae and Schizopathidae FIGURE 3 | Phylogenetic tree reconstruction based on all available (newly sequenced from this study and GenBank samples) antipatharian mitochondrial genomes using the 13 energy pathway genes based on a Bayesian inference (BI) and maximum likelihood (ML) analysis. The trees are rooted at a midpoint between Antipatharia and the hexacoral outgroups (See Figure S1 for outgroup relationships). BI posterior probabilities and ML bootstrap values shown at nodes (BI/ML). The arrow indicates the apparent loss of the COX1 intron at an internal node.
Frontiers in Marine Science | www.frontiersin.org suggest a complete evolutionary cyclic history of intron invasion, degradation and deletion has occurred within antipatharians, as has previously been shown in Actinaria (Emblem et al., 2014), and in Scleractinia and Corallimorpharia (Celis et al., 2017). This is further implied by the phylogenetic analysis which suggests the COX1 intron was present in the ancestral state due to its presence in Leiopathidae, which is sister to all other antipatharian families. It was then subject to a single loss in the branch leading to the Cladopathidae + Schizopathidae clade (see Figure 3), suggesting a close evolutionary relationship between these two families as shown by Brugler et al. (2013). Interestingly, species of Cladopathidae and Schizopathidae are the only antipatharians found in the abyssal zone >3,000 m (Molodtsova and Opresko, 2017). Whether the loss of the COX1 intron is an adaptation to the deep-sea environment warrants further investigation.
In contrast to the COX1 intron, the ND5 intron is present within all sequenced genomes. The insertion position is at the same location as found in all hexacoral orders (Chi and Johansen, 2017), referred to as position 717 (Emblem et al., 2011). Emblem et al. (2011) propose that the ND5 intron originally contained a HEG gene and was subject to horizonal transfer and a similar cyclic scenario as the COX1 intron. However, at one stage, the insertion of ND1 and ND3 resulted in the intron becoming a compulsory element, resulting in a vertical pattern of inheritance, supported by the conserved intron position. In contrast, three insertion sites have been identified for the COX1 intron across the hexacorals (Celis et al., 2017), suggesting a lack of positive selection.

Phylogenies and Interfamilial Relationships
Previous molecular studies show consistent congruent sister groupings between the following sets of antipatharian families: Cladopathidae + Schizopathidae Bo et al., 2018), Myriopathidae + Stylopathidae , Antipathidae + Aphanipathidae (Lapian et al., 2007;Bo et al., 2009aBo et al., , 2018Lapian, 2009;Wagner et al., 2010;Brugler et al., 2013). The phylogenies presented in this study are in agreement with these family groupings. Based on morphology, Opresko (1998) considered Leiopathidae to be a distinctive family and suggested it may warrant a higher level classification, while Brugler et al. (2013) suggested it may be the most primitive family. However, molecular studies have shown conflicting results with both Antipathidae and Leiopathidae being recovered as sister to the remaining antipatharian families Bo et al., 2018). As the ingroup topology is sensitive to the choice of outgroup and number of outgroup samples (Smith, 1994), we have chosen to include three samples from a close sister order, Actiniaria, and a single representative from each of the three remaining hexacoral orders, Corallimorpharia, Scleractinia and Zoantharia. This study presents the largest molecular character set to date (13,255 characters) and supports Leiopathidae as sister to all other families, reconciling the morphological/molecular positions. This study should be considered a preliminary assessment of the interfamilial relationships within Antipatharia until broader geographic and taxonomic representation can be achieved using full mitogenomes.
The species reported here as Sibopathes cf. macrospina (family Cladopathidae) clustered within a clade with three Parantipathes species (family Schizopathidae) with marginal genetic differences. Similar results were reported by Brugler et al. (2013) for material from the Gulf of Mexico and Western Atlantic. This suggests S. cf. macrospina may indeed belong to the genus Parantipathes. There are a number of morphological differences reported between S. cf. macrospina and the type species of the genus Sibopathes [S. gephura van Pesch, 1914], that include pinnulation pattern, form of spines and morphology of polyps (van Pesch, 1914;Opresko, 1993;unpubl. observations T.M.). Molecular studies including S. gephura are necessary to validate the phylogenetic position of the genus Sibopathes and resolve the apparent polyphyly in Cladopathidae and Schizopathidae. In addition, species in the genera Trissopathes and Cladopathes formed a well-supported clade with C. cf. plumosa and T. cf. tetracrada in particular showing almost no evolutionary distance between them in the nuclear phylogenies (Figures S4, S5). Further studies may conclude that they represent the same genus, however it is significant to note that the mitogenomes of C. cf. plumosa and T. grasshoffi are currently incomplete.
Polyphyly in Antipathidae has been noted in previous studies Bo et al., 2018) with the suggestion that a taxonomic revision is necessary. In particular the genera Antipathes and Stichopathes have been posited for revision and splitting (Bo et al., 2012a(Bo et al., , 2018. The present study found the single Antipathes sample A. cf. dichotoma, formed a 100% BS (BPP=1) supported clade with S. abyssicola, with little evolutionary distance between them. A similar lack of phylogenetic signal between the genera Antipathes and Stichopathes has been reported by previous studies (Bo et al., 2012a(Bo et al., , 2018Brugler et al., 2013). Bo et al. (2012a) highlights that the single-stem morphology (attributed to Stichopatheslike specimens), is an unreliable characteristic as it may have evolved independently in different taxa. This would account for the lack of grouping between the three Stichopathes species in this study.
With the current study, only one sample is representative of Aphanipathidae, Phanopathes sp., which forms a wellsupported clade with the antipathid, Stichopathes n. sp.. This is consistent with the findings of previous studies Bo et al., 2018), which have found Phanopathes species clustering with A. dichotoma and various Stichopathes species. With only one full genome representing Aphanipathidae, the current analysis is limited in scope, however the high support values for the Phanopathes sp. + Stichopathes n. sp. clade is noteworthy, although the large branch lengths indicate high genetic differentiation between these species. Bo et al. (2018) suggests that the name bearing type for the genus Antipathes and ultimately the family Antipathidae, A. dichotoma Pallas, 1766, is likely to be unrepresentative of the genus and therefore the family. Our study supports the need for taxonomic revision within Antipathidae and Aphanipathidae Bo et al., 2018), which would benefit from further sequencing of full mitogenomes representative of a greater range of diversity within these families.
The polyphyly and particularly short genetic distances found within two genera of the Schizopathidae samples; Bathypathes and Stauropathes, support the morphological similarities observed between Telopathes, Stauropathes and Bathypathes (MacIsaac et al., 2013). They form a tight clade in the current and previous studies (MacIsaac et al., 2013;Bo et al., 2018), warranting further morphological analysis to validate these genera.

Conclusions and Future Investigations
Antipatharia appears to have a stable gene order and content, with the exception of the absence of the COX1 intron in two families. The most variable gene region is found to be ND4, in contrast to the commonly recommended COX1 'barcode of life'. The variability in size range of IGRs, presence of multiple ORFs and HEG gene length variations across families, offer an avenue for future studies to investigate additional evolutionary scenarios and possible adaptations to metabolic functions. Phylogenetically we find that Leiopathidae is sister to all other antipatharian families. Families Antipathidae, Cladopathidae and Schizopathidae are polyphyletic, and several genera may be in need of taxonomic revision.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
The project was conceived by AA, RH, and CY and designed by CY with input from NB. Sample collection, initial taxonomic identification, and DNA extraction was carried out by AA and RH, with the exception of Sibopathes cf. macrospina, which was collected by AW. Final taxonomic identification was carried out by TM. DNA library preparation for sequencing was carried out by KH and NB. The pipeline assembly including data cleaning, genome assembly, annotation, and all downstream analysis was carried out by NB with advice from CY. The manuscript was written by NB with contributions from all authors. All authors edited the manuscript before submission.