Genetic Variation Within the HLA-DRA1 Gene Modulates Susceptibility to Type 1 Diabetes in HLA-DR3 Homozygotes

Type 1 diabetes (T1D) is an autoimmune disease that results from destruction of pancreatic islet β-cells by autoreactive T cells (1). These T cells form an inflammatory infiltrate (insulitis) in the islets of affected individuals. Clinical onset is preceded by a prodrome of β-cell autoimmunity marked by the appearance of autoantibodies. It is a non-Mendelian polygenic disorder involving interaction of multiple gene variants, environmental factors, and immunoregulatory dysfunction (1). The major genetic loci for both β-cell autoimmunity and T1D susceptibility are those encoding HLA-DR and -DQ, with a less significant contribution from genes encoding HLA class I (2,3) and HLA-DP (4). Many non-MHC loci are also associated with the disease, but their individual contributions to risk are relatively small (5). Although HLA genotyping can identify most persons at risk for T1D, only 1 in 15 (∼7%) of individuals with even the highest risk HLA genotype (HLA-DR3/4) become diabetic (6). Recent genetic studies have improved our ability to predict T1D in the general population (7), but additional genetic risk likely remains to be discovered (8). Understanding of the detailed genetic risk conferred by HLA-DR and -DQ remains incomplete.

Several “extended DR3” haplotypes containing DRB1*03:01-DQA1*05:01-DQB1*02:01, additional single nucleotide polymorphisms (SNPs), and specific HLA class I alleles are known to modify T1D risk (3,9,10). These extended haplotypes, although statistically significant, have not increased the accuracy of T1D prediction enough to be clinically useful and have not provided complete understanding of the mechanism of genetic risk. Additional data are clearly needed, but genetic heterogeneity and linkage disequilibrium in the highly polymorphic HLA region make identification of other T1D susceptibility loci difficult. Here, we report novel genetic risk elements in a cohort of individuals drawn from large genetic studies; all of whom carried the high-risk homozygous HLA-DR3 genotype to reduce HLA heterogeneity.

Case and control samples of DNA from individuals previously genotyped as homozygous for HLA-DR3 (DRB1*03:01-DQA1*05:01-DQB1*02:01) were selected from collections of three very large published studies. These included both T1D case subjects and control subjects without diabetes from the Type 1 Diabetes Genetics Consortium (T1DGC) (3,11), newly diagnosed T1D patients in the Swedish Better Diabetes Diagnosis (BDD) study (12), newborn children in Sweden screened in the Diabetes Prediction in Skåne (DiPiS) study (13), and children screened but not enrolled in The Environmental Determinants of Diabetes in the Young (TEDDY) study (14). No TEDDY data are included in this report. Numbers obtained from each study are shown in Table 1. Diabetes status was determined at the time of sample acquisition.

Each study cohort included approximately equal numbers of samples from males and females. HLA genotyping of the T1DGC samples was performed by one of the authors (J.A.N.) (15). HLA typing of the BDD and TEDDY samples was performed as previously described (12,14,16). All samples were originally acquired after human subject review. Reanalysis of samples in this report did not require additional approval, as they were completely de-identified.

T1DGC samples were supplied as DNA from immortalized B cells as previously described (15). Swedish samples were supplied as DNA extracted from dried blood spots using QIAamp 96 DNA Blood Kits (Qiagen, Venlo, the Netherlands).

