The R package “TCGAbiolinks” was used to download the log-transformed FPKM expression profiles and clinical information from the TCGA-LUAD dataset. A total of 497 tumor samples with both expression data and survival information were retained for the construction of the prognostic signature. The GSE41271 and GSE42127 bulk expression datasets were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/) and were used to validate the prognostic signature. The data processing standard of the GEO bulk expression dataset was as follows: the probes were converted to gene symbols according to the probe correspondence with the platform. If one probe corresponded to multiple genes, the probe was removed, and if multiple probes corresponded to the same symbol, the median value was taken. Single-cell RNA-seq data were obtained from the GSE131907 dataset of the GEO database. Fifteen primary LUAD samples from the GSE131907 dataset were used for the analyses. The clinical and transcriptomic data of the GSE126044 and GSE135222 cohorts, in which NSCLC patients were treated with the PD-1/PD-L1 blockade, were downloaded from the GEO database and used to evaluate the predictive efficacy of the prognostic signature. A total of 480 lung T_RM cell-specific genes were obtained from a previous publication (11).

The R package “limma” was used to identify the differentially expressed genes (DEGs) between LUAD and adjacent normal tissues, with thresholds set at |log2FC|≥1 and FDR<0.05. The DEGs intersecting with the lung T_RM cell-specific genes were regarded as differentially expressed (DE) T_RM cell-related genes and were chosen for subsequent analysis.

