Genetic Mapping of MAPK-Mediated Complex Traits Across

Summary

Wild yeast strains differ in phenotypes that are controlled by highly conserved signaling pathways. Yet it remains unknown how naturally occurring genetic variants influence signaling pathways in yeast. We have developed an approach to facilitate the mapping of genetic variants that underlie these phenotypic differences in a set of wild strain. Our mapping approach requires minimal strain engineering and enables the rapid isolation of mapping populations from any strain background. In particular, we have mapped genetic variants in twelve highly diverse yeast strains. Further, we demonstrate how the mapping results from these twelve strains can be used jointly to narrow the number of genetic variants identified to a set of putative causal variants. We identify genetic variants in genes with various roles in cell signaling. Our results illustrate the interplay of different signaling pathways and which signaling genes are more likely to contain variants of large phenotypic impact.

Introduction

Cellular survival is dependent on the ability to sense and respond to changing environmental conditions. Mitogen-activated protein kinase (MAPK) signaling cascades are ubiquitous in eukaryotic organisms and enable them to react to extracellular stimuli [1]. MAPK cascades are composed of three sequentially acting kinases [2], which upon sensing an extracellular stimulus trigger a cellular response by activating transcription factors and other regulatory proteins [2], [3]. The yeast Saccharomyces cerevisiae has long been used as a eukaryotic model system to elucidate the principles of these signaling pathways, as many of the core signaling proteins are conserved from yeast to human [3].

In yeast, MAPK pathways facilitate response to environmental cues such as the presence of mating pheromones and stresses such as cell wall damage and high osmolarity [4]. Adaptation to high osmolarity is conducted by the HOG (high osmolarity glycerol) MAPK pathway. The human homolog of the MAPK Hog1, p38α, not only mediates the response to hyperosmolarity as well [4], but also plays key roles in inflammation and cancer [5], [6]. Insults to the cell wall are sensed by the cell-wall integrity (CWI) pathway, which is anchored by the MAPK Slt2, homolog of the mammalian MAPK7 [7]. While these pathways have been scrutinized in great detail, it remains largely unknown how they are affected by genetic variation among different yeast isolates. Response and resistance to MAPK-activating stress conditions are highly variable among different yeast isolates [8], [9], [10]. While genes of the core MAPK cascades are highly conserved across species, upstream regulatory components, such as stress sensors, and downstream targets, such as transcription factors, exhibit high levels of divergence [11], [12]. Yet, it is unknown how genetic differences in such elements of the MAPK pathways contribute to the phenotypic differences between isolates of the same species. We sought to examine how sequence variation among yeast isolates results in quantitative differences in MAPK-dependent phenotypes.

Yeast is an ideal model system for the study of quantitative traits. Quantitative trait locus (QTL) mapping studies in S. cerevisiae have characterized the genetic architecture of complex traits, including global gene expression and resistance to small molecules [13], [14], [15], [16], as well as shed light on the sources of “missing heritability” [17]. We have previously developed a bulk-segregant-analysis approach (BSA) [18], [19], named X-QTL, which utilizes extremely large mapping populations for increased mapping power [20]. Such mapping populations are subjected to selection to enrich for alleles that contribute to a trait of interest. We have recently used this method to determine the genetic underpinnings of protein abundance variation [21]. The majority of mapping studies to date have interrogated a pair of strains at a time, with only some recent studies expanding to crosses among four parent strains [22], [23]. Isolation of individual recombinant haploids, as used in linkage mapping studies, is very labor intensive, and current X-QTL protocols require extensive strain engineering. Consequently, mapping studies have queried a fraction of the S. cerevisiae species-wide genetic and phenotypic variation and have left much of the genetic architecture of quantitative traits across the species unexplored.

Here, we introduce a novel X-QTL approach and use it to map variants that alter MAPK-dependent traits in a panel of diverse yeast isolates. Our method takes advantage of plasmid-borne fluorescent markers to leverage the precision and high throughput of fluorescence-activated cell sorting (FACS) to rapidly generate large mapping populations of recombinant haploids. We charted the genetic architectures of salt (sodium chloride) and caffeine tolerance in a set of strains that capture much of the species-wide genetic diversity. Salt tolerance is mediated by the HOG MAPK pathway, while exposure of yeast to caffeine results in the activation of the CWI MAPK pathway [7]. In addition to these two pathways, we leverage the ability of our new mapping approach to isolate mapping populations of either mating type to examine the effect of the MAPK mating pathway on growth. We have identified QTL that contain genes with various roles in MAPK signaling and other signaling pathways. Our results illustrate how genetic variation within S. cerevisiae shapes MAPK-dependent traits. Finally, we illustrate how a round-robin cross design can be leveraged to identify causal genetic variants.

Results

Fluorescence-based isolation of mapping populations

Our original X-QTL approach relied on the expression of the auxotrophic marker HIS3 under the control of the MATa-specific STE2 promoter for the isolation of MATa recombinant haploids [24]. However, most isolates of S. cerevisiae are generally prototrophs, lacking the auxotrophic mutations required for the original X-QTL approach. In addition, auxotrophic mutations can have large phenotypic effects [14]. For instance, the HIS3 deletion required under the original X-QTL approach results in sensitivity to various salts, which is not desirable when investigating natural variation in salt resistance [25]. Furthermore, the use of the STE2 promoter restricted X-QTL mapping populations to segregants of the MATa mating type. Here, we designed a novel approach that relies on mating-type specific expression of plasmid-borne fluorescent reporters, enabling the isolation of haploids of either mating type without the need to integrate the reporters into the genome (Fig. 1A). We placed the fluorescent markers mCitrine and mCherry under the control of the STE2 promoter and its MATα-specific equivalent, the STE3 promoter, respectively [24]. We combined these two fluorescent markers in a single construct and created a set of plasmids with drug resistances that permit facile introduction into any strain background (see Methods for details).

