1. Genetics and Genomics
Download icon

Population-scale proteome variation in human induced pluripotent stem cells

  1. Bogdan Andrei Mirauta
  2. Daniel D Seaton
  3. Dalila Bensaddek
  4. Alejandro Brenes
  5. Marc Jan Bonder
  6. Helena Kilpinen
  7. HipSci Consortium
  8. Oliver Stegle  Is a corresponding author
  9. Angus I Lamond  Is a corresponding author
  1. European Molecular Biology Laboratory, European Bioinformatics Institute, Wellcome Genome Campus, United Kingdom
  2. Centre for Gene Regulation & Expression, School of Life Sciences, University of Dundee, United Kingdom
  3. Wellcome Sanger Institute, United Kingdom
  4. King's College, United Kingdom
  5. European Molecular Biology Laboratory, Genome Biology Unit, Germany
  6. Division of Computational Genomics and Systems Genetic, German Cancer Research Center, Germany
Research Article
  • Cited 0
  • Views 778
  • Annotations
Cite this article as: eLife 2020;9:e57390 doi: 10.7554/eLife.57390

Abstract

Human disease phenotypes are driven primarily by alterations in protein expression and/or function. To date, relatively little is known about the variability of the human proteome in populations and how this relates to variability in mRNA expression and to disease loci. Here, we present the first comprehensive proteomic analysis of human induced pluripotent stem cells (iPSC), a key cell type for disease modelling, analysing 202 iPSC lines derived from 151 donors, with integrated transcriptome and genomic sequence data from the same lines. We characterised the major genetic and non-genetic determinants of proteome variation across iPSC lines and assessed key regulatory mechanisms affecting variation in protein abundance. We identified 654 protein quantitative trait loci (pQTLs) in iPSCs, including disease-linked variants in protein-coding sequences and variants with trans regulatory effects. These include pQTL linked to GWAS variants that cannot be detected at the mRNA level, highlighting the utility of dissecting pQTL at peptide level resolution.

Introduction

Induced pluripotent stem cells (iPSC) hold great promise for advancing basic research and biomedicine. By enabling the in vitro reconstitution of development and cell differentiation, iPS cells allow the investigation of mechanisms underlying development and the aetiology of many forms of genetic disease. To realize this potential, it is essential to characterise genetic and non-genetic sources of variability of molecular and cellular phenotypes in human iPSCs.

Recently, multiple reference panels of human iPSC lines have been established (Kilpinen et al., 2017; Panopoulos et al., 2017; Carcamo-Orive et al., 2017), providing a valuable resource for functional experiments in pluripotent cells. These cell lines, together with associated data, have enabled the characterisation of variability in iPSC transcriptomes, identifying genetic and non-genetic determinants of expression variation, including expression quantitative trait loci (eQTL) (Kilpinen et al., 2017; Rouhani et al., 2014; DeBoever et al., 2017) in cis.

While RNA-centric analyses are informative for studying gene regulatory mechanisms at the transcriptional level, most cellular phenotypes ultimately involve downstream mechanisms that are mediated by proteins. Several proteogenomics studies, primarily in cancer (Zhang et al., 2014; Mertins et al., 2016), have underlined the relevance of protein measurements to interpreting how genomic changes act at the phenotypic level. Moreover, recent evidence has shown that genetic alterations can have effects on RNA that are attenuated at the protein level (Gonçalves et al., 2017; Roumeliotis et al., 2017). Vice versa, the mapping of protein quantitative trait loci (pQTL), predominantly in lymphoblast cell lines (Battle et al., 2015; Stark et al., 2014; Wu et al., 2013) and for the plasma proteome (Sun et al., 2018; Yao et al., 2018; Liu et al., 2015; Johansson et al., 2013; Lourdusamy et al., 2012), has revealed genetic effects on protein traits that do not manifest at the RNA level. However, the extent of RNA-independent protein regulation is not yet understood, with previous analyses performed only at gene-level resolution and, in some cases, without comparing protein and RNA data from the same cellular material.

Here, we report the first comprehensive, population-scale, combined proteomic and transcriptomic analysis of human iPSC lines. Our data provide matched quantitative proteomic (Tandem Mass Tag Mass Spectrometry) and transcriptomic (RNA-Seq) profiles of 202 iPSC lines, derived from 151 donors from the HipSci project (Kilpinen et al., 2017). We identify both genetic and non-genetic effects associated with variability in protein expression between individuals and describe the first high-resolution pQTL map in human iPSCs, including loci not detected as eQTLs at the RNA level.

Results

A population reference proteome for human iPSCs

A set of 217 iPSC lines from the HipSci project (Kilpinen et al., 2017), derived from 163 distinct donors, was selected for protein analysis, using material from the identical batches of cells that were used for RNA-Seq and other assays (Materials and methods). Quantitative mass spectrometry was carried out in batches of 10 lines, using tandem mass tagging (TMT, Thompson et al., 2003), with one common reference sample shared across batches (Brenes et al., 2018) (Materials and methods). Collectively, we identified 255,015 distinct (unmodified) peptide sequences, corresponding to 16,773 protein groups (groups of protein isoforms with no discriminating peptides; hereon denoted proteins) with median sequence coverage of 46%.

After quality control, 202 lines (from 151 donors) with matched genotype, RNA-Seq and proteome data, were selected for further analysis (Figure 1a; Figure 1—figure supplement 1; Figure 1—source data 1). We identified 11,140 recurrently detected proteins, corresponding to 10,198 genes (detected in at least 30 lines; Materials and methods) and RNA expression for 12,363 protein-coding genes (population average TPM >1). Out of these, 9013 protein coding genes were detected both at the RNA and protein levels (Figure 1—source datas 2 and 3).

Figure 1 with 6 supplements see all
| Characterising variation in the iPSC proteome and transcriptome.

(a) Experimental design and assays considered in this study. Genotype, RNA-Seq and quantitative proteomics data were obtained from the same cell material of 202 iPSC lines derived from 151 unrelated donors. (b) Variance component analysis of RNA and protein abundances across genes, considering different technical and biological factors. Shown is the distribution of the fraction of variance explained by different factors (upper panel) across proteins, and the number of genes with substantial variance contribution for each factor (>20% contribution; lower panel). Also shown are the number of genes that retain greater than 20% contribution after adjusting for the effect of the corresponding RNA profiles on protein abundance (light blue; see Materials and methods). (c) Association of protein level with X chromosome inactivation (XCI) status across 110 female iPSC lines. Shown are lowess regression curves for 322 and 312 proteins respectively that were identified as significantly up (red) - and down (blue) - regulated with loss of XCI in female iPSC lines (lower panel; 10% FDR). Selected gene ontology enrichments for these sets of proteins are shown (right-hand panel; Materials and methods). XCI status was estimated as the average fraction of allele-specific expression for the inactive chromosome across all chromosome X genes (Materials and methods). (d) Scatter plot of the fraction of variance explained by donor at the RNA (x-axis) versus protein (y-axis) level. Encoded in colour is the fraction of variance explained by donor effects at the protein level after adjusting for the effect of the corresponding RNA profiles on protein abundance (Materials and methods).

Collectively, these data provide the most comprehensive analysis of the human iPSC proteome reported to date, and one of the most comprehensive proteomic datasets reported for any human primary, or derived, cell type (Supplementary file 1). When overlaying our data with the Human Proteome Map (Kim et al., 2014), the iPSC proteome was most similar to foetal and reproductive organs (Figure 1—figure supplement 2), which is consistent with the expected expression of pluripotency markers in these tissues (Kerr et al., 2008a; Kerr et al., 2008b). We also assessed differences between healthy donors and disease-bearing donors, identifying no systematic expression differences in the iPSC state (Supplementary file 2; Figure 1—figure supplement 3).

RNA and proteome variability

We assessed a range of factors to explain the variation in protein expression between iPSC lines. Leveraging our experimental design with data from two or more lines for 34% of donors, we assessed the effects of donor, alongside age, sex, and the contributions of technical and cell culture related factors (on 6009 genes; Figure 1b,d; Figure 1—source data 4; using a linear mixed model; Materials and methods). Overall, the fraction of variance explained by biological factors was lower for protein levels, compared to RNA variation, which points to higher assay noise and/or stochastic variability of protein abundance. Consistent with previous results using RNA data (Kilpinen et al., 2017), we identified donor genome (i.e. DNA sequence variation) as the most relevant factor, followed by culture medium (Figure 1b). Critically, however, significant donor effects remained after accounting for RNA variability (Figure 1b,d; Figure 1—figure supplement 4; Materials and methods). This indicates that (i) genetic differences between individual donors and experimental differences between culture conditions play an important role in causing the observed variation in proteome expression between the iPSC lines and (ii) post-transcriptional mechanisms also contribute to these effects. Notably, many of the proteins showing strong donor effects were previously identified as differentially expressed between reprogrammed iPS cells and embryonic stem cells (ESCs) (Phanstiel et al., 2011; Munoz et al., 2011), suggesting that some of these previously reported differences could be due to genetic variation between donors, rather than intrinsic differences between the iPSC and ESC cell types (Figure 1—figure supplement 5).

The sex of the donor affected proteome expression, including a subset of proteins uniquely encoded on the male-specific X chromosome. There was also a strong (i.e. >20%) gender-related effect on the expression of a subset of 88 proteins (Figure 1b), which are enriched for proteins encoded on the X chromosome (Odds ratio = 24.8, PV = 3×10−32, Fisher’s exact test). This reflected the partial erosion of X chromosome inactivation (XCI) observed in a subset of iPSC lines derived from female donors, as confirmed both by quantification of allele-specific expression and XIST expression in these lines (Figure 1—figure supplement 6). Incomplete XCI has been linked previously to poor iPSC differentiation and changes in RNA levels (Mekhoubad et al., 2012; Salomonis et al., 2016). However, our data provide the first opportunity to link the XCI status of 110 distinct female iPSC lines (as inferred from allele-specific expression; Materials and methods; Figure 1—figure supplement 6), with changes in the abundance of both proteins and RNAs. This identified 1374 genes for which either protein, or RNA levels, or both, showed changes correlated with XCI status (Figure 1c, FDR < 10%; Figure 1—source data 5). Further analysis indicated that XCI status preferentially impacts catabolic processes and mitochondria at the protein level, while this was not observed at the RNA level (Gene Ontology; Materials and methods; Figure 1c; Figure 1—source data 6). These data thus reveal an important effect of XCI status on global gene expression in iPSC lines from female donors, including specific effects at the protein level that are not detected by transcriptomic analysis.

Mapping cis genetic effects on protein abundance

Next, we mapped cis quantitative trait loci at both the RNA and protein levels, considering 8543 autosomal protein coding genes that were quantified at both levels (MAF >5%; within +/- 250 kb of the gene boundaries; using a linear mixed model; Materials and methods). The number of pQTLs identified was greatly increased by adapting the PEER-based adjustment, which was previously developed for mapping of eQTLs (Stegle et al., 2012), for use with proteomic data (Materials and methods; Figure 2—figure supplement 1). Across all autosomal genes, we report 654 genes with at least one cis pQTL and 3487 genes with a cis eQTL (FDR < 10% for both eQTL and pQTL mapping; lead variants only; Figure 2a; Figure 2—source datas 1 and 2). Among these, 273 genes were shared and had identical or correlated lead variants, whereas 82 genes showed evidence for an eQTL and pQTL with independent lead variants (LD-based criterion, r2 <0.1, Figure 2—source data 1). Genes with substantial donor components, as identified based on the variance component analysis (>20%; Figure 1b), were enriched for significant cis pQTL (215 genes out of 962; Odds ratio = 4.2; Figure 2—figure supplement 2).

