1. Genetics and Genomics
Download icon

Single cell functional genomics reveals the importance of mitochondria in cell-to-cell phenotypic variation

  1. Riddhiman Dhar
  2. Alsu M Missarova
  3. Ben Lehner  Is a corresponding author
  4. Lucas B Carey  Is a corresponding author
  1. The Barcelona Institute of Science and Technology, Spain
  2. Universitat Pompeu Fabra, Spain
  3. Indian Institute of Technology, India
  4. Institució Catalana de Recerca i Estudis Avançats (ICREA), Spain
Research Article
  • Cited 1
  • Views 1,893
  • Annotations
Cite this article as: eLife 2019;8:e38904 doi: 10.7554/eLife.38904

Abstract

Mutations frequently have outcomes that differ across individuals, even when these individuals are genetically identical and share a common environment. Moreover, individual microbial and mammalian cells can vary substantially in their proliferation rates, stress tolerance, and drug resistance, with important implications for the treatment of infections and cancer. To investigate the causes of cell-to-cell variation in proliferation, we used a high-throughput automated microscopy assay to quantify the impact of deleting >1500 genes in yeast. Mutations affecting mitochondria were particularly variable in their outcome. In both mutant and wild-type cells mitochondrial membrane potential – but not amount – varied substantially across individual cells and predicted cell-to-cell variation in proliferation, mutation outcome, stress tolerance, and resistance to a clinically used anti-fungal drug. These results suggest an important role for cell-to-cell variation in the state of an organelle in single cell phenotypic variation.

https://doi.org/10.7554/eLife.38904.001

Introduction

Isogenic populations often exhibit considerable phenotypic heterogeneity even in an identical environment. One common phenotypic variation that has been observed in isogenic populations of microbial and mammalian cells, including cancer cells is variation in proliferation rate (Sandler et al., 2015; Fridman et al., 2014; Levy et al., 2012; Yaakov et al., 2017; Ferrezuelo et al., 2012; Polymenis and Schmidt, 1999; Gupta et al., 2011; Kiviet et al., 2014; Roesch et al., 2010). Phenotypic variations that are often coupled with variation in proliferation rate are the abilities of an individual cell to survive stress and drug treatment (Levy et al., 2012; Yaakov et al., 2017). In this regard, the existence of ‘persister’ cells in microbial populations is well known and poses a significant challenge for antibiotic treatment (Fridman et al., 2014; Balaban et al., 2004; Kussell et al., 2005; Wakamoto et al., 2013; LaFleur et al., 2006; Bojsen et al., 2017). Similarly, individual cells in tumors have been shown to vary in their ability to survive anticancer drugs and can lead to drug-resistant populations (Shaffer et al., 2017; Sharma et al., 2010; Ramirez et al., 2016; Hata et al., 2016; Rego et al., 2017; Márquez-Jurado et al., 2018). Recent advances in single-cell techniques are revealing the extent of transcriptomic and metabolic differences among isogenic cells (Dey et al., 2015; Trapnell, 2015). The existence of such heterogeneity in gene expression in isogenic microbial and animal populations has been shown – to some extent – to underlie the variable outcome of mutations (Dickinson et al., 2016; Burga et al., 2011; Raj et al., 2010; Eldar et al., 2009; Horvitz and Sulston, 1980). Incomplete mutation penetrance and variable expressivity is also common in human disease (Cooper et al., 2013; Zlotogora, 2003; Giudicessi and Ackerman, 2013; Otterson et al., 1999).

Heterogeneity can arise due to stochastic fluctuations in biological processes taking place inside cells. This can happen due to the small numbers of molecules involved in processes such as transcription (Berg, 1978; Rigney, 1979; Cookson et al., 2010) or during stochastic partitioning of cellular components during cell division (Birky and Skavaril, 1984; Huh and Paulsson, 2011). Although genetic variation has been shown to influence proliferation heterogeneity (Ziv et al., 2013) and cell-to-cell variation in the expression level of single genes and levels of metabolite has been correlated with variation in proliferation rate and stress and drug resistance (Levy et al., 2012; Yaakov et al., 2017; Shaffer et al., 2017; Burga et al., 2011; Raj et al., 2010; Eldar et al., 2009; Li et al., 2018; Rotem et al., 2010; Battich et al., 2015), the true underlying causes of such phenotypic heterogeneity are poorly understood.

To identify genes and cellular processes involved in the generation of phenotypic heterogeneity we set up a high-throughput microscopy assay to quantify proliferation heterogeneity in a yeast population. Using this assay, we quantify the impact of deletion of >1500 genes on proliferation heterogeneity. We present evidence that the variation in mitochondrial membrane potential is an important determinant of phenotypic heterogeneity in individual cells. We also show that mitochondrial membrane potential impacts gene expression and stress tolerance and drug resistance in individual cells. Taken together, our work suggests an important role for an organelle in generating phenotypic heterogeneity across individual cells in a homogenous environment.

Results

Natural and lab yeast populations show proliferation heterogeneity

To investigate cell-to-cell variation in proliferation rates, we set up a high-throughput automated time-lapse microscopy assay that measures the proliferation rates of thousands of single-cells per plate as they grow into micro-colonies. The assay uses a microscope with laser-based autofocus for image acquisition and a liquid handling robot to minimize density-dependent effects on proliferation. The data obtained are highly reproducible with mode proliferation rate of a lab strain being 0.407 ± 0.011 h−1, (mean ±sd) during >2 years of data collection (n = 44 batches; Figure 1A).

High-throughput analysis of single cell proliferation rate heterogeneity.

(A) High throughput microscopy setup – log phase yeast cells were diluted onto conA coated microscopy plate using Biomek NX liquid handling system to have similar cell density across wells. Cells were observed using an ImageXpress Micro system. Images were processed using custom scripts and data for area of microcolony vs. time were obtained. The points in the area vs. time graph show actual data and the solid lines show lowess fits. Data collected from all fields of view in a well constitute a microcolony proliferation rate distribution for a strain. The common lab yeast strain BY4741 (WT) has ~10% slow proliferating sub-population. The density shows mean density and the shaded areas in grey represent ±1 s.d. value at each point. The dotted red line shows the expected proliferation distribution if it were normally distributed. (B) Natural strains of yeast (Ziv et al., 2013) also have slow proliferating sub-populations. Each point represents data for one strain. Solid lines show median value. (C) WT strain re-created the original proliferation distribution even after 20 generations of growth. The plot shows data from two replicate measurements.

https://doi.org/10.7554/eLife.38904.002

Laboratory strains of the budding yeast Saccharomyces cerevisiae showed substantial cell-to-cell variation in proliferation, with ~10% of cells forming a slow growing sub-population in defined growth medium (Figure 1A) (Levy et al., 2012; Ziv et al., 2013). This slow growing sub-fraction is not unique to laboratory strains but exists in all natural and clinical isolates that we tested (Figure 1B; Supplementary file 1) (Ziv et al., 2013). Growth of the culture for an additional 20 generations did not alter the proliferation rate distribution; the mixture of slow and fast proliferating cells is maintained (Figure 1C). Proliferation is therefore a stable heterogeneous phenotype within a population, with the amount of heterogeneity depending on the genetic background.

A genome-scale screen to identify genes that alter proliferation heterogeneity

The effect of individual gene deletions on population-level growth rate has been well studied (Giaever et al., 2002; Baryshnikova et al., 2010). Many deletions have been shown to reduce population growth rate and can do so in different ways. Deletions can uniformly affect fitness of all the cells or alternatively, can affect fitness of a sub-population whereas the rest of the population remains unaffected. Inter-individual variation in the outcome of mutations has been observed before in multicellular organisms (Dickinson et al., 2016; Burga et al., 2011; Raj et al., 2010) but its relative occurrence has not been systematically quantified.

We therefore used the automated microscopy assay to quantify proliferation rate heterogeneity in triplicate for 1600 gene deletion mutants (Supplementary file 2, including 1150 gene deletions previously reported as affecting growth rates [Giaever et al., 2002; Baryshnikova et al., 2010]). We obtained reproducible data (where at least two replicate measurements showed good agreement) for 1520 deletions, with 1112 of these reducing the population proliferation rate in our experiment (Mann-Whitney U test, FDR < 0.1).

Deletion strains with similar population proliferation rates often showed strikingly different degrees of intra-population heterogeneity (Fig, 2A-C). At the single cell level,~39% of all mutants with a significant reduction in population proliferation rate (1112 mutants) showed significantly higher variation in mutation outcome compared to the WT strain (Mann-Whitney U test, FDR < 0.1, after correcting for change in mode growth rate). Among these mutants,~13% had the same mode growth rate as the WT strain while showing higher variability (Mann-Whitney U test, FDR < 0.1). However, almost all mutants (1111 of 1112 mutants) had a subset of cells proliferating at the same rate as the bulk of the wild-type (WT) population (one sample Wilcoxon rank-sum test for overlap with bulk WT distribution differing from zero, FDR < 0.1; Figure 2D, Figure 2—figure supplement 1). Thus, a highly variable outcome is actually the normal outcome for proliferation rate at the single cell level when a non-essential gene is inactivated (Figure 2D, Figure 2—figure supplement 2A).

Figure 2 with 2 supplements see all
Single cell proliferation rate distributions for 1500 gene deletions.

(A) Mode growth rate (h−1) and % slow fraction for 1520 deletion strains. The points represent average values across replicates and the bars represent ±1 s.d. values. The colours show classification of mutants into different categories according to change in mode growth rate (see Materials and methods, FDR < 0.1) and change in % slow fraction (FDR < 0.1) compared to the wild-type (WT) strain. The table and pie chart show the number and proportion of strains in each group (colour coded). Replicate data for WT strain are shown by multiple black points. (B) Examples of growth distributions of mutants classified into different groups which are colour coded as in A. The distribution in dark grey shows WT growth distribution. (C) Coefficient of variation (CV) vs. mean growth rate for all strains. WT values are shown in black; mutants of genes that localize to mitochondrial envelope in red. The points represent average values across replicates and the bars represent ±1 s.d. values. (D) % of WT-like cells in all mutants showing variable mutation outcome. It was calculated for all mutants showing significant reduction in mean proliferation rate and had significant proportion of cells growing as fast as the bulk of the WT proliferation distribution (Wilcoxon rank-sum test). (E) Functional class enrichment (GOslim) analysis for different classification groups show significantly enriched functional classes (hypergeometric test, FDR < 0.1). P – Biological Process, F – Molecular Function, C- Cellular Component. Bars show % of genes in a particular group (colour coded) being present in that particular functional class.

