This article is mentioned in:

Introduction

Colorectal carcinoma (CRC) is one of the most common and lethal cancers worldwide. Although a trend of decline was observed in CRC due to screening, it still ranked the third in cancer mortality in the United States regardless of gender disparity in the year of 2014 and the third in cancer morbidity worldwide (1,2).

MicroRNAs (miRNAs) are non-coding RNAs participated in mRNA binding and silencing with approximate length of 22 nucleotides. With revelation of their prominent post-transcriptional regulation of miRNAs, there is accumulative evidence demonstrating that miRNAs play significant roles in numerous biological processes and diseases (3,4), and CRC is no exception. Various miRNAs were reported to be associated with CRC, a part of which are oncogenic miRNAs and the others are tumor suppressors (5,6).

MiR-1-3p, as the first numbered miRNA, undoubtedly correlated with multiple human diseases. MiR-1-3p is famous for its regulatory roles both in skeletal and cardiac muscle tissues, including myogenesis, muscle proliferation and differentiation (7). The remarkable roles of miR-1-3p in oncology are also noteworthy. It is widely proved that miR-1-3p functions as a notable tumor suppressor in various types of cancers such as prostate cancer (8), hepatocellular carcinoma (9), gallbladder cancer (10), non-small-lung cancer (11) and gastric cancer (12). Inverse expression of oncogenic factors and miR-1, including NRF2 (13), EGFR (14), MALAT1 (15), etc., indicates targeting relationships in tumorigenesis and development.

In CRC, the majority of existing researches unanimously agreed the tumor-suppressive function of miR-1-3p, but higher expression of miR-1-3p was still observed in brain metastases CRC than primary CRC (16). While revealing attenuated expression level of miR-1-3p and its clinical value in CRC, researchers were also exploring the potential pathways and coefficient target genes related to miR-1-3p. For instance, Wu et al (17), stressed that miR-1-3p could function as an effective diagnostic biomarker of CRC, with an AUC of 0.806. Furukawa et al (18), discovered that miR-1-3p might inhibit the migration of CRC cells by directly targeting NOTCH3, the latter can promote tumorigenesis and migration of CRC. However, despite of the promising regulatory roles miR-1-3p plays in CRC, the application of bioinformatics databases to further validate its function and potential clinical application is yet to be completed. Although several studies have found a part of target genes and pathways that associated with miR-1-3p, a more comprehensive map which extracts data in public databases is essential. Most importantly, since the inconformity in expression arose, a confirming meta-analysis is more persuasive to clarify the characteristics of miR-1-3p.

Hence, we demonstrated the clinical value of miR-1-3p in CRC by extracting information from Gene Expression Omnibus (GEO), ArrayExpress and existing literature, combining microarray data with previous studies. Furthermore, we attained the potential targets of miR-1-3p by gaining intersection of datasets that transfected miR-1-3p into CRC cells, online prediction databases and differentially expressed genes (DEGs) in The Cancer Genome Atlas (TCGA), and added confirmed targets in literature. Subsequently we conducted bioinformatics analysis based on aforementioned selected target genes to unravel the molecular mechanisms of miR-1-3p in CRC.

Materials and methods

Expression and diagnostic data of MiR-1-3p in CRC based on literature and microarray from GEO and ArrayExpress

The design of our investigation was shown in Fig. 1. We searched GEO and ArrayExpress databases with the following search strategy: (gut OR intestinal OR colorectal OR colonic OR rectal OR colon OR rectum OR colon OR colonic OR rectal OR rectum) AND (malignan* OR cancer OR tumor OR tumour OR neoplas* OR carcinoma). Entry type was filtered by ‘Series’ and the organism was restricted to Homo sapiens. The gene chips covering the miR-1-3p level between CRC and non-tumor controls were included for further analysis. The gene chips whose precision was less than 14 decimals were excluded. Samples with other kinds of cancer were eliminated (Fig. 2).

Figure 1. Flow chart of our investigation design.

Figure 2. Flow chart of retrieval in Gene Expression Omnibus (GEO) database.

As for literature, we retrieved PubMed, Science Direct, Google Scholar, Ovid, Wiley Online Library, EMBASE, Web of Science, Chong Qing VIP, CNKI, Wan Fang and China Biology Medicine Disc with the same searching strategy. In the aspect of expression data, studies which we could access the mean, standard deviation and the case numbers of CRC and non-tumorous group were included (Fig. 3). With regards of diagnostic value, we included studies that offered true positive (TP), false positive (FP), false negative (FN) and true negative (TN).

