Variation in the TLR10/TLR1/TLR6 locus is the major genetic determinant of interindividual difference in TLR1/2-mediated responses.

Toll-like receptor (TLR)-mediated innate immune responses are important in early host defense. Using a candidate gene approach, we previously identified genetic variation within TLR1 that is associated with hyper-responsiveness to a TLR1/2 agonist in vitro and with death and organ dysfunction in patients with sepsis. Here we report a genome-wide association study (GWAS) designed to identify genetic loci controlling whole blood cytokine responses to the TLR1/2 lipopeptide agonist, Pam(3)CSK(4) (N-palmitoyl-S-dipalmitoylglyceryl Cys-Ser-(Lys)(4)) ex vivo. We identified a very strong association (P<1 × 10(-27)) between genetic variation within the TLR10/1/6 locus on chromosome 4, and Pam(3)CSK(4)-induced cytokine responses. This was the predominant association explaining over 35% of the population variance for this phenotype. Notably, strong associations were observed within TLR10, suggesting that genetic variation in TLR10 may influence bacterial lipoprotein-induced responses. These findings establish the TLR10/1/6 locus as the dominant common genetic factor controlling interindividual variability in Pam(3)CSK(4)-induced whole blood responses in the healthy population.

Introduction

The innate immune system provides early recognition of microbial pathogens important to host defense. Toll like receptors (TLRs) play a key role in host defense, providing a mechanism to respond to highly conserved pathogen-associated molecular patterns (PAMPs).[1] In humans, there are ten unique TLR genes coding for receptors that initiate responses to PAMP ligands a robust inflammatory response. TLR2 heterodimerizes with TLR6, TLR1, and possibly TLR10, and these combinations facilitate the recognition of multiple distinct bacterial patterns diversifying innate immune sensing.[2-4] The importance of TLR2 in host defense has been well-established in mice where its deficiency has been associated with increased susceptibility to mycobacterial infection, pneumococcal meningitis, and sepsis due to Staphylococcus aureus and Listeria monocytogenes.[5-8] TLR1/2 and TLR2/6 heterodimers can discriminate the acylation of bacterial lipopeptides recognizing triacyl- and diacyl-lipopeptides respectively.[2,9-11] The synthetic triacyl lipopeptide N-palmitoyl-S-dipalmitoylglyceryl Cys-Ser-(Lys)4 (Pam3CSK4) and diacyl lipopeptide Fibroblast Stimulating Ligand-1 (FSL-1) derived from Mycoplasma salivarium have been shown to stimulate via TLR1/2 and TLR2/6 heterodimers.[2,12] Additionally, TLR2/6 heterodimers recognize peptidoglycan (PGN) and a yeast cell wall particle, Zymosan.[13,14]

A role for TLR1/2 and TLR2/6 in human disease has been suggested by candidate gene studies. We and others have demonstrated that there exists high inter-individual variability in terms of human leukocyte inflammatory responses to PAMPs[15,16] and that a portion of this variability is attributable to common genetic variants. Genetic variation in TLR2 has been shown to confer reduced responses to peptidoglycan and heat-killed S. aureus in vitro.[17] More recently, we have demonstrated that variants in TLR1 are highly associated with Pam3CSK4-induced whole blood cytokine production. We reported that common genetic variants in TLR1 conferred marked hyper-responsiveness to Pam3CSK4 and these same variants were associated with increased risk of organ dysfunction and death in septic shock.[15,18] Other studies have demonstrated associations between genetic variation in TLR1 with susceptibility to leprosy and tuberculosis.[19,20] These data support a role for TLR1/2-mediated responses in human disease. However, to date, our understanding of the role for genetic variation in TLR-mediated responses has been based on targeted candidate gene studies. Thus, in order to more comprehensively assess the genetic factors controlling TLR2-mediated responses in the healthy human population we undertook a genome wide association study to identify loci modifying Pam3CSK4-induced cytokine production in whole blood ex vivo.

Results

