Single cell profiling of γδ hepatosplenic T-cell lymphoma unravels tumor cell heterogeneity associated with disease progression

Clinical characteristics of the studied patient with HSTCL

The detailed clinical characteristics of the recruited patient with HSTCL were summarized in Table 1. This 17-years old female patient was diagnosed γδ HSTCL on May 18 2021 (Fig. 1A), initially presented with fever and fatigue, and subsequently developed severe hemorrhage due to thrombocytopenia (Table 1). CT scanning revealed splenomegaly and hepatomegaly (estimated spleen volume: 910.56 cm3 and liver volume: 1617.21 cm3 Fig. 1B left), while lymphadenopathy was absent (Fig. S1). Neoplastic lymphocytes were identified from bone marrow aspirate (Fig. 1C), which were large, round or irregular cells, with abundant cytoplasm but no or few particles in cytoplasm. The nuclei were round or irregular, depressed and distorted, with dense chromatin and 1–3 inconspicuous nucleoli, PAS glycogen staining was positive. The abnormal cells accounted for ~ 50% of the lymphocyte population. By flow cytometry, the neoplastic lymphocytes were found to be CD2 + , CD3 + , CD16 + , CD45 + , TCRγ/δ + , CD7P + , CD4-, CD8- and TCR α/β-, which were the typical immunophenotypes of γδ HSTCL [4, 6] Fig. 1D). TCR rearrangement assay further confirmed clonal rearrangement of TCR β (TCRB) and γ (TCRG) genes.

Table 1 The clinical characteristics of the studied patientFig. 1figure 1

Clinical presentations of the studied HSTCL patient. (A) Timeline of the patient diagnosis, chemotherapy administration and biopsies collected for single cell sequencing. (B) CT scan and 3D image reconstruction during the disease progression, the volume of the spleen and liver was measured and labeled. (C) Wright-staining of bone marrow aspirate, arrows indicate selected typical neoplastic T cells. (D) Immunophenotypes of the neoplastic γδ T cells identified by flow cytometry

After diagnosis, the patient began to receive chemotherapy hyper-CVAD-A regimen as the induction treatment (Fig. 1A, Table 2). During the first course of chemotherapy, the patient’s body temperature returned to the normal and the overall conditions were relatively stable. Just 17 days after the first chemotherapy completed, the patient was hospitalized again due to fever and thrombocytopenia. The patient was then received ICE regimen (Fig. 1A Table 2), which achieved similar therapeutic outcomes as hyper-CVAD-A. Separated by 32 days, another course of ICE regimen was administrated to the patient (Fig. 1A Table 2). During each chemotherapy period, the indexes of routine blood examinations such as red blood cell count, white blood cell count, platelet count could gradually returned to the normal range. However, the symptoms relapsed soon after the treatments were completed. During the second course of chemotherapy, CT examination revealed significantly reduced spleen volume (524.05 cm3) while liver volume stayed similar size (1613.09 cm3, Fig. 1B third column), but after the third course of treatment, the spleen and liver enlarged again to 1001.58 cm3 and 1645.19 cm3 (Fig. 1B. right). These results suggested HSTCL was very aggressive and refractory.

Table 2 Chemotherapy administrations for the patientEstablish a single cell atlas for HSTCL

To understand the cellular and molecular characteristics of HSTCL and how they can be changed during chemotherapies, we performed paired scRNA-seq and scTCR-seq for 4 samples collected from the patient. 2 pre-treatment samples were collected from bone marrow (BM1) and PBMC (PBMC1), respectively, and 2 were collected post-treatment (BM2, PBMC2) (Fig. 1A). After quality control and data pre-processing, 36,092 cells (8,295 BM1, 9,189 PBMC1, 5,730 BM2 and 12,878 PBMC2) were retained for downstream analysis. Followed by dimension reduction and Leiden clustering, the cells were visualized by uniform manifold approximation (UMAP) (Fig. 2A). 7 major clusters were identified from the cell population, and the cell type identities were assigned to each cluster by canonical markers from previous publication (Fig. 2B): T and NK cells (T_NK, marked by CD3D, CD3E, NKG7 and GNLY), B cells (B, marked by CD79A and MS4A1), Plasma cells (Plasma, marked by MZB1 and IGKC), Monocytes (Monocyte1 and Monocyte2, marked by CD14, LYZ and CST3), Erythroid and progenitor cells (Ery_Pro, marked by HBD and CA2) and malignant γδ T cells (GD_Tumor, marked by TRCD, TRDV1 and CD3D). We didn’t observe strong batch effects among different samples (Fig. S2). To further confirm the identity of the tumor cluster, we evaluated the expression of HSTCL markers previously identified by immunohistochemistry and flow cytometry. We confirmed that these GD_Tumor cells were CD2 + , CD3 + , CD4-, CD5-, CD7 + and CD8, CD38, CD56 (NCAM1) partial positive, aligned with the diagnostic guidelines (Fig. S3) [4, 6].

