Single-Cell Multiomics Reveals Clonal T-Cell Expansions and Exhaustion in Blastic Plasmacytoid Dendritic Cell Neoplasm

The immune system represents a major barrier to cancer progression, driving the evolution of immunoregulatory interactions between malignant cells and T-cells in the tumor environment. Blastic plasmacytoid dendritic cell neoplasm (BPDCN), a rare acute leukemia with plasmacytoid dendritic cell (pDC) differentiation, provides a unique opportunity to study these interactions. pDCs are key producers of interferon alpha (IFNA) that play an important role in T-cell activation at the interface between the innate and adaptive immune system. To assess how uncontrolled proliferation of malignant BPDCN cells affects the tumor environment, we catalog immune cell heterogeneity in the bone marrow (BM) of five healthy controls and five BPDCN patients by analyzing 52,803 single-cell transcriptomes, including 18,779 T-cells. We test computational techniques for robust cell type classification and find that T-cells in BPDCN patients consistently upregulate interferon alpha (IFNA) response and downregulate tumor necrosis factor alpha (TNFA) pathways. Integrating transcriptional data with T-cell receptor sequencing via shared barcodes reveals significant T-cell exhaustion in BPDCN that is positively correlated with T-cell clonotype expansion. By highlighting new mechanisms of T-cell exhaustion and immune evasion in BPDCN, our results demonstrate the value of single-cell multiomics to understand immune cell interactions in the tumor environment.

The immune system represents a major barrier to cancer progression, driving the evolution of immunoregulatory interactions between malignant cells and T-cells in the tumor environment. Blastic plasmacytoid dendritic cell neoplasm (BPDCN), a rare acute leukemia with plasmacytoid dendritic cell (pDC) differentiation, provides a unique opportunity to study these interactions. pDCs are key producers of interferon alpha (IFNA) that play an important role in T-cell activation at the interface between the innate and adaptive immune system. To assess how uncontrolled proliferation of malignant BPDCN cells affects the tumor environment, we catalog immune cell heterogeneity in the bone marrow (BM) of five healthy controls and five BPDCN patients by analyzing 52,803 single-cell transcriptomes, including 18,779 T-cells. We test computational techniques for robust cell type classification and find that T-cells in BPDCN patients consistently upregulate interferon alpha (IFNA) response and downregulate tumor necrosis factor alpha (TNFA) pathways. Integrating transcriptional data with T-cell receptor sequencing via shared barcodes reveals significant T-cell exhaustion in BPDCN that is positively INTRODUCTION Innovations in immuno-oncology, such as immune checkpoint blockade (ICB) therapy, have transformed cancer medicine. ICB and CAR-T cells have led to improved outcomes in various solid tumors and B-cell malignancies, respectively (1), however these methods have been less successful in myeloid leukemias than in other cancers. Various mechanisms have been proposed to explain the relative inefficacy of immunotherapies in myeloid leukemias, including expression of immunoregulatory molecules by malignant cells (2,3). The complexity of immune regulation in cancer is well illustrated by expression of interferon (IFN) related genes that influence the effectiveness of immunotherapy through immunostimulatory and immunosuppressive effects (4)(5)(6). Chronic interferon signaling has been associated with T-cell exhaustion in the setting of viral infections, but its role in cancer is controversial (7)(8)(9). Blastic plasmacytoid dendritic cell neoplasm (BPDCN) is an aggressive form of acute leukemia with few effective therapies that provides a unique opportunity to study IFN dysregulation in cancer. BPDCN is characterized by uncontrolled proliferation of transformed plasmacytoid dendritic cells (pDCs), specialized immune cells that link the innate and adaptive immune systems through the secretion of Type I interferons, including IFNA, particularly during viral infection (10). Studies of BPDCN up to this point have largely focused on the malignant pDC-like tumor cells, but few have focused on the T-cell response. An effective immune response relies on the interaction between healthy innate and adaptive immune systems, which is critical for immunotherapy and may be impacted by IFNA-producing pDCs.
Single-cell RNA-sequencing (scRNA-seq) has provided granular insights into the dynamics and phenotypes of immune cells. Specifically, scRNA-seq has been utilized to define subsets of exhausted T-cells in viral infection and cancer, including those that drive responses to ICB (11)(12)(13). Advances in other single-cell technologiessuch as single-cell methylation, chromatin accessibility, and mutation statushave allowed for multimodal data collection all from the same cell (14)(15)(16)(17). This type of analysis allows us to understand relationships between biological processes and heterogeneous cell types that cannot be studied using unimodal measurements alone. To study T-cell biology, scRNA-seq can be paired with sequencing of the T-cell receptor (TCR) aand b-chain variable regions. These regions can be enriched from single-cell transcriptomes while maintaining cell barcodes to integrate TCR sequencing data with gene expression profiles (18)(19)(20). Recently, TCR sequencing has been applied to solid and blood cancers to tie cell phenotypes to TCR properties, which has allowed for a deeper examination of T-cell subsets via multiple modalities (20)(21)(22). These studies also investigate the phenomenon of T-cell exhaustion, which occurs when cells enter a dysfunctional state in which they lose effector functions, such as proliferation and cytotoxicity, and gain immunoregulatory functions (7,23,24). Known drivers of Tcell exhaustion include continuous antigenic stimulation and chronic inflammation (7). Understanding the mechanism by which cells become exhausted is particularly important in the context of immunotherapy, as exhausted T-cells can regain effector functions via ICB therapy (24,25).
Here, we investigate the immune environment in five BPDCN bone marrow samples and five healthy controls using scRNA-seq and TCR sequencing. We used computational integration tools to create a unified healthy reference and applied multiple classification algorithms to annotate cell types in the BPDCN samples by consensus. These data show variable cell type expansions and a wide range of T-cell proportions between patients. Though we found heterogeneity in cell type proportions, we identified common signatures of upregulated IFNA response, downregulated TNFA signaling, and significantly increased CD8+ T-cell exhaustion. We applied a TCR sequencing approach termed T-cell Receptor Enrichment to linK clonotypes (TREK-seq) that we recently developed to identify expanded T-cell clonotypes in BPDCN samples and showed a correlation between CD8+ T-cell clone size and exhaustion scores. By providing a comprehensive map of cellular heterogeneity and T-cell transcriptomes/clonotypes in BPDCN patients, we lay the foundation for future development and evaluation of immunotherapies in this devastating cancer.

