Whole-genome mapping of 5' RNA ends in bacteria by tagged sequencing: a comprehensive view in Enterococcus faecalis.

Enterococcus faecalis is the third cause of nosocomial infections. To obtain the first snapshot of transcriptional organizations in this bacterium, we used a modified RNA-seq approach enabling to discriminate primary from processed 5' RNA ends. We also validated our approach by confirming known features in Escherichia coli. We mapped 559 transcription start sites (TSSs) and 352 processing sites (PSSs) in E. faecalis. A blind motif search retrieved canonical features of SigA- and SigN-dependent promoters preceding transcription start sites mapped. We discovered 85 novel putative regulatory RNAs, small- and antisense RNAs, and 72 transcriptional antisense organizations. Presented data constitute a significant insight into bacterial RNA landscapes and a step toward the inference of regulatory processes at transcriptional and post-transcriptional levels in a comprehensive manner.

INTRODUCTION

Enterococcus faecalis is a ubiquitous Gram-positive bacterium and one of the first colonizers of the human gastro-intestinal tract after birth. It belongs to the core-microbiota and lives in the guts during the entire human life, suggesting a contribution of the bacterium to intestinal homeostasis (Campeotto et al. 2007; Adlerberth and Wold 2009; Qin et al. 2010). In contrast to this potentially beneficial role, E. faecalis is also the third cause of nosocomial infections and may carry and transfer various antibiotic resistances to other bacterial species, making its presence in the medical environment a serious concern (Arias and Murray 2012). The opportunism of E. faecalis, i.e., the transition from commensalism to pathogenicity in response to environmental cues, underlines its capacity to adapt to and survive in harsh conditions. Thus, deciphering the regulatory pathways that enable E. faecalis to undergo the transition from commensalism to pathogeny is a key component in the understanding the dual lifestyle of this microorganism (Gilmore and Ferretti 2003). The V583 strain was one of the first discovered vancomycin-resistant clinical isolates of E. faecalis (Sahm et al. 1989). Its genome, a circular chromosome (3218 kbp) and three circular plasmids pTEF1 (66 kbp), pTEF2 (57.7 kbp), and pTEF3 (18 kbp), contains at least 3264 annotated protein-coding genes (Paulsen et al. 2003). Although partial transcriptomic analyses have been performed (Vebø et al. 2009, 2010; Aakra et al. 2010; Opsata et al. 2011), a comprehensive and dynamic view of the RNA landscape of V583 is missing.

Whole-transcriptome studies of prokaryotes via tiling arrays and RNA sequencing (RNA-seq) have unveiled a plethora of actively transcribed RNAs, and highly complex transcriptional organizations due to numerous promoters nested in open reading frames (ORFs), antisense (asRNAs) and small RNAs (sRNAs) genes (among other reviews Toledo-Arana and Solano 2010; Georg and Hess 2011). These global studies have been extremely valuable in determining transcribed regions in bacterial genomes and comparing RNA expression levels between different physiological states. Yet, their functional and regulatory insights remain incomplete as primary and processed RNAs cannot be distinguished and hence transcriptional (RNA synthesis) and post-transcriptional processes (RNA processing and stability) cannot be separated. The use of differential RNA-seq (dRNA-seq), an astute method that enriches an RNA population for primary transcripts, partially overcomes these limitations and gives access to the primary transcriptome (Albrecht et al. 2010; Bohn et al. 2010; Irnov et al. 2010; Sharma et al. 2010; Conway et al. 2014; Thomason et al. 2015). Yet, a major limitation of dRNA-seq is that all transcripts cannot be detected in a single experiment as they are degraded by a 5′-phosphate-dependent exonuclease, and thus information on post-transcriptional events is lost (Sharma et al. 2010). Global scale analysis of RNA stability has been performed in a few bacterial species, e.g., Bacillus cereus (Kristoffersen et al. 2012), Bacillus subtilis (Hambraeus et al. 2003), Escherichia coli (Selinger et al. 2003; Mohanty and Kushner 2006; Esquerré et al. 2013), Mycobaterium tuberculosis (Rustad et al. 2013), Lactococcus lactis (Redon et al. 2005), and Prochlorococcus (Steglich et al. 2010). These stability analyses have highlighted the broad and crucial contribution of RNA stability to gene expression reprogramming when bacteria face stresses, adapt to novel nutrient conditions, or grow at different rates. Those analyses consider transcribed regions as unique entities, where different sorts of RNA molecules can be present and cannot be seen.

