Overview of RNA-Seq transcriptomics analysis

To make the transcriptomics data comparable with previous proteomics data, we used the identical sampling conditions for transcriptomics as our previous proteomic analysis [26]. As described previously, Synechocystis was grown in BG11 supplemented with 0.20% (v/v) butanol and cell samples of both control and butanol treatment were collected by centrifugation (8,000 × g for 10 min at 4°C) at 24 h, 48 h and 72 h, corresponded to middle-exponential, exponential-stationary transition and stationary phases of the cell growth, respectively.

A total of 79.5-million raw sequencing reads was obtained from the RNA-seq transcriptomics analysis of six samples, with average reads of 13.2-million reads. After a two-step standard data filtering process, first to eliminate reads with low-quality bases (such as multiple N) and reads shorter than 20 bp, and then to eliminate sequence reads mapped to non-coding RNA of Synechocystis, a total of 27.5-million qualified mRNA-based sequence reads were identified (Table 1). The qualified sequence reads have an average genome mapping ratio of 66.4%. To assess the analytical reproducibility between biological replicates, we collected two biological replicates for butanol treated samples at 72 h, and plotted them using the normalized Reads Per Kilobase of Gene per Million Mapped Reads (RPKM) values, the result showed a correlation coefficient around 0.991 (Figure 1), indicating the overall good quality of RNA-sequencing based transcriptomics technology. The sequence reads matched to all 3189 coding genes in Synechocystis genome (data not shown), suggesting excellent sequencing depth and overall transcript coverage.

Potential gene targets related to butanol tolerance

One goal of the integrated OMICS analysis is to achieve a complete coverage of cellular molecules by using complementary techniques targeting different levels of information (i.e., RNA, protein or metabolites) [25]. In this study, our transcriptomic analysis also revealed new cellular responses which were not observed in the previous proteomic analysis [26]: i) Enhanced production of storage compounds: Polyhydroxyalkanoates (PHAs) are common carbon storage compounds that are accumulated during unbalanced growth conditions [43]. Two genes involved in PHA biosynthesis, slr1994 encoding a PHA-specific acetoacetyl-CoA reductase and slr1993 encoding a PHA-specific beta-ketothiolase were found up-regulated by butanol (Table 2). Cyanophycin is a non-ribosomally synthesized peptide, composed of arginine and aspartic acid, accumulates when cells are grown under all unbalanced nutrient conditions except nitrogen starvation, and has been considered as a primary nitrogen reserve compound in cyanobacteria [44]. Transcriptomic analysis showed that the key gene involved in cyanophycin synthesis, slr2002 encoding cyanophycin synthetase was up-regulated by butanol (Table 2). Although PHA and cyanophycin accumulation has been reported for many natural stress conditions, it may worth further investigation how these pathways respond to butanol stress; ii) Enhanced carotenoid biosynthesis: three genes involved in carotenoid biosynthesis were up-regulated: slr1254 encoding phytoene desaturase, slr0940 encoding zeta-carotene desaturase and slr0899 encoding cyanate lyase (Table 2). The results were consistent with the increased photosynthetic activity of Synechocystis upon butanol stress [26]. Carotenoid biosynthesis has been found up-regulated by strong light in Synechococcus PCC7942 [45], and in stress-tolerant mutants of Haematococcus pluvialis[46]. The results provided further evidences that the integrated OMICS approach could be advantageous in revealing global cellular responses.

Metabolomic signatures related to butanol response

GC-MS based metabolomic analysis was used to characterize the time-series metabolic responses of Synechocystis to butanol exposure, with unperturbed cultures as controls. Cell samples used for metabolomic analysis were collected at 24, 48 and 72 h, respectively, the identical time points of sampling for transcriptomic analysis. Three biological replicates were collected for each time point and treatment, thereby yielding a total of 18 samples. The analysis showed that a total of 73 metabolites were chemically identified with great confidence. Although more metabolites were detected in butanol-treated samples (70.4 ± 2.74) than the control samples (64.12 ± 4.01), the number of metabolites identified varied only slightly within control or treatment bins, implying an overall good analytical quality. To further assess the reproducibility of GC-MS metabolomics, we analyzed three technical replicates of one selected sample, and the results showed that most of the metabolites were identified in technical replicates (Date not shown).

