Genome-Wide Identification of m6A-Associated Single-Nucleotide Polymorphisms in Colorectal Cancer

Background

Colorectal cancer (CRC) is the third most deadly cancer worldwide, after breast and lung cancer.¹ Etiological studies have shown that genetic factors play an important role in CRC, and several genes were associated with CRC pathogenic, including MLH1, MSH2, PMS2, EPCAM, APC, and MUTYH.² However, these common genetic variants account for only 6% of CRC cases, which suggests a higher genetic complexity of CRC beyond these genes. Previous genome-wide association studies (GWAS) have identified several single-nucleotide polymorphisms (SNPs) related to CRC.^3,4 For instance, one study reported susceptibility variants at 8q23.3 (rs16892766) and 8q24.21 (rs6983267) that were associated with advanced-stage tumors and familial history of CRC.⁴

N6-methyladenosine (m6A) modification is a critical regulator of multiple cytopathological processes, including nuclear export, translation, splicing, and stability of mRNAs.⁵ Emerging evidence has shown that m6A modification plays a critical role in CRC; for example, METTL14 regulates m6A-dependent primary miR-375 processing to inhibit CRC progression,⁶ whereas METTL3 facilitates the progression of CRC via m6A-IGF2BP2/3-dependent mechanisms.^7,8 It was also confirmed that disease-associated genetic variants can influence m6A methylation by changing the RNA sequences or key flanking nucleotides of its target sites, suggesting that m6A-SNPs might affect mRNA stability, thereby contributing for the development of diseases.^9,10

Recently, m6A-SNPs have attracted considerable attention, and a number of prioritized SNPs have been identified by integrative analysis of cancer-related GWAS summary data, including in pancreatic, bladder, and gastric cancers.^11–13 The aim of the present study was to shed light and explore the potential contribution of m6A-SNPs in CRC pathogenesis.

Methods

Identification of CRC-Associated m6A-SNPs

The m6AVar database can provide host variants associated with m6A and may boost further mechanistic studies of genetic variants affecting m6A modifications.¹⁴ Therefore, CRC-associated m6A-SNPs were identified by integrating the data from GWAS and the m6Avar database. This study used publicly available GWAS data for CRC. The binary traits of GWAS datasets comprised 4562 CRC patients and 382,756 controls, relevant information was available from the link: ftp://share.sph.umich.edu/UKBB_SAIGE_HRC/PheCode_153_SAIGE_MACge20.txt.vcf.gz.The m6A-SNP list was downloaded from the m6Avar database (http://rmvar.renlab.org/download.html). Potential CRC-associated m6A-SNPs were identified through comparison of the SNPs in the GWAS datasets¹⁵ and the list of m6A-SNPs,¹⁶ which reached the genome-wide suggestive threshold (p < 0.001). Because SNP loci of genes encoding m6A regulators, including METTL3, METTL14, METTL16, WTAP, VIRMA, RBM15, FTO, ALKBH5, YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, IGF2BP1, IGF2BP2, IGF2BP3, HNRNPA2B1, and EIF3A were absent in the m6Var database, so these m6A-SNPs were identified, which reached the genome-wide suggestive threshold (p < 0.05).

Expression Quantitative Trait Loci (eQTLs) Analysis of CRC-Associated m6A-SNPs

