Uncommon mutational profiles of metastatic colorectal cancer detected during routine genotyping using next generation sequencing

RAS genotyping is mandatory to predict anti-EGFR monoclonal antibodies (mAbs) therapy resistance and BRAF genotyping is a relevant prognosis marker in patients with metastatic colorectal cancer. Although the role of hotspot mutations is well defined, the impact of uncommon mutations is still unknown. In this study, we aimed to discuss the potential utility of detecting uncommon RAS and BRAF mutation profiles with next-generation sequencing. A total of 779 FFPE samples from patients with metastatic colorectal cancer with valid NGS results were screened and 22 uncommon mutational profiles of KRAS, NRAS and BRAF genes were selected. In silico prediction of mutation impact was then assessed by 2 predictive scores and a structural protein modelling. Three samples carry a single KRAS non-hotspot mutation, one a single NRAS non-hotspot mutation, four a single BRAF non-hotspot mutation and fourteen carry several mutations. This in silico study shows that some non-hotspot RAS mutations seem to behave like hotspot mutations and warrant further examination to assess whether they should confer a resistance to anti-EGFR mAbs therapy for patients bearing these non-hotspot RAS mutations. For BRAF gene, non-V600E mutations may characterise a novel subtype of mCRC with better prognosis, potentially implying a modification of therapeutic strategy.


Results
We retrospectively collected data from 857 mCRC samples including 779 samples with valid NGS results. DNA quality was suitable for 91% of the samples for NGS and uncommon mutational profiles were reported in 22 (2.7% of total) samples.
The histological subtypes were: adenocarcinoma for 11 samples, liberkunhian adenocarcinoma for 9 samples, mucinous carcinoma for 1 sample, and ductal carcinoma for 1 sample. Twenty samples were from primary tumors and 2 were from liver metastases ( Table 1).
The range of coverage for rare KRAS, NRAS or BRAF mutations was 700 to 27 000x. Mutant allele fraction was between 1.8% and 40.6% for non-hotspot NRAS mutations, between 4.0% and 14.5% for non-hotspot KRAS mutations and between 7.9% and 46.2% for non-hotspot BRAF mutations. Thirty-five mutations were missense mutations, 5 were silent mutations and one was a stop mutation.
For each mutation, the PolyPhen-2 score, SIFT score, their interpretation, and protein domain impacted are described in Table 2. All KRAS and NRAS mutations observed in this study are located in the catalytic domain. The majority of these mutations are localized in the GTP binding site. For the BRAF gene, all mutations described in our study are found in the protein kinase domain. Predictions of the impact of each mutation are given in Table 2.

