Integrated single-cell analysis reveals the heterogeneity of tumor-ass

Introduction

Tumor-associated macrophages (TAMs) are one of the main cellular constituents in the colorectal cancer (CRC) microenvironment.^1–3 There are accumulating evidences that TAMs play prominent and multifaceted roles in modulating tumor biology.^1–5 The in-depth exploration of TAMs may benefit the development of TAM-associated therapy strategies.⁵

The invention of single-cell RNA sequencing (scRNA-seq) techniques may revolutionize our novelty knowledge of the CRC microenvironment by providing in-depth insights into its functional and molecular heterogeneity.^6–11 However, these single-cell studies are still limited by small sample sizes, mainly due to the availability and sequencing expenditure of fresh samples. It is difficult for individual single-cell studies with limited cell populations to define small macrophage subtypes and investigate the heterogeneity and plasticity of CRC TAMs. A large sample size will contribute to reveal the complexity of CRC TAMs by enriching the distinct subsets with a small cell subpopulation and by excluding individual variants. In recent years, the development of scRNA-seq processing techniques and strategies, including cell-cell communication evaluation, batch correction removal, and the developmental trajectory, will improve the analytical dimensions and accuracy of scRNA-seq data.^12–14 Although these previous studies have extensively explored the complexity of the CRC microenvironment,^6–11 there is still a lack of the systematic classification, functional analysis, gene signature definitions of CRC TAMs, unveiling their diverse roles in modulating tumor biology.

Therefore, we performed an integrated analysis study of three single-cell datasets with 151 primary CRC tissues and 62 adjacent normal tissues (GSE178341, GSE200997, and GSE231559). We established an accurate classification of TAMs subtypes through unsupervised clustering based on canonical markers, performed functional analysis and the evolutionary trajectory of the TAM subtypes, and explored the cell-to-cell crosstalk of fibroblasts with different TAM subtypes. Finally, we further performed an integrated single-cell atlas of five Asian studies with 52 primary CRC tissues and 27 adjacent normal tissues to confirm these findings of TAM characterization. Concomitantly, the M₂ polarization of TAMs might affect the therapeutic effect of immunotherapy.¹⁵ So we explored the relationship of different TAM subtypes with mismatch repair (MMR) status in the microarray dataset of GSE39582.¹⁶ This systemic study of two integrated single-cell cohorts might offer a profound understanding of the TAM subtypes to aid in the design of TAM-based treatment strategies.

Materials and Methods

Study Cohort

This study utilized eight scRNA-seq datasets (GSE178341, GSE200997, GSE231559, GSE132257, GSE132465, GSE144735, GSE221575, and GSE245552),^6–11 and two bulk sequencing datasets (GSE39582 and GSE17536).^16,17 The transcriptomic and clinical data of these datasets were obtained from the Gene Expression Omnibus (GEO) database. The use of these public datasets was exempt from approval based on the Measures for Ethical Review of Life Science and Medical Research Involving Human Subjects dated February 18, 2023, China. The integrated single-cell cohort of three American studies (GSE178341, GSE200997, and GSE231559) was adopted for the training cohort, including 151 primary CRC tissues and 62 adjacent normal tissues (Supplementary Table S1). The integrated single-cell cohort of five Asian studies (GSE132257, GSE132465, GSE144735, GSE221575, and GSE245552) with was used for the validation cohort, including 52 primary CRC tissues and 27 adjacent normal tissues (Supplementary Table S1). All scRNA-seq datasets were from the 10× Genomics platform.

Single-Cell Data Analysis Process

The single-cell data for each data set were integrated into an aggregate object using the Seurat package (version 5.0.1). The quality control criteria were as follows: Cells with detected genes fewer than 200 were discarded; Cells with detected genes greater than 6000 were further removed for excluding the possible doublets; Cells with >20% mitochondrial content were eliminated; Genes expressed in <3 cells were removed. After quality control, the remaining cells were retained for the subsequent analyses.^18,19

The Log-Normalize method with a scale factor of 10000 was used for the first normalization. The FindVariableFeatures function with the top 2000 variably features were retained for the subsequent analyses. Principal component analysis was applied to reduce the linear dimensionality of the integrated dataset. The ScaleData function was used to regress out the influence of highly expressed genes before principal component analysis. To correct batch effects between different samples, we employed the Harmony method to mitigate these batch effects by the harmony package (version 1.2.0). The FindNeighbors and FindClusters functions were utilized for preliminary clustering and annotation, utilizing the first 20 principal components and a resolution of 0.8. We employed the t-distributed stochastic neighbor embedding (t-SNE) method to visualize these identified clusters. Next, a similar procedure of sub-clustering analysis was applied for further characterizing subpopulations of myeloid cells.

The cell type of each cluster was annotated by its canonical marker genes. The FeaturePlot and Vlnplot functions were applied to depict marker genes for every cell type. We annotated seven major cell types based on their canonical marker genes, including B cells, T/natural killer (NK) cells, myeloid cells, epithelial cells, mast cells, endothelial cells, and fibroblasts (Supplementary Table S2).^18,19 We further annotated eight subsets of myeloid cells based on their canonical marker genes (Supplementary Table S2), including monocytes, neutrophils, CD1C⁺ dendritic cells (named as CD1C-DCs), LAMP3⁺ dendritic cells (named as LAMP3-DCs), LILRA4⁺ dendritic cells (named as LILRA4-DCs), CCL20⁺ TAMs (named as CCL20-TAMs), APOE⁺ TAMs (named as APOE-TAMs), and SLC40A1⁺ TAMs (named as SLC40A1-TAMs).^6–11,18,19

Functional Analysis

Gene Ontology (GO) analysis using the biological process category as the reference gene set was performed to explore the potential functions of these common co-upregulated genes in CCL20⁺ TAMs, APOE⁺ TAMs, and SLC40A1⁺ TAMs (Supplementary Table S3) by the clusterProfiler package (version 4.10.0).²⁰ An adjusted p-value of <0.05 for these enriched pathways was considered to be significant. We adopted the GSVA package (version 1.50.5) for gene set variation analysis (GSVA) on hallmark gene sets among three macrophage subsets (CCL20-TAMs, APOE-TAMs, and SLC40A1-TAMs).

