Somatic Mutation Profiles of MSI and MSS Colorectal Cancer Identified by Whole Exome Next Generation Sequencing and Bioinformatics Analysis

1Colorectal cancer is the third most common cancer with about 1 million cases worldwide. Over the last decades it has become clear that CRC evolves through multiple pathways and that these pathways can be roughly defined on the basis of molecular patterns such as the integrity of the mismatch repair system (MMR) or mutational and epigenetic patterns. Deficiency in the MMR is reflected in DNA microsatellite instability (MSI) which has also been associated with treatment outcome, but which needs to be further validated in additional clinical studies,,,,,.

2High-throughput Sanger sequencing studies on the other hand have shown that the mutation frequency of candidate cancer genes might be much higher than expected, and that the particular combination of mutations might influence the tumor's properties,,,,. With the development of next-generation sequencing (NGS) technologies the sequencing throughput has dramatically increased and the costs have decreased. In addition, and especially important for clinical settings, NGS can be applied to formalin-fixed and paraffin embedded FFPE tissue material as well as highly degraded DNA which is routinely prepared in pathology departments or found in ancient DNA,. Several studies have used NGS technologies for the identification of the underlying mutation in monogenetic diseases,. However, only a limited number of studies report on next-generation sequencing to identify new candidate cancer genes; one of the earliest studies examined cytogenetically normal acute myeloid leukemia, and breast cancer genomes,. In addition, studies on malignant melanoma and small-cell lung cancer cell lines have provided first insights into genomic alterations induced by ultraviolet light exposure or tobacco smoke,.

3To gain insight into the genomes of microsatellite stable and instable colorectal cancers and to identify functional relevant mutational patterns we used a hybridization based whole exome DNA capturing approach followed by 454 next generation sequencing. Applying stringent bioinformatics analyses, we narrowed down the amount of functionally significant somatic mutations in MSI to 359 and 45 in MSS cancers, thus highlighting specific mutation patterns depending on the microsatellite status. We were able to confirm our results by sequencing the exomes of four additional CRC cases (one MSI, three MSS) using a different enrichment and sequencing technology. Among these mutations are BRAF in the MSI cancer and KRAS and TP53 in the MSS cancer, further underscoring the validity of our selection approach. Further functional characterizations identified recurrent somatic mutations in BMPR1A, a protein which has been associated so far with juvenile polyposis syndrome, a cancer predisposition syndrome.

4We sequenced tumor and matching normal colon tissues from two patients with high grade adenocarcinoma of the colon (G3), patient 1 with a microsatellite instable and patient 2 with a microsatellite stable tumor (Table 1, Figure S1). For the determination of germline mutations we sequenced in addition to the tumor tissues from each patient adjacent not affected normal colonic tissue. Using Illumina sequencing and SNP arrays we determined that the tumor of patient 1 is copy number stable whereas patient 2 showed variations which we used for the re-evaluation of identified high stringency mutations (Table S3).

5We analyzed the complete exomes of more than 135,000 exons with single-read shotgun 454 sequencing (Figure 1, Figure S1, Figure S2, Table 2). To assess the effect of coverage depth on the sensitivity and specificity of sequence variant detection, genotype calls of the Affymetrix SNP array 6.0 were compared step-wise to the called nucleic acid positions and resulted in an accuracy of more than 99% (Figure S3). In addition to the SNP array, we used Sanger sequencing to confirm 23 selected mutations (Table S5, Figure S5).

6Searching for variants in coding regions we found 12,767 and 12,518 small nuclear variations in 6,428 and 6,205 genes, for tumor and normal respectively. Of these variants 1,428 for control and 2,404 for tumor have an average heterozygosity or minor allele frequency lower than 1% or have not been previously reported in dbSNP or the 1000 Genomes Project (Figure 2). Since indels (small insertions and deletions) at homopolymeric sites are a major source of sequencing errors of the 454 platform we ignored this type of alteration in our analyses.

