Abstract
Gastric cancer is one of the most common tumors of the digestive system. Here, analysis of the expression profiles of circular RNAs in advanced gastric adenocarcinoma and adjacent normal mucosa tissues revealed differential expression of 306 circular RNAs, among which 273 were predicted to exert regulatory effects on target microRNAs. The downstream pathway networks of circular RNA–microRNA were mapped and the node genes were identified. In particular, we found that the expression of hsa_circ_0058246 was elevated in tumor specimens of patients with poor clinical outcomes. Our collective findings indicate that circular RNAs play a critical role in gastric cancer tumorigenesis. Data from this study provide a new perspective on the molecular pathways underlying metastasis and recurrence of gastric cancer and highlight potential therapeutic targets that may contribute to more effective diagnosis and treatment of the disease.
Introduction
Gastric cancer, one of the most common tumors of the digestive system, is ranked the second leading cause of cancer-related deaths worldwide.1–3 Despite significant advances in early diagnostic techniques, numerous patients with gastric cancer are diagnosed at the advanced stages that develop into perivascular invasion, lymph node metastasis, and peritoneal metastasis.2,4 At present, the mainstay of gastric cancer treatment is radical surgical resection, 5 but the majority of patients require supplementation with chemotherapy, radiotherapy, and other treatment modes.5,6 Despite the range of therapeutic options available, overall prognosis of gastric cancer remains poor. 1 Gastric cancer is a multi-factorial and multi-step process that involves significant molecular participation and complex network regulation.
Recently, circular RNA (circRNA), rather than linear RNA molecules, have emerged as a common phenomenon in the regulation of gene expression in human cells.7,8 These RNA molecules play major roles in the development of several tumor types.9–11 At present, circRNA-related research is still in its infancy, and the functions and regulatory mechanisms of the members of known circRNA libraries are yet to be established. CircRNAs exert distinct effects, both as tumor suppressors and oncogenes.12,13 Hsa_circ_0000069 has been identified as a critical regulator in colorectal cancer, and thus presents a promising target in diagnosis and treatment of the disease. Knockdown of hsa_circ_0000069 inhibits colorectal cancer cell proliferation, migration, and invasion and induces G0/G1 phase arrest of the cell cycle. 10 A recent study showed that hsa_circ_002059 is significantly downregulated in gastric cancer and correlated with metastasis. 14 Expression of related circRNAs in laryngeal, liver, and esophageal cancer types has additionally been reported.15–17
Owing to rapid developments in gene sequencing technology, nucleic acid sequences in the gene database are accumulating considerably and the amount of biological information is expanding rapidly, providing further opportunities for data interpretation and potential breakthrough discoveries.
In this study, we analyzed the expression profiles of circRNAs in advanced gastric adenocarcinoma and adjacent normal mucosa tissues via use of biochips. The functions of differentially expressed circRNAs were further analyzed, and the target microRNAs (miRNAs) regulated by the circRNAs were predicted with a view to elucidating the regulatory pathways, mapping downstream pathway networks of interacting circRNA–miRNAs, and identifying the node genes that play critical roles in disease progression. Especially, we found that the expression of hsa_circ_0058246 was elevated in tumor specimens of patients with poor clinical outcomes. Our findings provide clues for the identification of key targets of posttranscriptional regulation and relevant individualized treatment strategies for gastric cancer.
Methods
Clinical specimens and ethical statement
Three pairs of III stage gastric cancer and paired paracarcinoma tissue specimens were obtained for circRNA microarray analysis. In addition, a total of 43 patients who underwent radical gastrectomy at Fudan University Shanghai Cancer Center between December 2014 and January 2016 were involved in this study. All selected patients were newly treated and did not receive preoperative radiotherapy, chemotherapy, or other treatments. All specimens were confirmed by pathology according to the guidelines of the International Union for Cancer Control/American Joint Committee on Cancer (UICC/AJCC) tumor–node–metastasis (TNM) staging criteria for gastric cancer (2010). Written informed consent was obtained from patients or family members. All the baseline characteristics of patients and patients’ follow-up data are shown in the Table 1.
The baseline characteristics of patients and patients’ follow-up data.
TNM: tumor–node–metastasis.
RNA extraction and purification
Total RNA was extracted and purified using the mirVana™ miRNA Isolation Kit (Ambion, Austin, Texas, USA) following the manufacturer’s instructions, and the RNA integrity number (RIN) was calculated using an Agilent Bioanalyzer 2100 (Agilent Technologies, Santa Clara, CA, USA).
RNA quality control
Total RNA was quantitated with a NanoDrop ND-2000 spectrophotometer and Agilent Bioanalyzer 2100 (Agilent Technologies) and subsequently used for biochip assay. This step was completed with the aid of Shanghai Biotechnology Co., Ltd (Shanghai, China).
RNA amplification and labeling
Ribosomal RNA (rRNA) was amplified and labeled using the Low Input Quick Amp WT Labeling Kit (Agilent Technologies) following the manufacturer’s instructions. Labeled complementary RNA (cRNA) was purified with an RNeasy Mini Kit (QIAGEN, Hilden, Germany). This step was completed with the aid of Shanghai Biotechnology Co., Ltd.
Chip hybridization
Each slide was hybridized with 1.65 µg Cy3-labeled cRNA using a Gene Expression Hybridization Kit (Agilent Technologies) in a hybridization oven (Agilent Technologies) according to the manufacturer’s instructions. After 17 h of hybridization, slides were washed in staining dishes (Thermo Shandon, Waltham, MA, USA) using a Gene Expression Wash Buffer Kit (Agilent Technologies), as recommended by the manufacturer. This step was completed with the aid of Shanghai Biotechnology Co., Ltd.
Data acquisition
Slides were scanned with an Agilent Microarray Scanner (Agilent Technologies) under default settings: dye channel: green, scan resolution = 3 µm, photomultiplier tube (PMT) gain: 100%, and dynamic range: 20 bit. Data were extracted with Feature Extraction software 10.7 (Agilent Technologies). Raw data were normalized with Quantile algorithm, LIMMA package in R. This step was completed with the help of Shanghai Biotechnology Co., Ltd.
Microarray data analysis
Gene ontology enrichment analysis of the target messenger RNAs of differentially expressed circRNAs
Gene ontology (GO) is a major bioinformatics initiative to integrate the expression of gene and gene product attributes. GO covers three domains: cellular component, molecular function, and biological process. Analyses were performed using DAVID (Database for Annotation, Visualization and Integrated Discovery; http://david.abcc.ncifcrf.gov/) software18–20 and Fisher’s exact test. We selected the standard as the difference in the number of genes on a term/GO ≥ 2, p < 0.05. The term/GO was drawn based on the enrichment factor value from the descending order of size, taking into account the first 30 results.
KEGG pathway enrichment analysis of target messenger RNAs of differentially expressed circRNAs
KEGG (Kyoto Encyclopedia of Genes and Genomes) is a database for systematic analysis of gene function, linking genome, and higher order information. This resource facilitates the study of gene and expression data as a whole network. Our study was performed using the KEGG database (http://www.genome.jp/kegg/)21,22 and Fisher’s exact test. We selected the standard as the difference in the number of genes on a term/pathway ≥ 2, p < 0.05. The term/pathway was drawn based on the enrichment factor value from the descending order of size, taking into account the first 30 results.
Prediction of circRNA–miRNA interactions
Prediction of interactions between circRNAs and miRNAs was based on TargetScan 6.0, which identifies miRNA targets and subsequently determines whether a given target is conserved or not across a given set of species.
CircRNA–miRNA co-expression network
To further elucidate the correlations between differentially expressed circRNAs and miRNAs in gastric cancer, we performed miRNA–circRNA interaction network analysis. The size of each node represents the number of putative miRNAs functionally linked to individual circRNAs. We identified all the different circRNAs and selected the top 30 circRNAs with differential expression as drawing objects.
Differentially expressed gene networks
In total, 306 differentially expressed genes were constructed using the STRING 23 (Search Tool for the Retrieval of Interacting Genes) database. A combined score was computed, which indicates a higher confidence when more than one type of information supports a given association.
Real-time polymerase chain reaction
Total RNA was extracted from cells with TRIzol reagent (Invitrogen, Carlsbad, CA, USA) and evaluated by measuring absorbance at 260 nm. Complementary DNA (cDNA) was synthesized with the ImProm-II Reverse Transcription System (Promega, Madison, WI, USA). Quantitative real-time polymerase chain reaction (qRT-PCR) was performed with strand-specific cDNA (0.5 µL), forward and reverse primers (0.5 µL), and SYBR Green Supermix (12.5 µL). qRT-PCR was performed using the 7000 Sequence Detection System (ABI, Foster, CA, USA). Relative expression levels were calculated using the 2−ΔΔCT (ΔCt = Ctgene−Ctcontrol) method.
Statistical analysis
SPSS software (version 17.0 for Mac OS X; SPSS, Chicago, IL, USA) was employed for statistical analysis. Values were calculated as mean ± standard deviation (SD) of three separate experiments. The Wilcoxon signed-rank test was used to analyze significant differences between circRNA expression levels. Differences or correlations between groups were calculated with Student’s t test. p < 0.05 was considered to be statistically significant.
Results
circRNA profile overview
High-throughput microarray was employed to detect circRNA expression in advanced gastric adenocarcinoma and adjacent normal mucosa tissues using the Circbase database including a total of 88,371 circRNAs. Original data were normalized and the Box plot (Figure 1(a)) used to display the overall characteristics of sample data distribution. Raw data of the chip were normalized and transformed into log2 values. In a two-dimensional coordinate system, a scatter plot (Figure 1(b)) was created and overall distribution of the two sets of data evaluated. We compared chip data between groups by screening the differences in fold change and performed the t test to determine the p value representing significant differences in the probe. The volcano plot (Figure 1(c)) shows the number and distribution of probes in the same plane that satisfy both the above conditions. Samples and genes were separately classified and gene expression in different samples was assessed using a heatmap (Figure 1(d)). Overall, expression levels of 165 circRNAs were decreased (fold change ≥ 2.0, p < 0.05) and those of 141 circRNAs increased (fold change ≥ 2.0, p < 0.05). Expression levels of hsa_circ_0023546, hsa_circ_0072040, and hsa_circ_0026674 in advanced gastric adenocarcinoma tissues were downregulated by 7.24-, 4.65-, and 4.56-fold, compared to adjacent normal mucosa tissues. Meanwhile, the levels of hsa_circ_0043637, hsa_circ_0043635, and hsa_circ_0043634 in advanced gastric adenocarcinoma tissues were upregulated by 64.41-, 63.24-, and 58.23-fold, compared to adjacent normal mucosa tissues (Table 2).
CircRNA microarray expression data between advanced gastric adenocarcinoma and adjacent normal mucosa tissues. (a) Normalized data were plotted as a Box plot to determine the overall characteristics of distribution. The distributions are almost the same after normalization. (b) Scatter plots of chip data assessing overall distribution of the two sets of data. Labeled red and blue portions indicate differentially expressed circRNA (two-fold or greater). (c) The volcano plot shows the number and distribution of probes that meet both conditions in the same plane. (d) The heatmap depicts an intuitive graphical representation of gene expression in different samples after hierarchical clustering of differentially expressed circRNAs.
CircRNA microarray expression data. The top 30 upregulated and 30 downregulated circRNAs, along with the degrees of differential expression, are listed.
CircRNA: circular RNA.
Real-time PCR validation of differentially expressed circRNAs
To verify microarray expression data, expression of a number of circRNAs (hsa_circ_0058246, hsa_circ_0091742, hsa_circ_0089548, hsa_circ_0003707, hsa_circ_0085553, hsa_circ_0006022, hsa_circ_0031979, and hsa_circ_0069131) was detected via real-time PCR in advanced gastric cancer and corresponding adjacent tissues from patients provided by Fudan University Shanghai Cancer Center (Figure 2). Both methods consistently and successfully revealed differential expression of circRNAs in tumors of advanced gastric cancer patients relative to the corresponding adjacent tissues, indicating that the biochip findings are highly accurate and reproducible.
Expression of circRNAs (hsa_circ_0058246, hsa_circ_0091742, hsa_circ_0089548, hsa_circ_0003707, hsa_circ_0085553, hsa_circ_0006022, hsa_circ_0031979, and hsa_circ_0069131) detected via real-time PCR using 20 pairs of tissues provided by the Fudan University Shanghai Cancer Center.
Bioinformatics analysis
Next, we predicted that the miRNAs sponge adsorbed by the circRNAs were differentially expressed in patients with advanced gastric adenocarcinoma and adjacent normal mucosa. Among the 306 circRNAs identified, 273 adsorbed related miRNAs and one circRNA could have one or more miRNA adsorption sites.
Using the corresponding miRNAs predicted using TargetScan, differentiated circRNAs were processed into interactions with miRNA (Supplementary Figure 1). We used the degree of interaction for quantification. To clearly observe the interactions between different circRNAs and miRNAs, we selected the top 30 upregulated and downregulated circRNAs and their associated miRNAs to map a network (Supplementary Figure 2). As evident from the Supplementary Figures 1 and 2, a single circRNA does not simply regulate multiple miRNAs but is also modulated by multiple circRNAs, constituting a complex regulatory network. We also use GO analysis and KEGG pathway annotation to determine the functions of differentially expressed circRNAs between advanced gastric adenocarcinoma and djacent normal mucosa tissues and related miRNAs (Figure 3(a)–(c)).
(a) GO enrichment analysis: The height and line graphs of the bar graph correspond to the number and percentage of genes, respectively. (b) GO enrichment analysis: Graph size represents the number of genes. The shape of the icon represents the three aspects of gene ontology, and the color represents the p value. (c) KEGG enrichment analysis: Graph size represents the number of genes, and graph color represents the p value.
To identify the node genes that play key roles in differentially expressed circRNAs, we used the STRING database, which integrates information from multiple databases and provides a combined score (Figure 4(a)). Consequently, the network relationships between differentially expressed circRNAs were mapped (Supplementary Figure 3). To determine the topological hub genes in the network, we selected genes with the highest degree. Overall, 29 genes were identified with degrees >20 (Figure 4(b)).
(a and b) To identify the topological hub genes in the network, genes with the highest degree were selected. In total, 29 genes with degrees >20 were identified. The height of the bar graph represents the level of the topological score.
We obtained follow-up data from patients who were enrolled in the study. Totally, we found 12 recurrence cases, including 3 cases of local recurrence, 4 cases of intraperitoneal extensive recurrence, 4 cases of liver metastases, and 1 case of distant metastasis (Table 1). We also compared the expression level of hsa_circ_0058246 in different clinicopathological stage group. We could not observe a statistically significant difference in the expression of hsa_circ_0058246 in the following groups: TNM stage group (IIIA vs IIIB vs IIIC) and T stage group (T2 vs T3 vs T4) (Figure 5(a) and (b)). Expression of hsa_circ_0058246 was significantly upregulated in group which have more than 16 lymph node invasions examined than the group less than 16 lymph nodes (Figure 5(c)). Interestingly, the patients who suffered recurrence of gastric cancer have a statistically significant increase in the expression of hsa_circ_0058246 (Figure 5(d)).
(a) The comparision of hsa_circ_0058246 expression among TNM stage groups (IIIA vs IIIB vs IIIC). (b) The comparision of hsa_circ_0058246 expression among T stage groups (T2 vs T3 vs T4). (c) The comparision of hsa_circ_0058246 expression between groups which have different lymph node invasions (more than 16 lymph node invasions vs less than 16 lymph node invasions; *p < 0.05). (d) The comparision of hsa_circ_0058246 expression between groups whose patients suffered recurrence or not (*p < 0.05).
Discussion
In the classical gene expression model, the information encoded by the genome is conveyed in the form of RNA molecules in individual cells. Each RNA molecule consists of a linear base sequence. The rapid development of gene sequencing technology has resulted in considerable accumulation of RNA sequence data, including a substantial fraction of spliced transcripts. 24 Traditional classical sequencing methods often ignore the “noises” of these genes. 25 At present, the number of circRNAs is more abundant than previously believed and their functions are additionally considered more important. 7 CircRNA is a type of special non-coding RNA with a closed-loop structure abundantly expressed in eukaryotic cells. These molecules are usually expressed in tissues or during developmental stages, and the majority of sequences are highly conserved.7,24 They are novel key molecules that play important roles in regulation. Owing to their closed-loop structure, circRNAs are not affected by RNA exonuclease and are more stable and less degradable. 26 CircRNAs have recently been developed as novel clinical diagnostic markers with obvious advantages. Compared with miRNAs, the mechanisms of circRNA formation and regulation are more complex and the regulatory network is three-dimensional involving combinations of DNA and RNA at multiple levels in gene transcription, translation, and other stages of target gene expression.
A recent hypothesis involving ceRNA (competing endogenous RNA) suggests a new mechanism of RNA interactions. 27 MiRNAs are known to cause gene silencing via binding to messenger RNA (mRNA), and ceRNAs regulate gene expression by competitively binding to miRNAs. ceRNAs interact with miRNAs through miRNA response elements (MREs), in turn, affecting gene silencing by miRNAs, suggesting the presence of an RNA in miRNA regulatory pathway with biological significance. 27 CircRNAs function as miRNA sponges and may compete with other endogenous RNAs. 28 Based on this hypothesis, we predicted that miRNAs with specific functions are differentially expressed. Our results showed that 273 of 306 differentially expressed circRNAs adsorbed one or more miRNAs and one circRNA could have one or more miRNA adsorption sites.
Here, we have presented a network map of circRNA–miRNA interactions. Specific miRNAs regulated by multiple circRNAs also play important roles in the regulation of tumor pathways. For example, in our study miR-4530 could be adsorbed by eight circRNAs. As one of the eight diagnostic index miRNAs, miR-4530 was reported to be differentially expressed in pancreatic or bile duct cancer. 29 Similarly, miR-3925-3p, regulated by six circRNAs, has been highlighted as one of the most prominent markers of esophageal adenocarcinoma. 30 Another study demonstrated downregulation of miR-663b expression in pancreatic cancer cells because of hypermethylation in its putative promoter region. Overexpression of miR-663b has been shown to inhibit cell proliferation, invasion, and migration, inducing apoptosis. Furthermore, miR-663b is regulated by the long non-coding RNA, HOX transcript antisense RNA (HOTAIR), and exerts a tumor suppressive effect by targeting insulin-like growth factor 2 (IGF2) in pancreatic cancer. 31 In our experiments, miR-663b was regulated by seven differentially expressed circRNAs. This complex regulatory network has a cascade effect, leading to altered expression of specific circRNAs, in turn, affecting numerous downstream targets and causing a butterfly effect. The data additionally provide a theoretical basis for optimizing individualized disease-targeted gene therapy.
In this study, GO analysis and KEGG pathway annotation were conducted to determine the functions of differentially expressed circRNAs between advanced gastric adenocarcinoma and adjacent normal mucosa tissues and related miRNAs. GO enrichment data revealed the target genes that are mainly involved in the regulation of key biological processes, cellular components, and molecular functions. In biological processes, such as cell migration, motility and locomotion were closely related to tumor cell migration. Cell structure and organization-related GO, such as anatomical structure development, cellular component organization, also constituted a high proportion of the target genes. In addition, cell metabolism–related protein and carbohydrate derivative binding were important functions.18–20 In KEGG pathway analysis, a number of important pathways closely associated with cancer were identified. Phosphoinositide 3-kinase (PI3K)-AKT is an important intracellular signal transduction pathway in the regulation of the cell cycle, shown to be directly related to the resting and proliferation stages, cancer progression, and longevity. 32 PI3K promotes phosphorylation and activates AKT, which localizes in the plasma membrane. 33 In many cancer types, the pathway is hyperactive, reducing apoptosis and promoting proliferation, and therefore presents an important target for cancer therapy. Another important pathway is Ras signaling, which regulates a number of key functions related to cell growth, differentiation, and apoptosis through actions on effector substrates, such as accelerating growth by shortening the cell cycle, prolonging lifespan by decreasing cell sensitivity to apoptotic signals, and inducing cell transformation. The Ras signaling pathway is reported to trigger abnormal cell proliferation leading to tumorigenesis, with at least 30% human cancers caused by Ras oncogene activation.34–37 Other pathways, such as epidermal growth factor receptor (EGFR) tyrosine kinase inhibitor resistance, extracellular matrix (ECM)-receptor interactions, and natural killer cell–mediated cytotoxicity, are additionally involved in regulatory processes in tumor cells.38–40
Differences in expression levels and roles of genes in tumor and normal tissues are two aspects of the same thing. In our experiments, 306 differentially expressed genes were introduced into the STRING database 23 and a network map of the different genes generated based on the combined score. We additionally focused on identifying the topological hub genes within the network. In total, 29 genes were identified with a degree greater than 20, most of which were closely related to the tumor. The differentially expressed genes were associated with interacting networks, indicative of a close relationship in terms of phenotypic regulation. Genes at the important nodes represented those with key regulatory roles. A few genes with high topological scores may be used as an example. In non–small cell lung cancer, cell adhesion and cell-to-cell signaling are mediated by PI3K/AKT, extracellular signal–regulated protein kinases 1 and 2 (ERK1/2), and nuclear factor kappa B (NF-κB) pathways. NOTCH1 has emerged as a driver connecting active signaling pathways, providing a potential therapeutic target for lung cancer transmission. 41 Recent studies have demonstrated that the human NOTCH1 gene is involved in metastasis and cancer stem cell maintenance. The level of NOTCH1 is significantly associated with TNM stage, metastasis, and triple-negative breast cancer, indicating that the protein is involved in metastasis and closely related to breast cancer stem cells. 42 FYN, a member of the Src family kinase (SFK), is overexpressed in various cancer types and associated with cell motility and proliferation. FYN enhances the expression of epithelial–mesenchymal transition (EMT)-related transcription factors, supporting the theory that EMT participates in FYN-induced breast carcinogenesis. Moreover, FYN is regulated by forkhead box protein O1 (FOXO1) and mediates fibroblast growth factor 2 (FGF2)-induced transcription of EMT-related genes through the PI3K/AKT and ERK/mitogen-activated protein kinase (MAPK) pathways. 43 LYN expression is reported to enhance cell proliferation, migration, and invasion. In lung adenocarcinoma cells, LYN is a useful prognostic marker. 44 Another earlier study demonstrated that LYN kinase is essential for chronic lymphocytic leukemia (CLL) progression. LYN plays a critical role in the formation of a microenvironment supporting leukemic growth. 45
In gastric cancer, our understanding of the expression and function of circRNA is still relatively weak. In one recent study, the expression of circRNAs (including 82 upregulated and 98 downregulated) in gastric cancer was significantly screened. Among them, circPVT1 had the most significant increase. The expression of circPVT1 in gastric cancer tissues and adjacent tissues was detected by qRT-PCR. The results showed that circPVT1 expression in gastric cancer tissues was significantly higher than that in adjacent tissues. Apart from promoting gastric cancer growth, functional study further confirms that circPVT1 is highly expressed in gastric cancer and is closely related to patient prognosis. 46 Another group of researchers found that hsa_circ_0000190 may be a novel non-invasive biomarker for the diagnosis of gastric cancer. They also found that the sensitivity and specificity of hsa_circ_0000190 were higher than carcinoembryonic antigen (CEA) and cancer antigen 19-9 (CA 19-9). 47 They also pointed out that hsa_circ_0000096 may be used as a potential novel biomarker for gastric cancer. Hsa_circ_0000096 could affect gastric cancer cell growth and migration through regulating cyclin D1, cell division protein kinase 6 (CDK6), matrix metalloproteinase 2 (MMP-2), and MMP-9. 48
In conclusion, the differences in circRNA expression between advanced gastric adenocarcinoma and adjacent normal mucosa tissues support an important role of these molecules in gastric cancer tumorigenesis. Data from the current study provide a new perspective on the mechanisms underlying gastric cancer metastasis and recurrence and highlight potential therapeutic targets that should contribute to improving diagnosis and treatment.
Footnotes
Declaration of conflicting interests
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Funding
The author(s) received no financial support for the research, authorship, and/or publication of this article.
References
Supplementary Material
Please find the following supplemental material available below.
For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.
For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.
