Driver gene alterations and activated signaling pathways toward malignant progression of gastrointestinal stromal tumors.

Mutually exclusive KIT and PDGFRA mutations are considered to be the earliest events in gastrointestinal stromal tumors (GIST), but insufficient for their malignant progression. Herein, we aimed to identify driver genes and signaling pathways relevant to GIST progression. We investigated genetic profiles of 707 driver genes, including mutations, gene fusions, copy number gain or loss, and gene expression for 65 clinical specimens of surgically dissected GIST, consisting of six metastatic tumors and 59 primary tumors from stomach, small intestine, rectum, and esophagus. Genetic alterations included oncogenic mutations and amplification-dependent expression enhancement for oncogenes (OG), and loss of heterozygosity (LOH) and expression reduction for tumor suppressor genes (TSG). We assigned activated OG and inactivated TSG to 27 signaling pathways, the activation of which was compared between malignant GIST (metastasis and high-risk GIST) and less malignant GIST (low- and very low-risk GIST). Integrative molecular profiling indicated that a greater incidence of genetic alterations of driver genes was detected in malignant GIST (96%, 22 of 23) than in less malignant GIST (73%, 24 of 33). Malignant GIST samples groups showed mutations, LOH, and aberrant expression dominantly in driver genes associated with signaling pathways of PI3K (PIK3CA, AKT1, and PTEN) and the cell cycle (RB1, CDK4, and CDKN1B). Additionally, we identified potential PI3K-related genes, the expression of which was upregulated (SNAI1 and TPX2) or downregulated (BANK1) in malignant GIST. Based on our observations, we propose that inhibition of PI3K pathway signals might potentially be an effective therapeutic strategy against malignant progression of GIST.

comprehensive cancer panel

copy number alteration

gene expression profiling

gastrointestinal stromal tumor

intestinal cells of Cajal

Integrative Genomics Viewer

short insertion and deletion

nonsense‐mediated decay

oncogene

single nucleotide variants

tyrosine kinase inhibitor

tumor mutation burden

tumor suppressor genes

Torrent Variant Caller

whole exome sequencing

INTRODUCTION

Gastrointestinal stromal tumors account for less than 1% of all gastrointestinal (GI) tumors and for approximately 5% of all sarcomas.1 However, they are the most common (80%) mesenchymal tumors of the GI tract. Worldwide annual incidences are largely consistent, ranging between 11 and 19.6 per million population.2 GIST originate from ICC or ICC‐like stem cell precursors.2 Age of onset ranges between teenage years to the 90s, with a peak onset age of 60 years.3 The most common organ affected is the stomach (50%‐60%) followed by small intestine (30%‐35%), rectum (5%), and esophagus (<1%).3 GIST often recur locally within the abdomen and/or are metastasized to the liver.2 Prognosis for GIST has been assessed by risk stratification schemes, including the modified NIH classification.4 This criterion, which considers tumor size, mitotic count, tumor location, and tumor rupture, is useful in identifying patients who might benefit from adjuvant therapy.5

Most GIST contain gain‐of‐function mutations in one of the two receptor tyrosine kinase genes, KIT (75%‐80%) and PDGFRA (5%‐10%), resulting in conformational changes of the respective proteins to constitutively activate downstream signaling pathways, including RAS/RAF/MAPK and PI3K/AKT/mTOR.2 The remaining 10%‐15% of GIST without KIT or PDGFRA mutations show different clinical and pathological features from the mutation‐carrying GIST, and these so‐called wild‐type GIST include neurofibromatosis‐type 1NF1, Carney triad, and Carney‐Stratakis syndrome.3 Discovery of TKI, including imatinib6 and sunitinib,7 targeting these oncogenic mutants, made a significant clinical impact on the drug treatment for GIST patients. Different types of mutations, including point mutations, deletions, and insertions can be found in the different exons of KIT and PDGFRA, and these mutations show type‐ and location‐specific relationships with risk stratification, clinical manifestations, and drug response. For example, GIST patients carrying KIT exon 11 deletion mutations show poor prognosis,8, 9 but are sensitive to imatinib.10 In KIT, oncogenic mutations found in exons other than exons 9 and 11 show primary resistance to imatinib, and some mutations in these regions arose from first‐line treatment as secondary resistant mutations.2, 11 Thus, second‐ or third‐line drug treatment is carried out for patients showing resistance to TKI.

Anderson et al and others have suggested that KIT/PDGFRA mutations are very common events in the early stage of GIST development, and are not sufficient for GIST progression, and that other genetic changes are required for clinical manifestation.12, 13, 14 Several chromosomal changes, including deletions in chromosome arms 1p, 13q, 14q, 15q, and 22q,14, 15, 16, 17 and gains in chromosomes 4 and 5,16, 17 have been associated with malignant progression of GIST. However, there are no consensus genetic alterations with mutations, amplification, or deletion for GIST development.13 Thus, we aimed to identify driver gene alterations and the subsequent signaling pathways that drive the progression of GIST in order to develop effective therapies, particularly for TKI‐resistant patients.

