An interaction map of circulating metabolites, immune gene networks, and their genetic regulation

RESEARCH Open Access An interaction map of circulating metabolites, immune gene networks, and their genetic regulation Artika P. Nath 1,2 , Scott C. Ritchie 2,3 , Sean G. Byars 3,4 , Liam G. Fearnley 3,4 , Aki S. Havulinna 5,6 , Anni Joensuu 5 , Antti J. Kangas 7 , Pasi Soininen 7,8 , Annika Wennerström 5 , Lili Milani 9 , Andres Metspalu 9 , Satu Männistö 5 , Peter Würtz 7,10 , Johannes Kettunen 5,7,8,11 , Emma Raitoharju 12 , Mika Kähönen 13 , Markus Juonala 14,15 , Aarno Palotie 6,16,17,18 , Mika Ala-Korpela 7,8,11,19,20 , Samuli Ripatti 6,21 , Terho Lehtimäki 12 , Gad Abraham 2,3,4 , Olli Raitakari 22,23 , Veikko Salomaa 5 , Markus Perola 5,6,9 and Michael Inouye 1,2,3,4* Abstract Background: Immunometabolism plays a central role in many cardiometabolic diseases. However, a robust map of immune-related gene networks in circulating human cells, their interactions with metabolites, and their genetic control is still lacking. Here, we integrate blood transcriptomic, metabolomic, and genomic profiles from two population-based cohorts (total N = 2168), including a subset of individuals with matched multi-omic data at 7-year follow-up. Results: We identify topologically replicable gene networks enrichedfordiverseimmunefunctions including cytotoxicity, viral response, B cell, platelet, neutrophil, and mast cell/basophil activity. These immune gene modules show complex patterns of association with 158 circulating metabolites, including lipoprotein subclasses, lipids, fatty acids, amino acids, small molecules, and CRP. Genome-wide scans for module expression quantitative trait loci (mQTLs) reveal five modules with mQTLs that have both cis and trans effects. The strongest mQTL is in ARHGEF3 (rs1354034) and affects a module enriched for platelet function, independent of platelet counts. Modules of mast cell/basophil and neutrophil function show temporally stable metabolite associations over 7-year follow-up, providing evidence that these modules and their constituent gene products may play central roles in metabolic inflammation. Furthermore, the strongest mQTL in ARHGEF3 also displays clear temporal stability, supporting widespread trans effects at this locus. Conclusions: This study provides a detailed map of natural variation at the blood immunometabolic interface and its genetic basis, and may facilitate subsequent studies to explain inter-individual variation in cardiometabolic disease. Background Over the past decade increasing evidence has implicated inflammation as a probable causal factor in metabolic and cardiovascular diseases. Consequently, research has begun to focus on the interplay between immunity and metabolism, or immunometabolism. While it is involved in diverse pathophysiologies, immunometabolism is particularly relevant to diseases of immense global health burden, such as type 2 diabetes (T2D) and atherosclerosis. For T2D, immune overactivation in adipose tissue has been implicated as a key driver [1, 2]. Studies have shown that macrophage infiltration and subsequent overexpression of proinflammatory cytokines, such as TNF-α, in adipose tissues is associated with insulin resistance [1, 2]. Moreover, evidence for metabolic inflammation has been shown in other tissues where, in blood, elevated glucose and free fatty acid levels potentiate IL1β-mediated destruction of pancreatic ß cells and subsequent T2D progression [3–5]. While circulating metabolites are known to be associated with cardiovascular disease [6], inflammation is an increasingly recognized factor in pathogenesis. In atherosclerosis, lipid-induced inflammatory response mechanisms have also been * Correspondence: [email protected] 1 Department of Microbiology and Immunology, The University of Melbourne, Parkville 3010, Victoria, Australia 2 Systems Genomics Lab, Baker Heart and Diabetes Institute, Melbourne, Victoria, Australia Full list of author information is available at the end of the article © The Author(s). 2017 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated. Nath et al. Genome Biology (2017) 18:146 DOI 10.1186/s13059-017-1279-y implicated in progression to myocardial infarction [7]. In atherogenic lesions, oxidized phospholipids are known to lead to a new macrophage phenotype [8], and cholesterol loading in macrophages promotes proinflammatory cytokine secretion [9]. Perhaps surprisingly, few large-scale studies have systematically assessed interactions between the human immune system and metabolites. Recent studies have investigated matched blood transcriptomic and metabolomic profiles to understand their interplay [10–16]. However, these studies had modest sample sizes and thus have not had the power to focus on the diverse range of immune processes that interact with circulating metabolites. Furthermore, even fewer have assessed effects of expression quantitative trait loci (eQTLs) on immune gene networks. A robust integrated map of immunometabolic relationships and their genetic regulation would provide a foundation for investigating the differential cardiometabolic disease susceptibility amongst individuals while also identifying key target interactions for mechanistic in vivo and in vitro follow-up. In this study, we present an integrated immunometabolic map using matched blood metabolomic and transcriptomic profiles from 2168 individuals from two population-based cohorts. We perform gene coexpression network discovery and cross-cohort replication to identify robust gene modules which encode immunerelated functions. Using a high-throughput quantitative NMR metabolomics platform that can separate lipids and lipoprotein sub-fractions as well as quantify a panel of polar metabolites, we identify significant interactions between immune gene modules and circulating metabolite measures. Genome-wide scans for QTLs affecting immune gene modules identify many cis and trans loci affecting module expression. Finally, we test the longterm stability of gene modules, their interactions with metabolite measures, and genetic control using a 7-year follow-up sampling of 333 individuals. Results and discussion Summary of cohorts and data We analyzed genome-wide genotype, whole blood transcriptomic, and serum metabolomics data from two population-based cohorts (“Methods”and Fig. 1). In DILGOM07, 240 males and 278 females aged 25–74 years were recruited (total N = 518). Data were available for a subset of 333 participants from DILGOM07 who were followed up after 7 years (DILGOM14). In YFS, relevant data were available for 755 males and 895 females aged 34–49 years (total N = 1650). DILGOM and YFS genotyping was performed using Illumina Human 610 and 670 arrays, respectively, with subsequent genotype imputation performed using IMPUTE2 [17] and the 1000 Genomes Phase I version 3 reference panel. For both cohorts, whole blood transcriptome profiling was performed using Illumina HT-12 arrays and serum metabolomics profiling was carried out using the same serum NMR metabolomics platform (Brainshake Ltd) [18]. Individuals on lipid-lowering medication and pregnant women were excluded from the metabolome analyses (“Methods”). Of the 159 metabolite measures analyzed, 148 were directly quantified and 11 derived (Additional file 1: Table S1). After filtering, matched transcriptome and metabolome data were available for 440 individuals in DILGOM07 and 216 of these individuals (DILGOM14) who were profiled at 7year follow-up. In YFS, 1575 individuals were available with similar data (see “Methods”for details). Robust immune gene coexpression networks from blood We first identified networks of tightly coexpressing genes in DILGOM07 and then used a permutation approach, NetRep [19], to statistically test replication patterns of density and connectivity for these networks in YFS. For module detection, we applied weighted gene coexpression network analysis (WGCNA) to all 35,422 probes in the DILGOM07 data, identifying a total of 40 modules of coexpressed genes (“Methods”). For each module, we used NetRep to calculate seven preservation statistics in the YFS, generate empirical null distributions for each of these test statistics, and calculate their corresponding Pvalues [19, 20]. A module was considered strongly preserved if the Pvalue was <0.001 for all seven preservation statistics (Bonferroni correction for 40 modules). Of the 40 DILGOM07 modules, 20 were strongly preserved in YFS (Additional file 2: Table S2). For each of the 20 replicated modules, we defined core gene probes, those which are most tightly coexpressed YFS DILGOM07 DILGOM14 Study Populations Molecular Data Matched TranscriptomeMetabolome Matched TranscriptomeGenotypes N=333 N=258 N=294 Metabolome Transcriptome Blood Samples Genotypes Discover Gene Modules in DILGOM07 Network analysis Replicate Gene Modules in YFS and DILGOM14 SNPs N=518 N=440 N=515 N=1650 N=1575 N=1400 Permutation statistics GO enrichment analysis Immune-linked gene networks Immune-metabolite association Genetic drivers of immune modules Analysis Strategy Transcriptome Fig. 1 The study design. GO Gene Ontology, SNP single nucleotide polymorphism Nath et al. Genome Biology (2017) 18:146 Page 2 of 15 and thus robust to clustering parameters, using a permutation test of module membership (“Methods”; Additional file 3: Table S3). To identify modules of putative immune function, we carried out Gene Ontology (GO) biological process enrichment analysis using GOrilla for the core genes of each replicated module [21]. Significant GO terms (false discovery rate (FDR) <0.05) were then summarized into representative terms based on semantic similarity using REVIGO [22] (Additional file 4: Figure S1). A module was considered immune-related if it was significantly enriched for GO terms “immune system processes” (GO:0002376) and/or “regulation of immune system processes”(GO:0002682) in the REVIGO output. Six out of 20 modules were enriched for at least one of these terms (Additional file 5: Table S4). We also identified two additional modules which were not enriched for any GO terms but have been previously linked to immune functions related to mast cell and basophil function [13] and platelet aggregation activity [23]. The eight modules encoded diverse immune functions, including cytotoxic, viral response, B cell, platelet, neutrophil, mast cell/basophil, and general immune-related functions. Each immune module’s gene content and putative biological function is summarized in Table 1. Immune module association analysis for eQTLs and metabolite measures For each gene module, we performed a genome-wide scan to identify module QTLs (mQTLs) that regulate expression. In DILGOM07 and YFS, the module eigengene was regressed on each SNP, then mQTL test statistics were combined in a meta-analysis (“Methods”). Significant mQTLs were further examined at individual gene expression levels. A genome-wide significance level (P value <5 × 10 −8 ) was used to identify mQTLs and significant trans effects on individual gene expression (Fig. 2 and Table 2). Leukocyte and platelet counts were available for YFS and were used to test the robustness of module associations with mQTLs and metabolite measures. Six modules showed statistically significant association with platelet or leukocyte counts (Pvalue <0.05) (Additional file 6: Table S5); however, adjustment for leukocyte counts did not affect mQTL nor modulemetabolite measure associations, with the exception of the platelet module (PM) and cytotoxic cell-like module Table 1 Immune module gene content and putative biological function based on GO terms (top three shown) and literature Module Size GO terms Literature-based immune-related function of genes Cytotoxic cell-like module (CCLM) 130 (115) Immune system process Defense response Immune response Cytotoxic effectors (GZMA,GZMB,GZMM,CTSW,PRF1 [66]); surface receptors (IL2RB,SLAMF6,CD8A,CD8B,CD2,CD247,KLRD1,KLRG1 [66–68]); T and NK cell differentiation (ID2 and EOMES [69,70]), activation (ZAP70 and CBLB [71,72]), and recruitment (CX3CR1,CCL5,CCL4L2 [73]) Viral response module (VRM) 95 (88) Response to virus Type I interferon signaling pathway Response to biotic stimulus Type I interferon-induced antiviral activity (IFITM1,IFIT1,IFIT2,IFIT3,IFIT5,IFI44, IFI44L,IFI6,MX1,ISG15,ISG20,HERC5 [74,75]); viral RNA degradation (OAS1, OAS2,OAS3,OASL,DDX60 [30]); type 1 interferon-signaling pathway (IRF9, STAT1,STAT2 [76,77]) B-cell activity module (BCM) 54 (49) Immune system process Immune response B cell activation B cell surface markers (CD79A,CD79B,CD22 [33,78]); B cell activation (BANK1, BTLA,CD40,TNFRSF13B,TNFRSF13C [79]), development (POU2AF1,BCL11A, RASGRP3 [80]), migration (CXCR5,CCR6 [80,81]), and their regulation (CD83, FCER2,FCRL5 [82]); antigen presentation (HLA-DOA,HLA-DOB [83]) Platelet module (PM) a 114 (106) Coagulation Blood coagulation Cell activation Platelet receptor signaling, activation, and coagulation (GP6,GP9,ITGA2B, ITGB3,ITGB5,MGLL,MPL,MMRN1,PTK2,VCL,THBS1,F13A1,VWF,[84]); regulating platelet activity (SEPT5,TSPAN9 [85,86]) Neutrophil module (NM) a 26 (26) Killing of cells of other organism Cell killing Response to fungus Anti-microbial, -fungal, and -viral activity (DEFA1,DEFA1B,DEFA3,DEFA4, ELANE,BPI,RNASE2,RNASE3 [87–90]); neutrophil-mediated activity (AZU1, LCN2,MPO,CEACAM6,CEACAM8,OLFM4 [90,91]) and its regulation (LCN2, CAMP,OLR1 [49,92,93]) Lipid-leukocyte module (LLM) a 13 (13) Mast cell and basophil function b Mast cell and basophil related immune response and allergic inflammation (FCER1A,HDC,GATA2,SLC45A3,CPA3,MS4A3 [13,94,95]) General immune module A (GIMA) 509 (482) Immune system process Defense response Regulation of response to stimulus These modules contain genes involved in a broad range of immune processes and their regulation such as signaling; cell death; defense response to stress, inflammation, and external stimuli; leukocyte activation, migration, and adhesion General immune module B (GIMB) 74 (69) Immune response-activating signal transduction Positive regulation of immune response Activation of immune response Size refers to the number of core genes in each module and the subset of these core genes with GO term annotations are listed in parentheses. Functions were assigned to each of these modules based on GO enrichments and literature-based searches for genes in the modules. a Modules previously reported to have immune related function b The LLM was not significantly enriched for any GO term Nath et al. Genome Biology (2017) 18:146 Page 3 of 15 (CCLM) discussed below (Additional file 7: Table S6). Since we did not have cell counts available for DILGOM07, all the associations between immune modules and metabolite measures discussed below, unless otherwise noted, have not been adjusted for cell counts. Cytotoxic cell-like module CCLM was associated with 24 metabolite measures, mainly consisting of fatty acids, intermediate density lipoproteins, and C-reactive protein (CRP; Fig. 3; Additional file 8: Table S7). Of these, the average degree Fig. 2 Module and expression QTL analysis. aManhattan plot of meta-analyzed Pvalues from the DILGOM/YFS module QTL analysis. The lead SNP and its closest genes are noted. Each significant mQTL locus is colored by module. The horizontal dashed line represents genome-wide (metaPvalue <5 × 10 −8 ) significance. b–dCircular plots summarizing the individual gene associations (meta-Pvalue <5 × 10 −8 ) for the lead module QTLs in the VRM, PM, and NM. Lead SNPs and cis genes are labeled outside the ring. PM platelet module, VRM viral response module, CCLM cytotoxic cell-like module, NM neutrophil module, BCM B-cell activity module Table 2 QTLs for immune gene modules Module Top SNP Chr Hg19 pos. (Mb) Allele (minor/major) MAF (avg) Pvalue DILGOM07 (effect size) Pvalue YFS (effect size) Meta Pvalue VRM rs182710579 4 19768086 G/T 0.012 2.01 × 10 –4 (0.05) 8.10 × 10 –6 (0.02) 9.23 × 10 –9 rs151234502 7 148950168 T/C 0.012 2.59 × 10 –1 (0.01) 5.31 × 10 –9 (0.03) 2.46 × 10 –8 rs147742798 11 70947761 T/C 0.016 1.51 × 10 –3 (0.04) 1.66 × 10 –6 (0.02) 9.43 × 10 –9 BCM rs2523489 6 31348878 T/C 0.186 1.42 × 10 –1 (0.005) 5.29 × 10 –8 (0.006) 6.27 × 10 –8 PM rs1354034 3 56849749 T/C 0.284 7.11 × 10 –14 (-0.02) 1.51 × 10 –16 (-0.008) 7.35 × 10 –28 rs28367734 6 3128657 A/G 0.108 5.40 × 10 –4 (0.02) 2.02 × 10 –5 (0.006) 5.44 × 10 –8 NM rs2485364 6 159512260 C/T 0.466 1.78 × 10 –3 (0.009) 6.05 × 10 –7 (0.004) 3.93 × 10 –9 rs13297295 9 131659724 C/T 0.085 4.26 × 10 –2 (0.009) 8.39 × 10 –11 (0.01) 3.93 × 10 –11 rs140929198 20 38555870 A/G 0.031 2.98 × 10 –2 (0.03) 8.47 × 10 –9 (0.01) 1.41 × 10 –9 GIMA rs2185366 8 131342722 T/C 0.421 2.0 × 10 –2 (0.007) 1.52 × 10 –6 (0.004) 1.05 × 10 –7 Modules: VRM viral response module, BCM B-cell activity module, PM platelet module, NM neutrophil module, GIMA general immune module A. MAF minor allele frequency Nath et al. Genome Biology (2017) 18:146 Page 4 of 15 of unsaturation in fatty acids was the most significant association (meta-Pvalue = 7.23 × 10 –7 ). The immunomodulatory effects of polyunsaturated fatty acids are well characterized; for example, omega-3 fatty acids have been shown to induce cytotoxicity in in vitro cancer cell lines as well as animal models of tumor incidence and growth [24, 25]. Adjustment of the associations between CCLM and metabolite measures for leukocyte counts resulted in the gain of 38 additional associations and loss of four (creatinine, ratio of polyunsaturated fatty acids to total fatty acids, very low density lipoprotein (VLDL) particle size, and CRP) existing associations (Additional file 7: Table S6). Varying proportions of leukocyte counts can be correlated with transcription-level variation in human blood [26] but not act to confound the latter's association with phenotypes. If this is the case, then adjusting for leukocyte count in the linear regression analysis can reduce noise and thus boost statistical power to detect an association, which may explain the additional associations noted with the CCLM module. CCLM had no significant mQTLs. IDL_P L_LDL_P VLDL_D XL_HDL_P XXL_VLDL_P ACACE ACE ALA ALB APOA1 APOB APOB/APOA1 IDL_L L_LDL_L LDL_D XL_HDL_L XXL_VLDL_L BOHBUT HDL_D IDL_PL L_LDL_PL XL_HDL_PL XXL_VLDL_PL CIT CREA CRP IDL_C L_LDL_C XL_HDL_C XXL_VLDL_C DHA DHA/FA IDL_CE L_LDL_CE XL_HDL_CE XXL_VLDL_CE EST_C IDL_FC L_LDL_FC XL_HDL_FC XXL_VLDL_FC FAW3 FAW3/FA FAW6 FAW6/FA FREE_C IDL_TG L_LDL_TG XL_HDL_TG XXL_VLDL_TG GLC GLN GLOL GLY GlycA L_HDL_P M_LDL_P XL_VLDL_P HDL_C HDL_TG HDL2_C HDL3_C HIS L_HDL_L M_LDL_L XL_VLDL_L ILE L_HDL_PL M_LDL_PL L P_LDLV_LX L_HDL_C M_LDL_C C_LDLV_LX L_HDL_CE M_LDL_CE XL_VLDL_CE LA LA/FA LAC LDL_C LDL_TG LEU L_HDL_FC M_LDL_FC XL_VLDL_FC MUFA MUFA/FA L_HDL_TG M_LDL_TG XL_VLDL_TG L_VLDL_P M_HDL_P S_LDL_P L_VLDL_L M_HDL_L S_LDL_L PC PHE PUFA PUFA/FA PYR L_VLDL_PL M_HDL_PL S_LDL_PL L_VLDL_C M_HDL_C S_LDL_C REMNANT_C L_VLDL_CE M_HDL_CE S_LDL_CE SERUM_C SERUM_TG SFA SFA/FA SM L_VLDL_FC M_HDL_FC S_LDL_FC TG_PG TOT_CHO TOT_FA TOT_PG TYR L_VLDL_TG M_HDL_TG S_LDL_TG UNSAT M_VLDL_P S_HDL_P VAL VLDL_C VLDL_TG M_VLDL_L S_HDL_L M_VLDL_PL S_HDL_PL M_VLDL_C S_HDL_C M_VLDL_CE S_HDL_CE M_VLDL_FC S_HDL_FC M_VLDL_TG S_HDL_TG S_VLDL_P S_VLDL_L S_VLDL_PL S_VLDL_C S_VLDL_CE S_VLDL_FC S_VLDL_TG XS_VLDL_P XS_VLDL_L XS_VLDL_PL XS_VLDL_C XS_VLDL_CE XS_VLDL_FC XS_VLDL_TG FAW6/FAW3 LLM (124) NM (122) GIMA (98) GIMB (82) PM (56) CCLM (24) BCM (14) VRM (8) FDR adjusted P-value < 1.25 x 10 -04 1.25 x 10 -03 6.25 x 10 -03 > 6.25x10 -03 Significance threshold 1.25 x 10 -04 Lipoprotein Subclasses HDL IDL LDL VLDL Lipoprotein particle size Apolipoproteins Cholesterol Glycerides & Phospholipids Fatty Acids Small Metabolites Amino Acids Inflammatory Markers Fig. 3 Metabolite measure associations with immune gene modules. Circular heatmap of associations between individual metabolite measure and the module eigengene of each module (colored by FDR-adjusted Pvalues). Concentric circles represent modules, with numbers in parentheses denoting total number of metabolite measures associated with that module at FDR-adjusted Pvalue <6.25 × 10 –3 .NM neutrophil module, LLM lipid leukocyte module, GIMA/GIMB general immune module A/B, PM platelet module, CCLM cytotoxic cell-like module, BCM B-cell activity module, VRM viral response module. See Additional file 1: Table S1 for full metabolite descriptions Nath et al. Genome Biology (2017) 18:146 Page 5 of 15 Viral response module Three genome-wide significant mQTLs were identified for the viral response module (VRM; Fig. 2a; Table 2). The strongest mQTL, rs182710579 (meta-P value = 9.22 × 10 –9 ), is within a known lincRNA locus (RP11608O21.1) (Additional file 4: Figure S2a). Rs182710579 was a trans eQTL for three genes in the VRM (Fig. 2b; Additional file 9: Table S8). The strongest association was seen with CCL2 (meta-Pvalue = 6.78 × 10 –12 ), a proinflammatory chemokine involved in leukocyte recruitment during viral infection [27, 28]. Also, adipocytederived CCL2 is known to play an important role in obesity-associated adipose tissue inflammation and insulin resistance [29]. The next strongest mQTL, rs151234502, resides within intron 4 of the relatively unstudied ZNF212, part of a zinc finger gene cluster at 7q36 (Additional file 4: Figure S2b). Rs151234502 modulated expression of 11 VRM genes in trans (Fig. 2b; Additional file 9: Table S8). The strongest association was with OAS2 (meta-Pvalue = 8.98 × 10 –10 ), an interferon-induced gene encoding an enzyme promoting RNase L-mediated cleavage of viral and cellular RNA [30]. The third mQTL, rs147742798, was an intergenic SNP located between SHANK2 and DHCR7 at 11q13.4 (Additional file 4: Figure S2c). Rs147742798 was a trans eQTL for two genes in the VRM, BST2 and PARP9 (Fig. 2b; Additional file 9: Table S8). BST2 encodes a trans-membrane protein with interferon-inducible antiviral function [31]. Studies have previously shown induction of fatty acid biosynthesis by a range of viruses [32]. VRM was associated with eight metabolite measures, including amino acids (alanine, phenylalanine), fatty acids (omega-6 fatty acids, polyunsaturated fatty acids, saturated fatty acids, and total fatty acids), and cholesterol esters in medium VLDL (Fig. 3; Additional file 8: Table S7). Consistent with its putative role in viral response, VRM was strongly associated with CRP (meta Pvalue = 2.38 × 10 –10 ). B-cell activity module The B-cell activity module (BCM) was associated with 14 metabolite measures including CRP, histidine, lactate, apolipoprotiens, and mainly the medium high-density lipoprotein (HDL) subclass of lipoproteins (Fig. 3; Additional file 8: Table S7). The strongest association was seen with CRP (meta-Pvalue = 2.65 × 10 –8 ). Histidine was the second strongest association. This is interesting given that histidine is a substrate for histamine, and both histamine release and B-cell activity are central parts of an allergic reaction. While no mQTLs for BCM reached genome-wide significance, there was some evidence in the YFS for the MHC class I locus (Fig. 2a and Table 2). The top signal was located between HLA-B/C and MICA (rs2523489, meta-Pvalue = 6.27 × 10 –8 ; Additional file 4: Figure S3). The HLA class I region is well known to be associated with autoimmune diseases, where the role of B cells is well recognized. Rs2523489 was a trans eQTL for CD79B (meta-Pvalue = 1.16 × 10 –9 ), a subunit of the antigen-binding B-cell receptor complex [33]. Platelet module PM had the strongest mQTL of any gene module, an intronic SNP of the ARHGEF3 gene at 3p14.3 (rs1354034; meta-Pvalue = 7.35 × 10 –28 , Fig. 2a; Table 2; Additional file 4: Figure S4a). ARHGEF3 encodes a Rho guanine nucleotide exchange factor, a catalyst of Rho GTPase conversion from inactive GDP-bound to active GTP-bound form. Rs1354034 was an eQTL for the majority of genes in the PM, all of which were in trans. An intergenic SNP, rs2836773 (meta-Pvalue = 5.4 × 10 –8 ), at the HLA locus was also identified as an mQTL for PM (Additional file 4: Figure S4b). The ARHGEF3 mQTL (rs1354034) exhibited a strong trans-regulatory effect and was associated with 61 PM genes (65 unique probes) (Fig. 2c; Additional file 9: Table S8). The top trans eQTL was ITGB3 (meta-Pvalue = 5.09 × 10 –42 ), a gene encoding the β 3 subunit of the heterodimeric integrin receptor (integrin α IIb β 3 ). This integrin receptor is most highly expressed on activated platelets and plays a key role in mediating platelet adhesion and aggregation upon binding to fibrinogen and Willebrand factor [34, 35]. Our data are consistent with previous observations of the diverse trans eQTL effects of rs1354034 [23], including the putative splice-QTL effects of rs1354034 on TPM4, a significant eGene in the PM. ARHGEF3 itself is of intense interest to platelet biology. It has previously been shown that silencing of ARHGEF3 in zebrafish prevents thrombocyte formation [36]. To test whether ARHGEF3 expression had an effect on PM genes, we regressed out ARHGEF3 levels and re-ran the eQTL analysis. Adjusting for ARHGEF3 did not attenuate the trans-associations of rs1354034, suggesting either independence of downstream function for ARHGEF3 and rs1354034 or post-transcriptional modification of ARHGEF3. Previous GWAS studies have shown rs1354034 is associated with platelet count and mean platelet volume [36]; however, perhaps due to power, we found no significant relationship between platelet counts and rs1354034 in YFS. While platelet counts were positively associated with the PM (β= 0.29; Pvalue = 8.23 × 10 –30 ; Additional file 6: Table S5), the association between rs1354034 and the PM was still highly significant when conditioning on platelet counts (β=−0.33; Pvalue = 1.40 × 10 –17 ). PM displayed diverse metabolic interactions and was associated with 55 metabolite measures, largely comprising of lipoprotein subclasses and fatty acids, as well as CRP (Fig. 3; Additional file 8: Table S7). Cholesterol Nath et al. Genome Biology (2017) 18:146 Page 6 of 15 esters in small HDL particles were most strongly associated with the PM (meta-Pvalue = 9.45 × 10 –20 ). HDL has been shown to exhibit antithrombotic properties by modulating platelet activation and aggregation and the coagulation pathway [37]. Also, various LDL subclasses of lipoproteins were associated with the PM, which is consistent with our understanding that LDL influences platelet activity. For example, LDL has been shown to influence platelet activity either by enhancing platelet responsiveness to aggregating stimuli or by inducing aggregation [38, 39]. Moreover, LDL-specific binding sites on platelets have also been reported [40, 41]. As noted above, the PM was associated with platelet counts, and adjustment for platelet counts in the YFS resulted in attenuation of approximately half of the weakest associations between PM and metabolite measures; however, the strongest were maintained (Additional file 7: Table S6). Association with VLDL particle size and three others were gained following the adjustment (Additional file 7: Table S6). Neutrophil module Three loci were identified as mQTLs for the neutrophil module (NM; Fig. 2a and Table 2). The top mQTL was intronic to LRRC8A at 9q34.11 (rs13297295; meta-P value = 3.93 × 10 –11 ; Additional file 4: Figure S5a). LRRC8A encodes a trans-membrane protein shown to play a role in Band T-cell development and T cell function [42, 43]. Two additional intergenic mQTLs were located at the TAGAP locus at 6q25.3 (rs2485364; meta-P value = 3.93 × 10 –9 ) and at 20q12 (rs140929198; meta-P value = 1.41 × 10 –9 ) (Additional file 4: Figure S5b, c). Rs13297295 was a strong trans regulator of NM and was an eQTL for eight NM genes (ten unique probes), in particular the major alpha defensins (DEFA1-DEFA4), the genes of highest centrality in the module (Fig. 2d; Additional file 9: Table S8). Rs13297295 was a cis-eQTL for another core NM gene, LCN2 (permuted meta-P value = 1 × 10 –4 ) (Fig. 2d; Additional file 9: Table S8). LCN2 is expressed in neutrophils and inducible by TLR activation, acting as an antimicrobial agent via sequestration of bacterial siderophores to prevent iron uptake [44–46]. LCN2’s role in acute phase response appears to be related to cardiovascular diseases, such as heart failure [47]. At the TAGAP locus, rs2485364 was a transeQTL for eight NM genes (ten probes) and was also a strong driver of LCN2 (meta-Pvalue = 9.11 × 10 –17 ) (Fig. 2d and Additional file 9: Table S8). Consistent with our findings, neutrophils from LCN2-deficient mice have been shown to have impaired chemotaxis and phagocytic capability and increased susceptibility to bacterial and yeast infections compared to wild type [48, 49]. This suggests a possible functional role of TAGAP variants in regulating neutrophil migration through LCN2. NM was associated with 121 circulating metabolite measures (~76% of all metabolite measures analyzed) as well as CRP (Fig. 3; Additional file 8: Table S7). The strongest is the previously reported association with inflammatory biomarker GlycA (meta-Pvalue = 2.68 × 10 – 25 ) [10]; however, NM’s association with various lipoprotein subclasses, particle sizes of lipoproteins, fatty acids, cholesterol, apolipoproteins, glycerides and phospholipids, amino acids, and other small molecules indicates it has a potentially major role in linking neutrophil function to metabolism. Lipid-leukocyte module Together with NM, the lipid-leukocyte module (LLM) showed extensive metabolic associations. Overall, 123 metabolite measures and CRP were associated with LLM, with the strongest being the ratio of triglycerides to phosphoglycerides (meta-Pvalue = 5.16 × 10 –138 ; Fig. 3; Additional file 8: Table S7). With the inclusion of the YFS, these findings strongly replicate previous associations between LLM and metabolite measures [14] as well as detecting additional associations. We also confirm the previous strong negative association between CRP and LLM (meta-Pvalue = 8.16 × 10 –20 ). Consistent with previous studies, no mQTLs were detected for LLM. General immune modules A and B No mQTLs were associated with general immune modules A and B (GIMA and GIMB); however, these modules were associated with 97 and 82 metabolite measures, respectively (Fig. 3; Additional file 8: Table S7). Cholesterol esters in small HDL and the mean diameter for VLDL particles exhibited the strongest associations with GIMA (meta-Pvalue = 1.56 × 10 –30 ) and GIMB (meta-Pvalue = 1.83 × 10 –15 ), respectively. The GIMA was also associated with omega-3 fatty acid levels (meta-Pvalue=4.1×10 –8 ) and CRP (meta-Pvalue = 5.7 × 10 –5 ) while GIMB was not, perhaps due to the subtly difference pathway enrichments for each module (Table 1). Other metabolite measures associated with these two modules include mainly the VLDL and HDL subclass of lipoproteins and fatty acids; due to their large size and heterogeneous composition, however, interpretation of metabolic relationships of GIMA and GIMB is limited. Long-term stability of interactions between metabolite measures, immune gene modules, and mQTLs The 216 individuals in both the DILGOM 2007 and 2014 follow-up allowed investigation of the long-term stability of immunometabolic and mQTL relationships. Across this seven-year period, the eight immune gene coexpression networks were strongly preserved (all preservation statistics’permutation Pvalues <0.001; Additional file 10: Table S9). The metabolite– Nath et al. Genome Biology (2017) 18:146 Page 7 of 15 metabolite correlation structure was also largely consistent between DIGOM07 and DILGOM14 (Additional file 4: Figure S6). Next, we examined how metabolite interactions with immune gene modules changed over the 7-year time period (“Methods”). The LLM–metabolite measure associations were the most consistent over time with 90 and 79 metabolite measures reaching significance in DILGOM07 and DILGOM14, respectively, of which 74 were significant at both time points (Fig. 4a; Additional file 11: Table S10). The direction and effect size of LLM– metabolite measure associations were largely maintained (Fig. 4b). For the neutrophil module, the pyruvate association was significantly maintained over time; however, there was some evidence that other expected associations with NM were stable over time, including GlycA (Additional file 11: Table S10). While no associations with metabolite measures were significantly maintained for the platelet module, rs1354034 was a temporally stable mQTL of PM (mQTL Pvalue = 4.87 × 10 –7 ). No other mQTLs reached significance for temporal stability. While we were powered to topologically replicate immune modules between time-points, power to detect module–metabolite associations and mQTLs was still limited, with only the strongest associations reaching significance. For the latter, the effect sizes for module associations were generally consistent between the time points (Additional file 4: Figure S7), with the exception of GIMB. With the particularly strong consistency of associations for the LLM and NM, it may be that smaller modules, which capture more defined transcriptional programs, are the most temporally stable in terms of their phenotype associations. However, given the robustness of these associations between independent cohorts, we anticipate that, as long-term omics follow-up of population-based cohorts increase in sample size, more of these discovered associations will become statistically significant over time. Conclusions This study has utilized over 2000 individuals to map the immuno-metabolic crosstalk operating in circulation. We have identified and characterized eight robust immune gene modules, their genetic control, and interactions with diverse metabolite measures, including many of clinical significance (e.g., triglycerides, HDL, LDL, branched-chain amino acids). Also, several significant metabolite measures identified here, particularly branched chain amino acids and fatty acids, have been previously shown to be predictive of cardiovascular events and the development of T2D [6, 50]. Furthermore, our findings are consistent with and build upon those of previous studies. In addition to five newly identified gene modules, their mQTLs and metabolite interactions, we have replicated the previously characterized LL module and confirm its association with lipoprotein IDL_P L_LDL_P VLDL_D XL_HDL_P XXL_VLDL_P ACACE ACE ALA ALB APOA1 APOB APOB/APOA1 IDL_L L_LDL_L LDL_D XL_HDL_L XXL_VLDL_L BOHBUT HDL_D IDL_PL L_LDL_PL XL_HDL_PL XXL_VLDL_PL CIT CREA CRP IDL_C L_LDL_C XL_HDL_C XXL_VLDL_C DHA DHA/FA IDL_CE L_LDL_CE XL_HDL_CE XXL_VLDL_CE EST_C IDL_FC L_LDL_FC XL_HDL_FC XXL_VLDL_FC FAW3 FAW3/FA FAW6 FREE_C IDL_TG L_LDL_TG XL_HDL_TG XXL_VLDL_TG GLC GLN GLOL GLY GP L_HDL_P M_LDL_P XL_VLDL_P HDL_C HDL_TG HDL2_C HDL3_C HIS L_HDL_L M_LDL_L XL_VLDL_L ILE L_HDL_PL M_LDL_PL L_HDL_C M_LDL_C L_HDL_CE M_LDL_CE XL_VLDL_CE LA LA/FA LAC LDL_C LDL_TG LEU L_HDL_FC M_LDL_FC XL_VLDL_FC MUFA MUFA/FA L_HDL_TG M_LDL_TG XL_VLDL_TG L_VLDL_P M_HDL_P S_LDL_P L_VLDL_L M_HDL_L S_LDL_L PC PHE PUFA PUFA/FA PYR L_VLDL_PL M_HDL_PL S_LDL_PL L_VLDL_C M_HDL_C S_LDL_C REMNANT_C L_VLDL_CE M_HDL_CE S_LDL_CE SERUM_C SERUM_TG SFA SFA/FA SM L_VLDL_FC M_HDL_FC S_LDL_FC TG_PG TOT_CHO TOT_FA TOT_PG TYR L_VLDL_TG M_HDL_TG S_LDL_TG UNSAT M_VLDL_P S_HDL_P VAL VLDL_C VLDL_TG M_VLDL_L S_HDL_L M_VLDL_PL S_HDL_PL M_VLDL_C S_HDL_C M_VLDL_CE S_HDL_CE M_VLDL_FC S_HDL_FC M_VLDL_TG S_HDL_TG S_VLDL_P S_VLDL_L S_VLDL_PL S_VLDL_C S_VLDL_CE S_VLDL_FC S_VLDL_TG XS_VLDL_P XS_VLDL_L XS_VLDL_PL XS_VLDL_C XS_VLDL_CE XS_VLDL_FC XS_VLDL_TG Apolipoproteins Cholesterol Small Metabolites Amino Acids Inflammatory Markers Lipoprotein particle size Glycerides & Phospholipids Fatty Acids HDL IDL LDL VLDL FDR adjusted P-value < 1.25 x 10 -04 1.25 x 10 -03 6.25 x 10 -03 > 6.25x10 -03 Significance threshold 1.25 x 10 -04 Lipoprotein Subclasses Lipid Leukocyte Module (LLM) DILGOM14 (N=79) DILGOM07 (N=90) XL_VLDL_PL C_LDL V _LX FAW6/FA −0.50 −0.25 0 0.25 0.50 −0.50 −0.25 0 0.25 0.50 PUFA/FA FAW6/FA L_HDL_FC L_HDL_FC L_HDL_CE XL_HDL_PL UnSat XL_HDL_FC XL_HDL_CE XL_HDL_C HDL2_C L_HDL_PL XL_HDL_P TG_PG S_VLDL_TG M_HDL_TG VLDL_TG M_VLDL_TG L_VLDL_PL L_VLDL_TG L_VLDL_C HDL_TG MUFA S_LDL_TG SFA XXL_VLDL_PL ApoB Remnant_C XS_VLDL_P S_VLDL_C M_VLDL_CE Beta Estimates DILGOM07 Beta Estimates DILGOM14 Both DILGOM07 and DILGOM14 DILGOM07 only DILGOM14 only Not significant B A Significant metabolite measure associations with LLM LA/FA XS_VLDL_L L_VLDL_CE Fig. 4 Temporally stable metabolite measure associations with the LLM. aCircular heatmap for association between each metabolite measure and the LLM. bComparison of the effect size estimates of metabolite measure association with LLM in DILGOM07 and DILGOM14 shows that the overall association patterns are consistent across the two time-points. Colors denote metabolites that are significantly associated with the LLM in DILGOM07 only (orange), DILGOM14 only (blue), and across both time-points (green). The grey dashed line is the x = y line Nath et al. Genome Biology (2017) 18:146 Page 8 of 15 subclasses, lipids, fatty acids, and amino acids [13, 14]. Associations between the core genes in the LL module and isoleucine, leucine, and various lipids were also identified independently in the KORA cohort [12]. Importantly, we have shown the long-term stability of LL and neutrophil module coexpression and interactions with metabolite measures, and we have greatly expanded the number of known biomarkers associated with the NM from one (GlycA) to 123 [10]. Our study has also expanded the widespread trans eQTL effects at the ARHGEF3 locus [23], shows them to be strongly maintained within individuals over time, and further identifies extensive interactions with lipoprotein measures that may be a consequence of these trans effects. Taken together, our analyses illustrate the rapidly growing body of evidence intimately linking the immunoinflammatory response to the blood metabolome. With finer-resolution maps of these interactions, new biomarkers of chronic and acute inflammatory states are likely to emerge. With in vivo and interventional studies, modulation of these metabolite–immune interactions through existing lipid-lowering medications, gut microbe effects, or dietary changes may provide new ways the immune system itself can be utilized to lessen the burden of cardiometabolic disease. Methods Study populations This study used data from two population-based cohorts, the Dietary, Lifestyle, and Genetic determinant of Obesity and Metabolic syndrome (DILGOM; N = 518) and the Cardiovascular Risk in Young Finns Study (YFS; N = 1650), which have been described in detail elsewhere [13, 51]. All subjects enrolled in these studies gave written informed consent. The DILGOM study is a subsample of the FINRISK 2007 cross-sectional population-based survey, which recruited a random sample of 10,000 individuals between 25 and 74 years of age, stratified by sex and 10-year age groups, from five study areas in Finland. All 6258 individuals who participated in the FINRISK 2007 baseline health examination were invited to attend the DILGOM study (N = 5024), 630 of whom underwent at least one of the genotyping, transcriptomics, or metabolomics profiling considered here. In 2014, a follow-up study was conducted, for which 3735 individuals from the original study re-participated. Samples collected in 2007 and 2014 are referred to as DILGOM07 and DILGOM14, respectively. The YFS is a longitudinal prospective cohort study that started in 1980, with follow-up studies carried out every 3 years, to monitor cardiovascular disease risk factors in children and adolescents from five major regions of Finland (Helsinki, Kuopio, Turku, Oulu, and Tampere). A total of 3596 children and adolescents in age groups 3, 6, 9, 12, 15, and 18 years participated in the baseline study; these children were randomly selected from the national public register and their details are described in [51]. In this current study, data collected from the 2011 follow-up study (participants aged 34, 37, 40, 43, 46, and 49 years) were analyzed. Sample collection Venous blood was collected following an overnight fast in all three studies. Samples were centrifuged and the resulting plasma and serum samples were aliquoted into separate tubes and stored at −70 °C for analyses. Protocols for the blood sampling, physiological measurements, and clinical survey questions were similar across the YFS and DILGOM studies and are described extensively in [13, 52]. Genotyping and imputation Whole blood genomic DNA obtained from both cohorts was genotyped using the Illumina 610-Quad SNP array for DILGOM07 (N = 555) [13] and a custom generated 670 K Illumina BeadChip array for YFS (N = 2443) [53]. The 670 K array shares 562,643 SNPs with the 610-quad array. The 670 K array removes poorly performing SNPs from the 610-quad array and improves copy number variation coverage [53]. Genotype calling was performed with the Illuminus clustering algorithm [54]. Quality control was as previously described in [13] and [53] for DILGOM and YFS, respectively. Genotypes were imputed to the 1000 Genomes Phase 1 version 3 reference panel using IMPUTE2 in both DILGOM and YFS [17]. Poorly imputed SNPs based on low call-rate (<0.90 for DILGOM, <0.95 for YFS), low-information score (<0.4), minor allele frequency <1%, and deviation from Hardy– Weinberg equilibrium (P<5×10 –6 ) were then removed. A total of 7,263,701 SNPs in DILGOM and 6,721,082 in YFS passed quality control, with 6,485,973 common between the two. A total of N = 518 samples in DILGOM and N = 2443 samples in YFS individuals passed quality control filters. Metabolomics profiling Metabolite concentrations for DILGOM07 (N = 4816), DILGOM14 (N = 1273), and YFS (N = 2046) were quantified from serum samples utilizing a high-throughput NMR metabolomics platform (Brainshake Ltd, Helsinki, Finland) [18, 55]. Details of the experimental protocol, including sample preparation, NMR spectroscopy, and metabolite identification, have been previously described in [13, 18]. A total of 159 metabolite measures were assessed, of which 148 were directly measured and 11 were derived (Additional file 1: Table S1). The 148 measures include the constituents of 14 lipoprotein subclasses (98 Nath et al. Genome Biology (2017) 18:146 Page 9 of 15