Screening for characteristic microRNAs between pre-invasive and invasive stages of cervical cancer

The aim of the present study was to investigate the characteristic microRNAs (miRNAs) expressed during the pre-invasive and invasive stages of cervical cancer. A gene expression profile (GSE7803) containing 21 invasive squamous cell cervical carcinoma samples, 10 normal squamous cervical epithelium samples and seven high-grade squamous intraepithelial cervical lesion samples, was obtained from the Gene Expression Omnibus. Differentially expressed genes (DEGs) were identified using significance analysis of microarray software, and a Gene Ontology (GO) enrichment analysis was conducted using the Database for Annotation, Visualization and Integrated Discovery. The miRNAs that interacted with the identified DEGs were selected, based on the TarBase v5.0 database. Regulatory networks were constructed from these selected miRNAs along with their corresponding target genes among the DEGs. The regulatory networks were visualized using Cytoscape. A total of 1,160 and 756 DEGs were identified in the pre-invasive and invasive stages of cervical cancer, respectively. The results of the GO enrichment demonstrated that the DEGs were predominantly involved in the immune response and the cell cycle, in the pre-invasive and invasive stages, respectively. Furthermore, a total of 18 and 26 characteristic miRNAs were screened in the pre-invasive and invasive stages, respectively. These miRNAs may be potential biomarkers and targets for the diagnosis and treatment of the different stages of cervical cancer.


Introduction
Cervical cancer is the most common gynecological malignancy, and is the second leading cause of cancer-associated mortality in females worldwide (1,2). One reason for the high levels of prevalence of this cancer is the lack of awareness and early detection approaches (3,4). Therefore, understanding the underlying molecular mechanisms of cervical cancer, and establishing more effective therapies are important areas of ongoing research.
The identification and characterization of key microRNAs (miRNAs) that participate in cervical cancer, is essential for determining the underlying mechanisms of this disease and establishing novel therapeutic strategies. miRNAs are 20-24 nt RNAs that are derived from distinct hairpin precursors in animals, plants and fungi, which bind to complementary sequences on target mRNAs (5,6). miRNAs regulate gene expression by cleaving target mRNAs, and by translational suppression at the post-transcriptional level (7). Previous studies have shown that miRNAs have important roles in various biological and metabolic processes, including cell growth, apoptosis, viral infection, differentiation, signal transduction and cancer (8)(9)(10)(11). Numerous studies have demonstrated that miRNAs are involved in the initiation and progression of cancer, and may be potential biomarkers for the diagnosis and prognosis of tumors, in addition to functioning as potential therapeutic targets (12)(13)(14). Therefore, it may be beneficial to identify novel miRNAs to act as diagnostic and therapeutic biomarkers, or therapeutic targets, in cervical cancer.
Recently, molecular network analysis technology, combined with gene expression profile data, has exhibited potential in a number of areas, including classification of diseases and the identification of novel therapeutic targets (15,16). In the present study a microarray dataset of healthy and malignant cervical samples was downloaded from the Gene Expression Omnibus (GEO) database. Differentially expressed genes (DEGs) were identified between these groups. Based on the TarBase v5.0 database, regulatory networks were constructed from selected miRNAs and their corresponding target genes from the identified DEGs. Key miRNAs, which may be used as potential biomarkers or therapeutic targets in cervical cancer, were subsequently identified.

