Identification of cell cycle as the critical pathway modulated by exosome-derived microRNAs in gallbladder carcinoma

Gallbladder cancer (GBC), the most common malignancy in the biliary tract, is highly lethal malignant due to seldomly specific symptoms in the early stage of GBC. This study aimed to identify exosome-derived miRNAs mediated competing endogenous RNAs (ceRNA) participant in GBC tumorigenesis. A total of 159 differentially expressed miRNAs (DEMs) was identified as exosome-derived miRNAs, contains 34 upregulated exo-DEMs and 125 downregulated exo-DEMs based on the expression profiles in GBC clinical samples downloaded from the Gene Expression Omnibus database with the R package. Among them, 2 up-regulated exo-DEMs, hsa-miR-125a-3p and hsa-miR-4647, and 5 down-regulated exo-DEMs, including hsa-miR-29c-5p, hsa-miR-145a-5p, hsa-miR-192-5p, hsa-miR-194-5p, and hsa-miR-338-3p, were associated with the survival of GBC patients. Results of the gene set enrichment analysis showed that the cell cycle-related pathways were activated in GBC tumor tissues, mainly including cell cycle, M phase, and cell cycle checkpoints. Furthermore, the dysregulated ceRNA network was constructed based on the lncRNA-miRNA-mRNA interactions using miRDB, TargetScan, miRTarBase, miRcode, and starBase v2.0., consisting of 27 lncRNAs, 6 prognostic exo-DEMs, and 176 mRNAs. Together with prognostic exo-DEMs, the STEAP3-AS1/hsa-miR-192-5p/MAD2L1 axis was identified, suggesting lncRNA STEAP3-AS1, might as a sponge of exosome-derived hsa-miR-192-5p, modulates cell cycle progression via affecting MAD2L1 expression in GBC tumorigenesis. In addition, the biological functions of genes in the ceRNA network were also annotated by Gene Ontology and Kyoto Encyclopedia of Genes and Genomes. Our study promotes exploration of the molecular mechanisms associated with tumorigenesis and provide potential targets for GBC diagnosis and treatment. Supplementary Information The online version contains supplementary material available at 10.1007/s12032-021-01594-8.

Exosomes, a type of extracellular vesicles (EVs) with diameters ranging from 40 to 100 nm, are widely released from many cell types into the extracellular space. Recently, mRNAs, lncRNAs, and miRNAs have been identified in exosomes. Dysregulation of exosome-derived ncRNAs is associated with tumor invasion, migration and metastasis, and drug resistance [22][23][24]. However, there are no studies specifically focused on the exosome-derived ncRNAs and their regulation of ceRNA networks in GBC. In the present study, we conducted a multi-step analysis using various R language packages on the expression profiles of clinical samples downloaded from the Gene Expression Omnibus (GEO) database to identify the differentially expressed mRNAs and ncRNAs in the GBC. The exosome-derived miRNAs associated with GBC patients' survival were identified, and the related ceRNA network was established based on the lncRNA-miRNA-mRNA interaction. The gene set enrichment analysis (GSEA) of the differentially expressed genes regulated by exosome-derived miRNAs was identified to investigate the downstream regulation role of the exosomederived miRNAs.

Differential expression analysis
The differentially expressed genes (DEGs), lncRNAs (DELs), and miRNAs (DEMs) between GBC tissues and normal tissues of the GEO datasets, were calculated using the "limma" package with voom method in R, respectively [25]. The False Discovery Rate (FDR)-adjusted P value < 0.05 by the Benjamini-Hochberg method and |log2 fold change (FC)|> 2 were set as cut-off criteria for DEMs and DELs, while as adjusted P value < 0.05 and |log2FC|> 1.5 were for DEGs. Visualization of the identified DEGs and DEMs including volcano plot and heatmap were performed with the "ggplot2" and "pheatmap" packages of R, respectively [26]. Then, DEMs and exosome miRNAs, which were obtained from the exosome databases, ExoR-Base and EVmiRNA, were intersected using the "VennDiagram" package of R, and the overlapped genes were identified as the exo-DEMs. The GBC tissues from GSE104165 were divided into the long survival group and short survival group based on the patients' survival. Then, the expression profiles of exo-DEMs were analyzed between the two groups, and the survival-related exo-DEMs were selected for subsequent analysis.

