The genome landscape of ESCC and its stage-associated variations

We performed genome sequencing on DNA from 104 ESCC tumors and matched adjacent normal tissues, including 51 stage I and 53 stage III cases (Additional file 1: Table S1) [18]. Five stage I and nine stage III tumors underwent WGS (median coverage of 65×); 46 stage I and 44 stage III tumors underwent WES (median coverage of 132×; Additional file 2: Figure S1). The average number was 3.9 coding mutations/Mb for WGS and 2.4 non-silent mutations/Mb for WES (Additional file 3: Table S2) [18]. Candidate non-silent mutations identified from 48 stage I and 48 stage III tumors were selected for TCS (at least 365×; Additional file 3: Table S2C) [18]. The validation rates were 97.8 % for identified single-nucleotide variations (SNVs) and 58 % for indels. Consequently, we obtained 10,330 somatic point mutations in total: of these, 65 % resulted in missense changes and 6 % resulted in nonsense changes. There were 184 alterations of splice sites and 90 small indels: of the indels, 84 % introduced frameshifts and 16 % were in-frame (Additional file 4: Table S3 and Additional file 5: Figure S2B). We found 73 (range 24–189) non-silent mutations per tumor in this cohort, and this rate is in line with published rates (Additional file 5: Figure S2C) [15–17, 20], underscoring the representative nature of our analysis.

Next, to identify possible differences in the overall genomic architecture between stage I and stage III tumors, we compared the spectrum of mutations and broad SCNAs in these two sets. We observed no clear difference between stage I and stage III tumors in terms of genome-wide somatic SNV counts, even when the groups were subdivided into stage IA, IB, IIIA, IIIB and IIIC (Additional file 6: Figure S3). However, we observed marked differences in the stage-associated pattern of SCNAs. Stage I tumors harbored fewer SCNAs than stage III tumors (p < 2.2 × 10⁻¹⁶, Fig. 1a) but there was no correlation within sub-stage IIIA, IIIB or IIIC (Additional file 6: Figure S3B). Genomic identification of significant targets in cancer (GISTIC) [21] analysis in the WGS set yielded universal deletions affecting 4p, 11p, 16p, 19p and 19q, and frequent gains of 3q, 5p, 7p, 7q, 8p, 8q, 12p, 14q, 18p, 20q, 21q, Xp and Xq (Additional file 6: Figure S3C). In particular, recurrent somatic amplifications at 8q (containing MYC and FAM84B) [22] were found to be enriched in stage I tumors; the deletions of 4p-q (containing VEGFC, FBXW7 and FAT1) [23, 24] and 5q (containing PTTG and MAML1) [25, 26] were particularly identified in stage III tumors (p < 0.05, Fisher’s exact test, Fig. 1b), suggesting that these alterations are associated with stage progression. Furthermore, copy-number analyses verified the amplifications of candidate genes located within these significantly altered regions in 36 atypical hyperplasia tissues (Fig. 1c). Thus, although stage I and stage III tumors of ESCC are genomically similar, our results reveal that the copy-number variations exhibit a pattern that is associated with the clinical stage of the tumor.

We also applied a modified GISTIC method to profile genome segments with copy number variations in the 14 tumors analyzed by WGS, which revealed 126 significantly altered regions (Additional file 7: Figure S4). Moreover, to identify genes affected by recurrent SCNAs, we manually inspected the 126 significantly altered regions using the Integrative Genomics Viewer (IGV). This approach identified recurrent focal CNAs, including one of the most amplified regions, 8q24.13-q24.21, which contains FAM84B (Additional file 8: Table S4) [27]. Amplification of this gene was found in 44 % (46 out of 104) of patients (Fig. 2a and Additional file 9: Figure S5) and this was further validated by fluorescence in situ hybridization (FISH) in atypical hyperplasia tissues and ESCC tumors (Fig. 2b). FAM84B is involved in the formation of DNA-repair complexes and little is known about its function in human cancers [27]. In our cohort of 104 patients, this gene was markedly highly expressed in 57 % of cases (the immunoreactivity score for FAM84B expression in tumors was at least double that of matched normal tissue, T_IRS/N_IRS > 2, Fig. 2c). To verify the specific abundance of this gene during the initiation of ESCC, we used immunohistochemistry on tissue microarray-sourced preclinical and ESCC samples. We observed a marked relative increase in FAM84B expression in atypical hyperplasia samples compared with that of normal tissues (Fig. 2d). In addition, we knocked down FAM84B in KYSE150 and TE-1 cells that have high levels of endogenous FAM84B and observed that FAM84B depletion attenuated cell growth, migration and invasion (Fig. 3). This finding suggests that FAM84B amplification and the resultant increased levels of FAM84B protein are associated with progression from normal squamous epithelium to ESCC and may be a potential diagnostic marker for susceptibility to ESCC.

