Heterologous Hsp90 promotes phenotypic diversity through network evolution

Biological processes in living cells are often carried out by gene networks in which signals and reactions are integrated through network hubs. Despite their functional importance, it remains unclear to what extent network hubs are evolvable and how alterations impact long-term evolution. We investigated these issues using heat shock protein 90 (Hsp90), a central hub of proteostasis networks. When native Hsp90 in Saccharomyces cerevisiae cells was replaced by the ortholog from hypersaline-tolerant Yarrowia lipolytica that diverged from S. cerevisiae about 270 million years ago, the cells exhibited improved growth in hypersaline environments but compromised growth in others, indicating functional divergence in Hsp90 between the two yeasts. Laboratory evolution shows that evolved Y. lipolytica-HSP90–carrying S. cerevisiae cells exhibit a wider range of phenotypic variation than cells carrying native Hsp90. Identified beneficial mutations are involved in multiple pathways and are often pleiotropic. Our results show that cells adapt to a heterologous Hsp90 by modifying different subnetworks, facilitating the evolution of phenotypic diversity inaccessible to wild-type cells.


Author summary
Biological processes in living cells are often carried out by gene networks. Hubs are highly connected network components important for integrating signal inputs and generating responsive functional outputs. Heat shock protein 90 (Hsp90), a versatile hub in the protein homeostasis network, is a molecular chaperone essential for cell viability in all tested eukaryotic cells. In yeast, about a quarter of the expressed proteins are profoundly influenced when Hsp90 activity is reduced. Despite its pivotal role, we found that the function of Hsp90 has diverged between two yeast species, Yarrowia lipolytica and Saccharomyces cerevisiae, which split about 270 million years ago. To understand the impacts and adaptive strategies in cells with an altered network hub, we conducted laboratory evolution experiments using a S. cerevisiae strain in which native Hsp90 is replaced by its counterpart in Y. lipolytica. We observed different fitness gain or loss under various stress conditions in individual evolved clones, suggesting that cells adapted via different evolutionary paths. Genome sequencing and mutation reconstitution experiments show that beneficial

Introduction
Biological functions are often considered the collective output of genes interacting coordinately in a network [1][2][3][4]. Although most genes in a network only have a few links, some "hub" genes are highly connected to network components and play a vital role in maintaining intra-network connectivity [2]. As a consequence, hub genes usually possess pleiotropic functions and have a profound influence on cell fitness [5,6]. Phenotypic variation was enhanced in yeast cells in which hub genes were disrupted, revealing a role for hub genes in network robustness against genetic or environmental perturbations [7]. Comparative genomics has shown that some hub genes can diverge significantly at the sequence level. However, little is known about to what extent the function of hub genes can change over long-term evolution or how changes in hub genes influence the evolutionary trajectory of an entire network.
The molecular chaperone heat shock protein 90 (Hsp90) represents a well-documented network hub. Hsp90 is present in bacteria, fungi, plants, and animals [8] and is essential for cell viability in all studied eukaryotic species [9][10][11][12][13]. Large-scale studies in S. cerevisiae indicate that Hsp90 is one of the central hubs in both physical and genetic interaction networks [14,15]. More than 10% of the whole yeast and human proteomes have been shown to be clients of Hsp90 [16]. Moreover, the yeast Hsp90-dependent proteome is enriched with other hubs and essential genes [17].
In general, Hsp90 acts at a late stage of protein folding. It assists partially folded clients to reach structural maturation [16], or it stabilizes the components of multi-subunit protein complexes and facilitates their assembly [18][19][20]. Hsp90 has a high degree of specificity for client recognition, which is attributed to its sophisticated co-chaperone system [21,22]. Interestingly, it has been found that the activity or stability of some Hsp90 client proteins can become Hsp90 independent after introducing mutations [23,24], suggesting that the link between Hsp90 and its clients may be changed over time.
Apart from these physiological functions, Hsp90 has been suggested to buffer genetic and nongenetic variations to ensure developmental stability in fungi, plants, and animals [25][26][27][28]. Organisms exhibit constant phenotypes under normal conditions, even if they host genetic or nongenetic variations. However, when Hsp90 activity is compromised or titrated by extreme environmental stresses, heterogeneous phenotypes are revealed in a population. This increased phenotypic diversity is hypothesized to help the population survive or adapt during drastic environmental changes [29].
To understand the implications of hub evolution, we replaced the Hsp90-coding gene of S. cerevisiae with orthologs from other yeast species and examined their functions. S. cerevisiae cells carrying Y. lipolytica-Hsp90 (Ylip-Hsp90) exhibited decreased fitness in most tested conditions but improved growth in high-salt medium, suggesting that the function of Ylip-Hsp90 has differentiated from its S. cerevisiae ortholog. We further evolved the Ylip-Hsp90-hosting strain to study how cells respond to a changed network hub. The evolved clones adopted various evolutionary trajectories, but none involved changing the heterologous Hsp90 per se. Furthermore, fitness measurements in different stress environments showed that these clones had evolved diverged phenotypes that differed from each other and also deviated from wild-type S. cerevisiae cells.

