Identification of High Clonality and Low Clonality Subtypes in GC by Clonality Analysis.

The 78 GC samples in discovery panel and their matching blood cells were subjected to WES. A mean coverage depth of 167X was achieved for GC tumor DNA, and 170X for matched normal blood cells. Eighty-five percent of exons were covered at 20X or greater depth for both tumor and normal samples (Fig. S1 and Dataset S2). We applied MuTect (13) and VarScan2.2.5 (14) to whole-exome data to detect somatic point and indel mutations. A total of 13,866 mutations were detected in this cohort, including 3,421 synonymous, 8,558 missense, 576 nonsense, 241 splice site, 967 frameshift, and 103 nonframeshift indels. The overall mutation frequency was 112.5 mutations (range of 1–1,486 mutations) with a median of 82.5 for nonsilent mutations per tumor (Dataset S3).

To gain insight into GC heterogeneity, we performed SNV clonality analysis using a recently reported algorithm SciClone (15), which estimates the clusters of clones in each sample. This analysis identified cluster numbers ranging from 1 to 8 in our cohort with a bimodal distribution dipping at a cluster number of 4. We thus define those GC cases with cluster numbers larger than 4 as high clonal or HiC (nine cases) (Fig. 1A and Dataset S4) and cases with cluster numbers equal to or fewer than 4 as low clonal or LoC (68 cases) (Fig. 1B and Dataset S4). Association analyses with clinical parameters revealed that the HiC subtype was associated with older disease onset age (Fig. 1C) (Pearson correlation coefficient = 0.27, P = 0.02). There were a total of 11,768 protein-coding subclonal SNVs (92.3% of all detected mutations, with a median of 106 mutations per sample) (Dataset S3). Gene-level mutation analysis revealed that the HiC subtype had a significantly higher frequency of TP53 mutations (Fig. 1D), whereas the LoC subtype was enriched for ARID1A mutations (Fig. 1E), both thought to be the driving events in GC tumorigenesis. Previous studies of leukemia (16) have reported that the higher clonal complexity was associated with poor survival, likely because of chemoresistance. Among our cohort of 78 GC patients, 32 (41.03%) died from GC with a median follow-up of 25.1 months (range from 0.20 to 78.95 months) (Dataset S1). Survival analysis showed that the HiC subtype was significantly associated with shorter survival compared with the LoC subtype (P = 0.02, log-rank test) (Fig. 1F). There were no significant differences in the treatment regimens between the HiC and LoC subtypes. This association remained significant after adjusting for age, sex, stage, Lauren subtype, TP53 mutation, and ARID1A mutation [adjusted hazard ratio (HR), 4.69; 95% confidence interval (CI), 1.62–13.6; P = 0.0043], indicating that the clonal subtype is an independent prognostic factor for GC.

• Download figure
• Open in new tab
• Download powerpoint

Fig. 1.

Identification of HiC and LoC subtypes by clonality analysis and their clinical outcomes. (A and B) SciClone analysis with cluster number 4 as cut-off identified 9 HiC and 68 LoC cases. Each peak infers one subclone by SciClone algorithm. (C) Association between clonality (y axis) and age at diagnosis (x axis). (D) Fractions of TP53 mutation carriers (y axis) in different clonality GC (x axis). (E) Fraction of ARID1A mutation carriers (y axis) in different clonality GC (x axis). (F) Kaplan–Meier survival curves for HiC and LoC cases.

The distinct clinical behavior of the HiC and LoC subtypes suggest they may have fundamental differences in tumorigenesis, as was recently proposed to be manifested by various nucleotide level mutational signatures (17, 18). Detailed mutational signature analysis of these two subtypes showed a significant difference in the percentages of C-to-G transition and C-to-T transversion mutations. HiC GCs had a significantly smaller fraction of C-to-T transversion mutation (45% in HiC vs. 52% in LoC) but greater fraction of C-to-G transition mutation (25% in HiC vs. 9% in LoC, P = 0.002). Specifically, the HiC subtype had a significantly higher fraction of C > G mutation in almost all contexts, except for the G in the 5′-adjacent to C (Dataset S4 A–C).

