PNPase is involved in the coordination of mRNA degradation and expression in stationary phase cells of Escherichia coli

mRNAs are very stable in stationary phase

The deletion strains of RNase R or PNPase (the rph-1 Δrnr and rph-1 Δpnp double mutant strains, respectively) were constructed in the E. coli K-12 MG1655 background referred as the rph-1 control strain. The E. coli rph-1 control strain and the rph-1 Δrnr double mutant displayed similar growth profiles in LB medium while the rph-1 Δpnp double mutant led to slower growth. Since this work aimed at comparing the role of PNPase and RNase R exoribonucleases during the stationary phase of growth, samplings for genome scale mRNA stabilities were performed at late stationary phase in standardized conditions, i.e. 3 h after the growth stopped (Additional file 1: Figure S1). Using our calculation and quality control method (see material and methods for details), 2856 half-lives were available for all the three strains. This dataset was retained for further analysis (Fig. 1a). Half-lives have been determined with classical linear model fitted on log concentration of mRNA (see example of a set of 50 randomly selected mRNAs in Additional file 2: Figure S2). The half-lives obtained (Fig. 1, median value > 13 min) were higher than previously published in the rph-1 control strain and the rph-1 Δpnp double mutant [2, 29, 34], however our growth conditions were different (prolonged stationary phase versus exponential phase), and this could justify in part the discrepancy in the results. We also hypothesized that a delay in the transcription arrest during the mRNA half-life measurement could lead to an increase in the determination of the general values of half-lives. Such a delay, probably due to RNA polymerase elongation activity in presence of rifampicin, was first demonstrated by Chen and coworkers [34] and modelled more recently [35]. We have subsequently re-estimated the half-lives taking into account the delay using a model composed of a stable baseline level during the delay followed by an exponential decay (as described by [35]). The boxplot obtained with the delay (Fig. 1b) is similar to the one obtained without delay, meaning that the half-life values are very similar with and without delay (see correlation coefficient close to 1 in Additional file 3: Figure S3).

The values of half-life were confirmed by Northern-blot experiments with 3 different mRNAs (osmB, ompA and nlpl) measured in the 3 strains. Although the absolute values differed, the two sets of mRNA half-lives were in the same range of magnitude (Additional file 4: Figure S4). A perfect match was even found for the half-life of osmB in the 3 strains.

The high mRNA half-lives obtained for the control strain and for the two exoribonuclease mutant strains were not related neither to the measurement method nor to the fitting estimation. These values were thus rather linked to the E. coli physiology in prolonged stationary phase. Here we demonstrated that mRNAs are very stable in non-growing cells. Stabilization of mRNAs at zero growth is in agreement with the recent demonstrations of mRNA stabilization at low growth rates [2, 34, 36]. Overall our results show an important stabilization of mRNAs associated to E. coli adaption to unfavorable growth conditions.

Only a few mRNAs are stabilized in the rph-1 Δrnr and rph-1 Δpnp double mutants

It was expected that mRNAs would be stabilized at the genome-wide level in the absence of the ribonucleases RNase R or PNPase. The half-lives in the rph-1 Δrnr and rph-1 Δpnp double mutants were compared to the rph-1 control strain at the global level (comparison of statistical distribution of data) and for individual mRNAs applying the statistical test defined in the material and methods section. In the rph-1 Δrnr double mutant, the median of the mRNA half-lives was 28.3 min while for the rph-1 control strain the median of the half-lives was 24.2 min (Fig. 1). Between rph-1 Δrnr and the rph-1 control there was no mRNA with a significant variation of stability suggesting only a small systematic increase in global mRNA stability in the rph-1 Δrnr double mutant. In the rph-1 Δpnp double mutant there is a global destabilization (see the lower mean and median half-lives compared to the rph-1 control in Fig. 1) but 214 mRNAs were found to be significantly stabilized. These stabilized mRNAs belonged to several functional categories. Functional enrichment analyses processed with AMIGO or KEGG Mapper revealed a significant enrichment of mRNAs involved in the central carbon metabolism. The functional categories “glycolysis” and “TCA cycle” were indeed enriched. More specifically, important genes of these processes were stabilized in the rph-1 Δpnp double mutant (pgi, pfkA, pgk, pck, acp, gltA, sdaA, sucA, sucC, ldhA, agp, talB and yccx). Interestingly, Bernstein et al. similarly reported the stabilization of pgi, pfkA, pgk, pckA, gltA and talB in a rph-1 Δpnp double mutant strain compared to a rph-1 strain in exponential phase [29] (data for the other genes mentioned here were not available). Recent papers also pointed out the link between PNPase and the regulation of the expression of genes related to the central carbon metabolism [33, 37].

