Systematic Evaluation of Genes and Genetic Variants Associated with Type 1 Diabetes Susceptibility

Ramesh Ram, Munish Mehta, Quang T. Nguyen, Irma Larma, Bernhard O. Boehm, Flemming Pociot, Patrick Concannon and Grant Morahan
J Immunol April 1, 2016, 196 (7) 3043-3053; DOI: https://doi.org/10.4049/jimmunol.1502056

Abstract

Genome-wide association studies have found >60 loci that confer genetic susceptibility to type 1 diabetes (T1D). Many of these are defined only by anonymous single nucleotide polymorphisms: the underlying causative genes, as well as the molecular bases by which they mediate susceptibility, are not known. Identification of how these variants affect the complex mechanisms contributing to the loss of tolerance is a challenge. In this study, we performed systematic analyses to characterize these variants. First, all known genes in strong linkage disequilibrium (r2 > 0.8) with the reported single nucleotide polymorphisms for each locus were tested for commonly occurring nonsynonymous variations. We found only a total of 22 candidate genes at 16 T1D loci with common nonsynonymous alleles. Next, we performed functional studies to examine the effect of non-HLA T1D risk alleles on regulating expression levels of genes in four different cell types: EBV-transformed B cell lines (resting and 6 h PMA stimulated) and purified CD4+ and CD8+ T cells. We mapped cis-acting expression quantitative trait loci and found 24 non-HLA loci that affected the expression of 31 transcripts significantly in at least one cell type. Additionally, we observed 25 loci that affected 38 transcripts in trans. In summary, our systems genetics analyses defined the effect of T1D risk alleles on levels of gene expression and provide novel insights into the complex genetics of T1D, suggesting that most of the T1D risk alleles mediate their effect by influencing expression of multiple nearby genes.

This article is featured in In This Issue, p.2911

Introduction

Type 1 diabetes (T1D) affects ∼30 million people worldwide (1). It is a complex autoimmune disease causing the destruction of pancreatic β cells. The largest genetic studies of T1D have been carried out by the Type 1 Diabetes Genetics Consortium (T1DGC) (2–4). These and other reports have now defined genetic variants associated with T1D in >60 different chromosomal regions (see Ref. 5 for review).

There is a need to identify the causative variants that are in linkage disequilibrium (LD) with the single nucleotide polymorphisms (SNPs) found by such association studies, and to define the molecular bases by which they contribute to disease susceptibility. The challenge of post–genome-wide association studies functional studies (6–8) is in finding ways to translate genetic associations into clinically useful information. The strong genetic association of the disease with HLA class II genes of the MHC is well established (9), but the identity of the genes associated with many of the non-HLA loci remains largely unknown, especially with respect to those associated SNPs located in noncoding regions of the genome (2, 5). Therefore, this study focuses on characterizing the non-HLA T1D risk loci.

In principle, most genetic variants could plausibly affect biological processes by changing amino acid residues in encoded proteins or by changing their levels of expression in particular tissues. Various DNA sequence repositories allow identification of commonly occurring nonsynonymous (missense) variations in genes, and amino acid substitution polymorphisms could be characterized for their potential to affect biological processes (10). Expression quantitative trait locus (eQTL) analyses can identify genes whose variation in expression is associated with specific SNP markers. For example, sequence variation in promoters or enhancer elements could result in differential cis regulation. Genetic variants can also regulate expression of genes at greater distances from, or on different chromosomes than, the regulatory element, that is, trans regulation (11). The mechanisms involved in trans regulation could include indirect genetic effects, for example, by means of variation in encoded proteins such as transcription factors, or by other effects, such as steric hindrance (11). Some loci could exert both cis and trans effects.

In the present study, we performed systems genetics (12) analyses of the 55 loci (2, 13–25) (Table I) showing highest evidence of association with T1D using data generated by the T1DGC (2) and Immunochip projects (13). Additionally, four new SNPs (rs6691977, rs4849135, rs2611215, and rs11954020) that showed strong associations (p < 5 × 10−8) with T1D in (13) were included in our study. SNPs at these loci were assessed for disease gene candidacy. Expression data of 47,323 high-quality transcripts (Illumina, HT-12 V4) were correlated with SNPs reported in T1D loci adjusting for confounding factors such as population structure.

Materials and Methods

Study samples

The T1DGC study has been described elsewhere, including phenotypic and extensive genetic characterization of >4000 affected sibling-pair families (3). Upon joining the T1DGC, family members provided blood samples. PBMCs were isolated and aliquots were used to provide DNA samples, to derive EBV-transformed B (EBV-B) lymphoblastoid cell lines (26, 27), and they were frozen for later use. EBV-B cells from 202 European subjects from the T1DGC family collection were examined in the present study. These samples consisted of 46 unaffected subjects and the rest were T1D cases. EBV-B cells were either unstimulated or treated with PMA (28) for 6 h (26, 27). PMA-stimulated samples consisted of 49 unaffected subjects. Cell lines were stimulated on a second occasion to provide a duplicate sample. SNPs were genotyped using the Immunochip (13) platform.

Frozen PBMC samples from 113 T1DGC family members were thawed, cultured overnight, stained, and separated into CD4+ and CD8+ T cell populations by flow sorting. Sufficient RNA was obtained from 102 CD4+ T cell samples and 84 CD8+ T cell samples to perform microarrays. Sex, HLA-DR, and autoantibody statuses of the affected subjects are summarized in Supplemental Table Ii.

Microarray analyses

After cell culture or flow sorting, RNA was extracted using TRIzol reagent (Invitrogen) following the manufacturer’s instructions. The RNA quantity was measured by a NanoDrop 1000 spectrophotometer (Thermo Scientific), and RNA quality was checked on an Agilent 2100 Bioanalyzer (Agilent Technologies). Samples with an RNA integrity number of ≥8 were biotin labeled using an Illumina TotalPrep RNA amplification kit (Ambion) as per the manufacturer’s instructions. The biotin-labeled samples were hybridized onto Illumina HumanHT-12 v4.0 expression beadchips and beadchips were scanned by a BeadArray Reader (Illumina) following the manufacturer’s instructions. Raw data were finally exported by GenomeStudio software (Illumina) for analysis.

Microarray and eQTL analysis