We employed samples from 360 healthy Caucasian subjects who had an average age of 35±14 years and were 39% male. Given that many innate immunity genes demonstrate population differences in allele frequencies including the genes coding for TLRs,[21] we performed principal components analysis (PCA) to address the possibility that there might exist population admixture within our genotyped subjects. PCA revealed that subjects who self-reported as Caucasian cluster with Caucasians from Utah (CEU) and the Toscani in Italia (TSI) populations from the HapMap3 collection[22] (Supplemental Figure 1). However, we did identify associations between eigenvalues from the first three principal components and TLR agonist-induced cytokine production and so these eigenvalues were used as covariates in the multiple linear regression models for the GWAS.

We used a genome wide association test adjusted for age, gender and the first 3 principal components, and identified 19 SNPs within the TLR10/1/6 locus on chromosome 4 that were associated with Pam3CSK4 induced IL-6 (Figure 1A), IL-1β, and TNF-α production in whole blood (Supplemental Figure 2) at a genome-wide level of significance (p = 1 × 10−8 – p = 1 × 10−27) (Table 1). No other loci achieved associations at a genome-wide level of significance including SNPs found in genes involved in TLR1/2 signaling such as TIRAP, IRAK4, and IRAK1 that we had anticipated a priori would be associated with the cytokine induced phenotypes (Table 2). Notably, all cytokine values obtained from the whole blood assay were normalized to a monocyte count obtained from the donor at the time of phlebotomy. In this way we mitigated the chances of identifying variation that merely affected the number of circulating monocytes.

Figure 1

The TLR10/1/6 locus is highly associated with Pam3CSK4-induced IL-6. A) Manhattan plot showing the primary association statistics for the Pam3CSK4–induced IL-6 concentration across all chromosomes for the 483,197 genotyped SNPs. Embedded quantile-quantile plot of −log10(P value) vs. the expected −log10(P value) for SNP associations with Pam3CSK4-induced IL-6 phenotype. B) Similar association statistics for FSL-induced IL-6 concentration.

Table 1

Top ranked associations with Pam3CSK4-induced IL-6[1]

SNP	P-value2	Gene	Left Gene	Right Gene	Allele	MAF	β3
rs4543123	1.35e-27	intergenic	TLR10	TLR1	C/T	0.23	+0.36
rs4833095	4.29e-26	TLR1	TLR10	TLR6	C/T	0.24	+0.35
rs11466640	1.85e-25	TLR10	KLF3	TLR1	C/T	0.18	+0.37
rs5743563	3.66e-25	TLR1	TLR10	TLR6	C/T	0.20	+0.36
rs11096956	5.57e-25	TLR10	KLF3	TLR1	G/T	0.22	+0.35
rs1873195	3.46e-19	FAM114A1	TLR6	TMEM156	A/G	0.22	+0.30
rs10008492	1.51e-17	intergenic	KLF3	TLR10	C/T	0.35	+0.26
rs4331786	1.72e-16	intergenic	KLF3	TLR10	A/G	0.31	+0.26
rs11096957	1.85e-16	TLR10	KLF3	TLR1	A/C	0.31	+0.26
rs10024216	7.25e-15	intergenic	KLF3	TLR10	A/G	0.30	+0.25
rs2890620	2.70e-13	intergenic	KLF3	TLR10	C/T	0.30	+0.24
rs6824105	1.45e-12	intergenic	KLF3	TLR10	A/G	0.29	+0.23
rs7660102	3.81e-11	intergenic	KLF3	TLR10	A/C	0.22	+0.22
rs4833103	5.87e-11	intergenic	TLR1	TLR6	A/C	0.49	+0.19
rs11944159	9.74e-10	FAM114A1	TLR6	TMEM156	C/T	0.27	−0.19
rs17429245	2.01e-09	intergenic	TLR6	FAM114A1	C/T	0.26	−0.19
rs6824001	2.53e-09	FAM114A1	TLR6	TMEM156	C/T	0.27	−0.18
rs13132956	8.93e-09	intergenic	KLF3	TLR10	C/T	0.20	+0.21
rs6833914	9.74e-09	intergenic	TLR6	FAM114A1	C/T	0.29	−0.18
rs4832792	1.05e-08	FAM114A1	TLR6	TMEM156	A/C	0.29	+0.18

