ReporterSeq measures how endogenous genes affect the expression of an exogenous reporter gene. It accomplishes this through pairing a genetic perturbation with barcode levels from a pathway-specific readout (Figure 1). Though ReporterSeq can be implemented with any encodable genetic perturbation (e.g. RNAi, CRISPR knockouts), we used CRISPRi (Gilbert et al., 2013) as a perturbation to lower expression of each gene and paired it with a barcoded mRNA encoding GFP downstream of an Hsf1-responsive synthetic promoter built upon a ‘crippled’ CYC1 promoter sequence derived from a 225 nucleotide fragment of the CYC1 promoter (Brandman et al., 2012; Guarente and Mason, 1983). This synthetic promoter, previously used in a fluorescence-based HSR screen (Brandman et al., 2012; Guarente and Mason, 1983), is sensitive to Hsf1 activity as well as proteins that may regulate the CYC1 sequence fragment (note that no transcription factors other than Hsf1 are known to bind to the sequence). To maximize dCas9-Mxi1 (a transcription inhibitor) and sgRNA expression, we expressed dCas9-Mxi1 and the ReporterSeq cassette on two separate high copy plasmids that were coexpressed in each cell. We verified that the genetic knockdowns could modulate the heat shock reporter by measuring the effects of known HSR regulators (Brandman et al., 2012) on reporter output (Figure 1—figure supplement 1A). We therefore concluded that CRISPRi knockdowns were effective in our system and capable of altering HSR reporter activity.

Figure 1 with 1 supplement see all

Download asset Open asset

The first step of ReporterSeq is to synthesize a diverse library of constructs, each containing one sgRNA and one barcode driven by a promoter of interest. The library is then introduced into cells such that each cell contains one construct. Cells are grown in the desired conditions and then harvested at the desired timepoints. RNA and DNA are harvested from the cells and the counts of each barcode are tallied using deep sequencing. The RNA counts for a given barcode are proportional to the total transcriptional output of all cells containing that barcode, while the DNA barcode counts are proportional to cell numbers. Thus, the RNA/DNA ratio of a given barcode reflects the activity of the reporter under the influence of the specific CRISPR knockdown corresponding to that barcode. This ratio is used to identify the effect of a given knockdown on the HSR in a single condition. In situations where the same pool of yeast is divided and exposed to multiple conditions, RNA counts between these samples can be directly compared to compute gene–stressor interactions (discussed below).

Our library targeted each gene in the genome, including genes producing non-coding RNAs, with up to 12 sgRNAs and paired multiple random barcodes with each sgRNA. Multiple sgRNAs per gene were necessary due to variability in efficacy of each knockdown and to minimize the influence of off-target effects from any one guide RNA. Using multiple barcodes for each sgRNA minimized the effects of barcode bias on stability of the reporter mRNA. We cloned a library of approximately 1 million random barcode-sgRNA pairs and performed paired-end high-throughput sequencing on the library to identify these pairs. Paired-end sequencing revealed that most genes were targeted in the library by all 12 designed sgRNAs (Figure 1—figure supplement 1B). Additionally, there was a wide distribution of barcodes associated with each sgRNA, with the mean representation being approximately 14 barcodes per sgRNA (Figure 1—figure supplement 1C). Essential genes were only weakly depleted from the library, reflecting that our screen measures the effect of knockdown and that many sgRNAs do not lead to effective knockdown (Figure 1—figure supplement 1D). Note that after multiple generations of growth, we expect that only one barcoded construct will be expressed in each cell even in the rare event that two constructs initially enter the same cell. Thus, the plasmid library was poised to measure the effect of each genetic knockdown on the output of the HSR reporter.

We first used the ReporterSeq library to measure the effect of each genetic knockdown on the HSR reporter in untreated, log-phase yeast. The full-genome library was transformed into >100,000 cells, resulting in a subset of the plasmid library expression in cells. RNA and DNA barcodes were sequenced and matched to their corresponding sgRNAs, revealing that 91% of genes were represented by six or more sgRNAs where both RNA and DNA counts were above 10 (Figure 1—figure supplement 1E) and 74% of these sgRNAs were represented by two or more barcodes (Figure 1—figure supplement 1F). To test whether the RNA/DNA ratio of barcodes reflected activation of the HSR under a specific perturbation, we assessed the effect of two predicted regulators of our synthetic HSR reporter (Figure 2A). sgRNAs targeting SSA2, an Hsp70 chaperone that binds Hsf1 and inhibits the HSR (Krakowiak et al., 2018; Zheng et al., 2016), resulted in a higher than average RNA/DNA ratio as expected. Conversely, sgRNAs targeting CYC1, the gene from which our synthetic reporter was derived and shares significant homology with, lowered the RNA/DNA ratio compared to average. Thus, the DNA/RNA ratio for a single barcode reflects the HSR reporter’s activity under a specific genomic perturbation. Note that since we do not explicitly measure knockdown efficiency of the CRISPR constructs, our results cannot rule out the effect of any gene on the HSR.