Figure 2 with 4 supplements see all
Human iPSC cis protein and RNA QTLs.

(a) Number of genes with a protein (blue) or RNA (green) cis QTL (FDR < 10%) and pairwise replication of genetic effects. Left: Number of genes with a pQTL, either with (dark blue) or without (light blue) replicated RNA effect. Right: Number of genes with an eQTL, either with (dark green) or without (light green) replicated protein effect. Replication defined by assessing nominal significance (PV <0.01) of QTL in the respective other layer. (b) Local Manhattan plots displaying negative log p-values (PV) from cis RNA (top) and protein (bottom) QTL mapping for PEX6. The dashed line and the grey box indicate the genomic positions of the lead QTL and of the gene. Boxplots show RNA and protein expression for different alleles at the pQTL lead variant rs11752813, a variant in LD (r2 = 1, 1000 Genomes European populations phase 3) with the Alzheimer risk variant rs1129187 (Jun et al., 2016) (OR 1.13). (c) Cumulative fraction of eQTLs with replicated protein effects as a function of the eQTL effect size (from highest to lowest). (d) Prediction of protein replication of eQTLs, considering features derived from gene annotations, eQTL, RNA and protein data. Predictions were obtained using a random forest model trained on the protein replication status of eQTL (as in a; Materials and methods). Left: Feature importance scores. Right: Precision-recall curve for the model, evaluated in independent test fractions. The model performance was assessed by random sampling of training/testing data with a 80/20 split, performed 50 times. Shown in red is the average precision-recall across all sampled training/test splits and in thin grey lines results of individual folds.

To identify DNA sequence variants with effects at both the RNA and protein levels, we considered the pairwise replication of pQTLs at the RNA level and vice versa (lead QTL variants, defining ‘replication’ as nominal PV <0.01 with consistent effect direction; Materials and methods), which is more sensitive than assessing overlapping QTL variants. This identified 473 pQTLs (72%) with replicated eQTL effects. Conversely, 893 eQTLs (26%) had replicated protein effects, with globally concordant effect size directions and distance from gene boundaries (Figure 2—figure supplement 3).

Lack of replication of eQTLs at the protein level could arise from a combination of technical and/or biological factors. We identified the eQTL effect size as a strong predictor for protein replication, with larger effects being associated with increased replication rates (Figure 2c). To systematically characterise the determinants of eQTL replication, we considered a random forest model trained to predict the protein replication status (Figure 2d). In addition to the eQTL effect size, this identified other predictive factors, including the protein coefficient of error (estimated from technical replicate samples; Materials and methods) and the protein coefficient of variation across lines (Figure 2d; Figure 2—source data 2).

To explore the physiological relevance of iPSC pQTL variants, we examined their overlap with variants identified in genome-wide association studies (GWAS). Specifically, we probed for QTLs that tag GWAS variants contained in the GWAS catalogue (MacArthur et al., 2017) (i.e. are in LD r2 >0.8), identifying 136 (of 654) pQTL signals that tag a known GWAS variant (Figure 2—source data 1). In addition, we assessed the statistical evidence for co-localisation of pQTL and GWAS signals for 51 studies for which full summary statistics were obtained (using eCAVIAR Hormozdiari et al., 2016; Materials and methods; Figure 2—source data 1), yielding 49 pQTLs with evidence of co-localisation (i.e cumulative co-localisation probability greater than 0.1). Among these, examples of pQTLs with corresponding effects at the RNA level include the variant rs7872034, a pQTL for SMC2 with co-localisation evidence for serous invasive ovarian cancer (Phelan et al., 2017), and the variant rs11752813, a pQTL for PEX6 and in LD with Alzheimer's disease in APOE e4+ carriers risk variant rs1129187 (Jun et al., 2016; Figure 2—figure supplement 4; Figure 2b).

Notably, for 33 pQTLs linked to GWAS variants, either via co-localisation or LD tagging, no replicated effect was identified at the RNA level, suggesting protein-specific regulation (Figure 2—source data 1). For example, rs11601507 has no RNA effects, and is associated with TRIM5 protein abundance and with coronary artery disease risk (van der Harst and Verweij, 2018; Figure 2—figure supplement 4). Such cases raise the question of the mechanisms by which these variants modulate protein abundance and, ultimately, phenotypic traits, as addressed below.

pQTL linked to isoform-specific transcript expression

To investigate the mechanisms that underlie discordant eQTLs and pQTLs in more detail, we performed transcript isoform and protein peptide QTL analyses. cis QTL mapping of 33,050 reference transcript isoforms (Zerbino et al., 2018) (quantified using Salmon Patro et al., 2017; Materials and methods) and 119,747 peptides identified 3810 genes with a transcript QTL (tQTL) and 566 genes with a peptide QTL (pepQTL), respectively (Figure 3—figure supplement 1, Materials and methods, Figure 3—source datas 1 and 2).

Transcript-level QTL mapping could explain the lack of protein effects for a small fraction of the 2594 eQTL without a replicated pQTL effect (Figure 2a). For 48 of these, the eQTL variant was identified as exclusively associated with abundance changes of non-coding transcript isoforms (nominal PV <0.01), which explains the absence of protein effects (Figure 3—figure supplement 2). Furthermore, when considering 1262 transcript QTL that neither replicate at the eQTL, nor at the pQTL level, in 45 instances we observed consistent replication when considering peptide QTL (Figure 3—figure supplement 2b).

Among 181 pQTLs without eQTL replication (Figure 2a), 61 had nominally significant transcript QTLs (PV <0.01; Figure 3a). For 12 of these, including a pQTL for MMAB (Figure 3b), we observed genetic effects with opposite directions on coding and non-coding transcript isoforms, which explains the lack of genetic effects when considering gene-level RNA abundance.

Figure 3 with 4 supplements see all
Putative mechanisms of pQTL from transcript isoform regulation and protein-altering variants.

(a) Categorisation of 654 pQTL into four classes according to their putative mechanism: gene expression effect (i.e. replicated at eQTL level), transcript-isoform specific effect (i.e. not replicated at eQTL level, but significant at transcript isoform level), protein-altering variant (i.e. at least one inframe variant in LD with lead pQTL variant) without expression effect at RNA level, and without any putative mechanism identified. (b) Example pQTL without eQTL replication (rs6663; gene MMAB), with a directional opposite effect on a coding and non-coding isoform (cyan: ENST00000540016; grey: ENST00000537496), resulting in no overall change in gene expression level. (c) The pQTL variant (rs1051061) is a protein-altering variant associated with VRK2 protein abundance (below), and lacks detectable effect on RNA expression. The pQTL signal is observed across 15 peptides spanning the VRK2 protein sequence (above, left). This variant is associated with schizophrenia risk, and is located at the kinase active site, proximal to the proton acceptor residue (above, right). The dashed line and the grey box indicate the genomic positions of the lead QTL and of the gene. (d) Enrichment of RNA-independent pQTL in different categories of predicted variant effects, using gene variants in high LD with pQTLs (proxy gene variants; r2 >0.8; within the cis gene boundaries). Enrichment calculated using Fisher’s exact test.

pQTL arising from protein-altering variants

Next, we set out to characterise further the remaining 120 pQTL without replication at either eQTL, or transcript QTL levels. When classifying the corresponding lead pQTL variants based on their predicted functional effect, we identified 24 inframe variants, a striking enrichment compared to pQTL with replicated RNA effects (3.8-fold enrichment; PV = 4×10−5, Fisher’s exact test; Figure 3c and d). These findings are in line with previous observations in lymphoblast cell lines (Battle et al., 2015). Of note, peptides containing protein-altering variants were excluded from the quantifications (Materials and methods), and the reported pQTL effects were observed for multiple peptides from the same proteins (Figure 3—figure supplement 3), providing further confidence in genuine regulatory effects. We assessed whether the 24 pQTL have effects at the RNA level (eQTL) in other cell types, and for 11 of these pQTL we did not find evidence of eQTL nominal significance in any of the 48 GTEx (PV <0.01/48; Battle et al., 2017; Figure 1—source data 1), which further points to RNA-independent mechanisms.

Inframe variants have the potential to affect protein function. We estimated whether a variant is likely to be deleterious to protein function using SIFT scores, which capture evolutionary conservation and amino acid similarity (Ng and Henikoff, 2003). This revealed a clear enrichment of the 24 RNA-independent pQTL that tag inframe variants, 10 of which have predicted deleterious effects (SIFT score <0.05), compared to four among all other pQTL (Odds ratio = 27.5, PV = 3.8×10−8, Fisher’s exact test; Figure 3d; Figure 1—source data 1). Putative effects of these variants on protein function include loss of enzymatic activity and disruption of protein structure. For example, the variant rs1051061 in VRK2 lies in a conserved sequence in the kinase domain, proximal to the proton acceptor residue, likely impacting kinase activity (Figure 3c). The identical variant has been identified as GWAS risk variant for schizophrenia (Yu et al., 2017) (OR 1.17), with the risk allele being associated with decreased protein abundance. The effect direction is consistent with previous studies that have linked decreased VRK2 expression to neurological disorders including schizophrenia (Azimi et al., 2018; Tesli et al., 2016).

These data show important roles of transcriptional regulation underlying cis pQTL effects, while also highlighting how isoform-specific effects, which are invisible to standard eQTL mapping approaches, can be detected at the protein level. For a substantial subset of pQTLs, we identified linked protein-altering variants, many with deleterious effects. Together with previous observations, these results suggest that proteomics information can aid understanding of pathogenic mechanisms of deleterious variants.

Proteome-wide effects of cis QTLs

Building on the compendium of cis pQTL identified here in iPS cells, we set out to characterise downstream proteome-wide changes. We mapped proteome-wide trans pQTL, considering 654 cis pQTL variants. This identified 51 cis-pQTL lead variants with trans effects on a total of 68 proteins (FDR < 10%; Figure 4—source data 1; Materials and methods). To rule out synthetic associations, we discarded associations with evidence for sequence similarity between cis and trans proteins, and we verified the consistency of the identified trans effects across multiple independent peptides (Materials and methods; Figure 3—figure supplement 3). The detected pairs of proteins with shared genetic regulation were strongly enriched for known protein-protein interactions (CORUM Ruepp et al., 2010, IntAct Orchard et al., 2014, StringDB Szklarczyk et al., 2017; Odds ratio = 9.1, PV = 1.5×10−10, Fisher’s exact test; Figure 4b). The cis and trans effects had similar effect directions and effect sizes, consistent with genetic effects mediated via stabilising protein-protein interactions (Figure 4c). This interpretation of our data in human iPSCs is consistent with the significant donor variance component we observed for many protein complexes (Figure 4d). It is also consistent with previous observations in an outbred mouse cross, showing that protein modules sharing genetic effects in trans are enriched in protein interactions (Ruepp et al., 2010), and identification of trans protein effects due to somatic aberrations in human cancer cell lines (Gonçalves et al., 2017; Roumeliotis et al., 2017). Importantly, our results generalise these previous observations to genetic effects of common variants that segregate in human populations.

