Multi-platform-based characterization of ferroptosis in human colorectal cancer.

Ferroptosis is a type of programmed cell death potentially playing an important role in colorectal cancer (CRC) development. However, comprehensive investigations toward ferroptosis in human CRC are lacking. Here, we performed multiple investigations on cancer and para-cancer tissues. We demonstrated that the changes of structural variation and chromatin accessibility in CRC were more associated with the altered mRNA expression of ferroptosis-related genes (FRGs), and the expression of CDKN2A, GPX4, ALOXE3, and LINC00336 was related to the overall survival rates. Subsequently, we revealed that CYBB and YAP1 were potentially the hub genes, and that HSF1 and STAT2 were potentially FRGs' upstream transcription factors. Finally, we depicted the crosstalk between ferroptosis and necrosis, autophagy, and apoptosis. Based on multi-dimensional analyses, we characterized ferroptosis, probable core genes, and the upstream regulators in human CRC. The findings here may improve our understanding of ferroptosis in CRC and provide new opportunities for clinical diagnosis and treatment.

Introduction

According to a report in 2020 from the International Agency for Research on Cancer, colorectal cancer (CRC) kills almost 930,000 people worldwide, making it the world’s third most deadly cancer after breast and lung. Due to the availability of early diagnosis and novel approaches of treatment, there is a significant improvement in incidence and death rates of CRC in recent years (Piawah and Venook, 2019; Ni et al., 2014). However, CRC still causes almost one million deaths a year globally and patients with advanced tumors frequently fail to get effective remission (International Agency for Research on Cancer, 2021).

Ferroptosis is a new form of programmed cell death discovered in 2012 (Dixon et al., 2012). Distinct from necrosis, autophagy, and apoptosis, ferroptosis is frequently caused by excessive lipid peroxidation (Xu et al., 2019; Chen et al., 2019; Bebber et al., 2020), with features of iron accumulation, plasma membrane deficiency, cytoplasm and organelles swelling, chromatin condensation, and mitochondrial disorders (Tang et al., 2021; Vabulas, 2021; Friedmann Angeli et al., 2014). Recently, some studies have revealed that ferroptosis is closely associated with the initiation and progression of CRC (Nie et al., 2018; Li et al., 2020). Multi-omics studies may reveal the molecular changes of ferroptosis at distinct dimensionalities, ranging from DNA, RNA, and protein to epigenetic modification, and improve our understanding of ferroptosis in human CRC. However, such studies are still lacking.

In this study, we performed multi-platform-based investigations on structural variants, DNA methylation, chromatin accessibility, proteome, and phosphoproteome of patients with CRC using our original datasets (n = 6–8). Integrated with the RNA-Seq datasets from The Cancer Genome Atlas (TCGA) (n = 538), we characterized the alterations of ferroptosis in patients with CRC at various molecular levels, from DNA, RNA, and protein to epigenetic modification. We depicted the network of ferroptosis and searched for the hub genes using the proteomic and phosphoproteomic datasets. Moreover, we explored the upstream transcription factors (TFs) potentially regulating the ferroptosis-related genes (FRGs). Finally, we investigated the crosstalk between ferroptosis and another three types of programmed cell deaths, including necrosis, autophagy, and apoptosis.

Results

The proteogenomic characterization of ferroptosis of patients with CRC

To understand the molecular changes of ferroptosis in human CRC, we primarily collected the cancerous and normal adjacent intestinal epithelium tissues from eight patients with CRC, and pooled two samples from the same stage (Table 1). Subsequently, we analyzed the molecular alterations of ferroptosis at multiple levels using our original datasets, including structural variation (n = 6), DNA methylation (n = 6), chromatin accessibility (n = 6), proteome (n = 8), and phosphoproteome (n = 8). We also analyzed the changes of ferroptosis at mRNA level using transcriptome datasets of healthy people (n = 51) and patients with CRC (n = 538) downloaded from TCGA (Figure 1A, Table 2). 159 known FRGs were collected from GeneCards and NCBI (Table S1).

Table 1

The clinical information of eight patients with CRC in our study

	I-II		III-IV
Age
≤65	0	(0.0%)	1	(12.5%)
>65	4	(50.0%)	3	(37.5%)
Types
Adenocarcinoma	4	(50.0%)	4	(50.0%)
Genders
Male	1	(12.5%)	3	(37.5%)
Female	3	(37.5%)	1	(12.5%)
T stages
T1	1	(12.5%)	0	(0.0%)
T2	0	(0.0%)	0	(0.0%)
T3	2	(25.0%)	1	(12.5%)
T4	1	(12.5%)	3	(37.5%)
N stages
N0	4	(50.0%)	0	(0.0%)
N1	0	(0.0%)	4	(50.0%)
N2	0	(0.0%)	0	(0.0%)
M stages
M0	4	(50.0%)	1	(12.5%)
M1	0	(0.0%)	3	(37.5%)
Locations
Ascending colon	2	(25.0%)	1	(12.5%)
Sigmoid colon	1	(12.5%)	3	(37.5%)
Appendix colon	1	(12.5%)	0	(0.0%)
Therapies
Yes	0	(0.0%)	0	(0.0%)
No	4	(50.0%)	4	(50.0%)

Figure 1

The proteogenomic alterations of FRGs in patients with CRC

(A) Schematic showing the overview of this study.

