Rates and relations of mitochondrial genome evolution across the Echinoidea, with special focus on the superfamily Odontophora

Abstract In order to better characterize the placement of genus Tripneustes, as a representative of the Toxopneustidae family within the broader sea urchin mitochondrial (MT) phylogeny, the complete MT genome of Tripneustes gratilla was generated and compared with all published echinoid MT genomes currently available on NCBI GenBank. The MT genome phylogeny supports the existence of the superfamily Odontophora (consisting of the families Strongylocentrotidae, Echinometridae, and Toxopneustidae). A relaxed molecular‐clock time calibration suggests a split between the three key Odontophore MT lineages occurred during the late Eocene/Oligocene. Major global oceanographic changes have been inferred during this time frame, potentially driving species diversification through environmental selection pressures. To test for signatures of selection acting on the mitochondria, the historical rate of gene evolution of individual MT genes was assessed through a branch‐site comparison of nonsynonymous to synonymous substitution ratios (ω). Models of positive selection and neutral evolution, as compared via a likelihood ratio test, show no evidence of strong historical positive selection on mitochondrial genes at the genesis of the Odontophora. However, while pairwise ω comparison revealed signatures of strong negative selection, relatively elevated ω values were observed within the Strongylocentrotus genus.


| INTRODUCTION
There are presently over 1,000 described species of sea urchins that collectively make up the class Echinoidea, the sea urchins, within the Phylum Echinodermata (WoRMS, 2017). Arising early on in the fossil record (Ordovician) as a distinct class, and lending themselves well to fossilization, echinoids have long been a hallmark of systematic research (Agassiz & Clark, 1907-1917Clark, 1912;Fell, 1974;Kroh & Smith, 2010;Littlewood & Smith, 1995;Smith, 1988;Smith & Savill, 2001). Following the publication of Theodore Mortensen's detailed alpha taxonomy of echinoids between 1928 and 1951, sea urchins became a tractable system to investigate questions of evolution and speciation (Mayr, 1954;Mortensen 1928-1954, Palumbi & Lessios, 2005. Within the Echinoidea resides the ominously named Toxopneustidae family, whose members include the highly venomous "flower urchin," Toxopneustes pileolus (Nakagawa et al., 1996), the much-studied Caribbean "green sea urchin," Lytechinus variegatus (Watts, McClintock, & Lawrence, 2013), and the Indo-Pacific "collector urchin," Tripneustes gratilla (Lawrence & Agatsuma, 2013), along with eight other recognized species (WoRMS, 2017). Though these species are well recognized and studied, their familial phylogeny is not as well established. Previous attempts at placement of the Toxopneustidae family within the greater sea urchin phylogenetic tree using morphology have suggested that the Toxopneustidae, Echinometridae, and the Strongylocentrotidae together form the superfamily Odontophora (Kroh & Smith, 2010). The complete MT genome of T. gratilla, reported here, represents the last genome needed in order to compare all MT family lineages within the proposed Odontophora superfamily.
Markers derived from MT genomes have long served as the molecular standard for species delineation, and population connectivity, as well as deeper evolutionary relationships (Avise et al., 1987;Ballard & Whitlock, 2004). Part of the appeal of the MT marker has been not only its ease of recovery during DNA extraction owing to the high copy numbers of MT genetic material per cell, but also the lack of recombination and a general assumption of freedom from strong positive selection (Grey, 1989). However, a growing number of studies have suggested natural selection on MT genes may not be as rare as previously assumed (Ballard & Whitlock, 2004;Bazin, Glémin, & Galtier, 2006;Doi, Suzuki, & Matsuura, 1999;Stojković et al., 2016).
Genetic variation of MT genes may be driven by thermal adaptation in ectothermic poikilotherms, such as sea urchins, as the thermal stability of transcribed proteins is critically important to function (Guderley & St-Pierre, 2002;Hazel, 1995). Indeed, periods of great oceanographic change, including significant global ocean cooling, have been linked to species divergence as a result of climate-driven selection pressures (Prothero & Berggren, 1992). Beyond thermal adaptation, cytonuclear incompatibility can also serve as a selection force on MT genes, which rely on nuclear molecular machinery. Cytonuclear incompatibility may arise through population divergence, as nuclear and MT genes evolve within a population, and can serve as a mechanism to restrict gene flow during secondary contact between populations (Burton & Barreto, 2012).
One approach to detect molecular signatures of selection across lineages is to quantify the proportion of nucleotide substitutions occurring at nonsynonymous codon sites compared to the proportion of synonymous substitutions, accounting for the degree of degeneracy at each codon site. This ratio of nonsynonymous to synonymous substitutions (dN/dS or ω) can serve as a highly conservative estimate of selection when comparing lineages. An assumption is made that the majority of synonymous substitutions are selectively neutral, while the occurrence of nonsynonymous substitutions are presumably selection driven. Values of ω close to 1 suggest neutrality, values much less than 1 suggest the action of purifying selection in removing amino acid substitutions, and values much greater than 1 pointing to positive selection on amino acid change (Kryazhimskiy & Plotkin, 2008;Nielsen & Yang, 2003). Extending this approach to inferred historical sequences of shared ancestors allows for an estimation of past instances of selection-driven divergence. This study does just that by estimating ω values for each branch of a phylogenetic tree generated from complete MT genomes of 14 taxa across the Echinoidea, including 10 from the Odontophora. Two branch-site tests of positive selection, a "strict" and "relaxed" variety, were performed on the branch giving rise to the Odontophora, in order to assess the influence of climate shifts driving lineage divergence through positive natural selection.
However, due to ambiguous phylogenetic signals across the T. hardwickii MT genome and a polytomous assignment of the species, it was dropped from the final analysis.
Full MT genome alignments were performed with MAUVE as implemented in Geneious v.6 (Darling, Mau, & Perna, 2010). Sequence length comparisons of the 15 MT genes were also performed in Geneious v.6. The two noncoding 12s-rRNA and 16s-rRNA sequences were aligned in T-Coffee using the ribosomal folding algorithm implemented in the rcoffee mode (Notredame, Higgins, & Heringa, 2000;Wilm, Higgins, & Notredame, 2008). Comparative alignments of the rRNA sequences were executed with the Muscle (Edgar, 2004) and ClustalW (Thompson, Higgins, & Gibson, 1994) algorithms, as implemented in Geneious v.6, under default settings. The 13 coding DNA sequences (CDS) were aligned using the amino acid alignment algorithm in Geneious v.6, using the echinoderm MT codon table. Best fit nucleotide substitution models for the initial phylogenetic analysis for each of the 15 gene alignments were determined with JModeltest 2.1.7 (Darriba, Taboada, Doallo, & Posada, 2012).
From the 15 gene alignment dataset, a phylogenetic tree was generated with MrBayes 3.2.3, which was run for 5,000,000 steps; with a 1,250,000 step burn-in (Huelsenbeck, Ronquist, Nielsen, & Bollback, 2001;Ronquist & Huelsenbeck, 2003), with A. lixula designated as an out group in order to root the tree. To confirm the robustness of the topology, a phylogenetic tree was also generated from the same data using RAxML v.8, which was run for 50,000 steps, with nonparametric bootstrapping enabled (Stamatakis, 2014). For the time calibration and gene evolution analysis a second Bayesian majority consensus tree was generated, under identical parameters as the prior tree, but including only the 13 CDS, not the 12S and 16S sequences. This CDS tree was generated to prevent confounding of codon substitution rate assessments with noncoding gene sequences.