| Trans effects on the iPSC proteome.

(a) Strategy for mapping trans genetic effects on protein abundance. Lead cis pQTL variants were considered for proteome-wide association analysis. (b) Enrichment of previously catalogued protein-protein interactions among significant trans pQTLs. Shown is the fraction of cis-trans gene pairs linked by a trans pQTL with evidence of protein-protein interactions (based on the union of CORUM, IntAct, and StringDB), as a function of the considered FDR threshold for trans pQTL discovery. The dashed lines correspond to FDR < 10%. Numbers indicate the number of trans pQTL identified for each FDR threshold. (c) Comparison of genetic effect sizes, in cis and trans, for significant (FDR < 10%) trans pQTLs. Red points indicate cis-trans pairs with evidence for protein-protein interactions defined as in b. (d) Left: Protein co-expression of protein complex subunits defined based on CORUM. Right: i) subunit with the most significant cis pQTL; ii) fraction of subunits in association with the cis pQTL at nominal significance (PV <0.01). iii) fraction of the average cluster protein expression level explained by donor effects. (eTrans regulation of the PEX26-PEX6-PEX1 complex. The variant rs11752813 (LD r2 = 1 with rs1129187) is associated in cis with changes in the RNA and protein abundance of PEX6 and in trans with changes in the protein abundance of PEX1 and PEX26.

In summary, the trans effects we detected appear to induce strong correlations across protein complex subunits (Figure 4d), whereby a variant associated in cis with one subunit was also associated in trans with other subunits. This is illustrated by PEX26-PEX6-PEX1, a protein complex involved in peroxisome biogenesis. As noted above, the underlying pQTL variant rs1129187 is associated in cis with an increase in both PEX6 RNA and protein abundance (Figure 2b) and is a known risk variant for Alzheimer's disease in APOE e4+ carriers (Jun et al., 2016). This cis pQTL in turn induces downstream associations on the remaining complex subunits, PEX26 and PEX1 (Figure 4e), suggesting that PEX6 acts as a limiting subunit of this complex in iPSCs. Thus, our results provide a potential biological mechanism underlying this risk variant, namely acting through changes in the abundance of the PEX26-PEX6-PEX1 complex. Notably, there is prior evidence for an implication of peroxisomal function in the development of Alzheimer’s disease and in other neurodegenerative processes (Lizard et al., 2012; Berger et al., 2016), providing further support for this hypothesis.

Discussion

Here, we report the first in-depth characterisation of the human iPSC proteome, connecting genetic variation to changes in RNA and protein levels. Beyond the relevance for iPS cell biology, this study, to our knowledge, provides the most detailed population-level analysis of parallel RNA/protein profiles in human cells. By quantifying genome-wide protein and transcript expression variation across more than 200 human iPSC lines, we identified both genetic and non-genetic mechanisms that underlie variation in both protein and RNA levels. We have mapped more than 600 cis protein quantitative trait loci (pQTLs) and analysed how these relate to cis eQTLs, how they impact other proteins in trans, and how pQTLs link to human disease variants.

The variance component analysis explained a lower overall fraction of variance in the protein data compared to RNA variation, which likely reflects larger technical effects and/or stochasticity in protein expression levels. Among the explainable fraction of variance, donor-specific genetic factors are a major contributor to the differences in protein expression detected across the iPSC lines. The corollary is that protein expression variation across iPS cells reflects genetic diversity in the human population. Consistent with this, we identified 654 common genetic variants associated with changes in protein abundance.

Globally, there were substantially fewer pQTLs than eQTLs, and while most pQTLs had effects of the same direction at eQTL, only 30% of eQTLs are nominally significant at the protein level. It is possible that technical factors resulting from the protein measurement methods may contribute, at least in part, towards attenuating the signal detected at the protein level. However, considering our data in light also of results from previous studies, some of which employed alternative technologies for protein detection to the MS methods used here, we suggest that the signal attenuation between eQTL and pQTL levels is not exclusively the result of limitations in protein measurements. Instead, many eQTLs may reflect variation in RNA abundance that does not cause significant changes in steady state protein levels.

By the systematic comparison of matched protein and RNA data, including detailed analysis of separate isoforms, we demonstrated that in order to fully understand the propagation of genetic effects to proteins, isoform-resolution protein and RNA phenotypes are indispensable. In particular, this approach identified additional RNA-dependent regulation that manifests in protein QTL, thereby improving the ability to identify genuine RNA-independent pQTL.

We showed that the pQTLs for which no corresponding changes in transcript levels were detected, are enriched in deleterious missense variants. This result suggests that the phenotypic effects of such variations may be exerted through protein abundance changes. Because most deleterious variants, and in particular pathogenic variants, are rare, larger sample sizes will be required to fully assess the protein components of this class of regulatory genetic effects.

Our study presents the first comprehensive map of pQTLs at peptide resolution, considering a total of 119,747 peptides from 8543 proteins for genetic analysis. This identified 566 peptide QTL, several of which were not detectable when considering whole protein expression levels, as illustrated with the variant rs12795503, pepQTL for gene CTTN (Figure 3—figure supplement 1). While we mapped fewer significant pepQTLs than pQTLs, peptide level analyses were shown here to overcome potential artefacts raised by protein quantification, in particular when mapping trans pQTL, and are invaluable in identifying isoform-specific effects.

Our data highlight the ability of protein-protein interactions to propagate genetic effects in human populations. A long-standing hypothesis has been that certain protein complexes may have a rate-limiting subunit that determines complex abundance, with any excess subunits produced being rapidly degraded (e.g. because of exposure of hydrophobic residues). This implies that cis eQTLs affecting the levels of rate-limiting subunits should also have effects in trans on the abundance of the whole complex, and on most, if not all, subunits therein. While trans genetic effects were previously reported to be mediated by protein interactions in high heterozygosity samples, that is outbred mice, (Chick et al., 2016) and for somatic aberrations in cancer cell lines (Gonçalves et al., 2017; Roumeliotis et al., 2017), to our knowledge, this study provides the first example that such effects act through common genetic variants in untransformed human cells. In the future, the approach we have taken here could be extended by mendelian randomisation-based approaches to formally assess the mediating role of the cis pQTL on protein complex members.

Understanding the mechanisms through which genetic variations act in the human population is of great relevance to characterising risk factors and susceptibility to disease. There is on-going interest in the potential for studying disease mechanisms using disease relevant tissues that are derived from panels of iPSCs (Cayo et al., 2017; Li et al., 2018; D'Aiuto et al., 2014; Schwartzentruber et al., 2018). This study provides important information showing how direct analysis of human iPSCs can advance our understanding of the genetic regulation of protein expression and how this influences cell phenotypes and disease mechanisms.

Materials and methods

Key resources table
Reagent type
(species) or
resource
DesignationSource or
reference
IdentifiersAdditional
information
Cell line (Homo-sapiens)iPSCwww.hipsci.orgRRID:SCR_003909
Software, algorithmMaxQuanthttps://www.maxquant.org/RRID:SCR_014485
Software, algorithmTrim Galorehttps://www.bioinformatics.babraham.ac.uk/projects/trim_galore/RRID:SCR_011847
Software, algorithmSTARhttps://www.ncbi.nlm.nih.gov/pmc/articles/PMC3530905/RRID:SCR_015899
Software, algorithmSalmonhttps://combine-lab.github.io/salmon/NA
Software, algorithmg:Profilerhttp://biit.cs.ut.ee/gprofiler/RRID:SCR_006809
Software, algorithmeCAVIARhttp://zarlab.cs.ucla.edu/tag/ecaviar/NA
Software, algorithmVEPhttps://www.ensembl.org/info/docs/tools/vep/index.htmlRRID:SCR_007931
Software, algorithmMutFunchttp://www.mutfunc.com/NA
Software, algorithmLimixhttps://github.com/limix/limixNA
Software, algorithmPeerhttp://www.sanger.ac.uk/resources/software/peer/RRID:SCR_009326
Software, algorithmScikit-learnhttp://scikit-learn.org/RRID:SCR_002577

Cell lines

Request a detailed protocol

As described in Kilpinen et al., 2017, all HipSci samples were collected from consented research volunteers recruited from the NIHR Cambridge BioResource, and iPSC were generated from fibroblasts by transduction with Sendai vectors. In brief, cells were cultured on fibroblasts (MEF-CF1) feeder layer with selected lines being transferred on Essential 8 (E8) medium. Pluripotency was assessed based on expression profiling, detection of pluripotency markers in culture and response to differentiation inducing conditions. Mycoplasma screening was performed using a standard PCR kit. Sample identity was confirmed using a Fluidigm genotyping assay containing 24 SNPs. The ID numbers and details for each cell line included in this analysis are listed in Figure 1—source data 1.

RNA-Seq data processing

Request a detailed protocol

Raw RNA-Seq data for 331 samples were obtained from the ENA project: ERP007111. CRAM files were merged on a sample level and were converted to FASTQ format. Sequencing reads were trimmed to remove adapters and low-quality bases (using Trim Galore! Dobin et al., 2013), followed by read alignment using STAR (version: 020201) (Liao et al., 2014), using the two-pass alignment mode and the default parameters as proposed by ENCODE (c.f. STAR manual). All alignments were relative to the GRCh37 reference genome, using ENSEMBL 75 as transcript annotation (Zerbino et al., 2018).

Samples with low-quality RNA-Seq were discarded, if they had less than 2 billion bases aligned, had less than 30% coding bases, or had a duplication rate higher than 75%. This resulted in 323 lines for analysis, of which 202 had matched proteome data.

Gene-level RNA expression was quantified from the STAR alignments using featureCounts (v1.6.0) (Robinson et al., 2010), which was applied to the primary alignments using the ‘-B’ and ‘-C’ options in stranded mode, using the ENSEMBL 75 GTF file. Quantifications per sample were merged into an expression table using the following normalisation steps. First, gene counts were normalized by gene length. Second, the counts for each sample were normalised by sequencing depth using the edgeR adjustment (McAlister et al., 2014).

Transcript isoform expression was quantified directly from the (unaligned) trimmed reads using Salmon (Zerbino et al., 2018) (version: 0.8.2), using the ‘--seqBias’, ‘--gcBias’ and ‘VBOpt’ options in ‘ISR’ mode to match our inward stranded sequencing reads. The transcript database was built on transcripts derived from ENSEMBL 75. The TPM values as returned by Salmon were combined into an expression table.

Quantitative proteomics data generation

Request a detailed protocol

Upon establishment in culture, lines were expanded for banking and molecular assays (including genomics, transcriptomics and proteomics). We selected 217 lines (banked frozen pellets) for in-depth proteomic analysis with Tandem Mass Tag Mass Spectrometry. A subset of 202 lines (112 normal and 90 disease; Figure 1—source data 1) with matched mRNA and protein data were considered for further analysis. Previous studies have shown that this sample size enables the identification of a large number of cis RNA and protein QTLs (Battle et al., 2015).

Sample preparation

Request a detailed protocol

For protein extraction, frozen iPSC cell pellets were washed with ice cold PBS and redissolved immediately in 200 μL of lysis buffer (8 M urea in 100 mM triethyl ammonium bicarbonate (TEAB) and mixed at room temperature for 15 min. DNA in the cell lysates was sheared using ultrasonication (6 × 20 s at 10°C). The proteins were reduced using tris-carboxyethylphosphine TCEP (25 mM) for 30 min at room temperature, then alkylated in the dark for 30 min using iodoacetamide (50 mM). Total protein was quantified using the fluorescence based EZQ assay (Life Technologies). The lysates were diluted four-fold with 100 mM TEAB for the first protease digestion with mass spectrometry grade lysyl endopeptidase, Lys-C (Wako, Japan), then diluted a further 2.5-fold before a second digestion with trypsin. Lys-C and trypsin were used at an enzyme to substrate ratio of 1:50 (w/w). Digestions were carried out for 12 hr at 37°C, then stopped by acidification with trifluoroacetic acid (TFA) to a final concentration of 1% (v:v). Peptides were desalted using C18 Sep-Pak cartridges (Waters) following manufacturer’s instructions and dried.

Tandem Mass Tag mass spectrometry analysis

Request a detailed protocol

For Tandem Mass Tag (TMT)-based quantification, the dried peptides were redissolved in 100 mM TEAB (50 µL) and their concentration was measured using a fluorescent assay (CBQCA) (Life Technologies). 100 µg of peptides, from each cell line to be compared, in 100 µL of TEAB were labelled with a different TMT tag (20 µg ml−1 in 40 µL acetonitrile) (Thermo Scientific), for 2 hr at room temperature. After incubation, the labelling reaction was quenched using 8 µl of 5% hydroxylamine (Pierce) for 30 min and the different cell lines/tags were mixed and dried in vacuo. TMT-ten plex was used to label ten IPSC lines and quantify them in parallel. In total, 24 TMT-ten plex experiments were performed, where one IPSC line (HPSI0314i-bubh_3) was chosen as a reference cell line and was kept constant in all TMT batches. The other nine quantification channels were used to label 9 different cell lines. No randomisation was applied in assigning the samples to batches.

The TMT samples were fractionated using off-line, high pH reverse phase chromatography: samples were loaded onto a 4.6 × 250 mm Xbridge BEH130 C18 column with 3.5 µm particles (Waters). Using a Dionex bioRS system, the samples were separated using a 25 min multistep gradient of solvents A (10 mM formate at pH 9) and B (10 mM ammonium formate pH 9 in 80% acetonitrile), at a flow rate of 1 ml/min. Peptides were separated into 48 fractions, which were consolidated into 24 fractions. The fractions were subsequently dried and the peptides redissolved in 5% formic acid and analysed by LC-MS.

5% of the material was analysed using an orbitrap fusion tribrid mass spectrometer (Thermo Scientific), equipped with a Dionex ultra high-pressure liquid chromatography system (nano RSLC). RP-LC was performed using a Dionex RSLC nano HPLC (Thermo Scientific). Peptides were injected onto a 75 μm × 2 cm PepMap-C18 pre-column and resolved on a 75 μm × 50 cm RP- C18 EASY-Spray temperature controlled integrated column-emitter (Thermo), using a 4-hr multistep gradient from 5% B to 35% B with a constant flow of 200 nL min−1. The mobile phases were: 2% ACN incorporating 0.1% FA (Solvent A) and 80% ACN incorporating 0.1% FA (Solvent B). The spray was initiated by applying 2.5 kV to the EASY-Spray emitter and the data were acquired under the control of Xcalibur software in a data-dependent mode using top speed and 4 s duration per cycle. The survey scan is acquired in the orbitrap covering the m/z range from 400 to 1400 Th, with a mass resolution of 120,000 and an automatic gain control (AGC) target of 2.0 e5 ions. The most intense ions were selected for fragmentation using CID in the ion trap with 30% CID collision energy and an isolation window of 1.6 Th. The AGC target was set to 1.0 e4 with a maximum injection time of 70 ms and a dynamic exclusion of 80 s.

During the MS3 analysis for more accurate TMT quantifications, five fragment ions were co-isolated using synchronous precursor selection, using a window of 2 Th and further fragmented using a HCD collision energy of 55% (McAlister et al., 2014). The fragments were then analysed in the orbitrap with a resolution of 60,000. The AGC target was set to 1.0 e5 and the maximum injection time was set to 105 ms.

Proteomics data processing

Request a detailed protocol

The TMT-labelled samples (24 batches of TMT-ten plex) were analysed using MaxQuant v. 1.6.0.13 (Cox et al., 2011; Schwanhäusser et al., 2011). Proteins and peptides were identified using the UniProt human reference proteome database (Swiss Prot + TrEMBL) release-2017_03, using the Andromeda search engine. Run parameters and the raw MaxQuant output have been deposited at PRIDE (PXD010557).

The following search parameters were used: reporter ion quantification, mass deviation of 6 ppm on the precursor and 0.5 Da on the fragment ions; Tryp/P for enzyme specificity; up to two missed cleavages, ‘match between runs’, ‘iBAQ’. Carbamidomethylation on cysteine was set as a fixed modification. Oxidation on methionine, pyro-glu conversion of N-terminal Gln, deamidation of asparagine and glutamine and acetylation at the protein N-terminus, were set as variable modifications (Cox and Mann, 2008; Cox et al., 2011; Schwanhäusser et al., 2011; Tyanova et al., 2016).

Peptides and protein groups, that is groups of protein isoforms identified by common peptides (for details see Brenes et al., 2018), were reported at FDR < 5%. The same FDR threshold was used for reporting the Peptide Spectrum Matches (PSM). We performed the FDR calculation on an extended set and removed the Razor Protein FDR calculation constraint (for more details see reference Ramasamy et al., 2013). In total, we identified 255,015 peptides detected in at least one sample (after removing reverse and contaminant peptides; on the 217 lines and 23 replicates of the reference line), which corresponds to 16,773 protein groups.

Quality control and quantification

Request a detailed protocol

Peptides that overlap non-synonymous variants may be incorrectly detected and could result in synthetic associations, similarly to the polymorphism-in-probe effects in microarrays (Schlaffner et al., 2017). We performed the following quality control and filtering steps. First, using PoGo; (Lippert et al., 2015) (applied with Gencode release 30 mapped on GRCh 37), we reconstructed the genomic loci of 250,171 peptides that could be assigned to one or multiple genomic locations. Next, we assessed the overlap of these peptides with common, non-synonymous variants in human populations (MAF >0.01 in 1000 Genomes European populations phase 3). To define non-synonymous variants, we considered the overlap with any transcript isoform that contains the analysed peptide obtained from Gencode release 30 mapped on GRCh 37. This resulted in 6273 peptides that overlap such polymorphisms, which were discarded from further analyses (Supplementary file 3). Peptides mapped to multiple locations in the genome were discarded if they overlap non-synonymous variants for at least one of these locations. Finally, we discarded 4844 peptides for which a genetic position was not identified.

We discarded 10 lines with fewer than 67,000 identified peptides (corresponding to 75% of the median number of peptides identified; Figure 1—figure supplement 1), resulting in a proteomics dataset consisting of 207 lines, 202 of which had matched RNA-Seq data and hence were considered for further analysis. In addition, the technical replicates for the included reference line in each TMT batch were retained to aid the normalisation of protein quantifications between batches; see below.

Protein group abundances were estimated using the remaining peptides as the sum of the intensities of individual peptides mapped to the protein group. Peptide abundance estimates were obtained from the intensity values reported in the ‘Peptides’ file from MaxQuant.

For downstream analysis, we considered the subset of peptides that were recurrently detected in at least 30 of the 202 lines. Similarly, we discarded all protein groups that did not contain at least one recurrent peptide. This resulted in a final dataset of 11,140 recurrent protein groups and 132,078 recurrent peptides (Figure 1—figure supplement 1), corresponding to 10,198 genes (Ensembl ID).

To adjust for technical effects during the acquisition of protein data in TMT batches, we scaled the abundance estimate for each feature (i.e. protein or peptide) as follows. For each feature and batch, we multiply the intensity with a scaling coefficient defined as the ratio between the median intensity across all lines and the median across the subset of lines within a given TMT batch. Next, we employed quantile normalisation for peptide and protein abundance estimates, by performing quantile normalisation of the feature distribution in each line relative to a normalisation reference line (the line with the highest number of total peptides detected). Briefly, for each line we assigned for each feature the value observed in the reference line corresponding to the rank of the value of that feature in the line to be normalised: y’pl = r [rank ypl ], where ypl are the intensity values for feature p and line l obtained after batch scaling, r is the sorted vector of intensities from the normalisation reference line, and y’pl is the normalised value. In order to evaluate how the data transformation tackles batch effects, we performed a PCA analysis of protein quantifications and compared the peptide coefficient of error before and after data transformation (Figure 1—figure supplement 1d and e).

Assessment of TMT ratio compression effects

Request a detailed protocol

We assessed quantitative compression of our proteomics data by examining changes in peptide abundance for those peptides that were discarded because they overlap non-synonymous variants, following Battle et al., 2015. The rationale behind this approach is that a non-synonymous variant in a peptide prevents detection of that peptide, as its sequence will not exist in the proteome reference. Thus, in samples heterozygous for the non-synonymous variants, the measured peptide abundance is expected to be half of that of samples homozygous for the reference variant. Although the data indicated that ratio compression effects can be noticed in our study (Figure 3—figure supplement 4), the QTL results show that protein measurements derived from TMT Mass spectrometry analyses are suitable for the detection of protein abundance changes.

Comparisons of iPSC proteome profiles to existing datasets

Request a detailed protocol

To compare our iPSC proteome dataset to the Human Proteome Map (HPM) (Kim et al., 2014; Figure 1—figure supplement 1d), we first mapped the RefSeq IDs of proteins quantified in the HPM to UniProt IDs. We considered the subset of 8333 proteins that were expressed in our iPSC dataset and in at least one HPM tissue, and for which IDs could be mapped. We then calculated spearman correlation coefficients between the aggregate iPSC proteome abundance profile (averaged across lines) and each HPM tissue.

Variance component analysis

Request a detailed protocol

In order to calculate the contribution of each factor k to variation in protein abundance, we fitted a random effects model as follows: y = μ + ∑k uk + ε ; uk ∼ N(0, σk2 ⋅ Mk ); ε ∼ N(0, σr2I); Mk [i, j] = {1 if fk [i] = fk [j]; 0 if fk [i ] ≠ fk [j]).

where y denotes the (N x 1) vector of normalised log protein abundances, uk are the random effects, Mk is the (N x N) covariance structure, σkis the standard deviation, and ε is the residual (i.i.d. noise). The random effect components are defined based on a categorical covariance function defined on covariates fk, that is the vector of observed values for factor k (e.g. fk [i] ∈ {′male′, ′female′} when k is the donor sex component). We considered donor identity, donor sex, donor age (treated as a categorical variable, with each group in age windows of 5 years), culture medium, TMT batch, and TMT channel as random effect components. In order to accurately estimate the donor variance component, we restricted this analysis to the set of lines from the subset of 51 donors for which two cell lines were assayed and on 6009 genes with TPM >1 and with proteins detected in this subset of lines. Analogous analyses were considered for RNA abundance, leaving out the TMT-specific, random effects.

In order to study the effect of different factors on protein abundance, after adjusting for the effects of RNA abundance on protein abundance, we also applied the variance decomposition analysis to protein abundance values after adjusting for RNA variation. Adjusted protein abundances were calculated by regressing out the effects of RNA abundance (i.e. gene-level quantifications of RNA), on protein abundance for each RNA-protein pair. To do this, we fitted a linear model between RNA and protein abundances across lines (using the Numpy function poly1d in Python), taking the model residuals as the adjusted protein abundance values. Variance decomposition models were then fitted as described above.

All variance component models were fitted using the LIMIX package (https://github.com/PMBio/limix; https://doi.org/10.1101/003905) (Reimand, 2016).

Quantification of X chromosome inactivation (XCI)

Request a detailed protocol

The X chromosome inactivation (XCI) status of female cell lines was quantified using allele-specific counts from RNA-Seq reads mapping to the X chromosome. These allele-specific counts were obtained for SNPs present in DBSNP using GATK ReadCounter with the command ‘GenomeAnalysisTK.jar -T ASEReadCounter -U ALLOW_N_CIGAR_READS --minMappingQuality 10 --minBaseQuality 2’. For known heterozygous SNPs in each line, the allele-specific fraction of expression was defined as the fraction of reads mapping to the less expressed allele (i.e. the allele-specific fraction was ≤0.5). The XCI status of each cell line was then defined as the mean of the allele-specific fractions across all heterozygous X chromosome SNPs with at least 20 overlapping reads in the corresponding RNA-Seq sample.

Gene ontology enrichment was performed against the 6335 genes included in the XCI analysis using g:Profiler (Ongen et al., 2016), and p-values were adjusted for multiple testing using the Benjamini-Hochberg FDR procedure.

QTL mapping of RNA and protein traits

Cis QTL mapping

Request a detailed protocol

We used PEER (Stegle et al., 2012) to account for unwanted variation and confounding factors. PEER was applied to log normalised protein abundance and log normalised gene TPM, considering the most highly expressed 10,000 proteins and genes, respectively. We fit 7 PEER factors for protein and 13 PEER factors for RNA, settings that were determined as the largest number of PEER factors that retain statistical independence of the inferred factors (R < 0.7; Figure 2—figure supplement 1).

For cis genetic analyses, we considered common variants (MAF >5%) in gene-proximal regions of 250 k upstream and downstream of gene transcription start and end sites (GRCh37). The chosen size of this cis analysis window is a compromise between comprehensiveness to detect distal regulatory elements, while managing the multiple testing burden. We used a linear mixed model implemented in LIMIX, thereby controlling for both population structure and repeat lines from the same donor, using kinship as a random effect component. The population structure random effect was accounted for using the realized relationship covariance,that is dot product of the genotype matrices. PEER factors were included as fixed effect covariates.

To adjust for multiple testing across cis variants for each gene, we fit an empirical null distribution using a permutation scheme combined with a parametric fit to the null distribution, similar to the approach taken in Fast QTL (Walter et al., 2015), Briefly, for each gene, we obtained p-values from 100 permutations of the tested variants. We then estimated an empirical null distribution by fitting a parametric Beta distribution to the obtained p-values. Using this null model, we estimated cis region adjusted p-values for QTL lead variants. For multiple testing adjustment across genes, we performed Benjamini-Hochberg adjustment.

For protein, peptide and transcript QTLs, as multiple of these features map to the same gene, we used Bonferroni adjustment to correct for feature multiplicity for each gene, followed by Benjamini-Hochberg adjustment, as performed for the gene-level eQTLs.

Trans QTL mapping

Request a detailed protocol

For trans QTLs analysis, we considered lead cis QTLs (FDR < 10%; 654 pQTLs) versus 11,140 recurrently detected proteins. Genome-wide adjustment for multiple testing was performed using Benjamini Hochberg (BH) across all tests (7⋅106 variants × proteins).

To rule out any potential artefact linked to the mis-mapping of cis protein peptides to the trans proteins, we aligned all trans peptide sequences to the cis protein sequences. We considered all peptides used for the quantification of proteins associated in trans and locally aligned them to the reference sequence of the proteins associated in cis (using pairwise2.align.localxs from the Biopython Project, with a gap penalty of 1). This identified a single trans association for which the peptides had less than two mismatches, which was excluded from the reported results.

Downstream analysis of QTL results

Cis eQTL and pQTL replication

Request a detailed protocol

To assess the replication of QTL across molecular layers, we considered the QTL detected in one layer and assessed the nominal significance in the other layer (PV <0.01), as well as the consistency of the effect directions in the second layer.

To identify the determinants of eQTL to pQTL replication, we trained a Random Forest model to the replication status of 3487 eQTL (Figure 2—source data 2). For each RNA-protein pair, we defined eight features: ‘eQTL effect size’, ‘protein coefficient of error’, ‘protein coefficient of variation’, ‘protein abundance’, ”SNP MAF’, ‘RNA abundance’, ‘number of peptides’, and ‘number of missing measurements’. The feature ‘protein coefficient of error’ was computed as the coefficient of variation of the reference line across the set of technical replicates (TMT batches). This model was fit in Python using the scikit-learn library v0.21.3, with the sklearn.ensembl.RandomForestClassifier model with n_estimators = 100 (i.e. the suggested default). The model was trained and tested 50 times, training on a random sample of 80% of the data, and tested on the remaining 20%.

Annotation of cis-trans protein pairs with protein-protein interactions

Request a detailed protocol

Protein-protein interactions were obtained from the union of CORUM (Ruepp et al., 2010), IntAct (Orchard et al., 2014) and protein-binding interactions from StringDB (Szklarczyk et al., 2017). For CORUM, we considered pairwise interactions between all protein complex subunits, discarded any isoform extension from the protein UniProt IDs, and intersect cis-trans protein pairs with the protein-protein interactions reference list.

Overlap with disease variants and GTEx eQTLs

Request a detailed protocol

Following the approach in Kilpinen et al., 2017, we defined proxy variants of each cis QTL as variants in high LD (r2 >0.8) based on the UK10K European reference panel (Buniello, 2019) and located in the same cis window. A QTL was defined as GWAS-tagging if at least one proxy variant was annotated in the NHGRI-EBI GWAS Catalog (download on 2019–03) (Staley et al., 2016). A pQTL was defined as replicating in GTEx eQTLs if it was mapped at nominal significance (PV <0.01/48) in any of the 48 tissues (Battle et al., 2017).

Complementary, we considered a more stringent criterion based on statistical co-localisation of GWAS signals with QTLs. Specifically, we used summary statistics from phenotypic traits mapped in 51 studies (McLaren et al., 2016), using eCAVIAR (Hormozdiari et al., 2016) to test for co-localisation. For each pQTL cis region and each GWAS trait and study we first intersected the variants assessed both in GWAS and in our study, resulting between 104 and 2 × 105 variants. We then selected the traits with at least one variant with PVGWAS <10−5 and genome wide-significant in our study. For each trait - QTL pair we performed the co-localisation of z-scores (computed as the ratio between effect size and effect size standard deviation).

Annotation of pQTLs with variant effects

Request a detailed protocol

For each pQTL lead variant, a set of proxy variants were identified based on LD (r2 >0.8). Each proxy variant located in the gene body for the corresponding gene was annotated based on its position and predicted effect, using Variant Effect Predictor (Eilbeck et al., 2005). The variants were grouped in parent categories using Sequence Ontology (2016 release; Wagih et al., 2018) as follows: 'inframe_variant' includes variants annotated as 'inframe_deletion', 'inframe_insertion', 'incomplete_terminal_codon_variant', 'stop_lost', 'stop_gained' and 'missense_variant'; 'splicing_variant' includes variants annotated as 'splice_acceptor_variant', 'splice_donor_variant' and 'splice_region_variant'; 'frameshift_variant' includes variants annotated as 'feature_elongation', 'feature_truncation' and 'frameshift_variant'. The set of inframe variants were further classified as deleterious or not according to their SIFT scores (Ng and Henikoff, 2003), as provided by MutFunc 69 (SIFT score <0.05 corresponds to a deleterious mutation). Each pQTL SNP with at least one proxy variant was annotated with the predicted effects of all its proxy variants. For each class of variants (inframe, splicing, etc.), enrichment in the set of pQTLs without any evidence of RNA mechanism (i.e. no eQTL or transcript isoform QTL) compared to all pQTLs was evaluated by Fisher’s exact test, restricted to pQTLs with at least one proxy variant in the gene body (Figure 3d).

References

  1. 1
  2. 2
  3. 3
  4. 4
  5. 5
  6. 6
    The NHGRI-EBI GWAS catalog of published genome-wide association studies, targeted arrays and summary statistics
    1. A Buniello
    (2019)
    Nucleic Acids Research 201947:D1005–D1012.
  7. 7
  8. 8
  9. 9
  10. 10
  11. 11
  12. 12
  13. 13
  14. 14
  15. 15
  16. 16
  17. 17
  18. 18
  19. 19
  20. 20
  21. 21
  22. 22
  23. 23
  24. 24
  25. 25
  26. 26
  27. 27
  28. 28
  29. 29
  30. 30
  31. 31
  32. 32
  33. 33
  34. 34
  35. 35
  36. 36
  37. 37
  38. 38
  39. 39
  40. 40
  41. 41
  42. 42
    Identification of 12 new susceptibility loci for different histotypes of epithelial ovarian Cancer
    1. CM Phelan
    2. KB Kuchenbaecker
    3. JP Tyrer
    4. SP Kar
    5. K Lawrenson
    6. SJ Winham
    7. J Dennis
    8. A Pirie
    9. MJ Riggan
    10. G Chornokur
    11. MA Earp
    12. PC Lyra
    13. JM Lee
    14. S Coetzee
    15. J Beesley
    16. L McGuffog
    17. P Soucy
    18. E Dicks
    19. A Lee
    20. D Barrowdale
    21. J Lecarpentier
    22. G Leslie
    23. CM Aalfs
    24. KKH Aben
    25. M Adams
    26. J Adlard
    27. IL Andrulis
    28. H Anton-Culver
    29. N Antonenkova
    30. G Aravantinos
    31. N Arnold
    32. BK Arun
    33. B Arver
    34. J Azzollini
    35. J Balmaña
    36. SN Banerjee
    37. L Barjhoux
    38. RB Barkardottir
    39. Y Bean
    40. MW Beckmann
    41. A Beeghly-Fadiel
    42. J Benitez
    43. M Bermisheva
    44. MQ Bernardini
    45. MJ Birrer
    46. L Bjorge
    47. A Black
    48. K Blankstein
    49. MJ Blok
    50. C Bodelon
    51. N Bogdanova
    52. A Bojesen
    53. B Bonanni
    54. Å Borg
    55. AR Bradbury
    56. JD Brenton
    57. C Brewer
    58. L Brinton
    59. P Broberg
    60. A Brooks-Wilson
    61. F Bruinsma
    62. J Brunet
    63. B Buecher
    64. R Butzow
    65. SS Buys
    66. T Caldes
    67. MA Caligo
    68. I Campbell
    69. R Cannioto
    70. ME Carney
    71. T Cescon
    72. SB Chan
    73. J Chang-Claude
    74. S Chanock
    75. XQ Chen
    76. YE Chiew
    77. J Chiquette
    78. WK Chung
    79. KBM Claes
    80. T Conner
    81. LS Cook
    82. J Cook
    83. DW Cramer
    84. JM Cunningham
    85. AA D'Aloisio
    86. MB Daly
    87. F Damiola
    88. SD Damirovna
    89. A Dansonka-Mieszkowska
    90. F Dao
    91. R Davidson
    92. A DeFazio
    93. C Delnatte
    94. KF Doheny
    95. O Diez
    96. YC Ding
    97. JA Doherty
    98. SM Domchek
    99. CM Dorfling
    100. T Dörk
    101. L Dossus
    102. M Duran
    103. M Dürst
    104. B Dworniczak
    105. D Eccles
    106. T Edwards
    107. R Eeles
    108. U Eilber
    109. B Ejlertsen
    110. AB Ekici
    111. S Ellis
    112. M Elvira
    113. KH Eng
    114. C Engel
    115. DG Evans
    116. PA Fasching
    117. S Ferguson
    118. SF Ferrer
    119. JM Flanagan
    120. ZC Fogarty
    121. RT Fortner
    122. F Fostira
    123. WD Foulkes
    124. G Fountzilas
    125. BL Fridley
    126. TM Friebel
    127. E Friedman
    128. D Frost
    129. PA Ganz
    130. J Garber
    131. MJ García
    132. V Garcia-Barberan
    133. A Gehrig
    134. A Gentry-Maharaj
    135. AM Gerdes
    136. GG Giles
    137. R Glasspool
    138. G Glendon
    139. AK Godwin
    140. DE Goldgar
    141. T Goranova
    142. M Gore
    143. MH Greene
    144. J Gronwald
    145. S Gruber
    146. E Hahnen
    147. CA Haiman
    148. N Håkansson
    149. U Hamann
    150. TVO Hansen
    151. PA Harrington
    152. HR Harris
    153. J Hauke
    154. A Hein
    155. A Henderson
    156. MAT Hildebrandt
    157. P Hillemanns
    158. S Hodgson
    159. CK Høgdall
    160. E Høgdall
    161. FBL Hogervorst
    162. H Holland
    163. MJ Hooning
    164. K Hosking
    165. RY Huang
    166. PJ Hulick
    167. J Hung
    168. DJ Hunter
    169. DG Huntsman
    170. T Huzarski
    171. EN Imyanitov
    172. C Isaacs
    173. ES Iversen
    174. L Izatt
    175. A Izquierdo
    176. A Jakubowska
    177. P James
    178. R Janavicius
    179. M Jernetz
    180. A Jensen
    181. UB Jensen
    182. EM John
    183. S Johnatty
    184. ME Jones
    185. P Kannisto
    186. BY Karlan
    187. A Karnezis
    188. K Kast
    189. CJ Kennedy
    190. E Khusnutdinova
    191. LA Kiemeney
    192. JI Kiiski
    193. SW Kim
    194. SK Kjaer
    195. M Köbel
    196. RK Kopperud
    197. TA Kruse
    198. J Kupryjanczyk
    199. A Kwong
    200. Y Laitman
    201. D Lambrechts
    202. N Larrañaga
    203. MC Larson
    204. C Lazaro
    205. ND Le
    206. L Le Marchand
    207. JW Lee
    208. SB Lele
    209. A Leminen
    210. D Leroux
    211. J Lester
    212. F Lesueur
    213. DA Levine
    214. D Liang
    215. C Liebrich
    216. J Lilyquist
    217. L Lipworth
    218. J Lissowska
    219. KH Lu
    220. J Lubinński
    221. C Luccarini
    222. L Lundvall
    223. PL Mai
    224. G Mendoza-Fandiño
    225. S Manoukian
    226. L Massuger
    227. T May
    228. S Mazoyer
    229. JN McAlpine
    230. V McGuire
    231. JR McLaughlin
    232. I McNeish
    233. H Meijers-Heijboer
    234. A Meindl
    235. U Menon
    236. AR Mensenkamp
    237. MA Merritt
    238. RL Milne
    239. G Mitchell
    240. F Modugno
    241. J Moes-Sosnowska
    242. M Moffitt
    243. M Montagna
    244. KB Moysich
    245. AM Mulligan
    246. J Musinsky
    247. KL Nathanson
    248. L Nedergaard
    249. RB Ness
    250. SL Neuhausen
    251. H Nevanlinna
    252. D Niederacher
    253. RL Nussbaum
    254. K Odunsi
    255. E Olah
    256. OI Olopade
    257. H Olsson
    258. C Olswold
    259. DM O'Malley
    260. KR Ong
    261. NC Onland-Moret
    262. N Orr
    263. S Orsulic
    264. A Osorio
    265. D Palli
    266. L Papi
    267. TW Park-Simon
    268. J Paul
    269. CL Pearce
    270. IS Pedersen
    271. PHM Peeters
    272. B Peissel
    273. A Peixoto
    274. T Pejovic
    275. LM Pelttari
    276. JB Permuth
    277. P Peterlongo
    278. L Pezzani
    279. G Pfeiler
    280. KA Phillips
    281. M Piedmonte
    282. MC Pike
    283. AM Piskorz
    284. SR Poblete
    285. T Pocza
    286. EM Poole
    287. B Poppe
    288. ME Porteous
    289. F Prieur
    290. D Prokofyeva
    291. E Pugh
    292. MA Pujana
    293. P Pujol
    294. P Radice
    295. J Rantala
    296. C Rappaport-Fuerhauser
    297. G Rennert
    298. K Rhiem
    299. P Rice
    300. A Richardson
    301. M Robson
    302. GC Rodriguez
    303. C Rodríguez-Antona
    304. J Romm
    305. MA Rookus
    306. MA Rossing
    307. JH Rothstein
    308. A Rudolph
    309. IB Runnebaum
    310. HB Salvesen
    311. DP Sandler
    312. MJ Schoemaker
    313. L Senter
    314. VW Setiawan
    315. G Severi
    316. P Sharma
    317. T Shelford
    318. N Siddiqui
    319. LE Side
    320. W Sieh
    321. CF Singer
    322. H Sobol
    323. H Song
    324. MC Southey
    325. AB Spurdle
    326. Z Stadler
    327. D Steinemann
    328. D Stoppa-Lyonnet
    329. LE Sucheston-Campbell
    330. G Sukiennicki
    331. R Sutphen
    332. C Sutter
    333. AJ Swerdlow
    334. CI Szabo
    335. L Szafron
    336. YY Tan
    337. JA Taylor
    338. MK Tea
    339. MR Teixeira
    340. SH Teo
    341. KL Terry
    342. PJ Thompson
    343. LCV Thomsen
    344. DL Thull
    345. L Tihomirova
    346. AV Tinker
    347. M Tischkowitz
    348. S Tognazzo
    349. AE Toland
    350. A Tone
    351. B Trabert
    352. RC Travis
    353. A Trichopoulou
    354. N Tung
    355. SS Tworoger
    356. AM van Altena
    357. D Van Den Berg
    358. AH van der Hout
    359. RB van der Luijt
    360. M Van Heetvelde
    361. E Van Nieuwenhuysen
    362. EJ van Rensburg
    363. A Vanderstichele
    364. R Varon-Mateeva
    365. A Vega
    366. DV Edwards
    367. I Vergote
    368. RA Vierkant
    369. J Vijai
    370. A Vratimos
    371. L Walker
    372. C Walsh
    373. D Wand
    374. S Wang-Gohrke
    375. B Wappenschmidt
    376. PM Webb
    377. CR Weinberg
    378. JN Weitzel
    379. N Wentzensen
    380. AS Whittemore
    381. JT Wijnen
    382. LR Wilkens
    383. A Wolk
    384. M Woo
    385. X Wu
    386. AH Wu
    387. H Yang
    388. D Yannoukakos
    389. A Ziogas
    390. KK Zorn
    391. SA Narod
    392. DF Easton
    393. CI Amos
    394. JM Schildkraut
    395. SJ Ramus
    396. L Ottini
    397. MT Goodman
    398. SK Park
    399. LE Kelemen
    400. HA Risch
    401. M Thomassen
    402. K Offit
    403. J Simard
    404. RK Schmutzler
    405. D Hazelett
    406. AN Monteiro
    407. FJ Couch
    408. A Berchuck
    409. G Chenevix-Trench
    410. EL Goode
    411. TA Sellers
    412. SA Gayther
    413. AC Antoniou
    414. PDP Pharoah
    415. AOCS study group
    416.  EMBRACE Study
    417. GEMO Study Collaborators
    418. HEBON Study
    419. KConFab Investigators
    420. OPAL study group
    (2017)
    Nature Genetics 49:680–691.
    https://doi.org/10.1038/ng.3826
  43. 43
  44. 44
  45. 45
  46. 46
  47. 47
  48. 48
  49. 49
  50. 50
  51. 51
  52. 52
  53. 53
  54. 54
  55. 55
  56. 56
  57. 57
  58. 58
  59. 59
  60. 60
  61. 61
  62. 62
  63. 63
  64. 64
  65. 65
  66. 66
  67. 67
  68. 68

Decision letter

  1. Stephen CJ Parker
    Reviewing Editor; University of Michigan, United States
  2. Patricia J Wittkopp
    Senior Editor; University of Michigan, United States
  3. Arushi Varshney
    Reviewer; University of Michigan, United States
  4. Roderic Guigó
    Reviewer; Center for Genomic Regulation, Spain

In the interests of transparency, eLife publishes the most substantive revision requests and the accompanying author responses.

Acceptance summary:

This exciting study quantifies the genetic architecture of gene expression (eQTL) and protein abundance (pQTL) regulation in iPSCs and will server is an important benchmark for the community. The rigorous comparison of these e/pQTL maps reveals differences that represent expected biological diversity, with a notable example in Figure 3B. Further, this work will serve as a foundation for many novel future studies, including GWAS colocalization and additional omics layers.

Decision letter after peer review:

Thank you for submitting your article "Population-scale proteome variation in human induced pluripotent stem cells" for consideration by eLife. Your article has been reviewed by three peer reviewers, including Stephen CJ Parker as the Reviewing Editor and Reviewer #1, and the evaluation has been overseen by Patricia Wittkopp as the Senior Editor. The following individuals involved in review of your submission have agreed to reveal their identity: Arushi Varshney (Reviewer #3); Roderic Guigó (Reviewer #4).

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.

Summary:

Here the authors report pQTL results on 202 iPSC lines from 151 donors which already had RNA-seq data, enabling direct population-scale comparisons of eQTL and pQTL. This manuscript is well-written, clear, and has impactful results. We think it would be a great resource for the field. This is excellent work, and the authors are to be commended on an exciting study. We have a few items of feedback:

Revisions:

1) It is unclear how the acquisition of RNA and proteomics data were related. For example, were the same batches prepped, split, and material frozen for later profiling. Or, were the lines grown up at separate times for each of the different profiling modalities? If the latter, how could this potentially RNA-protein confounding batch effect influence the results?

2) In subsection “pQTL arising from protein-altering variants”, paragraph three, the authors advocate using protein abundance instead of RNA to interpret pathogenic mechanisms of rare variants. However, QTL studies generally have low power to detect effects for SNPs with low MAF. In fact, in this study, the authors used MAF 5% or higher. So, how does one reconcile this assertion to use protein abundance to interpret rare variant effects with the massive sample size it would take to do so. Addressing this idea in the text would be helpful.

3) The authors should make the full summary scan results available so that the rest of the community can use them as a resource.

4) Discussion paragraph three, the word "significant" is missing a "ly" and we'd advise removing that word altogether. In general, it should be used when making a statistical comparison and the associated test and p-value should be provided. This is not mentioned anywhere in the sentence. So, either those results should be disclosed, or a different word without a statistical connotation should be used. The same issue is present in paragraph three of subsection “pQTL arising from protein-altering variants”.

