The 500FG cohort of healthy individuals of Western European ancestry comprises of 237 males and 296 females with age range of 18 to 75 years, being part of the Human Functional Genomics Project, HFGP (www.humanfunctionalgenomics.org ).

The 1M-scBloodNL cohort comprises of 120 individuals, 53 males and 67 females between 27 to 78 years. This cohort is part of the Lifelines DEEP cohort, a prospective population cohort of participants from the northern Netherlands (11).

Inflammatory protein concentrations were measured using Olink^® proteomics platform. In fact, supernatants were collected after 24h of Candida albicans stimulation and submitted to Olink Proteomics for analysis using the inflammatory panel assay of 92 analytes. Olink data are presented as Normalized Protein eXpression values (NPX, based on log2 scale). Immunoassays utilized by Olink are based on the Proximity Extension Assay (PEA) technology (14), which makes use of oligonucleotide-labeled antibodies binding to their respective protein. When the two antibodies are brought in proximity, a DNA polymerase target sequence is formed, which is subsequently quantified by quantitative real-time polymerase chain reaction (qPCR).

Filtering of Olink generated data was restricted to only Proteins as the samples passed Olink internal quality control across all proteins. We excluded all proteins which failed to be quantified in at least 85% of the samples, meaning all proteins with more than 15% samples missingness (NPX value below the protein-specific limit of detection (LOD) value) were excluded from downstream analysis. The remaining Proteins with NPX values below the LOD were replaced with protein-specific LOD values.

Following per-protein cleaning, 42 and 35 proteins were available for the 1M-scBlodNL and 500FG cohorts respectively ( Figure S1A ). Out of these proteins, 26 were common between both cohorts.

The protein distributions on log2 scale were not normally distributed ( Figure S1B ). We applied rank-based inverse normal transformation as implemented in the GenABEL R package (15), to transform the data to mimic Gaussian distribution ( Figure S1C ).

The procedures for genotyping, genetic data filtering and genotype imputation of the 500FG cohort had been previously described (6). Extracted DNA was genotyped using the commercially available SNP chip, Illumina HumanOmniExpressExome-8 v1.0. Following pre-imputation filtering steps for both markers and individuals, the remaining dataset SNP genotypes were imputed with GoNL as reference panel (16).

For Lifelines Deep cohort, genotyping and imputation was performed as previously described (17). Both the HumanCytoSNP-12 BeadChip and the ImmunoChip platforms (Illumina, San Diego, CA, USA) were used to genotype the isolated DNA. Independent markers quality control was performed for both platforms and subsequently merged into one dataset. After merging, genotype SNPs were imputed using IMPUTE2 (18) against the GoNL reference panel.

