Multi-omics profiling of primary small cell carcinoma of the esophagus reveals RB1 disruption and additional molecular subtypes

Primary small cell carcinoma of the esophagus (PSCCE) is a lethal neuroendocrine carcinoma. Previous studies proposed a genetic similarity between PSCCE and esophageal squamous cell carcinoma (ESCC) but provided little evidence for differences in clinical course and neuroendocrine differentiation. We perform whole-exome sequencing, RNA sequencing and immunohistochemistry profiling on 46 PSCCE cases. Integrated analyses enable the discovery of multiple mechanisms of RB1 disruption in 98% (45/46) of cases. The transcriptomic landscape of PSCCE closely resembles small cell lung cancer (SCLC) but differs from ESCC or esophageal adenocarcinoma (EAC). Distinct gene expression patterns regulated by ASCL1 and NEUROD1 define two molecular subtypes, PSCCE-A and PSCCE-N, which are highly similar to SCLC subtypes. A T cell excluded phenotype is widely observed in PSCCE. In conclusion, PSCCE has genomic alterations, transcriptome features and molecular subtyping highly similar to SCLC but distinct from ESCC or EAC. These observations are relevant to oncogenesis mechanisms and therapeutic vulnerability.

number of mutations or just of single nucleotide substitutions and not indels, for example. 3. Similarly, in L167-169, the reader is left to wonder how many TP53 were in total, assuming 52 stands for all nonsynonymous mutations. 4. L 107-109 "The mutation rate of PSCCE in this study was significantly higher than that reported by Wang et al12 (P = .0225) ,owing to the higher sequencing depth." -What was Wang mutation rate? This claim should be supported by subsampling the reads to similar sequencing depth and repeating the analysis. 5. Several methods are not cited appropriately. No reference for MutsigCV, MuTect2, CNVkit, mutational signatures (L536), RSEM, DESEQ, Consensus clustering and more. 6. L162-169 described the genes identified as significantly mutated. The L170-176 list other genes that were mutated -how were those identified? If by manual inspection, how many genes were manually inspected? What was the selection criteria? How many of those were mutated? 7. L182 -homozygous deletion of a gene is by definition mutually exclusive with a mutation in that gene. Possibly partial deletion? 8. L223,262 -differential expression analysis of PSCCE against what? Matched normal samples? Of which tissue? Other cancer types? 9. L232 seed should be see 10. Figure 3c -labels seem to be misplaced. 11. Figure 3g.g -what is the order of the genes? 12. L278 and elsewhere -consider replacing changeable by unstable. 13. L289 'integrated analysis revealed' -what type of analysis? How significant is the association? 14. L293, 294 -what test and significance level were used? 15. L308 -consistent with which of the findings described before? 16. L322, 324 -significant by which test? P-value? 17. L334-335 -where is this finding shown? 18. Figure 5d -what do N and T stand for? 19. L357 Projection on what? The PSCCE RNA-seq data? 20. L362 Deconvolution does not confirm signature projection, merely repeat the anlysis -those are two computational ways to show the same thing based on the same data. 21. L388 filed should have probably be field 22. L479 ethnic should probably be ethic 23. L536-L539 -Consider replacing the last two sentences of the paragraph. 24. L545 -RNA sequencing description is not detailed enough -read length and paired or non paired reads should be specified. 25. L548 -what makes a sequence read validated? How many of the reads were not validated? 26. In several places (L594, L516, L609 and elsewhere), p-value threshold is stated as 0.05, and though multiple comparisns were done, no correction is mentioned. State type of correction if performed, or perform correction if not performed. 27. L608 what filtering was done prior to differential expression analysis?
We greatly appreciate the reviewers' time and constructive comments. To address the concerns of the reviewers, we performed additional analyses and validations, and revised our manuscript accordingly. We believe that the inspiring comments from the reviewers helped us substantially improve our manuscript.
Please find our responses below: Reviewer #1 (Remarks to the Author):

