Differential Hsp90-dependent gene expression is strain-specific and common among yeast strains

Summary Enhanced phenotypic diversity increases a population’s likelihood of surviving catastrophic conditions. Hsp90, an essential molecular chaperone and a central network hub in eukaryotes, has been observed to suppress or enhance the effects of genetic variation on phenotypic diversity in response to environmental cues. Because many Hsp90-interacting genes are involved in signaling transduction pathways and transcriptional regulation, we tested how common Hsp90-dependent differential gene expression is in natural populations. Many genes exhibited Hsp90-dependent strain-specific differential expression in five diverse yeast strains. We further identified transcription factors (TFs) potentially contributing to variable expression. We found that on Hsp90 inhibition or environmental stress, activities or abundances of Hsp90-dependent TFs varied among strains, resulting in differential strain-specific expression of their target genes, which consequently led to phenotypic diversity. We provide evidence that individual strains can readily display specific Hsp90-dependent gene expression, suggesting that the evolutionary impacts of Hsp90 are widespread in nature.


Hsp90-dependent gene expression varies among different yeast strains
Hsp90 differentially influences transcription factors among yeast strains Differential expression is correlated with phenotypic change on Hsp90 inhibition Hsp90-mediated strainspecific regulation manifests under environmental stress

INTRODUCTION
Living organisms constantly face the challenges of environmental change. Although most organisms are equipped with sophisticated regulatory mechanisms to acclimate to such change, devastating situations occasionally arise that can wipe out most individuals of a given population, leaving only a few individuals with diverse phenotypes to survive. Darwin's evolutionary theory states that heritable phenotypic variation provides the foundation on which natural selection operates, allowing organisms to evolve. 1 Nonetheless, maintaining phenotypic diversity may have fitness costs under typical growth conditions, especially if divergent phenotypes deviate from optimal fitness. 2 It is crucial for a population to maintain high fitness so that it can compete with other populations under normal conditions, but it must also harbor enough phenotypic diversity to allow it to survive under stressful conditions.
A balance between the requirements of population fitness and phenotypic diversity can probably be achieved through cryptic genetic variation, i.e., allowing variations that do not cause phenotypic change under normal conditions, but that can do so when environments or genetic backgrounds are altered. 3,4 Because of its neutral phenotypic effect in normal conditions, a significant proportion of cryptic genetic variation accumulates in the populations probably through genetic drift, similar to other types of neutral mutations. However, under certain environments or genetic backgrounds in which the buffering effect is compromised, some cryptic variants can also be quickly selected and enriched if they provide beneficial effects.
Multiple mechanisms have been proposed to explain the on/off phenotypic effect of cryptic genetic variation. 5,6 One of the mechanisms is the network hubs that work as capacitors to buffer phenotypic variation. Computational modeling has shown that complex gene networks can buffer the effects of genetic and environmental perturbations. 7,8 Moreover, deleting the hub genes often leads to phenotypic diversity even in an isogenic yeast population, confirming the hub genes' capacitor function. 9 On the other hand, the buffering effect of network hubs on genetic variation is not systematically tested in natural populations.
Among those network hubs, heat shock protein 90 (Hsp90) represents an interesting capacitor candidate. Hsp90 is a highly connected network hub in eukaryotes. 10,11 In the budding yeast, Saccharomyces cerevisiae, there are two Hsp90 paralogs, Hsc82 and Hsp82. Both genes are among the top 25 protein-protein interaction network hubs that interact with more than 10% of the yeast proteome and also among the top 10 genetic interaction network hubs that interact with more than 17% of yeast genes. 12 Moreover, the primary function of Hsp90 is to assist in protein folding or to maintain protein stability. 13 Therefore, it is suggested that Hsp90 may be able to directly alleviate the effects of some mutations that affect protein conformation or stability when they occur in Hsp90 substrates (also called Hsp90 clients). 4 The role of Hsp90 in protein folding and gene network suggests that Hsp90 may have the ability to buffer genetic variation directly or indirectly. More interestingly, because Hsp90 is often recruited to perform other tasks (such as unfolding denatured proteins) under stress conditions, the cryptic variants have a chance to be unbuffered in stress conditions and reveal their phenotypic effect. 14,15 Many Hsp90 clients are involved in signaling transduction pathways and transcriptional regulation, 16 further raising the possibility that the Hsp90 buffering (or capacitor) effect can have a strong influence on cell physiology and organism development.
The buffering effect of Hsp90 has been observed in a wide variety of organisms including yeast, plants, flies, fish, mice, and humans, 14,15,[17][18][19] suggesting that it is a conserved mechanism. However, with few exceptions, the respective studies have focused solely on a single pathway or phenotype. In addition, a recent study in yeast cell morphology indicated that the effects of newly occurred mutations were more likely to be enhanced by Hsp90, 20 suggesting the mutations that can be buffered by Hsp90 may not be as common as what we previously assumed.
Other than its buffering effect, Hsp90 can also function as a potentiator to enhance rather than suppress mutation effects, i.e., phenotypic diversity is revealed under normal conditions but is reduced when Hsp90 is compromised. 21 This potentiating effect of Hsp90 has been implicated in cancer development. [22][23][24] However, just as for Hsp90-dependent buffering, little is known about Hsp90-dependent potentiation in natural populations, apart from the results of two yeast studies. 15,20 S. cerevisiae provides an ideal system for studying Hsp90-dependent buffering and potentiating events at the population level. Whole genomes of divergent yeast strains isolated from various ecological environments and geographic locations have been sequenced, and their population structures are well understood. [25][26][27] Moreover, standing genetic variation may have experienced natural selection, which would purge out the ones with deleterious phenotypic effects. 20 Therefore, the buffered genetic variation observed in populations will provide additional information that cannot be extracted from newly occurred mutations in laboratory evolution experiments.
Because many potential Hsp90 clients are transcription factors (TFs), 28 we assumed that the buffering and potentiating effects of Hsp90 could be revealed by differences in the transcriptional profiles of multiple strains. We examined transcriptomes in response to compromised Hsp90 among five highly diverged yeast strains. We observed several hundred Hsp90-dependent genes that were differentially expressed in each strain. Using comparative genomics and gene expression profiling, we identified the candidate TFs associated with strain-specific gene expression variation. We experimentally demonstrate that differential activities or abundances of those TFs displaying strain-specific Hsp90-dependent regulation could further lead to phenotypic variation among strains. 2) Strain-pair-specific genes: genes whose expression shows a strong differece in Hsp90 dependnecy between two strains strain GdA Figure 1. Hsp90-dependent transcriptome analysis of five diverse yeast strains (A) Identification of differentially expressed genes and strain-pair-specific genes. After mRNA levels had been measured in the presence or absence of the Hsp90 inhibitor GdA, we calculated fold-change (FC) in expression between the two conditions (GdA+/GdAÀ) for individual genes in each strain. Genes were classified as differentially expressed if the adjusted pvalue was <0.05 and FC was >2 or <0.5. Next, by comparing two different strains (A and B), we defined a gene as a strain-pair-specific gene if the ratio of its fold-change in expression between the two strains was greater than 2 or less than 0.5. These strain-pair-specific genes represent genes that respond very differently to Hsp90 inhibition between the two strains. The strain-pair-specific genes were further classified into buffered or potentiated genes based on differences in expression under normal and Hsp90-compromised conditions. Buffered genes  (Table S1). To ensure that the variation in expression profiles was not caused by different sensitivities to GdA in different strains, we assayed the Hsp90 activity by measuring the kinase activity of v-Src, a well-known Hsp90 client, 30 with or without GdA treatments ( Figure S2C and STAR Methods). The v-Src kinase activity was reduced to the background level in all tested strains on GdA treatment. Although we cannot rule out that the differences between the strains are below the detection limit of v-Src assays, this result indicates that the sensitivity difference to GdA is unlikely the primary reason for the observed expression differences.
We also noticed that the ML strain grew more slowly than the other strains when treated with 50 mM GdA (Figure S2D), raising the possibility that the expression difference may result from the secondary effect of slow cell growth. Therefore, we performed an additional RNA-seq using the ML cells treated with 25 mM GdA. Under this condition, the ML strain did not exhibit such slow growth, but its transcriptome profile is still more similar to the ML strain treated with 50 mM GdA than the other strains (Spearman's correlation, r = 0.83, 0.81, 0.74, 0.72, and 0.59 for 25 mM GdA-treated ML compared to 50 mM GdA-treated ML, NA, WA, SK, and LAB, respectively). Next, we compared the transcriptome profiles of 50 mM-GdA treatment with other stress conditions that slow down cell growth (Table S2). The correlation between the GdA treatment and other stress conditions is much lower than the correlation between different GdA-treated strains. These data suggest that the observed regulatory effects are not the result of slow growth in GdA-treated cells.
Moreover, we examined the expression of 133 genes known to be up-regulated under stress conditions (Table S3). 31 GdA-treated SK cells had the highest number of up-regulated stress-response genes (27/ 133 = 20%), raising the possibility that GdA-treated SK cells were slightly under stress. On the other hand, most of the stress response genes were not induced. We also compared our differentially expressed genes with the targets of Hsf1, a general TF of heat shock response. 32 The Hsf1 targets were not significantly enriched (Table S4). Together, our data suggest that the GdA treatment did not trigger a general stress response and observed gene expression patterns were specific to Hsp90 inhibition.
Hundreds to more than one thousand genes have changed their expression levels on GdA treatment (fold change >2 or <0.5, multiple testing correction adjusted pvalue <0.05, Figures 1A, S1B, S1C, and S3, and Table S5). Although the numbers of differentially expressed genes varied drastically among the strains, principal component analysis indicated that the major expression difference resulted from the environmental effect (i.e., with or without GdA) rather than the genetic background ( Figure S4A). Of interest, the lab strain showed the least differentially expressed genes, probably an outcome of long-time laboratory adaptation (Gu and Steinmetz, 2005).
Among the Hsp90-dependent differentially expressed genes, 298 genes were shared by at least four strains ( Figure S3). Gene ontology (GO) analysis showed that seven biological processes were enriched (multiple testing correction adjusted pvalue <0.05) in these common Hsp90-dependent genes (Table S6), which likely represent the general physiological effects of Hsp90 or ancestral cryptic genetic variation shared by multiple strains. For instance, Hsp90 and its co-chaperones have been shown to influence mitochondrial protein import, with this function being conserved in both yeast and mammalian cells, 33,34 perhaps explaining GO enrichment for cellular respiration and response to oxidative stress. We also examined the GO enrichment of differentially expressed genes in each strain to find the commonly shared and strain-specific responsive pathways. Five enriched GO terms were shared by more than three strains ( Figure 1B) that completely overlapped with the seven GO terms derived from the common Hsp90-dependent genes (Table S6). Nevertheless, even for commonly shared GO terms, we observed different levels of enrichment among strains ( Figure 1). These data indicate that the general effect of Hsp90 on individual cellular processes can be further adjusted depending on genetic background. Continued show a greater difference in expression between two strains under the Hsp90-compromised condition, whereas potentiated genes show a greater difference in expression between two strains under normal conditions. (B) GO analysis of genes differentially expressed in response to Hsp90 inhibition (+GdA/ÀGdA) for each strain. Only enriched GO terms have been plotted and circle size represents adjusted pvalues. GO terms enriched in more than three strains were classified as a common response GO term (shown in red).
(C) GO analysis of the strain-pair-specific genes for all strain pairs. Those genes exhibited more than a 2-fold difference in Hsp90-dependent expression between two strains. See also Figures S1, S3-S5, and Tables S3, S5 iScience Article In addition to common Hsp90 responses, we observed many strain-specific responses that were only enriched in one or two strains ( Figure 1B), suggesting that different strains can display unique Hsp90-dependent expression patterns. We used them to further dissect the underlying mechanisms and possible physiological effects.
Both Hsp90-buffering and potentiated events are observed in yeast populations Because we had already observed many strain-specific responses in gene expression, next we searched for differentially expressed genes that exhibited apparent differences in Hsp90 dependency between strain pairs. To do this, we compared fold-changes in mRNA abundance in cells of two given strains subjected to either of two conditions (i.e., GdA treatment or lacking GdA treatment). If a gene presented a between-strain ratio greater than 2 or less than 0.5, that gene was defined as a strain-pair-specific gene ( Figures 1A, S5, and Table S7). In later experiments, we used this strain-pair-specific gene information to facilitate the identification of the TFs that are influenced by Hsp90 in a strain-specific manner.
When we performed GO enrichment analysis on strain-pair-specific genes, we found that the enriched GO terms strongly overlapped with those from differentially expressed genes ( Figures 1B and 1C). This outcome indicates that many differentially expressed genes indeed exhibit significant differences in their Hsp90-dependencies among strains. The identified strain-pair-specific genes could be further divided into buffered and potentiated groups depending on whether the differential expression between two strains under GdA-treated conditions was greater or smaller than that under untreated conditions, respectively ( Figures 1A, S4B and Tables S7 and S8).
Genetic buffering is known to increase the complexity of the phenotype. [4][5][6] Using the transcriptome profiles to construct expression distance trees (see STAR Methods), it shows that the expression tree in the normal condition (ÀGdA) is highly similar to the genetic distance tree (Spearman's correlation, r = 0.71, Figures S1A and S1D). Again, the LAB and ML strains show the longest distance. The NA and SK strains have the second smallest distance (slightly longer than the distance between NA and WA). On the other hand, the expression tree in the Hsp90-inhibition condition (+GdA) is very different from the genetic distance tree (Spearman's correlation, r = À0.11, Figure S1E). These data suggest that the transcriptomic profiles diverge over time, similar to the genetic distance. However, the Hsp90 inhibition adds another layer of complexity to the divergence as predicted by the capacitor hypothesis. 4 Identification of TFs that contribute to variable Hsp90-dependent gene expression One of the major regulators of gene expression are TFs, which recognize specific DNA sequences (also called transcription factor binding sites or TFBSs) within promoter regions and then activate/suppress expression of the target genes. 35 Because genes involved in similar cellular processes may be regulated by the same TFs, we decided to examine if the strain-specific effects of Hsp90 on differential gene expression are mediated through specific TFs.
First, we analyzed the distribution of TFBSs among different strains. Although the information on TFBS location from chromatin immunoprecipitation (ChIP)-on-chip or ChIP-seq experiments are widely used in analyses of TF-based regulation, such experiments in yeast have typically been done under normal growth conditions and in laboratory strains. 32,[36][37][38] Inhibition of Hsp90 triggers physiological changes and leads to altered transcriptional regulation. 39 Moreover, several studies have shown that TF-bound targets can be dramatically different among strains and that these differences in TF binding are often associated with levels of sequence variation. 40,41 Thus, current experimentally validated TFBS data do not provide sufficient information to reveal the impacts of Hsp90 inhibition or how they differ among different strains. To gain a global picture of how TFBSs are distributed in our selected strains, we used a bioinformatics pipeline to predict potential TFBS via a genome-wide promoter scan based on known binding motifs ( Figure S6, also see STAR Methods). 42 Among all predicted TFBSs, only the ones that were conserved across strains (96%) were used for further analyses (Table S9). This analysis pipeline was highly sensitive and, inevitably, generated some false-positive signals. We used the ChIP-seq data from a previous study to validate our prediction. 32 On average, about 58% of them overlapped with our predicted candidates (Table S10, see also STAR Methods), indicating that our analysis pipeline is reliable. Moreover, we found that many of our predicted results could be validated by further experiments (see later sections).

OPEN ACCESS
iScience 26, 106635, May 19, 2023 5 iScience Article After constructing the TFBS map, we grouped genes sharing the same TFs, and then searched for TFs whose downstream targets exhibited differential Hsp90-dependent expression among strains. We hypothesized that if a TF or its upstream signal transduction components contain strain-specific genetic variations that are buffered or potentiated by Hsp90 directly or indirectly, such events may be revealed in the expression profiles of the TF target genes when Hsp90 is compromised. Among 187 analyzed TFs, we identified 21 TF whose targets showed significant differences in expression profile among our five tested strains (Figures S6, S7, and Table S11). Of interest, Hsp90-interacting TFs are significantly enriched in our candidate list (7 out of 21 in our list versus 23 out of 187 in the whole genome, hypergeometric test, p = 0.0011). Previously, our identification of strain-pair-specific genes allowed us to detect the differentially expressed genes that exhibited strong strain-specific Hsp90-dependencies ( Figure 1 and Table S7). We used this information as another criterion to further refine our list of TF candidates. If strain-pair-specific genes are enriched among the targets of a specific TF, it further supports that this TF may have strong impacts on strain-specific differential gene expression. Eleven TFs were observed to exhibit such enrichment in at least one strain pair ( Figure 2, Tables S11 and S12). Moreover, some of those TFs shared similar enrichment patterns (Figures 2 and S8). Next, we selected six of the eleven TFs (Msn2, Msn4, Com2, Rap1, Dot6, and Tod6), representing different expression patterns among strain pairs, for experimental validation.

Hsp90 inhibition leads to strain-specific upregulation of Msn4
Msn2 and Msn4 physically interact with Hsp90, 43 and they are functionally redundant paralogs that activate genes involved in stress responses. 44,45 We found that when Hsp90 was compromised, target genes of Msn2 and Msn4 in the SK and NA strains were significantly upregulated relative to in the ML strain (Figure 3A). However, because both Msn2 and Msn4 can bind to the same stress-responsive element (a DNA motif containing AGGGG) and share a large proportion of their targets, 46 it raises the possibility that only one of them is responsible for the observed changes in expression. Thus, the identification of The targets of these candidate TFs show significant differences in between-strain Hsp90-dependent expression. Dot size represents the Bonferroni-adjusted pvalue from the Wilcoxon rank-sum test (Table S11). TF names are labeled in red on the yaxis if they are Hsp90 clients. See also Figures S6-S8, and Tables S11 and S12.  To further establish the effect of Hsp90 on Msn2 and Msn4, we measured protein abundances of Msn2 and Msn4 in response to Hsp90 inhibition (GdA treatment). We found that, in the SK strain in which Msn2 and Msn4 targets were upregulated in Hsp90-compromised cells, Msn4 was also highly induced in both mRNA and protein levels ( Figure 3B and Table S3). In contrast, the ML strain maintained similar levels of Msn4 under GdA-treated and -untreated conditions. Moreover, Msn2 abundance was not affected by Hsp90 inhibition in either strain ( Figure 3B and Table S3). Because the functions of Msn2 and Msn4 are also regulated by post-translational modifications and nuclear translocation, 47 we performed one-hybrid reporter assays to directly measure the transcriptional activities of Msn2 and Msn4 (see STAR Methods). Consistent with our protein abundance data, Msn4 activity was increased in both the 25 mM and 50 mM GdA-treated SK strain relative to that of the GdA-treated ML strain ( Figures 3C and S9A). Msn2 only showed mild Hsp90-dependent differences in activity (<1.5-fold) between strains. To ensure the observed effect is really Hsp90dependent, we repeated the experiment using a different Hsp90 inhibitor, radicicol, and observed a similar SK-specific Msn4 induction ( Figure S9A). Together, these data suggest that Msn4 is the dominant factor contributing to the differential Hsp90-dependent expression of Msn2 and Msn4 targets. Besides, strong induction of Msn4 in GdA-treated SK cells might be the reason why more stress-response genes were up-regulated in the SK strain (Table S3).
Hsp90 affects Com2 activity without changing its protein abundance in the NA strain Com2 shares a similar core motif (AGGGG) with Msn2 and Msn4 (Table S11), but its function is not well characterized. 42,48 In our reporter enrichment analysis, Com2 exhibited a different pattern to that of Msn2/4 ( Figures 2 and S8), suggesting that Com2 targets were affected by Hsp90 via different regulatory mechanisms. Basal protein levels of Com2 varied among strains ( Figure 3E), but the alterations in protein abundance because of Hsp90 inhibition did not differ significantly among strains. This outcome indicates that specific upregulation of Com2 targets in the NA strain ( Figure 3D) was probably not mediated by altered Com2 abundance ( Figure 3E). However, Com2 transcriptional activity of the NA strain, but not the LAB or ML strains, increased significantly under GdA treatment ( Figure 3F). Consistent with the activity assay, the target genes of Com2 (YEF1, FMP33 and CYC7) also showed strong NA-specific induction in their mRNA levels under radicicol treatment when measured by quantitative PCR (q-PCR) ( Figure S9B and Table S13). Thus, post-translational regulation is probably involved in the Hsp90-dependent effect of Com2 observed here, unlike the case of Msn4.

Hsp90 influences the expression of ribosomal genes through multiple TFs
Because we noted that several enriched GO terms of the differentially expressed genes were related to ribosomal function ( Figures 1B and 1C), we next focused on Rap1, an essential TF involved in diverse cellular processes, including chromatin silencing, telomere length maintenance and transcriptional activation of (C) Msn4 activity is strongly enhanced in the SK strain, but not in the ML strain, upon Hsp90 inhibition (two-sided unpaired t-test with Welch's correction, p < 0.0001 for Msn4, p = 0.0618 for Msn2). The transcriptional activities of Msn2/4-lexA fusion proteins were measured by the one-hybrid assay in which the reporter (b-galactosidase) is transcriptionally regulated by the lexA fusion protein and b-galactosidase activity is used to represent relative TF activity.
(D) Expression profiles of Com2 targets reveal significant differences in Hsp90 dependency between the LAB, ML and NA strains (two-sided Wilcoxon rank-sum test) (  Figure S9, and Tables S11, S12, and S13.  Figure 4A).   (E) Dot6 and Tod6 are involved in Hsp90-dependent regulation of expression in the SK strain. A reporter construct lacking the Rap1 binding site (Sk-RPL40B-Rap1-BSD) was used in this experiment to measure the Rap1-independent effect. Promoter activity became almost Hsp90-independent (fold change was approached to 1) when either DOT6 or TOD6 was deleted from the SK strain (two-sided unpaired t-test with Welch's correction, p = 0.0056 between wild type and dot6Dmutant, p = 0.0026 between wild type and tod6D mutant).
(F) A model describing how Hsp90 influences the expression of ribosomal genes via Rap1, Dot6, and Tod6. The solid line represents physical interactions or direct transcriptional regulation, and the dashed lines represent the regulation indirectly inferred from geneticdata. *, adjusted pvalue <0.05; **, adjusted pvalue <0.01; ***, adjusted pvalue <0.001. See also Figure S9, and Tables S11, S12, and S13. We examined two known Rap1-regulated ribosomal genes, RPL40B and RPS26B, using q-PCR to determine if differential expression patterns were similar for individual genes (Table S13). Consistent with the overall Rap1 target data ( Figure 4A), both of these genes exhibited strong, intermediate, or weak Hsp90-dependent effects in the SK, WA, and LAB strains, respectively ( Figures 4B and S9C). A similar pattern of Hsp90 dependency was also observed for Rap1 protein abundance ( Figure 4C). This data is consistent with a previous observation that in an Hsp90 temperature-sensitive mutant, the Rap1 protein abundance was reduced on heat inactivation of Hsp90. 52 Together, these results indicate that Hsp90 can influence the same TF to different extents depending on strain background. We also constructed an SK-RPL40B promoter-lacZ reporter gene to investigate to what extent Rap1 contributes to RPL40B expression. When we deleted the Rap1 binding site (CCATCCGTGCCTC; located 273 base pairs upstream of the translation start site) from the SK-RPL40B promoter, we found that b-galactosidase activity was drastically reduced in the SK strain ( Figure 4D), confirming that Rap1 is a crucial activator of ribosomal gene expression.
Our list of six TFs contains a pair of paralogs, Dot6 and Tod6, which inhibit ribosomal gene expression in response to environmental change. 53 As for Msn2 and Msn4, Dot6 and Tod6 bind to similar binding motifs and share a significant proportion of their targets. However, Dot6 and Tod6 are controlled separately by the TORC1 and PKA pathways, respectively, indicating that these two TFs have diverged and are controlled by different upstream signals in response to diverse environments. 53 Besides, Hsp90 physically interacts with Dot6 but not Tod6. 12 To investigate the contribution of Dot6 and Tod6 to differential Hsp90-dependent gene expression, we performed promoter activity assays on the SK-RPL40B promoter-lacZ reporter in the SK strain. Our TFBS analysis indicated that the RPL40B promoter contains both Rap1 and Dot6/Tod6 binding sites (Table S11). We used a reporter construct lacking the Rap1 binding site to directly measure the effects of Dot6 and Tod6. When we deleted DOT6 or TOD6, expression of the SK-RPL40B promoter-lacZ reporter became almost Hsp90-independent in the SK strain ( Figure 4E). Thus, unlike Msn2 and Msn4, both Dot6 and Tod6 play a role in strain-specific Hsp90 dependency. Together, these findings reveal that differential Hsp90-dependent expression of ribosomal genes is controlled by a complex network of TFs that includes the transcription activator Rap1, as well as the repressors Dot6 and Tod6 ( Figure 4F).

Hsp90-dependent expression profiles are correlated with cellular responses to environmental change
Our experiments validated that Hsp90 exerts differential effects on several TFs depending on the strain background, leading to strain-specific differential expression of many target genes. However, such differences will only impact population evolution if they result in physiological changes in the cell. Msn2 and Msn4 are known to be involved in oxidative stress responses. 44 Because we observed Hsp90-dependent upregulation of Msn4 in the SK strain (Figure 3), we tested if such a change could lead to fitness differences between the SK and ML strains under oxidative stress. In the absence of H 2 O 2 , inhibition of Hsp90 had only a mild effect on the growth of ML and SK strains. However, when we challenged cells with H 2 O 2 , the ML strain displayed greater sensitivity to that oxidative stress than the SK strain ( Figure 5A).
Given that expression of ribosomal genes controlled by Rap1, Dot6, and Tod6 also exhibited different levels of Hsp90 dependency between the LAB and SK strains (Figure 4), we directly tested if these two strains had different resistances to rapamycin, a drug that represses the expression of ribosomal genes by inhibiting the TOR pathway. 54 Consistent with levels of Hsp90 dependency among its impacted genes, the SK strain indeed showed lower resistance to rapamycin when Hsp90 was compromised ( Figure 5B). Our results demonstrate a direct link between strain-specific differential expression and altered cell physiology, suggesting that the observed buffering or potentiating effects of Hsp90 have the potential to influence population adaptability.

Heat stress triggers yeast cells to manifest strain-specific phenotypic variation
The capacitor hypothesis predicts that Hsp90-buffered genetic variation allows cells to manifest phenotypic variation when faced with environmental challenges. 4 We already revealed several strain-specific phenotypic variations on inhibiting Hsp90 by GdA treatment. Could such variations also be manifested when cells encounter environmental stress, as predicted by the capacitor hypothesis? Previous studies have suggested that the massive amounts of misfolded proteins induced by acute environmental stress could titrate out the cytosolic activity of Hsp90 and cause Hsp90-compromised phenotypes. 14 iScience Article First, we measured Msn4 protein abundance in the conditions with or without heat stress. Although heat stress induces Msn4 activation, 56 we found that the ML strain did not significantly increase Msn4 protein abundance after 1 h of heat stress, whereas SK showed a strong induction ( Figure 6A). Next, we examined Msn4 transcriptional activity via one-hybrid reporter assays with or without heat stress. Consistent with the protein abundance data, we observed much higher induction in the SK strain than in the ML strain (Figure 6B). We also measured the expression of an Msn4 target, CTT1 (Table S13), using the CTT1 promoter-lacZ reporter. Consistent with the Msn4 activity data, CTT1 expression was much higher in the SK strain on heat shock ( Figure 6C). Moreover, when we deleted the Msn4 binding site in the CTT1 promoter (AAGGGGC; located 364 base pairs upstream of the translation start site), the induction of CTT1 was significantly reduced in the SK strain but not in the ML strain ( Figure S9D). These data provided evidence that Msn4 was the major contributor to the heat-induced strain-specific induction of CTT1. We also measured Rap1 protein abundance and used the RPL40B promoter-lacZ reporter to investigate how heat stress affected the expression of ribosomal genes. Similar to GdA treatment (Figure 4), Rap1 protein abundance was reduced in WA and SK strains but not changed in the LAB strain. Expression of the RPL40B promoter was also significantly reduced in the SK strain after heat stress ( Figures 6D and 6E).
Through our GdA-treatment experiments, we showed that the observed strain-specific patterns of differential expression are correlated with fitness differences under changing environments ( Figure 5). Extending that analysis, we challenged cells under a heat stress environment with the additional challenge of either oxidative stress or rapamycin treatment. Similar to our findings from GdA treatment, under heat stress the ML strain exhibited greater sensitivity to H 2 O 2 , whereas the LAB strain grew better than the SK strain in the presence of rapamycin (Figures 6F and 6G). Together, the phenotypes observed under heat stress were consistent with the ones revealed on Hsp90 inhibition (Figures 3, 4, 5, and 6), suggesting that heat stress can trigger cells to manifest strain-specific phenotypic variation, just as treatment with a specific Hsp90 inhibitor can. Nonetheless, our data could not completely rule out the possibility that heat stress triggered other Hsp90-independent responses, leading to observed strain-specific phenotypic differences. More experiments are required to determine whether Hsp90-dependent effects are the major contributor to heat-induced protection.

DISCUSSION
Network hubs have been suggested to buffer genetic and non-genetic variation. 7,8,57 However, the underlying molecular mechanisms are not well understood because many genes are involved in a Relative growth rate (+GdA/-GdA) * * *** ** Figure 5. Hsp90-dependent expression profiles are correlated with cellular responses to stress (A) The ML strain is more sensitive to the oxidative stress induced by H 2 O 2 treatments than the SK strain when Hsp90 is compromised (two-sided unpaired t-test with Welch's correction between strains, p = 0.0015 for YPD condition, p = 0.0003 for the oxidative condition). This result is consistent with the observation that Msn4 is highly induced by Hsp90 inhibition only in the SK strain.
(B) The LAB strain is more resistant to the TORC1 inhibitor, rapamycin, than the SK strain on Hsp90 inhibition (two-sided unpaired t-test with Welch's correction between strains, p = 0.0008 for YPD, p = 0.0075 for the rapamycin condition). This result is consistent with the observation that the SK strain shows a stronger Hsp90-dependent reduction in the expression of ribosomal genes than the LAB strain. Each dot represents a biological repeat data. **, adjusted pvalue <0.01; ***, adjusted pvalue <0.001.  Figure 6. Heat stress triggers yeast cells to manifest Hsp90-dependent strain-specific phenotypic variation (A) SK-Msn4, but not ML-Msn4, is induced after heat stress (two-sided unpaired t-test with Welch's correction between normal and heat stress conditions). Msn4 protein abundance was measured by western blot after shock (37 C) for 0 and 1 h. Individual dots in the bottom panel represent the quantified fold change from each repeat. (B) Msn4 activity is induced under heat stress in the SK strain, but not in the ML stain (two-sided unpaired t-test with Welch's correction between strains, p = 0.0006). Msn4 activity was measured by one-hybrid assay before and after heat shock (37 C) for 1 h. (C) Expression of CTT1, a target of Msn4, shows an induction pattern similar to that of Msn4 activity after heat shock (37 C) for 1 h (two-sided unpaired t-test with Welch's correction between strains, p = 0.0007). iScience 26, 106635, May 19, 2023 iScience Article hub-dependent subnetwork. Here, we showed that the network hub, Hsp90, can buffer (or potentiate) the gene expression differences via regulating Hsp90-dependent TFs in a strain-specific manner. Because the Hsp90-dependent network involves more than 1/3 of the yeast genome, 12 the potential buffering or potentiation targets of Hsp90 may cover a vast gene set. Instead of mapping them, we focus on understanding how the Hsp90-dependent genetic variation leads to phenotypic variation. By analyzing differential gene expression patterns among five diverse yeast strains, we identified 11 TF that are subject to Hsp90-dependent and strain-specific regulation (Figure 2). We performed experiments to further confirm the involvement of five TFs in the observed patterns of gene expression and predicted impacts on cellular physiology. Together with various enriched GO terms observed in differentially expressed genes, our data suggest that strain-specific Hsp90-dependent regulation is prevalent in natural yeast populations.

OPEN ACCESS
The TFs identified here can indicate which pathways/subnetworks have the genes carrying cryptic genetic variation. As a network hub and also a molecular chaperone, Hsp90 may regulate cryptic genetic variation through two different mechanisms. 58 In the direct buffering (or potentiating) model, the cryptic variation resides in the Hsp90 client and affects its folding or interacting ability. Hsp90 directly interacts with the variant-carrying protein to maintain its function or protein stability. Some disease-related mutations and oncogene mutants likely belong to this category. 18,22,23 If our candidate TFs physically interact with Hsp90 and carry such genetic variation, they might influence expression profiles through the direct buffering (or potentiating) mechanism ( Figure 7A). The indirect buffering (or potentiating) model suggests that the variant-carrying protein is not itself an Hsp90 client, but is a component of the Hsp90-dependent subnetwork. When Hsp90 is inhibited, the compromised client synergistically enhances the mutant effect, leading to a detectable phenotype. Hsp90-buffered gene expression arising from repression of endogenous retroviruses in mice represents a well-documented example of indirect buffering. 19 In our case, the casual mutations might occur in the TF or signal transduction gene downstream of an Hsp90 client and only partially compromise its function so normal gene expression is maintained. However, when Hsp90 inhibition reduces the function or stability of the client protein, the defect of the downstream mutations is multiplied, leading to a strain-specific expression pattern ( Figure 7B). Through indirect buffering (or potentiation), Hsp90 can impact a much broader range of proteins even though Hsp90 clients only represent around 16% of the yeast proteome. 12 Seven out of the 11 high-confident Hsp90-dependent TFs do not directly interact with Hsp90. These nonclient TFs can be regulated by several possible mechanisms. First, the non-client TFs may be the direct targets of the Hsp90-dependent client TFs. By examining the promoter regions of the seven non-client TFs, we identified the binding motifs of the Hsp90 client TFs (DOT6, MSN4, and GIS1) in five of the non-client TFs (Table S13). However, the general expression patterns of the client TF targets are not consistent with that of the motif-containing non-client TFs and their targets except for the DOT6-TOD6 pair (Table S14). Because Dot6 and Tod6 share many common targets, it remains unclear how significant such regulation is in this case. Second, the non-client TFs may be regulated by Hsp90 clients post-transcriptionally. As revealed by the phosphoGRID database, 59 all of our non-client TFs contain multiple post-translational modification sites although the modification enzymes are not identified in most cases (Table S14). Many protein modifiers are known to be Hsp90 clients and their activities are controlled by Hsp90. 60,61 In addition, our results have shown that one of our Hsp90-dependent TFs, Com2, is probably regulated by Figure 6. Continued (D) Rap1 protein abundance is significantly reduced after heat stress in SK, ML, and WA strains (Ordinary one-way ANOVA followed by Turkey's multiple comparison test). Rap1 protein level was measured by Western blot after heat shock (37 C) for 0, 1, and 4 h. Individual dots in the bottom panel represent the quantified fold change from each repeat. (E) Expression of RPL40B is reduced more in the SK strain than in the LAB strain under heat stress (two-sided unpaired t-test with Welch's correction between strains, p = 0.0162). b-Galactosidase assays were performed before and after heat shock for 1 h using the strains carrying the RPL40B promoter-lacZ reporter. (F) The SK strain grows better than the ML strain in the presence of oxidative stress at high temperatures (two-sided unpaired t-test with Welch's correction between strains, p = 0.0205 for YPD, p < 0.0001 for the H 2 O 2 condition). (G) The LAB strain grows better than the SK strain in the presence of rapamycin at high temperatures (two-sided unpaired t-test with Welch's correction between strains, p = 0.0319 for YPD, p = 0.0003 for the rapamycin condition). In F and G, we chose 35 C for high-temperature growth because some tested strains grew very slowly at 37 C. Each dot represents a biological repeat data. *, adjusted pvalue <0.05; **; adjusted pvalue <0.01, ***, adjusted pvalue <0.001. See also Figure S9, and iScience Article post-translational modifications ( Figure 3). We examined the one-degree interactors of the non-client TFs and identified 40 Hsp90 client proteins, including TF cofactors, Hsp90 cochaperones, protein modifiers, and ribosomal subunits (Table S14). Hsp90 may regulate the non-client TFs through these interacting Hsp90 client proteins.
Recent studies evidence that transcriptional rewiring is quite common during yeast evolution. [62][63][64] Of interest, several of our identified TFs were involved in these rewiring events. For instance, the regulation of ribosomal genes was dramatically rewired from Rap1 in S. cerevisiae to Tbf1 in Candida albicans. 65,66 Moreover, it has been suggested that Dot6 and Tod6 targeting was rewired among yeast species based on motif changes. 67 With regard to the regulatory pathway for responses to oxidative stress, the requisite TFs were rewired from Msn2/Msn4 in S. cerevisiae to Fkh2 in C. albicans. 66 It has been suggested that transcriptional rewiring often occurs through an intermediate state in which regulatory changes are selectively neutral to the cells. 62 Because Hsp90-buffered genetic variation allows cells to accumulate neutral regulatory changes, it raises the interesting possibility that Hsp90 may play a role in facilitating transcriptional rewiring.
Altered transcriptional regulation is thought to play a crucial role in the adaptive evolution of various organisms. [68][69][70] For example, a comparison of humans with other primates has revealed that TFs are enriched among genes exhibiting a human-specific increase in expression. 71 In threespine stickleback fish, the adaptive pelvic reduction observed in different natural populations is caused by recurrent mutations in Pitx1, a TF required for normal hindlimb development in vertebrates. 72 Repeated mutations in another TF, optix, were also found to be responsible for convergent evolution of wing pattern mimicry among different butterfly species. 73 In plants, genetic variations in two MADS-box TFs, SOC1 and FUL, alter plant growth form A B Figure 7. Two models showing how Hsp90 can buffer a variant-carrying protein in a transcriptional regulatory pathway (A) In the direct buffering model, the variant-carrying TF is an Hsp90 client and its function is stabilized by Hsp90. When the cell encounters stress, Hsp90 is recruited to clear misfolded proteins caused by stress, so the TF function is no longer buffered, allowing differential gene expression. (B) In the indirect buffering model, the TF function is regulated by protein X, which is an Hsp90 client. When X is fully functional, the variant-carrying TF or protein Y does not show obvious defects. However, when X is compromised after titration of Hsp90 (by stress-induced misfolded proteins), the combined effects of compromised X and variant-carrying TF (or Y) result in differential gene expression. iScience Article (annual or perennial) and longevity. 74 Adaptive evolution mediated by TF variation is equally common in microorganisms. In yeast, altered expression or activity of TFs has been shown to contribute to populationor species-level differences in morphology, environmental responses, or lifecycle. [75][76][77][78][79][80][81][82] Most of the observed changes in TFs and gene expression we report in our study are probably not the outcome of adaptive evolution. However, we have shown that TF activity or abundance is regulated by Hsp90 in a strain-specific manner, leading to Hsp90-dependent expression variations and also detectable phenotypic differences when the function of Hsp90 is perturbed by stress. These Hsp90-dependent variations provide a large repertoire of material on which natural selection can act when cells are challenged by changing environments.

Limitations of the study
Although we have identified several candidate pathways and TFs that exhibit different Hsp90-dependent behaviors among tested strains, we do not know whether the strain-specific Hsp90-dependent variants occur on those TFs or their upstream genes. More experiments are needed to map those variants and understand how Hsp90 exerts its buffering and potentiating effects.

STAR+METHODS
Detailed methods are provided in the online version of this paper and include the following:

ACKNOWLEDGMENTS
We thank Gianni Liti for the natural yeast strains and members of the Leu lab for helpful discussion and comments on the manuscript. We also thank John O'Brien for manuscript editing, and the IMB Genomics core and Academia Sinica Sequencing core for sequencing services. This work was supported by Academia Sinica of Taiwan

Lead contact
Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Jun-Yi Leu (jleu@imb.sinica.edu.tw).

Materials availability
Yeast strains described in this study are available from the lead contact.
This study did not generate new unique reagents.
Data and code availability

EXPERIMENTAL MODEL AND SUBJECT DETAILS Yeast strains
All S. cerevisiae strains used in this study are shown in Table S15. Primer pairs and the template DNA used for construction are listed in Table S16. We used polymerase chain reaction (PCR) to amplify DNA fragments where necessary, and PCR products were purified using a PCR cleanup kit or gel extraction kit (K-3037, Bioneer, Oakland, CA, USA). For genome editing, DNA fragments containing the selection marker flanked by 5'and 3' gene-specific homologous regions were transformed into cells to replace the endogenous region through homologous recombination. For TAP tagging, the 5 0 homologous region was amplified by the primers gene-F1 and gene-F1-TAP-R1, and the 3 0 homologous region was amplified by the primers gene-TAP-F3 and gene-R3. The TAP and HIS selection markers were amplified from the TAP collection. 98 Then, a second round of PCR was used to fuse these fragments. The final products were purified and transformed into yeasts using the lithium chloride method. 99 Strains in which final products had been correctly integrated were selected under a selection plate and checked by PCR and sequencing.
For the one-hybrid assay, the reporter system (including lacZ driven by the promoter with eight operator sites of lexA) was derived from the pSH18-34 plasmid. 100 To increase basal promoter activity for testing repressor activity, five operator sites for Gcn4-binding were inserted downstream of the lexA operators. The entire cassette was inserted into the HIS3 locus using the URA3 marker. The lexA-tagging for TFs was achieved using pSFS4A-lexA plasmid, a marker-free system derived from the pSFS2A plasmid whereby the drug-resistant gene SAT1 can be ''flipped out'' by the FLP gene. 101 Since the Candida albicans MAL2 promoter controlling the FLP gene is not well regulated in S. cerevisiae, we changed the CaMAL2 promoter to the ScGAL1 promoter using an In-Fusion kit (Takara Bio USA Inc., Mountain View, CA, USA) to generate pSFS4A plasmid. Next, the lexA-binding domain was inserted into the pSFS4A plasmid via the ApaI and XhoI sites to construct pSFS4A-lexA.
To insert lexA into the C-terminal domain of TFs, we inserted the 5 0 and 3 0 gene-specific homologous regions via the KpnI-ApaI and NotI-SacI sites on the pSFS4A-lexA plasmid, respectively. For the 5' region, the homologous 300-500 bp region before the translation stop site was amplified using the primer pair TF-kpnI-F1 and TF-ApaI-R1. For the 3' region, the homologous 300-500 bp region after the translation stop site was amplified using the primer pair TF-NotI-F3 and TF-SacI-R3. The resulting fragment was transformed into yeasts by a modified electroporation method. 102 Log-phase cells ( To flip out the drug-resistant marker SAT1, cells were cultured in 2% galactose medium overnight to express the FLP gene and remove SAT1 through recombination at the FLP recombination target sequence. For the promoter activity assay, the promoter was inserted into the pRS41H-lacZ plasmid using an In-Fusion kit (Takara Bio USA Inc.). The pRS41H-lacZ plasmid was constructed by inserting the lacZ fragment amplified from the pSH18-34 plasmid into the pRS41H plasmid. 100,103 The full-length promoters used for the activity assay were amplified using the primer pair promoter-infu-F and promoter-infu-R, and the promoter lacking TF binding sites was constructed by two-fragment PCR. Fragment I was amplified using primers promoter-infu-F and promoter-delTF-R, and fragment II was amplified using primers promoter-delTF-F and promoter-infu-R. The promoter fragment was ligated into the ApaI-linearized pRS41H-lacZ vector using an In-Fusion kit (Takara Bio USA Inc.). The sequence-checked plasmid was transformed into yeasts using the lithium chloride method. 99 For gene deletion, the 5 0 homologous region was amplified using primers gene-pF1 and gene-del-R1, and the 3 0 homologous region was amplified using primers gene-del-F2 and gene-dR2. The URA selection marker was amplified from the plasmid pUG72. 104

Construction of the genetic distance and expression distance trees
For calculating the genetic distance, we first downloaded the raw reads of each strain (LAB and ML strains from Song et al., 105 NA and SK strains from Strope et al., 106 and the WA strain was resequenced by this study) and mapped them to the reference genome (saccharomyces_cerevisiae_R64-2-1_20150113) by the BWA-MEM. 90 Then, we called the variants by GATK (version 3.8) and performed the hard filtering for SNPs and Indels. 107 The filtered variants were merged into the g.vcf format and then transformed into a fasta file by VCF-Kit. 91 The genetics distance was calculated by the R package, ape (function dist.dna in the ''TN93'' model), 92 and the neighbor-joining tree was constructed by BIONJ. 108 We took the average Transcripts Per Kilobase Million (TPM) value from three biological repeats divided by 1000 to calculate the Euclidean distance among strains by dist function in R, and then constructed the neighbor-joining tree by BIONJ. 108

Heatmap and clustering analysis
We used the function heatmap.2 in the R software (gplots package) for plotting the heatmap. To cluster the strains by the input data, we used the Euclidean method in the function dist to calculate the distance and clustered by completing the method in the hclust function in R.

V-src activity assay
Yeast cells were transformed with the Y316V-srcv5 plasmid so the v-Src protein abundance and activity could be measured directly. 89 The v-Src gene was driven by the GAL promoter. Since the GAL promoter was not induced in the WA strain in galactose-containing medium, we removed the WA strain from the To test whether our predicted targets are consistent with the targets identified by experiments, we used the the ChIP-seq data from a recent study. 32 Around 43% of the ChIP targets have the expected TF binding motif in the upstream 500-bp region from the translation start site (Table S10). Some of the ChIP-seq data did not have biological repeats so the data quality was relatively noisy. If we only considered the ChIP-seq data with multiple biological repeats, 58% of the ChIP-seq targets have the expected TF binding motif.

Identification of candidate TFs
To examine if TF targets exhibited differential Hsp90-dependent expression among strains, we performed nonparametric Kruskal-Wallis tests (KW test) and then adjusted the multiple tested p-values using the Bonferroni method. Four TFs-Aca1, Cep3, Cup9, and Yrm1-were excluded from this analysis because they exhibited fewer than three targets in at least one strain. The remaining 187 TFs were tested. Twenty-one of those 187 TFs passed the KW test (multiple testing correction adjusted p-value <0.05). Next, we examined if expression of TF targets differed between any pair of strains by performing two-sided nonparametric Wilcoxon rank-sum tests followed by p-value adjustment using the Bonferroni method. The multiple testing correction adjusted p-value <0.05 was used to define significance. We also tested these candidates using only the TF targets overlapped with the ChIP-seq data (Table S10). 32 With a few exceptions, most of the strain-pairs showed similar statistical results (Table S11).
To further identify high-confidence candidate TFs, we tested if the strain-pair-specific genes, i.e., genes exhibiting significant differences between a given pair of strains, are enriched among TF targets using one-sided two-proportional tests (proportion greater than the background) followed by adjustment of multiple tested p-values using the Bonferroni method. Eleven high-confidence Hsp90-dependent TFs passed the test (multiple testing correction adjusted p-value <0.05). To cluster the TFs based on the Wilcoxon ranksum test and enrichment test results, we labeled the TF by 0 if the adjusted p-value R0.05 in the Kruskal-Wallis test for a strain pair, by 1 if the adjusted p-value <0.05 in the Kruskal-Wallis test and the adjusted p-value R0.05 in the Wilcoxon rank-sum test, and by 2 if the adjusted p-value <0.05 in the Kruskal-Wallis and Wilcoxon rank-sum tests. These three groups represent non-Hsp90-dependent TFs, Hsp90dependent TFs, and high-confident Hsp90-dependent TFs, respectively. The data were then analyzed and plotted by R. The distance matrix is calculated by the Euclidean method and the cluster method is the complete method.

Gene ontology analysis
Biological processes as GO terms were assigned to each gene based on the go_slim_mapping.tab file downloaded from the SGD website (https://downloads.yeastgenome.org/curation/literature/). In total, 5748 genes have GO annotation, gene expression, and genomic sequence data for all analyzed strains. For each GO term, we conducted a hypergeometric test, followed by a multiple test correction using the Bonferroni method, to establish if genes are enriched for specific biological processes (multiple testing correction adjusted p-value <0.05).

Quantitative real-time PCR
We subjected 2 mg total RNA to cDNA synthesis using a HighCapacity cDNA Reverse Transcriptase Kit (Applied Biosystems, Waltham, MA, USA). The cDNA was diluted 16-fold and 5 ml was subjected to realtime quantitative PCR using gene-specific primers (listed in Table S16) and SYBR Green PCR master mix in an ABI-7000 sequence detection system (Applied Biosystems). The samples were measured by three technical repeats.

b-Galactosidase assays
To test the activity of TFs, we performed one-hybrid or promoter activity assays as described previously. 102 The activity of lacZ-encoded b-galactosidase was used to represent TF or promoter activity. For the b-galactosidase assay, log-phase cells cultured in the presence of DMSO or 50 mM GdA were harvested and resuspended in 400 ml Z buffer (0.06 M Na 2 HPO 4 , 0.04 M NaH 2 PO 4 , 0.01 M KCl, 1 mM Mg 2 SO 4 ・7H 2 O, 2.5% b-mercaptoethanol). Next, the cells in 100 ml suspension buffer were lysed by adding 15 ml 0.1% SDS and 30 ml chloroform and then vortexed for 15 seconds. We added 0.2 ml of 4 mg/ml ONPG (o-nitrophenyl-b-galactopyranoside) (N1127, Sigma-Aldrich) substrate into the suspension and vortexed for a further 15 seconds. Since b-galactosidase is activated at 37 C, we incubated the suspension at 37 C until the color of the suspension changed to yellow.