Research Article

Common Genetic Variants Modulate Pathogen-Sensing Responses in Human Dendritic Cells

See allHide authors and affiliations

Science  07 Mar 2014:
Vol. 343, Issue 6175, 1246980
DOI: 10.1126/science.1246980

Structured Abstract

Introduction

Variation in an individual’s response to environmental factors is likely to influence susceptibility to complex human diseases. The genetic basis of such variation is poorly understood. Here, we identify natural genetic variants that underlie variation in the host innate immune response to infection and analyze the mechanisms by which such variants alter these responses.

Embedded Image

Identifying the genetic basis of variability in the host response to pathogens. A cohort of 534 individuals donated blood for (a) genotyping of common DNA variants and (b) isolation of immune DCs. DCs were stimulated with viral and bacterial components, and the variability in individuals’ gene expression responses was mapped to specific DNA variants, which were then shown to affect binding of particular transcription factors.

Methods

We derived dendritic cells (DCs) from peripheral blood monocytes of healthy individuals (295 Caucasians, 122 African Americans, 117 East Asians) and stimulated them with Escherichia coli lipopolysaccharide (LPS), influenza virus, or the cytokine interferon-β (IFN-β) to generate 1598 transcriptional profiles. We genotyped each of these individuals at sites of common genetic variation and identified the genetic variants that best explain variation in gene expression and gene induction between individuals. We then tested mechanistic predictions from these associations using synthetic promoter constructs and genome engineering.

Results

We identified 264 loci containing genetic variants associated with variation in absolute gene expression in human DCs, of which 121 loci were associated with variation in the induction of gene expression by one or more stimuli. Fine-mapping identified candidate causal single-nucleotide polymorphisms (SNPs) associated with expression variance, and deeper functional experiments localized three of these SNPs to the binding sites of stimulus-activated transcription factors. We also identified a cis variant in the transcription factor, IRF7, associated in trans with the induction of a module of antiviral genes in response to influenza infection. Of the identified genetic variants, 35 were also associated with autoimmune or infectious disease loci found by genome-wide association studies.

Discussion

The genetic variants we uncover and the molecular basis for their action provide mechanistic explanations and principles for how the innate immune response to pathogens and cytokines varies across individuals. Our results also link disease-associated variants to specific immune pathways in DCs, which provides greater insight into mechanisms underlying complex human phenotypes. Extending our approach to many immune cell types and pathways will provide a global map linking human genetic variants to specific immunological processes.

Immune Variation

It is difficult to determine the mechanistic consequences of context-dependent genetic variants, some of which may be related to disease (see the Perspective by Gregersen). Two studies now report on the effects of stimulating immunological monocytes and dendritic cells with proteins that can elicit a response to bacterial or viral infection and assess the functional links between genetic variants and profiles of gene expression. M. N. Lee et al. (10.1126/science.1246980) analyzed the expression of more than 400 genes, in dendritic cells from 534 healthy subjects, which revealed how expression quantitative trait loci (eQTLs) affect gene expression within the interferon-β and the Toll-like receptor 3 and 4 pathways. Fairfax et al. (10.1126/science.1246949) performed a genome-wide analysis to show that many eQTLs affected monocyte gene expression in a stimulus- or time-specific manner.

Abstract

Little is known about how human genetic variation affects the responses to environmental stimuli in the context of complex diseases. Experimental and computational approaches were applied to determine the effects of genetic variation on the induction of pathogen-responsive genes in human dendritic cells. We identified 121 common genetic variants associated in cis with variation in expression responses to Escherichia coli lipopolysaccharide, influenza, or interferon-β (IFN-β). We localized and validated causal variants to binding sites of pathogen-activated STAT (signal transducer and activator of transcription) and IRF (IFN-regulatory factor) transcription factors. We also identified a common variant in IRF7 that is associated in trans with type I IFN induction in response to influenza infection. Our results reveal common alleles that explain interindividual variation in pathogen sensing and provide functional annotation for genetic variants that alter susceptibility to inflammatory diseases.

Susceptibility to complex diseases depends on both genetic predisposition and exposure to environmental factors, with interactions between the two (G × E interactions) likely contributing substantially to disease risk (1, 2). However, the extent and mechanisms by which common human genetic variants interact with the environment remain poorly explored and have been difficult to detect in clinical studies (1, 3). Genetic analysis of molecular traits, such as gene expression profiles, offers a promising way to dissect the molecular mechanisms underlying G × E interactions. Expression quantitative trait locus (eQTL) studies have been used to map genetic variants contributing to variation in gene expression, but have largely focused on steady-state expression in humans (4, 5)and thus have excluded G × E interactions. In model organisms, differences in growth conditions or treatment with various stimuli have revealed the existence of response eQTLs (reQTLs) (69), defined as QTLs associated with the change in expression after stimulation. Here, we sought to identify reQTLs in humans, to explain the mechanism by which the environment interacts with these variants, and to determine whether these variants are associated with immune-mediated diseases.

