High-throughput detection of RNA processing in bacteria.

BACKGROUND: Understanding the RNA processing of an organism's transcriptome is an essential but challenging step in understanding its biology. Here we investigate with unprecedented detail the transcriptome of Pseudomonas aeruginosa PAO1, a medically important and innately multi-drug resistant bacterium. We systematically mapped RNA cleavage and dephosphorylation sites that result in 5'-monophosphate terminated RNA (pRNA) using monophosphate RNA-Seq (pRNA-Seq). Transcriptional start sites (TSS) were also mapped using differential RNA-Seq (dRNA-Seq) and both datasets were compared to conventional RNA-Seq performed in a variety of growth conditions.
RESULTS: The pRNA-Seq library revealed known tRNA, rRNA and transfer-messenger RNA (tmRNA) processing sites, together with previously uncharacterized RNA cleavage events that were found disproportionately near the 5' ends of transcripts associated with basic bacterial functions such as oxidative phosphorylation and purine metabolism. The majority (97%) of the processed mRNAs were cleaved at precise codon positions within defined sequence motifs indicative of distinct endonucleolytic activities. The most abundant of these motifs corresponded closely to an E. coli RNase E site previously established in vitro. Using the dRNA-Seq library, we performed an operon analysis and predicted 3159 potential TSS. A correlation analysis uncovered 105 antiparallel pairs of TSS that were separated by 18 bp from each other and were centered on single palindromic TAT(A/T)ATA motifs (likely - 10 promoter elements), suggesting that, consistent with previous in vitro experimentation, these sites can initiate transcription bi-directionally and may thus provide a novel form of transcriptional regulation. TSS and RNA-Seq analysis allowed us to confirm expression of small non-coding RNAs (ncRNAs), many of which are differentially expressed in swarming and biofilm formation conditions.
CONCLUSIONS: This study uses pRNA-Seq, a method that provides a genome-wide survey of RNA processing, to study the bacterium Pseudomonas aeruginosa and discover extensive transcript processing not previously appreciated. We have also gained novel insight into RNA maturation and turnover as well as a potential novel form of transcription regulation. NOTE: All sequence data has been submitted to the NCBI sequence read archive. Accession numbers are as follows: [NCBI sequence read archive: SRX156386, SRX157659, SRX157660, SRX157661, SRX157683 and SRX158075]. The sequence data is viewable using Jbrowse on www.pseudomonas.com .

Background

Pseudomonas aeruginosa is a medically important γ-proteobacterium that is noted for causing opportunistic infections in hospitalized patients and chronic lung infections in cystic fibrosis patients [1]. A substantial cause of human morbidity and mortality, P. aeruginosa has been broadly studied due to its metabolic diversity and its ability to undergo substantial lifestyle changes that include biofilm formation, swarming motility and quorum sensing, adaptive responses to antibiotics, and complex virulence adaptations [1]. While P. aeruginosa PAO1 is the type strain for this model organism, detailed knowledge of transcriptional start sites (TSS) is currently lacking for this isolate. The post transcriptional modifications of RNA transcripts are largely unknown in Pseudomonas and are generally poorly studied for any living organism. In addition to enhancing our understanding of the basic biology of P. aeruginosa, the detailed mapping of TSS and subsequent RNA processing of transcripts involved in virulence, antimicrobial resistance, and essential cellular functions, will aid in understanding the regulation of pathogenesis and drug resistance, and facilitate the identification of promising drug targets.

A systematic inventory of TSS and RNA processing sites is fundamental to understanding a broad range of cellular processes. Determining the set of post-transcriptional modifications found in the transcriptome of an organism is an important and yet very challenging objective that can be partially addressed by RNA sequence-based analysis. After transcription, a series of highly regulated secondary modifications occur that result in the maturation of an RNA transcript. These processing steps strongly influence the overall lifetime of the RNA molecule and are instrumental in the functionality of many RNAs [2, 3]. A primary transcript contains a terminal 5′ triphosphate [4] (Fig. 1). In bacteria, selective removal of the 5′ triphosphate by the conserved pyrophosphatase RppH (YgdP in P. aeruginosa) leaves a 5′ monophosphate [5]. This 5′ monophosphate destabilizes mRNAs by making them more susceptible to degradation. An important complex in this regard is the multi-subunit degradosome, which contains the 5′ phosphate sensitive exonuclease/endonuclease RNase E at its core [3]. The presence of a 5′ monophosphate significantly enhances RNase E’s endonuclease activity [6]. Similarly endonucleases that cleave RNA so as to leave a 5′ phosphate can either activate RNA for degradation via pathways that include the degradosome, or can result in the production of stable RNAs essential for cellular function. Transfer RNAs and the 4.5S RNA transfer-messenger RNA (tmRNA) that rescues ribosomes found on broken messages are all endolytically cleaved by RNase P, producing a stable 5′ phosphate terminus [3]. RNase E/G, which performs endolytic cleavage at A/U rich ssRNA regions, plays a central role in the 5′ maturation of the 5S RNA as well as the maturation of the 3′ ends of tRNAs [3]. Concurrently, RNase III serves to cleave dsRNA and has a primary role in the maturation of the 16S and 23S rRNAs [7]. The complex interplay between these and other endo- and exo- nucleases presumably acts on numerous other unstudied RNAs within a cell, helping to regulate the maturation and lifetime of expressed RNA. Studying such processing with high-throughput methodologies provides a significant window into understanding global aspects of transcriptional regulation.

Fig. 1

RNA transcription and processing. a Transcription of RNA is initiated from a promoter sequence (indicated in red) within the genome. Ribonucleoside triphosphate polymerisation results in a 5′ triphosphate at the 5′ end of the nascent transcript and a 3′ hydroxyl at its 3′ terminus. b mRNAs can be internally cleaved by endonucleases to yield two RNA fragments or can be degraded by exonucleases from either their 5′ or 3′ termini. The activity of exonucleases is often triggered by the selective dephosphorylation of a terminal triphosphate to a monophosphate by a pyrophosphatase. c This study focuses on all RNA processing events that result in either a 5′ triphosphate (dRNA-Seq) or 5′ monophosphate (pRNA-Seq) and that simultaneously contain a terminal 3′ hydroxyl

RNA transcription is one of the initial steps in a complex regulatory cascade that enables cells to synthesize and regulate the expression of cellular factors in response to environmental changes. The highly conserved process of sigma factor dependent transcriptional initiation is of central importance in all bacteria [8]. Sigma factors form part of the RNA polymerase holoenzyme during transcriptional initiation and determine which promoters are active in specific cellular states [9]. P. aeruginosa PAO1 has 24 putative sigma factors, of which 14 have yet to have their DNA binding sites identified [10]. These sigma factors and their associated regulators are responsible for the correct transcriptional response to changing environmental conditions including low oxygen [11], limited iron [12, 13], and overall nutrient levels [14]. Identifying TSS within the genome and defining upstream sequence motifs helps to identify sigma-dependent promoters. To date, only 83 TSS have been annotated in the PAO1 strain [15, 16]. Wurtzel et al. [17] performed a key expansion of annotated TSS in P. aeruginosa strain PA14 by employing dRNA-Seq to find 2117 putative TSS. Notably however, the PA14 strain differs from the PAO1 strain [18] in that it contains an additional ~ 200 genes, is known to be more virulent [19], and has an estimated 5977 open reading frames versus the 5688 in PAO1 [1]. There is therefore considerable benefit in systematically defining and exploring TSS in PAO1.

The dRNA-Seq method was pioneered [20] to comprehensively map TSS found in a prokaryotic genome that are expressed in a given condition. This was accomplished by sequencing, in an orientation specific manner, RNA transcripts with triphosphates (PPP) at their 5′ ends. RNA-Seq technology has quickly become the new standard in transcriptome analysis and the data derived from these experiments has allowed us to view transcriptomes in unparalleled detail (see, e.g. [17, 21–24]). Single base pair resolution maps of transcriptional products derived from high throughput sequence data allow for gene by gene quantification of expression levels and novel gene discovery. However, 5′ degradation occurs quickly in bacterial RNA samples and it is difficult to tell where TSS are located based on standard RNA-Seq data. dRNA-Seq allows for the identification of TSS by sequencing only those transcripts that contain triphosphates at their 5′ ends [20], but this method cannot be used to examine further processing of transcripts after their synthesis by RNA polymerase.