Although TP53 mutations were commonly found in the HiC group, they were present only in the minor clones (mutation frequency < 15%) in all six positive HiC cases (Dataset S4 D and E). Conceivably, any of the minor clones can become a major clone after a single agent therapy, which has been illustrated in a study in leukemia (16). Thus, in principle, all of the subclones in a tumor should be eliminated for curative therapy. We used a drug–gene interaction database, DGIdb (19), to search for drugable targets for the clones present in the HiC tumors (Dataset S4F). As expected, many gene mutations and the clones harboring them do not have corresponding drugs available. However, for tumors like PGM71 and PGM32, targetable mutations were identified for almost all subclones [e.g., ALK, ABL2, SMAD4 (SMAD family member 4), FANCG, NRG1, KRAS (Kirsten rat sarcoma viral oncogene homolog), and others] (Dataset S4G). In these cases, combination chemo- and targeted therapies might be much more effective than those using single agents. For other HiC cases the targetable mutations were observed in only two or three clones, underscoring the need for drug development strategies that test multiple combinations (Dataset S4H).

Clonal Evolution of GC During Metastasis.

To define the type and extent of intratumoral heterogeneity as GC progresses to metastasis, we performed WGS of two cases (Pt1 and Pt2) where, for each case, we obtained the information by dissection and anatomical mapping of three different areas of primary tumor and two lymph node metastases. Pt1 had a total of 4,082 nonsilent mutations, including 2,809 missense, 143 nonsense, 132 splice site, 958 frameshift, and 40 in-frame indels. Pt2 had a lower mutation rate with 287 nonsilent mutations (253 missense, 14 nonsense, 10 splice site, 3 frameshift, and 7 in-frame). Consistent with the higher mutation rate in Pt1, mutations were found in several DNA-mismatch repair and chromatin remodeling genes, including MSH6, TGFBR2, KDM5A, and MLL4. In Pt2, the SETBP1 gene (p.S944N) was mutated, as it is in myeloid malignancies (20), together with amplification of AKAP9 and CDK6 (cyclin-dependent kinase 6). Intriguingly, mutations of the gene encoding the voltage-dependent calcium channel protein, CACNA1D p.V529G, was found in both Pt1 and Pt2. Copy number alteration analysis showed that a large proportion of the genome was altered in both Pt1 and Pt2, accounting for 516.6 Mb and 706.8 Mb, respectively. Pt2 exhibited consistent genomic amplification of 7q21.2–7q21.3, and Pt1 had amplifications across chr7 and chr8. Pt1 and Pt2 had a comparable number of structural variations with 3,812 and 3,390 breakpoints, respectively.

Because clonal evolution is closely linked to cancer metastasis, we applied phylogenetic tree analysis (SI Methods) to the five samples for each GC case and observed considerable divergence among the three primary tumor samples in each (Fig. 2 and Fig. S2). Interestingly, only one part of the primary tumor was closely related to the two highly similar lymph-node metastases, suggesting they were derived from a common clonal ancestor. The multiple clones uncovered within primary tumor and lymph-node metastases suggests the necessity of combined therapy that targets different genetic changes found within the patient’s combined tumor burden. We therefore used DGIdb (19) to search for potentially drugable targets in different parts of the same tumor and their resultant metastases. Results of the drug–sample interaction analysis suggest that Pt1 might benefit from combined targeting of PI3K, CDK7, and Notch (as well as immunotherapy) (Fig. 2), whereas combined AURKC and CDK6 inhibitors might benefit Pt2 (Fig. S2).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 2.

Clonal evolution in GC metastasis revealed by phylogenetic tree. WGS data from five samplings (three primary and two lymph node metastases) of one GC case were used in the similarity analysis of somatic substitutions. Branch and trunk lengths are proportional to the number of nonsynonymous mutations acquired on the corresponding branch or trunk. Driver mutations and copy number events were acquired by the indicated genes in the branches the arrows indicate. The potential drugs which target each tumor portion are shown as nodes with arrows pointing to the samples.

Validated NRGs-ERBBs Mutations in GC.