Significantly mutated genes

For most cancer types, genomic landscapes consist of many ‘hills’ (corresponding to genes that are altered in a low percentage of the tumors) and a small number of ‘mountains’ (genes that are altered in a higher percentage of tumors) [28]. This also holds true for ESCC, as revealed by our cohort (Additional file 10: Figure S6; Additional file 11: Table S5) [18]. Only 14 genes (TP53, TTN, MUC16, NOTCH1, FAT1, PIK3CA, CSMD3, PCLO, LRPIB, MLL2, EP300, SYNE1, SPTA1 and PKHD1L1) contained somatic mutations that were detected in at least 10 % of the samples whereas 85 genes were mutated in 5-10 % of the tumors. The remaining 1,309 genes were altered in 2–5 % of the cohort of 104.

We applied MutSigCV to identify significantly mutated genes (SMGs) associated with ESCC. This analysis led to the identification of eight SMGs with q < 0.2 (Additional file 12: Table S6A). In general, the number of SMGs in an individual cancer ranges from zero to four. Of the 104 tumors that we comprehensively characterized, 50 % (52 out of 104 cases) displayed mutations in two or more SMGs, 44 % (46 out of 104 cases) of cases harbored alterations in one SMG, and six cases did not have mutations in any of the SMGs. We observed 21 different combinations of SMGs and so there seems to be substantial variation in the drivers of oncogenesis. Four of the eight SMGs (TP53, NOTCH1, PIK3CA and FAT1) were mutated in more than 10 % of cases. Collectively, these contributed 86 % (158 of 184) of the non-silent mutations. The other four SMGs, each contributing relatively infrequently, were responsible for the remaining 14 % (26 of 184) of the non-silent mutations (Fig. 4).

Identical mutations suggest stage-associated NOTCH pathway alterations

To discern the SMGs and pathways relevant to tumorigenic capacity and tumor progression in more detail, we investigated the relationship between the prevalence of each SMG and tumor stage and conducted separate pathway analyses for stage I versus stage III tumors. Interestingly, a pronounced diversity of NOTCH1 mutations was observed in 35 % of stage I tumors but in only 8 % of stage III tumors (p < 0.0006, Fisher’s test). Likewise, NOTCH signaling, a fundamental signaling system comprising the NOTCH receptor (including NOTCH1, NOTCH2, NOTCH3 and NOTCH4), Delta and Serrate/Jagged (DSL) ligands and CSL DNA-binding proteins [28], was altered in 55 % of stage I tumors versus 32 % of stage III tumors (p < 0.02, Fisher’s test, Fig. 5a). This indicates that alteration of the NOTCH pathway may be early event in the development of a subgroup of ESCCs. NOTCH1 mutations are relatively common in head and neck squamous cell carcinoma (HNSCC), lung SCC and breast cancer, with 5 % to 15 % of tumors harboring protein-coding changes [28–30]. In our cohort, 22 somatic mutations were identified in NOTCH1 with a mutation frequency of 21 % (Additional file 12: Table S6B). Of the 22 non-silent mutations, eight are nonsense, ten are missense, three are frameshift indels and one is a splicing site mutation, and these are, in general, all predicted to be loss-of-function mutations (Fig. 5b). Notably, we found another six non-silent mutations including four missense and two frameshift indels, with a mutation frequency of 16.7 % in the 36 atypical hyperplasia tissues (Additional file 13: Table S7). Many of these missense mutations occurred at or near important domains such as the ligand-binding domain (EGF repeats). Moreover, the nonsense mutations observed in NOTCH1 generate a premature stop codon, resulting in a C-terminally truncated NOTCH1 protein lacking a PEST sequence (a sequence rich in proline, glutamic acid, serine and threonine) that is important for transcription activation. In addition, five out of nine mutations identified in NOTCH2/3 were truncating, and two stopgains were identified in RBPJ, one of the target genes of NOTCH. Thus, in contrast to T-cell acute lymphoblastic leukemia, chronic lymphoblastic leukemia and breast cancer, in which NOTCH1 serves as an oncogene [31], this pattern of mutations suggest that the NOTCH pathway has a tumor suppressing role in ESCC.

