Plasma proteomic associations with genetics and health in the UK Biobank

The Pharma Proteomics Project is a precompetitive biopharmaceutical consortium characterizing the plasma proteomic profiles of 54,219 UK Biobank participants. Here we provide a detailed summary of this initiative, including technical and biological validations, insights into proteomic disease signatures, and prediction modelling for various demographic and health indicators. We present comprehensive protein quantitative trait locus (pQTL) mapping of 2,923 proteins that identifies 14,287 primary genetic associations, of which 81% are previously undescribed, alongside ancestry-specific pQTL mapping in non-European individuals. The study provides an updated characterization of the genetic architecture of the plasma proteome, contextualized with projected pQTL discovery rates as sample sizes and proteomic assay coverages increase over time. We offer extensive insights into trans pQTLs across multiple biological domains, highlight genetic influences on ligand–receptor interactions and pathway perturbations across a diverse collection of cytokines and complement networks, and illustrate long-range epistatic effects of ABO blood group and FUT2 secretor status on proteins with gastrointestinal tissue-enriched expression. We demonstrate the utility of these data for drug discovery by extending the genetic proxied effects of protein targets, such as PCSK9, on additional endpoints, and disentangle specific genes and proteins perturbed at loci associated with COVID-19 susceptibility. This public–private partnership provides the scientific community with an open-access proteomics resource of considerable breadth and depth to help to elucidate the biological mechanisms underlying proteo-genomic discoveries and accelerate the development of biomarkers, predictive models and therapeutics1.

each rack in storage. 6,000 samples were preselected by consortium members; the remaining samples were selected from racks visited by consortium picking, with stratified sampling used to ensure a representative distribution through the UK Biobank cohort. Additionally, UK Biobank apply a pseudo-randomisation algorithm to maximise variation of phenotypes on each generated sample plate. Typically, 4 racks (24 samples per rack) are used to generate a rack of 96 picked samples. The algorithm works to ensure the racks are selected in an order that avoids clustering of racks from the same timepoint/sample collection clinic. This in turn results in pseudo-randomisation of other phenotypes like age and gender. Day of week of collection, selfreported participant ethnicity* and deprivation index for selected participants were confirmed as representative of cohort distributions. This created a total selection of 50,002 samples. For the second picking process a total of 7,000 samples were supplied. First, 1,020 pre-selected samples from the Consortium members were chosen. A further 3,637 samples were selected from participants attending the COVID-19 case-control imaging study. 1,270 participants formed this section of the study, who had been invited to attend UK Biobank for a repeat imaging assessment on the basis of a COVID-19 diagnosis from linked healthcare records or a positive home-based SARS-CoV-2 antibody lateral flow test provided by UKB, or were invited as matched controls of these COVID-19 cases. Full inclusion criteria are described in https://biobank.ndph.ox.ac.uk/showcase/showcase/docs/casecontrol_covidimaging.pdf.
For these 1,270 participants, samples from initial recruitment, pre-COVID and post-COVID imaging visits were included where possible. An additional 2,343 baseline samples were included in this picking process selected randomly as described for the first sample selection, optimised to include selection locations required for the COVID imaging and consortiumselected samples.
*We note that self-reported participant ethnicity was collected at each UKB Assessment Centre via touchscreen questionnaires, asking participants, "What is your ethnic group?" and "What is your ethnic background?" These questions were dropped from the touchscreen protocol on October 24 th , 2016, as noted in UKB Data Field 21000. We have retained the terms "ethnicity" and "ethnic background" when referring to the touchscreen questionnaire responses; however, all primary analyses described in this manuscript refer to genetic ancestry.

UKB sample handling
EDTA (9ml) vacutainers were collected and fractioned to 850µl aliquots of EDTA plasma, buffy coat and red cells as per the UK Biobank collection protocol described in 2 . Aliquots were stored in an automated -80°C sample archive. EDTA plasma from participants selected for the study were withdrawn from the automated sample archive in a quasi-randomised order 1 . Prior to processing, aliquots were stored in a -80°C freezer. Racks containing 85 aliquots were thawed (column 12 plus one additional random well left empty for controls and two spaces for blinded duplicates, Figure S2) and 60µl of EDTA plasma was pipetted to a PCR plate (P/N AB0800, Thermo Scientific™) using a TECAN freedomEVO with full sample tracking. Two blinded duplicates (using aliquots from the same plate) were added to each plate before plates were sealed with adhesive seals (P/N 4306311, Applied Biosystems™) and stored at -80°C.
Samples were shipped on dry ice to Olink Analysis Service in Sweden for analysis accompanied by an electronic sample manifest containing sample level information.

