Paired whole exome and transcriptome analyses for the Immunogenomic changes during concurrent chemoradiotherapy in esophageal squamous cell carcinoma

Background The immunogenomic changes triggered by concurrent chemoradiation therapy (CCRT), a standard neoadjuvant treatment for locally advanced esophageal squamous cell carcinoma (ESCC), are unknown. We aimed to analyze the early immunogenomic changes in ESCC induced by CCRT and to correlate them with clinical outcomes. Methods We collected biopsy samples from 40 patients with ESCC and the surgical candidates were treated with 5-fluorouracil (5-FU)/Cisplatin and concurrent radiation therapy. Endoscopic biopsy was performed before and after one treatment cycle of 5-FU/Cisplatin and 5 to 18 fractions of radiation. We analyzed immunogenomic changes using paired whole-exome sequencing (n = 29) and paired whole-transcriptome sequencing (WTS, n = 23). Multiplex immunohistochemistry (IHC) was conducted in four representative pair samples. Results Fourteen out of 23 WTS samples (60.8%) showed increased immune scores after CCRT, as calculated by ESTIMATE. The rate of progression-free survival was higher in patients with increased immune scores compared with the remaining patients (83.1% vs. 57.1%, p = 0.25). Tumor mutation burden and neoantigen load were significantly reduced after CCRT (p < 0.001). We observed no specific correlation with non-synonymous mutations and no changes in the single-nucleotide variant spectrum after CCRT. Post-CCRT samples were enriched in gene sets related to immune signaling pathways, such as interferon gamma signaling and CD28 co-stimulation. Multiplex IHC showed an incremental trend in the proportion of CD4 positive cells in cytokeratin positive region after CCRT. However, CD8, CD20, FOXP1, PD-L1 showed no definitive trend. Proportion of immune cells calculated by CIBERSORT, showed that significant increase in neutrophils after CCRT. Conclusions We have comprehensively analyzed the early immunogenomic changes induced in ESCC by CCRT and correlated them with clinical outcomes. Our results provide a potential basis for combining immunotherapy with CCRT for the treatment of ESCC. Electronic supplementary material The online version of this article (10.1186/s40425-019-0609-x) contains supplementary material, which is available to authorized users.


Background
Esophageal cancer, the sixth leading cause of cancer-related deaths, is classified into two main histological subtypes: squamous cell carcinoma (SCC) and adenocarcinoma. Despite the decreased incidence of esophageal squamous cell carcinoma (ESCC) in Western countries, SCC remains prevalent in Asia, Africa, and South America [1]. The combined modality approach of preoperative concurrent chemoradiation therapy (CCRT) and surgery has shown superior clinical outcomes compared with surgery or chemotherapy alone for the treatment of locally advanced ESCC [2][3][4]. Despite the promising effects of preoperative CCRT, the challenge of cancer relapse after curative resection in patients with ESCC needs to be addressed.
Immune checkpoints are the downregulators of the anti-tumor immune response and their inhibitions are the good treatment strategy as cancer immunotherapy [5,6]. The synergistic effects of CCRT and immunotherapy have been tested in many clinical trials. In a recent trial, patients with locally advanced non-small cell lung cancer that were treated with durvalumab, an inhibitor of programmed cell death-ligand 1 (PD-L1), after CCRT showed prolonged progression-free survival (PFS) compared with patients that only underwent CCRT [7]. However, the underlying immunogenomic changes induced by CCRT in the tumor and its microenvironment remain unknown. Analysis of such changes may offer additional insights into the synergistic effects of immunotherapy and CCRT and also provide potential predictive or prognostic biomarkers [8].

