Genetic analysis of the pathogenic molecular sub-phenotype interferon-alpha identifies multiple novel loci involved in systemic lupus erythematosus.

Systemic lupus erythematosus (SLE) is a chronic autoimmune disorder characterized by inflammation of multiple organ systems and dysregulated interferon responses. SLE is both genetically and phenotypically heterogeneous, greatly reducing the power of case-control studies in SLE. Elevated circulating interferon-alpha (IFN-α) is a stable, heritable trait in SLE, which has been implicated in primary disease pathogenesis. About 40-50% of patients have high IFN-α, and high levels correspond with clinical differences. To study genetic heterogeneity in SLE, we performed a case-case study comparing patients with high vs low IFN-α in over 1550 SLE cases, including genome-wide association study and replication cohorts. In meta-analysis, the top associations in European ancestry were protein kinase, cyclic GMP-dependent, type I (PRKG1) rs7897633 (P(Meta) = 2.75 × 10(-8)) and purine nucleoside phosphorylase (PNP) rs1049564 (P(Meta) = 1.24 × 10(-7)). We also found evidence for cross-ancestral background associations with the ankyrin repeat domain 44 (ANKRD44) and pleckstrin homology domain containing, family F member 2 gene (PLEKHF2) loci. These loci have not been previously identified in case-control SLE genetic studies. Bioinformatic analyses implicated these loci functionally in dendritic cells and natural killer cells, both of which are involved in IFN-α production in SLE. As case-control studies of heterogeneous diseases reach a limit of feasibility with respect to subject number and detectable effect size, the study of informative pathogenic sub-phenotypes becomes an attractive strategy for genetic discovery in complex disease.

Introduction

Systemic lupus erythematosus (SLE) is a systemic autoimmune disorder characterized by involvement of multiple organ systems including skin, musculoskeletal, renal and hematologic systems. The pathogenesis of SLE is driven by a combination of both genetic and environmental risk factors, which lead to an irreversible break in immunologic self-tolerance [1]. SLE is four times more common in African-Americans compared to European-Americans [1], and both immunologic and genetic differences are appreciated between SLE patients from these ancestral backgrounds [2-4]. Familial aggregation and monozygotic twin studies strongly support the idea that SLE has a genetic component. There is a 50% concordance between identical twins, while first-degree relatives of SLE cases have a twenty-fold higher risk of getting SLE [5]. Genetic studies in SLE in various world populations have identified numerous susceptibility loci, however these account for far less than half of the heritability of SLE [6-11], and most of the genes described have modest overall effect sizes (odds ratio (OR) ~1.5 to 1.2) [10, 11].

Further characterizing the heritability of SLE is challenging because of the large amount of genetic and phenotypic heterogeneity. Different genetic variations and molecular pathways may be of varying importance in different patients. Previous work from our group has shown that some of the established SLE-risk loci are characterized by strong subphenotype effects, which are much greater than the overall case-control effect size [12]. This heterogeneity between patients greatly reduces the power of overall case-control studies in SLE, and is a likely explanation for much of the “missing heritability” in this disease. Designing genetic studies for SLE focusing on pathogenic molecular sub-phenotypes should greatly increase our power to detect pathogenic loci.

Interferon alpha (IFN-α) is a molecular subphenotype which is central to the pathogenesis of SLE. IFN-α is a cytokine which works at the interface of the innate and adaptive immune systems, with the potential to break self-tolerance by activating antigen-presenting cells after the uptake of self-material [13]. Serum IFN-α is elevated in many SLE patients, and levels are stable over time[14-16]. Many lines of investigation support IFN-α as a primary causal factor in human SLE[17]. We have previously demonstrated familial aggregation of high IFN-α in SLE families [5], suggesting that high IFN-α is a heritable risk factor for SLE. Additionally, recombinant human IFN-α administered to humans as a therapy for chronic viral hepatitis and malignancy can induce de novo SLE in some cases. This IFN-α-induced SLE typically resolves after the IFN-α therapy is discontinued, which supports the idea that IFN-α is causal [18, 19]. Case-control genetic studies in SLE have demonstrated remarkable over-representation of genes involved in type I interferon (IFN) signaling, production and response [11]. We have shown that many of these SLE-risk loci in the IFN-α pathway are associated with increased IFN-α pathway activity in SLE patients [20-23], supporting the idea that these loci are gain-of-function in humans. High circulating levels of IFN-α correspond to particular clinical manifestations [24], and thus activation of this pathway contributes to both susceptibility and heterogeneity in SLE [25]. We suspect that heterogeneity in the molecular pathogenesis of SLE between patients is a major factor in the unexplained heritability of the disease to date. In this study, we directly address this heterogeneity by mapping the causal IFN-α molecular trait, which allowed for detection of novel genetic variations underlying SLE disease pathogenesis. Additionally, over-activity of the IFN-α pathway has been implicated in other autoimmune diseases such as Sjogren’s syndrome and inflammatory myositis [26, 27], and it is possible that these IFN-related loci underlie some of the genetic architecture of these conditions as well.

