Divergence in the Saccharomyces Species’ Heat Shock Response Is Indicative of Their Thermal Tolerance

Abstract The Saccharomyces species have diverged in their thermal growth profile. Both Saccharomyces cerevisiae and Saccharomyces paradoxus grow at temperatures well above the maximum growth temperature of Saccharomyces kudriavzevii and Saccharomyces uvarum but grow more poorly at lower temperatures. In response to thermal shifts, organisms activate a stress response that includes heat shock proteins involved in protein homeostasis and acquisition of thermal tolerance. To determine whether Saccharomyces species have diverged in their response to temperature, we measured changes in gene expression in response to a 12 °C increase or decrease in temperature for four Saccharomyces species and their six pairwise hybrids. To ensure coverage of subtelomeric gene families, we sequenced, assembled, and annotated a complete S. uvarum genome. In response to heat, the cryophilic species showed a stronger stress response than the thermophilic species, and the hybrids showed a mixture of parental responses that depended on the time point. After an initial strong response indicative of high thermal stress, hybrids with a thermophilic parent resolved their heat shock response to become similar to their thermophilic parent. Within the hybrids, only a small number of temperature-responsive genes showed consistent differences between alleles from the thermophilic and cryophilic species. Our results show that divergence in the heat shock response is mainly a consequence of a strain's thermal tolerance, suggesting that cellular factors that signal heat stress or resolve heat-induced changes are relevant to thermal divergence in the Saccharomyces species.


Introduction
Temperature is a universal life parameter and organisms vary widely in their temperature tolerance.Temperature also varies over time and most organisms exhibit widespread changes in gene expression in response to temperature shifts (Lindquist 1986;Hofmann et al. 2000).In yeast, most temperature-induced changes in expression are related to a general environmental stress response, characterized by largely transient changes in genes related to protein synthesis, cellular growth, and metabolism (Gasch et al. 2000).In response to heat, cells also activate heat shock proteins (HSPs), many of which are molecular chaperones that help disaggregate and refold denatured proteins (Richter et al. 2010).Whereas the heat shock response is present in most organisms, the threshold temperature for induction varies across organisms and is correlated with their optimal growth temperature (Lindquist 1986;Feder and Hofmann 1999).These and other observations implicate the heat shock response as a trait important to the ecology and evolution of most organisms (Tomanek 2010).
The heat shock response is an important mechanism of thermal tolerance (Richter et al. 2010).However, the fitness effects of the heat shock response are complex.In yeast, only a small fraction of proteins that are induced by heat are required to survive thermal stress (Gibney et al. 2013).The heat shock response also generates acquired thermotolerance whereby a mild heat stress increases the later survival of a more severe heat stress.Acquired stress tolerance can also be activated by and protective against a variety of stresses besides heat (Berry et al. 2011).Both a cell's stress response and its resistance to stress can depend on its growth rate, cell cycle phase, and metabolic state (Lu et al. 2009;Slavov et al. 2012).
The heat shock response has also been associated with thermal adaptation within and between species.Early work found that both the constitutive expression of HSPs and the temperature at which the heat shock response is activated are associated with thermotolerance (Feder and Hofmann 1999).Recent work has expanded upon these findings to implicate changes in HSP copy number and transposable element-mediated changes in HSP expression (Chen et al. 2018).In yeast, experimental evolution at high temperature resulted in constitutively activated HSPs, even at low temperatures (Satomura et al. 2013;Caspeta et al. 2016).Thus, there is strong support for investigating the role of the heat shock response in the evolution of thermotolerance.
The Saccharomyces yeast species exhibit substantial differences in their thermal growth profiles and provide an advantageous system to examine thermal adaptation (Gonçalves et al. 2011;Salvadó et al. 2011;Baker et al. 2019).Both Saccharomyces cerevisiae and Saccharomyces paradoxus can grow at 37°C, and, at higher temperatures, S. cerevisiae grows better than S. paradoxus.In comparison, Saccharomyces uvarum, Saccharomyces eubayanus, and Saccharomyces kudriavzevii cannot grow at 37°C, but they grow faster than S. cerevisiae at low temperatures.These thermal growth differences impact competitive growth differences and are particularly relevant to wine making (Williams et al. 2015;Alonso-del-Real et al. 2017).Two of the most thermally divergent species, S. cerevisiae and S. uvarum, also differ in their ability to survive heat shock treatment (Li 2018).
A number of studies have investigated the genetic basis of thermal divergence in the Saccharomyces species.Glycolytic enzymes have divergence in their activity and thermal stability (Gonçalves et al. 2011).The mitochondrial genome contributes to growth differences at high and low temperatures (Baker et al. 2019;Li et al. 2019).A genomewide screen identified eight genes involved in cell division and growth that provide S. cerevisiae with higher thermal tolerance relative to S. paradoxus but compromise growth at low temperature (Weiss et al. 2018;AlZaben et al. 2021).In S. cerevisiae and S. uvarum hybrids, there are abundant cis-acting differences in gene expression, but relatively few of these differences depend on the temperature at which the hybrid was grown (Li and Fay 2017).
Thermal growth differences among the Saccharomyces species are typically captured by monitoring growth after cells are shifted to a new temperature.The goal of this study was to examine the initial heat and cold shock response prior to any growth differences in order to determine whether and how these species have diverged in their temperature response.Species' hybrids show that thermal tolerance is dominant, raising the question of whether thermally protective alleles mediate their effects prior to or subsequent to temperature changes and the heat shock response.Our results identify a small number of genes consistent with protective effects and show that hybrids activate a stronger stress response than their thermophilic parent but resolve their response more rapidly than their cryophilic parent.