The Developmental Trajectory for Three Macrophage Subsets

We utilized the “monocle” package (version 2.30.0) to delineate the developmental trajectory of three macrophage subsets (CCL20-TAMs, APOE-TAMs, and SLC40A1-TAMs) in the training and validation cohorts, separately.²¹ The DDRTree method was applied to reduce dimensionality and construct the pseudo-temporal order.²¹

Cell-Cell Interaction Analysis

We utilized the CellChat package (version 1.6.1)¹³ to analyze and visualize the interactions between three TAM subsets (CCL20-TAMs, APOE-TAMs, and SLC40A1-TAMs) and other major cell types (CD4⁺ T cell, CD8⁺ T cell, mast cell, B cell, endothelial cell, fibroblasts, and epithelial cell) in the primary tumors tissues of the training and validation cohorts, separately. Putative ligands and receptors were determined according to the expression of each cell type. To accurately estimate the extent of cell-cell interaction, we also performed a random sampling of 1000 cells per population from each macrophage subset and other major cell types.

Determination of a Marker Gene Signature for Each Macrophage Subset

To determine a marker gene signature for each macrophage subset, we utilized the FindAllMarkers function to compare cells within a macrophage subset with other myeloid cell subsets. Differentially expressed genes (DEGs) were defined as having a adjusted p.value of <0.001 and a log2 fold-change value of >1.5, with detectable expression in more than 25% of the cells in the studied subset of both the training and validation cohorts. The top fifteen DEGs were determined by the mean log2 fold-change values of these DEGs in both the training and validation cohorts. The marker gene signatures for each macrophage subset were determined according to these top fifteen DEGs (Supplementary Table S4).

Correlations of Three Macrophage-Gene Signatures with Clinical Outcome

Two microarray datasets of GSE39582 and GSE17536 (Supplementary Table S1) were obtained from the GEO database.^15,16 The expression levels were further normalized and log2-transformed via the limma package. For genes with multiple probe sets, the median value was used. We employed the RobustRankAggreg package to mitigate batch differences between two microarray sets (GSE39582 & GSE17536).²² We respectively calculated the relative abundance of CCL20⁺ TAMs, APOE⁺ TAMs, and SLC40A1⁺ TAMs according to the mean values of their marker genes (Supplementary Table S4) in this integrated microarray cohort. The optimal cutoff values of CCL20⁺ TAM, APOE⁺ TAM, and SLC40A1⁺ TAM abundances were determined according to the MaxStat method.²³ In the microarray dataset of GSE39582,¹⁶ we explored the relationship of TAM subtypes with MMR status, including deficient-mismatch repair (dMMR) and proficient-mismatch repair (pMMR).

Multiplex Immunofluorescence Analysis

The paraffinem-bedded sections with a thickness of 4 μm for 20 CRC patients were adopted for the multiplex Immunofluorescence staining. The CaseViewer System (Danjier Company, China) at 400×magnification were used for obtaining the multispectral images. This study was authorized by the ethics committees of Liuzhou People’s Hospital (Approval No. KY2022-038-02), and conformed to the declaration of Helsinki.

CD68⁺LGALS9⁺ Cells

Multiplex Immunofluorescence staining was performed with the TSAPLus Fluorescent Staining Kit (Dual-Labeling Three-Color, Servicebio, China) according to the manufacturer’s instructions. The tissue sections were dewaxed, rehydrated, and subjected to heat-mediated antigen repair. The sections were incubated with bovine serum albumin at indoor temperature for 30 min, the following primary antibodies of anti-CD68 (rabbit, GB113150, dilution 1:2000, Servicebio, China) overnight at 4°C, the secondary antibody (horseradish peroxidase-labeled anti-Rabbit IgG) at room temperature for 50 min. Subsequently, the sections were processed with TSA plus working solution and were sequentially microwave heat-treated after each TSA treatment. The primary antibody of anti-LGALS9 (rabbit, 17938-1-AP, dilution 1:3000, Proteintech, China) was sequentially applied, followed by incubation with the secondary antibody and TSA treatment. Cell nuclei were stained with 4′-6′-diamidino-2-phenylindole (DAPI, Servicebio, China) after both antigens had been labeled.

CD163⁺STAT6⁺LGALS9⁺ Cells

Multiplex Immunofluorescence staining was performed using the TSAPLus Fluorescent Staining Kit (Three-Labeling Four-Color, Servicebio, China) according to the manufacturer’s instructions. The paraffinem-bedded sections were dewaxed, rehydrated, and subjected to heat-mediated antigen repair. The sections were incubated with Bovine Serum Albumin at indoor temperature for 30 min, the following primary antibodies of anti-CD163 (rabbit, GB113150, dilution 1:2000, Servicebio, China) overnight at 4°C, the secondary antibody (horseradish peroxidase-labeled anti-Rabbit IgG) at room temperature for 50 min. Subsequently, the sections were processed with TSA plus working solution and were sequentially microwave heat-treated after each TSA treatment. The primary antibodies of anti-STAT6 (rabbit, R380957, dilution 1:2000, Zen-bioscience, China) and anti-LGALS9 (rabbit, 17938-1-AP, dilution 1:3000, Proteintech, China) were sequentially applied, followed by incubation with the secondary antibody and TSA treatment. Cell nuclei were stained with DAPI (Servicebio, China) after three antigens had been labeled.

Statistical Analysis

We carried out the statistical analysis with the R program (version 4.2.2). The main analysis was recurrence-free and overall survival (RFS and OS). The Kaplan-Meier estimate and Log rank test were used to compare survival differences between the high- and low-infilation groups (CCL20-TAMs, APOE- TAMs, or SLC40A1-TAMs). We performed the multivariate Cox analysis to identify whether infiltration density (CCL20-TAMs, APOE-TAMs, or SLC40A1-TAMs) was independent of those clinicopathological variables (age, gender, and tumor stage), respectively. A p.value of < 0.05 was considered to be significant.

Results

The Integrated Single-Cell Analysis of Primary CRC Tumors and Adjacent Normal Tissues

