Single-cell transcriptomic and TCR repertoire profiles of TFHL

We first analyzed scRNA-seq data of 68,738 LN and 98,891 PB cells collected from 14 TFHL patients and controls (Fig. 1A; Table S1). Unsupervised clustering and annotation by canonical markers revealed four main cell types: CD4/CD8+ T, natural killer (NK; NK/gamma-delta [γδ] T for PB), B, and myeloid cells (Fig. 1B, Fig. S2A).

Subsequently, we analyzed LN scTCR-seq data to identify tumor cells and found unique TCR clonotypes in 93.7% of T cells (Fig. S2B, C). Major and minor tumor-cell clonotypes were defined based on clonality (Fig. 1C, Fig. S2D; Table S3). The median proportion of tumor cells in all LN MNCs was 14.1% (range 5.7–83.8%; Table S4). In 7 of 9 TFHL LNs, tumor cells were clones within a patient, suggesting single-TCR clonal expansion (Fig. 1C, Fig. S2D). Contrastingly, in AITL4, two clones with identical TCRβ chains and differing TCRα chains expanded (Fig. 1C; Table S3). Notably, 31.3% (2652 cells) of tumor cells expressed only one productive α or β chain, contrasting dogma that αβ T cells physiologically express a pair of productive α and β chains (Fig. 1C, Fig. S2D). In three TFHLs (AITL7, AITL8, and AITL12), the single-TCR clonotype was detected in most tumor cells (Fig. S2E). Moreover, in these TFHLs, FCM analysis revealed the expansion of abnormal CD4+ PD1high T-cell populations lacking cell-surface expression of TCRα and β chain-interacting CD3ε [23], which was not observed in TFHLs with physiological paired-chain TCRs (Fig. S2F, G, and data not shown). Tumor cell infiltration was also identified in 12 of 16 PB samples (0.24–74.0% in PB MNCs; Table S4), five of which showed expansion of single-chain TCR tumor cells with CD3ε- CD4+ PD1high T-cell populations by FCM, consistent with the LN data (Fig. S2F and S3A–C).

Genomic characterization of TFHL

We next analyzed somatic mutations using bulk WES data to construct genomic profiles of TFHL, finding 397 total somatic mutations in 14 TFHL patients. This included recurrent TET2 mutations in 11 (78.6%), G17V in 6 (42.9%), DNMT3A in 6 (42.9%), and IDH2 R172 in 2 (14.3%), comparable to previous reports [24,25,26,27] (Fig. S4A; Table S5). Furthermore, we analyzed CNVs and found that gain of chromosome (chr) 5 was most frequently observed (35.7%; 5/14), and gains of chr7/7q, 19, 21, or 22q were also detected in ≥2 of 14 cases (14.3%) with co-occurrence of chr5 and 21 gain in 1 case (7.1%), consistent with previous reports [28] (Fig. S4B).

Inter-patient and intra-tumor TFHL heterogeneity

To investigate tumor cell characteristics and heterogeneity, we sorted 16,358 tumor cells from scRNA-seq data into five clusters by unsupervised clustering (C0–C4; Fig. 2A, Fig. S5A), finding consistent inter-patient heterogeneity (Fig. S5B). Notably, TFH marker [4] expressions were highly heterogeneous between clusters (Fig. 2B, Fig. S5C; Table S6). PDCD1 and ICOS were globally expressed, while MME (encoding CD10) and other TFH markers were predominantly expressed in C0 and C3–4, respectively (Fig. 2B). Subsequently, we conducted gene set variation analysis (GSVA) [29] on each cell using TFH and non-TFH effector T cell (TEFF) signatures derived from bulk RNA-seq data (GSE58596 [30]). This revealed enrichment of the TFH signature in C3–4 and the TEFF signature in C0–2 (Fig. 2C). Moreover, cell-cycle scores revealed that C0, C1, and C3 were primarily in G1 stage, but C2 and C4 were exclusively in G2M or S stages, indicating active cell proliferation in C2 and C4 (Fig. S5D).

Fig. 2: Single-cell analysis of tumor cells from TFHL LN and PB samples.

