Genetic Risk Score Modelling for Disease Progression in New-Onset Type 1 Diabetes Patients: Increased Genetic Load of Islet-Expressed and Cytokine-Regulated Candidate Genes Predicts Poorer Glycemic Control
Genome-wide association studies (GWAS) have identified over 40 type 1 diabetes risk loci. The clinical impact of these loci on β-cell function during disease progression is unknown. We aimed at testing whether a genetic risk score could predict glycemic control and residual β-cell function in type 1 diabetes (T1D). As gene expression may represent an intermediate phenotype between genetic variation and disease, we hypothesized that genes within T1D loci which are expressed in islets and transcriptionally regulated by proinflammatory cytokines would be the best predictors of disease progression. Two-thirds of 46 GWAS candidate genes examined were expressed in human islets, and 11 of these significantly changed expression levels following exposure to proinflammatory cytokines (IL-1β + IFNγ + TNFα) for 48 h. Using the GWAS single nucleotide polymorphisms (SNPs) from each locus, we constructed a genetic risk score based on the cumulative number of risk alleles carried in children with newly diagnosed T1D. With each additional risk allele carried, HbA1c levels increased significantly within first year after diagnosis. Network and gene ontology (GO) analyses revealed that several of the 11 candidate genes have overlapping biological functions and interact in a common network. Our results may help predict disease progression in newly diagnosed children with T1D which can be exploited for optimizing treatment.
In type 1 diabetes (T1D) the pancreatic β-cells are destroyed by the immune system in a process involving the proinflammatory cytokines interleukin-1-β (IL-1β), interferon-γ (IFNγ), and tumor necrosis factor-α (TNFα) released from antigen-presenting cells and T-cells [1, 2]. Genome-wide association scans (GWAS) have identified more than 40 genomic regions that are associated with T1D risk  (http://www.t1dbase.org). Many of the GWAS candidate genes have annotated immune-cell functions and most of the genetic risk variants have therefore been suggested to modulate immune-regulatory pathways [4, 5]. However, recent studies have highlighted that a significant proportion of the candidate genes are also expressed in human islets suggesting functional effects in β-cells [6–8] and possibly involvement in inflammation- and immune-mediated β-cell killing mechanisms thereby potentially affecting disease progression after clinical onset . As most variants identified through GWAS contribute to only modest effects to disease risk, it is likely that a combination of variants will better capture effects of clinical relevance. In T1D, very few studies have analyzed the impact of multiple variants on disease prediction and progression [10–12], although candidate gene-focused studies have demonstrated association with parameters of disease progression [13–16].
In the current study, we aimed at investigating whether a combined genetic risk score of T1D risk variants can predict glycemic control and residual β-cell function as assessed by HbA1c and insulin dose-adjusted HbA1c (IDAA1c) during disease progression in children with newly diagnosed T1D. We exclusively included SNPs for candidate genes expressed and transcriptionally regulated by cytokines in the target tissue of T1D, that is, human islets, as we hypothesized that these qualify as the most directly involved predictors.
2. Research Design and Methods
2.1. Expression Profiling of Candidate Genes in Human Islets
Human pancreatic islet preparations from nine nondiabetic donors (aged 8–57 years; 6 males and 3 females) were obtained from a multicenter European Union-supported program on β-cell transplantation in diabetes. None had classical T1D-associated HLA-DR risk genotypes. The program was approved by central and local ethical committees. Islet preparation, cytokine stimulation (5000 U/mL TNFα + 750 U/mL IFNγ + 75 U/mL IL-1β for 48 h), and RNA extraction have been described previously . Relative gene expression of candidate genes was evaluated by TaqMan assays using the Low Density Array system on TaqMan 7900HT (Applied Biosystems). Target gene expression was normalized to the geometric mean of three housekeeping genes (GAPDH, 18S-RNA, and PPIA) and evaluated using the delta-delta Ct method . One of the identified genes (IL10) whose expression was modulated following cytokine treatment was only detected in three of the human islet preparations. Genes with Ct values < 37 were considered as expressed.
2.2. Study Populations from the Hvidoere Study Group (HSG) on Childhood Diabetes
The study population was collected through HSG and is described in . The cohort included in total 257 children (126 girls and 131 boys). Eighty-four percent of the patients were white Caucasian, and age at clinical diagnosis was years (mean ± SEM), BMI kg/m2, and HbA1c at the time of diagnosis. DKA (HCO3 ≤ 15 mmol/L and/or pH ≤ 7.30) was present in 20.7% of the cases at the time of diagnosis. Exclusion criteria were suspected non-T1D (type 2 diabetes, maturity-onset diabetes of the young (MODY), or secondary diabetes), decline of enrolment into the study by patients or parents, and patients initially treated outside of the centers for more than 5 days. The diagnosis of T1D was according to the World Health Organization criteria. The study was performed according to the criteria of the Helsinki II Declaration and was approved by the local ethic committee in each center. All patients, their parents, or guardians gave informed consent. In the current study, patients with missing values for genotyping and clinical outcome measures were excluded leaving a total of 182 patients with complete genotype profile and clinical characterization.
2.2.1. HbA1c and IDAA1c (Insulin Dose-Adjusted HbA1c)
HbA1c was analyzed centrally by ion-exchange high-performance liquid chromatography at onset and 1, 3, 6, 9, and 12 months after diagnosis. IDAA1c is defined as actual HbA1c + (4 × insulin dose (U/Kg/24 h)). A calculated IDAA1c ≤ 9 corresponds to an estimated maximal C-peptide level above 300 pmol/L and has been used to define clinical remission .
Genotyping of rs2290400/GSDMB, rs2327832/TNFAIP3, rs4948088/COBL, rs7202877/CTRB1, rs7804356/SKAP2, rs1990760/IFIH1, rs3184504/SH2B3, rs6897932/IL7R, rs3024505/IL10, rs3825932/CTSH, and rs689/INS was done using the KASPar system (KBioscience, Hoddesdon, UK). Typing of the HLA-class II DRB1 locus was performed by direct sequencing of exon 2 of DRB1 according to Immuno Histocompatibility Working Group. The HLA risk groups were defined as high risk (DRB1 03/04, 04/04), moderate risk (DRB1 03/03, 04/08), and low risk (all other DRB1 genotype combinations).
2.4. Gene Ontology Terms and Network Construction
We used PANTHER  to perform functional annotation of the 11 input candidate genes. The enrichment for gene ontology (GO) terms in the biological process category was identified based on binomial test. The human genome was used as the reference list. To construct protein networks on the 11 input candidate genes, the STRING network tool was used. STRING is a database of known and predicted protein interaction data from multiple sources including experiments, coexpression, and text mining. In total, STRING covers nearly 10,000,000 proteins from over 2,000 organisms (http://string-db.org). Network was built with a medium confidence score (0.400) and up to 10 interactors.
2.5. Statistical Analysis
A genetic risk score was calculated for each individual based on the cumulative number of risk alleles carried for the 11 SNPs and was used as a continuous variable to test for association with IDAA1c and HbA1c levels at 1, 3, 6, 9, and 12 months after T1D onset in linear regression models. The assigned risk alleles for CTSH and SKAP2 were opposite compared to risk of T1D due to regression analyses from individual SNP models. Regression models were adjusted for the covariates sex, age group (0–5, 5–10, and >10 years at diagnosis), and HLA risk groups. Forward stepwise regression models were selected from all SNPs and covariates. A value below 0.05 was considered statistically significant. All statistical analyses were performed in SAS version 9.2.
3.1. Cytokine-Induced Gene Expression in Human Islets
Gene expression may represent an intermediate phenotype between genetic variation and disease. We therefore first evaluated the expression of established/pinpointed T1D GWAS candidate genes in human islets left untreated or exposed to a combination of proinflammatory cytokines (IL-1β + IFNγ + TNFα) for 48 hrs to mimic disease. We found 31 out of 46 tested genes to be expressed. Of these, 11 significantly changed their expression level following cytokine treatment () (Table 1). Six candidate genes were upregulated by cytokines, TNFAIP3, IFIH1, GSDMB, IL7R, IL10, and SH2B3, whereas 5 genes were downregulated, COBL, CTRB1, CTSH, SKAP2, and INS (Figure 1). Comparable expression profiles of these genes were observed in a recently published human islet dataset .
3.2. Genetic Risk Score Modelling of Glycemic Control and β-Cell Function
A genetic risk score model was constructed from the GWAS-identified SNPs linked to the 11 genes identified above to investigate the cumulative effect of T1D-associated risk alleles on disease progression in new-onset T1D children. The risk allele distribution is described in Supplementary Table in Supplementary Material available online at http://dx.doi.org/10.1155/2016/9570424. HbA1c and IDAA1c (a surrogate marker for β-cell function ) levels were increased in carriers with risk allele numbers at and above the 75th percentile (corresponding to minimum 15 risk alleles) during disease progression (HbA1c: 3, 6, and 9 months after disease onset, , , and , resp.; and IDAA1c: 9 months, ) (Figures 2(a) and 2(b)).
We then performed a multiple linear regression analysis adjusted for age, sex, and HLA risk groups and found significantly increased HbA1c and IDAA1c levels with increasing genetic risk score (GRS) from 3–12 months following T1D onset (Table 2). The validity of including GRS in the regression analysis was tested by comparing the variance explained by the model (). This clearly showed that including GRS as explanatory factor improved the model (Supplementary Table ). These findings suggest that residual β-cell function declines faster following diagnosis in patients carrying increased genetic load of islet-expressed and cytokine-regulated candidate genes.
3.3. Network and GO Analyses of Candidate Genes
We next asked if any of the 11 candidate genes may interact with each other in a functional protein network which could explain their cumulative effects on disease progression. This was evaluated by the STRING network tool which constructed a network that contained 7 out of the 11 genes (Figure 3). Consistent with this, the functional annotation of these candidate genes based on GO analyses revealed significantly enriched GO terms in biological processes category (Table 3). The 11 candidate genes were found enriched for various immune-mediated processes including regulation of immune response () and immune system process (). These findings support that several of the 11 candidate genes act in common networks and pathways to affect disease risk and progression.
Recent GWAS have identified a large number of loci affecting T1D risk . In this study, we investigated the clinical relevance of a genetic risk score on markers of disease progression. An increased genetic risk score associated with increasing HbA1c and IDAA1c levels the first year after disease onset, indicating that a higher genetic load of islet-expressed candidate genes predicts poorer glycemic control and residual β-cell function, respectively. One additional risk allele resulted in a 0.15% point increase in HbA1c after 12 months and a corresponding 0.19% increase in IDAA1c corresponding to a calculated 4% lower stimulated C-peptide . Our cumulative genetic risk score assumes that each risk variant contributes with equal effects to the traits, which probably does not reflect the true underlying biology. An alternative approach would be to weight each variant by published effect sizes for T1D risk, which has been done in type 2 diabetes [22, 23]. We chose the unweighted cumulative score because disease risk and disease progression are different outcomes, which will likely not have identical effect sizes. This is underlined by our previous observation that there is no statistically significant association between HLA risk and T1D progression . This is also the reason why we did not include HLA risk genes in the risk score model.
The observed poorer glycemic control associated with higher genetic load might prove to be a valuable tool for prediction of disease progression. This should, however, be validated in independent cohorts. An advantage of our study is that the inclusion of variants in the genetic risk score was based on prior “biological” knowledge, as we strictly focused on islet-expressed and cytokine-regulated candidate genes. Because regulated gene expression is often highly dynamic due to positive and negative feedback mechanisms, that is, the expression of a specific gene might be increased at one time point but decreased in another and vice versa, we did not take into account in the risk score model whether genes were up- or downregulated by cytokines but simply focused on the fact that their expression level changed as we considered this most important. We may have missed genes that changed expression at different time points compared to those examined at the 48 hrs, and a more detailed time-course study in human islets would likely have allowed a greater number of SNPs to be included in the risk score and thus provide even more accurate predictions.
Interestingly, we found that 7 of the 11 investigated genes interact in a protein network and several of the genes also shared GO terms suggesting that they affect the same biological mechanisms within the β-cells. We hypothesize that the genes are modulating β-cell function in terms of insulin secretion and/or the regenerative capacity and/or regulate the vulnerability of the β-cell to immune-mediated destruction. Indeed for some of the candidate genes, functional studies in β-cells have been performed. Hence, we recently demonstrated that CTSH regulates insulin gene transcription and secretion and also has antiapoptotic properties in β-cells . Similarly, A20, the protein name of the gene product encoded by TNFAIP3, is an antiapoptotic protein that inhibits apoptosis induced by cytokines by blocking activation of the transcription factor NFκB . In conclusion, a cumulative genetic risk score comprising variants from 11 islet-expressed candidate genes predicted significantly poorer glycemic control and β-cell function during disease progression in new-onset T1D children. This knowledge might be useful to better predict disease progression after diagnosis with T1D.
Conflict of Interests
All authors declare no conflict of interests.
Caroline A. Brorsson and Lotte B. Nielsen contributed equally to this paper. Lotte B. Nielsen, Caroline A. Brorsson, and Joachim Størling designed the study, researched data, and wrote the paper. Simranjeet Kaur, Marie Louise Andersen, Lars Hansen, and Regine Bergholdt researched data, contributed to discussion, and reviewed and edited the paper. Henrik B. Mortensen, Flemming Pociot, and Joachim Størling reviewed and edited the paper and contributed to discussion.
This study was funded by the Danish Council for Independent Research, the Sehested-Hansen Foundation, the National Institute of Health (1 DP3 DK085678), Beckett Fonden, Savværksejer Jeppe Juhl og hustru Ovita Juhls Mindelegat, European Foundation for the Study of Diabetes (EFSD), Novo Nordisk A/S, and Juvenile Diabetes Research Foundation (JDRF). The authors thank the members of the Hvidoere Study Group on Childhood Diabetes who have contributed to the Hvidoere Remission Phase Study: Henk-Jan Aanstoot, M.D., Ph.D., Center for Pediatric and Adolescent Diabetes Care and Research, Rotterdam, Netherlands; Carine de Beaufort, M.D., Clinique Pédiatrique, Luxembourg; Francesco Chiarelli, Professor M.D., Clinica Pediatrica, Chieti, Italy; Knut Dahl-Jørgensen, Professor, M.D., Dr. Med. SCI and Hilde Bjørndalen Göthner, M.D., Ullevål University Hospital, Department of Paediatrics, Oslo, Norway; Thomas Danne, Professor, M.D., Kinderkrankenhaus auf der Bult, Hannover, Germany; Patrick Garandeau, M.D., Unité D’endocrinologie Diabetologie Infantile, Institut Saint Pierre, France; Stephen A. Greene, M.D., DC, University of Dundee, Scotland; Hilary Hoey, Professor, M.D., FRCPI, University of Dublin, National Children’s Hospital, Tallaght, Ireland; Reinhard W. Holl, Professor M.D., University of Ulm, Germany; Mirjana Kocova, Professor, M.D., Pediatric Clinic-Skopje, Republic of Macedonia; Pedro Martul, Professor M.D., Ph.D., Endocrinologia Pediatrica Hospital de Cruces, Spain; Nobuo Matsuura, Professor, M.D., Kitasato University School of Medicine, Japan; Henrik B. Mortensen, Professor, M.D., DMSc, Department of Pediatrics, Herlev Hospital, Faculty of Health and Medical Sciences, University of Copenhagen, Denmark, Kenneth J. Robertson, M.D., Royal Hospital for Sick Children, Yorkhill, Glasgow, Scotland; Eugen J. Schoenle, Professor, M.D., University Children’s Hospital, Zurich, Switzerland; Peter Swift, M.D., Leicester Royal Infirmary Childrens Hospital, Leicester, UK; Rosa Maria Tsou, M.D./Professor, Manuel Fontoura, Paediatric Department Oporto, Portugal; Maurizio Vanelli, Professor, M.D., Clinica Pediatrica, Centro di Diabetologia, University of Parma, Italy; Jan Åman, M.D., Ph.D., Örebro Medical Centre Hospital, Department of Paediatrics, Sweden.
Supplementary Table 1: The risk allele distribution of the 11 T1D candidate genes.
Supplementary Table 2.: Variance explained by regression models for HbA1c and IDAA1c with and without genetic risk score (GRS).
R. Bergholdt, C. Brorsson, A. Palleja et al., “Identification of novel type 1 diabetes candidate genes by integrating genome-wide association data, protein-protein interactions, and human pancreatic islet gene expression,” Diabetes, vol. 61, no. 4, pp. 954–962, 2012.View at: Publisher Site | Google Scholar
M. L. M. Andersen, M. A. Rasmussen, S. Pörksen et al., “Complex multi-block analysis identifies new immunologic and genetic disease progression patterns associated with the residual beta-cell function 1 year after diagnosis of type 1 diabetes,” PLoS ONE, vol. 8, no. 6, Article ID e64632, 2013.View at: Publisher Site | Google Scholar
L. B. Nielsen, H. B. Mortensen, F. Chiarelli et al., “Impact of IDDM2 on disease pathogenesis and progression in children with newly diagnosed type 1 diabetes: reduced insulin antibody titres and preserved beta cell function,” Diabetologia, vol. 49, no. 1, pp. 71–74, 2006.View at: Publisher Site | Google Scholar
A. Petrone, M. Spoletini, S. Zampetti et al., “The PTPN22 1858T gene variant in type 1 diabetes is associated with reduced residual beta-cell function and worse metabolic control,” Diabetes Care, vol. 31, pp. 1214–1218, 2008.View at: Google Scholar
T. Flyøel, C. Brorsson, L. B. Nielsen et al., “CTSH regulates β-cell function and disease progression in newly diagnosed type 1 diabetes patients,” Proceedings of the National Academy of Sciences of the United States of America, vol. 111, no. 28, pp. 10305–10310, 2014.View at: Publisher Site | Google Scholar
H. B. Mortensen, P. G. Swift, R. W. Holl et al., “Multinational study in children and adolescents with newly diagnosed type 1 diabetes: association of age, ketoacidosis, HLA status, and autoantibodies on residual beta-cell function and glycemic control 12 months after diagnosis,” Pediatric Diabetes, vol. 11, no. 4, pp. 218–226, 2010.View at: Publisher Site | Google Scholar
E. A. Andersson, K. H. Allin, C. H. Sandholt et al., “Genetic risk score of 46 type 2 diabetes risk variants associates with changes in plasma glucose and estimates of pancreatic β-cell function over 5 years of follow-up,” Diabetes, vol. 62, no. 10, pp. 3610–3617, 2013.View at: Publisher Site | Google Scholar