Plasma profiling using the Olink technology
The Olink technology uses Proximity Extension Assay, where a matched pair of antibodies labelled with unique complimentary oligonucleotides (proximity probes) bind to the respective target protein in a sample. As a result, the probes come into close proximity and hybridize to each one, enabling DNA amplification of the protein signal, which is quantified on a next generation sequencing read-out which has been described in detail previously 3,4 . Antibodies Samples were serially diluted to 1:1, 1:10, 1:100, 1:1000, 1:10000, transferred to 384-well plates consisting of four abundance blocks for each of the eight panels per 96 samples. In the immune-reaction, plasma samples were incubated overnight at 4℃ with the proximity probes.
Oligonucleotides in close proximity are extended and amplified using DNA polymerase and create a DNA sequence which is amplified by polymerase chain reaction (PCR 1), to create amplicons encoding protein assay information. All amplicons for each sample from the four abundance groups per panel are combined, resulting in one well of amplicons for each sample/panel. Four unique 96-index plates are added to every 4 sample plates, followed by a second PCR (PCR 2), to enable all samples in a plate to be combined in to one library per panel.
Each library is bead purified and libraries quality controlled using a Bioanalyzer. Four sample plates from four identical panels were pooled, denatured and sequenced on individual lanes on the Novaseq600 using S4 flow cells v1.5 (35 cycles) and 384-samples and 384-assays measured. Counts of known sequences were translated in to Normalized Protein eXpression (NPX) values within Olink's MyData Cloud Software.

Olink inbuilt assay quality control
The Olink workflow has an inbuilt quality control system, 3 engineered internal controls that are spiked into every sample, and each abundance block. The incubation control (Inc Ctrl) is a non-human assay, green fluorescent protein (GFP) and is used during data QC. The extension control consists of two paired oligonucleotides coupled to an antibody molecule with the DNA arms in close proximity and was used for normalization of the data. The amplification control (Amp Ctrl) consists of a synthetic double stranded DNA template, used for QC and to monitor the PCR steps in the protocol. In addition, each sample plate includes external controls in column 12. A negative control ran in triplicate is used to calculate the limit of detection (LOD) of each assay in every plate, and a plate control sample, consisting of a pooled plasma sample is run in triplicate to adjust the levels between plates. A duplicate pooled sample control is included to estimate precision within and between runs. In each of the four panels Cardiometabolic, Inflammation, Neurology and Oncology, 3 overlapping assays of IL6, IL8 (CXCL8), and TNF are included; in each of the four panels Cardiometabolic II, Inflammation II, Neurology II and Oncology II, 3 overlapping assays of LMOD1, SCRIB and IDO1 are included. Overlapping assays are used for quality control (QC) purposes to assess correlation.
Olink's internal QC assessment is performed at two levels; run QC and sample QC. For run QC, each abundance block per panel and sample plate should fulfil the mean absolute deviation (MAD) in both internal controls (Inc Ctrl and Amp Ctrl) which should not exceed 0.3 NPX, the deviation of sample QC level is allowed for up to 1/6 samples and in each panel the median of 90% assays in plate and negative controls should be in the accepted range from predefined values set during validation. The sample QC assesses all samples individually using the internal controls (Inc Ctrl and Amp Ctrl) which should be within ± 0.3 NPX from the plate median across the abundance block, in addition to the mean assay count for a sample may not be less 500 counts. Samples that do not fulfil these criteria will receive a sample warning for a given abundance block in the data set. Data from these assays should be treated with caution. Assays where the median of the negative Control triplicates deviate more than 5 SDs from predefined values set for each assay during validation will receive a QC warning in the results NPX file.
Data was generated according to Olink's standard procedures. Normalized Protein eXpression (NPX) is Olink's relative quantification unit on a log-2 scale. Data generation consists of normalization of matched counts of an assay to the extension control spiked into every sample, log-2 transformation of the data, and level adjustment using the plate control.

NPX calculation and normalization
Samples from the study were divided into two sets: i) Set 1 -UKB; and ii) Set 2 -COVID; depending on the time point they were randomly selected from the UKB population. Samples were randomly assigned to 96-well plates, and fully randomized within plates. Each plate contained: i) 87 samples from set 1 or set 2; ii) 1 empty well; iii) 2 Olink control samples used for quality control; iv) 3 Olink negative control samples used to compute the baseline assay level of each plate; and v) 3 Olink plate control samples used for normalization of protein expression ( Figure S2). All Olink samples were in column 12 of each plate while the remaining 87 samples + 1 empty well were randomized across columns 1-11 and rows A to H.  Calculation of Normalized Protein eXpression (NPX) values was performed stepwise as described below. Initially, we calculated the log2 ratio of counts of each assay of each sample to the counts of the extension control, and the assay-specific median value of the plate controls was subtracted. This provided us with plate normalized NPX values for both sets. For samples in batches 0-6, we subtracted the assay-specific plate median NPX value. 3 assays (PNLIPRP2, TDGF1 and FOLR3) with proven bimodal distribution driven by genotypes 4,5 were excluded from the last step, hence, remained plate-normalized. At this stage data was normalized within each batch. Next, we computed adjustment factors from the difference of the assay specific median NPX value of each batch to the reference batch (batch 1). The selection of the reference batch does not impact coefficients of variation as the across-batches normalization uses differences of assay-specific median NPX values between batches. These are fixed factors in NPX scale that aim at harmonizing median NPX values of assays between batches. Finally, adjustment factors from above were added to the NPX values of each batch of set 1. In summary, set 1 was normalized using a two-step approach of within-batch and across-batches intensity normalization. Samples of set 2 were normalized using reference (bridge) samples that were shared between the two sets. All plates that contained at least one sample of set 2 were assigned a randomly selected sample from the set 1. 93 samples with missing frequency <10% and representative of the NPX dynamic range were selected from batches 1-6 in set 1.
These samples were assigned to empty wells of the 93 plates from batch 7 of set 2. In this case, adjustment factors were computed from the assay-specific median of the pair-wise differences between set 1 and set 2. Adjustment factors were added to the NPX values of set 2. Intensity normalized NPX values for set 1 and bridge normalized NPX values for set 2 consisted of the final set of NPX values that was used for downstream analysis.
Within batch normalization centers data at NPX=0 by subtracting the plate-specific median per assay from all samples and assays in the same plate. Across batches normalization computes a set of adjustment factors as the difference of the assay-specific median NPX values of each batch. The first step accounts for plate-to-plate variation within batch and the second step accounts for batch-to-batch variations within the study. Both steps are shifts of an assay-specific fixed factor in NPX scale; namely, plate median in the first step and difference of assay-specific medians between batches in the second step. Normalization steps do not affect intra-plate CVs as the exact same factors are applied in both steps. Inter-plate CVs after the first step of normalization improve within-batch CV, while inter-plate CVs after the two-step normalization improve CV for the full dataset.

