Dicyemida and Orthonectida: Two Stories of Body Plan Simplification

Abstract

Two enigmatic groups of morphologically simple parasites of invertebrates, the Dicyemida (syn. Rhombozoa) and the Orthonectida, since the 19th century have been usually considered as two classes of the phylum Mesozoa. Early molecular evidence suggested their relationship within the Spiralia (=Lophotrochozoa), however, high rates of dicyemid and orthonectid sequence evolution led to contradicting phylogeny reconstructions. Genomic data for orthonectids revealed that they are highly simplified spiralians and possess a reduced set of genes involved in metazoan development and body patterning. Acquiring genomic data for dicyemids, however, remains a challenge due to complex genome rearrangements including chromatin diminution and generation of extrachromosomal circular DNAs, which are reported to occur during the development of somatic cells. We performed genomic sequencing of one species of Dicyema, and obtained transcriptomic data for two Dicyema spp. Homeodomain (homeobox) transcription factors, G-protein-coupled receptors, and many other protein families have undergone a massive reduction in dicyemids compared to other animals. There is also apparent reduction of the bilaterian gene complements encoding components of the neuromuscular systems. We constructed and analyzed a large dataset of predicted orthologous proteins from three species of Dicyema and a set of spiralian animals including the newly sequenced genome of the orthonectid Intoshia linei. Bayesian analyses recovered the orthonectid lineage within the Annelida. In contrast, dicyemids form a separate clade with weak affinity to the Rouphozoa (Platyhelminthes plus Gastrotricha) or (Entoprocta plus Cycliophora) suggesting that the historically proposed Mesozoa is a polyphyletic taxon. Thus, dramatic simplification of body plans in dicyemids and orthonectids, as well as their intricate life cycles that combine metagenesis and heterogony, evolved independently in these two lineages.

Introduction

In spite of more than one hundred years of studies, the evolutionary relationships of the Mesozoa are still elusive. The name of this taxon reflects the traditional view of mesozoans as organisms with intermediate organization between unicellular protozoans and multicellular metazoans (Van Beneden, 1876; Hyman, 1940). Indeed, the two groups of microscopic parasitic invertebrates, the Dicyemida, and Orthonectida, display a remarkably simple morphological organization and a nearly complete absence of tissues and organs (Malakhov, 1990). Adult dicyemids inhabit the renal sacs of cephalopod mollusks and consist of just about 40 somatic cells, lack recognized muscular, nervous, sensory cells, and the organs typical for eumetazoans (Furuya et al., 2004). Dicyemids do not have a morphologically recognized basal membrane (Czaker, 2000), and never develop “true” tissues throughout their complex life cycle (Furuya and Tsuneki, 2003). The trophic stage of orthonectids is a syncytial plasmodium, which resides inside the invertebrate host and generates ephemeral ciliated organisms that exit the host for reproduction (Slyusarev, 2008). These organisms are composed of several hundred somatic cells without anatomically recognized digestive, circulatory, or excretory systems. Before the discovery of muscular and nervous systems in the swimming stages of orthonectids (Slyusarev and Starunov, 2015), they were thought to have a planula-like organization and were grouped with dicyemids in the Mesozoa as multicellular animals with an incredibly simple body plan, perhaps – the simplest among all Metazoa, and comparable to placozoans.

Intricate life cycles of dicyemids and orthonectids exhibit the alternation of asexual and sexual generations, termed metagenesis. Ameiotic generative cells (agametes) develop inside the dicyemid axial cell and later produce the next vermiform generation possessing gametic cells that undergo self-fertilization. In orthonectids, agametes develop inside the parasitic plasmodium and produce the free-living diecious (or hermaphroditic) generation (Cheng, 1986; Slyusarev, 2008). The phenomenon of successive sexual parthenogenetic and amphimictic generations is termed heterogony. In this sense, orthonectids and dicyemids as well as parasitic flatworms combine metagenesis and heterogony in their life cycles. Particularly, trematode sporocysts and rediae that parasitize gastropod mollusks produce the next generation from ameiotic generative cells (Dobrovolskij and Ataev, 2003; Ataev, 2017). Similarities in life cycles for long sustained the hypothesis about close relationships of dicyemids and orthonectids with digenetic trematodes. On the other hand, intracellular localization of generative cells relates dicyemids and orthonectids with myxozoans rather than trematodes. Such intricate combination of traits makes life strategies in dicyemids and orthonectids unique among animals.

The phylogenetic affinity of dicyemids and orthonectids has been called into question on the grounds of morphology (Kozloff, 1990; Brusca and Brusca, 2003; Ruppert et al., 2004). Molecular data conclusively demonstrated that both dicyemids and orthonectids are in fact bilaterians (Katayama et al., 1995; Hanelt et al., 1996; Pawlowski et al., 1996; Aruga et al., 2007) and belong to the diverse clade of Lophotrochozoa (=Spiralia) (Kobayashi et al., 1999, 2009; Petrov et al., 2010; Suzuki et al., 2010; Mikhailov et al., 2016; Lu et al., 2017; Schiffer et al., 2018), thus implying that their simple organization evolved as the result of their parasitic lifestyle.

In molecular phylogenetic analyses, dicyemid and orthonectid lineages display extremely high levels of divergence, and their exact placement among the spiralians remains ambiguous and potentially prone to long branch attraction artifacts. Complicating the matter is the uncertainty in relationships between other spiralian taxa, including the Annelida, Mollusca, Nemertea, Brachiopoda, Entoprocta, and Bryozoa (Kocot, 2016). Recent phylogenomic analyses lead to conflicting conclusions regarding the mesozoan phylogeny. Lu et al. (2017) using a dataset of 348 orthologs (58,124 alignment positions) from 23 spiralian species, including an orthonectid and a dicyemid, report the monophyly of the Mesozoa either as a sister group to the Rouphozoa (Platyhelminthes + Gastrotricha) or within the Gastrotricha. Alternatively, Schiffer et al. (2018) using a dataset of 469 orthologs (190,027 alignment positions) from 29 spiralian species, including an orthonectid and two dicyemids, conclude that Orthonectida and Dicyemida evolved independently within the Lophotrochozoa, with the orthonectids exhibiting clear affinity to annelids, and dicyemids occupying an isolated position within Lophotrochozoa. Here, we obtained transcriptomic and genomic data for dicyemid species to resolve this contradiction.

The dicyemid genome is distinguished by uncommon features, such as the genome rearrangements during the life cycle and generation of circular DNAs (Noto et al., 2003), including those that encode mitochondrial proteins and rRNAs (Watanabe et al., 1999; Catalano et al., 2015). It is not yet established if the mitochondrial protein-coding genes are encoded only by small circular DNA molecules (Watanabe et al., 1999) or whether they are produced during the dicyemid development from a precursor mitochondrial DNA with a more typical metazoan organization (Awata et al., 2006). Using high-throughput genomic sequencing we sought to find any properties of dicyemid sequences that would reveal their genome organization. We also estimated the extent of gene losses due to the simplification of dicyemid morphological organization, and analyzed whether losses in particular gene families and regulatory pathways are the same or different compared to an orthonectid Intoshia linei.

Results and Discussion

Genomic Sequencing and Assembly of Dicyema sp.

Direct assembly of a dicyemid genome from whole DNA extracts using standard approaches is an extremelly challenging problem due to drastic genome rearrangements that occur in dicyemids during development. Previous studies have demonstrated that somatic cells of dicyemids undergo drastic genome rearrangements and chromatin elimination (Noto et al., 2003), and suggested that selective and whole genome amplification takes place at different stages of their development (Awata et al., 2006). Accordingly, the sequencing of whole DNA extracts from Dicyema sp. resulted in a highly fragmented assembly with uneven coverage and N50 of 942 bp, where the largest contig was only around 20 Kb. The total size of the assembly is 858 Mbp in nearly 1 million contigs over the length of 500 bp, and includes contaminating cephalopod sequences. Due to significant genetic difference between the dicyemid host Enteroctopus dofleini and the available genomic sequence of Octopus bimaculoides, the filtering of the assembly was performed at the level of predicted gene products. Only predictions identifiable by hits against the InterPro database were retained for the subsequent comparative analyses and filtered from the cephalopod contamination using the best hit approach with BLAST searches against the NCBI nr database. Out of 38,410 predictions with InterPro hits, 43% were discarded as contamination, resulting in 21,842 putative dicyemid genes with 71% complete and 12% fragmented universal eukaryotic orthologs evaluated by BUSCO (Table 1). Similar values are obtained for gene predictions after normalizing on the number of BUSCOs found in at least one filtered transcriptome: 76% complete and 12% fragmented. The total percentage of BUSCOs recovered by at least one sequencing library, including genomic and transcriptomic filtered data, approaches values seen in typical metazoan genomes: 91% complete and 3% fragmented. For all analyses in Sections 2.4–2.11 we used original genomic data on Dicyema sp., and the three transcriptomes, including the two originally obtained and the one of Dicyema japonicum available from the published source (Lu et al., 2017).

The dicyemid genes display miniaturization of spliceosomal introns – the median length of introns is 27 bp, and approximately two thirds of predicted introns are under the length of 30 bp (Figure 1). This agrees with an earlier survey that revealed extreme intron shortening in a set of 40 genes from D. japonicum (Ogino et al., 2010). The estimated intron density in Dicyema sp. is 4.9 introns/gene for predictions with intact start and stop codons, which is also similar to the 5.3 introns/gene reported for D. japonicum. Similar value of intron density is seen in the genome of orthonectid I. linei (Mikhailov et al., 2016). Notably, the orthonectid genes also harbor short spliceosomal introns, but the majority of its introns are longer than 30 bp, and the median size is 57 bp, considerably exceeding the intron lengths observed in dicyemid genes.

FIGURE 1

“Circular” Contigs in Genomic Assembly of Dicyema sp.

Using the genomic assembly we have identified 24,065 “circular” contigs (see section “Materials and Methods”). The distribution of circular contig lengths in the assembly is multimodal (Figure 2A). The first abundant pool of sequences is formed from contigs less than 500 bp. The second pool, which includes sequences of a length over 500 bp, consists of 3,220 contigs with the median length of 702 bp. The properties of the sequences in this pool (such as length and abundance) are consistent with previous data of DNA gel electrophoresis, EM and PCR experiments (Noto et al., 2003), which supports the conjecture that these sequences are circular DNA rather than direct repeats. “Short” circles (up to 500 bp length) were shown to possess 38.1% low complexity regions, while “long” circles – only 2.9%. This observation might suggest that a fraction of predicted short circles represents direct repeats. Following this rationale, we considered the two sub-pools separately in analyses.

FIGURE 2

The lengths of sequences from the second pool of circular contigs are distributed non-uniformly which is particularly evident within the 600–800 bp range (Figure 2B). The average distance between two adjacent peaks of this distribution is 10.44, which closely corresponds to the number of base pairs in one turn of B-DNA. Multimodal distribution was also observed (Kolmogorov–Smirnov test p-value is 0.999) when performing assembly with the varying k-mer size (55 or 77) and with another assembly method (Supplementary Figure S1). The presence of this pattern is unexpected, and presumably could be attributed to the greater stability of circles or tendency to circularize for molecules with an integer amount of turns of a relaxed form of DNA. A similar effect has also been observed in short (<200 bp) sequences as a result of rolling circle replication bias (Joffroy et al., 2018). This distribution can result from the random ligation of linear molecules cut from the genome as it leads to the reduction in DNA supercoiling. Alternatively, replicating mini-circular DNA molecules can be selected in length to reduce their supercoiling. Figure 3 shows that the coverage value for “long” circular contigs is not lower than for linear ones, which casts doubt on the proposed diminution of circular molecules during ontogenesis.

