Abstract

Hepatocellular carcinoma (HCC) is a primary cause of cancer-related death in the world. Despite the fact that there are many methods to treat HCC, the 5-year survival rate of HCC is still at a low level. Emodin can inhibit the growth of HCC cells in vitro and in vivo. However, the gene regulation of emodin in HCC has not been well studied. In our research, RNA sequencing technology was used to identify the differentially expressed genes (DEGs) in HepG2 cells induced by emodin. A total of 859 DEGs were identified, including 712 downregulated genes and 147 upregulated genes in HepG2 cells treated with emodin. We used DAVID for function and pathway enrichment analysis. The protein-protein interaction (PPI) network was constructed using STRING, and Cytoscape was used for module analysis. The enriched functions and pathways of the DEGs include positive regulation of apoptotic process, structural molecule activity and lipopolysaccharide binding, protein digestion and absorption, ECM-receptor interaction, complement and coagulation cascades, and MAPK signaling pathway. 25 hub genes were identified and pathway analysis revealed that these genes were mainly enriched in neuropeptide signaling pathway, inflammatory response, and positive regulation of cytosolic calcium ion concentration. Survival analysis showed that LPAR6, C5, SSTR5, GPR68, and P2RY4 may be involved in the molecular mechanisms of emodin therapy for HCC. A quantitative real-time PCR (qRT-PCR) assay showed that the mRNA levels of LPAR6, C5, SSTR5, GPR68, and P2RY4 were significantly decreased in HepG2 cells treated with emodin. In conclusion, the identified DEGs and hub genes in the present study provide new clues for further researches on the molecular mechanisms of emodin.

1. Introduction

Hepatocellular carcinoma (HCC) is a severe disease of the digestive system and is the sixth most common and second leading cause of death due to cancer. While the early stage of HCC can be successfully treated by surgical resection and liver transplantation, most HCC cases are diagnosed only at the late stage, when treatment options are more limited [13]. Targeted therapies using sorafenib, regorafenib, and lenvatinib, as well as immunotherapy with a PD-1 inhibitor (nivolumab), have been approved for the treatment of advanced-stage HCC [47]. However, there remains an urgent need for the identification and development of novel therapeutic agents and strategies for the treatment of HCC.

Emodin (1,3,8-trihydroxy-6-methylanthraquinone; Figure 1) is an active ingredient derived from Polygonum cuspidatum [8], Rheum palmatum [9], Cassia occidentalis [10], and Polygonum multiflorum [11] and has been used in China for many centuries. The role of emodin as an anticancer drug has been previously described. Modern pharmacological studies have revealed that emodin exhibits various biological activities, such as apoptosis-inducing and antiproliferative effects in breast cancer [12], pancreatic cancer [13], prostate cancer [14], gastric cancer [15], lung carcinoma [16], colon cancer [17], and HCC [18]. Furthermore, emodin can inhibit metastasis, invasion, and migration in HCC and breast cancer [19, 20]. Dong et al. showed that emodin induces apoptosis in human HCC [21]. Hsu et al. confirmed that emodin inhibited the growth of hepatoma cells [22]. Numerous studies have confirmed the therapeutic effect of emodin on liver cancer. Thus, it is necessary to identify the key genes associated with emodin in HepG2 cells by conducting comprehensive bioinformatic analysis.

High-throughput technologies, such as transcriptome, protein, metabolite, and RNA sequencing, are high accuracy tools that can be used to identify biomarkers for the treatment, diagnosis, and prognosis of various diseases [23]. RNA sequencing (RNA-seq) uses deep-sequencing technologies to provide precise information regarding transcription profiles. The use of RNA-seq in analyzing the effects of drug treatments presents significant advantages including the identification of differentially expressed genes (DEGs) associated with the drug. Network and functional enrichment analyses are also beneficial in understanding the molecular mechanisms underlying drug action.

