Differentiation of ncRNAs from small mRNAs in Escherichia coli O157:H7 EDL933 (EHEC) by combined RNAseq and RIBOseq – ryhB encodes the regulatory RNA RyhB and a peptide, RyhP

Abstract

Background

While NGS allows rapid global detection of transcripts, it remains difficult to distinguish ncRNAs from short mRNAs. To detect potentially translated RNAs, we developed an improved protocol for bacterial ribosomal footprinting (RIBOseq). This allowed distinguishing ncRNA from mRNA in EHEC. A high ratio of ribosomal footprints per transcript (ribosomal coverage value, RCV) is expected to indicate a translated RNA, while a low RCV should point to a non-translated RNA.

Results

Based on their low RCV, 150 novel non-translated EHEC transcripts were identified as putative ncRNAs, representing both antisense and intergenic transcripts, 74 of which had expressed homologs in E. coli MG1655. Bioinformatics analysis predicted statistically significant target regulons for 15 of the intergenic transcripts; experimental analysis revealed 4-fold or higher differential expression of 46 novel ncRNA in different growth media. Out of 329 annotated EHEC ncRNAs, 52 showed an RCV similar to protein-coding genes, of those, 16 had RIBOseq patterns matching annotated genes in other enterobacteriaceae, and 11 seem to possess a Shine-Dalgarno sequence, suggesting that such ncRNAs may encode small proteins instead of being solely non-coding. To support that the RIBOseq signals are reflecting translation, we tested the ribosomal-footprint covered ORF of ryhB and found a phenotype for the encoded peptide in iron-limiting condition.

Conclusion

Determination of the RCV is a useful approach for a rapid first-step differentiation between bacterial ncRNAs and small mRNAs. Further, many known ncRNAs may encode proteins as well.

Background

Bacterial RNA molecules consist of non-coding RNAs (ncRNAs including rRNAs and tRNAs), and protein-coding mRNAs. ncRNAs are encoded either in cis or in trans of coding genes and their size ranges from 50–500 nt [1, 2]. Cis-encoded ncRNA templates are localized opposite to the gene to be regulated and, accordingly, have full complementarity to the mRNA. Their expression leads to a negative or positive impact on the expression of the regulated gene [3–5]. This type of gene regulation has been exploited in applied molecular biology [6]. However, only few experimentally verified cis-encoded ncRNAs exist, in contrast to trans-encoded ncRNAs. Trans-encoded ncRNAs are usually found in intergenic regions and have a limited complementarity to the regulated gene. Recent research has led to the view that trans-encoded ncRNAs are involved in the regulation of almost all bacterial metabolic pathways (see [7], and references therein).

The number of annotated ncRNAs known from different bacterial species is rapidly increasing. For instance, 329 ncRNAs are annotated for E. coli O157:H7 str. EDL933 [2]. Around 80 of them have been experimentally verified in E. coli [8]. Numerous bioinformatic studies on E. coli K12 and other bacterial species predicted the number of ncRNAs to range between 100 and 1000 (e.g. [9–11]). As E. coli O157:H7 strain EDL933 (EHEC) contains a core genome of 4.1 Mb which is well conserved among all E. coli strains [12], many similar or identical ncRNAs are assumed to exist in EHEC.

In the past, ncRNAs have been predicted by different bioinformatics methods (see [13] for a review about ncRNA detection in bacteria). A commonly used tool in ncRNA-prediction is RNAz, which has been used to predict ncRNAs in Bordetella pertussis [14], Streptomyces coelicolor [15] and others. However, any such studies require experimental verification [13] of which next-generation sequencing is of prime interest for this task.

While experimental large scale screenings for ncRNAs, especially strand-specific transcriptome sequencing using NGS, are becoming more and more important (e.g. [16–18]), it is not possible to determine whether a transcript is translated, based solely on RNAseq (see, e.g. [19]). In order to distinguish “true” ncRNAs from translated short mRNAs, we modified the ribosomal profiling approach developed by Ingolia et al. for yeast [20] and applied this technique to E. coli O157:H7 strain EDL933. Ribosomal profiling, which is also termed ribosomal footprinting or RIBOseq, detects RNAs which are covered by ribosomes and which are, therefore, assumed to be involved in the process of translation. The RNA population which is covered by ribosomes is termed “translatome” [21] and bioinformatics tools are now available to analyze these novel data [22]. Combined with strand-specific RNA-sequencing, we suggest that this approach provides additional evidence to distinguish between non-coding RNAs and RNAs covered by ribosomes.

In the past, RNAs have been found which function as ncRNA (i.e. having a function as RNA molecule not based on encoding a peptide chain) and, at the same time, as mRNA (i.e. encoding a peptide chain). Therefore, those RNAs were either termed dual-functioning RNAs (dfRNAs [23]) or coding non-coding RNAs (cncRNAs [24]). The former name is now used for RNAs with any two different functions (e.g., base-pairing and protein binding [25]), the latter describes the fact that the DNA-encoded entity functions on the level of RNA (hence, non-coding) and additionally on the level of an peptide (i.e. coding). Less than ten examples of cncRNAs are known from prokaryotes, e.g., RNAIII, SgrS, SR1, PhrS, gdpS, irvA, and others [23, 24, 26, 27].

Methods

Microbial strain

Strain E. coli O157:H7 EDL933 was obtained from the Collection l’Institute de Pasteur (Paris) under the collection number CIP 106327 (= WS4202, Weihenstephan Microbial Strain Collection) and was used in all experiments. The strain was originally isolated from raw hamburger meat, first described in 1983 [28], originally sequenced in 2001 [12] and its sequence improved recently [29]. The genome of WS4202 was re-sequenced by us to check for laboratory derived changes (GenBank accession CP012802).

RIBOseq