Results

SNPs associated with IFN-α in the discovery cohort

We generated serum IFN-α data (using reporter cell assay described in Methods section to detect functional IFN-α activity) in the SLE cases who were genotyped in the SLEGEN consortium GWAS study for whom there was a serum sample available (n=400) [7]. This group was used as our discovery cohort. Re-analyzing the GWAS data in a case-case analysis stratified by high vs. low serum IFN-α, we found a number of strong associations (OR>2.0) with serum IFN-α activity. These included SNPs in the chromosome 7 open reading frame 57 (C7orf57), protein kinase, cGMP-dependent, type I (PRKG1), purine nucleoside phosphorylase (PNP), activating transcription factor 7 interacting protein (ATF7IP), and ankyrin repeat domain 44 (ANKRD44) loci (Supplemental Table 1, Figure 1a). Furthermore, a quantile-quantile (Q–Q) plot of the observed P values showed a deviation at the tail of the distribution from the null distribution (Figure 1b). We conducted a pathway analysis to identify canonical functional pathways that are enriched in the genes nearby these SNPs, and the results from this analysis are described later. This top 10 list did not share any SNPs or loci in common with the known case-control genetic associations with SLE, supporting the ability of this approach to detect novel associations. Many of the underlying genetic variations with SLE could impact particular pathways or subsets of this heterogeneous disease, and these genetic variants can be missed by large case-control SLE studies in which all patients are grouped together. We then planned to replicate all SNPs in the discovery phase with a p<10−4 (323 SNPs). In this replication list, there was one SNP which has been previously reported as a case-control association with SLE (ITGAM, rs1143678, enrichment p=0.044) [7, 28], and there were two loci on the replication list which had previously been associated with serum IFN-α levels in SLE patients (EFNA5, and ZKSCAN1/LAMTOR4, enrichment p=0.036) [29, 30] (Supplemental Table 2).

Figure 1

Top signals of association with increased serum IFN-α activity in SLE cases in the discovery phase. A) Manhattan plot shows top association signals by chromosome. B) Q–Q Plot showing association of SLE GWAS SNPs with serum IFN-α.

Validation of GWAS candidates with serum IFN-α activity in an independent cohort

The 323 top SNPs which had a p<10−4 were validated in both European-American and African-American SLE cases using an independent multi-ethnic replication cohort of 1165 SLE cases (see Supplemental Table 3 for the characteristics of the replication cohort). We used logistic regression analysis to test SNPs from the discovery cohort for association with serum IFN-α, and European- and African-American ancestral groups were analyzed separately. SNPs in the PRKG1 (rs7897633, rs7906944) and PNP (rs1049564) loci showed strong evidence for association (Table 1) in the European-American patients. In meta-analysis, both PRKG1 rs7897633 and PNP rs1049564 were associated with serum IFN-α in European ancestry with p-values that exceeded a Bonferroni correction for multiple comparisons (p<1.71×10−7, Table 1). Thus, the novel loci identified in the current study achieve genome-wide significance in the overall meta-analysis of discovery and replication sets, but none of the loci have been formally replicated. Table 2 shows a list of top SNPs in African-Americans. No significant SNP-SNP interactions were detected. Haplotype analysis was performed when evidence for association was observed for two nearby SNPs, but none of the haplotype models were superior to the individual SNP models of association. For the SNPs which demonstrated evidence for association in both European-American and African-American ancestral backgrounds, those with homogenous effects by Breslow-Day testing were analyzed in meta-analysis assuming a fixed-effect model. The two SNPs included in this cross-ancestral background meta-analysis were ANKRD44 rs4850410 [OR=0.64; 95%CI (0.48 – 0.84); PMeta =1.3 × 10−6] and pleckstrin homology domain containing, family F member 2 gene (PLEKHF2) rs297573 [OR=0.70; 95%CI (0.50 – 0.98); PMeta =1.2 × 10−4].