Figure 2 with 3 supplements see all

Download asset Open asset

To facilitate comparisons between many data points, we developed a scoring method called the neighborhood normalized score (NNS) that combined the information from the multiple barcodes and corresponding sgRNAs which target a single gene into a single score. Briefly, each barcode was assigned an NNS, representing how extreme the RNA/DNA ratio is for the barcode, relative to its ‘neighborhood’ of barcodes with similar levels of expression. The magnitude of the score for a given barcode corresponds to the number of standard deviations from the mean, where both the standard deviation and mean are computed based on this neighborhood. The sign of the score indicates the direction of the effect on the HSR (positive NNS indicate that the knockdown increases the HSR, while negative NNS indicate the knockdown decreases the HSR). Most barcodes have NNS near 0, reflecting the fact that most perturbations do not significantly affect the HSR. The NNS for each gene was tallied as the average barcode NNS over all sgRNAs and barcodes corresponding to the gene in question. (See the Methods and methods section for details on how the NNS is calculated.) To give two examples, SSA2 knockdown, had the fourth highest NNS at 7.8, indicating that the knockdown of SSA2 likely increased reporter activity (Figure 2B), with the corresponding barcodes being 7.8 standard deviations above the mean, on average. CYC1, conversely, had the third lowest NNS of −5.5, meaning that it likely decreased the reporter activity. All gene NNS are provided in supplementary data file 1 and fold changes for top and bottom 20 NNS hits (e.g. SIS1 = ~100% increase, ERV1 = ~30% reduction) are provided in supplementary data file 2. Individual sgRNA NNS as well as barcode NNS for the top scoring sgRNA are shown for CYC1, SSA2, and, another top hit, SIS1 (Figure 2—figure supplement 1). Individual sgRNA counts and NNS for all genes are provided in supplementary data file 3 (basal NNS) and supplementary data file 4 (all stress conditions).

We quantified the effect of each genetic knockdown on the basal HSR reporter activity using our scoring metric and used a one-sample t-test to quantify the reproducibility of each NNS in the four replicates that we performed (Figure 2C). CRISPRi knockdown of chaperones or cochaperones SSA2, SIS1 (NNS = 16.0), FES1 (NNS = 4.1), SSE1 (NNS = 5.3), HGH1 (NNS = 4.0), SSB2 (NNS = 3.6), and STI1 (NNS = 3.1) strongly activated the HSR. Indeed, the Hsp40 SIS1 was the strongest basal suppressor of the HSR. Additional strong HSR suppressors include the negative regulators of the Ras/cAMP pathway GPB2 (NNS = 3.8) and IRA2 (NNS = 4.1), the mitochondrial lyso-phosphatidylcholine acyltransferase TAZ1 (NNS = 6.0), the nonsense-mediated decay regulators UPF3 (NNS = 3.3) and NMD2 (NNS = 3.0) and the small ribosomal protein ASC1 (NNS = 9.3). LST8 (NNS = 13.1) was also a strong hit in the screen, likely because it shares a promoter with SIS1 and thus silenced by the same sgRNAs (Pincus et al., 2018). Our implementation of ReporterSeq cannot distinguish between adjacent genes that may be targeted by the same sgRNAs, a general shortcoming of CRISPRi. In addition to CYC1, several mitochondrial genes (ERV1 NNS = −5.8, RSM27 NNS = −3.9, and TIM12 NNS = −3.1) and translation genes (RPL25 NNS = −5.6, RPL12A NNS = −3.3, and RPS18A NNS = −3.1) strongly lowered the HSR when impaired, consistent with previous observations that knockout of mitochondrial genes lowers the HSR (Brandman et al., 2012). To test whether top hits in the screen were dependent on the Hsf1 binding sites in our reporter rather than the crippled CYC1 sequence, we compared the effect of deleting three of the top hits (SSA2, HSC82, and ASC1) on the HSR reporter versus a crippled CYC1 reporter without the Hsf1-specific binding sequence (Figure 2—figure supplement 2). All deletions activated the HSR reporter without activating the crippled CYC1 reporter, which has very low basal activity, suggesting that most hits in the screen are specific to Hsf1.

In order to identify processes that regulate the HSR, we analyzed which gene ontology (GO) categories (Ashburner et al., 2000; The Gene Ontology Consortium, 2019) were enriched for genes causing high or low activation of the HSR when knocked down. We assigned each GO category a score based on the Student's t-test between the scores of all genes and those of the specific category (supplementary data file 5). As with the interaction score, the sign indicates the direction of the effect of the knockdown on the reporter activity, and the magnitude of this score indicates the confidence. Many of the GO category scores followed trends largely consistent with previous studies, with chaperone, protein folding, and protein quality control complexes such as the ribosome-associated quality control complex (RQC) amongst the most highly scoring GO categories (Figure 2D). The highest scoring GO category was negative regulation of Ras signal transduction, indicating that high Ras activity leads to high HSR activity. This is consistent with a previous study showing an inverse correlation between the HSR and the Ras/PKA-inhibited Msn2/4 pathway (Brandman et al., 2012) and studies concluding that the PKA activity induces the HSR (Murshid et al., 2010; Verma et al., 2020) The lowest scoring categories included mitochondria and translation, consistent with results of a previous HSR screen (Brandman et al., 2012). Knockdowns of genes involved in vacuolar acidification, a process important for many cellular activities such as protein sorting and pH maintenance (Yamashiro et al., 1990) as well as regulation of the HSR (Triandafillou et al., 2018) also caused a significant downregulation of the HSR.