We did not observe any massive stabilization in the rph-1 Δrnr and rph-1 Δpnp double mutants, neither taking in consideration the results of mRNA stability all together nor at the individual scale.

mRNAs are extensively destabilized in the rph-1 Δpnp double mutant

Many mRNAs displayed decreased half-lives in the rph-1 Δpnp double mutant compared to the rph-1 control strain. In comparison only five mRNAs, namely osmB, cstA, uspB, ycgB and ychH, were destabilized in the rph-1 Δrnr double mutant. Deletion of pnp thus led to an overall loss of stability, suggesting that the PNPase could be involved in a protective mechanism of mRNAs.

Massive reorganization of genome-wide mRNA levels in the rph-1 Δpnp double mutant

Modification of RNase activity in the different RNase mutants is expected to directly change mRNA stabilities, therefore altering mRNA concentrations. However, many indirect effects on mRNA concentrations are also likely to occur, for instance when RNases targeted the stability of transcriptional regulators [32, 33]. We have notably found in this work that 70 transcriptional regulators (among the 162 identified according to RegulonDB) exhibited a modified stability in the rph-1 Δpnp double mutant (Additional file 5: Table S1). Therefore, in order to decipher the role of RNase R and PNPase, we have estimated the genome-scale mRNA concentrations by transcriptomic analysis and compared the data in the different strains. Raw data of steady-state transcriptome before addition of rifampicin were recomputed with a RMA-derived method (see material and methods for details). Dramatic changes were observed in the rph-1 Δpnp double mutant whereas only a minor transcriptional reorganization occurred in the rph-1 Δrnr double mutant (Fig. 2). Indeed, applying the statistical criteria as defined in material and methods, only 27 genes were differentially expressed between the rph-1 Δrnr double mutant and the rph-1 control strain whereas 2753 displayed significantly different levels between the rph-1 Δpnp double mutant and the rph-1 control. Most of the genes were up-regulated in the rph-1 Δpnp double mutant compared to the rph-1 control strain (Fig. 2). Our transcriptome results are in agreement with recently published RNA-Seq data showing that, in exponential phase, PNPase deletion led to more up-regulated transcript concentrations compared to RNase R deletion [33]. Such changes in mRNA levels in the rph-1 Δpnp double mutant were corroborated by the 4-fold increase of the total mRNA concentration and by the 2-fold increase of the total RNA (mRNA, rRNA and tRNA) content in this strain in comparison to the rph-1 control strain (Additional file 6: Table S2).

Overall, the deletion of PNPase led to a massive reorganization of the gene expression most probably due to indirect transcriptional regulations. The global up-regulation of mRNA concentrations observed in the rph-1 Δpnp double mutant is in the opposite direction to the down-regulation of mRNA stabilities also found in this strain.

PNPase affects the relationship between mRNA stability and concentration

