Telomere length is not a main factor for the development of islet autoimmunity and type 1 diabetes in the TEDDY study

The Environmental Determinants of Diabetes in the Young (TEDDY) study enrolled 8676 children, 3–4 months of age, born with HLA-susceptibility genotypes for islet autoimmunity (IA) and type 1 diabetes (T1D). Whole-genome sequencing (WGS) was performed in 1119 children in a nested case–control study design. Telomere length was estimated from WGS data using five tools: Computel, Telseq, Telomerecat, qMotif and Motif_counter. The estimated median telomere length was 5.10 kb (IQR 4.52–5.68 kb) using Computel. The age when the blood sample was drawn had a significant negative correlation with telomere length (P = 0.003). European children, particularly those from Finland (P = 0.041) and from Sweden (P = 0.001), had shorter telomeres than children from the U.S.A. Paternal age (P = 0.019) was positively associated with telomere length. First-degree relative status, presence of gestational diabetes in the mother, and maternal age did not have a significant impact on estimated telomere length. HLA-DR4/4 or HLA-DR4/X children had significantly longer telomeres compared to children with HLA-DR3/3 or HLA-DR3/9 haplogenotypes (P = 0.008). Estimated telomere length was not significantly different with respect to any IA (P = 0.377), IAA-first (P = 0.248), GADA-first (P = 0.248) or T1D (P = 0.861). These results suggest that telomere length has no major impact on the risk for IA, the first step to develop T1D. Nevertheless, telomere length was shorter in the T1D high prevalence populations, Finland and Sweden.

