Integrative Genomic Analysis Predicts Regulatory Role of N6-Methyladenosine-Associated SNPs for Adiposity

Genome-wide association studies have identified many susceptible loci to explore the genetic factors of adiposity. However, the specific mechanisms by which these SNPs (single nucleotide polymorphism), particularly in the non-coding region, are involved in the pathogenesis of adiposity remain unclear. Recently, genetic variation is thought to affect N6-methyladenosine (m6A) RNA modification, which is the most common post-transcriptional messenger RNA modification. In this study, we identified a large number of BMI (body mass index)-associated m6A-SNPs from published GWAS summary statistics through a public database and explored their potential mechanisms involved in the pathogenesis of adiposity. In summary, the integrative analysis detected 20,993 BMI-associated m6A-SNPs and 230 m6A-SNPs which reached the genome-wide suggestive threshold (5.0E-05), while 215 of them showed eQTL signals and 167 are the corresponding genes. The leading SNP rs8024 (C/A) was located next to the m6A modification site of 3′UTR of the IPO9 gene, which was predicted to affect the m6A modification site and regulate the expression of the IPO9 gene to participate in the pathogenesis of adiposity. This m6A-SNP/gene expression/adiposity triplets provide a new annotation for the pathogenic mechanism of adiposity risk loci identified by GWAS.


INTRODUCTION
Adiposity is considered to be both an independent disease and a clear risk factor that is closely related to the occurrence and death risk of non-communicable chronic diseases such as hypertension, cardiovascular and cerebrovascular disease, diabetes, and specific types of cancer, which has become one of the main sources of the burden of preventability worldwide (Apovian, 2016). Genome-wide association studies (GWAS) have identified a large number of potential risk loci and susceptible genes for exploring genetic factors of adiposity (Speakman et al., 2018). However, in addition to SNPs located in the coding regions of proteins which can directly affect amino acid sequences, many SNPs located in non-coding regions are also considered to affect epigenetic regulation and thus alter gene expression. FTO (fat mass and obesity associated gene) is the first gene associated with adiposity traits identified in the GWAS study (Frayling et al., 2007). Certain genetic variants of FTO gene appeared to be correlated with adiposity in human. Jia et al. (2011) revealed that N 6 -methyladenosine in nuclear RNA is a major substrate of FTO.
Similarly, our previous study also demonstrated that knockdown of m6A methylase METTL3 could promote the adipogenesis of bone marrow mesenchymal stem cells while inhibit osteogenesis, suggesting that m6A modification has an important role in adipogenesis and adiposity . m6A modification, first discovered in mammalian mRNAs in 1974, is considered to be the most abundant internal modification in eukaryotic mRNAs, and has been proven to be highly conserved in viruses, bacteria, yeast, plants, and vertebrates (Desrosiers et al., 1974;Perry and Kelley, 1974;Krug et al., 1976). m6A sites are mainly distributed in the stop codon regions, 3'UTR, 5'UTR and internal long exons of mRNA (Narayan et al., 1994). In recent years, the establishment of a sequencing method based on m6A antibody specific binding immunoprecipitation (MeRIP-Seq) proved that m6A modification was dynamically reversible, and made it possible to identify m6A modification sites across the whole-transcriptome (Dominissini et al., 2012;Meyer et al., 2012). m6A modification could alter the structure, alternative splicing, stability, and translation of mRNA, thereby regulates gene expression, and is widely involved in the regulation of stem cell pluripotency, cell differentiation, neural maturation, and embryo development (Fu et al., 2014;Roundtree et al., 2017;Zhao et al., 2017;Frye et al., 2018). Some genetic variations in specific regions may affect m6A methylation and subsequent biological processes by changing the RNA sequence of target sites or key flanking nucleotides. These genetic variations are known as m6A-SNP Lin et al., 2019). Recently, Xiao et al. found an association between m6A modification and gene expression homeostasis in the whole-transcriptome m6A methylation across 21 major fetal tissues, revealing that human SNP sites with eQTL effects were enriched near the m6A modified region. It was confirmed that the SNP near the m6A modification site could affect the m6A modification of RNA .
Therefore, m6A-SNP can be regarded as an important genetic functional variant, which can provide new clues for understanding the pathogenic molecular mechanisms of genetic variants. Till now, there are still no reports about the role of m6A-SNP in adiposity, and identification of m6A-SNP associated with adiposity could enrich the understanding of its genetic mechanism and the role of m6A modification in adiposity. Therefore, here we will use the public GWAS dataset to screen and explore the potential impact of m6A-SNP on adiposity.

