Substantial Hierarchical Reductions of Genetic and Morphological Traits in the Evolution of Rotiferan Parasites

Abstract

Within the last 800 million years, animals evolved a vast range of diversity of species exhibiting an enormous disparity of forms and lifestyles. The process involved an increase in complexity from life forms with few cell types to organisms with many hundreds of cell types. However, neither genome size nor number of protein-coding genes can explain these differences, and their biological basis remains elusive. Yet, recent studies suggest that the evolution of complexity is closely linked to the acquisition of a class of noncoding gene regulators called microRNAs. To test this hypothesis, we investigated the association between loss of organismal complexity and microRNAs in Syndermata, an invertebrate group including free-living wheel animals (Monogononta, Bdelloidea), epibiotic Seisonidea, and endoparasitic thorny-headed worms (Acanthocephala). Analyses of genomic, transcriptomic, and morphological data of altogether 25 syndermatan species revealed strong correlations of microRNA losses with reductions of protein-coding genes and morphological traits. The hierarchical pattern sums up to ∼85% loss of microRNAs and a ∼50% loss of conserved metazoan core genes (Benchmarking Universal Single-Copy Orthologs) on the lineage to thorny-headed worms. Extraordinarily reduced microRNA complements were confirmed by small RNA sequencing data. Endoparasitic Acanthocephala was additionally distinguished by the most morphological reductions of ancestral features, such as the digestive tract. Together, we observed that reductions of ∼400 protein-coding genes and 10 metazoan core genes tended to accompany the loss of single microRNA families. Furthermore, 4 microRNA families and 34 metazoan core genes appeared to be associated, on average, with the reduction of a single morphological trait.

Keywords: Rotifera; complexity; core genes; evolution; microRNAs; piRNAs.

Graphical Abstract

Fig. 1.

MicroRNA complements of Syndermata reflect stepwise losses from free-living to parasitic ancestors and unprecedented losses of 85% of conserved microRNAs in acanthocephalans. a) Banner-plot of 25 microRNA complements of included syndermatan species and two outgroups. White field means that no microRNA was found. Heatmap function refers to the paralogue number in each microRNA family (log e). Note that both outgroups do not belong to Syndermata and, hence, have none of the Rotifera microRNA families. Furthermore, S. kowalevskii is not a protostome and hence lacks protostome microRNA families. See supplementary file S1, Supplementary Material online “species microRNA families”. Arrows depict availability of small RNA sequencing data (pink, publicly available; green, novel). b) Schematic trees of the four syndermatan groups representing loss of usually conserved and gain of novel microRNA families. Branch lengths correspond to the number of gains and losses. c) Overview of microRNA precursor lengths in selected representatives of monogononts (B. plicatilis), bdelloids (A. vaga), the seisonid S. nebaliae, and the acanthocephalan P. laevis. d) Selected MIR-1 examples of the same species as in c) highlighting the length deviations in S. nebaliae. e) Alignment of the mature sequence of all MIR-1 genes in syndermatans (N = 54) highlights its very strong conservation in the pattern characterized for bilaterian mature miRNAs (Wheeler et al. 2009; Fromm et al. 2015) with strong conservation in the seed and 3′ CR and higher variability in positions 9 to 12 and 17 to 22. Species order corresponds to a); consensus is given in black.

Fig. 2.

Numbers of protein-coding genes and protein-coding core genes (BUSCO), but not genome size or assembly quality (N50), correlate with microRNA loss. a) Genome size, N50, number of protein-coding genes, and metazoan BUSCO completeness of 27 syndermatan species and two outgroup representatives (see supplementary file S1, Supplementary Material online “Syndermata genomes”). b) Correlation of protein-coding gene numbers with cumulative numbers of microRNA losses in representative syndermatans (B. koreanus, A. vaga, S. nebaliae, P. laevis) and two outgroup representatives (C. teleta, S. kowalevski). Gray area corresponds to 95% confidence interval. The asterisk (*) refers to a high number of paralogues in bdelloids due to genome duplication (see supplementary file S1, Supplementary Material online “correlations”). c) BUSCO completeness highlights duplicated genomes in bdelloids (see supplementary file S1, Supplementary Material online “Busco_syndermatan”). d) Correlation of metazoan BUSCO core genes with the loss of microRNA genes in the same representative syndermatans as in b) and tree giving BUSCO gene losses. Gray area in graph corresponds to 95% confidence interval. Branch lengths correspond to the number of losses.