www.nature.com/scientificreports/ Telomere length is, in part, inherited with a stronger influence on variation in telomere length from the father, defined by the stronger correlation in telomere length between father-son and father-daughter 3,4 . In addition, paternal age is positively correlated with longer telomeres 5 . At the population level, there is no difference in telomere length at birth between boys and girls 6 ; however, in young children 7 and young adults (around 30 years), telomeres are longer in women than in men, although the attrition rate is higher in women 8 . The preservation of telomeres in women may, therefore, occur during the first two decades of life 9 .
In healthy subjects, HLA-DRB1*04 alleles have been shown to be associated with shorter telomeres in CD4 + T-cells, although the difference in telomere length between HLA-DR4 + and HLA-DR4subjects was not identified at birth. Thus, during the first 20 years of life, telomere attrition may be accelerated in subjects with HLA-DR4 +10 .
Several mechanisms may contribute to the attrition of telomeres, but there is only one known mechanism in healthy individuals, beyond embryogenesis, that counteracts the shortening of telomeres. Some cell types, such as lymphocytes, are capable of activating telomerase, an enzyme that can elongate telomere sequences and thereby modulate cellular lifespan. Certain hormones can also stimulate telomerases. For example, estrogens have been shown to activate telomerase activity 11 while cortisol can decrease telomerase activity. In vitro experiments have shown that high concentrations of hydrocortisone (comparable to cortisol levels that may occur in vivo due to stress) can reduce telomerase activity by up to 50% 12 . Women with gestational diabetes (GDM), a condition of altered glucose metabolism and hormonal patterns during pregnancy, is associated with shorter telomeres in girls born to GDM mothers, but not in boys 13 .
Telomere length exhibits disease-specific patterns in different autoimmune diseases 2,14,15 . Shorter telomere length in leukocytes has been reported in subjects with type 1 diabetes (T1D) compared to controls, although shorter telomere length was not associated with the duration of T1D 16 . In T1D, all-cause mortality is also associated with a shorter telomere length 17 .
Here, we utilize the TEDDY nested case-control (NCC) cohort to estimate telomere length from wholegenome sequencing (WGS) data. Analyses of telomere length were performed to find out if there was a relation to the development of islet autoimmunity (IA), type 1 diabetes (T1D), or both, as well as its association with HLA-DR-DQ haplogenotypes. Primary end-points in TEDDY included any IA, micro insulin autoantibody (mIAA)-first, the 65 kDa glutamate decarboxylase autoantibody (GADA)-first, or T1D.
The first blood sample for HLA-screening was obtained either as cord blood in the maternity clinic or as a dry blood spot (DBS) on day three or four. If the child was eligible for TEDDY, the family was contacted by a study nurse and invited to participate in the 15-year follow-up study with blood sampling for determining islet autoantibody (GADA, islet antigen-2 antibody (IA-2A) and mIAA) status every three months between 3 and 48 months of child's age and every six months thereafter (while a child with islet autoantibody positivity remained on the three months visit schedule). HLA-haplogenotypes were confirmed in a second blood sample at the 9-month visit. For the present study, blood samples were taken at the median age of 0.82 years (interquartile range (IQR) 0.75-5.43 years) for extraction of DNA and subsequent sequencing.
All procedures performed in studies involving human participants were in accordance with the ethical standards of the institutional and/or national research committees and with the 1964 Helsinki declaration and its later amendments or comparable ethical standards. Study outcomes: IA and T1D. The primary outcome of the TEDDY study is the development of persistent confirmed IA assessed every three months. In the U.S.A., all sera were assayed at the Barbara Davis Center for Childhood Diabetes at the University of Colorado, Denver, CO; in Europe all sera were assayed at the University of Bristol, Bristol, U.K. GADA, IA-2A and mIAA were all measured using radio-binding assays [19][20][21] . IA was confirmed if an autoantibody was identified in a sample at both Reference Laboratories. All samples positive for IA and 5% of negative samples were re-tested for confirmation in both Reference Laboratories. www.nature.com/scientificreports/ Persistent autoimmunity was defined by the presence of confirmed IA (GADA, IA-2A or mIAA) on two or more consecutive visits. The date of persistent confirmed IA was defined as the first blood draw date at which the child was found positive for at least one islet autoantibody. As children may be born with maternal IA, positive results due to maternal transmission via placenta were excluded when defining the child's IA status. If the child was IA-positive at 3 or 6 months of age, the IA status of the mother was assessed at the 9-month visit in order to distinguish maternal autoantibody from IA in the child. If maternal autoantibody was present, the child was not considered persistently IA-positive unless the child had a negative sample prior to the first positive sample.
The threshold for positivity differed by the laboratory. Positivity was defined by Bristol (GADA ≥ 33 DK-U/ mL, IA-2A ≥ 5 DK-U/mL and mIAA > 0.95 Local Units) and Denver (GADA > 20 DK-U/mL, IA-2A > 5 DK-U/ mL and mIAA > 0.010 index). Both Reference Laboratories have shown high sensitivity and specificity as well as concordance for defining the presence of confirmed IA 10,21 .
HLA typing. HLA genotype screening was performed using either a dried blood spot (DBS) punch or a small volume blood specimen format 18,22 . After polymerase chain reaction (PCR) amplification of exon 2 of the HLA class II genes (DRB1, DQA1 or DQB1), alleles were identified either by direct sequencing, oligonucleotide probe hybridization or other genotyping techniques. Typing to certify specific HLA-DR-DQ haplogenotypes was determined for each clinical center. Confirmation of the HLA haplogenotypes was performed by the central HLA Reference Laboratory at Roche Molecular Systems, Oakland, CA on the eligible subjects using the 9-month sample.
Estimation of telomere length from the WGS data. WGS was performed on the subjects in two nested case-control (NCC) studies: one for IA and the other for T1D (Table 1). In the NCC study, the "cases" were those subjects that had developed any IA or T1D as of 30 June 2018 and "controls" were randomly selected from those autoantibody-negative or non-T1D children in the index case's risk set, matched for sex, clinical site and first-degree relative (FDR) status 23 . Matching was performed in a 1:1 format (one control was matched to each case).
Extraction of DNA was performed by the Center for Public Health Genomics at the University of Virginia, VA, U.S.A. The WGS data were generated on the Illumina HiSeq X Series platform with paired-end 2 × 150 bp reads with targeted 30× coverage by Macrogen, Inc, Rockville, MD, U.S.A. FastQC, MultiQC 24 and FastQ Screen 25 were used to identify low-quality sequences and contaminants in the raw sequence data. The raw sequences were aligned to the Genome Reference Consortium Build 37 (GRCh37d5) for qMotif and to the Genome Reference Consortium Build 38 (GRCh38DH) for Computel, Telseq, Telomerecat and Motif_counter. Post-alignment processing was performed by utilizing the University of Michigan's (UM) Docker-alignment pipeline (https:// github. com/ statg en/ docker-align ment) based on Burrows-Wheeler Aligner (BWA) 26 . Variant concordance analysis was performed by comparing the variants identified from the WGS data to the T1DExomeChip array (Illumina, CA) data to determine mislabeled, swapped or contaminated samples.
Estimation of the telomere length from the WGS data was performed using five different tools: Computel, Telseq, Telomerecat, qMotif and Motif_counter (Supplementary methods). Pairwise Spearman correlation analysis was used to compare performance of these tools against each other. All subsequent statistical analysis on telomere length was conducted using estimates from Computel which is based on aligning raw sequence data in FASTQ format to a telomeric reference sequence. Table 1. Length of telomeres (kb) in study subjects by the status of islet autoimmunity (IA) and type 1 diabetes (T1D) in the Environmental Determinants of Diabetes in the Young (TEDDY) study. Telomere length data were available on both case and matched control subjects in 389 IA pairs and 118 T1D pairs. Computel was used to estimate the telomere length from the whole-genome sequence (WGS) data. Telomere length is presented as mean ± standard deviation (S.D.) in kb.

