The training set consisted of gene expression data (HTSeq-FPKM) and associated clinical information from the TCGA database, including 458 LUAD patients with complete follow-up data and a follow-up duration longer than 30 days. Patients lacking comprehensive survival information were excluded. The validation set followed identical inclusion criteria and incorporated datasets GSE68465 (439 LUAD patients) and GSE72094 (386 LUAD patients) from the GEO database. Data from the TCGA-LUAD, GSE68465, and GSE72094 datasets were harmonized using the “Combat” algorithm from the R package “sva” (16). We have supplemented detailed steps and parameter settings for integrating data from different platforms using the “Combat” algorithm to ensure transparency and reproducibility of the data integration process. Specifically, we have supplemented detailed steps and parameter settings for integrating data from different platforms using the “Combat” algorithm. Specifically, to eliminate batch effects between the TCGA-LUAD, GSE68465, and GSE72094 datasets, we employed the “Combat” algorithm from the R package “sva” for data integration. First, we performed log2 transformation on the gene expression matrices of the three datasets to obtain approximately normally distributed expression values. Subsequently, we identified and retained genes common to all three datasets to construct a merged expression matrix containing all samples. When using the ComBat function, we set the parameter par.prior=TRUE to enable the empirical Bayesian method for enhancing the robustness of batch effect correction, while setting mean.only=FALSE to adjust both location and scale parameters simultaneously. To preserve biologically significant variation, we incorporated important clinical covariates (such as age, gender, and tumor stage) using the mod parameter, ensuring that the influence of these variables was retained during the batch correction process. After batch correction, we verified the effective removal of batch effects using principal component analysis (PCA) and t-SNE visualization methods, confirming that samples no longer clustered according to their original datasets but rather distributed based on their biological characteristics.

To further investigate the mechanisms of erlotinib resistance (ER), we incorporated resistance data from four GEO datasets of LUAD cell lines:

GSE75308-PC9: This dataset contains gene expression profiles of PC9 cell line and its erlotinib-resistant derivatives, generated using the Illumina HumanHT-12 V4.0 expression microarray platform. It comprises 12 samples (6 sensitive and 6 resistant strains).

GSE67051-PC9: This dataset includes gene expression data of PC9 cell line at different time points of erlotinib exposure, using the Affymetrix Human Gene 1.0 ST array platform. It contains 15 samples, documenting the dynamic process of resistance acquisition.

GSE67051-HCC827: Expression data of HCC827 cell line and its resistant derivatives from the same study, using the identical Affymetrix platform, comprising 12 samples.

GSE123031-PC9: This dataset contains RNA-seq data of PC9 and various acquired resistant PC9 variant cell lines, generated using the Illumina HiSeq 2500 platform, comprising 9 samples.

These cell line datasets were downloaded from the GEO database and processed using the same standardization procedures as the primary clinical datasets. For microarray data, we employed the RMA algorithm for background correction and normalization; for RNA-seq data, we applied the same HTSeq-FPKM normalization method used for TCGA data. Prior to integrative analysis, all cell line datasets underwent batch effect correction using the Combat method as described above.

A comprehensive flowchart of the study is depicted in Figure 1 .

Tissue samples were obtained from the Fourth Affiliated Hospital of Soochow University. Clinical data for these patients were retrieved from their medical records. This study was approved by the Ethics Committee of the same institution, and informed consent was obtained from all participants. Additionally, the study adhered to the ethical principles outlined in the Declaration of Helsinki, as published in the British Medical Journal (July 18, 1964).

The PC9 and THP-1 cell lines were acquired from the Shanghai Institutes for Biological Sciences (Shanghai, China). An erlotinib-resistant PC9 cell line (PC9-ER) was developed over six months by progressively increasing the erlotinib concentration from 0.1 μM to 40 μM using a stepwise incremental method. Both PC9 and PC9-ER cells were cultured in DMEM supplemented with 10% fetal bovine serum and antibiotics, maintained in a humidified atmosphere with 5% CO₂ at 37°C. THP-1 cell line was maintained in RPMI 1640 medium supplemented with 10% heat-inactivated fetal bovine serum, 2 mM L-glutamine, and 1% penicillin/streptomycin in a humidified incubator at 37°C with 5% CO₂.