Adaptation to a heterologous Hsp90 in laboratory evolution experiments
Alteration of a network hub often has a strong impact on cell fitness [5,6]. Cells may be selected to quickly restore the original functions of the hub, or the network may be rewired to accommodate or compensate for the altered hub. In the latter scenario, it may also open the possibility for cells to evolve novel phenotypes. To explore the possible evolutionary trajectories of an altered Hsp90 network, we evolved 12 independent haploid lines carrying Ylip-HSP90 and 12 control lines carrying Scer-HSC82 that were all derived from the same parental S. cerevisiae clone. These cell lines were grown in regular rich medium (Yeast extract-Peptone-Dextrose [YPD]) at 28˚C so that their major selective pressure is to restore the general growth defects caused by Ylip-HSP90. In each evolving culture, two isogenic strains carrying either the green or red fluorescent protein marker (Fig 2A) were mixed at an approximate initial ratio of 1:1. By monitoring the proportions of the two fluorescence-labeled populations at different time points, we could estimate when the first beneficial mutation was fixed in each culture (S3 Scer-HSC82 control lines were set up to evolve in rich media at 28˚C. Two isogenic subpopulations (red: DsRed-labeled subpopulation; green: GFP-labeled subpopulation) were mixed in an approximate ratio of 1:1 in the initial cultures. The subpopulation that obtained a beneficial mutation (noted by a star) would start expanding its frequency in the evolved culture. (B) The Ylip-HSP90 lines become fixed more quickly than the Scer-HSC82 lines. The mutation fixation of each culture was determined as the time when the frequency of a subpopulation was observed to be greater than 95% according to flow cytometry. The control line con6 and the evolving lines, evo1, evo2, and evo12 were excluded because of culture contamination during evolution. The frequency dynamics of each evolved line are shown in S3 Fig. ( Ylip-HSP90 hosting cells were initially under high selective pressure to resolve the incompatibility between the heterologous Ylip-Hsp90 and its interactive network in the S. cerevisiae background. In contrast, the Scer-HSC82 lines only experienced very low or no selective pressure, because the evolution experiments were conducted under a normal growth condition.  After evolving for around 2,200 generations, three Ylip-HSP90 (evo1, evo2, and evo12) lines and one Scer-HSC82 (con6) line were lost because of contamination. We further analyzed the remaining evolved cultures. The fitness of evolved Ylip-HSP90 populations in YPD at 28˚C improved by 16%-24% compared with their ancestor (Fig 2C). The degree of improvement was significantly higher than that observed in the evolved Scer-HSC82 populations, suggesting that the adaptive phenotype is likely to be Ylip-HSP90 specific. When the fitness of evolved Scer-HSC82 clones was compared with their ancestor, no significant difference was observed in four evolved lines (S2 Table), suggesting that both genetic drift and beneficial mutations contributed to the evolved phenotypes of Scer-HSC82 populations.
To characterize the evolved phenotype, we isolated three single colonies from each evolved culture and measured their fitness. The colonies with the highest fitness from each population were chosen to represent the lineage for further examination (denominated E3 to E11 for evolved Ylip-HSP90 clones and C1 to C12 for evolved Scer-HSC82 clones). It was previously reported that S. cerevisiae cells often converge to a diploid state during laboratory evolution [36]. We examined the genome content of evolved clones using flow cytometry and found that all of the evolved Scer-HSC82 clones had become diploid or triploid (S4A Fig). In contrast, most of the evolved Ylip-HSP90 clones remained haploid, except for E10, which had become diploid. To ensure that the ploidy expansion observed in individual clones is representative of the whole population, we also examined the evolved cultures using flow cytometry and observed the same results. Because all the cultures were evolved under the same selective regime, the difference in ploidy observed in evolved cultures further implicates functional divergence between Ylip-Hsp90 and Scer-Hsc82. To further test the effect of Hsp90 on ploidy expansion or maintenance, we generated diploid strains by mating haploid cells and compared their fitness with haploid strains. Haploid Ylip-HSP90-carrying cells exhibited higher fitness than diploid Ylip-HSP90-carrying cells (S4B Fig). This result provides the possible reason why most Ylip-HSP90 lines remained haploid after evolution.

Fitness improvement in the evolved clones is specific to Ylip-Hsp90
Because the evolved clones restored the fitness of Ylip-HSP90 cells in rich media, we wondered whether they would become sensitive to hypersaline stress like Scer-HSC82 cells. When evolved Ylip-HSP90 clones were tested on 0.2 M LiCl-containing plates, they displayed similar or higher fitness compared with ancestral Ylip-HSP90 strains (Fig 3A), indicating that the Y. lipolytica-specific hypersaline-resistant phenotype had been stably maintained in the Hsp90 network of the evolved cells.
The function of Ylip-Hsp90 might be enhanced if its expression level is increased or the protein is mutated in the evolved clones. We measured protein levels of Ylip-HSP90 using western blots. All evolved clones presented similar levels of the Hsp90 protein compared to the ancestral Ylip-HSP90 strains (S4C Fig). We also sequenced the Ylip-HSP90 gene of all evolved clones and found no mutation in any of them. Together, these data indicate that the evolved cells had not improved their fitness by modifying the Ylip-HSP90 gene.
The observed fitness improvement in evolved Ylip-HSP90 clones might be specific to the heterologous Hsp90 or simply be due to mutations that enhance general cell proliferation (e.g., loss of the costly genes or gain of glucose uptake genes). To differentiate between these two possibilities, we replaced the Ylip-HSP90 gene of the evolved and ancestral clones with Scer-HSC82 and measured their fitness in YPD at 28˚C. The extent of fitness improvements was drastically diminished when Ylip-Hsp90 was replaced ( Fig 3B and S2 Table), indicating that the beneficial effects of evolved mutations are dependent on Ylip-Hsp90.

Ylip-HSP90 clones
Because Hsp90 is an important molecular chaperone involved in maintaining protein homeostasis, we suspected that the reduced fitness of ancestral Ylip-HSP90 cells in normal conditions might reflect disturbed internal proteostasis caused by the compromised Hsp90 network. We used an Hsp104 aggregation assay to investigate this possibility [37]. Hsp104 is a disaggregase that localizes to protein aggregates upon heat stress and helps with the disaggregation and refolding process. Ancestral and evolved cells containing fluorescent Hsp104-BFP fusion protein were first grown at 25˚C and then shifted to 37˚C for high-temperature adaptation. By monitoring the fraction of cells containing Hsp104-BFP foci during heat treatment, we could assess how well protein homeostasis is maintained.
Ancestral Ylip-HSP90 cells contained more Hsp104-BFP foci than Scer-HSC82 cells before and during heat treatment, suggesting that Ylip-HSP90 cells indeed have a high endogenous and then shifted to 37˚C to induce heat adaptation. The fraction of cells containing Hsp104-BFP foci was counted at 0, 1, and 3 hours after the temperature shift. At least seven replicates were measured for each strain, and more than 75 cells were counted for each replicate. Asterisks indicate that the fraction of cells containing Hsp104-BFP foci is significantly lower compared with ancestral Ylip-HSP90 cells from the same time point (p-value < 0.01 after the Benjamini-Hochberg correction, N � 7, one-tailed Wilcoxon rank-sum test). (D) Evolved Ylip-HSP90 clones E3, E4, and E9 recover from elongated cell morphology, whereas E5 and E11 become even more elongated than the ancestor. The ratio of the long-versus short-axis diameter of each cell was measured, and more than 80 cells were examined for each strain. The distribution of each sample was compared to ancestral Ylip-HSP90 cells and asterisks indicate a p-value <0.01 after the Benjamini-Hochberg correction by two-tailed Wilcoxon rank-sum test. The numerical data used in panels (B, C, D) are included in S1 Data. See also S5 Fig Our results suggest that different evolved lines probably took different paths to fix the compromised Hsp90 network. In addition to attenuated protein homeostasis, ancestral Ylip-HSP90 cells also exhibited an elongated cell morphology that has previously been observed in cells lacking sufficient Hsp90 activity to stabilize its client proteins [28]. When we examined the cell morphology of evolved Ylip-HSP90 clones, we found that only E3, E4, and E9 cells recovered from the elongated cell morphology. E5 and E11 evolved to an even more elongated form, whereas the remaining clones were not significantly different than their ancestor ( Fig 3D and S5C Fig). These results reinforce that different adaptive paths had been taken by the individual evolved lines.

Independent Ylip-HSP90 clones evolved through different adaptive trajectories
Both the Hsp104 aggregation assay and our cell morphology measurements suggested that the phenotypes of evolved Ylip-HSP90 lines are diverged from each other and from their ancestor and Scer-HSC82 cells. To obtain a more complete picture of the evolved phenotypes, we subjected the evolved clones to two growth media and nine different stress conditions that challenge different aspects of protein homeostasis and then measured their fitness (see Materials and methods). Our results revealed a broad spectrum of adaptive phenotypes in the evolved clones ( Fig 4A and S2 Table) so that, even if some strains shared similar fitness levels for one or two conditions, their fitness levels differed under other conditions. Moreover, for many of the conditions we tested, the evolved clones exhibited higher fitness than Scer-HSC82 cells (S2 Table and Fig 4A). We also tested whether the observed fitness was correlated with the protein abundance of Hsp90 in evolved Ylip-HSP90 clones. No significant correlation was observed in any tested condition (S3 Table).
To obtain a more general overview of the evolved features in different clones, we performed hierarchical clustering and principal component analysis (PCA) according to their fitness measurements under the 11 tested conditions ( Table). Hierarchical clustering revealed that clones E4 and E9 were more similar to Scer-HSC82 cells than other evolved clones. The first principal component (explaining 33.8% of the variance) of our PCA showed that most clones evolved toward the Scer-HSC82 phenotype (S6A Fig). Nonetheless, they were not closely clustered. Phenotypic divergence between the evolved clones was further revealed by the second and third principal components (explaining 28.4% and 14.2% of the variance, respectively). These results suggest that Ylip-HSP90 lines have evolved diverged phenotypes that both differ from each other and from ancestor and Scer-HSC82 cells.
The observed phenotypic variation in the evolved Ylip-HSP90 clones suggests that after the heterologous Hsp90 was introduced, different Hsp90-related functions in individual evolving lines were modified to regain fitness. Under this hypothesis, high phenotypic variation should not occur in the evolved Scer-HSC82 clones (i.e., the control lines), because the network hub has not changed. We performed two types of analysis to address this issue.
We first measured the fitness of evolved Scer-HSC82 clones under the same 11 growth conditions (S6B Fig) and calculated the variance of fitness improvements among evolved Ylip-HSP90 and Scer-HSC82 clones in individual conditions. We then compared the variance between evolved Ylip-HSP90 and Scer-HSC82 groups. The data show that the evolved Ylip-HSP90 clones exhibited greater variance than the evolved Scer-HSC82 group (p-value = 0.039, one-tailed Student t test), indicating that responses of evolved Ylip-HSP90 clones to different stress conditions were more varied than evolved Scer-HSC82 clones (S6C Fig).
In the second analysis, we calculated the Pearson correlation distances between fitness levels under different conditions (see Materials and methods, S5 Table) to quantify the variation in   clones under 11 different growth conditions, many of which challenge different aspects of protein homeostasis. Cells were grown in liquid cultures and growth rates were measured by plate readers. Error bars are standard errors, N � 3. Unless specified otherwise, cells were grown at 28˚C. Fitness improvement was calculated by comparing the fitness of the evolved (or Scer-HSC82) clone with that of the ancestral Ylip-HSP90 strain. (B) Hierarchical clustering of the fitness values (i.e., growth rates) reveals divergent evolution of the Ylip-HSP90 clones. For many of the conditions, some evolved Ylip-HSP90 clones display enhanced fitness compared with Scer-HSC82 clones. Individual fitness values were divided by the mean of all fitness values in the same condition to normalize the fitness data to a similar scale between conditions. The color code represents the normalized fitness value (red: high fitness; blue: low fitness). (C) PCA of the fitness values shows that all Ylip-HSP90 clones evolved diverged phenotypes, scattered across the three principal component dimensions. Explanatory power is shown in brackets next to each principal component. (D) The evolved phenotypes of Ylip-HSP90 clones are more diverse than those of Scer-HSC82 clones. Pearson correlation distance (d = 1 − Pearson correlation coefficient) was used to measure whether the evolved clones displayed similar fitness trends under different conditions. A smaller distance indicates a higher similarity between the two conditions. In general, phenotypic distances for evolved Ylip-HSP90 clones were significantly greater than those of evolved Scer-HSC82 clones (p = 0.001, one-tailed Wilcoxon rank-sum test), suggesting that evolved Ylip-HSP90 clones have more divergent phenotypes. The color code represents the Pearson correlation distance value (cyan: high similarity; red: low similarity). The numerical data used in the figure are included in S1 Data. See also S6 Fig. AZC distribution patterns of individual fitness between each pair of conditions. If individual clones show similar ranking patterns for fitness between two conditions, Pearson correlation distances will be small. In contrast, large Pearson correlation distances indicate that individual evolved clones behave very differently in two conditions. Overall, pairwise comparisons revealed that the evolved Scer-HSC82 clones behaved similarly under different stress conditions, except for the complete supplement mixture (CSM) and azetidine-2-carboxylate (AZC) media ( Fig 4D and S5 Table). The stress condition-based phenotypic distance of evolved Scer-HSC82 clones is significantly smaller than that of evolved Ylip-HSP90 clones (S6D Fig), supporting our hypothesis that high phenotypic variation only occurred in evolved Ylip-HSP90 clones.

Most mutations in evolved Ylip-HSP90 clones are associated with various Hsp90-related functions
To explore the genetic basis of the evolved phenotypes, we sequenced the genomes of the evolved Ylip-HSP90 clones. On average, each evolved clone hosted eight protein-sequencechanging mutations (seven nonsynonymous and one insertion/deletion) and 14 total mutations (S7 Fig). The calculated mutation rate (5.45 × 10 −10 per bp per generation) is comparable with the estimated spontaneous mutation rates from previous studies [38,39]. Only a small group of genes were recurrently mutated in multiple evolved clones (i.e., BRE5 in three clones and BUD2, FMP30, SIR3, and HXK2 in two clones, S6 Table). We investigated the BRE5 gene further in the reconstitution experiments. In evolved clones, we found that large-scale copy number variation only occurred in mitochondrial genomes and ribosomal DNA, and that small-scale copy number variation was also rare outside subtelomeric regions (S7 Table).
When the gene functions of all identified protein-sequence-changing mutations were analyzed further (Table 1 and S8 Table), we found that most of them could be categorized into several Hsp90-related functions [14,40], including protein synthesis, folding and modification, intracellular trafficking, ribosome biogenesis, cellular conjugation, chromatin organization, and signaling, as well as mitochondrial physiology (Fig 5). Many of these genes have known interactions with Hsp90 or belong to the Hsp90-dependent proteome [17,41], raising the possibility that evolved yeasts might adapt by modifying different Hsp90 subnetworks.

Segregant analysis identifies beneficial mutations with strong fitness effects
Not all the mutations in the evolved clones are likely to contribute equally to the adaptive phenotypes. We performed segregant analysis to map the mutations with strong effects in three evolved clones (E4, E5 and E7). These clones were selected because E4 and E5 clones showed opposite trends in cell morphology, and E7 has the largest improvement in protein homeostasis (Fig 3). The evolved Ylip-HSP90 clones were backcrossed to the ancestral Ylip-HSP90 strain and their F1 spores (from 40-45 four-viable spore tetrads) were analyzed. In all three evolved clones, our genetic analysis indicated that multiple loci were involved in the evolved phenotypes (S9 Table). The F1 progeny with phenotypes similar to evolved or ancestral clones were grouped into good and bad spore pools, respectively (see Materials and methods). Genomic DNA was isolated from both pools and then subjected to Illumina sequencing.
Based on our computational simulation, two criteria were used to select the mutations with strong effects: (1) the mutation frequency in the good spore pool is greater than 70% or (2) the difference in mutation frequencies between good and bad spore pools is greater than 45%. Nine mutations-sec9-E310K and tap42-F72S from E4; izh3-S115 � ( � indicating a premature stop codon), bud4-P981fs (fs indicating a frameshift mutation), bre5-L433F, and mdm12-F148S from E5; and siz1-P401L, pbp1-N318fs, and nop53-E289_N290insE (ins indicating an insertion) from E7 (S10 Table)-were identified and used for reconstitution experiments to confirm their effects.

Individual mutations exert condition-specific effects
Our phenotypic assays had revealed that independent Ylip-HSP90 lines evolved diverse phenotypes (Fig 4). To understand how the evolved mutations contribute to the broad spectrum of phenotypes, we used the clustered regularly interspaced short palindromic repeats/CRISPRassociated protein 9 (CRISPR/Cas9) system to reintroduce these mutations into the ancestral Ylip-HSP90 strain and examined the fitness of these mutants under different stress conditions (Fig 6A, S8 Fig and S11 Table). Four major conclusions can be drawn from the phenotypic assays. First, individual mutations usually had dominant effects under at least one stress condition (e.g., in clone E7, the siz1, pbp1, and nop53 mutants exhibited strongest effects under conditions of 37˚C, macbecin II treatment, and doxycycline treatment, respectively), confirming that these mutations really contribute to the evolved phenotypes. Second, several mutations exhibited positive effects under one stress condition but negative effects under others (e.g., in clone E7, pbp1 mutants had a positive effect under macbecin II treatment but a negative effect under the MG132 stress condition), indicating that antagonistic pleiotropy is common among the evolved mutations. Third, in a few stress conditions (e.g., for clone E5 in medium containing macbecin II), the sum of individual effects exceeded the effect observed in the evolved clone, suggesting the existence of negative epistasis between mutations. Fourth, for some other conditions (e.g., clone E4 in medium containing MG132), the sum of all individual effects was still far below the effect observed in evolved clones. This outcome can be explained if there is positive epistasis between mutations. Alternatively, some mutations with weak individual effects but that are strongly epistatic to other mutations might have been missed by our mapping exercise. We selected one mutation from the E4 clone, cbk1-L461F, that exhibited a mutation frequency slightly lower than our segregant analysis criteria to test these possibilities. The  cbk1-L461F single mutation did not result in any significant effects under any tested stress condition ( Fig 6A). However, when cbk1-L461F was combined with sec9-E310K or tap42-F72S and tested in all conditions, it displayed positive epistatic effects with sec9-E310K and tap42-F72S in medium containing brefeldin A (BFA) and MG132, respectively (Fig 6B and  S11 Table).
In addition to its fitness improvement, our Hsp104 aggregation assay showed that the E7 clone has a strong aggregate clearance phenotype (Fig 3C). This allowed us to characterize the contribution of individual mutations to protein homeostasis. We found that all three single mutations (siz1, pbp1, and nop53) reduced the formation of Hsp104 foci at 37˚C, but pbp1-N318fs had the strongest effect (Fig 6C). When all three mutations were combined in a single clone, the further enhanced aggregate clearance phenotype was observed at the 3-hour time point. These results reveal a specific role of pbp1-N318fs in maintaining protein homeostasis under heat stress conditions.

Discussion
The Hsp90 network performs many crucial functions in eukaryotic cells and probably represents a robust system resilient to genetic perturbations [16,40]. Nevertheless, the ability of Ylip-Hsp90 to enhance hypersaline tolerance in S. cerevisiae cells and the compromised fitness of such cells under other stress conditions indicate that even a central hub like Hsp90 can change its functions or partners during evolution. This seeming discrepancy may be explained by a computational simulation study of complex signaling circuits that showed phenotypic robustness and variability are positively correlated [42]. Because a broad range of genotypes usually supports a robust phenotype, a greater diversity of phenotypes is accessible from among the genotypic neighborhood. This model implies that a robust system like the Hsp90 network can still be genetically adjusted to facilitate adaptation to novel environments. If the Hsp90 network is continually modified in a specific direction (e.g., to tolerate hypersaline stress), it may eventually lead to alterations in a network hub such as Ylip-Hsp90. Evolvability represents an organism's capacity to generate selectable phenotypic traits during evolution [43], so the level of evolvability may influence the likelihood of an organism surviving in or adapting to a new habitat. In developmental biology, large phenotypic changes or innovations are often associated with changes in the transcriptional network hubs [44][45][46]. Despite strong evidence supporting a link between network hubs and evolvability, it is difficult to study directly the effect of changed network hubs on evolutionary trajectories using natural populations. Our laboratory evolution experiments provide an alternative approach to address this issue. In the current study, we replaced the Hsp90-coding gene of S. cerevisiae with the ortholog from a distant species. Nonetheless, a similar evolutionary process may occur in nature under the following conditions: (1) a hybrid progeny carries incompatible gene pairs that compromise the function of a network hub. Genetic incompatibility has been observed to exist between closely related species and cause defects in different important cellular functions [47]. Although these types of hybrids have lower fitness and are likely to be eliminated in most cases, they may also create the chance for the organism to evolve new phenotypes, especially in the hybrid zone where hybrids are constantly generated. (2) Mutations occur in critical domains of a hub gene. Mutagenesis experiments in Hsp90 have shown that even a single amino acid change can alter the activity or interacting partners of the protein completely [48]. Such mutations will often result in the low fitness of the cell, as in the hybrid case. However, they may lead to the evolution of new phenotypes if compensatory mutations that restore fitness occur and are subsequently fixed by selection [49]. (3) A multifunctional gene is obtained through horizontal gene transfer (HGT). In recent years, HGT has been shown to occur frequently, even in eukaryotes [50,51]. A foreign multifunctional protein may disturb the original cellular network and open a new evolutionary path.
Network structures contribute to phenotypic robustness against genetic or environmental perturbations [52]. When network hubs are compromised, pre-existing cryptic genetic or nongenetic variation may be exposed that facilitates short-term population adaptation [7]. Our study reveals a different role played by network hubs in long-term evolution in a population without much genetic variation. S. cerevisiae cells immediately gained hypersaline tolerance like Y. lipolytica cells solely by replacing a single hub gene, HSP90, with the corresponding ortholog. This is similar to a previous observation that bacteria dramatically increased the liquid cultures and growth rates were measured by plate readers. Individual mutations exerted dominant effects in at least one condition, and antagonistic pleiotropy under different stress conditions was widely observed. Error bars are standard errors, N � 3. Asterisks indicate that the fitness values are significantly different from the ancestral Ylip-HSP90 strain (p < 0.05 after the Benjamini-Hochberg correction, two-tailed Wilcoxon ranksum test). See also S8 Fig. (B) Positive genetic interactions between cbk1-L461F and the other two mutations result in the fitness improvement not observed in single mutants when grown in media containing BFA or MG132. Error bars are standard errors, N � 4. Asterisks indicate that the fitness of double mutants is significantly greater than that of single mutants (p < 0.05 after the Benjamini-Hochberg correction, two-tailed Wilcoxon rank-sum test). (C) All three single mutations in the E7 clone reduce the formation of Hsp104 foci at 37˚C, but pbp1-N318fs has the strongest effect. Each measurement was compared with that of ancestral Ylip-HSP90 cells from the same time point, and asterisks indicate a p-value <0.01 after the Benjamini-Hochberg correction, according to one-tailed Wilcoxon rank-sum test. N � 14. Error bars are standard errors. The numerical data used in the figure are included in S1 Data. AZC, azetidine-2-carboxylate; BFA, brefeldin A; CRISPR/Cas9, clustered regularly interspaced short palindromic repeats/CRISPR-associated protein 9; CSM, complete supplement mixture; CYH, cycloheximide; Dox, doxycycline; EtOH, ethanol; GuHCl, guanidine hydrochloride; L461F, leucine 461 to phenylalanine; MacII, macbecin II; N318fs, asparagine 318 frameshift; Ylip, Yarrowia lipolytica; YPD, Yeast extract-Peptone-Dextrose. growth rate at lower temperatures with the introduction of psychrophilic chaperonins [53]. Moreover, the novel Hsp90 allowed individual clonal populations to evolve diverse phenotypes unobserved in cells hosting native Scer-Hsc82. These novel phenotypes may enable evolved cells to explore niches inaccessible to Scer-Hsc82-hosting cells and, in the long term, may broaden opportunities for population diversification and speciation.
Previous laboratory evolution experiments in bacteria showed that when essential ribosomal genes or translation factors were replaced by their orthologs from other bacteria or the ancestral sequence, expression levels of the introduced orthologous genes were often amplified to compensate for the incompatibility [54][55][56]. In our evolved clones, protein levels of Ylip-Hsp90 were not changed, suggesting that the cells were able to quickly gain other mutations to ameliorate the growth defect in the conditions under which they evolved. One major function of network hubs is to coordinate different subnetworks (or pathways) so that the cell can quickly respond to environmental changes. However, when a hub is compromised or altered (such as in the Ylip-Hsp90-hosting cells), any mutation that can improve cell fitness under the existing growth conditions will be selected, even if it has costs in other conditions or distorts the balance between different subnetworks. The evolved Ylip-HSP90 clones restored only some of the functions of Hsp90 despite that all clones had similar or higher fitness in rich medium compared to Scer-HSC82 cells (Figs 3C, 3D and 4).
By sequencing the evolved clones and reconstituting mutations, we confirmed that individual subnetworks of Hsp90 had been modified in different evolved lineages. As expected, some of the evolved mutations showed strong antagonistic pleiotropy, exhibiting improved cell fitness in rich medium but growth defects under other stress conditions. Because cells are constantly challenged by environmental fluctuations, it is likely that another layer of compensatory mutations can be selected to resolve the fitness conflict of pleiotropic mutations under different conditions [57]. These compensatory mutations may not exhibit adaptive effects by themselves but will be epistatic to the previous pleiotropic mutations. This hypothesis is indirectly supported by our reconstitution experiments, which showed that intergenic epistasis is common among the evolved mutations. After obtaining compensatory mutations, modified subnetworks are likely to stay in the lineage and to have a long-term influence on evolutionary paths.
The fact that no mutations occurred in the Ylip-HSP90 gene indicates a strong constraint on the hub, probably due to pleiotropic effects. Therefore, evolved cells resolved the fitness conflict through modifying the functions or structures of some subnetworks. Similar processes have been observed in the evolution of transcription networks. Although different species may share the same transcription factors, the downstream targets of the orthologs could have been rewired greatly [58,59]. Previous studies of digital organisms have suggested that deleterious mutations might act as stepping stones for the evolution of complex phenotypes inaccessible to the wild type [60,61]. Our study provides direct evidence that diverse phenotypes can evolve when an essential network hub is altered. Robust networks are often connected to a larger phenotypic space [42]. An initial adaptation to a changed or compromised network hub allows a cell to reshape its network structures in different directions, which further broadens its evolutionary potential (Fig 7).

Construction of the HSP90 replacement strains
All the S. cerevisiae strains were derived from R1158 (JYL1821) containing the tTA transactivator with the CMV promoter that allows us to down-regulate the expression level of the TetO 7 promoter-containing gene by adding doxycycline [62]. In S. cerevisiae, the essential function of Hsp90 is performed by two highly similar and functionally redundant genes, HSC82 and HSP82 (protein identity and similarity are 0.969 and 0.987, respectively) [9]. We first generated the HSP90 plasmid shuffling strain (JYL1827) in which both HSC82 and HSP82 were deleted, and the essential Hsp90 functions were performed by an HSC82-containing plasmid (pRS416-Scer-HSC82). The coding regions of HSC82 and its orthologs were PCR amplified from the genomic DNA or cDNA of S. cerevisiae, N. castellii, K. lactis, or Y. lipolytica and cloned into a plasmid (pRS413-pTetO 7 ) containing the kan R -TetO 7 -TATA CYC1 cassette from RP188 [63]. These ortholog-containing plasmids were then transformed into the shuffling strain to replace pRS416-Scer-HSC82. The ortholog-replacement strains were checked using PCR to confirm the loss of pRS416-Scer-HSC82 and are denominated as JYL5001 (Scer-HSC82), JYL1831 (Ncas-HSP90), JYL1833 (Klac-HSP90), and JYL1835 (Ylip-HSP90). In our experiments, the orthologs of Hsp90 were all carried by the plasmid so they could be easily shuffled to check whether the observed effects are ortholog specific.

Experimental evolution
For the evolution experiment, JYL5001 and JYL1835 were transformed with HindIII-linearized GFP (pGS62) and DsRED (pGS64) integration plasmid DNA [64] to generate green and red fluorescent protein-labeled ancestral strains. We mixed the two colored subpopulations of Ylip-HSP90 or Scer-HSC82 cells in an approximate ratio of 1:1 to establish starting cultures and then divided them into 12 evolving and 12 control lines. These cell lines were grown in YPD medium at 28˚C through a daily 1,000-fold dilution (about 10 generations) with an effective population size of 4.69 × 10 5 cells on the first 140 days and 2.79 × 10 6 cells after day 140. The effective population size was estimated using the formula, N e = N 0 × g, in which N 0 is the initial population size and g is the number of generations during a single cycle of growth [65]. We tracked the proportion of the green and red subpopulations in each lineage using the BD LSR II system (Becton Dickinson, Franklin Lakes, NJ). Fixation of the evolved culture is defined when either the green or red subpopulation reaches a population frequency of 95%.

Cell fitness assays
For the fitness spot assay on plates, around 10 3 cells were taken from overnight cell cultures growing in YPD at 28˚C and spotted onto the plates following 10-fold serial dilutions. Plates were then incubated at specified temperatures until individual colonies became visible. Colony size and cell density were used to estimate the cell fitness and survival rate for each condition. For the fitness assay in liquid, overnight cell cultures growing in YPD at 28˚C were inoculated into 120 μL fresh medium in 96-well tissue culture plates at an initial concentration of 0.1 OD 595 , and cell densities (OD 595 ) were measured every 10 minutes using Infinite 200 series plate readers (Tecan, Mannedorf, Switzerland) with a shaking cycle (1 minute shaking and 3 minutes standing) until the culture reached diauxic shift. Growth rates (OD 595 /hour) of the sample were calculated using the 10-point sliding window, and the maximal growth rate was used to represent the fitness value. Fitness improvement compared with ancestors was calculated as (F x /F a −1) × 100%, where F x is the fitness value of evolved or reconstituted clones and F a is the fitness value of the ancestral strain.
To characterize the evolved phenotypes, we measured the fitness of evolved clones under nine different stress conditions: AZC (cat. no. A0760, Sigma-Aldrich), a proline analog that can be incorporated into newly synthesized proteins and cause failure in folding or structural instability [66,67]; cycloheximide (CYH, cat. no. C7698, Sigma-Aldrich), which inhibits the 80S ribosome from translating mRNA into proteins [68]; BFA (cat. no. B6542, Sigma-Aldrich), which blocks protein trafficking and affects the central vacuolar secretory pathway, resulting in ER/Golgi membrane mixing [69] and an unfolded protein response [70]; high temperature (37˚C) and ethanol (EtOH), which induce protein misfolding and change the membrane protein composition [71]; doxycycline, which reduces the expression level of Hsp90 orthologs; GuHCl, which inhibits the disaggregase activity of Hsp104 and reduces cell thermotolerance [72]; MG132 (cat. no. C2211, Sigma-Aldrich), a compound that inhibits protein degradation by blocking proteasome activity [73]; and macbecin II (MacII, cat. no. NSC 330500, NIH, Bethesda, MD), which is a compound that inhibits the ATPase activity of Hsp90. CSM at 28˚C was a suboptimal condition. Apart from CSM, as well as MG132 that was added to synthetic medium (0.17% yeast nitrogen base without ammonium sulfate supplemented with 0.1% proline, amino acids, 2% glucose, and 0.003% SDS) [74], all these chemicals were added into YPD medium. Liquid fitness assays were used in all measurements.

Western blotting for measuring the abundance of Hsp90
We constructed pRS413-pTetO 7 -TAP-Scer-HSC82 and pRS413-pTetO 7 -TAP-Ylip-HSP90 plasmids, in which a tandem affinity purification (TAP) tag was fused to the N-terminus of the Hsp90 protein. These plasmids were then transformed into corresponding ancestral or evolved clones. Cell pellets were collected by centrifugation, washed once with ddH 2 O, and lysed using the NaOH-lysis method [75]. The anti-TAP-tag antibody (cat. no. CAB1001, Thermo Fisher Scientific, Waltham, MA) was used to detect TAP-Scer-Hsc82 or TAP-Ylip-Hsp90 proteins, and the anti-G6PDH antibody (cat. no. A9521, Sigma-Aldrich) was used to detect the G6PDH protein as the internal control.

V-src kinase assay
Ancestral Scer-HSC82 and Ylip-HSP90 strains were first transformed with the plasmid Y316v-srcv5 containing v-src with a C-terminal V5 tag under a galactose-regulated promoter [76]. Cells were then maintained on 2% glucose CSM-Ura plates to prevent the expression of v-src. Before the experiment, cells were grown to log phase in 2% glucose CSM-Ura liquid cultures, transferred to 2% raffinose CSM-Ura, and then grown to about 5 × 10 6 cells/mL. Cells were pelleted, resuspended in 2% galactose CSM-Ura, and grown for 6 hours to induce the expression of v-src. The whole experiment was carried out at 28˚C. Cell pellets were collected by centrifugation, washed once with water, and lysed with the NaOH-lysis method. Total protein extracts were mixed with the sample buffer (50 mM Tris, pH 6.8, 2% SDS, 10% glycerol, 0.004% bromophenol blue, and 2% 2-Mercaptoethanol) and boiled at 100˚C for 5 minutes. The anti-phospho-tyrosine antibody 4G10 (cat. no. 05-321, Merck Millipore, Burlington, MA) was used to detect phosphorylated tyrosine residues, and the anti-V5-Tag antibody (cat. no. MCA1360, BioRad, Hercules, CA) was used to detect the v-src protein.

Flow cytometry of DNA content
Flow cytometry was used to determine the DNA content of ancestral and evolved cells. Total 5 × 10 6 log-phase cells were harvested, resuspended in ice-cold fixation buffer (40% ethanol, 0.1 M sorbitol, and 5 mM EDTA), and then kept at −20˚C overnight. The fixed cells were then washed twice with 1 mL ddH 2 O and then washed with 1 mL PBS + 0.5% Triton X-100. After washing, the cells were treated with 0.5 mL 50 mM Tris-Cl (pH 8.0) with 150 μg/mL RNase A and incubated at 37˚C for 16-18 hours. SYTOX Green nucleic acid stain (Invitrogen, Carlsbad, CA) and 38 mM sodium citrate were mixed in a 1:800 ratio, and 300 μL of the mixture was added into the cell solution. The stained cells were diluted into 1 mL 0.1 M EDTA and sonicated for at least 3 minutes prior to flow cytometry. A total of 10,000 cells were scored for DNA content using the BD FACScan system (Becton Dickinson). The E10 clone was excluded from further examination because of its altered ploidy.

Hsp104 aggregation assay
A blue fluorescence protein (BFP) cassette was amplified from the pFA6a-link-yomTagBFP2-Kan plasmid and integrated into the genome to generate an Hsp104-BFP fusion protein in all tested strains [77]. Log-phase cells growing in YPD at 25˚C were transferred to a 37˚C water bath and cells were collected at 0, 1, 2, and 3 hours after the temperature shift. The samples were resuspended in 200 μL PBS and transferred to a Glass Bottom ViewPlate-96F (PerkinElmer, Waltham, MA) that was precoated with concanavalin A (cat. no. C2010, Sigma-Aldrich). Images of Hsp104-BFP foci were acquired using an ImageXpress Micro XL system (Molecular Devices, San Jose, CA) and analyzed using a custom-built module under MetaXpress High-Content Image Acquisition and Analysis Software (Molecular Devices).

Hierarchical clustering and PCA
For hierarchical clustering, we divided individual fitness values by the mean of all fitness values for the same stress condition to normalize fitness data to a similar scale between conditions. We applied a hierarchical clustering method in R (v3.4.0) [79] to group the evolved Ylip-Hsp90-hosting clones according to their fitness similarity between different conditions. The strain and condition dendrograms were constructed with the hclust function and the heatmap was drawn with the heatmap.2 function in the gplots R package [80]. For PCA, fitness values were used directly, i.e., without normalization, and the prcomp function in R was used to conduct the analysis.

Calculation of Pearson correlation distances
We used the growth rates of all evolved Scer-HSC82 and Ylip-HSP90 clones to conduct pairwise Pearson correlations between each condition. The Pearson correlation distance (d = 1 − Pearson's correlation coefficient) is used to represent the phenotypic distance between each pair of conditions. A large d means the fitness distributions of evolved clones have very different patterns in two conditions. We used the fviz_dist function of the factoextra package in R to draw the heatmap [81].

Whole genome sequencing analysis
Yeast genomic DNA was isolated using the previously described phenol/chloroform method with one more run of 200 μL of PCIA resuspension and 95% ethanol precipitation to remove the residual RNase A [82]. Isolated genomic DNA was fragmented using a Covaris microTUBE M220 Focused-ultrasonicator (Covaris, Worburn, MA) and DNA libraries were prepared using a Tru-Seq DNA PCR-Free HT Library Prep Kit (Illumina, San Diego, CA). Ancestral and each evolved Ylip-HSP90 clone was sequenced to an average depth of 100-fold coverage with a rapid-mode Illumina HiSeq 2500 system (250 bp paired-end reads) and the spore pools were sequenced to an average depth of 200-fold coverage with the NextSeq 500 system (150 bp paired-end reads).
For single nucleotide polymorphism (SNP) analysis, raw reads were mapped, and SNPs with coverage greater than 20-fold and a frequency higher than 35 (for the segregant analysis, the allele frequency threshold was lowered to 10 to include more low-frequency evolved mutations) were called using CLC bio software (CLC bio, Aarhus, Denmark). Evolved mutations were manually identified by comparing the mutation lists between evolved and ancestral Ylip-HSP90 clones, and poor-quality SNPs were excluded after manual checking using the Integrative Genome Viewer [83]. Two kinds of mutations were further excluded from our evolved mutation list (S4 Table): double-called SNPs due to overlap between two open reading frames, and SNPs (grx3-A196V, hop1-K510fs, isu1-I100V, and rsf2-C966G) shared by multiple evolved clones that were likely to be pre-existing SNPs in the ancestral population.
For copy number variation (CNV) analysis, raw reads were mapped via BWA-MEM and CNVs were analyzed via the SAMtools program [84]. The copy number of each gene was calculated by normalizing gene coverage by the median of whole genome coverage derived from the SAMtools bedcov module. Genes with copy number variation greater than 0.8 compared with the ancestral strain were selected and the fold-change is provided in S7 Table. F1 segregant analysis E4, E5, and E7 clones were backcrossed to the ancestral strain and the diploids were sporulated. Tetrads were dissected and only spores from four-viable spore tetrads were used. The fitness of F1 haploid progeny was examined under the following conditions: YPD at 28˚C, YPD at 37˚C, and YPD + GuHCl at 28˚C. The latter two conditions were included because they had a better resolution to separate the ancestral and evolved phenotypes. F1 haploid progenies with fitness similar to the evolved or ancestral clones were assigned into good or bad spore pools, respectively. A maximum of one good and one bad spore was selected from a single tetrad. If all four progeny from a single tetrad had fitness that deviated from the evolved and ancestral strains, then none were selected. Individual E7 spores behaved similarly under the high temperature and GuHCl conditions (r = 0.851, Pearson correlation). When cells behaved differently under two conditions, we used the more prominent phenotype of the evolved clone to select the spores (i.e., fitness at 37˚C for E4 and fitness in GuHCl-containing medium for E5). We picked 5-20 spores for each pool and employed Illumina sequencing to examine the mutant allele frequency (S9 Table).
The two frequency thresholds used to select the mutations with strong effects were determined based on computational simulations following the previous method [85]. The 90%, 95%, and 99% percentiles of 10,000 times of simulations with the coverage of 200 and pool size N = 5, 9, 10, 18, and 21 were listed in S10 Table. We selected the threshold frequency that approximated to N = 9 as the general criteria for all three analyzed evolved clones.

Gene ontology enrichment and network analysis
The genes containing protein sequence mutations from all evolved Ylip-HSP90 clones were subjected to Gene Ontology (GO) analysis using Yeast GO-Slim Process embedded in Saccharomyces Genome Database (SGD: https://www.yeastgenome.org/). GO categories with only one gene were excluded. Enrichment scores (Table 1 and S8 Table) were calculated by dividing the proportion of mutated genes classified into the indicated category by the proportion of those in the genome background. Significance of enrichment was analyzed using hypergeometric test (phyper module) in R [79], and p-values were adjusted using the Benjamini Hochberg method [86]. We excluded genes with uncharacterized function or dubious open reading frames, including YLR345W, YEL034C-A, YMR102C, YLR419W, YOL166C, and Q0182.
To understand the functional relationship between the mutated genes and Hsp90, we constructed the functional network between them. The genetic and physical interactions of Hsp90 were first constructed using STRING's website (https://string-db.org/), and then each interacting pair was confirmed using the information from the SGD (https://www.yeastgenome.org/). We manually incorporated the up-and down-regulated protein-level data from a proteomics analysis [17] to make the regulatory network more complete. Gene functional groups shown in Fig 5 were derived from enriched GO categories and gene-specific description in SGD. The interaction network was constructed with both S. cerevisiae Hsp90 coding genes and all genes that contain protein-sequence-changing mutations in evolved Ylip-HSP90 clones using the cytoscape software [87].
pRS426-gRNA-AB was linearized by BamHI and EcoRI digestion. Guide RNA was amplified by two-fragment fusion PCR using four primers: IgRNAAB-F, target gene gRNA-R, target gene gRNA-F, and IgRNAAB-R (S12 Table). Donor DNA sequences were around 300-600 bp in length and were PCR amplified from the evolved clones. The original plasmid carrying the high-specificity Cas9 nuclease, eSpCas9 (1.1) [90], was obtained from Addgene (plasmid number 71814, Cambridge, MA) and the Cas9 nuclease was further subcloned into the yeast pRS414 vector in our lab. MDM12 and BRE5 mutations were constructed using the wild-type Cas9 nuclease from Dr. George Church's group [88]. CRISPR-Cas9-reconstituted clones were screened by PCR with either Donor F or R primer and an allele-specific primer that has the mutation site at the last position of the 3 0 end and a mismatched nucleotide at the third last position to decrease the tendency of nonspecific annealing [91]. Fitness data were obtained by measuring at least three independent CRISPR-reconstituted clones. In contrast, Ylip-HSP90 replacement strains showed reduced fitness compared to the Scer-HSC82 strain under most tested conditions but similar fitness in media containing 150 mM MgCl 2 , 1.8 M sorbitol, or 0.5 M NaCl. Unless specified otherwise, cells were grown at 28˚C. Cells were grown in liquid cultures and growth rates were measured by plate readers. p-values were calculated by two-tailed Student t test with the Benjamini-Hochberg correction. � p < 0.05. �� p < 0.01. ��� p < 0.001. Data sets did not deviate from the normal distribution (Shapiro-Wilk test). Error bars are standard errors, N � 4. The numerical data used in the figure are included in S1 Data. See also S1

S2 Fig. The phenotypic differences between Ylip-Hsp90-and Scer-Hsc82-carrying cells
are not due to differences in the expression level or activity of these two orthologs. (A) The protein levels of Scer-Hsc82 and Ylip-Hsp90 are not significantly different between the ancestral Scer-HSC82 and Ylip-HSP90 strains (p = 0.071, z = −1.807, two-tailed Wilcoxon rank-sum test). Total cell protein was extracted and examined by western blot. The Hsp90 protein was N-terminally fused with a TAP tag and detected by the anti-TAP antibody. G6PDH was used as the internal control. The TAP-Hsp90/G6PDH ratios were normalized to that of the Scer-HSC82 strain, and the bottom panel shows quantitative data of the western blot. Error bar represents the standard error, N � 3. (B) Scer-Hsc82 and Ylip-Hsp90 exhibit similar v-src folding activities (p = 0.827, z = −0.218, two-tailed Wilcoxon rank-sum test). Log-phase cells were grown in galactose-containing medium for 6 hours to induce the expression of v-src. The same amount of total cell protein was loaded for each sample and examined by western blot. The phospho-tyrosine and v-src were detected by anti-phospho-tyrosine and anti-v-src antibodies, respectively. The phospho-tyrosine/v-src ratios were normalized to that of the Scer-HSC82 strain, and the right panel shows quantitative data of the western blot. Error bar represents the standard error, N = 3. (C) Yarrowia lipolytica cells have higher fitness than Saccharomyces cerevisiae cells in hypersaline conditions. Serially diluted cell cultures were spotted onto plates containing 1 M NaCl or 0.2 M LiCl and incubated at 25˚C until colonies became visible. The numerical data used in panels (A, B) are included in S1 Data. G6PDH, glucose-6-phosphate dehydrogenase; Scer, Saccharomyces cerevisiae; TAP, Tandem affinity purification; Ylip, Yarrowia lipolytica. Two isogenic subpopulations labeled with red or green fluorescence protein were mixed in an approximately 1:1 ratio in the initial cultures and allowed to evolve in YPD at 28˚C. The subpopulation that obtained a beneficial mutation would start expanding its frequency in the evolved culture. The proportions of subpopulations were monitored using flow cytometry. The y-axis shows the green subpopulation frequency and the x-axis shows the numbers of evolved generations. The control line con6 and evolving lines evo1, evo2, and evo12 were excluded because of culture contamination during evolution. The numerical data used in the figure are included in S1 Data. YPD, Yeast extract-Peptone-Dextrose. (EPS)

S4 Fig. The ploidy and Ylip-HSP90 protein expression levels of evolved clones. (A) Evolved
Scer-HSC82 clones all became diploid or triploid, but most of the evolved Ylip-HSP90 clones remained haploid. All ancestral strains were confirmed to be haploid by flow cytometry. (B) Diploid Ylip-HSP90-carrying cells exhibit lower fitness than haploid Ylip-HSP90-carrying cells (p = 0.003, z = 2.744, one-tailed Wilcoxon rank-sum test), while diploid and haploid Scer-HSC82-carrying cells show similar fitness. Cells were grown in liquid cultures, and growth rates were measured by plate readers. Error bar represents the standard error, N � 4. (C) The protein level of Ylip-Hsp90 in the evolved clones is not significantly different from that of the ancestral strains (two-tailed Wilcoxon rank-sum test). G and R indicate the ancestral strains carrying the green and red fluorescent proteins, respectively. Total cell protein was extracted and examined by western blot. The Hsp90 protein was N-terminally fused with a TAP tag and detected by the anti-TAP antibody. G6PDH was used as the internal control. The TAP-Hsp90/ G6PDH ratios were all normalized to that of the ancestral Ylip-HSP90 strains, and the bottom panel shows quantitative data of the western blot. Error bar represents the standard error, N � 3. The numerical data used in the figure are included in S1 Data. G6PDH, glucose-6-phosphate dehydrogenase; Scer, Saccharomyces cerevisiae; TAP, Tandem affinity purification; Ylip, Yarrowia lipolytica.  (XLSX) S1 Data. Excel spreadsheet containing, in separate sheets, the underlying numerical data of Figs 1B, 2B, 2C, 3B, 3C, 3D, 4A, 4B, 4C, 4D, 6A, 6B, 6C, S1, S3, S5B, S6A, S6D, S8, S2A,  S2B, S4A, S4B, S4C, S6B, S6C and S7. (XLSX) mapping analyses and Jiri Koubek for assistance in strain reconstruction. We thank the Sequencing Core of Academia Sinica for sequencing services, the Genomics, Bioinformatics, Imaging, and FACS Cores of IMB for technical assistance, and John O'Brien for manuscript editing. We thank Daniel Bolon, Didier Picard, and Gavin Sherlock for sharing plasmids.