Comprehensive profiling of T-cell exhaustion signatures and establishment of a prognostic model in lung adenocarcinoma through integrated RNA-sequencing analysis

Abstract

Background

T-cell exhaustion (TEX) in the tumor microenvironment causes immunotherapy resistance and poor prognosis.

Objective

We used bioinformatics to identify crucial TEX genes associated with the molecular classification and risk stratification of lung adenocarcinoma (LUAD).

Methods

Bulk RNA sequencing data of patients with LUAD were acquired from open sources. LUAD samples exhibited abnormal TEX gene expression, compared with normal samples. TEX gene-based prognostic signature was established and validated in both TCGA and GSE50081 datasets. Immune correlation and risk group-related functional analyses were also performed.

Results

Eight optimized TEX genes were identified using the LASSO algorithm: ERG, BTK, IKZF3, DCC, EML4, MET, LATS2, and LOX. Several crucial Kyoto encyclopedia of genes and genomes (KEGG) pathways were identified, such as T-cell receptor signaling, toll-like receptor signaling, leukocytes trans-endothelial migration, Fcγ R-mediated phagocytosis, and GnRH signaling. Eight TEX gene-based risk score models were established and validated. Patients with high-risk scores had worse prognosis (P < 0.001). A nomogram model comprising three independent clinical factors showed good predictive efficacy for survival rate in patients with LUAD. Correlation analysis revealed that the TEX signature significantly correlated with immune cell infiltration, tumor purity, stromal cells, estimate, and immunophenotype score.

Conclusion

TEX-derived risk score is a promising and effective prognostic factor that is closely correlated with the immune microenvironment and estimated score. TEX signature may be a useful clinical diagnostic tool for evaluating pre-immune efficacy in patients with LUAD.

Keywords

lung adenocarcinoma T-cell exhaustion prognosis estimate score

1 Introduction

Lung cancer remains a leading cause of cancer-related deaths worldwide, with an estimated 1.8 million deaths according to GLOBOCAN 2020 statistical data.¹ Lung adenocarcinoma (LUAD), the major pathological subtype of non-small cell lung cancer (NSCLC), is characterized by rapid progression and concomitant distal metastasis, which contribute to its suboptimal treatment response and diminished survival rates.² Recently, rapid progress in immunotherapy, particularly immune checkpoint inhibitors (ICIs), such as pembrolizumab and nivolumab, either alone or in combination with chemotherapy or targeted therapy, has markedly improved the prognosis of patients with NSCLC.^3,4 Notably, a significant challenge is that only a minority of patients with advanced disease benefit from ICIs, while the majority of initially responsive patients may later develop primary or secondary resistance.⁵ Therefore, there is a pressing need to advance our understanding of the key mechanisms governing immunotherapeutic efficacy and resistance, with the ultimate goal of enhancing patient response.

The tumor microenvironment (TME) significantly influences both cancer development and antitumor processes, potentially contributing to immunotherapy resistance. TME is composed of a diverse array of immune cells, stromal cells, and extracellular matrix molecules. Tumor-infiltrating lymphocytes (TILs) in the TME, primarily CD8+ T cells, have emerged as predictive biomarkers of immune efficacy. However, tumor-specific CD8+ T cells that are exposed to chronic antigen stimulation often enter a dysfunctional state known as “T-cell exhaustion” (TEX).⁶ Sustained stimulation of T-cell receptor (TCR) signaling and downstream Ca2+ signaling activates transcription factors of the NFAT, TOX, and NR4A families.^7,8 This signaling cascade leads to the upregulation of multiple inhibitory immune checkpoint proteins, such as PD-1, CTLA4, CD39, and LAG3.⁹ Additionally, exhausted T cells are characterized by progressive erosion of effector capabilities, memory recall deficits, metabolic dysregulation, impaired homeostatic self-renewal, and transcriptional and epigenetic reprogramming.^10–12 TEX is a reliable immunophenotype associated with decreased efficacy of ICIs and prognosis.^13,14 Investigating TEX within patient tumor tissues provides mechanistic insights for predicting immunotherapy responses and resistance, laying the foundation for enhancing the therapeutic effectiveness of combinatory approaches.

In this study, we screened for prognostic TEX molecules in LUAD and normal samples obtained from a public database. To create a predictive model for patient survival, we performed LASSO analysis and established a TEX risk score signature. Therefore, we developed a nomogram for predicting survival. Additionally, we explored the relationships between TEX risk scores and various factors, including immune infiltrating cells, gene set variation analysis (GSVA) pathway enrichment, and immune phenotypic variables such as tumor purity, stromal score, and estimated score.

2 Methods

2.1 Data source and preprocessing