We used dendritic cells (DCs) of the innate immune system as a model system for reQTL studies, with physiological and clinical relevance. DCs play a direct role in the host recognition of pathogens using specialized sensors that engage well-characterized signaling and transcriptional networks. For example, bacterial lipopolysaccharide (LPS) activates two distinct arms of the Toll-like receptor 4 (TLR4) pathway, whereas influenza infection primarily activates the RNA-sensing TLRs (for example, TLR3) and the RIG-I–like receptors (for example, RIG-I) (10). These, in turn, lead to the translocation of transcription factors (TFs) from the cytoplasm into the nucleus to induce the expression of immune genes, including interferon-β (IFN-β) secretion that engages the type I IFN response pathway to induce the expression of hundreds of antiviral effectors. Genetic studies have associated common variants near many genes in these pathways with risk of different inflammatory diseases (11, 12). DCs also play a critical role in the pathologic immune responses underlying inflammatory diseases (1113), also reflected in recent genome-wide association studies (GWAS) of several diseases (1417), especially inflammatory bowel disease (14).

Results

Assessing the Impact of Genetic Variation on Pathogen Sensing in Primary Human DCs

We developed an integrated experimental and computational pipeline (Fig. 1A) to identify variability in human DC responses and associate this variability with common genetic variants. First, we optimized a high-throughput protocol to isolate primary CD14+CD16lo monocytes from human blood samples (fig. S1, A to F), differentiate them into monocyte-derived DCs (MoDCs), and stimulate them with three immunostimulatory agents (Fig. 1B): Escherichia coli LPS, influenza virus, or IFN-β (a cytokine induced by LPS and influenza). Second, we collected genome-wide transcript profiles from resting and stimulated DCs from 30 healthy individuals and computationally identified a “signature” of 415 genes that would be informative of variation in the response in a larger cohort. Third, we generated 1598 transcriptional profiles [using an amplification-free platform suitable for small cell numbers (18)] from DCs isolated from 534 healthy individuals (295 Caucasians, 122 African Americans, and 117 East Asians) (table S1) in four states: resting, LPS-stimulated, influenza-infected, and IFN-β–stimulated. Finally, we associated expression variation with genetic variation to identify cis- and trans-eQTLs and reQTLs, and we investigated candidate reQTLs for their mechanism of action.

Fig. 1 A strategy to identify gene-by-environment interactions in the innate immune responses of primary human DCs.

(A) Strategy used to identify baseline and response eQTLs and reQTLs, consisting of five steps: (i) high-throughput isolation and stimulation of primary human MoDCs from 560 healthy individuals (dotted slices, male; solid-colored slices, female) collected as part of the PhenoGenetic cohort; (ii) whole-genome gene expression measurements in a subset of the cohort; (iii) selection of signature gene set, consisting of regulators and regulated genes; (iv) digital multiplex gene expression measurements of signature genes in the entire cohort; and (v) mapping of genetic variation to expression variation. GM-CSF, granulocyte-macrophage colony-stimulating factor; IL-4, interleukin-4. (B) Model of innate immune pathways activated by three stimuli demonstrating their downstream relationships. LPS from E. coli engages the TLR4 receptor; IFN-β engages the heterodimeric IFNAR receptor; influenza A/PR8 (ΔNS1) (“FLU”) engages the cytoplasmic TLR3 and RIG-I receptors. Receptor engagement activates signal transduction cascades that regulate expression of inflammatory genes, IFNs, and IFN-stimulated genes. IFNAR activation also occurs during LPS and FLU stimulations because LPS and FLU both induce IFN production, leading to activation of ISREs. JAK1, Janus kinase 1.

Variability of Pathogen-Sensing Responses Between Individuals

To assess the interindividual variability of DC responses, we first profiled genome-wide expression using microarrays in resting, LPS-stimulated, and influenza-infected MoDCs isolated from each of 30 healthy individuals (18 Caucasians, 6 Asians, and 6 African Americans). We found 1413 genes that were regulated in LPS- or influenza-treated cells [log2(fold change) >0.75 or <−1.5; false discovery rate (FDR) < 0.01; fig. S2A and table S2A; see table S2, B to I, for enriched biological processes] at 5 and 10 hours, respectively (time points selected to have maximal expression of induced clusters) (fig. S1, B and C).