7Our somatic variant detection strategy was designed to minimize false positive somatic variant calls rather than to determine zygosity. We used a two-step approach with two different stringency levels to detect variants in tumor and benign tissues similar to Pleasance et al.. In the first step, tumor variants were called under stringent criteria, which were determined by comparison to the SNP genotyping array data. The second step ascertained whether the tumor variant was germline or somatic. To keep the false negative rate in the benign tissue at less than 10%, we set the coverage cut-off at 5-fold, below which no conclusions were drawn regarding whether a variant was somatic or germline. If the coverage cut-off was met, a single read showing the same variant in benign and tumor tissue resulted in the categorization of the variant as germline.

8Using this strategy we identified 915 somatic non-synonymous mutations affecting 864 genes. The majority of somatic mutations were missense mutations (65%). However, many (7%) are located within untranslated regions of genes and might therefore result in altered expression or increased decay of mRNA species. Furthermore, approximately 0.5% of these mutations are found at splice sites and could influence splicing events, leading to an altered transcriptome structure. In addition, three somatic variants were identified in miRNA regions (Table S6). These mutations are of particular interest because miRNAs have been implicated as master regulators of tumor homeostasis. Analysis of the specific types of nucleic acid variations, including known and unknown variants, showed essentially the expected rates of nucleotide exchanges, as determined by calculations using dbSNP130 (Figure S4).

9Since not all of the 1,304 somatic mutations are likely to be pathologically relevant, we sought to identify those that probably destroy protein function or affect highly conserved amino acids and might therefore be functionally important. We used Polyphen and MutationTaster classification tools to predict the functional consequences of amino acid changes or frameshift mutations and found that 359 genes had at least one potentially destructive mutation (Table S1),.

10Of the potentially destructive somatic mutations, 309 were located in positions highly conserved in 44 different species, including opossum (Monodelphis domestica), chicken (Gallus gallus) and lamprey (Petromyzon marinus). Of these, 259 were located in genes expressed in the colon, of which 47 were repair, receptor, or kinase genes. Visualization of selected mutations on protein structures indicates that these nucleotides are on the protein surface, potentially resulting in disrupted protein-protein interactions.

11The Catalog of Somatic Mutations in Cancer (COSMIC) database is a comprehensive collection of cancer-related mutations. Approximately 39% of our mutated genes were already described in this database, and 13% were found by Wood et al. As did these previous databases, we found the BRAF p.V600E mutation in the MSI case and we identified KRAS and TP53 mutations in the MSS tumor. BRAF mutations are found in approximately 10% of CRCs, predominantly MSI and 30 to 35% of all patients with sporadic colorectal cancers carry somatic KRAS and TP53 mutations. These findings further demonstrate the sensitivity of our classification strategy.

12For the MSS colon cancer we identified 10,622 small nuclear variations. After the same filtering processes as for the MSI cancer using the dbSNP database and the data from the 1000 Genomes Project 1,288 variants remained which either had low prevalence or were unknown. Of these, 198 were somatic and 45 were predicted to alter gene function based on MutationTaster and Polyphen calculations (Figure 2B, Table S2),. In regard to copy number variations five of the 45 identified mutations are located within amplified regions, and, as expected, none in regions with deletions. The ratios of reads with reference sequence to mutated sequence are not exceeding ratios in copy number stable areas which supports the SNV-calling algorithm.

13In contrast to 1,304 somatic mutations in the MSI tumor we found 198 somatic mutations in the MSS tumor which demonstrates that the defective MMR system in MSI tumors results in a significant increase in mutation rates in colorectal cancer.

14Furthermore, looking at intersections between both cancer types we found BMPR1A, WDTC1 (WD and tetratricopeptode repeats 1) and EHD3 (EH-domain containing 3) mutated in both tumors. The selection was based on functional impairment with high probability in Polyphen and MutationTaster,. All mutations are located on the surface of the protein and are highly conserved. Since we found significant cancer-related pathways associated only with BMPR1A but not with WDTC1 or EHD3, and in addition germline mutations in BMPR1A are a known risk factor for juvenile polyposis syndrome, we chose BMPR1A for additional functional assays (Figure 3, Figure S5). Using reporter assays with wild type and mutated BMPR1A proteins we were able to show that the mutated proteins are strongly impaired in their signalling function and that stimulation with BMP2 results in a reduced maximum activity (Figure 3D).