To determine which genes affect the cell’s ability to activate the HSR in response to heat stress, we sought to quantify the interaction between the effects of heat stress and those of each genetic knockdown (a ‘stressor–gene’ interaction) on the output of the HSR reporter (Figure 3A), similar to how genetic interactions (‘gene–gene’ interactions) are conceptualized (Collins et al., 2007; Tong et al., 2001). Stressor–gene interactions are based on a comparison of the reporter output of the combined condition relative to that of the knockdown and stress independently. A stressor–gene interaction NNS of zero denotes that knockdown of the gene and the stressor itself have independent effects on the HSR. A positive NNS denotes that the combination of gene and stress resulted in a superadditive effect on the HSR (an ‘aggravating’ interaction), while a negative NNS denotes that the combination results in a subadditive effect on the HSR (a ‘blocking’ interaction). A positive or negative interaction NNS suggests that a stressor and gene are related in the context of the HSR. Traditional methods to compute interaction scores require three measurements: the phenotype of each individual perturbation and the phenotype of the double mutant. By assuming that (1) the majority of genetic knockdowns do not interact with the stressor and (2) short stressor treatments that impair cell growth (relative to the ~2.5 hr doubling time of unstressed yeast in the growth conditions used) do not significantly alter cell counts for a specific CRISPRi perturbation, we computed stressor–gene interactions by comparing RNA levels between stressed and unstressed conditions from the same yeast sample using the same scoring metric as with the untreated RNA/DNA ratios. This allowed us to compute the stressor–gene interactions without calculating the absolute induction of the gene (its RNA/DNA ratio), reducing noise arising from DNA normalization and increasing robustness by performing more comparisons than traditional methods for computing gene–gene interactions.

Figure 3

Download asset Open asset

We then examined which genes showed strong interactions with a heat stress of 20 min at 39°C. As an example of the interaction NNS in practice, we first examined SSA2 and CYC1. SSA2 had a blocking interaction with heat (NNS = −7.4) whereas CYC1 had an aggravating interaction with heat stress (NNS = 2.1; Figure 3B). This suggests (1) that heat stress and SIS1 depletion activate the HSR through common mechanisms and (2) that targeting both the CYC1 gene as well as the HSR reporter directly (because the reporter is built upon the CYC1 promoter sequence) do not reduce its inducibility by heat. Indeed, induction was higher than expected in knockdown of CYC1 (a positive interaction) and, generally, heat stress interaction NNS and the basal knockdown NNS were negatively correlated (r = −0.45), suggesting that knockdowns basaly increasing the HSR generally blocked HSR activation in heat and those decreasing basal HSR resulted in a bigger HSR increase during heat stress.

One knockdown target that blocked the response to heat less than would be predicted based on the general trend for strong HSR inducers was ASC1 (Figure 3C). This suggests that activation of the HSR by ASC1 may be through mechanisms independent of the other genes that strongly activated the HSR, such as SIS1 and SSA2. To test this, we measured the genetic interaction between combined genetic perturbations (deletions combined with CRISPRi) to determine if these genes belonged in a shared or independent pathway of activating the HSR (Figure 3D). As expected, SIS1 (encoding a cochaperone of Hsp70) and SSA2 (encoding an Hsp70) had a blocking interaction with each other, as the combination of the two perturbations was 65% of the expected GFP levels (quantified as log GFP levels). By contrast, ASC1 was independent of SSA2, with the combination of genetic perturbations having an almost completely additive effect on the reporter (GFP levels 94% of expected levels). This suggests that lowering Asc1 levels may activate the HSR through a separate mechanism than lowering Hsp70 levels.