Table 1

List of top replicated SNPs associated with IFN-α in European-Americans

Chromosome	Locus	SNP	SNP type	Associated allele/Minor allele	Odds Ratio (95% CI)	P-discovery	P-replication	PMeta
10	PRKG1	rs7897633	intron	C	0.59 (0.44 – 0.78)	1.07E-05	2.96E-04	2.75E-08
14	PNP	rs1049564	missense	T	2.08 (1.34 – 3.21)	1.32E-05	9.88E-04	1.24E-07
6	DLL	rs1028488	intergenic	A	0.51 (0.38 – 0.70)	8.50E-04	3.12E-05	2.21E-07
4	GRXCR1	rs6850606	intergenic	A	0.64 (0.50 – 0.83)	4.75E-04	5.88E-04	1.81E-06
1	CHIA	rs7411387	intron	C	1.61 (1.24 – 2.1)	1.23E-03	3.80E-04	3.07E-06
11	TMPRSS5	rs3934007	intergenic	T	1.55 (1.19 – 2.00)	4.86E-04	9.98E-04	3.12E-06
2	ANKRD44	rs1429411	intron	C	1.55 (1.20–2.00)	9.56E-04	8.4E-04	5.04E-06

Table 2

List of top SNPs associated with serum IFN-α in African-Americans

Chromosome	Locus	SNP	SNP type	Associated allele/ Minor allele	Odds Ratio (95% CI)	P-value
10	NRG3	rs1649949	intron	C	1.60 (1.20 – 2.15)	1.37E-03
2	ANKRD44	rs4850410	intron	T	0.64 (0.48 – 0.85)	1.69E-03
5	LOC729506	rs1666793	intergenic	C	1.5 (1.10 – 2.12)	1.10E-02
8	ASPH	rs7812327	intron	T	0.66 (0.48 – 0.93)	1.59E-02
20	PLCB4	rs2299676	intron	G	0.70 (0.50 – 0.95)	2.47E-02
5	FGF18	rs7711912	near 3′	A	1.45 (1.04 – 2.02)	2.90E-02
16	RBFOX1	rs4608354	intron	A	1.57 (1.03 – 2.40)	3.44E-02
8	PLEKHF2	rs297573	near 3′	C	0.70 (0.50 – 0.98)	3.83E-02
12	KCNA5	rs526654	near 3′	G	0.75 (0.57 – 1.00)	4.00E-02

Association of GWAS candidates with autoantibody subsets in the replication cohort

Because the presence of particular autoantibodies has been strongly associated with high IFN-α in SLE [15], we also tested the SNPs which were replicated from the GWAS study for association with SLE autoantibodies. Supplemental Tables 4 and 5 show the autoantibody associations observed in different ancestral groups in replication cohort. These include some of the SNPs which were associated with IFN-α, such as PRKG1 and PLEKHF2, as well as one SNP in a locus that we’ve previously found to be associated with autoantibodies in SLE (EFNA5) [29]. None of these serological associations withstood correction for multiple comparisons. While this does not exclude the possibility of an association, we did not find strong evidence for autoantibody associations with the loci examined in this study. Each of the two SNPs which demonstrated evidence for association with both IFN-α and a particular autoantibody were examined in detail to determine the relationships between genotype and the two phenotypes. In the case of PRKG1, the IFN and serological associations were independent (Supplemental Figure 1). For the PLEKHF2 rs297573 SNP, we found that the association with IFN-α was secondary to the serological association (Supplemental Figure 2).

Canonical pathway analysis of GWAS candidate SNPs

