Skip to main content
  • Research article
  • Open access
  • Published:

Bacterial transcriptome reorganization in thermal adaptive evolution

Abstract

Background

Evolution optimizes a living system at both the genome and transcriptome levels. Few studies have investigated transcriptome evolution, whereas many studies have explored genome evolution in experimentally evolved cells. However, a comprehensive understanding of evolutionary mechanisms requires knowledge of how evolution shapes gene expression. Here, we analyzed Escherichia coli strains acquired during long-term thermal adaptive evolution.

Results

Evolved and ancestor Escherichia coli cells were exponentially grown under normal and high temperatures for subsequent transcriptome analysis. We found that both the ancestor and evolved cells had comparable magnitudes of transcriptional change in response to heat shock, although the evolutionary progression of their expression patterns during exponential growth was different at either normal or high temperatures. We also identified inverse transcriptional changes that were mediated by differences in growth temperatures and genotypes, as well as negative epistasis between genotype—and heat shock-induced transcriptional changes. Principal component analysis revealed that transcriptome evolution neither approached the responsive state at the high temperature nor returned to the steady state at the regular temperature. We propose that the molecular mechanisms of thermal adaptive evolution involve the optimization of steady-state transcriptomes at high temperatures without disturbing the heat shock response.

Conclusions

Our results suggest that transcriptome evolution works to maintain steady-state gene expression during constrained differentiation at various evolutionary stages, while also maintaining responsiveness to environmental stimuli and transcriptome homeostasis.

Background

Evolutionary experimentation is a powerful tool for the exploration of how living organisms adapt to environmental change and is commonly applied in studies of evolution, typically focusing on changes in genome sequences. The results of these studies have provided experimental evidence to substantiate or revise numerous theories of evolution [1, 2]; the evolutionary mechanisms were generally explained by changes in DNA sequences and/or some particular genes [3–6].

Nevertheless, adaptive evolution also works to optimize cellular gene expression; thus, transcriptome evolution is a relevant area of study. In addition to studies analyzing the transcriptome response to stimuli such as heat shock and nutritional stress [7–11], recent studies report the effects of highly abstract phenomena on global transcriptome parameters. For example, studies showing gene transcriptional changes that were correlated with growth rates, performed in both yeast [12–15] and bacteria [16–18], as well as studies of the trade-off relationship for gene expression responsible for growth fitness and the stress response [18, 19]. Our previous studies identified transcriptome-wide growth-induced transcriptome reorganization [18] and demonstrated that a similar global coordination was also found for stochastic adaptation independent of regulatory mechanisms [20].

Global optimization of gene expression is critical for living cells to achieve novel adaptive states in new environments and must occur not only during transient adaptive responses but also during long-term adaptive evolution [19]. The pioneering transcriptome studies using experimental evolution provided detailed transcriptional changes for genes that are specifically involved in evolution [21, 22]. Studies of parallel evolution in different environments usually found transcriptome diversity and a correlation to genetic background [23–25]. Although a number of studies have reported the general features of transcriptome evolution [6, 23–27], it is unclear whether transcriptome evolution requires alteration of the genomic response to acquire an adaptive state.

As a pilot investigation, we compared transcriptomes of the heat shock response and thermal adaptation. We previously performed experiments in thermal adaptive evolution with a laboratory Escherichia coli strain. We clarified the mechanism of genome evolution, i.e., positive selection vs. neutral evolution [5]; however, how transcriptomes reorganize during the evolution of thermal adaptation remained unknown. The heat shock response is a critical mechanism that living cells use to tolerate high temperatures/thermal stress [28, 29]. In E. coli, this mechanism is largely dependent on feedback regulation from the sigma factor 32 gene (rpoH) and heat shock proteins such as GroEL/ES and DnaK [29, 30]. It is intriguing to know whether Escherichia coli cells alter this responsive machinery to reach adaptive states, or maintain heat shock responsivity using other expression patterns during evolution.

To address this question, we performed a microarray analysis with the representative ancestor (Anc) Escherichia coli cells, and three evolved cell populations (41B, 43B and 45 L), which were acquired at varied evolutionary periods and were well adapted to growth temperatures of 41.2, 43.2 and 44.8 °C, respectively. To distinguish the thermal stress response states and the states of adaptation to high temperatures, we obtained the transcriptomes of both the exponential growth phase at the different growth temperatures (designated as the steady states) and the heat shock response (designated as the responsive states). Multilevel analyses were conducted to investigate the evolution of thermal adaptive transcriptomes. Overall, our results found that long-term evolutionary adaptation to high temperatures was different for the transient response to thermal stress. We summarized these results by proposing a potential mechanism that illustrates transcriptome evolution accompanied by genome evolution.

Methods

Strains

