Abstract

Object. β-Elemene is an emerging antitumor Chinese medicine, but the exact mechanism of action of β-elemene in colorectal cancer (CRC) remains unclear. This study aimed to explore the mechanism of the lncRNA-miRNA-mRNA network in the process of β-elemene inhibiting CRC. Methods. RNA sequencing was performed on CRC cells from the control group (untreated) and the case group (β-elemene-treated). According to the sequencing data, we screened the differentially expressed (DE) lncRNAs, miRNAs, and mRNAs and then analyzed them by functional enrichment analyses. Through the lncRNA-miRNA-mRNA network, the key miRNAs and mRNAs involved in the process of β-elemene inhibiting CRC were further identified. Results. Totally, 607 upregulated and 599 downregulated DElncRNAs, 12 downregulated and 24 upregulated DEmiRNAs, and 3153 downregulated and 3248 upregulated DEmRNAs were identified. Through the lncRNA-miRNA-mRNA network, 3 miRNAs (miR-7109-3p, miR-4506, and miR-3182), 7 prognostic mRNAs (ALPG, DTX1, HOXD13, RIMS3, SLC16A8, SYT1, and TNNT1), and 2 key mRNAs (RIMS3 and SLC16A8) were determined to participate in the inhibitory mechanism of β-elemene in CRC. Conclusion. This study revealed for the first time that the lncRNA-miRNA-mRNA network is involved in the regulation of β-elemene in CRC, and these identified miRNAs and mRNAs could be new clinical prognostic biomarkers and therapeutic targets for CRC patients.

1. Background

Colorectal cancer (CRC) is a kind of disease with an increased morbidity year by year [1]. Many factors may contribute to the incidence of CRC, including genetics, environment, colon polyp, diet habits, and so on [2, 3]. 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 [3, 4]. 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 [5]. 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 [68]. 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 [9, 10]. 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 [11].

β-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 [12], breast cancer [13], colorectal cancer [14], glioma [15], osteosarcoma [16], and so on and decrease the side effects of radiotherapy and chemotherapy [17]. 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 [14]. 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 [18]. 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 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 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.

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)).

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.

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)).

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 , 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.

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)).

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 [19]. Cao et al. found through miR-206, lncRNA-RMRP increased bladder cancer development [20]. Besides, some researchers demonstrated lncRNA-CDC6 might act as a ceRNA to facilitate breast cancer development by sponging miRNA-215 [21]. 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 [22]. 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 [23]. 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 [24]. 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 [25]. 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 [26]. Nagy discovered miR-4506 was differentially expressed in adenoma compared to normal both in CRC tissue and plasma samples [27]. Moreover, it was found that the biological functions and potential mechanisms of miR-3182 were altered after Ganoderic acid treatment in CRC [28]. 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 [29]. 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 [30]. 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 [31]. 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 [32]. 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).