Tumor tissues and blood were collected from treatment-naive Chinese patients with oral cancer, who then completed the treatment and the follow-up for at least 12 months by their physician. Tumor tissues were obtained by precise laser micro-surgical resection, and blood samples were drawn through the upper arm vein. Ethics approval was obtained from the human Ethics Committee of Ninth People's Hospital (app. no. 2014011). Written informed consent was obtained from all patients. All patients were diagnosed with oral cancer without having any other type of cancer. Histological sections of lesions were examined by two oncopathologists, and the diagnosis was based on the World Health Organization (WHO) guidelines. Clinicopathological staging was determined according to the Union for International Cancer Control (UICC) TNM staging system v6. The patient characteristics are summarized in Table I. All patients received surgery combined with radiation therapy and/or chemotherapy. Each patient was followed-up each month after surgery. The clinical detection of recurrence was confirmed by histopathology. All patients experienced residual or recurrent disease.

The bisected tissue samples and blood were stored in liquid nitrogen. Thereafter, genomic DNA was isolated according to standard protocols as previously described (9–11). The genomic DNA libraries were constructed using protocols provided by Illumina (San Diego, CA, USA). Approximately 50 Mb of sequences from the whole exons and their flanking regions were enriched from the fragmented genomic DNA using the Agilent SureSelect DNA Capture array. The captured DNA libraries were processed using the Illumina HiSeq 2500 Genome Analyzer with an overall >10X coverage depth.