All SNPs meeting genome-wide significance (p<1e-8).

Adjusted for gender and eigenvalues from first 3 principal components.

Effect size and direction associated with copy number of minor allele (change in mean log10[IL6] with each copy of minor allele).

Table 2

Genes anticipated a priori to be associated with Pam3CSK4-induced IL-6 phenotype1

Gene1	Chromosome	Loci (Mb)2	SNP3	SNP Gene4	P-value
TLR2	4	154.84–154.85	rs2405432*	RNF175	0.17
NFKB1	4	103.54–103.76	rs2085548	intergenic	0.25
CD14	5	139.99–139.99	rs1583005	IK	0.02
MYD88	3	38.16–38.16	rs9825655	DLEC1	0.02
IRAK4	12	42.43–42.46	rs7972025	intergenic	0.04
LY96	8	75.07–75.10	rs10504553	intergenic	0.04
TIRAP	11	125.66–125.66	rs478309	FOXRED1	0.11

TLR and TLR signaling genes anticipated to be associated with the agonist-induced cytokine concentration.

For each gene, a window 50Kb from either end of the gene was included to select the most highly associated SNP.

SNP most highly associated within the gene range. Asterisk signifies the SNP was imputed.

Gene in which the SNP was located.

We next sought to identify loci associated with responses to TLR2/6 ligands FSL-1, PGN, and Zymosan in 167 subjects for whom we had measured whole blood responses to these ligands. We did not identify any associations reaching genome-wide significance (Figure 1B) and, notably, no SNPs within the TLR10/1/6 locus or TLR2 were even nominally associated (p>0.05) with responses to these ligands. Nonetheless, there were several moderately strong associations detected at other genomic loci with these cytokine responses ranging from p=1.55×10−6 (Zymosan-induced IL-6), p=3.30×10−6 (FSL-induced IL-6) to p=4.37×10−6 (PGN-induced IL-6). Since these analyses included fewer subjects than the GWAS of Pam3CSK4-induced responses we re-ran the GWAS of Pam3CSK4–induced IL-6 using only these 167 subjects. This analysis still identified multiple SNPs that were associated at a genome-wide level of significance (p<4.7×10−12) demonstrating that while statistical power for this sub-study may have been limiting, the associations with Pam3CSK4–induced responses are orders of magnitude stronger than any associations with TLR2/6 agonist-induced responses.

In order to identify SNPs within the TLR10/1/6 locus not directly genotyped by our platform that may be driving the observed associations with Pam3CSK4-induced cytokine production, we used imputation to infer missing genotypes on chromosome 4 using 1000 genomes NCBI Build 37[23] as a reference population. These imputed SNPs were tested for association with the Pam3CSK4-induced cytokine phenotypes. We observed a 222kb region across the TLR10/1/6 locus that was associated with Pam3CSK4-induced IL-6 at a genome wide level of significance (Figure 2). The SNP most highly associated with hypermorphic responses was rs67719080 (p=1×10−27), an intergenic SNP between TLR10 and TLR1. Of the SNPs that fell within genes, SNPs within TLR10 were most highly associated with hypermorphic cytokine responses (Figure 2). The most highly associated TLR10 coding SNP was rs4129009 (TLR102323A/G), a non-synonymous polymorphism that causes an amino acid change in the highly conserved Toll/Interleukin-1 receptor (TIR) domain. Individuals homozygous for the rare allele had increased IL-6 production consistent with a hypermorphic response (Figure 3). In addition to the TIR domain SNP, we also identified a missense SNP in TLR10, rs11096955 (I369L), near leucine-rich repeat 9 (LRR9: aa 349–368) of TLR10 that was strongly associated with hypermorphic responses to Pam3CSK4 (p=5.36×10−16).

Figure 2