Gene expression data for LUAD, specifically normalized to log₂ (FPKM + 1) expression values, were obtained from The TCGA (https://gdc-portal.nci.nih.gov/). Overall, 585 samples were used in this study. Gene expression profiling was performed using an Illumina HiSeq 2000 RNA Sequencing platform. Concurrently, a clinical survival set for the samples was also collected, allowing us to retain LUAD samples with survival and prognosis data as well as normal control samples. Our analysis included 501 LUAD and 58 normal samples from TCGA, which served as the training dataset.

For additional validation, we retrieved LUAD datasets from the GEO database (http://www.ncbi.nlm.nih.gov/geo/), as outlined in. Specifically, we obtained the GSE50081 dataset,¹⁵ which consists of 181 tumor tissue samples as well as prognostic data. Among these, 127 samples were LUAD tumor specimens that were used as the validation dataset. Gene expression profiling was performed using GPL570 Affymetrix Human Genome U133 Plus 2.0.

2.2 TEX score assessment

We initially extracted the published TEX factors from references¹⁶ and assessed the TEX scores for the samples based on whole-genome expression data from TCGA, utilizing R3.6.1 GSVA (http://www.bioconductor.org/packages/release/bioc/html/GSVA.html) Version 1.36.3.¹⁷ We then analyzed the disparity in TEX score distributions between LUAD and normal samples using the Kruskal–Wallis test.

Additionally, we integrated clinical prognostic information from LUAD tumor samples and employed the Xtile tool (https://en.freedownloadmanager.org/Windows-PC/X-tile-FREE.html) to identify the optimal TEX score threshold and stratify the samples into high and low TEX score groups. Subsequently, we assessed the correlation between the TEX score groups and prognosis using Kaplan–Meier analysis in the R3.6.1 language survival package (version 2.41-1, http://bioconductor.org/packages/survivalr/).

2.3 Screening differential expressed TEX genes between LUADs and normal samples

In TCGA profiles, we extracted the expression levels of TEX genes from the corresponding samples. Differential analysis was performed on the LUAD vs. CTRL comparison group using Limma (version 3.34.7, https://bioconductor.org/packages/release/bioc/html/limma.html).¹⁸ Significantly differentially expressed TEX genes were selected using a threshold of FDR value of < 0.05.

Furthermore, the COR function in R3.6.1 (https://stat.ethz.ch/R-manual/R-devel/library/stats/html/cor.test.html) was employed to calculate the correlation between the expression levels of differentially expressed genes. Genes showing significant correlations were visualized in a network using Cytoscape Version 3.9.0 (http://www.cytoscape.org/).

Finally, Gene Ontology (GO) biological processes and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were conducted for these TEX genes using DAVID¹⁹ (https://david.ncifcrf.gov/) with a significance threshold of p less than 0.05.

2.4 TEX gene-based prognostic model

For the screened TEX genes, a univariate Cox regression analysis was conducted using the survival package. Genes with P-values less than 0.05 were selected based on their significant associations with survival outcomes. Subsequently, the LASSO algorithm was applied to the target gene set using the lars package²⁰ (Version 1.2, https://cran.r-project.org/web/packages/lars/index.html) to perform regression analysis, thereby selecting an optimized combination of relevant genes.

We constructed a Risk Score (RS) model according to the LASSO coefficients and gene expression levels. The RS is calculated as follows:

\begin{aligned} Risk score = \frac{e^{sum (each gene's expression levels * corresponding coefficient)}}{e^{sum (each gene's mean expression levels * corresponding coefficient)}} \end{aligned}

Moreover, the RS was calculated for these samples in the training and validation sets. The entire sample was divided into high- and low-risk groups based on median scores. The correlation between the RS value and actual survival time was analyzed using a survival package.

2.5 Nomogram model for prognosis prediction

In the training dataset of TCGA tumor samples, we integrated clinical information and performed Cox regression analyses. Clinical factors that independently influenced the prognosis were screened based on a significance threshold.

Independent clinical characteristics, RS status, and risk factors were integrated to generate nomograms for 1-, 3-, and 5-year survival rate predictions using the rms package version 5.1-2 (https://cran.r-project.org/web/packages/rms/index.html). Subsequently, we calculated the C-index using R3.6.1, survcomp (version 1.34.0, http://www.bioconductor.org/packages/release/bioc/html/survcomp.html).²¹ The C-index is a measure used to rank the survival times of all individual pairs based on Harrell's C-statistics. Furthermore, we conducted a decision curve analysis for each factor and the nomogram combined model using the R3.6.1 rmda package Version 1.6²² (https://cran.r-project.org/web/packages/rmda/index.html) to assess the prognostic benefit of each factor.

2.6 Correlation analysis of risk group and immune

We assessed the proportion of immune cells in samples from the TCGA-LUAD dataset using the CIBERSORT tool (https://cibersort.stanford.edu/index.php).²³ Then, we utilized the “estimate” package (http://127.0.0.1:29606/library/estimate/html/estimateScore.html) to calculate ESTIMATE scores, immune scores, stromal scores, and tumor purity for LUAD samples in the TCGA dataset. the Kruskal–Wallis test was used to compare immune cell proportions and ESTIMATE scores among the different risk groups. Finally, we analyzed the correlation among immune cell infiltration, ESTIMATE scores, and TEX gene expression.

2.7 KEGG pathway analysis between differential risk groups

Using gene expression data from TCGA, we employed Gene Set Enrichment Analysis (GSEA, http://software.broadinstitute.org/gsea/index.jsp)²⁴ to identify the KEGG signaling pathways related to the risk groups. Based on the three key parameters (enrichment score, ES; normalized enrichment score, NES; nominal P-value), we ultimately selected a P-value less than 0.05 as the threshold to screen the significantly enriched KEGG pathways.

We downloaded all KEGG pathways and their corresponding genes from MSigDB²⁵ and conducted a quantitative analysis of these pathways using the GSVA package Version 1.36.3. Moreover, we used Kruskal–Wallis tests to assess the distribution differences of the pathways. Then, using the Cor function in R3.6.1, we calculated the correlation between the quantified pathways and the TEX gene expression levels screened from the previous model.

2.8 Real-time quantitative PCR (rt-qPCR)

Expression of the identified TEX genes was further validated in human cells using RT-qPCR. Human bronchial epithelial (BEAS-2B) and lung cancer cell lines (A549 and NCI-H1299) were purchased from the Cell Bank of the Chinese Academy of Sciences (Shanghai, China). BEAS-2B and A549 cells were maintained in Dulbecco's modified Eagle's medium (DMEM) containing 10% fetal bovine serum (FBS) and 1% penicillin/streptomycin, while NCI-H1299 cells were cultured in RPMI1640 medium supplemented with 10% FBS and 1% penicillin/streptomycin. All cells were incubated in an incubator with 5% CO₂ at 37 °C. After culturing to 80–90% confluence, total RNA was isolated from different cells using TRIzol reagent. The isolated total RNA was reverse-transcribed into cDNA using the PrimeScript 1st Strand cDNA Synthesis Kit (Takara, Japan). The RT-qPCR reaction was initiated at 95 °C for 10 min; followed by a total of 40 cycles at 95 °C for 15 s and 60 °C for 60 s. Glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as the reference gene, and the sequences of all primers are shown in Table 1. The relative mRNA expression levels of related genes were calculated by the 2⁻ $^{Δ Δ}$ ^CT method.

Table 1.
Primer sequences for RT-qPCR assay.

Primer Sequences (5′-3′)

GAPDH-hF TGACAACTTTGGTATCGTGGAAGG

GAPDH-hR AGGCAGGGATGATGTTCTGGAGAG

ERG-hF CGTGCCAGCAGATCCTACG

ERG-hR GGTGAGCCTCTGGAAGTCG

DCC-hF ACCCAAGCTGGCTTTTGTACT

DCC-hR TGTGACGGCATCAGAAGGTTC

MET-hF AGCAATGGGGAGTGTAAAGAGG

MET-hR CCCAGTCTTGTACTCAGCAAC

LATS2-hF ACTTTTCCTGCCACGACTTATTC

LATS2-hR GATGGCTGTTTTAACCCCTCA

LOX-hF CGGCGGAGGAAAACTGTCT

LOX-hR TCGGCTGGGTAAGAAATCTGA

Primer	Sequences (5′-3′)
GAPDH-hF	TGACAACTTTGGTATCGTGGAAGG
GAPDH-hR	AGGCAGGGATGATGTTCTGGAGAG
ERG-hF	CGTGCCAGCAGATCCTACG
ERG-hR	GGTGAGCCTCTGGAAGTCG
DCC-hF	ACCCAAGCTGGCTTTTGTACT
DCC-hR	TGTGACGGCATCAGAAGGTTC
MET-hF	AGCAATGGGGAGTGTAAAGAGG
MET-hR	CCCAGTCTTGTACTCAGCAAC
LATS2-hF	ACTTTTCCTGCCACGACTTATTC
LATS2-hR	GATGGCTGTTTTAACCCCTCA
LOX-hF	CGGCGGAGGAAAACTGTCT
LOX-hR	TCGGCTGGGTAAGAAATCTGA

3 Results

3.1 Screening of crucial TEX genes in LUAD versus normal samples

First, based on a comprehensive analysis of whole-genome expression profiles, we assessed the TEX scores of LUAD and normal samples. The results showed that the LUAD samples exhibited lower TEX scores (P < 0.001; Figure 1(a)). Considering the clinical survival data, we used the Xtile tool to define the optimal cut-off value. Finally, we selected 0.59 a threshold to classify the LUAD samples into two distinct groups: the low- and high-TEX score groups, which consisted of 322 samples and 179 samples, respectively. Kaplan–Meier curves showed that patients with high TEX scores exhibited better survival (P = 0.029, HR: 0.707 [0.517–0.967], Figure 1(b)).

Figure 1.

Analysis of crucial T-cell exhaustion (TEX)-related genes in samples with lung adenocarcinoma (LUAD), compared with control samples. (a) TEX score distribution in LUAD and control samples. (b) Kaplan–Meier curve displays the prognosis in two groups. (c) Aberrant expressed TEX genes in LUAD samples. (d) Network displays the expressed correlation of TEX gene. The red lines indicate positive expressed correlation, while green lines represent negative expressed correlation. (e) Kyoto encyclopedia of genes and genomes (KEGG) pathways related to crucial TEX genes.

Moreover, by retrieving previous literature, we identified 40 TEX-related genes and extracted their expression values from tumor samples. Using the limma package, we identified 31 differentially expressed TEX genes between LUADs and control samples (Figure 1(c)). Using the defined screening criteria (P < 0.05, absolute correlation values >0.3), we obtained a total of 183 gene pairs (Figure 1(d)).

Subsequently, we conducted a functional enrichment analysis and found 20 biological processes and eight KEGG signaling pathways associated with these TEX genes, including signal transduction, immune response, protein phosphorylation, pathways in cancer, Barr virus infection, Th17 cell differentiation, and transcriptional misregulation in cancer (Figure 1(e)).

3.2 Establishment and validation of TEX gene-based prognostic signature

After integrating the clinical prognostic data of these 31 TEX genes, we identified 16 genes that were significantly associated with survival outcomes (Figure 2(a)). Using the LASSO algorithm, we identified eight optimized TEX-related DEGs: ERG, BTK, IKZF3, DCC, EML4, MET, LATS2, and LOX (Figure 2(b)).

Figure 2.

Construction and validation of the T-cell exhaustion (TEX) risk model. (a) Forest map of univariate cox regression analysis. (b) LASSO parameter diagram. (c) Kaplan–Meier curves showed the correlation of eight TEX-related gene expression levels and survival outcomes. d-e. Kaplan–Meier curves and receiver operating character (ROC) curves display the accuracy of risk score for predicting survival in the cancer genome atlas (TCGA) (d) and GSE50081 cohort (e) Left: Kaplan–Meier curves associated with eight optimal TEX-associated gene prognostic models. Middle: Risk score distribution and actual survival time. Right: ROC curves of eight optimal TEX-related gene prognostic models. Numbers represent the specificity and sensitivity of the ROC curves.

Samples were classified into high- and low-expression groups. The association between TEX expression and patient survival was also assessed. Kaplan–Meier curves showed that patients with high expression levels of ERG, DCC, BTK, IKZF3, and EML4 had a better prognosis (Figure 2(c), P < 0.05). Moreover, patients with low expression of MET, LATS2, and LOX showed longer overall survival times (P < 0.05).

The prognostic prediction values of the TEX signature were verified using the TCGA training set and the GSE50081 validation dataset (Figure 2(d) and (e)). After calculating the risk score, the patients were divided into high- and low-risk groups. Furthermore, higher risk scores were clinically related to poor prognosis, and the results were consistent in both the TCGA and GSE50081 datasets (P < 0.001). The areas under the curve (AUCs) for the risk scores of the two datasets were 0.811 and 0.746, respectively.

Table 2.

Identification of independent clinical prognostic factors according to cox regression analysis.

	Uni-variable cox		Multi-variable cox
Clinical characteristics	HR (95% CI)	P value	HR (95% CI)	P value
Age (years, mean ± sd)	1.009 [0.994–1.024]	2.64E-01	–	–
Gender (male/female)	1.060 [0.792–1.418]	6.95E-01	–	–
Pathologic M (M0/M1)	2.111 [1.232–3.616]	5.36E-03	0.174 [0.0368–1.523]	2.74E-01
Pathologic N (N0/N1/N2/N3)	1.710 [1.443–2.027]	2.00E-10	0.812 [0.484–1.361]	4.29E-01
Pathologic T (T1/T2/T3/T4)	1.550 [1.289–1.863]	3.10E-06	1.105 [0.812–1.503]	5.24E-01
Pathologic stage (I/II/III/IV)	1.679 [1.463–1.928]	2.94E-14	2.223 [1.248–3.957]	6.66E-03
Tumor recurrence (yes/no)	2.392 [1.700–3.367]	2.44E-07	2.735 [1.794–4.169]	2.94E-06
Smoking history (current/reform/never)	0.765 [0.542–1.081]	1.29E-01	–	–
RS model status (high/low)	2.044 [1.508–2.771]	2.26E-06	2.049 [1.305–3.220]	1.85E-03

3.3 Construction of nomogram model for predicting survival rate of LUAD

To better utilize the TEX signature, we selected TCGA-LUAD for Cox regression analyses. Finally, we identified three independent prognostic clinical characteristics: pathological stage, tumor recurrence, and risk score (Table 2). The overall survival of patients with recurrence was significantly shorter than that of patients without recurrence (P = 2.438 × 10⁻⁷) and the overall survival worsened with an increase in the pathologic stage (Figure 3(a)). The risk score in patients with recurrence was significantly higher than that in patients without recurrence (P = 0.0299) and gradually increased with an increase in the pathologic stage (Figure 3(b)).

Figure 3.

Identification of prognostic-related clinical factors. (a) Kaplan–Meier curves for different recurrence status group and pathological stage. (b) Risk score distribution between differential disease status and differential pathological stage.

Subsequently, a nomogram integrating the three clinical factors was generated to predict the 1-, 3-, and 5-year survival rates of patients with LUAD (Figure 4(a)). The survival probability of the patients was evaluated using the point axis generated by adding each factor. Receiver operating characteristic (ROC) curves and calibration dots were used to assess efficacy. The AUCs were 0.756, 0.821, and 0.885. Calibration curves showed that the predicted survival and actual OS were consistent in the training cohort (Figure 4(a)). The decision curve revealed that this risk model had better clinical benefit than all other clinical factors (Figure 4(b)).

Figure 4.

Nomogram model for survival rate prediction in lung adenocarcinoma (LUAD). (a) Construction of nomogram model (up diagram) to predicted 1-, 3-, and 5-year survival probability and accuracy validation for actual survival outcomes (down diagram). Horizontal indicates predictive survival, while the vertical refers to actual prognosis. (b) Decision analysis curve of the LUAD samples.

3.4 Correlation analysis of risk signature and immunity

We employed CIBERSORT to evaluate immune cell proportions in TCGA samples and subsequently calculated the ESTIMATE scores. Thirteen immune cells exhibited significantly different distributions, including memory B cells, CD8+ T cells, plasma B cells, CD4+ memory resting T cells, resting NK cells, regulatory T cells (Tregs), M0-M2 macrophages, activated myeloid dendritic cells, activated or resting mast cells, and neutrophils (Figure 5(a)). Furthermore, the ESTIMATE score differed significantly between the two risk groups (Figure 5(b)). Subsequently, the correlations between the 13 immune cell proportions, four types of estimate scores, and eight TEX gene expression were analyzed (Figure 5(c)). We found that four types of estimate scores were significantly correlated with TEX gene expression, in addition to EML4 (P < 0.001). Additionally, the expression of BTK and IKZF3 was related to four types of estimated scores and 13 immune cell infiltrations (P < 0.05).

Figure 5.

Correlation analysis between risk score and immunity. (a) Immune cell infiltration assessment between the two risk score groups. (b) Distribution of four estimate score variate in the different risk groups. (c) Correlation between immune cells proportion, estimate score, and TEX genes expression.

3.5 Functional enrichment analysis

Based on the whole-genome gene expression levels in TCGA-LUAD, we used GSEA to identify KEGG pathways associated with the risk groups, and obtained 16 crucial pathways, such as B cell receptor-signaling, calcium signaling pathway, cell adhesion molecules, chemokine signaling, cytokine-cytokine receptor interaction, Fc epsilon RI signaling pathway, Fc-gamma R-mediated phagocytosis, GnRH signaling pathway, hematopoietic cell lineage, JAK-STAT signaling pathway, leukocyte trans-endothelial migration, natural killer cell-mediated cytotoxicity, neuroactive ligand-receptor interaction, TCR signaling pathway, toll-like receptor signal pathway, and vascular smooth muscle contraction. The details of several crucial pathways are shown in Figure 6.

Figure 6.

Functional enrichment analysis for T-cell exhaustion (TEX)-related genes. (a) Crucial signaling pathways associated with TEX signature. (b) Kyoto encyclopedia of genes and genomes (KEGG) score of signaling pathways display significant difference between two groups. (c) Network visualizes the correlation between sixteen KEGG signaling pathways and 8 TEX genes.

Furthermore, the KEGG scores of these pathways were assessed between the different risk groups, and 15 pathways exhibited significantly different quantified values between the two risk groups (Figure 6(b)). Next, we calculated the correlations between these pathways and the expression levels of eight model TEX-related genes. We retained relationships with P-values < 0.05 and absolute correlation values exceeding 0.3, resulting in a total of 71 connection pairs (Figure 6(c)).

3.6 Verification of the expression of the identified TEX-related genes by rt-qPCR

Four of the eight TEX-related genes (ERG, DCC, MET, and LOX) with the same expression trends as those shown in Figures 1(c) and 2(c) were chosen to determine their expression levels in lung cancer cells. Compared to control cells (BEAS-2B), the expression of ERG and DCC was significantly downregulated in lung cancer cells (A549 and NCI-H1299, P < 0.05); whereas the expression of MET and LOX was up-regulated in A549 and NCI-H1299 cells (P < 0.05, Figure 7). These results show that the expression trends of ERG, DCC, MET, and LOX were consistent with the bioinformatics analysis, implying the high relative reliability of our bioinformatics analysis.

Figure 7.

Expression of the identified T-cell exhaustion (TEX)-related genes (ERG, DCC, MET, and LOX) in the lung cancer cells. * P < 0.05, ** P < 0.01, compared with the BEAS-2B cells.

4 Discussion

LUAD, the most common type of lung cancer, accounts for ∼40% of all lung cancer cases and seriously threatens the health and life of patients.²⁶ Cancer is associated with TEX, which is a state of diminished function characterized by the progressive loss of T cell effector functions and self-renewal ability; as well as TEX is considered to be crucial pathway of resistance for cellular immunotherapy.¹³ However, the mechanism of action of TEX in LUAD progression remains unclear. Therefore, in this study, using TEX score assessment combined with bulk RNA sequencing data analysis, we established a TEX gene-based risk model that demonstrated remarkable and predictive utility for the prognosis of LUAD. Notably, patients with high-risk scores exhibited considerably worse prognoses, which was consistently observed in both the training and GSE50081 validation datasets. These abnormally expressed TEX genes were mainly related to 16 pathways, including TCR signaling, calcium signaling, chemokine signaling, and cell adhesion molecules. We identified three independent prognostic factors (pathological stage, tumor recurrence status, and risk score) to generate a nomogram model for survival rate prediction. Finally, the TEX score was correlated with the immune cell infiltration, estimated score, immune score, stromal score, and tumor purity.

In our study, eight optimized TEX-related DEGs were identified, including ERG, BTK, IKZF3, DCC, EML4, MET, LATS2, and LOX, which showed that ERG and DCC were significantly downregulated, whereas MET and LOX were upregulated in lung cancer cells. The endothelial ETS transcription factor ERG plays a key role in endothelial homeostasis, driving expression of lineage-specific genes, and repressing pro-inflammatory genes.^27,28 ERG actively promotes vascular development and angiogenesis via Wnt/β-catenin and angiopoietin-dependent Notch signaling.²⁹ It also acts as a repressor of pro-inflammatory genes, such as ICAM1 and IL89, and protects against endothelial-to-mesenchymal transition.³⁰ Additionally, ERG initiates a transcriptional regulatory network that controls the expression of crucial B cell lymphopoiesis genes, Ebf1 and Pax5.³¹ However, the aberrant expression of ERG has detrimental consequences and is associated with carcinogenesis. Recent studies have identified ERG as a transcriptional target of the oncogenic gene EVI1, with abnormal ERG activation being a requisite factor for both human and mouse EVI1 acute myeloid leukemia (AML). ERG suppression causes terminal differentiation of EVI1-related cancer cells, whereas ectopic ERG expression alleviates their dependence on EVI1.²⁸ In solid tumors, such as prostate cancer, TMPRSS2-ERG fusion is highly prevalent and often occurs in combination with the loss of PTEN, resulting in an aggressive invasive phenotype.³² Furthermore, TMPRSS2-ERG fusion, along with gain-of-function mutant p53, promotes prostate cancer through β-catenin activation and pyrimidine synthesis.³³ However, the precise involvement of ERG in LUAD remains unknown. In our study, we discovered that ERG was significantly associated with the estimated score, immune cell infiltration, and prognosis. These findings suggest that the ERG plays a crucial role in the progression of LUAD.

The netrin-1 guidance cue receptor, DCC, directs growing axons to appropriate targets and organizes fine neuronal connections throughout the life cycle by controlling target recognition, axon tree formation, and synaptic formation. Its role has been widely reported in psychopathology, including depressive and cognitive disorders.³⁴ A previous study by Li et al.³⁵ demonstrated that somatic mutations in DCC were associated with improved clinical outcomes in patients with melanoma treated with immune checkpoint blockade. Dysregulation of MET has been observed in NSCLC and its role as an oncogenic driver in NSCLC has been reported.³⁶ Nevertheless, our study confirmed the correlation between MET and TEX in LUAD. Recently, the up-regulation of LOX-1 in different tumors has demonstrated its involvement in cancer initiation, progression, and metastasis; as well as its roles in tumor spread and metastasis may be related to VEGF induction, HIF-1α activation, and MMP-9/MMP-2 expression.³⁷ Taken together, we speculated that DCC, MET and LOX closely related to TEX, may participate in the occurrence and progression of LUAD.

In addition, BTK functions as a major component of BCR signaling and is essential for B-cell proliferation.³⁸ The involvement of BTK in T-cell function has also been documented, where T-cell expression of BTK can lead to the activation of autoreactive T cells and worsen aplastic anemia.³⁹ Moreover, BTK holds promise as a predictive factor for lung adenocarcinoma and as a sign of TME alteration.⁴⁰ BTK inhibitors significantly reduce the proliferation of some EGFR-mutant lung cancers by blocking EGFR autophosphorylation.⁴¹ BTK is associated with resistance to EGFR inhibition in NSCLC by mediating stemness and EMT characteristics.⁴² IKZF3 encodes Aiolos, which can form heterodimers with IKZF/Ikaros to modulate lymphocyte differentiation.⁴³ IKZF1 and IKZF3 are vital transcription factors in multiple myeloma (MM). Notably, because of its substantial modulation of T-cell activity, increased IKZF3 expression in T-cells has been shown to be associated with improved clinical outcomes in patients with advanced MM receiving immunomodulatory treatment.⁴⁴ IKZF3 acted as a repressor of IL-2 expression. Treatment with the IKZF3 small-molecule inhibitor lenalidomide, either as a single dose or in combination with immune checkpoint blockade (ICB), partially rescued exhausted T cells within the TME.⁴⁵ In lung cancer, Aiolos expression predicts significantly shorter patient survival. Aiolos affects the protein structure of SHC1 and silences the anchorage reporter p66.⁴⁶ Aiolos overexpression promotes EMT and cancer stem cell-like characteristics in lung cancer.⁴⁷ Our study revealed that a unique IKZF3 transcriptome signature corresponds to favorable results in LUAD. IKZF3 expression was closely associated with the infiltration of CD4+ and CD8+ T cells and negatively correlated with the infiltration of M2 macrophages and neutrophils.

LATS2 is a major regulator of Hippo signaling, which controls organ size and cancer development. LATS2 serves as a classical tumor suppressor that stimulates mitochondrial fission and activates the JNK-Mff pathway.⁴⁸ Aberrant LATS2 expression in lung cancer corresponds to EGFR mutations, cisplatin response, and poor prognosis.^49,50 However, the functional roles of LATS1/2 in cancer immunity remain unclear. LATS1/2 deficiency can improve immunogenicity, resulting in tumor elimination by increasing the antitumor response. Mechanistically, tumor cells with LATS1/2 deletions trigger an interferon response through Toll-like receptors.⁵¹ EML4 is a microtubule-associated protein essential for ensuring proper chromosome congression during mitosis.⁵² Abnormal fusion events involving the EML4 and the kinase portions of ALK give rise to EML4-ALK fusion transcripts. This fusion gene was the primary translocation observed in NSCLC. The EML4-ALK fusion protein disrupts normal cellular signaling pathways, leading to excessive cell growth. Lung cancer harboring an oncogenic EML4-ALK fusion is responsive to ALK inhibitors.^53–55 Together with our findings, these reports suggest that BTK, IKZF3, EML4, and LATS2 are viable targets and prognostic markers in patients with LUAD.

However, this study has several limitations. First, the patient cohorts we examined were relatively limited. Despite the significant efficacy of the TEX signature in characterizing immunity and predicting patient outcomes, it is imperative to validate its effectiveness across a diverse range of patient cohorts in subsequent investigations. Second, our reliance on public databases for our analysis inevitably introduced the potential for deviation from real-world conditions into the prediction results. To enhance the usefulness of the model and improve the accuracy of immunotherapy predictions, it is crucial to gather additional data from patients with LUAD, as well as the need for associated clinical validation to continuously modify and improve the model according to actual values. Finally, a comprehensive understanding of the biological and molecular mechanisms governing these TEX genes requires dedicated biological experiments in vitro and in vivo.

5 Conclusion

In summary, we identified a TEX signature-based risk model for predicting survival in patients with LUAD. These TEX signature genes serve as pivotal tools for characterizing the immune environment of LUAD and show potential as predictive biomarkers and therapeutic targets. With the promise of further prospective validation, this signature has the potential to guide clinical decision-making for patients with LUAD, effectively addressing the growing demand for precision medicine.

Footnotes

Abbreviations

Acknowledgements

None.

Ethics approval and consent to participate

Ethical approval is not applicable because this is a bioinformatics analysis.

Consent for publication

Not applicable.

Author contributions

Conception and design of the research: YYZ; acquisition of data: JQC, PYJ; analysis and interpretation of data: LZL, HJY; statistical analysis: CXY, SZ; drafting the manuscript: YYZ; revision of manuscript for important intellectual content: SZ. All authors read and approved the final manuscript.

Funding

The authors received no financial support for the research, authorship, and/or publication of this article.

Declaration of conflicting interests

The authors declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.

Availability of data and materials

The datasets used and/or analyzed during the current study are available from the corresponding author upon reasonable request.

References

Sung

Ferlay

Siegel

, et al. Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J Clin 2021; 71: 209–249.

Mithoowani

Febbraro

. Non-small-cell lung cancer in 2022: A review for general practitioners in oncology. Current Oncology (Toronto, Ont) 2022; 29: 1828–1839.

Mok

TSK

Kudaba

, et al. Pembrolizumab versus chemotherapy for previously untreated, PD-L1-expressing, locally advanced or metastatic non-small-cell lung cancer (KEYNOTE-042): A randomised, open-label, controlled, phase 3 trial. Lancet (London, England) 2019; 393: 1819–1830.

Provencio

Nadal

González-Larriba

, et al. Perioperative nivolumab and chemotherapy in stage III non-small-cell lung cancer. N Engl J Med 2023; 389: 504–513.

Chaft

Shyr

Sepesi

, et al. Preoperative and postoperative systemic therapy for operable non-small-cell lung cancer. J Clin Oncol: Off J Am Soc Clin Oncol 2022; 40: 546–555.

McLane

Abdel-Hakeem

Wherry

. CD8T cell exhaustion during chronic viral infection and cancer. Annu Rev Immunol 2019; 37: 457–495.

Scott

Dündar

Zumbo

, et al. TOX Is a critical regulator of tumour-specific T cell differentiation. Nature 2019; 571: 270–274.

Liu

Wang

, et al. Genome-wide analysis identifies NR4A1 as a key mediator of T cell dysfunction. Nature 2019; 567: 525–529.

Blank

Haining

Held

, et al. Defining ‘T cell exhaustion’. Nat Rev Immunol 2019; 19: 665–674.

10.

Beltra

Manne

Abdel-Hakeem

, et al. Developmental relationships of four exhausted CD8(+) T cell subsets reveals underlying transcriptional and epigenetic landscape control mechanisms. Immunity 2020; 52: 825–841. e8.

11.

Franco

Jaccard

Romero

, et al. Metabolic and epigenetic regulation of T-cell exhaustion. Nat Metabol 2020; 2: 1001–1012.

12.

Khan

Giles

McDonald

, et al. TOX transcriptionally and epigenetically programs CD8(+) T cell exhaustion. Nature 2019; 571: 211–218.

13.

Chow

Perica

Klebanoff

, et al. Clinical implications of T cell exhaustion for cancer immunotherapy. Nat Rev Clin Oncol 2022; 19: 775–790.

14.

Lin

, et al. The identification of genes associated T-cell exhaustion and construction of prognostic signature to predict immunotherapy response in lung adenocarcinoma. Sci Rep 2023; 13: 13415.

15.

Zhu

Zheng

Liu

, et al. Significance of macrophage infiltration in the prognosis of lung adenocarcinoma patients evaluated by scRNA and bulkRNA analysis. Front Immunol 2022; 13: 1028440.

16.

Zhang

Chen

, et al. Pan-cancer landscape of T-cell exhaustion heterogeneity within the tumor microenvironment revealed a progressive roadmap of hierarchical dysfunction associated with prognosis and therapeutic efficacy. EBioMedicine 2022; 83: 104207.

17.

Choi

Jeon

Heo

, et al. Gene set enrichment analysis reveals that fucoidan induces type I IFN pathways in BMDC. Nutrients 2022; 14: 2242.

18.

Townsend

Daigle

BJ Jr

. GEOlimma: Differential expression analysis and feature selection using pre-existing microarray data. BMC Bioinf 2021; 22: 44.

19.

Sherman

Hao

Qiu

, et al. DAVID: A web server for functional enrichment analysis and functional annotation of gene lists (2021 update). Nucleic Acids Res 2022; 50: W216–Ww21.

20.

Begun

Kulinskaya

Ncube

. A double-cox model for non-proportional hazards survival analysis with frailty. Stat Med 2023; 42: 3114–3127.

21.

Shan

Chen

Jia

. Transcriptome analysis revealed a highly connected gene module associated with cirrhosis to hepatocellular carcinoma development. Front Genet 2019; 10: 305.

22.

Cao

, et al. A novel nomogram for predicting prognosis and tailoring local therapy decision for ductal carcinoma in situ after breast conserving surgery. J Clin Med 2022; 11: 5188.

23.

Guan

Jiao

Zhou

. Immune infiltration analysis with the CIBERSORT method in lung cancer. Dis Markers 2022; 2022: 3186427.

24.

Zito

Lualdi

Granata

, et al. Gene set enrichment analysis of interaction networks weighted by node centrality. Front Genet 2021; 12: 577623.

25.

Huang

, et al. Identification and validation of autophagy-related gene expression for predicting prognosis in patients with idiopathic pulmonary fibrosis. Front Immunol 2022; 13: 997138.

26.

Denisenko

Budkevich

Zhivotovsky

. Cell death-based treatment of lung adenocarcinoma. Cell Death Dis 2018; 9: 117.

27.

Caporarello

Lee

Pham

, et al. Dysfunctional ERG signaling drives pulmonary vascular aging and persistent fibrosis. Nat Commun 2022; 13: 4170.

28.

Schmoellerl

Barbosa

IAM

Minnich

, et al. EVI1 Drives leukemogenesis through aberrant ERG activation. Blood 2023; 141: 453–466.

29.

Birdsey

Shah

Dufton

, et al. The endothelial transcription factor ERG promotes vascular stability and growth through wnt/β-catenin signaling. Dev Cell 2015; 32: 82–96.

30.

Kharwadkar

Ulrich

Chu

, et al. ERG Functionally overlaps with other ets proteins in promoting TH9 cell expression of Il9 during allergic lung inflammation. J Immunol 2023; 210: 537–546.

31.

Coughlan

Hediyeh-Zadeh

, et al. An Erg-driven transcriptional program controls B cell lymphopoiesis. Nat Commun 2020; 11: 3013.

32.

Kobelyatskaya

Pudova

Snezhkina

, et al. Impact TMPRSS2-ERG molecular subtype on prostate cancer recurrence. Life (Basel, Switzerland) 2021; 11: 588.

33.

Ding

Blee

Zhang

J III

, et al. Gain-of-function mutant p53 together with ERG proto-oncogene drive prostate cancer by beta-catenin activation and pyrimidine synthesis. Nat Commun 2023; 14: 4671.

34.

Torres-Berro

Lopez

Bagot

, et al. DCC Confers susceptibility to depression-like behaviors in humans and mice and is regulated by miR-218. Biol Psychiatry 2017; 81: 306–315.

35.

Wang

Chen

, et al. Somatic mutations in DCC are associated with genomic instability and favourable outcomes in melanoma patients treated with immune checkpoint inhibitors. Br J Cancer 2022; 127: 1411–1423.

36.

Remon

Hendriks

LEL

Mountzios

, et al. MET alterations in NSCLC-current perspectives and future challenges. J Thorac Oncol 2023; 18: 419–435.

37.

Murdocca

De Masi

Pucci

, et al. LOX-1 and cancer: An indissoluble liaison. Cancer Gene Ther 2021; 28: 1088–1098.

38.

McDonald

Xanthopoulos

Kostareli

. The role of Bruton’s tyrosine kinase in the immune system and disease. Immunology 2021; 164: 722–736.

39.

Xia

Liu

Cao

, et al. T-cell expression of Bruton’s tyrosine kinase promotes autoreactive T-cell activation and exacerbates aplastic anemia. Cell Mol Immunol 2020; 17: 1042–1052.

40.

Wei

Qin

, et al. BTK has potential to be a prognostic factor for lung adenocarcinoma and an indicator for tumor microenvironment remodeling: A study based on TCGA data mining. Front Oncol 2020; 10: 424.

41.

Yeh

Chen

Satriyo

, et al. Bruton’s tyrosine kinase (BTK) mediates resistance to EGFR inhibition in non-small-cell lung carcinoma. Oncogenesis 2021; 10: 56.

42.

Giordano

Vaira

Cortinovis

, et al. p65BTK is a novel potential actionable target in KRAS-mutated/EGFR-wild type lung adenocarcinoma. J Exp Clin Cancer Res 2019; 38: 260.

43.

Cippitelli

Stabile

Kosta

, et al. Role of Aiolos and Ikaros in the Antitumor and Immunomodulatory Activity of IMiDs in Multiple Myeloma: Better to Lose Than to Find Them. Int J Mol Sci 2021; 22: 1103.

44.

Yang

Lin

, et al. Novel IKZF3 transcriptomic signature correlates with positive outcomes of skin cutaneous melanoma: A pan-cancer analysis. Front Genet 2022; 13: 1036402.

45.

Hay

Acklam

Turner

, et al. Abstract 5586: Demonstrating restoration of T cell function in exhausted T cells with IKZF3 small molecule inhibitor, lenalidomide. Cancer Res 2022; 82(12_Supplement): 5586.

46.

Wang

, et al. N160 of aiolos determines its DNA-binding activity. Anatomical Record (Hoboken, NJ : 2007) 2019; 302: 2014–2019.

47.

Hung

Kao

Huang

, et al. Overexpression of Aiolos promotes epithelial-mesenchymal transition and cancer stem cell-like properties in lung cancer cells. Sci Rep 2019; 9: 2991.

48.

Xie

Zhang

, et al. LATS2 Promotes apoptosis in non-small cell lung cancer A549 cells via triggering Mff-dependent mitochondrial fission and activating the JNK signaling pathway. Biomed Pharmacother 2019; 109: 679–689.

49.

Luo

Sit

Sihoe

, et al. Aberrant large tumor suppressor 2 (LATS2) gene expression correlates with EGFR mutation and survival in lung adenocarcinomas. Lung Cancer (Amsterdam, Netherlands) 2014; 85: 282–292.

50.

Luo

Kwok

Yang

, et al. Expression of large tumour suppressor (LATS) kinases modulates chemotherapy response in advanced non-small cell lung cancer. Transl Lung Cancer Res 2020; 9: 294–305.

51.

Zhu

Liu

, et al. WWC proteins mediate LATS1/2 activation by Hippo kinases and imply a tumor suppression strategy. Mol Cell 2022; 82: 1850–1864. e7.

52.

Lei

Shi

, et al. EML4-ALK Fusion gene in non-small cell lung cancer. Oncol Lett 2022; 24: 277.

53.

Cheon

Kwon

. Molecular anatomy of the EML4-ALK fusion protein for the development of novel anticancer drugs. Int J Mol Sci 2023; 24: 5821.

54.

Elshatlawy

Sampson

Clarke

, et al. EML4-ALK Biology and drug resistance in non-small cell lung cancer: A new phase of discoveries. Mol Oncol 2023; 17: 950–963.

55.

Zhang

Nagasaka

Zhu

, et al. Going beneath the tip of the iceberg. Identifying and understanding EML4-ALK variants and TP53 mutations to optimize treatment of ALK fusion positive (ALK+) NSCLC. Lung Cancer 2021; 158: 126–136.