Molecular inversion probe (MIP) capture was used to sequence regions of interest in all samples (17,18). The technique permits identification and quantification of SNPs in all genes and loci targeted for analysis (19). Variant frequency was enhanced by incorporation of unique molecular identifiers within each molecule captured by an MIP (20). MIP procedures were performed as previously described (17). Probes are single-stranded oligonucleotides ∼100 nucleotides in length; all have a shared backbone and also have sequences complementary to the target region. They were designed using algorithms developed by us (17) and synthesized commercially (Integrated DNA Technologies, Coralville, IA). We report here results obtained from two probes designed to target rs9268645, a known genome-wide association study T1D susceptibility locus in intron 1 of HLA-DRA (21). The probe sequences are as follows: 5′-AAGGCAGGGAATGGCTATCACTGNNNNAGATCGGAAGAGCACACGTGACTCGCCAAGCTGAAGNNNNNNNNNNAATTGCCACTATGGGTATGGGA-3′ for MIP-1 and 5′-GTGCTGTTGTAAAGATGTACTGTAGAAAAGNNNNAGATCGGAAGAGCACACGTGACTCGCCAAGCTGAAGNNNNNNNNNNAAGCAACAATTTACCCCAGG-3′ for MIP-2. Samples were sequenced using an Illumina NextSeq 500 (www. illumina.com). We called alleles only on those samples that had a total sequence read depth of ≥10 per locus (17).

Sequence data were processed as previously described (17). One of the two MIP probes (MIP-1) used for genotyping the HLA-DRA intron covered the three SNPs reported here, whereas the second probe (MIP-2) covered only two of the three SNPs. Samples were excluded from analysis if the MIP-2 data were discordant with the MIP-1 data.

Expression of HLA class II genes in defined tissues and cell types was extracted from published databases. Effects of each high-risk SNP allele on gene expression were identified using the database at https://www.gtexportal.org/. We also interrogated the 1000 Genomes Project (22) database for individuals expressing the three-SNP risk haplotype that we discovered (see below) using four populations (CEU, FIN, GBR, and TSI) that reflect the ancestral origins of our T1D cohort (https://hgdownload.cse.ucsc.edu/gbdb/hg19/1000Genomes/phase3/ALL.chr6.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.vcf.gz). Finally, we collected class II gene expression data from experiments using Epstein-Barr virus (EBV)-transformed lymphocytes available for 483 of the individuals in the 1000 Genomes Project (https://www.ebi.ac.uk/arrayexpress/files/E-GEUV-1/GD660.GeneQuantRPKM.txt.gz). Class II gene expression was then stratified according to three-SNP haplotypes of interest.