We quantified reproducibility in these responses by recalling 12 of the 30 donors 2 to 9 months after the first collection for MoDC isolations and profiling (table S1); we identified genes whose expression interindividual variance is significantly higher than their intraindividual variance based on the serial replicates [taking into account known covariates including gender, age, and population and unknown factors using surrogate variable analysis (19)]. Two hundred twenty-two of the 1413 regulated genes (16%) showed significantly higher (FDR < 0.1) inter- than intraindividual variability, either in their absolute expression or in the differential expression (stimulated/baseline level) to at least one stimulus (Fig. 2A and table S3). These results suggest that there is consistent variation in these traits that may have a genetic basis.

Fig. 2 Genome-wide expression profiles in MoDCs reveal response phenotypes.

(A) Coefficient of variation (CV) of gene expression between 30 different donors (“Interindividual CV”) plotted against CV of expression within 12 serial replicate samples (“Intraindividual CV”) for each differentially regulated (fold change >0.75 or <–1.5) gene after LPS or FLU stimulation. Yellow (up-regulated) and purple (down-regulated) circles represent genes with significant (moderated t test, FDR < 0.1) inter- versus intraindividual variation. (Right) log2(expression, microarray data) of CLEC4F in baseline, LPS-stimulated, and FLU-infected MoDCs from 30 donors and 12 replicates, demonstrating example of a gene that shows significant (FDR < 0.01) inter- versus intraindividual variation after LPS and FLU stimulations but not at baseline (fig. S2B). Standard error of replicate samples (n = 12) is shown for each sample. (B) Pie chart of 415 signature genes selected for Nanostring code set: 222 (49%) are regulated genes that showed significant (mixed model variance components test, permutation FDR < 0.1) inter- versus intraindividual variability; 61 (14%) are curated, regulated genes with a known function in the innate immune response; 76 (17%) are curated regulators in the TLR4, TLR3, RIG-I, and IFNAR pathways; 41 (9%) are control genes including low-variance genes, sex-specific genes, and nonexpressed genes; 28 (6%) are regulated genes that were reported in the regions of autoimmune and infectious disease GWAS; and 21 (5%) are regulated genes that showed significant inter- versus intrapopulation variability. (C) Gene expression heat map of the 415-gene signature in MoDCs from the microarray study (30 individuals) and the Nanostring study (534 individuals). Each row represents a gene; each column represents a donor sample at baseline, stimulated with LPS, infected with FLU, or stimulated with IFN-β. Rows were clustered by k-means clustering of Nanostring data set with major clusters (I, II, IIIa, IIIb, and IV) labeled. Between the two heat maps, each row was labeled with colored dashes corresponding to one of the six categories described in (B).

Expression Profiling of a Pathogen Response Signature

Mapping the genetic basis of interindividual variation in pathogen responses requires profiling of DC gene expression in a larger cohort; this poses a substantial technical challenge given the limited numbers of primary cells and multiple stimuli. Furthermore, the responses of DCs to virus, bacterial ligands, and IFN are limited in scope and do not encompass the entire genome. We therefore defined a 415-gene signature set (Fig. 2B and table S3) that could be monitored in small numbers of cells using a sensitive multiplex RNA detection system (18), allowing us to scale up our study. The signature consisted of (i) all 222 of 1413 regulated genes with greater inter- than intraindividual variability in the microarray study; (ii) all 24 regulated genes exhibiting greater inter- than intrapopulation variability (FDR < 0.1) (table S3) in the microarray study; (iii) 76 genes comprising the known components of the TLR4, TLR3, RIG-I, and IFNAR pathways (Fig. 1B) (10, 2022); (iv) 61 regulated genes that play key roles in the DC response (for example, IFNB1 and IFITM3) (2224); (v) 28 regulated genes residing in loci previously associated with autoimmune or infectious diseases (25); and (vi) 35 control genes including those that had among the lowest interindividual variance in the microarray profiles. This representative signature allowed us to monitor genes with high interindividual variability, as well as key components and responses of the pathogen-sensing pathway.

