Clinical characteristics of the study group and expanded deconvolution cohort

Baseline clinicopathological characteristics of BPDCN cases included in the current study are briefly summarized in Table 1. Clinical outcomes, when available, reflected previous dismal observations in BPDCN with a median progression-free and OS of 8 and 12 months, respectively. Following our previous study, RNA-seq data were deconvoluted according to single-cell DC and monocyte datasets (Supplementary Fig. 1), and previous observations were recapitulated/extended to new samples including differentially mutated genes (Supplementary Fig. 2). Additionally, mutational landscape, MutSigCV analysis, and immunohistochemical profiles were extended to new cases as described (Supplementary Tables 1, 2, and 3) [6]. Tumor mutational burden (TMB) was confirmed to be significantly higher (p = 0.0055) among C1-BPDCN and the set of differentially mutated genes was updated (Supplementary Fig. 3a, b).

Table 1 Baseline clinicopathological characteristics of the study group.Epigenetic profiling reveals significant deregulation of key regulatory pathways through loss of DNA methylation compared to dendritic cells

To assess epigenetic processes contributing to the malignant transformation from DCs to BPDCN and to allocate the entity within the spectrum of blood cells, we performed a principal component analysis (PCA) of BPDCN and various cell types of the peripheral blood (Fig. 1) [26]. Expectedly, we observed a clear segregation of blood cell types and BPDCN, but beyond this, a marked difference between BPDCN and previously profiled DC subsets became apparent (Fig. 1a) [27]. This was further reflected in significantly higher global mean DNA methylation levels in both early and mature DCs, signifying the extent of DNA methylation loss during transformation of DCs towards BPDCN (Fig. 1b). Enrichment analysis of differentially methylated regions (DMRs)(CpG FDRs <0.01, absolute difference above 0.3, enrichment analysis against HALLMARK and REACTOME) revealed significantly reduced DNA methylation in genes involved in extracellular matrix organization, collagen modulation, and neuronal systems. Oncogenic driver processes affected by these altered profiles included KRAS signaling (Fig. 1c, d).

Fig. 1: Epigenetic profiles of BPDCN and sorted hematopoietic cell populations.

a First and second principal components of the 10,000 most variable DNA methylation sites in BPDCN (C1 and C2) and various hematopoietic cell types (B B lymphocytes, Granulo granulocytes, HSC hematopoietic stem cells, iDC immature dendritic cells, mDC mature dendritic cells, Mono monocytes, NK natural killer cells; ellipses show 95% confidence intervals of multivariate normal distribution). b Average genome-wide DNA methylation level (beta values) of BPDCN and various cell types and of BPDCN cluster C1 and C2. Individual estimates are shown as dots and cell type-specific distributions are shown as box- and violin-plots; significant differences against BPDCN were assessed by unpaired Wilcoxon test and significant levels are indicated by asterisks (*p < 0.05, **p < 0.01, and ***p < 0.001). c Enrichment analysis of differentially methylated regions (DMRs) between BPDCN and dendritic cells (harmonic mean of the individual CpG FDRs <0.01, absolute difference above 0.3) against HALLMARK and REACTOME gene sets. Only significant gene sets are shown (FDR < 0.1). d Network enrichment against REACTOME for DMRs (as in c); fold changes of DMRs (BPDCN vs DC) are color-scaled (red: higher DNA methylation in BPDCN; blue: higher DNA methylation in DC) and gene sets are denoted by light-brown nodes.

BPDCN is characterized by a DNA methylation profile distinct from its related entities yet borderline cases exist

