CD4⁺CXCL13⁺ exhausted T cells drive immune microenvironment divergence in synchronous double primary lung adenocarcinoma with different degrees of invasiveness

Introduction

Double primary lung adenocarcinoma (DPLA) is defined as the presence of two lung adenocarcinomas of independent origin in the lung tissue of the same patient (1,2). In contrast to solitary lung adenocarcinoma, each lesion in DPLA may exhibit distinct stages of development and varying degrees of tumor invasion. Previous studies have indicated that the tumor microenvironment (TME) of solitary lung adenocarcinoma may demonstrate heterogeneity at different stages of development and degrees of tumor invasion, primarily reflected in differences in gene mutation spectra, immune cell infiltration patterns, and metabolic characteristics across lesions (3,4).

As a central immune component of the TME, tumor-infiltrating lymphocytes (TILs) play a critical role in regulating the dynamic balance between anti-tumor immune responses and immunosuppressive networks (5,6). In recent years, advances in single-cell sequencing and spatial multi-omics technologies have revealed the functional diversity of TILs within the TME, their interactions with tumor clonal evolution, and their impact on therapeutic responses (7). For instance, CD4⁺ TILs include subsets such as helper T cells (Th; e.g., Th1 and Th17), regulatory T cells (Tregs), and follicular helper T cells (Tfh), with their differentiation states being strictly regulated by tumor cell signals, metabolic microenvironments, and the extracellular matrix (ECM) (8). Extensive research has also been conducted on the role of CD8⁺ TILs in the development of lung tumors and their effects on immunotherapy response (9-11). Systematic investigations into the infiltration of TILs in multiple primary lung adenocarcinoma (MPLA) lesions could extend understanding of the tumor development mechanism (12). Moreover, a more in-depth exploration of the CD4⁺ and CD8⁺ TILs in the tumor immune microenvironment (TIME) could reveal sophisticated biomarkers, which hold promise in pinpointing patient groups that would respond positively to existing immune checkpoint blockade therapies, and would also facilitate the discovery of new targets for therapeutic intervention.

Research has reported that solitary lung adenocarcinoma lesions with different histological grades show distinct microenvironmental patterns (13). However, studies on the TME in synchronous DPLA (sDPLA) lesions are limited. From a clinical perspective, the management of sDPLA remains challenging. Individual lesions within the same patient often display heterogeneous invasiveness and biological behavior, yet current surgical and surveillance strategies are largely extrapolated from solitary lung adenocarcinoma. Thus, this study adopted a multi-omics approach to elucidate the characteristics of the TME, with a particular focus on TILs, in sDPLA with different degrees of infiltration. To exclude the interference of intrapulmonary metastasis, and investigate the characteristics and differences in the TIME of lung adenocarcinoma at different stages of development and degrees of tumor invasion, we conducted a comparative analysis between the invasive and non-invasive lesions in patients with sDPLA. We performed tumor driver gene detection, bulk RNA sequencing (RNA-seq), single-cell RNA-seq (scRNA-seq), multiplex immunofluorescence, and flow cytometry on all paired invasive and non-invasive sDPLA lesions. We present this article in accordance with the MDAR reporting checklist (available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-2025-1-1489/rc).

Methods

Patient cohort

We collected samples from a cohort of patients with sDPLA who underwent surgical resection at the Department of Thoracic Surgery, Ruijin Hospital, Shanghai Jiao Tong University School of Medicine between 2023 and 2024. Each patient had at least two distinct lesions pathologically confirmed as minimally invasive adenocarcinoma (MIA) and invasive adenocarcinoma (IAC), enabling lesion-level matched comparisons. Patients with any history of prior anti-cancer treatment were excluded from the study.

For the next-generation sequencing (NGS) genomic profiling analysis, a total of 14 paired lesions (seven MIA and seven IAC) from seven patients were analyzed. For the scRNA-seq, nine lesions from three patients, including three MIA, three IAC, and three matched normal lung (NL) tissues, were analyzed. For multiplex immunohistochemistry (mIHC), paraffin-embedded tissue sections from all three patients in the scRNA-seq cohort and one patient in the NGS cohort (P1) were used. For flow cytometry validation, fresh tumor and NL tissues were collected from the same three patients who underwent scRNA-seq profiling. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. This study was approved by the Ethics Committee of Ruijin Hospital, Shanghai Jiao Tong University School of Medicine [No. Ruijin-2021(219)]. Informed consent was taken from all the participants.

Single-cell suspension preparation

Under sterile conditions, the freshly collected tissues were washed twice with pre-cooled Roswell Park Memorial Institute (RPMI) 1640 medium containing 0.04% bovine serum albumin (BSA). The tissues were then mechanically fragmented into approximately 0.5 mm³ fragments using surgical scissors and transferred into freshly prepared enzymatic digestion solution. The digestion mixture contained RPMI 1640 (Corning Inc., Corning, NY, USA; cat. No. 10-040-CVR), 0.04% BSA (MACS, Miltenyi Biotec, Bergisch Gladbach, Germany; cat. No. 1000076), and 0.2% collagenase II (Gibco, Thermo Fisher Scientific, Waltham, MA, USA; cat. No. 17101015). The samples were incubated at 37 ℃ for 30–60 minutes with gentle inversion every 5–10 minutes.

The digested cell suspension was filtered through a 40-µm cell strainer (Falcon, Thermo Fisher Scientific; cat. No. 352340) 1–2 times. The filtrate was centrifuged at 300 ×g for 5 minutes at 4 ℃. The cell pellet was re-suspended in appropriate medium, mixed with an equal volume of red blood cell lysis buffer (Miltenyi Biotec; cat. No. 130-094-183), and incubated at 4 ℃ for 10 minutes. After centrifugation at 300 ×g for 5 minutes, the supernatant was discarded. The pellet was washed once with medium followed by another centrifugation at 300 ×g for 5 minutes, and the final supernatant was removed. Finally, the cells were re-suspended in 1 mL of RPMI 1640 (Conring Inc.; cat. No. 10-040-CVR) with 0.04% BSA. Single-cell suspension concentration and cell viability were then evaluated using the Luna-FL cell counter (Logos Biosystems, Anyang, Korea) or the trypan blue staining method.