(B) The differential SVs, (C) DNA methylations, (D) peaks of accessible chromatin, (E) mRNAs, (F) proteins, and (G) phosphorylations of FRGs in CRC. Venn analyses of (H) the different SVs, (I) DNA methylations, (J) accessible chromatin, and the differentially expressed mRNAs of FRGs in CRC. Venn analyses of (K) the different SVs, (L) DNA methylations, (M) accessible chromatin, and the survival-relevant FRGs in CRC, related to Tables S2 and S3.

Table 2

The clinical information of 538 patients with CRC from TCGA

	Type	Patients
Fustat	Alive	431 (80.1%)
	Dead	107 (19.9%)
Ages	≤65	232 (43.1%)
	>65	306 (56.9%)
Genders	Male	254 (47.2%)
	Female	284 (52.8%)
Stages	I	95 (17.7%)
	II	208 (38.7%)
	III	149 (27.6%)
	IV	86 (16.0%)
T stages	T1	17 (3.2%)
	T2	92 (17.1%)
	T3	372 (69.1%)
	T4	57 (10.6%)
M stages	M0	452 (84.0%)
	M1	86 (16.0%)
N stages	N0	313 (58.2%)
	N1	128 (23.8%)
	N2	97 (18.0%)

Consequently, 18 FRGs had structural variation in non-exon regions and nine had structural variation in exon regions (Figure 1B). The promoters of 21 FRGs were hyper-methylated and 44 FRGs were hypo-methylated (Figure 1C). Chromatin accessibility of 69 FRGs was increased and 13 FRGs was decreased (Figure 1D). The mRNAs of 65 FRGs were upregulated and 25 FRGs were downregulated (Figure 1E). The proteins of 11 FRGs were upregulated and four FRGs were downregulated (Figure 1F). The phosphorylation of one FRG was upregulated and three FRGs were downregulated (Figure 1G).

Since the transcription of mRNA might be affected by structural variation, DNA methylation, and chromatin accessibility, we screened genes with differential mRNA expression that might be caused by altered structural variation, DNA methylation, and chromatin accessibility. As shown in Figure 1E, we used the transcriptomic dataset downloaded from TCGA for analysis and found that 90 FRGS were differentially expressed in CRC tissues versus normal tissues. As a result, 17 FRGs with structural variants were differentially expressed (68%) (Figure 1H), 26 FRGs with different methylation were differentially expressed (40%) (Figure 1I), and 47 FRGs with distinct chromatin accessibility were differentially expressed (64.4%) (Figure 1J). These results indicated that structural variation and chromatin accessibility potentially affected FRG expression in human CRC the most.

Moreover, we searched for survival-related FRGs, which also had genomic or epigenetic changes in patients with CRC. Firstly, we obtained 1,598 survival-relevant genes through analyses of RNA-Seq datasets from TCGA (Table S2). Subsequently, 15 FRGs were identified as survival-related genes (Table S3). Among all these 15 FRGs, three genes harbored structural variants (20.0%) (Figure 1K) and seven genes had hyper or hypo DNA methylation (46.7%) (Figure 1L). In addition, 10 genes were found to have altered chromatin accessibility (66.7%) (Figure 1M). This result indicated that the altered chromatin accessibility of FRGs in patients with CRC might be more closely associated with patient survival.

The transcriptomic analysis revealed a significant change of ferroptosis in patients with CRC and uncovered an association between ferroptosis and apoptosis and autophagy

Nowadays, although ferroptosis has been demonstrated to be related to the development of CRC (Sui et al., 2018; Xia et al., 2020), most functional FRGs in CRC are still not clear, as are their interaction networks. Therefore, we analyzed the FRGs mRNAs using the transcriptomic datasets downloaded from TCGA. As a result, 90 FRGs were differentially expressed in CRC (56.6% of all known FRGs) (Figure 2A), suggesting a great change of ferroptosis in human CRC. Next, we analyzed the survival-related genes in all differentially expressed FRGs using the RNA-Seq datasets of patients with CRC downloaded from TCGA and uncovered that the mRNA expression of four FRGs was survival-associated, including CDKN2A, GPX4, ALOXE3, and LINC00336 (Figures 2B and 2C). We also used databases from the Gene Expression Omnibus (GEO) (GSE33113) to validate the mRNA expression of CDKN2A, GPX4, and ALOXE3 in 90 CRC samples and 6 healthy controls (Figure 2D). We checked CDKN2A protein expression in The Human Protein Atlas (HPA) and the Clinical Proteomic Tumor Analysis Consortium database (CPTAC), and discovered that CDKN2A was overexpressed in CRC tissues versus normal tissues (Figures 2E and 2F). Furthermore, we performed enrichment analysis on these differentially expressed FRGs and found some noteworthy pathways, such as negative regulation of extrinsic apoptotic signaling pathway, autolysosome, and p53 signaling pathway (Figures 2G–2J). These results suggested that the functions of ferroptosis were potentially relevant to apoptotic and autophagic pathways.

Figure 2

The transcriptome analysis demonstrated that the abnormal ferroptosis in CRC was associated with patient survival, and revealed a potentially regulatory role in apoptosis and autophagy

(A) Volcano plot of differentially expressed FRGs in CRC.

(B and C) The (B) survival analysis of four genes in patients with CRC and (C) their mRNA expression. “N” represented normal adjacent tissues and “T” represented tumor tissue. Statistical significance is calculated with Student’s t test (∗p < 0.05).