Discussion
In this study, we assessed the value of NGS-based testing to identify uncommon KRAS, NRAS and BRAF mutational profiles associated to in silico projection to evaluate their therapeutic and clinical implications in 779 samples of patients with mCRC. We identified 22 uncommon mutational profiles with a few missense variants that have not been previously reported in the literature, to our knowledge. Structural modelling of the observed missense variants in BRAF, KRAS and NRAS shows that most of the observed mutations can be accommodated in the protein structures without clearly adverse impacts on protein stability (Table 2). In a few cases, FoldX predicts the observed mutation to be highly destabilizing, primarily due to the introduction of interatomic clashes (e.g. KRAS Gly115Glu). While some such mutations may indeed lead to destabilization and inhibition of protein folding, it should be noted that the modelling procedure in FoldX does not consider backbone conformational changes. It is possible that some such mutations could be accommodated in the protein structure once backbone conformational changes are considered; an example is the RAS Gln61Leu mutation (not observed in this study) 14 . Despite this limitation, FoldX provides a quick, relatively accurate and parsimonious means to evaluate the impacts of mutations on protein stability.
We now briefly describe the results of modelling the previously unreported mutations. Modelling the Gly115Glu mutation in KRAS ( Fig. 1) suggests that this mutation has an adverse impact on protein stability. The three point mutations in KRAS that have not been previously reported (Lys88Asn, Cys118Tyr and Asp132Asn) are located on the protein surface and lead to minimal impacts on stability when KRAS is considered in isolation. The same is true of the NRAS Asp33Glu mutation, although this is located in the Switch 1 loop and is not directly involved in substrate binding. Examination of structures of RAS proteins in complex with RASGAP and SOS (PDB identifiers 1WQ1 and 1XD2, respectively) showed that most of these residues in KRAS and NRAS were not directly involved in the interaction with RASGAP and SOS. The exception is NRAS Asp33Glu, which may affect the interaction with SOS (Fig. 2) 15 . The BRAF Ser602Phe mutation leads to moderate stabilization relative to WT BRAF. This residue is located on the surface of the protein and is not in the BRAF dimer interface.
BRAF non-V600E mutations were detected in 1.43% of patients tested, which is consistent with previous published studies 16 . The predictive value of these mutations in the context of mCRC is indeterminate. However, exon 11 codon 469 BRAF mutations are located on the gene region coding for protein kinase function 11 . For exon 11 codon 466 mutations, the amino acid change is located within the glycine-rich loop in the kinase domain 17 . For samples #14 and #16, BRAF exon 15 codon 594 mutations are located on the gene region coding for protein kinase function and lead to impaired kinase activity, potentially conferring a favorable prognosis 18,19 . BRAF p.(Lys601Glu) (sample #21) is described as pathogenic but without published data on response to anti-EGFR mAbs therapy.
We now turn our attention to the 14 mutational profiles in our dataset with concomitant RAS and/or BRAF mutations. KRAS and BRAF mutations have been frequently described as mutually exclusive in CRC and concomitant KRAS and BRAF mutations are rare, occurring in less than 0.001% of cases 20 . In our study, in 10 samples with uncommon BRAF mutations, 6 are concomitant with a RAS mutation, and 1 BRAF hotspot mutation had a concomitant RAS mutation. The percentage of concomitant KRAS and BRAF mutations in our dataset is higher than previously described [21][22][23][24][25][26][27] . Since these studies only assessed BRAF codon 600 hotspot mutations, we may infer that hotspot KRAS mutations are more frequently associated with rare BRAF mutations than with hotspot mutations. This inference agrees with a previous study which found that patients with BRAF non-V600E mutations were more likely to have concomitant RAS mutations than patients with the BRAF V600E mutation 16  www.nature.com/scientificreports www.nature.com/scientificreports/ tumors have a different biology and natural history than KRAS or BRAF mutant tumors. We should illustrate the potential significance of this concomitant mutation on sample #12. Sample #12 bears a KRAS p.(Val14Ile) mutation described as pathogenic. Whether resistance to anti-EGFR mAbs is conferred by this mutation is unknown, however codon 14 belongs to the same domain as codons 12 and 13. As with mutations on codons 12 and 13, this mutation may be associated with a clinical resistance to anti-EGFR antibodies 2 . In this case, the knowledge of this rare mutation may be significant in guiding a therapeutic decision and to explain potential resistance to anti-EGFR antibodies. However, the presence of a BRAF non-V600 mutation could modify response to therapy and may even lead to potential treatment resistance. Even if first observations show that BRAF plays only a slight role in resistance to anti-EGFR mAbs, some BRAF non-hotspot mutations might contribute to reduced efficacy of anti-EGFR mAbs 18,28 .
In silico prediction of the functional effects of the observed mutations using both sequence-and structure-based approaches suggests possible biochemical mechanisms related to uncommon mutational profile and relevant in cancer. These in silico results warrant further in vivo examination to assess the relevance of detection of non-hotspot RAS mutations and their implication in resistance to anti-EGFR mAbs therapy. For the BRAF gene, non-V600E mutations may describe a novel subtype of mCRC with better prognosis, implying potentially different treatment management strategies. . The two oligo pools were hybridized to DNA samples. The specific hybridized targets were ligated, extended and PCR amplified with adaptors containing index with specific barcode sequences. Two complementary libraries were generated by targeting the forward and reverse DNA strands. The PCR-amplified amplicon libraries obtained were then purified using AMPure XP beads in order to remove non-specific products and reaction components.