Results
To examine divergence in heat and cold shock responses, gene expression was measured at four time points (0, 15, 30, and 60 min) after a shift from 25 °C to a low (12 °C) or high (37 °C) temperature for four Saccharomyces species and their six interspecific hybrids.Gene expression was measured using species-specific read counts from RNA sequencing and alignment to complete genomes of each species.Because only a draft genome with annotations was available for S. uvarum, we generated a complete chromosome-level assembly along with gene annotations (see Materials and Methods) and identified 5,121 genes with one-to-one orthologs in each of the four species.

Thermal Response in Parental Species
Across the four species, more genes showed altered expression in response to heat (4,639) than to cold (1,874) (false discovery rate [FDR] < 0.01, table 1), consistent with prior work in S. cerevisiae (Sahara et al. 2002;Schade et al. 2004).The average thermal response was similar in the two thermophilic species, S. cerevisiae and S. paradoxus, whereas it was initially smaller or slower in S. kudriavzevii in response to both temperature shifts and stronger in S. uvarum in response to heat (fig.1).
To examine species differences in their heat and cold shock response, we used 3,125 heat-responsive and 445 cold-responsive genes with a significant interaction between the species' alleles and time (FDR < 0.01, table 1).Using expression values normalized to the initial time point, we identified genes with large expression differences between the two thermophilic and the two cryophilic species.In response to heat, there were 197 genes whose average expression was 2-fold higher in the cryophilic compared with the thermophilic species.These genes were significantly enriched for involvement in de novo posttranslational protein folding (n = 8, P = 7.0e −5 , supplementary table S1, Supplementary Material online) and showed increased expression in response to heat and a larger increase in the cryophilic species (fig.2A).There were 50 genes whose average expression was 2-fold lower in the cryophilic species compared with thermophilic species.These genes were significantly enriched for involvement in cytoplasmic translation (n = 14, P = 3.5e −10 , supplementary table S1, Supplementary Material online) and showed lower expression in response to heat in the cryophilic species (fig.2B).In response to cold shock, there were only 11 genes with a 2-fold difference between the thermophilic and cryophilic species.Together, these results show that the cryophilic species exhibit a stronger heat shock response than the thermophilic species for a subset of genes.
Thermal Response in Species' Hybrids Gene expression differences between species can be caused by any number of cellular factors involved in the detection and response to thermal treatment.All the hybrids except S. kudriavzevii × S. uvarum (Sk × Su) grow at 37 °C (Li et al. 2019).We thus examined whether hybrids with either an S. cerevisiae or S. paradoxus parent adopted a heat response similar to their parents.
The hybrids showed a large but variable number of genes that responded to temperature (FDR < 0.01, table 1).Consistent with hybrid thermal growth profiles, hybrids with an S. cerevisiae parent had fewer genes that changed in response to heat or cold, and the S. cerevisiae × S. paradoxus (Sc × Sp) hybrid had the fewest heat-responsive NOTE.-Time, allele, and time * allele indicate the number of genes with a significant effect of time point, species' allele, and an interaction between time and species' allele, respectively.

