The Discovery cohort included 814 JIA cases, 658 local controls of self-reported non-Hispanic, European-American ancestry and 2400 “out-of-study” controls from the Molecular Genetics of Schizophrenia non-GAIN samples available from dbGAP. The demographic details are provided in Supplemental Table 1
. The cases were classified based on the ILAR revised criteria for JIA (17
). To control for potential population substructure, two principal components were identified and included in the logistic regression model as covariates. The genome-wide inflation factor was 1.04 and there was little evidence of systematic departure from expectation as shown in the Q-Q plot (Figure S1
). A total of 561,137 SNPs had MAF>0.05, no differential missingness between cases and controls, < 5% missing data and no evidence of HWE departures (P<0.01 in controls, P<0.0001 in cases). P-values of ≤ 5 × 10−8
were considered significant for genome-wide testing.
The results of the principal component adjusted analysis identified both HLA and non-HLA associations with JIA (). The strongest associations were within the MHC region and include HLA Class I and Class II loci () as well as non-HLA loci robust to linkage disequilibrium (manuscript in preparation). The strongest SNP associations outside the MHC region included novel loci as well as loci implicated in other autoimmune diseases. We have previously reported the overlap with autoimmune disease susceptibility loci in this dataset (14
). The shared autoimmune loci are labeled in blue on the Manhattan plot (). To extend association findings for JIA, ten SNPs, representing five regions not yet reported in JIA, were selected for replication testing in independent samples. These SNPs and the association results for the Discovery cohort are listed in and represent chromosomal regions 3q13, 4q31, 10q21, 15q26 and 16q12.
Figure 1 A. JIA genome-wide association results (−log10[p]) plotted on a genomic scale (Manhattan plot). The horizontal line represents a p-value of 5×10−8. This conservative threshold for genome-wide significance ignores linkage disequilibrium (more ...)
Considering the genome-wide association analysis in the Discovery cohort, evidence for association with JIA was found for the chromosome 3q13 region which includes the genes CD80, KTELC1 and C3orf1. Visual resolution of the association pattern across this region is provided in and indicates that the signal maps to a single linkage disequilibrium block. Replication was attempted for rs4688011 and rs6766899 based on the Discovery cohort statistical findings (OR=1.37, P=1.9×10−6 and OR=1.37, P=9.8×10−5, respectively). Association with JIA was consistent in the Replication cohort for rs4688011 (OR=1.18, P=0.0033) and was the strongest of all markers tested in meta-analyses (OR=1.23, P=3.6×10−7, ). Due to the disparate relative sample sizes of the cohorts, the weighted meta-analysis tends to be dominated by the UK cohort. Weighting each cohort equally modestly increased the evidence of association for both of these loci ().
Figure 2 Regional plots of JIA loci. Genotyped and imputed SNPs are plotted with their Discovery P-values (−log10[p]) as a function of genomic position (assembly hg18) within a 400 kb region surrounding the most significant SNP (purple diamond). Recombination (more ...)
Replication of JIA association findings*
Compelling evidence for genetic association was also found at 10q21, a region that has not been reported for any other autoimmune disease. Specifically, rs6479891 (OR=1.59, P=6.1×10−8) is located within the jumonji domain containing 1C (JMJD1C) gene and represents the single strongest finding outside the MHC region for the Discovery cohort. Analysis of the Replication cohort for SNPs in the 10q21.3 region revealed odds ratios that were consistent in direction, but of smaller magnitude and less significance () than the Discovery cohort. Notably, only one of the 5 SNPs, rs10995450, reached P<0.05 levels of statistical significance in the Replication cohort and noting that only the Discovery and UK cohorts were genotyped for this SNP. Weighting each cohort equally increased the evidence for association for all but rs10995450, which had comparable evidence under either approach.
An eQTL analysis was performed to assess the potential impact of the polymorphisms on expression levels of neighboring genes using two resources: 1) a dataset comprised of 68 JIA cases and 23 healthy controls from the GWAS and 2) an online dataset with GWAS and expression analysis from lymphoblastoid cells lines (LCL) collected from individuals with asthma (18
The SNPs that showed association in 3q13 were significantly associated with expression levels in the eQTL analyses (). Specifically, in the JIA/control dataset, both rs4688011 in C3orf1
and rs6766899 in CDGAP
were associated with the TMEM39A
probe set 222690_s_at (P=0.024 and P=0.0097, respectively). The order of genes on chromosome 3, and illustrated in , depicts TMEM39A
as located between CDGAP
. Interestingly, rs6766899 genotypes correlated with expression levels for CD80
(P=0.065, 1554519_at) as well. Considering the publically available online LCL dataset relative to our chromosome 3 results, we found that rs4447803, which is in LD with rs4688011 (r2
=0.54), was strongly associated with KTELC1
) with differences in expression levels comparable in magnitude to reports of KTELC1
expression in celiac disease (19
). (Note: for Celiac disease, the associated SNP is rs1599796, which is in almost complete LD with rs4688011 (r2
= 0.95 in CEU)). Neither C3orf1
expression data was available in the online LCL dataset.
Significant findings for eQTL analysis considering regions supported by genetic association*
Similarly, the SNPs that showed association with JIA in 10q21 were associated with certain expression levels (). Specifically, rs12411988 and rs6479891 were associated with expression levels measured by the JMJD1C
probe set 230007_at (P=0.008 and P=0.042, respectively; ) and other JIA-associated SNPs in the region were also associated with either JMJD1C
expression. These included SNPs in modest linkage disequilibrium (LD) with rs6479891, such as rs10761747 and rs10995450 (r2
=0.48 and 0.57, respectively). This is notable since in the online LCL dataset, both SNPs were associated with differential expression for the JMJD1C
probe set 241391_at (P=1.6×10−6
, respectively) (18
). In total, there were eight SNPs from the JMJD1C
region in the online dataset that were highly associated with expression of this gene. These SNPs were in high LD (r2
>0.77) with JIA associated SNPs and located within or near JMJD1C
which encompasses a region of over 150 kb (). Furthermore, one of the SNPs that relates to gene expression values (rs10761725) codes for a conservative substitution (Ser to Thr) in the JMJD1C protein, but is not available in the Affymetrix SNP dataset.
A third region of interest identified by genome-wide analysis was located at 4q31 and includes the gene encoding interleukin 15 (IL15). One SNP, rs4254850 showed evidence of association (OR=0.72; P=7.8×10−5), but not replication, whereas a second nearby SNP in high LD with rs425850 with more modest results in the Discovery cohort was supported by modest evidence of association in the Replication cohort (rs13139573, r2=0.84, Discovery cohort P=2.4×10−4; Replication cohort P=0.02). The regional association plot () highlights the statistical strength of the association in this region and the localization of the signal to the IL15 gene. No correlation with expression levels was found for IL15 loci in the eQTL analysis.
Replication was also attempted for two additional loci identified in the genome-wide analysis. An association with JIA was found for rs9302588 in the 16q12 region (OR=1.44, P=8.1×10−6, Discovery Cohort), but this finding lacked evidence for association in the Replication cohort. Similarly, rs12719740 on 15q26 was associated with JIA in the Discovery cohort (OR=1.45, P=6.8×10−8) but not supported by replication.
While the genetic architecture of JIA susceptibility has been informed by this and other studies, it still remains unclear how much common genetic variation accounts for the risk to JIA, as opposed to rare variants, epigenetics and gene-gene as well as gene-environment interactions. To address this, we used this GWAS data to estimate the fraction of the JIA risk that could be attributed to common SNP variation (minor allele frequency > 0.05) using a variance component threshold liability model (27
) assuming a disease prevalence ranging from 25 to 140 per 100,000 (28
). We found that common SNP variation accounts for an estimated one third of JIA susceptibility. This result holds even without the out-of study controls (). When partitioning this estimate into the individual chromosomes and extended MHC region, the extended MHC accounted for about ~8% while the remainder of the genome accounted for ~20% of the variation in JIA susceptibility.
Estimates of phenotypic variance in JIA susceptibility explained by genome-wide single nucleotide polymorphisms among unrelated individuals using a threshold model and the restricted maximum likelihood method*.