Exploring the Epigenetic Regulatory Role of m6A-Associated SNPs in Type 2 Diabetes Pathogenesis

Introduction

With the continuous increase in obesity and life expectancy in recent decades, the global incidence of diabetes has also risen rapidly.¹ Type 2 diabetes (T2D) is the most common type of diabetes, accounting for approximately 90% of all diabetes cases.² T2D is a chronic non-infectious disease characterized by impaired islet beta-cell function, insulin resistance and the resultant hyperglycemia. Besides the common environmental factors, cumulative evidence indicates that genetic factors also play a decisive role in the occurrence and progression of T2D. Clarifying the involvement of genetic variants in the development of T2D is critical to understand its pathogenesis. Indeed, many endeavors have been made to investigate T2D susceptibility genes over the past years.^3–5 Numerous susceptibility loci for T2D have been identified by a series of genome-wide association studies (GWAS).^6–8

More than 400 genomic regions have been identified as affecting T2D.⁹ For example, variations in hematopoietically expressed Homeobox (HHEX), the Solute Carrier Family 30 Member 8 (SLC30A8), and the transcription factor-7-like2 (TCF7L2) gene regions were the earliest identified susceptibility loci in GWAS related to T2D. Most of these risk genes were related to glucose metabolism and lipid synthesis.^10–12 Unfortunately, besides a small number of the published T2D-associated single nucleotide polymorphisms (SNPs) locating in protein coding regions, most of these SNPs are located in non-coding regions, making it more difficult to explain their biological effects.^13,14

N⁶-methyladenosine (m6A) is the methylation of the adenosine at the sixth position of the nitrogen, which is the most common post-transcriptional RNA modification in eukaryotes. With the landmark study that identified the first m6A demethylase fat mass and obesity-associated protein (FTO) in 2011 and the transcriptome-wide profiling of m6A by high-throughput sequencing,¹⁵ precise mapping of m6A distribution on transcript has been achieved. These m6A peaks have a conserved modification motif RRACH (R for A or G, H for A, U or C).¹⁶ In mRNA transcripts, the recognition of modification sites by m6A methyltransferase is highly dependent on these m6A-specific sequences, while genetic variation may change these m6A-specific recognition sequences on local gene transcripts, thus affecting the modification and recognition of m6A. Through the combination of m6A sequencing and whole genome sequencing (WGS) of large-scale cohort, Xiong and Zhang et al proved that genetic variation is an effective driver of m6A modification, respectively.^17,18 Recently, the correlations between m6A RNA methylation-related SNPs (m6A-SNPs) and hereditary diseases were widely explored.^19–23 In studies of Mo et al, they identified a large number of m6A-SNPs affecting a series of phenotypes related to cardiovascular diseases, such as blood pressure, blood lipids and coronary heart disease.^24–26 In the field of neurology, m6A-SNPs have been verified to be closely associated with stroke and multiple sclerosis.^27,28 Besides, m6A-SNPs could also participate in cancer development. m6A-SNP rs178184 was shown to affect the occurrence of colon cancer by regulating the expression of NOVA1 gene.²² Our previous study proved that m6A RNA methylation-related SNPs might be involved in adiposity, which had a close causal relationship with T2D.²⁹ These findings are crucial for a better understanding of the role of m6A modification in pathological changes caused by genetic variation. However, no study has been reported on the role of m6A-SNPs in T2D pathogenesis until now. In this study, we explored the association of m6A-SNPs with T2D and demonstrated the potential functions of these m6A-SNPs in T2D pathogenesis.

Materials and Methods

Determination of m6A-SNPs for T2D