FIG.
1.-Species' response to heat and cold treatment.Genes that significantly changed over time in response to heat or cold treatment were normalized to the zero time point and shaded areas show the 50%, 90%, and 95% quantiles (from darkest to lightest) of log2 changes in gene expression.Species are labeled with the first letter of their genus and species name: Sc, Sp, Sk, and Su.
Divergence in the Saccharomyces Species' Heat Shock Response GBE Genome Biol.Evol.15(11) https://doi.org/10.1093/gbe/evad207Advance Access publication 16 November 2023 genes of all the hybrids.To facilitate comparisons among hybrid gene expression, we used principal component analysis (PCA) of all genes regardless of which hybrids showed significant time-dependent changes.
For the heat shock response, the first component explained 58% of the variance and the remaining components each explained less than 6% of the variance.The first component was associated most strongly with time point, followed by strain and allele (analysis of variance [ANOVA] R 2 of 83%, 5.1%, and 0.9%, respectively).Using the first component, we found the hybrids differed from one another and displayed a mixture of parental phenotypes that changed over time (fig.3).At 15 min, S. uvarum showed the strongest response and S. kudriavzevii showed the weakest response among the parental species.At 15 min, the hybrids were all more similar to S. uvarum, including the Sc × Sp hybrid.At 30 min, the response of S. kudriavzevii and the S. cerevisiae × S. kudriavzevii (Sc × Sk) hybrid was noticeably stronger than at 15 min.Except for the Sc × Sk hybrid, all the hybrids with a thermophilic parent showed a reduced response at 30 compared with 15 min, and the S. cerevisiae and S. paradoxus hybrid response was weakest and similar to both of its parents.Finally, at 60 min, most of the hybrids showed a reduced response, similar to that of S. cerevisiae and S. paradoxus.The two exceptions to this pattern were the Sc × Sk and Sk × Su hybrids that retained their response level at or above S. kudriavzevii and S. uvarum.Together, these results show a strong heat shock response in the hybrids that is more rapidly resolved for hybrids with a thermophilic parent.
For the cold shock response, the first component accounted for 31% of the variance and was most strongly associated with time point (ANOVA R 2 = 77%).The overall pattern of the first component differed from the heat shock response (fig.3).Like S. kudriavzevii and S. uvarum, all the hybrids showed an increasing response over time, with the A B FIG. 2.-Genes that exhibit a stronger heat shock response in S. kudriavzevii and S. uvarum compared with S. cerevisiae and S. paradoxus.(A) Eight genes involved in de novo posttranslational protein folding with a 2-fold higher average response in the cryophilic species.(B) Fourteen genes involved in cytoplasmic translation with a 2-fold lower average response in the cryophilic species.60-min response being at or above the 30-min response.Across all the time points, the Sc × Sk hybrid showed the weakest response, similar to its parental species.Thus, unlike the heat shock response, we found no association between cryophilic species and the cold shock response.
An alternative way of examining dominance in the hybrid is to focus on those genes that differ among the parental species.We thus considered genes that were 2-fold higher (197) or 2-fold lower (50) in the cryophilic species and examined their expression in the hybrids.Normalized to the initial time point, the average heat response in the thermophilic species hybrid (Sc × Sp) compared with the cryophilic species hybrid (Sk × Su) often mimicked differences found in their parents (fig.4).Genes that were 2-fold or higher in the thermophilic species (Sc and Sp) were also higher in the thermophilic hybrid (Sc × Sp) compared with the cryophilic hybrid (Sk × Su).Similarly, genes that were 2-fold lower in the thermophilic parents were also lower in the thermophilic hybrid.In comparison, these same genes did not show the same pattern in comparisons of hybrids containing one thermophilic parent and one cryophilic parent (Sc × Sk, S. paradoxus × S. kudriavzevii [Sp × Sk], S. cerevisiae × S. uvarum [Sc × Su], and S. paradoxus × S uvarum [Sp × Su], fig.4).In comparison with the thermophilic (Sc × Sp) or cyrophilic (Sk × Su) hybrid, these hybrids displayed intermediate phenotypes (supplementary fig.S1, Supplementary Material online).These results show that similar to all genes, genes with a stronger heat shock response in the cryophilic species show the strongest response in the Sk × Su hybrid and the weakest response in the Sc × Sp hybrid.

Allele Differences in Species' Hybrids
Whereas each hybrid showed a different temperature response, within each hybrid both species' alleles experience the same cellular environment and enable cis-acting regulatory divergence to be quantified.We examined allele-specific expression differences within each hybrid separately.Compared with their parents, the hybrids showed fewer allele differences that changed over time.There was an average of 256 heat-responsive and 11 cold-responsive genes that exhibited significant interactions between the species' alleles and time (FDR < 0.01, table 1).Consistent with the number of genes with Divergence in the Saccharomyces Species' Heat Shock Response time-dependent changes in expression, the three hybrids with an S. cerevisiae parent had fewer genes with significant allele-time interactions in response to both the heat and cold treatment (table 1).Although the hybrid of the two thermophilic species, S. cerevisiae and S. paradoxus, showed the fewest allele-time interactions in response to heat, these two species are also the most closely related.
We next examined whether there was any consistent relationship between genes with allele differences and how they responded to heat.Only two of the hybrids (Sc × Sp and Sk × Su) showed an association between allele differences and heat response based on the number of genes with significant allele and time effects compared with all other genes (chi-squared test, Bonferroni-corrected P < 0.05, supplementary fig.S2, Supplementary Material online).Four hybrids (all but Sc × Sp and Sc × Sk) showed an association based on the number of genes with significant allele-time interactions.However, the major trend for this later association was that genes with allele-time interactions were more likely to increase in response to heat and not that alleles of thermophilic species were expressed higher than those from cryophilic species for either heat-induced or heat-repressed genes.
We also assessed whether there were any genes with a significant allele-time interaction that showed a consistent difference between the alleles of the thermophilic and cryophilic species.Most of the genes were significant in only one of the hybrids (775/1,095).Of the remaining, 208 were significant in at least two of the four hybrids between thermophilic and cryophilic species, and five genes were significant in all four of the hybrids (fig.5: TDH1, UTR4, MCR1, PER33, and DAP1).All of these five genes showed a consistent expression difference between the thermophilic and cryophilic species' alleles.

Thermal Response of Gene Families
In addition to the 5,121 one-to-one orthologs, we also examined 321 gene families where the number of family members varied across species or where gene conversion between family members obfuscated one-to-one ortholog assignment.The gene families include not only many duplicated ribosomal genes but also HSPs (HSC82/HSP82, SSB1/SSB2, and HSP31-34) and cell wall mannoproteins (DAN/TIR/PAU).To accommodate the variable number of family members, we used the summed expression of all family members within a species.For these gene families, the fraction that showed time dependence and allele-time interactions was similar to that of one-to-one orthologs (supplementary table S2, Supplementary Material online).Only eight of these genes showed significant allele-time interactions between two of the four thermophilic-cryophilic hybrids.One of these families encodes 2-deoxyglucose-6-phosphate phosphatases (Randez-Gil et al. 1995) and has two copies in S. cerevisiae and S. paradoxus (DOG1/DOG2), but only one copy in S. kudriavzevii and S. uvarum.Whereas the S. kudriavzevii and S. uvarum alleles were strongly induced upon heat shock, the S. cerevisiae and S. paradoxus alleles were not (supplementary fig.S3, Supplementary Material online).

