Recurrent mutations of chromatin remodeling genes and kinase receptors in pheochromocytomas and paragangliomas

Rodrigo A. Toledo, Yuejuan Qin, Zi-Ming Cheng, Qing Gao, Shintaro Iwata, Gustavo M. Silva, Manju L. Prasad, I. TolgayOcal, Sarika Rao, Neil Aronin, Marta Barontini, Jan Bruder, Robert L. Reddick, Yidong Chen, Ricardo C. T. Aguiar, Patricia L. M. Dahia

Supplementary Appendix

List of contents

  1. Supplementary Methods
  2. Supplementary Table 1 (Table S1): features of the 41 samples of the next-generation sequencing (NGS) cohort
  3. Supplementary Table 2 (Table S2): mutations identified in known pheochromocytoma/paraganglioma-related genes in the NGS cohort
  4. Supplementary Table 3 (Table S3): features of the validation cohort of 136 pheochromocytoma/paraganglioma samples
  5. Supplementary Figure 1 (Figure S1): Alternative views of Figure 1 displaying novel variants and known PPGL mutations in our cohort
  6. Supplementary Figure 2 (Figure S2): Histone 3.3 G34W mutation in samples of a patient with metastatic paraganglioma and giant cell tumor of bone
  7. Supplementary Figure 3 (Figure S3): Histone 3.3 and G34W mutant in silico analysis
  8. Supplementary Figure 4 (Figure S4): Immunohistochemistry of trimethylated Histone 3
  9. Supplementary Figure 5 (Figure S5): Expression profiling of pheochromocytomas and paragangliomas with a H3F3A mutation
  10. Supplementary Table 4 (Table S4):Pathway analysis of G34‐mutant vs. wtpheochromocytomas and paragangliomas from microarray expression
  11. Supplementary Figure 6 (Figure S6): Familial pheochromocytoma segregating with a MET gene mutation
  12. Supplementary Table 5 (Table S5): Main clinical featuresand MET gene variants found in a validation cohort of 136 pheochromocytomas and paragangliomas

Supplementary Methods

Samples: Pheochromocytomas and paraganglioma samples were collected through an international consortium for the study of Familial Pheochromocytoma between 1993 and 2013. We obtained written informed consent or anonymized samples stripped of all personal identifiers through a tumor repository approved by the University of Texas Health Science Center at San Antonio (UTHSCSA) Institutional Review Board (IRB). Diagnosis of pheochromocytoma and/or paraganglioma was confirmed by histology in every case. Forty three samples from 41 individuals were used for whole exome (n=40) or transcriptome sequencing (n=3). Of the 41 individuals, thirty had pheochromocytomas, seven had paragangliomas and four had both tumors. Germline samples were selected from patients at high-risk for inherited PPGL, i.e. with a family history of PPGLs, early onset disease (<30 years old at diagnosis), multiple tumors, including bilateral pheochromocytomas and/or co-occurrence of other tumors which have been previously associated with pheochromocytomas or paragangliomas (Table S1). Seven patients (17%) had metastatic PPGLs. A separate cohort of 136 pheochromocytoma or paraganglioma samples was used for verification of detected mutations (general clinical and genetic features of this group are shown in Table S3). Genetic screen for KIF1B, MAX, RET, SDHB, SDHC, SDHD, TMEM127 and VHL, and in a subset of samples, SDHA, were previously performed in these samples. Tumor specimens collected during surgery were immediately frozen and stored in liquid nitrogen until processing. In some cases, tumor was available as formalin-fixed paraffin-embedded material; in these samples only germline, but not tumor, DNA was used for exome sequencing. Tumor and blood samples were obtained and processed using conventional methods. A summary of the main clinical features of this cohort are shown in Table S1. In one of the cases, multiple tissues were obtained: three fresh-frozen tumors, consisting of right (196) and left (197) pheochromocytomas and a retroperitoneal (bladder, 369) paraganglioma resected three years later were used for whole exome sequencing, while formalin-fixed, paraffin embedded (FFPE) sections of these three tumors as well as a paraaortic paraganglioma, giant cell tumor (GCT) of the tibia and normal gallbladder excised due to cholelithiasis were used for targeted sequencing analysis and immunohistochemistry. In addition, FFPE sections were also obtained from a retroperitoneal paraganglioma, liver metastases, normal adjacent liver, tibial GCT, femoral GCT and lung GCT metastases from a previously reported patient(1). We also obtained four anonymous GCTs, two tibial, one femoral and one fibular, from the UTHSCSA Pathology Tumor bank (median age, 34 years old). Diagnoses were confirmed by histology in every case. Expression profiling data from 39 of the PPGL cases has been previously reported (GEO accessions GSE2841 and GSE19987) (2, 3). Tissue samples were processed for DNA (all cases) and RNA (fourteen samples) using conventional methods, as reported(2).