As expected, the presence of this construct resulted in green fluorescence in MATa strains and red fluorescence in MATα strains (Fig. 1B & C). Diploids carrying the marker construct were non-fluorescent, as neither promoter is active in diploid cells [24], [26]. After sporulation and germination, haploid progeny express the fluorescent marker corresponding to their mating type. We tested this aspect of the reporter construct by introducing it into the diploid hybrid between the well-characterized BY and RM strains and generating a pool of BY/RM recombinant haploids (see Methods for details). Flow cytometry revealed that this pool encompassed three cell populations: non-fluorescent cells and, in roughly equal proportions, green and red cells (Fig. 1D). We used fluorescence-activated cell sorting to isolate populations of green and red cells and then sequenced them in bulk to determine the BY and RM allele frequencies in each population. The BY parent strain contributes the MATa allele of the MAT locus, located on chromosome 3, while the MATα allele originates from the RM parent. As expected, the MATa allele was highly enriched in the green fluorescent recombinant haploids, while the MATα allele was enriched in the red population (Fig. 1E). This result demonstrates that the fluorescent markers can be used to isolate populations of MATa, as well as MATα, haploids.

In addition to the expected large allele frequency differences at the MAT locus, we identified several allele frequency skews that were highly concordant between the two populations (Fig. 1E). These skews correspond to previously identified growth QTL [17], [20]. For example, the skew on chromosome 12 is located at the gene HAP1; the BY allele of this gene results in partial loss of function and poorer growth compared to the RM allele, leading to the observed reduction in frequency of the BY allele. The agreement between the allele frequency distributions in the two mating type pools indicates that the isolated cell populations can be used to map QTL with excellent reproducibility in both mating types.

Phenotyping of wild strains and design of the round-robin cross

We set out to use our approach to dissect phenotypes mediated by different MAPK pathways. We used two stress conditions to query two different MAPK-signaling networks: high salt concentrations and caffeine. High salt conditions lead to activation of the HOG signaling cascade, while caffeine triggers the cell-wall integrity pathway.

We determined sodium chloride and caffeine tolerances for a collection of 65 S. cerevisiae isolates [27] and selected twelve diverse strains that spanned a range of resistances for our mapping studies (Fig. 2A, S1 Table). These twelve strains vary at approximately 80% of the common SNPs (allele frequency ≥0.05) segregating within the larger strain collection. We designed a round-robin crossing scheme in which each of the twelve strains was crossed to two other strains, for a total of twelve crosses (Fig. 2B). While the round-robin crossing scheme is the most efficient way to intercross a set of strains while maintaining equal contributions of each strain [28], it has only been used to map traits in a small number of organisms [29], [30]. Most mapping studies are based on crosses between resistant and sensitive strains, which implicitly assumes that resistance alleles are shared between resistant strains. By contrast, we randomly assigned the order of strains within the round-robin design without regard to strain phenotypes. As such, the crossing scheme included crosses between parent strains with different, as well as with similar, resistance phenotypes.

To assess the genetic complexity of salt and caffeine resistance, we isolated up to twenty individual segregants from each cross. We measured growth of these segregants and the parent strains under high salt and caffeine stress conditions (S1 Fig.). Three of the 24 phenotype distributions of the cross progeny were approximately bimodal, suggesting one large-effect determinant of resistance in these crosses. Other crosses exhibited phenotype distributions consistent with directional genetics and transgressive segregation [31]. In two crosses, the mean progeny phenotype was significantly lower than the midparent phenotype (S2 Table), indicative of interactions between the underlying causal loci [31], [32]. These findings illustrate the genetic complexity of the resistance phenotypes among wild isolates.

Mapping MAPK-dependent traits within the round-robin crosses

Next, we established a pipeline to map QTL between any pair of strains. Sequences of parent strains were used to generate lists of SNPs that segregated within a cross and could serve as genetic markers for QTL mapping. Mapping populations of both mating types were isolated and grown under permissive as well as selective conditions (YPD, YPD containing 0.5 M sodium chloride, 1 M sodium chloride, 15 mM caffeine or 20 mM caffeine). The resulting populations were bulk-sequenced to assess allele frequencies at SNP sites segregating in a particular cross. We used the MULTIPOOL software to calculate LOD scores testing whether allele frequencies at a given position in the genome in each selected population differed from those in the corresponding control population (S2 Fig.) [21], [33].

We used replicate BY/RM mapping experiments to determine the robustness and reproducibility of our mapping approach (see Materials & Methods for details). Allele frequency distributions of replicate experiments were extremely similar: comparisons between replicates did not exceed a LOD of 2.62 (S3A–S3B Fig.). At a LOD score threshold of five, QTL detection was highly reproducible, with 90.1% of QTL reproducibly detected in replicate experiments (S3C–S3D Fig., S3 Table). Reproducibility was also high between selections based on mapping populations of opposite mating types (88.3%).

We applied our mapping strategy to identify QTL that influence MAPK-dependent resistance traits in the twelve Round-robin crosses. We used the MATa and MATα selections for each cross as replicate mapping experiments (see S4A–S4B Fig. for combined LOD plots for each cross). We used a stringent criterion that a QTL must exceed a LOD threshold of five in both selections, and identified 155 QTL across the twelve crosses and the different conditions. Of QTL identified in only one of the two matched experiments, the majority had LOD scores close to the LOD threshold (median LOD 6.7 for unreplicated QTL vs. 12.3 for replicated QTL, Fig. 4A). Indeed, for 93.5% of unreplicated QTL the direction of the underlying allele frequency skews was the same in the matched MATa and MATα experiments (permutation based p<0.001). This suggests that many of these QTL stem from small effect variants that narrowly escaped detection in one of the two experiments.