We then generated 1598 transcriptional profiles, using the 415-gene signature, from MoDCs isolated from 534 healthy individuals (and 37 serially collected replicates) in up to four conditions: resting (528 individuals), LPS-stimulated (356), influenza-infected (342), and IFN-β–stimulated (284) (Fig. 2C; fig. S1, E and F; and table S1). The IFN-β stimulus allowed us to partition genes that were induced by both LPS and influenza (cluster III) into those induced secondary to type I IFN signaling (cluster IIIa) or by other mechanisms shared between these bacterial and viral sensing pathways (cluster IIIb), such as nuclear factor κB (NF-κB) or activating protein 1 (AP-1) activation (Figs. 1B and 2C). The signature expression profiles clustered similarly to those from the microarray analysis (Fig. 2C and fig. S2C) and captured most of the variance in the genome-wide profile [cross-validation, >0.99 ground-truth Pearson correlation coefficient; fig. S2D (26)], demonstrating the utility of the signature to capture the genome-wide response.

reQTLs That Modulate Cellular Responses to Pathogens

We mapped cis-eQTLs and cis-reQTLs by testing for association between common single-nucleotide polymorphisms (SNPs) (minor allele frequency >5%; genotyped with Illumina HumanOmniExpress BeadChip) and variation in either absolute transcript abundance (eQTL) or stimulation-induced change in transcript abundance (reQTL) in a nearby gene (in which the transcriptional start site or stop codon is within 1 Mb of the SNP). We pooled individuals from the three human populations together and included covariates and principal components to increase statistical power while avoiding systematic confounders such as population structure and expression heterogeneity. We identified 264 genes with cis-eQTLs in at least one condition (permutation FDR < 0.05); (Fig. 3A; fig. S3, A to C; and table S4, A to D) (27). Notably, 22 of 264 genes also associated with additional independent cis-eQTLs after conditioning on the top five most significant SNPs in each cis-eQTL region (table S6 and fig. S3D) (27), which reflects the complexity of the regulatory landscape around each gene.

We detected 121 cis-reQTLs (Fig. 3B and table S4, E to G; permutation FDR < 0.05; 91% internal reproducibility in at least two human populations, table S5), the subset of genetic variants that affect the induction of gene expression by LPS, influenza, or IFN-β. Of the 121 genes with cis-reQTLs, 7 associations were found only in the influenza condition (for example, IFNA21; Fig. 3C and fig. S3B; meta-analysis, P < 1 × 10−5; permutation FDR = 0.02 to 0.03); 15 in both the LPS and influenza conditions (for example, TEC; Fig. 3C and fig. S3B); and 57 in all three stimulation conditions (for example, ARL5B, SLFN5, and CLEC4F; Fig. 3, C to E, and fig. S3B), likely reflecting activation of the shared IFN-β pathway. We hypothesized that causal variants driving cis-reQTLs alter the sequence of genomic elements that respond to TFs downstream of the pathogen-sensing receptors (Fig. 3F), and represent gene-by-environment interactions.

Fig. 3 Association analysis reveals cis-eQTLs and cis-reQTLs.

(A and B) Manhattan plot of cis-eQTLs (A) (baseline expression) and cis-reQTLs (B) (LPS-, FLU-, and IFN-β–stimulated fold changes relative to baseline) showing −log10(P values) (left y axis) and R2 values (right y axis) for all cis-SNPs, displayed on the x axis with associated genes ordered by chromosomal location. (C) Box-whisker plots showing expression [left; log2(nCounts), y axis] or fold change [right; log2(fold), y axis] of DCBLD1, IFNA21, TEC, and ARL5B in resting, LPS-stimulated, FLU-infected, and IFN-β–stimulated MoDCs as a function of genotype of the respective cis-SNPs (x axis: rs27434, rs10964871, rs10938526, and rs11015435). African Americans, Asians, and Europeans in this order are displayed as separate box-whisker plots adjacent to each other in each condition. −Log10(P values) and β statistics are displayed in top right corners. (D and E) Allelic imbalance analysis of SLFN5 (D) and CLEC4F (E) in resting, LPS-stimulated, FLU-infected, and IFN-β–stimulated MoDCs, showing the ratio of gene expression between the major and minor alleles in heterozygote (rs11651240 for SLFN5, rs2075221 for CLEC4F) cDNA samples [n = 8 in (D); n = 9 in (E)] normalized to the ratio in the corresponding genomic DNA samples; significant deviation from 1.0 (dashed line) is consistent with allelic imbalance. Data are from one experiment representative of three (mean and SD shown). *P < 0.01, **P < 0.001, compared to unstimulated cells (Student’s t test). On the right panels, box-whisker plots showing expression [left; log2(nCounts), y axis] or fold change [right; log2(fold), y axis] of SLFN5 (D) and CLEC4F (E) in resting, LPS-stimulated, FLU-infected, and IFN-β–stimulated MoDCs as a function of genotype of the respective cis-SNPs: rs11867191 and rs2075221. (F) Schematic showing the different combinations of stimuli leading to significant cis-reQTLs, with the most significant examples listed. Specificity to conditions was defined with M value >0.9 taken as the inclusion criteria and M value <0.1 taken as the exclusion criteria for each condition.