The pooled format of ReporterSeq and very short sample collection time allowed us to measure the effect of each knockdown on the HSR in a time-resolved way with a precision that genome-wide approaches generally cannot achieve. We collected cells at 10, 20, 40, and 120 min after transition to 39°C in duplicate and isolated RNA to calculate stressor–gene interactions each time point (Figure 4A). CRISPR knockdown of two sample genes, SIS1 and SSA2, were not affected by heat shock duration, suggesting that knockdown efficiency does not change during the time course (Figure 4—figure supplement 1). Many genes had similar interaction NNS for the first and last timepoints. For example, knocking down SIS1 or TPS1 blocked the HSR throughout the time course (Figure 4B). Tps1 is required for synthesis of trehalose and has been demonstrated to be required for induction of the HSR under heat shock (Conlin and Nelson, 2007). Other interactions, however, greatly depended on the length of heat shock. Knockdown of the Ras gene, RAS2, caused an aggravation of the HSR at 10 and 20 min, had no effect at 40 min, and then subsequently blocked the HSR at 120 min. Knockdown of HSF1 had a strong blocking effect at 10 min that tapered to near zero by 40 min. Conversely, knockdown of HSP42, BTN2, or OPI10 had relatively little effect until 40 and 120 min, at which point the HSR became hyper-induced. These genes have been reported to be critical for survival in heat shocked cells, with Hsp42 and Btn2 sequestering toxic proteins (Haslbeck et al., 2004; Malinovska et al., 2012; Specht et al., 2011) and the human and Schizosaccharomyces pombe orthologs of Opi10 (‘Hikeshi’) facilitating the stress-induced import of Hsp70 into the nucleus (Kose et al., 2012; Song et al., 2015). The delay in the effects of these factors likely represents the time in which cellular damage accumulates under heat shock. Curiously, the nearly identical, highly expressed HSP70 paralogs SSA1 and SSA2 had strikingly different interaction profiles over time, with SSA2 having strong blocking interactions with all time points and SSA1 mostly not interacting with heat stress but having aggravating interactions at 40 and 120 min. These distinct profiles suggest differential regulation of the HSR by SSA1 and SSA2 over prolonged heat stress.

Figure 4 with 1 supplement see all

Download asset Open asset

To identify cellular processes that affected the early and late responses to heat shock differently, we analyzed which GO categories had differential effects over the time course (Figure 4C; Supplementary file 5). Like the aforementioned GO scores, these GO scores were calculated based on the Student’s t-test between the scores of all genes and those of the specific category. Knockdown of genes that are a part of the heat shock response, such as ‘mRNA export from nucleus in response to heat stress’, generally exacerbated the heat shock response at later time points, suggesting that if the response to heat is impaired, the stress response will be prolonged. Another mRNA regulator, the Lsm1-7-PatI complex, which degrades some mRNAs, also caused a stronger late-phase heat shock response. Additionally, genes involved in Ras signal transduction followed the same trend as RAS2 and were unusually high at early time points but low at later ones. Conversely, negative regulators of Ras showed the opposite trend. The lowest interaction scores found in later time points were proteasomal genes, while the highest scores were for ribosome biogenesis genes. ReporterSeq can thus identify cellular components that are critical for different phases of a stress response.

We next performed ReporterSeq on 11 additional stress conditions in duplicate, including glucose starvation and toxins (Supplementary file 6, Figure 5—figure supplement 1). These stressors affect cells through a variety of mechanisms and cause varying degrees of HSR reporter activation as measured by qPCR (Figure 5A). We applied these stressors for between 20 min and 3 hr. As expected, the pleiotropic proteotoxic stressors and amino acid analogs known to cause accumulation of misfolded protein (Jacobson et al., 2012; Trotter et al., 2002; Weids et al., 2016) robustly activated the HSR. The other stressors, including glucose starvation, and the ER stressors, DTT and tunicamycin, had milder effects on the HSR. We calculated interaction NNS for each gene in each stress condition (Figure 5B). Some of the stressors shared the same strong genetic interactions that we observed under heat shock conditions. For example, all stressors, with the exception of the proteasome inhibitor bortezomib, had a negative interaction NNS with the SIS1 and SSA2 knockdowns.

Figure 5 with 3 supplements see all

Download asset Open asset

Many genes had strong interactions in some stressors, but not in others. For example, knockdown of the Hsp90 co-chaperone STI1 (Chang et al., 1997; Nicolet and Craig, 1989) had aggravating interactions specifically with macbecin and canavanine, two stressors that have been demonstrated to reduce Hsp90 availability (Alford and Brandman, 2018), and unexpectedly, bortezomib (Figure 5C). Knockdown of GCN3, a gene involved in translation initiation regulation (Yang and Hinnebusch, 1996), had the strongest aggravating interactions in the screen, hyperinducing the response to hydrogen peroxide, and arsenite stress. Knockdown of SNF4, a Snf1 kinase subunit that plays a key role in glucose-repressed gene transcription (McCartney and Schmidt, 2001), had a strong aggravating interaction only with glucose starvation. Similar to what we observed in the heat shock time course, we saw distinct interaction profiles from the chaperone genes SSA1 and SSA2, with the SSA1 knockdown having generally weaker blocking interactions, as well as aggravating interactions with AZC and canavanine, two stressors that were strongly blocked by SSA2 knockdown. We verified the difference between SSA1 and SSA2 interactions with the HSR using gene knockouts and an integrated version of the HSR reporter used in a previous study (Alford and Brandman, 2018; Figure 5—figure supplement 2). Principal component analysis of chaperone gene NNS (chaperone list provided in Supplementary file 8) over all stressors revealed that the first two components scored strongest for SIS1/SSA2 and SSA1, respectively (Figure 5—figure supplement 3A). This suggests that the distinct interaction profiles of SIS1/SSA2 and SSA1 are characteristic of broad classes shared by other genes. The stressors that have the strongest effect on component one are the heat shock conditions and cycloheximide, indicative of the strong interactions these stressors had with SSA2/SIS1 (Figure 5—figure supplement 3A, B). Component 2, on the other hand, was driven most strongly with bortezomib, DTT, canavanine, macbecin, and tunicamycin. Exactly why chaperones like SSA1 had interactions with these stressors is unclear, but suggests a separate role in stress sensing and response from SSA2/SIS1.