Although emodin exhibits good clinical efficacy, its gene regulatory mechanisms in liver cancer cells have not been systematically elucidated. Therefore, it is necessary to measure the expression levels of DEGs in cancer cells after treatment with emodin and to systematically analyze the functions of these genes. To overcome the aforementioned issue, the transcriptomes of emodin-treated HepG2 cells were profiled using RNA-seq method. The DEGs induced by emodin treatment were then examined in more detail using a series of analysis tools. The hub genes were extracted, and their corresponding expression levels were compared. A series of survival analyses was then conducted to determine whether the hub genes are correlated with poor prognosis. The relationship between the hub genes and tumor progression in patients with HCC was analyzed. Finally, statistical analysis of functional DEGs was performed, and their potential possible contributions to the anticancer effects of emodin were discussed.

2. Materials and Methods

2.1. Drug

Emodin was purchased from the Chinese Medicine Center in Beijing, dissolved in dimethylsulfoxide (DMSO) at a concentration of 100 mM, and stored at −20°C. The compound was diluted in the appropriate medium to 25, 50, 75, and 100 μM immediately before use. The final concentration of DMSO was <0.1%.

2.2. Cell Culture

Human hepatocellular carcinoma HepG2 cells were purchased from cell bank of the institute of Biochemistry and Cell Biology, Shanghai Institutes for Biological Sciences, Chinese Academy of Science (http://www.cellbank.org.cn/). Cells were cultured in RPMI1640 medium (GIBCO, Grand Island, NY, USA), with 10% fetal bovine serum (FBS, Gemini, US), 100 U/ml penicillin, and 100 mg/ml streptomycin in a humidified atmosphere with 5% CO2 at 37°C (Thermo Fisher Science, MA, USA). The cells with 80% confluence were treated by emodin of different concentrations. In this study, we have pooled 3 biological repeats into one RNA-seq experiment to achieve the same amount of total RNA content. More specifically, at the cell culturing section, cancer cells in 6 different petri dishes were divided into two equal groups. One group was used as emodin treated group, and the other was used as control/untreated group.

2.3. Cell Proliferation Assay

Cell Counting Kit-8 assay was used to measure cell proliferation. HepG2 cells were seeded in 96-well plates at 5 × 103 cells/well and incubated at 37°C in complete medium for 24 h before being treated with increasing concentrations of emodin for up to 72 h. Then the cells were incubated with 10 μl CCK8 at 37°C for 2 h. Absorbance at 450 nm was determined using a microplate reader as recommended by the manufacturer. All experiments were performed in triplicate and repeated at least three times with essentially similar results. The inhibition rate was calculated as follows: inhibition rate (%) = [average OD value (control)-average OD value (medication)]/average OD value (control) × 100%. The IC50 value was calculated on the nonlinear regression fit method by the GraphPad Prism software.

2.4. RNA-Seq and Analysis

We have pooled 3 biological repeats into one RNA-seq experiment to achieve the same amount of total RNA content. More specifically, at the cell culturing section, cancer cells in 6 different petri dishes were divided into two equal groups. One group was used as emodin treated group, and the other was used as control/untreated group. Briefly, after being treated with emodin for 72 h, total RNA was extracted using Trizol reagent (Invitrogen, Carlsbad, CA, USA). We use Agilent 2100 Bio analyzer (Agilent RNA 6000 Nano Kit) to do the total RNA sample QC: RNA concentration, RIN value, 28S/18S, and the fragment length distribution. We filter the low quality reads (more than 20% of the bases qualities are lower than 10), reads with adaptors, and reads with unknown bases (N bases more than 5%) to get the clean reads. Then we assembled those clean reads into Unigenes, followed with Unigene functional annotation, SSR detection and calculate the Unigene expression levels and SNPs of each sample. Finally, we identify DEGs (differential expressed genes) between samples and do clustering analysis and functional annotations. We use internal software SOAPnuke to filter reads, as follows: (1) remove reads with adaptors; (2) remove reads in which unknown bases (N) are more than 10%; (3) remove low quality reads (we define the low quality read as the percentage of base whose quality is lesser than 15 and is greater than 50% in a read). After filtering, the remaining reads are called “Clean Reads” and stored in FASTQ format. We use HISAT (Hierarchical Indexing for Spliced Alignment of Transcripts) to do the mapping step because HISAT is much faster and sensitive and is a high accuracy analysis software [24]. We mapped clean reads to reference using Bowtie2 [25], and then calculated gene expression level with RSEM [26]. RSEM is a software package for estimating gene and isoform expression levels from RNA-Seq data. Then, we calculate Pearson correlation between all samples using cor, perform hierarchical clustering between all samples using hclust, perform PCA analysis with all samples using princomp, and draw the diagrams with ggplot2 with fuctions of R. We detect DEGs with DEGseq [27]. DEGs at each stage or site were used for further analyses of GO (gene ontology) Molecular Function and KEGG pathways using by the Database for Annotation, Visualization, and Integrated Discovery (DAVID: http://david.abcc.ncifcrf.gov).

2.5. PPI Network Construction and Module Analysis

The PPI network was predicted using Search Tool for the Retrieval of Interacting Genes (STRING; http://string-db.org) [28] online database. Analyzing the functional interactions between proteins may provide insights into the mechanisms of generation or development of diseases. In the present study, PPI network of DEGs was constructed using STRING database, and an interaction with a combined score >0.7 was considered statistically significant. Cytoscape is an open source bioinformatics software platform for visualizing molecular interaction networks [29]. The plug-in ClusterONE (version 1.0) of Cytoscape is an APP for clustering a given network based on protein-protein interaction networks [30]. The PPI networks were drawn using Cytoscape and the most significant module in the PPI networks was identified using ClusterONE. Subsequently, the KEGG and GO analyses for genes in this module were performed using DAVID.

2.6. Hub Genes Selection and Analysis

The hub genes were selected with degrees ≥10 [31]. A network of the genes and their coexpression genes was analyzed using cBioPortal (http://www.cbioportal.org) online platform [4, 14]. The overall survival and disease-free survival analyses of hub genes were performed using Kaplan-Meier curve in cBioPortal [32].

2.7. Quantitative Real-Time PCR

A quantitative real-time PCR (qRT-PCR) assay was developed for the detection and quantification of AXIN2, WNT5B, WNT3A, CATENIN, and GSK3B transcripts using ACTIN as an endogenous control. The primers used in this study were designed as follows: ACTIN forward 5′-CACCCAGCACAATGAAGATCAAGAT-3′; reverse 5′-CCAGTTTTTAAATCCTGAGTCAAGC-3′. P2RY4 forward 5′-CTGGACTGTTGGTTTGATGAGGA-3′; reverse 5′-CAGCGACAGCACATACAAGGT-3′; C5 forward 5′-GGAGTGACGGTGCTGGAGTTT-3′; reverse 5′-CCCTCGTGCCAAAGTGGATAA-3′; SSTR5 forward5′- CTACATTCTCAACCTGGCAGTGG-3′; reverse 5′-GCTCATGACTGTCAGGCAGAAGA-3′. LPAR6 forward 5′-GGTAAGCGTTAACAGCTCCCAC-3′; reverse 5′- CATTTCGGACTTTGAGGACGC-3′. GPR68 forward 5′-CAACTCCTCGATGAGCTGTACCA-3′; reverse 5′- AGGTAGCCGAAGTAGAGGGACA-3′. qRT-PCR was performed in a 20 μL mixture containing 2 μL of the cDNA preparation, 10 μL 2x SYBR Green Premix ExTaq (Takara), and 10 μM primer on an ABI 7500 Real-Time PCR System (Applied Biosystems, Grand Island, NY, USA). The PCR conditions were as follows:10 min at 95°C, followed by 40 cycles of 15s at 95°C, and 1 min at 60°C. Each sample was tested in triplicate. Threshold values were determined for each sample/primer pair, and the average and standard errors were calculated.

3. Results

3.1. Emodin Inhibits HepG2 Cell Growth in the Time- and Dose-Dependent Manner

We analyzed the effect of emodin on HepG2 cell proliferation. CCK8 assays were performed using the HCC cell line, HepG2, after treatment with different concentrations of emodin for 24, 48, and 72 h. As shown in Figure 2, emodin decreased cell viability in a dose- and time-dependent manner, and the 50% inhibitory concentration (IC50) observed at 72 h was 19.12 μM in HepG2 cells.

3.2. RNA-Seq Landscape

Treatment of HepG2 cells via emodin dose-response assays revealed that emodin had IC50 values of about 20 μM at 72 h. The whole transcriptomic profiles of HepG2 cells and of those treated with emodin were assessed at base-pair resolution via RNA sequencing. After removing the low-quality reads, the raw reads from each sample were mapped to a reference genome (Table 1).

3.3. Identification of DEGs

DEGs were screened and generated by comparing the emodin samples to the control samples. As shown in the column diagram (Figure 3(a)) and volcano plot (Figure 3(b)), the analysis identified 859 DEGs comprising 712 downregulated genes and 147 upregulated genes that satisfy the criteria (fold change) >3 and adj. p-value <0.001.

3.4. KEGG and GO Enrichment Analyses of DEGs

To determine the biological classifications of the DEGs, functional and pathway enrichment analyses were performed using DAVID, and the GO functional enrichments of the upregulated and downregulated genes were determined at a p-value of less than 0.05. GO analysis results showed that changes in biological processes (BP) of DEGs were significantly enriched in cellular heat acclimation, positive regulation of apoptotic process, response to lipopolysaccharide, negative regulation of inclusion body assembly, negative regulation of extrinsic apoptotic signaling pathway in absence of ligand, negative regulation of endopeptidase activity, extracellular matrix organization, homophilic cell adhesion via plasma membrane adhesion molecules, wound healing, and fibrinolysis (Table 2). Changes in cell component (CC) of DEGs were mainly enriched in extracellular space, extracellular region, proteinaceous extracellular matrix, integral component of membrane, integral component of plasma membrane, and myofibril (Table 2). Changes in molecular function (MF) were mainly enriched in ATPase activity, coupled, serine-type endopeptidase activity, serine-type peptidase activity, bradykinin receptor activity, protein domain specific binding, calcium ion binding, serine-type endopeptidase inhibitor activity, extracellular matrix structural constituent, structural molecule activity, and lipopolysaccharide binding (Table 2). KEGG pathway analysis revealed that the DEGs were mainly enriched in complement and coagulation cascades, amoebiasis, ECM-receptor interaction, protein digestion and absorption, platelet activation, MAPK signaling pathway, rheumatoid arthritis, estrogen signaling pathway, legionellosis, and metabolism of xenobiotics by cytochrome P450 (Table 2).

3.5. PPI Network Construction and Module Analysis

To analyze the functional contributions of the DEGs, protein–protein interaction (PPI) network analysis was performed using STRING and Cytoscape. The false discovery rate (FDR) for each p-value was calculated. In general, the terms with FDR values < 0.01 were considered significantly enriched. The PPI network and the most significant module of the DEGs were determined (Figures 4(a) and 4(b)). Functional and pathway enrichment analyses of the genes involved in the module were performed using DAVID. Results showed that the genes in the most significant module were predominantly enriched for terms associated with inflammatory response, neuropeptide signaling pathway, and positive regulation of cytosolic calcium ion concentration (Table 3).

3.6. Hub Gene Analysis

A total of 25 genes were identified as hub genes with degree ≥10. The 25 most significant genes showing significant interaction were OPRD1, AVP, BDKRB2, TAS2R4, KNG1, BDKRB1, AGT, PTGDR2, LPAR6, C5, OPRL1, ADRA2C, CCL16, OXER1, CORT, SSTR5, PYY, MCHR1, UTS2R, ANXA1, ADCY1, GPR68, PIK3R1, P2RY4, and HCAR2. The functional roles of these hub genes are shown in Table 4. Coexpression of these hub genes was determined using cBioPortal (Figure 5). To further analyze these hub genes, overall survival and disease-free survival analyses were performed using a Kaplan-Meier curve. HCC patients with C5 and LPAR6 alterations showed low overall survival (Figure 6(a)), and HCC patients with SSTR5, P2RY4, LPAR6, and GPR68 alterations showed low disease-free survival (Figure 6(b)). Results based on qRT-PCR revealed that the mRNA levels of LPAR6, C5, SSTR5, GPR68, and P2RY4 were significantly downregulated in HepG2 cells (Figure 7). These findings indicated that C5, SSTR5, P2RY4, LPAR6, and GPR68 play important roles in the molecular mechanisms involved in emodin therapy for HCC.

4. Discussion

Transcription is an important biological process that determines the proteome of the cells. In this study, the RNA profile of HepG2 cells was determined via RNA-seq to explore the functions of DEGs in cancer cells treated with emodin. A total of 859 DEGs were screened, including 712 downregulated genes and 147 upregulated genes in HepG2 cells treated with emodin. These DEGs were found to be associated with emodin-mediated inhibition of liver cancer proliferation. Functional and pathway enrichment analyses of the 859 DEGs were performed using DAVID. GO analysis results showed that the DEGs were mainly enriched in regulation of apoptotic process and extracellular matrix organization. KEGG pathway analysis revealed that the DEGs were mainly enriched in mitogen-activated protein kinase (MAPK) signaling pathway. MAPK is a family of protein kinases comprising the p38 MAPKs, the c-jun N-terminal kinases (JNKs), and the extracellular regulated kinases (ERKs) [36]. MAPKs play crucial roles in regulating various cellular processes, such as cell proliferation. Our findings indicated that emodin-mediated inhibition of HCC cell proliferation is associated with the MAPK signaling pathway. The DEGs induced by emodin treatment were then carefully identified using a suite of sequence analysis software packages. Weighted gene coexpression network analysis has been successfully used to identify coexpression modules and intramodular hub genes based on DEG expression data. Analyzing the functional interactions between proteins can provide insights into the mechanisms underlying the effects of emodin in HCC. PPI networks were generated using Cytoscape and the most significant module in the PPI networks was identified using ClusterONE. The hub genes of the most significant module were then selected based on the cutoff degree ≥10. A total of 25 genes were identified as hub genes: OPRD1, AVP, BDKRB2, TAS2R4, KNG1, BDKRB1, AGT, PTGDR2, LPAR6, C5, OPRL1, ADRA2C, CCL16, OXER1, CORT, SSTR5, PYY, MCHR1, UTS2R, ANXA1, ADCY1, GPR68, PIK3R1, P2RY4, and HCAR2. The functional and pathway enrichment analyses of the 25 hub genes were performed using DAVID. GO analysis results showed that changes in BP of hub genes were significantly enriched in positive regulation of cytosolic calcium ion concentration, G-protein coupled receptor signaling pathway, adenylate cyclase-inhibiting G-protein coupled receptor signaling pathway, neuropeptide signaling pathway, and inflammatory response. Changes in CC of hub genes were mainly enriched in integral component of plasma membrane and plasma membrane. Changes in MF were mainly enriched in neuropeptide binding and G-protein coupled receptor activity. KEGG pathway analysis revealed that the hub genes were mainly enriched in neuroactive ligand-receptor interaction. To further analyze these hub genes, overall survival and disease-free survival analyses were performed using a Kaplan-Meier curve. Results showed that C5 and LPAR6 are associated with the overall survival of HCC patients. SSTR5, P2RY4, LPAR6, and GPR68 are related to the disease-free survival of HCC patients. These five genes are closely related to prognosis of HCC patients. Verification of the results by qRT-PCR showed that emodin can downregulate the mRNA levels of LPAR6, C5, SSTR5, GPR68, and P2RY4 in HepG2 cells (Figure 7). Therefore, C5, SSTR5, P2RY4, LPAR6, and GPR68 are likely to play important roles in the molecular mechanisms involved in emodin therapy for HCC. In summary, the five hub genes can reflect the molecular mechanisms of emodin therapy for HCC and could serve as targets for emodin therapy for HCC.

Complement activation is regulated to provide a permanent source of complement mediators that can maintain the inflammatory microenvironment that favors tumor growth [13, 37]. Complement component 5 (C5) is the fifth component of the complement and can be cleaved into C5a and C5b by C5-convertase. It plays a critical role in cell killing and inflammatory processes [38]. Many studies have reported that a variety of tumor cell lines can produce the complement activation product C5a. Furthermore, it has been reported that plasma C5a levels are elevated in liver cancer and other tissue-specific cancers [1, 6, 15]. He et al. showed that C5 levels were upregulated in AFP(-) HBV-related HCC and that C5 is potentially strongly associated with the progression of AFP(-) HBV-related HCC [17]. In addition, tumor inflammatory microenvironments were found to contain the complement-activating components C3, C4, C5, C1q, and MAC in many cancer models [39]. Somatostatin receptor type 5 (SSTR5) is a receptor that can lead to somatostatin-mediated inhibition of the release of hormones and secretory proteins [40]. A previous study reported that SSTR5 levels are upregulated in advanced-stage HCC [23]. SSTR5 can bind to somatostatin analogues, such as octreotide, which can help determine the antiproliferative efficacy of somatostatin analogues [2]. In addition, a positive correlation has been reported between SSTR5 expression and tumor size [8]. Lysophosphatidic acid receptor 6 (LPAR6) is a G protein-coupled receptor that can bind to lysophosphatidic acid [41]. One study reported that LPAR6 is essential for maintaining the tumorigenic properties of HCC cells; patient data and the experimental evidence supported the claim that LPAR6 promotes tumorigenicity and growth in HCC by activating the protooncogene Pim-3 [33]. Emodin can downregulate the expression of C5, SSTR5, and LPAR6. So far, no studies have examined the expression levels of P2Y purinoceptor 4 (P2RY4) and G-protein coupled receptor 68 (GPR68) in HCC. P2RY4, a G-protein coupled receptor, is responsive to uridine nucleotides [18] and plays an important role in transporting chloride in the epithelium of the jejunum [10]. The role of P2RY4 in apoptosis and cell proliferation, based on the effect of extracellular nucleotides, has been studied in HCT8 and Caco 2 cells [7]. In ovarian cancer, GPR68 is coupled to the PLC/Ca2+ pathway via the Gq/11 protein [42]. GPR68 is nearly undetectable in the healthy pancreas but is highly expressed in pancreatic ductal adenocarcinoma (PDAC) [34]. In human ovarian cancer cells, GPR68 has been shown to promote adhesion of cells to the extracellular matrix [35]. Among the five hub genes, C5, SSTR5, and LPAR6 were found to be involved in HCC (Table 5). These three genes could serve as the targets of emodin therapy for HCC (Table 5). However, few studies have investigated the relationships among P2RY4, GPR68, and liver cancer, and elucidating these relationships could be a direction for future research.

5. Conclusion

The present study attempted to identify DEGs that may be involved in the molecular mechanisms of emodin therapy for HCC. A total of 859 DEGs and 5 hub genes were identified and may be regarded as targets of emodin therapy for HCC. However, further studies are needed to elucidate the molecular mechanisms and biological function of these genes involved in emodin therapy for HCC.

Abbreviations

HCC:Hepatocellular carcinoma
RNA-seq:RNA sequencing
DEGs:Differentially expressed genes
BP:Biological processes
CC:Cell component
MF:Molecular function
GO:Gene ontology
KEGG:Kyoto Encyclopedia of Genes and Genomes
DMSO:Dimethylsulfoxide
qRT-PCR:Quantitative real-time PCR
IC50:Inhibitory concentrations of 50%
MAPK:Mitogen-activated protein kinase
JNKs:c-jun N-terminal kinases
ERKs:Extracellular regulated kinases
C5:Component 5
SSTR5:Somatostatin receptor type 5
LPAR6:Lysophosphatidic acid receptor 6
P2RY4:P2Y purinoceptor 4
GPR68:Ovarian cancer G-protein coupled receptor 1
PDAC:Pancreatic ductal adenocarcinoma.

Data Availability

The data used to support the findings of this study are included within the article.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Authors’ Contributions

Rui-sheng Zhou, Qin-feng Sun, and Zeng Jie Ye conducted experiments, collected the data, and wrote the manuscript. Jian-wei Liu and Dai-Han Zhou conducted experiments. Ying Tang collected data, contributed to the discussion, and reviewed the manuscript. Xiong-Wen Wang revised the manuscript. Ying Tang is the guarantor of this work and, as such, has full access to all the data in the study and takes responsibility for the integrity of the data and the accuracy of the data analysis. Rui-sheng Zhou and Xiong-Wen Wang contributed equally to this work.

Acknowledgments

This work was supported by High-level University Construction of Guangzhou University of Chinese Medicine (A1-AFD018181A29 and A1-2606-18-414-016).