Identification of evolutionary mechanisms of myelomatous effusion by single-cell RNA sequencing

Multiple myeloma (MM) is characterized by the clonal proliferation of malignant plasma cells (PCs) and accounts for 10% of malignant hematopoietic diseases.¹ The development of protease inhibitors, immunomodulatory drugs, and monoclonal antibodies has increased survival rates.² However, MM is incurable, and disease progression results in drug resistance and poor outcomes.³ When a subclonal of PCs become independent of marrow microenvironment and egress out of bone marrow (BM), it is termed as extramedullary multiple myeloma (EMM). The overall incidence of EMM is 6% at the initial diagnosis and 6% to 20% during the follow-up.⁴ EMM can spread to skeletal tissue and extraosseous sites, including lymph nodes, liver, skin, central nervous system, pleura, and blood.⁵ Extraosseous MM is associated with worse outcomes.^6,7 

Myelomatous effusion (ME) is manifested with the presence of neoplastic PCs in bodily fluids, including pleural effusion and ascites. ME occurs in <1% of MM cases.^8,9 Although a few cases of remission with long-term survival have been reported,^10,11 the median survival of patients with ME is only 4 months.^12,13 The pathogenesis of ME has not been fully elucidated.¹⁴ The increased proliferation of PCs and the activation of antigen presentation, proteasomes, oxidative phosphorylation, glycolysis, and immune evasion in these cells are involved in the formation of EMM.¹⁵ Pleural and peritoneal lesions and hematogenous metastasis may be related to ME.¹⁶ However, the molecular mechanisms underlying ME need to be studied further.

This study included single-cell RNA sequencing (scRNA-seq) of malignant PCs as well as immune cells in 3 patients with MM who relapsed with ME. BM, peripheral blood (PB), and pleural effusion/ascites were collected. The transcriptome analysis of PCs in these tissues allowed for the identification of cell subclones involved in the migration to extramedullary sites. We determined the molecular mechanism of ME formation and the impact of target genes on myeloma invasiveness using in vitro experiments. Our findings underscore the importance of analyzing single-cell transcriptome profiles to identify PC clones able to migrate and survive in extramedullary sites.

Primary BM, PB, and pleural effusion/ascites samples were collected from 3 patients with MM with ME at the Department of Hematology of Jiangsu Province Hospital from January 2019 to December 2021. The main properties and quality statistics of the experiment are summarized in supplemental Materials. All procedures in this study were approved by the ethics committees of Jiangsu Province Hospital and in accordance with the Declaration of Helsinki. All the enrolled patients signed an informed consent.

Cells (1 × 10⁵ cells per mL) were harvested and suspended in equal volumes of phosphate-buffered saline. Cells were loaded onto a microfluidic chip, and scRNA-seq libraries were constructed using a GEXSCOPE Single-Cell RNA Library Kit¹⁷ (Singleron Biotechnologies, Nanjing, China). Individual libraries were diluted to 4 nM and pooled, and the pools were sequenced using an Illumina NovaSeq 6000 platform to produce 150 bp paired-end reads.