https://doi.org/10.7554/eLife.38904.004

Deletion of genes involved in mitochondrial function alter heterogeneity

To identify the determinants of this cell-to-cell variation in growth-rate and mutational impact we classified each of the deletions by how it affected both the mode and distribution of cellular proliferation rates (Supplementary file 2, Figure 2A,B). Approximately, 17% of the mutants showed no change in either mode proliferation rate or percentage of slow sub-population (in grey), whereas ~43% exhibited a change in mode proliferation rate but no change in slow fraction (in light blue). Interestingly, 48 mutants reduced the slow fraction without any change in mode proliferation rate (in red) and 97 mutants increased the slow fraction without altering the mode proliferation rate (in blue). In addition, there were 78 mutants that reduced both the slow fraction and the mode proliferation rate (in orange). Finally, 370 mutants reduced the mode growth rate but increased the slow sub-population (in magenta, Figure 2A). Across mutants, we observed a strong inverse relationship between mean growth rate and noise (co-efficient of variation) (Figure 2C), as has been observed for gene expression (Newman et al., 2006; Bar-Even et al., 2006).

To identify biological processes associated with changes in the slow growing sub-population, we performed a GO functional enrichment analysis on genes in these categories (FDR < 0.1). Deletions causing the largest increase in the fraction of slow proliferating cells were highly enriched for nuclear genes encoding mitochondrial proteins (Figure 2C,E). Among the mutants that increased the slow fraction but also reduced mode growth rate (Figure 2E, magenta), ~30% localized to mitochondria (~1.2 fold enrichment),~13% localized to the mitochondrial envelope (>1.6 fold enrichment) and ~4.6% were involved in cellular respiration (~2-fold enrichment). In particular, deletion of genes that localized to the mitochondrial envelope resulted in a large increase in slow fraction and noise (Figure 2E, Figure 2—figure supplement 2B, Supplementary file 2). Mutations that affect mitochondria, and in particular the mitochondria membrane, increase heterogeneity, suggesting that heterogeneity in proliferation might be associated with cell-to-cell variation in mitochondria. Among genes with other functional associations, deletion of genes with kinase activity (STE11, SNF1, NPR1, ATG1, HXK2), genes with cytoskeletal protein binding capability (HOF1, SRV2, SHE1) and genes associated with carbohydrate transport (SNF3, HXK2) reduced percentage of slow growing cells but did not alter the mode growth rate.

Mitochondrial membrane potential but not amount predicts slow growth

To further investigate the role of mitochondria in proliferation heterogeneity, we used the MitoTracker dye to quantify mitochondrial amount in WT cells and five deletion strains with very different proliferation distributions. Total mitochondria amount varied little across the strains (Figure 2—figure supplement 2C), ruling out cell-to-cell variation in segregation of the organelle as a driver of heterogeneity. However, signal from the mitochondria membrane potential dye TMRE varied substantially across WT, mutants and natural strains (Figure 3A,B). This suggested that the mitochondria membrane potential – but not their amount – might be driving proliferation heterogeneity.

Figure 3 with 2 supplements see all
Variation in mitochondria potential across single cells underlies proliferation heterogeneity.

(A) TMRE stain intensity (log transformed) measured by flow cytometry in WT and deletion mutants. (B) TMRE intensity in WT and natural isolates of S. cerevisiae strains. (C) WT cells were sorted by TMRE signal intensity into four bins HI, M1, M2 and LO with gates as shown (~5% of the population sorted in each bin) and growth rate distributions were measured using high throughput microscopy setup. HI bin was enriched for slow growing cells. (D) % of respiration deficient cells in each bin from WT strain. The columns represent the average values from 12 independent experiments and the bars show ±1 s.d. values. (E) Percentage of respiration deficient cells in WT and mutant strains is positively correlated with the percentage of slow growing cells. The blue dotted line represents y = x line. The error bars represent ±1 s.d. measured from at least two biological replicates for each strain. (F) Percentage of respiration deficient cells in UCC strains (Dimitrov et al., 2009) is strongly positively correlated with the percentage of slow growing cells. The blue dotted line represents y = x line. The error bars represent ±1 s.d. measured from at least two biological replicates for each strain.

https://doi.org/10.7554/eLife.38904.008

To determine if mitochondrial membrane potential is correlated with variation in growth rate within a population, we sorted wild-type cells according to their TMRE staining and measured the fraction of slow proliferating cells. The population with high TMRE was highly enriched for slowly proliferating cells (Figure 3C). This same fraction was also strongly enriched for respiration deficient cells (Figure 3D). This strong enrichment for slow proliferation and respiration deficiency was also observed for the cells with high TMRE in gene deletion strains and in natural isolates (Figure 3—figure supplement 1A–D).

We further quantified the TMRE signal and proliferation distribution in a set of twelve strains that differ only by naturally occurring polymorphisms known to affect mitochondrial function by altering mtDNA inheritance (Dimitrov et al., 2009). Across all datasets, the percentage of slowly proliferating cells showed a high correlation with the percentage of respiration-deficient cells (Figure 3E,F) as well as with the percentage of high TMRE cells (Figure 3—figure supplement 2A).

Although cell-to-cell variation in mitochondrial amount did not predict proliferation rate variation (Figure 3—figure supplement 2B), mtDNA copy number was substantially lower in the cells with high TMRE (Figure 4A; Figure 4—figure supplement 1A,B), suggesting a likely role of mtDNA copy number in defining mitochondrial membrane potential and ultimately, in generation of growth rate heterogeneity.

Figure 4 with 1 supplement see all
Reduction in mtDNA copy number causes slow growth.

(A) mitochondrial DNA (mtDNA) copy number in the sorted bins HI-LO from WT strain measured through quantitative PCR. Two columns show results from two independent experiments. The column represents average mtDNA copy number calculated based on five pairs of primers binding mtDNA and five pairs of primers binding nuDNA and three technical replicates for each of these primers. The bars show ±1 s.d. values. (B) Overexpression of Mip1 gene in WT strain led to significant increase in mtDNA copy number (C) Overexpression of Mip1 gene led to significant reduction in percentage of respiration deficient cells and in slow growing subpopulation in WT strain and tim11Δ mutant. Data are from at least four biological replicates. (D) Overexpression of MIP1 gene in WT strain reduced percentage of cells with high TMRE signal. (E) Percentage of slow growing sub-population was strongly correlated with mtDNA copy number in mutant strains. The dotted lines represent values for WT strain. The error bars represent ±1 s.d. values. (F) Pre-growing WT strain overnight in medium containing ethanol as sole carbon source (that required respiration) reduced percentage of slow growing sub-population by ~50% compared to pre-growth in medium containing glucose as the sole carbon source. Data are from six biological replicates.

https://doi.org/10.7554/eLife.38904.012

To establish a causal relationship between mtDNA copy number and slow growth, we introduced an extra copy of the mitochondrial DNA polymerase Mip1 (Genga et al., 1986), which increased mtDNA copy number 3-fold (Figure 4B). This reduced both the fraction of slow proliferating and respiration-deficient cells and the fraction of cells with high TMRE signal (Figure 4C,D), suggesting that variation in mtDNA copy number can be causal for variation in both mitochondrial membrane potential and proliferation. Consistent with an effect of mtDNA copy number on growth, knocking out of Mip1 gene led to complete loss of mtDNA and resulted in completely slow growing yeast population compared to WT (Figure 4—figure supplement 1C,D) (Genga et al., 1986). Furthermore, mtDNA copy number showed a strong negative correlation with the percentage of slow proliferating cells across mutants (Figure 4E, Figure 4—figure supplement 1E). Finally, forcing cells to respire by pre-growing them on ethanol as the sole carbon source prior to growth in glucose decreased the fraction of slowly proliferating cells (Figure 4F). Taken together, these results suggest that alterations in mitochondrial membrane potential, which can be caused by mtDNA copy number reduction below a threshold and other mechanisms is the underlying cause of slow growth in individual cells.

Variation in mitochondrial membrane potential predicts additional phenotypic heterogeneity including drug resistance

To systematically understand the physiological differences between sub-populations that vary by mitochondrial membrane potential, we analyzed the transcriptome by RNA sequencing. Cells with high TMRE have low expression of respiratory and proliferation-associated genes (Figure 5A, Figure 5—figure supplement 1A,B). Consistent with previous analyses of slow proliferating cells (Yaakov et al., 2017; van Dijk et al., 2015), they also exhibit a DNA damage response (Figure 5B) and signs of iron starvation (Figure 5C), which has previously been reported for respiration deficient cells (Veatch et al., 2009; Puig et al., 2005). Cells with intermediate TMRE have very similar proliferation distributions to cells with low potential (Figure 3C). However, their gene expression was substantially different (Figure 5—figure supplement 1C,D), including reduced expression of respiratory genes in cells with intermediate TMRE (Figure 5—figure supplement 1D).

Figure 5 with 12 supplements see all
Cell-to-cell variation in mitochondria potential predicts single cell drug resistance.

(A) Heatmap of expression of respiration genes in cells sorted by their TMRE signal intensity (bins HI-LO). (B) Heatmap of expression of DNA damage response and DNA repair genes. (C) Heatmap of expression of genes associated with iron deficiency (Puig et al., 2005Veatch et al., 2009). Data are from four independent experiments. (D) Sorted bins from WT cells were grown in a commonly used antifungal drug fluconazole and were observed under microscope for growth over 7 days. The images show growth of cells in bins HI, M1, M2 and LO in 50 μg/ml of fluconazole after 7 days. (E) Cells of HI bin showed significantly higher survival compared to other bins in both 50 μg/ml (three independent experiments) and 60 μg/ml fluconazole (four independent experiments). Cells were grown in liquid medium supplemented with fluconazole on microscopy plates and viability was calculated from microscopic observations over 7 days. Colonies showing growth rate above 0.02 h−1 after first time point were considered to be survivors. Error bars show ±1 s.d. values. (F) Percentage survival of high and low TMRE cells on fluconazole plates. High TMRE cells showed higher survival than low TMRE cells (Mann-Whitney U test). A substantial fraction of surviving high TMRE cells were respiration competent. The error bars represent ±1 s.d. values from six technical replicates for each bin. X-axis shows fluconazole concentrations used from three independent experiments. (G) From RNA sequencing data, cells from HI bin showed significantly higher expression of multidrug transporter PDR5 gene and its transcriptional activator PDR3 compared to cells from bins M1, M2 and LO. Results are from four independent experiments.

