Phenotypics and transcriptomic analysis of the 59A×S288c lineage
We phenotyped 44 segregants obtained from a cross between the laboratory strain S288c and the wine yeast derivative 59A (as described in Methods). Fermentations were performed in a synthetic medium simulating a grape must (SM425) and containing para-amino benzoate (PABA) to counteract the effect of the ABZ1
]. Our analysis were performed in more stringent conditions than that of Ambroset et al.
]: the amounts of ergosterol and oleic acid used were half those generally used (the final concentrations were 7.5 mg/l for ergosterol and 2.5 μl/l for oleic acid) and fermentations were performed at 24°C. We estimated the kinetic properties regarding the lag time and the fermentation rate at 2 stages of the fermentation: Rmax (maximal fermentation rate) and R70 (fermentation rate at 70% of fermentation). The fermentation kinetics of parental strains and of several segregants are shown in Figure . Population size and metabolite quantities were determined at the end of the fermentation.
Figure 1 Fermentation profiles for parental strains and for segregants. (A) Fermentation profile for the two parental strains and for 59A×S288c: 59A in blue, S288c in green, 59A×S288c in black. (B) Fermentation profile for 8 segregants. Sugar was (more ...)
We first checked that phenotypes affected by ABZ1
allele, such as Rmax, were corrected by the addition of PABA and independent of the allelic form (data not shown). Most of the phenotypes had a high heritability (80% to 97%), indicating that genetic variations had a major impact on overall variations. The dry weight was the only phenotype displaying a low heritability (>50%) and was disregarded. R70 values followed a bimodal distribution, suggesting that the phenotype is mainly controlled by one locus. The other phenotypes, such as Rmax and the amounts of metabolites followed continuous distributions, indicating a probable polygenic control (Additional file 1
Transcriptome profiles were obtained at 70% of fermentation progress (66 g/l CO2 released), corresponding to late stationary phase, 20 to 40 hours after the end of the growth phase (depending on the segregants). At this point in the fermentation, the yeast is subject to nutrient starvation and ethanol stress (8% v/v). We compared RNA abundance between segregants by Agilent mono-color labeling and hybridization on oligonuclotides microarrays as described in Methods. We identified 1610 genes differentially expressed between 59A and S288c with a log2 of fold-change (LogFC) of more than 0.7 in either direction (adjPv<0.01).
We assessed the correlation between the fermentation rate (R70) and gene expression and found large sets of genes positively correlated (347 genes involved in cell trafficking and nitrogen reutilization; Funspec [33
]) and negatively correlated (351 genes involved in oxidative stress and respiration; Funspec).
Clustering analysis revealed a partial disomy of chromosome 16 affecting genes expression
A clustering analysis of the whole gene expression dataset revealed significant associations (Pearson correlation coefficient ≥ 0,84) for 887 genes (Additional file 2
). Most of the gene clusters obtained were enriched in functional classes, such as mating, mitochondrial translation, ribosomal proteins, ergosterol synthesis and sulfur assimilation consistent with the known coordination of expression of these genes. Some clusters also corresponded to genes correlated with phenotypes (R70, Rmax). Other clusters included physically associated genes (ASP cluster, regions A B and C [13
], telomeric genes).
We identified a cluster of 37 genes displaying a higher expression in 13 segregants. These genes were physically linked and all were located in the first 373 kbp at the start of chromosome 16 (Additional file 2
, cluster “chr16 left-arm genes”). A careful mining of the CGH data (from the Affymetrix chip signals) revealed that these 13 segregants carried a 373 kbp duplicated region on the left arm of chromosome 16 (Figure A). No other chromosomal aberration was observed. This partial disomy resulted from the translocation of an arm of chromosome 16 onto the chromosome 8 originating from the wine yeast strain, in addition to the presence of a standard chromosome 16 (Figure B). This 8–16 translocation is known to be responsible for higher level of sulfite resistance and was found in many wine yeasts [10
]. We examined the impact of this partial disomy on transcriptome profiles by comparing the expression data for the 13 disomic segregants with those for the other 31 segregants. As expected, the genes of the duplicated region were more strongly expressed in segregants carrying the disomic region (155 of the 191 genes in the area concerned, adjPv<0.05), consistent with the doubling of gene copy number. We also observed differences in expression for 1305 other genes (adjPv<0.01) not located in this area. The upregulated genes include a high proportion of genes involved in methionine biosynthesis (17 genes, Pv=10-9
]) whereas the down-regulated genes included a high proportion of genes involve in the oxidative stress response (20 genes, Pv=10-6
-FunSpec) and aerobic respiration (37 genes, Pv=10-14
-FunSpec). The pattern of functional enrichment of the deregulated genes suggests that disomic strains are less affected by stresses of alcoholic fermentation conditions. Consistent with this hypothesis, we observed that the partial disomy had a major impact on fermentation kinetics (Figure C). Partial disomy was associated with significantly higher fermentation rates, a shorter duration of fermentation (t
) and an absence of sluggish fermentation profiles. Levels of acetic acid production were also lower for the disomic segregants (data not shown, t
-test Pv<0.01), supporting our hypothesis of lower levels of stress in these clones.
Figure 2 The 8–16 translocation and impact on fermentation. (A) CGH profile for the segregant 7a (in red) which has the partial disomy, in comparison to S288C CGH profile (in black). (B) Parental and daughter genotype: Crossing 59A which has the translocation (more ...)
Phenotypic and expression QTL linkages
We generated a new high-density genetic map (1.81 markers/10 kbp Additional file 3
) from the Affymetrix S98 microarrays genotyping data for the 44 segregants. Linkage analysis was performed by the interval mapping method, for QTLs and eQTLs.
Three phenotypic QTLs were identified with a false discovery rate (FDR) threshold of 0.05: one QTL explaining flocculation on chromosome 1 encompassing FLO1, another explaining clumpiness and encompassing AMN1 and a third, explaining fermentation rate (R70), on the left arm of the chromosome 16 (140 kbp), overlapping with the duplicated region and clearly triggered by the partial disomy. Limiting the population to the 31 segregants with no partial duplication of chromosome 16 eliminated the detection of the R70 QTL on chromosome 16. It also led to the detection of three QTLs for metabolites levels, mapping to chromosome 2 (208 kbp) for pyruvic acid and chromosome 14 (460 kbp) for succinic acid and glycerol.
Genome wide linkage analysis of gene expression led to the detection of 1063 eQTLs as summarize in Table . Dividing the FDR by two only slightly reduce the number of eQTLs, demonstrating the robustness of the data. A high amount (56%) of local-eQTLs (or cis-eQTLs, see Methods) were detected with higher average LOD than trans-eQTLs (Table ). Some of the genes differentially expressed between the parents displaying cis-eQTLs corresponded to genes absent from one of the parental strains. As expected, genes from the three new regions acquired by horizontal transfer in EC1118 [13
] and conserved in 59A formed three clusters (Additional file 2
), each controlled in cis-regulations. We did not consider these genes, nor the 111 genes missing from EC1118 [13
] for subsequent analyses, reducing the number of eQTLs to 967.
We found that 113 of the 1460 genes affected by the partial disomy displayed eQTLs, and 48 of these eQTLs mapped to the deleted or duplicated regions on chromosomes 8 and 16. Given that the partial disomy could affect the eQTL data, we adapted the analysis (analysis 1) to take this effect into account. Two new eQTL linkage analyses were performed as described in Figure A: (a) a search with a population limited to 31 non disomic segregants (analysis 2) (b) a filtering of the data for the 13 disomic segregants, eliminating the markers in the duplicated region and ignoring the 1460 genes affected by the partial disomy (analysis 3).
Figure 3 Analyses to by-pass the effect of partial disomy on transcriptome. (A) Principle of the analyses: Analysis 2 was the safest way to eliminate all the partial disomy effect and would have been applicable with a larger population. However, using a 31 segregants (more ...)
These two new data sets were compared with the first one on a Venn diagram Figure B and in Additional file 4
. Most of the eQTLs (547) were common to all the three data sets (group 123 on Figure B) and corresponded to eQTLs with high LOD-score (cis-eQTLs). In analysis 2, restricted to the 31 non disomic segregants, the detection power was lower: at lod4, the FDR was 0.13. Data set 3 overlapped strongly with data set 1, with 888 of 967 eQTLs common to these two data sets. This suggests that the partial disomy slightly decreased our ability to detect true eQTLs, but had only a limited effect on the robustness of the data. The eQTL analysis was performed for the union of the three data sets (1465 eQTLs), taking the “Venn” group for each eQTL into account. We, therefore, concluded that we could overcome the effect of the partial disomy.
Variations in expression networks are triggered by few loci resulting in the formation of hotspots
The locations of the phenotypic QTLs and eQTLs are summarized in Figure and Table . Many eQTLs overlapped in several hotspots and most were common to the three analyses. However, the disomy was responsible for two hotspots on chromosomes 8 and 16 (hotspot 6 and hotspot 10, in red Figure B). Furthermore, 113 of the eQTLs discovered in the two new analyses are overlapping in two hotspots on chromosome 2, at 208 kbp and 283 kbp (hotspots 1 and 2 in blue/green Figure A-B). Due to the genetic similarities between EC1118 and RM11-1a, some phenotypic QTLs and eQTLs are identical to those previously identified for the cross between BY4716 and RM11-1a (BxR population: hotspot 3, 4, 8 and 9, [25
Figure 4 QTL and eQTL localization. (A) Localization of the 1465 eQTLs in function of the corresponding transcript position. Colors correspond to the group in the Venn diagram Figure B. Spot in the first diagonal correspond to cis-eQTLs. (B) (more ...)
However, we also found some linkages not previously identified. Hotspot 5 on chromosome 4 was found to be enriched in genes involved in thiamine biosynthesis pathway and is probably controlled by variations of the THI3
regulator gene located in the hotspot area. Hotspot 7 on chromosome 10 brings together 18 genes with various functions (mitochondria, glucose sensing: SNF3
) and GRR1
seams to be the most relevant candidate gene for the control of this hotspot. It is indeed involved in glucose repression and amino acid sensing [36
We also found variations in several networks involving fewer genes but potentially involved in the wine-making process. The high-affinity sulfate permease gene, SUL2
, has a cis-eQTL (LOD 4.5) overlapping with trans-eQTLs for three other genes of the sulfur pathway MET2
has two SNPs in its promoter (but none in the consensus sequences for transcription factors binding sites, Yeastract [37
]) and five non-synonymous SNPs in its coding sequence. Changes in sulfate transport efficiency are clearly responsible for the higher expression level of the three other genes involved in sulfur amino-acid metabolism (Figure ).
Figure 5 Control of genes involved in sulfuric amino-acid synthesis by SUL2 expression. (A) Correlation between MET2, MET5 or MET32 expression and SUL2 expression. Orange spots correspond to the segregants with SUL2 from S288c. Blue spots correspond to the segregants (more ...)
Similarly, we observed a perturbation of the nitrogen uptake system that might have a significant impact on wine yeast performance [38
]. The membrane peptide transporter gene, PTR2
, and the arginine/alanine aminopeptidase gene, AAP1
, both have eQTLs mapping to CUP9
(chromosome 16) which encodes the known transcriptional repressor of PTR2
has two non-synonymous SNPs in its coding sequence and the 59A allele clearly increases the expression of PTR2
and, probably, also that of AAP1
(Figure ). Byrd et al.
] described the activation of PTR2
by ubiquitin-dependent Cup9p degradation in the presence of peptide [40
]. However, in the 59A strain, PTR2
has a frame-shift in its sequence and is, therefore, probably non-functional [41
]. The involvement of CUP9
in the control of AAP1
is consistent with a broader role of this regulator in peptide metabolism. The synthetic must used here did not contain peptides, but these molecules may have been generated by endogenous nitrogen metabolism. The mutation in 59A may act by modulating the repressive activity or degradation of Cup9p.
Figure 6 Correlation between PTR2 and AAP1 expression.CUP9 allelic origin is likely to control both genes. Red open circles correspond to disomic segregants, CUP9 being localized in the left arm of chromosome 16. Orange spots correspond the segregants with CUP9 (more ...)
Finally, a flocculation QTL was found on chromosome 1, mapping to FLO1
as previously described [24
] (Figure C). Flocculation phenotype was found linked to FLO1
expression level (Additional file 6
). Surprisingly, multiple QTL searches on the basis of flocculation phenotype did yield a locus containing FLO8
. However, a second QTL was found on chromosome 1, in a 45 kbp region (Additional file 6
). This region contains the gene OAF1
involved in fatty-acid and peroxisome biogenesis. It is possible that OAF1
mutations triggered variation of cell surface hydrophobicity that could impact flocculation [42
Variations in detoxification mechanisms are triggered by mutations in transcription factor genes
In this linkage analysis, we detected changes in the detoxification network. Several membrane transporters involved in drug export displayed eQTLs (Table ). Some of these genes were under local regulation (cis-eQTL): SVS1, required for resistance to vanadate, the two polyamine exporter, TPO1 and TPO2, and the ion transporters ALR2 and SSU1. The cis-eQTL of SSU1 is due to the 8–16 translocation. Other detoxification system genes (e.g. SNG1, PDR12, QDR2) displayed distant control. Several eQTLs associated variations in the expression of these transporters with variations of zinc-finger transcription factor genes: HAP1, PDR8, YRR1 and WAR1.
Genes involved in detoxification having an eQTL
gene, which displays a cis-eQTL, has three SNPs in its promoter area and three non synonymous SNPs in its coding sequence. The S288c allele of TPO1
is more strongly expressed than the wine yeast allele. As TPO1
is involved in octanoic acid resistance [3
], we investigated the effect of variations of TPO1
expression on octanoic acid resistance. We measured population growth in the synthetic must medium (pH3.3) supplemented with octanoic acid (0.2 mM). Octanoic acid resistance dependent principally on the TPO1
-test Pv<0.005, Figure ) but was also weakly correlated with TPO1
expression level (data not shown). Our data therefore suggest that the form of TPO1
encoded by the S288c allele is more effective, conferring a higher octanoic acid resistance.
Figure 7 TPO1 allele controls is own expression and octanoic acid resistance. (A) Box plot showing TPO1 allele effect over TPO1 expression level in the population. (B) Box plot comparing segregants resistance to octanoic acid by TPO1 allele (t-test Pv<0.005). (more ...) PDR11
all had an eQTL in hotspot 8 mapping to HAP1
(Table ). Another gene, PDR8
, encoding a pleiotropic drug resistance transcription factor lies in the same locus, 25
kbp away from HAP1
, and may control some of these genes (Additional file 5). Consistent with this hypothesis, Steyer et al.
], recently showed that high nerolidol production were associated with the S288c allele of PDR8
. They demonstrated that QDR2
regulation by PDR8
accounted the variation of nerolidol release into the medium. We assessed the impact of PDR8
allelic variation on gene expression by performing an allelic switch for PDR8
and constructed 59A strain expressing the PDR8
allele from S288c (59A PDR8
-S288c). We compared the transcriptome of this strain with that of wild type 59A in fermentation conditions (the same conditions as for the global analysis). We detected variation of expression for very few genes (adjPv<0.01, Table ). QDR2
was the gene most strongly upregulated by allelic replacement, with a LogFC of 2.24, consistent with the results of Steyer et al.
]. Among the other targets of PDR8
, we only found YLR179c
controlled by the allelic form of PDR8
. Both QDR2
had an eQTL mapping to PDR8
. In addition to a slight down regulation of PDR8
itself, we identified three other genes (YLR149c
) as upregulated in the constructed strain and potentially playing a role in detoxification. These genes have not previously been described as PDR8
-controlled. We did not detect the targets of PDR8
reported by Hikkel et al.
]. This was probably due to the specific conditions of alcoholic stress and nutrient starvation used in our study. Furthermore, we performed an allelic switch with a smaller impact than the over-expression of PDR8
used in other studies. Finally, only QDR2
were controlled by PDR8
in hotspot 8, suggesting that expressions of the other genes are triggered by HAP1
Gene expression modified by PDR8 allelic switch
The gene paralog of PDR8
, displays a cis-eQTL on chromosome 15 and has four SNPs in its promoter region and six non-synonymous SNPs in its coding sequence. SNG1
, one of the targets of YRR1
responsible for nitrosoguanidine resistance [44
], has an eQTL mapping to the position of YRR1
. Surprisingly, we found no linkages of this locus to other YRR1
target genes. Moreover, the inverse patterns of parental behavior and locus segregation for SNG1
expression indicates that other loci are involved in SNG1
control (Figure ).
Figure 8 YRR1 control over SNG1 expression. Orange spots correspond the segregants with YRR1 from S288c, blue spots correspond the segregants with YRR1 from 59A. Open circles correspond to disomic segregants, There is a permutation between parent and the allelic (more ...)
gene, encoding a plasma membrane ABC transporter responsible for organic acid efflux, had an eQTL mapping to its known transcription factor WAR1
(chromosome 13), an activator of weak acid stress response. The WAR1
allele of 59A has five non synonymous SNPs in its coding sequence and was associated with higher levels of PDR12
expression level (Figure A). In addition, WAR1
was itself one of the genes for which expression was affected by HAP1
variation (hotspot 8). The level of WAR1
expression had no significant impact on PDR12
transcript abundance. Instead PDR12
expression levels were influenced by the type of WAR1
allele (Figure B). PDR12
was the only known WAR1-
dependent gene with an eQTL mapping at WAR1
(even with a LODscore below the threshold) [45
]. As PDR12
has been implicated in sorbic acid resistance [46
], we investigated whether variation of the expression of this gene modulated resistance to this acid. We assessed the growth of segregants in the synthetic medium (pH3.3) supplemented with sorbic acid (0.5 mM, see Methods). Sorbic acid resistance appeared to be partly controlled by the WAR1
allele (Figure C, t
-test Pv<0.005) consistent with the effect of WAR1
allele on PDR12
expression. We then performed an allelic switch for WAR1
and constructed a 59A strain expressing the WAR1
allele from S288c (59A WAR1
-S288c). The growth profiles of the strains were monitored in the synthetic must medium supplemented with various amount of sorbic acid. The lag phase was longer for 59A WAR1
-S288c than for 59A wild-type (Figure D). This phenotypic difference was subtler than that observed with a knockout strain, but these results nevertheless confirm that sorbic acid resistance is linked to PDR12
expression and to WAR1
genetic variation. Variations of detoxification networks may thus be linked to variations of several transcription factor genes. Two of these linkages were validated by allelic switch experiments. Moreover, we demonstrated that the levels of resistance to octanoic acid and sorbic acid were clearly modulated by these genetic variations.
Figure 9 Control of sorbic acid resistance by WAR1 allele. (A) Box plot showing WAR1 allele effect over PDR12 expression level in the population. Values for parental strains are indicated. (B) Confrontation of WAR1 expression level and PTR12 expression level show (more ...)
As polymorphism may result from specific adaptation, we tried to infer the phylogeny of War1p and Yrr1p from nucleotidic sequence available for different strains in SGD. Interestingly; the two allelic versions (wine and laboratory) of these two genes diverge from their S. paradoxus
orthologue, and are clearly divergent but apparently in a different manner. For WAR1
, S288C and 59A alleles are located on different branches, whereas for YRR1
, the phylogeny indicates two clusters one containing palm wine and Asian alleles and the second cluster containing 59A wine, rum and bread alleles, as well as S288C alleles. However in this case S288C and 59A are also located on different branches of this cluster (Additional file 7
, MEGA 5 [47
]). In order to detect if these two genes have evolved neutrally, we performed a McDonald Kreitman test calculating the neutrality index (NI) [48
]. This test was applied to a set of 15 WAR1
alleles from strains isolated from various substrates and NI was 2.28 (Pv = 0.02), indicating a significant excess of non-neutral mutations. As the cluster of Asian alleles might inflate this test, we removed the sequence and performed the test again, which indicate a similar result (NI = 2.56 and Pv = 0.03). We performed the similar analysis on YRR1
alleles. We obtained NI = 2.51 (Pv = 0.001) for the whole set of alleles and NI = 3.93 (Pv = 0.002) after removal of the Asian and palm wine alleles. This suggests that WAR1
do not evolve neutrally: these two genes are either subject to the accumulation of slightly deleterious mutations that are eliminated by negative selection during speciation, or alternatively they present substantial diversity that might be associated with balanced selection resulting from specific adaptation of the strain to wine fermentation conditions.