Patient information and Quality control of samples
Among the 24 patients, 11 achieved pCR(54%). The baseline clinical characteristics including sex, age, smoking history, tumor location, tumor clinical T stage and N status were similar between two response groups(table S2).
WES was performed on 18 tumor samples, as six lacked blood control. Seven RNA samples showed significant degradation on electrophoresis, and one sample contained insufficient protein for spectrometry. These samples were not included in the experiment. In total, RNA sequencing and DIA spectrometry were performed on 17 and 23 samples, respectively(table S2).
Interferon signaling and RIG-I-like receptor signaling pathways were enriched in non-pCR tumors in the joint analysis of transcriptomics and proteomics
RNA-seq was performed in 17 of the 24 patients. Genes differentially expressed between pCR and non-pCR patients were analyzed using EdgeR and DESeq2 in R. The 85 genes at the intersection of the two results were identified as differential genes(p <0.005, log2FC>0.58 or <-0.58) (figure S1, table S3). GSEA analysis highlighted that the most significantly enriched pathways in the non-pCR group were the type I interferon signaling pathway and the defense response to viruses (adjusted p-value = 0.000). Conversely, the pCR patients showed significant enrichment in keratinization-related pathways (adjusted p-value < 0.001) (figure 1A and table S4). Furthermore, in non-pCR patients, the MDA-5 signaling pathway, an upstream regulator of interferon signaling, was enriched, as revealed by Gene Ontology (GO) enrichment analysis of the differential genes (adjusted p-value = 0.051).
DIA mass spectrometry was performed in 23 of the 24 patients, identifying 4080 proteins. Analysis using PCA and unsupervised hierarchical clustering with valid proteins revealed distinct protein expression patterns in patients with distinct treatment responses(figure 1S). Combined limma and t-test identified 385 differentially expressed proteins(p <0.05, log2FC>0.58 or <-0.58) between pCR and non-pCR patients(figure 1S). Gene ontology enrichment analysis was performed based on the differential proteins, and the top 30 enriched pathways were all derived from the non-pCR group(figure 1B). The top enriched pathways included RNA catabolic-related and neutrophil activation pathways.
In accordance with the transcriptomic enrichment results, the type I interferon signaling pathway was also among the significant pathways enriched in the non-pCR group(adjusted p-value =0.003). In addition, the dsRNA response pathway and Wnt signaling pathway enriched in proteomics corresponded with the MDA-5 signaling pathway and Wnt protein secretion pathway in transcriptomics, respectively (figure 1C)
To identify the consistently enriched pathways across different gene expression omics while also detecting pathways not apparent in any single-omic dataset, we integrated the transcriptomics and the proteomics data using ActivePathways package in R.11 In genes elevated in the non-pCR group, the joint analysis of ActivePathways highlighted 315 genes that were significantly enriched in 165 pathways(figure 1D). Approximately 20% of the enriched pathways were supported by data from both RNA and protein levels, and 24% of the pathways were identified by multiomics joint analysis. Notably, the type I interferon signaling and interferon regulatory pathways were frequently observed in these two categories. Additionally, our analysis highlighted the significant enrichment of the upstream pathway of interferon, the "RIG-I signaling pathway," in the non-pCR group(Qpathway=0.015).
Both RIG-I and MDA5 are members of retinoic acid-inducible gene I (RIG-I)-like receptors (RLRs). RLRs are RNA sensors localized in the cytosol and could induce type I interferon production in the detection of unusual nucleic acids.13 Our data support that the activation of both the RLR pathway and its downstream interferon signaling pathway contributes to neoadjuvant chemoradiotherapy treatment resistance in ESCC. It’s well established that the JAK-STAT pathway is the predominant downstream signal transduction pathway after interferon signal activation.14 In our differential expression gene analysis, STAT1 was significantly up-regulated in non-pCR group in both transcriptomics data(p-value=0.001,FC= 0.43) and proteomics data(p-value=0.03,FC= 0.48). Our findings are corroborated by previous research, which revealed that RIG-I-induced STAT1 activation leads to chemotherapy and radiotherapy resistance in breast cancer.15
Genetic characteristics of patients receiving NCRT with different response
A total of 18 biopsied samples were subjected to whole exome sequencing with the corresponding blood samples as controls, among which 9 samples(50%) were pCR after NCRT.
First, we evaluated the general genetic differences between the pCR and non-pCR groups (figure 2). Tumor mutation burden(TMB), described as mutations(SNVs or indels) per megabase, was not significantly different between groups(pCR median 10.32, SD3.63, IQR4.84; nonpCR median 9.36, SD4.75, IQR3.05; Mann-Whitney U test, p=0.840). The differences in other genetic variations(table S5), including total somatic CNA length(t-test p=0.608 for copy number gain, p=0.392 for copy number loss), somatic variation(SV) counts(t-test p=0.358), microsatellite instability (MSI)(t-test p=0.729), and tumor purity(t-test p=0.889), were not statistically significant between the two response groups. Gene mutation signatures were investigated via nucleotide mutation signature and the Catalog of Somatic Mutations in Cancer (COSMIC) database v3.16 the results revealed that the most common nucleotide mutations were C>T and T>C among all samples; there was no nucleotide mutational difference between the two groups(figure S2). The COSMIC signatures SBS25 and SBS9 were detected in all patients, and SBS9 was more frequently detected in the pCR group(figure S2). We further analyzed the differences in neoantigen load between the groups. Neither neoantigen load nor expressed neoantigen load was associated with the pathological response(figure S3).
Gene mutations
Next, we performed a significantly mutated gene (SMG) analysis of the two response groups. TP53 was the most mutated gene with a mutation frequency of 72% across all samples(figure 3A). FBXW7 was the only gene with significantly different mutation frequencies between the groups (Fisher exact test p-value=0.029), with five out of nine pCR patients carrying non-synonymous mutations, in contrast to no mutation in the non-pCR group. Mutation types included missense and frameshift mutations in the F-box, COG2319, and WD40 domains(figure 3B).
In a previous study, Gstalder et al. demonstrated that inactivation of FBXW7 in a melanoma model was associated with impaired tumor-intrinsic expression of the dsRNA sensors MDA5 and RIG-I, and thus decreased induction of type I interferon expression, which leads to anticancer therapy resistance.17 Given the observed differences in RLR and type I interferon signaling pathway expression levels between pCR and non-pCR patients, we sought to investigate if these expression variations were linked to FBXW7 mutation status. To do so, we conducted a differential gene expression and enrichment analysis, comparing samples with FBXW7 mutations to those with FBXW7 wild-type status using transcriptomic data. The results of GSEA revealed a notable reduction in the expression of the type I interferon signaling pathway in FBXW7 mutated samples (figure 3C and 3D). Meanwhile, GO enrichment analysis showed that the interferon upstream RLR pathways, including the MDA-5 signaling pathway(adjusted p-value =0.002) and the RIG-I signaling pathway(adjusted p-value=0.018), were both downregulated in FBXW7 mutated samples by GO enrichment analysis. RNA expression data revealed positive correlations between FBXW7 expression and key RLR pathway genes(Spearman correlation coefficient IFIH1 r=0.602, DDX58(RIG-I) r=0.400, TBK1 r=0.694), as well as type I interferon-induced JAK-STAT pathway key gene expression(Spearman correlation coefficient STAT1 r=0.449, JAK1 r=0.588, TYK2 r=0.306) (figure 3E) in 17 samples that underwent RNA-seq. Moreover, we analyzed the TCGA ESCC dataset (esophageal carcinoma, TCGA, Firehose Legacy) on cBioPortal,18 and verified the co-expression tendency of FBXW7 with key genes in the RLR and JAK-STAT pathways(measured by Spearman correlation coefficient, table S6).
In addition to its regulatory effect on RLR and interferon signaling, FBXW7 was also reported to contribute to DNA repair via nonhomologous end-joining and recovery of the cell cycle by promoting P53 degradation after DNA-damaging therapies. Inactivation of FBXW7 enhanced tumor radiosensitivity and chemosensitivity.19 The GO enrichment analysis in FBXW7 mutated samples in our cohort exhibited a significantly lowered expression in “DNA damage response” pathway(adjusted p-value=0.002).
Collectively, our data suggest that inactivation of FBXW7 in ESCC is sensitive to neoadjuvant chemoradiotherapy. It might exert an effect on the treatment response through the regulation of interferon-related pathways and DNA damage response pathways.
Copy number variations
Copy number variations(CNV) were investigated in all samples using GISTIC (figure S4). In the non-pCR group, the most significant copy number variations were 20q13.33 amplification(q-value=0.000) and 4q35.2 deletion(q-value=0.044). In the pCR group, 8q24.3 (q-value=0.000) was the top amplified cytoband, encompassing genes that encode a series of microRNAs and oncogenes such as PARP10. The most significant cytoband loss in the pCR group was 9p21.3 (q-value=0.001), which is a 917kb long chromosome band containing CDKN2A/B, MTAP, and a type I IFN gene cluster(figure 3F, 3G). Loss of chromosome 9p21.3 has been reported as the most common homozygous deletion in human cancers. While 9p21.3 deletion almost always involves CDKN2A, co-deletion of the IFN gene cluster only occurs in certain patients, which disrupts the tumor microenvironment and changes the tumor biological behaviors in these patients.20
In our pCR cohort, GISTIC results showed that five out of nine pCR patients carry 9p21.3 copy number variations, including CDKN2A deletion and co-deletion of IFN genes. We presume that the chromosome 9p21.3 IFN cluster deletion may partially contribute to the relatively lower expression of the type I IFN signaling pathway in the pCR group at both the transcriptomic and proteomic levels. In the ESCC TCGA cohort(esophageal carcinoma, TCGA, Firehose Legacy) from cBioPortal, 51 out of 96 patients carry 9p21.3 deletion, among which 13 patients had co-deletion of IFN genes. Initially, we found no significant difference in the expression level of the interferon signaling pathway between patients with and without non-IFN cluster 9p21.3 copy number deletion. However, when comparing patients with and without 9p21.3 deletions that included the IFN cluster, we observed a significant increase in the expression of RIG-I, IFN signaling, and JAK-STAT pathways in the copy number deletion group.(figure S4). 36 patients in the TCGA ESCA cohort had recorded response evaluation of radiotherapy, among which three patients were found to carry 9p21.3 IFN cluster deletion, and all of their radiotherapy responses were complete remission (CR).
The findings in the TCGA cohort reaffirmed our assumption that chromosome 9p21.3 copy number deletion with IFN-cluster is a sensitive genetic variance that enhances ESCC neoadjuvant chemoradiotherapy treatment response, and it might exert its pro-therapeutic effect by downregulating interferon signaling in the tumor cells.
Assessment of the IFN effect on ESCC cell line radiosensitivity
Previous studies have reported that constitutive elevation of type I IFN was seen in some cancers.21 And the chronic exposure to low dose type I IFN may promote cell resistance to DNA damage through elevated expression of a subset of IFN induced genes including STAT1, STAT2 and IRF9.22 To investigate whether chronic exposure to type I IFN promotes radiotherapy resistance of esophageal squamous cell cancer, we treated the ESCC cell line TE1 with low-dose IFNβ(5 IU/ml) continuously for 16 days. The treated cells and untreated control received irradiation of 0, 2, 4, 8Gy, and apoptotic percentages were measured 48h after radiation.
In comparison to the control group, TE1 cells treated with 5 IU/ml IFNβ for 16 days exhibited a significantly reduced number of apoptotic cells following exposure to 2 Gy(t-test p=0.031), 4 Gy(t-test p=0.018), and 8 Gy(t-test p=0.012) radiation(figure 4A-B), indicating increased radioresistance after low-dose IFNβ treatment. The results supported the hypothesis that the primary radioresistance observed in NCRT treatment may derive from the constitutive activation of type I IFN signaling in ESCC.
Knockdown of FBXW7 alters intracellular IFN levels
To evaluate the effect of FBXW7 mutation on the ESCC type I IFN expression levels, we used siRNA to knockdown FBXW7 in KYSE510 cell line. After transfection for 48 hours, the expression of FBXW7 was downregulated to 40%. The intracellular levels of both IFNα and IFNβ were significantly decreased in the FBXW7-siRNA group compared with the control group(IFNα p=0.000, IFNβ p=0.010)(figure 4C). This result supports our hypothesis that the loss-of-function mutation of FBXW7 could lower the constitutive expression level of type I IFN in ESCC, and may further alter the radiosensitivity through this effect.
Interferon stimulated genes(ISGs) as resistance biomarkers of ESCC patients treated with neoadjuvant chemoradiotherapy
After neoadjuvant chemoradiotherapy, patients with pathological complete remission have distinct gene expression profiles at the protein level compared with non-pCR patients. We aimed to define a protein biomarker panel that can identify NCRT-resistant tumors. Protein-protein interaction(PPI) analysis was performed using 385 differentially expressed proteins in the two response groups. Next, the genes in the PPI network were ranked based on their topological score calculated using the cytoHubba MCC method.23 45 hub genes were identified at a cutoff score of 50000. The hub genes were further selected using transcriptomic data, and the final protein biomarkers were defined as the intersection of PPI hub genes and transcriptomic differential expression genes(DESeq2, p-value<0.05)(table S7). The biomarker panel comprised 12 proteins: STAT1, EIF2AK2, MX1, BST2, TRIM21, SAMHD1, IFI44L, GBP1, PARP14, ISG15, HLA-B, and IFIT3. All proteins had relatively elevated expression in the non-pCR group; thus, the panel could identify ESCC patients with potential NCRT resistance. The corresponding genes of these proteins were all members of the interferon-stimulated gene family, and active interactions were observed within their protein-protein interaction (PPI) network(figure 5A).24
The 12-protein biomarker panel was tested in the 36-patient subset with a recorded radiotherapy response from TCGA ESCA cohort(table S8). GSEA demonstrated significant enrichment of the biomarkers in non-complete response(radiographic progressive disease, stable disease) patients(NES -1.852, adjusted p-value 0.002)(figure 5B), as expected. Stratified by high and low levels of the biomarker panel expression, Kaplan-Meier plotter (kmplot.com/analysis/) was used to assess the survival outcomes of ESCC patients from databases including TCGA,EGA, and GEO etc..12 Patients with higher biomarker expression level had significantly impaired overall survival(median OS, 18.9 months) compared with the low-expression subset(median OS, 42.1 months)(HR:2.29, 95% CI 1.00-5.23; log-rank p = 0.044)(figure 5C).
To further test the discriminative potential of the 12 NCRT-resistant biomarkers, we performed immunofluorescence staining of the biomarkers using FFPE slides from an independent cohort of NCRT-treated ESCC patients. Each biomarker was stained on pre-treatment(NCRT) biopsied primary tumor samples from 5 pCR patients and 5 non-pCR patients. Biomarker positive score(BPS) was defined as the percentage of tumor cells with positive biomarker staining at any intensity(the number of positive tumor cells divided by the total number of tumor cells). The BPS values of the 12 biomarkers are shown in figure 5D(p-values for t-test are listed in table S9). Ten out of the 12 biomarkers had significantly elevated expression in the tumor cells of the non-pCR patients using two-sample t-test, including IFIT3(p=0.000), BST2(p=0.000), SAMHD1(p=0.000), HLA-B(p=0.000), EIF2AK2(p=0.001), GBP1(p=0.002), IFI44L(p=0.003), MX1(p=0.006), STAT1(p=0.026), ISG15(p=0.038) (figure 5E and figure S5).