A UMAP plots of LN and PB tumor cell subclusters (top). Cells are shown separately for each clinical status (bottom). B The expression levels of known T follicular helper (TFH) cell markers for each cluster. Adjusted P values are calculated by the Wilcoxon rank-sum test with Bonferroni correction across all genes expressed by tumor cells and shown only when there is a significant difference. C Feature (left panel) and violin (right panel) plots of gene set variation analysis (GSVA) enrichment score for TFH (top) or CD4+ effector T cell (TEFF, bottom) signatures in tumor cells. In the violin plots, the dotted lines represent mean enrichment scores across all clusters. The boxplots show the median (center line), mean (center dot), interquartile range (box limits), and minimum to max values (whiskers) for each group. Adjusted P values are calculated by the pairwise Wilcoxon rank-sum test for each cluster against the mean value of all clusters. *P < 5.0 × 10−2, ****P < 1.0 × 10−4.

Next, we compared LN and PB tumor cells, revealing that C3–4 were predominant in LNs while C1–2 were more abundant in PB (Fig. S5E). Moreover, differentially expressed gene (DEG) analysis and GSVA revealed that LN tumor cells exhibited relatively higher expression of TFH markers whereas PB tumor cells showed high expression of cytotoxicity-associated genes [31] (Fig. S5F, G). FCM analysis confirmed significantly higher PD1 expression in LN tumor cells versus those of PB (Fig. S5H, I). Subsequent trajectory analysis using Monocle3 [32] inferred that LN tumor cells terminally differentiated starting from C0 through C3 to C4, whereas PB tumor cells differentiated from C0 to C2 via C1 (Fig. S5J).

We next performed DEG and gene ontology (GO) analyses between tumor cells with single-chain and paired-chain TCRs (Fig. S5K). Unexpectedly, both types had activated TCR and cytokine signaling pathways (Fig. S5L). Notably, genomic mutations known to activate TCR pathways (G17V, CTNNB1), or those in TCR signaling-related genes (LCK, KRAS) [33], were detected in 4 of 5 cases with expanded single-chain TCR clonotypes (Table S5).

Estimation of tumor evolution at single-cell resolution

To assess mutational evolution using scRNA-seq tumor cell data, we reanalyzed the initially detected 11 somatic mutations and focused on mutations with a total read coverage ≥100×, identifying seven mutations present in over 20 cells, of which only G17V was detected in multiple patients (Fig. S6A–C; Table S7). In TFHL patients where G17V was initially identified by WES, G17V was detected in 86.7% of tumor cells while 9.4% of the cells lacked mutant reads and 3.9% had no coverage (Fig. 3A). Considering the frequent allelic dropout in scRNA-seq [34], cells without mutant reads or no coverage status were categorized as “unknown”. Conversely, cells from TFHL patients without G17V mutations by WES were dubbed “G17V wild-type (WT)” (Fig. 3A). GSVA, DEG, and GO analyses revealed enrichment of TFH signatures and pathways related to the adaptive immune system and B-cell activation, including CD40LG, in G17V mutant cells compared to G17V WT cells, consistent with our previous report [18] (Fig. 3A, Fig. S6D, E). In AITL4LN/PB, with two distinctive tumor clones featuring identical TCRβ chains and different TCRα, LRRC41 mutations were detected only in tumor-clone 1 (Fig. S6F). Based on UMAP visualization, these cells with different clones exhibited different gene expression profiles, indicating clonal evolution across divergent genetic/transcriptomic directions (Fig. S6G).

Fig. 3: Estimation of genomic alternations at the single-cell level.

A Distribution of RHOA G17V (G17V) mutant cells (left panel) and comparison of GSVA enrichment score for the TFH signature between G17V mutant (red) and wild-type (WT, white) cells (right panel) in LN and PB tumor cells. In the violin plots of the right panel, adjusted P values are calculated for all tumor cells (left) and each tissue (right). “Unknown” cells had no mutant reads or no coverage. “WT” cells were from TFHL patients without G17V mutations by WES. MUT, mutant cells. B Distribution of tumor cells with chromosome (chr)5 gain (left panel) and comparison of GSVA enrichment score for the TFH signature between tumor cells with chr5 gain (pink) and those without (blue) (right panel). In the violin plots of the right panel, adjusted P values are calculated for all tumor cells (left) and each tissue (right). NS, not significant. C Phylogenetic trees based on copy number variation (CNV) patterns identified by inferCNV [35] in TFHL LN (left panel) and PB (right panel) samples with partial chr5 gain. All boxplots show the median (center line), mean (center dot), interquartile range (box limits), and minimum to max values (whiskers) for each group. *P < 5.0 × 10−2, ****P < 1.0 × 10−4.

