Comparative transcriptome analysis of obligately asexual and cyclically sexual rotifers reveals genes with putative functions in sexual reproduction, dormancy, and asexual egg production

RNA extraction and library preparation

OP and CP clones of Brachionus calyciflorus homozygous for the op and wild-type alleles, respectively, were generated from the same heterozygous mother by selfing [18]. Six B. calyciflorus females were added to 600 mL cultures of 5 × 10⁵ cells/mL Chlamydomonas reinhardtii in COMBO medium [25]. The cultures grew at 25 C and with 24 hour light exposure for ~6 days. Density of females was determined every day by removal of ~10 mL of culture and counting individual females under a microscope at low magnification. When cultures reached a density of ~25-40 females/mL, they were harvested for RNA. To confirm the presence of sexually-reproducing females in the CP cultures, a sample of females was used to determine rate of mixis induction as described previously [26], and the absence of males and resting eggs was confirmed for the OP cultures. The rotifer cultures were filtered through 20 μm Millipore nylon filters to remove algal cells, the rotifers were washed in COMBO medium, and rinsed into 15 mL conical tubes. Total RNA was isolated using Trizol reagent (Ambion, Life Technologies, Grand Island, NY) according to manufacturer’s instructions. RNA quality was determined by Experion Automated Electrophoresis (Bio-Rad Laboratories, Hercules, CA) according to manufacturer’s instructions before construction of mRNA-seq libraries using the mRNA-seq Sample Prep kit (Illumina, Inc., San Diego, CA). Seventy-six base pair single end reads were sequenced from each library on an Illumina Genome Analyzer IIx platform at the Iowa State University DNA Facility (http://www.dna.iastate.edu/nextgen.html).

Illumina sequence analysis

Analysis of Illumina sequence reads was performed in part using the Galaxy server (http://galaxyproject.org/) [27–29]. Adapter sequences were trimmed from the ends of the reads, and a blastn search revealed reads that matched C. reinhardtii transcripts [30]. Reads with an alignment length ≥ 20 base pairs matching an algal sequence were removed from the libraries. Transcriptome assembly was performed using the Tuxedo pipeline, where reads were mapped to a partial assembly of the B. calyciflorus genome [31]via Tophat [32], assembled into transcripts using Cufflinks (assembly for individual samples with 0.1 minimum isoform fraction and 0.1 pre-mRNA fraction) and Cuffmerge (combined assembly from all samples). Differential expression was determined between OP and CP samples with Cuffdiff using a minimum alignment count of 500 [33, 34]. Differential expression between OP and CP samples was also determined with edgeR using a table of counts constructed with the fragments per kilobase locus length per million reads mapped (FPKM) values determined by Cufflinks [34]. Gene ontology mapping and enrichment analyses were performed with BLAST2GO (http://www.blast2go.com/b2ghome). The transcriptome assembly projects have been deposited at GenBank under the accessions GACQ00000000 and GACL00000000. The versions described in this paper are the first versions, GACQ01000000 and GACL01000000. Assembled transcripts shorter than 200 base pairs are given in Additional files 1, 2, 3 and 4.

Library validation

OP and CP RNA samples used for Illumina library construction were treated with DNase I prior to first-strand cDNA synthesis using Superscript II reverse transcriptase (Invitrogen, Life Technologies, Grand Island, NY). Quantitative reverse transcriptase polymerase chain reaction (qRT-PCR) was performed with Platinum SYBR Green quantitative PCR (qPCR) Supermix (Invitrogen, Life Technologies, Grand Island, NY) and a 1:25 dilution of the cDNA. Primers designed for amplification of thirteen genes from cDNA in the qPCR reaction (200 nM final concentration) are listed in Additional file 5. Relative quantification of gene expression was determined for each gene, with actin transcript amplification used for normalization.

Phylogenetic analysis

For analysis of histone H2A, amino acid sequences of metazoan and fungal homologs were obtained from the National Center for Biotechnology Information (NCBI) database by keyword search and aligned with B. calyciflorus sequences by ClustalX version 2.1 [35]. Alignments were manually curated in MacClade (http://macclade.org/index.html) and phylogenetic analysis was performed using the Bioportal parallel computational resource [36]. Maximum likelihood analysis was performed using PhyML [37] with Whelan and Goldman (WAG) substitution model, best tree topology search, and one thousand bootstrap replicates.

Bdelloid sequence analysis

Transcript sequences for Adineta ricciae were obtained from GenBank (Accession numbers: HE687510 to HE716431). A. ricciae transcripts were used as queries against the assembled B. calyciflorus transcriptome using tblastx [38] with a Bit Score cutoff of 50. Gene ontology term enrichment was performed using BLAST2GO [39].

Mixis induction

The trigger for the transition from asexual to sexual reproduction in many species of monogononts is a secreted protein signal that accumulates in their environment [14]. Although little is known about the signal or the pathway(s) it activates to induce the transition, several lines of evidence indicate that it is mediated by steroid hormone signaling. First, isolation of the protein signal and sequencing of its N-terminus revealed that its sequence in the species B. plicatilis is similar to the mammalian steroidogenesis-inducing protein [14]. Second, exposure to hormones, including estradiol-17 and progesterone, resulted in increased production of sexual females in populations of B. plicatilis[21, 45]. Third, RNAi knockdown of the progesterone receptor in B. manjavacas reduced production of sexual females [20]. Intriguingly, there appears to be species-specificity in the response to steroid hormones, as sexual reproduction is not induced by progesterone treatment alone in B. calyciflorus cultures [21, 22].

Steroid hormone signaling is mediated by members of the nuclear receptor superfamily, a large family of proteins that are well-conserved across Metazoa [46]. Nuclear receptors have critical functions in developmental processes, including the induction of meiosis during oogenesis in mammals [47, 48]. Members of the nuclear receptor superfamily are transcription factors whose activity is modulated through direct interactions with steroid ligands. Nuclear receptors frequently act as homodimers or heterodimers when modulating transcription of target sequences via histone acetylation or deacetylation, including auto-regulating and cross-regulating the transcription of other nuclear receptors [48, 49].

We identified 73 genes in our transcriptome analysis with significant sequence identity to the DNA-binding domains of the nuclear receptor superfamily, of which five have significantly higher expression levels in CP populations (Additional file 12). Although we were unable to obtain phylogenetic resolution to confirm the orthology of the receptors (data not shown), these data are consistent with the role of steroid signaling in mixis induction.

We also identified differentially expressed genes with functions and GO term designations associated with steroid metabolic and signaling pathways (Additional file 9). CBP/p300 is a transcription factor with histone acetyltransferase activity that acts as a co-regulator of transcription with nuclear receptor superfamily members [50]. Protein Kinase A modulates nuclear receptor activity via phosphorylation and has roles during gametogenesis [51]. Bcl-X mediates apoptosis and is regulated through steroid hormone signaling pathways [52]. These genes provide additional candidates for roles in the transition from asexual to sexual reproduction in B. calyciflorus.

Meiosis, recombination, and gametogenesis

Genes with known roles in meiosis based on evidence in well-studied model systems (e.g. yeast, mouse, and nematodes) were previously annotated in monogononts [31]. The expression levels of these genes are shown in Additional file 11. The distribution of expression values for these genes as a group is significantly different in CP than in OP strains, with a greater median expression value in the CP strain (Table 2, Additional file 10). This difference in expression value distribution is also found in more specific categories of gene function within the meiotic gene inventory (Table 2, Additional files 10 and 11).

Histone H2A variants

Histone octamers comprised of the histones H2A, H2B, H3, and H4 are the primary protein component of chromatin. Histones are globular proteins that have a highly conserved core sequence and a C-terminal tail sequence that is post-translationally modified to alter chromatin structure and transcriptional activity. Canonical histones may be replaced by variants in a spatial or temporal manner in response to internal or external cues [56]. In particular, variants of histone H2A are associated with DNA repair and gametogenesis. Histone H2A.X localizes to DNA double-strand breaks, where it recruits components of repair machinery. This occurs when breaks are induced through exposure to endogenous or exogenous agents, as well as in meiosis, when double-strand breaks are created by the topoisomerase Spo11 (Fernandez, 2004). Another histone variant specific to mammals, H2A.bbd, localizes to transcriptionally active sites within mouse testes, where genes with meiotic expression patterns are located (Soboleva, 2011).

Bdelloid rotifers, which apparently reproduce solely through ameiotic parthenogenesis, possess three histone H2A variants, some with multiple identified alleles. These H2A sequences have long C-terminal tail sequences that do not correspond to the tails of canonical or known variants of H2A, including the well-conserved SQE/D sequence in the tail of H2A.X [57]. Although functionally uncharacterized, because of the importance of H2A variants to double-strand break repair in other systems, it is hypothesized that these novel variants may have evolved to aid the bdelloids in DNA repair following bouts of desiccation, which the rotifers can survive at any life stage. This same report identified three apparently canonical histone H2A sequences in the monogonont rotifer B. plicatilis[57].

We searched the transcriptome and partial genome assembly for B. calyciflorus to identify all histone sequences (Table 3, Additional file 11), and identified 13 histone H2A variants (Figure 3). Previous histone H2A sequence analysis was unable to resolve histone H2A relationships at deeper phylogenetic nodes, although the variant H2A.Z is monophyletic [57]; we found similar results through our analysis. Although H2A.Z was not identified in the previous analysis of monogonont H2A sequences [57], two of the B. calyciflorus sequences fall within this clade with strong bootstrap support, suggesting that loss of H2A.Z in rotifers may be limited to the bdelloid lineage. Six B. calyciflorus H2A sequences cluster with bdelloids, although there is only moderate bootstrap support (<70%) for this relationship. As is typically observed in H2A genes, no introns were found in 12 of the genes [58], however H2A11 contains three introns, and encodes a much longer C-terminal tail than the other H2A genes in B. calyciflorus. A long C-terminal tail has previously been described for macroH2A, a vertebrate-specific H2A variant [59].

Our transcriptome analysis indicated that three of the histone H2A sequences are differentially expressed between OP and CP populations, with higher expression in the CP strains (Table 3). Differential expression of these sequences suggests a role for these H2A genes in processes related to sexual reproduction in monogononts, perhaps similar to the functions performed by H2A.X or H2A.bbd in other systems.

In addition to these novel H2A sequences, we identified two histone interacting genes with significantly higher expression in the CP populations (Additional file 9). Tonsoku interacts with both histone H2A and H2B, and has a role in resolution of stalled replication forks and double strand breaks [60]. Pax-Interacting Protein 1 recruits the H3K4 methyltransferase complex to trimethylate specific sites and activate transcription [61], and is required for both embryonic and somatic tissue differentiation [62]. H3K4 trimethylation is also required during meiosis for pairing of homologs and progression through pachytene [63].

Dormancy

Further, enrichment for other terms likely relating to resting egg components/structure was identified, such as chitin binding (GO:0008061) and chitin metabolic process (GO:0006030). Asexually produced eggs are encased in one chitin-containing shell, while resting eggs have two shells, with the inner shell of the resting eggs—also composed of chitin—thicker than the asexual egg shell [43], providing a likely explanation for the enrichment of these terms in the differentially expressed genes. Interestingly, these GO terms were not identified as enriched in the previous resting egg transcriptome analyses [24], which focused on transcripts exclusively in the eggs. This discrepancy may indicate the maternal production of these transcripts during early stages of egg development, prior to egg extrusion.

In addition to significantly different distribution of expression values for dormancy genes as a group in CP samples (Table 2, Additional file 10), we identified several individual genes with significantly higher expression in the CP strains that were previously associated with dormancy in the monogonont B. plicatilis (Additional file 11). Aquaporin 3 is involved in water and small molecule transport [94]. Genes involved in lipid metabolism, including cathepsin L, fatty acid binding genes, and lipases may be involved in yolk processing that differs between resting and amictic eggs in B. calyciflorus[95, 96]. Heat shock proteins are molecular chaperones that prevent the aggregation of misfolded proteins and have been implicated in dormancy in rotifers and other systems [97]. Trehalose is a likely energy source for dormant organisms, suggesting a role for trehalose metabolic enzymes in monogononts [98, 99]. Late embryogenesis-abundant genes have been associated with embryonic development and dormancy in several organisms, including rotifers [100, 101].

Asexual reproduction

The number of genes with significantly higher expression in the OP strain was much lower than that for the CP strains (Additional file 9). This is expected given that asexually reproducing females cannot be morphologically distinguished from sexually reproducing females in the CP strains. These females were therefore also included in CP samples, which likely decreased the fold change in expression of any gene specific to asexual females and/or their eggs. However, several interesting genes were found to have significantly higher expression in the OP strains.

Gene expression in bdelloid rotifers

Using a published gene expression dataset containing 28,922 assembled transcripts from the bdelloid rotifer A. ricciae[44], we identified transcripts with significant sequence similarity to the B. calyciflorus sequences. Using tblastx with a bit score cutoff of 50, a total of 12,512 (~31%) of the B. calyciflorus sequences had at least one significant A. ricciae transcript hit (Additional file 14).

We further examined how the proportion of B. calyciflorus sequences with significant sequence similarity to A. ricciae transcripts changed with the difference in expression observed between CP and OP strains of B. calyciflorus. If genes with increased expression in the CP strain represent those that function in a process specific to sexual reproduction, it is likely that no identifiable homolog is present in bdelloids, as selection would no longer maintain the gene for that function. In contrast, genes that show no difference in expression likely have constitutive functions and would be more likely to be conserved in bdelloids. Finally, the presence of genes with higher expression in the OP strain in the A. ricciae transcripts would be more likely given that they may function in ameiotic egg production.

The distribution of A. ricciae transcripts across the observed CP/OP fold change in B. calyciflorus is shown in Figure 4 and Additional file 8. As expected, B. calyciflorus genes with little to no change in expression between CP and OP were the most likely to have at least one significant hit in the A. ricciae transcripts, with a decrease in the proportion of genes with corresponding bdelloid sequences when expression is higher in the CP strain. However, for genes with the greatest increase in expression in the CP strains, the proportion of sequences with at least one significant A. ricciae transcript hit increases. To examine which genes account for this increase in A. ricciae transcript hits, a GO term enrichment analysis was performed for the B. calyciflorus genes that had both a significant increase in expression in the CP strains and significant sequence similarity to at least one A. ricciae transcript (Table 5). The enriched GO terms represent genes with roles in responses to antioxidant stress. Although in monogononts dormancy and desiccation tolerance only occur in resting eggs, bdelloids are capable of entering dormancy at any life stage [15, 105]. Notably, the A. ricciae transcripts were assembled from a pool of samples obtained from bdelloids under normal conditions as well as various stages of desiccation recovery [44]; thus, expression of genes with functions in dormancy (Additional file 11) may be attributed to the inclusion of transcripts from the recovering samples. However, constitutive expression of genes associated with dormancy cannot be excluded: an extremely efficient antioxidant response has been observed in bdelloids after treatment with ionizing radiation, suggesting such proteins are present to allow for a more rapid response to environmental changes [106].

A decrease in the proportion of B. calyciflorus genes with significant sequence similarity to A. ricciae transcripts was also observed for genes with increased expression in the OP strain (Figure 4, Additional file 6). While this could suggest a larger trend for more highly regulated genes in monogononts having faster rates of evolution or species-specificity, it may also suggest that the genes associated with asexual reproduction have roles in the suppression of sexual processes; these genes would be unnecessary in the bdelloids, as they are not known to undergo the reproductive switch that occurs in monogononts.

Finally, we examined the identified genes in B. calyciflorus with known roles in mixis induction, meiosis, and gametogenesis that had significant sequence similarity to A. ricciae transcripts. Of the seventy-three nuclear receptor superfamily genes identified in B. calyciflorus, twenty-eight had at least one significant A. ricciae transcript hit (Additional file 12). Three of the receptors with A. ricciae hits are differentially expressed between CP and OP strains of B. calyciflorus. However, as mentioned above, the highly conserved DNA binding domain of nuclear receptors makes it difficult to differentiate members within the superfamily, so additional studies will be required to determine the roles of these genes in both monogononts and bdelloids.

For genes with potential roles in meiosis, seven of the histone H2A genes identified in B. calyciflorus have significant hits in the A. ricciae dataset (Table 3). No significant hits were found for H2A.Z, consistent with a previous study suggesting this variant is absent in bdelloids [57]. The three H2A genes that have significant differential expression in B. calyciflorus have at least one A. ricciae transcript with significant sequence similarity. Of the 122 monogonont genes with known roles in meiosis in model systems, 73 have A. ricciae transcripts, including 31 genes with roles in chromosome structure and interaction, 16 genes involved in recombination, and 26 genes involved in cell cycle progression (Additional file 11). CDC20A, which is differentially expressed in B. calyciflorus CP and OP strains, does not have a significant A. ricciae transcript hit, although A. ricciae transcripts were identified for two other copies, CDC20C and CDC20D. For “meiosis-specific” genes, no A. ricciae transcripts were identified with the exception of the MutS homologs Msh4 and Msh5. These genes act as a heterodimer in the resolution of the double-Holliday junction during meiotic recombination [107], and have also previously been identified in bdelloids via degenerate PCR (Schurko et al., unpublished).

Of the 39 genes with putative functions in male gametogenesis, twenty-two have significant hits in the A. ricciae transcripts (Additional file 13). The presence and expression of these genes in bdelloids combined with the lack of significant differential expression between CP and OP strains of B. calyciflorus suggests that roles for these genes are not exclusive to spermatogenesis in rotifers.