Phylogenomic analyses of mud dragons (Kinorhyncha)

Mud dragons (Kinorhyncha) are microscopic invertebrates, inhabiting marine sediments across the globe from intertidal to hadal depths. They are segmented, moulting animals like arthropods, but grouping with the unsegmented priapulans and loriciferans within Ecdysozoa. There are more than 300 species of kinorhynchs described within 31 genera and 11 families, however, their evolutionary relationships have so far only been investigated using morphology and a few molecular markers. Here we aim to resolve the relationships and classification of major clades within Kinorhyncha using transcriptomic data. In addition, we wish to revisit the position of three indistinctly segmented, aberrant genera in order to reconstruct the evolution of distinct seg- mentation within the group. We conducted a phylogenomic analysis of Kinorhyncha including 21 kinorhynch transcriptomes (of which 18 are new) representing 15 genera, and seven outgroups including priapulan, lor- iciferan, nematode and nematomorph transcriptomes. Results show a congruent and robust tree that supports the division of Kinorhyncha into two major clades: Cyclorhagida and Allomalorhagida. Cyclorhagida is composed of three subclades: Xenosomata, Kentrorhagata comb. nov. (including the aberrant Zelinkaderes ) and Echinorhagata. Allomalorhagida is composed of two subclades: Pycnophyidae and Anomoirhaga nom. nov. Anom- oirhaga nom. nov. accommodates the aberrant genera Cateria (previously nested within Cyclorhagida) and Franciscideres together with five additional genera . The distant and derived positions of the aberrant Zelinkaderes , Cateria and Franciscideres species suggest that their less distinct trunk segmentation evolved convergently, and that segmentation evolved among kinorhynch stem groups.


