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).

Open in a separate window

Fig. 4

A view of the genome landscape of ESCCs at various stages. a Stage-distributions of candidate driver mutations identified by MutSigCV significance analysis. The type of each mutation is shown for every sample, including the gene-specific total number of mutated samples (right); mutation subtypes are denoted by color. If multiple mutations were found in a gene in a single sample, only one is shown. The significance of the mutations in each gene is shown to the left by the false discovery rate (FDR) q value. The lower bars indicate smoking/drinking status and family history, respectively. The full list of mutated genes is given in Additional file 8: Table S4. b Percentage of ESCC patients harboring one or two SMGs in the 104 cohort. The ESCC patients harboring none of the identified SMGs are also shown

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.

We thank BGI for sequencing and analyzing data. This work was supported by the National Natural Science Foundation of China (81330063 & 81272189), the Key Project of Chinese Ministry of Education (NO213005A), Specialized Research Fund for the Doctoral Program of Higher Education (20121417110001), Research Project Supported by Shanxi Scholarship Council of China (2013–053 & 2015 key 3), the Innovative Team in Science & Technology of Shanxi (2013-23), the Program for the Outstanding Innovative Teams of Higher Learning Institutions of Shanxi(2015-313) to Y.P.C., and the National Natural Science Foundation of China (81272694 to X.L.C., 81201956 to J.L., 81402342 to L.Z.). The authors declare no competing financial interests.