The cellular growth rate controls overall mRNA turnover, and modulates either transcription or degradation rates of particular gene regulons

Published online 29 December 2015 Nucleic Acids Research, 2016, Vol. 44, No. 8 3643–3658 doi: 10.1093/nar/gkv1512 The cellular growth rate controls overall mRNA turnover, and modulates either transcription or degradation rates of particular gene regulons Jos´ eGarc ´ ıa-Mart´ ınez1,2,†, Lidia Delgado-Ramos3,4,†, Guillermo Ayala5, Vicent Pelechano6, Daniel A. Medina2,7, Fany Carrasco2,7,Ram´ on Gonz´ alez8, Eduardo Andr´ es-Le´ on3, Lars Steinmetz6,9,10, Jonas Warringer11, Sebasti´ an Ch´ avez3,4,* and Jos´ eE.P ´ erez-Ort´ ın2,7,* 1Departamento de Gen´ etica, Facultad de Ciencias Biol´ ogicas, Universitat de Val` encia. C/Dr. Moliner 50. E46100, Burjassot, Spain, 2ERI Biotecmed, Facultad de Ciencias Biol´ ogicas, Universitat de Valencia. C/Dr. Moliner 50. E46100, Burjassot, Spain, 3Instituto de Biomedicina de Sevilla (IBiS), Hospital Virgen del Roc´ ıo-CSIC-Universidad de Sevilla, C/Antonio Maura Montaner, E41013 Sevilla, 4Departamento de Gen´ etica, Universidad de Sevilla, Avenida de la Reina Mercedes s/n, E41012, Spain, 5Departamento de Estad´ ıstica e Investigaci´ on Operativa, Facultad de Matem´ aticas, Universitat de Val` encia. C/Dr. Moliner 50. E46100, Burjassot, Spain, 6European Molecular Biology Laboratory (EMBL), Genome Biology Unit, Meyerhofstrasse 1, 69117 Heidelberg, Germany, 7Departamento de Bioqu´ ımica y Biolog´ ıa Molecular, Facultad de Ciencias Biol´ ogicas, Universitat de Valencia. C/Dr. Moliner 50. E46100, Burjassot, Spain, 8Instituto de Ciencias de la Vid y del Vino (CSIC, Universidad de La Rioja, Gobierno de La Rioja), Finca La Grajera LO-20 Salida 13, Autov´ ıa del Camino de Santiago, E26007 Logro˜ no, Spain, 9Stanford University School of Medicine, Department of Genetics, Stanford, CA 94305, USA, 10Stanford Genome Technology Center, 3165 Porter Dr. Palo Alto, CA 94305, USA and 11Department of Chemistry and Molecular Biology, University of Gothenburg, Medicinaregatan 9 c, 40530 G¨ oteborg, Sweden Received November 04, 2015; Revised November 30, 2015; Accepted December 16, 2015 ABSTRACT We analyzed 80 different genomic experiments, and found a positive correlation between both RNA polymerase II transcription and mRNA degradation with growth rates in yeast. Thus, in spite of the marked variation in mRNA turnover, the total mRNA concentration remained approximately constant. Some genes, however, regulated their mRNA concentration by uncoupling mRNA stability from the transcription rate. Ribosome-related genes modulated their transcription rates to increase mRNA levels under fast growth. In contrast, mitochondria-related and stressinduced genes lowered mRNA levels by reducing mRNA stability or the transcription rate, respectively. We also detected these regulations within the heterogeneity of a wild-type cell population growing in optimal conditions. The transcriptomic analysis of sorted microcolonies confirmed that the growth rate dictates alternative expression programs by modulating transcription and mRNA decay. The regulation of overall mRNA turnover keeps a constant ratio between mRNA decay and the dilution of [mRNA] caused by cellular growth. This regulation minimizes the indiscriminate transmission of mRNAs from mother to daughter cells, and favors the response capacity of the latter to physiological signals and environmental changes. We also conclude that, by uncoupling mRNA synthesis from decay, cells control the mRNA abundance of those gene regulons that characterize fast and slow growth. INTRODUCTION Cells must adapt to changing environmental conditions to maintain fitness and to compete with other genotypes during the natural selection process. As fitness is the net growth of a genotype over time, one of the most important variables for single cell organisms during the course of adaptation is the growth rate (GR) of a population composed only of one genotype in relation to the growth rate of other isogenic populations. It is often assumed that higher population growth rates in microorganisms require higher protein synthesis rates (1–5). This is because proteins constitute a *To whom correspondence should be addressed. Tel: +34 963543467; Fax +34 963544635; Email: [email protected] Correspondence may also be addressed to Sebasti´ an Ch´ avez. Tel: +34 955923127; Fax +34 955 923 101; Email: [email protected] †These authors contributed equally to this work as first authors. C The Author(s) 2015. Published by Oxford University Press on behalf of Nucleic Acids Research. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact [email protected] 3644 Nucleic Acids Research, 2016, Vol. 44, No. 8 large fraction of dry mass in both prokaryotes (6) and eukaryotes (e.g. yeast (7)) and, therefore, their synthesis is the most energetically demanding process (8). This process is not exclusive of unicellular organisms and probably fastproliferating tumor cells are also governed by these rules (9,10). The translation machinery includes the most abundant noncoding RNAs: rRNA and tRNAs. Thus the eukaryotic RNA polymerases (RNA pol) devoted to the synthesis of rRNA and tRNA (RNA pol I and III) must increase their transcription rates (TRs, see the explanation about the acronyms used in M&M) in parallel to the GR (11). RNA pol II, however, transcribes a much larger set of mRNAs subjected to multiple regulatory influences, many of which lower in concentration with the GR (12). So although a large set of mRNAs, which are highly transcribed, is also devoted to ribosome biosynthesis and translation factors (8), the relationship between the GR and RNA pol II TR is not obviously predictable. Unlike rRNAs and tRNAs, which are stable molecules, mRNAs have a shorter half-life, and both their synthesis and decay significantly contribute to regulate their abundance (13). The RNA pol II transcription rate, especially for single cell organisms, is highly variable and a proxy of their physiology and metabolism (3). The influence of the cell cycle (14) and cell size (15)on the transcription rate has been investigated in the model yeast Saccharomyces cerevisiae. Cell cycle length and cell size are related to the GR (15). Transcriptome dependence on the yeast GR has been thoroughly investigated in D. Botstein’s laboratory (16–19). These authors have found that the mRNA amount/concentration (RA) of about one third of the genes varies in a GR-dependent manner. In fact based on an RA signature of 72 genes, it is possible to predict the instantaneous GR with reasonable accuracy (17). The [mRNA] of a given gene is controlled by the action of two opposite forces: synthesis by RNA pol II and degradation by several exonuclease pathways (20). In the exponential phase, the GR is a distinctive feature of any yeast strain and environmental condition. In this phase, most mRNAs are in a steady state (21), in which synthesis and degradation rates (DR) are identical. Synthesis is a zero-order reaction, whereas degradation is a first-order reaction that depends on a kinetic degradation constant (kd) and the mRNA concentration. Therefore, TR =DR =kd·RA where kdis a particular feature of each mRNA species and is inversely related with the mRNA half-life (HL): kd=ln2/HL. These simple equations allow any of the three kinetic parameters (RA, TR, or kd) to be derived from the other two (22). Given the appearance of several genome-wide techniques to measure the synthesis and degradation rates and the mRNA concentrations, in several model organisms, it is possible to quantify the respective contributions of transcription and degradation to the actual concentration for each individual mRNA (13,23), and also for the sum of them all. To do this, we used multiple data sets of different mutants and altered the key pathways in gene regulation and under different growth conditions obtained by ourselves and others to test the relationship of the three kinetic parameters with the growth rate at global and single gene levels. We found that TR and DR generally tended to increase when growing faster, and in such a way that the total mRNA abundance for most genes remained independent of the GR in a yeast cell. However, some genes, which corresponded to particular protein functions, showed increases and/or decreases exclusively on one side of the equilibrium, which provoked changes in the net abundance of specific mRNAs. With an increasing GR, the abundance of the translation-related mRNAs and the stress-induced mRNAs increased and decreased, respectively, due to transcription in both cases. The abundance of the mitochondriarelated mRNAs also decreased with an increasing growth rate, which led to a higher proportion of fermentation over respiration at a higher GR, even in the same high glucose medium. However, this was caused by increased mRNA destabilization, and not by reduced transcription. We confirmed the latter by a totally different experimental approach. We took advantage of the intrinsic heterogeneity of the GR within wild-type yeast populations under optimal environmental conditions (24)tosetupanew in-house developed procedure to allow the separation of yeast microcolonies according to their GR from a population of genetically identical cells. When comparing mRNA expression levels, we found that slowand fast-proliferating subpopulations exhibited differential transcriptomic signatures. This suggests the existence of distinct expression programs that are closely linked to the GR and are independent of the environment. The mitochondriaand respirationrelated genes were overexpressed in the cells from slower growing microcolonies. The PUF3 gene, which encodes an RNA-binding protein (RBP) that regulates the stability and translation of many mRNAs related with mitochondrial functions, was also overexpressed at a low growth rate in both the mutant collection and slow-growing microcolonies. Slower growing microcolonies were specifically enriched in long mRNA isoforms with Puf3-binding sites at their 3UTRs. We propose that the growth rate influences gene expression by acting on both sides of the [mRNA] equilibrium: synthesis and degradation. The importance of these two opposite processes is distinct for different classes of gene functions. In particular for respiration-related mRNAs, Puf3 exerts a stabilizing influence and promotes respiration at a low GR in a manner that is independent of external glucose concentrations. MATERIALS AND METHODS Meta-analysis of experimental data sets In the first meta-analysis we used several S. cerevisiae transcriptomic data sets (published and unpublished) obtained by us. These data included 42 data sets for nascent transcription rates and 21 data sets for mRNA amounts, both obtained by the GRO protocol (25). All these data sets were obtained from exponentially growing yeast populations. Differences in GRs were due to distinct culture conditions (carbon source, temperature), and also to the single deletions of individual genes. In 34 cases, the TR of a reference sample was simultaneously determined. The nascent TR and RA values for individual genes were measured as ratios with regard to a relevant reference sample: a wildtype (mostly BY4741, see Supplementary Table S1) strain Nucleic Acids Research, 2016, Vol. 44, No. 8 3645 at 28◦C in YPD (yeast extract 1%, peptone 1%, glucose 2%) medium. In this way, data can be compared across experiments. Individual values were summed to obtain the global estimates of all the genes. For the total TR (RNA pol I + II + III), the total radioactivity incorporated during a run-on experiment was used. A recently published set of 44 yeast mutants was also used, in which the authors employed uniform culture conditions: exponential phase cells grown at 30◦CinYPD(26) for which the mature synthesis rate (TR), RA and kdwere determined according to the cDTA protocol. Since different mutants have distinct cell volumes, the TR and RA values were corrected by dividing each TR and RA gene value by each mutant/wt fold of cell volume compared to its wildtype strain. The data for the different yeast mutants and their wild type were obtained from various sources (27,28). Nascent TR and mature mRNA synthesis rate represent different aspects of the same phenomenon (22). In this paper, for the sake of simplicity, we used the acronym TR to refer to the molecular process of RNA synthesis by any RNA polymerase. Most of the time, we refer only to the RNA pol II TR, otherwise we indicate what we refer to. In GRO experiments TR is calculated as nascent transcription rate but in the cDTA method (26) a mature mRNA synthesis rate is determined by in vivo RNA labeling. The term synthesis rate is used when talking about chemical equilibrium (as in the equation described in the Introduction section); where it refers to the rate of change of the mature mRNA concentration in the cytoplasm (22). The synthesis rate can be inferred from nascent TR data by assuming that a fixed percentage of nascent mRNA molecules reach the cytoplasm. Growth rate estimations For all our experimental conditions, we obtained the GR by growing 50 ml of yeast cultures in 250-ml flasks with shaking (190 rpm) at the desired temperature. Aliquots were taken every 30 min in the exponential phase and their OD600 (from 0.05 to 0.7) were measured. The GR (in h−1)inthe exponential phase was calculated from growth curves (Supplementary Table S1A). Similarly for the 44 yeast mutants used from Sun et al. (26), the GRs of 38 of these strains were calculated from the growth curves obtained by microcultivation in a Bioscreen C reader, as described elsewhere (29) (Supplementary Table S1B). The exponential phase GRs were extracted from growth curves and averaged over replicates (n =2) (Supplementary Table S1A), as previously described (30). Experiment set and data normalization We employed data sets from many different experiments obtained from various platforms. The expression values acquired from these experiments were normalized by means of quantile normalization (31), and were implemented into the function NormalizeBetweenArrays of the Limma R package (32). Correlation analysis between parameters and the growth rate: global values The global tendencies between different parameters and the GR for our data sets and for those of Sun et al. (26)werecalculated by plotting the total RNA pol II transcription, and also for the mRNA concentration and the mRNA degradation constant (kd). They were all represented as being relative to their wild-type or reference condition. The observed Pearson’s correlation coefficient, r, and the P-value of the statistically significant deviation from the null hypothesis of no correlation (r=0) were calculated. A multiple regression model was also applied (see statistical appendix in Supplementary Information). Enrichment analyses for gene categories For each sample in these experiments, we had the covariable GR and the corresponding expression profile per gene. The GR-gene expression covariation was quantified using a modified version of Pearson’s correlation (for details, see R. Tibshirani, G. Chu, B. Narasimhan and J. Li, 2011, SAM, Significance Analysis of Microarrays. R package version 2.0. http://CRAN.R-project.org/package= samr) between both values. A gene-to-gene analysis of the differential expression was performed by the Significance Analysis of Microarrays (SAM) method (33,34) with a false discovery rate of q <0.05. A gene set enrichment analysis was applied to previously detected sets of differentially expressed genes. A unilateral Fisher’s exact test was applied where the gene sets to be compared were the Gene Ontology groups. Analyses were run with R packages (R. Gentleman, with contributions from S. Falcon and D. Sarkar, 2014, Category Analysis. R package version 2.32.0 and M. Carlson, 2014), org.Sc.sgd.db: Genome wide annotation for Yeast. R package version 3.0.0). The whole R code used in this paper is found as Supplementary Material. A second enrichment analysis was run by the Gene Set Analysis (GSA) method (B. Efron and R. Tibshirani, 2010 GSA:, Gene set analysis. R package version 1.03. http://CRAN.R-project. org/package=GSA)(35). A detailed statistical protocol is given as a statistical appendix in Supplementary Information. Gene ontology representations The genes from the global experiment sets for all the parameters (TR, RA, and kd) from both data sources were ranked according to their correlation with the GR. The over-represented functional categories (Gene Ontology) obtained in the SAM and GSA analyses were reduced and visualized with the ReviGO web server ((36)http://revigo.irb. hr/). Calculating the respiratory quotient (RQ) index The protocol described in Quir´ os et al. (37) was used with the following modifications. RQ determinations were made in MiniBio bioreactors (250 ml nominal volume) (Applikon, Schiedam, The Netherlands) equipped with Peltierrefrigerated gas condensers, filled with YPD broth (200 ml). Cultures were sparged with air at a gas flow of 600 ml/h. 3646 Nucleic Acids Research, 2016, Vol. 44, No. 8 The gas flow was controlled with MFC17 mass flow controllers (Aalborg, Orangeburg, NY, USA), whose calibration was verified regularly by a soap bubble flowmeter. The instant CO2and O2concentrations in the exhaust gas were recorded every minute in a BlueInOne gas analyzer (Bluesens, Herten, Germany). The measured percentages of both gases were converted into instant production or consumption rates (expressed in moles of gas per liter of culture hourly) by taking into account the contribution of CO2production and O2consumption to the total gas flow. The accumulated production or consumption of each gas was determined as the integral over time of the instant values. The time point that corresponded to 0.15 moles of CO2produced per liter of culture was chosen as a way to standardize the RQ values for comparisons. At this time point, most of the initial sugar was already consumed, thus the corresponding RQ values (RQ 0.15) would reasonably summarize all the culture stages. RQ 0.15 was calculated as the ratio between the accumulated amounts (expressed in moles) of the produced CO2and the consumed O2for this time point. Cell microencapsulation and microcolony assay The microencapsulation of individual yeast cells in alginate microspheres (no more than one cell per microparticle) was performed using a Cellena microencapsulator (Ingeniatrics), as described elsewhere (38). Encapsulated cells were incubated together in YPD medium under standard culturing conditions. After ethanol fixation, microcolonies were analyzed and sorted using a BioSorter large-particle flow cytometer (Union Biometrica). A detailed protocol of the sorting procedures is described in (39). Transcriptomic analysis of microcolonies The total RNA extracted from the sorted microcolonies and from the non-sorted controls was sequenced following the 3T-fill protocol as previously described (40). The starting material used was 500–3000 ng of total RNA. The raw next generation sequencing reads to be first assessed for their quality with the FASTQC tool kit Bad quality reads (phred score <20) were trimmed. The genome sequence of S. cerevisiae S288c (version R64) and its annotations were retrieved from SGD (Saccharomyces Genome Database, http://yeastgenome.org/) and used for all the analyses. Raw reads were aligned against the reference genome using BWA/Bowtie2. We quantified the gene expression from the mapped reads using the HTSeq-count package (41) in order to obtain counts of mapped reads per gene (in the ‘intersection-nonempty’ mode). As the count process was dependent on sequencing depth, those samples with a small number of reads were removed, whereas those samples with a huge number of reads were down-sampled using the FastqSampler utility from the Shortread package (42). The statistical environment R (version 3.1; (43)), was used to perform the statistical analysis. Hierarchical clustering, boxplots, and Multidimensional Scaling (MDS) analyses were performed before and after normalization to measure the differences among samples. The edgeR package from R/BioConductor (44) was used to normalize and fit count data for the differential gene expression analysis. In order to avoid a bias in poorly expressed genes, those genes with less than five mapped reads per million in at least two samples were removed. An isoform-specific analysis was performed and an assignment was made to transcripts using GSNAP (version 2012–01–11) as previously described (40). RESULTS Transcription rate scales with the yeast growth rate To analyze the relationship between mRNA turnover and growth rate, we analyzed a set of 42 genome-wide transcription rate experiments (25,45–48)(P ´ erez-Ort´ ın, Medina and Ant´unez, unpublished), performed by the genomic run-on (GRO) method that provides nascent TR data (25). These experiments corresponded to the yeast strains (wild-types and mutants) that grow in different carbon sources (YPD, YPGal, YPRaff) and at distinct temperatures (from 23◦Cto 37◦C; see Supplementary Table S1A). In each experiment, we obtained individual mRNA amounts, nascent TR and kd for most protein-coding genes, as well as the total TR (for all three nuclear RNA polymerases). All the experiments were comparable because the data were normalized by cell numbers and cell volumes (48). We first analyzed the global RNA pol II nascent TR and RA data by summing all the protein-coding genes with at least 32 valid data points (5411 genes). Figure 1A shows a clear positive correlation between the total RNA pol II TR and the growth rate (in red). Within the analyzed GR range, an increase in the TR was approximately proportional to an increase in the GR. This dependence was not caused by the genes that encoded the translation machinery because the overall tendency remained almost identical when they were excluded (Figure 1A in blue). The total TR (including RNA pol I, II and III), which was also obtained in the runon experiment, correlated positively with the growth rate (Supplementary Figure S1). This was expected because approximately a 75% of the total TR came from RNA pol I + III transcription (8), which was done to transcribe the noncoding components of the translation machinery (rRNAs and tRNAs). Given the very diverse set of strains and growth conditions (see Supplementary Table S1A), they all revealed their own particular influence on the total TR. Thus the observed growth rate dependence of the TR must be regarded as remarkably robust. Nonetheless, all the data were obtained by a particular technique, GRO, which could be subject to technique-specific systematic biases. To independently confirm our results, we used an independent published data set based on metabolic mRNA labeling, which directly measures the mature mRNAs that appear in the cytoplasm. We analyzed the data that corresponded to 44 yeast mutant populations that grew exponentially at 30◦C in YPD from which the TR, RA and kdfor most of the protein-coding genes were extracted (26). To complement this data set, we measured the GRs of most of these mutants (38 strains, see Supplementary Table S1B) by micro-cultivation (30). As the TR should be measured as a change in mRNA concentration over time (22), and the original data were not corrected, we used the volume estimates of these mutants obtained by Nucleic Acids Research, 2016, Vol. 44, No. 8 3647 Figure 1. The mRNA turnover correlated positively with the yeast growth rate. (A) We plotted 34 of the 42 different growth variations (mutant strains, temperature or carbon source variations) for which we had relative transcription rates to a reference sample. A clear dependence of the total nascent transcription rate (TR) of RNA pol II is observed, as measured by GRO (23) with the GR of the culture in aerated flasks. Values are relative to the corresponding reference strain. Red symbols correspond to the average of all the data, whereas blue symbols correspond to the average of all the data except genes RP and RiBi. Pearson’s correlation coefficient (r) and the associated P-value are shown. (B–D). We used the data for the average total mRNA synthesis rate (TR), mRNA levels (RA) or mRNA degradation constant (kd) from 38 mutant strains (26) and the GR determined by micro-cultivation. The original TR and RA data were corrected by cell volume using data from (27,28). Once again, we observe a positive dependence of TR and in kd,which kept the mRNA levels approximately constant. Pearson’s correlation coefficients (r) and associated P-values are shown. (E) Multiple Regression Model shows that the RA can be predicted from the TR and the kdwith highly significant P-values but it cannot be predicted from GR (see Statistical Appendix in Supplementary Information for further explanation). 3648 Nucleic Acids Research, 2016, Vol. 44, No. 8 Jorgensen et al. (27) and Truong et al. (28) to correct them. The plots for the three parameters are shown in Figure 1B– D. Here we see that the TR steadily and proportionally increased with an increasing GR (Figure 1B), whereas mRNA stability decreased (kdincreased, Figure 1D). The reverse dependencies of the TR and mRNA stability on growth rate ensured that the steady-state mRNA levels, which reflected the joint consequence of mRNA synthesis and degradation rates, were largely independent of the GR (Figure 1C). A Multiple Regression Model confirms that RA can be predicted from the TR and from the kd, as expected, with highly significant P-values but it cannot be predicted from the GR (Figure 1E). Similar results for the kdand RA dependence on GR were obtained for our GRO experiments data set from a subset of experiments (see Supplementary Table S4). Together these results compellingly support that the RNA pol II transcription rates, measured as both the nascent or mature TR, were strongly associated with the growth rate. The degradation rate constant (kd,seeIntroduction section) also increased with a rising growth rate, and total [mRNA] remained constant. Thus the mRNA turnover rate of yeast cells correlated directly with the population growth rates. Moreover, the increase in mRNA turnover was quantitatively similar to that of the GR, which suggests that the proportionality between both should be maintained. Protein synthesis genes are more transcribed and respiratory gene transcripts are destabilized at fast growth The general tendency of the RNA pol II TR represents the average over all protein-coding genes. The particular tendencies of individual genes may diverge substantially from this average and can be studied from our original data sets or from that published by P. Cramer’s group (26). To investigate the dominant biological trends among individual genes for each TR, kdand RA, we ordered all the normalized data according to their correlation with the growth rate. Two different statistical analyses, the SAM and GSA methods (see M&M), were then applied to identify the gene functions (GO categories) enriched in the top (highest positive r) or the bottom (lowest positive r or highest negative r, depending on the plot) of each ranking list. Both methods and both data sets gave similar results. As the general TR and kdtendencies resulted in strong positive correlations with the GR (Figure 1), very few genes showed negative correlations. Thus for the TR, and kddescribed below, ‘top’ means the genes with the strongest positive correlation, whereas ‘bottom’ generally indicates the genes with the lowest positive (around zero) correlation. In contrast, the RA data sets are centered on zero given the general absence of a change withaGR(Figure1C). Therefore the ‘bottom’ group contains the genes that truly and negatively correlate with the growth rate. As all the data sets were normalized by their median, the found tendencies represent biases of the functional categories with regard to the typical behavior of the population. The SAM statistical analyses of TR-GR dependence found RiBi (Ribosome Biogenesis), and other proteinsynthesis-related genes, to be enriched among the genes whose nascent TRs correlated positively with the GRs (Figure 2A, Supplementary Table S2). In contrast, the TRs of the respirationand mitochondria-related genes tended to be the least positively correlated with GRs (Figure 2B). The statistical analyses of RA-GR dependence also found that the steady-state mRNA abundance of RiBiand the other protein synthesis-related genes tended to increase with rising growth rates (Figure 2C, Supplementary Table S2). In contrast, the steady-state mRNA abundance of the respirationand mitochondria-related genes diminished with faster growth (Figure 2D). We did not analyze the mRNA stability tendencies in this data set because very few data sets were available. Once again, the analysis of the TR-GR dependence in the alternative data set (26) found that the mRNA synthesis rate of the RiBiand amino acid biosynthesis-related genes correlated strongly and positively with growth rate, whereas the TRs of the cell cycle and regulation genes tended to decouple from it (see Supplementary Figure S2 and Table S3). Given the larger number of estimates in this case, we also analyzed kdGR dependence to find that the mRNA degradation of the RiBi and Glycosylation genes did not correlate with growth rate, whereas the mRNA degradation of the vacuoleand mitochondria-related genes tended to increase with rising GR (see Figure 2E–F and Supplementary Table S2). Other alternative searches between our TR and RA data and the TRandRAdataofSunet al. (26) were made by the GSA method. Similar results were obtained. D. Botstein’s lab made a systematic search for genes whose RA correlated either positively or negatively with the growth rate (16,17,49). In these studies, a signature of 72 genes was obtained, which strongly depended on the GR. We further analyzed the response of these genes and found that some 50 genes showed essentially identical (decreasing or increasing) tendencies of their RA and TR profiles (Supplementary Figure S3). This indicates that most growth rate-regulated transcripts are controlled manly at the synthesis level. However, we also found seven genes for which mRNA abundance and the TR rate data correlated inversely, which suggests a strong regulation at the stability level. In conclusion, apart from the general tendency of cells to increase the total mRNA synthesis rate with faster growth, the increased relative transcription of some functional classes of genes was particularly pronounced. As a result of the higher TRs, the steady-state mRNA abundance of these particular types of genes, unlike that observed for the overall transcriptome (Figure 1C), seemed to also increase with faster growth. In contrast, the steadystate mRNA concentration of other functional classes of genes actually decreased as a result of faster growth, caused apparently by a particular relative decrease in their mRNA stability (increase in kd). Ribosome biogenesis, mitochondria-related and stressinduced genes use different mechanisms to adjust mRNA levels to the growth rate We previously found that some GO categories differed from the general tendency of no change in mRNA abundance by displaying changes in the GR. To shed further light onto these opposing trends, we separately considered the ri- Nucleic Acids Research, 2016, Vol. 44, No. 8 3649 Figure 2. Search for the Gene Ontology categories that show enhanced dependence on the growth rate. GO enrichment searches for the highest positive and negative GR dependences on the nascent TR, mRNA levels and degradation rates are shown. Apart from the general tendency of the nTR and RA showninFigure1, certain functional categories are especially enriched among the highest positive correlations (top) with the GR, or among the genes with the lowest correlation (bottom). We show here (A–D) the analyses done with the data sets used in Figure 1for the nTR and RA, and from the data set of Sun et al. (26)forkd(E–F) using an SAM analysis and summarized by the ReViGO program. Similar results for the TR and RA were obtained with the data sets used in Figure 1B–C (26); see Supplementary Figure S3. 3650 Nucleic Acids Research, 2016, Vol. 44, No. 8 bosome biogenesis (GO:0042254) and mitochondrial ribosome (GO:0005761) functional categories. Figure 3A illustrates how the TR of the RiBi genes increased to result in higher mRNA abundance and thereby, presumably, maintained the ribosome numbers level with GR demands. This is the clearest example of a regulon controlled at the transcriptional level. RiBi mRNA stabilities did not change in a GR-dependent manner (Figure 3B). This behavior was similar in ribosomal protein (RP) genes (Supplementary Figure S4). In contrast, the mitochondrial ribosome genes showed no regulation at the transcriptional level (Figure 3C). Nevertheless, steady-state mRNA abundance diminished with the GR, which was achieved by a marked reduction in the stability of these mRNAs (Figure 3D). This is the clearest example of a post-transcriptional regulon (50). Given that the mRNA levels for the respiratoryand other mitochondria-related genes lowered with faster growth, and by assuming no additional compensatory regulation of the translation or post-translational levels, we expected respiration to decrease accordingly. To test this prediction, we selected a set of nine strains (mutants and wild types) from our data set (see Figure 1A and Supplementary Table S1A) and determined their respiratory quotient (fermentation/respiration, RQ) in a microfermentor (37). As seen in Supplementary Figure S5, the lower the mRNA levels of the genes related with respiration (GO:0009060) or mitochondrial ribosome (GO:0005761), the higher the RQs and, thus, the relative influence of respiration on metabolism. This shows that the reduction in mRNA levels of the respiratory and other mitochondrial genes with higher growth rates is physiologically relevant and serves to reduce respiration. This strongly supports the dependence of faster-growing yeast cells on the fermentative metabolism. Finally, we paid special attention to the frequently studied stress-induced genes of the so-called ‘environmental stress response’ (ESR) (51). These genes have been shown to be up-regulated in the RA in slow-growing mutants (12,16), which implies a negative correlation between their mRNA abundance and growth. It was our intention to establish whether this was indeed the case, and if so, whether or not it was achieved by a change in the TR or DR. Indeed the mRNA abundance of the stress-induced section of the ESR diminished at faster growth (Figure 3E). This decrease was achieved exclusively by lowering TRs, and no general change took place in mRNA degradation (Figure 3F). Thus the stress-induced part of the ESR exemplified a third regulatory strategy: a negative effect of growth rate on transcription rates. mRNA dependence on the growth rate can be detected within the heterogeneity of a wild-type cell population So far we analyzed the dependence of mRNA concentrations, and their synthesis and degradation rates on cellular growth rate across different genotypes and environmental conditions, including mutant strains with impaired proliferation capacity or wild-type cells that grow in different environments. Similarly, previous experiments that have described the influence of the GR on gene expression have also been based on the comparison of cells populations in different carbon sources or under chemostat-regulated conditions (18,52). These experimental designs cannot rule out the fact that the seeming dependence of TR, RA and kdon the GR might be a confounding effect of both depending on a shared third factor in the form of specific mutations and specific adverse environmental effects. The connection between the growth rate and RNA metabolism should happen at the single cell level. In order to rule out that any of those effects occurred only at the population level, we measured the GR of single cells as their ability to produce microcolonies, and we analyzed their GR transcriptomes. To do this, we analyzed the microcolonies generated after the encapsulation of the single cells that belonged to a clonal wild-type population in alginate microspheres (Supplementary Figure S6). Alginate encapsulation allows the free diffusion of nutrients, is fully compatible with yeast proliferation, and helps sort cells according to their GR potential (38). No decreased viability was detected after the encapsulation procedure (Supplementary Figure S7A). Encapsulated cells were incubated in rich liquid medium under standard culture conditions and were allowed to proliferate. At the different time points, samples were taken and observed by optical microscopy. The very different size observed in microcolonies indicated an evident variation in the proliferation capacity within this clonal population (Figure 4A). Quantification of the microcolony size by large particle flow cytometry confirmed an intrinsic heterogeneity that increased with incubation time (Figure 4B). We tested whether the different microcolony size was due to a distinct size of the microparticle itself or if it was a biased localization within the microparticle. We found no difference for either (Supplementary Figure S7B and S7C). Microcolony size heterogeneity could also be due to an unequal lag phase. The increasing statistical dispersion shown by microcolony size during the time-course experiments, by maintaining symmetrical distribution, was hardly compatible with this possibility (Figure 4B). Our results are in agreement with previous reports on proliferation heterogeneity conducted under optimal culture conditions in S. cerevisiae, detected with technical approaches that did not involve encapsulation (24). The combination of microencapsulation and large particle flow cytometry allows sister microcolonies to be physically separated by sorting (Supplementary Figure S6). We isolated the 10% top and bottom subpopulations of microcolonies according to their size (Supplementary Figure S6). These two subpopulations exhibited an average cell number of 32 and 2519, respectively, which involves apparent doubling times of approximately 162 and 76 min according to the incubation time under growing conditions (Figure 4C). This longer duplication time of the small microcolony cells agrees with their longer G1 phase and reflects their differential cell-cycle regulation (Delgado–Ramos, Mu˜ noz– Centeno and Ch´ avez, to be published elsewhere). The average cell number and apparent doubling time of non-sorted microcolonies gave intermediate values. However, when all the analyzed microcolonies were considered as a single population, the resulting population doubling time was comparable to the standard liquid cultures (Figure 4C). This simple procedure allowed us to isolate mRNA preparations from the two types of sorted colonies, and to Nucleic Acids Research, 2016, Vol. 44, No. 8 3651 Figure 3. Ribosome biogenesis, mitochondria-related and stress-induced genes use different mechanisms to adjust mRNA levels to the growth rate. Using the data from the Figure 1A data set, we extracted data for the RiBi regulon (GO: 0042254, panels Aand B), for the mitochondrial ribosome (GO: 0005761, panels Cand D), and for the induced part of the ESR (51)(panels Eand F). The genes in the Mitochondrial ribosome GO category (and similarly in the Puf3 regulon, see Figure 5B) seem to have lower mRNA levels (RA) with the growth rate (GR) of the culture, mainly by mRNA destabilization (increased degradation constant, kd). The mRNA levels of the RiBi genes (and the ribosomal protein genes, see Supplementary Figure S4) tend to rise with the GR due to increasing transcription (nTR). The induced part of the ESR conversely lowers, and RA reduces with decreasing transcription. Pearson’s correlation coefficients (r) and associated P-values are shown. 3658 Nucleic Acids Research, 2016, Vol. 44, No. 8 54. Lesage,P., Yang,X. and Carlson,M. (1996) Yeast SNF1 protein kinase interacts with SIP4, a C6 zinc cluster transcriptional activator: a new role for SNF1 in the glucose response. Mol. Cell. Biol.,16, 1921–1928. 55. Vincent,O. and Carlson,M. (1998) Sip4, a Snf1 kinase-dependent transcriptional activator, binds to the carbon source-responsive element of gluconeogenic genes. EMBO J.,17, 7002–7008. 56. Gasmi,N., Jacques,P.E., Klimova,N., Guo,X., Ricciardi,A., Robert,F. and Turcotte,B. (2014) The switch from fermentation to respiration in Saccharomyces cerevisiae is regulated by the Ert1 transcriptional activator/repressor. Genetics,198, 547–560. 57. Miller,M.A., Russo,J., Fischer,A.D., Lopez Leban,F.A. and Olivas,W.M. (2014) Carbon source-dependent alteration of Puf3p activity mediates rapid changes in the stabilities of mRNAs involved in mitochondrial function. Nucleic Acids Res.,42, 3954–3970. 58. Lee,C.D. and Tu,B.P. (2015) Glucose-regulated phosphorylation of the PUF protein Puf3 regulates the translational fate of its bound mRNAs and association with RNA granules. Cell Rep.,11, 1638–1650. 59. Gupta,I., Clauder-Munster,S., Klaus,B., Jarvelin,A.I., Aiyar,R.S., Benes,V., Wilkening,S., Huber,W., Pelechano,V. and Steinmetz,L.M. (2014) Alternative polyadenylation diversifies post-transcriptional regulation by selective RNA-protein interactions. Mol. Syst. Biol.,10, 719. 60. Riordan,D.P., Herschlag,D. and Brown,P.O. (2011) Identification of RNA recognition elements in the Saccharomyces cerevisiae transcriptome. Nucleic Acids Res.,39, 1501–1509. 61. Saint-Georges,Y., Garcia,M., Delaveau,T., Jourdren,L., Le Crom,S., Lemoine,S., Tanty,V., Devaux,F. and Jacq,C. (2008) Yeast mitochondrial biogenesis: a role for the PUF RNA-binding protein Puf3p in mRNA localization. PloS One,3, e2293. 62. Porrua,O. and Libri,D. (2015) Transcription termination and the control of the transcriptome: why, where and how to stop. Nat. Rev. Mol. Cell Biol.,16, 190–202. 63. Pelechano,V., Chavez,S. and Perez-Ortin,J.E. (2010) A complete set of nascent transcription rates for yeast genes. PloS One,5, e15442. 64. Keren,L., van Dijk,D., Weingarten-Gabbay,S., Davidi,D., Jona,G., Weinberger,A., Milo,R. and Segal,E. (2015) Noise in gene expression is coupled to growth rate. Genome Res.,25, 1893-1902 65. Shalem,O., Groisman,B., Choder,M., Dahan,O. and Pilpel,Y. (2011) Transcriptome kinetics is governed by a genome-wide coupling of mRNA production and degradation: a role for RNA Pol II. PLoS Genet.,7, e1002273. 66. Molin,C., Jauhiainen,A., Warringer,J., Nerman,O. and Sunnerhagen,P. (2009) mRNA stability changes precede changes in steady-state mRNA amounts during hyperosmotic stress. RNA,15, 600–614. 67. Perez-Ortin,J.E., Alepuz,P.M. and Moreno,J. (2007) Genomics and gene transcription kinetics in yeast. Trends Genet.,23, 250–257. 68. Gerber,A.P., Herschlag,D. and Brown,P.O. (2004) Extensive association of functionally and cytotopically related mRNAs with Puf family RNA-binding proteins in yeast. PLoS Biol.,2, E79. 69. Diaz-Ruiz,R., Rigoulet,M. and Devin,A. (2011) The Warburg and Crabtree effects: on the origin of cancer cell energy metabolism and of yeast glucose repression. Biochim. Biophys. Acta,1807, 568–576. 70. Lagunas,R. (1986) Misconceptions about the energy metabolism of Saccharomyces cerevisiae. Yeast,2, 221–228. 71. Dashko,S., Zhou,N., Compagno,C. and Piskur,J. (2014) Why, when, and how did yeast evolve alcoholic fermentation? FEMS Yeast Res., 14, 826–832. 72. Wang,J., Atolia,E., Hua,B., Savir,Y., Escalante-Chong,R. and Springer,M. (2015) Natural variation in preparation for nutrient depletion reveals a cost-benefit tradeoff. PLoS Biol.,13, e1002041. 73. Venturelli,O.S., Zuleta,I., Murray,R.M. and El-Samad,H. (2015) Population diversification in a yeast metabolic program promotes anticipation of environmental shifts. PLoS Biol.,13, e1002042. 74. Abernathy,E., Gilbertson,S., Alla,R. and Glaunsinger,B. (2015) Viral Nucleases induce an mRNA degradation-transcription feedback loop in mammalian cells. Cell Host Microbe,18, 243–253.