15The analyses presented so far have been based on 454 whole exome sequencings of one colorectal cancer patient for each microsatellite status. To further confirm that the increased amount of coding mutations in MSI cancers can be generalized and is not due to the technology used (‘array’ enrichment and 454 sequencing) we sequenced the exomes of four additional colorectal cancers, each with matching normal tissues. This time we used ‘in solution hybridization’ for capturing of DNA followed by SOLiD sequencing. After the same filtering procedures as described for the first two patients we again determined up to 8-fold higher mutation rates for the MSI colorectal cancer than for MSS cancers (Table 3). In this regard we found 532 non-synonymous somatic SNVs in the additional MSI CRC and only 65, 74 and 76 in the three MSS CRC cases. Thus, the differences in mutation rates are reproducible and independent of the sequencing technology used.

16Using next-generation sequencing, we sequenced the exomes of MSS and MSI colon cancer patients, with mean coverages of approximately 20-fold. We applied a two-sided classification algorithm to uncover functionally relevant mutations. Using this approach, we demonstrate for the first time that an array-capture NGS one-step work flow is a powerful tool for deep characterization of solid tumors and show that MSI tumors carry eight times more functional relevant mutations than MSS tumors (Figure 2).

17The functional impact of the somatic variations was predicted using two functional prediction algorithms, Polyphen and MutationTaster, and we found 359 somatic mutations for the MSI and 45 for the MSS cancer that are highly likely to cause functional impairment,. The heterogeneity of mutated genes suggests that not a specific gene per se but the affected pathway plays a major role for tumor development. In this regard we find a significant enrichment of mutations in cancer-related pathways such as cancer, cellular development and DNA replication, recombination and repair (Table S5). Interestingly, we find 50% of the most significant enriched pathways in the MSS cancer also as significantly enriched pathways for the MSI cancer, indicating that even though MSI cancers harbour an increased number of mutations both cancers might develop through overlapping pathomechanisms.

18Historically, microsatellite testing in colorectal cancers was the first predictive test for the identification of an underlying mismatch repair (MMR) mutation. Since more than 90% of hereditary non-polyposis colorectal cancers (HNPCC) show MSI, the microsatellite status has become a common diagnostic marker of MMR. In addition, survival advantages and therapeutic consequences have been reported for patients with MSI tumors,. Using ultradeep sequencing of one conserved region, UCR41, de Grassi and colleagues show that this region has higher mutation rates in HNPCC samples than in healthy controls and suggest that this might be used as sensitive molecular assay of genomic instability. We have extended their analyses to whole exomes and found 8-fold differences in the numbers of somatic mutation of MSI and MSS colorectal cancers. In comparison, a study by Greenman et al. which reported on the sequencing of 518 protein kinase genes in 210 diverse human cancers found an approximately 25-fold higher mutation rate for MMR-deficient cancers. However, these are extrapolations and in contrast to our study they included tumors from different origins and examined all somatic mutations irrespective of their functional relevance. With our sequencing approach we also detected a somatic MLH3 mutation in the MSI tumor, which, even though MLH3 mutations do not belong to the classical MMR mutations in CRC, might contribute to the microsatellite instability phenotype. Furthermore, the combination of MSI and BRAF mutation, as detected for the MSI tumor described, is most frequently found for CpG island methylation phenotype 1 (CIMP1) tumors which are associated with MLH1 promoter methylations,. The promoter methylation in turn is associated with gene silencing mechanisms which is suggestive as an explanation for the MSI status of the tumors. On the other side, it has been proposed that chromosomal instability (CIN) and CIMP represent two independent and inversely related mechanisms of instability. CIMP-negative cases are associated with p53 (71%) and KRAS (33%) mutations, but are rarely found with BRAF (2%) mutations. Since we found large copy number variations as well as KRAS and TP53 mutations in the MSS tumor analyzed this tumor is most likely CIMP-negative,.