Fine mapping of associations in TLR10/1/6 locus with imputed genotypes. Association statistics for the imputed SNPs on chromosome 4 versus the −log10(P value) of the Pam3CSK4-induced IL-6 phenotype with associated LD plot. Area of focus is the TLR 10/1/6 locus and highly associated SNPs are shown as black squares. The SNPs shown by rs number are the most highly associated coding SNPs in each gene. The lower plot shows all of chromosome 4 with the gray box representing cytoband 4p14.

Figure 3

Minor alleles in TLR1 and TLR10 are associated with hypermophic effects on Pam3CSK4-induced IL6. Coding SNPs for TLR1 (A, B) and TLR10 (C) most highly associated with Pam3CSK4-induced IL6 showing hypermorphic responses with the rare genotype.

Coding SNPs within TLR1 were also highly associated with the Pam3CSK4-induced cytokine phenotype including rs4833095 (TLR1742A/G) and rs5743618 (TLR11805G/T) but were not in high LD with the TLR10 coding SNP rs4129009 (Table 3) suggesting a distinct association. Notably, rs5743551 a SNP found 5′ to TLR1 that we have previously shown to be highly associated with death and organ dysfunction in sepsis was also highly associated (p=2.8×10−24). Finally, we also found a strong association with a non-synonymous variant in TLR6 (rs5743818, TLR61932T/G) and Pam3CSK4-induced responses (p=1.28×10−9). This SNP was not found to be in high linkage disequilibrium with the other most-highly associated coding SNPs in TLR1 (R2=0.11) and TLR10 (R2=0.08) (Table 3).

Table 3

Coding SNPs in TLR10/1/6 locus most-highly associated with Pam3CSK4-induced responses1

SNP	P-Value2	Gene	Alleles	Feature	LD3
rs4833095	1.15e-25	TLR1	C/T	Missense: S248N	1.0
rs4129009	5.04e-25	TLR10	A/G	Missense	0.66
rs4274855	5.04e-25	TLR10	A/G	utr-5′	0.66
rs5743566	9.71e-25	TLR1	C/G	utr-5′	0.73
rs5743565	9.71e-25	TLR1	A/G	utr-5′	0.73
rs9715841	1.61e-24	TLR10	C/T	utr-3′	0.61
rs10776482	1.61e-24	TLR10	C/T	coding-synonymous	0.61
rs10776483	1.61e-24	TLR10	C/T	coding-synonymous	0.61
rs11096956	1.61e-24	TLR10	G/T	coding-synonymous	0.61
rs5743618	7.10e-24	TLR1	G/T	Missense: I602S	0.86
rs11466661	5.36e-16	TLR10	A/C	utr-3′	0.61
rs11096955	5.36e-16	TLR10	A/C	Missense: L369I	0.61
rs11096957	5.36e-16	TLR10	A/C	Missense: N241H	0.61
rs5743818	1.28e-09	TLR6	G/T	coding-synonymous	0.11

Most highly associated coding SNPs to PAM3CSK4-induced IL-6.

Adjusted for age, gender, and eigenvalues from first 3 principal components.

Linkage disequilibrium (R2) between each SNP and the highest TLR1 coding SNP rs4833095

Discussion

In this genome-wide association study, we found that the TLR10/1/6 region on chromosome 4 is the dominant common genetic locus controlling inter-individual variation in responses to Pam3CSK4 in whole blood from healthy subjects ex vivo.

While the genes coding for TLRs are distributed throughout the genome, TLR10, TLR1, and TLR6 cluster at a locus on chromosome 4p14. Evidence suggests that this tandem arrangement arose from a gene duplication event.[24] Notably, all three of these genes have significant allelic heterogeneity with an abundance of rare variants that may indicate an influence of purifying selection.[24] In addition, there exist significant geographic differences in genetic variation between European populations within the TLR10/1/6 locus.[21] However, our principal components analysis shows that our subjects clustered with Caucasian populations in HapMap3 and our adjustment with principal components in the linear regression suggests that the association testing is not confounded by cryptic population substructure.