Among the total of 13,866 mutations, we also sought to identify clinically actionable mutations in GC. We analyzed the WES data of the 78 GCs with statistical approaches including MutSigCV (Methods) and identified 16 significantly (q value < 0.2) mutated genes. Of these, 13 were mutated in 5% or more of the tumors. The analysis identified previously reported gene mutations in GC: TP53, ARID1A, CDH1 [cadherin 1, type 1, E-cadherin (epithelial)], APC (adenomatous polyposis coli), RHOA, PIK3CA (phosphatidylinositol-4,5-bisphosphate 3-kinase, catalytic subunit alpha), SMAD4, MYC (v-myc avian myelocytomatosis viral oncogene homolog), and KRAS. Further pathway analysis of recurrent mutated genes revealed that mutated genes in our cohort belonged to the TP53 pathway, Wnt pathway, ERBB pathway, and HR pathway, and so forth (Dataset S5). Based on mutation frequency and novelty (Dataset S3), we selected 103 recurrently mutated genes in Wnt (Fig. S3), ERBB (Fig. S4A), and HR (Fig. S5) pathways for targeted sequencing verification in an additional 216 GC samples (Dataset S1) using the Ion-Torrent platform (Methods). Among the 103 genes, 86 (83%) were found have at least one mutation in the validation cohort (Dataset S3 and S6–S10).

The extensive follow-up data associated with our GC cohort allowed us to evaluate the clinical significance of many of the mutations identified in this and previous studies. We first evaluated if there are some genes specifically associated with different anatomical locations in GC. Using 78 WES samples, we identified 10, 3, and 4 genes, which showed moderate enrichment in antrum [i.e., CSMD3 (30.7%), TSHZ3 (15.4%), PCDHA11 (19.2%), ARHGAP28 (15.4%), DST (23.1%), SORCS1 (19.2%), MUC17 (15.4%), PCDH20 (19.2%), USP9X (19.2%), CACNA2D1 (15.4%)], body [i.e., FLG2 (15.4%), GRM1 (15.4%), SYNE2 (15.4%)], and cardia [i.e., CDH1 (23.1%), GFRAL (15.4%), SBNO1 (15.4%), STAB2 (15.4%)], respectively (P < 0.05). Seven (the underlined ones) of those genes remained significant when we included the extra 216 targeted sequencing cases. Intriguingly, none of these gene mutations were observed to significantly associate with same anatomical location in TCGA data, suggesting different genetic background (Russian vs. Chinese) or even diets might contribute to the mutation distribution. In addition to anatomical location, we also interrogated the association between gene mutation and Lauren subtypes of GC. We found that CDH1 somatic mutations (21 of 294 cases, 7.1%) were strongly enriched in the diffuse GC (Fig. S6A) and 8 of these 21 somatic CDH1 mutations were loss-of-function types (e.g., nonsense and frameshift mutations). RHOA mutations were also enriched in diffuse GC, consistent with recent reports (11, 12). Interestingly, Wnt pathway mutations showed a moderate enrichment among intestinal GC. Survival analysis showed that the CDH1-mutated GCs had shorter survival, whereas Wnt pathway mutated GCs had longer survival (log-rank P = 0.04) (Fig. S6B). This difference was not a result of their association with the two Lauren classifications because we did not observe significantly different survival of these latter two groups in our cohort. Interestingly, RHOA-mutated cases were not associated with shorter survivals similar to CDH1-mutated GCs, suggesting RHOA driven and CDH1 driven diffuse GCs are different clinically.