Data pre-processing and quality checking
The raw UKB-Olink data contained 58,776 samples and 54,221 individuals; taking out control samples, unprocessed samples and participants who have withdrawn from the study reduced the data to 58,353 samples and 54,219 individuals.
Outliers were identified by two approaches applied to each panel of proteins: (1) principal component analysis (PCA), and (2) examining the median and IQR of NPX across proteins by sample. The data used for PCA were the NPX data excluding: (1) control samples, (2) those who withdrew from the study/data were not processed, (3) data points missing NPX values, (4) those missing in covariates such as sex and sampling center. We removed data points with (1) a standardized PC1 (the component that captures the most variation) or PC2 (second largest component) value more than 5 standard deviations from the mean (which is zero in standardized PCA), or (2) a median NPX greater than 5 standard deviations from the mean median, or an IQR of NPX greater than 5 standard deviations from the mean IQR. We excluded outliers, data points with a QC or assay warning, as well as likely sample swaps where we excluded the sample across all panels if half or more panels were affected; the remaining data contained 56,695 samples and 52,790 individuals. Suspected sample swaps were identified using inconsistency from proteomic predicted sex and outliers from cis pQTLs where we summed the standardized squared residuals over all proteins for each individual and divided by the sum of squared protein levels -samples with the wrong genotype should have larger values than those with correct genotypes. We excluded 1 protein with high degree of QC failures (GLIPR1) from our main analyses. We did not remove data based on LOD (lower limit of detection). We did not perform further processing of the resulting NPX data. The attrition of data is summarized in Supplementary

Protein coefficients of variation (CVs)
Two sets of duplicate samples were provided both by Olink Proteomics ('Olink controls') and UK Biobank ('Blind Spike Duplications' / BSDs). We measured the intra-person variability of protein of individual on plate for duplicate samples by calculating the coefficient of variation (CV) of NPX, as recommended by Olink: Then, to summarize the CV for each protein, we took the median of the CVs across individuals and plates by protein, plotting the median CVs for each protein: Intra-person CVs of each Olink protein are illustrated in Figure S4 and shows that the CVs of the proteins range between 1.8% to 27.2%.

Detecting batch effects
To detect potential batch effects, we fit a simple random-effects model on each protein and calculated the percentage of variability attributable to batch. A high percentage suggested greater variability in NPX across batch. The model used was: The percentage of variability attributable to batch, by protein, is illustrated in Figure S5.
Overall, the proportion of variability attributable to batch was low. No notable evidence of batch effects was observed.

Detecting plate effects
We calculated two sets of CVs for each plate, plotting one series of CVs against the other, and identified whether any plate had high CVs on either axis of the plots. The two sets of CVs were defined as follows: Approach Inter-person versus inter-replicate variability for all plates is illustrated in Figure S6, and there is no obvious plate whose CV is out of range of others. Thus, we conclude that no evidence of plate effect is observed.

Figure S6. Plotting CV measuring intra-person vs inter-replicate variability for all plates.
We performed another analysis to examine potential plate effects as we did for batch effects; that is, we fit a simple random-effects model on each protein and calculated the percentage of variability attributable to plate. A high percentage suggested greater variability in NPX across plate, indicating the presence of plate effect. All except those labeled in the plot had inter-class correlation (percentage of variability attributable to plate) below 10%, and we did not think there was a systematic plate effect ( Figure S7). Figure S7. Percentage of variability attributable to plate by protein.