Review
Li and colleagues present a multi-omic analysis of small cell cancers of the esophagus. This descriptive study demonstrates that these tumors are more similar to small cancer of the lung than to EAC or ESCC-which is not a surprising conclusion yet still good to demonstrate. While this will ultimately make a useful contribution to the literature, some improvements are needed both in writing and in some of the analysis and discussion.
1. Of note, they contrast a prior paper with exomes of PSCCE-finding 80% p53 but only 27% RB1…. This is contrasted with their much higher assessment of RB1 loss-but from multiple ways of assessing-it would be useful to clearly compare the rates of somatic mutation with this prior paper.
Response: 2. There are some concerns regarding the quality of the copy number data. I have concerns that many of the most significant peaks are spurious-this may specifically be a problem with joint analysis of FFPE and fresh samples as copy-number is challenging with FFPE samples…. Are any peaks selective to the FFPE samples? There are not enough details on how these data are processed and what methods are used to properly clean the results. Many focal peaks are present but not discussed-the fear is that many are artifacts. Also, they are rather selective in describing the targets of the peaks-for example, they cite MYC as an amplified target but myc is outside the peak listed for 8q24.

Response:
We appreciated very much for this comment. CNV calling from FFPE WES is challenging, due to non-uniform distribution of coverage and fragmentation of FFPE DNA (Ref 1).
Following reviewer's suggestion, we thoroughly checked and revised our CNV analysis pipeline, data cleaning process and results. In the renewed CNV analysis pipeline, we cleaned our results mainly by the following measures: (1) PCR duplicates were marked with Picard (https://broadinstitute.github.io/picard/, 2.9.0) and removed from further analysis. This step was also performed but not mentioned in the initial Manuscript. This information was now added to method section.
(2) Two piled-up references, namely fresh-frozen reference and FFPE reference, were built from fresh-frozen normal samples and FFPE normal samples, respectively. Piled-up reference guarded against random fluctuation introduced during library preparation and sequencing process. Piled-up FFPE reference would also help to reduce non-uniform distribution of sequencing coverage in normal samples.
(3) CNV analysis was performed separately for fresh-frozen and FFPE samples. Specifically, fresh-frozen tumors were compared against fresh-frozen reference; FFPE tumors were compared against FFPE reference.
(4) Bins (each bin equals a single capture region of WES) whose log2 coverage ratio against reference < -15.0 against reference were highly likely due to poor capture or alignment (Ref 2). These segments were excluded from further analysis.
(5) As high coverage depth (> 300×) was achieved in present study, a justified resolution to one bin could by achieved (Ref 2). We allowed report of segments containing as few as one bin (compared to ≥ 5 bins in old version results) to decipher small SCNVs.
(6) Segments were filtered for CNVs that were commonly observed in healthy population. We built a reference set of CNVs in healthy population (DGV.refset) from Database of Genomic (7) The default thresholds (also the thresholds used in initial Manuscript) for calling a SCNV in CNVkit were relatively loose, mainly to guarantee sensitivity in low purity tumor samples.
With the aforementioned measures applied, we observed the following feature in the revised  c) We discovered more RB1 "exon deletions" by the revised SCNV calling pipeline. We found homozygous deletion events often affected only some exons (median: 11 exons, range: 1 to 25) but not full gene body of RB1. We use "exon deletions" to replace the imprecise "homozygous deletions of RB1" used in initial Manuscript.  We followed the reviewer's comment to search selective peaks in FFPE samples. Numerous GISTIC peaks appeared selective to FFPE samples ( Fig R2). For example, peak 8q24 harboring MYC were called from FFPE group but not observed in fresh-frozen group. However, MYC was also amplified in 3 fresh-frozen tumors. We propose SCNVs at 8q24 in the relatively small fresh-frozen group (n = 13) might fail to get statistically significant simply due to limited observation numbers.
Similarly, 13q14 peak harboring RB1 deletion were not called in fresh-frozen group, but RB1 deletions were observed in 6 fresh-frozen tumors. Amplification peak on 6p22 was selective to FFPE samples.
Genes in this peak, including E2F3 and SOX4, were amplified in none of 13 fresh-frozen samples. Indels affecting EP300 (n = 1) and CREBBP (n = 1) have not been previously described (COSMIC database v91, searched at 2020-05-06). We validated them using Sanger sequencing ( Fig R3). 4. They claim homozygous Rb1 deletion in 8 cases. The GISTIC results show only a borderline significant event-if it was a homozygous event in ~20% of tumors I would expect a much more significant result. How is 'homozygous' deletion called?

