A flow chart is outlined in Figure 1 to summarize the patient selection. From September 2002 to April 2021, a total of 2,445 patients with newly diagnosed DLBCL were included, with the last follow-up through August 31, 2021. Histological diagnosis was established based on the World Health Organization (WHO) classification reviewed by two experienced pathologists, LD and HY (12). Survival analysis was performed on 1,896 patients receiving R-CHOP-based immunochemotherapy. DNA and RNA sequencing were performed on 1,150 and 385 patients with available tumor and blood samples, respectively, for detection of genetic aberrations, tumor microenvironmental analysis, and gene set enrichment analysis (GSEA). This study was approved by the Shanghai Rui Jin Hospital Review Board with informed consent obtained in accordance with the Declaration of Helsinki.

DNA sequencing of 340 patients by whole-genome sequencing (WGS, n = 117) and whole-exome sequencing (WES, n = 223) were performed on available frozen or quality controlled FFPE tumor samples, as reported by our previous studies, along with detailed procedures for DNA sequencing (13–15). Targeted sequencing of 55 lymphoma associated genes were performed among 810 patients with available FFPE tumor samples based on the criteria as previously described (13–15). The mean depth of samples sequenced by WES and WGS was 120.25× (range 50–200×), with an average 97.65% (range 82.64–99.06%) of the target sequence being covered sufficiently deep for variant calling (≥10× coverage). The mean depth of samples sequenced by targeted sequencing was 1,351× (range 505–3,224×), with an average 88.57% (range 59.56–98.26%) of the target sequence being covered sufficiently deep for variant calling (200× coverage). Single nucleotide variations (SNVs) and indels were called by Genome Analysis Toolkit (GATK, v3.4) Haplotype Caller, and GATK Unified Genotyper and mapped to the genome location using the UCSC Genome Browser (http://genome.ucsc.edu) for annotation. The filtration of SNVs and indels was carried out by homemade pipeline with the software mentioned above. Variant allele frequency of mutations included should be over 5%. Mutations were filtered according to the rules listed below. Mutations were preserved if they met the following conditions: 1) sites reported as somatic mutations in our previous studies (13–15); 2) sites verified as somatic mutations by sequencing on paired blood samples; 3) sites commonly considered as hotspot mutations like MYD88L265P; 4) sites categorized into tier I and II variants according to the Guideline for Evidence-Based Categorization of Somatic Variants (16); 5) sites not included in the SNP database, or related to hematological malignancies (N >5) as reported in the COSMIC (the Catalogue of Somatic Mutations in Cancer). Mutations were excluded if they met the following conditions: 1) sites verified as germline mutations by sequencing on paired blood samples; 2) sites categorized into tier IV variants according to the Guideline for Evidence-Based Categorization of Somatic Variants; and 3) sites with extraordinarily high frequencies but never reported in previous studies. The mutation data among the 55 genes are listed in Supplementary Table 1 . Clinical and pathological features of patients with WGS/WES/targeted sequencing data are shown in Supplementary Table 2 . No statistically significant difference was found in the frequencies of mutations detected by WGS/WES/targeted sequencing called in the components of the cohort after balancing the baseline of clinical and pathological characteristics using propensity score matching, except DDX3X ( Supplementary Tables 3 , 4 ).

Based on the Gene Ontology database, 55 mutated genes were assigned to biological processes, namely, chromatin organization, immune response, cell-cycle/p53, and also oncogenic signaling pathways B-cell receptor (BCR)/NF-κB, JAK-STAT, PI3K-AKT, and Wnt. The list of genes involved in biological function is shown in Supplementary Table 5 .

RNA sequencing was performed on 385 patients with available frozen tumor samples. A total of 361 patients of them were from our previous studies (13–15) and 24 patients were newly analyzed, namely, 98 patients with WGS, 98 with WES, and 189 with targeted sequencing data, respectively. Clinical and pathological features of patients with RNA sequencing data are shown in Supplementary Tables 6 – 8 . RNA was extracted by Trizol and an RNeasy Mini Kit (Qiagen) from available qualified frozen tumor samples of 385 patients. RNA was purified by Ribo-Zero rRNA Removal Kits (Illumina). RNA concentration was assessed on NanoDrop and integrity by an Agilent 2100 Bioanalyzer. RNA libraries were generated with a TruSeq RNA Library Preparation Kit (Illumina) based on instructions of the manufacturer. The concentration and the quality of RNA libraries were controlled by Qubit and BioAnalyzer 2100 system. Paired-end sequencing was performed on Illumina HiSeq sequencer following Illumina-provided protocols. The read pairs were aligned to Refseq hg19 through Burrows‐Wheeler Aligner version 0.7.13‐r1126. The transcript counts table files were generated by the HTSeq (17). Potential false positive results were excluded by visual inspection. R package “sva” was applied to remove batch effect by r 4.0.3. The expressional data related are shown in Supplementary Table 9 . Cell of origin classification based on gene expression was performed using Lymph2Cx assay (18). Based on raw counts data of RNA sequencing, the infiltration of immune cells including 24 immune cell types were estimated by ImmuCellAI algorithm, a website tool (http://bioinfo.life.hust.edu.cn/web/ImmuCellAI/) (19). Pathway enrichment analysis was performed using GSEA v4.0.1 software as recommended by the GSEA team (http://www.broadinstitute.org/gsea). Pathways were considered of statistical significance with a P-value <0.05 and a false discovery rate <0.25.

Immunohistochemistry was performed on paraffin sections using antibodies against CD10, BCL6, MUM1, BCL2, and MYC by indirect immunoperoxidase method. Germinal center B-cell-like (GCB) and non-GCB phenotypes were determined with 30% cutoff values of CD10, BCL6, and MUM1, according to Han’s algorithm (20). BCL2/MYC double-expressors were defined by BCL-2 and MYC with cutoff values of 50 and 40%, respectively (12).

Fluorescence in-situ hybridization of BCL2, BCL6, and MYC rearrangements was performed on paraffin sections with 10% cutoff values. The results in detail are listed in Supplementary Table 10 .

DLBCL genotypes were identified as described by Lacy et al. using the 47 genes available among 1,150 patients with DNA sequencing data (R code version, https://github.com/ecsg-uoy/DLBCLGenomicSubtyping) (21).

Pearson’s χ ² test or Fisher’s exact test was used to analyze baseline characteristics of patients. Progression-free survival (PFS) was calculated from the date of diagnosis to the date when disease progression or relapse was recognized or the date of last follow-up. Overall survival (OS) was measured from the date of diagnosis to the date of death or the date of last follow-up. Survival functions were analyzed using the Kaplan–Meier method and compared by the log-rank test. Univariate hazard was analyzed using the Cox regression method and the significant variables were then kept in multivariate set. Normalized gene expression in two groups was analyzed using Mann–Whitney U test. All statistical analysis was performed by Statistical Package for the Social Sciences (SPSS) 26.0 software. All tests were two-sided, with statistical significance defined as P <0.05.

Among 2,445 patients with newly diagnosed DLBCL, 1,140 patients were >60 years old and 1,305 patients were ≤60 years old at diagnosis. The clinical characteristics of the patients are listed in Table 1 . Patients diagnosed at age >60 years were associated with advanced Ann Arbor stage (P <0.001), elevated serum lactate dehydrogenase (LDH) (P <0.001), poor performance status (P <0.001), multiple extranodal involvement (P = 0.001), high percentage of double expressor subtype (P = 0.012), as compared to those at age ≤60 years. No significant difference was observed between the two age groups, according to gender, cell of origin (Hans), cell of origin (Lymph2Cx), or double/triple-hit.

Among 1,896 patients receiving R-CHOP-based immunochemotherapy, the median follow-up time was 55.0 months (0.2–224.2 months). The 3-year PFS and OS rates of patients diagnosed at age >60 years were 60.8 and 71.0%, significantly lower than those of patients at age ≤60 years (72.1%, P <0.001; 83.2%, P <0.001, respectively) ( Figures 2A, B ). Using univariate analysis, the 5 factors of IPI were of great significance and were included in the multivariate analysis. In the Cox proportional-hazards model, independent prognostic factors of inferior PFS and OS were age (P both <0.001), along with Ann Arbor stage (P both <0.001), serum LDH (P both <0.001), performance status (P both <0.001), and multiple extranodal involvement (P both <0.05) ( Table 2 ). Similar results were observed based on age factor according to NCCN-IPI. The 3-year PFS and OS rates of patients >75 years were 51.2 and 57.6%, significantly lower than those of patients 61–75 years (62.0 and 72.8%), 41–60 years (72.5 and 83.4%), and ≤40 years (71.4 and 82.8%) ( Figures 2C , D ).

Oncogenic mutations closely related to age at diagnosis were analyzed in 1,150 patients, namely, 117 cases by WGS, 223 cases by WES, and 810 cases by targeted sequencing. A total of 55 genes related to the tumorigenesis of DLBCL were analyzed ( Figure 3A ). The association of oncogenic mutations and age at diagnosis among the patients were assessed by univariate logistic regression. Eight genes were significantly correlated with age ( Figure 3B ), namely, PIM1 (OR = 1.015, 95% CI = 1.005–1.026, P = 0.002), MYD88 (OR = 1.022, 95% CI = 1.011–1.033, P <0.001), BTG2 (OR = 1.014, 95% CI = 1.003–1.025, P = 0.015), CD79B (OR = 1.030, 95% CI = 1.016–1.044, P <0.001), TET2 (OR = 1.016, 95% CI = 1.003-1.030, P = 0.014), BTG1 (OR = 1.014, 95% CI = 1.001–1.028, P = 0.037), CREBBP (OR = 1.016, 95% CI = 1.002–1.031, P = 0.025), and TBL1XR1 (OR = 1.019, 95% CI = 1.002–1.036, P = 0.025). Meanwhile, the same results were observed on MYD88 ^L265P mutation alone (OR = 1.024, 95% CI = 1.011–1.037, P <0.001), or with CD79B mutation (OR = 1.038, 95% CI = 1.013–1.064, P = 0.002) ( Figure 3C ).

All 1,150 patients were genetically classified (21), namely, MYD88-like (191, 16.6%), SOCS1/SGK1-like (148, 12.9%), BCL2-like (222, 19.3%), and NEC-like subtype (589, 51.2%). MYD88-like subtype was significantly associated with age (OR = 1.024, 95% CI = 1.012–1.036, P <0.001) ( Figure 3D ). Based on biological functions, mutations involving BCR/NF-κB signaling (OR = 1.017, 95% CI = 1.009–1.026, P <0.001), B-cell differentiation (OR = 1.012, 95% CI = 1.003–1.020, P = 0.009), and histone acetylation (OR = 1.012, 95% CI = 1.002–1.023, P = 0.018) were significantly associated with age ( Figure 3E ). The total mutation load in the coding genome also showed positive correlation with age in 340 patients with WGS or WES data (P = 0.005, r = 0.152) ( Figure 3F ).

The abundance of immune cells in tumor microenvironment was evaluated using RNA sequencing data. The association between the abundance of immune cells and age was estimated using Spearman’s rank correlation. With abundance of CD4⁺ naïve T and CD8⁺ naïve T cells decreased at age (P = 0.006 and P = 0.006), and abundance of exhausted T cells and macrophages increased with age (P = 0.002 and P = 0.046) ( Figure 4A ).

As exhausted T cells and macrophages play essential roles in immune evasion, their relations with age-associated mutations and effects on tumor progression were further assessed. All the patients were subsequently divided into two groups according to the median of abundance of the two immune cells, respectively. CD79B mutations were found more frequently mutated in the exhausted T-high group and macrophage-high group (P = 0.009 and P = 0.030) ( Figure 4B ). As revealed by GSEA, BCR signaling pathway was upregulated according to CD79B mutations ( Supplementary Figure 1A ). IL10, one of the key cytokines associated with T-cell exhaustion and macrophage M2 polarization, was found significantly higher in CD79B mutation patients (P = 0.001) ( Supplementary Figure 1B ). The abundance of macrophages showed positive linear correlations with the abundance of exhausted T cells (P <0.001) ( Figure 4C ).

Correspondingly, GSEA analysis revealed negative regulation of T-cell immunity (T-cell differentiation, lymphocyte activation, and T-cell proliferation, and IL-2 production), upregulation of PD-1 signaling, inhibitory cytokines production and signaling, and also the recruitment of the inhibitory immune cells in exhausted T-high group, as compared to exhausted T-low patients ( Figure 4D ). IL10, but not TGFB1 (TGFβ), showed positive linear correlations with the abundance of exhausted T-cells (P <0.001) ( Supplementary Figure 2A ). When analyzing the expression levels of inhibitory receptors and their ligands and the abundance of exhausted T cells, inhibitory receptors, namely, LAG3, PDCD1 (PD-1), CD244 (2B4), HAVCR2 (TIM3), CTLA4, TIGIT, CD160, and BTLA, and ligands to inhibitory receptors, namely, TNFRSF14 (HVEM), LGALS9, PDCD1LG2 (PDL2), CD48, and CD274 (PDL1) showed positive linear correlations with the abundance of exhausted T cells (all P <0.05) ( Figure 4E and Supplementary Figure 2B ). Among them, the expression of CD244 and PDCD1 was increased with age (P = 0.016 and P = 0.031) ( Supplementary Figure 2C ).

According to the abundance of macrophages, GSEA analysis indicated that upregulation of macrophage recruitment and activation, immunoregulatory interactions of lymphoid and non-lymphoid cells, inflammation, collagen degradation, PD-1 signaling, angiogenesis, aging, and also inhibitory cytokine production and signaling in macrophage-high group, as compared to macrophage-low group ( Figure 5A ). Markers of macrophage M1 and M2 were further compared, showing that M2 markers, namely, CD163 (log₂foldchange = 2.693, P <0.001), MSR1 (log₂foldchange = 2.785, P <0.001), and MRC1 (log₂foldchange = 1.807, P <0.001) were significantly upregulated in macrophage-high group ( Figure 5B ). Meanwhile, cytokines IL10, CSF1, IL1B, IL6, and TGFB1, which were involved in M2 polarization, showed positive linear correlations with the abundance of macrophages (all P <0.001) ( Supplementary Figure 3A ). Since macrophages may contribute to T-cell exhaustion through the expression of ligands to inhibitory receptors, the relationships between the expression levels of inhibitory receptors and their ligands with the abundance of macrophages were analyzed. Among them, ligands to inhibitory receptors, namely, PDCD1LG2, CD274, LGALS9, TNFRSF14, PVR, and CD48, and also inhibitory receptors, namely, HAVCR2, CD244, LAG3, VSIR (VISTA), PDCD1, CTLA4, CD160, and ENTPD1 (CD39), showed positive linear correlations with the abundance of macrophages (all P <0.05) ( Figure 5C and Supplementary Figure 3B ). When analyzing the expression levels of chemokines and the abundance of macrophage, chemokines related to macrophage recruitment CCL8, CXCL16, CCL2, CCL18, CCL3, CCL4, CCL5, CXCL12, CXCL8, CX3CL1, CCL14, and CCL26 showed positive linear correlations with the abundance of macrophages (all P <0.05) ( Figure 5D ). Among these chemokines, the expression of CCL3 and CCL5 was increased with age (P = 0.024 and P = 0.045) ( Supplementary Figure 3C ).