To systematically evaluate which classes of genes are important for the HSR under different conditions, we analyzed which GO categories had strong interaction scores with each of the stress conditions (Figure 5D,E, GO categories provided in Supplementary file 7). We observed a wide diversity of stress-specific interactions. For example, knockdown of translation-related GO categories strongly aggravated the response to ethanol and glucose starvation, and to a lesser extent canavanine, AZC, and DTT. Translation genes had blocking interactions with macbecin and arsenite. Knockdown of mitochondrial genes strongly aggravated response to hydrogen peroxide and glucose starvation. Proteasome knockdown significantly blocked most responses, but had aggravating interactions with macbecin and bortezomib. ER-related genes had blocking interactions with DTT, tunicamycin, hydrogen peroxide, AZC, and canavanine. Enrichment scores for all GO categories are provided in Supplementary file 5. The diversity of genetic interactions observed between stressors and the HSR likely reflects effects of the knockdowns on the damage caused by the stressor as well as on stress-specific signaling to the HSR.

A prevalent model for HSR regulation is that the HSR senses changes in chaperone availability, becoming more active when chaperones become less available. In yeast, the Hsp70 and Hsp90 chaperones have both been proposed to have this role. Thus, lowering levels of either of these chaperones should cause constitutive HSR activation and a loss of sensitivity to stressors (i.e. a blocking interaction). Comparing HSR reporter mRNA induction in each stressor to interaction NNS for both SSA2 and HSC82 (one of two paralogs encoding Hsp90) revealed blocking interaction NNS with these chaperones that generally scaled with HSR induction strength, showing that both chaperones were important for responding to stressors (Figure 5—figure supplement 3C). Yet some interactions were specific to each chaperone. HSC82’s interaction with macbecin was stronger than with AZC, and SSA2’s interaction stronger with AZC, consistent with previous results (Alford and Brandman, 2018) and reflecting macbecin’s inhibition of Hsp90 and SSA2's role in sensing AZC stress.

Notably, cycloheximide had strong blocking interactions with chaperone genes, even though it does not substantially activate the HSR. This may occur through inhibition of translation, which lowers the load of misfolded proteins and requirements for chaperones in the cell and thus may block the HSR induction that occurs when chaperones are downregulated. Strikingly, plotting basal induction vs cycloheximide interactions revealed a negative correlation (correlation coefficient = −0.28), similar to the effect seen with heat shock (Figure 5—figure supplement 3D). Hierarchical clustering of stressors revealed that cycloheximide clustered with late heat shock (40 and 120 min; Figure 5F). Thus, many genetic perturbations that affect the HSR are dependent on translation, consistent with reports that the bulk of misfolded proteins in a cell are newly translated, nascent chains (Medicherla and Goldberg, 2008) and the dependence of the HSR on translation in heat stress (Feder et al., 2021). Translation-dependent HSR activation may explain blocking interactions observed in other stressors with low HSR induction, such as DTT, tunicamycin, and glucose starvation.

The knockdown with the strongest and most specific stress–gene interactions was the translation initiation factor, GCN3, displaying strong aggravating interactions with hydrogen peroxide and arsenite (Figure 5C). G is a component of eIF2B, the nucleotide exchange factor for eIF2. GCN3 has been shown to sense eIF2 phosphorylation and prevent nucleotide exchange of eIF2 thereby halting translation (Yang and Hinnebusch, 1996).

To further study the effects of GCN3 knockdown on the HSR, we used a single CRISPR sgRNA selected from the pool targeting GCN3. The HSR reporter encodes GFP, allowing us to measure its output by GFP fluorescence. We measured both GFP mRNA levels and fluorescence in a condition that was predicted to be aggravated by GCN3 knockdown, arsenite, as well as one that should be unaffected, AZC. GCN3 knockdown increased mRNA levels of the HSR reporter while reducing its protein levels under arsenite stress (Figure 6A), resulting in decreased translation efficiency (Figure 6B). As expected, no significant effect of GCN3 knockdown in both protein and mRNA levels is observed after treatment with AZC (Figure 6A, B). This suggests that GCN3 increases translation efficiency of the HSR reporter under arsenite stress.

Figure 6

Download asset Open asset

