Intersection of population variation and autoimmunity genetics in human T cell activation

Abstract

T lymphocyte activation by antigen conditions adaptive immune responses and immunopathologies, but we know little about its variation in humans and its genetic or environmental roots. We analyzed gene expression in CD4(+) T cells during unbiased activation or in T helper 17 (T(H)17) conditions from 348 healthy participants representing European, Asian, and African ancestries. We observed interindividual variability, most marked for cytokine transcripts, with clear biases on the basis of ancestry, and following patterns more complex than simple T(H)1/2/17 partitions. We identified 39 genetic loci specifically associated in cis with activated gene expression. We further fine-mapped and validated a single-base variant that modulates YY1 binding and the activity of an enhancer element controlling the autoimmune-associated IL2RA gene, affecting its activity in activated but not regulatory T cells. Thus, interindividual variability affects the fundamental immunologic process of T helper activation, with important connections to autoimmune disease.

Fig. 1. A systematic approach to characterize the variation in T cell activation and response

(A) Overview of study, from (1) cell purification and activation in four conditions, (2) genomewide or (3) signature expression profiling, (4) computational analysis to decompose expression variance and identify eQTLs and (5) functional validation of potential causal variants. (B) Time course expression profiles of expression in unbiased α3+28 cells (fold change relative to 0 hr baseline, top), and response to IFNβ (fold change relative to unbiased α3+28, middle) or to Th17 polarization cytokines (fold change relative to unbiased α3+28, bottom). Each row (gene) is normalized to its mean fold change. (C) Profiles of gene expression for seven gene clusters, for conditions outlined and color-coded as in A. A sample of representative genes in each cluster is listed. (D) Expression of 16 cytokines across 348 individuals and 5 conditions (black = unstimulated, other color codes as in A).

Fig. 2. Inter-individual variation in T cell response

(A) Relation between expression values (log2 scale) for pairs of Th defining cytokines (IL4, IL17F and IFNγ) in unbiased 4 hr α3+28 and 48 hr α3+28 conditions, each point is an individual sample. (B) Heatmap of 16 cytokines clustered by expression across 348 individuals in 48 hr α3+28 with demographic covariates race, age and sex. (C) Pairwise correlation between 16 cytokines in covariate-adjusted data during early activation (left), late activation (middle), and late activation with PCA-adjusted data to account for overall responsiveness (right).

Fig. 3. Sources of variation across T cell activation conditions

(A) Top 56 non-cytokine genes or (B) 16 key cytokines ranked by maximum repeatability across all conditions estimated from 8 recalled individuals. For each gene, repeatability (medium shading), contributions of cis heritability (dark shading) and of physiological covariates (age, gender, body-mass index, weight, height, diastolic, systolic, season and month; light shading) to repeatable variation are estimated (European cohort only). Cis heritability is defined as the proportion of phenotypic variance explained by variants within 1MB of the gene estimated using linear mixed models (29). Physiological variance explained is the proportion of variance explained by all known covariates including gender, age, weight, height, BMI, blood pressure estimated from multiple regression. Error bars indicate bootstrap standard errors around each estimate. (C) Percent difference of average population expression (median) from overall average (median) that show population differentiation in expression. (D) Fold change (effect size) in normalized expression for donors of European (left, N = 183) and Asian (right, N = 74) ancestry relative to donors of African ancestry (N=91), comparing 48 hrs α3+28 data adjusted for all PCs (x-axis) or for experimental and physiological covariates only (y-axis). Darker shading indicates more statistically significant population differentiated genes as determined by an omnibus F test.

Fig. 4. Common cis variation associated with gene expression during T cell activation and differentiation

(A) Manhattan plot of the significance of cis (within 1M of gene) associations for the 236 genes in each condition. Permutation significant (FDR < 0.01) associations are colored. (B–F) Individual plots of expression per genotype for illustrative SNPs, numbers in each panel indicate the significance (−log10(p-value), upper number) and effect size (beta, lower number) of association. (G) Per genotype plots of IL2RA expression in 48 hr α3+28, distinguished by population. Width of boxes representing each genotype is scaled by its frequency in the population.

Fig. 5. Trans-ancestry fine mapping identifies activation specific cis-eQTLs in GWAS regions

(A) In 48 hr α3+28/Th17, rs2863204 (MAF = 0.53, European) is best associated with IL23R expression and is independent from rs11209026 (R381Q). (B) At IL23R, a coding variant (rs11209026, R381Q, MAF = 0.06, European) and a regulatory variant (rs12095335) are associated with Crohn’s disease (CD). (C) Conditioning on rs2863204 recovers rs12095335 (secondary CD association) as a significant secondary association. (D) At the IL2RA locus, the best T1D association is a regulatory variant, rs12722495. (E) In 48 hr α3+28, rs12251836 (MAF = 0.41, European) is best associated with IL2RA expression and is independent from rs12722495. (F) Conditioning on rs12251836 recovers rs12722495 (primary T1D association). X-axis defines genomic intervals local to each gene and y-axis shows the −log10 p-value of association in GWAS (CD, T1D) and eQTL analysis (48 hr α3+28 or 48 hr α3+28/Th17).

Fig. 6. Validation of rs12251836, a causal variant in the IL2RA locus

(A) Genomic locations of five top associated SNPs in the IL2RA locus, together with mammalian sequence conservation and DNAase hypersensitivity sites (in Th1 cells). The ~200bp fragments centered around each of the variants used in enhancer assays are marked (I–V). (B) Enhancer activity in activated Jurkat cells for both alleles (as marked) of the 5 selected fragments cloned into a minimal promoter vector (firefly luciferase, all values normalized to co-transfected Renilla). (C) Enhancer activity of fragment IV, as in B, with or without PMA+Ionomycin (P+I) stimulation. (D) Protein binding to the two alleles of rs12251836 tested by EMSA with nuclear extract from Jurkat or human primary CD4⁺ T cells – these gels were run for different times and apparent sizes should not be compared. (E) Specific binding of YY1 to the rs12251836^T allele in nuclear extract of primary CD4⁺ T cells, with or without prior 6 hr activation with α3+28 (F) Top: predicted TF-binding sites around the two alleles of rs12251836. Bottom: pull-down assay of the three main predicted TFs (transfected into HEK293 cells, HA or FLAG tagged) with oligonucleotides encompassing both rs12251836 alleles or an irrelevant control (EBNA), and detected by immunoblotting for the relevant tags. (G) Trans-activating activity of over-expressed YY1 in Jurkat cells also transfected with the enhancer reporter plasmids used in B (with either allele at rs12251836). Data are representative of three or more independent experiments; error bars = SD.