Prediction of a miRNA‐mRNA functional synergistic network for cervical squamous cell carcinoma

Cervical squamous cell carcinoma (CSCC) accounts for a significant proportion of cervical cancer; thus, there is a need for novel and noninvasive diagnostic biomarkers for this malignancy. In this study, we performed integrated analysis of a dataset from the Gene Expression Omnibus database to identify differentially expressed genes (DEGs) and differentially expressed miRNAs (DEmiRNAs) between CSCC, cervical intraepithelial neoplasia (CIN) and healthy control subjects. We further established protein‐protein interaction and DEmiRNA‐target gene interaction networks, and performed functional annotation of the target genes of DEmiRNAs. In total, we identified 1375 DEGs and 19 DEmiRNAs in CIN versus normal control, and 2235 DEGs and 33 DEmiRNAs in CSCC versus CIN by integrated analysis. Our protein‐protein interaction network indicates that the common DEGs, Cyclin B/cyclin‐dependent kinase 1 (CDK1), CCND1, ESR1 and Aurora kinase A (AURKA), are the top four hub genes. P53 and prostate cancer were identified as significantly enriched signaling pathways of common DEGs and DEmiRNA targets, respectively. We validated that expression levels of three DEGs (TYMS, SASH1 and CDK1) and one DEmiRNA of hsa‐miR‐99a were altered in blood samples of patients with CSCC. In conclusion, a total of four DEGs (TYMS, SASH1, CDK1 and AURKA) and two DEmiRNAs (hsa‐miR‐21 and hsa‐miR‐99a) may be involved in the pathogenesis of CIN and the progression of CIN into CSCC. Of these, TYMS is predicted to be regulated by hsa‐miR‐99a and SASH1 to be regulated by hsa‐miR‐21.

exploring the mechanism in progression of CIN developing into CSCC.
miRNAs are a group of small, noncoding RNAs of [20][21][22] nucleotides that modulate about 60% of protein-coding genes [5]. miRNAs may play a key role in novel diagnostic, prognostic and therapeutic markers in clinical oncology. In addition, no curative therapy can be used for CSCC, and the current methods have only limited efficacy. Therefore, it is essential to identify accurate and credible biomarkers for CSCC diagnosis and treatment.
In our study, we intended to obtain more credible results than possible with individual study via integrated analysis. We performed functional annotation of differentially expressed genes (DEGs) and a CSCCspecific miRNA-target gene network to seek key DEGs and differentially expressed miRNAs (DEmiR-NAs) in CIN and biomarkers in the development of CIN into CSCC.

DEGs identification of CSCC and CIN
We searched gene expression datasets of CSCC and CIN from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo) [6]. Search keywords were ['cervical intraepithelial neoplasia' (MeSH Terms) OR cervical intraepithelial neoplasia (All Fields)] OR {['cervical' (MeSH Terms) OR cervical (All Fields)] AND ['carcinoma' (MeSH Terms) OR carcinoma (All Fields)]} AND 'gse' (Filter) [7]. The study types were limited to 'expression profiling by array'. Datasets that meet the following criteria would be included in our study: (a) the selected dataset must be genome-wide mRNA or miRNA transcriptome data; (b) these data were obtained from tissues of the CIN, CSCC and normal control (N) (no drug stimulation or transfection); and (c) normalized or raw datasets were considered in this study. Finally, three datasets of mRNA and miRNA were screened and included. The study methodologies conformed to the standards set by the Declaration of Helsinki.

Identification of DEGs and DEmiRNAs in CIN versus N and CIN versus CSCC
The method used for P-value consolidation is the inverse normal method in the LIMMA package and metaMA [8]. The adopted standard is false discovery rate (FDR) < 0.05, and all datasets have the same direction of different expression. Finally, the DEGs and DEmiRNAs of CIN versus N and CIN versus CSCC were obtained. Intersection of DEGs and DEmiRNAs was obtained from CIN versus N and CSCC versus CIN.