T1D controls (in kb)
Country  (Fig. 1). The T1DExomeChip is a custom genotyping array with more than 90,000 custom content SNPs added to the Infinium® CoreExome-24 v.1.1 BeadChip (Illumina, CA). Genotype calling and the following quality control steps were applied to the dataset: (1) individuals with low call rate (< 95%), or discordance with reported sex and prior genotyping were not considered in the analysis, (2) SNPs with low call rates (< 95%) were excluded, (3) SNPs from autosomal chromosomes with allele distributions strongly deviating from Hardy-Weinberg equilibrium in controls with the European ancestry (P < 10 -6 ) were discarded (except for chromosome 6 with the Hardy-Weinberg equilibrium in controls with the European ancestry (P < 10 -20 ) due to HLA eligibility requirements), (4) SNPs on chromosome X were excluded with heterozygosity > 1% in men, (5) SNPs on chromosome Y were excluded with genotypes > 2% in women, (6) Monomorphic SNPs were excluded, (7) SNPs were discarded for which the concordance rate among duplicates was < 99% and (8)  Statistical analyses. Multiple linear regression analysis was performed to examine the association between the telomere length and other factors including HLA-categories (HLA-DR3/3 or HLA-DR3/9 (as reference), HLA-DR3/4 and HLA-DR4/4, HLA-DR4/X), sex, parental age, mother´s gestational diabetes status, FDR status and the country of residence. The age at sample collection and the indicator of the long-distance protocol sample collection (yes vs. no) were also included in the model. Potential variation introduced by batches was taken account by a random intercept. Telomere length data were available on a total of 1119 TEDDY children. Primary analyses were focused on examining the association between telomere length with the risk of any IA and/or T1D. A conditional logistic regression model was used to consider the matching criteria in the NCC study for any IA or T1D, adjusting for the HLA-categories (HLA-DR3/3 or HLA-DR3/9 (as reference), HLA-DR3/4 and HLA-DR4/4, HLA-DR4/X). In the analysis, telomere length was age-adjusted by using the residuals from the linear regression of telomere length on the age at sample collection. Subgroup analyses were conducted with respect to the type of first appearing autoantibody (mIAA-first and GADA-first). The magnitude of the association was estimated by odds ratio (OR), with 95% confidence internal (CI).
Secondary analyses were performed to examine the effect of SNPs associated with telomere length on the risk of any IA, mIAA-first, GADA-first and T1D (Supplementary Table S2). A total of 236 SNPs with minor allele frequency (MAF) > 0.05 were included (Fig. 1). Cox proportional hazard models were used to study the Whole-genome sequencing (WGS) was performed on the subjects in the nested case-control (NCC) study for islet autoimmunity (IA) and type 1 diabetes (T1D). In the NCC study, the "cases" were those subjects that had developed any IA (n = 389) or T1D (n = 118) and "controls" were randomly selected from those autoantibody-negative or non-T1D children in the index case's risk set, matched for sex, clinical site and first-degree relative (FDR) status. Telomere length estimation was performed using five different tools. Single nucleotide polymorphisms (SNPs) were genotyped using the T1DExomeChip array. A total of 835 SNPs associated with telomere length were available on the T1DExomeChip and 236 of those had a minor allele frequency (MAF) > 0.05. SNP analysis was performed on 8093 TEDDY children. www.nature.com/scientificreports/ association between each SNP and the risk of any IA, mIAA-first, GADA-first and T1D, respectively, adjusting for sex, the country of residence, FDR, HLA-haplogenotype and accounting for population stratification (ancestral heterogeneity) (using the top three principal components estimated from the T1DExomeChip). Of the 8676 enrolled children, 583 were excluded: indeterminate autoantibody statuses (n = 14), HLA ineligible (n = 112), and no SNP data available (n = 457). A total of 7093 families had one child, 476 families had 2 children, and 16 families had 3 children who were included in the analysis. A robust variance estimate was used to account for the dependence within families in the Cox proportional hazard model 39 . The magnitude of each estimated association was defined by the Hazard ratio (HR), with 95% confidence intervals (CIs). To correct for multiple hypothesis testing, the false discovery rate adjusted P-values were calculated 40 . SNPs with nominal significance can be found in Supplementary Table S4 through S7. Statistical analyses were performed using Statistical Analysis Software (Version 9.4, SAS Institute, Cary, North Carolina, U.S.A.) and R (Version 4.0.2). For the multiple linear regression analysis, estimates with P < 0.05 were considered as statistically significant.