Materials and methods
Affymetrix microarray data. A gene expression profile generated by Zhai et al (17) was used in the present study, which was deposited in the GEO database (http://www.ncbi.nlm.nih. gov/geo/query/acc.cgi?acc= GSE7803). This gene expression profile is based on the GPL96 platform (Affymetrix Human Genome U133A Array). A total of 38 samples were available, including 21 invasive squamous cell cervical carcinoma (SCC) samples, ten normal squamous cervical epithelium (NE) samples and seven high-grade squamous intraepithelial cervical lesion (HSIL) samples.
Screening of DEGs. In order to identify the DEGs, the original GSE7803 dataset was converted into an identifiable expression form and was normalized. Probe sets were mapped to the National Centers of Biotechnology Information genes (http://www.ncbi.nlm.nih.gov). Probe sets that corresponded to numerous genes or to no genes were removed from subsequent analyses. For genes that corresponded with numerous probe sets and had a plurality of expression values, the expression values were averaged. Subsequently, the SAMR package (18) in R and a significance analysis of microarray (SAM) were used to identify the DEGs between the samples (19). SAM software is a practical tool used for detecting significantly expressed genes, and for controlling the proportion of falsely detected genes. In the present study, genes with a fold-change >1.2 and a false discovery rate (FDR) <0.05 were selected as DEGs. In addition, the identified DEGs were divided into two groups: DEGs from the NE and HSIL samples were considered pre-invasive DEGs, whereas DEGs from the HSIL and invasive SCC samples were considered invasive DEGs.
Functional enrichment analysis of DEGs. The Database for Annotation, Visualization and Integrated Discovery (DAVID; http://david.abcc.ncifcrf.gov/) is a web-accessible program that provides a comprehensive set of functional annotation tools, which may be used by investigators to understand the underlying biological functions of large lists of genes (20). The present study used DAVID to perform a Gene Ontology (GO) enrichment analysis of the identified DEGs. Based on hypergeometric distribution, GO terms were enriched, and numerous testing corrections were conducted using the Benjamini-Hochberg method (21). An FDR<0.05 was set as the cut-off value.
Construction of regulatory networks. TarBase is a database that contains a manually curated collection of experimentally supported miRNA targets from a animal, pant and viral species of central scientific interest (22). TarBase v5.0 is the updated and extended version of the TarBase database, with >1,300 experimentally supported miRNA-target interactions (MTIs). It contains 1,094 human MTIs between 285 miRNAs and 1,721 target genes.
In the present study, human miRNA target gene data were downloaded from the TarBase v5.0 database (http://diana. cslab.ece.ntua.gr/tarbase/). miRNAs that interacted with the identified DEGs were then selected. Subsequently, MTIs regulatory networks were constructed from these selected miRNAs and their corresponding target genes within the DEGs. The MTIs regulatory networks were visualized by Cytoscape (23). In addition, the MTIs regulatory networks were divided into two groups: The regulatory network constructed from the selected miRNAs and the pre-invasive DEGs was termed the pre-invasive regulatory network, whereas the regulatory network constructed from the selected miRNAs and the invasive DEGs was termed the invasive regulatory network.
Comparison of the regulatory networks. In order to determine the differences between the pre-invasive and invasive stages of cervical cancer, regulatory networks were constructed and compared. Regulatory networks may be characterized by topological properties, such as degree (24). Degree is defined as the number of edges per node, which indicates the number of interacting partners. The present study used Freeman's degree centrality to analyze the degree of the regulatory networks (25). Freeman's degree centrality consists of ingoing (in-degree) and outgoing degree (out-degree). In-degree refers to the number of links a node receives from other nodes, whereas out-degree refers to the number of links originating from a particular node.

Results
DEG analysis. The original GSE7803 dataset was downloaded from the GEO database, and the DEGs were identified using SAM. Genes with a fold-change >1.2 and an FDR <0.05 were classed as DEGs. A total of 1,160 pre-invasive and 756 invasive DEGs were identified. In addition, 2,001 DEGs were identified from the NE and invasive SCC samples.
GO analysis of DEGs. In order to study the DEGs that contributed to cervical cancer, a GO enrichment analysis for the pre-invasive and invasive DEGs was performed using DAVID software. The pre-invasive DEGs (e.g. PSMB10, POU2AF1, ST6GAL1, CLU, SERPING1 and APOL2) were predominantly involved in the immune response, such as the acute inflammatory response (FDR=3.06E-04; Table I). By contrast, the invasive DEGs (e.g. TTK, AURKA, BRCA2, PSMC3IP, CDK10 and TUBG1) were predominantly involved in the regulation of the cell cycle, such as Cell Cycle (FDR=1.25E-19; Table II).
Construction of regulatory networks. Based on human MTIs data, pre-invasive and invasive regulatory networks were constructed. The pre-invasive regulatory network consisted of 80 pairs of regulatory interactions between 18 miRNAs and 66 pre-invasive DEGs (Fig. 1). The invasive regulatory network consisted of 64 pairs of regulatory interactions between 26 miRNAs and 51 invasive DEGs (Fig. 2). The highest out-degree was observed in miR-124, in the pre-invasive as well as the invasive regulatory networks.
Comparisons between the regulatory networks. Based on the topological properties of the networks, the similarities and differences between the pre-invasive and invasive regulatory networks were identified. The invasive regulatory network (Fig. 2) consisted of many smaller sub-networks and the out-degree of miRNAs was decreased, compared with those in the pre-invasive regulatory network (Fig. 1). For example, there were 14 DEGs associated with miR-1, and 21 DEGs associated with miR-124 in the pre-invasive regulatory network (Fig. 1). However, only eight and nine DEGs were associated with miR-1 and miR-124 in the invasive regulatory network, respectively (Fig. 2).  A total of 10 common miRNAs were identified in the regulatory networks. Three miRNAs: miR-1, miR-124 and miR-16, had a degree change >5. In addition, there were eight miRNAs that were only detected in the pre-invasive regulatory network (Fig. 1), including miR-126 and miR-199a. By contrast, there were 16 miRNAs that were only detected in the invasive regulatory network (Fig. 2), including miR-127, miR-143, miR-17-5p, miR-26a, miR-29a, miR-34a and miR-375.

Discussion
Malignant transformation during tumor progression results from a series of genetic alterations (26). In order to gain a better understanding of the genetic changes that occur during the progression of cervical cancer, a gene expression profile (GSE7803) was analyzed using a bioinformatics approach. In the present study, a total of 756 invasive DEGs, 1,160 pre-invasive DEGs, and 2,001 DEGs from invasive SCC and NE samples, were identified. These findings are in accordance with those of previous studies, which have consistently shown that the expression of genes is markedly altered in invasive tumor cells, compared with that of noninvasive and normal cells (27,28). Furthermore, the results of a GO enrichment of the identified DEGs, indicated that the expression of key genes differs between the pre-invasive and invasive stages of cervical cancer.
Clusterin (CLU) was initially identified as a secreted glycoprotein that has a cytoprotective role. However, numerous intracellular CLU variants have recently been identified in diverse pathological conditions (29)(30)(31). Furthermore, recent studies have shown that CLU is involved in various biological functions, such as cell death, tumor progression and neurodegenerative disorders (32,33). A previous study used DNA microarray data to identify novel candidate molecular markers for cervical cancer diagnosis and therapy, and observed the downregulation of human C1 inhibitor (SERPING1) in invasive cervical carcinoma cells (34). In addition, a recent genomic study demonstrated that apolipoprotein L2 (APOL2) is markedly upregulated in cervical cancer (35). These findings, as well as the results of the present study, indicate that CLU, SERPING1 and APOL2 may have important roles in the progression of cervical cancer. Figure 1. Pre-invasive regulatory network. Regulatory network constructed by microRNAs and differentially expressed genes from normal squamous cervical epitheilum and high grade squamous intraepithelial cervical lesion samples. Red nodes represent target differentially expressed genes and green nodes represent microRNAs. Blue lines represent microRNA-target regulatory interactions (in-degree), and arrows indicate microRNA target differentially expressed genes (out-degree). The size of each green node represents the out-degree. As the out-degree increases, the associated green node becomes larger.
TTK has been shown to be associated with metastasis via chromosomal instability, in a previous study, which aimed to identify genes associated with the progression and metastasis of advanced cervical cancer following radiotherapy (36). Furthermore, genetic variants of Aurora A kinase (AURKA) have been shown to be associated with a radiotherapy-induced early adverse reaction in patients with cervical cancer (37). Previous studies have demonstrated that both BRCA1 and BRCA2 participate in a common DNA damage response pathway, and are involved in the activation of homologous recombination and double-strand break repair (38). By contrast, Narayan et al (39) reported the downregulation of BRCA1 in a small subset of patients with cervical cancer. These previous findings and the results of the present analysis suggest that TTK, AURKA and BRCA2 may participate in the progression of cervical cancer.
In order to obtain the upstream regulatory information of the DEGs, two regulatory networks were constructed based on the TarBase v5.0 database. These regulatory networks were then compared, and the common and specific miRNAs were identified. The miRNA with the highest out-degree was shown to be miR-124, in the pre-invasive as well as the invasive regulatory networks. miR-124 has previously been shown to be the most abundant miRNA expressed in neuronal cells (40). Furthermore, previous studies have shown that the upregulation of miR-124 induces neuronal differentiation of various tumor cell lines in mice (41)(42)(43). Wilting et al (44) previously demonstrated that the silencing of miR-124 expression, by methylation, inhibited the development of cervical carcinoma. These results suggest that miR-124 may be a potential therapeutic target for cervical cancer therapy.
Of the eight miRNAs specific to the pre-invasive regulatory network, miR-126 has previously been reported to be downregulated in cervical cancer tissues (45), and miR-199a has previously been suggested as a potential therapeutic target for cervical cancer therapy (46). miR-126 is a human miRNA that is expressed only in endothelial cells, throughout capillaries as well as in larger blood vessels (47), and acts upon Figure 2. Invasive regulatory network. Regulatory network constructed by microRNAs and differentially expressed genes from invasive squamous cell cervical carcinoma samples and high grade squamous intraepithelial cervical lesion samples. Red nodes represent target differentially expressed genes and green nodes represent microRNAs. Blue lines represent microRNA-target regulatory interactions (in-degree), and arrows indicate microRNA target differentially expressed genes (out-degree). The size of each green node represents the out-degree. As the out-degree increases, the green node becomes bigger. various transcripts in order to control angiogenesis (48). miR-126 has been identified as a tumor suppressor and as an oncogene, depending on the type of cancer involved. Inhibition of cancer progression by miR-126 is achieved through the negative control of proliferation, migration, invasion and cell survival. However, miR-126 may also support cancer progression through the promotion of blood vessel formation and inflammation at the site of activation (49). According to these previous findings and the results of the present study, miR-126 may be a potential biomarker for the diagnosis of cervical cancer, and a therapeutic target for the pre-invasive stage of this disease.
Of the 16 miRNAs specific to the invasive regulatory network, seven have been reported in previous studies and described as being upregulated or downregulated in cervical cancer. These include miR-127, miR-143, miR-17-5p, miR-26a, miR-29a, miR-34a and miR-375 (45,(50)(51)(52)(53)(54). Lee et al (46) demonstrated that the expression of miR-127 was significantly increased in patients with invasive squamous cell carcinoma, which had metastasized to the lymph nodes. The results of the present study are in accordance with those of previous studies, which indicate that miR-127 may be a marker for lymph node metastasis in invasive cervical cancer. miR-143 is highly conserved in vertebrates (55) and changes in miR-143 expression have frequently been implicated in cancer (56)(57)(58). Furthermore, the upregulation of miR-143 has previously been observed in a hepatocellular carcinoma model during tumor metastasis, through repression of FNDC38 (59). However, reduced expression of miR-143 has also been observed in a range of cancer stages, including at very early stages (60). The results of previous studies and of the present study indicate that miR-143 may be involved in tumor progression, and may be a candidate for RNA-targeted treatment of tumors (61). Wang et al (52) previously reported that miR-375 is downregulated in squamous cervical cancer, and inhibits cell migration and invasion by targeting the transcription factor, SP1. This finding indicates that deregulation of miR-375 may have an important role in the malignant transformation of cervical cancer cells. However, the elucidation of the underlying molecular mechanisms of miR-17-5p, miR-26a, miR-29a and miR-34a in the progression of cervical cancer, and the use of other miRNAs screened in the present study as biomarkers or therapeutic targets in cervical cancer require further investigation.
In conclusion, a total of 1,160 and 756 DEGs were identified in the pre-invasive and invasive stages of cervical cancer, respectively. The GO enrichment analysis demonstrated that the DEGs were primarily involved in the immune response and regulation of the cell cycle, in the pre-invasive and invasive stages, respectively. These findings indicate that the expression of key genes differs between the pre-invasive and invasive stages of cervical cancer progression. Based on the analysis of the regulatory networks, a total of 18 and 26 key miRNAs were screened in the pre-invasive and invasive stages, respectively. It is hypothesized that these miRNAs are involved in the malignant transformation of cervical cancer cells. In addition, these miRNAs may have a function as novel biomarkers in cervical cancer diagnosis and detection, and as therapeutic targets in this disease. Further studies in independent patient cohorts are required, in order to validate the potential roles of these miRNAs.