Data Information and Clinical Impact of KRAS Mutation in CRC
We downloaded the information from Colorectal Adenocarcinoma (TCGA, Firehose Legacy) Samples with mutation and CNA data (220 samples/patients) with complete follow-up profiles. There were colorectal carcinoma patients with the 44 % KRAS mutation, and mutation types included missense mutations, Amplification, deep deletion, and KRAS mRNA expression heatmap of Colorectal Adenocarcinoma was shown in Figure 1A.
We next investigated the influence of KRAS mutation on CRC progression and prognosis. Results showed that Alteration event frequency in COAD, and negative correlation with KRAS mutation of CRC patients (Figure 1B, 1C). Analysis of the relationship between the KRAS mutation status and disease prognosis showed that patients with the KRAS mutation had poorer prognosis on survival (Figure 1D), which indicated that KRAS mutation may contribute to CRC disease progression. Furthermore, the KEAP1 mutation status and disease prognosis showed that patients with the KEAP1 mutation (MGC9454) had poorer prognosis on survival in CRC progression (Figure 1E).
Identification of DEGs
To further investigate the KEAP1 activation in KRAS mutant colorectal cancer progression, we identified the DEGs among the transcriptional data and corresponding clinical data of CRC samples (GSE39582). The limma package for examining differential expression of RNA-Seq count data, was used according to GSE39582 by the user’s guide for GEOR2. Based on the following criterion: |fold change (FC)| >1; both the P-value and FDR <0.05, a total of 2489 genes were identified as DEGs of GSE39582 or GSE31084. The volcano plot, comparison of the expression of the DEGs in the CRC samples (GSE39582, GSE31084) is shown in Figure 1F, 1G.
Differential expression analysis and comparison of KRAS and KEAP1-related DEGs
To analyze the KRAS and KEAP1 related genes between the DEGs in the CRC samples (GSE39582), KRAS and KEAP1 related CRC genes were checked by the GeneCards database (https://www . genecards.org/). Considered that the significance of the role of these original key genes in HCT116 CRC cells may be different, thus, we also matched these genes with the predicted genes and targets were insected for DEGs in GSE31084 (HCT116 CRC cells with KRAS G13D Mutation), DEGs in GSE39582 and KEAP1 CRC in GeneCard were intersected to obtain the key genes. 3050 KRAS mutant related CRC genes and 1099 KEAP1 related CRC genes were downloaded and used for differentially expressed KRAS and KEAP1 related genes were intersected by Venn diagram. After taking the intersection by Venn diagram, 28 differentially expressed KRAS and KEAP1 related genes were finally identified (Figures 2A).
Based on RNA-seq results, a total of 28 KRAS and KEAP1-related differentially expressed genes were identified in the GEO (GSE39582). After screening the genes with differential expression in GSE39582, 13 up-regulated genes and 15 down-regulated genes with differential expression were determined and constructed heatmaps (Figure 2B) and then created a PCA diagram (Figure 2C) to obtain all the up or down-regulated DEGs among the two GEO datasets. Then, the interaction between KRAS and KEAP1-related 28 differentially expressed genes in the PPI network were constracted in STITCH (Figure 2D).
Pearson correlation of KRAS and KEAP1-related DEGs
In addition, through the Pearson correlation test, there were correlations between the KRAS and KEAP1-related DEGs at the mRNA levels (Figure 3A, 3B). Pearson correlation test between 20 KRAS and KEAP1-related DEGs, including MET, IGFBP3, GDF15, FOS, RNF43, CDKN1A, SPP1, CXCL1,GALNT12, CD44, TIMP3,SLC7A11, BCL2L11, EREG, CTNNB1, BCL2L1, IL18, CEACAM1, CDKN2B, HPSE.
GO and KEGG Analyses
In order to analyze the functional level, we submitted all the Pearson correlation of KRAS and KEAP1-related 20 DEGs were identified for were identified for GO Classification analysis, suggested that the significantly enriched in biological adhesion, biological regulation, cell killing, cellular component organization or biogenesis, cellular process, developmental process, establishment of localization, growth, immune system process, localization, locomotion, metabolic process and so on (Figure 3C, 3D).
Furthermore, in the KEGG Pathway Classification analysis, significantly enriched in Cancer: overview, Signal transduction, Immune system, Cancer: specific types, Infectious disease: viral, Endocrine system, Drug resistance: antineoplastic, Infectious disease: bacterial and so on (Figure 3E). In the KEGG enrichment analysis (Enrichment score>10) of Pearson correlation of KRAS and KEAP1-related 20 DEGs, significantly enriched in regulation of apoptotic signaling pathway, FoXO and PI3K-Akt signaling pathway, Colorectal cancer, MicroRNAs in cancer, and Pathway in cancer (Figure 3F).
Correlation Between the KRAS and KEAP1-related DEGs
In the heatmap-coefficient test, the prominent positive or negative correlations between the KRAS and KEAP1-related DEGs, such as SPP1 significant positive correlations with FOS or HPSE (Correlate score > 0.80). Moreover, SPP1, FOS and HPSE significant negative correlations with CEACAM1, CXCL1, SLC7A11, GDF15, CD44, CDKN1A, IL18, METI, GFBP3 (Correlate score < -0.80) (Figure 4A). Most important, CD44, CDKN1A, IL18, MET or CEACAM1, CXCL1, SLC7A11, GDF15, IGFBP3 showed the significant both positive and negative correlations with KRAS and KEAP1-related DEGs. Furthermore, Curve of samples in the GEO (GSE31084) indicated that the mRNA expression of 20 DEGs in the positive and negative correlations (Figure 4B).
The samples of the GEO (GSE31084) indicated the up-regulated and down-regulated genes with differential expression were determined and constructed heatmaps (Figure 4C, 4D). Then, the Boxplot for FPKM mRNA expression to obtain all the up or down-regulated DEGs among the samples of GSE31084 (Figure 4E), and a Hierarchical clustering showed that mRNA expression levels are correlated with 20 DEGs among multiple types of samples, and mRNA expression of DEGs is associated with five types of samples among CRC (Figure 4F).
The Relationship between hub genes and Clinical Features in COAD
Then, we constructed a KRAS and KEAP1-related gene prognostic signature comprising 20 key genes in the training cohort. Based on their median risk score, COAD patients were categorized into high-risk and low-risk groups.
Similarly, we performed analyses in the validation cohort and entire cohort. Figure 5 showed that KRAS and KEAP1-related key gene expression profiles in the validation cohort and entire cohort of TCGA-COAD. Furthermore, stratification analyses based on low expression of CDKN2A, SPP1, FOS, and BCL2L11 showed that the OS of low-risk patients was better than that of high-risk patients (Figure 5C, 5F, 5I, 5L). According to the KM survival curves, high-risk of CDKN2A, FOS, and BCL2L11 patients showed a lower OS than low-risk patients (Figure 5C, 5I, 5L), illustrate the distribution of risk scores and survival status of COAD patient. In the training cohort, the areas under the ROC curves (AUCs) for the
KRAS and KEAP1-related hub gene prognostic signature, we applied in the “ROC plotter” to validate our findings. The AUCs of CDKN2A, SPP1, FOS, BCL2L11 and HPSE were 0.88-1.00, respectively, shown in Figures 5B,5E, 5H, 5K.. The AUCs curve analysis indicated that the SPP1, FOS of high-risk patients was significantly worse than that of low-risk patients, which was significant correlations with KRAS and KEAP1-related DEGs in our correlation analysis.
Correlation between the key hub gene prognostic signature and clinical characteristics
To determine the mechanism by which KRAS and KEAP1-related hub gene affected the prognosis of COAD patients, we used the UALCAN online database to analyze 374 COAD samples and 41 normal samples in TCGA according to molecular subtypes and histological classification. Next, Figure 6A demonstrates the association of clinical characteristics with prognostic gene expression levels among the KRAS and KEAP1-related key hub genes (CDKN2A, SPP1, FOS, BCL2L11 and HPSE).
We further explored the value of the prognostic signature in the two risk groups stratified by different clinical characteristics (gender and stage). CDKN2A, SPP1, FOS, BCL2L11 expression was both significantly higher in COAD samples than in normal samples, gender and the COL1A1 expression did not differ significantly, and the COL1A2 expression and gender were not significantly correlated (Figure 6B). Additionally, CDKN2A, SPP1 and FOS expression significant correlations were found between stage I and stage III. The CDKN2A, SPP1 and FOS expression of stage I, II and stage III were significantly different, the CDKN2A, SPP1 and FOS expression of stage IV and stage III were significantly different. Furthermore, HPSE expression was both significantly lower in COAD samples than in normal samples, gender and the HPSE expression did not differ significantly, the HPSE expression of stage I, II and stage III, IV were no significantly different (Figure 6C). These results support that KRAS and FOS regulators are closely related to the Raf/MEK/ERK signaling pathway in colorectal cancer (Figure 6D).
Correlation Between Immune Infiltrates and hub genes in COAD
The HPA database was then used to further study the protein expression of CDKN2A, SPP1, FOS, BCL2L11 and HPSE at the translation level. Meanwhile, CDKN2A and FOS expression showed high expression in Adenocarcinoma and Mucinous adenocarcinoma than in normal (Figure 7A, 7E). Also, SPP1 expression decreased in Adenocarcinoma and Mucinous adenocarcinoma with a histological subtype dependent mannerh (Figure 7C).
To determine if hub genes have an effect on the immunological milieu of tumors, we examined the connection between CDKN2A, SPP1, FOS, BCL2L11 and HPSE expression and the degree of immune cell infiltration in TCGA-COAD. We indeed found FOS was strong association in multiple malignancies using the infiltration scores of six immune cell types (e.g., CD8 + T cells, CD8 + Tex, CD4 + T cells, B cells, neutrophils, macrophages, and dendritic cells) accessible in the Tumor Immune Estimation Resource (TIMER) database and obtained from TCGA (Figure 7F). Meanwhile, SPP1 was positively correlated with the extent of immune cell types (macrophage, neutrophil and dendritic cell) (Figure 7D).
The Key Genes and Targeted for the predicted Herbal Ingredients
In this study, an interesting gene of secreted phosphoprotein 1 (SPP1), which the Target Related Ingredients is Cyclopamine in the Herb database (Supply Table 1), and we found the Paper mined target genes of Cyclopamine including SRY-box transcription factor 9 (SOX9). As a member of SOX family, SOX9 can regulate cancer progression by affecting the transcription of diverse genes [6,7], and promote Oxaliplatin resistance in colorectal cancer [8]. Intestinal adenomas and colorectal adenocarcinomas from mouse models and patients demonstrate ectopic and elevated expression of SOX9. Based on SPP1, which was presented in DEGs of GSE39582, KRAS and KEAP1 related genes, thus, we considered it as the key target to SOX9, we directly predicted herbal ingredients such as Cyclopamine and honokiol in HERB (Supply Table 2). The chemical structure and DEGs of honokiol is shown in Figure 8A, 8B. The GO and KEGG analysis of honokiol, suggested that the pathways regulated by honokiol are the significantly enriched in Estrogen, GnRH signaling pathway and AGE-RAGE signaling pathway(Figure 8C, 8D). The chemical structure and Differentially Expressed Genes (DEGs) of cyclopamine is shown in Figure 8E, 8F. The GO and KEGG analysis suggested that the pathways regulated by cyclopamine, are the significantly enriched in hedgehog signaling pathway, Basel cell carcinoma, cAMP signaling pathway (Figure 8G, 8H).
Molecular docking studies of hit compounds in the cancer-associated target genes
The Glide module of the Schrodinger suite software was used for calculating the docking scores of propolis hit compounds cyclopamine and honokiol against the active sites of cancer-associated target genes SOX9. From Figure 9, it can be observed that cyclopamine or honokiol had the stabilized docking score against SPP1. Whereas, cyclopamine or honokiol exhibited the most stabilized interaction with SOX9.
The 2D and 3D interaction diagrams of honokiol in the active site of SOX9 (PBD ID 4euw) revealed the formation of two hydrogen bonds between 3 and 4′hydroxyl groups with DT2 and DC3 (Figure 9A), while cyclopamine in the active site of SOX9 revealed there was the formation of hydrogen bonds between 3′-hydroxyl groups (Figure 9D). The 2D and 3D interaction diagrams of honokiol in the active site of SPP1 (PBD ID 6j2p) revealed the formation of two hydrogen bonds between 1′-, 3′- and 7′-hydroxyl groups with ASP40, PHE62 and LYS58 (Figure 9B), while cyclopamine in the active site of SPP1 revealed there was the formation of hydrogen bonds between TRP45 and ASP40 (Figure 9C).
In Vitro Validation of Predicted Targets in CRC Cells
Cell death in HCT116 cells
To determine the cytotoxic effect of cyclopamine or honokiol on the proliferation of CRC cell lines, we treated colorectal cancer HCT116 cell lines with different concentrations of cyclopamine or honokiol for 12, 24, 36, 48 h and assessed cell viability via MTT assay. The results showed that cyclopamine or honokiol significantly reduced the viability of HCT116 CRC cell lines in a concentration- and time- dependent manner compared to control cells (Figure 10A, 10B).
To investigate whether cyclopamine or honokiol induces apoptosis as previously reported, the cell death was detected by flow cytometry analysis using Annexin V and PI staining. Results showed that cyclopamine or honokiol could induced apoptosis after 24 or 48 hours treatment. Interestingly, for short time treatment (12 hours) with cyclopamine or honokiol in HCT116 cells, a portion of the cells gathered in the area of Annexin V-/ PI+ (Figure 10C). This result suggested that cyclopamine may have induced or may induce the apoptosis type of cell death within 12 hours. To further investigate what kind of cell death was induced by cyclopamine or honokiol in this short time, cell death were tested together with α-hederin. This result indicated that cyclopamine induced cell death within 48 hours is apoptotic cell death, suggesting that apoptosis might be dominant for longer time treatment (Figure 10D). All these results indicated that cyclopamine -induced early cell death might be apoptosis while apoptosis might be induced upon longer time treatment.
Induces KEAP1\Nrf2 activation in CRC cells
In order to investigate the underlying mechanism of cyclopamine or honokiol induced cell death, RNA-sequencing analysis was employed. NFE2L2, is a well-known upstream element of the genes that are up-regulated in response to ROS [14]. Our results showed that cyclopamine significantly decreased expression of NFE2L2 in HCT116 cells at the mRNA level in a concentration-dependent manner compared to honokiol treated cells (Figure 11A). At the same time, we tested the protein levels of KEAP1 and NFE2L2 in HCT116 cells for 24 h after the treatment of cyclopamine or honokiol, and the protein levels of KEAP1 and NFE2L2 decreased (Figure 11B). This result suggested that KEAP1- NFE2L2 signals was inhibited when treated with cyclopamine which might play a critical role in inducing cell death of HCT116 cells, and cyclopamine induced activation of KEAP1\NFE2L2 signaling pathway, which cell death are characteristic features of apoptosis.
SOX9 and SPP1 activation in CRC cells
To investigate whether SOX9 or SPP1 was involved in cyclopamine -induced cell death, HCT116 cells were treated with cyclopamine or honokiol to check the genes of SOX9 and SPP1. As shown in Figure 10F, western blot assay showed that the cyclopamine involved the activation of SOX9 and SPP1 in HCT116 cells decreased within 24 h after treatment (Figure 11C).