Study flow and patients
Patients with esophageal cancer at clinical stages T1b-T4a and N0 or N+ were evaluated for surgery by a multidisciplinary team that included a medical oncologist, a radiation oncologist, a radiologist, and a thoracic surgeon. The surgical candidates (n = 40) underwent neoadjuvant chemoradiotherapy. Radiation therapy was delivered with a dose of 4300-4400 cGy in daily 210-215 cGy per fraction over 4 to 5 weeks. Concurrent chemotherapy of 5-FU (4000 mg/m 2 over 4 days) and cisplatin (60 mg/m 2 on day 1) were administered every 3 weeks for up to two cycles during RT. Endoscopic biopsies were performed at the time of initial diagnosis and 2 to 3 weeks after the start of preoperative CCRT (Fig. 1a). This genomic analysis study was approved by the institutional review board of the Samsung Medical Center (IRB no. SMC-2013-10-112) and written informed consent was obtained from all enrolled patients.

Isolation of genomic DNA and RNA
Genomic DNA (gDNA) and RNA were purified from cancer tissues using the AllPrep DNA/RNA Mini Kit (Qiagen, USA). gDNA from peripheral blood was extracted using the QIAamp DNA Blood Mini Kit (Qiagen). gDNA concentration and purity were measured using a NanoDrop 8000 UV-Vis Spectrometer (Thermo Scientific Inc., USA) and a Qubit 2.0 Fluorometer (Life Technologies Inc., USA). To measure DNA degradation, median DNA size and ΔCt values were measured using a 2200 TapeStation Instrument (Agilent Technologies, USA) and real-time polymerase chain reaction (PCR; Agilent Technologies), respectively. RNA concentration and purity were measured using the NanoDrop and Bioanalyzer (Agilent Technologies).

Whole-exome sequencing
High quality gDNA in each sample was sheared with an S220 ultra-sonicator (Covaris, USA) and used to construct a library with the SureSelect XT Human All Exon v5 and SureSelect XT reagent kit, HSQ (Agilent Technologies), according to the manufacturer's protocol. This kit is designed to enrich 335,756 exons from 21,058 genes, covering~71 Mb of the human genome. Enriched exome libraries were multiplexed and sequenced on the HiSeq 2500 platform (Illumina, USA). A paired-end DNA sequencing library was prepared via gDNA shearing, end-repair, A-tailing, paired-end adaptor ligation, and amplification. After hybridizing the library with bait sequences for 16 h, the captured library was purified and amplified with an indexing barcode tag, and library quality and quantity were assessed using a 2200 TapeStation Instrument and Qubit 2.0 Fluorometer, respectively. The exome library was sequenced using the 10-bp paired-end mode of the TruSeq Rapid PE Cluster kit and the TruSeq Rapid SBS kit (Illumina).