A pathway analysis of the networks enriched among the top SNPs in the discovery cohort was generated through the use of IPA (Ingenuity Systems, www.ingenuity.com). All SNPs from the discovery cohort with p<10−4 were included. The top canonical pathways related to IFN-α associated SNPs which pass a Benjamini-Hochberg false discovery rate of 0.05 are shown in Table 3. There was prominent representation of pathways associated with neural signaling and transmission and purine metabolism, and T cell signaling was also highlighted. Some of the key molecules defining these pathways were also some of the top validated serum IFN-α associated loci in our replication cohort, such as PNP and PRKG1. Networks enriched in our study included those with various cellular functions such as cell morphology, cellular assembly and organization (PRKG1), cellular development and cell-mediated immune response (PNP), (Supplemental Table 6).

Table 3

Top canonical pathways from IFN-α associated SNPs in initial discovery GWAS data

Canonical Pathways	Ratio	P value
Axonal Guidance Signaling	0.02	5.02E-04
Synaptic Long Term Depression	0.03	4.58E-03
Xanthine and Xanthosine Salvage	0.11	7.65E-03
Dopamine-DARPP32 Feedback in cAMP Signaling	0.03	8.37E-03
Guanine and Guanosine Salvage I	0.11	1.52E-02
Adenine and Adenosine Salvage I	0.11	1.52E-02
Antiproliferative Role of TOB in T Cell Signaling	0.08	1.67E-02
Cellular Effects of Sildenafil (Viagra)	0.03	1.73E-02
Caveolar-mediated Endocytosis Signaling	0.04	1.91E-02
Cardiac Î2-adrenergic Signaling	0.03	2.01E-02

Ratio and p value are calculated as described in the Methods section.

Genome-scale Integrated Analysis of gene Networks in Tissues

Because the top loci identified in this study were not classical type I IFN pathway genes, we used the GIANT software to query potential relevance of the gene products encoded by these loci in various immune cell subsets. Figure 2 shows the networks produced by the GIANT algorithm when the top hits from our study are used as the input data in the various immune cell subset datasets available for analysis. Networks with the highest density were observed in dendritic cells and natural killer cells, and low density networks were seen in T and B lymphocytes (Table 4). Similarly, the top associations with serum IFN-α generally demonstrated the greatest network strength in plasmacytoid dendritic cells and NK cells. These data support biological relevance of the transcripts in dendritic cells, which have been implicated as the major IFN-α producing cell type in SLE [31], and NK cells, which have been reported to play a critical cooperative role with dendritic cells in the production of IFN-α [32]. Additionally, when examining the other molecules functionally implicated in these networks, a number of SLE-associated molecules are observed in the network diagrams, including IL12, TLR7, and the JAK/STAT pathways.

Figure 2

Tissue specific analysis of gene networks in different immune cells. Networks demonstrate relationships between PNP, PRKG1, ANKRD44 and PLEKHF2 to other molecules in immune cells. Edges with weight (relative confidence) greater than 0.4 are shown. Each network diagram represents a different immune cell type as follows: A: B lymphocyte, B: Dendritic cell, C: Monocyte, D: Neutrophil, E: NK cell, F: T lymphocyte.

Table 4

Network density and network strength analysis for tissue specific gene networks in different immune cells

Cells	Network Density	Network Strength
		ANKRD44	PNP	PRKG1	PLEKHF2
B lymphocyte	0.12		1.1		5.5
Dendritic cell	0.58			9.3	27.4
Monocyte	0.18	0.4	3.7		3.4
Neutrophil	0.13	1.4	1.9		8.7
T lymphocyte	0.09	0.4	1.4		7.3
NK cell	0.52	11.7	1.4	17.0	24.6

Networks generated by the GIANT software program for each immune cell type. Network density and strength calculated as described in the Methods. Density is calculated for the overall network in the cell, and strength is calculated for each of the loci entered in the analysis.

Discussion

