Transcriptome Analyses of the Anti-Proliferative Effects of 20(S)-Ginsenoside Rh2 on HepG2 Cells

20(S)-ginsenoside Rh2 (Rh2), a well-known protopanaxadiol-type ginsenoside from Panax ginseng has especially gained attention for its anticancer activities on various types of human cancer cells. However, the molecular mechanism through which Rh2 promotes apoptosis in hepatocellular carcinoma (HePG2) cells is not known at the transcriptome level. Rh2 can specifically inhibit the proliferation of HePG2 in a dose- and time-dependent manner. Moreover, Rh2 can significantly increase the apoptosis which was related with an increase in protein expression levels of caspase-3, caspase-6, and poly (ADP-ribose) polymerase. Comparison of RNA-seq transcriptome profiles from control group and Rh2-treated group yielded a list of 2116 genes whose expression was significantly affected, which includes 971 up-regulated genes and 1145 down-regulated genes. The differentially expressed genes in p53 signaling pathway and DNA replication may have closely relationships to the cells apoptosis caused by Rh2 treatment. The results of qPCR validation showed that dynamic changes in mRNA, such as CDKN1A, CCND2, PMAIP1, GTSE1, and TP73.


INTRODUCTION
Ginseng (Panax ginseng C.A.Mey) has been known as the "King of Herbs" since ancient times in eastern Asia (Hossen et al., 2017). It is a perennial herb that grows in cool habitats and belongs to the Araliaceae family (Park and Cho, 2009). The root of ginseng been widely used in traditional medicine in eastern Asian countries to promote health (Yu et al., 2015). Ginsenosides are the main essential bioactive ingredients that have various pharmacological activities, including antiinflammatory, anti-allergic, anti-fatigue, anti-stress, and anti-cancer activities (Oh et al., 2014). More than 100 ginsenosides from Panax species have been identified and categorized into four types: protopanaxadiol, protopanaxatriol, oleanolic acid, and ocotillol-type (Nag et al., 2012). Ginsenoside Rh2 is a metabolite produced by Rg3, Rb1, Rb2, and Rc by human intestinal bacteria after absorption (Chi et al., 2005). Ginsenoside Rh2 is classified into 20(S) and 20(R) forms according to the different arrangement of the hydroxyl group at the C-20 position, and 20(S)-Rh2 was found to possess strong cytoxic activities to several cancer cells (Liu et al., 2010;Kim et al., 2017). 20(S)-Rh2 has a number of pharmacological properties, such as anti-inflammatory, anticancer, cardioprotection, and neuroprotective effects . Rh2 has especially gained attention for its broadspectrum anti-proliferative effect on a variety of human cancer cells (Oh et al., 1999;Kim et al., 2004;Qi et al., 2011;Guo et al., 2014;Shi et al., 2014;Qian et al., 2016;Yang et al., 2016;Wang et al., 2018). In addition, additional in vivo work showed that Rh2 also can significantly suppress the growth of uterine leiomyoma in SD rats (Zhu et al., 2015). The tumor growth inhibitory activity of Rh2 was also found in H22 tumor-bearing mice (Chen et al., 2017) and nude mice bearing human ovarian cancer cells (Nakata et al., 1998).
The molecular pathways affected by Rh2 in cancer cells are mainly through apoptosis-related pathways. for example, Rh2 induces apoptosis in colorectal cancer cells through activation of p53 (Li et al., 2011). in human leukemia Reh cells, Rh2 could induce apoptosis through the mitochondrial pathway (Xia et al., 2014). Rh2 also suppresses proliferation of hepatocellular carcinoma (HCC) cells by targeting EZH2 to regulate CDKN2A-2B gene cluster transcription (Li et al., 2017). Although apoptosis and other related pathways and genes have been identified in Rh2-treated cancer cells, the molecular mechanism through which Rh2 promotes cancer cell apoptosis is not known at the transcriptome level. To the best of our knowledge, this study is the first transcriptome analysis of the anti-proliferative effects of Rh2 on hepg2 cells. The results will aid the identification of functional genes and provide a more comprehensive understanding of the role of Rh2 in promoting apoptosis in human hematoma cells.