Results
Telomere lengths from the WGS data were calculated from 1119 TEDDY children (Fig. 1). Descriptive data on the estimation of the telomere lengths is shown by the status of IA and T1D in Table 1. The estimated median telomere length was 5.10 kb (IQR 4.52-5.68 kb) using Computel. The estimations of the telomere lengths were strongly correlated between Computel and Telseq (r = 0.97, P < 0.001) as well as between Computel and Motif_ counter (r = 0.93, P < 0.001) (Supplemental Fig. S1).
Association between telomere length and other phenotypes. The age of the child when the sample was drawn for WGS had a significant impact on the length of telomeres. There was a significant negative association between age of the children and telomere length (P = 0.003) ( Table 2). Children that carried HLA-DR4/4 or HLA-DR4/X had significantly longer telomeres compared to children with HLA-DR3/3 or HLA-DR3/9 haplogenotypes (P = 0.008). European children, particularly those from Finland (P = 0.041) and from Sweden (P = 0.001), had shorter telomeres than children from the U.S.A. Though the paternal age (in years) did not differ between the countries (Analysis of variance (ANOVA) P = 0.102) as shown in Supplementary Table S3. Paternal age (P = 0.019) and the female sex of the child (P = 0.069; although not significant) were positively associated with telomere length. FDR status, presence of GDM in the mother, and maternal age did not have a significant impact on estimated telomere length ( Table 2).
SNPs associated with telomere length and the risk of any IA, mIAA-first, GADA-first and T1D. We selected 236 SNPs previously reported to be associated with telomere length (see Material and methods). A total of 31 SNPs in 13 regions were nominally associated with risk of developing any IA, mIAA-first, GADA-first or T1D; however, none of these associations withstood correction for multiple hypothesis testing (Supplementary Table S4

Discussion
The key finding in this study was that children from high incidence countries for T1D i.e., Finland and especially from Sweden had shorter telomeres in white blood cells as compared to children from other TEDDY sites (Germany and three sites in the U.S.A.). The age when the blood sample was drawn had a significant negative correlation with telomere length, consistent with previous reports 41 . Interestingly, children that carried HLA-DR4/4 or HLA-DR4/X had significantly longer telomeres compared to children with HLA-DR3/3 or HLA-DR3/9 haplogenotypes. Also, paternal age was significantly correlated with telomere length. Children with HLA-DR3/4 had a significantly increased risk to develop any IA as compared to children with HLA-DR3/3 or HLA-DR3/9 haplogenotypes. Moreover, there were significantly increased risks for children with HLA-DR3/4 www.nature.com/scientificreports/ and HLA-DR4/4 or HLA-DR4/X to develop mIAA as the first autoantibody as compared to those with HLA-DR3/3 or HLA-DR3/9 haplogenotypes. In contrast, GADA as the first autoantibody was not associated with any of the HLA-DR-DQ haplogenotypes. To our knowledge, there has been no reports on the risk for IA, mIAA-first or GADA-first in relation to telomere length. The estimations of the telomeres were highly correlated between Computel, Telseq and Motif_ counter with R 2 values ranging from 0.85 to 0.94 (Supplemental Fig. 1). The finding of similar telomere lengths in young boys and girls confirms a previous study of telomere length taken at birth 6 . There has been a previous report on shorter telomere length in children with T1D 16 , a finding that was not reproduced here. This discrepancy could be explained by the older age of the participants in the previous report (17-48 years), compared with the TEDDY study in which the oldest children with T1D were 14.16 years. We could not reproduce the findings of shorter telomeres in children born to mothers with GDM, although this could be due to the older age (9-16 years) in the earlier study 13 . In addition, we detected a significant effect of paternal age but not maternal age on the length of telomeres 5 . The SNPs that had nominally significant association to our primary and secondary outcomes-any IA, mIAA-first, GADA-first or T1D-were mainly from genes coding for proteins involved in double DNA strand repair and maintenance of telomeres (MRE11 and RAD50) 42 , DNA replication and repair (BLM) 43 , genome stability (ATM) 44 or degradative endocytic trafficking (PIK3C3) 45 . The BLM gene has gene ontology including nucleic acid binding and ATPase activity also involved in DNA replication and repair.
The major strength of the present study is the use of both IA and T1D as endpoints. Additional strengths of this study are the inclusion of unrelated subjects from both the general population and FDRs without islet autoantibodies or T1D from the same countries in the same age ranges as the subjects with outcome events. One limitation of the present study is the relatively low number of subjects with IA, particularly those with GADAfirst or with T1D, raising concerns for statistical power of the association analysis. Nevertheless, the screening of more than 400,000 newborns in four different countries to identify those with an increased genetic risk for T1D, defined by high-risk HLA-DR-DQ haplogenotypes, followed with blood draws every three months for the first four years and every six months thereafter is a unique effort. The proportional hazard ratio identifies risk for time-to-event in this longitudinal study, thereby increasing the power compared with cross-sectional case-control study design.
We concluded that genes in the HLA-region particularly HLA-DR3/4 and HLA-DR4/4, HLA-DR4/X remain the most important genetic risk factors for IA, the first step towards T1D. Children with the high-risk genotypes for T1D from Finland and Sweden had shorter telomeres and telomere length is associated with HLA-DR4/4 or HLA-DR4/X haplogenotype as compared to HLA-DR3/3 or HLA-DR3/9. Moreover, paternal age was positively associated with telomere length.

Data availability
The datasets generated and analyzed during the current study will be made available in the NIDDK Central Repository at https:// repos itory. niddk. nih. gov/ studi es/ teddy/. The TEDDY Whole Genome Sequencing (WGS) data that support the findings of this study have been deposited in NCBI's database of Genotypes and Phenotypes (dbGaP) with the primary accession code phs001442.v3.p2.