Ribosomal footprinting was conducted according to Ingolia et al. [20], but was adapted to sequence bacterial footprints using strand-specific libraries obtained with the TruSeq Small RNA Sample Preparation Kit (Illumina, USA). Cells were grown in ten-fold diluted lysogeny broth (LB; 10 g/L peptone, 5 g/L yeast extract, 10 g/L NaCl) with shaking at 180 rpm. At the transition from late exponential to early stationary phase the cultures were supplemented with 170 μg/mL chloramphenicol to stall the ribosomes (about 6-times above the concentration at which trans-translation occurs [30]). After two minutes, cells were harvested by centrifugation at 6000 × g for 3 min at 4 °C. Pellets were resuspended in lysis buffer (20 mM Tris-Cl at pH8, 140 mM KCl, 1.5 mM MgCl₂, 170 μg/mL chloramphenicol, 1% v/v NP40; 1.5 mL per initial liter of culture) and the suspension was dripped into liquid nitrogen and stored at −80 °C. The cells were ground with pestle and mortar in liquid nitrogen and 2 g sterile sand for about 20 min. The powder was thawed on ice and centrifuged twice, first at 3000 × g at 4 °C for 5 min and next at 20,000 × g at 4 °C for 10 min. The supernatant was saved and A_260nm determined. After dilution to an A_260nm of 200, RNase I (Ambion AM2294) was added to the sample to a final concentration of 3 U/μL and the sample was gently rotated at room temperature (RT) for 1 h. Remaining intact ribosomes with protected mRNA-fragments (footprints) were enriched by gradient centrifugation. A sucrose gradient was prepared in gradient buffer (20 mM Tris-Cl at pH 8, 140 mM KCl, 5 mM MgCl₂, 170 μg/mL chloramphenicol, 0.5 mM DTT, 0.013% SYBR Gold). Nine different sucrose concentrations were prepared in 5% (w/v) steps ranging from 10 to 50% and 1.5 mL of each concentration was loaded to a centrifuge tube. Five hundred μL of the crude ribosome sample were loaded onto each gradient tube and centrifuged at 104,000 × g at 4 °C for 3 h. The layer containing the ribosomes was visualized using UV-light and the tube was pierced at the bottom to slowly release the gradient and the band containing intact 70S ribosomes was collected. To ensure that RNA which is not protected by ribosomes is fully digested, and to get a highly enriched ribosomal fraction, the procedure of RNase-digestion and gradient centrifugation was repeated: The ribosomal fraction was diluted 1:1 with gradient buffer (without SYBR Gold and sucrose) and was loaded on a sucrose gradient without the 10% sucrose layer. After centrifugation, complete 70S ribosomes were collected by slowly releasing the gradient as described above and frozen in liquid nitrogen. To obtain the protected ribosomal footprints, 1 mL Trizol was added to 200 μL of the ribosome suspension following the manual for Trizol extraction of RNA (life technologies, USA). The final footprint-RNA pellet was dissolved in RNase free water. To ensure no carry-over of genomic DNA fragments, DNase treatment was performed using the TURBO DNA-free Kit (Applied Biosystems, USA) according to the manual. For footprint size-selection, the crude RNA-preparation was loaded to a 15% denaturing polyacrylamide gel. An oligonucleotide of 28 bp was used as a marker which is about the size of a ribosomal footprint [31, 32]. After staining with SYBR Gold, the region of about 28 nt was excised from the gel. The RNA was extracted from the gel slice as described [20]. Results of pilot experiments showed that RNase I cuts the 5′ ends of the 16S rRNA producing a fragment of about the size expected for the footprints, contributing about 50% to the size-selected RNA fragments after sequencing. For this reason, these fragments were removed with oligonucleotides complementary to the 5′-end of the 16S rRNA using the MICROBExpress bacterial mRNA enrichment kit (life technologies, USA) following the manual. Furthermore, true footprints were found to be shorter than expected (see Results). Enriched footprint-RNAs were dephosphorylated using Antarctic phosphatase (10 units per 300 ng RNA, supplemented with 10 units Superase, 37 °C for 30 min). Footprints were recovered using the miRNeasy Mini Kit (Qiagen, Germany). Subsequent phosphorylation was carried out using T4 polynucleotide kinase (20 units supplemented with 10 units Superase, 37 °C for 60 min) and cleaned using the miRNeasy Mini Kit as before. Finally, the entire sample was processed with the TruSeq Small RNA Sample Preparation Kit (Illumina) according to the manual, using 11 PCR cycles, and was sequenced on an Illumina MiSeq.

Transcriptome sequencing

The same cultures used for ribosomal footprinting were also used for transcriptome sequencing (i.e., strand specific RNAseq). Fifty μL of the diluted cell extract with an A_260nm of 200 units (see above) were added to one 1 mL of Trizol and total RNA was isolated. Since 90–95% of the total RNA consists of ribosomal RNA [33], the Ribominus Transcriptome Isolation Kit (Yeast and Bacteria, Invitrogen, USA) was applied according to the manual and the RNA was precipitated with the help of glycogen and two volumes 100% ethanol. DNase treatment was performed as described above. One μg RNA was fragmented as described [34] and the RNA-fragments were precipitated with glycogen and 2.5 volumes 100% ethanol. For sequencing on an Illumina MiSeq, the fragments were resuspended in 25 μL RNase free water and further processed like the cleaned footprint-RNAs (see above).

Northern blots

RNA was isolated in the same manner and under the same conditions as for the NGS experiments. Northern blots were performed using the DIG Northern Starter kit (Roche, Switzerland). Primers to generate DIG (digoxygenin) labeled probes are listed in Additional file 1: Table S1. For preparation of the probes, electroblotting, crosslinking, hybridization and detection, the manufacturer’s protocol was followed, except that electroblotting was performed using polyacrylamide gels and that for crosslinking EDC (1-ethyl-3-(3-dimethylaminopropyl) carbodiimide) was used [35]. After exposure to CDP-Star (included in the DIG Northern Starter kit), luminescence activity of the hybridized probes was measured using an In-Vivo Imaging System (PerkinElmer, USA).

Competitive growth assays for the overexpression phenotype of RyhP