To assess the potential of DNA methylation analysis in the differential diagnosis between BPDCN and neighboring entities, we performed a comparative analysis of DNA methylation data between our cohort and the BEAT-AML cohort as well as a CMML cohort previously published by Palomo et al. This analysis revealed a separation of samples into entity-specific DNA methylation classes in both an unsupervised (PCA; Fig. 2a) as well as a supervised (PLS-DA; Fig. 2b) approach. Distinction from T-cell acute lymphoblastic leukemia (T-ALL) and malignant melanoma, which was added for comparison, as another entity driven by a UV-related mutational signature, was even more pronounced (Fig. 2j). Intriguingly, we observed four borderline cases (Fig. 2c; BPDCN_01; 15; 27; 37), exhibiting a relevant amount of overlap between both groups in terms of DNA methylation. In contrast to the overall distinction between BPDCN and AML, we found three of these cases to exhibit synchronous concurrent manifestations of AML and/or transformed CMML with pDC features and one of these three cases even presented with genetic features typically encountered in AML with pDC-features (RUNX1 mutation; Fig. 2c, d–i). The patient was allocated closest to the 95% CI cut-off for the AML definition and exhibited an atypical, more immature C2-BPDCN phenotype. Intriguingly, the RUNX1 mutation was observed at an approx. VAF of 50%, despite the obvious phenotypical heterogeneity with typical BPDCN, infiltrates alongside cells more typical of pDC-like AML emerging from an MDS, suggesting a partially shared clonal architecture. While RUNX1 mutations are considered rather rare in BPDCN we observed eleven cases, however, most of these (nine out of eleven) where these mutations were subclonal events and some cases even exhibited syn- and/or metachronous myeloid neoplasms other than BPDCN. Subsequent comparative pathway enrichment analysis against HALLMARK and REACTOME gene sets for most DMRs within gene-body regions revealed a significant enrichment across RHO GTPases, cell migration control, and leukemic stem cell maintenance (HSF1 activation) in BPDCN, whereas promotor regions in BPDCN compared to AML were methylated to a significantly higher degree in epigenetic and transcriptional regulation as well as TP53 regulation and cell cycle control. (Fig. 2m–o).

Fig. 2: BPDCN DNA methylation in comparison to AML, CMML, t-ALL, and melanoma.

a Visualization of the first and second principal components of the 10,000 most variable DNA methylation sites (ellipses show 95% confidence intervals of multivariate normal distribution). b Partial-least squares discriminant analysis (PLS-DA) of adjusted beta values. c First and second principal components of the comparison between BPDCN and AML (RUNX1 wild-type and mutated samples highlighted differently). The four BPDCN cases falling inside the 95% confidence interval of the AML data are labeled. d–j a prototypical borderline case with both typical BPDCN as well as AML with pDC-like features. d Morphology of the neoplastic infiltrate within the lymph node resembles acute leukemia with polymorphic blast-like cells of variable size (H&E, 40×). e Uniform expression of CD123 initially led to the inclusion of BPDCN into the differential diagnosis (CD123, 40×). f Further immunophenotypic work-up revealed several atypical features, reminiscent of AML with partial pDC phenotype, including variable expression of CD33 in a significant fraction of the malignant infiltrate (CD33, 40×), yet only partial expression of CD56 (g; CD56, 200×) and CD117 (h; CD117, 200×). i The bi-phenotypic character of the infiltrate is further underlined by a strong CD34 expression of a minor fraction of the blast-like cells alongside the vascular structures, resembling the pDC-like AML phenotype component, whereas the majority of blasts resemble phenotypically characteristic BPDCN cells. j Two-dimensional density plots of average CpG site DNA methylation in BPDCN vs AML, CMML, t-ALL, and melanoma (low density: orchid; high density: yellow/orange). k Average genome-wide DNA methylation level (beta values) of BPDCN, AML, CMML, t-ALL, and melanoma; for BPDCN subcluster estimates are shown as well. l DNA-methylation-based mitotic clock (epiCMIT) estimates for each entity and BPDCN subtypes. m Network enrichment against REACTOME for DMRs between BPDCN and AML; fold changes are color-scaled (red: higher DNA methylation in BPDCN; blue: higher DNA methylation in AML) and gene sets are denoted by light-brown nodes. n Pathway enrichment against REACTOME gene sets of gene-associated CpGs between BPDCN (blue) and AML (orange). o Pathway enrichment against REACTOME gene sets of promotor-associated CpGs between BPDCN (blue) and AML (orange). If not stated differently, differences between BPDCN and the four other entities were assessed by unpaired Wilcoxon test, and significant levels are indicated by asterisks (*p < 0.05, **p < 0.01, and ***p < 0.001).

