Single-cell RNA sequencing reveals intrinsic and extrinsic regulatory heterogeneity in yeast responding to stress.

Strains and growth conditions

All experiments were done in the BY4741 background. Unless noted, cells were grown in rich YPD medium in batch culture at 30 °C for at least seven generations to mid-log phase, at which point an aliquot was removed to serve as the unstressed sample. NaCl was added to a final concentration of 0.7 M in the remaining culture and cells were grown for 30 min. Unless otherwise noted for specific applications, cells were collected by brief centrifugation, decanted, and flash frozen in liquid nitrogen. Strains expressing tagged proteins were generated by integrating an mCherry- HIS3 cassette downstream of MSN2 in BY4741 strains from the GFP-tagged collection [ 108 ], which were verified to harbor the GFP– HIS3 cassette downstream of DOT6 , SFP1 , or HOG1 (generating strains AGY1328, AGY1329, and AGY1331, respectively). Single-cell sorting, library preparation, and sequencing Fluidigm’s C1 microfluidic platform was adapted to perform cDNA synthesis from single yeast cells. Flash frozen cells were resuspended in 1 mL of 1 M Sorbitol on ice, counted on a hemocytometer, and then diluted to approximately 4 × 10 5 cells per mL in a final volume of 200 μL. To generate partial spheroplasts that could easily lyse on the Fluidigm C1 microfluidic device, we titrated each sample with different amounts of zymolyase (0.025 U, 0.0125 U, 0.00625 U, and 0.003125 U) and incubated cells for 30 min at 37 °C. This was done because unstressed and stressed cells displayed different sensitivities to zymolyase digest. After incubation, cells were spun at 250 g for four min and resuspended in Sorbitol Wash buffer (0.455x C1 Cell Wash Buffer, 1 M sorbitol, 0.2 ug/μl BSA, 0.08(8) U/μL SUPERase RNAse Inhibitor). Samples with the maximal number of intact spheroplasts (compared on a Leica DMI 6000 inverted microscope) were diluted to a final concentration of 600 cells/μL, and 9 μL of these cells were mixed with Fluidigm Suspension reagent at final loading concentration of 275 cells/μL and loaded onto the primed C1 Chip designed to capture 5–10-μm cells, according to manufacturer instructions. The cell concentration in the loading mixture was crucial to maximize the number of wells capturing single cells inside the microfluidic device. Another modification was that 1 M Sorbitol was added to all wash buffers to prevent premature lysis. After cell loading, each chip was visually inspected and imaged to tabulate single-cell capture rates. Roughly 50% of wells contained a single cell, verified by imaging and manual inspection. This rate is lower than the normal capture rate because yeast cells are smaller and deform less. Spheroplasts were lysed in the Fluidigm instrument and cDNA was generated using Clontech reagents for Fluidigm C1 based on the single-cell RNA-seq protocol (cat # 635025). Finally, cDNA was harvested from each Fluidigm C1 chip and into a 96-well plate for storage at −20 °C. ERCC spike-in sequences (mix A) were added at 1:4 × 10 5 /μL of the concentration provided in the original product (Ambion catalog number 4456740). Before library preparation, cDNA from each cell was quantified on an AATI Fragment Analyzer. Using concentrations calculated from a smear analysis between 450 bp and 4,500 bp, cDNA from each cell was diluted with TE to approximately 0.2 ng/μL using the Mosquito X1 pipetting robot (TTP Labtech). Diluted cDNA served as the template for Nextera XT library generation following manufacturer protocol (Illumina catalog number FC-131-1096) with some modifications. Because we used TTP’s Mosquito HTS 16 channel pipetting robot (capable of accurately aliquoting volumes down to 50 nL), we were able to scale down the total volume of each Nextera XT library to four μL. More specifically, for each 400 nL of input template DNA, we added 400 nL Tagmentation mix and 800 nL Tagmentation buffer for a final volume of 1.6 nL. The Tagmentation reaction was incubated at 55 °C for 10 min. Neutralization was done by adding 400 nL Neutralization buffer to the above reaction and incubating 10 min at room temperature followed by the addition of primers at 400 nL each and NPM PCR master mix at 1,200 nL. The PCR step was run for 10 cycles. One μL of each library was combined to form two separate pools, one for unstressed cells and one for stressed cells. Two rounds of size selection were performed using Agencourt AMPure beads (Beckman Coulter catalog number A63882). One hundred ng of each pool was combined and sequenced in one run on three lanes of an Illumina HiSeq-2500 1T v4 sequencer for 150-bp paired-end sequencing. Reads generated across the three lanes were merged and demultiplexed using Illumina software bcl2fastq v1.8.4, allowing no mismatches and excluding the last position (eighth index base). Paired-end reads were mapped to the S288c S . cerevisiae genome R64-2-1 [ 42 , 109 ] with ERCC spike-in sequences added using BWA mem Version: 0.7.12-r1039 and default parameters [ 109 ]. Reads were processed with Picard tools Version: 1.98(1547) cleansam and AddOrReplaceReadGroups as required by downstream applications. Resulting bam files were sorted and indexed using Samtools Version 1.2. Paired-end fragments were deduplicated using the RemoveDuplicate function in Picardtools, and read counts mapped to genes were extracted using FeatureCounts Version 1.5.0. Sequenced wells were removed from the analysis if they had 0.2 [ 110 ]. Data were normalized by SCNorm [ 111 ] in R version 3.3.1; ERCC spike-in samples were not used in the normalization. Normalized read counts for each gene were logged and then centered by subtracting the mean log 2 (read counts) for that gene across all cells in the analysis. All data are available in the NIH GEO database under access number GSE102475 .

Statistical analysis of differential expression

Genes defined in iESR, RP, and RiBi clusters [ 24 ] are annotated in S8 Table . Individual cells with altered expression of defined gene groups (e.g., cells with low RP expression as in Fig 1C or high Rpn4/Hsf1 targets as in Fig 6 ) were determined using a two-tailed Welch t test, comparing the set of mean-centered log 2 (read count) values in that cell to the combined set of mean-centered log 2 (read count) values for all other unstressed (or stressed) cells, taking FDR 1e−5). All significant tests shown in Fig 3 remained significant ( p < 0.05) except for RP transcripts after stress ( p = 0.36).