Identify BMI-Associated m6A-Single Nucleotide Polymorphisms
By analysis of the BMI GWAS summary statistics with the m6A-SNP list, we identified the intersection of both, namely BMI-associated m6A-SNPs. Among them, BMI GWAS summary statistics is available on the GIANT website 1 , including a metaanalysis of previous GWAS of the GIANT consortium and GWAS of participants in UKBB, with a total sample size of 1 http://portals.broadinstitute.org/collaboration/giant/index.php/Main_Page approximately 0.7 million. The m6A-SNP list was downloaded from the m6Avar database 2 . The m6AVar database currently contained a high confidence of human m6A-SNP of 13,703 (miCLIP/PA-m6A-seq experiment), 54,222 medium (MeRIP-Seq experiment) and 245,076 (random genome prediction based on random forest algorithm) of low confidence. The 5.0E-05 was set as the suggestive threshold. In addition, we perform GO enrichment to explore the biological functions of these m6A-SNPs with an online website 3 (Zhou et al., 2019).

Functional Annotation of BMI-Associated m6A-SNPs
For the BMI-associated m6A-SNPs that have been identified, we used online tools to annotate them to further explore the mechanism by which they may play biological functions. The influence on the modification of nearby m6A site was obtained through m6AVar database query. One of the important potential pathways for these m6A-SNPs to play a biological role is affecting the post-transcriptional regulation of local genes. We have checked whether they have cis-eQTL effects in several different tissues and cells using the HaploReg browser (Ward and Kellis, 2016). In addition, we also queried about their other potential functions in transcription regulation, such as altering the binding capacity of transcription factors.

Prediction of m6A Modifications Near m6A-SNP
Sequence-based RNA adenosine methylation site predictor (SRAMP) is a powerful m6A modification prediction tool based on a random forest classifier using genomic or cDNA sequences as input, which can predict the confidence thresholds of m6A modifications at different loci 4 (Zhou et al., 2016). In order to determine the possibility of the identified m6A-SNPs affecting RNA modification, we used the online m6A modification prediction tool SRAMP with gene sequences of different genotypes as input to predict the changes of m6A modification.

Differential Gene Expression Analysis
We further queried three transcriptome-wide studies on adiposity in the GEO public database. To verify whether m6A-SNPs showing cis-eQTL signals could affect local gene expression in the pathogenesis of adiposity. GSE88837 includes the gene expression level of visceral adipose tissues extracted during abdominal surgeries on 30 lean and obese adolescent females. GSE109597 includes whole peripheral blood samples from 90 healthy controls and overweight participants. GSE70353 includes a large-scale population study of subcutaneous adipose tissue gene expression level in 770 men with different BMI values. BMI is a person's weight in kilograms divided by the square of height in meters. A high BMI can be an indicator of high body fatness, though there are slightly differences in BMI criteria for different genders and ages. Usually, the adults are regarded as  overweight if their BMI is between 25 and 30 kg/m 2 , and obese when it is greater than 30 kg/m 2 (Kuczmarski and Flegal, 2000). We used GEO2R online tool to calculate differentially expressed genes in adiposity (BMI > 30) and healthy (18.5 < BMI < 25) populations, and examined whether the BMI-associated genes (m6A-SNP with cis-eQTL) are differentially expressed between cases and controls. A significance level of p < 0.05 was used for differential expression analysis. Only m6A-SNPs with medium or high confidence and differential local gene expression were selected.

Identification of BMI-Associated m6A-SNPs
First, we analyzed 27381302 SNPs in BMI GWAS and 301529 m6A SNPs in the m6Avar database and identified 18568 m6A-SNPs on mRNA and 2425 m6A-SNPs related to m6A modifications on other RNAs (Figures 1, 2). These SNPs were not only located in the protein-coding gene regions, but also related to other RNA such as lincRNA, miRNA, and snoRNA. A large proportion of m6A-SNPs are distributed in exon regions (66%), including 54% in protein coding sequence and 12% in other exon regions. In addition, 29% of m6A-SNPs are distributed in 3 UTR, 5% of m6A-SNPs are distributed in 5 UTR, and a very small proportion of m6A-SNPs are distributed in the intron region (less than 1%) (Supplementary Figure S3).
For the 20993 unique m6A-SNPs found in the GWAS dataset, 713, 3559, and 16721 m6A-SNPs were belonged to the high, medium, and low confidence region categories, respectively. Among them, the high-confidence m6A-SNP included variants that disrupt m6A motif obtained from miCLIP/PA-seq. For the medium-confidence m6A-SNP, variants that changed the sequence features for m6A modification in MeRIP-seq were included. The low-confidence m6A-SNP included variants near the genome-wide prediction based on random forest algorithm of m6A sequences. We used 5.0E-05 as the suggestive threshold to screen the association between these 20,993 m6A-SNPs and adiposity. Of these SNPs, 230 appeared to be associated with adiposity (p < 5.0E-05) (Supplementary Table S1). Through GO enrichment analysis, we found that the gene functions of these m6A-SNPs were enriched in biological processes related to transcriptional regulation, such as transcription factor binding, histone binding and methyltransferase activity, suggesting that these genes may affect transcriptional regulation (Supplementary Figure S1).

