NRP1 is targeted by miR-130a and miR-130b, and is associated with multidrug resistance in epithelial ovarian cancer based on integrated gene network analysis.

Multidrug resistance (MDR) in epithelial ovarian cancer (EOC) remains a public health issue for women worldwide, and its molecular mechanisms remain to be fully elucidated. The present study aimed to predict the potential genes involved in MDR, and examine the mechanisms underlying MDR in EOC using bioinformatics techniques. In the present study, four public microarray datasets, including GSE41499, GSE33482, GSE15372 and GSE28739, available in Gene Expression Omnibus were downloaded, and 11 microRNAs (miRNA; miRs), including miR‑130a, miR‑214, let‑7i, miR‑125b, miR‑376c, miR‑199a, miR‑93, miR‑141, miR‑130b, miR‑193b* and miR‑200c, from previously published reports in PubMed were used to perform a comprehensive bioinformatics analysis through gene expression analysis, signaling pathway analysis, literature co‑occurrence and miRNA‑mRNA interaction networks. The results demonstrated that the expression of neuropilin 1 (NRP1) was upregulated, thereby acting as the most important hub gene in the integrated gene network. NRP1 was targeted by miR‑130a and miR‑130b at the binding site of chromosome 10: 33466864‑3466870, which was involved in the axon guidance signaling pathway. These results suggested that alteration of the gene expression levels of NRP1 expression may contribute to MDR in EOC. These data provide important information for further experimental investigations of the drug resistance‑associated functions of NRP1 in EOC.

Introduction

Ovarian cancer is one of the leading contributors to mortality rates in women worldwide, with a five-year survival rate of 30–45% (1–3). There were ~21,880 new cases diagnosed and 13,850 cases of associated mortality reported in the United States in 2010 (4). Epithelial ovarian cancer (EOC) is the most common histological type, comprising 80–90% of all ovarian cancer cases (5). Chemotherapy is an important treatment for EOC, and patients with advanced EOC have high initial response rates to chemotherapy (≥80%); however, 70–80% of patients eventually relapse, with a progression-free survival rate of 18 months (6,7). Although chemotherapy treatment regimens have improved in previous decades, the efficacy of chemotherapy for EOC has improved simultaneously. Multidrug resistance (MDR) is the predominant reason for chemotherapy failure and poor prognosis in patients with EOC (8).

The role of microRNAs (miRNAs; miRs), and their target genes associated with MDR, in EOC have been investigated in previous years (9). As a cluster of small non-coding molecular RNAs, miRNAs are important in cell differentiation, proliferation, apoptosis, organism growth and the development of human disease (10–12). Mature miRNA primarily targets the 3′-untranslated region (3′-UTR) of an mRNA strand with a complementary sequence, and then induces mRNA degradation or inhibits mRNA translation at the post-transcriptional level (13–15). The dysfunction of miRNAs is associated with the development and metastasis of human cancer (16,17).

A single miRNA is able to target multiple mRNAs and, similarly, a single mRNA can be targeted by multiple miRNAs, forming multiple complex gene expression regulatory networks (18). However, the mechanisms underlying MDR in EOC remain to be fully elucidated. As integrated network analysis has been applied in investigations of EOC as a promising technique (19–23), it may provide important information in investigating the molecular mechanisms underlying MDR in EOC. In the present study four published microarray datasets (GSE41499, GSE33482, GSE15372 and GSE28739) and 11 miRNAs (miR-130a, miR-214, let-7i, miR-125b, miR-376c, miR-199a, miR-93, miR-141, miR-130b, miR-193b* and miR-200c) were downloaded from public databases in order to perform a comprehensive bioinformatics analysis through gene expression analysis, signaling pathway analysis, literature co-occurrence and miRNA-mRNA interaction networks. The aim of the present study was to identify any potential genes, and obtain their bioinformatics information, highly associated with MDR in EOC.

Materials and methods

Microarray datasets/miRNAs/miRNA target genes