Discussion
The Saccharomyces species differ in their thermal growth profiles.These growth differences are rendered subsequent to any heat or cold shock response that occurs upon shift to a high or low temperature.Our study of gene expression in response to temperature shows that the Saccharomyces species and their hybrids differ substantially in their response to temperature, particularly heat.We find that major components of expression divergence are predicted by a species' thermal tolerance whereby sensitive species exhibit the strongest heat shock response.This pattern implies that the heat shock response is a signal of heat stress rather than FIG.4.-Hybrids recapitulate parental strain differences.Hybrid expression differences are shown in comparison to genes that differed two-fold between thermophilic (Sc & Sp) and cryophilic (Sk & Su) species.For the same genes, differences between the thermophilic (Sc × Sp) and cryophilic (Sk × Su) hybrid recapitulate many of the parental differences, whereas the two S. kudriavzevii hybrids (Sp × Sk and Sc × Sk) and the two S. uvarum hybrids (Sp × Su and Sc × Su) did not.Expression difference (log2) is the average response to heat across time points and species' alleles for each comparison (above and below the horizontal lines).Numbers at the top indicate the number of genes with differences greater than two-fold up/ down.a driver of thermotolerance.However, most hybrids with a thermophilic parent are not entirely protected and initially show a stronger stress response that is only later reduced to levels present in their thermophilic parents.Below, we discuss these results, deviations found in the Sc × Sk hybrid, and genes/pathways that may play a role in thermal divergence of the Saccharomyces species.
Multiple aspects of the heat shock response are associated with growth at high temperature.At 37°C, S. kudriavzevii, S. uvarum, and their hybrid do not grow, but there is also variation in growth among the other species and hybrids (Li et al. 2019).S. cerevisiae and the Sc × Sp hybrid produce the largest colonies.The two other S. cerevisiae hybrids produce slightly smaller colonies, with Sc × Su being larger than Sc × Sk. S. paradoxus colonies and its hybrids (Sp × Sk and Sp × Su) are of equivalent size and smaller than those of S. cerevisiae and its hybrids.These growth differences mirror multiple aspects of the heat shock response.First, fewer genes showed significant responses to heat and allele interactions in hybrids with an S. cerevisiae parent and the fewest in the Sc × Sp hybrid (table 1).Second, the PCA analysis showed that the strength of the heat shock response at the later time points corresponded to the heat sensitivity of the hybrids, with the Sc × Sk hybrid being one notable exception.Likewise, genes that were more strongly induced (protein folding) or repressed (ribosomal) in the cryophilic species were also more responsive to heat in the Sk × Su compared with the Sc × Sp hybrid (fig.4) and showed an intermediate response in the other hybrids (supplementary fig.S1, Supplementary Material online).
Our findings are similar to a prior study of Saccharomyces hybrids that found a correspondence between heat shock response and growth at high temperature (Kempf et al. 2017).Using Hsp104-foci as a marker for protein aggregation, they found that an Sc × Sp hybrid resolved Hsp104-foci more rapidly than the other S. cerevisiae hybrids.Additionally, they found that an Sc × Sk hybrid showed Hsp104-foci at a lower temperature, and these foci reverted the slowest among all the hybrids.
In our study, interpreting expression of the Sc × Sk hybrid is not straightforward.The Sc × Sk hybrid, like its parents, is a uracil auxotroph.Although we did not observe substantial growth defects in Chardonnay grape juice, prior to heat shock, the Sc × Sk hybrid and its parents showed higher expression of genes involved in uracil biosynthesis (URA1, URA2, and URA8) and lower expression of genes involved in lysine biosynthesis (LYS4, LYS9, LYS12, LYS20, and LYS21).Equivocally, the hybrid also showed expression differences in glutathione metabolism (supplementary fig.S4, Supplementary Material online).Glutathione metabolism is of interest since glutathione can alleviate temperatureinduced oxidative stress (Morano et al. 2012), and secretion of glutathione has been shown to elevate the maximum temperature at which populations of S. cerevisiae can grow in a density-dependent manner (Laman Trip and Youk 2020).In our experiments, glutathione synthetase (GSH1) was expressed at lower levels in S. kudriavzevii and the Sc × Sk hybrid compared with all other strains.The major glutathione-degrading enzyme ECM38 showed higher expression of the S. kudriavzevii allele in S. kudriavzevii and all of its hybrids.For both of these genes, the S. cerevisiae allele was more similar to S. kudriavzevii than S. paradoxus, making the Sc × Sk hybrid different from the Sp × Sk hybrid.Additionally, GLR1, which converts glutathione to its reduced form, showed greatly reduced expression in response to heat, but only in S. cerevisiae and for the S. cerevisiae allele in the Sc × Sk hybrid.The phospholipid hydroperoxide glutathione peroxidase (GPX2) uses glutathione to neutralize reactive oxygen species and showed low expression of the S. kudriavzevii allele in all the hybrids.Although uncertain as to the cause, the Sc × Sk hybrid showed both altered metabolism prior to heat shock and a heat shock response that stands out from the other hybrids.
We identified five genes that showed consistent allele differences between the thermophilic and cryophilic species.While none of them are directly linked to thermal tolerance, two are important for Erg11 activity and ergosterol biosynthesis.MCR1 encodes a mitochondrial NADH-cytochrome b5 reductase that is one of two electron donors to the cytochrome P450 Erg11 and potentially Erg5 and Erg1 (Lamb et al. 1999).DAP1 is a heme-binding protein that positively regulates Erg11 and Erg3 activity (Craven et al. 2007).In S. cerevisiae, inhibition of the later steps in ergosterol biosynthesis (Erg3-5) improves thermal tolerance, potentially through altered sterol composition of membrane lipids (Caspeta et al. 2014;Liu et al. 2017;Yang et al. 2023).However, the expression patterns of MCR1 and DAP1 may be incidental as they show opposite (compensatory) patterns in the thermophilic species.
Of the remaining three genes, TDH1 functions in glycolysis and gluconeogenesis, UTR4 is involved in the methionine salvage pathway, and PER33 is associated with the ER and nuclear pore complexes.One interesting facet of PER33 is that across many experimental conditions (SGD SPELL v2.0.3) (Hibbs et al. 2007), PER33 expression is most correlated with HCH1 (HSP regulator), CPR6 (chaperone activity with Hsp82), and CUR1 (protein quality control); and its human homolog TMEM33 is involved in the unfolded protein response (Sakabe et al. 2015).
Compared with heat shock, there were fewer genes that changed in response to cold treatment.This is consistent with prior work showing a small early response (<2 h) and a larger later response (>2 h) that is related to the environmental stress response (Schade et al. 2004).We found increased expression of genes previously described to be induced by cold treatment (NSR1, TIP1, BFR2, OLE1, and TIR1/2/4).Similar to prior studies of log phase expression in hybrids (Li and Fay 2017;Hovhannisyan et al. 2020;Timouma et al. 2021), we found an appreciable number of allele differences but very few allele differences that changed in response to cold treatment.
Relevant to both our heat and cold shock results, prior work has shown that gene expression in interspecific hybrids depends on which mitochondria they inherited (Solieri et al. 2008;Hewitt et al. 2020).The hybrids used in this study had mitochondria from their more thermophilic parent.The growth differences among hybrids depend on their mitochondria DNA and growth differences are mostly absent in petite hybrids lacking mitochondrial DNA (Li et al. 2019).It thus seems likely that the gene expression differences we observed may be quite different if the mitochondria were from the less thermophilic parent.
Gene families present a challenge to examining expression divergence when family members are closely related.Their analysis also depends on accurate annotation of complete genome sequences.For gene family members that could not be assigned into one-to-one orthologous groups, for example, varied in copy number among species, we used the sum of their expression values within each species.One interesting example is DOG1/2, where the singlecopy DOG1/2 alleles in S. kudriavzevii and S. uvarum showed a much stronger heat response than the two copies in S. cerevisiae and S. paradoxus (supplementary fig.S3, Supplementary Material online).Mutations in DOG1/2 confer 2-deoxyglucose resistance as Dog1/2 metabolizes derivatives of the glucose analog 2-deoxyglucose (Sanz et al. 1994).Although there are no direct links between DOG1/ 2 and thermotolerance, the top DOG1 phenome hit is sensitivity to fluvastatin, a statin that inhibits HMG-CoA reductase and ergosterol biosynthesis (Turco et al. 2023).
Among the other gene families, three are HSP families.Whereas SSB1/2 and HSC82/HSP82 are present as closely related duplicates in each of the four species (consistent with duplication prior to speciation accompanied by gene conversion within each species), the HSP31-34 family varies in copy number.HSP31 is present in all four species as one-to-one orthologs, increases in response to heat, and shows much higher expression in the two thermophilic species (supplementary fig.S3, Supplementary Material online).HSP32 is only present in S. cerevisiae, HSP33 is only present in S. cerevisiae and S. paradoxus, and HSP34 is likely to be a very recent duplicate as it is not present in any of the genomes we analyzed.Members of the HSP31 gene family have methylglyoxalase activity and redox-dependent chaperone activity that acts early in protein misfolding (Miller-Fleming et al. 2014;Amm et al. 2015;Tsai et al. 2015;Bankapalli et al. 2020).Deletions of HSP31 family members cause reduced thermotolerance, sensitivity to oxidative stress, altered mitochondrial homeostasis, and an increase in reactive oxygen species and oxidized glutathione (Miller-Fleming et al. 2014;Bankapalli et al. 2020).These studies make the HSP31 gene family a good candidate for functional studies of thermal divergence between species.
In summary, we find that thermal tolerance is indicative of the thermal response of the Saccharomyces species and their hybrids.What is the signal that turns on the heat shock response, and why is it stronger in the thermally sensitive species?The heat shock transcription factor, HSF1, and two stress response transcription factors, MSN2 and MSN4, are responsible for most of the heat shock response (Morano et al. 2012).Phosphorylation of Msn2/4 by growth-related signaling pathways (cAMP/PKA) blocks their import to the nucleus and inhibits their activation of the stress response.Hsf1 is also posttranslationally modified upon heat shock, and under the chaperone titration model, Hsf1 is negatively regulated by Hsp70, which is recruited away upon proteotoxic stress leading to Hsf1 activation of the heat shock response (Masser et al. 2020).However, heat shock also causes intracellular acidification and the formation of stress granules.Protein constituents of stress granules condense by phase separation in a temperature-and pH-dependent manner that may constitute an adaptive response to limit translation and protein misfolding (Glauninger et al. 2022).Two proteins, Pab1 and Ded1, are known to assemble into condensates in a temperature-and species-dependent manner such that the S. kudriavzevii protein does so at lower temperatures than that of S. cerevisiae (Iserman et al. 2020;Kik et al. 2023).Along with other thermal sensing proteins, there is growing evidence that intrinsically disordered regions can function as cellular sensors of temperature (Riback et al. 2017;Jung et al. 2020;Moses et al. 2023).Thus, further characterization of the Saccharomyces species' response to temperature could yield insight into the regulation of the heat shock response as well as to how it evolves.

