The GEO database is a high-throughput microarray and sequence functional genomic database (https://www.ncbi.nlm.nih). In this study, the GSE62504 dataset consisted of two afatinib-resistant HCC827 replicates and one parental HCC827 cell line (14). The GSE75037 dataset (15) included 83 lung adenocarcinoma (LUAD) samples and 83 matched adjacent lung samples, whereas the GSE18842 dataset contained 32 squamous cell carcinoma, 14 adenocarcinoma, and 45 adjacent lung tissues (16).

The original microarray data files of these three downloaded datasets were analyzed through GEO2R (https://www.ncbi.nlm.nih.gov/geo/geo2r/). This online tool can be used to compare two groups under the same experimental settings or more sets of samples (17). The P<0.05 and |fold change (FC)| > 1.0 were set as the cutoff standards to define the DEGs.

To perform functional annotation of DEGs, we applied annotation, visualization, and used a comprehensive database (DAVID, https://david.abccncifcrf.gov/) (18) to perform the Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses, specifically GO enrichment analysis and KEGG pathway analysis.

The gene expression array dataset of ONCOMINE (www.oncomine.org) is a publicly accessible, online cancer microarray database that helps facilitate research using genome-wide expression analyses (19). For DEGs, comparison between cancer specimens and normal control dataset analysis was performed.

The TIMER database was used to systematically analyze tumor-infiltrating immune cells (TIICs) in 32 cancer types using more than 10,000 samples from The Cancer Genome Atlas (TCGA) (https://cistrome.shinyapps.io/timer/) database (20). We initially used this database to assess differences in DEG expression levels in tumor types and explored the relationship between the expression of BIRC5 and the abundance of immune infiltrates in lung squamous cell carcinoma (LUSC) and LUAD by considering P<0.05 as the cut-off criterion for statistical significance.

The platform UALCAN is an omnibus and interactive web-based tool based on The Cancer Genome Atlas (TCGA) for deep analysis of gene expression using genomics data from 31 cancer types (21). The database UALCAN was used to analyze the correlation between BIRC5 mRNA transcriptional levels in LUSC and LUAD patients with different stage, gender, age, smoking habits, and lymph node metastasis. P<0.05 was regarded as indicating statistically significant results.

Kaplan-Meier plotter (https://www.kmplot.com) is an online survival analysis tool consisting of 10,461 cancer samples (including samples from 5143 breast cancer, 1816 ovarian cancer, 2437 lung cancer, and 1065 gastric cancer patients) that can facilitate evaluation of the impact of 54,675 gene pairs on overall survival OS (22). Based on the expression levels of DEGs, lung cancer patients were divided into a high expression group and a low expression group, the first progression(FP), OS and post progression survival (PPS) rate was further analyzed.

The LinkedOmics database contains multi-omics data and clinical data for 32 cancer types and 11,158 patients from the TCGA project (23). The “LinkInterpreter” module was used to derive biological insights and perform analysis of kinase targets, miRNA targets, and transcription factor targets for BIRC5 chemokines. Gene Set Enrichment Analysis (GSEA) was used to perform analyses with a minimum number of genes (size) of 3 and a simulation of 500 within the LUAD dataset. Results were analyzed statistically using the Pearson’s correlation test. The P value cutoff was 0.05.

GeneMANIA (https://www.genemania.org/) was used to predict the potential function of the BIRC5 gene and to predict the functions of specific genes (24). After obtained the kinase networks using LinkedOmics, we put those genes into GeneMANIA and understood their potential functions.

GSCALite, a bioinformatics platform for gene set cancer analysis, offered several type of analyses, including methylation analysis, cancer-related pathway analysis, miRNA network analysis, etc. (25). GSCALite provided pathway activity analysis and drug sensitivity analysis in our study with TCGA SKCM sample. The spearman correlation was used to explore the correlation between the gene expression and drug sensitivity.

Afatinib-sensitive NSCLC cell lines with mutated EGFR (DelE746-A750), HCC827 (Keygenbio, Nanjing, China) were cultured in RPMI-1640 medium (Gibco, USA) supplemented with 10% fetal bovine serum (Gemini, USA), penicillin (100 U/mL) and streptomycin (100 mg/mL). All cell lines were maintained in a humidified incubator at 37°C with 5% CO₂.

The HCC827 cells were cultured in a medium containing afatinib (Selleck Chemicals, Houston, TX, USA)at 1 nmol/L. When the cells showed viability similar to that of cells without afatinib, we gradually increased the concentration of afatinib until it reached 6 μmol/L. The entire induction period lasted approximately six months.

HCC827 and HCC827-AR cells (3×10³) were seeded in 96-well plates, grown overnight, and then treated with varying drug concentrations for 72 h; then, serum-free medium was replaced. Next, 10 μL of CCK8 (Beyotime, Shanghai, China) was added to each well for incubation at 37°C with 5% CO₂ for 2 h. The OD value was measured at 450 nm wavelength (Thermo, Waltham, MA, USA).

The cells were lysed with RIPA buffer (Beyotime, Beijing, China) in the presence of protease inhibitor PMSF (Beyotime, Beijing, China). Protein samples were quantified using BCA protein assay kit (Beyotime, Beijing, China). Equal amounts of proteins (25 μg) were separated on 15% SDS-polyacrylamide gels and transferred to polyvinylidene fluoride membranes (Beyotime, Beijing, China). Then, the membrane was blocked with defatted milk at room temperature for 2 h and incubated with anti-BIRC5 rabbit polyclonal antibody (1:500; Cell Signaling Technology, #2808S) and anti-GAPDH rabbit monoclonal antibody (1:6000; Cell Signaling Technology, #5174S) at 4°C overnight. Next, the membrane was incubated with fluorescence-labeled anti-rabbit IgG antibody (1:1000; Cell Signaling Technology, #7074) at room temperature for 2 h and then visualized using enhanced chemiluminescence (Thermo, Waltham, MA, USA).

GraphPad Prism 8 software (GraphPad Software, Inc., USA) statistical software was used to process the data, and P<0.05 was considered to define statistical significance.

To identify the potential genes conferring afatinib resistance in NSCLC, we investigated GSE62504 using |log2 (FC)|>1.0 and P <0.05 as thresholds to screen for differential genes ( Figure 1A ). A total of 1483 DEGs were screened, of which 700 DEGs were found to be significantly upregulated and 783 DEGs were downregulated ( Figure 1B ). To characterize the functions of these significantly upregulated DEGs, GO and KEGG analyses were then performed as previously described. The top five GO terms and enrichment pathways were determined, and these significantly upregulated DEGs were found to be highly associated with cell adhesion, inhibition of apoptosis, and inhibitor RNA transcription height ( Figure 1C ). As shown in Figure 1D , the DEGs that were simultaneously upregulated were also enriched in intercellular junctions and phagosomes, cell adhesion pathways, and PI3K-Akt signaling pathways.

Tumor gene heterogeneity is one of the determinants of drug resistance during tumor treatment (26, 27). To screen for genes related to the occurrence of afatinib resistance, the GSE75037 dataset consisting of lung cancer and adjacent normal tissues was further searched through the GEO database. A total of 984 DEGs were screened ( Figure 2A ), of which 312 genes were significantly upregulated and 672 genes were downregulated in cancer tissues ( Figure 2B ). Interference agents with mitotic arrest in cells promotes their death, which is a successful clinical strategy (28). The top GO terms and KEGG pathways of upregulated DEGs affected nuclear division and mitosis ( Figure 2C ) and the cell cycle ( Figure 2D ). Gene expression control mechanisms are complex and unique, and the integration of transcription and proteomic data provides additional information about gene expression control that cannot be obtained from data from a single source (29). To verify the reliability of the results, the GSE18842 datasets composed of data on lung cancer and normal tissues were retrieved. As the results show in Figure 3A , 990 DEGs were identified, among which 402 DEGs were up-regulated and 588 DEG were significantly down-regulated ( Figure 3B ). The GO terms were highly consistent with those for the GSE75037 dataset ( Figure 3C ), while the DEGs were enriched in cell cycle pathways ( Figure 3D ). Upon comparing the DEGs significantly upregulated in GSE62504, GSE18842, and GSE75037 ( Figure 4 ), six genes were identified, namely BIRC5, MCM4 (30), SPP1 (31), NMU (32), CTHRC1 (33), and UHRF1 (34) had reported on drug resistance. However, little research has explored BIRC5 in the context of afatinib resistance in lung cancer (35). Therefore, we chose BIRC5 as the potential target gene of interest in this study.

To determine the role of BIRC5 in lung cancer, we first evaluated its expression and diagnostic and prognostic value in patients with NSCLC. Oncomine data revealed that mRNA expression of BIRC5 was significantly higher in NSCLC tissues than in normal tissues ( Figures 5A, B ). According to the TCGA database, mRNA expression of BIRC5 in LUSC and LUAD tissues was higher than in normal tissues ( Figure 5C, D ). Using the human protein profile database, it was further verified that BIRC5 was positively expressed in NSCLC tissues and negatively expressed in normal lung tissues ( Figure 5E ). Then, we analyzed the transcription levels of BIRC5 by tumor stage, patients’ age, patients’ gender, smoking habits, and nodal metastasis status for LUSC and LUAD. Regardless of tumor stage, gender, age, smoking habits, and nodal metastasis status, BIRC5 transcription levels in tissues were significantly higher than in normal lung tissues ( Figure 6 ). Furthermore, we investigated the correlation between BIRC5 overexpression at the mRNA level and patient prognosis by plotting and comparing FP, OS and PPS of LUSC and LUAD patients with that of healthy individuals through Kaplan-Meier plotter ( Figure 7 ). BIRC5 overexpression was associated with worse FP (HR=3.13 (2.23-4.4), P<0.001), OS (HR=2.42 (1.9-3.09), P<0.001) in LUAD,which was negative in LUSC. Overall, the findings above imply that the mRNA expression of BIRC5 is remarkably correlated with LUAD patient survival rates, and BIRC5 expression may represent a promising biomarker for prediction of survival in LUAD patients. All of the above data indicate that the upregulation of BIRC5 expression levels promotes the development and progression of NSCLC.

Great progress has been made in immunotherapy for the treatment of advanced lung cancer (36). However, the role of BIRC5 in immune infiltrates in NSCLC is unknown. Using the TIMER database, we investigated the relationship between the transcriptional level of BIRC5 and immune infiltration. For LUSC, the levels of BIRC5 expression correlated positively with tumor purity (Cor=0.336, P=4.47e-14), while BIRC5 expression was negatively correlated with CD4+T cells (Cor=-0.262, P=6.53e-09, macrophages (Cor=-0.289, P=1.16e-10), neutrophils (Cor=-0.137, P=2.83e-03), and dendritic cells (Cor=-0.218, P=1.68e-06); there was no significant association between B cells and CD8+T cells. We also analyzed the distribution of immune infiltration in relation to LUAD mutated with BIRC5. The levels of BIRC5 expression correlated negatively with B cells (Cor=-0.225, P=5.95e-07), CD4+T cells (Cor=-0.185, P=4.05e-05), and dendritic cells (Cor=-0.093, P=4.08e-02), while there was no significant association between tumor purity, CD8+T cells, macrophages, neutrophils ( Figure 8 ). The relationship between BIRC5 expression level and immune markers was analyzed. As shown in Table 1 , the expression of BIRC5 correlated significantly with the expression of the marker genes of different subsets of T cells in LUAD, namely,B cell markers, CD19, CD79A;neutrophil markers CD66b (CEACAM8), CD11b (ITGAM), and CCR7;dendritic cell markers, HLA-DPB1, HLA-DQB1, HLA-DRA, HLA-DPA1, and BDCA-1 (CD1C); Th2 markers, STAT6 and STAT5A; Tfh markers, BCL6; Th17 markers, STAT3; Treg markers, STAT5B and TGFb (TGFB1); T cell exhaustion markers, LAG3 and GZMB. BIRC5 was negatively correlated with a variety of immune cells, suggesting that this gene is associated with immunotherapy resistance in NSCLC.

To gain insight into its biological significance, we examined BIRC5 co-expression and dark patterns, performed functional and enrichment analysis, and identified regulators of this gene in the LUAD cohort. The 5922 genes (red dots) showed positive correlations with BIRC5 and 7497 genes (green dots) showed negative correlations ( Figure 9A ). In addition, LinkedOmics was used to positively and negatively of the top 50 important genes associated with BIRC5 co-expression in LUAD ( Figures 9B, C ). Significant Gene Ontology (GO) term annotation indicated that BIRC5 correlated genes mainly participated in organelle fission, ribonucleoprotein complex biogenesis and mitotic cell cycle phase transition ( Figure 9D ). The KEGG pathway showed that BIRC5 related genes were primarily associated with mitotic cell cycle phase transition, double-strand break repair and negative regulation of cell cycle process ( Figure 9E ). Moreover, the significant top four genes were considered as the hub genes, namely CENPA (cor=9.214e-01, p=1.776e-158), CDC20 (cor=9.098e-01, p=1.422e-147), CCNB2 (cor=9.012e-01, p=2.107e-140), and NUF2 (cor=8.971e-01, p=3.370e-137); these were positively correlated with BIRC5 in LUAD ( Figure 9F ).

Due to the significance of BIRC5 in LUAD, we further explored BIRC5 networks of kinase, miRNA or transcription factor targets in LUAD. For kinase networks of BIRC5, the top most significant targets were involved in the kinase AURKB. We found that the top 5 most significant miRNA network targets were MIR-17-5P, MIR-507, MIR-24, MIR-138 and MIR-199A ( Table 2 ). Moreover, the transcription factor network target (V$E2F1_Q6, V$E2F4DP1_01, V$E2F_Q6, V$E2F_Q4, V$E2F_Q4_01) of BIRC5 was shown in Table 2 .

To determine whether BIRC5 functions in the drug resistance of cancer, GeneMANIA was employed. BIRC5 showed interactions with 20 proteins/genes. And the results showed that these genes were primarily involved in the regulation of microtubule associated complex, spindle, spindle assembly, chromosome, centromeric region, mitotic nuclear division, chromosomal region and mitotic spindle organization ( Figure 10A ).

To further explore the regulators of BIRC5 in expression, we performed an analysis of kinases, miRNAs and transcription factors of BIRC5 gene ( Table 2 ). For kinase networks of BIRC5, only one kinase target of BIRC5 was identified (Kinase_AURKB) according to the LinkedOmics database. To probe the function of AURKB in more details, we next constructed PPI networks of AURKB. And the results showed that these genes were primarily involved in the regulation of microtubule associated complex, spindle, spindle assembly and mitotic nuclear division ( Figure 10B ).

BIRC5 and the top four significant genes, namely CENPA, CDC20, CCNB2, and NUF2, were selected as the hub genes for genetic alteration, pathway, and drug sensitivity analyses. We also explored the role of hub genes in all known cancer related pathways, including Apoptosis, Cell Cycle, DNA Damage Response, EMT, Hormone AR, Hormone ER, PI3K/AKT, RAS/MAPK, RTK, and TSC/mTOR pathways. Results showed that BIRC5 is involved in the activation of Apoptosis, Cell Cycle, DNA Damage Response, EMT, and Hormone AR pathways ( Figure 11 ).

By drug sensitivity analysis, we evaluated the correlation of gene expression and IC50 using Genomics of Drug Sensitivity in Cancer (GDSC) of small molecules. As shown in Figure 12 , a high NUF2 level was resistant to 5 drugs or small molecules, whereas a high CENPA level was resistant to 2 drugs or small molecules. Moreover, high expression of BIRC5 was associated with resistance to trametinib, RDEA119, and selumetinib. These findings may provide support for the selection of drugs for targeted therapy in LUAD.

In the present study, stepwise exposure to increasing concentrations of afatinib was used to establish afatinib-resistant lung cancer cells. The half maximal inhibitory concentration (IC50) of afatinib treated for 72 h in HCC827 and HCC827-AR cells was measured by CCK8 ( Figure 13A ). To explore the role of BIRC5 in NSCLC cell resistance to afatinib, western blot analysis was performed. As shown in Figure 13B , BIRC5 expression was significantly higher in HCC827-AR cells than in HCC827 cells (P<0.001). Next, we tested the sensitivity of HCC827-AR cells to afatinib and BIRC5 inhibitor YM155 in vitro. HCC827-AR has a half-maximal inhibitory concentration of afatinib (IC50)=1–10 μM), YM155 (IC 50 = 1–10 nM) ( Figure 13C ). The results confirmed that YM155 is about 1000 times more potent than afatinib in vitro in afatinib-resistant cells. This finding provides potential new therapeutic strategies against lung adenocarcinoma cells with acquired drug resistance.