(D) The mRNA expression of three survival genes in 90 CRC samples and 6 normal controls as analyzed using the dataset from GEO (GSE33113). Data are represented as mean ± SD. Statistical significance is calculated with Student’s t test (∗p < 0.05, ∗∗p < 0.01).

(E and F) The (E) protein expression and (F) immunohistochemical detection of CDKN2A in normal and colorectal adenocarcinoma tissues. The immunohistochemical results were obtained from HPA. “Normal” represented normal adjacent tissues and “CRC” represented CRC tissues.

(G–I) The GO and (J) KEGG analysis of differentially expressed FRGs in CRC.

Proteomic investigations uncovered the network of ferroptosis in CRC and disclosed CYBB potentially as a hub gene regulating ferroptosis

Proteins are the direct units that carry out functions in cells. Subsequently, we analyzed ferroptosis in patients with CRC on views of proteome. After overlapping the 159 FRGs with differentially expressed proteins in CRC tissues versus normal adjacent tissues, we found that the expression of 15 FRGs was altered (Figure 3A). To explore the network of ferroptosis in human CRC, we disclosed the linking of these 15 differentially expressed proteins and their first neighbors using Cytoscape (Figure 3B). In order to further understand how these genes functioned in tumorigenesis, we performed enrichment analysis of the 15 FRGs and their first neighbors using Metascape database. The results showed that these genes were potentially involved in neutrophil degranulation, response to oxidative stress, reactive oxygen species metabolic process, detoxification of reactive oxygen species, etc (Figure 3C). This result suggested that ferroptosis in CRC might affect tumor immunity and oxidative stress. In addition, we searched for the top ten hub genes within the network and depicted their interaction using cytoHubba. As a result, ten high-scoring genes were selected, including CYBB, ITGAM, ITGB2, CYBA, ADAM10, PLD1, SLC2A3, ITGAV, PTPRJ, and SLC44A2 (Figures 3D and 3E). Interestingly, CYBB, which was one FRG, gained the highest score, suggesting that this gene was likely to play a core role in ferroptosis network. Therefore, we then focused on CYBB, exploring its expression and potential functions in CRC. We used data from TCGA to confirm the mRNA expression of CYBB in 620 CRC samples and 51 normal controls. Beyond that, we verified the increased protein expression of CYBB using the dataset from CPTAC (ntumor = 97, nnormal = 100) (p < 0.001) (Figures 3F and 3G). We further obtained the co-expressed proteins of CYBB from LinkedOmics and performed enrichment analysis to gain the insight of CYBB biological meaning in CRC. The top 50 positively correlated and negatively correlated proteins were shown in Figure S1. The KEGG analysis showed that the pathways related to complement and coagulation cascades, leukocyte transendothelial migration, and systemic lupus erythematosus were enriched (Figure 3H), indicating that the functions of CYBB might be highly associated with tumor immunity in CRC. Therefore, we examined the correlation between CYBB expression and immune score, stroma score, microenvironment score, and infiltration of immune cells using the datasets from Xcell database. As a result, the expression of CYBB was significantly associated with the immune score, stroma score, and microenvironment score, as well as the infiltration of many immune cells, such as neutrophils, natural killer cells, myeloid dendritic cells, macrophages, and B cells (Figure 3I). These results proclaimed a strong correlation between CYBB expression and the immune microenvironment, and suggested a possible relationship between ferroptosis and tumor immunity microenvironment in CRC.

Figure 3

Proteomic investigations proclaimed the network of ferroptosis in CRC and disclosed CYBB as a potentially hub gene in ferroptosis

(A) The heatmap showing protein levels of the 15 FRGs in normal adjacent tissues and CRC tissues. “N” represented normal adjacent tissues, and “C” represented CRC tissues.

(B) The network between FRGs and their first neighbors at protein level.

(C) The functional cluster network of FRGs and their first neighbors at protein level.

(D and E) The (D) scores and (E) network of the hub FRGs at protein level.

(F and G) The (F) mRNA expression and (G) protein expression of CYBB in CRC tissues and normal tissues. “Normal” represented normal adjacent tissues and “CRC” represented CRC tissues. Statistical significance is calculated with Student’s t test (∗∗∗p < 0.001), related to Figure S1.

(H) KEGG enrichment of CYBB co-expression proteins.

(I) Immune correlation analysis of CYBB.

Phosphoproteomic analysis revealed a potentially important role of decreased phosphorylation of YAP1 at Ser109 in ferroptosis network