Fig. 3.

Missing BUSCO genes in Pararotatoria are highly enriched for biological regulation GO terms. a) Venn diagram illustrating the overlap of missing metazoan BUSCO genes in the Pararotaria representatives S. nebaliae (Seisonidea) and Neoechinorhynchus agilis (Acanthocephala). Of the 954 protein-coding genes usually conserved in Metazoa, the seisonid and acanthocephalan collectively lack 280 genes. Additional 76 and 148 genes are absent in either the seisonidean or the acanthocephalan, respectively. b) Significantly enriched GO terms in the metazoan BUSCO genes missing in Pararotatoria. Enrichment of three Molecular Function GOs and one Cellular Component GO contrast with 81 enriched Biological Process GO terms. Note the high number of terms referring to regulation (pink arrows). c) Significantly enriched Biological Process GOs for missing genes in pararotatorians, whereby each dot represents one term. Color and weight of the dots correspond to the scale below the plot. Note the high number of regulatory terms (pink arrows).

Fig. 4.

a) Syndermatan character matrix combining novel characters with the ones compiled by Deline et al. (2018). Dark gray versus white in taxon columns indicates character presence versus absence. Column 2 indicates the source of respective characters: A—Deline et al. (2018), B—modified wording of character definition by Deline et al. (2018), C—Ahlrichs (1995), D—assessment by authors considering: ¹Ahlrichs (1997), ²Fontaneto and De Smet (2015), ³Nicholas and Mercer (1965), ⁴Herlyn (2021), ⁵Ahlrichs (1995), ⁶Amin (1987), ⁷Herlyn and Taraschewski (2017). b) Schematic trees of the four syndermatan groups representing loss of plesiomorphic (purple, top) and gain of apomorphic (green, bottom) morphological features. Branch lengths correspond to gains and losses. c) Morphological character loss (purple) and gain (green) in the extreme example of Acanthocephala. The eponymous attachment organ (proboscis) is inverted in the juvenile female and everted in the adult male. Drawings refer to Neoechinorhynchus rutili, a close relative of N. agilis. Modified from Steinsträsser (1936). d) Dependence of plesiomorphic morphological characters on loss of microRNA genes in representative syndermatans and two outgroup representatives (same species as in Fig. 2). e) Interrelations of metazoan BUSCO core genes and loss of plesiomorphic morphological characters in representative syndermatans and two outgroup representatives. Outgroup species are the same as in Fig. 2 and were set to zero morphological losses. The last common ancestor of extant Acanthocephala might still have possessed protonephridia. In any case, the two acanthocephalans studied (P. laevis and N. agilis) lack protonephridia.

Fig. 5.

piRNA retention in syndermatans. a) Length distribution of genome-mapped small noncoding RNAs, with putative piRNAs highlighted. b) Distribution of overlap between putative piRNAs on opposite strands, highlighting a peak at 10nt where present, indicative of ping-pong biogenesis. c) Cumulative frequency plot showing the tendency of putative piRNAs to be clustered together along the genome, suggestive of potential piRNA-generating loci or piRNA clusters. d) Conservation of proteins related to small noncoding RNA pathways. Annotation of the presence or absence of homologue genes based on best reciprocal blast searches is illustrated for each species accordingly (dark blue, multiple hits; light blue, one hit). e) The piRNAs/2 kb windows for each of the 2 kb windows with at least 1,000 piRNAs mapped. f) The clustering ratio, defined as the ratio of the fraction of the genome containing 90% of randomly shuffled piRNAs compared to the fraction containing 90% of true piRNA positions. A higher clustering ratio thus indicates that a smaller fraction of the genome contains 90% of the piRNAs than expected by chance. Avag, A. vaga; Bkor, B. koreanus; Nag, N. agilis; Plae, P. laevis; Sneb, S. nebaliae.