Fig. 2figure 2

A cell type atlas of γδ hepatosplenic T-cell lymphoma. (A) Cell type atlas of γδ HSTCL presented by UMAP. (B) Dotplot of the canonical marker genes used to annotate major cell types. Dot size represents % of cells of that cluster expressing the given gene, while color indicates the expression level of that cluster. (C) Barplot to show the cell type compositions in each sequenced sample. Significant composition changes are observed pre- and post-treatment. (D-G) Cell type % or absolute count during treatments. Dotted lines indicate chemotherapy-initiated date. (D) Tumor cell % estimated from histology slides (Wright staining) at different time points. Tumor cell % is sharply reduced during the treatments. (E) Absolute monocyte counts flatulate during the disease progression, with a drop after chemotherapy administration and soon increase to an unusual high level. (F) In clinical bone marrow examinations, erythroid lineage cell % (polychromatophilic (Poly-E) and orthochromatophilic (Ortho-E)) shows an increasing trend. (G) In lymphocyte subtype counting, B cell absolute counts drop to very low level after treatment

We observed unique cell type composition distribution in HSTCL and substantial composition changes post-treatment (Fig. 2C). Tumor cells were greatly reduced post-treatment in both PBMC (pre: 6.1% vs post: 1.1%) and BM (pre: 45.2% vs post: 18.5%), In Wright staining, tumor cell was also significantly decreased post-treatment (Fig. 2D). Unusual high percentage of monocytes was observed in PBMC samples, and became even higher post-treatment (pre: 40.2% and post: 65.0%, Fig. 2C), which was consistent with clinical blood examinations (Fig. 2E). Previous studies suggested that elevated monocyte count may associate with poor prognosis and identify high-risk patients in lymphomas [15, 16]. Ery_pro cell % were increased post-treatment in both BM (pre:14.2% vs post: 51.1%) and PBMC (pre:1.6% vs post: 5.0%). In Ery_pro cluster, most cells were erythroid and erythroid progenitor cells (Fig. S4). In clinical bone marrow examination, erythroid lineage cell % also showed an increasing trend (Fig. 2F). B cells reduced to a very low percentage post-treatment in both BM (pre: 5.2% vs post: 0.8%) and PBMC (pre: 15.7% vs post: 0.7%). In lymphocyte subtype counting by flow cytometry, B cells were significantly dropped after chemotherapy (Fig. 2G). These results suggested that the tumor microenvironment (TME) was dynamically reshaping as responses to the disease progression and treatment.

HSTCL malignant γδ T cells had unique clonal and molecular features

To investigate the clonal expansion of the malignant γδ T cells, we performed scTCR-seq in parallel with scRNA-seq. Among the 5,515 cells in the GD_Tumor cluster, 4,505 (81.7%) cells had recoverable γδ TCR information. The TCR repertoire of malignant cells was almost founded by a single clonotype, suggesting the tumor cells likely arose from a single ancestor cell and clonal expanded subsequently (Fig. 3A). For the 4,419 cells had assembled γ chains, 4,374 (99.0%) had the same γ chain sequence (TRGV3 + TRGJ1 + CATWDNLMENYYKKLF), for the 3,614 cells had assembled δ chains, 3,576 (99.0%) had the same δ chain sequences (TRDV1 + TRDJ1 + CALAVVVGGLGTDKLIF), and for the 3,494 cells had both γ chain and δ chain assembled, 3,445 cells had the above γ chain and δ chains sequences (Fig. 3A blue dots). As a comparison, we identified 559 cells in NK_T cluster with γδ TCR, suggesting they were normal γδ T cells from this patient. These normal γδ T cells had a more diverse TCR repertoire, ~ 60% of the cells were singletons (Fig. 3B).

Fig. 3figure 3

Molecular features of γδ hepatosplenic T-cell lymphoma. (A) TCR sequencing reveals the distribution of a single prominent clonotype in the malignant γδ T cells (Blue dots). (B) TCR sequencing shows that the normal γδ T cells have a diverse TCR repertoire, colors indicate different clonotypes. (C) Volcano plot to show the differentially expressed genes between malignant and normal γδ T cells, top 10 up and down genes by fold changes are labeled. (D) Heatmap to present the top 20 differentially expressed genes between malignant and normal γδ T cells. (E) KEGG pathway enrichment of the differentially expressed genes. Dot size represents the number of significantly changed genes in the pathway and color indicates adjusted p values