| RESULTS
The   Table 1  A. lixula, the outgroup, displays the split between Order Arbacioida and the Camarodonta (Figure 1 Averaging across all gene pairwise comparisons, the highest relative ω values were consistently observed between S. purpuratus and its three congenerics. Average pairwise ω is summarized in Figure 2. Testing specifically the Odontophore branch, no CDS had significant results of positive selection according to the "strict" branch-site test, with only three genes (ND1, ND5, and COX1) exhibiting differences greater than zero. Of those three, only ND1 showed a significant difference at the 5% level in the "relaxed" test for selection (Table 2).

| DISCUSSION
The placement of the Toxopneustidae family, represented here by the Tripneustes gratilla MT genome, as a closer sister group to the Strongylocentrotidae than the Echinometrid Heliocidaris crassispina contradicts the accepted species tree (Kroh & Smith, 2010). It must be specified that this relationship is based solely on MT sequences and is thus an insight into the evolutionary relationships of the   (Lee et al., 2004). The absence of regular recombination, generally assumed with mitochondria, can significantly inhibit the rates of fixation of advantageous mutants (Hill & Robertson, 1966).
This could explain the lack of a strong positive selection-driven signal. Only difference values that were greater than zero are reported. p-values are given for a 50:50 mix of point mass 0 and χ 1 2 distribution.
Urchin Genome Sequencing Consortium 2006). However, consistent pairwise ω values much below one in most all taxa comparisons suggests that the major force affecting MT sequences is purifying, or negative, selection.
The generation of a complete MT genome of Tripneustes gratilla, the first for a member of Toxopneustidae, has allowed for the confirmation of the existence of the Odontophora superfamily. This underscores the fact that MT genomes can serve as a powerful analytical unit when estimating species relationships, especially in systems that may otherwise lack publicly available nuclear sequences for comparisons. It is the hope of the author that researchers generating NGS data will view the isolation, proper annotation, and dissemination of captured MT genome sequences to be a worthwhile endeavor. Allen and two anonymous reviewers, for helpful comments and suggestions.