An integrative genomics approach identifies KDM4 as a modulator of trained immunity.

Innate immune cells are able to build memory characteristics via a process termed "trained immunity." Host factors that influence the magnitude of the individual trained immunity response remain largely unknown. Using an integrative genomics approach, our study aimed to prioritize and understand the role of specific genes in trained immunity responses. In vitro-induced trained immunity responses were assessed in two independent population-based cohorts of healthy individuals, the 300 Bacillus Calmette-Guérin (300BCG; n = 267) and 200 Functional Genomics (200FG; n = 110) cohorts from the Human Functional Genomics Project. Genetic loci that influence cytokine responses upon trained immunity were identified by conducting a meta-analysis of QTLs identified in the 300BCG and 200FG cohorts. From the identified QTL loci, we functionally validated the role of PI3K-Akt signaling pathway and two genes that belong to the family of Siglec receptors (Siglec-5 and Siglec-14). Furthermore, we identified the H3K9 histone demethylases of the KDM4 family as major regulators of trained immunity responses. These data pinpoint an important role of metabolic and epigenetic processes in the regulation of trained immunity responses, and these findings may open new avenues for vaccine design and therapeutic interventions.

Introduction

Upon recognition of pathogen‐associated molecular patterns, invading microorganisms induce a first line of immune defense consisting of production of cytokines, defensins, and reactive oxygen radicals, as well as activation of phagocytosis and intracellular killing by phagocytic cells [1]. These innate immune responses are rapid and efficient, yet they were believed to be nonspecific and lacking immunological memory. However, the ability of the innate immune system to adapt after previous insults has been recently described, representing a de facto nonspecific innate immune memory, also termed “trained immunity” [2]. Certain infections or vaccinations induce long‐term changes in the response of innate immune cells, which lead to an enhanced inflammatory response when these cells are stimulated by similar or unrelated secondary stimuli [3]. Microbial molecular patterns such as β‐glucan, a component of the Candida albicans cell wall, and the tuberculosis vaccine Bacillus Calmette‐Guérin (BCG), have been shown to induce trained immunity both in vitro and in vivo [4, 5, 6]. A full understanding of the molecular mechanisms that underlie this process is still ongoing, but studies have shown that the enhanced immune response is at least partially mediated through metabolic and epigenetic reprogramming of innate immune cells and their progenitors [7, 8]. In particular, the induction of trained immunity has been linked to genome‐wide changes in the distribution of histone H3 acetylation at lysine 27 (H3K27ac) and histone 3 trimethylation of lysine 4 (H3K4me3) at the promoter sites of proinflammatory genes [4, 5, 6, 9]. Accordingly, enzymes involved in adding and removing these histone modifications may be important drug targets to modulate trained immunity responses. However, very little is known regarding the role of epigenetic enzymes that mediate these processes, except for KDM5 histone demethylases that have been shown to be inhibited during induction of trained immunity responses [10]. Therefore, there is a need to further understand how epigenetic enzymes regulate trained immunity.

Previous BCG vaccination studies indicated considerable interindividual variation in the ability of innate immune cells to build memory responses [6, 11], which may influence the nonspecific protective effects of BCG vaccination and susceptibility to trained immunity‐mediated diseases, such as atherosclerosis. Both genetic and nongenetic factors may influence the magnitude of the trained immunity responses; however, a systematic genome‐wide study on the role of host genetics, and how it regulates trained immunity, has not yet been performed. Quantitative trait loci (QTL) studies will allow us to characterize genetic loci that have an impact on trained immunity and molecular pathways that underlie this process. Conducting a QTL study is challenging as it requires a large number of samples, and it also provides limited insight into how genes are linked to trained immunity responses. Alternatively, integrative strategies in which QTL studies are combined with follow‐up functional validation studies can help to prioritize candidate genes and to better understand their role in these responses.

