The Multi-Omics Analysis of Key Genes Regulating EGFR-TKI Resistance, Immune Infiltration, SCLC Transformation in EGFR-Mutant NSCLC

Introduction

Lung cancer is one of the most common malignant tumors in the world, and approximately 80% of lung cancer is NSCLC.¹ Although the efficacy of treatment for advanced NSCLC is still not satisfactory, the prognosis of NSCLC increased gradually during recent years because of the emerging new molecular targeting medicines.² The successful treatment of EGFR mutation by EGFR-TKI is almost the milestone of precision medicine in NSCLC. It could be used as a first-line treatment for advanced NSCLC harboring TKI-sensitive EGFR mutations. TKI is superior to chemotherapy in terms of progression-free survival and response rate. When combined, TKI and chemotherapy might further improve overall survival.³ However, patients would suffer from drug resistance after taking TKI for about 9–14 months.^4,5 Although some resistance-related alterations have already been found, such as T790M of 20th exon,⁶ MET amplification,⁷ transforming to SCLC,⁸ and so on, the resistance mechanisms of the 1st or 2nd generation of EGFR-TKI (gefitinib, afatinib, etc.) are still not completely clear. Therefore, it is urgently needed to estimate some new biomarkers to predict EGFR-TKIs resistance, and facilitate the best choice of treatment in clinical practice.

Cancer-associated immune cell infiltration in the tumor microenvironment could lead to tumor growth, metastasis, and therapy resistance. The immune and inflammatory cells include lymphocytes, neutrophils, monocytes, tumor-associated macrophages, and so on. EGFR-TKIs could regulate tumor-infiltrating immune cells of EGFR-mutant NSCLC, and tumor-infiltrating CD8+ lymphocytes significantly decreased after treatment.⁹ And such changes might associate with EGFR-TKI resistance, and provide clues for subsequent treatment like immune checkpoint inhibitors.

On the other hand, although NSCLC and SCLC have different pathological and genome characteristics, they might have a common cell source. It accounts for approximately 3–15% of resistance related to NSCLC transforming to SCLC.¹⁰ RB1 and TP53 loss,¹⁰ neuron-specific enolase (NSE),¹¹ and plasma gastrin-releasing peptide (ProGRP)¹² could be early predictive markers of NSCLC-to-SCLC transformation. And there might be other molecules driving the transformation.

Therefore, this study aims to explore whether there might be some DEGs linked to the resistance to both the 1st and 2nd generation of EGFR-TKIs, immune cell infiltration, and SCLC transformation by bioinformatics analysis and clinical checking samples. Moreover, we might find some mechanisms, and provide a new insight to overcome the resistance to EGFR-TKI. A study workflow is presented in Figure 1.

Figure 1 Study workflow.

Materials and Methods

Access Data

Gene Expression Omnibus (GEO: www.ncbi.nlm.nih.gov/geo) belongs to the National Center for Biotechnology Information (NCBI), an affiliation of the National Institutes of Health (NIH).¹³ And GEO is one of the largest public gene expression data resources now. We selected 3 profiles of gefitinib, erlotinib, and afatinib from the GEO database for analysis: GSE122005¹⁴ based on GPL570 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array, GSE31625¹⁵ based on GPL96 [HG-U133A] Affymetrix Human Genome U133A Array, and GSE62504¹⁶ based on GPL6947 Illumina HumanHT-12 V3.0 expression beadchip. GSE122005 (gefitinib) contained 3 resistant samples and 3 sensitive samples, GSE31625 (erlotinib) contained 28 samples and 18 sensitive samples, and GSE62504 (afatinib) contained 2 resistant samples and 1 sensitive samples.

Identification of DEGs

Firstly, we used GEO2R (GEO online tool) to search DEGs between EGFR-TKI-resistant tissue and sensitive tissue samples. The cut-off criteria for identifying DEGs were P<0.05 with the log FC ≥1 or ≤−1. And we used the Venn diagram to find common DEGs.

Functional Annotation of DEGs