5) The variance component analysis showed effects of culture medium, sex, age etc. Were these taken as covariates in the QTL analysis? Looking at effects of the culture medium, were the culture passage numbers comparable across lines and does that have an effect?

6) For the X chromosome inactivation (XCI) section, it is unclear how exactly the XCI status was quantified. Subsection “RNA and proteome variability” paragraph two and the legend to Figure 1 reference the Materials and methods section but a description of this analysis is missing there. These results are hard to interpret without methodological clarity.

7) Figure 1D – While the random forest model analysis is interesting, the authors should elaborate on their selection of these specific variables. Some other factors that would be informative to include would be MAF and RNA expression level.

8) Figure 2B and 3C – is the grey bar/shaded region in these plots the genomic location of the respective gene? If so, this should be specified in the legend.

9) Figure 3C and the corresponding text briefly describe a pQTL for the VRK2 gene where the pQTL SNP is also a GWAS SNP for schizophrenia risk. The authors should elaborate on the pQTL direction of effect with respect to GWAS risk. Is the variant the lead SNP for GWAS or what is the LD r2 with the lead SNP and do these signals colocalize in that case? Is there some evidence of this protein being relevant in schizophrenia related cellular mechanisms?

10) Figure 4E – This is an interesting example. What is the eQTL/pQTL and trans pQTL direction of effect with respect to the Alzheimer's GWAS risk allele? In this example, the SNP rs1129187 is associated with PEX6 mRNA expression and protein abundance, and also associated with PEX1 and PEX26 abundance. To directly test if the trans-association of the SNP with PEX1 and PEX26 is through the association with PEX6 (complex stability) and not through other mechanisms, have the authors tried to regress out the PEX6 abundance from the association between the SNP and PEX1 or PEX26 and check if the association disappears?