Response:
As described in response to comment 2, a segment with a log2 ratio against normal sample reference < -2.0 (equals to < 0.5 copy in a diploid genome) was considered "homozygous deletion".
Genes whose partial or all were involved in a deleted segment were then reported as "homozygous deletion".
RB1 exon deletions in PSCCE were short in length (median: 52.5kb, range: 0.2kb to 2828kb; described in detail in response to comment 2). These short exon deletions might have been partially masked by improper noisy-reducing parameter (joint segment size) setting in initial Manuscript.
We re-set "joint segment size" to its default value of zero to include small SCNVs. With the revised SCNV results and updated GISTIC parameter, fourteen RB1 deletions lead to a sharp 13q14 peak with q value of 6.91× 10 (Fig R2 in response to comment 2).
In the revised Manuscript, we updated the parameters in GISTIC analysis and updated results presented in Fig 1d and Supplementary Table 3. 5. In Figure 2a-how do you explain a homozygous deletion and a splicing abnormality in the same tumors? The authors need to evaluate a collection of 'normal' tissues some and ideally other cancer RNAseq data with their analytics to determine how often these splice forms are seen.

Response:
RB1 has an extraordinarily large gene body, with small exons separated by large introns. As was described in response to comment 2, most of the deletions observed in RB1 gene deleted only part of but not whole RB1 locus ("exon deletions"). The remaining exons of RB1 could still be transcribed into mRNA and joint together to produce splicing abnormalities observed in RNA-seq. Fig R4 gives two examples of exon deletions as genomic basis for splicing abnormality. We used RNA-seq data from 23 match normal esophageal samples sequenced together as normal control in our initial Manuscript. Here, following reviewer's advice, we provided Sashimi plots of 23 matched normal esophageal tissue sequenced in our study ( Fig R5). Additionally, we provided sashimi plots of 10 normal lung sample RNA-sequenced by our lab (unpublished data, Fig R6).
As was shown in Fig  Thus we confidently conclude that the splicing abnormalities described in PSCCE tumors were cancer-specific events. 6. Negative IHC for Rb1 is a difficult marker-can they benchmark their results in small cell cancers vs. other tumors with a larger dataset? Also, these results need to be read in a blinded fashion (e.g. a pathologist needs to read them without knowing which are small cell vs other) and verified by 2nd blinded reader.

Response:
Negative IHC marker bears the risk of producing false negative when IHC procedure fails. In awareness of this, the IHC of Rb1 showed as part of results in initial Manuscript was performed on a tissue microarray (TMA), with matched normal esophageal tissues alongside the corresponding tumor samples ( Fig R7). By checking positive Rb1 signal in nuclei of normal samples, we could confidently exclude an IHC failure. Results of Rb1 IHC in benchmarking cohort were summarized in Table R3 below. The proportion of Rb1 negative SCLC was highly similar to PSCCE. The proportion of Rb1 negative tumors in PSCCE were substantially higher than ESCC and LCNEC. Based on the above evidences, we could confidently conclude that our Rb1 IHC evaluation pipeline is trustworthy. We added the Rb1 IHC of matched normal samples in revised 7. Line 251-the statement that the 3 groups showed differential 'dependence' on oncogenic pathways is not correct as stated-dependence needs to be shown functionally and cannot be determined from gene set analyses Response: Thank you for pointing this out. We replace L251 in original Manuscript with the following "The three groups also showed differential activation of oncogenic signaling pathways".
8. 5a as a data slide on the notch pathway is a schematic but does not show real data to support this point.