Patient Samples
Bone marrow aspirate was collected from 5 patients with blastic plasmacytoid dendritic cell neoplasm (BPDCN) who consented to an excess sample banking and sequencing protocol that covered all study procedures and was approved by the Institutional Review Board (IRB) of the Dana-Farber/Harvard Cancer Consortium. Additionally, 5 samples from healthy donors were collected for use as a control in this study, following the same approved protocol or a protocol approved by the IRB of Mass General Brigham. Blast quantification for each diagnostic patient sample was measured using four distinct methods. Aspirates were tested by Giemsa stain for cell morphology, flow cytometric measurement of blast markers (CD45 lo SSC lo/mid , abnormal CD4/CD123/CD56 expression), and Rapid Heme Panel for targeted sequencing of leukemia-associated genes; a bone marrow core biopsy was tested by hematoxylin and eosin (H&E) stain and immunohistochemical (IHC) histological analysis. Details for each of these samples are located in Supplemental Table 1.

Dataset Processing
Raw sequencing data were processed using CellRanger software (27) to generate FASTQ files and count matrices. Ambient RNAs were estimated and removed from the datasets using SoupX with default parameters (28). Each dataset was filtered to retain cells with >= 1000 UMIs, >=500 genes expressed, and <30% of the reads mapping to the mitochondrial genome. Three of the healthy control bone marrow samples were originally integrated using the IntegrateData() function in Seurat v4.0.3 and clustered using the same software at the default resolution of 0.5 (29). Cell cycle genes ("cc.genes" within the Seurat software) were removed from the integration anchors to combine smaller cell cycle-driven clusters and to eliminate redundancy, though these genes were retained in the final dataset for further analyses. Further, a combination of high resolution sub-clustering with Seurat and analysis of key T-cell gene expression were used to improve granularity in the naive T-cell compartment. Two additional healthy samples from older individuals were further integrated into the original controls using Seurat's TransferData () function. The BPDCN samples were processed and clustered individually using Louvain clustering within Seurat with default options and a resolution of 0.5, though these clusterings were not used in downstream analyses.