scRNA-seq library construction

The freshly prepared single-cell suspension was adjusted to a concentration of 700–1,200 cells/µL. Library preparation and loading were performed using the MobiCube High-Throughput Single-Cell 3' Transcriptome Set, version 2.1 (cat. No. PN-S050200301), in accordance with the manufacturer’s instructions. The constructed libraries were sequenced on the Illumina Nova 6000 PE150 platform for high-throughput sequencing.

scRNA-seq data processing

The FASTQ files were processed and aligned to the GRCh38 human reference genome using MobiVision software (version 3.2) from 10× Genomics, with the unique molecular identifier (UMI) counts summarized for each barcode. The UMI count matrix was then analyzed using Seurat (14) (version 4.0.0) R package. To remove low-quality cells and likely multiplets, any cells that met the following filtering criteria were excluded: (I) number of detected genes <200; (II) total UMI count <1,000; (III) log₁₀ genes per UMI <0.7; (IV) percentage of mitochondrial UMI reads >10%; and/or (V) percentage of hemoglobin UMI reads >5%. Subsequently, the DoubletFinder package (version 2.0.3) was used to identify potential doublets. To obtain the normalized gene expression data, library size normalization was processed using the NormalizeData function. Gene expression data were normalized at the single-cell level using the global-scaling normalization approach implemented in Seurat (“LogNormalize”), in which expression values were normalized by total cellular counts, scaled by a factor of 10,000, and subsequently log-transformed. Highly variable genes were identified by selecting the top 2,000 features based on mean expression and dispersion using the FindVariableGenes function with the FastExpMean and FastLogVMR methods. Dimensionality reduction was conducted through principal component analysis (PCA) using the RunPCA function. Cell clustering was then performed using a graph-based algorithm implemented in the FindClusters function, grouping cells according to transcriptional similarity. For visualization, cells were projected into a two-dimensional space using Uniform Manifold Approximation and Projection (UMAP) via the RunUMAP function. Cluster-specific marker genes were identified using the FindAllMarkers function with the presto statistical test. Differentially expressed genes (DEGs) were selected using the function FindMarkers (test.use = presto). Bonferroni correction was applied for multiple testing adjustments of P values. A P value <0.05 and |log₂fold change| >0.58 was set as the threshold for significant differential expression. Gene Ontology (GO) enrichment and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses of the DEGs were performed using R (version 4.0.3) based on the hypergeometric distribution. The sequencing and bioinformatics analyses were performed by OE Biotech Co., Ltd. (Shanghai, China).

InferCNV analysis

Initial copy number variations (CNVs) for each region were estimated by inferCNV (version 1.0.4) (15) R package. CNV profiles were inferred from scRNA-seq data based on gene expression levels, using a threshold of 0.1. Genes were ordered according to their chromosomal positions, and smoothed expression values were obtained by calculating a moving average across a sliding window of 101 consecutive genes. The resulting expression values were then mean-centered to achieve a zero baseline. Epithelial cells were designated as putative malignant cells, whereas all other cell types were considered non-malignant reference populations. After noise reduction, tumor cells were identified according to CNV scores and hierarchical clustering patterns generated by inferCNV.

Cell clustering and annotation

Clustering was performed using the Louvain algorithm with the resolution parameters optimized to capture biologically meaningful subpopulations. The cell types were annotated based on canonical marker genes and referenced scRNA-seq atlases. Major immune cell types were identified including T/natural killer (NK) cluster CD4⁺ T cells (CD3D, CD3E, and CD4), CD8⁺ T cells (CD3D, CD3E, CD8A, and CD8B), NK cells (e.g., GNLY and PRF1), B cells (MS4A1 and CD79A), myeloid cells (CD68, LYZ, and FCER1A), and epithelial cell adhesion molecules (EpCAMs). The subclusters were assigned identities [e.g., CD4⁺ exhausted T (Tex), CD8⁺ effector T (Teff), naïve T, and NK bright cells] based on the expression of the marker genes, which are listed in Table S1.

Differential abundance and marker analysis of T/NK subpopulations

The proportions of the T/NK subclusters in each lesion were computed and compared between the paired MIA and IAC lesions. Differential abundance testing was performed using the propeller method implemented in the speckle package (version 0.0.3), and further evaluated with paired t-tests. The DEGs between the MIA and IAC lesions in each subcluster were identified using the Wilcoxon rank-sum test, with Bonferroni correction for multiple comparisons. The top markers were visualized via heatmaps and dot plots.

Trend analysis of T/NK subpopulations

To evaluate dynamic changes in T/NK composition across patients and conditions, we applied soft clustering using the Mfuzz package (version 2.58.0). The input data consisted of normalized subpopulation proportions across lesions. Mfuzz identified clusters with similar distribution trends, revealing possible differentiation or phenotypic transitions.

Gene set enrichment analysis (GSEA) of subpopulations

