Yuanqi Wang,Xiao Xiang,Huadong Chen,Luyao Zhou,Shuling Chen,Guopei Zhang,Xiaofei Liu,Xuxin Ren,Juncheng Liu,Ming Kuang,Juan Jiang,Jinbiao She,Zhichong Zhang,Ruidong Xue,Hong Jiang,Ji Wang,Sui Peng
Intratumoral erythroblastic islands restrain anti-tumor immunity in hepatoblastoma
hepatoblastoma,erythroid cells,erythroblastic island,immunosuppression
Summary Erythroblastic islands (EBIs) are the specialized structures for erythropoiesis, but they have never been found functional in tumors. As the most common pediatric liver malignancy, hepatoblastoma (HB) requires more effective and safer therapies to prevent progression and the lifelong impact of complications on young children. However, developing such therapies is impeded by a lack of comprehensive understanding of the tumor microenvironment. By single-cell RNA sequencing of 13 treatment-naive HB patients, we discover an immune landscape characterized by aberrant accumulation of EBIs, formed by VCAM1+ macrophages and erythroid cells, which is inversely correlated with survival of HB. Erythroid cells inhibit the function of dendritic cells (DCs) via the LGALS9/TIM3 axis, leading to impaired anti-tumor T cell immune responses. Encouragingly, TIM3 blockades relieve the inhibitory effect of erythroid cells on DCs. Our study provides an immune evasion mechanism mediated by intratumoral EBIs and proposes TIM3 as a promising therapeutic target for HB.
Intratumoral erythroblastic islands restrain anti-tumor immunity in hepatoblastoma

Erythroblastic islands (EBIs) are the specialized structures for erythropoiesis, but they have never been found functional in tumors. As the most common pediatric liver malignancy, hepatoblastoma (HB) requires more effective and safer therapies to prevent progression and the lifelong impact of complications on young children. However, developing such therapies is impeded by a lack of comprehensive understanding of the tumor microenvironment. By single-cell RNA sequencing of 13 treatment-naive HB patients, we discover an immune landscape characterized by aberrant accumulation of EBIs, formed by VCAM1+ macrophages and erythroid cells, which is inversely correlated with survival of HB. Erythroid cells inhibit the function of dendritic cells (DCs) via the LGALS9/TIM3 axis, leading to impaired anti-tumor T cell immune responses. Encouragingly, TIM3 blockades relieve the inhibitory effect of erythroid cells on DCs. Our study provides an immune evasion mechanism mediated by intratumoral EBIs and proposes TIM3 as a promising therapeutic target for HB.

Pediatric cancer is the leading cause of mortality in children., In recent decades, limited treatments for them have been discovered, especially for pediatric solid tumors. A combination of non-specific chemotherapy and surgery is still the standard of care. However, mortality and recurrence rates are high in patients with refractory and metastatic pediatric cancers. Notably, dose-dense chemotherapy brings a high incidence of severe adverse events, which could have permanent effects on the learning, development, and quality of life of pediatric patients., Therefore, novel therapies with improved effectiveness and safety are urgently needed for pediatric patients with solid tumors. Although immune checkpoint blockades (ICBs) have shown safety and clinical benefits in multiple adult cancers,, only a few clinical trials of ICBs (anti-PD-1, PD-L1, or CTLA-4) in pediatric cancers have been conducted. All of these trials had a low incidence of adverse events as always, but unfortunately had limited response rates.,, It is unclear whether the low effectiveness of ICBs for pediatric solid tumors is ascribable to a lack of tumor antigens or the immune-suppressive tumor microenvironment (TME) or both. Due to the embryonic origin and low environmental exposure to carcinogens, pediatric cancers exhibit a lower mutation burden than adult cancers., Moreover, pediatric cancers are characterized by less infiltration of T cells, which is regarded as one of the reasons underlying the low response rate to ICBs. A better understanding of the immune status of pediatric tumors can aid in the discovery of new molecular targets for immunotherapies other than PD-1/PD-L1 and CTLA-4. Hepatoblastoma (HB) is the most common pediatric liver cancer, particularly among children under the age of 3 years. However, studies of its TME have lagged behind those of other pediatric cancers. HB has been proposed to originate from the aberrant differentiation of hepatocyte progenitors in the liver, which undergoes several major changes during the early postnatal period., For example, the fetal liver has abundant hematopoietic progenitor cells, which diminish within 5 weeks of birth. Recent studies have described a brief immune landscape of HB, using immunohistochemistry (IHC) and bulk RNA sequencing,, and tumor cell heterogeneity by performing single-cell RNA sequencing (scRNA-seq)., But a comprehensive picture of TME at single-cell resolution is still missing. Furthermore, patients in these studies mostly received neoadjuvant chemotherapy, which might alter the TME and even induce tumor pathological type conversion., The immunological milieu of treatment-naive HB is crucial for the development of new immunotherapies that are safer and more effective than current chemotherapies for children. Here, we mapped the cellular landscape of TME in 13 treatment-naive HB patients at the single-cell resolution, revealing that erythroblastic islands (EBIs) composed of erythroid cells and VCAM1+ macrophages were significantly enriched in HB and inversely correlated with the long-term survival of patients. These intratumoral erythroid cells suppressed the function of dendritic cells (DCs) via the LGALS9/TIM3 axis, resulting in a decreased activated state of CD8+ T cells. Encouragingly, anti-TIM3 antibodies could significantly reverse the suppression of DCs, proposing TIM3 as a potential target of immunotherapies for HB.