Strains
A total of six diploid hybrid strains were made by crossing haploid strains of four Saccharomyces species (supplementary table S3, Supplementary Material online).The parental strains are derivatives of the following strains: T73 (S. cerevisiae), originally isolated from a vineyard in the Alicante wine region of Spain (Querol et al. 1992) and used as a commercial wine strain (Lallemand Inc., Montreal, Canada); N17 (S. paradoxus) originally isolated from oak exudate in Tatarstan (Glushakova et al. 2007) and a member of the European population (Liti et al. 2009); CR85 (S. kudriavzevii) isolated from an oak tree bark in Agudo, Ciudad Real, Spain (Lopes et al. 2010); and BMV58 (S. uvarum) a commercial wine strain (Lallemand Inc.), originally isolated from wine in the Utiel-Requena region of Spain (Henriques et al. 2021).Haploid derivatives of different mating types were obtained by replacing the HO gene with a dominant drug-resistant marker (Goldstein and McCusker 1999;Vorachek-Warren and McCusker 2004), sporulation, and selecting haploid progeny by mating-type polymerase chain reaction (PCR) (Huxley et al. 1990).Diploid hybrids with known mitochondrial types were generated by crossing parental strains with one parent lacking mitochondrial DNA (Li et al. 2019).

RNA Sequencing
Gene expression was measured by RNA sequencing of hybrids and their parental strains subjected to either a heat or cold shock.Two replicates of the six hybrids and the four parental strains were placed into 5-mL rich medium (YPD: 1% yeast extract, 2% peptone, and 2% dextrose) and grown overnight at 25°C.Out of each of these precultures, two sets (A and B) of four 2-mL tubes containing Chardonnay grape juice (Vintner's Reserve Chardonnay, 10L, Cultures for Health) at 25°C were inoculated at an OD 600 of 0.3 and placed at 25°C and 100 rpm for 5 h.For thermal shock, sample sets A were placed at 12°C and sets B at 37°C.Sample collection for gene expression profiling was done at four different times: 0, 15, 30, and 60 min after thermal treatment.The total number of samples was 160: ten strains, four time points, two temperatures, and two replicates.To reduce the number of sequencing libraries, cells from different tubes were mixed in order to bring the number of samples from 160 down to 64 (supplementary table S4, Supplementary Material online), with the criterion that none of the mixed strains shared the same genomic background (e.g., Sc × Sp and Sp × Su).RNA samples were enriched for mRNA with the NEBNext Poly(A) mRNA Magnetic Isolation Module (NEB: New England Biolabs, Ipswich, MA, USA), and library preparation was carried out by RNA fragmentation, cDNA synthesis, adaptor ligation, size selection, and PCR enrichment using the NEBNext Ultra Directional RNA Library Prep Kit (NEB).PCR enrichment step was performed with NEBNext Multiplex Oligos for Illumina (Dual Index Primers, NEB #E7600), using 64 different combinations of indexed primers.A pool of the 64 libraries was run twice on a NextSeq Sequencing System from Illumina (2 × 75 bp run), generating a total of 936 million paired reads.

Gene Expression
For read mapping, we selected representative genome assemblies from NCBI based on completeness and similarity to the strains used in this study.For S. cerevisiae, we selected DBVPG6765 (GCA_002057805.1)(Yue et al. 2017), a European wine strain with Chr15R elongated by a 65-kb horizontal gene transfer from Torulaspora microellipsoides (Marsit et al. 2015).For S. paradoxus, we used CBS432 (GCA_002079055.1)(Yue et al. 2017), a common reference for European strains (Liti et al. 2009).For S. kudriavzevii, we used CR85 (GCA_003327635.1),and for S. uvarum, we generated a genome assembly of YJF1449, a derivative of CBS7001 (described below).All assemblies contained complete chromosome-level contigs with no gaps except for CR85, which had a six-gene gap between PMT6 and ADE3, and an internal gap in ATG13, and DBVPG6765 and CBS432, which had their ribosomal RNA genes removed.
Reads were mapped to a combined reference of the four species' genomes using TopHat (v2.1.1)(Kim et al. 2013) with a minimum and maximum intron length of 40 and 1500, respectively.HTSeq-Count (v0.9.1) (Anders et al. 2015) was used to count unique reads mapping to gene features using a minimum alignment quality of 10.Over all the samples, 78% of read pairs were included yielding a median of 2.1 million counts per species-allele/sample, with a range of 82 thousand to 16 million after removing one sample with very few reads (S26, S. paradoxus, cold, 15 min, supplementary table S4, Supplementary Material online).The median number of species-specific counts per gene was 125 for the 5,121 one-to-one orthologs (see Orthology Assignments).

Genome Sequencing
Two S. uvarum strains were used for genome assemblies: YJF1449 and YJF4719, both derivatives of CBS 7001 (supplementary table S3, Supplementary Material online).High molecular weight DNA was extracted (Barbitoff et al. 2021) with some modifications.Strains were grown in 500-mL YPD overnight at 25 °C in a shaking incubator at 300 rpm.After pelleting (2,500 × g for 5 min) and washing twice with DNAse-free distilled water, pellets were resuspended in 10 mL sorbitol lysis buffer (1 M sorbitol, 0.05 M EDTA, and pH 8.0).Cells were treated with 100-µL zymolase (25 mg/mL, 20 T) and incubated for 1 h at 37 °C to make spheroplasts.Spheroplasts were pelleted and resuspended in 10-mL total lysis buffer (0.1 M NaCL, 0.01 M Tris-HCl, 0.025 M EDTA, and 0.5% w/v sodium dodecyl sulfate) and treated with 50-µL RNAse A (10 mg/mL) for 2 hours at 37 °C.Lysates were treated with 100-µL proteinase K (20 mg/mL) overnight at 50 °C.The next day, 15-mL 25:24:1 phenol:chloroform:isoamyl alcohol pH 8.0 (PCIaA) was added to each lysate and incubated for 10 min on a rotator.Lysates were added to phase separation tubes (50-mL conical tubes with 2-mL Dow Corning high vacuum grease prespun to the bottom of the tubes) and centrifuged at 4,500 × g for 5 min.The upper phase was transferred to new conical tubes where a second round of PCIaA extraction was performed.DNA was then precipitated with 1/10th volume 2 M NaCl followed by 2/3 volume of isopropanol and mixed by gentle inversions after each addition.After a 10-min incubation at room temperature, a glass hook was used to capture and place the DNA into 1.5-mL 70% ethanol, followed by 100% ethanol, and then dried.Precipitated DNA was rehydrated in 1-mL TE for 3 days.Size selection was performed using the Short Read Eliminator Kit (Circulomics).
Long-read DNA libraries were prepared using the Ligation Sequencing Kit (Oxford Nanopore SQK-LSK110, version GDE_9109_v110_revJ_10Nov2020) without the addition of lambda phage control DNA and were run through an R9.4.1 flow cell on a MinION.Each strain was sequenced on its own flow cell over which multiple rounds of sequencing libraries were run for 18-24 h each.Between sequencing runs, the flow cell was cleaned and refreshed using the Flow Cell Wash Kit (EXP-WSH004, version SFC_9120_V1_revB_ 08Dec2020).Flow cells were reused until less than 5% of the pores were deemed available after washing.Basecalling was performed using Guppy v.5.0.15 (Oxford Nanopore,dna_r9.4.1_450bps_sup.cfg configuration).
A second DNA extraction was performed for short-read libraries using a YeaStar Genomic DNA Kit (Zymo Research).Libraries were generated using the Nextera DNA Flex Library Prep Kit (Illumina) with unique dual indexes (IDT for Illumina) at 1/3 reaction volumes.Paired-end (2 × 100) sequencing was performed on a HiSeq (Illumina) by Novogene.
The filtered long reads were assembled with Canu v.2.2 (Koren et al. 2017) and then polished with all long reads that passed initial basecalling with Medaka v.1.5.0 (Oxford Nanopore) and then with five rounds of short-read sequences using Polypolish v.0.5.0 (Wick and Holt 2022).The resulting assemblies were oriented and named according to the chromosome numbers assigned in the S. uvarum genome assembly for strain CBS 7001 (Chen et al. 2022) (NCBI accession number JAIPTR010000000).Two modifications were made to the YJF1449 assembly.One of the YJF1449 contigs contained a concatenated repeat of the complete mitochondrial sequence.This contig was oriented and trimmed to a single iteration by aligning to an existing S. uvarum mitochondrion sequence (NCBI accession number KX657742).On chromosome XII, the YJF1449 assembly differed from that of YJF4719 and JAIPTR010000000 at the repeated rRNA locus.This region was thus reassembled and incorporated into YJF1449 chromosome XII using Gap5 v.1.2.14-r (Bonfield and Whitwham 2010).Reassembly was performed by extracting all Q10 and higher long reads in the region and assembling with Canu.Unlike the original assembly, the reassembled region showed even read coverage, similar to that obtained by mapping to the YJF4719 and JAIPTR010000000 assemblies, indicating that the reassembly was more congruent with the input read data.The number of reads, bases, and assembly statistics were tabulated for both strains (supplementary table S5, Supplementary Material online).

Genome Annotation
We generated genome annotations for each species.Our rationale was that we could use S. cerevisiae and S. paradoxus as controls to assess the quality and completeness of our annotation methods since they were previously accurately annotated using the LRSDAY pipeline (Yue et al. 2017;Yue and Liti 2018).
Augustus (v3.4.0) and YGAP were used to predict proteincoding genes in the nuclear genome.We were unable to install and generate LRSDAY predictions for S. kudriavzevii and S. uvarum due to software dependencies.Augustus was run with the saccharomyces_cerevisiae_S288C model, softmasking, hint parameters (extrinsic.M.RM.E.W.P.cfg), and with flags to report coding sequences with stop codons in gff3 format.We used RNA-sequencing alignments from TopHat (above) and scripts from Augustus to filter for unique hits (filterBam -uniq) and generate intron hints (bam2hints with max gene length = 10 kb).A preliminary assessment of Augustus predictions using intron hints did not show increased concordance with LRSDAY annotations of DBVPG6765, and so we only used RNA-sequencing hints for manual curation and not for Augustus gene predictions.To generate protein hints, we used Spaln2 (v2.4.1) (Iwata and Gotoh 2012) to map S. cerevisiae proteins (5,943 genes annotated in S288c version R64-3-1 after excluding dubious genes) to each genome (max gene size = 10 kb).Hits were filtered with the align2hints script from Augustus (max intron size 1,500) and yielded hints for 5,909 S. cerevisiae, 5,756 S. paradoxus, 5,555 S. kudriavzevii, and 5,471 S. uvarum proteins.For each predicted gene, IDs were designated by Augustus and a descriptive name was obtained by the best hit of the coding (DNA) sequence to all S288C coding genes using Blastn (-word_size 9, -evalue 1e-40, -strand plus), or by the best Blastp hit if no cDNA hit was returned.YGAP annotations were obtained online (Post-WGD species, No frameshift correction).GenBank output files were converted to gff using bioperl's bp_genbank2gff3 script.

Comparison of Gene Predictions and Manual Curation
Augustus, YGAP, and LRSDAY annotations for S. cerevisiae and S. paradoxus were compared.Most annotations (86-90%) were exact matches across all three sets, with the same start/stop and intron/coding-exon positions (supplementary table S6, Supplementary Material online).Most of the remaining annotations were exact matches in two out of the three sets and involved slight differences in start/stop and intron positions.To assess the YGAP and Augustus annotations, we used manual curation based on protein and RNA-sequencing hints (described above) and gene structure in S288c (SGD, R64-3-1) including gene and intron size.The majority of annotations with exact matches in two of the three sets were curated as correct and we thus found clear misannotations in each set including missing introns and merged or split genes.In general, Augustus did better at annotating genes with introns and YGAP did better at annotating small genes.We only identified 17 S. cerevisiae and 31 S. paradoxus annotations that were curated as correct in LRSDAY but incorrect in both YGAP and Augustus (supplementary table S6, Supplementary Material online).Based on these results, we curated all YGAP and Augustus annotations that did not match each other in S. kudriavzevii and S. uvarum.For each species, we generated a final set of annotations (supplementary table S6, Supplementary Material online) by combining our curation of nuclear proteins with annotations from tRNAscan and MFannot and repeats from RepeatMasker.For consistency in the gene expression analysis, we used these annotations rather than those from LRSDAY for S. cerevisiae and S. paradoxus.

Orthology Assignments
Hierarchical orthologous groups (HOGs) were defined by OrthoFinder (v2.5.4) (Emms and Kelly 2019).Hierarchical orthogroups are defined at each node in the species tree and OrthoFinder found 5,442 HOGs at the ancestral (N0) node, of which 5,270 had at least one ortholog in each species and 5,121 had one-to-one orthologs in all four species (supplementary table S7, Supplementary Material online).

Expression Analysis
Species-specific read counts were normalized using DESeq2 (Love et al. 2014).For each treatment (heat shock and cold shock), count data were fit to a model with interacting time and species-allele terms: counts ∼ time * allele.Genes with significant changes over time, between alleles, or the interaction between the two were found by dropping those respective terms and evaluation using the negative binomial likelihood ratio test with a FDR cutoff of 0.01 (Benjamini and Hochberg 1995).For 321 gene families that did not have one-to-one orthologs in all four species, we used the sum of the read counts across family members within each species.
For parental species differences, genes showing a significant interaction between allele and time were analyzed.For these genes, the average difference between the thermophilic and cryophilic species was obtained from the species and time point average after each gene was normalized to the initial time point.Gene ontology term analysis was conducted using SGD with a background set equal to all 5,121 genes in the dataset and Holm-Bonferroni correction for multiple tests.
PCA was used to compare the hybrid and parental species response to heat and cold treatment.Each gene's expression was normalized to the zero time point, and principal components were obtained from the centered but unscaled expression data obtained from DESeq2's variance stabilizing transformation function.
FIG. 3.-Hybrid and parental strains exhibit different heat and cold shock responses.Bars show the mean of the first principal component and whiskers show the standard error from replicates and species' alleles for the heat shock (A) and cold shock (B) response.Because each gene's expression was normalized to the zero time point, the first component represents deviations from this time point.Strains are represented by the first letter of the genus and species name: Sc, Sp, Sk, and Su.

FIG. 5 .
FIG.5.-Five genes exhibit consistent allele differences between thermophilic and cryophilic species.Lines show the mean and whiskers show the point range.