Whole exome and transcriptome sequencing and variant detection. Whole-exome sequencing was performed in 40 and RNA sequencing was performed in 3 PPGL samples. Twenty six were tumors, 13 were germline samples and we also sequenced four matched tumor-germline pairs. Three of the tumors were from the same patient. The four paired samples were previously described(4). For the remaining samples, libraries were prepared using Agilent SureSelect XT2 Target Enrichment System for Illumina Multiplexed Sequencing, version C.2, and exomes enriched using Agilent Sure Select 44Mbp or NimbleGen 44Mb SeqCap EZ Exome Kit v2 (at the Beijing Genomics Institute) or Agilent Sure Select XT2 Human All Exon V5 capture (at University of Texas at Austin Genomic Core facilities). Paired-end 90- or 100-bp reads of enriched exomes were generated in Illumina HiSeq2000. Image analysis and base calling were performed using Illumina Real-Time Analysis Pipeline version 1.12.4.2 with default parameters. Alignment was performed with Burrows-Wheeler Aligner (BWA 0.5.5) and NextGEne (v2.3.4, SoftGenetics, State College, PA, USA) to NCBI Build 37. The mean depth of coverage of exomes was 107×. For RNAseq, poly-A mRNA was enriched from total RNA using MicroPoly(A) Purist kit (Ambion Life Technologies, USA). RNA-seq libraries were then created for each sample using RNA ligation as described elsewhere(5). The libraries were sequenced on Illumina HiSeq2500 according to the manufacturer's instructions. Alignment was performed with TopHat and NextGEne to NCBI Build 37 and RefSeq using default parameters. The average read number of the three RNAseq samples was 140Mbp. Variant annotation was performed with VarScan2 and NextGENe. The VarScan2 parameters were set to require minimum coverage of 10, minimum phred base quality of 20, allele frequency of at least 10% for tumors and 25% for germline samples, and a p-value of <0.05(6, 7). Similar parameters were applied for NextGEne, and variants with an overall score (calculated based on read coverage, directional allele balance, homopolymer and mismatch allele) greater than 10 were selected(8). Synonymous changes, intronic variants located more than 6 nucleotides from the splice junctions, variants that coincided with the alignment of known pseudogenes and unmapped variants were removed. The estimated pathogenicity of missense variants was assessed by five different prediction programs (Polyphen2, SIFT, PhyloP, MutTaster and LRT); variants assigned scores within the ‘pathogenic’ range (>0.85) by at least two of these algorithms were kept for further analysis. We excluded variants reported with minor allele frequency (MAF) greater than 0.05 in the dbSNP135, 1000 Genome Project, NHLBI GO Exome Sequencing Project (ESP) and Exome Aggregation Consortium (EXAC). To prioritize relevant mutations, missense variants were analyzed by CRAVAT (cancer-related analysis of variants toolkit). We applied two tools of CRAVAT for this analysis: CHASM (cancer-specific high-throughput annotation of somatic mutations) was used for tumor samples, with each of the reference tumor cohorts input(9), while germline variants were analyzed by both CHASM and by VEST (variant effect scoring tool), designed specifically for inherited diseases(10). The driver score feature of each program was used to prioritize variants likely to be significantly involved in tumorigenesis: scores >0.6 for VEST (highest scores associated with pathogenic variants) and <0.4 (lowest scores associated with pathogenic scores) for CHASM were deemed significant. Finally, selected variants were annotated by Oncotator (v.1.5.3.0) and further filtered to keep variants that were absent in the reference databases listed above, present in COSMIC and retained ‘damaging’ or ‘deleterious’ status in PolyPhen2, Mutation Taster and SIFT. Mutations considered pathogenic by these criteria were validated visually using NextGEne viewer. A subset of these mutations was further verified by Sanger sequencing as described below. To examine mutation co-occurrence vs. mutual exclusivity of the detected variants, we calculated the log odds ratio, as reported(11). For these analyses, we compared both individual tumors and the full cohort for variants in 1) novel classes (i.e., chromatin-related and kinase), 2) known PPGL, 3) in both classes, or 4) in neither one. The significance (P value) of co-occurring or anti-co-occurring cancer gene mutations was determined using Fisher's exact test based on these numbers.