Next, we conducted single-cell level estimation of CNVs in LN tumor cells and classified them into subgroups based on their CNV patterns using inferCNV [35]. The estimated CNVs exhibited a high degree of consistency with paired bulk WES results (Fig. S4B and S7A). Similar to WES, the most frequently observed CNV in LN tumor cells was chr5 gain (6/9; 66.7%), comprising four samples detected by both scRNA-seq and WES, and two samples exclusively identified through scRNA-seq (AITL10LN and AITL12LN; Fig. S7A). In these two samples and AITL8LN, chr5 gain was identified only in a portion of tumor cells with the identical TCR clones, suggesting that it was acquired after TCR clonal expansion (Fig. S7A). Notably, chr5 gain was most prevalent in C3–4, showing enriched TFH phenotypes (C3–4) and active cell proliferation (C4), followed by C0 (Fig. 3B). Consistently, DEG analysis revealed that TFH-related genes, including CXCR5, CXCL13, IL21, and CD40LG, were significantly upregulated in chr5-gain tumor cells compared to those without and GSVA showed that the TFH signature enrichment was observed in total and LN tumor cells while being insignificant in PB tumor cells (Fig. 3B, Fig. S7B). Furthermore, the phylogenetic trees of partial chr5-gain samples diverged into two major branches based on the presence or absence of chr5 gain, plus the subsets with chr5 gain accompanied by other CNVs frequently detected in TFHL [28] (e.g., gains of chr7, 19, and 21) that evolved into C3–4 (Fig. 3C). In these samples, the TFH signature was significantly enriched in chr5-gain positive cells compared to non-gain (Fig. S7C). Contrastingly, in AITL10PB and AITL12PB, where partial chr5 gain was detected in LN tumor cells, fewer chr5-gain cells were detected in PB tumor cells compared to those of LNs (41.1 versus 7.8% and 11.0 versus 1.1%), and most of the tumor cells belonged to C0–2, suggesting that C0–2 without chr5 gain and TFH signatures favored infiltration into PB versus C3–4 (Fig. 3C, Fig. S7D). Indeed, CNV scores [36] were significantly higher in C3–4 and LN tumor cells than those in C0–2 and PB tumor cells, respectively, indicating the genomic complexity of C3–4 and LN tumor cells (Fig. S7E).

In summary, genomic aberrations, particularly G17V and chr5 gain, were closely correlated with distinct transcriptomic profiles in TFHL tumor cells at the single-cell resolution.

PLS3 as a novel, tumor-specific marker

To discover novel tumor-specific markers, we conducted DEG analysis using scRNA-seq data between LN/PB tumor cells and non-malignant MNCs or normal TFH cells (for LNs only). DEGs with significantly higher expression in tumor cells across ≥2 samples, as well as those upregulated in both LN and PB tumor cells, were extracted (Fig. 4A; Table S8). Finally, PLS3, IGFL2, and CA8 were identified as candidate genes (Fig. S8A; Table S9). FCM was performed to validate PLS3, revealing high expression in both LN and PB tumor cells versus non-malignant T cells (Fig. 4B). Furthermore, immunohistochemical staining revealed that PLS3 was strongly expressed on the tumor cell-surface in 17 of 35 PTCLs (10 of 20 angioimmunoblastic T-cell lymphomas, 2 of 5 nodal PTCLs with TFH phenotypes, and 5 of 10 PTCLs not otherwise specified), accounting for 48.6% PTCL samples examined, while significant expression was observed in only one of 168 (0.6%) mature B-cell lymphomas (Fig. 4C; Table S10). These findings suggest that PLS3 is a promising, tumor-specific marker in PTCL.

Fig. 4: Identification of novel tumor-specific cell markers.

A Overview of tumor-specific cell marker identification. For both LN (n = 9) and tumor cell-bearing PB (n = 12) samples, differentially expressed gene (DEG) analysis was performed to select genes whose expression was significantly elevated in tumor cells of two or more samples versus (vs.) all MNCs or normal TFH cells (for only LN tumor cells) for all samples. Next, genes specifically expressed in tumor cells were selected by manual inspection using UMAP plots. B Representative FCM plots of PLS3 expression in PD1bright CD4+ cells from PB of TFHL (AITL4, red) and PD1dim CD4+ cells from PB of a healthy donor (HD, blue). C Representative images of PLS3 expression by immunohistochemical staining of FFPE samples from TFHL (AITL4, left) and B-cell lymphoma (B03, right). White triangles indicate tumor cells expressing PLS3 on the cell membrane surface. Images were scanned using the NanoZoomer (2.0HT, Hamamatsu Photonics, Shizuoka, Japan) and acquired at ×10 and ×40 with the NDP.view.2 software (v2.9.29). FL, follicular lymphoma; G, grade of FL according to the 4th edition of the World Health Organization classification [1]; Scale bar, 250 μm (left) and 50 μm (right).