To profile the microenvironment landscape of CRC, we performed an integrated single-cell atlas of 151 primary tumors and 62 adjacent normal tissues in the training cohort of three USA studies (Figure 1A). After initial quality control, a total of 321243 cells from these 213 samples were included and integrated based on the harmony algorithm for batch effect removal. After unsupervised clustering, seven major cell types were identified according to their canonical marker genes (Figure 1B), including 100738 T/NK cells, 74410 B cells, 42043 myeloid cells, 3994 mast cells, 6862 endothelial cells, 84421 epithelial cells, and 8775 fibroblasts.

Figure 1 The integrated single-cell atlas of CRC in the training cohort. (A) The t-SNE plot of 321243 cells. (B) Dot plot showing the marker gene expression of seven major cell types. (C) Balloon plot showing the distribution of seven major cell types between primary tumors and adjacent normal tissues. (D) The t-SNE plot of 42043 myeloid cells. (E) Dot plot showing the marker gene expression of nine myeloid cell subsets. (F) Balloon plot showing the distribution of nine myeloid cell subsets between primary tumors and adjacent normal tissues.

We further performed an integrated single-cell atlas of 52 primary tissues and 27 adjacent normal tissues to confirm the microenvironment landscape of CRC in the validation cohort of four Asian studies. (Supplementary Figure S1A). A total of 202390 cells were retained for the subsequent analyses after quality control. After batch correction and clustering, we identified seven major cell types according to their canonical marker genes (Supplementary Figure S1B), including 64665 T/NK cells, 42373 B cells, 20997 myeloid cells, 2695 mast cells, 4467 endothelial cells, 50942 epithelial cells, and 16251 fibroblasts.

In the training cohort, myeloid cells were mainly enriched only in primary CRC tissues, accounting for 17.0% of all detectable cells in tumor tissues but only 3.3% of all detectable cells in adjacent normal tissues (all p<0.05, Figure 1C). We obtained equivalent findings for myeloid cells in the validation cohort (Supplementary Figure S1C).

Heterogeneity and Diversity of Myeloid Cells

We performed a subcluster analysis of 42043 myeloid cells to profile their diversity in the training cohort (Figure 1D). We identified nine myeloid cell subsets according to their canonical marker genes (Figure 1E) and DEGs, including one monocyte subset (6374 monocytes), one neutrophil subset (1822 neutrophils), three dendritic cell subsets (2010 CD1C⁺ dendritic cells, 1546 LAMP3⁺ dendritic cells, and 169 LILRA4⁺ dendritic cells), and three macrophage subsets (6252 CCL20⁺ TAMs, 9536 APOE⁺ TAMs, and 11214 SLC40A1⁺ TAMs) and one undetermined subset (3120 cells).

We further performed a subcluster analysis of 20997 myeloid cells to confirm their diversity (Supplementary Figure S1D) in the validation cohort. We also identified the same nine myeloid subsets according to their canonical marker genes (Supplementary Figure S1E) and DEGs, including one monocyte subset (2110 monocytes), one neutrophil subset (5654 neutrophils), three dendritic cell subsets (1359 CD1C⁺ dendritic cells, 356 LAMP3⁺ dendritic cells, and 460 LILRA4⁺ dendritic cells), three macrophage subsets (3655 CCL20⁺ TAMs, 3885 APOE⁺ TAMs, and 3000 SLC40A1⁺ TAMs) and one undetermined subset (518 cells).

The three TAM subtypes were heterogeneously distributed between primary CRC tumors and adjacent normal tissues in the training cohort (Figure 1F). Importantly, CCL20⁺ TAMs and APOE⁺ TAMs are tumor-specific macrophages in the training cohort, accounting for 15.8% and 24.3% of myeloid cells in tumor tissues but only 2.5% and 1.4% of the myeloid cells in adjacent normal tissues (all p<0.05). In the training cohort, SLC40A1⁺ TAMs were significantly enriched in both primary CRC tumors and adjacent normal tissues. We obtained equivalent findings for three TAM subsets in the validation cohort (Supplementary Figure S1F).

Characteristics of Three Macrophage Subsets

CCL20⁺ TAMs

We identified all co-upregulated genes of CCL20⁺ TAMs in the training and validation cohorts (Supplementary Table S3). We defined a robust gene signature of CCL20⁺ TAMs with these top fifteen upregulated genes (Supplementary Table S4). And CCL20⁺ TAMs expressed relatively higher levels of several chemokines (CCL20, CXCL3, CXCL2, and CCL3), which were known as marker genes of proinflammatory macrophages.

According to the optimal cutoff of CCL20⁺ TAMs (7.716163) in the merged microarray cohort, 756 patients were divided into low- (n=402) and high-infiltration (n=354) subgroups. The Kaplan-Meier curves demonstrated that CCL20⁺ TAM infiltration was not significantly associated with RFS (Log rank test: p=0.094, Supplementary Figure S2A) and OS (Log rank test: p=0.940, Supplementary Figure S2B) in CRC patients. The multivariate COX analysis revealed that CCL20⁺ TAM infiltration may be an independent factor for RFS (hazard ratio [HR]=0.747, 95% confidence interval [CI]=0.560–0.970; p=0.029, Supplementary Figure S2C) but not of OS (HR=0.942, 95% CI=0.734–1.208; p=0.637, Supplementary Figure S2C). GO analysis indicated that these common up-regulated genes of CCL20⁺ TAMs (Shown in Supplementary Table S3) was involved in many inflammation-related pathways (cell migration, chemotaxis and chemokines, Supplementary Figure S2D). The GSVA results suggested that CCL20⁺ TAMs were involved in several activated proinflammatory (inflammatory response) and M₁ polarization (IL-6/JAK/STAT3 and TNFA/NFKB) pathways in both the training (Figure 2A) and validation (Figure 2B) cohorts.

Figure 2 GVSA analysis among three TAM subsets in the training (A) and validation (B) cohorts.

In summary, this single-cell study revealed that CCL20⁺ TAMs might display with proinflammatory and anti-tumor properties.

APOE⁺ TAMs

We identified all co-upregulated genes of APOE⁺ TAMs in the training and validation cohorts (Supplementary Table S3). We defined a robust gene signature of APOE⁺ TAMs according to these top fifteen upregulated genes (Supplementary Table S4). APOE⁺ TAMs expressed higher levels of lipid metabolism genes (APOC1, APOE, LIPA, ACP5, FABP5, and GPNMB), which were the canonical genes of lipid-associated macrophages.