Introduction
Kinorhyncha, also known as mud dragons, are a phylum of marine, segmented microscopic invertebrates, ranging in size from 100 to 1000 µm, and belonging to Ecdysozoa (Aguinaldo et al., 1997). Ecdysozoans are moulting animals and include the clades Panarthropoda (with tardigrades, onychophorans and arthropods), Nematoida (nematodes and nematomorphs), and Scalidophora. Kinorhynchs belong to the Scalidophora, which also includes the unsegmented priapulans (penis worms) and loriciferans (girdle wearers). Even though Scalidophora is broadly accepted as a clade, merely based on morphology, its monophyly has been either questioned (Laumer et al., 2019), or only partially tested in broad molecular phylogenetic analyses with representatives of maximum two out of the three groups (e.g., Mallatt and Giribet, 2006;Dunn et al., 2008;Sørensen et al., 2008;Campbell et al., 2011;Laumer et al., 2015;Yamasaki et al., 2015, summarized in Giribet andEdgecombe (2017)). In most of these phylogenetic analyses kinorhynchs and priapulans appear as sister groups, while the position of loriciferans has varied. Segmentation has been assumed to evolve independently at least three times within bilaterians (in panarthropods, annelids and chordates) (e.g., Scholtz, 2020). Within Ecdysozoa, kinorhynchs are closely related to unsegmented groups and always found distantly related to panarthropods. Thus, their segmented body plan most likely represents a second independent evolutionary event within ecdysozoans. Most kinorhynchs show a similar body plan including a radial head bearing an introvert and a mouth cone with a terminal mouth, a neck, and a trunk divided into eleven articulated segments (Sørensen and Pardos, 2020) (Fig. 1). The general morphology differs in aberrant kinorhynchs, which have a much more elongated habitus, thin cuticle, and a less distinct external segmentation (Herranz et al., 2019b(Herranz et al., , 2021a (species marked with asterisks in Fig. 1). In order to understand the evolution of kinorhynch segmentation it is necessary to have a solid morphological and phylogenetic background. During the last decade, the morphology of aberrant kinorhynchs has been thoroughly investigated (Dal Zotto et al., 2013;Herranz et al., 2019aHerranz et al., , 2021aNeuhaus and Kegel, 2015;Yamasaki, 2019;Rucci et al., 2020); however, their phylogenetic position is still contentious, especially regarding the genus Cateria (see Sørensen et al., 2015). Accordingly, two possible scenarios emerge: i) Indistinctly segmented, aberrant forms are derived, which suggests that distinct segmentation would be synapomorphic for kinorhynchs and that aberrant traits therefore surged from modifications of a conserved, distinctly segmented pattern. ii) Aberrant kinorhynchs branch out as sister group to all remaining kinorhynchs, suggesting that their "less segmented" appearance is a plesiomorphic trait, and that distinct segmentation evolved progressively within the phylum.
Multiple attempts have been made to understand kinorhynch phylogeny using molecular data, yet, these studies either focused on a specific kinorhynch subgroup (Pycnophyidae in Sánchez et al., 2016) or suffered from a very limited taxon sampling (e.g., Dal Zotto et al., 2013;Yamasaki et al., 2013). The most complete kinorhynch phylogeny was based on a combined approach using morphological and molecular data, yet only including two molecular markers (18S and 28S rRNA) . This study established a new classification, dividing kinorhynchs into two major clades, Cyclorhagida and Allomalorhagida. Cyclorhagida contained the three subclades Echinorhagata (=Echinoderidae), Kentrorhagata, and Xenosomata (=Campyloderidae). Allomalorhagida included Dracoderes (a genus previously considered as part of Cyclorhagida) as sister group to a clade accommodating all traditional homalorhagids together with the aberrant Franciscideres and Gracilideres . Interestingly, other aberrant genera such as Cateria and Zelinkaderes remained within Cyclorhagida. However, the position of Cateria was based exclusively on morphological data.
In the present study, we carried out the first phylogenomic analysis of kinorhynchs. Our dataset includes transcriptomes from 28 species, of which 21 are kinorhynchs representing 15 genera. Of the kinorhynch transcriptomes, 18 are new, inclusive four species representing aberrant genera (Cateria styx, Franciscideres kalenesos, Zelinkaderes brightae and Zelinkaderes yong). This dataset also includes representatives of priapulans, loriciferans, nematodes and nematomorphs as outgroups. Our aims are: 1) Building a robust phylogeny of Kinorhyncha utilizing hundreds of genes from transcriptomes. 2) Resolving the position of aberrant kinorhynchs, especially of Cateria styx. 3) Determining if distinct segmentation appeared before or after the diversification of kinorhynchs based on the position of the aberrant forms.