The Escherichia coli strains analyzed in the present study were from a thermal adaptive evolution study performed previously [5]. The genetically engineered Escherichia coli strain DH1ΔleuB::gfpuv5-Km R was used for the evolution experiment, in which the cells were serial transferred with a temperature upshift, from 36.9 to 44.8 °C. The ancestor and evolved strains were designated as Anc, 41B, 43B and 45 L, respectively (Fig. 1a). The genome mutations were determined as previously reported [5]. These four cell populations were subjected to the microarray analysis.

Fig. 1
figure 1

Gene expression of the cells experienced thermal adaptive evolution. a Overview of thermal adaptive evolution. The ancestor (Anc) and evolved strains (41B, 43B, 45 L) used in the present study are indicated, according to our previous report [5]. b Growth rates of the ancestor and evolved strains. Exponential growth rates at both regular (37 °C) and evolution-stimulating (high) temperatures are shown as open and filled bars, respectively. Standard errors of three independent tests are indicated. c Dendrogram of gene expression. Gene expression patterns are clustered as indicated. Mean of the repeated microarray results is used. The three experimental conditions of heat shock, exponential growth at regular and evolutionary temperatures are indicated as hs, r and e, respectively. Black bars highlight the two categories of gene expression, designated as responsive and steady states, respectively

Cell culture and growth

The Escherichia coli cells were cultured in 5 mL of minimal medium M63, supplemented with 2 mM leucine and 25 μg/mL kanamycin, at corresponding temperatures. The culture conditions were exactly the same as previously reported [5]. Cells were counted using flow cytometers, either the FC500 (Beckman) or the FACSCanto II (Becton-Dickson). Cell concentrations were calculated as the ratio of gated particles representing the number of Escherichia coli cells carrying the reporter gene gfp (green fluorescent protein) and beads of known concentrations, as previously described [18, 31]. The growth rates were calculated using the initial and final cell concentrations and the culturing time, as described elsewhere [18, 20, 31]. The initial and final cell densities were maintained at approximately 1 × 104 and 2 × 108 cells/ml, respectively. The culturing times were varied from 14 to 23 h. Cells within the exponential growth phase were collected for microarray analysis.

Heat shock experiments

The condition for the heat shock experiment was a 5-min incubation following a temperature upshift from 36.9 °C to 44.8 °C, as previously described [11]. Exponentially growing cells at approximately 2 × 108 cells/mL were rapidly transferred to an adjacent water bath incubator (Personal-11EX; Taitec) set at 44.8 °C. Following a 5-min incubation at 44.8 °C, the cell culture was immediately poured into a cold phenol-ethanol solution to prepare the samples for microarray analysis. Each heat shock experiment was performed separately to enable precise control of the timing of the heat shock to ensure that the mRNA levels accurately reflected the stress response.

Sample preparation and microarray

Three independent cell cultures were used for the microarray analysis to acquire the mean expression under each culture condition. RNA sample preparation, microarray analysis with the Affymetrix GeneChip system, and data extraction were all based on the finite hybridization (FH) model [32, 33] and were performed as previously described [18, 20]. A single-nucleotide tilling array (high-density DNA microarray) that covers the entire Escherichia coli genome [34] was used. Thus, single-nucleotide substitution bias within ORF (open reading frame) regions could be disregarded in our evaluation of gene expression level. The results of 32 arrays for four bacterial strains in steady and responsive states were used for the analysis, which included triplicates of exponential growth at regular (36.9 °C) and evolution-stimulating (41.2, 43.2 or 44.8 °C) temperatures and duplicates of heat shock conditions.

Data normalization and annotations