Here, we use monophosphate RNA-Seq (pRNA-Seq), to study RNA processing in P. aeruginosa strain PAO1. In addition, we have used the differential RNA-Seq (dRNA-Seq) methodology of Sharma et al. [20] to characterize TSS and have conducted RNA-Seq inventories under four different growth conditions, in addition to selected additional downstream experiments, to provide a more complete picture of RNA expression in this organism. In addition to locating 1741 5′ monophosphate cleavage sites, we have also identified the sequence motifs corresponding to these sites, and were able to propose specific nucleases that might be responsible for some of the observed cleavage events. In addition we identified 3159 probable TSS in PAO1, significantly expanding our understanding of TSS in P. aeruginosa. A fraction of these TSS were found to be arranged in antiparallel pairs, implying that transcriptional initiation at either site might be conditionally dependent on the other. Through further downstream experiments, we demonstrated that certain small non-coding RNAs (ncRNAs) show significant changes in expression during swarming and biofilm formation, suggesting important roles for these RNAs in determining these complex adaptations. Collectively, these studies have heightened our understanding of transcription and RNA processing in the γ-proteobacteria, revealing layers of RNA processing complexity that were previously unexplored.

Results

Our transcriptome analysis consisted of three main facets: First, we sequenced the genome of our isolate of PAO1-UW using the Illumina methodology to confirm that it is the ref-seq isolate [1]. Our PAO1-UW isolate was nearly identical to the published reference sequence NC_002516, differing by only 11 single nucleotide polymorphisms (SNPs), each of which had been previously shown to be variable among laboratory strains of PAO1-UW [25], as well as 10 small indels (See Additional file 1: Table S1A and B respectively for a summary). The reference sequence genome was then used as a basis to map all reads from libraries derived from the transcriptome of PAO1-UW. Second, we analyzed the transcriptome for RNAs terminated with a 5′ monophosphate using our new high throughput method, pRNA-Seq. Third, we used dRNA-Seq to identify TSS and supplemented this information with an RNA-Seq analysis of transcription under four distinct growth conditions selected to simulate both laboratory and a range of infectious conditions. For all libraries, reads were then mapped to either the plus or minus strands (except for the RNA-Seq data, where information on strand orientation was not available) of the PAO1-UW genome using criteria summarized in the methods. We integrated data from the 5′ termini of the strand orientation sensitive dRNA-Seq and pRNA-Seq libraries into a detailed map of the PAO1-UW genome that identified both cleavage sites and TSS. RNA-Seq data was used to calculate read depth (see Methods) that again typically correlated well with both TSS and sites of major RNA processing as determined by dRNA-Seq and pRNA-Seq respectively. Statistics on the size and composition of each library can be viewed in Table 1. For online Jbrowse access to our data sets see www.pseudomonas.com

Table 1

Library growth conditions and summary of the number and percentage of reads mapped to the Pseudomonas aeruginosa PAO1 reference sequence and including (+ RNA) or excluding (− RNA) rRNA, tRNA and tmRNA genes

Sample		including rRNA		excluding rRNA
	Total Reads	Reads Mapped	Percent Mapped	Reads Mapped	Percent Mapped
pRNA-SeqLB Medium, 37 °C, OD600 = 0.7 (A06027)	272,983,632	149,142,710	55%	51,067,158	19%
dRNA-SeqLB Medium, 37 °C, OD600 = 0.7 (110817_SN865)	131,130,793	108,900,693	83%	29,801,108	23%
RNA-SeqSynthetic Cystic Fibrosis Medium, 37 °C, OD600 = 0.7 (A03674)	54,000,716	51,910,678	96%	16,003,536	30%
RNA-SeqArtifical Sputum Medium, 37 °C, OD600 not determined (A06026)	63,854,706	61,850,374	97%	6,802,731	11%
RNA-SeqLB Medium, 37 °C, OD600 = 0.7 (PA0004)	75,437,354	72,729,310	96%	30,068,678	40%
RNA-SeqLB Medium, 34 °C, OD600 = 0.7 (PA0001)	20,643,394	17,559,382	85%	3,352,489	16%

In all libraries but 110817_SN865, which consists of single end reads, the mapped reads are in proper pairs with a maximum insert size < 1000. A moderate drop in read quality (55% aligned reads) was observed with the dRNA-Seq library, which we attributed to the extra RNA manipulation steps required for library construction. We employed the same rigorous alignment quality thresholds with the dRNA-Seq library as with all other libraries to ensure that only high quality reads were mapped

Terminal monophosphate RNA data analyzed by pRNA-Seq revealed both expected and novel transcriptome processing sites

RNA processing sites are defined here as genomic locations with ≥100 reads “first bp coverage” from the pRNA-Seq library (i.e. oriented and aligned sequences whose 5′ most nucleotide was found at least 100 times; such sequences can and often had non-homogenous 3′ ends – see Methods). Processing sites were found that corresponded to previously-characterized processing sites of RNAs involved in translation. The ssrA RNA (tmRNA) was matured at precisely the residue expected based on data from other γ-proteobacteria (Fig. 2a) [26], as were the 12 ribosomal RNAs (4 copies of each of 5S, 12S and 23S RNA) [7]. However, the ribosomal RNAs also contained multiple internal cleavage sites. This was unexpected, as ribosomes are believed to be relatively stable structures within the cell, with core rRNAs that are resistant to the action of nucleases. The degradation process of ribosomes has, however, not been intensively studied [7], and it is currently unclear exactly what fraction of ribosomes are degrading during normal growth conditions. We examined the forty-one tRNAs known to be expressed in the pRNA-Seq library growth conditions (from our LB 37 °C RNA-Seq data) for cleavage sites. Isoleucine, alanine, serine and leucine tRNAs were observed to be processed to monophosphates at their 5′ termini, as would be expected based on the activity of RNase P. Overall a total of 13 tRNAs were cleaved at some point within the transcript, and by lowering the cut off threshold for cleavage site prediction (see Methods) from 100 reads to 50 reads, 5 additional tRNAs were revealed to be subject to cleavage. Transfer RNAs within an operon containing multiple tRNAs (i.e. Fig. 2b) showed evidence of not only 5′ processing but also 3′ processing and occasionally internal tRNA cleavage; this has been previously described as a mode of tRNA regulation in other bacteria, as well as in eukaryotes [27, 28]. Typically tRNA processing signals were more prominent towards the 5′ ends of the operons as would be expected from our use of a random primer reverse transcription (RT) step to generate cDNA after the initial ligation of a 5′ adapter sequence. Given that our RNA isolation method preferentially excluded tRNAs, it is unsurprising that they were not always present in large enough quantities to pass the step size threshold filtering steps for the pRNA-Seq library. Together, these data served as important biological controls for the overall pRNA-Seq methodology.

Fig. 2

Identification of precise RNA processing events within the ssrA transcript (tmRNA) and tRNA operon transcripts. a A histogram showing precise processing of the ssrA gene (by RNase P), with the 5′ ends of processed pRNA-Seq paired end reads aligning exactly at a single genomic location. Transcription for this gene is initiated 60-nt upstream of the location of RNase P cleavage. b Conventional RNA-Seq suggests the possibility that complex RNA processing occurs within the tRNA operon shown. pRNA-Seq indicates that a series of precise RNA cleavages occur downstream from a single strong initial transcriptional start

In addition to finding the expected evidence of RNA processing in certain transcripts, we observed evidence that unique cleavage events occurred in a diverse set of transcripts. Overall we identified 1741 5′ monophosphate RNA processing sites in the PAO1 transcriptome that met our statistical criteria for significance (Fig. 3). Five hundred putative cleavage sites not overlapping with pyrophosphatase sites were identified in 240 protein-coding mRNAs, which were among the most heavily transcribed in both the dRNA-Seq and RNA-Seq LB 37 °C libraries. Of the 240 most abundantly transcribed protein coding genes in the LB 37 °C RNA-Seq library (Additional file 1: Table S2), transcripts from 111 of these genes, or 46%, were also present in the dRNA-Seq library. Forty three (nearly 18%) of these genes were ribosomal proteins, which was by far the most highly represented Kyoto Encyclopedia of Genes and Genomes (KEGG; [29]) category among the transcripts. Three ncRNAs (PA4406.1, rnpB and crcZ) also showed evidence of cleavage. In addition, we found clear evidence of cleavage in transcripts for a wide range of proteins associated with critical cellular functions other than protein translation. For example, cleavage in PA3648/opr86 was detected, particularly in the central region of the gene. PA3648 encodes the only essential integral outer membrane protein in P. aeruginosa PAO1 – a protein critical for outer membrane biogenesis [30].

Fig. 3

