Identification of Key Circulating Exosomal microRNAs in Gastric Cancer

Exosomal miRNAs (EmiRs) can be used for prediction of gastric cancer (GC) development. Supposedly, both plasma and urinary microRNAs can also be potential biomarkers for screening, but the diagnostic values of EmiRs in blood and urine are not fully studied. We here collected both types of samples from GC patients and healthy individuals and conducted miRNA sequencing to identify key members of EmiRs in GC. The exosomes samples derived from blood and urine were collected from 3 healthy individuals and 7 GC patients. Differentially expressed miRNAs (DEmiRNAs) were acquired, ontology enrichment analysis and Protein-protein Interaction (PPI) enrichment analysis were performed. There were 8 DEmiRNAs in the serum and 3 DEmiRNAs in the urine. For GC patients, there were three up-regulated DEmiRNAs (hsa-miR-130b-3p, hsa-miR-151a-3p and hsa-miR-15b-3p) in the serum exosomes, and one up-regulated DEmiRNA (hsa-miR-1246) in the urinary exosomes. Using miRNA target prediction databases, we found 418 common targets of hsa-miR-15b-3p, 35 common targets of hsa-miR-151a-3p, 117 common targets of hsa-miR-130b-3p, and 357 common targets of hsa-miR-1246. Some commonly enriched ontology terms were found, including GO BP terms like cell surface receptor signaling pathway involved in cell-cell signaling, positive regulation of catabolic process, morphogenesis of an epithelium, and GO CC terms perinuclear region of cytoplasm. The PPI network show some key nodes, including TAOK1, CMTM6, SCN3A, WASF3, IGF1, CNOT7, GABRG1, PRKD1. Together, this study provided an integrative analysis of expression profile of key circulating exosomal microRNAs. Four key exosomal miRNAs (hsa-miR-130b-3p, hsa-miR-151a-3p and hsa-miR-15b-3p) and the interaction network or enrichments based on their targets (TAOK1, CMTM6, SCN3A, WASF3, IGF1, CNOT7, GABRG1, PRKD1) may provide a reference of the molecular mechanisms in the GC development.


INTRODUCTION
Gastric cancer (GC) is the fourth most common malignance worldwide (1). An exact detection based on GC biomarkers has a high clinical significance. For GC patients, the exosomal microRNA (miRNA or miR) may have potential clinical application value in diagnosis and status evaluation (2,3). MicroRNA is a short-chain noncoding RNA molecule (around 22 nucleotides in length). It regulates protein expression of a particular mRNA (by incomplete base pairing), causing inhibition of protein translation of its target genes. Exosomes are a subset of extracellular vesicles containing a variety of bioactive molecules including proteins, lipids, and RNAs. It can play a role in intercellular crosstalk. During tumor development, exosomes can facilitate construction of the specific microenvironment; besides, it strongly affects the occurrence, metastasis, and even drug resistance. Tumor secreted exosomes are sharply different among cancers, and between healthy individuals vs. cancer patients. Theoretically, the contents of exosomes include different tumor-related biomarkers. Commonly, previous studies observed the samples of plasma exosomal miRNAs (EmiRs) for prediction of GC development, progression, and treatment outcomes. Although the development of new biomarkers in blood tests has shown great potential, limited reports have applied multi-biomarkers based on the strategy of exosomal miRNA detection. Supposedly, urinary microRNA can also be a potential biomarker for tumor screening. For example, a lower urine level of miR-30a-5p was found in gastric cancer and colon carcinoma patients when compared to ovarian serous adenocarcinoma, and urinary miR-30a-5p from ovarian cancer patients was notably reduced following the surgical removal of ovarian serous adenocarcinoma (4). However, the diagnostic value of EmiRs in urine for GC is unknown, and that of plasma EmiRs is also not fully studied. We here collected both plasma samples and urine samples from GC patients and healthy control individuals for miRNA sequencing, and some key members were identified.

Exosomal miRNA Sequencing
The exosomes samples derived from blood and urine were collected from 3 healthy individuals and 7 GC patients (Stage I-II), including 6 blood samples (3 GC vs 3 controls) and 10 urine samples (7 GC vs 3 controls). The median age of the 7 patients was 52 years, and that of the healthy controls was 45 years. All GC patients were confirmed by histopathology. The frozen samples were thawed and centrifuged at 2000 g for 30 min to remove cells/debris. Next, 0.2 volumes of the Total Exosome Isolation reagent (Thermo Fisher, Inc.) were added, the solution was thoroughly mixed and incubated for 30 min at 4°C. Next, the samples were centrifuged at 10000 g for 10 min and the exosome pellet was resuspended in PBS. Next, the exosome-derived RNA was extracted from exosomes using the Total Exosome RNA and Protein Isolation Kit, the total amount of RNA was detected. Subsequently, the Illumina NextSeq 500 SE50 (20M) sequencing was performed. After sequencing, the FastQ-format data were stored and cleaned for further analysis.