Taxon sampling
Eighteen kinorhynch species from fifteen genera ( Fig. 1) were sampled for de novo transcriptome sequencing and combined with three publicly available transcriptomes of Echinoderes kohni (Varney et al., 2019), Echinoderes dujardinii (Laumer et al., 2015), and Pycnophyes sp. (Smythe et al., 2019). Each genus is represented by a single species, except for Echinoderes, Pycnophyes and Zelinkaderes that are represented by two to four species each. Details on the selected species, sampling locality and SRA accession number are summarized in Table 1. Samples were collected during several campaigns from 2015 to 2019 in Canada, Brazil, Italy, Korea and USA. Necessary collection permits were obtained from each country through individual applications as in Brazil (Sistema de Autorizaçao e informação em biodiversidade, SISBIO-47601) and South Korea (Department of World Cultural Heritage − 12568 (2018.04.27)) or covered under general collection permits from host research institutions. Specimens were extracted from the sediment using the "bubble and blot" method (Sørensen and Pardos, 2020) and subsequently isolated, identified, washed in autoclaved seawater and stored in cryotubes with RNAlater (Thermo Fisher) at − 80 • C. As outgroups we used raw available transcripts from two priapulans (Meiopriapulus fijiensis and Priapulus caudatus), one loriciferan (Armorloricus elegans), one nematomorph (Nectonema munidae) and three nematodes (Oncholaimidae sp., Caenorhabditis elegans and Caenorhabditis remanei). Predicted transcripts from available genomes for C. elegans and C. remanei were obtained from Table 1 List of the species included in the transcriptomic analysis and corresponding SRA accession numbers. Boldfaced species are those whose transcriptomes were generated for the present study, and deposited under BioProject accession number PRJNA728538. Abbreviations: G, genomic; NA, not available; T, transcriptomic. EnsemblMetazoa (http://metazoa.ensembl.org/).

Transcriptome sequencing, assembly and processing
RNA was extracted from single individuals using Clonetech SMART-Seq HT kit (Clonetech) following manufacturer's instructions. Each specimen was recovered from RNAlater, rinsed several times in ddH 2 O, cut with a sterilized micro scalpel, and quickly transferred to a buffer for lysis, cDNA synthesis and PCR amplification (16 cycles). cDNA quality, concentration and molecular weight was assessed using an Agilent 2100 Bioanalyzer system (Agilent Biosciences) with the High sensitivity DNA kit (50-7000 bp). Dual index libraries were prepared using a Nextera XT DNA kit (Illumina) with approximately 1 ng of cDNA as input, multiplexed in groups of 10 and sequenced in several lanes of an Illumina NovaSeq 6000 system (10 libraries per lane, 150 bp paired-end). Sequencing and demultiplexing was carried out by Genewiz (Leipzig, Germany). Raw reads had on average 55.6 million reads per taxon ranging from 39 to 78.4 million (Table S1). They were deposited in the NCBI sequence read archive (SRA) with BioProject accession number PRJNA728538 (Table S1). Raw reads from published transcriptomes had on average 28.7 million reads (range 11.5-72.4 million) and were processed in the same manner as the new sequence data (Table S1).
Paired-end reads were quality assessed with FASTQC v.0.11.8 (www. bioinformatics.babraham.ac.uk) and subsequently trimmed using TrimmGalore! v.0.6.5 (www.bioinformatics.babraham.ac.uk). Parameters were adjusted to analyse paired-end libraries and automatic detection of adapter sequences, otherwise default settings were used. Except for the genome-derived transcripts of C. elegans and C. remanei, all transcriptomes were de novo assembled using the transcriptome assembly pipeline of Agalma v.3 (Dunn et al., 2013). The pipeline does an initial quality check using FASTQC, Bowtie2 mapping of insert size (Langmead and Salzberg, 2012), assembly using Trinity v.2.9.1 (Grabherr et al., 2011) and translation with blast hits against SwissProt database. Completeness of the transcriptome assemblies was assessed with BUSCO v.3 (Benchmarking Universal Single-Copy Orthologs, Simão et al. (2015)) ran on gVolante v1.2.1 (Nishimura et al., 2017) using the metazoan database with default settings. Raw read values, quality and assembly statistics are compiled in Table S1. The assembly pipeline was run on the server of the Biocomputing Core Facility, Department of Biology, University of Copenhagen.

Phylogenetic analyses
Orthologue identification was carried out on the translated amino acid sequences in OrthoFinder v. 2.5.2 (Emms and Kelly, 2019) using the -M msa option for gene tree inference and disabling alignment trimming (-z option). After orthogroup identification, multiple sequence alignments were generated for 115,378 orthogroups using MAFFT v.7.453 (Katoh et al., 2002). Initial gene trees were generated with FastTree v.2.1.11 (Price et al., 2010) within OrthoFinder to build multi-copy gene trees for 54,907 of the alignments with at least four sequences (the minimum needed to build a phylogenetic tree). These orthogroups often contain various duplications of genes in different species, while most phylogenetic analyses require one gene sequence per species. In order to obtain single-copy sequences, we used PhyloTreePruner (Kocot et al., 2013) to identify the most inclusive single-copy subtree from the rooted multi-copy gene trees produced by OrthoFinder. After pruning, 15,976 single-copy alignments with at least four sequences remained, which were realigned with MAFFT. Single-copy gene trees were estimated with IQ-TREE2 (Minh et al., 2020) using automatic model selection of Mod-elFinder (Kalyaanamoorthy et al., 2017) and 1000 ultrafast bootstrap replicates (Hoang et al., 2018).
We analysed 15,976 gene trees in a coalescent-based species tree framework using ASTRAL-III v.5.15.4 (Zhang et al., 2018). We also used pruned alignments that had a high representation from the 28 species to generate concatenated matrices. This included a 70% occupancy-matrix (at least 20 taxa present, 382 gene regions; 212,043 amino acids) and an 80% occupancy-matrix (at least 22 taxa present, 171 gene regions; 90,124 amino acids), which were analysed in IQ-TREE2 using Model-Finder, partition merging and 1000 ultrafast bootstraps.
We also explored a recently developed extension of ASTRAL that uses multiple copies for each species. The 54,907 multi-copy gene trees from OrthoFinder were used in coalescent-based species tree estimation using ASTRAL-Pro v.1.1.5 . Coalescent-based species tree estimation depends on the number of gene trees (Mirarab et al., 2016;Molloy and Warnow, 2017) and we therefore summarized all available gene trees. In addition, we tested the robustness of the obtained topology by reanalysing the multi-copy dataset with ASTRAL-Pro after (1) removing the kinorhynch lower quality transcriptomes from previous studies that had <50% BUSCO genes (E. dujardinii, E. kohni, Pycnophyes sp.), and (2) removing the transcriptomes of the Nematoida and Loricifera to assess if their long branches influenced the internal topology of kinorhynchs. Computational analysis was carried out on the Danish National Supercomputer for Life Sciences Computerome 2.0.
A simplified workflow of all the phylogenetic analyses carried out in this study is included in Supplementary figure 1.

Results
Our dataset includes 28 terminals of which 21 are kinorhynchs representing all families, 15 out of the 31 existing genera, and seven outgroup representatives of nematodes, nematomorphs, loriciferans and priapulans. Assembly statistics, quality assessment values and completeness of each transcriptome are summarized in Table S1. Newly generated kinorhynch transcriptomes have BUSCO completeness values ranging from 69 to 96% compared to 13-48% from the existing ones.
Kinorhyncha was monophyletic in all analyses with priapulans consistently recovered as their sister group with full support, and loriciferans appearing as sister group of kinorhynchs and priapulans (support values: bs = 100, PP = 0.98) (Fig. 2). Our results show well supported trees and topology congruency between concatenation with 70% occupancy (382 gene regions) and coalescent-based analyses (15.976 gene trees) with single copy gene trees (  figure 6) yielded an identical topology to that of the two first analyses (Fig. 2, Supplementary figure 2) with slightly higher support values. Here we favour the topology derived from the 70% concatenation matrix (Fig. 2) since it is based on twice as many genes and is congruent with the single copy coalescent analysis (Supplementary figure 2).
Our results show with maximum support that Kinorhyncha is divided into two major clades Cyclorhagida and Allomalorhagida (agreeing with Sørensen et al., 2015). The first clade, except for the exclusion of C. styx, is consistent with Cyclorhagida sensu Sørensen et al. (2015) and recovers Campyloderes vanhoeffeni as sister group of two subclades, Kentrorhagata comb. nov. and Echinorhagata. The support of the relationship between the two subclades is moderate in the analyses of the 70% occupancy matrix (bs = 69, Fig. 2) and when all the single-copy gene trees are analysed (PP = 0.89, Fig. 2, Supplementary figure 2), but it increases when analysing the multi-copy gene trees (PP = 0.99, Supplementary  figure 3). The cyclorhagid subclade, Echinorhagata, has maximum support across all analyses and includes Echinoderes species mixed with M. macracanthus, thus recovering Echinoderes as paraphyletic. Its sister clade, Kentrorhagata comb. nov., also receives full support in all analyses, and accommodates all the kentrorhagid species, except for C. styx, making Kentrorhagata sensu Sørensen et al. (2015) polyphyletic. Within this subclade there are two additional clades with full support in all analyses: one including Centroderes spinosus as sister taxon to Zelinkaderes brightae and Tubulideres seminoli, and another one including Zelinkaderes yong, Semnoderes armiger, Sphenoderes neptunus and Antygomonas paulae. The two aberrant species of Zelinkaderes, nested within the Kentrorhagata comb. nov., are not recovered as sister taxa in any analysis. Among the analyses there are discrepancies in the relative position of Z. yong and S. armiger (Fig. 2, Supplementary figure 2 Fig. 2), and Z. yong as the sister group to all the previous.
The second major kinorhynch clade, Allomalorhagida, is composed of Cateria styx, Dracoderes abei, Paracentrophyes quadridentatus, Franciscideres kalenesos, Cristaphyes yushini, Pycnophyes ilyocryptus, Pycnophyes giganteus and Pycnophyes sp. Except for the inclusion of C. styx, for which no molecular data existed previously, all remaining lineages correspond with Allomalorhagida sensu Sørensen et al. (2015). Within Allomalorhagida there are two subclades, Pycnophyidae and Anomoirhaga nom. nov., with maximum support (Fig. 2). Pycnophyidae accommodates Pycnophyes spp. and C. yushini as sister groups, which supports the monophyly of the family Pycnophyidae sensu Sánchez et al. (2016). The second clade nests a mix of aberrant and non-aberrant kinorhynchs with C. styx branching out as sister group to the remaining taxa and D. abei branching out next, as sister group to F. kalenesos and P. quadridentatus. The two aberrant genera Cateria and Franciscideres are thereby not closest relatives within the clade. Anomoirhaga nom. nov. accommodates four genera represented in the present analyses, Cateria, Franciscideres, Paracentrophyes and Dracoderes (Fig. 2). We suggest that also Gracilideres, Mixtophyes and Neocentrophyes are assigned to Anomoirhaga nom. nov. (for further discussion see Section 4.2. below).
Based on the obtained results we propose a new kinorhynch classification (Table 2). Additional genera not included in the present dataset were classified according to Sørensen et al. (2015), which we find to be justified by the generally great congruence between the results of the two studies as well as the morphological considerations mentioned in Sørensen et al. (2015).

Discussion
Our phylogenomic analysis based on hundreds to thousands of gene regions provided a well-resolved, stable phylogenetic hypothesis for Kinorhyncha, and a necessary backbone to understand their evolution. The results largely overlap with the topology shown in Sørensen et al. (2015). The main differences are: i) Cateria styx is no longer a cyclorhagid, but is instead recovered within Allomalorhagida, nested in the new clade Anomoirhaga nom. nov.; ii) Due to the new position of C. styx, Kentrorhagata sensu Sørensen et al. (2015) is polyphyletic and is here redefined as Kentrorhagata comb. nov.; iii) Dracoderes is no longer recovered as sister group of all other allomalorhagids but nested within Anomoirhaga nom. nov.; iv) Zelinkaderes is not recovered as monophyletic.
Relevant congruencies between both studies are: i) Echinoderes is recovered as paraphyletic; ii) Campyloderes is confirmed as sister group of all other cyclorhagids; iii) Echinoderidae sensu  and Pycnophyidae sensu Sánchez et al. (2016) remain monophyletic.