RNA-Seq data analysis
The reads from the FASTQ files were mapped against the hg19 human reference genome using TopHat version 2.0.6 (http://ccb.jhu.edu/software/tophat/index.shtml). Raw read counts mapped to genes were measured using the BAM format file by HTSeq version 0.6. 1 [30]. A total of 18,161 coding genes were analyzed for transcript abundance and poorly expressed genes were eliminated based on the criteria of a maximum read count > 20 for all samples. Read counts were normalized using the Trimmed Mean of M-values normalization method. Differentially expressed genes were identified using the DESeq R package (www.huber.embl.de/users/anders/DESeq/). A gene set enrichment analysis (GSEA) [31] was conducted to analyze up-or down-regulated gene sets in specific groups of ESCC samples. Stromal and immune scores based on WTS were calculated using ESTIMATE [24]. Fractions of immune-associated cell types were calculated by CIBER-SORT using RNA-seq expression profiles [23]. The immune cytolytic activity (CYT) was measured by the geometric mean of GZMA and PRF1 expression values in TPM [13].

Cancer cell fraction measurement
Cancer cell fractions (CCF) were measured by PyClone, which de-convolves tumor sequences into sub-clones based on a hierarchical Bayesian clustering model [32]. Input data were generated from somatic single-nucleotide variants (SNVs) detected by MuTect and corresponding copy number variations. SNVs were grouped into clusters with similar CCF values.
Tumor mutation burden and prediction of candidate neoantigens TMB was measured by the number of somatic single nucleotide variants and indel mutations per megabase in the coding region [33]. Somatic single nucleotide variants include nonsynonymous as well as synonymous mutations. Non-coding alterations were not counted. In addition, known somatic alterations in COSMIC and truncations in tumor suppressor genes were excluded from the count.
Neoantigens were predicted using MuPeXI v.1.1.3 [34]. The three types of human leukocyte alleles (HLA-A, −B, and -C) were identified from the WTS data of each patient using seq2HLA [35]. Somatic mutations, gene expression counts, HLA types for each patient, and peptide lengths (8-11 mer) were provided as input for MuPeXI. Peptides with a half maximal inhibitory concentration (IC50) value ≤500 nM were considered to have a high binding affinity for the major histocompatibility complex (MHC). Expressed mutant peptide sequences with an IC50 value of ≤500 nM that showed binding affinity above normal were picked as candidate neoantigens.
Slides were scanned using the PerkinElmer Vectra 3.0 Automated Quantitative Pathology Imaging System (Perkin-Elmer, Waltham, MA), and images were analyzed using the inform 2.2 software and TIBCO Spotfire™ (Perkin-Elmer, Waltham, MA). Each cell was identified by detecting nuclear spectral elements (DAPI). All the immune cell populations from each panel were characterized and quantified using the cell segmentation tool by the InForm image analysis software. All cells in each slide were designated as positive or negative for each antibody, and the data were categorized and exported to an xls file for analysis. We used the Spotfire™ program after tissue and cell segmentation and expression intensity was compared and then judged based on the cut-off value. The numbers of CD4, CD20, FOXP3, PD-L1, CD8, and CK positive cells were counted in each slide.

Statistics
The two-sided t-test was used for comparisons of tumor purity, mutation burden, immune and stromal score, cytolytic score, fraction of immune cell from the pre-and post-CCRT samples. The Kaplan-Meier curves and log-rank test was used for the survival analysis. P-value less than 0.05 was considered as a statistically significant.

Baseline demographics
From a total of 40 study participants, 29 and 23 patients participated in paired DNA analysis and whole-transcriptome analysis, respectively (Fig. 1b). The pre-treatment clinical stages were IVA (n = 4), III (n = 21), and II (n = 4). Of the 29 patients that enrolled for paired genomic analysis, four patients did not undergo surgery due to disease progression (n = 2), refusal of surgery, and failure to follow-up. Surgical samples after CCRT (n = 25) showed post-neoadjuvant stages of IVA (n = 2), IIIB (n = 3), IIIA (n = 7), II (n = 1), I (n = 7) and 5 cases of pathologic complete response (pCR). The first biopsy was conducted at a median time of 4 days (range 1-14) before CCRT and the second biopsy was conducted at a median of 18 days (range 4-24) after CCRT. The resection margins were negative for all samples. Seventeen patients are currently enrolled in other clinical trial, which aims to study the effects of randomized adjuvant durvalumab treatment versus placebo (NCT02520453; Additional file 1: Table S1) [36].
Changes in the somatic mutation landscape and copy number alteration in pre-and post-CCRT samples We compared the somatic mutation landscape between pre-and post-CCRT samples, and analyzed representative genes related to cell cycle, histone modification, Hippo pathway, Notch pathway, and the PIK3CA pathway (Fig. 2a). Twenty-four pre-CCRT samples had missense, nonsense, or splicing mutations in the tumor suppressor gene, TP53, and these alterations were maintained in 11 samples after CCRT. Although nine pre-CCRT samples had a missense mutation in the nuclear factor erythroid 2 like 2 gene, NFE2L2, only two post-CCRT samples retained that mutation. SNV analysis showed that transition mutations, specifically C to T or G to A, were prominent in all samples. Comparison of copy number amplification and deletion regions between genomes of pre-CCRT and post-CCRT samples (Fig. 2b), showed that copy number at chromosome 7p14.1 was amplified only in post-CCRT samples. Additionally, the region harboring the cyclin D1 (CCND1) oncogene was amplified and the region harboring the tumor suppressor genes, Cyclin Dependent Kinase Inhibitor 2A (CDKN2A) and CDKN2B, was deleted in both pre-and post-CCRT samples.
Analysis of TMB, immune and stromal score profile, and immune cell composition in the tumor microenvironment in pre-and post-CCRT samples We found no difference in tumor purity, between pre-and post-CCRT samples (Fig. 2c). This allowed us to analyze the TMB and neoantigen load using WES and WTS data. We found that the TMB and neoantigen load were significantly lower in post-CCRT samples (p < 0.001) compared with pre-CCRT samples (Fig. 2d). Using ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data), we found that 14 post-CCRT samples (60.9%) showed an increased immune score compared with pre-CCRT samples. Of those 14 post-CCRT samples, 10 also showed a concurrent increase in the stromal score (Fig. 3a). Therefore, CCRT leads to increased immune and stromal scores in the tumor tissue (Figs. 3b and c). No difference in immune CYT was observed between pre-and post-CCRT samples (Fig. 3d). Using CIBERSORT, we analyzed the changes in the population of immune cell types before and after CCRT. We found that the numbers of resting natural killer (NK)/T-cells (p < 0.001), follicular helper T-cells (p = 0.049), and regulatory T-cells (p = 0.038) significantly decreased after CCRT. However, the population of neutrophils was significantly increased in the post-CCRT samples (p = 0.001; Fig. 3e). We found no difference in the abundance of major immune cells such as B-cells, CD4 T-cells, CD8 T-cells, and M0, M1, and M2 macrophages between pre-and post-CCRT cells.