To survey the TME in HB, we performed scRNA-seq on 23 samples from 13 treatment-naive HB patients as a discovery cohort, including 13 tumors, 7 paired non-tumor adjacent liver tissues, and 3 peripheral blood samples. Moreover, we collected paired tumor and adjacent liver samples from an additional 19 treatment-naive HB patients to perform multiplex immunofluorescence staining (mIF) and bulk RNA-sequencing validation. Detailed clinical and pathological information on all the enrolled patients is listed in Table S1. Other validation cohorts were created using the public sequencing data of HB and human fetal livers, to confirm our findings from the discovery cohort (Figure 1A). After discarding doublets and cells with low mRNA quality, we obtained the transcriptomes of 45,204 and 100,883 cells for the non-tumor adjacent liver tissues and tumor samples, respectively (Figure S1A; Table S2). To define the cellular composition in the TME, we performed a graph-based clustering analysis on these cells and identified eight major cell types (Figure 1B). Cell subclustering revealed 63 subcategories based on known signature genes or specific genes expressed by each cell type (Figure 1C). Transcriptional features of all subclusters are provided in Table S3. We first identified tumor cells based on single-cell copy-number variation (CNV) speculation. Copy-number profiles were homogeneous in different patients, with gains of 1q, 2q, and 15q and losses of 6q (Figure S1B). Next, we identified immunocytes, including myeloid cells, CD4+ T cells, CD8+ T cells, B cells, innate lymphoid cells (ILCs) (including natural killer [NK], ILC2, ILC3, and γδT cells), and progenitor cells. All these subsets were shared among tumor samples in similar proportions, in which myeloid cells accounted for the largest proportion, followed by ILCs (Figure 1D). Further exploring the distribution of these cells between tumor and adjacent liver samples, we discovered a lower infiltration of multiple immunocytes, especially T cells (Figure 1E). Antigen shortage may not be the key mechanism underlying the low infiltration of T cells, because HB tumor cells actually expressed high levels of tumor-associated antigens, including GPC3 and AFP,, at the single-cell level and bulk level (Figure S1C). Thus, suppressor cells may exist and play a pivotal role. A very unique feature of HB compared with adult tumors was the significant enrichment of hematopoietic progenitor cells in the TME, including hematopoietic stem cells/multipotent progenitors (HSCs/MPPs), megakaryocyte-erythroid-mast cell progenitors (MEMPs), granulocyte-myeloid cell progenitors (GMPs), erythroid cells (progenitors and precursors of red blood cells), and pre-progenitor B (pre-pro-B) and pro-B cells, all with high expression of signatures for progenitor cells resembling corresponding cells from fetal livers, and nearly no expression of CD45 (PTPRC) (Figures S1D and S1E). We then performed the pseudotime trajectory analysis and found that the HSC/MPP, MEMP, and GMP lineages were located at the beginning of differentiation (Figure S1F), confirming their progenitor characteristics. In accordance with the previous study, we also found that hepatic B cell lymphopoiesis remained in the liver after birth (Figure S1G). B cell receptor (BCR) sequencing data revealed that most B cells in HB were progenitor cells or immature B cells (Figure S1H), presenting a naive state of B lymphoid lineage. As a result, the TME of HB was fetal-like, characterized by less infiltration of adaptive immune cells, but enrichment of progenitor cells in the tumors.

Erythroid cells accounted for the major proportion of progenitors (Figure 2A), existing exclusively in tumors (Figure S2A). The only exceptions were two newborns (HB14, age 15 days, and HB07, age 21 days) who also had erythroid cells in their normal liver tissues. This was a normal physiological phenomenon, since erythroid cell hematopoiesis occurs in the liver during the fetal period and continues up to 5 weeks after birth. Nevertheless, the existence of erythroid cells in tumors of other patients, especially in those whose ages were far beyond 5 weeks, was abnormal. In addition, there was no significant correlation between erythroid cell proportion and patient age, further confirming that erythroid cell enrichment in HB was a tumor-related event rather than a delayed loss of hematopoiesis in the whole liver (Figure 2B). According to the expression pattern, we divided erythroid cells in HB into three stages (Figure S2B). Among them, early erythroid cells made up 8.6% (522/6,026) of erythroid cells in scRNA-seq and mainly contained erythroid progenitors due to the low expression of PTPRC and GYPA, including burst-forming unit-erythroid (BFU-E) and colony-forming unit-erythroid (CFU-E). Mid-to-late erythroid cells were mainly erythroid precursors, which were the main component of erythroid cells in HB. Mid-erythroid cells (57.4%, 3,456/6,026) mainly contained proerythroblasts (ProEs) and basophilic erythroblasts (BasoEs; PTPRC−GYPA+TFRC+ENG+). Late erythroid cells (34.0%, 2,048/6,026) mainly contained polychromatic erythroblasts (PolyEs) and orthochromatic erythroblasts (OrthoEs; PTPRC−GYPA+TFRC+ENG−). Thus, erythroid cells in HB were in a process of continuous dynamic change, from BFU-Es to OrthoEs. Flow cytometry analyses showed similar results in that the majority (∼85%) of erythroid cells in HB were erythroid precursors (CD71+CD235a+) (Figures S2C and S2D). Further Giemsa staining from fluorescence-activated cell sorting (FACS)-isolated erythroid cells revealed that these precursors are mainly ProEs, BasoEs, and PolyEs, while only a few cells entered the very late stage of erythropoiesis (OrthoEs < 10%) (Figure S2E). Flow cytometry analysis of paired HB tumor and adjacent liver samples confirmed a greatly higher number of CD71+CD235a+ erythroid cells inside tumors than in normal liver tissues (Figure 2C). Much higher levels of erythroid cell infiltration were also observed in HB tumors than in adjacent liver tissues in validation cohorts by immunohistochemical staining (Figure 2D). Bulk RNA sequencing on 15 pairs of HB tumor and adjacent liver tissues in our validation cohort also validated these results. GYPA and TFRC, two erythroid cell markers, were both found to be considerably more highly expressed in tumors (Figure 2E) as well as the marker gene set of erythroid cells (Figure S2F), indicating the higher infiltration of erythroid cells in HB tumors. Next, we compared intratumoral erythroid cells with conventional erythroid cells, which were found in the adjacent liver tissue of two newborn patients. To our surprise, erythroid cells in tumors were even more active in proliferation and differentiation, demonstrated by the upregulation of cell differentiation and metabolism-related pathways compared with conventional erythroid cells, which was further confirmed by flow cytometry analysis, cell-cycle analysis, and stemness analysis (Figures 2F and S2G–S2K), indicating the high proliferation of intratumoral erythroid cells. Furthermore, we calculated erythroid score in bulk RNA-sequencing datasets from our validation cohort, which was correlated with erythroid density in IHC staining and erythroid proportion in scRNA-seq (Figure S3A). Further long-term survival analysis revealed that HB patients with higher erythroid score had much shorter recurrence-free survival time in both Pavel’s cohort and Juan’s cohort (Figures 2G, S3B, and S3C). Similar results were observed in the progenitors of erythroid cells, including MEMPs and HSCs/MPPs, rather than other immunocytes (Figure S3D). Taken together, we identified intratumoral erythroid cells that were actively proliferating in HB tumors and correlated with long-term prognosis of HB patients.

