Identification of a Key Competing Endogenous RNA Axis, “COL5A1/miR-137-3p/FSTL1”, Associated With Gastric Cancer Progression by Integrated Bioinformatics Analyses and Experimental Verification

doi:10.21203/rs.3.rs-824329/v1

Download PDF

Research

Identification of a Key Competing Endogenous RNA Axis, “COL5A1/miR-137-3p/FSTL1”, Associated With Gastric Cancer Progression by Integrated Bioinformatics Analyses and Experimental Verification

https://doi.org/10.21203/rs.3.rs-824329/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Background

MicroRNAs (miRNAs) and their target genes have been shown to play an important role in gastric cancer (GC), but this role has not been fully clarified. Therefore, our goal was to find the key miRNA-mRNA regulatory network in GC by combining a variety of bioinformatics and experimental analyses.

Methods

Differentially expressed miRNAs (DEMs) and genes (DEGs) were screened from TCGA and GEO, respectively. Survival-related differentially expressed miRNAs (SRDEMs) were determined by Cox proportional hazards regression and lasso regression analyses. Differentially expressed target genes (DETGs) of SRDEMs were predicted by TargetScan and miRDB and overlapped with DEGs. We constructed a protein–protein interaction (PPI) network of DETGs and conducted weighted gene coexpression network analysis (WGCNA) to screen the hub genes. Then, qRT–PCR and western blotting were performed to detect the expression level, and a dual-luciferase reporter assay was conducted to verify the binding of miRNA and mRNAs. CCK-8, EdU, wound healing and Transwell assays were conducted to compare the proliferation, migration and invasion abilities of GC cells in the different groups.

Results

We identified 11 SRDEMs and 233 DETGs, from which we selected miR-137-3p and its target gene COL5A1 for further research because of their key roles in the results of the bioinformatics analyses. Then, we showed that miR-137-3p was significantly downregulated in GC and that overexpression of miR-137-3p suppressed the proliferation, invasion and migration of GC cells by targeting COL5A1. Furthermore, we found that COL5A1 could regulate the expression of FSTL1 and GC progression by sponging miR-137-3p. Finally, bioinformatics analyses showed that FSTL1 might promote GC progression by regulating the immune infiltration of GC.

Conclusions

miR-137-3p played a tumor-suppressive role in GC, and its target gene COL5A1 could competitively bind miR-137-3p to upregulate the expression of FSTL1, which affects immune infiltration.

Clinical Pharmacology

Cancer Biology

Oncology

Bioinformatics analysis

miR-137-3p

COL5A1

FSTL1

ceRNA

gastric cancer

Gastric cancer (GC) is one of the most common malignant tumors in the digestive tract, and its incidence and mortality rank fifth and fourth, respectively, among all cancers [1]. Patients with GC have no specific clinical symptoms in the early stage, and most patients are in advanced stages when they seek medical treatment, which leads to poor therapeutic effects and prognoses [2]. At present, there are few clinically targeted drugs for GC, and the effect is not satisfactory. Therefore, it is necessary to find new potential targets for the treatment of GC through effective strategies.

MicroRNAs (miRNAs) are an important type of gene expression regulator. Mature miRNAs are endogenous, single-stranded, non-coding RNAs that can be selectively embedded in the RNA-induced silencing complex (RISC), inhibiting protein translation by incomplete binding to target mRNAs or causing degradation of target mRNAs by almost complete binding to target mRNAs to regulate genes at the posttranscriptional level [3–4]. Previous studies have shown that the expression of a variety of miRNAs in GC is dysregulated and affects the prognosis of patients in terms of tumorigenesis, metastasis, and tolerance to chemotherapy and radiotherapy [5–8], demonstrating the excellent potential of miRNAs to become targets for the treatment of GC [9–10]. However, this kind of potential has not been well exploited because of individual differences, as well as differences in sample sizes or detection platforms and methods for miRNA microarray analyses [11–12]. Therefore, how to systematically and effectively identify miRNAs that are dysregulated and play important roles in GC has become a new problem.

Bioinformatics is applied to examine the collection, processing, storage, dissemination, analysis and interpretation of biological information. It can help to reveal the biological mysteries endowed by a large amount of complex biological data through the comprehensive utilization of biology, computer science and information technology [13]. With the gradual development of molecular biology technology, bioinformatics has become increasingly important in cancer research. Commonly used databases, such as The Cancer Genome Atlas (TCGA), Gene Expression Omnibus (GEO), Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO), are often used to study various cancers [14–17]. Therefore, bioinformatics technology is expected to become a method to solve the above problem. To date, many studies have reported the identification of miRNAs and their target genes that are important in the oncogenesis and development of GC by various bioinformatics analyses [18–19]. However, many of the previous studies only conducted bioinformatics analyses without experimental verification, or the analytic method and data source were too narrow. These deficiencies limit the application of bioinformatics technologies in screening out key miRNAs and their target genes in GC. Whether the miRNAs and their target genes they screened out have universality and important functions is unclear.

Based on the above situation, the aim of this study was to analyze and identify GC-related miRNAs and their target genes by integrating bioinformatics resources and carrying out experimental verification to identify potential therapeutic targets and to provide a new direction for the treatment of GC.

Tumor data download and differential expression analysis

We downloaded miRNA stem loop expression data, including 446 GC samples and 45 normal samples, from the TCGA database on August 26, 2020, and performed differential analysis using the edge package in R software (version 3.6.1) to find differentially expressed miRNAs (DEMs) [20]. Additionally, we downloaded GSE118916 dataset, including 15 GC samples and 15 paired normal samples, from the GEO database and performed differential analysis using the limma package in R software to obtain differentially expressed genes (DEGs) [21]. The cutoffs were |log₂FC| > 1.0 and FDR (or P_adj) < 0.05, where FC denotes the fold change and FDR denotes the false discovery rate.

Cox proportional hazards regression model based on DEMs

We used the survival package in R software to perform univariate Cox proportional hazards regression analysis of DEMs to evaluate the impact of a single miRNA on the survival of GC patients and screen survival-related differentially expressed miRNAs (SRDEMs). The screening cutoff was p < 0.05. Then, we used SRDEMs in the multivariate Cox proportional hazards regression analysis and the glmnet package for lasso regression analysis to eliminate some relatively unimportant miRNAs in SRDEMs. The remaining SRDEMs were used to construct a model that could predict the prognosis of GC patients, and the performance of this model was evaluated [22]. The risk scores for each patient were calculated using the following formula: miRNA risk score = β_miRNA1 × exp (miRNA1) + β_miRNA2 × exp (miRNA2) +... + β_miRNAn × exp (miRNAn), where β is the regression coefficient derived from the multivariate Cox proportional hazards regression model, and exp () is the expression level of SRDEMs [23]. We divided the patients into a high-risk group and a low-risk group according to the median risk score and compared the Kaplan–Meier survival curves of both groups. We then calculated the 1-year, 3-year and 5-year survival rates of both groups and plotted receiver operating characteristic (ROC) curves to test whether the prediction ability of the model was reliable [24].

Prediction of target genes of SRDEMs and functional enrichment analysis of DETGs

We queried the mature miRNAs of each SRDEM through miRBase on August 28, 2020 [25], predicted their target genes through TargetScan and miRDB [26–27], and regarded the target genes of mature miRNAs belonging to the same SRDEM as the target genes of that SRDEM. Then, we used the venndiagram package in R software to intersect the target genes of SRDEMs predicted by TargetScan and miRDB and then intersected these overlapping target genes with DEGs to obtain the differentially expressed target genes (DETGs) of each SRDEM. Next, we used the clusterprofiler package in R software for KEGG and GO enrichment analysis of DETGs (cutoff of p < 0.05).

Construction and analysis of PPI networks with DETGs

The Search Tool for the Retrieval of Interacting Genes (STRING) can help to understand protein-protein interactions (PPIs) by integrating a large number of known and predicted data between proteins [28]. To study the interaction between DETGs and obtain potential key genes, we constructed their PPI network through the STRING database on August 29, 2020. Genes with significant interactions were screened out based on a confidence score ≥ 0.4, and the filtered results were imported into Cytoscape software (version 3.8.0) for network visualization [29]. We used the cytoHubba plug-in to calculate the Maximal Clique Centrality (MCC) value of each node in the PPI network and selected the genes with the top 20 MCC values as the hub genes [30].

Weighted gene co-expression network analysis

