Genome-Wide Identification of N6-Methyladenosine Associated SNPs as Potential Functional Variants for Type 1 Diabetes

Objectives: N6-methyladenosine (m6A) is essential in the regulation of the immune system, but the role that its single nucleotide polymorphisms (SNPs) play in the pathogenesis of type 1 diabetes (T1D) remains unknown. This study demonstrated the association between genetic variants in m6A regulators and T1D risk based on a case-control study in a Chinese population.

Methods: The tagging SNPs in m6A regulators were genotyped in 1005 autoantibody-positive patients with T1D and 1257 controls using the Illumina Human OmniZhongHua-8 platform. Islet-specific autoantibodies were examined by radioimmunoprecipitation in all the patients. The mixed-meal glucose tolerance test was performed on 355 newly diagnosed patients to evaluate their residual islet function. The functional annotations for the identified SNPs were performed in silico. Using 102 samples from a whole-genome expression microarray, key signaling pathways associated with m6A regulators in T1D were comprehendingly evaluated.

Results: Under the additive model, we observed three tag SNPs in the noncoding region of the PRRC2A (rs2260051, rs3130623) and YTHDC2 (rs1862315) gene are associated with T1D risk. Although no association was found between these SNPs and islet function, patients carrying risk variants had a higher positive rate for ZnT8A, GADA, and IA-2A. Further analyses showed that rs2260051[T] was associated with increased expression of PRRC2A mRNA (P = 7.0E-13), and PRRC2A mRNA was significantly higher in peripheral blood mononuclear cell samples from patients with T1D compared to normal samples (P = 0.022). Enrichment analyses indicated that increased PRRC2A expression engages in the most significant hallmarks of cytokine-cytokine receptor interaction, cell adhesion and chemotaxis, and neurotransmitter regulation pathways. The potential role of increased PRRC2A in disrupting immune homeostasis is through the PI3K/AKT pathway and neuro-immune interactions.

Conclusion: This study found intronic variants in PRRC2A and YTHDC2 associated with T1D risk in a Chinese Han population. PRRC2A rs2260051[T] may be implicated in unbalanced immune homeostasis by affecting the expression of PRRC2A mRNA. These findings enriched our understanding of m6A regulators and their intronic SNPs that underlie the pathogenesis of T1D.

Introduction

Type 1 diabetes (T1D) is an autoimmune disease involving insulin-producing pancreatic islet beta cell prolonged destruction by autoreactive effector T lymphocytes (1). Although management with exogenous insulin changes T1D from an acutely fatal disease into a chronic disease, achieving glycemic control is still subject to significant limitations (2). The incidence of T1D continuously increases worldwide while varies a lot in different regions. In addition to environmental triggers, genetic factors likely play a vital role in the development of T1D. To date, nearly 60 risk loci have been identified in genome-wide association studies (GWAS), including MHC (major histocompatibility complex), INS, PTPN22 (protein tyrosine phosphatase N2), and so forth (3, 4). Most risk variants for T1D are located in noncoding regulatory regions (5) and are supposed to take part in altering gene expression. However, few studies have focused on potential functional variants that may affect N6-methyladenosine (m6A) methylation; the cellular mechanisms corresponding to their pathogenic processes are poorly understood.

N6-methyladenosine (m6A) is the most common internal modification of mRNA in eukaryotes. As a critical chemical mark, m6A regulates post-transcriptional gene expression by affecting RNA splicing, nuclear export, stability, and translation (6). The m6A RNA modification is reversible and can be dynamically modulated by a series of regulators. Among them, methyltransferase “writers”, m6A-binding proteins “readers”, and demethylase “erasers” are responsible for catalyzing, recognizing, and removing the m6A marks, respectively (7). Emerging studies have demonstrated that aberrant expression of m6A regulators such as METTL3, METTL14, and FTO can lead to dysregulated m6A modification, which is involved in the pathophysiology of various diseases, including cancers, immune disorders, viral infection, etc. (8–10). Given its importance in regulating gene expression and immune homeostasis, m6A modification could be involved in the pathogenesis of T1D.

Several single nucleotide polymorphisms (SNPs) located in m6A regulator genes have been identified as genetic variants related to different diseases, such as obesity (11), cancers (12), and autoimmune diseases (13). However, the association between genetic variants in m6A regulator genes and T1D risk remains elusive. To further investigate the potential pathogenesis of m6A modification in T1D, it is necessary to evaluate the functional variants in m6A modification.

Therefore, we conducted a case-control study of 1005 patients with T1D and 1257 healthy controls in a Chinese population to evaluate the association between tagging SNPs in m6A regulators and T1D risk. Next, we investigated the association between risk loci and clinical features. To verify whether these m6A regulator SNPs have potential regulatory effects on gene expression and whether these disturbances may lead to T1D, we performed eQTL and differential expression analyses using the GTEx portal and the GEO public dataset.

Materials and Methods

Study Subjects