Circular plot showing distribution of mapped RNA-Seq reads and 5′ monophosphate cleavage sites throughout the Pseudomonas aeruginosa PAO1 genome. The outer green and blue tracks represent reverse-strand and forward-strand genes, respectively. The third track from the outside is a heat map showing log10 first base-pair coverage (100 bp window) of RNA-Seq reads where a transition from yellow to green to blue correlates with increased transcription. The histogram with a grey background shows log10 coverage of 5′ monophosphate sites on both the forward and reverse strands. This is followed by a track containing rRNA (green), tRNA (blue) and ribosomal proteins (purple). The innermost track shows G-C skew (1000 bp window)

RNA cleavage sites within annotated genes were often clustered in the 5′ untranslated region of transcripts and were correlated with the reading frame position

Fig. 3 presents a global overview of the locations of cleavage sites within the genome. The RNA cleavage site locations within protein-coding genes were analyzed relative to the start and stop sites of annotated protein coding genes. This revealed that 31.6% (158/500) of cleavage sites fell upstream of the translation initiation site, with the remaining cleavage sites being nearly uniformly distributed across genes (Fig. 4a). We sought to determine whether there was a correlation between ribosomal binding site (RBS) location and cleavage site location. Prodigal software [31] was used to predict RBSs upstream of PAO1 protein coding sequences. Of the genes possessing a recognizable RBS and for which we had both TSS and RNA cleavage data, 2% (11) produced transcripts that were cleaved either within or immediately downstream of (i.e. 3′ to) the RBS. The open reading frames (ORFs) of such cleaved RNAs would presumably lack a RBS and consequently be expected to be poorly translated. More than half of the cleavage sites within gene ORFs were associated with a particular codon position, with 267 out of 500 (53.4%) of the ORF cleavage sites being located immediately following the first base in a codon, 5’-N1^N2 N3 (^: site of cleavage, Fig. 4b). This is likely due to the notable G + C bias in the genome influencing the location of G and C nucleotides in the cleavage site motif. Most of the cleavage sites located within genes were disproportionately located within transcripts for genes associated with specific functional categories as defined by the KEGG [29]. These included the basic bacterial functions of oxidative phosphorylation (54 cleavage sites, 14 genes, corrected p-value = 0.0000049) and purine metabolism (37 cleavage sites, 13 genes, corrected p-value = 0.0021), indicating a potentially important role in the regulation of core metabolism by the post-transcriptional secondary cleavage of RNA transcripts. Protein coding genes that contained cleavage sites tended to code for essential cellular functions. These genes encoded proteins that acted the toluene degradation pathway (3 genes, corrected p-value 0.004.4), and RNA polymerization (2 genes, corrected p-value 0.015).

Fig. 4

Relative cleavage site position within genes and reading frame dependent cleavage bias. a The relative distance that cleavage sites fall within protein coding genes was plotted in a histogram (see materials and methods). Cleavage sites are shown in turquoise, 5′ monophosphate sites occurring at the location of a TSS (i.e. potential 5′ pyrophosphatase sites) are shown in yellow. The remaining sites are distributed relatively evenly throughout the ORF. b The reading frames of cleavage products were determined for both positive (navy blue) and negative (red) stranded transcripts. Cleavage at site 1 occurs 5′ to nucleotide 1 (↓N1 N2 N3), cleavage at site 2 occurs 5′ to nucleotide 2 (5’-N1↓N2 N3) and cleavage at site 3 occurs 5′ to nucleotide 3 (5’-N1N2↓N3)

RNA cleavage patterns correlated with RNA cleavage motifs

We next searched for RNA sequence motifs associated with cleavage sites within annotated genes to determine whether specific nucleases might be involved in their creation. Processing sites found at locations meeting our peak height cut off criteria, and not overlapping with pyrophosphatase sites were aligned and examined for potential patterns associated with nuclease digestion. These potential patterns were inspected in a 10-nt window upstream and downstream of the dominant cleavage site which had the highest read coverage in the window. This alignment was used to explore the hypothesis that distinct RNA digestion patterns might be correlated with specific RNA sequence motifs. Strikingly, a single global motif dominated the dataset (Fig. 5a), [(A,C,g,u)(A,C,g,u)(G,a,c)(A,g,u)↓(A,c,u)(C,u)(A,c,g)(C,a,g,u)(C,a,g)], wherein nucleotides present in 10–30% of the sequences are depicted in lower case letters, nucleotides present in 31–64% of the sequences are depicted in uppercase letters and nucleotides present in 65% or more of the sequences are depicted in bold upper case letters. The position of the predominant cleavage site is indicated by the downward arrow (see also the graphical views of motifs in Fig. 5 and Additional file 1: Figure S2).

Fig. 5

5’ Monophosphate processing patterns and their corresponding sequence motifs. Cleavage sites are derived from genome locations with 100 reads or more 1st bp coverage from our pRNA-Seq library. Peak shape refers to the number of mapped transcripts surrounding a cleavage site. Peak shapes were categorized using k-mean clustering. Motifs were calculated from peak shape clusters using MEME [68]. The graph at the top of each panel shows normalized peak height. The WebLogo [69] at the bottom of each panel shows the sequence motif associated with each peak shape. a The global motif, which is derived from the entire pRNA-Seq dataset. b the “Sharp” peak shape motif, shows strong similarity to the RNase E motif in E. coli [38]. c the “Tail L” motif and (d) the “Tail R” motif

K-means clustering decomposed the RNA digestion patterns into five distinct classes (plus a ‘noisy’ class) based on peak shape that correlated with specific RNA sequence motifs (See Additional file 1: Figure S1). By peak shape we mean the shape of the pattern produced by the mapped 1st bp coverage of transcripts surrounding a cleavage site. For more in-depth analysis of motifs surrounding cleavage sites, we chose to analyze cleavage sites that fall within ORFs and are > = 10 nt downstream of another peak. Of the 383 cleavage sites analyzed, 97% fell into one of the categories described below (with the others falling into the ‘noisy’ class). For a list of all RNA sequence cleavage motifs, see Additional file 1: Table S3. The predominant ‘Sharp’ RNA cleavage pattern (182/383 cleavage sites) consisted of a single RNA cleavage event with little or no cleavage at adjacent nucleotides, and therefore resembles a sharp peak when viewed graphically (Fig. 5b). The Sharp sequence motif was consistent with the overall global motif, containing the nucleotides [(A,C,g,u)(A,C,u)(G,a,c)(A,g,u)↓(A,c,g)(C,u)(A,c,g)(C,a,g,u)(C,a,g)], where the downward arrow indicates the cleavage site. RNAs that contained this cleavage motif were located in genes that disproportionately belonged to the KEGG categories “oxidative phosphorylation” (number of genes = 12, number of cleavage sites = 54, corrected p-value = 0.000026) and “purine metabolism” (number of genes = 12, number of cleavage sites = 37, corrected p-value = 0.0011) (Table 2).

Table 2

Association of cleavage site motif types with KEGG functional terms

Test by peak shape	KEGG functional terms	p value	Number of genes in set with this function	Number of Cleavage Sites by KEGG Term
Global	Oxidative phosphorylation	4.90E-06	14	54
	Purine metabolism	2.08E-03	13	37
	Aminoacyl-tRNA biosynthesis	8.94E-03	7	16
	RNA polymerase	1.72E-02	3	19
	Pyrimidine metabolism	4.67E-02	8	26
	Protein export	4.93E-02	5	37
Sharp	Oxidative phosphorylation	2.60E-05	12	54
	Purine metabolism	1.13E-03	12	37
	RNA polymerase	9.13E-03	3	19
	Pyrimidine metabolism	1.14E-02	8	26
	Aminoacyl-tRNA biosynthesis	1.89E-02	6	16
Tail L	RNA polymerase	7.11E-05	3	19
	Protein export	1.24E-02	3	37
	Purine metabolism	1.39E-02	5	37
Tail R	RNA polymerase	4.52E-02	3	54
	Protein export	1.51E-02	2	37
Twin R	Oxidative phosphorylation	9.70E-04	5	54
	Protein export	1.03E-02	3	37
	Oxidative phosphorylation	2.52E-02	4	54
	Citrate cycle (TCA cycle)	4.94E-02	3	39

Certain cleavage site motifs are disproportionately located within genes that are associated with certain KEGG functional terms. The following table lists the KEGG terms associated with each cleavage site motif type and the p-values associated with each KEGG term. P-values were calculated from hypergeometric tests for KEGG pathways and Holm’s test was used to correct for multiple testing

