Physiological consequences of altered PNPase activity
To analyze the functionality of PNPase in vivo, we designed and cloned a pnp mutant strain of Rhodobacter sphaeroides 2.4.1. The KH-S1 RNA binding domains were removed and a stop codon was introduced at the end of the remaining coding sequence of pnp resulting in a truncated enzyme lacking those domains. The knockout was confirmed via selection on agar containing gentamicin and subsequent RNA sequencing analysis (Fig. 2a). Growth behavior of this strain differed from that of the wild type (Fig. 2b). When cultivated under microaerobic conditions, the growth rate was reduced, but both wild type and mutant finally reached the identical OD660. Under phototrophic conditions the mutant leaves exponential phase earlier than the wild type reaching a lower final OD660. A previous study revealed that reduced RNase E activity strongly impeded phototrophic growth of R. sphaeroides, while it had no effect on chemotrophic growth [7]. Moreover, the pnp mutant and the parental wild type strain vary in pigment composition (Fig. 2c+d). These differences are strongly dependent on the cultivation conditions: A significantly lower concentration of carotenoids and bacteriochlorophyll a was observed in the pnp mutant under microaerobic conditions (p-values < 0.05), while the pnp mutant exhibited repeatedly higher pigment concentrations under phototrophic conditions. However, this difference was statistically not significant.
In E. coli, Yersinia enterocolitica and Photorhabdus sp. PNPase plays an important role in the cold shock response due to selective degradation of mRNAs for cold shock proteins at the end of the acclimation phase to low temperature [43,44,45,46]. Based on this observation, we decided to test the R. sphaeroides strains for their ability to adapt to low and high temperatures. Wild type and pnp mutant cells were incubated at 4 °C or 42 °C on agar plates for 1 day and then shifted to an optimal temperature of 32 °C. In both cases growth of the pnp mutant was strongly impeded, while the wild type was able to grow at 42 °C and 4 °C (Fig. 2e). Also, in contrast to the wild type, the pnp mutant was not able to grow on malate minimal agar containing 300 μM tBOOH, while the wild type showed weak growth. Tertiary butyl-alcohol is representing organic peroxides that are produced e. g. during photo-oxidative stress.
Our results show that PNPase of R. sphaeroides is involved in cold adaptation as other bacterial PNPases and also is strongly impeded in its adaptation to heat. Whether the same molecular mechanisms are responsible for the phenotype as in other bacteria remains to be elucidated. This study for the first time analyses the function of PNPase in a phototrophic bacterium. The effect of PNPase on the bacteriochlorophyll levels and on carotenoid levels depends on growth conditions. Many genes are involved in the formation of photosynthetic complexes and it is not possible to correlate these phenotypic changes to specific changes of the transcriptome. We observed before that a temperature-sensitive variant of RNase E had little effect on growth under microaerobic conditions but strongly impeded phototrophic growth [7]. For the PNPase mutant we observed slower growth under both conditions, phototrophic growth was less affected, in contrast to the rne mutant.
PNPase modulates the transcriptome of R. sphaeroides
PNPase is an enzyme involved in many RNA processing reactions, and a global influence on the transcriptome can be expected as also shown for the Gram-positive S. pyogenes [31]. For the transcriptome analysis, three pre-cultures of the wild type and the pnp mutant strain of R. sphaeroides were inoculated with cells from three different single colonies. With each of these pre-cultures, three main cultures were inoculated (nine in total), grown under microaerobic conditions and later harvested during the exponential growth phase. All cultures initially derived from one colony in the first step were pooled. Total RNA was isolated and the DNA-free RNA was sequenced on an Illumina NextSeq 500 platform. The overall reproducibility within the replicates was fair, only one replicate obtained from the wild type strain showed some deviation to the other samples of the group (Supplementary Fig. S1, Additional file 1). In total 98% of the entire variance can be explained by the first two principal components.
Figure 3 shows the result of the DESeq2 analysis (version 1.26.0 [40];) and illustrates the log2-fold changes of the normalized read numbers in the pnp mutant versus the wild type strain (see Supplementary Table S3, Additional file 2). All transcripts with a log2-fold change ≤ − 1 or ≥ + 1 and an adjusted p-value ≤0.05 (Benjamini Hochberg algorithm) were considered to have a significant differential abundance within the two strains (coloured dots). We then decided to only keep those differentially expressed genes which have a basemean ≥50 (red dots) in order to further decrease the number of false positive hits. In total 334 transcripts met these strict criteria, 226 of them showed lower abundance in the pnp mutant strain and 108 showed higher abundance in the pnp mutant strain compared to the wild type. The most prominent differences were observed in the feature classes tRNA and rRNA: 94% of all tRNAs (51 out of 54) and 100% of all rRNAs (9 out of 9) showed a lower abundance in the pnp mutant strain (Fig. 3b). Altogether 37% of all non-coding RNAs, here merged of sRNAs and ncRNAs (including 6S, SRP RNA and tmRNA), were observed to have a differential abundance. Within the groups of RNAs with increased or decreased abundance, no distinct orthologous group of encoded proteins (COG) could be found to be prominent (Supplementary Fig. S2A + B, Additional file 1).
The transcriptome is directly affected by the action of RNases. Moreover, the RNA entity is modulated through secondary effects by the PNPase-mediated processing of sRNAs and mRNAs that code for regulatory elements, for example transcription factors. Thus, our transcriptome analysis reflects both direct and indirect PNPase dependent regulations and does not allow a distinction. In either case, our data emphasize the effect which PNPase has especially on stable RNAs (rRNA, tRNA). A similar effect was also observed in E. coli, although both rRNAs and tRNAs were more abundant in the pnp mutant despite a conducted rRNA depletion prior to RNA sequencing [47]. Further, Płociński et al. [48] demonstrated, that PNPase is involved in processing of ribosomal RNA and tmRNA in Mycobacterium smegmatis and M. tuberculosis.
We further validated these predictions using a different algorithm. An empirical Bayes approach integrated in the baySeq package (version 2.20.0 [41];) was used to identify differential expression (Supplementary Table S4, Additional file 2). The results of the two methods perfectly agree, since virtually all genes could be properly assigned. Every transcript (blue dot) with a log2-fold change ≤0 (pnp mutant/wild type) according to the DESeq2 analysis was also observed to be lower abundant in the pnp mutant according to the baySeq algorithm and vice versa (Fig. 3c). This includes every differently expressed gene which fulfills the strict criteria as mentioned above.
The RNA sequencing data was further used to investigate the cellular RNA 3′ elongation. All reads that could not be mapped end-to-end were instead mapped in very sensitive local mode with bowtie2 (version 2.2.6). To increase the quality of the analysis, all reads without a minimal matching sequence of 10 nt at the 5′ end were excluded. Only soft clipped sequences at the 3′ ends of the remaining reads were extracted with awk (version 4.1.3) (Supplementary Fig. S3A, Additional file 1). The overall results are similar for both strains: The lengths of elongated sequences are comparable, the majority of them (95%) is shorter than 39 nt in length. Further, the base frequency for each nucleotide position of the 3′ tail reveals an enrichment of guanine within the first 20 bases (Supplementary Fig. S3B + C + D, Additional file 1). A sequence motif which is related to a PNPase-dependent elongation could not be identified. Since both the lengths and base frequencies of the 3′ tails do not differ in between the analyzed strains, we conclude that the deletion of the KH-S1 domains does not have a major impact on the overall RNA 3′ elongation events in R. sphaeroides.
Levels of regulatory sRNAs are influenced by PNPase
An important effect of the PNPase on levels of small RNAs was reported: the enzyme does not only influence mRNA but also sRNA stability [49,50,51]. We were especially interested in those sRNAs that are derived from 5′ or 3′ UTRs and wanted to investigate the role of PNPase during the maturation process. For further analysis, we selected five sRNAs which showed a different pattern in the read coverage comparing pnp mutant and wild type. Two of them, CcsR1 and SorY, are known to have a regulatory function during the oxidative stress response in Rhodobacter sphaeroides [52, 53]. UpsM is processed from the mraZ 5′ UTR in a stress-dependent manner by RNase E [54]. The other two sRNAs have not been described so far and their function is still unknown. One is located in the intergenic region between RSP_1711 and rpsL and is derived from the rpsL 5′ UTR. The second one is derived from the 5′ UTR of RSP_6083. During the exponential growth phase, three of these sRNAs differed in abundance comparing the total RNA from the pnp mutant and the wild type strain (Fig. 4a+b). Moreover, processing products of the sRNAs IGR_1711_rpsL and 5′ UTR_6083 were prominently enriched in the pnp mutant. Interestingly, the abundance of the mature transcript of SorY and 5′ UTR RSP_6083 does not vary between the strains. To further evaluate the sRNA stability, we added rifampicin during the exponential phase and determined the RNA half-lives (Fig. 4c+d). CcsR1, SorY, IGR_1711_rpsL and 5′UTR_6083 are strongly stabilized in the mutant lacking PNPase, resulting in prolonged half-lives. In contrast to that, the half-life of UpsM drops form 12.2 min in the wild type to 4.0 min in the pnp mutant. The changed stabilities are in agreement with the observed sRNA levels during exponential phase (Fig. 4a). These observations highlight the role of PNPase during the maturation of sRNAs and in the homeostasis of their levels which has been described in E. coli. Cameron and De Lay [50] reported a stabilizing function by PNPase on Hfq-dependent sRNAs during the exponential growth phase in E. coli. They speculate, that PNPase may for example protect Hfq-bound sRNAs by degrading binding sites for other ribonucleases. Our data show a different trend in Rhodobacter sphaeroides, which suggests a mainly destabilizing effect on the selected sRNAs in this study. Even though CcsR1–4, UpsM and SorY are Hfq-dependent [52, 54, 55], CcsR1–4 and SorY are destabilized by PNPase and only UpsM fits the model proposed for E. coli. The two UTR-derived sRNAs are also destabilized by PNPase. The reason for these differences remains unidentified. Given this major influence of PNPase on the abundance of regulatory sRNAs, the pleiotropic effect of a pnp mutant becomes even more perspicuous.
Deletion of the KH-S1 domains of PNPase leads mainly to enriched RNA 3′ ends
As a 3′-to-5′ exoribunuclease, PNPase plays an important role in RNA turnover and decay from the 3′ end. Therefore, we analyzed how the RNA 3′ ends differ in abundance comparing the pnp mutant to the wild type. For this study we developed XPEAP, an analysis pipeline. It allows the detection of RNA 5′ or 3′ ends in prokaryotic NGS data and covers all relevant steps form data preprocessing to the final statistical analysis. As input the raw read files, three replicates of each strain, were used. After trimming and read alignment to the reference genome, READemption’s subcommand coverage was used to generate coverage files which contain the nucleotide positions of the 3′ end bases of each aligned read. The coverage files from the plus and minus strand were integrated in one data set. Subsequently all nucleotide positions that did not exhibit a coverage of ≥10 in at least one of the analyzed libraries were excluded. We further calculated the ratio of the 3′ end coverage and the full read coverage for every replicate and every position. To improve the signal noise ratio, only positions with a ratio higher than 0.05 were kept. A DESeq2 analysis was conducted based on this nucleotide-wise coverage files to detect differences within the two strains. All positions which showed a log2-fold change ≤ − 1 or ≥ + 1 and an adjusted p-value ≤0.05 (Benjamini Hochberg algorithm) were kept for further analysis, the other positions were rejected. Due to the fact that the 3′-to-5′ end processing is a dynamic process, it is supposed that in some cases several 3′ ends per RNA molecule will be detected. This is why we decided to merge all nucleotide positions within a range of 3 nt with BEDtools’ subcommand merge. The range of mapped positions which belong to one 3′ end is defined as distribution size. All resulting positions are regarded as true differential 3′ ends.
In total 1072 differential 3′ ends could be detected, the majority of them (around 68%) were mapped to a single position (Fig. 5a). By far most of all 3′ ends (82%) were strongly enriched in the pnp mutant strain showing a log2-fold change between + 5 and + 9 (Fig. 5b). These ends represent either termination sites or arise from cleavage by endoribonucleases and are further processed by PNPase in the wild type. The ends were strand-specifically assigned to the different feature classes according to their genomic position with BEDtools intersect. All 3′ ends that did not overlap with any feature were classified in the group UTR. This group may also contain additional sRNAs that have not yet been identified. The most prominent changes could be observed within the coding sequences and the untranslated regions: In both feature classes, around 10 times more RNA 3′ ends were enriched in the pnp mutant compared to the wild type (Fig. 5c). These 885 ends are interpreted as PNPase-degraded. To search for putative motifs, the 15 nt upstream sequence of every pnp enriched differential 3′ end was extracted with BEDtools getfasta. Only non-overlapping 3′ ends were kept for the following analysis to reduce bias during the motif analysis. No binding motif or consensus sequence could be found with MEME Suite (version 5.3.0).
We further detected 189 3′ ends which are enriched in the wild type. One possible explanation for this observation is, that these RNAs are degraded by PNPase from the 3′ end, but at the detected positions stable secondary structures prevent further RNA 3′-to-5′ end degradation in the wild type. For a test of this hypothesis all sequences within windows of 20 nt, 30 nt, 40 nt and 50 nt upstream of the detected wild type enriched 3′ ends were extracted as described above. As a control we selected sequences of the same length but located downstream of the 3′ ends (Supplementary Fig. S8A, Additional file 1). These RNA sequences are supposed to have no effect on the stop of PNPase decay since they are properly degraded. RNAfold (version 2.4.17) was used to compute the minimal folding energy (MFE) of every sequence. Independently of the window size, the distributions of both groups are highly similar and no shift to low MFE values was observable in the upstream sequences (Supplementary Fig. S8B, Additional file 1). But remarkably, the number of unstructured sequences (MFE = 0.0 kcal/mol) was three times higher downstream of the 3′ ends (n = 29, 13.8%) compared to the regions upstream of the 3′ ends (n = 9, 4.8%) (Supplementary Fig. S8C, Additional file 1). We therefore conclude that unstructured sequences may at least enhance degradation by PNPase while highly structured RNA sequences may facilitate a stop of decay. Moreover, other still unknown factors, e.g. binding proteins, are likely to influence degradation by PNPase. Using MEME, no recurring motif that could hint to conserved binding sites was detectable in the upstream sequences.
The group of tRNAs exhibited less differential 3′ ends in the pnp mutant strain, all with rather minor log2-fold changes (median − 1.3). This affects 32 of in total 54 tRNAs. We observed characteristic 3′ ends which show a clear edge in the wild type RNA coverage profile of the tRNAGly and tRNALeu (see Supplementary Fig. S9A, Additional file 1). The other tRNAs coverage profiles of the wild type harbour only minor edges but are stronger sloped than in the pnp mutant (see Supplementary Fig. S9B, Additional file 1). The reason for these differences still remains concealed. Nevertheless, the detection of those 3′ ends depicts the high sensitivity of XPEAP. The influence of the E. coli PNPase on tRNA maturation and degradation has been investigated intensively. PNPase is involved in the repair process of several tRNAs [56], although this enzyme does not affect the tRNA poly(A) tail length [57]. Furthermore, both PNPase and RNase II remove the Rho-independent terminator structure of the leuX tRNA [58]. In chloroplasts of A. thaliana, the PNPase activity is directly linked to the decay of tRNAs [59].
Moreover, our data suggest that PNPase-dependent degradation is not limited to only one site per gene. In many cases, more than one RNA 3′ end per gene was enriched in the mutant (Fig. 5d+e). This affects in particular the 23S rRNAs, but also several mRNAs and ncRNAs. This is not surprising, since often several RNase E cleavage sites per gene could be detected in Rhodobacter sphaeroides, but also in Salmonella enterica [6, 7]. After RNase E cleavage, the RNA fragments become potential new substrates for PNPase.
It should be mentioned that some of the 1072 differential 3′ ends may arise from a stop of the sequencing reaction after about 75 nucleotides in the single-end sequencing. 3′ ends at such a distance from the 5′ end account for only 4.8% of all differential 3′ ends. Moreover, this number also contains real 3′ ends as depicted in Fig. 4: as predicted, a 75 nt RNA derived from IGR 1711-rpsL occurs only in the mutant but not in the wild type.
Intersection analysis
Multiple ribonucleases are involved in RNA processing and turnover. In many organisms, the ribosomal RNA maturation requires an initial endonucleolytic cleavage of the long nascent precursor transcript by RNase III. This is followed by further enzymatic reactions, performed for example by RNase E, J, G and various other enzymes. PNPase can process those RNA species which are newly generated by endonucleases from the 3′ end. Furthermore, both PNPase and RNase E are part of the E. coli degradosome and work together in RNA degradation. In the R. capsulatus degradosom fraction PNPase activity could be detected, but only a small amount [11]. To get more insight into the interplay of RNase E, RNase III and PNPase we analyzed the correlation of 3′ ends from the wild type, the pnp mutant, an RNase III mutant and a strain with reduced RNase E activity. In the RNase III mutant strain, the rnc gene was removed by homologous recombination [24]. Since RNase E is an essential enzyme in Rhodobacter sphaeroides, the rne mutant strain was achieved by replacing the native rne by the gene of the thermosensitive RNase E from E. coli [54]. At 32 °C the enzymatic activity is already reduced and is even more reduced at 42 °C [7].
First, the differential 3′ ends were detected as described in the previous section (see Additional file 3 for the full list of detected RNA 3′ ends). For the comparison of a mutant to the corresponding wild type strain, we only analyzed data that was obtained from the very same sequencing chip. Next, the overlap between the detected 3′ ends was investigated with BEDtools’ subcommand window. Using this function, a window of 1 nt upstream and 1 nt downstream of every differential pnp mutant 3′ end was set prior to the intersection analysis to compensate a potential inaccuracy during 3′ end determination. Two hundred eighty-nine differential 3′ ends could be detected in the rnc mutant, more than half of them (166; 57%) were uniquely overlapping with PNPase dependent 3′ ends (Fig. 6a). Within this group, no 3′ end showed a reduced abundance in the pnp mutant and an increased abundance in the rnc mutant strain (Fig. 6b). By far the highest number of overlapping ends could be found in tRNAs (Fig. 6c). Taken together, 9.7% of all RNase III generated RNA 3′ ends are further processed by PNPase (Fig. 6a + b: quadrant rnc depleted, pnp enriched). Although in E. coli RNase III is involved in pnp mRNA processing [20, 60], we could not observe any effect on the pnp mRNA levels in the RNase III mutant strain in Rhodobacter sphaeroides.
A much greater number of total 3′ ends was identified to be RNase E-dependent: 5445 at 32 °C and 14,873 at 42 °C (Fig. 6d+g). At the permissive temperature, the majority of differential 3′ ends could be assigned to CDS and UTR (Fig. 6e+f). At the non-permissive temperature more than 2.5% of all differential pnp mutant RNA 3′ ends (279 in total, 271 unique pnp mutant ends) overlap with those differential 3′ ends identified in the rne mutant. Plotting the log2-fold changes reveals a specific pattern: 75% of all intersecting 3′ ends are pnp enriched and rne depleted (Fig. 6h), suggesting that these 3′ ends are generated by RNase E cleavage and further removed by PNPase. This group is mainly composed of ends located in CDS and UTRs, but also in several ncRNAs (Fig. 6i). Interestingly only 5.9% of all 3′ ends which are generated by RNase E are further degraded by PNPase. Such a low fraction of around 6% was also observed in S. pyogenes [31]. We hypothesize, that a portion of the RNase E generated 3′ ends can at least partly be trimmed by other 3′-to-5′ exonucleases and thus are not detected in the pnp mutant strain. Also, stable secondary structures at the newly generated 3′ ends may prevent PNPase degradation. At the permissive temperature, the pnp mRNA levels do not differ, whereas at 42 °C the pnp mRNA levels are more than doubled (log2-fold change = 1.75, comparison rne mutant versus wild type) strongly supporting an effect of RNase E on pnp stability.
We further elucidated the statistical significance of the observed overlapping 3′ ends. For every comparison of detected 3′ ends (pnp mutant versus rnc/rne mutant) BEDtools’ subcommand fisher was used to first compute the number of possible 3′ ends, taking into account the genome size and the individual distribution sizes of determined 3′ ends. Second, the number of overlaps and non-overlaps was computed and the Fisher’s exact test applied. For a possible number of 633,092 (pnp versus rnc), 748,506 (pnp versus rne 32 °C) and 667,153 (pnp versus rne 42 °C) RNA 3′ ends the resulting two-tail p-value are close or equal to 0 in each case (0, 7.3 × 10− 241, 2.9 × 10− 239). This strongly suggests a higher number of RNA 3′ ends at the same genomic positions than would be expected, if the given 3′ ends would be randomly distributed within the genome.
The total number of detected 3′ ends as well as the overlapping 3′ ends differ within the strains analyzed in this study. It cannot be excluded, that this observation may be influenced by a different number of uniquely aligned reads: The samples of the RNase E mutant and the corresponding wild type strain previously published [7] showed on average a four times higher number of uniquely aligned reads than the samples sequenced for this study (PNPase mutant, RNase III mutant and wild type strain). On the other hand, although Lécrivain et al. [61] used a somewhat different algorithm and parameters to detect differential 3′ ends in Streptococcus pyogenes, the number of identified 3′ ends (wild type enriched: 183; pnp mutant enriched: 1255) is strikingly similar to our findings in the Gram-negative organism R. sphaeroides (wild type enriched: 189; pnp mutant enriched: 885). This similarity strongly supports the reliability of the data obtained in this study and suggests valid differences in between the analyzed strains.
To illustrate these overlapping 3′ ends we chose two of the previously shown sRNAs. Several differential 3′ ends can be found in the sRNAs CcsR1–4, which are derived from one co-transcript by RNase E cleavage [52]. The read coverage depicts a strong enrichment of CcsR1–4 in the pnp mutant (Fig. 7a), which agrees with the Northern Blot data (Fig. 4a). Several pnp mutant enriched 3′ ends have been detected. Three of them overlap with those RNA 3′ ends which are depleted in the rne mutant in comparison to the wild type, even at 32 °C (Fig. 7a, highlighted in red color). The read coverage of the sRNA IGR_RSP_1711_rpsL reveals, that the first 75 bases of the RNA are enriched in the pnp mutant (Fig. 7b). This shorter fragment was also detected by Northern blot analysis (Fig. 4a). Comparing wild type and the RNase E deficient strain, a 3′ end depleted in the mutant could be found at the very same site. Despite some changes in overall abundance, no RNase III-dependent differential 3′ ends were detectable. For both of the two described sRNAs, a similar processing pattern is proposed. The initial transcript is first processed by RNase E, indicated by the depleted 3′ ends in the rne mutant. Followed by this cleavage, the remaining RNA molecule is further degraded by PNPase thus leading to enriched 3′ ends in the pnp mutant. This observation is in agreement with previous studies in Gram-negative (reviewed in [29, 30]) as well as in Gram-positive bacteria. In Streptococcus pyogenes for example, the endonuclease RNase Y generates new RNA 3′ ends, which are subsequently trimmed by PNPase [31]. PNPase seems not to be a component of the degradosome of the related species R. capsulatus [11], which most likely also accounts for R. sphaeroides. Even though this direct interaction between RNase E and PNPase may be lacking, our data provide strong evidence that stepwise RNA processing is mediated by these enzymes also in Rhodobacter.