1. Background
Colorectal cancer (CRC) is a kind of disease with an increased morbidity year by year [Reference Akimoto, Ugai and Zhong1]. Many factors may contribute to the incidence of CRC, including genetics, environment, colon polyp, diet habits, and so on [Reference Keum and Giovannucci2, Reference Sawicki, Ruszkowska, Danielewicz, Niedzwiedzka, Arlukowicz and Przybylowicz3]. Bloody stool is usually the earliest and commonest symptom of CRC. When it comes to the advanced stage, there will be anemia, acute peritonitis, tumor metastasis, and so on [Reference Sawicki, Ruszkowska, Danielewicz, Niedzwiedzka, Arlukowicz and Przybylowicz3, Reference John, George, Primrose and Fozard4]. In the recent 50 years, although much progress has been made in medical conditions, the 5 years survival rate of CRC still hovers around 50% due to its high recurrence [Reference Baidoun, Elshiwy and Elkeraie5]. At present, it is still important to explore new, safe, and effective treatment options for CRC patients.
In recent years, studies have shown that RNA with different lengths interacts with each other, especially among microRNA (miRNA), long noncoding RNA (lncRNA), and messenger RNA (mRNA), thus forming a mutual regulatory network of lncRNA-miRNA-mRNA that is important to show the interactions between RNA molecules and explain their functions [Reference Ye, Yang, Zhao, Song, Wang and Zheng6–Reference Zhou, Zhang and Sun8]. With relatively high specificity in cells and tissues and a wide existence in various human tissues from serum and stool, miRNAs are good noninvasive biological agents for the measurement of precancerous and cancerous lesions and therapeutic targets in many diseases [Reference Link and Kupcinskas9, Reference Pereira, Marques, Soares, Carreto and Santos10]. Thus, searching for miRNAs in CRC progression and its mRNA targets could be a method to find more effective diagnostic biomarkers and treatment targets for CRC patients. Previously, some researchers analyzed the correlation between abnormally expressed miRNAs and CRC by constructing a differential miRNA bioinformatics expression map and predicted that it would affect CRC by interfering with downstream signaling pathways [Reference Fonseca, Ramalhete and Mestre11].
β-Elemene is a bioactive compound extracted from natural plants. As an emerging antitumor Chinese medicine, β-elemene has been applied in medical therapies because it could boost the effects of other antitumor Chinese medicines on lung cancer [Reference Feng, Li and Liu12], breast cancer [Reference Pan, Wang and Huang13], colorectal cancer [Reference Chen, Li and Zhang14], glioma [Reference Zhang, Chen and Yao15], osteosarcoma [Reference Zhang and Guo16], and so on and decrease the side effects of radiotherapy and chemotherapy [Reference Zhan, Li, Guo, Du, Wu and Zhang17]. Previously, Chen P et al. found that by triggering ferroptosis and suppressing EMT, the β-elemene and cetuximab therapies were sensitive to KRAS-mutated CRC cells [Reference Chen, Li and Zhang14]. The other report indicates that β-elemene could suppress the occurrence and development of CRC in the nude mouse by affecting transplanted tumor proliferation, apoptosis, and autophagy process [Reference Wang18]. However, the exact mechanism between β-elemene and CRC has not been studied yet. Therefore, we intend to analyze the regulatory lncRNA-miRNA-mRNA mechanism of β-elemene-miRNA-mRNA in CRC by whole quasi-transcriptome sequencing.
This study plans to investigate the mechanism of the lncRNA-miRNA-mRNA network in the process of β-elemene inhibiting CRC and to further find promising therapeutic targets and prognostic biomarkers for CRC. Our findings will bring a new perspective to the CRC study and lay the therapeutic basis for the clinical treatment of CRC patients.
2. Materials and Methods
2.1. Cell Culture and Treatment
The HCT-116 cells were from the China National Collection of Authenticated Cell Culture (Shanghai, China). The cells were cultured in Dulbecco’s modified Eagle’s medium (DMEM) obtained from Carlsbad, CA, USA, with 2 ml glutamine and 10% Fetal bovine serum (FBS) obtained from Gibco, USA, in a humidified environment with 5% CO2 at 37°C. β-Elemene is available from Sigma. Cells treated with β-elemene were further used for RNA sequencing.
2.2. Small RNA Sequencing
We selected the control (no treatment) and case groups (β-elemene treatment) for RNA sequencing and pertinent analysis. Then, the RNA sequencing was performed by Origin-Biotech Inc. (Shanghai, China).
2.3. Identification of Differentially Expressed (DE) lncRNA, miRNA, and mRNA
Based on the RNA sequencing data, DElncRNA, DEmiRNA, and DEmRNA were identified by edgeR with a threshold of P < 0.05 and the absolute values of |log2 fold change (FC)|>1, respectively.
2.4. Functional Enrichment Analysis on the DElncRNA, DEmiRNA, and DEmRNA
Gene Ontology (GO) contains a set of related terms or concepts, and it describes the understanding of biology from three domains: cellular component (CC), biological process (BP), and molecular function (MF). The Kyoto Encyclopedia of Gene and Genome (KEGG) is a predictive calculation tool on the basis of relevant knowledge. Given a complete set of genes on a chromosome, it can anticipate the relevance of protein interaction networks in a variety of biological functions. Herein, the R package (v 3.5.1) and clusterProfiler were applied for this analysis.
2.5. Protein-Protein-Interaction (PPI) Network on the DElncRNA, DEmiRNA, and DEmRNA
Then, the Search Tool for the Retrieval of Interacting Genes (STRING v 11.5, https://cn.string-db.org/) tool was adopted to map PPI networks for all DElncRNA, DEmiRNA, and DEmRNA with a composite interaction score ≥0.4. Cytoscape was adopted to visualize the PPI networks. The connectivity among DElncRNA (n = 10), DEmiRNA (n = 3), and DEmRNA (n = 53) with high degrees was demonstrated by an independent PPI network.
2.6. The LASSO Analysis and Prognostic Signature Model Construction on the above 53 DEmRNAs
After the screening of the top 53 DEmRNAs from the sequencing data of β-elemene-treated CRC cells, we performed the least absolute shrinkage and selection operator (LASSO) regression model by the “glmnet” package of the R language. The relationship between the partial likelihood deviation and log (lambda, λ) was plotted, and the important parameters of the signature model were obtained. Moreover, we downloaded clinical samples of CRC from The Cancer and Genome Atlas (TCGA, https://tcga-data.nci.nih.gov/tcga) dataset and analyzed the prognostic values of mRNAs in clinical. First, CRC samples were classified into high-risk (n = 309) and low-risk (n = 310) groups in the risk score analysis, and the survival status of CRC samples corresponding to 15 different mRNAs was displayed. Next, the overall survival (OS) probability between two groups was compared by the Kaplan–Meier (KM, https://kmplot.com/analysis/) survival curve. Finally, the receiver operating characteristic (ROC) curve was used to compare the area under the curve (AUC) values of 1-, 3-, and 5-year survival rates of CRC, and the HR with 95% CI was calculated.
2.7. The Establishment of Prognostic Nomogram in CRC
To further verify the key mRNAs on the prognosis of CRC patients, the “forestplot” package was used to draw forest plots to display the results of univariate and multivariate Cox regression analyses on the above 15 mRNAs, and the corresponding P value and hazard ratio (HR) with 95% CI were calculated. Then, 7 mRNAs related to CRC prognosis were identified. Then, the corresponding nomogram was constructed by the “rms” package to predict the 1-, 3-, and 5-year survival rates with 7 prognostic mRNAs for CRC patients. The closer the dots are to the calibration curve, the higher the accuracy of this nomogram.
2.8. The Expression Verification and KM Analysis on the 7 Prognostic mRNAs
Based on the above findings, we detected their expressions in CRC tumor (n = 620) and normal (n = 10) samples in the TCGA database to identify the key mRNAs related to CRC. Then, the relations between the key mRNA expressions and OS and progression-free survival (PFS) curves were studied.
3. Results
3.1. DElncRNA, DEmiRNA, and DEmRNA from β-Elemene-Treated Samples
First, we analyzed the DElncRNA, DEmiRNA, and DEmRNA expression profiling data from β-elemene-treated samples. According to the RNA-seq analysis, 607 upregulated and 599 downregulated DElncRNAs were identified and displayed by volcano and heat maps (Figures 1(a) and 1(b)), and the top 5 upregulated and downregulated DElncRNAs are demonstrated in Table 1. As shown in Figures 1(c) and 1(d), we identified 12 downregulated and 24 upregulated DEmiRNAs, and the top 5 of each group are shown in Table 2. Besides, 3153 downregulated and 3248 upregulated DEmRNAs are displayed in Figures 1(e) and 1(f), and the top 5 DEmRNAs of them are also shown in Table 3.
FIGURE 1: Identification of DElncRNA, DEmiRNA, and DEmRNA induced by β-elemene treatment in CRC. (a)-(b) Volcano map and heat map of lncRNAs. Blue stands for downregulated DEGs, red stands for upregulated DEGs, and gray stands for no significant change in genes. (c)-(d) Volcano map and heat map of miRNAs. Green stands for downregulated DEGs, orange stands for upregulated DEGs, and gray stands for no significant change in genes. (e)-(f) Volcano map and heat map of mRNAs. Yellow stands for downregulated DEGs, purple stands for upregulated DEGs, and gray stands for no significant change in genes.
TABLE 1: Top 5 increased and decreased DElncRNAs.
TABLE 2: Top 5 increased and decreased DEmiRNAs.
TABLE 3: Top 5 increased and decreased DEmRNAs.
3.2. The Results of GO Term and KEGG Pathway Enrichment Analyses
In GO analysis, DElncRNAs were enriched in positive regulation of muscle hyperplasia, regulation of B cell activation, lung growth (BP), and mRNA binding involved in posttranscriptional gene silencing (MF, Figure 2(a)). DEmiRNAs were enriched in negative regulation of transcription by RNA polymerase II, positive regulation of gene expression (BP), transcription regulator complex, chromatin, heteromeric SMAD protein complex (CC), DNA-binding transcription factor binding, and sequence-specific DNA binding (MF, Figure 2(c)). DEmRNAs were enriched in positive regulation of phospholipase C activity (BP), basolateral plasma membrane (CC), clathrin-coated vesicle membrane, and RNA polymerase III complex (MF, Figure 2(e)). In the KEGG pathway analysis, DElncRNAs were related to the CCL18 signaling pathway, ncRNAs involved in STAT3 signaling in hepatocellular carcinoma, lncRNA mediated mechanisms of therapeutic resistance (Figure 2(b)), DEmiRNAs were related to the integrated breast cancer pathway, the TGF-beta signaling pathway, senescence, and autophagy in cancer (Figure 2(d)), and DEmRNAs were related to the ectoderm differentiation, the leptin signaling pathway, and so on (Figure 2(f)).
FIGURE 2: The GO and KEGG enrichment analyses on DElncRNAs, DEmiRNAs, and DEmRNAs. (a) The GO analysis of lncRNAs. (b) The KEGG pathway enrichment analysis of lncRNAs. (c) The GO analysis of miRNAs. (d) The KEGG pathway enrichment analysis of miRNAs. (e) The GO analysis of mRNAs. (f) The KEGG pathway enrichment analysis of mRNAs.
3.3. The lncRNA-miRNA-mRNA Coexpression Network Analysis
Next, we constructed a PPI network of the DERNAs with the Cytoscape tool (Figure 3(a)). According to the degrees in each group, DERNAs with high degrees (composite interaction score ≥0.4) in each group were selected for the construction of the lncRNA-miRNA-mRNA network, including 10 DElncRNAs, 3 DEmiRNAs, and 53 DEDmRNAs. As demonstrated in Figure 3(b), miR-7109-3p, miR-4506, and miR-3182 were closely related to ambient DElncRNAs and DEmRNAs, which indicated that miR-7109-3p, miR-4506, and miR-3182 had the potential to be clinical biomarkers or therapeutic targets in CRC.
FIGURE 3: The lncRNA-miRNA-mRNA coexpression network analysis. (a) The PPI network of the DERmNAs. (b) The lncRNA-miRNA-mRNA network of RNAs with top degrees in each category.
Nodes represent genes, and edges represent interactions between genes. The green nodes represent the lncRNA, the red nodes represent the miRNA, and the purple nodes represent the mRNA.
3.4. The LASSO Analysis and Prognostic Signature Model Construction of the 53 DEmRNAs
The LASSO regression analysis was first performed on the above 53 DEmRNAs to determine the optimal parameter number of the signature model (Figures 4(a) and 4(b), each line represented a DEmRNA). Then, CRC samples in the TCGA database were divided into high- (n = 309) and low-risk groups (n = 310), and it was found that the 15 mRNAs, M1AP, C1QTNF4, SYT1, SEPTIN5, GRID2IP, PAQR6, NIN, DTX1, TNFRSF11 A, HOXD13, ZNF20, ALPG, RIMS3, SLC16A8, and TNNT1, had the potential to be prognostic biomarkers in CRC (Figure 4(c)). In the KM survival curve, the median time of the two groups was 3.1 years, and the OS probability of the high-risk group was poor, and its HR was 3.175 (>1), indicating that the model was a risk model (Figure 4(d)). In addition, the AUC values in the ROC curve were all greater than 0.7, suggesting a good prognostic prediction ability of the model (Figure 4(e)).
FIGURE 4: The LASSO analysis and prognostic signature model construction on the 53 mRNAs. (a)-(b) The LASSO regression mode analysis on the 53 mRNAs. (c) The risk score evaluation of CRC patients (top). The corresponding survival time of each patient (middle). The Z-score of expression of 15 candidate prognostic mRNAs (bottom). (d) Prognostic analysis of high-risk and low-risk groups. (e) The ROC curve analysis.
3.5. The Identification of mRNAs with Prognostic Value in CRC
To investigate the mRNAs with prognostic values in CRC, we performed univariate and multivariate Cox regressions on the above 15 DEmRNAs. Under the premise of P < 0.05, it was found that ALPG, DTX1, HOXD13, RIMS3, SLC16A8, SYT1, and TNNT1 had a close relationship with the prognosis of CRC patients (Figures 5(a) and 5(b)). Subsequently, a nomogram on these genes was designed to predict the 1-, 3-, and 5-year survival rates of CRC patients (Figure 5(c)). As shown in Figure 5(d), the different dots were relatively close to the calibration curve, indicating that the prognostic model had good predictive ability.
FIGURE 5: The prognostic nomogram with prognostic mRNAs. (a)-(b) Univariate and multivariate Cox regression analyses. (c) The prognostic nomogram with prognostic mRNAs. (d) The dotted line is the ideal calibration curve of the nomogram.
3.6. RIMS3 and SLC16A8 Were the Key mRNAs Related to CRC in β-Elemene Treatment
To further explore the exact mRNAs related to CRC in β-elemene treatment, we detected the expressions of the 7 prognostic mRNAs in CRC tumor and normal samples. After detection, it was found that RIMS3 was significantly increased in normal samples, while SLC16A8 was substantially upregulated in tumor samples (Figures 6(a) and 6(d)), which indicated that RIMS3 was a suppressor gene while SLC16A8 was an oncogene in CRC. In the OS and PFS analyses on them, low expression of RIMS3 represented a poor survival rate, while low expression of SLC16A8 represented a relatively good survival rate (Figures 6(b), 6(c), 6(e), 6(f)).
FIGURE 6: RIMS3 and SLC16A8 were the key mRNAs related to CRC in β-elemene treatment. (a) RIMS3 significantly downregulated in CRC samples. (b)-(c) The Kaplan–Meier survival curve showing that lower expression levels of RIMS3 were significantly associated with shorter OS and PFS in CRC. (d) SLC16A8 upregulated in tumor samples. (e)-(f) The Kaplan–Meier survival curve showing that higher expression levels of RIMS3 were significantly associated with shorter OS and PFS in CRC. ∗P < 0.01, ****P < 0.0001.
4. Discussion
In recent years, the lncRNA-miRNA-mRNA network has been applied in the mechanism study of various diseases. For example, Li et al. found a novel lncRNA-miRNA-mRNA signature that predicted recurrence and disease-free survival in cervical cancer [Reference Li, Tian and Guo19]. Cao et al. found through miR-206, lncRNA-RMRP increased bladder cancer development [Reference Cao, Liu, Huang, Yue and Xi20]. Besides, some researchers demonstrated lncRNA-CDC6 might act as a ceRNA to facilitate breast cancer development by sponging miRNA-215 [Reference Kong, Duan and Sang21]. In this study, we also tried to investigate the exact mechanism in β-elemene inhibiting CRC progression from the aspect of the lncRNA-miRNA-mRNA network.
Through deep sequencing, high-throughput screening, and microarray technology, many dysregulated miRNAs were found in cancer cells. Due to the tissue specificity of miRNA regulation, miRNAs have much potential to be biomarkers in cancer diagnosis, treatment, and prognosis [Reference Xu, Geng, Yuan, Sun, Liu and Chen22]. In this study, we identified 1026 DElncRNAs, 6401 DEmRNAs, and 36 DEmiRNAs in β-elemene-treated HCT116 cells using RNA sequencing. Then, these DERNAs were analyzed by functional analyses and PPI networks. Enrichment analysis showed that DElncRNAs were associated with the CCL18 signaling pathway. Ruixue Yuan et al. showed that high CCL18 levels can be an independent biomarker for predicting better survival in CRC patients [Reference Yuan, Chen and He23]. The DEmiRNAs were related to the TGF-beta signaling pathway. TGF-beta signaling is one of the important cellular pathways that plays a key role in tissue maintenance. Nan et al. mentioned that LINC00941 plays an important role in metastatic CRC by activating the TGF-β/SMAD2/3 axis [Reference Wu, Jiang and Liu24]. DEmRNAs are related to the leptin signaling pathway. Leptin is an activator of cell proliferation, an antiapoptotic molecule, and an inducer of cancer stem cells in many cell types. Studies have shown that the stringent binding affinity of leptin/Ob-R and the overexpression of leptin/Ob-R and its targets in cancer cells make it a unique drug target for the prevention and treatment of CRC, especially in obese colorectal patients [Reference Zhou, Tian, Gong, Guo and Luo25]. In the PPI network, we selected the top RNAs in each group and observed their mutual relations, including miR-7109-3p, miR-4506, and miR-3182. At present, there are few reports on miR-7109-3p. One study highlights LINC00973-miR-7109-Siglec-15 function in the immune evasion of clear-cell renal cell carcinoma [Reference Karpov, Spirin and Zheltukhin26]. Nagy discovered miR-4506 was differentially expressed in adenoma compared to normal both in CRC tissue and plasma samples [Reference Nagy, Wichmann and Kalmar27]. Moreover, it was found that the biological functions and potential mechanisms of miR-3182 were altered after Ganoderic acid treatment in CRC [Reference Chen, Wan and Zeng28]. The 3 miRNAs all have the potential to be clinical biomarkers in CRC.
Besides, we focused on the mRNA associated with β-elemene and CRC. Through LASSO, risk score evaluation, and Cox analyses, 7 prognostic mRNAs, ALPG, DTX1, HOXD13, RIMS3, SLC16A8, SYT1, and TNNT1, were identified and applied to design the predictive nomogram for CRC prognosis, and the nomogram was proved to have a good prediction ability. In the previous research, few studies were conducted on these genes. In a study on TNNT1, researchers demonstrated that TNNT1, a target of miR-873, is related to CRC prognosis [Reference Chen, Wang and Wang29]. Herein, our findings show that the 7 mRNAs could be the prognostic biomarkers for CRC patients.
Furthermore, two key mRNAs related to CRC were identified by expression and KM analysis, namely, RIMS3 and SLC16A8. According to the public database, RIMS3 was a suppressor gene in CRC, while SLC16A8 was an oncogene in CRC. Regulating synaptic membrane exocytosis 3 (RIMS3) is supposed to enable transmembrane transporter binding activity and participate in the calcium ion-regulated exocytosis of neurotransmitters. Currently, it is suggested that RIMS3 is a gene related to autism. Functional investigations on RIMS3 variations like p.E177 A deepen our understanding of synaptic proteins’ function in autism pathogenesis [Reference Kumar, Sudi and Babatz30]. Additionally, solute carrier family 16 member 8 (SLC16A8), alias MCT3 and REMP, belongs to a family of proton-coupled monocarboxylate transporters regulating lactate transport across cell membranes [Reference Yoon, Donoso and Philp31]. The study of Klipfel L demonstrates the loss of transport function of the hypomorphic allele in the SLC16A8 gene and offers a methodological framework for the investigation of other SLC16A8 alleles linked to age-related macular degeneration [Reference Klipfel, Cordonnier and Thiebault32]. Currently, there are no reports about these two genes in CRC development. Based on our study, we conclude that RIMS3 and SLC16A8 are the key mRNAs involved in the β-elemene treatment in CRC progression, which certainly reveals the inhibitory mechanism of β-elemene in CRC. The study still has limitations. First, the expression profiles given by differential expression need to be validated with clinical samples. Second, further experimental studies are required to confirm the function of the lncRNA-miRNA-mRNA network in the inhibition of CRC progression by β-elemene.
To sum up, this study reveals the lncRNA-miRNA-mRNA network in the inhibitory function of β-elemene in CRC. These identified differentially expressed mRNAs (RIMS3 and SLC16A8) and miRNAs (miR-7109-3p, miR-4506, and miR-3182) all could be promising clinical biomarkers or targets in CRC diagnosis, treatment, and prognosis.
Data Availability
The datasets used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This study was funded by the Natural Science Foundation of Anhui Province (2008085QH388) and Clinical Research Project of the Anhui University of Chinese Medicine (2021efylc02).
Open access
