A three‑gene signature correlated with MAPK/ERK activation characterizes acquired resistance to EGFR‑tyrosine kinase inhibitors in non‑small cell lung cancer

Introduction

According to global cancer statistics data reported in 2024, lung cancer remains the leading cause of cancer-related mortality, accounting for 20% of all cancer-related mortalities worldwide (1). Lung cancer is classified into small cell lung cancer and non-small cell lung cancer (NSCLC) based on histopathological characteristics, with NSCLC accounting for 80% of lung cancer cases (2). Clinical statistics indicate that ≥70% of patients with NSCLC are diagnosed at an advanced stage (3), with the 5-year survival rate for patients with advanced NSCLC being <10% (4). The 2022 NCCN Clinical Practice Guidelines for NSCLC recommend targeted therapy as the preferred first-line treatment for advanced-stage NSCLC (5).

Epidermal growth factor receptor (EGFR) is a member of the erythroblastic leukemia viral oncogene B (ErbB) family. Upon ligand binding, the activated tyrosine kinase domain of EGFR triggers downstream signaling pathways, which subsequently promote tumor cell proliferation, differentiation and migration (6). Abnormally high EGFR expression is present in >60% of NSCLC tumor tissues (7). As a result, EGFR has become a key target for targeted therapy in advanced-stage NSCLC. Tyrosine kinase inhibitors (TKIs) have been developed and extensively applied in the treatment of advanced-stage NSCLC with EGFR activation (8). Among the EGFR-targeting TKIs, second-generation inhibitors are increasingly recognized for their pan-ErbB inhibitory activity, rendering them particularly effective in the treatment of uncommon EGFR mutations (9). Afatinib, as an irreversible inhibitor of the ErbB family, has demonstrated robust preclinical inhibitory activity against both wild-type and mutant EGFR (10). However, afatinib treatment leads to acquired resistance, leading to clinical treatment failure (11). Therefore, elucidating the molecular mechanisms underlying resistance to afatinib is key to addressing the clinical challenge of resistance to afatinib.

Numerous studies have revealed EGFR secondary mutations (12), oncogene amplification (13) and cellular phenotypic transitions (14), which contribute to the resistance to afatinib in NSCLC. However, 20–30% of drug-resistant cases remain to be elucidated, suggesting the involvement of unidentified mechanisms (15,16). Notably, tumors often undergo various adaptive responses during cancer treatment, which are key contributors to treatment failure and the development of resistance. Resistance caused by an adaptive response is a non-cell-autonomous mechanism by which tumor cells acquire transient drug resistance through the rapid reprogramming of signaling pathways, metabolism or phenotypic plasticity, without the genetic or epigenetic alterations (17,18). While this resistance may be reversed upon the cessation of the drug, it can also lay the groundwork for the subsequent development of stable resistance (19,20). Due to the dual role of short-term adaptive responses in the resistance to afatinib, longitudinal transcriptome profiling was conducted in the present study using A549 parental cells, afatinib-treated A549 cells, and afatinib-resistant A549-AR cells to identify resistance signatures and underlying signaling pathways throughout resistance development.

The aim of the present study was not to re-evaluate the canonical EGFR-targeted efficacy of afatinib, but to elucidate the cellular stress responses elicited by TKI exposure in NSCLC cells. Specifically, the present study investigated the mechanisms through which drug-induced adaptive responses, characterized by cytokine secretion, signaling reprogramming and stress-related protein upregulation, contribute to the initiation of resistance. To minimize confounding from constitutively active EGFR mutations, EGFR-wild-type NSCLC cell lines (A549 and H1299) were employed, thereby enabling the exploration of EGFR-independent adaptive mechanisms underlying TKI-induced resistance. Using RNA-sequencing (RNA-seq) analysis on both short-term afatinib-treated and stably resistant cells, the present study further delineated the dynamic transcriptional changes that accompany the transition from adaptive tolerance to stable resistance.

Materials and methods

Cell lines and culture

The A549 and H1299 cell lines, obtained from the National Collection of Authenticated Cell Cultures, were maintained in RPMI-1640 medium supplemented with 10% fetal bovine serum (GeminiBio LLC). The cells were incubated at 37°C in a humidified atmosphere with 5% CO2. The mutation profiles of EGFR, Ras and p53 in all NSCLC cell lines were confirmed according to the Catalogue Of Somatic Mutations In Cancer and American Type Culture Collection databases (Table SI). All cell lines were authenticated using short tandem repeat profiling and were routinely evaluated for mycoplasma contamination.

Drug preparation

Afatinib, osimertinib, selumetinib and SB431542 were purchased from Selleck Chemicals and dissolved in 100% DMSO at stock concentrations of 10 mM. All these drugs were stored at −20°C for long-term use.

Establishment of afatinib resistant cells