For selected T cell subtypes (e.g., CD4⁺ Tex and CD8⁺ Teff cells), the DEGs between the MIA and IAC lesions were subjected to GO enrichment analysis using the GSEABase package (version 1.44.0), which was used to load the gene set file that was downloaded and processed from the KEGG database (https://www.kegg.jp/) and GO database (https://geneontology.org/). A GSEA (16) was conducted to complete the GO and KEGG term enrichment analyses.

Cell communication analysis by CellChat

The cell communication analysis was performed using CellChat (17) (version 2.1.2) R package. First, we imported the normalized expression matrix to create the CellChat object using the create CellChat function. Second, the data were preprocessed using the identify OverExpressedGenes, identify OverExpressedInteractions, and projectData functions with the default parameters. The compute CommunProb, filter Communication (min.cells =10), and compute CommunProbPathway functions were then used to determine any potential ligand-receptor interactions. Finally, the cell communication network was aggregated using the aggregateNet function.

Monocle2 pseudotime analysis

Pseudotemporal trajectories were reconstructed using the Monocle2 package (version 2.9.0) (18). Raw count matrices were converted from Seurat objects into CellDataSet objects via the importCDS function. Genes used for trajectory ordering were identified with the differentialGeneTest function, selecting features with q-values below 0.01 that were informative for pseudotime progression. Dimensionality reduction was subsequently carried out using the reduceDimension function, and cells were ordered along the inferred trajectories using the orderCells function with default settings. Temporal changes in gene expression along pseudotime were visualized using the plot_genes_in_pseudotime function.

Flow cytometry validation

Fresh tumor tissues from the three patients in the scRNA-seq cohort were processed into single-cell suspensions as described above. The samples were incubated with Fc block and stained with fluorochrome-conjugated antibodies targeting CD3, CD4, CXCL13, and TOX2. The dead cells were excluded using Live/Dead Fixable dyes. The dells were analyzed using a BD FACSymphony A5 flow cytometer, and the data were processed using FlowJo (version 10.8.1). The CD4⁺ Tex cells were defined as CD4⁺CXCL13⁺TOX2⁺. Their frequencies were calculated as a proportion of the total CD4⁺ T cells and compared between the MIA and IAC lesions.

mIHC

Formalin-fixed paraffin-embedded (FFPE) sections from four patients (P1–P4) in the bulk RNA-seq cohort were used for mIHC. The tissue sections (4 µm) were deparaffinized, rehydrated, and subjected to heat-induced epitope retrieval. Sequential staining was performed using the Opal 7-Color immunohistochemistry (IHC) Kit (Akoya Biosciences, Marlborough, MA, USA), and antibodies against CD3, CD4, CXCL13, and TOX2. Each staining round was followed by microwave antigen retrieval and fluorophore deposition. The slides were scanned using Vectra Polaris (Akoya Biosciences), and the image analysis was performed using inForm software. The cell phenotypes were defined by the co-expression of markers, and densities of CD4⁺ Tex (CD4⁺CXCL13⁺TOX2⁺) and CD8⁺ Tex (CD8⁺PD1⁺TOX2⁺) cells were quantified per mm² and compared between the MIA and IAC lesions.

NGS-based genomic profiling using a 40-gene panel

Genomic DNA was extracted from the FFPE specimen blocks using the AmoyDx FFPE DNA Kit (cat. No. 20150079, AmoyDx^®, Amoy Diagnostics, Xiamen, China) in accordance with the manufacturer’s instructions. DNA concentration was quantified using a Qubit fluorometer (Thermo Fisher Scientific). Fragment size was evaluated using the Agilent 2100 Bioanalyzer with the DNA HS Kit (cat. No. 5067-1504/5067-1511, Agilent Technologies, Santa Clara, CA, USA). DNA was fragmented into 200–250 bp fragments using a Covaris LE220 (Woburn, MA, USA).

Targeted NGS was performed on each patient’s tumor DNA using the AmoyDx^® HANDLE Classic Panel (Amoy Diagnostics), which targets 40 cancer-related genes and assesses microsatellite instability status. A complete list of the genes included in this panel is provided in Table S2.

Public dataset comparison

To validate our findings, bulk RNA-seq data from the Gene Expression Omnibus (GEO) public database were downloaded and processed. The GEO dataset GSE27719 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE27719) with solitary adenocarcinoma lesions was evaluated to assess the heterogeneity of the T cell signatures between the solitary and synchronous cohorts. Single-sample GSEA (ssGSEA) scores for CD4⁺ Tex signatures were calculated and compared across subgroups stratified by histologic features.

Statistical analysis

All the statistical analyses were performed using R (version 4.2.1). For paired comparisons between the MIA and IAC lesions, including immune cell densities from mIHC, flow cytometric proportions, gene expression values, and ssGSEA scores, we used paired two-tailed Student’s t-tests. P values <0.05 were considered statistically significant unless otherwise specified. The data were visualized using ggplot2 (version 3.4.0) and ggpubr (version 0.6.0).

Results

Clinical and pathological features of the sDPLA cohort

The current study conducted a comprehensive sequencing analysis on a cohort of 10 individuals diagnosed with sDPLA (Table 1).

NGS was performed on 14 paired lesions (seven IAC and seven MIA) from a subset of seven patients to detect distinct genetic origins between the MIA and IAC lesions. Subsequently, a scRNA-seq analysis and multi-omics validation were performed on nine lesions (three IAC, three MIA, and three NL) from three patients. The age of the 10 patients ranged from 48 to 72 years (median, 59.9 years). The diameters of the invasive lung cancer lesions ranged from 0.8 to 2.0 cm (median, 1.5 cm). The diameters of the non-invasive lesions ranged from 0.6 to 1.1 cm (median, 0.9 cm). In terms of the distribution of the lesions, all the enrolled patients had lesions located in different lobes. All the patients had early-stage lung adenocarcinoma (before stage IA2). Of the 14 paired lesions from the NGS cohort, 13 harbored EGFR mutations. Additionally, in all the paired IAC and MIA lesions, the distribution patterns of gene mutations differed between the MIA and IAC lesions from the same patient, indicating the independent origins of the MIA and IAC lesions. The follow-up period for all patients ranged from 8 to 37 months (median, 19.2 months). As of the latest assessment, none of the patients exhibited any clinical signs or symptoms indicative of tumor recurrence or metastasis.

Characterization of the TME landscape and T cell subsets in sDPLA

