The single-cell datasets HRA000910 and HRA000363, which contain 13 samples, were downloaded from the National Genomics Data Centre (NGDC) Genome Sequence Archive (GSA). The raw data from each sample were subjected to sequencing library demultiplexed and aligned to the human reference genome GRCh38, followed by quantification of UMI counts using 10x Genomics Cell Ranger (v3.1.0). After cell identification using DropletUtils (v1.6.1), the quantified expression matrix was QC’d using the Seurat package (v5.1.0) for R software (v4.4.1). Cells with less than <15% mitochondrial gene content, together with a total gene count >300 and gene expression level between 500 and 15,000 and expressed in at least three cells were retained.

After quality control, the data were filtered for γδ T cells based on PTPRC, CD3D, CD3E, TRDV1, and TRDV2 gene expression levels greater than 1 using the subset function. The gene expression data were normalized using the NormalizeData function, which employs the default parameters of the Seurat package, the 2000 genes with the largest variations were selected using the FindVariableFeatures function, the gene expression was normalized using the ScaleData function and finally the data were linearly downscaled by PCA using the RunPCA function. The FindNeighbors function was performed by selecting 1:15 PCs based on the ElbowPlot function, and the FindClusters function was performed with a resolution of 0.6. For visual clustering, the RunUMAP function was used with the same number of PCs to generate the Uniform Manifold Approximation and Projection (UMAP) algorithm.

Marker genes for each cluster relative to all other clusters were determined using the FindAllMarkers, employing a Wilcoxon test with a p-value <0.05, Bonferroni correction test with p-adj < 0.05 and a differential expression threshold of 0.25, where selected marker genes were expressed in at least 25% of the target cell subpopulations. Each cluster was determined to specifically express the top 3 genes using the COSG package (v0.9.0). The method for gene marker identification based on cosine-based values. Cluster was labeled by manually matching typical cell marker genes with the algorithmically calculated genes that are characteristic of each cluster.

G1/S and G2/M phase signature genes were analyzed using the CellCycleScoring function to predict the temporal cell cycle phases of γδ T cell subsets.

Cell trajectories were analyzed using the monocle package (v2.32.0) for the Vδ1 and Vδ2 cell clusters, respectively. Seurat data were analyzed downstream for differential gene expression using a negative binomial distribution model with the newCellDataSet. Normalization was performed using the estimateSizeFactors and estimateDispersions functions, and feature selection was performed with the differentialGeneTest function to filter out 2000 genes that had a significant impact on cell trajectories. The reduceDimension function was employed to reduce the dimensionality of the data using the DDRTree algorithm, and the orderCells function was used to sort the cells according to their developmental trajectory. The BEAM function was used to model and analyze the branch point-dependent genes within the cell developmental trajectory to obtain the significance score for each gene.

Enrichment analyses were performed for each cell cluster using the clusterProfiler package (v4.12.6) and the compareCluster function to select the enrichGO Biological Processes (BP) GO database. P-values were corrected using the Benjamini & Hochberg (BH) method. Additionally, pathway enrichment analysis was performed for each cell cluster using the ReactomePA package (v1.48.0) with the compareCluster function to select the Reactome database.

We obtained fresh mononuclear cells from peripheral blood (PBMC) and pleural fluid (PFMC) of five TPE patients for single-cell sequencing (9), and used the Seurat R package (10) for quality control and screening of the raw sequencing data for γδ T cells (TRDC) ( Figures 1A, B ). Ultimately, 959 γδ T cells from peripheral blood and 1259 γδ T cells from pleural fluid were included in the study. The screened γδ T cells were linearly downscaled using principal component analysis (PCA), and the top 15 principal components were then selected for downscaled clustering based on the fragmentation map ( Figures 1C, D ). The results of the PCA demonstrated that there is a distinction in the transcriptomic profiles exhibited by PBMCs and Pleural Fluid Mononuclear Cells (PFMCs) in TPE patients.

The unsupervised clustering algorithm yielded a total of 7 γδ T cell subpopulations ranging from 0 to 6, and the top 10 gene heatmaps for each cluster were as follows ( Figure 2A ). Using the non-linear clustering UMAP method (11) to visualise the γδ T cell subpopulations, we classified γδ T cells into 7 clusters ( Figure 2B ). Consistent with previous PCA results, UMAP visualization results showed that PBMCs and PFMCs from TPE patients have distinct clusters ( Figure 2C ). Clusters C1, C4, C5, and C6 had a higher proportion of PBMCs with 36.1%, 17.5%, 16.4%, and 14%, respectively, while PFMCs were more predominant in clusters C0 and C3 with 70.5% and 16.5%, respectively, and cluster C2 is the same in both ( Figures 2D, E ).

According to the expression of γδ T cell signature molecules (TRDV1, TRDV2), C3 and C6 were identified as Vδ1 cells, while C0, C1, C2, C4, and C5 were classified as Vδ2 cells ( Figure 3A ). The scoring calculation of the cell cycle status (12) for each γδ T cell cluster showed that all seven γδ T cell subpopulations exhibited the same cell cycle status ( Figures 3B, C ).