The A549 and H1299 cells were exposed to serial concentrations of afatinib (2, 5, 7.5, 10 and 12 µM) at 50% confluency for 3 days, with each concentration repeated for a minimum of two treatment cycles. Subsequently, the established afatinib-resistant (AR) cells (A549-AR and H1299-AR) were validated using the Cell Counting Kit-8 (CCK-8) assay. Briefly, parental and AR cells were treated with 1, 2, 5 and 10 µM afatinib for 48 h to validate afatinib resistance, and with 1, 2, 5 and 10 µM osimertinib for 48 h to assess cross-resistance.

CCK-8 assay

The cells were digested and resuspended to a concentration of 1×105 cells/ml. A total of 30 µl of the resuspended cell suspension was added into 1 ml culture medium and subsequently seeded into 96-well plates. Following overnight incubation at 37°C, the medium was replaced with fresh medium containing SB431542 or selumetinib at 1 or 10 µM, and the cells were incubated at 37°C for an additional 48 h. Subsequently, 10 µl CCK-8 reagent (Selleck Chemicals) was added to each well and the plates were incubated for 2–4 h at 37°C. The absorbance was measured at 450 nm using a microplate reader (Agilent Technologies, Inc.).

Reverse transcription-quantitative polymerase chain reaction (RT-qPCR)

Total RNA was extracted from the A549 cells, A549 cells treated with 20 µM afatinib for 48 h and A549-AR cells using TRIzol® reagent (Thermo Fisher Scientific, Inc.) according to the manufacturer's protocol. First-strand complementary DNA (cDNA) was synthesized from 1 µg total RNA using the Hifair® III 1st Strand cDNA Synthesis SuperMix according to the manufacturer's protocol (Shanghai Yeasen Biotechnology Co., Ltd.). qPCR was performed using ChamQ Universal SYBR® qPCR Master Mix (Vazyme Biotech, Co., Ltd.) on a RT-PCR system (Roche Diagnostics). Each 20-µl reaction mixture contained 10 µl ChamQ SYBR® qPCR Master Mix, 0.4 µl of forward primer (10 µM), 0.4 µl of reverse primer (10 µM), 1 µl of cDNA template and 8.2 µl of nuclease-free water. The thermal cycling conditions were as follows: Initial denaturation at 95°C for 30 sec, followed by 40 cycles of denaturation at 95°C for 10 sec and annealing/extension at 60°C for 30 sec. Relative gene expression was calculated using the 2−ΔΔCq method (21), with GAPDH serving as the internal reference gene. The primer sequences were designed by author WX and are listed in Table SII.

Western blotting

The A549 and A549-AR cells were lysed in RIPA buffer (X-Blot Technology Co., Ltd.) containing protease and phosphatase inhibitors on ice for 30 min. Lysates were centrifuged at 12,000 × g for 15 min at 4°C and the supernatants were collected. Protein concentrations were quantified using a BCA protein assay kit (Thermo Fisher Scientific, Inc.) following the manufacturer's instructions. Equal amounts of protein (20 µg per sample) were resolved on 10% gels using SDS-PAGE and transferred to PVDF membranes (MilliporeSigma). The membranes were blocked with 5% bovine serum albumin in TBST (20 mM Tris-HCl, 150 mM NaCl, 0.1% Tween-20 and pH 7.5) for 1 h at room temperature and incubated overnight at 4°C with primary antibodies against ERK1/2 (cat. no. 4695; 1:1,000; Cell Signaling Technology, Inc.), anti-phospho-ERK1/2 (Thr202/Tyr204; cat. no. 4370; 1:1,000; Cell Signaling Technology, Inc.) and GAPDH as the reference (cat. no.60004-1-Ig; Proteintech Group, Inc.). After washing with 1X TBST containing 0.1% Tween-20 (v/v), the membranes were incubated with HRP-conjugated anti-rabbit IgG secondary antibody [cat. no. BE0106; 1:5,000; Baioyijie (Beijing) Technology Co., Ltd.] or HRP-conjugated anti-mouse IgG secondary antibody [1:5,000; cat. no. BE0104; Baioyijie (Beijing) Technology Co., Ltd.] for 1 h at room temperature. Bands were visualized using an enhanced chemiluminescence kit (X-Blot Technology Co., Ltd.) and captured with the ChemiDoc XRS+ imaging system (Bio-Rad Laboratories, Inc.). Band intensities were quantified using ImageJ software (version 1.54; National Institutes of Health). The p-ERK/ERK ratio was calculated for quantitative analysis, and GAPDH was used as the loading control.

RNA-seq and differential gene analysis

