Intratumoral erythroblastic islands restrain anti-tumor immunity in hepatoblastoma

Abstract

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.

Keywords: erythroblastic island; erythroid cells; hepatoblastoma; immunosuppression.

Graphical abstract

Figure 1

scRNA-seq reveals the fetal-like immune landscape of HB (A) Workflow of the experimental strategy. Left: the experimental design of single-cell RNA sequencing (scRNA-seq) as the discovery cohort. Right: validation experiments and corresponding cohorts. (B) Atlas of 146,087 single cells collected from 13 HB tumors and 7 normal adjacent liver tissues (NATs), displaying eight major cell types and 63 subsets. (C) Gene expression heatmap showing the expression levels of marker genes in each cluster. The columns show each subcluster, while the rows show each marker gene. The color of each square indicates the average expression level of each gene in each subcluster. (D) Bar chart showing the proportion of each cell type in 13 tumors and 7 adjacent liver tissues. (E) Bar chart showing the fraction of cells originating from tumors and adjacent liver tissues in each cell type (the dashed line shows the total cell fraction of tumor and adjacent liver tissues). See also Figure S1 and Tables S1, S2, and S3.

Figure 2

Intratumoral erythroid cells originated from fetal livers (A) Uniform manifold approximation and projection (UMAP) plot of identified cell types of erythroid-related lineage, labeled in different colors. (B) Scatterplot showing the correlation between the proportion of erythroid cells (divided by the total cell number) in HB tumors and patients’ age at the time of diagnosis (n = 13). r, correlation coefficient. Spearman’s correlation test. (C) Flow cytometry plots gating on CD71⁺CD235a⁺ erythroid cells from paired HB tumor (left) and normal adjacent liver tissues (right) showing the proportion of erythroid cells in CD45⁻ cells. CD71 (TFRC) and CD235a (GYPA) are both markers of erythroid cells. (D) Representative IHC staining images indicating GYPA⁺ erythroid cell (labeled by red arrows) infiltration into HB tumors and adjacent liver tissues. Scale bar, 50 μm. The right bar plot shows the erythroid cell density in paired tumor and adjacent liver tissues (n = 19). Paired t test; ∗∗∗p < 0.001. (E) The expression levels of two erythroid cell markers (GYPA and TFRC) based on bulk RNA-seq data from paired tumor and adjacent liver tissues of treatment-naive HB patients in the validation cohort (n = 15). Paired t test; ∗p < 0.05, ∗∗∗p < 0.001. (F) Volcano plot showing the differentially expressed genes between erythroid cells in tumors (red dots) and adjacent liver tissues (blue dots). The names of the most significant genes are indicated in the plot. (G) Kaplan-Meier curve showing the recurrence-free survival curves of HB patients characterized by either low (blue) or high (red) erythroid score, which was calculated in the bulk RNA-seq data from Pavel’s cohort. The numbers of patients and the risk classification are indicated. Significance was calculated using the log-rank test. (H) UMAP plot of 12 identified cell clusters from three HB peripheral blood samples in the discovery cohort. (I) Dot plot showing the expression of marker genes and transcription factors of erythropoiesis on erythroid cells from normal adjacent liver tissues (NAT), HB tumors, and peripheral blood. (J) Scatterplot showing the correlation between the proportion of erythroid cells (divided by the total cell number) in HB and patients’ serum hemoglobin level (n = 13). r, correlation coefficient. Spearman’s correlation test. (K) Dot plot showing the expression of erythroid lineage marker genes in erythroid cells at different stages from fetal livers and HB. (L) UMAP showing the integrated erythroid-lineage cells from fetal livers, bone marrow, and HB samples. See also Figures S2 and S3.

Figure 3

VCAM1⁺ macrophages form EBIs with erythroid cells (A) Circle plot showing the crosstalk count between erythroid cells and other immune cell types in HB. The band width is proportional to the interaction numbers. (B) UMAP plot of 13 identified myeloid cell clusters, labeled in different colors. (C) Boxplots showing the expression scores of phagocytosis-associated, M1 phenotype, and M2 phenotype signature genes among macrophage subsets. ANOVA, ∗∗∗∗p < 0.0001. (D) Dot plot showing the significantly (p < 0.05) enriched receptor-ligand interactions between macrophage subsets (purple) and erythroid cells (orange). The color gradient of the dots indicates the communication probability of selected ligand-receptor interactions. (E) Dot plot showing the expression levels of EBI central macrophage marker genes in each macrophage subset. (F) Multiplex immunofluorescence image of an HB tumor showing a representative EBI structure (arrow) formed by CD68⁺VCAM1⁺ macrophages and surrounding CD45⁻GYPA⁺ erythroid cells. Scale bar, 20 μm. (G) Bar plots showing the distance between erythroid cells and VCAM1⁺ macrophages (left) and the ratio of erythroid cells near VCAM1⁺ macrophages (distance ≤25 μm) in HB tumors and adjacent liver samples (right). Paired t test (n = 19), ∗∗p < 0.01, ∗∗∗p < 0.001. See also Figure S4.

Figure 4

