3.1. Identification of prognostic IRDEGs and functional analysis
We analysed DEGs in mRNA expression profiles between irinotecan-treated and untreated CRC cell lines using GSE145356 cohort in the GEO database. Total 548 IRDEGs were screened out, of which 233 were upregulated and 315 downregulated (Fig. 1A, B) under the threshold of FDR<0.05 and |log2 FC| > 1. Pathway enrichment analysis was conducted to explore the function of the IRDEGs. The results showed that the IRDEGs to be most enriched for terms related to tumour cell functions and tumour growth, such as cell adhesion, cellular component movement, tissue morphogenesis and the VEGFA-VEGFR2 signalling pathway (Fig. 1C and Supplementary Figure S1). These results indicated that genes related to irinotecan treatment probably play a pivotal role in malignant occurrence and progression of CRC. Subsequently, to screen the IREDGs related to prognosis of CRC patients, univariate Cox regression analysis was performed to find prognosis-related genes in TCGA database. We then intersected IREDGs with CRC prognosis-related genes (Fig. 1D), and 51 overlapping genes were selected for subsequent analyses (Fig. 1E, F). The details of gene expression, interaction and prognostic role among the IRDEGs were shown in the supplemental materials (Supplementary Figure S2).
3.2. Establishment of IRDEG-based subtypes and clinical validation
Based on the expression features of the 51 prognostic IRDEGs, we established irinotecan-related subtypes in CRC using unsupervised consensus clustering. According to the CDF graph and delta region graph (Fig. 2A-C), we separated patients in the CRC cohort from TCGA into two clusters (k value=2): Cluster A and Cluster B. Principal component analysis (PCA) was conducted to validate the dependability of the clustering (Fig. 2D). As the survival curve showed, patients in Cluster A had a significant survival advantage compared with patients in Cluster B (p<0.001) (Fig. 2E). We next investigated the relationship between IRDEG-based subtype and clinical information (Fig. 2F). As the heatmap showed, the N stage, M stage, TNM stage and survival status correlated significantly with the IRDEG-based subtype.
3.3. Construction of the iriscore model based on IRDEG-based subtypes and validation of its prognostic value
We further constructed the iriscore model by PCA based on IRDEG-based subtypes to explore the value of the signature genes in terms of prognosis. Then, we divided CRC patients into a high-iriscore group and a low-iriscore group based on the median value of iriscore. The low-iriscore patients had longer OS than patients with a high iriscore (P < 0.001), which indicated that iriscore had a significant prognostic value in CRC (Fig. 3A). The iriscore results for those who died was more elevated than that of those who survived (P<0.001), and the high-iriscore group had a higher proportion of patients who died (Fig. 3B, C). Areas under the curve (AUCs) were 0.693, 0.640, 0.701 and 0.753 at 1, 3, 5 and 10 years, respectively (Fig. 3D).
To assess the validity of the prognostic value of iriscore, we used a validation cohort (GSE39582) from the GEO database. Consistent with the results above, patients in the low-iriscore group presented a longer OS than those in the high-iriscore group (P<0.001) (Fig. 3E). The incidence of death was also correlated with iriscore strikingly (P<0.05) (Fig. 3F, G). AUCs of the GSE39582 cohort for the 3-, 5-, 10- and 13-year OS rates were 0.544, 0.560, 0.573 and 0.661, respectively (Fig. 3H). Then, subgroup analysis was conducted according to T stage. Based on Kaplan‒Meier curve analysis, patients with T1-2 in the low-iriscore group exhibited no significant differences compared with the high-iriscore group in TCGA CRC (Supplementary Figure S3) and GSE39582 (Supplementary Figure S5) cohorts. In contrast, patients with T3-4 in the low-iriscore group had better prognosis than patients in the high-iriscore group in both TCGA CRC (P<0.001) (Supplementary Figure S4) and GSE39582 (P=0.004) cohorts (Supplementary Figure S6). The similarity between the training dataset and validation dataset indicates that the prognostic value of iriscore is specific and sensitive.
3.4. Close correlation of the tumour mutation burden (TMB) with iriscosre
Recently, the tumour mutation burden (TMB), which means the number of somatic mutations after removal of germline mutations from tumor genomes, has been increasingly proved to have predictive value of tumour response to immunotherapy for patients with metastatic CRC[19]. Therefore, we confirmed the features of TMB in the high-iriscore and low-iriscore groups (Fig. 4A, B) and explored the relationship between TMB and iriscore (Fig. 4C). The TMB of the low-iriscore group presented higher compared to the high-iriscore (P<0.01) (Supplementary Figure S7), with the TMB correlating negatively with iriscore (R=-0.2, P<0.001) (Supplementary Figure S8). According to the results of survival analysis, we found no statistically significant difference in OS between the high-TMB group and low-TMB group (p>0.05) (Fig. 4D). To further investigate whether iriscore could be a marker to distinguish the survival probability in patients with various levels of TMB, patients were separated into high-iriscore with high TMB, high-iriscore with low TMB, low-iriscore with high TMB and low-iriscore with high TMB. The results indicated that patients with low-iriscore in the low TMB group had a significant survival advantage compared with the other groups; patients with high-iriscore in the high TMB group had the worst prognosis (Fig. 4E). The results above revealed that the constructed iriscore model was closely associated with TMB in CRC patients.
3.5. Iriscore as an efficient predictive tool of responses in chemotherapeutic and targeted therapy
To distinguish the treatment response of various iriscore level in patients with CRC, half the maximum inhibitory concentration (IC50) of anticancer drugs were downloaded from the Genomics of Drug Sensitivity in Cancer (GDSC) website to predict chemotherapeutic and targeted therapeutic responses. Compared with the low-iriscore group, the high-iriscore group had lower levels of IC50 in traditional chemotherapeutic drugs, including docetaxel, camptothecin, doxorubin, and etoposide. (P<0.05) (Fig. 5A-D). Besides, CRC patients with higher iriscore had lower levels of IC50 in targeted drugs such as elesclomol, axitinib, midostaurin, lenalidomide, pazopanib, vorinostat, nilotinib, and temsirolimus (P<0.05) (Fig. 5E-L). On the basis of results above, iriscore could be regarded as an efficient predictor of treatment responses in patients receiving chemotherapeutic and targeted therapy.
3.6. Immune landscapes of low- and high-iriscore CRC patients
The state of microsatellites is closely associated with the efficacy of immunotherapy. However, only approximately 15% of CRC patients are microsatellite instability-high (MSI-H), which is highly sensitive to immunotherapy[20]. The remainder are microsatellite stability (MSS) or microsatellite instability-low (MSI-L), which is resistant to immunotherapy. Hence, it is of great importance to identify MSI-H CRC patients. Here, we investigate the correlation between iriscore and the state of microsatellites. The proportion of patients in the low-iriscore group was 27%, which was much higher than that in the high-iriscore group at 7% (Fig. 6A). The iriscore in MSI-H CRC patients was significantly different from that in MSS or MSI-L CRC patients (P<0.01) (Fig. 6B), and the immune checkpoint TIM-3 was significantly different between the high- and low-iriscore group (P<0.05) (Fig. 6C). We next investigated the relationship between iriscore and immunotherapy. In the IMvigor210 cohort, the CR/PR (response) group had a higher iriscore than the SD/PD (non-response) group (P=0.0009) (Fig. 6D). In addition, differences in immunophenoscore (IPS) when using cytotoxic T lymphocyte-associated antigen-4 (CTLA-4)/programmed death-1 (PD-1) inhibitors between the low-iriscore and high-iriscore groups were assessed to further verify the value of iriscore in immunotherapy (Fig. 6E-H). The IPS of the low-iriscore group using PD-1 inhibitors was higher than that of the high-iriscore group (P<0.05) (Fig. 6F). The low-iriscore group also benefited from CTLA-4 inhibitors (P<0,01) (Fig. 6H). Nonetheless, there was no benefit when using PD-1 and CTLA-4 inhibitors at the same time (P>0.05) (Fig. 6G). Overall, iriscore displayed an obvious correlation with immune signalling pathways and immune cells (Fig. 6I, J). These results show that iriscore has a close association with the immunotherapy response and tumour immune microenvironment.
3.7. Single-RNA sequencing of CRC cells revealed the expression signature of crucial irinotecan-related regulators
We investigated the expression profiles of 10 crucial irinotecan-related regulators in CRC cells through scRNA-seq. As depicted in Fig. 7A and B, 18 cell clusters were identified by unsupervised clustering. Annotation of these clusters revealed the expression signature of the most important irinotecan regulators (Fig. 7C). FSTL3 and TMEM98 had a high expression level in tissue stem cells and fibroblasts in CRC, which may promote cancer progression through regulation of tumour cell proliferation. BCL10, TBC1D10A, CASZ1, CYFIP2, ECHDC1 and LUZP1 were found to be mainly expressed in T cells, macrophages and monocytes; hence, these regulators may influence the effect of immunotherapy by regulating the innate and acquired tumour immune microenvironment. Furthermore, BCL10, CYFIP2, ECHDC1, and LUZP1 were identified in B cells. ANKRD22 was only identified in innate immune cells, such as macrophages, monocytes, neutrophils and dendritic cells. PTP4A2 presented a high expression level in all cell clusters.
3.8. Identification of crucial irinotecan-related regulators and the transcriptional and immune microenvironment landscape
We used random forest analysis to screen the core irinotecan-related regulators further, and the results showed 18 crucial irinotecan-related regulators (Fig. 8A, B). The state of gene alteration, interaction and prognostic features among the 18 irinotecan-related regulators were illustrated in Fig. 8C. FSTL3, ENO2, RIMKLB, PTPRU and CREG2 were the five most important risk factors associated with prognosis in CRC (Fig. 8C), with FSTL3 being the most important. Therefore, FSTL3 was chosen for further analysis. Coexpression of irinotecan-related regulators in the CRC cohort from TCGA was shown (Fig. 8D). We found a positive correlation between FSTL3 expression levels and PTPRU expression levels (R = 0.52, P<0.0001). Considering that irinotecan-related regulators may influence the expression of transcription factors in CRC progression, we utilized a Sankey diagram to investigate the specific relationship between TFs and irinotecan-related regulators (Fig. 8E). The functions of irinotecan-related regulators in CRC may be reflected in the immune microenvironment. As the results showed, immune cell infiltration (Fig. 8F) and immune checkpoints (Fig. 8G) exhibited a close relationship with these irinotecan-related regulators. In addition, FSTL3 presented a strong association with TNFSF4 (R= 0.57, P<0.001), HAVCR2 (R= 0.57, P<0.001) and Macrophage. M2 (R= 0.37, P<0.001) (Supplementary Figure S9-11). Overall, the multiple traits of irinotecan-related regulators indicated that targeting these regulators at transcriptional and immune microenvironment levels might become potential therapeutic approaches for irinotecan resistance in CRC.
3.9. Validation of therapeutic agents targeting irinotecan-related regulators in CRC
We further explored the potential therapeutic drugs (DTP NCI-60) of irinotecan-related regulators and predicted 28 potential irinotecan-related agents, including telatinib (Fig. 9A). Telatinib is an selective small-molecule inhibitor of VEGFR2 VEGFR3 and PDGFR-β tyrosine kinases[21] and is widely used in solid tumour treatment. The three-dimensional structure of telatinib (Fig. 9B) and the molecular docking sites of the telatinib and FSTL3 proteins (Fig. 9C) were exhibited. FSTL3 can form a noncovalent bond with telatinib through asparagine at position 83 and aspartic acid at position 132 (Fig. 9C). The minimum binding energy obtained from the docking of FSTL3 with telatinib was calculated as -7.95 kcal/mol with Autodock. Since telatinib can bind to the FSTL3 protein, it could be served as a target drug in patients with a high level of FSTL3 and might have a potential synergistic effect with irinotecan. Next, potential targets of telatinib were confirmed by using Swiss Target Prediction. We performed the Gene Ontology (GO) and KEGG enrichment analyses of these targets from the database and the results indicated that these genes were mostly associated with cell proliferation, differentiation and metabolism, particularly “regulation of MAPK cascade”, “regulation of kinase activity”, and “PI3K-Akt signalling pathway” (Fig. 9D, E). Overall, therapeutic drugs targeting irinotecan-related regulators, especially telatinib, might act with irinotecan in CRC, as they target protooncogenes such as FSTL3.
3.10. FSTL3 is a core irinotecan-related regulator associated with prognostic and clinical features
K-M survival curves showed that patients with a higher FSTL3 expression had a poorer survival (Fig. 10A). As the ROC curves showed, the AUCs at 1, 3, 5 year was 0.718, 0.682, 0.674 respectively, which indicated the superior predictive ability of FSTL3 expression on prognosis in CRC (Fig. 10B). The distribution of the relative expression level of FSTL3 in CRC patients in the cohort from TCGA is provided in Fig. 9C. The expression level of FSTL3 was significantly higher in CRC patients who did not survive than in those who did survive (Fig. 10C). In addition, the expression level of FSTL3 was closely associated with TNM stage, survival status, and T stage (P<0.05) but not age or sex (P>0.05) in CRC (Fig. 10D-H). The predictive value of FSTL3 had a close relationship with TNM stage (P<0.05) but not age, sex, or T stage in CRC (P>0.05). In general, high expression levels of FSTL3 often predicted poor survival time in patients with T stage III-IV; however, in patients with T stage I-II, there was no significance of FSTL3 expression in the prognosis of CRC (Fig. 10L). Regardless of age or sex, the patients with high expression of FSTL3 had poor prognosis compared to the low expression group (Fig. 10I, J). Interestingly, there was no significant difference between the high- and low-FSTL3 groups with regard to TNM stage I-II or stage III-IV (Fig. 10K).
3.11. Construction and validation of a nomogram integrating FSTL3 and other independent prognostic factors
To determine whether FSTL3 could serve as an independent prognostic predictor, we screened the clinical characteristics (including age, sex, T stage, N stage, M stage and TNM stage) of patients in the CRC cohort using univariate and multivariate Cox regression analyses. The results revealed that age (P<0.05, HR=2.986), T stage (P<0.05, HR=2.113), M stage (P<0.05, HR=2.902) and FSTL3 expression (P<0.05, HR=1.365) were independent predictive factors in CRC patients (Fig. 11A). According to these results, we constructed a nomogram using these independent predictive factors to predict survival probability in CRC patients at 1, 3, and 5 years quantitively (Fig. 11B). The calibration curve was utilized to verify the reliability of the nomogram for predicting the OS probability at 1, 3, and 5 years of CRC patients. The slope of the calibration curve is close to 45°,which showed the reliability of the established nomogram (Fig. 11C-E). Next, we used DCA to evaluate the predictive performance of the clinical model and nomogram for OS time at 1, 3, and 5 years, and DCA showed a net benefit of FSTL3 expression on the prediction of survival possibilities compared to the clinical model including independent prognostic factors in CRC (Fig. 11F-H). Hence, FSTL3 expression level could be considered as an independent prognostic factor and the nomogram we established has a superior performance of prediction on OS.
3.12. FSTL3 might contribute to cancer progression by inhibiting ferroptosis
Ferroptosis is a unique cell death process that is mechanistically and morphologically different from other forms and plays an indispensable part in colorectal cancer progression. To study the mechanism by which FSTL3 promotes cancer progression at the level of cell death, we analysed correlations of FSTL3 and ferroptosis-related crucial regulators and metabolic pathways. The ferroptosis score of each sample in the CRC cohort from TCGA was calculated by utilizing the GSVA algorithm as described above. The ferroptosis score of the low-FSTL3 group was higher than that of the high-FSTL3 group (P<0.05), and there was a significant correlation between FSTL3 and ferroptosis (R=-0.13, P<0.05) (Fig. 12A, B). Next, we confirmed the close correlation between FSTL3 and ferroptosis-related hub genes, and the hub genes GPX4, HSPB1, PTGS2, FTH1, TFRC, SLC40A1, RGS4, ACSL4, TF and DHODH were statistically different between the low- and high-FSTL3 groups (P<0.05) (Fig. 12C). We further found that the FSTL3 expression level was closely associated with ferroptosis markers, such as RGS4 (R=0.59, P<0.001) and HSPB1 (R=0.31, P<0.001) (Fig. 12D-G). FSTL3 also showed a significant correlation with the ferroptosis-related metabolic pathway (Fig. 12H). Importantly, some crucial ferroptosis-related metabolic pathways, including lipid oxidation, regulation of lipid metabolism, reactive oxygen species (ROS) metabolism, phosphatidic acid metabolism, triglyceride acid metabolism, fatty acid metabolism, glutathione metabolism and glutathione peroxidase activity, were significantly different between the high- and low-FSTL3 groups (P<0.05).
3.13. Inhibition of FSTL3 and crucial ferroptosis defence proteins in CRC cells treated with SN38 and ferroptosis induction by irinotecan
At the transcriptional level, we identified relative mRNA expression of the 10 most important regulators in the two CRC cell lines treated with SN38 (Fig. 13A, B). The level of mRNA expression of FSTL3 in both LS180 and HCT116 cells was significantly inhibited after treatment with SN38 (P<0.05). To confirm the effect of ferroptosis induction by SN38, intracellular ROS were detected by flow cytometry analysis in HCT116 and LS180 cells treated with SN38 after 6, 12, 24, 48h, and the ROS level increased with the processing time in both cell lines (Fig. 13C, E). Moreover, SN38 induced lipid peroxidation in HCT116 and LS180 cells (Fig. 13D, F). Western blotting indicated that crucial ferroptosis defence proteins, including SLC7A11 and GPX4, were obviously inhibited in both cell lines treated with SN38, and an increasing inhibitory effect was observed with increasing concentrations of SN38 (Fig. 13G).