scRNA-seq was performed to characterize the TME in sDPLA, focusing on the T/NK cell populations and their subtypes across the NL, MIA, and IAC tissues. UMAP clustering (Figure 1A) revealed distinct T/NK cell populations in all three tissue types. To further explore the transcriptional signatures of these subpopulations, a heatmap was generated (Figure 1B) that highlighted the gene expression patterns of various T/NK cells. The marker genes of each T/NK subset are listed in Table S1. Notably, the CD4⁺ Tex and CD8⁺ Tex cells had unique gene expression profiles, with an elevated expression of Tex-related genes in the IAC lesions. Notably, the CD4⁺ Tex cells not only expressed CD4 together with canonical exhaustion markers such as PDCD1 and CTLA4, but also showed significantly increased expression of CXCL13, based on which, this subset was designated as CD4⁺ Tex CXCL13. The correlation analysis showed strong co-expression between Tex-related markers and other immune cell markers in the CD4⁺ Tex CXCL13 cells, which was more pronounced in the IAC than the MIA or NL tissues (Figure 1C). These findings suggest that CD4⁺ Tex CXCL13 cells in IAC play a significant role in immune evasion, potentially contributing to the suppression of anti-tumor immunity in the invasive TME.

Figure 1 Single-cell clustering, characterization, and trend analyses of T/NK cell subtypes in IAC, MIA, and NL tissues. (A) UMAP clustering of major T/NK cell populations (left) and detailed subpopulations (right) across NL, MIA, and IAC tissues. (B) Heatmap showing the gene expression patterns of various T/NK subpopulations across tissues, with a focus on the genes associated with exhaustion, activation, and inflammation in the CD4⁺ Tex CXCL13 and CD8⁺ Tex cells. The gene expression levels are scaled by row, where red indicates high expression and blue indicates low expression. (C) Correlation heatmap showing the relationship between the expression of the Tex-related genes across the T/NK cell subpopulations. (D) UMAP plot of T/NK cell distributions across the IAC, MIA, and NL tissues, highlighting the increased proportion of CD4⁺ Tex CXCL13 cells in the IAC tissues compared to the MIA and NL tissues. (E) Stacked bar plot showing the relative abundance of each T/NK cell subtype in the IAC, MIA, and NL tissues. (F) Mfuzz trend analysis of the T/NK cell subtypes in the IAC vs. MIA tissues, highlighting the increased infiltration of the CD4⁺ Tex CXCL13 cells in the IAC tissues. IAC, invasive adenocarcinoma; MIA, minimally invasive adenocarcinoma; NK, natural killer; NL, normal lung; Teff, effector T; Tem, memory T; Tex, exhausted T; Tprf, T proliferation; Treg, regulatory T cells; Trm, tissue-resident memory T; UMAP, Uniform Manifold Approximation and Projection.

A further analysis of the distribution of the T/NK subpopulations in the IAC, MIA, and NL tissues was conducted using UMAP plots (Figure 1D). The results confirmed that the CD4⁺ Tex CXCL13 cells were more prevalent in the IAC than the MIA and NL tissues, and showed a distinct enrichment in the invasive lesions. This finding was corroborated by a stacked bar plot (Figure 1E), which showed that the CD4⁺ Tex CXCL13 cells comprised a higher proportion of the T/NK cell population in the IAC tissues. Conversely, CD4⁺ and CD8⁺ memory T (Tem) cells were similarly distributed across all tissue groups, indicating that their abundance was less affected by the invasive characteristics of the tumor.

The Mfuzz trend analysis provided a deeper understanding of the temporal dynamics of the T/NK cell subtypes in the IAC and MIA tissues (Figure 1F). The analysis revealed that the CD4⁺ Tex CXCL13 cells exhibited a significant increasing trend in the IAC tissues compared to the MIA tissues, suggesting they may be more prominently recruited to invasive tumor sites. Conversely, the CD4⁺ Tem and NK T cells showed similar trends in both the IAC and MIA tissues, with no clear divergence between the two tissue types. This further underscores the potential role of the CD4⁺ Tex CXCL13 cells in immune suppression in IAC; the other T/NK subsets did not show the same degree of dynamic change.

Dynamic reprogramming of Tex cells in the transition from MIA to IAC

To investigate the dynamics of intercellular communication associated with immune infiltration across different stages of lung adenocarcinoma, we applied CellChat to predict ligand-receptor interactions among all the identified cell populations in the NL, MIA, and IAC tissues. As shown in Figure 2A, the interaction strength was higher in the IAC tissues than both the MIA and NL tissues, while the number of interactions in IAC tissues was greater than in the NL tissues but slightly lower than in the MIA tissues, reflecting the intrinsic heterogeneity of tumor progression. Moreover, network visualization further highlighted differences in cell-cell communication patterns between the MIA and IAC lesions (Figure 2B). In particular, the Tex cells, including both the CD4⁺ and CD8⁺ cell subsets, exhibited significantly increased interaction frequencies with tumor cells in the IAC group. The differential interaction analysis corroborated these findings (Figure 2C), revealing significantly enhanced communication between the Tex cells and multiple immune and stromal compartments in the IAC tissues compared with the MIA tissues. Building on our previous observation that Tex cells, especially the CD4⁺ subset, were enriched in the IAC tissues, we next examined the transcriptional profiles of the CD4⁺ and CD8⁺ Tex cells between the MIA and IAC tissues. The differential expression analysis revealed distinct transcriptional signatures between the two groups (Figure 2D). The KEGG pathway enrichment analysis demonstrated that the genes upregulated in the CD4⁺ Tex CXCL13 cells from the IAC tissues were significantly associated with key immune signaling pathways, such as JAK-STAT and PI3K-AKT signaling (Figure 2E). Similarly, the CD8⁺ Tex cells in IAC exhibited the activation of immune-regulatory and tumor-associated pathways, including ECM-receptor interaction and focal adhesion (Figure 2F). Collectively, these results suggest that the transition from MIA to IAC is characterized by enhanced infiltration and rewired intercellular interactions of Tex cells, which may foster the development of an immunosuppressive or remodeled TME.