Gene ontology and pathway enrichment analysis of DEGs
To identify the characteristic biological function and potential pathways of common DEGs, we performed gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment by using the online software GENECODIS3 [9]. All common DEGs were analyzed by GO and KEGG enrichment analysis using the R language (GSEABase package).

Protein-protein interaction network construction of common DEGs
Cytoscape 3.5.0 was used to search for all common DEGs based on the existing data of protein interaction in the String database. All common DEGs in CIN versus N and CIN versus CSCC groups were used to construct the protein-protein interaction (PPI) network [10][11][12].
DEmiRNA-target interaction network in CIN versus N and CIN versus CSCC miRTarBase (http://mirtarbase.mbc.nctu.edu.tw/php/index. php) [13] provides the latest and extensive experimentally validated miRNA-target interaction information. The potential target genes of miRNAs were predicted by miR-TarBase, which is an experimentally validated miRNA-target interactions database in CIN versus N and CIN versus CSCC, respectively.

Functional annotation of target mRNAs
Based on the earlier analysis results, we obtained seven target genes of DEmiRNAs. GO functional enrichment and KEGG functional enrichment analysis were performed on target genes using GENECODIS3 (http://genecodis.cnb.csic.es/a nalysis). FDR < 0.05 was defined as the criterion of statistical significance.

Quantitative RT-PCR confirmation
According to the results of GEO integrated analysis, we selected two DEmiRNAs (hsa-miR-99a and hsa-miR-21) and three targets including TYMS, SASH1 and Aurora kinase A (AURKA) in CIN versus CSCC as candidate genes. A total of 12 blood samples were collected from four N subjects, four patients diagnosed with CIN and four patients diagnosed with CSCC. Informed written consent was obtained from all participants, and research protocols were approved by the Ethics Committee of our hospital.

Differential expression analysis of genes in CIN and CSCC
After retrieving, we obtained three microarray studies of mRNA and three microarray studies of miRNA according to the inclusion criteria from the GEO database. The characteristics of the individual database for the integrated analysis are displayed in Table 1.
A total of 1375 DEGs were obtained with FDR < 0.05 in CIN compared with N, among which the expression level of 719 genes was increased and the expression level of 656 genes was decreased. Likewise, 2235 DEGs with 1213 up-regulated and 1022 down-regulated genes were obtained in CSCC compared with CIN. The top 20 DEGs in CIN versus N and CESC versus CIN are listed in Tables 2  and 3, respectively. The hierarchical clustering heatmap of the top 100 most significantly up-regulated or down-regulated genes in CIN versus N and CSCC versus CIN is shown in Fig. 1A,B. In total, 392 common DEGs were obtained by taking the intersection of CIN versus N and CSCC versus CIN (Fig. 1C).
Nineteen DEmiRNAs (7 up-regulated DEmiRNAs and 12 down-regulated DEmiRNAs) were obtained with P < 0.05 in CIN compared with N. Likewise, 33 DEmiRNAs with 18 up-regulated and 15 down-regulated genes were obtained in CSCC compared with CIN. The top 20 DEmiRNAs in CIN versus N and CESC versus CIN are listed in Tables 4 and 5, respectively. The hierarchical clustering heatmap of DEmiRNAs is shown in Fig. 1D,E. A total of six common DEmiRNAs were obtained by taking the intersection of CIN versus N and CSCC versus CIN (Fig. 1F).

PPI network and module analysis of common DEGs
To identify potential interactions between common DEGs, we constructed a PPI network. In total, 259 nodes (genes) and 886 edges were identified among DEGs in the results, which are shown in Fig. 3. Among them, the higher degree genes are CDK1 (de- . The expression levels of these genes in three databases are shown in Fig. 4.

DEmiRNA-target interaction network in CIN versus N and CIN versus CSCC
In total, 393 DEmiRNA-target interaction pairs were obtained. The DEmiRNA-target regulatory network was constructed based on these DEmiRNA-target interaction pairs, which consisted of 326 nodes and 393 edges (Fig. 5A). Based on the CIN-specific DEmiRNA-target interaction network, has-miR-24 (degree = 64), has-miR-149 (degree = 45) and has-miR-519d (degree = 38) were the top three DEmiR- NAs that regulated most DEGs. All of these DEmiR-NAs were down-regulated in CIN based on the GEO database. Likewise, 1396 DEmiRNA-target interaction pairs were obtained, which consisted of 828 nodes and 1396 edges (Fig. 5B). Based on the CIN-specific DEmiRNA-target interaction network, has-miR-26b (degree = 212), has-miR-16 (degree = 98) and has-let-7a (degree = 82) were the top three DEmiRNAs. All of these DEmiR-NAs were down-regulated in CIN based on the GEO database. Also, the intersection of seven DEmiRNAtarget pairs was obtained in the above two sections. Among them, two interaction pairs (hsa-miR-21-SASH1 and hsa-miR-99a-TYMS) were involved in CIN and the progression of CIN into CSCC (Fig. 6). Topological properties of the DEmiRNA-mRNA interaction pair of CIN versus N and CESC versus CIN are shown in Tables S2 and S3, respectively.  Table S4, among which pathway enrichment analysis for TYMS, CDK1, AURKA and SASH1 was further performed in SMPDB, BIOCYC, KEGG and REACTOME databases, which is shown in Table 6.

Quantitative RT-PCR confirmation
To indicate the results of integrated analysis, we selected two DEmiRNAs (hsa-miR-99a and hsa-miR-21) and three target genes including TYMS, SASH1 and AURKA in CIN versus CSCC. Based on the results of quantitative RT-PCR, the expression of hsa-miR-99a, hsa-miR-21, AURKA and SASH1 was down-regulated, whereas the expression of TYMS was up-regulated in CIN compared with CSCC. The expression of TYMS, SASH1 and hsa-miR-99a was consistent with the results of our integrated analysis (Fig. 8).

Discussion
CIN is a class of precancerous lesions of the CSCC, but the mechanisms of CIN developing into CSCC need to be elucidated more clearly [2]. In this study,     and TYMS) under the regulation of two DEmiRNAs (hsa-miR-21 and hsa-miR-99a) were associated with CIN and CSCC. CDK1 and AURKA are two major hallmarks of both CIN versus N and CIN versus CSCC. CDK1 has more than 70 regulatory targets, which play a vital role in the control of the cell cycle. In transcription and cell-cycle progression, various target substrates are directly phosphorylated by CDK1 in response to various stimuli [14]. As previously described, many patients with cancers had aberrant activation of CDKs and their modulators. Abnormal cell proliferation and genomic instability were caused by dysregulation of CDKs [15]. According to our integrated analysis, up-regulated CDK1 was modulated by hsa-miR-205 and hsa-miR-24 in the tissues of CIN versus N and by hsa-miR-195 and hsa-miR-497 in CIN versus CSCC, which have the same pattern in a previous study [16]. In the PPI analysis, CDK1 has the highest degree among the hub proteins. Also, CDK1 was enriched in the p53 signaling pathway in the KEGG analysis. Luo et al. [16] reported that CDK1 played a complicated role in regulating genetic networks involved in the progression of cervical cancer. Prognosis of advanced stage cervical cancer may be enhanced by new therapeutics targeting CDK1 or its related pathways. This was consistent with our results, which indicated that CDK1 might serve as a biomarker for CIN and the progression of CIN into CSCC.
AURKA is a member of a family of mitotic serine/ threonine kinases. During mitosis and meiosis, AURKA is correlated with crucial processes, whose appropriate function is integral for normal cell proliferation [17]. The first study related to the family of kinases in tumorigenesis reported that Aurora A and B were overexpressed in primary breast and colon tumor samples. Emerging studies found that AURKA was amplified or overexpressed in other tumor types, such as pancreatic, ovarian and hepatocellular tumors [18][19][20]. Nae-Fang Twu et al. [21] have shown that expression of AURKA was significantly increased in carcinoma and CIN 3 compared with the normal cervix. In our study, AURKA was up-regulated in CIN versus N and CIN versus CSCC, which showed the same pattern with the previous study [21]. In the PPI analysis, AURKA was among the top four hub proteins. In the GO and KEGG analysis, AURKA was enriched in the items of mitotic cell cycle, cell division and oocyte meiosis.
Recently, many studies reported that SASH1 played a crucial role in inhibiting various tumors. For example, He et al. [22] reported that the gene of SASH1 inhibited the metastatic progression of hepatocarcinoma cells via regulating the sonic hedgehog signaling pathway. Ren et al. [23] showed that up-regulation of SASH1 inhibited proliferation and migration of ovarian carcinoma cells. In line with these studies, another study [24] showed that SASH1 was down-regulated in cervical cancer, indicating that SASH1 might play a  negative role in cervical cancer. Up-regulated TYMS by hsa-miR-99a and down-regulated SASH1 by hsa-miR-21 were detected in CIN versus N and CIN versus CSCC in our study, indicating its important role in CIN, and up-regulated TYMS and down-regulated SASH1 might serve as a biomarker for CIN and the progression of CIN into CSCC. In addition, the expression of these two DEmRNAs in CIN in our integrated analysis was consistent with our quantitative RT-PCR results.
Recently, miRNAs have been described as potential diagnostic or prognostic markers for many cancers and can function as neoteric targets for cancer therapies, including cervical cancer [25,26]. As previously reported, miR-21-5p was expressed abnormally in patients with CSCC [27]. Likewise, other researchers have reported that miR-21 influenced tumorigenesis in the cervical squamous cell and served as an oncology miRNA (oncomiRNA) in cervical cancer [28][29][30]. Upregulated hsa-miR-21 was detected in our integrated analysis. There were 50 and 34 targets of hsa-miR-21 in the DEmiRNA-target interaction network in CIN versus N and CIN versus CSCC, respectively.
Many studies have indicated that miR-99a is correlated with tumor pathogen [31][32][33]. The miR-99a is down-regulated in human cancers, such as endometrioid endometrial carcinoma, suggesting that miR-99a may inhibit tumor progression [34]. However, the role of miR-99a in cervical cancer still needs to be elucidated. A previous study provided evidence that miR-99a was down-regulated in cervical cancer tissues [35]. Our results revealed that hsa-miR-99a was down-regulated in CIN versus N and CIN versus CSCC. There were 7 and 12 targets of hsa-miR-99a in the DEmiRNA-target interaction network in CIN versus N and CIN versus CSCC, respectively.
In conclusion, four DEGs (TYMS, SASH1, CDK1 and AURKA) and two DEmiRNAs (hsa-miR-21 and hsa-miR-99a) may be involved in the pathogenesis of CIN and the progression of CIN into CSCC, which might contribute to developing novel diagnostic and therapeutic strategies for early-stage CIN. Among them, TYMS was regulated by hsa-miR-99a and SASH1 was regulated by hsa-miR-21. Cell cycle (Homo sapiens) AURKA Gene expression (transcription) (Homo sapiens), metabolism of proteins (Homo sapiens) SASH1 - Fig. 8. Quantitative RT-PCR results of three DEGs (TYMS, SASH1 and AURKA) and two DEmiRNAs (hsa-miR-21 and hsa-miR-99a) in CIN versus N. The error bars indicate SD. The statistical test used to determine significance was one-way ANOVA. All of the data had biological duplicates, and each gene had three duplicates. *SASH1 has significant difference. performed the experiment. JJ and YD performed the data analyses. YD and XY contributed significantly in writing the manuscript. All authors read and approved the final manuscript.

Supporting information
Additional supporting information may be found online in the Supporting Information section at the end of the article. Table S1. GO (biological process, cellular component and molecular function) and KEGG items of common DEGs. Table S2. Topological properties of the DEmiRNA-mRNA interaction pair of CIN versus N. Table S3. Topological properties of the DEmiRNA-mRNA interaction pair of CESC versus CIN. Table S4. GO (biological process, cellular component and molecular function) and KEGG items of DEmiRNA targets.