Proteins below limit of detection (LOD)
The proportion of samples with measurements below LOD for proteins across the panel are summarised in Figure S8. Most of the protein analytes (2,154 out of 2,941) have the proportion of samples below LOD lower than 30%. 1,981 proteins have the proportion ranged between 0-10%; 92 and 81 proteins ranged between 10-20% and 20-30%, respectively. 787 proteins have 30% or more data falling below the lowest level of detection ( Figure S8a).
We also examined how many proteins had less than 5% of samples below LOD, overall and  Figure S8b).

Global proteomic analysis across sex, ancestry and technical factors
In addition to the analyses above, we also performed PCA analysis on the global proteomics data and investigated for any systematic clustering due to sex, ancestry and sampling centre.
We also visualize sex and ancestry effects by UMAP to investigate any extreme hidden clusters.

Sex effects
We plotted the global PCA and UMAP results, coloring each data point by sex of the individual; each data point represented a sample. We did not observe clustering of data points of either sex; no global sex effect was indicated.

Ancestry effects
The fact that there were six ancestry groups, in which those of European origin accounted for a vast majority of data points, made it difficult to detect any clustering visually. Thus, we separated those of European heritage from others and made two PCA and UMAP plots. We did not observe any clustering of data points in either plot indicating minimal ancestry effect.

Sampling center effects
To examine potential sampling center effects, plotting the PCA, we examined the distribution of first principal component (PC1) and second principal component (PC2) of each center ( Figure S11). The median of PC1 and 2 were within 2SD of the means of all centre medians for PC1 and 2. We also took an approach similar to as for batch effect investigation, i.e., fit a simple random-effects model to examine the proportion of variation attributable to sampling center: For PC1, the proportion of variability attributable to sampling center was 3.5%, and the corresponding figures for PC2 was 1.0%. Because the proportion of variability attributable to sampling center was low, we do not think there was strong evidence indicating the existence of sampling center effect.
Figure S11: Distribution of PCs 1 and 2 (the components that explain the largest share of variation) by UKB sampling centers. Sample numbers for each center are at the top. Each box plot presents the median, first and third quartiles, with upper and lower whiskers representing 1.5x inter-quartile range above and below the third and first quartiles respectively.
We did not see a particular center that is out of range of others, and hence we did not believe a sampling center effect was indicated.

Sample age effect
Sample age is defined as the difference between the time blood was taken and the protein measurement. Because sample age is a continuous variable (measured in days in its original scale), it may not be very effective to make a PCA plot to detect potential clustering of data as we did for sex or sampling age. Instead, we examine potential "sample age effect" using an approach like batch effect; that is, we fit a simple random-effects model on each protein and calculated the percentage of variability attributable to age of sample Overall the plots show that for most proteins, the percentage of variability in NPX attributable to age of sample is below < 10% except those labeled on the plots (Figure S12, Supplementary Table 3). Thus, we do not think a "sample age effect" is indicated globally.
We found higher sample age contributions to protein variabilities are also correlated to lower proportion of proteins with measures below lower limit of detection (Spearman's r=-0.18, p=6.8x10 -22 ), suggesting no evidence of sample-age susceptible proteins being less reliable or abundant than more sample-age stable ones. Storage time effects may vary protein to protein as each protein has its own unique stability, thus care is still needed to account for these effects in corresponding analyses.

Comparison of Olink proteins with independent assays in UKB
We compared seven protein measures acquired using the antibody-based Olink assay with seven corresponding protein measurements obtained using independent assays in UKB,  Table 4). Two sets of 3 proteins (CXCL8, IL6, TNF and IDO1, LMOD1, SCRIB) were assayed across   four protein panels: CXCL8, IL6, TNF on the first set of Olink protein panels and IDO1, LMOD1, SCRIB on the version II of the respective protein panels (Supplementary Table 3).

Correlation between same proteins measured across panels
We observed reasonably strong correlations between measurements across panels for each of the 3 proteins measured on each set of four protein panels (Extended Data Figure 2a)

Proteomic associations with age, sex and BMI
In total, we found 1,944, 2,092 and 2,348 associations between protein levels and age, sex, and   Figure 3b).
We also explored interaction effects between age, sex and BMI on protein levels in the same model. In total, we found 40 proteins levels with evidence of significant interactions (p<1.7x10 -5 ) between age, sex and BMI; 1,936 between age and sex; 677 between sex and BMI; and 828 between age and BMI (Supplementary Table 6). For example, we found the strongest interactions between age and sex for follitropin subunit beta (FSHB, p=2.6x10 -1113 ) and glycodelin (also known as progesterone-associated endometrial protein, PAEP, p= 6.5x10 -1421 ).
Follitropin is a key hormone in female reproductive health; the observed increase in follitropin levels at menopausal age in women only (Figure 1d) is a well-established hormonal change 16 .
Glycodelin is a glycoprotein expressed in mammary glands and endometrial tissues 17 . Levels of glycodelin decreased with age for women only, particularly before the age of menopause (~50 years), whilst for men, levels minimally increased with age (Figure 1d). After 55 years of age, levels of glycodelin minimally increased in women at a similar rate to men. These effects are consistent with the role of glycodelin in female reproductive tissues and their associated changes in hormone levels (such as progesterone) around menopause 17 -demonstrating that the proteomic assay used in this cohort can capture longitudinal physiological effects from cross-sectional measurements .