SLE is a highly heterogeneous disease, hence it is likely that certain genetic factors will be related to particular disease phenotypes and pathogenic pathways [1, 12, 33], and that genetic associations will not be shared between all SLE patients. We suspect this is a major factor in the unexplained heritability of the disease to date. Directly studying this heterogeneity by mapping a causal molecular trait in this study greatly enhanced our power to detect novel genetic variations underlying SLE disease pathogenesis. The top loci in our study have not been previously reported in other case-control studies of SLE, and were not top loci in the initial case-control analysis of the GWAS data set we used in this study [7]. Thus, our alternative strategy was capable of finding genetic variants associated with disease that are not readily apparent in case-control designs, supporting a complexity in genetic architecture that will require molecular subphenotyping to fully delineate.

PRKG1 (rs7897633, p=2.75×10−8) was the strongest association observed in our study. This gene codes for the soluble isoforms of the cyclic GMP-dependent protein kinase (Iα and Iβ), which are important components of signal transduction processes in diverse cell types [34]. Canonical pathway analysis revealed this gene was associated with pathways such as synaptic long term depression, Dopamine-DARPP32 Feedback in cAMP Signaling, and netrin signaling; pathways in which modulation of cGMP and cAMP plays an important role in signaling and function. GIANT analysis supported biological function for PRKG1 in dendritic cells and NK cells, two cell types which cooperate to generate IFN-α in SLE [32]. PRKG1 function was not as strongly supported in other immune cells such as T and B lymphocytes, which are not thought to be major IFN-α producing cells. It is not immediately clear how PRKG1 might impact type I IFN production or signaling, but mechanistic experiments directed at the IFN pathway in both dendritic and NK cells are likely to be informative.

The other strongly associated polymorphism in our study was a missense SNP in the PNP gene (rs1049564, p=1.24×10−7). PNP encodes the nucleoside phosphorylase enzyme, which is involved in purine metabolism. PNP together with adenosine deaminase (ADA), serve a key role in purine catabolism in the salvage pathway. Deficiency in this pathway leads to build up of elevated deoxy-GTP levels, resulting in T-cell toxicity and deficiency [35, 36].

Rare autosomal deficiency of the PNP gene results in a metabolic disorder characterized by defective T-cell and B-cell immunity as well as defective antibody responses [37, 38]. Interestingly, PNP-deficient patients have also been reported to develop autoimmune disorders, such as SLE, autoimmune hemolytic anemia, and idiopathic thrombocytopenic purpura [39]. The SNP we find associated is a common coding-change variant which does not cause complete deficiency, and whether this variant results in some change in enzyme function is not currently known. In silico bioinformatic analysis using PolyPhen2 and SIFT predicts this SNP as non-damaging, but an effect on enzyme activity would still be possible. There was strong representation of the purine metabolic pathway in our canonical pathway analysis, and PNP was the key molecule associated with this pathway. Some rare, highly penetrant variants in genes involved in nucleic acid metabolism have been associated with SLE, such as TREX1 and DNASE1L3 [40-42]. Given this precedent, PNP is a fascinating genetic association with SLE.

The non-European ancestral backgrounds studied were smaller, and did not allow for strong independent significance. Our discovery set was exclusively of European ancestry, and thus variants specific to other ancestral backgrounds could not have been discovered. Despite these limitations, we observed some interesting evidence for associations which were of similar effect in both European- and African-American ancestral backgrounds. Both ANKRD44 (rs4850410; PMeta =1.3 × 10−6 ) and PLEKHF2 (rs297573; PMeta =1.2 × 10−4 ) were associated with IFN-α in African-American and European ancestral backgrounds. PLEKHF2 is an endosome-associated protein responsible for modulating the structure and function of endosomes, as well as the endocytotic process [43]. PLEKHF2 can increase the activity of caspase 12, and a role for in ER-related apoptotic pathway has been suggested [44]. ENCODE ChIP-seq data demonstrate that rs297573 SNP downstream of PLEKHF2 which was associated with IFN-α in our study resides in NFκB transcription factor binding site. ANKRD44 has not been extensively studied, but it binds to the catalytic subunit of protein phosphatase 6 [45], which plays a role in cell cycle progression.