The role of exoribonuclease activities was thus complex, combining effects on both mRNA stability and mRNA concentration. This suggests a delicate balance between regulations of mRNA stability and concentration in cells. Until now, we have compared the behavior of genome-wide mRNA concentrations and half-lives between two strains. We focus now on the relationship between mRNA concentrations and half-lives within the same strain. The overall relationship between mRNA stability and concentration characterizes the in vivo activity of the degradation machinery (at the level of enzyme-substrate interaction) [42]. In order to determine the impact of PNPase or RNase mutations on the in vivo degradation activity, we have analyzed the stability - concentration relationship for all mRNAs of each strain using a genome-wide correlation analysis. We have combined on Fig. 3, the three plots of the degradation rate constant k (k = ln2/ half-life) versus the mRNA concentration obtained for the rph-1 control strain and the rph-1 Δrnr and rph-1 Δpnp double mutants, respectively. This figure allowed the visualization of a global negative relationship between degradation rate constant and the mRNA concentration for each individual strain but the level of correlation differed with the mRNA concentration. Three correlation phases were differentiated using a K-means clustering algorithm on the combined three datasets (see the correlation coefficients and p-values in the legend of Fig. 3). In phase I, for the lowest mRNA concentrations, the correlation was very strong; in phase II for intermediary mRNA concentrations, the correlation was less strong and in phase III for the highest mRNA concentrations, the correlation was not significant. Phases determined with combined data of the three strains were also relevant for individual strains (see correlation coefficient for individual strains in the three phases in Table 2). For each of the three strains, most of the genes showed a negative correlation between the degradation rate constant (k) and the mRNA concentration (85, 88 and 72% of genes were in phases I and II for the rph-1 control strain and the rph-1 Δrnr and rph-1 Δpnp double mutants, respectively) (Table 2). This result shows that in a strain mRNA stability and concentration are generally tightly linked in vivo. However in the rph-1 Δpnp double mutant, a more than 2-fold higher number of genes with uncorrelated degradation rate constant and concentration was observed (1144 genes in the rph-1 Δpnp double mutant were in the phase III against about 400 genes in the two other strains) (Table 2). Most of the genes in phase III in the rph-1 Δpnp double mutant (779) belonged to phases I or II in the rph-1 control strain. This was not the case in the rph-1 Δrnr double mutant since nearly all the genes found in phase III were also found in this phase in the rph-1 control.

	Data	PHASE I	PHASE II	PHASE III
rph-1 control	Number of genes	1889	834	476
	Correlation coefficient	−0.80	−0.60	0.27
	p-value	4.14E-108	1.07E-16	9.87E-09
rph-1 Δrnr double mutant	Number of genes	1957	750	369
	Correlation coefficient	−0.80	−0.69	0.04
	p-value	9.66E-107	1.88E-23	4.23E-01
rph-1 Δpnp double mutant	Number of genes	1150	1747	1144
	Correlation coefficient	−0.83	−0.70	0.28
	p-value	8.75E-121	1.03E-23	3.29E-09

mRNA concentration was thus negatively correlated to stability for a large majority of mRNAs (phase I and II) even in the rph-1 Δpnp double mutant which displayed an unexpected swing in mRNA concentrations and stabilities when compared to the rph-1 control strain. However, the number of mRNAs in phase III with unrelated stability and concentration values increased significantly in the rph-1 Δpnp double mutant. PNPase mutation thus affected significantly the relationship between mRNA concentration and stability unlike the RNase R mutation. This suggests that PNPase is a key enzyme of the RNA degradation machinery and the implication on the regulation of the in vivo kinetics of mRNA degradation will be presented more in details in the discussion section.

Strains and growth conditions