Sequence data in FASTQ files were quality-checked using FASTQC v0.11.4 (http://www.bioin-formatics.babraham.ac.uk/projects/fastqc/). The filtered sequence reads were mapped to the hg19 reference sequence using Burrows-Wheeler Aligner (BWA) v0.7.12 (http://source-forge.net/projects/bio-bwa/) and SAM files were generated. Local realignment around indels, duplicate removal and Base quality score recalibration were performed using the Genome Analysis Toolkit (GATK) v3.5 (https://www.broadinstitute.org/gatk/). Alignment SAM files were converted to BAM files using Samtools v1.3. Somatic variants were called using MuTect v1.1.7 (https://www.broadinstitute.org/gatk/download/auth?package=MuTect) with default parameters. High confidence variants were identified with the criteria of a minimum of 10 reads covering a site in the tumor tissue and 8 reads in the blood, and the maximum alternative allele frequency in the blood <0.3 of that in the tumor tissue. Somatic mutations which are defined as called mutations in the absence of corresponding reads in the blood DNA samples were annotated using the Annovar main package (http://annovar.openbioinformatics.org/en/latest/) to infer the locations of the variants within the genes. Somatic single nucleotide variants (SNVs) were then combined and searched in the dbSNP and COSMIC databases. The reference human genome hg19 and dbSNP132 were used to call and annotate the SNVs. The novel variant identified by WES was verified by Sanger sequencing. The primers used for PCR and sequencing were 5′-AAGTTGTCGTAGAGGCAGGC-3′ as the forward primer and 5′-GACCCCTGGCGGCAATG-3′ as the reverse primer. PCR was used for Sanger sequencing validation of the new variant. The primers used for PCR and sequencing were forward, 5′-AAGTTGTCGTAGAGGCAGGC-3′ and reverse, 5′-GACCCCTGGCGGCAATG-3′. PCR was performed on an ABI9700 PCR instrument (Applied Biosystems, Inc., Foster City, CA, USA). in a 30 µl reaction volume including 19 HotStarTaq buffer, 2.8 mM Mg²⁺, 0.1 U of HotStarTaq polymerase (Qiagen, Inc., Valencia, CA, USA), 2 ng of blood genomic DNA, 0.5 pmol of each primer, 0.5 mmol of dNTPs and 0.5 mmol of ddNTPs. Thermocycling was carried out at 94°C for 15 min, followed by 45 cycles at 94°C for 20 sec, 56°C for 30 sec, and 72°C for 1 min, with a final incubation at 72°C for 3 min.

First of all, the effects of the SNVs were evaluated using Mutation Taster (http://www.mutationtaster.org/) (12), which automatically yielded types of predictions, i.e., 'disease_causing_automatic', 'disease_causing', 'polymorphism' and 'polymorphism_automatic', and the probability P-value of the prediction. A P-value close to 1 indicates a high possibility of the prediction results. A transcript ID and a sequence covering each single-nucleotide polymorphism (SNP) were required as input.

Sorting Intolerant From Tolerant (SIFT; http://sift.jcvi.org/) uses a query sequence and multiple alignment information to predict roles of variants (13). The alleles with calculated probabilities <0.05 are predicted to be deleterious; probabilities ≥0.05 are tolerated (14). The protein FASTA AA sequence was used as the query sequence and novel mutated SNPs were analyzed.

Polymorphism Phenotyping v2.0 (PolyPhen-2; http://genetics.bwh.harvard.edu/pph2/) predicts the influence of variants through automatically selecting 8 sequence-based and 3 structure-based predictive features by an iterative greedy algorithm. If the classifier probability is <0.15, the allele is predicted to be benign, and if the probability is ≥0.15, the allele is predicted to be probably or possibly damaging. Protein FASTA AA sequence, amino acid position and amino acid variants are required (15).

PANTHER (http://www.pantherdb.org/tools/csnpScore-Form.jsp) predicts the functional consequences of substitutions on the protein by calculating the substitution position-specific evolutionary conservation (subPSEC) score derived from the Hidden Markov Model and alignment score (16). The SNPs with subPSEC ≤-3 were considered damaging.

SNPs&GO (http://snps.biofold.org/snps-and-go/) is based on SVM and functional information codified by gene ontology (GO) terms to predict the impact of protein variations. The present analysis was protein sequence-based. The prediction (either disease or neutral) and RI (probability of disease-related class) are given (17).

SNAP2.0 (https://rostlab.org/services/snap/) predicts the effect of mutations by identifying all variants based on neutral network method that utilizes protein information. It requires protein sequence as input and yields 80% accuracy (18).

I-Mutant v2.0 (http://folding.biofold.org/i-mutant/i-mutant2.0.html) is a neural network-based tool for the routine analysis of protein stability and alterations caused by the single-site mutations, which also provides the scores for the Gibbs free energy change (DDG) calculated with the FOLD-X (19). The FOLD-X analysis tool provides the comparison between wild-type and mutant models in the form of van der Waals clashes, which greatly influence the energy decomposition. The protein FASTA AA sequence was used as the query sequence.

The residues mutation was performed using the SwissPDB Viewer viewer and energy minimization for 3D structures was performed by the NOMAD-Ref server (http://lorentz.immstr.pasteur.fr/gromacs/minimization_submission.php)(20). This server uses Gromacs as the default force field for energy minimization based on the methods of steepest descent, conjugate gradient, and L-BFGS. We used the conjugate gradient method for optimizing the 3D structures. The deviation between the two structures is evaluated by their root-mean-square deviation (RMSD) values.

To examine the stability of the native and mutant modeled structures, the identification of the stabilizing residues is useful. We used the server SRide (http://sride.enzim.hu/) (21) for identifying the stabilizing residues in native protein and mutant models. Stabilizing residues were computed using parameters, such as long-range order, surrounding hydrophobicity, conservation score and stabilization center.

In total, tissue and blood samples from 5 patients. All 5 patients were male with a mean age of 44.8±2.28, ranging from 42 to 48 years (Table I). The patients with oral cancer were all HPV-negative and were all alcohol consumers, but were not tobacco users. All patients presented with advanced disease (stage IV)s.

Coding exons of approximately 20,000 protein coding genes were sequenced from DNA isolated from the blood and primary tumor of each patient. The mean depths of sequencing were 136.90±11.96 for blood DNA and 139.87±16.29 for tumor DNA. The details for each patient are provided in Table II. The coding regions of genomes of the 5 patients contained 253 somatic variants, of which 51 (20.2%) were predicted to be synonymous, 138 (54.5%) were missense and 5 (2%) nonsense. The average mutation rate/Mb was estimated to be 2.04±0.73 (Table III).

The TP53 gene was mutated in at least 4 samples and 10 genes in 2 or more samples (Fig. 1). Since NOTCH signaling has been suggested to play important roles in head and neck squamous cell carcinoma (2), a newly identified variant (chr19:15288426A>C) in the NOTCH3 gene, a homology of NOTCH1, and a known variant (chr6:32168732C>T) in the NOTCH4 gene, another homology of NOTCH1 attracted our attention. The variant (chr6:32168732C>T) had a frequency of 8.35E-06 in the database of ExAC (http://exac.broadin-stitute.org) and was a silent mutation, resulting in no amino acid substitution. We then confirmed the newly identified variant (chr19:15288426A>C) in the NOTCH3 gene using Sanger sequencing (Fig. 2). This nonsense SNP (nsSNP) has the same location with the known nsSNP of rs761506399, but has a different allele mutation. The nsSNP, rs761506399, is an A>G mutation and results in 1438 Asn changing into Ser, while the new variant is an A>C mutation and results in 1438 Asn changing into Thr.

We first evaluated the functional effect of this nsSNP by Mutation Taster, which was predicted to be 'disease causing' with a probability of 0.9997 and suggested to possibly affect protein features.

Subsequently, SIFT, which predicts tolerated and deleterious substitutions at each position of the query sequence with multiple alignment score (13), indicated that this nsSNP was tolerated (score, 0.24) (Table IV). PolyPhen-2, which predicts the functional importance of an allele replacement using a training set and naive Bayes classifier, indicated that this nsSNP was 'possibly damaging' with a score of 0.986. We also examined the probability of the variant causing a deleterious functional change by PANTHER, which is based on the Hidden Markov model and alignment score in a very low value of subPSEC, −2.0575. Another method utilized was SNAP2.0, which is a neutral network method combined with a training set. It designated that this nsSNP affects the function of the NOTCH3 protein. Finally, the mutants were analyzed using SNPs&GO, which is an SVM-based classifier and obtains quite accurate (79%) results. SNPs&GO predicted that this nsSNP was neutral (RI=5) (Table IV).

I-Mutant v2.0, which is a neural network based tool for the routine analysis of protein stability and alterations caused by the single-site mutations, indicated that it can decrease protein stability with a DDG value of −1.67 (Table IV).

Thus, the 3D structure of NOTCH3 or the domain of NOTCH3 was analyzed. The 3D structure of part of NOTCH3 was resolved as 4ZLP.pdb and the novel mutation was located within the sequence of 4ZLP.pdb. Subsequently, the mutation residues of NOTCH3 were anlayzed using the SwissPDB Viewer and energy minimization for 3D structures was performed using the NOMAD-Ref server (20) (Fig. 3). We used the conjugate gradient method for optimizing the 3D structures. The total energy for the native structure was −1,229.84 kcal/mol, and that for the mutant structure (N1438T) was −1,051.39 kcal/mol. The difference between the wild-type structure and the N1438T structure was evaluated by the RMSD value of 1.63.

The computed total solvent accessible surface area of the wild-type protein was 10,287 Ų and thesurface area excluding the solvent was 8,285 Ų, while the total solvent accessible surface area of theN1438T protein was 10,268 Ų, and the surface area excluding the solvent was 8,268 Ų. The mutant amino acid reduces the solvent accessible surface area by 19 Ų for N1438T, and solvent excluded surface area by 17 Ų in spite of having more number of atoms as compared to the wild-type. The reduction in surface area may have little effect on the solubility of the protein, and thus on its function. The mutant amino acid is also an AA without charge just like the wide type AA. The hydrophobicity slightly increases for domain with mutant AA (data not shown).

The SRide server (21) was used to identify the stabilizing residues of the native structure and mutant modeled structures with the conservation score threshold of 6, LRO threshold of 0.012, and surrounding hydrophobicity threshold of 20. Four stabilizing residues were identified in the native structure, which were CYS67, ASP70, GLY110 and PRO173. The single-mutant model (N1438T) has 3 identified stabilizing residues, which were CYS67, ASP70 and GLY110. This analysis revealed that one stabilizing residue in the N1438T-mutant model less than the native protein structure may decrease the stability of the N1438T mutant protein.