To give a more comprehensive description of the malignant γδ T cells and search for novel candidate genes relevant to HSTCL pathogenesis, we compared the malignant γδ T cells gene expressions with those of normal γδ T cells. We found a substantial number of genes were differentially expressed between the malignant and normal γδ T cells (Fig. 3C, Table S1). Among them, the top up-regulated genes in tumor cells were those related to cancer progression (eg, AREG, CCL3, PLEKHA5, VCAM1 and MCM7, Fig. 3D), which might be the targets of innovative treatments. Normal γδ T cells have strong cytotoxicity by expressing NK cell-associated genes like NKG2D, GNLY and a profile of cytotoxic molecules (eg, GZMA, GZMB, GZMM, and PRF1) [17,18,19]. However, these genes were mostly absent in the malignant γδ T cells (Fig. S5). On the contrary, killer cell immunoglobulin-like receptor 2DL2 (KIR2DL2) and KIR2DL3 were increased in the malignant γδ T cells, which encode the inhibitory receptors of the KIR family and reduce the cells’ cytotoxic function [20]. These results suggested that malignant γδ T cells differed from normal γδ T cells by lacking the cytotoxic modules. KEGG pathway enrichment analysis further showed that the top 1000 tumor up-regulated genes were significantly enrichment in pathways involving in cell proliferation and cancer development, such as Cell cycle, DNA replication, p53 signaling pathways (Adjusted p value = 6.99e-14, 5.28e-11 and 1.03e-03, respectively. Figure 4E).

Fig. 4figure 4

scRNA-seq reveals tumor cell heterogeneities. (A) UMAP of malignant γδ T cells, which are separated into two distinct sub-clusters. (B) Barplot to show that Tumor_1 cells are mostly from pre-treatment samples (BM1 and PBMC1) while Tumor_2 cells are dominant in post-treatment samples (BM2 and PBMC2). (C) Cell cycle phase distribution among Tumor_1 and Tumor_2. Tumor_1 has a higher portion of cells in S and G2M phases. (D) Dotplot of proliferation markers expression in Tumor_1 and Tumor_2. (E) Volcano plot to show the differentially expressed genes between Tumor_2 and Tumor_1 cells, top 10 up and down genes by fold changes are labeled. (F) Violin plot of selected differentially expressed genes between Tumor_1 and Tumor_2. (G) Violin and box plot of Cancer Hallmarks scores. Tumor_1 has higher scores in most cancer hallmarks. (H) Heatmap of the inferred CNVs in Tumor cells using normal immune cells as references. Red indicates CNV gains and blue indicates CNV loss. Green arrows indicate chromosome 7 and 8 duplications inferred in majority of the tumor cells. Yellow arrows indicate tumor sub-cluster specific CNVs. (I) Selected tumor sub-cluster specific CNVs. Left: a chromosome 4 partial duplication is enriched in Tumor_2. Right: a chromosome 19 partial deletion is enriched in Tumor_2. Blue indicates cells with specified CNVs

Transcriptional landscape revealed the heterogeneity of HSTCL malignant cells

We next investigated the heterogeneity of the malignant cells at single-cell level. The 4505 cells from GD_Tumor cluster with γδ TCR were taken out and re-clustered. Two distinct tumor clusters were identified from these cells (Fig. 4A), which had a strong correlation with the treatment status (Fig. 4B). Tumor_1 cluster consisted 3439 (75.3%) tumor cells, 98.0% of which were from pre-treatment samples (BM1 and PBMC1) and were almost absent in post-treatment samples. Tumor_2 cluster was relative smaller, included 1066 cells but dominating the post-treatment samples (BM2 and PBMC2), suggesting the chemotherapies may have a selection effect on the tumor cell populations.

We further did a comprehensive comparison between the two tumor clusters. From the expression of G2M and S-phase genes, we inferred the cell cycle phases for each tumor cell and observed different phase distribution between the two clusters (Figs. 4C, S6). 76.0% of the Tumor_1 cells were at S and G2M phases, while only 60.1% in the Tumor_2 cluster. We then surveyed the expression of proliferation associated genes (eg, MKI67, PCLAF, PCNA), Tumor_1 cells also had higher expression than Tumor_2 among these genes (Fig. 4D), together with the cell cycle analysis, suggesting that Tumor_1 cells were more actively proliferating than Tumor_2 cells.