Of related interest, 8 % of ESCC tumors (12 % of stage I tumors and 4 % of stage III tumors) harbored mutations in the F-box protein FBXW7 (Fig. 5a). FBXW7 is the substrate-recognition subunit of an SCF-type ubiquitin ligase complex that regulates the cell cycle by targeting many proto-oncoproteins and antiapoptotic molecules including NOTCH1, cyclin E, c-Myc, c-Jun and mTOR for ubiquitin-mediated degradation. It therefore acts as a tumor suppressor protein [32] and has been found to be mutated in various tumors [33, 34]. Of the eight mutations detected in our cohort, two nonsense and one frameshift deletion are inactivating mutations (Fig. 5b, bottom panel). Moreover, most of the mutations occur within the WD40 domain involved in substrate recognition; mutation of this site prevents the recognition of targets such as NOTCH1 for degradation. Taken together, we identified 31 mutations in NOTCH family genes. Most of these mutations result in truncated protein products or have deleterious effects, suggesting that NOTCH signaling is significantly disrupted in our sample set and may be one of the main mechanisms associated with the development of a subgroup of ESCCs.

In addition to the NOTCH pathway, we found that mutations in genes involved in several major metabolic pathways were significantly enriched in the stage I group. In particular, mutations occurred in genes associated with pyrimidine metabolism (31 % of stage I tumors versus 13 % of stage III), glycine/serine/threonine metabolism (16 % versus 2 %) and fructose and mannose metabolism (16 % versus 2 %). Conversely, the hedgehog (Hh) signaling pathway showed a significantly higher mutation frequency in the late stage III group (p < 0.05, Fisher’s test; Additional file 14: Table S8). Thus, despite there being no marked differences between the tumors at various stages at the level of individual genes, such differences seem to exist at the level of pathways. For other frequently or significantly mutated genes, we found no significant correlation between their mutation frequencies and tumor stage.

Samples and clinical data

We recruited tumor samples and adjacent normal tissues from 104 ESCC patients from the Han Chinese population who live in the Shanxi and Henan provinces, Taihang Mountain, North-Central China. All samples were obtained before treatment according to the guidelines of the local ethical committees (IRB of Shanxi Medical University, Approval No. 2009029, and the Ethics Committee of Henan Cancer Hospital, Approval No. 2009xjs12). This study was approved by the ethical committee of the Shanxi and Henan, China. All ESCC cases collected for this study were staged using the American Joint Committee on Cancer (AJCC) Cancer Staging Standards, 7th Edition (2010). This cohort includes 51 cases of stage I tumors and 53 cases of stage III tumors. Different subsets of patients were assayed on each platform: 14 tumors and matched normal samples from 14 patients, including 5 cases of stage I and 9 cases of stage III, had WGS data available (65×); 90 samples from 46 cases of stage I and 44 cases of stage III, had WES data available (132×); 96 of the 104 samples including 48 cases of stage I and 48 cases of stage III had TCS data available (365×) [18]. A detailed description of the clinical characteristics of this cohort is shown in Additional file 1: Table S1. A summary of the next-generation sequencing analyses in this study is presented in Additional file 2: Figure S1.

Sequencing

