High Histological Skeletonization as a Measure of Heterogeneity Is Associated with Poor Prognosis of Diffuse Large B-Cell Lymphoma

1. Introduction

Diffuse large B-cell lymphoma (DLBCL) is one of the most frequent lymphomas, accounting for approximately 25% of non-Hodgkin lymphomas [1].

DLBCL usually presents as a symptomatic lymph node enlargement in the neck or truck, but extranodal involment is also possible. In 30% of the cases, the patients manifest fever, night sweats, or weight loss, known as “B” symptoms [1]. The incidence of DLBCL is around 3-7 cases per 100,000 persons per year [2].

At presentation, in around 60% of the patients, the disease presents as advanced stage III-IV and elevated lactate dehydrogenase (LDH) in serum. Bone marrow is infiltrated in up to 30% of the cases. In the case of nodal DLBCL, secondary involvement of other organs such as the liver, kidney, lung, bone marrow, and central nervous system can occur [3,4,5,6].

From a histopathological point of view, DLBCL is characterized by a diffuse proliferation of large transformed B lymphocytes with prominent nucleoli and basophilic cytoplasm with abundant proliferation. The neoplastic lymphocytes resemble centroblasts or immunoblasts, and other morphologies, such as anaplastic, multinucleated, clear cells, epithelioid, and spindle cell shapes can be found [6,7,8,9,10,11].

At diagnosis, immunohistochemistry and/or flow cytometry are used to determine the immunophenotype at the protein level, which include B-cell markers (CD20, CD79a), CD10, BCL6, MUM1 (IRF4) for assessing the cell-of-origin Hans’classification, and BCL2, Ki67, MYC, CD5, and p53 for evaluating the proliferation status and prognostication [12,13,14,15,16,17,18,19]. Both the immunophenotype and cytological morphologies are useful in DLBCL diagnosis and prognosis [20].

Cytogenetic features are also analyzed despite no typical or diagnostic change. MYC, BCL2, and BCL6 gene rearrangements are routinely used to subclassify DLBCL into the DLBCL/high-grade B-cell lymphoma categories associated with a more aggressive clinical evolution [21,22,23,24,25,26,27,28,29,30].

The prediction of DLBCL prognosis using artificial intelligence and digital pathology has gained interest recently. For example, Jeong Honn Lee et al. recently retrospectively collected 251 slide images and performed deep learning and multi-modal prediction with contrastive learning and feature extraction that provided a performance with an area under the curve of 0.8 in the ROC analysis [31]. Other recent works include the work of Vrabac D et al. [32], Li Z. et al. [33], and Ma D. et al. [34] that also focused on the evaluation of histological images. However, the focus on topological features is less studied.

This study analyzed the topological characteristics of DLBCL using hematoxylin and eosin slides in a series of 153 patients, including 110 patients with DLBCL and 43 patients with reactive lymphoid tissue. The study found that DLBCL was characterized by lower skeletonization than reactive lymphoid tissue. In DLBCL, high skeletonization is associated with poor prognosis and enrichment of immuno-oncology markers.

2. Materials and Methods

2.1. Patients and Samples

A series of 153 patients were selected from the Department of Pathology, Tokai University, School of Medicine. The series included 110 cases of diffuse large B-cell lymphoma and 43 cases of reactive lymphoid tissue and reactive tonsils. The cases were diagnosed according to the current lymphoma classification [3,35,36,37,38,39]. The diagnosis was based on the clinicopathological characteristics of the patients, including the evaluation of the hematoxylin and eosin (H&E) stainings, immunophenotype, and molecular pathology data when required [37,38].

This study was conducted according to the guidelines of the Declaration of Helsinki and the World Medical Association for medical experimentation version 24. This study was approved by the Tokai University ethical committee (Institutional Review Board, IRB20/156).

2.2. Skeletonization and Area Measurement

Skeletonization reduces binary objects to 1 pixel-wide representations. The skeletonization technique is useful for feature extraction and represents an object (image)’s topology.