Methods
Library DNA concentrations were quantified using Qubit 3.0 Fluorometer (ThermoFisher Scientific Inc, Massachusetts, USA) and their quality was assessed on Fragment Analyzer (Advanced Analytical, Ankeny, USA) using the Standard Sensitivity NGS Fragment Analysis Kit (Advanced Analytical). PCR product sizes have to be around 260 base pairs in length. All 616 validated libraries were normalized to enable similar amplification and sequencing levels for each sample library in the same run. Sequencing was performed according to the manufacturer's instructions. All libraries were pooled before sequencing on the MiSeq instrument (Illumina). Sequencing data analysis was performed on Sophia DDM ® software (Sophia genetics, Saint Sulpice, Switzerland). Reference sequences NM_033360.2 for KRAS, NM_002525.4 for NRAS and NM_004333.5 for BRAF were used for alignment and variant calling. Fourteen samples had insufficient coverage to be interpretable and 602 samples had interpretable sequencing data.

Uncommon mutational profiles. Uncommon mutational profiles were defined as i) a concomitant KRAS
and NRAS hotspot mutations, or ii) a KRAS, NRAF or BRAF non-hotspot mutation (associated or not associated with other mutations). A threshold of 1% allele frequency has been reported clinically relevant for KRAS mutations linked with lack of response to anti-EGFR therapy 29 .
For samples with an uncommon mutational profile, mutational status of MAP2K1 and HRAS genes (available in our gene panel and implicated in the MAPkinase pathway) have been also identified. www.nature.com/scientificreports www.nature.com/scientificreports/ In silico prediction of mutation impact. We used PolyPhen-2 (Polymorphism Phenotyping) and SIFT (Sorting Intolerant from Tolerant) scores to predict mutation impact on the protein, as well as FoldX to model the observed variants in protein structures [30][31][32][33][34] . PolyPhen-2 score predicts the possible impact of an amino acid substitution on the structure and function of a human protein. PolyPhen combines amino acid composition analysis in multiple sequence alignments with information from solved protein structures (where available). Sequence composition is evaluated using the Position Specific Independent Counts (PSIC) tool, which calculates sequence  www.nature.com/scientificreports www.nature.com/scientificreports/ profiles 35 . Differences in the profiles calculated by PSIC for the different variants at a site are indicative of damaging impact. In cases where structural data are available, the impact of a mutation is assessed by considering physicochemical properties such as residue size and hydrophobicity, and the maintenance of contacts with ligands, metals or other interacting proteins. Sequence-and structure-based features are combined to produce predictions using a naïve Bayes classifier. PolyPhen-2 produces scores between 0 and 1 along with annotations of whether the mutation is predicted to be benign or damaging.
SIFT also predicts whether an amino acid substitution affects protein function. SIFT scores are used to predict the damaging effect of nucleotide substitutions and frame shifts (insertions/deletions) on protein function For clarity, only residues whose atoms contact Glu115 are shown. The mutation to Glu at this position causes severe atomic clashes, primarily with residues Ile84 and Arg123, although there is potential for hydrogen bond formation with Arg123. Clashes are indicated by colored discs drawn between atoms. The color and size of the disc reflects the severity of the clash, with wider, redder discs indicating the most severe clashes. www.nature.com/scientificreports www.nature.com/scientificreports/ based on the maintenance of amino acid composition in alignments of the target sequence with closely related sequences. The SIFT server assigns scores for each residue from 0 to 1, where mutations with a score of ≤0.05 are predicted to not be tolerated, and mutations with score >0.05 are predicted to be tolerated 36 .
SIFT and PolyPhen-2 scores were determined using Sophia DDM ® software (version 5.0.7). Structural modelling of observed missense mutations was carried out using FoldX version 4 33,34 . FoldX estimates the impact of point mutations on the folding energy or stability of the protein using rigid-backbone modelling and a classical forcefield whose parameters are trained to reproduce experimental observations of mutational impacts on folding energy. Starting from the wildtype structures for KRAS, BRAF and NRAS (PDB identifiers 4OBE, 5VAM and 5UHV, respectively), the missense mutations observed in these proteins were modeled using the FoldX BuildModel tool, and the structural and energetic impacts of each mutation were assessed as the predicted change in folding energy (∆∆G). Synonymous (silent) and stop mutations were not modeled. All FoldX modelling results are available at https://doi.org/10.5281/zenodo.1467311.
Protein domain impacted. KRAS and NRAS protein domains impacted by the different mutations described were identified using UCSF Chimera 37 . The COSMIC database was used to identify protein domains impacted by BRAF mutations 38 .