A fibroblast-associated signature predicts prognosis and immunotherapy in esophageal squamous cell cancer

ESCC scRNA-seq data was acquired in Gene Expression Omnibus (GEO) database (accession number GSE191756). We screened out single cells with any gene expressed in fewer than three cells or those expressing fewer than 250 genes. The percentage of rRNA and mitochondria was then calculated with the PercentageFeatureSet function in the Seurat R package (22). Consequently, 12118 cells were totally obtained for subsequent analysis.

We further collected transcript data, single-nucleotide variants (SNV), copy number variants (CNV), and corresponding clinical data of ESCC from The Cancer Genome Atlas (TCGA) database. Samples lacking outcome status or survival data were excluded and 94 ESCC samples were obtained, which were utilized for external validation. GSE53624↗ data with 119 tumor samples and 119 normal ones was used as the training cohort after abandoning samples without follow-up acquired from Gene Expression Omnibus (GEO) database (The clinical characteristics of both the training cohort and the test cohort were exhibited in Supplementary Tables). Based on the literature, ten cancer-associated pathways (Cell Cycle, NRF1, MYC, NOTCH, HIPPO, PI3K, TP53, PI3K, WNT, and TGF-Beta) were identified and analyzed about their gene expression profiles in our dataset.

The scRNA-seq data of ESCC was re-analyzed using the Seurat package, with the aim of providing a systematic characterization of the CAF signature. To start the data preprocessing, cells with less than 250 or more than 6000 expressed genes were removed, and the remaining expressed genes were log-normalized. Next, the FindIntegrationAnchors function was utilized. To reduce the dimensionality of the data, the t-distributed Stochastic Neighbor Embedding (tSNE) method was applied, utilizing a resolution of 0.1 and selecting 30 principal components. The tSNE method employed was non-linear in nature. To classify the single cells into various subgroups, we utilized the FindNeighbors and FindClusters functions (dim = 30 and resolution = 0.1). Additionally, we performed tSNE dimensional reduction using the RuntSNE function. Fibroblasts were annotated according to four marker genes, including PDGFRB, ACTA2, FAP, and NOTCH3. Subsequently, the fibroblasts were re-clustered using the FindClusters and FindNeighbors functions. To define the marker genes for each CAF cluster, we used the FindAllMarkers function with a comparison between different clusters (minpct = 0.35, logFC = 0.5, and adjust p-value < 0.05). We utilized the CopyKAT R package to analyze the CNV characteristics of the CAF clusters and distinguish them from tumor cells and normal ones. Finally, Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis was conducted on the marker genes using the clusterProfiler package (23).

Differentially expressed genes (DEGs) between normal and tumor tissue were identified using the limma package (24, 25), based on criteria of |log2(FoldChange)|>1 and a false discovery rate (FDR)<0.05. Next, the correlations between CAF clusters and DEGs were evaluated, and key CAF-related genes with p<0.01 and cor>0.4 were identified. To identify prognosis-related genes, the survival package was utilized to conduct univariate Cox regression analysis (26). The least absolute shrinkage and selection operator (lasso) was used to reduce the number of genes (27). Multivariate Cox regression analysis was conducted using the stepwise regression method to establish a CAF-based risk signature, which was calculated using the formula: 0.093*ANGPTL7 + 0.15*C6 + 0.121*CSRP1+-0.08*EXPH5 + 0.12*F2RL2 + 0.014*KCNMA1+-0.373*MAGEC3 + 0.143*MAMDC2+-0.188*SLC4A9. The patients were classified into low- and high-risk groups using zero-mean normalization. The predictive value of the risk signature was evaluated using the timeROC package to perform receiver operating characteristic curve (ROC) analysis. The results demonstrated that the risk signature had significant predictive value for patient prognosis. In summary, our analysis provides important insights into the molecular mechanisms underlying tumor development and highlights the potential of CAF-related genes as prognostic biomarkers for cancer patients.

After conducting univariate and multivariate Cox regression analyses based on the risk signature and clinicopathological features (24, 28), a novel nomogram was constructed to predict the prognosis of ESCC using variables with p<0.05 in the multivariate Cox model. The predictive accuracy of the model was evaluated by generating a calibration curve.