Cell line and Cell Culture
The HepG2 cell line was purchased from the Cell Bank of Type Culture Collection of the Chinese Academy of Sciences (Wuhan, China) and cultured in minimum essential medium supplemented with 10% fetal bovine serum, 100 U/ml penicillin, and 100 μg/ml streptomycin at 37°C in a humidified atmosphere of 95% air-5% CO 2 in an humidified CO 2 incubator (Heracell 150i, Thermo Fisher Scientific, Waltham, MA, USA).

Cell Viability
The cytotoxicity of Rh2 on HepG2 cells was determined using the MTT assay as described previously (Bai et al., 2019). Briefly, cells were seeded (1 × 10 4 cells/well) in 96-well plate for overnight. Then the cells were incubated with an increasing concentration of Rh2 (0-60 μM) for 24 and 48 h, respectively. The supernatant was removed via aspiration and 100 μl of MTT solution (0.5 mg/ ml) was added and cells were continuously cultured until 100 μl of MTT stop solution (10% SDS, 0.01 M hydrochloric acid) was added to each well to solubilize the formazan. The absorbance was measured at 550 nm using a microplate reader (Tecan Infinite M200 Pro, Männedorf, Switzerland).

Cell Morphology Analysis
To examine the effect of Rh2 on the morphology of HepG2 cells, cells (5 × 10 4 cells/well) were seeded on sterilized cover glasses in 6-well plate for overnight. Then the cells were incubated with Rh2 (20 or 50 μM) for 24 h and the cell morphology was observed with Olympic (Tokyo, Japan). The changes of microstructure were surveyed under scanning electron microscope (SEM). Treated cells were washed with PBS three times and fixed with 2.5% glutaraldehyde solution for 1 h at room temperature. Subsequently, cells were rinsed again with PBS three times, and then dehydrated using graded ethanol (30%, 50%, 70%, 90%, 100%). Cells were fixed onto a copper stub and coated with gold using a sputter coater under vacuum, and then the surface morphology of cells was characterized by SEM (FEI Quanta 450 FEG, Hillsboro, USA).

Annexin V-Fitc Apoptosis Assay
The cells for apoptosis were examined using Annexin V-FITC apoptosis detection kit. Cells (4 × 10 5 cells/well) were seeded in 6-well plate for overnight. Then the cells were incubated with Rh2 (20 or 50 μM) for 24 h. After treatment, the cells were collected by trypsinization with 0.25% EDTA and centrifuged at 1500 rpm for 3 min to remove the medium. Then the precipitation was washed with pre-cold PBS and resuspended in 100 μl 1 × binding buffer containing 5 μl FITC-Annexin V and 5 μl propidium iodide. After incubation in the dark for 15 min, the fluorescence intensity was evaluated under an Accuri C6 Plus flow cytometer (Becton, Dickinson and company, CA, USA).

RNA Extraction, library Construction, and Sequencing
HepG2 cells were seeded in 40 mm dishes at a density of 1 × 10 6 cells/well. Following treatment with 50 µM Rh2 for 24 h, the cells were washed twice with cold PBS and collected by centrifuged for 3 min. Total RNA of each sample was extracted with a Trizol reagent kit (Invitrogen). RNA quality was assessed on an Agilent 2100 Bioanalyzer (Agilent Technologies, Palo Alto, CA, USA) and checked by agarose gel electrophoresis. After enrichment of mRNA using Oligo(dT) beads, it was fragmented into short fragments in fragmentation buffer and reverse transcripted into cDNA with random primers. Then, DNA polymerase I, RNase H, dNTP, and buffer were used to synthesize the second-strand cDNA. The synthesized product were purified with a PCR purification kit (Qiagen, Venlo, The Netherlands), end repaired, poly(A) added, and ligated to Illumina sequencing adapters. After that, the ligated products were size selected by agarose gel electrophoresis, PCR amplified, and sequenced with Illumina HiSeq2500 (Gene Denovo Biotechnology Co. Guangzhou, China).

RNA-Seq Data Analysis
In order to ensure the quality of the data for the following analysis, the raw data of six samples were firstly processed through in-house perl scripts. In this step, clean data were obtained by removing reads containing adapter, reads containing ploy-N, and low quality reads from raw data. At the same time, Q20, Q30, and GC content of the clean data were calculated. All the downstream analyses were based on the clean data with high quality.

Differentially Expressed Genes
RNAs differential expression between control and Rh2 treated groups was analyzed using DESeq2 software (Love et al., 2014). We defined the genes with false discovery rate (FDR) < 0.05 and absolute fold change (FC) ≥ 2 as differentially expressed genes (DEGs).

