Whole-exome sequencing identifies somatic mutations associated with lung cancer metastasis to the brain.

BACKGROUND: Lung cancer is the most aggressive cancer, resulting in one-quarter of all cancer-related deaths, and its metastatic spread accounts for >70% of these deaths, especially metastasis to the brain. Metastasis-associated mutations are important biomarkers for metastasis prediction and outcome improvement.
METHODS: In this study, we applied whole-exome sequencing (WES) to identify potential metastasis-related mutations in 12 paired lung cancer and brain metastasis samples.
RESULTS: We identified 1,702 single nucleotide variants (SNVs) and 6,131 mutation events among 1,220 genes. Furthermore, we identified several lung cancer metastases associated genes (KMT2C, AHNAK2). A mean of 3.1 driver gene mutation events per tumor with the dN/dS (non-synonymous substitution rate/synonymous substitution rate) of 2.13 indicating a significant enrichment for cancer driver gene mutations. Mutation spectrum analysis found lung-brain metastasis samples have a more similar Ti/Tv (transition/transversion) profile with brain cancer in which C to T transitions are more frequent while lung cancer has more C to A transversion. We also found the most important tumor onset and metastasis pathways, such as chronic myeloid leukemia, ErbB signaling pathway, and glioma pathway. Finally, we identified a significant survival associated mutation gene ERF in both The Cancer Genome Atlas (TCGA) (P=0.01) and our dataset (P=0.012).
CONCLUSIONS: In summary, we conducted a pairwise lung-brain metastasis based exome-wide sequencing and identified some novel metastasis-related mutations which provided potential biomarkers for prognosis and targeted therapeutics. 2021 Annals of Translational Medicine. All rights reserved.

Introduction

Lung cancer (LC) is the leading cause of cancer death in both men and women, accounting for one-quarter of all cancer deaths (1). The 5-year survival rate has failed to improve significantly over the last 30 years and remains at a mere 19%, due to recurrence and metastasis. Metastasis accounts for about 90% of cancer-related deaths and is the inevitable outcome of most human tumors. The most common site of lung cancer metastasis is the brain, and about 50% of all lung cancers develop into brain metastasis (BM) during the disease process (2,3). It has been reported that the rate of brain metastasis from lung cancer has recently increased, placing a great burden on public health services.

The ‘seed-and-soil’ hypothesis, which is the most widely accepted hypothesis for the formation of metastasis, denotes that the growth of metastatic cancer cells depends on the intrinsic abilities of the cancer cells themselves (‘seeds’) and the target organ microenvironment (‘soil’) (4). The cancer cell population has multiple genetically heterogeneous subpopulations (5). Metastasis is a Darwinian natural selection process in which ‘seeds’, with distinct metastatic traits that enable them to obtain metastatic advantage, are selected from a genetically- and epigenetically-heterogeneous tumor cell subpopulation (6,7). The advantageous lung ‘seeds’ proliferate in brain ‘soil’ that provides a congenial ground and form metastatic brain tumors, whose genetic landscape is reshaped (8).

A large number of studies have attempted to weigh BM relevant factors in lung cancer, revealing factors that include young age (<60 years) (9,10), non-squamous cell carcinoma (10,11), and the presence of clinically bulky mediastinal lymph nodes (≥2 cm) (9) are associated with a high BM rate. However, other studies have reported conflicting results (12,13). Many candidate metastasis genes have also been found to be involved in metastasis through changes in gene expression levels (14). The expression levels of E-cadherin, N- cadherin, KIFC1, and FALZ may be used to identify patients at high risk of lung-brain metastasis (15,16). However, the molecular basis of metastatic gene expression remains largely unknown, and the genetic profiles of brain metastases from LC might give us a close insight into tumor initiation, dissemination, and local progression (14). It is universally acknowledged that the development of lung cancer determined by gene mutations and these oncogenes play as critical drivers in tumor metastasis. The identification of driver oncogenes and the development of targeted therapies have improved lung cancer patients’ outcomes. WES provide a convenient approach for understanding tumor mutation, mutational burden and further guide therapeutic decision making. To reveal the molecular mechanisms and the genetic alterations involved in metastasis of lung tumors to the brain, we carried out whole-exome sequencing (WES) of the primary tumors and their corresponding BMs from 12 patients with metastatic non-small-cell lung carcinoma. Our study might be instrumental for the identification of new genetic targets which may provide new therapeutic strategies for the design of drug intervention to improve the severity of the disease.