Genome-wide gene expression values from GenomeStudio (Illumina) for each of 47,323 probes were subjected to background correction using control probe profile, variance stabilizing transformation, and robust spline normalization using the lumi package (29) in R. We then removed from the analysis 95 transcripts that are ERCC spike-in controls (having gene symbols starting with ERCC). Four separate gene expression datasets were created. Upon examining initial principal component (PC) analysis (PCA) plots, batch effects were evident. For correction within each cell type, normalized expression data for each gene were centered by batch and centered again after merging batches. The batch correction was validated by PCA (Supplemental Fig. 1A, 1B), and pair plots of PCs 1–4 did not reveal any further batch effects. BLASTN software was used to identify probesets with unique sequences.

The data generated in this study have been deposited in the National Center for Biotechnology Information Gene Expression Omnibus (series accession no. GSE77350) and are available at http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE77350.

To assess association between SNP genotype and gene expression, the Matrix eQTL R (30) package was used. To adjust for unknown confounders in the expression, two correction methods were used and results were compared.

RUV-2 correction

First, association of all T1D SNPs with normalized uncorrected data was performed and p value association of every SNP–gene pair was obtained. For each SNP, the top 5000 associated genes ranked by p value were excluded and the rest were treated as empirical controls for RUV-2 correction (31, 32) using the naiveRandRUV method, with parameter k set to 20. After correction, the same SNP was tested against the corrected set and p value association of the SNP–gene pair was recorded. This procedure was repeated for all SNPs and finally Benjamini false discovery rate (FDR) correction was applied to the set of recorded nominal p values.

PCA correction

PCs were derived from individual whole-expression sets and tested against whole-genome Immunochip SNPs (200,000). The PCs that showed no or weak genome association (i.e., minimum SNP-PC association FDR p > 0.001) were chosen as unassociated PCs (33). These PCs were incrementally added in their order of precedence as covariates to assess SNP-gene associations with an aim to maximize the number of significant cis gene detections (at FDR p < 0.001) for the 77 T1D SNPs tested. Based on analysis shown in Supplemental Fig. 1E and 1F, the four gene expression datasets were corrected as follows: seven PCs (1–6 and 8) were removed from EBV-B basal cell line samples, three PCs (1, 4, and 9) were removed for PMA-stimulated EBV-B cell line samples, four PCs (1–4) were removed for CD4+ samples, and two PCs (1 and 2) were removed for CD8+ samples.

We compared numbers of cis- and trans-regulated genes detected in each cell type using two methods (Supplemental Table Iii). The RUV-2 method of correction yielded more significant results than did PCA methods.

Statistical analysis

Differential gene expression analysis was performed using the Limma package written for R (34). Transmission disequilibrium test (sibship) tests were performed using the software package UNPHASED (35, 36).

Enrichment analysis

Candidate gene names were converted to Entrez gene IDs and were analyzed using the DAVID (37, 38) function annotation tool (http://david.abcc.ncifcrf.gov/). Further pathway and network analyses were performed using GATHER (http://gather.genome.duke.edu) (39) and GENEMANIA (http://www.genemania.org) (40), respectively.

Results

Systematic evaluation of nonsynonymous SNPs in genes in T1D-associated regions

First, we searched for commonly occurring nonsynonymous SNPs (nsSNPs) in LD (r2 > 0.8) with the T1D SNPs (2, 13–25) in the 1000 genomes and HAPMAP (41) CEU datasets. All amino acid substitutions were subject to prediction of the effect of these changes, evaluated as benign, probably damaging, or possibly damaging by PolyPhen-2 (10). This search returned 25 nsSNPs in strong LD with only 16 of the 60 non-HLA T1D loci. These SNPs occurred in a total of 22 unique genes. The seven potentially damaging effects were found in two genes, SULT1A2 and GSDMB. Prediction status does not affect candidacy per se, so all genes listed in Table II should be evaluated in further studies. Additionally, none of the four SNPs recently discovered in Onengut-Gumuscu et al. (13) was in strong LD (r2 > 0.8) with any nsSNPs. Among the LD SNPs, there were three splice-region variants and one stop-gain variant (summarized in Supplemental Table Iiii).

View this table:
Table I. List of reported T1D SNPs located in 59 non-HLA T1D loci
View this table:
Table II. nsSNPs found in LD (r2 > 0.8) with 16 T1D loci

Next, we searched whether any nsSNPs showed better association with T1D than did the reported SNP itself. For this, we performed a transmission disequilibrium test (sibship) using UNPHASED (35, 36) on a dataset of 2676 nuclear families with unaffected parents and two or more affected siblings. Results are presented in Table II. Association p values for SNPs not included in the Immunochip genotyping were derived from Barrett et al. (2). At six T1D loci, the nsSNPs were the reported best SNPs. From those nsSNPs that were genotyped by Immunochip, rs7498665 associated with SH2B1 showed slightly better association than the reported rs4788084 (Δp = 0.1, where Δp = pns SNP/preported SNP). Two other ns SNPs (rs2305480 and rs229527) also showed very small (Δp > 0.1) improvement in association compared with the reported T1D SNP. Most of the T1D loci did not have associated nsSNPs in nearby genes.

Gene expression analyses

EBV-B cell lines were produced from blood samples obtained from T1DGC family members (3). RNA was extracted from 202 available EBV-B cell lines that were cultured under basal conditions and stimulated with PMA. We also purified CD4+ and CD8+ T cells from peripheral blood samples provided by 113 subjects. None of these subjects overlapped with the donors of the 202 EBV samples. After quality control, sufficient high-quality RNA to perform microarrays was obtained from 102 CD4+ T cell samples and 84 CD8+ T cell samples. The EBV-B cell samples were derived from both T1D cases and unaffected subjects. The unaffected controls included were first-degree relatives of the subset of case samples, and islet autoantibody status was not determined for these unaffected subjects. The details regarding the autoantibody, sex, and HLA-DR status of the affected subjects are summarized in Supplemental Table Ii. As expected, there were no significant differences in the gene expression between cases and unaffected subjects or between cases and unaffected first-degree relatives (Supplemental Fig. 1C, 1D), so all samples were used to search for eQTLs.

These RNA samples were hybridized to Illumina microarrays (HT-12v4). Data processing was carried out as described in Materials and Methods. Batch effects were corrected for each cell type by centering the normalized gene expression data by batch and centering again after merging batches. The batch correction was validated by PCA (Supplemental Fig. 1A, 1B). To eliminate probesets with potential cross-hybridization problems, a BLAST search of each probe sequence was carried out on a custom database of all 47,323 Illumina probeset sequences, and 38,500 probes that had a single hit were retained. In doing so, probes associated with two known T1D candidates RPS26 (due to sequence similarity with probes associated with RPS26 pseudogenes) and DEXI (due to sequence similarity with a probe associated with LOC653752) were removed. There were 95 ERCC spike-in controls in the probeset, which were also excluded from analysis. We also performed a search for SNPs within probeset coordinates and excluded any probes that contained SNPs from further analysis. We performed differential expression analysis of unstimulated EBV-B cells and after 6 h PMA stimulation. The negative log10 (adjusted p value) of each probe showing differential expression was plotted against the log2 fold change in a volcano plot (Fig. 1). Adjusted p < 0.0001 was selected as a cut-off for differential expression. A total of 1465 genes were differentially expressed at this threshold with at least a modest fold change (absolute log2 fold change > 0.3). Genes with the highest fold changes in expression included CCL3, CCL4, EGR1, EGR2, DUSP21, PIP4K2C, ILDR1, and IL9R.

FIGURE 1.
FIGURE 1.

Comparison of gene expression in EBV-B cells between basal and 6 h PMA-stimulated samples. Differences (log2 fold change) in gene expression are shown on the x-axis; y-axis shows −log10 (adjusted p values).

Parameters for systems genetics analyses

Genotypes of T1DGC subjects were previously determined (2–4, 13) at 77 SNPs in 55 of 60 T1D risk loci (Table I). Based on the risk allele’s code at each T1D SNP, an additive recode [0,1,2] was applied so that the risk allele’s effect on gene expression could be determined. Separate analyses were performed for each of the four expression sets (EBV-B basal, EBV-B 6 h PMA-stimulated, CD4+ and CD8+). For these analyses, we conservatively defined a cis transcript as being from a gene whose transcription start or end site was located within 1 Mbp from the T1D SNP. A trans-regulated transcript was defined as a gene located elsewhere in the genome. For each set, 3672 cis interactions pairs were tested, ∼2.9 million trans interactions pairs were tested, and FDR p value corrections were applied separately for cis and trans eQTLs. The Matrix eQTL R package (30) was used to perform these eQTL tests. Owing to unknown confounding factors that could limit the power of detecting significantly differentially expressed genes, we performed two methods of correction independently: 1) removing unwanted variation (RUV-2) (31, 32), and 2) adding genome-wide unassociated expression–derived PCs as covariates (described in Materials and Methods).