To validate cis-reQTLs, we quantified allele-specific expression for pan-stimulation–specific associations in heterozygote individuals (Fig. 3, D and E). This was feasible for a gene with an exonic SNP (CLEC4F rs2075221) that was the most significant reQTL SNP, and for a gene with an exonic SNP (SLFN5 rs11651240) in linkage disequilibrium with the most significant reQTL SNP (rs11867191, R2 = 0.501 CEPH, CEU). As predicted, transcripts derived from the major and minor alleles differed in expression after stimulation with LPS (SLFN5: P < 0.001, t test; CLEC4F: P < 0.01); influenza (SLFN5: P < 0.001; CLEC4F: P < 0.01); or IFN-β (SLFN5: P < 0.001; CLEC4F: P < 0.001) but not at baseline (Fig. 3, D and E).

Cis-reQTLs That Alter the Sequence of TF Binding Sites

We hypothesized that the causal variants underlying cis-reQTLs functionally alter the chromosomal binding sites of TFs activated by one or more stimuli. To fine map reQTLs, we performed a trans-ethnic meta-analysis (28) of imputed variants in each population (~10 M). We then examined whether the most highly associated meta-reQTLs overlap with TF binding sites identified in high-throughput human chromatin immunoprecipitation with massively parallel DNA sequencing (ChIP-seq) data sets (for example, ENCODE) or computationally predicted conserved regulatory elements (2931). We found substantial enrichment (table S7) of known binding sites for TFs from the signal transducer and activator of transcription (STAT) family (TFs that are activated downstream of type I IFN signaling). The most significant enrichment over background was found for STAT2 (116-fold, binomial P < 2.55 × 10−21) and STAT1 (126-fold, binomial P < 2.98 × 10−13) binding sites (derived from ChIP-seq in IFN-α–stimulated K562 cells) within cis-eQTLs after IFN-β stimulation.

Among the 57 genes (for example, SLFN5, CLEC4F, and ARL5B; Fig. 3F) with cis-reQTLs in all three stimuli, we observed that the most significant cis-reQTL in the (i) SLFN5 locus (rs11080327) lies in an ENCODE ChIP-seq signal (29) for STAT1 (Fig. 4, A to C); (ii) CLEC4F locus (rs35856355) alters a commonly occurring cytosine in a canonical IFN-stimulated response element (ISRE) that is a target of IFN-activated TFs (Fig. 4, B and C) (22, 3032); and (iii) ARL5B locus (rs2130531) changes a guanine in the canonical ISRE motif (Fig. 4, B and C). Additional cis-reQTLs that alter putative STAT binding sites based on ChIP-seq data or predicted ISRE motifs include rs10086852 (PTK2B), rs1981760 (NOD2), rs73023464 (C19ORF12), rs1331717 (IFI44), and rs12064196 (IFI44) (Fig. 4A). Because STAT TFs are activated downstream of the type I IFN receptor (IFNAR) and bind to ISRE motifs (Fig. 4A) (22), we hypothesized that the SNPs in these sites are likely causal SNPs within their respective cis-reQTL regions.

Fig. 4 Functional fine-mapping and mechanism of cis-reQTLs.

(A) Pathway diagram of signal transduction cascade downstream of IFNAR activation. Activation of receptor leads to downstream activation of JAK-STAT cascade, leading to posttranslational activation of STAT and IRF TFs. (B) LocusZoom plots showing the −log10(P values) of imputed cis-eQTLs (y axis) in the chromosomal regions (x axis) of SLFN5, CLEC4F, and ARL5B. The most significant imputed SNPs in each locus are labeled. (C) Schematic representation of alleles in the regions near the SLFN5, CLEC4F, and ARL5B genes that are in STAT2 ChIP-seq binding sites or that perturb ISRE motifs (SNPs are shown as vertical bars and in red letters). (D) EMSAs with 24- to 26-bp radiolabeled dsDNA probes—containing a known ISRE motif control, a mutated ISRE motif control, the CLEC4F rs35856355 major (C) or minor allele (A) sequences, or the ARL5B rs2130531 major (G) or minor allele (A) sequences—incubated with nuclear lysates from IFN-β–stimulated MoDCs. On the right, supershift assays with or without antibodies against IRF1, IRF9, and STAT2 (designated α-IRF-1, and so on) with the CLEC4F rs35856355 major (C) probe are shown. (E) Luciferase expression from reporter constructs transfected into HEK293 cells that were left unstimulated or were stimulated with IFN-β (1000 U/ml) for 21 hours. Sequences (150 to 200 bp) from the major and minor haplotypes of the SLFN5, CLEC4F, and ARL5B regions were subcloned 5′ of a minimal promoter and firefly luciferase gene. Firefly luciferase expression was normalized to Renilla luciferase expression expressed from cotransfected plasmid. (F) Fold change log2(IFN-β stim/unstim) of signature genes in wild-type (wt) HEK293 cells (rs11080327A/G), plotted against fold change in CRISPR-converted (rs11080327G/G) HEK293 cells. Data are from one experiment representative of three [mean and SD shown in (E)]. *P < 0.05, **P < 0.01, compared to unstimulated cells (Student’s t test).