Next, we asked how QTL are distributed and shared among the Round-robin crosses (Fig. 3 A & B, S4 Table). The salt QTL fell into 37 unique genomic regions, while the caffeine QTL comprised 23. There were 13 QTL regions that recurred between the conditions suggesting they maybe due to pleiotropic variants that affect both traits. We found an average of 6.2 salt resistance QTL (range 1 to 13) and 3.7 caffeine resistance QTL (range 0 to 7) per cross. Interestingly, the number of QTL detected per cross was independent of the genetic distance of the parent strains involved (divergence & salt QTL: Pearson correlation 0.154, divergence & caffeine QTL: 0.087). However, when the QTL of the two conditions were considered together per cross, this resulted in a weak correlation with the extent of parental divergence (0.257). We tested if the number of QTL per cross was a function of the phenotypic difference of the parent strains, since a previous study found that a higher number of QTL could be identified between yeast strains of similar phenotypes [34]. We found that this was the case for the salt resistance QTL (phenotypic difference & salt QTL: Pearson correlation −0.517), but not for caffeine resistance QTL (phenotypic difference & caffeine QTL: 0.395). While increasing the phenotypic similarity of yeast parent strains through the experimental fixation of large effect QTL can permit the identification of additional small effect QTL [35], [36], this trend is condition-specific for crosses between different wild isolates.

Because each genetic variant is introduced into the Round-robin scheme at least twice, each QTL is expected to be identified at least twice assuming 100% power and no genetic interactions. Contrary to this expectation, previous studies involving sets of interconnected crosses have found that the majority of QTL are ‘context-dependent’, that is their identification was dependent on a specific cross and not just on the parent strains involved [23], [34], [37]. Indeed, 33 of the 60 QTL regions we identified across the two conditions were found in only one cross (Fig. 4B). The overall lower LOD scores of these ‘context dependent’ QTL as a class partially explained their lack of detection in additional crosses (Fig. 4C; Pearson correlation 0.417, p = 2.6×10⁻⁶). In addition, interactions between the QTL and the genetic background likely contribute to context dependency [37].

Functional connections of QTL to signaling

We asked if the QTL mapping had resulted in any functional enrichment in the underlying genes. While genes within caffeine resistance QTL were only marginally enriched for the GO-term ‘Response to Stimulus’ (p = 0.077), genes found under salt resistance QTL were significantly enriched for ‘Cellular Response to Stimulus’ (p = 0.006) and ‘Kinase Activity’ (p = 0.002) (S5 Table). This suggests that our mapping experiments point towards genetic variants that impact cellular signaling pathways.

Several of the identified QTL tower among the genetic architectures of the two resistance traits (Fig. 3A & B, S4 Table). Tolerance of high salt concentrations is dominated by variation at the ENA locus on chromosome 4, which we detected in eleven out of twelve crosses. The ENA genes encode sodium efflux pumps that, in the presence of high salt, are transcriptionally activated by the HOG pathway [38]. European S. cerevisiae isolates carry an ENA variant that resulted from an introgression from the yeast species S. paradoxus and exists in varying copy numbers, referred to as ENA1 through ENA5. In contrast, non-European strains carry a single copy of the ancestral ENA6 gene [39]. Both the introgressed variant and its copy number have been shown to be associated with high salt tolerance [10], [40]. We determined the copy numbers of the two ENA variants within our parent strains and found that strains with high salt resistance indeed carried multiple copies of the introgressed ENA gene (S6 Table). Interestingly, the identified ENA QTL reflected differences beyond the two divergent ENA variants and copy number. We identified ENA QTL in crosses between parents with equal copy numbers of the introgressed variant (Cross 9), as well as between non-European isolates carrying alleles of the ENA6 gene that differed from each other by three missense mutations (Cross 12). Our findings corroborate the importance of the ENA locus in shaping the genetic architecture of salt resistance [10], [40] and further illustrate the multi-allelic nature of this locus.

The hallmarks of S. cerevisiae caffeine resistance architecture are the QTL containing TOR1 and TOR2 on chromosomes 10 and 11, respectively (Fig. 3B). In contrast to other eukaryotes, yeast carries two copies of TOR. As in higher organisms, TORC complexes respond to nutrient signals to regulate a plethora of cellular processes through two functionally distinct complexes. In yeast, Tor1p is specific to TORC1 complexes, while Tor2p is able to participate in both TORC1 and TORC2 complexes [41]. Caffeine inhibits Tor1p function [42], leading to the activation of the CWI pathway [43]. We identified a TOR1 QTL in eight of the twelve crosses, illustrating the large impact of TOR1 on the genetic architecture of caffeine resistance. TOR2 QTL were found in five crosses, and their identification likely speaks to Tor2p's ability to functionally compensate for Tor1p in TORC1 complexes. However, the QTL we detected as encompassing TOR genes were not mutually exclusive. In the four cases where we found both TOR1 and TOR2 QTL within the same cross, resistance alleles originated from the same or different parent strains an equal number of times. Furthermore, TOR negatively regulates transcription factors that are activated by the HOG pathway in response to stress (reviewed in [11]). Indeed, we uncovered salt resistance QTL containing TOR2 in three crosses (S4 Table).

In addition to these QTL shared among many of the Round-robin crosses, several other QTL had large effects in a smaller number of crosses and contain genes with well-established roles in cellular signaling. A caffeine resistance QTL shared between crosses 2 and 3 encompasses MSS4, which encodes the essential Phosphatidylinositol-4-phosphate 5-kinase required for activation of the CWI MAPK pathway [44]. Two additional large effect QTL highlight the interconnectivity of the cellular signaling network. These QTL on chromosomes 10 and 15 were pleiotropic and strongly influenced growth in the presence of both high salt and caffeine. The QTL on chromosome 10 contains CYR1, the adenylate cyclase essential for cAMP production and as such cAMP-PKA signaling [45]. Similar to TOR, the cAMP-PKA pathway responds to nutrient signals to modulate a myriad of cellular processes, including the negative regulation of stress response transcription factors (reviewed in [11]). The QTL on chromosome 15, on the other hand, includes a gene that represents a functional counterpart to this aspect of cAMP-PKA signaling: WHI2 encodes a phosphatase required for activation of the general stress response [46].