DLBCL images were transformed from hematoxylin and eosin (H&E) to a skeletonized image. H&E had previously been stained using an automated slide stainer (Tissue-Tek Prisma Plus, Sakura Finetek Japan Co., Ltd., Tokyo, Japan) as previously described [40]. H&E images were digitalized using a digital slide scanner and whole slide imaging system (NanoZoomer S360, #C13220-01, Hamamatsu Photonics K.K., Hamamatsu, Japan) and visualized using NDP.view2 image viewing software (#U12388-01, Hamamatsu Photonics K.K.). From each case, digitalized jpeg images were exported at 200x magnification and 150 dpi. For better image processing and performance, each image was split into image patches of 224× 224× 3 resolution using PhotoScape v3.7 (http://photoscape.org/; last accessed on June 20, 2026).

The skeletonization process extracted the centerline of an object while preserving the topology and Euler number (also known as the Euler characteristic). The process was initiated on grayscale image patches. The object of interest were dark threads against a light background. Skeletonization required a binary image with white foreground (1) and black background (0) pixels. The complement of the original grayscale image was used. Then, the image was binarized, and skeletonization was performed. In this study, pruning was not performed. Finally, the area of the skeleton was measured for each image patch.

When required only for visualization, the skeleton was overlayed on the original image, and small spurs that appeared on the skeleton were pruned.

2.3. Computer Code

Computer vision analysis was performed using Matlab (R2026a Update 3 (26.1.0.3276743), May 27 2026). This study used the ‘bwskel’ function. Morphological operations on binary images using ‘bwmorph’ were also tested, such as the ‘skel’ function that created comparable skeleton images as the ‘bwskel’ function. However, this study finally used the ‘bwskel’ function for the final analysis. The skeletonization codes for one image and a folder with multiple images are shown in Table 1. In this study, small spurs were not pruned.

2.4. Immunohistochemistry

Immunophenotype was performed using conventional immunohistochemical techniques using a Leica Bond Max II automated immunohistochemistry stainer following the manufacturer’s instructions (Leica Biosystems K.K., Tokyo, Japan). Immunohistochemistry reagents included Bond Dewax solution (#AR0084), Bond epitope retrieval solution 2 (#AR0087), primary antibody diluent (#AR9352), and Bond polymer refine detection kit (#DS9800). Slides were mounted using a fully automated glass coverslipper (Leica CV5030, Leica Biosystems).

All primary antibodies were purchased from Novocastra (Leica Biosystems), including CD20 (clone L26, #CD20-L26-L-CE), CD79a (clone JCB117, #CD79A-599-L-CE), CD3 (clone LN10, #CD3-565-L-CE), CD5 (clone 4C7, #CD5-4C7-L-CE), BCL2 (#BCL-2-L-CE), Ki67 (clone K2, ACK02), and MYC (anti-c-Myc antibody [Y69], ab32072, Abcam K.K., Tokyo, Japan).

The cell of origin characterization using the Hans’ algorithm [12] used the combination of CD10 (clone 56C6, #CD10-270-L-CE), BCL6 (clone LN22, # BCL-6-564-L), and MUM1 (clone EAU32, #MUM1-L); more than 30% of positivity by neoplastic B lymphocytes was considered positive for calculating Hans’ algorithm.

Immune microenvironment images were retrieved from our previous publications [41,42,43,44], and immunohistochemistry was reanalyzed. The primary antibodies that targeted the immune microenvironment, immune checkpoint and cell cycle were the following: Ki67 (Cell proliferation, clone MM1, Novocastra, Leica), IL10 ( Immuno-oncology, #LS-B7432, Lifespan Bioscience), PD-L1 (CD274) (Immuno-oncology, clone E1J2, Cell Signaling), CSF1R (Immuno-oncology, clone FER216, CNIO), CD163 (Tumor-associated macrophages, clone 10D6, Novocastra, Leica), CASP8 (Active subunit p18, clone 11B6, Novocastra), TNFAIP8 (Apoptosis, #14559-MM0, Sino Biological), LMO2 ( Hematopoietic development, clone 299B, CNIO), MYC (Proto-oncogene, clone Y69, Abcam), MDM2 (Proto-oncogene, IF2, Invitrogen), CDK6 (Cell cycle, clone 98D, CNIO), E2F1 (Cell cycle, clone Agro368V, CNIO), and TP53 (Cell regulation, DO-7, Novocastra, Leica).

2.5. In Situ Hybridization

Characterization of the gene rearrangement and translocation status was performed using DNA fluorescent in situ hybridization (FISH) break-apart probes targeting MYC, BCL2, and BCL6 genes (Vysis, Abbott Japan Co., Ltd., Tokyo, Japan).

RNA in situ hybridization to detect infection by Epstein–Barr virus (EBV)-targeted Epstein Barr virus encoded RNA (EBER) using Bond in situ hybridization probes (#PB0589, Leica Biosystems).

2.6. Gene Expression

RNA was extracted from formalin-fixed paraffin-embedded (FFPE) tissue sections of 10 micrometer thickness using the RNeasy FFPE Kit (50) (#73504, Qiagen). Before RNA extraction, the histology was revised, and cases with at least 80% of neoplastic areas were selected. After quality control, the RNA was used to quantify gene expression levels using the nCounter® PanCancer Profiling Panel (Bruker Spatial Biology, Bruker Corporation). The tissue sections were outsourced to Celgene Corporation (New Jersey, USA) for this analysis. This panel included 730 immune-oncology genes and 40 housekeeping genes. The analysis was the standard, including the calculation of log2 fold change and -log10 P values to create a volcano plot and later the heatmap, and hierarchical clustering using the one minus pearson correlation and average linkage method.

2.7. Statistical Analysis

Comparison between groups used nonparametric tests, including Mann–Whitney U test for 2-independent samples test and Kruskal-Wallis H test for several (≥3) independent samples. Crosstabulation was used to correlate qualitative variables with the Pearson chi-square test and Fisher’s exact test. Overall survival was calculated from the time of diagnosis and biopsy to the date of death or last follow-up using the Kaplan–Meier and log-rank tests. A p-value of less than 0.05 was considered statistically significant.

3. Results

3.1. Evaluation of Skeletonization in Test Images

Skeletonization reduces binary objects to 1 pixel-wide representations. The skeletoniza-tion technique is useful for feature extraction and represents an object (image)’s topology. Skeletonization area measurement was tested on a set of 8 images to evaluate the usefulness of the ‘bwskel’ function. The area values from images 1-8 were as follows: 0, 40, 100, 215, 526, 857, 1786, and 2603, respectively. Therefore, the function was measured correctly (Figure 1).

3.2. Evaluation of the Skeletonization Area in DLBCL

3.2.1. Difference Between DLBCL and Reactive Lymphoid Tissue

The skeletonization technique is useful for feature extraction and represents an object’s topology. The DLBCL skeleton area ranged from 3396.12 to 8176.94, with a mean of 5665.20 ± 1012.82 and median of 5591.65. In comparison to reactive lymphoid tissue, DLBCL was characterized by lower skeletonization area, 5665.20 ± 1012.82 vs 6341.16 ± 548.23 respectively (p < 0.001) (Figure 2).

3.2.2. Correlation Between Skeletonization and Overall Survival of the DLBCL Patients

The skeletonization area correlated with the overall survival of patients with DLBCL and the cases were classified into 2 groups: low skeletonization (45/110, 40.9%) and high skeletonization (65/110, 59.1%). Overall survival analysis using the Kaplan-Meier and log rank test showed that high skeletonizaton cases were associated with worse overall survival than cases with low skeletonization (log rank (Mantel-Cox) P value = 0.002). High skeletonization patients had a hazard risk of 2.484 (95% CI for Exp(B) = 1.362–4.531) (P = 0.003) (Figure 3).

Examples of reactive, high, and low skeletonization in DLBCL are shown in Figure 4, Figure 5, Figure 6, Figure 7 and Figure 8.

3.2.3. Correlation Between Skeletonization and Clinicopathological Features of DLBCL Patients

The skeletonization was correlated with several clinical and pathological variables. High skeletonization correacted with higher EBER positivity (P < 0.001), lower CD5 positivity (P = 0.006), lower E2F1 percentage (P < 0.001), lower BCL2 percentage (P = 0.010), lower ISY1 percentage (P = 0.012), and lower TNFAIP8 (P = 0.014) (Table 2, Table 3 and Table 4 and Figure 9).

3.2.4. Gene Expression Analysis

Gene expression analysis using a pan-cancer immuno-oncology panel of 730 genes was available in a subset of 30 cases, including 18 cases with high skeletonization and 12 of low skeletonization. Conventional biostatistics was performed, including the calculation of differentially expressed P value and fold change. The analysis showed that 166 immuno-oncology genes were upregulated in the high skeletonization group and 8 genes upregulated in the low skeletonization group.

The first 10 most significant in the high skeletonization group were C4BPA (complement), MUC1 (CD molecules), IL13RA1 (chemokines, T-cell functions), RRAD (cell functions), CXCL1 (chemokine, regulation), PLAUR (CD molecule), IL22 (anti-inflammatory cytokine), TLR3 (innate immune response, Toll-like receptor), and PVR (regulation of immune response)

The 8 genes upregulated in the low skeletonization group were BLNK (humoral immune response), TNFRSF13B (TNF superfamily), CD3EAP (adaptive immune response), ILF3 (chemokine), CD1C (T-cell function), CD37 (adaptive immune response), POU2F2 (humoral immune response), and TLR9 (innate immune response, toll-like receptor) (Figure 10).

4. Discussion

DLBCL is the most common histological subtype of non-Hodgkin lymphoma (NHL) and accounts for approximately 35 of NHL [2,3,9]. The molecular pathogenesis is complex and a multistep process initially initiated by the neoplastic transformation of germinal center or post-germinal B lymphocytes. Most of the steps of neoplastic transformation are still unknown. One of the most known oncogenic events is the rearrangement of MYC, BCL2, and BCL6 genes. This study included the study of rearrangements in most of the cases. However, the rearrangement status did not correlate with skeletonization.

The diagnostic category of DLBCL was previously used to describe a group of lymphomas with heterogeneous morphology, genetics, and biological behavior. There are several distinct large B-cell lymphoma variants, including primary mediastinal large B-cell lymphoma (PMBL), T-cell-rich large B-cell lymphoma, intravascular large B-cell lymphoma, lymphomatoid granulomatosis, Epstein-Barr virus (EBV) positive DLBCL, DLBCL associated with chronic inflammation, fibrin-associated large B-cell lymphoma, primary large B-cell lymphoma of immune-privileged sites (central nervous system, vitreoretinal, and testis), fluid overload-associated large B-cell lymphoma, DLBCL associated with IRF4 rearrangements, and ALK-positive large B-cell lymphoma. This study of 110 cases of DLBCL did not differentiate with all these subtypes but used the DLBCL “not-otherwise specified” (NOS) diagnostic category, and the DLBCL cases were classified according to their primary site as nodal (plus spleen), Waldeyer’s ring, gastrointestinal, and other extranodal. When correlation was made with the skeletonization level, no differences were found. However, skelotonization correlated with the presence of Epstein–Barr virus (EBV): lows, 4/44 (9.1%) vs high 24/64 (37.5%) (P< 0.001). Therefore, our data indicate that in a subset of cases with virus infection, the level of skeletonization (understood as heterogeneity as well) depends on pathogenic mechanisms associated with EBV. Notwithstanding, 72.5% of high skeletonization cases was not associated with EBV. Therefore, other mechanisms may be associated.

Gene expression has been extensively studied in DLBCL and led to the identification of three main types: germinal center B-cell-like (GCB), activated B-cell-like (ABC), and not-otherwise-specified type 3 [45,46,47]. In this study, we used the surrogate immunohistochemical-based cell-of-origin classification of Hans that differentiates 2 groups, GCB and non-GCB. The non-GCB DLBCL is associated with poor prognosis [12,13,48]. We speculated that the high skeletonization group associated with poor prognosis would also correlate with non-GCB subtype. Nevertheless, that was not the case. Other factors may be associated with the histological subtype.

We correlated skeletonization level with other markers of the immune microenvironment and cell cycle and found that high skeletonization also correlated with lower CD5, E2F1, BCL2, ISY1, and TNFAIP8 protein levels (all P values < 0.05). We previously described that high immunohistochemical expression of TNFAIP8 is associated with poor prognosis in DLBCL [49]. To deepen the pathological background of high skeletonization and its association with the prognosis of the DLBCL patients, we performed gene expression analysis using a pan cancer immune profiling panel. We found that 166 immuno-oncology genes were upregulated in the high group and only 8 genes were up-regulated in the low skeletonization group. Therefore, the high skeletonization enriches immuno-oncology markers. This association is important as immunotherapy in DLBCL has gained interest recently [50,51,52,53].

5. Limitations

This study has the limitation of the number of cases. Larger series of cases will be needed to confirm our findings.

6. Conclusions

DLBCL is characterized by lower skeletonization than reactive lymphoid tissue. In DLBCL, high skeletonization is associated with poor prognosis and enrichment of immuno-oncology markers.