Cis-reQTLs That Affect Differential Binding of Stimulus-Activated TFs

To experimentally validate our predicted regulatory mechanisms, we determined whether IFN-β–activated TFs bind differentially at these SNPs. Using radiolabeled 24- to 26-bp (base pair) double-stranded DNA (dsDNA) probes in electrophoretic mobility shift assays (EMSAs), we found that probes encompassing the major alleles (rs35856355C and rs2130531G, which are associated with increased expression in CLEC4F and ARL5B), but not the minor alleles, shifted after incubation with IFN-β–stimulated MoDC nuclear lysates (Fig. 4D). Consistently, nonradiolabeled major, but not minor, allele probes competed for binding to the IFN-β–stimulated nuclear factors (fig. S4A). By incubating with antibodies to TFs known to bind the ISRE element, we found that IRF1 (IFN-regulatory factor 1) supershifted the low–molecular weight band and STAT2 and IRF9 supershifted the higher–molecular weight band, consistent with the known complex of these latter two proteins (32) (Fig. 4D). These results suggested that the SNPs in CLEC4F and ARL5B alter the binding of IFN-activated IRF1, STAT2, and IRF9 TFs.

To directly test the functional effects of the SNPs, we cloned 150- to 200-bp genomic regions surrounding rs11080327 (SLFN5), rs35856355 (CLEC4F), and rs2130531 (ARL5B) upstream of a luciferase reporter driven by a minimal promoter (Fig. 4E). For each region, we created two constructs that differ only at the respective SNPs. Because almost all cell types respond to IFN-β stimulation, we transfected the reporter constructs into human embryonic kidney 293 (HEK293) cells and stimulated the cells with IFN-β. Consistent with our hypothesis, the constructs containing the major alleles of rs11080327 (SLFN5), rs35856355 (CLEC4F), and rs2130531 (ARL5B) induced 11.6-, 2.1-, and 1.5-fold, respectively, more luciferase production than the cognate constructs containing the minor alleles (t test, P < 0.001, P < 0.001, and P < 0.01, respectively) (Fig. 4E). Furthermore, mutation of motif-conserving nucleotides across the CLEC4F ISRE motif reduced induction by IFN-β, whereas mutation of nonconserved nucleotides did not (fig. S4B), consistent with the ISRE consensus motif.

Finally, to further test (33) the functionality of rs11080327—the SLFN5 variant with the strongest functional effect—in its native chromosome, we directly edited the genome using the CRISPR system (34, 35), converting the heterozygous rs11080327A/G in HEK293 cells to a homozygous rs11080327G/G (Fig. 4F and fig. S4C). We then stimulated both native and CRISPR-converted cells with IFN-β. The fold induction of SLFN5 decreased from 3.64 in heterozygous cells to 1.03 in the converted rs11080327 (G) homozygotes, whereas the induction of other genes was unaffected (Fig. 4F and fig. S4C, independent CRISPR construct).

The enrichment, EMSA, reporter assays, and directed mutagenesis results suggest that we have identified causal variants in some of these cis-reQTL regions and demonstrate that the reQTLs occur because of differential binding of stimulus-activated TFs.

A Trans-eQTL That Also Associates with IRF7 as a Cis-eQTL