To identify the T2D-associated SNPs that probably affect RNA m6A modification, we retrieved a recently released large-scale T2D GWAS summary statistics among 81412 T2D cases and 370832 controls.⁶ The data is publicly available on the DIAGRAM website (http://diagram-consortium.org/index.html). We identified the potential m6A-SNPs by comparing the m6A-SNP list in m6AVar database with T2D summary statistics.³⁰ The identified m6A-SNPs were annotated, and a preliminary screening was performed with a threshold of 5.0e-05 for subsequent analysis. GO enrichment analysis of m6A-SNPs passing threshold was performed through an online website metascape (http://metascape.org/gp/index.html).³¹

Integrative Analysis of T2D-Associated m6A-SNPs

It has been confirmed that m6A modification can affect the mRNA splicing, nuclear export, stability and degradation, thus affecting local gene transcription, which might be an important mechanism for m6A-SNPs to involve the pathogenesis of T2D. To investigate whether these genetic variants could lead to changes in gene transcription, we retrieved their expression quantitative trait loci (eQTL) effects through the HaploReg browser (http://pubs.broadinstitute.org/mammals/haploreg/haploreg.php).³² Further, we also searched for several other potential regulatory mechanisms of these SNPs on transcriptional regulation, such as altering protein binding and changing motif. The results are shown in Table 1 and Supplementary Table 1.

Table 1 The 20 Most Important m6A-SNPs Related to T2D

Prediction of m6A Modification Near m6A-SNPs

Among the T2D-associated m6A-SNPs identified above, we predicted the potential m6A modification site on the transcript by inputting the reference and altered sequence to an online m6A modification prediction tool, SRAMP (http://www.cuilab.cn/sramp/).³³ SRAMP is a website-based computational predictor for identifying mammalian RNA m6A modification sites. SRAMP includes three random forest classifiers and uses genomic sequence or cDNA sequence as input. To identify whether m6A-SNP could affect the nearby m6A modifications, the altered sequence and reference sequence were input into SRAMP, and the confidence of m6A modifications was predicted.

Differential Expression Analysis

For m6A-SNPs that showed eQTL signal, we further tried to determine whether the related gene expression of these T2D-associated m6A-SNPs genes in T2D was changed. Based on the transcriptomic data available in the gene expression omnibus (GEO) database, we displayed the expression level of the identified T2D-associated genes (m6A-SNPs showed eQTL signal). Three datasets GSE23343, GSE38642 and GSE78721 containing standardized T2D-associated gene expression signals were retrieved and downloaded from GEO database (http://www.ncbi.nlm.nih.gov/geo). Samples in GSE23343 were collected from 7 subjects with normal blood glucose levels and 10 patients with T2D. GSE38642 contained the gene expression profile in 63 islets from cadaver donors (including 9 diabetic and 54 non-diabetic). And GSE78721 contained gene expression level of adipose tissues from 63 patients suffering from T2D and 68 age/BMI matched normal glucose tolerance controls. Considering that multiple tests are used among the three datasets, we use the Bonferroni adjustment to set a more stringent threshold, p < 0.0167 (0.05/3) was set as a threshold for defining differential genes. Also, to prove that these gene expression differences were mediated by m6A-SNPs instead of the T2D phenotype, we downloaded the transcriptome data of tissues that are not closely related to the T2D phenotype as control sets, including enteroendocrine cells (GSE132831) and peripheral blood monocytes (GSE156061). We also explored the expression of these differential genes in control sets.

For single-cell RNA sequencing (scRNA-seq) data analysis, we downloaded the islet scRNA-seq data of high-fat diet-induced diabetic mice and healthy controls in GSE162512, and used the Seurat package to load the hdf5 format data for subsequent process.³⁴ Based on the expression of classic cell markers, we defined the cell types and compared the expression of HLA-DQB1 (H2-Ab1) in different cell populations.

Quantitative Reverse Transcription Polymerase Chain Reaction (qRT-PCR)

A total of 8 db/db mice and 8 heterozygous mice at the age of 8 weeks were used. Male db/db mice were used as animal model of T2D, and the heterozygous mice were used as control. The islets of the mouse pancreas were isolated as previously described.³⁵ In brief, the pancreatic tissue was cut into small pieces of 1~2 mm and digested in 1 mg/mL collagenase XI (Sigma) for 20 minutes. Ficoll solution was used to purify the islets by centrifugation. Total RNA of the obtained islets was extracted using Trizol reagent (Invitrogen), and then reversely transcribed to obtain stable cDNA using PrimeScript^TM RT reagent Kit with gDNA Eraser (TaKaRa Bio). The qRT-PCR was performed using SYBR Premix Ex Taq II (TaKaRa Bio) in Quant Studio^TM 3 real-time fluorescent quantitative PCR instrument (ThermoFisher Scientific). Glyceraldehyde 3-phosphate dehydrogenase (Gapdh) was used as an internal reference to normalize the gene expression. The result was calculated using the 2^−ΔΔCt method and expressed as a multiple change relative to Gapdh. The primer sequences were summarized in Supplementary Table 2. The animal studies were approved by the Animal Care and Ethics Committee of West China Hospital of Stomatology, Sichuan University. The animal experiments complied with the ARRIVE guidelines and were carried out in accordance with the National Institutes of Health guide for the care and use of Laboratory animals (NIH Publications No. 8023, revised 1978).

Results

Determination of T2D-Associated m6A-SNPs

We identified 15,124 m6A-SNPs and 16,420 potentially affected m6A sites (Figure 1). For these 15,124 T2D-associated m6A-SNPs identified, 11820, 2950 and 354 m6A-SNPs belong to the low, medium and high confidence categories, respectively. p < 5.0e-05 was defined as the significance threshold and 71 appeared to be associated with T2D (p < 5.0e-05) (Figure 2A). GO enrichment analysis showed that these 71 genes were mainly enriched in antigen processing and presentation related biological processes (Figure 2B). Among these m6A-SNPs, 56, 12 and 3 belong to the low, medium and high confidence intervals, respectively. 21 SNPs are non-synonymous, 20 are synonymous mutations, and 28 are located at untranslated regions (UTR) (Supplementary Table 1). As for other potential regulatory mechanisms of these m6A-SNPs, we found that 52 altered the motifs and 19 influenced the binding of regulatory proteins in different cell types from the ENCODE transcription factor ChIP-seq datasets (Supplementary Table 1).

Figure 1 Design and main implementation steps of this study.

Figure 2 Identification of T2D related m6A-SNPs. (A) The identified T2D related m6A-SNPs with genome-wide significance are displayed in the Manhattan plot. P < 5.0e−05 was set as the suggestive threshold, and p < 5.0e-08 was set as the genome-wide threshold. (B) GO enrichment analysis of the T2D related m6A-SNPs (p < 5.0e-05). (C) Gene expression differences between T2D patients and controls in islet (GSE38642), adipose (GSE78721), and liver tissue (GSE23343).

Identification of Differentially Expressed Genes

For the 56 T2D-associated m6A-SNPs with eQTL signals, we analyzed whether the related gene expression changed in published datasets about the gene expression of the tissues from T2D patients and controls. The etiology of T2D is complicated, and the abnormal function of organs related to glucose metabolism may lead to blood glucose homeostasis disorders. For instance, liver dysfunction will lead to a decline in glycogen synthesis capacity, resulting in an imbalance in blood glucose homeostasis.^36,37 And the excessive accumulation of adipose tissue, especially visceral adipose tissue, can lead to chronic inflammation and promote the progression of T2D.^38,39 Also, the dysfunction and loss of β-cells that secrete insulin in the endocrine pancreas can also cause T2D.⁴⁰ Therefore, the transcriptomic data of these three tissues were selected in the GEO database to validate the expression change of the m6A-SNPs related genes in T2D patients. We found that 11 of these genes were differentially expressed in at least one dataset (Supplementary Table 1). Then, the expression level of representative genes with medium or high confidence were selected to be shown (Figure 2C). Islet beta-cells can secrete insulin to regulate blood glucose, and impaired beta-cell function can lead to hyperglycemia caused by insufficient insulin secretion. In GSE38642, HLA-DQB1 (p = 0.00365124) and KCNJ11 (p = 0.00538383) were found to be differentially expressed in islet tissue. Insulin resistance caused by obesity also leads to T2D. PPARA (p = 0.00839) and MFHAS1 (p = 0.00067) were found to be abnormally expressed in GSE78721 (adipose tissue). In addition, the liver can respond to blood glucose concentration and regulate the secretion of insulin or glucagon, thereby regulating the homeostasis of glucose. In GSE23343 (liver tissue), RREB1 (p = 0.0107675) was differentially expressed. Therefore, we found that 11 m6A-SNP-related genes were differentially expressed between T2D patients and healthy controls. The abnormal expression of these genes may be involved in the pathogenesis of T2D. In the control sets (enteroendocrine cells and peripheral blood monocytes), we found that 8 of 11 were also differentially expressed in at least one control set (Supplementary Table 3). Therefore, we believe that the differential expression of these genes is caused by m6A modification and not related to T2D phenotype.

To validate whether the top m6A-SNP-related genes showed expression changes during T2D pathogenesis, we detected the expression of these genes in the pancreatic islets of db/db mice and heterozygous control by qRT-PCR. The results showed that the expression of 8 genes was found to be significantly different between the diabetic group and control. Although the expression of Mfhas1, Ppara and Rreb1 did not show a significant difference, their average expression levels were the same as the trend of these genes in GEO datasets (Supplementary Figure 1). In short, the results of qRT-PCR proved that the expression of these m6A-SNP-related genes changed under the diabetic condition.

Integrative Analysis of T2D-Associated m6A-SNPs

To further elucidate the potential functional mechanisms of the identified 71 T2D-associated m6A-SNPs, we further investigated whether the 71 T2D-associated m6A-SNPs were related to the expression level of local genes and the underlying functional mechanism. In all, 56 T2D-related m6A-SNPs displayed eQTL signals related to local genes expression. Among these m6A-SNPs, 19 of them might change the binding of transcription factors in different cell lines, and 52 of them might change motifs (Table 1 and Supplementary Table 1). From a transcriptional influence standpoint, the leading SNP rs4993986 was located at the 3ʹ-UTR of HLA-DQB1 with DNaseI hypersensitivity and ENCODE transcription factor clusters binding signals. RNA-binding protein (RBP) analysis revealed that the RBP PABPC1 showed a binding region to rs4993986 (Figure 3A). According to the results of SRAMP prediction, the original m6A modification peak disappeared (labeled in red) when the altered gene sequence was input (Figure 3B), and the RNA secondary structure prediction showed that rs4993986 (C>G) was located nearby the m6A modification site (Figure 3C). In short, m6A-SNPs may affect the expression of local genes at the transcriptional level through the above-mentioned multiple pathways.

Figure 3 Integrated analysis of the leading m6A-SNP rs4993986. (A) The protein coding region of rs4993986 in HLA-DQB1 gene sequence shows high transcriptional activity and DNaseI hypersensitivity. RNA-binding protein analysis revealed that rs4993986 exhibited a binding region to the RBP PABPC1. (B) Prediction of m6A modification on the HLA-DQB1 transcript (ENST00000434651.6) on the SRAMP website. (C) The secondary structure of HLA-DQB1 transcript and the position of rs4993986 near the m6A modification site. (D) The expression of HLA-DQB1 (H2-Ab1) was identified in the mouse pancreas scRNA-seq data. Violin plot shows that HLA-DQB1 expression is increased in diabetic islet acinar cells. (E) HLA-DQB1 expression of different cell types in islet of diabetes and control group.

Furthermore, in the scRNA-seq dataset of islet from diabetic mice induced by high-fat diets, we found that HLA-DQB1 was highly expressed in two cell populations (macrophages and acinar cells). The expression of HLA-DQB1 was increased in the acinar cells of diabetic group (Figure 3D and E), which is consistent with the up-regulation of HLA-DQB1 expression in diabetic patient islet.

Discussion

It has been found that more than 170 RNA modifications are involved in physiological processes and disease progression,⁴¹ among which m6A modification is considered to be the most abundant modification in mRNA.⁴² m6A modification has been reported to exert an important role in blood glucose homeostasis. Conditional knockout of the m6A methyltransferase Mettl14 led to decreased insulin secretion in β-cells and blood glucose intolerance.^43,44 In addition, the hepatocyte-specific knockout of Mettl3, which is also a key component of the m6A methyltransferase complex, improved the mouse hepatic insulin sensitivity and inhibited the fatty acid metabolism.⁴⁵ However, no studies have yet found the contribution of m6A-SNPs to the pathogenesis of T2D. In the current study, a great amount of T2D-associated m6A-SNPs have been identified, which are verified by other published datasets at genes expression level.

m6A modification of RNA can be recognized by different m6A readers, which may improve the stability of RNA transcripts (eg IGFBPs),⁴⁶ or promote the degradation of transcripts (eg YTHDF2).⁴⁷ RNA methylation has been reported to regulate the entire mRNA metabolism process, including RNA stability, splicing, nuclear transport, translation, and degradation, thereby affecting the key biological functions in the regulation of different cellular processes.^48,49 Therefore, the expression of genes corresponding to m6A-SNPs may be affected. Experimental epitranscriptomic study conducted by Tian and his colleagues proved that overexpression of the rs8100241[A] (an m6A-SNP, G>A) allele significantly increased the level of ANKLE1 m6A and promoted the expression of ANKLE1 protein compared with that of rs8100241[G] allele.⁵⁰ To validate whether the screened m6A-SNP-related genes were involved in the diabetes pathogenesis, we retrieved the GEO database and selected three tissues (liver, adipose, and pancreatic islets) from T2D patients to explore the expression of m6A-SNP-related genes under diabetic status. 11 m6A-SNP related genes were differentially expressed between T2D patients and the controls in at least one dataset. Also in diabetic mice, the expression trend of the m6A-SNPs related genes was verified by qRT-PCR. HLA-DQB1 encoded by the leading m6A-SNP rs4993986 was significantly upregulated in islets of diabetic patients and islet acinar cells in diabetic mice. With a minor allele frequency of 0.3691, rs4993986 was considered as a common SNP. Human leukocyte antigen (HLA) class II molecules are the gene loci found to be most highly associated with autoimmune diabetes mellitus.⁵¹ Increased expression of HLA-DQB1 results in beta cell failure mediated by autoimmunity.⁵² The structural changes of the HLA-DQB1 gene also participated in T2D pathogenesis,⁵³ indicating that genetic variations in HLA-DQB1 locus were closely related to the occurrence of T2D. To determine the cell-type specific contribution of these m6A-SNPs, we further analyzed the scRNA-seq data of islet tissue in diabetic mice, and found that HLA-DQB1 was abnormally expressed in acinar cells under diabetic condition, which might participate in the progression of diabetes.

In order to avoid a large number of false-positive results, traditional GWAS usually sets a strict significance threshold (usually p < 5.0e-08) to screen for susceptibility loci.⁵⁴ With the traditional significance threshold, the moderately correlated variations found in GWAS data can be ignored, and further research is needed to determine if they are involved in the disease.⁵⁵ Thus, a relatively loose threshold of p < 5.0e-05 was set to avoid missing some valuable SNPs in this study.²⁶ Exploring the potential functions of identified m6A-SNPs as multifunctional variants is essential for the study of T2D pathogenesis. In addition to affecting m6A modification, genetic variation could also alter gene expression by affecting transcriptional regulation factor binding, changing motif and RNA secondary structure. Therefore, the biological role of m6A-SNP should be considered as a multifunctional variation. By querying the HaploReg browser, we found several other potential regulatory mechanisms for these m6A-SNPs to affect gene transcription. For instance, rs4993986 was located in the DNaseI hypersensitivity cluster region, which may affect the binding of four proteins and alter six motifs. rs4993986 had a binding region with RBP PABPC1, which could be recruited to m6A modification sites by the m6A reader protein IGFBP2, and increased the transcript stability.⁵⁶ Therefore, rs4993986 should be considered as a multifunctional variant involved in the pathogenesis of T2D. In addition, the possible m6A methylation site is predicted by a sequence of the HLA-DQB1 transcript on the SRAMP website. A m6A modification peak with a very high confidence at the 3ʹ end of the HLA-DQB1 transcript near the SNP rs4993986 is predicted, suggesting that rs4993986 is likely to affect the modification of m6A methylation of this site.

In conclusion, we made the first endeavor to explore the relationship between T2D and m6A-SNPs in this study. 71 m6A-SNPs were potentially associated with T2D and the potential functions were identified. For example, rs4993986 has been shown to be located near a m6A modification at the 3ʹ end of the HLA-DQB1 transcript, which may participate in the pathogenesis of T2D by altering HLA-DQB1 expression. Gene therapy targeting HLA-DQB1 locus and m6A demethylation enzymes may alleviate the progression of T2D. Also, there are still several limitations in the current study. We speculated that m6A-SNPs participated in T2D pathogenesis by affecting local gene expression. However, m6A modification has also been proved to affect translation process.⁵⁷ Therefore, the validation of the screened m6A-SNPs at protein expression level has yet to be performed. As the exploration of RNA epigenetic regulatory networks increases, more m6A modifications in non-coding regions will be detected. Genetic variations near these sites may also be considered as m6A-SNPs to affect the pathogenesis of T2D, thus more m6A-SNPs that may affect non-coding regions remain to be identified. Further study targeting these m6A-SNPs to alleviate the development of T2D may become a new direction.