The BPDCN genome is characterized by an exceptional degree of DNA methylation loss and epigenetic signs of mitotic stress

Global DNA methylation is commonly diminished in malignant compared to healthy cells, being stably maintained in the latter. In this line, BPDCN revealed the highest methylation loss compared to AML, CMML, T-ALL, and melanoma irrespective of the C1 or C2 subtype (Wilcoxon rank-sum test, p < 0.001) (Fig. 2k) [8, 28,29,30,31]. Intriguingly, the aggressive clinical nature of BPDCN was further reflected in our epiCMIT analysis, as a DNA methylation-based mitotic clock, which recapitulates the proliferative history of a given tumor sample [32]. epiCMIT predicted unveiled an accelerated mitotic history in BPDCN compared to CMML and even AML, matched only by T-ALL (Fig. 2l). Although an independent prognostic impact of a high epiCMIT score was observed for a wide range of blood cancers, we observed no such trend in our cohort, plausibly attributable to the limited sample size.

DNA methylation of tumor suppressor genes is highly deregulated in BPDCN compared to CMML and AML

In keeping with a substantially deregulated DNA methylation profile in BPDCN and global loss of DNA methylation, we observed promotor regions of tumor suppressor genes to be methylated to an exceptionally high degree, compared to CMML and AML (Fig. 3a). A significantly lower level of gene-body associated methylation in the same genes was identified (Fig. 3b), indicating a substantial oncogenic impact of deregulated DNA methylation in BPDCN (manually selected candidate TSGs, Fig. 3c, for exhaustive information on deregulated TSGs see Supplementary Fig. 2). Strikingly, we hereby observe an increasingly deregulated DNA methylation profile from CMML to BPDCN in parallel to the aggressiveness of the entity’s clinical behavior.

Fig. 3: DNA methylation of tumor suppressor genes in BPDCN, AML, and CMML.

a Average promoter DNA methylation (beta values) of tumor suppressor genes in BPDCN, AML, and CMML. b Average gene body DNA methylation (beta values) of tumor suppressor genes in BPDCN, AML, and CMML. c Average promoter and gene body DNA methylation in selected tumor suppressor genes (for a complete representation of significantly divergent TSGs see Supplementary Fig. 2). Differences between BPDCN and the AML/CMML were assessed by unpaired Wilcoxon test and significant levels are indicated by asterisks (*p < 0.05, **p < 0.01, and ***p < 0.001).

DNA methylation patterns and gene expression signatures differentiate C1 and C2 subtypes in BPDCN and shape a JAK/STAT-driven profile in C2-BPDCN

Building on transcriptional BPDCN subtype classification, we performed a differential DNA methylation analysis. Alongside differentially mutated genes in BPDCN, we discovered 114 probes, which were differentially methylated between C1 and C2 (p < 0.0001; 10,136 probes with p < 0.01), corresponding to a relatively similar methylome, in keeping with our above PCA/PLS-DA (depicted as scales beta values in Fig. 4a; beta values see Supplementary Table 4).

Fig. 4: Differential DNA methylation, expression, and mutation patterns in BPDCN subtypes C1 and C2.

a Top part of the heatmap shows mutational patterns of 14 genes previously identified as significantly enriched between the two clusters [6]. The bottom part shows beta values (DNA methylation levels) of CpG sites where an absolute mean difference above 0.25 was observed between C1 and C2. b Volcano plot of beta values showing log2 fold-changes and p values with gene annotations for significantly different CpGs (p < 0.0001). c Gene set enrichment analysis results of DNA methylation data against REACTOME gene sets (p < 0.05). d Volcano plot of expression profiles showing log2 fold-changes and p values with gene annotations for significant differentially expressed tumor suppressors genes, oncogenes, and genes involved in cell adhesion/cell cycle (p < 0.0001). e Gene set enrichment analysis of RNA-seq data against REACTOME gene sets (padj < 0.001, absolute enrichment >0.3). f Transcription factors with significantly different inferred activity (p < 0.05) in C1 (light blue) and C2 (dark blue). g Pathway activities in C1 (light blue) and C2 (dark blue) with significantly different pathway activities shown in red font (p < 0.05). h, i FFPE-ATAC-seq estimates of chromatin accessibility for CDK6 and STAT5A for four BPDCN cases belonging to subcluster C1. CDK6 is located on the minus strand and STAT5B on the plus strand, respectively.