The interactive relationships of the lung T_RM cell-related genes were acquired from the STRING database (https://www.string-db.org/), and a protein-protein interaction (PPI) was constructed based on this information.

The R package “clusterProfiler” was applied for the functional annotation of the DE lung T_RM cell-related genes, with the p-value cutoff set at 0.05. Functional enrichment analysis was performed to predict the potential biological functions of these genes.

Univariate Cox regression analysis was used to determine the hazard ratios (HR) and prognostic significance. Genes with p values < 0.05 were prognosis-associated genes. Least Absolute Shrinkage and Selection Operator (LASSO) regression analysis was applied to further identify key prognostic factors, and a risk score model for predicting survival was constructed by weighting the expression of each key prognostic gene with LASSO regression coefficients (“exp” represents the expression level of the genes, and “coef” represents the Cox regression coefficient):

The patients were divided into high-risk and low-risk groups based on the median risk score. The “Kaplan-Meier” method was used to generate survival curves for the prognostic analysis, and the “log-rank” test was used to evaluate the significance of the differences in overall survival between groups. The receiver operating characteristic (ROC) curve was used to assess the predictive efficacy of the prognostic models. The R package “timeROC” was used to visualize the “area under the curve” (AUC). Univariate and multivariate Cox regression analyses were performed to evaluate the independent predictive value of the prognostic model.

The “ESTIMATE” algorithm was used to calculate the immunity score, stroma score, and tumor purity for each tumor sample, and then the “Wilcoxon” test was subsequently used to compare the differences in the immunity score, stroma score, and tumor purity among different subgroups of samples. The correlations between the risk score and the immunity score, stroma score, and tumor purity were calculated using Spearman analysis. Single-sample gene set enrichment analysis (ssGSEA) was used to evaluate the relative abundance of each infiltrating cell in the TME. The gene sets of the 28 types of immune cells used in the analysis were obtained from a previous publication (21). The R packages “GSVA” and “GSEABase” were used to compare the differences in biological pathways and immune functions.

The “calcPhenotype” function of the R package “oncoPredict” was used to assess the IC50 values of the samples for the drugs. The correlation coefficients between the risk score, the expression of genes included in the prognostic model, and the drug IC50 values were calculated using Spearman analysis.

The R package “Seurat” (version 4.1.0) was used for quality control of the scRNA-seq data. To exclude some low-quality cells and genes expressed at low levels, we set the thresholds as follows (1): each gene was expressed in at least three cells (2); the number of features per cell was between 500 and 6000, and the number of counts per cell was between 1000 and 20,000; and (3) the number of mitochondrial and erythrocyte genes was less than 20% of the total number of genes in each cell. Next, the “NormalizeData” function was used for normalization, and the “FindVariableFeatures” function was used to identify highly variable genes on the basis of their average expression values (greater than 0.1 and less than 3) and dispersion (greater than 0.5). The R package “Harmony” was used to perform batch correction between the samples to avoid batch effects interfering with downstream analysis. The data were then scale transformed and downscaled via principal component analysis (PCA), and the top 50 principal components were selected for downstream analysis and visualized via the “RunTSNE” function.

Malignant cells in which at least two model genes were detected were selected for subsequent analysis. After standardization, normalization, identification of highly variable genes, removal of batch effects and PCA, the first 50 principal components were selected at a resolution of 0.1. Three subtypes of tumor cells were subsequently identified by clustering and grouping again. The marker genes of each subtype of tumor cells were identified via the “FindAllMarkers” function (avg_log2fc > 0.25, p_val_adj < 0.05). The CellScore was calculated based on the genes included in the prognostic model via the “AddModuleScore” function of the “Seurat” package. The malignant cells were divided into high and low groups based on the median cell score.

The R package “monocle2” was used to conduct the trajectory analysis of the tumor cells. Different states reflect the internal transformation of tumor cells. The R package “CellChat” was used to analyze the communication between tumor cells and other cells.

The expression levels of the genes included in the lung T_RM cell-related prognostic model were validated at the protein level using immunohistochemical staining images from the Human Protein Atlas database (https://www.proteinatlas.org/). The staining intensity levels of each gene in normal lung tissues and LUAD tissues were observed and compared.

Lung adenocarcinoma samples were collected from the pathology department of Shandong Provincial Hospital from 2017 to 2021. Written informed consent was obtained from all participants. Tumor tissues were obtained from excised biopsies, fixed in formalin and embedded in paraffin (FFPE) for histological evaluation. After paraffin wax removal and rehydration, the sections were placed in citrate antigen retrieval solution and boiled for 15 minutes for antigen retrieval. An endogenous peroxidase blocker was then added to block the endogenous peroxidase activity in the sections. After incubation at room temperature for 30 min, 50 µL of goat serum working solution was added to each sample, which was subsequently incubated at 37°C for 20 min to block nonspecific staining. The sections were subsequently incubated with a primary antibody (rabbit anti-thymidylate synthase antibody, 1:100, ab108995, Abcam) for 1 h at 37°C. After 3 × 5-minute washes with PBS, the sections were incubated with a biotinylated secondary antibody at room temperature for 30 min, followed by subsequent washes (3 × 5 min in PBS). The sections were subsequently dried with absorbent paper and incubated with 50 µL horseradish peroxidase-labeled streptavidin for 20 min at 37°C. The sections were then rinsed with PBS for 3 × 5 min each. After immunostaining, the sections were visualized using an MBMbio Intelligence 400 scanner according to the manufacturer’s protocol. The slides were independently examined by two experienced pathologists according to the WHO criteria. The expression levels of each gene were characterized using a scoring system. The staining intensity was graded into four levels: 0, no positive staining (negative), 1 point for light yellow (weakly positive), 2 points for brownish yellow (positive); and 3, dark brown (strongly positive). The percentage of positive cells was also classified into four levels: 1 point was given when it was ≤25%, 2 points when it ranged from 26% to 50%, 3 points when it was between 51% and 75%, and 4 points when it was >75%. The final scoring results were obtained by multiplying the scores of the above two items. Based on the results, the samples were divided into four grades: negative expression (0 points), low expression (1–4 points), moderate expression (5–8 points) and high expression (9–12 points).

The human cell line, H1395, was purchased from the National Laboratory Cell Resource Sharing Platform (Beijing, China) at the beginning of this study, with STR authentications. H1395 cells were cultured in RPMI 1640 medium supplemented with 10% fetal bovine serum (FBS) and 100 U/mL penicillin/streptomycin (Invitrogen, Carlsbad, CA, USA) at 37°C in a humidified incubator with 5% CO₂.

The siRNA oligo sequences (5’-3’) against TYMS mRNAs (si-TYMS-1#: sense, GGGAUUCUCCACCAGAGAATT; antisense, UUCUCUGGUGGAGAAUCCCTT; si-TYMS-2#: sense, CCAACUGCAAAGAGUGAUUTT; antisense, AAUCACUCUUUGCAGUUGGTT) were synthesized by GenePharma Co. (Shanghai, China). H1395 cells were transfected with the siRNAs using Omifection-R (OMIGET, China) siRNA transfection reagent according to the manufacturer’s instructions when the cells reached a confluence of 60–80% confluence. The successful knockdown of TYMS expression was confirmed by quantitative RT-PCR (qRT-PCR) and western blotting 48 h post-transfection. Scramble siRNAs (sense: 5’-UUCUCCGAACGUGUCACGUTT-3’; antisense: 5’-ACGUGACACGUUCGGAGAATT-3’) were used as negative controls.

The total RNA of the cell lines was isolated using the Total RNA Isolation kit (TRIcom Reagent) of GenStone Biotech and then reverse-transcribed into cDNA using TransScript First-Strand cDNA Synthesis SuperMix (TransGen Biotech, China) according to the manufacturer’s instructions. Next, qRT-PCR was performed using the FastStart Universal SYBR Green Master (ROX) (Roche, Germany) on an ABI-7500 Fast system (Applied Biosystems). ALU was used as the endogenous reference gene for the cultured cell lines. Each sample was analyzed quantitatively in six replicates. The relative expression levels of these genes were determined using the ΔΔCt method. The differences in target gene expression between different groups were analyzed using the Kruskal-Wallis test and plotted using GraphPad Prism 10.1.2. P < 0.05 was considered statistically significant (***indicates p < 0.001). The primer sequences are shown in Table 1 .

Total protein was extracted from cultured cells using RIPA buffer. Primary polyclonal antibodies against TYMS (15047-1-AP, ProteinTech) and β-Actin (66009-1-Ig, ProteinTech) were used at dilutions of 1:3,000 and 1:20,000, respectively. The signals were visualized using an enhanced chemiluminescence kit (Millipore) and an Alpha Imager system.

The long-term dynamic proliferation of the H1395 cells was observed using a long-term dynamic observation platform (IncuCyte, Essen, MI, USA). The cells were seeded into 96-well plates (3000 cells per well, six wells per group) and cultured for 120 h to generate proliferation curves. The cells were photographed every 24 h on the platform and analyzed using IncuCyte ZOOM software (Essen, Ann Arbor, MI, USA).

All analyses were performed using R software (version 4.4.2). For significance analysis between various values (such as expression levels, infiltration ratios, and various eigenvalues), the Wilcoxon rank-sum test was applied to compare the differences between two groups of samples, and the Kruskal-Wallis test was used to compare the differences between multiple groups of samples. For the plot presentation, ns indicates p > 0.05, * indicates p < 0.05, ** indicates p < 0.01, *** indicates p < 0.001, and **** indicates p < 0.0001. Survival curves for the prognostic analysis were generated using the Kaplan-Meier method, and the significance of the differences was determined using the log-rank test.

A flow chart of the study is shown in Figure 1 . To assess whether the expression of lung T_RM-cell-specific genes affects tumorigenesis and tumor progression in LUAD, differential expression analysis was performed in LUAD and adjacent normal tissues. First, 1002 downregulated genes and 741 upregulated genes in tumor tissues were screened ( Figures 2A, B ; Supplementary Table 1 ), including 130 lung T_RM cell-specific genes ( Figure 2C ; Supplementary Table 2 ). Protein-protein interaction network (PPI) analysis results revealed extensive interactions among the DE lung T_RM cell-specific genes, and the node connectivity of the RRM2, CDK1, CCNA2 and EXO1 genes was relatively high, which may indicate that these genes play a dominant role in the regulatory network ( Figure 2D ).

The results of functional enrichment analysis revealed that the DE lung T_RM cell-specific genes were significantly enriched in biological processes such as cell cycle regulation and chromosome segregation and were significantly associated with functions such as MHC class II molecule receptor activity, antigen binding, and immune receptor activity ( Figures 2E1–E4 ), suggesting that these genes are related to T_RM cells.

To investigate the clinical value of the DE lung T_RM cell-specific genes, a lung T_RM cell-related prognostic signature was constructed and validated. First, a univariate Cox regression analysis was performed. There were 62 genes associated with overall patient survival, and the top 20 genes with the greatest significance are shown in Figure 3A ; Supplementary Table 3 . The KM curves of the top six genes with the lowest p-values are presented in ( Figures 3B1–B6 ). Least absolute shrinkage and selection operator (LASSO) regression analysis was subsequently conducted to further investigate the clinical significance of these genes. The trajectory of each independent variable was obtained ( Figure 3C ), and as the lambda gradually increased, the number of independent variable coefficients gradually decreased to zero ( Figure 3C ). Ten-fold cross-validation was used to build the model, and the confidence intervals for each lambda value are shown in Figure 3D . Fourteen genes were identified when the model was optimized. Therefore, we selected the 14 genes for the subsequent analyses and constructed a risk score model based on their coefficients and expression levels of the 14 genes ( Figure 3E ). The formula for calculating the risk-score model is as follows:

Score = SLC16A3 * (0.128) + ARHGAP11A * (0.031) + PTTG1 * (0.020) +DTL * (0.025) + GPRIN1 * (0.044) + EXO1 * (0.018) + GAPDH * (0.128) + TYMS * (0.039) + DAPK2 * (-0.023) + CCL20 * (0.038) + HLA-DQA1 * (-0.076) +ADAM12 * (0.010) + ALOX5AP * (-0.040) + OASL * (0.014).

Using the 14-gene risk score model, the samples in the TCGA-LUAD training cohort were divided into high- and low-risk groups according to the median risk score. Overall survival analysis revealed that the OS of patients in the high-risk group was significantly lower than that of patients in the low-risk group in both the training cohort (TCGA-LUAD) ( Figure 4A ) and the two validation cohorts: GSE41271 ( Figure 4D ) and GSE42127 ( Figure 4G ). The ROC curve revealed that the AUCs of the patients at 1, 3, and 5 years were relatively high (0.688, 0.698, and 0.648, respectively) in the training cohort ( Figure 4B ). The AUCs of patients at 1, 3, and 5 years were 0.649, 0.638, and 0.646, respectively, in validation cohort GSE41271 ( Figure 4E ). The AUCs of the patients at 1, 3, and 5 years were 0.867, 0.662, and 0.672, respectively, in validation cohort GSE42127 ( Figure 4H ). To test whether the risk score model was an independent prognostic factor for LUAD patients, we performed univariate and multivariate Cox regression analyses via the “coxph()” function in the R package “survival”. In all the training and validation cohorts, the risk score was an independent prognostic factor among other clinical features, such as age, sex, and tumor stage ( Figures 4C, F, I ). These results demonstrated that the 14-gene prognostic signature based on the DE lung T_RM cell-specific genes had strong prognostic efficacy with high robustness and generalizability.

The associations between the T_RM-related prognostic signature and patients’ clinicopathological features were further analyzed. The results revealed that the proportions of patients aged <60 years, with advanced-stage disease, and with a history of smoking was significantly greater in the high-risk score group than in the low-risk score group (p < 0.05). The proportions of patients of different sexes, ALK rearrangements, EGFR mutations, and KRAS mutations did not significantly differ between the high- and low-risk score groups ( Supplementary Figure 1 ). Patients aged <60 years, with advanced-stage disease, male sex, and a history of smoking had significantly higher risk scores than the other groups of patients, and there was no significant difference in the risk scores for patients with EGFR mutations or KRAS mutations ( Supplementary Figure 2 ).

To further explore the functions of lung-specific T_RM cells in the TME of LUAD, gene set enrichment analysis (GSEA) and immune cell infiltration analysis were conducted in the high- and low-risk score groups of patients. The results revealed that signaling pathways, such as P53, B-cell receptor, and MAPK, were significantly activated in the low-risk score group ( Figure 5A ; Supplementary Table 4 ), and immune-related biological processes, such as T-cell activation, proliferation, and B-cell activation, were also significantly activated in the low-risk score group ( Figure 5B ; Supplementary Table 4 ). Further analysis of immune cell infiltration revealed that the infiltration of immune cells, such as activated B cells, activated CD8⁺T cells, central memory CD4⁺T cells, central memory CD8⁺T cells, and effector memory CD8⁺T cells, was significantly greater in the low-risk score group ( Figure 5C ). These results demonstrate that patients in the low-risk score group had stronger antitumor immunity and greater infiltration of T_RM cells, which may be the reason for their longer survival time. The ESTIMATE, immunity, and stroma scores were significantly greater in the low-risk score group, whereas the tumor purity was significantly greater in the high-risk score group ( Figures 5D1–D4 ). Next, the expression levels of the immune checkpoint genes were compared between the high- and low-risk score groups. The results revealed that the expression levels of several immune checkpoint genes, including CD276 and LAG3, were significantly different between the two groups ( Figure 5E ). These findings suggest the possibility of exploring novel targets for immunotherapy.

A total of 51935 cells, including 4827 B lymphocytes, 635 endothelial cells, 10998 epithelial cells, 1764 fibroblasts, 1735 MAST cells, 9098 myeloid cells, 22878 T/NK cells, and 27578 cells, were detected in the GSE131907 scRNA-seq cohort ( Supplementary Figure 3 ). The PCA results revealed that there was a significant batch effect between samples ( Supplementary Figures 4A, B ), and the batch effect between samples was removed via the R package “Harmony” ( Supplementary Figures 4C, D ). The distribution of different cell types was determined via UMAP analysis ( Supplementary Figure 4E ), and heterogeneity in the distribution of cells among the samples was detected ( Supplementary Figure 4F ).

A total of 3906 malignant tumor cells with at least two model genes detected were extracted for subsequent analyses. These malignant cells were renormalized and clustered, and three subtypes of malignant cells were identified ( Figures 6A, B1–B3 ). These cell subtypes were defined according to the genes that were highly expressed in the clusters, and these three subtypes were named IFI27 ⁺Mal, FBXO2 ⁺Mal and HMGB2 ⁺Mal ( Figures 6A, B1–B3 ). With the “FindAllMarkers” function, we identified the marker genes of each cell subtype ( Supplementary Table 5 ), and the top 5 marker genes of each cell subtype are shown in Figure 6C . FBXO2 ⁺Mal highly expressed genes such as CXCL14, TNNC2 and ASS1, which were significantly enriched in biological processes such as the regulation of the apoptosis signaling pathway and peptidase activity ( Supplementary Figure 5A ). HMGB2 ⁺Mal highly expressed genes such as STMN1, TUBA1B and UBE2C, which were significantly enriched in biological processes such as the regulation of cell adhesion, leukocyte migration and leukocyte chemotaxis ( Supplementary Figure 5B ). IFI27 ⁺Mal highly expressed genes such as SFTPA2, SFTPA1 and SCGB3A1, which were significantly enriched in biological processes such as the regulation of cell adhesion, leukocyte migration and leukocyte chemotaxis ( Supplementary Figure 5C ).

The CellScore of each malignant tumor cell line was calculated via the “AddModuleScore” function ( Figure 6D ), and the malignant tumor cells were divided into high- and low-CellScore groups ( Figure 6F ). Among the three cell subtypes, FBXO2 ⁺Mal and HMGB2 ⁺Mal had higher CellScores ( Figure 6E ), and the CellGroup was high ( Figure 6G ). GSEA of the cells in high- and low-CellScore groups revealed that immune-related biological processes, such as T-cell migration and the B-cell receptor signaling pathway, were also significantly activated in the low-CellScore group ( Supplementary Figure 6 ; Supplementary Table 6 ).

Trajectory analysis of the extracted malignant epithelial cells revealed three differentiation states ( Figure 7A ). In the trajectory from State1 to State2 cells, the IFI27 ⁺Mal subpopulation decreased significantly, whereas the HMGB2 ⁺Mal subpopulation increased significantly ( Figure 7D ). In the State1 to State3 cell trajectories, the proportion of the FBXO2 ⁺Mal subpopulation increased, but the HMGB2 ⁺Mal subpopulation also increased ( Figure 7D ). In the trajectory from State1 to State2, there was no significant increase in the CellScore. However, in the trajectory from State1 to State3, there was a significant increase in the CellScore ( Figures 7B, C ) and an increase in the proportion of high-cell groups ( Figures 7E, F ). This suggests that the malignancy of the tumor cells increased from low to high in the trajectory.

To further analyze the differences in physiological activity between high- and low-neoplastic populations, cell-to-cell communication was analyzed via the “CellChat” package. Extensive cellular communication was observed between the cell populations ( Supplementary Figure 7A ). High neoplastic cells were more likely to be outgoing signaling-dominant senders than low neoplastic cells ( Supplementary Figure 7B ), and different cell populations were found to have outgoing signaling patterns in different biological pathways ( Supplementary Figures 7C, D ). Compared with low neoplastic patients, high neoplastic patients exhibited specific cellular communication in the CSF and KIT signaling pathways ( Supplementary Figures 7E, F ).

The IC50 values of the drugs in the training cohort were predicted using the R package “oncoPredict” via the use of the drug information from the GDSC database combined with the expression profiles of the training set. Spearman correlation analysis was performed between the prognostic signature and the log2(IC50) value for each drug ( Supplementary Table 7 ). Patients in the high-risk group had a poorer prognosis; therefore, the top six drugs with the most significant negative correlations were selected according to the absolute values of the correlation coefficients. The six drugs used were AZD6738_ 1917, BI.2536_1086, docetaxel_1007, docetaxel_1819, MK.1775_179, and paclitaxel_1080 (p<0.05). The log2(IC50) values of these six drugs were lower in the high-risk score group than those in the low-risk score group and had greater sensitivity ( Figures 8A1–A6 ). The Spearman correlation coefficients between the risk scores and drug log2 (IC50) values were also calculated, and the top 50 drugs were selected for display, which revealed that there was a correlation between gene expression levels and most of the drug log2 (IC50) values ( Figure 8B ).

To explore the predictive efficacy of the lung T_RM-related prognostic model for patients receiving anti-PD1/PD-L1 immunotherapy, two immunotherapeutic cohorts, GSE126044 and GSE135222, were used for the prognostic analysis. The results revealed that patients in the low-risk score group had a superior overall survival status compared with patients in the high-risk score group in both cohorts ( Figures 8C, E ). Patients with low-risk scores had a greater response rate ( Figure 8D ) and a greater progression-free rate to anti-PD1/PD-L1 immunotherapy ( Figure 8F ).

To validate whether the protein expression levels of the genes involved in the lung T_RM cell-related prognostic model were consistent with the RNA expression levels, immunohistochemical staining images were obtained from the Human Protein Atlas database (https://www.proteinatlas.org/). The results of immunohistochemical staining for SLC16A3, ARHGAP11A, PTTG1, GPRIN1 and TYMS were greater in LUAD tissues than in normal lung tissues, which was consistent with the RNA expression levels of these genes ( Figures 9A–F ). Immunohistochemical staining for HLA-DQA1, ALOX5AP and OASL was lower in LUAD tissues than in normal lung tissues, which was consistent with the RNA expression levels of these genes ( Figures 9A, G–J ). The LUAD proteome expression data and corresponding clinical information were obtained from the supplementary data of the study by Xu et al. (22). Survival analysis based on the proteome data revealed that SLC16A3, TYMS, ALOX5AP and OASL were risk factors for patient prognosis, whereas HLA-DQA1 was a protective factor ( Figures 10A–E ). These findings validated the functions of key genes in the T_RM cell-related prognostic signature.

Previous results revealed that TYMS was a hub gene in the T_RM cell-related prognostic signature in the PPI analysis and the construction of the risk score model ( Figures 2D , 3E ; Supplementary Table 3 ). However, the function of TYMS in LUAD has rarely been investigated. To further investigate the clinical and biological roles of TYMS in LUAD, we performed immunohistochemical staining experiments using LUAD samples and assessed the proliferation of LUAD cell lines. Images of LUAD samples with negative, low, moderate, and high TYMS expression, magnified 40 × 10 times under a light microscope, are presented in ( Figures 11A1–A4 ). The clinical samples included in this study were from a total of 30 patients with LUAD, 10 of whom were positive for TYMS and 20 of whom were negative. Survival analysis revealed that patients with negative TYMS staining in their tumor tissues had longer overall survival times than patients with positive TYMS staining (p=0.038, Figure 11B ). Patients with M1 tumors had a higher TYMS-positive staining rate than those with M0 tumors ( Figure 11C ), and patients with clinical stage IV tumors had a higher TYMS-positive staining rate than those with stages II and III tumors ( Figure 11D ). These findings indicated that TYMS is a risk factor for LUAD patients, which is consistent with our previous findings ( Figure 3E ; Supplementary Table 3 ).

In the cell proliferation experiment, first, TYMS was successfully knocked down through the siRNA oligos, as evaluated by qRT-PCR ( Figure 11E ) and western blotting ( Figures 11F, G ). The results of long-term dynamic observation experiments using the IncuCyte platform revealed that the proliferative capacity of H1395 cells was significantly impaired in the TYMS-knockdown groups compared with the control groups ( Figures 11H, I ), suggesting that TYMS may promote LUAD cell proliferation. These findings validate the role of TYMS in enhancing the growth of H1395 cells, suggesting its potential in promoting LUAD progression.

To further verify the link between tissue-resident memory T cells infiltration and the expression of the T_RM cell-related prognostic predictors, Spearman’s correlation analysis was performed between the expression levels of the marker genes of T_RM cells: CD8, CD69, CD103, and the 14 genes included in the T_RM cell-related prognostic signature. The analytical results were consistent in the three datasets: TCGA ( Figure 12A ), GSE41271 ( Figure 12B ), and GSE42127 ( Figure 12C ). Interestingly, the expression levels of genes that were significant risk factors for LUAD patient prognosis, such as GAPDH, GPRIN1 and EXO1, were negatively correlated with the expression levels of T_RM cell marker genes ( Figures 3E , 12A–C ). In contrast, the genes that are protective factors for the prognosis of LUAD, such as HLA-DQA1 and ALOX5AP, were positively correlated with the three T_RM cell markers ( Figures 3E , 12A–C ). Moreover, the risk score of the LUAD patients was negatively correlated with the important T_RM cell marker gene CD69 in all three cohorts ( Figures 12A–C ). These results confirmed again that the prognostic predictors involved in the T_RM cell-related signature could affect patient prognosis by regulating the infiltration of T_RM cells.