Proteomic associations with health burden and prevalent diseases
We investigated the effects of smoking, medication use and the top 20 most prevalent conditions on protein levels (Supplementary Table 2  We note that all proteomic associations with health burden and disease have not been adjusted for the potentially complex influences of disease severity or time between diagnosis and blood collection; thus, these observational findings should be interpreted carefully, and downstream validation studies are strongly recommended.

Proteomic associations with renal and liver function markers
We investigated proteomic associations with estimated glomerular filtration rate (eGFR) -a measure of renal function, and liver function enzymes alanine transaminase (ALT) and aspartate transaminase (AST) (Methods). We identified 1,815 proteins significantly associated with eGFR (p<1.7x10 -5 ). As anticipated, CST3 was most strongly associated with eGFR, (Supplementary Table 7), owing to its use in the calculation of eGFR values 20 via the CST3 immuno-turbidimetric assay (Supplementary Table 4). We also found 16 assayed proteins that have previously been associated with end stage renal disease 21 to be significantly associated with reduced eGFR (max p=1.8x10 -17 ). For liver function, 2,016 and 1,843 proteins were associated with ALT and AST, respectively, and 18 of the top 20 associations with ALT/AST were enzymes (Supplementary Table 7). The most strongly associated nonenzymatic protein for both ALT and AST was KRT18, which is abundantly expressed in the liver and released in liver cell death, and serves a biomarker for liver disease 22 . We found that We trained LASSO models with 10-fold cross-validation to determine whether plasma proteomics can be used to infer demographics (age, sex, BMI), renal (eGFR) and liver (ALT, AST) functions as well as blood groups (Methods). Evaluating the models on 20% random held-out data, we found that proteomic data alone provided good predictive performance

Consequences of high impact pQTLs
Following annotation of the primary pQTLs, we identified 37 cis pQTLs annotated as potential high-impact variants, with 33 (89%) associating with lower cis protein levels, and all 15 of the primary cis pQTLs coding for start codon lost/stop codon gained, associating with lower corresponding protein levels (Supplementary Table 12). 24 trans pQTLs SNPs were also annotated as potential high-impact (Supplementary Table 12

List of previous pQTL studies
To evaluate whether the pQTLs in the discovery set were novel, we used a list of published GWAS with proteomics (http://www.metabolomix.com/a- studies were excluded if they tested or reported SNPs in candidate genes or only tested/reported SNPs in cis genes. 4. Appropriate p-value adjustment (genome-wide threshold) was included and reported; no false discovery rate corrections or studies that fail to report p-value adjustment methods.
5. Phenotypes tested must be protein based: Studies from flow cytometry or other cell imaging modalities were excluded.
6. Studies must report all proteins in a panel. Studies that generated proteomic data on a platform but only reported selected proteins from that panel were not included.
Following this curation, thirty-four studies were included (Supplementary Table 32). Using a p-value threshold of 1.7x10 -11 , we identified the sentinel variants and associated protein(s) from the previously published studies and queried those against our discovery pQTLs. If a previously associated sentinel variant-protein pair fell within a 1 Mb window of the discovery set pQTL sentinel variant for the same protein and had an r 2 >0.8 with any significant SNPs in the region, it was considered a replication.

Identification, fine mapping, and colocalization of independent signalsfurther details
We identified 29,420 independent pQTL signals with SuSiE regression of individual-level protein levels on genotype dosages and confirmed statistical independence through multiple linear regression models (Supplementary Table 16 We also performed fine-mapping with SuSiE to narrow down 95% credible sets of causal variants for each pQTL signal (Supplementary Table 16  To explore whether the cis pQTL signals for a given protein also had cis or trans effects on other proteins, we performed pairwise colocalization analyses using the coloc with SuSiE framework. Of 10,750 independent cis signals, 11% (1,233) were colocalized with at least one signal from another protein (Supplementary Table 18 However, we also observed instances where colocalized cis pQTLs had effects on two nearby genes that otherwise appeared unrelated such as SPP1/PKD2, FURIN/FES and TIMP1/CFP, which could be explained by non-coding variants in enhancers that regulate more than one gene.

Trends of pQTL associations with increasing sample size and proteins assayed
We observed an initial increase in detectable cis pQTLs at sample sizes below 5,000 before slowly plateauing as the number of cis pQTLs trended towards the number of proteins tested - Mean variance explained by cis associations quickly plateaued beyond samples sizes >5,000 whilst the mean variance explained by trans associations continued to slowly increase and drive most of the increase in mean total variance explained at sample sizes >5,000 (Figure 2f).
We found the mean proportion of variance explained by independent pQTLs increased the most at sample sizes less than 5,000 (Figure 2f). Mean variance explained by cis associations quickly plateaued beyond samples sizes >5,000 whilst the mean variance explained by trans associations continued to slowly increase and drive most of the increase in mean variance explained at sample sizes >5,000 (Figure 2f).
Overall, the rate of increase in the number of pQTL associations with increasing proteins measured slowed beyond the most abundant 1,000 assayed proteins (Methods, Figure 2g).
This is consistent with more lower abundance proteins being increasingly sampled and the reduced pQTL detection at lower dilutions (lower expected abundance) observed in this study.
Although the yield of new pQTLs is starting to decrease as newer larger assay platforms begin to measure more low abundant proteins (% of proteins in 1:1 dilution section = 76% in Olink II panels vs 57% in Olink I panels, Supplementary Table 3), we still anticipate considerable additional associations to be gained, especially at large sample sizes, beyond the proteins measured here (Figure 2g), before the saturation point of all detectable proteins present in plasma.
We also found a shift towards an increasing number of genomic regions harboring associations with multiple proteins with larger sample sizes, indicating greater detectability of pleiotropic loci at increased study sizes (Extended Data Figure 9a). Furthermore, we found a slightly sublinear increase in the number of proteins with at least one interacting protein in any trans loci as sample size increased (Extended Data Figure 9b) -suggesting additional trans target interacting loci for other proteins can be found with larger studies. However, the proportion of trans loci containing at least one interacting protein with the tested protein decreased at a slowing rate with sample size (Extended Data Figure 9c) -suggesting increased detections of trans loci are not driven by direct protein-protein interactions.

Sensitivity analyses of pQTLs
We also explored, a priori, the impact of blood cell composition, BMI, seasonal and fasting time before blood collection on pQTL effects (Supplementary Table 24, Figure S13).

Effects of blood cell counts
Most primary associations from the discovery analysis (84.0% [12,007/14,287], including 99.4% significantly associated with 444 proteins in the discovery analysis, and 95 proteins in the blood-cell sensitivity analysis. The ARHGEF3 locus is established to be highly pleiotropic and known to associate with platelet counts 29,30 . A previous plasma pQTL study suggested that the observed pleiotropy at ARHGEF3 may be driven by genetically determined increases in platelet counts and related sequelae that cause proteins to be secreted into plasma during sample handling and preparation 29 . We further tested this hypothesis through a formal mediation analysis of 349 variant -protein associations for platelet counts at the ARHGEF3 locus using individual participant data (Methods). After correction for multiple testing (Bonferroni correction for 349 variant -protein associations and 9 blood cell phenotypes; p=1.25x10 -5 ), 79.08% of the associations were found to be fully mediated by platelet counts, in line with the observed decrease in pleiotropy at ARHGEF3 variant rs1354034 after adjusting for blood-cell composition.
We cross-referenced the 2,415 proteins with pQTLs for enrichment in various tissues including blood reported by Uhlen et.al 31 (Supplementary Table 24). Within the list of 2,415 proteins, 3 proteins (CLEC4C, IFNL1 and MAP1LC3B2) were enriched in blood and 34 other proteins were enriched in more than one tissue including blood (Supplementary Table 33). The list of these proteins along with the tissues and the specific distributions are detailed in Supplementary Table 24.

Effects of BMI
In the BMI-adjusted analysis, all but 1.8% [264/14,287; 259 trans associations, 5 cis associations] of the primary associations from the discovery analysis remained significant (p<1.7x10 -11 ). Of these 259 trans associations, only one association (leptin levels -rs56094641) fell below a significance threshold of p<0.05 (p=0.15). This association with leptin levels was driven by an intronic variant (rs56094641) in FTO, an established obesity associated locus, suggesting that leptin association with this variant is mediated by obesity 32 .

Effects of season and amount of time fasted at blood collection
The majority (99.4% [14,207/14,287]) of sentinel pQTLs identified by the discovery analysis remained genome-wide significant (p<1.7×10 -11 ) after adjustment for season and participantreported fasting time at blood collection. P-values for the 80 associations (2 cis, 78 trans pQTLs) that were no longer genome-wide significant after accounting for multiple testing ranged from 5.0×10 -11 to 1.7×10 -11 , suggesting minimal impact of season and/or fasting time on variant associations with protein levels.

Proteomic variance contributions by blood cell composition, BMI, seasonal and fasting time
We additionally explored the effects of the aforementioned factors along with demographics contributing to the proteomic variance explained (Figure S14, Supplementary Table 25). We first created a reference model regressing the demographic covariates used in GWAS (age, gender, age*gender, age 2 , age2*gender, ancestry and 20 genetic PCs) on protein levels. We then calculated the difference in the variance explained by the reference model and a multivariable model including the base covariates and each of the non-genetic factors separately to evaluate the percentage of variance in protein levels explained by the non-genetic factors alone. On average, demographic covariates explain 3.7% of the variance in plasma protein levels, while blood cell covariates, BMI, fasting time, and season explain 4%, 1.5%, 0.06%, and 0.03% respectively. Figure S13. Comparison of effect size and p-value changes before and after sensitivity analyses. P-values derived from REGENIE regression GWAS (two-sided, unadjusted).

Assay Target Specificity category Enhanced tissues CLEC1B
Group enriched blood,liver CLEC4C Tissue enriched blood CLEC4D Group enriched blood,bone marrow,lymphoid tissue CLEC6A Group enriched blood,heart muscle,lung,lymphoid tissue CSF2 Tissue enhanced blood,lung,pancreas DRAXIN Tissue enhanced blood,brain FCER2 Group enriched blood,lymphoid tissue FCRL1 Group enriched blood,lymphoid tissue FCRL2 Group enriched blood,intestine,lymphoid tissue FCRL3 Tissue enhanced blood,lymphoid tissue FCRLB Tissue enhanced blood,lymphoid tissue FGF23 Tissue enhanced blood,heart muscle,liver,urinary bladder FOLR3 Tissue enhanced blood,bone marrow,lymphoid tissue GFRAL Tissue enhanced adipose tissue,blood,liver HBQ1 Group enriched blood,bone marrow,brain,vagina IFNL1 Tissue enhanced blood IL27 Tissue enhanced blood,brain,liver KIR3DL1 Tissue enhanced blood,bone marrow,lymphoid tissue LAIR2 Group enriched blood,lymphoid tissue,pituitary gland,placenta MAP1LC3B2 Tissue enriched blood MCEMP1 Group enriched blood,bone marrow,lung,lymphoid tissue OTOA Group enriched blood,lymphoid tissue,testis PADI4 Group enriched blood,bone marrow,lymphoid tissue PGLYRP1 Tissue enhanced blood,bone marrow RNASE3 Group enriched blood,bone marrow SCT Group enriched blood,intestine SHD Group enriched blood,brain,heart muscle SIGLEC5 Group enriched blood,bone marrow SIT1 Group enriched blood,lymphoid tissue TCL1A Group enriched blood,lymphoid tissue TEX101 Group enriched blood,testis TNFRSF13B Tissue enhanced blood,lymphoid tissue,skeletal muscle TNFRSF13C Tissue enhanced blood,intestine,lymphoid tissue TNFRSF9 Group enriched blood,lymphoid tissue TNFSF11 Tissue enhanced blood,lymphoid tissue TREML2 Tissue enhanced blood,bone marrow,placenta VSTM1 Group enriched blood,bone marrow Figure S14. Proportion of variance explained by various a priori defined sensitivity factors. Each box plot presents the median, first and third quartiles, with upper and lower whiskers representing 1.5x inter-quartile range above and below the third and first quartiles respectively; n=2,940 independent protein analytes tested.

Co-localization with expression QTLs
Colocalized pairs of pQTL and eQTL signals were largely concordant with respect to their effects on circulating proteins and tissue gene expression levels, with 82% (7,728/9,385) sharing the same direction of effect overall. We observed the lowest directional concordance rates in tissues from the brain (Extended Data Figure 10a), including the cerebellar hemisphere (64%; 85/133) and cerebellum (68%; 100/148), and the highest in the liver (90%; 146/141), which could potentially be explained by factors affecting access to circulation such as the blood brain barrier. One example of a gene with tissue-specific discordant direction was ADAM23, where 5 cis pQTL signals colocalized with ADAM23 eQTLs across 12 different tissues, showing discordant directions of effect for the cerebellum (Extended Data Figure   10b).

ABO blood group and FUT2 secretor epistasis effects provide insights into GI pathophysiology
We observed pleiotropic associations at the ABO blood group and fucosyltransferase 2 (FUT2) loci on chromosomes 9 and 19 respectively. The FUT2 enzyme facilitates expression of ABH antigens on red cells of corresponding blood groups in mucal and gastro-intestinal (GI) secretions. Approximately 20% of white Europeans are homozygous for deletion of the FUT2 functional secretor allele (rs601338, Trp154Ter), leading to truncation and inactivation of the enzyme and non-secretion of the blood group antigens 33 . The FUT2 deletion has been associated with cholestatic and gastrointestinal conditions [34][35][36] . This led us to explore the biologically informed hypothesis that FUT2 secretor status modifies the effect of blood group antigen expression on protein levels, serving as an example of long-range gene-by-gene interaction.
We found the proportion of FUT2 non-secretor carriers (25.4%) in UKB-PPP to be in line with reported prevalence and we did not observe any evidence of dependencies between ABO blood group genotypes and FUT2 secretor status (c 2 p=0.58). At a multiple testing corrected threshold of p<1.7x10 -5 , 432 proteins were associated with ABO blood groups and 225 proteins were associated with secretor status (Supplementary Table 27). We found significant interactions between blood group and secretor status for 55 proteins.
Some of the top interactions include MUC2, FAM3D and ALPI (Figure 4a). MUC2 levels were lower in secretors compared to non-secretors and similar between blood groups in nonsecretors only, whilst in secretors, we found O blood group had lower MUC2 levels vs A and in turn AB/B blood groups (Figure 4a left). For FAM3D, we found higher levels in secretors, with differences between blood groups only seen in secretors where the O blood group had lower levels vs non-O, and levels were higher in A than AB and in turn B blood groups (Figure   4a middle). For ALPI (intestinal alkaline phosphatase), we replicated and extend on the only previous reported such interaction effect seen in a smaller targeted study for alkaline phosphatase (ALP) in a Japanese cohort 37 -we found higher levels in secretors and differences between blood groups in secretors only where O blood groups had higher levels vs non-O, and levels were lower in A than AB and in turn B blood groups (Figure 4a right) We found significant gene expression enrichments for proteins with significant interaction effects across multiple human gastrointestinal tissues 45 , including duodenum, small intestine, colon, rectum, and pancreas -reinforcing with the role of FUT2 in GI secretions (Figure 4b   left). Enrichment in the intestine was also observed in orthologous genes in a mouse tissue expression data 46 (Figure 4b right), indicating a degree of conservation between these two species.
Our results provide evidence of blood group and secretor interaction in the modulation of proteomic concentrations, which may underlie potential susceptibility mechanisms to various FUT2/ABO associated GI conditions (such as IBD) and infections.

Proteomic insights into pathways underlying COVID-19 associated loci
The COVID-19 pandemic continues to accelerate research into the mechanisms and pathways influencing risk of COVID-19 infections and potential target candidates for drug compounds.
Here we integrated pQTL data with the largest GWAS meta-analysis of reported and hospitalized COVID-19 cases conducted to date (https://www.covid19hg.org/results/r7/) using multi-trait colocalization under the HyPrColoc framework 47 .
For six of the COVID-19 hospitalization loci, we found high posterior probability of colocalization (PP>0.9) with pQTLs for two proteins with expression in the lungs: surfactant protein D (SFTPD) and lysosome-associated membrane glycoprotein 3 (LAMP3)

Insights into inflammasome pathways
Inflammasomes are multimeric protein complexes that mediate innate immune responses, primarily through the activation of CASP1 and subsequent cleavage, activation, and noncanonical secretion of pro-inflammatory cytokines IL-18 and IL-1b 56,57 . Rare, protein altering variants in inflammasome components are known to cause many inherited autoinflammatory conditions 58 . The causal relationship between genetic alterations in the inflammasome and autoinflammation has been clinically validated by their successful treatment with anti-IL-1b therapies 59 .
Previous proteo-genomics studies using the SomaScan aptamer platform reported inconsistent trans pQTL results at the NLRP12 locus. The INTERVAL 11 study identified this as a pleotropic hot spot, whereas AGES and a cross-platform study did not 23,60 , leading to speculation that pQTL associations at this locus may be platform dependent and/or the result of inter-study differences in sample handling. The INTERVAL study identified more than 300 trans pQTLs within a 5Kb window of the NLRP12 locus. Using similar filtering criteria, we identified 44 trans pQTLs within the same window. While head-to-head comparisons between these studies are difficult due to greater statistical power of PPP and broader protein screening of INTERVAL, these results suggest the UKB-PPP findings are more consistent with the typical pleotropic loci identified in other studies.
In contrast to other loci such as ARHGEF3 (see Effects of blood cell counts section), trans pQTL results at the NLRP12 locus were robust to cell blood count sensitivity analyses. This supports the hypotheses that these signals result from inflammasome driven biological effects, rather than artificial increases in protein concentrations caused by sample handling related cell lysis.
Furthermore, we found a missense variant in NLRC5 (rs74439742, Pro191Gln) to be the lead trans pQTL associated with reduced levels of all MHC class I proteins tested: HLA-A, HLA-E and B2M, as well as BTN3A2. This is highly consistent with NLRC5 being a key transcriptional activator of MHC class I genes 61  Taken together, these results indicate that -in addition to known, rare, highly penetrant, disease-causing variants -common forms of genetic variability play a more subtle, but significant, role in inflammasome-mediated innate immune responses.

PCSK9 MR with extended instruments expands mirroring of clinical trial pharmacological effects on cholesterol and indicated diseases
The causal effects of PCSK9 levels on LDL and total cholesterol have been well established through various orthogonal means, with several randomized clinical trials demonstrating the efficacy of PCSK9 inhibitors on cholesterol levels and cardiovascular events [66][67][68][69] . Leveraging multiple cis pQTLs as genetic instruments to proxy directly for the effect of PCSK9 levels, we employed Mendelian randomization to examine causal effects of PCSK9 levels on lipids (HDL, LDL and total cholesterol), cardiovascular outcomes (coronary heart disease (CHD), myocardial infarction (MI)) and ischaemic stroke (IS: large-artery (IS-LA) and small-vessel We also compared the gain in precision for the PCSK9 genetic instrument and MR analyses between the current study and three previous pQTL studies [73][74][75]  on the lipid, CHD and stroke outcomes also became more significant in the current study when compared with instruments from previous pQTL studies, reflecting higher precision in the genetic instrument with increasing sample sizes (Extended Data Figure 12b). We note there is potential for a small degree of sample overlap between UKB instruments for PCSK9 and outcomes through random sampling of overlapping populations, but we expect this impact to be minimal.