To gain insight into the effect of differential DNA methylation of both subtypes, we assigned differentially methylated probes to their respective genes (Fig. 4b). GST (Fig. 4c) revealed enrichment in elevated DNA methylation levels for interleukin signaling genes and prominent members of JAK/STAT signaling (STAT5B) in C1-BPDCN, whereas C2-BPDCN samples exhibit significantly pronounced DNA methylation of posttranslational modifications and metabolic processes (e.g., vitamin and heme metabolism).

Next, we performed differential gene expression analysis (Fig. 4d, e). Hereby, we observed an induction of innate and adaptive immunological processes alongside an upregulation of extracellular matrix interactions in C2 cases. Further, we found several prominent candidate genes, including STAT5A, CDK6, CCR4, CCND2, and FOXO1 to be expressed at significantly higher levels in C2 cases. Corresponding to higher gene expressions in C2, CDK6 (p = 5.3 × 10−4, log2 fold-change 1.52) and STAT5A (p = 5.5 × 10−3, log2 fold-change 1.19) had promotor-associated sites that were significantly higher methylated in C1-BPDCN leading to relative transcriptional inactivation in C1. Correlation of DNA methylation at promotor and gene body sites and relative RNA-seq-derived gene expression data yielded a substantial number of significant correlations. Findings regarding the enrichment of gene promotor methylation in particular TSGs are only partially recapitulated on the transcriptional level (Supplementary Table 5 and Supplementary Fig. 4a). Various statistically significant correlations between DNA-methylation at either type of site and gene expression, regardless of potential confounding factors, including copy number variants and mutations, affecting prominent TSGs and oncogenes are identified through both a focused and genome-wide approach (see Supplementary Fig. 4b, c). Significantly correlated candidates include BCL2, ALK, GNAQ, and RUNX1.

To further focus our observations on potential therapeutic applicability, we inferred TF activities from bulk RNA-seq data employing CollecTRI as a resource for TF pathway activity inference with PROGENy. Hereby we observed the transcriptional mirror image of the divergence in DNA methylation profiles between C1 and C2-BPDCN. In particular, we observed significantly higher activity of NFkB (driven by FOXC1, NFKB1, and NFKB; p < 0.05) and more strikingly JAK-STAT (predominantly driven by STAT3, STAT1, and STAT5A; p < 0.05) associated TF in C2 and an EZH2 dependence in C1-BPDCN (Fig. 4f, g). We supplemented these observations by FFPE-ATAC-seq of four C2-BPDCN, which revealed substantial chromatin accessibility, in keeping with our epigenetic and transcriptional findings (Fig. 4h, i).

Tumor microenvironment by MethylCIBERSORT and immunohistochemistry reveals distinct immunological subtypes correlated with tumor genomics