We previously described a method that enables us to differentially tag 5′ ends of primary and processed RNAs (Fouquier d'Hérouël et al. 2011). In the present work, we have coupled this method to RNA-seq, yielding novel insights into the bacterial RNA landscape where the primary and the processed transcriptomes (ppRNome) are unveiled within a single experiment. We have sorted transcription start sites (TSSs) and processing sites (PSSs) and validated the method by reproducing known results for E. coli. The presented data provide a first snapshot of the RNA landscape of the human pathogen E. faecalis.

RESULTS AND DISCUSSION

Global view of the E. faecalis RNA landscape

Bacterial native (or primary) transcripts undergo cleavage that can be maturation or degradation processes (Rochat et al. 2013). Without the ability to identify and discriminate primary from processed transcripts, we have only partial information of gene expression control at the genome scale.

We have developed a method combining our previously introduced 5′ RNA end tagging (Fouquier d'Hérouël et al. 2011) with deep sequencing technologies and termed it “tagRNA-seq” to visualize the primary and processed transcripts of E. faecalis in a genome-wide manner (Fig. 1; Supplemental Table SA and “The ppRNome browser” website, see “Visualization of Results” in Materials and Methods).

FIGURE 1.

Three examples of 5′ RNA ends viewed by tagRNA-seq on the E. faecalis V583 chromosome from the “ppRNome” browser. Below the line “Annotated genes,” coordinates are those of the chromosome. The location of tags detected is indicated by the black vertical lines and the red arrows. TSS-tags are shown in the upper line, PSS-tags in the line below. “RNA levels” show the RNA signal detected; in red from St, in blue from Rt. Accurate values obtained for TSS- and PSS-tag counts and RNA levels are provided in Supplemental Table SA. (A) Transcription start site mapped at 769663/-5 for ef0809. This TSS could be easily predicted from the signal coverage. (B) TSS mapped at 1992494 for ef2071. This TSS is internal to the signal provided by the transcription of ef2072 and would be difficult to predict. (C) Processing site mapped at 1,121,951/-53 for the RNA RnpB. This PSS is a dozen nucleotides downstream from the previously mapped TSS (see “Processing sites”).

TagRNA-seq was performed on total RNAs extracted from bacteria grown in static (S) and respiratory (R) conditions, providing transcriptomes coined “St” and “Rt,” respectively (Supplemental Section S2, Table S1). In parallel to these, and as control, three other RNA libraries from E. faecalis and one from E. coli were sequenced on different next-generation sequencing platforms (see Materials and Methods; Supplemental Section S2). Globally, St and Rt show that significant transcription occurs in a limited portion of the E. faecalis genome. Out of the 3.34 Mbp long genome, 1.65 Mbp appears to be transcribed in each condition (coverage >2×), including 90 kbp due to antisense transcription and 470 kbp made up by nonannotated and/or noncoding portions, i.e., 5′- and 3′-untranslated regions (UTRs), unannotated ORFs, and as- and sRNAs (see below). These data are in line with previous reports highlighting that the information provided by annotations of bacterial genomes on their gene content remains incomplete (Toledo-Arana et al. 2009; Albrecht et al. 2010; Irnov et al. 2010; Sharma et al. 2010; Mitschke et al. 2011; Wurtzel et al. 2012; Conway et al. 2014; Thomason et al. 2015).

The calling of TSS and PSS from tagRNA-seq

We compared deep-sequencing data obtained with tagged and untagged RNA libraries prepared from E. faecalis grown in S conditions. Predictions of 5′ edges based on sequence coverage signals from both libraries show good agreement, indicating that the tagging procedure does not alter the allocated coverage nor the location of putative transcription starts (Supplemental Section S3, Figure S3; the ppRNome browser). Several TSSs mapped previously by other methods were retrieved by tagRNA-seq at near-identical locations (±2 bp) attesting to the reliability of the method. For example, we find the TSSs of sodA (ef0463), coding for the superoxide dismutase, ptb (ef1663), coding for the phosphotransbutyrylase, fsrB/D (ef1821), coding for the cysteine protease-like processing enzyme FsrB and the autoinducing propeptide FsrD of the fsr system, a homolog of the accessory gene regulator (agr) of Staphylococcus aureus, and gelE (ef1818), coding for a gelatinase (Qin et al. 2000, 2001; Ward et al. 2000; Nakayama et al. 2006; Verneuil et al. 2006) (see below; Supplemental Table SB).

In the ideal case, the procedure should identify unambiguously TSSs and PSSs. In practice, a fraction of 5′ ends attached to a TSS-tag were also ligated to a PSS-tag. Indeed, in vivo, 5′ triphosphate RNA ends are enzymatically converted to monophosphate, often as a first step of RNA degradation (Bail and Kiledjian 2009, and references therein). Therefore, a fraction of TSSs are expected to be associated with the PSS-tag. This effect may be further strengthened by spontaneous hydrolysis of 5′ triphosphate RNA ends during the ligation step of the PSS-adaptor and preceding RNA treatments, generating 5′ ends opened for ligation. On the other hand, the first step of the tagging procedure using the T4 RNA ligase is certainly not complete and acts with different efficiency on different RNA molecules (Zhuang et al. 2012; Raabe et al. 2014). Therefore, at the second ligation step, 5′ monophosphate ends (i.e., PSSs) that have escaped the first tag can be ligated to the TSS-adaptor and appear as false TSSs. For each 5′ RNA end mapped in this study, Figure 2A presents the number of each tag counted. The distribution of 5′ termini extends continuously between the two axes and hence does not give an immediate way to distinguish TSSs from PSSs. However, the distribution can be sorted by additional arguments, paying the price of discarding information on a fraction of mapped positions. (1) PSSs (i.e., 5′ monophosphate groups) for which the first ligation step was partial and also tagged with the TSS-tag sequence at the second step, should not give more TSS-tag counts than PSS-tag counts since the enzyme should act with the same efficiency on the same RNA end at each step. Therefore, points (i.e., 5′ RNA ends) above the diagonal may be either TSSs or partially ligated PSSs, but 5′ RNA termini falling below the diagonal in Figure 2A should be TSSs. Obviously, such a cutoff eliminates true TSSs that would exist in vivo mainly as 5′ monophosphate ends. (2) In accordance with the previous argument, all other TSSs known from the literature fall below the diagonal with one exception (Fig. 2B), the ncRNA Ref25C (RNA in 25C), which we discuss in more detail in Supplemental Section S4. (3) We considered separately 5′ edges of transcribed regions that feature an absence of detectable expression upstream and should therefore be a signature of a TSS. As expected for those selected RNA ends, and in accordance with the two first arguments, a clear distribution below the diagonal appears (Fig. 2C). (4) A motif search in DNA regions between 0–30 and 20–40 nt upstream of 5′ RNA ends located below the diagonal shows that >80% of them contain at least one canonical sequence featuring a promoter region (−10 and/or −35 boxes). In contrast, the same search performed for 5′ RNA ends above the diagonal does not retrieve any sequence reminiscent of a canonical promoter region (see below). The presence of promoter motifs in one area delineated by the diagonal is a very strong argument in favor of the location of true TSSs <45° in the plot presented in Figure 2A. Considering these rules and in order to err on the side of caution, in this work we will only consider points (i.e., 5′ RNA ends) below 30° as “TSSs,” and >60° as “PSSs”; for points in between, 5′ RNA ends cannot be assigned with certainty and will be considered as undetermined. Compared with other single nucleotide resolution RNA-seq methods, tagRNA-seq provides, for the first time, an accurate mapping of TSSs buried in transcribed regions and of RNA cleavage sites at a comprehensive scale in a single view, without requiring comparison between transcriptomes (Sharma et al. 2010; Nicolas et al. 2012; Wurtzel et al. 2012) (Fig. 1; the “ppRNome” browser; Supplemental Table SA).

FIGURE 2.

Scatter plot showing TSS-tag counts versus PSS-tag counts. (A) At each position of the genome. (B) At genomic locations within 2 bp of previously experimentally mapped transcription start sites. (C) At genomic locations within 2 bp of 5′ RNA edges of transcribed regions (see Materials and Methods; Supplemental Section S3). About 80% of predicted 5′ RNA ends fall below the diagonal.

Transcription start sites in E. faecalis

Within the area <30° in Figure 2A, we mapped a total of 559 TSSs on the V583 E. faecalis genome, combining both St and Rt (Supplemental Table SA,SB). A total of 327 TSSs were common to both transcriptomes. Among candidates classified as TSSs in St but not in Rt, 49 were classified as inconclusive due to a location between the 30° and 60° lines in Figure 2A, 1 was classified as PSS and 27 were inconclusive due to a weak (TSS-tag + PSS-tag) signal (i.e., <3.2× per million of tagged reads aligned). For the corresponding candidates in the Rt conditions, those numbers are respectively, 36 between the 30° and 60° lines, 3 classified as PSSs and 116 were inconclusive due to low signal in St (Table 1).

TABLE 1.

Summary of the numbers of TSSs reported in R and S conditions

Motif detection and promoter features in the E. faecalis genome

Up to date, <50 TSSs have been experimentally characterized in E. faecalis (Fouquier d'Hérouël et al. 2011, and references therein). In order to better define promoter regions in this species, we took advantage of our 5′ ends mapping and performed a blind search for common sequences nested in DNA regions preceding RNA extremities using the MEME suite (Bailey et al. 2009). By doing so, this search also enabled us to challenge our classification of 5′ RNA ends based on the tagging method as presented in Figure 2A. We defined four groups of DNA regions: two groups below the diagonal, one from 0° to 30° (called as TSSs), a second from 30° to 45° (called as undetermined, but expected to contain mainly TSSs), and two groups above the diagonal, one from 45° to 60° (called as undetermined but with a few TSSs, e.g., Ref25C) (see Supplemental Section S4), and a second group from 60° to 90° (called as PSSs). DNA sequences used as input for MEME and a detailed list of the motifs discovered are presented in Supplemental Table SC. For groups below 30°, the analysis reveals motifs with strong statistical significance (E-values below 10−30) and consensus sequences: Within the region [−30 … 0] and centered around position −9.7 ± 2.6, we found GnTATAAT, the canonical −10 box; in the [−40 … −20] region, the motif TTGACAA was found centered at −31.5 ± 2.3, the canonical −35 box. The −10 box appears with a high frequency (83.5%) and ends 5–9 bp from the 5′ RNA ends mapped. The −35 box was found in 20.6% of input sequences. At least 90% of the sequences where a −35 box is detected also have a canonical −10 box. Boxes defined as −10 and −35 are spaced by a 16–22 bp long sequence. Thus, the most significant motifs discovered correspond to the canonical −10 (TATAAT) and −35 (TTGACA) sequences of promoters recognized by the vegetative RNA polymerase loaded with the transcription initiation factor SigA (RpoD, σA or σ70) in the most studied bacteria E. coli and B. subtilis (Harley and Reynolds 1987; Helmann 1995). The presence and the location of −10 and −35 boxes on DNA regions upstream of 5′ RNA ends falling in the area defined by the angle between 0° and 30° in Figure 2A, reinforces our previous conclusion that these RNA extremities are TSSs. In line with this conclusion, for features with an angle between 0° and 45°, −10 and −35 canonical boxes are still the most frequently found motif but the numbers fall to 80.8% and 14.7%, respectively, which indicates that the density of true TSSs is indeed higher for signals corresponding to a low angle in the plot (≤30°). Within the two groups of sequences above 45°, the most significant motif discovered is AACGA/TAC/G A/G, found in <10% of sequences. To our knowledge, this purine-rich motif does not resemble any canonical sequence of bacterial promoter described previously. One might speculate that this sequence represents a frequent RNA motif targeted by an endoribonuclease, but further experiments will be required to confirm this hypothesis. Nonetheless, this observation reinforces our conclusion that the majority of TSSs do not locate >45° in Figure 2A.

In addition to SigA, three other sigma factors have been predicted in E. faecalis V583, SigH (Ef0049, the heat-shock factor), SigV (Ef3180, an “extracytoplasmic” factor) and SigN (Ef0782, a σ54-like factor) (Paulsen et al. 2003). ORFs coding for SigH and SigV are not expressed in S and R growth conditions (Supplemental Table SD; the ppRNome browser), hence we did not expect to find TSSs whose promoter regions would carry consensus sequences recognized by either one of these factors. In contrast, the sigN encoding sequence is transcribed and we sought manually for the consensus sequence of SigN-dependent promoter ahead of TSSs mapped (−24/−12; TTGCCACNNNNNTTGCT) (Buck et al. 2000; Héchard et al. 2001; Iyer and Hancock 2012). Only six corresponding locations were found across the whole genome: upstream ORFs coding for components of phosphor-sugar transfer systems (PTS), ef0019, ef1012, ef1017, ef1954, ef3210, and fabF-2 coding for an enzyme involved in fatty acid and biotin metabolism. Out of those six locations, the TSS for ef1012 is detected and a tag signal below our selection threshold is found upstream of ef1017.

Processing sites in E. faecalis

PSS-tags are found ∼50% more abundant than the number of total TSS-tags detected (Supplemental Table SA, Section S2). In contrast to TSS-tags that appear with a discrete distribution at 5′ edges or nested within transcribed regions, PSS-tags, in addition to colocalize with TSS-tags, also tend to spread out over RNA signals. Although we cannot rule out experimental RNA breaks, such a distribution of PSS-tags is expected as they label any type of 5′ monophosphate RNA ends, including processing sites, degradation products and hydrolyzed 5′ triphosphate ends. To pinpoint major PSSs within the transcriptome, we only considered 5′ ends located within the area delineated by the 60° angle in Figure 2A and above our acceptance threshold in both St and Rt. Ignoring rRNA and tRNA loci we mapped a total of 352 PSSs candidates (Supplemental Table SE).

Up to now, most of bacterial transcriptomic studies have focused on TSSs, RNA levels and the discovery of unannotated genes (e.g., Toledo-Arana et al. 2009; Sharma et al. 2010; Nicolas et al. 2012). In addition to these aspects, tagRNA-seq allows to visualize RNA processing sites and shows that the “processed RNA landscape” is an important part of the total transcriptome that has been less studied up to now. For example, the well-known ubiquitous sRNA RnpB, the ribozyme element of RNase P (Frank and Pace 1998), provides an illustration of the information accessible via tagRNA-seq. We previously mapped the rnpB TSS at location 1,121,939 in the E. faecalis V583 chromosome (Fouquier d'Hérouël et al. 2011), which is not detected by tagRNA-seq, most likely due to the higher amplification of the signal via the RACE-derivative method compared with the SOLiD procedure. The functional RnpB molecule, also termed M1, originates from a series of maturation processes conserved across the three domains of life that we may reasonably speculate to also operate in E. faecalis due to the high degree of structural and functional conservation of RnpB (Li et al. 1998; Mann et al. 2003; Griffiths-Jones et al. 2005, and references therein). TagRNA-seq data enables us to map locations 1,121,951/-53 with high tag counts corresponding to PSSs (Supplemental Tables SA, SE; the ppRNome browser). The RnpB upstream-most 5′ end predicted in the Rfam database allocates a position at 1,121,944 in the chromosome (Griffiths-Jones et al. 2005), a location spaced by 4 and 7 nt from the TSS and PSS we have mapped, respectively.

Further experiments will be necessary to shed light on the details of the processed transcriptome and its complex organization. Nevertheless, to our knowledge, this is the first study mapping PSSs at a global scale in bacteria.

Transcription start sites and processing sites in E. coli

Unlike in E. faecalis, transcription start sites in E. coli have been extensively studied and TSSs have been mapped in genome-wide manner in at least three separate contributions: Mendoza-Vargas et al. (2009) reported TSSs for ∼1000 (∼23%) of the ∼4500 ORFs in E. coli MG1655, Conway et al. (2014) mapped with high accuracy about 2100 promoters, and Thomason et al. (2015) detected over 14,000 TSSs using dRNA-seq.

In order to challenge the tagRNA-seq method and our analysis, we applied the same procedure for the E. coli transcriptome as we applied for E. faecalis with a significance threshold set to five reads. We were thus able to retrieve 397 TSSs in the U00096.3 reference genome (Supplemental Table SB). This lower number compared with E. faecalis can be explained by the smaller number of reads obtained from this sequencing experiment (see Supplemental Material, Section S2) and the 40% larger genome of E. coli, out of which ∼33% (1.55 Mbp) appear to be transcribed in our experiment (coverage higher than 2×). Out of those 397 TSSs, 314 (79%) were found within 2 bp of a TSS reported in at least one of the three contributions described above (113 in Mendoza-Vargas et al. [2009], 218 in Conway et al. [2014], and 285 in Thomason et al. [2015]).

In contrast, relatively few RNA processing sites (PSSs) have been mapped with single nucleotide accuracy in E. coli in standard growth conditions. Supplemental Section S5 and Table S4 provides 16 examples of PSSs reported in the literature and how they appear in our results of E. coli: (i) 11 PSSs are clearly recovered and fall in the area >60°, albeit three carry tag counts below the chosen threshold of five reads; (ii) five PSSs reported elsewhere are found within the area called “undetermined” (Fig. 2). These examples support that the tagRNA-seq method enable us to map PSSs within the bacterial RNA landscape.

Nonannotated genes, small RNAs, and particular transcriptional organization

Up to now, transcriptomic studies in E. faecalis have used microarrays designed to examine expression of annotated ORFs (Aakra et al. 2005; Solheim et al. 2007; Makhzami et al. 2008; Vebø et al. 2009, 2010; Abrantes et al. 2011; Mehmeti et al. 2011; Vesić and Kristich 2013), or custom-made tiling arrays containing a limited number of intergenic regions (IGRs) to search for sRNAs (Shioya et al. 2011). Although informative, these approaches provide partial information on the bacterial transcriptome compared with RNA-seq methods (Sittka et al. 2008; Rasmussen et al. 2009; Toledo-Arana et al. 2009; Sharma et al. 2010; Chao et al. 2012; Nicolas et al. 2012). We took advantage of our 5′ RNA end mapping for a detailed transcriptional analysis, looking for previously nonannotated genes in the genome of E. faecalis V583. Among other transcripts, sRNAs were primarily identified as stand-alone signals, whose length is up to 500 nt, located in “empty” regions (i.e., nonannotated regions), or transcripts antisense to annotated ORFs. In addition to the sRNAs previously identified (Fouquier d'Hérouël et al. 2011; Shioya et al. 2011), we unveiled a total of 85 novel sRNAs (Fig. 3; Supplemental Table SF). Continuing on our previous nomenclature (Fouquier d'Hérouël et al. 2011), these new sRNAs were named from “Ref47” to “Ref120” when present in the chromosome, and for sRNAs encoded by plasmids pTEF1 and pTEF2, from “RefA1” to “RefA9” and “RefB1” to “RefB12,” respectively. Among the 45 sRNAs previously reported in Fouquier d'Hérouël et al. (2011), 18 are confirmed by tagRNA-seq. For the remaining 27, we detect some low and inconclusive signals for about half of them, and do not detect anything for the rest. Four among the 11 unnamed sRNAs reported in Shioya et al. (2011) were confirmed and named Ref57 (IGR ef0408-9 in the chromosome), Ref78 (IGR ef1368-9), RefA9 (IGR efa0080-1 in pTEF1), and RefB12 (IGR efb0062-63 in pTEF2) (Supplemental Table SF). Additionally, comparison with Table S3 in the supplementary material of Shioya et al. (2011) confirmed six sRNAs in our list (Ref81, Ref115, RefA3, Ref93A,B, RefA8, and RefB11), as well as the antisense to RefA7 observed in our experiment but unnamed due to absence of tags, and the long 3′ UTR of ef3249.

FIGURE 3.

Global view of sRNAs and antisense organizations currently known in E. feacalis V583 (the chromosome and plasmids pTEF1 and pTEF2). All the sRNAs from Supplemental Table SF are represented, in “gray” on the forward strand and “orange” on the reverse strand. The inner plot visually describes with vertical blue lines the location and importance of antisense organizations detected (see the ppRNome browser for details). On the chromosome, the pathogenicity island (red) and other mobile genetic elements are annotated on the chromosome, i.e., efaC1/C2 (dark green), vancomycin resistance region (orange), and the six prophages (bright green) (Lepage et al. 2006; Matos et al. 2013).

sRNAs have important regulatory functions. Generally, regulatory RNAs not embedded in a transcriptional antisense organization (stand-alone) modulate the activity of proteins or affect translation (up or down) by pairing to mRNAs; a class of sRNAs also named “trans-acting sRNAs” (Repoila and Darfeuille 2009; Waters and Storz 2009). Although not functionally characterized so far, many sRNAs found in E. faecalis are likely trans-acting regulators, e.g., Ref50, Ref52, Ref72, Ref78, Ref79, Ref95, Ref102, RefA1, RefA4 (Supplemental Table SF). Some sRNAs have been shown to carry a dual function since they can exert their regulatory role via different mechanisms or can also encode peptides (Wadler and Vanderpool 2007; Loh et al. 2009; Livny and Waldor 2010; Jørgensen et al. 2012; Sayed et al. 2012). Some of the sRNAs predicted here may also encode for peptides as previously predicted for other sRNAs in E. faecalis (Fouquier d'Hérouël et al. 2011). Another category of sRNAs are antisense RNAs (asRNAs), transcribed from the complementary DNA strand of genes, and thereby forming transcriptional antisense organizations. As regulatory consequence, the expression of an asRNA can potentially impact the transcription initiation efficacy of the opposite gene (promoter interference), may provoke premature arrest of transcription elongation, and/or possibly modulate the translation and the stability of the cognate RNA (Georg and Hess 2011; Brantl 2012; Sesto et al. 2013). Many of the novel Ref sRNAs appear to form antisense transcriptional organizations, e.g., Ref89 and RefB4 are antisense to sRNAs Ref90, RefB5, respectively; Ref94, Ref114, and Ref115 are antisense to transcripts bearing ORFs ef2025, ef3087, and ef3088, respectively (Supplemental Table SF).

Also, long 3′ UTRs have been reported in several bacterial species and in a few cases their involvement in RNA-mediated regulation has been demonstrated (Sittka et al. 2008; Chao et al. 2012). In Supplemental Section S6, we present several cases found in the genome of V583.

In addition, antisense transcriptional organization also results from overlapping mRNAs and may involve coding sequences as well as 5′ or 3′ UTRs (e.g., Rasmussen et al. 2009; Toledo-Arana et al. 2009; Sharma et al. 2010; Nicolas et al. 2012; Wurtzel et al. 2012). For instance, the 5′ UTR ef0282 (fabI) overlaps the 5′ UTR ef0283 (fab-F1); ORFs ef0479, and ef0480 are embedded in a long opposite transcript originating 3000 bp upstream; and the transcript that contains ef0522–ef0523 in an operon is antisense to a transcript carrying ef0524. Similar examples are observed on plasmids pTEF1 and pTEF2 (the ppRNome browser and Supplemental Table SF). However, one of the most striking antisense organizations was found in the region spanning from ef2298 to ef2324 (Fig. 4). It is clearly visible in all of the E. faecalis transcriptomes, regardless of the growth condition, tagging, or sequencing protocol. It encompasses ∼22 kbp on the chromosome and involves two transcribed regions of 16 and 17 kbp long that overlap by 11.5 kbp. In the positive direction, the transcribed RNA originates 265 bp upstream of ef2304, the unique predicted ORF contained within this 16 kbp long RNA, and would code for a putative transcriptional regulator (Paulsen et al. 2003). In addition, this RNA is antisense to ef2312 and ef2314 that code for the DNA topoisomerase III (TopB-2) and a putative bacteriocin, respectively. On the negative strand, the second RNA originates 225 bp upstream of ef2308, and carries ef2298 and ef2299 (Fig. 4). These later ORFs encode the two-component regulatory system VanRB/SB, a vital element for E. faecalis V583 to resist to vancomycin, a major clinical antibiotic against Gram-positive infections (Huycke et al. 1998; Arias and Murray 2012). Experimental validation will be required, but it is tempting to speculate that this antisense regulation may control vancomycin resistance in E. faecalis V583.

FIGURE 4.

Long antisense organization in the chromosome of E. faecalis V583. The transcriptional antisense organization encompasses 22 kbp. Note that only a single ORF, ef2304, is predicted in the transcript originated from the positive DNA strand at coordinate 2221569. ORFs, ef2312 and ef2314 encode the DNA topoisomerase III and a putative bacteriocin; they are not transcribed in the growth conditions used but are “covered” by the antisense RNA (17 kb). vanRB and vanSB, encoding the two-component regulatory system of the vancomycin resistance locus, are contained at the end of the transcript originated at coordinate 2233229 from the minus DNA strand. The two antisense RNAs overlap by 11.6 kb. Coordinates mapped for TSSs of each corresponding transcript are noted at the corresponding location. The color boxes denoted by “Rt,” “St,” “KTH,” “KTHr,” and “IlluminaSt” applied to RNA levels of corresponding transcriptomes. RNA levels shown are normalized.

Differential gene expression in static and respiratory growth

We performed a standard differential expression analysis between the R and S conditions for E. faecalis (see Materials and Methods; Supplemental Table SD). Annotation and gene ontology information show that most of the differentially expressed genes in Supplemental Table SD are involved in the central metabolism associated with the consumption of sugars and amino acids as sources of carbon and nitrogen, mainly centered on glycolysis, pyruvate metabolism, and citrate cycle (TCA) (http://www.genome.jp/kegg/pathway.html and Fig. 5). For instance, operons celAB, ef1017–ef1020, and ef2959-ef2961 have a higher expression in S compared with R conditions, and are predicted to code for phosphor-sugar transfer systems (PTSs) of cellobiose and ribose uptake, and degrading enzymes, respectively. These transporter systems allow the entry and the phosphorylation of specific sugars that are catabolized via the glycolytic pathway. In R, compared with S conditions, an increased expression is detected for operons ef0097–ef0100 that code for the subunits of the serine-dehydratase (Sdh) that converts serine into pyruvate (or the reverse reaction), the end product of glycolysis, and ef1657–ef1663 coding for enzymes that feed or make part of the TCA cycle and glycolysis, and participate to amino acid metabolism (valine and serine) from pyruvate. We also observed an increased expression in R for operon glpFOK (ef1927–ef1929) coding for enzymes enabling glycerol catabolism via glycolysis in aerobic conditions (Bizzini et al. 2010).

FIGURE 5.

Metabolic pathways deduced from functions assigned to RNAs differentially expressed in S and R growth conditions. Transcriptomes St and Rt both show significant expression of genes involved in central metabolism. Differential expression analysis indicates that groups of enzymes specifically involved in sugar uptake systems (PTSs) and nucleotide bases metabolism pathways are induced in static growth conditions (S). Similarly, in respiration (R), the formation of acetyl-CoA from pyruvate, formate, and ethanol is enhanced, as well as the citrate cycle (TCA) and amino acids metabolism. Genes with up-regulated expression in each condition are represented by thick arrows in the corresponding plot. The names on the arrows correspond to gene names as given in Supplemental Table SD, which also contains their assigned functions from the KEGG database.

The results of this standard differential expression analysis are in line with previous reports on metabolic and transcriptomic studies performed on E. faecalis and closely related species such as Lactococcus lactis (Garrigues et al. 1997; Redon et al. 2005; Pedersen et al. 2008). Under S conditions (anaerobic growth) and, in the growth medium used in this study, E. faecalis adopts fermentation, a physiological state characterized by a high rate of sugar consumption in comparison to respiration, and employs the glycolysis pathway for glucose catabolism, ending up in pyruvate formation (Garrigues et al. 1997; Pfeiffer et al. 2001). In E. faecalis, the respiratory chain is complete but not functional as genes encoding enzymes required for heme synthesis are absent. However, if oxygen and heme are supplied as in R growth conditions, E. faecalis switches from homolactic fermentation to respiration and under respiration, pyruvate conversion is not exclusively catalyzed by Ldh-1, but also by Als (Ef1213), Pdh (Ef1353/Ef1354), Pfl (Ef1612/Ef1613), and AdhE (Ef0900), ending in the production of various products (acetate, ethanol, formate), and regenerate NAD+ from NADH (Winstedt et al. 2000; Yamamoto et al. 2006; Pedersen et al. 2008).

In summary, with the caveats of high expected false negative rate when comparing few samples, it can be checked that all the reported variations in Figure 5 and Supplemental Table SD are consistent with previously reported features of E. faecalis biology with and without respiration.

CONCLUSION

In this work, we have introduced a new method to distinguish primary and processed RNA ends and achieved the first RNA-seq transcriptome of E. faecalis. The discovery of numerous sRNAs and antisense organizations in E. faecalis transcriptomes highlights, as in many other species, the importance of RNA-dependent regulatory processes. The association of the RNA-seq method with the differential labeling of 5′ RNA ends, enabled us to provide the two “faces” of a bacterial RNA landscape. We mapped 559 TSSs and predicted promoter motifs at the genome-wide scale in a species where <50 were previously known, and 352 major PSSs, providing a first snapshot of a bacterial processed RNA landscape. As TSS- and PSS-tags hallmark transcription initiation and processing, the next step in the exploitation of our technology allowing identification of 5′ ends of primary and processed transcripts will be to perform quantitative studies in order to pinpoint the contribution of RNA synthesis and RNA stability in gene expression reprogramming accompanying physiological adaptation. This study constitutes a significant advance in the understanding of the organization and the expression of the genetic information of the human pathogen E. faecalis, and a key improvement of the functional analysis of bacterial transcriptomes.

MATERIALS AND METHODS

Bacterial growth and RNA preparation

E. faecalis V583 (VE14002 in our laboratory collection) was grown in brain–heart infusion (BHI) medium at 37°C in static (S) or respiratory (R) conditions. Briefly, one fresh bacterial colony was used to inoculate 5 mL BHI and grown at 37°C in S conditions for 20 h. A 1/500 dilution was then made in prewarmed BHI medium to obtain the test culture. In S growth conditions, the culture was not shaken. In R growth, hemin (SIGMA) was added to the test culture at 10 μM final concentration and culture was shaken at 200 rpm. Bacterial growth was followed by measuring the optical density at 600 nm (OD600). Total RNAs were prepared from bacterial cultures grown to an optical density (OD600) ranging between 0.7 and 0.85, as previously described (Fouquier d'Hérouël et al. 2011). In the course of this work, we discovered that our laboratory strain did not contain the plasmid pTEF3 (Paulsen et al. 2003); although we used the appellation of “V583” throughout the text, data presented are those obtained for our strain VE14002.

E. coli strain MG1655 was grown in LB medium at 37°C under agitation (200 rpm) until an OD600 of 0.5. Bacteria were pelleted and total RNA prepared as previously described (Fouquier d'Hérouël et al. 2011).

RNA tagging and sequencing

In bacteria, transcriptional start sites (TSSs, or “+1”) are characterized by the presence of a 5′ triphosphate group. In contrast, 5′ RNA ends created by endonucleolytic cleavages (PSSs) are 5′ monophosphate. We exploit this chemical difference by labeling differentially mono- and triphosphate 5′ RNA ends with two short RNA oligonucleotides, the “tags” (Fouquier d'Hérouël et al. 2011).

5′ RNA ends were differentially labeled with two short and different RNA oligonucleotides (tags) (Fouquier d'Hérouël et al. 2011; Supplemental Section S1). Briefly, primary transcripts contain 5′ ends with a triphosphate group which is brought by the first nucleotide triphosphate used by the RNA polymerase to initiate RNA synthesis at TSSs. In contrast, RNA processing events generate, at cleavage sites (PSSs), 5′ ends with monophosphate groups. RNAs with PSS and hydrolyzed 5′ triphosphate RNA ends were tagged by a first ligation step with the PSS-RNA adaptor (PSS-tag). Subsequently, RNAs were treated with the tobacco alkaline phosphatase (TAP) to transform triphosphate groups into monophosphate groups and were then tagged by a second ligation step with a TSS-RNA adaptor (TSS-tag). TSS- and PSS-tag sequences were adapted to RNA-seq in such a way that they cannot be mistaken with any regions of the V583 reference genome (Supplemental Section S1).

In order to account for variations in total number of reads and to be able to compare experiments, RNA levels are reported normalized to the total number of reads mapped, as commonly done in RNA-seq (Robinson and Oshlack 2010). Additionally, the ligation procedure introduces a new variability in the experiment that is corrected for by normalizing the number of tagged reads mapped at a given position to the total number of tagged reads mapped for the entire V583 genome (Supplemental Table SA).

Two RNA libraries were obtained from total RNAs prepared from E. faecalis grown in S and R conditions. They were tagged with PSS- and TSS-tags according to our 5′ RNA end discriminative method, treated according to SOLiD manufacturer's protocols for sequencing (Applied Biosystems, Life Technologies Corporation), and sequenced on a SOLiD 5500 platform (MetaGenoPolis, INRA). The corresponding transcriptomes were named “St” and “Rt,” respectively (Supplemental Table SA).

Additionally, as control, three other RNA libraries prepared from two independent growths in S conditions were sequenced. In one experiment, the bacterial culture was grown at Karolinska Institute, Sweden, as previously described in Fouquier d'Hérouël et al. (2011), and sequenced on a SOLiD v3 platform (Viiki) without applying our tagging procedure. Two libraries, denoted as “KTHr” and “KTH” respectively, were prepared from this experiment. For one of them, ribosomal RNAs (rRNAs) were removed using Ambion MICROBExpress Bacterial mRNA Enrichment Kit; in the other, rRNAs were retained. In the second S growth culture, bacteria were grown at INRA and RNAs where prepared as described here. Sequencing was done on a Hi-seq platform (IMAGiF, CNRS) following the Illumina Trueseq protocol, resulting in the “IlluminaSt” transcriptome. Due to technicalities in the ligation, the tagging in this sample was ineffective leading to a sample that here has been used only as standard RNA-seq.

A single RNA sample was prepared from E. coli total RNA, tagged using the same RNA adaptors (TSS- and PSS-tags) and sequenced on the SOLiD Wildfire platform (MetaGenomPolis, INRA). The resulting transcriptome was named “Coli.”

Alignment and coverage

Reads were aligned to the E. faecalis v583 and E. coli K12 substrain MG1655 reference genomes (respectively, GenBank Accession IDs [GenBank:AE016830.1] (chromosome), [GenBank:AE016831.1] (pTEF2), and [GenBank:AE016833.1] (pTEF1); [GenBank:U00096.3]) using Bowtie 1.0.0 (Langmead et al. 2009) with default options, but allowing for multiple matches (-a best command line option). The coverage is calculated by counting the number of reads mapped at each position on the genome for each strand. In case of multiple matches, the number of matches correspondingly divides the contribution to the read count. In cases mentioned explicitly in the text where repeated regions are excluded from the analysis; those multiply matched reads are ignored in the count. Similarly, when rRNAs and tRNAs are excluded, we impose a zero coverage over the corresponding regions. In order to reduce the effect of fragment bias (reads that are not uniformly distributed within the transcripts they represent) (Roberts et al. 2011), we define a quantity called “coverage density” similar to the coverage, except that reads mapped so that they start at the same genomic position are counted only once. The resulting signal is thus less sensitive to the specific amplification of the different fragments at the cost of losing dynamic range. The coverage density signal has the useful and exploitable feature that edges of an expressed region are always staircase-shaped. We use the coverage density signal as a means to predict transcript edges from our RNA-seq data. Supplemental Section S2 compiles the raw sequencing output for the various transcriptomes performed.

Gene expression level

We calculated gene expression levels of annotated genes of the E. faecalis genome and performed differential expression analysis between R and S growth condition using Cuffdiff from the Cufflinks suite v2.1.1 (Trapnell et al. 2010). Cuffdiff was run on the Bowtie output files with the command line option -u --library-type fr-secondstrand using the genome of E. faecalis V583 and its annotation. Regions corresponding to rRNA in the annotation were masked using the -M option.

Comparing Rt and St, we obtained a list of 31 ORFs that displayed significant variation of RNA levels between the two tagged RNA-seq (Supplemental Table SD). We assessed the false positive rate by comparing to the results including in the analysis data from the “IlluminaSt” transcriptome in the S growth condition. This improves the prediction quality of some called genes (decreases P-values for 16 ORFs), adds 16 more called genes to the list (moves their P-values below threshold), and still calls 23 ORFs in the list. For the nine ORFs that are removed from the list, the P-value, on the other hand, increases above threshold. Results from the analysis are available in Supplemental Table SD and discussed in the Results and Discussion section.

Predictions of transcription start sites

Starting from the coverage density signal, we developed an iterative algorithm to detect transcribed regions, filtering out signals of low quality originating from sequencing errors or misalignments. The algorithm is inspired by the edge thinning operation in image processing (Davies and Plummer 1981). All regions where the signal is greater than a given but arbitrary confidence threshold are marked as “strong” signal. The signal in the immediate vicinity of this strong signal is recursively annexed to the strong signal region. All signals not marked as strong are discarded (Supplemental Section S3). The algorithm discriminates low signals within transcribed regions and eliminates those likely caused by noise. The orientation of the aligned reads and the edges of signals enable us to assign TSSs.

Detection of transcription starts and processing sites using 5′ tags

The addition of tags allows to readily map 5′ ends of RNA molecules and to discriminate primary transcripts (ligated to TSS-tags) from processed transcripts (ligated to PSS-tags). Prior to alignment, reads are sorted according to tag sequences or their absence and, when present, tag sequences are removed from reads, leaving only sequences from bacterial RNAs. Both operations are performed simultaneously with Flexbar—Flexible Barcode and adapter removal for sequencing platforms—v2.4 (Dodt et al. 2012) allowing for up to two mismatches in the 13 nt of the tags (command line parameters --barcode-trim-end LEFT --barcode-threshold 1.6 --barcode-unassigned --barcode-min-overlap 9 --min-read-length 35). After alignment, reads with tags are classified into TSS or PSS candidates according to the rules described in the Results and Discussion section “The calling of TSS and PSS from tagRNA-seq.”

For transcriptomes obtained from the Illumina Hi-seq and SOLiD Wildfire, standard removal of 3′ sequencing adapters is performed using Flexbar in an additional preprocessing step. This step is not needed for SOLiD v3 and 5500 where the insert size is typically much longer than read length (Innocenti and Aurell 2013).

As a 5′ RNA end can be tagged by both TSS-tag and PSS-tag sequences (see below), we considered a 5′ end to be present in the sequenced RNA population when at a given location, the sum of TSS-tags and PSS-tags have counts of at least 3.2× per million of tagged reads aligned in E. faecalis transcriptomes. Such a threshold corresponds to a total of 5 tags (TSS- + PSS-tags) detected at the concerned position for the St transcriptome and seven tags for Rt. For the “Coli” transcriptome, this threshold was kept at 5 tags, or 4.06× per million of tagged reads aligned.

On one hand, a careful examination of the transcriptomes reveals many instances of one or two isolated reads in isolated positions. It is reasonable to assume that many of these are noise and thus setting the limit above this level was chosen to eliminate them. On the other hand, as described in detail in Supplemental Section S7, there seems to be no natural threshold in the data and a lower threshold simply leads to more candidates. The threshold values given above are simply one reasonable choice.

It is well known that transcription initiation at a transcription start site is not always initiated with single nucleotide accuracy (Schlüter et al. 2010; Sharma et al. 2010; Cortes et al. 2013; Morton et al. 2014). To take this into account, when locations distant by 4 bp or less from each other have mapped reads with TSS-tags and at least one of them is classified as a TSS candidate, the multiple tag signals are grouped in a single region that encompasses all those locations (Supplemental Tables SA, SB). The location with the highest tag signal is taken to be the most probable location of the TSS, and the total tag signal for the region is taken as the sum of the signals at all locations in the group. Although the length of such a “TSS region” can reach 6 bp or more in rare cases, many TSSs are detected with single nucleotide resolution (Supplemental Figs S8a,b in Section S8). Furthermore, an analysis of the average signal around retrieved TSSs shows that most of the signal concentrates within a region of ±2 bp around the most probable location (Supplemental Fig. S8c in Section S8).

As much less is known about the accuracy of the different ribonucleases, PSSs are reported in the text as point location on the genome and neighboring nucleotides with tag signal classified as PSSs are counted as different PSS sites.

Motif detection

We performed unbiased de novo motif search using MEME v4.9.1 (Bailey et al. 2009) upstream of genomic positions with TSSs and PSSs. The search was limited to the 10 most significant motifs with a width between 4 and 8 bp (command line arguments -nmotifs 10 -minw 4 -maxw 8). Short DNA sequences were extracted from the reference genome between [20–40 bp] and [0–30 bp] upstream of locations of interest and classified according to their ratios of PSS to TSS tag signals (as described in the Results and Discussion section “The calling of TSS and PSS from tagRNA-seq”). Those sequences were used as input to MEME without any filtering. When the feature was a TSS region as defined in the previous section, the most probable location was used as reference position for the sequence. The list of input sequences and results of the analysis are available Supplemental Table SC.

Visualization of results

All coverage information and tag signals resulting from our experiments and analysis can be visualized in a user-friendly and interactive manner online at the address http://ebio.u-psud.fr/eBIO_BDD.php (select the entry called “The ppRNome browser”). The visualization uses the Genome Browser (GBrowse) (Supplemental Section S9; Stein et al. 2002). The data presented online are also available in a numerical format in Supplemental Tables SA and SD for the Rt and St transcriptomes.

DATA DEPOSITION

The data are available in NCBI SRA repository, reference number PRJNA272574.

SUPPLEMENTAL MATERIAL

Supplemental material is available for this article.