Identification of differentially expressed miRNAs in early-stage cervical cancer with lymph node metastasis across The Cancer Genome Atlas datasets.

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.

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.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.2 Some studies have shown that miRNAs are associated with the development and prognosis of CESC.3–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.6–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.10

miRNAs are involved in regulating 30% of human genes.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).12–17 Additionally, miRNAs are stable in tissue, blood, and urine.18–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.3,5,21–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 (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.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.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.26,27 The protein–protein interaction (PPI) networks were analyzed through the STRING 10.5 database and Molecular Complex Detection (MCODE) 1.5.1.28,29 The network was further visualized using Cytoscape 3.6.1.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 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; 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

Abbreviations: CESC, cervical cancer; LNM, lymph node metastasis; NA, not available.

Figure 1

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 (Supplementary table S1). As shown in 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

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 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 Figure 4 and Table 2.

Figure 3

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

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

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 (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

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 (Figure 5, Table 4). In addition, the KEGG pathways were mostly enriched in cancer pathways, including the MAPK, cAMP, PI3K/Akt, mTOR, and estrogen signaling pathways (Figure 6, 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 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 (Figure 8). Additionally, the ion, cation, metal ion, and zinc ion binding components were significant in the MF network (Figure 9).

Figure 5

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-508hsa-miR-509-2hsa-miR-526b
GO:0008360	Regulation of cell shape	28	2.42×10−4	GNA13, DMTN, SHROOM3, WASF3, HEXB, FERMT2, ATP10A, RHOQ, etc	hsa-miR-508hsa-miR-509-2hsa-miR-526b
GO:0001764	Neuron migration	22	6.93×10−4	DCC, TUBB2B, SOX1, NDN, ASTN1, FZD3, SPOCK1, NTN1, etc	hsa-miR-508hsa-miR-509-2hsa-miR-526b
GO:0008654	Phospholipid biosynthetic process	12	8.80×10−4	PPARD, AGPAT5, LPGAT1, SH3GLB1, SERAC1, HEXB, MBOAT2, PISD, etc	hsa-miR-508hsa-miR-509-2hsa-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-508hsa-miR-509-2hsa-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-508hsa-miR-509-2hsa-miR-526b
GO:0030522	Intracellular receptor signaling pathway	11	2.15×10−3	NOTCH2, PPARD, NCOA2, THRB, NR1D2, ESRRB, RORB, NR2F2, etc	hsa-miR-508hsa-miR-509-2hsa-miR-526b
GO:0007010	Cytoskeleton organization	27	4.47×10−3	DMTN, TUBB2B, WASF3, ABLIM3, NEDD9, RHOU, DMD, STRIP2, etc	hsa-miR-508hsa-miR-509-2hsa-miR-526b
CC
GO:0030864	Cortical actin cytoskeleton	12	1.46×10−3	EEF1A1, PPP1R9A, SHROOM2, NF2, MED28, SHROOM4, LASP1, FLOT2, etc	hsa-miR-508hsa-miR-509-2hsa-miR-526b
GO:0014069	Postsynaptic density	31	1.57×10−3	NETO2, DMTN, LZTS1, SYNDIG1, CNN3, CPEB3, BMPR2, SPOCK1, etc	hsa-miR-508hsa-miR-509-2hsa-miR-526b
GO:0055038	Recycling endosome membrane	12	1.79×10−3	SCAMP1, RAP2B, RAB11FIP4, RAP2A, RAB11FIP2, RAP2C, WASH1, SLC26A7, etc	hsa-miR-508hsa-miR-509-2hsa-miR-526b
GO:0000932	Cytoplasmic mRNA processing body	17	1.83×10−3	NANOS2, CNOT2, BTBD6, APOBEC3H, CNOT7, APOBEC3F, CAPRIN1, PATL1, etc	hsa-miR-508hsa-miR-509-2hsa-miR-526b
GO:0043235	Receptor complex	23	2.93×10−3	EGFR, RNMT, IL23R, OLR1, LDLR, RXRA, TGFBR1, LIFR, SMAD3, etc	hsa-miR-508hsa-miR-509-2hsa-miR-526b
GO:0000139	Golgi membrane	75	4.42×10−3	SLC9A8, OSBP, SEC24A, NDST1, VAPB, FAM20B, CNGB1, RHOU, etc	hsa-miR-508hsa-miR-509-2hsa-miR-526b
GO:0031901	Early endosome membrane	21	4.68×10−3	EGFR, MMGT1, SYNDIG1, WASH1, PML, EEA1, CFTR, WLS, etc	hsa-miR-508hsa-miR-509-2hsa-miR-526b
GO:0005884	Actin filament	14	5.93×10−3	MYO5A, DMTN, MYO1B, RHOQ, ACKR2, ACTN2, MYO9B, GNG12, etc	hsa-miR-508hsa-miR-509-2hsa-miR-526b
GO:0015629	Actin cytoskeleton	33	5.99×10−3	MTSS1, DMTN, CNN3, SHROOM4, ONECUT2, KNTC1, PDLIM2, ARHGAP35, etc	hsa-miR-508hsa-miR-509-2hsa-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-508hsa-miR-509-2hsa-miR-526b
GO:0005254	Chloride channel activity	13	3.95×10−3	GABRG1, SLC1A4, GABRG2, GABRG3, GABRA4, CLIC4, SLC26A7, CLIC5, etc	hsa-miR-508hsa-miR-509-2hsa-miR-526b
GO:0005112	Notch binding	7	4.85×10−3	NOV, NOTCH1, HIF1AN, DNER, AAK1, JAG1, GALNT11	hsa-miR-508hsa-miR-509-2hsa-miR-526b
GO:0003779	Actin binding	41	4.86×10−3	MYO5A, LRRC10, DMTN, SHROOM2, CNN3, WASH1, WASF3, IMPACT, etc	hsa-miR-508hsa-miR-509-2hsa-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-508hsa-miR-509-2hsa-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-508hsa-miR-509-2hsa-miR-526b

Abbreviations: BP, biological process; CC, cellular component; GO, gene ontology; MF, molecular function.

Figure 6

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-508hsa-miR-509-2hsa-miR-526b
hsa04010	MAPK signaling pathway	37	7.52×10−3	FGF7, PPP3R1, PPM1A, PPP3R2, NFKB1, GNG12, TNFRSF1A, KRAS, etc	hsa-miR-508hsa-miR-509-2hsa-miR-526b
hsa04024	cAMP signaling pathway	29	1.68×10−2	ADCY1, ADCY7, ATP1B4, OXTR, NFKB1, GRIN3B, GABBR2, CNGB1, etc	hsa-miR-508hsa-miR-509-2hsa-miR-526b
hsa04068	FoxO signaling pathway	21	2.35×10−2	EGFR, IRS2, SGK1, RBL2, NLK, TGFBR1, SMAD4, SMAD3 etc.	hsa-miR-508hsa-miR-509-2hsa-miR-526b
hsa04151	PI3K–Akt signaling pathway	44	2.94×10−2	CSH2, FGF7, PHLPP2, TLR2, NFKB1, GNG12, IL7R, PTEN, etc	hsa-miR-508hsa-miR-509-2hsa-miR-526b
hsa04150	mTOR signaling pathway	11	4.07×10−2	PRKCA, MAPK1, RPS6KA3, TSC1, IGF1, CAB39, RICTOR, PRKAA2, etc	hsa-miR-508hsa-miR-509-2hsa-miR-526b

Abbreviation: KEGG, Kyoto Encyclopedia of Genes and Genomes.

Figure 7

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

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

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 (Figure 10).

Figure 10

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.2 LNM in the early stage affects the prognosis of CESC.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.31 In gastric cancer, hsa-miR-508 inhibits the invasion of gastric cancer cells by targeting S-phase kinase-associated protein 2.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.33 Moreover, miRNA microarray profiling showed that hsa-miR-508 expression is downregulated in lymph node metastatic tissues.34 In addition, Zhai et al indicated that hsa-miR-508 inhibits renal cell cancer cell invasion and migration.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.36–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.41 Additionally, low expression levels of hsa-miR-509 were associated with poor outcomes in pancreatic cancer patients.39 hsa-miR-526b can inhibit tumor progression by targeting MMP1.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.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α.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.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.13,46,47 Zhang et al found that the MAPK pathway was associated with LNM in head and neck cancer.48 Research has shown that the Akt/ mTOR pathway is activated in medullary thyroid carcinoma with lymph node metastases.49 In the coexpression regulatory network of target genes, multiple studies have shown that a TP53 gene mutation increases the risk of CESC.50–52 MMP1 has been associated with LNM in esophageal squamous cell cancer,53 lung cancer,54 and gastric cancer.55 Liu et al reported that MMP1 promotes esophageal squamous cell carcinoma cell migration and invasion via the PI3K/ Akt pathway.53 In CESC, NOTCH1 has been associated with LNM, and high expression of NOTCH1 indicates a poor prognosis in patients.56 At the same time, Park et al found that NOTCH1 could be a molecular marker of LNM in papillary thyroid cancer.57 In pancreatic and colorectal cancer, SMAD4 was found to be associated with LNM and survival.58,59 Another gene in the network, NFKB1, was reported to increase the risk of LNM in gastric cancer and oral cancer.60,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.

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

Abbreviations: ID, miRNA ID; FC, fold change; FDR, false discovery rate.