FIGURE 3

Long circular contigs are predominantly not similar in nucleotide sequences. Only 15% of them have at least one fairly similar contig, and only three families of contigs unite more than 10 members (Figure 4).

FIGURE 4

Two independent motif detection methods (Bailey and Elkan, 1994; Rubanov et al., 2016) have been applied to the circles of length 600–800 bp with a coverage logarithm of over 3 (2,031 sequences). In 1,871 sequences (92.12% of sequences in the analysis) common motifs have been found (E-value: 4.8e-82, see Figure 5A).

FIGURE 5

At a p-value < 10⁻⁵, the most common motif occurs on average once every 874 bp in “long” circles and every 21,863 bp throughout the entire assembly (statistical significance of the difference provided by the chi-squared criterion: p-value < 0.001). The search for highly conserved sequences in various subsets of genome sequences has demonstrated that less common motifs with high information content can also be found in circles (Figure 5B–E).

The search for conserved domains in circular contigs recovered only domains of mtDNA-encoded proteins (10 conserved domains, 13 contigs including paralogs). These sequences are presumably transcribed as they are also found in the RNA-seq data (blastn search, E-value < 1e-30).

Mitochondrial DNA of Dicyema sp.

Genomic data confirm the localization of mitochondrial genes of dicyemids on circular DNA molecules (Watanabe et al., 1999; Catalano et al., 2015). The search for mtDNA genes in the genomic data found 21 circles with length varying from 344 to 1605 bp. The following gene sequences were found: cox1-3, cob, nad1-5, atp6, rrnL, rrnS, trnH, trnI, trnK, trnL1, trnN, trnP, trnQ, trnR, trnS2, and trnY (Figure 6). In earlier studies the dicyemid mitochondrial contigs were found to carry either one protein coding gene (Watanabe et al., 1999) or a protein coding gene and a tRNA gene (Robertson et al., 2018). We found one circle that contains two genes – cox2 and rrnS, and three circles that contain two tRNA genes each. Protein identity between mitochondrial predictions for Dicyema sp. and the earlier published D. japonicum (Robertson et al., 2018) varies from 39% (nad2) to 75% (cox1). The majority of mitochondrial genes can also be found in the transcriptomic data, except for atp6 and nad5. The mtDNA circles also contain the motif described above (Figure 5A) (p-value < 10⁻⁵).

FIGURE 6

The nad2 and atp6 genes were found in two different variants in the genomic data. Two paralogs of nad2 with lengths of 215 and 252 amino acids have 42% identity at the amino acid level. Two paralogs of atp6 with lengths of 117 and 149 amino acids have 89% identity at the amino acid level, and share two long deletions with other dicyemids. These deletions are specific for dicyemids and are not found in other taxa including Orthonectida. Both of dicyemid deletions are located outside of the transmembrane helices – the first one with the length of 16 amino acids is located in the region facing the mitochondrial matrix and the second one with the length of 17 amino acids is located in the region facing the intermembrane space, according to the alignment of atp6.

We predicted 11 mitochondrial tRNA genes in Dicyema sp. including two paralogs of glutamine tRNA gene (Supplementary Figure S2). Both dicyemid glutamine tRNAs have similar secondary structures and lack a T-arm. Dicyemid arginine tRNA also lacks a T-arm and lysine tRNA lacks a D-arm. Other mitochondrial tRNAs maintain the typical clover leaf structure, although several tRNA genes have single nucleotide insertions and/or non-complementary pairs in stems. Experimental evidence is needed to confirm all the predicted tRNA genes, as well as decisions whether numerous not listed tRNA-like sequences with p-value below the threshold are functional genes.

Read mapping to the genomic assembly revealed no read pairs that would facilitate mtDNA scaffolding. Whenever one read from a pair would map to the circular mitochondrial contig, the other would map to the same contig or have a sequence of low complexity. Thus, our genomic data fails to confirm the hypothesized existence of an unprocessed mtDNA precursor, which would generate the mtDNA circles (Awata et al., 2006).

The presence of common sequence motifs in circles with mtDNA genes and without them seems to be surprising. It can be interpreted as a consequence of a similar mechanism of generation and maintaining of circles irrespective of their function.

The partitioning of mtDNA into circular molecules is a rare feature for the animal mitochondrial genomes (Odintsova and Yurina, 2005; Burger et al., 2012; Kolesnikov and Gerasimov, 2012; Smith and Keeling, 2015; Lavrov and Pett, 2016, for review). In bilaterians, the mtDNA is fragmented into a large number of mini-chromosomes in the cyst-forming nematodes Globodera spp. (Armstrong et al., 2000; Gibson et al., 2007) and sucking lice (Shao et al., 2009). Notably, the mitochondrial DNAs from orthonectids Intoshia linei, Intoshia variabili, and Rhopalura ophiocomae retain typical structure for metazoans and encode the full set of mitochondrial genes on a single circular molecule (Robertson et al., 2018; Bondarenko et al., 2019). The reason why the mitochondrial Dicyema spp. genome is fragmented is unknown. Earlier, the fragmentation of the mitochondrial genome of sucking lice was considered (Shao et al., 2009) as an adaptation to the high rate of molecular evolution, which is even more characteristic of Dicyema spp. It is possible that under conditions of high mutagenesis, a set of uncorrupted genes is easier to assemble from individual than concatenated molecules.

Analysis of mitochondrial DNA suggests an explanation of the multiple observed circular contigs. For searches with the tblastx algorithm we used proteins from the annotated mitochondrial contigs as the query and all 3,220 “long” circle contigs as the database. The searches returned many circular contigs that encode highly diverged genes cox3, nad2, and nad4 (Figure 7). The cox3 homolog is largely diverged, while nad2 and nad4 contain stop codons and frame shifts. These contigs therefore represent mitochondrial pseudogenes.

FIGURE 7

Kiełbasa et al., 2011) (threshold E-value for local alignment 0.008).

Previous publications and our new data confirm the presence of two unusual features of the Dicyema genome. First, mitochondrial genes in Dicyema are not located on a single long DNA molecule as in most animals, but are partitioned into smaller circular molecules. The second interesting feature of Dicyema is the presence of thousands of non-coding circular DNA sequences. Both types of circular DNA molecules fall in a similar range of size and coverage in DNA assembly and bear a common set of similar 12–20 bp DNA patterns, which might be hypothetical signal sequences. We assume that all circular DNAs in Dicyema may have a common origin, although experimental evidence is necessary. We speculate that the presence of multiple mtDNA mini-rings instead of one long molecule might have produced serious problems in mitochondrial division. This requires special mechanisms to correct distribution of multiple minicircular DNA molecules upon mitochondrion division so that both descendants would obtain a complete set of genes. Specific signal patterns like the ones we observe could be used to support circular mtDNA duplication, their protection against elimination or their correct distribution between descendent mitochondria. When such mechanisms are established, it is possible that rings carrying mutated (pseudo)genes or other selfish non-coding DNA circular elements acquire similar signal sequences that ensure their preservation in a similar way as with parasitic mobile genetic elements.

Homeobox Transcription Factors

Homeodomain (homeobox) transcription factors are crucial regulators of animal development that play central roles in tissue differentiation and axial body patterning. Bilaterian genomes encode from over 300 to around 60 homeobox genes. The genome of orthonectid I. linei was found to possess one of the smallest repertoires of homeoboxes (Mikhailov et al., 2016), which matches the reduced complexity of their organization. To determine how the extreme simplification of body plan seen in dicyemids relates to their homeobox gene content we searched for these genes in the dicyemid genomic and transcriptomic data. For analyses with HMMER, we used gene predictions coming from the genomic assembly of Dicyema sp. (PRJNA527259; designated as Dicyema sp. 1) and transcriptome assemblies of Dicyema sp. (SRR827581; designated as Dicyema sp. 2) and Dicyema japonicum (DRR057371). HMMER searches using the homeodomain profile identified 38, 39, and 55 homeoboxes in the three dicyemids after filtering out contaminating cephalopod sequences. Phylogenetic inference suggests that dicyemid homeoboxes form up to 39 families, and each dicyemid was found to contain 31 or 34 families (Table 2). A high level of sequence divergence complicates classification of the dicyemid homeoboxes. Although most dicyemid sequences could be assigned to one of the homeobox classes (Zhong and Holland, 2011), their attribution to known families is inconclusive.

Phylogenetic inference with PRD class homeoboxes reveals six dicyemid families, including the previously identified orthologs of the Pax6 and Otx (Aruga et al., 2007; Kobayashi et al., 2009). Five dicyemid sequence groups were found among the LIM homeoboxes, two of which are grouped with the Lhx6/8 and Islet family sequences. An additional LIM class homeobox of the Lhx2/9 family was found in D. japonicum, but could neither be confirmed by data from the other dicyemids nor discarded as contamination. Another five dicyemid gene groups belong to the POU class homeoboxes, but branch outside of any known families. The dicyemids form at least 5 TALE class sequence groups, with one group branching within the Pbx family. An additional single member of the TALE Tgif family was found only in the genomic data. Four SINE class families were found among the dicyemid sequences, with one grouping with the Six3/6 family. The dicyemids also possess a zinc finger homeobox, and a group of Onecut family sequences, which can be subdivided into 2 dicyemid-specific families.

Reconstructions with the ANTP class homeoboxes recover 8 dicyemid sequence groups (Figure 8). Three of these groups fall within the central Hox sequences. One of the dicyemid central Hox groups corresponds to orthologs of DoxC – a dicyemid member of the spiralian Lox5 family (Kobayashi et al., 1999), which was shown to have an expression pattern consistent with defining anterior–posterior boundaries in the developing dicyemids (Aruga et al., 2007). The analysis suggests that dicyemids possess another member of the Lox5 family – all three dicyemids were found to encode a paralog of the DoxC. The paralog displays greater sequence divergence, but similar to other members of the family retains a Lox5-specific motif flanking the C-terminus of the homeodomain (de Rosa et al., 1999). The third dicyemid group within the central Hox sequences is found outside the Lox5 family and tends to group with the Lox2/Lox4 families, but beyond that does not lend itself to classification. A single posterior type Hox gene was found in D. japonicum, but once again could not be verified using other dicyemid sequences or rejected as contamination by BLAST searches. The dicyemid ANTP class homeoboxes also include members of the Evx, Dlx, and Nk2 families, and a conspicuous group of Hox-like sequences (Figure 8). Sequences within the dicyemid Hox-like assemblage share a common ancestor and retain a YPWM motif, which is essential for binding Hox cofactors (Prince et al., 2008), but this group is too divergent to be classified with any family, and is placed with the longest branch of Hox-like genes – the ParaHox Cdx family.

FIGURE 8