The highest percentage of effector Vδ2 cell cluster (C1) (36.1%), along with elevated expression of cytotoxicity-related genes (PRF1, GZMA, GZMB, CX3CR1, and LAMP1) and effector transcription factors (ID2, TBX21, and ZEB2) were found in PBMCs from TPE patients. The highest percentage (70.5%) of tissue-resident Vδ2 cell cluster (C0) was found in PFMCs, with high expression of tissue-resident related genes (ITGAE and CD69) ( 13), and its high expression of IL7R was speculated to be a possible tissue-resident memory Vδ2 cell subpopulation. Central memory Vδ2 cells (C2) was found in both PBMCs and PFMCs, expressing high levels of CD27, IL7R and mitochondrial oxidative phosphorylation genes (ATP5E, ATP5D, etc.) (14), as well as cytotoxicity-related genes (PRF1, GZMA). In PFMCs, there was a higher percentage (16.5%) of naive Vδ1 cell cluster (C3), characterized by elevated expression of resting stemness-related genes (TCF7, SELL, LEF1, BACH2, and IKZF2) ( 15) and almost no expression of effector-related genes. In PBMC, the INPP4B⁺ Vδ2 cell cluster (C4) accounted for 17.5% and was characterized by high expression of the inositol phosphatase INPP4B, along with elevated expression of metallothionein-related genes (MT1E and MT2A) (16) ( Figures 3D, E ). The MT⁺ Vδ2 cell cluster (C5) accounted for (16.4%) in PBMCs, which was characterized by high expression of metallothionein-related genes (MT1E, MT2A, MT1G, and MT1X) ( 17, 18). The NK-like Vδ1 cell cluster (C6) accounted for (14%) in PBMCs, which was characterized by high expression of NK-related genes (GNLY, KLRC2, and KLRF1) ( 19, 20), as well as high expression of cytotoxicity-associated genes (PRF1, GZMA, GZMB, CX3CR1, and LAMP1) and effector transcription factors (TBX21 and ZEB2) ( Figures 3D, E ). Taken together, these results reflect the phenotypic and functional heterogeneity of γδ T cells in PBMCs and PFMCs of TPE patients.

We analyzed the differentiation trajectories of Vδ2 cell subpopulations using pseduotime analysis (21) and found that the tissue-resident Vδ2 cell cluster (C0) was predominantly located in the lower left branch, whereas the effector Vδ2 cell cluster (C1) and the central memory Vδ2 cell cluster (C2) were located in the upper left and right branches, respectively ( Figures 4A–C ). Based on pseduotime calculations, it was speculated that the tissue-resident Vδ2 cell cluster (C0) mainly differentiated into the effector Vδ2 cell cluster (C1) and the central memory Vδ2 cell cluster (C2), whereas the INPP4B⁺ Vδ2 cell cluster (C4) and the MT⁺ Vδ2 cell cluster (C5) were mainly located in the intermediate transition state ( Figures 4A–C ). In contrast, the proposed time series analysis of the differentiation trajectory of the Vδ1 cell cluster showed that there was no branching of the differentiation trajectory, and the naive Vδ1 cell cluster (C3) differentiated into the NK-like Vδ1 cell cluster (C6) ( Figures 4D, E ). We speculate that the Vδ2 cluster in situ in the lungs of TPE patients stimulated by TPE infection has two main directions of differentiation, one into the peripheral blood to differentiate into a Vδ2 cluster with effector functions, and the other to become a central memory Vδ2 cluster to cope with TPE in the long term. Overall, these results reveal the differentiation trajectory of γδ T cell subsets.

Effector Vδ2 cells are the mainstay of the TPE immune response, and further study of the functions and characteristic genes of these effector Vδ2 cells can enhance our understanding of the TPE immune response of γδ T cells. Enrichment analysis of biological processes in the GO database22 revealed that a cluster of effector Vδ2 cells exhibated high activation and cytotoxic functions ( Figure 5A ; Supplementary Material 3 ). Enrichment analysis of the Reactome database23 revealed that a cluster of effector Vδ2 cells had specific FCGR-activated signaling pathways ( Figure 5B ; Supplementary Material 3 ). We further analyzed the characteristic genes of each cluster using the COSG package24, in which the characteristic gene of the effector Vδ2 cell subpopulation was FCGR3A, and the violin plots also showed that FCGR3A was highly expressed in the effector Vδ2 cell subset ( Figures 5C, D ). The expression of FCGR3A increased with the differentiation of Vδ2 cells and reached the highest level in the effector Vδ2 cell subset, and at the same time, the expression of FCGR3A with the other cytotoxicity-related genes (PRF1, GZMB and GZMH) and effector transcription factors (TBX21 and ZEB2) had similar pseudotime expression patterns ( Figure 5E ). We used a regression algorithm to search for pseudotime differential genes in the effector Vδ2 cell cluster, and the heatmap showed FCGR3A as a branching-dependent gene in the effector Vδ2 cell subset (q-value < 0.001) ( Figure 5F ; Supplementary Material 1 ).

In conclusion, CD16 (encoded by the FCGR3A gene) was specifically expressed in a subset of effector Vδ2 cells in the peripheral blood of TPE patients and correlated with cytotoxic and effector transcription factor expression.