Targeted amplicon-based next-generationsequencing. DNA from fresh-frozen or FFPE from the various tumor and nontumor tissues from two cases carrying a H3F3A G34W mutation, as well as three unrelated PPGLs, and four sporadic GCTs were used to quantify the variant allele representation at ultra-deep sequencing using IlluminaTruSeq Custom Amplicon with slight modification. A similar approach was used to detect the METc.2416G>A; p.V806M variant in samples from three individuals with familial pheochromocytoma without a known driver mutation, including germline DNA from two affected siblings and FFPE from a nonpheochromocytoma tissue section (bladder) from their affected father, along with unrelated samples. Primers were designed to span a 170 bp flanking the H3F3A-G34W mutationand 150bp surrounding the MET-V806M variant, and contained adaptors to allow for subsequent sequencing on the IlluminaMiseq, with interspersed ‘1N’ base between the gene-specific sequence and the adaptors to add sequence diversity given the small targeted genomic region. Each amplicon was purified using AMPure XP beads (Beckman Genomics) and used as template for a second round of PCR using indexed primers from Illumina FC-131-2001-Nextera® XT Index Kit v2 Set. Indexed PCR products were again purified with AMPure XP beads and quantified using a Qubit (Life Technologies). Libraries were multiplexed at equimolar amounts with 5%PhiX using and run on the IlluminaMiSeq150 cycle V3 kit. Sequencing reads in fastq format were aligned to the human reference genome hg19 using Burrows-Wheeler Aligner (BWA)(12), followed by sorting and mpileupwithSAMtools(13). mpileup files were used on VarScan2(6) for annotation. The H3F3A amplicon hadan average read depth of 67700x (±47969SD), as shown in Fig.1b. The MET target fragment had average read counts of 135,602±4,227 reads (Suppl. FigSb). Variant allele counts below the quality threshold (1%) were considered negative.

Sanger sequencing.Primers are available upon request. Sequencing was processed at Beckman Genomics and analyzed with Mutation Surveyor (Softgenetics, PA), as previously reported(4). The Mutation Quantifier tool of Mutation Surveyor was used to measure frequency, expressed as % of the mutant allele in the tumor samples. This measurement is based on the intensity ratio (degree of drop in the reference allele peak and the gain of the variant allele at that position) of the variant of interest in relation to the intensity of the neighboring peaks of the same nucleotide as the variant:, where ‘i’ represents the number of neighboring peaks of the same nucleotide as the variant allele. Peaks immediately preceding or after the target variant are not included in this calculation as those can be affected by the mutation.