Gene ontology (GO) covers three aspects: biological processes, molecular functions, and cellular components. Kyoto Encyclopedia of Genes and Genomes (KEGG) is a database resource for understanding high-level functions and utilities of the biological system. GO analysis of DEGs was performed by the DAVID v6.8 online tool (https://david.ncifcrf.gov/).¹⁷ Moreover, we used KEGG rest API (https://www.kegg.jp/kegg/rest/keggapi.html) to obtain the latest KEGG pathway gene annotation and used R package clusterProfiler for enrichment analysis.

Construction of PPI Network and Hub-Gene Extraction

A network could be built to show relationships among proteins by STRING database (https://string-db.org/).¹⁸ We used it to construct a PPI (protein–protein interaction) network for DEGs, used Cytoscape to visualize protein network, and extracted key modules in PPI network by MCODE.¹⁹ The gene that meets degree over 10 would be defined as the hub gene.

Hub Genes Expression Validation

We performed some expression analyses of hub genes by GEPIA (http://gepia.cancer-pku.cn/index.html).²⁰ To further confirm results, we analyzed the protein expression of hub genes in normal lung and lung adenocarcinoma tissue microarrays by The Human Protein Atlas (http://www.proteinatlas.org/).²¹ The tissue microarray was stained by antibody HPA027541. LUAD patients’ tumor RNA-seq expression profiles, paired normal tissue samples, and clinical information were downloaded from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/). And EGFR wild or mutant samples were determined according to the prediction results of mutect2 software. Statistical analyses were performed by using R-4.0.4, and figures were implemented by R package ggplot2.

Prognostic Survival Analyses

We used GEPIA to perform survival analysis of hub genes. Furthermore, EGFR-mutant LUAD patients were divided into four groups according to the expression level of SPP1. Patients’ survival status and expression heatmap were drawn by R package ggrisk. The DSS (disease-specific survival) Kaplan–Meier analysis with log rank test was used to compare differences of survival curves by R package survival and survminer. Time-dependent ROC analysis was performed to compare the predictive accuracy of SPP1 by R package timeROC.

Analysis of Immune Cell Infiltration

Analysis was performed in TIMER (https://cistrome.shinyapps.io/timer/)²² of the correlations of SPP1 and immune cells infiltration level in LUAD. We divided LUAD patients in TCGA into EGFR mutant/wild-type and SPP1 high/low expression (median cutoff) groups. We utilized R package immunedeconv to find differential immune infiltrations between groups by integrating xCell algorithm.²³ And we used ssGSEA to calculate the immune infiltration score for 28 kinds of immune cells of erlotinib-resistant and sensitive samples in GSE31625.²⁴ The methods of anti-tumor immunity (cells regulate anti-cancer immune responses) and pro-tumor suppression (cells regulate immune suppressive functions) were referred from Jia et al’s recent publication.²⁵ The relationship of erlotinib-resistant and sensitive samples’ anti-tumor immunity and pro-tumor suppression was estimated by Pearson’s correlation analysis. To better understand the relationship among various immune cells, we also calculated Pearson’s correlation coefficients among 28 immune cells to build the network. Finally, a box-plot was constructed to compare the difference of 28 immune cells between erlotinib-resistant and sensitive subgroups. All above immune analyses of GSE31625 were achieved by using R package GSVA.²⁶

Kaplan–Meier Plotter Database Analysis

We performed survival analysis to estimate the relationship between immune cell infiltration and the expression level of SPP1 by KM Plotter (http://kmplot.com).²⁷ LUAD patients were divided into decreased and enriched subgroups according to the difference of immune cell infiltration, and Kaplan–Meier survival analysis was used to evaluate the effect of SPP1 on patients’ prognosis. Statistical significance was defined as P value <0.05.

Co-Expression of SPP1, Relationship Between SPP1 and TMB/MSI

To further understand biological significance of SPP1 in LUAD, we performed the enrichment analysis of co-expressed genes of SPP1 in KEGG and Reactome pathway by GSEA in LinkedOmics (http://www.linkedomics.org/admin.php).²⁸ We specifically estimated the expression relationship between SPP1 and PD-L1 (CD274). RNA-seq data of 33 kinds of tumor were downloaded from TCGA GDC data portal website (https://portal.gdc.cancer.gov/). For the methods of calculating tumor mutational burden (TMB) and microsatellite instability (MSI) we refer to 2 publications.^29,30 R 4.0.4 was used for statistical analysis with a significant difference determined at P value <0.05.

Hub Genes of SCLC, TKIs Resistance and miRNA Network

We found SCLC profile to obtain DEGs between SCLC and normal tissues in GEO database: SCLC profile GSE149507 based on GPL23270 [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array, and GSE149507 contains 18 pairs of SCLC and adjacent lung tissues. And we used a Venn diagram to find co-expressed genes between SCLC and previous TKI resistance profiles. Thus, we might have genes related to SCLC transformation of TKI resistance. Additionally, we downloaded SCLC and NSCLC microRNA (miRNA) list, and got the common miRNAs of these two diseases from HMDD (the Human microRNA Disease Database) (http://www.cuilab.cn/hmdd/), which is a database that curated experiment-supported evidence for human miRNA and disease associations.³¹ And we predicted the interaction between transformation-related genes and common miRNAs by miRWalk (http://mirwalk.umm.uni-heidelberg.de/).³²

Methylation and Prognostic Model of Hub Genes

The expression and DNA methylation of SPP1 in LUAD were analyzed by DNMIVD: DNA methylation interactive visualization database (http://119.3.41.228/dnmivd/index/). We used the prognostic model in DNMIVD to get CpGs related to 5 hub genes, and construct univariate, multivariate proportional hazards regression model. Kaplan–Meier plot was drawn by dividing patients into high-risk or low-risk groups based on the results of multiple proportional hazards regression model.³³

Drug Sensitivity Analysis of SPP1

We downloaded data (Compound activity: DTP NCI-60, RNA: RNA-seq) from Cell Miner (https://discover.nci.nih.gov/cellminer/home.do).³⁴ The NCI-60 is a panel of 60 diverse human cancer cell lines used to screen over 100,000 chemical compounds and natural products. Drug activity levels expressed as 50% growth-inhibitory levels (GI50) are determined by the Developmental Therapeutics Program (DTP). It defined the drug activity Z scores above 0 or below 0 indicating that cells are sensitive or resistant to some drugs, respectively. Cell Miner data analysis was performed by using R-4.0.4 to draw the scatter plot between SPP1 expression level and drug activity Z scores, and each one point means one of 60 various human cancer cell lines.

Real-Time Quantitative PCR

Total RNA was extracted and isolated from thirteen EGFR-mutant (19Del/21L858R) LUAD samples, eight of them associated with paired paracancerous tissue samples, by using the TRIzol method. Afterwards, RNA was reverse-transcribed into complementary DNA (cDNA) using the PrimerScript RT reagent kit (TaKaRa). The transcription level of SPP1 was measured by qRT-PCR according to the assay protocol (Bio-Rad CFX ConnectTM Real-Time Detection System). Glyceraldehyde-phosphate dehydrogenase (GAPDH) was used as internal control for mRNA normalization, and relative expression of target gene was quantified by using the 2-ΔΔCT method. The reference sample was the average of all paracancerous tissue samples. The unpaired t-test was used for intergroup comparisons. Results were considered statistically significant at P value <0.05.

Results

Identification of DEGs of TKI Resistance

In GEO database, we found three EGFR-TKI-resistant versus sensitive profiles, and DEGs of three profiles were analyzed by GEO2R. A Venn diagram was employed to get common DEGs, and in total 107 DEGs are identified in Figure 2A.

Figure 2 (A) Venn diagram, common DEGs of GSE122005, GSE31625, GSE62504. (B) a. GO enrichment; b. KEGG pathway annotation. (C) PPI network of DEGs was constructed in Cytoscape; the greater the degree of correlation between genes, the larger the circle and darker the colors. Abbreviations: DEGs, differentially expressed genes; GO, gene ontology; KEGG, Kyoto Encyclopedia of Genes and Genomes; BP, biological processes; CC, cellular components; MF, molecular functions; PPI, protein–protein interaction.

GO and KEGG Pathway Enrichment of DEGs

The results showed that DEGs were mainly enriched in cell adhesion, extracellular matrix organization, positive regulation of GTPase activity, negative regulation of cell proliferation, and cell−cell signaling for biological processes. Cellular components were significantly enriched in extracellular exosome, plasma membrane, extracellular space, extracellular region, and integral components of plasma membrane, and molecular functions were significantly enriched in protein binding (Figure 2B-a). And KEGG pathway analysis of DEGs showed that they were associated with cell adhesion molecules, focal adhesion, microRNAs in cancer, PI3K-Akt signaling, T cell receptor signaling, and so on (Figure 2B-b).

DEGs PPI Network

PPI network of DEGs was constructed by STRING database (Supplementary Figure 1), and visualized in Cytoscape. MCODE extracted four key modules (Supplementary Figure 2); the greater the degree of correlation among genes, the larger the circle and darker the colors (Figure 2C). Finally, we found 11 highly connected genes (hub genes): RAB25, CD44, SPP1, ITGB4, ESRP1, FYN, CDH1, ERBB2, CLDN7, ZEB1, and MUC1. The expression of ZEB1 and SPP1 was up-regulated in three profiles. On the other hand, RAB25, ERBB2, ESRP1, ITGB4, MUC1, CLDN7, and CDH1 were down-regulated in three profiles. CD44 and FYN were both decreased in erlotinib- and gefitinib-resistant profiles and increased in afatinib-resistant profile.

Hub Gene Expression and Survival Analyses

Next, we analyzed the expression level of hub genes on GEPIA. The expression level of RAB25, ERBB2, SPP1, CDH1, ESRP1, and ITGB4 was significantly higher, and that of of ZEB1 and FYN was significantly lower in lung adenocarcinoma compared to normal tissue, respectively (Figure 3A). And lower expression of SPP1 and ITGB4 related to better overall survival both in LUAD and LUSC patients compared to higher expression. On the other hand, the simultaneously high expressions of SPP1 and ZEB1 were associated with poor overall survival in LUAD and LUSC patients compared to low expression (Figure 3B). Therefore, we thought that SPP1 was of great significance in the survival of NSCLC.

Figure 3 (A) The expression level of hub genes in lung adenocarcinoma and normal tissues in GEPIA (*P < 0.01). (B) Survival analysis. Expression level of SPP1 and ITGB4 was related to the overall survival of LUAD and LUSC patients (P < 0.05), expression level of SPP1+ZEB1 was related to overall survival of LUAD and LUSC patients (P < 0.05). Abbreviation: GEPIA, gene expression profiling interactive analysis.

The Expression of SPP1 in Pan-Cancers and Prognostic Potential

Thus, we investigated the expression of SPP1 in pan-cancers. The expression of SPP1 was significantly higher in multiple kinds of tumor tissues compared with their cognate normal tissues (Figure 4A). Based on the staining of SPP1 protein in tissue microarray, SPP1 was mainly located in cytoplasm. And the staining was moderate in lung adenocarcinoma cells, and low in normal lung tissue cells according to the results of HPA database (Figure 4B). And we compared SPP1 expression in EGFR-mutant, wild-type of LUAD, and normal samples from TCGA. The expression level of SPP1 was significantly higher in the EGFR-mutant LUAD compared to the other two types (Figure 4C). According to the expression level of SPP1, we divided EGFR-mutant LUAD patients into four groups (Figure 5A). Among them, the lowest expression group had the longest DSS. The median survival time of 25–50% and 50–75% group was 3.8 years and 3.2 years, respectively (Figure 5B). And the AUC of the 3-year ROC is 0.762 (Figure 5C). Therefore, the expression level of SPP1 might be a prognostic factor in EGFR-mutant lung adenocarcinoma.

Figure 4 (A) Validation of the SPP1 expression in multiple tumors. (B) Validation of the SPP1 expression by the Human Protein Atlas database. a. SPP1 of immunohistochemistry in normal lung tissue; b. SPP1 of immunohistochemistry of LUAD tissue. (C) SPP1 expression in EGFR-mutant/wild LUAD and normal groups from TCGA. P value significant codes: *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001.

Figure 5 (A) Upper panel: The scatter plot of SPP1 expression from low to high, different colors represent different expression groups; Middle panel: The survival time and survival status of patients; Lower panel: Heatmap of SPP1 in EGFR-mutant LUAD. The abscissas of the upper, middle, and lower panels all represent same samples. (B) Kaplan–Meier survival analysis of SPP1. (C) Time-dependent ROC analysis of SPP1.

Immune Landscape of EGFR-Mutant/Wild and TKI-Resistant Groups

We further analyzed the degree of immune cell infiltration related to prognostic hub gene expression of lung adenocarcinoma by TIMER database. There was a positive correlation between the expression of SPP1 and the infiltration of macrophage, neutrophil, and dendritic cells (Figure 6A). xCell results showed that the abundance of CD8+ T cells and CD4+ Th2 was lower in LUAD patients with EGFR mutation compared to wild-type (Figure 6B). There were more macrophage M1 and CD4+ Th2 cells, and fewer CD8+ T cells in LUAD patients with high expression of SPP1 compared to low expression (Figure 6C). These results might suggest that SPP1 could regulate tumor immune tolerance and immune escape. Immune infiltration landscape is shown in Figure 7A. Two categories of immunity ability (anti-tumor immunity and pro-tumor suppression) were positively associated (R=0.5438, P<0.001) in Figure 7B. The correlation heatmap (Figure 7C) revealed that effector memory CD4 T cells, NK T cells, MDSC (myeloid-derived suppressor cells), and Th2 T cells were moderately or highly correlated. We compared the differences between each kind of immune cells and erlotinib-sensitive and resistant samples in Figure 7D. We found fewer CD8+ T cells and natural killer cells, and more CD4+ T cells in the resistant group compared to the sensitive group.

Figure 6 (A) The relationship of SPP1 expression level and tumor immune infiltration level. (B) The immune scores distribution of EGFR-mutant and wild-type in LUAD. Blue color represents EGFR-mutant group, yellow represents EGFR-wild group. (C) The immune scores distribution of SPP1 high and low expression in LUAD. Blue color represents SPP1 high group, yellow represents SPP1 low group. P value significant codes: *P < 0.05, **P < 0.01, ****P < 0.0001.

Figure 7 Immune landscape of GSE31625. (A) Heatmap of two subtypes from the GSE31625 based on ssGSEA scores calculated from 28 immune cells. (B) Correlation between anti-tumor immunity and pro-tumor immune suppressive functions; R: coefficient of Pearson’s correlation; the shaded area represents 95% confidence interval. (C) Correlation of 28 immune cells. (D) The relative immune infiltration scores of 28 immune cells between erlotinib-sensitive and resistant groups. P value significant codes: *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001.

Prognostic Analysis of SPP1 with Immune Cells in LUAD

We performed prognosis analyses based on the expression of SPP1 in related immune cell-enriched or depleted LUAD subgroups (Figure 8A). LUAD patients with high expression of SPP1 and decreased infiltration of basophils, CD4+ memory T cells, CD8+ memory T-cells, and regulatory T-cells had poor prognosis. And LUAD patients enriched in eosinophils, mesenchymal stem cells, natural killer T-cells, and type 1 T-helper cells with high expression of SPP1 had poor prognosis too. What is more, LUAD patients with high expression of SPP1 had significantly poorer prognosis with whatever infiltration level of B-cells and macrophages. These results indicated that SPP1 might be an effective prognosis marker of LUAD with different immune infiltration.

Figure 8 (A) A forest plot showed the prognostic value of SPP1 expression according to different immune cell subgroups in LUAD patients. (B) Correlation between PD-L1 and SPP1. (C) Correlation analysis of TMB and SPP1. (D) Correlation analysis of MSI and SPP1. The horizontal axis in the figure represents the correlation coefficient; the ordinate is different tumors. The size of dot represents the correlation coefficient, and colors represent the significance of P value. (E) A Venn diagram of SCLC and three TKI-resistant profiles. (F) RNA-miRNA network of five hub genes. Abbreviations: TMB, tumor mutational burden; MSI, microsatellite instability; SCLC, small cell lung cancer.

The Co-Expression of SPP1 and TMB/MSI in Pan-Cancer

The positive correlation between PD-L1 (CD274) and SPP1 is explored in Figure 8B. The KEGG pathway of SPP1 co-expressed genes mainly promoted tumor cell proliferation in osteoclast differentiation, ECM-receptor interaction, HIF-1 signaling pathway, glycolysis, DNA replication, NF-kappa B signaling pathway, and cytokine–cytokine receptor interaction (Supplementary Figure 3A). And Reactome pathway analysis revealed that co-expressed genes promoted extracellular matrix organization, immunoregulatory interactions between a lymphoid and a non-lymphoid cell, neutrophil degranulation, synthesis and replication of DNA, and so on, shown in Supplementary Figure 3B. These enrichment results showed that SPP1 expression is related to tumor cell proliferation and immunity. Moreover, we found that SPP1 was not only closely related to TMB and MSI in LUAD but also related to various tumors by the analysis of SPP1 expression and TMB/MSI. The expression of SPP1 was associated with high TMB in DLBC, THYM, SARC, ACC, and PRAD (Figure 8C). On the other hand, SPP1 was negatively correlated to MSI in LUAD and LUSC. High correlation of SPP1 and MSI was observed in COAD, ACC, SARC, READ, and TGCT (Figure 8D).

SCLC Transformation Co-Expression RNA-miRNA Network

SCLC transformation is also one of the reasons for the resistance to EGFR-TKI. We compared the SCLC and TKI-resistant DEGs, and 42 genes were co-expressed in both profiles. CD44, MUC1, ESRP1, SPP1, and ZEB1 were hub genes of both SCLC and TKI resistance (Figure 8E). And we found that hsa-miR-495-3p, hsa-miR-24-3p, hsa-miR-181a-5p, and hsa-miR-125a-3p were related to both SCLC and NSCLC by using disease-miRNAs network (Supplementary Figure 4). Meanwhile, 5 hub genes were predicted to be target genes of these miRNAs (Figure 8F). A RNA-miRNA regulatory network might induce the TKI resistance by SCLC transformation.

Establishment and Evaluation of Methylation Prognostic Risk Models

SPP1 expression and DNA methylation were compared between normal and LUAD patients. The expression of SPP1 was higher, and promoter methylation level was lower in LUAD compared with normal tissue (Figure 9A). High methylation of SPP1 promoter represented decreased gene expression. A Cox proportional hazards model including CpGs of 5 hub genes (CD44, MUC1, ESRP1, SPP1, and ZEB1) was used. The univariate and multivariate analysis revealed that SPP1 CpG cg00088885 and CD44 CpG cg20971158 in LUAD were significantly associated with OS (Figure 9B and C). Z-scores of partial hazard matched to the prognostic classifier and patient survival status. We found that patients with lower Z-scores generally had better outcome than those with higher scores (Figure 9D-a). According to multivariate proportional hazards regression model, patients were divided into high-risk or low-risk groups, and the low-risk group had better survival benefit (Figure 9D-b).

Figure 9 (A) a. Boxplot of expression of SPP1 in LUAD; b. Boxplot of DNA methylation of SPP1 in LUAD. (B) Univariate Cox analysis of hub gene CpGs. (C) Multivariate Cox analysis of hub gene CpGs. (D) a. Z-score distribution of the prognostic classifier and patients’ survival status; b. Kaplan–Meier survival analysis for LUAD patients based on multivariate proportional hazards regression model. (E) The high expression of SPP1 in EGFR-mutant lung adenocarcinoma compared with paracancerous tissues by qRT-PCR.

TKIs Sensitivity Associated with SPP1

Relying on Cell Miner data, we acquired the SPP1 expression level in 60 diverse human cancer cell lines. The SPP1 expression of the NCI-60 cell lines and their Z scores for gefitinib, afatinib, and osimertinib are shown in Supplementary Figure 5. Most of cancer cell lines had Z scores below 0, meaning they are resistant to drugs with SPP1 high expression. Not only lung cancer but also other kinds of cancers such as breast cancer, glioblastoma, colon cancer, and melanoma were resistant to TKIs when SPP1 gene expressed.

The Expression Levels of SPP1 in EGFR-Mutant Lung Adenocarcinoma and Paracancerous Tissue

To examine the expression of SPP1 in EGFR-mutant LUAD, thirteen LUAD tissues and eight paracancerous tissues were collected. qRT-PCR analysis showed that SPP1 expression was significantly increased in EGFR-mutant LUAD tissues compared with adjacent normal tissues (P=0.0188) (Figure 9E).

Discussion

NSCLC patients harboring EGFR mutations developed resistance after taking TKI therapy for a period of time. The mechanisms of resistance were associated with T790M mutation, MET amplification, RET amplification, and so on. However, multiple unknown mechanisms still existed and needed further exploration. And resistance was also related to intratumoral clonal heterogeneity, dynamic changes of tumor cells, and tumor microenvironments. Although EGFR mutation could drive carcinogenesis of NSCLC, it is equally important to estimate the meaning of the copy number changes of some hub genes.³⁵ On the other hand, the relevant scores of immune infiltration induced by these hub genes might be helpful for the guidance of immunotherapy after TKI resistance.

In this study, we found that the expression level of two hub genes, ZEB1 and SPP1, increased in resistant samples compared with sensitive samples. It was interesting that ZEB1 was down-regulated in LUAD, but up-regulated in resistant samples. And the expression of SPP1 was not only higher in EGFR-mutant LUAD than that of wild-type, but also highly expressed in the 1st and 2nd generation TKI-resistant samples. Recent studies have pointed out that abnormally high SPP1 expression could activate FAK/AKT and ERK signaling pathways by up-regulating αVβ3 to promote cell proliferation for EGFR-TKI resistance in EGFR-mutant NSCLC.³⁶ Moreover, we found that SPP1 could be a prognostic factor in EGFR-mutant patients. Thus, it suggests that both ZEB1 and SPP1 might be related to resistance to EGFR-TKIs.

SPP1 encodes the protein, secreted phosphoprotein 1 (SPP1), also known as osteopontin.³⁷ Now, SPP1 has been found to abnormally express in a variety of cancers, and induces drug resistance, progression, recurrence, and metastasis in breast, ovarian, and colon cancer.^38–40 However, related resistance genes could be identified in both resistant cell lines and wild-type cells. Similarly, some kinds of TP53 mutation could be in found in NSCLC with or without EGFR mutation. The TP53 and EGFR co-mutations were associated with shorter progression-free survival for the 1st or 2nd generation EGFR-TKIs treatment compared with TP53 wild-type. It was a negative prognostic factor in advanced NSCLC.⁴¹ When EGFR mutation is coupled with high expression of SPP1, we might predict the prognosis of advanced NSCLC patients more accurately.

Interestingly, highly expressed SPP1 could induce lung cancer immune escape by up-regulating the expression of PD-L1, and regulating the polarization of tumor-associated macrophages.⁴² Although PD-L1 expression increased due to oncogenic signaling pathway activation in some EGFR-mutant lung adenocarcinomas, it resulted in unsatisfactory prognosis due to the lack of functional tumor-infiltrating lymphocytes (TILs) or “cold” TIME (tumor immune micro-environments).⁴³ This cold TIME needs the regulation of immune microenvironments to enhance the response. In EGFR-TKI-resistant patients, if SPP1 is highly expressed, it would imply high infiltration of CD4+ T cells in TIME based on our results. Thus, the TIME is transformed into PD-L1+ with more TILs. This type of TIME is also described as an “inflamed” TIME.

The Car-T treatment also introduced T cells exogenously to regulate the infiltration of TILs in TIME to achieve anti-tumor effect,⁴⁴ like the regulation of SPP1 in immune cell infiltration. This modification could strengthen patients’ response to immune checkpoint inhibitor therapy. What is more, CD4 is a reliable biomarker of ICIs, and it correlates to promoting the activity of cytotoxic T lymphocytes and establishing anti-tumor immunity. Specific dendritic cells transmitted signals from CD4+ T cells to CD8+ T cells to amplify the cytotoxic response.⁴⁵ The poor efficacy of ICIs in lung cancer patients with EGFR mutation also might be related to poor antigen presentation with low CD4+ T cells,⁴⁶ and low expression of SPP1.

In the IMPower150 trial, the patients with EGFR-sensitive mutation who were treated with atezolizumab plus bevacizumab and chemotherapy had better results in terms of efficacy than the anti-vascular combined with chemotherapy groups, with a median OS of 29.4 months vs 18.1 months (HR=0.6).⁴⁷ Vascular endothelial growth factor-targeted therapy could “normalize” tumor blood vessels, introduce more immune cells to the TIME, and regulate the tumor microenvironment to enhance the immune response.⁴⁸ For EGFR-TKI-resistant patients, Lam et al have shown that the combination of atezolizumab, bevacizumab, pemetrexed, and carboplatin achieved promising efficacy in metastatic EGFR-mutant NSCLC after TKI failure. The median PFS of patients reached 9.4 months, while the median OS did not reach and 72.5% of patients survived more than 1 year. After the treatment failure, many patients rechallenged EGFR-TKIs, and could still benefit from it, which means achieving targeted drug resensitization.⁴⁹

Therefore, when SPP1 is highly expressed, it could regulate the immune microenvironment of TKI-resistant patients like anti-vascular drugs, and such patients might benefit from the immunotherapy. The expression level of SPP1 should be a basic test for subsequent treatment in patients resistant to EGFR-TKI to determine who would be suitable for immunotherapy. Overall, SPP1 might be viewed as a new biomarker for TIME and ICIs, and might serve as a prognostic biomarker associated with immune cell infiltration in lung cancer.

Lung cancer stem cells have been proven to be one of the key factors for EGFR-TKI resistance, and TKI resistance can induce a cancer stem cell (CSC) phenotype.^50,51 Moreover, as for small cell transformation in TKI-resistant patients under the selective pressure of TKIs, multidrug-resistant cancer stem cells were retained and enriched. Thus, alveolar type II cell as the common precursor could differentiate into NSCLC firstly, and transform into SCLC under selective pressure.⁸ This transformation is due to the plasticity of the cell line, and the mechanism is not fully clear.⁵²

We obtained a regulation network by using TKI-resistant and SCLC co-expressed hub genes. Among them CD44 not only binds to SPP1 ligands,⁵³ but also functions as a recognized tumor stem cell marker.⁵⁴ And ZEB1 could induce epithelial mesenchymal transition (EMT), and ESRP1 and CD44 contribute to early pathogenesis and metastatic potential in lung cancer.⁵⁵ Some studies found that knocking down SPP1 significantly reduced stemness characteristics in pancreatic cancer cells.⁵⁶ All of these hub genes were involved in the stemness regulation. Therefore, they might be involved in the small cell transformation in TKI resistance. However, not all TKI resistance samples were related to SCLC transformation. In these three TKI-resistant GEO profiles, only the afatinib group contained the co-increasing of SPP1 and CD44, indicating that this group might have the potential of SCLC transformation. This inference would be investigated in our future work by performing experiments on transformation to SCLC.

Conclusion

In summary, the resistance to EGFR-TKI is a critical problem after the initiation of targeting therapy. Several molecular mechanisms are associated with the resistance. Our results suggested that the up-regulation of SPP1 might induce resistance to the 1st and 2nd generation EGFR-TKIs, and SPP1 might be viewed as an independent prognostic biomarker of TKI treatment. Moreover, the high expression of SPP1 might promote the anti-tumor immune cell infiltration, and ICIs therapy could be considered in this situation. At the same time, SPP1 and CD44 DNA methylation might be independent prognostic factors in LUAD patients. Moreover, 5 hub genes might promote transformation of LUAD to SCLC by regulating the tumor stem cells in EGFR-TKI resistance.