Multiplex immunohistochemistry results in pre-and post-CCRT samples
We found a strong correlation between outcomes from gene expression profile and multiplex IHC. Multiplex IHC for CD8 was significantly correlated with CYT score (p = 0.012). Especially PRF1 (Perforin 1) expression among CYT genes was significant (p = 0.001) while GRZA (Granzyme A) was not significant (p = 0.105). IFN-r was also correlated with CD8 (p = 0.104), and PD-L1 expression was negatively correlated with CD8 (p = 0.022) (Fig. 4a). Multiplex IHC is conducted in four patients available for both pre-and post-CCRT samples (Fig. 4b, Additional file 2: Table S2). Interestingly, all the samples showed an incremental trend in CD4 cell proportion after the CCRT. The proportion of cell expressing CD8, CD20, FOXP3, PD-L1 showed no definitive trend. (Fig. 4c). TMB was also highly correlated with CD8, but it was not significant (r = 0.69, p = 0.059) (data not shown).

Tissue characteristics and survival analysis of patients with an increased immune score after CCRT
The 14 post-CCRT samples that showed increased immune scores also showed an increased abundance of neutrophils and a decrease in the number of follicular helper T-cells, regulatory T-cells, and resting NK/T-cells compared with the baseline pre-CCRT samples. We found that pre-CCRT samples with a higher proportion of M0 macrophages and lower resting mast cells were likely to show an increase in immune score after CCRT (Additional file 3: Table S3).
Survival analysis showed extended PFS and overall survival (OS) in patient subgroups with increased immune scores (12-month PFS rate: 83.1% vs. 57.1%, p = 0.248; 12-month OS rate: 92.3% vs. 85.7%, p = 0.702) compared with the remaining subgroups ( Fig. 5a and b). To validate the strength of these potential prognostic factors, we analyzed patient survival using pCR, a well-known prognostic factor for ESCC [37]. Disease free survival (DFS) and OS were extended in patients that showed pCR after CCRT (12-month DFS rate: 100.0% vs. 62.2%, p = 0.210; 12-month OS rate: 100.0% vs. 82.2%, p = 0.465), which was similar to the difference in survival rates analyzed according to the change of immune score (Figs. 5c and d).