For WGS, genomic DNA extracted from 5 stage I tumors, 9 stage III tumors and matched normal tissues were randomly fragmented and purified. The WGS libraries were constructed and subjected to WGS on an Illumina HiSeq 2000. At least 65× target depth and 30-fold haploid coverage for tumors and normal samples were achieved in all samples. A detailed description is presented in Additional file 3: Table S2.

For WES, the qualified genomic DNAs from 46 tumors from stage I, 44 tumors from stage III and matched normal tissues were randomly fragmented, amplified by ligation-mediated PCR (LM-PCR), purified and hybridized to the NimbleGen SeqCap EZ exome (44 M, Roche Company, Indianapolis, IN, USA) array for enrichment. Each captured library was subjected to an Illumina HiSeq 2000 platform for high-throughput sequencing. The mean coverage achieved was 130× in the tumor and 133× in the normal tissues. A detailed description is presented in Additional file 3: Table S2B. The Agilent SureSelect in Solution is described elsewhere [44].

Illumina sequencing analysis pipeline

For detection of somatic point mutations, sequencing reads from an Illumina HiSeq 2000 sequencer were aligned to the human reference genome (hg19) sequence using Burrows-Wheeler Aligner (BWA) [45]. After removing duplicated reads (redundant information produced by PCR) using the SAMtools, an in-house cancer sequencing analysis pipeline (CSAP) was used to identify somatic mutations [15, 46]. Several recently published studies have compared the somatic mutation calling tools (e.g. SomaticSniper, Virmid, Strelka, MuTect and VarScan2) and found that VarScan2 outperformed all other tools in detecting high-quality somatic mutations (high coverage and allele frequency) [47–49]. Therefore, we detected mutations by VarScan2.2.5 with previously determined parameters [50]. High-confident SNVs were annotated with ANNOVAR and used in follow-up analyses [51]. Briefly, we used an in-house pipeline for statistic sequencing error rates, base calling accuracy and ATCG base contents of raw sequencing data. Sequence data per lane needed to meet the following quality controls: i.sequencing error rates ≤ 2.5 %; ii. base calling accuracy, measured by the Phred quality score (Q score), Q20 ≥ 80 %, Q30 ≥ 75 %; iii. GC or AT separation rate ≤ 0.4 %. We applied SOAPnuke to remove adapter and filter low quality reads [52]. Reads were kept if they met the following criteria: i. adaptor rate ≤ 10 %; ii. N rate of every single read ≤ 10 %; iii. low quality base rate (base quality < 5) of every single read ≤ 50 %. After removing reads containing sequencing adaptors and low-quality reads, the high-quality single-end reads were aligned to the NCBI human reference genome (hg19) using BWA by default parameters. Each sample should meet the following criteria: i. average sequencing depth of WGS ≥ 30×, average sequencing depth of WES ≥ 100×; ii. Mapping rate ≥ 95 %; iii. Mismatch rate ≤ 10 %. High-confident somatic SNVs were called when they met the following criteria: i. both the tumor and normal samples should be covered sufficiently (≥10×) at the genomic position; ii. the variants should be supported by at least 10 % of the total reads in the tumor compared with less than 2 % in normal tissues; iii. the variants should be supported by at least three reads in the tumor; iv. distance between adjacent somatic SNV distance should be over 10 bp; v. mapping qualities of reads supporting mutant alleles in the tumor should be significantly higher than 30 (Wilcoxon rank sum test, p < 0.2); vi.base qualities of reads supporting mutant alleles in the tumor should be significantly higher than 20 (Wilcoxon rank sum test, p < 0.05); vii.mutations should not be enriched within 5 bp 5′ or 3′ of read end (Wilcoxon rank sum test, p < 0.1); viii. the changes of mutant allele frequency between tumor and normal should be statistically significant (Fisher’s exact test, p < 0.05); ix. reads supporting mutations should not be significantly enriched within either forward or reverse genomic strand (Fisher’s exact test, p < 0.0001).

Detection of small indels

The indel calling step was performed by a GATK SomaticIndelDetector with default parameters [53]. The high-confident indels were identified by an in-house pipeline and further annotated with ANNOVAR [51] as either germline or somatic based on whether any evidence for the event at the same locus was observed in the normal data. High-confident somatic insertions and deletions (indels) were called through the following steps: i. candidate somatic indels were predicted with the GATK SomaticIndel Detector with default parameters; ii.for each predicted somatic indel, local realignment was performed with combined normal and tumor BAM files; iii. high-confident somatic indels were defined after filtering germline events.

