Associations of GWAS-Supported Non-MHC Genes with Autoimmune Thyroiditis in Patients with Type 1 Diabetes

Introduction

Type 1 diabetes (T1D) is characterized by a severe autoimmune pancreatic β-cell destruction, usually leading to absolute insulin deficiency.¹ In clinical practice, this disorder is closely associated with autoimmune thyroid diseases (AITD), another organ-specific T-cell mediated disease. Epidemiological data have shown that T1D and AITD frequently coexist in the same individual, who is also classified as a variant of autoimmune polyglandular syndrome type 3 (APS3v).² A higher prevalence of AITD has been reported in T1D compared with those without T1D. In different studies, AITD occurs in 17% to 35% of T1D patients.^2–7 A more recent study declared that lifetime risk of AITD in longstanding T1D is about 41% and is greater in females.⁸

The underlying mechanisms for the concurrent development of AITD in T1D patients are still unclear. A possible explanation lies in the common genetic background in conjunction with environmental factors.⁹ In the past years, a genome-wide association study (GWAS) in Caucasian population identified multiple signals with the HLA region, as well as 5 non-MHC genes (PHTF1, PTPN22, MAGI3, BCL2L15, and QRFPR) associated with the joint risk for AITD and T1D. Notably, the novel QRFPR gene was unique to the APS3v phenotype that had not been suggested by previous studies of T1D or AITD.¹⁰ However, the allele frequency distribution of variant in QRFPR was quite different from different populations. According to the 1000 Genomes Browser (http://www.ncbi.nlm.nih.gov/variation/tools/1000genomes/), unlike Europeans with minor allele “A”, East Asians have a lower allele frequency of allele “G”. Because of possible genetic heterogeneity among different ethnic groups, it is essential to replicate this association in Chinese population.

Moreover, female is an independent risk factor for a large number of autoimmune diseases. It has been reported that interactive effects of sex and gene determinants may synergistically influence the risk of autoimmune diseases.^11–14 Here, we examined the potential associations between polymorphisms in GWAS-supported non-MHC genes and risk of AITD in patients with T1D in Chinese Han population. The interactive effects of genotypes and sex were also analyzed.

Materials and Methods

Study Subjects and Data Collection

A case–control study was designed. We consecutively recruited unrelated Chinese Han patients with diagnosed T1D from the Department of Endocrinology and Metabolism in Zhongshan People’s Hospital between October 2017 and March 2020. T1D was diagnosed based on the American Diabetes Association criteria.¹ Those diagnosed with T1D over the age of 30 were defined as late-onset T1D. Of all T1D patients without known thyroid disorder were screened for AITD. AITD was defined as the presence of elevated thyroid antiperoxidase antibody (TPOAb) and/or antithyroglobulin antibody (TgAb) values and of the typical pattern of the thyroid ultrasound.¹⁵ All of the subjects with AITD were newly diagnosed and previously untreated. Participants were then divided into two groups: cases were defined as T1D patients with AITD and controls were defined as T1D patients without AITD.

Trained physicians collected clinical and biochemical characteristics. Blood samples were taken after 12h overnight fasting for laboratory investigations. Fasting blood glucose, fasting C-peptide, glycated hemoglobin, glutamic acid decarboxylase autoantibody (GADA), lipid concentrations, and creatinine were monitored.

This study was approved by the Ethical Review Board of Zhongshan People’s Hospital, Zhongshan, China (Approval number K2017-06). The study was performed in accordance with the ethical standards laid down in the 1964 Declaration of Helsinki. Participants ≥18 years old or a parent/guardian of participants <18 years old provided written informed consent.

Genotyping

Variants were selected based on a GWAS for co-occurrence of T1D and AITD in Caucasians.¹⁰ The GWAS associated 5 non-MHC genes (PHTF1, PTPN22, MAGI3, BCL2L15, and QRFPR) with co-occurrence of T1D and AITD. Among them, rs2476601 in PTPN22, which reached genome wide significance (P < 5×10⁻⁸) in Caucasians, showed no polymorphism in Chinese population. Therefore, we instead included rs1217407 in PTPN22 that close to the genome-wide significance (P = 1.87×10⁻⁶). As a result, a total of 5 single-nucleotide polymorphisms (SNPs), 4 located on chromosome 1p13 (PHTF1 rs1111695, PTPN22 rs1217407, MAGI3 rs2153977, and BCL2L15 rs2358994) and 1 located on chromosome 4q27 (QRFPR rs7679475), were included in this study. Supplementary Table 1 shows details of the selected SNPs.

Genotyping was conducted by Improved Multiple Ligase Detection Reaction (iMLDR)¹⁶ according to the standard procedures, with technical assistance from Center for Genetic and Genomic Analysis, Genesky Biotechnologies Inc. (Shanghai, China). Genotyping was successfully performed in all samples. For quality control, 5% of the samples were randomly selected for reproducibility assessments and the results were 100% concordant.

Statistical Analyses

Departures from Hardy–Weinberg equilibrium were tested using chi-square tests. Categorical variables were presented as count (percentage) and continuous variables as median (interquartile range). Chi-square tests or Fisher’s exact was conducted for categorical data as appropriate, while Mann–Whitney U-tests for continuous variables. Univariate and multivariate logistic regression were performed to investigate the association of each individual polymorphism with risk of AITD in T1D patients. Variables with a marginal significant difference (P < 0.1) in baseline comparison or with an important biological association (age, sex, and BMI) were included as covariates in the multivariate model. Sex-genotype interaction, linkage disequilibrium (LD) analysis, and haplotype analysis were analyzed by the SNPStats web tool.¹⁷ Haplotypes with frequencies <0.01 were excluded. Multiple testing corrections were carried out with the Bonferroni method. A P value <0.01 (0.05/5 for five SNPs) was considered as statistically significant. Statistical analysis was performed by R version 3.5.3.

In silico Analysis

To provide function annotation for QRFPR polymorphism rs7679475, we performed in silico analysis using Haploreg v4.1¹⁸ (https://pubs.broadinstitute.org/mammals/haploreg/haploreg.php) and UCSC Genome Browser (http://genome.ucsc.edu/).

Results

Baseline Characteristics

In total, 489 T1D patients were enrolled in the study, in which 159 (32.5%) patients with AITD. The baseline characteristics of the study population are summarized in Table 1. Patients with AITD were older (P = 0.004) and had a higher percentage of females (71.7% vs 59.4%, P = 0.011), higher proportion of late-onset T1D (52.8% vs 42.4%, P = 0.031), higher fasting C-peptide level (P = 0.010), higher positivity to GADA (65.4% vs 46.4%, P < 0.001), lower low-density lipoprotein (LDL-C) level (P = 0.011), and higher high-density lipoprotein (HDL-C) level (P < 0.001) compared with patients without AITD. No statistically significant difference in other factors between the two groups was observed.

Table 1 Baseline Characteristic of Total Population, T1D Patients with AITD (Cases), and T1D Patients without AITD (Controls)

Genotype and Allele Distribution

Genotype distributions were not deviated from Hardy–Weinberg equilibrium for all SNPs (Supplementary Table 1). Genotype frequencies of the five investigated variants in T1D cases with AITD and T1D controls without AITD are shown in Table 2. QRFPR rs7679475G allele frequency was 39.3% in cases with AITD and 30.2% in controls (P = 0.004). The distribution of genotype of rs7679475 differed significantly between case and control groups (A/A, G/A, and G/G were 37.1%, 47.2%, and 15.7%, respectively, in cases vs 46.7%, 46.4%, and 7.0%, respectively, in controls; P = 0.005). Univariate logistic regression suggested that G/G genotype of rs7679475 conferred risk for AITD in T1D individuals (codominant model: OR 2.84; 95% CI 1.49–5.39; P = 0.001; recessive model: OR 2.49; 95% CI 1.36–4.54; P = 0.003). Upon adjusted for age, sex, and BMI, the results were remained significant. Moreover, after further adjustment for variables with a marginal significant difference (P < 0.1) in baseline comparison (late-onset T1D, fasting C-peptide, positivity to GADA, LDL-C, and HDL-C), rs7679475G was still associated with AITD in patients with T1D in codominant model (G/G vs A/A, OR 2.93; 95% CI 1.44–5.96; P = 0.003) and dominant model (G/A-G/G vs A/A, OR 1.81; 95% CI 1.17–2.79; P = 0.007). In recessive model, though not reaching the Bonferroni-corrected threshold, there was a suggestive association between rs7679475 and AITD in T1D patients (G/G vs A/A-G/A, OR 2.28; 95% CI 1.17–4.43; P = 0.015).

Table 2 Genotype and Allele Frequencies of the Five Selected SNPs in T1D Patients with and without AITD

However, none of the other SNPs (rs1111695, rs1217407, rs2153977, and rs2358994) was significantly associated with the presence of AITD in T1D patients, and multivariate adjustment did not change these results.

Sex Interaction Analysis

Since female is widely acknowledged to be a risk factor for a large amount of autoimmune diseases, we performed an interaction analysis between rs7679475 variant and sex. Figure 1 shows the effect of rs7679475 on the risk of AITD in T1D individuals stratified by sex after adjusting for multiple confounders. Compared with rs7679475 A/A genotype in men, T1D women with G/A had a 5.11-fold (95% CI 2.23–11.71) increased risk of AITD, while males with G/G had an OR of 6.43 (95% CI 1.79–23.12). The highest risk of AITD in T1D patients was found in females with G/G genotype (OR 11.63; 95% CI 3.94–34.34). Moreover, interaction between rs7679475 and sex was significant (P _interaction = 0.005).

Figure 1 Adjusted odd ratio (OR) of autoimmune thyroid diseases (AITD) in type 1 diabetes (T1D) patients according rs7679475 and sex, male with A/A genotype as a reference (adjusted for age, sex, BMI, proportion of late-onset T1D, fasting C-peptide, glutamic acid decarboxylase autoantibody, low-density lipoprotein, and high-density lipoprotein). The interaction between rs7679475 and sex was significant (P interaction = 0.005).

Haplotype Analysis

Four of the five target SNPs are located in chromosome 1p13 locus. LD analysis showed that MAGI3 rs2153977, PHTF1 rs1111695, and PTPN22 rs1217407 were in an LD block (Supplementary Figure 1). In haplotype analysis, after adjusting for potential confounders, no significant association between any haplotype and the risk of AITD in patients with T1D (Supplementary Table 2) was observed (P > 0.05).

Functional Annotation

We carried out in silico analysis to evaluate the potential function of the significant QRFPR polymorphism. Rs7679475 is located in the 11859bp upstream of QRFPR and in H3K27ac histone modification marked region, a marker of enhancer/promoter-associated regulatory element (Figure 2). Further bioinformatics analysis by HaploReg v4.1 indicated that rs7679475 could change the binding affinity of regulatory motifs Brachyury_2, PPAR_1, and SIX5_known1 (Table 3).

Table 3 Regulatory Motifs Altered for Rs7679475 Based on HaploReg V4.1

Figure 2 Evidence from UCSC data for regulatory function of rs7679475. The plot is presented on chromosome 4q27 (NCBI GRCh37). Tracks (from top to bottom) in the picture are genome base position, chromosome bands, UCSC genes, protein interactions from curated databases and text-mining, human mRNAs from GenBank, gene expression, H3K27Ac Mark, DNaseI Hypersensitivity Clusters, ENCODE transcription factor ChIP-seq clusters, ENCODE chromatin state segmentation by HMM, and simple nucleotide polymorphisms (dbSNP 151).

Discussion

In the current study, we comprehensively evaluated the associations of five GWAS-supported non-MHC genes polymorphisms and AITD in patients with T1D in a Chinese Han population. The main finding was that QRFPR rs7679475G was associated with an increased risk of AITD in T1D individuals. Moreover, we observed significant interaction between QRFPR genotype and sex.

Thyroid dysfunction and diabetes mellitus, two closely associated disorders, share common underlying pathophysiological mechanisms.¹⁹ Shared genetic susceptibility might contribute to the co-occurrence of AITD and T1D. In 2015, a GWAS in Caucasian reported that 5 non-MHC genes were associated with the co-occurrence of AITD and T1D at approximate genome-wide significance.¹⁰ A subsequent replicate study in a Japanese population only found a significant association between rs7679475 and Hashimoto’s thyroiditis.²⁰ Concordantly, our study demonstrated that rs7679475 was an independent risk factor for AITD in those diagnosed T1D in Chinese population, but did not show significant association with other SNPs that identified by the Caucasian GWAS. The inconsistency may be partly explained by genetic heterogeneity in different ethnic origins and regions. According to the 1000 Genomes Browser, the allele frequencies of rs7679475 are quite different between Chinese and European (CHS: A = 0.6476, G = 0.3524 and CEU: A = 0.4091, G = 0.5909), suggesting that ethnic differences may exert different influences on the associations observed in our study.

The QRFPR gene (also known as GPR103 gene) encodes a G-protein-coupled receptor that is expressed in brain, urinary bladder, thyroid gland, human pancreatic β-cells, and so on.²¹ Recent data indicate that QRFPR and its ligand QRFP are key regulators of glucose homeostasis²² by increasing insulin sensitivity and insulin secretion,^23,24 as well as preventing pancreatic β-cells death and apoptosis.²⁵ Besides, a locus on chromosome 4q27 that is approximately 625kb downstream the QRFPR gene was reported to be associated with a number of autoimmune diseases.^26–29 Therefore, it is possible that QRFPR or its signaling pathway may be involved in the development of T1D plus AITD.

According to bioinformatics analysis, rs7679475 is located in the sequence upstream of QRFPR. The variant is enriched for enhancer/promoter-related histone marker H3K27ac. It also can change the binding of regulatory motifs (rachyury_2, PPAR_1, and SIX5_known1). All above results provided evidence that the polymorphism might fall within putative functional regulatory regions. Therefore, we hypothesized that the rs7679475G allele is likely to change QRFPR expression in thyroid and pancreatic β-cells, influencing immunity and β-cells apoptosis, thus, predisposing AITD and T1D. Future studies are required to confirm the function of the rs7679475 on QRFPR expression since we cannot fully exclude the possibility that this SNP is not itself responsible for the observed association, but instead is in LD with a causative SNP.

Female has been reported to be strongly related to a large number of autoimmune diseases.¹¹ Our results demonstrated a clear difference in the magnitude of AITD risk in T1D across rs769475 genotypes according to sex, and revealed a possible gene–sex interaction. But the mechanisms behind this are still unclear. A previous study in female rats supported a role for estradiol in the regulation of hypothalamic prepro-QRFP mRNA expression.³⁰ Nevertheless, to the best of our knowledge, there is no research to date investigating interactions between sex hormones and QRFP or QRFPR expression in peripheral organs.

Several limitations should be considered when interpreting our findings. First, our study is hospital-based design, which might induce selection bias. Second, previous studies have reported that the GADA may influence the development or progression of AITD,^31,32 which were consistent with our results. However, adjustment of GADA did not appreciably change the result for rs7679475 and AITD. Still, we lacked other T1D antibodies for every participant; therefore, we could not determine potential impacts from the existence of other T1D antibodies. Finally, this study was conducted in Chinese Han patients, so our findings may not be generalizable to other ancestry groups. The replication of these results in multi-ethnic larger samples and the research of potential molecular mechanisms underlying the genetic associations are warranted.

Conclusion

In summary, our results suggested that in Chinese Han patients with T1D, QRFPR rs7679475 polymorphism may be associated with risk for AITD and this effect may be modulated by sex. Further work is required to understand the particular mechanisms underlying the observed effects.