Whole genome sequencing analysis identifies recurrent structural alterations in esophageal squamous cell carcinoma

Clinical samples

Tumor and normal samples were obtained from 20 patients in Kindai University hospital, Osaka, Japan. All the patients agreed to participate in the study and provided written informed consent following ICGC guidelines. The study was approved by the Institutional Review Board at Kindai University Hospital and RIKEN (approval number 25-031) and Personal history of ESCC in these patients were unavailable. All except one patients (OK047) had not received any cancer treatment before the sample collection. OK047 received cisplatin-based chemotherapy before sample collection. Tumor tissues in esophagus were collected by biopsy, and histologically confirmed as ESCC. The patients were treated with neoadjuvant chemotherapy after collecting the samples. The clinico-pathological data are available in the Table 1 and Table S1.

Whole genome sequencing

We performed WGS of the 20 pairs of matched tumor and normal samples. The tumor DNA was extracted from the ESCC samples, and normal DNA was from the lymphocytes in blood. The libraries were prepared using TruSeq Nano DNA Library Prep Kit (Illumina) following the manufacturer’s protocol. Paired-end sequencing of 101- or 126-bp reads was performed using HiSeq2000/2500. The Fig. S1 shows the schematic representation of the WGS analysis pipeline performed in this study. Sequence reads were mapped to the human reference genome GRCh37 using BWA. We removed PCR duplicates using Picard tool (http://broadinstitute.github.io/picard/).

Somatic mutation calling and mutation signature profiling

Somatic single nucleotide variations (SNVs) and short insertions/deletions (INDELs) were called as previously described (Fujimoto et al., 2016). Functional annotation of the detected SNVs and INDELs was performed with Annovar (Wang, Li & Hakonarson, 2010). We applied dNdScv method (Martincorena et al., 2017) to search for genes with significant recurrent mutations (q-value < 0.05). We further summarized and visualized the annotated variants using the MAFtools package (Mayakonda et al., 2018) of the R software (https://www.r-project.org/). For the detection of mutational signatures in 20 ESCC, we used the SignatureAnalyzer (https://software.broadinstitute.org/cancer/cga/msp). Identified mutational signatures were compared with the COSMIC mutational signatures (version 2) using the cosine similarity scoring.

Structural variation calling

Somatic SVs were called by merging calls of two software: in-house pipeline (Fujimoto et al., 2016) and Genomon2 structural variation (SV) detection tool (https://github.com/Genomon-Project/GenomonSV). In the Genomon2 SV detection, The SVs were called using minimum junction number 2, maximum control variant read pair 10 and minimum overhang size 50. The SVs were then filtered using parameters minimum allele frequency 0.07, maximum control variant read pair 1, control depth threshold 10 and minimum overhang size 100. The inversion size threshold set to 1,000 and the simple repeats were removed. Here, SVs were categorized into four classes based on the mapping information for a read pair. The four classes are intrachromosomal deletion, inversion and tandem duplication, and interchromosomal translocation, respectively.

The breakage-fusion-bridge (BFB) events were identified based on the information of fold-back inversion and loss of telomere (Hermetz et al., 2014). We implement the following criteria in order to infer BFB: (i) Inversion is single inversion (either forward or reverse) i.e., without any reciprocal partner, (ii) Inversions must have copy number change versus the adjacent position, and (iii) The two ends of the fold-back inversion must be separated by <20 kb. To detect kataegis, we used Maftools (Mayakonda et al., 2018) package, which is defined by mutation clusters of six or more consecutive mutations localized in a small region with an average inter-variant distance of less than or equal to 1 kb. Chromothripsis was inferred in ESCC samples based on the criteria provided by a previous study (Korbel & Campbell, 2013). Chromothripsis was identified in samples which show clustering of SV breakpoints, usually more than 10 breakpoints within 50 kb, with regular oscillation of copy number states.

Copy number alteration calling

Somatic CNAs were called by analyzing read depth of matched tumor and normal using the Varscan2 software (Koboldt et al., 2012). The thresholds used to call the CNAs using CopyNumber function were p-value 0.001, minimum segment size 100 and maximum segment size 1000. The output of the above step was filtered using the copyCaller function. The raw CNAs were segmented by the circular binary segmentation method implemented in the R package DNAcopy. The GISTIC2.0 algorithm was used to identify the significant recurrent copy number amplified and deleted regions (Mermel et al., 2011).

Identification of druggable genes

To identify the druggable genes across the ESCC samples, we used online gene-drug interaction database (DGIdb) (Griffith et al., 2013). We used the genetically altered genes as target data in order to determine the druggability of the genes. The target genes were detected by different method such as SNV/INDEL analysis, SV analysis and CNA analysis in this study. The database provides two options to examine gene-drug interactions either by gene or by drug names. In this case, we identified gene-drug interactions by providing the gene names. The genes that have at least one interaction with drug target was considered as druggable gene.

Whole genome sequencing of ESCC samples

To identify the mutational events and driver genes that contributing to the development of ESCC in Japanese population, we performed WGS of 20 pairs of tumor and matched blood samples. The samples were collected by biopsy from the individuals before neo-adjuvant chemotherapy (except OK047). Among the 20 samples, 10 samples responded well to the therapy while rest 10 samples showed poor response. The average genome coverage was 43. 4 × for the tumor samples and 34. 3 × for blood samples, after removal of polymerase chain reaction (PCR) duplicates. The WGS data was computationally analyzed to call somatic alterations of the following types: single nucleotide variations (SNVs), small insertions and deletions (INDELs), structural variations (SVs), and copy number alterations (CNAs). In all, WGS analysis identified 104,534 somatic SNVs, 10,523 somatic INDELs, and 2,641 somatic SVs in the 20 tumors (Table S2). In addition, WGS analysis detected 10 significant copy number altered regions in ESCC.

Recurrent coding mutations in ESCC

We investigated the somatic SNVs, short INDELs in the protein-coding regions and their splice sites in the 20 ESCCs. Our WGS analysis identified recurrently mutated genes, including previously known esophageal cancer associated oncogenes and tumor-suppressor genes. We evaluated the somatic alterations by dNdSCV (Martincorena et al., 2017) method in order to find out the significantly recurring mutations across the ESCC genomes. We identified significant mutations in TP53 and ZNF750 genes in ESCC consistent with findings by previous studies (Fig. 1 and Fig. S2). TP53 mutations were most frequent and found in 55% of the samples followed by ZNF750 (15%). In addition, other genes that recurrently mutated in ESCC include FAT1 (10%), PTCH1 (10%), EP300 (5%), FAT2 (5%), FBXW7 (5%), KMT2D (5%), NFE2L2 (5%), NOTCH1 (5%), PIK3CA (5%), RB1 (5%), RIPK4 (5%) and TP63 (5%). These genes were previously reported in ESCC by different studies (Zhang et al., 2015; Sawada et al., 2016; Chang et al., 2017; Qin et al., 2016; Li et al., 2018; TCGA, 2017).

The mutational signatures of ESCC

In order to understand mutational mechanisms of ESCC in Japan, we analyzed the mutational signatures of the 20 ESCC tumors. A Bayesian variant of the non-negative matrix-factorization method was applied to trinucleotide substitution patterns and extracted six mutational signatures (Figs. 2A–2F). The mutational signatures in ESCC were then compared with the signatures of the Catalogue of Somatic Mutations in Cancer (COSMIC) database Alexandrov et al., 2013 (Table S3).

The identified Signature W1 was highly similar to the COSMIC signature 18 (cosine similarity 0.933), but the biological aetiology of this signature is not known. Signature W2 was characterized by C>A and C>T mutations and found in almost 85% of the ESCC patients (Fig. 2G). Signature W2 was similar to the COSMIC signature 1 (cosine similarity 0.850), which is an age-dependent signature. Signature W4 mainly represented by T>C mutation, was similar to COSMIC signature 5 and 16 (cosine similarity 0.868 and 0.9, respectively. Signature W5 displayed high similarity with COSMIC signature 2 and 13 (cosine similarity 0.811 and 0.835, respectively), which was characterized by C>T and C>G mutations. COSMIC signatures 2 and 13 were assigned to the hyperactivity of the APOBEC family enzyme, cytidine deaminases (Alexandrov et al., 2013). One of our samples, OK101, was basaloid squamous cell carcinoma of esophagus, a rare type of malignancy, and a hypermutator in our cohort (Fig. 2G). Somatic SNVs of this sample were dominated by the Signature W5. Signature W6 was characterized by CpTpT-to-CpGpT mutations and highly similar to COSMIC signature 17 (cosine similarity 0.945). COSMIC signature 17 is often found in esophagus and stomach cancer, but its aetiology is unknown. Signature W3 which was defined by C>A, C>T and T>C mutations appeared to have low similarity with any of the COSMIC signatures (all cosine similarity <0.7).

Since, previous studies did not find signatures associated with smoking alone in ESCC despite smoking being a major risk factor of this cancer (Sawada et al., 2016; Zhang et al., 2015). We examined the association of the mutation signatures with clinical features in ESCC. We found signature W4 was significantly elevated in the smoking patients compared to non-smoking patients (P = 0.01474, t-test) (Fig. 2H). Mutations of W4 may be caused by the carcinogenic chemical in tobacco smoke. Signature W4 shared similarity with COSMIC signature 5 and 16. Smoking associated mutational signatures were found higher in the squamous cell carcinomas of lung and head & neck (Wang et al., 2019). However, there was no difference in total mutation burden between the smoker and non-smoker group (P = 0.7, t-test) unlike liver cancer. In liver cancer patients, similar signature showed higher mutation rate in smoker group than the nonsmokers (Alexandrov et al., 2016). The other signatures showed no association with smoking, gender, response to chemotherapy and alcohol drinking status (P > 0.05).

Structural variations in ESCC

WGS analysis detected a total of 2,641 SVs in the 20 ESCC samples with an average of 132 SVs per tumor. The number of SVs varied, ranging from 0 to 515 across the 20 ESCC cases which shows the heterogeneous nature of tumor genome (Fig. S3). In particular, compared to other ESCC cases OK007 and OK008 had high number of SVs which affected most of their chromosomes, indicating genomic instability in these samples (Figs. 3A–3L and Fig. S4). The deletions were the most abundant type of SVs across the samples (Fig. 3M). We found 1,090 (41.27%) deletions followed by 793 (30.03%) inversions, 391 (14.80%) translocations and 367 (13.90%) tandem duplications, respectively.

We found numerous cancer associated genes affected by SVs in ESCC (Fig. 3N and Table S4). Consistently with a report on Chinese ESCC (Chang et al., 2017), LRP1B (30%) and TTC28 (30%) were the most commonly affected genes in Japanese ESCC. The tumor suppressor gene LRP1B was mostly affected by deletions, whereas TTC28 was by interchromosomal translocations. We were also able to identify SDK1, a novel gene, affected by SVs in 25% of the ESCC samples. SDK1 is a cell adhesion molecule that plays an active role in cancer development. Somatic mutations in this gene was reported in adrenocortical carcinoma (Juhlin et al., 2015). Recurrent SVs in CSMD1, WWOX, ERC1, PDE4D and SHANK2 were also identified in five tumor samples. The CUB and Sushi multiple domains 1 (CSMD1) is a tumor suppressor gene reported to be associated with poor prognosis in many cancer types including breast cancer, gastric cancer, head and neck squamous cell carcinoma (HNSCC), and hepatocellular carcinoma (Deng et al., 2012; Zhang et al., 2019; Jung et al., 2018). Deletion of WWOX, a tumor suppressor gene, is frequent in esophageal adenocarcinoma (32%) and stomach adenocarcinoma (30.2%) and also observed in other human cancer types such as colon adenocarcinoma, bladder urothelial carcinoma and lung adenocarcinoma (Hussain et al., 2019). Structural rearrangements of ERC1 was previously reported in Chinese ESCC (Chang et al., 2017). ERC1 was also found as a prognostic biomarker in HNSCC (Szczepanski et al., 2013). Previously, a genome wide association study in Chinese Han population identified that SNP rs10052657 in PDE4D on 5q11 was associated with ESCC risk (Wu et al., 2011). Homozygous deletion of PDE4D was also identified in breast, lung and gastric cancers which established it as a tumor-promoting gene (Lin et al., 2013a; Lin et al., 2013b).

Breakage-Fusion-Bridge (BFB) is a mechanism supported by previous studies in cancer (Cheng et al., 2016; Yang et al., 2017). The BFB event is characterized by a special type of structural rearrangements called ‘fold-back’ inversion. Fold-back inversion can be defined as somatic structural variants with single inverted breakpoints exhibiting copy-number changes (Campbell et al., 2010). We implemented these information to identify BFB in each ESCC genome. In total, we detected 101 fold-back inversions across the 20 ESCCs, of which chromosome 11 appeared to have highest fold-back inversions (30) (Fig. S5A). BFB event was present in total 14 ESCC cases (70%) in this study. Moreover, fold-back inversions were observed on chromosome 11 around amplification of CCND1 locus (69455873-69469242) in eight ESCC cases (Fig. S5B). Notably, our CNA analysis identified amplification of CCND1 in fourteen patients, and 57% of those amplifications (8/14) were caused as a result of BFB. In addition, oncogenes such as FGFR1 (1/20) (Fig. S5C) and EGFR (1/20) were also detected in the amplified regions which were affected by BFB event. In all, this analysis presented an important insight of BFB in ESCC, and targeting the amplified oncogenes in therapies will benefit the ESCC patients in future.

Kataegis loci, which are localized hyper-mutation clusters, were identified in 30% of the ESCC patients (6/20) (Fig. S6A and Table S5). In total 11 kataegis loci were detected in six cases, of which four kataegis had SVs in their close vicinity. Previous studies reported kataegis in ESCC (Chang et al., 2017; Cheng et al., 2016), and in breast cancer it was observed in more than 50% of the cases (Nik-Zainal et al., 2012; D’Antonio et al., 2016). Furthermore, chromothripsis (Korbel & Campbell, 2013), a phenomenon that affects chromosomes with more than ten structural rearrangements with regular oscillation of copy number changes, was found in one ESCC patient (OK008) (Fig. S6B) on chromosome 12 and 14. Chromothripsis lead structural rearrangements have been reported in many cancers usually in low frequency, however, more than 40% chromothripsis were found in glioblastomas and lung adenocarcinomas (Korbel & Campbell, 2013; Cortes-Ciriano et al., 2020).

Somatic copy number alterations in ESCC

We called the somatic CNAs from the whole genome of the 20 pairs of matched tumor and normal samples to investigate CNAs in ESCC. GISTIC2.0 was then used to identify the recurrently amplified and deleted regions (Mermel et al., 2011). We identified four frequent amplified regions (3q26.33, 8p11.23, 11q13.3 and 14q21.1) and six deleted regions (1q21.1, 4q35.2, 5q13.2, 9p21.3, 10p12.33 and 21p11.2) (Figs. 4A and 4B, Tables S6 and S7).

The copy number amplification of 3q26.33 was found in 63% of the samples and included important caner driver genes PIK3CA, SOX2, FGF12 and TP63. Notably, SVs in TP63 was also found in 15% of the ESCC samples in this study. It was observed that 11q13.3 was the most frequently amplified region (74%) in our dataset, consistently with the observation in Chinese ESCC (Ying et al., 2012). 11q13.3 gain harbored many important genes such as CCND1, FGF3, FGF4 and FGF19 which established this region as a prominent target in ESCC. Amplification of 8p11.23 involved FGFR1, which plays an active role in cell growth and differentiation. Earlier studies reported FGFR1 as a potential drug target in many human cancers (Von Loga et al., 2015; Chang et al., 2014; Lin et al., 2014).

Copy number deletion of 9p21.3 was found in 79% of the samples and it was the most common deletion in the ESCC. This region contains CDKN2A, an essential regulator of cell cycle (Fig. 4C). 10p12.33 was deleted in 79% of the samples and harbored gene MRC1, a M2 macrophage antigen known to be associated with tumor development, invasion, metastasis and angiogenesis (Weber et al., 2016; Fang et al., 2017). Deletion of 4q35.2 was found in 63% of the samples. This region includes the tumor suppressor FAT1, which was recurrently deleted in colorectal cancers, glioblastoma and HNSCC (Morris et al., 2013). Notably, recurrent SNVs/INDELs in FAT1 was observed in our study as well. Overall, this study detected many copy number altered peaks and important ESCC-associated driver genes such as FGFR1, PIK3CA, CCND1, CDKN2A and MRC1, which have the potential to be used as therapeutic targets in future (Du et al., 2017; Lin et al., 2014; Padhi et al., 2017; Weber et al., 2016; Weber et al., 2014).

Assessment of druggable genes in ESCC

In order to examine the druggability of the genes identified by the CNA analysis, we analyzed the genes with the drug-gene interaction database (Griffith et al., 2013; Cotto et al., 2018). We found that 67 genes out of the 442 genes had at least one drug target that interact with it (Table S8). HTR1A, PIK3CA, FGFR1, ADRB3, HPIK3R1, MTNR1A, CDKN2A, CCND1, CDK7, and ANO1 were the top druggable genes that showed large interactions with drugs. Notably, amplification of CCND1 (74%), PIK3CA (63%), FGFR1 (37%) and deletion of CDKN2A (79%) and CDK7 (74%) genes were observed in a considerable amount of samples in our analysis. In all, CNA analysis was able to identify potentially druggable genes including CCND1, FGFR1, PIK3CA and CDKN2A which might be used for the treatment of ESCC in future.

The druggability of genes detected by our SNV analysis was also examined. We identified 212 genes out of the 1,111 genes that have at least one interaction with drug (Table S9). We identified PIK3CA, TP53, BRAF, NOTCH1, FGFR2, F11, MTOR, LRP2, and RB1 were the most prominent druggable genes with large number of drug-interactions. On the other hand, 149 genes out of 958 have the druggable property which were identified by SV analysis here (Table S10). LRP1B, PDE4D, WWOX, GPHN, KCNB2, FHIT, CDKN2A, TP53 and BRCA2 were the important druggable genes determined with large number of interactions with drugs. LRP1B (30%), PDE4D (25%), WWOX (25%), GPHN (20%) and KCNB2 (20%) were frequently affected by SVs in ESCC.

We also combined the analysis of druggable genes identified by all the mutational events SNVs, SVs and CNAs in ESCC. It was observed that some of the genes such as TP53, NOTCH1, RB1, JAK2, ACAN and SCN9A were commonly altered by both SNVs and SVs in ESCC (Fig. 5 and Table S11). Genes, for example, PIK3CA, F11, EPHB3 and TLR3 were commonly altered by SNVs and CNAs. We also detected some genes such as CDKN2A, ANO1 and NDUFB5 were altered by both SVs and CNAs. However, we found no druggable gene that was altered commonly by all the three mutation patterns in our analysis.