Genome-wide analysis of 3′-untranslated regions supports the existence of post-transcriptional regulons controlling gene expression in trypanosomes

Conserved structural RNA elements in 3′-UTRs of mRNAs encoding metabolically interrelated proteins of T. cruzi

Using a bioinformatic approach, we focused our work on the detection of potential structural cis-elements located in noncoding regions of functionally related transcripts. Towards this goal, we carried out a global analysis using T. cruzi genomic data (El-Sayed et al., 2005). Since there is no information available for T. cruzi RNA-seq reads, noncoding sequences have been inferred from average lengths of 5′- and 3′-UTRs of T. cruzi published transcripts (Brandao & Jiang, 2009; Campos et al., 2008) and extracted from TriTrypDB (http://tritrypdb.org/tritrypdb/) (Methods). These sequences were classified into functional categories in accordance to the KEGG pathway database (Kanehisa & Goto, 2000) (http://www.genome.jp/kegg/pathway.html). We next generated lists of putative 3′ noncoding regions for each KEGG class containing genomic sequences resembling 3′-UTRs. Allelic copies identified in the hybrid TcVI CL Brener genome having similar 3′-UTRs, 80% identity or higher, were filtered to reduce redundancy (see File S1 and Methods for details). As a result, we categorized the T. cruzi proteins within 80 groups including 1814 genes, but only those classes having at least 10 sequences were used in this paper. Thus, we limited our search to 53 categories termed tcr00010 to tcr04650 (see Table 1 for descriptions) which enclose 1617 total genes.

Linear motifs are reportedly difficult to define, especially in repeat-rich and atypical TriTryp genomes which contain pyrimidine-rich elements in the intergenic regions (El-Sayed et al., 2005; Hendriks & Matthews, 2007). Thus, we used the CMfinder software (Yao, Weinberg & Ruzzo, 2006) (http://bio.cs.washington.edu/yzizhen/CMfinder/) for structural RNA motif prediction in the putative 3′-UTR sequences of each group. Covariance models are RNA motif models that represent both the sequence and structure binding preferences of RBPs. We chose the best top ranked motif provided by the program. Therefore, 53 new RNA structural elements were identified and termed according to the number of the KEGG pathway from which the motifs were obtained: e.g., m00010 is the motif derived from the tcr00010 dataset (Glycolysis/Gluconeogenesis), m00020 from tcr00020 (Citrate cycle), etc. Figure 1 illustrates the motif discovery pipeline used (Fig. 1A) and a pie chart distribution of the metabolic groups having at least 10 genes used as the input data (Fig. 1B).

Figure 2 shows the RNA structures for the predicted motifs. Structured elements had a length ranging from 28 nts (tcr00240, Pyrimidine metabolism) to 87 nts (tcr03010, Ribosome). Nearly all the consensus motifs fold as a predicted stem-loop structure, with an average hairpin length of 15 bp and a loop ranging from 3 to 18 nts, giving rise to loops of a median length of 4 nts. Based on the logo representation, some motifs were classified according to their nucleotide composition. File S2 shows the consensus sequence, secondary structure in bracket notation and sequence logo of all the candidate RNA elements.

Evaluating the significance of the motif enrichment by randomization tests

Next, we further analyzed the specific-enrichment of the RNA elements in the KEGG groups. Consequently the motif representation was calculated as the percentage of element-containing sequences over the total number of sequences in each category (detailed under Methods). Overall, 79% of the groups have specific RNA elements. Accordingly, 42 out of 53 KEGG categories encompassed conserved structural motifs statistically enriched in their 3′-UTRs in comparison with control groups using random 3′-UTR datasets (Z-test, FDR 10%) (Fig. 3 and Table S1). For example, the RNA motif m00030 that was originally discovered in the 3′-UTRs of genes from the Pentose phosphate pathway (tcr00030) was detected in 69.7% of the sequences (23 of 33 mRNAs) but only in 26% ± 6% against 50 random searches.

These results reveal that ≈ 80% of the metabolic pathways analyzed contain specific signature elements in their noncoding regions different to what could be expected by chance in a group of random sequences. To see if these RNA motifs are specifically located at 3′-UTRs, we repeated our searches using the elements identified in the 3′-UTRs as queries against a 5′-end dataset. This dataset is composed of 350-nt fragments of the 5′-end of T. cruzi transcripts, grouped according to KEGG pathways (Methods). We next calculate the relative frequencies of elements (number of hits found in each group divided by its sequence length). Virtually all the motifs were noticeably over-represented in the 3′-UTRs compared with the 5′-end subset (Fig. S1), indicating a preferred 3′-UTR localization. See File S3 for a complete list of element-containing genes.

Constraints in motif orientation and position

RNA-binding elements need to be oriented in a particular direction relative to the coding sequence in order to carry out their regulatory function (Elemento, Slonim & Tavazoie, 2007). Thus, motif sequences are biologically significant when located on one strand, but not on the other. To evaluate sequence orientation preference, we compared the motif occurrences on the transcribed strand (genuine transcript) against the information acquired by its occurrences on the non transcribed strand (anti-sense transcript). The coverage test showed that most of the motifs previously described (38 out of 42) have an orientation bias, with a particular sense orientation (χ² test, FDR 5%) (Table S2 and Fig. S2). This property is consistent with RNA regulatory motifs located in a precise orientation with respect to the coding sequence. A definitive list containing the motif’s representation of the 38 significant RNA elements and their P-values are shown in Table 2.

Because there is no data available for the full T. cruzi transcriptome, we searched a trypanosomal transcript database for the presence of Expressed Sequence Tag records that could match the predicted 3′-UTR used in this work. About 10% of the total motif-containing KEGG sequences (107 out of 1090 non-redundant genes) gave a positive blast hit in the database considering an alignment of more than 150 nts with >95% of identity. Remarkably, 78.5% (84 out of 107) of these sequences contain the candidate element fully aligned against the EST transcript sequence. Table S3 shows a chart with the BLAST output of the KEGG genes indicating the motif position in the predicted 3′-UTR. Next, we sought to determine the percentage of sequence-containing motifs within the EST for each of the 53 KEGG metabolic group. As shown in Table 1 and Fig. S3, most of the KEGG groups analyzed here (43 out of 53) contain more than 50% of the sequences having the entire RNA candidate element aligned within the ESTdb hit; thus reinforcing the idea that the motifs are located within UTRs.

In this regard, we found that the majority of the motifs that were identified in 3′-UTRs have a position bias. If we split the 350-nt sequence into two segments, almost 60% of the elements were preferentially located proximal to the CDS. The proportion of RNA motifs located in the vicinity of the ORF increase to 75% when we restricted the study to those elements having the most significant e-values (E < 10⁻⁹) (Table S4). This observation has also been reported for some regulatory elements (Elemento, Slonim & Tavazoie, 2007; Mao, Najafabadi & Salavati, 2009). In Fig. 4A we show the distribution of RNA element localizations within the 3′-UTRs. However, this distribution is dependent on how we bin data. Contrarily, the probability density plot in Fig. 4B, which produces a smoother representation of the histogram, is independent of bin size. Both Figs. 4A and 4B show a declining trend for the localization of the identified RNA motifs as we move away from the stop codon. We are aware that this declining trend (observed at both ends in Fig. 4B) may reflect an artifactual boundary effect. However, because this effect is noticeably more pronounced in the distal part of the molecule (300–350 nts) we interpret this to suggest that the optimal localization of these motifs could be proximal to the coding region.

Altogether, these results suggest that most of the sequences selected for this study could be considered as 3′ noncoding regions, with the RNA motifs most probably making part of the 3′-UTR and located at a particular distance from the stop codon, increasing the probability that the reported motifs have biological functions.

Motif representation in other metabolic pathways indicates that most of the RNA elements belong to a specific group

We next performed an all-against-all comparison of the presence of each motif in the 3′-UTRs of all KEGG categories, to evaluate the distribution of the previously identified elements in different metabolic groups. In this strategy, each RNA element is used to search into all individual KEGG T. cruzi datasets (e.g. the m00010 element is separately searched against dataset tcr00010 to tcr04650) (see Table S5). However this examination is not valid if KEGG groups share a high number of genes. Under this scenario, a significant motif density in two different groups can simply be the result of the presence of a high number of shared genes, and not of the conservation of RNA structure in the absence of sequence similarity. To avoid this situation, we first clustered the initial KEGG dataset into 23 disjoint (non-overlapped) sets (see Fig. S4 and Methods). From this set, we selected 19 non-overlapping groups that have significantly enriched RNA elements. Figure 5A shows a box plot chart containing the 3′-UTR motifs in the x-axis and the percentage of motif-containing mRNAs from each KEGG dataset in the y-axis. It is noteworthy that 15 out of 19 elements (78.9% of total cases) were statistically more abundant in mRNAs of the KEGG dataset from which the motif was initially identified. For example, the m00330 has a motif representation of ≈ 65% in the KEGG tcr00330 dataset, but less than 40% in all the remaining groups. Arrows above each motif represent the 15 significant motifs that are specifically enriched in their classes (group-specific motifs). Only in four cases the motif ranked at the second position: m00280, m0970, m03050 and m03010 (see below).

Additionally, data were also visualized by a heatmap analysis using R package (Fig. 5B). The heatmap plot shows that [1] as anticipated, most of 3′ motifs were specifically enriched in the group from which they were derived (a white/yellow color indicates high correlation, red color indicates low correlation); [2] some elements displayed lower abundances: m00480 and m04120; and [3] other elements are widely enriched and have a high abundance in 3′-UTR mRNAs of metabolically unconnected groups: m00030, m00280, m04130 and m03010, being this last motif the one having the most widespread distribution. This evidence indicates that structured 3′-UTR elements can constitute a signal network, being mostly present in groups of mRNAs that belong to a given metabolic pathway (group-specific elements) or being part of a wide-ranging distribution. The dendogram depicted in Fig. 5B separates motifs into three main clusters. Cluster I contains 11 motifs involved in specific cellular processes. In addition, four elements from Cluster II are specifically enriched in the group from which they were discovered. However, it is noteworthy that these RNA elements (m00030, m00620, mm03020 and m04130) also present an important enrichment in the 3′-UTR of other groups. In this context, we speculate that certain cis-acting signals may be shared among different but related metabolic categories, i.e. m00030 (Penthose phosphate pathway) and m00620 (Pyruvate metabolism) thus organizing a combinatorial network of RNA-RBP interactions. Finally, Cluster III contains the global element m03010. The Top 2 RNA motifs extensively distributed, m00280 and m03010, also presented the highest abundances in our validation strategy using randomized datasets. This produces a high representation of these motifs among unconnected metabolic groups, which hampers the analysis of their biological relevance. Altogether, the estimation of motif density among groups that did not share any transcripts (19 non-overlapped datasets) showed that at least 15 group-specific elements were more abundant in their corresponding category compared with all remaining groups. Nevertheless, a few elements such as m00030 and m04130 were also enriched in other datasets. A similar analysis was performed using 27 metabolic groups sharing less than 25% of the genes and comparable results were obtained [85% of the RNA elements were group-specific (data not shown)].

Validating predicted motifs using independent genomic and gene ontology data sets

In previous sections we found that a number of RNA motifs are specifically enriched in groups of functionally interrelated coding genes. To obtain an independent validation of the representation of these motifs in functional categories, we performed a reverse validation in which we searched the complete T. cruzi genome with each of these 53 motifs, using the cmsearch algorithm. We next analyzed the positive hits for each motif to see if there was a significant enrichment in annotated functional categories in each group. To do this we fetched the results of these searches into the DAVID server (NIAID, National Institutes of Health) (Huang et al., 2007) to evaluate functional enrichment in each of the 53 lists of genome motif-containing targets. The results show that 73.6% of the putative elements (39 out of 53) were successfully used to predict related gene targets with the same biological function of the category from which the motif was originated (results available in Table S6).

To cross-validate the motifs, we used a dataset grouped according to the Gene Ontology project (http://www.geneontology.org/). We analyzed 93 groups of 3′ downstream sequences having at least 10 sequences each (see File S4). We found that 35 out of 38 motifs (92%) were over-represented in at least one GO category compared with random searches. As expected, 11 out of 35 elements (31%) were enriched in groups of closely related function: m00270, m00310, m00480, m00562, m03010, m03050, m03420, m03430, m03440, m04120 and m04130 (Table S7). These results show that ≈ 74% of the motifs analyzed here were successfully validated against a full search in the T. cruzi genome and a third part of the elements could also be confirmed using a different dataset grouped by ontology characters.

Co-occurrence of specific RBP recognition sites and regulon’s signals

When we inspected the list of mRNAs co-immunoprecipitated in TcDHH1-containing granules (Holetz et al., 2010) (cellular structures that appear to be more related to a RNA degradative rather than stabilizing process), we could not detect over-representation of any elements, with the exception of m00250 (Ala, Asp and Glu metabolism) (data not shown). This observation suggests that aggregation of transcripts in TcDDH1 granules may not be guided by specific cis-acting signals but for a general recruitment’s mechanism, still unknown.

We next searched for co-occurrence of experimental RNA-binding elements (previously described in our laboratory) and structured motifs identified in this study. To this end, we examined the targets of UBP1 and RBP3 -containing the recognition motifs UBP1m and RBP3m (Noe, De Gaudenzi & Frasch, 2008)- and counted the number of targets containing (or not) any of the candidate elements. The data obtained revealed that the RBP3 binding motif was enriched in the Ribosome KEGG group (tcr03010 dataset) (P < 0.05, χ² test) (see Table 3). When UBP1 mRNA target hits were analyzed, we found that the UBP1 binding motif was enriched in six KEGG groups (m00230, m03420, m00020, m00010, m00785, m03450) compared with the entire T. cruzi genome representation (χ² test, Bonferroni correction). Therefore, for seven categories a significant number of genes have a co-occurrence of both sequences (one specific RBP binding site and one predicted KEGG element) in their 3′-UTRs, suggesting that these particular groups could be coordinately regulated by specific trans-acting factors and thus defining a post-transcriptional regulon.

RNA motifs are differentially expressed during parasite development and stress response

We used T. cruzi microarray data provided by the Tarleton laboratory (Minning et al., 2009) to investigate the motif representation among clusters of developmentally co-expressed genes. Firstly, we used the coXpress v1.3 program (http://coxpress.sf.net) to obtain 74 clusters of co-expressed genes. For each sequence in a cluster, we obtained 350 nts downstream of the stop codon of the CDS as annotated in TriTrypDB to obtain sequences resembling 3′-UTR (see Methods). Secondly, predicted 3′-UTRs of transcripts included in these groups were utilized to analyze the motif density for each putative RNA element. To analyze the statistical significance of the enrichment of these elements, we compared the experimental data against random distributions (P < 0.001, Z-test). Interestingly, 18 out of the 53 RNA motifs were statistically enriched in 20 groups of developmentally regulated (co-expressed) genes (see Fig. 6 and File S5). From these 18 elements, 13 (72%) belong to the 38 statistically significant candidates listed in Table 2: m00310, m00380, m00410, m00480, m00500, m00562, m00670, m00900, m00970, m03430, m03440, m04070 and m04120 (see Table 4 for a detailed list of developmentally regulated motifs).

Using a similar strategy we next investigated the transcriptional response of T. cruzi epimastigotes submitted to hyperosmotic stress. Using the data from Li et al. (2011), we analyzed the expression profiles of co-expressed genes and identified 33 sets of similarly regulated genes (File S6). As before, we analyzed the statistical significance of these element’s enrichments comparing the experimental data against random distributions (P < 0.001, Z-test). When epimastigote cells were subjected to hyperosmotic stress during a time-course experiment, nine RNA motifs were statistically over-represented in 10 co-regulated gene clusters (CI, CII, CIII, CVIII, CIX, CX, CXI, CXII, CXXII and CXXIII) compared to random searches (P < 0.001, Z-test). From these nine elements, six (m00310, m00510, m00562, m00670, m00785 and m03440) belong to the 38 statistically significant motifs (Table 5).

There are four RNA motifs that are over-represented in four clusters of genes differentially up-regulated by stress. These elements are: m00310 (in CII and CXI), m00562 (in CII, CIII and CVIII), m00670 (in CII) and m00785 (in CVIII). On the other hand, two structured motifs were enriched in two clusters of down-regulated genes: m00510 (in CXII) and m03440 (in CXII and CXXIII). These results suggest that a number of candidate RNA elements identified in this work may be involved in the post-transcriptional regulation of a variety of genes whose expression changes significantly during the parasite’s life-cycle or upon stress.

Databases

The T. cruzi database (T. cruzi CL Brener genomic sequence Release 5.1) utilized in this work was obtained from TriTrypDB Database Resource (www.tritrypdb.org) (Agüero et al., 2006). 5′ upstream genomic sequences and 3′ downstream genomic sequences were obtained using the TriTrypDB sequence retrieval tool. A length of 80 nts upstream to the CDS was used to obtain sequences resembling the 5′-UTR, while 350 nts downstream to CDS were used for 3′-UTR, in agreement to previously reported data from trypanosomes (Campos et al., 2008). The 5′-end of transcripts was estimated as a 350-nt sequence encoding the predicted 5′-UTR followed by the first 270 nts of the coding region (in order to cover the same total length of the 3′-UTRs).

Trypanosomatid genes were grouped by the KEGG pathway database (Kanehisa & Goto, 2000). EST databases (T. cruzi filtered) were downloaded from NCBI and BLASTn searches were performed with the following parameters “−F F −W 7 −E 1e−5 −S 1 −b 1 −v 1 −m 8”. DAVID Functional Annotation Chart tool (Huang et al., 2007) (http://david.abcc.ncifcrf.gov/) was used to categorize and compare the different gene lists against a T. cruzi background (using Fisher’s exact test and Benjamini-Hochberg correction).

RNA motif elucidation

All computational analyses were performed using free software, available in the public domain and compiled for a LINUX environment (Ubuntu 9.10 distribution). Using CMfinder stand-alone software version 0.2 (Yao, Weinberg & Ruzzo, 2006) (http://bio.cs.washington.edu/yzizhen/CMfinder/), common elements were identified in 53 different datasets containing at least 10 sequences, using the following parameters “−s1 −f 0.6 −c10”. As a control, motifs were also obtained using random sets. The 53 candidate motifs obtained were used to build stochastic context-free grammar (SCFG) model with the Infernal program (Nawrocki, Kolbe & Eddy, 2009). These models are representations of the RNA secondary structure and were used as queries in more refined searches (cmsearch algorithm). To validate the usage of models, a full test comparison between cmsearch and CMfinder’s output was performed to determine the parameters of true and false positive rates. Thus, each model was used to assess the element coverage within the 3′ sequences in the KEGG dataset (m00010 vs. tcr00010, m00020 vs. tcr00020, etc). We found that cmsearch tool gave a very good performance, being able to find most of the originally predicted CMfinder elements (mean “sensitivity” 0.91), giving very few false positives (mean “1-specificity” of 0.13) with a mean accuracy of 0.90.

Randomization test

Random datasets were constructed using a custom Perl script (version 5.8). Fifty different groups containing 3′-UTRs were obtained by randomly shuffling the original dataset of 1814 sequences from the KEGG repertoire and used in searches with CMfinder and Infernal programs. Each KEGG motif was used to search with Infernal program (Nawrocki, Kolbe & Eddy, 2009) against 50 random groups containing a similar number of sequences.

Statistical analysis

All statistical tests were corrected for multiple-testing, using the Benjamini-Hochberg False Discovery Rate (FDR) procedure. To test statistical significance of detected motifs, we first evaluated the specific enrichment of a given element. We defined specific enrichment as the mean difference of sequences having the motif in the original set, compared to sets of non-related sequences. For this, we generated a distribution of specific-motif enrichment by searching a given motif into 50 groups of non-related 3′-UTRs (randomly generated groups). The motif coverage (i.e., the number of sequences having a motif divided the total number of sequences in the group) was normally distributed, allowing us to calculate a Z-score for each putative motif [(specific-motif coverage — mean coverage in random groups)/standard deviation]. As the expected motif coverage was sensitive to the size of the group, we actually generated distributions for different group sizes (N = 5 to N = 200) and evaluated each motif significance with the corresponding distribution (according to its size). To determine if the specific-enrichments are significant, we generated a null distribution of motif-specificity, from sets of non-related sets of UTRs. A second Z-score was calculated and used to determine whether a particular specific-enrichment is expected by chance or not. Because this test was more restrictive, and to avoid losing too much sensitivity, we set an FDR threshold of 0.1 (i.e., 10% of false positives are expected). Motif candidates were further filtered based on a motif coverage test. The motif coverage was calculated for each KEGG group using Infernal algorithm [using parameter “−toponly −E 1” (Nawrocki, Kolbe & Eddy, 2009)]. Then, we tested if the motif coverage in a given set of UTRs is significantly higher compared to the set of complementary strands (i.e., non-transcribed strands where no motifs are expected). For motif detection in non-transcribed regions, the parameter “−bottomonly” was applied. Statistical significance of the coverage difference between UTRs and non-transcribed sequences (sense vs. anti-sense) was assessed by a χ² test (FDR 5%).

Sequence logos and secondary structures

The motif logo was constructed using WebLogo (http://weblogo.berkeley.edu). RNAfold tool (http://rna.tbi.univie.ac.at/) (Hofacker, 2003) was used to plot the secondary structure of the 53 predicted RNA motifs.

Gene filtering

Many T. cruzi genes are present in many copies. To avoid such sequence bias in motif detection procedures, we filtered out paralogous genes whose UTRs have >80% of sequence identity (BLAST alignment with at least 280 of the 350 nucleotide matches). In addition, truncated UTRs (<350 nts) and UTRs having a stretch of 5 or more ambiguous nucleotides (NNNNN) were discarded. A total of 106 genes were eliminated after applying these filters using custom Perl scripts.

Generation of disjoint (non-overlapping) KEGG groups

In order to unbiasedly test if some motifs are present in more than one KEGG group, we first defined non-overlapped KEGG groups (i.e., groups that have no genes in common). A measure of overlap between a pair of KEGG groups was calculated as the number of shared genes divided by the number of total distinct genes considering both groups (set intersection/set union). Figure S4 shows a dendogram of KEGG groups after a complete-linkage clustering. All branches starting at a height of 1 are disjoint (non-overlapped), while groups joining at lower heights have increasing degrees of overlap. To make a dataset of non-overlapped groups, we simply selected one KEGG group of each branch.

Clustering

Hierarchical clustering was carried out using CoXpress 1.3 software (http://coxpress.sf.net). Briefly, 5268 genes obtained from T. cruzi microarray data previously published by the Tarleton lab (Minning et al., 2009) were clustered into 74 profiles using the following parameters: h = 0.05, r = 0.95. Stress data was extracted from Table S2 of Li et al. (2011) collected by the Docampo lab. Concisely, 1468 genes affected after hyperosmotic stress were clustered into 33 profiles using the following parameters: h = 0.01, r = 0.95. RNA motifs were used to search against clustered and 50 random groups using cmsearch algorithm and statistical significance was calculated by Z-test.

Abbreviation List

ARE

AU-rich element

mRNP

messenger ribonucleoprotein

RBP

RNA-binding protein

RRM

RNA-recognition motif

UTR

untranslated region