Single cell transcriptome data containing 2 reactive lymph node tissue samples and 2 DLBCL lymph node biopsy tissue samples were obtained from the heiDATA database (https://heidata.uni-heidelberg.de) (Supplementary Table S1). Single cell samples were prepared and Single-cell RNA sequencing (scRNA-seq) as follows: single cell suspensions, synthetic complementary DNA and single cell libraries were prepared using Chromium Single Cell v2 3ʹ kits (10x Genomics) according to the manufacturer’s instructions. Each was sequenced on a single NextSeq 550 lane (Illumina). The data were aligned to the hg38 reference genome with Cell Ranger (v2.1, 10x Genomics) using “mkfastq” and “count” commands and default parameters. The results of the Cell Ranger analysis contained the count values of unique molecular identifiers assigned to each gene in each of the cells for each individual sample using all mapped reads (Supplementary Table S2).

The R package Seurat (v4.0.2) (Butler et al., 2018) was used to perform quality control. Gene counts per cell, UMI counts per cell, and percentages of mitochondrial and ribosomal transcripts were calculated using the functions of the Seurat package. Genes expressed in three or fewer cells were excluded from downstream analysis. Before further analysis, libraries with >5% of mitochondrial transcripts, libraries with UMI numbers indicating an abnormal range of potential doublets, and libraries with less than 200 genes were screened out. After removing low-quality cells, we analyzed scRNA-seq profiles of 11,729 cells with an average sequencing depth of approximately 1,400 genes per cell.

The canonical correlation analysis (CCA) and mutual nearest neighbor (MNN) algorithms in the R package Seurat (v4.0.2) (Butler et al., 2018) were used for sample whole and correction of batch effects. After identifying the different cell types, the subsetdata function was used to split the dataset into subsets of different cell types.

We used Seurat (v4.0.2) (Butler et al., 2018) to perform clustering analysis of cells. Data was normalized to log scale using the “NormalizeData” function with a default scale parameter of 10,000. “FindVariableFeatures” function was used to identify highly variable genes with parameters for “selection.method = vst, nfeatures = 2000”. We standardized the data with the “ScaleData” function. These variable genes were used as input for PCA using the “RunPCA” function. The first 20 principal components (PCs) and a resolution of 0.5 were used for clustering using “FindClusters”. Uniform manifold approximation and projection for dimension reduction (UMAP) was used for two-dimensional representation of first 20 PCs with “RunUMAP”. We used the “FindAllMarkers” or “FindMarkers” function to determine the marker genes of each cluster relative to all other clusters or to a specific cluster. The selected parameters of marker genes were detected in at least 25% of the cells in the target cluster, under p value of Wilcoxon test <0.05 and the differential expression threshold of 0.25 log fold change. FeaturePlot, DotPlot, VlnPlot and DoHeatmap were used for visualization of gene expression levels. We labeled the obtained clusters as T cells, B cells, NK cells, Dendritic cells (DC) and monocytes by known classical markers (T cells: CD3D, CD3E, CD3G, TRAC; B cells: MS4A1, CD79A; NK cells: NKG7, GNLY; DC: IRF7, IRF8; monocytes: LYZ, CD68.).

Because DLBCL1 contains significantly more cells than DLBCL2, in order to perform a systematic analysis of intercellular communication, we re-clustered DLBCL1 for annotation and used the R package CellChat (v1.1.3) (Jin et al., 2021) to explore the expression of ligand-receptor pairs.

Branching developmental trajectories of CD8⁺ T cell subpopulations were calculated using the R package Monocle 2 (v2.16.0) (Qiu et al., 2017). Monocle introduces the strategy of ordering single cells in pseudo-time, by taking advantage of the asynchronous progression of individual cells in these processes and aligning them along trajectories corresponding to biological processes, such as cell differentiation.

After annotation of each cell type by characterization of cell type marker genes, we used the SCENIC package (v1.2.4) (Aibar et al., 2017) to analyze the enriched transcription factors in cell subpopulations. The input matrix is a normalized expression matrix, output by Seurat.

Hallmark gene sets were downloaded from the MSigdb (Molecular Signatures Database) database and Gene Set Variation Analysis (GSVA) was performed using the R package GSVA to determine the molecular characteristics of different cell subpopulations. Gene-cell matrices are converted into gene set-cell matrices and GSVA scores are calculated for sets with at least 5 detected genes; all other parameters are default.

RNA-seq data and clinical information of 481 DLBCL patients were downloaded from The Cancer Genome Atlas (TCGA) database (https://cancergenome.nih.gov/) for screening prognostic genes and developing prognostic models. RNA sequencing data and clinical information for 420 DLBCL patients from the external validation cohort GSE10846 dataset were obtained from the GEO database. Data for the IMvigor210 immunotherapy cohort were obtained from the website http://research-pub.gene. com/IMvigor210CoreBiology. Extracted CD8⁺ T_EXH subpopulation-related genes obtained from differential gene expression analysis were used to construct prognostic models. In the TCGA cohort, univariate Cox regression analysis was performed using the R package Survival to screen prognosis-related genes (p < 0.05). Lasso regression analysis was performed using the R package glment to further screen prognosis-related genes, and finally six prognosis-related genes were obtained by multivariate Cox regression analysis for the construction of the prognostic risk model. The risk score of each patient was calculated as follows: Risk score = ∑ j = 1 n ( βj   ×   expGj ) where β is the regression coefficient obtained by multivariate Cox regression analysis and expG is the prognostic gene expression level. Based on the median risk scores obtained from the prognostic model, the DLBCL samples were divided into high-risk and low-risk groups, and survival differences between the different risk subgroups were compared by Kaplan-Meier curves. We plotted time-dependent subject operating characteristic (ROC) curves with 1, 3 and 5 years as the defined points, calculated the corresponding area under the ROC curve to assess the predictive power of the risk model, and verified whether the risk score was an independent prognostic indicator for DLBCL by Cox regression analysis. The GSE10846 cohort was used as an independent external validation cohort to verify the efficacy of the prognostic model.

ESTIMATE is an algorithm that uses expression data to estimate stromal and immune cells in malignant tumor tissues, allowing estimation of stromal and immune scores for each DLBCL sample (Yoshihara et al., 2013). The deconvolution algorithm CIBERSORT is a method for characterizing cell composition from gene expression profiles of complex tissues, allowing inference of the relative content of immune cells from large amounts of tumor transcriptome data (Newman et al., 2015). Gene set enrichment analysis was performed using GSEA software (v4.1.0) to identify pathways that are predominantly enriched between high- and low-risk groups. Significantly enriched gene sets were screened with a threshold of p < 0.05. To validate the predictive power of prognostic models for immunotherapy response, the IMvigor210 immunotherapy cohort was used to assess differences in response to PD-L1 treatment in patients in different risk groups. Spearman correlation analysis was used to characterize the correlation between immune checkpoint genes and risk scores.

All statistical analyses were performed in R (v4.0.5). Comparisons between groups were performed using the Wilcoxon test and t-test. Correlations were analyzed by using Spearmans correlation. Survival curves were compared using log-rank test. Statistical significance was accepted for p < 0.05. *p < 0.05, **p < 0.01, ***p < 0.001.

In this study, single-cell transcriptomic data obtained from 10x Genomics sequencing were used to investigate the cellular diversity and molecular features in DLBCL tissues. After data quality control and filtering, 11,729 cells were obtained for subsequent analysis. After normalization of gene expression data, descending and clustering were performed using principal component analysis and UMAP, respectively. Twelve cell subpopulations were obtained by dimensionality reduction and clustering (Figure 1A), and these cells were assigned to five different cell types using known marker genes (Figures 1B,D): B cells (marker genes: MS4A1 and CD79A), T cells (marker genes: CD3D, CD3E, CD3G and TRAC), NK cells (marker genes: GNLY and NKG7), DC cells (marker genes: IR7 and IR8), monocytes (marker genes: LYZ and CD68). Notably, B cells and T cells are the major cell subsets of DLBCL (Figure 1C).

To investigate the transcriptomic heterogeneity of malignant B cells in DLBCL tissues, we re-clustered the B cells and identified 13 cell subpopulations. (Figure 2A). To further distinguish malignant B cells from non-malignant B cells, we took advantage of the fact that the malignant B cell population expresses only one type of immunoglobulin light chain, i.e. κ or λ light chains. The ratio of light chains per B cell (κ/λ) was calculated based on the expression of the genes IGKC (encoding a constant portion of the κ light chain) and IGLC2 (λ light chain). Malignant lymph nodes contain malignant B cells that uniformly express κ light chains, whereas reactive lymph node samples contain only non-malignant B cells (Figure 2B). We then re-clustered the malignant B cells and obtained eight malignant B cell subpopulations (Figure 2C), which showed a high degree of heterogeneity. SCENIC analysis identified EGR1, FOS and STAT1 as potential transcription factors (Figure 2D). Gene differential expression analysis revealed different transcriptional profiles among malignant B cell subpopulations: subpopulation 0 showed high expression levels of the malignancy-promoting factors S100A6 and LY6E, subpopulation 1 showed high expression levels of the tumor suppressor BTG1 and TXNIP, subpopulation 2 showed high expression levels of the immune-related genes CD74 and HLA-DRA, subpopulation 3 and subpopulation 4 showed high expression levels of cell proliferation genes MCM3, H2AFY, PCN, MKI67, TK1, subpopulation 5 showed high expression levels of metabolism-related genes FABP5, LDHA, ENO1, and subpopulation 6 showed high expression levels of cell cycle-related genes CENPF, CCNB1, CDC20 (Figure 2E). GSVA analysis showed different molecular signatures among malignant B-cell subpopulations: interferon response-dominant signature (subpopulation 0), cell proliferation-dominant signature (subpopulation 3 and subpopulation 4), metabolism-dominant signature (subpopulation 5), and hypoxia-dominant signature (subpopulation 7) (Figure 2F). In conclusion, these results reveal a high degree of inter-tumor heterogeneity in DLBCL.

Tumor-infiltrating immune cells are highly heterogeneous and play an important role in tumor cell immune evasion and response to immunotherapy. To investigate the transcriptomic heterogeneity of T cells in DLBCL tissues, we re-clustered T cells and identified 13 T cell subpopulations (Figures 3A–D). The T cell subpopulations were annotated by differentially expressed marker genes as: CD4⁻CD8⁻Navie T (IL7R,SELL, CCR7 and LEF1, subpopulations: 0, 3, 4 and 10), CD4⁺ T_H (CD4 and TRAC, subpopulations: 1, 2 and 5), CD8⁺ T_TOX (CD8A, GZMK and NKG7, subpopulation: 6), CD4⁺ T_REG (FOXP3,TIGIT, ICOS and CTLA4, subpopulations: 7 and 9), CD8⁺Navie T (CD8A, SELL and IL7R, subpopulation: 8), T_PRO (MKI67 and TOP2A, subpopulation: 11), CD8⁺ T_EXH (CD8A, GZMA, NKG7, LAG3 and HAVCR2, subpopulation: 12). To understand the state transitions between CD8⁺ T cell subtypes, we used Monocle2 to construct potential developmental trajectories of T cells. Developmental trajectories inferred from expression data or marker genes suggest (Supplementary Figures S1A, B) that CD8⁺ T cells have two differentiation pathways: cytotoxic CD8⁺ T cells (CD8⁺ T_TOX) and exhausted CD8⁺ T cells (CD8⁺ T_EXH). GSVA analysis revealed different signaling pathway enrichment among subpopulations: WNT and TGF signaling (CD4⁺ T_H), TGF and TNF signaling (CD8⁺ T_TOX), IL6/STAT3, IL2/STAT5 and KRAS signaling (CD4⁺ T_REG), and interferon response (CD8⁺ T_EXH) (Supplementary Figure S1C). SCENIC analysis identified SREBF2, RAD21, IRF7 as potential transcription factors in different T cell subpopulations (Supplementary Figure S1D). Taken together, our single-cell analyses reveal that CD4⁺ T_REG are highly immunosuppressive and CD8⁺ T_EXH highly express exhaustion markers such as LAG3, TIGIT and HAVCR2.

To explore the interactions between cells in the DLBCL microenvironment, we used CellChat to infer and analyze intercellular communication networks. Dimension reduction, clustering and cell type annotation of sample DLBCL1 identified 13 cell subpopulations containing 9 malignant B cell subpopulations (MB1-9), 3 T cell subpopulations (T_REG, T_TOX, Naive T), and 1 DC cell subpopulation (DC). CellChat analysis revealed complex interactions between malignant B cell subpopulations and with other cell subpopulations, and 22 important pathways between 13 cell subpopulations were detected in DLBCL tissues, with the MIF signaling pathway being the prominent incoming and outgoing signaling mode (Figures 4A–D). Network centrality analysis of the inferred MIF signaling network showed that malignant B cell subpopulations (MB-2, MB-7) are the major senders and DCs are the major receivers of the MIF signaling pathway (Figure 4E). Notably, among all known ligand-receptor pairs, MIF signaling was predominantly dominated by the MIF ligand and its multimeric CD74/CXCR4 receptor (Figure 4F). CellChat uses a pattern recognition approach based on non-negative matrix decomposition to identify global communication patterns as well as key signals in different cell groups (i.e. pattern recognition modules). The output of this analysis is a set of the so-called communication patterns that connect cell groups with signaling pathways either in the context of outgoing signaling (i.e. treating cells as sources) or incoming signaling (i.e. treating cells as targets). The application of this pattern recognition module revealed three patterns of the outgoing signal and three patterns of the incoming signal (Figures 5A,B). The outgoing signaling of all malignant B cells is characterized by pattern #1, which includes the MHC-II, MIF, MHC-I, CD22, CD45 and other pathways, the outgoing signaling of T cells is characterized by pattern #2, which represents the ADGRE5, LCK, IFN-II, VCAM, PECAM1 and other pathways, and the outgoing signaling of DC is characterized by pattern #3, which includes the APP, BAFF, ICAM and other pathways. On the other hand, the communication patterns of target cells show that incoming malignant B cell signaling is dominated by patterns #1, which includes signaling pathways such as CD22, CD45, CD70, BAFF, IFN-II, etc. Incoming T cell signaling is characterized by two patterns #2 and #3, driven by pathways such as MHC-I, LCK, VCAM, ICAM, etc., while incoming DC signaling is also characterized by patterns #3. These results suggest that different cell types in the same tissue have different signaling networks and the pattern of malignant B-cell communication is homogeneous.

The CD8⁺ T_EXH subpopulation-related genes obtained from the differential analysis were extracted for the construction of the prognostic model. Nineteen genes were obtained by univariate Cox regression analysis and lasso regression analysis (Figures 6A,B), and finally six prognosis-related genes for model construction were obtained using multivariate Cox regression analysis (GABRA3, HOXC8, RTN4R,CRLF1, BIRC3, REXO5). Using the regression coefficients for each of the above 6 prognostic genes, we constructed a prognostic model for DLBCL patients and calculated the risk score according to the following formula: risk score = (2.201 × GABRA3 expression level) + (−0.719 × HOXC8 expression level) + (−0.765 × RTN4R expression level) + (0.545 × CRLF1 expression level) + (−0.013 × BIRC3 expression level) + (−0.226 × REXO5 expression level). Using the median value of the risk score as the threshold, we divided DLBCL patients into low-risk and high-risk groups. Survival analysis showed that patients in the high-risk group had a poorer prognosis (p < 0.001) (Figure 6C), with an area under the ROC curve of 0.83, 0.80 and 0.80 for 1-year, 3-years and 5-years OS, respectively (Figure 6D). In the external validation cohort, survival analysis also showed a poorer prognosis for patients in the high-risk group (p < 0.001) (Figure 6E), with an area under the ROC curve of 0.71, 0.70 and 0.63 for 1-year, 3-years and 5-years OS, respectively (Figure 6F). The results of univariate and multivariate Cox regression analyses indicated that risk score was an independent prognostic factor (Figures 6G,H).

We performed gene set enrichment analysis (GSEA) to identify potential biological processes between high- and low-risk groups. The results showed that pathways such as nitrogen metabolism, oxidative phosphorylation, ribosomes, Alzheimer’s disease, and Parkinson’s disease were enriched in the high-risk group, and pathways such as extracellular matrix receptor interactions, focal adhesion, gap linkage, pathways in cancer, and regulation of the actin cytoskeleton were enriched in the low-risk group (Figures 7A,B).

We used the ESTIMATE algorithm to assess the TME immune and stromal abundance in the different risk groups, and the results showed that the high-risk group had higher levels of immune and stromal component abundance (Supplementary Figures S2A–C). We also analyzed the proportion of 22 types of immune infiltrating cells among different risk groups in 481 DLBCL samples using the CIBERSORT algorithm (Supplementary Figure S2D), and the results showed that seven types of immune infiltrating cells were associated with risk scores: resting CD4 memory T cells, activated CD4 memory T cells, regulatory T cells, γδ T cells, and M0, M1, and M2 macrophages (Figure 8A). Correlation analysis of risk scores with immune checkpoint genes showed that most of the immune checkpoint gene expression levels were positively correlated with risk scores (Figure 8B). In addition, higher risk scores in the IMvigor210 immunotherapy cohort were associated with anti-PD-L1 treatment response (Figure 8C).