Comprehensive analyses reveal TKI-induced remodeling of the tumor immune microenvironment in EGFR/ALK-positive non-small-cell lung cancer.

Tyrosine kinase inhibitors (TKI) play a pivotal role in the treatment of non-small-cell lung cancer (NSCLC) with mutations in epidermal growth factor receptor (EGFR) and rearrangements in anaplastic lymphoma kinase (ALK). However, the influences of TKIs on the tumor immune microenvironment (TIM), especially dynamic changes of responders, have not yet been fully elucidated. Therefore, RNA sequencing and whole-exome sequencing were performed on EGFR/ALK-positive NSCLC samples before and after TKI treatment. In combination with neoantigen and mutational-load estimations, xCell and single-sample gene set enrichment analysis (ssGSEA) were used to assess tumor immune-cell infiltration and activity. Furthermore, weighted-gene correlation network analysis and the bottleneck method were used to identify the hub genes that affected treatment-related immune responses. We found that TKI treatment remodeled the TIM in treatment-responsive samples. Profound increases in the rate of anti-tumor cell infiltration and cytotoxicity was observed following TKI treatment, while antigen presentation was limited in ALK-rearranged samples. However, no significant change in anti-tumor cell infiltration or cytotoxicity was found between pre-treatment and post-progression samples. Subsequently, we found that neurofilament heavy (NEFH) mutations were enriched in samples after TKI treatment and were associated with reduced neutrophil infiltration. The cytotoxicity of EGFR-mutant NSCLCs with co-driver TP53 mutation and ALK-rearranged samples with wild-type TP53 seems to be more easily induced by TKI. Finally, the immune-associated score generated by hub genes was positively correlated with immune infiltration, immune activation, and a favorable prognosis. In conclusion, the dynamic changes in the TIM provide clues to drug selection and timing for TKI-immunotherapy combinations.

Introduction

Lung cancer is the leading cause of cancer-related death worldwide, and non-small-cell lung cancer (NSCLC) is the major histological subtype of lung cancer.[1] Nearly two-thirds of NSCLC patients have carcinogenic driver mutations, approximately half of which have treatable target lesions. These expand the scope of treatment options and may help improve survival and safety over conventional chemotherapy.[2] At present, most of the targetable oncogenic mutations in receptor tyrosine kinases occur in lung adenocarcinomas (ADCs). These include epidermal growth factor receptor (EGFR) mutations and anaplastic lymphoma kinase (ALK) rearrangements.[3]

Although the therapeutic effects of tyrosine kinases inhibitors (TKIs) have improved the prognosis of NSCLC patients, the responses to TKIs are usually incomplete and temporary.[2] There has been strong interest in determining whether immune-checkpoint inhibitors may produce more complete and longer-term efficacious responses in NSCLC patients. However, findings on the efficacies of immune-checkpoint inhibitors either after TKI treatments or in combination with TKI treatment have been disappointing.[4,5] Several explanations for these poor responses to immunotherapy in NSCLC patients harboring EGFR mutations or ALK rearrangements have been proposed – such as variable programmed death-ligand 1 (PD-L1) expression,[6] low mutational load, and an inactive tumor immune microenvironment (TIM) characterized by a lack of infiltration of cytotoxic T cells[4,7] – all of which may negatively affect responses to PD-1/PD-L1 inhibitors. In addition, TKIs may inhibit the expression of PD-L1 up-regulated by EGFR or ALK activation,[8,9] thus indicating that TKIs and anti-PD-1 antibody may have a similar effect on blocking the interaction between PD-1 and PD-L1, but not a synergistic effect.[10,11] Therefore, there is a need to further elucidate the effects of TKI treatments on the function of the immune system. Previous studies showed that ALK inhibitor crizotinib and anti-EGFR monoclonal antibody cetuximab induce immunogenic cell death in NSCLC and are linked to increased T-lymphocyte infiltration.[12,13] However, evidence has suggested that EGFR-TKIs decrease T-cell infiltration in post-TKI-resistant samples with EGFR mutations[14] and that there is no significant difference in the abundance of CD8-positive (CD8+) T cells between naïve-treatment and post-TKI-resistant mouse samples with ALK fusion.[15] Although there are a limited number of studies on TKI-induced TIM reprogramming of TKI-resistant samples, the immunological effects of TKIs in treatment-responsive patients remain widely unknown.