Figure 2 Cell-cell communication landscape and functional enrichment analysis of CD4⁺ Tex CXCL13 and CD8⁺ Tex cells in MIA and IAC. (A) Bar plots showing the total number of predicted ligand-receptor interactions (left) and cumulative interaction strength (right) among cell populations in NL, MIA, and IAC groups as inferred by CellChat. (B) Network diagrams illustrating the number of interactions between the cell subtypes in the MIA (left) and IAC (right) tissues. Each node represents a specific cell subtype, and edges represent inferred ligand-receptor interactions. (C) Differential interaction network between IAC and MIA, highlighting changes in the number of cell-cell interactions across cell subtypes. Red and blue edges indicate increased or decreased interactions in IAC relative to MIA, respectively. (D) Heatmap of DEGs in the CD4⁺ and CD8⁺ Tex cells. (E) KEGG pathway enrichment analysis of CD4⁺ CXCL13⁺ Tex cells revealing distinct transcriptional signatures and pathway activation in the IAC vs. MIA tissues. (F) KEGG pathway enrichment analysis of CD8⁺ Tex cells revealing distinct transcriptional signatures and pathway activation in the IAC vs. MIA tissues. DEG, differentially expressed gene; IAC, invasive adenocarcinoma; KEGG, Kyoto Encyclopedia of Genes and Genomes; MIA, minimally invasive adenocarcinoma; NK, natural killer; NL, normal lung; Teff, effector T; Tem, memory T; Tex, exhausted T; Tprf, T proliferation; Treg, regulatory T cells; Trm, tissue-resident memory T.

Progressive differentiation and transcriptional programs of CD4⁺ Tex cells

Based on our initial clustering and subtype annotation, which revealed an enrichment of CD4⁺ Tex CXCL13 cells in the IAC compared to the MIA lesions, we next performed a pseudotime trajectory analysis to investigate the differentiation dynamics of the CD4⁺ T cells (Figure 3A). Using Monocle 2, we reconstructed a branched developmental trajectory of the CD4⁺ T cells from both lesion types, which could be divided into five distinct states (Figure 3B, states 1 to 5). These states reflected a progressive differentiation process from early to terminally differentiated cells, with activated and exhausted-like subpopulations—such as CD4⁺ Tem and CXCL13⁺ Tex-like cells—predominantly occupying states 4 and 5. Notably, the CD4⁺ T cells from the IAC lesions were more enriched in these terminal states, suggesting a more advanced differentiation landscape compared to the MIA microenvironment (Figure 3B-3E).

Figure 3 Pseudotime trajectory and gene module analysis of CD4⁺ T cell differentiation in MIA and IAC lesions of sDPLA. (A) Monocle 2-based pseudotime trajectory reconstruction of CD4⁺ T cells revealed a branched structure with continuous differentiation from the early to terminal states. The major CD4⁺ T cell subsets were annotated. (B) The cells along the trajectory were classified into five pseudotime states (states 1 to 5), representing distinct differentiation phases. (C) Distribution of CD4⁺ T cells across the five states in the NL, MIA, and IAC lesions, showing a relative increase in terminal states (states 4 and 5) in IAC. (D) Pseudotime coloring showed differentiation progression and its lesion-specific distribution, with IAC cells skewed toward later pseudotime values. (E) Representative terminal CD4⁺ T cell subclusters (e.g., memory T- and Tex-like cells) mapped onto the trajectory, highlighting their enrichment in terminal branches of IAC lesions. (F) A branched pseudotime heatmap of dynamically regulated genes identified three (modules 1–3) of four gene modules that were preferentially upregulated in the terminal states (state 4/5). (G) KEGG enrichment analysis of modules 1–3 genes. IAC, invasive adenocarcinoma; KEGG, Kyoto Encyclopedia of Genes and Genomes; MIA, minimally invasive adenocarcinoma; NK, natural killer; NL, normal lung; sDPLA, synchronous double primary lung adenocarcinoma; Teff, effector T; Tem, memory T; Tex, exhausted T; Tprf, T proliferation; Treg, regulatory T cells.

To further investigate the transcriptional programs driving the development of terminal CD4⁺ T cells, we identified genes dynamically regulated along pseudotime and performed branched heatmap clustering, revealing three (modules 1, 2, and 3) of four gene modules that were preferentially upregulated in the terminal states (Figure 3F, states 4 and 5). These modules represent distinct gene expression programs associated with terminal CD4⁺ T cell differentiation; specifically, module 1 was enriched for pathways such as focal adhesion, MAPK, and sphingolipid signaling, which may support structural and migratory adaptation in tissue-localized effector or Tex-like cells (Figure 3G); module 2 was associated with antigen processing and presentation, TCR signaling, and Th1/Th17 lineage programs, indicating continued immune activation and polarization; module 3 was most strongly enriched in pathways related to T cell exhaustion and immune regulation, including PD-1 signaling, cytokine-cytokine receptor interaction, and apoptosis-related signaling, reflecting a transcriptional signature of functional exhaustion. These results indicate that multiple exhaustion-associated transcriptional programs are involved in the terminal differentiation of CD4⁺ T cells in IAC lesions. The co-existence of distinct gene modules highlights the functional heterogeneity of this differentiation trajectory, which may contribute to shaping the immunosuppressive microenvironment of IAC.

Multiplex immunofluorescence analysis of T cell subtypes in MIA and IAC lesions

To validate the transcriptomic findings at the protein level, we performed mIHC on FFPE tissue sections from four patients with paired MIA and IAC lesions. A panel of T cell markers was used to define and quantify major T cell subtypes based on co-localized marker expression, including CD8⁺ T cells, CD4⁺ T cells, CD8⁺/CD4⁺ tissue-resident memory T (Trm) cells, and Tex cells. The cell densities (cells/mm²) of each subtype were computed from five high-power fields per lesion.