According to the optimal cutoff point of APOE⁺ TAMs (8.032915) in the merged microarray cohort, 756 patients were divided into low- (n=446) and high-infiltration (n=310) subgroups. Kaplan-Meier curves revealed that the high infiltration of APOE⁺ TAMs was significantly associated with poor RFS (Log rank test: p=0.014, Figure 3A) and OS (Log rank test: p=0.042, Figure 3B) in CRC. Multivariate COX model revealed that the infiltration of APOE⁺ TAMs was an independent factor for RFS (HR=1.420, 95% CI=1.096–1.840; p=0.008, Figure 3C) and OS (HR=1.336, 95% CI=1.041–1.716; p=0.023, Figure 3C).

Figure 3 Survival analysis and Functional analysis of APOE+ TAMs. (A) The Kaplan‒Meier curve of recurrence-free survival. (B) The Kaplan‒Meier curve of overall survival. (C) The COX regression analysis of recurrence-free and overall survival. (D) GO analysis.

The GO analysis indicated that these up-regulated genes of APOE⁺ TAMs (Supplementary Table S3) were involved in many lipid-metabolism pathways, including lipid, lipoprotein, sterol, cholesterol, triglyceride, acylglycerol, and steroid metabolic processes (Figure 3D). The GSVA results suggested that APOE⁺ TAMs were involved in several lipid-metabolism pathways (fatty acid metabolism, bile acid metabolism, and adipogenesis) in both the training (Figure 2A) and validation (Figure 2B) cohorts. APOE⁺ TAMs utilized the limited nutrients within the CRC microenvironment by activating metabolism-related pathways and facilitated their anti-infammatory properties by enhancing oxidative phosphorylation.

We observed that tumor-associated fibroblasts secreted the high levels of collagen (COL6A1, COL6A2, COL6A3, COL1A1, COL1A2, and COL4A2) to mainly interact with the highly expressed receptor of CD44 on APOE⁺ TAMs via the COLLAGEN pathway in both the training (Figure 4A–C) and validation (Supplementary Figure S3A–C) cohorts. These findings indicated that the fibroblast-derived collagen pathways might mainly interact with APOE⁺ TAMs to promote the formation of the extracellular matrix in the CRC microenvironment. This might lead to a desmoplastic niche, thus promoting tumor progression.^24–27

Figure 4 Cell communication analysis between different TAM subsets and other cell types in the training cohort. (A) Number of interactions between different cell types. (B) Strength of interactions between different cell types. (C) Dot plot of cell communications between other cell types and TAM subtypes.

In summary, this single-cell study revealed that APOE⁺ TAMs might display with lipid-metabolism and pro-tumor properties. APOE⁺ TAMs might foster a desmoplastic niche by the fibroblast-derived collagen pathways.

SLC40A1⁺ TAMs

We identified all co-upregulated genes of SLC40A1⁺ TAMs in the training and validation cohorts (Shown in Supplementary Table S3). We defined a robust gene signature of SLC40A1⁺ TAMs according to these top fifteen upregulated genes (Shown in Supplementary Table S4). Among five immunosuppressive ligands (PD-L1, PD-L2, LGALS9, FCL1, and PVR), SLC40A1⁺ TAMs only expressed the high level of LGALS9 in the training (Figure 5A) and validation (Figure 5B) cohorts. The multiplex Immunofluorescence staining found that CRC tissues were infiltrated with plenty of LGALS9⁺ TAMs (LGALS9⁺CD68⁺ cells, Figure 5C). LGALS9⁺ TAMs were associated with tumor progression, immune suppression, T cell exhaustion, and immunotherapy resistance.^28–31

Figure 5 Immunosuppressive ligand expression among three TAM subsets. (A) The immunosuppressive ligand expression of three TAM subsets in the training cohort. (B) The immunosuppressive ligand expression of three TAM subsets in the validation cohort. (C) Multiplex immunofluorescence staining of CD68+LGALS9+ cells in CRC tissues.

According to the optimal cutoff point of SLC40A1⁺ TAMs (7.294983), 756 patients were divided into low- (n=636) and high-infiltration (n=120) subgroups in the merged microarray cohort. Kaplan-Meier curves revealed that the high infiltration of SLC40A1⁺ TAMs was significantly associated with poor RFS (Log rank test: p<0.001, Figure 6A) and OS (Log rank test: p=0.010, Figure 6B). Multivariate COX analysis demonstrated that SLC40A1⁺ TAM infiltration was an independent factor for RFS (HR=1.621, 95% CI=1.188–2.214; p=0.002, Figure 6C) and OS (HR=1.500, 95% CI=1.096–2.055; p=0.011, Figure 6C).

Figure 6 Survival analysis and Functional analysis of SLC40A1+ TAMs. (A) The Kaplan‒Meier curve of recurrence-free survival. (B) The Kaplan‒Meier curve of overall survival. (C) The COX regression analysis of recurrence-free and overall survival. (D) GO analysis.

GO analysis revealed that these up-regulated genes of SLC40A1⁺ TAMs (Shown in Supplementary Table S3) were involved in several immunity-related pathways of B cell-mediated immunity, humoral immune response, immunoglobulin-mediated immune response, and circulating immunoglobulin-mediated humoral immune response (Figure 6D). GSVA results suggested that SLC40A1⁺ TAMs were associated with the inhibition of the inflammatory response, interferon-alpha/interferon-gamma response, and M₁ polarization (IL-6/JAK/STAT3 and TNFA/NFKB) pathways in both the training (Figure 2A) and validation (Figure 2B) cohorts.

We revealed that tumor-associated fibroblasts could highly express MIF, which bound to the CD74 receptor on SLC40A1⁺ TAMs in the training (Figure 4A–C) and validation (Supplementary Figure S3A–C) cohorts. The previous studies revealed that the MIF pathways^32,33 were associated with an immunosuppressive microenvironment. And SLC40A1⁺ TAMs highly expressed the immunosuppressive molecule of LGALS9, which bound to the CD45 receptor on CD8⁺ T cells in both training (Supplementary Figure S4A) and validation (Supplementary Figure S4B) cohorts.

In summary, these results demonstrated that SLC40A1⁺ TAMs might exhibit an M₂ phenotype with immunosuppressive properties.

