text:
Simultaneous quantification of mRNA and protein in single cells reveals post-transcriptional effects of genetic variation
Trans-acting DNA variants may specifically affect mRNA or protein levels of genes located throughout the genome. However, prior work compared trans-acting loci mapped in separate studies, many of which had limited statistical power. Here, we developed a CRISPR-based system for simultaneous quantification of mRNA and protein of a given gene via dual fluorescent reporters in single, live cells of the yeast Saccharomyces cerevisiae. In large populations of recombinant cells from a cross between two genetically divergent strains, we mapped 86 trans-acting loci affecting the expression of ten genes. Less than 20% of these loci had concordant effects on mRNA and protein of the same gene. Most loci influenced protein but not mRNA of a given gene. One locus harbored a premature stop variant in the YAK1 kinase gene that had specific effects on protein or mRNA of dozens of genes. These results demonstrate complex, post-transcriptional genetic effects on gene expression.
Phenotypic variation in genetically complex traits is shaped by DNA variants throughout the genome. Many variants that influence complex traits alter gene expression (Albert and Kruglyak, 2015; Maurano et al., 2012). Understanding this regulatory genetic variation is critical for our knowledge of genetic disease, evolutionary change, and for applications of biology in agriculture and industry.
Regulatory genetic variation is typically studied by measuring mRNA abundance in genetically diverse populations to map regions of the genome that contain variants that cause variation in mRNA levels (expression quantitative trait loci, eQTLs). Many genes are influenced by cis-acting eQTLs that are caused by DNA variants in cis-regulatory elements. The proximity of cis-eQTLs to the genes they affect has aided their identification in numerous species (GTEx Consortium, 2017; Albert et al., 2018; Brem et al., 2002; Cheung et al., 2005; Clément-Ziza et al., 2014; Hasin-Brumshtein et al., 2014; Higgins et al., 2018; Ka et al., 2013; Kita et al., 2017; Morley et al., 2004; West et al., 2007). However, most genetic variation in gene expression arises from the combined effects of many trans-acting variants throughout the genome (Albert et al., 2018; Grundberg et al., 2012; Wright et al., 2014). Trans-acting variants affect the activity or abundance of diffusible factors that in turn alter the expression of other genes that can be located anywhere in the genome. Therefore, the identification of trans-eQTLs requires a large number of tests for statistical association between genotypes and mRNA levels. This, combined with the fact that trans-eQTLs have smaller effects than cis-eQTLs, has hampered the identification of individual trans-eQTLs in human association studies. By contrast, experimental crosses in organisms such as yeast (Albert et al., 2018; Brem et al., 2002; Brion et al., 2020; Clément-Ziza et al., 2014; Emerson et al., 2010; Thompson and Cubillos, 2017), plants (Fu et al., 2013; West et al., 2007; Zhang et al., 2011), worms (Snoek et al., 2017; Viñuela et al., 2010), flies (Coolon et al., 2014; Everett et al., 2020; Huang et al., 2015; McManus et al., 2010), and mouse (Gerrits et al., 2009; Hasin-Brumshtein et al., 2016 ) have quantified the global contribution of cis and trans-acting variation and identified thousands of trans-eQTLs.
Post-transcriptional mechanisms can shape protein levels independently of mRNA abundance (Buccitelli and Selbach, 2020; Csárdi et al., 2015; McCarthy, 1998), raising important questions about how DNA variation shapes protein levels compared to mRNA abundance. It is unclear whether most eQTLs, which by definition affect mRNA, also affect the protein abundance of the same genes, such that they also result in a protein-QTL (pQTL) (Damerval et al., 1994). It is also unclear to what extent there are protein-specific pQTLs that influence protein abundance without affecting the mRNA of the same gene. Recent studies have begun to address these questions (Battle et al., 2015; Cenik et al., 2015; Chick et al., 2016; Foss et al., 2011; Foss et al., 2007; Ghazalpour et al., 2011; Großbach et al., 2019; Mirauta et al., 2020; Picotti et al., 2013; Sun et al., 2018; Wu et al., 2013; Yao et al., 2018) and reported a wide range of observations.
Across several species, half or more of cis-eQTLs were reported to also affect protein abundance, suggesting high levels of agreement between mRNA and protein (Albert et al., 2018; Albert et al., 2014b; Hause et al., 2014; Skelly et al., 2013). However, other work found many cis-eQTLs whose effects on protein were smaller than those on mRNA, suggesting the existence of mechanisms that buffer protein levels against variation in mRNA abundance (Battle et al., 2015; Ghazalpour et al., 2011; Mirauta et al., 2020). Trans-acting QTLs were reported to show even greater differences between mRNA and protein. In mouse crosses and some yeast crosses, (Mirauta et al., 2020) trans-eQTLs and trans-pQTLs mapped for the same genes showed very little overlap, suggesting that while many trans-eQTLs are buffered at the protein level, there are also many trans-pQTLs that affect proteins via mechanisms that do not affect mRNA (Chick et al., 2016; Foss et al., 2011; Foss et al., 2007; Ghazalpour et al., 2011; Williams et al., 2016). While more recent work in yeast showed somewhat higher agreement between trans-eQTLs and trans-pQTLs (Albert et al., 2018; Großbach et al., 2019), yeast crosses also revealed instances of 'discordant' trans-acting loci that affect both mRNA and protein of the same gene, but in opposite directions (Albert et al., 2018; Albert et al., 2014b; Foss et al., 2011; Großbach et al., 2019). In plants, many trans-eQTLs appear to be buffered at the protein level, with few protein-specific pQTLs (Fu et al., 2009). In human populations, about half of the identified pQTLs could be accounted for by underlying eQTLs, although the exact fraction varied between studies (Emilsson et al., 2018; Hause et al., 2014; Sun et al., 2018; Yao et al., 2018). In yeast and humans, genetic effects on translation as measured by ribosome profiling were similar to those on mRNA, suggesting that translation cannot account for protein-specific pQTLs (Albert et al., 2014a; Battle et al., 2015; Cenik et al., 2015). Instead, protein-specific effects may arise from protein degradation, especially for proteins that are part of protein complexes (Chick et al., 2016; Großbach et al., 2019; He et al., 2020). Together, this literature suggests that genetic variants can independently affect the different layers of gene expression regulation.
However, there are substantial caveats to this conclusion. First, small sample sizes of dozens to a few hundred individuals were used in most studies to date. Most eQTLs and pQTLs change the expression of the gene they affect by 10% or less (Albert et al., 2018). Detecting such small-effect loci requires high statistical power resulting from the analysis of large numbers of individuals (Albert et al., 2018; Albert et al., 2014b; Bloom et al., 2013; Ehrenreich et al., 2010; Visscher et al., 2017). With limited sample size, real loci can be missed. For example, a locus that truly affects both mRNA and protein but happens to pass the detection threshold in mRNA but not in protein data could be incorrectly inferred to be buffered at the protein level. Thus, low statistical power can inflate discrepancies between eQTLs and pQTLs.
Second, experimental differences between studies may also inflate discrepancies between eQTLs and pQTLs. In nearly all cases, eQTLs and pQTLs were compared between experiments conducted in different laboratories, sometimes using different experimental designs, and, in the case of some human studies, comparing eQTLs and pQTLs mapped in different tissues. These comparisons are likely to be confounded by environmental influences, which can drastically alter the effects of regulatory variants. For example, stimulation of human immune cells can reveal new cis-eQTLs and trans-eQTLs (Fairfax et al., 2014; Lee et al., 2014). In yeast, many trans-eQTLs are specific to a specific culture medium (Lewis et al., 2014; Smith and Kruglyak, 2008), and some trans-eQTLs switch the sign of their effect in different media (Smith and Kruglyak, 2008). Even when experiments aim to recreate the same conditions, subtle differences (for example, in the precise stage of cell growth, medium composition, temperature, or handling of the isolated mRNA or proteins) could influence measures of gene expression (Gallego Romero et al., 2014; Gibson et al., 2007; Hayeshi et al., 2008), potentially inflating the discrepancy between eQTLs and pQTLs.
To date, no study of regulatory variation has measured mRNA and protein with high statistical power and in the same samples. This challenge is not easy to overcome, because combined application of RNA-sequencing (RNA-seq) and mass spectrometry in thousands of matched samples remains prohibitively expensive. While recent studies in yeast (Großbach et al., 2019) mice (Chick et al., 2016; Williams et al., 2016), and human iPSC lines (Mirauta et al., 2020) used RNA-seq and mass spectrometry on the same cultures or tissue samples, their sample sizes of less than 200 individuals limited QTL detection. As a result, the role of genetic variation on post-transcriptional processes remains unclear, especially for trans-acting variation.
Here, we addressed this challenge by developing a system for quantifying mRNA and protein from the same gene simultaneously, in the same, live, single yeast cells using two fluorescent reporters. We reasoned that such an approach would equalize all environmental confounders and most of the technical biases that could obscure the relationship between eQTLs and pQTLs. Our assay is sensitive enough to be used in fluorescence-activated cell sorting (FACS), permitting the use of a well-powered bulk segregant QTL mapping approach (Michelmore et al., 1991) that leverages millions of genetically different yeast cells to yield high statistical power (Albert et al., 2014b; Parts et al., 2014). Genetic mapping across ten genes revealed 86 trans-acting loci. Although our approach was explicitly designed to maximize the chance of detecting similar genetic effects on mRNA and protein, the majority of the identified loci affected only mRNA or protein for a given gene, or had discordant effects on mRNA and protein. These results demonstrate considerable differences in the genetic basis of variation in mRNA vs protein abundance.
A reporter system for quantifying mRNA and protein in single, live cells
We designed a dual reporter system for the simultaneous quantification of mRNA production and protein abundance of a given gene in single, live cells. In this system, protein abundance is measured via a fluorescent green fluorescent protein (GFP) tag fused to the C-terminus of the given protein of interest (Huh et al., 2003). To measure mRNA, we reasoned that a clustered regularly interspaced short palindromic repeats (CRISPR) guide RNA (gRNA) (Doudna and Charpentier, 2014) produced in equal molarity with the mRNA of interest would be able to drive proportional expression of a reporter gene via CRISPR activation (Gilbert et al., 2014; Konermann et al., 2015). To implement this idea, we created a gRNA tag located in the 3'UTR of the gene, downstream of the sequence encoding GFP (Figure 1A). After transcription of the mRNA along with this tag, the gRNA is released from the mRNA by two flanking self-cleaving ribozymes (Hammerhead, Hh; and Hepatitis Delta Virus, HDV) (Gao and Zhao, 2014). Because gRNA cleavage separates the mRNA from its poly-adenylated (polyA) tail, we added a synthetic polyA tail between the GFP tag and the Hh ribozyme (Gao and Zhao, 2014). Once released, the gRNA directs a catalytically deactivated CRISPR associated enzyme (dCas9) fused to a VP64 activation domain (dCas9-VP64) to drive the expression of a red fluorescent mCherry gene integrated in the genome (Farzadfard et al., 2013). After gRNA release, stability and half-life of the mRNA no longer affects gRNA abundance, such that mCherry expression primarily reports on mRNA production.
The reporter system is implemented as two cassettes (Figure 1A). The 'GFP-gRNA tag' cassette is added at the 3' end of the gene of interest. A second cassette, which we call the 'CRISPR reporter', comprises the remaining components: dCAS9-VP64 and the mCherry gene under the control of an inactive CYC1 promoter fragment. This promoter contains one recognition sequence that, when targeted by the gRNA and dCas9-VP64, drives mCherry expression (Farzadfard et al., 2013). The two cassettes are stored on two plasmids that can be used to easily construct strains for quantification of mRNA and protein of any gene of interest (Figure 1—figure supplement 1).
We tested the reporter system in diploid BY strains tagged at two genes with different expression levels: the highly expressed TDH3, and GPD1, which has an average expression level compared to other genes in the genome. Both genes gave green and red fluorescent signals in a plate-reader (Figure 1B). A strain carrying the CRISPR reporter and TDH3 tagged with GFP but no gRNA produced no mCherry fluorescence, demonstrating that the gRNA is required for driving mCherry expression (Figure 1B & C). Quantitative real-time reverse-transcription PCR (qPCR) confirmed expression of the gRNA and the mRNA (Figure 1D). Absence of qPCR signal from primers that spanned the ribozyme cut sites in cDNA confirmed that the ribozymes cleaved the mRNA (Figure 1D and Figure 1—figure supplement 2).
mCherry fluorescence provides a quantitative readout of mRNA production
To characterize the quantitative response of our reporter system to a range of gene expression levels, we used the synthetic Z3EV system, which allows quantitative regulation of transcription via the concentration of estradiol in the culture medium (McIsaac et al., 2013). We cloned the Z3EV promoter upstream of a GFP-gRNA sequence (Figure 2A) in a strain that also contained the CRISPR reporter and grew this strain in a range of estradiol concentrations. Along with the expected increase in green fluorescence (McIsaac et al., 2013), red fluorescence increased as a monotonic function of estradiol concentration (Figure 2B). Similar results were observed in the RM11-1a strain, which has a different genetic background than BY (Figure 2—figure supplement 1). Thus, mCherry provides a quantitative readout of the expression of the tagged gene.
While green fluorescence continued to increase throughout the tested estradiol range, red fluorescence ceased to increase at concentrations of more than 4 nM estradiol (Figure 2B). qPCR quantification of the gRNA showed that mCherry fluorescence followed gRNA abundance (Figure 2C), confirming that the mCherry reporter gene is quantitatively regulated by gRNA abundance. gRNA abundance was linearly related to GFP mRNA and GFP fluorescence at lower doses of estradiol but stopped increasing at higher doses (Figure 2D & E). This suggests that mCherry production is limited by gRNA availability at high expression levels. Increasing the concentration of dCas9 proteins or binding sites for the gRNA had no effect on the mCherry expression plateau (Figure 2—figure supplement 2 and Figure 2—figure supplement 3).
The linear relationship between mCherry fluorescence and mRNA abundance of the tagged gene was present up to an expression level that corresponded to half of the abundance of ACT1 mRNA, which we had used as a reference gene in qPCR (Figure 2D). In previous RNA-seq data (Albert et al., 2018), 95% of S. cerevisiae genes fall below this threshold (Figure 2—figure supplement 4, Supplementary file 1). While differences in growth media (SC here and YNB in the RNA-seq data), as well as different qPCR efficiencies for ACT1 compared to GFP RNA (Figure 1—figure supplement 2), make this comparison imperfect, it suggests that the majority of genes can be quantified by our mRNA reporter. For lowly expressed genes, the GFP tag does not provide a strong enough signal to enable protein quantification (Huh et al., 2003; Newman et al., 2006, Figure 2—figure supplement 4). Based on these results, we concluded that our dual reporter system can be used to simultaneously measure mRNA and protein of more than half of the genes in the S. cerevisiae genome.
Simultaneous mapping of genetic variation affecting mRNA and protein levels
Our reporter system quantifies mRNA production and protein abundance at the same time, in the same live cells, exposed to the same environment. These features enable mapping of the genetic basis for variation in mRNA and protein levels, free from environmental or experimental confounders. We selected ten genes for genetic mapping (Figure 5—source data 1), based on several criteria. Five genes (ARO8, BMH2, GPD1, MTD1, UGP1) had previously been reported to have multiple differences between their respective eQTLs (Albert et al., 2018) and pQTLs (Albert et al., 2014b). Three genes (CYC1, OLE1, TPO1) had shown high agreement between their eQTLs and pQTLs. The remaining two genes (CTS1 and RPS10A) had low protein abundance based on GFP-tag quantification (Huh et al., 2003) compared to their mRNA levels (Albert et al., 2018).
To identify genetic loci affecting mRNA production and protein abundance, we used the strains BY4741 (BY), a reference strain frequently used in laboratory experiments, and RM11-1a (RM), a vineyard isolate closely related to European strains used in wine-making. These two strains differ at 47,754 variants in the yeast genome. We engineered RM to carry the CRISPR reporter inserted at the NPR2 gene and a synthetic genetic array (SGA) marker for selection of MATa haploid strains (Tong and Boone, 2007) at the neighboring CAN1 gene. We engineered a series of BY strains, each carrying one gene tagged with the GFP-gRNA tag (Figure 3). We crossed these BY strains to the RM strain and obtained populations of recombinant haploid progeny carrying both the tagged gene and the CRISPR reporter. Flow cytometry detected a range of GFP and mCherry signals from single cells (Figure 3).
To study the relationship between mRNA and protein among single cells, we examined the cell-to-cell correlation between mCherry and GFP fluorescence in our genetically heterogeneous populations during exponential growth (Figure 3—figure supplement 1A). After correcting for cell size (Figure 3—figure supplement 2), mCherry and GFP were positively correlated for all tested genes (Figure 3—figure supplement 1B). The strength of the correlation varied from gene to gene. Lower correlations between mCherry and GFP were observed for the genes with many differences between published eQTLs and pQTLs compared to those with more concordant eQTLs and pQTLs. Thus, different genes are influenced by mRNA-specific or protein-specific variation to different degrees.
Fluorescence-based readouts of mRNA and protein quantification in single cells enabled the use of bulk segregant analysis, a genetic mapping paradigm that gains statistical power from the analysis of millions of cells (Albert et al., 2014b; Ehrenreich et al., 2010). In each of the segregating populations, we used FACS to collect four subpopulations of 10,000 cells with high or low GFP or mCherry fluorescence, respectively, at a cutoff of 3–5% (Figure 3). In prior work, similarly stringent selection provided high power for QTL mapping (Albert et al., 2014b; Ehrenreich et al., 2010; Metzger and Wittkopp, 2019; Parts et al., 2014).
To gauge the heritability of gene expression among single cells, we measured fluorescence between the high and low populations after 13 generations of growth. In almost all cases, the sorted populations showed significant (T-test, p<10−5) differences in fluorescence, confirming the presence of genetic variants affecting mRNA and protein levels (Figure 3—figure supplement 3).
To map QTLs, we performed pooled whole-genome sequencing of all collected populations, computed the allele frequency of each DNA variant in each population, and calculated the difference in allele frequency (∆AF) between high and low populations along the genome. A significant ∆AF at a locus indicated the presence of one or more genetic variants affecting protein abundance (GFP) or mRNA production (mCherry, Figure 4—source data 1). We call loci mapped in this paper 'mRNA-QTLs' and 'protein-QTLs', to distinguish them from published 'eQTLs' and 'pQTLs'. QTL mapping was performed in two to six biological replicates for all but one gene (RPS10A). Because any allele frequency differences among replicate populations sorted on the same parameters (e.g. two high GFP populations for the same gene) represent false positives, we used the replicate data to estimate false discovery rates. We chose a significance threshold (logarithms of the odds; 'LOD'=4.5) corresponding to a false discovery rate of 7% (Figure 4—figure supplement 1). Between replicates, 76% of the protein-QTLs and 78% of the mRNA-QTLs were reproducible at genome-wide significance (Figure 4A).
Figure 4—source data 1
- https://cdn.elifesciences.org/articles/60645/elife-60645-fig4-data1-v2.xlsx
Figure 4—source data 2
- https://cdn.elifesciences.org/articles/60645/elife-60645-fig4-data2-v2.xlsx
Figure 4—source data 3
- https://cdn.elifesciences.org/articles/60645/elife-60645-fig4-data3-v2.xlsx
Across the ten genes, we detected 78 protein-QTLs and 44 mRNA-QTLs (Figure 4—source data 2 and Figure 4—source data 3). By design, all detected loci were trans-acting, and most were located on a different chromosome than the tagged gene. One locus located at ~450 kb on chromosome XIV affected mCherry levels in the same direction in all ten genes. This region was also observed in a control experiment, in which mCherry was expressed constitutively using an ACT1 promoter, and without a gRNA present (Figure 4—figure supplement 2). This region harbors the MKT1 gene, which carries a variant affecting a variety of traits (Deutschbauer and Davis, 2005; Fay, 2013). While the highly pleiotropic MKT1 locus may truly affect all ten genes we tested, it could also affect mCherry fluorescence via mCherry maturation or degradation, independently of any tagged gene. We excluded this region from further analyses.
The number of protein-QTLs per gene identified here (median = 7) agrees well with results from a previous study using the same mapping strategy (median = 8 for the same genes; Albert et al., 2014b), confirming that individual proteins are influenced by multiple trans-acting loci. The effects of individual protein-QTLs and pQTLs showed a positive correlation across studies (Pearson r = 0.73, p-value<10−15, Figure 4B). The number of mRNA-QTLs per gene in our study (median = 3 after removing the MKT1 locus) was lower than th from a previous study using RNA sequencing in 1012 segregants (median = 8 for the same genes; Albert et al., 2018). This difference could be due to using our reporter in single cells with high stochastic variation compared to RNA-seq in individually grown segregant cultures in the earlier study (see Discussion). However, while the mRNA-QTLs detected by our reporter primarily reflect influences on mRNA production, the RNA-seq-based eQTLs may reflect effects on transcription as well as mRNA degradation, which our system was not designed to capture. The effects of mRNA-QTLs were significantly correlated between studies (r = 0.44, p-value=5 × 10−6, Figure 4C). Some of the QTLs we detected harbored variants known to affect gene expression. For example, a region at ~650 kb on chromosome XII that contains the gene HAP1 affected protein abundance and/or mRNA production of GPD1, CYC1, OLE1, and TPO1 (Figure 4A). In the BY strain, the HAP1 coding sequence is interrupted by a transposon insertion, which alters the expression of thousands of mRNAs in trans (Albert et al., 2018; Brem et al., 2002). Overall, these agreements with previous analyses confirmed the reliability of our new reporter as a means for mapping the genetic basis of gene expression variation.
We detected several QTLs that were not shared with prior work and vice versa (Figure 4B-C). Most of these QTLs tended to have small-effect sizes, suggesting that they could have been missed due to incomplete power in either study. Alternatively, these QTLs may reflect experimental differences between studies, such as different growth media. For example, we observed a new, strong protein-QTL affecting Aro8 on chromosome XIV. The regulation of Aro8 expression by amino acid levels (Iraqui et al., 1998) suggests that this QTL could be due to the synthetic complete medium used here vs. YNB medium in earlier work.
Differences between mRNA-QTLs and protein-QTLs
Genetic mapping using our reporter enabled us to compare mRNA-QTLs and protein-QTLs, free from environmental or experimental confounders. We classified 86 loci based on the presence and effect direction of their respective mRNA-QTLs and/or protein-QTLs (Figure 5A and Figure 5—figure supplement 1, Figure 5—source data 2).
Of these 86 loci, 16 affected mRNA and protein of a given gene in the same direction. Such loci are expected for variants that alter a gene's mRNA production such that, in the absence of other effects, they also result in a concordant effect on protein abundance. A majority of the loci corresponded to protein-QTLs that did not overlap an mRNA-QTL. These 52 protein-specific QTLs may arise from variants that affect translation or protein degradation, without an effect on mRNA production.
There were eleven mRNA-QTLs that did not overlap with a protein-QTL and seven loci where mRNA-QTLs and protein-QTLs overlapped but had discordant effects. These two categories may occur when protein abundance and mRNA production of the same gene are regulated separately, through two different trans-acting pathways. These two pathways could be affected by two distinct but genetically linked causal variants at the same locus, or by a single variant with distinct pleiotropic effects on the two pathways. Alternatively, buffering mechanisms (Battle et al., 2015; Großbach et al., 2019) may compensate for changes in mRNA production perfectly (resulting in an mRNA-specific QTL) or may overcompensate (resulting in a discordant QTL pair) (Figure 5—figure supplement 1).
Genes differed widely in the complexity and specificity of trans-acting loci that influenced their expression. For example, four genes (BMH2, GPD1, UGP1, and CTS1) were each influenced by multiple loci, none of which affected mRNA and protein levels in the same direction. By contrast, most of the loci influencing CYC1 had concordant effects on mRNA and protein (Figure 5A).
While more than 73% of loci were specific for mRNA or protein, this difference might be inflated by loci that are truly concordant, but at which either the mRNA-QTL or the protein-QTL narrowly failed to meet the significance threshold. To bypass this potential limitation, we compared effect sizes, expressed as ∆AF, at significant mRNA-QTLs or protein-QTLs, irrespective of the significance of the locus in the other data (Figure 5B). When considering all loci, we observed a significant, positive correlation between mRNA and protein effect sizes (r = 0.41, p-value=8.4 × 10−5, Figure 5B). This overall correlation was almost exclusively driven by the concordant QTL pairs, whose effect sizes showed a strong correlation (r = 0.88 p-value=9 × 10−6). In sharp contrast, neither protein-specific QTLs (r = 0.2, p-value=0.23) nor mRNA-specific QTLs (r = −0.05, p-value=0.9) had correlated effects across the two data types, as expected if these loci specifically affected only mRNA or only protein. Considerable differences between mRNA-QTLs and protein-QTLs were also observed when simply considering effect directions. Overall, only 64% of QTLs affected mRNA and protein in the same direction. While this was more than the 50% expected by chance (binomial test p-value=0.006), it left 36% of loci with discordant effects. Protein-specific QTLs showed similar directional agreement (63%) at lower significance (p-value=0.04), while only 55% of mRNA-specific QTLs had an effect in the same direction in the protein data, which was not significantly different from chance (p-value=0.5). Together, these results are consistent with the existence of many QTLs that specifically affect mRNA production or protein abundance.
Several loci were shared across the ten genes. Even these shared loci differed in the specificity of their effects on mRNA or protein. For example, the locus containing the HAP1 gene had strong, concordant effects on both mRNA and protein for CYC1 and OLE1, but affected only the protein abundance of UGP1, and had significant but discordant effects on mRNA and protein for GPD1. Overall, these results revealed complex trans-acting influences on gene expression, in which genes were affected by different sets of multiple loci, with different degrees of mRNA or protein specificity.
A premature stop mutation in YAK1 affects gene expression post-transcriptionally
The causal variants in most trans-acting loci are unknown, limiting our understanding of the underlying mechanisms. In particular, very few causal variants with specific trans effects on protein abundance are known (Chick et al., 2016; Hause et al., 2014). We noticed a region at ~155 kb on chromosome X that affected the protein abundance but not mRNA production of ARO8, BMH2, and especially GPD1 (Figure 4A). This region spanned about 20 kb and contained 15 genes and 99 sequence variants. To identify the causal variant, we systematically divided this region into four tiles, swapped alleles in each tile using double-cut CRISPR-swap, an efficient scarless genome editing strategy (Lutz et al., 2019), and quantified the effect of these swaps on Gpd1-GFP fluorescence (Figure 6A-D).
This strategy, followed by analysis of our segregant population sequencing data, pinpointed a single G→A variant at 148,659 bp in the YAK1 gene as the causal variant. While this variant is present in neither the BY nor RM reference genomes (Figure 6E and Figure 6—figure supplement 1), our sequence data showed it to be present in all BY derivatives we used from the GFP collection (specifically, strains tagged at ARO8, BMH2, GPD1, MTD1, and UGP1; Figure 6—figure supplement 1; Huh et al., 2003). We observed this variant in two additional strains we genotyped from the GFP collection (FAA4 and YMR315W) and all four strains we genotyped from the tandem affinity purification (TAP)-tagged collection (PGM1, NOT5, EMI2, and TUB1) (Ghaemmaghami et al., 2003). This variant was not present in a BY4741 strain that we obtained from the ATCC stock center (#201388), suggesting that the YAK1 variant arose very recently in the specific BY4741 strain used to construct both the GFP and TAP-tagged collections. YAK1 encodes a protein kinase involved in signal transduction in response to starvation and stress, indirectly regulating the transcription of genes involved in various pathways (Figure 6E). The causal variant changes the 578th codon (glutamine) to a premature stop codon that is predicted to disrupt translation of the Yak1 kinase domain (Figure 6E).
The YAK1Q578* variant led to a diminution of Gpd1-GFP fluorescence, suggesting a decrease of Gpd1-GFP protein abundance (Figure 6D). While YAK1 may control transcription of genes in the glycerol biosynthesis pathway (Lee et al., 2008; Rep et al., 2000), which includes GPD1, our QTL results suggested no link between the variant and GPD1-GFP mRNA level. Consistent with a protein-specific trans-effect on GPD1, deletion of YAK1 in a strain in which GPD1 was tagged with GFP-gRNA caused a reduction of green fluorescence but had no detectable effect on mCherry fluorescence (Figure 6F). Further, qPCR indicated no difference in the level of GPD1-GFP mRNA in YAK1Q578* or yak1Δ compared to matched YAK1wt (Figure 6F).
We explored the genome-wide effects of the YAK1 variant by comparing mRNA and protein abundance between strains carrying the YAK1Q578* and the BY YAK1 wildtype allele using RNA-seq and mass spectrometry (Figure 7A, Figure 7—source data 1, Figure 7—source data 2 and Figure 7—source data 3). Among 5755 quantified mRNA transcripts, 262 were up-regulated and 310 down-regulated in the presence of YAK1Q578* (Benjamini-Hochberg (BH) adjusted p-value<0.05) (Benjamini and Hochberg, 1995). The variant reduced the abundance of 82 of 2590 quantified proteins, and increased another 82 proteins (BH adjusted p-value<0.05). By comparing mass spectrometry and RNA-seq results, we classified genes as affected only at the mRNA level (58 genes up, and 118 genes down-regulated), only at the protein level (60 genes up, and 50 genes down-regulated), or at both mRNA and protein (15 genes up, and 27 genes down-regulated). There was a strong enrichment for genes involved in cytoplasmic translation (q-value <10−10) among the genes with reduced mRNA abundance, which is consistent with the role of Yak1 as a regulator of transcription of ribosomal genes through Crf1 phosphorylation (Martin et al., 2004; Figure 7—figure supplement 1, Figure 7—source data 2). Genes up-regulated at the mRNA level showed an enrichment in amino acid biosynthesis (q-value = 0.001). The most differentially expressed genes included known targets of the YAK1 pathway (Figure 7A–B). Gpd1 protein was strongly reduced (BH adjusted p-value<0.004), with a non-significant effect on GPD1 mRNA (adjusted p-value=0.10) (Figure 7A).
Figure 7—source data 1
- https://cdn.elifesciences.org/articles/60645/elife-60645-fig7-data1-v2.xlsx
Figure 7—source data 2
- https://cdn.elifesciences.org/articles/60645/elife-60645-fig7-data2-v2.xlsx
Figure 7—source data 3
- https://cdn.elifesciences.org/articles/60645/elife-60645-fig7-data3-v2.xlsx
Figure 7—source data 4
- https://cdn.elifesciences.org/articles/60645/elife-60645-fig7-data4-v2.xlsx
Finally, we investigated if the YAK1Q578* mutation affects other phenotypes. As YAK1 and GPD1 are involved in osmotic stress resistance (Lee et al., 2008), we grew strains carrying YAK1wt, YAK1Q578*, and yak1Δ, in a range of sodium chloride concentrations (Figure 6—figure supplement 2A). While this osmotic stress reduced growth, strains with YAK1Q578* and yak1Δ had a higher growth rate than wild-type, consistent with the role of Yak1 in triggering cell cycle arrest in response to stress. Gpd1-GFP abundance increased with stronger osmotic stress in YAK1wt and yak1Δ, with consistently lower expression of Gpd1-GFP in yak1Δ (Figure 6—figure supplement 2B–C).
We developed a fluorescence-based dual reporter system for the simultaneous quantification of mRNA and protein from a given gene in single, live cells. This system enabled the use of a statistically powerful mapping strategy to identify genetic loci that affected mRNA production or protein abundance in trans. Because mRNA and protein were quantified in exactly the same condition, mRNA-QTLs and protein-QTLs could be compared without environmental confounding. This reporter system and mapping strategy were explicitly designed to maximize the chance of detecting concordant effects on mRNA and protein. Nevertheless, of the 86 trans-acting loci that contributed to variation in the expression of ten genes, 84% did not have concordant effects on mRNA production and protein abundance. This result demonstrates the importance of variants that act on specific layers of gene expression.
The genes ARO8, BMH2, GPD1, MTD1, and UGP1, which we had selected for high discrepancy between previous mRNA-QTLs and protein-QTLs, showed many QTLs (89%) that were not shared between mRNA and protein in our data. In contrast, CYC1, OLE1, and TPO1, which we had selected for higher agreement between published QTLs, showed fewer discrepant QTLs in our data, although even for these genes the majority of QTLs was not shared between mRNA and protein (58%). These three genes showed fewer QTLs overall and all had one locus with strong effect size (Figure 4A; the HAP1 locus for CYC1 and OLE1, and the IRA2 locus (Brem et al., 2002; Smith and Kruglyak, 2008) for TPO1). Based on these results, strong trans-acting loci may be more likely to cause concordant effects on mRNA and protein, while loci with smaller effects could be more likely to be specific to mRNA or protein.
While half of the mRNA-QTLs we detected had concordant effects on protein (16 out of 34), most protein-QTLs had no effects on mRNA (52 out of 75), in line with observations from prior work (Albert et al., 2014b; Albert et al., 2018; Foss et al., 2011). That 70% of our protein-QTLs had protein-specific effects suggests that the causal variants underlying many of these loci affect post-transcriptional processes, i.e. mRNA stability, translation, or protein degradation. The indirect nature of our mCherry reporter and its lower signal intensity compared to GFP fluorescence are potential sources of measurement noise, which could have reduced detection power for mRNA-QTLs compared to protein-QTLs. However, our analyses of the magnitudes and directions of effects on mRNA and protein, which did not require loci in the other data to meet a significance threshold, also suggested that many protein-QTLs specifically influence protein.
We found seven loci that had discordant effects on mRNA production and protein abundance of the same gene. For example, at the HAP1 locus, the BY allele increased Gpd1 protein abundance but decreased GPD1 mRNA production, as had been seen when comparing QTLs across different studies (Albert et al., 2018). The highly pleiotropic effects of HAP1 on mRNA and protein levels of many genes (Albert et al., 2014b; Albert et al., 2018; Smith and Kruglyak, 2008) reinforces the hypothesis that HAP1 alleles influence Gpd1 protein abundance and mRNA production via two separate trans-acting mechanisms.
QTLs with discordant effects on mRNA and protein, as well as mRNA-specific QTLs, may be caused by buffering of mRNA variation at the protein level. A well-studied example of this phenomenon is the regulation of expression of genes that encode members of a protein complex, in which excess protein molecules that cannot be incorporated in the complex are rapidly degraded (Taggart et al., 2020). Among the genes we investigated, RPS10A encodes a part of the ribosome small subunit complex. RPS10A showed the highest number of mRNA-specific QTLs, possibly because Rps10a is subject to buffering mechanisms.
The nonsense mutation (Q578*) we identified in the YAK1 gene provides an informative example of the complexity with which trans-acting variants can shape the transcriptome and the proteome. YAK1Q578* changed protein abundance for many genes more strongly than mRNA abundance, but also affected mRNA but not protein for many other genes. Thus, the consequences of this mutation span mechanisms that affect mRNA as well as protein-specific processes. A reduction of ribosomal gene transcription may account for some of these observations by reducing the translation of multiple genes.
The YAK1Q578* variant likely arose as a new mutation in the BY4741 ancestor of the GFP and TAP-tagged collections. Its relatively large effect and rarity in the global yeast population are consistent with population genetic expectations (Eyre-Walker, 2010; Gibson, 2012) and observations (Bloom et al., 2019; Fournier et al., 2019) for a deleterious variant that may have drifted to fixation in this specific background, as has been suggested for many causal variants in natural yeast isolates (Warringer et al., 2011). Alternatively, faster growth of a strain carrying the YAK1Q578* variant during osmotic stress (Figure 5—figure supplement 1) may have contributed to adaptive fixation of this variant in this specific isolate of BY4741. While the large effect of YAK1Q578* aided our ability to fine-map it (Rockman, 2012), we suspect that its diverse, mRNA-specific as well as protein-specific mechanistic consequences may be representative of more common trans-acting variants with smaller effects.
To simultaneously quantify mRNA and protein and eliminate potential environmental confounders in expression QTL mapping, we developed a system in which a gRNA drives CRISPR activation of a fluorescent reporter gene in proportion to a given mRNA of interest. Standard methods for mRNA quantification require lysis of cell cultures or tissues, constraining sample throughput and statistical power for mapping regulatory variation. Single-cell RNA-seq (Picelli, 2017; Tang et al., 2009) and in situ fluorescent hybridization techniques (Buxbaum et al., 2014; Player et al., 2001; Rouhanifard et al., 2018) are improving rapidly, including in yeast (Gasch et al., 2017; Li and Neuert, 2019; Nadal-Ribelles et al., 2019; Wadsworth et al., 2017). However, these approaches still have throughput that is orders of magnitude below that available using FACS. Further, they involve cell lysis or fixation, precluding bulk segregant approaches that rely on growing cells after sorting. By contrast, our reporter allows quantification of mRNA production of a given gene within millions of single, live cells by flow cytometry. Because mCherry production in our system amplifies mRNA abundance signals, it is able to quantify genes with mRNA levels that would likely be hard to detect by FACS using in situ hybridization methods. The readout of our system is driven by a gRNA after it detaches from its mRNA. Therefore, the resulting signal is independent of the fate of the mRNA after gRNA release. Given the hammerhead ribozyme has a rate constant for self-cleavage of 1.5 per minute (Wurmthaler et al., 2018), gRNA abundance is not expected to reflect the half-life and stability of most yeast mRNAs, which have a median half-life of 3.6 min (Chan et al., 2018). By contrast, standard methods usually used in eQTL mapping quantify mRNA at steady state, which may explain some of the differences we observed between our mRNA-QTLs and known eQTLs identified by RNA-seq.
Future versions of our reporter could be improved in several ways. The mCherry used here has a maturation time of about 40 min (Merzlyak et al., 2007), which limits the temporal resolution at which we can observe dynamic expression changes. Fluorescent proteins with faster maturation times and shorter half-lives, could enable precise investigation of rapid temporal change in mRNA production. Using brighter fluorescent proteins or multiple copies of mCherry and its gRNA-targeted promoter could increase fluorescence and increase mRNA detection further. Finally, we observed that beyond a certain mRNA level, the abundance of gRNA no longer follows that of the tagged mRNA. Nevertheless, we estimated that our CRISPR-based reporter can be used to quantify the mRNA production of most S. cerevisiae genes. Like other reporters that tag the 3' end of genes, our system decouples the native 3' UTR from the gene of interest. 3' UTRs can influence gene expression (Shalem et al., 2015), and genetic trans-effects that involve mechanisms that target the 3' UTR would be missed by our system. However, loss of the native 3' UTR is not expected to inflate differences between mRNA-QTLs and protein-QTLs. Instead, 3' UTR removal may reduce differences between mRNA-QTLs and protein-QTLs because it eliminates an element of the mRNA that could be a target for mRNA-specific processes. Finally, CRISPR activation has been demonstrated in many organisms (Long et al., 2015; Maeder et al., 2013; Park et al., 2017), suggesting that similar reporters for RNA production could be developed in other species.
Our reporter system for quantifying mRNA and protein of a given gene in the same live, single cells combined with a mapping strategy with high statistical power was deliberately designed to minimize technical or environmental confounders that may have inflated differences between the genetics of mRNA and protein levels in earlier work. Yet, fewer than 20% of the detected trans-acting loci we detected had concordant effects on mRNA and protein levels, providing strong support for the existence of discordant trans effects on mRNA vs proteins. Whether cis effects show similar degrees of specificity will be a topic of future work. The fact that the majority of trans-QTLs identified here were protein-specific suggests that protein abundance is under more complex genetic control than mRNA abundance.
Materials and methods
Yeast strains
Request a detailed protocolWe used 160 yeast strains, 12 of which were obtained from other laboratories, including six from the GFP collection, and 148 that were built for this study (complete list in Supplementary file 2). All strains are based on two distinct genetic backgrounds: BY4741 (BY), which is closely related to the commonly used laboratory strain S288c, and RM11-1a (RM), a haploid offspring of a wild isolate from a vineyard, which is closely related to European strains used in wine-making. Both strains carried auxotrophic markers, and RM had been engineered earlier to facilitate laboratory usage (BY: his3∆one leu2∆0 met15∆0 ura3∆0; RM: can1Δ::STE2pr-URA3 leu2Δ0 HIS3(S288C allele) ura3Δ ho::HYG AMN1(BY allele); Supplementary file 2). Most strains were built using conventional yeast transformation (Gietz and Schiestl, 2007) and DNA integration based on homologous recombination. Integrated DNA fragments were produced by PCR (Phusion Hot Start Flex NEB M0535L, following manufacturer protocol, annealing temperature: 57°C, 36 cycles, final volume: 50 µl) and gel purified (Monarch DNA Gel Extraction Kit, NEB T1020L), with primers carrying 40 to 60 bp overhanging homologous sequence as required. All primers are available in Supplementary file 3. For transformation, fresh cells from colonies on agar plates were grown in YPD media (10 g/l yeast extract, 20 g/l peptone, 20 g/l glucose) overnight at 30°C. The next day, 1 ml of the culture was inoculated in an Erlenmeyer flask containing 50 ml of YPD and grown under shaking at 30°C for 3 hr to reach the late log phase. Cells were harvested by centrifugation and washed once with pure sterile water and twice with transformation buffer 1 (10 mM TrisHCl at pH8, 1 mM EDTA, 0.1 M lithium acetate). We resuspended the cells in 100 µl of transformation buffer 1, added 50 µg of denatured salmon sperm carrier DNA (Sigma #D7656) and 1 µg of the DNA fragment to be integrated, and incubated for 30 min at 30°C. Alternatively, when transforming a replicative plasmid, we used 0.1 µg of plasmid DNA and skipped this first incubation. We added 700 µl of transformation buffer 2 (10 mM TrisHCl at pH8, 1 mM EDTA, 0.1 M lithium acetate, 40% PEG 3350) and performed a second incubation for 30 min at 30°C. A heat shock was induced by incubating the cells at 42°C for 15 min. The transformed cells were then washed twice with sterile water. If the selective marker for the transformation was an antibiotic resistance gene, the cells were resuspended in 1 ml of YPD, allowed to recover for 2 hr at 30°C, and spread on a YPD plate (2% agar) containing the antibiotic (200 ng/l G418, 100 ng/l nourseothricin sulfate/CloNAT, or 300 ng/l hygromycin B). Alternatively, if the transformation was based on complementation of an auxotrophy, the cells were resuspended in 1 ml of sterile water and spread on a plate containing minimal media lacking the corresponding amino acid or nucleotide (YNB or Synthetic Complete (SC): 6.7 g/l yeast nitrogen base (VWR 97064–162), 20 g/l glucose, with or without 1.56 g/l SC -arginine -histidine -uracil -leucine (Sunrise science 1342–030), complemented as needed with amino acids: 50 mg/l histidine, 100 mg/l leucine, 200 mg/l uracil, 80 mg/l tryptophan). After two to three days of incubation at 30°C, colonies were streaked on a fresh plate containing the same selection media to purify clones arising from single, transformed cells. DNA integration in the correct location was confirmed by PCR (Taq DNA Polymerase NEB M0267L, following manufacturer protocol, annealing temperature: 50°C, 35 cycles, final volume: 25 µl, primers in Supplementary file 3). To store the constructed strains, we regrew the validated colony on a new selection media plate overnight at 30°C, scraped multiple colonies, resuspended the cells in 1.4 ml of YPD containing 20% glycerol in a 2 ml screw cap cryotube and froze them at −80°C.
We constructed seven plasmids: three plasmids that do not replicate in yeast and that carry the GFP-gRNA tag, the CRISPR reporter, and Z3EV system, respectively (Figure 1—figure supplement 1), and four yeast-replicating plasmids to investigate the quantitative properties of our reporter (Figure 2—figure supplement 3). These plasmids were constructed through multiple rounds of cloning using DNA fragments from yeast DNA or plasmids acquired from Addgene (kind gifts from John McCusker: Addgene #35121–22, from Michael Nick Boddy: Addgene #41030, from Benjamin Glick: Addgene #25444, from Timothy Lu: Addgene #64381, #64389, #49013; complete list of plasmids in Supplementary file 4). Plasmids were assembled using Gibson assembly (NEBuilder HiFi DNA Assembly Cloning Kit, NEB E5520S). Fragments were either PCR amplified with a least 15 bp overlap at each end (Phusion Hot Start Flex NEB M0535L, manufacturer protocol, annealing temperature: 57°C, 36 cycles, final volume: 50 µl, primers in Supplementary file 3) or obtained by restriction digestion of already existing plasmids (also shown in Supplementary file 3).
The fragment encoding the gRNA tag, containing the two ribozymes and the gRNA sequence itself, was purchased as a 212 bp double-stranded DNA oligo from IDT (we used the 'C3' gRNA from Farzadfard et al., 2013, as it was reported to provide the highest reporter gene expression). The synthetic polyA tail following the GFP sequence (Figure 1A) was introduced by using a PCR primer containing 45 thymines in its overhang sequence (primer OFA0038 in Supplementary file 3). Fragments for assembly were purified using agarose electrophoresis and gel extraction (Monarch DNA Gel Extraction Kit, NEB T1020L). For assembly, the given fragments were mixed at equi-molar amounts of 0.2–0.5 pM in 10 µl. Assembly was done by addition of 10 µl of NEBuilder HiFi DNA Assembly Master Mix and incubation at 50°C for 60 min. From this reaction, 2 µl of the final products were transformed into E. coli competent cells (10-beta Competent E.coli, NEB C3019I) through an incubation of 30 min on ice and a heat shock of 30 s at 42°C. Transformed cells were spread on LB plates (10 g/l peptone, 5 g/l yeast extract, 10 g/l sodium chloride, 2% agar) containing 100 mg/l ampicillin and grown overnight at 37°C. After cloning, the final plasmids were extracted (Plasmid Miniprep Kit, Zymo Research D4036) and verified by restriction enzyme digestion or PCR (Taq DNA Polymerase NEB M0267L, 30 cycles, 25 µl final volume, primers in Supplementary file 3). We also verified by Sanger sequencing that the gRNA tag in the plasmid had no mutation. To store the plasmids, the host bacteria were grown in LB with ampicillin overnight at 37°C and 1 ml of the culture was mixed with 0.4 ml of a sterile solution containing 60% water and 40% glycerol. The cells were stored at −80°C. The three plasmids containing the different parts of the reporter are available on Addgene (ID #157656, #157658, and #157659) along with their full DNA sequence.
Plate reader-based fluorescence measurements
Request a detailed protocolYeast fluorescence was measured in 24 hr time courses during overnight growth in a BioTek Synergy H1 plate-reader (BioTek Instruments). Fresh cells from agar plates were inoculated in 100 µl of minimal YNB media (Figure 1) or SC media (all other figures) containing any complements necessary for growth of auxotrophic strains, at an initial optical density at wavelength 600 nm (OD600) of 0.05 in a 96-well flat bottom plate (Costar #3370). We conducted the experiments in Figure 1 in the presence of G418 (which requires the use of YNB + glutamate) because, at the time, we were concerned about potential loss of the reporter construct from the genome due to recombination between its flanking sequences. This concern later turned out to be unfounded. The plates were sealed with a Breathe Easy membrane (Diversified Biotech). Cells were grown in the plate reader at 30°C and with circular agitation in between fluorescence acquisition. During each acquisition, performed every 15 min, we recorded OD600, GFP fluorescence (read from bottom, excitation 488 nm, emission 520 nm, 10 consecutive reads averaged, gain set to 'extended') and mCherry fluorescence (read from bottom, excitation 502 nm, emission 532 nm, 50 consecutive reads averaged, gain set to a value of 150). We took 97 measurements during 24 hr of growth, unless individual runs were manually terminated early.
Raw measurements of OD600 and fluorescence were processed using R version 3.5.1 (https://www.r-project.org/, scripts and raw data available at the github repository at Brion, 2020; https://github.com/BrionChristian/Simultaneous_RNA_protein_QTLs; copy archived at swh:1:rev:74628305678f30d5b126745e5c6cd6cec4091c12). 'Blank' values from wells with no cells were subtracted from OD600 and fluorescence measurements of wells that had been inoculated with cells. OD600 was log-transformed and manually inspected to identify the late log phase, i.e. a time point about 3/4 into the exponential growth phase. This stage was identified separately for each well, and usually corresponded to an OD600 of 0.1–0.3. We extracted the OD600 and fluorescence measurements at the five time points centered on our selected time point. The mCherry and GFP fluorescence ratios were calculated as the ratio between the fluorescence and the OD600 at these five time points (example in Figure 1B–CB), allowing us to estimate fluorescence while correcting for culture density. Focusing on the late log phase allowed measurements at higher cell density to provide more robust fluorescence reads. Growth rates were estimated as the slope of a linear fit of the log of OD600 over time.
RNA quantification by qPCR
Cell harvest
Request a detailed protocolWe quantified mRNA and gRNA abundance by quantitative real-time reverse-transcription PCR of RNA extracted from exponentially growing cells. Cells were grown in either 50 ml of medium (YNB with auxotrophic complements, results shown in Figure 1D) in shaking Erlenmeyer flasks or in 1.2 ml of media (SC with auxotrophic complements and estradiol (Sigma, E1024-1G), Figure 2C–EC and Figure 6F) in a shaking 2 ml 96-deep-well plate. The OD600 was monitored to identify the second half of the exponential growth phase (corresponding to an OD600 of 0.35–0.45 OD600 in flasks, and 0.20–0.30 in the deep-well plates). At this point, GFP and mCherry fluorescence ratios were recorded in a BioTek Synergy H1 plate reader (BioTek Instruments). Cells were then harvested immediately. Cells were washed with sterile water through either short centrifugation using 5 ml of culture from flasks, or vacuum-filtration through a 96-well filter plate (Analytical Sales 96110–10) using the entire remaining 1 ml of culture from the deep-well plate. Cells were then immediately flash-frozen in either isopropanol at −80°C (pellet from flask) or liquid nitrogen (filter plate) and stored at −80°C until RNA extraction.
RNA extraction from flasks
Request a detailed protocolTo extract the RNA from cells grown in flasks, we used the ZR Quick-RNA Kit (Zymo Research R1054). Frozen cell pellets were resuspended in 800 μl RNA Lysis Buffer one from the kit and transferred to a ZR BashingBead Lysis Tube. The cells were shaken in a mini-bead beater (BioSpec Products) for ten cycles of one minute in the beater, one minute on ice. Cell debris and beads were centrifuged for one minute at full speed and 400 µl of supernatants were transferred into Zymo-Spin IIIC Columns. The columns were centrifuged for one minute, and 400 µl 100% ethanol was added to the flow-through. After mixing, the flow-throughs were transferred into Zymo-Spin IIC Columns and centrifuged for 1 min to bind the RNA and DNA to the columns. The columns were washed with 400 μl RNA Wash Buffer from the kit. DNA was digested in columns by adding a mixture of 5 μl DNase I and 75 μl DNA Digestion Buffer from DNase I Set kit (Zymo Research E1010) followed by a 15 min incubation at room temperature. The columns were then washed three times with 400 μl RNA Prep Buffer, 700 μl RNA Wash Buffer, and 400 μl RNA Wash Buffer. RNA was eluted in 50 µl DNase/RNase-free water, quantified using Qubit RNA BR or HS Assay Kit (Thermo Fisher Scientific Q10210 or Q52852), and stored at −20°C.
RNA extraction from 96-well plates
Request a detailed protocolTo extract the RNA from cells grown in 96-well plates, we used the ZR RNA in-plate extraction kit (ZR-96 Quick-RNA Kit, Zymo Research R1052), which followed the same protocol as the flask RNA extraction above, with a few minor differences. Bead-beating was done in an Axygen 1.1 ml plate (P-DW-11-C-S) with 250 µl of acid washed 425–600 µm beads (Sigma G8722) per well, sealed with an Axymat rubber plate seal (AM-2ML-RD-S). RNA purified from 200 µl of the resulting supernatant. DNA digestion and washing steps were done on Silicon-A 96-well plates from the kit. The RNA was eluted in 30 µl of DNase/RNase-free water, quantified, and stored at −20°C.
Reverse-transcription and qPCR
Request a detailed protocolRNA was reverse-transcribed using the GoScript RT kit (Promega A5000) following the kit protocol. We performed negative controls, no-enzyme and no-primer, which generated no qPCR signals. Quantitative PCRs were done in a 96-well plate (Bio-Rad HSP9645) using GoTaq qPCR kit (Promega A6001). Plates were sealed using a microseal 'B' Adhesive Seal (Bio-Rad MSB1001) and the reaction progress was recorded during 40 cycles using a C1000Touch plate reader (Bio-Rad). We quantified four different parts of the tag cDNA (GFP, Hh ribozyme cleavage, gRNA, and HDV ribozyme cleavage, Figure 1D), as well as ACT1 cDNA as a reference gene. Primer sequences are in Supplementary file 3. The primers were tested and calibrated by running qPCR measurements on nuclear DNA extracts at a range of known input concentrations (Figure 1—figure supplement 2).
Segregant populations
Request a detailed protocolBY strains (BY4741 background) carrying a given GFP-gRNA-tag and RM (YFA0198) carrying the CRISPR reporter and the SGA marker were mixed for crossing on a plate with medium that allows only hybrids to grow (SC agar -leucine -histidine). Growing cells were streaked on the same medium, and a single hybrid colony was kept for storage and for generating the segregant population. For sporulation, hybrid strains were incubated in sporulation medium (2.5 g/l yeast extract, 2.5 g/l glucose, 15 g/l potassium acetate, 200 mg/l uracil, 100 mg/l methionine) at room temperature under vertical rotation in a glass tube for seven days. After verifying sporulation under a light microscope, 1 ml of medium containing the tetrads was pelleted (13,000 rpm for 5 min) and resuspended in 300 µl of sterile water containing about 15 µg of zymolyase. The resulting ascii were digested at 30°C for 30 min with agitation. Spores were separated by vortexing for about 15 s, and 700 µl of pure sterile water was added to the tube. We spread 250 µl of this spore suspension on a plate containing segregant selection media (SC agar, 50 mg/l canavanine, -uracil -leucine) allowing growth of haploid segregants carrying the following three alleles: (1) cells with mating type MATa, selected via the SGA marker with URA3 under control of the STE2 promoter, which resulted in a ura+ phenotype only in MATa cells, (2) the SGA marker integrated at the CAN1 gene (whose deletion conferred canavanine resistance), which also selected for the CRISPR reporter that we had integrated at NPR2, the gene next to CAN1, (3) the given gene of interest tagged with the GFP-gRNA tag and LEU2 selectable marker. We used SC because this medium complemented a growth defect of strains with an NPR2 deletion in YNB medium. After three days of incubation at 30°C, segregants were harvested by scraping the entire plate in 10 ml of sterile water. Cells were centrifuged, resuspended in 3 ml of segregant selection media, and incubated at 30°C for 1.5 hr. To store these genetically diverse segregant populations, 1 ml of the culture was mixed with 0.4 ml of a sterile solution containing 60% water and 40% glycerol in a 2 ml screw cap cryo tube and frozen at −80°C.
Cell sorting for QTL mapping
Request a detailed protocolOne day before cell sorting, the segregant population was thawed from the −80°C stock, mixed well, and 8 µl of culture were used to inoculate 5 ml of segregant selection media. The cells were reactivated with an overnight growth at 30°C under shaking. The next day, 1 ml of the growing culture was transferred to a new tube containing 4 ml of segregant selection medium (see above) and grown for an additional two hours before cell sorting, roughly corresponding to the middle of the exponential growth phase.
Cell sorting was performed on a BD FACSAria II P0287 (BSL2) instrument at the University of Minnesota Flow Cytometry Resource (UFCR), without compensation. Cells were gated to exclude doublet and cellular fragments. To focus on cells in approximately the same stage of the cell cycle, an additional gate selected cells in a narrow range of cell size as gauged by the area of the forward scatter signal (FSC). From the cells within this gate, we sorted five populations per experiment, each comprising 10,000 cells: (1) a control population from the entire gate without fluorescence selection, (2) the 3% of cells with the lowest GFP fluorescence, (3) the 3% of cells with the highest GFP fluorescence, (4) the 3% of cells with the lowest mCherry fluorescence, and (5) the 3% of cells with the highest mCherry fluorescence. Each population was collected into 1 ml of segregant selection medium. After overnight growth at 30°C, 0.9 ml of culture was mixed with 0.4 ml of a sterile solution containing 60% water and 40% glycerol, and frozen at −80°C until sequencing. The remaining 0.1 ml were inoculated into 0.9 ml of segregant selection medium and grown for 3 hr before analyzing the population using flow cytometry (see below). In total, we obtained 125 sorted populations from 25 experiments across the ten tagged genes, with 1 to 6 biological replicates per gene, as well as the untagged population (Figure 5—source data 1). Sorting was done in four batches on different dates. Biological replicates were performed as independent sporulations of the stored diploid hybrids, and thus represent independent populations sorted in separate experiments.
Flow cytometry
Request a detailed protocolSingle-cell fluorescence analysis was performed using cultures in the late log growth phase. We used a BD Fortessa X-30 H0081 flow cytometer at UFCR equipped with blue and yellow lasers and 505LP and 595LP filters to measure green (GFP) and red (mCherry) fluorescence, respectively. We did not use compensation. Forward scatter (FSC), side scatter (SSC), GFP, and mCherry fluorescence were recorded for 50,000 cells, excluding doublets and cellular debris. The voltaic gains were set as follows: 490 for FSC, 280 for SSC, 500 for GFP, and 600 for mCherry. We monitored for possible cross-contamination from cells retained in the instrument using strains expressing either only GFP or mCherry, and observed no cross-contamination. Recorded data on. fsc files were analysed using R and the flowCore package (Hahne et al., 2009). Raw data and scripts are accessible on github (https://github.com/BrionChristian/Simultaneous_RNA_protein_QTLs). The data were filtered to discard outlier cells based on unusual FSC and SSC signals. We used the fluorescence data from the sorted populations to determine correlations between red and green fluorescence, as well as heritability (Figure 3—figure supplement 1 & Figure 3—figure supplement 2). For these analyses, fluorescence values were corrected for cell size (FSC) by calculating the residuals of a loess regression of fluorescence on FSC. Loess regression avoided the need to assume a specific mathematical relationship between the two parameters (Figure 3—figure supplement 3).
DNA extraction and sequencing
Request a detailed protocolDNA extraction for whole-genome sequencing was performed in 96-well plate format using E-Z-96 Tissue DNA kits (Omega D1196-01). The stored, sorted populations were thawed, mixed, and 450 µl transferred into a 2 ml 96-deep-well plate containing 1 ml of segregant selection medium for an overnight growth at 30°C. The plate was centrifuged for 5 min at 3700 rpm, and the supernatant was removed by quick inversion of the plate. Then, 800 µl of Buffer Y1 (182 g/l sorbitol, 0.5 M EDTA, pH 8, 14.3 mM β-mercaptoethanol, 50 mg/l zymolyase 100T) were added to the pellets, and the cells were resuspended and incubated for 2 hr at 37°C. The spheroplasts were centrifuged and the supernatant discarded. The pellets were resuspended in 200 µl of TL buffer and 25 µl of OB Protease Solution from the kit and incubated overnight at 56°C. The next day, RNA was denatured by addition of 5 µl of RNAse A (20 mg/ml) and incubated at room temperature for 5 min. After addition of 450 µl of BL Buffer from the kit, the mixture was transferred onto a E-Z 96 column DNA plate and centrifuged at 3700 rpm for 3 min. The columns were washed once with 500 µl of HBC Buffer and three times with 600 µl of DNA Wash Buffer from the kit. After an additional centrifugation to dry the column, the DNA was eluted in 100 µl of pure sterile water, quantified using Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific Q32854) and stored at 4°C for library preparation the next day.
Library preparation for Illumina sequencing was performed using Nextera DNA Library Prep kit (Illumina) with modifications. The tagmentation was done on 5 ng of DNA using 4 µl of Tagment DNA buffer ('TD' in the kit) and 0.25 µl of Tagment DNA enzyme (corresponding to a 20-fold dilution of 'TDE1' from the kit) and incubating for 10 min at 55°C. Fragments were amplified with index primers (8 Nextera primers i5 and 12 Nextera primers i7, for up to 96 possible multiplex combinations) on 10 µl tagmented DNA by adding 1 µl of each primer solution (10 µM), 5 µl of 10X ExTaq buffer and 0.375 µl of ExTaq polymerase (Takara) and water to a final volume of 50 µl. The amplification was run for 17 PCR cycles (95°C denaturation, 62°C annealing, 72°C elongation). 10 µl of each reaction were pooled for multiplexing and run on a 2% agarose gel. DNA that migrated between the 400 and 600 bp was extracted using Monarch DNA Gel Extraction Kit (NEB T1020L). The pooled library DNA concentration was determined using Qubit dsDNA BR Assay Kit (Thermo Fisher Scientific Q32853), and submitted for sequencing. Sequencing was performed at the University of Minnesota Genomics Core (UMGC). Our 125 populations were processed in four batches extracted and sequenced at different times. Two were sequenced using an Illumina HiSeq 2500 (high-output mode; 50 bp paired-end) and two were sequenced using an Illumina NextSeq 500 (mid-output mode, 75 bp paired-end). Read coverage ranged from 5-fold to 24-fold coverage of the genome (median: 13-fold). The reads are available on NCBI SRA via BioProject PRJNA644804.
QTL mapping
Request a detailed protocolFor each sorted and sequenced population, reads were filtered (MAPQ ≥30) and aligned to the S. cerevisiae reference genome (version sacCer3, corresponding to BY strain) using BWA (Li and Durbin, 2009, command: mem -t). We used samtools (Li et al., 2009, command: view -q 30) to generate bam files and collapse PCR duplicates using the rmdup command. We used 18,871 variants previously identified as polymorphic and reliable between RM and BY (Bloom et al., 2013; Ehrenreich et al., 2010) (list available on github: https://github.com/BrionChristian/Simultaneous_RNA_protein_QTLs, samtools: mpileup -vu -t INFO/AD -l), generating vcf files with coverage and allelic read counts at each position for each population.
The vcf files were processed in R to identify bulk segregant analysis QTLs using code adapted from Albert et al., 2014b (available on github: https://github.com/BrionChristian/Simultaneous_RNA_protein_QTLs). Briefly, for plotting the results, the allele frequency of the reference (that is, BY) allele was calculated at each position in each population. Random counting noise was smoothed using loess regression, and the allele frequency of a given 'low' fluorescence population subtracted from its matched 'high' fluorescence population to generate ΔAF. A deflection from zero indicated the presence of a QTL. To identify significant QTLs, we used an R script that implemented the MULTIPOOL algorithm (Edwards and Gifford, 2012), which calculates LOD score based on ΔAF and depth of read coverage in bins along the genome. We used MULTIPOOL output to call QTLs as peaks exceeding a given significance threshold (see below), along with confidence intervals for the peak location corresponding to a 2-LOD drop from the peak LOD value. We applied the MULTIPOOL algorithm using the following parameters: bp per centiMorgan: 2,200; bin size: 100 bp, effective pool size: 1000. We excluded variants with extreme allele frequencies of >0.9 or <0.1. We initially set a permissive detection threshold of LOD >3.0 to identify a set of candidate QTLs, which we then integrated across replicates (507 QTLs, Figure 4—source data 1). A second, more stringent, threshold of LOD >4.5 was then applied to retain only significant QTLs based on our estimated false discovery rate (FDR).
To estimate FDR, we applied the multipool QTL detection algorithm to pairs of populations sorted into the same gates in different replicates. Any 'QTLs' in such comparisons must be due to technical or biological noise. We restricted these analyses to replicates sequenced on the same instrument, resulting in 80 inter-replicate comparisons. From these data, we calculated the FDR as a function of the significance threshold (thr): FDRthr = (NrepQTLthr/Nrep) / (NfluoQTLthr/Nfluo), where NrepQTLthr is the number of false 'QTLs' from comparing the same gate across replicate populations at a LOD score threshold of thr, Nrep is the number of such inter-replicate comparisons (Nrep = 80), NfluoQTLthr is the number of fluorescence-QTLs at a LOD threshold of thr, and Nfluo is the number of high vs. low fluorescence comparisons (Nfluo = 48; the untagged experiment was excluded). At a significance threshold of LOD = 4.5, the estimated FDR was 7.3% (Figure 4—figure supplement 1), which we used to call significant QTLs. For some overlap analyses (see below), we used a threshold of LOD = 3.0, which corresponded to an FDR of 13%.
To call significant QTLs across replicates, we first scanned each replicate for QTLs at a permissive threshold of LOD > 3.0. Second, at each resulting QTL peak position, we averaged ΔAF and LOD scores across all available replicates without applying a LOD filter to each replicate. Third, we collapsed groups of overlapping QTLs, which we defined as QTLs whose peaks were within 75,000 bp of each other in the different replicates. For each group of these overlapping QTLs, we averaged the LOD scores, the ΔAFs, the peak positions, and the location confidence intervals to form one merged QTL. Of the resulting merged QTLs, we retained those that exceeded our stringent significance threshold of LOD ≥ 4.5.
To gauge reproducibility of these significant QTLs, we counted the number of replicates in which a given QTL had been detected at the permissive LOD >3.0, using the same definition of positional overlap as above. The majority (74%) of significant QTLs were shared across all the corresponding replicates. Two tagged genes had more than two replicates (GPD1 and UGP1). For these genes, requiring all replicates to be significant is conservative. Therefore, we also estimated the average reproducibility of all mRNA-QTLs or all protein-QTLs by calculating the average fraction of replicates that had a QTL at a given merged QTL:
Here, NshareQTLij is the number of replicates for which the QTL i is detected for the tagged gene j at LOD > 3, and Nrepj is the number of replicates performed for the tagged gene j. Note that if only a single replicate has a QTL at a given merged QTL, this fraction takes on a value of zero because in such a case, there is no overlap among replicates at this QTL. The observed fraction_overlap was 0.76 for the protein-QTLs and 0.78 for the mRNA-QTLs.
As a measure of effect size for a QTL, we used ΔAF. Briefly, in the presence of a QTL, there are two subpopulations of segregants that carry the two QTL alleles. The effect of the QTL is given by the difference in mean between these two subpopulations, assuming identical standard deviations for the two alleles. During FACS of the combined population with both alleles, cells that exceed the threshold for inclusion in the high population are enriched for the allele that increases expression, and vice versa for the low population. The ratio of cells that carry the two alleles in each extreme population determines the allele frequency at the locus in that population, and subtracting these allele frequencies for the two populations gives the ΔAF metric. For QTLs with larger effect sizes, the two subpopulations are further apart in phenotype. This increases the degree of genotype enrichment in the two extreme populations, as reflected in a larger ΔAF.
Comparison of mRNA-QTLs and protein-QTLs
Request a detailed protocolTo compare mRNA-QTLs and protein-QTLs of the same gene, we first considered all merged QTLs that exceeded a permissive threshold of LOD > 3.0 (after merging replicates as described above). We considered an mRNA-QTL and a protein-QTL for the same gene with overlapping confidence intervals as a QTL pair across mRNA and protein. We manually curated the result of this overlap analysis for six cases; after curation, QTLs on chromosomes XV (ARO8), VIII (MTD1), XIII (CYC1) and XIII (RPS10A) were considered to be pairs, and QTLs on chromosome V (GPD1) and XIV (MTD1) were considered to be mRNA-specific.
From this initial set, we retained those QTL pairs at which the given mRNA and/or protein QTL met a more stringent LOD score of >4.5 (FDR = 7.2%). Applying this higher threshold only after the more permissive overlap analysis allowed us to consider QTL pairs even if one of the paired QTLs did not pass the strong significance threshold of LOD >4.5. As an example, we considered overlapping QTLs on chromosome XI that affected OLE1 expression (mRNA-QTL LOD = 4.4, protein-QTL LOD = 15.5) to be a pair even if the mRNA LOD score was below the stringent significance threshold. In such cases, we deemed it more conservative to assume that the weaker QTL exists but narrowly failed to reach significance than to declare the stronger QTL as specific for mRNA or protein. We discarded all QTLs located between 350 and 550 kb on chromosome XIV, as this region may affect mCherry fluorescence independently of the tagged gene.
We distinguished four types of QTLs (Figure 5A and Figure 5—figure supplement 1). The shared QTL pairs either had similar effects on mRNA and protein abundance (16 QTL pairs, defined as having the same sign of ΔAF), or discordant effects on mRNA and protein (7 QTL pairs, different sign of ΔAF). All QTLs that were not part of a pair were considered to be specific (11 mRNA-specific QTLs, 52 protein-specific QTLs).
Finally, we conducted an analysis of mRNA or protein QTL effect sizes and directions that avoided having to define potentially paired QTLs as significant or not. For each mRNA-QTL (or protein-QTL), we extracted the ΔAF from the protein-QTL (or mRNA-QTL) data at the same exact position, irrespective of significance in the other data. We used these values to compute correlations of effects and to examine shared directionality of effects between mRNA-QTLs and protein-QTLs (Figure 5B, Figure 5—source data 2).
Allelic engineering for YAK1 fine-mapping
Request a detailed protocolTo obtain strains with scarless allelic swaps in haploids, we used a strategy based on double-cut CRISPR-swap (Lutz et al., 2019). We flanked each of the four tiles to be switched by two resistance markers (hphMX and KanMX) using our regular yeast transformation protocol (see above). The yeast were then transformed with 100 ng of CRISPR-Swap plasmid (pFA0055-gCASS5a, Addgene plasmid # 131774) and 1 µg of DNA repair template amplified either from BY, BY GPD1-GFP, or RM. The transformed cells were spread on SC -leucine plates, selecting for the presence of the plasmid expressing CAS9 and a gRNA targeting and cleaving a sequence present in both of the two resistance cassettes. We used strains in which GPD1 was tagged with GFP but not the gRNA tag, as the gRNA in our tag would likely have directed CAS9 to cleave the mCherry promoter. Cleavage of both cassettes resulted in the region in between the resistance cassettes to be replaced by the repair template. Transformed clones were screened for the double loss of antibiotic resistance to identify those with successful editing.
We introduced the 148,659 G→A variant, which we had detected through sequence analysis (see below), in YAK1 by single-cut CRISPR-swap (Lutz et al., 2019). We replaced the YAK1 sequence with a hphMX resistance cassette insertion to create yak1Δ::hphMX. We then delivered the CRISPR-Swap plasmid along with a repair template DNA produced by fusion PCR to carry either the G or A allele at the variant position (primers OFA0874 to OFA0881 in Supplementary file 3). Five clones of each allele, (YAK1wt and YAK1Q578*) were confirmed by Sanger sequencing (primers OFA0883 and OFA0882 in Supplementary file 3).
Sequence analyses to identify the YAK1 and other new DNA variants
Request a detailed protocolTo search for new variants in our populations that were not known to be present in the BY and RM strains, we used the sequencing reads from each selected segregant population. In each of 125 populations, we analyzed bam files after collapsing PCR duplicates. We applied samtools mpileup (--min-BQ 0) and bcftools call (-vc), either locally in the region affecting Gpd1-GFP (-r chrX:146000–150000) or on the whole genome to generate vcf files containing variant information. The vcf files were merged in R to generate matrices of polymorphic positions along with their allele frequency and coverage. The allele frequencies of the YAK1 polymorphisms were plotted along the genome for each population (Figure 6—figure supplement 1). We excluded all 47,754 previously known BY/RM variants from the whole-genome polymorphism matrices, and also removed variants with a bcftools quality score below 30. Among the 7624 remaining variants sites, 5822 were identified in only one or two populations and were deemed to be sequencing errors. Only one variant was shared across most (71 out of 75) of the populations created from strains from the GFP collection (strains tagged at ARO8, BMH2, GPD1, MTD1, and UGP1; Figure 5—source data 1) and absent in all other populations. This variant was the 148,659 G-to-A SNV in YAK1.
YAK1 genotyping
Request a detailed protocolThe region containing the YAK1Q578* variant was PCR-amplified from genomic DNA isolated from strains carrying FAA4-GFP and YMR315W-GFP (GFP collection), PGM1-TAP, NOT5-TAP, EMI2-TAP, TUB1-TAP (TAP-tag collection), BY4741 (ATCC 201388), YLK1879 (a BY strain from the Kruglyak lab) and YLK1950 (an RM strain from the Kruglyak lab) using primers OFA0883 and OFA0882 (Supplementary file 3). The resulting PCR product was Sanger sequenced using OFA0883 to genotype the YAK1 variant. The YAK1Q578* variant was observed only in the strains obtained from the GFP and TAP-tag collections (Ghaemmaghami et al., 2003; Huh et al., 2003).
Differential expression analysis by RNA-seq and mass spectrometry
Cell harvest
Request a detailed protocolWe quantified RNA and protein from five biological replicates (different clones obtained after CRISPR-Swap) of YAK1wt and YAK1Q578*. For each of these 10 strains, fresh colonies were used to inoculate 5 ml of SC medium (completed with uracil, arginine, histidine, and leucine) and the culture was grown overnight at 30°C. The next day, 50 ml of SC media in an Erlenmeyer flask was inoculated with the overnight culture to an initial OD600 of 0.05 and were grown under shaking at 30°C. When the OD600 reached 0.35–0.45 (late log phase), four aliquots of 7 ml of the culture was transferred to a falcon tube, centrifuged, and washed with 1 ml of sterile PBS buffer at 30°C (Phosphate Buffered Saline, pH 7.5). The pellets were immediately frozen in liquid nitrogen. For each strain, one cell pellet was used for RNA-seq and another for protein mass spectrometry.
RNA extraction and sequencing
Request a detailed protocolRNA extraction and library preparation were conducted as described in Lutz et al., 2019. Briefly, RNA extraction was done in two batches that each contained equal numbers of clones from the two groups (first batch: clones 1 and 2 of BY YAK1wt and clones 1 and 2 of BY YAK1Q578*, second batch: clones 3, 4, and 5 of BY YAK1wt and clones 3, 4, and 5 of BY YAK1Q578*). We used ZR Fungal/Bacterial RNA mini-prep kit with DNase I digestion (Zymo Research R2014), following the kit manual. RNA was eluted in 50 µl of RNase/DNase free water, quantified, and checked for integrity on an Agilent 2200 TapeStation. All RNA Integrity Numbers were higher than 9.5 and all concentrations were above 120 ng/µl. The RNA samples were stored at −20°C until use. Poly-A RNA selection was done using 550 ng of total RNA using NEBNext Poly(A) mRNA Magnetic Isolation Module (NEB E7490L), processing all samples in one batch. We prepared the library using NEB Ultra II Directional RNA library kit for Illumina (NEB E7760) used dual index primers (from NEBNext Multiplex Oligos for Illumina, NEB E7600S) for multiplexing, and amplified the library for ten cycles. The libraries were quantified using Qubit dsDNA HS Assay Kit (Thermo Fisher Scientific Q32854) and pooled at equal mass for sequencing using 75 bp single-end reads on Illumina NextSeq 550 at UMGC. The sequencing reads are available on NCBI GEO via BioProject PRJNA644804. Using the trimmomatic software (Bolger et al., 2014), reads were trimmed of adapters and low quality bases and filtered to be at least 36 bp long. Reads were then pseudo-aligned to the S. cerevisiae transcriptome (Ensembl build 93) and counted using kallisto (Bray et al., 2016). We used RSeQC to calculate Transcript Integrity Numbers (TIN) which provided an estimation of alignment quality for each gene of each sample. We excluded any gene with at least one read count of zero or at least one TIN of zero across the ten samples. After filtering, 5755 genes remained for analysis. Differential expression analysis was performed using the DESeq2 R package (Love et al., 2014), using the extraction batch information as covariate. DESeq2 provided, for each gene, the log2-fold-change (YAK1Q578* vs. YAK1wt) and the p-value adjusted for multiple-testing using the Benjamini-Hochberg method (Benjamini and Hochberg, 1995; Figure 7—source data 1).
Protein extraction and mass spectrometry
Request a detailed protocolProtein extraction and quantification using mass spectrometry was performed by the Center for Mass Spectrometry and Proteomics at the University of Minnesota. Briefly, cells from the pellets were lysed by sonication (30%, 7 s in Branson Digital Sonifier 250) in protein extraction buffer (7 M urea, 2 M thiourea, 0.4 M triethylammonium bicarbonate pH 8.5, 20% acetonitrile and 4 mM tris(2-carboxyethyl)phosphine). The proteins were extracted using pressure cycling (Barocycler Pressure Biosciences NEP2320, 60 cycles of 20 s at 20 kpsi and 10 s at 0 kpsi), purified by centrifugation in 8 mM iodoacetamide (10 min at 12000 rpm), and quantify by Bradford assay. For each sample, 40 µg of extracted proteins were fragmented by trypsin digestion (Promega Sequencing Grade Modified Trypsin V5111), and labeled by tandem mass tag isotopes (TMT10plex Isobaric Label Reagent Set, Thermo Scientific 90110). All the tagged samples were pooled. Peptides were separated first based on hydrophobicity through two consecutive liquid chromatographies, followed by separation based on mass per charge after ionization in the first mass spectrometry step. The second mass spectrometry step after high-energy collision-induced dissociation allowed for the identification of the peptide and the quantification of the TMTs. MS-MS analysis was conducted on an Orbitrap Fusion Tribrid instrument (Thermo Scientific). Database searches were performed using Proteome Discoverer software, and post-processing and differential expression analysis was done using Scaffold 4.9. Differential expression statistics per protein were computed on mean peptide abundances after inter-sample normalization. Normalized abundances for each of the 2590 detected proteins are provided in Figure 7—source data 3. We adjusted p-values provided by Scaffold using the Benjamini-Hochberg method.
Protein and mRNA data comparison
Request a detailed protocolTo compare the effect of the YAK1 variant on mRNA and protein abundance, we examined the log2 fold-change of the mRNA abundance (from DESeq2) and protein abundance (from Scaffold) for the 2577 genes present in both datasets. We considered genes that were differentially expressed in mRNA or protein (adjusted p-value<0.05) and that showed no evidence of difference in the other quantity (raw p-value>0.05), to be specifically affected at the mRNA or protein level (Figure 7A, Figure 7—source data 4). In various categories of differentially expressed genes, we looked for gene ontology enrichment using GOrilla (Eden et al., 2009) with the list of 2577 genes detected in both mass spectrometry and RNA-seq as the background set (Figure 7—source data 2).
Data availability
Raw DNA reads from bulk segregant mapping are available via the NCBI BioProject PRJNA644804. Transcriptome sequencing data is available at GEO under accession GSE155998. Source Data files are available for Figures 4, 5, and 7.
NCBI Gene Expression OmnibusID GSE155998. Differential RNA abundance between S. cerevisiae yeast BY4741 strains with YAK1_WT and YAK1_Q578*.
NCBI BioProjectID PRJNA644804. Simultaneous mRNA and protein quantification in single cells reveals post-transcriptional effects of genetic variation.
Genetic influences on translation in yeastPLOS Genetics 10:e1004692.https://doi.org/10.1371/journal.pgen.1004692
The role of regulatory variation in complex traits and diseaseNature Reviews Genetics 16:197–212.https://doi.org/10.1038/nrg3891
Controlling the false discovery rate: a practical and powerful approach to multiple testingJournal of the Royal Statistical Society: Series B 57:289–300.https://doi.org/10.1111/j.2517-6161.1995.tb02031.x
Trimmomatic: a flexible trimmer for illumina sequence dataBioinformatics 30:2114–2120.https://doi.org/10.1093/bioinformatics/btu170
Near-optimal probabilistic RNA-seq quantificationNature Biotechnology 34:525–527.https://doi.org/10.1038/nbt.3519
Pervasive phenotypic impact of a large nonrecombining introgressed region in yeastMol Biol Evol Msaa 101:msaa101.https://doi.org/10.1093/molbev/msaa101
mRNAs, proteins and the emerging principles of gene expression controlNature Reviews Genetics 21:630–644.https://doi.org/10.1038/s41576-020-0258-4
Tempo and mode of regulatory evolution in DrosophilaGenome Research 24:797–808.https://doi.org/10.1101/gr.163014.113
Quantitative trait loci underlying gene product variation: a novel perspective for analyzing regulation of genome expressionGenetics 137:289–301.
Quantitative trait loci mapped to single-nucleotide resolution in yeastNature Genetics 37:1333–1340.https://doi.org/10.1038/ng1674
High-resolution genetic mapping with pooled sequencingBMC Bioinformatics 13 Suppl 6:S8.https://doi.org/10.1186/1471-2105-13-S6-S8
Natural selection on Cis and trans regulation in yeastsGenome Research 20:826–836.https://doi.org/10.1101/gr.101576.109
Gene expression networks in the Drosophila genetic reference panelGenome Research 30:485–496.https://doi.org/10.1101/gr.257592.119
Tunable and multifunctional eukaryotic transcription factors based on CRISPR/CasACS Synthetic Biology 2:604–613.https://doi.org/10.1021/sb400081r
The molecular basis of phenotypic variation in yeastCurrent Opinion in Genetics & Development 23:672–677.https://doi.org/10.1016/j.gde.2013.10.005
Genetic basis of proteome variation in yeastNature Genetics 39:1369–1375.https://doi.org/10.1038/ng.2007.22
System-wide molecular evidence for phenotypic buffering in ArabidopsisNature Genetics 41:166–167.https://doi.org/10.1038/ng.308
RNA sequencing reveals the complex regulatory network in the maize kernelNature Communications 4:2832.https://doi.org/10.1038/ncomms3832
Self-processing of ribozyme-flanked RNAs into guide RNAs in vitro and in vivo for CRISPR-mediated genome editingJournal of Integrative Plant Biology 56:343–349.https://doi.org/10.1111/jipb.12152
Yeast responses to stresses associated with industrial brewery handlingFEMS Microbiology Reviews 31:535–569.https://doi.org/10.1111/j.1574-6976.2007.00076.x
Rare and common variants: twenty argumentsNature Reviews Genetics 13:135–145.https://doi.org/10.1038/nrg3118
Mapping Cis- and trans-regulatory effects across multiple tissues in twinsNature Genetics 44:1084–1089.https://doi.org/10.1038/ng.2394
Identification and validation of genetic variants that influence transcription factor and cell signaling protein levelsThe American Journal of Human Genetics 95:194–208.https://doi.org/10.1016/j.ajhg.2014.07.005
Comparison of drug transporter gene expression and functionality in Caco-2 cells from 10 different laboratoriesEuropean Journal of Pharmaceutical Sciences 35:383–396.https://doi.org/10.1016/j.ejps.2008.08.004
The sequence alignment/Map format and SAMtoolsBioinformatics 25:2078–2079.https://doi.org/10.1093/bioinformatics/btp352
CRISPR RNA-guided activation of endogenous human genesNature Methods 10:977–979.https://doi.org/10.1038/nmeth.2598
Posttranscriptional control of gene expression in yeastMicrobiology and Molecular Biology Reviews 62:1492–1553.https://doi.org/10.1128/MMBR.62.4.1492-1553.1998
Regulatory divergence in Drosophila revealed by mRNA-seqGenome Research 20:816–825.https://doi.org/10.1101/gr.102491.109
Single-copy gene detection using branched DNA (bDNA)The Journal of Histochemistry and Cytochemistry 49:603–611.https://doi.org/10.1177/002215540104900507
ClampFISH detects individual nucleic-acid molecules using click chemistry based amplificationNature Biotechnology 37:nbt.4286.https://doi.org/10.1038/nbt.4286
mRNA-Seq whole-transcriptome analysis of a single cellNature Methods 6:377–382.https://doi.org/10.1038/nmeth.1315
Book16 High-Throughput Strain Construction and Systematic Synthetic Lethal Screening in Saccharomy cescerevisiaeIn: Stansfield I, Stark M. J, editors. Methods in Microbiology, Yeast Gene Analysis. Academic Press. pp. 369–707.https://doi.org/10.1016/S0580-9517(06)36016-3
10 years of GWAS discovery: biology, function, and translationThe American Journal of Human Genetics 101:5–22.https://doi.org/10.1016/j.ajhg.2017.06.005
mRNA detection in budding yeast with single fluorophoresNucleic Acids Research 45:e141.https://doi.org/10.1093/nar/gkx568
Heritability and genomics of gene expression in peripheral bloodNature Genetics 46:430–437.https://doi.org/10.1038/ng.2951
Decision letter
Patricia J WittkoppSenior and Reviewing Editor; University of Michigan, United States
David GreshamReviewer; New York University, United States
In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.
Acceptance summary:
We are excited to publish this work because it presents a novel, creative strategy for studying variation in RNA and protein expression in parallel. It then uses this strategy to simultaneously quantify RNA and protein expression in genetically variable cells and examines the genetic basis for this expression variation. Consistent with prior work, they find different loci contributing to variation in RNA and protein expression, but this work greatly strengthens this observation because the new methodology removes many other possible sources of differences between the two.
Decision letter after peer review:
Thank you for submitting your article "Simultaneous quantification of mRNA and protein in single cells reveals post-transcriptional effects of genetic variation" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by Patricia Wittkopp as the Senior and Reviewing Editor. The following individual involved in review of your submission has agreed to reveal their identity: David Gresham (Reviewer #1).
The reviewers have discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
We would like to draw your attention to changes in our revision policy that we have made in response to COVID-19 (https://elifesciences.org/articles/57162). Specifically, we are asking editors to accept without delay manuscripts, like yours, that they judge can stand as eLife papers without additional data, even if they feel that they would make the manuscript stronger. Thus the revisions requested below only address clarity and presentation.
This paper presents a creative, novel strategy for quantifying expression of RNA and protein of the same genes in the same cells. While some concerns were raised about the mRNA level plateau that will need to be addressed, the methods and analyses were generally felt to be sound. Using this clever method, levels of mRNA and protein were used to simultaneously map QTL affecting mRNA and protein expression levels. The primary finding was that different QTL are affecting these two levels of gene expression, which is consistent with prior work. This consistency lead to a range of opinions about the novelty of the work.
In a revision, it will be critical to make more clear the new biological insight provided by the work, beyond the experimental design. Specifically, you should more clearly state what was learned about eQTLs and pQTLs with these new approaches that goes beyond what was already known. Clarifying expectations would also help: Do you expect to simultaneously have the same eQTLs and pQTLs? What are the molecular mechanisms responsible for different mRNA and pQTLs? Authors claim that is not statistical, however large effect trans-QTLs are conserved across e-pQTLs. Either way, how different trans-factors can interact with target genes and be responsible for expression and protein differences?
We have left the full set of reviewer comments below in this case because they contain detailed, specific, and complementary comments and reflect this range of opinions.
Reviewer #1:
In this paper, Brion et al. performed expression QTL mapping in a budding yeast cross using two methods for quantifying gene expression at the level of transcript abundance and protein abundance. Previous eQTL studies have been performed separately for transcript abundance and protein abundance using separate experimental procedures. These studies have found that there is a largely unique set of trans QTL that underlie interindividual variation in mRNA and protein abundance. However, these comparative analyses have a number of limitations that the authors outline. To rigorously address the extent to which these loci overlap the authors developed a method that enables simultaneous detection of protein and mRNA levels in single cells. Protein expression is detected using GFP fusions. mRNA expression is quantified using a novel, and clever, CRISPR-based system that uses self-cleaving ribozymes to generate a guide RNA that directs a dCas9 to drive expression of mCherry. Using this approach the authors performed bulk segregant mapping in a cross between two yeast strains to map loci that underlie differences in expression of mRNA and protein. Consistent with prior studies, they find little overlap between eQTL affecting the two levels of expression. They follow up on one locus, YAK1, in detail using RNAseq and mass spec.
This is an interesting study that presents a very clever method of quantifying mRNA and protein expression simultaneously. The experiments and analysis are well-performed. Although the result is consistent with previous studies, it greatly strengthens those observations and thus is an important addition to the field. Prior to publication the authors should address the following:
– The authors state that "most genetic variation in gene expression arises from trans-acting variants…" and cite several studies from model organisms that support this. However, my understanding is that this is not the case in human studies. The authors should expand on their explanation as to why they haven't been found in humans, beyond the fact that they can occur anywhere in the genome.
– The authors refer to the poor correlation between mRNA and protein levels, but should mention the issues identified with this comparison in Csardi et al., 2015 https://pubmed.ncbi.nlm.nih.gov/25950722/
– Figure 1 states the cells were grown in YNB + glutamate. Why was this media used for the study?
– Figure 3—figure supplement 3 aims to show the heritability of the phenotype of sorted cells. This would be more convincing if scatter plots were presented that show that after the sorting there is no correlation between the green and red signal, or if there is, an explanation as to why.
– It would be easier to read the paper if after the method has been introduced the text referred to mRNA and protein expression rather than the color of the reporter.
– The authors refer to the protein QTL as "post-transcriptional", but I think that term makes one think of mRNA stability, which is not what is measured. Perhaps "translational control" would be more correct and clearer.
– What is the rationale for showing δ-allele frequencies in Figure 4A? Showing LOD scores in this figure would make more sense to me for reporting QTL.
– It is not intuitive to me that δ-allele frequencies correspond to effect sizes and this should be more clearly explained.
– It is not clear to me if there is any compensation used in the FACS analysis. Is there any correlation between the green and red channels when controls containing only green or only red are used?
– Why is a loess regression used to correct for cell size? Wouldn't dividing by FSC be simpler and have the desired effect of normalizing for cell size?
Reviewer #2:
In this work, the authors used a novel approach to simultaneously measure mRNA and protein levels and identify trans-acting genetic variants affecting both or either of both. The authors present a fast and accurate novel strategy, taking advantage of a CRISPR strategy, and FACS analysis in a highly used yeast cross. Unfortunately, this work does not provide novel genetic, neither biological insights and conclusions are similar to previous findings. I would recommend the authors try to extend their work to a higher number of genes/environments and obtained novel biological insights that could explain how trans-acting factors, including post-transcriptional and post-translation modifications, impact mRNA, protein and traits. Furthermore, a time course would be interesting to analyse and determine how trans-acting variants affect mRNA and protein levels over time and whether they correlate or not.
Reviewer #3:
In this manuscript, Brion Lutz and Albert detail the simultaneous mapping of trans-acting QTLs influencing mRNA production and steady-state protein level differences between yeast strains BY4741 and RM-11 for ten genes. While much work has previously been done evaluating mRNA-QTLs in this system and others, the simultaneous comparison of transcription and protein QTLs, even for a relatively small set of genes, is significant and timely. The authors used a clever experimental system that converts transcription rates of test genes into mCherry expression by CRISPR activation. Notably, the authors find most QTLs affect protein levels, but not mRNA levels, and that mRNA-QTLs do not always have the same effects on protein levels. They also report genome-wide comparisons of a YAK1 kinase allele-swap on steady-state protein and mRNA levels. Overall, the manuscript describes important results that suggest more study of pQTLs is needed to better model the impact of transcription regulatory variation. While generally solid, the authors should further evaluate the effects of the mCherry plateau and examine protein-protein interactions for each of the ten genes to evaluate how these might contribute to the disparities between RNA-QTLs and protein-QTLs. They are also missing important citations that should be included for scholarly accuracy. Numbered issues (chronological):
1) Introduction: the authors should include Emerson et al., 2010 for yeast, and references of similar work in Drosophila, including McManus et al., 2010; Coolon et al., 2014 and Huang et al., 2015 and potentially others.
2) Introduction paragraph three: Cenik et al., 2015 also used ribosome profiling and mass spectrometry to examine post-transcriptional variation amont human LCLs and should be cited.
3) One caveat to the experimental approach is the plateau in mCherry measurements at high transcription rates. The authors argue that this plateau does not affect most genes, by comparing qPCR signal from inducible GFP and ACT1. However, this seems problematic because the PCR efficiency of GFP and ACT1 may not be identical. Thus it's possible the mRNA plateau occurs at a lower level. If this were the case, one would expect to see fewer mRNA-QTLs for highly expressed genes than for low-expression genes. Is there any relationship between expression levels (RNA-seq) and number of mRNA-QTLs among the 10 genes tested? If high mRNA genes are artifactually missing mRNA-QTLs, it's possible more pQTLs also affect mRNA levels.
4) A related complication of the assay system is that the tested genes do not have native 3' UTRs or cleavage and polyadenylation sites (CPS). The "strength" of CPS, has an important impact on mRNA levels and may influence both mRNA stability and transcription (See Shalem et al. PMID: 25875337). The 3' UTR issue should be discussed.
5) "A majority of the loci corresponded to protein-QTLs that did not overlap an mRNA-QTL". Were these non-matching QTLs more common in genes with high mRNA levels (i.e. could this reflect the mCherry plateau)?
6) Figure 5 shows variation among the number of pQTLs and mRNA-QTLs for the 10 studied genes. In the Discussion, the authors note that some p-QTLs could result by changes in protein complex members. Is there any relationship between the number of pQTLs for genes and the number of protein-protein interactions each gene has?
7) The YAK1 mutation was identified as altering protein (but not mRNA) levels of ARO8, BMH2 and GPD1. The genome-wide RNA-seq and Mass-spec comparison of YAK1 alleles appears to only show this affect for GPD1. Were the other genes also significantly altered in mRNA or protein levels in this orthogonal experiment?
8) The paper ends with the suggestion that protein abundance is under more complex genetic control than mRNA abundance. While this seems very likely in general (before the results presented in this paper) due to protein turnover and post-translational modifications, I think the authors should reiterate at the end of the paper that this is for trans-QTLs specifically. I expect that cis-QTLs have more consistent influences on mRNA and protein levels.
https://doi.org/10.7554/eLife.60645.sa1Author response
This paper presents a creative, novel strategy for quantifying expression of RNA and protein of the same genes in the same cells. While some concerns were raised about the mRNA level plateau that will need to be addressed, the methods and analyses were generally felt to be sound. Using this clever method, levels of mRNA and protein were used to simultaneously map QTL affecting mRNA and protein expression levels. The primary finding was that different QTL are affecting these two levels of gene expression, which is consistent with prior work. This consistency lead to a range of opinions about the novelty of the work.
In a revision, it will be critical to make more clear the new biological insight provided by the work, beyond the experimental design. Specifically, you should more clearly state what was learned about eQTLs and pQTLs with these new approaches that goes beyond what was already known. Clarifying expectations would also help: Do you expect to simultaneously have the same eQTLs and pQTLs? What are the molecular mechanisms responsible for different mRNA and pQTLs? Authors claim that is not statistical, however large effect trans-QTLs are conserved across e-pQTLs. Either way, how different trans-factors can interact with target genes and be responsible for expression and protein differences?
1) The plateau in mRNA levels has been addressed as follows:
– As requested by reviewer #2, we examined the robustness of our chosen threshold of half the abundance of ACT1 mRNA by comparing the abundance of ACT1 we had used to that in another RNA-Seq dataset. In both datasets, the expression of the vast majority of genes are below this threshold, which is consistent with our prior interpretation of these data in the paper. We have also reworded the manuscript to make clear that we view this "threshold" as a broad guideline rather than a hard, numerically precise threshold.
– As requested by reviewer #3, we asked whether mRNA abundance as measured by external RNA-Seq data correlates with the number of mRNA-QTLs or protein-specific QTLs we detected. We found no such relationship, suggesting that the mCherry plateau does not greatly influence our results.
2) Our manuscript presents several novel advances:
First, and most obviously, we introduce a novel method. We anticipate that the reporter system we describe will find uses well beyond QTL mapping. To name just one example, at conferences where we have presented this work, we had discussions with bioengineers who were considering our method for real-time monitoring of the production of compounds or peptides in yeast bioreactors, in cases where large protein tags such as GFP would interfere with compound production. While detailed discussion of these possibilities is out of the scope of the current manuscript, we are confident that our method will be useful in many applications.
Regarding a perceived lack of biological insight, we believe that this was largely an issue with our previous presentation of the background literature. We have now addressed this with a completely rewritten Introduction (see point 3 below) to highlight the lack of consensus in the field on the relationship between eQTLs and pQTLs. We started this study believing this lack of consensus is mostly due to technical limitations of prior work (small sample size and experimental variability). We explicitly designed our system to be as free from these issues as possible. Thus, a possible outcome of our work was that in our well-controlled system, most discrepancies between eQTLs and pQTLs would disappear. The fact that these discrepancies did not disappear is an important and novel result.
Finally, we want to stress that our manuscript includes fine-mapping of a trans-acting pQTL to the causal variant. Fine-mapping causal variants is not a trivial task, making this fine-mapped variant a novel result in its own right.
3) To help readers appreciate the novelty of our results, we have followed the editor's suggestion to clarify expectations regarding the discrepancies between eQTLs and pQTLs. To do so, we rewrote the Introduction entirely. It now includes an expanded review of the eQTL/pQTL literature in multiple species (human, mouse, plants, and yeast). The new Introduction details discrepancies between the results reported by different studies. Briefly, different studies report different degrees to which a typical eQTL is also a pQTL for the same gene. There is also a range of reports on what fraction of pQTLs arises from protein-specific processes. We see no consensus in the literature on these questions, and our Introduction now makes this point more forcefully.
We go on to describe in more detail than before exactly how the technical shortcomings of earlier work (i.e., small sample size and / or comparisons across studies) may have inflated discrepancies between eQTLs and pQTLs. We believe that this rewritten Introduction makes it clear why the work we describe here was necessary, and why our results matter: although our system controlled for environmental factors and had high statistical power, we still found many discrepancies between eQTLs and pQTLs.
4) Regarding the comment "Authors claim that [it] is not statistical, however large effect trans-QTLs are conserved across e-pQTLs.": This is a fair point, and we discuss this below in more detail in our response to reviewer #2. To briefly restate here, we cannot formally rule out that incomplete power to detect mRNA-QTLs may have contributed to some apparent protein-specific QTLs, just like we cannot rule out the inverse scenario for mRNA-specific QTLs. Our manuscript explicitly acknowledges this possibility, includes analyses that aim to address limited statistical power, and provides evidence that while incomplete power undoubtedly leads to some false calls of specific QTLs, it is unlikely to account for all of them.
5) Regarding the comments "the molecular mechanisms responsible for different mRNA and pQTLs" and "how different trans-factors can interact with target genes and be responsible for expression and protein differences?". To our knowledge, the only protein-specific trans-pQTLs that have been reported in some detail are effects on protein complex stoichiometry (although we are not aware of experimental confirmation of these causal relationships). This mechanism is discussed in the revised Discussion section. We agree that experimental dissection of additional protein-specific trans-pQTLs would help to elucidate the mechanistic nature of eQTL and pQTL; however, we believe this is out the scope of the present paper, given the considerable challenges in identifying the causal genes and variants in QTL regions.
Reviewer #1:
This is an interesting study that presents a very clever method of quantifying mRNA and protein expression simultaneously. The experiments and analysis are well-performed. Although the result is consistent with previous studies, it greatly strengthens those observations and thus is an important addition to the field. Prior to publication the authors should address the following:
– The authors state that "most genetic variation in gene expression arises from trans-acting variants…" and cite several studies from model organisms that support this. However, my understanding is that this is not the case in human studies. The authors should expand on their explanation as to why they haven't been found in humans, beyond the fact that they can occur anywhere in the genome.
It is indeed true that individual trans-eQTLs (let alone causal variants in these eQTL regions) are harder to find in human populations than in model organisms. However, as a group, trans-eQTLs contribute more genetic variation to gene expression than cis-eQTLs in humans, just like in model organisms. These human results are not based on the effects of individual trans-eQTLs, but on genome-wide aggregate estimates from mixed linear models. These findings are described in two studies of human populations Grundberg et al., and Wright et al., which we cite in our Introduction.
To clarify these points, we have reworded this section as part of our extensively rewritten Introduction.
– The authors refer to the poor correlation between mRNA and protein levels, but should mention the issues identified with this comparison in Csardi et al., 2015 https://pubmed.ncbi.nlm.nih.gov/25950722/
We have added this reference as part of the rewritten Introduction.
– Figure 1 states the cells were grown in YNB + glutamate. Why was this media used for the study?
This medium was used only in the experiments in Figure 1. We had conducted these experiments in the presence of G418 (which requires the use of YNB + glutamate) because, at the time, we were concerned about potential loss of the reporter construct from the genome due to recombination between its flanking sequences. This concern later turned out to be unfounded. In the remainder of this study, we used synthetic complete medium, as indicated in the Materials and methods and the figure legends throughout the paper. We have also added more explanation for our choice of media to the Materials and methods section.
– Figure 3—figure supplement 3 aims to show the heritability of the phenotype of sorted cells. This would be more convincing if scatter plots were presented that show that after the sorting there is no correlation between the green and red signal, or if there is, an explanation as to why.
This is an interesting point. First, we'd like to emphasize that these experiments, in which we measured the response to phenotypic selection after a few generations of growth, are informed by the breeder's equation in quantitative genetics. This equation relates the strength of selection that is applied to a population (defined as the difference in mean phenotype between selected and unselected individuals) to the phenotypic response to selection (defined as the shift in mean of the population after one generation). The breeder's equation is usually applied to a single generation of sexually reproducing diploid organisms, making its use in our haploid segregants (which undergo multiple mitotic cell divisions but no meiosis after FACS) imperfect. Nevertheless, it provides a framework for interpreting the shifts in mean phenotypes that we reported.
Second, our selection scheme would not necessarily be expected to result in a strong reduction in correlation between mRNA and protein. Specifically, the effects of each individual locus are small compared to the nongenetic cell-to-cell noise in gene expression, which can reach up to 100-fold variation even in isogenic yeast cell populations. This means that our selections are most likely not strong enough to result in pure populations, in which every cell carries all alleles that increase (or decrease) protein or mRNA at all or even most loci across the genome. Instead, we fully expect that our sorted extreme populations still contain alleles for high and low expression at many QTLs. This expectation is confirmed by the fact that complete fixation of opposite alleles in high vs low populations was almost never observed in this work, consistent with our previous bulk-segregant studies.
Therefore, each of our extreme populations still segregates at multiple QTLs. Because a portion of these QTLs do influence both mRNA and protein, these loci are expected to create a correlation between mRNA and protein even after selection. To confirm this, we computed the correlations between red and green fluorescence in each population, and compared their magnitude before and after sorting:
In Author response image 1, the y-axis shows the difference between the mCherry/GFP correlation coefficients measured after FACS and growth of an unsorted population ("all_pop") or the fluorescence-based sorted population ("sorted_pop"). The numbered populations are: 1=low mCherry, 2=high mCherry, 3=low GFP, 4=high GFP. There is one subplot per gene, with lines connecting the four populations from a given biological replicate. The "NA" population refers to untagged cells (as explained in Figure 3—figure supplement 3).
In most populations, there is not a consistent difference in correlation coefficients between cells during FACS and after growth. The only consistent exception is the gene CYC1, where the post-growth cells do show reduced correlation. This gene is affected by a strong QTL, such that FACS can create relatively "pure" populations of cells that all share one of the alleles at this locus. Note that even here, the other QTLs affecting this gene do still segregate and create correlations among single cells. We hope that these explanations address this comment. We have chosen to keep our presentation of the results in Figure 3—figure supplement 3 unchanged.
– It would be easier to read the paper if after the method has been introduced the text referred to mRNA and protein expression rather than the color of the reporter.
This is a great point, and one that we had discussed extensively while preparing our manuscript. As the reviewer will appreciate, there is a tradeoff here between slightly greater ease of reading and more exact language. Specifically, while we certainly assume (and believe to have shown) that "red" mCherry corresponds to mRNA and "green" GFP to protein, we ultimately can only say with certainty that we measure fluorescence signals during flow cytometry. Therefore, when describing flow cytometry, we had opted to refer to the observed quantities (i.e. fluorescence) rather than the molecules they report on. For example, we do feel it is more appropriate to say that sorting was performed based on "mCherry", rather than, say "mRNA". However, note that as the text progresses to presenting the QTL results, it already has a transition in terminology to "mRNA-QTLs" and "protein-QTLs", as opposed to "mCherry-QTLs" and "GFP-QTLs".
We have reworded some instances of where this terminology may have been confusing (e.g. Figure 4 legend).
– The authors refer to the protein QTL as "post-transcriptional", but I think that term makes one think of mRNA stability, which is not what is measured. Perhaps "translational control" would be more correct and clearer.
Our use of "post-transcriptional" was intended to be understood very literally, to refer to any process that may affect protein abundance after transcription. As the reviewer correctly notes, because our mCherry system stops reporting on the fate of the mRNA after the gRNA has been released, it indeed cannot report on mRNA degradation (other than for some mRNAs that are degraded before gRNA release is complete). However, GFP-tag signal is sensitive to mRNA degradation, to the extent that mRNA degradation limits the time available for translation and therefore the number of protein molecules that can be generated from a given mRNA molecule. Further, we feel that "translational control" is too specific, since it omits possible influences of protein degradation on protein abundance. For these reasons, we think that the term "post-transcriptional" is appropriately specific [because it only excludes the one process (transcription) that will have the strongest influence on mCherry signal] but also sufficiently broad (because it captures any mechanism that can shape protein abundance downstream of transcription). Therefore, we have decided to keep this term in the paper. We have added a brief statement to the Discussion that makes clear that we take "post-transcriptional" to mean "mRNA stability, translation, or protein degradation"
– What is the rationale for showing δ-allele frequencies in Figure 4A? Showing LOD scores in this figure would make more sense to me for reporting QTL.
The rationale for showing δ-AF is four-fold. First, LOD scores would not show the direction of effect at a given QTL (i.e., whether the BY or the RM allele increases expression), while δ-AF does show this. Second, the relationship between δ-AF and LOD scores is not linear: as δ-AF approaches fixation, the LOD score rises quickly. From a visualization standpoint, this creates the issue that strong QTLs, if shown at full scale, "blow out" the display, such that weaker loci become very hard to see. These two points are the main reason we chose to use δ-AF in this figure.
Third, LOD scores estimated by the "Multipool" algorithm we use are sensitive to the depth of sequencing coverage. At higher coverage, LOD scores increase for a given effects size, because deeper sequencing provides more certainty to the model that the allele frequencies in the two populations are not the same. We acknowledge that this is not a major concern for our present work, because sequencing coverage is fairly similar across our samples. Fourth, the δ-AFs are a more "raw" display of the data than LOD scores; these allele frequencies are directly derived from the actual measurements (i.e. allele counts in Illumina data). Together, these points led us to use δ-AF as the more appropriate choice for displaying the data visually.
Note that our supplementary tables do of course report LOD scores, along with many other characteristics of each QTL, such that interested readers have full access to this information.
– It is not intuitive to me that δ-allele frequencies correspond to effect sizes and this should be more clearly explained.
We illustrate the relationship between effect size and δ AF in Author response image 2, which is based on simulated bulk segregant mapping experiments. We have also added a brief explanation of these considerations to the Materials and methods, at the end of the "QTL mapping" section.
Consider a population of cells made up of two subpopulations that each carry one of two alleles at a QTL that affects the given trait (left panel in Author response image 2). The effect of the QTL is given by the difference in mean between the two subpopulations. When we use FACS on the combined population, we use the indicated thresholds for picking "high" vs "low" cells. Note that in the high population defined by the "high" threshold, the fraction of cells with the BY allele is lower than the cells with the RM allele, and vice versa for the low population. The ratio of BY to RM alleles in the extreme populations determines the allele frequencies, and subtracting these from each other gives the δ AF metric that we use as a measure of effect size.
If we increase the effect size of the QTL, the two subpopulations are being pushed further apart, and the degree of genotype enrichment in the two extreme populations increases. This increases δ AF, as shown on the right for a range of simulated effect sizes. Note that the relationship between effect size and δ AF is not linear, especially for strong QTLs, but that δ AF rises monotonically as a function of effect size.
– It is not clear to me if there is any compensation used in the FACS analysis. Is there any correlation between the green and red channels when controls containing only green or only red are used?
We did not use compensation during FACS. This information has been added to the Materials and methods. During our pilot studies for this work, we did not notice interference between the green and red fluorescence signals.
– Why is a loess regression used to correct for cell size? Wouldn't dividing by FSC be simpler and have the desired effect of normalizing for cell size?
We chose LOESS regression over simple division, because LOESS allows for relationships between FSC and fluorescence that are not strictly linear with a constant slope across the entire range of the data. In practice, there is almost no difference between correlations between red and green fluorescence using LOESS vs division, see Author response image 3.
Reviewer #2:
In this work, the authors used a novel approach to simultaneously measure mRNA and protein levels and identify trans-acting genetic variants affecting both or either of both. The authors present a fast and accurate novel strategy, taking advantage of a CRISPR strategy, and FACS analysis in a highly used yeast cross. Unfortunately, this work does not provide novel genetic, neither biological insights and conclusions are similar to previous findings. I would recommend the authors try to extend their work to a higher number of genes/environments and obtained novel biological insights that could explain how trans-acting factors, including post-transcriptional and post-translation modifications, impact mRNA, protein and traits. Furthermore, a time course would be interesting to analyse and determine how trans-acting variants affect mRNA and protein levels over time and whether they correlate or not.
We thank the reviewer for their positive assessment of our method. While we look forward to applying our reporter system to additional proteins and questions in future work, we consider the scale of experiments suggested by the reviewer to be out of the scope of the present manuscript. In addition, we politely disagree that our work has not provided novel insights, but do acknowledge that our presentation of the background literature for this work may have obscured some novelty, as also pointed out by the editor.
To address this, we have worked to improve our presentation of the background literature as well as our new insights in the paper. Specifically, we have rewritten the entire Introduction to provide a more detailed review of existing comparisons between eQTLs and pQTLs from several species (human, mouse, plants, and, of course, yeast). In our opinion, the literature is far from settled on how well eQTLs and pQTLs agree. For example, while some studies have shown that most cis-acting eQTLs are also pQTLs, others have reported systematic buffering of the effects of cis-eQTLs such that they have smaller, if any, effects on the proteins of the same gene. Likewise, for trans-eQTLs, some studies report nearly complete lack of overlap (indeed, Foss et al., 2011 even pointed out that trans-QTLs may show less overlap than expected by chance), while others focus on relatively strong correlations of effect sizes to indicate that most trans-acting variants have consistent effects. Whether or not eQTL/pQTL discrepancies exist is an important and open question. Among other implications for the biology of gene expression, this question raises the practical issue of how informative existing studies of mRNA expression variation are, including many large-scale projects in human genetics, if mRNA expression changes do not correlate with changes in protein expression.
Our revised Introduction now stresses more explicitly technical considerations in earlier work related to 1) low statistical power and 2) experimental differences between studies across which eQTLs and pQTL are almost always compared. These two issues have the potential to create the impression of discrepant eQTLs and pQTLs, even if eQTLs and pQTLs were consistent in reality.
Our new approach here was explicitly designed to directly address both of these issues. We aimed to create a system that maximizes the chance of detecting consistency between eQTLs and pQTLs. The fact that we nevertheless detected considerable discrepancies between eQTLs and pQTLs is an important result. In addition to our rewritten Introduction, we have also emphasized these points more in the Discussion section.
On a more technical level, we are convinced that our reporter system could find uses outside of mapping regulatory genetic variation. Finally, we stress that fine-mapping QTLs to the causal variant remains a difficult challenge. That we were able to accomplish it here for the YAK1 gene adds further novelty to our study.
Reviewer #3:
In this manuscript, Brion Lutz and Albert detail the simultaneous mapping of trans-acting QTLs influencing mRNA production and steady-state protein level differences between yeast strains BY4741 and RM-11 for ten genes. While much work has previously been done evaluating mRNA-QTLs in this system and others, the simultaneous comparison of transcription and protein QTLs, even for a relatively small set of genes, is significant and timely. The authors used a clever experimental system that converts transcription rates of test genes into mCherry expression by CRISPR activation. Notably, the authors find most QTLs affect protein levels, but not mRNA levels, and that mRNA-QTLs do not always have the same effects on protein levels. They also report genome-wide comparisons of a YAK1 kinase allele-swap on steady-state protein and mRNA levels. Overall, the manuscript describes important results that suggest more study of pQTLs is needed to better model the impact of transcription regulatory variation. While generally solid, the authors should further evaluate the effects of the mCherry plateau and examine protein-protein interactions for each of the ten genes to evaluate how these might contribute to the disparities between RNA-QTLs and protein-QTLs. They are also missing important citations that should be included for scholarly accuracy. Numbered issues (chronological):
1) Introduction: the authors should include Emerson et al., 2010 for yeast, and references of similar work in Drosophila, including McManus et al., 2010; Coolon et al., 2014 and Huang et al., 2015 and potentially others.
We thank the reviewer for these suggestions. As part of our revised Introduction, we have added the requested references, along with others. Some of the references suggested by the reviewer examined F1 hybrids and parent strains, but did not perform linkage analysis in segregants (which was the intended topic of this sentence) and in some cases examined divergent species rather than strains of the same species. For these reasons, we had not cited these papers before. To accommodate these references in the revision, we have changed the sentence in the manuscript as follows:
"By contrast, experimental crosses in organisms such as yeast (Albert et al., 2018; Brem et al., 2002; Brion et al., 2020; Clément-Ziza et al., 2014; Emerson et al., 2010; Thompson and Cubillos, 2017), plants (Fu et al., 2013; West et al., 2007; Zhang et al., 2011), worms (Snoek et al., 2017; Viñuela et al., 2010), flies (Coolon et al., 2014; Everett et al., 2020; Huang et al., 2015; McManus et al., 2010), and mouse (Gerrits et al., 2009; Hasin-Brumshtein et al., 2016) have quantified the global contribution of cis and trans-acting variation and identified thousands of trans-eQTLs."
2) Introduction paragraph three: Cenik et al., 2015 also used ribosome profiling and mass spectrometry to examine post-transcriptional variation amont human LCLs and should be cited.
We have added this reference, as well as Wu et al., 2013 as an early application for mass spectrometry for pQTL mapping in human cells.
3) One caveat to the experimental approach is the plateau in mCherry measurements at high transcription rates. The authors argue that this plateau does not affect most genes, by comparing qPCR signal from inducible GFP and ACT1. However, this seems problematic because the PCR efficiency of GFP and ACT1 may not be identical. Thus it's possible the mRNA plateau occurs at a lower level. If this were the case, one would expect to see fewer mRNA-QTLs for highly expressed genes than for low-expression genes. Is there any relationship between expression levels (RNA-seq) and number of mRNA-QTLs among the 10 genes tested? If high mRNA genes are artifactually missing mRNA-QTLs, it's possible more pQTLs also affect mRNA levels.
We have conducted the analysis reviewer suggested, and found no relationship between a gene's mRNA expression in the Albert, Bloom et al., 2018 data and the number of mRNA-QTLs (or any other type of QTL). Specifically, while the correlation between the number mRNA-QTLs and mRNA abundance did have a negative coefficient (rho = -0.15), the correlation was not significant (p = 0.67). Thus, this analysis suggests that the mCherry plateau does not cause systematic underestimation of the number of mRNA-QTLs in our data.
However, we agree with the reviewer's point that the comparison of ACT1 qPCR and RNA-seq is not without problems, as also noted by reviewer 2 above. This is true even though we did correct for different PCR efficiencies using the calibration curves shown in Figure 1—figure supplement 2.
In the revised manuscript, we explicitly acknowledge this issue and have softened our language describing how many genes can likely be analyzed with our reporter system.
4) A related complication of the assay system is that the tested genes do not have native 3' UTRs or cleavage and polyadenylation sites (CPS). The "strength" of CPS, has an important impact on mRNA levels and may influence both mRNA stability and transcription (See Shalem et al. PMID: 25875337). The 3' UTR issue should be discussed.
We agree, and have added a discussion of the possible consequence of 3' UTR removal to the Discussion section:
"Like other reporters that tag the 3' end of genes, our system decouples the native 3' UTR from the gene of interest. 3' UTRs can influence gene expression, and genetic trans-effects that involve mechanisms that target the 3' UTR would be missed by our system. However, loss of the native 3' UTR is not expected to inflate differences between mRNA-QTLs and protein-QTLs. Instead, 3' UTR removal may reduce differences between mRNA-QTLs and protein-QTLs because it eliminates an element of the mRNA that could be a target for mRNA-specific processes"
5) "A majority of the loci corresponded to protein-QTLs that did not overlap an mRNA-QTL". Were these non-matching QTLs more common in genes with high mRNA levels (i.e. could this reflect the mCherry plateau)?
Across the ten genes, the correlation between the number of protein-specific QTLs and mRNA abundance had a rho value of -0.02 with a p-value of 0.96. We do not think that the mCherry plateau is responsible for these protein-specific QTLs.
6) Figure 5 shows variation among the number of pQTLs and mRNA-QTLs for the 10 studied genes. In the Discussion, the authors note that some p-QTLs could result by changes in protein complex members. Is there any relationship between the number of pQTLs for genes and the number of protein-protein interactions each gene has?
To address this question, we extracted physical protein/protein interactions from the bioGRID database (https://thebiogrid.org/). Genes with more of these interactions tended to have a higher number of protein-QTLs (rho = 0.51) as well as protein-specific QTLs (rho = 0.55) but neither of these results was statistically significant (p-values: 0.13 and 0.10, respectively). We suspect that ten genes is likely not enough to test this question rigorously. Future analyses with more genes may well strengthen the significance of this relationship. We thank the reviewer for this excellent suggestion, which we will keep in mind for future work.
7) The YAK1 mutation was identified as altering protein (but not mRNA) levels of ARO8, BMH2 and GPD1. The genome-wide RNA-seq and Mass-spec comparison of YAK1 alleles appears to only show this affect for GPD1. Were the other genes also significantly altered in mRNA or protein levels in this orthogonal experiment?
The results are shown in Author response table 1, and can also be found in the source data for Figure 5—figure supplement 2 and Figure 7—figure supplement 3:
As the reviewer correctly spotted, ARO8 and BMH2 were not significantly different in either RNA-seq or mass spectrometry. We believe that this is likely due to the small effect of the YAK1 variant on the protein levels of these two genes, which is much smaller than its effect on Gpd1. In Author response table 1, note the big differences among these genes in terms of their deltaAF metric (which is related to the effect size of a QTL; see our response to reviewer 1 above) as well as their LOD scores, which show that GPD1 has a much stronger protein-QTL than the other two genes. Our mass-spectrometry dataset with n = 5 in each group is likely not powered to detect the small effects at the other two genes.
8) The paper ends with the suggestion that protein abundance is under more complex genetic control than mRNA abundance. While this seems very likely in general (before the results presented in this paper) due to protein turnover and post-translational modifications, I think the authors should reiterate at the end of the paper that this is for trans-QTLs specifically. I expect that cis-QTLs have more consistent influences on mRNA and protein levels.
We have clarified in the final paragraph of the paper that our results exclusively cover trans-QTLs.
https://doi.org/10.7554/eLife.60645.sa2Article and author information
Author details
National Institute of General Medical Sciences (R35-GM124676)
- Frank Wolfgang Albert
Alfred P. Sloan Foundation (FG-2018- 10408)
- Frank Wolfgang Albert
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
We thank Scott McIsaac and Leonid Kruglyak for yeast strains, Joshua Bloom for implementing the MULTIPOOL algorithm in R, and Mahlon Collins for the mCherry/GFP protein fusion flow cytometry data. We thank the University of Minnesota Genomics Center, the University of Minnesota's Center for Mass Spectrometry and Proteomics, and the University of Minnesota's Flow Cytometry Resource. We thank members of the Albert lab for critical feedback on the manuscript.
Senior and Reviewing Editor
- Patricia J Wittkopp, University of Michigan, United States
- David Gresham, New York University, United States
Version history
- Received: July 2, 2020
- Accepted: November 14, 2020
- Accepted Manuscript published: November 16, 2020 (version 1)
- Version of Record published: December 1, 2020 (version 2)
© 2020, Brion et al.
This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.
- 8,230
- Page views
- 554
- Downloads
- 23
- Citations
Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.
Download links
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Further reading
- Evolutionary Biology
- Genetics and Genomics
Maintaining germline genome integrity is essential and enormously complex. Although many proteins are involved in DNA replication, proofreading, and repair, mutator alleles have largely eluded detection in mammals. DNA replication and repair proteins often recognize sequence motifs or excise lesions at specific nucleotides. Thus, we might expect that the spectrum of de novo mutations – the frequencies of C>T, A>G, etc. – will differ between genomes that harbor either a mutator or wild-type allele. Previously, we used quantitative trait locus mapping to discover candidate mutator alleles in the DNA repair gene Mutyh that increased the C>A germline mutation rate in a family of inbred mice known as the BXDs (Sasani et al., 2022, Ashbrook et al., 2021). In this study we developed a new method to detect alleles associated with mutation spectrum variation and applied it to mutation data from the BXDs. We discovered an additional C>A mutator locus on chromosome 6 that overlaps Ogg1, a DNA glycosylase involved in the same base-excision repair network as Mutyh (David et al., 2007). Its effect depends on the presence of a mutator allele near Mutyh, and BXDs with mutator alleles at both loci have greater numbers of C>A mutations than those with mutator alleles at either locus alone. Our new methods for analyzing mutation spectra reveal evidence of epistasis between germline mutator alleles and may be applicable to mutation data from humans and other model organisms.
- Evolutionary Biology
- Genetics and Genomics
Although gene expression divergence has long been postulated to be the primary driver of human evolution, identifying the genes and genetic variants underlying uniquely human traits has proven to be quite challenging. Theory suggests that cell-type-specific cis-regulatory variants may fuel evolutionary adaptation due to the specificity of their effects. These variants can precisely tune the expression of a single gene in a single cell-type, avoiding the potentially deleterious consequences of trans-acting changes and non-cell type-specific changes that can impact many genes and cell types, respectively. It has recently become possible to quantify human-specific cis-acting regulatory divergence by measuring allele-specific expression in human-chimpanzee hybrid cells—the product of fusing induced pluripotent stem (iPS) cells of each species in vitro. However, these cis-regulatory changes have only been explored in a limited number of cell types. Here, we quantify human-chimpanzee cis-regulatory divergence in gene expression and chromatin accessibility across six cell types, enabling the identification of highly cell-type-specific cis-regulatory changes. We find that cell-type-specific genes and regulatory elements evolve faster than those shared across cell types, suggesting an important role for genes with cell-type-specific expression in human evolution. Furthermore, we identify several instances of lineage-specific natural selection that may have played key roles in specific cell types, such as coordinated changes in the cis-regulation of dozens of genes involved in neuronal firing in motor neurons. Finally, using novel metrics and a machine learning model, we identify genetic variants that likely alter chromatin accessibility and transcription factor binding, leading to neuron-specific changes in the expression of the neurodevelopmentally important genes FABP7 and GAD1. Overall, our results demonstrate that integrative analysis of cis-regulatory divergence in chromatin accessibility and gene expression across cell types is a promising approach to identify the specific genes and genetic variants that make us human.