Among the SNPs within TLR1 showing the strongest associations in our study were several that have been previously associated with susceptibility to leprosy (rs5743618)[25], risk for prostate cancer and placental malaria (rs4833095).[26,27] These findings are consistent with the assertion that functional responses mediated by TLR1/2 heterodimers might drive important biologic responses and alter risk for disease. We were more surprised to find strong associations with coding SNPs within TLR10 as there is no known ligand specific for TLR10 and it is not known that TLR10 ligation actually generates an intracellular response.[4,28] These findings suggest that SNPs within TLR10 may contribute to associations between disease susceptibility and the TLR10/1/6 locus.

The most highly associated non-synonymous SNP in TLR10, rs4129009 causes an amino acid change in the TIR domain of the intracellular portion of the protein. The TIR domain is critical for intracellular signaling in other TLR family members.[29,30] A recent study has shown that a chimeric receptor containing the extracellular domain of TLR10 and the intracellular domain of TLR1 (including the TIR domain) induced a cellular response to Pam3CSK4 comparable to wild-type TLR1.[4] This study suggests that the extracellular portion of TLR10 recognizes Pam3CSK4 but that the intracellular portion of TLR10 does not translate this recognition event to an intracellular signal. Our study shows that individuals homozygous for the rare allele of rs4129009 in TLR10 have increased cytokine responses suggesting that this genetic alteration of the TIR domain may result in a functionally active TLR10 molecule. Of note, this SNP has previously been reported to be associated with decreased risk of atopic asthma.[31] In addition to this SNP in the TIR domain, we identified another highly associated missense SNP in TLR10, rs11096955 (I369L), near LRR9, that could alter ligand binding. In order to best identify whether the TLR10 signal is an independent association, future research should be aimed at other racial groups where haplotype blocks in these region are smaller. Future work will need to more finely delineate whether SNPs in TLR10 or TLR1 (or both) are causally responsible for the associations observed. However, due to moderate LD, conditional regression analysis adjusting for the top SNPs in this analysis was underpowered to detect independent associations.

The importance of genetic variation in TLR genes and downstream TLR signaling genes is highlighted by candidate gene studies that have demonstrated associations between variants in these genes and diseases for which host defense and inflammation is pathologic. With respect to genes encoding the TLR1/2 heterodimer, functional polymorphisms within the TLR10/1/6 locus and TLR2 have been associated with altered susceptibility to the mycobacterial infections of leprosy and tuberculosis.[19,20,32] A TLR1 polymorphism (rs5743618, Ser602Ile) that mediates higher levels of signaling and cell surface expression[15,19] is associated with protection from recurrent urinary tract infection and pyelonephritis.[33] In sepsis, where severe infection leads to overwhelming inflammation and end-organ dysfunction, a TLR 1 polymorphism (rs5743551) associated with marked hyper-responsiveness has been associated with risk of death and organ dysfunction and sepsis induced acute lung injury.[15,18] Outside of infectious diseases, polymorphisms within the TLR10/1/6 locus have been variably associated with prostate cancer, non-Hodgkin lymphoma, Crohn’s disease, asthma, and chronic sarcoidosis.[26,31,34-40] Our findings that the TLR10/1/6 locus explains a large portion of population variance in TLR1/2-mediated responses in vitro provides additional support for the importance of this locus in human disease.

Several previous reports have demonstrated associations between disease risk and genetic variation in TLRs and genes of the TLR intracellular signaling pathway including TLR2, TIRAP, IRAK4, and IRAK1.[41-43] In spite of these previous findings, we detected only a nominally significant association with variants in some TLR-related genes (Table 2). It should be noted that this study was designed to have adequate statistical power to detect associations with common genetic factors (MAF >5%). This study is inadequately powered for detection of associations with rare genetic variants (MAF<1%) and, therefore, we cannot exclude the possibility that rare variants within these or other genes may also play a role in modulating these effects. Nonetheless, our findings suggest that common genetic variation in TLR pathway genes outside of the TLR10/1/6 locus play only a minor role in modifying TLR1/2 responses in the Caucasian population.