For the production of the peptide RyhP encoded in RyhB, two versions of the corresponding ORF (named P1 and P2) were cloned onto pBAD/Myc-His C (Invitrogen). Similarly, two versions of this ORF with either the second or the third codon changed into stop codons to terminate translation were used as negative controls (named T2 and T3). For cloning, primer pairs (for primer see Additional file 1: Table S1) were hybridized forming RyhP-coding dsDNA fragments. The pBAD was opened by NcoI and BglII in restriction buffer NEB3.1 (NEB) and was subsequently column cleaned (Genelute PCR Clean-Up Kit, Sigma-Aldrich). RyhP-DNA fragments and pBAD were ligated (T4 ligase, NEB) and transformed in E. coli TOP10. After sequencing (eurofins), verified plasmids were transformed in E. coli O157:H7 EDL933. EHEC strains (containing either P1, P2, T2 or T3) were grown overnight in LB medium with a final concentration of 120 μg/ml ampicillin. The cell was density measured and both strains were mixed 50:50. Minimal Medium (MM) M9 without any iron added [36], but supplemented with a final concentration of 120 μg/ml ampicillin and 0.2% arabinose (for induction), was inoculated 1:1000 using the mixture and incubated 24 h at 37 °C with shaking at 150 rpm. Of both, the initial mixture and of the MM-culture, the plasmids were isolated and Sanger sequenced using the primer pBAD-C-R. The peak heights of the two nucleotides changed to form the stop codon in T2 or T3 were measured in comparison to the P variants, and the mean CI was calculated according to CI = (T(out) · P(in))/(P(out) · T(in)) [37] of P1 against, T2, P1 against T3 and P2 against T3. Given are mean and the standard deviations of three biological independent experiments.

Bioinformatics procedures

NGS mapping and evaluation

Raw data were deposited at the Gene Expression Omnibus [GEO: GSE94984]. Illumina output files (FASTQ files in Illumina format) were converted to plain FASTQ using FastQ Groomer [38] in Galaxy [38, 39]. The FASTQ files were mapped to the reference genome (NC_002655) using Bowtie2 [40] with default settings, except for a changed seed length of 19 nt and zero mismatches permitted within the seed in the Illumina data due to the short length of the footprints. Visualization of the data was carried out using our own NGS-Viewer [41] or BamView [42] implemented in Artemis 15.0.0 [43].

The number of reads was normalized to reads per kilobase per million mapped reads (RPKM) [44]. Using this method, the number of reads is normalized both with respect to the sequencing depth and the length of a given transcript. For determination of counts and RPKM values, BAM files were imported into R (R Development Team [45]) using Rsamtools [46]. For further processing, the Bioconductor [47] packages GenomicRanges [48] and IRanges were used [49]. The locations of the 16S rRNA and 23S rRNA are given by the RNT file from RefSeq [50]. findOverlaps of IRanges [49] was used to determine the remaining reads overlapping a 16S or 23S rRNA gene on the same strand. Reads from these rRNA-genes were excluded from further analysis as most rRNA had been removed using the Ribominus kit, as described above. countOverlaps can also determine the number of reads overlapping a gene on the same strand (counts). Using these counts, RPKM values were generated. For the value “million mapped reads”, the number of reads mapped to the genome, less the remaining reads overlapping a 16S or 23S rRNA gene, were used. Pearson correlation was calculated using Excel and Spearman rank correlation according to Wessa [51].

RCV thresholds

To distinguish between translated and non-translated for a given RNA, the ribosomal coverage value (i.e., reads of ribosomal footprints per reads of mRNA) was examined [52]. A negative control set contains the RCVs of tRNAs (“untranslated”). Sixteen phage encoded tRNAs, one tRNA annotated as a pseudogene, and one tRNA containing less than 20 reads in the combined transcriptome data set were disregarded since phage tRNAs sometimes have unusual properties [53, 54]. The RCVs of the tRNAs were transformed to ln(RCV), abbreviated LRCV. A density function $f\overset{}{\hat{\mkern6mu} }$ _LRCV-tRNA(x), with x = LRCV, was estimated by a kernel density estimation with Gaussian kernels and bandwidth selection according to Scott’s rule [55], furthermore a normal distribution was fitted as well for comparison. This was also conducted for the annotated genes (i.e., “translated” set), excluding zero RCVs (261 genes). To test the hypothesis “the RCV of the RNA belongs to the tRNA distribution”, we used the estimated tRNA LRCV distribution to compute a P value for an observed ncRNA with LRCV x as

$$ Pval\ (x)={\displaystyle {\int}_x^{+\infty }{\mathrm{f}}_{\mathrm{LRCV}\hbox{-} \mathrm{tRNA}}\ (x) dx,} $$

where we numerically evaluate the density function. For example, the hypothesis will be rejected for α = 0.05 for any x ≥ −1.816817 which corresponds to an RCV of 0.162542. Similar for α = 0.01 we obtain an RCV of 0.354859. For α = 0.05 we reject 52 of 115 annotated ncRNAs to be not translated, and for α = 0.01 we reject 63.

Since the interpretation of the results depends on the assumed distribution, we also used, at least for tRNAs, a fit of the normal distribution. The tails of the normal distribution tend to zero faster than before, which results in different P values. For example, for α = 0.05 a corresponding RCV of 0.646079 is obtained and for α = 0.01 the bound for the RCV is 0.928702. However, the normal distribution has no good fit (not shown) and is henceforth excluded.

In a similar way as for the tRNAs, we can use the gene distribution to test the hypothesis “the RCV of the RNA belongs to the mRNA distribution” by using the RCV of all annotated genes (aORFs) as a negative control set. In this case, the P value is computed by

$$ Pval\ (x)={\displaystyle {\int}_{-\infty}^x{\mathrm{f}}_{\mathrm{LRCV}\hbox{-} \mathrm{aORF}}(x) dx.} $$