As shown in Figure 4A,4B, the densities of most T cell subtypes—including the total CD8⁺ and CD4⁺ T cells, Trm cells, and CD8⁺ Tex cells—did not differ significantly between the MIA and IAC lesions (all P>0.10). Conversely, the CD4⁺ Tex CXCL13 cells showed a notable trend toward enrichment in the IAC lesions (P=0.08), which suggests that this subset may preferentially infiltrate more invasive foci. To further assess the relative abundance of the CD4⁺ Tex CXCL13 cells, we quantified their proportion among the total CD4⁺ T cells. This analysis also demonstrated a consistent trend of higher CD4⁺ Tex cell ratios in the IAC tissues compared to the MIA tissues (P=0.10; Figure 4C). These findings support the results of our scRNA-seq and flow cytometry analyses, highlighting CD4⁺ Tex CXCL13 cells as a potentially distinct immunological feature associated with more invasive lesions in multifocal lung adenocarcinoma.

Figure 4 Multiplex immunofluorescence analysis of co-localized positive T cell subtypes in MIA and IAC lesions. (A) Representative mIHC staining images of key T cell subtypes in IAC and MIA lesions, magnification: 200x. The co-expression of canonical markers was used to define T cell subsets, including CD8⁺ T cells (CD3⁺CD8⁺), CD4⁺ T cells (CD3⁺CD4⁺), Trm cells (CD8⁺/CD4⁺CD45RO⁺CD103⁺), and Tex cells (CD8⁺/CD4⁺TIM3⁺). (B) Quantitative analysis of cell density (cells/mm²) for each T cell subtype in paired MIA and IAC samples (n=4). (C) Boxplots showing paired comparisons of the CD4⁺ Tex cell ratio (%) among total CD4⁺ T cells (left) and the CD4⁺ Tex cell density (right) between the MIA and IAC lesions. IAC, invasive adenocarcinoma; MIA, minimally invasive adenocarcinoma; mIHC, multiplex immunohistochemistry; Tex, exhausted T; Trm, tissue-resident memory T.

Flow cytometry confirms divergent T cell infiltration patterns in IAC and MIA

To assess the immune microenvironmental differences between the IAC and MIA lesions, we performed flow cytometric profiling of tumor-infiltrating T cells using paired IAC and MIA samples from three patients (n=3 pairs). NL tissues adjacent to the tumors were also analyzed as a control. As illustrated in Figure 5A, the T cell populations were first gated into CD4⁺ and CD8⁺ subsets. Within these compartments, we further characterized the Tex cells, defined by the co-expression of PD-1 and TIM-3, and Trm cells, defined by the co-expression of CD103 and CD45RO. A quantitative analysis of the T cell subset composition revealed overall comparable distributions across the IAC, MIA, and NL samples (Figure 5B). However, one notable difference was identified in the exhausted CD4⁺ T cells—the proportion of CD4⁺ Tex cells (among the total CD4⁺ T cells) was significantly higher in the IAC lesions than the matched MIA lesions (Figure 5C, P=0.040). This suggests that IAC lesions harbor a more pronounced immunosuppressive microenvironment, potentially driven by chronic antigen exposure or immune evasion mechanisms during tumor progression. No statistically significant differences were observed in the frequencies of the other T cell subsets, including the CD8⁺ Tex and Trm cell populations, between the IAC and MIA lesions. These findings highlight CD4⁺ Tex CXCL13 enrichment as a potential immunological signature of invasive transformation in early-stage lung adenocarcinoma.

Figure 5 Flow cytometric analysis of tumor-infiltrating T cell subsets in IAC and MIA lesions of synchronous MPLA. (A) Representative flow cytometry gating strategies used to identify major tumor-infiltrating T cell subsets in paired IAC, MIA, and NL tissues. CD4⁺ and CD8⁺ T cells were gated based on CD4 and CD8 expression (top row). Tex cells were defined by the co-expression of PD-1 and TIM-3 (middle row), and Trm cells were defined by the co-expression of CD103 and CD45RO (bottom row). (B) Quantitative comparison of the relative frequency of the T cell subsets in the CD4⁺ and CD8⁺ T cell compartments across the IAC, MIA, and NL samples. Data are presented as a percentage of each subset among the total CD4⁺ or CD8⁺ T cells, with the mean ± standard deviation shown for three matched patient cases. (C) Boxplot of CD4⁺ Tex cells between the IAC and MIA groups. Each dot represents one patient sample (n=3 pairs). FITC, fluorescein isothiocyanate; IAC, invasive adenocarcinoma; MIA, minimally invasive adenocarcinoma; MPLA, multiple primary lung adenocarcinoma; NL, normal lung; Teff, effector T; Tem, memory T; Tex, exhausted T; Treg, regulatory T cells; Trm, tissue-resident memory T.

Discussion

Advancements in imaging technology (e.g., high-resolution computed tomography) and the widespread adoption of early screening have increased the detection rate of multiple primary lung cancer (MPLC) patients. Epidemiological studies have reported that the incidence of MPLC has risen significantly from 0.2–0.5% to 8% in recent years (19). Data from the Japanese Lung Cancer Registry (n=18,978) further revealed that among 9,689 stage I non-small cell lung cancer patients, the proportion of MPLC cases showed a year-by-year increasing trend (20). Research data from the Guangdong Lung Cancer Institute [2005–2013] also confirmed that the proportion of MPLC among lung cancer patients has significantly increased compared to previous reports (21). Among MPLC patients, DPLA is a common pathological type. A study involving 482 lung cancer patients revealed that the incidence of sDPLA in primary lung cancer patients was 1.6% (22), significantly higher than the 0.67% incidence of triple primary lung cancer (21). In terms of the distribution of the different pathological types, a study by Lyu et al. (23) showed that among MPLC cases, lung adenocarcinoma accounted for 70.78%, a figure significantly higher than the 2.25% observed for squamous cell carcinoma. Further, a large-scale study based on the United States Surveillance, Epidemiology, and End Results (SEER) database (n=1,419) further established the dominant role of DPLA in MPLC (24). However, systematic research on patients with sDPLA, which is representative of sDPLA, remains relatively scarce. In light of this, this study conducted multi-omics detection and analyses of lesion tissues from a cohort of DPLA patients at our hospital to examine the complexity and diversity of the immune context of the TME. It also innovatively explored the characteristics and differences in the immune microenvironment of sDPLA lesions at different developmental stages and with varying degrees of invasiveness.