Differential gene expression analysis of PPGLs. Gene expression data generated using the Affymetrix U1332.0 platform were previously reported (GEO accessions GSE2841 and GSE199877)(2, 3). Normalized data from the three tumors with H3F3A G34W mutation and from 36 tumors without this mutation but with a detectable mutation in other pheochromocytoma/paraganglioma susceptibility gene were used for unsupervised clustering using dCHIP, and applying the following parameters: 1 <average SD/mean across categories <1000 and were present in >=50% of samples at signal level >= 20.00, as previously reported(2). In addition, we analyzed categorized samples by gene set enrichment analysis (GSEA). The two sample categories (G34W and WT) were processed using GSEA v2.0, with gene set as permutation type, 1000 permutations, and log2 ratio of classes as metric for ranking genes(14) (Suppl. Table 2). DAVID Bioinformatics Resources 6.7, NIAID/NIH toolset was used for further annotation of the pre-ranked list of differentially expressed genes identified in GSEA based on gene ontology terms, protein-protein interactions and protein functional domains(15). An enrichment score and associated p value and were generated. Publicly available data from the glioma expression data (GSE34824(16)) were retrieved from the Gene Expression Omnibus ( and analyzed in GSEA and DAVID using the same parameters as those applied to the PPGL samples.

Structure Modeling.Structural predictions of histone 3.3 variants were performed using the I-TASSER server(17). I-TASSER is a fully automated tool that selected the nucleosome PDB structure 1KX5, determined by X-ray crystallography at 1.9Å resolution, as the highest significance template for predicting the final models. We used the full length amino acid sequence for human histone H3F3A (ENST00000366813) and computationally modeled the wild-type (WT) and the G34W mutant. The C-scores (confidence score, range -5 to 2), TM-scores (template modelling) and RMSD values are shown in Suppl. Figure S2a for the WT and G34W structural models provided by I-TASSER. The scores indicate a high quality of the predictions and correct topologies of the models based on the posterior alignment with the template. The model with the highest score was selected for further analysis. 3D structure visualization and structure superimposition were performed using Pymol (DeLano Scientific). Electrostatic surface potential was calculated using Adaptive Poisson-Boltzmann Solver (APBS) within Pymol and hydrophobicity analysis was based on the scale defined by Eisenberg et al(18).

Western Blots. Whole cell lysates from tumors or cell lines were prepared as previously described(3, 19). Fifty-μg(or 5μg, for histone analysis) whole cell lysates were run on 12% SDS-polyacrylamide gels, transferred to PVDF membranes (Immobilon, Millipore) and hybridized with the antibodies listed above according to the manufacturer’s instructions. Filters were developed with a chemiluminescence assay (West Pico, Pierce) and exposed to X-Ray film (Kodak).Proteins were detected with the following antibodies: H3K27me3, H3K36me3, Histone H3, MYCN, MERTK, phospho ERK (all from Cell Signaling Technology) or β-actin (Sigma-Aldrich, #A2228). PVDF membranes were stripped with RestoreTMPlus Western Blot Stripping Buffer (Thermo Scientific, #46430) and re-probed with Histone H3 and β-actin antibodies for loading control.

References

1.Iwata S, Yonemoto T, Ishii T, Araki A, Hagiwara Y, Tatezaki S-i. Multicentric Giant Cell Tumor of Bone and Paraganglioma: A Case Report. The Journal of Bone and Joint Surgery. 2013;3:e23.

2.Dahia PL, Ross KN, Wright ME, Hayashida CY, Santagata S, Barontini M, et al. A HIF1alpha regulatory loop links hypoxia and mitochondrial signals in pheochromocytomas. PLoS Genet. 2005;1:72-80.

3.Qin Y, Yao L, King EE, Buddavarapu K, Lenci RE, Chocron ES, et al. Germline mutations in TMEM127 confer susceptibility to pheochromocytoma. Nat Genet. 2010;42:229-33.

4.Toledo RA, Qin Y, Srikantan S, Morales NP, Li Q, Deng Y, et al. In vivo and in vitro oncogenic effects of HIF2A mutations in pheochromocytomas and paragangliomas. Endocr Relat Cancer. 2013;20:349-59.

5.Podnar J, Deiderick H, Huerta G, Hunicke-Smith S. Next-Generation Sequencing RNA-Seq Library Construction. Current protocols in molecular biology / edited by Frederick M Ausubel [et al]. 2014;106:4 21 1-4 19.

6.Koboldt DC, Larson DE, Wilson RK. Using VarScan 2 for Germline Variant Calling and Somatic Mutation Detection. Current protocols in bioinformatics / editoral board, Andreas D Baxevanis [et al]. 2013;44:15 4 1- 4 7.

7.Koboldt DC, Zhang Q, Larson DE, Shen D, McLellan MD, Lin L, et al. VarScan 2: somatic mutation and copy number alteration discovery in cancer by exome sequencing. Genome Res. 2012;22:568-76.

8.Miyake N, Elcioglu NH, Iida A, Isguven P, Dai J, Murakami N, et al. PAPSS2 mutations cause autosomal recessive brachyolmia. J Med Genet. 2012;49:533-8.

9.Douville C, Carter H, Kim R, Niknafs N, Diekhans M, Stenson PD, et al. CRAVAT: cancer-related analysis of variants toolkit. Bioinformatics. 2013;29:647-8.

10.Carter H, Douville C, Stenson PD, Cooper DN, Karchin R. Identifying Mendelian disease genes with the variant effect scoring tool. BMC Genomics. 2013;14 Suppl 3:S3.

11.Ciriello G, Cerami E, Sander C, Schultz N. Mutual exclusivity analysis identifies oncogenic network modules. Genome Res. 2012;22:398-406.

12.Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009;25:1754-60.

13.Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics. 2009;25:2078-9.

14.Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102:15545-50.

15.Huang da W, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4:44-57.

16.Schwartzentruber J, Korshunov A, Liu XY, Jones DT, Pfaff E, Jacob K, et al. Driver mutations in histone H3.3 and chromatin remodelling genes in paediatric glioblastoma. Nature. 2012;482:226-31.

17.Yang J, Yan R, Roy A, Xu D, Poisson J, Zhang Y. The I-TASSER Suite: protein structure and function prediction. Nat Methods. 2015;12:7-8.

18.Eisenberg D, Schwarz E, Komaromy M, Wall R. Analysis of membrane and surface protein sequences with the hydrophobic moment plot. J Mol Biol. 1984;179:125-42.

19.Duncan EM, Muratore-Schroeder TL, Cook RG, Garcia BA, Shabanowitz J, Hunt DF, et al. Cathepsin L Proteolytically Processes Histone H3 During Mouse Embryonic Stem CellDifferentiation. Cell. 2008;135:284-94.

1