The survey of dicyemid genes suggests that overall they possess fewer homeoboxes than the orthonectid I. linei and their sequences are also markedly more diverged. Unlike the orthonectid, no ParaHox or anterior Hox families could be readily identified in the dicyemid data. Reduction of homeobox transcription factors in dicyemids is consistent with extreme simplification of their body plan. Unexpectedly, the dicyemids also experience several lineage-specific expansions of homeoboxes, notably the duplication of central Hox gene DoxC, which opposes the general trend of regulatory gene loss.

Basement Membrane

The basement membrane is a structure that enables the compartmentalization of cells to form tissues and organs. It is present in the majority of metazoans, with exception of sponges, placozoans, and acoelomorphs. The reported loss of a morphologically recognized basement membrane in dicyemids would indicate unprecedented simplification in this animal group. Even though this topic has been studied (Czaker, 2000), it is still unclear whether dicyemids have a basement membrane during any of their life cycle stages. The basement membrane consists of a set of “basement membrane toolkit” proteins, but the most important are collagen IV and laminin (Fidler et al., 2017). Both laminin and type IV collagen are multi-domain proteins that include specific domains (LamNT for laminin and C4 in the case of collagen type IV) and non-specific domains (EGF-like and other). The BLAST and Pfam searches showed that these domains of canonical molecules forming basement membrane are absent from the sequenced genome of Dicyema sp., therefore supporting the proposed secondary loss of this trait in dicyemids. The apparent absence of the recognized basement membrane is parallel with a reported loss of muscular and nervous systems in these animals. Indeed, in bilaterians, the basement membrane supports the maintainance of the muscular and nervous system architecture, their development and compartmentalization, and supporting growth factor signaling gradients among other functions.

The complete life cycle of dicyemids is not entirely understood, and more complex structures of transitional obscure life forms of these organisms are not excluded. An unknown stage can potentially exist between the infusorioform larvae that exits the host and the vermiform embryos found in cephalopods. The lack or reduced representation of genes encoding key elements of the basement membrane or other mediators of organ formation further supports the idea that dicyemids are secondarily simplified to an outstanding state.

Membrane Receptor Proteins

Cell surface membrane receptors act in cell signaling and allow communication between the cell and the extracellular space. Their diversity reflects the complexity of the organism and its ability to respond to different external signals. The number of genes encoding receptor proteins in dicyemids is exceptionally low. We found only two PF00001 domain hits corresponding to the 7 pass transmembrane receptor proteins of rhodopsin family in Dicyema sp. This family of G-protein-coupled receptors (GPCRs) is ubiquitously present and abundant in metazoans and contains tens to hundreds of members in different species. The minimum number of the rhodopsin family genes (six per genome) is detected in the sponge Amphimedon queenslandica; even the genome of the simplified orthonectid I. linei contains 32 genes of the rhodopsin family. The actual specificity of these GPCRs proteins is unknown, although their BLAST search shows best similarity to the rhodopsin family neuropeptide receptors from other animals. Four proteins from another GPCR 7 pass transmembrane receptor family – secretin family (PF00002) were predicted in the Dicyema sp. data. This is fewer than in most metazoans yet some flatworms have even fewer (Zamanian et al., 2011), and the Orthonectida have no such proteins. We found one putative metabotropic glutamate receptor with a PF00003 domain. Curiously, this metabotropic glutamate receptor also contains a (LIVBP)-like domain that is characteristic of ionotropic glutamate receptors. Two ionotropic glutamate receptors (iGluRs) that are ligand-gated ion channels activated by the neurotransmitter glutamate with Lig_chan (PF00060) domain were identified in Dicyema sp. One of them with a PF10613 (Lig_chan-Glu_bd) and another with PF01094 (ANF_receptor). Thus, both distinct types of glutamate receptors (ionotropic and metabotropic types) are present in Dicyema sp. It is well known that glutamate is often associated with non-neuronal signaling and is highly abundant in some animals that lack nervous systems (such as sponges and Trichoplax). Previously, we reported that iGluRs are absent in the genome of orthonectid (Mikhailov et al., 2016). Glutamate receptors are also found in plants and many other eukaryotes outside Metazoa (Turano et al., 2001).

Another big group of ionotropic receptors is the Cys-loop ligand-gated ion channel superfamily that is composed of nicotinic acetylcholine, GABA-A, GABA-A-ρ, glycine, 5-HT3, and zinc-activated (ZAC) receptors. We found 8 genes for this superfamily in Dicyema sp., identified by the specific transmembrane region domain (PF02932) and the ligand binding domain (PF02931). All these receptors are predicted to be nicotinic acetylcholine-like receptors.

Ion Channels

Despite the reported absence of muscles and neurons, tetrameric ion channels that are often associated with cellular electrical excitability are present in Dicyema sp. in numbers similar to the orthonectid I. linei (33 and 36 sequences with PF00520, and 11 and 9 with PF07885 in Dicyema and the orthonectid, respectively). Although unlike Orthonectida no signatures for voltage-gated sodium ion channel (Na_trans_assoc PF06512) were detected in Dicyema, Pfam analysis (for Ca_chan_IQ PF08763) and reciprocal BLAST searches indicates the presence of voltage-gated calcium ion channels in this animal group. The presence of such channels together with tetrameric potassium ion channels implies that electrical excitability in the form of action potentials might exist in dicyemid cells. Figure 9 provides a hypothetical schema of the intercellular communication and an analog of the neuromuscular junction in dicyemids. This structure may potentially assemble from key predicted proteins typical to many other metazoans.

FIGURE 9

Genes Encoding Putative Contractile/Muscular Elements

“True” muscle cells are absent in dicyemids and detection of the muscle-specific genes in these animals is interesting. Most of the core muscle proteins, including a type II myosin heavy chain (MyHC) motor protein was already present in unicellular eukaryotes before the origin of multicellular animals (Steinmetz et al., 2012). At the same time, the troponin complex appears to be a universal innovation of bilaterians. Troponin is a complex of three proteins (troponin C, troponin I, and troponin T). These proteins are detected in the dicyemid data by BLAST search, and the troponin domain PF00992 is found by Pfam search. The troponin complex is characteristic of skeletal and cardiac muscles, but not for smooth muscles. It appears that throughout the radical simplification in dicyemids that resulted in massive gene loss (including most of genes encoding the extracellular matrix ECM molecules) and in the absence of specialized muscle cells the troponins remain essential. The presence of troponins relates dicyemids to all other bilaterians with one remarkable exception – the orthonectid. In contrast to dicyemids and other bilaterians, the genome of orthonectid I. linei has no troponins despite having specialized muscles. Morphological data suggest that muscles in I. linei are similar to smooth muscles, so troponin was likely lost in I. linei, and its absence is a derived feature. At the same time another bilaterian hallmark – the myogenic regulatory factor (Myogenic Basic domain PF01586) – is present in the genome of I. linei, but was not detected in dicyemids. These findings support the mosaic evolution of many bilaterian traits, supporting the possibility of independent simplifications in these two parasitic lineages.

Gap Junctions and Adhesion Molecules

Gap junctions are a distinct type of intercellular communication channels. In Metazoa, the gap junction proteins belong to two unrelated families – connexins and pannexins (also known as innexins). The connexins are only found in chordates, while the pannexin family is widespread in invertebrates. The presence of gap junctions and innexin/pannexins in dicyemids was demonstrated earlier by transmission electron microscopy (TEM) (Furuya et al., 1997) and molecular cloning (Suzuki et al., 2010). BLAST and Pfam searches with our dicyemid data detected 21 hits with the innexin/pannexin-specific Pfam domain (PF00876) and no connexins. The number of dicyemid pannexins is similar to other invertebrates (25 in Caenorhabditis elegans, 13 in Drosophila melanogaster). It appears that unlike the highly reduced chemical signaling, direct intercellular communication via gap junctions is conserved in dicyemids.

Other hallmarks of multicellularity – the adhesion molecules and adherens junctions are retained in dicyemids and were demonstrated in these organisms earlier by TEM (Furuya et al., 1997). The universal metazoan proteins Integrin alpha and Integrin beta are detected in dicyemids in single copies; immunoglobulin domain is present in 6 sequences and Cadherin in 18 copies.

Axon Guidance Molecules and Their Receptors

The simplicity of the nervous system in Orthonectida is associated with a reduction of genes encoding components of axon guidance and synapse formation (Mikhailov et al., 2016). Dicyemids are presumably entirely deprived of the nervous system and follow the same trend of gene loss. Both animal groups lack genes encoding semaphorins, important neuronal pathfinding signaling molecules, and their receptors (plexins). Genes potentially involved in the nervous system development, such as Netrin, Ephrins, and Ephrin receptors are present in Orthonectida but were not identified in Dicyemids. Interestingly, the fasciclin domain (PF02469) is absent in the genome of I. linei, but we found its three orthologs in Dicyema. Fasciclin (FAS1 domain) is a cell adhesion domain found in neural cell adhesion molecules involved in axonal guidance in insects (Grenningloh and Goodman, 1992).

Peroxisome

The proteins and Pfam domains specific to peroxisome organelles, found in most metazoans, are absent from the dicyemid data. The peroxisomal proteins PEX3, PEX10, PEX12, and PEX19, mandatory for peroxisome function are apparently missing. Failure to detect these genes unequivocally suggests the absence of the organelle. Eight Pfam domains (PF01756, PF04088, PF04614, PF04882, PF05648, PF07163, PF09262, and PF12634) linked to peroxisome in the GO database¹ were not detected in Dicyema spp. In this respect, dicyemids are similar to Orthonectida and parasitic flatworms (Tsai et al., 2013).

Phylogenetic Analyses

To clarify the relationships of the two mesozoan groups, Orthonectida and Dicyemida, we used the sequenced transcriptomes of two unidentified species of Dicyema. We included the gene predictions of the orthonectid I. linei (Mikhailov et al., 2016) in the set of orthologous genes based on the dataset published by Struck et al. (2014). Given the high uncertainty in phylogenetic affinities of orthonectids and dicyemids, we extended taxonomic sampling by adding 30 spiralian taxa from available transcriptomes (see section “Materials and Methods”). Although the data broadly covers the spiralian diversity, several taxonomic groups are still missing or underrepresented in the complement of sequenced genes. To minimize missing data, we merged closely related species within several operational taxonomic units (OTUs) and produced the final matrix with 73 OTUs (69 OTUs for spiralian species) and 87,610 aa positions from 452 individual protein alignments. The proportion of missing data in the concatenated alignment is 40%.

Highly divergent sequences of mesozoans pose a formidable challenge for inference methods due to the confounding effect of long branched taxa on phylogenetic reconstructions. A recognized approach to tackle the long branch attraction (LBA) problem is to use a site-heterogeneous model of sequence evolution (Rodríguez-Ezpeleta et al., 2007). In the Bayesian tree constructed with PhyloBayes (Lartillot et al., 2013) under the site-heterogeneous CAT-GTR model, the dicyemid and orthonectid lineages form the longest branches, yet they do not group thus contradicting monophyly of the Mesozoa (Figure 10). We recovered the position of the orthonectid Intoshia within annelids with the posterior probability of 1.0. Specifically, the orthonectid forms a branch of the Pleistoannelida that comprises the annelid groups Errantia and Sedentaria (Weigert et al., 2014), while Owenia, Magelona, Chaetopteridae, and Phascolosoma (Sipuncula) occupy more basal positions in the annelid subtree.