https://doi.org/10.7554/eLife.38904.015

In bacteria (Ihssen and Egli, 2004) and yeast (Levy et al., 2012; Brauer et al., 2008; Lu et al., 2009), slow growing cells can have increased stress resistance. We therefore tested whether the cells with high TMRE in a population are more resistant to acute heat stress. However, cells with high TMRE were more sensitive to heat shock as well as to oxidative stress (Figure 5—figure supplement 2A,B) and expressed some stress-response genes at lower levels (Figure 5—figure supplement 2C).

Slow growing microbes and cancer cells often have increased drug resistance (Fridman et al., 2014; Yaakov et al., 2017; Wakamoto et al., 2013; Ramirez et al., 2016; Levin-Reisman et al., 2017; Moore et al., 2012). Moreover, in several species of fungi, complete loss of mitochondria is associated with elevated tolerance to some antifungal drugs (Hallstrom and Moye-Rowley, 2000; Brun et al., 2004; Zhang and Moye-Rowley, 2001) and respiration deficient strains are often isolated from drug-treated patients (Bouchara et al., 2000; Ferrari et al., 2011; Peng et al., 2012).

We tested therefore whether the high TMRE cells differed in their sensitivity to a clinically used antifungal drug, fluconazole. We found that cells with high TMRE were ~5–7 fold more resistant to high concentrations of fluconazole (Figure 5D,E; Figure 5—figure supplement 3A). The cells surviving fluconazole treatment included a sub-fraction able to respire (Figure 5F). Cell-to-cell variation in mitochondrial membrane potential is therefore also an important predictor of cell-to-cell variation in drug resistance.

Fluconazole targets the cytochrome P450 14α-sterol demethylase enzyme (ERG11), resulting in depletion of ergosterol, a key component of the yeast cell membrane (Hitchcock et al., 1990; Ghannoum and Rice, 1999). Survival in fluconazole has been previously reported to depend on the multidrug transporter PDR5 (Brun et al., 2004; Parsons et al., 2004; Miranda et al., 2010). High TMRE cells had significantly higher level of PDR5 expression (Figure 5G; Figure 5—figure supplement 3B) as well as higher expression of the ergosterol biosynthesis pathway (Figure 5—figure supplement 3D). Consistent with previous work (Katzmann et al., 1994; Delaveau et al., 1994), the elevated expression of PDR5 in the high membrane potential cells was dependent on the PDR3 transcription factor (Figure 5—figure supplement 3C). Thus, the increased resistance in fluconazole of high TMRE cells is likely to be mediated, at least in part, by increased expression of a multidrug transporter.

Discussion

In summary, we have shown here that mitochondrial membrane potential – but not the amount of mitochondria – varies substantially across individual yeast cells and that this is associated with cell-to-cell heterogeneity in proliferation, mutation outcome, and stress and drug resistance. Laboratory strains of yeast have long been known to generate respiratory deficient ‘petite’ colonies at quite high frequency (Dimitrov et al., 2009; Nagley and Linnane, 1970; Evans et al., 1985). However, a slow growing sub-population of cells was observed in all the laboratory, natural, and clinical strains that we tested (Figure 1A,B).

Although, mitochondrial genes showed the strongest enrichment for an increased slow fraction in our gene deletion screen, other causes of slow growth will, of course, also exist. For example, deletions of genes associated with chromosome segregation and nucleus organization also affected heterogeneity but had no apparent relation to mitochondrial function. Environmental conditions, both during growth prior to microscopy, and during microscopy, are likely to affect heterogeneity; pre-growth in ethanol reduces proliferation heterogeneity (Figure 4F), presumably due to selection for cells with well-functioning mitochondria and high mtDNA copy number.

Previous theoretical studies have proposed that variability in the partitioning of cellular components could lead to heterogeneity (Birky and Skavaril, 1984; Huh and Paulsson, 2011). However, prior experimental work on the fidelity of mitochondria inheritance has shown it to be high, suggesting it is likely to be of little phenotypic consequence for single cells (Jajoo et al., 2016). In contrast, we have shown here that cell-to-cell variation in the state of the organelle can be high and predicts phenotypic variation among single cells. Although variation in mitochondrial membrane potential has been shown to exist in isogenic yeast populations (Fehrmann et al., 2013), here we have linked such variations to heterogeneity in growth and drug resistance phenotypes. Variation in mitochondrial membrane potential was related to variation in mtDNA copy number in individual cells, but we do not currently know if this is the only – or even the most common – cause of variation in the state of the organelle across single cells. Future work will be required to track down the upstream, proximal causes of this cell-to-cell variation in organelle functional state. The list of gene deletions that alter growth heterogeneity that we have reported here provide a rich resource for this future work.

In addition to the TMRE method we used, three different methods have now been used to identify slow-proliferating yeast within a population: GFP tagged TSL1, GFP tagged HSP12, and FitFlow, a direct method (Levy et al., 2012; Li et al., 2018; van Dijk et al., 2015). These methods agree on some aspects and disagree on others. For example, early work on proliferation heterogeneity in yeast has correlated expression noise in a stress response gene with proliferation heterogeneity (Levy et al., 2012) with this ultimately linked to the variation in the activity of the Ras/cAMP/protein kinase A (PKA) pathway and its target transcription factors MSN2 and MSN4 in individual yeast cells leading to variation in expression of downstream stress response genes (Li et al., 2018). Some of our results are identical to that of Li et al. (Li et al., 2018); we also find that pde2 and ira2 cells have a reduced fraction of slow growing microcolonies. Others differ. Importantly, TSL1 mRNA expression is lower in the high TMRE sub-population, as are the MSN2/MSN4 targets CTT1, DDR2 and HSP12.

Other studies have suggested a role for the DNA damage response in the generation and or maintenance of proliferation heterogeneity (Yaakov et al., 2017; van Dijk et al., 2015). Results from our screen show that the slow growing sub-populations in an isogenic yeast population are primarily driven by loss of mitochondrial function in lab, natural and clinically-isolated yeast strains. Loss of mitochondrial function leads to a loss of respiration capability and induced DNA damage response, consistent with previous reports (Yaakov et al., 2017; van Dijk et al., 2015). However, while it is tempting to speculate that the mitochondria defects in the high TMRE sub-population cause the DNA damage seen in other studies, the partial conflict in the mRNA-sequencing data suggest that this is not the case.

Taken together, results across all four studies suggest that there are multiple types of slow-proliferating cells. It will be interesting to determine the full extent of the heterogeneity in cell-state and stress resistance within isogenic populations.

Prior work on the causes of variation in proliferation rates, stress and drug resistance, and mutation outcome across individuals and individual cells has focused on fluctuations in gene expression as causative influences (Levy et al., 2012; Yaakov et al., 2017; Shaffer et al., 2017; Burga et al., 2011; Raj et al., 2010; Eldar et al., 2009; Li et al., 2018; Rotem et al., 2010; Battich et al., 2015). Here we have shown that, in yeast, variation in an organelle is strongly associated with heterogeneity in gene expression across single cells. In animals, inherited and somatic genetic variation in the mitochondrial genome can act as an important modifier of phenotypic variation (Haag-Liautard et al., 2008; Kujoth et al., 2007; Tyynismaa and Suomalainen, 2009; Kauppila et al., 2016). Recent work has also revealed substantial variation in mtDNA copy number across human tumors (Reznik et al., 2016). Moreover, in mammalian cells, mitochondrial variability has been suggested to be an important influence on cell-to-cell variation in gene expression and splicing (Johnston et al., 2012; Guantes et al., 2015; Guantes et al., 2016) and to influence variability in cell death by modulating apoptotic gene expression (Márquez-Jurado et al., 2018). Taken together, these results suggest important roles for cellular organelles, in general, and mitochondria, in particular, in the generation of heterogeneity among individual cells. In future work, therefore, it will be important to test the extent to which cell-to-cell variation in the state of mitochondria and other organelles also contributes to variable phenotypic outcomes, mutation effects, and drug resistance in human cells, including in cancer.

Materials and methods

High-throughput microscopy assay

Our high-throughput microscopy assay was inspired by the microcolony growth measurement assay by Levy et al. (2012). 96 strains were grown from glycerol stocks in a 96-well plate containing Synthetic Complete medium (0.67% Yeast Nitrogen Base without amino acids and 0.079% Complete Synthetic Supplement (ForMedium, UK)) with 2% glucose (SCD) for 24 hr at 30°C. The cells were diluted 1:50 in fresh medium, grown for 20 hr and diluted again 1:50 in fresh medium. Finally, cells were grown for 4 hr, cell densities were determined by OD at 600 nm in a Tecan plate reader and then were diluted to another plate containing SCD or appropriate medium required for microscopy experiment using a Biomek NX (Beckman Coulter) liquid handling robot, capable of pipetting variable volumes of cells across wells in a 96 plate, to a target density of ~17000 cells/μl. This minimized any possible bias due to variability in cell densities among strains. A final 5-fold dilution was done by pipetting 80 μl cells onto a pre-coated 96-well microscopy plate containing 320 μl of SCD. The microscopy plate was then sealed with LightCycler 480 sealing foils (Roche), cells were spun at 450 rpm for 2 min and taken for microscopy observations.

Microscopy plates (96-well glass bottom, MGB096-1-2-LG-L, Brooks Life Science Systems) were coated with 200 μl sterile solution of 200 μg/ml concanavalin A (type IV, Sigma) at 37°C for 16–18 hr. The solutions were then pipetted out and the plates were washed twice with sterile milli-q water. Plates were dried at 4°C for at least 24 hr. Imaging was performed using an ImageXpress Micro (Molecular Devices) microscope, with laser autofocusing, at an interval of 90 min for up to 12 hr. The microscope chamber was maintained at 30°C.

Image processing