Cis-acting eQTL (cis-eQTL) analysis can be used to evaluate the potential function of the m6A-SNPs that showed cis-eQTL signals in transcription regulation, such as altering protein binding, changing motifs, and affecting deoxyribonuclease.¹⁷ Therefore, cis-eQTL analysis was performed to investigate which CRC-associated m6A-SNPs could affect the RNA modification using the HaploReg browser (https://www.encodeproject.org/software/haploreg/).

Differential Expression Analysis

The corresponding genes of the identified eQTL m6A-SNPs were further evaluated according to differential expression among CRC patients and controls.¹⁸ Hence, three microarray gene expression CRC datasets (GSE89076, GSE21510, and GSE32323) publicly available in the Gene Expression Omnibus database (http://www.ncbi.nlm.nih.gov/geo) were used. Then, we used GEO2R online tool (https://www.ncbi.nlm.nih.gov/geo/geo2r/) to examine whether the CRC-associated genes (m6A-SNP with cis-eQTL) are differentially expressed between CRC and controls. A significance level of p < 0.05 was used for differential expression analysis.

Results

CRC-Associated m6A-SNPs

A total of 402 m6A-SNPs were extracted from the raw GWAS data by comparing the SNPs identified from the GWAS datasets and the m6A-SNPs from the m6AVar database (Figure 1). Next, the SNPs within genes encoding m6A regulators, including METTL3, METTL14, METTL16, WTAP, VIRMA, RBM15, FTO, ALKBH5, YTHDC1, YTHDC2, YTHDF1, YTHDF2, YTHDF3, IGF2BP1, IGF2BP2, IGF2BP3, HNRNPA2B1, and EIF3A, were collected as these SNP loci were absent in the m6Var database. This analysis revealed that rs112126539 (in IGF2BP1) could be a CRC-associated m6A-SNP (p = 0.018) (Table 1).

Table 1 Rs112126539 (in IGF2BP1) Could Be a CRC-Associated m6A-SNP

Figure 1 Flow chart of study design and analysis.

eQTL Analysis

eQTL analysis was performed on the 402 m6A-SNPs, which revealed that 98 of these m6A-SNPs had eQTL signals, with 76 and 22 m6A-SNPs having lost and gained modification functions, respectively (Table S1). Moreover, the rs2037844 polymorphism showed a stronger eQTL signals (p = 3.60×10⁻⁸) as compared to the other SNPs (the second was rs2957748; p = 1.17×10⁻⁴) (Figure 2A). Additionally, 98 CRC-associated m6A-SNPs showed eQTL signals displaying a unique distribution pattern, with most found in the intron and exon regions, and with few being in the coding sequences, 3ʹ/5ʹ-untranslated regions, or stop codon regions (Figure 2B).

Figure 2 Manhattan plot of genome-wide identified CRC-associated m6A-SNPs (A) and 98 CRC-associated m6A-SNPs showed eQTL signals displaying a unique distribution pattern (B).

Differential Expression Analysis

Three microarray datasets containing gene expression signals in CRC were analyzed: GSE89076, comprising paired normal and tumor tissues from 275 CRC patients; GSE21510, comprising 148 microarray datasets obtained from CRC tissue samples; and GSE32323, comprising 17 pairs of cancer and non-cancerous tissues from CRC patients.

Among the 98 CRC-associated m6A-SNPs that showed eQTL signals, three genes harboring m6A-SNPs (rs178184 in neuro-oncological ventral antigen 1 [NOVA1], rs35782901 in hydroxytryptamine receptor 4 (HTR4), and rs60571683 in solute carrier organic anion transporter family member 1B3 (SLCO1B3) were found to be differentially expressed in at least one CRC dataset (p < 0.05). NOVA1 (p = 3.41×10⁻¹¹) and HTR4 (p = 5.56×10⁻⁷) were significantly downregulated in GSE89076, whereas SLCO1B3 was significantly overexpressed in GSE32323 (p = 3.27×10⁻⁵), GSE21510 (p = 1.09×10⁻⁶), and GSE89076 (p = 7.63×10⁻⁶) (Figure 3 and Table 2). Therefore, these findings suggest that these three m6A-SNPs (rs178184, rs35782901, and rs60571683) may alter the expression of NOVA1, HTR4, and SLCO1B3, and subsequently impact on CRC pathogenesis.

Table 2 Three Genes of the m6A-SNPs (rs178184, rs35782901, rs60571683) Were Differentially Expressed Between Controls and CRC in at Least One Data Set

Figure 3 Expression levels of selected genes were displayed among controls and CRC of GSE32323, GSE21510, and GSE89076 datasets.

Discussion

Given the evidence that m6A contributes to CRC,^19,20 the relationship between candidate SNPs in 20 m6A regulators and the risk of CRC was also investigated.¹⁸ However, a broader analysis method, integrating independent GWAS summary statistics with eQTL data to identify potential functional genetic variants of CRC-associated m6A-SNPs, has more research value.^21,22 Therefore, we used this method combining raw data and summary statistics of genome-wide association studies with eQTL and differential expression data to screen potential candidate CRC-associated m6A-SNPs.

In this study, three CRC-associated m6A-SNPs, specifically rs178184 in NOVA1, rs35782901 in HTR4, and rs60571683 in SLCO1B3, were found to be associated with altered gene expression in CRC. As we known, m6A-SNPs not only affect gene expression level, also involve in disease progression by influencing the ratio between different RNA isoforms and the translation level of their protein products. However, we only explored the potential effect of m6A-SNPs on gene expression level in this study, because other data are currently publicly limited.

The role of NOVA1, HTR4, and SLCO1B3 polymorphisms has not been previously explored in CRC, and genetic variants found near m6A sites have more possibilities to influence the pathogenesis of CRC. We found the m6A-SNP rs178184 is located in the NOVA1 coding gene on chromosome 14, which is essential for growth and invasion-related signaling in cancer cells and is a master regulator of alternative splicing. NOVA1 is a crucial regulator of alternative splicing in pancreatic beta cells,²³ acts as an oncogene in the development of melanoma,²⁴ and regulates hTERT splicing and promotes cell growth in non-small cell lung cancer.²⁵ In CRC, Hong et al showed that NOVA1 is involved in cancer progression, suggesting that NOVA1 might be a valuable prognostic biomarker and a target for CRC treatment. The m6A-SNP rs35782901 is located in the HTR4 coding gene on chromosome 5. Studies have shown that HTR4 variants are associated with chronic obstructive pulmonary disease.²⁶ Moreover, 5-HTR4 was found predominantly in high-grade tumors, with 5-HTR4 inhibition reducing the proliferative activity of androgen-independent prostate cancer cell lines.²⁷ The m6A-SNP rs60571683 is located in the SLCO1B3 coding gene on chromosome 12. SLCO1B3 is a functional transporter that is normally expressed in the liver but that was also detected in different cancers and reported to be involved in cancer.²⁸ For example, SLCO1B3 inhibits tumorigenesis and progression of breast cancer.²⁹ Genotypic variants of SLCO1B3 affect docetaxel pharmacokinetics.³⁰ In addition, the SLCO1B3 GG/AA haplotype is associated with impaired testosterone transport and improved survival in patients with prostate cancer.³¹

Conclusion

This study reports the first comprehensive analysis of GWAS raw data and summary statistics combined with eQTL and differential gene expression data to identify candidate CRC-associated m6A-SNPs. Three m6A-SNPs (rs178184, rs35782901, and rs60571683) were found to be potentially associated with CRC, as demonstrated by their high eQTL signals and altered the expression of their coding genes (NOVA1, HTR4, and SLCO1B3, respectively). However, this study has certain limitations. For example, the identified m6A-SNPs have not been validated in tissue samples. Despite these limitations, the finding will provide the opportunity of further research to elucidate the practical impact of m6A-SNPs on the pathogenesis of CRC.