Identification of Novel Type 1 Diabetes Candidate Genes by Integrating Genome-Wide Association Data, Protein-Protein Interactions, and Human Pancreatic Islet Gene Expression

A total of 395 positional candidate genes were identified from non-HLA type 1 diabetes–associated linkage disequilibrium (LD) regions from GWAS. The LD intervals were calculated based on the HapMap CEU founders data in snpMatrix (http://www.bioconductor.org/packages/release/bioc/html/snpMatrix.html) using different D´ and r² thresholds, modified from (1,2). The positional candidate genes were used as input into the STRING database (http://string-db.org/) to build protein networks (Fig. 1). In STRING, the functional associations are derived from four sources: genomic context, high-throughput experiments, conserved coexpression, and text mining, while the physical interactions are derived from only the high-throughput experiments. Protein networks were built running the Markov CLustering (MCL) algorithm using the functional and physical interaction scores from STRING with an inflation parameter of 4.5 and forcing all the nodes to be connected. MCL is a fast, scalable algorithm for unsupervised clustering of graphs and networks, which has been used by many groups to build protein complexes or families based on protein-protein interactions (20,21). To build meaningful networks, only the protein networks that contained input proteins encoded by genes from at least two different type 1 diabetes–associated loci were considered. We obtained 11 networks of proteins that interact functionally (networks 1–11) and 6 networks of proteins that interact physically (networks 12–17). These networks were produced in STRING and recolored in Inkscape (www.inkscape.org). The edges that connected proteins that physically interact were highlighted in the networks that were differentially regulated (Figs. 2 and 3). All of the 235 proteins contained in the resulting 17 networks (Figs. 2 and 3 and Supplementary Fig. 1) were, if possible, analyzed for gene expression levels in human pancreatic islets.

Human pancreatic islet preparations were obtained from a multicenter European Union–supported program on β-cell transplantation in diabetes. The program has been approved by central and local ethics committees. Human islets were obtained from eight nondiabetic donors (aged 8–57 years; six males and two females). None had classical type 1 diabetes–associated HLA-DR risk genotypes (DR1/3, DR11/52, DR4/5, DR2/5, DR6/13, DR4/6, DR4/5, and one donor missing HLA-DR genotype). Islet preparation, cytokine stimulation, including tumor necrosis factor-α (TNF-α; 5,000 units/mL), interferon-γ (IFN-γ; 750 units/mL), and interleukin (IL)-1β (75 units/mL) for 48 h, and RNA extraction have been described previously (7,22,23).

Relative expression levels of selected genes were evaluated by use of TaqMan assays. The Low Density Array system (Applied Biosystems) containing assays for the individual genes as well as housekeeping genes was used on TaqMan 7900HT (Applied Biosystems). Of the 235 genes contained in the 17 protein-protein interaction networks, 222 could be examined by a TaqMan gene expression assay. Expression levels of genes were normalized against the geometric mean of three different housekeeping genes, GAPDH, 18S-RNA, and PPIA, and evaluated using the ΔΔCt (cycle threshold) method (24). The relative expression levels for unstimulated versus cytokine-stimulated islet preparations were compared using paired t tests. P < 0.05 was considered statistically significant. Transcripts with Ct values <38 were considered to be expressed.

INS-1 cells were maintained in RPMI 1640 medium (11 mmol/L glucose) supplemented with 10% FBS, 100 units/mL penicillin, and 100 μg/mL streptomycin (all from Invitrogen). In addition, the media contained 50 μmol/L β-mercaptoethanol. For mRNA purification, 100,000 cells were seeded in quadruplicates per condition in 48-well dishes. mRNA was extracted by RNeasy kit according to the manufacturer’s protocol (QIAGEN). cDNA was prepared from total RNA as described by the manufacturer (Applied Biosystems). Relative expression levels of target genes (Cd83, Il17rd, Il27ra, Cxcr7, Myo1b, Plcg2, Ifngr1, and Traf3ip2), the input gene (Ccr7), and the housekeeping genes (Ppia and Gapdh) were evaluated by the use of commercially available rat TaqMan gene expression assays (Applied Biosystems). Expression levels of genes were normalized against the geometric mean of Ppia and Gapdh and evaluated using the ΔΔCt method.

To evaluate whether the networks were statistically enriched for cytokine-regulated genes, we compared the number of regulated genes after false discovery rate (FDR) correction within each network with the total number of genes that were regulated by cytokines on the Affymetrix Human Genome U133 Plus 2.0 array as a reference. For these experiments, total RNA from four human islet preparations provided through the Juvenile Diabetes Research Foundation (JDRF) Islet Distribution Program (JDRF award 6-2005-1178) was used. The islets were treated with IL-1β (1 ng/mL), IFN-γ (20 ng/mL), and TNF-α (8 ng/mL) for 48 h before RNA extraction. Although not identical to the conditions used for transcript profiling of network genes, we consider this design to be adequate for the (cytokine-induced) gene enrichment analysis. The gene expression was normalized using the robust multiarray analysis method, and probes were annotated using an updated probe set definition (25). Of the 17,491 genes analyzed on the array, 154 were significantly regulated (P < 0.05, adjusted for multiple testing by FDR) (26). Enrichment scores for significantly regulated genes within the networks compared with the Affymetrix microarray were calculated by Fisher exact test.

To evaluate whether the networks included noninput genes that harbored type 1 diabetes–associated SNPs, we used the meta-analysis GWAS data generated by the Type 1 Diabetes Genetics Consortium (T1DGC) (1). The SNPs in the GWAS dataset were first mapped to Ensembl gene IDs (mapping based on National Center for Biotechnology Information build 36) and then mapped to proteins within the 17 networks. We included only the minimal promoter and 3′-untranslated region (i.e., 2 kilobases [kb] up- and downstream, respectively, from the coding region) to reduce overlap between genes and, thus, SNPs that mapped onto more than one gene. SNPs located >2 kb from a gene or in gene deserts were not considered in the present analysis. Networks 8 and 9 did not contain any genes outside the type 1 diabetes–associated regions, and these networks were therefore not analyzed further.

To assess whether the individual networks were significantly enriched for specific GO terms or KEGG pathways, we used FatiGO (27). For each of the networks 1 to 11, which are based on functional interactions, the set of member proteins was compared with the entire set of proteins from those networks (in total, 147 proteins). Likewise, every set of proteins from networks 12 to 17 was compared with the whole set of proteins that belong to the networks based on only experimental interactions (in total, 88 proteins).

To evaluate whether the number of type 1 diabetes input proteins observed in the networks was higher than expected by chance, we generated 100 random networks with the same size and topological characteristics for each of the 17 original networks. The number of input proteins observed in the random networks was compared with the number in the original networks, and if an equal or higher number of input proteins was observed in <5% of the random networks, we considered the original network to have an excess of type 1 diabetes input proteins.

To evaluate the specificity of the regulations that we observed in pancreatic islets, we searched the Gene Expression Omnibus database for human gene expression datasets suitable for comparison with our expression results. The approach was to search for datasets that were related to the treatment conditions of our study (i.e., inflamed or cytokine-exposed tissues/cells) but in tissues/cells that are unrelated to the target tissue in diabetes. We identified the datasets GDS3158 (rheumatoid arthritis cartilage), GDS3011 (IL-1–exposed keratinocytes), GDS2939 (cytokine-exposed thyroidea cells), GDS2516 (IFN-exposed endothelial cells), and GDS2478 (inflammatory acne), available from http://www.ncbi.nlm.nih.gov/gds, and we compared how many of the genes within the significantly regulated networks were regulated with a P value of <0.1 in the five datasets.

All statistical analyses were performed in R (www.r-project.org).

In this study, summarized in Fig. 1, we used all 395 positional candidate genes within the LD intervals of 44 non-HLA type 1 diabetes–associated regions identified from GWAS (1) as input genes (Table 1 and Supplementary Table 1) to construct biological protein networks using the STRING database (18). This database includes both direct (physical) and indirect (functional) protein interactions. Our hypothesis was that minor susceptibility genes cooperate in networks to confer type 1 diabetes risk. Hence, two or more minor susceptibility genes functioning in the same network or pathway may synergize to affect phenotype. Therefore, when constructing the networks, only networks with at least two input genes were considered. The HLA region genes (177 genes in the LD interval of the major histocompatibility complex [MHC] region, position: 29.91–33.50 Mb on chromosome 6) were not included as input genes because the strong LD and genetic contribution from this region would likely mask all other effects. We also did not allow HLA/MHC region genes as secondary interaction partners in the constructed networks. By this approach, we identified 17 protein networks (Figs. 2 and 3 and Supplementary Fig. 1) containing in total 235 nodes, of which 41 were input genes (i.e., encoded by genes in type 1 diabetes–associated GWAS regions) (Table 2 and Supplementary Table 2).

Gene expression may represent an intermediate phenotype between the genetic variation and the disease process. Therefore, we investigated whether cytokines affected the identified type 1 diabetes protein networks at the expression level. We evaluated mRNA expression by real-time PCR of 222 of the 235 genes in the 17 networks between untreated (control) and cytokine-exposed (disease) human pancreatic islets from eight donors. A total of 27 genes were differentially regulated after cytokine exposure of the islets (P < 0.05) (Table 2 and Supplementary Table 2). We then evaluated whether any of the 17 networks was statistically enriched for cytokine-regulated genes. For this, we used Affymetrix microarray data from human islets stimulated with cytokines to get a value of global cytokine-induced gene regulation. A total of 154 of the 17,491 genes that were interrogated on the Affymetrix array were significantly regulated by cytokines (corrected P < 0.05). The number of regulated genes (FDR corrected) in each of the 17 networks was then compared with this array background distribution using Fisher exact test. Three networks (networks 3, 4, and 16) were significantly enriched for cytokine-regulated genes (P values: 3.262 × 10⁻⁶, 0.006694, and 0.0449, respectively, see 2research design and methods) (Table 2). Network 3 included 11 differentially regulated genes, whereas cytokines modulated the expression of 4 genes in network 4 and 7 genes in network 16 (Table 2). In addition to the enrichment of cytokine-regulated genes in the networks, our expressional analysis also confirmed basal expression, defined as Ct values <38, of the majority of the genes in the 3 networks in human islets, highlighting that all of these networks can be functional in islets. For network 3, the fraction of expressed genes was 60% (26 of 43); for network 4, it was 78% (14 of 18); and for network 16, it was 86% (38 of 44) (Table 2).

To identify nominally associated SNPs (P < 0.05) in network genes, we used the complete genetic information from the GWAS meta-analysis (1). SNPs and their corresponding P values were mapped to all the genes in the 17 networks except for the genes used as input. Networks 3, 4, and 16, which were enriched for cytokine-regulated genes, also contained several genes with nominally associated SNPs (Figs. 2 and 3 and Table 2). It is interesting that four of the genes (IL17RD, CD83, IFNGR1, and TRAF3IP2) in network 3 that were regulated by cytokines also contained nominally associated SNPs. In a similar manner, IL27RA in network 4 was both differentially regulated and harbored an associated SNP. In network 16, this was true for three genes (PLCG2, MYO1B, and CXCR7). None of these eight genes have previously been identified as candidates to be involved in type 1 diabetes.

These new candidate genes have often been associated with roles in immune cells, but our results suggest that they are also expressed in human islets and potentially in the β-cells. Therefore, we examined the expression of the novel candidate genes in the rat β-cell line INS-1. We found that all the genes, with the exception of Il27ra,were both expressed and regulated by cytokines in INS-1 cells (Fig. 4). Both direction and fold changes for these genes were comparable with our findings in human islets (Table 2). We also evaluated the expression of Ccr7 in INS-1 cells because this gene was the only input gene in these networks that changed expression upon cytokine treatment in human islets. We found Ccr7 also to be expressed and significantly regulated by cytokines in the rat β-cell line (3.95 ± 1.70-fold induction by IL-1β compared with untreated control, P < 0.05, n = 4). Together, these findings verify that IL17RD, CD83, IFNGR1, TRAF3IP2, PLCG2, MYO1B, and CXCR7 are expressed in β-cells and therefore indicate that they may play a functional role in these cells (Table 3).

Type 1 diabetes is caused by an immune-mediated destruction of the insulin-producing β-cells, and current interpretation of GWAS data supports the involvement of classical immunoregulatory pathways, such as modulation of the IL-2 pathway, cytokine signaling, and changes in subsets of T cells (1,28). However, it is obvious from most models of type 1 diabetes that this immunodysregulation does not per se lead to type 1 diabetes development unless initial β-cell damage has taken place (2). Knowledge of the β-cell and its role in the pathogenesis is therefore essential in the fight against diabetes. Proinflammatory cytokines, such as IL-1β, IFN-γ, and TNF-α, have been proposed to be mediators of β-cell dysfunction in type 1 diabetes (19,29).

In the current study, we have constructed protein interaction networks from all genes located in non-HLA loci associated with type 1 diabetes. The protein network approach is considered more flexible and powerful than sorting genes into classical pathways or functional categories because the results of these methods might be limited to a priori knowledge (17,30,31).

We applied transcript profiling of the protein networks to score the importance of the networks that were enriched for genes affected by proinflammatory cytokines in human islets. In addition to the identification of three significantly regulated networks, we demonstrated that the majority of the genes within the networks were expressed in islets. The high percentage of confirmed islet expression strengthens the idea of using protein networks to add a biological context to candidate genes identified through GWAS.

The three highlighted networks included eight regulated genes that contained SNPs that were significantly associated with type 1 diabetes after removing the input genes. Seven of these were also expressed and regulated by cytokines specifically in β-cells. The overlap between the genetic evidence and expressional support from studies of human pancreatic islets and INS-1 β-cells indicates a prominent role for the networks in β-cells for type 1 diabetes pathogenesis. We examined if any of the most significant SNPs would give rise to changes in coding sequences, but all SNPs were found to be silent. However, the SNP in TRAF3IP2 is in perfect LD (D′ = 1, r² = 1) with rs1043730, which changes the amino acid at position 323 from His to Gln. Also, the SNP in the candidate IL17RD is in moderate LD (D′ = 1, r² = 0.56) with rs6780995, which leads to an amino acid substitution at position 255 from Thr to Met. The exact effects of these amino acid substitutions in TRAF3IP2 and IL17RD are difficult to predict, but they both are predicted to be benign (PolyPhen-2, http://genetics.bwh.harvard.edu/pph2/). It is interesting that TRAF3IP2 is suggested to be a key link between IL-17–mediated immune response and nuclear factor-κB (NF-κB)–controlled innate immunity (32). TRAF3IP2 is recruited to the IL-17 receptor after ligand binding with subsequent signaling via TRAF3 and TRAF6 and activation of the NF-κB and mitogen-activated protein kinase pathways. Coding variants in the TRAF3IP2 gene have been associated with psoriasis in recent GWAS and also with altered binding of TRAF3IP2 to TRAF6 and likely subsequent altered activation of the NF-κB pathway (33). Such alterations would have phenotypic effects in the β-cells and, thus, would potentially affect β-cell function and/or the destruction process in type 1 diabetes. Future studies should be designed to shed more light on this.

In network 3, 11 genes were differentially regulated (Table 2). Of particular interest, disruption of the input gene CCR7, which encodes a chemokine receptor and was upregulated by cytokines in islets and β-cells, modulates both streptozotocin (STZ)-induced and spontaneous diabetes in the nonobese diabetic (NOD) mouse (34,35). Our data support the hypothesis that chemokine (C-X-C motif) ligand (CXCL)-10, and also CXCL9 and CXCL11, which were all strongly induced by cytokines, accelerates the autoimmune process by enhancing the attraction of antigen-specific lymphocytes to the pancreatic islets and by directly causing β-cell damage (36). Also CD83, which harbored SNPs nominally associated with type 1 diabetes, is involved in antigen presentation and T-cell activation (37). The cytokine IL-17 promotes diabetes development in NOD mice (38,39) and potentiates human islet apoptosis in vitro induced by IL-1β plus IFN-γ and IFN-γ plus TNF-α (40). Of note, CD4 T cells in new-onset type 1 diabetic patients secrete IL-17 (40). Our findings emphasize that several proteins of the IL-17 signaling pathway are expressed in islets, regulated by cytokines, and/or harbor disease-associated SNPs and, thus, support a central role for IL-17 in β-cell destruction in type 1 diabetes. This is further supported by a recent study demonstrating that reversal of type 1 diabetes in NOD mice with anti-CD3 and IL-1 blockade is accompanied by a reduced intrapancreatic expression of IL-17 (41). Thus, this network highlights a prominent role of islet-expressed chemokines and cytokines for the pathogenesis. Consistent with this, analyses of KEGG and GO terms revealed that network 3 harbored several significant terms, including “cytokine-cytokine receptor interaction” (KEGG), “immune response” (GO), “nucleus” (GO), and “receptor binding” (GO) (Supplementary Table 3).

In network 4, the gene programmed cell death 1 ligand 2 (PDCD1LG2) plays an essential role in immune tolerance (42). Neutralization of IL-15 attenuates diabetes in mice induced by multiple low-dose STZ (43). The expression of IL-27 receptor A (IL27RA) was suppressed by cytokines in human islets, but Il27ra was not expressed in INS-1 β-cells. The ligand for IL27RA, IL-27, potentiates subdiabetogenic STZ treatments by inducing β-cell damage and diabetes (44). This network also contained IL2RA, which was not regulated by cytokines in pancreatic islets. However, haplotypes of IL2RA are associated with type 1 diabetes (45), and the expression level of IL2RA on CD4⁺ memory T cells correlates with these haplotypes (46). In combination with significant GO terms, such as “cell activation” and “leukocyte activation” (Supplementary Table 3) for network 4, these findings support a role of this network and ILs for inflammation and type 1 diabetes.

Regarding candidate genes from network 16, studies have suggested that CXCR7 and PLCG2 play a role in immune-mediated destruction of β-cells (47,48). The two input genes, KCNH7 and SH2B3, interacted with GRB2, which was highly expressed in islets (data not shown) and is central in the network. Of interest, GRB2 is involved in the insulin signaling pathway through interaction with insulin receptor substrate 1 (49)—a pathway that is also essential in β-cells for proper insulin biosynthesis and release (50).

For 14 of the 17 networks, we observed a significant excess of input proteins (P < 0.05) compared with randomly generated networks with the same size and topological characteristics. Two of the remaining 3 networks (networks 1 and 6), were not enriched for cytokine-regulated genes. Network 16 was also not significantly enriched for input proteins, which is not surprising based on the low number of input proteins (two) compared with the large size of the network. However, this network is entirely based on physical interactions and contains a significant excess of differentially regulated genes in the network, which validates a possible functional role in type 1 diabetes pathogenesis.

To assess whether spurious signals accounted for the observed transcriptional profiles, we compared the expression data of the cytokine-regulated networks with other human gene expression datasets from inflamed or cytokine-exposed tissue of unknown relevance for type 1 diabetes, such as rheumatoid arthritis cartilage, IL-1–exposed keratinocytes, cytokine-exposed thyroidea cells, IFN-exposed endothelial cells, and inflammatory acne available from Gene Expression Omnibus Data Sets (http://www.ncbi.nlm.nih.gov/gds). Although there was limited overlap between these datasets, the comparison clearly showed that there are other conditions that closely resemble the inflammatory process that we mimic with the experimental conditions in the pancreatic islets. This strengthens the relevance of these networks in the inflammation that takes place within the islets at the initiation of type 1 diabetes. There also appear to be pathways that are more specific to the process in islets, such as the IL-17 pathway, which were not regulated in any of the five datasets but highly regulated in the human islets after cytokine stimulation.

In conclusion, the current study captured information not obtainable by conventional means of GWAS and facilitated a systems-based understanding of type 1 diabetes. We identified novel type 1 diabetes susceptibility genes that may play a particularly relevant role in complex interactions for type 1 diabetes in human islets. The observations support that the β-cells play an active role in the inflammatory and metabolic processes leading to type 1 diabetes. On the basis of the current findings, novel drugs and treatment strategies may be developed to treat or prevent type 1 diabetes.

Specific parts of the study were funded by the National Institutes of Health (1-DP3-DK-085678-01), the JDRF (33-2008-391), the Sehested Hansen Foundation, and Novo Nordisk A/S. This research uses resources provided by the T1DGC, a collaborative clinical study sponsored by the National Institute of Diabetes and Digestive and Kidney Diseases, the National Institute of Allergy and Infectious Diseases, the National Human Genome Research Institute, the National Institute of Child Health and Human Development, and the JDRF and supported by U01-DK-062418. The work carried out in this study was in part supported by the Novo Nordisk Foundation Center for Protein Research.

No other potential conflicts of interest relevant to this article were reported.

R.B., C.B., J.S., and F.P. researched data, contributed to discussion, and wrote, reviewed, and edited the manuscript. A.P. and L.J.J. researched data, contributed to discussion, and reviewed and edited the manuscript. L.A.B. and T.F. researched data. C.H.B.-B. and K.S.F. researched data and contributed to discussion. F.P. is the guarantor of this work and, as such, had full access to all the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis.

The authors thank Jens Høiriis Nielsen, University of Copenhagen, for providing human islet preparations and Nils Billestrup, University of Copenhagen, for RNA from human islets provided through the JDRF Islet Distribution Program (JDRF award 6-2005-1178). Furthermore, the authors thank Oliver Burren, University of Cambridge, and the T1DGC for access to the GWAS data.