Meta-analysis of shared genetic architecture across ten pediatric autoimmune diseases

Abstract

Genome-wide association studies (GWASs) have identified hundreds of susceptibility genes, including shared associations across clinically distinct autoimmune diseases. We performed an inverse χ(2) meta-analysis across ten pediatric-age-of-onset autoimmune diseases (pAIDs) in a case-control study including more than 6,035 cases and 10,718 shared population-based controls. We identified 27 genome-wide significant loci associated with one or more pAIDs, mapping to in silico-replicated autoimmune-associated genes (including IL2RA) and new candidate loci with established immunoregulatory functions such as ADGRL2, TENM3, ANKRD30A, ADCY7 and CD40LG. The pAID-associated single-nucleotide polymorphisms (SNPs) were functionally enriched for deoxyribonuclease (DNase)-hypersensitivity sites, expression quantitative trait loci (eQTLs), microRNA (miRNA)-binding sites and coding variants. We also identified biologically correlated, pAID-associated candidate gene sets on the basis of immune cell expression profiling and found evidence of genetic sharing. Network and protein-interaction analyses demonstrated converging roles for the signaling pathways of type 1, 2 and 17 helper T cells (TH1, TH2 and TH17), JAK-STAT, interferon and interleukin in multiple autoimmune diseases.

Figure 1

The ten pAID case cohorts and top pAID-association loci identified. (a) Percentage and relative contribution of cases for the ten pediatric autoimmune diseases studied. (b) Top pAID-association signals identified by inverse χ² meta-analysis. The top 27 loci (where at least one lead SNP reached genome-wide significance: P_META < 5 × 10⁻⁸) are annotated with the candidate gene symbol. (c) Novel and established pAID-association loci. Top left: rs706778 (chr10p15.1) is a known DNase I peak and an intronic SNP in IL2RA and was associated with THY, AS, PSOR, CEL, T1D and JIA. Top right: rs755374 (chr5q33.3) is an intergenic SNP upstream of IL12B and was associated with AS, CEL, UC and CD. Bottom left: rs2807264 (chrXq26.3), mapping near CD40LG, was associated with CEL, UC and CD, and chr15q22.33 (rs72743477), also mapping to an intronic position in SMAD3, was associated with UC, CD and AS. Bottom right: SNPs are colored according to pairwise LD (r²) with respect to the most strongly associated lead SNP in the locus. Associated pAIDs are indicated at the upper left. pAID associations are color-coded according to the key in each plot. (d) Pleiotropic candidate genes have pleiotropic effect sizes and directions across pAIDs. Although a few pleiotropic SNPs had consistent effect directions across diseases (e.g., IL21), for many loci (e.g., PTPN22 and CLEC16A), the candidate SNP had variable effect directions across diseases. The radii of the wedges correspond to the absolute values of the Z-scores (beta/s.e.) for each pAID, and the color indicates whether the SNP is protective (green) or risk-associated (red) for each disease.

Figure 2

Pleiotropic loci with heterogeneous effect directions across pAIDs. (a) Disease-specific Z-scores (beta/s.e.) for each SNP identified as having different effect directions across the ten pAIDs and as detailed in the figure. Circles (color-coded by disease as in key) denote diseases where the indicated SNP had an opposite effect compared with that of the group of pAIDs identified as sharing the lead association on the basis of results of the model search (black triangles). (b) Clustering of pAIDs across the lead loci on the basis of disease-specific effect sizes. Agglomerative hierarchical clustering across ten pAIDs on the basis of normalized directional Z-scores (beta/s.e.) resulting from logistic regression analysis in each disease for the 27 lead loci based on those disease combinations identified by the model search analysis as producing the strongest association-test statistics.

Figure 3

Integrated annotation of pAID-association loci using existing predictive and experimental data sets. (a) Biological, functional and literature annotations for the 27 loci reaching genome-wide significance in meta-analysis. Loci (identified by the lead SNPs and candidate genes) are organized by column; the colors in the table denote the associated pAIDs, functional annotations are presented at the top of the table, and the color bar at the bottom represents the meta-analysis P_meta values (according to key at right). For each locus, the lead SNP and proxy SNPs (r² > 0.8) were included in the annotation protocol (Online Methods). (b) Distribution and enrichment of experimental and predicted annotations for the top 27 GWS SNPs. The annotation frequencies were used to calculate the relative enrichment of pAID SNPs (blue bars) as compared with that of 10,000 random 100-SNP sets drawn from the genome in each annotation category. CpG, CpG islands; DNase, DNase-hypersensitivity I sites; gad, known genetic association; gerp_phast, conserved positions; mir, miRNAs; sift_pp, functional mutations in SIFT; tfbs, TF–binding sites.

Figure 4

Tissue-specific gene set enrichment analysis (TGSEA) of pediatric and adult autoimmune data sets identifies autoimmune-associated gene expression patterns across immune cells and tissues. (a) Expression enrichment of autoimmune-associated genes across human tissues. Distribution of TGSEA enrichment score (ES) values across 126 tissues for pAID-associated genes (center) either with (circles, top curve) or without (triangles, bottom curve) the extended MHC. Results for the pAID gene set are compared with those obtained for known genes associated with CD (left) and schizophrenia (right). Tissue and cell types are classified as immune (red) or non-immune (blue) and are ranked left to right on the basis of the magnitude of the ES test statistic. (b) Enrichment of pAID-associated gene expression across diverse murine immune cell types. Distribution of pAID-associated gene ES values across murine immune cell types either including (red) or excluding the genes within the MHC (gold); results are compared with those for genes associated with CD, schizophrenia (Schizo, turquoise), LDL cholesterol (LDL, magenta) or body mass index (BMI, blue) abstracted from the National Human Genome Research Institute (NHGRI) GWAS Catalog. (c) Hierarchical clustering based on the expression of pleiotropic candidate genes associated with three or more autoimmune diseases across the murine immune cells. Boxes outlined in black denote gene clusters enriched for specific disease associations discussed in the text. An enlarged version of c is presented in Supplementary Figure 8.

Figure 5

Genetic variants shared across the ten pAIDs reveal autoimmune disease networks. (a) Quantification of pAID genetic sharing by GPS test including SNPs within the extended MHC. Correlation plot of results of the pairwise pAID GPS test; the color intensity and the size of each circle are proportional to the strength of the correlation as the negative base ten logarithms of the GPS test P values (color-coded numbers in squares). (b) Quantification of autoimmune disease genetic sharing by locus-specific pairwise sharing. Undirected weighted network graph depicting results from the LPS test. Edge size represents the magnitude of the LPS test statistic; labeled nodes for each of the 17 autoimmune diseases are positioned on the basis of a force-directed layout. Edges represent significant pairs after Bonferroni adjustment (P_adj < 0.05). (c) Protein-protein interaction network analysis of the top pAID-associated protein candidates in STRING; action view of protein interactions observed across the top 46 GWM (P < 1 × 10⁻⁶) signals, of which 44 could be mapped to corresponding proteins. Views were generated on the basis of results for known and predicted protein interactions produced by the STRING DB Homo sapiens database. The plots shown are results of the ‘action’ view, where the molecular actions (stimulatory, repressive or binding) are illustrated by arrows.