FIGURE 10

The same analysis placed the dicyemid lineage near the base of a group uniting the Rouphozoa (Platyhelminthes, Gastrotricha) and Entoprocta + Cycliophora. However, the position of dicyemids in Bayesian inference is unstable. In about one-third of trees dicyemids were recovered as a sister group to the clade uniting Annelida, Nemertea, Lophophorata (Brachiopoda + Phoronida + Bryozoa), and Mollusca. In about 10% of trees the dicyemids branch off at the base of this group plus (Platyhelminthes + Gastrotricha) plus (Entoprocta + Cycliophora) (Figure 11, green branch). The grouping of Intoshia linei and Pleistoannelida has been observed in all summed trees. However, the exact position of the orthonectids relative to pleistoannelids is less certain in our analyses. The basal placement of the orthonectids is observed in 50% of trees, and the orthonectids were recovered as a sister group of Sedentaria or Errantia in 38 and 11% of trees, respectively.

FIGURE 11

The consensus Bayesian tree was obtained from four independent chains. The majority of bipartitions are shared across chains, but convergence on a single topology was not observed. Topologies in each chain uniquely reflect the concurrent hypotheses of spiralian relationships (Kocot, 2016; Kocot et al., 2017). All four chains of our Bayesian analysis are consistent in several major assemblages, including the Rouphozoa (Platyhelminthes, Gastrotricha), Gnathifera (Gnathostomulida, Micrognathozoa, and Syndermata), the sister relationship of Rouphozoa and Entoprocta + Cycliophora, and the basal position of Gnathifera relative to other spiralians. Importantly, in all topologies, the orthonectid is nested within the Annelida, and the dicyemid lineage is inferred sister to the assemblage of Platyhelminthes, Gastrotricha, Entoprocta, and Cycliophora (to the inclusion of Bryozoa in some topologies, Figure 11).

Alternative groupings obtained in our analysis include the Lophophorata (Brachiopoda, Phoronida, and Bryozoa) versus Polyzoa (Entoprocta, Cycliophora, and Bryozoa), and Vermizoa (Annelida, Nemertea) versus Nemertea + Mollusca (Figure 11). A comparison of topologies across chains based on site-wise likelihoods computed with PhyloBayes (the “sitelogl” option of the PhyloBayes readpb_mpi) under the CAT-GTR model and the approximately unbiased (AU) test (Shimodaira, 2002) show that the difference in likelihoods of topologies in chains 1–3 is not significant, but the topology likelihood in chain 4 is significantly lower (p-value = 0.01). Chain 2 (Figure 11) converges on a topology identical to the consensus four-chains topology (Figure 10) but its likelihood is lower than in chain 1 (non-significantly). Excluding chain 4 that failed the AU test and constructing the consensus with the three remaining chains does not affect the topology itself but only node supports due to eliminating the effects of non-monophyletic Lophophorata in chain 4 (Supplementary Figure S3).

The best scoring topology supports the monophyletic Lophophorata, the grouping of Annelida and Nemertea, and also the monophyly of macrodasyid and chaetonotid gastrotrichs, which frequently find themselves separate in our analyses (Figure 10). Maximum likelihood (ML) analyses of the same dataset with RAxML (Stamatakis, 2014) and IQ-TREE (Nguyen et al., 2015) produce a different view on the phylogeny of Mesozoa. The dicyemids and the orthonectid form a monophyletic group in ML trees with maximal support. In the RAxML analysis under the GTR model the monophyletic Mesozoa branch off with chaetonotid gastrotrichs, similarly to the result obtained in the recent phylogenomic analysis (Lu et al., 2017), but the support of the group is minimal (Supplementary Figure S4). In the IQ-TREE run under the C60 profile mixture model the monophyletic Mesozoa are found at the base of the Rouphozoa, again with weak support (64% of ultrafast bootstrap replicates) (see Supplementary Figure S5).

Although ML analyses show disagreement with the result of Bayesian inference, modeling of site-heterogeneity by the IQ-TREE profile mixture model does shed light on some spurious cases in spiralian relationships. The divergent annelid Myzostoma is correctly grouped with other annelids in the IQ-TREE analysis, in contrast with the RAxML tree where it forms a clade with long branches of the Rouphozoa, Gnathifera, and Mesozoa. The clustering of Rouphozoa and Gnathifera referred to as the Platyzoa, receives maximal support in the RAxML analysis but was previously shown to be artefactual (Struck et al., 2014; Laumer et al., 2015). This grouping is not inferred by both the IQ-TREE and Bayesian analyses. In contrast with IQ-TREE and PhyloBayes, the RAxML tree supports monophyletic Polyzoa, a group uniting Entoprocta, Cycliophora, and Bryozoa, which was also suggested to be erroneous and caused by the compositional bias in amino acid sequences (Nesnidal et al., 2013).

To test for expected LBA effects, particularly to exclude the possibility of the orthonectid being attracted to annelids by the divergent Myzostoma, we conducted additional analyses excluding each of the long branched lineages. Additional datasets were generated by removing Myzostoma, Myzostoma and both dicyemids, Myzostoma and Intoshia. Bayesian analyses of additional datasets recovered placement of Intoshia within annelids in the absence of Myzostoma (Supplementary Figures S6, S7). The position of dicyemids is also unaffected by the exclusion of other long-branched taxa – the dicyemids occupy a basal position within the Lophotrochozoa after the divergence of Gnatifera in all analyses of the additional datasets (Supplementary Figures S6, S8).

We also tested our dataset for the effects of compositional heterogeneity by discarding highly heterogeneous alignments and utilizing the data recoding approach (Susko and Roger, 2007). Bayesian inference with a concatenate of 150 protein alignments retained after discarding highly compositionally heterogeneous alignments from the original dataset recovers the same groupings of the mesozoan taxa as the analysis of the full dataset. The orthonectid is nested within the annelid clade (1.0 posterior probability) and the dicyemids branch off at the base of the Rouphozoa + Entoprocta + Cycliophora clade with weak support (0.46 posterior probability) (Supplementary Figure S9). Similarly, inference with the Dayhoff-recoded alignment groups the orthonectid with annelids, while leaving the position of the dicyemids uncertain within the Lophotrochozoa (Supplementary Figure S10). Remarkably, the PhyloBayes run with recoded data shows adequate convergence between chains (maxdiff = 0.17) and infers the monophyletic Gastrotricha. Several conventional groupings, such as the Rouphozoa, are not recovered. Consistent with the proposed artefactual nature of the grouping of Bryozoa and Entoprocta due to compositional heterogeneity (Nesnidal et al., 2013), both test datasets support the monophyletic Lophophorata, whereas the alternative Polyzoa was frequently observed for the complete and non-recoded data. The lack of support for monophyletic Rouphozoa in the analysis with the recoded dataset was similarly obtained in a recent study of the spiralian phylogeny, which aimed at counteracting the impact of compositional heterogeneity (Marlétaz et al., 2019).

Schiffer et al. (2018) selected proteins that support annelid monophyly as an approach to verify the orthonectid position. We also selected 111 protein alignments that contain the annelid signal but with a different method, and used those for Bayesian inference with the PhyloBayes program. In contrast to other Bayesian runs, the consensus presents a stable topology (maxdiff value 0.15). In this tree, the orthonectid I. linei stabilizes inside the annelids [posterior probability (PP) 1.0], whereas the species of Dicyema are not attracted to annelids (Figure 12). The position of Dicyema remains uncertain within the lineage of long-branched taxa (Platyhelminthes, Gastrotricha, Entoprocta, Cycliophora). Lophophorata and Gastrotricha are reconstructed with PP 1.0 (as in case of the Dayhoff-recoded dataset mentioned above and the non-recoded full dataset in chain 1 that reaches the highest likelihood). Bayesian topologies obtained in chain 1 (Figure 11) and both the sub-sampled datasets of 111 proteins with the annelid signal and the 150 proteins with low compositional heterogeneity (Supplementary Figure S9) reconstruct the sister relationship of annelids with nemertines. Marlétaz et al. (2019) also report the grouping of annelids and nemertines, with the inclusion of Platyhelminthes.

FIGURE 12

The lack of convergence in most PhyloBayes analyses precludes strong assertions regarding problematic areas of the spiralian tree. Nevertheless, some clades are reconstructed consistently. We do not observe monophyly of the Mesozoa in any of the chains, in contrast to the recent study by Lu et al. (2017) and in agreement with the evidence from Schiffer et al. (2018), nor do we observe their direct relations with Platyhelminthes.

The orthonectid I. linei occupies a stable position within the annelid part of the tree. Its placement is among the major conflicts between ML and Bayesian topologies, which likely indicates the impact of a more complex CAT-GTR model in the presence of long branches of highly divergent lineages like orthonectids and dicyemids. Noteworthy, polyphyletic Mesozoa and the proposed affinity of orthonectids to annelids was also recovered in Schiffer et al. (2018) in Bayesian analyses of a dataset with the less extensive representation of the lophotrochozoan diversity.

The orthonectids share with annelids certain morphological features: the presence of microvillar cuticle, metameric muscles, gonochory (Slyusarev, 2008), and the dorsal ganglion in adult specimens (Slyusarev and Starunov, 2015). Cases of dramatic morphological reduction in annelids are known in archiannelids (Andrade et al., 2015) and lobatocerebrids (Laumer et al., 2015), and especially in dwarf males of the echiurid Bonnelia, dinophidid Dinophilus gyrociliatus, spionid Scolelepis laonicola (Vortsepneva et al., 2008), and siboglinids (Worsaae and Rouse, 2010). Adaptations of orthonectids that had led to the complete loss of coelomic cavity, gonad wall, chaetae, gastral system, nephridia, trochophore larva, and spiral cleavage further demonstrate the extent of morphological regress in the evolution of annelids.

The dicyemid lineage in our analyses exhibits affinity to the Rouphozoa clade, in congruence with Lu et al. (2017), with the intercalation of Entoprocta or Polyzoa, which were not included in analyses by Lu et al. (2017). Schiffer et al. (2018) report an uncertain position of dicyemids at the base of the Lophotrochozoa. We did not recover the position of the dicyemids as part of the Platyhelminthes. A previous analysis of innexin genes (Suzuki et al., 2010) also rejects the kinship of the dicyemids and Platyhelminthes. Furthermore, rhabditophoran platyhelminths are known to possess a unique non-standard mitochondrial genetic code, which was shown to be not the case for the studied mesozoans (Telford et al., 2000; Schiffer et al., 2018). Dicyemids were placed within Spiralia in various taxonomic contexts in molecular phylogenetic studies (Pawlowski et al., 1996; Petrov et al., 2010; Suzuki et al., 2010; Lu et al., 2017) but the interpretation of their body plan remains enigmatic. Being among the simplest known bilaterians, they yet possess multiciliated epithelia, which is not a primitive trait and suggests secondary evolutionary regress, and do not display evident synapomorphies with other animal phyla. Dicyemids might represent a relict lineage of lophotrochozoan animals with no direct relatives that had survived to the present days.