Subsequently, we analyzed the phosphoproteomic changes of the FRGs in patients with CRC. As a result, the phosphorylation of four FRGs was altered (Figure 4A). We then described the network between the four FRGs and their connected first neighbors using Cytoscape (Figure 4B). Next, we performed enrichment analysis of these genes using Metascape database to further understand the biological processes that these genes participated in. As shown in Figure 4C, these genes were mainly associated with processes such as epithelial cell development, signaling by Hippo, and apoptotic signaling pathway. Additionally, we analyzed the top ten hub genes in the network using cytoHubba, and found that YAP1, EGFR, TJP1, MSN, SLC9A3R1, CTNND1, CREB1, TJP2, IQGAP1, and SMARCA4 were screened (Figures 4D and 4E). Remarkably, the gene gained the highest score, YAP1, was also one known FRG. In our phosphoproteomic examination, we discovered that the phosphorylation of YAP1 at Ser109 was downregulated in patients with CRC(Figure 4A). Besides, we found that YAP1 and HSPB1 were both altered at the protein and phosphorylation levels (Figure 4F). These results suggested that YAP1 Ser109 might play an important role in ferroptosis. Therefore, we verified the phosphorylation level of YAP1 at Ser109 in 97 patients with CRC using the datasets from UALCAN database. As shown in Figure 4G, the phosphorylation of YAP1 at Ser109 was reduced in cancerous tissues compared to normal tissues. Additionally, to investigate the function of YAP1 phosphorylation in CRC, we further explored the phosphorylation co-expressed with YAP1 Ser109. The top 50 phosphorylations positively or negatively correlated with YAP1 Ser109 were shown in the heatmap (Figure S2). We analyzed all the expression-correlated phosphorylation of YAP1 Ser109 using Gene Set Enrichment Analysis, and discovered that YAP1 Ser109 potentially participated in processes named actin cytoskeleton reorganization and G0 to G1 transition in CRC (Figure 4H).

Figure 4

Phosphoproteomic analysis revealed a potentially important role of YAP1 S109 in ferroptosis network

(A) The heatmap showing phosphorylation level of the four FRGs in normal adjacent tissues and CRC tissues. “N” represented normal adjacent tissues, and “C” represented CRC tissues.

(B) The network between FRGs and their first neighbors at phosphorylation level.

(C) The functional cluster network of FRGs and their first neighbors at phosphorylation level.

(D–F) The (D) scores and (E) network of the hub FRGs at phosphorylation level. Venn analyses of (F) the differentially expressed proteins and phosphorylations of FRGs in CRC.

(G) The phosphorylation expression of YAP1 at Ser109 in CRC tissues versus normal tissues. “Normal” represented normal adjacent tissues and “CRC” represented CRC tissues, related to Figure S2.

(H) KEGG enrichment of co-expressed phosphorylation of YAP1 Ser109 in CRC.

A multilayer investigation identified HSF1 and STAT2 as the potentially upstream TFs regulating ferroptosis in patients with CRC

As is well known, transcription factors regulate the expression of target genes. Therefore, we queried 90 differentially expressed FRGs using hTFtarget database and obtained 162 experimentally confirmed TFs of the FRGs in CRC tissues (Table S4). Transcription factors control gene transcription by recognizing binding motifs on DNA, thus we firstly searched for the DNA-binding motifs of the 162 TFs using our original assay for transposase-accessible chromatin with high-throughput sequencing datasets. As a result, the binding motifs of SP1, SP2, SP3, ZBTB7B, KLF5, GLIS1, ESRRA, NKX2-5, and USF2 were only opened in para-cancerous tissues, while the binding motifs of JUN, SOX4, STAT2, SPIB, JUND, GATA4, HSF1, ASCL2, KLF3, REST, GATA2, and ELF1 were specifically opened in cancerous tissues (Figures 5A and 5B). Among them, SP1 and JUND potentially regulated 95.6% and 56.7% of all the differentially expressed FRGs (Figures 5C and 5D). Next, we analyzed the protein expression and phosphorylation of the 162 TFs in CRC tissues versus normal tissues using our original data. The results showed that five TFs that possibly targeting FRGs in CRC were changed at the protein level, including YY1, CDX2, HSF1, NFKB2, and STAT2 (Figure 5E), and two TFs were changed at the phosphorylation level, including CBX3 and CREB1 (Figure 5F). Noticeably, YY1 potentially regulated more than 50% of expression-altered FRGs (Figure 5G). CBX3 potentially targeted more than 10% of differentially expressed FRGs in CRC (Figure 5H).

Figure 5

Multi-platform-based analyses disclosed HSF1 and STAT2 as the possibly upstream TFs regulating ferroptosis network in patients with CRC

(A and B) The binding motifs of the TFs probably regulated the differentially expressed FRGs in (A) normal and (B) CRC tissues, related to Table S4.

(C) The proportion of the genes that potentially targeted by the TFs selected in (A) in all differentially expressed FRGs.

(D–H) The (D) proportion of the genes that potentially targeted by the TFs selected in (B) in all differentially expressed FRGs. The potential FRG TFs were altered on (E) the protein and (F) phosphorylation level, and (G–H) the proportion of their potentially targeting expression-changed FRGs in all differentially expressed FRGs. Data are represented as mean ± SD. Statistical significance is calculated with Student’s t test (∗p < 0.05, ∗∗p < 0.01), related to Table S4.

(I) The Venn plot showed the intersection of TFs obtained from (B) and (E), related to Figure S3 and S4.

(J–N) The (J) mRNA expression and (K–L) protein expression of STAT2 and HSF1, and (M−N) the immunohistochemical detection of these two genes in normal and CRC tissues. The immunohistochemical results were obtained from HPA database. “Normal” represented normal adjacent tissues and “CRC” represented CRC tissues.