For the latter function, we obtained the bounds 0.532837 and 0.197320 for α = 0.05 and α = 0.01, respectively. Thus, all RNAs above those values might be considered mRNAs.

Examination of known and novel ncRNAs

Escherichia coli O157:H7 EDL933 (genbank accession AE005174) contains 329 known ncRNAs (Rfam database, April, 30^th 2014 [56]). All ncRNAs which should naturally have ribosomal footprints (e.g., are leader peptides, riboswitches (several contain a translatable ORF [57]), occur within genes on the same strand, or tmRNA) were excluded from the analysis, as well as rRNAs and tRNAs. Thus, the excluded RNAs are 5S_rRNA (8x), ALIL (19x), Alpha_RBS, C4, Cobalamin, cspA (4x), DnaX, FMN, greA, His_leader, IS009 (3x), IS102 (2x), iscRS, isrC (2x), isrK (2x), JUMPstart (3x), Lambda_thermo (2x), Leu_leader, Lysine, Mg_sensor, mini-ykkC, MOCO_RNA_motif, nuoG, Phe_leader (2x), PK-G12rRNA (7x), QUAD_2, rimP, rncO, rnk_leader, rne5, ROSE_2, S15, SECIS (3x), SgrS, ssrA (tmRNA), sok (10x), SSU_rRNA_archaea (14x), STnc40, STnc50, STnc370, t44/ttf, Thr_leader, TPP (3x), tRNAs (99x), tRNA-Sec, Trp_leader, and yybP-ykoY. The remaining 116 RNAs were grouped in translated, non-translated and undecided according to their RCV. Translated ncRNAs were three-frame translated and proteins sequences were searched against the non-redundant database “nr” of genbank using blastp [58]. Cases in which the ORFs of the ncRNA generated a single hit to the database were excluded since a false annotation of the hit is likely for those.

In order to provide an initial in silico characterization of the putative function for the novel intergenically-encoded ncRNAs, we used CopraRNA [59, 60] and examined the functional enrichments returned for the predictions. CopraRNA was called with default parameters for each set of putative ncRNA homologs. To find ncRNA homologs for the CopraRNA prediction, GotohScan (v1.3 stable) [61] was run with an e value threshold of 10⁻² against the set of genomes listed in the Additional file 2: Table S2. The highest scoring homolog (i.e. having the lowest e value) for each organism was retained, if more than one GotohScan hit was present.

Ka/Ks ratio

The most likely ORF encoding a peptide was chosen according to the RIBOseq data. Homologs were searched using NCBI Web BLAST in the database nr using blastn. Hits with the highest e value but still achieving 100% coverage and displaying no gaps in the alignment were chosen (Additional file 3: Table S3). Gene pairs were examined using the KaKs_Calculator 2.0 [62] providing a number of algorithms which are compared and evaluated.

Shine-Dalgarno prediction

For any novel ncRNA with a significant blastp hit (e value ≤ 10⁻³, see above), a start codon (ATG, GTG, TTG) of the respective frame was searched closest to the start position of the ncRNA (except sgrS for which the start codon position is known, but ATG in E. coli K12 corresponds to ATT in EHEC, a rare but possible start codon; see Discussion). The maximum distance allowed between the ncRNA start coordinate and proposed start codon was ±30 bp. The region upstream of the putative start codon was examined for the presence of a Shine-Dalgarno sequence (optimum taAGGAGGt) according to [63] and [64]. A Shine-Dalgarno motif was assumed to be present at a ΔG° threshold of ≤ −2.9 kcal/mol (according to [63]) to allow weak Shine-Dalgarno sequences to be reported since even leaderless mRNAs exist [65].

For global examinations, we used PRODIGAL bins of the Shine-Dalgarno sequence and their distance to the start codon (Additional file 4: File S1) according to Hyatt et al. [66]. Bins without genes were omitted, and bins containing less than 100 genes were combined to superbins: S0, S2-3-4, S6, S7-8-9-12, S13, S14-15, S16, S18-19-20, S22, and S23-24-26-27 containing 629, 115, 116, 133, 1095, 664, 1191, 145, 687, and 327 genes, respectively.

Results and discussion

Sequencing statistics and footprint size

Two biologically independent replicates were used to assay reproducibility (Additional file 5: Figure S1). The numbers of footprint reads per gene of both RIBOseq replicates have a Pearson correlation of 0.86 and a Spearman rank correlation of 0.92, which was found to be slightly less compared to other NGS experiments [17, 67]. Nevertheless, the data sets were combined to increase the overall sequencing depth. In summary, 32.0 million transcriptome reads and 20.6 million translatome reads could be mapped to the EHEC genome (NC_002655; see Additional file 6: Table S4). Interestingly, the percentage of tRNA, an RNA species not translated, in both experiments was quite different. In the transcriptome, tRNAs contributed 31% of the library, whereas in the footprint libraries, tRNAs contributed only 0.3%. Such a difference is expected, since in the transcriptome sequencing, the tRNAs are processed together with the total RNA isolated. In contrast, in translatome sequencing, only translated RNAs are sequenced since the RNase digestion will destroy any RNA outside the ribosomes, including most tRNAs. However, some tRNAs might be trapped in the ribosomes and are recorded despite the RNase treatment. Thus, we reasoned that tRNAs would represent the best maximum background value for any carry-over of a non-translated RNA in the translatome sequencing.