Next, we investigated the origin of intratumoral erythroid cells. Erythroid cells have been found in adult hepatocellular carcinoma (HCC) tumors and were considered to be first induced in spleen due to tumor-induced anemia and then recruited into the tumor through the peripheral blood.,, To clarify whether erythroid cells in HB were also recruited from peripheral blood, we explored the cellular composition of peripheral blood mononuclear cells (PBMCs) from three HB patients. The results showed that PBMCs contained only a small amount of erythroid cells (Figure 2H), which rarely expressed the transcription factors (TFs) for erythroid cell differentiation (GATA1 and KLF1) (Figure 2I). Moreover, there was no correlation between erythroid cell proportion in tumors and patients’ serum hemoglobin levels (Figure 2J), indicating that erythroid cells in HB tumors were neither induced by anemia nor recruited from peripheral blood. Thus, we speculated that erythroid cells in HB tumors had different origins from adult HCC. Several studies suggested that HB tumor cells may derive from fetal hepatic precursor cells and reproduce the process of fetal liver development., We found that erythroid cells from HB and fetal livers both highly expressed fetal liver-specific erythroid TFs (MYC, HDAC2, HMGA1) and rarely expressed bone marrow-specific erythroid genes (ZNF385D, CXCL8, KLF2, TNFRSF12A, PRTN3) (Figure 2K). To investigate if erythroid cells in HB were more similar to erythroid cells from fetal liver, we calculated the gene signature scores of fetal liver-derived erythroid cells and bone marrow-derived erythroid cells. The results showed that erythroid cells in HB had a higher score for fetal liver-derived erythroid signature (Figure S3E). To further investigate the similarities between erythroid cells from HB and fetal livers, we used Harmony to integrate erythroid cells and related progenitors from our data with the public scRNA-seq dataset from fetal liver and bone marrow. The intermixing of these cells in a single-cell atlas indicated transcriptomic similarities between erythroid cells from HB tumors and fetal liver-derived erythroid cells, rather than bone marrow (Figure 2L). Collectively, the results suggested that erythroid cells in HB were most likely to originate in the fetal liver along with HB tumor development, rather than being recruited from peripheral organs or bone marrow. Two key questions remained to be explored. How could these intratumoral erythroid cells persist for such a long time and not diminish along with their counterparts in normal liver tissues? How do these intratumoral erythroid cells restrain anti-tumor T cell responses?

To find the answers to these questions, we first performed intercellular cross-talk analysis, revealing that erythroid cells had the closest interactions with macrophages and DCs (Figure 3A). To confirm the specific cell clusters that interacted with erythroid cells, we subclustered myeloid cells based on their marker genes (Figure 3B). Among five macrophage subclusters, VCAM1+ macrophage (cluster 1) had the highest score for phagocytosis, while CXCL10+ macrophages (cluster 5) and TREM2+ macrophages (cluster 4) had a higher score for M1 or M2 phenotype, respectively (Figure 3C). CXCL10+ macrophages showed their anti-tumor function by recruiting T cells, as previous studies reported (Figure S4A). TREM2+ macrophages highly expressed SPP1 (Figure S4B), which was associated with worse prognosis in multiple cancer types., Ligand-receptor (L-R) analysis suggested that the specific interactions between macrophage subclusters and erythroid cells were mainly related to cellular adhesion. Therein, VCAM1+ macrophages had the strongest cross talk with erythroid cells via VCAM1/α4β1 interaction (Figure 3D), which is the main L-R pair in EBIs. EBIs are the specific niche for erythroid hematopoiesis and are composed of a central macrophage surrounded by erythroid cells. The central macrophages in EBIs provide adhesion molecules and growth factors to promote erythroid cell proliferation and survival, which were all highly expressed in VCAM1+ macrophages (Figure 3E)., Meanwhile, the higher phagocytosis-associated gene signature expression in VCAM1+ macrophages indicated its role in erythroid nucleus phagocytosis (Figure 3C). Importantly, mIF staining of HB showed the island structure of an EBI in situ, in which a VCAM1+ macrophage in the center was surrounded by several erythroid cells (Figure 3F). Spatial analysis also showed a closer spatial co-location between both cells in HB tumors than in adjacent livers (Figure 3G). Together, the above results indicated that VCAM1+ macrophages in HB tumors were EBI central macrophages and formed the specific structure of EBIs in situ, which provided the basic niche for supporting the persistence and proliferation of erythroid cells in HB tumors.