Impact of strain history and source on QTL identified

While trait variation among S. cerevisiae isolates is great in comparison to related fungal species [10], the evidence for adaptive alleles is scant [9], [47]. Rather than being strongly shaped by their source environment, yeast phenotypes have been proposed to be a result of population history [10]. The ENA locus illustrates this pattern. The introgression of the S. paradoxus ENA variant is widespread in European isolates regardless of their source environment. Our parent strains included European clinical and vineyard isolates, yet these ecological niches had no bearing on the ENA QTL we identified in crosses involving these strains. We wondered if alleles, denoted by the strain they derived from, tracked with source environment for any of the other QTL. We examined QTL detected in non-overlapping crosses (S7 Table), which excludes QTL likely due to single instances of loss- or gain-of-function mutations, and found two QTL where the grouping of alleles was unlikely to be a result of shared population history. Caffeine resistance alleles that gave rise to TOR2 QTL belonged to wine and clinical strains from a diverse range of geographic sources, while the corresponding sensitive alleles stemmed from natural isolates, such as oak strains. A salt resistance QTL on chromosome 15, containing YGK3, which encodes a stress response controlling kinase, exhibited a similar pattern: a European and an American clinical strain contributed resistance alleles, while diverse natural isolates supplied sensitivity alleles. Further studies are needed to clarify whether these two cases are coincidental or whether they represent cases of environmental adaptation.

Identification of MAT-dependent growth QTL

Yeast mating, akin to the stress responses examined above, is regulated by a MAPK signaling pathway. In addition to providing us with independent selection experiments, our ability to isolate mapping populations of either mating type allowed us to search for variants that influence growth through a genetic interaction with the mating type locus. Fitness differences between the two mating types have been reported [48], [49], and we have shown that the MAT locus influences the differential expression of numerous genes [13]. In order to detect MAT-dependent growth QTL, we directly compared the allele frequencies resulting from corresponding MATa and MATα selection experiments.

For the round-robin crosses, we first compared the permissive YPD selections for each mating type. Among the twelve crosses, we identified six significant QTL (S8 Table). One of these QTL appeared in two non-neighboring crosses and encompassed the mating pathway regulator GPA1 located on chromosome 8. Polymorphisms within Gpa1, which negatively regulates the mating pathway, can impair its function and hence result in costly activation of the pathway in the absence of mating pheromone [50]. Importantly, the GPA1 allele selected for in one of our QTL experiments (Fig. 5) contains a known fitness-associated variant [50], [51]. We previously uncovered a genetic interaction between the MAT locus and GPA1 that affects the expression of other mating pathway genes [13]. While this interaction did not result in a growth QTL in the BY/RM cross, our identification of MAT-dependent GPA1 QTL suggests that this interaction can have an effect on growth in other genetic backgrounds. In addition to the GPA1 QTL, we identified a QTL containing SST2, which encodes the GTPase-activating protein for Gpa1 [52].

Next, we asked if mating type had any influence on the allele frequency distributions arising under the two MAPK-inducing stress conditions. We found 19 instances of MAT-dependent stress resistance QTL (S8 Table). While it is difficult to determine whether these QTL represent an instance of MAT-dependence or failure to replicate a weak QTL, 14 of these QTL were identified in independent experiments using different dosages of the same stress condition or in both stress conditions. Among the latter group, we found that the GPA1 and SST2 mating type dependent QTL persisted in the presence of either stress condition. The mating and HOG MAPK pathways share the MAPKKK Ste11, and several studies have investigated how signaling specificity is achieved despite this shared signaling component [53], [54]. In particular, faithful signaling along either pathway could be maintained by either mutual inhibition [53] or pathway insulation [54]. Our results suggest that the baseline fitness cost of the mating pathway [50] remains unaltered in the face of stressors that activate other MAPK cascades, which is consistent with the insulation of the mating MAPK from other MAPK signaling cascades [54].

Utilizing the Round-robin design to identify causative variants

The Round-robin design enables a strategy to identify candidate causative variants that is not available using pair-wise crosses. Specifically, the segregation of individual variants in relation to the presence or absence of QTL in multi-parental crosses can be leveraged to reduce the search space of candidate causative variants [55], [56].

Many sequence variants occur in several of the Round-robin crosses and, assuming additive effects independent of genetic background, should have similar consequences in these different crosses. We generated de novo assemblies of the parent strain genomes [57] and cataloged the non-synonymous coding variants in each QTL. We determined whether a given variant segregated within each cross by examining if the respective parent strains carried different alleles of the variant. For each QTL, we then assessed to what extent the segregation of a particular variant in each cross was associated with the detection of the QTL in the crosses (S4 Table). We used the maximum LOD scores within a particular QTL interval for each cross as a quantitative measure of QTL detection. This association analysis resulted in a 17-fold reduction in the number of variants to be considered as potentially causative for each QTL. We confirmed causality of individual sequence variants in two QTL to illustrate how this approach can aid in the identification of causal variants.

We first focused on the QTL on chromosome 15 that contained the stress response regulator WHI2 [46] and strongly influenced the resistance to both salt and caffeine in crosses 7 and 8 (Fig. 6A). Both crosses shared CLIB219 as a parent strain, and in both sets of experiments, the CLIB219 allele was strongly selected against, indicating that it leads to stress sensitivity. As this QTL did not appear in this form in other crosses besides the two involving CLIB219, the causative variant likely was private to CLIB219. Four variants out of 64 variants in the region were private to CLIB219 (Fig. 6B). One of the variants was a frame-shift mutation in WHI2. Allele replacements of the CLIB219 allele of WHI2 resulted in reduced growth in the presence of high salt and caffeine and thus confirmed that it was indeed the causal variant (Fig. 6C).