Raw reads from scRNA-seq were processed to generate gene expression matrixes using CeleScope version 1.1.7 (https://github.com/singleron-RD/CeleScope). Low-quality reads were removed using CeleScope, and poly-A tails and adapter sequences were removed using Cutadapt2 version 1.17. Cell barcodes and unique molecular identifiers (UMIs) were extracted. High-quality reads were mapped to the reference genome GRCH38 (Ensembl version 92) using STAR version 2.6.1a.¹⁸ UMI counts and gene counts were grouped using featureCounts¹⁹ version 2.0.1 to generate expression matrixes.

Cells with a range from 200 to 5000 genes, less than 30 000 UMI counts, and less than 50% mitochondrial content were retained; 37 731 cells passed quality control. Six samples from GEXSCOPE (MM01-ME, MM01-PB, MM01-BM, MM03-BM, MM03-ME, and MM03-PB) and 2 samples from 10X Genomics (MM02-BM and MM02-ME) were integrated using Harmony version 1.0. Dimensionality reduction, unsupervised clustering, and data were visualized using Seurat²⁰ version 3.1.2. Gene expression was normalized using the Normalize Data function of the Seurat package, in which the number of UMIs of each gene was divided by the total UMIs of each cell, multiplied by 10 000, and log-transformed (ln [UMI-per-10 000 + 1]). Highly variable genes were identified using the Find Variable Genes function and analyzed using principal component analysis. Clustering with 20 principal components at a resolution of −0.8 was performed via graph-based clustering. Cells were visualized in a 2-dimensional space using a uniform manifold approximation and projection (UMAP) algorithm and the RunUMAP function of the Seurat package.

The cell types in each cluster were identified using canonical markers from literature data and the reference SynEcoSys database (Singleron Biotechnologies, Nanjing, China). Cell doublets were estimated based on the expression pattern of canonical cell markers. Clusters enriched with multiple cell type–specific markers were excluded from the analysis. Subtype analysis of major cell types is summarized in supplemental Materials.

Differentially expressed genes (DEGs) were determined using the FindAllMarkers function of the Seurat package (Wilcoxon rank-sum test with Bonferroni correction for multiple comparisons). For computing DEGs, genes expressed in at least 10% of PCs in either of the 2 samples being compared were selected for analysis. Cut-off criteria included the absolute value of difference in gene expression on a natural log scale of at least 0.25 and adjusted P value < .05. Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes enrichment analyses were performed using clusterProfiler6²¹ version 3.16.1.

Survival analysis was performed using the survminer package in R. Gene expression and survival data from the The Cancer Genome Atlas and MMRF CoMMpass data sets were downloaded from the Genomic Data Commons platform. The MMRF CoMMpass data set contains 859 RNA sequences from CD138⁺ PCs in BM samples from patients with MM. Data on overall survival or progression-free survival were available for 787 patients. The samples were divided into 2 groups (high and low gene expression) based on median values. Survival was assessed using the Kaplan-Meier method and compared using the 2-sided log-rank test.

Single-cell regulatory network analysis was performed using pySCENIC²² version 0.11.0 in Python, following the protocol described in the single-cell regulatory network inference and clustering (SCENIC) workflow. Coexpression gene regulatory networks were generated using the “pyscenic grn” command and the GRNBoost2 method. AUCell analysis was performed using the “pyscenic aucell” function with the following parameters: rank_threshold = 5000, auc_threshold = 0.05, and nes_threshold = 3.

The generation of LILRB4 knockout (KO) NCI-H929 cells and LILRB4-overexpressing (OE) NCI-H929 cells and procedures of cell culture, cell proliferation and viability analysis, western blot, real-time quantitative polymerase chain reaction, cell cycle analysis, transendothelial migration assay, tube formation in human umbilical vein endothelial cells (HUVECs), and flow cytometry can be found in the supplemental Materials.

We evaluated 3 patients with MM who relapsed with ME. The disease course is summarized in Figure 1A, and patient characteristics are shown in Table 1. Malignant PCs in MEs were detected using flow cytometry. BM, PB, and ME were collected from each patient and analyzed using scRNA-seq. After quality control, the single-cell transcriptomes of 37 731 cells were selected (supplemental Figure 1A). In addition, 17 599 malignant and 20 132 nonmalignant cells were distinguished using 3 complementary approaches. After the integration, single-cell transcriptomes were visualized in a two-dimensional map using the UMAP method for unsupervised dimensionality reduction (Figure 1B). The BM of all patients contained PCs, immune cells, myeloid cells, and erythroblasts (Figure 1C), whereas the pleural and peritoneal effusion contained abundant PCs and few immune components. The heatmap of the top 10 DEGs for each cell type is shown in supplemental Figure 1B.

Then, we investigated the immune microenvironment within or outside the BM. The predominant immune cell types were T cells (Figure 1D), indicating humoral immunodeficiency in patients with myeloma. The predominant T-cell cluster (with ∼3546 cells) contained the following 4 cell subsets: naive T cells, CD8⁺ effector T cells, γδ T cells, and regulatory T cells. MEs contained a higher percentage of cytotoxic T cells, whereas PB contained a higher proportion of naive T cells (Figure 1D; supplemental Figure 2A). Of note, there was a population of Treg cells in the effusion group, indicating that reduced immune response might contribute to the formation of ME (Figure 1E).

We analyzed different components of T cells in PB and ME. We found that T cells in the ME highly expressed proliferation protein (eg, MKI67), proteins involved in the regulation of the microtubule filament system (eg, STMN1), and mobility protein (eg, HMGB2), indicating that T cells in the ME were proliferating (supplemental Figure 2B). CD8⁺ effector T cells differed in the expression of cell cycle proteins (eg, TOP2A) and were associated with cell proliferation (Figure 1F; supplemental Figure 2C). We observed higher cytotoxic scores of the T cells in the BM than those in extramedullary lesions, whereas there were no statistically significant differences in exhaustion scores between that within and outside the BM (supplemental Figure 2D). We observed the same phenomenon in CD8⁺ effector T cells (supplemental Figure 2E).

Next, we determined the scRNA-seq profile of the BM from 3 patients with MM for matching extramedullary involvement from these patients (ME and PB). scRNA-seq was not performed in 1 PB sample because the number of cells was small. A total of 17 599 malignant cells clustered based on the tumor of origin (Figure 2A). We removed immunoglobulin genes from the clustering analysis as the high expression of clones limited the assessment of tumor heterogeneity. Different clusters of malignant PCs were identified in the BM, ME, and PB (Figure 2B). Some clusters mainly derived from single samples, whereas others included samples from both the extramedullary site and BM. Every cluster contained unique signature genes, suggesting the intertumoral heterogeneity of PCs. DEGs in each cluster are shown as a heatmap (Figure 2C). There was considerable variability between the patients and their sites (BM vs ME). Clusters 1, 4, and 6 contained PCs originated from ME samples, expressing cellular stress response protein (SGK1; cluster 1), MT1X (cluster 4), and GNB2L1 (cluster 6) gene that played a role in myeloma cell survival, migration, and contact with the extracellular matrix. The myeloma cells in BM samples were clusters 2, 5, and 7 expressing the regulator of G-protein signaling family RGS1 (cluster 2), cyclin family gene CCND2 (cluster 5), and the inflammatory responses gene ANXA1 (cluster 7) that regulated the cell cycle and promoted the growth of myeloma cells in the BM. Although PC clusters were site- and patient-specific, myeloma cells in the BM and extramedullary sites shared the same subclusters, indicating the presence of a transition state of extramedullary migration. MM01 presented PC clusters 3, 8, 10, and 11, whereas MM03 presented PC clusters 9 and 11, which existed both inside and outside the BM; these clusters contained genes encoding cell cycle regulation proteins (eg, CCND2), proliferation marker proteins (eg, MKI67), chemokine receptors (eg, CXCR4), extracellular matrix proteins (eg, MMP9), and high-mobility group proteins (eg, HMGB2), respectively. Clusters 3, 8, and 10 contained genes encoding major histocompatibility complex (MHC) class II molecules, such as HLA-DRA and HLA-DRB1 (supplemental Table 1).

Next, we performed differential gene expression analysis of PCs from the BM and extramedullary site. The expression of the metallothionein family (eg, MT1X and MT2A) and SGK1 increased in extramedullary PCs, whereas gene encoding S100 family of proteins (eg, S100A8 and S100A9) and hemoglobin (eg, HBA1, HBA2, and HBB) were upregulated in PCs of the BM (Figure 2D).

Signaling pathways enriched in malignant PCs from the BM and extramedullary sites were analyzed. Extramedullary disease upregulated pathways associated with cellular response to oxidative stress, protein kinase regulator activity, and reactive oxygen species metabolic process, and downregulated pathways associated with protein folding and endoplasmic reticulum (ER) response (supplemental Table 2).

The malignant PCs of MM01 were heterogeneous, which made up of clusters 3, 5, 7, 8, 10, and 11. Kyoto Encyclopedia of Genes and Genomes analysis showed that the actin cytoskeleton pathway and leukocyte transendothelial migration were upregulated in cluster 11 (supplemental Figure 3A). GO enrichment analysis indicated that the focal adhesion and cell-substrate conjunction were downregulated in the cellular component of cluster 10 (supplemental Figure 3B). These 2 clusters (10 and 11) were mainly present in the BM and PB but rarely in the ME. Cluster 8 not only existed in the BM but also occupied the majority portion in the PB and ME. Proteasome and antigen presentation signal pathway were upregulated in this cluster (supplemental Figure 3C).

MM02 presented PC clusters 3, 4, 6, and 7. Clusters 3 (13 cells in the ME of MM02) and 7 were observed in both the BM and ME. GO analysis indicated that cell cycle control and oxidative phosphorylation pathways were enriched in cluster 3, and ER and unfolded protein response (UPR) pathways were enriched in cluster 7 (supplemental Figure 4A-B). Clusters 4 and 6 were found in the ME and were associated with the upregulation of reactive oxygen species pathways (supplemental Figure 4C-D).

MM03 showed PC clusters 1, 2, 8, 9, and 11. Clusters 1 (8 cells in the BM of MM03) and 2 (2 cells in the ME of MM03) were found in the BM and ME and were associated with the upregulation of ER pathways (supplemental Figure 5A-B). In addition, cell cycle control and oxidative phosphorylation pathways were upregulated in cluster 9 (supplemental Figure 5C).

We performed hotspot analysis and found that the gene modules 3 and 7 were highly expressed in intramedullary and extramedullary PC clusters, suggesting that these modules were associated with cell viability and migration (Figure 3A). The gene sets in these modules were present in PCs both from extramedullary lesions and BM (Figure 3B). We observed a higher expression of modules 3 and 7 in the PCs of the BM compared to PB and ME, according to UCell Gene Set Scoring²³ (supplemental Figure 6A). We speculated that the gene set was a key factor in the survival and proliferation of myeloma cells migrating into the pleural space, and we defined cell clusters expressing this gene set as EMICs. Hub genes in module 3 (TOP2A, NUSAP1, and TPX2) correlated with nuclear division, whereas hub genes in module 7 (MKI67 and HMGB1) were implicated in cell proliferation (supplemental Table 3). We further subdivided the functions of the genes and found that the functions of most of them were related to the cell cycle, P53 signal pathway, and other cellular functions (Figure 3C), suggesting that the improvement in these functions could help tumor cells adapt to the extramedullary environment. We subsequently explored the expression of modules 3 and 7 in GSE66291 and GSE39683 data sets. Single sample gene set enrichment analysis was performed on 8 newly diagnosed multiple myeloma (NDMM) cases and 8 secondary plasma cell leukemia (sPCL) cases from GSE66291 and 7 NDMM cases and 7 sPCL cases from GSE39683. We observed higher enrichment scores of modules 3 and 7 in the sPCL cases, indicating that the EMICs played an important role in the extramedullary migration of PCs (supplemental Figure 6B).

There were 68 genes in modules 3 and 7. Overlapped with 50 genes of patients with relapsed/refractory multiple myeloma from the DisGeNET database, we suggested that ANP32B and PCNA might be the biomarkers of EMICs (supplemental Figure 6C).

We analyzed the relationship between gene expression and prognosis in 2 MM cohorts from the MMRF CoMMpass data set (n = 787). The high expression of genes in these modules was associated with worse overall survival in patients with MM (Figure 3D), whereas the other modules had no significant effect on prognosis (supplemental Figure 7A).

Then, we considered how EMICs relate to genetic heterogeneity and tumor extramedullary evolution. We assessed genetic differences between EMICs and other tumor cells. EMICs expressed microtubule-associated genes (TUBA1B and NUSAP1), proliferative genes (TTK and BIRC5), and cell cycle–associated genes (CCNB1 and FAM111A) (supplemental Table 4). In addition, we identified the molecular pathways in which the differential genes were expressed. The cell cycle pathway was upregulated, whereas protein processing in ER pathway was downregulated in EMICs (supplemental Figure 7B). Genetic changes in extramedullary clones and intramedullary clones indicated they are substantially heterogeneous.

Furthermore, considering that both intramedullary and extramedullary cells possess EMICs, we further examined the genetic differences between them. The expression of leukocyte immunoglobulin–like (LIL) family proteins LILRB4 and long noncoding RNAs (MALAT1 and NEAT1) were upregulated in extramedullary EMICs (supplemental Table 5). Intramedullary EMICs had the potential to migrate extramedullary after acquiring some extra genetic alterations.

Finally, we investigated the different transcription-factor expression levels between EMICs and other tumor cells. Cells that highly expressed gene signatures of modules 3 and 7 exhibited specific expression of transcription factors essential for the survival and growth of MM tumors such as E2F1, YY1, and SMAD1 (supplemental Figure 7C). These 3 genes were associated with transforming growth factor β and NF-κB signaling pathways, suggesting that the NF-κB pathway played a crucial role in myeloma migration to extramedullary sites. Given that transcription factors were important in defining cell identity, we then compared the expression levels of transcription factors of intramedullary and extramedullary EMICs and found high activities of RUNX2 and antigen presentation (MHC class II) pathway in extramedullary EMICs. Upregulated transcription factors, such as PRDM1, STAT1, and IRF family, related to the master regulator of MHC class II gene expression–MHC Class II transactivator, might reduce T-cell–mediated cytotoxicity (supplemental Figure 8A). These results suggested that EMICs were involved in tumor cell proliferation and migration and immune evasion.

Gene signatures associated with myeloma cell migration were analyzed. The genes LILRB4, H3F3B, ARPC2, TMEM173, LAP3, and CD63 were upregulated, whereas FOS, SLC1A5, BCMA, DUSP1, P4HB, and JUN were downregulated at extramedullary sites (Figure 4A).

LILRB4, an immune inhibitory receptor that contains putative immunoreceptor tyrosine–based inhibitory motifs, was more highly expressed in the ME than in the BM (Wilcoxon, P < 2.2e-16; Figure 4B). The analysis of the MMRF CoMMpass data set showed that the high expression of LILRB4 or ARPC2 was associated with poor survival, whereas H3F3B, TMEM173, and LAP3 had no significant effect on prognosis (Figure 4C; supplemental Figure 8B). Furthermore, the LILRB4 gene was significantly upregulated in extramedullary EMICs (supplemental Figure 8C).

To verify the role of LILRB4 on ME development, we used the CRISPR-Cas9 strategy to knock out LILRB4 (LILRB4-KO), and overexpress LILRB4 (LILRB4-OE) via lentiviral transfection in MM cell line (NCI-H929). Western blot (Figure 4D) and real-time polymerase chain reaction (Figure 4E) were performed to verify the overexpression or KO of LILRB4. Then, we detected the level of apoptotic marker cleavage of PARP and caspase 3 in the LILRB4-OE and LILRB4-KO cell lines. The cleaved PARP and caspase 3 were upregulated in LILRB4-KO cells and downregulated in LILRB4-OE cells (Figure 4F).

The effect of LILRB4 on cell cycle progression was also assessed. LILRB4-OE induced arrest in the G1 phase, whereas the percentage of cells in the G1 phase was strikingly decreased in LILRB4-KO cells (Figure 4G-H). There was no significant difference in cell proliferation between LILRB4-OE/KO and control H929 cells using the CCK-8 assay, indicating that LILRB4 did not induce cell proliferation (Figure 5A-B). The effect of LILRB4 on myeloma cell migration was evaluated using a transwell chamber (Figure 5C) and tube formation in HUVECs. LILRB4 overexpression increased cell migration (Figure 5D-E) and promoted the tubular structure of HUVECs (supplemental Figure 9A-B). Meanwhile, we observed decreased transendothelial migration (Figure 5F) and disrupted the tubular structure of HUVECs in LILRB4-KO (supplemental Figure 9C). Moreover, LILRB4 changed the messenger RNA expression of angiopoietin 1, integrin αv, and integrin β3 (Figure 5G-I).

EMM is a malignancy that presents as a soft tissue mass of monoclonal PCs. Pleural and peritoneal involvement, including ME, is rare, accounting for less than 1% of MM cases.²⁴ ME was assessed in case reports^25,26 and a systematic review.¹² Nonetheless, the molecular mechanisms underlying this entity are incompletely understood. The use of scRNA-seq analysis allows for determining the transcriptional profile of immune cells in the BM and extramedullary sites.

RNA-seq analysis of the immune microenvironment showed that the genes implicated in cell proliferation, migration, and cell cycle progression were highly expressed in effector T cells from MEs. The expansion of effector T-cell subtypes does not produce a cytotoxic effect on tumor.²⁷ Genes associated with T-cell exhaustion, including PD-1, CTLA-4, 2B4, and CD160, were expressed in T cells from MM,²⁸ especially in extramedullary sites,¹⁵ indicating that the dysregulation of T-cell immune surveillance might play a role in ME formation.

EMM has a more aggressive clinical course than myeloma.²⁹ ScRNA-seq helps to assess the evolution of malignant PCs.³⁰ We found that clusters 3, 8, 9, 10, and 11, which were both present at BM and extramedullary sites, were associated with PC migration and shared proliferation marker genes (eg, MKI67) and cell migration–related genes (CXCR^31,32 and HMGB^33,34 family), which were related to poor prognosis in MM. Consistent with a previous study,¹⁵ these clusters were associated with the proteasome, cell cycle control, antigen presentation, and oxidative phosphorylation pathways.

Cancer stem–like cells (CSCs) contribute to cancer progression and metastasis.^35-37 However, CSCs in MM have not been thoroughly defined.³⁸ ME often occurs during relapse, and disease progression and metastasis are driven by clonal selection.³⁹ We identified a PC subtype in the BM, PB, and ME through hotspot analysis, and this subtype survived in extramedullary sites. In addition, this subtype expressed genes related to cell cycle control, chromatin remodeling, microtubule dynamics, and other cellular functions. Cell cycle dysregulation⁴⁰ and higher proliferation indices⁴¹ contribute to drug resistance and poor outcomes. In addition, we determined that the expression of modules 3 and 7 was higher in the BM of sPCL than that of NDMM using 2 other databases (GSE39683 and GSE66291). On the basis of data from a recent study, there was significant tumor heterogeneity between BM and extramedullary sites in MM.³ We speculated that EMICs, similar to CSCs, were a high-risk PC clone and more difficult to remove.⁴² It was necessary to identify biological markers for predicting the potential of developing ME. Two potential genes were filtered out by overlapping the signature of relapsed/refractory multiple myeloma and EMICs. ANP32B, a potent survival factor for exogenous cytokine–dependent myeloma cells, upregulated the expression of the 2 major antiapoptotic proteins in myeloma cells, Mcl-1 and Bcl-2.⁴³ PCNA was an essential factor for DNA replication and repair. Its expression on MM cells, interacting with the NKp44 receptor, induced immuno-compromise by inhibiting NK cell function.⁴⁴ 

The transcription factors E2F1, SMAD1, and YY1 were upregulated in EMICs. E2F1, a gene encoding DNA-binding protein, was indirectly activated by cyclin D1 and induced the growth of MM.^45,SMAD1 was a determinant of drug resistance in MM cells and an unfavorable prognostic marker,⁴⁶ whereas a higher nuclear YY1 expression was associated with disease progression.⁴⁷ 

EMICs promoted tumor cell survival in extramedullary sites, indicating that EMICs were the prerequisite for extramedullary metastasis. Thus, we speculated that based on the presence of EMICs, further genetic alterations contributed to extramedullary infiltration. A set of genes was consistently expressed in PCs in the BM and extramedullary sites. LILRB4, an inhibitory receptor of immune checkpoint pathways,⁴⁸ was upregulated in ME and was known to suppress the immune response in solid tumors⁴⁹ and promote the infiltration of myeloid leukemia cells.⁵⁰ Furthermore, LILRB4 induced nonsmall lung cancer invasion by increasing epithelial-mesenchymal transition and the expression of vascular endothelial growth factor A.⁵¹ The KO of LILRB4 in MM cells decreased cell migration, and LILRB4 overexpression had the opposite effect. These results suggested that LILRB4 was involved in myeloma extramedullary migration. Besides this, we also detected the association between high LILRB4 expression and poor outcomes. Of interest, we found that BCMA, an immunotherapy target for MM,⁵² was downregulated in extramedullary sites. The efficacy of anti-BCMA chimeric antigen receptor T-cell therapy for EMM is limited.^53,54 

In summary, the transcriptome analysis of PCs and the immune microenvironment at different physiological sites revealed that patients with EMM had a unique gene expression profile, including gene sets associated with cell migration and survival.

This study has limitations. Firstly, the number of ME samples was small because the incidence of this entity was low. Notwithstanding, this study helped elucidate the mechanisms underlying MM cell proliferation and invasion. Through the analysis of intratumor heterogeneity, we believe that the development of ME needed a subpopulation of proliferating stem–like cells, whereas invasion-related genes such as LILRB4 mediated the migration of PCs.

The authors thank Xiaming Kong for helping with data analysis of single-cell sequencing.

This study was supported by National Natural Science Foundation of China (81720108002), CSCO Research Foundation (Y-Roche2019/2-0090), Science Foundation Project of Ili & Jiangsu Joint Institute of Health (yl2021ms04), and Jiangsu Province Hospital (the First Affiliated Hospital with Nanjing Medical University) Clinical Capacity Enhancement Project (JSPH-MB-2021-14 and JSPH-MB-2021-9).

Contribution: X.Q., J.L., Z.S., and J.J. contributed to the conception and design; Y.L., Y.C., and F.L. contributed to the acquisition of data; X.Q., Z.S., and J.J. contributed to the analysis and interpretation of data; X.Q., Z.S., and J.J. contributed to the writing, review, and/or revision of the manuscript; X.Q. and J.L. supervised the study; and all authors read and approved the final manuscript.

Conflict-of-interest disclosure: The authors declare no competing financial interests.

Correspondence: Jianyong Li, Department of Hematology, The First Affiliated Hospital of Nanjing Medical University, Jiangsu Province Hospital, 300 Guangzhou Rd, Nanjing, 210029, China; e-mail: lijianyonglm@126.com; and Xiaoyan Qu, Department of Hematology, The First Affiliated Hospital of Nanjing Medical University, Jiangsu Province Hospital, 300 Guangzhou Rd, Nanjing, 210029, China; e-mail: quxiaoyan@jsph.org.cn.