Weighted gene co-expression network analysis (WGCNA) was used to analyze the coding genes in GSE118916 [31]. After setting the soft threshold to 12, the co-expression network met the scale-free distribution. We used a dynamic tree-cutting algorithm (module size = 30) to merge genes with similar expression patterns into the same module. Then, we selected three modules containing the most DEGs and DETGs for further analysis. In these three modules, we screened the top 250 weight edges and input them into Cytoscape for network visualization.

Cell culture

The human gastric epithelial cell line GES-1 and GC cell lines AGS, BGC823, HGC27 and SGC7901 were purchased from KeyGEN BioTECH (China). All cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM, Gibco, USA) supplemented with 10% fetal bovine serum (FBS, Biological Industries, Israel), 100 KU/L penicillin and 100 mg/L streptomycin. Cells were maintained in an incubator at 37°C with 5% CO₂.

Plasmid construction

We constructed three kinds of COL5A1 knockdown plasmids based on the pLKO.1-TRC plasmid and packaged them as lentiviruses. In addition, based on the pCDH plasmid, we constructed plasmids overexpressing the COL5A1 3'UTR, FSTL1 3'UTR and FSTL1. COL5A1 wt/mut, FSTL1 wt/mut firefly luciferase plasmid and Renilla luciferase plasmid were purchased from Shanghai Genechem Co., Ltd (China).

Transfection

Cells were seeded into 6-well plates and cultured for 24 hours. When the density reached 50–60%, Lipofectamine 3000 reagent (Invitrogen, USA) was used to transfect the plasmids and miRNA mimics or inhibitor (Biomics, China) according to the instructions.

RNA extraction and quantitative real-time PCR (qRT–PCR)

Total RNA was extracted using an RNA extraction kit (Promega, China). CDNAs of COL5A1 and FSTL1 were synthesized using a reverse transcription kit (Yesen, China) and detected using a qRT–PCR kit (Yesen, China). Reverse transcription and qRT–PCR detection of miR-137-3p were performed with an All-in-One miRNA qRT–PCR Detection Kit (GeneCopoeia, USA). All qRT-PCRs were performed on a Roche Lightcycler 480 instrument. With GAPDH or U6 as internal parameters, the relative COL5A1, FSTL1 and miR-137-3p levels were calculated by the 2 ^−ΔΔCt formula. All primers are shown in Table S1.

Western blotting

Cell proteins were obtained by RIPA buffer and quantified using the BCA method. Then, the proteins were separated by sodium dodecyl sulfate-polyacrylamide gel electrophoresis (SDS–PAGE) and transferred to polyvinylidene fluoride membranes (Millipore, USA). After blocking the membrane with 5% skim milk for 1 hour at room temperature, the membrane was immunoblotted with primary antibodies overnight and secondary antibodies for 1 hour and visualized on an Odyssey CLx Infrared Imaging System (LI-COR Biosciences, NE, USA). COL5A1, FSTL1 and GAPDH antibodies were purchased from Bioworld (China).

Cell Counting Kit-8 (CCK-8) assay

Cells were seeded at 5000 cells per well and cultured in 96-well plates for 24 hours. Then, 10 µL of CCK-8 reagent (Yesen, China) per well was added and incubated at 37°C for 2 hours in the next 4 days. The relative proliferation was calculated based on the absorbance of each well.

Ethynyldeoxyuridine (EdU) proliferation assay

Cells were seeded at 2000 cells per well and cultured in 96-well plates for 24 hours. EdU reagent was added and incubated for 2 hours according to the instructions provide with the kFluor488Click-iT kit (KeyGen BioTECH, China). EdU-positive cells were observed and photographed under a fluorescence microscope.

Transwell migration/invasion assay

Transwell chambers (Corning, USA) coated with Matrigel (BD Biosciences, USA) were used for the invasion assay, and chambers without Matrigel were used for the migration assay. After seeding 2×10⁵ cells per chamber, the upper wells were supplemented with DMEM without serum, and the lower wells were supplemented with DMEM containing 10% FBS. The cells were incubated for 48 hours at 37°C with 5% CO₂. Then, the chambers were fixed with 4% paraformaldehyde for 30 minutes and stained with crystal violet dye for 20 minutes. Finally, the cells were observed and photographed under an inverted microscope after removing nonmigrating or noninvading cells.

Wound healing assay

Cells were seeded at 3×10⁵ cells per well in 6-well plates and cultured for 24 hours. When the cells had completely covered each well, a 10-µL sterile tip was used to scratch the cell surface to create a wound. Pictures were obtained under an inverted microscope at 0 and 48 hours to calculate the width of the wound.

Dual-luciferase reporter assay

Cells were seeded at 5×10⁴ cells per well in 24-well plates and cultured for 24 hours. The dual-luciferase reporter assay was carried out according to the instructions provided with the Dual-Luciferase Reporter Assay Kit (Vazyme, China). Cell lysis and luciferase detection were performed 48 hours after transfection. The firefly luciferase activity was normalized to the Renilla luciferase activity. Data are expressed as the percent of luciferase activity in control cells (100%).

ESTIMATE and CIBESORT

We downloaded the mRNA expression data of 375 GC samples and 32 normal samples from the TCGA database on August 26, 2020. The proportion of stromal and immune components in the tumor microenvironment (TME) of each sample was estimated with the estimate package in R software and displayed in three forms: stromal score, immune score, and ESTIMATE score, which were positively correlated with the proportion of stromal components, immune components and their sum, respectively. The higher the scores, the greater are the proportion of corresponding components in the TME [32]. In addition, we estimated the content of 22 types of immune cells in all GC samples using the CIBESORT package in R software. Only 358 tumor samples with p < 0.05 were selected for analysis [33].

Statistical analysis

GraphPad Prism 8.0 (San Diego, CA) and R software were used to analyze the data, and a t-test was employed to evaluate differences between the different groups. p < 0.05 was considered statistically significant: *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001, #p > 0.05.

Cox proportional hazards regression model of DEMs

The flow chart of this study is shown in Fig. 1a. By analyzing the TCGA data, we identified 242 DEMs with statistical significance, which consisted of 178 upregulated miRNAs and 64 downregulated miRNAs (Fig. 1b). To determine which miRNAs significantly affected the survival of GC patients, we identified 18 SRDEMs by univariate Cox proportional hazards regression analysis (p < 0.05, Table 1). Then, we used these SRDEMs for multivariate Cox proportional hazards regression analysis. Lasso regression analysis was utilized to eliminate some relatively unimportant miRNAs in SRDEMs, and 11 SRDEMs (hsa-miR-328, hsa-miR-549a, hsa-miR-708, hsa-miR-217, hsa-miR-371a, hsa-miR-7-2, hsa-miR-675, hsa-miR-137, hsa-miR-548v, hsa-miR-2115, hsa-miR-3923) were screened out to construct a prognostic model (Fig. 1c-d). The risk score was calculated as follows: miRNA risk score = (0.07279 × hsa-miR-328) + (-0.16586 × hsa-miR-549a) + (0.07277 × hsa-miR-708) + (0.06800 × hsa-miR-217) + (0.06526 × hsa-miR-371a) + (-0.14550 × hsa-miR-7-2) + (0.04510 × hsa-miR-675) + (0.05038 × hsa-miR-137) + (0.16243 × hsa-miR-548v) + (-0.02891 × hsa-miR-2115) + (0.20189 × hsa-miR-3923). The forest plot of multivariate Cox analysis is shown in Fig. 1e, and the nomogram of the model is shown in Fig. 1f. The calibration charts for predicting the 3- and 5-year survival rates of GC patients from TCGA based on the nomogram are shown in Fig. 1g and Fig. 1h, respectively. The 3- and 5-year survival rates predicted by the nomogram were basically consistent with the actual 3- and 5-year survival rates of GC patients. Kaplan–Meier survival analysis showed a significant difference in survival between the high-risk and low-risk groups (p < 0.001, Fig. 1i). The areas under the 1-year, 3-year and 5-year ROC curves were 0.654, 0.787 and 0.735, respectively, indicating that the model could effectively predict the prognosis of GC patients (Fig. 1j).

Table 1

Univariate Cox regression analysis of the 18 miRNAs associated with survival in GC patients.
miRNA	HR	p-value
hsa-miR-7-2	0.837016	0.000874
hsa-miR-328	1.263305	0.002276
hsa-miR-3923	1.244871	0.002407
hsa-miR-675	1.095968	0.002585
hsa-miR-7-3	0.864100	0.006473
hsa-miR-549a	0.840086	0.008898
hsa-miR-125a	1.271646	0.009330
hsa-miR-708	1.158467	0.017433
hsa-miR-217	1.124465	0.017771
hsa-miR-137	1.137458	0.018146
hsa-miR-100	1.120973	0.018551
hsa-miR-2115	0.822270	0.028516
hsa-miR-548v	1.136731	0.029136
hsa-miR-6511b-1	1.201655	0.034219
hsa-miR-187	1.067162	0.034778
hsa-miR-145	1.096116	0.042586
hsa-miR-371a	1.113365	0.044637
hsa-miR-216a	1.114984	0.048947