Pairwise Pearson correlation analysis using the ‘corrr” R package was performed on the normalized protein abundances after adjusting for age and sex. Based on Pearson’s correlations for each pair of proteins, co-expression protein networks were reconstructed (using significant correlation coefficient threshold of 0.7, (absolute(r) > 0.70), a cut-off denoting a very strong strength of association.

The association analysis of genotype-phenotype correlation was carried out using two main approaches. In the first method, the univariate approach was performed using the linear regression in PLINK (19). The pQTL analysis was conducted independently for both cohorts, that is one analysis for the 1M-scBloodNL and another for 500FG cohorts. In the second method, multivariate test of association based on canonical correlation analysis (CCA) (20), was conducted. CCA extracts the linear combination of highly corrected traits that explain the largest possible amount of the covariation between genetic variants and all traits (Proteins in this case). To control for potential confounding factors, we adjusted for covariates such as age and sex on normalized protein abundances and, regressed the residuals of each protein and protein network on SNP genotypes.

Summary statistics from both primary studies were utilized to perform meta-analyses.

Association results for both the 500FG and 1M-scBloodNL cohorts obtained from the univariate approach were combined using the weighted sum fixed-effect model as implemented in the METAL software program (21). The multivariate approach implemented in PLINK does not compute beta estimates and standard errors. Therefore, the meta-analysis of the multivariate P-values was carried out using sum of z (Stouffer’s) method as implemented in the metap R package (22).

Colocalization analysis of cell-type-specific cis-eQTL and pQTL signals was conducted using Bayesian colocalization method which is implemented in the coloc package in R (23). We retrieved the genome-wide cell-type-specific cis-eQTL summary statistics from a previously published study using the 1M-scBloodNL cohort (13). Additionally, we performed trans-eQTL mapping only for the top pQTLs to ascertain their trans-eQTL effect. The eQTLs from this study were identified using PBMCs in an unstimulated condition as well as after 3h and 24h in vitro-stimulations with three different pathogen stimulations, namely Candida albicans, M. tuberculosis and P. aeruginosa. Cis-eQTL was defined as SNP-gene distance of 100kb window and FDR < 0.05. However, for overlap comparison with pQTLs, we re-tested the relevant SNPs to ensure that SNPs outside the defined cis-region are included.

From our pQTL mapping results, we selected the index variants (SNPs with the smallest P value) for each protein and extracted all variants within a window size of 1Mb around the index pQTL variant for further analysis. The default prior probability (1x10^-6) that a random variant in the region is causal to both traits was applied. A posterior probability (PP4 >= 0.75) is considered as strong evidence of colocalization. LocusCompareR, being an R package was used for the visualization of results (24).

To identify SNPs associated with the concentration of specific proteins induced upon 24h Candida albicans stimulations, pQTL analysis was performed in two independent population-based cohorts. A total of 445 participants (500FG (N= 325) and 1M-scBloodNL (N=120) cohort) were studied for whom both genotype and Olink protein abundances were measured upon ex vivo 24h Candida albicans stimulation of their PBMCs ( Figure 1 ). After genotype imputation and quality control, 4,095,761 SNPs and 26 inflammatory proteins remained that were common between both cohorts and that were used as input for the pQTL analysis. The protein data for both cohorts were generated in different batches. Therefore, pQTL analysis was run in each cohort separately, after which results ( Figure S2 ) were integrated using meta-analysis to increase statistical power. We identified a genome-wide significant cis-acting pQTL variant rs484915 (P value = 1.81x10^-8) on chromosome 11 correlating with MMP-1 production ( Figure 2A ; Table 1 ). In addition, the other top SNPs correlating with the remaining 25 proteins were trans-acting pQTLs with strong suggestive associations (P value > 5.0 x10^-8 to 5 x10^-6). For example, the second most significant hit aside the genome-wide significant cQTL, is an intronic SNP rs62205465 residing in the ZNF133 locus and exhibited an association strength closed to the genome-wide significant threshold (P value = 5.80 x10^-8) with MCP-3 concentrations. Figure 2A illustrates the association results of all the proteins analyzed and their corresponding top SNPs (lowest P-value).

Previously, we had conducted a genome-wide eQTL analysis per major cell type using scRNA-seq data from unstimulated and 3h and 24h pathogen-stimulated (Candida albicans, M. tuberculosis, P. aeruginosa) PBMCs in the 120 individuals from the 1M-scBloodNL cohort (13). This data was used to interrogate whether the same pQTL SNP could also affect gene expression levels to explore the degree of association between pQTLs and eQTLs. Note that the pQTL data captures the bulk secretion of proteins coming from PBMCs that were stimulated for 24h with Candida albicans, whereas the eQTL data was collected for each cell type separately. Out of the 24 unique top pQTL variants with a suggestive or genome-wide significant association (P value from 1.81 x 10^-8 to 5.89 x 10^-6) identified from meta-analyzed results of the univariate approach ( Table 1 ), we could only overlay 20 with the scRNA-seq derived eQTL data of the 1M-scBloodNL cohort, as eQTL data for the remaining 4 SNPs (rs36067904, rs13033376, rs62129298, and rs7018706) were not available. Seven out of the 20 tested SNPs showed an eQTL effect upon 24h Candida albicans stimulation in at least one cell type (FDR < 0.05) ( Figure S3 ). Among these seven, two pQTLs showed genome-wide significant association with gene expression levels specifically in monocytes. The first cis-acting pQTL, rs484915 showed association with both MMP-1 protein concentrations and gene expression ( Figure 2B ). Of note, only this pQTL among the seven was associated with the corresponding gene of the protein. The second is trans-acting SNP rs10276582 associated with CCL4 protein production affected AMZ1 gene expression levels ( Figure 2C ). This observation suggests AMZ1 may regulate CCL4 protein levels depending on an individual’s genotype at SNP rs10276582 (or any other SNP in high LD).

To elucidate whether the effects of the significant eQTLs were detectable only in Candida stimulations after 24h, we further assessed their impact in unstimulated condition, 3h Candida stimulations, as well as in Mycobacterium tuberculosis stimulations after 24h. We observed that the effects of the cis-eQTLs were less pronounced before stimulation and after 3h Candida stimulation, which suggests temporal regulation of gene expression following Candida stimulation. Nearly equal strength of association was identified in the case of M. tuberculosis stimulations after 24h for only the trans-acting variant SNP rs10276582 ( Figure 2C ).

Next, to uncover the potential mechanisms underlying the observed pQTLs, we tested whether SNPs impacting protein concentrations at specific loci are also the same causal regulatory eQTLs through colocalization analysis. We identified strong evidence of colocalization for both pQTLs and eQTLs. For instance, in the MMP-1 locus ( Figure 3A ), the posterior probability (PP.H4) was 0.998 ( Figure 3B ). Also, the trans-acting pQTL SNP rs10276582 located in the AMZ1 locus on chromosome 7 ( Figure 3C ), showed significant colocalization (PP.H4 = 0.996) with the eQTL effect. eQTL analysis revealed association of this variant with three different cis genes, GNA12 (P = 7.1 x 10^-1), TTYH3 (P = 1.7 x 10^-2) and AMZ1 (P = 7.31 x 10^-12) ( Figure 3D ), indicating the presence of multiple causal genes at this locus. However, both GNA12 and TTYH3 exhibited weaker strength of association than AMZ1, which showed statistically significant association. Thus, pinpointing AMZ1- which is predicted to belong to a large metallopeptidase protein family and interacts with cell receptors and growth factors, as the potential causal gene in this genomic region. This observation demonstrates that cis-regulation of gene expression levels maybe involved in the mechanisms by which distal variants impact protein expression.

It is possible that some of the genetic variants may affect multiple proteins, so we wanted to explore whether there is strong correlation between protein concentrations upon stimulation. To explore the patterns underlying Candida-induced protein production, we performed correlation analyses. We observed mostly significant positive pairwise correlation (with the exception of MCP-1 and MCP-2) between the 26 proteins in the 500FG cohort, which contains the largest number of individuals from which samples were collected ( Figure 4A ). However, in the 1M-scBloodNL cohort, divergent patterns of correlation strengths were observed, including weaker and negative correlations between CSF-1, IL-12B and IFN-gamma, contrasting the positive correlation observed in the 500FG cohort ( Figure 4B ). To understand what might underlie this observation, we performed the Fligner-Killen test to evaluate whether there is significant donor variation of these proteins between the two cohorts. While we observe significant difference for CSF-1 (Fligner-Killeen:med chi-squared = 19.893, p-value = 8.19 x 10^-6) and IL-12B (Fligner-Killeen:med chi-squared = 5.0642, p-value = 0.02442), there was no evidence to suggest that the variance in IFN-gamma concentrations significantly differ between both cohorts (Fligner-Killeen:med chi-squared = 0.3097, p-value = 0.5779).

In the 500FG cohort, the strongly pairwise correlated proteins acting as protein network consisted of 24 out of the 26 proteins common between both cohorts ( Figure 4C ). For example, the correlation coefficient between uPA and OSM was as high as 0.91, while CD40 exhibited an approximately 0.9 between CD5 and CD6. Next, we sought to more accurately infer the protein network by replicating the analysis in the 1M-scBloodNL cohort. While MCP-1 exhibited strong negative correlation with some proteins (such as OSM and uPA) in the 500FG cohort, this pattern of correlation was hidden in the network obtained from the 1M-scBloodNL cohort ( Figure 4D ). A rather weaker and positive correlation was manifested between MCP-1 and OSM (0.26) and uPA (0.24) and thus, MCP-1 was excluded from multivariate genetic association analysis as shared genetic architecture is likely to underlie or trigger the observed strong positive correlations.

Taking advantage of the tightly co-regulated pattern of the proteins, we sought to identify SNPs impacting protein network upon Candida stimulations aside the univariate association analysis. After independent analyses in both cohorts, we identified a genome-wide significant hit in the 500FG cohort as well as several suggestive associations. The genome-wide significant SNP was rs938662 (P value = 4.37 x 10^-8), residing in the PRKCE locus ( Figure S6A ). This suggest that PRKCE locus exhibit a pleotropic role as the top SNP in this genomic region is associated with multiple highly correlated proteins. There was no evidence of inflation of the test statistics as genomic inflation factor (Λ) was computed to be 1.01 ( Figure S6B ). In the case of the 1M-scBloodNL dataset, no statistically significant SNP was found to be associated with proteins ( Figure S6C ), and corresponding genomic inflation factor (Λ) was 0.99, indicating lack of inflation of test statistics ( Figure S6D ).

To identify true or robust association signals, we synthesized summary statistics from both cohorts through meta-analysis. Even though no loci reached genome-wide significance, we identified strong suggestive associations ( Figure 5A ). The top signal identified was rs484915 (P value = 1.05x10^-7), located at the MMP1 locus and has previously been shown to alter expression levels of cis-genes in multiple tissues and also blood protein concentrations (26). In the univariate analysis, we found the same SNP rs484915 to be associated with MMP-1 with slightly much stronger association based on P-value (1.81x10^-8). Interestingly, five proteins, namely MMP-1, CXCL5, CCL20, CXCL1 and CXCL10 among the proteins forming the network ( Figure 5B ), contributed with relatively stronger weights or correlation coefficients to the observed association result of SNP rs484915 to the protein network. This observation from the multivariate approach demonstrates the pleiotropic effect of SNP rs484915 which cannot be captured directly via univariate analysis and also tease apart the main proteins whose expression levels are being regulated.

Next, we sought to investigate whether the joint analysis of multiple correlated proteins with genetic variants offers some advantage over univariate analysis. To achieve this aim, we directly compared the P-values or distribution of P-values of genetic variations identified using both approaches and restricted the analysis to only the largest cohort, 500FG. We reasoned that such analysis will be more credible to be conducted in a specific cohort due to the variation seen in proteins’ correlation structure between both cohorts. As expected, we observed stronger associations emanating from the multivariate analysis compared to the univariate manner ( Figure 5C ). The top 6 independent strong suggestive loci (P < 9 x10^-7) identified via the multivariate approach were in all cases showing a much stronger association when compared to the strength of association (P values) of each protein analyzed separately ( Figure 5D ). Furthermore, to characterize the identified pQTLs, we targeted the top 6 strongest pQTL variants (rs938662, rs1501565, rs188465730, rs1020993, and rs11902595), given the lack statistically significant SNPs to elucidate their effect on cell-type specific eQTLs. We found 3 out of 6 to be an eQTL as well in at least one cell type albeit one with nominally significant effect. The overlapped SNPs showed weak association with various cis-genes. For instance, the strongest effect was observed for intronic SNP rs11902595 on chromosome 2 which was nominally (P= 0.013) correlated with a lincRNA (RP11-191L7.1) in CD8T cells.

We further evaluated the overlap and strength of association between our Candida 24h stimulated pQTL associations and previously reported highly powered (21,758 participants) unstimulated pQTL study published under the SCALLOP consortium (5). We found only 8 proteins (MMP-1, CSF-1, CXCL1, EN-RAGE, CCL4, MCP-1, CD40, CCL20) common between the 90 cardiovascular proteins measured in the previous study and the 26 proteins measured in the inflammatory panel used in our study. Using nominal significance P-value < 0.05, the percentage of shared Candida 24h stimulated pQTLs vs pQTLs identified using the SCALLOP consortium data, ranged from 4.5% to 5.6%. Of these, the top SNPs among the shared variants show nominal association with 6 proteins after stimulation. However, cis-acting genome-wide significant pQTL variant rs484915 (P = 1.81 x 10^-8) correlating with MMP-1 upon stimulation showed very strong association (P = 1.87 x 10^-220) in the SCALLOP consortium data ( Figure 6A ). We further observed one trans-acting pQTL variant rs3014874 exhibiting suggestive association with EN-RAGE protein (5.87 x 10^-4), but exceeded the genome-wide significant threshold with P = 1.64 x 10^-29 in the SCALLOP consortium data ( Figure 6B ). Evidence from previous studies indicate that this downstream variant (rs3014874) located on chromosome 1 affects the expression levels of cis-genes in blood as well as different tissues ( Figure 6C ). Given that eQTL effects are dependent on context and tissue being studied, we investigated the effect on this SNP on nearby genes in different cell types after candida stimulation. Among all the cis-genes, SNP rs3014874 showed association with only S100A9 gene ( Figure 6D ), making it the likely causal gene. This observation further highlights the role of trans-regulatory network in determining the abundance of circulating plasma proteins in blood.