Images were processed using custom scripts written in perl. Yeast cells were identified by juxtaposition of bright and dark pixels (10,36). A pixel was considered ‘bright’ if its intensity exceeded mean +2.2 s.d. value and a pixel was considered ‘dark’ if its intensity was below mean-2.2 s.d. value. In addition, Sobel’s edge detection algorithm (Sobel and Feldman, 1968) was applied for identifying yeast cell boundaries with sharp changes in pixel intensity. Clustering was used to identify the microcolonies and the centroid position for each microcolony was calculated. Microcolonies were tracked through centroid tracking over time. Sudden increase or decrease in centroid number in a time series indicated a failure in image acquisition or image processing and such images were discarded from the analysis. To differentiate cells from cellular debris, residuals of concanavalin A coating etc., two filtering steps were used. First, only objects that were bigger than 50 pixels at the start of observation were considered. Second, the object had to increase its size to greater than 2-fold at the end of observation. Whether neighbouring colonies touch each other at any point in time during microscopy observation was also checked. If they did, they were tracked only up to the time they touched each other. Similar to the findings of Ziv et al. (Ziv et al., 2013), lag phase during growth was observed in some microcolonies in our experiments. In addition, growth slowdown near the end of observation for some microcolonies was also observed, possibly due to nutrient limitation (Figure 5—figure supplement 4).

To calculate growth rate, linear regression on natural log-transformed area vs. time for three consecutive time points was performed and only fits with R2 ≥0.9 were considered. This was repeated using a three-point moving window over all time points. Of all these regressions, the maximum value was chosen as the microcolony growth rate to avoid biases because of slow down of growth during lag and/or due to possible substrate limitation near the end of observation.

Screening of deletion mutants, classification and functional enrichment analysis

Growth distributions for deletion mutants were measured in three independent experiments on different days. To calculate reproducibility of growth rate between replicates, mean growth rate was allowed to vary up to 0.05 h−1 and then Kolmogorov-Smirnov distance (K-S distance) (Justel et al., 1997) was calculated between all replicates after shifting one of them (through addition/subtraction) by difference in mean growth rates. Three replicates were considered as three nodes in a graph with K-S distance between them as the edge weight. If the K-S distance between two replicates exceeded 0.1, no edge was drawn between those two nodes. The sub-graph where the maximum number of nodes was connected to each other directly and via shortest possible distance was considered as reproducible replicates. Only mutants with at least two reproducible replicates were considered in our analysis. The number of reproducible replicates for each mutant is given in Supplementary file 2. The proliferation distributions for all mutants are shown in Supplementary file 4.

To calculate slow fraction from a proliferation distribution, first, a cumulative distribution function (cdf) was calculated with density being calculated at an interval of 0.01 h−1. The cdf function was then scanned for maximum slope using a window of 5 points. At the point with maximum slope, a line with the maximum slope was fitted and the points that deviated from the fitted line by >0.02 h−1 were considered as the edges of the main subpopulation. In the next step, if the left sub-population was bigger than the right sub-population, the percentage of slow fraction was calculated as (% left sub-population-% right sub-population) and the percentage of fast sub-population was set to zero. If the right sub-population was bigger than the left sub-population, the percentage of fast fraction was calculated as (% right sub-population-% left sub-population) and the percentage of slow sub-population was set to zero.

If the mode of a growth distribution is reduced (compared to WT) and the growth rate of the slow fraction is not reduced, the main sub-population growth distribution is likely to overlap with and mask a slow growing sub-population. To avoid such scenarios, all the reproducible growth distributions for the WT strain were collected and the mode growth rate was computationally reduced in steps of 0.01 h−1 without changing the growth rate of the slow sub-population. The percentage slow fraction was calculated at each step. As expected, reduction in mode growth rate without moving the slow fraction led to a reduction in % of slow fraction (Figure 2—figure supplement 2A).

Mutants with altered mode proliferation rate compared to WT strain were identified through Mann-Whitney U test (FDR < 0.1). Mutants with altered slow fraction were identified by Mann-Whitney U test (FDR < 0.1) after correcting for any change in mode growth rate (Figure 2—figure supplement 2A). Comparison between replicate measurements of mode growth rate and percentage of slow fraction was done (Figure 5—figure supplement 5A). The mean proliferation rate of mutant strains obtained in our assay was comparable with published values (Figure 5—figure supplement 5B).

We used GOslim gene annotation (Gene Ontology Consortium, 2015Gene Ontology Consortium, 2018) for functional class enrichment analysis and we performed a hypergeometric test as follows. Let us assume that in a group ‘g’ from screening (for example, the group with increased slow fraction but no change in mode growth rate), out of total Ng genes, Xg genes are associated with function f according to GOslim annotation. Let us also assume that out of total N genes screened in our data, X genes belong to the functional class f according to GOslim annotation. Thus, the probability that the group ‘g’ contains more number of genes of functional class f than expected by chance alone is given by p=i=XgNg(Xi)(NXNgi)(NNg), which gives the p-value. A further multiple testing correction was done using Benjamini-Hochberg procedure with FDR<0.1.

Quantification of incomplete penetrance

Incomplete penetrance was calculated for all mutants that showed significant reduction in mean proliferation rate compared to the WT strain (Mann-Whitney U test, FDR < 0.1). For each of these mutants, replicate proliferation distributions were compared with replicate proliferation distributions of WT strain that were reproducible across the screening experiment. Average overlap of the mutant proliferation distributions with the bulk sub-population of each of the WT proliferation distribution was calculated. For WT strain, to calculate bulk sub-population from a proliferation distribution, first, a cumulative distribution function (cdf) was calculated with density being calculated at an interval of 0.01 h−1. The cdf function was then scanned for maximum slope using a window of 5 points. At the point with maximum slope, a line with the maximum slope was fitted and the point that deviated from the fitted line by >−0.02 h−1 was considered as the edge of the bulk sub-population. Thus, for each mutant, this led to a distribution of a percentage of cells showing WT-like proliferation (Figure 2—figure supplement 1). In the next step, it was tested whether the distribution of percentage of WT-like cells was significantly different from zero (Wilcoxon rank-sum test for one sample) and an FDR correction for multiple testing was performed (FDR < 0.1).

Mitotracker green and TMRE staining

To perform mitotracker green (MitoTracker Green FM, Molecular Probes, Thermo Fisher Scientific) staining, cells were centrifuged at maximum speed for 2 min and washed twice with buffer containing 10 mM HEPES (pH 7.4) and 5% glucose. Cells were then re-suspended in the same buffer and Mitotracker Green (10 μM stock dissolved in DMSO) was added to a final conc. of 100 nM. Cells were incubated for 20 mins at 30°C, washed twice with PBS (pH 7.4) and quantified by flow cytometry (LSR Fortessa, BD Biosciences).

TMRE (Tetramethylrhodamine, Ethyl Ester, Perchlorate) (Molecular Probes, Thermo Fisher Scientific) is a positively charged dye that accumulates inside mitochondria depending on the mitochondria transmembrane potential generated due to transfer of protons across mitochondrial membrane resulting in net negative charge inside the mitochondria (Crowley et al., 2016). For TMRE staining, cells were grown as in pre-growth step in microcolony assay, precipitated, washed twice with PBS, were re-suspended in PBS, and TMRE was added to a final conc. of 100 nM from a 10 mM stock dissolved in DMSO. Cells were incubated at 30°C for 30 min, were washed twice with PBS and were analysed by flow cytometry or were sorted. There was a gap of 15–20 min between the end of staining and beginning of flow cytometry experiments due to the time required for cleaning, priming and setting up of flow cytometry machine parameters. Day-to-day variations were observed in measurement of TMRE distributions.

Cell sorting and growth measurement of sorted bins

Cells were sorted by TMRE signal into four bins HI, M1, M2, LO (Supplementary file 5, Figure 3C) in an Aria II SORP cell sorter (BD Biosciences). For growth rate measurement, stress resistance measurement, and mitochondrial DNA quantitation by qPCR in the sorted bins, 100,000 cells per bin were sorted at room temperature into 1.5 ml tubes pre-filled with 600 μl of PBS. After sorting, 200 μl of YPD was added to each tube, cells were centrifuged for 5 min at maximum speed at room temperature, and the supernatant was thrown away. Cells were re-suspended in 600 μl of PBS before proceeding for subsequent experiments. For heat shock experiments, 100 μl of sorted cells were put into PCR tubes and were subjected to heat shock in a PCR machine, put on ice for 1 min before measurement of growth and viability. For RNA sequencing experiments, 750,000 cells per bin were sorted in three 1.5 ml tubes, each pre-filled with 800 μl of PBS. After sorting, 200 μl of YPD was added to each tube, centrifuged at maximum speed for 5 min, and supernatants were discarded. The cell pellets were gently washed twice using PBS. Total RNA was isolated using MasterPure yeast RNA isolation kit (Epicentre) following manufacturer’s protocol. Cells from four sorted bins were grown in SCD medium up to 48 hr and their growth distributions were remeasured using high-throughput microscopy assay (Figure 5—figure supplement 6A).

To determine percentage of respiration deficient cells, cells were plated on plates containing Synthetic Complete medium with 3% glycerol and 0.1% glucose (SCDG) solidified with 1.5% agar to a target density of ~100–150 colonies per plate. After 5–7 days, number of small and big colonies were counted and the percentage of respiration deficient cells were determined as - percentage of respiration deficient cells = Number of small colonies(Number of small+big colonies) ×100

Respiration deficient cells showed almost no mtDNA and remained slow growing even after seven days of growth in Synthetic Complete medium with 2% glucose (SCD) (Figure 5—figure supplement 6B,C).

Cells were also tested for switching of respiration capability (Figure 5—figure supplement 6D,E). Equal number of sorted cells from each bin was plated onto SCDG plates (with 0.1% glucose and 3% glycerol as the carbon sources) and SCG plates (with 3% glycerol as the carbon source) and allowed to grow for 5 days. Percentage of cells regaining respiration capability was calculated as (No. of respiration capable cells on SCDG plate – No. of colonies on SCG plate)/Total no. of colonies on SCDG plate.

Growth rate switching in microcolonies