Conclusion

We confirm that orthonectids are extremely simplified annelids and do not form a monophyletic group with dicyemids. Mesozoa is a polyphyletic taxon. Dramatic simplification of their body plan, as well as the metagenetic life cycle, evolved independently in the two lineages. Many conserved bilaterian genes are absent in the genomes of Dicyemida and Orthonectida. At the same time, the pattern of their loss and presence is different, which supports the conclusion that these animal groups are not close relatives and have simplified independently. Analyses of genes related to the basement membrane, neuronal and muscular systems expose the extreme simplicity of dicyemids. Intriguingly, dicyemids lack muscle cells and the genetic factors of muscle cell differentiation but possess the troponin complex specific for striated muscles. Taken together with detection of a relatively big set of nicotinic acetylcholine receptors often associated with neuromuscular signaling and the presence of voltage-gated ion channels, this fact urges reevaluation of the traditional view that dicyemids completely lost the neuro-muscular system. Appealing is to experimentally check if some contractility and movements could be induced in dicyemids by signal molecules such as acetylcholine or glutamate, and for the presence of electrical excitability in the form of propagated calcium action potential in their cells. Small circular extrachromosomal molecules are present in total DNA extracts of dicyemids. Mitochondrial rRNA, tRNA, protein-coding genes and pseudogenes are located on circular molecules. There are short nucleotide sequence motifs confined specifically to circular DNAs in Dicyema sp.

Materials and Methods

Biological Material, Genome and Transcriptome Sequencing

The original live material on Dicyema sp. 1 was collected at the Vostok marine biological station of the Institute for Marine Biology of the Russian Academy of Sciences (the Vostok Bay of the Sea of Japan, Vladivostok, Russia) from dissected kidneys of the giant Pacific octopus E. dofleini. Live dicyemids were rinsed individually in filtered marine water and fixed in the RNAlater stabilization solution (Ambion). Total DNA was isolated from tissue samples by Diatom DNA Prep (IsoGene). The sequencing of dicyemid genomic data was performed with an Illumina HiSeq2000 system, generating 140 million paired-end reads.

Total RNA was isolated by TRIzol kit (Invitrogen) and further used for ds cDNA synthesis using the SMART approach (Zhu et al., 2001). SMART-prepared amplified cDNA was then normalized using the DSN normalization method (Zhulidov et al., 2004). Normalization included cDNA denaturation/reassociation, treatment by the duplex-specific nuclease (Shagin et al., 2002), and PCR amplification of the normalized fraction (8 PCR cycles: 95°C for 7 s; 65°C for 20 s; 72°C for 3 min). Normalized cDNA libraries were sequenced using the Roche 454 sequencing technology, producing about 480,000 reads with an average length of 444.6 bases.

Specimens of Dicyema sp. 2 was collected at the Friday Harbor Laboratories (Friday Harbor, WA, United States) from circulatory system and kidneys of the octopus E. dofleini. All individual animals were washed 3–5 times in 0.2 μm filtered seawater. Then RNA was extracted from individual animals and processed as described elsewhere (Moroz and Kohn, 2013) for Illumina HiSeq 2000 sequencing.

The sequences are deposited in the NCBI: BioProject PRJNA527259 (Dicyema sp. 1) and SRA SRP021079 (Dicyema sp. 2).

Assembly and Filtering of Dicyemid Sequences

The reads obtained from the DNA library for Dicyema sp. 1 were trimmed for adapters with Trimmomatic (Bolger et al., 2014) and assembled by SPAdes (Nurk et al., 2013) using k-mer values of 21, 33, 55, and 77. We also performed genome assembly with the Newbler GS De Novo Assembler software (v. 2.9) (using 1/10 of all reads) as a control to our method of circular contigs identification. Gene prediction was performed with Augustus (Stanke and Waack, 2003) after constructing a training set of 200 dicyemid sequences identified in the genomic assembly. The predicted genes were queried against the InterPro database (Finn et al., 2017) with InterProScan (Jones et al., 2014) and genes with InterPro hits were screened for cephalopod sequences with BLAST (Altschul et al., 1997) searches against the NCBI nr database. Predictions producing best hit with cephalopod sequences were discarded from the gene set. Completeness estimates were performed with BUSCO (Waterhouse et al., 2017) using the eukaryota_odb9 ortholog set (Zdobnov et al., 2017). HMMER (Eddy, 2011) searches were carried out with Pfam (Finn et al., 2016) Homeodomain (PF00046) and Homeobox_KN (PF05920) profiles to identify homeobox transcription factors in the data. Phylogeny reconstructions for homeobox sequences were performed with IQ-TREE (Nguyen et al., 2015) using the LG + C20 + G4 model of sequence evolution or with PhyloBayes (Lartillot et al., 2013) using the LG +CAT + G4 model.

The reads obtained from the cDNA library for Dicyema sp. 1 were trimmed for adapters, non-coding RNA, low-quality and low-complexity sequences with the SeqClean software (Dana-Farber Cancer Institute²), and about 430,000 reads were retained. Data was further assembled with the original 454 Newbler GS De Novo Assembler software (v. 2.9) utilizing flowgram quality data and settings that maximize contig overlap. The “-urt” option was invoked to improve contigging in low depth portions of the assembly. Fusions of transcripts that can potentially occur with low-depth assembly extensions in densely packed genomes are subsequently eliminated in our experimental design by alignment filtering at the supermatrix construction step. The obtained assembly contained 19,641,638 bases, and 22,082 isotigs of average size 889 bases, N50 size of 1,081, and the largest isotig size of 9,199. Protein coding regions were predicted using TransDecoder (Haas et al., 2013) with settings to maximize the sensitivity of capturing ORFs regardless of the predicted coding likelihood score by accounting for homology to known proteins in the Pfam (Finn et al., 2014) and UniProtKB/Swiss-Prot (UniProt Consortium, 2015) curated databases. Coding region prediction with TransDecoder was set to the minimal predicted protein length of 80 aa. The predicted proteome contained 15,227 unique coding regions.

The second dicyemid transcriptome sequenced using the Illumina platform was assembled with Trinity (Grabherr et al., 2011). Before assembly the reads were processed with the SeqClean software, and the prediction of coding regions was performed by TransDecoder, similarly to the transcriptome of the first dicyemid.

The transcriptomes of dicyemids were derived from samples contaminated with their cephalopod host. Therefore, we paid special attention to avoid mixing dicyemid and cephalopod sequences in the phylogenetic analysis. The transcriptomes of dicyemids were first screened for cephalopod sequences by performing BLAST (Altschul et al., 1997) searches against the NCBI RefSeq database (O’Leary et al., 2016). Two dicyemid transcriptomes were processed independently. In the first step of decontamination we filtered out proteins having best hit in RefSeq belonging to prokaryotes (and having at least 50% identity). This lead to rejection of only 35 proteins for Dicyema sp. 1 and 363 proteins for Dicyema sp. 2. In the second step we removed all dubious proteins if their local alignment score with any cephalopod protein higher than in all the other considered species (with the same query protein). Sequences with best hits to cephalopods were discarded from the transcriptomes if the sequence identity exceeded 70%. For the third filtering step we queried the proteins of O. bimaculoides combined with several transcriptomes of Octopus vulgaris (NCBI BioProject PRJNA79361 and Sequence Read Archive entries SRR331946, SRR1507221) against a custom database containing 9 metazoan proteomes (4 molluscs, 2 annelids, a brachiopod Lingula anatina, an ecdysozoan Limulus polyphemus, and a deuterostome Danio rerio) and the dicyemid transcriptome, and inspected dicyemid sequences that produced hits with the highest match to the cephalopod queries among the 10 metazoans. All dubious sequences (hits with at least 80% identity) captured by this method were discarded from the dicyemid transcriptomes as potential cephalopod contamination.

Search for “Circular” Contigs, Signals, and Mitochondrial Sequences

The contigs constructed from shotgun fragments display special characteristics emerging from the genome assembly algorithms based on De Bruijn graph of k-mers. This approach results in “circular” contigs starting and ending with the same k-mer. After assembly, terminal repeats equal in length to the k-mer were cut off. Contigs analyzed in sections “Circular Contigs in Genomic Assembly of Dicyema sp.” and “Mitochondrial DNA of Dicyema sp.,” and NCBI submission data have been cleaned off the terminal repeats. In this study, a contig was considered “circular” if it had terminal direct repeats ≥ 77 nt in length (k77). The length distribution of contigs assembled by different methods (Newbler and SPAdes) was compared with the two-sample Kolmogorov–Smirnov test implemented in the SciPy package in Python 3. Here the null hypothesis is that contig lengths come from the same distribution. High p-values in this case reflect high probabilities of this hypothesis. Low complexity regions were detected with the DUST algorithm from the MEME Suite (Bailey and Elkan, 1994) with standard settings. MEME and ChIPMunk (Kulakovskiy et al., 2010) tools with the default parameters were applied to the task of finding specific motifs. The reverse lookup for the signal presence was done via FIMO (from the MEME Suite) with the p-value threshold of 10⁻⁴. Moreover, highly conserved elements of circles in dicyemids were found utilizing the technique borrowed from (Rubanov et al., 2016). The method identifies highly conserved DNA elements on the base of the identification of dense subgraphs in a specially built multipartite graph (whose parts correspond to genomes). Specifically, the algorithm does not rely on genome alignments, no pre-identified perfectly conserved elements; instead, it performs a fast search for pairs of words (in different genomes) of maximum length with the difference below the specified edit distance. Such pair defines an edge whose weight equals the maximum (or total) length of words assigned to its ends. The graph composed of these edges is then compacted by merging some of its edges and vertices. The dense subgraphs are identified by a cellular automaton-like algorithm; each subgraph defines a cluster composed of similar inextensible words from different genomes (Rubanov et al., 2016).

HMMER3 package (Eddy, 2011) along with the Pfam-A database were used to find the circles containing protein-coding sequences, whereas an additional verification step was performed in BLAST. The search itself was conducted through the database composed of six-frame translated circular sequences. The search for genes coding for mitochondrial proteins was conducted with BLAST using mitochondrial protein-coding gene sequences from flatworms as queries, MITOS (Bernt et al., 2013) and HMMER3 using HMM profiles from the Pfam-A database. Mitochondrial rrnS genes in dicyemids are highly diverged and poorly detected with BLAST. Their detection was conducted with HMMER3 with HMM profiles preliminarily generated from the set of 140 rrnS alignments from other organisms (140 species of bilaterians, cnidarians, and placozoans). All findings were verified using blastp or blastn with nr NCBI database. It was proposed that the dicyemid small mitochondrial circular DNA molecules are generated from the usual long multigene mitochondrial DNA (Awata et al., 2006). If such long mtDNA exists together with mitochondrial mini-circles we can expect the cases when one read from the sequencing library corresponds to a particular mitochondrial mini-circle while its pair read maps elsewhere. Blastn with minimal word size was used to map raw paired end reads to circular contigs coding for mitochondrial genes to search for hypothetical high-molecular-weight mtDNA. Reads pair analysis was conducted after that in order to find the reads whose pair does not map to the initial circular contig. Mitochondrial tRNA secondary structures were predicted using the MiTFi program (Jühling et al., 2012).

