Systematic and functional identification of small non-coding RNAs associated with exogenous biofuel stress in cyanobacterium Synechocystis sp. PCC 6803

Strains, culture, and stress conditions for sRNA samples

Synechocystis sp. PCC 6803 was grown in BG11 medium (pH 7.5) under a light intensity of approximately 50 μmol photons/m² s¹ in an illuminating incubator (HNY-211B Illuminating Shaker, Honour, China) at 130 rpm and 30 °C with a starting cell density of OD₇₃₀ = 0.1 [21, 23–27]. Cell density was measured with a UV-1750 spectrophotometer (Shimadzu, Japan). For growth and stress treatment, 10 mL fresh cells at OD₇₃₀ approximately 0.5 was collected by centrifugation and inoculated into 50 mL BG11 liquid medium in a 250-mL flask. Ethanol 1.5% (v/v), butanol 0.2% (v/v), hexane 0.8% (v/v), and salt 4% (w/v) were added at the beginning of cultivation. The nitrogen starvation condition was established by excluding NaNO₃ from the BG11 medium. Two milliliter of culture was sampled and measured at OD₇₃₀ every 12 h. Finally, a total of 18 culture samples including six conditions at three time points (i.e., 24, 48, and 72 h) were collected for RNA preparation.

RNA preparation and cDNA synthesis

Approximately 10 mg of cell pellets were frozen in liquid nitrogen immediately after centrifugation at 8000×g for 10 min at 4 °C, and cell walls were broken by liquid nitrogen mortar grinding. Cell pellets were re-suspended in TRIzol reagent (Ambion, Austin, TX) and mixed well by vortexing. Total RNA extraction was achieved using a miRNeasy Mini Kit (Qiagen, Valencia, CA). Contaminating DNA in RNA samples was removed with DNase I according to the instructions for the miRNeasy Mini Kit (Qiagen, Valencia, CA). The RNA quality and quantity were determined using an Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA) and subjected to complementary DNA (cDNA) synthesis. The RNA integrity number (RIN) of every RNA sample used for sequencing was more than 7.0. To enrich small RNA for the sRNA-seq analysis, the pool of total RNAs was size-selected, and transcripts smaller than 250 nucleotides (nt) was used to prepare cDNA libraries. For each sample, 500 ng size-fractionated sRNAs were subjected to cDNA synthesis using a NuGEN OvationW Prokaryotic sRNA-Seq System according to the manufacturer’s protocol (NuGEN, San Carlos, CA). The resulting double-stranded cDNA was purified using the MinElute Reaction Cleanup Kit (Qiagen, Valencia, CA).

Library preparation for sRNA and sequencing

The double-stranded cDNA was subjected to library preparation using the Illumina TruSeq™ RNA Sample Preparation Kit (Illumina, San Diego, CA), through a four-step protocol including end repairing, addition of adenylate 3′ ends, adapter ligation, and cDNA template enrichment. The amplification program were: 98 °C 30 s; 98 °C 10 s, 60 °C 30 s, 72 °C 30 s for 15 cycles; 72 °C for 5 min, and hold at 4 °C. To determine the quality of the libraries, a QubitW 2.0 Fluorometer and Qubit™ dsDNA HS (Invitrogen, Grand Island, NY) were first used to determine the DNA concentration of the libraries, and a FlashGel DNA Cassette (Lonza, USA) or Agilent Technologies 2100 Bioanalyzer (Agilent, Santa Clara, CA) was used to determine the product size of the libraries, with good libraries typically 250–300 bp. The product was used directly for cluster generation using an Illumina Solexa Sequencer according to the manufacturer’s instructions.

RNA 2 × 100 bp paired-end sequencing was performed using the standard protocol for the Illumina Genome Analyzer IIx. The cDNA library for each sample was loaded onto a single lane of an Illumina flow cell. The image deconvolution and calculation of quality values were performed using a Goat module (Firecrest v.1.4.0 and Bustard v.1.4.0 programs) with Illumina pipeline v.1.4. Sequenced reads were generated by base calling using the Illumina standard pipeline.

sRNA data analysis