To test whether microcolonies switch from fast to slow or slow to fast growth rate, microcolony growth rates at all time points of tracking were calculated. Growth rate of a microcolony at a time point was calculated using linear regression of ln(area) with time including one preceding time point and one subsequent time point and only fits with R2 ≥0.9 were considered. In case of two cutoff values (say ‘c1’ as the lower cutoff and ‘c2’ as the higher cutoff) for growth rate for determining slow and fast growing microcolonies, if a microcolony showed growth rate above ‘c2’ for at least three consecutive time points and afterwards showed growth rate below ‘c1’ in at least three consecutive time points and the time of switching from fast to slow growth was not within last three time points of observation (to avoid incorrect classification due to slowdown in growth at later time points), the colony was classified to be switching from fast to slow growth (Figure 5—figure supplements 7, 8 and 9). Similarly, if a microcolony showed growth rate below ‘c1’ for at least three consecutive time points and afterwards showed growth rate above ‘c2’ in at least three consecutive time points and the time of switching from slow to fast growth was not within first 3 time points of observation (to avoid incorrect classification due to the lag phase), the colony was classified to be switching from slow to fast growth (Figure 5—figure supplements 7, 8 and 9). In case of a single cutoff value, the same criteria were applied as above with ‘c1’ being equal to ‘c2’. This was applied to identify switching in cells of TMRE sorted bins HI, M1, M2 and LO in the WT strain. Various cut-off values for identifying switching were tested to check switching of microcolony growth rates across different ranges of growth rates (Figure 5—figure supplements 7, 8 and 9).

For calculating switching in unsorted WT and deletion strains, the criteria for classification were slightly modified since there were fewer time points of observation. Specifically, if a microcolony showed growth rate above ‘c2’ for at least two consecutive time points and afterwards showed growth rate below ‘c1’ in at least two consecutive time points and the time of switching from fast to slow growth was not within last 2 time points of observation (to avoid incorrect classification due to slowdown in growth at later time points), the colony was classified to be switching from fast to slow growth (Figure 5—figure supplement 7). Similarly, if a microcolony showed growth rate below ‘c1’ for at least two consecutive time points and afterwards showed growth rate above ‘c2’ in at least two consecutive time points and the time of switching from slow to fast growth was not within first 2 time points of observation (to avoid incorrect classification due to lag), the colony was classified to be switching from slow to fast growth (Figure 5—figure supplement 7).

Switching of cells from high to low membrane potential

To test for switching between high and low TMRE states, cells with high, medium or low TMRE were grown for 48 hr after sorting and their mitochondrial membrane potential values were remeasured (Figure 5—figure supplement 10). Cells from the HI TMRE bin consisted of 99% of high TMRE cells and 1% low TMRE cells (impurity). As measured by time-lapse microscopy (Figure 5—figure supplement 11), high TMRE cells in the HI bin consisted of slow-growing (80%) and fast growing (20%) cells. Low TMRE cells in the HI bin (impurity) were assumed to be fast growing, as this gave a conservative estimate for the number of cells switching from high to low TMRE state in 24 hr. The mean growth rates of the slow and fast-growing cells (geometric mean) were 0.20 h−1 and 0.37 h−1 respectively (Figure 5—figure supplement 11). The percentages of cells switching from low TMRE to high TMRE state in 24 hr was estimated as the average switching rate of cells from M1, M2 and LO bins to high TMRE state and it was estimated to be ~10% in 24 hr (Figure 5—figure supplement 10). Given these data, assuming exponential growth for all sub-populations over 24 hr and in the absence of switching from high to low TMRE state,

%oflowTMREcells=0.01×e0.37×240.1×(0.01×e0.37×24)0.80×e0.20×24+0.19×e0.37×24+0.01×e0.37×24×100% =71.877.1997.20+1365.49+71.87× 100% =4.21%

This estimate is 8-fold lower than the observed % of low TMRE cells (34.7%).

For growth of cells from 24 hr to 48 hr, assuming exponential growth for 24 hr and no switching from high to low TMRE state,

%oflowTMREcellsafter48hr=0.347×e0.37×240.1×(0.347×e0.37×24)0.80×0.653×e0.20×24+0.20×0.653×e0.37×24+0.347×e0.37×24×100%  =2493.81249.463.47+1437.36+2493.81× 100% =56.2%

Again, this estimate is substantially lower than the observed % of low TMRE cells (76.5%). Taken together, these results suggest considerable switching from high to low TMRE state.

Measurement of mtDNA copy number by qPCR

To determine mtDNA copy number per cell using quantitative PCR (qPCR), five primer pairs specific to nuclear DNA (ACT1, ALG9, KRE11, TAF10, COX9) and five primer pairs specific to mitochondrial DNA (COX1, ATP6, COX3, ATP9, tRNA – primer picked around tQ(UUG)Q gene) were used (Supplementary file 3). A standard curve for each of primer was made, using six concentrations of genomic DNA serially diluted from the highest concentration by 4-fold at each step. Absolute quantification of DNA copy number was performed using the standard curve. Three technical replicates for each primer and for each sample were set up totaling 30 reactions per sample. To compare mtDNA copy number across sorted bins, nuclear DNA and mtDNA copy numbers in all bins were normalized by the respective values for LO bin. Two sample t-test was used to check whether the normalized value for nuclear DNA differs significantly from the normalized value for mtDNA and a p-value was calculated using a two-sample t-test. Mean mtDNA copy number per cell was calculated by the ratio of mtDNA to nuclear DNA and standard deviations were calculated by taking error propagation models into account.

To overexpress the MIP1 gene, the MIP1 gene under the control of the native promoter (930 bp upstream and 262 bp downstream, total insert length - 4957 bp) was cloned into pRS413 plasmid and then transformed into NEB 10B electrocompetent E. coli cells. The plasmid with the verified construct was then isolated and transformed into yeast cells.

RNA sequencing experiment and data analysis

Isolated total RNA (using MasterPure yeast RNA isolation kit (Epicentre)) was checked and quantified using bioanalyzer. 200 ng of total RNA for each sample was taken and was mixed with 4 μl of 1:1000 dilution of ERCC spike-in mix1 (Thermo Fisher Scientific). Sequencing was done in Illumina HiSeq with paired end 2 × 50 bp reads. Quality of the sequenced reads was checked using FastQC (Andrews, 2016) and then the reads were mapped to reference yeast transcriptome (R64-1-1 reference cdna sequence from Ensembl [Hubbard et al., 2002]) using bowtie2 (Langmead and Salzberg, 2012). Mapping statistics was calculated using a custom script where only read pairs mapping concordantly and uniquely to the reference sequence were considered. The data were normalized using ERCC spike-in reads as controls using RUVg method in R package RUVSeq (Risso et al., 2014). Correlation between replicates were checked through distance heatmap and PCA analysis (Figure 5—figure supplement 12A,B), using R package DESeq2 (Love et al., 2014). Differentially expressed genes were identified using package DESeq2. Functional enrichment analysis on sets of differentially expressed genes was done using a hypergeometric test as described above with multiple testing correction (FDR < 0.1) (Benjamini-Hochberg method) with GOslim gene annotations.

Reconstruction of single mutants

Gene deletion mutants were remade in the WT strain using sequence specific homologous recombination. First, the deletion cassette from the appropriate deletion strain from the collection was amplified using primers such that the amplified region contained the deletion cassette with KanMX marker and 50–300 bp of overhang on either side of the cassette. Particular care was taken to avoid neighbouring genes from being amplified. The PCR product was transformed into competent yeast cells (prepared using lithium acetate and PLI – made by mixing 1 ml water, 1 ml 1M lithium acetate and 8 ml 50% PEG3350) and colonies were selected on G418 plates. Two verified clones for each mutant were picked for experiments. Some of the mutants associated with mitochondrial function were found to be compensated in the deletion collection (Figure 5—figure supplement 12C,D). Beyond the initial high-throughput measurement of proliferation distribution of deletion mutants, all experiments were performed with freshly made deletion mutants.

Long-term microscopy-based growth measurements for measurement of stress tolerance and drug resistance

To observe growth of yeast cells in drug (fluconazole dissolved in DMSO, stock conc. 5 mg/ml) over 7 days, yeast cells were imaged under the microscope every ~24 hr. To have a bigger part of a well imaged and to increase the number of data points, 48 fields of view per well were imaged in these experiments. Yeast cells and microcolonies were identified as above. Before tracking the microcolonies over time, the images for all fields of view in a well were merged which allowed tracking of microcolonies even if the plate was positioned slightly differently in the microscope at different time points. Microcolonies were tracked over time and growth rates were calculated. A growth rate of 0.02 h−1 after the first time point was taken as cut-off for survival on fluconazole, as most colonies showed initial growth but then stopped growing. Percentage survival in heat shock and hydrogen peroxide treatment was calculated as the ratio of the number of colonies showing growth under stressed condition compared to the total number of colonies showing growth under unstressed condition.

Measurement of respiration capability in drug resistant cells

To test whether the cells that survive fluconazole treatment can still respire, the drug resistance of the sorted sub-populations from the HI and LO bins were measured on agar plates after 15 days of growth in SCD medium supplemented with fluconazole (9.5 or 10 µg/ml) and solidified with 1.5% agar. This assay needed lower concentrations of fluconazole compared to the microscopy-based assay, as only the colonies that divided multiple times were visible on the plate. Sorted cells from bins HI and LO were plated directly after sorting onto the drug plates (5–6 replicates per bin), onto plates without any drug as well as onto SCDG plates to calculate the percentage of cells capable of respiration. Cells were counted after 15 days and 40–50 colonies from each plate were randomly picked and checked for respiration capability by plating onto plates containing 3% glycerol as the carbon source.

Data availability

RNA-sequencing data that support the findings of this study have been deposited in NCBI GEO with the accession code GSE104343. Microscopy images have been submitted to openmicroscopy.org. The raw microcolony growth data for the WT and mutant strains are available at https://github.com/lehner-lab/MicroscopyCode-Dhar_et_al/tree/master/Microscopy_screen_processed_data.

Code availability

Custom codes for analysing microscopy images are available at https://github.com/lehner-lab/MicroscopyCode-Dhar_et_al (Dhar and Faure, 2019; copy archived at https://github.com/elifesciences-publications/MicroscopyCode-Dhar_et_al). 