Calcium-related genes were sourced from GeneCards, selecting those with a relevance score ≥ 8. Univariate Cox survival analysis and Kaplan-Meier (KM) survival analysis were performed on the TCGA-LUAD, GSE68465, and GSE72094 datasets using the “survival” package. A stringent P-value threshold of < 0.001 was applied to ensure result reliability, leading to the identification of 33 CRPGs.

Unsupervised clustering was conducted using the “ConsensusClusterPlus” R software package to identify novel subtypes (17). Subtypes were analyzed using the “survival” and “survminer” packages, employing KM survival analysis and log-rank tests. Mesenchymal and immune cell presence in malignant tissues was estimated using the ESTIMATE algorithm, with tumor purity calculated for different molecular subpopulations (18). Immune cell abundance across molecular subpopulations was assessed using MCP-Counter (19). The GSVA algorithm was employed to investigate significant TME characteristic differences (20, 21). Additionally, HLA gene expression levels were compared between subtypes.

LASSO-Cox regression analysis was utilized to refine the 33 CRPGs, eliminating collinearity with TCGA LUAD patients serving as the training set. This process led to the development of a calcium-related risk score prognostic signature, calculated by multiplying the β (Coef) values by the expression levels of CRPGs. The risk score formula used was: Risk score = (β1*CRPG1 + β2*CRPG2 + β3*CRPG3 +… + βn*CRPGn), where β signifies the CRPG coefficient (22, 23). The median risk score was used as a cut-off threshold to classify patients into high- and low-score groups. KM survival analyses were carried out to evaluate these groups further. The accuracy of the risk model was measured by calculating the area under the curve (AUC) for predicting 1-, 2-, and 3-year OS outcomes.

Validation sets GSE68465 and GSE72094 were similarly analyzed. Univariate and multivariate regressions, visualized through forest plots, as well as KM survival analyses, were conducted to verify the validity of the signatures. Stratified analyses of clinical variables, including age, gender, stage, and T and N classifications, were performed to compare OS between the high- and low-risk groups.

To assess immune infiltration in patients, seven algorithms were employed: MCPcounter (19), CIBERSORT (24), xCell (25), TIMER (26), EPIC (27), Cibersort-ABS (28), and QUANTISEQ (29). These algorithms were used to compare the high-risk and low-risk groups, identifying differences in immune infiltration. Pearson correlations determined the relationship between risk scores and immune cell content. Normalized enrichment scores (NES) were calculated using the Hallmark gene set with the ‘GSEA’ package to compute NES and false discovery rate (FDR). Gene set enrichment analysis (GSEA) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were performed on the Hallmark gene set for both high-risk and low-risk groups.