Prediction of target genes of SRDEMs and functional enrichment analysis of DETGs

Since miRNAs generally function by targeting genes, we used two independent online tools (TargetScan and miRDB) to predict the target genes of 11 SRDEMs. As shown in Figure S1A-K, after intersection of the prediction results from two websites, 114, 408, 496, 240, 274, 695, 279, 851, 497, 1303 and 196 target genes of hsa-miR-328, hsa-miR-549a, hsa-miR-708, hsa-miR-217, hsa-miR-371a, hsa-miR-7-2, hsa-miR-675, hsa-miR-137, hsa-miR-548v, hsa-miR-2115 and hsa-miR-3923 were determined, respectively. In addition, analysis of the GSE118916 dataset identified 1080 DEGs with statistical significance, including 474 upregulated genes and 606 downregulated genes (Fig. 2a). By intersecting 1080 DEGs with the predicted target genes, 6, 20, 38, 13, 13, 31, 23, 56, 26, 82 and 13 DETGs of these SRDEMs were obtained, with a total of 321 target relationships and 233 DETGs (some DETGs were targeted by more than one DETG, Table 2, Figure S2). To further reveal the potential biological functions of these 233 DETGs, we performed KEGG and GO enrichment analyses. The results of the KEGG analysis showed that these DETGs were mainly enriched in the PI3K-Akt signaling pathway, focal adhesion, proteoglycan in cancer, extracellular matrix and other cancer-related pathways (Figure S3A). The results of the GO analysis showed that these DETGs were mainly involved in biological processes such as cell proliferation, cell migration, extracellular matrix and angiogenesis (Figure S3B). These results suggested that these DETGs might play an important role in the development of GC.

Table 2

Two hundred and thirty-three DETGs shared by the target genes of 11 SRDEMs and DEGs
SRDEM	DETGs
hsa-miR−328	ESRP1, IFIT1, MICALL1, MLLT11, MN1, PRDM16
hsa-miR−549a	ATP11A, C1R, CDH11, EPPK1, FUT9, GLUL, GPD1L, LAPTM5, NIPAL1, OCLN, PPP2R3A, PRELID2, RAI14, SGK1, SH3BGRL2, SOX21, ST3GAL6, SYNJ2BP, TFCP2L1, THBS2
hsa-miR−708	AADAC, ABHD2, AHCYL2, ANK3, AQP4, BCL6, CAP2, CCL28, CCND2, CDH2, CMTM4, COL10A1, EMP1, EPB41L3, ETNK1, FGD4, FUT2, GABRP, GPR155, IFIT1, KIAA1958, MAGI3, MAP1B, NLRP1, PAFAH2, PLSCR1, PTGER3, REPS2, RGN, RHOBTB1, SCRG1, SIK2, STS, TMEM161B, UGT8, XK, ZEB1, ZNF462
hsa-miR-217	ABCC9, ANK3, ANLN, ATP1B1, CALD1, CORO1C, ESRRG, FN1, HOMER2, SEMA3A, SLC4A4, TMEM246, UBL3
hsa-miR-371a	ABHD2, CALU, ELOVL6, HAND2, MSRB3, NEUROD1, OCLN, SGPP2, SH3BGRL2, SLC16A7, SYNJ2BP, TMEM185B, TOB1
hsa-miR-7-2	BUB1, CALU, COL2A1, CTSB, DSP, EHF, EPB41L3, ESRRG, GALNT3, GATA6, GATM, GJA1, GJC1, GREM1, HOMER2, IL20RA, KLF4, MAP1B, NTN4, PDE4B, PDZRN4, PRKCB, PTGER3, SGK1, SH3BGRL2, SKAP2, SLC16A7, SLC4A4, SPP1, SYNJ2BP, TCF4
hsa-miR-675	ABCC9, ANKRD22, ATP11A, C6, CDH11, CDH13, COBLL1, COL5A2, COLCA1, CORO1C, DMRTA1, GC, ISPD, MYBL1, PIGR, PIP5K1B, SERPINF1, SLC12A2, SMOC2, STS, TFEC, UGT8, USP53
hsa-miR-137	ABCC9, ABHD2, AHCYL2, AJUBA, ANGPTL2, ATP1B1, CNTN3, COL5A1, CXADR, DSP, DUSP10, DUSP4, EHBP1L1, ELOVL6, ESRRG, FAT3, FBXO32, FGD4, FGL2, FRMD6, FSTL1, FUT9, GFPT1, GJC1, GPX7, GULP1, KIAA1958, KLF15, KLF4, LBH, LCP2, MPC1, MSRB3, NEUROD1, NIPAL1, NT5DC2, PDLIM3, PKDCC, PLEKHO2, RCN3, RFTN1, RHOBTB1, RRAGD, SERPINA3, SIK2, SIPA1L2, SLC12A2, TBC1D1, TCF4, THBS4, TMEM125, TMEM56, TRPS1, TWIST1, YBX3, ZNF385B
hsa-miR-548v	ABHD2, APOBEC1, ATP1B1, CCND2, CDH2, CHN1, CLDN1, CMTM4, FUT9, GALNT3, GATA4, GPSM2, HOXB7, MAP4K4, MRAP2, NAT1, NNT, PRKACB, PRSS8, PTPRZ1, RAB27B, SLC16A7, SLC4A4, SYNC, TMEM92, TNFSF13B
hsa-miR-2115	ABCC5, ADAM28, ALCAM, ALDH6A1, ARSD, ATP11A, BCAS1, BEX5, C6orf58, CA8, CCL28, CD44, CDS1, CENPF, CEP170, CKS1B, CNTN3, COL11A1, CORO1C, CRISPLD1, CSGALNACT2, CXADR, CXCL9, DAB2, DSP, DUSP10, EHF, EMILIN2, EMP1, FAM83F, FAT3, FGF7, FOXA1, FSTL1, FUT8, FUT9, GALNT6, GNB4, GPD1L, GPR155, GPSM2, GSTA1, GSTA2, HOOK1, HOXB6, IGF1, KL, LRRC17, NEURL1B, NID2, NQO1, NSUN7, PDE4B, PEA15, PLAU, PMP22, PRELID2, PRKCB, PTPN3, PTPRN2, RAB11FIP2, RAB27B, RAP1GAP, RHOBTB1, SEMA3A, SGK1, SH3GL2, SIK2, SIPA1L2, SLC16A7, SMOC2, SPARC, SULF1, SYNJ2BP, TFCP2L1, TMEM161B, TMEM229A, TRPS1, TSPAN12, TWIST1, UGT8, ZEB1
hsa-miR-3923	ATP1B1, DSC2, GKN1, GRIA4, KL, LIFR, MSR1, NEUROD1, ODAM, PIP5K1B, PTGS2, SAMD13, SOSTDC1
Bold represents genes targeted by more than 1 SRDEM.

Identification of hub genes in DETGs

A SRDEM can target multiple DETGs, and a DETG can also be targeted by multiple SRDEMs, indicating that there are correlations between these DETGs. To clarify their interactions, we input all DETGs into the STRING database to construct the PPI network (Fig. 2b), which contained 148 nodes and 355 edges (several discrete nodes were eliminated). Each node in the network represents a gene, and the connection between nodes symbolizes the interactions between proteins encoded by corresponding genes. The MCC value can be used to evaluate the importance of each node and help to find the hub genes in the network [30]. Therefore, we used the cytoHubba plug-in to calculate the MCC value of each node. The nodes with the top 30 MCC values are shown in Fig. 2c. Then, we regarded the genes with the top 20 MCC values as the hub genes in the network: FN1, CD44, IGF1, SPARC, CDH2, SPP1, GJA1, COL5A2, TWIST1, COL5A1, COL2A1, ZEB1, FSTL1, PTGS2, KLF4, COL11A1, OCLN, PLAU, THBS2, and GATA4. The PPI network of these 20 hub genes is shown in Fig. 2d.

WGCNA identification of the hub module