Expression data sets from biological replicates were shown as mean values. The raw data sets were subjected to global normalization, resulting in a common median value of zero (logarithmic value) in all data sets, as previously described [18]. The data sets of a total of 4383 genes were used in the analysis. Both the normalized expression data sets and the raw CEL files were deposited into the NCBI Gene Expression Omnibus database under the GEO Series accession number GSE61749. The full datasets of gene names, gene categories [35], and gene regulations (TF) were downloaded from GenoBase, Japan (https://dbarchive.biosciencedbc.jp/jp/genobase-v6/desc.html) and RegulonDB v8.0 [36] (http://regulondb.ccg.unam.mx), as described previously [11, 18].

Computational analysis

Binomial tests were performed to evaluate the statistical significance of extracted gene groups. These statistical analyses were carried out using free software packages available from the Broad Institute (http://www.broadinstitute.org). All statistical tests and computational analyses, except for the gene set enrichment analysis (GSEA), were performed using R [37]. Gene sets enrichment analysis (GSEA) [38] was performed as previously described [11, 18]. The gene categories and gene regulatory links (TF) comprising more than 15 genes were used for the enrichment analysis and bimodal tests. Principal component analysis (PCA) [39, 40], which classifies expression patterns according to gene expression level variance between the culturing conditions, was performed as described previously [18, 20].

Results and discussion

Overview of transcriptome evolution for thermal adaptation

To identify changes in gene expression for thermal adaptation, four Escherichia coli strains were acquired during experimental evolution that we performed previously [5]. We used an ancestor (Anc) population and evolved cell populations 41B, 43B and 45 L, which were adapted to the temperatures 41.2, 43.2 and 44.8 °C, respectively (Fig. 1a). Their growth rates were examined at both the regular (_r) and evolutionary high (_e) temperatures. The results showed that the evolved cells demonstrated highly recovered growth fitness at their corresponding evolutionary temperatures (Fig. 1b), consistent with our previous report [5]. These growing cells were collected for microarray analysis, such that growth rates corresponded to transcriptomes at the exponential growing phases analyzed in this study (Additional file 1: Figure S1). In addition, we acquired the transcriptomes of cells undergoing a heat shock response (_hs, a transient responsive state induced by a temperature increase from 36.9 to 44.8 °C) to identify whether thermal adaptation changed responsivity to a temperature increase (Additional file 1: Figure S1).

An overview dendrogram based on hierarchical clustering analysis clearly shows that global gene expression patterns could be divided into two main clusters, the responsive and the steady states, which corresponded to the transcriptomes during heat shock response and the exponential growing phases, respectively (Fig. 1c). These results suggest that the transcriptome reorganization that occurs during thermal adaptive evolution differs from the reorganization that occurs during a heat shock response. In addition, the transcriptomes of cells in the early and the late evolutionary processes were roughly separated, within the steady-state cluster. This tendency suggests that genetic changes (i.e., genome mutations) played a role in the transcriptional changes during the evolutionary process at 43.2 °C (from 41B to 43B). According to the genome sequence analysis performed previously [5], a comparable number of genome mutations responsible for transcriptional regulation were newly accumulated in each strain (summarized in Additional file 1: Table S1). Overall, these results suggest that global transcriptome optimization might underlie transcriptome evolution.

Equivalent magnitudes of responsive changes

Because the heat shock response is a common reaction to a transient increase in temperature in living cells, we evaluated whether thermal adaptation reduced the magnitude of the general transcriptional changes associated with the heat shock response. Comparison of gene expression during the regular temperature and at heat shock showed that, according to the correlation coefficients, heat shock transcriptional reorganization occurred in the evolved strains, with comparable magnitudes to that found in the ancestor strain (Fig. 2a), which was consistent with the hierarchical clustering results (Fig. 1c). This similarity in the response to heat shock was confirmed by the fact that the expression of the major heat shock genes showed equivalent induction in all strains (Additional file 1: Figure S2A). These results suggest that the evolved strains maintain responsivity and sensitivity to a temperature increase, regardless of long-term evolution under high temperatures.

Fig. 2
figure 2

Changes in gene expression at heat shock responsive states. a. Scatter plots of gene expression in response to heat shock. Individual steady-state expression at regular temperature (_r) is plotted against its heat shock expression (_hs). The strain names and their correlation coefficients are indicated. The direction of evolution is indicated with the arrowed line for reference. b. Heat maps of enriched gene groups showing differentiated expression. The results of GSEA based on two types of annotations, gene category and gene regulation (TF), are shown in the left and right panels, respectively. The statistical significance (FDR q value) is indicated in the heat map. Vivid colors represent high significance in the directions of either up-regulated (yellow) or down-regulated (blue) genes

The enrichment analysis consistently found similar patterns of transcriptional change in most gene categories among all strains, although the genes in the categories of partial information and carrier showed significant differences between Anc and 45 L (Fig. 2b, left). A similarity in the patterns of transcriptional change was also observed for transcriptional regulation. Despite the small number of regulatory networks identified, such as the gene networks regulated by iscR and fhlA, most gene regulation patterns maintained identical directional changes among all strains (Fig. 2b, right). We note that such comparable changes in gene expression were also found in the transcriptional networks controlled by the global regulators harboring nonsynonymous substitutions or deletion mutations (e.g., lrp and oxyR; Additional file 1: Table S1). These results indicate that thermal adaptive evolution did not disturb the mechanism of heat shock response and that the evolved competence at high temperatures might have been acquired using other pathways.

Differentiated transcriptomes at steady states

Because thermal adaptive evolution did not alter the sensitivity of the transcriptome to heat shock, we subsequently analyzed whether and how thermal adaptive evolution disturbed the steady-state transcriptomes at regular and evolution-stimulating temperatures. Comparison of the ancestor and the evolved steady-state transcriptomes showed that longer evolution proceeded larger differentiation (i.e., lower correlation) of the transcriptome from the ancestor at the regular temperature (Fig. 3a). These results suggest that thermal adaptive evolution optimized the global gene expression patterns, not only at evolution-stimulating high temperatures but also at the regular temperature, which did not simulate an evolutionary condition. These results are supported by the fact that the evolved cells all grew faster (improved growth fitness) than the ancestor strain at the regular temperature (Fig. 1b). In addition, comparison of the steady-state transcriptomes at the different growth temperatures showed that transcriptome reorganization for thermal adaptation corresponds with the evolutionary process; i.e., the correlations of gene expression in the steady states at the regular and evolution-stimulating temperatures became weaker in the order of 41B, 43B and 45 L (Fig. 3b), in congruence with the evolutionary process that actually occurred (Fig. 1a).

Fig. 3
figure 3

Changes in gene expression at steady states. a Scatter plots of gene expression at steady states under the regular temperature. The gene expression of the ancestor strain at regular temperature (Anc_r) is plotted against gene expression of the evolved strains at regular temperature (41B_r, 43B_r and 45L_r), which represents a comparison of genotypes (genome mutations). Correlation coefficients are indicated. b Scatter plots of gene expression at steady states under the evolution-stimulating temperatures. The gene expression patterns of the evolved strains at regular (41B_r, 43B_r and 45L_r) and evolutionary-stimulating (41B_e, 43B_e and 45L_e) temperatures are plotted individually, representing a comparison of the growth temperatures. Correlation coefficients are indicated. c Heat maps of enriched gene groups of differentiated expression. The results of GSEA based on gene category are shown. The left and right panels represent comparisons of genotypes and growth temperatures, respectively, corresponding to (a) and (b) panels, respectively. The statistical significance (FDR q value) is indicated in the heat map. Vivid colors represent high significance in the directions of either up-regulated (yellow) or down-regulated (blue) genes. The strains are indicated. The labeled black bars (C1–C4) located at the right of the heat maps indicate the manually categorized clusters showing distinguished changing patterns

Gene set enrichment analysis (GSEA) showed that transcriptome reorganization for thermal adaptation did not randomly occur but followed directional changes (Fig. 3c). The changes in gene expression compared in Fig. 3a represent the transcriptome differences derived from genotypes (i.e., the genomes harboring all the mutations fixed during the thermal adaptive evolution) and similar differences in Fig. 3b are derived from the growth temperatures. Gene function enrichment analysis was performed on the two types of differences, the genotype- and the growth temperature-mediated changes. Roughly, the changes in gene expression mediated by the genotypes were inversely correlated to changes mediated by the growth temperatures (Fig. 3c). The resultant patterns could be manually divided into four clusters (C1–C4). C1 showed significant inverse correlations; C3 and C4 showed gradual directional changes. Such directional changes were also found in regulatory networks (Additional file 1: Figure S3). The inverse transcriptional changes reflected a homeostasis of the transcriptome, which is supported by the finding that the expression of genes responding to a temperature increase (i.e., heat shock genes) returned to normal levels in all evolved cells, regardless of the high temperatures (Additional file 1: Figure S2B).

Thermal adaptive transcriptomes between the regular and the heat shock states

To identify the evolutionary direction of the thermal adaptive transcriptome, principal component analysis (PCA) using all data sets was performed. Analysis based on either the individual data sets (Additional file 1: Figure S1B) or the averaged data sets (Additional file 1: Figure S1A) resulted in the same conclusions (Fig. 4, Additional file 1: Figure S4). Overall, transcriptome reorganization could be roughly interpreted in four main PCs (representing approximately 72 % and 85 % of the total variance in Fig. 4a and Additional file 1: Figure S4A, respectively). Although the magnitude of the changes in gene expression in response to heat shock was approximately consistent in all strains (Fig. 2), the PC1/PC2 space showed that the difference between the transcriptome at evolution-stimulating temperatures and that at the heat shock temperature became smaller when evolution proceeded (Fig. 4a, left). In particular, PC1 represented the adaptivity of the strains in the corresponding conditions because PC1 was correlated with the growth rates (Additional file 1: Figure S5), as previously reported [14–16, 18]. The thermal adaptive expression patterns were located between the steady expression at regular temperature and the responsive expression to heat shock (Fig. 4a, left). We assumed that the adaptive transcriptome was detached from the responsive states and approached the novel steady states. In addition, PC3 and PC4 were likely related to the evolutionary process and to variation in genotypes, respectively (Fig. 4a, right). As shown in Fig. 4a, 43B and 45 L were located at the opposite sides of the Anc–41B zone in PC4. These results suggest that the genetic changes were distinct for thermal adaptation in the later evolutionary phase.

Fig. 4
figure 4

Transcriptome reorganization for thermal adaptation. a Principal component analysis based on individual microarrays. The main four components (P1–P4) are shown. Color variation filled in the circles represents differences among the strains (Anc, 41B, 43B and 45 L), as indicated. The letters hs, r, and e, which are indicated within the circles, represent gene expression conditions, heat shock responsive states, and steady states at regular and evolution-stimulating (high) temperatures, respectively. The weights of each PC are indicated. Pale gray and pink lines illustrate the two zones of the steady expression at the regular temperature and the responsive expression, respectively. b Gene categories and gene regulation significantly contributed to the main PCs. The top 5 % of the genes (439 genes) weighted on each PC were subjected to analysis. Color bars indicate the statistical significance in log-scaled p values obtained using binomial tests with Bonferroni corrections. Asterisks indicate either non-synonymous single-nucleotide substitutions or InDel mutations that occurred during evolution

Functional enrichment analysis was performed for the genes positively and negatively loaded on the four main PCs (top 5 % each, 439 genes in each PC). The enrichment gene analysis (p < 0.001) failed to identify any essential functions (Fig. 4b, Additional file 1: Figure S4). The result suggests that transcriptome adaptation is not simply achieved by a central function of gene activity but is largely mediated by genes with diverse functions. The enrichment of transcriptional networks analysis resulted in largely varied transcription factor regulatory links (TF, p < 0.001) (Fig. 4b). The results showed that only a few regulatory networks were enriched in the most weighted components, PC1 and PC2, indicating that transcriptional networks participated as hubs in the transcriptome, and considerably contributed to thermal adaptation. In comparison, more regulatory links were enriched in the lower weighted components, PC3 and PC4. Differentially expressed regulatory networks were most abundantly identified in PC4, which represented variation in genotypes. This is consistent with the fact that the mutated regulators were largely identified in PC4 and, furthermore, demonstrated that genome mutations (lrp, rpoD and oxyR) did contribute to transcriptional changes at regulatory levels. In addition, of the mutated genes, most were nonsynonymous or InDel (insertion or deletion) mutations, found at considerable frequencies in PC1 and PC4 (Additional file 1: Figure S6). These findings indicate that genome evolution shaped the transcriptomes at PC1, which correlates with growth direction, and at PC4, which correlates with genotype.

Negative epistasis in transcriptome reorganization

Considering our findings of inverse correlations (Fig. 3) and novel adaptive steady states (Fig. 4), we hypothesized that the contributions from genome mutations and temperature increase were negatively correlated. In other words, together, these factors cancel any effect on transcriptome reorganization. Because transcriptional changes caused by heat shock and genome reduction had a negative epistasis relationship [11], we next evaluated whether such a negative epistasis existed in the transcriptomes during thermal adaptive evolution. The transcriptional changes triggered by genome mutations that had accumulated from the ancestor strain (genotype-mediated changes, ΔG) showed low or no correlation to that induced by temperature increase in the ancestor strain (heat stock-induced change, ΔHSA) (Fig. 5a), which was different from a previous study that reported a significant positive correlation between genotype-mediated and heat shock-induced transcriptional changes [11]. As the process of evolution progressed, the weak correlation between the two contributors, genotype and heat shock, seemed to turn from positive to negative. These results indicated that the genome mutations and heat shock contributed to transcriptome reorganization nearly independently, although their interaction also had an effect on evolution.

Fig. 5
figure 5

Epistasis in transcriptome between heat shock and genome mutation. a Scatter plots of the changes in gene expression mediated by heat shock and genotype. The heat shock-induced changes in gene expression (ΔHSA) are plotted against the genotype-mediated changes in gene expression (ΔG) at the regular temperature. Both changes are based on the ancestor (Anc). Correlation coefficients (r) and p-values (p) are indicated. b Negative epistasis. The additive changes in gene expression (ΔHSA + ΔG) are plotted against the simultaneous changes (Δ(G + HSA)), which represent the changes in gene expression between the evolved strains at the heat shock temperature (41B_hs, 43B_hs and 45L_hs) and the ancestor under the regular temperature (Anc_r). The red line indicates linear regression. The values with slopes subtracted from 1 are indicated, representing the magnitudes of the negative epistasis. Gene expression is in log-scale

Despite the weak correlations between genotype and heat shock, a negative epistasis in transcriptome reorganization was clearly detected (Fig. 5b). An approximate 20–30 % reduction in transcriptional changes was estimated to be the result of the simultaneous occurrence of genome mutations and heat shock, compared to their additive occurrence. This result is highly consistent with a previous finding that reported a ~30 % reduction in transcriptional changes [2]. The negative epistasis was most significant for genes that were highly weighted in the main PCs, e.g., a ~50 % reduction in overlapping genes from PCs1–4 (Additional file 1: Table S2). This negative epistasis explains the fact that changes in gene expression were largely inversely related to genotype and temperature (Fig. 3c). Taken together, the changes in gene expression attributed to genetic alteration tended to cancel the changes attributed to the environmental perturbation. Such annulment might enable the cells to maintain the cellular homeostasis of the transcriptome. Although the underlying mechanism was not clearly resolved, our results suggest that adaptive evolution might involve a process of repressing gene expression changes triggered by multiple factors.

A potential path of transcriptome evolution for thermal adaptive

Considering that transition from a heat shock state to a regular steady state generally explained the negative feedback mechanism [28, 29], which was regulated by sigma factor 32 (rpoH) and molecular chaperones, we proposed an evolutionary trajectory for thermal adaption of the transcriptome (Fig. 6). The cell population order corresponds to the evolutionary process (Fig. 1a). The responsivity (Y-axis, the gray circles at the peaks of the distributions) to a transient temperature increase was maintained normally, independent of long-term evolution (Fig. 2). Genome evolution (genome mutations) caused a differentiation of gene expression at the steady states in the order of genome evolution, even under the common regular temperature (X-axis, the hidden gray circles) (Fig. 3). The adaptive transcriptome transitions followed relaxation curves (the distributions in red and ivory), resulting from the feedback mechanism to the heat shock response. Normally, as the result of a heat shock response and negative epistasis (Fig. 5), thermal adaptive evolution of the transcriptome occurs as an extended relaxation process (Fig. 6, the red arrowed line), in which the steady states (Fig. 6, yellow circles) gradually return from the responsive states (Fig. 6, gray circles), bound for a relaxed state different from that of the steady states at the regular temperature (Fig. 6, hidden gray circles). This trajectory represents thermal adaptive evolution with a preference for fewer transcriptional changes between the transient and long-term temperature increase, which was verified by the increasing positive correlations, from 0.48 to 0.72 in 41B to 45 L, between the heat shock-induced changes in gene expression and the growth temperature-mediated changes in gene expression (Additional file 1: Figure S7). The resultant triangles (broken black lines, as an example for 41B) that formed as a result of the changes initiated by heat shock and genome mutations separately and simultaneously represent the cancellation effect on the transcriptome (negative epistasis, Fig. 5b), and this effect possibly became stronger during the course of evolution. We hypothesize that unknown global feedback mechanisms might have played an essential role in transcriptome evolution for thermal adaptation. In summary, thermal adaptive evolution is a process that not only drives the genome from positive to neutral evolution [5] but also optimizes the transcriptome for a balance between responsivity and adaptivity.

Fig. 6
figure 6

Schematic drawing of the proposed evolutionary trajectory. Changes in gene expression in response to a temperature increase are illustrated with transparent distributions, representing the well-known feedback regulatory mechanism in the heat shock response. The responsive states at heat shock are indicated with circles located at the peaks of the distributions indicated with the names of the ancestor and evolved cells. The equivalent heights of the distributions (i.e., the vertical distances from the X-axis to the peaks) indicate that the heat shock response (ΔHS) of both the ancestor and the evolved strains remained in common. The steady states of the evolved cells at the regular temperature are shown with circles located on the axis for genome evolution and are shadowed by the transparent distributions. The order of these strains (i.e., the distances from the starting point of the X-axis) represents the changes in gene expression (ΔG) mediated by the genotypes (i.e., genome mutations). The steady states of the ancestor and evolved cells at the corresponding growth temperatures are highlighted with bright yellow circles. These steady states (circles) are placed at the sides opposite to the regular states along the distributions, suggesting that the thermal adaptive states were distinct to the regular steady states. In addition, they are set on the tails of the distributions, which implies that the steady states under high temperatures were relaxed from the heat shock stress to the corresponding thermal adaptive states. The broken triangle represents an example (41B) of the negative epistasis at transcriptome, as described in the main text. The solid red line, linking up the four steady states at the corresponding evolutionary temperatures, indicates the trajectory of transcriptome evolution along with genome evolution

Conclusions

Thermal adaptation in Escherichia coli cells differentiated steady-state transcriptomes under both normal and evolution-stimulating high temperatures; nevertheless, the heat shock response was maintained with a high sensitivity to thermal stress, as found in the ancestor population. Thermal adaptive evolution directed the transcriptomes to novel steady states different from regular steady states or responsive states. These findings indicate that long-term evolution does not alter existing response machineries but rather adjusts gene expression homeostasis to adaptive steady states.

Availability of supporting data

The raw data is available at the GEO database under:

http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE61749.

References

  1. Elena SF, Lenski RE. Evolution experiments with microorganisms: the dynamics and genetic bases of adaptation. Nat Rev Genet. 2003;4(6):457–69.

    Article  CAS  PubMed  Google Scholar 

  2. Barrick JE, Lenski RE. Genome dynamics during experimental evolution. Nat Rev Genet. 2013;14(12):827–39.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  3. Philippe N, Crozat E, Lenski RE, Schneider D. Evolution of global regulatory networks during a long-term experiment with Escherichia coli. Bioessays. 2007;29(9):846–60.

    Article  PubMed  Google Scholar 

  4. Barrick JE, Yu DS, Yoon SH, Jeong H, Oh TK, Schneider D, et al. Genome evolution and adaptation in a long-term experiment with Escherichia coli. Nature. 2009;461(7268):1243–7.

    Article  CAS  PubMed  Google Scholar 

  5. Kishimoto T, Iijima L, Tatsumi M, Ono N, Oyake A, Hashimoto T, et al. Transition from positive to neutral in mutation fixation along with continuing rising fitness in thermal adaptive evolution. PLoS Genet. 2010;6(10):e1001164.

    Article  PubMed Central  PubMed  Google Scholar 

  6. Wang L, Spira B, Zhou Z, Feng L, Maharjan RP, Li X, et al. Divergence involving global regulatory gene mutations in an Escherichia coli population evolving under phosphate limitation. Genome Biol Evol. 2010;2:478–87.

    Article  PubMed Central  PubMed  Google Scholar 

  7. Patten CL, Kirchhof MG, Schertzberg MR, Morton RA, Schellhorn HE. Microarray analysis of RpoS-mediated gene expression in Escherichia coli K-12. Mol Genet Genomics. 2004;272(5):580–91.

    Article  CAS  PubMed  Google Scholar 

  8. Durfee T, Hansen AM, Zhi H, Blattner FR, Jin DJ. Transcription profiling of the stringent response in Escherichia coli. J Bacteriol. 2008;190(3):1084–96.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  9. Brockmann-Gretza O, Kalinowski J. Global gene expression during stringent response in Corynebacterium glutamicum in presence and absence of the rel gene encoding (p)ppGpp synthase. BMC Genomics. 2006;7:230.

    Article  PubMed Central  PubMed  Google Scholar 

  10. Jozefczuk S, Klie S, Catchpole G, Szymanski J, Cuadros-Inostroza A, Steinhauser D, et al. Metabolomic and transcriptomic stress response of Escherichia coli. Mol Syst Biol. 2010;6:364.

    Article  PubMed Central  PubMed  Google Scholar 

  11. Ying BW, Seno S, Kaneko F, Matsuda H, Yomo T. Multilevel comparative analysis of the contributions of genome reduction and heat shock to the Escherichia coli transcriptome. BMC Genomics. 2013;14:25.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  12. Regenberg B, Grotkjaer T, Winther O, Fausboll A, Akesson M, Bro C, et al. Growth-rate regulated genes have profound impact on interpretation of transcriptome profiling in Saccharomyces cerevisiae. Genome Biol. 2006;7(11):R107.

    Article  PubMed Central  PubMed  Google Scholar 

  13. Castrillo JI, Zeef LA, Hoyle DC, Zhang N, Hayes A, Gardner DC, et al. Growth control of the eukaryote cell: a systems biology study in yeast. J Biol. 2007;6(2):4.

    Article  PubMed Central  PubMed  Google Scholar 

  14. Brauer MJ, Huttenhower C, Airoldi EM, Rosenstein R, Matese JC, Gresham D, et al. Coordination of growth rate, cell cycle, stress response, and metabolic activity in yeast. Mol Biol Cell. 2008;19(1):352–67.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  15. Fazio A, Jewett MC, Daran-Lapujade P, Mustacchi R, Usaite R, Pronk JT, et al. Transcription factor control of growth rate dependent genes in Saccharomyces cerevisiae: a three factor design. BMC Genomics. 2008;9:341.

    Article  PubMed Central  PubMed  Google Scholar 

  16. Nahku R, Valgepea K, Lahtvee PJ, Erm S, Abner K, Adamberg K, et al. Specific growth rate dependent transcriptome profiling of Escherichia coli K12 MG1655 in accelerostat cultures. J Biotechnol. 2010;145(1):60–5.

    Article  CAS  PubMed  Google Scholar 

  17. You C, Okano H, Hui S, Zhang Z, Kim M, Gunderson CW, et al. Coordination of bacterial proteome with metabolism by cyclic AMP signalling. Nature. 2013;500(7462):301–6.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  18. Matsumoto Y, Murakami Y, Tsuru S, Ying BW, Yomo T. Growth rate-coordinated transcriptome reorganization in bacteria. BMC Genomics. 2013;14:808.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  19. Lopez-Maury L, Marguerat S, Bahler J. Tuning gene expression to changing environments: from rapid responses to evolutionary adaptation. Nat Rev Genet. 2008;9(8):583–93.

    Article  CAS  PubMed  Google Scholar 

  20. Murakami Y, Matsumoto Y, Tsuru S, Ying BW, Yomo T. Global coordination in adaptation to gene rewiring. Nucleic Acids Res. 2015;43(2):1304–16.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  21. Cooper TF, Rozen DE, Lenski RE. Parallel changes in gene expression after 20,000 generations of evolution in Escherichia coli. Proc Natl Acad Sci U S A. 2003;100(3):1072–7.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  22. Riehle MM, Bennett AF, Lenski RE, Long AD. Evolutionary changes in heat-inducible gene expression in lines of Escherichia coli adapted to high temperature. Physiol Genomics. 2003;14(1):47–58.

    Article  CAS  PubMed  Google Scholar 

  23. Fong SS, Joyce AR, Palsson BO. Parallel adaptive evolution cultures of Escherichia coli lead to convergent growth phenotypes with different gene expression states. Genome Res. 2005;15(10):1365–72.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  24. Le Gac M, Cooper TF, Cruveiller S, Medigue C, Schneider D. Evolutionary history and genetic parallelism affect correlated responses to evolution. Mol Ecol. 2013;22(12):3292–303.

    Article  PubMed  Google Scholar 

  25. Puentes-Tellez PE, Kovacs AT, Kuipers OP, van Elsas JD. Comparative genomics and transcriptomics analysis of experimentally evolved Escherichia coli MC1000 in complex environments. Environ Microbiol. 2014;16(3):856–70.

    Article  CAS  PubMed  Google Scholar 

  26. Vijayendran C, Barsch A, Friehs K, Niehaus K, Becker A, Flaschel E. Perceiving molecular evolution processes in Escherichia coli by comprehensive metabolite and gene expression profiling. Genome Biol. 2008;9(4):R72.

    Article  PubMed Central  PubMed  Google Scholar 

  27. Pelosi L, Kuhn L, Guetta D, Garin J, Geiselmann J, Lenski RE, et al. Parallel changes in global protein profiles during long-term experimental evolution in Escherichia coli. Genetics. 2006;173(4):1851–69.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  28. Guisbert E, Yura T, Rhodius VA, Gross CA. Convergence of molecular, modeling, and systems approaches for an understanding of the Escherichia coli heat shock response. Microbiol Mol Biol Rev. 2008;72(3):545–54.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  29. Yura T, Nakahigashi K. Regulation of the heat-shock response. Curr Opin Microbiol. 1999;2(2):153–8.

    Article  CAS  PubMed  Google Scholar 

  30. Hartl FU, Martin J. Molecular chaperones in cellular protein folding. Curr Opin Struct Biol. 1995;5(1):92–102.

    Article  CAS  PubMed  Google Scholar 

  31. Ying BW, Tsuru S, Seno S, Matsuda H, Yomo T. Gene expression scaled by distance to the genome replication site. Mol Biosyst. 2014;10(3):375–9.

    Article  CAS  PubMed  Google Scholar 

  32. Ono N, Suzuki S, Furusawa C, Agata T, Kashiwagi A, Shimizu H, et al. An improved physico-chemical model of hybridization on high-density oligonucleotide microarrays. Bioinformatics. 2008;24(10):1278–85.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  33. Ono N, Suzuki S, Furusawa C, Shimizu H, Yomo T. Development of a physical model-based algorithm for the detection of single-nucleotide substitutions by using tiling microarrays. PLoS One. 2013;8(1):e54571.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  34. Suzuki S, Ono N, Furusawa C, Kashiwagi A, Yomo T. Experimental optimization of probe length to increase the sequence specificity of high-density oligonucleotide microarrays. BMC Genomics. 2007;8:373.

    Article  PubMed Central  PubMed  Google Scholar 

  35. Riley M, Abe T, Arnaud MB, Berlyn MK, Blattner FR, Chaudhuri RR, et al. Escherichia coli K-12: a cooperatively developed annotation snapshot--2005. Nucleic Acids Res. 2006;34(1):1–9.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  36. Salgado H, Peralta-Gil M, Gama-Castro S, Santos-Zavaleta A, Muniz-Rascado L, Garcia-Sotelo JS, et al. RegulonDB v8.0: omics data sets, evolutionary conservation, regulatory phrases, cross-validated gold standards and more. Nucleic Acids Res. 2013;41(Database issue):D203–213.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  37. Ihaka R, Gentleman R. R: A language for data analysis and graphics. J Comput Graph Stat. 1996;5(3):299–314.

    Google Scholar 

  38. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102(43):15545–50.

    Article  PubMed Central  CAS  PubMed  Google Scholar 

  39. Venables WN, Ripley BD. Modern Applied Statistics with S. New York: Springer; 2002

  40. Mardia KV, Kent JT, Bibby JM. Multivariate Analysis. London: Academic Press; 1979.

    Google Scholar 

Download references

Acknowledgements

We thank Junko Asada for technical assistance. This work was partially supported by JST KAKENHI No. 23231061 (to TY) and in part by a research grant from the Institute for Fermentation, Osaka, Japan (to BWY).

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Tetsuya Yomo.

Additional information

Competing interests

The authors declare that there is no competing interest.

Authors’ contributions

BWY and TY conceived the research. BWY designed the experiments. KK, SS and BWY performed the experiments. BWY, YM and TY analyzed the data. KK and TK provided experimental materials. SS, NO, CF and TY provided analytical tools. BWY wrote the paper. All authors read and approved the final manuscript.

Additional file

Additional file 1:

Supplementary figures and tables. (PDF 320 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Ying, BW., Matsumoto, Y., Kitahara, K. et al. Bacterial transcriptome reorganization in thermal adaptive evolution. BMC Genomics 16, 802 (2015). https://doi.org/10.1186/s12864-015-1999-x

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-015-1999-x

Keywords