11) Wrong figures are referenced in some places in the manuscript. Figure 3D is referenced before 3B etc.

12) Were there genes for which significant eQTL and also pQTL associations were identified but the variants were independent (low LD r2?)

13) It was confusing that the authors do not clearly distinguish between the variant affecting the phenotype of the gene (transcript or protein expression) and the affected gene. They write "we report 654 genes with a cis pQTL and 3487 genes with a cis eQTL. I assume that will in general find multiple p/eQTLs for a given gene (althought these number do not appear to be reported). These are thus the numbers of p/eGenes, but not of p/eQTLs. However, when they set to investigate replication of p/eQTLs, the numbers correspond to p/eGenes. The authors equate the numbers of QTLs with the number of affected genes. This part could be more clear.

14) Figure 1B. It looks to me that the fraction of the total variance explained by the factors the authors use in their model is much larger for transcriptomics than for proteomics data. I suggest the authors to report this number. If I am correct, this would mean that proteomics data behaves "somehow" more stochastically than transcriptomics data, maybe reflecting technical issues. It maybe also linked to the lack of replication of eQLTs at the protein level.

15) I understand the rationale of using 250Kb for eQTL analysis, since much of the regulation of gene expression is likely to reside in the promoter region. However, I do not see a biological rationale for using the same window for pQTL analysis. I understand that using the same window maybe the only way of making meaningful comparisons between eQTLs and pQTLs, and I think that this is ok. By using the same window the authors are implicitly assessing to what extent variation affecting gene expression also affects protein expression, that is the genetic variation in which the impact on protein expression is mediated by the impact on gene expression. Maybe the authors should acknowledge this.