The four microarray datasets, GSE41499, GSE33482, GSE15372 and GSE28739, were downloaded from the Gene Expression Omnibus database (GEO; http://www.ncbi.nlm.nih.gov/geo/) (28). All four microarray datasets satisfied the following four criteria: i) Data contained information regarding EOC genome-wide RNA expression; ii) data provided a comparison of EOC samples between chemotherapy resistance and chemotherapy sensitivity; iii) a minimum of three samples were included in each group; iv) each sample provided detailed information on chemotherapy resistance or sensitivity in EOC (Table I).

Table I

Characteristics of the datasets selected for the present study.

GEO accession	Reference	Chip	Experimental design	Probes	Chemotherapy response
					Sensitivity	Resistance
GSE41499	(24)	U133A	Paired cell lines PEO1 and PEO4	22,277	4	4
GSE33482	(25)	Agilent_HumanGenome	Paired cell lines A2780 and A2780cis	41,078	6	6
GSE15372	(26)	U133 Plus 2.0	Paired cell lines A2780 and Round5 A2780	54,675	5	5
GSE28739	(27)	Agilent_Human1Av2	Unpaired, tissues	16,096	20	30

GEO, Gene Expression Omnibus.

A total of 11 miRNAs (miR-130a, miR-214, let-7i, miR-125b, miR-376c, miR-199a, miR-93, miR-141, miR-130b, miR-193b* and miR-200c) were investigated in PubMed (http://www.ncbi.nlm.nih.gov/pubmed/; Table II), which were those reported to be highly associated with MDR in EOC (29–39). The above-mentioned 11 miRNAs were entered into TargetScan (http://genes.mit.edu/tscan/targetscan2003.html) and PicTar (http://pictar.mdc-berlin.de/) (40,41), respectively, to determine the relative miRNA target genes.

Table II

Characteristics of the miRNAs selected for the present study

miRNA	Target gene	Expression in EOC chemotherapy tissue/cell line	Reference
miR-130a	M-CSF	Downregulation	Sorrentino et al (29)
miR-214	PTEN	Upregulation	Yang et al (30)
let-7i	/	Downregulation	Yang et al (31)
miR-125b	Bak1	Upregulation	Kong et al (32)
miR-376c	ALK7	Upregulation	Ye et al (33)
miR-199a	CD44	Downregulation	Cheng et al (34)
miR-93	PTEN	Upregulation	Fu et al (35)
miR-141	KEAP1	Upregulation	van Jaarsveld et al (36)
miR-130b	CSF-1	Downregulation	Yang et al (37)
miR-193b*	/	Upregulation	Ziliak et al (38)
miR-200c	TUBB3	Upregulation	Prislei et al (39)

miR/miRNA, microRNA; EOC, epithelial ovarian cancer; M-CSF, macrophage colony-stimulating factor; PTEN, phosphatase and tensin homolog; Bak1, B cell lymphoma 2-antagonist/killer 1; KEAP1, Kelch-like ECH-associated protein 1; TUBB3, tubulin β3 class III;/, unavailable.

Gene expression analysis

The Benjamini-Hochberg (BH) (42) method was used to analyze gene expression, which was performed in GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r/) (43), a web tool allowing users to perform R statistical analysis without command line expertise. An adjusted P-value of P<0.05 was used as the screening criterion for statistically significantly expressed genes. The fold change (FC) method (44) was also used to estimate gene expression. When logFC<0, the expression of the genes was downregulated, whereas the expression of the genes was upregulated when logFC>0.

Pathway enrichment analysis

Genetic pathway enrichment analysis was performed in the Kyoto Encyclopedia of Genes and Genomes (KEGG) database using the Database for Annotation, Visualization and Integrated Discovery (DAVID) tool (http://david.abcc.ncifcrf.gov/) (45,46). Another web server, DIANA-miRPath (http://www.microrna.gr/miRPathv2) (47) was also used, as it is specifically designed for miRNA-targeted pathway analysis based on interaction levels. Fisher's exact probability method was used to determine the significant difference of pathway enrichment analysis, with P<0.05 indicating significance.

Co-occurrence analysis

Text mining methods from the literature and disease levels were combined to screen for MDR-associated genes in EOC, which were performed in the COREMINE (http://www.coremine.com/medical/#search) and IPAD (http://bioinfo.hsc.unt.edu/ipad/) (48) databases, respectively. The corresponding genes and the exact keywords 'drug resistance', 'drug resistance, multiple' and 'drug resistance, neoplasm' were input in COREMINE for co-occurrence analysis, and the disease information between differentially expressed genes and EOC were mined in the IPAD database. CytoScape2.6.1 software (49) was used to construct a graph of the association between genes and MDR.

Integrated gene network analysis

Integrated gene network analysis, based on miRNAs and their target genes, was performed using Pajek software (50). The topological characteristics of the integrated gene network comprised degree centralization (DC), betweenness centralization (BC), closeness centralization (CC) and clustering coefficient (CC′), which were calculated using the Pajek software. Degree of node indicates the number of adjacent nodes or connected edges each node has. The higher the number of neighbors (nodes and edges) a node has, the more importance it has in the network. Therefore, the node is also called a hub node (51). Correspondingly, the gene in the position of the hub node was termed the hub gene. The binding sites of the miRNA-target interactions were finally analyzed in StarBase (http://star-base.sysu.edu.cn/) (52), which was designed for deciphering miRNA-target interactions, including miRNA-mRNA interaction networks from large-scale CLIP-Seq data.

Results

Gene expression and miRNA target genes

Using the BH method in GEO2R, a total of 5,003 significantly expressed genes were obtained from GSE41499, 3,372 from GSE33482, 2,029 from GSE15372 and 267 from GSE28739. Among these, 2,505 genes were upregulated and 2,498 were downregulated in GSE41499, 1,487 genes were upregulated and 1,885 genes were downregulated in GSE33482, 798 genes were upregulated and 1,231 genes were downregulated in GSE15372, and 180 genes were upregulated and 87 genes were downregulated in GSE28739, respectively.

The present study also obtained 47,077 target genes using TargetScan and 1,675 target genes using PicTar, based on the previously mentioned 11 miRNAs (miR-130a, miR-214, let-7i, miR-125b, miR-376c, miR-199a, miR-93, miR-141, miR-130b, miR-193b* and miR-200c).

Genetic pathway enrichment analysis were performed in the KEGG database using the DAVID tool, based on upregulated genes and downregulated genes from the four microarray datasets (GSE41499, GSE33482, GSE15372 and GSE28739). A total of 11 upregulated signaling pathways were enriched in the KEGG database, including the mitogen-activated protein kinase (MAPK) signaling pathway, ubiquitin-mediated proteolysis, axon guidance, focal adhesion, neurotrophin signaling pathway, pathways in cancer, renal cell carcinoma, citrate cycle, terpenoid backbone biosynthesis, mismatch repair and Huntington's disease. In addition, seven downregulated signaling pathways were identified, including glycerolipid metabolism, pentose phosphate pathway, fructose and mannose metabolism, glutathione metabolism, proteasome, p53 signaling pathway and lysosome. The corresponding upregulated and downregulated genes are shown in Tables III and IV, respectively.

Table III

Common upregulated pathways and their corresponding genes.

Upregulated pathway	Upregulated gene
MAPK signaling pathway	FLNA, FLNC, MAPK8IP2, NF1
Ubiquitin mediated proteolysis	CDC16, CDC34, CUL4A, NEDD4, UBE2D1
Axon guidance	CXCR4, DPYSL2, FYN, NRP1, SLIT2
Focal adhesion	CAV1, FLNA, FLNC, FYN, PIK3CA
Neurotrophin signaling pathway	BDNF, PIK3CA, SHC3
Pathways in cancer	CSF2RA, EPAS1, FGF18, FGF5, FZD2, FZD7, PIK3CA, TGFB2
Renal cell carcinoma	EPAS1, PIK3CA, TGFB2
Citrate cycle	ACLY, ACO1, ACO2, FH, IDH3B, IDH3G, PDHA1, SDHA
Terpenoid backbone biosynthesis	ACAT2, HMGCR, IDI1, MVD
Mismatch repair	EXO1, RFC2, RFC4, RFC5
Huntington's disease	CLTB, NDUFB3, NDUFB8, NDUFS6, NDUFV1, SDHA

FLNA, filamin Aα; FLNC, filamin Cγ; MAPK8IP2, mitogen-activated protein kinase 8-interacting protein 2; NF1, neurofibromin 1; CDC16/34, cell division cycle 16/34; CUL4A, cullin 4A; NEDD4, neural precursor cell expressed developmentally downregulated 4 E3 ubiquitin protein ligase; UBE2D1, ubiquitin-conjugating enzyme E2D 1; CXCR4, chemokine (C-X-C motif) receptor 4; DPYSL2, dihydropyrimidinase-like 2; NRP1, neuropilin 1; SLIT2, slit homolog 2 (Drosophila); CAV1, caveolin 1; PIK3CA, phosphatidylinositol-4,5-bisphosphate 3-kinase (catalytic subunit α); BDNF, brain-derived neurotrophic factor; CSF2RA, colony-stimulating factor 2 receptor α; EPAS1, endothelial PAS domain protein 1; FGF18/5, fibroblast growth factor 18/5; FZD2/7, frizzled class receptor 2/7; TGFB2, transforming growth factor β2; ACLY, ATP citrate lyase; ACO1/2, aconitase 1/2; FH, fumarate hydratase; IDH3B/G, isocitrate dehydrogenase 3 (NAD+) β/γ; PDHA1, pyruvate dehydrogenase (lipoamide) α1; SDHA, succinate dehydrogenase complex, subunit A; ACAT2, acetyl-CoA acetyltransferase 2; HMGCR, 3-hydroxy-3-methylglutaryl-CoA reductase; IDI1, isopentenyl-diphosphate delta isomerase 1; MVD, mevalonate (diphospho) decarboxylase; EXO1, exonuclease 1; RFC2/4/5, replication factor C2/4/5; CLTB, clathrin light chain B; NDUFB3/8/6, NADH dehydrogenase (ubiquinone) 1β subcomplex 3/8/6; NDUFV1, NADH dehydrogenase (ubiquinone) flavoprotein 1.

Table IV

Common downregulated pathways and their corresponding genes.

Downregulated pathway	Downregulated gene
Glycerolipid metabolism	AGPAT3, PNPLA3, PPAP2C
Pentose phosphate pathway	PFKL, PGD, PGM1
Fructose and mannose metabolism	PFKL
Glutathione metabolism	ALDH9A1, HIBCH
Proteasome	ACACA, ALDH9A1, HIBCH
P53 signaling pathway	CDKN2C, MCM6, SKP2, TTK
Lysosome	FAS, RRM2

AGPAT3, 1-acylglycerol-3-phosphate O-acyltransferase 3; PNPLA3, patatin-like phospholipase domain containing 3; PPAP2C, phosphatidic acid phosphatase type 2C; PFKL, phosphofructokinase liver; PGD, phosphogluconate dehydrogenase; PGM1, phosphoglucomutase 1; ALDH9A1, aldehyde dehydrogenase 9 family member A1; HIBCH, 3-hydroxyisobutyryl-CoA hydrolase; ACACA, acetyl-CoA carboxylase α; CDKN2C, cyclin-dependent kinase inhibitor 2C; MCM6, minichromosome maintenance complex component 6; SKP2, S-phase kinase-associated protein 2; RRM2, ribonucleotide reductase M2.

The involvement of 11 significantly expressed genes in cancer drug resistance determined in the COREMINE database, on the basis of literature co-occurrence (Fig. 1), whereas eight genes associated with MDR in ovarian cancer were identified in the IPAD database. Among these, five genes, including aconitase 1 (ACO1), brain-derived neurotrophic factor (BDNF), chemokine (C-X-C motif) receptor 4 (CXCR4), 3-hydroxy-3-methylglutaryl-CoA reductase (HMGCR) and neuropilin 1 (NRP1) were upregulated, and three genes [FAS, cyclin-dependent kinase inhibitor (CDKN)2C) and S-phase kinase-associated protein 2 (SKP2)] were downregulated. In addition, three important signaling pathways were identified in the IPAD database, including the MAPK signaling pathway, the P53 signaling pathway and axon guidance (Table V).

Figure 1

Co-occurrence analysis of significantly expressed genes in multidrug resistance. Red circles represent upregulated genes and grey circles downregulated genes. The lines represent the co-occurrence between nodes. CXCR4, chemokine (C-X-C motif) receptor 4; CDKN2C, cyclin-dependent kinase inhibitor 2C; BDNF, brain-derived neurotrophic factor; ACO1, aconitase 1; SKP2, S-phase kinase-associated protein 2; PIK3CA, phosphatidylinositol-4,5-bisphosphate 3-kinase (catalytic subunit α); NRP1, neuropilin 1; HMGCR, 3-hydroxy-3-methylglutaryl-CoA reductase; FLNA, filamin Aα.

Table V

Enrichment analysis in the IPAD database.

Signaling pathway	Gene	Expression status	AE	RE	MJI	P
Axon guidance	CXCR4, NRP1	Upregulation	2	35.89	0.21	0.04
Axon guidance	NRP1	Upregulation	1	8.77	0.10	0.15
MAPK signaling pathway	BDNF	Upregulation	1	8.54	0.10	0.15
P53 signaling pathway	FAS	Downregulation	1	57.17	0.17	0.04

AE, absolute enrichment value; RE, relative enrichment value; MJI, mean Jaccard index; MAPK, mitogen-activated protein kinase; CXCR4, chemokine (C-X-C motif) receptor 4; NRP1, neuropilin 1; BDNF, brain-derived neurotrophic factor.

The above-mentioned eight genes (ACO1, BDNF, CXCR4, HMGCR, NRP1, CDKN2C, FAS and SKP2) interacted with the 47,077 target genes from TargetScan and 1,675 target genes from PicTar. The corresponding text file was subsequently converted into a format file (.net), to enable it to be recognized by the Pajek software. Following these steps: File>Network>Read and Draw>Net work>Layout>Circular>using Permutation commands, the integrated gene network based on miRNAs and their target genes (data from TargetScan) was constructed using Pajek software. In the integrated gene network, the upregulated gene, NRP1, was found to represent the most important hub gene (Fig. 2). Only NRP1, targeted by miR-130a and miR-130b, was identified in PicTar.

Figure 2

Integrated gene network. Red nodes represented microRNAs/genes. NRP1 was identified as the most important hub gene in the network, and was targeted by miR-130a and miR-130b (circled). The lines indicate the miRNA interacting with their target genes. miR, microRNA; ACO1, aconitase 1; BDNF, brain-derived neurotrophic factor; HMGCR, 3-hydroxy-3-methylglutaryl-CoA reductase; NRP1, neuropilin 1.

The topological characteristics of the integrated gene network, including DC, BC, CC and CC′ were 1.50, 0, 0.44 and 0, respectively. The BC and CC′ values of each node were 0, which indicated that no clustering phenomenon existed in the whole or local network. The calculated results of DC and CC are shown in Table VI. The topological analysis also demonstrated that the hub gene, NRP1, exhibited a higher DC and CC, reflecting the orthocenter of the integrated gene network.

Table VI

Topological characteristics of the integrated gene network.

Node	DC	CC
let-7i	1	0.342857
miR-125b	2	0.342857
miR-130a	1	0.363636
miR-130b	1	0.363636
miR-141	5	0.631579
miR-200c	2	0.363636
miR-214	1	0.363636
miR-376c	3	0.521739
ACO1	2	0.428571
BDNF	3	0.461538
FAS	4	0.500000
HMGCR	2	0.428571
NRP1	5	0.545455

DC, degree centralization; CC, closeness centralization; miR, microRNA; ACO1, aconitase 1; BDNF, brain-derived neurotrophic factor; HMGCR, 3-hydroxy-3-methylglutaryl-CoA reductase; NRP1, neuropilin 1.

The binding sites between miR-130a and NRP1, and miR-130b and NRP1 were identified in StarBase. The two binding sites were found to be located on chromosome 10: 33466864-33466870, determined by CLIP-Seq datasets in the deepView genome browser (Fig. 3). Axon guidance (Fig. 4) was identified based on miRNA signaling pathway analysis in DIANA-miRPath, which involved the hub NRP1 gene.

Figure 3

Illustrative screen shot from the deepView browser. The binding site for miR-130a and NRP1, and the binding site for miR-130b and NRP1, were located on chromosome 10 at position 33466864-33466870. miR, microRNA; NRP1, neuropilin 1.

Figure 4

Axon guidance pathway involving neuropilin 1. Axon guidance acts as one of multiple developmental events in the mitogen-activated protein kinase signaling pathway (circled) (56).

Discussion

MDR is the predominant form of tumor resistance to chemotherapy, and was first reported by Biedler and Riehm in 1970 (53). MDR in EOC may involve complex interactions among genes and other molecules, including miRNAs, proteins or transcription factors, in numerous biological processes. In organisms, these molecular interactions occur between genes, between miRNAs and genes, and between miRNAs, forming complex networks to regulate gene expression. Following the construction of networks for miRNA target genes and transcription factor target genes in ovarian cancer, a previous study identified miR-16, cyclin E1, CDKN1A and E2F transcription factor 1 as hub genes, all of which may be potential biomarkers for ovarian cancer (54). However, there has been little focus on MDR in EOC on the basis of integrated gene network analysis of miRNAs and their target genes. In the present study, a comprehensive bioinformatics analysis was performed through gene expression, pathway enrichment, co-occurrence, and interaction networks, based on published microarray datasets and miRNAs from public databases. This was performed in order to predict the potential molecular mechanism underlying MDR in EOC. The results of the present study suggested that NRP1 was the hub gene of the network, targeted by miR-130a and miR-130b. Therefore, axon guidance involving NRP1 may be a novel signaling pathway for MDR in EOC.

NRP family members have been an area of intense investigation in the field of cancer research over the last 10 years (55,56). NRP receptors are expressed in several types of tumor and endothelial cell, and interact with various soluble molecules, including vascular endothelial growth factor (VEGF), integrin, c-Met and transforming growth factor receptors to modulate the progression of cancer (57). NRP1 encodes one of two NRPs, which contain specific protein domains, allowing NRP1 to be involved in various signaling pathways, which control cell migration (58). NRPs bind several ligands and various types of co-receptors, including VEGF. VEGF has the ability to promote cancer stemness and renewal by directly affecting cancer stem cells via NRP1 in an autocrine loop, and deletion of NRP1 in normal cells prevents tumor initiation (59). It has been reported that immunoreactivity to NRP1 is observed in the vessels of normal tissue samples adjacent to cancer tissue samples, as well as in 98–100% of carcinoma, and the inhibition of NRP1 signaling results in defective angiogenesis and recapitulated the effects of anti-VEGF treatment (60,61). It was also confirmed that knockdown of NRP1 in regulatory T cells may delay or eliminate oncogeny in mouse models of several types of human cancer (62). Disorders of the expression of NRP1 is widely observed in human cancer, including breast cancer, colorectal cancer and chronic lymphocytic leukemia (63–65). NRP1 in plasmacytoid dendritic cells and T regulatory cells is also reported to be a promising therapeutic target for the treatment of cancer (66).

Previously, the expression of NRP1 was found to be upregulated in EOC, and high expression levels of NRP1 enhanced proliferation and were associated with ovarian malignancy, making NRP1 a potential drug targeting candidate for the treatment of EOC (67). These results were concordant with those of the present study. However, the association between high expression levels of NRP1 and MDR in EOC remains to be elucidated and requires further experimental verification in the future.

The results of the present study also suggested that axon guidance was involved in the alteration of the expression of NRP1. The axon guidance signaling pathway represents a pivotal stage in the development of the neuronal network. Axons are guided along specific signaling pathways by various attractive and repulsive guidance molecules, including netrins, slits, ephrins and semaphorins, a number of which have been implicated in human cancers (68–70). However, the role of axon guidance in MDR in EOC remains to be fully elucidated. The present study hypothesized that axon guidance is important in the formation of MDR in EOC, as it acted as one of multiple developmental events in the MAPK signaling pathway (Fig. 4) (71). The latter mediates cisplatin-induced apoptosis and triggers DNA damage and drug resistance in EOC (72,73). Further investigations are required in order to investigate these hypotheses.

Bioinformatics analysis based on genes and miRNAs has been used for investigations on chemotherapeutic response and prognosis in EOC (74–76). MDR in EOC is often accompanied by alterations in gene expression levels and dysfunction of miRNAs. Alterations in gene expression and dysfunction of miRNAs often affects signaling pathways involved in cell proliferation, adhesion, migration, invasion, apoptosis, drug resistance and survival (77–79). Gene expression analysis is one of the basic methods of bioinformatics analysis, and is regularly used to identify dysregulated genes serving as molecular markers for EOC. As high-throughput gene expression data are available in the GEO, four microarray datasets (GSE41499, GSE33482, GSE15372 and GSE28739) were downloaded for gene expression analyses in the present study, which were performed using an R-based web application called GEO2R with the BH method. In addition, 11 miRNAs were mined from previous literature, all of which were revealed to be associated with MDR in EOC. The miRNA target genes in the present study were mined from TargetScan and PicTar. TargetScan and PicTar represent the first and second generation of miRNA target prediction algorithms, respectively, and have high positive predictive values (40,41). The similarity between the two target prediction algorithms is that the seed sequence of miRNA is complementary with the 3′-UTR of mRNA, characterized by a thermodynamic dimer of miRNA target gene (40,41). Furthermore, PicTar was developed based on the design of the first generation of prediction algorithms, breaking the limitation of cross-species conservation. The above characteristic form the basic of integrated network analysis.

COREMINE is a web tool used for literature mining, which performs automated analysis of titles and abstracts by extracting experimental and theoretical biomedical data to create a gene to gene co-citation network (44). IPAD is a web-based database and tool used for mining gene function based on the enrichment analysis of multiple genomic or proteomic data/sources. The present study used a unique approach with that combined literature co-occurrence and disease enrichment. Pathway enrichment combined with network analysis was a further novel approach used in the present study. This was performed to identify the potential genes involved in drug resistance from putative pathways and manually drawn networks. All theoretical conclusions of from the present study require experimental validation and clinical cohorts in the future.

In conclusion, based on a comprehensive bioinformatics analysis using gene expression analysis, signaling pathway analysis, literature co-occurrence and miRNA-mRNA interaction networks, the present study demonstrated that NRP1 is targeted by miR-130a and miR-130b, and may contribute to MDR in EOC. The binding sites of the miRNAs were found to be located on chromosome 10: 33466864-33466870, and all located in axon guidance, a potential pathway associated with MDR in EOC. To the best of our knowledge, the present study is the first to report a gene function of NRP1 associated with MDR in EOC. The findings provide important information for further experimental investigations on the MDR-associated functions of NRP1 in EOC.