Since erythroid cells also had a strong interaction with DCs, we then set out to explore the particular L-R pairs between erythroid cells and DCs in HB, finding that DCs mainly displayed immunosuppressive interactions with erythroid cells. Among these interactions, the LGALS9/HAVCR2 axis had the highest communication probability for both conventional DC1 (cDC1) and cDC2 (Figure 4A). Next, we asked whether other known immunosuppressive L-R pairs, such as PD-1/PD-L1 or CD86/CTLA-4, might participate. Interestingly, both cDC1 and cDC2 showed low expression of common immunosuppressive receptors, except for HAVCR2 (Figure 4B). At the single-cell level, the immune checkpoint ligand LGALS9 was significantly more highly expressed in HB-specific erythroid cells than in fetal liver-derived erythroid cells (Figures S5A and S5B). Its receptor HAVCR2, also named TIM3,, was broadly expressed in both cDC1 and cDC2 (Figure S5C). Flow cytometry analysis confirmed that most intratumoral erythroid cells exhibited high expression of LGALS9 compared with their counterparts in adjacent liver samples (Figure 4C). On the other hand, we also confirmed that 37.6% of intratumoral DCs expressed a high level of TIM3 (Figure S5C). mIF staining showed that TIM3+CD11c+ DCs resided near LGALS9+ erythroid cells (Figure 4D). Spatial analysis showed that more TIM3+ DCs were located near LGALS9+ erythroid cells in tumors than in paired adjacent liver tissues (Figures 4E and S5E). Cumulative distributions of cross-type distances showed that the spatial organization of TIM3+ DCs and LGALS9+ erythroid cells in the erythroid-high patient represented a clustered pattern (the observed G(r) curve was above the theoretical G(r) curve). In sharp contrast, the spatial organization in the erythroid-low patient and adjacent liver represented a dispersed or random pattern, respectively (Figure 4F). Thus, erythroid cells in HB tumors had stronger cross talk with DCs via the LGALS9/TIM3 axis. LGALS9 was reported to inhibit antigen processing and presentation by bonding to its receptor TIM3 on DCs., To further confirm if erythroid cells had a suppressive capability on DCs, we divided patients into the erythroid-high and the erythroid-low groups based on the median of erythroid cell proportion in tumor samples, which was significantly correlated with the density of erythroid cells in the tumor based on the IHC analysis (Table S1; Figure S5F). We found that DCs from the erythroid-high group displayed lower expression of major histocompatibility complex (MHC)-II and antigen presentation-related genes than DCs from the erythroid-low group and adjacent livers (Figures 4G and S5G), indicating an impaired capability to present antigen. Flow cytometry analysis of DCs from the validation cohort confirmed that intratumoral DCs from the erythroid-high group had lower expression of CD80, CD40, CD86, and HLA-DR (Figures 4H and S6). Further, differentially expressed gene (DEG) analysis revealed that both cDC1s and cDC2s from the erythroid-low group enriched immune activation-related genes (Figure 4I). Conversely, DCs from the erythroid-high group enriched only the cellular physiological function-related pathways (Figure 4J). To exclude the intersample DC gene expression heterogeneity, we compared the fluorescence intensity of functional markers (CD86 and HLA-DR) on TIM3+ DCs that were proximal or distant to LGALS9+ erythroid cells and LGALS9− erythroid cells on the same tumor section. Results showed that TIM3+ DCs near LGALS9+ erythroid cells had a lower intensity of HLA-DR and CD86 than TIM3+ DCs distant from LGALS9+ erythroid cells. In addition, we calculated the ratio of surface marker intensity, dividing marker intensity of TIM3+ DCs that were distant from erythroid cells by marker intensity of those proximal to erythroid cells. Ratios from several patients were compared, demonstrating that LGALS9+ erythroid cells had a stronger inhibitory effect on DCs than LGALS9- erythroid cells (Figure S5H). All of the findings indicated that DCs in HB with higher erythroid cell infiltration displayed a phenotype consistent with an impaired ability to present antigens for activating T cells. Next, we performed the cell co-culture experiment to validate the suppressive function of erythroid cells on DCs in vitro (Figure 5A). We sorted Galectin-9+ and Galectin-9− erythroid cells from fresh HB tumor samples (Figure S7A) and generated monocyte-derived DCs from the peripheral blood of healthy donors, which also expressed TIM3 (Figure S7B). Co-culture assays demonstrated that Galectin-9+ erythroid cells suppressed the expression of antigen presentation-related molecules (HLA-DR, CD40, CD86, and CD80) of DCs, consistent with the above results. Conversely, Galectin-9− erythroid cells had no significant effect on DCs. To investigate the therapeutic potential of targeting the LGALS9/TIM3 axis on erythroid cells in HB, we added anti-TIM3 antibodies into the co-culture system, demonstrating that blocking TIM3 could relieve the inhibitory effect of Galectin-9+ erythroid cells on DCs (Figures 5B and 5C). Thus, our data confirmed that erythroid cells in HB had an immunosuppressive effect on the antigen presentation of DCs via the LGALS9/TIM3 axis. Blocking TIM3 could be a potential treatment strategy to restore the function of intratumoral DCs.

According to the above results, erythroid cells inhibited the antigen presentation of DCs, which might consequently suppress the activation of T cells. To explore the phenotype of T cells in patients with higher erythroid cell infiltration, we subclustered T cells and ILCs into 17 clusters according to their specific genes and T cell receptor (TCR) sequencing results (Figures 6A and S8A). In patients with high erythroid cell infiltration, NK cells were enriched, but CD8+ T cells were significantly reduced, especially CXCL13+CD8+ T cells and GZMK+CD8+ T cells (Figure 6B). To explore the potential function of these CD8+ T cells, we performed gene signature score calculation and pseudotime analysis. Both GZMK+CD8+ T cells and CX3CR1+CD8+ T cells, with higher cytotoxicity scores and clone expansion, were annotated as effector CD8+ T cells (Figures S8B and S8C). CXCL13+CD8+ T cells highly expressed exhaustion-related genes (PDCD1, TIGIT, LAG3) and had modest cytotoxicity scores (Figure S8D), representing exhausted CD8+ T cells., The inferred developmental trajectory based on pseudotime analysis validated the above annotations. Naive T cells and ZNF683+ tissue-resident T cells were located at the beginning of trajectory, while CX3CR1+ effector-memory T cells were positioned at the opposite end of CXCL13+ exhausting T cells, with GZMK+ effector T cells located in between (Figure S8E). These findings indicated that HB tumors with higher erythroid cell infiltration had reduced infiltration of cytotoxic T cells. Next, we compared the CD8+ T cell phenotype in different groups of erythroid cell infiltration. CD8+ T cells from the erythroid-high group had lower effector-memory scores and lower proliferation scores (Figure 6C), as well as smaller amount of clonal T cells (Figure 6D), indicating reduced cytotoxicity and activation. Moreover, interactions necessary for T cell activation between DCs and CD8+ T cells were diminished in the erythroid-high group (Figure 6E), demonstrating our suggestion that erythroid cells inhibited DCs from presenting antigen to CD8+ T cells, leading to CD8+ T cell activation failure. To further validate the above results, we performed flow cytometry analysis in CD8+ T cells from patients with different erythroid cell infiltration levels. In contrast to the erythroid-low group, CD8+ T cells in the erythroid-high group had lower expression of cytotoxic markers, including IFNγ, TNF-α, and granzyme B (Figure 6F), indicating that these cells were in a low cytotoxic state. Thus, all these data demonstrated that HB with higher erythroid cell infiltration had an immunosuppressive microenvironment, characterized as less infiltration of effector CD8+ T cells with impaired cytotoxicity.