The score plot of principal component analysis (PCA) was applied to evaluate the similarities and differences between the 18 metabolomic profiles (Figure 2). The score plot revealed the following features: i) the samples with or without butanol treatment at different time points were distinctly separated, suggesting significant metabolic differences between samples; ii) for the control samples, metabolic changes along the time courses were relatively small, as showed by the clustering patterns of 9 samples; and iii) when compared with controls, significant metabolic changes were observed for butanol-treated samples, especially for samples with 48 and 72 h butanol treatments. One of the butanol-treated biological replicates was slightly different from other two biological replicates at 48 h and 72 h, probably due to the fact the long-term butanol treatment has caused significant cell aggregation [26], which increased the sample heterogeneity. Nevertheless, the overall similar response patterns can still be observed in these replicate samples according to their position in the score plot (Figure 2). Using a cutoff ratio of 1.5 fold between butanol-treated and control samples, and change in at least 5 out of 9 replicate ratios in any time point, we determined 46 metabolites were differentially regulated, in which 35, 41 and 38 metabolites were detected in 24, 48 and 72 h, respectively (Table 5). Pattern analysis showed the 48 metabolites can be divided into at least 6 clusters according to their changes along the treatment time courses. For example, Cluster I included 7 metabolites up-regulated in all three time points, while Cluster II included 7 metabolites up-regulated only in 48 and 72 h after butanol exposure (Table 5).

Metabolomic analysis has identified several metabolites induced by butanol treatment, including 3-phosphoglycerate (3-PG) and glycerol 1-phosphate induced significantly in all three time points, serine induced at 24 and 48 h, and glycine induced at 48 and 72 h after butanol exposure, respectively. The findings were consistent with early studies which showed 3-phosphoglycerate is increasingly withdrawn from the Calvin cycle in S. elongatus PCC 7942 under iron limitation stress [47]. In addition, phosphoglycerate kinase that catalyzes the production of 3-phosphoglycerate from 1,3-bisphosphoglycerate was also found induced in Anabaena sp. PCC7120 under arsenic stress [48]. Moreover, early study has shown that the intracellular levels of organic acids (glyceric, glycolic and glyoxylic acids) and amino acids (glycine and serine) were elevated in salt-treated Anabaena sp. PCC 7120 as compared to those in the control cells [49]. The results suggested that these metabolites could be important part of metabolic responses to both butanol and general environmental stresses.

Previous proteomic study found that a common stress response of Synechocystis under various environmental perturbations, irrespective of amplitude and duration, is the activation of atypical pathways for the acquisition of carbon and nitrogen from urea and arginine, as evidenced by the significant up-regulation of urease that converts urea into CO₂ and ammonia, under most conditions [50]. Our metabolomic analysis showed that urea was induced by butanol, especially at 48 and 72 h. Previous proteomic analysis showed that cyanophycinase, involved in the breakdown of cyanophycin, a storage molecule for excess carbon and nitrogen, into arginine and aspartic acid, was moderately up-regulated under several conditions [50]. Arginine and aspartic acid can be further converted to glutamate and succinate, respectively [51]. Metabolomic analysis showed that aspartic acid was significantly induced at all three time points, and succinic acid and L-glutamic acid were both induced at 48 and 72 h by butanol treatment. These results implied that a similar up-regulated degradation of cyanophycin may also occur under butanol stress.