The remaining cleaved RNAs sorted into RNA digestion patterns that showed asymmetries in their cleavage patterns (Fig. 5c and d). When viewed graphically, the second most abundant peak shape (58/383 cleavage sites) had a shoulder immediately 5′ to the dominant peak (Fig. 5c) and was named “Tail L”. This motif (Fig. 5c) contained the nucleotides [(A,C)(A,c,g,u)(G,a,u)(A,g,u)↓(C,U,a)(C,u)(A,c,g)(A,C,g,u)(C,a,g)] and was consistent with either two adjacent cut sites or an initial cleavage event followed by the removal of an additional nucleotide in the 5′- > 3′ direction. The motif for this cluster shared properties with the predominant Sharp motif but was notably lacking in sequence conservation at the + 2 and + 3 position. The Tail L motif also frequently (> 30%) had a U at the + 1 position that was absent in the predominant Sharp motif. Transcripts containing this cleavage motif included those from genes belonging to the KEGG categories “RNA polymerase” (number of genes = 3, number of cleavage sites = 19, corrected p-value = 0.000071), “protein export” (number of genes = 3, number of cleavage sites = 37, corrected p-value = 0.012) and “purine metabolism” (number of genes = 5, number of cleavage sites = 37, corrected p-value = 0.014).

The third most abundant peak shape (54/383 cleavage sites) had a shoulder immediately 3′ of its main peak and was named “Tail R” (Fig. 5d). This motif contained the nucleotides [(A,C,u)(A,C,u)(G,a,c,u)(A,G)↓(A,c,g,u)(U,c)(C,a,g,u)(A,C,g)(C,a,g)] and was quite different from the Sharp and Tail L motifs (Fig. 5c). RNAs containing this motif included transcripts from genes belonging to the KEGG categories “RNA polymerase” (number of genes = 2, number of cleavage sites = 19, corrected p-value = 0.0096) and “protein export” (number of genes = 3, number of cleavage sites = 37, corrected p-value = 0.010) among others (Additional file 1: Table S3, Additional RNA cleavage patterns were also identified – see Additional file 1: Figure S2 for details).

A subset of the pRNA-Seq-determined cleavage sites were found to localize exactly with our determined TSS locations. In total 131 sites (92 in coding genes + 39 in rRNA) met this criterion and were thus likely to be due to dephosphorylation of transcripts whereby the triphosphate was removed from the 5′ end of the RNA molecule and a 5′ monophosphate remained (Fig. 1). One of the steps during preparation of the dRNA-Seq TSS libraries was the removal of 5′ monophosphate-containing RNA using terminator-5′-phosphate-dependent exonuclease. It was possible that some TSS might have been identified as false positives due to the incomplete digestion of 5′ monophosphate RNA during TSS library construction. However, if indeed removal of 5′ monophosphate RNAs was only partially complete, then we would have expected that the 5′ termini of tRNA and tmRNA, which are known to possess 5′ monophosphates, would be identified in our TSS (dRNA-Seq) library data. Since this was not observed, we can tentatively conclude that either the desphosphorylation sites identified were biologically significant, or that the secondary structures with the 5′ termini of these RNAs might prevent dephosphorylation. Genes with transcripts that possessed dephosphorylation sites were not significantly associated with any KEGG terms. Due to the low number of dephosphorylation sites identified, motif analysis of the downstream sequence was inconclusive, as was an attempt to determine potentially conserved RNA secondary structure at these sites.

TSS prediction from dRNA-Seq identified promoter regions and novel sets of potentially co-expressed genes

We predicted a total of 3159 TSS from which RNA was actively initiated at 37 °C in LB media. These sites were associated with 2030 genes (in many cases, more than one TSS was associated with a single gene) that represented 36% of the strain PAO1 genome. Strikingly, just 54% (1695) of these TSS lay outside of ORFs. The remaining 1467 TSS lay within ORFs implying a potential regulatory role for such transcripts. We compared our predicted TSS to a set of 51 previously described TSS from strain PAO1 grown under various conditions [15], and for these genes 44 (86%) of our dRNA-Seq TSS lay within ±3 nt of previously published TSS (see Additional file 1: Figure S3).

TSS correlations between plus and minus strands

A correlation analysis for TSS was performed to explore the hypothesis that adjacent TSS might be related to a particular biological process. This analysis revealed a novel and previously unsuspected correlation between back to back TSS pairs. These antiparallel TSS pairs, based on opposite strands, were separated by 18 bp (Fig. 6a) and often contained a palindromic “TAT(T/A)ATA” motifs that could be − 10 promoter elements. These elements are equidistant between the two antiparallel TSS (Fig. 6b). When a similar correlation analysis was performed for pairs of TSS on the same strand, no significant correlation between peaks was found (Fig. 6a). A total of 105 antiparallel transcriptional pairs were found, with 59 of these sites being found between genes that were also in an antiparallel arrangement. Interestingly, 14 of these sites would be predicted to express antisense RNA for an otherwise parallel arrangement of genes.

Fig. 6

Correlation between TSS occurring in the antiparallel orientation and a palindromic sequence motif. a TSS show a significant correlation with TSS in the opposite direction that are 18-bp upstream (black curve, n = -18). No statistically significant correlation was found between TSS having the same strand orientation (red curve). The correlation function defined as: C(n) = Sum[x(l)y(l + n)], l = 1..G-n, with G being the genome size, x(l) plus strand TSS counts, y(l) either minus strand TSS counts (black curve) or plus strand TSS counts (red curve). The maximum value of the correlation function was normalized to one for each curve. b A “TATNATA” motif occurs between antiparallel TSS. Of the 105 antiparallel TSS found that where 18-bp appart, 3 match the motif exactly, 27 have one mismatch, 25 have two mismatches and 13 have three mismatches

PAO1 transcriptional promoter map

The dRNA-Seq data analyzed here was used to develop a promoter map for strain PAO1. Promoters were predicted for all (1612) primary TSS (TSS in an intergenic region and on same strand as downstream gene), including 111 primary antisense TSS (TSS in an intergenic region but on the opposite strand of the closest gene). Promoters can be viewed on Jbrowse at www.pseudomonas.com. As it has been shown that many virulence factors share a common transcription factor [32], analyses were conducted to identify putative novel virulence factors based on the similarity of promoter motifs to those upstream of known virulence factors. Three potential novel virulence factors that share promoter motifs with the gene flhA were identified (See Supplementary Material, Additional file 1: Figure S4). In addition, we aimed to determine whether novel binding sites for the known sigma factor RpoN could be identified based on sequence motif similarity. Thirty two putative novel RpoN binding sites were identified within promoter regions of known genes. Four of these genes were predicted to be regulated by RpoN, but binding sites had not been described, while 25 were upstream of genes where RpoN involvement in transcription had not previously been hypothesized (See Supplementary Material).

Chromosomal gene position affects transcription

Recently, it was reported that for short-read Illumina sequencing of bacterial genomes, sequence reads near the chromosomal origin of replication are more frequent than sequence reads distal to the origin. This is thought to be due to the nature of circular chromosome replication, such that there is a higher copy number of genes/sequences near the origin where DNA replication is initiated. Such read frequencies can even aid in the identification of genome rearrangements [33]. We examined whether RNA-Seq sequencing reads would have a corresponding bias towards higher frequency around the origin vs. the terminus of replication. RNA-Seq sequence reads did indeed show a decrease in frequency in a region near the known terminus (Fig. 3) that had been previously shown to also have reduced transposon mutagenesis frequency [34]. The fold-change of read density was also calculated in 0.5 Mbp increments along the genome when compared to the region with the lowest read density (Additional file 1: Table S6). In all RNA-Seq libraries, the region with the lowest read density lay between 2 and 2.5 Mbp, which is the location of the terminus of replication as revealed by a G-C skew plot (Fig. 3). There was an array of rRNA genes located between 5 and 5.5 Mbp, which raised the fold-change of read density values to very high levels. In addition, the highly transcribed tmRNA gene was located at 1.1 Mbp, thus increasing fold-change values in this region. The regions with the highest read density were proximal to the origin of replication, having implications for the analysis of gene expression and illustrating the importance of gene location in impacting its expression.

Confirmation, and additional functional analysis, of small non-coding (nc) RNAs

Estimates of the number ncRNAs vary widely and the small size of ncRNA genes can make TSS mapping challenging. Gomez-Lozano et al. [35] utilized an automated analysis to identify 513 novel ncRNA transcripts and verified the expression of several previously-known ncRNAs. Another study by Wurtzel et al. [17] on strain PA14 found 165 ncRNAs in total. Here we adopted a conservative approach by manually analyzing and reporting only on those 26 novel ncRNA transcripts that were clearly expressed based on our selection criteria [H(dRNA)cut_off = 500] for our RNA-Seq datasets (Table 3); these transcripts were subsequently confirmed by RT-qPCR (NB. four additional putative ncRNAs identified from the RNA-Seq data could not be confirmed by RT-qPCR.). As controls, 5 of the 63 previously-identified ncRNAs [16, 36], (P1, rsmY, prrF1, prrF2 and prrH) were PCR confirmed here (Table 3). Thirty of the 31 ncRNAs verified through RT-qPCR were differentially expressed (Table 3, Additional file 1: Figure S5) under one or both conditions of biofilm growth or swarming motility, which represent complex adaptive lifestyles in Pseudomonas.