Odds ratios (ORs) were calculated using the statistical R package “epitools” (https://cran.r-project.org/web/packages/epitools/epitools.pdf). Two-sided P values were calculated using Fisher exact test in epitools to analyze 2 × 2 tables. The χ² statistic was used to analyze higher-order tables with Bonferroni adjustment for multiple comparisons. For evaluating published class II RNA expression data from EBV-transformed lymphoblasts, RPKM (reads per kilobase of transcript per million mapped reads) values were stratified by haplotype and tested for significance by two-way ANOVA (GraphPad Prism), corrected for multiple comparisons by the Holm-Šídák method, with α = 0.05.

We discovered a previously unreported haplotype of three SNPs, rs3135394, rs9268645, and rs3129877, located in intron-1 of the HLA-DRA1 gene (Supplementary Fig. 1). This “tri-SNP” haplotype was identified because the MIPs designed to interrogate rs9268645 revealed the presence of the other SNPs. As shown in Table 2, this tri-SNP haplotype is strongly associated with T1D risk in DR3/3 homozygous individuals. The reference and risk alleles for each SNP are shown in Table 3. We designate the higher-T1D-risk T1D risk haplotype as “010” to designate reference (0), alternate (1), and reference (0) alleles, respectively. Table 2 shows the ORs and P values for the homozygous HLA-DR3 sample in the current study (for people with 010/010 [N = 115], 101/010 [N = 245], or 101/101 [N = 650] haplotypes) relative to the lower-risk 101/101 haplotype. The 010 haplotype is strongly associated with T1D in both the T1DGC and Swedish DR3/3 cohorts; the overall OR for the combined data set was 4.65, P = 1.69 × 10⁻¹³, well below the genome-wide significance threshold of 5 × 10⁻⁸.

In our combined DR3/3 cohort, only 46 of 2,066 (2.2%) chromosomes had haplotypes other than 010 or 101. The frequencies of these SNPs and haplotypes in our HLA-selected T1D samples are quite different from those in the general population, including the CEU, FIN, TSI, and GBR cohorts documented in the 1000 Genomes Project and the cohort in the Swedish population reported in the SweGen Project (23). In those cohorts, non-101 genotypes predominate (Table 4). In contrast, in the targeted cohort of HLA-DR3/3 homozygotes, 91% (Swedish) and 81% (T1DGC) of samples had either the 101/010 or 101/101 genotype (Table 4).

Because the Swedish case and control subjects are largely derived from different studies, we tested for possible geographic stratification but found no significant difference in tri-SNP haplotype as a function of geography (Supplementary Table 1). In addition, risk association remained highly significant (P = 1.17 × 10⁻⁹) when data were reanalyzed restricting the Swedish cases to the Skåne region, from which all the control samples were obtained.

The 010 haplotype appears to confer significant T1D risk compared with the 101 haplotype. OR for the 010/010 genotype is 7.09, relative to the 101/101 genotype, whereas the heterozygous genotype is intermediate (OR 3.42 relative to 101/101). The 101-containing haplotypes (especially 101/101) are significantly protective for T1D compared with the 010 haplotypes. These results effectively stratify the high-risk DR3/3 individuals into two groups: those who carry the 010 haplotype and those who do not (Table 2).

The expression quantitative trait loci (eQTL) database at https://www.gtexportal.org/home/ groups gene expression in numerous tissues and cell types by the allele at any given SNP. Using this eQTL database, we found that RNA expression changes have been associated with the risk alleles of each of our three SNPs and that they statistically significantly affect expression of HLA genes such as DQB1 (data not shown).

The tri-SNP haplotype cannot be studied using the GTEx database because it contains data only for individual SNPs. To study the haplotypes in silico, we adapted a strategy developed previously to investigate susceptibility to multiple sclerosis (24) using both sequence and expression data available through the 1000 Genomes Project. We collected data on HLA class II gene RNA expression in EBV-transformed B lymphocytes and stratified these results by tri-SNP haplotype (N = 93 individuals carried at least one copy of 101 or 010). As noted above, the 101/101 haplotype is rare in the general population, and only two such individuals are included in the 1000 Genomes data set. Nonetheless, despite this limitation, we found that for people carrying the putative protective homozygous haplotype (101/101), class II gene expression in general was lower, and significantly so for HLA-DQB1, compared with that observed in individuals homozygous for the risk haplotype (010/010 [N = 60]) (Supplementary Fig. 2). In the 101/010 heterozygotes (N = 31), the expression of class II genes was intermediate between those values seen in either homozygote.

We report a novel approach to analyzing genetic susceptibility to T1D by minimizing HLA heterogeneity. By studying a population of case and control subjects homozygous for HLA-DR3, we were able to identify a previously unrecognized susceptibility locus within intron 1 of HLA-DRA. The discovery is novel and highly statistically significant. It surpasses the standard rigorous metric of genome-wide significance (5 × 10⁻⁸) despite the relatively small sample size and the lack of polymorphism of the DR-α chain encoded by this locus. This result is consistent with discoveries reported in the multiple sclerosis literature (24), suggesting that occult, HLA-specific susceptibility loci may be important. Our data suggest a possible explanation for the fact that most persons with HLA-DR3 (>90% even for the highest-risk DR3/4 haplotype (6)) do not become diabetic. Notably, 75% of the Swedish DR3/3 samples in this study are homozygous for the relatively protective 101 haplotype.

The current study cohort is limited to people with homozygous “DR3” haplotypes. The predictive tri-SNP haplotype occurs within this cohort, but whether it exists, or is diabetes predictive for other high-risk haplotypes, such as “DR4,” is not known. HLA-DRA1 is nearly invariant, but differences in noncoding regions of HLA-DRA1 do occur. Planned future studies include determining whether the tri-SNP risk haplotype is present in other T1D risk haplotypes.

Studies of RNA expression by others indicate that each SNP in the tri-SNP haplotype has the characteristics of an eQTL that could influence the level of transcription of HLA and other immunologically relevant genes. Our in silico finding is preliminary, but if it can be substantiated directly by in vitro studies, it could contribute to our understanding of how the HLA-DR associations in T1D work. The data suggest that modification of class II gene expression could be a T1D risk mechanism.

In summary, we report that a three-SNP haplotype in intron 1 of HLA-DRA1 appears to modify the risk of T1D in individuals homozygous for the HLA-DR3 haplotype. This intronic SNP haplotype may well function as an eQTL, affecting expression of HLA class II and perhaps other genes. The data also suggest that additional genetic information relevant to susceptibility to T1D remains to be discovered.

BDD Study Group Investigators and Co-Authors. The investigators in the BDD Study Group are as follows: Martina Persson, Department of Medicine, Clinical Epidemiology, Karolinska University Hospital, Stockholm, Sweden; Helena Elding Larsson, Department of Clinical Sciences, Lund University/CRC, Skåne University Hospital, Malmö, Sweden; Gun Forsander, Department of Pediatrics, Institute for Clinical Sciences, Sahlgrenska Academy, University of Gothenburg and the Queen Silvia Children’s Hospital, Sahlgrenska University Hospital, Gothenburg, Sweden; Sten-Anders Ivarsson, Department of Clinical Sciences, Lund University/CRC, Skåne University Hospital, Malmö, Sweden; Johnny Ludvigsson, Department of Clinical and Experimental Medicine, Linköping University, Linköping, Sweden; Claude Marcus, Department of Clinical Science, Karolinska Institutet, Stockholm, Sweden; and Annelie Carlsson, BDD coordinator, Department of Clinical Sciences, Lund University, Skåne University Hospital, Lund, Sweden.

Acknowledgments. The authors thank Gregory Martin (Children’s Hospital Oakland Research Institute, Oakland, CA) and Anita Ramelius (Department of Clinical Sciences, Lund University/CRC, Skåne University Hospital, Malmö, Sweden) for their excellent technical assistance.

Funding. This work utilizes resources provided by the T1DGC, a collaborative clinical study sponsored by the National Institute of Diabetes and Digestive and Kidney Diseases, National Institute of Allergy and Infectious Diseases, National Human Genome Research Institute, National Institute of Child Health and Human Development, and JDRF and supported by U01 DK062418. This work was supported in part by American Diabetes Association grant 1-19-ICTS-076 (to J.P.M.), a grant from the UMass Memorial Healthcare JL Stock Fund (to J.P.M.), and National Institutes of Health grant AI 39095 (to J.P.M.).

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

Author Contributions. Ö.A. oversaw all MIP sequencing and data analysis and cowrote the manuscript. J.A.N. provided the T1DGC samples, codesigned the study, and cowrote the manuscript. J.A.B. had overall responsibility for MIP sequencing and data analysis and cowrote the manuscript. Å.L. provided the Swedish samples and cowrote the manuscript. P.M. performed the MIP sequencing. A.A.S. performed the DNA extraction of the Swedish BDD samples. F.B. performed the eQTL analyses of the published 1000 Genomes data. E.P.B. oversaw all data analyses, originated and supervised the eQTL analyses, and cowrote the manuscript. J.P.M. had overall responsibility for the entire project and cowrote the manuscript. The BDD Study Group Investigators (see appendix below) oversaw the collection and transmission of the BDD DNA sample set and gave permission for its use. J.P.M. 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.

Prior Presentation. Parts of this study were presented in abstract form at the Immunology of Diabetes Society Congress 2018, London, U.K., 25–29 October 2018.