Functional Annotation of BMI-Associated m6A-SNPs
To further explore the potential functional mechanisms of the 230 m6A-SNPs associated with adiposity, we investigated whether they were related to local gene expression level. In total, 215 BMI-associated m6A-SNPs tested (p < 0.05) showed cis-eQTL signals with 167 corresponding genes ( Supplementary  Table S1). Interestingly, by retrieving these m6A-SNPs in genotype-tissue expression (GTEx) database (Battle et al., 2017), we found that m6A-SNPs cause loss-of-function of m6A are more inclined to reduce local gene expression, and those gain-of-functions seem to have the reversed eQTL effect FIGURE 3 | (A) Regional association plots of the rs8024 locus. The regional plot was plotted from published data of the BMI GWAS using the locus Zoom website tool (http://locuszoom.sph.umich.edu//). SNP rs8024 (purple) had a high degree of linkage disequilibrium (r 2 > 0.8) with other high-risk SNP sites associated with BMI (p = 7.10E-39). (B) Integrative analysis of the potential function of SNP rs8024 in adiposity by querying USCS browser. rs8024 is located at the 3 UTR of IPO9. This region shows high transcription level and DNaseI hypersensitivity. RIP-chip GeneST from ENCODE/SUNY Albany data showed that the RNA-binding protein IGF2BP1 might have a potential interaction with rs8024-containing transcripts. Figure S3). However, the mechanism of m6A modification affecting gene expression is very complicated, the impact of loss/gain of function of the m6A-SNPs on local gene expression is still remained to be explored. At the same time, to explore other possible mechanisms of these SNPs affecting transcriptional regulation, we queried the Haploreg browser and found that 20 of them changed the binding of one or more proteins in different cell types from the ENCODE transcription factor ChIP-seq dataset, and some of these transcriptional factors have been proved to be related with adiposity. For example, rs1802036 is related to the binding of transcriptional factor GATA1, which is considered to be a specific gene of brown adipose tissue (BAT). The increased expression of GATA1 is related to BAT formation, which may affect fat accumulation and lead to adiposity (Lu et al., 2012), suggesting that these SNPs may act as multifunctional variants to affect the transcriptional regulation of local genes ( Table 1 and Supplementary Table S1).

Differential Expression Analysis
For the above 167 m6A-SNPs showing cis-eQTL signals, we compared the expression level of local genes in three transcriptome studies of obese population. In these three datasets, we found that 88 of the 167 genes were differentially expressed in at least one study (p < 0.05) (Supplementary Table S1).
For example, the established BMI-gene ADPGK expression is gradually downregulated in healthy, overweight, and adiposity populations (GSE109597, p = 0.00708315). We believe that m6A-SNP may be involved in the occurrence of adiposity by affecting the expression level of corresponding genes.

DISCUSSION
More than 100 chemical modifications have been detected on mRNA, of which N 6 -adenine (m6A) was one of the most extensive epigenetic modifications. During biological process such as cell differentiation, embryonic development, and stress, certain mRNAs can undergo modifications including N 1 -adenine methylation, N 5 -cytosine methylation, pseudouracil, and N 6adenine methylation (Boccaletto et al., 2018). Together, they form an apparently modified transcriptome for post-transcriptional regulation of mRNAs, enabling precise timing control of the process of mRNA translation into proteins, especially m6A modification that can regulate a series of biological processes in cells by regulating mRNA metabolism and translation. Recent evidence suggested that m6A methylation could be involved in regulating blood glucose and lipid metabolism, and that abnormalities in m6A modification were increasingly linked to a variety of metabolic diseases (Peng et al., 2019;Wu et al., 2017;Zong et al., 2019). Nutritional challenges or dietary supplements have shown an effect on m6A methylation level, altering gene expression and phenotype. At the same time, some studies have demonstrated that the reduction of the overall level of m6A modification obtained by knocking down m6A methylase METTL3 could promote adipogenic differentiation (Yao et al., 2019). However, there is no report on the role of m6A-SNP in adiposity.
In this study, we found several m6A-SNP loci with high confidence that may be involved in adiposity. We used rs8024 located in the 3 UTR region of the IPO9 gene as an example for subsequent analysis. A recent study has shown that the IPO9 could interact with the JNK1/2 and promote their nuclear translocation (Zehorai and Seger, 2019). JNK1/2 has been proved to be an inhibitor of adipogenesis (Engin, 2017). Genetic variation in IPO9 may disrupt the nuclear localization of JNK1/2 and promote adipogenesis, thus leading to adiposity. rs8024 was located on the 3 UTR of the IPO9 gene exon on chromosome FIGURE 5 | Expression levels of selected genes were displayed among adiposity and healthy populations in GSE88837 (visceral adipose tissue), GSE109597 (peripheral blood sample), and GSE70353 (subcutaneous adipose tissue).
1, and had a high degree of linkage disequilibrium (r 2 > 0.8) with other high-risk SNP sites nearby (Figure 3). By querying the m6AVar database, we found that AGO2, EIF4A3, PTBP1, IGF2BP2, and IGF2BP3 were binding to the rs8024-related RNA region in IPO9 transcript in the CLIPdb and starBase2 databases. Among them, the IGF2BP family has recently been proven to be an m6A reader, which was believed to stabilize mRNA by binding to the m6A site on the mRNA. Knockdown of IGF2BP brought about a down-regulation of target gene expression (Huang et al., 2018). Therefore, rs8024 might lead to the loss of nearby m6A modifications, making the corresponding m6A reader unable to recognize this site and reducing gene expression. On the UCSC genome browser, we also found that the transcript of IPO9 had a binding site to the RNA-binding protein IGF2BP1 (Figure 3).
In addition, possible m6A methylation sites IPO9 transcript sequence was predicted on the SRAMP website. After input of the IPO9 reference sequence, we found a highly convincible m6A modified predicted peak near rs8024, but this predicted peak disappeared after entering the altered sequence of IPO9, which also proved that rs8024 would affect m6A modification (Figure 4). To explore the mechanism in which rs8024 affects adiposity, we performed a phenotypic study on GWAS Atlas 5 and found that rs8024 was also related to also related to daytime napping (p = 3.21E-20) and hypertension (p = 8.65E-12) (Supplementary Figure S2 and Supplementary Table S2).
Except for affecting m6A modification to change gene expression, m6A-SNP may also have other ways to participate in the regulation of gene expression. By querying the Haploreg database, we found that rs8024 was also located in the DNaseI hypersensitivity cluster region and would affect the binding of three transcription factors. Therefore, rs8024 should be considered a multifunctional variant involved in adiposity. As for expression level, it was found that mRNA expression of more than half of the m6A-SNP-associated genes was statistically different in at least one dataset (Figure 5 and Supplementary  Table S1). In conclusion, we demonstrated that the identified m6A-SNPs may affect adiposity through giving rise to abnormal gene expression.

CONCLUSION
In this study, we identified a total of 230 m6A-SNPs that may be associated with adiposity, and explored their potential functions using public databases. However, further research is still needed to verify how these m6A-SNPs affect m6A modification and its practical impact on the pathogenesis of adiposity. It is known that besides influencing global gene expression level, m6A-SNPs could also involve in disease progression by affecting the ratio between different RNA isoforms and the translation level of their protein products. However, such data are currently publicly limited, so only the potential impact of m6A-SNPs on gene expression level was explored in this study. Moreover, with the development of different detection m6A technologies, such as antibody-free m6A sequencing and nanopore sequencing, more accurate, disease-specific m6A sites would be identified. Genetic variants found near these sites have more possibilities to influence the pathogenesis of human diseases, so there are still more potential m6A-SNPs needed to be determined.

DATA AVAILABILITY STATEMENT
The datasets generated for this study can be found in the Data of BMI GWAS were publicly available at GIANT consortium (http://portals.broadinstitute.org/collaboration/giant/index.php/ Main_Page). The m6ASNPs data can be downloaded from the m6AVar database (http://m6avar.renlab.org/). The three public gene expression data sets, GSE88837, GSE109597, and GSE70353 can be downloaded from GEO database (https://www.ncbi.nlm. nih.gov/geo/).

AUTHOR CONTRIBUTIONS
WL designed the study and wrote the manuscript. HX generated the figures and table. QY checked and polished the language. SZ designed the study and reviewed the manuscript.