Prognostic value of genomic mutation signature associated with immune microenvironment in southern Chinese patients with esophageal squamous cell carcinoma

We retrospectively identified patients who underwent surgical resection of the esophagus and were histologically confirmed ESCC at Fudan University Shanghai Cancer Center (FUSCC) and Sun Yet-Sen University Cancer Center (SYSUCC) from 2007 to 2008. After obtaining informed consent, cases were interviewed to collect information on demographics and clinical information. Exclusion criteria included missing baseline clinicopathological features, mixed histology and incomplete follow-up data. This study was conducted in accordance with the provisions of the Declaration of Helsinki and was approved by Institutional Review Boards of FUSCC and SYSUCC.

Data on somatic mutation, and/or the mRNA expression data, and the corresponding clinicopathological characteristics of ESCC patients from the ICGC database (https://dcc.icgc.org/↗) and The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/↗) were acquired to make comparative analysis. Clinical characteristics included age, TNM stage, survival status, and survival time.

Tumor tissues obtained during surgery were formalin-fixed, paraffin-embedded (FFPE), and stored at -130℃ until used. The specimens were chosen for this study based on two criteria: (a) histological diagnosis of ESCC independently confirmed by two experienced pathologists at FUSCC and SYSUCC; (b) availability of purity tumor tissue (at least > 75%). To address the FFPE-DNA damage, mitigation strategies were used in the process of sample quality control, DNA repair treatments, and analytical sample preparation, as suggested by previous studies [16, 17].

Genomic DNA was prepared for paired-end library construction using a QIAamp DNA FFPE Tissue Kit and QIAamp DNA Blood Midi Kit (Qiagen, Hilden, Germany). FFPE-induced sequencing artifacts were controlled mainly through the following aspects. Firstly, for DNA extraction, a FFPE sample underwent a pathologist review to ensure each sample at least had the area of one square centimeter, nucleated cellularity of 80% and tumor content of 20%. And specialized FFPE sample reagents were utilized for DNA extraction. Secondly, data quality was inspected and controlled by examining sequencing coverage and uniformity. The variation of CT/GA was eliminated if the ratio of transition and transversion was too high.

The genomic alterations were examined using the YuanSuTM panel (OrigiMed, Shanghai). This panel covers 637 genes important in the oncogenesis of solid tumors. The gene list is available upon request. The targeted regions of YuanSuTM panel were captured and sequenced at a mean depth of 1000× with an Illumina NovaSeq 6000. The details of genomic alteration detection including copy number variations (CNVs) were listed in the supplemental materials.

TMB was estimated by counting the coding somatic mutations, including single-nucleotide variants (SNVs) and insertions or deletions (INDELs), per megabase of the sequence examined in each patient. After calculating the TMB of each patient in this method, we divided the patients into high- and low-TMB groups according to the median or the quantile value for subsequent analyses.

The top 50 most frequently mutated genes were selected and inputted into the GeneMANIA website (http://genemania.org↗) to generate the PPI network. The “org. Hs.eg.db,” “ggplot2,” “clusterProfiler,” and “enrichplot” packages were utilized for analysis of the Kyoto Gene and Genome Encyclopedia (KEGG) and gene ontology (GO) pathways; FDR < 0.05 was considered to be statistically significant. Immune-related gene set enrichment analysis (GSEA) of top 200 mutated genes was used to visualize the different biological processes in different groups. Gene sets with |NES|> 1, NOM p < 0.05, and FDR q < 0.25 were considered to be statistically significant.

We randomized our ESCC cohort in a 7:3 ratio, with 70% of the dataset as the training set and 30% as the validation set. In addition, ICGC ESCC dataset was used as our independent validation sets. Clinical characteristics and gene mutation status were used as candidate features. The important features were selected by minimal depth method and predicted in the validation group. To evaluate the model, receiver operating characteristic curves were drawn and the areas under the curves (AUCs) were calculated. To evaluate the clinical significance of the selected features, we further conducted Cox regression analysis. A nomogram was constructed to show the results of Cox analysis. The R packages used in this procedure included randomForestSRC, glmnet, dplyr, survival, surviminer, ggplot2, forestplot, survcopm, and prodlim.

Immune cell abundance (immune score), stromal cell infiltrating level (stromal score), and tumor purity (ESTIMATE score) were estimated via the Estimation of Stromal and Immune cells in Malignant Tumor tissues using Expression (ESTIMATE) algorithm [18]. Potential response of immune checkpoint blockade therapy was estimated with Tumor Immune Dysfunction and Exclusion (TIDE) algorithm (http://tide.dfci.harvard.edu/)↗ [19]. Using Tumor Immune Estimation Resource 2.0 (TIMER 2.0, http://timer.cistrome.org/)↗ [20], a comprehensive analysis of immune infiltration in the tumor microenvironment of ESCC was carried out. MCPcounter algorithm was used to estimate the relative proportions of 22 immune cells in ESCC. The infiltration of 22 immune cells was indicated by enrichment scores, which were calculated by single sample GSEA using Gene Set Variation Analysis R package.

Two different FFPE blocks of tumor tissues were obtained from each patient and spotted in adjacent 2.5-mm cores as per the tissue microarray map to capture the tumor heterogeneity. Briefly, the regions of interest were marked by experienced pathologists, and then a recipient paraffin block was made, incubated, cooled, fixed, and finally stored with paraffin. The details of the construction process were listed in the supplemental materials. A panel of two markers was used to explore the macrophage infiltration in ESCC, which included (a) marker of epithelial cells and tumor cells: cytokeration; and (b) marker of macrophages: CD68. A panel of two markers was further used to explore the infiltration of different types of macrophages, which included (a) marker of pro-inflammatory M1 macrophages, CD86; and (b) marker of anti-inflammatory M2 macrophages, CD163. The details of the staining process were listed in the supplemental materials.

All statistical analyses were conducted using SPSS 24.0, R software (version 4.0.2), and GraphPad Prism 8.0. Student’s t test was used to compare continuous variables, and the Chi-square test or Fisher’s exact test was used to compare categorical data. Fisher’s exact test was used to analyze the correlations between TMB and gene mutations that co-occurred or were exclusive of each other. Cox regression analysis was used to determine risk mutated genes. The clinical effects of the prognostic model were assessed using the log-rank test and Kaplan–Meier method. A nomogram model was then constructed, and predictive performance of the nomogram was estimated using C-index and calibration plot.

Primary tumor tissue samples were obtained from 102 patients with histologically confirmed ESCC. After discarding samples that failed to pass the quality check, a final cohort of 92 patients was available for analysis, with 32 from FUSCC and 60 from SYSUCC. All the patients had TGS data and CNV data. Figure 1 demonstrates the study design. The detailed clinicopathological information was shown in Supplemental Table 1.

Among the 92 ESCC samples with qualified TGS data, 1,979 somatic mutations were identified, comprising 1,799 SNVs and 168 INDELs. These tumors harbored a median of 22 nonsynonymous SNVs and 2 INDELs. As demonstrated in Fig. 2A, the most prominent mutated genes were TP53 (found in 94.6% of tumors), followed by KMT2C (27.2%), KMT2D (27.2%), LRP1B (27.2%) and FAT1 (25%). As for the mutation sites, ALOX12B (c.1565C > T), SLX4 (c.2786C > T), LRIG1 (c.746A > G), and SPEN (c.6915_6917del) were the most frequent (6.5%), followed by RELA (c.1033 + 6del) at 5.4%. Some of these mutation sites have never been reported before (Fig. 2B). To assess somatic CNVs, reported oncogenes and tumor suppressor genes observed in the top rank (residual q<1×10-4) GISTIC peaks were examined (Fig. 2C). FGF19 was the gene most frequently affected by somatic CNVs, with alterations in 94.6% of patient samples. Frequent amplifications were observed in FGF3 (43%), FGF4 (42%), FGF19 (42%), TP53 (39%), and CCND1 (35%). Frequent deletions were observed in CDKN2A (42%) and CDKN2B (36%). Frequent losses were observed in BIRC3 (45%) and YAP1 (38%), and frequent gains in EGFR (45%). The mean TMB was 8.2 mutations per megabase (mut/Mb), and median TMB was 7.5 mut/Mb. After the 92 samples were divided into four groups by TMB quantiles, gender was the only one among the clinicopathological features to be significantly related to TMB (p < 0.001, Fig. 2A). A total of 12 mutated genes were significantly associated with TMB (p < 0.05), including ZNF750, TP53, TET1, ROS1, LRP2, KMT2D, GLI2, AR, FAT1, SPEN, KMT2C, and DPYD (Supplemental Fig. 1).

We integrated the top 200 most frequent mutated genes into pathway analysis to identify the functionally aberrant pathways in ESCC. As a result, we identified a highly interconnected network of aberrations targeting the cell cycle regulators, epigenetics regulators, PI3K/AKT signaling pathway, and NOTCH signaling pathway (Fig. 3). Several tyrosine kinases (that is, FGF3, FGF4, FGF12, and FGF19) and receptors (that is, EGFR and ERBB2) exhibited CN gains, as well as their downstream signal transducer PIK3CA and AKT3, together contributed to the proliferation of ESCC. The loss of CDKN2A and CDKN2B, together with the gain of CCND1, CCND3, CDK4, and CDK6, resulted in the disruption of cell cycle restriction.

To determine whether the mutated genes functionally interacted with each other and were involved in tumorigenesis, the top 50 frequently mutated genes were selected to perform a PPI network-based analysis and functional pathway analysis. The PPI networks are visualized in Fig. 4A. These genes and their encoded proteins were involved in positive regulation of cell cycle process, regulation of protein kinase B signaling, intrinsic apoptotic signaling pathway, and regulation of NOTCH signaling pathway, which was similar to the results in the above aberrant pathway analysis. The co-occurrence and exclusion relationship of the mutated genes in ESCC is demonstrated in Fig. 4B. Among the top 50 frequently mutated genes, KMT2C and TET1 had the most frequent co-occurrence relationship with other 13 genes in mutation, hinting the important role of epigenetic regulators in ESCC pathogenesis and development. KEGG and GO pathway analysis demonstrated enrichment of mutated genes in several signaling pathways involved in malignancy, such as the PI3K/AKT pathway, cell cycle, and NOTCH signaling pathway, all of which were also enriched in TCGA-ESCC cohort (Fig. 4C and D)

Random survival forest model was applied to determine the mutated genes of most significance to the survival of ESCC patients. According to the variable importance and minimal depth[21], 19 variables were finally selected, including 17 mutated genes (BRCA2, NOTCH2, ARHGEF17, NCOA2, PIK3CG, ZNF703, NF1, MECOM, ATM, FAM135B, SERPINB4, ZNF750, PARP2, NOTCH1, RECQL4, EPHB4, and TET1), age and stage (Fig. 5A). Among the 17 genes, NOTCH2, ARHGEF17, NCOA2, and ZNF750 were significantly associated with the OS of patients in our cohort (Supplemental Fig. 2). Based on the model, each patient was scored. The median risk score value was used to divide the samples into high-risk and low-risk groups both in our cohort and in the ICGC cohort. The survival curves can be significantly separated between the two groups in both cohorts (p < 0.001, Fig. 5B and C). To validate this model, receiver operating characteristic curves were drawn, and the AUCs for 36 months, 24 months, and 12 months were 0.890, 0.861, and 0.780, respectively in our validation set (Fig. 5D) and 0.695, 0.627, and 0.628, respectively in the ICGC validation set (Fig. 5E). Then, a predictive nomogram based on these risk factors was constructed to clearly show the survival rate (Fig. 6A). Calibration curve indicated that the observed and predicted values were consistent in predicting OS in our cohort and the ICGC cohort (Fig. 6B and C).

Of note, TMB was not selected under this screening method. We further examined the association between TMB and the OS of patients in our cohort. The survival curves twined together in the patients grouped by either the median value or the quantile values of TMB (Supplemental Fig.). 3

To further reveal the potential relationship between the risk score and immune microenvironment in ESCC, the TCGA-ESCC cohort with both mutation data and mRNA expression data was analyzed. We used GSEA to visualize the enriched immune-related biological processes in different risk groups classified by the mutated genes-related risk model in the TCGA-ESCC cohort. Our results indicate that patients in the low-risk group were prone to show associations with Th17 cell differentiation, regulation of Th17 cell differentiation, IL-1 receptor binding, cytokine receptor interaction, natural killer cell mediated cytotoxicity, JAK-STAT signaling pathway, B cell receptor signaling pathway, and T cell receptor signaling pathway (Fig. 7A). Compared to the low-risk group, the high-risk group was enriched in fewer immune-related pathways. Using the ESTIMATE algorithm, we found that the EstimateScore was significantly lower in the low-risk score group (Fig. 7B). The relationship between the risk score and the immune checkpoint expression level was then explored. The risk score was in a significant correlation with several stimulatory immune checkpoints (p < 0.05), including TNFSF4, ITGB2, CXCL10, CXCL9, and BTN3A1 (Fig. 7C). And the risk score tended to have a correlation with many stimulatory or inhibitory immune checkpoints with p < 0.1 (Supplemental Fig. 4). To predict the response to immunotherapy, the TIDE algorithm was used to divide the patients into responders and non-responders (Fig. 7D). The key markers in the TIDE algorithm (CD274, IFNG, and TAMM2) exhibited significantly different levels in the patients with different risk scores (Fig. 7E). Taken together, these results suggest that the 17-mutated gene-related risk model might play an additional role in hinting the immune status in ESCC.

The potential relationship between the risk score and infiltrating immune cells was further explored in the TCGA-ESCC patients grouped by the median of the risk score based on the 17-mutated gene-related model. Figure 8A demonstrates the proportion of different infiltrating immune cells in each sample. Among the 22 types of immune cells, the proportion of macrophage ranked first, followed by T cells CD4 + memory resting, M2 macrophage and M1 macrophage (Supplemental Fig. 5). Compared to patients in low-risk group, patients in high-risk group were estimated to harbor a significantly higher proportion of M2 macrophage and a lower proportion of activated dendritic cells (Fig. 8B). As the proportion of macrophage ranked first, we then unveil the profile of macrophage infiltration in ESCC in our cohort. Macrophage was marked and classified as pro-inflammatory M1 macrophage and anti-inflammatory M2 macrophage by immunofluorescence double staining. Macrophage was detected in ESCC tissue from patients in both low-risk group (n = 14) and high-risk group (n = 12), with no significant difference in macrophage level (Fig. 8C and E). M1 macrophage was observed in a similar level in the two groups; however, patients in the high-risk group had a significantly higher level of M2 macrophage than patients in the low-risk group, hinting the potential role of M2 macrophage instead of M1 macrophage in the regulation of ESCC immune microenvironment (Fig. 8D and E).

The authors declare that they have no competing interests.

This study was conducted in accordance with the provisions of the Declaration of Helsinki and was approved by Institutional Review Boards of Fudan University Shanghai Cancer Center and Sun Yet-Sen University Cancer Center.

Not applicable.

Datasets could be downloaded directly from the indicated websites or are available from the corresponding author on reasonable request.