Differentially expressed genes were identified from the two clusters (Fig. 4E and F, Table S2). NK cell receptors (NKR) were frequently expressed in normal γδ T cells. We found the kill receptor KLRC1 (NKG2A) was expressed in 62.0% Tumor_1 cells but only 3.3% in Tumor_2.KLRK1 (NKG2D), another NKR highly expressed in normal γδ T cells, was almost absent in Tumor_1 but presented in 25.8% Tumor_2 cells. Unexpectedly, we observed high expression level of CD19, a marker for B cell, in Tumor_1. Other B cell marker such as CD79, CD20 were not identified in Tumor_1 (Fig. S7), suggesting the CD19 expression was not caused by doublet. Therefore, Tumor_1 had a CD3, CD19 double positive phenotype, which is rarely reported in previous studies [21, 22]. In Tumor_2, the tumor cells lost the CD19 expression. Among the genes highly expressed in Tumor_2, many were associated with tumor progression and potential drug resistance. IL32 was reported to accelerates the proliferation of cutaneous T-cell lymphoma cell lines through mitogen-activated protein kinase (MAPK) and NF-κB signaling pathways [23]. TOX2 is a novel tumor driver, which promotes in Natural Killer/T-Cell Lymphoma cell growth and enhances ability of colony formation, as well as protects cell viability under adverse condition [24]. AIF1 and AKAP12 were reported to correlated with various chemo resistance [25, 26]. Interestingly, we also found CD38 was highly expressed in Tumor_2 cells, which is a marker of multiple myeloma. These results suggested that Tumor_2 might obtained unique survival features under the therapeutic conditions.

To describe different phenotypic characteristics of the malignant cells, we calculated the 10 cancer hallmarks scores curated by Cancer Hallmark Genes (CHG) database [27]. By permutation test, we identified 9 out of the 10 hallmark scores were significantly different between the two tumor clusters (Fig. 4G, Table S3). Except “Inducing Angiogenesis”, all other 8 hallmarks were significantly higher in Tumor_1 cluster, indicating Tumor_1 had stronger cancer related signaling pathways. Since Tumor_2 cells were mostly from post-treatment samples, the decreased hallmarks scores could be a response to therapeutic pressure. Tumor_2 might be a less aggressive tumor subtype, however, may have unique features to tolerant the chemotherapy.

Chromosomal aberrations are frequently observed in HSTCL [28, 29]. To determine if the transcriptional heterogeneity was potentially associated with chromosomal abnormalities, we inferred large scale copy number variations (CNV) at single cell resolution using inferCNV R package. We identified gains of chromosome 7 and 8 in majority of the tumor cells (Fig. 4H green arrows), which are the most common chromosomal abnormalities in HSTCL, presented in up to 63% and 50% of cases, respectively, and frequently cooccurring [28, 29]. The two tumor clusters also showed cluster-specific CNVs. For example, a portion of cells in Tumor_2 cluster demonstrated a CNV gain in chromosome 4 and a loss in chromosome 19 (Fig. 4H yellow arrows and Fig. 4I). These findings suggested that some chromosomal aberrations might have occurred during the evolution of the tumor cells under therapeutic pressures and may contribute to the survival of the tumor cells.

Tumor microenvironment heterogeneity associated with disease progression

To characterize the features of TME in this HSTCL patient, we performed analysis for the immune cells in the single cell profiles, with a focus on T cells and monocytes. We first re-clustered the cells in T_NK cluster, and grouped them into 4 major sub-clusters. Based on the expression signatures, these clusters can be annotated as CD4 T cells (CD4_T), CD8 memory T cells (CD8_Mem_T), CD8 naïve T cells (CD8_Naive_T) and NK cells (NK) (Figs. 5A and B, S8). In the post-treatment samples, the CD8 memory T cell % were increased (BM1 21.9% vs BM2 40.6% and PBMC1 21.7% vs PBMC2 35.3%, Fig. 5C), indicating these cells populated during disease progression. scTCR-seq results showed that the cells in CD4_T and CD8_Naive_T were mostly singletons (71.5% and 73.0%, respectively). Meanwhile, CD8_Mem_T cluster had more clonal expanded cells, only 40.6% were singletons, consisting with their cytotoxic and memory features (Fig. 5D). Cells in the post-treatment samples were more clonal expanded than the pre-treatment samples (Fig. 5E). In addition, we observed multiple clonotypes were already detected in pre-treatment samples then further expanded in the post-treatment samples (Fig. 5F), suggesting they were potential effector cells, expanding during the disease progression to perform anti-tumor functions.