WGCNA can be used to analyze the gene expression patterns of multiple samples and can cluster genes with similar expression patterns [31]. Therefore, we performed WGCNA to analyze the GSE118916 dataset to determine a class of genes playing an important role in GC. With different power values ranging from 1 to 20, we determined that the best soft threshold was 12 (Fig. 2e), and then 32 different modules were generated according to the independence and average connectivity of the network. The co-expression degree of genes in the same module was high, while the co-expression degree of genes from different modules was low (Fig. 2f). Among them, the top 5 modules containing the most DEGs and DETGs are shown in Table S2. Moreover, we constructed weighted gene co-expression networks for the top 3 modules (bisque4, blue, dark gray) using the top 250 weight edges in each module and input them into Cytoscape for network visualization (Fig. 2g, Figure S4A-B). We found that although the bisque4 module contained the most DEGs and DETGs, it contained only 4 hub genes, while the dark gray module (Fig. 2g) had fewer DEGs and DETGs but contained 8 hub genes. Thus, we identified this module as the hub module. Interestingly, we found that collagen genes (COL genes) accounted for a very high proportion of the hub genes in the PPI network (COL5A2, COL5A1, COL2A1, COL11A1), while COL genes also frequently appeared in the dark gray module, in which only COL5A1 was located in the center, suggesting an important role of COL5A1. Bioinformatics analyses showed that COL5A1 was highly expressed in GC, while previous prediction analyses showed that the mature miRNA regulating COL5A1 was only miR-137-3p, which was expressed at low levels in GC, suggesting that miR-137-3p and COL5A1 were involved in the progression of GC. Therefore, we selected miR-137-3p and COL5A1 for further experimental verification.

MiR-137-3p suppresses the proliferation, invasion and migration of GC cells

First, we verified the expression of miR-137-3p and COL5A1 in GC cells. By comparing the expression levels of miR-137-3p (Fig. 3a) and COL5A1 (Fig. 3b-c) in GES-1 and 4 GC cell lines (AGS, BGC823, HGC27, SGC7901), we obtained results consistent with the bioinformatics analyses; that is, miR-137-3p was expressed at low levels in GC, while COL5A1 was highly expressed. Then, we transfected miR-137-3p inhibitor into AGS with high expression of miR-137-3p and miR-137-3p mimics into SGC7901 with low expression of miR-137-3p (Figure S5A). Using the CCK-8 (Fig. 3d) and EdU assays (Fig. 3e), we confirmed that the cell proliferation of AGS was enhanced after transfection of the miR-137-3p inhibitor, while the cell proliferation of SGC7901 was decreased after transfection of miR-137-3p mimics. We also found that the migration and invasion ability of AGS cells were enhanced after transfection with the miR-137-3p inhibitor, while the opposite results were observed in SGC7901 cells transfected with miR-137-3p mimics (Fig. 3f-g), suggesting that miR-137-3p played a tumor-suppressive role in GC.

To further clarify whether miR-137-3p could target and regulate the expression of COL5A1, we detected the effect of changing the expression level of miR-137-3p on the expression level of COL5A1. Compared with the control group, the expression of COL5A1 increased after transfection of the miR-137-3p inhibitor, while it decreased significantly after transfection of the miR-137-3p mimics, indicating that COL5A1 was regulated by miR-137-3p (Fig. 3h-i). Then, to verify whether miR-137-3p could directly bind to the COL5A1 3'UTR, we carried out a dual-luciferase reporter assay. First, we designed COL5A1-wt/mut firefly luciferase plasmid (Fig. 3j) for the binding site (as shown in Table S3, we selected the one with the highest score) between miR-137-3p and the COL5A1 3'UTR, cotransfected them with the NC or miR-137-3p mimics into 239 T cells, and detected the luciferase activity after 48 hours. Compared with the control group, the luciferase activity in the COL5A1-wt and miR-137-3p mimic cotransfected groups decreased significantly, while the luciferase activity in the COL5A1-mut and miR-137-3p mimic cotransfected groups did not obviously change (Fig. 3k). These results indicated that COL5A1 was the target gene of miR-137-3p.

COL5A1 can reversely sponge miR-137-3p and upregulate the expression of FSTL1 through a ceRNA mechanism

Studies have shown that when there are multiple miRNA binding sites (also known as miRNA response elements, MREs), the mRNAs of target genes have the potential to become miRNA sponges [34–35]. When searching for binding sites between COL5A1 mRNA and miR-137-3p, we found that the 3'UTR of COL5A1 mRNA had up to three high-quality MREs that could bind to miR-137-3p (Table S3), suggesting that COL5A1 might not only be the target gene of miR-137-3p but also the sponge of miR-137-3p, thus regulating the expression of some key genes.

Therefore, we further queried the genes that might be regulated by COL5A1 mRNA via the competing endogenous RNA (ceRNA) mechanism in the starBase database (3717 genes in total, Table S4) [36]. These genes were then intersected with the 20 hub genes mentioned above. Seven genes (FSTL1, CD44, FN1, KLF4, OCLN, PLAU, THBS2, Fig. 4a) may have been regulated by COL5A1 mRNA, of which only FSTL1 and KLF4 were the target genes of miR-137-3p. Then, we used the starBase database to query the expression correlations between COL5A1 and FSTL1 (Fig. 4b) or KLF4 (Figure S5B) in GC. We found that only the expression of FSTL1 and COL5A1 was significantly positively correlated (r = 0.805, p < 0.0001). Thus, we speculated that COL5A1 mRNA might regulate the expression of FSTL1 by competitive binding to miR-137-3p through a ceRNA mechanism.

To verify this hypothesis, we detected the expression level of FSTL1 in GES-1 and 4 GC cell lines. As expected, we found that the expression of FSTL1 was highly consistent with COL5A1; that is, it was also high in BGC823/SGC7901 cells with high COL5A1 expression and low in AGS cells with low COL5A1 expression (Fig. 4c-d). Moreover, the expression of FSTL1 increased after transfection of miR-137-3p inhibitor in AGS and decreased after transfection of miR-137-3p mimics in SGC7901, indicating that the expression of FSTL1 was also regulated by miR-137-3p (Fig. 4e-f).

Next, to verify whether miR-137-3p could directly bind to the FSTL1 3'UTR, we again carried out a dual-luciferase reporter assay. We designed FSTL1-wt/mut firefly luciferase plasmid for the binding site of miR-137-3p and the FSTL1 3'UTR (Fig. 4g, Table S3) and cotransfected them with the NC or miR-137-3p mimics into 239 T cells. The luciferase activity was detected 48 hours later. Compared with the control group, the luciferase activity in the FSTL1-wt and miR-137-3p mimic cotransfected group was significantly decreased, while that in the FSTL1-mut and miR-137-3p mimic cotransfected group was not significantly changed (Fig. 4h). These results suggested that FSTL1 was the target gene of miR-137-3p.