16) Related to the above. eQTLs tend to cluster around the TSS. Do they observe the same clustering for pQTLs? What is the comparative distribution of p/eQTLs along the tested region?

https://doi.org/10.7554/eLife.57390.sa1

Author response

Revisions:

1) It is unclear how the acquisition of RNA and proteomics data were related. For example, were the same batches prepped, split, and material frozen for later profiling. Or, were the lines grown up at separate times for each of the different profiling modalities? If the latter, how could this potentially RNA-protein confounding batch effect influence the results?

The lines were expanded once for both banking of material and all molecular assays, including RNA and proteomics quantification. Hence, the same batch of cells was split and frozen material was distributed and used for both RNA and proteomics processing. We have clarified this in the main text (Results section), and included further details about the processing steps in the Materials and methods section.

2) In subsection “pQTL arising from protein-altering variants”, paragraph three, the authors advocate using protein abundance instead of RNA to interpret pathogenic mechanisms of rare variants. However, QTL studies generally have low power to detect effects for SNPs with low MAF. In fact, in this study, the authors used MAF 5% or higher. So, how does one reconcile this assertion to use protein abundance to interpret rare variant effects with the massive sample size it would take to do so. Addressing this idea in the text would be helpful.

We agree that this section was not fully clear. We intended to advocate the value of proteome level data in addition to RNA abundance data and to state that our results demonstrate that the proteomics data add another important level of information for annotating specific pathogenic variants. We agree, however, that the sample size is a concern, which limits the frequency spectrum that can be studied using our data (we consider MAF>5% for our analyses). We have reworded implicated sections accordingly.