Figure 3. Flow chart of retrieval in online publications databases.

Statistical analysis

The data of each individual gene chip was converted to log2 scale. The number, the mean and the standard deviation of control group and experimental group were calculated based on each single gene chip.

We utilized Stata 14 (StataCorp LLC, College Station, Texas, USA) to display the expression level of miR-1-3p between neoplastic tissues and non-tumorous tissues in forest plots. Standardized mean difference (SMD) and 95% confidence interval (Cl) were calculated to indicate the expression difference. I2>50% or P<0.05 was considered as heterogeneous. Fixed effect model was applied at first and then random effect model was applied when evident heterogeneity still remained. In order to test the reliability of our analysis, we also assessed the publication bias of the study by funnel plots, in which better symmetry suggested less publication bias. Sensitivity analysis was carried out to explore the source of heterogeneity.

IBM SPSS v23 (IBM Corp., Armonk, NY, USA) was applied to perform diagnostic trial based on individual gene chip to acquire TP, FP, FN and TN. We carried out diagnostic meta-analysis using Stata 14. Summarized receiver operating characteristic curve (SROC) were conducted.

Collection of GEO data and literature about intervening experiments

In order to find the promising target genes of miR-1-3p, we searched GEO database again about the experiments which intervened CRC samples or cell lines with miR-1-3p by knockout or transfection, with the following search strategy: (malignan* OR cancer OR tumor OR tumour OR neoplas* OR carcinoma) AND (miR-1 OR miRNA-1 OR microRNA-1 OR miR1 OR miRNA1 OR microRNA1 OR ‘miR 1’ OR ‘miRNA 1’ OR ‘microRNA 1’ OR miR-1-3p OR miRNA-1-3p OR microRNA-1-3p OR miR-1-1 OR miR-1-2 OR miR1-1 OR miR1-2). The search range was set as ‘All Fields’. Gene chips after miR-1-3p intervention which contained mRNAs expression were included.

Meanwhile, in order to collect the target genes of miR-1-3p, we searched PubMed, Science Direct, Google Scholar, Ovid, Wiley Online Library, EMBASE, Web of Science, Chong Qing VIP, CNKI, Wan Fang and China Biology Medicine Disc using the same key words as in GEO database. The target genes of miR-1-3p which were verified via experiments were correspondingly recorded.

Obtainment of prospective target genes of MiR-1-3p

In GSE29760, we considered it significant when log2FC of down-regulated genes was less than −1.

We downloaded data at the entry of colorectal cancer from TCGA. There are 459 patients in total contributing 480 tumor samples and 41 non-neoplasm samples. DESeq data R package was performed to analyze the DEGs. The up-regulated genes with log2FC>1 and P<0.05 were considered significant.

As for prediction tools, miRWalk 2.0 was applied to predict the potential target genes of miR-1-3p based on 12 online databases, namely miRWalk; Microt4; miRanda; miRBridge; miRDB; miRMap; miRNAMap; PicTar; PITA; RNA22; RNAhybrid and Targetscan. For the purpose of acquiring the possible target genes of miR-1-3p as complete as possible, genes were selected on condition that they were recorded no less than 2 prediction databases. The genes selected above were cross-referenced with the DEGs in GEO and TCGA. We also added the confirmed target genes of miR-1-3p in CRC from publication databases. The final genes were considered possible target genes of miR-1-3p.

Analyses in silico via DAVID, cytoscape and STRING