GSEA analysis of the survival-related exo-DEMs target genes
To elucidate the molecular mechanism of the exo-DEMs associated with survival, the interactions of miRNA-target genes were obtained from the reliable online miRNA-mRNA databases, including miRDB, TargetScan, and miRTar-Base, using the "multiMiR" package in R (http:// multi mir. org/). Interactions of miRNA-target genes were obtained based on experimental verification of luciferase reporter assay. Then, the differential expression analysis of these target genes between GBC tissues and normal tissues from GSE74048 was performed, and the possible pathways and functions were predicted using the GSEA method with the Kyoto Encyclopedia of Genes and Genomes (KEGG) and Reactome databases. NOM P value < 0.05 was considered statistically significant. The results were visualized with the "clusterProfile" package in R [27].

Construction of the survival-related exo-DEMs mediated ceRNA network
LncRNA-miRNA interaction pairs were predicted using the online databases, miRcode (http:// www. mirco de. org/ index. php) and starBase v2.0 (http:// starb ase. sysu. edu. cn/ starb ase2/). Then integrated with the miRNA-mRNA interactions to establish a dysregulated lncRNA-miRNA-mRNA ceRNA network and visualized using Cytoscape software. Furthermore, the genes and the lncRNAs in ceRNA networks were intersected with DEGs and DELs, respectively, using Venn diagram, to identify the DEGs and DELs in the ceRNA networks. The gene expression was analyzed in GBC and normal tissue using the paired Student's t-test. The twosided P < 0.05 was considered statistically significant.

Functional enrichment analysis and protein-protein interaction of the ceRNA network
To better understand the biological functions of the integrated ceRNA network, Gene Ontology (GO) covering biological processes (BP), molecular functions (MF), and cellular components (CC), enrichment analyses were performed. The whole human genome was set as the background, and functional categories with adjusted P < 0.01 were considered significant. The results were visualized with the "PerformanceAnalytics" (https:// cran.r-proje ct. org/ web/ packa ges/ Perfo rmanc eAnal ytics/) package of R. Furthermore, the protein interaction network was mapped with the Search Tool for the Retrieval of Interacting Genes/ Proteins (STRING) website (http:// string-db. org/) and visualized using Cytoscape software [28].

Identification of DEMs, DELs, and DEGs in gallbladder carcinoma
A total of 843 DEGs and 895 DELs between GBC tissues and the matched peri-carcinomatous tissues were obtained after analyzing expression profiles from GSE74048 datasets. Among them, 246 mRNAs and 716 lncRNAs were significantly upregulated; while 597 mRNAs and 179 lncRNAs were significantly downregulated, respectively. Meanwhile, 204 DEMs, including 71 up-regulated and 133 down-regulated miRNAs, between GBC tissues and normal tissues were obtained from GSE104165 dataset, respectively. The volcano plots and heatmaps illustrated the significant differences and distribution of the fold change in DEMs, DELs, and DEGs (Fig. 1a, b).

The exosome DEMs associated with the GBC patients' survival
To screen the exosome-related DEMs in GBC compared with normal tissues, Venn diagram analysis was used to obtain the intersection between DEMs and miRNAs from exosomes. Subsequently, a total of 159 overlapped DEMs were identified, containing 34 upregulated exo-DEMs and 125 downregulated exo-DEMs (Fig. 1c). Seven survival associated exo-DEMs were obtained through the differential expression analysis between the long survival group and the short survival group (Fig. 1d), including 2 up-regulated exo-DEMs, hsa-miR-125a-3p and hsa-miR-4647, and 5 downregulated exo-DEMs, hsa-miR-29c-5p, hsa-miR-145a-5p, hsa-miR-192-5p, hsa-miR-194-5p, and hsa-miR-338-3p. Results showed that the 2 up-regulated exo-DEMs were up-regulated in the short survival group, while the 5 downregulated exo-DEMs were up-regulated in the long survival group (Fig. 1e).

GSEA analysis of the survival-related exo-DEMs target genes
To elucidate the molecular mechanism of the exo-DEMs with patients' survival, GSEA analysis of target genes was performed. Result showed that the GBC tumor tissues were mainly enriched in KEGG pathways including cell cycle (hsa04110, P = 0.002), protein processing in endoplasmic reticulum (hsa04141, P = 0.002), and carbon metabolism (hsa01200, P = 0.004), while as the normal tissues were mainly enriched in proteoglycans in cancer (hsa05205, P = 0.001), and serotonergic synapse (hsa04726, P = 0.004) ( Fig. 2a; SI 1). That indicated the genes in cell cycle, protein processing in endoplasmic reticulum, and carbon metabolism were activated in GBC tumor tissues, and the genes in pathways of proteoglycans in cancer, and serotonergic synapse were suppressed (Fig. 2b). Reactome pathway enrichment also indicated that the cell cycle-related pathways were activated in GBC tumor tissues, including cell cycle

Survival-related exo-DEMs mediated ceRNA network
To better understand the biological regulation roles of the exo-DEMs with patients' survival, we constructed the dysregulated ceRNA network based on the lncRNA-miRNA-mRNA interactions. The regulatory relationship between 176 mRNAs and 6 of survival-related exo-DEMs with patients' survival was found based on experimental verification of luciferase reporter assay (SI 3). Then, 27 lncR-NAs were predicted to interact with the exo-DEMs using the online databases, miRcode, and starBase v2.0. Finally, integrated the relationship of miRNAs and mRNAs, miRNAs and lncRNAs, a ceRNA network consisting of 27 lncRNAs, 6 miRNAs, and 176 mRNAs was constructed (Fig. 3a). In addition, we calculated the connection degree of each gene by topology to clarify its importance in the ceRNA network. mRNA (CDH2), lncRNA (KCNQ1OT1), and miRNAs (hsa-miR-145-5p) were considered the most important genes among the lncRNAs, miRNAs, and mRNAs, respectively (SI 4).

Functional enrichment and expression analysis of genes in the ceRNA network
GO enrichment and KEGG pathway analysis were used to investigate the mechanisms associated with the ceRNA network with the threshold set as adjusted P < 0.01. The significantly enriched GO terms were illustrated and available in    Fig. 1 Identification of DEL, exosomal DEMs and DEGs between GBC and normal tissues. a Volcano plot of DEL, DEMs, and DEGs for datasets from GEO. X-axis: log2 fold change; Y-axis: − log 10 (P value) for each gene; vertical-dotted lines: fold change ≥ 2 or ≤ 2 for DELs and DEMs, fold change ≥ 1.5 or ≤ 1.5 for DEGs; horizontal-dotted line: the significance cut off (adjusted P value = 0.05). The red dot represents up-regulated genes, and the blue dot represents down-regulated genes. b Gene expression heatmap of DEL, DEMs, and DEGs, respectively. c Venn diagrams of the overlapping genes between exosome-derived miRNAs with up-regulated DEMs and down-regulated DEMs, respectively. d Box plots that represent the expression levels of exo-DEMs between tumor tissues and normal tissues. e Box plots that represent the expression levels of exo-DEMs between the tissues from long survival and short survival patients. **** indicates P < 0.0001, *** indicates P < 0.001, ** indicates P < 0.01 Protein-protein interaction of the gene in the ceRNA networks shown that the high-scoring proteins were obtained including CD44, MYC, CCND1, ESR1, IL6, and VEGFA, which suggesting that they might jointly regulate the occurrence and development of tumor (Fig. 4). Expression analysis indicated that five lncRNAs and ten mRNAs were differentially expressed between GBC tumor tissues and normal tissues, including seven mRNAs down-regulated in tumor tissues, while the other three mRNAs and all the five lncRNAs up-regulated in tumor tissues (Fig. 5).

Cell cycle as a critical pathway in GBC tumorigenesis
Based on the integrated ceRNA network and GSEA analysis, the STEAP3-AS1/hsa-miR-192-5p/MAD2L1 axis, which participants in cell cycle pathway, was identified as critical in GBC tumorigenesis. Besides, as shown in Fig. 6, we also found that mainly members in this pathway were differentially expressed, including 17 genes were up-regulated in GBC tumor tissues, while three genes were down-regulated.

Discussion
As dysregulation of exosomal ncRNAs has been highlighted for critical functions in tumor invasion, migration and metastasis, and drug resistance. With the purpose of identifying the exosomal miRNAs mediated ceRNA participant in GBC tumorigenesis and to investigate the potential biomarkers for better detection and therapy, we analyzed the gene expression profiles of the GEO datasets and established a lncRNA-miRNA-mRNA network mediated by prognostic exo-DEMs at the transcriptomewide level to provide a useful foundation for the regulatory function research in GBC.
After that, GSEA analysis was performed to investigate the molecular mechanism of the prognostic exo-DEMs. Results indicated that the enriched KEGG and Reactome pathways which activated in GBC tumor tissues were mainly associated with the cell cycle-related pathways including cell cycle, M phase, and cell cycle checkpoints. Furthermore, the dysregulated ceRNA network based on the lncRNA-miRNA-mRNA interactions was constructed, consisting of 27 lncRNAs, 6 prognostic exo-DEMs, and 176 mRNAs. GO enrichment and KEGG pathway analysis also showed that the genes of the ceRNA network were mainly enriched in cell cycle-related GO terms and KEGG pathways. Before the reproductive cycle of replicating DNA, cells enter a phase called G1 during which they release plenty of signals that contribute to cell division and cell fate. miRNAs regulate genes involved in checkpoints, DNA repair, and cell cycle. Dysregulated miRNAs can lead to the loss of cell cycle regulation and finally contribute to pathological situations, such as cancer [38][39][40]. Expression analysis indicated that 5 lncRNAs and 10 mRNAs in the ceRNA network were differentially expressed between GBC tumor tissues and normal tissues. Together with prognostic exo-DEMs, the STEAP3-AS1/ hsa-miR-192-5p/MAD2L1 axis was identified. Due to the main members in the cell cycle pathway were differentially expressed including MAD2L1, we speculated that the STEAP3-AS1/hsa-miR-192-5p/MAD2L1 axis was critical in GBC tumorigenesis through regulation of cell cycle. A recent study has found that STEAP3-AS1 downregulation could increase the expression of cyclin-dependent kinase inhibitor 1C (CDKN1C) by STEAP3 upregulation and modulate the cell cycle progression in CRC [41]. STEAP3-AS1 was also identified as a novel, potential prognostic biomarker, and therapeutic target for tongue squamous cell carcinoma (TSCC) [42]. STEAP3-AS1, together with HOTAIR and SOX21-AS1, as the exo-lncRNA signature, was proven to be an independent prognostic factor in glioblastoma (GBM) [43]. MAD2L1, MAD2 mitotic arrest deficient-like 1, is a critical mitotic checkpoint gene. Overexpressing MAD2L1 has shown the effect on cancer cell proliferation, migration, invasion, and cell cycle arrest [44][45][46]. Activation of miR-192/215 by p53 includes a set of genes which regulate DNA synthesis and G1 and G2 checkpoints, including CDC7, MAD2L1, and CUL5, to induce cell arrest and reduce tumor cell growth [47][48][49][50].
In our study, we found that STEAP3-AS1 and MAD2L1 were up-regulated in GBC tumor tissues, while as hsa-miR-192-5p was down-regulated, suggesting lncRNA STEAP3-AS1, might as a sponge of exosome-derived hsa-miR-192-5p, modulates cell cycle progression via affecting MAD2L1 expression in GBC tumorigenesis. Further experimental studies are required to verify this hypothesis. These findings suggest that the ceRNA network have significant effects on the pathogenesis and prognosis of GBC.

Conclusion
This study provides a better understanding of the lncRNA-miRNA-mRNA regulatory mechanism in GBC mediated by exosome-derived miRNAs. In future research, we will explore the prognostic genes and their potential function of the ceRNA axis based on the present study.

Conflict of interest
The authors declare that they have no competing interests.
Ethical approval Not applicable.

Consent for publication Not applicable.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.