We then overlapped all the TFs screened in the above analysis, and found that the TFs HSF1 and STAT2 were selected (Figure 5I). In our proteomic investigation, the protein expression of HSF1 and STAT2 was both upregulated in CRC. Meanwhile, the binding motifs of HSF1 and STAT2 on DNA were also more opened in CRC tissues than in normal tissues. This result suggested that these two TFs were more likely to be functional TFs in CRC. Subsequently, to clarify the FGRs targeted by these two TFs, we queried the target genes of HSF1 and STAT2 using hTFtarget database. As a result, the FRGs HMGB1 and TP53 were found to be targeted by HSF1 and STAT2, respectively. Furthermore, we discovered that the mRNA levels and protein levels of HMGB1 and TP53 were both upregulated in CRC tissues compared to normal tissues, further indicating that HMGB1 and TP53 were potentially downstream genes of HSF1 and STAT2 in CRC cells (Figures S3 and S4). Meanwhile, using TCGA data, we found that HSF1 and STAT2 mRNA levels were also higher in CRC tissues than in normal tissues (ntumor = 620, nnormal = 51) (Figure 5J). Moreover, we analyzed their protein expression in 97 patients with CRC and 100 healthy people using datasets from CPTAC database. As a result, HSF1 and STAT2 exhibited increased expression in CRC tissues versus normal tissues (Figures 5K–5L). Additionally, we evaluated the expression of HSF1 and STAT2 in CRC tissues and normal tissues by immunohistochemical detection from HPA database and also found high expression of the two proteins in CRC (Figures 5M−5N).

The crosstalk between ferroptosis and apoptosis, necrosis, and autophagy

Some studies have proclaimed that there is probable crosstalk between ferroptosis and necrosis, autophagy, and apoptosis in CRC (Guo et al., 2018; Hong et al., 2017). However, little is known about the genes or pathways that connect ferroptosis with other programmed deaths. Therefore, we primarily downloaded genes involved in necrosis, autophagy, and apoptosis from the Harmonizome database (Tables S5, S6, and S7) and then investigated the protein changes of the genes involved in these four types of cell deaths in CRC. As a consequence, the protein expression of 22 necrosis-related genes was upregulated and 13 was downregulated in CRC tissues versus para-cancerous tissues; the protein expression of 25 autophagy-related genes was upregulated and nine was downregulated in CRC tissues compared to para-cancerous tissues; the protein expression of 129 apoptosis-related genes was upregulated and 79 was downregulated in CRC tissues versus para-cancerous tissues; the protein expression of 11 FRGs was upregulated and four proteins was downregulated in CRC tissues versus normal tissues (Figure 6A). Subsequently, we uncovered the phosphorylation changes of the proteins of the four types of cell deaths in CRC. The results revealed that one necrosis-related protein had increased phosphorylation and 11 proteins had decreased phosphorylation in CRC tissues versus para-cancerous tissues; six autophagy-related proteins had increased phosphorylation and 24 proteins had decreased phosphorylation in CRC tissues versus para-cancerous tissues; 40 apoptosis-related proteins had increased phosphorylation and 78 proteins had decreased phosphorylation in CRC tissues versus para-cancerous tissues; one ferroptosis-related protein had increased phosphorylation and three had decreased phosphorylation in CRC tissues versus para-cancerous tissues (Figure 6B). Notably, after overlapping all the differentially expressed proteins and phosphorylation selected above, we discovered that HSPB1 and RIPK1 participated in all these four cell deaths simultaneously, indicating that they might be the hub genes linking these four types of cell programmed deaths together (Figures 6C and 6D). Then, we analyzed the protein and phosphorylation levels of HSPB1 and RIPK1 using our own proteomic and phosphoproteomic datasets, and found that the protein content of HSPB1 was higher in CRC samples versus normal controls (Figure 6E), and the phosphorylation levels of HSPB1 and RIPK1 were both downregulated in CRC tissues compared to para-cancerous tissues (Figures 6F and 6G). Subsequently, to further understand the crosstalk between the four types of deaths in CRC, we constructed networks between each two types of cell deaths using the differentially expressed proteins and phosphorylations. Shown in Figure 6H, ferroptosis linked to necrosis through the proteins TFRC, RIPK1, and HSPB1; ferroptosis linked to autophagy through RIPK1, HSPB1, HSF1, YAP1, and VDAC2 (Figure 6I); ferroptosis linked to apoptosis through FH, GOT1, C12orf5, TFRC, HSPB1, and CYBB (Figures 6J and 6K). To our best knowledge, this is the first complete picture that shows the connectivity between ferroptosis and the other three programmed deaths in patients with CRC.

Figure 6

The crosstalk between ferroptosis and apoptosis, necrosis, and autophagy

(A–D) The differentially expressed (A) proteins, and (B) phosphorylations of necrosis, autophagy, apoptosis and ferroptosis in CRC patients, related to Tables S5, S6, and S7. Venn analyses of the differentially expressed (C) proteins and (D) phosphorylations of necrosis, autophagy, apoptosis and ferroptosis, related to Tables S5, S6, and S7.

(E–K) The (E–G) protein or phosphorylation levels of HSPB1 and RIPK1 in normal and CRC tissues. Data are represented as mean ± SD (H) ferroptosis and necrosis, (I) ferroptosis and autophagy, (J) ferroptosis and apoptosis, and (K) ferroptosis and apoptosis networks based on the differentially expressed proteins and phosphorylations.

Discussion

