CircCCDC134 was significantly upregulated in CC tissues
circRNA-Seq was used to analyse the dysregulated circRNA expression between 4 pericarcinomatous and 4 CC tissues, which was supported by GENESEED (Guangzhou, China). We found hundreds of circRNAs that varied between the CC tissue and adjacent normal tissue. We detected 21 upregulated circRNAs with a 2-fold change (p < 0.05) (Fig. 1A and B), and one of the new molecules that attracted our attention was circCCDC134 (hsa_circ_0008806, 580 bp). Next, the characteristics of circCCDC134 in CC were analysed. Primers for circCCDC134 were designed for the qRT–PCR experiments (Supplementary Fig. 1 A and B). The circular structure of circCCDC134 was confirmed by amplification with divergent primers, followed by Sanger sequencing, which showed that circCCDC134 was generated from exon 2 to exon 6 of the CCDC134 gene through back splicing (Fig. 1C). Oligo (dT) and random primers were employed to identify whether circCCDC134 contains poly-A tails by qRT–PCR. If the expression of the target gene cannot be detected in the reverse transcribed cDNA of Oligo (dT), the RNA has no poly-A tail and is a circular RNA. We verified the circular feature of circCCDC134 and conducted qRT–PCR using oligo (dT) and random primers, and the results suggested that compared with random primers, circCCDC134 expression was noticeably reduced, while CCDC134 expression did not change when using oligo (dT) primers, revealing that circCCDC134 contained no poly-A tail (Fig. 1D). To further detect the features of circCCDC134, both circCCDC134 and the linear mRNA CCDC134 were treated with RNase R and the transcription inhibitor actinomycin D. Ribonuclease R is only able to degrade linear RNA but does not affect circRNA. The results showed that circCCDC134 was more resistant to both RNase R (Fig. 1E and F) and actinomycin D (Fig. 1G), suggesting that circCCDC134 was more stable than linear CCDC134. For the localization of circRNA, FISH and RNA fractionation assays were further performed. The cytoplasmic and nuclear RNA of ME180 and SiHA cell lines were separated and extracted, and qRT–PCR was performed to determine circCCDC134 expression in the cytoplasm and nucleus by RNA fractionation assays. The fluorescence signal (Fig. 1H) and the expression of circCCDC134 (Fig. 1I and Supplementary Fig. 1C) were found in both the nucleus and cytoplasm in the ME180 and SiHA cell lines. qRT–PCR assays were used to explore the expression of circCCDC134 in CC tissues and cells. The qRT–PCR results revealed circCCDC134 overexpression in the tumour and metastatic tissues (Fig. 1J) and CC cells. Among the CC cells, the expression of circCCDC134 was higher in the ME180 and SiHA cells and lower in the HCC94 cells than in the other CC cells (Fig. 1K).
m6A modification of circCCDC134 enhanced its stability
First, we explored the potential mechanism of circCCDC134 upregulation in CC. Previous research has shown that m6A modification plays a pivotal role in posttranscriptional regulation and the biogenesis of circRNAs. We analysed the m6A modification of circCCDC134 using circPrimer (Fig. 2A) and the SRAMP prediction server (http://www.cuilab.cn/sramp/) (Fig. 2B), and there were many m6A modification sites of circCCDC134. The motif analysis of the circCCDC134 methylation site based on the very high confidence of SRAMP is shown in Fig. 2C.
To explore this further, we first performed RNA pull-down assays and a mass spectrometry analysis (ChIRP-MS) to screen circCCDC134-interacting proteins and comprehensively identify RNA-binding proteins. The results of ChIRP-MS (Supplementary Fig. 1D and E) and silver staining (Fig. 2D) showed that circCCDC134 cooperated with multiple proteins, including ALKBH5 and YTHDF2. Previous studies demonstrated that a-ketoglutarate-dependent dioxygenase AlkB homologue 5 (ALKBH5) (a demethylase) can remove m6A methylation from its target RNAs and lead to lower levels of m6A [18]; the reader protein YTH N6-methyladenosine RNA binding protein 2 (YTHDF2) can selectively bind m6A-modified RNAs, recruit them to mRNA decay sites and then control target RNA stability [19]. To gain insight into the functions of circCCDC134-binding proteins, we implemented Protein–Protein Interaction Networks (PPI) and Gene Ontology (GO) based on mass spectrometry data. The PPI (Supplementary Fig. 1 F) and GO (Fig. 2E) results revealed that ALKBH5 and YTHDF2 mediated m6A modification, which could affect RNA stability. Therefore, we propose that circCCDC134 stability in CC is tuned by ALKBH5-mediated m6A modification in a YTHDF2-dependent manner, increasing its expression. To further identify the ALKBH5 protein and YTHDF2 protein interacting with circCCDC134, we performed circRNA pull-down and RIP assays. The results of the circCCDC134 pull-down assay showed that the ALKBH5 protein and YTHDF2 protein could bind circCCDC134-MS2 (Fig. 2F and G). We performed RIP assays using anti-ALKBH5 or anti-YTHDF2 antibodies, and the results demonstrated that both ALKBH5 and YTHDF2 could interact with circCCDC134 (Fig. 2H and I). In order to quantify m6A level of circCCDC134, m6A RNA immunoprecipitation (MeRIP) was performed. Gene-specific m6A qPCR to detect the m6A methylation status of circCCDC134 was performed and the results indicated that an reduction of m6A methylation in the ALKBH5 overexpression group of SiHA (Supplementary 1 G) and ME180 (Supplementary 1 H) cell lines. Next, we further confirmed the effect of ALKBH5 expression on circCCDC134 expression. We further analysed the correlation between ALKBH5 and circCCDC134 expression in CC. The qPCR and IHC results showed that the expression of ALKBH5 was low in CC (Fig. 2J and K), and patients with high ALKBH5 expression had much longer overall survival (Supplementary Fig. 1I). Next, circCCDC134 expression was significantly negatively correlated with ALKBH5 expression in CC (Fig. 2L) based on the qPCR results. Moreover, the qRT–PCR and actinomycin D assays showed that the expression (Fig. 2M) and stability (Fig. 2N) of circCCDC134 were decreased with the overexpression of ALKBH5. Taken together, our data suggest that the ALKBH5-mediated m6A methylation of circCCDC134 is responsible for circCCDC134 upregulation by increasing its RNA stability via YTHDF2.
circCCDC134 facilitates CC cell proliferation and metastasis in vitro and in vivo
To gain insight into the function of circCCDC134, we characterized the oncogenic phenotypes in ME180 and SiHA cells with circCCDC134 silencing (Si-circCCDC134-1 and Si-circCCDC134-2) and HCC94 cells with circCCDC134 overexpression (EX-circCCDC134). The qPCR results showed that Si-circCCDC134-1 and Si-circCCDC134-2 could significantly inhibit circCCDC134 expression in CC cells (ME180 and SiHA cells) (Fig. 3A), and we chose Si-circCCDC134-1 and Si-circCCDC134-2 for the cell function experiments. We investigated the role of circCCDC134 in the proliferation of CC cells by a colony formation assay. To assess the influence of circCCDC134 on cell migration and invasion, Transwell assays were used to detect the cell migration and invasion capacity. The colony formation and Transwell assays revealed that the knockdown of circCCDC134 impaired the proliferation (Fig. 3B), migration (Fig. 3D), and invasion (Fig. 3E) of CC cells. Additionally, the overexpression of circCCDC134 significantly promoted CC cell proliferation (Fig. 3C), migration and invasion (Fig. 3F). To further evaluate the effect of circCCDC134 on CC tumour growth and metastasis in vivo, we established xenograft growth and lung metastasis models. Notably, the xenografts (Fig. 3G) formed by ME180 cells bearing sh-circCCDC134 exhibited significantly slower growth rates (Fig. 3H) and lower tumour weights (Fig. 3I) than the scrambled controls. Furthermore, the volumes of the metastatic nodes in the lung in the circCCDC134 knockdown group were apparently smaller than those in the control group (Fig. 3J), and the number of metastatic nodes was significantly decreased in the circCCDC134 silencing group (Fig. 3K). Taken together, these findings suggest that circCCDC134 acts as an oncogene in CC metastasis by promoting cellular proliferation, migration and invasion.
circCCDC134 dysregulation contributes to the abnormal expression of genes in CC
To clarify the abnormal expression of genes related to CC metastasis and circCCDC134 disorder, mRNA-seq of CC (Fig. 4A), mRNA-seq of metastasis tissues (Fig. 4B) and mRNA-seq of si-circCCDC134/ME180 cells (ME180 cells transfected with si-circCCDC134-1) (Fig. 4C and D) were carried out. Furthermore, we implemented Kyoto Encyclopedia of Genes and Genomes (KEGG) (Fig. 4E) and gene set enrichment analysis (GSEA) (Fig. 4F and G) based on the si-circCCDC134 mRNA-seq data. The results revealed that the upregulated circCCDC134 expression was closely associated with the TNF signalling pathway. It has been acknowledged that the TNF signalling pathway [20] and epithelial-mesenchymal transition (EMT) [21] contribute to the metastatic progression of CC. To identify potential TNF and EMT mRNAs regulated by circCCDC134, we analysed the mRNA-seq data of CC tissues, metastatic tissues and the si-circCCDC134 group. Regarding TNF-related genes, by combining the TNF signalling pathway in the si-circCCDC134, CC and metastasis mRNA-seq data results, 5 genes were found, including SIK1, IER2, EFNA1, KLF2 and PER1 (Fig. 4H and J). Regarding EMT-related genes, 200 EMT-related genes were downloaded from the Molecular Signature database v7.1 (http://www.broad.mit.edu/gsea/msigdb/) [22]. By combining the EMT-related gene data, si-circCCDC134 mRNA-seq data, CC mRNA-seq data and metastasis mRNA-seq data, 4 genes were found, including DCN, PCOLCE, PMP22 and COL5A3 (Fig. 4I and J). Next, we explored the expression and survival analysis of these genes in CC based on TCGA databases, and the analysis revealed that EFNA1, KLF2, PER1 and COL5A3 expression was closely associated with CC overall survival (Supplementary Fig. 2). In addition, the qRT–PCR analysis confirmed that the expression of EFNA1 was decreased in ME180 and SiHA cells and increased in HCC94 cells that were induced by circCCDC134. Regarding KLF2, PER1 and COL5A3, their expression was increased in ME180 and SiHA cells and decreased in HCC94 cells induced by circCCDC134 (Fig. 4K and L). These results demonstrate that circCCDC134 dysregulation contributes to the abnormal expression of TNF-related mRNAs and EMT-related mRNAs.
circCCDC134 stimulates HIF1A transcription by recruiting p65 and acting as a miR-503-5p sponge
Many studies have demonstrated that circRNAs play regulatory roles by binding proteins and sponging miRNAs. In the above experiments, ChIRP-MS was performed to determine whether circCCDC134 could interact with binding proteins, and the results showed that circCCDC134 cooperated with multiple nuclear proteins. Among the binding proteins, we found that circCCDC134 could interact with p65, which is a very important transcription factor that promotes CC metastasis [23] and interacts with the TNF signalling pathway [24] and EMT [25]. We explored the possibility that circCCDC134 interacts with p65, and the secondary structure of circCCDC134 was assembled in Mfold (version 2.3) following the methods described by Du William W et al. [26]. The molecular simulation results supported that circCCDC134 could perfectly dock with p65 (Fig. 5A). The results of the circRNA pull-down (Fig. 5B) and RIP (Fig. 5C) assays demonstrated that p65 could interact with circCCDC134. The western blot results showed that there was no significant change in the p65 protein levels upon circCCDC134 knockdown, indicating that circCCDC134 does not regulate the expression of the P65 protein (Fig. 5D). Furthermore, ChIP-seq and ChIP–qPCR assays were used to identify target genes regulated by p65. As shown in Fig. 5E and F, p65 binds the upstream 5 k region of 7.34% of genes that contain the gene promoter region based on the ChIP-seq results. Furthermore, AnimalTFDB 3.0, which is a comprehensive resource for the annotation and prediction of animal transcription factors, was used to analyse the key target genes of p65 [27]. By combining AnimalTFDB 3.0 data of p65, si-circCCDC134 mRNA-seq data, CC mRNA-seq data and P65 ChIP-seq data, HIF1A was found (Fig. 5G and H). We tested whether circCCDC134 regulates HIF1A transcription by enhancing p65, and the results were as follows: 1. p65 was enriched in the region − 1200 to − 400 bp from the transcription start site of the HIF1A promoter (Fig. 5I); 2. p65 enrichment at − 1200 to − 400 bp from the transcription start site of the HIF1A promoter was reduced after interference with p65 and circCCDC134 expression (Fig. 5J and K). Moreover, the qPCR results showed that the expression of HIF1A was rescued by transfection with p65 and si-circCCDC134 (Fig. 5L). Immunohistochemistry showed that the expression of HIF1A was positively correlated with the expression of circCCDC134 in the xenograft growth and lung metastasis models (Fig. 5M and N). These results suggest that circCCDC134 interacts with p65 to enhance HIF1A transcription.
Another key role of circRNAs is to sponge miRNAs and regulate downstream genes. To explore the miRNAs that could be sponged by circCCDC134, miRNA-seq data from CC tissues and metastatic tissues (Fig. 6A and B), the circBank database (http://www.circbank.cn/) and the circRNA database (https://circinteractome.nia.nih.gov/) were used. Among these data, only hsa-miR-503-5p was low in CC tissues and may be sponged by circCCDC134 (Fig. 6C). The qRT–PCR results revealed that the expression of miR-503-5p was low in the tumour and metastatic tissues (Fig. 6D) and CC cells (Fig. 6E). The circRNA pull-down assays showed that miR-503-5p could specifically bind circCCDC134 (Fig. 6F). The RNA FISH assays revealed that circCCDC134 and miR-503-5p colocalized in the cytoplasm (Fig. 6G). The RIP assay also illustrated that circCCDC134 and miR-503-5p enrichment was increased in the Ago2 group compared that in the IgG group (Fig. 6H). We investigated the role of miR-503-5p in the proliferation of CC cells by a colony formation assay. To assess the influence of miR-503-5p on cell migration and invasion, Transwell assays were used to detect the cell migration and invasion capacity. The results showed that the transfection of the miR-503-5p mimic resulted in a decrease in proliferation, migration and invasion ability in ME180 cells, and the function of miR-503-5p was rescued by the reintroduction of circCCDC134 (Supplementary Fig. 3A). Furthermore, we sought to identify the specific target gene of miR-503-5p via the miRDB, miRTarBase and TargetScan databases. According to the miRDB, miRTarBase and TargetScan database analysis, 90 target genes were found (Supplementary Fig. 3B). Combined with the si-circCCDC134 mRNA-seq data, 14 target genes were analysed (Supplementary Fig. 3C).
Next, we explored the OS (overall survival) and RFS (relapse-free survival) associated with these 14 genes in CC based on TCGA databases. The analysis revealed that MYB (Supplementary Fig. 3D), CD2AP, KPNA3, NUFIP2, and TLL1 (Supplementary Figs. 4 and 5) expression was closely related to CC OS and RFS. Since MYB is known to be associated with CC migration and invasion [26], we chose MYB for further study. The luciferase reporter gene assay consists of cloning both the wild-type and mutated forms of the 3’UTR of the miRNA predicted mRNA target downstream of a firefly luciferase reporter. The luciferase reporter gene assay showed that luciferase activity was significantly inhibited in HEK-293 T cells that were cotransfected with the miR-503-5p mimic and wild-type luciferase reporter (circCCDC134-WT and MYB-WT) compared with those in the control group (Fig. 6I). Moreover, the qPCR results showed that the expression of MYB was rescued by transfection with sh-circCCDC134 and anti-miR-503-5p (Supplementary Fig. 3E) or transfection with circCCDC134 and miR-503-5p mimics. These findings indicate that circCCDC134 can act as a sponge for miR-503-5p to upregulate MYB in CC cells.
HIF1A as a downstream target of MYB
In previous studies, MYB was shown to be a very important transcription factor. To explore the downstream target of MYB, the JASPAR database (https://jaspar.genereg.net/) was used. Interestingly, HIF1A was an important MYB target gene (Fig. 7A). The motif analysis of the HIF1A transcription promoter site based on the high confidence of JASPAR is shown in Fig. 7B. The ChIP results showed that 1. MYB was enriched in the region − 1200 to − 200 bp from the transcription start site of the HIF1A promoter (Fig. 7C), and 2. MYB enrichment in − 1200 to − 200 bp from the transcription start site of the HIF1A promoter was reduced after MYB interference (Fig. 7D). Moreover, the qPCR results showed that the expression of HIF1A was rescued by transfection with MYB (Fig. 7E). These results demonstrate that MYB could stimulate HIF1A transcription.