MATERIALS AND METHODS

Clinical samples

Tumor tissue samples dissected from surgical specimens, along with whole blood samples were obtained from Shizuoka Cancer Center under protocols approved by the institutional review board at Shizuoka Cancer Center (authorization number: 25‐33).18 Written informed consent was obtained from all patients enrolled in the study. All experiments using clinical samples were carried out in accordance with the approved guidelines.

Sequencing analysis

Isolation and characterization of genomic DNA for WES and CCP has been described previously.19, 20 Exome library for WES was constructed using an Ion AmpliSeq Exome RDY Kit (Thermo Fisher Scientific). Briefly, 100 ng DNA was used in the target amplification under the following conditions: 99°C for 2 minutes, followed by 10 cycles at 99°C for 15 seconds and 60°C for 16 minutes, and a final hold at 10°C. Amplicons were ligated with Ion Torrent Proton adapters (Thermo Fisher Scientific) at 22°C for 30 minutes followed by 72°C for 10 minutes, and the library was purified with Agencourt AMPure XP beads (Beckman Coulter). This exome library supplied 293 903 amplicons that cover 57.7 Mb of the human genome, comprising 34.8 Mb exons of 18 835 genes registered in RefSeq. The constructed library was quantified using quantitative PCR, and DNA was sequenced using a semiconductor DNA sequencer (Ion Proton Sequencer; Thermo Fisher Scientific).

For CCP, the targeted DNA library comprising 1.6 Mb exon regions and splice sites of 409 genes was constructed using an Ion AmpliSeq Comprehensive Cancer Panel Kit (Thermo Fisher Scientific) with the Ion Library Equalizer Kit (Thermo Fisher Scientific). Briefly, 40 ng DNA was used for multiplex PCR amplification with four separate primer pools. The amplicons were treated with FuPa reagent (Thermo Fisher Scientific) and ligated to a uniquely barcoded adapter. After purification using Agencourt AMPure XP beads, the constructed library from each primer pool was quantified using the Ion Library Quantification Kit (Thermo Fisher Scientific) and pooled together. Template preparation was carried out using the Ion Chef System (Thermo Fisher Scientific) and Ion PI Hi‐Q Chef Kit (Thermo Fisher Scientific). DNA was sequenced using a semiconductor DNA sequencer (Ion Proton Sequencer).