The correlation between the risk signature and the tumor immune microenvironment (TIME) was comprehensively assessed using several algorithms, including CIBERSORT, EPIC, MCPCOUNTER, and TIMER (29). Stromal scores, immune scores, and estimate scores (stromal scores + immune scores) were calculated using the “estimate” R package to evaluate differences in the tumor microenvironment of patients (30, 31). Besides, the proportions of 22 immune cell subtypes were estimated using the CIBERSORT algorithm based on the GSE53624↗ cohort. The correlation between genes comprising the signature and immune score were further explored to illuminate the great impact those genes exert on immune-related functions.

Anti-PD-1 or anti-PD-L1 checkpoint inhibition therapy has gained increasing attention as a crucial component of immunotherapy. Transcriptomic data as well as corresponding clinical data from patients who received anti-PD-L1 therapy from the IMvigor210 cohort were collected to evaluate the performance of the risk signature in predicting responsiveness to immunotherapy (immune checkpoint blocks). Additionally, transcriptomic data from the GSE78220↗ cohort, which included melanoma patients who received anti-PD-1 checkpoint inhibition therapy before treatment, were downloaded.

To further probe into the heterogeneity of ESCC, all ESCC patients were separated into different clusters according to the expression of CAF-related genes with the R package ‘ConsensusClusterPlus’ (22). Differences in survival, TIME, and immune checkpoints were evaluated among subgroups using the same methodology as previously employed. The immune landscape of ESCC patients based on different clusters was demonstrated in the form of heatmap.

Total RNA was isolated from ESCC cells or tissues using TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, USA). The complementary DNA (cDNA) was synthesized as per the manufacturer’s instructions, utilizing the RevertAid™ First Strand cDNA Synthesis Kit (Thermo Fisher Scientific). qRT-PCR was performed with SYBR Green PCR kit (Takara Bio, Otsu, Japan) on a StepOne Real-Time PCR system (Thermo Fisher Scientific). The relative gene expression levels were quantified by employing the 2-△△CT method.

All statistical analyses were performed using R software (version 4.1.0). The Wilcoxon test was used for comparing two groups, while Spearman or Pearson correlation was used for correlation matrices. Survival differences through K-M curves were assessed using the Log-rank test, where statistical significance was defined as p-value < 0.05.

The flow-process diagram of our study was depicted in Figure 1. After initial screening, 18024 cells were totally obtained based on the scRNA-seq data. As was shown in Figure S1, the details of data preprocessing were demonstrated. After conducting log-normalization and dimensionality reduction, 32 subpopulations were obtained (Figure 2A). As presented in Figure 2B, six CAF populations were further identified with four marker genes (FDGFRB, FAP, ACTA2, and NOTCH3). The proportions of the six clusters in each cohort were then calculated and the results were illustrated in histograms (Figure 2C). Moreover, KEGG analysis illuminated that the DEGs (which were obtained using R package ‘FindVariableFeatures’) were significantly enriched in various pathways, including tight junction, complement and coagulation, focal adhesion and so on (Figure 2D). Additionally, distributional differences between tumor and normal cells in the six CAF clusters were presented in Figure 2E. In addition, the expression of TOP 5 DEGs were respectively exhibited in heatmap (Figure 2F), bubble diagram (Figure 2G), and volcano plot (Figure 2H).

To probe into the correlations between tumor progression and the CAF clusters, we explored the features of ten tumor-related pathways based on the six CAF clusters. GSVA scores of these divergent pathways were estimated based on various CAF clusters, and the results were depicted in Figure 3A. Significant differences about the ratio of malignant cells were obtained in CAF_0 and CAF_4, where the malignant cells accounted for a few proportions. By contrast, the ratio of malignant cells was remarkably higher among CAF_1, CAF_2, CAF_3, and CAF_5 (Figure 3B). Besides, slight differences were identified after performing the GSVA analysis of these tumor-related pathways between non-malignant and malignant cells in each CAF cluster (Figures 3C-F). (GSVA scores analysis based on CAF_0 CAF_1 was shown in Figures S2A, B)