To evaluate the TME, we conducted a MethylCIBERSORT analysis on all 45 high-quality genome-wide DNA methylation profiles. From inferred relative abundance of T-cell subpopulations, B-cells, natural killer (NK)-cells neutrophils, monocytes, eosinophils, and stromal cells two TME classes (IC1 and IC2) were predicted. The smaller class (IC1) exhibited a depletion in monocytes, B- and NK-cells, an enrichment in Tregs, and a trend towards higher counts of neutrophils and cytotoxic T cells (Fig. 5a). Additional deconvolution revealed a significantly higher proportion of keratinocytes in skin samples compared to non-skin samples (p = 0.036). These observations were subsequently validated via IHC for tumor-infiltrating T cells and monocytes with a significant correlation between relative distributions of cell populations (Fig. 5b–g; monocytes: Pearson correlation coefficient = 0.6437, p = 1.84 × 10−6 and T-cell populations: Pearson correlation coefficient = 0.3809 p = 0.0098). The correlation of TME with mutational profiles (restricted to the eight genes mutated in at least 30% of BPDCN cases) is presented in Supplementary Fig. 5. We observed enrichment in mutations affecting CDH11, ERBB2, ASXL1, EP300, KMT2C, JAK2, SMAD2, NOTCH1, and DNMT3A (Fig. 5h, Supplementary Fig. 6a) and a trend towards a higher TMB (Supplementary Fig. 6b) in IC1. Further, IC1 patients had significantly shorter progression-free survival (p = 0.044) and a trend towards inferior OS (p = 0.14) (Fig. 5i, j). Similar characteristics in terms of TME composition were observed by immunohistochemistry and MethylCibersort with enrichment of T-cells in IC1 and monocytes in IC2 (Supplementary Fig. 6c–f). epiCMIT was applied to estimate the history of proliferative stress/DNAm age in both immunological clusters and revealed significantly higher proliferative stress in IC2 compared to IC1, which resembles DNA methylation-based pre-aging in this subgroup (p = 0.038; Supplementary Fig. 6g). A partial overlap between TME classes IC1/2 and transcriptional clusters C1/C2 with a strong enrichment of IC2 cases in C1 BPDCN patients was observed. An integrated summary of immunohistochemical and molecular features including cluster allocation is provided in Supplementary Table 3. Of particular clinical relevance and in keeping with the impact on survival we identify substantial differences in terms of overall response rates between C1/C2 (66% vs 31%) and IC1/IC2 (12% vs 62%) (Fig. 5f, g). Previous reports identified a correlation between DNA methylation and chronological age [33, 34]. We observed a trend towards higher levels of tumor-infiltrating T-cells in C2-BPDCN samples despite the overall higher TMB in C1 patients, regardless of patient age or location of the tumor sample by both MethylCIBERSORT and IHC (Supplementary Fig. 7a, b).

Fig. 5: Tumor immune composition by MethylCIBERSORT identifies clusters of divergent immunogenicity.

a DNA methylation data were deconvoluted according to immune cell populations (MethylCIBERSORT). This revealed two different types of BPDCN (Treg and CD14 driven; named IC1 and IC2, respectively) presenting with significantly differing immune cell subsets in regard to the markers CD14 (monocytes/macrophages), CD19 (B-cells), CD4 (T-helper cells) as well as fibroblasts, NK-cells and T-regulatory cells (T-regs). b–e A borderline BPDCN/AML pDC-like case analyzed by MethylCIBERSORT and a comparative immunohistochemical assessment of the tumor microenvironment is presented. b H&E staining reveals a cutaneous infiltrate covered by an intact epidermis. c Giemsa staining reveals small blastoid cells with partly roundish occasionally monocytoid nuclei, small nucleoli, and weakly basophilic cytoplasm with increased mitotic activity. d Staining for myeloid peroxidase reveals expected negativity in the malignant infiltrate alongside a few positive, tumor-infiltrating myeloid cells. e However, CD14-Expression highlights both a typical negative BPDCN population, as well as a relevant monocytoid population, including few tumor-infiltrating monocytes alongside a larger subgroup of malignant cells. f Mosaic plot visualizing the overlap between transcriptional phenotype C1/C2 and IC1/IC2. g Overall, complete and partial response rates according to C1/C2 and IC1/IC2. h Oncoplot displays mutational patterns of 9 genes that were found to be more differentially mutated between IC1 and IC2. Additionally, the heatmap illustrates TME cell proportions for each individual sample. i, j Progression-free (PFS) and Overall survival (OS) analysis for patients with available clinical follow-up according to IC1 vs IC2 identifies a significant inferior prognostic impact for the IC1 subtype regarding PFS accompanied by a trend towards inferior OS.

DNA methylation-based clusters cannot be fully recapitulated in transcriptional BPDCN subtypes