We next attempted to detect cases where variation in gene expression is explained by distant genetic variants acting in trans. We noted that cis-eQTLs for regulator genes tend to have smaller effect sizes and less significant P values than those for regulated genes (fig. S3E). To decrease the multiple testing burden of detecting trans-associations, we restricted testing SNPs local to genes on our gene signature set. A number of these are cis-eQTLs associated with the expression of genes encoding regulators (for example, TFs) in pathways that we stimulated (Fig. 5A). Of these, rs12805435 was the most significant, associating with the expression of the master antiviral TF IRF7 in cis only after influenza, LPS, or IFN-β stimulation (table S4, A to D). This SNP associated with both the expression and induction of seven additional genes in trans (that is, the transcriptional start site or stop codon is located >1 Mb from the SNP) only after influenza infection: NMI, IFNA4, IFNA10, IFNA13, IFNA21, IFNA17, and IFNA5 (P < 6.26 × 10−5 in expression, P < 1.68 × 10−5 in differential expression) (Fig. 5B and tables S8 to 10).

Fig. 5 Trans-reQTL association at the IRF7 cis-regulatory locus.

(A) Diagram showing selected components of TLR4, TLR3, RIG-I, and IFNAR pathways. Components with significant cis-eQTLs (permutation FDR < 0.05) are shown in black (or red if they also have a trans-eQTL); components that do not have significant cis-eQTLs are gray. (B) Manhattan plot showing the trans-association of rs12805435 to all 415 genes on signature gene set in baseline, LPS-stimulated, FLU-infected, and IFN-β–stimulated conditions. Trans-reQTL to NMI and cluster of IFN-α genes (IFNA4, IFNA10, IFNA13, IFNA17, and IFNA21) are labeled. (C) Expression [log2(nCounts)] of 415 signature genes in FLU-infected MoDCs overexpressing IRF7, plotted against expression in FLU-infected MoDCs overexpressing eGFP control. (D) Expression [log2(nCounts)] of 415 signature genes in HEK293 cells overexpressing IRF7, plotted against expression in HEK293 cells overexpressing eGFP control. Right, expression of genes in cells with eGFP overexpression versus cells without cDNA.

To test whether the trans-associated genes are indeed targets of IRF7, we infected MoDCs that overexpressed IRF7 or a control gene with influenza virus, and found that IRF7 overexpression induced the expression of IFNA2, IFNA13, and IFNA14 in influenza-infected cells (Fig. 5C). In addition, in HEK293 cells, which normally express very low levels of IRF7, we overexpressed IRF7 and observed induction of NMI, IFNA4, IFNA10, IFNA13, IFNA14, and IFNA21 (Fig. 5D), suggesting that IRF7 is sufficient to drive downstream expression. We have thus identified a stimulus-specific cis-eQTL associated with IRF7 expression, which is also a trans-reQTL that underlies the variability of IFN induction in response to influenza infection in humans.

DC reQTLs That Overlap with Autoimmune and Infectious Disease SNPs from GWAS

We examined the subset of loci associated in GWAS with inflammatory disorders. We first extended an analysis of Crohn’s disease in which loci are enriched for genes specifically expressed in DCs (12, 14) to multiple sclerosis, celiac disease, psoriasis, and leprosy (fig. S5A). We found that genes nearest to the susceptibility loci for these diseases were enriched not only in DC-specific genes but also in genes induced by LPS and/or influenza (fig. S5B), which suggests that some of the GWAS loci modulate the expression of the corresponding genes in activated DCs.

Supporting this hypothesis, 15 cis-reQTLs and 23 cis-eQTLs that were not significant at baseline but only became significant after stimulation are the same SNPs previously identified in GWAS of autoimmune and infectious diseases, including Crohn’s disease, multiple sclerosis, celiac disease, psoriasis, and leprosy (25) (Fig. 6, A and B, and table S11). These include NOD2 with leprosy (rs9302752), IRF7 with systemic lupus erythematosus (SLE) (rs4963128), TRAF1 with rheumatoid arthritis (rs881375), and CREM with Crohn’s disease (rs12242110) and ulcerative colitis (rs4246905). For example, rs9302752—a SNP previously associated with susceptibility to leprosy (17)—was associated in our study with the absolute and fold expression of NOD2 under IFN-β stimulation (P = 3.49 × 10−25) but not at baseline (Fig. 6A and table S11); NOD2 plays a known role in pathogen sensing and possibly mycobacterial immunity (36). Similarly, rs4963128—a variant associated with SLE (37)—was associated in our study with the expression of IRF7 after IFN-β stimulation (P = 1.10 × 10−16) but not at baseline (Fig. 6A), in line with the importance of type I IFN responses in SLE pathogenesis (38). We note that rs4963128 is on the same haplotype (R2 = 0.69, D′ = 0.94) as the IRF7 SNP rs12805435 that is associated with the trans-reQTL effect described above. These results suggest a role for innate immune pathogen-sensing pathways in the pathogenesis of these inflammatory disorders.