In this study, we employed scRNA-seq to analyze the immune landscape of DPLA by examining the immune infiltration in MIA and IAC lesions. Our findings revealed significant differences in immune cell populations, particularly in the T/NK cell compartments, highlighting a marked enrichment of CD4⁺ Tex CXCL13 cells in the IAC tissues compared to the MIA and NL tissues. This suggests that CD4⁺ Tex CXCL13 cells may contribute to the immune suppression observed in invasive tumors, a key feature of cancer progression. This study represents the first comprehensive multi-omics analysis of the immune microenvironment in DPLA lesions. By using paired lesions from the same patient, we effectively controlled for inter-individual genetic and treatment variability, thereby minimizing confounding factors. This design enabled a more precise investigation into the immune alterations occurring as MIA progresses to invasive IAC. Our findings offer novel insights into the immune evolutionary mechanisms of early-stage lung adenocarcinoma progression and provide a solid foundation for the development of targeted immunotherapeutic strategies (25).

Further, our research was the first to show that in dual-primary lesions with different infiltration degrees from the same patient, the CD4⁺ Tex CXCL13 level in the IAC group was higher than that in the MIA group. However, an analysis of solitary IAC and MIA lesions derived from different patients in the GEO public database (cohort GSE27719) revealed a different trend (Figure S1; Wilcoxon test, P=0.19). Taken together with our findings, this suggests a potential difference in the evolutionary pattern of the immune microenvironment during the progression from non-invasive to invasive stages between solitary adenocarcinomas and DPLA.

We also observed a significant enrichment of CD4⁺ Tex CXCL13 cells in IAC tissues compared to MIA and NL tissues. These CD4⁺ Tex CXCL13 cells exhibited high expression of immune checkpoint molecules such as PD-1 and TIM-3, and demonstrated impaired cytokine production, characteristics commonly associated with T cell exhaustion in chronic infections and cancer (26,27). The accumulation of these cells in IAC lesions suggests their potential role in immune evasion and tumor progression. These findings align with previous studies highlighting the importance of CD4⁺ Tex CXCL13 cells in the TME and their implications for cancer immunotherapy.

However, several limitations in this study need to be acknowledged. First, the small sample size (three patients with six paired lesions) might limit the generalizability of our findings, especially given the heterogeneity of the sDPLA microenvironment. Studies with larger cohorts need to be conducted to validate these results. Second, scRNA-seq captures transcriptional states but lacks functional validation. Although the CD4⁺ Tex CXCL13 cells were identified based on exhaustion-related genes, direct functional validation using in vitro co-culture systems was not performed. The rarity of sDPLA and the limited availability of fresh paired samples, together with technical challenges in establishing organoids from early-stage MIA lesions, constrained the feasibility of such experiments within the current study. Finally, technical constraints of scRNA-seq—such as dropout events and the loss of spatial context—may affect data interpretation. The integration of spatial transcriptomics or mIHC data in future studies will help address these gaps.

The identification of CD4⁺ Tex CXCL13 cells as a key immune subset in IAC highlights the need to further investigate their molecular mechanisms and functional roles. Future studies could explore their differentiation dynamics and plasticity using approaches such as pseudotime analysis and gene perturbation models to uncover the pathways driving exhaustion and suppression. Additionally, examining their interactions with CD8⁺ T cells, myeloid cells, and tumor-associated macrophages may provide deeper insights into immune regulation in the IAC microenvironment.

Conclusions

This study shed light on the complex immune landscape of DPLA and highlights CD4⁺ Tex CXCL13 cells as a potential target for therapeutic intervention. While our findings provide valuable insights, future studies with larger cohorts, integrated functional analyses, and multimodal technologies are required to further elucidate the role of immune cell subtypes in tumor progression and to develop effective immunotherapeutic strategies tailored to the unique immune microenvironment of DPLA.

Acknowledgments

We acknowledge the patients for their consent to use the tumor tissue for the study.

Footnote

Reporting Checklist: The authors have completed the MDAR reporting checklist. Available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-2025-1-1489/rc

Data Sharing Statement: Available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-2025-1-1489/dss

Peer Review File: Available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-2025-1-1489/prf

Funding: This study was supported by the National Natural Science Foundation of China (Nos. 82072557, 81871882, and 82372855), the National Key Research and Development Program of China (No. 2021YFC2500900), the Program of Shanghai Academic Research Leader from Science and Technology Commission of Shanghai Municipality (No. 20XD1402300), the Novel Interdisciplinary Research Project from Shanghai Municipal Health Commission (No. 2022JC023), the Interdisciplinary Program of Shanghai Jiao Tong University (No. YG2023ZD04), and the Shanghai Sailing Program (No. 21YF1427100).