The number of nucleotides which are protected by the ribosomes, i.e., the size of the footprints, was reported to be 28 nt in prokaryotes as well as in eukaryotes [20, 31, 32, 34, 68, 69]. Additionally, other studies using ribosome profiling in eukaryotes were able to determine the ribosome position of the footprints at sub-codon resolution (e.g. [70, 71]). The situation is quite different in bacteria: In one of the first studies in bacteria, Li et al. [72] determined the footprint size to range between 25 and 40 nt. Based on these results, O’Connor et al. [73] suggested that the footprint size may vary due to different progression rates of the ribosome. However, the enzyme used to obtain the bacterial ribosomal footprints in these studies was micrococcal nuclease which is known to prefer sites rich in adenylate, deoxyadenylate or thymidylate, which explains the varying length of the footprints [72]. In our study, after sequencing E. coli ribosomal footprints, the major peak of fragment sizes was observed at 23 nt, even despite the size-selection targeting 28 nt. We believe that RNase I, which we used, is a better choice [74, 75]. We also tested a number of commercially available RNases and mixtures of endo- and exo-cutting enzymes and received a consistent footprint size of about 23 nt and not 28 nt (unpublished data). The observed value of 23 nt may be explained by the different size of prokaryotic and eukaryotic ribosomes. Klinge et al. [76] estimated the mass of ribosomes to be 3.3 MDa for the eukaryotic and 2.5 MDa for prokaryotic, respectively. Assuming a roughly proportional scaling between the mass of the ribosome and its diameter suggest a bacterial footprint size of about 23 nt.

Putative novel ncRNAs with low ribosomal coverage

The ribosome coverage value (RCV) gives the ratio of RPKM footprints over RPKM transcriptome. ncRNAs should have low RCVs. The RCV is similar to the “translational efficiency” applied for eukaryotes [77] to determine the translatability of a given mRNA. The RCV varied between zero (for 261 annotated genes) and a maximum value of nearly 39 for an annotated gene. Low or zero RCVs for annotated genes can be explained by the internal status of the cells controlling translation independent of transcription. For instance, some mRNAs are blocked by riboswitches or bound by ncRNA (e.g. [78]). We examined the genes with zero reads in some detail. This group contains about 3-times more phage associated genes compared to all genes (36% versus 13%). The genes are shorter compared to all (about half the size) and a larger fraction is annotated as hypothetical (50% compared to 30% in the annotation NC_002655). We looked for transcription under any of 11 different growth conditions [17] and found transcription for less than 20% of those genes under any condition. However, the other genes might be activated in specific circumstances not tested yet. This is corroborated by our findings that some genes were induced when EHEC was grown in co-culture with amoeba (unpublished results), but are not activated in any other condition of the published data set [17].

To analyze the data for novel ncRNAs, the transcriptome data was analyzed for contiguous transcription patterns (no gaps allowed) containing at least 20 transcriptome reads which do not correspond to an annotated gene (i.e., in a distance of more than 100 nt to a same-strand annotated ORF of a gene). Start and end of the novel ncRNAs were defined as the first and last nt of the contiguous read pattern. The chosen value of 20 reads was applied independently of any length restriction. For a 100-bp transcript in our dataset this approximately corresponds to an RPKM of 20, which is about 200-times above background level for transcriptome sequencing [17].

Each novel transcript was analyzed for its RCV to determine whether it is potentially translated. As a negative control, we chose tRNAs which have RCVs in a range between 0.000173 and 0.094843. While the RCVs are small for tRNAs, the ratio between the highest and lowest RCV of the tRNAs is about 500-fold. We surmised that tRNA abundance might correlate either to the RCV or to the codon usage of EHEC (which correlates with tRNA abundance). However, no relationship was found (not shown) and the reasons for the difference in RCV remain unknown. For convenience, the RCV is shown as ln(RCV) (=LRCV) in Fig. 1. Figure 1a shows a histogram of the LRCV of tRNAs together with an estimated density function $\hat{f}$ _LRCV (x) obtained by a kernel density estimation (blue line). Next, the LRCV distribution of the annotated genes is shown in Fig. 1b (green line). Finally, Fig. 1c shows the LRCV of all annotated ncRNAs (red line; less those known to be translated; see Table 1). To determine, whether the RCV of a given RNA belongs either to the tRNA distribution group or the gene distribution group, we determined the lower and upper limit of the RCV corresponding to a probability of error of 99% (α = 0.01), respectively (see Methods). Below the RCV threshold 0.197 a transcript is considered to be untranslated and above 0.355 it is considered to be a candidate for translation. Thus, a transcript is qualified as a putative novel ncRNA only, if its RCV was below the lower threshold.

Fig. 1

Logarithmic (ln) ribosomal coverage (LRCV) of tRNAs, annotated genes, annotated ncRNAs and a merger of the former. a Histogram of the LRCVs (X-axis) of the tRNAs together with either the estimated density function (blue curve). The density of the individual tRNAs is shown as little blue bars on top of the X-axis. b LRCV histogram as before, but of the annotated genes and their estimated density function (green). c LRCV histogram as before, but of the known ncRNAs (see Table 1) together with their estimated density function (red). d A combination of the estimated density functions for the tRNAs (blue), the annotated genes (green) and the ncRNAs (red) of the former panels, shown a substantial overlap between the annotated genes and the ncRNAs supposedly non-coding

Table 1 Transcriptome and translatome profiles of 115 ncRNAs known from E. coli O157:H7 EDL933

Using the RCV limits mentioned in the methods section (i.e., RCV <0.197), 150 putative ncRNAs were discovered of which three examples are shown in Fig. 2. All novel ncRNA candidates are listed in Table 2, including the read counts, RPKM values and RCV values for each transcript. The putative novel ncRNAs range between 27 and 268 nt with an average size of 77 nt. One (ncR3609372) had a match in the Rfam database [56] as being a tRNA. We analyzed these transcripts to see whether they contained a potentially protein coding ORF. Of the 150 identified transcripts, 44 do not contain any ORF at all and only a minority of 6 candidates contains a putative ORF coding for more than 30 amino acids, indicating that most transcripts identified are truly non-coding. This agrees with the fact that all RCVs are below the threshold for translation. The RPKM-transcriptome values of the novel ncRNA transcripts range between 8 and 8857, the average being 198 (Table 2).

Fig. 2