We present the following article in accordance with the STREGA reporting checklist (available at http://dx.doi.org/10.21037/atm-21-1555).

Methods

Patients and specimens

The pairwise lung-brain tumor samples and adjacent histologically normal tissue samples from 12 patients were collected at the Sun Yat-sen Memorial Hospital from 2010 to 2015. All procedures performed in this study involving human participants were in accordance with the Declaration of Helsinki (as revised in 2013). The Ethics Committee at Sun Yat-sen Memorial Hospital approved the utilization of samples. All informed consents were taken from 12 LC patients. A total of 12 paired samples were subjected to hematoxylin and eosin (HE) staining, then the normal cells were isolated from the tumor cells through histopathological examination. The normal cells were regarded as matched normal controls. The total amount of DNA extracted from the archived formalin-fixed paraffin-embedded (FFPE) samples of tumor tissue was up to standard and qualified.

Next-generation sequencing, variant calling, and annotation

The genomic DNAs were exacted and sonicated to an average size of 200 bp. The targeted DNA fragments were captured and pulldown and exon-wide libraries were created using the Roche SeqCap EZ Exome V3 (Roche Diagnostics, Basel, Switzerland) and TruePrep DNA Library Prep Kit V2 for Illumina (#TD501, Vazyme, Nanjing, China), and paired-end sequence data were generated using Illumina HiSeq machines (Illumina, Inc., San Diego, CA, USA). The sequence data were aligned to the human reference genome (NCBI build 37) using BWA (Burrows-Wheeler Aligner) (17) and sorted and removed polymerase chain reaction (PCR) duplication using GATK 4.1.2.0 (https://gatk.broadinstitute.org/) (18). Somatic mutation calling was performed using Mutect1, Mutect2 (18), and VarDict (19). Somatic mutations existing in at least 2 of the results of the 3 software were selected as high confident mutations and to be involved in the further bioinformatics and bio-statistical analysis. Copy number variants (CNVs) from WES data were detected by CNVKIT (Genome-Wide Copy Number Detection and Visualization from Targeted DNA Sequencing) (20). The dN/dS (non-synonymous substitution rate/synonymous substitution rate) ratios of each domain were calculated via The DiversiTools according to Xia’s description (21). The GISTIC2.0 (Genomic Identification of Significant Targets in Cancer) algorithm was used to identify regions of the genome that are significantly amplified or deleted across a set of samples (22). Somatic variants were annotated by Ensembl Variant Effect Predictor (23). Transition (Ti) and transversion (Tv) ratios were applied to measure the selection in cancer genomes and to show mutation characteristics between different cancer types.

Spectrum, signatures of somatic mutations, and subclonal architecture inference

Non-negative matrix factorization and model-selection were applied to delineate mutational processes underlying genome-wide single nucleotide variants (SNVs) and to identify the major mutational signatures (24). In the cohort, 1 sample was considered as a strong association with 1 mutational signature if the proportion of the contribution >20% with MutationalPatterns (version 1.10) (25) and deconstructSigs (26). The sub-clonal architecture of tumors was inferred by sciClone (27) and clonevol (28).

Survival analysis, enrichment analysis

Pathway analysis was based on the Database for Annotation, Visualization and Integrated Discovery (DAVID) bioinformatics (29), and WebGestalt (30) with significant Benjamini adjusted P value (P<0.05). The co-mutation profile was prepared with R package ComplexHeatmap (31). The Cancer Genome Atlas (TCGA) mutation and survival data were downloaded from the Genomic Data Commons (GDC) database (https://portal.gdc.cancer.gov/exploration). In the validation study to ERF, we downloaded the expression and survival data of ERF mutation and expression from TCGA project, and Cox-regression was conducted to binary gene expression data dichotomized by median expression level for TCGA dataset.

Statistical analysis

Cox regression was applied for survival analysis between mutation and overall survival (OS) time and the Kaplan-Meier plot were used to show the difference between survival time among different groups. Since our research involved a small sample size of participants, all the clinically-related statistical analysis was considered to be significant when P<0.05 without multiple correction testing. Data is shown as the mean ± standard deviation (SD), and a two-tailed P value <0.05 was considered statistically significant.

Results

Landscape of somatic mutations in primary LC and BM

We collected and quantified DNA from the original 12 non-small cell lung cancer (NSCLC) patient FFPE samples and the matched BM samples. The average coverage depth for the tumor cells and normal cells were 194× and 120×, respectively. Detailed clinicopathological information is summarized in Table S1. We identified 1,702 SNVs and 6,131 mutation events in 1,220 genes from 12 paired LC and BM, including the most frequent LC driver gene mutations such as TP53, EGFR, BRCA1, BRCA2, and BRAF. We identified a mean of 3.1 driver gene mutation events per tumor with the dN/dS of 2.13 which is slightly higher than non-metastatic LC samples in TCGA, indicating a significant enrichment for the cancer driver gene mutations. We did not find any difference in the dN/dS ratio between primary tumor (dN/dS =2.20), BM tumor (dN/dS =2.06), and shared mutations between LC and BM (dN/dS =2.25).

We identified several LC metastases-associated genes (KMT2C, AHNAK2, PDE4DIP, ANKRD36C, and BAGE2), and the mutations of these genes showed distribution diversity among the LC and BM samples (). Mutations of KMT2C were found in 25% of LC samples; however, the mutation frequency in BM was up to 50%, indicating the positive selection of KMT2C mutations during metastasis. According to TCGA dataset AHNAK2 have significant enrichment in LC, with a mutation ratio of 18.8% in LC while 9.98% in pan-cancer (P=7.2×10−9, Chi-square test). However, the mutation frequency of AHNAK2 in our dataset is as high as 26.9% which is 1.43 times that of the LC population in TCGA dataset (P=0.02, Chi-square test). We also observed all EGFR mutations were shared between LC and BM, suggesting that EGFR mutations are drivers and likely to be an early event before BM.

Figure 1

The mutational landscape and patterns of genetic divergence in paired primary lung cancer and brain metastases. (A) Top mutated genes in LC and brain cancer samples; (B) tumor-specific mutations and lung-brain shared mutations; (C) drive mutation genes in tumor specific mutations and shared mutations; (D) oncoprint of lung and brain mutations in all the samples. LC, lung cancer.

We found more somatic mutations in BM lesions (median 71, range 23–180) than in LC lesions (median 48.5, range 13–187), while the difference was not statistically significant (P=0.069, Student’s t-test) (Figure S1A,B). High correlation between TMB (tumor mutation burden) of LC and TMB of BM were confirmed by Pearson coefficient 0.65 (P=0.02) (Figure S1C,D), indicating the TMB of BM can be estimated by that of primary LC when brain tissue is not available, so as to screen patients who will most likely benefit from programmed death-ligand 1 (PD-L1) immunotherapy. A sum of 18.2% (0.5–35.9%) of all mutations were shared between LC and BM, clearly suggesting a common ancestral truncal clone with 30.0% (9.3–60.8%) LC-specific and 51.8% (19.9–79.7%) BM-specific, respectively (). Although metastases had more private mutations than the primary tumor, they were not enriched for the pan-cancer driver genes (32) (). This suggested that few additional private genomic driver genes were required for metastasis when the primary cancer was already advanced.

In order to provide more landscape for the mutations identified in our study, we conducted a pathway analysis to the most frequently mutated genes (mutation frequency >5%) (Table S2). We found the most important tumor pathways were chronic myeloid leukemia (P=0.002), ErbB signaling pathway (P=0.0014), and glioma pathway (P=0.05). Keyword enrichment indicated important metabolic abnormalities for the lung-metastasis cancers including EGF-like domain and tyrosine-specific phosphatase (Table S3).

Copy number variations (CNVs)

To further explore BM-related molecular events, genomic CNVs were analyzed: 8q21.2, 6p22.1, 12p13.33, and 5q35.3 were the most common chromosome deleted regions in both LC and BM, and 8q24.13 were the most common regions with gain copy numbers in both LC and BM (Figure S2A). Loss of 6p22.1, which harbors HLA-A, HLA-G, and HLA-H, was most frequent in both the LC and BM of participant 9 and in BM of participant 11 (Figure S2B). Interestingly, these samples also had high TMB (Figure S1C). This may be due to the loss of HLA function associated with higher overall mutation burden and a larger fraction of HLA-binding neoantigens (33). The recurrent deletion of HLA was detected as the early events, indicating the important role of the immune system in LCBM, and that these patients may benefit from immunotherapy.

The differences include gains of chromosomes 7q35 and loss of 7q22.1 and 7q36.3, which were more frequent in metastasis samples, and gains of chromosomes 11q13.2 and losses of 7q11.23 and 2q13, which were less frequent. Most recurrent CNV regions were shared in LC and BM samples, indicating that CNVs are early molecular events in tumorigenesis and metastasis.

Spectrum and signatures

To determine the relationships between the mutational spectra and tumor organ sites, we analyzed the spectra of LC and BM from our study and primary lung adenocarcinoma (LUAD), lung squamous cell carcinoma (LUSC), low grade glioma (LGG), and glioblastoma multiforme (GBM) from TCGA dataset. The C to T substitution was the most common transition in our LC and BM samples, which was much closer to primary brain cancer (LGG and GBM) but significantly different from primary lung cancer (LUAD and LUSC) which has higher C to A transversion (). These evidences were consistent with our hypothesis that the mutations identified in our study have a higher probability of being associated with BM. The mutational spectrum of LC and BM samples from the same individuals are more similar to each other than those from different patients, implying that different mutational processes were involved during the development of metastasis between the different patients ().

Figure 2

Spectrum and signatures. (A) Mutation spectrum of lung tumor specific mutations, brain tumor specific mutations, shared mutations, and TCGA lung and brain mutations; (B) mutation spectrum of all the samples in our cohort; (C) mutation signatures of mutations for lung tumor samples and brain tumor samples. TCGA, The Cancer Genome Atlas.

We further analyzed mutational signatures in BM and LC and signatures 1, 3, and 4, which have been linked to aging, BRCA1/2 mutations, and smoking, respectively, were identified as dominant in either BM or LC samples (34) (). There was no significant difference between the signature levels of LC and BM tissues (Wilcoxon rank sum test, P>0.05), indicating the change of mutational signatures happened before metastasis and may not lead to their great difference between the two groups.

Clonal evolution during the development of LCBM

Phylogenetic trees give clear overviews of the order of mutation events, allowing the track of emergence and movement of clones from LC to BM (). Phylogenetic trees of the 9 participants showed that mutations on the trunk were probably earlier genetic alteration events, followed by mutations on the branch which occurred later during tumorigenesis and BM development. Clonal evolution analyses revealed that LC and BM tumors had the same evolutionary process in 3 participants (P02, P07, and P10), LC tumors harbored a cluster of LC-private clones in another 3 participants (P05 and P08), and BM tumors harbored clones that were nonexistent in matched LC tumors in 3 other participants (P06, P13, and P15), indicating the mutations on BM-private clones may contribute to metastatic progression.

Figure 3

Phylogenetic trees indicating 4 types of evolution. (A) LC tumors and BM tumors had the same evolutionary process; (B) BM tumors harbor a cluster of BM-private clones; (C) LC tumors harbor a cluster of LC-private clones; (D) LC and BM tumors harbor a cluster of LC-private clones and a cluster of BM-private clones respectively. LC, lung cancer; BM, brain metastasis.

OS by genotype

In order to identify independent prediction factors for outcomes, we conducted survival analysis on several potential factors. Participants with aberrations of 18 genes in LC (Figure S3) and 15 genes in BM (Figure S4) had significantly worse OS than those without these aberrations (P<0.05). Of these genes, we identified a significant survival-associated mutation gene, ERF, which was confirmed by both TCGA (P=0.01) and our dataset (P=0.012) (). Additionally, in order to show the prognostic roles of ERF, we also found that high expression of ERF genes in TCGA was a significant risk factor for the OS time [hazard ratio (HR) =1.46, P<1.2×10−22, ]. Taken together, our findings revealed an important role for ERF in prognostic prediction of LC.

Figure 4

Overall survival analysis. (A) Survival analysis of several potential factors; (B) survival analysis between ERF mutation and OS times in TCGA; (C) survival analysis between ERF expression and OS times in TCGA; (D) survival analysis between ERF mutation and OS times in our cohort. *, P<0.05; **, P<0.01. OS, overall survival; TCGA, The Cancer Genome Atlas.

Multivariate analysis demonstrated gender (P=2.02×10−119), smoking status (P=1.21×10−269), metastatic tumor size (P=0), and the ratio of shared mutations in lung and brain cancers (P=0.019) were significantly associated with OS time, while no significant associations were found in drinking status (P=0.996), the number of metastatic tumors (P=0.746), mutation numbers of primary tumor (P=0.840), or metastatic tumor (P=0.248) ().

Discussion

The metastatic cascade involves multiple steps, including invasion, entry into the circulation from the primary tumor, systemic dissemination, arrest and extravasation in secondary organs, settlement into latency, reactivation, outgrowth, and potential seeding of tertiary metastasis (35). Genetic and epigenetic changes accumulating among primary tumor cells and metastases may contribute to these multiple steps of the metastatic cascade (36,37). Whole-exome sequencing is a convenient to acquire comprehensive characterization of tumor mutations, from which neoantigens, tumor mutational burden (TMB), and clonality can be obtained. Even though there are previous study characterized the genomic difference between primary tumors and paired brain metastases in lung cancer patients by whole-exome sequencing, this study has larger sample size in comparison to previous study; survival analysis between mutation and overall survival (OS) time was conducted (38).

Spinal metastasis, hepatic metastasis, bone metastasis, pleural dissemination is the important modes of metastasis and a poor prognostic factor of lung cancer. Tang using SCIMET tool to value timing of lung cancer dissemination and found the metastatic seeding happened approximately 2.74 years before clinical detection and 87.5% non-lymph node metastases were mainly seeded by the primary tumors (39). Hence, it is necessary to collect a well-defined cohort of matched primary tumors and BM and perform comparative deep sequencing analyses to acquire some biomarkers of metastasis. However, most patients with brain metastasis of lung cancer (LCBM), which were in the late stage at diagnosis, are typically treated with palliative approaches such as chemotherapy, targeted therapy, and whole-brain radiotherapy instead of neurosurgical resection. Thus, researchers rarely have the opportunity to investigate matched primary-metastatic tumors in the mutation status of analyzed genes between tumor sites on a large-scale (40).

In this study, we collected 12 paired LC and BM samples and identified some genes associated with LCBM. The gene KMT2C (lysine-specific methyltransferase 2C, also known as MLL3), which belongs to the mixed-lineage leukemia (MLL) family of histone methyltransferases, was the most commonly mutated gene among our samples. Recent studies have revealed frequent mutations of KMT2C in several epithelial and myeloid cells, and it has been identified as a putative tumor suppressor (41,42). It has been proposed that in metastatic spread, circulating tumor cell (CTC) populations in the blood of carcinoma patients contain cells with the clonal capacity to initiate metastatic growth in distant organs, thus having similarity with the hematologic tumor. The involvement of KMT2C was originally identified as oncogenic fusions in leukemia (43), and the most enriched pathway in our study was chronic myeloid leukemia. Aljohani et al. reported that the Keap1-Nrf2-ARE pathway was mutated in NSCLC patients that had metastasized to the brain and in CTC according to WGS (44). Recent reports have demonstrated that metastatic brain tumors of NSCLC and LUAC patients carried higher EGFR mutations rate than those without metastases (45,46). Our results further emphasized the association between the BM of LC and leukemia. The gene AHNAK nucleoprotein 2 (AHNAK2) is a prognostic marker and an oncogenic protein for clear cell renal cell carcinoma and hypoxic upregulation of AHNAK2 supports epithelial-mesenchymal transition (EMT) and cancer cell stemness (47). Cancer cells acquire characteristics of self-renewal, motility, and invasiveness, traits that facilitate metastatic dissemination during EMT (48). That is, the driver gene mutation of AHNAK2 may promote metastatic colonization of the lung to the brain by supporting EMT.

Our findings, depicting that chronic myeloid leukemia and the ErbB signaling pathway were mutated in the majority of LCBM patients, supports our hypothesis that mutations in these pathways may indeed provide a survival advantage to these cells and help them reach distant sites. Of note, glioma pathways have also been identified. The Ti/Tv profile showed our mutation profile was much closer to brain cancer mutation profile as brain cancer has high C to T transitions were more frequent, while LC had a higher frequency of C to A transversion. These evidences strongly indicate that brain tumor-related events are involved in the process of LCBM.

Multivariate analysis demonstrated that the high ratio of shared mutations were associated with better prognosis (49,50). This implies that patients might have had a preferable prognosis when the tendency of metastatic cancer mutation was more inclined to primary cancer, that is the BM sample might not have evolved from the primary cancer but rather they had a shared antecedent. Limited inter-tumor heterogeneity between LC and BM within the same patient results in effective postoperative chemotherapy and radiotherapy in this circumstance.

Overall, we revealed the genomic difference between metastatic and primary tissues, and identified several genes associated with LCBM. While further molecular biology studies to validate the role of identified candidates and a larger-scale of LCBM samples would be required to confirm our findings, this study provides a preliminary evidence of the genetic evolution of LC metastases and furthers our advantage in revealing the therapeutic vulnerabilities of LC metastatic tumors.

The article’s supplementary files as