Differential miRNAs Analysis
First, miRNA expression was analyzed by miRExpress (http:// mirexpress.mbc.nctu.edu.tw), which works based on the miRNA database (http://www.mirbase.org) for quantitative analysis of miRNA. After the second-generation sequencing raw data in FastQ format are imported, the software compares the data with the known miRNA sequence in miRBase and then acquires the read count of the same sequence fragments as Counts, which is used to measure the miRNA levels in a sample. According to the grouping information, the sample expression data was imported into the edgeR R package to compare the expression differences between samples. The Fold Change and p-value were obtained. Then, according to the threshold of p-value<0.05 and Log (Foldchange) >2 or <-2, the DEmiRNAs was acquired. The volcanic plots were produced to present DEmiRNAs using the ggplot2 R package. And a principal component analysis (PCA) was performed using the ggfortify R package based on all detected miRNA levels.

Expression Analysis of DEmiRNAs in TCGA and EVmiRNA
Exosomal miRNAs that are up-regulated in the urine or blood exosomes of tumor patients have potential clinical applications, so the differential miRNAs obtained from screening were subjected to expression analysis. From the Cancer Genome Atlas (https://portal.gdc.cancer.gov/) and EVmiRNA (http:// bioinfo.life.hust.edu.cn/EVmiRNA#!/) databases, the raw count of tumor and exosomal miRNA expression data were obtained. And box plots were drawn by ggplot2 R software.

Prognostic Analysis of DEmiRNAs in TCGA
Corresponding clinical information were obtained from The TCGA dataset. High and low expression of miRNAs was defined with the median. For the Kaplan-Meier curve, the p value and the hazard ratio (HR) with 95% confidence interval (CI) are obtained by logrank test and univariate Cox proportional hazard regression. p<0.05 was considered statistically significant. And p < 0.05 was considered as statistically significant.

Target Prediction of the DEmiRNAs
First, we focused on the up-regulated DEmiRNAs which may be encapsulated in circulating exosomes. The mRNAs potentially targeted by these miRNAs were predicted using both miRDB and Targetscan databases. When the intersection of two databases regarding the targets of each up-regulated DE-miR was acquired, the overlap set of mRNAs were further used for enrichment analysis.

Predicted Target Gene Enrichment Analysis
The Metascape gene list analysis online tool was used for enrichment analysis of the predicted target genes. All targets were imported into the tool, gene annotation was performed, and pathway and process enrichment analysis were carried out with the following ontology sources: KEGG Pathway, GO Biological Process, GO Molecular Function and GO Cell Component. The Homo Sapiens background was used for enrichment analysis of all predicted target genes. Terms with a p-value < 0.01, a minimum count of 3, and an enrichment factor > 1.5 (the ratio between the observed counts and the counts expected by chance) were collected and grouped into clusters based on their according DEmiRNAs. The statistically enriched terms were presented by accumulative hypergeometric p-values and enrichment factors, and then hierarchically clustered into a tree based on Kappa-statistical similarities.

Protein-Protein Interaction Analysis
Protein-protein interaction analysis was carried out with STRING and BioGrid. And the Molecular Complex Detection (MCODE) was applied to identify densely connected network components. The MCODE networks identified for individual gene lists were gathered and drawn by the Metascape (https:// metascape.org/).

Network Hub Gene Expression Analysis
For the analysis of hub genes in the resulting MCODE subnetwork, GEPIA (http://gepia.cancer-pku.cn/) tools were used to analyze their expression profiles in the STAD from TCGA. And p < 0.05 was considered as statistically significant.

Identification of GC Associated DEmiRNAs
Principal component analysis were presented in Figure 1, which suggested that there is some heterogeneity in the expression profiles of blood and urine exosomal miRNAs in patients. As Table 1 and Figure 2 shown, there were 8 DEmiRNAs in the serum and 3 DEmiRNAs in the urine. The expression profiles of all differentially expressed miRNAs are shown in the Heatmap. For GC patients, there were 3 up-regulated DEmiRNAs (hsa-miR-130b-3p, hsa-miR-151a-3p and hsa-miR-15b-3p) in the serum exosomes, and one DEmiRNA (hsa-miR-1246) enriched in the urinary exosomes. Theoretically, these up-regulated miRNAs may be transferred to tissues or cells and target specific RNAs, which finally promotes the GC development.

Prognostic of DEmiRNAs in TCGA
As shown in Figure 4A, the results of univariate Cox analysis regarding these four miRNAs revealed that only one serum exosomal hsa-mir-15b-3p was identified to be associated with GC patient prognosis (HR = 0.814, P < 0.05). None of the other three were related, hsa-mir-130b-3p (HR = 0.977, P > 0.05), hsamir-151a-3p (HR = 0.835, P > 0.05), hsa-mir-1246 (HR = 0.804, P > 0.05). While the results of prognostic analysis in Figure 4B were different, high hsa-mir-151a-3p expression was considered to be associated with poor survival outcomes in GC patients (P < 0.05), while the other three miRNAs were not.

Predicted Target Gene Enrichment Analysis
As Figure 6 and Table 2 shown, the enrichment analysis was conducted based on above targets. Some commonly enriched ontology terms were found, including GO BP terms like cell surface receptor signaling pathway involved in cell-cell signaling, positive regulation of catabolic process, morphogenesis of an epithelium, and GO CC terms perinuclear region of cytoplasm.

PPI Network
Applying the STRING and BioGrid data, the PPI network was generated. Meanwhile, each MCODE network was drawn assigned by a unique color ( Figures 7B, D) or by the matched DEmiRNAs ( Figures 7A, C). Some genes were located at key nodes, such as TAOK1, CMTM6, SCN3A, WASF3, IGF1,    CNOT7, GABRG1 and PRKD1. And the information of these key genes from mcode screening was exhibited in Table S1.

Network Hub Gene Expression Analysis
As shown in Figure 8, the expression levels of eight key genes screened by MCODE in STAD tissues of TCGA were analyzed. The results showed that the expression of TAOK1, CMTM6 and CNOT7 were significantly lower in the adjacent cancers, while the expression of WASF3 was significantly higher in the adjacent cancers. SCN3A, IGF1, GABRG1 and PRKD1 were not significantly expressed.
There is a study has reported that miR-130b-3p was upregulated in GC tissues, and miR-130b-3p promoted survival, metastasis and angiogenesis of GC cells as well as enhanced tumor formation and angiogenesis in GC in vivo (14). Another one has identified exo-miR-15b-3p/DYNLT1/ Caspase-3/Caspase-9 axis do promotes GC development and malignant transformation. It means that serum exo-miR-15b-3p may be a potential GC diagnosis and prognosis biomarker, which  can be used in precise targeted GC therapy. Studies in vitro revealed that elevated serum miR-1246 was tumor-derived by being packaged into exosomes with the help of ELAVL1 (15). And in a other area, exosomal miR-1246 expressions in serum could differentiate GC patients with TNM stage I from healthy controls and patients with benign diseases (11). Furthermore, there are no correlation between hsa-miR-130b-3p, hsa-miR-151a-3p, hsa-miR-15b-3p and hsa-miR-1246 expression and   A literature search of important hub genes revealed that CKLF like Marvel transmembrane domain 6 (CMTM6) was involved in gene epigenetic regulation and tumorigenesis, and the combined cmtm6 and PD-L1 testing could be used as an indicator to determine the prognosis of patients with gastric cancer (16,17). Whereas Wiskott Aldrich syndrome protein family member 3 (WASF3) is required for tumor invasion and metastasis. Reports have indicated that WASF3 expression is associated with poor prognosis and is a potential prognostic factor for gastric cancer patients, therefore, targeting WASF3 is a novel potential therapeutic strategy for gastric cancer (18,19). In addition to this, it has also been documented that mir-218 can inhibit the proliferation, migration and EMT of gastric cancer cells SGC7901 by targeting WASF3 (20). Protein kinase D3 (PRKD3) promotes cancer cell proliferation, growth, migration, and invasion in various tumor types. A growing body of data supports that PRKD3 is a promising therapeutic target for the treatment of cancer (21).
In the PPI network, most of the key nodes were promising candidates that have been fully surveyed, such as TAOK1, CMTM6, SCN3A, WASF3, IGF1, CNOT7, GABRG1 and PRKD1.
Still, this work has some limitations. First, the published studies about hsa-miR-130b-3p, hsa-miR-151a-3p, hsa-miR-15b-3p and hsa-miR-1246 are very few, and known conclusions implied that they are tumor suppressors in all possibility. Again, many targets of our up-regulated DEmiRNAs have been reported to be oncogenic factors or increased in GC (as mentioned above). Their conclusions sharply contradict ours, and the deep reasons need to be further explored. Second, we failed to acquire any common DE-miR in blood and urine, which may be due to that only 3 or 7 samples were included in each group. Theoretically, there may be some enriched circulating exosomal miRNAs in both blood and urine, and more samples should be accumulated in future.

CONCLUSIONS
In conclusion, this study provides an integrative analysis of expression profile of key circulating exosomal microRNAs. Four key miRNAs were found: hsa-miR-130b-3p, hsa-miR-151a-3p, hsa-miR-15b-3p and hsa-miR-1246. These risk exosomal microRNAs, as well as the corresponding interaction network or enrichments based on their targets (such as TAOK1, CMTM6, SCN3A, WASF3, IGF1, CNOT7, GABRG1, PRKD1) may provide a better understanding of the molecular mechanisms in the GC development.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by The Ethics Committee of Renji Hospital (Shanghai, China). The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
Conception and design: XQ. Administrative support: DC. Provision of study materials or patients: HW. Collection and assembly of data: FX. Data analysis and interpretation: XQ.
Manuscript writing: All authors. All authors contributed to the article and approved the submitted version.