Furthermore, the ssGSEA score of the marker genes (the TOP 5 DEGs obtained in Figures 2F-H) were analyzed in each CAF cluster based on the GSE53624↗ cohort to illustrate the relationships between the CAF clusters and crucial clinicopathological characteristics. Interestingly, tumor samples were found having higher scores compared with normal ones only in CAF_0 and CAF_3, while among the other clusters, normal samples gained significantly higher scores (Figure S2C). In addition, ESCC samples of GSE53624↗ cohort were divided into high-and-low score groups according to the optimal cut-off value with survminer R package. In the CAF_2, CAF_4, and CAF_5 clusters, samples in the low-CAF score subgroup shared a more favorable prognosis compared with those in high-CAF subgroup. CAF_0, CAF_1, and CAF_3, however, were identified not associated with the prognosis of ESCC (Figure S2D).

Firstly, DEGs were screened out between normal and tumor samples to establish a risk signature. As depicted in Figure 4A, 17080 DEGs were totally obtained, with 7556 down-regulated and 9524 up-regulated DEGs. Among them, a total of 642 genes were identified significantly related with those prognosis-related CAF clusters (Figure 4B). After univariate Cox regression analysis, the prognosis of each gene was evaluated, with 8 genes being identified related to protective factors and 18 genes correlated with risk values. Lasso Cox regression analysis was then performed to reduce the number of genes (Figure 4C). Furthermore, the stepwise regression method was utilized to develop the risk signature after performing multivariate Cox regression analysis. The signature was composed with nine genes (Figures 4F, G), namely complement C6, MAM domain containing 2 (MAMDC2), cysteine- and glycine-rich protein 1 (CSRP1), coagulation factor II thrombin receptor like 2 (F2RL2), angiopoietin like 7 (ANGPTL7), potassium calcium-activated channel subfamily M alpha 1 (KCNMA1), exophilin 5 (EXPH5), solute carrier family 4, sodium bicarbonate cotransporter, member 9 (SLC4A9), and MAGE family member C3 (MAGEC3). And the risk model formula is as follows: RiskScore=“0.093*ANGPTL7 + 0.15*C6 + 0.121*CSRP1+-0.08*EXPH5 + 0.12*F2RL2 + 0.014*KCNMA1+-0.373*MAGEC3 + 0.143*MAMDC2+-0.188*SLC4A9”. The risk score of each sample was calculated using z-mean normalization, and the patients were then separated into high-and-low-risk groups. The Kaplan-Meier survival analysis exhibited that patients in high-risk groups encountered with worse prognosis compared with those in low-risk groups both in GSE53624↗ (Figure 4D) and TCGA cohorts (Figure 4E). Additionally, both the GSE53624↗ and TCGA cohorts exhibited satisfying AUC values of the model, revealing that the predictive power of the signature was excellent. The last, the distribution of patient survival status, risk score, and expression of hub genes in GEO and TCGA cohorts were depicted in Figure S3.

To improve the accuracy of our predictive model, we incorporated clinicopathological characteristics and the risk score through univariate and multivariate Cox regression analyses. Our multivariate analysis revealed that the risk signature was the most significant independent prognostic factor for ESCC, with a p-value of less than 0.001 (Figures 5A, B). We have developed a new nomogram that incorporates T-stage, N-stage, and the risk score (shown in Figure 5C). Through calibration plot analysis, this nomogram was found to have strong predictive power for actual survival outcomes (Figure 5D). The TimeROC analysis in the TCGA cohort has confirmed that the area under the curve (AUC) of both the nomogram and risk score outperformed other indicators (Figure 5E).

To explore the fundamental functions those DEGs play in initiation and progression of ESCC, Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) analysis were conducted. As shown in Figures 6A–E, the up-regulated genes were most enriched in base excision repair, cell cycle, and DNA replication, while the down-regulated genes were significantly correlated with arachidonic acid metabolism, calcium signaling pathway, and histidine metabolism. Likewise, the results of GO analysis were presented in Figures 6B-D. Furthermore, Gene Set Enrichment Analysis was conducted based on the nine genes involved in the risk signature. The results illustrated that 7 pathways were remarkably associated with these nine genes (Figures 6F, G). Interestingly, olfactory transduction was positively correlated with the genes except for EXPH5, MAGEC3, and SLC4A9, which were identified to have protective values in ESCC, indicating that olfactory transduction might suppress the immigration and progression of ESCC.