Response:
We appreciate this comment from Reviewer 1. To better present dysregulation of Notch pathway  9. In the methods-it states that both fresh and FFPE-in the text of the paper, this breakdown is necessary to detail in more depth in terms of the breakdown of types. Also, there needs to be evaluation for possible systematic biases in results of fresh vs. FFPE results.

Response:
Following reviewer's suggestion, we evaluated differences between fresh-frozen and FFPE DNA results in the following aspects: (1) The proportion of C>T in FFPE samples (53.1%) was higher than that in fresh-frozen samples  (4) As described in response to comment 2, we found by performing data cleaning and applying stringent threshold to pick out most prominent amplification and deletions. SCNV results were highly consistent with qPCR validation in both fresh-frozen and FFPE samples.
"Selective" peaks in FFPE samples were mainly due to larger sample size but not systemic bias.
In summary, we did find certain differences between FFPE sample and fresh-frozen sample, all had been previously reported. Some influence was introduced by FFPE samples, such as mutational burden and enrichment of C>T substitution in mutational signature analysis. However, as discussed above, these influences were mild.
Sequencing and subsequent analysis of FFPE samples are challenging. For a rare cancer such as PSCCE, archival FFPE samples were probably the only choice for any sizable sequencing study. FFPE samples were proved to have good performance in discovering driver mutations and actionable events (Ref 1-3). We also found that by stringent filtering and proper interpretation, FFPE WES sequencing provided trustworthy and valuable information on molecular characteristics of PSCCE.
In the revised Manuscript, we discussed the above aspects in corresponding sections. We also provided sample type information in Supplementary Table 1 for readers to comparing parameters of interest between two sample types. 10. The writing needs correction/editing for more proper English grammar.

Response:
We appreciate this comment from the Reviewer. Our revised Manuscript has been edited by American Journal Experts for more proper English grammar.
Reviewer #2 (Remarks to the Author): Li et al. perform a molecular characterization, at a substantially deeper genomic level, and for the first time at the transcriptomic level, of esophageal small cell carcinomas (PSCCE). These tumors present similarities to small cell lung carcinomas (SCLC), including ubiquitous loss of TP53 and RB1 loss, alterations in NOTCH family members and in epigenetic regulators, common molecular subtypes and strong immune suppression. These data are novel, informative, and of potentially high impact, suggesting that the molecular features described in SCLCs may be extendable to other neuroendocrine small cell carcinomas. The analysis appears to be well-performed, and the text requires only minor grammatical editing.
I have a few comments and suggestions regarding the manuscript.
Most importantly, the authors should clearly state the statistical tests used for each analysis in the text, and especially in the figures. As one example, please state the statistical tests used to compare the differentially altered signaling pathways within the different esophageal cancer histologic types (Figure 3b). It is also unclear what the comparator is in the DEG analysis of PSCCE (Figure 3a). Is it normal tissue or the tumors from the other esophageal cancer histologic types? Please address throughout.