Table 3

Small RNA species detected reliably by RNA-Seq and confirmed by RT-qPCR. Differential expression in biofilms and during swarming motility

Genomic Coordinates	Name	Size bp	Identity Gomez-Lozano, et al.	Identity Wurtzel, et al.	Complementarity (potential binding sites within other transcripts)	Fold change in biofilms	Fold change in swarming motility
143,349–143,517	PA0123.1	169	pant15	Not identified	None	1	10
326,875–327,066	PA0290.1	192	pant37	PA14sr_012	pilW, plcH	−9	1
334,491–334,686	PA0296.1	196	Not identified	P1	rne	−2	1
354,527–354,742	PA0314.1	216	pant42	Not identified	None	−5	1
586,867–586,990	PA0527.11	124	rsmY	rsmY	None	4	11
720,091–720,345	PA0667.1	255	Not identified	Not identified	PA3505, PA2897, PA0690	−3	1
798,865–799,255	PA0730.1	391	pant80	PA14sr_122	None	−3	44
883,307–883,582	PA0805.1	276	pant89	PA14sr_119PA14sr_120	None	−5	−5
1,045,414–1,045,733	PA0958.1	320	pant103	PA14sr_112	None	−6	2
1,182,820–1,183,057	PA1091.1	238	pant119	Not identified	PA0588	3	5
1,254,432–1,254,698	PA1156.1	267	pant125	PA14sr_105	PA1123, phuR	1	1
2,761,459–2,761,704	PA2461.1	246	pant225,	PA14sr_076,	PA2460, PA2458	−8	−3
2,761,599–2,761,911	PA2461.2	313	pant226	PA14sr_077	PA2460, PA2458	−5	1
2,964,898–2,965,137	PA2461.3	240	pant233	Not identified	PA5134	−2	3
2,977,373–2,977,611	PA2633.1	239	pant235	PA14sr_067	PA3672, recJ, nrdG, PA3522, PA3949, PA5325	5	8
3,312,577–3,312,693	PA2952.1	117	Not identified	PA14sr_061	PA4629	−2	2
3,545,572–3,545,872	PA3159.1	301	pant292	no PA14 ortholog	None	−2	3
3,697,226–3,697,433	PA3299.1	208	Not identified	Not identified	PA0690, cyoB	2	1
3,930,282–3,930,639	PA3514.1	358	pant326	no PA14 ortholog	tagF1	10	1
4,012,653–4,012,735	PA3580.1	83	pant337	Not identified	None	−4	1
4,536,541–4,536,848	PA4055.1	308	pant373	Not identified	PA2728, mfd, chpA, PA3641	−3	−2
5,080,450–5,080,630	PA4539.1	181	pant415	Not identified	wzz	2	10
5,207,898–5,208,463	PA4639.1	566	pant428	PA14sr_139	PA5156, PA2502, PA4510, aruI, PA0475, PA0558, PA1025	5	n/a
5,224,568–5,224,795	PA4656.1	228	pant430	Not identified	PA2038, PA3517, PA2152, pslE, PA2472, PA2750	−4	3
5,283,960–5,284,110	PA4704.11	151	prrF1	prrF1	prrF2	12	217
5,284,172–5,284,319	PA4704.21	148	prrF2	prrF2	prrF1	20	141
5,283,960–5,284,319	prrH 1	360	prrH	prrH	prrF1, prrF2	13	163
5,309,047–5,309,325	PA4726.3	279	pant439	PA14sr_141	ispA, hepA, PA2018, PA3461	−5	1
5,718,503–5,718,753	PA5078.1	251	pant465	Not identified	PA0312, kds	−6	3
5,973,539–5,973,681	PA5304.1	143	pant487	Not identified	None	1	5
5,986,186–5,986,427	PA5316.2	242	pant488	PA14sr_154	PA1302, PA2933, gcp, PA0241, PA0364, pilJ, hsiC2, PA2325, PA3037, rnhB, pchF, recD, algP	−7	−2

1RsmY and PrrF1–2, included as controls, were previously known to be up-regulated during the biofilm mode of growth [24]. The two small RNAs named PrrF1 and PrrF2 are adjacent and very similar in sequence and have been reported to also form a single transcript named PrrH

2n/a indicates that no consistent amplification or expression of that transcript was detectable under the given adaptive condition in three biological replicates. Expression was observed under mid logarithmic growth conditions

Four other putative sRNAs (68,836–69,271, 707,395–707,685, 830,970–831,031, and 99,801–100,048) could not be confirmed by RT-qPCR

Discussion

Here we introduce a technical advance in sequencing that characterizes post-transcriptional processing events that we call pRNA-Seq. In addition, we performed dRNA-Seq to investigate transcriptional start sites. This marriage of existing methodology with novel techniques provided new insights into the prokaryotic transcriptome and its regulation. We have provided a comprehensive map of cleavage and TSS sites for strain PAO1 grown under standard laboratory conditions (37 °C in LB media in the logarithmic phase of growth), characterized the associated promoters, identified RNA cleavage motifs, and further characterized ncRNAs.

Frequent RNA cleavage events as likely intrinsic check-points for RNA degradation

The continual turnover of RNA transcripts in a bacterial cell is a highly dynamic process that has been fine tuned by evolution. Our analysis of RNA post-transcriptional processing revealed three important steps involved in the eventual destruction of RNA species. Each process has a significant ramification for understanding the regulation of RNA turnover in the γ-proteobacteria.

Of the 240 most abundantly transcribed protein-coding genes, the transcripts from 111 (46%) were found as cleaved products in the dRNA-Seq library (Additional file 1: Table S2, bolded rows indicate transcripts also found in the pRNA-Seq library). Therefore, there appears to be a partial correspondence between the level of expression and the level of transcript cleavage. The observation that only about half of the most abundantly transcribed genes were processed, and that this occurred at specific sites, is consistent with the hypothesis that transcript processing reflects a regulatory process that serves to control RNA decay rate for only some of these abundant RNAs [37]. Intriguingly, nearly 18% of the most abundantly transcribed genes were ribosomal proteins, many of which were apparently processed, supporting the concept that the production of the protein synthesis machinery in actively growing cells is at least partly modulated by 5′ monophosphate generating cleavage events.

Our data indicate that cleavage sites in transcripts from protein-coding genes often occur in the 5′ untranslated regions of transcripts. We propose that this might be a post-transcriptional regulatory mechanism leading to modulation of translation. In addition, 1.5% of genes (19/1238) that possessed recognizable RBS produced transcripts that were cleaved such that their RBS were removed, even though the associated gene remained intact, implying a further level of translational suppression.

Remarkably, nearly all of the 500 cleavage sites from protein-coding transcripts in our data set corresponded quite closely to a sequence motif that, based on studies in E. coli, would be predicted to be cleaved by RNase E: [(G,A)(C,A)N(G)(G,U,A) ↓ (A,U)(C,U)N(C,A)(C,A)] [38]. Our data provides evidence in vivo for the widespread distribution and functionality of this cleavage motif, suggesting that RNase E plays an instrumental role in regulating the processing of an unprecedented number of cellular RNAs in Pseudomonas and, by extrapolation, the eubacteria. P. aeruginosa RNase E has 64% amino acid sequence identity to that of E. coli and the conservation of RNase E cleavage patterns across the γ-proteobacteria appears likely given the essential function of this enzyme [3, 7, 39].

P. aeruginosa has 66.6% G + C in its genome, therefore the third codon positions in this organism tend to be occupied by guanine and cytosine. Our analysis indicated that cleavage sites tended to be located between the first and second codon positions in annotated protein coding genes. Thus most of the cleavage sequence motifs showed a preponderance of G and C residues in their third codon positions (the − 2, + 2, + 5, etc. positions within the cleavage sequence motifs) (Fig. 5b and Additional file 1: Figure S2A). For example, the G residue at position − 2, which is known to confer rapid RNase E cleavage in E. coli [38], was found to be conserved in 65% of both the Sharp and Tail L motifs. The notable exception to this trend was the Tail R motif, where a U residue was conserved in more than 65% of the sequences at the + 2 location downstream of the cleavage site, at a 3rd codon position (see Fig. 5d). As this motif did not match any previously described motifs for known nucleases, the nuclease performing this cleavage, while possibly unique, is currently unknown. Sequence motifs such as Tail L, that showed RNA cleavage patterns 5′ of the cleavage site identified for the predominant Sharp cleavage site (Fig. 4c), likely resulted from a two-step cleavage process where a primary cleavage event would serve to recruit subsequent nuclease complexes that would ultimately degrade RNA in the 5′- > 3′ direction (consistent with a degradosome type of activity) [3].