Cell Classification
Cells in each BPDCN sample were classified with four classification methods using the integrated healthy control samples as a reference: random forest (30,31), cellHarmony (32), Seurat's TransferData() function (33), and scPred with default parameters (34). Reference input for three of the algorithms was the integrated Seurat object and feature selection was performed separately as a part of each method. cellHarmony instead used a gene expression matrix of the cells by top 50 marker genes for each cluster (as defined by Seurat), the full expression matrix, and a table to cell type classifications for the reference. Given highly consistent classifications between the four methods, we selected cellHarmony-defined cell type labels for all downstream analyses. Classified cells from BPDCN patients were projected into the UMAP space of the integrated controls using the MapQuery() function in Seurat v4.0.3.

Gene Set Enrichment Analysis (GSEA)
Gene set enrichment was performed separately on cells from each cell type and BPDCN sample relative to the healthy controls. Log fold change values for every gene were extracted from the cellHarmony output (~/input/cellHarmony/OtherFiles), sorted in decreasing order, and used as input for a custom GSEA function (35). This function uses two GSEA implementations in R, gage (36) and fgsea (37) and only reports the results that are significantly enriched with both methods. Hallmark gene sets from the GSEA website were used for this analysis (https://www.gseamsigdb.org/gsea/msigdb/collections.jsp). Pathways that were significantly enriched in at least one comparison were plotted using the Pheatmap R package.

Pathway Score Quantification
To score individual cells for gene signatures, we combined all 52,803 cells from the dataset into one Seurat object, read in a curated list of signatures [IFNA response from the GSEA Hallmark gene set, TNFA signaling by NFKB from the Hallmark gene set, T-cell exhaustion from Penter et al. (21)] and applied the function AddModuleScore to the Seurat object using signature genes as features. Statistical significance for each cell type was calculated by comparing the median scores of normal samples (n = 5) to the median scores of BPDCN samples (n = 5) using the R function wilcox.test and default parameters. The number of biological replicates was less than five if a cell type was not detected in one of the donors (for example, n = 4 for BPDCN CD8+ Memory T-cells).

T-Cell Receptor Sequencing
T-cell receptor sequencing was performed on both control and BPDCN bone marrow samples using a modified protocol originally developed for Seq-Well (18) that we term TREK-seq. The modifications to the original protocol are as follows: in the TCR enrichment master mix, we added PartialRead1 and PartialTSO primers at a final concentration of 1.25 µM each. For amplification of TCR transcripts following enrichment, we used the same primers at a final concentration of 0.4 µM each. For the final PCR, we used UPS2-N70x and 10X_SI-PCR_P5 primers at a final concentration of 0.2 µM each to add the Illumina P5 and P7 sequences. The libraries were sequenced using a 150 cycle kit on the Illumina MiSeq and loaded at a final concentration of 10 pM. 28 cycles were used for Read 1, which read the cell barcode and UMI. 150 cycles were used for Index 1, which read the TCR region. TCRa and TCRb-specific custom sequencing primers were used for index 1 at a final concentration of 2.5 µM. We aimed for a MiSeq cluster density of roughly 450k/mm2. Primer sequences can be found in Supplemental Table 2.

WAT3R Computational Pipeline
To analyze TREK-seq data, raw sequencing data were demultiplexed with bcl2fastq (v2.20.0), and the resulting FASTQ files were reformatted to join the cell barcode, UMI and corresponding TCR sequence in the same file. For this analysis, we developed and applied the bioinformatics pipeline WAT3R, or Workflow for Association of T-cell receptors from 3' single-cell RNA-seq (WAT3R) (https://github.com/mainciburu/ WAT3R) (38). Briefly, we first corrected cellular barcodes allowing one mismatch with the 10x Single Cell 3' v3 whitelist. UMI correction was also performed by clustering together UMIs with one mismatch. Every corrected barcode and UMI sequence was then added to the corresponding FASTQ read header. Next, we applied a quality filter to remove every read with an average score < 25. To account for barcode swapping, TCR sequences with identical barcode and UMI were subjected to clustering using an identity threshold of 0.9. The subsequent analysis was carried out exclusively with clusters representing at least 50% of the reads with identical barcodes and UMI and doubling the number of reads from the second most abundant cluster. Next, a consensus sequence was built to summarize the TCR sequences in each of those clusters. We required a minimum of 3 sequences and allowed for a maximum error rate of 0.5 and gap frequency of 0.5 per position. Consensus sequences were aligned to the VDJ segments reference available at the IMGT database. IgBLAST with the default parameters was used for the alignment. Finally, the V, D and J calls and CDR3 sequence with higher UMI counts, for both TRA and TRB genes, were assigned to each cell barcode. For downstream analysis, only TRB variable regions were used.

Statistics
P-values less than or equal to 0.05 were considered statistically significant. Cell population proportions between younger and older healthy donors were assessed for significant differences using a chi-squared test. For the GSEA heatmap ( Figure 3A), only significantly enriched pathways (P-value ≤ 0.05) in at least one sample were visualized. For the heatmaps of genes in the IFN Alpha Response and TNFA Signaling via NFKB hallmark gene sets ( Figure 4), P-values were calculated via t-test for each comparison of reference (BM1-5) to each BPDCN sample per cell type and adjusted using the Bonferroni method. Genes that were significant (P-value ≤ 0.05) following correction and had an expression value of 0.5 or greater (log normalized expression) in at least one BPDCN or reference sample were retained for the heatmap. A Wilcoxon signed-rank test was used for pathway score and exhaustion score quantification, a Kruskal-Wallis rank sum test was used to compare normalized T-cell clonotype sizes, and a Spearman ranked correlation analysis was used to test the relationship between clonotype size and exhaustion scores.

Data and Software Availability
All data used in this paper are publicly available in the Gene Expression Omnibus (GEO), accession number GSE189431. The scripts used to generate the results in this paper can be found at the project GitHub repository: https://github.com/ EDePasquale/BPDCN.

Identification of Cell Populations in Healthy Bone Marrow
To map the cellular diversity in healthy bone marrow samples, we performed scRNA-seq using the 10x Genomics platform (27). We profiled cells from three younger healthy donors (age 31-45) and two older donors (age 74 and 75) (Supplemental Table 1). The two older donors were included as age-appropriate controls for BPDCN, which has been reported to have a median age of 65 (39). Unique cellular barcodes were used to assign transcripts to cells and individual mRNA molecules were quantified using Unique Molecular Identifiers (UMI). Across all five healthy donors, we retained 25,726 cells following quality control and filtering.
Cells from these healthy donors (BM 1-5) were integrated using the Seurat package (IntegrateData and TransferData functions) to remove batch effects and generate a unified UMAP projection of the data ( Figure 1A and Supplemental Figure 1). We initially identified 16 unique clusters based on gene expression that were interrogated for expression of canonical marker genes and gene signatures of blood cell types ( Figure 1B). Additional sub clustering was performed in the lymphocytic compartment, and new T-cell classifications were made based on expression of CD8A, CD4, ITGB1, and CCR7 gene expression. Specifically, CD4+ T-cells and CD8+ T-cells could both be split into Naive (ITGB1-CCR7+) and Memory (ITGB1+ CCR7-/+) subsets. The final reference was comprised of 17 clusters representing known hematopoietic cell types, including HSCs and progenitors, myeloid, erythroid, and lymphoid cells ( Figure 1C). All cell types were represented in each donor at broadly comparable levels, though variation exists as expected, particularly between the younger and older donors (P < 0.05, chi-squared test) ( Figure 1D) (40).

Single-Cell Profiling of BPDCN Bone Marrow
Approximately half of BPDCN patients present with skin tumors only, whereas the other half present with bone marrow involvement, i.e. detection of malignant blasts in the bone marrow by conventional clinical assays, including Giemsa stain, flow cytometry, H&E, and Rapid Heme Panel (41). These assays often indicate divergent bone marrow involvement due to inherent differences in the assays and due to sampling. For example, bone marrow aspirate samples, as used in this study, can be variably diluted with peripheral blood elements depending on the quality of the biopsy and the number of passes taken. Further, our single-cell sequencing pipeline involves density centrifugation, cryopreservation and microfluidics that could influence cell type proportions. A comparison of bone marrow involvement estimates is included in Supplemental Table 1. To understand the proportional and transcriptional changes that occur in the bone marrow in the presence of BPDCN tumor cells, we selected five patients with 1.5-97% disease involvement at diagnosis. We performed scRNA-seq on BPDCN bone marrow samples using the same protocol as in the healthy samples. From this sequencing, we acquired 27,077 cells across all BPDCN donors following quality control and filtering.
We next wanted to classify cell types in the five BPDCN samples (BPDCN 1-5) using the healthy controls cells as a reference. To this end, we tested four computational classification methods that utilize different underlying algorithms to assign cell type labels to cells: random forest, cellHarmony, Seurat TransferData, and scPred. The resultant cell type predictions were concordant across all four methods in most cell types, particularly in the erythroid and lymphoid lineages, lending confidence to the classification methods ( Figure 2A). Classifications from the cellHarmony algorithm were chosen as the cell type labels for downstream analyses (Supplemental Table 3). Some variability was noted in the proportion of cells labeled as pDC, conventional dendritic cells (cDC), and Pro-B cells (ProB) in BPDCN 4 and 5; we hypothesized that these could be malignant cells that align poorly to any healthy population, which would be consistent with the clinical annotation of high malignant cell content in these bone marrows (Supplemental Table 1). To visualize shifts in gene expression programs and differential cell type proportions, we projected the cells from each BPDCN patient onto the same UMAP space as the healthy reference, colored by cellHarmony labels ( Figure 2B). Indeed, we observed that the pDC and cDC labeled cells in BPDCN 4 and 5 were shifted toward the HSC/Prog and ProB populations, suggesting that the malignant tumor cells in these samples were transcriptionally different from their healthy counterparts. For this study, we were primarily interested in the role of T-cells in response to and as a potential treatment for BPDCN, therefore we next focused on alterations in cells of the lymphoid lineage.

T-Cell Populations Show Differential Regulation of IFNA and TNFA Pathways
We used Gene Set Enrichment Analysis (GSEA) to identify pathways that were significantly enriched in the BPDCN samples. Using the Hallmark gene sets and gene expression changes derived via cellHarmony, we generated a heatmap containing all up-or down-regulated pathways enriched in any BPDCN sample compared to the healthy controls ( Figure 3A).

Memory T-cells, 3/5 in CD4+ Memory T-cells, 3/3 in CD8+
Naive T-cells, and 3/3 in CD4+ Naive T-cells, along with 3/5 samples in NK cells. IFNA is a marker of immune activation and pDCs produce IFNA, furthering our interest in this pathway. Using a different statistical framework, we next scored all individual CD8+ Memory T-cells for their expression of genes in the IFNA response gene set. We found a significant increase in BPDCN 1-4 relative to control samples, confirming increased IFNA response gene expression in BPDCN T-cells (P = 0.0159, Wilcoxon signed rank test, Figure 3B). BPDCN 5 did not have cells classified as CD8+ Memory T-cells and was not tested. Increased IFNA response gene scores were also observed in CD4+ Memory T-cells and NK cells (data not shown). We found that tumor necrosis factor alpha (TNFA) signaling was downregulated in many of the T-cell and NK populations: 4/4 CD8+ Memory T-cell, 3/5 CD4+ Memory T-cell, and 2/3 CD4+ Naive T-cell, and 3/5 of the NK cell samples with the exception of CD8+ Naive T-cells (0/3) ( Figure 3A). Scoring individual CD8+ Memory T-cells for the TNFA gene signature also showed significant reductions in BPDCN (P = 0.0159, Wilcoxon signed rank test, Figure 3C). These results demonstrate that T/NK-cells in BPDCN exhibit consistent gene expression changes compared to healthy control cells, including increased response to the pDCrelated cytokine IFNA and decreased TNFA signaling. We next investigated which genes in the IFNA response gene set were driving the enrichment in BPDCN T-cells. We generated a heatmap of genes that were significantly different from control expression levels in at least one sample/cell type pair and further filtered to those genes with normalized expression values of at least 0.5 in one pair (P ≤ 0.05, t-test, Bonferroni adjusted, Figure 4A). We also generated heatmaps with the full list of genes in each gene set (Supplemental Figure 2). Related to the IFNA response, genes that were upregulated in BPDCN compared to control T-cells in the majority of samples included: IRF1, PSMB9, PSME2, CD74, ISG15, ISG20, CD47, LY6E, PSMB8, and SELL ( Figure 4A). The increased expression of IRF1, PSMB9, and PSME2 were notable in all 5 BPDCN samples across lymphoid cell types, with the exception of BPDCN 4's CD4+ Naive T-cells due to low cell numbers. In contrast, CD74, ISG15, ISG20, CD47, LY6E, PSMB8, and SELL were most notably increased from control levels in CD8+ Memory T-cells and NK cells, and in BPDCN 1-3 specifically. These genes are involved in T-cell proliferation/ differentiation, antigen presentation, and inflammation, consistent with the view that IFNA signaling plays a critical role in immune activation. Related to TNFA signaling, genes that were downregulated in BPDCN T and NK cells compared to control samples are as follows: NFKBIA, TNFAIP3, DUSP1, DUSP2, FOS, JUN, ZFP36, BTG2, and CD69 ( Figure 4B). These genes were downregulated relative to controls in all T-cell and NK cell subsets in BPDCN 1-4, with the largest expression differences being isolated to the CD4+ and CD8+ Memory T-cell and NK cell clusters. These genes are involved in immune regulation via inflammatory response, proliferation and T-cell activation. Decreased expression of TNFA genes is indicative of a lowered immune response, particularly in relation to T-cell activation. Overall, our results suggest an altered balance between inflammatory pathways in T-cells of BPDCN patients, shifting towards IFNA at the expense of TNFA signaling (42). The changes we observe in IFNA and TNFA target genes are likely to affect the function of the adaptive immune system in BPDCN.

Expanded T-Cell Receptor Clonotypes in BPDCN Correlate With Signatures of Exhaustion
T-cell antigen specificity is determined by the TCR sequence, which consists of the aand b-chains encoded by TRAV and TRBV genes that are recombined during T-cell development. We reasoned that sequencing the variable region of the TCR may help to elucidate the biology of T-cells in BPDCN. We performed TCR sequencing on healthy control and BPDCN samples using TREK-seq, a protocol for paired transcript and TCR capture in 10x Genomics data (see Material and Methods). We developed and applied the Workflow for Association of T-cell receptors from 3' single-cell RNA-seq (WAT3R) computational pipeline to detect groups of cells expressing the same TCR sequence, i.e. Tcell clonotypes (38). We identified 10,870 T-cell clonotypes in our dataset, each supported by multiple sequencing reads. Using this pipeline, we found that TCR sequences were mainly detected in T-cell types in healthy controls, which is in line with the reported specificity highlighted in the WAT3R publication and lends confidence to TREK-seq protocol ( Figures 5A-E). In BPDCN, we mapped most of the TCR sequences to T-cell types ( Figures 5F-I), though more variability is seen in these samples than the healthy controls. Specifically, in BM 1-5, 94.7% of TCR sequences were detected in T-cells ( Figure 6A), and we observed clonal expansions in the CD8+ Memory T-cell compartment, which was expected ( Figure 6B). In BPDCN 1-4, 86.7% of TCR sequences mapped to T-cell types, though 3.9% mapped to NK cells and 7.8% mapped to likely tumor populations (pDC, cDC, and ProB). NK cells with TCR results may be a misclassified subset of NKT cells not prevalent enough in the healthy controls to warrant a separate cluster during classification. We observed TCR transcripts in 12.6% of pDC/cDC/ProB cells in BPDCN 4. This finding is of uncertain significance in this biological context as BPDCN blasts do not typically express the cell-surface proteins (e.g., CD3) associated with functional TCR complexes (43). TCR sequences from BPDCN 5 could not be successfully retrieved with the TREK-seq, likely due to the low number of T-cells in the sample, and it was therefore excluded from all downstream analyses.
Next, we ranked clonotypes in each sample by their normalized size (number of cells with the same TCR over all cells in which a TCR was detected). For BM 1-5, we found that none of the control samples were dominated by a single clone or handful of clones, but that the distribution of clonotypes was largely uniform (Supplemental Figures 3A-E). BPDCN samples showed heterogeneity in the distribution of clonotypes: while BPDCN 2-4 do not exhibit a single prevalent clone, similar to the healthy controls, we detected a clone that comprises 30.9% of all TCRs in BPDCN 1 (Supplemental Figures 3F-I). In this sample, one-third of the cells within the largest clone are labeled as NK cells (Supplemental Figure 3F), leaving open the possibility that the NK-labeled population with TCR results are NKT cells. We next evaluated T-cell exhaustion by scoring each T-cell for the expression of exhaustion-associated genes (21). We overlaid cell type label, exhaustion scores, and normalized clone size parameters onto UMAP plots for each BPDCN sample containing T-cells (Figure 7). CD8+ Memory T-cell clusters were enriched for exhausted cells and expanded clonotypes in all patients, which was confirmed by quantification of signature scores and clone sizes in cell types ( Figure 8A). Overall, exhaustion scores in CD8+ Memory Tcells in BPDCN 1-4 were significantly higher than controls (P = 0.0159, Wilcoxon signed-rank test), though the range of scores among individual cells within each sample was high ( Figure 8B). These results are consistent with the pathway analysis above, as IFNA-mediated immune stimulation and low levels of TNFA signaling were previously associated with T-cell exhaustion (7,8,44).
To explore the relationship between exhaustion and clonotype sizes, we correlated these two parameters at the single-cell level. We found that CD8+ Memory T-cell clonotypes of various sizes can express T-cell exhaustion signature genes, including some clonotypes that make up a small proportion of the total T-cell population (Figures 8C-F). Spearman ranked correlation analysis revealed a positive relationship between clonotype size and exhaustion signature expression: rho of 0.67, 0.61, 0.63, and 0.4 for BPDCN 1-4, respectively (Supplemental Figure 4). These results suggest that clonal expansion and exhaustion of CD8+ Memory T-cell clonotypes may be functionally linked in the BPDCN tumor environment.

DISCUSSION
To study the immune environment in BPDCN, we utilized paired single-cell sequencing of the transcriptome and T-cell receptors in healthy controls and patient samples. We identified 17 transcriptionally unique cell clusters in an integrated dataset of bone marrow from five younger and older healthy donors. We used this dataset to label cell types in five BPDCN patient bone marrow samples using four distinct classification algorithms, revealing heterogeneity in cell type proportions that correspond to clinical assessment of bone marrow involvement. An analysis of differentially regulated gene sets revealed that IFNA response genes were significantly increased in T-cells of the BPDCN samples and TNFA signaling genes were decreased. A deeper examination of expression within each gene set identified a subset of genes driving this enrichment, with the associated functions of T-cell proliferation, activation, and antigen presentation. TREK-seq showed T-cell specific TCR expression and uniformly distributed clone sizes in healthy controls; both of these results were more variable in BPDCN patient samples. Further, we saw increased T-cell exhaustion in CD8+ Memory Tcells that positively correlated with increasing clone sizes. These results demonstrate the utility of single-cell multiomics by establishing a resource of the BPDCN tumor environment that provides new insights relevant to immuno-oncology. Our results contribute to the body of research surrounding the role of the adaptive immune response in cancer and the relationship between T-cell receptor clonality and exhaustion. In a healthy immune system, activated pDCs produce large quantities of IFNA to stimulate activation of T-cells (10). In this study, we observe increased expression of IFNA response genes in T-cells relative to T-cells in the healthy controls, initially suggesting higher levels of IFNA production by pDC-like tumor cells. However, previous research on BPDCN has shown that pDC-like tumor cells produce lower levels of IFNA relative to pDCs from healthy individuals (45,46). Further, little research has been reported describing the role of TNFA signaling in BPDCN, though studies have shown increased TNFA signaling in the plasma of patients with acute myeloid leukemia (47,48). While the directionality of differential expression in these gene sets in T-cells is in contrast to what has been observed in the adjacent BPDCN cells, it is consistent with our findings of increased T-cell exhaustion. We hypothesize that, despite the  (BPDCN 1-4), separated by column. UMAPs are colored by cluster (first row); T-cell exhaustion score, with red indicating high exhaustion and gray indicating low exhaustion (second row); and normalized clone size, percentage of all cells in the dataset that share the same TCR sequence, with blue indicating high clone size and gray indicating low clone size (third row).
tumor cells exhibiting lower IFNA production at the individual cell level, abnormal accumulation of pDC-like tumor cells in BPDCN may lead to increased IFNA production and chronic Tcell activation, eventually leading to T-cell exhaustion and consequent TNFA downregulation.
Exhaustion of cytotoxic T-cells has been a major hurdle in establishing effective immunotherapy treatments for blood cancers (24,25,49). Reactivation of T-cell effector function in exhausted T-cells is the goal of ICB, but aberrant IFNA and IFNG signaling have been shown to inhibit this process (50,51). While the prevalence of exhaustion has been studied in other cancers including acute myeloid leukemia, its role in BPDCN remains unknown. In this study, we establish that BPDCN patients have heightened T-cell exhaustion relative to controls. We also find a positive association between T-cell exhaustion and T-cell clonal expansion, consistent with other cancers (21,52). Our findings highlight a potential mechanism of T-cell exhaustion via persistent IFNA signaling that might be targeted to restore anti-tumor immunity.
For this study, we selected five BPDCN patients with varying levels of clinical bone marrow involvement, though more samples will be needed to ensure that the full spectrum of patient heterogeneity is captured. The quantification of malignant blasts in the bone marrow is an important clinical feature at the time of BPDCN diagnosis. Independent bone marrow aspirates are evaluated for abnormal cell morphology (Giemsa stain), surface phenotypes (flow cytometry), and genetics (targeted sequencing, karyotyping). This can be supplemented by H&E histological analysis of a bone marrow core biopsy. Significant differences between these assays illustrate the heterogeneity and complexity of BPDCN. Across the samples we analyzed, the proportion of dendritic cells by scRNA-seq is generally in agreement with the overall clinical evaluation of malignant blasts, suggesting that single-cell sequencing may provide additional support when defining tumor involvement. Future studies should evaluate limitations introduced by sampling bias and cell processing; further opportunities could come from integration of genetic and phenotypic analysis with single-cell gene expression. Multimodal single-cell technologies enable investigation of gene regulation at a high resolution through multiple modalities representing interacting processes within a cell. In addition to our measurements of gene expression and TCR sequences, other multimodal technologies could elucidate mechanisms of regulation in BPDCN. For example, the application of singlecell genotyping methods in future work would help to classify and characterize the transcriptomes of malignant tumor cells (2,16). Inclusion of additional patient samples and complementary single-cell measurements would strengthen our initial findings and uncover new results that further illuminate T-cell biology in BPDCN.
In summary, we apply scRNA-seq and TCR sequencing with computational techniques to catalog the cellular heterogeneity in BPDCN patient samples. Multiomics technology allows us to gain a deeper understanding of immune cell dynamics by assessing the diversity of immune cell states via scRNA-seq and the expansion of T-cell clonotypes through TREK-seq. Our results suggest that the balance between IFNA and TNFA signaling is disrupted in BPDCN, potentially leading to T-cell clonotype expansion and exhaustion. The discovery of mechanisms by which BPDCN cells evade immune destruction will lead to the development of new cancer therapies that leverage tumor-reactive T-cells.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: GEO database, accession number GSE189431.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by Institutional Review Board (IRB) of the Dana-Farber/Harvard Cancer Consortium and IRB of Mass General Brigham. The patients/participants provided their written informed consent to participate in this study.