In this study, we described the molecular changes of ferroptosis in CRC at multiple levels, including structural variation, DNA methylation, chromatin accessibility, transcriptome, proteome, and phosphoproteome. We revealed the survival-related mRNAs, hub proteins, and phosphorylation within the ferroptosis network, and disclosed the potential upstream TFs of ferroptosis using multi-omics analysis. In addition, we delineated the crosstalk between ferroptosis and other three programmed cell deaths in patients with CRC, including necrosis, autophagy, and apoptosis. Taken together, our study characterized ferroptosis network in patients with CRC at various molecular levels and provided novel potential drug targets for CRC treatment.

Ferroptosis is a new form of cell death caused by iron-dependent lipid peroxidation, the abnormity of which probably leads to various human diseases (Wu et al., 2019; Dixon et al., 2012). On the one hand, ferroptosis is considered to suppress tumor development, whose activation may kill cancer cells (Badgley et al., 2020; Lei et al., 2020). On the other hand, ferroptosis has also been found to promote tumor growth (Dai et al., 2020). So far, more than 150 genes have been found participating in ferroptosis, including TP53, GPX4, SLC7A11, ACSL4, NOX, and NRF2 (Xu et al., 2021; Kang et al., 2019). In our study, we discovered that the mRNA expression of TP53 was elevated in CRC tissues versus normal tissues potentially regulated via STAT2. In previous studies of CRC, TP53 has been demonstrated to antagonize ferroptosis by forming the DPP4-p53 complex (Zhang et al., 2018). Therefore, the high expression of STAT2 and TP53 may be associated with inhibition of ferroptosis in CRC, thus promoting tumor initiation and progression. In conclusion, genes involved in ferroptosis may provide new therapeutic strategies for the treatment of CRC in clinic.

In view of the high recurrence rate and individual heterogeneity of patients with CRC, prognostic biomarkers can be used to predict the high-risk patients to help patients obtain appropriate treatment (Ahluwalia et al., 2019; Sveen et al., 2020). However, there are still no molecular prognostic biomarkers applied in the clinical diagnosis of CRC. In our study, we found four prognosis-relevant FRGs, including CDKN2A, GPX4, ALOXE3, and LINC00336. In future studies, we may expand the sample to verify the values of these four genes as prognosis factors and co-use the four genes as a signature, testing their efficacy in the prediction of CRC prognosis.

Nowadays, several key genes involved in programmed cell death have been developed as therapeutic targets for cancer treatment, such as BCL2L1, BCL2, BCL2L2, MCL1, and BCL2L10 (Barman et al., 2018). This indicates that programmed cell death may be a promising target for tumor therapy. In this study, we uncovered several genes as core regulators of ferroptosis network, including CYBB, YAP1, HSF1, STAT2, HSPB1, and RIPK1. These genes have potentials as drug targets functioned on ferroptosis. For example, some studies have shown that the protein expression of HSPB1 (also called human HSP27) was increased in CRC (Hung et al., 2020; Yu et al., 2010). There were considerable evidences that the increased expression of HSP27 in cells enhances the tumorigenicity of CRC (Lee et al., 2018). HSPB1 has been identified as a negative regulator of ferroptosis in various cancer cell lines, including HeLa, U2OS, and LNCaP (Sun et al., 2015). Consistent with the previous study, we found that HSPB1 was higher expressed in CRC tissues versus normal tissues. Therefore, the elevated expression of HSPB1 in CRC may inhibit ferroptosis, or even other types of cell deaths, thus promoting tumor initiation and progression. To sum up, HSPB1 and other hub FRGs highlighted in our study may be used as drug targets for CRC treatment.

In conclusion, our study provides a comprehensive molecular landscape of ferroptosis in CRC, uncovers key FRGs from distinct dimensionalities, and discloses potential regulators of the ferroptosis network. The findings in our study may provide a new perspective for understanding ferroptosis in tumorigenesis and lay a theoretical basis for diagnosis and treatment of CRC.

Limitations of the study

The number of samples used for omics analysis is relatively small, with only six to eight patients of colorectal adenocarcinoma. To some extent, the data may lack universality. High-throughput findings based on large samples or small samples may not be identical. Through validation by large-scale samples from public databases, we preliminarily confirmed our findings. However, the validation by the authors' own cohort was still absent. The discovered molecular changes may also vary at different levels in different populations due to differences in race, age, region, or pathological typing, etc.

In this study, our goal was to uncover a group of genes with molecular changes that might be involved in CRC’s ferroptosis process using a high-throughput approach. Whereas, the highlighted genes were not taken histological verification or functional analysis. Molecular alterations do not always lead to functional changes; we cannot be certain that the genes discovered through differential expression examination, survival analysis, enrichment investigation, and algorithmic calculations play a role in ferroptosis in colon cancer cells. From this study, only the potential panoramic molecular context for ferroptosis in colon cancer can be inferred. Histopathological or functional verification is still required in the future to clarify the genes’ roles in ferroptosis of colon cancer.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Yong Dai (daiyong22@aliyun.com or dai.yong@szhospital.com).

Materials availability

This study did not generate new unique reagents.

Experimental model and subject details

We obtained eight samples of colorectal adenocarcinoma from Shenzhen People’s Hospital. The samples were used for previous study (Zhang et al., 2021a, 2021b). Seven patients were older than 65, and one was younger than 65 (Four are male and four are female). All patients signed informed consent. This work was approved by the ethics committee of Shenzhen People's Hospital (LL-KY-2019213).

Method details

Protein extraction and trypsin digestion

