Predicting cellular growth from gene expression signatures

Abstract

Maintaining balanced growth in a changing environment is a fundamental systems-level challenge for cellular physiology, particularly in microorganisms. While the complete set of regulatory and functional pathways supporting growth and cellular proliferation are not yet known, portions of them are well understood. In particular, cellular proliferation is governed by mechanisms that are highly conserved from unicellular to multicellular organisms, and the disruption of these processes in metazoans is a major factor in the development of cancer. In this paper, we develop statistical methodology to identify quantitative aspects of the regulatory mechanisms underlying cellular proliferation in Saccharomyces cerevisiae. We find that the expression levels of a small set of genes can be exploited to predict the instantaneous growth rate of any cellular culture with high accuracy. The predictions obtained in this fashion are robust to changing biological conditions, experimental methods, and technological platforms. The proposed model is also effective in predicting growth rates for the related yeast Saccharomyces bayanus and the highly diverged yeast Schizosaccharomyces pombe, suggesting that the underlying regulatory signature is conserved across a wide range of unicellular evolution. We investigate the biological significance of the gene expression signature that the predictions are based upon from multiple perspectives: by perturbing the regulatory network through the Ras/PKA pathway, observing strong upregulation of growth rate even in the absence of appropriate nutrients, and discovering putative transcription factor binding sites, observing enrichment in growth-correlated genes. More broadly, the proposed methodology enables biological insights about growth at an instantaneous time scale, inaccessible by direct experimental methods. Data and tools enabling others to apply our methods are available at http://function.princeton.edu/growthrate.

Figure 1. Growth phases of a typical cellular culture.

Figure 2. Schematic of a chemostat.

In the chemostat, cells are grown in liquid media . A tank contains a large supply of nutrient containing high concentrations of all growth factors, but a limited concentration (S₀) of the controlling growth factor. The nutrient flows continuously into a growth tube of limited capacity, where the culture grows. The dynamic behavior of the density of the culture (X) and of the concentration of the controlling nutrient (S) in the growth tube is summarized with a system of Michaelis-Menten differential equations. The desired growth rate is attained by manually limiting the concentration of the controlling growth factor in the nutrient provided to the cells.

Figure 3. Representative genes responding to growth rate, specific nutrients, or unsystematically in our chemostat-derived training data.

Our statistical model of growth rate regulation is based on expression data collected from 36 chemostats at six growth rates (0.05 hr⁻¹ through 0.3 hr⁻¹) under six nutrient limitations (Glucose, Nitrogen, Phosphate, Sulfur, Leucine, and Uracil) as described in . By employing the genes responding strongly, consistently, and only to changes in growth rate (and not specific nutrients) as growth-specific genes, we can apply our model to predict relative growth rates in new expression data. Gene expression in our original 36 conditions fell into three main categories as shown here. (A) Genes strongly up- or down-regulated in response to changes in growth rate, independent of limiting nutrient. The most statistically significant members of this set became our growth-specific calibration genes for application of the linear model to other expression data. (B) A subset of conditions highlighting genes with expression levels showing some correlation with growth rate, but with a strong nutrient-specific component. This represents a sizeable portion of the genome (∼25%), with positively growth-correlated genes enriched mainly for ribosomal function and negatively correlated genes enriched for oxidative metabolism. (C) A subset of conditions highlighting genes showing a non-systematic or negligible change in gene expression. Unresponsive genes were enriched for a variety of cellular processes not expected to show a strong relationship with growth, e.g. transcription, DNA metabolism and packaging, secretion, and many others.

Figure 4. Predicted growth rates for S. cerevisiae gene expression datasets.

Our model of the growth rate transcriptional response can be used to predict the growth rate of a cellular culture from gene expression data, robust to the originating biological conditions, growth regime, and experimental platform. Here, we apply the model to three selected data sets to infer relative and absolute growth rates. (A) A brief (<30 s) heat pulse was administered to a steady state chemostat culture immediately before time zero, and gene expression was assayed with an expression time course (see Figure S1 and Table S1). The relative growth rates inferred from this data show an abrupt departure from steady state growth, followed by a return to steady state (including a brief regulatory overshoot). Our predictions monitor these changes in growth rate at an instantaneous time scale (<5 m) inaccessible by standard experimental assays for growth rate. (B) Predicted growth rates for a portion of the environmental stress response data , assaying the response to a 30–37°C heat shock. Our model captures the cessation and resumption of growth induced by the stress, even for a batch culture in which the growth rate is not fixed a priori. (C) A collection of 24 chemostats were run at four growth rates (0.05 hr⁻¹ through 0.2 hr⁻¹) and limited on six different nitrogen sources. Using only expression data from each condition, our model predicts accurate relative growth rates. However, when provided with the known growth rate for a single condition, the model is additionally able to infer absolute growth rates for all other data sets sharing that condition's mRNA reference channel. Note that the actual growth rate is measured empirically and thus deviates slightly from an ideal straight line due to technical variation in the growth equipment.

Figure 5. Assessment of accuracy and outlier detection during growth rate inference.