Differentially expressed gene analysis and gene set enrichment analysis
Sixty genes that satisfied pre-defined criteria (> 2-fold change and adjusted p < 0.01) were identified by differentially expressed gene analysis in pre-and post-CCRT tissue samples (Fig. 6a and Additional file 4: Table S4). We found that the expression of cell cycle-related tumor suppressor genes, CCND2 and CDKN1A, was increased after CCRT. Gene set enrichment analysis showed that several immune-related gene sets, such as those involved in interferon gamma signaling, cytokine signaling, adaptive immune system, innate immune system, PD1 signaling, T-cell receptor signaling, and CD28 co-stimulation were enriched in post-CCRT samples (Figs. 6b and c).

Correlation of genomic profile with treatment outcome
To analyze the genomic profile of the ESCC tissue samples based on their surgical pathology, we compared pCR (n = 5) and non-pCR (n = 20) subgroups using pre-andpost-CCRT samples. We observed no significant difference in immune and stromal scores between the two subgroups in pre-CCRT samples. However, the pCR subgroup showed significantly lower immune CYT compared with the non-pCR subgroup (p = 0.011) in pre-CCRT samples. We found a prominent decrease in CYT in the non-pCR group (p = 0.047) after CCRT, but no change in CYT in the pCR group compared with pre-CCRT samples (Fig. 7a). We also did not find any predictive somatic mutation markers in the pCR subgroup. NFE2L2, a gene frequently found to be mutated in our pre-CCRT samples, mostly disappeared after CCRT (Fig. 3a). CCF estimation analysis showed that NFE2L2 p.D15E was the only unique variant in pCR samples. However, the proportion of cancer cells with specific mutation failed to specifically correlate with the pathologic response (Fig. 7b). Analysis of changes in the population of immune cell types showed a significant increase in activated dendritic cells after CCRT in the pCR subgroup compared with pre-CCRT samples. Neutrophils were significantly increased in number in both pCR and non-pCR subgroups after CCRT compared with pre-CCRT samples ( Fig. 7c and Additional file 5: Table S5). However, analysis of individual gene expression patterns revealed no specific pattern in the pCR and non-pCR subgroups before and after CCRT (Additional file 6: Table S6 and Additional file 7: Fig. S1).
In addition, we looked into the difference in immune profiles between samples from inoperable patients (n = 2) and the patients who received surgery (n = 25). In pre-CCRT samples from inoperable patients showed significantly lower stromal and CYT scores compared with patients who received surgery (P = 0.0014 and P = 0.012, respectively). A similar pattern in immune scores was shown despite not significant (Additional file 8: Fig. S2). Looking into the immune cell compositions between two groups, a fraction of resting NK cells and regulatory T cells in pre-CCRT samples were significantly higher in inoperable patients compared to the patients who conducted surgery (P = 0.037 and P = 0.005, respectively). There was no significant difference in TMB and neoantigen load between both groups.