Three examples of novel ncRNAs detected using transcriptome and translatome analysis. A genomic area is visualized in Artemis 15.0.0 [43]. In the lower part of the panels, the genome (shown as grey lines) is visualized in a six-frame translation mode. Numbers given between the grey lines indicate the genome coordinates. On top of the forward strand are three reading frames and on the reverse DNA strand are three further reading frames. Each reading frame represented is visible by the indicated stop codons (vertical black bars). Annotated genes are shown in their respective reading frame (turquoise arrows) and also on the DNA strand itself (white arrows). The gene name is written below each arrow. Any protein-coding ORF must be at least located between two black bars, with the downstream stop codon being the translational stop. In the upper part of the panels, the DNA is indicated by a thin black line and the sequencing reads matching to the forward or reverse strand are shown above or below this line. The sequencing reads from the footprint (yellow line) and transcriptome (blue line) sequencing are shown as coverage plot, respectively. The pink shaded area in the coverage plot corresponds to the novel ncRNAs, which are drawn in by red arrows. Novel ncRNAs were identified by their very low RCV, thus, hardly any footprint reads (in yellow) but a number of transcriptome reads (in blue; see Table 2). Known ncRNAs are indicated on the DNA by a bright green arrow. Since ncRNAs supposedly do not contain a protein-coding ORF, these genes are only shown on the DNA. a ncR3665651. b ncR3690952. c ncR1085800

Table 2 Novel non-coding RNA (ncRNA) candidates (150 in total) based on transcriptome sequencing and ribosomal profiling. ncRNAs are identified by their start position on the genome given in the name (abbreviated as ncR#)

Presence of novel ncRNAs in E. coli K12

In E. coli O157:H7 EDL933, 329 ncRNAs have been annotated [2], but various bioinformatic studies suggest the existence of up to 1000 ncRNAs in E. coli (e.g. [8–11]) and probably in other bacteria as well (e.g. [19, 79]). Our current study presents even under a single growth condition 150 new ncRNA candidates. For comparison, we determined the presence of corresponding regions in the E. coli K12 strain MG1655. We found 102 of 150 novel ncRNAs regions present in MG1655. Next, we searched data of prokaryotes having both, transcriptome and translatome data of the same experiment. Only a single study was published by the Weissman group of MG1655 grown in MOPS glucose medium [80]. In addition, the ArrayExpress database contains a further dataset of MG1655 grown in LB (E-MTAB-2903). In MOPS medium with glucose at OD 0.3 and in LB medium at an OD of about 0.5, 43 and 66 of the 102 putative ncRNAs were found to be transcribed in MG1655, respectively. Combining both datasets confirmed transcription (without translation) of 74 of the 102 ncRNAs under either condition in E. coli MG1655 (Additional file 7: Table S5).

Detection of ncRNAs by Northern blots

To verify the existence of at least some annotated ncRNAs, Northern blot analysis was conducted for five of the annotated ncRNAs of different length and strength. Three were verified, namely ffs, sraJ, and STnc100_4 (Table 1 and Fig. 3). We then chose seven exemplary novel ncRNAs to be confirmed using Northern blots. However, of the novel ncRNAs only the two transcripts with the highest RPKM in the transcriptome of 8857 and 6404 could be verified as sum signal since they are indistinguishable on the basis of Northern blots (Fig. 3). Obviously, Northern blots have a certain detection limit. Under the conditions applied in this study, any RNA required an RPKM value of about 2000 to be detectable. RNAs transcribed at lower levels were not detected via hybridization. A sufficiently high number of RNA molecules are needed to generate a signal passing the detection threshold, a problem also common to microarrays [81, 82].

Fig. 3

Detection of novel and annotated ncRNAs by Northern blots. Since ncRNAs do not have defined ends like, e.g., ORFs which have start and stop codons, their actual length may differ somewhat from the expected length (compare to Table 1). The contrast of the bands has been adjusted by gamma correction using digital image processing for better visibility. a ncR1085800 and ncR1481381. Both ncRNAs are indistinguishable by their sequence. b STnc100_4. c Bacteria_small SRP/ffs. d GlmZ_SraJ_2

Putative functions and differential expression of the novel ncRNAs

To examine putative functions of the novel ncRNA candidates, we used the sRNA target prediction tool CopraRNA [59, 60]. For 14 of the 150 novel ncRNAs, a significant functional enrichment was found (Table 2). The targets include a diversity of metabolic and regulatory functions within the cell, e.g., synthesis pathways of amino acids and vitamins; but also respiratory functions and oxidation of components etc.

Interestingly, 121 of the novel ncRNAs were found to be expressed (i.e., ≥ 10 RPKM, which is ≥ 100-fold above background) in the data of a former study [17], when grown in eleven different growth conditions for at least one condition. Forty-six novel ncRNAs revealed 4-fold differential expression in at least one other condition when compared to plain LB. Example data are given in Table 3, the full data set can be found in the Additional file 7: Table S5. Combining these findings of CopraRNA predictions (14), Rfam match (1), expression (121), and regulation (46) suggests that at least 126 out of 150 putative ncRNAs are not just a random by-product of pervasive transcriptional activity, but might fulfill specific functions in the cell.

Table 3 Expression of exemplary novel ncRNAs under 11 different growth conditions (MM, minimal medium)

Evidence for translation of annotated ncRNAs

To our own surprise, a significant number of annotated ncRNAs had high RCVs indicating translation, which we examined further. Table 1 shows the known ncRNAs which i) are independent from protein coding genes (i.e., are not leader peptides or riboswitches, etc.), ii) are not ribosomal RNA or iii) do not encode tRNAs. The remaining 115 annotated ncRNAs were categorized according to their RCV (Fig. 1c; Additional file 8: Table S6). As expected for ncRNAs, 52 of these ncRNAs are not translated and have a low RCV (RCV ≤ 0.16). This indicates transcription but no translation. Surprisingly, we identified 52 ncRNAs with an RCV higher than 0.355 (α = 0.01) which we used as lower limit for considering a transcript to be translated (Additional file 9: Figure S2). For both cases, an ncRNA example with low (csrB) and high (arcZ) RCV is shown in Fig. 4. Eleven ncRNAs fall in an RCV range above the upper limit for untranslated and below the lower limit for translated RNAs and, thus, their translation status (i.e., either untranslated or weakly translated) is difficult to assess. In summary, the ncRNAs were divided into three groups with different ribosome coverage: low RCV similar to untranslated RNAs (52 or 45.5%), such of ambiguous nature (11 or 9%), and those with high RCV similar to translated genes (52 or 45.5%). Clearly, the RCV threshold at which an RNA is considered to be translated depends on the assumed distribution fitted to the tRNA values (see Methods). In any case, different thresholds only alter the region of uncertainty, but do not invalidate our principal finding that quite a number of annotated ncRNAs appear to be associated with ribosomes. Normally, translation is considered the main cause for ribosome binding of an RNA in RIBOseq experiments [83].