A total of 119 peaks corresponded exactly to a TSS found in our dRNA-Seq library. As mentioned above, we believe that these are mRNAs that have been matured by a pyrophosphatase and are consequently bona fide entities within P. aeruginosa cell. tRNAs and tmRNA (both require 5′ monophosphates to be biologically active) were present in the dRNA-Seq library, but not in the TSS library. It is possible that some TSS were identified incorrectly due to incomplete enzyme digestion during library preparation. However, since we do not see any tRNAs or tmRNA in the TSS libraries, we can conclude that enzyme digestion (and elimination of transcripts possessing 5′ monophosphates) from the TSS library was successful. Therefore, we can also conclude that the overlapping sites from the dRNA and pRNA-Seq libraries exist within the cell, and that such transcripts are processed by pyrophosphatases after transcription. Triphosphates at the 5′ end of E. coli transcripts can be removed by the pyrophosphatase RppH; thus it seems possible that the homolog YgdP in P. aeruginosa (67% identity to E. coli RppH), may be responsible for this activity in Pseudomonas.

RNA cleavage resulting in RNA with a 5′ monophosphate was associated with a surprisingly diverse class of genes that had not been previously thought to be subject to such regulation and was found at notable levels in a surprisingly high number of protein-coding genes. Genes of note included PA3648/opr86 (an Omp85 homolog), which encodes the only known essential integral outer membrane protein in P. aeruginosa, involved in outer membrane biogenesis. This gene and many other coding genes were not known to be subject to RNA-cleavage based regulation. Notably, the impacted protein-encoding genes tended to code for essential cellular functions. Our analysis represents a starting point for more in-depth characterization of the important role of RNA cleavage in transcriptional regulation and overall cell stability.

A comprehensive TSS profile highlighted correlated antiparallel transcription and alternative promoters

Our observation of 105 back-to-back TSS spaced by 18 bp and often containing a palindromic − 10 A/T motif at the precise center of the back-to-back TSS has at least two mechanistic explanations. First, an individual polymerase holoenzyme complex could bind to either motif and form an open-form transcriptional bubble in the palindromic region through sigma factor specific interactions. Such complexes might transiently prevent transcription in the opposite direction and/or provide a competitive mechanism for transcriptional initiation. This first model is supported by numerous lines of biochemical and structural evidence (see, for example, [40-47]). Second, we propose an additional model where transcription is potentially initiated by an RNA polymerase dimer. In this model the spacing of the polymerase active sites on the dimer would be responsible for the 18 bp spacing observed in our data, while the palindromic − 10 A/T motif would facilitate the opening of the DNA duplex so as to allow either unit of the dimer to compete for transcriptional initiation (Additional file 1: Figure S6). This model would therefore predict approximately equal transcriptional initiation in either direction for a fully palindromic site and biased transcriptional initiation for an asymmetric site. This second model is consistent with the finding that RNA polymerase dimers have previously been observed during the purification of bacterial RNA polymerase. Furthermore, the authors of [48] also proposed a model of and RNA polymerase dimer binding to DNA based on the purification of E. coli RNA polymerase in association with DNA probes.

Previously, dRNA-Seq studies have been shown to be an accurate method for determining prokaryotic TSS in Helicobacter pylori, Anabaena sp. PCC7120, Trichodesmium erythraeum IMS101 and Burkholderia cenocepacia J2315 [20, 49–51]. In this study, we have effectively expanded the dataset of probable TSS in PAO1 by a factor of thirty. Our variation from published data (78% of our TSS lie within ±2 nt of the 51 previously described strain PAO1 TSS [15] (Additional file 1: Figure S3)) is slightly larger than that observed by the H. pylori group who pioneered dRNA-Seq methodology [20] (87% of TSS within ±2 nt of published TSS). This might reflect the possibility that multiple promoters are used to transcribe the same gene, experimental errors, or variations in growth conditions in the previously published studies cf. our study (e.g., different OD600 at harvest) resulting in the use of different TSS. P. aeruginosa, with its notably larger genome and number of transcriptional regulators cf. H. pylori, might have a more complex transcriptome. Overall, our predicted TSS showed a mean deviation of 5.0 bases from previously published TSS. Most of this variation was due to the presence of 6 outliers that had differences from published TSS ranging from 8 to 66 nt. Interestingly 3 of the 6 outliers were involved in transcriptional regulation, and may have multiple promoters enabling different transcriptional hierarchies depending on the growth conditions.

Wurtzel et al. [17] greatly expanded the catalogue of annotated TSS in P. aeruginosa strain PA14 by employing a 5′ triphosphate transcript mapping strategy. Strains PA14 and PAO1 are highly similar organisms, with genomes differing by only roughly 200 genes [52]. In both the PA14 study and the present study, the bacteria were grown in identical media, and very similar TSS mapping methods were used. The number of TSS identified in PAO1 in the current study was 3159, while Wurtzel et al. [17] identified 2117, a difference of more than 1000. The difference likely lies in the 25–100 fold number of total reads mapped to the non-rRNA regions in the genome in these studies (29,801,000 reads in our 5′ triphosphate library vs. 218,000 reads in Wurtzel et al.’s 37 degree 5′ triphosphate library, and 1,262,000 reads in their 28 degree 5′ triphosphate library [17]). This greater sequencing depth enabled more stringent cutoffs to be employed here than in the PA14 study (a threshold of 500 reads mapping to a single genomic location, cf. a threshold of 5 reads in PA14) [17].

In many bacterial RNA-Seq studies, it appears to be fairly common to find transcripts generated from the DNA strands opposite to those on which ORFs are located. This has been termed “antisense transcription” [20, 24]. In our study, antisense transcripts were 7% of the primary transcriptome, as predicted by mapping TSS locations. Such transcripts were found to be 12% of the transcriptome of strain PA14 [24], while in H. pylori, they represented 27% of the transcriptome [20]. The large difference between these organisms might reflect different cutoff thresholds used for determining what qualifies as a TSS. It could, however, reflect the biological differences that exist between P. aeruginosa and H. pylori, or reflect the fact that the P. aeruginosa genome is 3-fold larger than that of H. pylori.

A RNA-Seq study of P. aeruginosa PA14 focused on the changes in gene expression between cells grown in planktonic vs. biofilm conditions [24]. This group interpreted the first base at the 5′ end of RNA-Seq read pileups as the TSS. They identified a total of 3389 putative TSS (1054 of which were present under more than one culture condition). Their study demonstrated a high degree of reproducibility between biological replicates as well as between different culture conditions, confirming that RNA-Seq can be used to detect expression of genes encoding essential proteins and proteins involved in housekeeping functions, as well as those genes that are environment-specific. However some differences were evident between this prior study and ours. Dotsch et al. [24] reported that 75% of TSS are upstream of start codons in strain PA14 based on conventional RNA-Seq, while here we determined by dRNA-Seq that only 55% of TSS are located upstream of start codons in strain PAO1. This might be due to differences in methodology, and we note that the dRNA-Seq method demonstrated similar proportions of transcripts (49%) beginning upstream of start codons in the distantly related ε-proteobacterium Helicobacter pylori [20]. The use of standard RNA-Seq data would make discerning the 5′ ends of reads that form peaks in the middle of actively transcribed genes very difficult due to the cDNA fragmentation process inherent in library construction and the nucleolytic degradation that occurs rapidly with labile prokaryotic RNA. In contrast, the dRNA-Seq methodology enables the identification of the 5′ ends of transcripts regardless of where they lie within a gene, and the method is not sensitive to nucleolytic degradation that might otherwise obscure TSS signals within genes. Similarly, the fragmentation process that occurs during library construction will not obscure the 5′ ends of the cDNA molecules. Therefore, dRNA-Seq data does not share the same compounding set of interpretation issues as standard RNA-Seq data.