Taxonomic Expansion of Alignments

The starting set of orthologous genes used in this work is based on a dataset for phylogenetic reconstructions within Spiralia assembled by Struck et al. (2014) that was later expanded with sequences of orthonectid I. linei (Mikhailov et al., 2016). The base set of orthologs contained 469 alignments with a total of 62 spiralian species and four ecdysozoan species. To extend the taxonomic sampling of Spiralia and minimize the missing data in the dataset we obtained predicted proteins from several genomic projects accessible through public databases and collected transcriptomic data from the NCBI Sequence Read Archive. The annotations for the genomes of Clonorchis sinensis, Echinococcus granulosus, L. anatina, O. bimaculoides, Priapulus caudatus were obtained from the GenBank database, and the proteins of Adineta vaga were obtained from the Genoscope database. The NCBI Sequence Read Archive was used to extract raw sequence data of another 31 spiralian species (see Supplementary Table S1).

The assemblies of the SRA transcriptome data were performed with Trinity (Grabherr et al., 2011) after cleaning the reads with SeqClean (Dana-Farber Cancer Institute³) from adapter sequences using the UniVec_Core database⁴ and filtering ribosomal RNA sequences using a database of eukaryotic rRNAs. The prediction of proteins in the assembled transcripts was performed with TransDecoder (Haas et al., 2013), which was assisted with searches against the Pfam (Finn et al., 2014) and UniProtKB/Swiss-Prot (UniProt Consortium, 2015) databases.

The addition of proteins from the newly assembled data to orthologous groups featured in the base set of alignments was performed using the procedure for mapping genes to existing orthologous groups (Fischer et al., 2011) of the OrthoMCL database (Chen et al., 2006). The genes from the initial dataset and novel transcriptomic and genomic data were assigned to orthologous groups of OrthoMCL-DB, and the genes within the same orthologous group were extracted and aligned together using MUSCLE (Edgar, 2004). When more than one sequence per organism was assigned to the same group of orthologs, only the sequence scoring highest against the orthologous group in the initial dataset was selected for the alignment.

Phylogenetic Analyses

The concatenation of individual gene alignments was performed with Scafos (Roure et al., 2007) using the option to construct chimeric sequences for several closely related taxa. The following 15 chimeric taxa were constructed for the analysis: Aplysia californica + Biomphalaria glabrata, Brachionus plicatilis + B. manjavacas, Chiton olivaceus + Chaetopleura apiculata, Clonorchis sinensis + Opisthorchis viverrini, Dugesia japonica + Dugesia ryukyuensis, Echinococcus granulosus + Echinococcus multilocularis, Echinorhynchus gadi + Echinorhynchus truttae, Euprymna scolopes + Idiosepius paradoxus, Lepadella patella + Lecane inermis, Pedicellina sp. + P. cernua, Protodrilloides symbioticus + P. chaetifer, Schistosoma mansoni + S. japonicum, Spiochaetopterus sp. + Chaetopterus variopedatus, Stenostomum leucops + Stenostomum sthenum, Symbion pandora + S. americanus. Another ten species that were present in the starting set of alignments were removed due to poor representation in the final alignment: Alcyonidium diaphanum, Fasciola gigantica, Flustra foliacea, Lumbricus rubellus, Philodina roseola, Rotatoria rotatoria, Spirometra erinacei, Stylochoplana maculata, Taenia solium, Turbanella ambronensis. The final number of operational taxonomic units featured in the analysis is 73. Before concatenation, the alignments were trimmed with TrimAl (Capella-Gutierrez et al., 2009) to remove poorly aligned regions. The trimming was performed with a gap threshold of 0.5 and a similarity threshold of 0.001. After the removal of invariant positions, the length of the concatenated alignment totaled 87,610 positions, with 40% missing data. Compositional heterogeneity in the alignment partitions (i.e., individual protein alignments after masking) was evaluated using the relative composition frequency variability (RCFV) metric (Zhong et al., 2011). The RCFV values were calculated using BaCoCa (Kuck and Struck, 2014). The low compositional heterogeneity dataset was generated by discarding 302 partitions (referred to in the paper simply as protein alignments) with RCFV value exceeding 0.115.

The phylogenetic reconstructions were performed with PhyloBayes-MPI 1.7 (Lartillot et al., 2013), RAxML (Stamatakis, 2014), and IQ-TREE (Nguyen et al., 2015). The RAxML analysis was carried out utilizing the complete analysis function (-f a) with 150 rapid bootstrap replicates and the PROTCATGTR model of evolution. The IQ-TREE analysis was performed using the LG + C60 + F + G4 evolutionary model, and node support was calculated using the ultrafast bootstrap approximation (Minh et al., 2013) with 1,000 replicates. The Bayesian inference with PhyloBayes was carried out using the CAT + GTR + Γ4 model, and the analyses were run with four chains. For the main dataset, the majority rule consensus tree was reconstructed after 30,000 cycles using one out of ten cycles with a 60% burn-in. PhyloBayes analyses of the additional datasets were conducted similarly to the main dataset; the consensus trees were reconstructed after 5,000 or 15,000 cycles with a 50% burn-in. Analysis of the recoded alignment was performed with PhyloBayes utilizing the recode option and the Dayhoff recoding scheme with six amino acid groups (Dayhoff et al., 1978). The convergence of the chains was assessed by comparing bipartitions using the pbcomp utility from the PhyloBayes package.

Comparison of topologies in the four chains of the Bayesian inference of the main dataset was performed using the CONSEL program (Shimodaira and Hasegawa, 2001) and the “sitelogl” option of the PhyloBayes readpb_mpi program. The site-specific marginal log likelihoods were computed for each chain across 10 data points sampled over 2,000 cycles after a 20,000 cycle burn-in.

Alignment partitions (i.e., individual protein alignments after masking) with the strong annelid signal were selected as follows. In a protein alignment we define two sets of sequences – G₁ (ingroup), and G₂ (outgroup). Only alignment positions containing no more than a half of missing data (gaps or X’s) in each of the two sets are considered. For each such position i-value q(i) is determined as the maximum of frequency differences of each amino acid in this position from G₁ and G₂. Missing data is ignored. Maximum q(i) value is 1 when G₁ consists only of one character, and G₂ does not contain this character. Under any q(i) ½ there exists a character a(i) observed in more than a half of taxa from G₁ but much less frequently in G₂ [frequency difference is q(i) ½]. In the phylogenetic context, when G₁ + G₂ constitute a monophyletic clade, and G₁ is a narrower natural clade, high q(i) values can be interpreted as presence of a synapomorphy against G₂. Notably, in this analysis q(i) values are used only to select partitions but not for alignment editing or positions removal. In our case of detecting the annelid signal, G₁ contained all annelids except the orthonectid, and G₂ – all non-annelid taxa except dicyemids in order to obtain q(i) estimates unbiased with respect to the lineages under study.

References

• 1

  AltschulS. F.MaddenT. L.SchafferA. A.ZhangJ.ZhangZ.MillerW.et al (1997). Gapped BLAST and PSI-BLAST: a new generation of protein database search programs.Nucleic Acids Res.253389–3402. 10.1093/nar/25.17.3389

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 2

  AndradeS. C.NovoM.KawauchiG. Y.WorsaaeK.PleijelF.GiribetG.et al (2015). Articulating “Archiannelids”: phylogenomics and annelid relationships, with emphasis on meiofaunal taxa.Mol. Biol. Evol.322860–2875. 10.1093/molbev/msv157

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 3

  ArmstrongM. R.BlokV. C.PhillipsM. S. (2000). A multipartite mitochondrial genome in the potato cyst nematode Globodera pallida.Genetics154181–192.

  • Google Scholar
  • View reference in article

• 4

  ArugaJ.OdakaY. S.KamiyaA.FuruyaH. (2007). Dicyema Pax6 and Zic: tool-kit genes in a highly simplified bilaterian.BMC Evol. Biol.7:201. 10.1186/1471-2148-7-201

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 5

  AtaevG. L. (2017). Reproduction of Trematode Parthenites: An Overview of the Main Theories.St. Petersburg: Nauka, 87.

  • Google Scholar
  • View reference in article

• 6

  AwataH.NotoT.EndohH. (2006). Peculiar behavior of distinct chromosomal DNA elements during and after development in the dicyemid mesozoan Dicyema japonicum.Chromosome Res.14817–830. 10.1007/s10577-006-1084-z

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 7

  BaileyT. L.ElkanC. (1994). Fitting a mixture model by expectation maximization to discover motifs in biopolymers.Proc. Int. Conf. Intell. Syst. Mol. Biol.228–36.

  • Pubmed Abstract
  • Google Scholar
  • View reference in article

• 8

  BerntM.DonathA.JühlingF.ExternbrinkF.FlorentzC.FritzschG.et al (2013). MITOS: improved de novo metazoan mitochondrial genome annotation.Mol. Phylogenet. Evol.69313–319. 10.1016/j.ympev.2012.08.023

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 9

  BolgerA. M.LohseM.UsadelB. (2014). Trimmomatic: a flexible trimmer for Illumina sequence data.Bioinformatics302114–2120. 10.1093/bioinformatics/btu170

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 10

  BondarenkoN.BondarenkoA.StarunovV.SlyusarevG. (2019). Comparative analysis of the mitochondrial genomes of Orthonectida: insights into the evolution of an invertebrate parasite species.Mol. Genet. Genomics .

  • Pubmed Abstract
  • Google Scholar
  • View reference in article

• 11

  BruscaR. C.BruscaG. J. (2003). Invertebrates.Sunderland, MA: Sinauer Associates.

  • Google Scholar
  • View reference in article

• 12

  BurgerG.JacksonC. J.WallerR. F. (2012). “Unusual mitochondrial genomes and genes,” in Organelle Genetics, ed.BullerwellC. (Berlin: Springer), 41–77. 10.1007/978-3-642-22380-8_3

  • CrossRef
  • Google Scholar
  • View reference in article

• 13

  Capella-GutierrezS.Silla-MartinezJ. M.GabaldonT. (2009). trimAl: a tool for automated alignment trimming in large-scale phylogenetic analyses.Bioinformatics251972–1973. 10.1093/bioinformatics/btp348

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 14

  CatalanoS. R.WhittingtonI. D.DonnellanS. C.BertozziT.GillandersB. M. (2015). First comparative insight into the architecture of COI mitochondrial minicircle molecules of dicyemids reveals marked inter-species variation.Parasitology1421066–1079. 10.1017/S0031182015000384

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 15

  ChenF.MackeyA. J.StoeckertC. J.Jr.RoosD. S. (2006). OrthoMCL-DB: querying a comprehensive multi-species collection of ortholog groups.Nucleic Acids Res.34D363–D368. 10.1093/nar/gkj123

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 16

  ChengT. C. (1986). General Parasitology, 2nd Edn.Cambridge, MA: Academic Press. 827.

  • Google Scholar
  • View reference in article

