Single cell functional genomics reveals the importance of mitochondria in cell-to-cell variation in proliferation, drug resistance and mutation outcome

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 developed a high-throughput automated microscopy assay and used it to quantify the impact of deleting >1,500 genes in yeast. Mutations affecting mitochondria were particularly variable in their outcome. In both mutant and wild-type cells mitochondria state – but not content – 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.


Introduction 39 40
Isogenic populations often exhibit considerable phenotypic heterogeneity even in an identical 41 environment. One common phenotypic variation that has been observed in isogenic populations 42 of microbial and mammalian cells, including cancer cells is variation in proliferation rate [1][2][3][4][5][6][7][8][9]. 43 Phenotypic variations that are often coupled with variation in proliferation rate are the abilities of 44 an individual cell to survive stress and drug treatment [3,4]. In this regard, the existence of 45 variable expressivity is also common in human disease [28][29][30][31]. 53 54 Heterogeneity can arise due to stochastic fluctuations in biological processes taking place inside 55 cells. This can happen due to the small numbers of molecules involved in processes such as 56 transcription [32][33][34] or during stochastic partitioning of cellular components during cell division 57 [35,36]. Although cell-to-cell variation in the expression level of single genes has been 58 correlated with variation in proliferation rate and stress and drug resistance [3,4,15,[24][25][26]37,38], 59 the true underlying causes of such phenotypic heterogeneity are poorly understood. 60 61 To identify genes and cellular processes involved in the generation of phenotypic heterogeneity 62 we set up a high-throughput microscopy assay to quantify proliferation heterogeneity in a yeast 63 population. Using this assay, we quantify the impact of deletion of >1,500 genes on proliferation 64 heterogeneity. We present evidence that the variation in mitochondria state is an important 65 determinant of phenotypic heterogeneity in individual cells. We also show that mitochondria 66 state impacts gene expression and stress and drug resistance in individual cells. Taken together, 67 our work suggests an important role for an organelle in generating phenotypic heterogeneity 68 across individual cells in a homogenous environment. 69

Natural and lab yeast populations show proliferation heterogeneity 72 73
To investigate cell-to-cell variation in proliferation rates, we developed a high-throughput 74 automated time-lapse microscopy assay that measures the proliferation rates of thousands of 75 single-cells per plate as they grow into micro-colonies. The assay uses a microscope with laser-76 based autofocus for image acquisition and a liquid handling robot to minimize density-dependent 77 effects on proliferation. The data obtained are highly reproducible with mode proliferation rate of 78 a lab strain being 0.407±0.011 h -1 , (mean±sd) during >2 years of data collection (n=44 batches; 79 Fig 1A). individual variation in the outcome of mutations has been observed before in multicellular 97 organisms [23-25] but its relative occurrence has not been systematically quantified. 98

99
We therefore used the automated microscopy assay to quantify proliferation rate heterogeneity in 100 triplicate for 1,600 gene deletion mutants (Supplementary table 2, including 1,150 gene deletions  101 previously reported as affecting growth rates) [40,41]. We obtained reproducible data (where at 102 least two replicate measurements showed good agreement) for 1,520 deletions, with 1,112 of 103 these reducing the population proliferation rate in our experiment (Mann-Whitney U test, 104 FDR<0.1). 105 106 Deletion strains with similar population proliferation rates often showed strikingly different 107 degrees of intra-population heterogeneity ( Fig, 2A-C). At the single cell level, ~39% of all 108 mutants with a significant reduction in population proliferation rate (1112 mutants) showed 109 significantly higher variation in mutation outcome compared to the WT strain. Among these 110 mutants, ~13% had the same mode growth rate as the WT strain while showing higher variability. 111 However, almost all mutants (1111 or 1112) had a subset of cells proliferating at the same rate as 112 the bulk of the wild-type (WT) population (one sample Wilcoxon rank-sum test for overlap with 113 bulk WT distribution differing from zero, FDR<0.1; Supplementary fig. 1, Fig. 2D). Thus, a 114 highly variable outcome is actually the normal outcome for proliferation rate at the single cell 115 level when a non-essential gene is inactivated (Fig. 2D, Supplementary fig. 2A). To identify the determinants of this cell-to-cell variation in growth-rate and mutational impact 120 we classified each of the deletions by how it affected both the mode and distribution of cellular 121 proliferation rates (Supplementary table 2 Fig. 2A,B). Approximately 17% of the mutants 122 showed no change in either mode proliferation rate or percentage of slow sub-population (in 123 grey), whereas ~43% exhibited a change in mode proliferation rate but no change in slow 124 fraction (in light blue). Interestingly, 48 mutants reduced the slow fraction without any change in 125 mode proliferation rate (in red) and 97 mutants increased the slow fraction without altering the 126 mode proliferation rate (in blue). In addition, there were 78 mutants that reduced both the slow 127 fraction and the mode proliferation rate (in orange). Finally, 370 mutants reduced the mode 128 growth rate but increased the slow sub-population ( Fig. 2A). Across mutants, we observed a 129 strong inverse relationship between mean growth rate and noise (co-efficient of variation) (Fig.  130 2C), as has been observed for gene expression [42,43]. 131