To determine whether Gcn3-dependent translation in arsenite was a general phenomenon or specific to the HSR reporter, we constructed two additional reporters with identical CDS and 3’ UTR as our HSR reporter, but regulated by different stress-dependent promoters. The first reporter includes the full-length promoter sequence of HSP12, a stress-responsive gene induced by a variety of stressors that includes an heat shock element (HSE) and seven STRE sequences. In a second reporter, we removed the HSE in this reporter by truncating 805 bp while preserving the STRE sequences. GCN3 knockdown decreased the translation efficiency of both versions of the pHSP12 reporters under arsenite stress (Figure 6B), while translation efficiency remained unchanged during AZC treatment (Figure 6B). No effects of GCN3 knockdown on translation efficiency were observed for any reporters in unstressed conditions (Figure 6B). These findings suggest a general role for GCN3 in increasing translation efficiency of stress-responsive genes under arsenite stress.

All experiments were performed with BY4741 yeast. The ssa2 knockout (yJP385) was made using NAT replacement of the SSA2 ORF and selected with nourseothricin. Yeast were grown at 30°C (unless otherwise indicated) in synthetic defined (SD) media with the appropriate dropout. SD media used for growth of yeast cultures contained: 2% w/v dextrose (Thermo Fisher Scientific, Waltham, MA), 13.4 g/L Yeast Nitrogen Base without Amino Acids (BD Biosciences, San Jose, CA), 0.03 g/L l-isoleucine (Sigma-Aldrich, St. Louis, MO), 0.15 g/L l-valine (Sigma-Aldrich), 0.04 g/L adenine hemisulfate (Sigma-Aldrich), 0.02 g/L l-arginine (Sigma-Aldrich), 0.03 g/L l-lysine (Sigma-Aldrich), 0.05 g/L l-phenylalanine (Sigma-Aldrich), 0.2 g/L l-threonine (Sigma-Aldrich), 0.03 g/L l-tyrosine (Sigma-Aldrich), 0.018 g/L l-histidine (Sigma-Aldrich), 0.09 g/L l-leucine (Sigma-Aldrich), 0.018 g/L l-methionine (Sigma-Aldrich), 0.036 g/L l-tryptophan (Sigma-Aldrich), and 0.018 g/L uracil (Sigma-Aldrich).

Approximately 12 sgRNAs were designed to target the 5’ end of each protein-coding and non-coding gene in the genome (between −750 and +50 nucleotides from the start of the ORF with priority toward sgRNAs closer to the ORF start) based on the Saccharomyces Genome Database (Cherry et al., 2012). A pool of oligonucleotides was ordered (CustomArray, Inc, Bothell, WA) containing these sequences and constant flanking sequences for the purposes of cloning (Supplementary file 10). A 2µ -URA plasmid was constructed from three DNA fragments using HiFi DNA assembly Master Mix (New England Biolabs, Ipswich, MA) for cloning. The first fragment was a diverse sample of inserts amplified via PCR using KAPA HiFi polymerase (Roche Basel, Switzerland), producing a sequence that contained random barcodes associated with a random sgRNAs from the library. This first fragment was generated in three steps. Step (1a) generation of a constant region using the parental plasmid as template, (1b) generation of random sgRNA sequences using the custom pool of oligonucleotides containing the sgRNA sequences, with flanking sequence to that of product of step 1a in its 5’ end. Step (2) PCR product using products of steps 1a and 1b to generate a constant region containing the random sgRNAs. Step (3) addition of random barcodes via PCR using product of step two as template, generating random barcodes in the 5’ end – which corresponds to the 3’ UTR of GFP – associated with random sgRNA sequences. The second fragment was the parent plasmid digested at two locations using the FastDigest enzymes BglI and BglII (Thermo Fisher, Sunnyvale, CA). The final fragment was another region of the plasmid amplified via PCR, which was necessary to connect the previous two pieces. These three fragments were assembled in a 100 µL reaction containing 700 ng of the first fragment, 600 ng of the second fragment, and 400 ng of the final fragment. After a 1 hr incubation at 50°C, the reaction was concentrated to 20 µL using a DNA clean and concentrate column (Zymo Research, Irvine, CA) and electroporated in five separate cuvettes containing 4 µL of the plasmid and 40 µL of ElectroTen-Blue electrocompetent cells (Agilent, Santa Clara, CA). After a 1 hr outgrowth in LB at 37°C, the plasmid was then grown in 1 L of LB containing carbenicillin for 2–3 days until it reached saturation. The plasmid was then purified via maxiprep (Qiagen, Hilden, Germany). A diagram of the plasmid library construction strategy and sequences used is shown in supplemental data files library_creation_diagram, library_primer_sequences, and library_sequences found in Supplementary file 11.

Two micrograms of plasmid was cut with the FastDigest PstI restriction enzyme (Thermo Fisher), and the 900 bp piece was gel purified using a gel DNA extraction kit (Zymo). This piece was then recirculized in a 20 µL reaction containing 60 ng of gel purified DNA and 1 µL T4 DNA Ligase (Thermo Fisher) in buffer. Cutting and recircularizing was necessary to increase efficiency of later sequencing steps for unknown reasons. The ligation product was then amplified with primers containing the appropriate overhangs for high-throughput sequencing and then gel purified. This piece was then sequenced through paired-end sequencing using custom primers on an Illumina (San Diego, CA) sequencing machine (Miseq or Nextseq).