• 17

  CzakerR. (2000). Extracellular matrix (ECM) components in a very primitive multicellular animal, the dicyemid mesozoan Kantharella antarctica.Anat. Rec.25952–59. 10.1002/(SICI)1097-0185(20000501)259:1<52::AID-AR6<3.0.CO;2-J

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 18

  DayhoffM. O.SchwartzR. M.OrcuttB. C. (1978). “A model of evolutionary change in proteins,” in Atlas of Protein Sequence and StructureVol. 5ed.DayhoffM. O. (Washington, D.C: National Biomedical Research Foundation), 345–352.

  • Google Scholar
  • View reference in article

• 19

  de RosaR.GrenierJ. K.AndreevaT.CookC. E.AdoutteA.AkamM.et al (1999). Hox genes in brachiopods and priapulids and protostome evolution.Nature399772–776. 10.1038/21631

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 20

  DobrovolskijA. A.AtaevG. L. (2003). “Thr nature of reproduction of digenea rediae and sprosycts,” in Taxonomy, Ecology and Evolution of Metazoan Parasites, edsCombesC.JourdaneJ.Ducreux-ModatA.PagèsJ. R. (Perpignan: University of Perpignan Press), 273–290.

  • Google Scholar
  • View reference in article

• 21

  EddyS. R. (2011). Accelerated profile HMM searches.PLoS Comput Biol.7:e1002195. 10.1371/journal.pcbi.1002195

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 22

  EdgarR. C. (2004). MUSCLE: multiple sequence alignment with high accuracy and high throughput.Nucleic Acids Res.321792–1797. 10.1093/nar/gkh340

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 23

  FidlerA. L.DarrisC. E.ChetyrkinS. V.PedchenkoV. K.BoudkoS. P.BrownK. L.et al (2017). Collagen IV and basement membrane at the evolutionary dawn of metazoan tissues.Elife6:e24176. 10.7554/eLife.24176

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 24

  FinnR. D.AttwoodT. K.BabbittP. C.BatemanA.BorkP.BridgeA. J.et al (2017). InterPro in 2017-beyond protein family and domain annotations.Nucleic Acids Res.45D190–D199. 10.1093/nar/gkw1107

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 25

  FinnR. D.BatemanA.ClementsJ.CoggillP.EberhardtR. Y.EddyS. R.et al (2014). Pfam: the protein families database.Nucleic Acids Res.42D222–D230. 10.1093/nar/gkt1223

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 26

  FinnR. D.CoggillP.EberhardtR. Y.EddyS. R.MistryJ.MitchellA. L.et al (2016). The Pfam protein families database: towards a more sustainable future.Nucleic Acids Res.44D279–D285. 10.1093/nar/gkv1344

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 27

  FischerS.BrunkB. P.ChenF.GaoX.HarbO. S.IodiceJ. B.et al (2011). Using OrthoMCL to assign proteins to OrthoMCL-DB groups or to cluster proteomes into new ortholog groups.Curr. Protoc. Bioinform.611–19. 10.1002/0471250953.bi0612s35

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 28

  FuruyaH.HochbergF. G.TsunekiK. (2004). Cell number and cellular composition in infusoriform larvae of dicyemid mesozoans (phylum Dicyemida).Zool. Sci.20877–889. 10.2108/zsj.21.877

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 29

  FuruyaH.TsunekiK. (2003). Biology of dicyemid mesozoans.Zool. Sci.20519–532. 10.2108/zsj.20.519

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 30

  FuruyaH.TsunekiK.KoshidaY. (1997). Fine structure of dicyemid mesozoans, with special reference to cell junctions.J. Morphol.231297–305. 10.1002/(SICI)1097-4687(199703)231:3<297::AID-JMOR8<3.0.CO;2-8

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 31

  GibsonT.BlokV. C.DowtonM. (2007). Sequence and characterization of six mitochondrial subgenomes from Globodera rostochiensis: multipartite structure is conserved among close Nematode relatives.J. Mol. Evol.65308–315. 10.1007/s00239-007-9007-y

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 32

  GrabherrM. G.HaasB. J.YassourM.LevinJ. Z.ThompsonD. A.AmitI.et al (2011). Full-length transcriptome assembly from RNA-Seq data without a reference genome.Nat. Biotechnol.29644–652. 10.1038/nbt.1883

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 33

  GrenninglohG.GoodmanC. S. (1992). Pathway recognition by neuronal growth cones: genetic analysis of neural cell adhesion molecules in Drosophila.Curr. Opin. Neurobiol.242–47. 10.1016/0959-4388(92)90160-m

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 34

  HaasB. J.PapanicolaouA.YassourM.GrabherrM.BloodP. D.BowdenJ.et al (2013). De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis.Nat. Protoc.81494–1512. 10.1038/nprot.2013.084

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 35

  HaneltB.Van SchyndelD.AdemaC. M.LewisL. A.LokerE. S. (1996). The phylogenetic position of Rhopalura ophiocomae (Orthonectida) based on 18S ribosomal DNA sequence analysis.Mol. Biol. Evol.131187–1191. 10.1093/oxfordjournals.molbev.a025683

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 36

  HymanL. H. (1940). The Invertebrates: Protozoa Through Ctenophra, Vol. 1. New York, NY: McGraw Hill

  • Google Scholar
  • View reference in article

• 37

  JoffroyB.UcaY. O.PrešernD.DoyeJ. P. K.SchmidtT. L. (2018). Rolling circle amplification shows a sinusoidal template length-dependent amplification bias.Nucleic Acids Res.46538–545. 10.1093/nar/gkx1238

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 38

  JonesP.BinnsD.ChangH. Y.FraserM.LiW.McAnullaC.et al (2014). InterProScan 5: genome-scale protein function classification.Bioinformatics301236–1240. 10.1093/bioinformatics/btu031

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 39

  JühlingF.PutzJ.BerntM.DonathA.MiddendorfM.FlorentzC.et al (2012). Improved systematic tRNA gene annotation allows new insights into the evolution of mitochondrial tRNA structures and into the mechanisms of mitochondrial genome rearrangements.Nucleic Acids Res.402833–2845. 10.1093/nar/gkr1131

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 40

  KatayamaT.WadaH.FuruyaH.SatohN.YamamotoM. (1995). Phylogenetic position of the dicyemid mesozoa inferred from 18S rDNA sequences.Biol. Bull.18981–90. 10.2307/1542458

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 41

  KiełbasaS. M.WanR.SatoK.HortonP.FrithM. C. (2011). Adaptive seeds tame genomic sequence comparison.Genome Res.21487–493. 10.1101/gr.113985.110

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 42

  KobayashiM.FuruyaH.HollandP. W. (1999). Dicyemids are higher animals.Nature401:762. 10.1038/44513

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 43

  KobayashiM.FuruyaH.WadaH. (2009). Molecular markers comparing the extremely simple body plan of dicyemids to that of lophotrochozoans: insight from the expression patterns of Hox, Otx, and brachyury.Evol. Dev.11582–589. 10.1111/j.1525-142X.2009.00364.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 44

  KocotK. M. (2016). On 20 years of Lophotrochozoa.Org. Divers. Evol.16329–343. 10.1007/s13127-015-0261-3

  • CrossRef
  • Google Scholar
  • View reference in article

• 45

  KocotK. M.StruckT. H.MerkelJ.WaitsD. S.TodtC.BrannockP. M.et al (2017). Phylogenomics of Lophotrochozoa with consideration of systematic error.Syst Biol.66256–282. 10.1093/sysbio/syw079

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 46

  KolesnikovA. A.GerasimovE. S. (2012). Diversity of mitochondrial genome organization.Biochemistry (Moscow)771424–1435. 10.1134/S0006297912130020

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 47

  KozloffE. N. (1990). Invertebrates.Philadelphia: Saunders College Publishing.

  • Google Scholar
  • View reference in article

• 48

  KuckP.StruckT. H. (2014). BaCoCa–a heuristic software tool for the parallel assessment of sequence biases in hundreds of gene and taxon partitions.Mol. Phylogenet. Evol.7094–98. 10.1016/j.ympev.2013.09.011

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 49

  KulakovskiyI. V.BoevaV. A.FavorovA. V.MakeevV. J. (2010). Deep and wide digging for binding motifs in ChIP-Seq data.Bioinformatics262622–2623. 10.1093/bioinformatics/btq488

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 50

  LartillotN.RodrigueN.StubbsD.RicherJ. (2013). PhyloBayes MPI: phylogenetic reconstruction with infinite mixtures of profiles in a parallel environment.Syst. Biol.62611–615. 10.1093/sysbio/syt022

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 51

  LaumerC. E.BekkoucheN.KerblA.GoetzF.NevesR. C.SorensenM. V.et al (2015). Spiralian phylogeny informs the evolution of microscopic lineages.Curr. Biol.252000–2006. 10.1016/j.cub.2015.06.068

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 52

  LavrovD. V.PettW. (2016). Animal mitochondrial DNA as we do not know it: mt-genome organization and evolution in nonbilaterian lineages.Genome Biol. Evol.82896–2913. 10.1093/gbe/evw195

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 53

  LuT. M.KandaM.SatohN.FuruyaH. (2017). The phylogenetic position of dicyemid mesozoans offers insights into spiralian evolution.Zool. Lett.3:6. 10.1186/s40851-017-0068-5

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 54

  MalakhovV. V. (1990). Mysterious Group of Marine Invertebrates: Tricholpax, Orthonectida, Dicyemida, Porifera.Moscow: Moscow University Press, 144.

  • Google Scholar
  • View reference in article

• 55

  MarlétazF.PeijnenburgK. T. C. A.GotoN.SatohN.RokhsarD. S. (2019). A new spiralian phylogeny places the enigmatic arrow worms among gnathiferans.Curr. Biol.29312–318. 10.1016/j.cub.2018.11.042

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 56

  MikhailovK. V.SlyusarevG. S.NikitinM. A.LogachevaM. D.PeninA. A.AleoshinV. V.et al (2016). The genome of Intoshia linei affirms orthonectids as highly simplified spiralians.Curr. Biol.261768–1774. 10.1016/j.cub.2016.05.007

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 57

  MinhB. Q.NguyenM. A.von HaeselerA. (2013). Ultrafast approximation for phylogenetic bootstrap.Mol. Biol. Evol.301188–1195. 10.1093/molbev/mst024

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 58

  MorozL. L.KohnA. B. (2013). Single-neuron transcriptome and methylome sequencing for epigenomic analysis of aging.Methods Mol. Biol.1048323–352. 10.1007/978-1-62703-556-9_21

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 59

  NesnidalM. P.HelmkampfM.MeyerA.WitekA.BruchhausI.EbersbergerI.et al (2013). New phylogenomic data support the monophyly of lophophorata and an ectoproct-phoronid clade and indicate that polyzoa and kryptrochozoa are caused by systematic bias.BMC Evol. Biol.13:253. 10.1186/1471-2148-13-253

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 60

  NguyenL. T.SchmidtH. A.von HaeselerA.MinhB. Q. (2015). IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies.Mol. Biol. Evol.32268–274. 10.1093/molbev/msu300

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 61

  NotoT.YazakiK.EndohH. (2003). Developmentally regulated extrachromosomal circular DNA formation in the mesozoan Dicyema japonicum.Chromosoma111359–368. 10.1007/s00412-002-0216-2

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 62

  NurkS.BankevichA.AntipovD.GurevichA. A.KorobeynikovA.LapidusA.et al (2013). Assembling single-cell genomes and mini-metagenomes from chimeric MDA products.J. Comput. Biol.20714–737. 10.1089/cmb.2013.0084

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 63

  OdintsovaM. S.YurinaN. P. (2005). Genomics and evolution of cellular organelles.Russ. J. Genet.41957–967. 10.1007/s11177-005-0187-5

  • CrossRef
  • Google Scholar
  • View reference in article