3) The authors should make the full summary scan results available so that the rest of the community can use them as a resource.

We agree that these data will be useful. We now provide the complete summary statistics for all eQTL, pQTL, transcript QTL and peptide QTL analyses we have conducted. These high-volume data are not suitable for including as supplementary tables and hence have been deposited on figshare. We have referenced them in the Data availability section.

4) Discussion paragraph three, the word "significant" is missing a "ly" and we'd advise removing that word altogether. In general, it should be used when making a statistical comparison and the associated test and p-value should be provided. This is not mentioned anywhere in the sentence. So, either those results should be disclosed, or a different word without a statistical connotation should be used. The same issue is present in paragraph three of subsection “pQTL arising from protein-altering variants”.

We thank the reviewers for highlighting this issue. We have fixed the corresponding sections and avoid statistical connotation in these specific instances, which we feel are not necessary for these specific statements made.

5) The variance component analysis showed effects of culture medium, sex, age etc. Were these taken as covariates in the QTL analysis? Looking at effects of the culture medium, were the culture passage numbers comparable across lines and does that have an effect?

For the QTL analyses we used PEER factors (Stegle et al., 2012) to adjust for batch and other confounding factors. We have experimented with alternative strategies, including to explicitly account for these additional covariates in the analysis. However, the inclusion of PEER factors performed best in terms of maximizing the number of discoveries. In general, we observed that these factors tend to capture broad covariates such as batch or culture media and hence there tends to be no significant benefit of including them separately. Regarding passage number – we agree that this an interesting covariate, and in particular would have added value to the variance component analysis. Unfortunately, the records on passage number were not fully complete for all lines and hence we could not consider this in our study.

6) For the X chromosome inactivation (XCI) section, it is unclear how exactly the XCI status was quantified. Subsection “RNA and proteome variability” paragraph two and the legend to Figure 1 reference the Materials and methods section but a description of this analysis is missing there. These results are hard to interpret without methodological clarity.

We thank the reviewers and apologize for leaving out this detail from the Materials and methods section. We have now added this in the Materials and methods section.

7) Figure 1D – While the random forest model analysis is interesting, the authors should elaborate on their selection of these specific variables. Some other factors that would be informative to include would be MAF and RNA expression level.

We thank the reviewer for these suggestions. We have extended the analysis and the revised Figure 2D now depicts the relevance of MFA and RNA abundance.

8) Figure 2B and 3C – is the grey bar/shaded region in these plots the genomic location of the respective gene? If so, this should be specified in the legend.

The grey box does indeed indicate the genomic location of the respective gene. We have extended the figure caption accordingly.

9) Figure 3C and the corresponding text briefly describe a pQTL for the VRK2 gene where the pQTL SNP is also a GWAS SNP for schizophrenia risk. The authors should elaborate on the pQTL direction of effect with respect to GWAS risk. Is the variant the lead SNP for GWAS or what is the LD r2 with the lead SNP and do these signals colocalize in that case? Is there some evidence of this protein being relevant in schizophrenia related cellular mechanisms?

Thank you for raising this point. The pQTL lead variant is indeed identical to a reported GWAS risk variant for schizophrenia (Yu et al., 2017). Unfortunately, we could not obtain access to the summary statistics, and hence we have not conducted a formal colocalization test for this locus. However, the effect size direction of this pQTL is consistent with prior evidence for the disease relevance of VRK2. Briefly, the risk allele (OR 1.17) is associated with decreased protein expression. Several previous studies have implicated VRK2 in schizophrenia and have linked downregulation of VRK2 RNA abundance to schizophrenia and other neurological disorders (Azimi et al., 2018, Tesli et al., 2016). We have extended the main text accordingly to mention this (subsection “pQTL arising from protein-altering variants”).

10) Figure 4E – This is an interesting example. What is the eQTL/pQTL and trans pQTL direction of effect with respect to the Alzheimer's GWAS risk allele? In this example, the SNP rs1129187 is associated with PEX6 mRNA expression and protein abundance, and also associated with PEX1 and PEX26 abundance. To directly test if the trans-association of the SNP with PEX1 and PEX26 is through the association with PEX6 (complex stability) and not through other mechanisms, have the authors tried to regress out the PEX6 abundance from the association between the SNP and PEX1 or PEX26 and check if the association disappears?

The risk allele at rs1129187 is associated with increased abundance of both PEX6 RNA and protein level, and is also positively associated with the abundance of the remaining complex subunits PEX26 and PEX1. We have extended the text to include this information.

Regarding the mediating role of PEX6 in the trans association, two lines of evidence indicate that PEX6 is mediating this QTL effect. First, all the complex members are correlated (e.g. r=0.42 for PEX6 and PEX1). Second, we have conducted a regression-based analysis to compare the evidence for a genetic effect on the downstream targets before and after adjusting for PEX6 expression. This results in decreased correlation (r=0.06 vs r=0.29 for PEX1 and r=0.36 versus 0.57 for PEX26), which is consistent with the hypothesized mediation. Nevertheless, this result remains a hypothesis at this stage and hence we have toned down this claim in the main text.

11) Wrong figures are referenced in some places in the manuscript. Figure 3D is referenced before 3B etc.

We thank the reviewer for highlighting these referencing errors. We have carefully revised and checked all references in the paper.

12) Were there genes for which significant eQTL and also pQTL associations were identified but the variants were independent (low LD r2?)

This is indeed a nice addition to our results as presented. In total, we identify 82 genes with significant eQTL and pQTL with independent lead variants (r2<0.1). The LD (r2) between the reported lead variants is now included in the Figure 1—source data 2, 3. We have added a statement to the main text (subsection “Mapping cis genetic effects on protein abundance”).