Integrated transcriptomic and metabolomic analysis has been proposed as a powerful tool to build the relationship between information elements (i.e., genes/transcripts) and functional elements (i.e., metabolites) in cells [25, 52, 53]. In one recent study, integrated transcriptomic and metabolomic approach was used to determine the infection mechanism of Rhodococcus fascians into Arabidopsis thaliana. The transcriptomic analysis showed a significant impact of infection on the primary metabolism of the host, which was then confirmed by subsequent metabolite analysis, for example, invertase transcripts and activities strongly enhanced upon infection, may related to the increase in the hexose:sucrose ratio [54]. In another study to compare the aerobic and anaerobic fermentations of Zymomonas mobilis, researchers found that greater amounts of end products such as acetate, lactate and acetoin were detected under aerobic conditions, while no change in terms of gene expression was found between aerobic and anaerobic conditions in the early exponential growth phase [55], implying the importance to applying integrated technology in uncovering related molecular mechanism. In this study, although only small number of metabolites can be chemically classified in Synechocystis, the metabolomic analysis found increased abundances of aspartic acid and serine, which was consistent with the induction of slr0550 encoding dihydrodipicolinate synthase involved in aspartate pathway, and sll0455 encoding homoserine dehydrogenase involved in serine pathway, respectively (Table 5). In addition, increased abundance of glutamic acid inside the cells was correlated with up-regulation of sll1883 encoding bifunctional ornithine acetyltransferase/N-acetylglutamate synthase protein, sll0461 encoding gamma-glutamyl phosphate reductase, slr0288 encoding glutamate--ammonia ligase, and slr1898 encoding acetylglutamate kinase that are involved in metabolism of glutamate family amino acids (Table 5). Moreover, metabolomic analysis showed the increased abundances of intermediates in the glycolysis pathway, such as glucose-6-P and 3-PG, consistent with the induction of two key genes, slr0752 encoding phosphopyruvate hydratase and sll0745 encoding 6-phosphofructokinase in the glycolysis pathway. Consistent with this result, up-regulation of glycolysis has been reported for various microbes under stress condition [56, 57]. In a recent ¹³C-based flux analysis, a thermophilic ethanol-tolerant Geobacillus thermoglucosidasius M10EXG was found to prefer glycolysis, the pentose phosphate pathway and the TCA cycle for glucose metabolism [58]. On the other hand, for some of differentially regulated metabolites identified, such as urea and cyanophycin, no change was observed for their functionally-related genes in the transcriptomic datasets, which may be due to multiple factors, such as the snapshot nature of the analysis and the different stability of RNA molecules [25]. Nevertheless, the results further demonstrated that transcriptomic and metabolomic technologies could be complementary to each other, allowing better decipherment of cellular responses of Synechocystis under butanol stress.

Validation of potential tolerance targets

Three genes, sll0690, slr0947 and slr1295 which were found induced by butanol exposure at all three time points (i.e. 24, 48 and 72 h) (Table 2), were selected for construction of knockout mutants and for validation of their involvement in butanol resistance. sll0690 encoding a probable transcription regulator, was up-regulated 5–6 folds, slr0947 encoding an OmpR-type DNA-binding response regulator, was up-regulated 2.4-5.5 folds, and slr1295 encoding an iron transport system substrate-binding protein was up-regulated 1.8-5.9 folds by butanol, respectively. Two corresponding proteins of the genes, Slr0947 and Slr1295, were identified in our previous proteomic analysis, in which they were also slightly up-regulated 1.16-1.55 and 1.16-1.57 folds after butanol treatment for 48 h, respectively [26]. After confirmed by PCR and sequencing, the mutants were grown in parallel with wild type Synechocystis in both normal BG11 medium and the BG11 medium supplemented with 0.25% (v/v) butanol. Comparative analysis showed that although there is no visible difference in terms of growth patterns between the wild type and all three mutants in the regular BG11 medium (Figure 3A), gene disruption of sll0690, slr0947 and slr1295 led to increased butanol sensitivity, suggesting they were involved in butanol resistance (Figure 3B). Currently little is known how these genes are involved in butanol tolerance, although early studies have found that the slr1295 gene product, a periplasm-located component of an iron transporter, has a function in protecting photosystem (PS) II [59] and was induced under salt-stress condition [60]; and the slr0947 gene was involved in the regulation of the coupling of phycobilisomes to photosynthetic reaction centers, and reduction of the copy number of slr0947 resulted in decreased efficiency of energy transfer from phycobilisomes to photosystem II relative to photosystem I [61].

Bacterial growth conditions and butanol treatment

Synechocystis sp. PCC 6803 was grown in BG11 medium (pH 7.5) as described previously [26, 27]. Butanol of 0.20% (v/v) was added at the beginning of cultivation. Cells were collected by centrifugation at 8,000 × g for 10 min at 4°C.

RNA preparation and cDNA synthesis