We performed Gene Ontology (GO) analysis via DAVID 6.8 (https://david-d.ncifcrf.gov/) based upon aforementioned prospective target genes. The processed GO categories were biological process (BP), cellular component (CC) as well as molecular function (MF), by which we could further understand the enrichment of the prospective target genes. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis was applied to explore the significant pathways based on the prospective target genes. BiNGO and EnrichmentMap plug-in components in Cytoscape version 3.5.0 were utilized to visualize the network of GO terms and KEGG pathways. Protein-protein interaction (PPI) network was built by STRING (http://www.string-db.org) to show the connections among the possible target genes of miR-1-3p. Disconnected nodes were hidden in the network.

Results

Expression of MiR-1-3p in GEO

There were 7 gene chips from GEO and 1 from ArrayExpress included, with 589 cases in total. No literature met our inclusion criteria. The meta-analysis showed the down-regulation of miR-1-3p in CRC compared to control group. The combined SMD was −0.51 with 95% confidence interval (CI) of −0.68 to −0.33 and I2=88.0% P=0.000 using fixed effect model (Fig. 4A). In random effect model, the SMD changed to −0.63 with 95%CI of −1.19 to −0.07 and I2=88.0% P=0.000 (Fig. 4B). The SMD in tissue subgroup was −0.67 (95%CI: −0.89, −0.46) and I2=89.4% (Fig. 4C). Funnel plot was carried out to test publication bias of our study (Fig. 4D, P>0.05). Sensitivity analysis showed that GSE28364 and GSE71008 might be the source of heterogeneity (Fig. 4E). The SMD after excluding these two gene chips was −1.06 (95%CI: −1.29, −0.82) and I2=43.8% (Fig. 4F).

Figure 4. Expression of miR-1-3p between colorectal carcinoma (CRC) and non-tumorous group. (A) forest plot with fixed effect model; (B) forest plot with random effect model; (C) subgroup analysis of tissue; (D) funnel plot. (E) sensitivity analysis; (F) forest plot after removal of GSE28364 and GSE71008.

Diagnostic value of MiR-1-3p based on GEO gene chips and literature

There were 10 studies included, comprising of 9 microarray datasets and 1 publication from Wu et al (17), containing 645 cases altogether. The combined sensitivity, specificity, positive likelihood ratio, negative likelihood ratio, diagnostic score and odds ratio were 0.74 (95%CI: 0.48, 0.90), 0.75 (95%CI: 0.35, 0.94), 2.94 (95%CI: 1.01, 8.55), 0.34 (95%CI: 0.19, 0.60), 2.15 (95%CI: 1.06, 3.23) and 8.57 (95%CI: 2.89, 25.36) (Fig. 5). The SROC showed the area under the curve (AUC) was 0.81 (Fig. 6). Deek's plot showed P=0.93, which suggested there was no publication bias (Fig. 6D). We also plotted Bivariate Boxplot (Fig. 6A), Likelihood Matrix (Fig. 6C) and Fagan plot (Fig. 7). To summarize, the down-regulation of miR-1-3p showed diagnostic potential.

Figure 5. Forest plots to display the (A) combined sensitivity, specificity, positive likelihood ratio, (B) negative likelihood ratio, (C) diagnostic score and odds ratio.

Figure 6. (A) Bivariate Boxplot, (B) SROC curve, (C) Likelihood Matrix and (D) Deek's plot.

Figure 7. Fagan plot.

Selection of possible target genes of MiR-1-3p

There were 29,834 genes downloaded from miRWalk2.0 with filtration of at least 2 sum scores. In GSE29760, 759 were selected. In TCGA, there were 2111 DEGs. We obtained 24 overlapping genes (NEB, SERPIND1, IL1RAP, SPARC, GNG4, FOXA2, SLC25A22, CSF2, F5, SULT1C2, C2orf61, CREG2, ODAM, PPM1H, TDGF1, E2F5, SH3TC2, SDR16C5, PF4, IL20RA, GPR143, HOXC11, IL1A and FAM57A) (Fig. 8). Then we added the 6 confirmed target genes of miR-1-3p in CRC from literature (NOTCH3, LIM, LASP1, PIK3CA, TWIST1 and GATA4).

Figure 8. Venn of promising target genes of miR-1-3p.

Bioinformatics analysis of the potential target genes of MiR-1-3p

With regard to GO annotation in DAVID, the target genes were closely related to positive regulation of transcription from RNA polymerase II promoter in biological process (BP) (Figs. 9 and 10, P=0.01), extracellular region in cellular component (CC) (Fig. 9, P=3.46E-04), transcription factor binding in molecule function (MF) (Figs. 9 and 11, P=0.01). Moreover, KEGG pathway analysis (Figs. 9 and 12) revealed that cytokine-cytokine receptor interaction (P=9.20E-04), thyroid hormone signaling pathway (P=0.02) and JAK-STAT signaling pathway (P=0.03) were significant. PPI network (Fig. 13) was comprised of GNG4, SPARC, PF4, F5, SERPIND1, IL20RA, CSF2, IL1A and IL1RAP, with enrichment P-value of 0.000427.

Figure 9. The content of Gene Ontology (GO) analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis.

Figure 10. GO biological process (BP) network.

Figure 11. GO cellular component (CC) network.

Figure 12. KEGG pathway network.

Figure 13. Protein-protein interaction (PPI) network.

Discussion

In our present study, we indicated lower expression level of miR-1-3p in CRC than in non-cancerous tissues, which indicated moderate diagnostic accuracy with AUC of 0.81. Bioinformatics analyses based on 30 promising target genes highlighted positive regulation of transcription from RNA polymerase II promoter, extracellular region and transcription factor binding in GO analysis; cytokine-cytokine receptor interaction in KEGG pathway analysis; PF4 in PPI analysis.

At present, the down-regulation of miR-1-3p in CRC and other types of cancer has been observed extensively. Nevertheless, exceptional circumstance still exists. As mentioned in Introduction, we noticed Li et al (16), indicated over-expression of miR-1-3p in CRC with brain metastasis. Because only 4 patients were enrolled in this study, more experiments both in vitro and in vivo are required to verify that whether miR-1-3p possesses contrary functions in CRC. Concerning how to collect data and demonstrate the reduced expression, most studies utilized microarray profiling and quantitative real-time polymerase chain reaction (qRT-PCR) to extract expression profile from tissue samples. However, the sample collected by a single research is limited, and there is no study yet concentrating on collecting existing cases together and excavating the massive materials in public gene expression databases. Thus the meta-analysis we conducted containing literature and microarray data is more reliable in proving the weakened expression of miR-1-3p in CRC. Combining results of diagnostic tests, we can deduct that the down-regulation of miR-1-3p might serve as a more powerful diagnostic target when grouping together with other possible biomarkers. The suppressed expression also presented that miR-1-3p is highly possible to function as a tumor suppressor.

In terms of GO analysis, positive regulation of transcription from RNA polymerase II (Pol II) promoter was listed as the most pivotal biological process. Pol II participated in catalysis of DNA transcription, and studies have found the degradation of Pol II in cancer cells induced by drugs could lead to cytotoxic effects (19). In CRC with microsatellite instability, TAF7, which indirectly impedes Pol II catalyzing by inhibiting the necessary transcriptional factor for Pol II to initiate, was found mutated (20). As a result, we infer that by targeting the genes which can promote the transcription of Pol II promoter, miR-1-3p could negatively regulate CRC proliferation and perform antitumor efficacy.

JAK-STAT signaling pathway ranked the third significant pathway in KEGG pathway analysis. Associated with malignancies and immune diseases, JAK-STAT signaling pathway assists cells to respond to cytokines (21,22), which coincides with ‘cytokine-cytokine receptor interaction’, the most significant signaling pathway in our analysis. In CRC, Slattery et al (23), proved JAK-STAT signaling pathway was intimately connected to CRC risk and survival in their case-control studies. Wang et al (24) blocked JAK-STAT signaling pathway by intervening its up-stream modulator, which contributed to inhibiting CRC progression and improving CRC survival. In addition, non-coding RNAs, especially miRNAs, were reported to have regulatory associations with JAK-STAT signaling pathway (25). For instance, miR-23a and miR-23b was discovered to target at JAK-STAT signaling pathway in prostate cancer by decreasing the expression of IL-6R (26). Therefore, it is deducted that miR-1-3p might be a roadblock in JAK-STAT signaling pathway by silencing the target genes involved in it, and confront CRC by inducing apoptosis and reduce proliferation.

PF4 was the hub gene with most interactions with other genes in PPI network. It mainly participates in blood coagulating, wound repair and inflammatory (27). With regard of neoplasms, PF4 was observed having different roles in various cancers. In animal breast cancer model, PF4 could induce cancer cell apoptosis with the assistance of rapamycin (28). However, contrary results appeared in lung cancer. Pucci et al (29) considered PF4 as an oncogene for it bettered the microenvironment of tumor and promoted its proliferation. Jian et al suggested that PF4 was a protective factor in lung cancer to impede metastasis (30). In CRC, three types of results emerged. Abbasciano et al (31), didn't notice any different expression of PF4 in large bowel carcinoma; Peterson et al (32), claimed PF4 was up-regulated in CRC and its high expression owned diagnostic potential; while Maione et al (33), applied an analogue of PF4 to inhibit CRC development. Apparently, we still have vast space to explore the functions and mechanisms of PF4, in malignancies particularly.

Although we have delineate the map of comprehensive expression and functional annotation with existing data, more experiments were still needed to further validate and better apply the clinical value of miR-1-3p in CRC. Clinicopathological variables that are statistically related to miR-1-3p expression level and survival statistics also require summarizing and analyzing to judge whether miR-1-3p is related to prognosis and other indexes in CRC.

To summarize, miR-1-3p is down-regulated in CRC and is likely to suppress CRC via multiple biological approaches, which indicated diagnostic potential and tumor suppressor efficacy.

References