- 1The First Affiliated Hospital of Chengdu Medical College, School of Clinical Medicine, Chengdu Medical College, Chengdu, China
- 2National Health Commission (NHC), Key Laboratory of Nuclear Technology Medical Transformation, Mianyang Central Hospital, School of Medicine, University of Electronic Science and Technology of China, Mianyang, China
Purpose: To characterize the entire profile of m6A modifications and differential expression patterns for circRNAs in colorectal cancer (CRC).
Methods: First, High-throughput MeRIP-sequencing and RNA-sequencing was used to determine the difference in m6A methylome and expression of circRNA between CRC tissues and tumor-adjacent normal control (NC) tissues. Then, GO and KEGG analysis detected pathways involved in differentially methylated and differentially expressed circRNAs (DEGs). The correlations between m6A status and expression level were calculated using a Pearson correlation analysis. Next, the networks of circRNA-miRNA-mRNA were visualized using the Target Scan and miRanda software. Finally, We describe the relationship of distance between the m6A peak and internal ribosome entry site (IRES) and protein coding potential of circRNAs.
Results: A total of 4340 m6A peaks of circRNAs in CRC tissue and 3216 m6A peaks of circRNAs in NC tissues were detected. A total of 2561 m6A circRNAs in CRC tissues and 2129 m6A circRNAs in NC tissues were detected. Pathway analysis detected that differentially methylated and expressed circRNAs were closely related to cancer. The conjoint analysis of MeRIP-seq and RNA-seq data discovered 30 circRNAs with differentially m6A methylated and synchronously differential expression. RT-qPCR showned circRNAs (has_circ_0032821, has_circ_0019079, has_circ_0093688) were upregulated and circRNAs (hsa_circ_0026782, hsa_circ_0108457) were downregulated in CRC. In the ceRNA network, the 10 hyper-up circRNAs were shown to be associated with 19 miRNAs and regulate 16 mRNAs, 14 hypo-down circRNAs were associated with 30 miRNAs and regulated 27 mRNAs. There was no significant correlation between the level of m6A and the expression of circRNAs. The distance between the m6A peak and IRES was not significantly related to the protein coding potential of circRNAs.
Conclusion: Our study found that there were significant differences in the m6A methylation patterns of circRNAs between CRC and NC tissues. M6A methylation may affect circRNA-miRNA-mRNA co-expression in CRC and further affect the regulation of cancer-related target genes.
Introduction
The latest epidemiological data in 2020 show that colorectal cancer (CRC) is the third most common cancer in the world, with more than 1.93 million new cases, accounting for 9.7% of the world’s newly diagnosed cancers (1). Treatment of CRC includes surgery, radiation therapy, and chemotherapy (2). However, approximately one-quarter of patients have developed liver metastases at the time of first diagnosis, and patients with advanced CRC have a poor 5-year survival and quality of life (3, 4). This highlights the need for a better understanding of the underlying pathogenic mechanisms that promote the development of CRC to identify and treat CRC in the early stage.
N6-Methyladenosine (m6A), methylated adenosine at the N6 position, is the most abundant internal modification of RNAs in eukaryotes (5, 6). Methylation is a reversible epigenetic modification that affects the fate of modified RNAs (7). Methylation involes in RNAs behaviors regulation, such as pre-mRNA splicing, polyadenylation, regulation of RNA stability and long noncoding RNAs biological functions (8–10). In recent years, the function of methylation has been extensively studied in the progression of various cancers, including leukemia, glioma, breast cancer and liver cancer (11–14). M6A modification plays a double role in human carcinogenesis and suppression. Methyltransferases acting as “writers” and demethylases acting as an “erasers” are important for maintaining balanced Methylation. M6A reader proteins specifically recognize m6A transcripts and further regulate gene expression and tumor development (15).
Circular RNAs (circRNAs) characterized by a covalently closed loop produced via back-splicing are a novel class of non-coding RNAs. The connection between the 5’cap and the 3’end make them more stable and not susceptible to degradation (16, 17). CircRNAs can accumulate at high levels in cells and can also be detected in the sera and exosomes (18). Therefore, circRNAs may be a viable candidate as tumor biomarkers. In the past, circRNAs were considered a byproduct but recent evidence suggest that they play a key role in various biological functions, such as acting as microRNA sponges (19, 20), transcriptional regulators (21–23), and translational intermediates (24–26). Multiple studies have shown that circRNAs disorders are involved in cancer progression and tumor chemotherapy resistance (19–23, 27, 28).
M6A modifications are also widespread in both circRNAs and mRNAs. M6A modifications are also widely present in circular RNAs and may affect tumorigenesis and development through various mechanisms. Current research shows that activation of YAP1 by N6-Methyladenosine-Modified circCPSF6 Drives Malignancy in Hepatocellular Carcinoma (29). CircMETTL3, upregulated in an m6A-dependent manner, promotes breast cancer progression. However, the mechanism of cicRna methylation involved in CRC is still in the preliminary exploratory stage (30). Here, we attempted to characterize the patterns of m6A modification in circRNAs from CRC tissues and tumor-adjacent normal control tissues (NC). At the same time, we analyze the relationship between m6A modification and the expression and encoding potential of circRNAs in CRC. We expect that this study can lay a foundation for exploring the oncogenic mechanism of CRC.
Materials & methods
Tissue samples
Five colorectal cancer specimens and adjacent normal tissue were obtained at the Chengdu Medical College from March 2017 to June 2018. None of these patients received radiation or chemotherapy before the specimens were collected. The collected specimens were frozen sectioned and stained with hemAtoxylin-eosin to confirm that the collected specimens were CRC tissues and corresponding normal tissues adjacent to the cancer. This study has been approved by the Ethics Committee of the First Affiliated Hospital of Chengdu Medical College. The clinicopathological data of 5 patients with CRC are the same as the article (doi: 10.3389/fcell.2021.760912) published by our team in Front Cell Dev Biol.
Extraction of total RNA
RNA was extracted using TRIzol Reagent (Invitrogen, Carlsbad, CA, USA) according to the manufacturer’s instructions. RNA concentration was determined using a NanoDrop ND-1000 spectrometer (Thermo, Waltham, MA, USA), and RNA integrity was evaluated by denaturing gel electrophoresis.
MeRIP-seq & RNA-seq
High-throughput RNA sequencing was performed by Cloud-Seq Biotech (Shanghai, China). M6A RNA immunoprecipitation was performed using the GenSeqTM m6A -MeRIP Kit (GenSeq, Beijing, China) according to the manufacturer’s instructions. The input samples without immunoprecipitation and the m6A IP samples were used as templates for the NEBNext® Ultra II Directional RNA Library Prep Kit (New England Biolabs, Inc., USA). These libraries were controlled for quality and quantified using the BioAnalyzer 2100 system (Agilent Technologies, Inc., Palo Alto, CA, USA). Then, library sequencing was performed on an Illumina HiSeq instrument with 150bp paired end reads.
Data analysis
After 3’ adaptor-trimming and the removal of any low-quality reads by the cutadapt software (v1.9.3), the clean reads were aligned to the reference genome using STAR software (v2.5.1b), and the circRNAs were detected and identified using DCC software (v0.4.4). The data were normalized with edgeR(v3.16.5) software and differentially expressed circRNAs with statistical significant expression changes (p ≤ 0.05 and fold change≥2) were identified. MACS (v1.4.2) software was used to identify methylated sites on the RNAs (m6A peaks) and differentially methylated sites were identified using diffReps software (v1.55.3) (P ≤ 0.05 and fold change≥2). Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis were performed by DAVID database.
Validation of circRNAs expression by RT-qPCR
Total RNA was extracted from 8 pairs of CRC and NC tissues using Tatal RNA Extration KIT (Solarbio, Beijing, China) according to the manufacturer’s instructions. The concentration and purity of the RNA was determined using Nanodrop 2000 microspectrophotometer (Thermo Scientific, New York, USA). The reversing transcription was performed according to the PrimeScript RT reagent kit with gDNA Eraser (Takara Biotechnology Co., Ltd, Beijing, China). Real-time quantitative PCR (RT-qPCR) was performed using TB Green™ Premix Ex Taq™ II (Takara Biotechnology Co., Ltd, Beijing, China). All RT-qPCR analyses were conducted in triplicate and the average value was calculated. The expression values of genes were normalized GAPDH and then 2−ΔΔCt transformed for the gene expression.
M6A correlation, ceRNA network, and coding potential prediction
The differences in the m6A fold enrichment for the total circRNAs populations between the CRC and NC groups were used to create a correlation map linking M6A status and expression level. We selected 10 hypermethylated circRNAs with upregulated expression and 14 hypo-methylated circRNAs with downregulated expression to establish the circRNA-miRNA-mRNA network. The regulatory relationships between the circRNA-miRNA-mRNA were visualized using the Target Scan (v8.0) and miRanda (v3.3a) software (31). Then, The correlations between m6A status and expression level in each of the groups were calculated using a Pearson correlation analysis. we used LGC software to predict the coding potential of circular RNA (32).
Results
General characteristics of circRNA m6A modification patterns in CRC and NC tissues
A total of 4340 m6A peaks of circRNAs in CRC tissue and 3216 m6A peaks of circRNAs in NC tissues were detected. A total of 2561 m6A circRNAs in CRC tissues and 2129 m6A circRNAs in NC tissues were detected. There were 2638 overlapping m6A peaks of circRNAs and 1886 overlapping m6A circRNAs between the two groups. Compared to NC tissues, CRC tissues had 1702 unique m6A peaks of circRNAs and 765 m6A circRNAs. These results showed that there was a significant difference in the overall m6A modification pattern between CRC and NC tissues (Figures 1A, B).
Figure 1 General characteristics of circRNA m6A modification patterns. (A) Venn diagram of the number of m6A peaks of circRNAs in CRC and NC groups; (B) Venn diagram of the number of m6A circRNAs in CRC and NC groups; (C) The distribution of the types of total circRNAs, m6A circRNAs and non-m6A circRNAs in CRC group; (D) The distribution of the types of total circRNAs, m6A circRNAs and non-m6A circRNAs in NC group; (E) Distribution of the number of circRNAs (y axis) plotted against the number of m6A peaks in each circRNA (x axis) for CRC and NC group.
CircRNAs can be classified into various types based on the type of sequences they contain. Therefore, we analyzed the distribution of the types of circRNAs represented in the total circRNAs, m6A circRNAs and non-m6A circRNAs groups from both the CRC and NC tissues. This analysis revealed that exonic circRNAs account for 80% of the total circRNAs, 77% of the m6A circRNAs, and 84% of the non-m6A circRNAs in CRC tissues (Figure 1C). Similarly, exonic circRNAs accounted for the highest proportion of transcripts in all three groups from NC tissues at approximately 80% (Figure 1D). This indicated that the exonic circRNAs were the most prevalent version of these circular RNAs, which was consistent with previous studies. In addition, it was found that most circRNAs had the unique m6A modified peak both CRC and NC tissues. A relatively small number of m6A modified genes contain two or more peaks (Figure 1E).
The majority of the total circRNAs from both CRC and NC groups contained two or three exons. The most common of the m6A circRNAs group and the non-m6A circRNAs group are circular RNA transcripts containing two exons (Figures 2A, B). The exon length analysis showed that circRNAs containing more than 6 exons were the longest with an average length of more than 1000bp, followed by single exon circRNAs (Figure 2C). In addition, the exon length of the m6A circRNAs was significantly longer than non-m6A circRNAs (Figure 2D).
Figure 2 Distribution patterns for the exons in the circRNAs. (A) Percentage (y-axis) of exon number (x-axis) of circRNA (left), m6A -circRNA and non-m6A circRNA (right) in CRC group; (B) Percentage (y-axis) of exon number (x-axis) of circRNA (left), m6A -circRNA and non-m6A circRNA (right) in NC group; (C) The distribution of exon lengths (y axis) for each of the input circRNAs (left), m6A circRNAs (middle), and non-m6A circRNAs (right) were plotted against the number of exons (x axis) each circRNA spans. For all box plots, the lower edge of the box represents the first quartile and the upper edge represents the third quartile. The horizontal line inside the box represents the median and the whiskers identify the farthest data points within a 1.5 x interquartile range (IQR); (D) Comparison of exon size (nt) between m6A circRNAs and non-m6A circRNAs.
Distribution of differentially methylated m6A circRNAs
There were 336 differential m6A peaks of circRNAs and 247 differential m6A circRNAs were found in CRC relative to the NC group. Meanwhile, There were 130 hyper-methylated circRNAs and 117 hypo-methylated circRNAs were found in CRC (Figure 3A). The top 10 circRNAs of m6A hypermethylated and hypomethylated with the highest fold-change values are shown in Table 1. The hyper-methylated circRNAs were predominantly located on chromosomes 6, 7 and 10, while the hypo-methylated circRNAs were predominantly located on chromosomes 1, 2, and 12 (Figure 3B). Moreover, The length of hyper-methylated and hypo-methylated m6A circRNAs is mainly enriched in ≤10000 bp (Figure 3C).
Figure 3 General characteristics of differentially methylated m6A circRNAs in CRC. (A) The number of differentially m6A peaks and differentially m6A circRNAs in CRC; (B) The chromosome origins for genes of these differentially methylated circRNAs; (C) The length of the differentially methylated circRNAs; (D) Schematic representation of the exons included in the hsa-circ-0019079 (circKIF20B) and has-circ-0108457 (circSETBP1) circularized transcripts.
We confirmed the back-splicing of hsa_circ_0019079 of hyper-methylated circRNAs and has_circ_0108457 of hypo-methylated circRNAs using CIRI software. Has_circ_0019079 was located on chromosome 10 and spliced by 6 exons from KIF20B. Has_circ_0108457 was located on chromosome 18 and spliced by a single exon from SETBP1 (Figure 3D). Interestingly, both KIF20B and SETBP1 were involved in the development of CRC (33, 34).
Differentially m6A modification circRNAs are involved in important biological pathways
GO and KEGG pathway analyses were performed to evaluate the biological significance of differentially m6A circRNAs. GO analysis shows that hyper-methylated circRNAs were mainly concentrated in cellular macromolecule metabolic processes, nuclear part and nucleic acid binding (Figure 4A). While hypo-methylated circRNAs were mainly enriched in myofibril assembly, contractile fiber and Ras guanyl-nucleotide exchange factor activity (Figure 4B). KEGG pathway analysis revealed that the hyper-methylated circRNAs were mainly involved in DNA replication, RNA transport, ribosome biogenesis in eukaryotes, and protein processing in the endoplasmic reticulum (Figure 4C). The hypo-methylated circRNAs were enriched in pathways of cancers, cGMP-PKG, and tight junction signaling pathways (Figure 4D).
Figure 4 Pathway analysis of differentially m6A circRNAs. (A) GO analysis of hyper-methylated m6A circRNAs; (B) GO analysis of hypo-methylated m6A circRNAs; (C) KEGG analysis of hyper-methylated m6A circRNAs; (D) KEGG analysis of hypo-methylated m6A circRNAs.
Differentially expressed circRNA profiles in CRC
RNA-seq was used to detect differentially expressed circRNAs (fold-change ≥2 and p-value <0.05) in CRC and NC tissues, which showed that there was a significant difference in the expression of circRNAs (Figure 5A). The volcano plot shows that 877 circRNAs were differentially expressed, including 522 downregulated and 355 upregulated circRNAs in CRC (Figure 5B). The top 10 upregulated and downregulated genes are listed in Table 2.
Figure 5 Differentially expressed circRNAs in CRC. (A, B) A heatmap hierarchical clustering and Volcano plot were used to visulize the differentially expressed circRNAs in each of the groups; (C) The RT-qPCR was used to detect the expression of 5 circRNAs in CRC and NC tissues; (D) GO enrichment analysis for Upregulated and Downregulated circRNAs; (E) KEGG pathway analysis for Upregulated and Downregulated circRNAs.
Among these differentially expressed circRNAs, we used RT-qPCR to detect the expression of 5 circRNAs in 8 pairs of CRC and NC tissues. The results of RT-qPCR showed the same expression trend with the RNA-Seq results. Compared with NC tissue, CircRNAs (has-circ-0032821, has-circ-0019079, has-circ-0093688) were upregulated, while circRNAs (hsa_circ_0026782, hsa_circ_0108457) were downregulated in CRC (Figure 5C).
GO analysis revealed that the upregulated circRNAs were involved in the nucleic acid metabolic process, intracellular part and nucleotide binding. The downregulated circRNAs were involved in cytoskeleton organization, intracellular part and nucleotide binding (Figure 5D). KEGG pathway analysis showed that the upregulated circRNAs were closely related to base excision repair, ribosome biogenesis in eukaryotes, and protein processing in the endoplasmic reticulum, while the downregulated circRNAs were closely related to cGMP-PKG signaling pathway and adherens junctions (Figure 5E).
Association between m6A methylation and circRNA expression
Methylation sequencing data showed that the m6A peaks in the circRNAs from CRC tissues was more than that of the NC tissues. Here, we analyzed the methylation correlation between the two groups. The results showed that the m6A level of circRNAs in the CRC was positively correlated with the m6A level of circRNAs in the NC group (Spearman’ s rho = 0.63, p <0.05) (Figure 6A).
Figure 6 Association between m6A methylation and circRNAs expression. (A) Scatter plot showing the linear correlation between circRNAs expression and m6A methylation in CRC and NC group; (B) Venn diagram showing the relationship between m6A modifications and the expression of circRNAs; (C) Scatter plot of the correlation between m6A levels and circRNAs expression in CRC group; (D) Scatter plot of the correlation between m6A levels and circRNAs expression in NC group; (E) Cumulative distribution of circRNAs expression for m6A circRNAs (red) and non-m6A circRNAs (blue).
Combining the results of the methylation sequencing and RNA sequencing, we identified 130 hypermethylated circRNAs, 10 of which were upregulated and 5 were downregulated. A total of 117 circRNAs were hypomethylated, 14 of which were downregulated and one was upregulated (Figure 6B). Interestingly, the methylation level of one of the circRNAs had both up-regulated and down-regulated m6A sites. To analyze the correlation between circRNA methylation and expression levels, we constructed a correlation graph using the fold enrichment of circRNA m6A methylation and expression values described by FPKM. The results indicate that there was a statistically significant positive correlation between methylation and expression of circRNAs in CRC and NC samples (Figures 6C, D).
To further analyze the m6A effects on circRNAs expression, we divided all of the circRNAs into m6A and non- m6A circRNAs groups. We then calculated the log two-fold change values for these circRNAs and generated a cumulative curve. There was no significant difference between m6A and non-m6A circRNAs (Figure 6E).
Co-expression network of circRNA-miRNA-mRNA in CRC
CircRNAs are generally considered as sponges for miRNAs to fine-tune the miRNA-mRNA regulatory network. We selected 10 hyper-up circRNAs and 14 hypo-down circRNAs to establish the circRNA-miRNA-mRNA network. The 10 hyper-up and 14 hypo-down circRNAs are listed in Table 3 and Table 4. In the ceRNA network, the 10 hyper-up circRNAs were shown to be associated with 19 miRNAs and regulate 16 mRNAs (Figure 7A). Similarly, 14 hypo-down circRNAs were associated with 30 miRNAs and regulated 27 mRNAs (Figure 7B). In addition, there was a positive correlation between circRNA and mRNA expression. when the expression of the circRNA was upregulated the expression of the target mRNA was upregulated (Figure 7). In these mRNAs, we found that 40 were associated with tumors and 28 of these were specifically associated with CRC through the GEPIA database. This suggests that m6A may regulate CRC related genes through the circRNA-miRNA-mRNA co-expression network (35–39).
Figure 7 The circRNA-miRNA-mRNA networks in CRC. (A) CeRNA analysis for hyper-methylated upregulated circRNAs; (B) CeRNA analysis for hypo-methylated downregulated circRNAs.
Relationship between circRNAs with coding potential and methylation
Previously, it was thought that circRNAs could not encode any functional products, but recent studies have shown that circRNAs also possess some coding potential with a handful even encoding putative proteins (24–26). In this study, we analyzed the protein-coding potential of circRNAs using LGC software (32). The results showed that 850 of the 7990 (10.64%) circRNAs have a certain protein coding potential. In the 2894 m6A circRNAs, 575 had protein-coding potential, accounting for 19.87%. In the 224 differentially methylated circRNAs, 101 have protein-coding potential, accounting for 45.09%. Of the 30 circRNAs exhibiting both differential methylation and expression, nine were shown to have protein-coding potential, accounting for 30.00% (Figure 8A). Compared with the total circRNA, there are significantly more circRNAs with protein coding potential in the m6A group. As a result, m6A may increase the coding potential of circRNAs.
Figure 8 CircRNAs with coding potential predominantly exhibit N6-adenosine methylation. (A) The percentages of circRNAs with protein-coding potential in total circRNAs, m6A circRNAs, differentially m6A circRNAs, and differentially m6A& differentially expressed circRNAs were compared; The distance between the methylation peak and the ORFs for each of the circRNAs was calculated and then the mean (B) or median (C) value of their coding potential scores were used to stratify the circRNAs into low and high probability groups which were then used to create the density function diagram.
CircRNAs can be translated with an internal ribosome entry site (IRES) (40, 41). And m6A can act as an IRES to drive circRNAs translation in a cap-independent manner with even a single m6A site being sufficient to initiate translation (42). Some non-coding RNAs are reported to have coding potential from 5′ open reading frames (ORFs) (43). This suggests that the ORFs commonly found in circRNAs may be translated by internal m6A sites. Therefore, we evaluated whether the distance between the m6A peak and ORFs could act as a deciding factor for their translation. Based on the mean or median value of their coding potential score, the circRNAs were divided into low and high groups and the density function diagram was drawn (Figures 8B, C). The results showed that the distance of high coding potential score was farther from IRES.
Discussion
In this study, we managed to profile circRNAs that were differentially methylated and expressed in CRC and NC tissues. These results showed that there was a significant difference in m6A methylation and expression of circRNAs between the CRC and NC tissues. There are several reported circRNAs arising from the exonic regions that contribute to cancer progression in CRC (44, 45). Our study revealed that exonic circRNAs account for 80% of the total circRNAs, 77% of the m6A circRNAs, and 84% of the non-m6A circRNAs in CRC tissues. The m6A methylation levels in the CRC group were significantly higher than NC group. There was a significant positive correlation between methylation and expression levels. In addition, m6A modification can affect the expression of circRNAs altering their fine-tuning of the miRNA-mRNA regulatory axis, and thus affecting the expression of tumor-related target mRNAs. We found that methylation affected the coding potential of the circRNAs but this effect was not related to the distance between the m6A peak and ORFs.
CircRNAs occur in the nucleus, and most circRNAs containing introns are confined to the nucleus but most of the exon circRNAs are localized to the cytoplasm (16). The m6A binding protein YTHDC1 regulates the export of methylated mRNA and can also regulate the transfer of m6A circRNAs from the nucleus to the cytoplasm (46). M6A modified circNSUN2 is exported to the cytoplasm by YTHDC1 and functions to improve the stability of HMGA2 mRNA and promote liver metastasis in CRC (44). Our results show that the methylation level of circRNAs in CRC tissues is significantly increased and that this methylation is mainly concentrated in circRNAs composed of long exons. some circRNAs arising from the exonic regions contribute to the development and progression of this cancer (44, 45). The exon circRNAs with methylation may be worthy of further study in colorectal cancer.
Circ_0032821 was significantly upregulated in human GC tumors and cells. its’ expression induced cell proliferation, EMT, migration, invasion and autophagy inhibition in human GC cells through activating MEK1/ERK1/2 signaling pathway (46). Circ_0032821 was also highly expressed in OXA-resistant GC cells and contributes to oxaliplatin (OXA) resistance of gastric cancer cells by regulating SOX9 via miR-515-5p (47). Hsa_circ_0026782 (circITGA7) and its linear host gene ITGA7 are both significantly downregulated in CRC tissues and cell lines. These decreased expression levels correlated with CRC progression. They play a suppressor in CRC. CircITGA7 inhibits the proliferation and metastasis of CRC cells by suppressing the Ras signalling pathway and promoting the transcription of ITGA7 (22). In another study found that circITGA7 sponges miR-3187-3p to upregulate ASXL1, suppressing colorectal cancer proliferation (48). ITGA7 also plays an important tumorigenic function and acts as a suppress gene in breast cancer (49). In this study, we found that circ_0032821 were upregulated and hsa_circ_0026782 were downregulated in CRC. The expression trend of those circRNAs is consistent. In addition, We verified that hsa_circ_0019079 (circKIF20B) were upregulated and has_circ_0108457 (circSETBP1) were downregulated in CRC tissues. There is no research on hsa-circ-0019079 and has-circ-0108457. However, both of KIF20B and SETBP1 are involved in the development of CRC (33, 34). Their role in colorectal cancer is worthy of further study.
By cross analyzing the m6A-Seq and RNA-seq data, we identified 30 circRNAs with differential methylation and expression. 33.3% (10/30) of the host genes (KIF20B, BRCA1, NAV1, SIPA1L2, FMN2, MMP28, SETBP1, EIF4E3, MICU3, and RUNX1T1) were associated with CRC. The expressions of NAV1, SIPA1L2, FMN2, MMP28, SETBP1, EIF4E3, MICU3, and RUNX1T1 have shown to be decreased in CRC samples and reduced expression levels of SIPA1L2, MMP28, and EIF4E3 are associated with lower overall survival rates in CRC patients. KIF20B is upregulated in CRC promoting the migration and invasion of CRC (33). BRCA1 is associated with CRC and breast cancer (50, 51). MMP28 is involved in the occurrence and metastasis of gastric cancer and CRC (36, 52). Therefore, those circRNAs with differential methylation and expression patterns deserve further study, which may help to clarify the molecular function underlying the occurrence and development of CRC.
When we established a circRNA-miRNA-mRNA network by integrating matched circRNAs, miRNAs, and mRNAs expression profiles we were able to demonstrate a positive correlation between the expression of circRNAs and mRNA. We found that 30 differentially methylated circRNAs were associated with 40 tumor-related mRNAs and 28 of these were associated with CRC. For example, circRNAs hsa_circ_0106593 and hsa_circ_000292 serve as sponges for hsa-miR-5000-5p and hsa-miR-6842-3p, respectively, regulating the expression of the target gene AJUBA. AJUBA promotes the growth of CRC by inhibiting apoptosis (30). circRNA hsa_circ_008487 binds to hsa-miR-608 and hsa-miR-6812-5p, which bind to CCL21, which further inhibits the migration and invasion of CRC cells (31). circRNAs hsa_circ_0019079 and chr1:35558995-35578782+ positively regulate the expression of SLC31A1, hsa_circ_0114420, and chr12:10572963-10588009- which can positively regulate the expression of MMP28. The expression of ZFPM2 decreases with decreasing hsa_circ_0112394 expression and SLC31A1, MMP28, and ZFPM2 are associated with the occurrence, progression, and prognosis of CRC (35, 36, 39). M6A is a post-transcriptional modification which regulates circRNAs expression. We speculate that m6A could be involved in the occurrence and development of CRC by affecting the circRNA-miRNA-mRNA co-expression networks. This may provide a theory for the mechanism underlying circRNA activity in CRC. Furthermore, regulating m6A modifications may be a future strategy for the treatment of CRC.
CircRNAs are defined as non-coding RNAs. however, it has recently been shown that some circRNAs actually encode proteins, and that these proteins are involved in the occurrence, development and drug resistance of many tumors (22–24). The 5’ and 3’ untranslated regions (UTRs) are essential elements for canonical cap-dependent translation in eukaryotic cells. Due to the lack of the 5 ‘cap and 3’ end, the translation of circRNA can only be initiated in a cap-independent manner. IRES- and m6A- mediated cap-independent translation initiation are accepted as important mechanisms for circRNA translation (40–42). IRESs are sequences located in the 5’ UTR of mRNAs that directly recruit ribosomes to initiate translation (53). M6A in the 5’ UTR can directly bind to eukaryotic initiation factor 3 (eIF3) prompting translation (54). M6A driven translation is an alternative mechanism often employed under stress conditions (55). Interfering with the m6A modification level in RNAs can influence translation efficiency (42, 56) and a single m6A site is sufficient to initiate circRNA translation via eIF4G2 and m6A reader YTHDF3 (42). In this study, we predicted the protein-coding potential of all the circRNAs and found that methylation increases the coding potential of circRNA transcripts. There were significantly more circRNAs with protein-coding potential in the m6A circRNA group and the proportion (45.09%) of potentially coding transcripts was highest in the differentially methylated circRNAs. This may be because these short RNA elements containing m6A sites have IRES-like activity initiating the translation of the circRNAs (42). We further analyzed if the distance between the methylation site and ORFs may predict greater protein coding potential. However, further analysis of the density function diagram indicated that methylation affected the coding potential of the circRNAs, but this effect was not related to the distance between the methylation peak and the ORFs.
Our study found that there were significant differences in the m6A methylation patterns of circRNAs between CRC and NC tissues and that the level of m6A methylation could affect the expression level and increase the coding potential of circRNAs. Bioinformatic analyses showed that the host genes of differentially methylated circRNAs were associated with tumor development-related pathways. In addition, we showed that m6A methylation may affect circRNA-miRNA-mRNA co-expression in CRC and further affect the regulation of cancer-related target genes. However, the differentially methylated circRNAs contribute to the occurrence and development of CRC needs further study. In our follow-up research, I will target the circRNAs selected in this experiment, such as hsa_circ_0093688, hsa_circ_0019079, hsa_circ_0108457, hsa_circ_0043136. Methyltransferase METTL3/METTL14, demethylase FTO/ALKBH5, methylation reader YTHDF/IGF2BP2 will be used to verify the methylation expression level of circRNAs, and functional exploration will be carried out based on KEGG and GO results. The network map of circRNA-miRNA-mRNA was used to find the regulated target genes for mechanism research.
Conclusions
This study found that there were significant differences in the m6A methylation patterns of circRNAs between CRC and NC tissues. M6A methylation may affect circRNA-miRNA-mRNA co-expression in CRC and further affect the regulation of cancer-related target genes.
Data availability statement
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE190388.
Ethics statement
The studies involving human participants were reviewed and approved by The First Affiliated Hospital of Chengdu Medical College. The patients/participants provided their written informed consent to participate in this study. Written informed consent was obtained from the individual(s) for the publication of any potentially identifiable images or data included in this article.
Author contributions
YZ and FH were involved in the conception and design of the study. QG and GXJ acquired the data and drafted the manuscript. QG collected the samples and conducted data analysis. YZ critically revised the manuscript for intellectual content, approved the final version. All authors read and approved the final manuscript.
Funding
The present study was supported by the scientific research project of Sichuan Provincial Health Commission (grant No. 19PJ190), Mianyang Central Hospital hospital-level project (grant No. 2021YJRC-002).
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.
References
1. Bray F, Ferlay J, Soerjomataram I, Siegel RL, Torre LA, Jemal A, et al. Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA CancerJ Clin (2018) 68:394–424. doi: 10.3322/caac.21492
2. Kopetz S, Chang GJ, Overman MJ, Eng C, Sargent DJ, Larson DW, et al. Improved survival in metastatic colorectal cancer is associated with adoption of hepatic resection and improved chemotherapy. J Clin Oncol (2009) 27(22):3677–83. doi: 10.1200/JCO.2008.20.5278
3. Arnold M, Sierra MS, Laversanne M, Soerjomataram I, Jemal A, Bray F, et al. Global patterns and trends in colorectal cancer incidence and mortality. Gut (2017) 66(4):683–91. doi: 10.1136/gutjnl-2015-310912
4. Watanabe T, Kobunai T, Yamamoto Y, Matsuda K, Ishihara S, Nozawa K, et al. Gene expression signature and response to the use of leucovorin, fluorouracil and oxaliplatin in colorectal cancer patients. Clin Transl Oncol (2011) 13(6):419–25. doi: 10.1007/s12094-011-0676-z
5. Zhou C, Molinie B, Daneshvar K, Pondick JV, Wang J, Van Wittenberghe N, et al. Genome-wide maps of m6A circRNAs identify widespread and cell-type-specific methylation patterns that are distinct from mRNAs. Cell Rep (2017) 20(9):2262–76. doi: 10.1016/j.celrep.2017.08.027
6. Zheng Y, Nie P, Peng D, He Z, Liu M, Xie Y, et al. m6AVar: A database of functional variants involved in m6A modification. Nucleic Acids Res (2018) 46(D1):D139–45. doi: 10.1093/nar/gkx895
7. Wu B, Li L, Huang Y, Ma J, Min J. Readers, writers and erasers of N(6)-methylated adenosine modification. Curr Opin Struct Biol (2017) 47:67–76. doi: 10.1016/j.sbi.2017.05.011
8. Wang X, Lu Z, Gomez A, Hon GC, Yue Y, Han D, et al. N6-methyladenosine-dependent regulation of messenger RNA stability. Nature (2014) 505(7481):117–20. doi: 10.1038/nature12730
9. Zhou J, Wan J, Gao X, Zhang X, Jaffrey SR, Qian SB, et al. Dynamic m(6)A mRNA methylation directs translational control of heat shock response. Nature (2015) 526(7574):591–4. doi: 10.1038/nature15377
10. Xiao W, Adhikari S, Dahal U, Chen YS, Hao YJ, Sun BF, et al. Nuclear m(6)A reader YTHDC1 regulates mRNA splicing. Mol Cell (2016) 61(4):507–19. doi: 10.1016/j.molcel.2016.01.012
11. Feng M, Xie X, Han G, Zhang T, Li Y, Li Y, et al. YBX1 is required for maintaining myeloid leukemia cell survival by regulating BCL2 stability in an m6A-dependent manner. Blood (2021) 138(1):71–85. doi: 10.1182/blood.2020009676
12. Du J, Ji H, Ma S, Jin J, Mi S, Hou K, et al. m6A regulator-mediated methylation modification patterns and characteristics of immunity and stemness in low-grade glioma. Brief Bioinform (2021) 22(5):bbab013. doi: 10.1093/bib/bbab013
13. Einstein JM, Perelis M, Chaim IA, Meena JK, Nussbacher JK, Tankka AT, et al. Inhibition of YTHDF2 triggers proteotoxic cell death in MYC-driven breast cancer. Mol Cell (2021) 81(15):3048–64.e9. doi: 10.1016/j.molcel.2021.06.014
14. Zhang C, Huang S, Zhuang H, Ruan S, Zhou Z, Huang K, et al. YTHDF2 promotes the liver cancer stem cell phenotype and cancer metastasis by regulating OCT4 expression via m6A RNA methylation. Oncogene (2020) 39(23):4507–18. doi: 10.1038/s41388-020-1303-7
15. Jiang X, Liu B, Nie Z, Duan L, Xiong Q, Jin Z, et al. The role of m6A modification in the biological functions and diseases. Signal Transduct Target Ther (2021) 6(1):74. doi: 10.1038/s41392-020-00450-x
16. Li Z, Huang C, Bao C, Chen L, Lin M, Wang X, et al. Exon-intron circular RNAs regulate transcription in the nucleus. Nat Struct Mol Biol (2015) 22(3):256–64. doi: 10.1038/nsmb.2959
17. Qu S, Zhong Y, Shang R, Zhang X, Song W, Kjems J, et al. The emerging landscape of circular RNA in life processes. RNA Biol (2017) 14(8):992–9. doi: 10.1080/15476286.2016.1220473
18. Qu S, Zhong Y, Shang R, Zhang X, Song W, Kjems J, et al. Circular RNA is enriched and stable in exosomes: a promising biomarker for cancer diagnosis. Cell Res (2015) 25(8):981–4. doi: 10.1038/cr.2015.82
19. Zhang X, Wang S, Wang H, Cao J, Huang X, Chen Z, et al. Circular RNA circNRIP1 acts as a microRNA-149-5p sponge to promote gastric cancer progression via the AKT1/mTOR pathway. Mol Cancer (2019) 18(1):20. doi: 10.1186/s12943-018-0935-5
20. Han D, Li J, Wang H, Su X, Hou J, Gu Y, et al. Circular RNA circMTO1 acts as the sponge of microRNA-9 to suppress hepatocellular carcinoma progression. Hepatology (2017) 66(4):1151–64. doi: 10.1002/hep.29270
21. Abdelmohsen K, Panda AC, Munk R, Grammatikakis I, Dudekula DB, De S, et al. Identification of HuR target circular RNAs uncovers suppression of PABPN1 translation by CircPABPN1. RNA Biol (2017) 14(3):361–9. doi: 10.1080/15476286.2017.1279788
22. Li X, Wang J, Zhang C, Lin C, Zhang J, Zhang W, et al. Circular RNA circITGA7 inhibits colorectal cancer growth and metastasis by modulating the ras pathway and upregulating transcription of its host gene ITGA7. J Pathol (2018) 246(2):166–79. doi: 10.1002/path.5125
23. Chen L, Yang X, Zhao J, Xiong M, Almaraihah R, Chen Z, et al. Circ_0008532 promotes bladder cancer progression by regulation of the miR-155-5p/miR-330-5p/MTGR1 axis. J Exp Clin Cancer Res (2020) 39(1):94. doi: 10.1186/s13046-020-01592-0
24. Chen CY, Sarnow P. Initiation of protein synthesis by the eukaryotic translational apparatus on circular RNAs. Science (1995) 268(5209):415–7. doi: 10.1126/science.7536344
25. Pan Z, Cai J, Lin J, Zhou H, Peng J, Liang J, et al. A novel protein encoded by circFNDC3B inhibits tumor progression and EMT through regulating snail in colon cancer. Mol Cancer (2020) 19(1):71. doi: 10.1186/s12943-020-01179-5
26. Zheng X, Chen L, Zhou Y, Wang Q, Zheng Z, Xu B, et al. A novel protein encoded by a circular RNA circPPP1R12A promotes tumor pathogenesis and metastasis of colon cancer via hippo-YAP signaling. Mol Cancer (2019) 18(1):47. doi: 10.1186/s12943-019-1010-6
27. Huang X, Li Z, Zhang Q, Wang W, Li B, Wang L, et al. Circular RNA AKT3 upregulates PIK3R1 to enhance cisplatin resistance in gastric cancer via miR-198 suppression. Mol Cancer (2019) 18(1):71. doi: 10.1186/s12943-019-0969-3
28. Hon KW, Ab-Mutalib NS, Abdullah NMA, Jamal R, Abu N. Extracellular vesicle-derived circular RNAs confers chemoresistance in colorectal cancer. Sci Rep (2019) 9(1):16497. doi: 10.1038/s41598-019-53063-y
29. Chen Y, Ling Z, Cai X, Xu Y, Lv Z, Man D, et al. Activation of YAP1 by N6-Methyladenosine-Modified circCPSF6 drives malignancy in hepatocellular carcinoma. Cancer Res (2022) 82(4):599–614. doi: 10.1158/0008-5472.CAN-21-1628
30. Li Z, Yang HY, Dai XY, Zhang X, Huang YZ, Shi L, et al. CircMETTL3, upregulated in a m6A-dependent manner, promotes breast cancer progression. Int J Biol Sci (2021) 17(5):1178–90. doi: 10.7150/ijbs.57783
31. Smoot ME, Ono K, Ruscheinski J, Wang PL, Ideker T. Cytoscape 2.8: new features for data integration and network visualization. Bioinformatics (2011) 27(3):431–2. doi: 10.1093/bioinformatics/bfq675
32. Wang G, Yin H, Li B, Yu C, Wang F, Xu X, et al. Characterization and identification of long non-coding RNAs based on feature relationship. Bioinformatics (2019) 35(17):2949–56. doi: 10.1093/bioinformatics/btz008
33. Lin WF, Lin XL, Fu SW, Yang L, Tang CT, Gao YJ, et al. Pseudopod-associated protein KIF20B promotes Gli1-induced epithelial-mesenchymal transition modulated by pseudopodial actin dynamic in human colorectal cancer. Mol Carcinog (2018) 57(7):911–25. doi: 10.1002/mc.22812
34. Torrejón B, Cristóbal I, Caramés C, Prieto-Potín I, Chamizo C, Santos A, et al. Analysis of potential alterations affecting SETBP1 as a novel contributing mechanism to inhibit PP2A in colorectal cancer patients. World J Surg (2018) 42(11):3771–8. doi: 10.1007/s00268-018-4684-9
35. Jia H, Song L, Cong Q, Wang J, Xu H, Chu Y, et al. The LIM protein AJUBA promotes colorectal cancer cell survival through suppression of JAK1/STAT1/IFIT2 network. Oncogene (2017) 36(19):2655–66. doi: 10.1038/onc.2016.418
36. Bister VO, Salmela MT, Karjalainen-Lindsberg ML, Uria J, Lohi J, Puolakkainen P, et al. Differential expression of three matrix metalloproteinases, MMP-19, MMP-26, and MMP-28, in normal and inflamed intestine and colon cancer. Dig Dis Sci (2004) 49(4):653–61. doi: 10.1023/B:DDAS.0000026314.12474.17
37. Zhang Z, He T, Huang L, Ouyang Y, Li J, Huang Y, et al. Two precision medicine predictive tools for six malignant solid tumors: from gene-based research to clinical application. J Transl Med (2019) 17(1):405. doi: 10.1186/s12967-019-02151-8
38. Rong Y, Chen X, Fan D, Lin X, Zheng X, Zhou C, et al. Influence of CCL21 on the invasion and metastasis of colorectal cancer. Zhonghua Wei Chang Wai Ke Za Zhi (2017) 20(11):1300–5.
39. Barresi V, Trovato-Salinaro A, Spampinato G, Musso N, Castorina S, Rizzarelli E, et al. Transcriptome analysis of copper homeostasis genes reveals coordinated upregulation of SLC31A1,SCO1, and COX11 in colorectal cancer. FEBS Open Bio (2016) 6(8):794–806. doi: 10.1002/2211-5463.12060
40. Pamudurti NR, Bartok O, Jens M, Ashwal-Fluss R, Stottmeister C, Ruhe L, et al. Translation of CircRNAs. Mol Cell (2017) 66(1):9–21.e7. doi: 10.1016/j.molcel.2017.02.021
41. Yang Y, Wang Z. IRES-mediated cap-independent translation, a path leading to hidden proteome. J Mol Cell Biol (2019) 11(10):911–9. doi: 10.1093/jmcb/mjz091
42. Yang Y, Fan X, Mao M, Song X, Wu P, Zhang Y, et al. Extensive translation of circular RNAs driven by N6-methyladenosine. Cell Res (2017) 27(5):626–41. doi: 10.1038/cr.2017.31
43. Ingolia NT, Brar GA, Stern-Ginossar N, Harris MS, Talhouarne GJ, Jackson SE, et al. Ribosome profiling reveals pervasive translation outside of annotated protein-coding genes. Cell Rep (2014) 8(5):1365–79. doi: 10.1016/j.celrep.2014.07.045
44. Chen RX, Chen X, Xia LP, Zhang JX, Pan ZZ, Ma XD, et al. N6-methyladenosine modification of circNSUN2 facilitates cytoplasmic export and stabilizes HMGA2 to promote colorectal liver metastasis. Nat Commun (2019) 10(1):4695. doi: 10.1038/s41467-019-12651-2
45. Hsiao KY, Lin YC, Gupta SK, Chang N, Yen L, Sun HS, et al. Noncoding effects of circular RNA CCDC66 promote colon cancer growth and metastasis. Cancer Res (2017) 77(9):2339–50. doi: 10.1158/0008-5472.CAN-16-1883
46. Roundtree IA, Luo GZ, Zhang Z, Wang X, Zhou T, Cui Y, et al. YTHDC1 mediates nuclear export of N6-methyladenosine methylated mRNAs. Elife (2017) 6:e31311. doi: 10.7554/eLife.31311
47. Jiang Y, Zhang Y, Chu F, Xu L, Wu H. Circ_0032821 acts as an oncogene in cell proliferation, metastasis and autophagy in human gastric cancer cells in vitro and in vivo through activating MEK1/ERK1/2 signaling pathway. Cancer Cell Int (2020) 20:74. doi: 10.1186/s12935-020-1151-0
48. Yang G, Zhang T, Ye J, Yang J, Chen C, Cai S, et al. Circ-ITGA7 sponges miR-3187-3p to upregulate ASXL1, suppressing colorectal cancer proliferation. Cancer Manag Res (2019) 11:6499–509. doi: 10.2147/CMAR.S203137
49. Bhandari A, Xia E, Zhou Y, Guan Y, Xiang J, Kong L, et al. ITGA7 functions as a tumor suppressor and regulates migration and invasion in breast cancer. Cancer Manag Res (2018) 10:969–76. doi: 10.2147/CMAR.S160379
50. Garcia JM, Rodriguez R, Dominguez G, Silva JM, Provencio M, Silva J, et al. Prognostic significance of the allelic loss of the BRCA1 gene in colorectal cancer. Gut (2003) 52(12):1756–63. doi: 10.1136/gut.52.12.1756
51. Waks AG, Cohen O, Kochupurakkal B, Kim D, Dunn CE, Buendia Buendia J, et al. Reversion and non-reversion mechanisms of resistance to PARP inhibitor or platinum chemotherapy in BRCA1/2-mutant metastatic breast cancer. Ann Oncol (2020) 31(5):590–8. doi: 10.1016/j.annonc.2020.02.008
52. Li Y, Sun Q, Jiang M, Li S, Zhang J, Xu Z, et al. KLF9 suppresses gastric cancer cell invasion and metastasis through transcriptional inhibition of MMP28. FASEB J (2019) 33(7):7915–28. doi: 10.1096/fj.201802531R
53. Godet AC, David F, Hantelys F, Tatin F, Lacazette E, Garmy-Susini B, et al. IRES trans-acting factors, key actors of the stress response. Int J Mol Sci (2019) 20(4):924. doi: 10.3390/ijms20040924
54. Meyer KD, Patil DP, Zhou J, Zinoviev A, Skabkin MA, Elemento O, et al. 5' UTR m(6)A promotes cap-independent translation. Cell (2015) 163(4):999–1010. doi: 10.1016/j.cell.2015.10.012
55. Zhou J, Wan J, Shu XE, Mao Y, Liu XM, Yuan X, et al. N6-methyladenosine guides mRNA alternative translation during integrated stress response. Mol Cell (2018) 69(4):636–47.e7. doi: 10.1016/j.molcel.2018.01.019
Keywords: colorectal cancer, MeRIP sequencing, M6A modification, CircRNA, RNA-Seq
Citation: He F, Guo Q, Jiang G-x and Zhou Y (2022) Comprehensive analysis of m6A circRNAs identified in colorectal cancer by MeRIP sequencing. Front. Oncol. 12:927810. doi: 10.3389/fonc.2022.927810
Received: 25 April 2022; Accepted: 01 August 2022;
Published: 22 August 2022.
Edited by:
Norfilza M Mokhtar, National University of Malaysia, MalaysiaReviewed by:
Dawei Chen, University of Liège, BelgiumJustin Roberts, University of South Alabama, United States
Copyright © 2022 He, Guo, Jiang and Zhou. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Yan Zhou, enFsdnp5MzE5QDE2My5jb20=
†These authors have contributed equally to this work