Fig. 6 Autoimmune and infectious disease–associated SNPs are cis-eQTLs and cis-reQTLs.

(A) Expression [log2(nCounts)] of NOD2 in resting and IFN-β–stimulated MoDCs from 184 Caucasians as a function of genotype of the leprosy GWAS SNP, rs9302752 (left). (Right) expression [log2(nCounts)] of IRF7 in resting and IFN-β–stimulated MoDCs from 184 Caucasians as a function of genotype of the SLE GWAS SNP, rs4963128. (B) Plot showing overlap of genome-wide significant (P < 5 × 10−8) GWAS SNPs with cis-eQTLs and reQTLs in MoDCs, with clinical phenotypes connected to corresponding gene expression phenotypes by lines. Orange circles represent cis-reQTLs (P < 10−7); yellow circles represent stimulus-specific cis-eQTLs (P < 10−7).

Discussion

Although genetic association studies have identified alleles that confer disease risk, little is known about how these genetic variants contribute to disease through their effects on specific biological processes and their interaction with environmental stimuli. We addressed this question by quantifying the DC response to pathogens in a set of genotyped individuals, and then leveraged our understanding of these pathways to explain the mechanisms underlying the observed genotype-environment-phenotype interactions.

Many of the QTL associations we identified are only detectable in the presence of specific stimuli, which underscores the need to activate cells to capture additional genotype-phenotype relations (3941). Furthermore, few eQTL studies have definitively identified the causal variants underlying the associations. By measuring variability in hundreds of individuals, applying stimuli that partially overlap in their downstream pathways, and leveraging genomic data sets such as the 1000 Genomes Project and ENCODE (29, 42), we (i) pinpointed causal variants in reQTL regions and (ii) identified the stimulus-activated TFs that bind differentially at these SNPs, which explain some of the G × E interactions. Our data set can thus be used to explore mechanisms of G × E interactions (43) and are consistent with deoxyribonuclease I sensitivity QTL (44) and ChIP-seq QTL (45) studies that showed that differential TF binding between individuals is pervasive in resting cells.

The reQTLs we identified provide genetic explanations for interindividual variation in innate immune responses. This is best exemplified by the trans-reQTL in the IRF7 locus that regulates the type I IFN antiviral response. Our study reveals the effects of this trans-reQTL on target genes (antiviral IFN module) in the context of a particular cell type (DCs) and in response to specific ligands (influenza). The changes in this immune response are, in turn, likely to affect organismal phenotypes that are driven by the IFN module, including susceptibility to viral infections and autoimmune diseases like SLE.

Overall, our high-throughput experimental pipeline and integrative analysis of primary human DCs reveals abundant gene-by-environment interactions, points to the effects of disease variants on pathogen detection, and motivates extending our approach to a wide variety of immune cell types and stimulatory conditions to better explain the impact of human genetic variation on the immune response.

Supplementary Materials

www.sciencemag.org/content/343/6175/1246980/suppl/DC1

Materials and Methods

Supplementary Text

Figs. S1 to S5

Tables S1 to S11

References (4672)

  • * These authors contributed equally to this work.

References and Notes

  1. Acknowledgments: We thank all study participants in the PhenoGenetic project for their contributions. We are grateful to Q. Sievers and C. Wu for advice about preparing primary cells and V. Agarwala and D. A. Landau for valuable discussions. This work was supported by the National Institute of General Medical Sciences, NIH, grant RC2 GM093080 (C.O.B. and N.H.); the National Institute of Allergy and Infectious Diseases grant U19 AI082630 and NIH Director’s New Innovator Award DP2 OD002230 (N.H.); the National Human Genome Research Institute grant P50 HG006193, NIH Director’s Pioneer Award DP1 CA174427, and Howard Hughes Medical Institute (A.R.); and NIH grant R01 HG004037 (M.K.). F.Z. is supported by an NIH Director’s Pioneer Award (DP1 MH100706), an NIH Transformative R01 grant (R01 DK097768), the Keck, McKnight, Merkin, Vallee, Damon Runyon, Searle Scholars, Klingenstein, and Simons Foundations, Bob Metcalfe, and Jane Pauley. M.N.L. is supported by the NIH Medical Scientist Training Program fellowship (T32 GM007753). T.R. is supported by NIH fellowship F32 AG043267. K.S. is supported by NIH grant T32 HG002295. Gene expression data are deposited in National Center for Biotechnology Information’s Gene Expression Omnibus under accession no. GSE53166.
View Abstract

Navigate This Article