18S rRNA variability maps reveal three highly divergent, conserved motifs within Rotifera
BMC Ecology and Evolution volume 21, Article number: 118 (2021)
18S rRNA is a major component of the small subunit of the eukaryotic ribosome and an important phylogenetic marker for many groups, often to the point of being the only marker available for some. A core structure across eukaryotes exists for this molecule that can help to inform about its evolution in different groups. Using an alignment of 18S rDNA for Rotifera as traditionally recognized (=Bdelloidea, Monogononta, and Seisonacea, but not Acanthocephala), I fitted sequences for three exemplar species (Adineta vaga, Brachionus plicatilis, and Seison nebaliae, respectively) to the core structure and used these maps to reveal patterns of evolution for the remainder of this diverse group of microscopic animals.
The obtained variability maps of the 18S rRNA molecule revealed a pattern of high diversity among the three major rotifer clades coupled with strong conservation within each of bdelloids and monogononts. A majority of individual sites (ca. 60%) were constant even across rotifers as a whole with variable sites showing only intermediate rates of evolution. Although the three structural maps each showed good agreement with the inferred core structure for eukaryotic 18S rRNA and so were highly similar to one another at the secondary and tertiary levels, the overall pattern is of three highly distinct, but conserved motifs within the group at the primary sequence level. A novel finding was that of a variably expressed deletion at the 3' end of the V3 hypervariable region among some bdelloid species that occasionally extended into and included the pseudoknot structure following this region as well as the central “square” of the 18S rRNA molecule. Compared to other groups, levels of variation and rates of evolution for 18S rRNA in Rotifera roughly matched those for Gastropoda and Acanthocephala, despite increasing evidence for the latter being a clade within Rotifera.
The lack of comparative data for comparable groups makes interpretation of the results (i.e., very low variation within each of the three major rotifer clades, but high variation between them) and their potential novelty difficult. However, these findings in combination with the high morphological diversity within rotifers potentially help to explain why no clear consensus has been reached to date with regard to the phylogenetic relationships among the major groups.
Together with numerous ribosomal proteins, 18S rRNA forms a major component of the small subunit of the eukaryotic ribosome. The single stranded RNA molecule itself has a characteristic and complicated secondary structure [see 1, 2], whereby it repeatedly folds back upon itself, with the resultant base-pairing or lack thereof creating stems and loops, respectively. Two or more stem-loop regions can also combine to form one of 14 different types of three-dimensional pseudoknot , thereby contributing to the tertiary structure of the molecule. In addition, eukaryotes share nine homologous regions in the molecule (and which are also present in the prokaryotic homologue 16S rRNA) that are especially variable and are labelled as the hypervariable regions V1—V9 . (The eukaryotic V6 region, however, is noticeably less variable compared to the other regions  and to such a degree that it is not counted among the hypervariable regions by some authors (e.g., )). The molecule itself is encoded by the 18S rDNA gene, a distinction that I will maintain throughout this paper although the terms rRNA are rDNA are often used interchangeably in the literature.
The phylogenetic utility of 18S rDNA arguably stems in part historically from practical considerations. Together with 5.8S rDNA, 28S rDNA, and two internal and two external transcribed spacers, it forms an array that is tandemly repeated throughout the eukaryotic genome (e.g., ca. 300 + copies clustered across five chromosomes in humans ). The multicopy nature of the gene made it easier to extract and amplify in the pre-PCR era and concerted evolution meant that the many copies are virtually identical [6, 7], thereby sidestepping questions of paralogy. In addition, 18S rDNA as a gene is found universally among eukaryotes and possesses conserved flanking regions that facilitated primer design. However, its practical utility was augmented by the broad phylogenetic information content yielded by its structural characteristics, with the slower, more conserved stem regions (due to the constraints of the base pairing) providing resolution deeper in the tree to compliment the more recent information provided by the faster, less conserved loops and especially by the hypervariable regions. Indeed, much of the basis for deep phylogenies within Metazoa and beyond derive from phylogenetic analyses of 18S rDNA or other rDNA molecules more generally, both at the sequence and, more recently, the meta-sequence (i.e., structural) levels (e.g., [8,9,10]). At the other end of the spectrum, certain hypervariable regions, or parts thereof, have been promoted as possible species barcoding regions in diatoms (e.g., [11, 12]) and other “protists” (e.g., ). Although the early promise of 18S rDNA as “the” phylogenetic marker has not been realized, it remains one of the most widely sequenced genes across all organismal groups, especially in a phylogenetic context .
Indeed, 18S rDNA is often one of only a few, if not the only, phylogenetic marker sequenced for a given group. A case in point is Rotifera, a historically recognized phylum of approximately 2000 named species of microscopic animals  for which the only comprehensive molecular phylogeny to date encompasses 53 species (plus numerous outgroup species) sequenced for up to four markers including 18S rDNA as part of a total-evidence analysis with morphological characters . Although countless studies based on either morphological or molecular data confirm the monophyly of the three major rotifer clades–Bdelloidea (ca. 388 species), Monogononta (ca. 1623 species), and Seisonacea (four species; more commonly referred to as Seisonidae)–their relationships to one another remain unclear [see 15], in part because of their highly distinct natures. Bdelloids are obligate asexuals (i.e., only parthenogenetic females are known) that can also undergo anhydrobiosis and whose genome has been shaped in part by (ancient) gene exchange. Monogononts, by contrast, are facultative asexuals that undergo cyclic parthenogenesis, possess only a single gonad and comprise many species presenting dwarf males to various degrees. Finally, seisonids are obligate sexuals with no male dwarfism and that live as ectoparasites or commensally on different species of the crustacean genus Nebalia, sometimes with different seisonid species living on different body parts of the same host [17, 18]. Even the application of molecular phylogenetics has failed to firmly resolve the relationships of these taxa to one another .
My goal in this paper is not so much to resolve the phylogenetic relationships within Rotifera, nor to address the question as to whether Acanthocephala nest within it (see “Methods”), but rather to examine rates of molecular evolution of the 18S rDNA gene within this traditional grouping based on newly obtained sequences combined with those present in GenBank. In so doing, however, my analyses reveal a pattern of highly distinct, yet conserved motifs among the major clades that might explain why there has been so little consensus about their interrelationships to date [see 15].
Examination of the uncorrected, average pairwise distances across groups determined using PAUP* 4.0a166  revealed a pattern (Table 1A) whereby each of the taxonomic groups Bdelloidea, Monogononta, Seisonacea were all highly distinct from one another as well as from the platyhelminth outgroup Calicophoron calicophorum (> 14% divergence), but showed little internal variation (< 3%). This pattern is also reflected indirectly in the increase in the percentage of gaps in the Rotifera alignment compared to either of the Bdelloidea or Monogononta alignments (Table 2). In addition, sequences from monogonont species were, on average, slightly more similar to those of the outgroup than they were to the remaining rotifer clades. The use of corrected distances calculated using a GTR model of evolution revealed the same pattern (Table 1B), with average pairwise distances being on a par with (within clade comparison) or slightly higher than (between clade comparisons) the uncorrected distances.
The inferred structural maps of the three exemplar rotifer species—Adineta vaga (Bdelloidea), Brachionus plicatilis (Monogononta), and Seison nebaliae (Seisonacea) (Fig. 1; see also Additional file 1)—each showed a good fit to the eukaryotic core structure for 18S rRNA proposed by Van de Peer and colleagues [20,21,22]. Unusual, however, is that most bdelloids (23 of 29 fully informative species; potentially the incompletely sequenced species Otostephanos jolantae and Zelinkiella synaptae as well) possess a unique and variably expressed deletion that spans 92 bps comprising the last 14 bps of region V3 and extends into the non-hypervariable region beyond this (Fig. 2). Six motifs of different lengths are present including the absence of the deletion and no species expresses the deletion in its full length (maximum deletion length is 68 bp). Although the deletion begins with a set of four nucleotides displaying very high relative rates of evolution (Fig. 1A), the patterns of the motifs suggest that the deletion originates from its 3′ end in a non-hypervariable region of the 18S rRNA molecule that is otherwise virtually constant in its sequence composition across Rotifera and encompasses most, if not all, of the pseudoknot following the V3 region in the eukaryotic core structure. The deletion is also variably expressed among species within each of the genera Adineta, Embata, Mnobia, Philodina and Rotaria (but not Abrochtha, Dissotrocha, or Habrotrocha) as well as in the species Dissotrocha aculeata, Dissotrocha macrostyla, Philodina citrina, and Philodina megalotrocha, with second GenBank sequences for each of these species (accession numbers JX494743, JX494745, JX494740, and JX494741, respectively) possessing even longer deletions. Altogether, this variability at both the genus and species levels together with the apparent presence of a more restricted deletion in the monogonont species Lindia tecusa and Lindia torulosa (Fig. 2) would suggest the convergent evolution of the deletion motifs barring any sequencing or identification errors for these GenBank sequences.
Average relative TIGER values across all sites were high, indicating slow rates of evolution (Table 3, Fig. 3), largely because > 60% of all sites (i.e., > 1000 bp) with sufficient coverage for the TIGER analyses within a given data set (Bdelloidea, Monogononta, and all Rotifera) were constant. The more homogenous bdelloid and monogonont data sets showed even greater proportions of constant sites (85.5% and 75.5%, respectively) and each also presented slower average relative rates across sites than did the entire rotifer data set. However, even the variable sites alone, which would include the hypervariable regions, were not unduly fast, with average TIGER values for Bdelloidea and Monogononta being slightly less than 0.5 and only decreasing to around 0.4 for all Rotifera (Table 3). Indeed, the rates for variable sites tended to cluster around these values such that variable sites that evolved extremely slowly (TIGER rate > 0.6) were all but absent as were those that evolved very rapidly (TIGER rate < 0.4) with the possible exception of across Rotifera as a whole.
There was no difference in the inferred relative rates of evolution between paired sites in each of the three data sets (column 2 in Table 4) and together these sites (“stems”) evolved slightly slower than unpaired sites (“loops”) for Bdelloidea only (column 3 in Table 4). Highly significant differences in relative rates were present between regions inferred as being within hypervariable versus non-hypervariable regions, regardless whether the former regions were pooled for the analysis or analysed individually (columns 4 to 6, respectively, in Table 4). The average relative rates for most hypervariable regions were faster than those for the pooled non-hypervariable regions (Fig. 4), with V6 being the notable, expected exception in addition to several hypervariable regions for bdelloids. In addition, the rates in the Rotifera data set were generally faster than those in either the bdelloid or monogonont data sets, reflecting in part the sequence differences between the motifs in the latter two. For the latter data sets, bdelloids generally presented slower rates than did monogononts with the exception of hypervariable region V3.
Substitution hotspots were present along the entire molecule (Fig. 5), but more common and with higher relative rates in the hypervariable regions (with the exception of V6). However, most of these hotspots did not express unduly high relative rates of evolution. Relative rates within the hypervariable regions could also differ as exemplified by region V4, where the 5′ half of the region displays noticeably higher rates than the 3′ end where two pseudoknots are located. Again, rates of evolution at any given position were generally the highest for Rotifera and lowest for Bdelloidea among the three data sets.
Indels were numerous in each of the three exemplar species when aligned against one another in the Rotifera data set, and disproportionately so in the hypervariable regions, but were very short on average (usually < 2 bp long; Table 5). Only hypervariable region V7 showed a tendency toward having longer indels despite it being noticeably shorter in rotifers than in the outgroup species C. calicophorum (as well as in Daphnia pulex, and Loricera foveata; see Methods). Region V6 is noteworthy insofar as no indels were inferred between the three major rotifer clades (Table 5).
Altogether, my examination of the evolution of 18S rDNA in rotifers revealed a clear pattern whereby the high morphological disparity among the three major clades is matched by their molecular disparity for this molecule. The high number of very short indels together with the extremely restricted sequence variation within each of the two largest clades (Bdelloidea and Monogononta) also indicates that the disparity derives mostly from substitutions, and then predominantly in the hypervariable regions. However, despite their name, even the hypervariable regions are largely conserved within the three major rotifer clades. By way of comparison, as well as to underscore the degree of sequence conservation within the major clades, the analogous average uncorrected pairwise distance for an alignment of 69 18S rDNA GenBank sequences for Acanthocephala (see Additional file 2) that includes all its four major subgroups within the clade is 13.1%. Although this value approaches that of the average pairwise divergence among the three rotifer 18S rDNA motifs (from 17.0 to 22.2% from Table 1 or 18.7% ± 1.5% among the three exemplar species only), it is an order of magnitude higher than that for within each of the major rotifer clades (< 3%; compare Table 1), a taxonomic status that arguably also applies for Acanthocephala (see Methods). Similarly, TIGER rates of evolution for Acanthocephala (Table 3) are comparable to those for Rotifera as a whole, rather than to those for either Bdelloidea or Monogononta.
This extreme pattern in Rotifera is also easily visualized through a maximum-likelihood phylogeny of the sequences examined in this study together with those of Acanthocephala for context (Fig. 6; Additional file 3). Here it is the pattern that is of chief interest rather than the relationships per se, with Bdelloidea and Seisonacea forming sister taxa at the ends of extended branches to the exclusion of Monogononta, which shows much less divergence from the rotifer common ancestor. Monogononta also displays extremely reduced molecular divergence between its members, despite being by far the best represented of the three major clades. Finally, Acanthocephala again display as much internal molecular divergence in the phylogeny as for all true rotifers.
Moreover, the variability among the rotifer clades is restricted to at most 40% of the 18S rRNA molecule, with the remaining sites being constant across all Rotifera. This value coincidentally matches the proportion of the molecule that comprises hypervariable regions. However, even though the hypervariable regions are indeed significantly more variable than the non-hypervariable ones (both in the number of variable sites and their rates; Table 6), variable sites are found in both regions. Although the species sampling in this study was necessarily restricted, the taxonomic diversity of Rotifera was well represented insofar as exemplars from all major taxonomic subdivisions within each of Rotifera, Bdelloidea, and Monogononta (see Methods) were present in the data set. It could be argued that the number of constant sites is slightly overestimated insofar as gaps between the major clades were often preferred to substitutions. However, the increase in length compared to length of the entire molecular is negligible (about 10%; see Table 2) and this problem would apply chiefly to the entire rotifer data set. In addition, the high proportion of constant sites across Rotifera matches that inferred across Gastropoda by Weigand et al.  and for Acanthocephala (Table 3). Thus, the general patterns observed here, especially the highly distinct motifs for the major clades, are likely to be accurate. In addition, preliminary results from 28S rDNA and MT-CO1 confirm this general pattern of highly distinct motifs (results not shown), although not to the same extent as seen for 18S rDNA.
The lack of comparative data for other, comparable clades prohibits a good assessment of the novelty of the observed patterns across eukaryotes and the loss of readily available structural data through the demise of the European Ribosomal Database  is strongly felt in this regard. As mentioned, the high proportion of constant sites is comparable to that seen across Gastropoda, although the latter show a more even distribution of rates across the entire spectrum from slow to fast than was the case for Rotifera, albeit with a strong spike at intermediate rates . Again, there is also a strong similarity in many parameters between Rotifera as a whole and Acanthocephala. Although the variability map of Van de Peer et al.  shows comparatively few constant sites and many very fast ones, it is at the level of all Eukaryota and so hardly comparable in terms of taxonomic breadth. The latter, however, does indicate that the pseudoknot following the V3 region is relatively conserved across eukaryotes, making its variable deletion within Bdelloidea all that more unusual. Even more extraordinary is that the deletion in both Rotaria neptunoida and Rotaria tardigrada includes the central “square” of the 18S rRNA molecule from which its three main arms originate and thereby could affect the tertiary structure of the entire molecule in some unknown and potentially severe fashion.
Regardless of the potential patterns in other taxonomic groups, one explanation for the pattern observed here is that the three rotifer crown groups are of relatively recent origin and, with the exception of Seisonacea, have undergone rapid adaptive radiations. If true, this scenario would imply a high degree of morphological plasticity in Monogononta in particular given both the higher number of species as well as morphological diversity across the group compared to Bdelloidea. What remains unclear, however, is why monogonont sequences, on average, remain more similar to that of a relatively distant platyhelminth outgroup (Table 1) instead of to the remaining rotifers or why this clade as a whole also does not subtend an extended branch like Bdelloidea and Seisonacea (and even Acanthocephala) as would be expected for a recent radiation within an otherwise ancient group. Of particular interest in this general context would be the sequencing of the remaining seisonid species. Given the lack of any obvious widespread dispersal abilities in seisonids, perhaps in concert with the hypothesis that they are among the oldest of the rotifer clades , the apparently exclusive association between them and their Nebalia hosts could be ancient, which agrees with the extended branch leading to S. nebaliae. Less clear, however, is whether the individual associations are also ancient or of more recent origin, especially given that different seisonid species can be found on the same host species if not the same host individual [17, 18]. Additionally, or alternatively, it could be that the three clades have independently reduced their rates of molecular evolution. However, it is not clear what the mechanism behind this would be and why the same process has not occurred in Acanthocephala, especially given that this taxon does indeed appear to nest within Rotifera [also 15, 16, 26–28]. Possible explanations for the latter discrepancy could lie with the endoparasitic lifestyle of all acanthocephalans as compared to the free-living true rotifers (with the exception of Seisonacea) or that the crown group is simply older and so shows more within-group molecular diversity.
Problematic in testing these hypotheses is that divergence time estimates for and within Rotifera are all but absent. Apart from a few reports of subfossilized rotifers from Holocene peat deposits (e.g., [29, 30]), the only other known rotifer fossils comprise contracted bdelloid specimens or their theca encased in amber, the oldest pieces of which have been dated to 35–40 Ma ago [31, 32]. Although the theca in particular have been assigned to the extant genus Habrotrocha , the contracted nature of the specimens and the simplistic features of the theca make it difficult to determine their species identity precisely and thus caution is perhaps advised as to whether or not either set of specimens belong to crown group Bdelloidea. Molecular based studies place the divergence between Rotifera and Platyhelminthes from anywhere between 492 to 1160 Ma ago (best estimate, 824 Ma ago; www.timetree.org ), with the origin of Rotifera being necessarily more recent than this, especially given that Platyhelminthes are likely not the immediate sister group of Rotifera, even among extant taxa [34, 35]. As such, it is unknown how old Rotifera are as a group as well as what the ages of its three major crown groups are, information that is needed to better understand the pattern of evolution of 18S rDNA witnessed in this paper, and possibly of other markers as well.
Nevertheless, the observed pattern potentially explains the severe problems in reconstructing the phylogenetic history of Rotifera at higher taxonomic levels. As summarized by Fontaneto and de Smet , about the only points of consensus in this context are the monophyly of each of the three major rotifer clades. Indeed, morphological and molecular analyses paint different pictures of higher-level rotifer phylogeny. Whereas the former tend to support Bdelloidea and Monogononta as sister taxa (= Eurotatoria) within a monophyletic Rotifera, the latter usually cluster Bdelloidea and Seisonacea together (as herein) within a paraphyletic Rotifera because of the inclusion of Acanthocephala, sometimes as the immediate sister group to Seisonacea. However, the pattern of evolution in 18S rDNA presented in this paper suggests that traditional sequence-based phylogenetic analyses of Rotifera could potentially be compromised by artefacts known to arise from long-branch attraction [see 36], which could be especially severe in this case given the lengths of the branches involved. In fact, the major rotifer clades are so distinct from one another molecularly that even analyses of MT-CO1 support the monophyly of each of Bdelloidea and Monogononta with high bootstrap support (results not shown), although MT-CO1 (or at least that part obtained using the Folmer  primers) normally loses phylogenetic signal above the genus level so rapidly that its use is typically restricted to species barcoding . Instead, inferences based on rare genomic changes [see 39] or other types of markers, including meta-sequence features [see 40] might be more reliable for unravelling higher-level relationships within Rotifera. Interestingly, results from an analysis of mitochondrial gene order  do match those from traditional sequence analyses, despite the limited taxon sampling in that study as well as the methodological problems that are known for analyses of gene order and that have prevented these data from playing a more prominent role in phylogenetic analyses to date [see 41, 42].
Fortunately, any potential artefacts caused by long-branch attraction appear to be limited to the inferences of the relationships among the three major clades, rather than the relationships within each of them where the branches are shorter (see Fig. 6). At these less inclusive levels, the inferred relationships in the 18S rDNA tree (Additional file 3) tend to reflect accepted taxonomic groups down to the genus level within Rotifera, especially within Monogononta (see The Rotifera World Catalog (www.rotifera.hausdernatur.at), despite the high proportion of constant sites. As such, analogous to the case with missing data (see ), it would appear that the limited number of variable sites are both sufficient in number and do not evolve unduly rapidly to provide good resolution. Nevertheless, the rate variation among sites that is present requires that some correction for rate heterogeneity is used [see also 20, 23] and even this might be insufficient at higher taxonomic levels within Rotifera.
Largely complete 18S rDNA sequences (one per species; see Table S1, Additional file 2) were either compiled from GenBank or newly generated within the working group (87 sequences, all from Monogononta), in part for other studies (e.g., [44, 45]). All previously unpublished sequences have been deposited in GenBank under the accession numbers MT522624–MT522695 and MT542324. The newly derived sequences were obtained following the protocols outlined in Kimpel et al.  and Wilke et al. . All taxonomic names were verified against The Rotifera World Catalog (accessed on April 10, 2020). In total, 198 rotifer sequences were used, including 162 monogonont sequences, 35 bdelloid sequences, and one seisonid sequence. Although this amounts to only roughly 10% of the described species diversity, all three major rotifer clades (Bdelloidea, Monogononta, and Seisonacea) were represented as were both major clades within Monogononta (Gnesiotrocha and Pseudotrocha; ) and all three within Bdelloidea (Adinetida, Philodinavida, and Philodinida; ); Seisonacea comprises only a single clade of four described species in two genera . Although there is mounting molecular evidence that Acanthocephala (ca. 1330 species) nests within Rotifera, possibly as the sister taxon to Seisonacea (e.g., [16, 26,27,28]), I have restricted my primary analyses to “true” rotifers as traditionally recognized only. The complete 18S rRNA sequence for the flatworm Calicophoron calicophorum (Platyhelminthes: Trematoda: Digenea; GenBank accession L06566) was added to the data set as the most closely related outgroup to Rotifera for which a structural map was available (from the now defunct European Ribosomal Database; bioinformatics.psb.ugent.be/webtools/rRNA/secmodel/Ccal_SSU.html; ).
The data set was aligned initially using the default settings of MUSCLE v3.8.31  as implemented in SeaView v5.0  and then improved by eye to correct for obvious errors with reference to the structural map for C. calicophorum whenever possible. Automated alignment programs that account for secondary structure, including LocARNA , MAFFT , or RNAsalsa , could not be used because they cannot accommodate the pseudoknots present within the eukaryotic 18S rRNA tertiary structure (see Analyses below). The alignment, however, was relatively trivial insofar as areas of disagreement were usually localized to several discrete and obvious indels of various sizes, usually in the hypervariable regions. From the complete data set, two subsets comprising each of bdelloids and monogononts only were also constructed. The aligned length of the final data set (available as Additional file 4) and subsets of it (Bdelloidea, Monogononta, and all Rotifera including Seisonacea) can be found in Table 2.
Fitting the rotifer 18S rRNA molecules to the core structure of the eukaryotic 18S rRNA molecule proposed by Van de Peer and colleagues [20,21,22] is difficult because the latter contains three pseudoknots—one immediately following region V3 and two within region V4 [1, 2, 53]—as well as some other nested base pairings, each of which represent computationally hard problems [54, 55]. Although numerous programs exist that can compute pseudoknots or can fit a sequence to a given (partially resolved) structure based on minimizing free energies, none exist to my knowledge that can do both. In addition, it is likely unwise to attempt to constrain the structures of the hypervariable regions based on known structures from distantly related taxa as is the case here with the platyhelminth outgroup. Thus, to maintain the highest degree of accuracy as well as to incorporate as much information from the core structure as possible, I employed a two-step, ad hoc procedure to transfer the eukaryotic core structure to three representative rotifer rRNA sequences.
First, I obtained a conservative model of the core structure by aligning the complete 18S rRNA sequences for C. calicophorum, the water flea Daphnia pulex (Crustacea: Branchiopoda: Cladocera; GenBank accession AF014011), and the ground beetle Loricera foveata (Insecta: Coleoptera: Carabidae; AF012503), which together represent the most closely related outgroups to Rotifera for which structural maps are available from the European Ribosomal Database. The initial alignment again used MUSCLE, with the subsequent, manual refinement of it focussing exclusively on the more conserved, non-hypervariable regions with the aid of the structural maps. Again, for these relatively constant regions, the alignment procedure was relatively trivial. In the non-hypervariable regions, I used a conservative approach insofar as only those positions that were consistently paired versus unpaired across these three relatively distantly related reference sequences were constrained as such for the structural analyses of the rotifer sequences (denoted using paired parentheses and an x, respectively). The structure of all remaining positions was unconstrained (denoted using a period). All hypervariable regions as well as all pseudoknots and other nested pairings were constrained to be unpaired, thereby forming unresolved, extended loops initially. This generalized structural constraint was then applied to each of three rotifer exemplar species—Adineta vaga (for Bdelliodea), Brachionus plicatilis (for Monogononta), and Seison nebaliae (for all Rotifera)—using RNAfold v2.4.11  to obtain the basic core structure. However, the constraint could be overruled by the program insofar as paired positions were not strictly enforced for these analyses (i.e., the switch—enforceConstraint was not set). Nucleotide positions unique to the rotifer sequences with respect to the core backbone (i.e., insertions) were left unconstrained (non-hypervariable regions) or constrained to be unpaired (hypervariable regions).
Second, the structures of the hypervariable regions in each exemplar rotifer species were resolved individually using RNAfold in the absence of any constraints apart from “bracketing” each region so that it began with a stem of at least two paired positions to provide some basal, structural context. Where possible, the brackets were obtained from the real sequences immediately adjacent to the hypervariable region according to the core backbone (regions V2, V4, V5, and V9). For regions V1, V6, V7, and V8, however, the brackets used obtained from the start of the hypervariable region, whereas an artificial bracket needed to be appended to either end of the sequences for region V3 and which was later removed. For region V4, only the 5′ region before the pseudoknots could be inferred in this fashion. The inferred structures for the hypervariable regions were then spliced into the core backbones determined in the previous step where these regions were forced to be unpaired and formed extended loops. Finally, the structures of the three pseudoknots and nested base pairings that RNAfold is unable to infer  were added by hand according to their homologous conserved positions and structures in the outgroup sequences.
Altogether, this procedure yielded structures that each strongly resembled the eukaryotic core structure and, within the context of the topological constraints, was based on an objective criterion for comparison (i.e., minimum free energies ) for those portions of the structure that were truly computed. The computational difficulties presented by the pseudoknots (e.g., RNAfold could not accurately estimate the structures of any of the three outgroup species in and around the regions of the pseudoknots, even given their own, complete structures as constraints) were also obviated in this manner.
Site-specific rates of molecular evolution were determined using the TIGER algorithm  as implemented in perlTiger.pl v1.0, with flanking gaps and Ns being ignored and ambiguous base calls being resolved into their component nucleotides. TIGER estimates relative rates of evolution using the congruence of site patterns between nucleotide positions as a proxy. The expectation is that slowly evolving positions will retain more phylogenetic signal and so, together with constant positions, will display more congruence with the remainder of the alignment. Conversely, positions that evolve rapidly should conflict with many others in the alignment. Thus, because TIGER determines relative rates in the absence of any phylogenetic information, it avoids some of the inherent circularity involved in determining absolute rates of evolution insofar as some a priori estimate of these rates (e.g., through a model of evolution) is needed to determine the molecular branch lengths from which the final rates are derived. Relative rates of evolution were determined for the entire data set as well as for the bdelloid and monogonont subsets of it. All analyses were restricted to the rotifer sequences (i.e., excluding C. calicophorum). Rates were only analysed for those sites where information was present for at least 15% of the species in the data set, a level chosen such that relative rates could still be calculated over the full rotifer data set for positions that comprised bdelloid sequences only (e.g., for bdelloid-specific insertions), but then still required at least 30 bdelloid species to be represented.
Non-parametric statistical analyses of the data were performed using the macOS version of PAST v4.02  to test whether or not each of (1) inferred paired sites, (2) sites that were inferred to be in stems versus loops, and (3) non-hypervariable versus hypervariable sites evolved at the same relative rates. The last comparison was tested both between pooled hypervariable and non-hypervariable regions (Mann–Whitney U) as well as among all individual hypervariable regions and the pooled non-hypervariable regions (Kruskal–Wallis). The nominal alpha value for each test was 0.05 and no correction for multiple comparisons was used.
Availability of data and materials
The DNA sequence data supporting the results of this article are available in the GenBank® repository (http://www.ncbi.nlm.nih.gov) under accession numbers MT522624–MT522695 and MT542324 (see also Tables S1 and S2 in Additional file 2). perlTiger.pl is freely available at www.uol.de/systematik-evolutionsbiologie/programme. All other data generated or analysed for this study are included in this published article or its associated Additional files.
Degrees of freedom
- et al.:
General time reversible
In other words
Apple® Macintosh® operating system
Mitochondrially encoded cytochrome c oxidase I
Phylogenetic Analyses Using Parsimony
Polymerase chain reaction
Randomized A(x)ccelerated Maximum Likelihood
Ribosomal deoxyribonucleic acid
Ribosomal ribonucleic acid
Tree Independent Generation of Evolutionary Rates
Visualization Applet for RNA
Gamma distribution for modeling rate heterogeneity
Neefs JM, Van de Peer Y, De Rijk P, Chapelle S, De Wachter R. Compilation of small ribosomal subunit RNA structures. Nucleic Acids Res. 1993;21(13):3025–49.
Gillespie JJ, Johnston JS, Cannone JJ, Gutell RR. Characteristics of the nuclear (18S, 5.8S, 28S and 5S) and mitochondrial (12S and 16S) rRNA genes of Apis mellifera (Insecta: Hymenoptera): structure, organization, and retrotransposable elements. Insect Mol Biol. 2006;15(5):657–86.
Pleij CWA. Pseudoknots: a new motif in the RNA game. Trends Biochem Sci. 1990;15(4):143–7.
Ki J-S. Hypervariable regions (V1–V9) of the dinoflagellate 18S rRNA using a large dataset for marker considerations. J Appl Phycol. 2012;24(5):1035–43.
McStay B. Nucleolar organizer regions: genomic ‘dark matter’ requiring illumination. Genes Dev. 2016;30:1598–610.
Brown DD, Wensink PC, Jordan E. A comparison of the ribosomal DNA’s of Xenopus laevis and Xenopus mulleri: the evolution of tandem genes. J Mol Biol. 1972;63(1):57–73.
Zimmer EA, Martin SL, Beverley SM, Kan YW, Wilson AC. Rapid duplication and loss of genes coding for the alpha chains of hemoglobin. Proc Natl Acad Sci U S A. 1980;77(4):2158–62.
Field KG, Olsen GJ, Lane DJ, Giovannoni SJ, Ghiselin MT, Raff EC, Pace NR, Raff RA. Molecular phylogeny of the animal kingdom. Science. 1988;239(4841):748–53.
Hug LA, Baker BJ, Anantharaman K, Brown CT, Probst AJ, Castelle CJ, Butterfield CN, Hernsdorf AW, Amano Y, Ise K, et al. A new view of the tree of life. Nature Microbiol. 2016;1:16048.
Woese CR. Bacterial evolution. Microbiol Rev. 1987;51(2):221–71.
Zimmermann J, Jahn R, Gemeinholzer B. Barcoding diatoms: evaluation of the V4 subregion on the 18S rRNA gene, including new primers and protocols. Org Divers Evol. 2011;11:173–92.
Rimet F, Chaumeil P, Keck F, Kermarrec L, Vasselon V, Kahlert M, Franc A, Bouchez A. R-Syst::diatom: an open-access and curated barcode database for diatoms and freshwater monitoring. Database. 2016;2016:21.
Pawlowski J, Audic S, Adl S, Bass D, Belbahri L, Berney C, Bowser SS, Cepicka I, Decelle J, Dunthorn M, et al. CBOL protist working group: barcoding eukaryotic richness beyond the animal, plant, and fungal kingdoms. PLoS Biol. 2012;10(11):e1001419.
Stocsits RR, Letsch H, Meusemann K, von Reumoent BM, Misof B, Hertel J, Tafer H, Stadler PF. RNA in phylogenetic reconstruction. In: Wägele WJ, Bartolomaeus T, editors. Deep metazoan phylogeny: the backbone of the tree of life: new insights from analyses of molecules, morphology, and theory of data analysis. Berlin: Walter de Gruyter GmbH; 2014. p. 531–8.
Fontaneto D, Smet WH. Rotifera. In: Schmidt-Rhaesa A, editor. Gastrotricha, Cycloneuralia and Gnathifera: Gastrotricha and Gnathifera, vol. 3. Berlin: Walter de Gruyter GmbH; 2015. p. 217–300.
Sørensen MV, Giribet G. A modern approach to rotiferan phylogeny: combining morphological and molecular data. Mol Phylogenet Evol. 2006;40:585–608.
Ricci C, Melone G, Sotgia C. Old and new data on Seisonidea (Rotifera). Hydrobiologia. 1993;255–256(1):495–511.
Ahlrichs WH, Riemann O. Seisonidae. In: Schmidt-Rhaesa A, editor. Handbook of Zoology. Berlin: De Gruyter; 2019. p. 55–85.
Swofford DL. PAUP*. Phylogenetic analysis using parsimony (*and other methods) Version 4. Sunderland: Sinauer Associates; 2002.
Van de Peer Y, de Wachter R. Evolutionary relationships among the eukaryotic crown taxa taking into account site-to-site rate variation in 18S rRNA. J Mol Evol. 1997;45(6):619–30.
Van de Peer Y, De Rijk P, Wuyts J, Winkelmans T, De Wachter R. The European small subunit ribosomal RNA database. Nucleic Acids Res. 2000;28(1):175–6.
Wuyts J, Van de Peer Y, Wachter RD. Distribution of substitution rates and location of insertion sites in the tertiary structure of ribosomal RNA. Nucleic Acids Res. 2001;29(24):5017–28.
Weigand AM, Dinapoli A, Klussmann-Kolb A. 18S rRNA variability map for Gastropoda. J Molluscan Stud. 2012;78(1):151–6.
Wuyts J, Perriere G, Van de Peer Y. The European ribosomal RNA database. Nucleic Acids Res. 2004;32(1):D101–3.
Van de Peer Y, Jansen J, De Rijk P, De Wachter R. Database on the strcuture of small ribosomal subunit RNA. Nucleic Acids Res. 1997;25(1):111–6.
Wey-Fabrizius AR, Herlyn H, Rieger B, Rosenkranz D, Witek A, Welch DB, Ebersberger I, Hankeln T. Transcriptome data reveal Syndermatan relationships and suggest the evolution of endoparasitism in Acanthocephala via an epizoic stage. PLoS ONE. 2014;9(2):e88618.
Struck TH, Wey-Fabrizius AR, Golombek A, Hering L, Weigert A, Bleidorn C, Klebow S, Iakovenko N, Hausdorf B, Petersen M, et al. Platyzoan paraphyly based on phylogenomic data supports a noncoelomate ancestry of Spiralia. Mol Biol Evol. 2014;31(7):1833–49.
Sielaff M, Schmidt H, Struck TH, Rosenkranz D, Mark Welch DB, Hankeln T, Herlyn H. Phylogeny of Syndermata (syn. Rotifera): mitochondrial gene oder verifies epizoic Seisonidea as sister to endoparasitic Acanthocephala within monophyletic Hemirotifera. Mol Phylogenet Evol. 2016;96(1):79–92.
Van Geel B. A palaeoecological study of Holocene peat bog sections in Germany and The Netherlands, based on the analysis of pollen, spores and macro- and micro- scopic remains of fungi, algae, cormophytes and animals. Rev Palaeobot Palynol. 1977;25(1):1–120.
Warner BG, Chengalath R. Holocene fossil Habrotrocha angusticollis (Bdelloidea: Rotifera) in North America. J Paleolimnol. 1988;1:141–7.
Poinar GO Jr, Ricci C. Bdelloid rotifers in Dominican amber: evidence for parthenogenetic continuity. Cell Mol Life Sci. 1992;48(4):408–10.
Waggoner BM, Poinar GO Jr. Fossil habrotrochid rotifers in Dominican amber. Experientia. 1993;49:354–7.
Hedges SB, Marin J, Suleski M, Paymer M, Kumar S. Tree of life reveals clock-like speciation and diversification. Mol Biol Evol. 2015;32(4):835–45.
Giribet G. New animal phylogeny: future challenges for animal phylogeny in the age of phylogenomics. Org Divers Evol. 2016;16(2):419–26.
Telford MJ, Budd GE, Philippe H. Phylogenomic insights into animal evolution. Curr Biol. 2015;25:R876–87.
Bergsten J. A review of long-branch attraction. Cladistics. 2005;21(2):163–93.
Folmer O, Black M, Hoeh W, Lutz R, Vrijenhoek R. DNA primers for amplification of mitochondrial cytochrome c oxidase subunit I from diverse metazoan invertebrates. Mol Mar Biol Biotechnol. 1994;3(5):294–9.
Rubinoff D, Holland BS. Between two extremes: mitochondrial DNA is neither the panacea nor the nemesis of phylogenetic and taxonomic inference. Syst Biol. 2005;54(6):952–61.
Rokas A, Holland PW. Rare genomic changes as a tool for phylogenetics. Trends Ecol Evol. 2000;15(11):454–9.
Donath A, Stadler PF. Molecular morphology: higher order characters derivable from sequence information. In: Wägele WJ, Bartolomaeus T, editors. Deep metazoan phylogeny: the backbone of the tree of life: new insights from analyses of molecules, morphology, and theory of data analysis. Berlin: Walter de Gruyter GmbH; 2014. p. 549–62.
Hu F, Lin YC, Tang J. MLGO: phylogeny reconstruction and ancestral inference from gene-order data. BMC Bioinf. 2014;15:354.
Bernt M, Merkle D, Middendorf M, Schierwater B, Schlegel M, Stadler PF. Computational methods for the analysis of mitochondrial genome rearrangements. In: Wägele WJ, Bartolomaeus T, editors. Deep metazoan phylogeny: the backbone of the tree of life: new insights from analyses of molecules, morphology, and theory of data analysis. Berlin: Walter de Gruyter GmbH; 2014. p. 515–30.
Wiens JJ. Missing data, incomplete taxa, and phylogenetic accuracy. Syst Biol. 2003;52(4):528–38.
Wilke T, Ahlrichs WH, Bininda-Emonds ORP. A comprehensive and integrative re-description of Synchaeta tremula (Muller, 1786) and the newly rediscovered Synchaeta tremuloida Pourriot, 1965 (Rotifera: Synchaetidae). Zootaxa. 2017;4276(4):503–18.
Wilke T, Ahlrichs WH, Bininda-Emonds ORP. The evolution of Synchaetidae (Rotifera: Monogononta) with a focus on Synchaeta: an integrative approach combining molecular and morphological data. J Zool Syst Evol Res. 2020;58:2.
Kimpel D, Gockel J, Gerlach G, Bininda-Emonds ORP. Population structuring in the monogonont rotifer Synchaeta pectinata: high genetic divergence on a small geographical scale. Freshw Biol. 2015;60(7):1364–78.
Melone G, Ricci C. Rotatory apparatus in Bdelloids. Hydrobiologia. 1995;313(314):91–8.
Edgar RC. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 2004;32(5):1792–7.
Gouy M, Guindon S, Gascuel O. SeaView version 4: a multiplatform graphical user interface for sequence alignment and phylogenetic tree building. Mol Biol Evol. 2010;27(2):221–4.
Raden M, Ali SM, Alkhnbashi OS, Busch A, Costa F, Davis JA, Eggenhofer F, Gelhausen R, Georg J, Heyne S, et al. Freiburg RNA tools: a central online resource for RNA-focused research and teaching. Nucleic Acids Res. 2018;46(W1):W25–9.
Katoh K, Standley DM. MAFFT multiple sequence alignment software version 7: improvements in performance and usability. Mol Biol Evol. 2013;30(4):772–80.
Stocsits RR, Letsch H, Hertel J, Misof B, Stadler PF. Accurate and efficient reconstruction of deep phylogenies from structured RNAs. Nucleic Acids Res. 2009;37(18):6184–93.
Wuyts J, De Rijk P, Van de Peer Y, Pison G, Rousseeuw P, De Wachter R. Comparative analysis of more than 3000 sequences reveals the existence of two pseudoknots in area V4 of eukaryotic small subunit ribosomal RNA. Nucleic Acids Res. 2000;28(23):4698–708.
Akutsu T. Dynamic programming algorithms for RNA secondary structure prediction with pseudoknots. Discrete Appl Math. 2000;104:45–62.
Lyngsø RB, Pedersen CNS. RNA pseudoknot prediction in energy-based models. J Comput Biol. 2000;7(3–4):409–27.
Lorenz R, Bernhart SH, HönerzuSiederdissen C, Tafer H, Flamm C, Stadler PF, Hofacker IL. ViennaRNA package 2.0. Algorithms Mol Biol. 2011;6:26.
Cummins CA, McInerney JO. A method for inferring the rate of evolution of homologous characters that can potentially improve phylogenetic inference, resolve deep divergence and correct systematic biases. Syst Biol. 2011;60(6):833–44.
Hammer Ø, Harper DAT, Ryan PD. PAST: Paleontological statistics software package for education and data analysis. Palaeontol Electron. 2001;4(1):9.
Darty K, Denise A, Ponty Y. VARNA: interactive drawing and editing of the RNA secondary structure. Bioinformatics. 2009;25(15):1974–5.
Stamatakis A. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics. 2014;30(9):1312–3.
Stamatakis A, Hoover P, Rougemont J. A rapid bootstrap algorithm for the RAxML web-servers. Syst Biol. 2008;75(5):758–71.
Stamatakis A. Phylogenetic models of rate heterogeneity: a high performance computing perspective. In: IPDPS2006: 2006; Rhodos, Greece.
Felsenstein J. Confidence limits on phylogenies: an approach using the bootstrap. Evolution. 1985;39:783–91.
I thank Wilko Ahlrichs for identifying all the rotifer species that were sequenced within the group as well as Fabian Deister, Julia Hawkins, Isabelle von Gösselin, and Tanja Wilke for performing the sequencing. Wilko Ahlrichs and Michael Raupach helped to focus and correct my thoughts regarding rotifers and 18S rRNA, respectively, and Davorka Radovčić motivated me to write perlTiger.pl, which is freely available at www.uol.de/systematik-evolutionsbiologie/programme. Peter Foster, Diego Fontaneto, and an anonymous reviewer provided additional comments in review to help improve the manuscript.
Open Access funding enabled and organized by Projekt DEAL. This work was supported in part by the Deutsche Forschungsgemeinschaft (BI 825/7–1), which played no role in the design of the study; in the collection, analysis, or interpretation of data; or in writing the manuscript.
Ethics approval and consent to participate
Consent for publication
The author declares that he has no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Fasta-formatted sequence file of 18S rDNA for the three exemplar rotifer species together with structural information in Vienna format.
List of rotifer (Table S1) and acanthocephalan (Table S2) 18S rDNA sequences used in the study.
Result files from the maximum-likelihood analysis of the alignment in Additional file 2.
Fasta-formatted sequence alignment of 18S rDNA for all species examined in this study.
About this article
Cite this article
Bininda-Emonds, O.R.P. 18S rRNA variability maps reveal three highly divergent, conserved motifs within Rotifera. BMC Ecol Evo 21, 118 (2021). https://doi.org/10.1186/s12862-021-01845-2
- Rate of evolution