Discussion
Recent advances in modern sequencing and analytical tool-kits have resulted in the evolution of immunogenomics as a critical component of cancer immunotherapy. Immunogenomic studies help to understand the mechanisms that control therapy response and resistance to immune response shaped by the tumor and its microenvironment. In this study, we have comprehensively analyzed the changes in the immunogenomic profile of ESCC in response to CCRT.
Using our sequencing dataset, we analyzed the genomic profiles of patients with ESCC to identify subgroups that could benefit from CCRT. We found that patient subgroups with increased immune scores after CCRT showed favorable survival outcomes. However, the elevated immune score could be primarily due to an Previous studies [38][39][40] have suggested that the molecular signature of NFE2L2 (NRF2), a master transcriptional regulator of stress response, serves as a predictive marker for esophageal tumor response to CCRT. A gain-of-function NRF2 mutation confers resistance to therapy in ESCC cells [39]. However, our results showed no significant correlation between NFE2L2 missense mutations and the pCR of the study population (pCR rate of 14.3% in patients with NFE2L2 mutations vs. 20.0% in the remaining patients). It is possible that the missense mutations observed in our study population failed to affect NFE2L2 function and further analysis is needed to study the functional alteration of the mutated forms of NFE2L2 in the ESCC tissue samples.
The tumor microenvironment contains diverse immune cell types in addition to tumors cells and its nature and composition change over time with treatment [41,42]. From our previous study, it is known that PD-L1 expression elevated in ESCC samples that received preoperative CCRT compared to the CCRT naïve sample [43]. These findings reflect that CCRT induces immune checkpoint protein expression in tumor which we could also expect alteration in immune cell composition in tumor microenvironment lead by upregulation of PD-L1. We found an enrichment of neutrophils in the tumor microenvironment after CCRT (Fig. 3e). However, we failed to see an increase in the number of adaptive immune cells needed for anti-tumor immune response, such as activated CD8 T-cells and dendritic cells, after CCRT. CYT, which is A B C Fig. 7 Comparison of changes in the immune, stromal, and immune cytolytic activity scores and immune cell fractions between pre-CCRT and post-CCRT samples based on their surgical pathologic profile. a, Immune and stromal scores, and immune cytolytic (CYT) activity in pathologic complete response (pCR), non-pCR, and no-operations (Op) samples. b, Somatic mutation changes in NFE2L2 using cancer cell fraction measurements based on the surgical pathologic profile.c, Changes in immune cell fractions based on the surgical pathologic profile indicative of activated CD8 T-cells, was also not elevated after CCRT (Fig. 3d). Radiation is known to induce an inflammation response by damaging tumor endothelial cells and triggering inflammatory cytokine signaling (via Interleukin 1 and the tumor necrosis factor) [44,45]. Radiation also induces the recruitment of circulating immune cells and increases antigen exposure to initiate an adaptive immune response [46]. Our results on neutrophil enrichment after CCRT suggested that chemoradiation, unlike radiation only, induce non-specific inflammation rather than an adaptive immune response in the tumor microenvironment of ESCC.
High TMB and neoantigen load in a tumor can generate T-cell responses that recognize and eradicate tumor cells [47]. Clinical trials have shown that high TMB increases the efficacy of immune checkpoint blockades in cancer immunotherapy [16]. We found that TMB and neoantigen load were significantly reduced (p < 0.001) in ESCC samples after CCRT. Radiotherapy is known to induce antigen presentation [46]. However, our results showed that underlying tumor mutation burden, which hypothetically shows positive correlation to tumor antigen presentation, were lowered after chemoradiation.
Based on our results that CCRT decreases the TMB and induces non-specific inflammation in ESCC cells, it is possible that the immune cell priming and reinvigoration induced by immune checkpoint inhibitors can have a greater impact when immunotherapy is combined before or at the time of initiating CCRT. Patients with non-small cell lung cancer (NSCLC) that undergo consolidation therapy with anti-PD-L1 inhibitor (durvalumab) following completion of CCRT show increased survival [7]. Concurrent durvalumab treatment with CCRT, starting with CCRT not after CCRT, is currently being tested for improved synergistic effect and survival benefit compared to previous result in same NSCLC population (NCT03519971) [48]. Also, current clinical trials are testing combination therapies in which immunotherapy is initiated one or two weeks prior to CCRT [49][50][51]. Results from these clinical trials may validate our hypothesis.
One of the challenges that we faced during this study was the selection of an optimal timepoint for the second biopsy to best identify the immunogenomic changes caused by CCRT. We had initially designed the study to compare tissue samples from patients with ESCC at the time of diagnosis and post-surgical samples from patients that had undergone CCRT. However, we were unable to perform genomic analyses on samples that were exposed to the entire CCRT regimen due to extensive tumor necrosis. Therefore, we conducted the second biopsy within 2 to 3 weeks after the initiation of CCRT. Due to the reason, small number of samples were available for the multiplex IHC which weakens the representative value of its result and showed inconsistent result in some result such as trend of CD4-cells before and after the CCRT. Future studies should confirm whether our timepoints were optimal for the evaluation of immunogenomic changes in the tumor and its microenvironment.
In conclusion, our study is the first to demonstrate the underlying genomic changes caused by CCRT. It will provide a basis for further genomic studies in patients that undergo CCRT and guide clinical trials that test combination therapy treatments of immunotherapy and CCRT.