In the present study, we first assessed the impact of age and sex on in vitro trained immunity responses using two independent population‐based cohorts of healthy individuals with Western European ancestry, the 300BCG and 200 Functional Genomics (200FG) cohorts from the Human Functional Genomics Project (http://www.humanfunctionalgenomics.org). Furthermore, we assessed the impact of genetic variation on in vitro‐induced trained immunity by performing a meta‐analysis of independent QTLs identified in the 300BCG and 200FG cohorts. We identified several suggestive, independent genetic loci that influence trained immunity (P < 1 × 10–5). However, the significance of these loci cannot be ascertained solely genetically due to the possibility of false‐positive findings, and functional validation is needed to confirm their biological relevance. In order to assess the importance of certain genes from these QTL loci, we prioritized three candidate genes for functional validation as important mediators of trained immunity: one that codes for a known epigenetic regulator that belongs to the KDM4 family of histone demethylases, and two that belong to the family of sialic‐acid‐binding Ig‐like lectins (Siglecs) that code for the Siglec‐5 and Siglec‐14 receptor, which are key regulators of immune‐cell signaling. In addition, we found that QTLs that influence cytokine levels upon induction of trained immunity with β‐glucan are significantly enriched for genes involved in the PI3K‐Akt signaling and the ECM‐receptor interaction pathway. Together, these findings demonstrate that QTL studies combined with functional validation help prioritize candidate genes, such as Sigelc‐5, Siglec‐14, and KDM4, which have an important role in trained immunity, and therefore contribute to a better understanding of trained immunity. Subsequently, a better understanding of these factors and pathways will help design new immunotherapeutic strategies against diseases in which trained immunity plays a role.

Results

Variability in human trained immunity responses

To investigate the interindividual variability in trained immunity responses, we assessed the in vitro induction of trained immunity in adherent monocytes in two independent cohorts of 267 and 110 healthy individuals from the 300BCG and 200FG cohort, respectively. Monocytes were stimulated for 24 h with RPMI control medium, β‐glucan, or BCG to mimic a primary infection [12, 13]. To simulate a secondary infection, cells were re‐stimulated on day 6 with either LPS or RPMI control medium for 24 h (Fig. 1A) and pro‐inflammatory cytokine responses were measured. Similar to previous studies, stimulation with β‐glucan or BCG‐induced trained immunity, as characterized by an enhanced production of IL‐6 and TNF‐α upon re‐stimulation with LPS as compared to RPMI control (Fig. 1B; 300BCG cohort, and Supporting Information Fig. S1 and S2, 200FG cohort). A strong correlation was observed between the fold change in IL‐6 and TNF‐α production for both β‐glucan and BCG‐induced trained immunity responses (Spearman ρ = 0.73, P < 1 × 10–4 and Spearman ρ = 0.66, P < 1 × 10–4, respectively; Supporting Information Fig. S1). Furthermore, we observed a high degree of inter‐individual variability in both β‐glucan and BCG‐induced trained immunity responses in both cohorts (Fig. 1C; Supporting Information Fig. S1). To understand the cause of the observed interindividual variation in trained immunity, we first tested the effect of age and sex. Age showed no correlation with induction of trained immunity (Fig. 1D), whereas male individuals showed higher trained immunity responses, as assessed by measuring TNF‐α production (Fig. 1E) and IL‐6 (data not shown) induced by either β‐glucan or BCG.

Figure 1

Inter‐individual variation in trained immunity responses. (A) Schematic outline of the in vitro training experiments. Human monocytes were incubated with culture medium (negative control) or different stimuli for 24 h (in 10% pooled human serum). After washing and a 5 day‐rest in culture medium, cells were re‐stimulated for 24 h with culture medium or LPS (10 ng/mL) and levels of IL‐6 and TNF‐α were assessed by ELISA in supernatants. (B) Log‐transformed fold changes of cytokine production in trained versus non‐trained monocytes in the 300BCG cohort. Human monocytes were incubated for 24 h with either culture medium or β‐glucan (2 μg/mL), or BCG (5 μg/mL). On day 6, cells were re‐stimulated with culture medium or LPS (10 ng/mL) for 24 h and levels of IL‐6 and TNF‐α were assessed using ELISA (C) Log‐transformed fold changes of cytokine levels in the 300BCG cohort. Green/red color corresponds to an increase/decrease in cytokine production, respectively, indicating interindividual variability in cytokine responses. Impact of (D) age and (E) sex on trained immunity responses in the 300BCG cohort. Correlations between cytokine responses and age were calculated using Spearman's rank correlation. Data were analyzed using the Wilcoxon signed‐rank test, 300BCG cohort n = 267 (B–E).

Genome‐wide QTL mapping combined with functional validation pinpoints SIGLECs and KDM4 as important mediators of trained immunity

Using an integrative genomics approach, in which QTL studies were combined with functional investigation, our aim was to prioritize and understand the role of specific genes in trained immunity responses. We first systematically investigated the impact of genetic variation on in vitro trained immunity responses to identify genetic variants that influence the magnitude of trained immunity. To increase power, we conducted a meta‐analysis of QTLs identified in two independent population‐based cohorts, the 300BCG (n = 267) and 200FG (n = 110) cohort. We first tested for correlation between single nucleotide polymorphisms (SNPs) and the magnitude of individual trained immunity response in each cohort (QTL mapping) after adjusting for the effect of age and sex. Next, we performed meta‐analysis by using the summary statistics of QTLs identified in the two independent cohorts, which revealed 20 independent genetic loci (trained‐immunity QTLs) that showed suggestive association with trained immunity responses (P < 1 × 10–5; Fig. 2 and Table 1). From these QTL loci, eight were associated with β‐glucan‐induced trained immunity responses, and 12 were associated with BCG‐induced responses. All of the QTLs were located in noncoding regions of the genome, suggesting that gene regulation induced by trained immunity‐QTLs is an important source of variation in these responses.

Figure 2

Results from meta‐analysis of trained immunity QTLs using the 300BCG and 200FG cohort.

Manhattan plots showing the QTLs influencing β‐glucan/BCG induced (A) IL‐6 and (B) TNF‐α levels. Y‐axis represents –log10 P‐values of SNPs. X‐axis shows chromosomal positions. The dashed line indicates the suggestive P threshold for association (<1 × 10–5). Depicted QTLs showed a P‐value of significance <0.05, and showed the same direction between the two cohorts (300BCG cohort n = 267, 200FG cohort n = 110).

Table 1

Twenty independent loci associated with trained immunity responses upon β‐glucan or BCG stimulation in monocytes (P < 1 × 10–5) and prioritized genes in a window of 1 MB around the QTLs are presented

Chr	Base pair	SNP	Minor allele	MAF	Z score	P‐value	Training	Cytokine	Gene(s)
1	233794603	rs6694698	A	0.23	4.95	7.42 × 10–7	BCG	IL‐6	KCNK1 a
1	114677948	rs529989	C	0.44	3.447	6.31 × 10–6	BCG	TNF‐α	SYT6 a , e
1	35253011	rs480749 c	C	0.34	4.486	7.24 × 10–6	β‐glucan	TNF‐α	GJB3 a , RP1‐34M23.5 a , GJA4 b , GJB5 b , e
2	33951513	rs10196645	T	0.33	5.072	3.93 × 10–7	β‐glucan	IL‐6	MYADML a , AC009499.1 a
2	108771986	rs855011	T	0.44	−4.715	2.42 × 10–6	β‐glucan	IL‐6	SULT1C3 a , AC019100.3 a
3	102588665	rs2590267	A	0.25	−4.495	6.94 × 10–6	BCG	IL‐6	ZPLD1 a , e , U6 a
4	187352601	rs72714110	G	0.30	4.542	5.57 × 10–6	BCG	TNF‐α	RP11‐215A19.1 a , LOC285441 a
4	41655258	rs7677930	C	0.20	−4.598	4.27 × 10–6	β‐glucan	TNF‐α	LIMCH1 a
5	73178550	rs2973528 d	T	0.41	4.569	4.90 × 10–6	BCG	IL‐6	RGNEF a , RP11‐428C6.1 a , ANKRA2 b
8	73736455	rs893891	C	0.36	−4.442	8.89 × 10–6	BCG	TNF‐α	KCNB2 a
9	118168921	rs4145473	T	0.23	−4.438	9.08 × 10–6	BCG	TNF‐α	DEC1 a
10	121516017	rs11199085	A	0.24	3.625	1.67 × 10–7	β‐glucan	TNF‐α	INPP5F a , e ,MCMBP b , e ,BAG3 b , e
10	26041292	rs2442853	C	0.10	4.102	9.86 × 10–6	β‐glucan	TNF‐α	GPR158 a , RN5S306 a
11	134426461	rs13377274	G	0.15	−4.421	9.83 × 10–6	BCG	IL‐6	LOC283177 a , AP004550.1 a
11	123141616	rs2156448	T	0.32	4.307	5.46 × 10–6	BCG	TNF‐α	CLMP a , e , U8 a
12	87299084	rs17371306	T	0.14	4.797	1.61 × 10–6	BCG	IL‐6	MGAT4C a , RP11‐202H2.1 a
13	30665705	rs598032	C	0.30	−4.508	6.56 × 10–6	β‐glucan	IL‐6	LINC00365 a , KATNAL1 a , e
16	71389182	rs8054321	A	0.39	−4.443	8.89 × 10–6	β‐glucan	TNF‐α	CALB2 a , ZNF23 b , FTSJD1 b , PHLPP2 b , e , HP b , e
19	52298815	rs11084122	A	0.27	4.48	7.47 × 10–6	BCG	IL‐6	FPR3 a , b , e , FPR2 b , ZNF613 b , SIGLEC14 b , e , LINC00085 b , SIGLEC5 b , e , ZNF577 b
20	23969193	rs2424607	T	0.12	−4.434	9.24 × 10–6	BCG	IL‐6	GGTLC1 a

Abbreviations: Chr, chromosome; MAF, minor allele frequency.

Closest annotated genes to QTL SNP.

SNP associated with trained immunity showed an eQTL effect based on publicly available eQTL datasets [14].

SNP proxy (r 2 ≥ 0.8), rs1764390, is a missense variant within GJA4 gene.

SNP proxy (r 2 ≥ 0.8), rs2973568, is a synonymous variant within RGNEF gene.

Gene displays a different transcriptional and/or epigenetic profile upon in vitro training with β‐glucan [15] or in vivo BCG vaccination [6].

Next, we prioritized genes within a window of 1 MB around the 20 trained‐immunity QTLs with a potential important role in trained immunity by following three approaches. First, we made use of the largest expression‐QTL (eQTL) study in blood through the eQTL Consortium [14]. Second, we tested whether any of the 20 trained‐immunity QTLs were in linkage disequilibrium (LD) with variants that alter the protein‐coding sequence of genes (both synonymous and nonsynonymous), and third, we explored whether the closest annotated genes to the identified trained‐immunity QTLs and the 16 eQTL genes from these loci (Table 1) showed epigenetic changes upon the induction of trained immunity using two previously generated datasets of epigenetic changes in human monocytes that were trained in vitro by β‐glucan (GEO dataset: GSE85245) [14] or in vivo by BCG vaccination (GEO dataset: GSE104149) [6].

A strong association among the 20 trained‐immunity QTL loci was found on chromosome 19 at rs11084122 (P < 7.47 × 10–6) that influences the fold change of IL‐6 production upon BCG‐induced trained immunity responses. This SNP showed an effect on the expression of SIGLEC5 and SIGLEC14 genes (Table 1). SIGLEC5 and SIGLEC14 code for the paired Siglec‐5 and‐14 receptors belonging to the family of Siglecs. Given that Siglecs are key regulators of immune‐cell signaling and are considered attractive therapeutic targets [15], we prioritized these genes for further functional studies on their role in trained immunity.

The role of SIGLECs genes as regulators of trained immunity responses

Siglec‐5 is an inhibitory receptor that dampens the immune response [16] and was found to display reduced levels of the activation mark H3K4me3 upon β‐glucan training as compared to control (Supporting Information Table S1 and Fig. S3A). Siglec‐14 has a ligand‐binding domain almost identical to Siglec‐5 but, in contrast to Siglec‐5, leads to the activation of the inflammatory response [17, 18]. In line with a heightened inflammatory profile in trained cells, we observed enhanced levels of H3K4me3 and H3K27ac in β‐glucan‐trained monocytes at SIGLEC14 (at ‐5491 bp and ‐20,552 bp for H3K37ac, and 1605 bp for H3K4me3) 24 h after stimulation as compared to RPMI control cells, indicative of increased transcription (Fig. 3A and B; and Supporting Information Table S1 and Fig. S3A).

Figure 3

The role of . (A) Heatmap of the changes in H3K27ac in β‐glucan‐trained monocytes compared to control at 4 h, 1 day and 6 days after start of incubation at the sites of suggestive genes (P < 0.05, fold change >1.5, n = 2 donors, 1 experiment) (GEO: GSE85245). (B) Average H3K27ac signal at enhancer sites of SIGLEC14 in β‐glucan‐trained monocytes and control (n = 2). (C and D) Monocytes from seven donors from an independent BCG vaccination study (GEO: GSE104149) were analyzed by (C) ChIP‐seq and (D) RNA‐seq before BCG vaccination (baseline) and 1 month after BCG (1 experiment). Responders (n = 3) and nonresponders (n = 4) were defined based on the level of the trained immunity response (1 experiment). (E) Volunteers from the 300BCG cohort were vaccinated with BCG (n = 278) and PBMCs were isolated and stimulated ex vivo with heat‐killed Staphylococcus aureus before vaccination and 90 days after vaccination. Fold increases (compared to baseline) of IL‐1β were taken as a measurement of trained immunity responses. (F) Regional association plot showing that SNPs in a window of 250 kb around SIGLEC14 influence in vivo trained immunity responses induced by the BCG vaccine in volunteers from the 300BCG cohort (P < 1 × 10–4). Data in (C) and (D) are shown as mean + SEM (independent t‐test).

Next, we made use of previously published epigenetic data in which human healthy volunteers were vaccinated with BCG and subsequently exposed to an attenuated yellow fever virus vaccine strain [6]. In this study, volunteers were divided into a group of responders (n = 3), who showed low levels of viremia (CT > 36, indicative of a strong trained immunity response) and non‐responders (n = 4) who showed high levels of viremia (CT < 36). Responders showed an increase in H3K27ac signal at SIGLEC14 1 month after BCG vaccination as compared to nonresponders (P = 0.05) (Fig. 3C). A similar pattern, although not statistically significant, was observed for Siglec‐14 expression levels (Fig. 3D). In contrast, no profound differences in the fold change of H3K27ac levels at SIGLEC5 upon BCG vaccination were observed between responders and nonresponders (Supporting Information Fig. S3B). These findings suggest that increased activity of Siglec‐14, which leads to enhanced inflammatory responses, is an important factor influencing the magnitude of trained immunity responses.

In addition to our in vitro model to assess trained immunity responses, we investigated whether genetic polymorphisms in SIGLEC14 affect trained immunity responses induced in vivo by BCG vaccination. To this end, we stimulated PBMCs isolated from healthy volunteers (n = 278) ex vivo with Staphylococcus aureus before and 3 months after BCG vaccination (Fig. 3E). As depicted in Fig. 3F, SNPs in a window of 250 kb around SIGLEC14 strongly influence the fold change of IL‐1β production upon ex vivo exposure to S. aureus (P < 1 × 10–4). This observation further supports an important role of SIGLEC14 in the induction of trained immunity responses.

Histone demethylase KDM4 family members mediate the induction of trained immunity through changes in H3K9me3 levels

Given that the induction of trained immunity is mediated by epigenetic modifications, epigenetic enzymes may be another attractive therapeutic target to modulate trained immunity responses. Of note, promising drug targets that regulate chromatin structure are the histone demethylases, which catalyze the removal of methyl marks from histone lysine residues [19]. Therefore, to identify genes that encode for epigenetic enzymes that potentially impact the induction of trained immunity, we extracted trained‐immunity QTLs at a P threshold of < 1 × 10–4 within a window of 250 kb around known genes encoding for epigenetic enzymes (Supporting Information Table S2). Within the family of histone demethylases, we found that members of the KDM4 family showed the strongest association with trained immunity responses (Fig. 4A; Supporting Information Table S2), suggesting that KDM4 demethylases play a role in these responses. To further verify their role in trained immunity, we performed functional validation using in vitro follow‐up experiments as described below.

Figure 4

The role of the histone lysine demethylase (KDM) subfamily 4 in trained immunity. (A) Heatmap of the P‐values of association between SNPs mapped to genes coding for epigenetic regulators of the KDM4 subfamily and the magnitude of cytokine production upon trained immunity. Related to Supporting Information Table S2. (B) Monocytes were stimulated for 24 h with culture medium, β‐glucan or BCG in the presence or absence of the KDM4 inhibitor JIB‐04 (100nM). On day 7, after restimulation with LPS for 24 h, IL‐6 and TNF‐α were measured in the supernatants. (C) Effect of JIB‐04 (100 nM) on ROS production by BCG. Monocytes were trained for 24 h with BCG and normalized cell numbers were restimulated using zymosan. The amount of luminescence was measured for 1 h. RLU/s, relative light units per second. (D) Monocytes were trained with β‐glucan in the presence or absence of the KDM4 inhibitor, and lactate production was assessed on day 6. (E) Human monocytes were stimulated for 24 h with culture medium or β‐glucan in the presence or absence of JIB‐04 (100 nM). Chromatin was fixated on day 6 and ChIP‐qPCR was performed. H3K9me3 was determined at promoters of IL6, TNFA, and IL1B. (F) Human monocytes were stimulated for 24 h with culture medium or GM‐CSF (5 ng/mL). On day 6, cells were counted and 5 × 105 cells were restimulated with control medium or LPS (10 ng/mL). After 24 h, levels of IL‐6 and TNF‐α were assessed. (G) Monocytes were trained with GM‐CSF in the presence or absence of the KDM4 inhibitor JIB‐04 (100 nM). On day 7, after restimulation with LPS for 24 h, IL‐6 and TNF‐α were measured in the supernatants. (H) Human monocytes were incubated for 24 h with culture medium or GM‐CSF (5 ng/mL). Chromatin was fixated on day 6 and ChIP‐qPCR was performed. H3K9me3 was determined at promoters of IL6, TNFA, IL1B, MTOR, HK2, PFKP, and PFKFB3. Data are presented as mean + SEM, n = 10 (B), n = 7 (C), n = 6 (D), n = 3 (E), n = 8 (F), n = 10 (G), n = 6 (H) from four (B, C, F, and G), two (D and E), and three (H) independent experiments. *P < 0.05, **P < 0.005 (Wilcoxon signed‐rank test).

To validate the role of KDM4 demethylases in trained immunity responses, we first inhibited KDM4 proteins by the small molecule JIB‐04 [20]. We observed significantly decreased trained immunity responses induced by either β‐glucan or BCG (Fig. 4B). To obtain further insight into the effect of KDM4 on the functional phenotype of trained monocytes, we measured the production of ROS by trained monocytes upon stimulation with opsonized zymosan. We observed that ROS production in BCG‐trained monocytes was significantly reduced when monocytes were trained while KDM4 activity was inhibited (Fig. 4C). We next assessed the effect of the presence of JIB‐04 during training in the first 24 h on glycolysis via measuring lactate production on day 7. Production of lactate was significantly reduced when monocytes were trained by β‐glucan in the presence of JIB‐04, indicating a role of KDM4 as a modifier of the genes important for the induction of glycolysis in trained monocytes (Fig. 4D).

Given that KDM4 proteins are known for their ability to promote gene transcription by removing the repressive histone modification H3 lysine 9 trimethylation [21], we hypothesized that inhibition of trained immunity induction in the presence of the KDM4 inhibitor might be due to increased levels of the repressor mark H3K9me3. In agreement with our hypothesis, ChIP‐qPCR experiments revealed elevated H3K9me3 levels in monocytes trained in the presence of the inhibitor JIB‐04 (Fig. 4E). Together, these data indicate the importance of changes in histone trimethylation at H3K9 in trained immunity as well as a key role of KDM4 proteins, which are important for demethylation of H3K9me3 [22].

In addition to exogenous microbial ligands (β–glucan and BCG) that we used in our experiments described above to investigate the broad role of H3K9me3 and KDM4 demethylases in trained immunity, we assessed their role in a model in which human monocytes are stimulated with an endogenous ligand known to mediate trained immunity in mice: the GM‐CSF [23]. Training of monocytes with GM‐CSF enhanced cytokine production when monocytes were re‐stimulated with LPS or other PRR ligands (Pam3Cys, Poly I:C) (Fig. 4F; Supporting Information Fig. S4 and S5). Next, we investigated whether modulating KDM4 would abrogate the training effect induced by GM‐CSF. Similar to our findings observed upon β‐glucan and BCG‐induced trained immunity, incubation of monocytes with the KDM4 inhibitor JIB‐04 significantly inhibited the GM‐CSF‐induced trained immunity response, indicating the importance of KDM4 demethylases in trained immunity induced by an endogenous stimulus (Fig. 4G). To study whether GM‐CSF training is associated with changes in histone trimethylation at H3K9, we determined levels of H3K9me3 by ChIP‐qPCR at day 6 after GM‐CSF induced training. We observed decreased levels of H3K9me3 at the promoter sites of pro‐inflammatory genes and at genes involved in glycolysis upon GM‐CSF training. This finding was in agreement with our observations upon β‐glucan‐induced training (Fig. 4H). Lastly, we tested whether genetic variation around KDM4 genes regulates in vivo trained immunity responses induced by BCG vaccination. For this, we mapped and extracted trained‐immunity QTLs within a window of 250kb around KDM4 genes and found that they strongly influence in vivo trained immunity responses induced by BCG vaccination (P < 1 × 10–5; Supporting Information Table S3). These data further support the role of KDM4 demethylases in the induction of trained immunity responses.

Trained‐immunity QTLs are enriched for genes involved in the PI3K‐Akt and ECM‐receptor interaction pathways

Next, we investigated the molecular pathways underlying trained immunity responses at a genetic level. For this, we performed a pathway enrichment analysis using genes from QTL loci that show a P < 1 × 10–3 [24]. Pathway enrichment analysis revealed that trained‐immunity QTLs influencing β‐glucan‐induced trained immunity are significantly enriched for genes involved in the PI3K‐Akt signaling pathway (P = 2.83×10–5, FDR = 0.01; Fig. 5A and B; Supporting Information Table S4 and S5). Given that this pathway has been previously shown to mediate a metabolic switch from oxidative phosphorylation to glycolysis, thereby increasing lactate production [10, 12], we measured levels of lactate in the supernatant of trained monocytes in the 200FG cohort. We observed increased lactate production in trained monocytes (Fig. 5C), indicative of an increased glycolytic rate. Furthermore, we found a positive association between the concentration of lactate released by monocytes and an increase in cytokine production capacity (Fig. 5D), indicating that the metabolic switch to glycolysis is an important determinant of the magnitude of individual trained immunity responses. In addition, we found that QTLs influencing β‐glucan‐induced trained immunity were enriched for the relaxin signaling pathway (P = 2.35 × 10–4, FDR 0.02). As for the QTLs influencing BCG‐induced trained immunity, they were found to be significantly enriched for the ECM‐receptor interaction pathway (Fig. 5A and B). Lastly, genes involved in TLR signaling and NK mediated cytotoxicity pathways appeared to be important for BCG‐induced trained immunity responses.

Figure 5

(A) Pathway enrichment analysis of genes from QTL loci associated with trained immunity responses (P < 1 × 10–3). KEGG was used as a database resource of the pathways. (B) Selection of genes associated with their designated pathway (an entire list of genes involved in these pathways can be found in Supporting Information Table S5). (C) Fold change in lactate production in trained versus nontrained cells (200FG cohort, mean + SEM, n = 113, ***P < 0.0001, Wilcoxon signed‐rank test). (D) Spearman correlation plot showing the relationship between the magnitude of trained immunity responses and fold change in lactate production (n = 113).

Discussion

In this study, we used an integrative genomics approach to reveal host genetics factors and pathways that influence the induction of trained immunity responses, with a specific interest on epigenetic regulators that mediate these processes. In line with previous findings in BCG‐vaccination studies [6], we found a significant interindividual variability in trained immunity responses in vitro. Trained immunity responses were more pronounced in men compared to women, which may reflect the different heterologous effects of vaccines in boys versus girls [25]. Of note, age did not seem to influence trained immunity responses in our cohorts. Interestingly, a previous study showed intact innate immune responses of monocytes in the elderly [26], suggesting that targeting the induction of trained immunity might be a potent preventive measure against infections in the elderly, who are known to respond poorly to classical vaccination [27]. A recent study in Greece demonstrated protective heterologous effects against infections by BCG vaccination in the elderly [28], and this avenue of improving host defense in the elderly warrants future investigations.

We performed a meta‐analysis of QTLs identified in the 300BCG and 200FG cohorts on a genome‐wide scale and prioritized novel genes that belong to the family of Siglecs, Siglec‐5 and Siglec‐14, and KDM4 demethylases to understand better their role in trained immunity responses. The Siglec family of receptors recognizes sialic acids and plays a significant role in immunity [29]. Most Siglec receptors, such as Siglec‐5, are immunosuppressive [16], although a limited number of Siglecs, including Siglec‐14, induces the activation of an inflammatory response [17, 18]. Interestingly, differential expression of Siglec‐14 has been shown to play a role in the incidence of premature delivery in Group B Streptococcus (GBS)‐positive mothers [18], COPD exacerbation [30] and susceptibility to tuberculosis [31]. In this study, we show that polymorphisms in SIGLEC14 influence the induction of in vitro as well as in vivo induced trained immunity responses. Siglec‐14 has been identified as a positive regulator of NLRP3 inflammasome activation, thereby increasing IL‐1β production in human macrophages [32]. Given the crucial role of the NLRP3 inflammasome and IL‐1β production in trained immunity [6, 8, 33], we speculate that increased IL‐1β production as a result of increased Siglec‐14 expression is a possible mechanism underlying trained immunity responses, which needs further exploration. Furthermore, we found that the induction of trained immunity was associated with epigenetic modifications at the promoter region of SIGLEC14, indicative of increased transcription of Siclec‐14 upon BCG vaccination. How epigenetic changes are selectively directed to specific positions in the genome needs to be further investigated. However, evidence indicates the involvement of dynamic loop‐mediated regulation, topologically associated domains (TADs) and immune gene‐priming long noncoding RNAs (IPLs) [34].

Furthermore, we provided functional evidence for a key role of the KDM4 family of demethylases in the induction of trained immunity via modulation of methylation at H3K9. We showed that inhibition of KDM4 impairs glycolysis in β‐glucan‐induced trained immunity. This observation suggests a link between the induction of epigenetic changes and alterations in cell metabolism as it has been previously shown [35]. KDM4C selectively interacts with HIF‐1α and functions as a co‐activator to stimulate transcription of HIF‐1 target genes, which encode proteins required for metabolic reprogramming [36]. HIF‐1α is also an essential transcription factor in β‐glucan‐induced trained immunity as it is responsible for the metabolic shift induced by β‐glucan. Of note, a recent study on trained immunity induced by a Western‐type diet in mice revealed increased expression of KDM4A‐D in myeloid progenitor cells in trained mice compared to control ones [33]. It is important to mention that KDM4 demethylases and Siglecs are considered attractive therapeutic targets in the field of oncology and other immune cell‐mediated diseases [37, 38]. Therefore, targeting these factors warrants further investigation as it may benefit the prevention and treatment of inflammatory diseases in which trained immunity plays an important role, such as atherosclerosis and gout.

In addition, trained‐immunity QTLs influencing β‐glucan‐induced trained immunity responses found to be significantly enriched for genes involved in the PI3K‐Akt signaling pathway (P = 2.83×10–5, FDR = 0.01) followed by the relaxin signaling pathway (P = 2.35 × 10–4, FDR 0.02). The importance of PI3K‐Akt signaling pathway in β‐glucan induced trained immunity is supported by a previous in vitro study that demonstrated increased phosphorylation of Akt in β‐glucan trained cells as compared to control ones [12]. In addition, inhibition of this pathway using the PI3K inhibitor wortmannin abrogated the trained immunity phenotype [12]. Interestingly, the relaxin signaling pathway has been shown to increase PI3K and Akt activity in THP‐1 cells [39]. The BCG‐induced trained‐immunity QTLs were found to be significantly enriched for the ECM‐receptor interaction pathway followed by the TLR signaling pathway and NK‐mediated cytotoxicity. The enrichment in the ECM‐receptor interaction pathway is in line with recent data that highlight the importance of the ECM in the innate immune response to microbial pathogens [40].

Due to the relatively limited size of our cohorts, we followed an integrative genomics approach in which we used different molecular datasets to prioritize genes and pathways that may play an important role in trained immunity responses. Of the prioritized genes, we provided genetic and functional evidence for Siglec‐5, Siglec‐14, and KDM4 as important mediators of trained immunity responses. In addition, due to the limited size of our cohorts, we cannot exclude the possibility that additional trained‐immunity QTLs were missed. Future studies should include larger cohorts to increase statistical power and detect more subtle effects and robust associations. Such large cohorts in combination with an integrative approach would allow gaining a more thorough understanding of the mechanisms that regulate trained immunity responses.

In summary, we present a framework in which we integrated genetic, epigenetic, and functional validation data to shed light on the regulation of trained immunity responses. This data integration revealed a novel role of SIGLECs and KDM4 genes on trained immunity responses. By studying the processes underlying the heterogeneity of trained immunity responses, the present study prioritizes and identifies important candidate genes and pathways that are implicated in trained immunity that could be novel targets for diagnostic and therapeutic purposes.

Materials and methods

Experimental model and subject details

The study was performed in two independent cohorts of healthy individuals of Western European descent from the Human Functional Genomics Project (http://www.humanfunctionalgenomics.org). The 300BCG cohort consists out of 325 adults from the Netherlands (44% males and 56% females, age range 18–71 year). The 200FG comprises 241 volunteers between 23 and 73 years of age and consisted of 23% females and 77% males.

Monocyte isolation and in vitro training experiments

Venous blood from the cubital vein of volunteers from the 200FG and 300BCG cohorts was drawn into 10 mL EDTA tubes (Monoject).  Peripheral blood mononuclear cell (PBMC) isolation was performed by density centrifugation of blood diluted 1:1 in pyrogen‐free PBS over Ficoll‐Paque (GE healthcare, UK). Cells were washed twice in PBS and suspended in RPMI culture medium (Roswell Park Memorial Institute medium, Invitrogen, CA, USA) supplemented with 50 μg/mL gentamicin (Centrafarm), 2 mM glutamax (Gibco), and 1mM pyruvate (Gibco). Next, the cells were counted in a Coulter counter (Beckman Coulter, Pasadena) and the concentration was adjusted to 5 × 106 cells/mL. A total of 5 × 105 PBMCs in a 100 μL volume were added to flat‐bottom 96‐well plates (Corning, NY, USA) and cells were let to adhere for 1 h at 37°C. Afterward cells were washed thrice with warm PBS to increase the purity of the adherent monocytes. Monocytes were trained according to the in vitro training protocol as described previously [13, 41]. Briefly, monocytes were incubated either with culture medium only as a negative control, 2 μg/mL of β‐1,3‐(d)‐glucan or 5 μg/mL BCG (300BCG cohort: BCG‐Bulgaria, Intervax, 200FG cohort: BCG Denmark 1331), for 24 h at 37°C (all in the presence of 10% pooled human serum). Cells were washed once with 200 μL warm PBS and incubated for 5 days in culture medium with 10% human pooled serum. Medium was changed once on day 3. On day 6, cells were re‐stimulated for 24 h with either 200 μL RPMI or LPS 10 ng/mL (serotype 055: B5; Sigma). Where indicated, monocytes were restimulated with 10 μg/mL Pam3Cys (EMC microcollections, Tübingen, Germany; L2000) or 50 μg/mL Poly I:C (Invivogen). After 24 h of stimulation, supernatants were collected and stored at ‐20°C until assayed (300BCG: n = 267, 200FG: n = 110).

For validation experiments, Percoll isolation of monocytes was performed. A total of 150–200 × 106 PBMCs were layered on top of a hyperosmotic Percoll solution (48.5% Percoll (Sigma), 41.5% sterile H2O, 0.16 M filter‐sterilized NaCl) and centrifuged for 15 min at 580 × g. Next, the interphase layer was isolated, cells were washed with cold PBS and resuspended in RPMI. Percoll‐isolated monocytes were let to adhere to polystyrene flat bottom plates for 1 h at 37°C; subsequently cells were washed with warm PBS to yield maximal purity. In some experiments, monocytes were isolated from PBMCs by negative selection using a MACS system (Miltenyi Biotech, Bergisch Gladbach, Germany), according to the manufacturer's protocol (>96% purity of monocytes). Where indicated, cells were pre‐incubated (before stimulation) for 1 hour with 100nM JIB‐04 (Tocris). The concentration used was the highest non‐toxic concentration demonstrated in pilot experiments and did not result in differences in cell counts as compared to the control condition (data not shown).

Cytokine and lactate measurements

Cytokine production was measured in supernatants using commercial ELISA kits for human IL‐6 and IL‐1β (Sanquin, Amsterdam, the Netherlands) and TNF‐α (R&D systems, MN, USA), in accordance with the manufacturer's instructions. Levels of lactate were determined using a Lactate Fluorometric Assay Kit (Biovision, CA, USA).

Cytotoxicity assay

Cell survival was analyzed using the CytoTox 96® Non‐Radioactive Cytotoxicity Assay (Promega, Leiden, the Netherlands), in accordance with the manufacturer's instructions. This assay measures lactate dehydrogenase (LDH) in a cell‐free supernatant. As LDH is a stable cytosolic enzyme only released into the supernatant after cell lysis, it is used to detect cell death.

ROS measurements

The production of ROS was determined by a luminal‐enhanced luminescence assay. Upon detachment and counting of monocytes, a total of 1 × 105 monocytes in a 200 μL volume were added to white 96‐well plates (Corning). Cells were left either unstimulated or stimulated with 1 mg/mL serum‐opsonized zymosan (from S. cerevisiae, Sigma). Next, 145 g/ml luminol was added, and chemiluminescence was measured every 142 s for 1 h. All measurements were performed at least in duplicate.

Assessment of histone marks

For ChIP‐qPCR experiments shown in Fig. 5, monocytes were trained in vitro as described above. On day 6, cells were detached from the plate with Versene and fixed in methanol‐free 1% formaldehyde. Next, cells were sonicated and IP was performed using antibodies against H3K9me3 (Diagenode, Seraing, Belgium). DNA was processed further for qPCR analysis using the SYBR green method. Samples were analyzed by a comparative Ct method according to the manufacturer's instructions. GAPDH was used as a negative control and ZNF as a positive control in experiments assessing H3K9me3. Primer sequences are listed in Supporting Information Table S6.

DNA isolation

DNA was isolated from EDTA venous blood using the Gentra Pure Gene Blood kit, in accordance with the manufacturer's instructions (Qiagen, Venlo, the Netherlands).

Genotyping, quality control, and imputation of the 300BCG cohort

DNA samples of individuals (n = 325) were genotyped using the commercially available SNP chip, Infinium Global Screening Array MD version 1.0 from Illumina. Opticall 0.7.0 with default settings was used for genotype calling [42]. Genetic variants with a call rate ≤ 0.01 were excluded, as were variants with a Hardy‐Weinberg equilibrium (HWE) ≤ 0.0001, and minor allele frequency (MAF) ≤ 0.001. The strands and variant identifiers were aligned to the 1000 Genome reference panel using Genotype Harmonizer [43]. One sample was excluded from the pre‐imputed dataset due to high relatedness (PI_HAT = 0.99). As measure of relatedness, we used the estimates of PI_HAT (proportion identity‐by‐decent) as calculated using plink version 1.90b [44]. We then imputed the samples on the Michigan imputation server using the human reference consortium (HRC r1.1 2016) as a reference panel. Data were phased using Eagle version 2.3, and we filtered out genetic variants with imputation quality score R 2 < 0.3 [45]. We identified and excluded 17 genetic outliers. Of these outliers, five individuals showed high relatedness in pair with other individuals (0.4 < PI_HAT < 0.5), and 12 samples were removed as ethnic outliers based on multidimensional scaling analysis (MDS; Supporting Information Fig. S6A). We selected 4,296,841 SNPs with MAF 5% for follow‐up QTL mapping.

Genotyping, quality control, and imputation of the 200FG cohort

DNA samples of individuals (n = 241) were genotyped using the commercially available SNP chip, Illumina HumanOmniExpressExome‐8 v1.0. Opticall 0.7.0 with default settings was used for genotype calling [42]. Quality control steps and imputation of the genetic data were done similarly as for the 300BCG cohort. We identified and excluded seven samples as ethnic outliers based on (MDS) (Supporting Information Fig. S6B). We selected 4,296,378 SNPs with MAF 5% for follow‐up QTL mapping.

QTL mapping

In vitro trained immunity responses in the 300BCG cohort

Both genotype and cytokine data on trained immunity responses was obtained for a total of 267 individuals. Three samples were excluded due to medication use (of which one was identified as a genetic outlier), and one sample due to onset of type 1 diabetes during the study. Few cytokine measurements upon low β‐glucan stimulation or, under the detection limit (<39 pg/mL for TNF‐α and <46.88 pg/mL for IL‐6) were set to missing. First, the fold change of cytokine production between trained and nontrained cells was taken as a measurement for the magnitude of the trained immunity response. To check for normality of the fold change of cytokine production, we followed a visual inspection of the fold change of cytokine production (Supporting Information Fig. S7) using both raw and log‐transformed data. Fold change of cytokine production followed non‐Gaussian distribution before data transformation. The individual trained immunity response was measured as the fold change in cytokine production in trained monocytes as compared to nontrained monocytes (RPMI control). Following quality check for cytokine distribution and after excluding genetic outliers, we mapped the log‐transformed fold changes of cytokine production to genotype data using a linear regression model with age and sex as covariates to correct the distributions of fold change of cytokine production. R‐package Matrix‐eQTL was used for cytokine QTL mapping [46]. We used a cutoff of P < 1 × 10–5 to identify suggestive QTL associations affecting trained immunity responses.

In vivo trained immunity responses in the 300BCG cohort

Individuals from the 300BCG cohort were vaccinated with 0.1 mL of BCG (BCG vaccine strain Bulgaria; Intervax, Canada). PBMCs were isolated as described above and stimulated ex vivo with 5 × 106 CFU/mL heat‐killed Staphylococcus aureus before vaccination, 14 days and 90 days after vaccination. IL‐1β production was measured after 24 h in supernatants and the fold change in cytokine production (after vaccination compared to baseline) was used as a measurement of the magnitude of the trained immunity response. Both genotype and cytokine data on trained immunity responses was obtained for a total of 296 individuals. Genetic outliers and outliers due to medication use and onset of type 1 diabetes during the study as described above were excluded. In addition, 18 evening vaccinated individuals were excluded, resulting in 278 samples before QTL mapping. The magnitude of the trained immunity response was assessed by the fold change in cytokine production (IL‐1β) after BCG vaccination as compared to levels before vaccination. The fold change of cytokine production was then log‐transformed and mapped to genotype data using a linear regression model with age and sex as covariates. R‐package Matrix‐eQTL was used for cytokine QTL mapping [46].

In vitro trained immunity responses using the 200FG cohort

Both genotype and cytokine data were obtained for a total of 110 samples. Following quality check for cytokine distribution (Supporting Information Fig. S8) and after excluding genetic outliers, both genotype and data on trained immunity responses were obtained for a total of 103 individuals. The same QTL mapping method used in 300BCG cohort was applied to 200FG samples.

Meta‐analysis of trained immunity QTLs between 200FG and 300BCG cohort

To jointly test the effect of genetic variation on the magnitude of trained immunity, we used the following meta‐analysis software, METAL. We tested for heterogeneity to decide whether observed effect sizes are homogeneous across samples and excluded QTLs with a significant P‐value of heterogeneity < 0.05.

Pathway enrichment analysis

To investigate whether specific pathways are over‐represented for the genes in a window of 500 kb around the trained‐immunity QTLs (P < 1 × 10–3), we performed an over‐representation analysis (ORA) using the web‐based gene set enrichment analysis toolkit WebGestalt [24]. Results were reported using KEGG as a database, and Benjamini–Hochberg was used to control the false discovery rate (FDR). Adjusted P‐values < 0.05 were considered as being significant.

Statistics

Statistical analyses and visualization were performed in R version 3.3.2. Unsupervised hierarchical clustering was performed using Spearman's correlation as the measure of similarity. R‐package Matrix‐eQTL, was used for QTL mapping, where linear model was applied with age and sex included as co‐variables. In vitro training experiments for functional validation studies were analyzed using a Wilcoxon signed‐rank test performed on GraphPad Prism 5.0 software (GraphPad). The data are expressed as mean + SEM unless mentioned otherwise. P‐values ≤ 0.05 were considered significant. In figures, asterisks denote statistical significance ∗ P < 0.05; ∗∗ P < 0.01; ∗∗∗ P < 0.001.

Conflict of interest

The authors declare no commercial or financial conflict of interest.

Ethics approval and consent to participate

Written informed consent was received from participants prior to inclusion in the study. Cohort studies were approved by the Ethical Committee of Radboud University Nijmegen, the Netherlands (2011/399) (200FG cohort) and (no. NL58553.091.16) (300BCG cohort). Experiments were conducted according to the principles expressed in the Declaration of Helsinki.

Author contributions

L.A.B.J., M.G.N., and S.J.C.F.M.M. were associated with study conceptualization; S.J.C.F.M.M., Y.L., V.A.C.M. K., M.O., R.J.R., O.B.B., V.P.M., S.P.S., R.A.G., V.M., J.H.P., C.H., S.K., L.G., V.K., L.C.J.d.B., and B.N. were associated with the methodology. S.J.C.F.M.M., Y.L., J.H.P., and V.K. were associated with the investigations; S.J.C.F.M.M. wrote the original draft; Y.L., V.M. V.A.C.M.K., L.C.J.d.B., J.H.P., N.P.R., R.v.C., R.J.X., C.W., V.P.M., L.A.B.J., and M.G.N. reviewed and edited the final manuscript; L.A.B.J. and M.G.N. supervised the study.

Availability of data and materials

The datasets used and/or analyzed during the current study are available from the corresponding author on reasonable request.

Peer review

The peer review history for this article is available at https://publons.com/publon/10.1002/eji.202149577

bacillus Calmette‐Guérin

peripheral blood mononuclear cell

quantitative trait locus

Supporting Information

Click here for additional data file.

Supporting Information

Click here for additional data file.