Lianmin Ye1, Wumin Jin2. 1. Department of Intensive Care, The First Affiliated Hospital of Wenzhou Medical University, Wenzhou, China. 2. Department of Reproductive Medicine Centre, The First Affiliated Hospital of Wenzhou Medical University, Wenzhou, China.
Abstract
BACKGROUND: Gastric cancer (GC) is one of the common digestive malignancies worldwide and causes a severe public health issue. So far, the underlying mechanisms of GC are largely unclear. Thus, we aim to identify the long non-coding RNA (lncRNA)-associated competing endogenous RNA (ceRNA) for GC. METHODS: TCGA database was downloaded and used for the identification of differentially expressed (DE) lncRNAs, miRNAs, and mRNAs, respectively. Then, the ceRNA network was constructed via multiple online datasets and approaches. In addition, various in vitro assays were carried out to validate the effect of certain hub lncRNAs. RESULTS: We constructed a ceRNA network, including 76 lncRNAs, 18 miRNAs, and 159 mRNAs, which involved multiple critical pathways. Next, univariate and multivariate analysis demonstrated 11 lncRNAs, including LINC02731, MIR99AHG, INHBA-AS1, CCDC144NL-AS1, VLDLR-AS1, LIFR-AS1, A2M-AS1, LINC01537, and LINC00702, and were associated with OS, and nine of those lncRNAs were considered as hub lncRNAs involved in the sub-ceRNA network. The in vitro assay indicated two lncRNAs, INHBA-AS1 and CCDC144NL-AS1, which were positively related to the GC aggressive features, including proliferation, invasion, and migration. CONCLUSIONS: We identified nine hub lncRNAs and the associated ceRNA network related to the prognosis of GC, and then validated two out of them as promising oncogenes in GC.
BACKGROUND: Gastric cancer (GC) is one of the common digestive malignancies worldwide and causes a severe public health issue. So far, the underlying mechanisms of GC are largely unclear. Thus, we aim to identify the long non-coding RNA (lncRNA)-associated competing endogenous RNA (ceRNA) for GC. METHODS: TCGA database was downloaded and used for the identification of differentially expressed (DE) lncRNAs, miRNAs, and mRNAs, respectively. Then, the ceRNA network was constructed via multiple online datasets and approaches. In addition, various in vitro assays were carried out to validate the effect of certain hub lncRNAs. RESULTS: We constructed a ceRNA network, including 76 lncRNAs, 18 miRNAs, and 159 mRNAs, which involved multiple critical pathways. Next, univariate and multivariate analysis demonstrated 11 lncRNAs, including LINC02731, MIR99AHG, INHBA-AS1, CCDC144NL-AS1, VLDLR-AS1, LIFR-AS1, A2M-AS1, LINC01537, and LINC00702, and were associated with OS, and nine of those lncRNAs were considered as hub lncRNAs involved in the sub-ceRNA network. The in vitro assay indicated two lncRNAs, INHBA-AS1 and CCDC144NL-AS1, which were positively related to the GC aggressive features, including proliferation, invasion, and migration. CONCLUSIONS: We identified nine hub lncRNAs and the associated ceRNA network related to the prognosis of GC, and then validated two out of them as promising oncogenes in GC.
Gastric cancer (GC) is one of the most frequent digestive system cancers and is the second cause of cancer mortality worldwide by 2018.
The cases in China account for more than 40% of the total number of GC worldwide due to a high incidence rate and a large population.
Moreover, the GC patients were more likely to be in the advanced stage when diagnosed because of non‐early specific symptoms. Unfortunately, the late diagnosis can significantly affect the 5‐year survival rate. In the past decade, due to the advancement of treatment and medicine, the GC prognosis has improved, but it is still not satisfied due to the relatively short disease‐free survival duration.
Thus, it is challenging and necessary to explore the underlying mechanisms of GC and identify novel biomarkers or treatment targets.On the contrary, the long non‐coding RNA (lncRNA) is a well‐known member of the non‐coding RNA family in the past decade, which is RNA with a length of over 200 nt.
,
In the past decade, the accumulating knowledge indicates the important role of the aberrant expression of lncRNA in GC.
,
LncRNA can function as sequence‐specific recruitment of proteins, competing for endogenous RNA (ceRNA) regulation, and molecular scaffolding of protein complexes, in which the ceRNA regulation is widely investigated nowadays. It is hypothesized that lncRNA can modulate the miRNA‐regulated mRNA expression by competitively binding miRNAs through endogenous molecular sponges.
This regulatory mechanism interprets the roles of lncRNA in various cancers, including GC.
,
,
,
However, the lncRNA‐associated ceRNA networks are far from clear.To further explore the role of specific lncRNA‐miRNA‐mRNA axis in GC, we first construct the ceRNA network via the online database. In addition, the hub lncRNAs with sub‐ceRNA networks related to prognosis were identified. To confirm the reliability and validity of the results, hub lncRNAs were validated in vitro. Overall, the present study aimed to establish a critical ceRNA network and identified novel diagnostic/therapeutic targets.
MATERIAL AND METHODS
Data resources and differential expression analysis
We downloaded the RNA sequence data with log2 (fpkm + 1) transformed (lncRNA and mRNA, level 3; Illumina HiSeq RNA‐Seq platform), miRNA sequence data (Illumina HiSeq miRNA‐Seq platform), and clinical information from the Xena dataset (https://xenabrowser.net/datapages/?dataset=TCGA‐STAD.htseq_fpkm‐uq.tsv&host=https%3A%2F%2Fgdc.xenahubs.net&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443; https://xenabrowser.net/datapages/?dataset=TCGA‐BRCA.htseq_fpkm‐uq.tsv&host=https%3A%2F%2Fgdc.xenahubs.net&removeHub=https%3A%2F%2Fxena.treehouse.gi.ucsc.edu%3A443). Annotation information for the RNA data was provided by the Ensemble database derived from “biomaRt” package.
The differentially expressed (DE) lncRNAs, miRNAs, and mRNAs between 32 normal samples and 372 cancer samples were identified by the “limma” package.
We set the |log2(Fold change)|>1 and p‐value <0.05 as the thresholds to identify DE genes. DE lncRNAs, miRNAs, and mRNAs were presented as a volcano plot. The study was approved by the Ethics Committee of The First Affiliated Hospital of Wenzhou Medical University. Informed patient consent was not required as the results showed are based upon the data generated by the TCGA database.
ceRNA network construction
To construct a ceRNA network, we extracted the interactions between mRNA and miRNA/lncRNA via multiple online datasets. For the interaction between miRNA and mRNA, we used the “miRNAtap” package (https://bioconductor.org/packages/release/bioc/html/miRNAtap.html), which consisted of commonly used five reliable and online datasets, including Pictar (https://pictar.mdc‐berlin.de/), DIANA,
Targetscan,
miranda,
and mirdb.
Only when the interaction was identified in at least three datasets, we considered it to be potential interactions. In contrast, the lncRNA‐miRNA interactions were predicted through the miRcode.For ceRNA construction, we submitted the above interactions and expression dataset to the “GDCRNATools” package (http://bioconductor.org/packages/devel/bioc/vignettes/GDCRNATools/inst/doc/GDCRNATools.html), which can further filter the interaction based on those two criteria, (1) expression of lncRNA and mRNA must be positively correlated, and (2) those common miRNAs should play similar roles in regulating the expression of lncRNA and mRNA. Those two factors were indicated via the Pearson correlation and regulation similarity.
Prognosis‐related lncRNA identification
To identify the potential prognosis‐related lncRNA, we grouped the samples into high‐ and low‐expression subgroups based on the mean expression. Then, univariate Cox regression analysis was used to find the prognosis‐related lncRNA. Then, the lncRNA achieving statistical significance (p < 0.05) in the univariate analysis was submitted into multivariate Cox analysis adjusting with age, gender, TNM stage, and histological grade. In addition, the KM plot with log‐rank test was carried out to further validate the prognosis‐related lncRNA.
Gene sets enrichment analysis
To reveal the function of the lncRNA‐associated ceRNA network, the DE mRNAs derived from ceRNA were subjected to gene sets enrichment analysis based on gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and REACTOME annotation with the clusterProfiler package.
GO is a structured standard biological model established by the GO consortium, including biological processes (BP), molecular functions (MF), and cellular components (CC).
KEGG is widely used as a reference for integrating large‐scale molecular datasets generated by sequencing and high‐throughput experimental technologies.
RACTOME is an open‐source, open access, manually curated, and peer‐reviewed pathway database (https://reactome.org/). Gene sets with a p‐value <0.05 were considered significant.
Cell line culture, RNA extraction, and real‐time PCR
GC cell line, MKN45, was purchased from the Cell Bank of the Chinese Academy of Sciences. Cells were incubated at 37℃ in a humidified atmosphere containing 5% CO2 and cultured in Dulbecco's Modified Eagle's Medium (DMEM) (ThermoFisher Scientific) combined with 10% fetal bovine serum (FBS).The primers were designed and purchased from GenePharma (Table 1). The total RNA was extracted using the ReliaPrep RNA Cell Miniprep System (Promega). It was first reversely transcribed into complementary DNA through the iScript Reverse Transcription Kit (Bio‐Rad), then, SYBR Green Supermix (Bio‐Rad) was employed to perform quantitative‐PCR (qPCR), which was run in Bio‐Rad CFX 96 PCR instrument (Bio‐Rad). The 2−ΔΔCT method was used to evaluate the relative gene expression. The GAPDH and U6 were used as the internal control for lncRNA/mRNA and miRNA, respectively.
TABLE 1
Forward and reverse primer sequences used for quantitative PCR
Forward and reverse primer sequences used for quantitative PCRAbbreviations: F, forward; R, reverse.
siRNA knockdown in cell line
siRNA designed and chemically synthesized by GenePharma. The sequences of siRNA against INHBA‐AS1 were: 5′‐GUCUCAUGACCACAGCUAAtt‐3′ (Sense) and 5′‐UUAGCUGUGGUCAUGAGACct‐3′ (Anti‐Sense), and siRNA against CCDC144NL‐AS1 was 5′‐UAGGUAGAUGGUGGAAUGAtt‐3′ (Sense) and 5′‐UCAUUCCACCAUCUACCUAtg‐3′ (Anti‐Sense). Before the transfection, MKN45 cells were cultured with antibiotic‐free DMEM for 24 h in advance. Transfections of the siRNAs were manipulated via the lipofectamine RNAiMAX (Invitrogen) and the same protocols were carried out following the manufacturer's instruction in all those two lncRNAs.
Proliferation assay
To evaluate the proliferation speed, 3 × 103 cells were seeded into a 96‐well plate per well (Greiner bio‐one), and then a time‐series assay every 2 days was carried out in triplicate. We used the MTS CellTiter 96 One Solution Cell Proliferation Assay (Promega) to measure proliferation. With 1, 3, 5, and 7 days, the absorbance at 490 nm was measured using the microtiter plate spectrophotometer (Benchmark Plus, Bio‐Rad) according to the manufacturer's protocol. Subsequently, proliferation was normalized based on the absorbance of the first day and calculated by the changes between the readings.
Transwell migration and invasion assay
For cell migration and invasion assay, 24‐well transwell inserts with a pore size of 0.8 mm (Corning) were used. After siRNA transfection, 1 × 105 cells were seeded in the upper chamber, and 650 µl of complete medium was added to the lower chamber as a chemo‐attractant. The 100‐µl Matrigel (Corning) with a concentration of 20 mg/ml was pre‐coated above the insert for invasion assay. The cells were allowed to migrate or invade toward the chamber for 12 and 16 h, respectively. The migrated and invaded cells below the membrane were fixed with 4% paraformaldehyde, stained with DAPI (Beyotime Biotechnology), and quantified from microscopic fields.
Statistical analysis
The student's t test was used to compare the continuous variables. The chi‐squared test was employed for the categorical variables. p values have two tails and only when it is less than 0.05 was considered significant. All the figures and statistical analysis were carried out via R software (version 3.6.1).
RESULTS
Identification of differently expressed lncRNA, miRNA, and mRNA
There were 1485 DE lncRNAs (1260 up‐regulation and 225 down‐regulation), 312 DE miRNAs (290 up‐regulation 12 down‐regulation), and 4260 DE mRNAs (2347 up‐regulation and 1913 down‐regulation) identified (Figure 1).
FIGURE 1
Volcano plots showing up and downregulated (A) lncRNA, (B) miRNA, and (C) mRNA. The red dots represent high expression lncRNA/miRNA/mRNA with LogFC ≥1 and FDR <0.05, and the blue dots represent low expression of lncRNA/miRNA/mRNA with LogFC ≤ −1 and FDR <0.05
Volcano plots showing up and downregulated (A) lncRNA, (B) miRNA, and (C) mRNA. The red dots represent high expression lncRNA/miRNA/mRNA with LogFC ≥1 and FDR <0.05, and the blue dots represent low expression of lncRNA/miRNA/mRNA with LogFC ≤ −1 and FDR <0.05
ceRNA network construction and function analysis
To construct the ceRNA, we first extracted the interaction between DE miRNAs and DE lncRNAs, and DE miRNAs and DE mRNAs as well via multiple online datasets. Finally, there were 6082 pairs of lncRNA‐miRNA interactions predicted, and 938 interactions between miRNA and mRNA were identified. Based on those interactions, we established 24847 potential lncRNA‐miRNA‐mRNA axes consisting of 892 lncRNAs, 18 miRNAs, and 278 mRNAs preliminarily.Then, the ceRNA network was conducted via the “GDCRNATools” package with “gdcCEAnalysis” function, which utilizes the Pearson correlation and regulation pattern to further determine the promising ceRNA. Finally, a ceRNA network (909 edges and 253 nodes), including 76 lncRNA, 18 miRNA, and 159 mRNA, was constructed with Pearson correlation coefficient ≥0.5, Pearson correlation p‐value >0.05, and regulation similarity ≠ 0 (Figure 2).
FIGURE 2
The ceRNA network including lncRNAs, miRNAs, and mRNAs. Red and green for all nodes represent up and downregulated directions between normal and cancer tissues, respectively. The node of V shape: miRNA; the node of triangle shape: lncRNA; the node of the cycle: protein coding gene; connecting line of red: lncRNA‐miRNA; connecting line of blue: miRNA‐mRNA
The ceRNA network including lncRNAs, miRNAs, and mRNAs. Red and green for all nodes represent up and downregulated directions between normal and cancer tissues, respectively. The node of V shape: miRNA; the node of triangle shape: lncRNA; the node of the cycle: protein coding gene; connecting line of red: lncRNA‐miRNA; connecting line of blue: miRNA‐mRNANext, we performed gene sets enrichment analysis to understand the potential biological effect of this ceRNA network (159 DE mRNA). We first divided the DE mRNA into two groups, including up‐ (n = 33) and downregulated (n = 126) genes. Then, we employed GO, KEGG, and REACTOME datasets and identified 20 significant GO terms, 1 KEGG, and 22 REACTOME for upregulated genes. In meantime, there were 220 significant GO terms and 5 REACTOME for downregulated genes. Then, we displayed the top 20 significant gene sets in Figure 3.
FIGURE 3
GO terms, REACTOME, and KEGG interpretation for functions of (A) up and (B) downregulated mRNAs derived from ceRNA network in GC. BP, biological pathway; CC, cellular component; MF, molecular function
GO terms, REACTOME, and KEGG interpretation for functions of (A) up and (B) downregulated mRNAs derived from ceRNA network in GC. BP, biological pathway; CC, cellular component; MF, molecular function
Survival‐associated lncRNA and mRNA identification
To identify the potential prognosis‐related lncRNAs, we utilized the univariate Cox analysis to filter 76 lncRNAs derived from the above ceRNA network. Then, the lncRNAs and mRNAs with p < 0.05 were further subjected to multivariable Cox analysis with the adjustment of age, gender, histological grade, and TNM stage. Then, there were 11 lncRNAs associated with the overall survival (OS) (Table 2). Intriguingly, they were all negatively related to OS. The mean expression of lncRNA was utilized as a cut‐off value to determine the high‐ and low‐expression groups. Then, the log‐rank test was applied to validate the relationship between the OS and the lncRNA (Figure 4).
TABLE 2
Univariate and multivariate Cox analysis for lncRNA and mRNA for overall survival of GC
RNA
Symbol
Univariable cox analysis
Multivariable cox analysis
HR
95% CI lower
95% CI upper
p
HR
95% CI lower
95% CI upper
p
lncRNA
LINC02731
1.527
1.109
2.103
0.036
1.574
1.139
2.175
0.007
lncRNA
MIR99AHG
1.447
1.052
1.992
0.008
1.433
1.036
1.984
0.005
lncRNA
INHBA‐AS1
1.419
1.034
1.949
0.024
1.420
1.033
1.953
0.007
lncRNA
LINC02613
1.480
1.078
2.032
0.023
1.550
1.126
2.133
0.003
lncRNA
CCDC144NL‐AS1
1.439
1.046
1.981
0.018
1.445
1.047
1.994
0.006
lncRNA
VLDLR‐AS1
1.379
1.004
1.896
0.030
1.502
1.090
2.069
0.004
lncRNA
LINC01497
1.596
1.159
2.198
0.045
1.509
1.086
2.097
0.005
lncRNA
LIFR‐AS1
1.466
1.066
2.016
0.005
1.446
1.046
1.998
0.006
lncRNA
A2M‐AS1
1.472
1.070
2.024
0.044
1.437
1.042
1.980
0.004
lncRNA
LINC01537
1.502
1.093
2.064
0.002
1.580
1.142
2.185
0.003
lncRNA
LINC00702
1.457
1.057
2.007
0.025
1.410
1.019
1.950
0.006
mRNA
NOVA1
1.789
1.294
2.473
0.010
1.856
1.340
2.569
0.006
mRNA
NPAS3
1.442
1.050
1.981
0.023
1.579
1.145
2.178
0.005
mRNA
CACNB2
1.412
1.023
1.949
0.030
1.377
0.996
1.902
0.008
mRNA
PDE7B
1.543
1.118
2.129
0.015
1.474
1.064
2.042
0.007
mRNA
CACNA2D3
1.467
1.051
2.047
0.025
1.368
0.980
1.911
0.006
mRNA
COL21A1
1.488
1.057
2.096
0.047
1.536
1.076
2.192
0.004
mRNA
SLC35F1
1.376
1.002
1.891
0.004
1.410
1.025
1.939
0.008
mRNA
PKIA
1.436
1.036
1.991
0.019
1.401
1.009
1.945
0.007
mRNA
KLF9
1.385
1.008
1.904
0.017
1.327
0.957
1.839
0.006
mRNA
PKNOX2
1.616
1.153
2.265
0.012
1.788
1.267
2.522
0.003
mRNA
COL5A2
1.387
1.009
1.906
0.021
1.490
1.072
2.071
0.014
mRNA
MATN3
1.651
1.197
2.277
0.000
1.695
1.222
2.351
0.009
mRNA
SNAP25
1.436
1.046
1.971
0.024
1.359
0.986
1.874
0.008
FIGURE 4
Kaplan‐Meier survival analysis for the correlation of DE lncRNAs with overall survival of the GC patients. Patients with expression ≥ mean expression were considered as high expression and otherwise as low expression
Univariate and multivariate Cox analysis for lncRNA and mRNA for overall survival of GCKaplan‐Meier survival analysis for the correlation of DE lncRNAs with overall survival of the GC patients. Patients with expression ≥ mean expression were considered as high expression and otherwise as low expression
Reconstruction and function analysis of hub lncRNA‐associated ceRNA network
We assumed those 11 lncRNAs play critical roles in the GC‐related ceRNA network. Thus, we extracted the corresponding 59 mRNAs that interacted with these 11 lncRNAs derived from the ceRNA network. Then, the univariate and multivariate Cox analyses were employed to filter the prognosis‐related mRNAs with the same procedures as above. There were 13 out of 59 mRNAs showing the independent relationship with OS, which was also validated via KM plot (Table 2 and Figure 5).
FIGURE 5
Kaplan‐Meier survival analysis for the correlation of DE mRNAs with overall survival of the GC patients. Patients with expression ≥ mean expression were considered as high expression and otherwise as low expression
Kaplan‐Meier survival analysis for the correlation of DE mRNAs with overall survival of the GC patients. Patients with expression ≥ mean expression were considered as high expression and otherwise as low expressionNext, a sub ceRNA network (24 edges and 27 nodes), including 9 lncRNAs, 5 miRNAs, and 13 mRNAs, was reconstructed (Figure 6A). A total of 9 lncRNA were considered as hub lncRNA, including LINC02731, MIR99AHG, INHBA‐AS1, CCDC144NL‐AS1, VLDLR‐AS1, LIFR‐AS1, A2M‐AS1, LINC01537, and LINC00702. To reveal the potential biological function of this sub‐network, the functional enrichment analysis was carried out, and it found 8 significant GO terms, 8 KEGG, and 30 REACTOME sets, such as MAPK signaling pathway (Figure 6B).
FIGURE 6
The nine‐hub lncRNAs‐associated sub ceRNA network. (A) The sub‐network including lncRNAs, miRNAs, and mRNAs. Blue circle, red triangle, and yellow V shape represent miRNA, mRNA, and lncRNA, respectively. (B) Functional gene sets derived from GO terms, REACTOME, and KEGG
The nine‐hub lncRNAs‐associated sub ceRNA network. (A) The sub‐network including lncRNAs, miRNAs, and mRNAs. Blue circle, red triangle, and yellow V shape represent miRNA, mRNA, and lncRNA, respectively. (B) Functional gene sets derived from GO terms, REACTOME, and KEGG
INHBA‐AS1 and CCDC144NL‐AS1 are potential oncogenes in GC
To validate the above result, two lncRNAs, INHBA‐AS1 and CCDC144NL‐AS1, were selected for further investigation. Subsequently, to validate the findings, siRNA‐mediated silencing of lncRNA was measured by RT‐qPCR, and the knockdown efficiencies of INHBA‐AS1 and CCDC144NL‐AS1 were significant in the MKN45 cell line (Figure 7A).
FIGURE 7
The expression of INHBA‐AS1/hsa‐miR‐98/COL5A2 and CCDC144NL‐AS1/hsa‐miR‐128‐1/MATN3 axis and in vitro functional assay. (A) The expression of INHBA‐AS1 and CCDC144NL‐AS1 with siRNA knockdown. (B) Proliferation assay with INHBA‐AS1 and CCDC144NL‐AS1 knockout, respectively. (C) Migration assay with INHBA‐AS1 and CCDC144NL‐AS1 knockout, respectively. (D) Invasion assay with INHBA‐AS1 and CCDC144NL‐AS1 knockout, respectively. (E) The expression of hsa‐miR‐98 and hsa‐miR‐128‐1 with INHBA‐AS1 and CCDC144NL‐AS1 knockout, respectively. (F) The expression of COL5A2 and MATN3 with INHBA‐AS1 and CCDC144NL‐AS1 knockout, respectively. *** Indicated p < 0.001; ** indicated p < 0.01; * indicated p < 0.05
The expression of INHBA‐AS1/hsa‐miR‐98/COL5A2 and CCDC144NL‐AS1/hsa‐miR‐128‐1/MATN3 axis and in vitro functional assay. (A) The expression of INHBA‐AS1 and CCDC144NL‐AS1 with siRNA knockdown. (B) Proliferation assay with INHBA‐AS1 and CCDC144NL‐AS1 knockout, respectively. (C) Migration assay with INHBA‐AS1 and CCDC144NL‐AS1 knockout, respectively. (D) Invasion assay with INHBA‐AS1 and CCDC144NL‐AS1 knockout, respectively. (E) The expression of hsa‐miR‐98 and hsa‐miR‐128‐1 with INHBA‐AS1 and CCDC144NL‐AS1 knockout, respectively. (F) The expression of COL5A2 and MATN3 with INHBA‐AS1 and CCDC144NL‐AS1 knockout, respectively. *** Indicated p < 0.001; ** indicated p < 0.01; * indicated p < 0.05Next, the result indicated that the knockdown of the two lncRNAs both suppressed cell proliferation as determined by MTS assays (Figure 7B). The results of migration and invasion assay indicated that MKN45 cell line with INHBA‐AS1 and CCDC144NL‐AS1 knockdown significantly less migrated and invaded than their counterpart (Figure 7C,D). Based on the ceRNA axis, we selected hsa‐miR‐98 and hsa‐miR‐128‐1 to verify the relationships between lncRNA and miRNA. Using RT‐qPCR, we observed that the knockdown of INHBA‐AS1 and CCDC144NL‐AS1 significantly increased the hsa‐miR‐98 and hsa‐miR‐128‐1 expression level, respectively (Figure 7E). Accordingly, the expression of COL5A2 and MATN3, which are corresponding lncRNA‐related mRNAs, showed a significant decrease compared to the controls (Figure 7F). These results indicated that INHBA‐AS1 and CCDC144NL‐AS1 might have an oncogenic function and act as ceRNA to sponge miRNAs in GC.
DISCUSSION
So far, the GC is one of the top‐ranking digestive cancers and has become a worldwide public concern. Thus, it is important to investigate the potential biomarkers and therapeutic targets. In our study, we identified an lncRNA‐associated ceRNA network involving GC tumorigenesis, which was based on the analysis of gene expression data obtained from the TCGA databases. Then, we identified nine hub lncRNAs accompanied with the sub ceRNA network related to OS. Among those nine lncRNAs, we validated the two of them, INHBA‐AS1 and CCDC144NL‐AS1, in vitro and found they were promising oncogene in GC.As mentioned above, lncRNA can influence the expression of mRNA via competitively binding to shared miRNA, which is defined as ceRNA and may play a critical role in the regulation of cancer development and progression, including GC.
For instance, LINC00152 regulated GACAT3 via miR‐103, and both are positively associated with poor clinicopathological characteristics in colorectal cancer.
For GC, lncRNA LINC01133 can inhibit GC progression by sponging hsa‐miR‐106a‐3p and then influence the APC expression.
In addition, lncRNA PTENP1 can regulate PTEN expression via binding to miR‐106b and miR‐93 in GC.
Our study identified a ceRNA network in GC involved in upregulation of MET activates PTK2 signaling, MET promotes cell motility and non‐integrin membrane‐ECM interactions. The MET activating the PTK2 signaling is related to MET receptor activating the focal adhesion kinase FAK1, which plays crucial role in focal adhesions (FAs). Specifically, FAs are large macromolecular complexes of integrins that mediate cell‐ECMs interactions and facilitate the metastatic process in cancer.
,
Previous studies identified that FAs is strongly associated with metastasis and lower survival rates.
,
,
,
Moreover, FAs can impact various tumor behaviors, such as migration, invasion, and proliferation.
Then, MET promotes cell motility, which may contribute to GC progression.
Non‐integrin membrane‐ECM interactions, such as dystroglycan and 37/67 laminin receptor, is found to be related to various epithelial cancers.Subsequently, we further identified nine‐hub lncRNAs, including LINC02731, MIR99AHG, INHBA‐AS1, CCDC144NL‐AS1, VLDLR‐AS1, LIFR‐AS1, A2M‐AS1, LINC01537, and LINC00702. Those hub lncRNAs‐associated ceRNA subnetwork is involved in actin filament binding and MAPK signaling pathway. A filament is a form of dense meshwork generated by lamellipodia, which facilitates cellular movement and plays anessential role in tumor cell metastasis.
MAPK signaling pathway is involved in various promoting‐cancer mechanisms, such as anti‐drug, inflammation, and immune evasion.
,
,
,
In terms of individual lncRNAs, most of them were related to the development and progression of various cancers in other studies. For instance, MIR100HG has been validated as an oncogene in the development of myeloid leukemia in vitro.
In addition, it was positively related to worse prognosis in GC via datasets other than TCGA.
LncRNA INHBA‐AS1 can promote multiple invasion features, including cell growth, migration, and invasion in oral squamous cell carcinoma, which targets on hsa‐miR‐143‐3p.
The INHBA‐AS1 in GC plasma was overexpressed compared to it in controls without further function assay.
Knockdown of lncRNA CCDC144NL‐AS1 attenuated migration and invasion in endometrial stromal cells.
The expression of VLDLR‐AS1 was independently related to the worse prognosis in thymoma.
The LIFR‐AS1/hsa‐miR‐29a/TNFAIP3 axis played an effect on the resistance of photodynamic therapy in colorectal cancer.
High expression of LIFR‐AS1 was correlated with poor survival in GC.
Upregulated A2M‐AS1 was associated with invasion and migration in breast cancer.
Besides, LINC00702 enhanced the progression of ovarian cancer through increased EZH2 expression.
Then, LINC00702/has‐miR‐4652‐3p/ZEBI axis can promote the progression of malignant meningioma through activating the Wnt/β‐catenin pathway.
Taken together, most of those nine‐hub lncRNAs were promising tumor‐promoting genes in diverse cancer and were worthwhile for further investigation in GC.Then, to validate our findings, INHBA‐AS1 and CCDC144NL‐AS1 and related axis were further verified in vitro and showed the promoting influence on proliferation, migration, and invasion. This indicated that two lncRNAs were promising oncogenic genes in GC. In terms of the related mRNA, INHBA‐AS1‐regulated COL5A2 and CCDC144NL‐AS1‐regulated MATN3 are related to GC prognosis.
,
MATN3 is a member of the Matrilin protein family, a noncollagenous extracellular matrix, which is associated with diverse cancers.
,
,
Specifically, it can induce the expression of MMP1, MMP3, MMP13, pro‐inflammatory cytokines, iNOS, and COX2, indicating MATN3 can regulate extracellular matrix degradation.
The COL5A2, collagen‐type V alpha 2 chain, encodes an alpha chain for one of the low abundances fibrillar collagens. It plays a critical role in the pathological process in multiple cancers including colorectal cancer, ovarian cancer, and bladder cancer.
,
Moreover, COL5A2 was strongly correlated with cell‐extracellular matrix organization, vascularization, and EMTs process function, and those functions were known to be involved in cancer invasion and metastasis.
Those findings may partially explain the functional effect of INHBA‐AS1 and CCDC144NL‐AS1 in vitro.So far, there are a few studies that are similar to ours. One study applied GEO dataset and paired GC/non‐tumorous tissues to identify DE lncRNAs and miRNAs, respectively. Then, TarBase and miRcode were used to establish a lncRNA‐miRNA‐mRNA network. Subsequently, one pair of ceRNA was validated in vitro without functional assays.
Another one used a small number of poorly differentiated GC and normal tissues to conducted numerous DE lncRNAs and selected one lncRNA, LINC02535, for further study.
Specifically, DE genes related to LINC02535 were filtered and used to conduct functional and protein—protein interaction analysis. LINC2535 alone was positively associated with cell proliferation, migration, invasion, and wound healing and negatively related to cell apoptosis via in vitro assays. Besides, one study utilized the TCGA GC dataset to identify DE lncRNAs, and selected LINC01234 for further validation.
Then, LINC01234 was proved to be positively associated with poor clinical characteristics in GC patients. Besides, they identified potential functions of LINC01234‐ and LINCO1234‐related network, including transcription factor (TF)‐lncRNA regulation, miRNA‐lncRNA relationship, as well as lncRNA‐RNA‐binding proteins interactions, via bioinformatics analysis; however, there was no functional assay. In terms of our study, we also used the TCGA GC data, which is the most well‐known and comprehensive dataset so far. Then, we utilized not only TarBase and miRcode but also other well‐known online tools to predict lncRNA‐miRNA and miRNA‐mRNA relationships and established the ceRNA network. Except for those, we further validated two pairs of ceRNAs in vitro with functional assays. Taken together, our study applied the latest and comprehensive dataset and well design methods to conduct the updated critical ceRNA network in GC. This may compensate for the shortage in this field.There are several limitations to our study. First, we only employed TCGA dataset, a frequently used online comprehensive cancer database. Second, although we combined a well‐designed bioinformatics study and in vitro validation, there was no in‐depth laboratory evidence, for example, a dual‐luciferase reporter assay, and mice model. Third, there is no clinical result in the present study. Taken together, a few vital experiments accompanied by the prospective stududies will be helpful to further validate our findings in the future.
CONCLUSIONS
So far, the role of the ceRNA network in GC is far from understood. In the present study, we established a promising lncRNA‐miRNA‐mRNA triple ceRNA network and identified a ceRNA subnetwork with nine‐hub lncRNAs involved in the prognosis of GC patients. Then, we validated two lncRNAs, INHBA‐AS1, and CCDC144NL‐AS1, accompanied their corresponding miRNAs and mRNAs, as potential oncogenic roles in GC. These findings need to be further confirmed in the future.
CONFLICT OF INTEREST
The authors declare no conflict of interest in preparing this article.
ETHICAL APPROVAL
The study was approved by the Ethics Committee of The First Affiliated Hospital of Wenzhou Medical University. Informed patient consent was not required as the results shown are based upon the data generated by the TCGA database.
Authors: W-Q Ma; J Chen; W Fang; X-Q Yang; A Zhu; D Zhang; H-L Zhong; B Yang; Z Luo Journal: Eur Rev Med Pharmacol Sci Date: 2020-02 Impact factor: 3.507
Authors: Jacques Ferlay; Isabelle Soerjomataram; Rajesh Dikshit; Sultan Eser; Colin Mathers; Marise Rebelo; Donald Maxwell Parkin; David Forman; Freddie Bray Journal: Int J Cancer Date: 2014-10-09 Impact factor: 7.396