Next, we turned to a QTL with a substantially more complex allelic pattern. The TOR1-containing caffeine resistance QTL on chromosome 8 was identified in eight crosses. Specific alleles were selected for or against depending on the particular cross, which was a clear indication that multiple alleles were present at this QTL (Fig. 6D). While no single variant perfectly recapitulated the pattern of LOD scores (Fig. 6E), two SNPs in TOR1 and NTA1 scored the highest. Closer examination of the segregation pattern of these SNPs revealed that they were common to two caffeine sensitive strains, M22 and YJM981. Allele replacements of the TOR1^3875A variant in two strain backgrounds confirmed that it contributes to caffeine sensitivity (Fig. 6F). These results demonstrate how the Round-robin crossing scheme can be leveraged to focus the list of candidate causative variants.

Discussion

Here, we have demonstrated the ability of a novel QTL mapping approach to query a large proportion of S. cerevisiae genetic diversity. We employed this approach to examine the genetic architectures of three MAPK-dependent traits: resistance to high salt or caffeine and growth differences due to mating type. The relationship of the mapped genetic variants to MAPK signaling pathways ranged from upstream modulators to core regulators and downstream targets.

Our modified X-QTL approach makes it possible to rapidly generate large mapping populations without extensive strain engineering. Quantitative trait mapping studies in yeast have uncovered the genetic architectures of traits central to biology, such as gene expression and protein-level variation [14], [21], [58], as well as traits of medical and industrial importance [15], [59]. Yet, these studies have only scratched the surface of the phenotypic variation of S. cerevisiae strains [10], [22], a space that is sure to expand with the continued discovery of additional yeast isolates [60], [61]. Our method promises to be of great utility in future yeast QTL mapping studies. Here, we have chiefly used our approach to isolate large pools of segregants for BSA-based QTL mapping, but we also illustrate its utility in isolating individual segregants that can be used in traditional linkage-based QTL mapping studies. In contrast to previous methods [20], [62], our method can be employed in nearly any strain background with minimal strain engineering and enables positive selection for haploids of either mating type. Furthermore, the fluorescent mating type markers are amenable to combination with fluorescent readouts of cell physiology or gene expression [20], [21]. Such reporter combinations will make it feasible to isolate mapping populations and select for phenotypes of interest simultaneously.

We have genetically dissected stress resistances mediated by two MAPK pathways, as well as interrogated the genetic contribution to the baseline cost of the mating pathway. We studied twelve diverse yeast strains that represent a large proportion of the genetic and phenotypic diversity of S. cerevisiae. Thus the QTL identified capture much of the species-wide genetic variation that influences the phenotypic consequences of these MAPK pathways. We identified QTL that encompass genes connected to MAPK signaling at multiple levels (S4 Table). Several QTL contained genes whose protein products are responsible for the first step of the MAPK signaling cascades: osmosensors (MSB2 & SHO1) and CWI sensors (WSC3 & MTL1) responsible for detecting the respective perturbations. Moving along the signaling cascades, we found several key regulators such as MSS4 and GPA1. Furthermore, one of the caffeine tolerance QTL contained SLT2, which encodes the MAPK at the center of the CWI response. Arriving at the output of the signaling pathways, we identified several genes transcriptionally activated in response to stress. In addition to the ENA genes, we found GPG1, which is activated by the HOG pathway to increase glycerol synthesis to counteract changes in osmolarity [63]. Finally, although the resistance phenotypes of the parent strains were not correlated, thirteen QTL occurred under both conditions. The overlap between QTL illustrates the crosstalk between the two stress responsive pathways [64], [65].

In addition to genes directly associated with MAPK pathways, we found genes that highlight the interconnectivity of the cellular signaling network. Variants in CYR1 and TOR2, core members of the cAMP-PKA and TOR pathways respectively, were tied to both stress responses. These pathways sense the internal availability of nutrients and repress the general stress response [11]. The TOR pathway in particular regulates a myriad of functions that mediate cell growth [66]. In contrast to higher eukaryotes, yeast carries two copies of the TOR gene. Our study links numerous alleles of TOR1 and TOR2 to differential growth under stress conditions. The two Tor proteins are redundant in regards to their participation in TORC1 complexes, but only Tor2p can partake in TORC2 complexes [41]. Despite this difference, the two TOR genes are diverging at the same rate as the majority of duplicated gene pairs [67]. However, these studies were based on individual strains. Our QTL results point towards functional variation in both Tor isoforms and raise questions about how they jointly regulate cell growth in different yeast isolates. Are the roles of the two proteins in TORC1 complexes truly equivalent, or are the TOR alleles indicative of functional tradeoffs? Similarly, it bears determination whether stress sensitive alleles of the TOR2 QTL, which were found in natural rather than ‘domesticated’ isolates, are beneficial under conditions approximating their source environment.

The identification of WHI2 further echoes the interplay between internal nutrient and external stress sensing pathways. While WHI2 is required for the activation of stress response pathways [46], its loss leads to prolonged growth under nutrient limited conditions [68]. Interestingly, WHI2 loss-of-function alleles were found to be driver mutations in experimental evolution experiments [69], [70]. Variants that abrogate Whi2 function were also discovered as recurrent secondary site mutations in yeast strains with specific gene deletions [71]. As such, whi2 mutants not only illustrate how laboratory evolution experiments can recapitulate the biology of wild isolates, but also suggest that the concomitant loss of stress responsiveness is adaptive under permissive growth conditions [69], [70]. Indeed, the WHI2 allele responsible for stress sensitivity was advantageous for growth under the permissive YPD condition.