Drug sensitivity analysis was carried out using the GSCA database, based on gene expression data (30). GSCALite provided a collection of 750 small molecule drugs from GDSC and CTRP, using gene expression data to identify differentially expressed genes associated with these drugs in samples with various risks of drug response. Among these differential genes, the top 150 that were significantly up-regulated or down-regulated were selected as signature-related markers. The CMap_gene_signatures. RData file, obtained from the database website (https://www.pmgenomics.ca/bhklab/sites/default/files/downloads), containing 1,288 compound-related features, was used to calculate the matching score. The analysis methodology adhered to those described in previous publications (31, 32).

We applied the train function from the caret package to train multiple machine learning models and used the explain function from the DALEX package to interpret these models. The predict function was used to evaluate the accuracy of each model and generate ROC curves, while the variable_importance function from the DALEX package was employed to calculate the importance of variables within the models. Reverse cumulative residual distribution plots and residual boxplots were generated to assess model performance. The top 10 most important genes for each model (as shown in the variable importance ranking plot) were identified. The root mean square error (RMSE) was used as the loss function, indicating the degree to which excluding a variable impacts the predicted values of the response variable, where larger RMSE values imply greater variable importance.

MR analyses were conducted using the R package TwoSampleMR V0.5.6.24 (33). LUAD data, identified by code ieu-a-984, were sourced from the Open GWAS IEU website (https://gwas.mrcieu.ac.uk/) as the outcome, while cis-eQTL data from eQTLGen served as the exposure (34). Screening criteria for the exposure file included a p-value threshold of 5e-8, a correlation coefficient of 0.001, a distance of 10,000 kb, and a minor allele frequency greater than 0.01, focusing on cis-loci within 1 MB upstream and downstream of the gene center. After data loading, outcome harmonization was performed using built-in functions. MR estimates for each SNP were calculated using the Wald ratio method, while the inverse-variance weighted (IVW) method was applied when multiple SNPs were present to compute a weighted mean of the ratio estimates, with weights based on the inverse variance of the ratio estimates. Horizontal pleiotropy was assessed by evaluating whether the MR Egger intercept significantly differed from 0 when at least three SNPs were analyzed. Cochran’s Q tests assessed heterogeneity among Wald ratios (35). Significance was determined using FDR-corrected p-values, with an FDR threshold set at less than 0.05. Statistically significant findings in replication studies were recognized when nominal p-values were below 0.05.

To assess the horizontal pleiotropy of potential drug targets and possible side effects, we performed PheWAS using the AstraZeneca PheWAS Portal (https://azphewas.com/). Multiple testing corrections were applied, and a significance threshold of 1E-8 was established, as recommended by the AstraZeneca PheWAS Portal, to minimize the risk of false positives (36).

Experiments were executed following the manufacturer’s protocols. Total RNA from cells or tissues was extracted using the TRIzol reagent (Invitrogen, Carlsbad, CA). Reverse transcription was conducted with the PrimeScript RT kit (Takara, Cat: RR036A, KeyGEN) to synthesize cDNA, which was subsequently subjected to qRT-PCR analysis. Each qRT-PCR trial was performed in triplicate, and results were normalized to β-actin using the 2^-ΔΔCt method. The primer sequences utilized for qRT-PCR analysis were as follows: ACTIN forward, 5′-GTCATTCCAAATATGAGATGCGT-3′; ACTIN reverse, 5′-GCATTACATAATTTACACGAAAGCA-3′; TNNC1 forward, 5′-GTCTGACCTCTTCCGCATGT-3′; TNNC1 reverse, 5′-ATGAGCTCCTCGATGTCGTC-3′.

Total protein was extracted from the cells using RIPA lysis buffer containing protease inhibitors. Protein concentration was determined using a BCA Protein Assay Kit. Equal amounts of protein (30 μg) were separated by 12% SDS-PAGE and transferred to PVDF membranes. The membranes were blocked with 5% non-fat milk in TBST for 1 hour at room temperature and then incubated with primary antibodies against TNNC1 (1:1000) and HA-tag (1:2000) overnight at 4°C. After washing with TBST, the membranes were incubated with HRP-conjugated secondary antibodies (1:5000) for 1 hour at room temperature. Protein bands were visualized using an enhanced chemiluminescence detection system. GAPDH served as a loading control.

Small interfering RNA (siRNA) targeting TNNC1 (RiboBio, Guangzhou, China) was transfected into PC9 or PC9-ER cells using Lipofectamine 3000 (Invitrogen, Carlsbad, CA, USA), following the manufacturer’s instructions. Forty-eight hours post-transfection, cells were harvested for qRT-PCR and other subsequent experiments. The siRNA sequences were as follows: si1-TNNC1 sense sequence, 5′-CGGUAGAGCAGCUGACAGA-3′; si2-TNNC1 antisense sequence, 5′-AAGAUAAUGCUGCAGGCUA-3′; si-NC sense sequence, 5′-UAACGACGCGACGACGUAAtt-3′; si-NC antisense sequence, 5′-UUACGUCGUCGCGUCGUUAtt-3′.

The TNNC1 overexpression system was constructed using lentiviral vectors. The full-length coding sequence of human TNNC1 was amplified by PCR and cloned into a lentiviral expression vector with an HA-tag. The recombinant plasmid was verified by DNA sequencing. Empty vector was used as a control. Lentiviral particles were produced in HEK293T cells by co-transfection with packaging plasmids. PC9 cells were infected with lentivirus containing either TNNC1-HA or empty vector control, and stable cell lines were established by puromycin selection (2 μg/mL) for 2 weeks.

Cells were seeded at a density of 4 × 10² cells per well in a 6-well plate. After overnight stabilization, cells were treated with erlotinib and incubated at 37°C until colonies formed and matured, with media changes every three days. After 10 days, colonies were fixed with 4% paraformaldehyde for 10 minutes, stained with crystal violet for 30 minutes, washed with PBS, and imaged.

PC9 and PC9-ER cells transfected with si-TNNC1 or si-NC were seeded in 96-well plates at a density of 2 × 10³ cells per well and incubated for 6–8 hours to establish baseline values. They were then exposed to various concentrations of erlotinib for 72 hours. Following treatment, 10 μl of CCK-8 solution was added to each well and incubated for 2 hours. Cell viability was assessed by measuring absorbance at 450 nm. For M2 macrophage polarization, PMA-differentiated THP-1 cells (M0) were stimulated with 20 ng/ml IL-4 and 20 ng/ml IL-13 for 48 hours.

THP-1 cells were differentiated into macrophages using phorbol 12-myristate 13-acetate (PMA). The differentiation of THP-1 monocytes into macrophage-like cells was induced by PMA. Briefly, THP-1 cells were seeded at a density of 5 × 10⁵ cells/ml in complete medium containing 50 ng/ml PMA for 24 hours to induce initial differentiation into M0 macrophages. After PMA treatment, the medium was replaced with fresh complete medium without PMA, and cells were rested for an additional 24 hours. For M2 macrophage polarization, PMA-differentiated THP-1 cells (M0) were stimulated with 20 ng/ml IL-4 and 20 ng/ml IL-13 for 48 hours.

Differentiated cells were washed twice with PBS and then incubated with PBS containing 2 mM EDTA at 37°C for 10 minutes to detach cells from culture plates. After collection, cells were centrifuged at 300 × g for 5 minutes, and the supernatant was discarded. Cells were resuspended in 100 μl of flow cytometry buffer (PBS containing 0.5% BSA and 2 mM EDTA). Fc receptor blocking reagent was added and incubated at 4°C for 15 minutes to reduce non-specific binding. Fluorochrome-conjugated antibodies were then added: anti-human CD206-PE (1:50 dilution) and anti-human CD163-APC (1:50 dilution), and incubated at 4°C in the dark for 30 minutes. Following incubation, cells were washed twice with flow cytometry buffer and resuspended in 300 μl of flow cytometry buffer. Cell analysis was performed using a flow cytometer (BD FACSCalibur), with a minimum of 10,000 events collected per sample. Data were analyzed using FlowJo software (Tree Star Inc.), and M2 macrophage polarization was assessed by measuring the percentage of CD206 and CD163 positive cells and their mean fluorescence intensity (MFI). Appropriate isotype control antibodies were used as negative controls.

Data analysis was conducted using R4.1 software. T-tests were employed for normally distributed data, while the Wilcoxon rank-sum test was used for non-normally distributed data. The Kruskal-Wallis test and one-way analysis of variance (ANOVA) served as non-parametric and parametric methods, respectively, for intergroup comparisons. Univariate Cox survival analysis and Kaplan-Meier survival analysis were conducted using the “survival” package. Data analysis and graphical representation were facilitated by GraphPad Prism 9. All experiments were independently repeated at least three times. A statistical significance level of P < 0.05 was considered.

A comprehensive flowchart of the study is depicted in Figure 1 .

To ensure high sample quality, a rigorous screening process was implemented, selecting only samples with complete survival data and a follow-up duration exceeding 30 days. The R package “sva” and the “Combat” algorithm were utilized to minimize non-biological variations across the three datasets, culminating in the inclusion of 1,283 patients (TCGA-LUAD=458, GSE68465 = 439, GSE72094 = 386). Univariate Cox survival analysis (P < 0.001) and KM survival analysis (P < 0.001) conducted on 392 calcium-related genes with a relevance score of ≥ 8 led to the identification of 33 CRPGs ( Figure 2A ; Supplementary Table S1 ).

Unsupervised clustering was applied to 1,283 LUAD tissues to identify novel subtypes based on the expression profiles of 33 CRPGs. Through rigorous analysis using cluster consensus, a consensus matrix, and a delta area plot, four clusters (C1, C2, C3, C4) were identified ( Figures 2B-D ; Supplementary Table S2 ). Principal component analysis (PCA) indicated significant transcriptional variations of CRPGs across these clusters ( Figure 2E ). All 33 CRPGs showed markedly different expression levels among the clusters (P < 0.0001) ( Figure 2F ). KM survival analyses and log-rank tests revealed substantial survival differences between most cluster pairs (P < 0.005, Figure 2G ). The Chi-Squared test demonstrated diverse survival rates among clusters, with C3 exhibiting the highest survival rate at 73%, followed by C1 and C4 (58% and 63%, respectively), while C2 had the lowest survival rate at 45% (P = 5.57E-11, Figure 2H ).

Further analysis of the subtypes revealed that C3 had the highest HLA gene expression and lowest tumor purity, indicative of stronger immune responses, whereas C2 displayed the highest tumor purity and the lowest immune cell infiltration, correlating with its poor prognosis ( Figures 3A-D ). GSVA and limma differential analysis uncovered significant pathway variations among the subtypes: C3 was enriched in immune response and T cell activation pathways, while C2 was characterized by enrichment in cell cycle and DNA replication pathways ( Figures 3E-J ). These differences in molecular characteristics likely explain the observed prognostic disparities and provide insight into potential precision treatment strategies for each subtype.

These findings emphasize the role of the 33 CRPGs in differentiating patients into four biologically distinct clusters, underscoring their potential prognostic and therapeutic significance. The molecular characteristic differences revealed by subtype analysis provided the basis for our prognostic model construction. Next, we developed a prognostic scoring model based on these 33 CRPGs and validated its predictive value in multiple independent cohorts.

Utilizing the TCGA-LUAD dataset as the training set, we applied LASSO-Cox regression analysis to streamline the list of 33 CRPGs, removing collinearity and optimizing prognostic characteristics ( Figures 4A, B ). The process led to the development of a risk score prognostic signature, calculated by multiplying the β (Coef) values by the CRPG expression levels, resulting in the formula: CRPGscore = (0.034 × S100P expression) + (0.109 × S100A10 expression) + (-0.031 × SLC34A2 expression) + (0.050 × CAMK2N1 expression) + (0.160 × PKP2 expression) + (-0.038 × CFTR expression) + (-0.006 × TNNC1 expression) + (0.113 × MAPK10 expression) + (0.132 × CDKN3 expression) ( Figure 4C ; Supplementary Table S3 ).

The biological roles of each CRPG within the signature further underscore their importance in shaping the prognostic model. For instance, S100P and S100A10, both calcium-binding proteins, are involved in regulating cell proliferation, differentiation, and migration. S100P is overexpressed in multiple cancers and linked to poor prognosis in LUAD, while S100A10 contributes to membrane repair and drug resistance in tumor progression. Contrasting with these findings, SLC34A2, a sodium-phosphate cotransporter downregulated in LUAD, functions as a tumor suppressor gene by maintaining calcium-phosphate balance. Similarly, downregulated TNNC1, a novel finding in LUAD, is associated with extended survival, suggesting potential tumor-suppressive functions. On the other hand, overexpressed genes like PKP2 and CDKN3 may promote tumor progression by affecting calcium-dependent cell-cell adhesion and cell cycle regulation, respectively. The inclusion of CAMK2N1 and MAPK10, involved in calcium-driven stress response and apoptosis, and CFTR, an immune microenvironment regulator, highlights the multifaceted roles of calcium signaling in LUAD progression. These insights provide a mechanistic basis for the prognostic utility of the CRPG signature.

To assess the predictive value of the signature, a random effects model was employed to conduct a prognostic meta-analysis across the GEO cohorts (GSE72094 and GSE68465) and the TCGA-LUAD cohort. This analysis indicated that the prognostic signature is a significant risk factor across datasets (HR: 3.47, 95% CI: 2.34-5.15, Weight = 100.0%). No significant heterogeneity was observed (I² = 64%, t = 0.0751, P = 0.06) ( Figure 4D ). Using the median risk score from the TCGA-LUAD scores as the cut-off value, patients across all three cohorts were classified into high- and low-score groups. The low-score group consistently showed significantly better OS compared to the high-score group in TCGA-LUAD (P < 0.001, Figure 4E ), GSE72094 (P < 0.001, Figure 4F ), and GSE68465 (P = 0.001, Figure 4G ). ROC analysis was performed to evaluate the sensitivity and specificity of the signature for prognosis, using the AUC as the performance metric. The AUC values for the TCGA-LUAD training set at 1, 2, and 3 years were 0.713, 0.684, and 0.696, respectively ( Figure 4H ). For the GSE72094 test set, the AUC values were 0.699, 0.694, and 0.711, whereas the GSE68465 test set yielded AUC values of 0.699, 0.709, and 0.691 ( Figures 4I, J ). Moreover, the reliability of the ROC in predicting outcomes in LUAD patients was confirmed ( Figures 4K-M ).

These results highlight the robustness and reliability of the constructed prognostic signature in accurately predicting the prognosis of LUAD across multiple datasets.

Based on the aforementioned analysis of the expression pattern and functional characterization of CRPGs in LUAD, we further evaluated the value of the constructed prognostic features for application in different clinical subgroups. Cox regression analysis was performed to evaluate the signature as an independent predictor in LUAD patients. The univariate Cox regression analysis showed that OS in LUAD was strongly associated with stage, T stage, N stage, and risk score ( Figure 5A ). Multivariate Cox regression analysis further confirmed that the risk score was a significant independent prognostic factor for LUAD patients (P < 0.001) ( Figure 5B ). An examination of the relationship between the signature and clinicopathological variables revealed that the CRPGscore was significantly higher in patients with stage III-IV compared to stage I-II (P = 0.0089, Figure 5C ), higher in N1–2 compared to N0 (P = 0.0015, Figure 5D ), and higher in males compared to females (P = 0.0099, Figure 5E ). These findings suggest that a higher CRPGscore is linked with increased malignancy in LUAD. Stratification of LUAD patients based on clinicopathological variables was conducted to investigate the prognostic value for OS. Stratification analyses were performed based on age (P < 0.001; P < 0.001; Figures 5F, G ), sex (P < 0.001; P < 0.001; Figures 5H, I ), N stage (P < 0.001; P = 0.005; Figures 5J, K ), stage (P < 0.001; P = 0.002; Figures 5L, M ), and T stage (P < 0.001; P = 0.026; Figures 5N, O ). Consistently across all stratified analyses, the high-score group exhibited significantly shorter OS than the low-score group, indicating that the signature’s prognostic value was not confounded by conventional clinical factors in LUAD patients.

We investigated the relationship between our gene signature and the tumor immune microenvironment in LUAD using seven algorithms (MCPcounter, CIBERSORT, xCell, TIMER, EPIC, CIBERSORT-ABS, and QUANTISEQ). The signature exhibited significant associations with various immune cell types, such as Tregs, T cells, NK cells, neutrophils, myeloid dendritic cells, monocytes, mast cells, and B cells, with notably higher infiltrating levels in the low-risk group than in the high-risk group ( Figures 6A, B ). Most CRPGs (except PKP2) were significantly correlated with immune-related scores, and S100P showed a strong positive relationship with TumorPurity ( Figure 6C ). Single-cell analyses further demonstrated enriched S100P expression in tumor-associated cells ( Figure 6D ).

Hallmark gene set-based GSEA and KEGG analyses revealed distinct pathway enrichments among high- and low-risk patients, indicating differing biological features ( Supplementary Figure S1 ).

To explore the clinical relevance of the signature in anti-cancer therapy, we performed Connectivity Map (CMap) analysis using the top 150 differentially expressed genes in the high-risk group. Five candidate compounds (TTNPB, W.13, iloprost, NS.398, and Gly.His.Lys) were identified for their possible inhibitory effects on tumor-promoting mechanisms ( Figure 7A ). Notably, previous research has documented the inhibitory effects of NS.398 (37), TTNPB (38) and iloprost (39) on LUAD. Further integration of the GDSC and CTRP databases revealed that TNNC1, S100A10, PKP2, CAMK2N1, and S100P exhibited positive correlations with chemotherapy resistance ( Figures 7B, C ). Notably, TNNC1, S100A10, PKP2, CAMK2N1, S100P, and CDKN3 in the CTRP dataset were linked to resistance across multiple drugs ( Figure 7C ). Collectively, these results underscore the potential clinical utility of these genes for predicting treatment outcomes.

The TCGA-LUAD dataset was randomly divided into two subsets using the createDataPartition function from the caret package, designating one as the training set and the other as the test set. Diagnostic models were developed using nine CRPGs. Residual analysis demonstrated strong predictive capabilities for these models ( Figure 8A ), with minimal prediction errors evident in the inverse cumulative distribution of residuals across 11 methods ( Figure 8B ). ROC curves based on these residuals indicated excellent model performance ( Figure 8C ). Among the nine CRPGs, TNNC1 consistently emerged as a critical tumor suppressor gene in LUAD across all algorithmic models ( Figure 8D ).

Building on this discovery, MR analysis employing 33 CRPGs further confirmed the tumor-suppressive role of TNNC1. The analysis revealed that increased TNNC1 expression significantly reduces the risk of LUAD (OR = 0.737; 95% CI, 0.550–0.988; P = 0.041) ( Supplementary Figure S2A ). Additionally, a PheWAS analysis conducted via the Phenoscanner platform simulated potential drug-related side effects associated with TNNC1, and no significant adverse effects were observed ( Supplementary Figure S2B and Supplementary Table S4 ). These results not only reinforce TNNC1’s protective role in LUAD but also highlight its potential as a novel therapeutic target.

Subsequent validation with the TCGA-LUAD dataset confirmed that TNNC1 expression was significantly lower in LUAD tissues compared to adjacent normal lung tissues ( Figures 8E, F , P < 0.001). ROC curves supported TNNC1’s robust differentiation between LUAD and normal lung tissues (AUC: 0.997, 95% CI: 0.994-0.999) ( Figure 8G ), and the Hosmer-Lemeshow test confirmed the model’s goodness of fit (P = 0.998, Figure 8H ). Further validation using external datasets GSE31547 and GSE40791 affirmed that TNNC1 expression is significantly lower in LUAD tissues compared to adjacent normal tissues ( Figures 8I, J , P < 0.001). Kaplan-Meier survival analysis showed that patients with higher TNNC1 expression had significantly better overall survival (P = 0.006, Figure 8K ), progression-free interval (P = 0.015, Figure 8L ), and disease-free interval (P = 0.021, Figure 8M ) than those with lower expression. These findings highlight TNNC1’s downregulation in LUAD tissues and its close association with patient prognosis, suggesting its potential role as a prognostic biomarker.

Resistance to EGFR-TKIs is a common challenge in treating late-stage LUAD, often leading to relapse (40). Understanding the mechanisms behind this resistance is crucial for improving chemotherapy outcomes. TCGA-LUAD patient data revealed that TNNC1 expression was significantly higher in patients classified in the progressive disease (PD) category compared to those with a complete response (CR) to EGFR-TKI treatment ( Figure 9A ). Further, an examination of erlotinib resistance (ER) datasets from LUAD cell lines (GSE75308-PC9, GSE67051-PC9, GSE67051-HCC827, GSE123031-PC9) consistently showed an increase in TNNC1 expression in erlotinib-resistant LUAD cells compared to sensitive ones ( Figure 9B ). Elevated TNNC1 expression in the ER group, as opposed to patients untreated with erlotinib, underscores its potential role as a key gene related to erlotinib resistance ( Figures 9C, D ). To explore TNNC1’s impact on erlotinib resistance, a PC9-ER cell model was developed using a dose-escalation approach with the PC9 cell line. The half-maximal inhibitory concentration (IC50) of erlotinib required for a 50% reduction in cell viability was assessed using the CCK-8 assay and a dose-response curve ( Figure 9E ). TNNC1 mRNA levels were found to be higher in PC9-ER cells compared to PC9 cells ( Figure 9F ). To investigate further, PC9-ER cells were transfected with TNNC1-specific siRNA (si1/2-TNNC1) to knock down TNNC1 expression, and transfection efficiency was confirmed 48 hours post-transfection via RT-qPCR ( Figure 9G ). Following TNNC1 knockdown, the IC50 value of erlotinib in PC9-ER cells decreased compared with the control ( Figure 9H ). The CCK-8 assay results showed that downregulation of TNNC1 significantly reduced the IC₅₀ value of erlotinib in these cells ( Figure 9I ). In addition, TNNC1 downregulation markedly inhibited PC9-ER cell proliferation as indicated by CCK-8 proliferation and colony formation assays ( Figures 9J, K ). To establish a direct causal relationship between TNNC1 overexpression and erlotinib resistance, we overexpressed TNNC1 in the erlotinib-sensitive PC9 cell line. Western blot confirmed successful expression of HA-tagged TNNC1 in PC9 cells ( Supplementary Figure S3A ). Cell viability assays demonstrated that TNNC1-overexpressing PC9 cells treated with 5μM erlotinib exhibited significantly higher viability at day 3 compared to vector control cells under the same treatment (P < 0.05) ( Supplementary Figure S3B ). These findings collectively suggest that TNNC1 enhances LUAD cell resistance to erlotinib.

Tumor-associated macrophage M2 polarization is implicated in tumor progression and chemotherapy resistance (41). This study assessed whether TNNC1 affects LUAD biological functions through M2 tumor-associated macrophage (TAM) polarization. xCell analysis showed that patients with high TNNC1 expression experienced more pronounced immune cell infiltration than those with low expression ( Figure 9L ), including significant enrichment of total macrophages and M2 polarization ( Figure 9M ). Immune-related genes, such as immune stimulatory and inhibitory genes, chemokines, and human leukocyte antigens were generally upregulated in the high TNNC1 expression group ( Figure 9N ). These data suggest a strong association between TNNC1 and M2 TAM polarization. The relationship between TNNC1 expression in M2 cells and LUAD patient survival was explored. Using the TIMER2.0 database to integrate M2 cell infiltration levels and TNNC1 transcriptional expression, LUAD patients were categorized into four groups. Results indicated that patients with a TNNC1+ & M2+ status had significantly poorer survival outcomes ( Figures 9O, P ). It was hypothesized that TNNC1 could influence macrophage M2 polarization. To evaluate this, THP-1 monocytes were treated with IL-3 to induce M2 polarization, then transfected with TNNC1-specific siRNA. IL-3 treatment raised levels of M2 markers, including CD206, Ym1, and Arg1, while TNNC1 knockdown reversed these increases ( Figure 9Q ). To directly confirm the role of TNNC1 in promoting macrophage M2 polarization, we overexpressed TNNC1 in THP1 cells and assessed its effects on M2 phenotype markers. Western blot confirmed successful expression of HA-tagged TNNC1 in THP1 cells ( Supplementary Figure S1C ). Flow cytometry analysis demonstrated that TNNC1 overexpression significantly increased the percentage of CD163+/CD206+ THP1 cells under IL-4/IL-13 stimulation compared to control cells (P < 0.01) ( Supplementary Figure S1C ). qPCR analysis further revealed that TNNC1 overexpression significantly upregulated mRNA expression of M2-related markers including Arg1 (P < 0.001), IL-10, TGF-β, CCL17, and CCL22 (all P < 0.05) ( Supplementary Figure S1D ). Taken together, these findings provide direct evidence that TNNC1 positively regulates macrophage M2 polarization.

In conclusion, TNNC1 may contribute to erlotinib resistance in LUAD cells by promoting macrophage M2 polarization.