Herein, using the whole-exome sequencing (WES) and RNA sequencing (RNA-seq) of the biopsies of pre-TKI and post-TKI responders, we present results assessing TKI-induced TIM remodeling and findings regarding differential patterns of immune responses between EGFR-mutant and ALK-rearranged samples via several methods, such as computing the signature score of immune cells and tumor-neoantigen estimation. We also analyzed and compared our data with public datasets to identify immune-related genes underlying TKI-induced dynamic changes of the TIM in response to TKI treatments.

Materials and methods

Data availability

Our ADC datasets generated and analyzed during the current study are not publicly available but are available from the corresponding author on reasonable request. Data from three publicly available datasets were incorporated into our study. TCGA level three RNA-seq data and clinical information from patients were acquired from the UCSC Xena website (https://xenabrowser.net/). Additionally, microarray data from 34 lung cancer samples with EGFR mutations from the GSE11969 dataset were downloaded from the Gene Expression Omnibus (GEO) (http://www.ncbi.nlm.nih.gov/geo). RNA expression data of paired pre- and post-TKI-resistant samples were available for the Osimertinib-treated cohort.[16]

Patients

We enrolled 23 patients with histologically confirmed EGFR mutations or ALK rearrangements in advanced lung ADCs as defined by the American Joint Committee on Cancer guidelines. Fifteen paired pre- and post-TKI samples with EGFR mutations and eight paired ALK-rearranged tissues were obtained from patients undergoing TKI treatments. With the help of computed tomography (CT) and the Response Evaluation Criteria in Solid Tumors (RECIST 1.1),[17] three post-treatment samples were identified as progressive disease (PD), and the other 20 post-treatment samples were defined as partial response (PR). In total, 45 samples (including three samples from P3 patient and four samples from P4 patient) were obtained for RNA-seq analysis, and 42 samples were used for WES analysis. Patient information is summarized in Supplementary Table S1. This study was approved by the ethics committee of NangFang Hospital (NFEC-2019-265), and informed consent was obtained from all patients.

RNA sequencing (RNA-seq)

Total RNA from fresh frozen tissues was extracted with Trizol according to the manufacturer’s recommendations. Sequencing libraries were generated using a KAPA Stranded RNA-Seq Kit (KAPA, USA) following the manufacturer’s instructions, and index codes were added to attribute sequences to each sample. Libraries were pooled and paired-end sequencing (2 × 150 bp reads) was performed on an Illumina X-ten PE150 platform.

Whole-exome sequencing (WES)

Fresh frozen tissues from tumor samples were used for genomic DNA extraction via a DNeasy Blood and Tissue Kit (QIAGEN) following the manufacturer’s instructions. Peripheral blood DNA extracted from individual patients was used for germline exome sequencing. Extracted tumor-genomic DNA was fragmented into 300–350 bp by sonication (Covaris, Woburn, MA). Sequencing libraries were prepared with a KAPA Hyper Prep kit (KAPA Biosystems) with optimized protocols. Libraries were then subjected to PCR amplification and purification before targeted enrichment. The enriched libraries were sequenced on Hiseq 4000 NGS platforms (Illumina) to cover depths at 200x.

Further analysis of RNA-seq data

Differential gene expression analysis was conducted using DEseq2 between pre-treatment and post-treatment samples. Gene Ontology (GO) enrichment analysis of differentially expressed genes (DEGs; adjusted p-value < 0.05) was subsequently performed by using the R package, “clusterProfiler”.[18] Gene set enrichment analysis (GSEA) was then performed on all pre-ranked genes, and enrichment scores were generated for Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways. GO terms and GSEA enrichment scores with false discovery rates (FDRs) less than 0.05 were considered to be significantly enriched. Single-sample GSEA (ssGSEA) was used to generate infiltration scores for 16 immune cell types and cytotoxic cells for each tumor based on previously published gene sets[19] using the R package, “gsva”.[20] Immune scores, microenvironment scores, and stroma scores were processed using xCell.[21]

Data processing of weighted-gene correlation network analysis (WGCNA) followed the Horvath Lab UCLA protocol (https://horvath.genetics.ucla.edu/html/CoexpressionNetwork/Rpackages/WGCNA). Then, the STRING network of genes in each color module was imported into Cytoscape 3.7, and hub genes were identified by the bottleneck method provided in the cytoHubba plugin.

Tumor purity estimation

ESTIMATE scores were calculated via the “ESTIMATE” R package, and the purity of each sample was evaluated using the formula described by Yoshihara et al.[22]

Tumor-neoantigen estimation

We used a QBRC pipeline for neoantigen-calling with parameters that were set according to recommended values (https://github.com/tianshilu/QBRC-Neoantigen-Pipeline).

Statistical analysis

All statistical tests were performed using R software, version 3.6.3. Data are expressed as mean ± standard deviation (SD). Fisher’s exact tests were used for categorical data when appropriate. For comparisons of two groups, statistical significance was estimated via Student’s t-tests. Pearson’s correlation and Spearman’s rank correlation coefficient was used to evaluate the correlation between two parameters when appropriate. Log-rank tests, Kaplan–Meier curves, and smooth HR curves for survival analysis were performed using the R package, “survival”. We used the “surv_cutpoint” function in the “survminer” package to define the cutoff value of NEFH and IAS. All heatmaps were generated by the “pheatmap” function (https://github.com/raivokolde/pheatmap). OncoPrint was used to visualize multiple genomic mutations and was constructed by the package, “maftools” (https://github.com/PoisonAlien/maftools). All p-values were based on a two-sided hypothesis, and a p-value of < 0.05 was considered to be statistically significant.

Results

Transcriptional landscape of EGFR-mutant and ALK-rearranged tissue responses to TKI treatments

We enrolled 23 advanced lung ADC patients with EGFR mutations (EGFR-m, n = 15) or ALK rearrangements (ALK-r, n = 8) who received EGFR- or ALK-TKI treatments (Supplemental Table S1). Paired analysis of pre-TKI treatment and partial-response samples was performed to explore transcriptomic changes. DEGs were screened with a criterion of an FDR ≤ 5%. We identified a total of 2,867 DEGs (1,442 up-regulated and 1,425 down-regulated DEGs; Figure 1(a,b); Supplementary Table S2) that were shared in EGFR-m and ALK-r patients in response to TKI treatments. GO enrichment analysis of shared up-regulated DEGs indicated that their potential functions were significantly related to immune responses, while down-regulated DEGs were enriched in processes regulating cell-cell adhesion and junctions (Figure 1(c)). Further analysis of unique DEGs in EGFR-m (2,787 DEGs) or ALK-r (1,538 DEGs) samples show that, compared with those of ALK-r samples, there were more DEGs related to biological processes involving immunity in EGFR-m samples (Figure 1(a,b,d,e)). Moreover, as determined by GSEA of KEGG-pathway gene sets, some immune-related and cell-cycle-related pathways were found to be significantly regulated by TKI treatments (Figure 1(f,g)). Collectively, the above results suggest that TKI treatments may change the tumor microenvironment of TKI responders with EGFR mutations or ALK rearrangements.

Figure 1.

The transcriptional landscape of NSCLCs with EGFR mutations or ALK rearrangements in response to TKI treatments

TKI treatments remodel the TIM of samples with EGFR mutations or ALK rearrangements

To further verify the impact of TKI treatments on the TIM, general immune estimation scores were generated via a gene signature-based method, xCell, which converts gene expression profiles to enrichment scores.[21] In the EGFR-m samples, immune scores, microenvironment scores, and stroma scores were significantly increased after TKI treatment (Figure 2(a)). Similarly, immune scores and microenvironment scores were significantly increased in ALK-r post-treatment samples (Figure 2(b)). However, no significant change was found between pre-treatment and post-progression samples (Supplementary Figure 1(a)).

Figure 2.

TKI treatment remodels the tumor immune microenvironment

Next, we performed ssGSEA to estimate the infiltration levels of 16 cell types in samples before and after TKI treatments. In addition, infiltration of cytotoxic cells was computed by genes overexpressed in CD8 + T cells, gamma-delta T cells, and natural killer (NK) cells to represent the overall anti-tumor activity,[19] which was dramatically increased in response to TKI treatments (Figure 2(c)) but the difference was not statistically significant in TKI-resistant samples (Supplementary Figure 1(b)). Further exploration revealed that the enhanced anti-tumor response in EGFR-m was mainly due to the increased infiltration of CD8 + T cells, B cells, and so on; in ALK-r samples, NK CD56dim cells – characterized by their nonspecific-antigen-dependent cell-killing ability[23] – may have had a greater impact on the enhanced anti-tumor response (Figure 2(c)). Moreover, the expression levels of gamma-interferon-signaling and immune-checkpoint genes (e.g., STAT5B, PDCD1, BTLA, and CD27) were up-regulated in response to TKI treatments (Figure 2(d)). However, some genes related to antigen presentation were down-regulated (e.g., CALR, CANX, and PDIA3) (Figure 2(d)), which may explain the limited T-cell infiltration in responders after ALK-TKIs treatment.[24] Findings from many recent preclinical and clinical studies have supported the long-standing hypothesis that tumors induce adaptive immune responses and that the antigens that drive effective T-cell response are tumor-specific mutant peptides (neoantigens), which are generated from somatically mutated genes.[25] Therefore, RNA-seq and WES were performed to estimate the burden of tumor neoantigens and mutations. Consistent with our previous results, the reduced tendency of neoantigens (Figure 2(e)) and mutational load (Figure 2(f)) may have induced weakened antigen presentation and subsequent limited anti-tumor responses of T cells in ALK-r samples during TKI treatments. Collectively, these findings suggest that the anti-tumor responses of immune cells in EGFR-m and ALK-r responders were enhanced after TKI treatments. EGFR-m samples were dominated by adaptive immune responses, while changes in innate immune responses of ALK-r samples may be more critical under the pressure of TKI treatments.

Somatic mutational landscape between pre-treatment and post-treatment samples

TKI resistance in patients with EGFR-m or ALK-r is mainly caused by acquired mutational changes, such as the EGFR T790M mutation[1] or ALK E1384K mutation;[26] however, somatic mutational changes in post-responding samples have still not been widely recognized. To this end, we investigated somatic mutations in samples before and after TKI treatments and identified that neurofilament heavy (NEFH) was mutated in multiple individuals in response to TKI treatments (Figure 3(a–c)). These two variants (p.E658_E659del and p.E658_K665del) led to the in-frame deletion of two or eight amino acids in the motif characterized by repeats of the tripeptide K-S-P, which may affect the phosphorylation and function of NEFH.[27] A previous study established that NEFH is a cancer driver because of its functional impact on KRAS and MTORC1 signaling.[28] We conducted survival analyses of LUAD patients with the EGFR 19Del/L858R mutation or EML4-ALK fusion in the TCGA dataset; our findings suggested that patients with high NEFH expression had worse prognoses than patients with low NEFH expression (Figure 3(d,e)). Survival analyses of patients with EGFR or ALK mutations in the Pan-Cancer dataset confirmed that NEFH was a cancer-promoting factor (Supplementary Figure 2(a,b)). Consistent with our hypothesis NEFH mutants with functional limitations with EGFR or ALK mutations in the Pan-Cancer dataset showed better progression-free survival and overall survival than NEFH wildtype patients (Supplementary Figure 2(c,d)).

Figure 3.

Somatic mutational landscape between pre-treatment and post-treatment samples

NEFH mutations are associated with reduced neutrophil infiltration

To investigate the relationship between immunity and NEFH mutations, we calculated the relative infiltration scores of immune cells by subtracting baseline scores before TKI treatment from the scores after TKI treatments. We found that there were lower abundances of neutrophils in NEFH-mutant than in NEFH-wildtype post-treatment samples (P = .011; Figure 3(f)), which may aid in inhibiting tumor progression.[29] Moreover, Spearman’s correlation analysis yielded a positive correlation between NEFH expression and neutrophil infiltration in the TCGA LUAD samples with the EGFR 19Del/L858R mutation or EML4-ALK fusion (P = .056, Rho = 0.29; Figure 3(g)). Further differential gene expression analysis revealed that the expression of genes related to neutrophil chemotaxis, such as CXCL2 and CXCL8,[30] was significantly reduced in NEFH-mutant samples, compared with matched samples before treatment (Figure 3(h)). However, the causal effect of NEFH mutation on neutrophil infiltration needs more research.

Taken together, these findings demonstrate that functionally restricted mutations of NEFH were enriched in samples after rapid TKI-induced responses and may have led to reduced neutrophil infiltration and a better prognosis. However, no significant changes in relative levels of cytotoxic cells were found, which may indicate that NEFH mutations are not the main cause of increased cytotoxic cell infiltration after TKI treatments (Figure 3(f)).

Identification of hub genes involved in TKI-induced immune remodeling via WGCNA

To further explore the mechanisms of TKI-induced cytotoxic cell infiltration, we performed WGCNA to identify relevant intramodular hub genes,[31] using the top-5000 variation genes in RNA-seq samples to build the co-expression network. This analysis identified nine distinct co-expression modules that corresponded to clusters of correlated genes (Figure 4(a)). Moreover, we analyzed the relationships between the module eigengene and sample clinical characteristics (i.e., pre- vs. post-TKI treatments, immune scores, stroma scores, microenvironment scores, and cytotoxic cell-infiltration scores) to identify co-expression modules associated with the above sample traits. We found that the blue module was strongly positively correlated with TKI treatments, immune scores, microenvironment scores, and cytotoxic cell-infiltration scores, while the turquoise module was negatively associated with TKI treatments and immune scores (Figure 4(b)). Therefore, we focused on the blue and turquoise modules and determined the top-10 genes from each module (Figure 4(c,d)) via the bottleneck method.[32] The immune-associated score (IAS) (i.e., mean expression of the hub genes in the blue module minus the mean expression of the hub genes in the turquoise module) was generated (Figure 4(e)). Furthermore, Spearman’s correlation analysis confirmed that the immune-associated score had a strong correlation with cytotoxic-cell infiltration levels, microenvironment scores (Figure 4(f)), and many immune-cell types (Supplementary Figure 3).

Figure 4.

WGCNA-identified hub genes involved in TKI-induced immune remodeling

Immune-associated score-based subtypes are related to immune infiltration, immune activity, and patient prognosis in the TCGA and GSE11969 cohorts

Next, we validated these hub genes and the IAS in samples with EGFR mutations or ALK fusion in the TCGA LUAD cohort. We found that the IAS was correlated with cytotoxic cell infiltration and immune scores (Figure 5(a,b)). Further analysis confirmed that there was a strong positive correlation between the IAS and many immune cell types (e.g., activated CD8 + T cells) in the TIM (Supplementary Figure 4). Previous studies have demonstrated that there is a direct link between an activated immune microenvironment and favorable clinical outcomes in patients with various forms of cancer.[33,34] To evaluate the prognostic value of the IAS, we performed Kaplan–Meier curves and smooth hazard ratio (HR) curves of progression-free survival (PFS) and overall survival (OS) of EGFR/ALK-positive patients in the TCGA LUAD cohort. Kaplan–Meier curves showed that under the optimal cutoff point, the high-score group had a higher PFS and OS compared to those of the low-score group (Figure 5(c,d)). Similarly, the analysis of patients with EGFR mutations in the GSE11969 dataset showed that the IAS was significantly associated with the immune remodeling process and patient OS (Figure 5(e,f)).

Figure 5.

Immune-associated score-based subtypes are related to immune infiltration, immune activity, and patient prognosis in the TCGA and GSE11969 cohorts

Discussion

NSCLCs with EGFR mutations or ALK rearrangements are generally inactive tumors with a low tumor mutational burden (TMB)[35] and a lack of T-cell infiltration.[4] Immunotherapy with immune-checkpoint inhibitors has been approved as a treatment for NSCLC.[36] It is unlikely that any single immunotherapy will be capable of changing the outcomes of NSCLCs harboring EGFR mutations or ALK rearrangements.[6,37,38] Although there have been many clinical trials attempting to explore a combination of TKIs and immunotherapy, the results have shown that the clinical benefit is far lower than expected.[5,39] At present, the data supporting the front-line use of these drugs in EGFR-mutant PD-L1-positive NSCLC patients have been encouraging, but few studies have provided conclusive evidence supporting the combination of TKIs and immunotherapy as front-line treatment for these patients.[40] Since immune status has a key impact on the efficacy of immunotherapy, further elucidation of the effects of TKI treatments on the immune microenvironment is needed. Unfortunately, previous studies have shown that after the development of TKI resistance, TKI cannot effectively improve T-cell infiltration,[14,15] and there is a lack of research on the remodeling of the TIM in post-TKI-responsive samples.

In the present study, compared with baseline parameters, TKI treatments significantly improved immune-cell infiltration and cytotoxicity in TKI-responsive samples but not in TKI-resistant samples, which suggests that the combined application of immune-checkpoint inhibitors may be more effective before TKI resistance. Previous studies have shown that the lower expression level of PD-L1 in EGFR-mutant NSCLCs results in poor inferior response to PD-1/L1 blockade.[4,7] In our study, TKI treatment did not increase the expression of CD274 (the gene encoding PD-L1 protein), suggesting that the combination of TKI and anti-PD-L1 may not have a significant synergistic effect on tumor treatment. However, a variety of genes associated with IFN-gamma signaling and immune-checkpoint processes were up-regulated in response to TKI treatments. Further analysis revealed that PDCD1 (the gene encoding PD-1 protein) and other inhibitory-checkpoint genes (e.g., BTLA) were up-regulated following TKI treatments. This finding may explain the poor effects of PD-1 inhibitors combined with TKIs (immune suppression of multiple immune checkpoints) and may provide evidence for targeting subsequent combination therapies. Furthermore, in ALK post-response samples in the present study, the reduction of antigen-presentation-related genes (e.g., CALR, CANX, and PDIA3), tumor neoantigens, and mutational load may have limited the activation of antigen-specific immune cells such as CD8 + T cells.

Since changes in the genome may also affect immune-cell infiltration, we determined that oncogene NEFH mutations that may have restricted functions were enriched in TKI-treated samples, as determined via WES analysis. We found that NEFH mutations were related to a good prognosis and decreased neutrophils, but were not related to increased cytotoxicity. The complexity of the genomic landscape of EGFR-mutant LUADs was unexpected, and many co-driver gene mutations such as TP53 were related to a poor prognosis.[41] We tried to determine the relationship between TP53 mutation and the changes in TIM after TKI treatment. Surprisingly, we found that the increase in cytotoxicity of EGFR-mutant LUADs with TP53 co-mutation was more significant than that of patients with TP53 wildtype in response to TKI (Supplementary Figure 5(a)), which may have resulted from the loss of immunosuppressive TP53 co-mutation after TKI treatment.[42] However, the more obvious cytotoxicity up-regulation was found in ALK-rearranged LUAD with TP53 wildtype (Supplementary Figure 5(b)). These findings revealed that TP53 mutation may be a predictor of the effectiveness of TKI-immunotherapy combinations for NSCLCs, but further studies are needed.

Next, we identified hub genes related to the TKI-induced increase in immune infiltration via WGCNA and the bottleneck method, from which we determined the IAS. Interestingly, we found that the IAS composed of the hub genes driven by TKI-induced immune remodeling was positively correlated with immune infiltration. Based on the IAS, we divided the samples of the public dataset into high-score and low-score groups, and these two groups of samples showed differential prognoses and immune infiltration patterns, the findings of which deserve further investigation.

Taken together, we revealed the complexity of TKI-induced immune remodeling in TKI-responsive samples with EGFR mutations or ALK rearrangements, in which changes in immune-related mRNA and NEFH mutations may also play an important role. Nevertheless, the present study has a limitation of not fully considering the issue of intratumoral heterogeneity, which was reported to be a potential cause for the difference in response to TKI.[41,43] Moreover, the relationship between TIM and TKI should be validated by using a large-sized sample and some other in vitro or in vivo experimental studies. Despite these limitations, the present study provides clues for combination therapy of TKI and immunotherapy.

Click here for additional data file.