132
To identify biological processes associated with changes in the slow growing sub-population, we 133 performed a GO functional enrichment analysis on genes in these categories (FDR<0.1). 134 Deletions causing the largest increase in the fraction of slow proliferating cells were highly 135 enriched for nuclear genes encoding mitochondrial proteins (Fig 2C,E). Among the mutants that 136 increased the slow fraction but also reduced mode growth rate (Fig. 2E, magenta), ~30% 137 localized to mitochondria (~1.2-fold enrichment), ~13% localized to the mitochondrial envelope 138 (>1.6-fold enrichment) and ~4.6% were involved in cellular respiration (~2-fold enrichment). In 139 particular, deletion of genes that localized to the mitochondrial envelope resulted in a large  fig. 2C), ruling out cell-to-cell variation in segregation of the organelle as 151 a driver of heterogeneity. However, signal from the mitochondria membrane potential dye 152 TMRE varied substantially across WT, mutants and natural strains (Fig 3A,B). This suggested 153 that the state of the mitochondriabut not their contentmight be driving proliferation 154

heterogeneity. 155 156
To determine if mitochondria state is correlated with variation in growth rate within a population 157 we sorted wild-type cells according to their TMRE staining and measured the fraction of slow 158 proliferating cells. The population with high TMRE was highly enriched for slowly proliferating 159 cells (Fig. 3C). This same fraction was also strongly enriched for respiration deficient cells (  In summary, we have shown here that mitochondria statebut not contentvaries substantially 246 and reversibly across individual yeast cells and that this is associated with cell-to-cell 247 heterogeneity in proliferation, mutation outcome, and stress and drug resistance. Laboratory 248 strains of yeast have long been known to generate respiratory deficient 'petite' colonies at quite 249 high frequency [44,67,68]. However, a slow growing sub-population of cells was observed in all 250 the laboratory, natural, and clinical strains that we tested (Fig. 1A,B). We propose that 'petite' 251 colonies are an irreversible extreme phenotype that is generated as part of the more general 252 variation in mitochondria state across single cells that we have identified here. if the left sub-population was bigger than the right sub-population, the percentage of slow 355 fraction was calculated as (% left sub-population-% right sub-population) and the percentage of 356 fast sub-population was set to zero. If the right sub-population was bigger than the left sub-357 population, the percentage of fast fraction was calculated as (% right sub-population-% left sub-358 population) and the percentage of slow sub-population was set to zero. performed a hypergeometric test as follows. Let us assume that in a group 'g' from screening 378 (for example, the group with increased slow fraction but no change in mode growth rate), out of 379 total N g genes, X g genes are associated with function f according to GOslim annotation. Let us 380 also assume that out of total N genes screened in our data, X genes belong to the functional class 381 f according to GOslim annotation. Thus, the probability that the group 'g' contains more number 382 of genes of functional class f than expected by chance alone is given by p=∑ To test for switching of microcolony growth from slow to fast state or from fast to slow state, 441 microcolony growth rates were calculated for each time point using linear regression with data 442 from 3 time-points and with R 2 >=0.9. Microcolonies for which growth rate changed from over 443 0.3 h -1 to less than 0.3 h -1 over time were counted as microcolonies switching from fast to slow 444 state. For fast to slow state transition, the transition must happen before last 4 time-points so as 445 to discount slowdown due to nutrient deprivation. Microcolonies for which growth rate changed 446 from less than 0.3 h -1 to more than 0.3 h -1 over time were counted as microcolonies switching 447 from slow to fast state. For this switching, the transition must happen after first 4 time-points so 448 as not to count lag to log phase switching.

RNA sequencing experiment and data analysis 472
Isolated total RNA (using MasterPure yeast RNA isolation kit (Epicentre)) was checked and 473 quantified using bioanalyzer. 200ng of total RNA for each sample was taken and was mixed with 474 4μl of 1:1000 dilution of ERCC spike-in mix1 (Thermo Fisher Scientific). Sequencing was done 475 in Illumina HiSeq with paired end 2×50bp reads. Quality of the sequenced reads was checked 476 using FastQC [82] and then the reads were mapped to reference yeast transcriptome (R64-1-1 477 reference cdna sequence from Ensembl [83]) using bowtie2 [84]. Mapping statistics was 478 calculated using a custom script where only read pairs mapping concordantly and uniquely to the 479 reference sequence were considered. The data were normalized using ERCC spike-in reads as 480 controls using RUVg method in R package RUVSeq [85]. Correlation between replicates were 481 checked through distance heatmap and PCA analysis (Supplementary fig. 12A,B), using R 482 package DESeq2 [86]. Differentially expressed genes were identified using package DESeq2. 483 Functional enrichment analysis on sets of differentially expressed genes was done using a 484 hypergeometric test as described above with multiple testing correction (FDR<0.1) (Benjamini-485 Hochberg method) with GOslim gene annotations. 486 487

Reconstruction of single mutants 488
Gene deletion mutants were remade in the WT strain using sequence specific homologous 489 recombination. First, the deletion cassette from the appropriate deletion strain from the collection 490 was amplified using primers such that the amplified region contained the deletion cassette with 491 KanMX marker and 50-300bp of overhang on either side of the cassette. Particular care was 492 taken to avoid neighbouring genes from being amplified. The PCR product was transformed into 493 competent yeast cells (prepared using lithium acetate and PLImade by mixing 1ml water, 1ml 494 1M lithium acetate and 8ml 50% PEG3350) and colonies were selected on G418 plates. Two 495 verified clones for each mutant were picked for experiments. Some of the mutants associated 496 with mitochondrial function were found to be compensated in the deletion collection 497

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

Data availability 528
RNA-sequencing data that support the findings of this study have been deposited in NCBI GEO 529 with the accession code GSE104343 and reviewer token "mtabqsqebncphaz".