Developmental Trajectory of Three Macrophage Subsets

In the training cohort, pseudo-time analysis revealed that a transition from CCL20⁺ TAMs to SLC40A1⁺ TAMs, with APOE⁺ TAMs located at the intermediate position during M₂ polarization (Supplementary Figure S5A and B). We observed that the marker gene expression of APOE and SLC40A1 gradually increased along this trajectory from CCL20⁺ TAMs to SLC40A1⁺ TAMs (Supplementary Figure S5C). However, we observed that the marker gene expression of CCL20 declined during this trajectory. We obtained equivalent findings for pseudo-time analysis in the validation cohort (Supplementary Figure S5D–F). These results indicated that APOE⁺ TAMs might be the possible transitional position during M₂ polarization.

As shown in Figure 7A, the single-cell data revealed that the expression level of STAT6 (the key molecule of STAT6 pathway), CD163/MRC1 (M₂-macrophage markers), and LGALS9 (an immunosuppressive ligand) gradually increased from CCL20⁺ TAMs to SLC40A1⁺ TAMs. The multiple immunofluorescence staining demonstrated that CD163⁺ TAMs expressed the high level of STAT6 and LGALS9 proteins in CRC tissues (Figure 7B–E). So we speculated that the STAT6 pathway might play a crucial role in M₂ polarization and the expression of LGALS9.^34–36

Figure 7 The expression of STAT6, M2-macrophage markers (CD163 and MRC1), and LGALS9 among different TAM subsets. (A) The transcription level of STAT6, M2-macrophage markers (CD163 and MRC1), and LGALS9 among different TAM subsets. Multiplex immunofluorescence staining results of CD163+ cells (B), STAT6+ cells (C), LGALS9+ cells (D), CD163+STAT6+LGALS9+ cells (E) in primary CRC tissues.

The Relationship of TAM Subtypes and Anti-PD-1 Monotherapy

In the microarray dataset of GSE39582 with MMR data, dMMR CRC had the higher infiltration abundance of CCL20⁺ TAMs (7.873 vs 7.647, p<0.001) and APOE⁺ TAMs (8.254 vs 7.824, p<0.001) than pMMR CRC (Supplementary Figure S6A). But dMMR CRC had similar infiltration abundance of SLC40A1⁺ TAMs (6.994 vs 6.900, p>0.050) with pMMR CRC (Supplementary Figure S6A). According to the previous cutoff of CCL20⁺ and APOE⁺ TAMs, CRC patients were divided into four groups, including the CCL20-TAM^highAPOE-TAM^high, CCL20-TAM^lowAPOE-TAM^high, CCL20-TAM^highAPOE-TAM^low, and CCL20-TAM^lowAPOE-TAM^low groups (Supplementary Figure S6B). And the CCL20-TAM^highAPOE-TAM^high group had the highest prevalence of dMMR phenotype (27.9% vs 16.4% vs 10.6% vs 7.7%, p<0.001). The previous studies demonstrated that dMMR CRC with locally advanced stage was highly sensitive to anti-PD-1 monotherapy.^37–39 The highest prevalence of dMMR phenotype (27.9%) occurred in the CCL20-TAM^highAPOE-TAM^high group, which further supported that the high infiltration of CCL20⁺ and APOE⁺ TAMs could help predict patients’ response to anti-PD-1 immunotherapy.

Discussion

In this single-cell study, we identified three main TAM subtypes, including CCL20⁺ TAMs with proinflammatory and anti-tumor properties, APOE⁺ TAMs with lipid-metabolism and pro-tumor properties, and SLC40A1⁺ TAMs with immunosuppressive and pro-tumor properties. TAMs are traditionally divided into two main subtypes based on their polarization states, including M₁-macrophages with a pro-inflammation phenotype and M₂-macrophages with an immunosuppressive phenotype.^1–5 However, another previous study showed that TAMs could not be fully characterized by the simple dichotomy of M₁– and M₂-macrophages.^1–5 The single-cell technique has emerged as a powerful tool for elucidating the heterogeneity of the CRC microenvironment.^6–11 The invention of scRNA-seq sequencing techniques has revolutionized our novelty knowledge of the CRC microenvironment by offering a detailed understanding of its complexity at the single-cell level.^6–11 However, there is still a notable gap in an accurate subtype system for TAMs. We performed an integrated analysis of three scRNA-seq studies, categorized CRC TAMs into three distinct subsets (CCL20⁺ TAM, APOE⁺ TAM, and SLC40A1⁺ TAM) through unsupervised clustering analysis in the training cohort of three American studies (GSE178341, GSE200997, and GSE231559), and subsequently conducted the GO and GSVA analysis of three TAM subsets to explore their functional properties. And we further obtained the equivalent findings of three TAM subsets in the validation cohort of five Asian studies (GSE132257, GSE132465, GSE144735, GSE221575, and GSE245552). We found the similar TAM subtypes between North American and Asian populations, which implicated that this TAM subtypes had excellent stability. Moreover, our study firstly found that CRC patients with enriched CCL20⁺ and APOE⁺ TAMs were characterized by the highest prevalence of dMMR (27.9%) and might achieve more benefit from the anti-PD-1 immunotherapy.

Among three macrophage subsets, APOE⁺ TAMs might be characterized by a lipid-metabolism phenotype. Our single-cell analysis revealed that APOE⁺ TAMs with the increased expression of lipid metabolism-related genes (APOE, APOC1, LIPA, ACP5, FABP5, and GPNMB) might regulate the process of lipid storage, synthesis, localization and metabolism. Both our single-cell study and the previous study⁴⁰ revealed that APOE⁺ TAMs might facilitate the anti-infammatory property by enhancing oxidative phosphorylation, and lead to a poor prognosis. Some previous studies also revealed that lipid-associated macrophages might promote immune escape and tumor growth by utilizing the limited nutrients of tumor-derived lipids and modulating lipid metabolism.^41–43 In our single-cell study, APOE⁺ TAMs also expressed high levels of SPP1 in both the training (Figure 1E) and validation (Supplementary Figure S1E) cohorts. Both our single-cell study and the previous study⁴⁴ revealed that the high infiltration of SPP1⁺ TAMs exhibited shorter progress-free survival, and the high prevalence of high microsatellite-instability status in CRC. Pseudo-time analysis revealed that APOE⁺ TAMs might be the possible transitional position between CCL20⁺ TAMs and SLC40A1⁺ TAMs, indicating lipid-associated macrophages might be intermediate state during M₂ polarization. It might be a new immunotherapy strategy for the inhibition of M₂ polarization by targeting the TAM-related lipid metabolism reprogramming.