Under the standard growth conditions used, only 34.5% of genes in PAO1 were found to have an upstream TSS, including genes found in predicted operons. In contrast, Toledo-Arana et al. report that under all conditions studied, the firmicute Listeria monocytogenes transcribes at least 98% of its genes [53]. This difference likely reflects the large metabolic diversity/flexibility evident in P. aeruginosa [1]. Its repertoire of 24 sigma factors [10] and nearly 10% of genes involved in transcriptional regulation, is large for a bacterium, and the diverse conditions under which P. aeruginosa lives requires the complex interplay of genes expressed under different conditions to ensure survival and competitiveness.

RNA-Seq data revealed novel genes and layers of transcriptional complexity

In addition to the novel layers of transcriptional complexity revealed by our dRNA-Seq analysis, this study detected other impacts on transcription. For example, data from the RNA-Seq library analysis indicated a trend towards a reduction in the level of transcription around the terminus of replication (Fig. 3). This is in concordance with Illumina DNA sequencing [33] and transposon mutagenesis studies [34] that also recovered decreased numbers of sequences and mutants respectively from this region of the genome. This phenomenon is likely related to the manner in which the genome replicates, since at any given time in an actively growing cell culture there would be more DNA present at the origin of replication than at the terminus. However, here it was demonstrated that this bias is also detectable at the transcriptional level with more transcripts evident from genes that are present in regions that are origin-proximal than terminus-proximal, likely due to a combination of gene dosage and the higher tendency for relaxation of supercoiling at the origin, which would impact gene expression. Genomic rearrangements occur more commonly in a symmetrical fashion around the terminus and origin, rather than between the terminus and origin regions [54], since such symmetrical rearrangements would conserve existing levels of gene expression by enabling genes to maintain the same distance from the origin of replication, and therefore the same copy number. In contrast it can be anticipated that there would be fitness costs to the organism if origin-proximal genes were relocated near to the terminus as supported by previous studies of detected genome rearrangements in P. aeruginosa [54].

Small, non-coding RNA (ncRNA, sRNA) transcripts are emerging as a major mechanism for regulating translational expression in bacteria and were also investigated here. At the time of the current study there were 140 annotated ncRNAs in PAO1, based on the Pseudomonas genome database [16] and the Rfam website [55]. All but 63 of these are rRNAs and tRNAs. Recent P. aeruginosa transcriptome investigations have sought to more thoroughly annotate ncRNAs. A recent RNA-Seq study [35] on strain PAO1 grown in LB at 37 °C into exponential and early stationary phase suggested an additional 513 ncRNAs in PAO1. However this was based on low-stringency, automated computational methods. The existence of several of these ncRNAs was verified by Northern blot. Conversely, others [17] identify 165 novel ncRNAs in their transcriptome analysis of P. aeruginosa PA14. Many of the ncRNAs defined in that study employed a very low threshold for TSS prediction (minimum number of 5′ triphosphate library reads aligned to a single genome position = 5), as well as a small number of total RNA-Seq reads spanning the region of the predicted ncRNA. Therefore, these ncRNAs predicted by others should be further verified using additional sequencing data, or an additional method such as Northern blotting or PCR, to confirm their existence and expression under different conditions, as well as looking at their potential regulatory roles in the cell. Our current investigation utilized high stringency methods and confirmatory RT-PCR to define 31 reliable ncRNAs that were produced under the investigated conditions, and not all of these were identified in the above-mentioned RNA-Seq studies. Intriguingly 30 of these were found to be dysregulated under the adaptive lifestyle conditions of biofilm formation and/or swarming motility. This indicates that it is important to consider different growth conditions when confirming ncRNA species and implies that translational regulation mediated through ncRNAs may be an important element in determining adaptive lifestyle changes. Certainly this is true of the ncRNA crc [56] and rsmYZ [57] and we have preliminary evidence implicating the importance of prrF1,2 and phrS as critical regulatory elements in one or both of these lifestyle changes [58].

Conclusions

Here, we use pRNA-Seq to identify RNA transcript 5′ monophosphate cleavage/processing sites in a genome-wide manner. Cleavage sites were predominately located between the first and second codon positions within protein-coding genes. Further examination of cleavage sites has revealed that they can be classified into five distinct categories, based on their cleavage peak shape and associated sequence motifs. The cleaved transcripts occur in genes associated with specific KEGG categories, but a much wider set of categories and genes was observed than initially anticipated. We also identified a correlation between TSS that lie ~ 18 bp apart on opposite strands of the transcriptome. These sites are separated by a distinct motif, and transcription may be initiated here by RNA polymerase dimers. This combination of pRNA-Seq, dRNA-Seq and RNA-Seq thus provided us with a more extensive view of the transcriptome.

Due to the diversity of lifestyles and complex adaptations that Pseudomonas can undertake, it would be necessary to perform similar analyses, under these conditions, to those described here in order to define the transcriptional complexity that underpins diversity in this organism. Nevertheless, this study has provided a new window into a previously unappreciated level of complexity in RNA processing in a bacterial transcriptome. It involves a greater extent of transcript cleavage than previously anticipated, evidently occurring in a regulated fashion through enzymatically-controlled processes. However, this is only the start, as we expand, in the future, our understanding of the role and significance of RNA processing events in maintaining a dynamic, flexible and robust bacterial transcriptome.

Methods

DNA extraction, genomic library construction and sequencing