RNA was extracted from A549 cells, A549 cells treated with 20 µM afatinib for 48 h and A549-AR cells using TRIzol reagent (Thermo Fisher Scientific, Inc.). The extracted RNA samples were submitted to Genewiz, Inc., for RNA sequencing as a commercial sequencing service, and subsequent RNA quality assessment, library preparation and sequencing were performed by the company. RNA quality and integrity were assessed using NanoDrop spectrophotometry (Thermo Fisher Scientific, Inc.) and RNA Quality Score analysis. Libraries were prepared using the VAHTS® Universal V8 RNA-seq Library Prep Kit for Illumina (cat. no. NR605; Vazyme Biotech Co., Ltd.). Paired-end sequencing with a read length of 150 bp (PE150) was performed (all other methods details are proprietary to Genewiz, Inc.). The RNA sequencing data were processed and analyzed in RStudio (version 2025.9.2.418; Posit Software), with differential gene expression calculated using the ‘DESeq2’ package (version 1.50.2; Bioconductor). The differential gene expression analysis was performed among the following comparisons: i) Treated A549 cells vs. A549 cells; ii) A549- AR cells vs. A549 cells; and iii) A549-AR cells vs. treated A549 cells. Gene Ontology (GO) enrichment analysis was performed for biological process, molecular function and cellular component categories (http://geneontology.org/). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was also performed (https://www.genome.jp/kegg/).

Time-series analysis

Time-series gene expression analysis was performed using the ‘maSigPro’ package (version 1.82.0; Bioconductor). Expression data were organized by time points and treatment conditions. Model fitting was performed using ‘p.vector’ with Q set to 0.05 and counts set to false, followed by backward stepwise regression using ‘T.fit’. Significant genes were identified using ‘get.siggenes’, applying an R-squared threshold of 0.7 and selecting all condition variables. Genes with adjusted P-values <0.05 were retained for visualization and downstream functional analysis.

Survival analysis for genes

Overall survival (OS) and disease-free survival (DFS) analyses were performed using the Gene Expression Profiling Interactive Analysis 2 (GEPIA2) web tool (http://gepia2.cancer-pku.cn/#survival), which is based on data from The Cancer Genome Atlas (TCGA; http://portal.gdc.cancer.gov/) and the Genotype-Tissue Expression databases (GTEx; http://gtexportal.org/home/). Survival curves were generated based on gene expression levels, with the group cut-off set to the quartile method to define the low-expression group as patients in the lowest 25% of gene expression and the high-expression group as patients in the highest 25% of gene expression. Patients with gene expression levels in the middle 50% were not included. For Kaplan-Meier curves showing crossover, P-values were estimated using the two-stage procedure as previously described (22,23) and implemented in the TSHRC package (version 0.1–6; http://cran.r-project.org/package=TSHRC). As individual-level raw survival data and full analysis settings were not available from GEPIA, the two-stage analysis was performed on survival coordinates digitized from the Kaplan-Meier plots using WebPlotDigitizer (version 5.2; Automeris; http://automeris.io/WebPlotDigitizer/); therefore, the resulting two-stage P-values should be interpreted as approximate. Hazard ratios (HRs) shown in the GEPIA2 survival plots were used as summary effect estimates over the full follow-up period.

Correlation analysis

ERK signaling-related genes, including EGFR, KRAS, NRAS, HRAS, BRAF, Raf1, MAP2K1, MAP2K2, MAPK3, MAPK1, FOS, FOSL1, FOSL2, Jun, JunB, JunD, dual-specificity phosphatase (DUSP)6, EGR1, CCND1, SPRY2, ETV4 and ETV5, were selected and entered into the first gene signature input in GEPIA2. The three-gene resistance signature (ABLIM3, HTR1D and HSPA1A) was entered into the second gene signature input. Pearson's correlation analysis was subsequently performed using TCGA-lung adenocarcinoma tumor dataset.

Meta-analysis

The HRs were obtained from the GEPIA2 survival plots. The corresponding 95% confidence intervals (CIs) were calculated based on the reported HRs and P-values. Subsequently, the HRs and their 95% CIs for differentially expressed genes or cancer types were pooled to calculate the overall effect estimates, applying a random-effects model. An HR >1 with the lower bound of the 95% CIs >1 indicated that the high expression level of the gene was significantly associated with reduced survival. A meta-analysis was performed using STATA software (version 16.0; StataCorp LP; http://www.stata.com/).

Statistical analysis

Data are presented as the mean ± standard error of the mean from at least three independent replicates. All statistical analyses were performed using GraphPad Prism software (version 10.6.1; Dotmatics). Statistical comparisons between two groups were performed using unpaired two-tailed Student's t-tests. When multiple treatment groups were compared, homogeneity of variance was first evaluated using the Brown-Forsythe (or Levene's) test; if the equal-variance assumption was met, one-way ANOVA followed by Tukey's or Dunnett's multiple-comparisons test was used and the corresponding adjusted P-values were reported. Specifically, Dunnett's multiple-comparisons test was applied for many-to-one comparisons in which all groups were compared with a single pre-specified reference group, whereas Tukey's multiple-comparisons test was used for analyses involving multiple pairwise comparisons among groups. If variances were unequal, Welch's one-way ANOVA was applied; post hoc multiple comparisons were performed using the Games-Howell test. When data did not meet parametric assumptions, the Kruskal-Wallis's test followed by Dunn's multiple-comparisons test was used. P<0.05 was considered to indicate a statistically significant difference.

Results

Cells with acquired resistance exhibit cross resistance to TKIs

To investigate the molecular mechanisms underlying acquired resistance to afatinib, AR cell lines were established. Initially, the optimal concentration of afatinib was determined based on cytotoxic effects. A range of afatinib concentrations was administered to the H1299 and A549 cells, revealing that a dose of 2 µM reduced cell viability by 20% (Fig. 1A and B). Therefore, the induction of acquired resistance was initiated using 2 µM afatinib, as this concentration provided sufficient selective pressure while still allowing surviving cells to recover and gradually adapt during subsequent treatment cycles. The H1299-AR cells exhibited a significantly reduced sensitivity to afatinib, as evidenced by the significantly lower inhibition of cell viability compared with the parental H1299 cells under identical treatment conditions (P<0.01; Fig. 1C). A similar resistance phenotype was observed in the A549-AR cells (P<0.01; Fig. 1D). Furthermore, both the H1299-AR and A549-AR cells demonstrated cross-resistance to osimertinib treatment, as shown by significantly higher cell viability at 5 and 10 µM in H1299-AR cells and at 1, 2, 5 and 10 µM in A549-AR cells compared with their parental cells (Fig. 1E and F). These findings suggested that the acquired resistance in NSCLC cells is not limited to afatinib, but extends to third-generation EGFR-TKIs, such as osimertinib. This suggested that the underlying mechanism may be independent of the EGFR mutation status and instead may be associated with adaptive responses contributing to acquired resistance.

Figure 1. Cells with acquired resistance to afatinib exhibit cross resistance to tyrosine kinase inhibitors. (A and B) Relative cell viability of (A) H1299 and (B) A549 cells treated with serial concentrations of afatinib (n=3). (C) Relative cell viability of H1299 and H1299-AR cells treated with serial concentrations of afatinib (n=3). (D) Relative cell viability of A549 and A549-AR cells treated with serial concentrations of afatinib (n=3). (E) Relative cell viability of H1299 and H1299-AR cells treated with serial concentrations of osimertinib (n=3). (F) Relative cell viability of A549 and A549-AR cells treated with serial concentrations of osimertinib (n=3). **P<0.01, compared with parental cells. AR, afatinib-resistant.

Transcriptome dynamics and functional enrichment reveal molecular transition from drug sensitivity to resistance

To investigate the adaptive responses driving the transition from drug sensitivity to stable acquired resistance, the transcriptome profiling of sensitive, treated and resistant cells was conducted (Fig. 2A). Hierarchical clustering analysis revealed distinct segregation among the three groups, with the treated cells exhibiting an intermediate transcriptional state, indicative of a stepwise progression toward resistance (Fig. 2B). Gene expression trajectories revealed a variety of dynamic patterns, such as high-low-high, high-high-low, low-high-high, low-high-low and low-low-high, emphasizing the heterogeneous nature of transcriptional reprogramming during resistance acquisition (Fig. 2C). The low-high-high pattern observed in cluster 3 was characterized by constitutively upregulated genes that potentially contribute to afatinib acquired resistance. To further investigate the functional significance of these persistently upregulated genes, GO analysis was conducted on cluster 3.

Figure 2. Transcriptome dynamics and functional enrichment reveal molecular transition from drug sensitivity to resistance. (A) Flowchart of RNA-sequencing analysis in A549 sensitive cells, A549 cells treated with afatinib for 48 h and A549-AR cells. (B) Heatmap of the significant genes after Time-series analysis of gene expression using ‘maSigPro’ package (n=3). (C) Time-course analysis of differential gene expression across sensitive, treated and resistant states. (D-F) GO (D) molecular function, (E) biological process and (F) cell component analysis for genes in cluster 3. GO, Gene Ontology; AR, afatinib resistant.

GO molecular function analysis revealed significant enrichment in ‘cadherin binding’, ‘phosphatase regulator activity’ and ‘protein phosphatase regulator activity’ (P<0.05; Fig. 2D), suggesting the enhanced regulation of kinase-phosphatase signaling and cell-cell adhesion, which may contribute to the resistant phenotype. GO biological process analysis identified pronounced enrichment in RNA-related functions, including ‘regulation of mRNA metabolic process’, ‘regulation of mRNA processing’ and ‘RNA splicing’ (Fig. 2E), implicating elevated transcriptional and post-transcriptional activity in the maintenance of resistance. Furthermore, cellular component analysis indicated significant enrichment in ‘spindle’, ‘focal adhesion’ and ‘cell-substrate junction’ (P<0.05; Fig. 2F), consistent with increased mitotic activity and dynamic cytoskeletal remodeling. These functional analyses suggested that resistance to afatinib is associated with enhanced kinase-phosphatase signaling, cell adhesion, and increased transcriptional and post-transcriptional activity.

Identification and characterization of differentially expressed genes during the transition from drug sensitivity to resistance

To identify the constitutively upregulated gene signature contributing to acquired resistance to afatinib, genes with a ≥2-fold increase in expression were selected in both treated and AR cells. As presented in the volcano plots, the treated cells displayed a moderate number of differentially expressed genes, while the resistant cells exhibited a more extensive transcriptional shift, reflecting a progressive reprogramming of gene expression (Fig. 3A). Venn diagram analysis further revealed that a subset of upregulated (n=63) and downregulated (n=219) genes was shared between the treated and resistant groups, while the majority were uniquely enriched in resistant cells (Fig. 3B), suggesting a gradual development toward a resistant state. Hierarchical clustering demonstrated clear segregation among the three groups, with treated cells occupying an intermediate transcriptional position, consistent with a gradual acquisition of resistance (Fig. 3C). Furthermore, representative gene trajectories revealed distinct dynamic expression patterns, including persistent upregulation of resistance-associated genes such as JUN, HSPA1A and SYTL2 (Fig. 3D), indicating a time-dependent shift in gene expression. Therefore, these findings suggested that resistance to afatinib develops through a gradual, multi-phase alteration in gene expression.

Figure 3. Identification and characterization of differentially expressed genes during the transition from drug sensitivity to resistance. (A) Volcano plot of the differential genes with 2-fold changes in treated and resistant cells compared with sensitive cells, n=3. (B) Overlapping of the upregulated and downregulated genes with 2-fold changes. (C) Heatmap images for the 63 continuously upregulated genes or 219 downregulated genes in the treated and resistant groups (n=3). (D) Comparison of the top 30 upregulated genes in treated or resistant cells versus sensitive cells (n=3). Statistical significance was assessed by one-way ANOVA followed by Dunnett's multiple-comparisons test versus the control; *P<0.05, ***P<0.001. ns, not significant; FPKM, fragments per kilobase of transcript per million mapped reads.

Validation of afatinib resistance-associated gene signature and its prognostic relevance across patient cohorts

To validate the gene signatures associated with resistance to afatinib and assess their clinical relevance, RT-qPCR was performed on selected candidates. Several genes, including BCAT1, GNG2, HSPA1A, ABLIM3, IL31RA, ROR1, SLCO2B1, HTR1D, and TBXAS1, were consistently upregulated in the resistant cells compared with those in both the parental and afatinib-treated cells, consistent with the RNA-seq findings (Fig. 4A). HR analysis revealed that multiple genes, such as ABLIM3, HTR1D and HSPA1A, were significantly associated with both poor OS and DFS in clinical lung cancer datasets (Figs. 4B and C, S1 and S2). Furthermore, stage-specific expression analysis revealed the progressive upregulation of ABLIM3, HTR1D and HSPA1A from stage I to IV tumors, suggesting their involvement in tumor progression (Fig. 4D-F). The three-gene signature stratified patients into distinct risk groups, with individuals in the high-signature group demonstrating significantly shorter OS (HR=1.7; two-stage P<0.01) and DFS (HR=2.1; P<0.01) times (Fig. 4G and H). Notably, this signature retained prognostic significance for both OS and DFS across multiple cancer types, including liver, lung and pancreatic cancer (Figs. 4I and S3). Collectively, these findings identified ABLIM3, HTR1D and HSPA1A as an afatinib resistance-associated gene signature with prognostic value, underscoring their potential roles in tumor progression and patient risk stratification.

Figure 4. Validation of the afatinib resistance-associated gene signature and its prognostic relevance across patient cohorts. (A) Reverse transcription-quantitative PCR analysis of candidate gene expression in drug-sensitive (A549), treated A549, and resistant (A549-AR) cells (n=3). Statistical significance was assessed by one-way ANOVA followed by Dunnett's multiple-comparisons test versus the drug-sensitive group. (B and C) Forest plots summarizing HRs with 95% CIs for (B) OS and (C) DFS in patients with lung cancer with high expression versus low expression of the indicated genes. (D-F) Relative expression levels of (D) ABLIM3, (E) HTR1D and (F) HSPA1A in different stages of lung cancer. (G) OS and (H) DFS of patients with lung cancer with high expression of the 3-gene signature compared with those with low expression. For OS, Kaplan-Meier curves exhibiting crossover were compared over the full follow-up period using the two-stage procedure. For DFS, Kaplan-Meier curves were compared over the full follow-up period using the log-rank P-value. (I) OS and DFS of patients with various cancer types, comparing high versus low expression levels of ABLIM3, HTR1D and HSPA1A. *P<0.05; ***P<0.001. HTR1D, 5-hydroxytryptamine receptor 1D; ABLIM3, actin-binding LIM protein family member 3; HSPA1A, heat shock protein family A member 1A; OS, overall survival; DFS, disease-free survival; HR, hazard ratio; CI, confidence interval; DL, DerSimonian-Laird; Pr, probability; ns, not significant.

Pathway enrichment reveals MAPK/ERK signaling as a major contributor to resistance to afatinib

To elucidate the activated signaling pathways in adaptive response that also contribute to the stable resistance to afatinib, KEGG pathway enrichment analysis was conducted on upregulated genes from both the treated and resistant cells. In the treated cells, the ‘MAPK signaling pathway’ was identified as the most significantly enriched (P<0.01), whereas the resistant cells demonstrated an additional enrichment of the ‘PI3K-AKT signaling pathway’ (P<0.01), the ‘Ras signaling pathway’ (P<0.01) and ‘TGF-β’ (P<0.01) signaling (Fig. 5A and B), indicating the progressive and coordinated activation of pathways during resistance acquisition. To validate the activation of MAPK/ERK signaling in drug-treated and AR cells, the phosphorylation levels of ERK1/2 were examined using western blotting analysis. The results revealed that the p-ERK/ERK ratio was elevated by 1.7-fold in the treated A549 cells compared with that in the parental cells (P<0.05; Fig. 5C). Similarly, a 1.3-fold increase in the p-ERK/ERK ratio was observed in the AR cells (P<0.01; Fig. 5C), indicating the sustained activation of the MAPK/ERK pathway in the development of resistance to TKI. Furthermore, a correlation analysis demonstrated a significant positive correlation between the MAPK/ERK pathway activity and the expression levels of the three-gene signature (P<0.01; r=0.42; Fig. 5D). Collectively, these findings indicated that the MAPK/ERK signaling pathway was activated in the development of resistance to afatinib and was positively correlated with the identified three-gene signature. Due to previous reports implicating the role of TGF-β in TKI resistance across multiple malignancies, including lung and gastric cancer (22–24), TGF-β signaling was also selected for further investigation using specific inhibitors. To functionally assess the contribution of these pathways, A549-AR cells were treated with afatinib alone or in combination with either the TGF-β receptor inhibitor, SB431542 or the MEK inhibitor, selumetinib. Co-treatment with SB431542 led to no additional effect compared with afatinib alone, while co-treatment with 10 µM selumetinib significantly reduced the viability of resistant cells compared with afatinib treatment alone (P<0.01; Fig. 5E). These results demonstrated that MAPK/ERK signaling serves a functional role in mediating resistance to afatinib. Notably, the inhibition of MAPK/ERK signaling in H1299-AR cells restored sensitivity to afatinib, resulting in a >20% decrease in cell viability compared with afatinib monotherapy (P<0.01; Fig. 5F), further supporting its role as a major contributor to the resistant phenotype.

Figure 5. Pathway enrichment reveals MAPK/ERK signaling as a major contributor to resistance to afatinib. (A and B) KEGG enrichment of the upregulated genes in the (A) treated and (B) resistant groups. (C) Representative western blot images of ERK1/2 and p-ERK1/2, and quantification of the p-ERK/ERK ratio in A549 cells treated with 20 µM afatinib or in A549-AR cells. *P<0.05, **P<0.01, compared with the untreated group or A549 cells. (D) Correlation of the three-gene signature with ERK1/2 signaling in TCGA lung cancer data. (E) Relative cell viability of A549-AR cells treated with 10 µM afatinib in combination with either 1 or 10 µM SB431542 or selumetinib (n=3). Statistical significance was assessed by one-way ANOVA followed by Dunnett's multiple-comparisons test versus afatinib (10 µM). Adjusted P-values (Dunnett) are reported. (F) Relative cell viability of H1299-AR cells treated with 10 µM afatinib in combination with either 1 or 10 µM SB431542 or selumetinib (n=3). Statistical significance was assessed by one-way ANOVA followed by Dunnett's multiple-comparisons test versus the afatinib (10 µM) group. Adjusted P-values (Dunnett) are reported. ***P<0.001 and ns, no significance, compared with afatinib-treated group. AR, afatinib-resistant; p, phosphorylated; KEGG, Kyoto Encyclopedia of Genes and Genomes; TPM, transcript per million; TCGA, The Cancer Genome Atlas.

Inhibition of ERK signaling also reverses resistance to osimertinib

The aforementioned results showed that AR cells also exhibit resistance to osimertinib, suggesting that drug-induced adaptive responses may drive a broader, cross-resistant phenotype. Since the pharmacological inhibition of ERK signaling was sufficient to restore afatinib sensitivity, it was hypothesized that targeting this pathway may similarly reverse resistance to osimertinib. To validate this hypothesis, AR cells were treated with osimertinib in combination with either SB431542 or selumetinib. SB431542 alone had no inhibitory effect compared with the untreated control group and its combination with osimertinib did not enhance inhibition compared with osimertinib alone (Fig. 6A). By contrast, selumetinib monotherapy modestly reduced cell viability, while its combination with osimertinib produced a marked synergistic effect, resulting in an additional 20% reduction in viability compared with osimertinib treatment alone (P<0.01; Fig. 6A). Furthermore, selumetinib independently attenuated osimertinib resistance, as evidenced by a 10% greater suppression of cell viability compared with osimertinib monotherapy (P<0.01; Fig. 6B). These findings further support the role of sustained MAPK/ERK signaling in mediating cross-resistance and highlight the potential of targeted inhibition of ERK to restore sensitivity to EGFR-TKIs.

Figure 6. Inhibition of ERK signaling also reverses resistance to osimertinib. (A and B) Relative cell viability of (A) A549-AR and (B) H1299-AR cells treated with 10 µM osimertinib in combination with either 1 or 10 µM SB431542 or selumetinib (n=3). Statistical significance was assessed by one-way ANOVA followed by Tukey's multiple-comparisons test for comparisons between the osimertinib-alone group and the corresponding combination groups, and between 10 µM selumetinib alone and 10 µM selumetinib plus osimertinib. Adjusted P-values are reported. *P<0.05, **P<0.01 and ***P<0.001; ns, no significance. AR, afatinib-resistant.

Discussion

The present study demonstrated that acquired resistance to afatinib also conferred cross-resistance to the third-generation EGFR inhibitor osimertinib, indicating that adaptive cellular responses, rather than secondary EGFR mutations, may serve a dominant role in the development of resistance. Mechanistically, resistance was associated with transcriptional reprogramming and the activation of the MAPK/ERK signaling pathway. A three-gene signature (ABLIM3, HTR1D and HSPA1A) was identified as a potential biomarker of resistance to afatinib and was positively correlated with MAPK/ERK pathway activity. Notably, the pharmacological inhibition of MAPK/ERK signaling using selumetinib effectively re-sensitized resistant cells to both afatinib and osimertinib. These findings highlighted MAPK/ERK pathway activation as a major contributor to resistance to EGFR-TKI and suggested that its targeted inhibition may represent a promising therapeutic strategy for NSCLC.

Since afatinib targets both wild-type and mutant EGFR, the present study employed EGFR-wild-type NSCLC models to specifically investigate EGFR-independent adaptive responses to TKI exposure. This approach minimizes the confounding effects of constitutively active EGFR mutations and enables a clearer characterization of stress-induced transcriptional programs that may precede or promote genetic resistance. By focusing on this model, the present study aimed to capture the early, non-mutational stage of resistance development, an adaptive process increasingly recognized as a major contributor to TKI treatment failure and a potential therapeutic target distinct from mutation-driven mechanisms.

Transcriptome analysis revealed that resistance of NSCLC cells to afatinib develops through a gradual transition from a drug-sensitive state to a stable resistant phenotype, beginning with an acute drug-induced stress response followed by persistent transcriptional reprogramming, consistent with previously reported heterogeneity in resistance mechanisms (16,25). The functional enrichment of resistance-associated gene clusters highlighted the aberrant regulation of kinase-phosphatase signaling, the upregulation of RNA metabolic activity, and the reorganization of cytoskeletal and adhesion structures, molecular features of cellular plasticity and therapeutic evasion. Among the oncogenic pathways, MAPK/ERK signaling emerged as a principal driver, with KEGG analysis indicating its early activation following exposure to afatinib. This activation was accompanied by the coordinated involvement of the PI3K/AKT, Ras and TGF-β pathways, which collectively contribute to the maintenance of the resistant phenotype (26–28). Notably, the modulation of the PI3K/AKT cascade has also been implicated in reversing the resistance of NSCLC cells to osimertinib, underscoring its interplay with MAPK/ERK signaling in acquired resistance to TKIs (29). Notably, pharmacological inhibition of ERK1/2 with selumetinib effectively reversed the resistance to both afatinib and osimertinib, supporting the key involvement of MAPK/ERK signaling in both the development and maintenance of acquired resistance to TKIs (30). These findings emphasize the complexity of dynamic, multi-pathway cooperation during resistance acquisition, suggesting that future therapeutic strategies may require the combined targeting of multiple signaling axes. Similar combined approaches, such as use of β-elemene in combination with hyperthermia, have been reported to enhance antitumor efficacy in NSCLC models (31), supporting the rationale for multi-target intervention.

The present study systematically investigated the role of the MAPK/ERK pathway in the development of resistance to EGFR-TKIs. It was demonstrated that the persistent activation of this pathway not only drives acquired resistance to the second-generation EGFR inhibitor, afatinib, but also mediates cross-resistance to the third-generation inhibitor, osimertinib, consistent with the findings of previous studies (32,33). Persistent ERK phosphorylation appears to promote stable signaling reprogramming independent of EGFR mutations, reflecting adaptive plasticity distinct from classical resistance mechanisms such as T790M or C797S mutations (34). This EGFR-independent resistance mode is consistent with previously reported bypass mechanisms involving the aberrant activation of signaling pathways, such as Src family kinase/focal adhesion kinase (35) or KRAS amplification (36).

Pharmacological validation confirmed the key role of MAPK/ERK signaling, as ERK inhibition restored the sensitivity to afatinib and osimertinib. These findings underscore MAPK/ERK activation as a central mediator of adaptive resistance. Furthermore, the present study identified a three-gene resistance signature (ABLIM3, HTR1D and HSPA1A), which was consistently upregulated in resistant cells and positively correlated with MAPK/ERK activity. Notably, HSPA1A can enhance DUSP1 expression, a negative regulator of MAPK signaling (37), complementing previous research demonstrating that the loss of DUSP1 or DUSP6 leads to ERK activation and drug resistance (34). Therefore, these results delineate a mechanistic association between stress-induced transcriptional reprogramming and sustained ERK signaling in TKI-acquired resistance.

At a broader conceptual level, resistance to anticancer drugs is increasingly understood as an adaptive process within a complex pathological ecosystem rather than a phenomenon driven solely by genetic mutations (38). Tumor evolution involves dynamic interactions among cancer cells, stromal and immune components, and metabolic constraints that collectively shape therapeutic outcomes (39). From this ecological and evolutionary perspective, the transcriptional plasticity and signaling rewiring observed in AR cells likely reflect adaptive reprogramming under selective TKI pressure, maintaining the equilibrium between sensitive and resistant subpopulations (40,41). Integrating ecological and evolutionary principles with molecular insights may provide a more comprehensive understanding of acquired resistance and lead to the development of adaptive therapeutic strategies.

Notably, the ABLIM3-HTR1D-HSPA1A signature has demonstrated robust prognostic value across multiple solid tumors, including colorectal carcinoma (42), hepatocellular carcinoma (43) and pancreatic cancer, highlighting its potential as a pan-cancer biomarker for resistance-associated risk stratification. Compared with previously proposed single-gene markers such as neurotrophic tyrosine receptor kinase 2 (44) or AXL (45), this multi-gene signature provides greater predictive accuracy and broader clinical relevance. Building on these findings, the three-gene signature also holds translational potential as a clinically accessible biomarker. Among the three genes, HSPA1A, a member of the HSP70 family, is detectable in the circulation and is associated with tumor burden and treatment response (46,47). Although ABLIM3 has been detected at low abundance in plasma proteomic datasets (48) and HTR1D currently lacks circulating evidence, their co-expression patterns suggest possible applicability for liquid biopsy-based monitoring. Future studies are warranted to determine whether longitudinal changes in circulating or exosomal HSPA1A, ABLIM3 and HTR1D levels reflect emerging resistance or treatment response in patients with NSCLC receiving EGFR-TKIs. Validation in xenograft and patient-derived xenograft models with serial plasma sampling will be essential to confirm the in vivo detectability and prognostic relevance of this gene signature. Such efforts may enable the real-time, non-invasive tracking of adaptive resistance and accelerate the translation of these mechanistic insights into clinical practice.

Despite the comprehensive elucidation of the molecular mechanisms underlying the development of acquired resistance to afatinib in the present study, several limitations should be acknowledged. First, the findings presented herein are primarily based on in vitro models of NSCLC and lack in vivo validation using xenografts or patient-derived xenografts models. This limitation may affect the translational relevance of the results. Second, while the MAPK/ERK signaling pathway was identified as a major contributor to the resistant phenotype, the upstream regulatory mechanisms governing its activation remain incompletely understood. In particular, the involvement of alternative receptor tyrosine kinases, adaptor proteins or cellular stress-related pathways in sustaining ERK signaling requires further investigation. Furthermore, synergy metrics (for example, Chou-Talalay) were not included in the present study, as the primary objective was to delineate single-agent-induced adaptive responses. Future studies are thus warranted to incorporate systematic combination index analyses to quantitatively define synergistic or additive effects between EGFR-TKIs and MEK inhibitors. Lastly, the prognostic evaluation of ABLIM3, HTR1D and HSPA1A was based on transcriptomic datasets rather than clinical specimens. Thus, the future validation of their protein expression in NSCLC tissues using immunohistochemistry will be essential to confirm their prognostic significance and clinical applicability. In addition, Kaplan-Meier curves with crossover, particularly the DFS curve of HSPA1A and the OS curve of the three-gene signature in lung cancer, should be interpreted with caution. Although a two-stage analysis was additionally performed for these curves, it was based on reconstituted data derived from digitized Kaplan-Meier plots rather than original individual-level survival data.

Supplementary Material

Supporting Data

Supporting Data

Acknowledgements

The authors would like to thank Dr Ying Sun, a postdoctoral researcher at the School of Medicine, Shanghai Jiao Tong University (Shanghai, China), for providing help in operating RStudio.

Funding

The present study was funded by the National Natural Science Foundation of China (grant no. 82204673), the Budget Research Project of Shanghai University of Traditional Chinese Medicine (grant no. 2021LK018) and the Priority Incubation Program (grant no. A1-U24-205-030306).

Availability of data and materials

The data generated in the present study may be found in the Genome Sequence Archive in the National Genomics Data Center, China National Center for Bioinformation/Beijing Institute of Genomics, Chinese Academy of Sciences under accession number GSA-Human: HRA012241 or at the following URL: https://ngdc.cncb.ac.cn/gsa-human/browse/HRA012241. The remaining data generated in the present study may be requested from the corresponding author.

Authors' contributions

CQ participated in data acquisition, analysis and manuscript revision. WZ and DT conceived the study concept, and performed data interpretation and revision of the manuscript. YY and BL contributed to the acquisition of data and revised the manuscript. ZH and YZ contributed to the interpretation of data and revised the manuscript. TY conceived the study concept, drafted the manuscript and revised the manuscript. WX conceived the study concept, contributed to data acquisition, analysis and interpretation, drafted the manuscript and revised the manuscript. All authors have read and approved the final manuscript. CQ and WX confirm the authenticity of all the raw data.

Ethics approval and consent to participate

Not applicable.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References