Although ERBB2 has been previously recognized as a therapeutic target for GC because of its gene amplification and mutation (confirmed in our study), our analysis revealed that another member of the ERBB pathway, ERBB4, was significantly mutated in both our discovery and validation cohorts. Mutation of ERBB4 or HER4 has recently been shown to be an oncogenic driver in melanoma (21), although this has not been previously described for GC. Mutated ERBB4 activates both ERBB4 and PI3K-AKT signal transduction and can be inhibited by the dual-tyrosine kinase inhibitor Lapatinib, for which efficacy is presently being evaluated in phase III clinical trials in ERBB4-mutated melanoma (clinicaltrials.gov/ct2/show/NCT01264081). In our discovery cohort, ERBB4 was mutated in seven samples. In the validation cohort, another 13 samples had ERBB4 mutations. Among these mutations, ERBB4 p.R50C, was also observed in melanoma (Fig. 3A). Unlike ERBB2, which has no known ligand, ERBB4 binds to NRGs, particularly NRG1, forming homo- and heterodimers, and leading to the subsequent activation of kinase activity. Previous studies have shown that overexpression of NRG1 leads to the activation of ERBBs and consequent Lapatinib-sensitivity (22). However, no somatic mutation of NRG1 has been reported previously. In our discovery cohort, NRG1 was mutated in five samples. In the validation cohort, NRG1 mutations were observed in another 16 cases. Taking these data together, we find that NRG1 and ERBB4 were mutated in 34 GC samples (11.6%), with 3 samples having mutations in both genes and 31 samples having mutations in either NRG1 or ERBB4 (Fig. 3). Detailed analysis showed that 7 of 21 (33%) ERBB4 mutations occurred in the kinase domain and 4 (20%) in the receptor domains, suggesting these mutations might impact receptor ligand interaction and kinase activity. Interrogation of the 289 GC cases in the TCGA database also revealed ERBB4 mutations present in both the receptor domain (i.e., p.R106H/C) and the kinase domain (i.e., p.V744L and p.774N/G). We also found NRG1 mutations in the TCGA cohort, including the mutations (i.e., p.A221T and p.A225P) in the EGF-like domain (predicted as the functional domain to bind ERBB4), as well as in regions adjacent to the EGF-like domain (i.e., p.E223G, p.R224Q, and p.S226P). These findings suggest that the EGF-like domain is a hotspot mutation for NRG1. Mutations were also present in other NRGs and ERBB family proteins (Dataset S11). Mutual exclusivity analysis was conducted for ERBB and NRG gene members in the combined Tianjin and TCGA cohorts. In 1,000-times permutation, the ERBBs and NRGs’ mutations were mutually exclusive of each other (P = 0.02).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 3.

Frequent NRG1 and ERBB4 gene mutations in GC. (A) Distribution of coding mutations in NRG1 gene. Red: mutations present in Tianjin cohort; green: mutations present in TCGA cohort. (B) Three-dimensional structure for NRG1 protein and mutation locations.

GCs with BRCA2 Mutations Show a Trend for Longer Survival.

Adjuvant chemotherapy as the standard of care for GC has limited success and there is no reliable marker to predict response. Among our 78 WES cohort of GC cases, 5 had BRCA2 mutations (interestingly, there were no BRCA1 mutations). Correspondingly, in the 216 targeted sequencing cohort cases, another 12 had BRCA2 mutations (10 somatic and 2 germ line) (Dataset S3 and S5). A higher frequency (28 of 289 cases) of BRCA2 mutations was also found in the TCGA Russian GC cohort (Fig. 4A). This finding is consistent with epidemiological studies showing that GC is the third major cancer type, joining breast and ovarian cancer, with BRCA2 mutations (23). We hypothesized that the BRCA2 mutation might be clinically important for GC because of accumulating evidence demonstrates that such mutations are associated with improved response to platinum-based chemotherapy and longer survival in ovarian cancer (24). BRCA2 mutant carrier patients are marginally correlated with longer survival in Tianjin cohort (P = 0.07, HR = 0.29, 95% CI: 0.07–1.19). A similar nonsignificant trend was also observed in TCGA cohort (P = 0.23, HR = 0.41, 95% CI: 0.10–1.77). When these two cohorts were pooled to increase the sample size to 569 (and BRCA2 mutation carriers to 45 cases), mutation carriers were significantly correlated with longer survival (log-rank P = 0.03) (Fig. 4B). After adjusting for age, stage, and different population, BRCA2 mutation remained a significant factor predicting longer survival (P = 0.05, HR= 0.37, 95% CI: 0.13–0.96) (Fig. 4C).

• Download figure
• Open in new tab
• Download powerpoint

Fig. 4.

BRCA2 mutation carriers have longer survival. (A) Distribution of coding mutations in BRCA2 gene. Red: mutations in Tianjin cohort; green: mutations in TCGA cohort. The number of total mutations is more than the number of patients is because some patient has multiple BRCA2 mutations. (B) Kaplan–Meier survival curves for BRCA2 wild-type and BRCA2 mutated GCs. (C) Forest plot of multivariable Cox-regression model including clinical and BRCA2 mutations.