Fig. 5figure 5

T cell subtypes and dynamics in tumor microenvironment. (A and B) T cell subtypes are identified by canonical markers and unique cell type features. (C) T cell subtype composition changes during the disease progression. CD8 memory T cells are increased post-treatment. (D) Clonal size distribution among the T cells. Clonal expanded cells are enriched in CD8 memory T cells. (E) Barplot to show the clonal size distribution in each sequenced sample, post-treatment samples have less singletons, suggesting lower diversity (F) Pie plot to show the clonotypes with clonal size >  = 5. Pie size represents the size of the clonotype, color indicate the cell origins. (G) Selected differentially expressed genes between post- and pre-treatment cells in each T cell subtype

We next analyzed differentially expressed genes in T-cell subsets pre- and post-treatment (Fig. 5G, Tables S4-S6). We observed noticeable changes in the transcriptomes of these cells, involving genes that related to tumor invasion, metastasis, and therapeutic resistance. For example, CD69 was increased post-treatment in all three T cell subsets. CD69 was demonstrated to regulate T-cell migration and retention in tissues, playing an important role in inducing the exhaustion of tumor-infiltrating T cells [30]. ANXA1 was increased in CD8 memory T cell post-treatment, which has been reported to inhibit the anti-tumor immunity and support the formation of an immunosuppressed tumor microenvironment that promotes tumor growth and metastasis [31].

We observed substantial increase of monocytes in the PBMC samples (Fig. 2C), suggesting monocytes may also play a role in the disease progression. Here, we further evaluated the heterogeneity of the monocyte population. We identified two major monocyte subtypes, the classic monocytes (Classic_Mono) and the CD16 + Monocytes (CD16_Mono) (Fig. 6A and B). The CD16 + nonclassical monocytes might have a function to scavenge tumor cells and debris due to their patrolling activities [32]. It accounted 25.4% and 58.5% of the monocyte population in BM1 and PBMC1, respectively, but decreased to 9.9% and 26.3% post-treatment. These findings demonstrated the actively reforming TME that may influence tumor development and therapeutic outcomes.

Fig. 6figure 6

Monocyte subtypes and composition changes in tumor microenvironment. (A) UMAP presentation of monocyte subtypes. (B) Heatmap to present the canonical marker genes used to identify the classical monocytes and the CD16 + nonclassical monocytes. (C) Barplot to show the monocyte subtype composition changes pre- and post-treatment. CD16 + nonclassical monocytes decreased post-treatment

Rewired Tumor-microenvironment interaction networks during treatment

As we observed substantial changes in TME during the disease progression, we further performed cell–cell interaction analysis to explore the communication dynamics between tumor cells and other TME cells. We revealed significantly decreased cell–cell interactions in post-treatment cells, for both total interaction counts and interaction strengths (Fig. 7A and B). The total interaction counts were 1487 in BM1 cells, decreased to 913 in BM2 cells. Among them, tumor cell involving interactions decreased from 709 to 418, contributing to 50.7% of the decrease (Fig. S9). We also observed that several signaling pathways were affected by the interaction network rewiring (Fig. 7C). For example, cell–cell interactions in MHC-1, MIF, GALECTIN and ICAM signaling pathways became weaker in post-treatment cells. We further evaluated each ligand-receptor pair between tumor and other major immune cells (Fig. 7D). As examples, MIF-(CD74 + CXCR4) interaction was decreased between tumor and most immune cell types post-treatment. MIF/CXCR4 axis were reported to induce an aggressive phenotype by inducing proliferation, adhesion, migration, and invasion of the colon cancer cells [33]. ADGRE5 (CD97)-CD55 interaction were also decreased post-treatment, an interaction has been shown to contributing to aggressiveness in multiple cancers [34,35,36]. Together, our results suggested that the interaction between malignant γδ T cells and the TME were actively reshaping during the disease progression. The Tumor_2 cells had cut down the communications with the TME, however might acquire survival advantages.

Fig. 7figure 7

Rewired cell–cell interaction networks during disease progression. (A) Overview of the cell–cell interactions inferred from ligand-receptor expressions. Post-treatment samples have less interactions than the pre-treatment samples. (B) Scatter plot of the incoming and outgoing signals for each cell type. Tumor cell signals are reduced post-treatment. (C) Heatmap to compare the signaling patterns pre- and post-treatment. (D) Bubble plot to show the increasing and decreasing ligand-receptor interactions between tumor cells and major immune cells. Bubble sizes indicate p values and colors indicate average expression of ligand and receptor in the given interaction

留言 (0)

沒有登入
gif