In summary, our study shows that genetic variation within the TLR10/1/6 locus is the major common genetic factor explaining inter-individual variation in TLR1/2-mediated cytokine responses to Pam3CSK4 in vitro. We find that the mostly highly-associated SNPs fall within TLR10 and that some of these SNPs are located in or near important functional domains (TIR domain and LRR9) of TLR10 suggesting that this receptor might have functional relevance. Overall, this study supports ongoing efforts to understand the importance of this locus to human diseases involving innate immunity.

Materials and Methods

Study Subjects

We used DNA samples and innate immune response phenotypes collected from 360 healthy Caucasian volunteers recruited from the Seattle metropolitan area from whom written informed consent was obtained. This was approved by the University of Washington Human Subjects Committee. This population has been previously described by our group.[15]

Cytokine Assays

Innate immune responses were measured in whole blood ex vivo as previously described.[16] We exposed whole blood collected from each subject to Pam3CSK4 (360 subjects at 100ng/ml), FSL-1, PGN (167 subjects at 100ng/ml) and Zymosan (179 subjects at 100μg/ml) for six hours, supernatants were collected, and production of Tumor Necrosis Factor-alpha (TNF-α), Interleukin-1 beta (IL-1β), Interleukin-6 (IL-6), Interleukin-8 (IL-8), Interleukin-10 (IL-10), Granulocyte colony stimulating factor (G-CSF), Interleukin 1-receptor antagonist (IL-1ra), monocyte-chemotactic protein-1 (MCP-1) was measured by cytometric bead-based immunoassay (Luminex™). A complete blood count with differential cell counts was obtained at the time of blood sampling for the stimulation assays and cytokine concentrations were normalized to monocyte counts.

Genotyping and Imputation

Genomic DNA was genotyped using the Illumina™ Human 1M Beadchip array. In addition, we imputed genotypes on chromosome 4 not present on the array with the BEAGLE software package version 3.3[44] using EUR genotypes from 1000 Genomes[23] as a reference.

Quality Control

Quality control was performed as described by Anderson et al.[45] We assessed for discordance between reported sex and genotype-determined sex, excess autosomal heterozygosity, excess relatedness (identity by descent of > 0.1875), and population substructure using principal components analysis (PCA) and removed 14 subjects resulting in a total of 346 subjects. All subjects had a genotype call rate of over 97%. The 561,491 SNPs were filtered to remove all SNPs with a minor allele frequency (MAF) <0.05, Hardy-Weinberg equilibrium p<0.001, or a call rate ≤ 0.90 resulting in 493,197 SNPs that were used for association testing. Imputed SNPs for chromosome 4 were filtered for an allelic R2 of 0.85.

Data analysis

We tested for associations between genome-wide genotypes and log10-transformed, monocyte normalized, cytokine values by multiple linear regression assuming additive effects. Subjects and SNPs passing QC filtering were tested for association with Pam3CSK4-induced, monocyte-normalized, whole blood cytokine production adjusting for covariates including age, gender and eigenvalues from the first three principal components generated by PCA clustering subjects with samples from HapMap3 (Release 3, NCBI build 36).[22] Correcting for multiple tests, we considered a p< 1 × 10−8 to be indicative of genome-wide significance. We assigned p-values to TLR signaling genes anticipated a priori to be associated with the cytokine phenotype by choosing the p-value of the highest SNP within a 50kB range from the 5′ and 3′ end of the gene. All above analyses were performed and linkage disequilibrium calculated using the Golden Helix™ software package.

Supplemental Figure 1. Study subjects cluster with European populations by principal component analysis. X-Y scatter plot of the first two principal components. Healthy subjects are shown in light blue circles and the HAPMAP3 racial phenotypes are in the other colored squares.

Supplemental Figure 2. Dominant peak of association shared across different Pam3CSK4-induced cytokine phenotypes. Manhattan plot showing the primary association statistics for the Pam3CSK4–induced IL-1β concentration (A) and TNF-α concentration (B) across all chromosomes for the 483,197 genotyped SNPs. Quantile-quantile plot of −log10(P value) vs. the expected −log10(P value) for SNP associations with Pam3CSK4-induced IL-1β (C) and TNF-α (D) phenotype.