The double E. coli mutant strains with rnr or pnp deletion (rph-1 ∆rnr and rph-1 ∆pnp, respectively) in the MG1655 background (F⁻ λ⁻rph-1) were respectively obtained by phage transduction from donor strains in the MG1693 background (F⁻ λ⁻rph-1) [12, 32, 51]. P1-mediated transductions were performed as previously described [52]. 50 μL of an overnight culture of the donor strains were used to inoculate 5 mL of LB with 0.4% glucose and 10 mM CaCl₂. After 30 min of growth at 37 °C, phage P1 lysate (100 μL) was added. After 2 h incubation at 37 °C, cell lysis was performed by chloroform (100 μl) addition. Cells were centrifuged at 13500 rpm for 10 min, the supernatant (transducing lysate) transferred to a new tube with 30 μl chloroform and kept at 4 °C. An overnight culture of the recipient strain was centrifuged and the pellet resuspended in 2.5 ml LB with 10 mM MgSO₄ and 5 mM CaCl₂. Cells (100 μl) were incubated with 10 to 200 μl of transducing lysate during 30 min at 30 °C. 300 μl of sodium citrate 1 M were added and the cells were incubated another hour at 30 °C. Infected cells were centrifuged and resuspended in 100 μL of the supernatant. Transductants were selected by plating on LB agar plates containing kanamycin (50 μg·ml^− 1) for the rph-1 ∆rnr strain and streptomycin (50 μg·ml^− 1) for the rph-1 ∆pnp strain. The rph-1 control strain and the rph-1 ∆rnr and rph-1 ∆pnp double mutants were grown in baffled flasks in LB medium at 37 °C and 180 rpm. Initial pH was set to 7. All cultures were inoculated at a low OD of 0.1 to maintain cells in exponential phase of growth for many generations prior to entry into stationary phase to dilute out any RNAs and proteins provided by the inoculum. Inoculation was performed from overnight pre-cultures performed in similar conditions. Biomass was estimated from absorbance at 600 nm. Each culture was repeated three times to provide independent biological replicates.

RNA extraction protocol

Three hours after the stationary phase of growth was reached, samples were collected for transcriptomic analysis; this time point was also the reference time point (T0) for the half-life determination procedure. Subsequently, rifampicin (500 μg.mL^− 1) was added to inhibit the initiation of transcription, and cells were harvested at three different time points after this addition. Cultures were performed in triplicate. Samples were taken at either at 0′, 2′, 8′ and 15′, or 0′, 5′, 10′ and 18′ or 0′, 3′, 12′ and 30′ min after rifampicin addition. Samples were immediately frozen in liquid nitrogen upon collection. After thawing and centrifugation steps, total RNA was extracted with TRIZOL® Reagent (Ambion) according to the manufacturer’s instructions. DNA contamination was eliminated with Turbo DNase kit (Ambion). Total RNA concentration and integrity were measured using a Nanodrop® spectrophotometer and Agilent BioAnalyzer, respectively. Total RNA extraction profiles were checked to be similar for all the three tested strains.

Microarray procedures