Among three macrophage subsets, our study demonstrated SLC40A1⁺ TAMs were associated with the inhibition of several immune and inflammatory response pathways. And high infiltration of SLC40A1⁺ TAMs was significantly related with worse survival in the integrated microarray cohort of CRC. The previous study revealed that SLC40A1⁺ TAMs were associated with a poor prognosis for cancer by inhibiting a key pro-inflammatory regulator of IL-1B.⁴⁵ In addition, SLC40A1⁺ TAMs were correlated with the invasion and migration of M₂-macrophages and a poor prognosis in cancer patients.^46,47 SLC40A1⁺ TAMs are characterized by increased expression of complement protein C1Q encoding genes (C1QC, C1QB, and C1QA). Recent evidences suggested that C1QC⁺ TAMs might suppress the inflammatory response and promote immunosuppression, leading to a poor prognosis for CRC patients.^48–50 In our study, Selenoprotein P (SEPP1) and Folate Receptor 2 (FOLR-2) were proved to be highly expressed in SLC40A1⁺ TAMs. Recent evidences revealed that SEPP1⁺ TAMs might contribute to promoting M₂ macrophage polarization, anti-inflammatory activation, and tumor recurrence.⁵¹ And FOLR2⁺ TAMs might contribute to an immunosuppressive microenvironment by inducing regulatory T cells activity.⁵² In our study, the sing-cell analysis revealed that SLC40A1⁺ TAMs highly expressed the immunosuppressive ligand of LGALS9. The findings of both single-cell analysis and multiple immunofluorescence staining supported that high expression of LGALS9 in SLC40A1⁺ TAMs might be derived from M₂ polarization by the STAT6 pathway. In summary, SLC40A1⁺ TAMs might exhibit pro-tumor and immunosuppression effects in CRC.

Tumor-associated fibroblasts and macrophages are the most abundant cell components within the microenvironment. Considering the pivotal role of tumor-associated fibroblasts and macrophages, we conducted cell‒cell interaction analysis to explore the potential mechanisms. The analysis revealed that several tumor-associated pathways between fibroblasts and TAMs. In our study, fibroblasts expressed the high levels of collagen (COL6A1, COL6A2, COL6A3, COL4A2, COL1A1, and COL1A2) to mainly interact with the CD44 receptor of APOE⁺ TAMs via the collagen pathways. These results indicated that fibroblast-derived collagen might interact with myeloid cells to prompt the formation of a desmoplastic microenvironment by collagen pathways, leading to unfavorable conditions for immunosuppression.^24–27 This analysis also revealed that fibroblasts secreted MIF, which bound to the receptors of CD74 on SLC40A1⁺ TAMs to promote immune escape. Notably, the fibroblast-macrophage interaction via the MIF axis might promote tumor tumor proliferation, inhibit tumor apoptosis, and drive suppressor cell-mediated immunosuppression.^32,33 Taken together, these results also revealed that SLC40A1⁺ TAMs might play a significant immunosuppressive role in CRC.

Our study has several limitations. Firstly, the classification method for macrophage subtypes in our study was rough and might not fully elucidate the diversity of three TAM subsets. Secondly, the underlying mechanisms of APOE⁺ and SLC40A1⁺ TAMs are not yet fully understood. The further studies are necessary to explore these potential mechanisms which may provide new insights to aid in the development of TAM-related treatment strategies in CRC.

Conclusions

In our study, we integrated the scRNA-seq and bulk RNA-seq data to construct a precise single-cell atlas of TAMs in the CRC microenvironment, providing a profound understanding of the TAM subtypes. Our findings provided the comprehensive and in-depth information about the biological functions of APOE⁺ and SLC40A1⁺ TAMs, with an emphasis on the interplay of APOE⁺ and SLC40A1⁺ TAMs with fibroblasts and T cells. And CRC with enriched CCL20⁺ and APOE⁺ TAMs were characterized by high prevalence of dMMR, which might obtain more benefit from the anti-PD-1 monotherapy. These findings could provide new sights to aid in the development of treatment strategies and to overcome unfavorable conditions in CRC patients.

Abbreviations

TAM, tumor-associated macrophages; CRC, colorectal cancer; scRNA-seq, single-cell RNA sequencing; t-SNE, t-distributed stochastic neighbor embedding; GEO, Gene Expression Omnibus; NK, natural killer; DEG, Differentially expressed gene; GO, gene ontology; GSVA, gene set variation analysis; RFS, recurrence-free survival; OS, overall survival; SEPP1, Selenoprotein P; FOLR-2, Folate Receptor 2.

Data Sharing Statement

The original data are included in our manuscript. Further inquiries should be directed to the corresponding author (Jing Li).

Consent for Publication

All the authors approved the submission of the final manuscript.

Funding

This study was supported by the Science and Technology Plan Project of Liuzhou (No.2022CAC0102), and the Specific Research Project for Research Bases and Talents in Guangxi (No. guikeAD21220090).

Disclosure

Guozeng Xu and Binglan Fang are co-first authors for this study. The authors report no conflicts of interest in this work.

References

1. Munir MT, Kay MK, Kang MH, et al. Tumor-associated macrophages as multifaceted regulators of breast tumor growth. Int J Mol Sci. 2021;22(12):6526. doi:10.3390/ijms22126526

2. Chen D, Zhang X, Li Z, Zhu B. Metabolic regulatory crosstalk between tumor microenvironment and tumor-associated macrophages. Theranostics. 2021;11(3):1016–1030. doi:10.7150/thno.51777

3. Kumari N, Choi SH. Tumor-associated macrophages in cancer: recent advancements in cancer nanoimmunotherapies. J Exp Clin Cancer Res. 2022;41(1):68. doi:10.1186/s13046-022-02272-x

4. Wang H, Tian T, Zhang J. Tumor-associated macrophages in colorectal cancer: from mechanism to therapy and prognosis. Int J Mol Sci. 2021;22(16):8470. doi:10.3390/ijms22168470