References

  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
    Elevated levels of petite formation in strains of saccharomyces cerevisiae restored to respiratory competence. I. Association of both high and moderate frequencies of petite mutant formation with the presence of aberrant mitochondrial DNA
    1. RJ Evans
    2. KM Oakley
    3. GD Clark-Walker
    (1985)
    Genetics 111:389–402.
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
    GO subset guide
    1. Gene Ontology Consortium
    (2018)
    Accessed December 2, 2014.
  29. 29
    A nuclear mutant of Saccharomyces cerevisiae deficient in mitochondrial DNA replication and polymerase activity
    1. A Genga
    2. L Bianchi
    3. F Foury
    (1986)
    The Journal of Biological Chemistry 261:9328–9332.
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
    Isolation and genetic characterization of cell-lineage mutants of the nematode caenorhabditis elegans
    1. HR Horvitz
    2. JE Sulston
    (1980)
    Genetics 96:435–454.
  41. 41
  42. 42
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47
  48. 48
  49. 49
  50. 50
  51. 51
  52. 52
  53. 53
  54. 54
  55. 55
  56. 56
  57. 57
  58. 58
  59. 59
  60. 60
  61. 61
  62. 62
  63. 63
  64. 64
  65. 65
  66. 66
  67. 67
  68. 68
  69. 69
  70. 70
  71. 71
  72. 72
  73. 73
  74. 74
  75. 75
  76. 76
  77. 77
  78. 78
  79. 79
  80. 80
    A 3x3 isotropic gradient operator for image processing
    1. I Sobel
    2. G Feldman
    (1968)
    Stanford Artificial Intelligence Project.
  81. 81
  82. 82
  83. 83
  84. 84
  85. 85
  86. 86
  87. 87
  88. 88
  89. 89

Decision letter

  1. Naama Barkai
    Senior and Reviewing Editor; Weizmann Institute of Science, Israel

In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.

Thank you for submitting your article "Single cell functional genomics reveals the importance of mitochondria in cell-to-cell phenotypic variation" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by Naama Barkai as the Senior and Reviewing Editor. The reviewers have opted to remain anonymous.

Please address all comments in the individual reviews below. In particular, the reviewers were concerned that the experiments showing transition from 'slow' to 'fast' state are important but are not very convincing and should be addressed. in this context, please note the paper from Gilles Charvin (Fehrmann et al., 2013). They used microfluidics to look at mitochondrial membrane potential in individual yeast cells as they aged. This paper should have been cited in the current manuscript, and it sets the bar both technically (for the microfluidic analysis of membrane potential and growth in individual cells) and conceptually (they found that decline of membrane potential is age-independent and leads to clonal senescence – i.e., is not reversible). The phenomenon studied in the current manuscript is different because it has to do with high-potential/low-respiration cells, but the Charvin study sets the precedent in the literature at least for irreversible mitochondrial-state changes.

Reviewer #1:

Dhar et al. address an important topic of broad interest: cellular heterogeneity and its relationship to proliferation and survival. Given that very little is known about the causes of heterogeneity, their approach of doing a genetic screen for yeast mutants that alter the distribution of growth rates within a population is sensible, and the scale of the experiment (>1500 mutants) is impressive. However, despite my enthusiasm for the research question and the general approach, I do have some major concerns:

1) The authors appear to have reinvented the micro-colony growth assay of Levy et al. (Levy et al., 2010) and should do a better job of presenting the similarities and differences between the two methods. Of particular concern is that they appear to have cell behaviors that were not seen by the other group – they mention lag (subsection “Image processing” and subsection “Cell sorting and growth measurement of sorted bins”, end of last paragraph) and nutrient deprivation, neither of which was seen by the other group (Levy et al., 2010 and Ziv et al., 2013) in medium with sufficient glucose.

2) The finding of a role for mitochondria in growth heterogeneity is interesting, but petites are well known to decrease growth rate and increase stress resistance. So the key advance here would be the finding that mitochondrial states switch back and forth, rather than petites just being a slow-growing dead end. Levy et al., 2010, had already demonstrated switching of growth states, and that single micro-colonies could give rise to the complete distribution of growth rates. So the question is whether the switching is driven by mitochondrial-state changes. Unfortunately, the experiments on switching are not entirely convincing. The authors took bins of cells sorted by TMRE (mitochondrial membrane potential) and assayed their growth. They found that high-TMRE cells tended to have slower growth rates and that high TMRE correlated with respiration deficiency. These results would be consistent with switching only in the direction of becoming petite. But the authors argue that switching can go in both directions because sorted bins could repopulate the full distributions of TMRE staining and growth. However, noise in TMRE staining or other technical reasons could mean that the sorted bins were contaminated to some extent with respiration-proficient cells, which then might be responsible for the repopulation of the distributions. A more direct tracking experiment is necessary.

3) The manuscript does not offer much mechanistic insight into the changes in mitochondrial state or mtDNA copy number. What is causing these changes? Why do the changes persist for the number of cell generations that they do? Is the mitochondrial-state phenomenon tied to, or independent from, the heterogeneity in gene expression seen in Levy et al., 2010?

Reviewer #2:

In this manuscript, Dhar et al. studied heterogeneity of growth between yeast cells by measuring single-colony growth rates of a large collection of yeast strains (gene deletions and natural strains). The originality was to track over time the size of colonies that emerge from single cells, deriving individual growth rates of these initial founder cells. The distribution of single-cell growth rates differed between strains and the authors found marked alterations of these distribution for mutants of nuclear genes involved in mitochondrial functions. They sorted cells based on TMRE staining, and observed that high signal correlated with i) slow growth ii) low mtDNA iii) low competence for respiration iv) peculiar transcriptomic signatures v) sensitivity to heat and oxidative stress but vi) higher survival to fluconazole. Overexpression of MIP1 increased mtDNA content and reduced both the fraction of slow growing cells and the fraction of respiratory-deficient cells.

General Assessment:

This is a very sound study. The amount of work is impressive and the complementarity of the experimental approaches is remarkable. Conclusions that heterogeneity in mitochondrial 'state' affects phenotypic heterogeneity between cells is very important and opens novel avenues for investigating cell-cell biological variation, which is a hot topic.

1) The raw data of colony growth assays is not made available. This is not acceptable, as others may later want to explore the dataset in various ways. The raw data should be uploaded in a public repository in two forms:

- Time-course data of the size of each colony, for all colonies and all yeast strains. This is important if others want to explore time-series for lag phase, saturations or other traits that may not be fully correlated to the growth rates computed by the authors.

- The scalar values of growth-rates of individual colonies, so that others can use these values directly.

2) GO analysis (Figure 2E) highlighted respiration based on many genes but little enrichment. Other terms were better enriched, although among fewer genes, especially for the 'red' set. The authors should mention in text what are the kinases, cytoskeletal binding proteins and carbohydrate transporters involved in maintaining the% slow fraction.

Reviewer #3:

The described research aims to identify the underlying mechanisms of investigating proliferation heterogeneity in isogenic cell populations. Through a quantitative microcolony assay, the authors demonstrate in isogenic cultures of yeast maintain heterogenous and stable rates of cell division that are influenced by genetic variation. By screening the yeast deletion collection, the authors discover that mitochondrial activity influences proliferation heterogeneity, and that mitochondrial membrane potential and respiratory dysfunction are features of the "slow fraction" of these cell populations. This work highlights the importance of mitochondrial dynamics as a master regulator of cells and how these dynamics may influence heterogeneity.

S. cerevisiae is known to throw off petite cells that lack complete mtDNAs, in a genetic background dependent manner (these petites can also propagate in liquid cultures so multiple serial transfers can alter the initial number of petites in the inoculums which may influence the overall observed heterogeneity differences between strains). The "slow fraction" could represent these non-respiring petite cells, however the authors provide evidence that the "fast" to "slow" transitions are reversible (i.e. that the High TMRE, respiratory deficient cells could give rise to Low TMRE cells (subsection “Mitochondria state but not content predicts slow growth”, last paragraph, Figure 5—figure supplement 2Figure). Have these "switched" cells regained respiratory capability? To decide if cells have switched between slow and fast growth, the authors rely on a threshold value of growth rates (as opposed to a percent change). It's not clear why the 0.3 h-1 threshold was chosen. As I understood, a microcolony that changes from 0.1 to 0.29 h-1 would be scored as slow growing while a microcolony that changes from 0.29 to 0.31 h-1 would be scored as one that switched between slow and fast growth. Can the authors clarify?

• Yeast cells growing in microplates are likely to be oxygen starved (even if being shaken) and near transition states in metabolism, so it makes sense that changes in mitochondrial "states" are being observed. A discussion of how environmental conditions could influence heterogeneity would help to put this into context.

[Editors' note: further revisions were requested prior to acceptance, as described below.]

Thank you for submitting your article "Single cell functional genomics reveals the importance of mitochondria in cell-to-cell phenotypic variation" for consideration by eLife. Your revised article has been reviewed by two peer reviewers, and the evaluation has been overseen by Naama Barkai as the Senior and Reviewing Editor. The reviewers have opted to remain anonymous.

The reviewers remained supportive of your study and appreciated most of the revisions. However, the evidences supporting the claim that cells switch from a low-growth to a fast-growing mode remained inconclusive. Since this is not the essence of your results, we are happy to accept the paper. however, we ask that you will back off this claim and move Figure 5 to the supplementary material.

The reviewers also felt that the discussion of previous findings, and their relation to your work, will improve the manuscript. please try to do that as well. see below specific comments.

Reviewer #1:

In this revised manuscript, Dhar et al. present new analyses aiming to address concerns about their inferences of switching from a slow-growing, respiration-deficient cell state to a fast-growing, respiration-capable cell state. I remain enthusiastic about the topic of the manuscript and I remain impressed by the scale of the mutant analysis, but I also remain concerned that the authors' analyses of switching rely on indirect arguments that are not entirely convincing. Specifically:

1) They argue (subsection “Switching of cells from high to low membrane potential”) that TMRE-high to TMRE-low switching must take place to explain the observation that 1% of cells are TMRE-low immediately after sorting into the TMRE-high bin (a low rate of contamination), whereas 34.7% of cells are TMRE-low 24 hours later. Their argument relies on a simple model in which TMRE-high cells consist of two subpopulations (80% slow-growing and 20% fast-growing), TMRE-low cells are fast-growing, and the average slow and fast growth rates represent the corresponding subpopulations. This model has some problems. First, the arithmetic mean growth rate might not be appropriate to model changes in subpopulation sizes over time. If each cell chooses a growth rate at random from its distribution then what matters is the geometric mean, which is heavily influenced by low values. Alternatively, if growth rate is relatively stably inherited through time then what matters is the fastest members of a subpopulation. Either way, arithmetic mean could be a very poor indicator of subpopulation growth. Second, even if we ignore the problem with the arithmetic mean and take the authors' argument and apply it to the transition between the 24-hr time point and the 48-hr time point, then their calculation actually predicts quite well the 48-hr subpopulation sizes from the 24-hr starting point, so TMRE-high to TMRE-low switching is not required to explain this transition. Third, the authors' simple model ignores lag, which would likely further under-represent the slow-growing (high-TMRE) fraction, and might actually be why the transition from 0 to 24 hours requires additional explanation whereas the transition from 24 to 48 hours does not.

2) The authors have verified that their analyses of growth rate switching in micro-colonies are robust to changes in cutoffs between "fast" and "slow" growth, which is good, but what these switches mean is insufficiently discussed. In particular, fast-to-slow switches such as those shown in Figure 5D necessarily involve concerted changes in most cells in a micro-colony because a single cell changing from fast to slow in the middle of micro-colony growth would not appreciably change the growth rate of the micro-colony. The authors should discuss (and perhaps experimentally investigate) what such switches at the level of the micro-colony mean.