Cateria is an allomalorhagid
Since the description of the first Cateria species, more than 60 years ago by Gerlach (1956), the genus has been quite enigmatic due to several morphological characteristics (e.g., very elongated habitus, thin cuticle, indistinct neck, elongated scalids) that led to consider the genus as aberrant (Higgins, 1968). Cateria was originally assigned to Cyclorhagida based solely on morphological characters (Higgins, 1968). However, the description of the new aberrant genera Franciscideres and Gracilideres, and their unequivocal assignment to Allomalorhagida (Family Franciscideridae) based on molecular data (Dal Zotto et al., 2013;Yamasaki, 2013 (note that in the latter Gracilideres is named as 'undescribed' or 'new genus')), prompted further investigations into Cateria's phylogenetic position. The first kinorhynch phylogeny including molecular as well as morphological data , positioned Cateria within Cyclorhagida only based on morphological data. However, Sørensen et al. (2015) discussed the contentious position of Cateria, suggesting its potential affinity with Franciscideridae and the need of obtaining molecular data for the genus. Our phylogenomic study clearly supports Cateria as an allomalorhagid taxon, and recovers C. styx as sister group of the clade composed of D. abei, P. quadridentatus and F. kalenesos. Due to the congruency and robustness of our results, also confirming the suspicions of Dal Zotto et al. (2013), Yamasaki et al. (2013) and Sørensen et al. (2015), we find it justified to now consider Cateria as an allomalorhagid genus nested within Anomoirhaga nom. nov.
Recent morphological data likewise supports a closer relationship between Cateria and Franciscideres, rather than between Cateria and the cyclorhagid, aberrant genus Zelinkaderes (Herranz et al., 2021a,b). This also agrees with the minor neuro-and myoanatomical differences described between Cateria and Franciscideres, compared with those in Zelinkaderes (Herranz et al., 2021a,b), indicating an independent origin of their aberrant habitus.
From a morphological point of view, it is difficult to find similarities among such different genera. Obvious resemblances are related to the aberrant appearance in C. styx and F. kalenesos, thoroughly discussed in previous studies (Dal Zotto et al., 2013;Sørensen et al., 2015;Herranz et al., 2019bHerranz et al., , 2021a. However, within Anomoirhaga nom. nov., F. kalenesos and P. quadridentatus appear as sister taxa, and as most closely related with D. abei, suggesting the aberrant Cateria and Franciscideres as distantly related within the clade. There are no unambiguous morphological synapomorphies that support all anomoirhagid taxa. However, some characters might be of systematic relevance, such as the presence of incomplete divisions on segment 1, presence of alternatingly displaced middorsal spines, and their thin cuticle. Partial episternal (ventrolateral) divisions at the anterior margin of segment 1 were originally described as a genus diagnostic character for Paracentrophyes species (Higgins, 1983) but recently, partially differentiated ventrolateral divisions were also reported from Dracoderes nidhug (Thomsen et al., 2013). Also, species of Cateria show ventrolateral divisions (Higgins, 1968), which extend from anterior to posterior margins in C. gerlachi (Neuhaus and Kegel, 2015) while they are present only in the posterior 2/3 of the segment in C. styx (Herranz et al., 2019b). F. kalenesos and G. mawatarii have no indication of longitudinal divisions on segment 1 (Dal Zotto et al., 2013;Yamasaki, 2019;Rucci et al., 2020), but the morphological similarities between these two genera are so evident Yamasaki, 2019) that it makes sense to consider Gracilideres (not included in the present analysis) as part of Anomoirhaga nom. nov., and probably as sister group of Franciscideres. Two additional allomalorhagid genera not included in this analysis are Mixtophyes and Neocentrophyes. Both have been suggested as close relatives to Paracentrophyes based on morphology , thus for now, we tentatively consider both genera as part of Anomoirhaga nom. nov. (Table 2). Neither Mixtophyes nor Neocentrophyes show any full or partial differentiation of the sternal plate of segment 1 (Higgins, 1969;Sánchez et al., 2014). This makes it difficult to predict whether differentiation of this sternal plate is autapomorphic for Anomoirhaga nom. nov., and subsequently lost within the clade, or if the subdivision of the sternal plate has occurred independentlymaybe even multiple timeswithin Anomoirhaga nom. nov. It is in any case worth noticing that most taxa of the clade show a tendency towards a ventral plate differentiation of segment 1.
The alternating position of dorsal spines has been one of the key characters to identify Dracoderes species (Higgins and Shirayama, 1990;Sørensen et al., 2012). However, recent studies also found slight alternation of the middorsal spines in F. kalenesos, G. mawatarii and C. styx (Yamasaki, 2019;Herranz et al., 2019b;Rucci et al., 2020). This character has so far been reported from these four genera only, since species of Paracentrophyes, Mixtophyes and Neocentrophyes have middorsal cuticular processes instead of middorsal spines (Higgins, 1969(Higgins, , 1983Sánchez et al., 2014). Thus, the alternation of dorsal spines could potentially be an autapomorphy for Anomoirhaga nom. nov.
Another character shared by most anomoirhagid species is the presence of thin cuticle in adults. Although cuticle thickness can be a very subjective character and measurements usually are unavailable, it is indisputable that adults of Franciscideres, Cateria and Gracilideres have conspicuously thin cuticle. Likewise, the cuticle of Paracentrophyes, Mixtophyes and Neocentrophyes is thinner than in the heavily armoured pycnophyids, sister group to Anomoirhaga nom. nov. Conversely, species of Dracoderes usually show thick, rigid cuticle. It is also noteworthy that a weak cuticularisation is not exclusive to the lineages contained in Anomoirhaga nom. nov., but also present in some cyclorhagid lineages such as Zelinkaderes and Triodontoderes.

Relationships within Kentrorhagata comb. nov.
All our analyses support the monophyly of Kentrorhagata comb. nov. However, morphologically, the topology within the clade is not completely meaningful. Our results recovered the two Zelinkaderes species (Z. brightae and Z. yong) as part of two different clades, where Z. brightae is sister taxon of T. seminoli, and Z. yong is sister to a clade composed of Semnoderes, Antygomonas and Sphenoderes species (Fig. 2). We do not see this as an indication of potential paraphyly of Zelinkaderes  Sørensen and Thormar, 2010Xenosomata Zelinka, 1907Campyloderidae Remane, 1929Campyloderes Zelinka, 1907Ryuguderes Yamasaki, 2016 though. The genus is morphologically well-supported, and it accommodates very similar species that mostly differ in spine and tube patterns (see, e.g., Higgins, 1990;Altenburger et al., 2015;Herranz et al., 2021b). Thus, we see the apparent paraphyly of Zelinkaderes as result of the limited taxon sampling of kentrorhagid terminals.
All analyses support Antygomonas paulae as sister taxon to Sphenoderes neptunus, and the topology derived from the 70% concatenation matrix, the ASTRAL single-copy species tree, and the multi-copy species tree without long-branched outgroups (Fig. 2, Supplementary figure 2 and Supplementary figure 6) support Semnoderes armiger as their sister taxon. The conspicuous morphological similarities and potential synapomorphies shared between these three genera have already been pointed out in several studies (Sørensen et al., 2010;Herranz et al., 2014;Sørensen and Landers, 2018). Thus, based on their potentially synapomorphic modifications in segment 1, closing mechanisms of neck and anterior trunk segment (see papers cited above), and the obtained tree topology (Fig. 2), we find it justified to reassign Antygomonas to Semnoderidae.

Midterminal spine as a cyclorhagid character only?
The presence of a midterminal spine in adults have so far been considered an exclusively cyclorhagid character present in all Kentrorhagata and Xenosomata but absent in Echinorhagata . The new position of C. styx, previously considered a cyclorhagid and having a midterminal spine, within Allomalorhagida contradicts that this character is a cyclorhagid autapomorphy. However, while the motile midterminal spine of the kentrorhagids has a pair of associated longitudinal muscles (Müller and Schmidt-Rhaesa, 2003) and a serotonin-like immunoreactive loop (Herranz et al., 2013), these structures are absent in C. styx (Herranz et al., 2021a,b). Thus, the midterminal spine of C. styx shows closer morphological resemblance to non-motile middorsal spines, suggesting it to be a middorsal spine displaced to a more terminal position. This is supported by the fact that C. styx lacks a middorsal spine on segment 11, opposite to most kentrorhagids, which show both midterminal and middorsal spines on this segment. This evidence all together suggests that the midterminal spine of C. styx is not homologous with the midterminal spine of kentrorhagids.
Other allomalorhagids with midterminal structures include Paracentrophyes species that show a midterminal process in juveniles and adults (Neuhaus, 1995). However, this structure is a cuticular, nonarticulated elongation of the trunk cuticle, and therefore considerably different from the articulated, motile kentrorhagid midterminal spine. Thus, the motile midterminal spines observed in cyclorhagids appear to be unique. Myo-and neuroanatomical information is still unavailable for species of Xenosomata, but such data could clarify whether the cyclorhagid midterminal spines evolve twice, at the branches leading to Xenosomata and Kentrorhagata comb. nov., respectively, or if midterminal spines are autapomorphic for all cyclorhagids and secondarily lost in Echinorhagata.

Evolution of segmentation in kinorhynchs
According to the most parsimonious interpretation of our results, segmentation is a synapomorphic character for the crown group kinorhynchs, whereas the aberrant worm-like body plans evolved at least three times independently within the phylum: at least once in cyclorhagids (Zelinkaderes) and twice in allomalorhagids (Cateria and Franciscideres). Additionally, morphological studies have shown that despite the weak external segmentation in the aberrant forms, both musculature and nervous system follow a segmental-like pattern, although not always showing a one-to-one correspondence with the cuticular divisions (Herranz et al., 2021a,b). If we hereby assume that segmentation evolved progressively among stem group kinorhynchs, we should search for traces of this transition among extinct lineages. Although their kinorhynch affinity is still questioned, the only potential stem kinorhynchs, Eokinorhynchus rarus and Zhongpingscolex qinensis (Zhang et al., 2015;Shao et al., 2020), show a vermiform habitus with annuli and no clear indications of segmentation. Only E. rarus might show indications of plates within annuli (Shao et al., 2020). More promising, yet undescribed, kinorhynch-like fossils from the Qingjiang biota (Fu et al., 2019) seem to show more resemblance to extant kinorhynchs, with signs of trunk segments, supporting the appearance of segmentation before kinorhynch diversification.

Conclusions
We present a new kinorhynch phylogeny based for the first time on a transcriptomic dataset composed of 21 in-group species (15 genera) and seven outgroup species (four phyla). Our results show a well-supported and robust topology where aberrant forms, represented by Cateria, Franciscideres and Zelinkaderes species, are congruently located in distant parts of the tree, indicating that their worm-like appearances evolved convergently. This suggests that segmentation is synapomorphic for kinorhynchs and probably evolved along the kinorhynch stem lineage.
The phylogeny resulted in a new kinorhynch classification where Cyclorhagida is composed of the clades Echinorhagata, Kentrorhagata comb. nov., and Xenosomata; and Allomalorhagida is composed of the clades Pycnophyidae and Anomoirhaga nom. nov. Pycnophyidae accommodates Pycnophyes and all genera described by Sánchez et al. (2016). Anomoirhaga nom. nov. accommodates very morphologically different genera such as Dracoderes, Franciscideres, Paracentrophyes, Cateria (the latter previously nested within Cyclorhagida), and tentatively also Gracilideres, Mixtophyes and Neocentrophyes. We encourage that detailed morphological investigations are conducted in the future to identify autapomorphies for the group.

Data availability statement
Raw data generated for this article is available on the NCBI Sequence Read Archive (SRA) with BioProject accession number PRJNA728538. Input and output files of all phylogenetic analyses are available on

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.
Center for Marine Studies, Federal University of Paraná (CEM/UFPR) for offering lab space and student support with field collections in Pontal du Sul in 2015 and 2019; and the staff at the Smithsonian Marine Station at Fort Pierce, Florida for their support with field collections in 2016. Kevin M. Kocot and Rebecca M. Varney are acknowledged for suggesting protocols and guidance. Matteo Dal Zotto is acknowledged for providing the SEM image of Pycnophyes giganteus used in Fig. 2. We thank the Faculty of SCIENCE and the Department of Biology at the University of Copenhagen for computational resources at Computerome2 and Pallas.

Funding
This project received funding from the European Union's Horizon 2020 research and innovation programme, under the Marie Skłodowska-Curie grant agreement No 797140 to MH. Sampling in South Korea and Brazil was funded by the Carlsberg Foundation to MVS (CF17-0054 and CF2013_01_0035).