Venous blood from the cubital vein of healthy volunteers was drawn in 10 ml EDTA Monoject tubes (Medtronic, Dublin) after obtaining written informed consent. PBMC isolation was performed as previously described (9). In short, the PBMC fraction was obtained using density centrifugation of EDTA blood diluted 1:1 in pyrogen-free saline over Ficoll-Paque (Pharmacia Biotech, Uppsala). Cells were then washed twice in saline and suspended in RPMI 1640 medium supplemented with gentamycin (10 mg/mL), L-glutamine (10 mM) and pyruvate (10 mM). PBMCs were counted in a Coulter counter (Beckman Coulter, Pasadena) and re-suspended in a concentration of 5 x 10⁶ cells/mL. A total of 5 x 10⁵ PBMCs were added in 100 μL to round-bottom 96-well plates (Greiner) and incubated with 100 μL of stimulus (heat-killed Candida albicans yeast of strain ATCC MYA-3573, UC 820, 1 x 10⁶/mL) or RPMI 1640 medium. After 24 hours, the supernatants were collected and stored at −20°C until assayed using the proximity extension assay (PEA) that is being used from Olink Proteomics (https://www.olink.com/).

C. albicans ATCC MYA-3573 (UC 820) was grown overnight to generate yeast cells in Sabouraud dextrose broth at 29°C, with shaking. Cells were harvested by centrifugation, washed twice with PBS, and resuspended in culture medium (RPMI 1640 Dutch modification). C. albicans yeast was heat killed for 30 min at 95°C.

All bulk RNA-seq data from PBMCs was previously generated (10) in 70 individuals from the GoNL cohort (15) This data was generated from PBMCs that were stimulated for 24 h with Candida or remained unstimulated (RPMI 1640 control condition). Principal components analysis (PCA) of the gene expression data generated from the Candida- and RPMI 1640-stimulated PBMCs showed that two samples were possibly mislabelled, which were subsequently removed before differential expression (DE) analysis ( Supplementary Figure 4 ). The differentially expressed genes upon stimulation were identified through DESeq2 package from Bioconductor using 72 Candida-stimulated- and 75 RPMI 1640-stimulated samples (19). The statistically significant threshold (FDR < 0.05) was applied.

Supernatant samples of PBMCs stimulated with C. albicans yeast and with RPMI 1640 control medium were analyzed using the proximity extension assay (PEA) that is being used from Olink Proteomics (20). We selected the Olink^® inflammatory panel (Olink Bioscience AB, Uppsala, Sweden). This panel includes 92 inflammatory proteins. To note, due to technical issues identified in the brain-derived neurotrophic factor (BDNF) assay from the Olink proteomics, we were unable to report any data for this analyte, and thus, we reported in this study only 91 inflammatory proteins ( Table S1 ). Olink proteomics delivers the protein measurements as Normalized Protein Expression (NPX) values, which are Cq values normalized by the subtraction of values for extension control, as well as inter-plate control. The scale is shifted using a correction factor (normal background noise) and reported in the Log₂ scale. We performed normality test on both raw and log-transformed data using Shapiro-Wilk normality test, and a P > 0.05 was used as a threshold for normal distribution. The majority of proteins followed a non-Gaussian distribution. The distribution of all protein measurements generated from Candida-stimulated PBMCs from 500FG and Lifelines Deep cohort are shown in Supplementary Figures 5A, B respectively. In addition to QC for protein distribution, Olink Proteomics performed QC per sample during which samples that deviate less than 0.3 NPX from the median pass the quality control.

To assess the correlation structure between pairs of inflammatory proteins measured in at least 85% of samples from Candida-stimulated PBMCs isolated from the 500FG cohort, unsupervised clustering was performed using Spearman’s correlation as the measure of similarity. All statistical analyses were performed in R version 3.2.2.

To test which proteins are differentially expressed in Candida-stimulated PBMC samples (n = 364) versus RPMI 1640-stimulated PBMC samples (n = 24) in the 500FG cohort, we performed a DE analysis using a linear model with age and sex as covariates. Twenty-six out of 364 Candida-stimulated PBMC samples did not pass QC per sample performed by OLINK proteomics as described above ( Figure 1 ). All samples stimulated with RPMI 1640 control medium passed quality control. After removing outliers, 337 Candida-stimulated samples (one sample out of the 338 samples was removed due to missing phenotype information) and 24 RPMI 1640-stimulated samples used for DE analysis. The total number of proteins that were measured in at least 85% of all samples and tested for DE analysis was 36 proteins. For DE analysis, we used the R package Limma adjusted to protein data, which is originally being used for the analysis of gene expression data (21). Limma uses an empirical Bayes method to moderate the standard errors of the estimated log-fold changes. Proteins with an FDR < 0.05 and a logarithm of fold change (logFC) > 2 compared to RPMI 1640-stimulated samples were considered statistically significant. Normalized protein measurements (NPX values) for all samples used in the DE analysis are provided in Table S2 .

We performed a meta-analysis of Candida-induced pQTLs by using the summary statistics of the pQTL datasets that were mapped in the two independent population-based cohorts, the 500FG and Lifelines Deep cohort. Lack of protein measurements for all available individuals restricted us to select 364 individuals from 500FG and 63 samples from Lifelines Deep cohort for whom both genotype and protein data were available ( Figure 1 and Supplementary Figure 1 ). In addition to the outliers that did not pass the QC from Olink Proteomics, we removed genetic outliers before QTL mapping as described above. In total, we removed 14 genetic outliers from the 500FG cohort and 11 genetic outliers from the Lifelines Deep cohort, resulting in 326 samples from 500FG cohort and 52 samples from the Lifelines Deep cohort ( Supplementary Figure 1 ). To correct the protein distributions for QTL mapping, we used a linear model adjusting for age and sex. NPX values of protein levels were used, and the correlation between protein abundance and genotype was tested independently for the two cohorts using the R-package Matrix-eQTL (22).

Next, we performed meta-analysis using the summary statistics of pQTLs that influence proteins that were measured in at least 85% of the samples in both cohorts. HGF protein was one of the proteins that was not measured sufficiently in the 500FG cohort in contrast to the Lifelines Deep cohort ( Table S1 ), and, hence, the number of proteins that used for meta-analysis was 35 proteins. To jointly test the effect of genetic variation in protein levels measured in both 500FG and Lifelines cohort, we used METAL (http://www.sph.umich.edu/csg/abecasis/metal/) using 35 proteins from 326 samples of the 500FG cohort and 52 samples of the Lifelines Deep cohort. Finally, we removed pQTLs from meta-analysis with opposite direction and those that show significant heterogeneity (P < 0.05) between the two cohorts (23). We considered P < 5 x 10^-8 to be the threshold for genome-wide significant pQTLs.

Both genotype and gene expression data upon Candida stimulation were available from 62 individuals of the GoNL cohort (15). The stimulation experiments and the RNA sequencing analysis of this dataset has been previously described (10). There were 32, 22 and 8 individuals carrying AA, GA and GG genotypes respectively at rs10114707 locus. We combined GA and GG in one group and performed DE analysis between two genotype groups (GA+GG versus AA) using the expression data of cis-genes (n = 30) extracted from the 5 GWAS loci (in a window of 1MB around each GWAS locus), which were identified in the meta-analysis of protein-QTLs. For DE analysis using gene counts, we used the DESeq2 package (19). Summary statistics of the DE analysis is shown in Table S3 .

To test all reported genome-wide significant pQTLs for pleiotropy, we initially inspected all genome-wide significant trans-pQTLs whether they show association with multiple distinct proteins. To correct for the multiple testing burden for all genome-wide pQTLs, we next inspected all genome-wide significant trans-pQTLs that had at least two associations with distinct proteins at P < 0.05/(5*35) = 0.0003. This cut-off represents a conservative approach to the multiple testing burden for all identified, independent GWAS SNPs (n = 5) with all protein traits (n = 35). The resulting association matrix was then visualized as a heatmap based on the negative log2 of the P values of association.

To identify genetic variants that are associated with survival, we ran a within-case genome-wide association study using our survival cohort consisted of 171 candidaemia patients, of which 115 were identified as non-survivors and 56 as survivors. After quality control per SNP and sample, we tested for association with survival of ~5 million SNPs using 104 non-survivors and 52 survivors. Age and sex were included as covariates in an additive model using a frequentist association test from SNPTEST v.2.5.4-beta3 (24). The quantile-quantile (QQ) plot that shows the distribution of the observed negative log₁₀ P values for the whole-genome SNPs against the theoretical distribution of expected negative log₁₀ P values indicated no or little evidence of population stratification ( Supplementary Figure 6 ). The genomic inflation factor based on median was equal to λ = 1.03.

To systematically study the inflammatory responses in the context of C. albicans infection, we first obtained PBMCs from 364 individuals of European origin from the 500FG cohort. We next measured the abundance of 91 inflammatory proteins in the supernatant of PBMCs in response to C. albicans yeast stimulation ( Figure 1 ). Of these proteins, 36 were detected in at least 85% of the samples ( Table S1 ), and all proteins showed significant inter-individual differences to C. albicans stimulation compared to RPMI 1640 medium stimulation ( Supplementary Figure 7 ). After removing samples that failed the QC (n = 24) and a sample with no age and sex information, we compared the abundance of the 36 proteins between stimulated (n = 337) and un-stimulated (using RPMI 1640 medium as control, n = 24) conditions by performing DE protein analysis using a linear model with age and sex as covariates. We identified significant up-regulation of 25 proteins out of 36 measured in at least 85% of PBMC samples that showed logFC > 2 compared to RPMI 1640-stimulated samples (adjusted P value < 0.05) ( Supplementary Figure 8A and Table S4 ). Next, we compared the correlation structure among the inflammatory proteins released from PBMCs by performing unsupervised clustering of the protein responses. Clustering revealed that the majority of proteins showed a strong positive correlation (P < 0.01) ( Supplementary Figure 8B , Table S5 ), with the exception of IL-8 and MCP-1, indicating the coordinated inflammatory protein response upon Candida stimulation.

Next, we investigated whether host genetic variation affects the inter-individual differences in inflammatory responses to C. albicans stimulation. For this, we used the genome-wide SNP genotype data and protein measurements of Candida-stimulated PBMCs isolated from two independent population-based cohorts consisting of 364 individuals from 500FG cohort and 63 individuals from Lifelines Deep cohort. Upon quality control (see Materials and Methods), we obtained 326 samples from 500FG cohort and 52 samples from the Lifelines Deep cohort with a total of 35 inflammatory proteins that were measured in at least 85% of the samples in both cohorts. For pQTL mapping, we selected SNPs that showed a minor allele frequency (MAF) >= 5% and passed other quality filters (see Materials and Methods). Using the protein and genotype data, we mapped pQTLs using each cohort independently and performed meta-analysis using the summary statistics of pQTLs identified in each cohort (n = 378). Meta-analysis revealed five independent trans-pQTLs that reached genome-wide significance level (P < 5 x 10^-8, P for heterogeneity > 0.05) ( Table 1 and Figure 2 ). These include one independent pQTL for IL-8 (chromosome 9), one for MCP-2 (chromosome 4), one for MMP-1 (chromosome 9) and two for CCL3 (chromosome 9 and 15).

In these five independent loci, we prioritized potential causal genes following two approaches. First, we tested whether genes that are located within 1MB window around the genome-wide significant SNPs were differentially expressed in response to C. abicans yeast stimulation after 24 hours in PBMCs using our previously published Candida-induced transcriptomic dataset generated from the Dutch population-based cohort, GoNL ( Table S6 ) (10). Secondly, we made use of the largest cis- and trans-expression quantitative trait locus (eQTL) study in blood from a total of 31,684 individuals through the eQTLGen Consortium (25) to test if the genome-wide SNPs affect gene expression in whole blood ( Table 1 ) or in human monocytes (26). Our strategy prioritized several interesting genes including CD38 and GADD45G. Briefly, CD38 gene is prioritized from the rs591999268 locus at chromosome 4 based on its significant differential expression upon 24-hour Candida stimulation (adjusted P 9.75 x 10^-50, log2 of fold change = 3.25). This finding suggests that CD38 may affect protein levels upon C. albicans stimulation. CD38 is a glycoprotein widely expressed in cells from the immune system, and animal studies indicate that it confers protection against infection by several bacterial and parasitic pathogens (27), suggesting its role on regulation of inflammatory protein response. As a second example, growth arrest and DNA damage inducible gamma (GADD45G) gene from the rs10114707 locus at chromosome 9 was prioritized based on two levels of evidence. Firstly, GADD45G is significantly differentially expressed in response to C. albicans at 24 hours (adjusted P 4.61 x 10^-97, log2 of fold change = 2.63). Secondly, we tested if GADD45G gene is differentially expressed among individuals carrying different genotypes (AA, GA and GG) at rs10114707, which could explain how SNP rs1011477 can affect the levels of inflammatory proteins. For this, we stratified individuals from the GoNL cohort based on their genotype at rs10114707 locus and tested whether cis-genes in a window of 1MB around each GWAS locus is differentially expressed. We observed that individuals heterozygous or homozygous for the minor allele G at rs10114707 tend to show a decrease in GADD45G expression compared to individuals homozygous for allele A ( Figure 3 and Table S3 ). Such stimulation specific expression QTLs will be helpful to identify potential regulators of inflammatory protein response.

It is possible that the same genetic locus may regulate multiple proteins (pleiotropy). To test the presence of such pleiotropic effects, we extracted all possible associations of genome-wide significant pQTLs with the 35 proteins in the same allelic direction between the two cohorts. We observed that all genome-wide significant pQTLs influence multiple distinct proteins ( Figure 4A and Table S7 ), indicating the pleiotropic effects of the pQTLs on the proteins. To account for multiple testing, we inspected all genome-wide significant trans-pQTLs that had at least two associations with distinct proteins at P = 0.05/(5*35) = 0.0003. This cut-off represents a conservative approach to correct for multiple testing for all identified GWAS SNPs (n = 5) with all protein traits (n = 35). SNP rs59199268 at chromosome 4 and rs55642284 at chromosome 15 showed only a single association with MCP2 and CCL3, respectively. We observed three loci on chromosome 9 at genetic variants rs10114707, rs12343816 and rs658817 that influence four, two and another two proteins respectively at P < 0.0003 in the same allelic direction ( Table S7 ). These pleiotropic effects suggest that these proteins are co-regulated in inflammation induced by C. albicans. For instance, individuals carrying the minor allele G at rs10114707 seem to upregulate the levels of VEGFA, TNF, MMP-1, and CCL20 ( Figures 4B–E ). We prioritized GADD45G as a potential causal gene at this locus as discussed above. GADD45 proteins are known to be induced by a number of cytokines and by bacterial lipopolysaccharide, have an essential role in the differentiation of myeloid cells as well as in the function of granulocytes and macrophages (28, 29). Given this important immune-regulatory function it is plausible that GADD45G could be one of the core regulators of multiple inflammatory proteins. However, more studies are needed to explain the specific mechanism by which these pleiotropic effects occur.

Next, we aimed to investigate whether genetic variation that influences inflammatory proteins (pQTLs) in blood circulation may determine patient survival to candidaemia. To identify genetic variants that are associated with survival in candidaemia patients, we performed a within-case GWAS analysis using 156 individuals, of which 104 were identified as non-survivors and 52 as survivors. After QC per SNP and sample, we tested 5450489 variants (MAF > 1%) for association with survival using an additive model with age and gender as covariates. No genome-wide significant associations were observed (P < 5 x 10^-8), which can be explained due to the small number of patients with reported survival outcome ( Figure 5 ). To test whether genetic variants that are associated with survival influence inflammatory proteins in response to Candida, we extracted the within-case GWAS summary statistics of survival and overlaid this with our trans-pQTLs at a P threshold 1 x 10^-4. A trans-pQTL, rs7651677, at chromosome 3 that influences uPA (P 8.29 x 10^-6, z score 4.46), also known as urokinase, showed a significant association with survival (P 3.52 x 10^-5 and OR 3), indicating that this locus may determine survival of candidaemia patients through modulating uPA. In this regard, it is interesting to note that the serum suPAR (soluble urokinase plasminogen activator receptor) levels are considered as a biomarker to assess the risk in sepsis (30).