The QTL results provide insights into the types of genes that carry variants shaping MAPK-associated phenotypes of different yeast isolates. The core elements of the HOG and CWI signaling cascades (HOG: SSK1, SSK2, STE11, PBS2, and HOG1; CWI: RHO1, PKC1, BCK1, MKK1, MKK2, and SLT2) are highly conserved across species. Analysis of the strain sequences revealed that despite their species-level conservation, all of these core proteins, with the exception of Rho1p, contain instances of non-synonymous coding variants. Yet, we only identified QTL at SLT2, indicating that the coding variants in the other core genes may be of little functional consequence.

In contrast to the central cascade, proteins responsible for sensing perturbations in osmolarity or cell wall integrity diverge rapidly between yeast and related fungal species [12], [72]. Divergent sensor proteins in other species are thought to enable the response to different habitats. For instance, other fungi use the HOG pathway to respond to a broader range of stresses than S. cerevisiae (reviewed in [11]). Consequently, we wondered if variation in sensor genes underlies the observed stress resistance variation among different S. cerevisiae isolates. We found stress sensors at two QTL for both salt and caffeine resistance, but these QTL had low LOD scores indicative of small effect variants. This is likely explained by the large degree of redundancy among the sensor proteins [73], [74]. Osmosensors are encoded by a set of five genes, and CWI sensors are encoded by six genes. While sensor redundancy allows for species-level fine-tuning of the stress responses, any single change is unlikely to have large phenotypic consequences. The same holds true for the large number of transcriptional targets of the MAPK pathways, with the exception of special effectors such as the introgressed ENA variants and GPG1 [75]. In contrast, many genes associated with large effects, such as CYR1, GPA1, MSS4 and TOR1/2, have functions essential to the MAPK and other signaling pathways. Together, these findings illustrate how genetic variants shape the phenotypic space of MAPK stress resistance traits within S. cerevisiae. The global picture that emerges from our results is that genes thought to underlie phenotypic differences between fungal species due to their high levels of divergence contribute little to trait variation among S. cerevisiae isolates.

The QTL we identified fall into several different classes in regards to their appearance among the twelve round-robin crosses. 33 of the 60 grouped QTL were ‘context dependent’, meaning they were detected in a single cross. The preponderance of ‘context dependent’ QTL was also notable in studies involving smaller panels of pairwise crosses [23], [34]. Context-dependency has been attributed to interactions between QTL and the genetic background. In particular, large effect variants can mask the presence of small effect variants [35], [36]. For salt resistance, the negative correlation between parent strain phenotype differences and the number of QTL identified was largely driven by the large effect of the ENA locus. For example, Cross 7 between two sensitive parent strains lacking the ENA introgression permitted the detection of a large number of small effect QTL. Of the twelve QTL detected in exactly two crosses, seven appeared in crosses involving the same parent strain, and as such can be classified as ‘strain dependent’ [34]. This set of QTL likely results from single alleles private to the shared parent strain [23]. Beyond this simple case, deducing the number of underlying alleles and their frequency becomes difficult. For instance, of the eleven QTL appearing in three crosses, nine were found in a pair of neighboring crosses plus one additional cross. Several models could explain this pattern, and we attempted to use the segregation pattern of the sequence variants within each QTL to determine which alleles best explained the observed QTL pattern. While this approach resulted in a significant reduction in the number of variants to consider, more than twelve crosses will be necessary to distinguish between the possible allele patterns with certainty. Going forward, our fluorescent markers allow facile creation of segregant pools from many additional crosses, as well as large mapping panels of individual segregants [17]. In conjunction with new genome engineering technologies, these approaches will enable even more complete genetic dissection of MAPK signaling and other complex traits.

Materials and Methods

Plasmid construction

Plasmids with drug markers were constructed in two steps. We first generated pRS416 and pRS426 Gateway expression vectors suitable for one and three fragment Gateway cloning. As described previously, the pRS vectors were digested with Sma1 and the Gateway cassettes introduced by ligation [76]. The drug markers were then introduced from a set of plasmids kindly provided by Mikko Taipale and Susan Lindquist (Whitehead Institute, Cambridge, MA). These Gateway drug markers plasmids were digested with BglII and XhoI in order to create fragments with the MX drug resistance markers alone. PCR was used to generate fragments of the pRS416 and pRS426 Gateway vector backbones with Xho1 and BglII cut sites and without the URA3 marker. PCR products were gel extracted and digest with Xho1 and BglII. The vector backbone with the two different Gatewat cassettes and the three different drug resistance markers where then combined by ligation.

The STE2 and STE3 promoters were cloned upstream of the genes encoding the fluorophores mCitrine and mCherry using overlap extension PCR. Promoter fragments based on Tong et al. [24] were PCR amplified from genomic DNA, while fluorophores were amplified from plasmids kindly provided by Gregory Lang (Lewis-Sigler Institute, Princeton University, Princeton, NJ). The resulting fusions, as well as the CYC1 terminator, were cloned into different Gateway entry vectors to facilitate their final combination in the three fragment Gateway plasmids.

Plasmids created for this study are listed in S9 Table, while primers used are listed in S10 Table.

Strain construction

The fluorescent mating type markers were tested in BY strains kindly provided by the Botstein lab (BY 12045 MATa Δura3, BY 12046 MATα Δura3). The BY/RM diploid was published previously (yLK1993 [17]). Other strains were previously described in Schacherer et al. [27]. To facilitate the generation of diploids for the round-robin cross, the KanMX cassette used to delete HO was replaced by the HphMX cassette in MATx strains [77]. MATa and MATx parent strains were mated in YPD and subsequently streaked to YPD+G418+Hygromycin to select for diploids. Putative diploids were tested for their ability to sporulate and to determine their ploidy [78]. The p41Nat plasmid carrying the fluorescent mating type markers was introduced by standard yeast transformation [79].