The cancer tissues and paracancer tissues of eight CRC patients from colon segment were ground into cell powder by liquid nitrogen. Then, the mixture was next crushed by the ultrasonic equipment (Scientz) with lysis buffer [8 M urea and 1% protease inhibitor cocktail III (Merck Millipore, 156535140)]. The preceding steps were performed on ice. After that, the crushing liquid was centrifuged at 12000 g for 10 min at 4°C, and the supernatant was collected to determine the protein concentration using the BCA kit (Beyotime P0011-1) as directed by the manufacturer.

Following that, trichloroacetic acid (TCA) was added to the ultrasonic lysate at a final concentration of 20%. After storage at 4°C for 2 h, the mixture was centrifuged at 4500 g for 5 min at 4°C. Then, the supernatant was removed, and the precipitation was washed with pre-cooled acetone 2-3 times. Then triethylammonium bicarbonate (TEAB) was added to the precipitation, bringing the final concentration to 200 mM. After centrifugation, trypsin was injected in a 1:50 ratio for incubation. After that, dithiothreitol (DTT) was employed to the mixture staying at 56°C for 30 min (final concentration = 5 mM). Finally, the mixture was added with Iodoacetamide (IAA) (final concentration = 11 mM) and sit at room temperature in the dark for 15 min.

Liquid chromatography mass spectrometry (LC-MS)/MS and protein annotation

The LC-MS/MS was carried out according to a standard methodology. The quadrupole time of flight mass spectrometer named timsTOF Pro used in this experiment is equipped a modified nano-electrospray ion source (Bruker Daltonics), and the running mode is parallel accumulation serial fragmentation (PASEF). The 75-m i.d., 15-cm long and 100-m i.d., 25-cm long handmade inverted-phase analytical columns for proteomic and phosphoproteomic investigations, respectively, contain the tested material dissolved in 0.1% formic acid (solvent A). For proteomic examination, the concentration of solvent B (98% acetonitrile and 0.1% formic acid) was increased from 6% to 24% for 70 min. The concentration of solvent B then increased from 24% to 35% for 14 min. After that, the concentration of solvent B reached up to 80% for 3 min and then maintained at 80% for additional 3 min. For phosphoproteomic examination, the concentration of solvent B (98% acetonitrile with 0.1% formic acid) increased from 2% to 22% within 50 min. Then it soared from 22% to 35% in 2 min, and continued to rise to 90% within 3 min and then maintained at that level for another 5 min. All the analyses were performed on the nanoElute UPLC system with a constant flow rate of 450 nL/min. The peptides were transferred to a glass capillary tube at a 1.6 kV electrospray voltage, and then processed by timsTOF Pro mass spectrometry. The precursor and debris were detected and analyzed using TOF, and MS/MS scanning was carried out in the range of 100 m/z-1700 m/z. The precursors with charge states ranging from 0 to 5 were chosen for segmentation, and PASEF-MS/MS scans were collected 10 times per cycle with a dynamic elimination time of 30 s.

The results of the MS/MS data were analyzed using the Maxquant software package (v1.6.6.0) (Cox and Mann, 2008), and Homo_sapiens_9606_SP_201911-15 (20380 sequences) was used as the reference, and a reverse library was included to calculate the false positive rate caused by random matching. Additionally, to mitigate the influence of contaminated proteins, typical polluted libraries were added to the database. The minimum peptide length was set to 7 and the maximum modification peptide quantity was set at 5. In the initial and principal search, precursor ions had a mass tolerance of 20 PPM, while fragment ions had a mass tolerance of 0.02 Da. The fixed modification was set to the cysteine alkylation, and the variable modification was set to methionine oxidation. In phosphoproteomic analysis, carbamidomethyl on Cys was set as the fixed modification, and oxidation of Met, acetylation of N-termini of proteins, and phosphorylation of Ser, Thr, Tyr residues were considered as variable modifications. Finally, the threshold of peptide-spectrum match (PSM) and false discovery rate (FDR) were kept within 1%.

The UniProt-GOA database was used to annotate GO terms using the InterProScan (Zdobnov and Apweiler, 2001). Besides, the identified proteins were annotated using KAAS v.2.0 (KEGG online website) for KEGG enrichment, and KEGG mapper V2.5 was used to match the proteins to KEGG pathways (Kanehisa et al., 2012). Finally, WoLFPSOR v.0.2 software was used to do the subcellular enrichment.

Assay for transposase-accessible chromatin with highthroughput sequencing (ATAC-Seq)

50,000 single cells from the digested tissues were broken down with the lysis buffer (0.1% IGEPAL CA-630, 3 mM MgCl2, 10 mM NaCl, 10 mM Tris-HCl, pH 7.5). Then, the centrifugation was performed at 500 g for 10 min at 4°C to remove the supernatant. The precipitation was incubated in the reaction solution (2.5 μl of Tn5 transposase and 1 × TD buffer in 50 μl total) at 37°C for 30 min. A MinElute Reaction Cleanup Kit (QIAGEN, 51306) was used to purify DNA, and the TruePrep™DNA Library Prep Kit V2 of Illumina (Vazyme Biotech, TD501/TD502/TD503) was utilized to build DNA libraries. The quality of library was checked using the real time PCR tests. The lengths of inserted fragments were determined using the HS 2100 Bioanalyzer (Agilent). Based on 150-bp paired-end reads, the high-quality library was sequenced using the Illumina sequencing platform HiSeq X Ten. The initial results were saved in FASTQ format, which included base sequences and quality control information. Clean reads were recovered after removing low-quality and adaptor-contaminated reads, and Bowtie2 was used to match them to the matching reference genome (Langmead and Salzberg, 2012). The data was visualized using the Integrative Genomics Viewer (IGV), and MACS2 was used to call peaks within sections of the open genome.