Dysfunctional CD8+ and regulatory T-cell expansion

Subsequently, subclustering of non-malignant LN T cells highlighted 11 clusters (T0–10). Based on canonical markers, these clusters were classified as conventional CD4+ or CD8+ T, encompassing naïve (TN, T0 and T6), central memory (TCM, T1 and T7), TEFF (T2 and T8), regulatory T (TREG, T3), TFH (T4), proliferating CD4+ T (CD4 TPRO, T5), dysfunctional CD8+ T (CD8 TDYS, T9), and proliferating dysfunctional CD8+ T (CD8 TPRO/DYS, T10) (Fig. 5A, Fig. S9A–E). T9 and T10 expressed “exhaustion” or “dysfunctional” markers such as PDCD1, HAVCR2 (encoding TIM-3), TIGIT, and LAG3 [37, 38] (Fig. S9C, D). T5 was subclassified into two FOXP3+IL2RA+ clusters (T5 FOXP3-1 and T5 FOXP3-2) and one cluster positive for the TFH marker (T5 PDCD1) (Fig. S9F), suggesting that T5 was a mixture of proliferating TREG and TFH. These manual annotations were further validated using SingleR [39, 40] (Fig. S9G).

Fig. 5: Subclustering of non-malignant T cells from LNs.

A UMAP plots of non-malignant LN T cell subclusters. Cells are shown separately for each clinical status (right panels). TCM, central memory T cell; TDYS, dysfunctional T cell; TEFF, effector T cell; TN, naïve T cell; TPRO, proliferative T cell; TPRO/DYS, proliferative dysfunctional T cell; TREG, regulatory T cell. B Comparison of proportions of each cluster in non-malignant MNCs per sample. The boxplots show the median (center line), interquartile range (box limits), minimum to max values (whiskers), and samples (dots) for each group. P values are shown only for significant differences. *P < 5.0 × 10−2. C Clone size (i) and overlap analysis (ii) of TCRs in non-malignant LN T cells. The number of cells expressing each clonotype was defined as clone size and illustrated for each cell. The “repOverlap” function of Immunarch measured the TCR sharing between each cluster, analyzed and illustrated for TFHL LNs (upper right) and HLNs (lower left), respectively.

The proportions of CD4 TPRO and CD8 TDYS/TPRO/DYS significantly increased in TFHL LNs versus HLNs, whereas those of CD4/CD8 TN and CD4 TCM decreased as previously reported [41] (Fig. 5B). Additionally, the proportion of TREG was higher in RR TFHL LNs than HLNs (Fig. 5B, Fig. S9H). Moreover, compared with HLNs, TREG in RR TFHL LNs exhibited higher expressions of genes associated with TREG inhibitory functions, such as FOXP3, IL2RA, TNFRSF4, TNFRSF9, and LAG3 [42], suggesting TREG activation in RR TFHL (Fig. S9I; Table S11). Additionally, non-malignant PB T cell subclustering demonstrated the expansion of CD8 TDYS/TPRO/DYS and upregulation of activated markers in TREG from TFHL (Fig. S10A–F; Supplementary Note 1).

We next analyzed the TCR clonality of non-malignant LN T cells and identified oligo-clonal expansion of TCRs in CD8 TEFF, TDYS, and TPRO/DYS of TFHL (Fig. 5C(i), Fig. S10G). Indicative of a common cellular origin, these cells partly shared identical TCRs between clusters, an overlap not observed in HLNs (Fig. 5C(ii)). Moreover, trajectory analysis revealed that a portion of CD8 TEFF diverged and differentiated into CD8 TDYS and TPRO/DYS in TFHL, consistent with TCR analysis (Fig. S10H). Notably, TCR commonality was observed between LN and PB CD8 TPRO/DYS of TFHL, suggesting that enrichment of immunoevasive phenotypes in PB potentially reflects the LN immune environment (Fig. S10I).

Increased immunosuppressive myeloid cells

Subclustering LN myeloid cells distinguished seven clusters (M0–6; Fig. 6A). Employing canonical markers and SingleR analysis [43, 44], we annotated each cluster as classical and intermediate monocytes (Class Mono, M0; Inter Mono, M1, respectively), complement component 1Q (C1Q)-positive macrophages (C1Q+ Mφ, M2), CD1C-positive type 2 conventional dendritic cells (CD1C+ cDC2, M3), CLEC9A-positive type 1 cDCs (CLEC9A+ cDC1, M4), plasmacytoid DCs (pDC, M5), and LAMP3-positive cDCs (LAMP3+ cDC, M6) (Fig. 6A, Fig. S11A–D).