Binary raw data derived from the semiconductor DNA sequencer were converted using Torrent Suite software (Thermo Fisher Scientific) into sequence reads that were mapped to the reference genome (UCSC hg19). At this step, sequence data derived from tumor and blood samples were saved as BAM files. Then, somatic mutation calling was applied separately for the WES and CCP datasets. For WES, two BAM files were uploaded to the Ion Reporter system and analyzed concurrently with AmpliSeq exome tumor‐normal pair workflow (ver. 4.4; Thermo Fisher Scientific) with a Custom Hotspot file that specifies somatic and pathogenic mutations registered in COSMIC and ClinVar, respectively. The list of identified mutations was processed by in‐house scripts to remove false‐positive calls, including sequencer‐derived errors. Mutations fulfilling at least one of the following criteria were discarded as false‐positive: (i) quality score <60; (ii) depth of coverage <20; (iii) variant read observed in one strand only; (iv) clipped sequence length <100 (avg_clipped_length <100); (v) variant located on either sequence end (avg_pos_as_fraction <0.05); and (vi) mutation matches one on an in‐house false‐positive list. Parameters specified in criteria (iv) and (v) were calculated by bam‐readcount with option "−q 1" (ver. 0.8.0) (https://github.com/genome/bam-readcount). For CCP, tumor CCP sequence reads were compared with blood WES reads to identify somatic mutations. First, variant calling for tumor samples was carried out using TVC (ver. 4.4; Thermo Fisher Scientific), and then blood data were analyzed by TVC with custom hotspot, which specifies all detected mutations in a tumor sample to determine whether these mutations were observed in the blood sample. Because a different variant caller and higher depth were used compared to those used for WES, the following criteria were used to identify unreliable mutations: (i) quality score <50; (ii) depth of coverage <20; or (iii) mutation matches one on an in‐house false‐positive list.

Effects of mutations were predicted using SnpEff.21 Nonsynonymous mutations, including substitutions, insertions or deletions at coding regions and splice‐sites were visually confirmed on the IGV,22 and subsequently validated by Sanger sequencing. PCR and Sanger sequencing were carried out as described previously,23 and primer sequences are shown in Table S1. Based on the number of WES reads, CNA was calculated according to the saasCNV method.24 We defined gain and loss in the case of copy number as ≥2.5 and <1.5, respectively. WES data were applied to estimate tumor purity using an in silico method.25 Data of SNV and INDEL from WES and CCP are listed in Data S1. CNA data from WES are listed in Data S2. For fusion gene analysis, total RNA was used as a template to prepare cDNA, and subjected to the Ion Proton System for detecting fusion transcripts from a panel of 491 fusion genes, as previously reported.26

Gene expression analysis

Total RNA was isolated and subjected to microarray analysis as described previously18 using SurePrint G3 Human Gene Expression 8 × 60K v2 Microarray (Agilent Technologies). RNA samples with RNA integrity number ≥5.9 were used for microarray analysis. Microarray analysis was carried out in accordance with the MIAME guidelines.27 Data analysis was carried out using GeneSpring GX (Agilent Technologies), Subio platform (Subio), and Microsoft Excel. Probes to be analyzed were selected according to the reference genome sequence, hg19, obtained from the UCSC Genome Browser.28 Raw signal intensity values were log‐transformed and normalized to the 75th percentile. Microarray data for mRNA expression are available through the NCBI database under accession http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE136755.

Statistical analysis

We used Fisher’s exact test for comparison of two datasets. Welch’s t test was carried out in the assumed normal distribution. P‐values <0.05 were considered statistically significant. We used the Z‐score, which indicates the number of standard deviations away from the mean of expression, to predict significant changes in gene expression.

RESULTS

Patient and tumor sample characteristics

A total of 65 surgically resected GIST tumor samples were obtained from 64 patients, which included six metastatic and 59 primary tumors (Table S2). One patient provided primary and matched metastatic tumor samples resected at a different time period. According to the risk criteria,4 primary GIST tumor samples were divided into four groups with high (17 cases), intermediate (9 cases), low (22 cases) and very low (11 cases) risk of progression. Nine patients were treated with imatinib or sunitinib before surgery. With a combination of WES, targeted sequencing of 409 cancer‐associated genes using the Ion AmpliSeq CCP, and Sanger sequencing (Tables S2 and S3), we identified oncogenic KIT and PDGFRA mutations in 57 (88%) and six (9%) of the 65 GIST samples, respectively. Among the KIT driver mutations, short deletions in exon 11 were observed more frequently (P = 8.5 × 10−3) in the metastatic/high‐risk groups (15 of 23 tumors) than in the other risk groups (15 of 42 tumors), as reported previously.8, 9 In addition to the drivers, secondary KIT mutations,2, 11 including V654A in exon 13, T670I in exon 14, and N822Y in exon 17, were identified in three samples. CNA data by WES showed copy number gain of KIT in eight samples, of which a sample was from KIT/PDGFRA wild type (Figure S1). Three samples showed copy neutral loss of heterozygosity (cnLOH) in KIT. GEP showed that KIT mRNA expression levels were lower in samples harboring PDGFRA mutations than in those with KIT mutations (P = 1.8 × 10−2), whereas PDGFRA mRNA levels were lower in samples harboring KIT mutations than in samples with PDGFRA mutations (P = 4.0 × 10−4). These observations were consistent with a previous report.17 No correlation (P = 2.8 × 10−1) was observed in KIT expression levels between samples with (copy number ≥2.5) or without (copy number = 2) KIT copy number gain. Even excluding samples harboring PDGFRA mutations, this conclusion was unchanged (P = 6.0 × 10−1).

Summary of WES and CCP

Mean depths of coverage for blood and tumor tissue sequences in WES were 131 (range, 94‐164) and 126 (range, 93‐126), respectively. In CCP, mean depth of tumor tissue sequences was 1183 (range, 851‐1438). Total numbers of nonsynonymous mutations in WES and CCP were 1084 (mean, 16.7) and 213 (mean, 3.3), respectively. Among the 1084 mutations from WES, 113 mutations were derived from 409 CCP genes, and 85 out of 113 mutations (75.2%) were also detected in CCP. In contrast, 39.9% of mutations detected in CCP (85 out of 213) matched mutations detected by WES. Size of TMB, defined as the total number of synonymous/nonsynonymous mutations per megabase obtained from WES, ranged from 0.06 to 1.75 with a median of 0.77 (Figure S1), indicating no TMB‐high (TMB ≥20) samples in this set of GIST samples.19

In addition to the GIST‐initiating mutations of KIT and PDGFRA, we used WES and CCP data to investigate other genetic alterations against GIST progression, including SNV and INDEL of the 707 driver genes (Table S4 and Data S1) that were selected referring to several publications29, 30, 31 and databases (COSMIC,32 https://cancer.sanger.ac.uk/cosmic; Oncomine Comprehensive Assay v3, OCAv3, https://www.thermofisher.com/jp/ja/home/clinical/preclinical-companion-diagnostic-development/oncomine-oncology/oncomine-cancer-research-panel-workflow.html). Of these, 320 and 301 genes were defined as OG and TSG, respectively, and the remaining 86 genes had both characteristics. Among the 1084 mutations detected by WES, 128 mutations were derived from the 707 driver genes, and their presence was further confirmed by IGV analysis. We carried out CCP to enhance detection rate of driver gene mutations, as CCP gives more depth of coverage. Among the 213 mutations detected by CCP, 158 mutations were derived from the driver genes, 76 mutations of which were also detected by WES. For the 82 mutations uniquely detected by CCP, our validation by IGV confirmed the presence of only four mutations, including FANCD2 (no. 8), KIT (nos 12 and 14), and KMT2A (no. 64) (Data S1B). Among the rest of the mutations, 42 out of 78 mutations (53.8%) were recurrently found mainly in the regions containing homopolymers and repetitive sequences, and primer‐ends, suggesting false positivity. Excluding KIT and PDGFRA mutations, 68 driver gene mutations, consisting of 21 mutations defined as OG, 39 mutations defined as TSG, and eight mutations defined as OG/TSG, were used for further analysis.

Oncogene mutations

In the case of the 406 OG, 29 mutations in 23 genes, consisting of 21 and eight mutations from 320 OG and 86 OG/TSG, respectively, were observed in 24 of the 65 samples (Figure 1). Among them, PIK3CA mutations, including G106R (observed in a metastatic sample, no. 1) and R88Q (observed in paired primary and metastatic samples, nos 6.1 and 6.2), were the only mutations that have been identified as activating mutations.33, 34 Of the remaining 26 mutations, MUC4, HDAC1, and MUC16 mutations were found in multiple samples, but the mutation patterns in each gene were different. Additionally, we also carried out fusion gene analysis and identified the oncogenic COL1A1‐PDGFB fusion transcript35 in an intermediate‐risk sample (no. 33) (Figure 1). Taken together, the observations in the two patients with either metastasis or high‐risk GIST suggests that PIK3CA mutations are possible driver alterations for malignant progression of GIST.

Figure 1

Driver gene mutations in 65 gastrointestinal stromal tumor (GIST) samples. Identified mutations of 707 driver genes, including oncogenes (OG) and tumor suppressor genes (TSG), and a fusion gene (F) are indicated. TSG mutations accompanied by loss of heterozygosity (LOH) are shown in red squares with copy number (CN) <1.5, and black squares with copy neutral LOH (cnLOH)

Tumor suppressor gene dysfunction by LOH

Oncogenes initiate carcinogenesis when mutations dominantly occur within a single copy of the gene, whereas TSG are required to follow Knudson’s ‘two‐hit hypothesis’ to recessively inactivate the gene.36 In addition to the classical TSG inactivation, even partial TSG inactivation, as a result of haploinsufficiency, or dominant‐negative TSG mutations could contribute to tumorigenesis.36, 37, 38 TSG inactivation is also due to epigenetic mechanisms of gene silencing (eg, hypermethylation of CpG islands located in the promoter region).39 Thus, we investigated TSG dysfunction considering the following categories: (i) LOH; (ii) deleterious mutations with expression reduction; (iii) haploinsufficiency/dominant‐negative mutations; (iv) copy number loss accompanied by expression reduction; and (v) expression reduction without copy number loss.

First, we searched LOH‐related TSG mutations among the 47 mutations, consisting of 39 and eight mutations from 301 TSG and 86 OG/TSG, respectively, using a combination of mutation and CNA data, and identified 15 mutations in 12 genes accompanied by copy number loss (copy number <1.5) (Figure 1, Data S2). Despite a lack of information on carcinogenesis, nonsense and frameshift mutations are predicted to be deleterious as a result of disruption of protein structure. According to the American College of Medical Genetics and Genomics (ACMG) standards and guidelines,40 certain types of variant, including nonsense and frameshift, disrupt gene function by leading to a complete absence of the gene product by lack of transcription or NMD of an altered transcript. However, ACMG states that we must carefully interpret truncating variants to consider pathogenesis if the predicted stop codon occurs in the last exon or in the last 50 base pairs of the penultimate exon, in which NMD would not be predicted. Considering the locations, all the four nonsense mutations found in RB1, MGA, and CBL (Table 1) were predicted to be nonfunctional by NMD. However, it is unclear whether the missense mutations play tumor‐suppressive roles. In addition, the effect of copy number loss on mRNA expression was assessed by GEP data. As summarized in Table 1, all of the seven deleterious mutations found in RB1, PTEN, TSC1, MGA, and CBL showed reduced expression levels (Z‐score <0), which were designated as ‘very likely’ for TSG inactivation. One of the two RB1 mutations was copy neutral LOH. Eight missense mutations from seven genes, including NOTCH3, PMS2, NF2, TNFRSF14, TCF7L1, MAX, and NF1, were designated as ‘possible’. Inactivation of NF1 by LOH or other factors could have led to GIST in neurofibromatosis41 in patient no. 31, in which no KIT and PDGFRA mutations were identified. LOH mutations were observed in five metastatic, four high‐risk, one intermediate, and two low‐risk tumor samples. This indicated a higher rate of LOH presence in the metastatic/high‐risk GIST samples than in the other samples (P = 2.6 × 10−3). In comparison with GIST derived from the stomach and other tissues, metastatic/high‐risk GIST showed a higher proportion in tissues (65%, 13 of 20) other than the stomach (22%, 10 of 45; P = 1.6 × 10−3). Accordingly, the presence of LOH was higher in GIST from other tissues (35%, 7 of 20) than in GIST from the stomach (11%, 5 of 45; P = 3.6 × 10−2).

Table 1

Loss of heterozygosity (LOH)‐driven tumor suppressor gene (TSG) inactivation

Case no.	Meta/Primary	Risk	Gene	Cytoband	Mutation						Copy no.	Expression			Purity	Prediction of TSG inactivationb
					Type	Allele	Site (coding DNA)	Site (protein)	Exon (exon position/no. of total exons)	Total length of protein (no. of amino acid residues)		Signal	Z‐score	Reduction (Z‐score <0)a
1	Meta	NA	RB1	13q14.2	NS	hetero	c.649C > T	p.Q217*	7/27	928	1.44	1.12	−2.88	Yes	0.58	Very likely
1	Meta	NA	NOTCH3	19p13.12	MS	hetero	c.3658C > T	p.R1220W	22/33	2321	1.03	2.01	0.08	No	0.58	Possible
2	Meta	NA	PMS2	7p22.1	MS	hetero	c.958T > C	p.F320L	9/15	862	1.42	‐0.52	−1.02	Yes	0.22	Possible
3	Meta	NA	PTEN	10q23.31	F	hetero	c.223_224delCA	p.H75fs	4/9	403	1.34	0.32	−2.58	Yes	0.51	Very likely
3	Meta	NA	NF2	22q12.2	MS	hetero	c.1414C > A	p.L472I	13/16	595	1.22	3.41	−0.74	Yes	0.51	Possible
5	Meta	NA	TNFRSF14	1p36.32	MS	hetero	c.17A > G	p.D6G	1/8	283	1.14	4.78	0.25	No	0.53	Possible
6.1	Meta	NA	TSC1	9q34.13	F	hetero	c.1360_1373delGTCACTCTAAGTGA	p.V454fs	14/23	1164	1.19	0.01	−3.67	Yes	0.29	Very likely
6.2	P	H	TSC1	9q34.13	F	hetero	c.1360_1373delGTCACTCTAAGTGA	p.V454fs	14/23	1164	0.86	0.47	−2.89	Yes	0.49	Very likely
6.2	P	H	TCF7L1	2p11.2	MS	hetero	c.887C > T	p.A296V	8/12	588	1.45	0.84	−0.57	Yes	0.49	Possible
8	P	H	MAX	14q23.3	MS	homo	c.108G > T	p.R36S	3/5	160	1.07	3.56	0.01	No	0.65	Possible
11	P	H	NOTCH3	19p13.12	MS	hetero	c.5707A > G	p.M1903V	31/33	2321	1.47	3.05	0.86	No	0.29	Possible
18	P	H	RB1	13q14.2	NS	homo	c.1324G > T	p.G442*	13/27	928	2 (cnLOH)	1.85	−1.64	Yes	0.56	Very likely
31	P	Int	NF1	17q11.2	MS	hetero	c.237A > T	p.L79F	3/58	2839	1.4	0.26	−2.14	Yes	0.19	Possible
50	P	L	MGA	15q15.1	NS	hetero	c.5654C > A	p.S1885*	17/24	3065	0.82	1.38	−1.17	Yes	0.34	Very likely
52	P	L	CBL	11q23.3	NS	hetero	c.891T > A	p.C297*	6/16	906	1.46	−0.67	−1.24	Yes	0.46	Very likely

Abbreviations: cnLOH, copy neutral LOH; F, frameshift; H, high; hetero, heterozygous; homo, homozygous; Int, intermediate; L, low; Meta, metastasis; MS, missense; NA, not available; NS, nonsense; P, primary.

Expression level was deduced to be ‘reduction’ when Z‐score was <0.

TSG inactivation was predicted to be ‘Very likely’ with nonsense or frameshift mutations, and ‘Possible’ with missense mutation.

Translation termination due to the presence of a stop codon.

Tumor suppressor gene dysfunction by other alterations

Among the five TSG harboring deleterious mutations without copy number loss (copy number ≥1.5), a frameshift mutation in RBM10 showed reduced mRNA expression level (Table S5). This result indicates that mutation of RBM10, an RNA‐binding protein and splicing regulator, has a tumor‐suppressive role. Genes with the remainder of the four deleterious mutations, along with 27 genes harboring the missense and in‐frame deletion mutations, can disrupt tumor‐suppressive function of protein by haploinsufficiency or dominant‐negative effect. However, according to the ACMG standards and guidelines42 described previously, the nonsense USP44 mutation was predicted to maintain its function as a result of its location on the penultimate exon. Thus, considering the types of mutations and locations, four deleterious mutations in RBM10, CHD2, MGA, and MAML3 were designated as ‘very likely’ for TSG inactivation. However, these mutations may not directly contribute to tumor growth and progression of GIST, commonly referred to as passenger mutations.

Next, we integrated TSG inactivation by genes carrying copy number loss with expression reduction. This corresponds to TSG inactivation as a result of copy number loss in one of the alleles and expression reduction, due to epigenetic changes (eg, promoter methylation) in the other allele. TSG showing copy number loss (copy number <1.5) and expression reduction (Z‐score ≤−2.5) were observed in 55 genes with a total of 66 alterations (Figure S2). Four of these cases were identified as LOH‐related TSG (Table 1). In GIST, deletions in chromosome arms 1p, 13q, 14q, 15q, and 22q are frequently observed,14, 15, 16, 17 and were also observed in our samples (Figure S3). Particularly, our CNA data showed preferential copy number loss at 13q in the metastatic/high‐risk GIST samples. Among the genes, expression reduction of RB1 was observed in multiple cases, including in two metastatic and one high‐risk tumors (Figure S2). This observation was consistent with the previous report that deletion at the RB1 locus frequently occurred in recurrent or metastatic GIST.43

Last, we integrated TSG inactivation by genes showing expression reduction without copy number loss. This corresponds to TSG inactivation as a result of epigenetic modifications in both alleles. TSG showing expression reduction (Z‐score ≤−2.5) without copy number loss (copy number ≥1.5) were observed in 112 genes with a total of 150 alterations (Figure S4). Seven genes, including ASXL2, ARID1B, EXT1, CREB3L1, FANCF, NPRL3, and SMARCE1, were downregulated in multiple samples of either metastatic or high‐risk tumors. TGFBR2, RUNX1TX, CDKN1B, and CDH1 showed expression reduction in multiple samples from either metastatic or high‐risk tumors along with low‐risk tumor samples.

Oncogene amplification

Oncogenes are activated by mutations, amplification, and chromosomal rearrangements, causing either an alteration in oncogene structure or an increase in or deregulation of its expression.44 In addition to mutation and fusion gene identification, amplification of OG were investigated. Our previous report showed that amplified genes are not always upregulated.18 Thus, we integrated GEP data with CNA data to assess OG amplification. Among the 406 OG, 98 genes showed copy number gain (copy number ≥2.5) accompanied by expression enhancement (Z‐score ≥1.5; Figure S5). Copy number gains of chromosomes 4 and 5 are frequently observed in GIST,16, 17 which were also shown in our samples (Figure S3). Additionally, our data showed gains of chromosome 20q in the metastatic/high‐risk GIST samples, in which expression levels of PLCG1, ZNF217, and GNAS were enhanced in the metastatic or high‐risk tumor samples. Other OG, which were identified in multiple tumor samples from independent metastatic/high‐risk patients, included SKP2, HOXA9, EZH2, CDK4, HMGA2, and FZD10.

Integration of activated OG and inactivated TSG

We summarized the activated OG and inactivated TSG in each tumor sample (Figure 2). Based on the types of alterations in OG and TSG, their effects on driver potential were classified into two types, including higher potential (‘very likely’) and lower potential (‘possible’), as follows. The ‘very likely’ alterations were: (i) OG (PIK3CA) mutations known as activating mutations (Figure 1); (ii) OG amplified with increased expression (Figure S5); (iii) nonsense or frameshift TSG mutations with LOH (Table 1); (iv) frameshift TSG (RBM10) mutation with reduced expression (Table S5); (v) frameshift TSG mutations without reduced expression (Table S5); (vi) TSG with copy number loss and reduced expression (Figure S2); (vii) TSG with reduced expression (Figure S4); (viii) oncogenic fusion transcripts (Figure 1). The ‘possible’ alterations were: (i) OG mutations with unknown function (Figure 1); (ii) mutations found in genes defined as OG/TSG (Figure 1); (iii) missense TSG mutations with LOH (Table 1); (iv) missense TSG mutations (Figure 1). ‘Very likely’ driver gene alterations with showed higher incidence in the metastatic/high‐risk GIST samples (96%, 22 of 23), including metastatic (100%, 6 of 6) and high‐risk (94%, 16 of 17) tumors, than in the other risk GIST samples (73%, 24 of 33) of low‐ and very low‐risk tumors (P = 3.6 × 10−2). In comparison with GIST derived from the stomach and other tissues, proportion of metastatic/high‐risk GIST was higher in tissues (65%, 13 of 20) other than the stomach (28%, 10 of 36; P = 1.1 × 10−2). Accordingly, ‘very likely’ driver gene alterations were shown to be greater in GIST from other tissues (100%, 10 of 10) than from stomach GIST (72%, 26 of 36; P = 9.5 × 10−3).

Figure 2

Summary of genetic alterations and expression modulations of driver genes. Alterations of driver potential were assigned and predicted by two different potential levels, including ‘Very likely’ and ‘Possible’ levels, which are shown in red and blue, respectively. LOH, loss of heterozygosity; TSG, tumor suppressor gene

The only metastatic/high‐risk GIST sample in which no drivers were identified was a high‐risk sample (no. 9). This sample remained without driver alterations even considering lower potential (‘possible’). Besides the 707 driver genes, we investigated other genes with driver potential in sample no. 9 by comparing GEP data between the metastatic/high‐risk GIST and the other risk GIST (low‐ and very low‐risk) samples (Figure 3A). To identify potential OG, 14 genes showing copy number gain (copy number ≥2.5) and expression enhancement (Z‐score ≥1.0) in sample no. 9 were extracted from the genes significantly upregulated (fold change ≥2, P < .05) in the metastatic/high‐risk GIST samples (Figure 3B). The 14 upregulated genes were located on chromosomes 4q or 20q, where copy number was gained in sample no. 9 (Figure S3). As a potential TSG, only BANK1 was extracted as a gene showing expression reduction in sample no. 9 among the genes significantly downregulated (fold change ≤0.5, P < .05) in the metastatic/high‐risk GIST samples (Figure 3B). For all 15 potential driver genes, no significant difference in expression level was observed between the metastatic and high‐risk tumors (P > .05).

Figure 3

Identification of oncogenic or tumor‐suppressive candidate genes. A, Volcano plot showing differentially expressed genes between two groups with different degrees of malignant progression, consisting of metastasis (M) and high‐risk (H), and low (L) and very low (vL) risk. All circles represent 25 434 microarray probes corresponding to mRNA; closed and open circles for probes showing copy number (CN) ≥2.5 in sample no. 9 with expression levels of Z‐score ≥1.0 and Z‐score <1.0, respectively. Closed triangle indicates probe showing expression reduction with Z‐score ≤−2.5 in sample no. 9. Gray circles indicate other probes. Log2 fold change in the H/M versus L/vL is represented on the x‐axis. The y‐axis shows the log10 of the Q value. A Q value of 0.05 and a fold change of 2 are indicated by horizontal and vertical dotted lines, respectively. B, Z‐score expression levels of 15 genes in 65 GIST tumor samples. Probes corresponding to these genes were selected as copy number ≥2.5 and Z‐score ≥1.0 in sample no. 9 (closed circles in Figure 3A) or as Z‐score ≤−2.5 in sample no. 9 (closed triangle in Figure 3A) among the probes showing differential expression in the M/H group. C, Expression and copy number levels of 15 genes. In each panel, the open red circle represents sample no.9 and closed red circles represent other samples that belong to the M/H group. Gray squares and blue circles represent the intermediate‐risk (Int) and L/vL groups, respectively.

Signaling pathway activation

Understanding of the genes along with pathways altered in tumor samples is essential to identify potential therapeutic options and vulnerabilities. According to the published information,30, 31 we determined a pathway related to each driver gene (Table S4). Then, we assigned driver genes with alterations as ‘very likely’ driver potential (Figure 2) to 27 pathways. Genetic alterations in the pathways of HIPPO, WNT, PI3K, NOTCH, cell cycle, and RAS‐MAPK were observed in 10 or more samples (Figure 4A). By comparing the alteration frequency, PI3K and cell cycle pathways were significantly (P = 2.1 × 10−2) altered in the metastatic/high‐risk GIST samples (metastatic and high‐risk tumors) than in the other risk GIST samples (low‐ and very low‐risk tumors (Figure 4B). Despite statistical nonsignificance, alteration in the chromatin SWI/SNF complex pathway was found in five samples, all of which were derived from the metastatic/high‐risk tumor samples. In all the samples, the presence of deletion mutations in KIT exon 11 was not significantly related to alterations in PI3K and cell cycle pathways (P = 3.8 × 10−1, P = 7.4 × 10−1, respectively). Also, in the metastatic/high‐risk tumor samples, no correlations between them were found (P = 4.0 × 10−1, P = 1.8 × 10−1, respectively). The PI3K pathway has been reported to be activated as a result of imatinib secondary resistance.42, 45, 46 In the present study, we observed PI3K pathway alteration even in tumor samples derived from patients not treated with neoadjuvant imatinib. In the metastatic/high‐risk tumor samples, the alterations in PI3K and cell cycle pathways were not significantly related to neoadjuvant imatinib (P = 4.1 × 10−1, P = 4.2 × 10−1, respectively).

Figure 4

Curated pathways. A, Genetic alterations and expression modulations found in the 27 pathways for each gastrointestinal stromal tumor (GIST) sample. On the right, frequency of pathway involvement was compared between two groups with different degrees of malignant progression, consisting of metastasis (M) and high‐risk (H), and low (L) and very low (vL) risk. B, Pathway members and interactions in the PI3K and cell cycle pathways. LOH, loss of heterozygosity; TSG, tumor suppressor gene

DISCUSSION

KIT or PDGFRA mutations activate most GIST (~90%). However, their prognosis and development are independent of the types of mutation. Thus, these mutations are believed to be early events in GIST development, which involves additional genetic alterations for malignant progression.12, 13, 14 Saponara et al reported that metastatic GIST showed frequent copy number loss at regions located on TSG, although no shared oncogenic mutations were observed except in KIT.13 In the present study, despite the fact that there were no common mutations in driver genes of GIST with the maximum frequency of 5%, we showed that PI3K and cell cycle signaling pathways were involved in GIST malignant progression.

KIT and PDGFRA oncogenic mutations activate downstream signaling pathways, including RAS/MAPK, PI3K/AKT, and STAT3.2 Saponara et al reported frequent inactivation of tumor‐suppressive PTEN, which is involved in the PI3K pathway, as a result of LOH or copy number loss with expression reduction in metastatic GIST.13 In addition to LOH of PTEN, we dominantly identified genetic alterations of PI3K pathway genes in malignant GIST samples with metastatic or high‐risk group, including oncogenic activation by PIK3CA mutations and copy number gain with expression enhancement in AKT1 and RICTOR, and tumor‐suppressive inactivation by LOH or expression reduction of TSC1/2, and NPRL3 (Figure S6). The PI3K pathway has been reported to be activated as a result of imatinib secondary resistance, including GIST.42, 45, 46 In the present study, we observed PI3K pathway alteration even in tumor samples derived from patients who were not treated with neoadjuvant imatinib, suggesting that PI3K activation was irrelevant to imatinib treatment in metastatic and high‐risk GIST. Also, no correlation was observed between KIT exon 11 deletions and PI3K activation.

Our observations indicated that the cell cycle was another pathway showing genetic alterations significantly observed in metastatic and high‐risk GIST. Reduced mRNA expression of tumor‐suppressive CDKN2A gene leads to upregulation of E2F1, increasing cell proliferation to drive poor prognosis for GIST.47 Saponara et al reported that copy number losses of CDKN2A and CDKN2B were the most frequent in metastatic GIST.13 We also observed significant losses of CDKN2A and CDKN2B in metastatic and high‐risk GIST compared with the other risk GIST (P = 6.4 × 10−3). However, as a result of insufficient decreased mRNA levels, these genes were not identified as altered TSG. Inactivating mutations of other cell cycle pathway genes, including RB1 and TP53, occur in high‐risk GIST, showing activation of the cell cycle pathway.48, 49 We also observed inactivation of RB1, but not of TP53. CDKN1B, another TSG, was downregulated in two samples from each metastasis and high‐risk GIST group. As CDKN1B is suppressed by AKT,2 the cell cycle pathway can cooperate with the PI3K pathway.

We identified 15 potential driver genes from a GIST sample (no. 9) without alterations in the 707 known driver genes. Analysis of associated pathways showed that SNAI1, TPX2, and BANK1 were involved in the PI3K pathway. SNAI1, a zinc‐finger transcription factor, promoted cell migration through downregulation of SERPINB5 (Maspin) and subsequent activation of PI3K/AKT‐dependent Rac1 during prostate cancer progression.50 TPX2, targeting protein for Xklp2, which is a microtubule‐associated protein, suppressed proliferation through repressing the PI3K/AKT signaling pathway in breast cancer cells.51 Downregulation of BANK1 promoted CD40‐dependent AKT activation in B cells.52, 53 Collectively, in sample no. 9, SNAI1 and TPX2 were activated by amplification‐dependent expression enhancement, and BANK1 was inactivated by expression reduction. Among other metastatic/high‐risk GIST, 12 and nine samples were upregulated (Z‐score ≥1.0) without copy number gain for SNAI1 and TPX2, respectively (Figure 5A). Two more metastatic/high‐risk GIST showed reduced expression (Z‐score ≤−2.5) of BANK1. Based on activated AKT (Figure 5B), insertion of these results into the alteration data identified that 20 (87%) of the 23 metastatic/high‐risk samples were involved in the PI3K pathway, suggesting that its activation drives malignant progression of GIST. In fact, the PI3K pathway has been implicated in metastasis for various types of tumors (review in ref54).

Figure 5

Genetic alterations and expression modulations of PI3K pathway‐related genes. A, Expression levels of SNAI1, TPX2, and BANK1. Expression data were added to the data on the PI3K pathway shown in Figure 4A. B, Overview of the PI3K pathway with involvement of SNAI1, TPX2, and BANK1. Interactions between genes were deduced based on published articles.50, 51, 53

In summary, using multi‐omics analysis, we identified driver gene alterations and subsequent signaling pathways. Our data indicate that the PI3K and cell cycle pathways play important roles in GIST malignant progression, which can be of significance for prognosis and treatment of GIST. Particularly, we propose that the development of PI3K inhibitors55 is of potential benefit for patients with metastatic/high‐risk GIST.

DISCLOSURE

Authors declare no conflicts of interest for this article.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.