Yeast were first transformed with a 2-micron pMET17-dCas9-MxiI/LEU2 plasmid (adapted from Gilbert et al., 2013) using standard methods. Five milliliters of this yeast was grown overnight in SD -Leu media to saturation. This yeast was then back diluted into a 1 L YPD culture to approximately ~0.1 OD₆₀₀. This culture was then grown at 30°C with shaking until the OD₆₀₀ reached ~0.6. The yeast were transformed with 600 ng of the library plasmid using standard procedures. However, instead of plating the yeast after the heat shock, they were grown in 1 L of SD -Ura -Leu for 3 days. The yeast were then either rediluted for experiments or frozen in 1 mL aliquots in media with 30% glycerol at −80°C for later use.

Yeast from a saturated library sample were grown overnight to log phase in SD -Ura/Leu. These yeast were then split into a number of cultures of at least 10 mL. Yeast were then either left untreated or stressed under the various conditions for the indicated times (Supplementary file 6). The yeast were then vacuum filtered onto nitrocellulose paper and flash frozen in liquid nitrogen. All replicates (see Supplementary file 1) are biological replicates, which came from separate transformations of the library into the yeast, separate treatments, and sample preparations. Some basal samples included technical replicates (independent sequencing sample preparation of the same yeast sample), and these were averaged into one biological replicate.

RNA was extracted through acid–phenol extraction. The frozen yeast pellets were first resuspended in 2.7 mL AE Buffer (50 mM sodium acetate pH 5.5, 10 mM EDTA), 0.3 mL 10% SDS, and 3 mL acid phenol:chloroform (Thermo Fisher). This mixture was heated to 65°C and shaken at 1400 rpm for 10 min followed by 5 min on ice. The solution was then spun at 13,000 g for 15 min. The supernatant was then added to 3 mL of chloroform and spun at 15,000 g for 5 min. The supernatant was then purified through a standard isopropanol extraction.

We then used NextFlex oligo-dT beads (PerkinElmer, Waltham, MA) to select polyadenylated mRNA from 100 µg of the extracted RNA. This sample was then subjected to DNAse treatment using Turbo DNase (Thermo Fisher) treatment and then purified using AMPure XP beads (Beckman Coulter, Indianapolis, IN). The resulting RNA was reverse transcribed using Multiscribe reverse transcriptase and a gene-specific primer (Thermo Fisher), and the cDNA was purified using alkaline lysis of the RNA followed by another round of AMPure bead purification.

To obtain the yeast DNA samples, frozen yeast was thawed and miniprepped using a Zymoprep Yeast Plasmid Miniprep Kit (Zymo Research).

The purified cDNA and the purified DNA were then amplified at the barcode region using primers with the appropriate Illumina sequencing ends and indexing barcodes to differentiate each sample. These resulting PCR products were then gel purified, mixed proportionately with each other, and then sequenced using a Illumina next-generation sequencing machine (MiSeq, NextSeq, or NovaSeq). The total number of barcode reads matching sgRNAs for each condition is listed in Supplementary file 9.

Basal neighborhood normalized scores (NNS) and interaction NNS were analyzed as pairs of conditions. For the case of the basal NNS, RNA levels and DNA levels were compared, whereas for the gene-stressors interaction NNS, stressor RNA was compared to unstressed RNA. While choosing the sgRNAs for library construction, an sgRNA mapped to two adjacent genes if the sgRNA targeted within 750 nucleotides upstream of both ORFs. For each pair of conditions, $A$ and $B$, and each barcode, $x$, we assign an NNS, $score(x,A,B)$, which captures how extreme (positive or negative) the read count is for barcode $x$ in condition $A$ versus condition $B$. These scores allow us to compare, in a principled fashion, the discrepancies between counts in different conditions, enabling a convenient aggregation to form a score for each gene or gene cluster. Crucially, the scores take into account (1) the higher relative amount of noise in small counts versus larger counts and (2) the fact that some conditions yield higher read counts across most barcodes.

