Targeted genomic analysis reveals widespread autoimmune disease association with regulatory variants in the TNF superfamily cytokine signalling network

Abstract

Background: Tumour necrosis factor (TNF) superfamily cytokines and their receptors regulate diverse immune system functions through a common set of signalling pathways. Genetic variants in and expression of individual TNF superfamily cytokines, receptors and signalling proteins have been associated with autoimmune and inflammatory diseases, but their interconnected biology has been largely unexplored.

Methods: We took a hypothesis-driven approach using available genome-wide datasets to identify genetic variants regulating gene expression in the TNF superfamily cytokine signalling network and the association of these variants with autoimmune and autoinflammatory disease. Using paired gene expression and genetic data, we identified genetic variants associated with gene expression, expression quantitative trait loci (eQTLs), in four peripheral blood cell subsets. We then examined whether eQTLs were dependent on gene expression level or the presence of active enhancer chromatin marks. Using these eQTLs as genetic markers of the TNF superfamily signalling network, we performed targeted gene set association analysis in eight autoimmune and autoinflammatory disease genome-wide association studies.

Results: Comparison of TNF superfamily network gene expression and regulatory variants across four leucocyte subsets revealed patterns that differed between cell types. eQTLs for genes in this network were not dependent on absolute gene expression levels and were not enriched for chromatin marks of active enhancers. By examining autoimmune disease risk variants among our eQTLs, we found that risk alleles can be associated with either increased or decreased expression of co-stimulatory TNF superfamily cytokines, receptors or downstream signalling molecules. Gene set disease association analysis revealed that eQTLs for genes in the TNF superfamily pathway were associated with six of the eight autoimmune and autoinflammatory diseases examined, demonstrating associations beyond single genome-wide significant hits.

Conclusions: This systematic analysis of the influence of regulatory genetic variants in the TNF superfamily network reveals widespread and diverse roles for these cytokines in susceptibility to a number of immune-mediated diseases.

Keywords: Autoimmunity; Autoinflammation; GWAS; Gene set analysis; Genetics; Genomics; TNF superfamily; eQTL.

Fig. 1

Flow chart of analyses. Flow chart demonstrates how results from each analysis feed into the next. Datasets analysed are listed in blue italics

Fig. 2

Expression of TNFSF-related genes differs across leukocyte subsets. Expression of TNFSF-related genes was measured across four cell subsets from five healthy controls by microarray. Expression values are hierarchically clustered. Cell types are coloured blue (CD4⁺ T cells), purple (CD8⁺ T cells), green (CD14⁺ monocytes) and red (CD16⁺ neutrophils). Genes are grouped by function and coloured yellow (TNFSF member ligands), orange (TNFRSF member receptors) and black (adaptors and signalling molecules in TNFSF signalling network)

Fig. 3

TNFSF-related genes are under extensive genetic regulation. a Normalised FAS expression in each cell subset is plotted against rs4406737 genotype. Association p values are indicated for eQTLs with FDR < 0.1. b Normalised TNFRSF10B and TNFRSF10C expression in monocytes and neutrophils is plotted against rs7009522 genotype. NS not significant. c TNFSF-related genes with a significant cis-eQTL (FDR < 0.1) in any cell type were extracted. For each gene, the SNP most significantly associated with expression in each cell type was extracted (best cis-eQTL) and the FDR corresponding to its p value in that dataset calculated. The plot depicts hierarchical clustering of the negative logarithm of these FDRs. Colours are as in Fig. 2. d The number of SNPs in cis-eQTL models after variable selection. Numbers around the periphery and grey shades indicate the number of eQTL SNPs remaining as predictors in the model for each gene. “n” indicates the number of genes with cis-eQTLs represented by each pie chart

Fig. 4

eQTLs are independent of total magnitude of gene expression and not preferentially associated with active enhancer marks. a Expression of TNFSF and TNFRSF members was measured by the NanoString nCounter analysis system. Each point represents average expression over eight or more individuals, including healthy controls and individuals with IBD. Genes are listed left to right in order of decreasing number of cell types in which an eQTL was detected. b Average H3K27ac ChIP-seq or input DNA sequencing counts per million intersecting TNFSF-related eQTL SNPs in the same cell type are plotted. c Average H3K27ac ChIP-seq counts per million intersecting eQTL SNPs in the same cell type are compared with a random distribution of intersections created by selecting the same number of SNPs from the cis genomic regions used for eQTL search. Error bars represent the standard deviation of 10,000 iterations of random selection. d eQTL chi-squared scores from the strongest eQTL SNP for each TNFSF-related gene (regardless of whether the association passed our eQTL significance threshold) are compared with H3K27ac ChIP-seq or input DNA sequencing counts per million at the same SNPs. Spearman correlation coefficient (rho) and correlation p values (p) are indicated for H3K27ac ChIP-seq counts per million versus eQTL score

Fig. 5

Immune-mediated disease risk alleles can either increase or decrease TNFSF-related gene expression. a Autoimmune and autoinflammatory disease GWAS hits tagged by TNFSF-related cis-eQTLs were identified and the directions of effect of the risk alleles on expression are plotted. Disease-associated SNPs that are eQTL SNPs in multiple cell types or are associated with multiple diseases are counted only once. “n” indicates the number of SNPs in each slice of the pie. b TNFSF14 expression is plotted against MS-associated SNP rs1077667 genotype. Allele A is protective. P values are provided for eQTLs with FDR < 0.1; NS indicates not significant. c Disease-associated eQTL SNPs depicted in a are plotted by eQTL cell type and disease association. Effect directions are coloured as in a. SNPs associated with more than one disease are plotted once per disease

Fig. 6

Genes regulated by disease-associated eQTL SNPs differ across diseases. Within each disease, permutation-based p values for gene-level disease association were calculated by combining the strongest significant cis-eQTL SNPs in each cell type. As in Table 1, proxy eQTL SNPs in each genetic dataset were filtered for relative independence before computation of a permutation-based disease association p value. P values were corrected within each disease by the Benjamini–Hochberg FDR method. The heatmap represents the negative logarithm of these corrected values such that genes marked in white-to-red shades show disease association, FDR < 0.1. Grey indicates no data available because the GWAS dataset did not include SNPs that tagged eQTLs for this gene with LD r² ≥ 0.8. Gene colours on the left side correspond to TNFSF ligands (yellow), TNFRSF receptors (orange) and adaptors/signalling molecules (black) as in Fig. 2. Genes are presented in order of decreasing association with any disease