19The sequencing and analysis strategy we have presented might be the basis for future classification tools for colorectal cancers because it may allow a parallel detection of an increased mutation frequency in MSI tumors as well as the detection of the underlying MMR defect. In addition, we were able to detect mutations of genes frequently associated with certain subtypes of colorectal cancers such as BRAF, KRAS and TP53. Within our high priority genes, encompassing all genes which pass all selection filters, BMPR1A stands out as mutated in both cases. The overall structure reveals that the mutated amino acids are all located at the C-terminal intracellular helix bundle at the protein kinase domain and suggests that protein-protein interactions are destroyed. Germline BMPR1A mutations predispose to juvenile polyposis syndrome; however, our findings indicate that also somatic mutations might play an important role in sporadic colorectal cancer development,,,,.

20Besides these mutations we have also identified several mutated cancer drug targets or genes that are associated with treatment outcome, including BRAF, KRAS, FGFR2 and MTOR, which might help to choose optimal drug combinations. As such similar targeted re-sequencing approaches and bioinformatics filtering strategies might become a gold standard for individually tailored colorectal cancer treatment in the future.

21The study has been approved by the Ethical Committee of the Medical University of Graz. For new samples patients have given their written informed consent. For old samples (15 years old) no informed consent was available, therefore all samples and medical data used in this study have been irreversibly anonymized.

22Patient 1 had a high-grade (G3) adenocarcinoma of the proximal colon, staged pT-3C, pN-0, pM-X, pR-0, microsatellite instable (Figure 1A). In addition, this case was selected because of its chromosomal stability, as determined using genome-wide next generation sequencing (NGS) (Figure 1B). Patient 2 had a high-grade (G3) adenocarcinoma of the proximal colon, staged pT-4B, pN-2, pM-X, microsatellite stable.

23Human tissue obtained during surgery was snap-frozen in liquid nitrogen. Cryosections (3 µm thick) were prepared and stained with haematoxylin and eosin to evaluate tumor cell content. Dissections were performed under the microscope to achieve a tumor cell content of >80%. DNA isolation was performed using the QIAamp DNA Mini Kit (Qiagen), according to the manufacturer's instructions.

24Genomic DNA of both tissues was subjected to whole exome sequence capture using Roche/NimbleGen's 2.1M Human Exome Array. This array is based on build 36.3 of the human genome sequence, and captures the coding regions of 16,755 NCBI RefSeq genes (approximately 180,000 coding exons) as well as 493 miRNA regions. Tumor and normal tissue DNA were subjected to whole exome sequence capturing according to the manufacturer's protocol. DNA was sheared by nebulization to fragment sizes below 800bp, cleaned (Zymo Research) and end-polished using T4 DNA Polymerase and T4 Polynucleotide Kinase. Linker adapters pSel3 (5′ – CTCGAG AAT TCT GGA TCC TC – 3′) and pSel4-P (5′ – Phos/GAG GAT CCA GAA TTC TCG AGT T – 3′) were ligated and size selection was performed using AMPure DNA Purification Beads (Agencourt). Quality was controlled with the Bioanalyzer system. LM-PCR was performed with LMPCR3 primers (5′ – CUC GAG AAU UCU GGA UCC UC – 3′) before the library was used for hybridization at 42°C for 72h. The arrays were washed two times at 47.5°C, two times at room temperature and two times at 42°C with washing buffers as recommended. Bound genomic DNA was eluted with 125 mM NaOH for 10 min at room temperature and amplified by LM-PCR using primers LMPCR3. Captured amplified samples were subjected to quantitative PCR to measure the relative enrichment.

25The enriched Nimblegen DNA was used to construct single-stranded Genome Sequencer FLX (454/Roche) libraries. After emulsion PCRs sequencing primers were annealed to the template and beads were incubated with Bst DNA polymerase, apyrase, and single-stranded binding protein. Pyrosequencing was performed on a 70×75 mm picotiter plate in 13 separate sequencing runs. After default raw data processing, a resequencing trimming filter was used to increase the data output. (Parameters used: doValleyFilterTrimBack = false, vfBadFlowThreshold = 6, vfLastFlowToTest = 168, errorQscoreWindowTrim = 0.01).