We now describe and motivate our definition of $score(x,A,B)$. Let $A\left(x\right)$ denote the logarithm of the count of barcode $x$ under condition $A$, and let $B\left(x\right)$ denote the logarithms of the counts for barcode $x$ under conditions $B$, scaled by a factor $c_{A,B}$ that compensates for the possibility that some conditions yield systematically higher counts than others. After this scaling, for most barcodes, $x\text{'}$, it is the case that $A\left(x\text{'}\right)\approx B\left(x\text{'}\right)$.

Let $d_{A,B}\left(x\right)=\frac{A\left(x\right)-B\left(x\right)}{A\left(x\right)+B\left(x\right)}$ denote the normalized discrepancy in these values, and note that this quantity is bounded between ±1. The score is this discrepancy, $d_{A,B}\left(x\right)$ after the median value of this is subtracted, normalized by a factor $S_{A\left(x\right),B\left(x\right)}$ that takes into account how much variance is typical in this quantity, given the values of $A\left(x\right)$ and $B\left(x\right)$. Note that we expect $d_{A,B}\left(x\right)$ to have a larger variance due to random noise in the experiment and measurement when $A\left(x\right)$ and $B\left(x\right)$ correspond to small counts. This scaling factor $S_{A\left(x\right),B\left(x\right)}$ is defined to be an outlier-robust analog of the standard deviation of $d_{A,B}\left(x\text{'}\right)$ computed across the set of barcodes $x\text{'}$ for which ${A\left(x\text{'}\right)+B\left(x\text{'}\right)\approx A\left(x\right)+B\left(x\right).}_{}$

The following pseudocode formalizes the above calculation of the NNS for each barcode in each pair of conditions:

• Calculate barcode: NNS

• Input; Barcode x, read counts for all barcodes in conditions A and B, and scaling factor c(A,B).

• Output: score(x,A,B)

• Define parameters Cohort_Size = 10,000 and Clip_Param = 100.

• Let A(x) denote the logarithm of 1 + read_count for barcode x in condition A, and let B(x) denote the logarithm of 1 + read_count for barcode x in condition B scaled by a constant c(A,B) defined as the ratio of the median of the logarithms of the top 1000 counts in condition A, to the median of the logarithms of the top 1000 counts in condition B.

• Define the set of barcodes

where delta is the smallest value such that |X| ≥ Cohort_Size.

• For each x' in X, let d(x') = (A(x)-B(x))/(A(x)+B(x))

• Define S = truncatedStandardDeviation(d(X),Clip_Param), which is the standard deviation of the multiset of values of d(x') taken across all elements x' of X, after the largest and smallest Clip_Param values are removed. (E.g. if Clip_Param = 100, the largest and smallest 100 values are removed prior to computing the standard deviation.)

• Return score(x,A,B)=[d(x)-median(d(X))]/S.

• If the sum of the read counts in conditions A and B is less than 10, we discard that score, and do not use it when computing gene scores of GO scores.

Basal NNS (four replicates) and stress interaction NNS (two replicates) were computed as the arithmetic mean of replicates.

GO scores were calculated based on the above calculated NNS after saturation of the top and bottom 0.5% NNS (to prevent strong outliers from dominating a category score). For each GO category, the distribution of NNS for genes within that category and those of the entire distribution (scores for all genes) were compared using a Student’s t-test. The GO score was the natural log of the reciprocal of the resultant p-value from the Student’s t-test. The GO score sign indicates in whether the mean of the gene scores for that gene category was above or below zero.

Hierarchical clustering was performed using Cluster 3.0 (de Hoon et al., 2004), using uncentered correlation similarity metric and complete linkage clustering method.

Mitochondrially encoded genes were removed from supplementary data file one and subsequent analysis. Clustering and GO enrichment was performed on all genes except tRNAs and retrotransposons. Excluded gene list is provided in Supplementary file 12.

Fold changes for each gene were by averaging the log₂ fold changes for each barcode count and subtracting the median fold change of the distribution. Counts less than five were excluded.

MATLAB computer programs and raw data to compute gene scores and GO category enrichment are provided in Supplementary file 12.

Yeast were grown overnight from a saturated culture to log phase (OD < 0.6), such that they were growing in log phase for at least 14 hr. These yeast were then measured on an Accuri flow cytometer at 10,000 yeast per sample, under the fast flow setting. Each flow cytometry data point is a biological replicate, in which the steps starting from the wild-type (or in some cases knockout) yeast strain are repeated separately for each replicate, including plasmid transformations and treatment conditions.

RNA was extracted, and treated with DNase according to the protocol described above. One microgram of RNA was then reverse transcribed using Multiscribe reverse transcriptase (Thermo Fisher) and random hexamers in a 25 µL reaction. One microliter of the cDNA was then mixed with Luna qPCR master mix (New England Biolabs) and primers (to a final concentration 100 nM), which amplified either GFP or ACT1, the reference gene. These mixtures were measured in triplicate for each sample-primer combination with a ViiA 7 qPCR machine (Thermo Fisher), and the median C_T value was used in analysis. Each qPCR data point is a biological replicate, in which the steps starting from the wild-type yeast strain are repeated separately for each replicate, including plasmid transformations and treatment conditions.

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

We thank L Persson, C Sitron, J Giafaglione, J Park, S Alfonso, Drs. Z Davis, Z Jaafar, R Rohatgi, M Krasnow, J Weissman, L Gilbert, and C Jan for helpful suggestions throughout the project and in preparation of the manuscript. Funding was provided by the Lucille P. Markey Basic Biomedical Research Fellowship and NIH 5 T32 GM007276 to BDA, and NSF Awards 1813049 and 1704417 and ONR YIP award N00014-18-1-2295 to GV, and R01GM115968 to OB.

1. Received: March 30, 2020
2. Accepted: June 11, 2021
3. Version of Record published: July 5, 2021 (version 1)

© 2021, Alford 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.