Approximately 10 mg of cell pellets were frozen by liquid nitrogen immediately after centrifugation and cell walls were broken with mechanical cracking at low temperature. Cell pellets were then resuspended in Trizol reagent (Ambion, Austin, TX) and mixed well by vortex. 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 instruction in the miRNeasy Mini Manual (Qiagen, Valencia, CA). The RNA quality and quantity were determined using Agilent 2100 Bioanalyzer (Agilent, Santa Clara, CA) and subjected to cDNA synthesis. The RNA integrity number (RIN) of every RNA sample used for sequencing was more than 8.0. For each sample, 500 ng total RNA were subjected to cDNA synthesis using a NuGEN Ovation® Prokaryotic RNA-Seq System according to manufacturer’s protocol (NuGEN, San Carlos, CA). The resulting double-stranded cDNA was purified using the MinElute Reaction Cleanup Kit (Qiagen, Valencia, CA).

RNA-seq library preparation

The double-stranded cDNA obtained was subjected to library preparation using the Illumina TruSeqTM RNA Sample Preparation Kit (Illumina, San Diego, CA), through a four-step protocol of end repairing, adenylate 3’ ends adding, adapter ligation, and cDNA template enrichment. To determine the quality of libraries, a Qubit® 2.0 Fluorometer and Qubit™ dsDNA HS (Invitrogen, Grand Island, NY) were first used to determine the DNA concentration of the libraries, and then 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 around 300 bp.

Next-generation sequencing

RNA 2 × 100 bp paired-end sequencing was performed using Illumina’s Solexa Genome Analyzer II using the standard protocol. The cDNA library of each sample was loaded to a single lane of an Illumina flow cell. The image deconvolution and calculation of quality value were performed using Goat module (Firecrest v.1.4.0 and Bustard v.1.4.0 programs) of Illumina pipeline v.1.4.

Transcriptomics data analysis

Sequence reads were pre-processed using FASTX Toolkit (v. 0.0.13) to remove low-quality bases, and reads shorter than 20 bp. The qualified sequence reads were then mapped to non-coding RNA (ncRNA) sequences using Bowtie (v. 2.0.0) with default settings. Genome sequences (including ncRNA sequences) and annotation information of Synechocystis sp. PCC 6803 were downloaded from NCBI (Downloaded on April 22, 2012) [27]. Reads that mapped to ncRNA sequences were excluded from further analysis in this study. For paired-end Illumina reads, both pairs were removed if either pair mapped to rRNA. Remaining reads were mapped to the Synechocystis genome using Bowtie (v. 2.0.0) with the default parameters. For gene expression determination, we performed a standard calculation of Reads Per Kilobase of Gene per Million Mapped Reads (RPKM) [63]. We performed comparative transcriptome analysis for all three time points (i.e., 24, 48 and 72 h). To identify the reliable gene targets related to butanol tolerance, only the genes with 1.5-fold induction by butanol at all three time points were regarded as up-regulated genes.

Quantitative real-time RT-PCR analysis

The identical RNA samples used for transcriptomics analysis as described above were used for RT-qPCR analysis. cDNAs were synthesized using RevertAidTM Reverse Transcriptase (Fermentas, Glen Burnie, MD). The qPCR reaction was carried out in 20 μl reactions containing 10 μl of SYBR® Green PCR Master Mix (Applied Biosystems, Foster City, CA), and 2 μl of each PCR primer at 2 mM, employing the StepOnePlus™ Real-Time PCR System (Applied Biosystems, Foster City, CA), under the following condition: 50°C for 2 min and 95°C for 10 min, followed by 40 cycles of 95°C for 15 s and 60°C for 1 min. Quantification of gene expression was determined according to standard process of RT-PCR which used serial dilutions of known concentration of chromosome DNA as template to make a standard curve. The rnpB gene (6803 s01) encoding RNase P subunit B was used as an internal control according to the previous publication [64]. Three technical replicates were performed for each gene. Data analysis was carried out using the StepOnePlus analytical software (Applied Biosystems, Foster City, CA). Data was presented as ratios of the amount of normalized transcript in the treatment to that from the control. The gene ID and their related primer sequences used for real-time RT-PCR analysis were listed in Additional file 1: Table S1.

GC-MS based metabolomics analysis