We used the NatMX cassette from pUG74 [80] to delete URA3 gene in the haploid round-robin parent strains in order to make them amenable to allele replacement protocols. Attempts to make allele replacements according to Erdeniz et al. were unsuccessful [81]. We modified the approach described by Stuckey et al. to make allele replacements [82]. For the WHI2 allele replacements, the entire WHI2 ORF was first deleted using the pCORE construct [83] and then replaced by transformation with PCR-amplified WHI2 ORFs plus approximately 250 bases up- and downstream from the different strains. For the TOR1 allele replacements, we first integrated K. lactis URA3 right at position 3875 of TOR1. We then removed the K. lactis URA3 integration according to Storici et al. using complementary oligos spanning the integration site and carrying the two different TOR1 alleles (TOR1 3875G and TOR1 3875A). Primers used in allele replacement experiments are listed in S10 Table.

Strain phenotyping

We phenotyped a set of strains previously described by Schacherer et al. [27], as well as sets of random segregants stemming from the round-robin crosses. The strains were grown overnight at 30°C in 96-well plates containing YPD media according to standard yeast culture protocols [78]. The strains were then pinned in quadruplicate onto agar plates using a Singer RoToR HDA. Strains were phenotyped on YPD plates (the permissive control condition), as well as YPD plates containing 0.5 M sodium chloride, 1 M sodium chloride, 15 mM caffeine or 20 mM caffeine. Both sodium chloride and caffeine were autoclaved along with the other ingredients needed to make standard YPD plates. After strain pinning, plates were grown at 30°C in large plastic containers to avoid excessive evaporation. YPD plates were scanned after two days of growth, while stress condition plates were scanned once growth reached a level equivalent to that observed on the YPD plates [17]. Images of the plates were processed and analyzed to extract strain growth information as previously described [17]. Growth under the stress conditions was adjusted for growth under the permissive stress condition by calculating the ratio of the two growth measurements. We used a modified version of the epistasis test described by Lynch and Walsh [32] as detailed in Brem and Kruglyak [31].

Spore isolation

We modified a published random spore prep protocol [78] in order to enrich for spores prior to FACS. Diploids were grown in 3 ml YPD overnight at 30°C. Cultures were spun down and washed with water. Cell pellets were resuspended in 1 ml sporulation media containing 25 mg/L ClonNAT or 200 mg/ml G418 to maintain the selection for fluorescent mating type marker plasmids. Half of the cell mixture was then added to 2.5 ml drug containing sporulation media. Sporulation cultures were incubated for four to six days at room temperature depending on the sporulation efficiency of a particular cross. 250 µl of the sporulation cultures was spun down for one minute in an eppendorf tube. After removal of supernatant the cell pellets were resuspended in 100 µl water. 17 µl of the cell suspensions was then added to new eppendorf tubes containing 3 µl of β-glucoronidase. After a brief vortex tubes were incubated in an Eppendorf Thermomixer at 30°C for two hours while shaking at 900 rpm. Digestion and opening of ascii after this two hour incubation was confirmed by light microscopy. Digestions were quenched by the addition of 100 µl water. In order to break apart digested ascii, tubes were vortexed for two minutes and then spun down for one minute. The supernatant was removed, another 100 µl was added and the tubes were vortexed for another two minutes. As hydrophobic spores stick to the sides of the eppendorf tubes all liquid is removed after the second vortex. Tubes were washed three times by addition of 1 ml water and inverting the tubes fully three times. After the washes, 1 ml 0.01% NP-40 (igepal) was added added to the tubes. Tubes were the sonicated using a tip sonicator for 1 min at power level 2.

In order to isolate individual spores or to assess the amount of spores isolated, dilutions of the spore prep (for most crosses 25 µl spore mix in 250 µl H2O was used) on the plates containing 25 mg/L ClonNAT or 200 mg/ml G418. Individual haploids were picked after two days of growth at 30°C using a fluorescence-equipped dissection microscope. For cell sorting 250 µl of the samples was plated on four plates containing the appropriate drug to select for the marker plasmid.

Selection and FACS

Spore preps were harvested after 2 days of growth at 30°C. Cells were washed off plates using 5 ml water and collected in 50 ml conical tubes. Cells were sonicated for 1 min at power level 2 and then diluted to an OD of 0.4 for FACS in 1xPBS. We used a BD FACSVantage SE to sort 500,000 green and red fluorescent cells each into vials containing 2 ml YPD with 100 µg/ml ampicilin to prevent bacterial contamination during the cell sorting. The sorted cells were added to 10 ml YPD containing 100 µg/ml ampicilin in glass tubes and cultured for 6 hours at 30°C. Cultures were spun down for 5 minutes at 3,000 rpm and after removal of supernatant cells were resuspended in 950 µl water. 200 µl of cells suspension were then plated on 5 different conditions: YPD, YPD with 0.5 M NaCl, YPD with 1 M NaCl, YPD with 15 mM caffeine, and YPD with 20 mM caffeine. Plates were grown for 2 days at 30°C. Plates were made by adding the appropriate amount of NaCl or caffeine to the standard YPD recipe prior to autoclaving.

DNA isolation and library prep

We assessed to what extent segregant populations grew under the selective conditions and proceeded with conditions that resulted in a reduced yet adequate amount of growth compared to the permissive YPD condition [20]. Cells were harvested by washing plates with 5 ml water and collected in 15 ml conical tubes. Cells were spun down and stored at −80°C after removal of supernatant. DNA was isolated from the harvested cells according to Lee et al. [84]. Cell pellets of approximately 50 µl were used for each DNA prep. DNA was isolated from parent strains in the same manner except that parent strains were grown overnight in 3 ml YPD at 30°C.

Sequencing libraries were generated using Epicentre Nextera DNA Sample Prep kits or Illumina Nextera Sample Prepartation kits as described previously [17]. Residual salts and other contaminants can cause Nextera library preps to fail. We cleaned the input DNA using the Zymo Research Genomic DNA Clean & Concentrator kits if the initial library prep failed. Sample libraries were indexed using primers designed according to Meyer and Kircher [85] and sequenced in groups of 25 per lane on an Illumina HiSeq sequencer.