Target capture sequencing

To provide high-confident mutations, non-synonymous mutations identified in the WGS and WES sets were selected for TCS. Briefly, non-synonymous SNVs (6873) and indels (125) in coding regions identified from 96 of the 104 samples were designed on a Nimblegen customized capture array (Roche Company, Indianapolis, IN, USA). Genomic DNAs from 48 stage I tumors, 48 stage III tumors and matched normal tissues were fragmented and libraries were constructed following the same method as the exome-capture experiment. For SNVs and small indels, we tiled ~200 bp targets across the variant of interest, including a minimum buffer of 100 bp in each direction. After library preparation and hybridization, sequencing was performed on an Illumina HiSeq 2000 platform. Each tumor was sequenced to at least a depth of 300× and SNVs were called with the same pipeline except that variant allele frequency was decreased to 5 % to guarantee that low-frequency mutations were retained. Somatic indels were manually inspected across 104 normal samples to remove germline mutations. The mean coverage achieved was 365× in TCS with 354× in tumors and 375× in the normal samples. A detailed description is presented in Additional file 3: Table S2C. Validation lanes were aligned to the reference sequence and BAM files created in the same manner as described above.

DNA copy number analysis

We performed SegSeq [54], a widely used method to identify copy number variation (CNVs) by comparing a tumor sample with a matched normal sample, to infer somatic CNVs in ESCC genomes based on WGS reads. The resulting copy number segments were mapped to individual genes to determine gene-level copy numbers and copy gain/loss statuses. Copy numbers of ≤ 1.5 were considered to indicate deletions and ≥ 2.5 were considered to be amplifications. To infer recurrently amplified or deleted genomic regions, we re-implemented the GISTIC [21] algorithm using copy numbers in 1 kb windows as markers instead of SNP array probes. G-scores were calculated for genomic and gene-coding regions based on the frequency and amplitude of amplification or deletion of each gene. A significant CNV region was defined as having an amplification or deletion with G-score > 0.1, corresponding to a p-value threshold of 0.05 from permutation-derived null distribution.

Identification of significantly mutated genes

For the identification of significantly mutated genes, we applied a novel analytical methodology, mutation significance with covariates (MutSigCV) [10], to avoid the false-positive findings detected by the standard significance analysis method (MutSig1.0) [55]. MutSigCV corrects for variation by using patient-specific mutation frequency and spectrum and gene-specific background mutation rates, incorporating expression levels and replication time. MutSigCV is freely available for non-commercial use [56].

Pathway enrichment analysis

We performed the pathway enrichment analysis using the Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.7 by examining the distribution of the non-synonymously mutated genes identified within the KEGG database [57] as described previously [18]. Significantly altered pathways were determined by p-values calculated based on hypergeometric distribution with Benjamini correction.

Cell lines

KYSE150 and TE1 ESCC cell lines that were tested and found to be free of mycoplasma contamination were used in this study. 293 T cells were used as a packaging cell line to produce virus. All cells were grown in DMEM/F12 media at 37 °C in 5 % CO₂. For the functional analysis of FAM84B, the ESCC lines KYSE150 and TE-1 with high endogenous expression levels were used for knockdown experiments.

Knockdown of FAM84B

Knockdown experiments were performed in at least two ESCC lines with high endogenous FAM84B expression. Two independent shRNAs were cloned into the pLKO.1-puro vector (Addgene, Cambridge, MA, USA) as described previously [58]. A non-specific targeting shRNA was also cloned into the pLKO.1-puro vector to be used as a scrambled control (SCR). shRNA knockdown efficiency was determined by western blot analysis for FAM84B proteins using an anti-FAM84B antibody (Proteintech, Chicago, USA). Relative expression was normalized to the β-actin expression level.

qPCR copy number validated analysis