The presence of tumor antigens is the foundation of anti-tumor immunity. Although HB is reported to have a limited tumor mutation burden, our study found high expression of tumor-associated antigens, such as GPC3, which hold the potential to induce anti-tumor T cells. Therefore, antigen shortage cannot explain the low response rate of ICBs in HB patients. In the present study, we portrayed a comprehensive single-cell landscape of treatment-naive HB tumors, revealing an immunosuppressive TME, which was mediated by intratumoral EBIs composed of VCAM1+ macrophages and erythroid cells. Surprisingly, our data supported that these intratumoral erythroid cells probably originate in the fetal period and persist in the liver along with HB development. Meanwhile, erythroid cells inhibited the antigen presentation function of DCs via the LGALS9/TIM3 axis, leading to the failure of T cell activation and a worse prognosis for HB. The above evidence demonstrated a pivotal role of intratumoral EBIs in inducing immune evasion and restraining anti-tumor immunity in HB (Figure 7). EBIs are located in the bone marrow, spleen, and fetal liver, but have never been found to have immune-related function in tumors. We observed the specific enrichment of EBIs in treatment-naive HB tumors and investigated the origin of intratumoral erythroid cells. Recent studies found erythroid cells in advanced adult cancers. However, those erythroid cells are quite different from the erythroid cells in HB. Erythroid cells in adult tumors are recruited from the spleen via extramedullary erythropoiesis, while only CD45+ erythroid cells show immunosuppressive function. Moreover, no evidence of EBIs was found in adult tumors.,, In our study, erythroid cells had nearly no expression of CD45, which was different from adult cancers, indicating their different origin. We speculated that erythroid cells persist in the liver during the embryonic period based on the following evidence. Erythroid cells were absent in adjacent liver tissues of patients over 5 weeks of age. There was no significant correlation between the proportion of erythroid cells in tumor and the patients’ serum hemoglobin levels or age. Moreover, the comparison between the transcriptomic profiles of erythroid cells from HB and fetal livers showed a high similarity, which was similar to a recent scRNA-seq study. Unlike in that study, all of our tumor samples were obtained before chemotherapy, which cut down the influence of chemotherapy on hematopoietic function. The presence of EBI structures with VCAM1+ macrophages and surrounding erythroid cells on mIF staining in situ indicated that HB tumors had the essential niches for erythroid cell maturation and differentiation. Our findings also provided key evidence that HB tumor cells are derived from fetal hepatic precursor cells indirectly., Due to their distinct origin and cell type, erythroid cells in HB should have an immunosuppressive mechanism different from that of adult cancers but might be similar among neonates. A previous study proposed that erythroid cells in neonates suppressed immunity in a paracrine manner through the secretion of immunosuppressive factors, which were, however, not identified., We found that erythroid cells could suppress DCs by secreting LGALS9, which is a member of the β-galactoside-binding protein families and binds to its receptor, TIM3. The expression of LGALS9 on erythroid cells and the LGALS9/TIM3 interaction between erythroid cells and DCs were validated by flow cytometry, mIF, and in vitro functional assays. Further experiments were merited to confirm this immunosuppression mechanism in vivo. Targeting the interaction between erythroid cells and DCs could be an underestimated strategy for developing a safe and effective treatment. Our results revealed that the addition of anti-TIM3 antibodies effectively reversed the suppression of LGALS9+ erythroid cells on DCs in vitro, suggesting that TIM3 antibodies could be a promising neoadjuvant therapy for HB. Several clinical trials using anti-TIM3 monotherapy or a combination of anti-TIM3 and anti-PD1 in adult cancer patients are currently underway. In addition, we observed a correlation between erythroid cell infiltration level and prognosis of HB patients, suggesting that erythroid cell infiltration could be a potential biomarker for predicting HB prognosis. In conclusion, this study provided a single-cell atlas of the immune microenvironment in HB, identifying the function of EBIs in mediating an immunosuppressive TME. Our results may assist in the development of safer and more effective immunotherapies for young children with HB.

There are several limitations of this study. First, the number of patient samples used for scRNA-seq analysis was relatively low. It is necessary to confirm the enrichment and function of erythroid cells in a larger cohort. Also, an extended follow-up period is required to obtain sufficient survival data in our cohort. Second, in vivo models are necessary to further investigate the immune-suppressive mechanism of erythroid cells. Third, we suppose that erythroid cells in HB originate in the fetal period. Additional validation in animal models is needed for stronger evidence. Unfortunately, there are no applicable HB models in immature mice yet. Finally, in vivo studies and translational clinical trials need to be conducted to validate and observe the treatment effectiveness of anti-TIM3 antibodies in HB patients.



Requests for resources, reagents and further information regarding this manuscript should be addressed and fulfilled by the lead contact, Sui Peng (pengsui@mail.sysu.edu.cn).

Reagents, resources and materials employed in this study are detailed in the key resources table. Materials and reagents used in this study and generated at our laboratory are available upon request.

A total of 32 patients from the Department of Pediatric Surgery, First Affiliated Hospital of Sun Yat-sen University were enrolled in this study, who were all pathologically diagnosed as HB and had no previous chemotherapy history. Fresh tumor samples (2 × 2 × 1cm), paired 7 non-tumor adjacent liver samples (at least 2 cm from matched tumor tissues), and 3 peripheral blood samples were collected as the discovery cohort, from 13 treated-naïve HB patients. All the samples were obtained at the time of surgery or biopsy. The detailed clinical and pathological information of the Discovery cohort was shown in Table S1. Cell suspensions of pairwise tumor and normal liver samples from 4 patients in the Discovery cohort were frozen to perform flow cytometry further as Validation cohort 1. Frozen tumor tissues (n = 20) and paired normal liver tissues (n = 15) were collected from 20 patients as Validation cohort 2. Formalin-fixed paraffin-embedded (FFPE) tissue blocks of pairwise tumors and normal liver samples were collected from 19 patients as Validation cohort 3. The available clinical features of Validation cohorts are summarized in Table S1. Additionally, we collected the whole peripheral blood from two healthy donors, to isolate monocytes and perform cell co-culture experiments. All patients and donors participating in this study provided written informed consent according to the ethical guidelines of the Institutional Review Board of First Affiliated Hospital of Sun Yat-sen University (No. [2021]492).