All transcripts with FDR p < 0.05 for each T1D SNP were followed up with enrichment analysis using the DAVID bioinformatics resource (37, 38). Additional pathway and network analyses were performed using GATHER (39) and GENEMANIA (40), respectively. The results from these analyses are summarized in Tables III–IV and are described below. Box plots of eQTL associations can be accessed online through our Web resource (42) where we compare effects explained by raw normalized gene expression against RUV-2– and PCA-corrected gene expression sets. A screenshot of the user interface is shown in Fig. 2.

View this table:
Table III. cis Genes associated with 24 T1D SNPs with minimum FDR p < 0.001
View this table:
Table IV. Transgenes associated with 25 T1D SNPs with minimum FDR p < 0.001
FIGURE 2.
FIGURE 2.

Screen capture of the Web interface for browsing box plots and gene networks (available at: http://www.sysgen.org/T1DGCSysGen/).

Effect of T1D-associated non-HLA SNPs on neighboring gene expression in EBV-B cell lines

We examined cis genes in EBV-B basal cell line samples at various FDR p value thresholds. At p < 0.001, 15 T1D SNPs were associated with differences in expression of 20 genes (Table III). Using lower thresholds of adjusted p values (<0.05), an additional 13 T1D SNPs affected the expression of a further 20 genes (Supplemental Table IIi). Hence, 28 T1D SNPs were found to be associated with changes in a total of 40 significant cis genes. Of these, three SNPs (rs10877012, rs4788084, and rs2290400) showed strong cis effects with multiple nearby genes that were either up- or downregulated by the corresponding risk allele. In testing the four newly discovered T1D SNPs (13), we observed that the risk allele associated with rs2611215 reduced expression of TMEM192 (FDR p = 0.008) (Supplemental Table IIi).

Next, we tested 6 h PMA-stimulated EBV-B cell line samples. Results confirmed the cis effects associated with 22 of 40 candidate genes identified in unstimulated EBV-B cells (at minimum FDR p < 0.05) and the effect directions were consistent. IFNGR1, SUOX, SPNS1, and UBASH3A were among genes that showed regulatory effects in basal cells but not after PMA stimulation. Additionally, 17 T1D SNP genotypes significantly regulated the expression of 16 new candidate genes (FDR p < 0.05). Of these, genes INO80B and LYRM2 were detected highly significant at FDR p < 0.001 (Table III). The expression of candidate genes IKZF1 and TSFM showed decreased association with their corresponding T1D SNPs after stimulation, compared with basal condition (refer to Ref. 42). The rest of these results are presented in Supplemental Table IIi.

In summary, 31 T1D SNPs affected the expression of a total of 38 candidate cis genes, 22 of which had shown evidence of cis effects in unstimulated EBV-B cells, whereas the remaining 9 showed association after PMA stimulation, thus suggesting genes that may play a role after immune activation.

Effect of T1D-associated non-HLA SNPs on neighboring gene expression in CD4+ and CD8+ T cells

Tests of CD4+ T cell samples revealed that 16 T1D SNP genotypes regulated the expression of 20 genes significantly. Of these genes, 11 (SMARCE1, LOC728734, SUOX, FAM119B, C16ORF75, GSDMB, IKZF1, ADCY3, ORMDL3, SKAP2, and IKZF3) were found to be cis regulated in both EBV-B and CD4+ T cells by the same T1D SNPs (Table III). In particular, the risk allele of rs2290400 (T) affected nearby genes ORMDL3, GSDMB, and IKZF3 similar to that observed in EBV-B cells. The effect directions between the cell types for the 11 shared genes were consistent, except for gene C16ORF75 where the risk allele increased expression in CD4+ cells but decreased it in EBV-B cells (42). We also noted that expression of candidate gene SUOX showed a clear increase in the significance of association (i.e., lower p value) with T1D risk allele rs705704 (T) in the CD4+ cells compared with EBV-B cells. Additionally, there were nine newly identified candidate genes associated with nine T1D SNPs. Five of these SNPs had shown cis effects in the EBV-B cells but had affected a different set of genes. Among these nine new candidate genes, CLECL1 was the most significantly associated (Table III). The cis genes detected at lower FDR thresholds of 0.01 and 0.05 are presented in the Supplemental Table IIi. These results suggest that the effects of the T1D risk SNPs on gene expression vary between cell types.

Finally, we performed analyses of the CD8+ T cell samples and identified 17 T1D SNP genotypes that regulated the expression of 19 genes across all samples tested. Excepting ADCY3, 10 candidate genes were found cis regulated in EBV-B cells, CD4+ T cells, and CD8+ cells. Thirteen of the 19 candidate genes were cis regulated in both CD4+ and CD8+ T cells and the effect directions were consistent. The remaining six that were not differentially regulated in EBV-B cells or in CD4+ cells were associated with six T1D SNPs in CD8+ cells. Of these, T1D SNP rs2292239 regulated the expression of candidate gene ERBB3 most significantly (FDR p < 0.001) (Table III). The rest of the results are presented in Supplemental Table IIi.

In summary, 24 T1D SNP genotypes regulated the expression of 31 candidate genes highly significantly at FDR p < 0.001 (Table III). Using lower FDR-adjusted p value thresholds (p < 0.05), 43 T1D SNP genotypes regulated the expression of 71 candidate genes. Using an even less stringent suggestive threshold of nominal unadjusted p < 0.001 for evidence of cis effect, we could define up to 85 candidate genes that were affected by 50 T1D SNPs in the four cell types tested.

T1D-associated SNPs associated with changes in expression of distant genes

Next, we investigated whether T1D loci showed trans regulatory effects. After performing ∼2.9 million tests for each cell type and appropriate statistical correction, we identified 38 genes that were highly significantly associated with 25 T1D SNPs at FDR p < 0.001 (Table IV). Five of these SNPs (rs1534422, rs1990760, rs11571291, rs9585056, and rs425105) did not show any cis effect on nearby genes in the cell types tested. Trans-regulated genes shared between B and T cells were detected at only one T1D locus (defined by T1D SNP rs705704) and the effect direction was consistent. Except for ZMYM5, GRAMD1B, and LOC389386, all significant trans genes were detected in the EBV-B cells. Upon characterizing the function of 38 trans genes in DAVID (37, 38), we identified two clusters: CD276, ST6GAL1, CCL5, and IRF8 were associated with immune response, and a further two genes (ID2 and IRF8) were associated with immune system and hemopoietic (lymphoid) organ development. Eight T1D SNPs (Table IV, highlighted in bold) showed highly significant cis as well as trans regulatory interactions in one or more cell types tested, suggesting coregulation between cis and trans genes. We describe tests for meaningful relationships between these genes in the next section.

In summary, in addition to the loci that affected genes in cis, we could identify five loci that exclusively affected genes in trans. Of the T1D loci that were not associated with expression changes in any of the four cell types, three loci contained nsSNPs defined in Table II. The trans regulatory effects detected at lower threshold levels are presented in Supplemental Tables IIii and III.

Enrichment analysis of genes associated with T1D susceptibility

We investigated the function of the genes whose expression was changed by individual risk SNPs. The DAVID enrichment analysis software (37, 38) tests whether sets of genes are enriched for terminology referenced by UniProt Protein Information Resource keywords, Gene Ontology (GO), and KEGG pathways. First, we performed analysis to explore for enrichment between the highly significant (FDR p < 0.001) cis and trans gene candidates for the eight T1D SNPs highlighted in Table IV. For three of these SNPs, the candidate genes shared a common keyword (Table V). Second, using the list of 86 candidate genes derived from Tables II–IV, we performed pathway and enrichment analysis using GATHER (39) and we report results obtained with high confidence (unadjusted p < 0.001) in Table VI. In these results, we found that the cytokine–cytokine receptor interaction pathway received the highest significance. Third, we performed network analysis using GENEMANIA (40) for the same list of 86 candidate genes. The significant functional findings are presented in Table VI. The full GENEMANIA report can be accessed online (http://www.sysgen.org/T1DGCSysGen/genemania.pdf). Finally, we analyzed the list of cis and trans genes detected at FDR p < 0.05 for every T1D SNP separately. We identified 21 enrichment terms (excluding GO cellular component terms) that were significantly enriched at Benjamini p < 0.05 for 10 T1D SNPs. These results are summarized in Table VII and below.

View this table:
Table V. Enrichment terms shared between cis genes and trans genes in association with three T1D SNPs
View this table:
Table VI. Network and pathway analysis of the list of candidate genes identified in Tables II–IV using GATHER and GENEMANIA (unadjusted p < 0.001)
View this table:
Table VII. Significant enrichment terms found using the DAVID bioinformatics resource

The term “lectin” was highly enriched for the T1D locus defined by rs10466829 because it affected expression of five c-type lectin genes (CLEC1A, CLEC2B, CLEC2D, CLECL1, and CD69) in the cell types tested. The T1D locus defined by rs17696736 was highly enriched for response to virus and anti-viral defense due to changes in expression of seven trans genes (EIF2AK2, IFI16, IFNGR1, MX1, MX2, PLSCR1, and STAT1). Furthermore, genes MX1 and MX2 are also known inflammatory and immune response genes. Additionally, the T1D SNP rs416603 showed significant enrichment for IL10–anti-inflammatory signaling pathway and intestinal immune network IgA production pathway through its regulation of three genes (IL10, IL10RA, and STAT5A). We also noted that two risk SNPs (rs2476601 and rs679574) showed association in trans with genes in the MHC (HLA-F, HLA-G, HLA-H, and DRB4), which gave positive enrichment for terms such as “antigen processing and presentation.” These results provide insights into the functions of genes whose expression is affected by the T1D loci.

Validation of trans regulatory gene interactions

To confirm our results, we searched using the blood eQTL browser (43) for the trans regulatory associations we identified at significance threshold FDR p < 0.05. Because not all T1D SNPs may be present in this browser, we allowed a 100 Kb window for the search of the expression SNP. Two trans genes were validated: UBE2L6 (EBV-B with/without PMA) associated with rs3184504 and STAT1 (EBV-B basal) with rs17696736. Second, we searched in the trans regulatory interactions reported by Fairfax et al. (44) and validated a further three gene interactions reported in their study: LOC728823, IP6K2, and LOC389386 all associated with the T1D SNP rs705704. Although many cis gene effects were clearly defined from our datasets, validating trans genes poses a challenge warranting further investigation.

Discussion

Our results provide a potential molecular basis for disease association at 46 of the 59 identified T1D loci (Table I). Sixteen of these loci contained nsSNPs in strong LD with the T1D SNP. Thirty-six of the loci showed cis effects on 75 nearby genes. The remainder showed statistically significant trans regulatory interactions that were substantiated by significant enrichment results (Tables V–VII). These candidate genes can be the focus for further studies. For example, a systems genetics study (45) into candidate gene CTSH, whose expression was affected by T1D SNP rs3825932, supported its product as a novel therapeutic target.

Onengut-Gumuscu et al. (13) recently confirmed several previously reported T1D-associated SNPs (2, 5) in addition to the identification of four additional new T1D risk SNPs of which one SNP (rs2611215) had high significance (p = 1.817 × 10−11) whereas p values of the rest only just exceeded the significance threshold (p < 5 × 10−8). This study found that the associated SNPs localized to enhancer sequences active in thymus, T and B cells, and CD34+ stem cells. Of the four new T1D-associated SNPs (13) we were able to establish likely candidacy for rs2611215 as TMEM192.

An important conclusion from our study is that the cell type was important in characterizing T1D SNP function, that is, eQTLs are cell type–specific. For example, the candidate gene ERBB3 was highly significantly cis regulated in CD8+ T cells but its variation effect was largely undetectable in other cell types. The risk allele associated with rs4788084 reduced expression of candidate gene TUFM exclusively in the CD4+ cells. Similarly, CLECL1 did not show any effect in EBV-B cell lines but showed highly significant effects in both T cell types tested. Among the weakly detected effects, there was evidence that suggested the risk allele associated with rs231727 reduced expression in cis of a well-known candidate (CTLA4) exclusively in the CD8+ cells (unadjusted p = 0.0003, FDR p = 0.04) (Supplemental Table IIi). Our CD4+/CD8+ cell type data also assisted in mapping candidate genes at otherwise anonymous T1D SNPs; the most significant of these candidates included SLC11A1 (rs3731865), C6Orf173 (rs9388489), and C10orf59 (rs10509540).

Sixteen transcripts (12 in cis, 4 in trans) were significantly associated with T1D SNPs in both EBV-B and the T cell types tested. Of these, a novel uncharacterized cis transcript LOC728734 (nuclear pore complex interacting protein family, member B8) was identified to be associated with T1D SNP rs4788084 (chromosome [Chr ] 16p11.2) where the risk allele decreased expression in all four cell types. The effect directions of cis and trans regulation by T1D SNPs on genes detected across multiple cell types were found consistent for all SNPs except C16ORF75. We also noted that probes associated with candidate genes DEXI and RPS26 also showed a strong cis regulatory effect in association with T1D risk SNPs rs12708716 and rs705704, respectively, in one or more cell types. However, due to quality control procedures relevant to cross-hybridization problems described in the previous section, these probes were excluded from further analysis. nsSNPs may also affect gene expression in trans. We found two examples of these: rs1990760 (Chr 2q24.2) in IFIH1 also affected the expression of LOC643997 in trans; similarly, rs2304256 (Chr 19p13.2) in TYK2 also affected the expression of ZNF280D in trans.

Pathway analysis identified the cytokine–cytokine receptor interaction pathway with highest confidence. The sulfur metabolism pathway also scored high significance because two genes SUOX (cis) and SULT1A2 (nonsynonymous) involved in this pathway were identified as candidates in this study. It is also well known that sulfur plays an important role in insulin production (for review, see Ref. 46). Furthermore, DAVID enrichment analysis of locus-specific cis and trans transcript perturbations revealed significant enrichment of 48 category terms in 15 of the T1D regions at FDR p < 0.05. Among the best enriched terms were “response to virus,” “acetylation,” “lectin,” and “IL10-anti-inflamatory pathway.” From the enrichment analysis for genes associated with each T1D SNP, upon examination T1D risk SNP rs17696736 (Chr 12q24.12) was notably associated with response to virus and antiviral defense due to trans genes that are involved in proinflamatory response (such as MX1 and MX2) in the Salmonella infection pathway (KEGG pathway 05132). In contrast, chemokine gene CCL5 was highly significantly associated with diabetes loci associated with T1D SNP rs425105 (Chr 19q13.32). These results support evidence found in a recent work (47) suggesting that Salmonella and chemokine vaccines can prove clinically useful in diabetes management and prevention.

In conclusion, our results confirm systems genetics (12) as a powerful tool for investigating the genetic architecture of complex diseases such as T1D. Many genes were identified whose expression levels were influenced by SNPs associated with T1D susceptibility. These nsSNPs, cis-regulated genes, and trans-regulated genes we identified are important candidates for further investigation. So that other researchers can extend the work reported in the present study, we have implemented a Web interface (42) allowing users to browse box plots for the eQTL interactions reported.

Disclosures

The authors have no financial conflicts of interest.

Acknowledgments

This research utilized resources provided by the Type 1 Diabetes Genetics Consortium, 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 by the Juvenile Diabetes Research Foundation International.

Footnotes

  • This work was supported by National Health and Medical Research Council of Australia Program Grants 53000400 and 37612600, the Diabetes Research Foundation (Western Australia), and National Institute of Diabetes and Digestive and Kidney Diseases Grant 1DP3DK085678. R.R. is supported by the MACA Ride to Conquer Cancer in association with the Harry Perkins Institute of Medical Research. B.O.B. is supported by the Deutsche Forschungsgemeinschaft and by a grant from the Boehringer Ingelheim Ulm University BioCenter.

  • The data presented in this article have been submitted to the National Center for Biotechnology Information Gene Expression Omnibus (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE77350) under accession number GSE77350.

  • The online version of this article contains supplemental material.

  • Abbreviations used in this article:

    Chr
    chromosome
    EBV-B
    EBV-transformed B
    eQTL
    expression qualitative trait locus
    FDR
    false discovery rate
    GO
    Gene Ontology
    LD
    linkage disequilibrium
    nsSNP
    nonsynonymous SNP
    PC
    principal component
    PCA
    principal component analysis
    SNP
    single nucleotide polymorphism
    T1D
    type 1 diabetes
    T1DGC
    Type 1 Diabetes Genetics Consortium.

  • Received September 24, 2015.
  • Accepted January 25, 2016.

References

  1. ↵
    1. Danaei G.,
    2. M. M. Finucane,
    3. Y. Lu,
    4. G. M. Singh,
    5. M. J. Cowan,
    6. C. J. Paciorek,
    7. J. K. Lin,
    8. F. Farzadfar,
    9. Y. H. Khang,
    10. G. A. Stevens,
    11. et al,
    12. Global Burden of Metabolic Risk Factors of Chronic Diseases Collaborating Group (Blood Glucose)
    . 2011. National, regional, and global trends in fasting plasma glucose and diabetes prevalence since 1980: systematic analysis of health examination surveys and epidemiological studies with 370 country-years and 2·7 million participants. Lancet 378: 31–40.
    OpenUrlCrossRefPubMed
  2. ↵
    1. Barrett J. C.,
    2. D. G. Clayton,
    3. P. Concannon,
    4. B. Akolkar,
    5. J. D. Cooper,
    6. H. A. Erlich,
    7. C. Julier,
    8. G. Morahan,
    9. J. Nerup,
    10. C. Nierras,
    11. et al
    . 2009. Genome-wide association study and meta-analysis find that over 40 loci affect risk of type 1 diabetes. Nat. Genet. 41: 703–707.
    OpenUrlCrossRefPubMed
  3. ↵
    1. Rich S. S.,
    2. B. Akolkar,
    3. P. Concannon,
    4. H. Erlich,
    5. J. E. Hilner,
    6. C. Julier,
    7. G. Morahan,
    8. J. Nerup,
    9. C. Nierras,
    10. F. Pociot,
    11. J. A. Todd
    . 2009. Overview of the Type I Diabetes Genetics Consortium. Genes Immun. 10(Suppl. 1): S1–S4.
    OpenUrlCrossRefPubMed
  4. ↵
    1. Morahan G.,
    2. M. Mehta,
    3. I. James,
    4. W.-M. Chen,
    5. B. Akolkar,
    6. H. A. Erlich,
    7. J. E. Hilner,
    8. C. Julier,
    9. J. Nerup,
    10. C. Nierras,
    11. et al
    . 2011. Tests for genetic interactions in type 1 diabetes: linkage and stratification analyses of 4,422 affected sib-pairs. Diabetes 60: 1030–1040.
    OpenUrlAbstract/FREE Full Text
  5. ↵
    1. Morahan G.
    2012. Insights into type 1 diabetes provided by genetic analyses. Curr. Opin. Endocrinol. Diabetes Obes. 19: 263–270.
    OpenUrlCrossRefPubMed
  6. ↵
    2010. On beyond GWAS. Nat. Genet. 42: 551.
    OpenUrlCrossRefPubMed
    1. Freedman M. L.,
    2. A. N. Monteiro,
    3. S. A. Gayther,
    4. G. A. Coetzee,
    5. A. Risch,
    6. C. Plass,
    7. G. Casey,
    8. M. De Biasi,
    9. C. Carlson,
    10. D. Duggan,
    11. et al
    . 2011. Principles for the post-GWAS functional characterization of cancer risk loci. Nat. Genet. 43: 513–518.
    OpenUrlCrossRefPubMed
  8. ↵
    1. Sacks D. B.,
    2. M. Arnold,
    3. G. L. Bakris,
    4. D. E. Bruns,
    5. A. R. Horvath,
    6. M. S. Kirkman,
    7. A. Lernmark,
    8. B. E. Metzger,
    9. D. M. Nathan,
    10. National Academy of Clinical Biochemistry,
    11. Evidence-Based Laboratory Medicine Committee of the American Association for Clinical Chemistry
    . 2011. Guidelines and recommendations for laboratory analysis in the diagnosis and management of diabetes mellitus. Diabetes Care 34: e61–e99.
    OpenUrlAbstract/FREE Full Text
  9. ↵
    1. Morahan G.,
    2. M. Varney
    . 2010. The genetics of type 1 diabetes. In The HLA Complex in Biology and Medicine: A Resource Book, 1st Ed. Mehra N. K., G. Kaur, J. McCluskey, F. T. Christiansen, F. H. Claas, eds. Jaypee Brothers Medical Publishers, New Delhi, India, p. 205–218.
  10. ↵
    1. Adzhubei I. A.,
    2. S. Schmidt,
    3. L. Peshkin,
    4. V. E. Ramensky,
    5. A. Gerasimova,
    6. P. Bork,
    7. A. S. Kondrashov,
    8. S. R. Sunyaev
    . 2010. A method and server for predicting damaging missense mutations. Nat. Methods 7: 248–249.
    OpenUrlCrossRefPubMed
  11. ↵
    1. Veyrieras J. B.,
    2. S. Kudaravalli,
    3. S. Y. Kim,
    4. E. T. Dermitzakis,
    5. Y. Gilad,
    6. M. Stephens,
    7. J. K. Pritchard
    . 2008. High-resolution mapping of expression-QTLs yields insight into human gene regulation. PLoS Genet. 4: e1000214.
    OpenUrlCrossRefPubMed
  12. ↵
    1. Morahan G.,
    2. R.W. Williams
    . 2007. Systems genetics: the next generation in genetics research? Novartis Found. Symp. 281: 181–188.
    OpenUrlPubMed
  13. ↵
    1. Onengut-Gumuscu S.,
    2. W.-M. Chen,
    3. O. Burren,
    4. N. J. Cooper,
    5. A. R. Quinlan,
    6. J. C. Mychaleckyj,
    7. E. Farber,
    8. J. K. Bonnie,
    9. M. Szpak,
    10. E. Schofield,
    11. et al
    . 2015. Fine mapping of type 1 diabetes susceptibility loci and evidence for colocalization of causal variants with lymphoid gene enhancers. Nat. Genet. 47: 381–386.
    OpenUrlCrossRefPubMed
    1. Cooper J. D.,
    2. D. J. Smyth,
    3. A. M. Smiles,
    4. V. Plagnol,
    5. N. M. Walker,
    6. J. E. Allen,
    7. K. Downes,
    8. J. C. Barrett,
    9. B. C. Healy,
    10. J. C. Mychaleckyj,
    11. et al
    . 2008. Meta-analysis of genome-wide association study data identifies additional type 1 diabetes risk loci. Nat. Genet. 40: 1399–1401.
    OpenUrlCrossRefPubMed
    1. Todd J. A.,
    2. N. M. Walker,
    3. J. D. Cooper,
    4. D. J. Smyth,
    5. K. Downes,
    6. V. Plagnol,
    7. R. Bailey,
    8. S. Nejentsev,
    9. S. F. Field,
    10. F. Payne,
    11. et al
    . 2007. Robust associations of four new chromosome regions from genome-wide analyses of type 1 diabetes. Nat. Genet. 39: 857–864.
    OpenUrlCrossRefPubMed
    1. Smyth D. J.,
    2. V. Plagnol,
    3. N. M. Walker,
    4. J. D. Cooper,
    5. K. Downes,
    6. J. H. Yang,
    7. J. M. Howson,
    8. H. Stevens,
    9. R. McManus,
    10. C. Wijmenga,
    11. et al
    . 2008. Shared and distinct genetic variants in type 1 diabetes and celiac disease. N. Engl. J. Med. 359: 2767–2777.
    OpenUrlCrossRefPubMed
    1. Bradfield J. P.,
    2. H. Q. Qu,
    3. K. Wang,
    4. H. Zhang,
    5. P. M. Sleiman,
    6. C. E. Kim,
    7. F. D. Mentch,
    8. H. Qiu,
    9. J. T. Glessner,
    10. K. A. Thomas,
    11. et al
    . 2011. A genome-wide meta-analysis of six type 1 diabetes cohorts identifies multiple associated loci. PLoS Genet. 7: e1002293.
    OpenUrlCrossRefPubMed
    1. Yang J. H.,
    2. K. Downes,
    3. J. M. Howson,
    4. S. Nutland,
    5. H. E. Stevens,
    6. N. M. Walker,
    7. J. A. Todd
    . 2011. Evidence of association with type 1 diabetes in the SLC11A1 gene region. BMC Med. Genet. 12: 59.
    OpenUrlPubMed
    1. Wellcome Trust Case Control Consortium
    . 2007. Genome-wide association study of 14,000 cases of seven common diseases and 3,000 shared controls. Nature 447: 661–678.
    OpenUrlCrossRefPubMed
    1. Swafford A. D.,
    2. J. M. Howson,
    3. L. J. Davison,
    4. C. Wallace,
    5. D. J. Smyth,
    6. H. Schuilenburg,
    7. M. Maisuria-Armer,
    8. T. Mistry,
    9. M. J. Lenardo,
    10. J. A. Todd
    . 2011. An allele of IKZF1 (Ikaros) conferring susceptibility to childhood acute lymphoblastic leukemia protects against type 1 diabetes. Diabetes 60: 1041–1044.
    OpenUrlAbstract/FREE Full Text
    1. Evangelou M.,
    2. D. J. Smyth,
    3. M. D. Fortune,
    4. O. S. Burren,
    5. N. M. Walker,
    6. H. Guo,
    7. S. Onengut-Gumuscu,
    8. W. M. Chen,
    9. P. Concannon,
    10. S. S. Rich,
    11. et al
    . 2014. A method for gene-based pathway analysis using genomewide association study summary statistics reveals nine new type 1 diabetes associations. Genet. Epidemiol. 38: 661–670.
    OpenUrlCrossRefPubMed
    1. Bailey R.,
    2. J. D. Cooper,
    3. L. Zeitels,
    4. D. J. Smyth,
    5. J. H. Yang,
    6. N. M. Walker,
    7. E. Hyppönen,
    8. D. B. Dunger,
    9. E. Ramos-Lopez,
    10. K. Badenhoop,
    11. et al
    . 2007. Association of the vitamin D metabolism gene CYP27B1 with type 1 diabetes. Diabetes 56: 2616–2621.
    OpenUrlAbstract/FREE Full Text
    1. Huang J.,
    2. D. Ellinghaus,
    3. A. Franke,
    4. B. Howie,
    5. Y. Li
    . 2012. 1000 Genomes-based imputation identifies novel and refined associations for the Wellcome Trust Case Control Consortium phase 1 data. Eur. J. Hum. Genet. 20: 801–805.
    OpenUrlCrossRefPubMed
    1. Heinig M.,
    2. E. Petretto,
    3. C. Wallace,
    4. L. Bottolo,
    5. M. Rotival,
    6. H. Lu,
    7. Y. Li,
    8. R. Sarwar,
    9. S. R. Langley,
    10. A. Bauerfeind,
    11. et al
    . 2010. A trans-acting locus regulates an anti-viral expression network and type 1 diabetes risk. Nature 467: 460–464.
    OpenUrlCrossRefPubMed
  25. ↵
    1. Wallace C.,
    2. D. J. Smyth,
    3. M. Maisuria-Armer,
    4. N. M. Walker,
    5. J. A. Todd,
    6. D. G. Clayton
    . 2010. The imprinted DLK1-MEG3 gene region on chromosome 14q32.2 alters susceptibility to type 1 diabetes. Nat. Genet. 42: 68–71.
    OpenUrlCrossRefPubMed
  26. ↵
    1. Burster T.,
    2. B. O. Boehm
    . 2010. Processing and presentation of (pro)-insulin in the MHC class II pathway: the generation of antigen-based immunomodulators in the context of type 1 diabetes mellitus. Diabetes Metab. Res. Rev. 26: 227–238.
    OpenUrlCrossRefPubMed
  27. ↵
    1. Rosinger S.,
    2. S. Nutland,
    3. E. Mickelson,
    4. M. D. Varney,
    5. B. O. Boehm,
    6. G. J. Olsem,
    7. J. A. Hansen,
    8. I. Nicholson,
    9. J. E. Hilner,
    10. L. H. Perdue,
    11. et al,
    12. T1DGC
    . 2010. Collection and processing of whole blood for transformation of peripheral blood mononuclear cells and extraction of DNA: the Type 1 Diabetes Genetics Consortium. Clin. Trials 7(1, Suppl.): S65–S74.
    OpenUrlAbstract/FREE Full Text
  28. ↵
    1. Nomura N.,
    2. M. Nomura,
    3. K. Sugiyama,
    4. J. Hamada
    . 2007. Phorbol 12-myristate 13-acetate (PMA)-induced migration of glioblastoma cells is mediated via p38MAPK/Hsp27 pathway. Biochem. Pharmacol. 74: 690–701.
    OpenUrlCrossRefPubMed
  29. ↵
    1. Du P.,
    2. W. A. Kibbe,
    3. S. M. Lin
    . 2008. lumi: a pipeline for processing Illumina microarray. Bioinformatics 24: 1547–1548.
    OpenUrlAbstract/FREE Full Text
  30. ↵
    1. Shabalin A. A.
    2012. Matrix eQTL: ultra fast eQTL analysis via large matrix operations. Bioinformatics 28: 1353–1358.
    OpenUrlAbstract/FREE Full Text
  31. ↵
    1. Risso D.,
    2. J. Ngai,
    3. T. P. Speed,
    4. S. Dudoit
    . 2014. Normalization of RNA-seq data using factor analysis of control genes or samples. Nat. Biotechnol. 32: 896–902.
    OpenUrlCrossRefPubMed
  32. ↵
    1. Jacob L.,
    2. J. A. Gagnon-Bartsch,
    3. T. P. Speed
    . 2016. Correcting gene expression data when neither the unwanted variation nor the factor of interest are observed. Biostatistics 17: 16–28.
    OpenUrlAbstract/FREE Full Text
  33. ↵
    1. Flutre T.,
    2. X. Wen,
    3. J. Pritchard,
    4. M. Stephens
    . 2013. A statistical framework for joint eQTL analysis in multiple tissues. PLoS Genet. 9: e1003486.
    OpenUrlCrossRefPubMed
  34. ↵
    1. Smyth G. K.
    2004. Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat. Appl. Genet. Mol. Biol. 3: 1–25.
    OpenUrlCrossRefPubMed
  35. ↵
    1. Dudbridge F.
    2008. Likelihood-based association analysis for nuclear families and unrelated subjects with missing genotype data. Hum. Hered. 66: 87–98.
    OpenUrlCrossRefPubMed
  36. ↵
    1. Dudbridge F.,
    2. P. A. Holmans,
    3. S. G. Wilson
    . 2011. A flexible model for association analysis in sibships with missing genotype data. Ann. Hum. Genet. 75: 428–438.
    OpenUrlCrossRefPubMed
  37. ↵
    1. Huang W.,
    2. B. T. Sherman,
    3. R. A. Lempicki
    . 2009. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat. Protoc. 4: 44–57.
    OpenUrlCrossRefPubMed
  38. ↵
    1. Huang W.,
    2. B. T. Sherman,
    3. R. A. Lempicki
    . 2009. Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 37: 1–13.
    OpenUrlAbstract/FREE Full Text
  39. ↵
    1. Chang J. T.,
    2. J. R. Nevins
    . 2006. GATHER: a systems approach to interpreting genomic signatures. Bioinformatics 22: 2926–2933.
    OpenUrlAbstract/FREE Full Text
  40. ↵
    1. Warde-Farley D.,
    2. S. L. Donaldson,
    3. O. Comes,
    4. K. Zuberi,
    5. R. Badrawi,
    6. P. Chao,
    7. M. Franz,
    8. C. Grouios,
    9. F. Kazi,
    10. C. T. Lopes,
    11. et al
    . 2010. The GeneMANIA prediction server: biological network integration for gene prioritization and predicting gene function. Nucleic Acids Res. 38: W214-20.
    OpenUrlPubMed
  41. ↵
    International HapMap Project. Available at: http://www.hapmap.org/hapmart.html.en. Accessed: November 27, 2015.
  42. ↵
    T1DGC Systems Genetics Project. Available at: http://www.sysgen.org/T1DGCSysGen/. Accessed: September 23, 2015.
  43. ↵
    1. Westra H. J.,
    2. M. J. Peters,
    3. T. Esko,
    4. H. Yaghootkar,
    5. C. Schurmann,
    6. J. Kettunen,
    7. M. W. Christiansen,
    8. B. P. Fairfax,
    9. K. Schramm,
    10. J. E. Powell,
    11. et al
    . 2013. Systematic identification of trans eQTLs as putative drivers of known disease associations. Nat. Genet. 45: 1238–1243.
    OpenUrlCrossRefPubMed
  44. ↵
    1. Fairfax B. P.,
    2. S. Makino,
    3. J. Radhakrishnan,
    4. K. Plant,
    5. S. Leslie,
    6. A. Dilthey,
    7. P. Ellis,
    8. C. Langford,
    9. F. O. Vannberg,
    10. J. C. Knight
    . 2012. Genetics of gene expression in primary immune cells identifies cell type-specific master regulators and roles of HLA alleles. Nat. Genet. 44: 502–510.
    OpenUrlCrossRefPubMed
  45. ↵
    1. Fløyel T.,
    2. C. Brorsson,
    3. L. B. Nielsen,
    4. M. Miani,
    5. C. H. Bang-Berthelsen,
    6. M. Friedrichsen,
    7. A. J. Overgaard,
    8. L. A. Berchtold,
    9. A. Wiberg,
    10. P. Poulsen,
    11. et al
    . 2014. CTSH regulates β-cell function and disease progression in newly diagnosed type 1 diabetes patients. Proc. Natl. Acad. Sci. USA 111: 10305–10310.
    OpenUrlAbstract/FREE Full Text
  46. ↵
    1. Wallace J. L.,
    2. R. Wang
    . 2015. Hydrogen sulfide-based therapeutics: exploiting a unique but ubiquitous gasotransmitter. Nat. Rev. Drug Discov. 14: 329–345.
    OpenUrlCrossRefPubMed
  47. ↵
    1. Husseiny M. I.,
    2. J. Rawson,
    3. A. Kaye,
    4. I. Nair,
    5. I. Todorov,
    6. M. Hensel,
    7. F. Kandeel,
    8. K. Ferreri
    . 2014. An oral vaccine for type 1 diabetes based on live attenuated Salmonella. Vaccine 32: 2300–2307.
    OpenUrlCrossRefPubMed
Back to top

In this issue

Print
Download PDF
Article Alerts
Email Article
Citation Tools

Jump to section