Our initial discovery cohort showed association of two SNPs in the C7orf57 locus in chromosome 7 with serum IFN-α activity; however, this locus failed to replicate. One of the possible reasons for lack of replication could be that this locus was related to some unique feature of the discovery cohort that was not present in the replication cohort. One previous GWAS study of amyotrophic lateral sclerosis which used a very similar Illumina genotyping platform found evidence for association between these two SNPs and ALS which then failed to replicate in an additional independent replication cohort [46]. It is possible that some peculiarity of the earlier Illumina genotyping platform made it more likely for these SNPs to be spuriously associated, although this locus was not associated in the original SLEGEN GWAS case-control study [7]. We used an entirely different rtPCR-based genotyping method for our replication cohort to eliminate potential platform-related biases.

As referenced above, discovery methods followed by replication in non-European ancestral backgrounds would be an important next step to this work. It is likely that some polymorphisms will be ancestry-specific, and will not be evident until a discovery strategy is used in that particular ancestral background. This would be especially important for African-Americans who have a higher incidence of SLE and more severe clinical manifestations [1]. African-American SLE cases also have higher levels of serum IFN-α activity [15], which could be one factor related to the increased incidence and severity of the disease. Our findings could have pharmacogenomic implications, as therapeutics targeting the IFN-α pathway are currently in development for SLE. Knowledge of the functional genetic factors underlying IFN-α dysregulation in a given patient could be useful in individualizing therapy with these agents.

Materials and Methods

Samples and Genotyping

Discovery cohort

Genome wide association study (GWAS) data from 755 SLE cases were obtained from multiple study centers as part of the international consortium for Systemic Lupus Erythematosus Genetics (SLEGEN) [7]. The cohort studied by the SLEGEN consortium for GWAS in SLE consists of unrelated women of self-reported European ancestry and has been described in detail [7]. Out of 755 SLE cases; 400 cases had serum available for IFN-α analysis and were included in the discovery GWAS phase. Samples were genotyped at 317,000 SNPs on the Illumina Infinium HumanHap300 genotyping Beadchip (Illumina Inc., San Diego, CA, USA). SNPs that failed the Hardy-Weinberg equilibrium test (p<0.001) were excluded, as were SNPs with a genotyping success rate less than 95% or with a minor allele frequency less than 0.05, resulting in 291,943 SNPs that were used in the analysis.

Replication cohort

The independent multi-ethnic replication cohort of 1165 SLE patients was obtained from the Lupus Family Registry and Repository (LFRR) at the Oklahoma Medical Research Foundation and consisted of the following self-reported ancestral backgrounds: 715 European-Americans and 450 African- Americans. We incorporated 238 Hispanic/Native-American, and 40 Asian-American SLE cases in addition to 1165 SLE cases in the principal components analysis to determine population stratification in the replication cohort. Association analyses were not performed in the Hispanic-American/Native American and Asian-American cases due to the small number of subjects. Clinical characteristics and demographic details for the patients in the replication cohort are summarized in Supplemental table 3. Informed consent was obtained from all patients in both cohorts included in this study, and the study was approved by the institutional review boards at the respective institutions. We followed up IFN-α associated SNPs which had a p-value of 1×10−4 or less from the initial discovery GWAS analysis. SNPs that failed SNP assay design were excluded, resulting in 323 SNPs which were genotyped in the replication cohort. A separate panel of 334 ancestry-informative markers (AIMs) [47] was also genotyped in the replication cohort. SNPs were genotyped in genomic DNA using the Fluidigm Biomark microfluidic qPCR system. All DNA samples were pre-amplified using the SNPtype™ primers from the genotyping assays, according to the manufacturer’s protocol. PCR data were analyzed using the BioMark SNP Genotyping Analysis software version 3 to obtain genotype calls. Scatter plots were all reviewed individually for quality, and SNPs that deviated significantly from the expected Hardy-Weinberg proportions (P < 0.001) or with less than 95% genotyping success were excluded from the analysis.

Measurement of serum IFN-α activity