Gene Ontology and Pathway Enrichment Analysis
Gene Ontology (GO) is a widely used ontology in the bioinformatics field, which offers a semantic vocabulary standard to define and describe the functions of genes and proteins, and can be updated and applied in any organism (Myhre et al., 2006). Hierarchical system was used in GO for gene classification and putting identical genes at the same level. The vocabulary of genes and proteins involved in GO covers three aspects of biology: molecular function (MF), cellular component (CC), and biological process (BP). Terms is the basic units of GO, which were also known as GO classes. Each GO term has a human readable name and a GO ID, and belongs to one of the three sub-ontologies.
Genes often interact with other genes when performing their functions in BP. We often classify genes involved in same biological functions into pathways, which are helpful for us to further understand the biological functions of genes. Kyoto Encyclopedia of Genes and Genomes (KEGG) is a pathwayrelated database resource for understanding high-level functions and utilities of the biological systems (Lin and Feng, 2012). KEGG pathway enrichment analysis in DEGs comparing with the whole genome background can identify DEGs into significantly enriched metabolic pathways or signaling pathways.
In the present study, we use the online DAVID functional annotation clustering tool (https://david.ncifcrf.gov) by choosing the default option to determine the significantly represented GO categories and KEGG enriched pathways for DEGs of HepG2 cells response to Rh2 treatment. The calculated p value was corrected with FDR ≤0.05 as a threshold under the FDR correction. At the same time, GO terms and KEGG Pathways of DEGs meeting this condition were defined as significantly enriched GO terms and KEGG pathways.

Real Time Quantitative RT-PCR
HepG2 cells were seeded in 40 mm dishes at a density of 1 × 10 6 cells/well. Following treatment with 50 µM Rh2 for 24 h, the RNA extraction was carried out as described in the section of RNA extraction, library construction, and sequencing. A 20 µl PCR reaction system was consist with 2 µl cDNA, 10 µl SYBR mixture, 1 µl forward primer, 1 µl reverse primer, and 6 µl deionized water. After mixing, the PCR reaction was performed using CFX-96™ Real-Time instrument (Bio-Rad, CA, USA). The GAPDH gene was used as a house keeping gene to normalize the expression level of the test genes, and the relative gene expression level was analyzed using the 2 −ΔΔCT method. All of the samples were analyzed in triplicate. Primers were synthesized by Sangon Biotech (Shanghai, China) and were listed in Table 1.

Western Blot Analysis
HepG2 cells were seeded in 40 mm dishes at a density of 1 × 10 6 cells/well. Following treatment with 20 or 50 µM Rh2 for 24 h, the cells were washed twice with cold PBS and centrifuged for 3 min, and the supernatant was removed. The cells were lysed in 200 µl of lysis buffer and protein content was measured with Bradford reagent using bovine serum albumin as a standard. Equal amounts of protein were resolved on 10%-15% SDSpolyacrylamide gels for 2 h at 120 V and then transferred to PVDF membranes. The membranes were blocked for 2 h using 5% milk powder in Tris-buffered saline containing Tween 20 (1 × TTBS). After incubation, the membranes were washed

Statistical Analysis
The results have been represented as the means standard deviation (SD). The significance was analyzed using SPSS 20.0 package (SPSS Inc., Chicago, IL, USA). Differences among samples were compared using a Duncan's multiple range test and p values less than 5% were considered to be statistically different.

Rh2 Inhibits HepG2 Cell Proliferation
HCC is the most malignant primary human liver cancer and has a poor outcome after combined surgical treatment, radiotherapy, and chemotherapy (Lin, 2014). Although our understanding of the mechanism of HCC carcinogenesis has improved, the prognosis of HCC remains poor. Hence, there is an urgent need to develop effective therapies for HCC (Shang et al., 2016). There is growing interest in developing more effective therapeutic agents from natural resources.
To confirm the effects of Rh2 on HepG2 human hepatocarcinoma cells, we performed an MTT assay (Figure 1) and found that Rh2 inhibited HepG2 cell growth in a concentration-dependent manner with an IC 50 of 51.97 and 45.46 μM at 24 and 48 h, respectively, which is consistent with the literature (Chen et al., 2015;Shi et al., 2016). To confirm the cytotoxic effects of Rh2 on HepG2, cell morphology was visualized under an inverted microscope and by SEM. As shown in Figure 2A, the untreated cells were spread out and flattened, whereas the Rh2-treated cells were round, shrunken, and showed altered adherence. Moreover, cell numbers were reduced. As shown in Figure 2B, control cells had a normal cell architecture  with a distinct cytoskeleton, whereas the Rh2-treated cells had atypical morphological changes, including cell shrinkage and incomplete membranes.

Rh2 Induces Apoptosis in HepG2 Cells
Apoptosis, also known as programmed cell death, is an important regulatory mechanism by which the body maintains cell structure and normal tissue development (Bennetts and Pierce, 2010). Inducing tumor cell apoptosis is an important way to inhibit tumor growth. Rh2 inhibits HCC cells in many ways, including inhibiting cell proliferation, inducing cell differentiation, promoting cell apoptosis, inhibiting cell invasion and metastasis, reducing drug resistance, and improving immunity (Li et al., 2017). Although the anti-HCC cell activity of Rh2 has been extensively reviewed, only a few studies have investigated the molecular responses in HCC HepG2 cells (Li et al., 2011;Cheong et al., 2015;Chen et al., 2016;Shi et al., 2016). In this study, we investigated the anti-cancer effects of Rh2 on HepG2 cells and the molecular mechanism by which Rh2 promotes HepG2 cell apoptosis at the transcriptome level.
To confirm the increased apoptotic effects of Rh2 on HepG2 cells, the cell lines were treated with 20 and 50 μM Rh2 for 24 h. After treatment, the cells were examined for apoptosis via flow cytometry using Annexin V-FITC/PI double staining ( Figure  3A). The percentage apoptosis increased from 0.057 in the control group to 46.5% and 66.07% (Figure 3B) at 20 and 50 μM, respectively. These results indicate that the anti-proliferative effects of Rh2 on HepG2 cells are mostly cytostatic, rather than cytotoxic. Caspases are a family of cysteine proteases and the primary drivers of apoptotic cell death, and play a critical role in the transduction of apoptosis signals (Vaculova and Zhivotovsky, 2008). Among the major caspases, caspases 3 and 6 are executioners of tumor growth and progression. PARP also plays important roles in apoptosis, the DNA damage response, and cell cycle regulation (Smulson et al., 2000). Interestingly, apoptosis is a genetically regulated cell suicide program mediated by activation of effector caspases 3 and 6. The effects of Rh2 on the expressions of apoptosis-related proteins were determined by assessing the expression levels of caspases 3 and 6 and PARP protein in HepG2 cells via Western blotting. As shown in Figure 4, treatment with 50 μM Rh2 dramatically increased caspase 3 and 6 expression in HepG2 cells but had less of an effect on cleaved-PARP protein. These results indicate that Rh2 promotes the apoptosis of hepatoma cells by upregulating the expression of caspase family proteins.

Data Generation and Comparison With RNA-Seq Analyses
To identify what genes were regulated via Rh2 treatment, and further elucidate the molecular mechanism by which Rh2 inhibits HepG2 cell proliferation, transcriptome analyses were carried out on Rh2-treated and control HepG2 cell samples. Illumina HiSeq2500 was used for sequencing. As shown in Table 2, there were 70,625,504 (99.77%), 60,206,338 (99.78%), and 66,661,250 (99.79%) clean reads in control groups CK-1 to -3 and 60,916,030 (99.81%), 62,125,098 (99.78%), and 67,314,974 (99.76%) in Rh2-treated groups 1-3, respectively. The Q30 proportion of all six samples exceeded 93%, indicating that the quality of the clean reads obtained was high, and the results met the requirement for the next step in bioinformatics analyses. The rates of clean reads for all six samples exceeded 97%, indicating that the utilization rate of sequencing data was high.

DEG Response to Rh2 Treatment
To identify the DEG response to Rh2 treatment, the differential expression multiple between different samples was determined from the level of gene expression. Using FC ≥ 2 and FDR < 0.05  as criteria, 2116 DEGs were identified between the control and Rh2-treated groups. Next, the DEG clusters were analyzed in a heat map ( Figure 5A). DEGs identified in biological replicates clustered together, indicating good reproducibility of the treatments. The DEG statistics indicated 971 upregulated genes and 1145 downregulated genes ( Figure 5B).

GO Enrichment Analyses of DEGs
The DEGs in the CK vs. Rh2 groups were analyzed for GO functional enrichment. Each group of DEGs was annotated into three classifications: MF, CC, and BP. The top 25 functional enriched classes of each group are shown in Figure 6. In MF, the DEGs were mainly enriched for protein dimerization activity, electron carrier activity, heme binding, and heparin binding ( Figure 6A). In CC, the DEGs were mainly enriched in the extracellular region part, microsome, endoplasmic reticulum, and vesicular fraction ( Figure 6B). In BP, the DEGs were mainly enriched in oxidation reduction, response to wounding, wound healing, and response to cAMP ( Figure 6C).

KEGG Pathway Enrichment Analyses of DEGs
KEGG pathway analyses of the DEGs in CK vs. Rh2 were performed online at OmicShare. Pathways and all of the DEGs clustered into six categories: organismal systems, metabolism, human diseases, environmental information processing, cellular process, and genetic information processing ( Figure 7A). The environmental information processing category includes signal transduction, signaling molecules and interaction, and membrane transport. The cellular processing category includes transport and catabolism, cell growth and death, cellular community-eukaryotes, and cell motility. A scatterplot was used to show the top signaling pathways in the KEGG enrichment analyses. As shown in Figure 7B, the apoptosis-related signaling pathways, including the p53 signaling pathway and DNA replication, were significantly enriched. The DEGs in these two pathway categories may have close relationships to the cell apoptosis caused by Rh2 treatment, highlighting the modulation of Rh2 in tumorigenesis

DEGs in Apoptosis and Cell Cycle-Related Pathways
Transcription factor p53 is regulated by its phosphorylation (Zhang et al., 2002). p53 is the upstream regulator of caspase family proteins in oxidative stress and DNA damage induced by various stress conditions (Vakifahmetoglu et al., 2006). The function of p53 is to induce mitochondria to release cytochrome C, which directly improves caspase 9 activity (Schuler et al., 2000). The activation of caspase 3 is promoted by the activation    of caspase 9 (Kaushal et al., 2004). The p53-mediated apoptotic signaling pathway is important for cancer treatment with multiple chemotherapies (Zhan et al., 2014). Recent studies have shown that various natural products can induce apoptosis, which is closely associated with the p53 pathway (Laux et al., 2004;Yeruva et al., 2005;Volate et al., 2010). Table 3 lists the DEGs in KEGG pathways related to cell apoptosis and the cell cycle. Fifteen DEGs were annotated in the p53 signaling pathway: six were upregulated and nine were downregulated.

Differential Expression of Representative DEGs
Based on the enrichment results of GO and KEGG, 15 DEGs belonging to the p53 signaling pathway were selected as representative DEGs for a heat map and are listed in Table 4. The color represents the relative expression of the genes. Red represents upregulation and green represents downregulation. As shown in Figure 8 and

RT-qPCR Verification of DEGs
RT-qPCR analyses were used to verify the accuracy of the transcriptome sequencing data. Eleven DEGs in the p53 signaling pathway were selected for qPCR analyses, including PMAIP1, SESN1, TP73, GTSE1, CCNE1, TP53I3, CDKN1A, CCND3, CCND2, BBC3, and DDB2. As shown in Figure 9, the DEGs PMAIP1, CDKN1A, CCND2, and BBC3 were significantly upregulated, while SESN1, TP73, GTSE1, CCNE1, TP53I3, and DDB2 were significantly downregulated on qPCR analyses. The expression of CCND3 was slightly downregulated. This differs slightly from the selection criteria (FC ≥ 2, FDR< 0.05) for DEGs. In general, the relative expression of the selected DEGs showed similar tendencies to the RNA-Seq results (Tables 3 and 4), indicating that the sequencing results were reliable.
FIGURE 9 | RT-qPCR verification of the selected DEGs. The mRNA expression levels of the test DEGs in the CK vs. Rh2 groups were quantified using RT-qPCR and the housekeeping gene GAPDH was used to normalize the relative expression level. Within the same gene, *p < 0.05 vs. control group.

CONClUSIONS
We demonstrated that Rh2 can inhibit proliferation and induce apoptosis in HepG2 cells. Investigations of the possible mechanism of the anti-proliferation effect revealed that Rh2 induced marked apoptosis that was related to changes in mRNAs in the p53 signaling pathway. Additional studies are currently underway to investigate its overall in vivo anticancer effect and potential mechanisms.

DATA AVAIlABIlITY STATEMENT
The data generated for this study can be found in the Sequence Read Archive, accession number PRJNA548682.