Bioinformatics analyses revealed the involvement of m6A modification in cervical cancer progression
To better understand whether and how m6A regulators contribute to cervical cancer progression, we first identified 9 m6A writers (WTAP, ZC3H13, METTL3, METTL14, METTL16, VIRMA, RBM15B, RBM15, and CBLL1), 15 m6A readers (FMR1, hnRNPA2B1, hnRNPC, YTHDF1/2/3, YTHDC1/2, LRPPRC, IGF2BP1, IGF2BP2, IGF2BP3, RBMX, EIF3A, and ELAVL1), and 2 m6A erasers (FTO and ALKBH5) that were previously reported to have pathogenic roles in other human cancer types [7, 8]. Using the TCGA CESC database, we identified copy number variations (CNVs) and somatic mutations in the m6A-related genes because amplification or deletions that affect copy number can also affect expression of the corresponding gene. This analysis showed that 21 m6A regulators harbored CNVs in cervical cancer (Additional file 2: Fig. S1). Among these genes, 13 of the 21 carried mutations at a frequency of 1–4%, while the remaining 8 genes had mutations at a frequency of < 1% (Additional file 2: Fig. S2a). Analysis of mutation co-occurrence indicated that several of these genes also carried mutations that significantly co-occurred with mutations in other m6A regulators (Additional file 2: Fig. S2b). RNA-seq data in the TCGA database further showed that these m6A regulators were differentially expressed in cervical cancer samples compared with the expression in adjacent normal tissues (Additional file 2: Fig. S2c), which suggested that CNVs and other somatic mutations could lead to dysregulation of the m6A-related genes in cervical cancer.
Continuing our analysis of possible relationships between m6A-related genes and cervical cancer in public datasets, we next focused on the potential clinical significance of these genes. Survival analyses indicated that 15 CNV-carrying genes had altered levels of expression that were significantly correlated with cervical cancer patient prognosis (Additional file 2: Fig. S3), whereas 6 of the genes showed no significant value as candidate prognostic indicators. Subsequent gene interaction network analysis indicated that 3 genes (ZC3H13, METTL14, and CBLL1) appeared to serve as network hubs, suggesting that dysregulation of m6A modification by these genes provided major contributions to the development and/or progression of cervical cancer (Additional file 2: Fig. S4a). Correlation analyses also revealed a strong association among these 21 m6A regulators (Additional file 2: Fig. S4b). Notably, unsupervised clustering of the 21 m6A-related gene expression patterns in different patients identified three distinct patterns of m6A modification (i.e., clusters A-C) associated with cervical cancer (Additional file 2: Fig. S4c, d). Survival analyses showed that the prognosis of patients with cluster B pattern generally had a better prognosis than patients exhibiting cluster A or C patterns (Additional file 2: Fig. S4e).
To determine the biological significance of these m6A modification patterns, we conducted GSVA to elucidate differentially-enriched KEGG pathways and GO terms in these clusters that are related to the m6A modification patterns identified in cervical cancer samples. The results showed that clusters A and C were enriched for canonical cancer signaling pathways and processes, such as the VEGF, mTOR, ERBB, MAPK, Wnt, TGF-β, hedgehog, and Notch pathways, as well as tight junctions, adherens junctions, cell cycle, non-homologous end joining, mitotic sister chromatid cohesion, mRNA processing, and related GO terms compared with samples exhibiting cluster B patterns (Additional file 2: Fig. S5), further indicating the involvement of m6A modification in cervical cancer progression.
CENPK expression is correlated with aberrant m6A modification and tumorigenic gene expression in cervical cancer
We compared the differential expression data from these 21 genes with the m6A methylation data available in the m6A-Atlas database and the data obtained in our previous study of cervical cancer patients [23] to identify specific genes that could be dysregulated by m6A modification in cervical cancer. First, we performed differential expression analyses between m6A cluster B and C patterns, which identified 3628 differentially-expressed genes in different m6A clusters [fold-change ≥ 2.0 (geneset A); P < 0.05]. Furthermore, experimental evidence in the m6A-Atlas database (http://180.208.58.66/m6A-Atlas/) suggested the participation of 7617 genes with an m6A modification in HeLa cells (geneset B). Our previous study [23] also showed a possible correlation between 151 differentially-expressed genes and cervical cancer cell progression [fold-change ≥ 2.0 (geneset C); P < 0.05]. A total of 40 overlapping genes were obtained from the intersection of the 3 genesets (Additional file 2: Fig. S6a).
Subsequent GSVA revealed that pathways involving CENPK had the greatest overlap with m6A-modified clusters associated with cervical cancer. Specifically, CENPK was enriched for genes in pathways, such as Wnt signaling, DNA damage repair signaling, the cell cycle, and DNA replication (Additional file 2: Fig. S6b). CENPK expression was elevated in pan-cancer datasets compared to that in corresponding normal tissues based on TCGA (Additional file 2: Fig. S6c). Interestingly, GSEA showed that CENPK also participated in the modulation of platinum drug resistance more prominently than radioresistance (Additional file 2: Fig. S6d). In addition, correlation analyses showed a positive association between CENPK and stemness markers, including EPCAM, CD133, SOX2, and OCT4 (Additional file 2: Fig. S6e). Taken together, these results indicated a strong correlation between CENPK expression and the dysregulation of m6A modifications in cervical cancer.
ZC3H13-mediated m6A methylation upregulated CENPK mRNA to activate pro-tumorigenic functions
We further examined whether CENPK expression was regulated by m6A RNA methylation based on the results of our bioinformatics analyses. The above mentioned survival analysis showed 9 m6A regulators (EIF3A, hnRNPC, LRPPRC, RBM15B, VIRMA, WTAP, YTHDF2, YTHDF3, and ZC3H13) with the potential to confer poor patient survival (Additional file 2: Fig. S3). Indeed, CENPK expression was correlated with CBLL1, ELAVL1, FMR1, METTL3, METTL14, RBM15, WTAP, YTHDC1, YTHDC2, and ZC3H13 expression (Fig. 1a, b). Thus, ZC3H13 and WTAP were correlated with poor patient survival and CENPK expression, and of these two candidates, CENPK had a more prominent correlation with ZC3H13 expression than with WTAP expression. Because CENPK has a potential oncogenic role in cervical cancer, we speculated that ZC3H13 might be responsible for m6A methylation of CENPK mRNA. To test this possibility, we generated HeLa and SiHa cells with ZC3H13 knockdown by siRNA (si-ZC3H13), which downregulated CENPK expression (Fig. 1c). In addition, MeRIP assays confirmed that CENPK mRNA was enriched by anti-m6A antibody in HeLa and SiHa cells, while m6A modification of CENPK mRNA was reduced in si-ZC3H13 cells (Fig. 1d). Moreover, DART-seq data in the RMVar database (http://rmvar.renlab.org) identified a site in the 3’ UTR of CENPK mRNA on chr5:65518395(-) as an m6A modification site (Fig. 1e). Luciferase reporter assays suggested that even though ZC3H13 suppression impaired the transcription of wild-type CENPK, non-synonymous cytosine-to-adenosine conversion mutations at this m6A site abolished the downregulation (Fig. 1f). These results thus indicated that CENPK expression was modulated by ZC3H13-associated m6A modification.
Based on these findings, we further investigated the effects of ZC3H13 on CENPK-mediated downstream signaling. Considering its overlap with cervical cancer-enriched pathways (vide supra), we focused on Wnt and p53 signaling. TOP/FOP luciferase reporter and dual-luciferase assays revealed that ZC3H13 knockdown decreased Wnt signaling activity and elevated p53 signaling activity in HeLa and SiHa cells, but this effect was abolished by overexpressing CENPK (ov-CENPK) (Additional file 2: Fig. S7a, b). Indeed, tumorsphere and colony formation, and Transwell, EdU, and immunofluorescence staining assays showed that ZC3H13 knockdown had a suppressive effect on cervical cancer stemness, chemoresistance, metastasis, and cell proliferation that was reversed by ov-CENPK in HeLa and SiHa cells (Additional file 2: Fig. S7c-i). These results collectively demonstrated that ZC3H13 functioned as a regulator of CENPK expression through m6A RNA methylation and that these genes functioned together in facilitating cervical cancer progression.
CENPK expression was correlated with cervical cancer pathology
To further establish the central role of CENPK in cervical cancer, we conducted IHC analysis of 119 cervical cancer samples and 35 adjacent normal tissues. The results showed elevation of CENPK expression in cervical cancer compared with adjacent normal tissues (Fig. 2a). Moreover, IHC staining of CENPK protein indicated a positive association with Ki67 protein levels (Fig. 2b) and was positively correlated with cancer recurrence (Additional file 1: Table S3).
Survival analysis suggested that cervical cancer patients with high CENPK expression exhibited poor overall survival (Fig. 2c). Because CENPK expression was correlated with cancer recurrence, we further investigated the relationship between CENPK expression and recurrence-free survival in cervical cancer patients. Survival analysis of cervical cancer patients confirmed that high CENPK expression was a predictor of poor recurrence-free survival, which is consistent with correlation analysis (Fig. 2d). Quantitative analysis of IHC staining, followed by univariate and multivariate COX hazard analyses, indicated that CENPK expression was an independent and unfavorable prognostic indicator for overall and recurrence-free survival in cervical cancer patients (Fig. 2e, f). A nomogram showing the role of CENPK, age, and T classification in predicting cancer recurrence in cervical cancer patients is shown in Fig. 2g. When stratified by age or T classification, CENPK expression remained correlated with overall and recurrence-free survival (Additional file 2: Fig. S8a, b). Consistent with these results, correlation analyses suggested a positive association between β-catenin, Ki67, and CENPK expression (Additional file 2: Fig. S8c). Moreover, ZC3H13 expression was positively associated with EPCAM, CD133, β-catenin, Ki67, N-cadherin, and CCND1 expression in the TCGA database (Additional file 2: Fig. S8d). Taken together, these data showed a clear role for CENPK in cervical cancer development and indicated the clinical value as a potential biomarker.
Silencing CENPK suppressed cervical cancer stemness, chemoresistance, metastasis, and proliferation
Given this validation of CENPK contribution to cervical cancer, we sought to determine how CENPK affected cell stemness, chemoresistance, metastasis, and proliferation. To this end, CENPK was silenced by shRNA (sh-CENPK) or siRNA (si-CENPK) in HeLa and SiHa cell lines. Tumorsphere formation assays and immunofluorescence staining of cervical cancer stem cell markers (CD133 and CD44) indicated that CENPK knockdown reduced stemness of HeLa and SiHa cells (Fig. 3a, b). Clonogenic assays revealed that CENPK suppression impaired cisplatin and carboplatin resistance (Fig. 3c), Transwell assays showed impaired migration and invasion (Fig. 3d), and MTT and colony formation assays revealed inhibited cell growth by sh-CENPK or si-CENPK in HeLa and SiHa cells (Fig. 3e, f). Moreover, immunofluorescence staining of γ-H2AX (Ser139) showed that downregulation of CENPK led to enhanced DNA damage in cervical cancer cells treated with platinum-based drugs (Fig. 3g). Consistent with the bioinformatics results, EdU assays suggested that CENPK knockdown decreased the percentage of cells in the S phase and inhibited DNA replication in HeLa and SiHa cells (Fig. 3h). In addition, the expression of DNA damage repair-associated proteins (p53 and p21) was enhanced, while cell stemness, epithelial-mesenchymal transition (EMT), and DNA replication-associated proteins (c-Myc, Vimentin, CCND1, and c-Jun) were suppressed in si-CENPK cells (Fig. 3i; Additional file 2: Fig. S9). These results indicated that high CENPK expression stimulated cervical cancer stemness, chemoresistance, metastasis, and proliferation.
We established cervical cancer models in BALB/c-nu mice harboring wild-type or stably-silenced CENPK to determine how CENPK affected cancer progression in vivo. We confirmed that CENPK expression was decreased in tumors treated with CENPK-targeted shRNA (Fig. 4a). Mice with downregulated CENPK in HeLa and SiHa cells exhibited decreased tumor formation, fewer lung metastases, and slower growth rates than the corresponding controls (Fig. 4b-f). IHC staining also showed lower Ki67 expression in sh-CENPK tumors compared to scramble shRNA control tumors (Fig. 4g). Moreover, CENPK silencing prolonged the overall survival time of model mice, and CENPK suppression led to synergistic effects with chemotherapy to further extend the mouse survival time compared with that of mice receiving chemotherapy only (Fig. 4h).
CENPK activated Wnt signaling and inactivated p53 signaling via SOX6 in cervical cancer
We next determined the specific mechanisms by which CENPK exerted oncogenic effects. A review of the BioGRID database (https://downloads.thebiogrid.org/BioGRID) suggested that SOX6 potentially interacted with CENPK, and our above bioinformatics analysis indicated that CENPK modulated Wnt signaling (Additional file 2: Fig. S6). Together with TOP/FOP-flash assays, which confirmed CENPK promoted Wnt signaling (Fig. 5a), and previous studies that showed SOX6 activated p53 signaling and inactivated Wnt signaling [27] while also suppressing cervical cancer progression [28], we selected SOX6 as a candidate protein interaction partner of CENPK. Co-IP assays verified that CENPK interacted with SOX6 (Fig. 5b), and immunofluorescent staining confirmed the co-localization of CENPK and SOX6 in HeLa and SiHa cells (Fig. 5c).
SOX6 is known to inhibit the transcription of target genes regulated by β-catenin via direct interaction with β-catenin to inactivate Wnt signaling [27]. Moreover, TOP/FOP-flash assays suggested that SOX6 knockdown had regulatory functions in CENPK-mediated Wnt signaling (Fig. 5a). These results thus implied that CENPK-SOX6 interaction might affect interactions between SOX6 and β-catenin. Co-IP analyses subsequently revealed an interaction between SOX6 and β-catenin, and this SOX6 interaction with β-catenin was enhanced in si-CENPK cells with decreased CENPK-SOX6 binding (Fig. 5d). Moreover, immunofluorescence staining and cell fractionation assays suggested that CENPK knockdown inhibited β-catenin expression and nuclear translocation; however, no impact on SOX6 localization or expression by CENPK was detected (Fig. 5e).
We next investigated the effects of CENPK on p53 signaling in light of the report that SOX6 activates p53 signaling by increasing protein stability via transcriptional suppression of c-Myc in HeLa cells [28]. Dual-luciferase assays showed increased p53 signaling in CENPK knockdown cells, whereas SOX6 silencing abolished this effect (Fig. 5f). ChIP assays indicated that SOX6 binding to the c-Myc transcription regulatory region was increased in HeLa and SiHa cells (Fig. 5g). Furthermore, cycloheximide chase assays, Co-IP analyses, and cell fractionation assays all confirmed that CENPK knockdown led to increased p53 stability and inhibition of p53 ubiquitination and nuclear export, all of which were reversed by SOX6 knockdown (Fig. 5h-j). These data suggested that SOX6 mediated the effects of CENPK on Wnt and p53 signaling.
CENPK promoted tumorigenic functions of cervical cancer cells via Wnt and p53 signaling
We then investigated how CENPK function in cervical cancer was mediated by Wnt and p53 signaling. Immunofluorescence staining, tumorsphere formation, clonogenic, MTT, and EdU assays collectively showed that β-catenin overexpression or p53 knockdown reversed the inhibitory effects of CENPK knockdown on cell stemness, chemoresistance, migration, invasion, and proliferation in HeLa and SiHa cells (Additional file 2: Fig. S10a-f). Immunofluorescence staining and Western blotting further showed that β-catenin overexpression or p53 knockdown abolished the diverse effects of CENPK knockdown on DNA damage repair, cell stemness, EMT, and DNA replication-associated gene expression [i.e., γ-H2AX (Ser139), p53, p21, c-Myc, Vimentin, CCND1, and c-Jun] (Additional file 2: Fig. S10g, h). These results thus showed that the pro-tumorigenic effects of CENPK were mediated by Wnt and p53 downstream regulatory targets.