Fig. 4

Visualization of ribosomal footprints and transcript reads mapping to annotated ncRNAs as coverage plots. A genomic area is visualized in Artemis 15.0.0 [43]. In the lower part of the panels, the genome (shown as grey lines) is visualized in a six-frame translation mode. Numbers given between the grey lines indicate the genome coordinates. On top of the forward strand are three reading frames and on the reverse DNA strand are three further reading frames. Each reading frame represented is visible by the indicated stop codons (vertical black bars). Annotated genes are shown in their respective reading frame (turquoise arrows) and also on the DNA strand itself (white arrows). The gene name is written below each arrow. Any protein-coding ORF must be at least located between two black bars, with the downstream stop codon being the translational stop. In the upper part of the panels, the DNA is indicated by a thin black line and the sequencing reads matching to the forward or reverse strand are shown above or below this line. The sequencing reads from the footprint (yellow) and transcriptome (blue) sequencing are shown as filled coverage plots, respectively. The known ncRNAs are indicated on the DNA by a bright green arrow. Since ncRNAs supposedly do not contain a protein-coding ORF, these genes are only shown on the DNA. a csrB: Very few footprint reads are seen for CsrB, indicating that this ncRNA is not translated. b arcZ: In contrast, ArcZ is covered with many footprints and a number of transcript reads are found. All further examples are shown in Additional file 9: Figure S2

We analyzed the potential ORFs of the 52 ncRNAs covered by ribosomes for their annotation status in other organisms using blastp [58]. Twenty were found to contain ORFs which achieve blastp-hits to multiple genes annotated in other enterobacteria (e value 10⁻³ or lower), mainly in other Escherichia coli strains. From these, 15 are annotated as hypothetical proteins, two belong to toxin-antitoxin systems, one encodes a conserved domain of phage origin and the remaining two are membrane proteins (Additional file 8: Table S6).

Correlation of translation with Shine-Dalgarno sequences

The presence or absence of a Shine-Dalgarno sequence in proper distance to the start codon can be an indicator for a translational start [66]. A strong Shine-Dalgarno sequence should correspond to a high RCV. On a global scale, i.e. taking average values of all genes with comparable Shine-Dalgarno sequences, such a correlation was found (Additional file 4: File S1). However, predictions are unreliable for single genes. Since several genes exist which either have no Shine-Dalgarno or are completely leaderless [65], a missing Shine-Dalgarno is not necessarily an indication for absent translation. We then searched for the presence of a Shine-Dalgarno sequence for those 20 ncRNAs which have a blastp hit. A start codon in reasonable distance to the start coordinate of the ncRNA was selected (see Methods) and a possible Shine-Dalgarno sequence was determined according to Ma et al. [63], also including weak Shine-Dalgarno sequences. In 11 of 20 cases, a putative Shine-Dalgarno sequence was found (Table 1, Additional file 8: Table S6). The Shine-Dalgarno sequences were also determined as above according to Hyatt et al. [66] (see Additional file 4: File S1), but this method is more stringent and misses some of the weaker sequences (Additional file 8: Table S6). The observation that 11 out of 20 translated ncRNAs with blastP hit (i.e., 55%) have Shine-Dalgarno sequences compares well to about 57% annotated genes possessing such a sequence in E. coli K12 [64].

Why are ncRNAs covered with ribosomes?

Translational profiling showed that 52 annotated ncRNAs have high RCVs. High RCVs may occur due to incomplete digestion of free RNA. Therefore, we had performed two rounds of RNase I digestion and sucrose density gradient centrifugation for ribosomal enrichment, which makes this assumption very unlikely. Most ncRNAs are reported in the Rfam database to bind Hfq and regulate via antisense pairing to their target genes; some ncRNAs are of completely unknown function, and few are involved in toxin-antitoxin interactions. We consider it unlikely that the high numbers of footprints are false-positives in all cases. While the phenomenon of “translated ncRNAs” is highly discussed for eukaryotes [70, 71, 84–89], this observation has, to our knowledge, only rarely been reported for bacteria, i.e. SgrS/SgrT or the “ncRNA” C0343 ([90]; see below, [91]).

In any case, the ribosomal “coverage” of tRNAs (median RCV 0.03), taken as background in this study, is far below the high ribosomal coverage of some ncRNAs. Finally, another explanation for ribosomal coverage of ncRNAs is regulatory functions performed by interaction of the ncRNA with the ribosomes and, thereby, causing accidental carry-over. However, ribosome-interacting ncRNAs are a minority according to Guttman et al. [86].

RNAs functioning as both ncRNA and mRNA?

A few ncRNAs which are also translated have been suggested to exist in bacteria and are termed coding non-coding RNAs (cncRNAs) [24]. sgrS/sgrT is the only known example for E. coli K12 [90]. In EHEC EDL933, the ATG start codon used by E. coli K12 is mutated to ATT. In addition, the Shine-Dalgarno sequence has changed from AAGGGGGT in K12 to AAGGAGGT in EDL933, the very best category S27 of Hyatt et al. [66]. Since a strong SD sequence compensates a weak start codon [63], and sgrS has an RCV of 1.55 (Table 1), and ATT is known to be a (very rare) start codon in E. coli [92–94], we hypothesize that EHEC synthesizes SgrT using the uncommon start codon ATT. Interestingly, the ORF encoding for SgrT gave a Ka/Ks ratio below 1, i.e. 0.15 with a P value of about 0.002. Unfortunately, most ORFs found covered with footprints proved to be too short for any meaningful Ka/Ks analysis (data not shown). Only one other footprint-covered ORF of the ncRNA MicA gave significant results. This ORF had a Ka/Ks ratio of about 0.35 with a P value of about 0.018 (Additional file S3: Table S3).

