Shuoling Chen1,2†
Zunnan Huang- 1Key Laboratory of Big Data Mining and Precision Drug Design of Guangdong Medical University, Research Platform Service Management Center, Guangdong Medical University, Dongguan, China
- 2The Second School of Clinical Medicine, Guangdong Medical University, Dongguan, China
- 3Key Laboratory for Research and Development of Natural Drugs of Guangdong Province, School of Pharmacy, Guangdong Medical University, Dongguan, China
- 4Institute of Marine Biomedical Research, Guangdong Medical University, Zhanjiang, China
Background: miRNAs and genes can serve as biomarkers for the prognosis and therapy of cervical tumors whose metastasis into lymph nodes is closely associated with disease progression and poor prognosis.
Methods: R software and Bioconductor packages were employed to identify differentially expressed miRNAs (DEMs) from The Cancer Genome Atlas (TCGA) database. GEO2R detected differentially expressed genes (DEGs) in the GSE7410 dataset originating from the Gene Expression Omnibus (GEO). A Cox proportional hazard regression model was established to select prognostic miRNA biomarkers. Online tools such as TargetScan and miRDB predicted target genes, and overlapping DEGs and target genes were defined as consensus genes. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment and Gene Ontology (GO) function annotations were performed to discern the potential functions of consensus genes. STRING and Cytoscape screened key genes and constructed a regulatory network.
Results: A combination of four miRNAs (down-regulated miR-502 and miR-145, up-regulated miR-142 and miR-33b) was identified as an independent prognostic signature of cervical cancer. A total of 94 consensus genes were significantly enriched in 7 KEGG pathways and 19 GO function annotations including the cAMP signaling pathway, the plasma membrane, integral components of the plasma membrane, cell adhesion, etc. The module analysis suggested that CXCL12, IGF1, PTPRC CDH5, RAD51B, REV3L, and WDHD1 are key genes that significantly correlate with cervical cancer lymph node metastasis.
Conclusions: This study demonstrates that a four-miRNA signature can be a prognostic biomarker, and seven key genes are significantly associated with lymph node metastasis in cervical cancer patients. These miRNAs and key genes have the potential to be therapeutic targets for cervical cancer. Among them, two miRNAs (miR-502 and miR-33b) and two key genes (PTPRC and CDH5) were first reported to be potential novel biomarkers for cervical cancer. The current study further characterizes the progression of lymph node metastasis and mechanism of cervical tumors; therefore, it provides a novel diagnostic indicator and therapeutic targets for future clinical treatments.
Introduction
Worldwide, cervical cancer is the fourth most common female malignancy with an extremely high morbidity and mortality rate (Cohen et al., 2019). According to the World Health Organization (WHO), there were 569,874 new cases of cervical tumors in 2018, accounting for 3.2% of all new cancer cases in that same year (Bray et al., 2018). Most cervical tumors are attributed to recurrent HPV infections. Aberrant expression of oncoproteins encoded by HPV genetic material, such as E6 and E7, partially leads to epigenetic instability, which affects the carcinogenesis and metastasis of cervical cancer (Doorbar, 2006).
MicroRNAs (miRNAs) are associated with the development of a wide range of cancers, including cervical cancer (Park et al., 2017). miRNAs are a group of non-coding, single-stranded RNA molecules that are approximately 22 nucleotides (nt) in length and encoded by endogenous genes (Bartel, 2004; Chen and Kang, 2015). Abnormally expressed miRNAs regulate various biological processes such as apoptosis, proliferation, and metabolism (Zhang et al., 2007; Lee and Dutta, 2009). Moreover, miRNA dysregulation plays a critical role in cancer progression and metastasis in multiple cancers, including cervical cancer (Gómez-Gómez et al., 2013; Huang et al., 2017). For example, Yao et al. (Yao et al., 2018) demonstrated that decreased expression of HPGD by miR-146b-3p induced proliferation, migration, and anchorage-independent growth of cervical cancer cells by activating the STAT3 and AKT signaling pathways. This indicated that miRNAs may be clinically applicable as potential biomarkers and therapeutic targets.
Tumor metastasis is the leading cause of death for most cervical cancer patients. It generally involves complex processes such as extracellular matrix degradation, lymphangiogenesis, angiogenesis, clonal growth at secondary sites, etc. (Zhang et al., 2015). The metastasis and invasion of tumor cells through blood and lymph nodes are crucial processes in the progression of cervical cancer (Dai et al., 2017). Newly formed lymphatic vessels are comprised of endothelial cells that are not tightly connected. Therefore, cervical cancer cells can easily invade this endothelial layer and metastasize to lymph nodes (Zhang et al., 2015). In addition, the poor prognosis of cervical cancer patients is correlated with high invasiveness and diffuse lymph node metastasis (Wang et al., 2016). Lymph node metastasis is the main metastatic pathway and the most critical factor in the prognosis and recurrence of cervical cancer cases (Huang and Fang, 2018). However, a detailed understanding of the specific signatures and molecular mechanisms of lymph node metastasis is lacking. Therefore, the search for a reliable lymphatic signature is important in determining malignant metastasis, and it also provides information that can guide the clinical treatment of cervical tumor patients.
Many have proposed conducting a large-scale systematic analysis of miRNAs, genes, and clinical data using bioinformatics to further characterize the functions of miRNAs and genes in certain disease states, clarify their potential as disease-related signatures, and discover new disease biomarkers and drug targets (Liang et al., 2017). Researchers can obtain tumor data from public databases and conduct differential expression, survival, and prognosis analyses as well as target gene prediction and functional characterization using R language (Sepulveda, 2020), TargetScan (Agarwal et al., 2015), miRDB (Wong and Wang, 2015), DAVID (Huang et al., 2009), etc. With these tools, tumor biomarkers can be screened, and their mechanism of action can be further elucidated. For example, Liang et al. (Liang et al., 2017) identified a three-miRNA signature (miR-145, miR-200c, and miR-218-1) that is a prognostic factor of cervical tumors by conducting Cox univariate and multivariate analyses on differentially expressed miRNAs screened from clinical samples in TCGA database.
In this study, we conducted a multi-step analysis using various R language packages on clinical samples downloaded from the TCGA (Tomczak et al., 2015) and GEO (Barrett et al., 2007) databases to identify DEMs and DEGs. A Cox proportional hazard regression model was then established to determine potential prognostic biomarkers from the available DEMs. Subsequently, the target genes of the miRNAs biomarkers were predicted by the TargetScan and miRDB online tools and the consensus genes were further determined based on overlap between DEGs and these target genes. Finally, MCODE (Bader and Hogue, 2003) software in Cytoscape (Shannon et al., 2003) was used to identify key genes related to lymph node metastasis caused by early-stage cervical tumors. Together, the prognostic miRNAs and key genes mined in this study will provide new insights in elucidating the molecular mechanisms of cervical cancer and contribute to finding new therapeutic targets and prognostic biomarkers for cervical cancer patients.
Materials and Methods
Data Preparation and Differential Expression Analysis
The expression profiles and clinical data for 310 samples were obtained from TCGA (https://portal.gdc.cancer.gov/) database on July 27, 2018. These samples included 307 primary solid cervical cancer tissue samples and three normal tissue samples. Even though this study used many more cancer samples than normal samples, previous studies by Yang et al. (2019), Zeng et al. (2018), etc. have demonstrated that appropriate prognostic models with diagnostic indicators and clinically significant therapeutic targets can be obtained with similar unbalanced sample sets. DEMs were analyzed using the edgeR (Robinson et al., 2010), gplots (Li S. et al., 2020), and limma (Ritchie et al., 2015) R language packages (Version 3.5.1) according to the screening criteria of |log2FC| > 1 and Padj < 0.05. The GSE7410 expression profile data from 24 samples were obtained from the GEO (https://www.ncbi.nlm.nih.gov/geo/) database. The samples included 19 early-stage cervical tumor tissues with lymph node metastasis and five healthy cervical tissues. DEGs were analyzed using GEO2R based on the screening criteria of |log2FC| > 1 and Padj < 0.05.
Establishment of a Cox Proportional Hazard Regression Model
The Cox proportional hazard regression model (George et al., 2014) was established to analyze the association between DEMs and overall survival. The survival package was used in the univariate and multivariate Cox analyses of DEMs (Liang et al., 2017). miRNAs with P < 0.05 as calculated by the univariate Cox analysis were considered strongly correlated with overall patient survival. The multivariate Cox analysis used stepwise regression to screen a prognostic model based on the Akaike information criterion (AIC) value (Vrieze, 2012), where the model with the smallest AIC value contains the smallest number of miRNAs that best predict cervical cancer patient prognosis. miRNAs with P < 0.05 as calculated by the multivariate Cox analysis were considered as independent prognostic factors.
Prognostic Model Construction
Using the results of the multivariate Cox analysis, the risk score was calculated as follows: Risk Score = Exp (miRNA1) × β1 + Exp (miRNA2) × β2 +… + Exp (miRNAn) × βn. The patients were divided into low- and high-risk groups based on their median risk score. A risk score curve was plotted to demonstrate the risk score differences according to the model. A survival status map was plotted to demonstrate the survival status of each sample. A heatmap was plotted to demonstrate the expression level differences of the four prognostic miRNAs in the low- and high-risk groups. A survival curve was plotted to demonstrate the 5-year survival in high- and low-risk groups. And an ROC curve of the model was constructed to determine its predictive ability. A model with an AUC value larger than 0.7 possesses a strong prediction function.
Target Gene Prediction
TargetScan and miRDB predicted the target genes of miRNAs obtained from the multivariate Cox analysis. Target genes that were predicted by both tools were considered target genes of those miRNAs. And the consensus genes were further obtained from the overlap between DEGs related to lymph node metastasis and candidate target genes from the prognosis miRNAs.
KEGG Pathways and GO Function Annotation Enrichment Analysis
To further understand the underlying biological significance of DEMs and DEGs, DAVID (https://david.ncifcrf.gov/) produced KEGG pathway (Kanehisa and Goto, 2000) enrichment and GO functional annotations (Ashburner et al., 2000) of the consensus genes using P < 0.05 as a demarcation criterion. The pathways or annotations with the smallest P-value or the largest count value were considered as key KEGG pathways and GO function annotations. GO function annotations include three parts: biological process (BP), cellular component (CC), and molecular function (MF).
PPI Network Construction and Module Analysis
STRING (http://string-db.org) (Szklarczyk et al., 2017) analyzed the consensus genes and obtained protein interaction data. Proteins with a minimum required interaction score greater than or equal to 0.400 were selected to construct the PPI (protein-protein interaction) network, and nodes with network interruption were hidden. The combined score was imported into Cytoscape software (Version 3.7.1, https://cytoscape.org/). The Molecular Complex Detection (MCODE) plug-in in Cytoscape calculated the MCODE score and selected the significant modules of key genes using the screening criteria of Degree Cut-off = 2, Haircut on, Node Score Cut-off = 0.2, k-core = 2, and Max. Depth = 100. Moreover, the logFC value of genes in the interaction network was also imported into Cytoscape to show the up/down regulation status.
Network Visualization
Regulatory relationship data were imported into Cytoscape along with miRNAs, key genes, consensus genes, key KEGG pathways, and GO functional annotations. The visualization network was performed with Cytoscape to explore the potential regulatory relationship among them.
Analysis Procedure
Figure 1 shows the analysis procedure of the data mining processes that was used to screen tumor biomarkers and key genes in this study. It was based on the extensive use of R language and various online analytical tools. We first obtained DEMs and DEGs using R language and GEO2R, respectively, by analyzing cervical cancer-related miRNAs and gene expression profiles downloaded from TCGA and GEO databases. A prognostic model was constructed via Cox analysis to detect key miRNAs to form a best-available prognostic signature whose target genes were predicted by TargetScan and miRDB. Subsequently, we performed KEGG pathway and GO annotation analyses using DAVID to clarify the function of the consensus genes from the overlap between DEGs and the predicted target genes of key miRNAs. Furthermore, STRING and MCODE selected key genes among consensus genes. Finally, we constructed a miRNAs-Genes-Pathways and Annotations network to elucidate potential regulatory mechanisms. The approach described here is a feasible protocol to identify potential tumor biomarkers and key genes.
Figure 1 Analysis procedure of the data mining process used to screen tumor biomarkers and key genes in this study. It includes specific bioinformatics methods, data processing tools, and partial research results.
Results
miRNA and Gene Differential Expression Analyses
A total of 110 DEMs were obtained after analyzing miRNA expression profiles from TCGA with R language using Padj < 0.05 and |log2FC| > 1 as screening criteria. Among them, 64 miRNAs were significantly down-regulated, and 46 miRNAs were significantly up-regulated (Table S1). At the same time, 1840 DEGs related to early cervical cancer lymph node metastasis were identified by analyzing the GSE7410 expression profile with GEO2R using Padj < 0.05 and |log2FC| > 1 as screening criteria. Among them, 1,298 genes were significantly down-regulated, and 542 genes were significantly up-regulated (Table S2). The volcano map illustrates the significant differences and distribution of the fold change in DEMs and DEGs (Figure 2).
Figure 2 Volcano plot of DEMs (A) and DEGs (B). The abscissa represents the log2 transformation value of the differential expression fold change between the cervical cancer samples and the normal samples. The larger the | logFC | value is, the greater the fold change is. The ordinate represents the -log10 transformation value of the FDR value. The larger the -log10 transformation value is, the more significant the difference is. Green dots represent significantly down-regulated miRNAs or genes. Red dots represent significantly up-regulated miRNAs or genes.
miRNA-Based Signature Identification as a Prognostic Biomarker
A total of 15 miRNAs related to patient survival were obtained from 110 DEMs using univariate Cox analysis (P < 0.05) (Table 1). Four miRNAs related to patient prognosis (miR-502, miR-142, miR-145, and miR-33b) were further screened from the above 15 miRNAs by multivariate Cox analysis (Table 2). Among them, miR-502 and miR-145 were down-regulated, and miR-142 and miR-33b were up-regulated in cervical cancer tissues. The multivariate Cox analysis demonstrated that these four miRNAs could be used as independent prognostic factors in cervical cancer (P < 0.05).
Table 1 Univariate analysis of cervical cancer patients.
Table 2 Multivariate analysis of cervical cancer patients.
Application of a Cox Proportional Hazard Regression Model in Disease Prognosis
The prognostic model consisting of four miRNAs was constructed by multivariate cox analysis. The formula of the model was defined as follows: Risk score = (Exp (miR-502) × (−0.676) + Exp (miR-145) × (−0.290) + Exp (miR-142) × (−0.297) + Exp (miR-33b) × (−0.209)). The survival risk score based on the four miRNAs was calculated by the model formula and divided patients into low- (154 samples) and high-risk (153 samples) groups according to their median risk score. The risk score curve demonstrates that the risk score of individuals in the low-risk group is small and consistent while the risk score of individuals in the high-risk group is larger and rises significantly (Figure 3A top). The survival status map demonstrates that individuals in the low-risk group have low mortality rates, and individuals in the high-risk group have a short survival time (most survived for less than 5 years) (Figure 3A middle). The heatmap demonstrates that low miRNA expression levels are associated with high-risk scores, and high miRNA expression levels are associated with low-risk scores (Figure 3A bottom). These all demonstrate the reliability of the prognostic model in classifying high- and low-risk groups. And the survival curve demonstrates that the survival rate in the low-risk group is notably higher than in the high-risk group (P = 3E-06; low-risk group: 80.3%, 95% CI = 71.3–90.4%; high-risk group: 51.6%; 95% CI = 41.6–64.0%) (Figure 3B). The time-dependent ROC curve analysis indicates that this prognostic model has a high level of credibility, sensitivity, and specificity with an AUC value of 0.766 (Figure 3C).
Figure 3 Construction of a prognostic model based on a 4-miRNA signature. (A) From top to bottom: the risk score curve, survival status map, and expression heatmap between the low- and high-risk groups. The color bar shows the relative miRNA expression value with red indicating high expression and green indicating low expression. (B) Survival curve for the low- and high-risk groups. (C) The ROC curve for survival predictions.
Target Gene Prediction for the Four Prognostic miRNAs
The prediction results combined with TargetScan and miRDB identified 267 target genes of miR-502, 423 target genes of miR-145, 370 target genes of miR-142, and 279 target genes of miR-33b (Figure S1). After excluding 111 repetitive target genes co-regulated by multiple miRNAs, 1,228 target genes of the four prognostic miRNAs were obtained. Moreover, 94 consensus genes (Table 3) were identified from the overlap between 1840 DEGs (Table S3) and 1228 target genes (Table S4).
Table 3 Consensus genes of four prognostic miRNAs discovered in the overlap between predicted target genes and DEGs related to cervical cancer lymph node metastasis.
KEGG Pathway and GO Function Annotation Details of Consensus Genes
Enrichment analyses produced 7 KEGG pathways and 19 GO function annotations from the consensus genes using the screening criteria of P < 0.05 (Figure 4). Consensus genes enriched in each pathway or function annotations were showed in Table S5. KEGG pathways involving insulin secretion, cAMP signaling, adrenergic signaling in cardiomyocytes, leukocyte transendothelial migration, etc. were mainly enriched. Among them, the insulin secretion pathway had the smallest P-value, and this pathway is associated with five consensus genes (P = 1.05E-03). The cAMP signaling pathway had the largest number (count = 6) of involved consensus genes, and its P-value was 3.87E-03. BP GO annotations involved in neuron migration, positive regulation of the G2/M mitotic cell cycle transition, extracellular matrix disassembly, cell adhesion, etc. were mainly enriched. The annotation with the smallest P-value in the BP involved neuron migration (P = 2.07E-03), and it correlated with five consensus genes. The largest number of related consensus genes involved cell adhesion (count = 7), and its P-value was 3.03E-02. And GO annotations in CC involving the plasma membrane, actin cytoskeleton, postsynaptic density, etc. were mainly enriched. The annotation with the smallest P-value in the CC was an integral component of the plasma membrane (P = 6.92E-03), and it involved 15 consensus genes. The largest number of linked consensus genes involved the plasma membrane (count = 28), and its P-value was 4.33E-02. Furthermore, GO annotations in MF involving heparan sulfate 6-O-sulfotransferase activity were mainly enriched and had a P-value of 1.98E-02 and two associated consensus genes.
Figure 4 Bar graph illustrating the enrichment analysis. The abscissa represents the number of consensus genes involved in KEGG pathways or GO function annotations. The ordinate represents items of the primary KEGG pathways or GO function annotations.
Identification of Key Modules Using a Cluster Analysis of Protein-Protein Interactions
The interaction network presents protein-protein interactions, the strength of the interactions within the protein modules (Figure 5A), and up/down regulation conditions of the involved genes (Figure 5B). It consists of 52 nodes and 59 edges. Two important modules with MCODE scores greater than 2.0 were screened by MCODE and included seven key genes: CXCL12, IGF1, PTPRC, and CDH5 (Figure 5C); RAD51B, REV3L, and WDHD1 (Figure 5D). Among them, CXCL12, IGF1, CDH5, RAD51B, and REV3L were significantly down-regulated while WDHD1 and PTPRC were notably up-regulated. It is worth to note that we also performed topology parameters analysis of PPI and obtained the same seven key genes above (Figure S2). The readers who are interested in this process can read the supplementary TopologyParametersAnalysis.doc.
Figure 5 Diagram of the protein-protein interaction network. The dark and light shading of the lines between solid circles indicates high and low interaction relationships, respectively, which is represented by the Combined Score. (A) The color of solid circles representing protein targets is distinguished by the MCODE score value. The blue circles are protein targets with an MCODE score value of 0, and these proteins did not participate in the construction of the key module. The dark and light orange circles are protein targets with a high MCODE score value of 3.000 and 2.000, respectively, and these proteins participated in the construction of the key module. (B) The color of solid circles representing protein targets is distinguished by the logFC value. The green circles are protein targets with a low logFC value and represent significantly down-regulated genes. The red circles are protein targets with a high logFC value and represent significantly up-regulated genes. (C) A module was constructed using the four protein targets with MCODE score values of 3.000. (D) A module was constructed using the three protein targets with MCODE score values of 2.000.
Visualization Network of MiRNAs-Genes-Pathways and Annotations
The visualization network demonstrated that seven key genes are regulated by four miRNAs and are involved in seven key KEGG pathways and GO functional annotations (Figure 6). WDHD1 and CDH5 are regulated by miR-502; RAD51B and REV3L are regulated by miR-145; REV3L and IGF1 are regulated by miR-142; and CXCL12 and PTPRC are regulated by miR-33b. The enrichment analysis of the consensus genes showed that there was no key gene enrichment in the key KEGG pathways, while CXCL12, IGF1, CDH5, and PTPRC were involved in the key GO functional annotations. Specifically, CXCL12 was enriched in neuron migration and cell adhesion of BP. CDH5 was enriched in cell adhesion of BP. IGF1, together with CDH5 and PTPRC, were enriched in the plasma membrane of CC. PTPRC was enriched in the integral component of the plasma membrane of CC.
Figure 6 MiRNAs-Genes-Pathways and Annotations visualization network denoting the relationships between miRNAs, consensus genes, key genes, key KEGG pathways, and GO functional annotations. Yellow solid circles represent miRNAs, blue solid circles represent consensus genes, red solid circles represent key genes, and green solid circles represent KEGG pathways and GO functional annotations.
Discussion
Cervical carcinoma is one of the most common malignancies in females. Approximately 90% of cervical cancers occur in low- and middle-income countries. Lymph node metastasis and recurrence are the main manifestations in cervical cancer patients with a poor prognosis (Cohen et al., 2019). miRNAs collectively regulate thousands of human cancer-related, protein-coding genes and regulate many important biological processes that facilitate cancer development (i.e. cell growth, invasion, and apoptosis) (Lu et al., 2005). Therefore, miRNAs have become a hotspot of tumor research in the past decade. To identify novel and reliable prognostic biomarkers and important regulatory genes related to cervical cancer and lymph node metastasis, this study first identified 110 DEMs and 1840 DEGs, separately, from expression profiling and clinical data uploaded into TCGA and GEO databases, respectively. Next, a cancer prognosis model based on four prognostic miRNAs (miR-502, miR-145, miR-142, and miR-33b) was established using univariate Cox, multivariate Cox, and survival analyses. Subsequently, 1,228 target genes were predicted using TargetScan and miRDB, and 94 consensus genes were obtained from the overlap of DEGs and predicted target genes. Finally, seven key genes related to lymph node metastasis (i.e. CXCL12, IGF1, PTPRC, CDH5, RAD51B, REV3L, and WDHD1) were identified using STRING and Cytoscape.
According to the prognostic model, the expression levels of miR-502 and miR-145 were down-regulated, and the expression levels of miR-142 and miR-33b were up-regulated in cervical cancer. Among these four miRNAs, miR-145 and miR-142 have been reported that they were related to cervical cancer in previous experimental studies.
Previous studies have shown that the down-regulation of miR-145 is closely related to cervical cancer and its lymph node metastasis. Azizmohammadi et al. (2017) used qRT-PCR and a multivariate Cox analysis to demonstrate that the expression of miR-145 is reduced in cervical cancer tissues. They also showed that reduced expression of miR-145 is related to lymph node metastasis (P = 0.02), advanced Federation International of Gynecology and Obstetrics (FIGO) stage (P = 0.007), and vascular invasion (P = 0.026), which confirms miR-145’s potential as a prognostic biomarker for the early detection of cervical cancer (Azizmohammadi et al., 2017). Ma et al. (2019) demonstrated that miR-145 is also down-regulated in cervical tumor cells, and up-regulation of miR-145 reduces cell proliferation by directly suppressing its target gene, FSCN1. Shi et al. (2012) also showed miR-145 is decreased in cervical cancer cells, and increasing miR-145 expression enhances chemo-sensitivity and inhibits invasion and migration by enhancing p53. The above experimental results are consistent with our prediction that the expression of miR-145 is suppressed in cervical cancer cells.
While our prognostic model predicted that miR-142 expression is increased, existing experimental studies illustrated that miR-142 expression is decreased in cervical cancer cells. Li et al. (2017) revealed that when compared with normal tissue, miR-142 expression is lower in cervical cancer cells and correlates with a poor prognosis. Deng et al. (2015) demonstrated that miR-142-3p is down-regulated in cervical tumors. The overregulation of miR-142-3p inhibits the expression of its target gene, FZD7, and further halts the proliferation and invasion of cervical cancer (Deng et al., 2015). Xia et al. (2018) reported that Metformin, an anti-cancer drug, up-regulates miR-142-3p expression in cervical cancer cells. They also showed that Metformin inhibits the invasion and migration of tumor cells by decreasing the sponge effect of MALAT1, up-regulating the expression of miR-142, and down-regulating the expression of the target gene, HMGA2 (Xia et al., 2018). These three studies showed that miR-142 is a tumor suppressor gene. Thus, the decreased expression of miR-142 in cervical cancer cells is contrary to our results from the data mining calculation that shows an increase in its expression. We do not know the reason(s) for this contradiction, but it should be clarified in the future.
Although no experiment has demonstrated a correlation between miR-502 and miR-33b in cervical cancer, studies have shown that the down-regulation of miR-502 and up-regulation of miR-33b are involved in other types of cancer. Sun et al. (2014) demonstrated that miR-502, which is down-regulated in breast cancer cells, suppresses early apoptosis by targeting TRAF2 and restraining the NF-κB signaling pathway. Furthermore, Li et al. (2009) demonstrated that the activation of NF-κB results in a lower tumor grade, larger tumor volume, higher invasiveness, and increased metastasis in cervical cancer tissues. Therefore, miR-502 may affect cervical cancer lymph node metastasis by participating in the NF-κB signaling pathway. Zhang et al. (2019) showed that the up-regulation of miR-33b inhibits the Wnt/β-catenin signaling pathway by decreasing ZEB1 expression and promoting endometriosis. Although endometriosis is a common and benign disease, it has similar characteristics to malignancies including cell proliferation, invasion, metastasis, and recurrence. Also, Ramachandran et al. (2012) showed that abnormal activation of the Wnt/β-catenin pathway is common in cervical tumors, which may enhance proliferation and prevent apoptosis of cervical cancer cells. Therefore, up-regulation of miR-33b may affect the proliferation and apoptosis of cervical cancer cells by inhibiting the Wnt/β-catenin signaling pathway. To summarize, miR-502 and miR-33b might be involved in cervical cancer formation through distinctive ways as explained above, but such speculation needs to be experimentally validated.
To date, one computational paper exploring the relationship between cervical cancer and the expression of miRNA through data mining has been published. Liang et al. (2017) constructed a three-miRNA signature containing miR-145, miR-200c, and miR-218-1 by processing miRNA data from TCGA database. Their study demonstrated that the expression of miR-145 is significantly decreased in cervical tumor tissues. Furthermore, they found that the expression levels of miR-142 and miR-33b are up-regulated (Liang et al., 2017). This study provided results consistent with ours regarding the expression levels of miR-145, miR-142, and miR-33b. It is worth noting that although the expression of miR-142 was experimentally shown to be down-regulated in cervical cancer cells, as we mentioned previously, this earlier computational study circumstantially confirms the results of our research.
According to our predictions, among the seven key genes involved in early-stage cervical cancer lymph node metastasis, CXCL12, IGF1, RAD51B, REV3L, and CDH5 are down-regulated, and WDHD1 and PTPRC are up-regulated. Table 4 compares the expression of these seven key genes in cervical cancer, other cancers, and lymph node metastasis as reported in previous experimental studies and our calculated results. The decreased expression of CXCL12 in cervical cancer cells and its role in lymph node metastasis was confirmed in previous experiments (Yadav et al., 2016), and the decreased expression of IGF1 and the increased expression of WDHD1 were also experimentally validated in cervical cancer cells but only in lymph node metastasis of other cancers (Serrano et al., 2006; Kümmel et al., 2007; Huang et al., 2008; Zhou et al., 2016; Liu et al., 2019). Interestingly, the expression of RAD51B is also decreased in cervical cancer cells, but experiments showed that its expression is increased in other cancer lymph node metastases (Cheng et al., 2016; Hang et al., 2016). In addition, REV3L expression is increased in cervical cancer cells and lymph node metastases of other cancers in the laboratory, which is inconsistent with our prediction results (Yang et al., 2015; Zhu et al., 2016). Finally, PTPRC and CDH5 have not been previously reported to be associated with cervical cancer or its lymph node metastasis. While experiments demonstrate that PTPRC expression is increased in other cancers and their lymph node metastases (Collette et al., 2007; Camacho et al., 2018), CDH5 is up-regulated in other cancers but down-regulated in lymph node metastasis of colorectal cancer (Tacconi et al., 2015; Hung et al., 2016; Higuchi et al., 2017).
Table 4 The expression of key genes reported from previous experimental studies.
CXCL12 is the ligand for the G-protein coupled receptor-like chemokine (C-X-C motif) receptors 4 and 7. It affects many cellular processes such as immune monitoring, inflammatory response, tumor growth, and metastasis (Colamussi et al., 2001; Yadav et al., 2016). Yadav et al. (2016) demonstrated that CXCL12 expression is absent in cervical cancer. They also illustrated that CXCL12 silencing enables cells to evade apoptosis and leads to the progression and metastasis of cervical cancer (Yadav et al., 2016). Meanwhile, Müller et al. (2001) reported that CXCR4 is usually highly expressed in tumor tissues, and the CXCR4/CXCL12 axis could prevent breast tumor lymph-node metastasis and lung metastasis when CXCR4 is neutralized. Thus, CXCL12 can act as a tumor suppressor in lymph node metastasis of cervical and other cancers, and it might inhibit tumor progression through the CXCR4/CXCL12 axis.
IGF1 is a cytokine that mediates cell growth and development (Macaulay, 1992). Serrano et al. (2006) reported low serum IGF1 expression in cervical cancer cells, and reduced expression of IGF1 is associated with an increased risk of cervical cancer. Huang et al. (2008) showed that lower levels of IGF1 can effectively predict survival in patients with cervical cancer. Furthermore, they found that a combination of increased carcinoembryonic antigen (CEA) levels and decreased IGF1 levels is significantly associated with an increased risk of death and could accurately predict patients with a poor prognosis. Though no experiment is available to confirm that IGF1 is involved in cervical cancer lymph node metastasis, Kümmel et al. (2007) reported that in breast tumor patients, plasma IGF1 expression is increased after dose-intensified chemotherapy, and they showed a more prominent increase in IGF1 expression in patients with positive lymph nodes than other patients. WDHD1, also known as AND1, is involved in signal transduction, pre-mRNA processing, replication, transcription, cytoskeleton assembly, chromosome assembly, etc. (Park et al., 2012). Zhou et al. (2016) showed that WDHD1 is up-regulated in primary human keratinocytes and spontaneously immortalized human foreskin keratinocytes cells expressing oncogene E7 in HPV-induced carcinogenesis. They also showed that WDHD1 can increase E7-induced G1 checkpoint abolition and duplication. They further demonstrated that the polyploidy ratio of cells expressing E7 is significantly reduced after knocking down WDHD1 with siRNA (Zhou et al., 2016). Although the role of WDHD1 in cervical cancer lymph node metastasis has not been explained, Liu et al. (Liu et al., 2019) found that miR-494 can inhibit EMT and lymph node metastasis of cholangiocarcinoma (CCA) cells by targeting overexpressed WDHD1. In summary, IGF1 and WDHD1 are likely connected with cervical cancer lymph node metastasis.
RAD51B is one of the RAD51 gene families involved in homologous recombination-mediated DNA repair (Thacker, 2005; Nagathihalli and Nagaraju, 2011). Hang et al. (2016) demonstrated that RAD51B is a tumor suppressor gene in cervical cancer cells, which is consistent with our prediction. They further showed that the miRNA binding sites of RAD51B genetic variants in cervical cancer cells may increase tumor susceptibility, and RAD51B is vital in gauging the cervical cancer risk of individuals and improving the effectiveness of preventive intervention (Hang et al., 2016). Currently, there are no reports that demonstrate RAD51B affects cervical cancer lymph node metastasis, but Cheng et al. (2016) showed that the overexpression of RAD51B in gastric cancer cells is significantly associated with lymph node metastasis (P = 0.001), advanced stage (P = 0.009), and invasive differentiation (P = 0.022), and it may act as a potential signature for early detection and poor prognosis. RAD51B expression in gastric cancer with lymph node metastasis is contrary to our prediction: it tended to represent the particularity of RAD51B expression in cervical cancer lymph node metastasis. The specific reasons for this contradiction need to be further clarified.
REV3L encodes the protein representing the catalytic sub-unit of Polζ, and inhibiting REV3L expression enables cancer cells to tolerate DNA damage and stunted growth (Knobel et al., 2011). Yang et al. (2015) showed that overexpression of REV3L promotes proliferation and colony formation and inhibits cervical cancer cell sensitivity to cisplatin. Thus, REV3L could be a potential therapeutic target for cervical cancer treatment. Zhu et al. (2016) also showed that REV3L is significantly up-regulated in esophageal squamous cell carcinoma tissues, and it positively correlates with lymph node metastasis (P <0.05) and clinical stage (P < 0.05). Additionally, overexpression of REV3L increases the expression levels of cyclin D1 and survivin, which work together to promote the growth and invasion of esophageal cancer cells (Zhu et al., 2016). These results demonstrate that REV3L is closely linked to cervical cancer lymph node metastasis. However, these REV3L experimental results are inconsistent with the inferences obtained from our data mining analysis. The specific reasons behind this contradiction need to be clarified.
PTPRC, also known as CD45, is a key regulator of cell growth, differentiation, mitosis, and carcinogenic transformation (Rheinländer et al., 2018). Camacho et al. (2018) showed that PTPRC is significantly overexpressed in head and neck squamous cell cancer cells, and tumor samples overexpressing PTPRC have significantly higher tumor-infiltrating lymphocytes (TIL) scores than tumor samples expressing low levels of PTPRC, leading to a poorer prognosis. Collette et al. (2007) noted that the expression of PTPRC plays a critical role in determining the signal transduction and proliferation response of human myeloma cells to growth factors such as IL-6 and IGF1. IL-6 and IGF1 separately induced CD45+ and CD45- myeloma cell colony formation through the MAPK/ERK signaling pathway in which CD45 is critical for myeloma proliferation (Collette et al., 2007). In addition, many experimental studies reported that the MAPK/ERK signaling pathway is crucial in cervical cancer formation. For example, Li et al. (2018) demonstrated that MAPK/ERK signaling pathway activation promotes cervical cancer cell proliferation. Tao et al. (2018) illustrated that miR-497 acts as a tumor suppressor by blocking the MAPK/ERK signaling pathway in cervical cancer cells. Furthermore, Wang et al. (2020) showed that the anti-cancer drug sclareol targets the MAPK/ERK signaling pathway and induces cervical cancer cell apoptosis and cell cycle arrest. Therefore, considering these observations, we can speculate that PTPRC and IGF1 may negatively regulate the occurrence and development of cervical cancer and its lymph node metastasis by affecting the MAPK/ERK signaling pathway.
CDH5, also known as VE-cadherin, represses endothelial cell apoptosis and participates in endothelial cell growth contact inhibition (Cavallaro et al., 2006). Abnormal CDH5 was found in cancer cells. Hung et al. (2016) and Higuchi et al. (2017) demonstrated that CDH5 is over-expressed in lung cancer and gastric cancer, respectively. However, Tacconi et al. (2015) demonstrated that CDH5 is down-regulated in colorectal cancer, and is negatively associated with its lymph node metastasis. They found that up-regulated VEGFC/VEGFR3 reduces CDH5 expression, enhances permeability, and increases trans-endothelial migration. Thus, it promotes lymphatic vessel density and colorectal cancer lymph node metastasis (Tacconi et al., 2015). Furthermore, Chaudary et al. (2011) demonstrated that suppressing the over-expression of VEGFC/VEGFR3 in cervical cancer cells inhibits proliferation and cervical cancer lymph node metastasis. We suspect that the over-expression of VEGFC/VEGFR3 restrains CDH5 in cervical cancer cells and further promotes lymph node metastasis. But this assumption needs to be further studied.
Table 3 shows that eight targeting relationships are available between four prognosis miRNAs and seven key genes. Specifically, miR-502 targets WDHD1 and CDH5, miR-145 targets REV3L and RAD51B, miR-142 targets REV3L and IGF1, and miR-33b targets CXCL12 and PTPRC. Among these relationships, miR-145-REV3L and miR-142-IGF1 have been reported previously. Chen et al. (2019) demonstrated that miR-145 can directly regulate the expression of REV3L in esophageal squamous cell carcinoma (ESCC) cells using a dual luciferase reporter assay and Western blot analysis. Also, low levels of miR-145 augmented REV3L mRNA and protein expression in ESCC. Xiong et al. (2018) showed that eight out of 12 target gene prediction programs including TargetScan and miRDB predicted the miR-142-IGF1 relationship. Subsequently, they suggested that decitabine exerts its therapeutic effect on hepatocellular carcinoma by inhibiting miR-142 DNA methylation, which enhances miRNA expression and further down-regulates the target genes of miR-142-5p such as IGF1 (Xiong et al., 2018). To summarize, the above studies demonstrated target/regulatory relationships between miR-145-REV3L and miR-142-IGF1. The remaining six regulatory relationships have yet to be reported, which might lay the foundation for future research on cervical and other cancers with lymph node metastasis.
In the key KEGG pathways and GO functional annotations of key gene enrichment, it had been reported that cell adhesion and integral component of the plasma membrane were closely related to tumor proliferation and metastasis. For example, Mierke (2008) showed that the integrity of the endothelial barrier is preserved by a complex balance of cell adhesion factors. Once the integrity is compromised, tumor cells can metastasize through blood vessels or lymphatic channels (Mierke, 2008). More importantly, Zheng et al. (2017) suggested that the CXCR4/CXCL12 axis plays a role in reducing the adhesion ability of colon cancer cells by regulating the Akt and IGF1R signaling pathways. Therefore, the key gene CXCL12, enriched in cell adhesion, may lead to cervical cancer lymph node metastasis by affecting cell-cell adhesion. In addition, protein components of the cell membrane are involved in cell signal transduction, cell interaction, and other important steps in the process of cancer cell metastasis (Lund et al., 2009). Most of the key molecules involved in the signaling pathway are distributed in the cell membrane. PTPRC, as a transmembrane tyrosine phosphatase, is expressed in all leukocytes, and it is involved in lymphocyte immunity against tumor cells (Tchilian et al., 2001). Therefore, we speculate that PTPRC, which is enriched in an integral component of the plasma membrane, may be an important signal molecule in cervical cancer metastasis. Moreover, it may be involved in a new approach to studying the mechanism of cervical cancer lymph node metastasis.
Conclusion
By data-mining differentially expressed miRNAs and genes along with other clinical information in a multi-step analysis, we obtained a prognostic model of cervical cancer with lymph node metastasis containing four miRNAs and seven genes. We showed that four miRNAs (miR-502, miR-145, miR-142, and miR-33b) are independent and common prognostic biomarkers for patients with cervical cancer, and seven proteins (CXCL12, IGF1, PTPRC CDH5, RAD51B, REV3L, and WDHD1) are key genes significantly related to lymph node metastasis. Among them, miR-145, miR-142, CXCL12, IGF1, and WDHD1 have been confirmed, while miR-502, miR-33b, PTPRC, and CDH5 are reported, here, for the first time. Also, the expression levels and/or roles of miR-142, miR-33b, RAD51B, REV3L, and CDH5 in cervical cancer lymph node metastasis need further clarification. In summary, our study may improve our understanding of the progression and lymph node metastasis mechanism of cervical cancer, and, therefore, provide a novel diagnostic indicator and therapeutic targets for future clinical treatments.
Data Availability Statement
Publicly available datasets were analyzed in this study. These data can be found here: https://portal.gdc.cancer.gov/, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE7410.
Author Contributions
Conceptualization: SC, CG, and ZH. Methodology: SC, CG, and YW. Software: SC, CG, and YW. Validation: SC, YW, and ZH. Formal analysis: SC, CG, and ZH. Investigation: SC and YW. Data curation: SC and CG. Writing—original draft preparation: SC and CG. Writing—review and editing: ZH. Visualization: SC and CG. Supervision: ZH. Project administration: ZH. Funding acquisition: ZH. All authors have read and agreed to the published version of the manuscript.
Funding
This work was supported by the National Natural Science Foundation of China (31770774), the Provincial Major Project of Basic or Applied Research in Natural Science, Guangdong Provincial, Education Department (2016KZDXM038) and the higher education reform project of Guangdong Province (2019268).
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.
Acknowledgments
We thank Wordvice for their help in revising the English grammar.
Supplementary Material
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphar.2020.00544/full#supplementary-material
Abbreviations
AUC, area under the curve; BP, biological process; CC, cellular component; CDH5, cadherin 5; CEA, carcinoembryonic antigen; CI, confidence interval; Coef, partial regression coefficient; CXCL12, C-X-C motif chemokine ligand 12; CXCR4, C-X-C motif chemokine receptor 4; CXCR7, C-X-C motif chemokine receptor 7; DAVID, Database for Annotation, Visualization, and Integrated Discovery; DEGs, differentially expressed genes; DEMs, differentially expressed miRNAs; EMT, epithelial-mesenchymal transition; Exp (Coef), relative risk; FDR, false discovery rate; GDF15, growth differentiation factor 15; GEO, Gene Expression Omnibus; GO, gene ontology; HMGA2, high-mobility group AT-hook 2; HR, hazard ratio; IGF1, insulin like growth factor 1; IGF1R, insulin like growth factor 1 receptor; KEGG, Kyoto Encyclopedia of Genes and Genome; LNM, lymph node metastasis; MALAT1, metastasis-associated lung adenocarcinoma transcript-1; MCODE, Molecular Complex Detection; MF, molecular function; Pol ζ, DNA translation synthesis polymerase ζ; PPI, protein-protein interaction; PTPRC, protein tyrosine phosphatase receptor type C; qRT-PCR, quantitative real-time polymerase chain reaction (reverse transcription-quantitative polymerase chain reaction (RT-qPCR)); RAD51B, RAD51 paralog B; REV3L, recovery protein 3 like; Se (Coef), standard error of regression coefficient; siRNA, small interfering RNA; STRING, Search Tool for the Retrieval of Interacting Genes; TCGA, The Cancer Genome Atlas; TNM, tumor node metastasis; WDHD1, WD repeat and HMG-box DNA binding protein 1.
References
Agarwal, V., Bell, G. W., Nam, J.-W., Bartel, D. P. (2015). Predicting effective microRNA target sites in mammalian mRNAs. eLife 4, e05005. doi: 10.7554/eLife.05005
Ashburner, M., Ball, C. A., Blake, J. A., Botstein, D., Butler, H., Cherry, J. M., et al. (2000). Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat. Genet. 25 (1), 25–29. doi: 10.1038/75556
Azizmohammadi, S., Safari, A., Azizmohammadi, S., Kaghazian, M., Sadrkhanlo, M., Yahaghi, E., et al. (2017). Molecular identification of miR-145 and miR-9 expression level as prognostic biomarkers for early-stage cervical cancer detection. J. Assoc. Physicians 110 (1), 11–15. doi: 10.1093/qjmed/hcw101
Bader, G. D., Hogue, C. W. V. (2003). An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinf. 4, 2. doi: 10.1186/1471-2105-4-2
Barrett, T., Troup, D. B., Wilhite, S. E., Ledoux, P., Rudnev, D., Evangelista, C., et al. (2007). NCBI GEO: mining tens of millions of expression profiles–database and tools update. Nucleic Acids Res. 35 (Database issue), D760–D765. doi: 10.1093/nar/gkl887
Bartel, D. P. (2004). MicroRNAs: genomics, biogenesis, mechanism, and function. Cell 116 (2), 281–297. doi: 10.1016/s0092-8674(04)00045-5
Bray, F., Ferlay, J., Soerjomataram, I., Siegel, R. L., Torre, L. A., Jemal, A. (2018). Global cancer statistics 2018: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. Cancer J. Clin. 68 (6), 394–424. doi: 10.3322/caac.21492
Camacho, M., Agüero, A., Sumarroca, A., López, L., Pavón, M.-Á., Avilés-Jurado, F.-X., et al. (2018). Prognostic value of CD45 transcriptional expression in head and neck cancer. Eur. Arch. Oto-Rhino-Laryngol 275 (1), 225–232. doi: 10.1007/s00405-017-4806-2
Cavallaro, U., Liebner, S., Dejana, E. (2006). Endothelial cadherins and tumor angiogenesis. Exp. Cell Res. 312 (5), 659–667. doi: 10.1016/j.yexcr.2005.09.019
Chaudary, N., Milosevic, M., Hill, R. P. (2011). Suppression of vascular endothelial growth factor receptor 3 (VEGFR3) and vascular endothelial growth factor C (VEGFC) inhibits hypoxia-induced lymph node metastases in cervix cancer. Gynecol. Oncol. 123 (2), 393–400. doi: 10.1016/j.ygyno.2011.07.006
Chen, L., Kang, C. (2015). miRNA interventions serve as ‘magic bullets’ in the reversal of glioblastoma hallmarks. Oncotarget 6 (36), 38628–38642. doi: 10.18632/oncotarget.5926
Chen, Q., Hou, J., Wu, Z., Zhao, J., Ma, D. (2019). miR-145 Regulates the sensitivity of esophageal squamous cell carcinoma cells to 5-FU via targeting REV3L. Pathol. Res. Pract. 215 (7), 152427. doi: 10.1016/j.prp.2019.04.019
Cheng, Y., Yang, B., Xi, Y., Chen, X. (2016). RAD51B as a potential biomarker for early detection and poor prognostic evaluation contributes to tumorigenesis of gastric cancer. J. Int. Soc. Oncodevelopmental Biol. Med. 37 (11), 14969–14978. doi: 10.1007/s13277-016-5340-3
Cohen, P. A., Jhingran, A., Oaknin, A., Denny, L. (2019). Cervical cancer. Lancet (London England) 393 (10167), 169–182. doi: 10.1016/S0140-6736(18)32470-X
Colamussi, M. L., Secchiero, P., Gonelli, A., Marchisio, M., Zauli, G., Capitani, S. (2001). Stromal derived factor-1 alpha (SDF-1 alpha) induces CD4+ T cell apoptosis via the functional up-regulation of the Fas (CD95)/Fas ligand (CD95L) pathway. J. Leukocyte Biol. 69 (2), 263–270.
Collette, M., Descamps, G., Pellat-Deceunynck, C., Bataille, R., Amiot, M. (2007). Crucial role of phosphatase CD45 in determining signaling and proliferation of human myeloma cells. Eur. Cytokine Network 18 (3), 120–126. doi: 10.1684/ecn.2007.0095
Dai, Y., Tong, R., Guo, H., Yu, T., Wang, C. (2017). Association of CXCR4, CCR7, VEGF-C and VEGF-D expression with lymph node metastasis in patients with cervical cancer. Eur. J. Obstetrics Gynecol. Reprod. Biol. 214, 178–183. doi: 10.1016/j.ejogrb.2017.04.043
Deng, B., Zhang, Y., Zhang, S., Wen, F., Miao, Y., Guo, K. (2015). MicroRNA-142-3p inhibits cell proliferation and invasion of cervical cancer cells by targeting FZD7. J. Int. Soc. Oncodevelopmental Biol. Med. 36 (10), 8065–8073. doi: 10.1007/s13277-015-3483-2
Doorbar, J. (2006). Molecular biology of human papillomavirus infection and cervical cancer. Clin. Sci. (London Engl: 1979) 110 (5), 525–541. doi: 10.1042/CS20050369
Gómez-Gómez, Y., Organista-Nava, J., Gariglio, P. (2013). Deregulation of the miRNAs expression in cervical cancer: human papillomavirus implications. BioMed. Res. Int. 2013, 407052. doi: 10.1155/2013/407052
George, B., Seals, S., Aban, I. (2014). Survival analysis and regression models. J. Nuclear Cardiol. 21 (4), 686–694. doi: 10.1007/s12350-014-9908-2
Hang, D., Zhou, W., Jia, M., Wang, L., Zhou, J., Yin, Y., et al. (2016). Genetic variants within microRNA-binding site of RAD51B are associated with risk of cervical cancer in Chinese women. Cancer Med. 5 (9), 2596–2601. doi: 10.1002/cam4.797
Higuchi, K., Inokuchi, M., Takagi, Y., Ishikawa, T., Otsuki, S., Uetake, H., et al. (2017). Cadherin 5 expression correlates with poor survival in human gastric cancer. J. Clin. Pathol. 70 (3), 217–221. doi: 10.1136/jclinpath-2016-203640
Huang, B.-X., Fang, F. (2018). Progress in the Study of Lymph Node Metastasis in Early-stage Cervical Cancer. Curr. Med. Sci. 38 (4), 567–574. doi: 10.1007/s11596-018-1915-0
Huang, Y. F., Shen, M. R., Hsu, K. F., Cheng, Y. M., Chou, C. Y. (2008). Clinical implications of insulin-like growth factor 1 system in early-stage cervical cancer. Br. J. Cancer 99 (7), 1096–1102. doi: 10.1038/sj.bjc.6604661
Huang, D. W., Sherman, B. T., Lempicki, R. A. (2009). Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat. Protoc. 4 (1), 44–57. doi: 10.1038/nprot.2008.211
Huang, Z., Zhu, D., Wu, L., He, M., Zhou, X., Zhang, L., et al. (2017). Six Serum-Based miRNAs as Potential Diagnostic Biomarkers for Gastric Cancer. Cancer Epidemiol. Biomarkers Prev. 26 (2), 188–196. doi: 10.1158/1055-9965.EPI-16-0607
Hung, M.-S., Chen, I. C., Lung, J.-H., Lin, P.-Y., Li, Y.-C., Tsai, Y.-H. (2016). Epidermal Growth Factor Receptor Mutation Enhances Expression of Cadherin-5 in Lung Cancer Cells. PloS One 11 (6), e0158395. doi: 10.1371/journal.pone.0158395
Kümmel, S., Eggemann, H., Lüftner, D., Gebauer, N., Bühler, H., Schaller, G., et al. (2007). Significant changes in circulating plasma levels of IGF1 and IGFBP3 after conventional or dose-intensified adjuvant treatment of breast cancer patients with one to three positive lymph nodes. Int. J. Biol. Markers 22 (3), 186–193.
Kanehisa, M., Goto, S. (2000). KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 28 (1), 27–30. doi: 10.1093/nar/28.1.27
Knobel, P. A., Kotov, I. N., Felley-Bosco, E., Stahel, R. A., Marti, T. M. (2011). Inhibition of REV3 expression induces persistent DNA damage and growth arrest in cancer cells. Neoplasia (N.Y.) 13 (10), 961–970. doi: 10.1593/neo.11828
Lee, Y. S., Dutta, A. (2009). MicroRNAs in cancer. Annu. Rev. Pathol. 4, 199–227. doi: 10.1146/annurev.pathol.4.110807.092222
Li, J., Jia, H., Xie, L., Wang, X., Wang, X., He, H., et al. (2009). Association of constitutive nuclear factor-kappaB activation with aggressive aspects and poor prognosis in cervical cancer. Int. J. Gynecol. Cancer 19 (8), 1421–1426. doi: 10.1111/IGC.0b013e3181b70445
Li, M., Li, B. Y., Xia, H., Jiang, L. L. (2017). Expression of microRNA-142-3p in cervical cancer and its correlation with prognosis. Eur. Rev. Med. Pharmacol. Sci. 21 (10), 2346–2350.
Li, S., Ma, Y.-M., Zheng, P.-S., Zhang, P. (2018). GDF15 promotes the proliferation of cervical cancer cells by phosphorylating AKT1 and Erk1/2 through the receptor ErbB2. J. Exp. Clin. Cancer Res. 37 (1), 80. doi: 10.1186/s13046-018-0744-0
Li, S., Cui, Z., Zhao, Y., Ma, S., Sun, Y., Li, H., et al. (2020). Candidate lncRNA-microRNA-mRNA networks in predicting non-small cell lung cancer and related prognosis analysis. J. Cancer Res. Clin. Oncol. 146 (4), 883–896. doi: 10.1007/s00432-020-03161-6
Liang, B., Li, Y., Wang, T. (2017). A three miRNAs signature predicts survival in cervical cancer using bioinformatics analysis. Sci. Rep. 7 (1), 5624. doi: 10.1038/s41598-017-06032-2
Liu, B., Hu, Y., Qin, L., Peng, X.-B., Huang, Y.-X. (2019). MicroRNA-494-dependent WDHDI inhibition suppresses epithelial-mesenchymal transition, tumor growth and metastasis in cholangiocarcinoma. Digestive Liver Dis. 51 (3), 397–411. doi: 10.1016/j.dld.2018.08.021
Lu, J., Getz, G., Miska, E. A., Alvarez-Saavedra, E., Lamb, J., Peck, D., et al. (2005). MicroRNA expression profiles classify human cancers. Nature 435 (7043), 834–838. doi: 10.1038/nature03702
Lund, R., Leth-Larsen, R., Jensen, O. N., Ditzel, H. J. (2009). Efficient isolation and quantitative proteomic analysis of cancer cell plasma membrane proteins for identification of metastasis-associated cell surface markers. J. Proteome Res. 8 (6), 3078–3090. doi: 10.1021/pr801091k
Müller, A., Homey, B., Soto, H., Ge, N., Catron, D., Buchanan, M. E., et al. (2001). Involvement of chemokine receptors in breast cancer metastasis. Nature 410 (6824), 50–56. doi: 10.1038/35065016
Ma, L., Li, L.-L. (2019). miR-145 Contributes to the Progression of Cervical Carcinoma by Directly Regulating FSCN1. Cell Transplant. 28 (9-10), 1299–1305. doi: 10.1177/0963689719861063
Macaulay, V. M. (1992). Insulin-like growth factors and cancer. Br. J. Cancer 65 (3), 311–320. doi: 10.1038/bjc.1992.65
Mierke, C. T. (2008). Role of the endothelium during tumor cell metastasis: is the endothelium a barrier or a promoter for cell invasion and metastasis? J. Biophysics (Hindawi Publishing Corporation Online) 2008, 183516. doi: 10.1155/2008/183516
Nagathihalli, N. S., Nagaraju, G. (2011). RAD51 as a potential biomarker and therapeutic target for pancreatic cancer. Biochim. Biophys. Acta 1816 (2), 209–218. doi: 10.1016/j.bbcan.2011.07.004
Park, S.-Y., Im, J.-S., Park, S.-R., Kim, S.-E., Wang, H.-J., Lee, J.-K. (2012). Mimosine arrests the cell cycle prior to the onset of DNA replication by preventing the binding of human Ctf4/And-1 to chromatin via Hif-1α activation in HeLa cells. Cell Cycle (Georgetown Tex.) 11 (4), 761–766. doi: 10.4161/cc.11.4.19209
Park, S., Eom, K., Kim, J., Bang, H., Wang, H.-Y., Ahn, S., et al. (2017). MiR-9, miR-21, and miR-155 as potential biomarkers for HPV positive and negative cervical cancer. BMC Cancer 17 (1), 658. doi: 10.1186/s12885-017-3642-5
Ramachandran, I., Thavathiru, E., Ramalingam, S., Natarajan, G., Mills, W. K., Benbrook, D. M., et al. (2012). Wnt inhibitory factor 1 induces apoptosis and inhibits cervical cancer growth, invasion and angiogenesis in vivo. Oncogene 31 (22), 2725–2737. doi: 10.1038/onc.2011.455
Rheinländer, A., Schraven, B., Bommhardt, U. (2018). CD45 in human physiology and clinical medicine. Immunol. Lett. 196, 22–32. doi: 10.1016/j.imlet.2018.01.009
Ritchie, M. E., Phipson, B., Wu, D., Hu, Y., Law, C. W., Shi, W., et al. (2015). limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 43 (7), e47. doi: 10.1093/nar/gkv007
Robinson, M. D., McCarthy, D. J., Smyth, G. K. (2010). edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinf. (Oxford England) 26 (1), 139–140. doi: 10.1093/bioinformatics/btp616
Sepulveda, J. L. (2020). Using R and Bioconductor in Clinical Genomics and Transcriptomics. J. Mol. Diagn. 22 (1), 3–20. doi: 10.1016/j.jmoldx.2019.08.006
Serrano, M. L., Romero, A., Cendales, R., Sánchez-Gómez, M., Bravo, M. M. (2006). Serum levels of insulin-like growth factor-I and -II and insulin-like growth factor binding protein 3 in women with squamous intraepithelial lesions and cervical cancer. Rev. Del Instituto Nacional Salud 26 (2), 258–268.
Shannon, P., Markiel, A., Ozier, O., Baliga, N. S., Wang, J. T., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13 (11), 2498–2504. doi: 10.1101/gr.1239303
Shi, M., Du, L., Liu, D., Qian, L., Hu, M., Yu, M., et al. (2012). Glucocorticoid regulation of a novel HPV-E6-p53-miR-145 pathway modulates invasion and therapy resistance of cervical cancer cells. J. Pathol. 228 (2), 148–157. doi: 10.1002/path.3997
Sun, L.-L., Wang, J., Zhao, Z.-J., Liu, N., Wang, A.-L., Ren, H.-Y., et al. (2014). Suppressive role of miR-502-5p in breast cancer via downregulation of TRAF2. Oncol. Rep. 31 (5), 2085–2092. doi: 10.3892/or.2014.3105
Szklarczyk, D., Morris, J. H., Cook, H., Kuhn, M., Wyder, S., Simonovic, M., et al. (2017). The STRING database in 2017: quality-controlled protein-protein association networks, made broadly accessible. Nucleic Acids Res. 45 (D1), D362–D368. doi: 10.1093/nar/gkw937
Tacconi, C., Correale, C., Gandelli, A., Spinelli, A., Dejana, E., D’Alessio, S., et al. (2015). Vascular endothelial growth factor C disrupts the endothelial lymphatic barrier to promote colorectal cancer invasion. Gastroenterology 148 (7), 1438–51.e8. doi: 10.1053/j.gastro.2015.03.005
Tao, L., Zhang, C.-Y., Guo, L., Li, X., Han, N.-N., Zhou, Q., et al. (2018). MicroRNA-497 accelerates apoptosis while inhibiting proliferation, migration, and invasion through negative regulation of the MAPK/ERK signaling pathway via RAF-1. J. Cell. Physiol. 233 (10), 6578–6588. doi: 10.1002/jcp.26272
Tchilian, E. Z., Wallace, D. L., Wells, R. S., Flower, D. R., Morgan, G., Beverley, P. C. (2001). A deletion in the gene encoding the CD45 antigen in a patient with SCID. J. Immunol. (Baltimore Md.: 1950) 166 (2), 1308–1313. doi: 10.4049/jimmunol.166.2.1308
Thacker, J. (2005). The RAD51 gene family, genetic instability and cancer. Cancer Lett. 219 (2), 125–135. doi: 10.1016/j.canlet.2004.08.018
Tomczak, K., Czerwińska, P., Wiznerowicz, M. (2015). The Cancer Genome Atlas (TCGA): an immeasurable source of knowledge. Contemporary oncology (Poznan, Poland) 19, 1A, A68–A77. doi: 10.5114/wo.2014.47136
Vrieze, S. I. (2012). Model selection and psychological theory: a discussion of the differences between the Akaike information criterion (AIC) and the Bayesian information criterion (BIC). psychol. Methods 17 (2), 228–243. doi: 10.1037/a0027127
Wang, L., Zhou, P., Gao, X.-M., Luo, M.-L., Meng, M.-J., He, Y., et al. (2020). Naturally Occurring Sclareol Diterpene Augments the Chemosensitivity of Human Hela Cervical Cancer Cells by Inducing Mitochondrial Mediated Programmed Cell Death, S-Phase Cell Cycle Arrest and Targeting Mitogen-Activated Protein Kinase (MAPK)/Extracellular-Signal-Regulated Kinase (ERK) Signaling Pathway. Int. Med. J. Exp. Clin. Res. 26, e920248. doi: 10.12659/MSM.920248
Wang, W., Chu, H. J., Shang, C. L., Gong, X., Liu, T. Y., Zhao, Y. H., et al. (2016). Long-Term Oncological Outcomes After Laparoscopic Versus Abdominal Radical Hysterectomy in Stage IA2 to IIA2 Cervical Cancer: A Matched Cohort Study. Int. J. Gynecol. Cancer 26 (7), 1264–1273. doi: 10.1097/IGC.0000000000000749
Wong, N., Wang, X. (2015). miRDB: an online resource for microRNA target prediction and functional annotations. Nucleic Acids Res. 43 (Database issue), D146–D152. doi: 10.1093/nar/gku1104
Xia, C., Liang, S., He, Z., Zhu, X., Chen, R., Chen, J. (2018). Metformin, a first-line drug for type 2 diabetes mellitus, disrupts the MALAT1/miR-142-3p sponge to decrease invasion and migration in cervical cancer cells. Eur. J. Pharmacol. 830, 59–67. doi: 10.1016/j.ejphar.2018.04.027
Xiong, D.-D., Dang, Y.-W., Lin, P., Wen, D.-Y., He, R.-Q., Luo, D.-Z., et al. (2018). A circRNA-miRNA-mRNA network identification for exploring underlying pathogenesis and therapy strategy of hepatocellular carcinoma. J. Trans. Med. 16 (1), 220. doi: 10.1186/s12967-018-1593-5
Yadav, S. S., Prasad, S. B., Prasad, C. B., Pandey, L. K., Pradhan, S., Singh, S., et al. (2016). CXCL12 is a key regulator in tumor microenvironment of cervical cancer: an in vitro study. Clin. Exp. Metastasis 33 (5), 431–439. doi: 10.1007/s10585-016-9787-9
Yang, L., Shi, T., Liu, F., Ren, C., Wang, Z., Li, Y., et al. (2015). REV3L, a promising target in regulating the chemosensitivity of cervical cancer cells. PloS One 10 (3), e0120334. doi: 10.1371/journal.pone.0120334
Yang, G., Zhang, Y., Yang, J. (2019). A Five-microRNA Signature as Prognostic Biomarker in Colorectal Cancer by Bioinformatics Analysis. Front. Oncol. 9, 1207. doi: 10.3389/fonc.2019.01207
Yao, S., Xu, J., Zhao, K., Song, P., Yan, Q., Fan, W., et al. (2018). Down-regulation of HPGD by miR-146b-3p promotes cervical cancer cell proliferation, migration and anchorage-independent growth through activation of STAT3 and AKT pathways. Cell Death Dis. 9 (11), 1055. doi: 10.1038/s41419-018-1059-y
Zeng, Y., Wang, K.-X., Xu, H., Hong, Y. (2018). Integrative miRNA analysis identifies hsa-miR-3154, hsa-miR-7-3, and hsa-miR-600 as potential prognostic biomarker for cervical cancer. J. Cell. Biochem. 119 (2), 1558–1566. doi: 10.1002/jcb.26315
Zhang, B., Pan, X., Cobb, G. P., Anderson, T. A. (2007). microRNAs as oncogenes and tumor suppressors. Dev. Biol. 302 (1), 1–12. doi: 10.1016/j.ydbio.2006.08.028.
Zhang, Z.-Q., Han, Y.-Z., Nian, Q., Chen, G., Cui, S.-Q., Wang, X.-Y. (2015). Tumor Invasiveness, Not Lymphangiogenesis, Is Correlated with Lymph Node Metastasis and Unfavorable Prognosis in Young Breast Cancer Patients (≤35 Years). PloS One 10 (12), e0144376. doi: 10.1371/journal.pone.0144376
Zhang, H., Li, G., Sheng, X., Zhang, S. (2019). Upregulation of miR−33b promotes endometriosis via inhibition of Wnt/β−catenin signaling and ZEB1 expression. Mol. Med. Rep. 19 (3), 2144–2152. doi: 10.3892/mmr.2019.9870
Zheng, F., Zhang, Z., Flamini, V., Jiang, W. G., Cui, Y. (2017). The Axis of CXCR4/SDF-1 Plays a Role in Colon Cancer Cell Adhesion Through Regulation of the AKT and IGF1R Signalling Pathways. Anticancer Res. 37 (8), 4361–4369. doi: 10.21873/anticanres.11830
Zhou, Y., Zhang, Q., Gao, G., Zhang, X., Liu, Y., Yuan, S., et al. (2016). Role of WDHD1 in Human Papillomavirus-Mediated Oncogenesis Identified by Transcriptional Profiling of E7-Expressing Cells. J. Virol. 90 (13), 6071–6084. doi: 10.1128/JVI.00513-16
Keywords: cervical cancer, miRNA, key gene, prognostic signature, lymph node metastasis
Citation: Chen S, Gao C, Wu Y and Huang Z (2020) Identification of Prognostic miRNA Signature and Lymph Node Metastasis-Related Key Genes in Cervical Cancer. Front. Pharmacol. 11:544. doi: 10.3389/fphar.2020.00544
Received: 27 February 2020; Accepted: 09 April 2020;
Published: 08 May 2020.
Edited by:
Feng Zhu, Zhejiang University, ChinaReviewed by:
Bifang He, Guizhou University, ChinaKangyu Lin, University of Texas MD Anderson Cancer Center, United States
Copyright © 2020 Chen, Gao, Wu and Huang. 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: Zunnan Huang, em5faHVhbmdAeWFob28uY29t
†These authors have contributed equally to this work