Pseudomonas aeruginosa strain PAO1 was grown in Luria Broth (LB) medium at 37 °C (to OD600 = ~ 0.7 at ~ 200 rpm. DNA was extracted using a protocol modified from Cheng and Jiang [51]. Briefly, cells were collected by centrifugation at 4 °C, washed twice with STE buffer (100 mM NaCl, 10 mM TRIS-HCl, 1 mM EDTA, pH = 8.0), then resuspended in TE buffer (pH = 8.0). Cells were lysed by adding phenol and vortexing for 60 s. Chloroform phenol extractions were performed to extract DNA. DNA was precipitated with ethanol and sodium acetate. DNA was sheared and size fractionated on a 10% SDS-poly acrylamide gel electrophoresis (PAGE). The 190–210 bp region was excised, and DNA eluted and purified using a QIAquick purification kit (Qiagen). Libraries were constructed using the Illumina Genome Analyzer protocol and 50 bp paired-end sequence reads were obtained using an Illumina Genome Analyzer II according to the manufacturer’s instructions.

SNP/indel analysis

Whole genome alignments were performed against the reference strain PAO1 (NC_002516) genome using 4 different tools: Bowtie [59], BWA [60], mrsFAST [61] and SSAHA2 [62]. All default parameters were used with the exception of minimum and maximum insert size specifications of 50 and 1000 for Bowtie and SSAHA2, kmer = 13 and skip = 2 for SSAHA2, and an edit distance of 3 for mrsFAST. After read alignments, single nucleotide polymorphisms (SNPs) were identified using Samtools [63] version 1.12 and were filtered for SNP quality scores greater than 90, read depth greater than 50 and percentage of non-reference bases greater than 90%. All heterozygous calls were removed since only a single allele is expected for haploid genomes. Most predictions overlapped across the 4 different alignments for this highly filtered set of SNPs, with the exception of BWA calling an insertion in place of 2 consecutive SNPs.

RNA extraction

P. aeruginosa strain PAO1 was grown in LB medium at 37 °C (libraries A06027, 110817_SN865 and PA0004) or 34 °C (library PA0001), synthetic cystic fibrosis medium (SCFM) [64], and artificial sputum medium (ASM) [65] after inoculation from an LB culture grown overnight at 37 °C at ~ 200 rpm. For growth temperatures, sample cultures OD600 and rates of shaking, see Table 1. RNA was extracted using Qiagen’s RNA Protect Bacteria Reagent and RNeasy Midi kit (Qiagen) using the manufacturer’s protocol except that all centrifugation steps were carried out at 4 °C. A DNase digestion step was carried out on the columns as per the manufacturer’s protocol. After elution from the columns, the RNA was further purified using Trizol (Invitrogen) following the manufacturer’s directions. rRNA depletion was performed twice using Ambion’s MICROBExpress kit. The above protocol was used for cells grown in all media with the following exceptions: A second DNase digestion was conducted for the LB 34 °C and the SCFM samples. The OD600 for the samples grown in ASM media was not determined due to biofilm formation in this medium; these samples were grown for 48 h before harvesting. ASM cultures were stabilized by adding an equal volume of RNA Later (Ambion), incubated at 23 °C for 10 min, then centrifuged for 30 min at 3200 x g, 4 °C. The pellet was resuspended in Sputasol (Oxoid) in order to break down the biofilm structure and incubated at 37 °C for 25 min with shaking. Two volumes of RNA protect were added, then the mixture was incubated at 23 °C for 5 min and centrifuged for 30 min at 3200 x g, 4 °C. The supernatant was removed, and the pellet was further extracted using the RNeasy kit beginning at step #6 of the manufacturer’s protocol.

pRNA-Seq library construction, sequencing and read mapping

The pRNA-Seq library construction process is depicted in Additional file 1: Figure S7. The cleavage site enriched library A06027 (LB 37 °C) was prepared by first ligating a 17 base long adenylated DNA oligonucleotide (17.71; Additional file 1: Table S7) onto the 3′ end of the 2× DNase treated RNA using T4 RNA ligase. A second 17 base long ribonucleotide (17.50, see Additional file 1: Table S7) was then ligated onto the 5′ end of the RNA using T4 RNA ligase. These RNA/DNA hybrid molecules with known adapter sequences at the 3′ and 5′ ends were then reverse transcribed using 16 base oligonucleotide 16.16 (see Additional file 1: Table S7) as a primer and Superscript III reverse transcriptase (Invitrogen). The resulting cDNA was PCR amplified using oligonucleotides 16.16 (16 bases) and 17.53 (17 bases) (Additional file 1: Table S7) and resolved on a 6% denaturing-PAGE gel. A gel fragment containing cDNA fragments ranging in size from 600 bp - 10 kb was excised to remove any small fragments. The excised DNA was eluted and again PCR amplified using 16.16 and 17.53 to enrich for full-length cDNAs. The cDNA was then sheared and again size fractionated via SDS-PAGE. The 190–210 bp area (which corresponded to the desired library fragment size) was excised, eluted and purified with a QIAquick purification kit (Qiagen). Libraries were constructed using the Illumina Genome Analyzer protocol and paired-end 75 bp sequence reads were obtained using an Illumina Genome Analyzer II according to the manufacturer’s instructions. All reads were from the A06027 library were checked for passage of Illumina quality standards, then converted into FASTQ format. An in-house Perl script was used to trim off the 5′ ends of any reads from library A06027 that had the 17 bp adapter attached from the library synthesis process. The script (provided in Supplementary Information) allowed for up to three mismatches. The trimmed reads were then aligned using Bowtie [59] as for standard RNA-Seq reads. Samtools [63] was then used to separate the individual reads into two strand-specific files.

dRNA-Seq library construction, sequencing and read mapping

The RNA used to construct library SN865 (LB 37 °C) was prepared as described previously [20] by Vertis Biotechnologie (Germany). The library was a single end Illumina library, and 50 bp strand-specific sequence reads were obtained using an Illumina Hi-Seq 2000 machine according to the manufacturer’s instructions. All reads were quality checked according to Illumina standards and then converted into FASTQ format. Reads were aligned to the PAO1 genome using Bowtie [59] except that the -X 1000 command was not used as this was a single-end library. Samtools [63] was also used to separate the reads into two strand-specific alignment files.

RNA-Seq library construction, sequencing and read mapping

The RNA used to construct libraries PA0001 (LB 34 °C), PA0004 (LB 37 °C), A03674 (SCFM) and A06026 (ASM) was reverse transcribed using random hexamers and the SuperscriptTM Double Stranded cDNA synthesis kit (Invitrogen). cDNA was sheared and size fractionated using SDS-PAGE. The 190-210 bp area was excised, eluted and purified with a QIAquick purification kit (Qiagen). Libraries were constructed using the Illumina Genome Analyzer protocol and paired-end 50 bp sequence reads were obtained using an Illumina Genome Analyzer II as per the manufacturer’s instructions. All reads were checked for passage of Illumina quality standards. Reads organized into FASTQ files were aligned using Bowtie [59] to the PAO1 genome using –X 1000 such that only mate pairs were reported if separated by less than 1000 bp. All other settings were the defaults. Once aligned, Samtools [63] was used to remove duplicates and select for reads that were aligned in proper pairs. The number of reads aligned to RNA genes were summarized using coverageBed (part of the bedtools software package; [66]). Reads per kb/million reads (RPKM) was calculated as a measure of expression of all genes individually under all conditions [67] using the formula: RPKM = number of mapped reads / total number of reads / gene length × 1000,000,000.

Cutoff identification: Peak height threshold determination for TSS and RNA cleavage site predictions

Please see Supplemental Information for a detailed description of how peak height cutoffs were determined for the dRNA-Seq and pRNA-Seq data.

Prediction of TSS and cleavage sites

Once mapped, mate pairs of the 5′ monophosphate reads that were from nuclease-cleaved ends were discarded. Remaining reads from both the pRNA-Seq and dRNA-Seq libraries were trimmed of adapter sequence up to the first base pair. Genomic locations with coverage over the statistically determined cutoff (100 reads for cleavage sites/pRNA-Seq, 500 reads for TSS/dRNA-Seq as per the threshold rationale above) were marked as a peak, which represents possible cleavage sites or TSS, respectively. Regions encoding rRNA or tmRNA and their respective upstream intergenic regions were subjected to a higher cutoff threshold since the coverage in these regions was disproportionally higher. For the pRNA-Seq library, the peak height cut-off on both strands in these areas was 2000. For the dRNA-Seq library, the cutoff on the + strand was 30,000 and on the – strand was 2000. If multiple peaks were observed within 5 nucleotides of one another, only the location with the highest coverage was used as a predicted cleavage site. The TSS peaks were categorized according to their locations. Peaks in an intergenic region and on the same strand as the closest downstream gene = primary. Peaks within gene boundaries and on the same strand as the gene = internal. Peaks within gene boundaries and on the opposite strand from the gene = antisense. Peaks within an intergenic region = primary antisense and peaks within an area where 2 gene boundaries overlap = internal genes overlap.

Cleavage site motif identification and function analysis

Cleavage site motifs were calculated using MEME software [68] with default parameters [68-86], except for utilizing the option to use the sense strand only, using sequence spanning from − 10 to + 10 bp of cleavage sites. The motif length was set incrementally from 6 bp to the full length of the sequence (21 bp). The motif search did not include the reverse complement of the extracted sequence and sites within 10 bp downstream of one another were removed from the analysis, since the shared subsequence would be confusing in identifying the motif. Once an initial motif was found for all cleavage sites, sites were binned into different peak shape categories for further analysis. Peak shape refers to the number of mapped transcripts surrounding a cleavage site. A window of ±5 bp was used to profile peak shapes. Peak shapes were subsequently categorized using k-means clustering with the R software package http://www.r-project.org/). The parameter k was estimated to be around 15 by plotting within-group variance for a number of the clusters. Therefore an initial k value of 20 was used and clusters with similar profiles were merged for subsequent analysis. An over-representation gene function analysis was calculated using hypergeometric tests based on KEGG pathway categories [29] and Holm’s test was used to correct for multiple testing.

Novel regulon member identification, RpoN binding site identification, Ribosomal binding site identification, small ncRNA identification and Functional analysis of ncRNAs by RT-PCR methods are described in Supplementary Information.

Supplementary results, figures and methods. Figure S1. Measures of cleavage site motif category sequence member distance within and between categories. Figure S2A and B. Additional cleavage site motifs. Figure S3. While our data is more comprehensive, 86% of previously published transcription start sites lie within ±3 nt of our dRNA-Seq set. Figure S4. A Promoter motif associated with several virulence factors and two hypothetical proteins is highly similar to the RpoN binding motif. Figure S5. Small non-coding RNA gene expression in biofilms (A) and cells undergoing swarming motility (B) compared to expression in planktonic cells. Figure S6. Proposed polymerase dimer binding to a palindromic antiparallel TSS. Figure S7. Step by step construction of the pRNA-Seq library. Figure S8. Peak height threshold determination for generating cutoffs for site predictions (Transcription start sites (TSS) and RNA cleavage sites). Table S1A. SNPs found in our sequenced isolate of PAO1 cf. the reference strain. Table S1B. Indels found in our sequenced isolate of PAO1 vs. the reference strain. Table S2. The 240 most abundant transcripts from the RNA-Seq study of cells grown in LB at 37 °C. Table S3. Sequence motifs associated with cleavage sites. Table S4. p-values for promoters matching motif associated with virulence factors.Table S5A. Summary of predicted RpoN binding sites and adjacent downstream genes in Pseudomonas aeruginosa PAO1. Table S5B. Function classes associated with genes downstream of predicted RpoN binding sites. Table S6. Fold change in RNA-Seq read density throughout the PAO1 genome. Table S7. Sequences of the oligonucleotides used in the preparation of the pRNA-Seq library. Table S8. Fitted slopes and intercepts for data from H = 1..100XX. Fit to G(H) = H0XS. Table S9. List of primers utilized for confirmation of sRNAs by RT-qPCR. The Perl script used for 5′ end trimming on library. (DOCX 1536 kb)