Intratumoral erythroid cells suppress the anti-tumor function of DCs via the LGALS9/TIM3 axis (A) Dot plot showing the selected receptor-ligand interactions between erythroid cells (orange) and cDCs (purple). The color gradient of the dots indicates the communication probability of ligand-receptor interactions. (B) Heatmap showing the expression level of selected immunosuppressive ligands on erythroid cells and corresponding receptors on cDC1 or cDC2. (C) Left: flow cytometry plot showing the Galectin-9 (LGALS9) expression in erythroid cells from paired HB tumors (the blue line) and adjacent liver samples (the red line) in the validation cohort. Right: dot plot showing the statistical data. Paired t test (n = 4), ∗p < 0.05. (D) Representative mIF staining image of HB tumors indicating the interaction between erythroid cells (LGALS9⁺GYPA⁺, labeled by red arrowheads) and DCs (CD11c⁺TIM3⁺, labeled by orange arrowheads). Scale bar, 20 μm. (E) Graph showing the number of TIM3⁺ DCs near each LGALS9⁺ erythroid cell (distance ≤25 μm) in HB tumors and adjacent liver sections, calculated by dividing the density of TIM3⁺ DCs near LGALS9⁺ erythroid cells by the density of LGALS9⁺ erythroid cells. Paired t test (n = 19), ∗p < 0.05. (F) Cumulative distributions of cross-type distance from erythroid cells (LGALS9⁺GYPA⁺) to the nearest DCs (CD11c⁺TIM3⁺) evaluated in an area with radius r. The theoretical curve, G_theo(r), with its confidence envelopes illustrates random sample distribution. (G) Boxplots showing the expression of antigen presentation- and MHC II molecule-related gene signatures in DCs from HB tumors with high- and low-erythroid cell filtration and adjacent liver tissues. Wilcoxon test, ∗∗p < 0.01, ∗∗∗p < 0.001, ∗∗∗∗p < 0.0001; ns, not significant. (H) Flow cytometry density plots showing the expression of antigen presentation-related genes in DCs from HB tumors with different erythroid cell infiltration levels in the validation cohort (left). Right: bar plots show mean ± SD of the mean fluorescence intensity (MFI) for depicted antigen presentation-related markers (n = 3). Unpaired t test, ∗p < 0.05; ns, not significant. (I) Volcano plots showing the differentially expressed genes between intratumoral DCs in patients from the erythroid-high group (red dots) and the erythroid-low group (blue dots). The names of the most significant genes are indicated in the plots. (J) Bar plot showing the enriched GO pathways in DCs from patients with high-erythroid cell infiltration, compared with DCs from the erythroid-low group. See also Figures S5 and S6.

Figure 5

LGALS9⁺ erythroid cells inhibited antigen presentation of DCs in vitro (A) The flowchart of functional in vitro cell co-culture experiments. Monocytes were sorted from peripheral blood mononuclear cells and then induced into immature DCs (imDCs). CD45⁻CD71⁺CD235⁺ erythroid cells were sorted from fresh HB tumor samples and sorted into Galectin-9⁺ or Galectin-9⁻ subgroups. Then imDCs were co-cultured with Galectin-9⁺ erythroid cells or Galectin-9⁻ erythroid cells for 48 h, and maturation cocktail was added during co-culture. One of the Galectin-9⁺ erythroid cell groups was treated with anti-TIM3 antibodies. (B) Density plots showing the expression of the depicted antigen presentation-related markers in each co-culture group, labeled with different colors. (C) Graphs showing mean ± SD of the geometric mean fluorescence intensity (gMFI) for the depicted antigen presentation-related markers (n = 3). Unpaired t test; ∗p < 0.05, ∗∗p < 0.01. See also Figure S7.

Figure 6

CD8⁺ T cells in HB with higher erythroid cell infiltration showed low-activated phenotype (A) UMAP plot of 17 identified T cell and ILC clusters. (B) Heatmap displaying the relative enrichment of each T cell and ILC subtype from the erythroid-high group and the erythroid-low group, calculated by the R_o/e score. Chi-square test, ∗p < 0.05, ∗∗p < 0.01, ∗∗∗p < 0.001. (C) Violin plot showing the expression scores of effector memory and proliferation-associated signature genes in CD8⁺ T cells from the erythroid-high group and the erythroid-low group. Wilcoxon test, ∗∗p < 0.01, ∗∗∗∗p < 0.0001. (D) Heatmap displaying the relative enrichment of each CD8⁺ T cell clonotype in HB tumors with high- and low-erythroid cell filtration, calculated by the R_o/e score. Chi-square test, ∗∗∗p < 0.001. (E) Dot plot showing the communication probability of CD80⁻CD28 and CD86⁻CD28 intercellular cross talk between cDCs (orange) and CD8⁺ T cells (purple) from different erythroid cell groups. (F) Flow cytometry plots showing the expression of cytotoxicity-related genes in CD8⁺ T cells from HB tumors with different infiltration levels in the validation cohort. Percentages were calculated over target populations. See also Figure S8.

Figure 7

Schematic diagram of EBI-induced immunosuppressive TME in HB VCAM1⁺ macrophages form EBIs with erythroid cells in HB, which inhibit the antigen presentation function of DCs via the LGALS9/TIM3 axis, leading to the inactivated CD8⁺ T cell anti-tumor response. Targeting TIM3 could be a promising therapy to relieve the immunosuppressive microenvironment.