Global coordination in adaptation to gene rewiring

Gene rewiring is a common evolutionary phenomenon in nature that may lead to extinction for living organisms. Recent studies on synthetic biology demonstrate that cells can survive genetic rewiring. This survival (adaptation) is often linked to the stochastic expression of rewired genes with random transcriptional changes. However, the probability of adaptation and the underlying common principles are not clear. We performed a systematic survey of an assortment of gene-rewired Escherichia coli strains to address these questions. Three different cell fates, designated good survivors, poor survivors and failures, were observed when the strains starved. Large fluctuations in the expression of the rewired gene were commonly observed with increasing cell size, but these changes were insufficient for adaptation. Cooperative reorganizations in the corresponding operon and genome-wide gene expression largely contributed to the final success. Transcriptome reorganizations that generally showed high-dimensional dynamic changes were restricted within a one-dimensional trajectory for adaptation to gene rewiring, indicating a general path directed toward cellular plasticity for a successful cell fate. This finding of global coordination supports a mechanism of stochastic adaptation and provides novel insights into the design and application of complex genetic or metabolic networks.


INTRODUCTION
Gene rewiring, or an alteration of preexisting genetic connections, is a common evolutionary phenomenon (1)(2)(3). In nature, rewiring is caused by mutations in promoters or by the transposition of genomic sequences to other loci and is considered a means of conferring novel gene regulation and/or essentialities in evolution (4)(5)(6)(7). Such genetic alterations may lead to either survival or extinction for the organism, as the alterations could in turn lead to disturbances in existing regulatory patterns. Nevertheless, increasing experimental studies indicate that cells can survive genetic rewiring with a high probability. For example, newly connected gene networks have been shown to be tolerated by bacteria (8), and successful adaptation has been observed in bacterial and yeast cells undergoing gene rewiring in response to external changes (9)(10)(11).
Successful cell fates were largely explained by stochastic switching of rewired genes (10,(12)(13)(14). That is, the stochasticity of gene expression (15)(16)(17)(18)(19) occasionally produces fit cells in response to unforeseen environments (i.e. stochastic adaptation, Supplementary Figure S1A) (10,13,(20)(21). Stochastic adaptation does not require specific regulation, but it largely relies on cell diversity within a genetically identical population (22,23). Therefore, stochastic adaptation is considered a universal survival strategy in living cells (20,24). Despite previous experimental results and a proposed mechanism for adaptation, how often the stochastic expression of rewired genes results in adaptation and whether stochasticity in rewired genes is solely sufficient for a successful cell fate are not known.
Current studies identified common features of global changes in living organisms, which prompted us to consider a general mechanism for adaptation via transcriptome reorganization under cellular plasticity in rewired cells. Pioneering transcriptome studies commonly observed random patterns in global gene expression, which indicated that cellular plasticity contributed to the adaptation of rewired cells (11,25). In comparison, common patterns of transcriptome reorganization in response to external stress were largely reported in native cells (i.e. wild-type strains). For example, the correlation between growth rate and gene expression was universal regardless of variations in external perturbations in yeast and bacteria (26)(27)(28)(29). Recent studies on cellular networks illustrated the general mechanisms of network dose compensation (30), metabolic restriction (31) and the coordination between proteome and metabolic signals (32). In comparison, whether and how the global reorganization of gene expression contributes to stochastic adaptation is not clear.
Previous studies using several genetic structures (9)(10)(11) were insufficient to reach a general conclusion of the probability of stochastic adaptation or any common patterns in transcriptome reorganization. Therefore, systematic surveys of rewiring cells with varied cell fates were required to determine common properties and investigate organizational principles of living systems (33). This study performed a systematic survey of gene rewiring under a defined stress condition to identify inherent survival approaches of cells and provide novel insights into cellular plasticity. Genetic disruptions were introduced by rewiring the structural genes within the His operon (34)(35)(36) to a monostable synthetic gene circuit at other chromosomal sites in Escherichia coli (37). Because the rewired gene essential for histidine biosynthesis was controlled by a foreign promoter, the stochastic switching of the rewired gene could provide the opportunity for survival from histidine starvation (10). This study first reported three different cell fates that were mediated by stochastic switching of gene expression, which were designated good survivors, poor survivors, and failures. The stochasticity of the rewired genes and reorganization in global gene expression were evaluated to identify similar and distinct features that corresponded to final cell fates. In conclusion, the stochasticity of the rewired gene itself was insufficient to reach a successful cell fate, but the coordinated reorganization of global gene expression within a restrictive dimension was essential. These results indicate a certain level of universal discipline during stochastic adaptation, which is represented by cellular plasticity in rewired cells.

Strains
The OSU11 and OSU12-geneX E. coli strains ( Figure 1A) were constructed as previously described (10,37). Both strains carry a core genetic circuit, P tetA -gfpuv5 and P trcdsred.t4-tetR-PEM7-zeo, located in their genomes at the galK and intC sites, respectively. PtetA-gfpuv5 consists of the promoter PtetA and the mutated gfp gene (gfpuv5) encoding green fluorescent protein (GFP), and P trc -dsred.t4-tetR-PEM7-zeo is comprised of the promoter P trc , dsred.t4 encoding red fluorescent protein (RFP), tetR encoding the repressor protein and the independent expression unit of the zeocin resistance gene, zeo, with its promoter PEM7. The correlation between the co-expressed gfp gene and the downstream rewired genes were previously demonstrated (10,37).

Cell culture and growth rate
Escherichia coli cells were cultured in minimal medium in the presence of 1 mM histidine at 37 • C for several passages to reach a constant growth rate, as previously described (10). Single colonies from the resultant populations were isolated and stored at −80 • C until analysis. All cell cultures were inoculated using the single colony derived glycerol stocks. Precultures were incubated with 100 M isopropyl ␤-D-1-thiogalactopyranoside (IPTG) to induce full expression of the TetR repressor, which was reported by RFP expression. Exponentially growing cells were subsequently transferred to fresh media containing 50 or 100 nM doxycycline hydrochloride (Dox) to induce the expression of the downstream rewired genes, as reported by GFP expression. Following the induction of the rewired genes, the cultures were subjected to both histidine-depleted (-His) and histidine-rich (+His) conditions for additional analysis. Under the conditions of +His and -His (>10 h) conditions, the initial and final cell concentrations were controlled at ∼10 3 and ∼10 7 cells/ml, respectively. Under the conditions of -His (10 min) and -His (2 h) conditions, the cell concentrations of the precultures (before transfer) were maintained at ∼10 7 cells/ml. The cells that were harvested 10 min or 2 h after transfer were counted to confirm their final concentrations, which were often close to their initial concentrations. The growth rate (h −1 ) was evaluated by the speed with which the number of offspring increases during the exponential growth phase, and calculated according to the following commonly used (10,29,38) formula: growth rate (μ t ) = ln(C t /C 0 )/t, where C t , C 0 and t represent the final and initial cell concentrations (cells/ml) and the culture time from initial to final time points (h), respectively.

Histidine depletion
Aliquots of 350 l of exponentially growing cells (10 6 -10 7 cells/ml) were harvested by centrifugation at 8000 rpm for 1 min at 37 • C using spin columns (Ultrafree-MC Centrifugal Filter Units, 0.22 m; Millipore). After discarding the flow-through fraction, the cells were washed with 350 l of the identical medium without histidine. After an additional centrifugation, the cells were suspended in 350 l of fresh, histidine-free medium. The concentration of the cell suspension was determined using flow cytometry, and cells were then inoculated at 10 6 -10 7 cells/ml (slow growth strains) or 10 4 -10 6 cells/ml (other strains) in 5 mL of the identical histidine-free medium.

Flow cytometry
The cell concentration, gene expression (fluorescence intensity), and cell volume of E. coli cells were determined using a flow cytometer (FACSCanto TM II; Becton Dickinson) equipped with a 488-nm argon laser, a 515-545-nm emission filter (GFP), and a 563-589-nm emission filter (RFP). The following PMT voltage settings were applied: forward scatter (FSC), 450; side scatter (SSC), 400; GFP, 500; and RFP, 600. The flow rate for the sample measurements was set to low. Each calculation was performed using 10 000 collected cells. Cell samples were mixed with fluorescent beads (Fluoresbrite YG Microspheres, Calibration Grade 3.00 m; Polysciences) to calculate the cell concentration. Daily detection accuracy was monitored using eight-peak beads (SPHERO TM Rainbow Calibration Particles (eight peaks), 3.0-3.4 m)) for data calibration.

FCM analysis
The data obtained by flow cytometry were converted to the fcs2.0 format and analyzed by custom-designed scripts written in R (39) using the free packages flowCore and flowViz, which are available from the Bioconductor web site (http: //www.bioconductor.org/). Fluorescent beads loaded for the calculation of cell concentration, cell debris and systematic errors resulting from events that occurred at the bottom or top of the instrument's range were eliminated as previously described (10,14). The protein abundance of the rewired genes, which are represented by GFP bias, was calculated by dividing the green fluorescence value (GFP FI) by the red fluorescence value (RFP FI). The relative cell size was evaluated by the FSC value.

Microscopic observation
Cell samples for FCM analysis were applied to a microscopic observation chamber. A volume of 1 l of exponentially growing cells in the presence or absence of histidine was placed on a glass slide and covered with a cover slip. For membrane staining, 1 l of a 200 mg/ml solution of FM 4-64 (Molecular Probes) was added. Fluorescent images were acquired using a fluorescence microscope (Eclipse TE2000-E, Nikon) and recorded using an EMCCD camera (iXon, Andor). Filter sets of 465-495 nm excitation, 515-555 nm emission and 530-550 nm excitation, >575 nm emission were used to measure fluorescence intensity from GFP and FM4-64, respectively. The area of cells was calculated by summing the number of fluorescent pixels using ImageJ software (http://rsbweb.nih.gov/ij/). The relative cell size was calculated by measuring at least 45 cells for each condition.

Microarrays and expression data normalization
Three biological replicates were performed for each condition, resulting in a total of 90 arrays in this study (Supplementary Figure S8). Total RNA was prepared and microarray analysis was performed using an Affymetrix GeneChip R system as described elsewhere (29,40). A highdensity DNA microarray was utilized with the Affymetrix GeneChip system, and data extraction was performed based on the finite hybridization model (41,42) as previously described (29,40). The raw expression data sets were subjected to normalization (43), resulting in a common mean value (logarithmic) in all data sets. To avoid potential noise caused by the small values, normalized expression data with values less than −1.5 were removed, resulting in a total of 3398 genes for subsequent analysis. The averaged expression value for each gene over the three biological replicates was used for computational analyses. Both the normalized expression data sets and the raw CEL files were deposited in the NCBI Gene Expression Omnibus database under the GEO Series accession number (http://www.ncbi. nlm.nih.gov/geo/query/acc.cgi?acc=GSE55719).

Computational analysis
The Bioconductor software package RankProd (44), which is based on the rank product method (45), was employed to identify differential gene expression (DEGs) caused by histidine depletion. Binomial tests were performed to evaluate the significance of the extracted gene groups 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 (39). GSEA (46), which was used to identify the gene groups with significant changes, and PCA (47,48), which was used to classify expression patterns according to gene expression level variance, were performed as previously described (29). PCA converts the high-dimensional data set (i.e. 3398 genes representing 3398 dimensions) into low-dimensional space (e.g. three dimensions). K-means clustering was performed as previously described (29). The gene expression levels (logarithmic scale) obtained under histidine-rich conditions were subtracted from those acquired under histidine-depleted conditions. Consequently, two transcriptional change values (−His (2 h) and −His) and a base value of zero (+His) were acquired for each gene in each strain. −His (10 min) data were eliminated because the growth rate contained substantial noise caused by the short time period being measured. K-means clustering analysis was performed on this data set, which comprised 22 × 3398 values. PCA was performed on three transcriptional change values (−His (10 min), −His (2 h) and −His), which also comprised 22 × 3398 values.

Annotation of gene function and regulation
The entire data set of gene names and categories was obtained from GenoBase, Japan (http://ecoli.aist-nara.ac.jp/ gb6/Download.html). Transcriptional network information was obtained from RegulonDB v8.0 (49) (http://regulondb. ccg.unam.mx). The transcriptional networks comprising more than 15 regulated genes controlled by a regulator were used in the analysis. MultiFun annotation was performed according to the GenProtEC database (50) (http: //genprotec.mbl.edu). The MultiFun classification was applied to annotations comprising >15 genes. The gene categories were in accordance with a previous report (51). GO term annotations for E. coli strain K-12 were obtained from the Gene Ontology database (52,53) (http:// www.geneontology.org). GO terms (biological processes) that comprised fewer than 15 or greater than 1000 genes were excluded from the analysis.

Success and failure in stochastic adaptation
An assortment of rewired E. coli strains, i.e. an OSU12 series ( Figure 1A), was employed to evaluate whether cells with rewired genetic structures were able to survive starvation (i.e. the generality of adaptation to disturbed regulation in response to histidine depletion, which was previously demonstrated using a single strain (10)). Each strain carried a gene from the native His operon rewired to a monostable synthetic gene circuit ( Figure 1A) at another chromosomal location, including a previously reported strain (10). Therefore, these structural genes (hisG, hisD, hisC, hisB, hisA, hisF  and hisI), which are essential for histidine biosynthesis (Figure 1A, I), were no longer subject to their native regulation by the histidine-sensing His operon (34-36) but were controlled by the foreign P tet promoter in response to the chemical inducer doxycycline ( Figure 1A, II). The expression of the rewired genes was reported by the coexpression of green fluorescent protein (GFP). The native strain (i.e. OSU11, Figure 1A, III), which carried the native His operon, was utilized as a control, as previously reported (10,37). Both success and failure were observed among a total of eight strains under histidine starvation. The growth rates of the rewired strains were comparable to the growth rate of the native strain under histidine-rich conditions (Figure 1B, filled), verifying that the genetic reconstruction and the diversity of the rewired genes slightly disturbed cell growth under histidine-rich conditions. However, significant growth decreases, which were accompanied by genespecific effects, were observed under histidine-depleted conditions ( Figure 1B, open). The strains with rewired hisB, hisC or hisF exhibited relatively high growth fitness that was similar to that of the native strain, indicating that these rewired strains successfully achieved adaptive states that were comparable to the native regulation. Conversely, cell growth was undetectable in the hisD-or hisG-rewired strains, providing the first cases of failure in adaptation to gene rewiring. Additionally, largely suppressed growth was observed in the hisA-or hisI-rewired strains, although these strains were sustainable. Such diverse cell fates were commonly observed in the presence of 50-200 nM doxycycline, which determined the basal expression levels of the rewired genes (Supplementary Figure S2). These results provided the first experimental evidence of differentiated cell fates mediated by the stochastic switching of rewired genes under identical conditions.

Cell fate decision was attributed to fluctuation in rewired genes
The fluctuating expression of rewired genes was clearly demonstrated in the population distributions of GFP bias, which represents the specific expression of the rewired gene  (GFP FI/RFP FI, where RFP was constitutively induced; Figure 2A, Supplementary Figure S3). In the good survivors (hisB, hisC and hisF), the GFP bias exhibited a strong preference for induced levels in response to histidine depletion, with high statistical significance ( Figure 2B, asterisks, Supplementary Figure S3), regardless of the primary mean expression levels ( Figure 2B, black). In contrast, such changes were undetected in the poor survivors (hisA and hisI) or in the failures (hisD and hisG), although their primary mean expression levels were relatively higher than that in the good survivors of hisC and hisF.
Notably, the variations in expression within the population were largely related to the final fate, and these results were consistent with the properties of stochastic adaptation. The good survivors exhibited a relatively high degree of cellto-cell variation within the initial populations compared to the other groups ( Figure 2C, filled). The poor survivors exhibited a relatively small initial variation, which tended to increase the variation for adaptation ( Figure 2C, asterisks). A positive correlation was found between the initial variation in expression (i.e. standard deviation) and the final adaptiveness (i.e. growth rate) among the six survivors (Figure 2D, left). In addition, a negative correlation between variation and adaptation in histidine depleted conditions ( Figure 2D, right) was also found. The selection of adaptive cells may have occurred in the populations exhibiting a large initial variation, such that the variation finally became small in the good survivors. Conversely, an increase in variation most likely occurred in the populations of poor survivors to help breed adaptive cells by chance.

Changes in cell size compensated for stochastic adaptation
In addition to the fluctuating gene expression (Figure 2), the cell size, as a typical indicator of adaptivity, was also observed. Cell size influences the fluctuated abundance of gene product, which is known as the system dilution effect (54). Therefore, the means and variations of cell sizes within populations were evaluated. Cell size distributions demonstrated the fluctuation of cell size (Supplementary Figure  S4), as well as what found in gene expression (Figure 2A,  Supplementary Figure S3). The survivors but not the failures tended to increase their cell size in response to histidine depletion ( Figure 3A). The increased cell size, specifically in cell length, was quantitatively verified by microscopic observation (Supplementary Figure S5). This filamentation strategy, which was commonly observed in wild-type cells as a stress response (55)(56)(57), was assumed to compensate the insufficient fluctuation in gene expression and to facilitate future propagation or generate fluctuations in abundance of gene products by system dilution effect. The strains with large variations in cell size, such as the poor-surviving hisIrewired strain ( Figure 3B), may still be in the process of randomly searching for better adaptations.
The slightly increased cell size detected in most strains was correlated with a slightly decreased growth rate under histidine-depleted conditions could be seen ( Figure 3C, Supplementary Figure S5C). We reasoned that most cells in the good survivor populations had completely adapted to the stress, such that their cell sizes were slightly increased with only a small variation. Conversely, both adaptive and maladaptive cells were present in the poor survivor populations, and thus the average cell size was large and exhibited a large variation ( Figure 3B). The phenotypic plasticity at the morphological level might compensate for the insufficient fluctuations in rewired gene expression, particularly in the poor survivors. The results suggested that the mean expression of the rewired gene was not the only determining factor for adaptation, and the variation of expression within populations and/or the changes in cell size facilitated random searching for an adaptive state.

Successful cell fate required the cooperative expression of the His operon
The expression of all rewired genes showed fluctuations to a certain extent (Figures 2 and 3), but the final cell fates (i.e. growth fitness) of these rewired strains in histidine-depleted conditions were dissimilar ( Figure 1B). The presence of cooperative gene regulation was further investigated using transcriptome analyses. Six survivors of the same growing conditions that were used for FCM analysis were examined. The expression patterns between histidine-rich and histidine-depleted conditions largely differentiated in the rewired strains compared with the native strain ( Figure 4A, dot plots). Several differentially expressed genes (DEGs, determined by the rank product method) were identified in the rewired strains; however, only a few overlapped with those in the native strain especially in poor survivors ( Figure 4A, Venn diagrams). Among all DEGs, only seven genes were shared by the survivors (Figure 4B). In addition, gene set enrichment analysis (GSEA), which identified the gene groups and/or networks of significant transcriptional changes in response to histidine depletion, presented large variability (Supplementary Figure S7). These results reflected random patterns of expression at the individual gene or group level in the rewired strains, supported by the previously reports (11,25). Additional analyses of the transcriptome changes in the early responsive phase were performed ( Figure 5A) because of this failure to obtain common regulatory features in the survivors growing in the exponential phase ( Figure 4). Overall, the induced expression of both the rewired gene and the other structural genes of the His operon was detected in the survivors in the absence of histidine and was comparable to the expression by the native strain ( Figure 5B, Supplementary Figure S9). Such cooperated upregulation of all components in the His operon might be necessary for sufficient biosynthesis of histidine. In particular, the structural genes of the His operon were highly upregulated 2 h after histidine depletion in most strains ( Figure 5B, blue), following a 10-min delay for significant changes ( Figure 5B, black). Induced expression of the rewired genes in the corresponding rewired strains, either in early responsive (2 h) or later exponential phases (>10 h), were crucial for the successful cell fate ( Figure 5B, black asterisks), which was consistent to FCM analyses ( Figure 2). Strains that failed to maintain high expression levels at 2 h ( Figure 5B, red asterisks) finally succumbed to maladaptation. Success appeared to largely depend on the original location of the rewired gene within the His operon. The strains carrying genes (hisG and hisD) that were rewired from loci close to the highly conserved 'attenuator' (35)(36)58) failed to maintain the induced expression of these essential rewired genes. The chromosomal locations of genes may limit cellular plasticity, which is consistent with previous findings on the relationship be-tween chromosome structure and plasticity in gene expression (59).

Cooperative global transcriptional reorganization guided a success path
Principal component analysis (PCA) using all data sets was performed to obtain a conceptual abstract of global cooperation for success in stochastic adaptation. Highdimensional transcriptome reorganizations were primarily restricted within a narrow space and composed of any two PCs of the three main PCs ( Figure 6A). In particular, a highly significant correlation was found between the two primary PCs (PC1 and PC2, except for failures), which accounted for ∼54% of the total changes ( Figure 6A, left). This negative correlation illustrated a trade-off in gene expression for a homeostatic transcriptome, and a singledirectional trajectory allowed for survival ( Figure 6A, red dashed line). The large varied distribution of transcriptome changes at 2 h ( Figure 6A, 2 h circles) reflected the manner of stochastic adaptation. However, survivors occurred along this trajectory, and the failures dropped off. Moreover, normalized distances on this PC1/PC2correlated trajectory showed that growth rates correlated to adaptiveness. The distance from the vertical projection of the circles on the trajectory to a defined initial point on the trajectory was calculated and normalized within one unit ( Figure 6B). Data sets of the 10 min bins were excluded because the growth rates that were estimated during this 10 min were unreliable. A positive correlation with high significance was observed between normalized distances and growth rates ( Figure 6C), which indicated that transcriptome reorganizations cooperated with the rewired genes to direct successful adaptation. Growth-correlated genes as determined by K-means clustering clearly demonstrated that the genes positively and negatively loaded on the PC1/PC2 trajectory mostly comprised clusters that positively and negatively correlated to the growth rates, respectively (Supplementary Figure S8A and B), which strongly supports the finding of the PC1/PC2 trajectory as a success path.
Gene enrichment analysis identified gene functions related to mobility/transport and transposition in the positive and negative loadings, respectively ( Figure 6D, Gene category, MultiFun, GO term). K-means clustering analysis confirmed that the enriched gene functions concentrated in nonessential processes (Supplementary Figure S8C). These results were quite different from studies on regular adaptation by native regulations, which often reported that central biological functions, such as translation, were enriched in the growth-correlated gene clusters (27,29). In addition, the enriched transcriptional networks in the positive load ( Figure 6D, TF, positive) partially overlapped with the networks that were identified in the native response to general stresses (29,40) (e.g. rpoD). However, the enrichment of gene regulations failed to identify any regulatory networks with significant changes in the negative loading ( Figure 6D, TF, negative), which is consistent with the enriched gene functions. These results indicate that the stochastic adaptation of these rewired strains was also attributable to cooperative reorganization of the expression of those nonessential genes (e.g. phage/IS in common), which is quite different from the well-known stress response mechanisms that are tightly regulated by the corresponding genes (e.g. rpoS for general starvation stress). This finding raised an intriguing question whether it was the cooperative transcriptional changes in the nonessential genes that potentially assisted the rewired genes of stochastic switching to stay at the proper expression level. In summary, the single-dimension PC1/PC2 trajectory that was highly abstracted from the multidimensional transcriptome reorganization guided success in stochastic adaptation and the contribution of the cooperative expression of both basal-related and nonessential genes. The population analysis based on the reporter gene gfp (Figure 2) and the transcriptome analysis evaluating the genome-wide expression patterns ( Figure 6) drew comparable conclusions, in which successful cell fates that are mediated by directional changes in gene expression can be achieved in the absence of specific regulation.

DISCUSSION
The present study systematically investigated the mechanism of stochastic adaptation that is often proposed to occur via successful fluctuations that correspond to the rewired genes. Previous studies have demonstrated that cells possess the ability to cope with the challenge presented by a particular gene being rewired (9)(10)(11). That is, cells were induced to search for the proper states independent of the specific regulatory mechanism by the stochastic switching of a certain gene (i.e. the rewired genes) essential for cell growth (10,(12)(13)20). In this study, the genes originally controlled by the well-known His operon were systematically subjected to rewiring. By considering growth recovery ( Figure 1B), three different cell fates were identified: highly adaptive (hisB, hisC and hisF), poorly adaptive (hisA and hisI), and maladaptive (hisG and hisD). This survey provides the first experimental evidence regarding the probability (five of seven strains) that a stochastic strategy will succeed under a certain defined condition (i.e. histidine depletion).
It is intriguing that the strains whose rewired genes were regulated in the same way and were involved in an identical pathway had the different cell fates. First, the success of stochastic adaptation likely depends on the order of either regulation or metabolism. Because hisG and hisD are the first genes expressed from the His operon ( Figure 1A), their expression might require highly precise regulation in contrast to other structural genes, which might be loosely controlled under native regulation. In addition to this locus priority in the His operon, the function of genes that are involved in the metabolic pathway may play a role. The initial and final steps of histidine biosynthesis are controlled by hisG and hisD, respectively (Figure 7, red). If the order of gene expression and the order in which a protein plays a role indicate the accuracy of the genetic and metabolic controls (60), then hisD and hisG must respond in a highly efficient or rapid manner. Transcriptome reorganizations of these two rewired strains may fail to reach the PC1/PC2 trajectory ( Figure 6) within the time limit for growth recovery, because random searching for proper expression patterns was timeconsuming. Second, the surrounding metabolic pathways (61,62) may influence the adaptivity. The strains carrying the rewired genes participating in the reactions upstream of the bifurcation for purine metabolism (61,62) were poorly adapted (Figure 7, blue), whereas those located downstream were highly adapted (Figure 7, green), indicating that the key metabolic pathway may contribute to the stochastic switching-mediated fate decision and supporting the conclusion that a global cooperative reorganization of gene expression is essential for the success of stochastic adaptation. Third, the hisF strain showed a less significant change in GFP bias ( Figure 2) and maintained its cell size ( Figure  3), compared with the other two good survivors. We assume that the dosage of hisF gene product that was required for cell growth might be lower than that for other strains, because this strain showed greater growth fitness than the other rewired strains under all of the induced expression conditions (Supplementary Figure S2). The dose-to-fitness relationship of hisF might be different from other genes, most likely because this gene works together with hisH (63) in histidine biosynthesis. Because the partner gene hisH was maintained under native regulation, the rewiring of hisF might be less lethal than the rewiring of other genes that function alone; therefore, perhaps less change is required.
As a first and intriguing finding, the cooperative transcriptional reorganization (Figures 4-6) that occurred globally within a limited space, in addition to the fluctuated expression of the rewired genes and cell size increase, was important for successful stochastic adaptation ( Figure 6C). The poor survivors showed a higher degree of global reorganization than the good survivors (Figures 4 and 6), but they all remained along the PC1/PC2 trajectory. This global coordinative change was assumed to occur partially by chance because it was only the phage-and IS-related genes, but not regulators that were significantly enriched in the negative load of the success path ( Figure 6D). This nonessentiality in enriched gene functions agreed well with the fact that large cell-to-cell variations of GFP bias were observed in the survivors (Figure 2).
A perspective model is proposed based on these results. The PC1/PC2 trajectory represents a success path ( Figure  8, red bold line) in response to environmental perturbations (e.g. histidine depletion). This path guides the success of stochastic adaptation in gene rewiring in the same direction that is regulated by regular adaptation via native regulations. Different cell fates are drawn when the rewired strains carrying damaged regulations encounter unfavorable environments, which is due to the cell-to-cell variation in gene expression and cell size within the initial populations (Figure 8, light gray circle). Cells that luckily switched to the induced expression level of the essential rewired genes (Figure 8, green) in cooperation with synchronized transcriptome reorganization following the success path turned to alive (Figure 8, red circles). Cells of the transcriptome reorganization left this success path (Figure 8, gray arrow) and tended to succumb to death despite stochastic switching of the rewired gene (Figure 8, open circle).
It is unclear whether all 3398 genes had randomly searched for the proper expression levels; nevertheless, the global reorganization occurred directionally along a onedimensional path. This success path is most likely restricted or shaped by the metabolome, because the cell fates mediated by a stochastic strategy were well categorized from a metabolic perspective ( Figure 7). As gene rewiring is thought to be highly likely in natural evolution (1), this onedimensional path may potentially provide a direction for evolutionary changes at the level of the transcriptome. The distance along this defined path may reveal the evolutionary potential for preexisting gene networks to form novel regulatory connections that are adaptive to unforeseen environments. The observed path was admittedly limited to the context of histidine starvation, but this first successful attempt to search for a general path among 3398 potential dimensions strongly indicates that a general design principle underlies cellular plasticity. This finding illustrates the mechanism of stochastic adaptation and offers novel insights into the systematic design and application of complex gene networks.

SUPPLEMENTARY DATA
Supplementary Data are available at NAR Online.