Genome sequence and annotation information for Synechocystis were downloaded from NCBI (ftp://ftp.ncbi.nlm.nih.gov/genomes). The sRNA sequence reads were pre-processed using an NGS QC Toolkit (v. 2.3) to remove low-quality bases and adapter sequences [55]. Reads after QC were aligned to the Synechocystis genome using the Burrows-Wheeler Alignment tool software (v. 0.7.10) [56] with perfect match parameters. As the length of some small RNAs might be shorter than 100 nucleotides, we applied a strategy of 50 cycles of read trimming and re-mapping to detect bacterial sRNAs between 50 and 100 bp [7]. Briefly, we re-extracted the reads that did not match the Synechocystis genome from the aligned SAM files and trimmed one low-quality base from 3′ or 5′ end of these reads (observed by FastQC software). We then re-mapped these trimmed reads to the Synechocystis genome. If reads still did not match, the entire process was repeated until the reads matched the Synechocystis genome. Reads shorter than 50 bp that could align were discarded. Finally, we use Samtools [57] software (v. 0.1.19) to merge all original SAM files from the same sample.

After sRNA pair-end reads were mapped to the Synechocystis genome, Bedtools (v. 2.20.1) [58] was used to calculate read mapping statistics for each sample from BAM files generated with Samtools. The coverage of each nucleotide was calculated by counting the number of reads mapped at corresponding nucleotide positions in the genome. To normalize the sRNA expression level in different samples, we removed nucleotide coverage deep in the rRNA operon and tRNA gene regions, summed all the nucleotide coverage in the remaining genome regions as total mapped bases, and normalized them to 100,000,000 bases. This created a reads base for all samples, which corresponded to an approximately 25× sequence depth for the whole Synechocystis genome.

The sRNAs were identified using a multi-tiered approach. We first searched for enriched regions of sRNA expression and then estimated their 3′ and 5′ positions through a manual correction. For highly transcribed sRNA, we defined s _i as the coverage depth at nucleotide i in the Synechocystis genome and then set a real expression level for a given sRNA at each nucleotide of at least 50× sequencing depth. We looked for the first location in s _i  > 50 representing the start site of an sRNA and then determined whether s _i+1 > 50, until s _i+j  < 50 is the end site of a sRNA. For low-transcribed sRNAs, we used previously reported sRNA gene candidates as references and only retained sRNA with an obvious reads coverage reduction 50 bp upstream or downstream of the adjacent region. To obtain a robust sRNA mapping with a low false positive rate, especially for the induced sRNAs, we retained the sRNAs that were repeatedly observed in at least in three samples across all 18 samples. Finally, manual correction of sRNA boundaries was conducted to identify a point of max rapid coverage decline, considered as the end of the sRNA. We used an R script based on the core code from Kopf [10] to produce a PDF file to display distribution of sRNA reads in the Synechocystis genome (Additional file 2: Figure S1) and four plasmids including pSYSM, pSYSA, pSYSG, and pSYSMX (Additional file 16: Figure S10, Additional file 17: Figure S11, Additional file 18: Figure S12, Additional file 19: Figure S13, respectively). Two adjacent sRNAs located on the same strand shorter than 50 bp with the same expression trend were merged as a single sRNA. The sRNAs shorter than 50 bp were discarded. By determining the 5′ and 3′ ends and inspecting the locations, these potential sRNA genes were classified into trans-encoded sRNAs, cis-antisense sRNAs and UTRs of mRNAs based on their location and were annotated as ncRNA (nc), asRNA (as) and UTR (U), respectively.

To overcome challenge that potential identified sRNAs due to reads matched to multiple locations [59], repetitive regions in the Synechocystis genome were searching by BLAST software. Fragments with identity >80%, E value <1e ⁻⁵ and length >50 bp were considered a repeat region. The IS was predict using Isfinder software [60]. The sRNA-seq data was compared with the Nr database and Rfam database with a 1e ⁻¹⁰ cut-off. Candidate ORFs and RBSs were predicted by Glimmer [61] and RBSfinder [62], respectively. Rho-independent terminators in Synechocystis were searched for using TransTermHP [63] with standard settings. RNA secondary structure and ∆ G analyses of sRNAs were performed using RNAfold software [45].

Quantification and statistical analysis of sRNAs

For all comparisons, the aligned total base counts were normalized to obtain relative levels of expression. We defined s _i as the coverage depth at nucleotide i in an sRNA, summed coverage depth at each nucleotide for an sRNA and divided by its length. The sRNA expression level was calculated as: \(\mathop \sum \nolimits_{i}^{i + j} {{{\text{s}}_{i} } \mathord{\left/ {\vphantom {{{\text{s}}_{i} } j}} \right. \kern-0pt} j}\). Differentially expressed sRNAs were identified using the R package from DESeq software [46] with identified sRNAs as input. For comparison, the resulting p values were adjusted using Benjamini and Hochberg’s approach for controlling the false discovery rate. The sRNAs with fold change >1.5 and adjusted p values <0.05 were assigned as significantly differentially expressed. Correspondence analysis was conducted using R software. The annular heatmap representation was conducted using the OmicCircos package [64] in R software.

Construction of sRNA regulatory networks

For construction of sRNA regulatory networks, each trans-encoded sRNA target prediction was first performed with CopraRNA software [65] using sRNA homologs from six strains of Synechocystis 6803 genomes (NC_000911, NC_017038, NC_017039, NC_017052, NC_017277, and NC_020286) as inputs. The top 100 predictions obtained from CopraRNA with a free-energy cut-off of −10 kcal/mol were retained to remove potential false positive targets [34]. Later, a Pearson correlation analysis with a set of paired RNA samples in previous studies was used to further improve the accuracy of target prediction. Only correlation <−0.4 or >0.4 and p values <0.05 by Fisher’s exact test in two-sided analysis were kept to reduce false positive target prediction. The approach is a modification of previous reports where only anti-correlated relationships were retained for sRNA regulatory network construction [66] as positive regulation of sRNA also found in Synechocystis [67]. To assess sRNAs-regulated metabolic pathways, we performed functional enrichment analyses for each trans-encoded sRNA (see “Functional enrichment analysis” section). Only pathways containing at least two target genes with hypergeometric test p values <0.05 were considered enriched metabolic pathways potentially regulated by sRNAs (Additional file 20: Table S7). Finally, we generated an sRNA regulatory network after assembling all significantly enriched target results. A display of the sRNA regulatory network was conducted with Cytoscape software [68].

Functional enrichment analysis

N is the total number of genes with KEGG pathway annotation information, M is the number of genes with a given KEGG pathway annotation, n is the number of sRNA target genes with all KEGG pathway annotation information and m is the number of the sRNA target genes with a given KEGG pathway annotation. KEGG pathways with p values less than 0.05 were considered enriched response pathways.

qRT-PCR analysis and two-step RT-PCR analysis

The RNA samples used in sRNA sequencing and qRT-PCR were prepared from identical cultures, and qRT-PCR analysis was performed as previously described [21]. Quantification of sRNA expression was determined according to a standard process of qRT-PCR that used serial dilutions of known concentrations of chromosomal DNA as a template to construct a standard curve. A total of 20 sRNAs were selected for verification and 16S rRNA was used as an internal control. Three technical replicates and three biological replicates were analyzed for each sRNA. The data analysis was carried out using the StepOnePlus analytical software (Applied Biosystems, Foster City, CA). Briefly, the amount of relative gene transcript was normalized by that of 16S rRNA in each sample, and the data presented were ratios of the amount of normalized transcripts in the treatment between the stress-treated and normal conditions.

In two-step RT-PCR analysis, a specific set of primer pairs located in the inner boundary of sRNA was designed for each sRNA (Additional file 9: Table S4). As Additional file 21: Figure S14 shows, two opposite RT primers that sit inside the sRNA region were separately added to approximately 200 ng of total RNA from first-strand cDNA synthesis at 40 °C for 1 h using the RevertAid™ First Strand Synthesis Kit (Fermentas, USA). Both nested primers of the sRNAs were added to 1 μL of the two separate first-strand cDNA products for amplification using Dream Taq Green PCR Master Mix (Thermo Scientific). The PCR cycling was as follows: 95 °C 2 min, then 95 °C 30 s, 55 °C 30 s or 57 °C 30 s, 72 °C 15 s for 35 cycles, followed by a 5-min final extension at 72 °C. The PCR products were separated on 3.0% agarose gels.

Construction and analysis of sRNA mutant strains

Gene expression vector pJA2, kindly provided by Prof. Paul Hudson (KTH Royal Institute of Technology of Sweden), was used to overexpress the sRNAs. All sRNAs were cloned under the control of the psbA2 promoter. Briefly, the pJA2 backbone was amplified by PCR, treated with DpnI and digested with BamHI and XbaI to create cohesive ends. The sRNA sequence was PCR-amplified using primers pJA2-sRNA-F and pJA2-sRNA-R and cloned into the BamHI/XbaI sites of pJA2, resulting in the recombinant plasmid pJA2-sRNA (details in Additional file 22: Figure S15). The plasmid was introduced into the WT by electro-transformation as previously described [69]. Positive clones grew on BG-11 agar plates with 10 μg/mL kanamycin and were confirmed by colony PCR analysis. The sRNA overexpression strains were designated WT/pJA2-sRNA.

Phenotypic analysis

Synechocystis and sRNA mutant strains were grown under the same culture condition with sRNA sampling starting at a cell density of OD₇₃₀ = 0.1. For biofuel treatment, 2.0–2.5% (v/v) ethanol, 0.20–0.25% (v/v) butanol, or 0.8–1.0% (v/v) hexane was added at the beginning of the cultivation. For salt treatment, BG11 with 5% NaCl (w/v) was prepared and sterilized. Next, 3–4% (w/v) NaCl of BG11 was prepared by adjusting the ratio between normal and 5% NaCl (w/v) BG11 at the beginning of cultivation. For nitrogen starvation treatment, fresh cells at the same (logarithmic) phase were collected by centrifugation at 1500×g at 4 °C and washed twice with nitrogen depletion BG11 medium. Cells were inoculated into nitrogen-depleted BG11 liquid medium in flasks. Cell density was measured on a UV-1750 spectrophotometer (Shimadzu, Japan) at OD₇₃₀. Growth experiments were repeated at least five times to confirm growth patterns.

Flow cytometric analysis

Flow cytometric analysis was performed using a fluorescence-activated cell sorting (FACS) Calibur cytometer (Becton–Dickinson) equipped with a 488-nm blue laser with the following settings: forward scatter (FSC), E00 log; side scatter, 400 V. Control and stress-treated cells were harvested at 48 h, washed twice with phosphate buffer (pH 7.2), and re-suspended in the same buffer to a final OD₅₈₀ of 0.3 (approximately 1.5 × 10⁷ cells/mL¹). A total of 5 × 10⁴ cells were used for each analysis according to the published method [70]. Chlorophyll fluorescence was detected by the FL3 channel with a 670/LP filter. The data analysis was conducted using CellQuest software, version 3.1 (Becton–Dickinson).