Cell classifications and gene clustering

Data were clustered with Pagoda [ 54 ] using default parameters, and clusters enriched for known cell-cycle regulators [ 65 ] or glycolysis transcripts were manually identified. We were unable to identify known cell-cycle markers that peaked in M phase, either in the Pagoda-clustered data or using known M-phase markers [ 65 ]. To classify cells according to cell-cycle phase, cells were organized by hierarchical clustering [ 113 ] based on the centroid (median) of mean-centered log 2 (transcript abundance) for transcripts in each cell-cycle phase ( S3 Table ) and manually sorted and classified based on overlapping peaks of each vector ( Fig 5A ). Cells that showed no expression peak in any of the vectors were scored as “unclassified”.

Identification of TF target expression differences

Compiled lists of TF targets were taken from [ 93 ]. We added to this an additional list of Msn2 targets, taken as genes with stress-dependent Msn2 binding within the 800 bp upstream region and whose normal induction during peroxide stress required MSN2/MSN4 [ 114 ] ( S7 Table ) and genes whose repression requires the DOT6/TOD6 repressors (defined here as genes with a 1.5X repression defect in two biological replicates of wild-type BY4741 and dot6Δtod6Δ cells responding to 0.7M NaCl for 30 min [ 53 ] ( S7 Table ). In total, we scored 623 overlapping sets of TF targets defined in various datasets [ 114 – 117 ] or summarized from published ChIP studies [ 118 ]. We identified TF targets with concerted expression changes in two ways. First, we identified the distribution tails in each cell, identifying all mRNAs in that cell whose mean-centered log 2 (read count) was ≥1.0 or ≥2.0 (i.e., two times or four times higher than the population mean). We then used the hypergeometric test to identify sets of TF targets enriched on either list, taking the lower of the two p -values for that TF–cell comparison. Comparable tests were done for the sets of genes whose relative abundance was ≤−1.0 or ≤−2.0 in each cell. TFs with −log 10 ( p -value) > 4 were taken as significant (equivalent to a cell-based Bonferroni correction, p = 0.05 / 623 tests = approximately 1e-4). We focused on TFs whose targets were enriched at the distribution tails in at least four cells, which is unlikely to occur by chance. Two sets of TF targets were removed because their targets heavily overlapped with ESR targets and their enrichment was not significant when those overlapping mRNAs were removed ( S6 Table ). In a second approach, we used Welch t tests to compare the mean-centered log 2 (read count) values of each set of TF targets compared to all other measured mRNAs in that cell, again taking −log 10 ( p -value) > 4 as significant and focusing on TFs identified in at least three cells. Finally, for Fig 6E and 6F , colored boxplots indicate TF targets whose relative abundances in the denoted cell were significantly different from all other stressed cells (FDR < 0.053). sm-FISH BY4741 was grown as described above and collected for fixation and processing as previously described [ 119 ]. FISH probes were designed using the Biosearch Technologies Stellaris Designer with either Quasar 670, CAL Fluor Red 610, or Quasar 570 dye, using 33 probes for PPT1 (Quasar 670) and RLP7 (Quasar 570) and 48 probes for SES1 (Quasar 570). PPT1 and RLP7 were measured in the same cells, and SES1 was measured in a separate experiment as a control. Images were acquired as z-stacks every 0.2 mm with an epifluorescent Nikon Eclipse-TI inverted microscope using a 100x Nikon Plan Apo oil immersion objective and Clara CCD camera (Andor DR328G, South Windsor, Connecticut, United States of America). Quasar 670 emission was visualized at 700 nm upon excitation at 620 nm (Chroma 49006_Nikon ET-Cy5 filter cube, Chroma Technologies, Bellows Falls, Vermont, USA). Quasar 570 emission was visualized at 605 nm upon excitation at 545 nm (Chroma 49004_Nikon ET-Cy3 filter cube). PPT1 and RLP7 transcripts were counted manually, while SES1 mRNA was counted by semiautomated transcript detection and counting in MATLAB using scripts adapted from [ 120 ].

Fixed-cell microscopy

Cells were grown as described above and fixed with 3.7% formaldehyde for 15 min either before or at indicated times after NaCl addition. Cells were washed two times with 0.1 M potassium phosphate buffer pH 7.5, stained 5 min with 1 μg/mL DAPI (Thermo Scientific Pierce, 62247), and additionally washed two times with 0.1 M potassium phosphate buffer pH 7.5. Cells were imaged on an epifluorescent Nikon Eclipse-TI inverted microscope using a 100x Nikon Plan Apo oil immersion objective. GFP emission was visualized at 535 nm upon excitation at 470 nm (Chroma 49002_Nikon ETGFP filter cube, Chroma Technologies, Bellows Falls, VT, USA). mCherry emission was visualized at 620 nm upon excitation at 545 nm (Chroma 96364_Nikon Et-DSRed filter cube). DAPI emission was visualized at 460 nm upon excitation at 350 nm (Chroma 49000_Nikon ETDAPI filter cube). Nuclear to cytoplasmic intensity values were calculated with customized CellProfiler scripts [ 121 ]. The fraction of cells with nuclear factor before stress was calculated by identifying nuclear masks in Cell Profiler (i.e., identifiable nuclear objects) in that channel that overlapped with DAPI masks and manually correcting miscalls.

Live-cell microscopy

Yeast strains AGY1328 (Dot6-GFP, Msn2-mCherry) and AGY1331 (Sfp1-GFP, Msn2-mCherry) were grown at 30 °C with shaking to OD 600 0.4–0.5 in low-fluorescence yeast medium (LFM) [ 119 ]; the fraction of cells with nuclear Msn2/Dot6/Sfp1 was very similar in fixed, unstressed cells growing in LFM versus YPD, not shown. Cells were imaged in a Focht Chamber System 2 (FCS2) (Bioptechs, Inc; Butler, Pennsylvania, USA) with temperature maintained at 30 °C. Cells were loaded into the chamber and adhered to a 40-mm round glass coverslip. Briefly, the coverslip was prepared by incubating with concanavalin A solution (two mg/ml in water) for two min at room temperature. The concanavalin A was aspirated, and 350 μl of cell culture added. Cells were allowed two min at room temperature to adhere to the concanavalin A before the gasket was washed once with 350 μl of fresh LFM media. A round 0.2-μM gasket contained the concanavalin A plus cell solution. The entire assembly, including coverslip, cells, and gasket, was assembled with the microaqueduct slide, upper gasket, and locking base—as per manufacturer instructions—to create an incubated perfusion chamber. This assembly was transferred to an epifluorescent Nikon Eclipse-TI inverted microscope for time-lapse imaging. Temperature was maintained at 30 °C using the Bioptechs temperature controller. Steady perfusion with LFM or LFM + 0.7 M NaCl (initiated at 81 min) was maintained utilizing microperfusion pumps feeding media through the FCS2 perfusion ports. Once loaded onto the microscope, cells and temperature were monitored for at least one hr to ensure stable temperature and robust growth. During each timecourse experiment, a TI-S-ER motorized stage with encoders (Nikon MEC56100) and PerfectFocus system (Nikon Instruments, Melville, New York, USA) were used to monitor specific stage positions and maintain focus throughout the timecourse. A Clara CCD Camera (Andor DR328G, South Windsor, Connecticut, USA) was used for imaging. GFP emission was visualized at 535 nm (50-nm bandwidth) upon excitation at 470 nm (40-nm bandwidth; Chroma 49002_Nikon ETGFP filter cube, Chroma Technologies, Bellows Falls, Vermont, USA). mCherry emission was visualized at 620 nm (60-nm bandwidth) upon excitation at 545 nm (30-nm bandwidth; Chroma 6364_Nikon Et-DSRed filter cube). Single-cell traces of nuclear localization were extracted from fluorescence images using custom Fiji [ 122 ] and Matlab scripts. Fiji was used to threshold and identify individual cells. Individual cell traces were constructed using a modified Matlab particle tracking algorithm [ 123 ] and then manually validated and corrected. Localization of individual transcription factors was quantified using a previously published localization score based on the difference between the mean intensity of the top 5% of pixels in the cell and the mean intensity of the other 95% of pixels in the cell [ 124 ]. Peaks of nuclear TF localization were identified using findpeaks2 (Matlab File Exchange).

Supporting information S1 Table

Statistics on sequencing features. (XLSX) Click here for additional data file.

S2 Table

Transcripts outside RP-splines and iESR-fit splines. iESR, induced-Environmental Stress Response; RP, ribosomal protein. (XLSX) Click here for additional data file.

S3 Table

Cell-cycle expression vectors and gene lists used for classification. (XLSX) Click here for additional data file.

S4 Table

Oscope-identified clusters of cycling genes. (XLSX) Click here for additional data file.

S5 Table

Wavecrest-identified transcripts associated with RP mRNAs. RP, ribosomal protein. (XLSX) Click here for additional data file. S6 Table Specific TF-target gene sets identified in Fig 6 analysis. TF, transcription factor. (XLSX) Click here for additional data file.

S7 Table

Lists of all TF targets considered for Fig 5 analysis. TF, transcription factor. (XLSX) Click here for additional data file.

S8 Table

Lists of genes included in iESR, RiBi, and RP clusters used in the analysis. iESR, induced-Environmental Stress Response; RiBi, ribosome biogenesis; RP, ribosomal protein. (XLSX) Click here for additional data file.

S9 Table

Summary of all cell classifications. (XLSX) Click here for additional data file.

S10 Table

FISH counts plotted in Fig 4B . Counts from unstressed and stressed cells are recorded for different transcripts on different tabs of the file. PPT1 and RLP7 were measured in the same cells; SES1 was measured alone in a different set of cells. FISH, fluorescence in situ hybridization. (XLSX) Click here for additional data file.

S11 Table

Nuclear versus cytosolic TF localization plotted in Fig 7A . Nuclear/cytosolic ratios from Cell Profiler are shown for each of two factors as measured in a single strain. Each strain data is shown on a separate tab. TF, transcription factor. (XLSX) Click here for additional data file.

S12 Table

Ratios from Fig 7C , as described in the text. (XLSX) Click here for additional data file. S1 Fig Mean iESR and RiBi transcript levels in cells with low RP mRNA. For each cell plotted, the median of the mean-centered log 2 (read count) values for the group of RP transcripts was plotted against the median of the mean-centered log 2 (read count) values for (A) iESR transcripts and (B) RiBi transcripts in that cell. Red points represent the cells in which iESR transcripts were concertedly high ( Fig 1C , asterisks). The RiBi-transcript detection rate (percentage of transcripts measured) for each cell is represented above the corresponding points in B, color-coded according to the key. R 2 is shown for all points and excluding red points in which iESR transcripts were relatively high. iESR, induced-Environmental Stress Response; RiBi, ribosome biogenesis; RP, ribosomal protein. (TIF) Click here for additional data file. S2 Fig iESR, RP, and RiBi abundance by cell-cycle phase. Boxplots (without whiskers) of (A) RP, (B) iESR, (C) RiBi mRNAs for cells classified in different cell-cycle phases. Each boxplot represents the distribution of relative mRNA abundances, as defined in Methods, within each cell. Significance was assessed by Welch t test on the pooled RP, iESR, or RiBi genes from cells within a given phase compared to the pooled set of those transcripts from all other cells; unstressed and stressed cells were analyzed separately. Dark blue boxplots in (A) represent significant groups, with FDR listed below. RP expression was slightly, but highly statistically significantly, higher among cells in G1 phase, particularly for subsets of cells. A subset of unstressed cells in S-phase showed particularly low mean-centered log 2 abundance of RP transcripts. In stressed cells, those in S-phase had particularly tight distribution of RP abundance, showing in effect weaker RP repression during stress than cells in other phases. There were no significant differences across cell-cycle phases for expression of iESR or RiBi transcripts. iESR, induced-Environmental Stress Response; RiBi, ribosome biogenesis; RP, ribosomal protein. (TIF) Click here for additional data file. S3 Fig No clear evidence for the YMC. (A) Mean-centered log 2 (read count) values of several groups of transcripts, in unstressed cells with low RPs, other unstressed cells, and stressed cells, as described in Fig 1C . Cells are ordered as in Fig 1 , except that unstressed, low-RP cells are displayed first. Gene groups include: 149 mRNAs annotated as part of the ESR RP cluster, 28 genes identified by Pagoda clustering that are heavily enriched for glycolysis transcripts, 23 and 40 YMC-related mRNAs identified by Tu et al . enriched for mitochondrial or peroxisomal (POX) functions; also shown are 82 transcripts encoding mitochondrial ribosomal proteins. (B) Boxplots (without whiskers) of the distribution of mean-centered log 2 (read counts) for the group of Pagoda-identified glycolysis transcripts, in individual cells ordered as in Fig 1C . Cells with statistically different expression of glycolysis mRNAs (FDR < 0.05) are highlighted in dark blue (unstressed and stressed cells were analyzed separately). Cells with lower RP expression from Fig 1C are indicated with arrows and asterisks, as described in Fig 1C . Cells with low-RP expression had no statistically significant difference in glycolysis-mRNA expression. ESR, Environmental Stress Response; POX, peroxisomal; RP, ribosomal protein; YMC, yeast metabolic cycle. (TIF) Click here for additional data file. S4 Fig Correlation between scRNA-seq and bulk RNA-seq for specific gene sets. Correlation between bulk RNA-seq log 2 (fold-change) after 0.7 M NaCl and average of single cell measurements. Bulk data represent the average of three biological replicates. R 2 is shown for each plot. Correlations were substantially higher when transcripts with no reads in a given cell were scored as zero reads (rather than treated as missing values), which supports the notion that missing measurements are influenced biological variation in transcript presence per cell, rather than purely noise specific to scRNA-seq data. RNA-seq, RNA sequencing; scRNA-seq, single-cell RNA sequencing. (TIF) Click here for additional data file. S5 Fig Genes with more upstream Msn2 sites do not show significant differences in pre-stress activation. The number of Msn2 binding sites (CCCCT or CCCCC) within 800bp upstream of the gene starts was identified in Msn2 targets defined by Huebert et al (2012). The distribution of mean-centered log 2 (read counts) in unstressed cells for genes with at least five upstream elements (left) was not statistically significantly different from genes with one or two upstream elements (right) (Welch t test comparing the two groups in each cell). Cells are ordered as in Fig 1C . None of the cells showed concerted differences in Msn2 target expression before stress (FDR < 0.05, Welch t test, see Methods ). There was no relationship between the median relative log 2 read counts and cells with low RPs. (TIF) Click here for additional data file.

S1 Data

Representative timecourse of Msn2-mCherry and Sfp1-GFP before and after addition of 0.7 M NaCl. The pre-stress phase is indicated by a white scale bar, which turns red upon NaCl addition. We were unable to call troughs of nuclear Sfp1-GFP before stress. Apparent Sfp1-GFP nuclear depletion, without nuclear Msn2-mCherry during the unstressed phase, is seen in the cell at the far upper left corner; the nucleus is clearly in focus judged by the nuclear Msn2-mCherry concentration upon NaCl treatment. GFP, green fluorescent protein. (AVI) Click here for additional data file.

S2 Data

Representative timecourse of Msn2–mCherry and Dot6–GFP before and after addition of 0.7M NaCl. The pre-stress phase is indicated by a white scale bar, which turns red upon NaCl addition. GFP, green fluorescent protein (AVI) Click here for additional data file.