Not all former entities named as ncRNA in the past, however, are cncRNAs. For instance, C0343 had formerly been described as ncRNA, but contains an ORF and yields an RCV of 2.49 in our study (not shown). This validates Washietl et al. [91] who shows that C0343 encodes a short 57-aa protein. Consequently, this entity was possibly falsely labelled as ncRNA and it had been removed from the Rfam database. However, a former study described 72 novel intergenic small protein-coding genes of EHEC [83]. We found six instances in which the locus of a novel protein-coding gene overlaps fully or partly with the locus of one of the ncRNAs (Additional file 8: Table S6), which also hints towards cncRNAs.

In any case, we suggest being cautious in labeling any ribosome covered “ncRNA” of E. coli found in this study as cncRNA since further experimental evidence is needed. Based on our current results, we conclude that ribosome covered ncRNAs may represent a mixture of misannotated short mRNAs, ncRNAs with a regulatory function including potential ribosomal binding, and cncRNAs translated indeed. To corroborate this hypothesis about additional cncRNAs and to confirm the existence of novel peptides from so called “non-coding” RNAs as indicated by ribosomal footprints, we tested the footprint-covered ORF of ryhB for a phenotype (see below).

ryhB supposedly is a novel cncRNA, encoding the RNA RyhB and a phenotype-causing peptide, RyhP

Closer examination of footprint signals for several ncRNAs revealed possible ORFs which encode novel peptides. We chose ryhB for further examination, since the encoded RNA-molecule RyhB has a well-known function in iron homeostasis for many bacteria [95–97]. Accordingly, we expected iron-limiting to be the most likely condition in which a phenotype for this novel peptide might be found. Thus, we picked the best matching ORF according to the RIBOseq data, coding for the nona-peptide MAHIASSIT (Fig. 5; start codon ATT) and named it ryh B-encoded peptide, RyhP, in the following. This ORF was introduced on a high-copy arabinose-inducible plasmid in EHEC wild type. In cloning, we omitted all non-coding parts of ryhB, to limit any effect the expressed (m)RNA-fragment might have (sequence P1). To even further reduce the possibility that the expressed RNA and not the peptide itself causes the phenotype, we changed all codons of the ORF such that the same peptide is produced, but the underlying RNA sequence differs maximally from the wild type sequence (P2). This strategy prevents the RNA made hybridizing with any natural target RNAs [e.g., 99]. Two negative controls were created, either with the second (T2) or third codon (T3) changed into a stop codon, terminating RyhP translation prematurely. Competitive indices (CI) under RyhB-inducing condition (i.e. low-iron) showed a significant advantage of the strain possessing the RyhP-producing plasmid over those strains containing a plasmid with stop codons in the RyhP-ORF (Table 4).

Fig. 5

Visualization of individual ribosomal footprints mapping to rhyB. The genomic area is visualized in Artemis 15.0.0 [43]. In the lower part of the panels, the genome (shown as grey lines) is visualized in a six-frame translation mode. Numbers given between the grey lines indicate the genome coordinates. On top of the forward strand are three reading frames and on the reverse DNA strand are three further reading frames. Each reading frame represented is visible by the indicated stop codons (vertical black bars). Annotated genes are shown in their respective reading frame (turquoise arrows) and also on the DNA strand itself (white arrows). The gene name is written below each arrow. In the upper part of the panels, the DNA is indicated by a thin black line and the footprint reads (blue) matching to the forward or reverse strand are shown above or below this line. The shaded areas indicate ryhB (pink), the coding ORF RyhP (green) and a putative weak Shine-Dalgarno sequence (brown; ggagaa)

Table 4 Competitive index values (CI) for EHEC strains possessing a wild-type like ORF encoding RyhP (P1 or P2) or an ORF with a premature stop codon (T2 or T3) plusminus their standard deviations (Std)

RyhB folds when not bound to its regulated target RNA (Fig. 6) and this, assumedly, makes the coding ORF unavailable for translation. However, ribosomes are able to resolve secondary structures of mRNAs [98]. Furthermore, RyhP has a weak putative Shine-Dalgarno motif (i.e., ggagaa) upstream. Upon binding a target mRNA like sodA [99], the RNA structure opens and the Shine-Dalgarno sequence is set free (Fig. 6). If this opening facilitates ribosomal binding for translation initiation of the RyhB RNA, and subsequent progression of ribosomes along the RNA, must remain open.

Fig. 6

Overview of the secondary structures formed by RyhB for the molecule on its own (top) and after binding to a target RNA, like sodA (bottom). Structures are taken from [99]. Individual bases have been highlighted. Underlined, putative Shine-Dalgarno sequence; green, start codon; violet/orange, individual codons along the frame; red, stop codon, bold; bases involved in hybridization to the sodA-target

Conclusion

In the past, very short proteins or peptides were excluded from annotation and believed to be unlikely. Some short mRNAs could have been labeled as ncRNA solely on this presumption. However, more and more small proteins are being discovered. For instance, a number of small genes have been described for E. coli in recent years. These genes were hard to detect because they appear to be membrane proteins and are induced under stress conditions only [100, 101]. In another study, we confirmed the existence of 72 novel and short protein-coding genes in the EHEC genome, some which were verified by proteome data [83]. Similar findings have been made by other groups (see, e.g., [102–104]), and future research could confirm the existence of more of these proteins similar to studies conducted in eukaryotic ribosomal profiling [70, 105–107].