Whole genome bisulfite sequencing (WGBS)

After preparing high-quality DNA with the SMRTbell Express Template Prep Kit 2.0 (Pacific Bioscience, 100-938-900), the DNA library was built using 1 g DNA. C-T conversion was achieved using the EZ DNA MethylationGoldTM Kit (ZYMO research, D5005). The preliminary DNA was then quantified using a Qubit 2.0 fluorometer (Qubit), and the range of library insert size was determined using a 2100 Bioanalyzer (Agilent). The library's effective concentration (>10 nM) was then determined using real-time PCR tests, and the verified library was sequenced on the HiSeq X10 sequencing instrument (Illumina) based on paired-end 150-bp reads.

Long-read whole genome sequencing

The DNA was broken into fragments of around 20 kb using a g-TUBETM (Covaris, 520079). Following that, magnetic beads were added (Pacifi Bioscience, 100-317-100). After adding stem-loop junctions to DNA fragments, a fresh round of screening was performed. Next, DNA and primers were annealed, and DNA polymerase was used for DNA extension (Pacifi Bioscience, 101-731-100). The quality of the DNA was then assessed using an Agilent 2100 Bioanalyzer, and paired-end 150-bp reads were generated by the Hiseq X10 sequencing platform (Illumina).

TCGA

The TCGA database was used to obtain the Level 3 RNA-Seq datasets, which included information from 51 healthy people and 647 CRC patients. A total of 538 samples were included in the study after the subjects with uncertain survival duration were removed. Based on the Ensembl database documentation, Ensembl IDs were converted to gene symbols. Following that, the datasets of healthy controls and CRC patients were retrieved from the symbol matrix. The survival analyses were performed using R package “survival” v.2.4.2. The codes for differential gene selection and survival curve drawing have been uploaded as supplemental information (Related to Data S1).

Other databases

The list of FRGs was obtained from Genecards and NCBI (Related to Table S1). The potential TFs for FRGs were collected using hTFtarget. The DNA-binding domains of TFs were obtained from The Human Transcription Factors. The enrichment analysis was performed on Metascape.

Quantification and statistical analysis

The statistic significance of differential events was calculated using T tests. The significance of enrichment results was calculated by fisher’s exact tests. Data management and statistical analysis were performed using GraphPad Prism 7 software. Data are reported as mean +/- standard deviation (SD). Significance was set to p < 0.05 (∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001).

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Chemicals, peptides, and recombinant proteins
urea	Sigma-Aldrich	V900119-500G
protease inhibitor cocktail III	Merck Millipore	156535140
trichloroacetic acid	Sigma-Aldrich	T4885-2KG
acetone	Hangzhou Hannuo	N/A
TEAB	Sigma-Aldrich	T7408-500mL
trypsin	Promega	SRT0201
dithiothreitol	Sigma-Aldrich	D9163-25G
Iodoacetamide	Sigma-Aldrich	V900335-5G
acetonitrile	ThermoFisherScientific	204433
formic acid	Fluka	A117-50
Tn5 transposase	vazyme	TD501-01
TD buffer	vazyme	TD501-01
Critical commercial assays
BCA kit	Beyotime	P0011-1
MinElute Reaction Cleanup Kit	QIAGEN	51306
TruePrep™DNA Library Prep Kit V2	Vazyme Biotech	TD501/TD502/TD503
SMRTbell Express Template Prep Kit 2.0	Pacific Bioscience	100-938-900
EZ DNA MethylationGold™ Kit	ZYMO research	D5005
MBSpure beads	Pacifi Bioscience	100-317-100
DNA polymerase	Pacifi Bioscience	101-731-100
Deposited data
Level 3 RNA-Seq	TCGA	https://portal.gdc.cancer.gov/
Proteome data	This paper	PXD021314
Phosphoproteome data	This paper	PXD021318
ATAC-Seq; WGBS; Long-read whole genomesequencing	This paper	PRJNA693028
Software and algorithms
Maxquant software package v1.6.6.0	Cox and Mann, 2008	https://www.maxquant.org/
InterProScan	Zdobnov and Apweiler, 2001	https://www.ebi.ac.uk/interpro/
KAAS v.2.0	N/A	https://www.genome.jp/tools/kaas/
KEGG mapper V2.5	Kanehisa et al., 2012	https://www.genome.jp/k egg/mapper/
WoLFPSOR v.0.2	N/A	https://wolfpsort.hgc.jp/
R package “survival” v.2.4.2	N/A	https://cran.r-project.org/ web/packages/survival/index.html
Bowtie2	Langmead and Salzberg, 2012	http://bowtie-bio.sourceforge.net/bowtie2/index.sht ml
Integrative Genomics Viewer	N/A	http://software.broadinstitute.org/software/igv/home
MACS2	N/A	https://pypi.org/project/MACS2/
GraphPad Prism 7.00	GraphPad Software	https://www.graphpad.com/