ELISA methods for the measurement of type I IFN in human sera have been complicated by low sensitivity and low specificity [48]. We used a well-documented sensitive and reproducible reporter cell assay was used to generate IFN-α activity data from patient sera [5, 49]. The reporter cells in this assay (WISH cells, ATCC #CCL-25) measure the ability of patient sera to cause IFN-induced gene expression. These cells are an epithelial-derived cell line that is highly responsive to IFN-α. Cells are incubated with patient serum for 6 hours. Then real-time PCR is used to quantify three canonical IFN-α-induced transcripts in the WISH cell lysates (IFIT1, MX1 and PKR). Sera from healthy unrelated controls (n=200) were tested to establish a normal value for the assay. Results from patient samples are expressed as the number of standard deviation (SD) above the mean of healthy unrelated control sera. The sum of the number of SD above healthy controls for the three transcripts is used as the quantitative output from the assay, representing a serum IFN-α activity score. This assay has been extremely informative in SLE and other autoimmune diseases[15, 27, 50].

Statistical Analysis

Control for Population Structure

To account for potential differences in admixture or population structure within self-reported ancestral backgrounds in the discovery and replication cohort, we performed a principal component analysis (PCA) using the GWAS SNPs and 334 independent AIMs SNPs respectively. PCA in the discovery cohort was carried out on all of the GWAS SNPs that passed quality control thresholds. This cohort is composed of SLE cases with self-reported European ancestry, and as shown in the principal components analysis plots (Figure 3A), there were no major population outliers. As expected, cases with varying proportions of Northern and Southern European ancestry were included in the study, and some cases cluster with the Ashkenazi Hap Map reference population, suggesting Jewish ancestry (Figure 3B). PCA of the AIMs genotyped in the replication cohort revealed that the principal component (PC) 1 obtained in this analysis provided a strong separation between subjects of self-reported African-American ancestry and the non-African ancestral backgrounds, while PCs 2 and 3 provided a separation between subjects of self-reported Asian-, Hispanic-, and European-American ancestry (Figure 3C and 3D). Self-reported Hispanic-American/Native-American ancestry subjects were largely overlapping in this analysis, and are considered together in this analysis. Association analyses were not performed in the Hispanic-American/Native American (n=238) and Asian-American cases (n=40), due to the small number of subjects. These subjects were included in the principal component analysis of the AIMs to assist with the determination of population structure. Correction for population structure within the discovery and replication cohort was done using the first three PCs as covariates in the logistic regression association analyses. PCA analysis was performed using Cluster 3.0 software [51].

Figure 3

Principal component analyses to detect population structure. A. and B. show principal components from all SNPs studied in the SLEGEN GWAS data set. All studied subjects are included, and reference populations from HapMap 3 samples are also included. Each circle represents an individual sample. PC = principal component, SLE=SLEGEN samples, CEU=Utah residents with Northern and Western European ancestry from the CEPH collection, CHB=Han Chinese in Beijing, China, JPT=Japanese in Tokyo, Japan, TSI=Toscani in Italy, YRI=Yoruba in Ibadan, Nigeria, AJ=Ashkenazi Jewish. C. and D. show the principal components derived from the AIMs in the replication cohort. Each symbol represents an individual sample, and colors represent self-reported ancestry.

Association Analyses

Logistic regression analysis was used to detect associations between the SNPs and serum IFN-α in both stages of the study. IFN-α activity was studied as a categorical trait because the trait distribution is highly skewed, such that log transformation does not result in a normal distribution and the highly skewed data did not allow for linear modeling in a QTL analysis. We used a binning strategy that has been highly informative in previous large scale studies and multivariate analyses of the serum IFN trait in SLE [4, 15, 52], in which subjects with a value greater than 2 standard deviations above the mean of healthy controls are binned as high IFN-α, and the rest are binned as low IFN-α. Using this binning strategy prevents high outlying values from exerting an inordinate amount of influence in the model. In the discovery cohort, 88 were categorized as high IFN-α, and 322 were categorized as low IFN-α. Logistic regression analysis was carried out using PLINK v.1.07 software [53]. The first three principal components from the PCA of the GWAS SNPs were used as covariates in the logistic regression to control for population structure in the discovery cohort. In the replication cohort, each self-reported ancestral background was analyzed separately, and the first three principal components were included as covariates to correct for population structure and admixture. Regression analysis was also performed to detect any potential associations between the presence of particular autoantibodies and SNPs in the replication cohort, because autoantibodies have been associated with high IFN-α in SLE patients [15]. In the replication cohort, we used the Benjamini-Hochberg procedure to control the false discovery rate at 0.05, and the SNPs which passed this threshold were considered for meta-analysis. The p-value threshold used for significance in the overall meta-analysis corrects for the number of SNPs which were analyzed for association in the initial GWAS discovery analysis, controlling the family-wise error rate at the 0.05 level. For SNPs that demonstrated a homogenous effect across the discovery and replication sets by Breslow-Day testing, meta-analysis was performed using the weighted Z-score method [54] using R 2.11.1 statistical analysis software (www.r-project.org). For statistical correction for multiple comparisons, we applied a Bonferroni correction to the meta-analysis results using the number of SNPs that passed quality control in the discovery GWAS (n=291,943), resulting in a threshold p-value for this study of P<1.71 × 10−7. In the cross-ancestral background analysis, SNPs that demonstrated a homogenous effect across both ancestral backgrounds were meta-analyzed using the same weighted Z-score method [54], assuming a fixed effect model. Enrichment p-values were calculated using a Fisher’s exact test with the following parameters: for the SNP-wise calculation for the SLE-associated ITGAM SNP, the number of possible confirmed SLE-risk SNPs in European ancestry was estimated at 40, and the number of SNPs that passed quality control in our GWAS screen was used as the denominator to establish the null proportion. The observed proportion was one SNP out of the 323 SNPs in our replication list. For the locus-wise comparison for genes associated with circulating IFN in SLE, we estimated 18 loci which have been previously associated with IFN in SLE, and 20,000 as the number of human gene loci to estimate the null proportion. The observed proportion was 2 loci out of the 277 loci represented by the 323 SNPs in the replication list.

Canonical pathway analysis

From the initial discovery GWAS data, IFN-α associated SNPs (n=323) with a p-value of 1×10−4 or less were analyzed further using Ingenuity Pathway Analysis (IPA) (Ingenuity® Systems, www.ingenuity.com) to identify the top canonical pathways related to IFN-α associated SNPs. SNPs were attributed to the nearby gene, and the genes were then compared to curated functional attribution lists organized by canonical pathway function. The magnitude of over-representation of a particular canonical pathway in the gene list from our study was calculated as the ratio of the number of molecules from our data set that map to the pathway divided by the total number of reference molecules in that pathway in the IPA database (a list of genes belonging to major cannoncial pathways is curated in IPA based on published literature). Statistical significance was determined using the Fisher’s exact test, comparing the observed ratio of genes in a particular pathway to the null expectation (that the genes would assort proportionally across all IPA pathways), to estimate the probability that the observed over-representation of the particular pathway would arise by chance.

Genome-scale Integrated Analysis of gene Networks in Tissues (GIANT)

The top genes from the replication cohort were queried using the GIANT software program to determine likely functional relationships of these genes in various types of immune cells. GIANT is a public, web-based software program that uses tissue-specific gene expression databases to predict tissue-specific gene interactions (http://giant.princeton.edu/about/). 145 tissues/cell-types are available to be queried, including major immune cell subsets. The software generates functional networks based on the genes queried via the integration of thousands of publicly available gene expression datasets, sequence data, transcription factor binding sites, and protein-protein interaction data to generate gene-association matrix. Bayesian weights derived from the gold-standard tissue-specific data sets are then applied, and networks are generated for each tissue queried which illustrate the most probable functional relationships between the queried genes and other molecules in that particular tissue. Network relationship confidence (edge weight) was set at a minimum of 0.4 for our analyses. After the networks were generated, we calculated overall network density and network strength of each of our study genes in each immune cell subset network. Network density (D) was defined as a ratio of the number of edges (E) to the number of possible edges, given by the binomial coefficient , giving D=2E/N (N-1); where N=number of nodes. In these weighted networks, we calculated strength as the sum of a node’s edge weights.

Prediction of the impact of coding-change SNPs

Prediction of consequences on protein structure and/or function of non-synonymous single nucleotide polymorphisms were evaluated using the prediction programs SIFT (Sorting Intolerant From Tolerant, http://blocks.fhcrc.org/sift/SIFT.html) and PolyPhen (Polymorphism Phenotyping, http://genetics.bwh.harvard.edu/pph/). These two programs use algorithms to determine the likelihood that a particular coding-change polymorphism impacts protein-folding based upon local protein structure as well as the particular amino acid substitution.