13) It was confusing that the authors do not clearly distinguish between the variant affecting the phenotype of the gene (transcript or protein expression) and the affected gene. They write "we report 654 genes with a cis pQTL and 3487 genes with a cis eQTL. I assume that will in general find multiple p/eQTLs for a given gene (althought these number do not appear to be reported). These are thus the numbers of p/eGenes, but not of p/eQTLs. However, when they set to investigate replication of p/eQTLs, the numbers correspond to p/eGenes. The authors equate the numbers of QTLs with the number of affected genes. This part could be more clear.

We have carefully revised the manuscript to clarify the number of eGenes versus distinct eQTL/pQTL. Note that given the moderate sample size of our study, we have limited our analysis to lead e/pQTL variants and hence in most cases there is a one to one mapping of these terms. Nevertheless, we agree that correct terminology is important.

14) Figure 1B. It looks to me that the fraction of the total variance explained by the factors the authors use in their model is much larger for transcriptomics than for proteomics data. I suggest the authors to report this number. If I am correct, this would mean that proteomics data behaves "somehow" more stochastically than transcriptomics data, maybe reflecting technical issues. It maybe also linked to the lack of replication of eQLTs at the protein level.

We agree that these differences are quite striking. The most likely explanation is higher assay noise reflecting the lower sensitivity of quantitative proteomics technologies compared to state-of-the-art deep RNA sequencing. This is indeed also the most likely explanation of the reduced mapping power. We have added a note in the main text (subsection “RNA and proteome variability”) as well as the Discussion section.

15) I understand the rationale of using 250Kb for eQTL analysis, since much of the regulation of gene expression is likely to reside in the promoter region. However, I do not see a biological rationale for using the same window for pQTL analysis. I understand that using the same window maybe the only way of making meaningful comparisons between eQTLs and pQTLs, and I think that this is ok. By using the same window the authors are implicitly assessing to what extent variation affecting gene expression also affects protein expression, that is the genetic variation in which the impact on protein expression is mediated by the impact on gene expression. Maybe the authors should acknowledge this.

We agree that these choices deserve some justification. Besides the expected biological mechanisms, e.g. promoter regulation or proximal enhancers, the choice is primarily motivated by tradeoffs in mapping power to identify genetic effects. Larger testing regions entail a multiple testing burden thereby prohibiting the ability to find eQTL/pQTL variants that are proximal to the TSS. There is no generic recipe for these tradeoffs and our choice can be considered a middle ground when comparing to previous pQTL studies (e.g. 20kb around the gene boundaries in Battle et al., 2015 and 1Mb around the TSS in Sun et al., 2018). We acknowledge that a smaller window may mean that we miss interesting distal effects. We have included a brief justification in the Materials and methods section.

16) Related to the above. eQTLs tend to cluster around the TSS. Do they observe the same clustering for pQTLs? What is the comparative distribution of p/eQTLs along the tested region?

Yes, this locational clustering in the vicinity of the TSS is expected, and indeed prior studies have reported this for eQTL (Kilpinen et al., 2017). The distribution of pQTLs is similar to the one observed for eQTLs. While initially we did not report these results, we follow the reviewer's suggestion and now include them (Figure 2—figure supplement figure 3).

https://doi.org/10.7554/eLife.57390.sa2

Article and author information

Author details

  1. Bogdan Andrei Mirauta

    European Molecular Biology Laboratory, European Bioinformatics Institute, Wellcome Genome Campus, Hinxton, United Kingdom
    Contribution
    Conceptualization, Data curation, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing
    Contributed equally with
    Daniel D Seaton and Dalila Bensaddek
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-9733-292X
  2. Daniel D Seaton

    European Molecular Biology Laboratory, European Bioinformatics Institute, Wellcome Genome Campus, Hinxton, United Kingdom
    Present address
    Human Genetics, GSK R&D, Stevenage, SG1 2NY, United Kingdom
    Contribution
    Conceptualization, Software, Formal analysis, Validation, Investigation, Visualization, Methodology, Writing - original draft, Writing - review and editing
    Contributed equally with
    Bogdan Andrei Mirauta and Dalila Bensaddek
    Competing interests
    No competing interests declared
  3. Dalila Bensaddek

    Centre for Gene Regulation & Expression, School of Life Sciences, University of Dundee, Dundee, United Kingdom
    Present address
    Bioscience Core Laboratories, King Abdullah University of Science and Technology, Building 2, Office 2271, Thuwal, 23955-6900, Saudi Arabia
    Contribution
    Resources, Data curation, Investigation, Writing - original draft
    Contributed equally with
    Bogdan Andrei Mirauta and Daniel D Seaton
    Competing interests
    No competing interests declared
  4. Alejandro Brenes

    Centre for Gene Regulation & Expression, School of Life Sciences, University of Dundee, Dundee, United Kingdom
    Contribution
    Resources, Data curation, Investigation, Writing - original draft, Writing - review and editing
    Competing interests
    No competing interests declared
  5. Marc Jan Bonder

    European Molecular Biology Laboratory, European Bioinformatics Institute, Wellcome Genome Campus, Hinxton, United Kingdom
    Contribution
    Data curation, Software, Investigation, Writing - original draft, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8431-3180
  6. Helena Kilpinen

    European Molecular Biology Laboratory, European Bioinformatics Institute, Wellcome Genome Campus, Hinxton, United Kingdom
    Present address
    Helsinki Institute of Life Science, University of Helsinki, Helsinki, Finland
    Contribution
    Investigation, Writing - review and editing
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-6692-6154
  7. HipSci Consortium

    1. European Molecular Biology Laboratory, European Bioinformatics Institute, Wellcome Genome Campus, Hinxton, United Kingdom
    2. Centre for Gene Regulation & Expression, School of Life Sciences, University of Dundee, Dundee, United Kingdom
    3. Wellcome Sanger Institute, Hinxton, Cambridge, United Kingdom
    4. King's College, London, United Kingdom
    Contribution
    Resources, Funding acquisition, Conceptualization
    Competing interests
    No competing interests declared
    1. Chukwuma A Agu
    2. Alex Alderton
    3. Petr Danecek
    4. Rachel Denton
    5. Richard Durbin
    6. Daniel J Gaffney
    7. Angela Goncalves
    8. Reena Halai
    9. Sarah Harper
    10. Christopher M Kirton
    11. Anja Kolb-Kokocinski
    12. Andreas Leha
    13. Shane A McCarthy
    14. Yasin Memari
    15. Minal Patel
    16. Ewan Birney
    17. Francesco Paolo Casale
    18. Laura Clarke
    19. Peter W Harrison
    20. Helena Kilpinen
    21. Ian Streeter
    22. Davide Denovi
    23. Oliver Stegle
    24. Angus I Lamond
    25. Ruta Meleckyte
    26. Natalie Moens
    27. Fiona M Watt
    28. Willem H Ouwehand
    29. Philip Beales
  8. Oliver Stegle

    1. European Molecular Biology Laboratory, European Bioinformatics Institute, Wellcome Genome Campus, Hinxton, United Kingdom
    2. European Molecular Biology Laboratory, Genome Biology Unit, Heidelberg, Germany
    3. Division of Computational Genomics and Systems Genetic, German Cancer Research Center, Heidelberg, Germany
    Contribution
    Conceptualization, Formal analysis, Supervision, Funding acquisition, Writing - original draft, Project administration, Writing - review and editing
    Contributed equally with
    Angus I Lamond
    For correspondence
    oliver.stegle@embl.de
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-8818-7193
  9. Angus I Lamond

    Centre for Gene Regulation & Expression, School of Life Sciences, University of Dundee, Dundee, United Kingdom
    Contribution
    Conceptualization, Resources, Supervision, Funding acquisition, Writing - original draft, Project administration, Writing - review and editing
    Contributed equally with
    Oliver Stegle
    For correspondence
    a.i.lamond@dundee.ac.uk
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0001-6204-6045

Funding

Wellcome Trust Strategic Award and UK Medical Research Council (WT098503)

  • Bogdan Andrei Mirauta
  • Dalila Bensaddek
  • Helena Kilpinen

Wellcome Trust Strategic Award (105024/Z/14/Z)

  • Bogdan Andrei Mirauta

EMBL Interdisciplinary Postdoctoral (EIPOD) programme under Marie Sklodowska-Curie Actions COFUND (grant number 291772)

  • Daniel D Seaton
  • Marc Jan Bonder

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Senior Editor

  1. Patricia J Wittkopp, University of Michigan, United States

Reviewing Editor

  1. Stephen CJ Parker, University of Michigan, United States

Reviewers

  1. Arushi Varshney, University of Michigan, United States
  2. Roderic Guigó, Center for Genomic Regulation, Spain

Publication history

  1. Received: March 30, 2020
  2. Accepted: August 8, 2020
  3. Accepted Manuscript published: August 10, 2020 (version 1)
  4. Accepted Manuscript updated: August 12, 2020 (version 2)
  5. Version of Record published: August 25, 2020 (version 3)

Copyright

© 2020, Mirauta 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.

Metrics

  • 778
    Page views
  • 107
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Download citations (links to download the citations from this article in formats compatible with various reference manager tools)

Open citations (links to open the citations from this article in various online reference manager services)

Further reading

    1. Developmental Biology
    2. Genetics and Genomics
    Meltem Weger et al.
    Research Article

    The glucose-sensing Mondo pathway regulates expression of metabolic genes in mammals. Here, we characterized its function in the zebrafish and revealed an unexpected role of this pathway in vertebrate embryonic development. We showed that knockdown of mondoa impaired the early morphogenetic movement of epiboly in zebrafish embryos and caused microtubule defects. Expression of genes in the terpenoid backbone and sterol biosynthesis pathways upstream of pregnenolone synthesis was coordinately downregulated in these embryos, including the most downregulated gene nsdhl. Loss of Nsdhl function likewise impaired epiboly, similar to MondoA loss of function. Both epiboly and microtubule defects were partially restored by pregnenolone treatment. Maternal-zygotic mutants of mondoa showed perturbed epiboly with low penetrance and compensatory changes in the expression of terpenoid/sterol/steroid metabolism genes. Collectively, our results show a novel role for MondoA in the regulation of early vertebrate development, connecting glucose, cholesterol and steroid hormone metabolism with early embryonic cell movements.

    1. Genetics and Genomics
    2. Microbiology and Infectious Disease
    Jeziel Dener Damasceno et al.
    Research Article Updated

    DNA replication is needed to duplicate a cell’s genome in S phase and segregate it during cell division. Previous work in Leishmania detected DNA replication initiation at just a single region in each chromosome, an organisation predicted to be insufficient for complete genome duplication within S phase. Here, we show that acetylated histone H3 (AcH3), base J and a kinetochore factor co-localise in each chromosome at only a single locus, which corresponds with previously mapped DNA replication initiation regions and is demarcated by localised G/T skew and G4 patterns. In addition, we describe previously undetected subtelomeric DNA replication in G2/M and G1-phase-enriched cells. Finally, we show that subtelomeric DNA replication, unlike chromosome-internal DNA replication, is sensitive to hydroxyurea and dependent on 9-1-1 activity. These findings indicate that Leishmania’s genome duplication programme employs subtelomeric DNA replication initiation, possibly extending beyond S phase, to support predominantly chromosome-internal DNA replication initiation within S phase.