• 64

  OginoK.TsunekiK.FuruyaH. (2010). Unique genome of dicyemid mesozoan: highly shortened spliceosomal introns in conservative exon/intron structure.Gene44970–76. 10.1016/j.gene.2009.09.002

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 65

  O’LearyN. A.WrightM. W.BristerJ. R.CiufoS.HaddadD.McVeighR.et al (2016). Reference sequence (RefSeq) database at NCBI: current status, taxonomic expansion, and functional annotation.Nucleic Acids Res.44D733–D745. 10.1093/nar/gkv1189

  • Pubmed Abstract
  • CrossRef
  • Google Scholar
  • View reference in article

• 66

  PawlowskiJ.Montoya-BurgosJ. I.FahrniJ. F.WuestJ.ZaninettiL. (1996). Origin of the Mesozoa inferred from 18S rRNA gene sequences.Mol. Biol. Evol.131128–1132. 10.1093/oxfordjournals.molbev.a025675

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 67

  PetrovN. B.AleshinV. V.PegovaA. N.OfitserovM. V.SlyusarevG. S. (2010). New insight into the phylogeny of Mesozoa: evidence from the 18S and 28S rRNA genes.Moscow Univers. Biol. Sci. Bull.65167–169. 10.3103/S0096392510040127

  • CrossRef
  • Google Scholar

• 68

  PrinceF.PrinceF.KatsuyamaT.OshimaY.PlazaS.Resendez-PerezD.et al (2008). The YPWM motif links antennapedia to the basal transcriptional machinery.Development1351669–1679. 10.1242/dev.018028

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 69

  RobertsonH. E.SchifferP. H.TelfordM. J. (2018). The mitochondrial genomes of the mesozoans Intoshia linei, Dicyema sp. and Dicyema japonicum.Parasitol. Open4:e16. 10.1017/pao.2018

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 70

  Rodríguez-EzpeletaN.BrinkmannH.RoureB.LartillotN.LangB. F.PhilippeH. (2007). Detecting and overcoming systematic errors in genome-scale phylogenies.Syst. Biol.56389–399. 10.1080/10635150701397643

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 71

  RoureB.Rodriguez-EzpeletaN.PhilippeH. (2007). SCaFoS: a tool for selection, concatenation and fusion of sequences for phylogenomics.BMC Evol. Biol.7(Suppl. 1):S2. 10.1186/1471-2148-7-S1-S2

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 72

  RubanovL. I.SeliverstovA. V.ZverkovO. A.LyubetskyV. A. (2016). A method for identification of highly conserved elements and evolutionary analysis of superphylum Alveolata.BMC Bioinformatics17:385. 10.1186/s12859-016-1257-5

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 73

  RuppertE. E.FoxR. S.BarnesR. D. (2004). Invertebrate Zoology: A Functional Evolutionary Approach, Seventh Edn.Boston, MA: Brooks/Cole Thompson Learning.

  • Google Scholar

• 74

  SchifferP.RobertsonH.TelfordM. J. (2018). Orthonectids are highly degenerate annelid worms.Curr. Biol.281970.e3–1974.e3. 10.1016/j.cub.2018.04.088

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 75

  ShaginD. A.RebrikovD. V.KozhemyakoV. B.AltshulerI. M.ShcheglovA. S.ZhulidovP. A.et al (2002). A novel method for SNP detection using a new duplex-specific nuclease from crab hepatopancreas.Genome Res.121935–1942. 10.1101/gr.547002

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 76

  ShaoR.KirknessF.BarkerS. C. (2009). The single mitochondrial chromosome typical of animals has evolved into 18 minichromosomes in the human body louse, Pediculus humanus.Genome Res.19904–912. 10.1101/gr.083188.108

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 77

  ShimodairaH. (2002). An approximately unbiased test of phylogenetic tree selection.Syst. Biol.51492–508. 10.1080/10635150290069913

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 78

  ShimodairaH.HasegawaM. (2001). CONSEL: for assessing the confidence of phylogenetic tree selection.Bioinformatics171246–1247. 10.1093/bioinformatics/17.12.1246

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 79

  SlyusarevG. S. (2008). Phylum Orthonectida: morphology, biology, and relationships to other multicellular animals.Zh. Obshch. Biol.69403–427.

  • Pubmed Abstract
  • Google Scholar

• 80

  SlyusarevG. S.StarunovV. V. (2015). The structure of the muscular and nervous systems of the female Intoshia linei (Orthonectida).Org. Divers. Evol.1665–71. 10.1007/s13127-015-0246-2

  • CrossRef
  • Google Scholar

• 81

  SmithD. R.KeelingP. J. (2015). Mitochondrial and plastid genome architecture: reoccurring themes, but significant differences at the extremes.Proc. Natl. Acad. Sci. U.S.A.11210177–10184. 10.1073/pnas.1422049112

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 82

  StamatakisA. (2014). RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies.Bioinformatics301312–1313. 10.1093/bioinformatics/btu033

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 83

  StankeM.WaackS. (2003). Gene prediction with a hidden Markov model and a new intron submodel.Bioinformatics19(Suppl. 2), ii215–ii225. 10.1093/bioinformatics/btg1080

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 84

  SteinmetzP. R.KrausJ. E.LarrouxC.HammelJ. U.Amon-HassenzahlA.HoulistonE.et al (2012). Independent evolution of striated muscles in cnidarians and bilaterians.Nature487231–234. 10.1038/nature11180

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 85

  StruckT. H.Wey-FabriziusA. R.GolombekA.HeringL.WeigertA.BleidornC.et al (2014). Platyzoan paraphyly based on phylogenomic data supports a noncoelomate ancestry of spiralia.Mol. Biol. Evol.311833–1849. 10.1093/molbev/msu143

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 86

  SuskoE.RogerA. J. (2007). On reduced amino acid alphabets for phylogenetic inference.Mol. Biol. Evol.242139–2150. 10.1093/molbev/msm144

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 87

  SuzukiT. G.OginoK.TsunekiK.FuruyaH. (2010). Phylogenetic analysis of dicyemid mesozoans (phylum Dicyemida) from innexin amino acid sequences: dicyemids are not related to Platyhelminthes.J. Parasitol.96614–625. 10.1645/GE-2305.1

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 88

  TelfordM. J.HerniouE. A.RussellR. B.LittlewoodD. T. (2000). Changes in mitochondrial genetic codes as phylogenetic characters: two examples from the flatworms.Proc. Natl. Acad. Sci. U.S.A.9711359–11364. 10.1073/pnas.97.21.11359

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 89

  TsaiI. J.ZarowieckiM.HolroydN.GarciarrubioA.Sánchez-FloresA.BrooksK. L.et al (2013). The genomes of four tapeworm species reveal adaptations to parasitism.Nature49657–63. 10.1038/nature12031

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 90

  TuranoF. J.PantaG. R.AllardM. W.van BerkumP. (2001). The putative glutamate receptors from plants are related to two superfamilies of animal neurotransmitter receptors via distinct evolutionary mechanisms.Mol. Biol. Evol.181417–1420. 10.1093/oxfordjournals.molbev.a003926

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 91

  UniProt Consortium (2015). UniProt: a hub for protein information.Nucleic Acids Res.43D204–D212. 10.1093/nar/gku989

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 92

  Van BenedenE. (1876). Recherches sur les dicyémides, survivants actuels d’un embranchement des mésozoaires.Bull. Acad. Belg. Cl. Sci.411160–1205.

  • Google Scholar

• 93

  VortsepnevaE.TzetlinA.PurschkeG.MugueN.Haß-CordesE.ZhadanA. (2008). The parasitic polychaete known as Asetocalamyzas laonicola (Calamyzidae) is in fact the dwarf male of the spionid Scolelepis laonicola (comb. nov.).Invert. Biol.127403–416. 10.1111/j.1744-7410.2008.00137.x

  • CrossRef
  • Google Scholar

• 94

  WatanabeK. I.BesshoY.KawasakiM.HoriH. (1999). Mitochondrial genes are found on minicircle DNA molecules in the mesozoan animal Dicyema.J. Mol. Biol.286645–650. 10.1006/jmbi.1998.2523

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 95

  WaterhouseR. M.SeppeyM.SimãoF. A.ManniM.IoannidisP.KlioutchnikovG.et al (2017). BUSCO applications from quality assessments to gene prediction and phylogenomics.Mol. Biol. Evol.10.1093/molbev/msx319 [Epub ahead of print].

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 96

  WeigertA.HelmC.MeyerM.NickelB.ArendtD.HausdorfB.et al (2014). Illuminating the base of the annelid tree using transcriptomics.Mol. Biol. Evol.311391–1401. 10.1093/molbev/msu080

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 97

  WorsaaeK.RouseG. W. (2010). The simplicity of males: dwarf males of four species of Osedax (Siboglinidae, Annelida) investigated by confocal laser scanning microscopy.J. Morphol.271127–142. 10.1002/jmor.10786

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 98

  ZamanianM.KimberM. J.McVeighP.CarlsonS. A.MauleA. G.DayT. A. (2011). The repertoire of G protein-coupled receptors in the human parasite Schistosoma mansoni and the model organism Schmidtea mediterranea.BMC Genomics12:596. 10.1186/1471-2164-12-596

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 99

  ZdobnovE. M.TegenfeldtF.KuznetsovD.WaterhouseR. M.SimãoF. A.IoannidisP.et al (2017). OrthoDB v9.1: cataloging evolutionary and functional annotations for animal, fungal, plant, archaeal, bacterial and viral orthologs.Nucleic Acids Res.45D744–D749. 10.1093/nar/gkw1119

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 100

  ZhongM.HansenB.NesnidalM.GolombekA.HalanychK. M.StruckT. H. (2011). Detecting the symplesiomorphy trap: a multigene phylogenetic analysis of terebelliform annelids.BMC Evol. Biol.11:369. 10.1186/1471-2148-11-369

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 101

  ZhongY. F.HollandP. W. (2011). HomeoDB2: functional expansion of a comparative homeobox gene database for evolutionary developmental biology.Evol. Dev.13567–568. 10.1111/j.1525-142X.2011.00513.x

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 102

  ZhuY. Y.MachlederE. M.ChenchikA.LiR.SiebertP. D. (2001). Reverse transcriptase template switching: a SMART approach for full-length cDNA library construction.Biotechniques30892–897. 10.2144/01304pf02

  • Pubmed Abstract
  • CrossRef
  • Google Scholar

• 103

  ZhulidovP. A.BogdanovaE. A.ShcheglovA. S.VagnerL. L.KhaspekovG. L.Kozhemyakoet al (2004). Simple cDNA normalization using kamchatka crab duplex-specific nuclease.Nucleic Acids Res.32:e37. 10.1093/nar/gnh031

  • Pubmed Abstract
  • CrossRef
  • Google Scholar