Response:
We appreciated your comment. We clarified the statistical tests by stating the following in the revised Manuscript (newly added part is underlined): (1) By enrichment analysis, we discovered several pathways that were largely distinguishable between PSCCE and esophageal cancers (Fig 3b). …… (In corresponding Method section) Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment was performed using R package "clusterprofiler". Enrichment with p < 0.01 and Benjamini-Hochberg corrected q value < 0.05 were considered significant. Pathways that was significantly enriched for upregulated DEGs but not for downregulated DEGs were considered upregulated; similarly, pathways which were specifically enriched for downregulated DEGs were considered downregulated.
(2) We compared the gene expressions of 38 tumors against 23 matched normal esophageal samples to identify differentially expressed genes…… (In corresponding Method section) Differentially expressed genes (DEGs) were identified using R package "DESeq2" 62 , by comparing 38 tumors against 23 matched normal esophageal samples sequenced together in the present study.
I would suggest a brief discussion of the potential significance of the enrichment of mutational signature E1 in PSCCE, especially in light of different etiologic drivers from SCLC (less strongly associated with tobacco exposure). What does this say about oncogenesis? Response: Signature E1 is highly similar to COSMIC mutational Signature 1 (v2, 2015), whose etiology is spontaneous deamination of 5-methylcytosine. We propose that the enrichment of signature E1 may provide clue of oncogenesis in the following aspects: (1) Endogenous and spontaneous deamination of 5-methlcytosine is the major mutagenesis mechanism in PSCCE.
(2) The differences of tobacco-related signature may reflect an organ-specific vulnerability to particular mutagens. The esophagus was not directly exposed to tobacco. Although smoking has been defined as risk factor for both ESCC ( It would be interesting to see if the expression of the transcription factors defining the other two SCLC subtypes, POU2F3 and YAP1, are differentially expressed among the PSCCE samples. Figure  4A seems to show some samples with low ASCL1 and NEUROD1 expression and it would be interesting to comment if these may express POU2F3 or YAP1. I am sure the authors have looked at this, and if the answer is no, this could just be very briefly stated.

Response:
We appreciate this comment from Dr. Rudin. Whether a POU2F3 or YAP1-dominated subtype existed in PSCCE was also of our concern.
We identified 5 sample lower for both ASCL1 and NEUROD1 (both < 10.0 TPM), which were listed in Table R4. POU2F3 in these five tumors were also very low. YAP1 expressions values were slightly higher. We found that it might be due to a higher organ-specific baseline of YAP1. We compared YAP1 and POU2F3 expression to ASCL1 and NEUROD1, and compared the ratios with SCLC-P and SCLC-Y samples. SCLC-Y subtype had both high YAP1/ASCL1 ratios and YAP1/NEUROD1 ratios (Fig R10). However, we did not observe any PSCCE had such outstanding expression of YAP1 (that is, 38 PSCCEs in the present study had relatively high level of either ASCL1 or NEUROD1). Moreover, YAP1-dominated subtype was reported to be associated with intact Rb protein, while Rb was completely abolished in all 38 RNA-sequenced tumors. We propose that the evidences obtained from the present study was insufficient to claim a YAP1-regulated subtype. "POU2F3 and YAP1 levels in the 38 RNA-sequenced PSCCEs were relatively low compared to ASCL1 and NEUROD1 and did not appear to be a selective master regulator (Supplementary Fig 5).
Given the ASCL1-or NEUROD1-regulated expression patterns observed in our 38 RNA-sequenced PSCCEs, the evidence obtained in the present study was insufficient to confirm a POU2F3-or YAP1-dominated subtype." We also discussed the potential cause of the paucity of POU2F3-or YAP1-dominated subtype in the Discussion section, by stating: "Although we observed shared ASCL1-and NEUROD1-regulated subtypes, we did not confirm any In SCLC, NEUROD1 high tumors exhibit milder neuroendocrine features than ASCL1 high counterparts (while still maintaining expression). Is this observation reproduced in PSCCE? Response:

Yes, this observation was reproduced in PSCCE.
We observed that PSCCE-N subtype had significantly lower level of neuroendocrine marker DDC and GPR than PSCCE-A subtype (Fig R11a, Wilcoxon  In terms of neuroendocrine markers widely used for diagnosis of both SCLC and PSCCE (Fig R11b), we found PSCCE-A subtype had significantly higher level of NCAM1 (also known as CD56), while CHGA (Chromagranin A) and SYP (synaptophysin) were comparable between two subtypes.

Response:
We appreciate this suggestion. We provide a new Supplementary table 13 same with Table R5 below, comparing clinical characteristics between two subtypes.
Major clinical characteristics, including age of diagnosis, gender, smoking, alcohol drinking, TNM stages and adjuvant therapies, showed no significant difference between two subtypes of PSCCE. Does multivariate analysis suggest that PSCCE subtype and NOTCH mutations are independent prognostic factors?

Response:
No, neither molecular subtype nor NOTCH status was independent prognostic factor by multivariate analysis in our cohort of 38 patients with RNA-seq data.
In univariate analysis (Table R6), we discovered that male gender, smoking, AJCC TNM 8 th stage IV, no adjuvant chemotherapy and subtype of PSCCE-N were significantly associated with poorer prognosis in 38 PSCCE with RNA-seq data. NOTCH1-4 somatic mutation (SNVs and indels) was no longer associated with overall survival in RNA-seq cohort, perhaps due to a further limited sample size.  2. Clustering which puts samples from different datasets in different clusters is very likely to be due to the inherent difficulty in combining datasets. Quantile normalization does not remove all technical differences. This is not discussed at all. There are ways to convince the reader this is not merely a technical effect, e.g. that normal samples of the same tissue from different datasets cluster together, or that tumors of the same type from different datasets cluster together. The fact that the Wang PSCCE data is not included is also strange.

Response:
We agree with the Reviewer. Cross-study comparison and clustering of RNA-seq data are challenging. Technical differences could be introduced in nearly every step of RNA sequencing, including RNA extraction, library preparation, sequencing platform and post-sequencing analysis (Ref  When preparing gene expression matrix, we found huge heterogeneity in "normal esophageal samples". We discovered that the majority of normal samples in TCGA esophageal study (Ref 5) were actually gastric mucosa, as they highly express gastric marker progastricsin PGC but lack squamous cell markers KRT14 and TP63 (Table R8). Thus, in an ideally batch-effect-free scenario, only a small fraction of TCGA-ESCA normal samples would cluster with our normal samples. We corrected expression matrix using quantile normalization, batch-effect correction algorithm ComBat and BUScorrect, respectively. We evaluate differences between datasets in each corrected matrix by dimension reduction using uniform manifold approximation and projection (UMAP, Ref 6). In brief, when performing batch-effect correction, each study was considered a batch. Gene with average expression in the upper 50% (n=9262) of corrected expression matrix were kept and log2 transformed.  As is shown in Fig R12a, raw expression matrix showed strong batch effect, as samples from each study tightly clustered together and no mixing of esophageal normal samples was observed. Quantile normalized (Fig RX12b), ComBat ( Fig R12c) and BUScorrect (Fig R12d) all reduced batch effect, by clustering normal esophageal samples from TCGA (red dots) and the present study (light blue dots) together. In all scenario, most TCGA ESCA study normal samples (those suspected as gastric) clustered with esophageal adenocarcinoma, perhaps by their shared glandular features.
Then we compared downstream analyses of each corrected matrix. As is shown in Fig R13a, three principal groups were reproduced using all three correction methods. We compared identity of each sample and found 98.3% (351/357) and 98.0% (350/357) samples pertained their group identity when using ComBat and BUScorrect correction methods, respectively. All 118 sample categorized into Group NE still clustered together as Group NE when using ComBat and BUScorrect correction. Both ComBat and BUScorrect corrected matrices produced similar heatmaps (Fig R13b) with three main groups corresponding to Group NE, Group Adeno and Group Squamous in the Manuscript. We also performed gene set enrichment analysis (GSEA) corrected matrices. We found the differential enrichment of pathways which we mentioned in the Manuscript, could also be reproduced from expression matrices corrected by ComBat and BUScorrect (Fig R14). In summary, we found that (a) quantile normalization did reduce batch-effect (but not fully eliminate for sure); (b) three widely used batch-effect minimization method produced highly similar results; (c) conclusions from quantile normalization could be cross-validated with genomic or pathological features not affected by RNA-seq batch-effect. As was shown and discussed above, we are confident that clustering results shown in Fig 3d was not merely technical differences. If required, we would be happy to provide the full response to comment 2 as a "not merely 3. The author did not identify PSCCE subtypes corresponding to two subtypes of SCLC by unsupervised analysis. Are those the least frequent subtypes? It may be worth to try supervised analysis to see if the marker genes of those subtypes are expressed in some tumors.

Response:
Yes, the other two subtypes, SCLC-P and SCLC-Y, were less frequently observed. SCLC-P subtype was estimated to comprise 12%-18% of all SCLC (Ref 1 We added a paragraph discussing expressions of POU2F3 and YAP1 to the results section, concluding "POU2F3 and YAP1 levels in the 38 RNA-sequenced PSCCEs were relatively low compared to ASCL1 and NEUROD1 and did not appear to be a selective master regulator…the evidence obtained in the present study was insufficient to confirm a POU2F3-or YAP1-dominated subtype". We also  4. Is the NOTCH status different between the subtypes? Response: No, Notch status showed no significant difference between subtypes.
We also checked expression level of Notch receptors (NOTCH1-4, Fig R16a), main effector (HES1), inhibitory effector (HES6) and inhibitory ligands (DLL3 and DLK1, Fig R16b). We observed multiple lines of evidence showing insufficient T cell infiltration in both PSCCE and SCLC, as were shown in initial Manuscript. Some signatures reported in initial Manuscript were signatures for lymphocytes (including both B cells and T cells, "Module 4" and "Module 5" in initial Supplementary Fig 5a). Besides T cells, we also observed insufficient infiltration of macrophages (Fig R17a), which was another commonly observed population in TME. Some populations, such as B cells, appear to show inconsistent results in some analyses (Fig R17b). For the conciseness and preciseness of the Manuscript, we did not discuss them thoroughly.

We chose to focus on T cells as: (1) low T cells and CD8 T cells abundance in PSCCE and
In the revised Manuscript, we clarified that multiple immune populations were checked. We listed all source of signatures in revised Supplementary Table 14. We provided the raw scores for all immune signatures in a new Source Data 4 to give readers clues on immune infiltrates not described in main text. We also discussed limitation of our computational method in Discussion section to remind readers of need for future studies. Impact of organ-of-origin and histology was examined using two-way ANOVA (degree of freedom: Organ-of-origin = 1; histology = 1), with P value and F value shown below corresponding boxplot. Scores between each pair of cancers was compared using Wilcoxon's rank-sum test. * P < 0.05; ** P < 0.01; *** P < 0.001.
Nonsignificant P values were shown on plot. 6. L524 requiring that a gene is expressed in at least half of the tumors to call it expressed for the purpose of candidate driver gene is unnecessarily too strict. A driver gene in 45% of tumors may still be extremely important. 10% is more reasonable, and possibly a higher expression level threshold, as driver genes are likely to be highly expressed in tumors where they are amplified.

Response:
Following the Reviewer's suggestion, we re-defined "expressed" as "having expression levels ≥ 10.0 TPM in at least 10% of 38 tumors with RNA-seq data". Two candidate genes, DNAH9 and ODF1 were filtered by the old expression filter. Their expression levels still not satisfied the new expression filter of "≥ 10.0 TPM in at least 10% tumors" (Table R9). On the other hand, all SMGs and genes with significant mutational cluster reported in initial Manuscript had relatively high expression and were not affected by the new expression filter (Table   R10). We appreciate this suggestion. We replaced "insufficient immune cell infiltration" by "T cell exclusion" in the abstract, the tumor microenvironment section and discussion of the revised Manuscript.
2. In mutation identification analysis, it will be informative to add the total number of mutations, not just the nonsynonymous. Synonymous mutations might also have an effect, especially if identified in the same gene in multiple tumors. The reader is left to wonder if the 4713 in L113 is the total number of mutations or just of single nucleotide substitutions and not indels, for example.

Response:
We appreciate this comment. We updated our description by stating (newly added part is 6. L162-169 described the genes identified as significantly mutated. The L170-176 list other genes that were mutated -how were those identified? If by manual inspection, how many genes were manually inspected? What was the selection criteria? How many of those were mutated? Response: By MutsigCV and mutational cluster analysis, we identified TP53, RB1, NOTCH1, EP300 and FBXW7 as significantly mutated genes. Detailed mutation of TP53 was described in L167-L169. RB1 and Notch pathway suffered dysregulation by multiple mechanism and were described as two independent section. We manually inspected 97 genes that were mutated in at least three tumors and attempted to identify cancer-associated genes (the COSMIC Cancer Gene Census, ref 1), or genes within same oncogenic pathway (according to ref 2). We found 22 cancer genes with established role in cancer.
Histone modifier, including EP300, CREBBP, KMT2D and KDM6A, comprised a large fraction of the 22 genes list. They were also reported to be frequently mutated in both SCLC and esophageal cancers.
We also found frequent mutation of atypical cadherins, including FAT1, FAT3 and FAT4 and presented them as part of Fig 1e. For the conciseness of the manuscript, we did not describe number of mutation in details for each gene. Instead, full mutation list were provided for readers to search for genes of interest.
We updated our description in the revised Manuscript by stating: "We sought genes with 7. L182 -homozygous deletion of a gene is by definition mutually exclusive with a mutation in that gene. Possibly partial deletion?

Response:
Yes, 13 out of 14 RB1 homozygous deletions observed in PSCCE are partial deletion, affecting a median of 11 exons (range: 1 to 25) but not the whole RB1 locus and having a median length of 52.5kb (range: 0.2kb to 2828kb). We used "exon deletions" to replace the expression "homozygous deletion of RB1" in the revised Manuscript. Table R11 below listed all the exon deletions observed in RB1 (more than initially reported because we revised our CNV calling pipeline following comment from Reviewer #1). In differentially expressed gene (DEG) analysis, we compared the 38 tumors against 23 matched normal esophageal samples sequenced together in our study.
We clarified this process in the revised Manuscript, by stating (newly added part is underlined): "We next turned to transcriptomic landscape of PSCCE…We compared gene expressions of 38 PSCCE tumors against 23 match normal esophageal samples to identify significantly differentially expressed genes (DEGs)." The corresponding Method section has also been updated to: "Differentially expressed genes (DEGs) were identified using R package "DESeq2", by comparing 38 tumors against 23 matched normal esophageal samples sequenced together in the present study." 9. L232 seed should be see

Response:
Thank you for pointing this out. We corrected the word in revised Manuscript.

Response:
Yes, these labels belonged to Fig 3f and were misplaced there. We removed these labels in the revised Manuscript.
The whole list of genes are as following: Similarly, we estimated FDR corresponding to P value of each tests described in L609 (identification of genes of PSCCE subtypes) through BH procedure. We required a signature gene to have P value < 0.05 and log 2 FC >= 2 and q value < 0.1. We found in signature genes of PSCCE-A subtype, 6 out of 115 were rejected; in signature genes of PSCCE-N, 3 out of 36 were rejected. The rejected genes were listed in Table R13. None of them was involved in downstream analysis.  Table S2, related to filter comparing variant reads distribution in tumor and normal) of all raw Mutect2 calls before any of the seven filters was applied. We found that, even by very conservative BH procedure estimation, a p value of 0.05 corresponded to an FDR of only 0.0931 (< 0.1). By applying a series of seven filters, FDR in final list of somatic mutations would be well below 0.1.
We clarified the correction processes applied to L594, L609 and L516, and updated the Method section. Benjamini-Hochberg FDR q values were added to Supplementary tables 2, 10 and 12.