All chemicals used for metabolome isolation and GC/MS analysis were obtained from Sigma-Aldrich (Taufkirchen, Germany). Cells were collected from control and butanol-treated (0.2% v/v) cultures at 24, 48 and 72 h, respectively. Three biological replicates were established for each sample, and every sample was analyzed three times. For each sample, cells from 5 to 20 mL culture, equivalent to 10⁸ cells mL^-1, were collected by centrifugation at 8000 × g for 10 min at 4°C (Eppendorf 5430R, Hamburg, Germany). The cell pellets were frozen in liquid nitrogen and then stored at −80°C before use. i) Metabolome extraction: cells were re-suspended in 1 mL cold 10:3:1 (v/v/v) methanol: chloroform: H₂O solution (MCW), and frozen in liquid nitrogen and thawed for five times. Supernatants were collected by centrifugation at 14,000 × g for 3 min at 4°C (Eppendorf 5430R, Hamburg, Germany). To normalize variations across samples, an internal standard (IS) solution (100 μg/mL U- 13C-sorbitol, 5 μL) was added to 100 μL supernatant in a 1.5-mL microtube before it was dried by vacuum centrifugation for 2–3 h (4°C). ii) Sample derivatization: derivatization was conducted according to the two-stage technique by Roessner et al. (2001) [65]. The samples were dissolved in 10 μL methoxyamine hydrochloride (40 mg/mL in pyridine) and shaken at 30°C for 90 min, then were added with 90 μL N-methyl- N -(trimethylsilyl) trifluoroacetamide (MSTFA) and incubated at 37°C for 30 min to trimethylsilylate the polar functional groups. The derivate samples were collected by centrifugation at 14,000 × g for 3 min before GC/MS analysis. iii) GC-MS analysis: sample analysis was performed on a GC-MS system-GC 7890 coupled to an MSD 5975 (Agilent Technologies, Inc., Santa Clara, CA, USA) equipped with a HP-5MS capillary column (30 m × 250 mm id). 2 μL derivatized sample was injected in splitless mode at 230°C injector temperature. The GC was operated at constant flow of 1 mL/min helium. The temperature program started isocratic at 45°C for 2 min, followed by temperature ramping of 5°C/ min to a final temperature of 280°C, and then held constant for additional 2 min. The range of mass scan was m/z 38–650. iv) Data processing and statistical analysis: The mass fragmentation spectrum was analyzed using the Automated Mass Spectral Deconvolution and Identification System (AMDIS) [66] to identify the compounds by matching the data with Fiehn Library [67] and the mass spectral library of the National Institute of Standards and Technology (NIST). Peak areas of all identified metabolites were normalized against the internal standard and the acquired relative abundances for each identified metabolite were used for future data analysis.

All metabolomic profile data was first normalized by internal control and cell numbers, and then subjected to Principal Component Analysis using software SIMCA-P 11.5 [68]. Differentially regulated metabolites were determined using a threshold of fold change greater than 1.5 between butanol-treated samples and controls. For each time point, three biological replicates of butanol-treated samples were compared with three biological replicates of control, generating 9 ratios. For each ratio, r > 1.5 was assigned as “ + 1”, r < −1.5 as “-1”, and −1.5 < r < 1.5 as “0”. The sums of the nine ratios for each metabolite at any time point were provided in Table 5.

Construction and analysis of knockout mutants

A fusion PCR based method was employed for the construction of gene knockout fragments [69]. Briefly, for the gene target selected, three sets of primers were designed to amplify a linear DNA fragment containing the chloramphenicol resistance cassette (amplified from a plasmid pACYC184) with two flanking arms of DNA upstream and downstream of the targeted gene. The linear fused PCR amplicon was used directly for transformation into Synechocystis by natural transformation. The chloramphenicol-resistant transformants were obtained and passed several times on fresh BG11 plates supplemented with 10 μg/ml chloramphenicol to achieve complete chromosome segregation. Three genes, sll0690, slr0947 and slr1295 that have been found differentially regulated by butanol exposure, were selected for construction of gene knockout mutants. PCR primers for mutant construction and validation were listed in Additional file 1: Table S1. Full segregation for sll0690 and slr1295 genes was confirmed by PCR. For Δslr0947 mutant, we found that it contained trace amount of original wild-type band in the DNA gels even after more than ten passages, it may worth further investigation whether slr0947 is a lethal gene for the condition. Comparative growth analysis of the wild type 6803 and the mutants were performed in 100-mL flasks each with 10 mL BG11 medium with or without 0.25% (v/v) butanol. Cultivation conditions are the same as described above. Growth analysis was performed in biological triplicates.