Fig. 6: Subclustering of myeloid cells from LNs.

A UMAP plots of LN myeloid cell subclusters. Cells are shown separately for each clinical status (right panels). C1Q+ Mφ, complement component C1q positive macrophage; CD1C+ cDC2, CD1C-positive type2 conventional dendritic cell; Class Mono, classical monocyte; CLEC9A+ cDC1, CLEC9A-positive type1 cDC; Inter Mono, intermediate monocyte; LAMP3+ cDC, LAMP3-positive cDC; pDC, plasmacytoid DC. B Trajectory inference by Monocle2 for cDCs of TFHL, color-coded by cluster (left of the left panel) and pseudo-time (right of the left panel). Cells are shown separately for each cluster (right panel). C Comparison of proportions of each cluster in non-malignant MNCs per sample. The boxplots show the median (center line), interquartile range (box limits), minimum to max values (whiskers), and samples (dots) for each group. P values are shown only when there is a significant difference. *P  < 5.0 × 10−2.

Notably, C1Q+ Mφ exhibited high C1Q-family gene expression (C1QA, C1QB, and C1QC), CD163 expression, and ligands for T-cell immune checkpoints, including LGALS9 (binds to TIM-3) and NECTIN2 (binds to TIGIT) (Fig. S11C), resembling tumor-associated macrophages with immunosuppressive signatures reported in multiple cancers [43, 45]. Trajectory inference showed a continuous spectrum from Class/Inter Mono to C1Q+ Mφ (Fig. S11E). LAMP3+ cDCs were characterized by high immune regulatory gene expression, including CD274 (encoding PD-L1) and PDCD1LG2 (encoding PD-L2) [46], consistent with “mregDCs,” a subset of DCs co-expressing maturation and immunoregulatory genes identified in various cancers [43, 46] (Fig. S11B, C and F). LAMP3+ cDCs have been proposed to develop from either cDC1s or cDC2s [43, 46]. Using a previously reported scoring system [43, 46], we inferred that most LAMP3+ cDCs in TFHL originated from cDC2s, unlike many cancers [43, 46] (Fig. S11G). Trajectory analyses further supported this, revealing a continuous connection from cDC2 to LAMP3+ cDCs (Fig. 6B). Next, we compared the proportions of each cluster and observed a significant increase in Inter Mono, C1Q+ Mφ, and LAMP3+ cDCs in RR TFHL LNs than in HLNs (Fig. 6C). Finally, DEG analysis showed the expression level of CCL17, a ligand for CCR4 [47], was significantly higher in LAMP3+ cDCs of RR TFHL LNs than HLNs (Fig. S11H, I).

Taken together, a distinct increase in immunoregulatory myeloid cells, particularly LAMP3+ cDCs, was observed in RR TFHL, which could support an immune-evasive TME.

Exhausted phenotype and clonal expansion of B cells

Subclustering analysis of LN B cells revealed 11 clusters: naïve B cell (NBC, B0), memory B cell (MBC, B1), FCRL4-positive MBC (FCRL4+ MBC, B2), pre-germinal center B cell (preGCB, B3), light-zone GCB (GCB [LZ], B4), dark-zone GCB (GCB [DZ], B5), plasmablast (PBL, B6), and plasma cell (PC, B7–10) (Fig. 7A, Fig. S12A–C). FCRL4+ MBCs displayed “exhausted” or “atypical” MBC phenotypes, as observed in chronic inflammatory or immunodeficiency conditions [48,49,50], distinguished by high expression of TBX21 (encoding T-bet), B cell inhibitory receptor genes (FCRL4 and SIGLEC6 [49]), and homing receptor genes (CXCR3 and ITGAX [48, 49]) as well as low expression of CR2 (encoding CD21) and CD27 (Fig. S12D; Table S12). Compared to HLNs, the proportions of each cluster varied between TFHL patients (Fig. S12A). While statistically insignificant, the proportion of FCRL4+ MBCs was markedly increased (>20% of B cells) in 2 (AITL9LN and AITL12LN) out of 3 RR TFHLs (Fig. S12A). Moreover, DEG analysis highlighted the greater upregulation of genes related to exhausted MBC phenotypes [48,