PBMCs were isolated from whole blood of two healthy donors using Ficoll-Paque PREMIUM density centrifugation. Then, monocytes were isolated from PBMCs using EasySep Human CD14 Positive Selection Kit II (STEMCELL, cat#17858) according to the manufacturer’s protocol. The isolated monocytes (3×106 cells) were cultured in 3 mL of complete RPMI 1640 medium, supplemented with 20 ng/ml IL-4 and 30 ng/ml GM-CSF, and incubated at 37°C and 5% CO₂ in 6-well plates for 4 days. On day 5, only adherent cells (immature DCs) were harvested. The viability, purity, and expression of TIM3 of immature DCs were determined by flow cytometry (Figure S7B). To generate erythroid cells, tumors were first minced into approximately 1–3 mm3 fragments and then digested into single-cell suspension using Tumor Dissociation Kit (cat# 130-095-929, Miltenyi Biotec). Erythroid cells were sorted using the FACS method and them co-cultured with DC in complete RPMI 1640 medium. The gating strategy for Galectin-9+ erythroid cells, and Galectin-9- erythroid cells was CD45−CD235a+CD71+Galectin-9+, and CD45−CD235a+CD71+Galectin-9-, respectively (Figure S7A). Flow cytometer antibodies are listed in the KEY resource table. Complete RPMI 1640 medium consists of RPMI-1640 (Gibco, cat#C11875500BT) supplemented with 10% fetal bovine serum (Pricella, cat#164210-50), 100 U/ml penicillin and 100 mg/mL streptomycin (Gibco, cat#15140-122).

Fresh tumor tissues and normal liver tissues from HB patients were digested into single-cell suspension using Tumor Dissociation Kit (cat# 130-095-929, Miltenyi Biotec) and Multi-tissue Dissociation Kit (cat# 130-110-201, Miltenyi Biotec) according to the manufacturer’s instructions. Briefly, the tissues were cut into 1-2mm pieces after washing with PBS. The minced tissues were digested in the dissociation enzyme for 30 min at 37°C. The cell suspension was filtered and collected through a 70 μm magnet-activated cell sorting (MACS) SmartStrainer (cat# 130-098-462, Miltenyi Biotec) and a 30 μm MACS SmartStrainer (cat# 130-098-458, Miltenyi Biotec). Peripheral blood mononuclear cells (PBMC) were isolated from whole blood using Ficoll-Paque PREMIUM (cat# 17544203, Cytiva) density centrifugation to get cell suspension. After centrifuging the cell suspension for 6 min at 400×g, red blood cells were lysed using RBC Lysis buffer (cat#00430054, eBioscience) for 5 min at 4°C. Cells were then washed twice with PBS and stained by AOPI for testing the concentration of single-cell suspensions and cell viability. The concentration of cell suspensions was then adjusted to ∼700-1,300 cells/μl to prepare the library construction. Redundant cells were frozen at a defined concentration in FBS containing 10% DMSO, in the case of flow cytometry analysis.

Prepared cells were loaded and converted to barcoded scRNA-seq libraries following the manufacturer’s instructions of “Chromium Next GEM Single Cell 5′ Reagent Kits to capture cells per library” (10x Genomics, v2 chemistry). Libraries were sequenced on Novaseq 6000 (Illumina, San Diego, CA). Raw sequencing reads were aligned to the GRCh38 human reference genome using CellRanger (version 5.0.1) to generate a raw gene expression matrix for each sample. To abrogate the ambient RNA profile present in all droplets, we performed ambient RNA decontamination using R package SoupX v1.5.2. Cells that met the following criterion were regarded as low-quality cells and were removed: (1) total expressed genes <200, (2) proportions of expressed mitochondrial genes of more than 25%. DoubletFinder was applied to each sequencing library to remove potential doublets. The expected doublet rate was set according to the Poisson doublet formation rate, which was calculated based on the cell concentration. In the end, we retained a total of 45,204, 100,883, and 25,752 cells for the non-tumor adjacent liver tissues, tumor samples, and PBMC datasets, respectively. The number of high-quality cells in each sample could be found in Table S2.

Unsupervised clustering was based on the normalized and scaled data using the Seurat R package (version 4.1.0). Briefly, we used the “LogNormalize” method to normalize the feature expression measurements for each cell by the total expression, and then multiply it by a scale factor of 10,000 followed by log-transforming the result. Based on the normalized gene expression matrix, 2,000 highly variable genes were identified by using the ‘FindVariableFeatures’ function with the ‘vst’ method. ‘ScaleData’ function was then used to scale the expression of each gene. Next, we used the function ‘FindClusters’ on 25 PCs with a resolution of 0.2 to perform the first-round cluster and annotated each cluster with known markers. We firstly identified eight major cell types without batch-correction, including 4 immune cell types (ILC/T cells, B cells, myeloid cells, and progenitors) and 5 non-immune cell types (hepatocytes, cholangiocytes, endothelial cells, fibroblasts, and malignant cells). For non-malignant cell sub-clustering, we performed batch effect correction by the Harmony package (version 0.1.0) of R based on the top 30 PCA components identified. The clustering analysis was performed based on the integrated joint embedding produced by Harmony with the Louvain algorithm after computing a shared nearest-neighbor graph with the Louvain algorithm. The original expression data were used in further differential analysis, pseudotime analysis, etc. Finally, to detect cluster-specific expressed genes, the clusters were compared pairwise using the original expression data in “RNA” assay via the Seurat ‘‘FindAllMarkers’’ function to test genes with >0.25-fold difference (log-scale) on average between the two groups of cells and detectable expression in more than 25% of cells in either of the two populations. The marker genes of each cluster were shown in Table S3.

To quantify the tissue distribution of each cell type between tumor and normal liver samples (or between erythroid-high and erythroid-low groups), the observed and expected cell numbers in each cell type were compared according to the following formula as Zhang, described previously:where the expected cell numbers for each combination of cell clusters and tissues are obtained from the Chi-square test.

The chromosomal copy-number variations (CNVs) profile of single cells was inferred using the inferCNV (version 1.8.1) R package. Immune cells, including myeloid cells and T cells, were used as a reference to define a baseline of normal karyotype. The clonal large-scale chromosomal copy number variations of malignant cells were estimated by sorting the analyzed genes by their chromosomal location and applying a moving average to the relative expression values, with a sliding window of 100 genes within each chromosome, as previously described.,

Differentially expressed genes (DEG) analysis between two groups, including erythroid cells from HB tumors and normal liver samples, or DCs from the erythroid cell-high group and erythroid cell-low group, was carried out using the original expression data in “RNA” assay via the ‘FindMarkers’ function in the Seurat package based on Wilcoxon Rank-Sum test. All DEGs (adjusted P threshold of 0.05) were used as input into the ClusterProfiler package v.4.0.5, for gene ontology functional enrichment. Results were visualized using the barplot. The metabolism-related pathways were analyzed via function ‘sc.metabolism.Seurat’ (method = “VISION”, imetabolism.type = “KEGG”) from scMetabolism package (version 0.2.1).

To compare the proliferation activity of erythroid cells, we calculated an S score per cell using the CellCycleScoring function with S cell cycle genes as Seurat R package previously described. The CytoTRACE algorithm was also performed to characterize cell stemness. R package CytoTRACE (version 0.3.3) was applied to calculate the CytoTRACE scores for erythroid cells. CytoTRACE scores range from 0 to 1, while higher scores indicate higher stemness (less differentiation) and vice versa.

Immunocytes gene expression signatures were scored using the “AddModuleScore” function of the Seurat package. The list of functional phenotype genes were shown in Table S5, including fetal liver-derived erythroid cells signature gene score and bone marrow-derived erythroid cells signature gene score for erythroid cells; effector-memory score, exhaustion score, and proliferation score for T cells; M1(pro-inflammatory), M2 (anti-inflammatory), and phagocytosis for macrophages and antigen-presentation score, MHC-II score for DCs. Meanwhile, GOBP dendritic cell antigen processing and_presentation, GOBP MHC class II biosynthetic process, and GOCC MHC class II protein complex genesets were also calculated by same method in DCs. Wilcoxon rank-sum test was used to indicate the statistical significance between each group.

To construct a developmental trajectory of progenitor cells and CD8+ T cells, we applied the Monocle 3 (version 1.0.0) in our expression matrix. Firstly, the data of all the progenitors and corresponding downstream immunocytes were preprocessed by ‘preprocess_cds’ with parameter num_dim = 50 and method = ‘PCA’. Then dimension reduction was performed by function ‘reduce_dimension’ with default parameters. Next, the Seurat UMAP of all the cells was merged and embedded with reduced cds data by function ‘Embeddings’. Finally, the trajectories of cells were constructed using function ‘learn_graph’ and estimated pseudotime using function ‘order_cells’. The HSC/MPP cluster was selected as the root when estimating the pseudotime.

Potential cell-cell interactions between erythroid cells and different immunocytes were identified using CellChat v1.1.3. All the CellChatDB databases were selected to perform cell-cell communication analysis. For the main analyses, the core functions “computeCommunProb”, “computeCommunProbPathway” and “aggregateNet” were applied using standard parameters and fixed randomization seeds. The potential interaction strength between each two cell subsets was predicted based on the expression of ligand-receptor pairs. And the aggregated cell-cell communication network was calculated by counting the number of links. Finally, the function “netVisual_bubble” was used to show all the significant interactions (L-R pairs) from selected cell groups, for example, erythroid cells and macrophages, erythroid cells and DCs.

TCR data was analyzed using scRepertoire (version 1.2.0). The contig_annotation data of different samples were combined into a single list object using the ‘combineTCR’ function. The combined TCR contig list file was then integrated with the corresponding Seurat object of the scRNA-seq data using the function ‘combineExpression’. Only the cells with both TCR and scRNA-seq data were kept for downstream clonotype analysis. The clonotype was defined according to the genes comprising the TCR and the nucleotide sequence of the CDR3 region. The frequency of each clonotype in each patient was then calculated as clone count. Based on the clone count of each clonotype, we divided the clonotypes into 5 subgroups: Hyperextension (clone count >25), Large (10< clone count ≤25), Medium (5< clone count ≤10), Small (1 < clone count ≤ 5), and Single (clone count = 1). The information of BCR for each cell was matched with its mRNA expression data via the unique cellular barcode. We defined the Ig classes according to the c-gene information of each heavy chain. IGHG, IGHA, and IGHE were defined as class-switched Ig (csIg). The composition of Ig classes of each B cell subcluster was calculated to validate their annotation.

Total RNA was collected and then isolated into poly (A) mRNA using Beads (Invitrogen, Carlsbad, CA) with oligo (dT). Fragmentation buffer was added for interrupting mRNA to short fragments. A random hexamer-primer was used to synthesize the first-strand cDNA with the short fragments as templates. The second-strand cDNA was synthesized using buffer, dNTPs, RNase H, and DNA polymerase I, respectively. Short fragments were purified with a QIAQuick PCR extraction kit (Qiagen, Hilden, Germany) and resolved with EB buffer for end reparation and adding poly (A). Then, the short fragments were connected with sequencing adaptors. For amplification with PCR, suitable fragments were selected as templates, concerning the result of agarose gel electrophoresis. At last, RNA libraries were constructed and loaded on the NovaSeq sequencing platform (Illumina), and 150 bp paired-end reads were generated. Qualified reads were obtained after removing raw reads with adapters or of low quality and then aligned to the human genome (hg19) by HISAT (version 2.1.0) with default parameters. RSeQC (version 2.3.4) was used to measure gene expression abundance as reads per kilobase per million mapped reads (RPKM). The P-value of GYPA and TFRC expression between adjacent liver and tumor samples was calculated by paired t-test.

For cell co-culture assays, a fixed number of immature DCs were seeded into 96-well plates individually or together with Galectin-9-/Galectin-9+ erythroid cells at a ratio of 1:2 for 48 hours and then stimulated for another 24 hours with maturation cocktail (2000 IU/mL IL-6, 400 IU/mL IL-1β, 2000 IU/mL TNF-α, and 2 μg/mL PGE₂). In some experiments, cells were treated with 10μg/ml blocking antibodies against TIM3 (Biolegend, cat#345009). After co-culture for 48h together, cells were collected for further flow cytometry analysis.

Flow cytometry and cell sorting were operated by CYTEK Aurora flow cytometry and BD FACSMelody™ Cell Sorter, respectively. For data analysis, Flow Jo v10.6.2 software (Tree Star) was used. Fluorochrome-conjugated antibodies used for the staining are listed in KEY resources table. For frozen tumor cells resuscitation, cells were incubated with 100μg/ml DNase I at 37°C for 15 min to reduce cell adhesion. Target cells were re-suspended in Stain Buffer (BD Bioscience, Cat# 554657). For cell surface staining, cells were incubated for 20 min at RT with fluorochrome-labeled antibodies at a dilution ratio of 1:100 with staining buffer. After washing and fixing, cells were immediately analyzed. For intracellular staining, cells were first stained with antibodies against cell surface markers and then fixed, permeabilized with BD Cytofix/Cytoperm (BD Bioscience, Cat#554722) at 4°C for 20 mins, and stained with corresponding antibody cocktails. For cytotoxic markers analysis of T cells, cells were activated with Leukocyte Activation Cocktail (BD Bioscience, Cat# 550583) for 4 hours first. After sorting CD45-CD71+CD235a+ cells, 700μl of cell suspension containing 105 sorted cells was seeded in a 24-well plate with poly-lysine-coated slides for 24 hours, an then centrifuged for 5 minutes at 500 rpm. The slides were then fixed in methanol for 10 minutes and stained with Giemsa according to the manufacturer’s recommendations (Solarbio, Cat#G1015). Images were acquaired by an Olympus BX63 system at 100× magnification.

Formalin-fixed paraffin-embedded (FFPE) blocks of HB samples were cut into 4 μm slides, deparaffinized in xylene, rehydrated through a decreasing ethanol series, and then fixed in NBF (10% neutral buffered formalin). Slides were stained to enable the simultaneous visualization of four markers on the same slide, using a PANO 7-plex IHC kit (Panovue, Cat# 0004100100). The primary antibodies in the panel showing the interaction of erythroid cell with Macrophage (Panel 1) were: anti-GYPA (Abcam, Cat# ab129024, 1:10000), anti-CD68 (CST, Cat# 76437S, 1:1000), anti-VCAM (Abcam, Cat# ab134047, 1:500) and anti-CD45 (CST, Cat# 13917S, 1:500). The primary antibodies in the panel showing the interaction of erythroid cell with DC (Panel 2) were: anti-Galectin-9 (CST, Cat# 54330S, 1:500), anti-GYPA (Abcam, Cat# ab129024, 1:500), anti-CD11c (Abcam, Cat# ab52632, 1:200), anti-TIM3 (Abcam, Cat# 45208S, 1:100). The primary antibodies in the panel showing the function of DC near or distinct to erythroid cells (Panel 3) were: anti-Galectin-9 (CST, Cat# 54330S, 1:500), anti-GYPA (Abcam, Cat# ab129024, 1:500), anti-CD11c (Abcam, Cat# ab52632, 1:200), anti-TIM3 (Abcam, Cat# 45208S, 1:200), anti-CD86 (CST, Cat# 91882, 1:2000), anti-HLA-DR (Abcam, Cat# ab92511, 1:4000). At each staining cycle, antigen retrieval was performed through microwave-heated treatment in EDTA solution. After blocking proteins for 10 min, corresponding primary antibodies were sequentially incubated at 37°C for 30 min. Then HRP-conjugated secondary antibody and tyramide signal amplification (TSA) with Opal were incubated. Corresponding staining cycle orders and antibodies/fluorescent dye combinations were shown in the KEY resource table. After each TSA operation, the Ab TSA complex was removed by microwave treatment at each cycle of staining. Finally, all slides were stained with DAPI (Panovue, Cat# 0012100500) for 5 min and enclosed with Mounting Medium (Panovue, Cat# 0022001010).

All the stained slides were scanned using the TissueFAXS platform (TissueGnostics, Vienna, Austria) at 20× magnification, which captures the fluorescent spectra at 20-nm wavelength intervals from 420 to 720 nm with identical exposure time; and the scans were combined to build a single stack image. Images of single-stained slides were applied to extract the spectrum of each fluorescein, to establish the spectral libraries which were used to unmix the multispectral images with the StrataQuest software (TissueGnostics, Vienna, Austria) for further analysis. The cut-off value for positivity was determined according to the staining pattern and intensities of all images for each primary antibody. Spatial co-localization analysis of multispectral imaging data was conducted on reconstructed whole slides without introducing edge effects and loss of data. The StrataQuest software (TissueGnostics, Vienna, Austria) was used to calculate the distance between interested cells and the number of mutual neighbors (at the maximum distance of 25μm) for each pair of cell phenotypes. To adjust the density of erythroid cells, we calculated the density of TIM3+ DCs surrounding LGALS9+ erythroid cells within a radius of 25μm, and then divided them by the density of LGALS9+ erythroid cells in tumor. So that, the number of TIM3+ DCs near each LGALS9+ erythroid cell was obtained. Furthermore, cumulative distributions of cross-type distance analysis was performed to describe the spatial organization between erythroid cells and DCs. The function G(r), which is the cumulative distribution of the distance from a typical random cell X to its nearest cell Y, is estimated by R package “spatstat” (version 3.0.2). The argument r is the radius of the area in which G(r) is assessed. The function G(r) and its derived curves can be used to distinguish whether the distribution of cell X and Y was clustered or dispersed spatial patterns when compared with a random point process model (theoretical curve). When the observed G(r) curve was above the theoretical G(r) curve, the spatial organization of the cell phenotype pair may be clustered. Otherwise, the spatial organization may be dispersed. Fluorescence intensity comparison on the same tumor section was performed to show the spatial relationship between LGALS9+ erythroid cells and TIM3+ DCs. The fluorescence intensity of functional markers (CD86 and HLA-DR) on TIM3+ DC, which were proximal or distant to LGALS9+ erythroid cells and LGALS9- erythroid cells, were compared separately. In addition, we calculated the ratio of surface marker intensity, dividing marker intensity of TIM3+ DCs that were distant to erythroid cells by marker intensity of those proximal to erythroid cells. According to the intensity ratio, we can know the extent of DC function inhibition when the DCs were proximal to erythroid cells. Ratios from several patients were compared by paired t-test.

For immunohistochemistry examination, FFPE blocks were sliced into 4m-thick slides. After deparaffinization, the slides were blocked with 10% normal goat serum for 30 minutes, then incubated with the anti-GYPA antibody (Abcam, Cat# ab129024) overnight in a humidified chamber, and subsequently probed with the secondary antibody (DAKO DAB kit). The slides were scanned with a KF-PRO Slide Scanner (Kfbio, China), and images were processed and analyzed with Qupath (version 0.3.0).

The bulk RNA sequencing data from the Affymetrix microarray data from Pavel’s cohort were used to evaluate the prognostic correlation of different immunocyte infiltration. The infiltration level of each cell type was calculated as single-sample gene set enrichment analysis (ssGSEA) score using the gsva function from GSVA R package (version 1.43.1) based on gene sets of each cell type, which were top 300 DEG compared with all the other cell types according to the order of avg_log2FC (Table S4). The signature genes of LGALS9+ erythroid cells were defined based on the combination of canonical marker genes of erythroid cells (GYPA and TFRC) and LGALS9+ erythroid cell cluster specific makers (Table S4). Patients were grouped into high and low expression groups using the cutoff value calculated by the ROC curve for 5-year recurrence-free survival in Pavel’s cohort (survivalROC package, version 1.0.3). The cut-off value was established according to the best combination of sensitivity and false-positive rate (Youden’s index) based on the ROC curves. Then, Kaplan-Meier survival curves were plotted using the R function ggsurvplot to compare the recurrence-free survival between patients from erythroid-high and erythroid-low groups. Meanwhile, univariate Cox regression analysis was performed to calculate the relationship between each celltype score and recurrence-free survival in Pavel’s cohort.

Statistical analyses were performed using R (version 4.1.2). Comparisons between groups were conducted using the Chi-square test or Fisher’s exact test for categorical variables and Student’s t-test or the Wilcoxon rank-sum test for continuous variables. For the experimental data, GraphPad Prism 8 (GraphPad Software) was used to perform statistical analyses and graphics production. A value of P < 0.05 was considered statistically significant.