5. Song J, Xiao T, Li M, Jia Q. Tumor-associated macrophages: potential therapeutic targets and diagnostic markers in cancer. Pathol Res Pract. 2023;249:154739. doi:10.1016/j.prp.2023.154739

6. Pelka K, Hofree M, Chen JH, et al. Spatially organized multicellular immune hubs in human colorectal cancer. Cell. 2021;184(18):4734–4752. doi:10.1016/j.cell.2021.08.003

7. Khaliq AM, Erdogan C, Kurt Z, et al. Refining colorectal cancer classification and clinical stratification through a single-cell atlas. Genome Biol. 2022;23(1):113. doi:10.1186/s13059-022-02677-z

8. Hsu WH, LaBella KA, Lin Y, et al. Oncogenic KRAS drives lipofibrogenesis to promote angiogenesis and colon cancer progression. Cancer Discov. 2023;13(12):2652–2673. doi:10.1158/2159-8290.CD-22-1467

9. Lee HO, Hong Y, Etlioglu HE, et al. Lineage-dependent gene expression programs influence the immune landscape of colorectal cancer. Nat Genet. 2020;52(6):594–600. doi:10.1038/s41588-020-0636-z

10. Ji L, Fu G, Huang M, et al. scRNA-seq of colorectal cancer shows regional immune atlas with the function of CD20(+) B cells. Cancer Lett. 2024;584:216664. doi:10.1016/j.canlet.2024.216664

11. Liu X, Wang X, Yang Q, et al. Th17 cells secrete TWEAK to trigger epithelial-mesenchymal transition and promote colorectal cancer liver metastasis. Cancer Res. 2024;84(8):1352–1371. doi:10.1158/0008-5472.CAN-23-2123

12. Korsunsky I, Millard N, Fan J, et al. Fast, sensitive and accurate integration of single-cell data with Harmony. Nat Methods. 2019;16(12):1289–1296. doi:10.1038/s41592-019-0619-0

13. Jin S, Guerrero-Juarez CF, Zhang L, et al. Inference and analysis of cell-cell communication using CellChat. Nat Commun. 2021;12(1):1088. doi:10.1038/s41467-021-21246

14. Qiu X, Mao Q, Tang Y, et al. Reversed graph embedding resolves complex single-cell trajectories. Nat Methods. 2017;14(10):979–982. doi:10.1038/nmeth.4402

15. Xiang X, Wang J, Lu D, Xu X. Targeting tumor-associated macrophages to synergize tumor immunotherapy. Signal Transduct Target Ther. 2021;6(1):75. doi:10.1038/s41392-021-00484-9

16. Marisa L, de Reyniès A, Duval A, et al. Gene expression classification of colon cancer into molecular subtypes: characterization, validation, and prognostic value. PLoS Med. 2013;10(5):e1001453. doi:10.1371/journal.pmed

17. Smith JJ, Deane NG, Wu F, et al. Experimentally derived metastasis gene expression profile predicts recurrence and death in patients with colon cancer. Gastroenterology. 2010;138(3):958–968. doi:10.1053/j.gastro.2009.11.005

18. Xie Z, Niu L, Zheng G, et al. Single-cell analysis unveils activation of mast cells in colorectal cancer microenvironment. Cell Biosci. 2023;13(1):217. doi:10.1186/s13578-023-01144-x

19. Zhang L, Li Z, Skrzypczynska KM, et al. Single-cell analyses inform mechanisms of myeloid-targeted therapies in colon cancer. Cell. 2020;181(2):442–459. doi:10.1016/j.cell.2020.03.048

20. Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–287. doi:10.1089/omi.2011.0118

21. Trapnell C, Cacchiarelli D, Grimsby J, et al. The dynamics and regulators of cell fate decisions are revealed by pseudo-temporal ordering of single cells. Nat Biotechnol. 2014;32(4):381–386. doi:10.1038/nbt.2859

22. Kolde R, Laur S, Adler P, Vilo J. Robust rank aggregation for gene list integration and meta-analysis. Bioinformatics. 2012;28(4):573–580. doi:10.1093/bioinformatics/btr709

23. Ogłuszka M, Orzechowska M, Jędroszka D, Witas P, Bednarek AK. Evaluate cutpoints: adaptable continuous data distribution system for determining survival in Kaplan-Meier estimator. Comput Methods Programs Biomed. 2019;177:133–139. doi:10.1016/j.cmpb.2019.05.023

24. Xu Y, Ying L, Lang JK, Hinz B, Zhao R. Modeling mechanical activation of macrophages during pulmonary fibro-genesis for targeted anti-fibrosis therapy. Sci Adv. 2024;10(13):eadj9559. doi:10.1126/sciadv.adj9559

25. Fabre T, Barron AMS, Christensen SM, et al. Identification of a broadly fibrogenic macrophage subset induced by type 3 inflammation. Sci Immunol. 2023;8(82):eadd8945. doi:10.1126/sciimmunol.add8945

26. Tharp KM, Kersten K, Maller O, et al. Tumor-associated macrophages restrict CD8⁺ T cell function through collagen deposition and metabolic reprogramming of the breast cancer microenvironment. Nat Cancer. 2024;5(7):1045–1062. doi:10.1038/s43018-024-00775-4

27. LaRue MM, Parker S, Puccini J, Cammer M, Kimmelman AC, Bar-Sagi D. Metabolic reprogramming of tumor-associated macrophages by collagen turnover promotes fibrosis in pancreatic cancer. Proc Natl Acad Sci USA. 2022;119(16):e2119168119. doi:10.1073/pnas.2119168119

28. Pan Z, Bao L, Lu X, et al. IL2RA+VSIG4+ tumor-associated macrophage is a key subpopulation of the immunosuppressive microenvironment in anaplastic thyroid cancer. Biochim Biophys Acta Mol Basis Dis. 2023;1869(1):166591. doi:10.1016/j.bbadis.2022.166591

29. Daley D, Mani VR, Mohan N, et al. Dectin 1 activation on macrophages by galectin 9 promotes pancreatic carcinoma and peritumoral immune tolerance. Nat Med. 2017;23(5):556–567. doi:10.1038/nm.4314