A double-stranded cDNA synthesis kit (InvitroGen) was used to produce cDNA from 2 μg aliquots of total RNA. Aliquots of 1 μg of cDNA were labeled using the one color DNA labeling kit and 2 μg of labeled cDNA were hybridized onto E. coli K-12 gene expression arrays (Nimblegen, Roche) for 17 h at 42 °C according to the manufacturers’ instructions. Arrays were washed and then scanned with a MS200 Microarray Scanner (Nimblegen, Roche). The images were analyzed with DEVA 1.2.1 software. Only raw data were used for further analyses. All array procedures were performed by the GeT-Biopuces platform (http://get.genotoul.fr).

Transcriptome and genome-wide mRNA half-life measurements

For transcriptomic analysis, raw probe intensities were processed and analyzed with the R computing environment using the affy [53] and limma [54] packages of Bioconductor. Raw data were submitted to a RMA-based background correction [55]. After background correction, intra-replicate quantile normalization was performed for each strain. A set of probes in the background for which the ranks were roughly invariant across all nine arrays was selected. The median value of the invariant probe set intensities in each condition was used as a scaling factor for normalization between the three strains. After normalization, the intensity of a transcript was calculated by a RMA-summarization procedure [55] within each condition. Intensity values were multiplied by the total RNA extraction yield (in μg total RNA per mg of dry cell weight) to provide the mRNA concentration value in arbitrary units per mg of dry cell weight. RNA extraction yields were 15.6 ± 3.0, 17.2 ± 2.3 and 35.9 ± 7.3 μg RNA per mg of dry cell weight for growth of the rph-1 control strain and the rph-1 ∆rnr and rph-1 ∆pnp double mutants, respectively. The multiplication step by the total RNA extraction yield allowed to take into account differences in total RNA content (rRNA, tRNA and mRNA) per cell weight between the strains. Differences in mRNA concentration were evaluated with a modified t-test in conjunction with an empirical Bayes method [54]. The p-values were adjusted for multiple testing by the “BH” False Discovery Rate (FDR) method [56]. A p-value threshold of 1% and a fold change higher than 2 or lower than 0.5 were used for significance of differences in mRNA concentration.

For mRNA half-life determinations, twelve arrays (three reference T0 samples, and nine time points after addition of rifampicin) were used. Only the normalization between arrays according to the invariant probeset intensities was performed. In each array, transcript-specific intensity was computed as the median value of the 16 targeting probe intensities. The linear regression coefficient, k, of ln(mRNA) versus time (12 points) and its associated coefficient of variation (standard error of slope/estimation of slope) were calculated for each mRNA species. The determination of k was considered as reliable only if the associated coefficient of variation was below 30%. The linear regression coefficient k corresponding to the degradation rate constant was inversely proportional to the mRNA half-life t_1/2, \( k=\frac{\ln 2}{{\mathrm{t}}_{1/2}} \). The statistical significance of differences in half-life was evaluated using the probability value of interaction between time and growth rate in a global model of linear regression. A statistical threshold of 10% was used for adjusted p-values by the “BH” FDR method [56].

For functional analyses, GO and KEGG enrichment analyses were performed on the differentially expressed genes as well as on the genes with differential stabilities. Enrichment significance was set with a cut-off of 5% for the associated p-value.

ncRNAs were quantified by RNA-Seq technology in the rph-1 Δpnp double mutant and the rph-1 control strain

RNA samples (20 μg) were sent to Vertis Biotechnologie AG, Germany, for library preparation and sequencing of libraries using an Illumina HiSeq platform (single end, 50-bp read length). For the library preparation Vertis Biotechnologie AG depleted the ribosomal RNA molecules from the total RNA preparations using the MICROBExpress Bacterial mRNA Enrichment Kit (Ambion). The rRNA depleted RNAs were then fragmented with RNase III and the 5’PPP structures were removed using RNA 5’ Polyphosphatase (Epicentre). Afterwards, the RNA fragments were poly(A)-tailed using poly(A) polymerase and a RNA adapter was ligated to the 5′-phosphate of the RNA fragments. First-strand cDNA synthesis was performed using an oligo(dT)-adapter primer and M-MLV reverse transcriptase. The resulting cDNA was PCR-amplified to about 30 ng/μl using a high fidelity DNA polymerase and sequenced. The reads were mapped against E. coli genome (NC_000913 downloaded from NCBI genome database). The ncRNAs expression was quantified using Artemis [57] and the Log2 Fold change was calculated for each of the ncRNAs.

Confirmation of mRNA half-life data by northern blot

Three genes osmB, ompA and nlpl were selected for additional calculation of their associated mRNA half-life using Northern blot method. Northern blots were performed using the total RNA samples extracted for microarrays. Briefly, 10–30 μg of total RNA were fractionated under denaturing conditions in 1.2% agarose formaldehyde gel in MOPS buffer. RNAs were transferred onto Hybond-N+ membrane (GE Healthcare) and cross-linked by UV irradiation using a UVC 500 apparatus (Amersham Biosciences). Membranes were hybridized overnight with radiolabeled specific probes in PerfectHyb Plus Hybridization Buffer (Sigma Aldrich) at 68 °C. Specific probes (primers are given in Additional file 4: Figure S4) were obtained by in vitro transcription reactions with PCR DNA templates carrying a T7 promoter sequence through the use of 32P-UTP and T7 RNA polymerase (Promega). Radiolabeled probes were purified on G25 Microspin columns (GE Healthcare).

Availability of supportive data

The data discussed in this publication have been deposited in NCBI’s Gene Expression Omnibus [58] and are accessible through GEO SuperSeries accession number GSE116652 and the GEO Series GSE60107 for microarrays and RNA-Seq data, respectively. Other supportive data are included as additional files.