Conflicts of Interest: All authors have completed the ICMJE uniform disclosure form (available at https://tlcr.amegroups.com/article/view/10.21037/tlcr-2025-1-1489/coif). The authors have no conflicts of interest to declare.

Ethical Statement: The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved. This study was conducted in accordance with the Declaration of Helsinki and its subsequent amendments. This study involves human participants and this study was approved by the Ethics Committee of Ruijin Hospital, Shanghai Jiao Tong University School of Medicine [No. Ruijin-2021(219)]. Informed consent was taken from all the participants.

Open Access Statement: This is an Open Access article distributed in accordance with the Creative Commons Attribution-NonCommercial-NoDerivs 4.0 International License (CC BY-NC-ND 4.0), which permits the non-commercial replication and distribution of the article with the strict proviso that no changes or edits are made and the original work is properly cited (including links to both the formal publication through the relevant DOI and the license). See: https://creativecommons.org/licenses/by-nc-nd/4.0/.

References

1. Kanemura H, Yokoyama T, Nakajima R, et al. The Tumor Immune Microenvironment Is Associated With Recurrence in Early-Stage Lung Adenocarcinoma. JTO Clin Res Rep 2024;5:100658. [Crossref] [PubMed]
2. Chen K, Liu A, Wang C, et al. Multidisciplinary expert consensus on diagnosis and treatment of multiple lung cancers. Med 2025;6:100643. [Crossref] [PubMed]
3. Won J, Chung YS, Jo SY, et al. Genomic classification of intrapulmonary metastasis and multiple primary lung cancer. Clin Transl Med 2025;15:e70463. [Crossref] [PubMed]
4. Takemura C, Yoshida T, Yoshida Y, et al. Unveiling the molecular and clinical risk landscape of second primary lung cancer in resected non-small cell lung cancer. Lung Cancer 2025;208:108750. [Crossref] [PubMed]
5. Kinoshita F, Takada K, Wakasu S, et al. Granzyme B (GZMB)-Positive Tumor-Infiltrating Lymphocytes in Lung Adenocarcinoma: Significance as a Prognostic Factor and Association with Immunosuppressive Proteins. Ann Surg Oncol 2023;30:7579-89. [Crossref] [PubMed]
6. Wang Y, Ding Q, Wei J. Current advances and future directions of combined ICIs and TILs in solid tumors. Cancer Lett 2026;636:218145. [Crossref] [PubMed]
7. Ge W, Wu W. Influencing Factors and Significance of Tumor-associated Macrophage Polarization in Tumor Microenvironment. Zhongguo Fei Ai Za Zhi 2023;26:228-37. [Crossref] [PubMed]
8. Imianowski CJ, Chen Q, Workman CJ, et al. Regulatory T cells in the tumour microenvironment. Nat Rev Cancer 2025;25:703-22. [Crossref] [PubMed]
9. Gueguen P, Metoikidou C, Dupic T, et al. Contribution of resident and circulating precursors to tumor-infiltrating CD8(+) T cell populations in lung cancer. Sci Immunol 2021;6:eabd5778. [Crossref] [PubMed]
10. Zhang L, Chen Y, Wang H, et al. Massive PD-L1 and CD8 double positive TILs characterize an immunosuppressive microenvironment with high mutational burden in lung cancer. J Immunother Cancer 2021;9:e002356. [Crossref] [PubMed]
11. Ma K, Xu Y, Cheng H, et al. T cell-based cancer immunotherapy: opportunities and challenges. Sci Bull (Beijing) 2025;70:1872-90. [Crossref] [PubMed]
12. Torri M, Sandell A, Al-Samadi A. The prognostic value of tumor-infiltrating lymphocytes in head and neck squamous cell carcinoma: A systematic review and meta-analysis. Biomed Pharmacother 2024;180:117544. [Crossref] [PubMed]
13. Bischoff P, Trinks A, Obermayer B, et al. Single-cell RNA sequencing reveals distinct tumor microenvironmental patterns in lung adenocarcinoma. Oncogene 2021;40:6748-58. [Crossref] [PubMed]
14. Hao Y, Hao S, Andersen-Nissen E, et al. Integrated analysis of multimodal single-cell data. Cell 2021;184:3573-3587.e29. [Crossref] [PubMed]
15. Puram SV, Tirosh I, Parikh AS, et al. Single-Cell Transcriptomic Analysis of Primary and Metastatic Tumor Ecosystems in Head and Neck Cancer. Cell 2017;171:1611-1624.e24. [Crossref] [PubMed]
16. Subramanian A, Tamayo P, Mootha VK, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A 2005;102:15545-50. [Crossref] [PubMed]
17. Jin S, Guerrero-Juarez CF, Zhang L, et al. Inference and analysis of cell-cell communication using CellChat. Nat Commun 2021;12:1088. [Crossref] [PubMed]
18. Trapnell C, Cacchiarelli D, Grimsby J, et al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nat Biotechnol 2014;32:381-6. [Crossref] [PubMed]
19. Tan L, Yin J. Diagnosis and Treatment for Multiple Primary Lung Cancer. Zhongguo Fei Ai Za Zhi 2018;21:185-9. [Crossref] [PubMed]
20. Shintani Y, Okami J, Ito H, et al. Clinical features and outcomes of patients with stage I multiple primary lung cancers. Cancer Sci 2021;112:1924-35. [Crossref] [PubMed]
21. Li F, Zhong WZ, Niu FY, et al. Multiple primary malignancies involving lung cancer. BMC Cancer 2015;15:696. [Crossref] [PubMed]
22. Shankar PS. Synchronous dual primary lung carcinomas. J Am Geriatr Soc 1981;29:328-30. [Crossref] [PubMed]
23. Lyu M, Luo L, Zhou L, et al. Emerging trends in the coexistence of primary lung Cancer and hematologic malignancy: a comprehensive analysis of clinicopathological features and genetic abnormalities. Cancer Cell Int 2024;24:84. [Crossref] [PubMed]
24. Wang G, Zeng Y, Zheng H, et al. Prognostic analysis and clinical characteristics of dual primary lung cancer: a population study based on surveillance, epidemiology, and end results (SEER) database. Gen Thorac Cardiovasc Surg 2022;70:740-9. [Crossref] [PubMed]
25. Liu J, Gu M, Xue Y, et al. Clinical significance of PD-L1 expression and CD8-positive tumor-infiltrating lymphocytes in patients with cavitary lung adenocarcinoma. Oncologie (Tech Science Press) 2021; [Crossref]
26. McManus DT, Valanparambil RM, Medina CB, et al. An early precursor CD8(+) T cell that adapts to acute or chronic viral infection. Nature 2025;640:772-81. [Crossref] [PubMed]
27. Galluzzi L, Smith KN, Liston A, et al. The diversity of CD8(+) T cell dysfunction in cancer and viral infection. Nat Rev Immunol 2025;25:662-79. [Crossref] [PubMed]

(English Language Editor: L. Huleatt)