To explore the importance of NTRK3 mutations in LUAD, we investigated the correlation between NTRK3 mutations and the outcome of ICI treatment in patients. First, a discovery immunotherapy cohort from Samstein  et al, which consisted of clinical and mutation data for patients receiving ICI treatment, was collected (Samstein et al., 2019). Then, after excluding five patients without mutation data (n=5), we divided the ICI treatment cohort with mutation data (n=266) into NTRK3-MT (mutant-type) and NTRK3-WT (wild-type) groups according to the nonsynonymous somatic mutation status of NTRK3 and then used Kaplan-Meier survival curves for analysis. Moreover, an R/Bioconductor package called TCGAbiolinks (Colaprico et al., 2016) was employed to download somatic mutation and OS data (n=494) from the Genomic Data Commons portal (https://portal.gdc.cancer.gov/) and the LUAD cohort from TCGA. For the latter, Kaplan-Meier survival curves were used to show the differences in OS (data from TCGA-LUAD) and disease-free survival (DFS) between patients with NTRK3-MT and NTRK3-WT. The DFS data (n=430) downloaded through cBioPortal (Cerami et al., 2012). In addition, we analyzed the association of OS with TKI-sensitive gene mutations in the ICI treatment cohort.

All samples with somatic mutations (n=266) reported by Samstein  et al. were analyzed by targeted next-generation sequencing (NGS) and evaluated with the Memorial Sloan Kettering Cancer Center-Integrated Mutation Profiling of Actionable Cancer Targets (MSK-IMPACT) test. NAL(neoantigen load) data (n=500) from TCGA-LUAD have been reported (Gu et al., 2016). In addition, we analyzed whole-exome sequencing (WES) and drug sensitivity response data for LUAD cell lines in the Genomics of Drug Sensitivity in Cancer (GDSC) database (Yang et al., 2013). Consistent with approaches used for published data (Chalmers et al., 2017), nonsynonymous mutations in TCGA-LUAD were used as raw mutation data and divided by 38 Mb to quantify the TMB. For the ICI treatment and TCGA-LUAD cohorts, the R package ComplexHeatmap (Gu et al., 2016) was applied to visualize mutations in the genes with the top 20 mutation rates and the clinical characteristics associated with these mutations. The R package Maftools (Mayakonda et al., 2018) was used to visualize NTRK3 mutation sites.

A previous study showed the predictive relationship between immune cell-related genes and chemosensitivity (Newman et al., 2015). Therefore, we used CIBERSORT (Hao et al., 2018) (http://cibersort.stanford.edu/) to examine the gene expression data (RNA-seq with the Illumina HiSeq platform) for the LUAD cohort from TCGA (n=505) using TCGAbiolinks. Then, we compared the infiltration of 22 types of immune cells between NTRK3-MT and NTRK3-WT LUAD. We also examined differences between NTRK3-MT and NTRK3-WT LUAD with regard to expression of immune-related gene at the mRNA level. As reported in some studies (Thorsson et al., 2019), immune-related genes, the expression level was quantified as fragments per kilobase of exon model per million mapped fragments (FKPM) values and log2 transformed, along with their functional classification and immune-related scores, were obtained. Data for LUAD cell lines were downloaded from GDSC and compared to determine the difference in drug sensitivity between NTRK3-MT and NTRK3-WT LUAD.

Broad GDAC Firehose (http://gdac.broadinstitute.org/) was used to download Affymetrix SNP 6.0 microarray data for TCGA-LUAD (hg19; germline/potential false-positive calls were removed) and the GISTIC2.0 module of GenePattern was utilized to evaluate downloaded copy number alteration (CNA) segments (Reich et al., 2006) (https://cloud.Genepattern.org/gp/pages/index.jsf). Default parameters were used (except, for example, the confidence level was set to 0.99, and the X chromosome was excluded before analysis). The R package Maftools was used to visualize the GISTIC2.0 CNA analysis (Mayakonda et al., 2018).

Differentially expressed RNAs (raw count) in the LUAD cohort from TCGA derived from TCGAbiolinks were identified using the R package edgeR (Robinson et al., 2010). The R package clusterProfiler (Yu et al., 2012) was used for gene annotation enrichment analysis. A P value of less than 0.05 was considered to indicate a significant difference in Gene Ontology (GO) terms, Kyoto Encyclopedia of Genes and Genomes (KEGG), and Reactome analyses. The data used for gene set enrichment analysis (GSEA) were obtained from Broad Institute Molecular Signatures Database (MSigDB) (Subramanian et al., 2005); gene sets of DNA damage response (DDR)–related pathways were obtained from the Broad Institute MSigDB collection (Subramanian et al., 2005). The DDR gene set was used to evaluate the number of nonsynonymous mutations in the immunotherapy (n=265) and TCGA-LUAD cohorts (n=514) and GDSC-LUAD cell lines and to identify differences in the number of nonsynonymous mutations in the DDR pathway between NTRK3-MT and NTRK3-WT LUAD.

Univariate Cox regression analysis was performed to identify the prognostic role of NTRK3 mutations and other common TKI-sensitive gene mutations in ICI treatment. Correlations between the NTRK3 status and the TMB, NAL, abundance of immune cells, expression of immune-related genes, age, smoking history (pack years) and number of gene mutations in the DDR pathway were assessed using the Mann-Whitney U test. Fisher’s exact test was applied to assess differences in the mutation status of the top 20 mutated genes in the immunotherapy cohort, sex and smoking history between patients with NTRK3-MT and NTRK3-WT LUAD. Fisher’s exact test was also employed to evaluate differences in sex, race, ethnicity, smoking history and clinical stage between patients with NTRK3-MT and NTRK3-WT in the LUAD cohort from TCGA. Survival curves were generated using the Kaplan-Meier method with the log-rank test. A P value of less than 0.05 was considered to indicate a significant difference. All statistical tests were two-sided. Statistical and visual analyses were carried out by using R software (version 3.6.1). The R package ggpurb was used to generate boxplots (Kassambara, 2018). CNAs with false discovery rate (FDR) of 5% as a cutoff criterion were visualized.

To investigate the potential role of key gene mutations in immunotherapy, we first collected a discovery immunotherapy cohort (LUAD, n=271) from Samstein et al, comprising clinical and mutation data for patients receiving ICI treatment; mutation data of the LUAD cohort from TCGA were also downloaded. Our bioinformatic analysis workflow is shown in Supplementary Figures S1A–C . The Venn diagram depicted in Figure S1D displays the sample size of different data sets. As the schematic representation shows, screening revealed found that 16 of the 266 patients with mutation data in the above cohort had the NTRK3-MT genotype ( Supplementary Figure S1A ). Furthermore, we estimated the ability of NTRK3-MT to predict an OS benefit in both the ICI-treated and TCGA-LUAD cohorts and to predict a DFS benefit in the latter ( Figure 1A ). Univariate survival analyses indicated that the OS of the two groups, namely, NTRK3-MT and NTRK3-WT, differed significantly in the ICI treatment cohort (n =266; P=0.01; HR = 0.3; 95% CI, 0.17 to 0.54). However, Kaplan-Meier results indicated no significant difference in OS (n =494; P=0.415; HR = 1.28; 95% CI, 0.66 to 2.47) or DFS (n =430; P=0.533; HR = 1.18; 95% CI, 0.67 to 2.08) between the NTRK3-MT and NTRK3-WT groups in the non-ICI-treated TCGA-LUAD cohort. In addition, we conducted univariate Cox regression analysis to demonstrate that compared to other common TKI-sensitive gene mutations (BRAF, EGFR, KRAS, PIK3CA), only NTRK3-MT [hazard ratios (HR)=0.30, 95% confidence interval (CI)=0.11–0.80, P <0.05] correlated with a good prognosis to ICI treatment in the ICI-treated cohort ( Figure 1B ). Although we also found that the TMB score (HR=0.97, 95% CI=0.95–0.99, P <0.05) may serve as a biomarker for a good prognosis, its HR in in the ICI-treated cohort was very close to 1, which suggested that it is not a reliable predictor in this context. To date, it remains unclear whether the TMB can identify patients who may benefit from ICI treatment (Ferrara et al., 2017; Chalmers et al., 2017). In addition, survival analysis for TKI-sensitive gene mutations (BRAF, EGFR, KRAS, PIK3CA) showed no significant difference between their wild-type and mutant types in the ICI-treated cohort ( Figure 1C ). Overall, the above results clearly indicate the predominant role of NTRK3-MT in immunotherapy, which is not affected by TKI-sensitive gene mutations.

We then explored whether the NTRK3 mutation status is associated with clinical characteristics. Figures 2A, B show the top 20 most significantly mutated genes and the clinical characteristics in the ICI treatment and TCGA-LUAD cohorts. No significant differences were observed in the clinical characteristics of NTRK3-MT and NTRK3-WT patients treated with ICIs except for the TMB and OS. The results for the dataset TCGA-LUAD revealed that NAL, TMB and smoking history correlated significantly with the NTRK3 mutation status but that the age, sex, stage, race, and ethnicity did not. Interestingly, we also found that the mutation rate of Tp53 was higher in patients with NTRK3-MT, consistent with prior reports that Tp53 mutations are associated with enhanced antitumor immunity in LUAD (Li et al., 2020). Overall, these results highlight a potential role of NTRK3-MT as a predictive biomarker for ICI treatment. Furthermore, we annotated each mutation in NTRK3 in both the LUAD-MSKCC panel and the TCGA-LUAD in a lollipop plot ( Figure 2C ). The data from both showed that most mutations in NTRK3 occur in the protein tyrosine kinase, mainly in the immunoglobulin I-set domain and fifth domain (immunoglobulin−like) of the Trk receptors TrkA, TrkB, and tropomyosin-related kinase C (TrkC). In addition, more NTRK3 mutations in the fourth domain (immunoglobulin−like) of the Trk receptors TrkA, TrkB, and TrkC and the leucine-rich repeat C-terminal domain and catalytic domain of the protein tyrosine kinase TrkC were found in the TCGA-LUAD than in the LUAD-MSKCC panel. For the former, differences in somatic CNAs between the NTRK3 mutation statuses were assessed using GISTIC2.0. As illustrated in Figure 2D , amplifications on chromosomes 1, 8, 12, and 20 were enriched in the NTRK3-MT group.

We further evaluated differences in antitumor immunity and immunogenicity between the NTRK3-MT and NTRK3-WT groups in the LUAD cohort from TCGA. Some immune activation-related genes, such as interferon genes (IFNA1) and interleukin genes (IL1A), were significantly upregulated in the NTRK3-MT group, whereas the immunosuppression-related gene KIR2DL1 was downregulated ( Figure 3A ). In addition, stimulatory immunomodulators such as chemokines (CCL5, CXCL10, and CXCL9), cytolytic activity-associated genes (GZMA and PRF1) and immune checkpoint biomarkers (CTLA4, LAG3, PDCD1, and TIGIT) exhibited higher expression in the NTRK3-MT group than in the NTRK3-WT group ( Figure 3B ). Then, the TMB was evaluated in the ICI-treated LUAD and TCGA-LUAD cohorts, and as expected, it was significantly higher in the NTRK3-MT groups ( Figure 3C ), consistent with previous studies showing that tumor immunity is positively associated with the TMB. Additionally, a higher NAL correlated with NTRK3 mutation in the LUAD cohort from TCGA ( Figure 3C ).

Moreover, we investigated immune cell-related genes in the LUAD cohort from TCGA and found that genes that correlated with antitumor immunity (NEIL3, GAL, FOSL1, and BASP1) were more highly enriched in the NTRK3-MT group ( Figure 3D ). Infiltration of 22 types of immune cells was compared in the NTRK3-MT and NTRK3-WT groups of the LUAD cohort from TCGA ( Figure 3E ) and these immune cells were generally more abundant in the NTRK3-MT group.

Collectively, these results confirm that NTRK3 mutations correlate with enhanced antitumor immunity and immunogenicity in LUAD and again indicate the possible predictive value of NTRK3-MT for immunotherapy.

To further explore whether NTRK3 mutations are involved in important tumor-related biological pathways and processes, we compared levels of signaling pathway enrichment between the NTRK3-MT and NTRK3-WT groups ( Figure 4 ). Pathways correlating with tumorigenesis and development, such as cell cycle arrest and apoptosis, negative regulation of NOTCH4 signaling, regulation of RAS by GAPs and positive regulation of insulin receptor signaling were more highly enriched in the NTRK3-MT group. In addition, signaling pathways related to the cholesterol biosynthetic process and very long-chain fatty acid biosynthetic process pathways differed significantly in the NTRK3-MT group. However, oncogenic pathways such as ATF6 (ATF6-alpha) activate chaperones, GPCR ligand binding and RHO GTPases activate NADPH oxidases were downregulated in the NTRK3-MT group. Collectively, NTRK3-MT is associated with these pathways, reminiscent of its role in tumorigenesis.

Reports have revealed the association of genomic instability and increased immunogenicity with the DDR pathway. To identify whether NTRK3 mutations function in DDR pathways, we evaluated alterations in various DDR pathways, such as the MMR, base excision repair (BER), nonhomologous end-joining (NHEJ), double-strand break (DSB) repair, single-strand break (SSB) repair, nucleotide excision repair (NER), and homologous recombination repair (HRR) pathways. More genomic alterations in these DDR-related pathways were observed in the NTRK3-MT group than in the NTRK3-WT group in both the ICI-treated LUAD and TCGA-LUAD cohorts ( Figures 5A, B ). We further compared the response to anticancer drugs between the NTRK3-MT and NTRK3-WT groups ( Figure 5C and Supplementary Figure S2 ) and found statistically significant differences for the three drugs.