(A) We performed an out-of-sample cross-validation of our model by randomly sub-sampling 24 of the 36 training expression arrays 1,000 times. We refit our linear model in each random sample, calculated bootstrapped null distributions for all gene parameters, and found sets of the most significant growth-specific genes. These were then used to infer growth rates for the 12 held-out conditions, providing an estimate of the accuracy of the model's growth rate predictions. (B) When predicting the growth rate of a new collection of expression data, our model excludes any calibration gene with an expression level outside the inner fence (1.5 times the inter-quartile range below or above the first or third quartiles). This improves predicted growth rate accuracy while also calling out genes potentially responding to specific non-growth stimuli under some biological condition. For example, in the mild heat shock time course, two of the six outliers are known heat shock genes (HSP26 and HSP78). The other four (YLR327C, MOH1, YBL048W, and TMA10) are uncharacterized genes, suggesting potential roles in the response to heat shock.

Figure 6. Predicted growth rates for S. bayanus and S. pombe expression datasets.

By examining genes orthologous to our ∼70 S. cerevisiae growth-specific calibration genes, we successfully applied our model to predict growth rates in S. bayanus (∼50 orthologous growth-specific genes, ∼20 M years diverged) and S. pombe (∼75 growth-specific genes due to one-to-many mappings, ∼1B years diverged). (A) Predicted growth rates for S. bayanus undergoing the diauxic shift from fermentative to respiratory growth (Table S3). As observed for the S. cerevisiae diauxic shift in , growth pauses as glucose is exhausted and resumes as the yeast begins consuming ethanol. (B) Predicted growth rates for S. bayanus exposed to a 25–37 C heat shock (Table S3). In contrast to Figure 4B, in which S. cerevisiae is observed to recover from a 37 C heat shock, the less-thermotolerant S. bayanus is predicted to halt growth at high temperatures. (C) Predicted growth rates for S. pombe wild-type and rad3Δ time courses, grown normally and exposed to hydroxyurea (HU, an inhibitor of DNA synthesis and thus growth) . Despite the wide evolutionary divergence between S. pombe and our S. cerevisiae training data, predicted growth rates are in substantial agreement with expected biology. Each time course begins with low growth in a synchronized culture. When the synchronization block is released, cells begin growing, wild-type more efficiently than the rad3Δ mutant. Exposure to HU decreases growth over time, and this effect is exacerbated by RAD3 deletion. While the S. cerevisiae RAD3 ortholog MEC1 is essential, knockouts of the MEC1 pathway members SOD1 and LYS7 have been previously observed to induce HU sensitivity .

Figure 7. Differences in growth characteristics of a metabolically cycling culture compared to cells synchronously undergoing the cell division cycle.

We predict periodic bursts of growth during the oxidative phase of the metabolic cycle as described by . Conversely, we observe essentially no variation in growth in cultures synchronously undergoing the cell division cycle, which has been shown to primarily occupy the reductive phase of the metabolic cycle . (A) In cells undergoing metabolic cycling, growth rates are predicted to peak during the oxidative phase of the cycle, where also observes strong upregulation of translational and ribosomal genes. (B) The predicted growth rate for the alpha-factor synchronized cell cycle is essentially constant, after an initial release from the synchronization block. (C) Predicted rates for the alpha-factor synchronized cell cycle also show an initial resumption of growth after alpha-factor block followed by relatively constant growth rate. Taken together, these observations support the claim that growth rate regulation is not specific to any one cell cycle phase. This also agrees with the fact that rapidly growing (and thus fermenting) S. cerevisiae does not partition metabolism into discrete stages, a phenomenon only occurring when reductive metabolism is hindered by nutrient limitation or other stresses.

Figure 8. Perturbations and potential transcriptional regulators of the growth rate response.

(A) Predicted growth rates for gal1Δ cells shifted to glucose, to galactose, and to galactose with a constitutively active RAS2G19V allele. On glucose, rapid growth is induced within ∼40 m; growth on galactose falls to low levels within ∼40 m, as it cannot be metabolized by this mutant. However, when glucose sensing is emulated by artificial activation of the Ras/PKA pathway, the transcriptional regulatory network attempts to induce rapid growth within ∼60–80 m despite the unavailability of appropriate nutrients. This disconnect between actual and perceived cellular state leads to cell death within 4–6 hours and suggests that nutrient sensing (as opposed to metabolic activity or internal cellular state) is responsible for a large portion of the transcriptional growth rate response. (B) Regulatory binding sites enriched in growth up- and down-regulated genes. We clustered the yeast genome by degree of growth rate response, yielding ten clusters with average responses ranging from −12.0 (strongly downregulated with increasing growth rate) to 8.6 (strongly upregulated). The FIRE program predicted 10 regulatory motifs in the upstream flanks and 3′ UTRs of the most up- and down-regulated clusters. These included the known stress-responsive MSN2/4 binding sites in downregulated genes, the ribosomal regulators RAP1 and PUF4 in upregulated genes, and INO4 sites in upregulated genes (possibly corresponding to its role in the stress response and fatty acid biosynthesis . We also identified five additional putative growth regulatory sites for which the binding factor is not yet known.