A total of 1005 patients with T1D and 1257 healthy controls were recruited from the First Affiliated Hospital of Nanjing Medical University between January 2008 and December 2016. All patients with T1D had been diagnosed according to the American Diabetes Association criteria and had a positive test for at least one pancreatic islet-specific autoantibody (ZnT8A, GADA, IA-2A, IAA) (14). Those with the clinical features of latent autoimmune diabetes in adults were excluded. The inclusion criteria for the non-diabetic healthy controls have been described previously (15). Briefly, the control subjects (422 female, 835 male; mean age 41.96 ± 16.36 years) were enrolled from healthy individuals in Jiangsu province without the medical history of diabetes or evident autoimmune diseases. Written informed consent was obtained from all participants or their guardians. Each participant was interviewed to collect demographic information and peripheral blood samples. All subjects were unrelated members of a Chinese Han population. Their clinical characteristics are shown in Supplementary Table 1. This study was approved by the Ethics Committee of the First Affiliated Hospital of Nanjing Medical University.

Selection of m6A Regulator Tag SNPs and Genotyping

The m6A RNA methylation regulators consist of 10 writers (METTL3/14/16, WTAP, VIRMA, ZC3H13, CBLL1, RBM15/15B, ZCCHC4), two erasers (FTO, ALKBH5), and 16 readers (YTHDC1/2, YTHDF1/2/3, HNRNPC/G, HNRNPA2B1, IGF2BP1/2/3, FMR1, PRRC2A, LRPPRC, ELAVL1, ZFP217) (6, 9, 16). A total of 21,912 SNPs in the m6A regulator gene and within ±10 kb flanking regions were retrieved from the 1000 Genomes project in the Chinese Han population (CHB) data (https://www.internationalgenome.org/). Genomic DNA extracted from peripheral blood lymphocytes was genotyped using the Illumina Human OmniZhongHua-8 platform, providing 787,224 qualified SNPs. Then we extracted m6A regulator variants from it based on quality control and imputation protocols (3). First, SNPs with minor allele frequency < 0.05 or a Hardy-Weinberg equilibrium (HWE) < 0.05 were excluded. Next, Haploview 4.2 software was used to select tag SNPs with an r² threshold > 0.4. Finally, 39 SNPs with call rates > 95% were selected. Detailed information on these SNPs is presented in Supplementary Table 2. In addition, variants were excluded from subsequent association analyses with differences in effect size compared to the European ancestry meta-analyzed UCSD T1D-GWAS (5).

Pancreatic Islet-Specific Antibodies Detection

Autoantibodies to ZnT8, GAD65, and IA-2 were detected separately using protein A radio-binding assays (1). Serum IAA was only tested within two weeks after diagnosis using an ELISA kit (Biomerica, USA) to avoid interference of injected insulin.

Mixed-Meal Tolerance Test

A total of 355 eligible newly diagnosed T1D patients underwent a mixed-meal tolerance test (MMTT) to assess their residual β cell function. Plasma C-peptide concentrations at baseline, 60 min, and 120 min after an MMTT were measured using an automated chemiluminescence assay (Roche, Switzerland). The 2-hour C-peptide area under the curve (AUC) was calculated with the trapezoidal method.

Functional Annotation for m6A Regulator Loci

The GTEx v8 Portal (http://www.gtexportal.org/home/) was searched to conduct the expression quantitative trait loci (eQTL) analyses. Then the functional annotations for the identified SNPs were performed with HaploReg v4.1 (https://pubs.broadinstitute.org/mammals/haploreg/haploreg.php). Expression profiles were downloaded from the GEO database [GSE9006 (17) and GSE72492 (18)]. We compared the expression data of PRRC2A and YTHDC2 in peripheral blood mononuclear cell (PBMC) samples between 79 T1D and 23 healthy controls. Additionally, the YTHDC2 mRNA expression in pancreatic samples was compared between 10 T1D and 7 controls. Then gene expression data of PBMC samples from patients with T1D were divided into high- and low- groups according to the median PRRC2A expression level. The differentially expressed genes (DEGs) between high and low PRRC2A expression groups were identified using the limma package (19), and the thresholds were |log2FoldChange| > 0.5 and P value < 0.05. GO and KEGG pathway analysis was performed on the DEGs between high and low PRRC2A expression groups with the clusterProfiler R package (20). A false discovery rate (FDR) cutoff of 0.1 was used for statistical significance. To elucidate the differences between high- and low- PRRC2A groups, gene set enrichment analysis (GSEA) was performed using the molecular signature database (MSigDB) collections: c2.all.v7.2.symbols, c7.all.v7.2.symbols by R package clusterProfiler (20, 21). Normalized enrichment scores (|NES| > 1), adjusted P-values < 0.05 and FDR q value < 0.25 were considered statistically significant.

Statistical Analysis

Departures from the HWE in control subjects were verified by the goodness-of-fit χ2 test. OR (95% CI) and P values were derived from logistic regression analyses with adjustment for sex and the ten principal components under the assumption of an additive genetic model. The principal components are described in a previous report (3). Multivariate logistic regression analyses for islet function were also performed to assess its effects on T1D susceptibility. The P values of eQTL were obtained from the linear regression model. Differences in the distribution of categorical variables were evaluated with the χ2 test, and the continuous variables were compared with Welch’s t-test. In cases of more than two groups, one-way ANOVA or the Kruskal-Wallis test was used. All analyses were performed using R 4.0.5 and plink1.07. Two-tailed P values below 0.05 were considered statistically significant.

Results

SNPs in PRRC2A and YTHDC2 Genes Are Associated With T1D Risk

Among all the 39 selected SNPs, three tag SNPs, namely, rs2260051 (A>T), rs3130623 (C>T), and rs1862315 (A>C) were significantly associated with T1D risk under the addictive model (P = 7.87E-05 for rs2260051, P = 3.72E-13 for rs3130623, and P = 0.045 for rs1862315). The results were further validated in an independent UCSD T1D-GWAS (5) dataset containing 73,011 individuals. Among these, rs2260051 and rs3130623 reached the threshold of GWAS significance (Table 1). In addition, we calculated the cumulative risk score across the three SNPs. As shown in Table 2, individuals who carry 3-4 risk alleles are associated with higher incidences of T1D (OR = 1.83, 95% confidence interval (CI): 1.52-2.21, P = 1.98E-10). A stronger trend was found for individuals with 5 or 6 risk alleles (OR = 4.92, 95% CI: 2.07-13.03, P = 5.65E-04). It suggests a significant association between cumulative risk score in PRRC2A and YTHDC2 gene with T1D risk in a Chinese Han population (P_trend = 1.93E-12).

TABLE 1

Table 1 Associations between the m6A regulator SNPs in PRRC2A, YTHDC2 gene, and T1D risk.

TABLE 2

Table 2 Associations between the m6A regulator SNPs cumulative risk score and T1D risk in a Chinese population.

The Significant Association Between m6A Regulator SNPs and Islet Autoimmunity

To explore the potential effects induced by m6A regulator SNPs, we also evaluated the association between the genotype distribution of the three SNPs and the positivity of islet-specific autoantibodies in 1005 T1D patients. As shown in Table 3, after adjusting for sex and the ten principal components, ZnT8A, GADA, and IA-2A were all positively associated with rs2260051 and rs3130623. Patients with T1D who carried a high cumulative risk score had a higher frequency of islet autoantibody positivity. However, there was no significant association between YTHDC2 rs1862315 and islet autoantibody positivity (all P > 0.05).

TABLE 3

Table 3 Associations of the m6A regulator SNPs with islet autoantibody positivity.

No Significant Correlation Between m6A Regulator SNPs and Residual C-Peptide

To assess the association between m6A regulator SNPs and residual islet function, we conducted an MMTT on 355 newly diagnosed patients with T1D. Generalized linear model analyses were performed to determine statistical significance and beta coefficients for age at onset, disease duration, body mass index, positivity of islet-specific autoantibodies, and residual islet β-cell function (fasting C-peptide, random C-peptide, or area under the C-peptide curve). As is shown in Table 4, there were no significant associations between m6A regulator SNPs and residual islet β-cell function (P > 0.05).

TABLE 4

Table 4 Associations of the m6A regulator SNPs with residual islet function in newly diagnosed patients with T1D.

Rs2260051 Was Associated With Increased Expression of PRRC2A

Functional explorations using HaploReg annotated these SNPs to evaluate how they confer disease susceptibility. Rs1862315 acted as a promoter for H3K4me3, H3K9ac, and H3K27ac in different immune cell subsets, including T cell, B cell, monocyte, etc. The other two SNPs, rs2260051 and rs3130623, are also located in the enhancer histone marks for H3K4me1 and H3K27ac. (Supplementary Table 3) We also found that the risk allele of rs2260051 was associated with increased expression of PRRC2A in the whole blood sample (P = 7.0E-13) (Figure 1A), while rs1862315 was associated with increased expression of YTHDC2 both in whole blood (P = 2.8E-36) (Figure 1B) and pancreas (P = 3.5E-9) (Figure 1C). Moreover, using transcriptome dataset from the GEO, we found that the mRNA expression of PRRC2A was increased in PBMC samples from patients with T1D relative to the controls (P = 0.022) (Figure 1D). However, there was no significant difference in YTHDC2 expression between T1D and healthy, neither from the PBMC samples (Figure 1E) nor the pancreatic samples (Figure 1F). These results provide evidence that the three m6A regulator SNPs are in the functional promoter or enhancer region of the PRRC2A and YTHDC2 genes. Furthermore, rs2260051 may induce T1D by regulating the expression of PRRC2A.

FIGURE 1

Figure 1 Effects of variants on gene expression of PRRC2A and YTHDC2. (A) The effect allele of rs2260051 was significantly associated with increased expression of PRRC2A in whole blood samples. The minor allele of rs1862315 was significantly associated with increased expression of YTHDC2 in (B) whole blood and (C) pancreas. (D) The expression levels of PRRC2A mRNA were significantly higher in the peripheral blood mononuclear cell (PBMC) samples from type 1 diabetes (T1D) than the healthy control. No significant differences in YTHDC2 between T1D and healthy control from (E) PBMC or (F) pancreatic samples.

Increased Expression of PRRC2A Associated With Disturbed Immune Homeostasis

A total of 800 DEGs (668 upregulated and 132 downregulated, Supplementary Table 4) were identified between the high and low PRRC2A expression groups (Figure 2A). To predict the functions of the interactive genes of PRRC2A, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were also performed. We found that cytokine-cytokine receptor interaction involved in immune response, cell chemotaxis pathways involved in cell adhesion, and trans-synaptic signaling involved in neurotransmitter regulation were significantly enriched by PRRC2A-associated DEGs (Figure 2B). Moreover, gene set enrichment analysis (GSEA) between high and low PRRC2A expression groups revealed significant differences (FDR < 0.25, P < 0.05, and |NES| >1) in the enrichment of MsigDB collection. Figure 2C and Supplementary Table 5 show that calcium signaling, PI3K/AKT signaling, cytokine-cytokine receptor interaction, and inflammatory response pathways were significantly enriched in T1D with high PRRC2A expression. In addition, neuroactive ligand-receptor interaction, and neurotransmitter release cycle pathways (acetylcholine, norepinephrine, serotonin) were significantly enriched in T1D with high PRRC2A expression (Figure 2D and Supplementary Table 5). Downregulation of Treg and Th2 were also enriched in the PRRC2A high-expression group. In contrast, upregulation of effector memory CD4+ T, CD8+ T cells, NK cells, and plasma cells were differentially enriched in PRRC2A high-expression phenotypes (Figure 2E and Supplementary Table 5).

FIGURE 2

Figure 2 The identification and functional enrichment of differentially expressed genes based on PRRC2A mRNA expression in type 1 diabetes. (A) Volcano plot of the differentially expressed genes (DEGs) between groups exhibiting high and low levels of PRRC2A mRNA. (B) Significantly enriched Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway terms of PRRC2A associated DEGs. (C) Gene set enrichment analysis (GSEA) indicated that high expression level of PRRC2A was correlated with calcium signaling, PI3K/AKT signaling, cytokine-cytokine receptor interaction, and inflammatory response pathways. (D) GSEA indicated that high expression level of PRRC2A was correlated with neuroactive ligand-receptor interaction, and neurotransmitter (acetylcholine, norepinephrine, serotonin) release cycle pathways. (E) GSEA indicated that high expression level of PRRC2A was correlated with downregulation of Treg and Th2, and upregulation of effector memory CD4+ T, CD8+ T cells, NK cells, and plasma cells.

Discussion

Previous studies have established the important role of m6A modification in modulating gene expression that regulates T cell differentiation, proliferation, and cytokine production (22, 23). Taking the importance of m6A regulators in shaping balanced immune homeostasis into account, allelic variants in m6A regulators could be genetic risk factors for autoimmune disease. Polymorphisms in the ALKBH5 (24) and METTL21B (25) genes have been identified as genetic variants associated with autoimmune thyroid disease and multiple sclerosis. However, little is found in the literature on m6A regulator gene susceptibility in T1D. In this study, we genotyped 21,912 SNPs in 28 m6A regulator genes and evaluated their association with T1D risk in a case-control study. We found that three SNPs within the intronic region of PRRC2A and YTHDC2 were significantly associated with a higher risk for T1D.

The proline-rich coiled-coil 2A (PRRC2A) gene, also known as BAT2 (26), is located in the human major histocompatibility complex (MHC) class III region. It has been newly recognized as an m6A reader in controlling neural development via regulating the stability of the Olig2 mRNA (27). Several reports have shown that microsatellite repeats and missense polymorphisms in the PRRC2A gene conferred increased risk in various immune-related diseases, including type 1 diabetes, rheumatoid arthritis, coeliac disease, lupus nephritis (28–31). Masao et al. first reported that BAT2 microsatellite alleles are associated with the age-at-onset of insulin-dependent diabetes mellitus (28). The type 1 Diabetes Genetics Consortium (T1DGC) study also detected seven missense SNPs in the BAT2 gene associated with type 1 diabetes (32). This study further found two independent intronic variants in the PRRC2A gene, rs2260051, and rs3130623, associated with T1D risk. Both loci were identified in the Chinese population and validated in an independent European ancestry UCSD T1D-GWAS dataset (5). Over 95% of the variants in high LD were in the intronic region, making it difficult to elucidate their functional roles in the pathogenesis of human disease (33). Based on eQTL analysis, we found that the rs2260051 may influence the mRNA expression level of PRRC2A. Interestingly, rs3130623 is in moderate linkage disequilibrium with another PRRC2A variant, rs2736157, which was recently reported to be associated with neuromyelitis optica and multiple sclerosis in a Chinses cohort (34). Thus, these PRRC2A variants may provide insight into the pathogenesis of dysregulated m6A modification in autoimmune diseases.

As an N6-methyladenosine reader, YTHDC2 can selectively bind RNA with m6A modification. It plays a significant role in enhancing the translation efficiency of its targets and decreasing their mRNA abundance (35). Several studies have reported the critical functions of YTHDC2 during spermatogenesis and oncogenesis (12, 36). Daniele et al. also reported YTHDC2 gene could be involved in pancreatic cancer susceptibility (37). In this study, we identified rs1862315, in the promotor of YTHDC2, was significantly associated with the susceptibility of T1D. Further functional annotation in silico demonstrated that rs1862315 upregulates YTHDC2 in the whole blood and pancreas sample. However, no significant difference was found in YTHDC2 expression between T1D and healthy controls.

We also showed that a genetic score based on the cumulative effects of variants is associated with increased risk for T1D. It has been demonstrated that T1D is a polygenic disease, and the combination of multiple SNPs may reveal a greater risk than single SNP (38). In addition, multiple variants of a single gene locus can have an additive effect on affecting gene expression (39). We and others previously found that a higher genetic risk score (GRS) was associated with lower fasting C-peptide levels at diagnosis (3). It also resulted in a sharper decline in residual β-cell function following the T1D diagnosis (40). However, this study did not observe any association between risk alleles in PRRC2A or YTHDC2 and lower residual islet function. This might indicate that other SNPs involved in the GRS model but not PRRC2A or YTHDC2 risk alleles affect β-cell function.

Islet autoantibodies are important markers of islet autoimmunity, indicating disease progression and diagnosis. However, there is considerable interindividual variability in the progression rate of islet beta cell destruction (41). Genetic susceptibility is likely to play a significant role in it. Several studies have identified correlations between immune-related variants and the positive rate of islet autoantibodies, such as IL27-IA2A, IFIH1-IA2A, IL2RA-GADA, IL1B-ZnT8A, and CTLA4-IA2A (15, 42, 43). Our results confirm an association between PRRC2A risk alleles and islet autoantibody positivity. It indicates that PRRC2A risk variants might regulate autoimmunity and have the potential to be biomarkers of T1D progression.

In our findings, the risk variant of rs2260051 was significantly associated with enhanced expression of PRRC2A. We conduct bioinformatic analyses to investigate the function of PRRC2A in T1D using PBMC samples from new-onset patients with T1D. On the one hand, KEGG and GSEA analysis found cytokine-cytokine receptor interaction and PI3K/AKT signaling pathway were significantly enriched in T1D with high PRRC2A expression. Prior studies have noted the importance of the PI3K/AKT signaling pathway in regulating immunity. During the adaptive immune response, TCR or BCR stimulation activates the PI3K/AKT pathway, which takes part in the regulation of lymphocyte differentiation in turn (44). The GSEA results also suggested that higher expression of PRRC2A may enhance the number or function of pathogenic immune cells (effector memory T cells, NK cells, and plasma cells) while limiting the regulatory T cells. These findings suggested that increased PRRC2A in T1D was positively correlated with immune response and lymphocyte differentiation. On the other hand, neurotransmitter regulation was enriched by PRRC2A-associated DEGs. Recently studies revealed the role of neurotransmitters in the modulation of the immune system. Immune cells can synthesize and release neurotransmitters such as acetylcholine and catecholamines. The receptors of these neuromodulators have also been identified on monocytes, lymphocytes, etc. In addition, cytokine receptors are also expressed on sensory neurons (45). These findings are in accord with the present study, which indicates that increased expression of PRRC2A may be involved in imbalanced immune homeostasis in T1D through the PI3K/AKT pathway and neuro-immune interactions.

Our study was the first to report an association between the variants in m6A regulator genes and T1D risk. These results should be interpreted with caution. Although rs2260051 and rs3130623 in the PRRC2A gene reached the threshold of GWAS significance in an independent validation, the significance of YTHDC2 rs1862315 was uncorrected for multiple comparisons. Although we elaborated the functional annotations of the m6A regulator SNPs based on the GTEx Portal and RNA expression profiling from the GEO database, they could not be directly validated due to the absence of experiments. Further studies are needed to clarify the biological mechanisms.

Overall, this study found that m6A regulator SNPs in the functional elements of PRRC2A and YTHDC2 gene, including rs2260051, rs3130623, and rs1862315, were associated with increased risk for T1D. Of these variants, rs2260051 and rs3130623 showed significant association with islet autoantibody positivity in T1D individuals. For rs2260051, the risk allele was associated with increased expression of PRRC2A. Transcriptomic analyses confirmed that higher PRRC2A expression correlates with disturbed immune homeostasis. These findings might inspire further research to investigate the roles of PRRC2A and m6A regulator SNPs in the development of T1D.

Data Availability Statement

The genotype data presented in the study are deposited in the GWAS Catalog repository, accession number GCST008377, GCST90014023. The RNA microarray data presented in the study are deposited in the GEO repository, accession number GSE9006, GSE72492. Further inquiries can be directed to the corresponding author.

Ethics Statement

The studies involving human participants were reviewed and approved by Ethics Committee of the First Affiliated Hospital of Nanjing Medical University. Written informed consent to participate in this study was provided by the participants’ legal guardian/next of kin.

Author Contributions

YC performed data analysis and drafted the initial manuscript. MS validated the results and helped with the data visualization. CJ and YH performed data analysis. YS, LJ, and YQ helped prepare for data collection and the literature search. YG and QF helped with the manuscript review and editing. HC contributed to the laboratory testing. KX worked on the methodology and supervision. TY designed the study and is the guarantor of this work. As such, TY has full access to all the study data and takes responsibility for the integrity of the data and the accuracy of the data analysis. All authors contributed to the article and approved the submitted version.

Funding

This study was supported by the National Natural Science Foundation of China [Grant number: 81900708, 81830023, 81530026, 82000747], the National Key Research and Development Program of China [Grant number:2016YFC1305004], and the Postgraduate Research & Practice Innovation Program of Jiangsu Province [Grant number: JX10213850]. The sponsors had no role in the study design; the collection, analysis, and interpretation of data; the writing of the report; or the decision to submit the article for publication.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We thank Professor Perrin C. White (University of Texas Southwestern Medical Center), Professor Kyle J. Gaulton, Maike Sander (University of California San Diego), Professor Linda Yip (Stanford University), and their colleagues for kindly sharing their data. We thank Dr. Jianming Zeng (University of Macau) and all the members of his bioinformatics team, trainees, for generously sharing their experience and codes.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fendo.2022.913345/full#supplementary-material

References

1. Liu J, Bian L, Ji L, Chen Y, Chen H, Gu Y, et al. The Heterogeneity of Islet Autoantibodies and the Progression of Islet Failure in Type 1 Diabetic Patients. Sci China Life Sci (2016) 59:930–9. doi: 10.1007/s11427-016-5052-3

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Warshauer JT, Bluestone JA, Anderson MS. New Frontiers in the Treatment of Type 1 Diabetes. Cell Metab (2020) 31:46–61. doi: 10.1016/j.cmet.2019.11.017

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Zhu M, Xu K, Chen Y, Gu Y, Zhang M, Luo F, et al. Identification of Novel T1D Risk Loci and Their Association With Age and Islet Function at Diagnosis in Autoantibody-Positive T1D Individuals: Based on a Two-Stage Genome-Wide Association Study. Diabetes Care (2019) 42:1414–21. doi: 10.2337/dc18-2023

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Pociot F, Lernmark Å. Genetic Risk Factors for Type 1 Diabetes. Lancet (2016) 387:2331–9. doi: 10.1016/S0140-6736(16)30582-7

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Chiou J, Geusz RJ, Okino M-L, Han JY, Miller M, Melton R, et al. Interpreting Type 1 Diabetes Risk With Genetics and Single-Cell Epigenomics. Nature (2021) 594:398–402. doi: 10.1038/s41586-021-03552-w

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Zaccara S, Ries RJ, Jaffrey SR. Reading, Writing and Erasing mRNA Methylation. Nat Rev Mol Cell Biol (2019) 20:608–24. doi: 10.1038/s41580-019-0168-5

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Meyer KD, Jaffrey SR. Rethinking M6a Readers, Writers, and Erasers. Annu Rev Cell Dev Biol (2017) 33:319–42. doi: 10.1146/annurev-cellbio-100616-060758

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Deng X, Su R, Weng H, Huang H, Li Z, Chen J. RNA N-Methyladenosine Modification in Cancers: Current Status and Perspectives. Cell Res (2018) 28:507–17. doi: 10.1038/s41422-018-0034-6

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Shi H, Wei J, He C. Where, When, and How: Context-Dependent Functions of RNA Methylation Writers, Readers, and Erasers. Mol Cell (2019) 74:640–50. doi: 10.1016/j.molcel.2019.04.025

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Tong J, Cao G, Zhang T, Sefik E, Amezcua Vesely MC, Broughton JP, et al. M6a mRNA Methylation Sustains Treg Suppressive Functions. Cell Res (2018) 28:253–6. doi: 10.1038/cr.2018.7

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Dina C, Meyre D, Gallina S, Durand E, Körner A, Jacobson P, et al. Variation in FTO Contributes to Childhood Obesity and Severe Adult Obesity. Nat Genet (2007) 39:724–6. doi: 10.1038/ng2048

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Liu J, Wang D, Zhou J, Wang L, Zhang N, Zhou L, et al. N6-Methyladenosine Reader YTHDC2 and Eraser FTO may Determine Hepatocellular Carcinoma Prognoses After Transarterial Chemoembolization. Arch Toxicol (2021) 95:1621–9. doi: 10.1007/s00204-021-03021-3

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Wang Y, Li L, Li J, Zhao B, Huang G, Li X, et al. The Emerging Role of M6a Modification in Regulating the Immune System and Autoimmune Diseases. Front Cell Dev Biol (2021) 9:755691. doi: 10.3389/fcell.2021.755691

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Chiang JL, Kirkman MS, Laffel LM, Peters AL, A. Type 1 Diabetes Sourcebook. Type 1 Diabetes Through the Life Span: A Position Statement of the American Diabetes Association. Diabetes Care (2014) 37:2034–54. doi: 10.2337/dc14-1140

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Chen Y, Chen S, Gu Y, Feng Y, Shi Y, Fu Q, et al. CTLA-4 +49 G/A, a Functional T1D Risk SNP, Affects CTLA-4 Level in Treg Subsets and IA-2A Positivity, But Not Beta-Cell Function. Sci Rep (2018) 8:10074. doi: 10.1038/s41598-018-28423-9

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Lee D-F, Walsh MJ, Aguiló F. ZNF217/ZFP217 Meets Chromatin and RNA. Trends Biochem Sci (2016) 41:986–8. doi: 10.1016/j.tibs.2016.07.013

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Kaizer EC, Glaser CL, Chaussabel D, Banchereau J, Pascual V, White PC. Gene Expression in Peripheral Blood Mononuclear Cells From Children With Diabetes. J Clin Endocrinol Metab (2007) 92:3705–11. doi: 10.1210/jc.2007-0979

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Yip L, Fuhlbrigge R, Alkhataybeh R, Fathman CG. Gene Expression Analysis of the Pre-Diabetic Pancreas to Identify Pathogenic Mechanisms and Biomarkers of Type 1 Diabetes. Front Endocrinol (Lausanne) (2020) 11:609271. doi: 10.3389/fendo.2020.609271

PubMed Abstract | CrossRef Full Text | Google Scholar

19. Ritchie ME, Phipson B, Wu D, Hu Y, Law CW, Shi W, et al. Limma Powers Differential Expression Analyses for RNA-Sequencing and Microarray Studies. Nucleic Acids Res (2015) 43:e47. doi: 10.1093/nar/gkv007

PubMed Abstract | CrossRef Full Text | Google Scholar

20. Yu G, Wang LG, Han Y, He QY. Clusterprofiler: An R Package for Comparing Biological Themes Among Gene Clusters. Omics J Integr Biol (2012) 16:284–7. doi: 10.1089/omi.2011.0118

CrossRef Full Text | Google Scholar

21. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene Set Enrichment Analysis: A Knowledge-Based Approach for Interpreting Genome-Wide Expression Profiles. Proc Natl Acad Sci USA (2005) 102:15545–50. doi: 10.1073/pnas.0506580102

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Li H-B, Tong J, Zhu S, Batista PJ, Duffy EE, Zhao J, et al. mA mRNA Methylation Controls T Cell Homeostasis by Targeting the IL-7/STAT5/SOCS Pathways. Nature (2017) 548:338–42. doi: 10.1038/nature23450

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Shulman Z, Stern-Ginossar N. The RNA Modification N6-Methyladenosine as a Novel Regulator of the Immune System. Nat Immunol (2020) 21:501–12. doi: 10.1038/s41590-020-0650-4

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Song R-H, Zhao J, Gao C-Q, Qin Q, Zhang J-A. Inclusion of ALKBH5 as a Candidate Gene for the Susceptibility of Autoimmune Thyroid Disease. Adv Med Sci (2021) 66:351–8. doi: 10.1016/j.advms.2021.07.006

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Mo X-B, Lei S-F, Qian Q-Y, Guo Y-F, Zhang Y-H, Zhang H. Integrative Analysis Revealed Potential Causal Genetic and Epigenetic Factors for Multiple Sclerosis. J Neurol (2019) 266:2699–709. doi: 10.1007/s00415-019-09476-w

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Jin G, Zhu M, Yin R, Shen W, Liu J, Sun J, et al. Low-Frequency Coding Variants at 6p21.33 and 20q11.21 are Associated With Lung Cancer Risk in Chinese Populations. Am J Hum Genet (2015) 96:832–40. doi: 10.1016/j.ajhg.2015.03.009

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Wu R, Li A, Sun B, Sun JG, Zhang J, Zhang T, et al. A Novel M(6)A Reader Prrc2a Controls Oligodendroglial Specification and Myelination. Cell Res (2019) 29:23–41. doi: 10.1038/s41422-018-0113-8

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Hashimoto M, Nakamura N, Obayashi H, Kimura F, Moriwaki A, Hasegawa G, et al. Genetic Contribution of the BAT2 Gene Microsatellite Polymorphism to the Age-at-Onset of Insulin-Dependent Diabetes Mellitus. Hum Genet (1999) 105:197–9. doi: 10.1007/s004399900100

PubMed Abstract | CrossRef Full Text | Google Scholar

29. Goudey B, Abraham G, Kikianty E, Wang Q, Rawlinson D, Shi F, et al. Interactions Within the MHC Contribute to the Genetic Architecture of Celiac Disease. PLoS One (2017) 12:e0172826. doi: 10.1371/journal.pone.0172826

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Xu R, Li Q, Liu R, Shen J, Li M, Zhao M, et al. Association Analysis of the MHC in Lupus Nephritis. J Am Soc Nephrol JASN (2017) 28:3383–94. doi: 10.1681/ASN.2016121331

CrossRef Full Text | Google Scholar

31. Singal DP, Li J, Lei K. Genetics of Rheumatoid Arthritis (RA): Two Separate Regions in the Major Histocompatibility Complex Contribute to Susceptibility to RA. Immunol Lett (1999) 69:301–6. doi: 10.1016/S0165-2478(99)00108-X

PubMed Abstract | CrossRef Full Text | Google Scholar

32. He C, Hamon S, Li D, Barral-Rodriguez S, Ott J, Diabetes Genetics C. MHC Fine Mapping of Human Type 1 Diabetes Using the T1DGC Data. Diabetes Obes Metab (2009) 11(Suppl 1):53–9. doi: 10.1111/j.1463-1326.2008.01003.x

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Ward LD, Kellis M. Interpreting Noncoding Genetic Variation in Complex Traits and Human Disease. Nat Biotechnol (2012) 30:1095–106. doi: 10.1038/nbt.2422

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Zhang J, Chen MJ, Zhao GX, Li HF, Wu L, Xu YF, et al. Common Genetic Variants in PRRC2A Are Associated With Both Neuromyelitis Optica Spectrum Disorder and Multiple Sclerosis in Han Chinese Population. J Neurol (2021) 268:506–15. doi: 10.1007/s00415-020-10184-z

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Yang N, Ying P, Tian J, Wang X, Mei S, Zou D, et al. Genetic Variants in M6a Modification Genes are Associated With Esophageal Squamous-Cell Carcinoma in the Chinese Population. Carcinogenesis (2020) 41:761–8. doi: 10.1093/carcin/bgaa012

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Hsu PJ, Zhu Y, Ma H, Guo Y, Shi X, Liu Y, et al. Ythdc2 Is an N-Methyladenosine Binding Protein That Regulates Mammalian Spermatogenesis. Cell Res (2017) 27:1115–27. doi: 10.1038/cr.2017.99

PubMed Abstract | CrossRef Full Text | Google Scholar

37. Fanale D, Iovanna JL, Calvo EL, Berthezene P, Belleau P, Dagorn JC, et al. Germline Copy Number Variation in the YTHDC2 Gene: Does it Have a Role in Finding a Novel Potential Molecular Target Involved in Pancreatic Adenocarcinoma Susceptibility? Expert Opin Ther Targets (2014) 18:841–50. doi: 10.1517/14728222.2014.920324

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Ilonen J, Lempainen J, Veijola R. The Heterogeneous Pathogenesis of Type 1 Diabetes Mellitus. Nat Rev Endocrinol (2019) 15:635–50. doi: 10.1038/s41574-019-0254-y

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Broekema RV, Bakker OB, Jonkers IH. A Practical View of Fine-Mapping and Gene Prioritization in the Post-Genome-Wide Association Era. Open Biol (2020) 10:190221. doi: 10.1098/rsob.190221

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Williams MD, Bacher R, Perry DJ, Grace CR, McGrail KM, Posgai AL, et al. Genetic Composition and Autoantibody Titers Model the Probability of Detecting C-Peptide Following Type 1 Diabetes Diagnosis. Diabetes (2021) 70:932–43. doi: 10.2337/db20-0937

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Lempainen J, Laine AP, Hammais A, Toppari J, Simell O, Veijola R, et al. Non-HLA Gene Effects on the Disease Process of Type 1 Diabetes: From HLA Susceptibility to Overt Disease. J Autoimmun (2015) 61:45–53. doi: 10.1016/j.jaut.2015.05.005

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Plagnol V, Howson JMM, Smyth DJ, Walker N, Hafler JP, Wallace C, et al. Genome-Wide Association Analysis of Autoantibody Positivity in Type 1 Diabetes Cases. PLoS Genet (2011) 7:e1002216. doi: 10.1371/journal.pgen.1002216

PubMed Abstract | CrossRef Full Text | Google Scholar

43. Li J, Sun X, Luo S, Lin J, Xiao Y, Yu H, et al. The Positivity Rate of IA-2A and ZnT8A in the Chinese Han Population With Type 1 Diabetes Mellitus: Association With Rs1143627 and Rs1143643 Polymorphisms in the IL1B Gene. Front Pharmacol (2021) 12:729890. doi: 10.3389/fphar.2021.729890

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Manning BD, Toker A. AKT/PKB Signaling: Navigating the Network. Cell (2017) 169:381–405. doi: 10.1016/j.cell.2017.04.001

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Pavlov VA, Chavan SS, Tracey KJ. Molecular and Functional Neuroscience in Immunity. Annu Rev Immunol (2018) 36:783–812. doi: 10.1146/annurev-immunol-042617-053158

PubMed Abstract | CrossRef Full Text | Google Scholar