Abstract Background and aim Previous studies have suggested that lymph node metastasis (LNM) in early-stage cervical cancer (CESC) may affect the prognosis of patients and the outcomes of subsequent adjuvant therapy. However, research focused on miRNA expression in early-stage CESC patients with LNM remains limited. Therefore, it is necessary to identify prognostic miRNAs and determine their molecular mechanisms. Methods We evaluated the differentially expressed genes in early-stage CESC patients with LNM compared to patients without LNM and evaluated the prognostic significance of these differentially expressed genes by analyzing a public dataset from The Cancer Genome Atlas. Potential molecular mechanisms were investigated by gene ontology, the Kyoto Encyclopedia of Genes and Genomes, and protein–protein interaction network analyses. Results According to the The Cancer Genome Atlas data, hsa-miR-508, hsa-miR-509-2, and hsa-miR-526b expression levels were significantly lower in early-stage CESC patients with LNM than in patients without LNM. A multivariate analysis suggested that three miRNAs were prognostic factors for CESC (P<0.05). The target genes were identified to be involved in the MAPK, cAMP, PI3K/Akt, mTOR, and estrogen cancer signaling pathways. Protein–protein interaction network analysis showed that TP53, MMP1, NOTCH1, SMAD4, and NFKB1 were the most significant hub proteins. Conclusion Our results indicate that hsa-miR-508, hsa-miR-509-2, and hsa-miR-526b may be potential diagnostic biomarkers for early-stage CESC with LNM, and serve as prognostic predictors for patients with CESC. Keywords: lymph nodes, miRNAs, prognosis, uterine cervical neoplasms Introduction Cervical cancer (CESC) is one of the leading etiologies of malignant tumor deaths in woman aged 20–39 years. According to the National Center for Health Statistics, in 2017, there were 12,820 estimated new cases of CESC and 4,210 estimated deaths caused by CESC in America.[31]^1 In fact, CESC deaths in developing countries account for 90% of the world’s total death toll, causing severe public health problems and financial burden.[32]^2 Some studies have shown that miRNAs are associated with the development and prognosis of CESC.[33]^3^–[34]^5 Previous studies have reported that lymph node metastasis (LNM) in early-stage CESC influences the patient’s survival and outcomes of subsequent adjuvant therapy, but there is limited research on the relationship between LNM and miRNAs in early-stage CESC.[35]^6^–[36]^9 Therefore, it is meaningful to study the molecular mechanisms of CESC and identify miRNA markers for individualized treatment plans and prognostic prediction of CESC. miRNAs are small noncoding RNAs of ~22 nucleotides in length that play regulatory roles in RNA silencing.[37]^10 miRNAs are involved in regulating 30% of human genes.[38]^11 Through the regulation of target genes, miRNAs play a role in cell cycle regulation, cell proliferation, apoptosis, drug resistance, and other biological processes (BPs).[39]^12^–[40]^17 Additionally, miRNAs are stable in tissue, blood, and urine.[41]^18^–[42]^20 Therefore, miRNAs show promise for use in the diagnosis and treatment of CESC. The Cancer Genome Atlas (TCGA) contains genomic and proteomic data from more than 20 tumor types, with the benefits of a large sample size, a unified testing platform, and complete follow-up information; moreover, an increasing number of studies utilize data from this database.[43]^3^,[44]^5^,[45]^21^–[46]^23 High-throughput dataset analysis may provide a useful method for the identification of hub genes and patterns in CESC. In this study, we analyzed the differences in miRNA expression patterns between early-stage CESC patients with LNM [LNM(+)] and without LNM [LNM(−)]. Additionally, we analyzed the prognostic value of these miRNAs (hsa-miR-508, hsa-miR-509-2, and hsa-miR-526b) and their possible mechanisms of action, which may offer new insight for improving CESC treatment and patient survival. Based on our results, hsa-miR-508, hsa-miR-509-2, and hsa-miR-526b may be potential diagnostic biomarkers and prognostic indicators of CESC. Materials and methods CESC patients in TCGA The information and miRNA data for 307 patients with CESC were obtained from the TCGA database ([47]https://cancergenome.nih.gov/). The exclusion criteria were as follows: 1) information on LNM was not provided or not applicable; 2) the clinical stage was not stage IA2–IIA (International Federation of Gynecology and Obstetrics 2009 criteria); and 3) the clinical data were incomplete. A total of 145 patients with early-stage CESC in the training group were used to analyze differentially expressed genes (DEGs). In addition, 275 CESC patients in the validation group were used to evaluate the prognosis of CESC. miRNA expression data analyses The miRNA expression analyses were performed using the “edgeR” package in R.[48]^24 From a total of 1,046 miRNAs, 422 miRNAs were detected in at least two cohorts, and miRNAs with no ID on miRBase21 or with a mean expression <3 raw counts were removed. Each individual miRNA was log2-transformed for further analysis. All DEGs were considered significant if the absolute fold changes were more than 1.5 and the P-value was <0.05. Prediction of miRNA target genes To date, miRWalk2.0 is the only comprehensive archive that can be accessed for free, providing the largest available predictive and experimental miRNA–target interaction data and offering novel and unique features.[49]^25 Target genes of three miRNAs were predicted based on the miRWalk database, and overlapped genes among the miRanda, PITA, and Targetscan databases were selected. A P-value<0.05 represented statistical significance. GO, KEGG, and PPI network analyses Gene ontology (GO) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses are the common methods for studying the functions and pathways associated with large-scale genome data or transcriptional data. GO and KEGG enrichment analyses were performed using BiNGO 3.0.3 and DAVID 6.8.[50]^26^,[51]^27 The protein–protein interaction (PPI) networks were analyzed through the STRING 10.5 database and Molecular Complex Detection (MCODE) 1.5.1.[52]^28^,[53]^29 The network was further visualized using Cytoscape 3.6.1.[54]^30 Statistical analyses The data are shown as the mean ± SD or as the median. Chi-square tests were used to assess differences in race, tumor grade, and tumor status between the LNM(+) and LNM(−) cohorts. Other data, such as age and smoking history, were evaluated using independent-sample t-tests. Kaplan–Meier survival analyses and univariate and multivariate Cox proportional risk regression analyses were conducted to compare miRNA levels (high- and low-expression levels) and prognosis. The data were analyzed using GraphPad Prism 5.0 (La Jolla, CA, USA) and R software, and a P-value <0.05 was considered statistically significant. Results The effect of lymph node status on the survival of CESC patients This study involved 145 patients with early-stage CESC, including 32 LNM(+) patients and 113 LNM(−) patients. The clinical data of these patients are presented in [55]Table 1. Indeed, there were no significant differences in age, race, smoking history, or tumor grade between the LNM(+) and LNM(−) patients. Survival analysis showed that the survival time of the LNM(+) group was shorter than that of the LNM(−) group (P<0.05; [56]Figure 1). Table 1. Clinical characteristics of patients with early-stage CESC Category Cohort LNM+ (n=32) Cohort LNM− (n=113) P-value Age (years) 0.580  Mean ± SD 47.13±13.14 46.04±12.89  Median 46 45 Race 0.122  Asian 2 10  White 21 86  Black or African American 7 7  Native Hawaiian or other Pacific islander 0 1  American Indian or Alaska native 0 2 Smoking history 0.422  Mean ± SD 1.75±1.29 1.54±1.22  Median 2 1 Tumor grade 0.911  G1 + G2 17 61  G3 + G4 14 48  GX 0 3  NA 1 1 Tumor status 0.001  With tumor 11 12  Tumor free 18 88  NA 0 4  Unknown 3 9 [57]Open in a new tab Abbreviations: CESC, cervical cancer; LNM, lymph node metastasis; NA, not available. Figure 1. Figure 1 [58]Open in a new tab Survival analysis of early-stage LNM(+) and LNM(−) cervical cancer patients. Note: Kaplan–Meier survival curves show that LNM(+) was associated with poor prognosis in early-stage CESC patients. Abbreviations: CESC, cervical cancer; LNM, lymph node metastasis. Differentially expressed miRNAs in LNM(+) and LNM(−) patients Comprehensive miRNA profiles of early-stage CESC were generated for LNM(+) and LNM(−) patients. In total, 75 miR-NAs were differentially regulated and had a P-value <0.05 and an absolute fold change of more than 1.5 ([59]Supplementary table S1). As shown in [60]Figure 2, of these 75 miRNAs, 24 miRNAs (32.0%) were downregulated in LNM(+) patients, while the remaining 51 (68.0%) miRNAs were upregulated in LNM(+) patients. Seventy-five miRNAs were selected for further analyses. Figure 2. [61]Figure 2 [62]Open in a new tab Differentially expressed miRNAs in LNM(+) patients. Notes: Volcano plots of miRNAs were constructed using fold changes and P-values. The vertical lines correspond to a 1.5-fold change in expression (up or down), and the horizontal line represents P=0.05. The green points on the plot represent the 24 differentially expressed miRNAs that were downregulated in LNM(+) patients, and the red points on the plot represent the 51 differentially expressed miRNAs upregulated in LNM(+) patients. Abbreviation: LNM, lymph node metastasis. Identification of miRNAs correlated with the survival of CESC patients To identify the association between miRNAs and the survival of CESC patients, the overall survival times of 275 CESC patients in the TCGA dataset were analyzed by the Kaplan–Meier method and log-rank tests. The results showed that patients with high hsa-miR-508 or hsa-miR-526b levels had a longer overall survival time than did those with low hsa-miR-508 level or hsa-miR-526b levels (log-rank test: P=0.006, 0.003), as shown in [63]Figure 3. The correlation of hsa-miR-509-2 level with overall survival time in the 275 CESC patients was not significant. The miRNA expression data are shown in [64]Figure 4 and [65]Table 2. Figure 3. [66]Figure 3 [67]Open in a new tab Kaplan–Meier survival curves of the expression of the three miRNAs of interest in CESC. Notes: Kaplan–Meier survival curves show that low miRNA expression was associated with poor survival in CESC patients. Patients were divided into high and low levels according to the median of each miRNA. (A) hsa-miR-508; (B) hsa-miR-509-2; (C) hsa-miR-526b. Abbreviation: CESC, cervical cancer. Figure 4. [68]Figure 4 [69]Open in a new tab Differential expression of the three miRNAs of interest in LNM(+) and LNM(−) patients. Notes: (A) The expression levels of hsa-miR-508 in LNM(+) and LNM(−) patients. (B) The expression levels of hsa-miR-509-2 in LNM(+) and LNM(−) patients. (C) The expression levels of hsa-miR-526b in LNM(+) and LNM(−) patients. The expression level values were calculated using log2-transformed values. Abbreviation: LNM, lymph node metastasis. Table 2. Three downregulated miRNAs in early-stage CESC patients with lymph node metastases No. ID log2(FC) P-value FDR Regulation 1 hsa-mir-508 −1.57 2.92×10^−4 4.56×10^−3 Down 2 hsa-mir-509-2 −1.74 1.47×10^−3 1.62×10^−2 Down 3 hsa-mir-526b −1.67 8.89×10^−3 5.89×10^−2 Down [70]Open in a new tab Abbreviations: FC, fold change; FDR, false discovery rate ID, miRNA identification. Univariate and multivariate Cox regression analyses in CESC patients The clinical features of age, smoking history, clinical stage, tumor grade, and metastasis were considered in the univariate and multivariate Cox regression analyses to analyze the prognosis of CESC patients. The univariate analysis suggested that there was a significant association between overall survival time and clinical stage (HR: 2.329, 95% CI: 1.423–3.812, P=0.001), hsa-miR-508 (HR: 0.511, 95% CI: 0.315–0.828, P=0.006), and hsa-miR-526b (HR: 0.484, 95% CI: 0.294–0.795, P=0.004). No significant association was observed between overall survival time and age, smoking history, tumor grade, or metastasis. The multivariate analysis suggested that clinical stage (HR: 2.368, 95% CI: 1.414–3.967, P=0.001), hsa-miR-508 (HR: 0.496, 95% CI: 0.298–0.826, P=0.007), hsa-miR-509-2 (HR: 0.596, 95% CI: 0.363–0.980, P=0.041), and hsa-miR-526b (HR: 0.591, 95% CI: 0.350–0.996, P=0.048) were prognostic factors for CESC ([71]Table 3). Table 3. Univariate and multivariate Cox regression analyses in CESC patients Category Univariable HR (95% CI) P-value Multivariable HR (95% CI) P-value Age (>60 years vs <60 years) 1.688 (0.994–2.799) 0.053 Smoking history (>3 years vs <3 years) 0.582 (0.276–1.226) 0.154 Clinical stage (III + IV vs I + II) 2.329 (1.423–3.812) 0.001 2.368 (1.414–3.967) 0.001 Tumor grade (G3 + G4 vs G1 + G2) 0.798 (0.482–1.320) 0.379 Metastasis (M1 vs M0) 2.450 (0.889–6.749) 0.083 hsa-miR-508 (high vs low) 0.511 (0.315–0.828) 0.006 0.496 (0.298–0.826) 0.007 hsa-miR-509-2 (high vs low) 0.655 (0.409–1.050) 0.079 0.596 (0.363–0.980) 0.041 hsa-miR-526b (high vs low) 0.484 (0.294–0.795) 0.004 0.591 (0.350–0.996) 0.048 [72]Open in a new tab Abbreviation: CESC, cervical cancer. Functional implication of prognostic miRNA signatures To further analyze the mechanisms of these miRNAs (hsa-miR-508, hsa-miR-509-2, and hsa-miR-526b) in CESC, we used miRWalk2.0, which identifies overlapping genes between miRanda, PITA, and the Targetscan databases to predict target genes. Then, we performed functional enrichment analysis with GO and KEGG using BiNGO and DAVID. The GO BP terms were remarkably enriched in epithelial cell development, regulation of cell shape, neuronal migration, and so on. For the GO cellular component (CC) terms, the target genes were concentrated in the cortical actin cytoskeleton, postsynaptic density, and so on. The molecular function (MF) category was focused on steroid hormone receptor activity, NOTCH binding, and so on ([73]Figure 5, [74]Table 4). In addition, the KEGG pathways were mostly enriched in cancer pathways, including the MAPK, cAMP, PI3K/Akt, mTOR, and estrogen signaling pathways ([75]Figure 6, [76]Table 5). The GO network was illustrated as an interaction network using the Cytoscape plug-in BiNGO, which maps the main functional categories of a given gene set in the GO hierarchy. As shown in [77]Figure 7, the regulation of cellular and primary metabolic processes was important in the BP network. In the CC network, the cell and intracellular components were in core positions ([78]Figure 8). Additionally, the ion, cation, metal ion, and zinc ion binding components were significant in the MF network ([79]Figure 9). Figure 5. [80]Figure 5 [81]Open in a new tab GO terms for the target genes of three miRNAs using DAVID. Abbreviation: GO, gene ontology. Table 4. GO terms enriched in target genes GO ID GO term Count P-value Genes miRNA BP GO:0002064 Epithelial cell development 8 1.52×10^−5 B4GALT1, SHROOM3, HYDIN, ONECUT1, SLC9A4, ONECUT2, TP63, SLC4A5 hsa-miR-508 hsa-miR-509-2 hsa-miR-526b GO:0008360 Regulation of cell shape 28 2.42×10^−4 GNA13, DMTN, SHROOM3, WASF3, HEXB, FERMT2, ATP10A, RHOQ, etc hsa-miR-508 hsa-miR-509-2 hsa-miR-526b GO:0001764 Neuron migration 22 6.93×10^−4 DCC, TUBB2B, SOX1, NDN, ASTN1, FZD3, SPOCK1, NTN1, etc hsa-miR-508 hsa-miR-509-2 hsa-miR-526b GO:0008654 Phospholipid biosynthetic process 12 8.80×10^−4 PPARD, AGPAT5, LPGAT1, SH3GLB1, SERAC1, HEXB, MBOAT2, PISD, etc hsa-miR-508 hsa-miR-509-2 hsa-miR-526b GO:0030335 Positive regulation of cell migration 31 2.05×10^−3 CGA, ONECUT1, WASH1, ONECUT2, LRRC15, SEMA5A, ITGAV, SEMA3F, etc hsa-miR-508 hsa-miR-509-2 hsa-miR-526b GO:0043401 Steroid hormone- mediated signaling pathway 14 2.06×10^−3 PPARD, THRB, ESRRB, RXRA, PAQR8, RORB, NR2C2, NR1H2, etc hsa-miR-508 hsa-miR-509-2 hsa-miR-526b GO:0030522 Intracellular receptor signaling pathway 11 2.15×10^−3 NOTCH2, PPARD, NCOA2, THRB, NR1D2, ESRRB, RORB, NR2F2, etc hsa-miR-508 hsa-miR-509-2 hsa-miR-526b GO:0007010 Cytoskeleton organization 27 4.47×10^−3 DMTN, TUBB2B, WASF3, ABLIM3, NEDD9, RHOU, DMD, STRIP2, etc hsa-miR-508 hsa-miR-509-2 hsa-miR-526b CC GO:0030864 Cortical actin cytoskeleton 12 1.46×10^−3 EEF1A1, PPP1R9A, SHROOM2, NF2, MED28, SHROOM4, LASP1, FLOT2, etc hsa-miR-508 hsa-miR-509-2 hsa-miR-526b GO:0014069 Postsynaptic density 31 1.57×10^−3 NETO2, DMTN, LZTS1, SYNDIG1, CNN3, CPEB3, BMPR2, SPOCK1, etc hsa-miR-508 hsa-miR-509-2 hsa-miR-526b GO:0055038 Recycling endosome membrane 12 1.79×10^−3 SCAMP1, RAP2B, RAB11FIP4, RAP2A, RAB11FIP2, RAP2C, WASH1, SLC26A7, etc hsa-miR-508 hsa-miR-509-2 hsa-miR-526b GO:0000932 Cytoplasmic mRNA processing body 17 1.83×10^−3 NANOS2, CNOT2, BTBD6, APOBEC3H, CNOT7, APOBEC3F, CAPRIN1, PATL1, etc hsa-miR-508 hsa-miR-509-2 hsa-miR-526b GO:0043235 Receptor complex 23 2.93×10^−3 EGFR, RNMT, IL23R, OLR1, LDLR, RXRA, TGFBR1, LIFR, SMAD3, etc hsa-miR-508 hsa-miR-509-2 hsa-miR-526b GO:0000139 Golgi membrane 75 4.42×10^−3 SLC9A8, OSBP, SEC24A, NDST1, VAPB, FAM20B, CNGB1, RHOU, etc hsa-miR-508 hsa-miR-509-2 hsa-miR-526b GO:0031901 Early endosome membrane 21 4.68×10^−3 EGFR, MMGT1, SYNDIG1, WASH1, PML, EEA1, CFTR, WLS, etc hsa-miR-508 hsa-miR-509-2 hsa-miR-526b GO:0005884 Actin filament 14 5.93×10^−3 MYO5A, DMTN, MYO1B, RHOQ, ACKR2, ACTN2, MYO9B, GNG12, etc hsa-miR-508 hsa-miR-509-2 hsa-miR-526b GO:0015629 Actin cytoskeleton 33 5.99×10^−3 MTSS1, DMTN, CNN3, SHROOM4, ONECUT2, KNTC1, PDLIM2, ARHGAP35, etc hsa-miR-508 hsa-miR-509-2 hsa-miR-526b MF GO:0003707 Steroid hormone receptor activity 14 1.82×10^−3 PPARD, THRB, ESRRB, RXRA, PAQR8, RORB, NR2C2, NR1H2, etc hsa-miR-508 hsa-miR-509-2 hsa-miR-526b GO:0005254 Chloride channel activity 13 3.95×10^−3 GABRG1, SLC1A4, GABRG2, GABRG3, GABRA4, CLIC4, SLC26A7, CLIC5, etc hsa-miR-508 hsa-miR-509-2 hsa-miR-526b GO:0005112 Notch binding 7 4.85×10^−3 NOV, NOTCH1, HIF1AN, DNER, AAK1, JAG1, GALNT11 hsa-miR-508 hsa-miR-509-2 hsa-miR-526b GO:0003779 Actin binding 41 4.86×10^−3 MYO5A, LRRC10, DMTN, SHROOM2, CNN3, WASH1, WASF3, IMPACT, etc hsa-miR-508 hsa-miR-509-2 hsa-miR-526b GO:0004879 RNA polymerase II transcription factor activity, ligand-activated sequence-specific DNA binding 10 5.30×10^−3 NR1H2, PPARD, NR1D2, ESRRB, RXRA, RORB, NR2F2, NR1H4, etc hsa-miR-508 hsa-miR-509-2 hsa-miR-526b GO:0000980 RNA polymerase II distal enhancer sequence-specific DNA binding 14 7.12×10–3 CDX2, ESRRB, NFKB1, MBD2, GZF1, NONO, ARX, FLI1, etc hsa-miR-508 hsa-miR-509-2 hsa-miR-526b [82]Open in a new tab Abbreviations: BP, biological process; CC, cellular component; GO, gene ontology; MF, molecular function. Figure 6. [83]Figure 6 [84]Open in a new tab KEGG pathways for the target genes of three miRNAs using DAVID. Abbreviations: HCM, hypertrophic cardiomyopathy; KEGG, Kyoto Encyclopedia of Genes and Genomes. Table 5. KEGG pathways enriched in target genes KEGG ID KEGG term Count P-value Genes miRNA hsa05200 Pathways in cancer 61 7.27×10^−5 GNA13, F2RL3, ADCY1, PPARD, FGF7, ADCY7, PTGS2, STAT5B, etc hsa-miR-508 hsa-miR-509-2 hsa-miR-526b hsa04010 MAPK signaling pathway 37 7.52×10^−3 FGF7, PPP3R1, PPM1A, PPP3R2, NFKB1, GNG12, TNFRSF1A, KRAS, etc hsa-miR-508 hsa-miR-509-2 hsa-miR-526b hsa04024 cAMP signaling pathway 29 1.68×10^−2 ADCY1, ADCY7, ATP1B4, OXTR, NFKB1, GRIN3B, GABBR2, CNGB1, etc hsa-miR-508 hsa-miR-509-2 hsa-miR-526b hsa04068 FoxO signaling pathway 21 2.35×10^−2 EGFR, IRS2, SGK1, RBL2, NLK, TGFBR1, SMAD4, SMAD3 etc. hsa-miR-508 hsa-miR-509-2 hsa-miR-526b hsa04151 PI3K–Akt signaling pathway 44 2.94×10^−2 CSH2, FGF7, PHLPP2, TLR2, NFKB1, GNG12, IL7R, PTEN, etc hsa-miR-508 hsa-miR-509-2 hsa-miR-526b hsa04150 mTOR signaling pathway 11 4.07×10^−2 PRKCA, MAPK1, RPS6KA3, TSC1, IGF1, CAB39, RICTOR, PRKAA2, etc hsa-miR-508 hsa-miR-509-2 hsa-miR-526b [85]Open in a new tab Abbreviation: KEGG, Kyoto Encyclopedia of Genes and Genomes. Figure 7. [86]Figure 7 [87]Open in a new tab The BP network showed as an interaction network using the Cytoscape plug-in BiNGO. Notes: Color depth represented the degree of enrichment of GO terms. The significance of enrichment is represented by yellow color (P<0.05). Abbreviations: BP, biological process; GO, gene ontology. Figure 8. [88]Figure 8 [89]Open in a new tab The CC network showed as an interaction network using the Cytoscape plug-in BiNGO. Notes: Color depth represented the degree of enrichment of GO terms. The significance of enrichment is represented by yellow color (P<0.05). Abbreviations: CC, cellular component; GO, gene ontology. Figure 9. [90]Figure 9 [91]Open in a new tab The MF network showed as an interaction network using the Cytoscape plug-in BiNGO. Notes: Color depth represented the degree of enrichment of GO terms. The significance of enrichment is represented by yellow color (P<0.05). Abbreviations: GO, gene ontology; MF, molecular function. Selection of hub proteins The PPI network of target genes was analyzed according to the STRING database, including 956 nodes and 3,448 edges. Then, we used the MCODE plug-in in Cytoscape software to build the core module, which consist of 29 nodes and 95 edges (MCODE score: 7.000). The highly connected node with the highest height is considered the hub protein of the network. A total of 29 proteins were selected as the most significant hub proteins, including TP53 (a target of hsa-miR-508, hsa-miR-509-2, and hsa-miR-526b), MMP1 (a target of hsa-miR-509-2 and hsa-miR-526b), NOTCH1 (a target of hsa-miR-509-2), SMAD4 (a target of hsa-miR-508 and hsa-miR-526b), and NFKB1 (a target of hsa-miR-508), which might be associated with LNM in early-stage CESC ([92]Figure 10). Figure 10. Figure 10 [93]Open in a new tab PPI network of hub proteins. Notes: The MCODE was conducted to screen modules of the PPI network with a degree cutoff =2, node score cutoff =0.2, k-core =2, and maximum depth =100. Abbreviation: PPI, protein–protein interaction. Discussion CESC is one of the major primary gynecological tumors that affect women’s health. In developing countries, the 5-year overall survival rate is only 46%, which is far lower than that of developed countries.[94]^2 LNM in the early stage affects the prognosis of CESC.[95]^6 Therefore, it is necessary to identify prognostic miRNAs and determine their molecular mechanisms. In this study, we found that hsa-miR-508, hsa-miR-509-2, and hsa-miR-526b were downregulated in LNM(+) patients and could effectively evaluate the prognosis of CESC. Moreover, we explained the potential mechanisms via GO and KEGG enrichment and PPI network analyses of target genes. As far as we know, there is limited research on the role of these three miRNAs (hsa-miR-508, hsa-miR-509-2, and hsa-miR-526b) in CESC. We first demonstrated that three miRNAs were downregulated in early-stage CESC LNM(+) patients and influenced the prognosis of CESC through bio-informatics analysis. In a recent study, Liu et al reported that hsa-miR-508 is an independent, protective prognostic factor of glioma, and functional tests confirmed that overexpression of hsa-miR-508 inhibits the migration of glioma cells, which was similar to the results of our study.[96]^31 In gastric cancer, hsa-miR-508 inhibits the invasion of gastric cancer cells by targeting S-phase kinase-associated protein 2.[97]^32 Huang et al also found that hsa-miR-508 expression is downregulated in gastric carcinogenesis cell lines and plays a major role as a tumor suppressor gene in gastric cancer.[98]^33 Moreover, miRNA microarray profiling showed that hsa-miR-508 expression is downregulated in lymph node metastatic tissues.[99]^34 In addition, Zhai et al indicated that hsa-miR-508 inhibits renal cell cancer cell invasion and migration.[100]^35 hsa-miR-509 serves as a cancer suppressor gene in various cancers, such as lung cancer, breast cancer, renal cell cancer, and pancreatic cancer.[101]^36^–[102]^40 hsa-miR-509 can inhibit tumor cell invasion and migration through target genes. Using bioinformatics analysis, hsa-miR-509 was found to be downregulated in breast cancer and brain metastasis compared to primary cancer and was shown to play a role in suppressing brain metastasis by regulating the RhoC-TNF-α network in both in vivo and in vitro experiments.[103]^41 Additionally, low expression levels of hsa-miR-509 were associated with poor outcomes in pancreatic cancer patients.[104]^39 hsa-miR-526b can inhibit tumor progression by targeting MMP1.[105]^42 Meanwhile, Zhang et al found that hsa-miR-526b is downregulated in lung cancer, and low expression was associated with poor clinical outcomes through the p53/p21 pathway.[106]^43 In colon cancer, hsa-miR-526b was reported to be downregulated in metastatic tissue and to inhibit cell proliferation and invasion by regulating the hub gene HIF-1α.[107]^44 Liu et al observed that low hsa-miR-526b expression correlates with poor prognosis in hepatocellular carcinoma and enhanced lung metastasis in mouse lungs via the ERK pathway.[108]^45 To investigate the mechanism of LNM in early-stage CESC, we conducted KEGG pathway and GO annotation enrichment analyses and PPI network analyses on three miRNA target genes. Pathway enrichment analysis showed that the target genes were mainly enriched in pathways in cancer, such as MAPK, cAMP, PI3K/Akt, and mTOR pathways. The MAPK pathway has been reported to play a role in the tumor progression of CESC.[109]^13^,[110]^46^,[111]^47 Zhang et al found that the MAPK pathway was associated with LNM in head and neck cancer.[112]^48 Research has shown that the Akt/ mTOR pathway is activated in medullary thyroid carcinoma with lymph node metastases.[113]^49 In the coexpression regulatory network of target genes, multiple studies have shown that a TP53 gene mutation increases the risk of CESC.[114]^50^–[115]^52 MMP1 has been associated with LNM in esophageal squamous cell cancer,[116]^53 lung cancer,[117]^54 and gastric cancer.[118]^55 Liu et al reported that MMP1 promotes esophageal squamous cell carcinoma cell migration and invasion via the PI3K/ Akt pathway.[119]^53 In CESC, NOTCH1 has been associated with LNM, and high expression of NOTCH1 indicates a poor prognosis in patients.[120]^56 At the same time, Park et al found that NOTCH1 could be a molecular marker of LNM in papillary thyroid cancer.[121]^57 In pancreatic and colorectal cancer, SMAD4 was found to be associated with LNM and survival.[122]^58^,[123]^59 Another gene in the network, NFKB1, was reported to increase the risk of LNM in gastric cancer and oral cancer.[124]^60^,[125]^61 These results suggest that the target genes of these three miRNAs play a significant role in LNM in early-stage CESC. However, there are some deficiencies in our research at present. Although the three miRNAs were downregulated in LNM(+) patients according to the bioinformatics analyses, the original samples need additional validation. Functional tests need to be conducted, and the mechanism of LNM in early-stage CESC needs further study. Our group is currently pursuing some of these research goals. Conclusion Using bioinformatics analysis, we have shown for the first time that three miRNAs were downregulated in early-stage CESC LNM(+) patients and are potential prognostic indicators for CESC. In addition, five central genes (TP53, MMP1, NOTCH1, SMAD4, and NFKB1) were found to be directly or indirectly involved in LNM in early-stage CESC. Taken together, the results of the present study suggest that hsa-miR-508, hsa-miR-509-2, and hsa-miR-526b may be potential new diagnostic and prognostic markers for CESC in the future. Supplementary material Table S1. Differentially expressed microRNAs in LNM(+) and LNM(−) patients No. ID log2(FC) P-value FDR Regulation in LNM(+) 1 hsa-mir-599 −5.611 2.83E–05 8.53E–04 Down 2 hsa-mir-499 −4.191 4.60E–06 3.08E–04 Down 3 hsa-mir-767 −2.967 1.43E–03 1.62E–02 Down 4 hsa-mir-551b −2.847 5.13E–06 3.08E–04 Down 5 hsa-mir-105-2 −2.463 1.38E–02 8.09E–02 Down 6 hsa-mir-105-1 −2.453 1.09E–02 6.66E–02 Down 7 hsa-mir-885 −2.388 3.05E–04 4.59E–03 Down 8 hsa-mir-891a −1.962 1.06E–02 6.54E–02 Down 9 hsa-mir-449a −1.897 4.71E–03 4.16E–02 Down 10 hsa-mir-520a −1.869 5.67E–03 4.87E–02 Down 11 hsa-mir-509-1 −1.801 4.55E–04 6.18E–03 Down 12 hsa-mir-506 −1.798 1.23E–02 7.42E–02 Down 13 hsa-mir-509-2 −1.741 1.47E–03 1.62E–02 Down 14 hsa-mir-509-3 −1.697 1.08E–03 1.34E–02 Down 15 hsa-mir-526b −1.666 8.89E–03 5.89E–02 Down 16 hsa-mir-508 −1.569 2.92E–04 4.56E–03 Down 17 hsa-mir-514-3 −1.526 3.19E–03 2.99E–02 Down 18 hsa-mir-514-2 −1.358 8.12E–03 5.78E–02 Down 19 hsa-mir-514-1 −1.331 9.10E–03 5.89E–02 Down 20 hsa-mir-1248 −0.902 1.60E–02 9.13E–02 Down 21 hsa-mir-181b-2 −0.811 8.98E–03 5.89E–02 Down 22 hsa-mir-3651 −0.721 4.59E–02 2.33E–01 Down 23 hsa-mir-582 −0.714 4.90E–02 2.34E–01 Down 24 hsa-mir-425 −0.694 1.43E–03 1.62E–02 Down 25 hsa-mir-548v 0.615 4.64E–02 2.33E–01 Up 26 hsa-mir-3677 0.628 9.46E–03 6.03E–02 Up 27 hsa-mir-1228 0.632 6.25E–03 5.06E–02 Up 28 hsa-mir-889 0.633 1.57E–02 9.03E–02 Up 29 hsa-mir-379 0.635 8.37E–03 5.78E–02 Up 30 hsa-mir-1274b 0.635 3.30E–02 1.78E–01 Up 31 hsa-mir-2277 0.644 8.34E–03 5.78E–02 Up 32 hsa-mir-337 0.666 8.26E–03 5.78E–02 Up 33 hsa-mir-299 0.714 6.80E–03 5.30E–02 Up 34 hsa-mir-382 0.735 5.90E–03 4.87E–02 Up 35 hsa-mir-138-2 0.751 3.35E–02 1.78E–01 Up 36 hsa-mir-409 0.761 2.11E–03 2.12E–02 Up 37 hsa-mir-411 0.794 1.95E–03 2.00E–02 Up 38 hsa-mir-127 0.815 4.16E–04 5.84E–03 Up 39 hsa-mir-543 0.877 1.02E–02 6.39E–02 Up 40 hsa-mir-654 0.881 1.89E–03 1.99E–02 Up 41 hsa-mir-3127 0.897 2.43E–05 8.52E–04 Up 42 hsa-mir-539 0.917 1.43E–03 1.62E–02 Up 43 hsa-mir-149 0.919 2.61E–04 4.39E–03 Up 44 hsa-mir-412 0.920 1.58E–03 1.71E–02 Up 45 hsa-mir-487a 0.920 7.78E–03 5.78E–02 Up 46 hsa-mir-129-1 0.927 7.72E–03 5.78E–02 Up 47 hsa-mir-129-2 0.940 3.19E–03 2.99E–02 Up 48 hsa-mir-487b 0.946 4.82E–04 6.35E–03 Up 49 hsa-mir-1266 0.956 8.46E–04 1.08E–02 Up 50 hsa-mir-494 0.980 2.67E–03 2.61E–02 Up 51 hsa-mir-758 1.016 3.54E–04 5.13E–03 Up 52 hsa-mir-615 1.020 4.74E–03 4.16E–02 Up 53 hsa-mir-154 1.029 1.55E–04 2.71E–03 Up 54 hsa-mir-369 1.036 5.57E–05 1.17E–03 Up 55 hsa-mir-548b 1.070 2.92E–04 4.56E–03 Up 56 hsa-mir-655 1.071 1.33E–04 2.55E–03 Up 57 hsa-mir-370 1.088 1.52E–04 2.71E–03 Up 58 hsa-mir-376c 1.141 2.82E–05 8.53E–04 Up 59 hsa-mir-493 1.164 3.47E–05 8.63E–04 Up 60 hsa-mir-134 1.167 1.26E–05 5.81E–04 Up 61 hsa-mir-377 1.180 3.48E–05 8.63E–04 Up 62 hsa-mir-410 1.183 3.11E–05 8.63E–04 Up 63 hsa-mir-495 1.221 6.34E–06 3.33E–04 Up 64 hsa-mir-137 1.245 8.76E–03 5.89E–02 Up 65 hsa-mir-432 1.250 5.50E–05 1.17E–03 Up 66 hsa-mir-937 1.256 1.38E–05 5.81E–04 Up 67 hsa-mir-219-2 1.262 7.88E–03 5.78E–02 Up 68 hsa-mir-485 1.278 4.86E–05 1.14E–03 Up 69 hsa-mir-433 1.336 1.07E–04 2.14E–03 Up 70 hsa-mir-376a-1 1.432 2.34E–05 8.52E–04 Up 71 hsa-mir-376b 1.552 4.70E–07 4.94E–05 Up 72 hsa-mir-380 1.691 2.57E–06 2.17E–04 Up 73 hsa-mir-431 1.739 2.74E–08 3.85E–06 Up 74 hsa-mir-323 2.175 5.30E–09 2.23E–06 Up 75 hsa-mir-376a-2 2.222 1.26E–08 2.64E–06 Up [126]Open in a new tab Abbreviations: ID, miRNA ID; FC, fold change; FDR, false discovery rate. Acknowledgments