To further confirm that the COL5A1 3'UTR could regulate the expression of FSTL1 by competitive binding to miR-137-3p through a ceRNA mechanism, we designed plasmids overexpressing the COL5A1 3'UTR (p-COL5A1 3'UTR) and plasmids overexpressing the FSTL1 3'UTR (p-FSTL1 3'UTR) and cotransfected them into 293T cells with FSTL1-wt or COL5A1-wt luciferase plasmid and miR-137-3p mimics, respectively. As shown in Fig. 4i, compared with the control group (pCDH/FSTL1-wt/miR-137-3p mimics cotransfected group), the luciferase activity was significantly increased in the p-COL5A1 3'UTR/FSTL1-wt/miR-137-3p mimics cotransfected group, indicating that the COL5A1 3'UTR might sponge miR-137-3p and reduce its binding to FSTL1-wt, resulting in increased luciferase activity. However, the luciferase activity in the p-FSTL1 3'UTR/COL5A1-wt/miR-137-3p mimics cotransfected group showed no significant change compared with the control group (pCDH/COL5A1-wt/miR-13p mimics cotransfected group), indicating that the regulation of COL5A1 on FSTL1 was unidirectional; that is, FSTL1 3'UTR could not regulate the expression of COL5A1 in turn. The qRT–PCR and western blotting results confirmed the above hypothesis; that is, after overexpression of the COL5A1 3'UTR in AGS, the mRNA and protein levels of FSTL1 were significantly upregulated compared with those of the control group (Fig. 4j-k), while after overexpression of the FSTL1 3'UTR, the mRNA and protein levels of COL5A1 were not significantly changed (Fig. 4l-m). We also detected the expression level of miR-137-3p in the above groups. As shown in Fig. 4n, after overexpression of the COL5A1 3'UTR in AGS, the expression level of miR-137-3p was significantly downregulated, while overexpression of the FSTL1 3'UTR had no significant effect on miR-137-3p. Then, we designed 3 shRNA sequences for COL5A1 and confirmed that the knockdown effect of the COL5A1-sh1 sequence was optimal in BGC823 and SGC7901 cells with high expression of COL5A1 (Figure S5C-D). We also found that after knockdown of COL5A1, the expression of miR-137-3p was significantly upregulated (Fig. 4o). These results suggested that FSTL1 mRNA could not regulate COL5A1 by competitive binding to miR-137-3p. In contrast, COL5A1 mRNA could regulate FSTL1 expression by competitive binding to miR-137-3p through a ceRNA mechanism.

Knockdown of COL5A1 inhibits the proliferation, migration and invasion of GC cells, which can be rescued by miR-137-3p inhibitor or overexpression of FSTL1

To further study the effect of COL5A1 regulation of FSTL1 expression on the function of GC cells, we used COL5A1-knockdown cell lines for follow-up experiments. In addition, we designed a plasmid overexpressing FSTL1. Verification of the overexpression efficiency is shown in Figure S5E-F. The mRNA and protein levels of FSTL1 were significantly increased after transfection of AGS with the plasmid overexpressing FSTL1. Next, we designed four different experimental groups in BGC823 and SGC7901 cells that highly expressed COL5A1 and FSTL1: the control group, COL5A1-knockdown group, and transfection with miR-137-3p inhibitor or plasmid overexpressing FSTL1 in the COL5A1-knockdown group. First, we detected the cell proliferation of these four groups by CCK-8 and EdU experiments. As shown in Fig. 5a-b, compared with the control group, cell proliferation of SGC7901 decreased after knockdown of COL5A1, while it was rescued after transfection of miR-137-3p inhibitor or plasmid overexpressing FSTL1. Then, we performed healing and Transwell assays to detect the cell migration and invasion abilities in the 4 groups. We found that compared with the control group, the migration and invasion ability of SGC7901 decreased significantly after knockdown of COL5A1, and this effect could also be rescued by transfection of miR-137-3p inhibitor or plasmid overexpressing FSTL1 (Fig. 5c-d). Simultaneously, we also detected the expression level of FSTL1 in the four groups. The results showed that the expression level of FSTL1 decreased significantly after knockdown of COL5A1, while it was rescued or even higher than that in the control group after transfection of miR-137-3p inhibitor or overexpression of FSTL1 (Fig. 5e-f). We repeated the above experiments in BGC823 cells and obtained similar results (Figure S6). These findings indicated that COL5A1 could regulate the expression of FSTL1 by competitive binding to miR-137-3p to promote the proliferation, migration and invasion of GC.

FSTL1 is related to immune infiltration in the TME of GC patients

We were interested in how FSTL1 regulates the progression of GC. By consulting the literature, we found that COL5A1 is related to immune infiltration in the TME of GC patients [37], so we wondered whether FSTL1 was also involved. Therefore, we evaluated the stromal score, immune score and ESTIMATE score of GC patients from TCGA by the ESTIMATE algorithm and divided them into high and low groups (50% each). Then, we compared the expression levels of COL5A1 and FSTL1 in the two groups. As shown in Fig. 6a-f, patients with high stromal, immune and ESTIMATE scores also had high expression levels of COL5A1 and FSTL1. Figure 6g-l shows the correlations between the expression of COL5A1 or FSTL1 and those scores. The expression of COL5A1 or FSTL1 was significantly positively correlated with the three kinds of scores, and the correlations between the expression of FSTL1 and those scores were significantly higher than that of COL5A1, suggesting that FSTL1 might affect tumor immune infiltration in GC patients. Next, we divided GC patients from TCGA into high/low groups (50% each) according to the expression of FSTL1. The content of 22 types of immune cells in the TME of each GC patient was estimated by the CIBESORT algorithm. As shown in Fig. 6m, there were differences in the contents of 9 types of immune cells between the two groups. We then combined the expression of FSTL1 with the content of various immune cells for correlation analysis, from which we found that 11 types of immune cells were correlated with the expression of FSTL1 (Fig. 6n-p, Figure S7). Finally, we intersected these 11 types of immune cells and the 9 types of immune cells mentioned above, identifying 8 types of immune cells related to FSTL1 expression (Fig. 6q), namely, naive B cells, activated memory CD4 + T cells, follicular helper T cells, monocytes, M0 macrophages, M2 macrophages, resting mast cells and eosinophils. These results suggested that FSTL1 might promote the progression of GC by affecting tumor immunity.

Bioinformatics analysis can effectively reveal the potential mechanism in the tumorigenesis and development of GC and provide new targets and ideas for the treatment of GC. In the past few years, many studies have reported the screening of key miRNAs and their target genes in GC through various bioinformatics analyses. However, with the development and updating of bioinformatics, previous methods have become outdated and not generalizable to some extent because many studies conducted only a preliminary verification. For example, Wenqian Qi et al. constructed a prognosis prediction model of GC patients composed of three miRNAs (miR-126-3p, miR-143-5p and miR-1275), and their expression was verified in clinical samples [19]. It is worth noting that most miRNAs function by targeting coding genes, while focusing only on miRNAs and ignoring their interaction with target genes will lead to an overamplification of the functions of some miRNAs [3–4]. Combined analysis of miRNAs and their target genes can avoid this issue. In this study, we screened DEMs from TCGA and constructed a prognostic model including miR-137-3p, which could be used to evaluate the prognosis of GC patients. Then, we predicted the target genes of 11 SRDEMs in the model and intersected them with DEGs, identifying some hub target genes through enrichment analysis, the PPI network and WGCNA. Subsequently, we backtracked the miRNAs regulating hub genes to achieve the combined analysis of miRNAs and their target genes. Using this method, we found that miR-137-3p might play an important role in GC, and we carried out experimental verification and mechanistic exploration.

MiR-137 is located on chromosome 1p21.3 and in the non-coding gene AK094607 [38], the mature miRNA of which is miR-137-3p. The other mature miRNA, miR-137-5p, has not been verified [27]. Therefore, we selected miR-137-3p for further research. MiR-137-3p (i.e., miR-137) plays a tumor-suppressive role in a variety of cancers. Yajun Wang et al. reported the function of miR-137-3p in cell proliferation, apoptosis, invasion, angiogenesis and the treatment of glioma [39]. Feng Din et al. reported that miR-137-3p plays a tumor-suppressive role in pancreatic cancer by targeting MRGBP [40]. Meiying Luo et al. reported that miR-137-3p regulates ferroptosis by targeting the glutamine transporter SLC1A5 in melanoma cells [41]. In GC, many studies have also reported the antitumor effect of miR-137-3p. Qiaoyan Gu et al. revealed that the expression of miR-137-3p was related to tumor differentiation, N stage, TNM stage and overall survival rate through analyses of 154 clinical samples of GC, indicating that the expression of miR-137-3p was an independent prognostic index for GC patients [42]. Yan Cheng et al. reported that miR-137-3p was expressed at low levels in GC and inhibited activation of the PI3K/Akt signaling pathway by targeting cyclooxygenase-2 (COX-2) [43]. Ruta Steponaitiene et al. confirmed that epigenetic silencing of miR-137-3p is a common event in GC, and hypermethylation of its promoter leads to its low expression [44]. Studies have also reported that long non-coding RNAs (lncRNAs) or circular RNAs (circRNAs) sponge miR-137-3p through a ceRNA mechanism in GC [45–46]. However, no one has reported that mRNA can regulate miR-137-3p through a ceRNA mechanism, nor has it been reported that miR-137-3p can target COL5A1 or FSTL1. In this study, we explain and supplement the phenomenon in which miR-137-3p is inhibited in GC and its tumor-suppressive role: the low expression of miR-137-3p in GC leads to high expression of its target gene COL5A1, which in turn can sponge miR-137-3p through a ceRNA mechanism to form a positive feedback regulation to further reduce miR-137-3p in GC.

A miRNA can regulate multiple mRNAs, and a mRNA can also be regulated by multiple miRNAs. It was found that when one or several target mRNAs were downregulated, the number of free miRNAs increased. These miRNAs could silence other target genes, resulting in a decrease in the abundance of all target genes [47]. In view of this phenomenon, Pandolfi et al. proposed the ceRNA hypothesis based on the experimental results of their group and other researchers [48]. The ceRNA hypothesis holds that all RNA transcripts, including pseudogenes, lncRNAs, circRNAs or mRNAs, regulate each other through their common MREs. The more MREs they share, the greater is the degree of communication. In recent years, the ceRNA hypothesis has been confirmed to play a crucial role in various tumors, including GC. For example, Xianzi Yang et al. reported that LINC01133 could regulate the expression of APC and the Wnt/β-catenin pathway by sponging miR-106a-3p, thus inhibiting the progression of GC [49]. Xing Zhang et al. reported that circNRIP1 could regulate the AKT1/mTOR pathway by sponging miR-149-5p, thereby promoting the progression of GC [50]. Research on ceRNAs in recent years has focused on non-coding RNAs, such as lncRNAs and circRNAs [51–52]. However, when proposed by Pandolfi et al. in 2011, research on ceRNAs focused more on mRNAs competing for the same miRNA or even multiple miRNAs; that is, the combination of miRNA and mRNA will not only lead to the silencing of target genes but also reduce the effective concentration and actual activity of miRNAs. In this way, the expression level of one mRNA will affect the expression level and physiological activity of other mRNAs. In the present study, our focus returned to mRNAs. From the starBase database, we learned that COL5A1 and FSTL1 shared the binding sites of 82 miRNAs, including miR-137-3p, and there were significant positive correlations between their expression levels in 31 types of tumors in TCGA (the correlation coefficient was greater than 0.5 in 26 types of tumors). These results suggest that COL5A1 and FSTL1 may have a similar co-expression relationship in other types of tumors [36].

Collagen is the most abundant protein in mammals and the main component of the extracellular matrix (ECM) in the TME [53–54]. Abnormalities in collagen will affect the tumorigenesis and development of a variety of tumors as well as the infiltration of immune cells in the TME [37, 55–56]. COL genes participate in the interaction between the ECM and receptors and the adhesion pathway, playing an important role in tumor invasion and metastasis. Therefore, the prognostic value of multiple COL genes for various tumors has been confirmed [57–60]. Among the hub genes in this paper, COL genes accounted for a large proportion. In the WGCNA results, COL genes were mostly concentrated in the same module, suggesting that they had the same expression pattern and might play the same role in GC. In the present study, we only selected COL5A1 from COL genes for experimental verification, while other COL genes are worthy of further attention and in-depth research. Previous studies have reported that COL5A1 promotes the progression of various tumors, including GC, but they were limited to functional analysis of the COL5A1 protein [37, 61–63]. However, we also studied the function of COL5A1 mRNA and analyzed the tumor-promoting effect of COL5A1 from another perspective; that is, it regulates the expression of FSTL1 by competitively binding miR-137-3p and promotes the progression of GC as a ceRNA.

Follistatin like-1 (FSTL1) is a secretory glycoprotein that is involved in a variety of signaling pathways and biological processes, including the regulation of angiogenesis and the immune response. However, the role of FSTL1 in cancer has been controversial [64]. FSTL1 is downregulated in ovarian cancer, endometrioid cancer and kidney carcinoma, and a decrease in FSTL1 is related to a poor prognosis [65–66]. In contrast, in esophageal squamous cell carcinoma, gastric cancer, colorectal cancer, and breast cancer, increased expression of FSTL1 is associated with a poor prognosis [67–70]. The TME contains nutrients, chemokines, immune cells and other noncancerous components. It is the main location for tumor cells to interact with the host immune system [71]. In recent years, immunotherapy has rapidly developed into a promising and effective therapy for GC patients [72–73]. Different immune cell subsets in the TME have different effects on tumor progression and treatment [74]. Previous studies have shown that FSTL1 is closely related to the inflammatory response [64]. Chie Kudo-Saito et al. reported that FSTL1 inhibits antitumor immunity and promotes tumor invasion as well as bone metastasis [75]. In subsequent studies, they demonstrated that FSTL1 is a determinant of immune dysfunction mediated by mesenchymal stem cells (MSCs) and immunomodulatory cells. Blocking FSTL1 could significantly inhibit tumor progression and metastasis in a mouse tumor model with increased MSCs. The expression of DIP2A (the receptor of FSTL1) in tumor cells is the key to FSTL1-induced immune resistance [76]. Our results showed that the expression of FSTL1 was closely related to immunity in the TME of GC patients. Specifically, FSTL1 was positively correlated with the stromal, immune and ESTIMATE scores, suggesting that the expression of FSTL1 was negatively correlated with the purity of gastric cancer, while a low tumor purity often indicates a poor prognosis [77–78].

In addition, we also found that FSTL1 was significantly associated with a variety of immune cells in the TME, and macrophages attracted our attention. Tumor-associated macrophages (TAMs) are an important type of immunosuppressive cell. Generally, monocytes and macrophages (M0 macrophages) can be polarized into classically activated macrophages (M1 macrophages) or alternatively activated macrophages (M2 macrophages) under different stimuli, and their functions are almost antagonistic to each other [79–81]. M1 macrophages promote the immune response and normalize irregular tumor vascular networks so that cancer cells are sensitive to chemotherapy and radiotherapy [82]. In contrast, TAMs are usually similar to M2 macrophages, which play an immunosuppressive role in promoting tumor progression, metastasis and chemoresistance [83–86]. Our results showed that the expression of FSTL1 was negatively correlated with the proportion of M0 macrophages (Fig. 6n) and positively correlated with the proportion of M2 macrophages (Fig. 6o) and monocytes (Fig. 6p). Therefore, we speculated that FSTL1 could recruit monocytes and macrophages or promote the polarization of monocytes and M0 macrophages into M2 macrophages. This phenomenon might be related to the fact that FSTL1 is a secretory glycoprotein and the interaction between FSTL1 and its receptor DIP2A [76]. This hypothesis and the specific underlying mechanism require further study and verification.

Through a variety of bioinformatics analyses and experiments, we confirmed that miR-137-3p played a tumor-suppressive role in GC, and its target gene COL5A1 could reversely sponge miR-137-3p to relieve its targeted inhibition of FSTL1, which might promote the progression of GC by affecting immune infiltration in the TME. These results may provide new ideas for the treatment and future research of GC.

GC: Gastric cancer; miRNAs: MicroRNAs; RISC: RNA-induced silencing complex; TCGA: The Cancer Genome Atlas; GEO: Gene Expression Omnibus; KEGG: Kyoto Encyclopedia of Genes and Genomes; GO: Gene Ontology; DEMs: Differentially expressed miRNAs; DEGs: Differentially expressed genes; FC: Fold change; FDR: False discovery rate; SRDEMs: Survival related differentially expressed miRNAs; ROC: Receiver operating characteristic; DETGs: Differentially expressed target genes; STRING: Search Tool for the Retrieval of Interacting Genes; PPI: Protein-protein interaction; MCC: Maximal Clique Centrality; WGCNA: Weighted gene co-expression network analysis; DMEM: Dulbecco’s modifed Eagle medium; FBS: Fetal bovine serum; wt: Wild type; mut: Mutant; qRT-PCR: Quantitative real-time PCR; SDS–PAG: Sodium dodecyl sulfate-polyacrylamide gel electrophoresis; CCK-8: Cell Counting Kit-8; EdU: Ethynyldeoxyuridine; TME: Tumor microenvironment; COL genes: Collagen genes; MREs: MiRNA response elements; ceRNA: Competing endogenous RNA; lncRNAs: Long noncoding RNAs; circRNAs: Circular RNAs; ECM: Extracellular matrix; MSCs: Mesenchymal stem cells; TAMs: Tumor associated macrophages

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Availability of data and materials

All data in our study are available upon request.

Competing interests

The authors declare that they have no competing interests.

Funding

This work was supported by the National Natural Science Foundation of China (No.81602568, No.81773130, No.82172833), China Postdoctoral Science Foundation (No.2018M643009), the National Natural Science Foundation of Hunan Province (No.2019JJ50906, No.2019JJ80030, No.2018JJ2599, No.2020JJ4863, No.2020JJ4857) and Hunan Cancer Hospital Climb Plan (No.2020NSFC-A004).

Authors' contributions

Conceptualization: Changwei Lin; Data Curation: Ming Yang, Liang Li; Bioinformatics analysis: Ming Yang, Zhixing Lu; Methodology: Changwei Lin, Runliu Wu; Software: Ming Yang, Zhixing Lu; Supervision: Xiaorong Li, Gui Hu, Yi Zhang, Changwei Lin, Lihua Huang; Validation: Ming Yang, Zhixing Lu, Liang Li, Ma Min, Fei Long, Jing Chou, Ni Gong, Kaiyan Yang; Visualization: Ming Yang, Liang Li; Writing & Editing: Changwei Lin, Ming Yang. The authors read and approved the final manuscript.

Acknowledgements

We thank all individuals who participated in this work.

Sung H, Ferlay J, Siegel RL, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin. 2021;71(3):209–49.
Smyth EC, Nilsson M, Grabsch HI, van Grieken NC, Lordick F. Gastric cancer. Lancet. 2020;396(10251):635–48.
Bartel DP. MicroRNAs: target recognition and regulatory functions. Cell. 2009;136(2):215–33.
Filipowicz W, Bhattacharyya SN, Sonenberg N. Mechanisms of post-transcriptional regulation by microRNAs: are the answers in sight? Nat Rev Genet. 2008;9(2):102–14.
Zhang Z, Li Z, Li Y, Zang A. MicroRNA and signaling pathways in gastric cancer. Cancer Gene Ther. 2014;21(8):305–16.
Zhang Z, Pi J, Zou D, et al. microRNA arm-imbalance in part from complementary targets mediated decay promotes gastric cancer progression. Nat Commun. 2019;10(1):4397.
Ebrahimi S, Hashemy SI. MicroRNA-mediated redox regulation modulates therapy resistance in cancer cells: clinical perspectives. Cell Oncol (Dordr). 2019;42(2):131–41.
Weidle UH, Birzele F, Nopora A. microRNAs Promoting Growth of Gastric Cancer Xenografts and Correlation to Clinical Prognosis. Cancer Genomics Proteomics. 2021;18(1):1–15.
Shin VY, Chu KM. MiRNA as potential biomarkers and therapeutic targets for gastric cancer. World J Gastroenterol. 2014;20(30):10432–9.
So JBY, Kapoor R, Zhu F, et al. Development and validation of a serum microRNA biomarker panel for detecting gastric cancer in a high-risk population. Gut. 2021;70(5):829–37.
Hwang J, Min BH, Jang J, et al. MicroRNA Expression Profiles in Gastric Carcinogenesis. Sci Rep. 2018;8(1):14393.
He X, Zou K. MiRNA-96-5p contributed to the proliferation of gastric cancer cells by targeting FOXO3. J Biochem. 2020;167(1):101–8.
Luscombe NM, Greenbaum D, Gerstein M. What is bioinformatics? A proposed definition and overview of the field. Methods Inf Med. 2001;40(4):346–58.
Tomczak K, Czerwińska P, Wiznerowicz M. The Cancer Genome Atlas (TCGA): an immeasurable source of knowledge. Contemp Oncol (Pozn). 2015;19(1A):A68–77.
Barrett T, Troup DB, Wilhite SE, et al. NCBI GEO: mining tens of millions of expression profiles–database and tools update. Nucleic Acids Res. 2007;35(Database issue):D760–5.
Kanehisa M, Goto S. KEGG: kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30.
Ashburner M, Ball CA, Blake JA, et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat Genet. 2000;25(1):25–9.
Zhang Z, Dong Y, Hua J, et al. A five-miRNA signature predicts survival in gastric cancer using bioinformatics analysis. Gene. 2019;699:125–34.
Qi W, Zhang Q. Development and clinical validation of a 3-miRNA signature to predict prognosis of gastric cancer. PeerJ. 2021;9:e10462. Published 2021 Feb 3.
Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40.
Ritchie ME, Phipson B, Wu D, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47.
Ahmed FE, Vos PW, Holbert D. Modeling survival in colon cancer: a methodological review. Mol Cancer. 2007;6:15.
Chen W, Gao C, Liu Y, Wen Y, Hong X, Huang Z. Bioinformatics Analysis of Prognostic miRNA Signature and Potential Critical Genes in Colon Cancer. Front Genet. 2020;11:478.
Heagerty PJ, Zheng Y. Survival model predictive accuracy and ROC curves. Biometrics. 2005;61(1):92–105.
Kozomara A, Birgaoanu M, Griffiths-Jones S. miRBase: from microRNA sequences to function. Nucleic Acids Res. 2019;47(D1):D155–62.
Agarwal V, Bell GW, Nam JW, Bartel DP. Predicting effective microRNA target sites in mammalian mRNAs. Elife. 2015;4:e05005.
Wong N, Wang X. miRDB: an online resource for microRNA target prediction and functional annotations. Nucleic Acids Res. 2015;43(Database issue):D146–52.
Szklarczyk D, Morris JH, Cook H, et al. The STRING database in 2017: quality-controlled protein-protein association networks, made broadly accessible. Nucleic Acids Res. 2017;45(D1):D362–8.
Shannon P, Markiel A, Ozier O, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504.
Chin CH, Chen SH, Wu HH, Ho CW, Ko MT, Lin CY. cytoHubba: identifying hub objects and sub-networks from complex interactome. BMC Syst Biol. 2014;8(Suppl 4(Suppl 4):11.
Langfelder P, Horvath S. WGCNA: an R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559.
Yoshihara K, Shahmoradgoli M, Martínez E, et al. Inferring tumour purity and stromal and immune cell admixture from expression data. Nat Commun. 2013;4:2612.
Newman AM, Liu CL, Green MR, et al. Robust enumeration of cell subsets from tissue expression profiles. Nat Methods. 2015;12(5):453–7.
Bosson AD, Zamudio JR, Sharp PA. Endogenous miRNA and target concentrations determine susceptibility to potential ceRNA competition. Mol Cell. 2014;56(3):347–59.
Denzler R, McGeary SE, Title AC, Agarwal V, Bartel DP, Stoffel M. Impact of MicroRNA Levels, Target-Site Complementarity, and Cooperativity on Competing Endogenous RNA-Regulated Gene Expression. Mol Cell. 2016;64(3):565–79.
Li JH, Liu S, Zhou H, Qu LH, Yang JH. starBase v2.0: decoding miRNA-ceRNA, miRNA-ncRNA and protein-RNA interaction networks from large-scale CLIP-Seq data. Nucleic Acids Res. 2014;42(Database issue):D92–7.
Wei Z, Chen L, Meng L, Han W, Huang L, Xu A. LncRNA HOTAIR promotes the growth and metastasis of gastric cancer by sponging miR-1277-5p and upregulating COL5A1. Gastric Cancer. 2020;23(6):1018–32.
Yin J, Lin J, Luo X, et al. miR-137: a new player in schizophrenia. Int J Mol Sci. 2014;15(2):3262–71.
Wang Y, Chen R, Zhou X, et al. miR-137: A Novel Therapeutic Target for Human Glioma. Mol Ther Nucleic Acids. 2020;21:614–22.
Ding F, Zhang S, Gao S, et al. MiR-137 functions as a tumor suppressor in pancreatic cancer by targeting MRGBP. J Cell Biochem. 2018;119(6):4799–807.
Luo M, Wu L, Zhang K, et al. miR-137 regulates ferroptosis by targeting glutamine transporter SLC1A5 in melanoma. Cell Death Differ. 2018;25(8):1457–72.
Gu Q, Zhang J, Hu H, Tan YE, Shi S, Nian Y. Clinical Significance of MiR-137 Expression in Patients with Gastric Cancer After Radical Gastrectomy. PLoS One. 2015;10(11):e0142377.
Cheng Y, Li Y, Liu D, Zhang R, Zhang J. miR-137 effects on gastric carcinogenesis are mediated by targeting Cox-2-activated PI3K/AKT signaling pathway. FEBS Lett. 2014;588(17):3274–81.
Steponaitiene R, Kupcinskas J, Langner C, et al. Epigenetic silencing of miR-137 is a frequent event in gastric carcinogenesis. Mol Carcinog. 2016;55(4):376–86.
Li W, Duan J, Shi W, Lei L, Lv P. Long Non-Coding RNA NCK1-AS1 Serves an Oncogenic Role in Gastric Cancer by Regulating miR-137/NUP43 Axis. Onco Targets Ther. 2020;13:9929–39.
Lu Y, Li L, Li L, Wu G, Liu G. Circular RNA circHECTD1 prevents Diosbulbin-B-sensitivity via miR-137/PBX3 axis in gastric cancer. Cancer Cell Int. 2021;21(1):264.
Arvey A, Larsson E, Sander C, Leslie CS, Marks DS. Target mRNA abundance dilutes microRNA and siRNA activity. Mol Syst Biol. 2010;6:363.
Salmena L, Poliseno L, Tay Y, Kats L, Pandolfi PP. A ceRNA hypothesis: the Rosetta Stone of a hidden RNA language? Cell. 2011;146(3):353–8.
Yang XZ, Cheng TT, He QJ, et al. LINC01133 as ceRNA inhibits gastric cancer progression by sponging miR-106a-3p to regulate APC expression and the Wnt/β-catenin pathway. Mol Cancer. 2018;17(1):126.
Zhang X, Wang S, Wang H, et al. Circular RNA circNRIP1 acts as a microRNA-149-5p sponge to promote gastric cancer progression via the AKT1/mTOR pathway. Mol Cancer. 2019;18(1):20.
Ala U. Competing Endogenous RNAs, Non-Coding RNAs and Diseases: An Intertwined Story. Cells. 2020;9(7):1574.
Chan JJ, Tay Y. Noncoding RNA:RNA Regulatory Networks in Cancer. Int J Mol Sci. 2018;19(5):1310.
Ricard-Blum S. The collagen family. Cold Spring Harb Perspect Biol. 2011;3(1):a004978.
Lu P, Weaver VM, Werb Z. The extracellular matrix: a dynamic niche in cancer progression. J Cell Biol. 2012;196(4):395–406.
Xu S, Xu H, Wang W, et al. The role of collagen in cancer: from bench to bedside. J Transl Med. 2019;17(1):309.
Martins Cavaco AC, Dâmaso S, Casimiro S, Costa L. Collagen biology making inroads into prognosis and treatment of cancer progression and metastasis. Cancer Metastasis Rev. 2020;39(3):603–23.
Gao X, Zhong S, Tong Y, et al. Alteration and prognostic values of collagen gene expression in patients with gastric cancer under different treatments. Pathol Res Pract. 2020;216(3):152831.
Kauppila S, Stenbäck F, Risteli J, Jukkola A, Risteli L. Aberrant type I and type III collagen gene expression in human breast cancer in vivo. J Pathol. 1998;186(3):262–8.
Kita Y, Mimori K, Tanaka F, et al. Clinical significance of LAMB3 and COL7A1 mRNA in esophageal squamous cell carcinoma. Eur J Surg Oncol. 2009;35(1):52–8.
Wu YH, Chang TH, Huang YF, Huang HD, Chou CY. COL11A1 promotes tumor progression and predicts poor clinical outcome in ovarian cancer. Oncogene. 2014;33(26):3432–40.
Liu W, Wei H, Gao Z, et al. COL5A1 may contribute the metastasis of lung adenocarcinoma. Gene. 2018;665:57–66.
Feng G, Ma HM, Huang HB, et al. Overexpression of COL5A1 promotes tumor progression and metastasis and correlates with poor survival of patients with clear cell renal cell carcinoma. Cancer Manag Res. 2019;11:1263–74.
Ren W, Zhang Y, Zhang L, Lin Q, Zhang J, Xu G. Overexpression of collagen type V α1 chain in human breast invasive ductal carcinoma is mediated by TGF-β1. Int J Oncol. 2018;52(5):1694–704.
Mattiotti A, Prakash S, Barnett P, van den Hoff MJB. Follistatin-like 1 in development and human diseases. Cell Mol Life Sci. 2018;75(13):2339–54.
Chan QK, Ngan HY, Ip PP, Liu VW, Xue WC, Cheung AN. Tumor suppressor effect of follistatin-like 1 in ovarian and endometrial carcinogenesis: a differential expression and functional analysis. Carcinogenesis. 2009;30(1):114–21.
Liu Y, Han X, Yu Y, et al. A genetic polymorphism affects the risk and prognosis of renal cell carcinoma: association with follistatin-like protein 1 expression. Sci Rep. 2016;6:26689.
Lau MC, Ng KY, Wong TL, et al. FSTL1 Promotes Metastasis and Chemoresistance in Esophageal Squamous Cell Carcinoma through NFκB-BMP Signaling Cross-talk. Cancer Res. 2017;77(21):5886–99.
Peng X, Wang P, Li S, Jiang Y, Wu C. Follistatin–like protein 1 knockdown elicits human gastric cancer cell apoptosis via a STAT6–dependent pathway. Oncol Rep. 2019;42(6):2806–13.
Gu C, Wang X, Long T, et al. FSTL1 interacts with VIM and promotes colorectal cancer metastasis via activating the focal adhesion signalling pathway. Cell Death Dis. 2018;9(6):654.
Cheng S, Huang Y, Lou C, He Y, Zhang Y, Zhang Q. FSTL1 enhances chemoresistance and maintains stemness in breast cancer cells via integrin β3/Wnt signaling under miR-137 regulation. Cancer Biol Ther. 2019;20(3):328–37.
Nagarsheth N, Wicha MS, Zou W. Chemokines in the cancer microenvironment and their relevance in cancer immunotherapy. Nat Rev Immunol. 2017;17(9):559–72.
Vrána D, Matzenauer M, Neoral Č, et al. From Tumor Immunology to Immunotherapy in Gastric and Esophageal Cancer. Int J Mol Sci. 2018;20(1):13.
Coutzac C, Pernot S, Chaput N, Zaanan A. Immunotherapy in advanced gastric cancer, is it the future? Crit Rev Oncol Hematol. 2019;133:25–32.
Fridman WH, Galon J, Dieu-Nosjean MC, et al. Immune infiltration in human cancer: prognostic significance and disease control. Curr Top Microbiol Immunol. 2011;344:1–24.
Kudo-Saito C, Fuwa T, Murakami K, Kawakami Y. Targeting FSTL1 prevents tumor bone metastasis and consequent immune dysfunction. Cancer Res. 2013;73(20):6185–93.
Kudo-Saito C, Ishida A, Shouya Y, et al. Blocking the FSTL1-DIP2A Axis Improves Anti-tumor Immunity. Cell Rep. 2018;24(7):1790–801.
Aran D, Sirota M, Butte AJ. Systematic pan-cancer analysis of tumour purity. Nat Commun. 2015;6:8971.
Zhang C, Cheng W, Ren X, et al. Tumor Purity as an Underlying Key Factor in Glioma. Clin Cancer Res. 2017;23(20):6279–91.
Murray PJ. Macrophage Polarization. Annu Rev Physiol. 2017;79:541–66.
Piccolo V, Curina A, Genua M, et al. Opposing macrophage polarization programs show extensive epigenomic and transcriptional cross-talk. Nat Immunol. 2017;18(5):530–40.
Yang L, Zhang Y. Tumor-associated macrophages: from basic research to clinical application. J Hematol Oncol. 2017;10(1):58.
Jarosz-Biej M, Kamińska N, Matuszczak S, et al. M1-like macrophages change tumor blood vessels and microenvironment in murine melanoma. PLoS One. 2018;13(1):e0191012.
Caux C, Ramos RN, Prendergast GC, Bendriss-Vermare N, Ménétrier-Caux C. A Milestone Review on How Macrophages Affect Tumor Growth. Cancer Res. 2016;76(22):6439–42.
Chen Y, Zhang S, Wang Q, Zhang X. Tumor-recruited M2 macrophages promote gastric and breast cancer metastasis via M2 macrophage-secreted CHI3L1 protein. J Hematol Oncol. 2017;10(1):36.
Sousa S, Brion R, Lintunen M, et al. Human breast cancer cells educate macrophages toward the M2 activation status. Breast Cancer Res. 2015;17(1):101.
De Palma M, Lewis CE. Cancer: Macrophages limit chemotherapy. Nature. 2011;472(7343):303–4.

Download PDF

Version 1

posted

You are reading this latest preprint version

Identification of a Key Competing Endogenous RNA Axis, “COL5A1/miR-137-3p/FSTL1”, Associated With Gastric Cancer Progression by Integrated Bioinformatics Analyses and Experimental Verification

Status:

Version 1

Abstract

Background

Methods

Results

Conclusions

Figures

Background

Materials And Methods

Tumor data download and differential expression analysis

Cox proportional hazards regression model based on DEMs

Prediction of target genes of SRDEMs and functional enrichment analysis of DETGs

Construction and analysis of PPI networks with DETGs

Weighted gene co-expression network analysis

Cell culture

Plasmid construction

Transfection

RNA extraction and quantitative real-time PCR (qRT–PCR)

Western blotting

Cell Counting Kit-8 (CCK-8) assay

Ethynyldeoxyuridine (EdU) proliferation assay

Transwell migration/invasion assay

Wound healing assay

Dual-luciferase reporter assay

ESTIMATE and CIBESORT

Statistical analysis

Results

Cox proportional hazards regression model of DEMs

Prediction of target genes of SRDEMs and functional enrichment analysis of DETGs

Identification of hub genes in DETGs

WGCNA identification of the hub module

MiR-137-3p suppresses the proliferation, invasion and migration of GC cells

COL5A1 can reversely sponge miR-137-3p and upregulate the expression of FSTL1 through a ceRNA mechanism

FSTL1 is related to immune infiltration in the TME of GC patients

Discussion

Conclusion

Abbreviations

Declarations

References

Supplementary Files

Status:

Version 1