The most 5600 variable CpG probes were subject to discovering DNA methylation subtypes using K-means clustering. The optimal number of clusters was determined using the average silhouette method and gap statistics. Both methods agreed on the optimal number of clusters (n = 2). Upon unsupervised cluster analysis, we observed no statistically significant recapitulation of transcriptional BPDCN subtypes C1 and C2. Some overlap was, however, found in shared genomic features (EP300 and ATRX mutations in C1 and MethC1), consistent with higher TMB in C1 and MethC1. A significantly more pronounced overlap was apparent upon correlation analysis between DNA methylation-based cluster allocation and immunological clusters according to our methyCIBERSORT-derived immunoclusters (Supplementary Fig. 8a, b).

Mutational drivers and promotor status of epigenetic regulators shape the proliferative fate and DNA methylation profile of BPDCN

To provide genomic context for DNA methylation profiles obtained within this study, we performed complementary WES on all patients, who were not part of our previous molecular BPDCN study (Fig. 6a). WES confirmed a higher mutational load in C1 BPDCN (Supplementary Fig. 3b), further supporting a non-mutational mechanism driving C2-BPDCN such as the deregulated DNA methylation and transcriptional profile outlined above. The pathophysiological role of epigenetic features within the C2 cluster is supported by more frequent DNMT3A (14% vs 36%) alterations bordering on statistical significance in this limited cohort and the significant mutational enrichment for splicing gene SRSF2 (0% vs 16%; Supplementary Fig. 3a).

Fig. 6: Impact of significant driver mutations and DNA methylation status of epigenetic regulators on proliferative history and global CGI DNA methylation.

a Oncoplot displaying putative driver genes inferred by MutSigCV analysis in addition to tumor mutational burden (TMB; upper bar plot), −log10p values (left bars), and the number of samples harboring mutations in a given gene (right bar). In total 6065 bioinformatically deleterious mutations, in 3970 genes were observed. Of these indels and SNVs, 4999 were missense (82.4%), 321 nonsense (5.3%), 12 non-stop (0.2%), and 487 indel mutations (8.0%). Mutation types are color-coded, and covariates are shown below for each sample (covariate “Other” refers to samples with tissue affected other than skin or bone marrow). b Analysis of which significant mutational driver alterations may confer a higher or lower proliferative capacity to BPDCN cells. Statistically significant associations are delineated in red and blue, respectively. c Comparative analysis regarding proliferative capacity in selected genes’ mutant and wild-type status. d A heatmap depicting the DNA methylation status of gene promoter regions across a set of preselected epigenetic regulators associated directly and indirectly with global DNA methylation, as described [48]. Average global DNA methylation and CpG islands (CGI) are estimated per patient.

Further, we wondered whether specific mutational drivers were associated with the proliferative history of a given case, thereby integrating mutational and epigenetic datasets to identify genetic alterations, which could confer a selective advantage as mitotic accelerators. We observed a trend towards a proliferative increase in cases harboring ARID1A, ZRSR2, JAK3, CIC, and CREBBP mutations and a (significant) decrease in mitotic activity in cases with mutant PPARGC1A, CDH1, SI, RUNX1, SF3B1, MSH2 and SETD2 (Fig. 6b, c). In summary, our findings suggest, that the proliferative potential of BPDCN is determined by mutational driver events. To unravel additional determinants, that shape the globally demethylated landscape in BPDCN, we investigated the promotor DNA methylation status of a set of genes associated with epigenetic regulation (Fig. 6d). Hypermethylation of the promotor region (promoter CGI methylation >0.2) of DNMT1, WT1, MYC, and TET1 was hereby observed in 42.2%, 73.2%, 40.0%, and 53.3% of patients, respectively. Further, we observed a significant correlation between promotor methylation and overall CGI methylation for 22 of the 24 genes (Spearman’s rank correlation, padj < 0.05) signifying a substantial impact on the epigenetic landscape in BPDCN. Only DNMT1 and KMD2B did not show a significant correlation. Of note, a substantial overlap between patients harboring WT1 and MYC promoter methylation was observed in 16/45 cases (35.6%). Given that MYC constitutes a known WT1 target, this hints at a coupled mechanism in BPDCN pathogenesis [35] (Fig. 6b, c).