Allele frequency measurements and QTL mapping

We created a set of genetic markers for each cross by aligning parent strain sequences to the sacCer3 reference genome and we then used custom scripts written in the programming language R to identify SNPs that differed between any set of parent strains. We then determined allele frequencies at these variant sites. We subtracted the YPD allele frequencies at each variant position from those obtained under the different growth conditions to determine the allele frequency skews that are due to a specific condition rather than permissive growth. To determine allele frequency differences due to mating type, allele frequencies from MATα samples were subtracted from those of corresponding MATa samples.

We converted allele frequency skews into LOD scores using the MULTIPOOL method [33] with the following parameters: -t, -n 1000, -m contrast, -c 2200, -r 100. MULTIPOOL outputs were used to generate genome-wide LOD score plots [21]. LOD plots were generated in R. All the subsequent analyses were conducted using the statistical programming language R. We called QTL at positions that exceeded a LOD threshold of 5 and QTL regions using LOD drop intervals of 2, as described previously [21]. In order to combine QTL across replicate experiments we determined which QTL had overlapping mapping intervals. We then used the mean of replicated QTL intervals to identify QTL that overlapped between different round-robin crosses. We combined QTL that were detected at different concentrations of the same conditions (0.5 M NaCl & 1 M NaCl, 15 mM Caffeine & 20 mM Caffeine) to provide a conservative estimate of the QTL LOD scores and regions.

We used a set of replicate BY/RM mapping experiments to determine the false discovery rate and reproducibility of our mapping approach. We generated allele frequency measurements for four independent experiments using two independent transformants of the BY/RM diploid with the mating type marker plasmid. We isolated mapping populations of either mating type and used the same selection conditions as for the round-robin crosses (YPD, +0.5 M NaCl, +1 M NaCl, +15 mM Caffeine, +20 mM Caffeine). The YPD 20 mM Caffeine condition resulted in very low growth of BY/RM segregants and we did not proceed further with this condition. Given the two mating type mapping populations and four independent experiments, we generated eight selections for each condition. We compared the allele frequencies resulting from technical replicate experiments stemming from one transformant (16 comparisons) and found no allele frequency fluctuations that resulted in LOD scores higher than 2.62. Analysis of biological replicates, i.e. comparing experiments based on different transformants, resulted in six peaks with LOD scores higher than our threshold of LOD 5. Yet, five of these peaks were actually the same region identified multiple times in different replicate pairs suggesting that they were due to a genetic difference between the two BY/RM transformants and not due to random allele fluctuations. The remaining peak tagged the mating type locus on chromosome 3 pointing towards a difference in the mating type selection rather than a random allele frequency fluctuation.

In order to asses the reproducibility of the mapping approach we determined condition-specific QTL for individual experiments and then asked to what extent QTL where found in pairs of replicates. To keep this assessment of reproducibility consistent with the round-robin mapping experiment, we did not ask whether QTL were identified multiple times across all replicates. Such a measure of reproducibility would be inflated in comparison to experiments with only two replicates, such as the experiments with the round-robin crosses. Rather we tested if QTL were reproducibly identified in designated pairs of experiments, either technical replicates performed with the same strain and the same condition or the experiments of opposite mating types of one replicate. At a LOD threshold of 5, 90.1% of QTL were reproducibly identified between technical replicates of the same mapping experiment. Reproducibility was 88.3% using mapping populations of opposite mating type. For the experiments with technical replicates, the directions of the allele effects of the QTL identified in only one experiment were consistent among replicates in 100% of the cases. This suggests that these QTL only narrowly missed the LOD threshold in replicate experiments. Indeed, their LOD scores were close to the threshold of LOD 5 (median LOD 5.9 vs. 9.84 for replicated QTL). For the round-robin experiments across all conditions, 74.4% of QTL (311 of 418 individually identified QTL, with one instance of one peak overlapping two in the second dataset, resulting in 155 jointly identified QTL) were identified in both of the mapping experiments. Since QTL that narrowly miss the LOD threshold are classified as not being reproducible, we asked to what extent allele frequency directions, in the control case BY versus RM, agreed over the identified region. We determined the mean allele frequencies for the region underlying a QTL and compared them across the pairs of replicate experiments. For the BY/RM experiments, the sign of the allele frequencies agreed for all eleven QTL that had missed replication. In the case of the round-robin crosses, allele frequency directions agreed for 93.5% of non-replicated QTL. This result was robust to permutation of the allele frequencies (p<0.001).

GO Analyses

We used VEP [86] to identify genes with non-synonymous coding sequence variants within the combined QTL intervals. GO analyses were performed using the tools provided by the Saccharomyces Genome Database [87].

Testing for association between variants and LOD scores

We used Cortex [57] to generate reference-assisted genome assemblies of the round-robin parent strains. We then used VCFtools [88] to subset the VCF file generated by Cortex to the individual regions of the QTL. We used VEP [86] to identify non-synonymous coding sequence variants within QTL intervals. We then generated a table for each QTL in the statistical programming language R that contained the replicate-averaged maximum LOD scores for each cross for the given QTL region and notations how each variant segregated among the crosses. For each variant, we then calculated the squared correlation between its segregation pattern and the pattern of LOD scores. We determined the distribution of for all the variants of a given region and used a p value of 0.01 to call significant variants. The coefficients of determination were visualized along QTL intervals using the UCSC Genome Browser [89].

Supporting Information

Attachment 1

Attachment 2

Attachment 3

Attachment 4

Attachment 5

Attachment 6

Attachment 7

Attachment 8

Attachment 9

Attachment 10

Attachment 11

Attachment 12

Attachment 13

Attachment 14