The copy numbers of the genes of interest were assessed in frozen atypical hyperplasia samples using genomic qPCR (TaqMan, Applied Biosystems, Foster, CA, USA) in triplicate. Prevalidated primers for the relevant genes were obtained from Applied Biosystems (accession numbers Hs02758348_Cn for MYC, Hs00655850_Cn for FAM84B, Hs01654625_Cn for VEGFC, Hs05965356_Cn for FBXW7 and Hs00703603_Cn for FAT1). RNase P (RPPH1 gene; Life Technologies, Shanghai, China, 4403328) was used as a diploid control. Data were analyzed using the comparative (delta-Ct) Ct method. An inferred copy number of < 0.3 was considered to indicate a homozygous deletion.

Fluorescence in situ hybridization analysis

To evaluate the amplification of FAM84B, we performed FISH. Tumor and matched normal tissues of interesting ESCC cases were cut into pieces in PBS, swollen in 65 mmol/L KCl for 5 mins at 37 °C, followed by fixation in cold acetic acid/methanol for 5 mins at 4 °C, then dropped onto slides. For interphase FISH analysis, slides were stained with Cytocell enumeration probes against chromosome 2 or FAM84B conjugated with FITC or Cy3.5 (Rainbow Scientific, Windsor, CT, YSA). Probes against chromosome 2 were used as controls. Staining was carried out according to the manufacturer’s protocol. FISH samples were viewed with a fully automated, upright Zeiss Axio-ImagerZ.1 microscope with a × 20 objective and DAPI, FITC and Rhodamine filter cubes. Images were produced using the AxioCamMRm CCD camera and Axiovision version 4.5 software suite. p-values were calculated using a two-sample test for equality of proportions with continuity correction.

Immunohistochemistry and tissue microarray resource

FAM84B was immunohistochemically stained with an anti-FAM84B antibody (Proteintech, Chicago, USA) as described previously [18]. A human ESCC tissue array (OD-CT-DgEso01-001) purchased from Shanghai Outdo Biotech (Shanghai, China) was used to detect the expression level of FAM84B. OD-CT-DgEso01-001 array contains 10 normal esophagus tissues, 22 atypical hyperplasia tissues and 58 ESCC cases (2cores/case; Additional file 16: Table S9). FAM84B immunohistochemistries were performed on tissue microarraysas previously described [59], using the polyclonal anti-FAM84B antibody (Additional file 17: Figure S8). Briefly, sections were incubated with the specific antibody at a 1:500 dilution for 14 h at 4 °C, followed by detection using the PV8000 (Zhongshan, Beijing, China) and DAB detection kit (Maixin, Fuzhou, China), producing a dark brown precipitate. Slides were counterstained with hematoxylin. All images were captured at × 100. Cytoplasm expression of FAM84B was quantified using AperioCytoplasma 2.0 software. Statistic analyses were performed using Graphpad Prism 5.0.

Immunoblotting

Immunoblotting was performed as previously described [19] using anti-FAM84B (Proteintech, Chicago, USA). Antibody binding was detected using horseradish peroxidase-labeled anti-mouse (Sigma, Santa clara,USA) or anti-rabbit (Cell Signaling, Boston,USA) antibodies and chemiluminescence was measured using a LAS4000 device chemiluminescence System (Sagecreation, Beijing, China). Equal protein loading was confirmed with antibodies against β-actin (Transgen, Beijing, China).

MTT and migration/invasion assays

To assess cell viability, MTT assays were performed as previously described [59]. Each experiment consisted of four replications and at least three independent experiments were carried out. Migration and invasion assays were performed in 16-well CIM plates in an xCELLigence RTCA DP system (ACEA Biosciences, San Diego, USA) using BD matrigel basement membrane matrix for real-time cell migration analysis as described previously [60, 61]. At least three independent experiments were carried out; for each independent experiment, five duplicates were performed for each group.

Statistical analysis

Experiments were done in triplicate and data were presented as mean ± SD. Student’s t-test was used for statistical analysis, and data from more than two groups were analyzed by one-way analysis of variance (ANOVA) in SPSS Statistics 19.0 followed by LSD-t test. Results were considered significant when p < 0.05.