Irinotecan (CPT11) is one of the most effective drugs for treating colon cancer, but its severe side effects limit its application. Recently, a traditional Chinese herbal preparation, named PHY906, has been proved to be effective for improving therapeutic effect and reducing side effects of CPT11. The aim of the present study was to provide novel insight to understand the molecular mechanism underlying PHY906-CPT11 intervention of colon cancer. Based on the GSE25192 dataset, for different three treatments (PHY906, CPT11, and PHY906-CPT11), we screened out differentially expressed genes (DEGs) and constructed a co-expression network by weighted gene co-expression network analysis (WGCNA) to identify hub genes. The key genes of the three treatments were obtained by merging the DEGs and hub genes. For the PHY906-CPT11 treatment, a total of 18 key genes including Eif4e, Prr15, Anxa2, Ddx5, Tardbp, Skint5, Prss12 and Hnrnpa3, were identified. The results of functional enrichment analysis indicated that the key genes associated with PHY906-CPT11 treatment were mainly enriched in ‘superoxide anion generation’ and ‘complement and coagulation cascades’. Finally, we validated the key genes by Gene Expression Profiling Interactive Analysis (GEPIA) and RT-PCR analysis, the results indicated that EIF4E, PRR15, ANXA2, HNRNPA3, NCF1, C3AR1, PFDN2, RGS10, GNG11, and TMSB4X might play an important role in the treatment of colon cancer with PHY906-CPT11. In conclusion, a total of 18 key genes were identified in the present study. These genes showed strong correlation with PHY906-CPT11 treatment in colon cancer, which may help elucidate the underlying molecular mechanism of PHY906-CPT11 treatment in colon cancer.
Colon cancer is a common malignant tumor and the second leading cause of mortality worldwide in 2018 . With the improvement in people’s living standards and changes in diet and other aspects, the incidence of colon cancer has increased year by year. Colon cancer development is related to many factors, such as age, sex, race/ethnicity, genetic makeup, dietary factors, increased BMI, red meat intake, low physical activity, long-term cigarette smoking (30–40 years), low vegetables and fruit consumption . Studies showed that intestinal flora plays an important role in the occurrence and development of colon cancer . Intestinal flora imbalance will cause damage to the intestinal barrier and impaired innate immune function, leading to inflammation and infection. There is mounting evidence that the intestinal flora can be used as biomarkers for cancer diagnosis and prognostication [4–6], which indicated that the intestinal flora will become an important component of future cancer prevention and therapy. In recent years, several studies have shown that colon cancer is closely related to inflammatory bowel disease (IBD) , and the relative risk of colon cancers in IBD patients is approximately two- to three-times greater than that in the general population. Additionally, a number of studies have reported the importance of single inflammatory genes in colon cancer [8,9]. For instance, a previous study conducted by Jiang and colleagues identified 14 inflammatory genes with high prognostic significance, and found a novel role of the inflammation–associated network in colon cancer . Based on these, anti-inflammatory drugs have become a new focus of researchers seeking effective drugs for colon cancer.
The therapy of colon cancer includes surgical resection, radiotherapy, chemotherapy, and others. Chemotherapy is an important treatment for colon cancer in all stages. Irinotecan (CPT11), a DNA topoisomerase I inhibitor, does not hinder the binding of topoisomerase I, but converts this ribozyme into a substance that is harmful to DNA. CPT11 is considered as one of the most effective drugs for treating colon cancer . However, the clinical application of CPT11 is often accompanied by a variety of side effects, including neutropenia, bone marrow suppression, and delayed diarrhea [12–14].
Traditional Chinese medicine (TCM) is one of the unique cultural treasures of China, which has gained growing recognition around the world in recent years. PHY906, one of the TCM, is composed of four main herbs, Scutelleria baicalensis Georgi, Paeonia lactiflora Pall., Glycyrrhiza uralensis Fisch., and Ziziphus jujuba Mill. PHY906 has been used for over 1800 years for treating gastrointestinal symptoms. Zou and colleagues  found that PHY906 exerts a protective effect in IBD by down-regulated production of pro-inflammatory cytokines and suppression of the NF-κB signaling pathway. Several clinical studies have demonstrated that TCM, as an alternative treatment for cancer, can reduce the incidence of bone marrow suppression and gastrointestinal reactions in chemotherapy [16–18]. In addition to mitigate the side effects of chemotherapy, the combination of chemotherapy and TCM also improves the chemosensitivity of diseases. Several Phase I/II clinical trials with encouraging results have been conducted in the United States to explore the toxicity and clinical efficacy of PHY906. PHY906 was shown to provide cytoprotective effects without dampening the anti-tumor activity of CPT11 in colon cancer patients under chemotherapy. A study demonstrated that PHY906 significantly amplifies the effects of CPT-11 in the tumor tissue; combination of PHY906 and CPT11 (PHY906-CPT11) enhanced anti-tumor activity by creating a unique tissue-specific response . However, few studies focused on the mechanism of anti-tumor modulation of the PHY906-CPT11 in chemotherapy.
Currently, microarray technology has been widely applied to elucidate the tumor progression in different conditions. Weighted gene co-expression network analysis (WGCNA) is a new tool that is used to identify co-expressed modules and hub genes. It has been reported that WGCNA can be useful to identify candidate biomarkers or therapeutic targets of various diseases, such as Alzheimer’s Disease , psychiatric disorders , glioblastoma , and colon cancer . At present, it has been widely used to compare differentially expressed genes (DEGs) and help investigating the interactions among genes in different modules . Besides, it is also used to identify rules for predicting survival of patients by investigating the relationship between clinical traits and tissue microarray data .
In the present study, we aimed to investigate molecular mechanism of different treatments in colon cancer. We analyzed the microarray datasets of colon cancer treated with PHY906, CPT11, PHY906-CPT11, and Phosphate Buffered Saline (PBS), respectively. DEGs from three different treatments were identified, then WGCNA was constructed to screen hub genes associated with the three treatments. DEGs and hub genes associated with different treatments were intersected and key genes which were significantly associated with the treatment of colon cancer were obtained finally. The functional annotation of genes (DEGs, hub genes, and key genes) were performed to understand the functions of these genes.
Materials and methods
The gene expression profiles of GSE25192 submitted by Wang et al. (2011)  were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo). The dataset comprised four groups: ten tumor-bearing mice treated with PBS (‘PBS’), nine tumor-bearing mice treated with PHY906 (‘PHY906’), ten tumor-bearing mice treated with CPT11 (‘CPT11’), and nine tumor-bearing mice treated with the combination of PHY-906 and CPT-11 (‘PHY906-CPT11’).
Screening and functional enrichment analysis of DEGs
The raw file data were preprocessed in ‘limma’ R package , including background correction and expression normalization. Under the threshold of FDR < 0.05 and |log2FC| ≥ 1, the ‘limma’ R package was used to screen out the DEGs between PHY906 and PBS, CPT11 and PBS, PHY906-CPT11 and PBS. Then the ‘pheatmap’ R package was used to construct the expression heatmap while the ‘ggpubr’ and ‘ggthemes’ R packages were used to construct the volcano plot of acquired DEGs. Venn diagram was performed by using the online tool jvenn (http://jvenn.toulouse.inra.fr/app/example.html) to overlap the DEGs in PHY906, CPT11, and PHY906-CPT11.
Construction of WGCNA and key module identification
The data of GSE25192 were used for cluster analysis to detect outliers using standardized connectivity (Z.K) method suggested by WGCNA authors , with the threshold Z.K score < −2. Then, the expression data were used as input datasets of WGCNA, and construct the sample dendrogram and trait heatmap after one outlier sample were excluded, and the similar gene expression profiles were divided into different gene modules. To ensure that the network can obey the scale-free criteria without affecting the network connectivity, the present study selected β = 8 as the soft-thresholding.
After determining the soft-thresholding, the network was developed. The adjacency matrix was transformed into a topological overlapping matrix to construct a network, and the gene dendrogram and module colors were established by using the degree of dissimilarity. To further analyze the module, we calculated the dissimilarity of the module eigengenes (MEs), hierarchically clustered the modules, and merged similar modules. In order to evaluate the co-expression similarity among different modules, characteristic adjacency degrees of correlation among each module were calculated. Interaction of different co-expression modules was assessed by flashClust function, and the heatmap was constructed.
After the modules were identified, the ME was summarized by the first principal component of the module expression levels. Module–trait relationships were estimated using the correlation between MEs and clinical treatment, which allowed efficient identification of the relevant modules. To evaluate the correlation strength, we calculated the module significance (MS), which is defined as the average absolute gene significance (GS) of all the genes in the module. The GS was measured as the log10 transformation of the P-value (lgP) in the linear regression between gene expression and clinical information. In the WGCNA, modules with the highest MS score among all modules are usually defined as the key module and selected for further analysis.
Transcription factors prediction and identification of hub genes for the key module
Enrichr (http://amp.pharm.mssm.edu/Enrichr/) is a comprehensive web-based tool that contains 180184 annotated gene sets from 102 gene set libraries. The gene information for the key module was imported into the Enrichr to obtain the interaction between transcription factors (TFs) and their target genes. To reduce the chance of identifying false-positives, we extracted only TFs with consensus target ChEA and ENCODE geneset libraries and with P<0.05 as determined by the Fisher exact test. We then used the Cytoscape 3.4.0 software (Cytoscape Consortium, SanDiego, CA, U.S.A.) to visualize the TF-target gene regulatory networks.
Hub genes were screened out using module connectivity, measured by MM ≥ median and clinical trait relationships, measured by GS ≥ median. Then we exported the full weighted network of the key module. Finally, all the obtained genes’ network from the full weighted network of the key module (if there are too many edges, we filtered the edges according to the weight to obtain the suitable number of edges for visualization by cytoscape) were extracted, and a subnetwork of the full weight network of the key module was finally obtained. Furthermore, the subnetworks were then constructed with the plug-in MCODE in the Cytoscape (degree cut-off ≥ 2, node score cut-off ≥ 0.2, K-core ≥ 2, and max depth = 100). Genes identified in the MCODE subnetworks and exhibiting high weight in the co-expression network were chosen as the hub genes for further analysis. The significance of the gene expression in different treatments was examined by the Wilcoxon’s test (P<0.05 was considered as significant differences).
Identification of key genes by merging DEGs with hub genes
The Venn map was constructed by website (http://bioinformatics.psb.ugent.be/webtools/Venn/) for intersection of the obtained data.
Validation of the key genes
An online tool, Gene Expression Profiling Interactive Analysis (GEPIA, http://gepia.cancer-pku.cn/), was used to analyze differential expression of the key genes and visualize their expression levels in different stages based on TCGA datasets. The overall survival analysis for each key gene was also performed using GEPIA, and the log-rank tests were used to measure the statistically significance (P<0.05).
Cell culture and treatment
Human colon cancer cell line SW480 was seeded in six-well plates (Corning, U.S.A.) at a density of 5 × 104 cells/well. The cells were cultured in Roswell Park Memorial Institute (RPMI)-1640 medium and incubated at 37°C with 5% CO2. A total of 10% fetal bovine serum (FBS), 100 U/ml penicillin, and 100 μg/ml streptomycin were contained in the culture media.
The cells were cultured to approximately 85% density and divided into four groups. Four groups of SW480 cells were respectively treated with PBS, PHY906 (500 μg/ml), CPT11 (5 μg/ml), and PHY906-CPT11 (500 μg/ml-5 μg/ml) for 24 h.
RNA isolation and qRT-PCR
Total RNA was extracted from SW480 that was treated with different treatments respectively by using the TRIzol reagent. The total RNA was used for reverse-transcription into cDNA using a specific kit for reverse transcription (Promega) according to the manufacturer-provided protocol. The qRT-PCR approach was used for amplification and determination of the mRNA expression level with the SYBR Green PCR Master Mix. The mRNA expression levels relative to the controls (GAPDH expression) were obtained using the 2−ΔΔCT method . Primer sequences used in the present study were provided in Table 1.
|Gene symbol .||Primer forward (5′–3′) .||Primer reverse (5′–3′) .|
|Gene symbol .||Primer forward (5′–3′) .||Primer reverse (5′–3′) .|
GO function and KEGG pathway enrichment analyses of the DEGs, key module, hub genes and key genes
To understand the biological meaning of the DEGs, key module, hub genes and key genes, clusterProfiler (R package) was used for performing GO enrichment and KEGG pathway analyses. P-value <0.05 was regarded as significant.
Screening and GO function and KEGG pathway enrichment analyses of DEGs
The volcano plots for all genes and the heatmap of the genes with the most significant differential expression were shown in Figure 1. In total, three DEGs were screened out (three up-regulated genes Exosc1, Slitrk2, and Afmid) between PBS and PHY906 (Figure 1A,B), 131 DEGs were screened out (74 up-regulated and 61 down-regulated) between PBS and CPT11 (Figure 1D,F). Genes such as HBB-B1, MYH4, and HBA-A2 were the most down-regulated while Sepp1, Afmid, and Calcb were the most up-regulated in PBS vs. CPT11 comparison. Moreover, 268 DEGs were identified (186 up-regulated and 82 down-regulated) between PBS and PHY906-CPT11, among which Capg, Metrnl, and Calcb were the most up-regulated while Hbb-B1, Hba-A2, and Prss2 were the most down-regulated (Figure 1H,I). Additionally, we found there was only one common gene (Slitrk2) among the DEGs from PHY906, CPT11, and PHY906-CPT11 (Supplementary Figure S1). As shown in Figure 1C, GO analysis revealed that the DEGs between PBS and PHY906 were enriched in different biological processes (BPs) among which tryptophan metabolic process, indolalkylamine metabolic process and cellular biogenic amine catabolic process were the most significantly enriched. In the category of cellular component (GO-CC), these genes were significantly enriched in the functional terms of nuclear exosome, exosome, and exoribonuclease complex. The results also showed that hydrolase activity and acting on carbon–nitrogen (but not peptide) bonds were the most enriched terms in the category of molecular function (GO-MF). For KEGG pathway enrichment analysis, the result showed that these genes were significantly enriched in the functional terms of glyoxylate and dicarboxylate metabolism, and tryptophan metabolism.
Analysis of DEGs in dataset GSE25192
The results of the GO function analysis of the DEGs between PBS and CPT11 are listed in Figure 1F. In the category of GO-BP, the genes were significantly enriched in the functional terms of interferon (IFN)-γ-mediated signaling pathway and response to IFN-γ, while in the category of GO-CC, the genes were significantly enriched in the functional terms of lysosomal membrane and lytic vacuole membrane. In the category of GO-MF, these genes were mainly involved in endopeptidase activity, GTPase activity, and so on. There was no enrichment of such genes in KEGG pathway analysis.
The GO function and KEGG pathway enrichment analyses of the DEGs between PBS and PHY906-CPT11 were further performed (Figure 1I). In the category of GO-BP, the genes were significantly enriched in the functional terms of positive regulation of response to external stimulus. The category of GO-CC showed that these DEGs were significantly enriched in collagen-containing extracellular matrix and extracellular matrix. In the category of GO-MF, the DEGs were mainly enriched in chemokine activity, calcium-dependent protein binding, enzyme inhibitor activity, and receptor ligand activity. The results of KEGG analysis showed that these DEGs were mainly enriched in lysosome and apoptosis.
Construction of WGCNA and key module identification
After detecting the outliers using Z.K method with the threshold Z.K score < −2, the GSM618744 was excluded (Figure 2A). Next, we performed the analysis of network topology for thresholding powers from 1 to 30 and identified the relatively balanced scale independence and mean connectivity of the WGCNA. As the lowest power for the scale-free topology fit index on 0.85, power value 8 was selected to produce a hierarchical clustering tree (Figure 2B). Then we used the WGCNA to assign genes with similar expression patterns into one module, and 12 gene modules were obtained (Figure 2C). Interaction relationships of the 12 modules were analyzed, and the network heatmap was plotted (Figure 2D). The results revealed that each module was independent from each other, suggesting a high level of independence among the modules and relative independence of genes expression in each module. WGCNA was then used to correlate each module with the four different treatments (including PBS, PHY906, CPT11, and PHY906-CPT11) for the colon tumors GSE25192 dataset by calculating the MS for each module–trait correlation (Figure 2E). After screening for strong correlations between all modules and PHY906, we found the ME in the red module exhibited a higher correlation with PHY906 than other modules (r = 0.58 and P=2e-04). Meanwhile, we also found that the pink module had the strongest relationship with CPT11. The black module was positively correlated with the PHY906-CPT11, while the greenyellow module was negatively correlated with the PHY906-CPT11. To identify genes associated with different chemotherapeutic agents of colon cancer, the red module, pink module, greenyellow module, and black module were selected for subsequent research, respectively.
Functional enrichment analysis, TFs prediction, and identification of hub genes for the red module associated with PHY906
Scatter plots of GS vs. module membership for the red module was shown in the Figure 3A. Then a total of nine TFs were identified in the red module, which included E2F1, PML, SIN3A, NANOG, USF1, FOXA1, GATA2, TCF3, and NFE2L2. The TF-target gene regulatory network was displayed in Figure 3B. The results of the significant functions and pathways in the red module were presented in Figure 3C. It was found that the genes in this module were associated with regulation of cell morphogenesis involved in differentiation, epithelial cell development, extracellular matrix organization, and so forth. There was no enrichment of this module in KEGG pathway analysis. Afterward, the most significant clusters in the red module was shown in Figure 3D. The results showed that there were 36 hub genes in the module. The boxplots demonstrating the correlations between the tumor treatment and the hub genes were shown in Figure 3E. Except for Cox17, Pck2, Psmc6, Sms, Ckmt1, Gtf2b and 1110005a03rik, the expression status of hub genes in the red module were negatively correlated with PHY906 treatment. According to Figure 3F, the results of GO enrichment analysis revealed the hub genes in red module were significantly enriched in the functional items of substrate adhesion-dependent cell spreading, regulation of substrate adhesion-dependent cell spreading, and cell–substrate adhesion in the category of GO-BP. The results showed the hub genes were significantly enriched in the functional items of basement membrane, cortical cytoskeleton, and actin cytoskeleton in the category of GO-CC. In the category of GO-MF, extracellular matrix structural constituent, nuclear hormone receptor binding, and hormone receptor binding were the representative enrichment terms. For KEGG pathway enrichment analysis, the hub genes in the red module were mainly related to focal adhesion, arginine and proline metabolism, and amoebiasis.
Analysis of the red module negatively associated with PHY906
Functional enrichment analysis, TFs prediction, and identification of hub genes for the pink module associated with CPT11
For the pink module, scatter plots of GS vs. module membership was shown in the Figure 4A. The TF-target gene regulatory network was displayed in Figure 4B, and we identified eight TFs in the module, including MYC, PRARD, RFX5, USF2, STAT3, IRF1, RUNX1, and IRF8. Next, the results of the functional enrichment analysis were listed in Figure 4C. The results showed that the genes in the pink module were significantly enriched in the terms of cytokine-mediated signaling pathway, response to virus, and response IFN-γ in the category of GO-BP, while the category of GO-CC showed that the genes in the pink module were significantly enriched in symbiont-containing vacuole, extracellular membrane-bounded organelle, and host cell cytoplasm. Moreover, phosphoprotein binding, CCR chemokine binding, and protein phosphorylated amino acid binding, and GTP binding were the representative enrichment terms in the category of GO-MF. In KEGG pathway analysis, the genes were significantly involved in several pathway, including influenza A, measles, hepatitis, and JAK-STAT signaling pathway, and so forth. For the pink module, we identified 35 hub genes (Figure 4D). As shown in Figure 4E, all the hub genes were significant in distinguishing CPT11 and other treatments for colon cancers (P<0.05), among which Lrs1, Pls1, Ank2, Adam10, Grpel2, Emp2, Adil, Pafah2, Hey1 and Aif1l were up-regulated while Bloclsl, Psmb8, Ostc, Rtp4, Gm12250, Mrpl3, Gm12250, Mrpl3, Gmppb, Lfitm3, Gbp2 and Epstl1 were down-regulated. The results of the GO analysis of the hub genes were shown in Figure 4F, the results revealed the hub genes were significantly enriched in the functional items of response to IFN-β, symbiont process, and regulation of viral genome replication in the category of GO-BP. Then, in the category of GO-CC, the hub genes had significant enrichment in function of symbiont-containing vacuole, extracellular membrane-bounded organelle, and host cell cytoplasm. Moreover, the hub genes were mainly involved in the functional terms of threonine-type endopeptidase activity in the category of GO-MF.
Analysis of the pink module negatively associated with CPT11
Functional enrichment analysis, TFs prediction, and identification of hub genes for the black module associated with PHY906-CPT11
For the black module, scatter plots of GS vs. module membership was shown in the Figure 5A. There was only one TF identified in the black module. It was SPI1 with 14 target genes. The TF-target gene regulatory network was displayed in Figure 5B. Then, the results of the functional enrichment analysis were listed in Figure 5C. We found genes in the black module most significantly enriched in the functional terms of positive regulation of chemotaxis in the category of GO-BP. In the category of GO-CC, collagen-containing extracellular matrix and extracellular matrix were the representative enrichment terms while peptidase regulator activity, endopeptidase inhibitor activity, endopeptidase regulator activity, and peptidase regulator activity were the most enriched terms in the category of GO-MF. The KEGG analysis revealed the genes in the black module were significantly enriched in the functional items of complement and coagulation cascades. In total 38 hub genes were identified in the black module (Figure 5D). The boxplots demonstrating the correlations between the tumor treatment and the hub genes was shown in Figure 5E. The expression status of Rbm25, Phc1, 3632451o06rik, Dag1, and Prkd were negatively correlated with PHY906-CPT11 treatment, while others had positive correlation. The results of functional enrichment analysis of hub genes in the black module were depicted in Figure 5F. We found the hub genes mainly enriched in positive regulation of chemotaxis, regulation of interleukin-6 production, and interleukin-6 production in the category of GO-BP, and significantly enriched the functional items of Golgi cisterna, Golgi cis cisterna, and Golgi stack in the category of GO-CC. In the category of GO-MF, the genes were mainly involved in peptidoglycan muralytic. In addition, KEGG pathway analysis revealed the hub genes were mainly enriched in the following pathways: complement and coagulation cascades, Staphylococcus aureus infection, and osteoclast differentiation.
Analysis of the black module positively associated with PHY906-CPT11
Functional enrichment analysis, TFs prediction, and identification of hub genes for the greenyellow module associated with PHY906-CPT11
For the greenyellow module, Scatter plots of GS vs. module membership was shown in the Figure 6A. A total of 18 TFs were identified in the greenyellow module which included IRF8, NFIC, ZBTB33, RUNX1, USF2, REST, TAF1, YY1, ZNF384, MAX, MYC, USF1, NELFE, E2F6, NRF1, KLF4, E2F1 and TRIM28. The TF-target gene regulatory network was displayed in Figure 6B. The results of the functional enrichment analysis were listed in Figure 6C. In the category of GO-BP, we found that the genes in black module were significantly enriched in the functional items of RNA splicing, while the genes in black module were mainly enriched in DNA replication factor A complex, SMN–Sm protein complex, and spliceosomal complex in the category of GO-CC. In the results of KEGG analysis, we found the genes in the black module were significantly enriched in spliceosome. Next, we identified 22 hub genes in the greenyellow module (Figure 6D). As shown in Figure 6E, all the hub genes were significant in distinguishing PHY906-CPT11 from other treatments for colon cancers (P<0.05). Except for Skint5, Clic1, Anxa2, and Prss12, the expression status of hub genes in the greenyellow module were negatively correlated with PHY906-CPT11 treatment. The results of functional annotation in the category of GO-CC showed that the hub genes in the greenyellow module were mainly enriched in spliceosomal complex, catalytic step 2 spliceosome, extracellular exosome, and extracellular vesicle. Then, we found that mRNA binding, repressing TF binding, and calcium-dependent protein binding were the representative enrichment in the category of GO-BP. Furthermore, for KEGG analysis, the hub genes in the greenyellow module were significantly enriched in terms of spliceosome (Figure 6F).
Analysis of the greenyellow module positively associated with PHY906-CPT11
Merging of DEGs with hub genes and their functional enrichment analysis
In order to further reveal the key genes involved in the treatment of colon cancer with different drugs (PHY906, CPT11, and PHY906-CPT11), we merged the acquired DEGs and the hub genes obtained through WGCNA in different groups, respectively. For PHY906 group, we merged the DEGs with the hub genes associated with PHY906. Then, we found no key genes after the merging. Next, we merged the DEGs acquired from PBS and CPT11 with the hub genes which for the pink module associated with CPT11 (Figure 7A). In total, three key genes were identified which included Ligp1, Gm12250, and Gbp2. As shown in Figure 7B, we performed functional enrichment analysis for the three genes. The enrichment in several GO-BP terms such as defense response to protozoan, response to protozoan, cellular response to IFN-β, response to IFN-β, and defense response to bacterium were detected. In the category of GO-CC, symbiont-containing vacuole membrane, symbiont-containing vacuole, extracellular membrane-bounded organelle, and host cell cytoplasm were the significantly enriched GO term. In the category of GO-MF, GTPase activity, GTP binding, purine bibonucleoside binding and purine nucleoside binding were the most enriched items. As for KEGG pathway analysis, these genes were mostly associated with NOD-like receptor signaling pathway.
Identification and functional enrichment analysis of key genes
For the PHY906-CPT11 group, after merging the DEGs and the hub genes, we finally identified 18 key genes including Eif4e, Prr15, Anxa2, Ddx5, Tardbp, Skint5, Prss12, Hnrnpa3, Pfdn2, Bcl2a1c, Ncf1, Rgs10, C3ar1, Gng11, Tyrobp, Tmsb4x, C1qc, and Lyz2 (Figure 7C). The results of functional enrichment analysis of these key genes are shown in Figure 7D. In the category of GO-BP, superoxide anion generation was the most significantly enriched GO term. The key genes were significantly enriched in the functional items of rough endoplasmic reticulum, Golgi cis cisterna, NADPH oxidase complex, and G-protein β/γ-subunit complex in the category of GO-CC. In the category of GO-MF, we found the key genes had significant enrichment in function of peptidoglycan muralytic activity, complement binding, G-protein β-subunit binding, and actin monomer binding. Furthermore, the results of KEGG analysis revealed the key genes were significantly enriched in terms of complement and coagulation cascades, Staphylococcus aureus infection, and osteoclast differentiation.
Validation of the key genes
Although we finally obtained 21 key genes, 5 of them were the genes of murine. Therefore, we only performed the validation for 16 genes which included GPB2, EIF4E, PRR15, ANXA2, DDX5, TARDBP, PRSS12, HNRNPA3, PFDN2, NCF1, RGS10, C3AR1, GNG11, TYROBP, TMSB4X, and C1QC. For GPB2, a key gene associated with CPT11 treatment, its expression status was significantly different between normal and tumor samples whether in COAD or READ (Supplementary Figure S2A,D). However, expression status of GPB2 was not related to both the tumor stages and overall survival in colon cancer which included COAD and READ (Supplementary Figure S2B,C,E,F). As for the key genes associated with PHY906-CPT11 treatment, we found that EIF4E, PRR15, ANXA2, PRSS12, HNRNPA3, PFDN2, RGS10, GNG11, and TMSB4X showed obviously different expression levels between normal and tumor samples in colon cancer (Supplementary Figures S3 and S4). Besides, as shown in Supplementary Figures S5 and S6, the expression of EIF4E, TYROBP, and PFDN2 were significantly related to the tumor stages in COAD. Next, the results for overall survival analysis was shown in Supplementary Figures S7 and S8, suggesting these 15 key genes might not correlate with OS of colon cancer patients.
To verify the results of the bioinformatics analyses, we conducted qRT-PCR to detect expression of the key genes in SW480 cells after treating with different treatments. As shown in Figure 8, the expression of EIF4E, HNRNPA3, and PFDN2 in SW480 cells of PHY906-CPT11 group were significantly lower than those of other groups (control, PHY906, and CPT11 groups). In the meanwhile, compared with other treatments (PBS or PHY906 or CPT11) the expression of ANXA2, RGS10, GNG11, C3AR1, and TMSB4X in SW480 cells were significantly increased after treating with PHY906-CPT11. Besides, we also found that the expression of GBP2 in CPT11 and PHY906-CPT11 groups were both significantly higher than that in control and PHY906 groups. Collectively, the results of qRT-PCR were similar with the validation results described as above. The EIF4E, HNRNPA3, PFDN2, ANXA2, RGS10, GNG11, C3AR1, and TMSB4X genes might play a role in improving the therapeutic index of CPT11.
RT-PCR analysis of 16 key genes
Colon cancer is one of the most common causes of cancer deaths worldwide, and its incidence has become higher and higher in recent years. CPT11 is one of the most effective chemotherapeutic drugs for colon cancer patients, but its severe side effects remain unresolved. PHY906 has been commonly used as an adjuvant agent for the treatment of various cancers, especially for colon cancer . Several researches revealed that PHY906 not only enhances the anti-tumor efficacy of CTP11 but also alleviates chemotherapy-induced side effects (diarrhea) in the treatment of advanced colon cancer .
In order to provide novel insight to understand the molecular mechanism of different drugs in treating colon cancer, the microarray data from the colon tumors treated with PHY906, CPT11, and PHY906-CPT11 were further analyzed. Initially, we found that only three DEGs were screened out between PBS and PHY906; this might be because PHY906, as an adjuvant in cancer treatment, has no significant effect when treating colon cancer alone . Moreover, 131 DEGs were identified between PBS and CPT11, and 268 DEGs were identified between PBS and PHY906-CPT11, which indicated that the CPT11 and the PHY906-CPT11 administration both induced significant alterations in gene expression levels of colon tumors when compared with PBS treatment. Then, we applied WGCNA to identify the key modules and hub genes in therapy of colon cancer by R and identified four key modules that were most significantly associated with γ three different treatments. Since TFs play important regulatory roles in human diseases, detailed analyses of these regulators will allow us to better understand the molecular mechanisms underlying colon cancer therapy; so, we also constructed TF-gene regulation networks within the key module. The results showed us that the red module mostly correlated with PHY906 treatment in colon cancer, and nine TFs (E2F1, PML, SIN3A, NANOG, USF1, FOXA1, GATA2, TCF3, and NFE2L2) were significantly associated with the genes within red module. Several studies have shown that E2F1 has potential for mediating multiple cancer hallmarks, such as sustaining proliferative signaling, metabolic reprogramming, and tissue invasion and metastasis [32,33]. NANOG is involved in the self-renewal of embryonic stem cells; a previous study indicated that NANOG gene is abnormally overexpressed during the development of malignant phenotype of tumor cells . FOXA1 plays a regulatory role in the evolution and development of organisms, and the recent studies revealed that FOXA1 is associated with a variety of cancers, such as prostate, breast, and gastric [35–37]. After filtering for GS and MM value, we eventually obtained 36 hub genes (Loc100041932, 1300001i01rik, Kdelr2, Lasp1, 1110005a03rik, Llgl2, Copa, Orai1, Bcar1, Actn4, Col4a2, Jmjd5, Arrb1, Cry2, Gtf2b, Tmem64, Zfp319, Myadm, Sparc, Ehmt1, Ctgf, Litaf, Psmc6, Sms, Ckmt1, Oat, Zyx, Col18a1, Cox17, Enpp5, Prickle1, Pck2, Slc44a4, Ppap2b, Fn1, and Tpp1) and their expression were significant in distinguishing PHY906-CPT11 and other treatments for colon cancers. Many studies have reported that some of these 36 genes such as Lasp1, Llgl2, Bcar1, Actn4, Jmjd5, Sparc, Ctgf, Litaf, Ckmt1 and Pck2, are cancer-related genes which function in tumorigenesis and prognosis. As for other hub genes, there are few reports implicating them in cancers. LASP1 is highly expressed in various malignant tumors such as colon cancer and prostate cancer, and is closely related to tumorigenesis, invasion, and metastasis [38,39]. BCAR1 is an inflammatory gene that has been shown to be a novel biomarker for colorectal cancer prognosis . A large number of clinical studies suggested that the changes of ACTN4 expression are related to the aggressiveness, invasiveness, and metastasis in some tumors [40,41]. It has been reported that JMJD5 was important for sustaining cell migration and invasion, and could be a potential oncogene for colon carcinogenesis . SPARC is highly expressed in many tumors, and a research showed that its expression is related with metastasis and recurrence of colon cancer [43,44]. CTGF has shown to be associated with progression, tumorigenesis, and prognosis of various cancers [45,46]. To explore biological functions for the hub genes, we performed the functional enrichment analysis, and the results showed these hub genes were mainly enriched in substrate adhesion-dependent cell spreading and focal adhesion. Focal adhesion is demonstrated to be an essential pathway for BPs including wound healing and tumor metastasis . Furthermore, we merged the DEGs and hub genes which associated with PHY906 treatment and found no common gene. This may be because PHY906 does not have a particularly obvious therapeutic effect on colon cancer, which is similar to the previous DEGs screening results.
For the CPT11 treatment, we identified the pink module as the most significantly related module. Different from the red module associated with PHY906 treatment, in addition to the enrichment in response IFN-γ, the genes in the pink modules were also significantly enriched in cytokine-mediated signaling pathway, and response to virus. In pathway analysis, genes in the pink module was found to be significantly related to JAK-STAT signaling pathway and herpes simplex virus 1 (HSV-1) infection. JAK-STAT pathway is widely involved in cell proliferation, differentiation, apoptosis, and inflammation. The activation of JAK-STAT pathway promotes the occurrence and development of various diseases, including various inflammatory diseases, lymphoma, leukemia, and solid tumors [48,49]. HSV-1 has reportedly been shown to be the most effective oncolytic virus, which can stimulate the body to produce a strong anti-tumor immune response. A previous study demonstrated that HSV-1 has a good clinical effect on colon cancer . TF-gene regulation networks within the pink module had been conducted, in total eight TFs were identified which including MYC, PRARD, RFX5, USF2, STAT3, IRF1, RUNX1, and IRF8. Except for the PRARD, other TFs were previously reported to be related to cancer progression. An increasing body of evidence supports that the cytoplasmic TF STAT3 is activated constitutively in a variety of cancers wherein it significantly affects the growth of tumors and facilitates metastasis [51,52]. USF2 is observed in various human cancers, plays important roles in embryogenesis, metabolism, and cancer development . A previous study indicated that USF2 is associated with tumor grade and inversely with survival in Stage II colon cancers . IRF1 and IRF8 are member of the IFN regulatory factor family; both of them were proven to be suppressor of colon cancer [55,56]. Fijneman and colleagues found a set of RUNX1 target genes in mouse colon indicative of changes in gut homeostasis, that is, genes involved in inflammation and intestinal metabolism, which are known to increase the risk of tumor development . Then, we obtained 35 hub genes within the pink module, and there were 23 cancer-related genes (Lsm4, Arpp19, Rtp4, Ifitm3, Flywch1, Ppia, Bre, Psmb8, Ank2, Mrpl3, Hey1, Impa2, Emp2, Irs1, Psmb9, Epsti1, Ccl5, Bloc1s1, Gbp2, Adam10, Adi1, Aif1l, and Dusp1), which functioned in tumorigenesis and prognosis. FLYWCH1 has the potential to become a biomarker and therapeutic target against metastatic colon cancer . Some studies demonstrated that PPIA expression was up-regulated in all 17 types of cancer assessed when compared with normal tissues . The expression level of IRS1 was found to be higher in colon cancer cells in comparison with the normal cell . Compared with the hub genes associated with PHY906 treatment, there were more cancer-related genes in the hub genes associated with CPT11 treatment. This result may explain why the anti-cancer effect of CPT11 was significantly stronger than that of PHY906. However, the results of the functional enrichment analysis of the 35 hub genes in the pink module seemed to have no obvious correlation with cancer progression. In addition, the overlap of the DEGs with these hub genes allowed the identification of a common list of three key genes which included Ligp1, Gm12250 and Gbp2. Similarly to our results, a previous study demonstrated that GBP2 might be involved in the progression of colon cancer . GO and KEGG pathway enrichment analysis showed that the key genes related to CPT11 treatment were significantly enriched in defense response to protozoan, cellular response to interferon-beta and NOD-like receptor signaling pathway. NOD-like receptor is an inflammatory immune receptor which is connected with many organelles in the cell. In recent years, it is generally believed that inflammation is closely related to tumors, especially chronic inflammation, and even the main driving factor of many tumors [62,63]. Therefore, the treatment of colon cancer with CPT11 is potentially regulated by LIGP1, GM12250 and GBP2, and exerts anti-cancer effects through mediating the NOD-like signaling pathway.
In order to clarify the anti-cancer mechanism of PHY906-CPT11 treatment more clearly, we constructed the WGCNA and screened out the two modules with the most positive and negative correlations. Black module was most positively correlated with PHY906-CPT11 treatment. The genes in the black module were significantly enriched in positive regulation of chemotaxis, which could relate to the cancer invasion and metastasis. Then, a total of 38 hub genes were obtained from the black module. We found that 30 of these genes which included Nckap1l, Ncf1, Pla2g15, Unc93b1, Lrrc33, Rgs10, Gmfg, Rbm25, C3ar1, Ccdc88b, Fam105a, Gng11, Phc1, Rasa4, Dag1, Tyrobp, Tmsb4x, Cfp, Prkd1, Gp49a, Trem2, C1qc, Lyz2, Igfbp7, Rab30, Cxcl16, Egfl7, Ube2r2, Pltp, and Cfb, showed a relation with cancer. Wang and colleagues first assessed GMFG expression in colorectal cancer cell lines and tissue specimens and found the targeting of GMFG expression may suppress colorectal cancer progression . It has been reported that CCDC88B has a critical function in colon inflammation and the pathogenesis of IBD , and a previous study proved that patients with IBD have a higher probability of developing colorectal cancer than other people . Strong expression of the CXCL16 protein in primary colon tumor tissue was reported to be correlated to a good prognosis . EGFL7 plays a vital role in controlling vascular angiogenesis during embryogenesis, organogenesis, and maintaining skeletal homeostasis, its dysregulation has been frequently found in several types of cancers [68–70]. Accumulating evidence suggested that EGFL7 plays a crucial role in cancer biology by modulating tumor angiogenesis, metastasis, and invasion . Afterward, we identified the yellowgreen module as the key module which most negatively correlated with PHY906-CPT11 treatment. A total of 18 TFs (IRF8, NFIC, ZBTB33, RUNX1, USF2, REST, TAF1, YY1, ZNF384, MAX, MYC, USF1, NELFE, E2F6, NRF1, KLF4, E2F1, and TRIM28) were identified in this module. Interestingly, we found that E2F1 was also identified in the red module associated with PHY906, while IRF8, RUNX1, and USF2 were also identified in the pink module associated with CPT11. Our study found that genes in the yellowgreen module were significantly enriched in terms of spliceosome. Recently, increasing studies have found that spliceosome has a very important position in the occurrence and treatment of cancer . In the yellowgreen module, a total of 22 hub genes were identified, which included Eif4e, Prr15, Clic1, Anxa2, Ddx5, Tardbp, Pde4dip, Crhr2, Scd2, Skint5, Naa35, Slc25a3, Prss12, Snrpd3, Gmnn, Aifm1, Hnrnpa3, Rcc2, Pfdn2, Rdh9, Dph5, and Fam100b. We found that 17 hub genes were cancer-related genes, which were involved in tumorigenesis and prognosis. Wang and colleagues  first demonstrated that CLIC1 regulates the migration and invasion of colon cancer. Several studies found Anxa2 overexpression in different types of tumors including colon cancer . A new research found that ANXA2 mRNA is up-regulated at all stages of colon cancer and ANXA2 protein levels associate with high probability of invasion and distant metastasis . In a previous study, CRHR2 was identified as a gene which contributing to reversal of colorectal cancer cell resistance . Consistent with the result of functional enrichment analysis for the yellowgreen module, the hub genes were also significantly enriched in terms of spliceosome. Compared with the previous two treatment, we found that there were significantly more cancer-related genes involved in PHY906-CPT11 treatment. This is consistent with the results of the comparison of the current treatment effects of the three different treatments, the best in PHY906-CPT11, the second in CPT11, and the worst in PHY906. Finally, we merged the DEGs and the hub genes obtained through WGCNA to find key genes associated with PHY906-CPT11 treatment. A total of 18 key genes were finally obtained, which included EIF4E, PRR15, ANXA2, DDX5, TARDBP, SKINT5, PRSS12, HNRNPA3, PFDN2, BCL2A1C, NCF1, RGS10, C3AR1, GNG11, TYROBP, TMSB4X, C1QC, and LYZ2. Except for SKINT5, other key genes are cancer-related genes and most of these cancer related genes are associated with progression, tumorigenesis, and prognosis of colon cancer. We performed a functional enrichment analysis of key genes, and the results showed that these genes were mainly enriched in the functional terms related to superoxide anion generation, complement and coagulation cascades, and Staphylococcus aureus infection. Some studies found an imbalance between the generation of reactive oxygen species (ROS) and the ability of a cancer cell to repair the resulting oxidative damage leads to a state of oxidative stress that can culminate in apoptosis . Complement can mediate immune response and inflammatory response. Complement coagulation cascades is involved in many physiological and pathological processes, including inflammatory processes. Increasing evidence showed that variation in the type and number of intestinal flora plays an important role in the occurrence and development of colorectal cancer [78,79]. Intestinal flora imbalance might induce bacterial infections and inflammation, which can lead to cancer. Therefore, we speculated that these key genes may be the important therapeutic targets for inducing colon cell apoptosis by producing ROS, and inhibiting bacteria and inflammation.
Furthermore, we verified the expression of key genes in TCGA datasets using GEPIA. Before validation of the key genes, we excluded five murine genes which included IIGP1, GM12250, SKINT5, BCL2A1C and LYZ2. The dysregulation of GBP2, EIF4E, PRR15, ANXA2, PRSS12, HNRNPA3, PFDN2, RGS10, GNG11, and TMSB4X were found in the samples of colon cancer, suggesting that these genes are related to progression of colon cancer. In the analysis for relation of gene expression and tumor stages, we found only three of the key genes were significantly related to tumor stages of COAD. Unfortunately, the expression of all key genes we identified were not related to overall survival of colon cancer patients. Besides, we analyzed the expression of 16 key genes by qPCR in SW480 cells treated with different drugs (PHY906, CPT11 and PHY906-CPT11). The results showed that the GBP2, EIF4E, PRR15, ANXA2, HNRNPA3, NCF1, C3AR1, PFDN2, RGS10, GNG11, and TMSB4X were significantly dysregulated in SW480 cells after treating with PHY906-CPT11, which was similar with the above results. Thus, we speculated that PHY906-CPT11 exerts therapeutic role in colon cancer by regulating the expression of GBP2, EIF4E, PRR15, ANXA2, HNRNPA3, NCF1, C3AR1, PFDN2, RGS10, GNG11, and TMSB4X. Previous evidence has indicated that these genes play regulatory roles with other proteins in several tumors including colon cancer [75,80,81].
However, there are some limitations in the present study. First, small sample size limits the statistical power to identify the key genes, but we validated the key genes based on GEPIA and PCR analysis to improve the credibility of our results. Besides, further studies about the genes obtained in the present study are required for elucidating the underlying mechanism of PHY906-CPT11 in colon cancer treatment. Since the NF-κB plays an important role in colon cancer progression, its involvement in PHY906-CPT11 treatment in colon cancer should also be elucidated by in vivo studies in future works.
Nowadays, researchers have increasingly revealed that combination of chemotherapy and TCM would improve the survival of patients, but most of the molecular mechanisms are not clear. In the present study, we applied bioinformatics method to screen out the genes involved in therapeutic process of colon cancer. We identified a total of 18 key genes (EIF4E, PRR15, ANXA2, DDX5, TARDBP, SKINT5, PRSS12, HNRNPA3, PFDN2, BCL2A1C, NCF1, RGS10, C3AR1, GNG11, TYROBP, TMSB4X, C1QC, and LYZ2) associated with PHY906-CPT11 therapy, and found that the expression of EIF4E, PRR15, ANXA2, HNRNPA3, NCF1, C3AR1, PFDN2, RGS10, GNG11, and TMSB4X in SW480 cells was affected by PHY906-CPT11 treatment. Functional enrichment suggested that these key genes might induce cell apoptosis by producing ROS, and inhibit bacteria and inflammation to play a therapeutic role in colon cancer. Although our investigation is of a preliminary nature and further research is needed to verify these findings, our work might provide direction for specific molecular mechanisms of combination of chemotherapy and TCM in future studies.
All data generated or analyzed during the present study are included in this published article.
The authors declare that there are no competing interests associated with the manuscript.
The authors declare that there are no sources of funding to be acknowledged.
S.X. analyzed the data and wrote the manuscript. K.H. contributed to design of the work. Y.W. and F.W. obtained the data from database. T.S. and Q.L. performed the WGCNA analysis and constructed the co-expression network. All the authors read and approved the final version of the manuscript.
body mass index
C-C chemokine receptor
differentially expressed gene
false discovery rate
Gene Expression Omnibus
Gene Expression Profiling Interactive Analysis
herpes simplex virus 1
inflammatory bowel disease
Kyoto Encyclopedia of Genes and Genomes
Molecular Complex Detection
Nucleotide binding oligomerization domain
phosphate buffered saline
quantitative reverse transcription polymerase chain reaction
reactive oxygen species
reverse transcription polymerase chain reaction
The Cancer Genome Atlas
traditional Chinese medicine
weighted gene co-expression network analysis