26For the sequencing we performed 13 Genome Sequencer FLX runs, which produced over 558 million bases and 1.43 million reads per run. Reads were aligned to the human reference genome, NCBI build 36 (http://hgdownload.cse.ucsc.edu/goldenPath/hg18/), using GS Reference Mapper Version 2.0.0.12 (Roche). The best matches in the genome were used as the location for the reads with multiple matches. Only unique reads with a minimum length of 50 bp were used for further analysis (see run statistics Tab. 2).

27Detection of variants was performed with the GS Reference Mapper Version 2.0.0.12 (Roche). Redundant reads were subtracted before variant callings. Only the HCDiff (high confidence differences) of the GS Mapper software were used as basis of variant detection. HCDiff callings presume at least three reads with the variant with both forward and reverse reads included; alternatively the quality scores at the variable positions must be over 20 (or over 30 if a homopolymer of five or more bases is involved). As additional quality criteria we used only variants with a coverage >10× of high quality reads.

28Enrichments and SOLiD library preparation were performed according to Agilent's SureSelect Target Enrichment protocol for the Applied Biosystems SOLiD system. In brief, whole genomic DNA was sheared and end repaired. For adapter ligations 30x excess of the adapters were used. Size selections for 150–200 bp DNA fragments were performed followed by a nick-translation and amplification step with Platinum polymerase (Invitrogen) and Pfu-Polymerase (Fermentas). For hybrid selection the libraries were adjusted to 500 ng in 3.4 µl volume and added to the SureSelect Block solutions. Hybridizations were performed for 24 h at 65°C, hybrids were extracted with 500 ng M-280 streptavidin Dynabeads (Invitrogen) and finally eluted with 50 µl Elution buffer. After amplification with Platinum polymerase the libraries were quantified by qPCR and DNA concentration was titrated to achieve a fraction of 10–20% monoclonal template beads in the emulsion PCR using in total 0.7 to 1 billion beads. Successive bead enrichment and deposition of 130 million beads per quarter slide (quad) was followed by standard 50 bp fragment runs. Each of the four patient samples was analyzed on a single quad.

29Mapping was performed with the Bioscope alignment pipeline using the seed & extend algorithm with a mismatch penalty score −2.0. Single Nucleotide Variants (SNV) were called with the DiBayes algorithm integrated in the Bioscope package.

30Tissue materials were genotyped on the Affymetrix 6.0 array, according to the manufacturer's protocol. Array positions with a quality score (p-value) <0.1 were used for comparison with the sequencing data. Sequencing data positions were used if their coverage exceeded 3-fold. This generated 46,000 and 49,000 positions for tumor and benign tissue, respectively, that were eligible for comparison. To determine false positive and false negative rates, we set the array data as standard and distinguished between reference call and SNP call dependence on the array data.

31For the detection of somatic variants, a bimodal strategy was applied with tumor variants called under stringent criteria, whereas variants in control tissue were called using less stringent criteria: A minimum threshold for reads was set with variants of 15% of all reads at a given position in tumor. Less stringent criteria were used for calling control tissue variants with a minimum of one variant read and a minimum coverage cutoff of 5. For coverages above 30-fold one variant read was accepted.

32Follow-up confirmation of identified SNVs was performed on an ABI 3730 (Applied Biosystems) capillary sequencing instrument following standard procedures.

33Preparation of single read libraries and sequencing were performed using the Solexa sequencing platform (GenomeAnalyzer IIx, Illumina) following the manufacturer's instructions. Image Analysis and base calling were performed using Firecrest 1.9.5_14 and Bustard 1.9.5_14 and reads were aligned to the human genome (NCBI36) using Bowtie 0.9.7.1. Copy number analysis was done in R using the DNAcopy package. In short, DNA read frequencies were determined for bins of 50 Kb. The log2 frequency ratio of corresponding bins was calculated for tumor versus normal tissue. Median of ratios was centered to zero experiment wise. Log ratios were smoothed by DNAcopy using default values and copy number variation was detected by DNAcopy using a threshold of two standard deviations.

34Genotyping on the Affymetrix 6.0 array was performed according to the manufacturer's protocol. Regions of copy number gain and loss were determined by paired and analysis using the Hidden Markov Model (HMM) of the Partek Genomics Suite software (Partek Inc, St.Louis, MO) with default parameter settings. For paired analysis, copy number values were generated by comparing tumor and benign tissue profiles from the same patient.

35For each tissue, variations were annotated using the gene models generated by Ensembl (ensembl 54.36, www.ensembl.org). All variations were mapped to all transcript models, which led to multiple annotations for several loci. For instance, a variant can lead to an amino acid change in one transcript and appear in the UTR of another. Comparison of tumor and benign tissue variants to dbSNP130 and 1000 Genomes Project data was carried out with the subset of dbSNP130 and 1000 Genomes Project positions with minor allele frequencies or average heterozygosity >0.01. Variants were subjected to many comparisons with external data sources. Most data sources are integrated in the UCSC genome browser (http://hgdownload.cse.ucsc.edu/goldenPath/hg18/database/) or were derived from websites like the 1000 Genomes Project (ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/release/2009_04/) the gene ontology data (http://archive.geneontology.org/full/2009-10-01/go_200910-termdb.obo-xml.gz), the cosmic database version 46 (http://www.sanger.ac.uk/genetics/CGP/cosmic/) or the cancer gene census database (http://www.sanger.ac.uk/genetics/CGP/census/). Functional classifications were performed using Polyphen and MutationTaster classification tools (http://genetics.bwh.harvard.edu/pph/, http://neurocore.charite.de/MutationTaster/),. Base conservation among 44 species was tested using the phyloP track of the UCSC browser (http://genome.cshlp.org/content/early/2009/10/26/gr.097857.109.abstract). Bases were considered highly conserved if their conservation score was greater or equal 2.0 (0.975 quantil of all conservation scores). For gene expression healthy control samples from http://www.ncbi.nlm.nih.gov/sites/GDSbrowser?acc=GDS2609 were used. We calculated genewise mean expression values across all samples and used the first quartile as threshold to determine gene expression. Pathway analyses were performed with the ingenuity pathway analysis tool (http://www.ingenuity.com). All new data from this study has been deposited at NCBI dbSNP database of genetic variation (user-name MPIMGCancerogenomics). Accession numbers are included as Table S4.

36Models for BMPR1A were obtained from SwissModel and ModBase. Very similar models except for some loops (total Calpha rmsd 0.66) were rendered in PyMol,.

37The activity of the wildtype mouse protein (wtmBMPR1A) and its mutants was determined by measuring induced Luciferase activity in the transiently transfected pre-myoblastic mouse cell line C2C12 (ATCC). WtmBMPR1A was amplified from mouse cDNA and cloned in the expression vector pCS2+. Both mutations W487R and E502G were inserted by Quikchange mutagenesis (Stratagene) using the following primer pairs:

38mBmpr1a_W487R_fwd caatcgtgtctaaccgcCggaacagcgatgaatg;

39mBmpr1a_W487R_rev cattcatcgctgttccGgcggttagacacgattg and

40mBmpr1a_E502G_fwd gttttgaagctaatgtcagGatgttgggcccataatc;

41mBmpr1a_E502G_rev gattatgggcccaacatCctgacattagcttcaaaac.

42C2C12 cells were cultured in DMEM glucose 4,5 g/L with 10% FCS were co-transfected with each Bmpr1a expression construct, a Smad Binding Element (SBE) luciferase construct and the normalization vector pRL-Tk (Promega Corporation, Madison, WI, USA) using Turbofect (Fermentas GmbH, St. Leon-Rot, Germany). Luciferase activity was determined as described previously.