After conducting an investigation into the landscape of immune and stromal cell infiltrations in both low- and high-risk groups, Figure 7A demonstrated that patients in the high-risk group exhibit higher proportions of immune and stromal cell infiltrations when compared to those in the low-risk group. Besides, the immune cells proportions between the low-and-high-risk groups were estimated using the CIBERSORT algorithm. It was found that the high-risk group had higher proportions of resting memory CD4 T cells, Macrophages (M2), and resting mast cells, while naive B cells were more enriched in low-risk group (Figure 7B). Figure 7C exhibited the results of immune-related functions differences between high-and-low-risk groups.

Additionally, the relationship between risk genes and immunity was probed into. On the one hand, Figures 7D–H demonstrated that the protective genes (including EXPH5, MAGEC3, and SLC4A9) were found negatively correlated with stromal score, immune score, and estimate score. On the other hand, the risk genes (including PR1, F2RL2, KCNMA1, and MAMDC2) were identified positively correlated with divergent immune cells (Figures 7F, I). Finally, the reciprocal communication between the 75 immune-related genes and the nine model genes were displayed in Figure 7G.

Under the circumstances that T-cell immunotherapy has gained great achievements in recent years, we performed the assessment of prognostic value of our signature for immune-checkpoint therapy in GSE78220↗ and IMvigor210 cohorts. Divergent degrees of responsiveness of anti-PD-L1 receptor blockers were identified in the 348 patients from the IMvigor210 cohort, including partial response (PR), complete response (CR), progressive disease (PD), and stable disease (SD). As depicted in Figures 8A-C, patients in the high-risk group accounted for more proportions in PD/SD, and had worse prognosis than those in the low-risk group. Besides, SD/PD patients tended to gain higher risk scores than CR/PR patients. However, significant survival differences were identified neither in Stage I+II nor in Stage III+IV patients between the different risk subgroups (Figures 8D, E). To validate our findings, the GSE78220↗ cohort was enrolled for further analysis. Corresponding with the results from IMvigor210, patients who achieved partial or complete response had lower risk scores and were less likely to be in the high-risk group (Figures 8F-H).

Moreover, unsupervised consensus clustering was conducted to explore molecular subtypes based on the expression of CAF-related genes comprising the risk signature. With k = 3 deemed as the optimal clustering stability, patients in GSE53624↗ cohort were grouped into three clusters (Figure 9A). The ridge plot exhibited the distribution of various clusters (Figure 9B). Besides, as presented in Sankey diagram (Figure 9D), Cluster 1(C1) and Cluster 3(C3) made up the low-risk group while the high-risk group was comprised of Cluster 2 (C2) and Cluster 3 (C3). Subsequent survival analysis illustrated those patients in the C1 group had the most favorable prognosis, while patients in the C3 group had the worst clinical outcomes (Figure 9C). The immune landscape based on different clusters were shown in heatmap (Figure 9E), which indicated that C2 cluster beard the highest immune cell infiltrations. The TME scores of varying clusters were then calculated, revealing that C2 cluster had the highest immune, stromal, and estimate score as well as lower tumor purity than the other clusters. (Figures 9F-I). After applying the immune checkpoints inhibitors analysis, it was identified that C2 cluster was significantly correlated with BTLA, CTLA4, CD48, and so on, suggesting that patients involved in the C2 group might benefit from immunotherapy, especially anti-PD-1 receptor blockers.

After comparing the efficacy of various chemotherapeutic agents across different clusters, we discovered that patients belonging to cluster 2 (C2) exhibited elevated IC50 values for chemotherapeutic medications such as Bosutinib, Gefitinib, and AICAR. Additionally, patients in cluster 1 (C1) were observed to be more receptive to AMG.706, IPA.3, and the like, while those in cluster 3 (C3) demonstrated poorer response rates to the majority of chemotherapeutic treatments (Figure 10).

To explore potential ESCC cancer risk-related genes, four genes involved in the risk signature were selected for further validation in ESCC patients. As demonstrated in Figure 11, F2RL2 exhibited elevated expression levels in tumors, whereas SLC4A9, EXPH5, and MAGEC3 exhibited significantly reduced expression levels in tumors. These distinctions align with our bioinformatic findings, suggesting that these genes may serve as innovative biomarkers for early ESCC diagnosis.