3) The conclusion that some cells regain ability to respire is an interesting piece of this work but it is not clear exactly what the authors mean by it. In Figure 5E they present the "percentage of colonies that regained capability to respire" but the calculation involves two different plates that were not replicas, so the ability to "regain" respiration was a more indirect inference than the text suggests.

Because of the above concerns, my view remains that more direct evidence of switching at the level of single cells (from high-TMRE to low-TMRE, from slow growth to fast growth, and from respiration-deficient to respiration capable) would substantially strengthen this paper.

Reviewer #3:

I completely agree that the resubmission only indirectly addresses our concerns about slow to fast growth. Even with the new statistical analysis, some quick growing contaminants could explain what they conclude are transitions. A simple replica plating experiment would have shown that their non-respiring cells can convert back to respiring cells. Given that they were previously asked for alternate experimentation to demonstrate the slow to fast transition, they need to back off this claim. (and minimize and move Figure 5 to supplemental).

In general, the authors could do better at discussing their findings in the context of previously published work.

Their proposal (in the Discussion) that petites are part of a mitochondrial heterogeneity continuum is not substantiated if they are unable to convincingly show reversibility.

In the Abstract, they should change the word "developed" a high-throughput method to "used".

The authors now refer to mitochondrial membrane potential in places though there are still several places where mitochondrial "state" is used and could/should be corrected (e.g. Subsection “Mitochondrial membrane potential but not amount predicts slow growth”, fourth and fifth paragraphs, subsection “Variation in mitochondrial state predicts additional phenotypic heterogeneity including drug resistance”, first paragraph and especially in the first paragraph of the Discussion).

https://doi.org/10.7554/eLife.38904.042

Author response

Please address all comments in the individual reviews below. In particular, the reviewers were concerned that the experiments showing transition from 'slow' to 'fast' state are important but are not very convincing and should be addressed. in this context, please note the paper from Gilles Charvin (Fehrmann et al., 2013). They used microfluidics to look at mitochondrial membrane potential in individual yeast cells as they aged. This paper should have been cited in the current manuscript, and it sets the bar both technically (for the microfluidic analysis of membrane potential and growth in individual cells) and conceptually (they found that decline of membrane potential is age-independent and leads to clonal senescence – i.e., is not reversible). The phenomenon studied in the current manuscript is different because it has to do with high-potential/low-respiration cells, but the Charvin study sets the precedent in the literature at least for irreversible mitochondrial-state changes.

We would like to thank the editor and the three reviewers for their constructive suggestions that helped improve our manuscript.

We have now added results of further analysis on ‘slow’ to ‘fast’ switching (Figure 5B-D and Figure 5—figure supplements 3-7). In addition, we also show transition of cells from high to low TMRE state and vice versa (Figure 5A and Figure 5—figure supplement 10). We have now cited the paper from the group of Gilles Charvin (Fehrmann et al., 2013) in the main text (subsection “Mitochondrial membrane potential but not amount predicts slow growth”, last paragraph).

We address the specific issues raised by each of the reviewers below.

Reviewer #1:

Dhar et al. address an important topic of broad interest: cellular heterogeneity and its relationship to proliferation and survival. Given that very little is known about the causes of heterogeneity, their approach of doing a genetic screen for yeast mutants that alter the distribution of growth rates within a population is sensible, and the scale of the experiment (>1500 mutants) is impressive. However, despite my enthusiasm for the research question and the general approach, I do have some major concerns:

1) The authors appear to have reinvented the micro-colony growth assay of Levy et al. (Levy et al., 2010) and should do a better job of presenting the similarities and differences between the two methods. Of particular concern is that they appear to have cell behaviors that were not seen by the other group – they mention lag (subsection “Image processing” and subsection “Cell sorting and growth measurement of sorted bins”, end of last paragraph) and nutrient deprivation, neither of which was seen by the other group (Levy et al., 2010 and Ziv et al., 2013) in medium with sufficient glucose.

We have now added more details to the subsection “Image processing” and show examples of lag and growth slowdown towards the end of the observation in Figure 5—figure supplement 4. We also note that lag was observed in the microcolony growth assay by Ziv et al. [Ziv et al., 2013].

2) The finding of a role for mitochondria in growth heterogeneity is interesting, but petites are well known to decrease growth rate and increase stress resistance. So the key advance here would be the finding that mitochondrial states switch back and forth, rather than petites just being a slow-growing dead end. Levy et al., 2010, had already demonstrated switching of growth states, and that single micro-colonies could give rise to the complete distribution of growth rates. So the question is whether the switching is driven by mitochondrial-state changes. Unfortunately, the experiments on switching are not entirely convincing. The authors took bins of cells sorted by TMRE (mitochondrial membrane potential) and assayed their growth. They found that high-TMRE cells tended to have slower growth rates and that high TMRE correlated with respiration deficiency. These results would be consistent with switching only in the direction of becoming petite. But the authors argue that switching can go in both directions because sorted bins could repopulate the full distributions of TMRE staining and growth. However, noise in TMRE staining or other technical reasons could mean that the sorted bins were contaminated to some extent with respiration-proficient cells, which then might be responsible for the repopulation of the distributions. A more direct tracking experiment is necessary.

We show that the contaminating cells in the HI bin cannot explain the percentage of low TMRE cells after 24 hours of growth (Figure 5A and Figure 5—figure supplement 10). We describe the calculations in the Materials and methods section (subsection “Switching of cells from high to low membrane potential”).

3) The manuscript does not offer much mechanistic insight into the changes in mitochondrial state or mtDNA copy number. What is causing these changes? Why do the changes persist for the number of cell generations that they do? Is the mitochondrial-state phenomenon tied to, or independent from, the heterogeneity in gene expression seen in Levy et al., 2010?

The underlying mechanism that drives changes in mitochondrial state and mtDNA copy number is indeed not clear. However, in Figure 4 we show that increased expression of the mitochondrial DNA polymerase Mip1 increased mtDNA copy number 3-fold (Figure 4B) and reduced both the fraction of slow proliferating and respiration-deficient cells and the fraction of cells with high TMRE signal (Figure 4C, D). This is consistent with defects in mtDNA replication causing the mitochondrial state and mtDNA copy number variation. The mitochondrial state phenomenon is likely to be different the from the heterogeneity in gene expression seen in Levy et al. 2010, as unlike Levy et al. 2010, we observed lower expression of TSL1 gene in our slow growing cells (Figure 6—figure supplement 2C).

Reviewer #2:

[…] 1) The raw data of colony growth assays is not made available. This is not acceptable, as others may later want to explore the dataset in various ways. The raw data should be uploaded in a public repository in two forms:

- Time-course data of the size of each colony, for all colonies and all yeast strains. This is important if others want to explore time-series for lag phase, saturations or other traits that may not be fully correlated to the growth rates computed by the authors.

- The scalar values of growth-rates of individual colonies, so that others can use these values directly.

Data is available at https://github.com/lehner-lab/MicroscopyData_Dhar_et_al

2) GO analysis (Figure 2E) highlighted respiration based on many genes but little enrichment. Other terms were better enriched, although among fewer genes, especially for the 'red' set. The authors should mention in text what are the kinases, cytoskeletal binding proteins and carbohydrate transporters involved in maintaining the% slow fraction.

We have now mentioned some of these genes in the main text (subsection “Deletion of genes involved in mitochondrial function alter heterogeneity”, last paragraph).

Reviewer #3:

[…] • S. cerevisiae is known to throw off petite cells that lack complete mtDNAs, in a genetic background dependent manner (these petites can also propagate in liquid cultures so multiple serial transfers can alter the initial number of petites in the inoculums which may influence the overall observed heterogeneity differences between strains). The "slow fraction" could represent these non-respiring petite cells, however the authors provide evidence that the "fast" to "slow" transitions are reversible (i.e. that the High TMRE, respiratory deficient cells could give rise to Low TMRE cells (subsection “Mitochondria state but not content predicts slow growth”, last paragraph, Figure 5—figure supplement 2). Have these "switched" cells regained respiratory capability? To decide if cells have switched between slow and fast growth, the authors rely on a threshold value of growth rates (as opposed to a percent change). It's not clear why the 0.3 h-1 threshold was chosen. As I understood, a microcolony that changes from 0.1 to 0.29 h-1 would be scored as slow growing while a microcolony that changes from 0.29 to 0.31 h-1 would be scored as one that switched between slow and fast growth. Can the authors clarify?

We have now addressed this issue in detail and we show that the results are robust to the cut off values that we use (Figure 5A-C and Figure 5—figure supplements 2 and 4). In addition, we have also added an analysis on bin switching where we divide the range of entire growth rate into 6 bins and calculate the percentage of cells switching from one bin to another (Figure 5—figure supplement 5). We also show that some cells lose respiration capability transiently and are able to regain this capability (Figure 5E).

• Yeast cells growing in microplates are likely to be oxygen starved (even if being shaken) and near transition states in metabolism, so it makes sense that changes in mitochondrial "states" are being observed. A discussion of how environmental conditions could influence heterogeneity would help to put this into context.

We have mentioned this now in the Discussion section (second paragraph).

[Editors' note: further revisions were requested prior to acceptance, as described below.]

The reviewers remained supportive of your study and appreciated most of the revisions. However, the evidences supporting the claim that cells switch from a low-growth to a fast-growing mode remained inconclusive. Since this is not the essence of your results, we are happy to accept the paper. however, we ask that you will back off this claim and move Figure 5 to the supplementary material.

The reviewers also felt that the discussion of previous findings, and their relation to your work, will improve the manuscript. please try to do that as well. see below specific comments.

We thank the editor and the reviewers again for their constructive suggestions.

We have moved Figure 5 to the supplementary section (Figure 5—figure supplement 6, 7) and have minimized the discussion on switching in the manuscript.

We have also added a part discussing the relevance of our work in the context of already published work in the Discussion.

Reviewer #1:

In this revised manuscript, Dhar et al. present new analyses aiming to address concerns about their inferences of switching from a slow-growing, respiration-deficient cell state to a fast-growing, respiration-capable cell state. I remain enthusiastic about the topic of the manuscript and I remain impressed by the scale of the mutant analysis, but I also remain concerned that the authors' analyses of switching rely on indirect arguments that are not entirely convincing. Specifically:

1) They argue (subsection “Switching of cells from high to low membrane potential”) that TMRE-high to TMRE-low switching must take place to explain the observation that 1% of cells are TMRE-low immediately after sorting into the TMRE-high bin (a low rate of contamination), whereas 34.7% of cells are TMRE-low 24 hours later. Their argument relies on a simple model in which TMRE-high cells consist of two subpopulations (80% slow-growing and 20% fast-growing), TMRE-low cells are fast-growing, and the average slow and fast growth rates represent the corresponding subpopulations. This model has some problems. First, the arithmetic mean growth rate might not be appropriate to model changes in subpopulation sizes over time. If each cell chooses a growth rate at random from its distribution then what matters is the geometric mean, which is heavily influenced by low values. Alternatively, if growth rate is relatively stably inherited through time then what matters is the fastest members of a subpopulation. Either way, arithmetic mean could be a very poor indicator of subpopulation growth. Second, even if we ignore the problem with the arithmetic mean and take the authors' argument and apply it to the transition between the 24-hr time point and the 48-hr time point, then their calculation actually predicts quite well the 48-hr subpopulation sizes from the 24-hr starting point, so TMRE-high to TMRE-low switching is not required to explain this transition. Third, the authors' simple model ignores lag, which would likely further under-represent the slow-growing (high-TMRE) fraction, and might actually be why the transition from 0 to 24 hours requires additional explanation whereas the transition from 24 to 48 hours does not.

We agree with the reviewer that our model is quite simple and we could make it more complex by considering many more sub-populations with different growth rates. However, this model gives us a conservative estimate for switching of cells from high to low TMRE state.

The geometric mean growth rates of slow and fast sub-populations are 0.20 and 0.37 h-1 (as compared to arithmetic mean growth rate of 0.21 and 0.38 h-1 respectively). This does not substantially change our calculations of expected% of low TMRE cells without any switching from high to low TMRE state (subsection “Switching of cells from high to low membrane potential”). We did not consider any lag phase as it is difficult to determine the duration of the lag phase from our results and it seems that the lag phase is more pronounced in fast growing cells (Author response image 1). Nevertheless, including lag phase in our calculation would reduce the time for exponential growth for slow and fast-growing cells, which means there would be less time for fast growing low TMRE cells to take over the population. For example, in our calculation in the aforementioned subsection, if we consider 2 hours of lag phase for both slow and fast-growing cells, the expected% of low TMRE cells comes out to be 4.11% which is lower than what we expect without considering lag phase.

We also show calculations for expected% of low TMRE cells after 48 hours and again, the expected% of low TMRE cells is substantially lower than the observed% of low TMRE cells.

2) The authors have verified that their analyses of growth rate switching in micro-colonies are robust to changes in cutoffs between "fast" and "slow" growth, which is good, but what these switches mean is insufficiently discussed. In particular, fast-to-slow switches such as those shown in Figure 5D necessarily involve concerted changes in most cells in a micro-colony because a single cell changing from fast to slow in the middle of micro-colony growth would not appreciably change the growth rate of the micro-colony. The authors should discuss (and perhaps experimentally investigate) what such switches at the level of the micro-colony mean.

Indeed, we do not have the mechanistic understanding of the switching at the microcolony level and how many cells in a microcolony are switching their growth rates is not clear from our data. Hence, we have moved these results to the supplement (Figure 5—figure supplement 7, 8, 9) and have minimized discussion of switching in the manuscript.

3) The conclusion that some cells regain ability to respire is an interesting piece of this work but it is not clear exactly what the authors mean by it. In Figure 5E they present the "percentage of colonies that regained capability to respire" but the calculation involves two different plates that were not replicas, so the ability to "regain" respiration was a more indirect inference than the text suggests.

Because of the above concerns, my view remains that more direct evidence of switching at the level of single cells (from high-TMRE to low-TMRE, from slow growth to fast growth, and from respiration-deficient to respiration capable) would substantially strengthen this paper.

The replica plating experiment for assessing the ability to regain respiration required us to know the time of switching which we were not able to determine through this experiment. Thus, we have removed this claim from the results part and have moved Figure 5 to the supplement (Figure 5—figure supplement 6, 7).

Reviewer #3:

I completely agree that the resubmission only indirectly addresses our concerns about slow to fast growth. Even with the new statistical analysis, some quick growing contaminants could explain what they conclude are transitions. A simple replica plating experiment would have shown that their non-respiring cells can convert back to respiring cells. Given that they were previously asked for alternate experimentation to demonstrate the slow to fast transition, they need to back off this claim. (and minimize and move Figure 5 to supplemental).

We have moved Figure 5 to the supplementary material (Figure 5—figure supplement 6, 7) and also have minimized discussion on this topic in the manuscript.

In general, the authors could do better at discussing their findings in the context of previously published work.

We have now more extensively discussed our findings in the context of earlier work in the Discussion.

Their proposal (in the Discussion) that petites are part of a mitochondrial heterogeneity continuum is not substantiated if they are unable to convincingly show reversibility.

We have removed this proposal from the main text.

In the Abstract, they should change the word "developed" a high-throughput method to "used".

Done.

The authors now refer to mitochondrial membrane potential in places though there are still several places where mitochondrial "state" is used and could/should be corrected (e.g. Subsection “Mitochondrial membrane potential but not amount predicts slow growth”, fourth and fifth paragraphs, subsection “Variation in mitochondrial state predicts additional phenotypic heterogeneity including drug resistance”, first paragraph and especially in the first paragraph of the Discussion).

Changed now.

https://doi.org/10.7554/eLife.38904.043

Article and author information

Author details

  1. Riddhiman Dhar

    1. Systems Biology Program, Centre for Genomic Regulation (CRG), The Barcelona Institute of Science and Technology, Barcelona, Spain
    2. Universitat Pompeu Fabra, Barcelona, Spain
    3. Department of Biotechnology, Indian Institute of Technology, Kharagpur, India
    Contribution
    Conceptualization, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing—original draft, Writing—review and editing, set up high-throughput microcopy assays, developed pipelines for image processing, performed high-throughput microscopic proliferation measurements of the deletion mutants and analyzed the data, performed the remaining experiments, analyzed the data, and wrote the manuscript
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-4642-0492
  2. Alsu M Missarova

    Department of Experimental and Health Sciences, Universitat Pompeu Fabra, Barcelona, Spain
    Contribution
    Conceptualization, Formal analysis, Investigation, Writing—review and editing, performed high-throughput microscopic proliferation measurements of the deletion mutants, analyzed the data and contributed to writing the manuscript
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-9472-2095
  3. Ben Lehner

    1. Systems Biology Program, Centre for Genomic Regulation (CRG), The Barcelona Institute of Science and Technology, Barcelona, Spain
    2. Universitat Pompeu Fabra, Barcelona, Spain
    3. Institució Catalana de Recerca i Estudis Avançats (ICREA), Barcelona, Spain
    Contribution
    Conceptualization, Supervision, Funding acquisition, Writing—original draft, Writing—review and editing
    Contributed equally with
    Lucas B Carey
    For correspondence
    lehner.ben@gmail.com
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8817-1124
  4. Lucas B Carey

    Department of Experimental and Health Sciences, Universitat Pompeu Fabra, Barcelona, Spain
    Present address
    Center for Quantitative Biology, Peking-Tsinghua Center for Life Sciences, and the Academy for Advanced Interdisciplinary Studies, Peking University, Beijing, China
    Contribution
    Conceptualization, Supervision, Funding acquisition, Writing—original draft, Writing—review and editing
    Contributed equally with
    Ben Lehner
    For correspondence
    lucas.carey@pku.edu.cn
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-7245-6379

Funding

H2020 European Research Council (616434)

  • Ben Lehner

AXA Research Fund

  • Ben Lehner

Ministerio de Economía y Competitividad (BFU2011-26206)

  • Ben Lehner

Fondation Bettencourt Schueller

  • Ben Lehner

Ministerio de Economía y Competitividad (BFU2015-68351-P))

  • Lucas B Carey

Agència de Gestió d’Ajuts Universitaris i de Recerca

  • Ben Lehner
  • Lucas B Carey

Schweizerischer Nationalfonds zur Förderung der Wissenschaftlichen Forschung

  • Riddhiman Dhar

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

Work in the lab of BL was supported by a European Research Council Consolidator grant (616434), the Spanish Ministry of Economy and Competitiveness (BFU2017-89488-P and SEV-2012–0208), the AXA Research Fund, the Bettencourt Schueller Foundation, Agència de Gestió d’Ajuts Universitaris i de Recerca (AGAUR SGR 1322), the EMBL Partnership, and the Generalitat/CERCA program. Work in the lab of LBC was supported by MINECO (BFU2015-68351-P) and AGAUR (2014 SGR 0974) and the Unidad de Excelencia Maria de Maeztu (MDM-2014–0370). The authors thank Dr. Raul Gomez Riera and the CRG Advanced Light Microscopy Unit for assistance with high-throughput microscopy, the CRG/UPF FACS facility for help with flow cytometry, the CRG Genomics Unit for RNA-seq and DNA sequencing experiments, and the UPF Genomics facility for help with automated liquid handling. RD was partially supported by the Swiss National Science Foundation Early Postdoc Mobility fellowship.

Senior and Reviewing Editor

  1. Naama Barkai, Weizmann Institute of Science, Israel

Publication history

  1. Received: June 4, 2018
  2. Accepted: January 13, 2019
  3. Accepted Manuscript published: January 14, 2019 (version 1)
  4. Version of Record published: February 7, 2019 (version 2)

Copyright

© 2019, Dhar et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 1,893
    Page views
  • 314
    Downloads
  • 1
    Citations

Article citation count generated by polling the highest count across the following sources: PubMed Central, Crossref, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Genetics and Genomics
    2. Neuroscience
    Raquel Barajas-Azpeleta, Kausik Si
    Insight
    1. Developmental Biology
    2. Genetics and Genomics
    Samuel Andrew Malone et al.
    Research Article