30. Yan ZX, Dong Y, Qiao N, et al. Cholesterol efflux from C1QB-expressing macrophages is associated with resistance to chimeric antigen receptor T cell therapy in primary refractory diffuse large B cell lymphoma. Nat Commun. 2024;15(1):5183. doi:10.1038/s41467-024-4949

31. Jiang J, Xu Y, Chen D, et al. Pan-cancer analysis of immune checkpoint receptors and ligands in various cells in the tumor immune microenvironment. Aging. 2024;16(15):11683–11728. doi:10.18632/aging.206053

32. Huang WC, Kuo KT, Wang CH, Yeh CT, Wang Y. Cisplatin resistant lung cancer cells promoted M₂ polarization of tumor-associated macrophages via the Src/CD155/MIF functional pathway. J Exp Clin Cancer Res. 2019;38(1):180. doi:10.1186/s13046-019-1166-3

33. de Azevedo RA, Shoshan E, Whang S, et al. MIF inhibition as a strategy for overcoming resistance to immune checkpoint blockade therapy in melanoma. Oncoimmunology. 2020;9(1):1846915. doi:10.1080/2162402X.2020.1846915

34. Li Y, Shen Z, Chai Z, et al. Targeting MS4A4A on tumour-associated macrophages restores CD8+ T-cell-mediated antitumour immunity. Gut. 2023;72:2307–2320. doi:10.1136/gutjnl-2022-329147

35. Zhou S, Zhao T, Chen X, et al. Runx1 deficiency promotes M₂ macrophage polarization through enhancing STAT6 phosphorylation. Inflammation. 2023;46:2241–2253. doi:10.1007/s10753-023-01874-7

36. Li Y, Sheng Q, Zhang C, et al. STAT6 up-regulation amplifies M2 macrophage anti-inflammatory capacity through mesenchymal stem cells. Int Immunopharmacol. 2021;91:107266. doi:10.1016/j.intimp.2020.107266

37. Cercek A, Lumish M, Sinopoli J, et al. PD-1 blockade in mismatch repair-deficient, locally advanced rectal cancer. N Engl J Med. 2022;386(25):2363–2376. doi:10.1056/NEJMoa2201445

38. Chen G, Jin Y, Guan WL, et al. Neoadjuvant PD-1 blockade with sintilimab in mismatch-repair deficient, locally advanced rectal cancer: an open-label, single-centre Phase 2 study. Lancet Gastroenterol Hepatol. 2023;8(5):422–431. doi:10.1016/S2468-1253(22)00439-3

39. Wang QX, Xiao BY, Cheng Y, et al. Anti-PD-1-based immunotherapy as curative-intent treatment in dMMR/MSI-H rectal cancer: a multicentre cohort study. Eur J Cancer. 2022;174:176–184. doi:10.1016/j.ejca.2022.07.016

40. Liu J, Gao M, Yang Z, et al. Macrophages and metabolic reprograming in the tumor microenvironment. Front Oncol. 2022;12:795159. doi:10.3389/fonc.2022.795159

41. O’Neill LA, Kishton RJ, Rathmell J. A guide to immunometabolism for immunologists. Nat Rev Immunol. 2016;16:553–565. doi:10.1038/nri.2016.70

42. Mills EL, Kelly B, O’Neill LAJ. Mitochondria are the powerhouses of immunity. Nat Immunol. 2017;18(5):488–498. doi:10.1038/ni.3704

43. Di Conza G, Tsai CH, Gallart-Ayala H, et al. Tumor-induced reshuffling of lipid composition on the endoplasmic reticulum membrane sustains macrophage survival and pro-tumorigenic activity. Nat Immunol. 2021;22(11):1403–1415. doi:10.1038/s41590-021-01047-4

44. Qi J, Sun H, Zhang Y, et al. Single-cell and spatial analysis reveal interaction of FAP+ fibroblasts and SPP1+ macrophages in colorectal cancer. Nat Commun. 2022;13:1742. doi:10.1038/s41467-022-29366-6

45. Zhang Q, He Y, Luo N, et al. Landscape and dynamics of single immune cells in hepatocellular carcinoma. Cell. 2019;179(4):829–845. doi:10.1016/j.cell.2019.10.003

46. Xue Y, Cheng Z, Liao Y, Chen X. Role of exosome-mediated molecules SNORD91A and SLC40A1 in M2 macrophage polarization and prognosis of ESCC. Discov Oncol. 2023;14(1):177. doi:10.1007/s12672-023-00797-x

47. Zhang X, Sun Y, Ma Y, et al. Tumor-associated M₂ macrophages in the immune microenvironment influence the progression of renal clear cell carcinoma by regulating M₂ macrophage-associated genes. Front Oncol. 2023;13:1157861. doi:10.3389/fonc.2023.1157861

48. Benoit ME, Clarke EV, Morgado P, Fraser DA, Tenner AJ. Complement protein C1q directs macrophage polarization and limits inflammasome activity during the uptake of apoptotic cells. J Immunol. 2012;188(11):5682–5693. doi:10.4049/jimmunol.1103760

49. Deng H, Chen Y, Liu Y, Liu L, Xu R. Complement C1QC as a potential prognostic marker and therapeutic target in colon carcinoma based on single-cell RNA sequencing and immunohistochemical analysis. Bosn J Basic Med Sci. 2022;22(6):912–922. doi:10.17305/bjbms.2022.7309

50. Peng Z, Ren Z, Tong Z, Zhu Y, Zhu Y, Hu K. Interactions between MFAP5⁺ fibroblasts and tumor-infiltrating myeloid cells shape the malignant microenvironment of colorectal cancer. J Transl Med. 2023;21(1):405. doi:10.1186/s12967-023-04281-6

51. Barrett CW, Reddy VK, Short SP, et al. Selenoprotein P influences colitis-induced tumorigenesis by mediating stemness and oxidative damage. J Clin Invest. 2015;125(7):2646–2660. doi:10.1172/JCI76099

52. Xiang C, Zhang M, Shang Z, et al. Single-cell profiling reveals the trajectory of FOLR2-expressing tumor-associated macrophages to regulatory T cells in the progression of lung adenocarcinoma. Cell Death Dis. 2023;14(8):493. doi:10.1038/s41419-023-06021-6