Spatiotemporal Immune Landscape of Colorectal Cancer Liver Metastasis at Single-Cell Level

Liver metastasis remains a major hurdle to long-lasting survival of patients with colorectal cancer (1, 2), which can be partly explained by the highly dynamic spreading routes of cancer cells (3, 4). Outside the cancer cells, the tumor microenvironment (TME) of liver metastasis harbors a highly immunosuppressive phenotype (5), induces a systemic loss of antigen-specific T lymphocytes (6), and drives the spread of tumors (7). It still remains largely unknown how the immune cells spatially orchestrate colorectal cancer liver metastasis (CRLM) progression and whether the metastatic cellular microenvironment differs from the primary ones (8–13). There is hence a pressing need for recording, identifying, and quantifying cellular or even spatial landscape of CRLM to refresh our understanding of metastasis biology.

Neoadjuvant chemotherapy (NAC) refers to chemical drugs that are administered prior to surgical removal of a tumor. Paradoxical results have been seen for the effect of NAC in CRLM. In some reports, NAC was shown to prolong progression-free survival (14) and overall survival (15). However, especially in the patients with low-risk resectable CRLM, NAC was unexpectedly associated with an unimproved overall survival (16). Another independent group also reported that NAC did not prolong disease-free interval (17). Those seemingly contradictory data raise a critical but poorly understood question of whether NAC could favor antitumor response and improve outcome of patients with CRLM. Thus, dissecting the precise effect of NAC, especially TME alterations, is crucial for understanding the key mechanisms of NAC and designing novel therapeutic strategies.

Here, using single-cell RNA sequencing (scRNA-seq) and spatial transcriptomics (ST) of 97 samples from 24 patients, we aim to explore the immune atlas of CRLM. We revealed that immune microenvironment showed dynamic cellular and spatial changes from primary to metastatic sites. Following NAC treatment, the immune phenotypes underwent antitumor remodeling in responsive patients but shifted toward more immunosuppression in nonresponsive patients, where the immunosuppressive cells were reprogrammed at metastatic sites. Our findings highlight the favorable effect of efficient NAC treatment in patients with resectable CRLM and allow for the data-driven design of novel therapeutic combinations such as NAC plus immunotherapy in selected patients.

To define the single-cell landscape of CRLM, we applied scRNA-seq and ST to quantify CD45⁺ cell dynamics (Fig. 1A; Methods) using paired samples of colorectal cancer, adjacent colon, liver metastasis, and adjacent liver, lymph nodes along colons, and peripheral blood mononuclear cells (PBMC). In total, 24 patients with resectable CRLM were recruited based on strict criteria (Supplementary Table S1). Resectability of the tumors has passed the multidisciplinary team (MDT) discussion composed by senior clinicians as described previously (18). In detail, 89 samples from 20 patients underwent scRNA-seq and 8 samples from 4 patients were sequenced by ST, among which 54 samples from 13 patients were treatment-naïve, 30 samples from 6 patients were treated with neoadjuvant chemotherapy with partial response (PR), and 13 samples from 5 patients were treated with neoadjuvant chemotherapy with progressive disease (PD) or stable disease (SD; Fig. 1A). All the tumors were microsatellite-stable.

After two rounds of quality control and doublet removal of scRNA-seq data (see Methods), 79,703 cells from untreated patients (n = 11), 36,284 cells from NAC-PD/SD patients (n = 5), and 62,643 cells from NAC-PR patients (n = 4) were used for further analysis. On average, each sample contained 2,007 cells, and 1,064 genes (median) were quantified in each cell, with no obvious batch effects across different samples (Supplementary Fig. S1A). We integrated all 178,630 CD45⁺ cells, performed the clustering analysis, and defined the main cell types using the SingleR (19) and cell marker genes (Supplementary Table S2; Methods). Largely consistent with previous scRNA-seq of colon cancer (8, 20), we identified myeloid cells, CD8⁺ T cells, CD4⁺ T cells, natural killer (NK) cells, and B cells from the CRLM samples (Fig. 1B and C; Supplementary Fig. S1A–S1C). Of note, T-cell proportion was obviously distinct across different tissue types, suggesting the tissue heterogeneity of the immune microenvironment (Fig. 1D). Both colorectal cancer and paired liver metastasis showed ∼10% regulatory T cells (Treg) of all CD45⁺ cells, whereas such cells were scarce in the adjacent normal tissues (Supplementary Fig. S1D and S1E). This indicated a suppressed intratumor immunity, which was previously ascribed to the reduced immunotherapeutic efficacy in liver metastasis (21). Likewise, we integrated the SingleR (19) and manual marker-based annotation to further define immune cell subsets (Fig. 1E; Supplementary Fig. S1F; Supplementary Table S2; see Methods) for further in-depth analysis.

In parallel, ST data of eight samples (paired colorectal cancer and liver metastasis) from two untreated patients and two NAC-PR patients (Fig. 1F and G; Supplementary Fig. S1G and S1H) were analyzed. Unsupervised clustering analysis classified the samples into distinct regions such as tumor and fibroblast regions (Fig. 1F). To integrate the scRNA-seq and ST data, we used Seurat (22) to quantify the main immune subpopulations. In untreated samples, we observed the immune cell pattern consistent with scRNA-seq, such as CD8⁺ T cells' remarkable enrichment in both liver metastasis and colorectal cancer (Fig. 1G). By contrast, in NAC-treated samples, we observed smaller tumor regions and irregular distribution of immune cells. Especially in liver metastasis, only a small fraction of spots was clustered into cancer cells, indicating the fundamental remodeling of cellular compartment after NAC treatment. Taken together, the results highlighted the spatiotemporal dynamic immune cell landscape of patients with CRLM and the huge remodeling after NAC. The single-cell CRLM atlas is now available online for exploring the immune diversity in primary and metastatic sites (http://www.cancerdiversity.asia/scCRLM).

Given the heterogeneous cellular compositions across patients, we utilized a bias-corrected and accelerated (BCa) bootstrap algorithm to estimate the cell proportion that uses bootstrap resampling for 1,000 times to estimate the SEs associated with cell-type proportion estimations (23). We filtered low-abundant immune cells, clustered the immune cell proportion, and scaled them across tissues (Fig. 2A; Supplementary Fig. S2A). In general, immune cell abundance could be classified into six distinct patterns in untreated samples (adjacent colorectal high, colorectal cancer high, adjacent liver high, LM high, colorectal high, and liver high). Adjacent colorectal samples showed the highest proportions of AIM2⁺ memory B cells, TCL1A⁺-naïve B cells, and CD4⁺-naïve T cells, whereas colorectal cancer samples were mostly enriched with FOXP3⁺ Tregs. MAIT cells, NK cells, and FGFBP2⁺ GZMB⁺ CD8⁺ T cells were enriched in adjacent normal liver tissues, analogous to our previous reports in primary liver cancer (24, 25). These results indicated an organ-specific distribution of certain immune cells. Interestingly, PDCD1⁺ CD4⁺ T cells and CTLA4⁺ CD8⁺ T cells were relatively enriched in primary tumors (2.24% and 1.49% of CD45⁺ cells) and liver metastases (1.17% and 3.02% of CD45⁺ cells), respectively, which is in accordance with other computational predictions (26). The co-occurrence of these cells in the TME indicated a common tumor education on immune cells. In contrast, some immunosuppressive cells were specifically present in liver metastasis, among which a sharp increase of SPP1⁺ macrophages and MRC1⁺ CCL18⁺ macrophages was observed (Fig. 2B). We asked the origin of such liver metastasis–specific macrophages and found that especially the MRC1⁺ CCL18⁺ macrophages shared the Kupffer cell signature (27), indicating their potential origin from hepatic Kupffer cells (Supplementary Fig. S2B; Methods). Neutrophils, which were reported to function as potential protumor players (28), were also enriched in liver metastasis (Fig. 2B). We further used single-sample gene set enrichment analysis (ssGSEA), a method proven to be robust for ST data (29), to quantify the spatial distribution of the immune cell subsets. Again, MRC1⁺ CCL18⁺ macrophages and SPP1⁺ macrophages were enriched in treatment-naïve liver metastasis (Fig. 2C). Those data were partly consistent with the scRNA-seq data, further confirming the potentially suppressive TME at liver metastasis.

We next sought to validate our findings in another three independent CRLM cohorts. First, we used multiplex IHC (mIHC) to quantify spatial distribution of selected immune cell subsets in additional 27 patients with CRLM (n = 18, untreated; n = 9, NAC-PR) in the same hospital (Fig. 2D–F; Supplementary Table S1). Among the untreated samples, we observed the significantly increased infiltration of FOXP3⁺ Tregs, MRC1⁺ CCL18⁺ macrophages, and SPP1⁺ macrophages in colorectal cancer and liver metastasis, largely consistent with our scRNA-seq data (Fig. 2F). Second, we utilized ssGSEA to score the immune cells in a bulk microarray data set of CRLM (ref. 30; n = 133; see Methods). Because of lack of adjacent tissue data, we could only check the cell abundance between primary and metastatic tumors. In accordance with our scRNA-seq data, MRC1⁺ CCL18⁺ macrophages, SPP1⁺ macrophages, and neutrophils were all significantly higher in metastatic tumors in this validation cohort (Supplementary Fig. S2C). Third, we explored the prognostic value of scRNA-seq defined immune cell subpopulations in The Cancer Genome Atlas (TCGA) colorectal cancer cohort (ref. 31; n = 430, Fig. 2G; Supplementary Fig. S2D). The high scoring of MRC1⁺ CCL18⁺ macrophages and SPP1⁺ macrophages in primary tumors both predicted a significantly worse prognosis (Fig. 2G; Supplementary Fig. S2D). To summarize, our single-cell and spatial immune landscape of CRLM could be validated in independent cohorts, indicating the robustness of our data. These observations also led us to hypothesize that specific macrophage subpopulations may play fundamental roles in protumor niche formation in CRLM.

Given the enrichment of macrophage subpopulations in liver metastasis, we hypothesized that those liver metastasis–enriched macrophages might be functionally distinct across different tissues. We hence performed the ligand–receptor-based cancer–immune cross-talk analysis and compared between primary and metastatic sites (Fig. 3A; Methods). Strikingly, SPP1⁺ macrophages and MRC1⁺ CCL18⁺ macrophages ranked high among all differential cross-talking cell types between liver metastasis and colorectal cancer. This result suggested that liver metastatic cancer cells may preferentially reprogram macrophages and induce their specific functional status (Fig. 3A), which is likely due to the intrinsic difference between primary and metastatic cancer cells (Supplementary Fig. S3A). Of note, metastatic tumor cells in liver metastasis showed preferential expression of ligand CD47, which is an important “don't-eat-me” checkpoint (32, 33), and therefore could potentially recruit or activate MRC1⁺ CCL18⁺ macrophages through the corresponding receptor SIRPA (Fig. 3B). We further utilized mIHC to visualize those proteins and directly observed the cross-talk between CD47⁺ tumor cells and SIRPA⁺ macrophages (Fig. 3C). These ligand–receptor pairs specifically enriched in liver metastasis provided clues for targeted therapy against liver metastasis.

To pinpoint the difference imprinted by unique TME, we performed differential gene-expression analysis between colorectal cancer and liver metastasis macrophages. Notably, liver metastasis–enriched MRC1⁺ CCL18⁺ macrophages highly expressed a broad spectrum of key molecules fundamental for macrophage polarization (Fig. 3D; Supplementary Table S3). APOE, which functions as an anti-inflammatory and pro-M2 conversion protein (34), was one of the significantly differentially expressed genes. Other M2 polarization–related genes such as MARCO (35) were also significantly upregulated in liver metastasis–enriched MRC1⁺ CCL18⁺ macrophages (Fig. 3E). By contrast, we found colorectal cancer–enriched MRC1⁺ CCL18⁺ macrophages showed higher expression of inflammatory cytokines including TNF, IL1B, CCL3, and CCL4. Such shift of molecular profiles was also observed in SPP1⁺ macrophages, indicating that macrophages shared common functional changes from primary to metastatic sites (Supplementary Fig. S3B and S3C; Supplementary Table S3). Then, pathway enrichment analyses showed (Fig. 3F) a strong enrichment of metabolic pathways in MRC1⁺ CCL18⁺ macrophages from both colorectal cancer and liver metastasis, highlighting the fundamental role of metabolism in determining those macrophage functions. Notably, we found that MRC1⁺ CCL18⁺ macrophages in colorectal cancer were enriched with oxidative phosphorylation, whereas their counterparts in liver metastasis were dominated by amino acid metabolism. These data highlighted that metabolic regulation might mediate the phenotypic and functional shift of macrophages in response to distinct microenvironmental cues.

Mining metabolic activity at single-cell resolution remains challenging, mostly due to lack of an easy-to-use computational pipeline. Although Locasale Lab recently reported their ssGSEA-based workflow (36), performing such analysis still requires a strong computational background and is time-consuming. Here, we present an R package, scMetabolism, which integrates dropout gene imputation, metabolism quantification, and data visualization. We demonstrated that this workflow shows good compatibility with Seurat scRNA-seq toolkit (22, 37) and is broadly applicable to scRNA-seq data sets. The R package of scMetabolism is now available on Github (https://github.com/wu-yc/scMetabolism).

The scMetabolism package consists of three serial steps (Fig. 4A). First, scMetabolism supports the imputation of the dropout value based on ALRA (38). Second, scMetabolism quantifies the metabolic pathway activity based on VISION (39), AUCell (40), and ssGSEA (41). We integrated reported gene sets (36) and manually reviewed gene sets from KEGG (42, 43) and REACTOME (43, 44) to generate a list of high-quality metabolic gene sets. Third, scMetabolism provides three different ways, including box plot, uniform manifold approximation and projection (UMAP)/t-distributed stochastic neighbor embedding plot, and dot plot, to visualize the metabolic diversity of single cells (see Methods). We also tested the accuracy (comparison with published methods; ref. 36), consistency (correlation between different methods such as VISION and AUCell), and stability (dropout simulation) of scMetabolism on four different data sets, including PBMCs, melanoma, head and neck squamous cell carcinoma, and testis (45–47), to test the performance of this pipeline (Supplementary Fig. S4A–S4E; see Methods). Finally, we developed an easy-to-use web server, scMetabolism online (http://www.cancerdiversity.asia/scMetabolism/), for biologists to easily quantify metabolic activity of single cells.

To understand the metabolic landscape of myeloid cells in CRLM, we first computed the score of all 76 active metabolic pathways in suppressive macrophages and monocytes. Among them, we selected top 50% variable metabolic scores and ranked myeloid cells (excluding myeloid cells with low abundance) according to their average metabolic score (Fig. 4B). Strikingly, liver metastasis–infiltrated MRC1⁺ CCL18⁺ macrophages showed the highest metabolic activities among all myeloid cells, indicating that they were extremely energetic, which could be linked to their unique functions at metastatic sites.

Further unsupervised clustering of all metabolism-associated genes (see Methods) of MRC1⁺ CCL18⁺ macrophages observed a strong metabolic preference associated with surrounding tissues (Fig. 4C), identifying 42 metabolism pathways potentially upregulated in colorectal cancer and liver metastasis–infiltrated MRC1⁺ CCL18⁺ macrophages. Among them, phenylalanine metabolism, which aids in generating tyrosine (48), was more highly expressed in liver metastasis–enriched MRC1⁺ CCL18⁺ macrophages. We subsequently closely looked at the differentially expressed metabolic genes and found metabolic genes specifically upregulated in primary and metastatic infiltrated MRC1⁺ CCL18⁺ macrophages, which were partly the target genes of known drugs (fetched from Enrichr; Fig. 4C; Methods; ref. 49). GSTO1, an enzyme in metabolism of xenobiotics by cytochrome P450 (50), was significantly higher in colorectal cancer and liver metastasis–infiltrated MRC1⁺ CCL18⁺ macrophages and associated with differentiation (Fig. 4D and E). IL4I1 and MIF, genes essential for macrophage activation (51–53), also showed similar features (Fig. 4D and E). We validated those metabolic enzymes by using mIHC and observed largely consistent results (Fig. 4F; Supplementary Fig. S4F and S4G).

We next asked whether such highly metabolic activated states of macrophages were conserved across species. We hence constructed the mouse model of CRLM (Fig. 4G; Methods) and applied scRNA-seq on the primary/metastatic tumors (Supplementary Fig. S4H and S4I). As a result, similar macrophage subpopulations, such as Mrc1⁺ macrophages and Spp1⁺ macrophages, were identified in the mouse CRLM model (Fig. 4H). However, we did not observe the significantly differential proportion of these macrophages between primary and metastatic sites, possibly due to the fact that primary tumors were established subcutaneously rather than in the colon (Supplementary Fig. S4J). Further metabolic quantification showed the featured metabolic states enriched in liver metastasis and colorectal cancer–infiltrated macrophages (Fig. 4I). We checked the top three differentially expressed pathways (liver metastasis vs. colorectal cancer in the mouse model), and some of them were conserved in human and mouse. For example, oxidative phosphorylation and phenylalanine metabolism were in consistency between human and mouse macrophages. The expression of MIF, the key enzyme of phenylalanine metabolism, is also consistent with human CRLM (Supplementary Fig. S4K). In summary, macrophages in liver metastasis harbored a sharp increase of metabolic activity such as phenylalanine metabolism. Inhibiting such metabolism activity might mobilize an antitumor immune response to control hepatic metastases in advanced colorectal cancer.

To illustrate the association of metabolic reprogramming with phenotypic changes of macrophages, we computed the classic phenotypic scores of macrophages such as antigen processing and presentation, inflammation, angiogenesis, immune suppression, phagocytosis, interleukin signaling pathway, and complement (Fig. 4J; Supplementary Fig. S4L and S4M; Supplementary Table S4). Indeed, within liver metastasis, MRC1⁺ CCL18⁺ macrophages harbored significantly increased scores of antigen processing and presentation and complement activity (Fig. 4J; Supplementary Fig. S4N), and SPP1⁺ macrophages also showed similar trend (Supplementary Fig. S4L). Certain genes like HLA-A (Supplementary Fig. S4M), one of the MHC class I molecules (54), were upregulated in liver metastasis–enriched MRC1⁺ CCL18⁺ macrophages. We further analyzed the correlation between the 42 tumor-high metabolic pathways and phenotypic scores of MRC1⁺ CCL18⁺ macrophages. Interestingly, several macrophage-specific metabolic pathways, including pentose and glucuronate interconversions, showed strong correlations with antigen processing and presentation in liver metastasis but weak in colorectal cancer (Fig. 4K; Supplementary Fig. S4O). These data collectively suggested the potentially essential role of metabolism in determining MRC1⁺ CCL18⁺ macrophages' phenotypes.

To explore the effect of NAC on antitumor immunity, we performed the clustering analysis of all CD45⁺ cells using TooManyCells (55). Interestingly, we observed strong difference of cell clustering between colorectal cancer and liver metastasis (Fig. 5A; Supplementary Fig. S5A). To understand which cell types were changed, we computed the bootstrapped cell proportion and performed the differential analysis (Fig. 5B and C; see Methods). Of note, the impact of NAC on the TME was highly cell type–specific and context-dependent (Fig. 5B and C). Generally, in PR tumors, NAC could downregulate immunosuppressive cells and activate cytotoxic cells in tumors (Fig. 5B). In liver metastasis, responsive NAC downregulated the SPP1⁺ macrophages but upregulated cytotoxic T cells such as GZMK⁺ CD8⁺ T cells and XCL1⁺ CD8⁺ T cells (Fig. 5B and D). Particularly, several MRC1⁺ CCL18⁺ macrophages and dozens of SPP1⁺ macrophages were left in NAC-PR tumors, suggesting a potential role of NAC in inhibiting macrophage lineages. However, in PD/SD tumors, we observed distinct immune cell changes form PR tumors. For example, an increased level of SPP1⁺ macrophages and MRC1⁺ CCL18⁺ macrophages was infiltrated in liver metastasis, whereas the cytotoxic immune cells (i.e., FGFBP2⁺ GZMB⁺ CD8⁺ T cells) were downregulated (Fig. 5C, E, and F). Those opposite results might be explained by the organ-specific and treatment efficacy–related remodeling of cancer cells. For example, metastatic PD/SD cancer cells showed the specific enrichment of DNA repair and fatty acid metabolism, whereas metastatic PR cancer cells were more associated with adipogenesis (Supplementary Fig. S5B). Alternatively, NAC treatment may also directly and differentially regulate individual immune subsets. In addition to tumors, NAC upregulated the proportion of GZMB⁺ CD8 T cells in adjacent normal liver and increased the proportion of FGFBP2⁺ GZMB⁺ CD8⁺ T cells in PBMC (Supplementary Fig. S5C and S5D) of PR patients. As for PD/SD patients, a reduced level of FGFBP2⁺ GZMB⁺ CD8⁺ T cells was also observed in adjacent colon and liver (Supplementary Fig. S5D). Altogether, the results demonstrated that effective NAC not only reprogrammed intratumor immune balance but also activated systemic antitumor immune responses, providing evidence to support the potential use of NAC in resectable CRLM.

Considering that NAC mostly affected myeloid populations, especially macrophages, we hence asked whether the trajectory of macrophages was edited by NAC. We split the UMAP distribution of myeloid cells according to receiving NAC treatment response and identified a rich diversity of myeloid phenotypes between PD/SD and PR samples (Fig. 6A). In untreated samples, we could observe the developing trajectory of myeloid cells, whereas myeloid cells in the PR group were enriched with a monocyte phenotype, a progenitor state of macrophages (56), indicating an early differentiation status (Fig. 6A). The bootstrap cell proportion differential analysis confirmed the lack of terminal status of macrophages in NAC-PR tumors, such as a significant downregulation of MRC1⁺ CCL18⁺ macrophages and SPP1⁺ macrophages (Fig. 6B; Supplementary Fig. S6A). By contrast, the NAC-PD/SD tumors harbored a significantly increased level of MRC1⁺ CCL18⁺ macrophages and SPP1⁺ macrophages compared with treatment-naïve tumors (Fig. 6B; Supplementary Fig. S6A). We next utilized mIHC to check those macrophages in another independent cohort (n = 18, untreated; n = 9, NAC-PR; Fig. 6C; Supplementary Table S1), confirming the lower infiltration of MRC1⁺ CCL18⁺ macrophages and SPP1⁺ macrophages in NAC-PR liver metastases (Fig. 6D; Supplementary Fig. S6B). Trajectory inferring analysis using dynverse (57) further identified a clear developing trajectory from S100A9⁺ monocytes to MRC1⁺ CCL18⁺ macrophages, indicating the opposite polarization of untreated versus NAC-PR cells (Fig. 6E). By focusing on MRC1 and S100A9 expression, which represents MRC1⁺ CCL18⁺ macrophages and monocytes, respectively, we again demonstrated that the chemotherapy-responsive cells exhibited monocyte-enriched phenotypes. On the contrary, in the PD/SD group, the trajectory reprogramming of macrophages was not observed (Supplementary Fig. S6C). The results indicate that the trajectory blocking of macrophages by chemotherapy is exclusive in responsive tumors.

Given the limited number of MRC1⁺ CCL18⁺ and SPP1⁺ macrophages in metastatic tumors, we performed the pathway enrichment analysis of all myeloid cells between untreated and NAC-treated samples and found that NAC induced specific functional changes of myeloid cells in the PR and PD/SD groups (Fig. 6F; Supplementary Fig. S6D). Liver metastasis macrophages in the PR group showed the enrichment of ribosome gene signature, whereas those in the PD/SD group were associated with complement and coagulation cascades. Then, the scMetabolism pipeline, focusing on metabolic changes, showed that NAC considerably blocked myeloid metabolism and the PR group showed more downregulated metabolism (Fig. 6G). The glycolysis and gluconeogenesis were consistently downregulated in macrophages from the PD/SD and PR groups in both colorectal cancer and liver metastasis. By contrast, macrophages in the PD/SD group were also associated with specific metabolic pathways such as inositol phosphate metabolism upregulation in colorectal cancer and liver metastasis. Together, NAC exerted specific metabolic reprogramming of myeloid lineage, which may extend our understanding of the NAC-induced antitumor mechanisms and informing potential treatment for patients with resectable CRLM.

A significant proportion of tumors within an organ are multifocal tumors, which harbor diverse genetic, molecular, and immunologic features (58–62), as demonstrated by the recently reported heterogeneous TME of multifocal primary liver cancer (60, 63, 64). Whether multifocal liver metastasis harbors a diverse immune microenvironment and distinct response to chemotherapy remains largely unclear. To generally capture intermetastatic immune heterogeneity, we performed the unsupervised clustering of all CD45⁺ cells of four patients with multifocal liver metastasis (P5, P7, P9, and P10). Interestingly, multifocal metastases showed distinct immune profiles (Fig. 7A). As for the untreated patient (P10), liver metastasis site 2 was in the same cluster of adjacent liver, whereas liver metastasis site 1 was not similar to any other samples. Another three NAC-treated multifocal liver metastases were also featured with intermetastatic heterogeneous immune profiles. Next, bootstrap cell proportion analysis showed that, in untreated samples, immunosuppressive cells were significantly and differentially infiltrated between different liver metastasis sites, where the MRC1⁺ CCL18⁺ macrophages were differentially infiltrated in liver metastasis of P10 (Fig. 7B). They also harbored distinct intermetastatic phenotypic/metabolic activities such as antigen processing and presentation (Supplementary Fig. S7A and S7B). Thus, such considerable intermetastatic heterogeneity may result in differential responses to chemotherapy and other treatments, which is well documented in clinical practice (65, 66).

Here, we defined the cellular and spatial immune landscape of CRLM by using the state-of-the-art high-throughput scRNA-seq and ST. It is known that the genomic divergence between the primary and metastatic colorectal cancer cells is relatively low (67–69). However, our data revealed that the immune microenvironment has undergone extensive spatiotemporal remodeling in liver metastasis with a strong enrichment of immunosuppressive cells (Fig. 7C), consistent with previous bulk sample–based studies (5, 26, 70). In particular, MRC1⁺ CCL18⁺ macrophages showed a terminally differentiated state, exhibited a metabolically energetic phenotype, and were sensitive to neoadjuvant chemotherapy. Our study not only unravels a previously unrecognized diversity of CRLM immune ecosystem, but also provides a comprehensive resource for the cancer research community.

By applying scMetabolism, a novel R package for quantifying the single-cell metabolic activity on our CRLM data, we surprisingly found that the suppressive MRC1⁺ CCL18⁺ macrophages displayed the highest metabolic activity, especially in the metastatic tumors, suggesting that intratumor immunosuppressive cells were highly energetic. However, after chemotherapy, metabolism of such immunosuppressive cells was significantly slowed down. For example, alanine, aspartate, and glutamate metabolism of MRC1⁺ CCL18⁺ macrophage cells, which is essential for macrophage development and fate decision (71), was downregulated by NAC. Our scMetabolism pipeline not only reveals the metabolic landscape of CRLM immune microenvironment, but also opens new computational opportunities for discovering functional metabolic checkpoints and pathways for other scRNA-seq data sets.

Compared with NAC-associated immune activation (increased infiltration of NK cells and oligoclonal expansion of T cells) in advanced ovarian cancer (72), we surprisingly observed that the cellular states of macrophages were reprogrammed in NAC-responsive patients. In PR samples, there existed only a small fraction of terminally suppressive macrophages such as MRC1⁺ CCL18⁺ macrophages and SPP1⁺ macrophages, whereas an opposite trend was observed in PD/SD samples. Mechanically, one possibility is that such immune reprogramming was induced by the chemotherapeutic modulation of cancer cells, and another group has also reported the chemomodulation effect by utilizing an ex vivo culture platform (73). Accompanied with the reduction of immunosuppressive cells, significant elevation of cytotoxic cells, such as XCL1⁺ CD8⁺ T cells, was noted. Therefore, our data showed that NAC powerfully restored intratumor immune balance in chemotherapy-responsive patients, decoding additional mechanisms that favored neoadjuvant chemotherapy in those patients.

Metastatic multifocal tumors harbored diverse immune microenvironmental composition. Even inside the same tissue of the same patient, the suppressive immune cells were characterized by differential abundance and distinct functional states. As for multifocal hepatocellular carcinoma, we recently reported that distinct lesions induced heterogeneous T-cell receptor expansions, which may impose immune selection (63). Another independent group also found that, in a patient with high-grade serous ovarian cancer with multifocal lesions (62), the progressive metastases excluded immune cells whereas the regressive lesions were infiltrated by CD8 T cells. This could partly explain the diverse drug responses of multifocal tumors of the same patients in the clinical practice. Our data will allow novel design of optimal rational combination therapies for multifocal CRLM such as combining NAC, metabolism checkpoint inhibitors, and immunotherapy.

Our study did have some limitations. Due to the limited patient number fitting the inclusion criteria (chemotherapy regimen and sample resectability), only four patients were recruited in our ST cohort. Second, the ST technique has not reached the single-cell resolution. Current 10X Visium platform merely supported an estimate of 1 to 10 cells in each spot (74). A more precise spatial map, such as Slide-seq and DBiT-seq (75–77), will be essential for better understanding how to therapeutically target the spatiotemporal heterogeneity of this complex disease. Third, it still remains unclear how the chemotherapy induces the functional changes of macrophages in metastatic tumors. Further experimental validation is required to validate that such state shift of macrophages is due to altered differentiation or population change (i.e., macrophage depletion or monocyte recruitment).

In summary, we present the spatial and cellular atlas of CRLM and pinpoint the functional impacts of neoadjuvant chemotherapy on TME. Our scMetabolism pipeline not only deepens our understanding into the immunobiology of CRLM, but also enables the unbiased discovery of active metabolic states in other complex tissues. These findings provide unique insights into the biology of cancer spreading and raise the possibility to target metabolism pathways in metastasis.

We collected samples of pathologically diagnosed CRLM from Zhongshan Hospital, Fudan University, with written informed consent and approval from the Institutional Review Board–approved protocols (B2020-348R). The study was conducted in accordance with ethical guidelines (Declaration of Helsinki). In detail, whether the tumors are resectable has all passed the MDT discussion composed by senior clinicians. In detail, 89 samples from 20 patients (n = 11, untreated; n = 5, PD/SD; n = 4, PR) underwent scRNA-seq, and eight samples from four patients (n = 2, untreated; n = 2, PR) were sequenced by ST. Another 104 samples from 27 patients were analyzed by mIHC (n = 17, untreated colorectal cancer and adjacent colon; n = 18, untreated liver metastasis and adjacent liver; n = 9, treated PR colorectal cancer and adjacent colon; n = 8, treated PR liver metastasis and adjacent liver).

Six-week-old female BALB/c mice were ordered from Shanghai Model Organisms Center maintained under specific pathogen-free housing with a maximum of five mice per cage. The experiments were performed following the institutional guidelines strictly and were approved by the Institutional Animal Care and Use Committee of the Shanghai Model Organisms Center (2019-0011). CT26 cells (5 × 10⁵) were intrasplenically injected at day 0, and a splenectomy was performed 30 minutes later. We then subcutaneously injected CT26 cells (3 × 10⁶) at day 0, as previously described (6). Samples of primary tumors (subcutaneous) and metastatic tumors were collected at day 20. Four mice underwent liver metastasis and one mouse underwent colorectal peritoneal metastases. Those primary tumors (n = 4) and metastatic tumors (n = 5) were then sequenced by scRNA-seq.

We cut the fresh tissue into small pieces and placed them in RPMI 1640 medium (Gibco; 11875093) for 1 hour at 37°C. In detail, 1 mg/mL collagenase IV (Gibco; 17104019) and 0.4 mg/mL hyaluronidase were used for digestion. To dissociate the tumor tissues, we used the gentleMACS tumor dissociation kit (Miltenyi Biotec; 130095929). We then suspended and filtered cells through a 400-μm mesh sieve. According to the manufacturer's instructions, flow cytometry was performed on a BD LSRFortessa cell analyzer (BD Biosciences) and analyzed using FlowJo software (V9.3.2), as we described previously (25).

Following tissue digestion, samples were sequenced by 10x Chromium single-cell platform (10x Genomics, 30 v3 chemistry, Rev C). In detail, Gel Bead-In EMulsion (GEM) generation and barcoding was first performed. Then, post GEM-RT cleanup and cDNA amplification were conducted. Finally, the 3′ gene-expression library construction was performed following the manufacturer's instructions. Completed libraries were sequenced on NovaSeq (Illumina). The CellRanger (V3) pipeline was used for demultiplexing, barcode processing, alignment, and initial clustering of the raw scRNA-seq profiles (25).

After resection, samples were collected, and a tissue optimization experiment was performed (10x Genomics, Visium Spatial Tissue Optimization, Rev A). In detail, tissue staining and imaging were first performed. Then, tissue permeabilization and fluorescent cDNA synthesis were performed, and RNA quality was checked (RIN > 7.0). Finally, tissue was removed, and slide imaging was performed. Thereafter, Visium Spatial Gene Expression was conducted (10x Genomics, Visium Spatial Gene-Expression Reagent Kits, Rev B). First, tissues were fixed and stained with imaging. Second, tissue permeabilization with different times (3, 6, 12, 18, 24, and 30 minutes) was performed, and tissue permeabilization with 12 minutes results in maximum fluorescence signal in both tumor regions and adjacent normal regions in the liver and colon. Next, reverse transcription, second-strand synthesis and denaturation, cDNA amplification and QC, Visium spatial gene-expression library construction, and ST sequencing were performed, respectively, following the manufacturer's instructions. This generated data sets containing 3,826 (ST-P1, liver), 4,658 (ST-P2, liver), 3,695 (ST-P3, liver), 3,721 (ST-P4, liver), 3,313 (ST-P1, colon), 4,174 (ST-P2, colon), 4,007 (ST-P3, colon), and 3,902 (ST-P4, colon) spots.

We performed the fluorescent dyes by using the Osteopontin/SPP1 rabbit anti-human antibody (Abcam; catalog no. ab214050), Ki-67 rabbit anti-human antibody (IBP; catalog no. IR098), CD68 mouse anti-human antibody (Maxim; catalog no. Kit-0026), CCL18 rabbit anti-human antibody (Sino Biological; catalog no. 10502-T16, RRID:AB_2860192), GranzymeK mouse anti-human antibody (ProteinTech; catalog no. 67272-1-Ig, RRID:AB_2882541), GZMB rabbit anti-human antibody (Biolynx; catalog no. BX50024), CD56 rabbit anti-human antibody (IBP; catalog no. IR040), CTLA4 mouse anti-human antibody (Origene; catalog no. TA810299), CD45 mouse anti-human antibody (BioLegend; catalog no. 304002, RRID:AB_2661811), CD15 mouse anti-human antibody (Maxim; catalog no. MAB-0779), DAPI (BioLegend; catalog no. 422801), MIF polyclonal rabbit anti-human antibody (Abcam, ab65869), GSTO1 rabbit anti-human antibody (Abcam, ab129106), IL4I1 rabbit anti-human antibody (Abcam, ab222102), SIRPA mouse anti-human antibody (OriGene, TA807751), CD47 rabbit anti-human antibody (Abcam, ab226837), and cytokeratin Pan (Pan Ck) mouse anti-human antibody (LBP Medicine, IHC-M067). Following the manufacturer's instructions (PerkinElmer, Opal Kit), we scanned the slides by using the PerkinElmer Vectra3 platform and quantified the results by using PerkinElmer Vectra3 platform as previously described (25, 78).

Raw fastq files were first processed by CellRanger V3. Raw sequencing reads were then aligned to human genome (GRCh38, ENSEMBL). We then used Seurat (V3.2.2; refs. 22, 37) to process the unique molecular identifier (UMI) count matrix and integrates all cells according to sample ID by using Harmony (79). We performed the clustering analysis of all cells and only reserved the CD45⁺ clusters. To define the main immune cell types, we then split the cells according to commonly used cell markers (CD4 for CD4 T cells, FOXP3 for Tregs, CD8A for CD8 T cells, SLC4A10 for MAIT cells, CD79A for B cells, MZB1 for plasma cells, LYZ for myeloid cells, and KLRC1 for NK cells; refs. 25, 80).

After labeling the main cell types, we performed clustering analysis of the split main cell populations. As for the CD4⁺ cells, we found eight clusters. Cluster 1 was SELL⁺ and was hence defined as CD4⁺-naïve T cells (81). CCL4, CCL5, and IFNG were highly expressed in cluster 2 (CCL5⁺ CCL4⁺ IFNG⁺ CD4⁺ T cells). Cluster 3 expressed CCL5 but not CCL4, so this cluster was classified as CCL5⁺ CCL4⁻ CD4⁺ T cells. Given clusters 2 and 3 were similar, we directly combined them into a same cluster and labeled it as CCL5⁺ CD4⁺ T cells. Cluster 4 expressed TOX2, TIGIT, and PDCD1, and we defined it as PDCD1⁺ CD4⁺ T cells. Cluster 5 expressed KLRB1 and CCR6, and we defined it as CCR6⁺ Th17 cells (82). Cluster 6 expressed LTB, and we defined it as LTB⁺ CD4⁺ T cells. Cluster 7 expressed HSP genes such as HSPA1A and HSP90AA1, and we defined it as HSP⁺ CD4⁺ T cells. Cluster 8 expressed FOXP3, TIGIT, and CTLA4, and we defined it as FOXP3⁺ Tregs.

As for the CD8⁺ cells, we found 10 clusters. Cluster 1 was GZMK⁺; hence, it was defined as GZMK⁺ CD8⁺ T cells. Cluster 2 expressed CD52; we hence defined it as CD52⁺ CD8⁺ T cells. Cluster 3 expressed SLC4A10, DPP4, and KLRB1 (83), and this cluster was defined as MAIT cells. Cluster 4 expressed CTLA4, TNFRSF9, HAVCR2, and TIGIT, so this cluster was classified as CTLA4⁺ CD8⁺ T cells (exhausted T cells; ref. 84). Cluster 5 expressed GIMAP family genes such as GIMAP7, GIMAP4, and GIMAP1. So we defined it as GIMAP7⁺ CD8⁺ T cells. Cluster 6 expressed SELL and CCR7, and we defined it as naïve CD8⁺ T cells (81). Cluster 7 expressed FGFBP2 and GZMB, and we defined it as FGFBP2⁺ GZMB⁺ CD8⁺ T cells (effector T cells; ref. 85). Cluster 9 expressed XCL1, and we defined it as XCL1⁺ CD8⁺ T cells. Cluster 10 expressed HSP genes such as HSPA6, HSPA1A, HSPD1, and HSPH1, and we defined it as HSP⁺ CD8⁺ T cells. In addition, clustering analysis also generates a cluster with high expression of TRDC but low expression of the CD8A and CD4 genes. We hence define it as TRDC⁺ γδ T cells (86).

As for the NK cells, we found four clusters. Cluster 1 expressed GZMK and XCL1 and was hence named as GZMK⁺ NK cells. Cluster 2 expressed GZMK, IL7R, and GPR183, and was hence defined as GZMK⁺ IL7R⁺ GPR183⁺ NK cells. Cluster 3 expressed GZMB and GNLY and was hence named as GZMB⁺ NK cells. Cluster 4 expressed GZMB and MYOM2. Thus, it was defined as GZMB⁺ MYOM2⁺ NK cells. Given that clusters 1 and 2 both expressed GZMK, we hence combined them into the same cell cluster and labeled it as GZMK⁺ NK cells. Given that clusters 3 and 4 both expressed GZMB, we combined them into the same cell cluster and labeled it as GZMB⁺ NK cells.

As for the myeloid cells, there existed nine clusters. Cluster 1 was MKI67⁺ and was hence defined as MKI67⁺ proliferation macrophages. S100A8, S100A9, and S100A12 were highly expressed in cluster 2 (S100A8⁺ monocytes; ref. 8). Cluster 3 expressed CLEC10A, and this cluster was defined as cDC2 cells (87). Cluster 4 expressed FCGR3B, so this cluster was classified as neutrophils (88). Cluster 5 expressed HSP genes such as HSPA1B, HSPH1, HSPA1A, and HSPB1, and we defined it as HSP⁺ monocytes. Cluster 6 expressed MRC1 and CCL18, and we defined it as MRC1⁺ CCL18⁺ M2-like macrophages (89). Cluster 7 expressed CLEC9A, and we defined it as CLEC9A⁺ cDC1 (8). Cluster 8 expressed SPP1 and we defined it as SPP1⁺ macrophages (8). Cluster 9 expressed CXCL10, and we defined it as CXCL10⁺ M1-like macrophages (89).

As for the B cells, cluster 1 was MKI67⁺ and was hence defined as MKI67⁺ proliferation B cells. Cluster 2 expressed AIM2, and this cluster was defined as AIM2⁺ memory B cells. Cluster 3 expressed TCL1A and SELL, so this cluster was classified as TCL1A⁺-naïve B cells (90). Cluster 4 expressed EGR1, which is essential for B-cell activation (91), and was labeled as EGR1⁺-activated B cells. Cluster 5 expressed HSP genes such as HSPA1A, HSP90AA1, HSPE1, and HSPH1, and we defined it as HSPA1A⁺ B cells. Cluster 6 expressed a broad range of IGHA genes such as IGHA2 and IGHA1. We hence labeled it as IGHA⁺ plasma B cells. Cluster 6 showed high expression of IGHG1, IGHG2, IGHG4, and IGHG3. This indicates that such cluster is IGHG⁺ plasma B cells.

As for quality control of scRNA-seq data, we set the cutoff min.cells as 3 and min.features as 500 in the CreateSeuratObject function in Seurat (V3.2.2). To remove the potential doublets, we perform two rounds of doublet removal by using doubletFinder (V2.0.3; ref. 92). We first run doubletFinder on all cells and then perform similar analysis on subcellular populations. When performing clustering analysis, we remove the cluster with high expression of mitochondrial genes (93).

To visualize our data, we use Seurat (V3.2.2) (22), ggplot2 (V3.3.2), dittoSeq (V1.1.10) (94), ggrepel (V0.8.2), ggpubr (V0.4.0), ggsignif (V0.6.0), pheatmap (V1.0.12), circlize (V0.4.10), and cowplot (V1.1.0).

To perform the differential cell-type composition analysis, we utilized the scDC algorithm (23) to conduct the bias-corrected and accelerated bootstrap analysis for cell-type proportion comparison. We used the scDC_noClustering function to perform the bootstrap resampling cell-type composition analysis (calCI = TRUE, calCI_method = BCa, nboot = 1,000). To find out the tissue-specific cellular types, we used the median value of bootstrap resampled cellular composition, filtered the immune cells with low composition, and performed the unsupervised clustering analysis.

We used Space Ranger V1.2.1 to process raw fastq files. We next used Seurat (V3.2.2; ref. 22) to process the Space Ranger output files. In detail, we used SCTransform (95) to normalize data, RunPCA to perform dimension reduction, FindNeighbors and FindClusters to cluster the ST spots, and RunUMAP to visualize the data. Given that each spot contains several cells (74), we next combined the scRNA-seq data and predicted main cell types by using Seurat (95). As for subpopulation annotation, we found that Seurat did not generate uniform results, which might be due to the overlapping cell markers among subcellular populations. We hence used ssGSEA (96) to score the subcell types, which has been proved to be robust for ST data (29).

We fetched bulk RNA-seq data of patients with colorectal cancer from the TCGA (https://portal.gdc.cancer.gov/). We next used ssGSEA (96) to estimate the immune cell infiltration level by using specific cell markers (avg_logFC > 0.25, pct.1 > 0.5, and pct.2 < 0.5). We next used survival (V3.2-7) and survminer (V0.4.8) to generate the Kaplan–Meier plot, and the P value was calculated through the log-rank test. To validate the cell infiltration level, we used another CRLM microarray data set (ref. 30; https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE41568) and utilized ssGSEA to estimate the immune cell level.

To infer the differentiation trajectory of intratumor myeloid cells, we first used monocle (V2.14.0; ref. 97) to infer the pseudotime of each cell (method = “DDRTress,” ordering_genes = marker genes). To compare the myeloid cell lineage/clustering difference between untreated samples and NAC-treated samples, we used dyno (V0.1.2, method = “slingshot”; refs. 57, 98) to analyze and visualize.

We used CellPhoneDB V2.0 to infer the ligand–receptor cross-talk between single cells (99). As for the differential tumor–immune cell cross-talk analysis (colorectal cancer vs. liver metastasis), we separately computed the tumor–immune cell cross-talk in colorectal cancer and liver metastasis. We selected the significantly differential cross-talk ligand–receptor pairs and used circlize (V0.4.10) to visualize the results.

We used FindMarkers function to perform differential gene-expression analysis of different cellular populations and set P_adj = 0.05 as the cut-off value. We next performed the pathway enrichment analysis of those differentially expressed genes by using clusterProfiler V3.19.0 (refs. 43, 100; enricher function, KEGG gene sets, or cancer hallmark gene sets from msigdb). We set P = 0.05 as the cut-off value.

To compare the cellular clustering difference between untreated samples and NAC-treated samples, we separately performed clustering analysis on colorectal cancer and liver metastasis samples by using TooManyCells (55). In detail, we used too-many-cells make-tree function and set draw-collection as “PieRing.”

scMetabolism was designed to easily quantify single-cell metabolic activity by using single-line command. Before metabolism gene signature scoring, users can choose whether to impute their scRNA-seq matrix (38). The core function of scMetabolism is to quantify metabolic pathway gene sets. We combined the published gene sets (36) and manually reviewed gene sets from the KEGG database (42, 43) and the REACTOME database (43, 44) to generate the list of metabolic gene sets. As for the quantification method, scMetabolism currently supports VISION (39), AUCell (40), and ssGSEA (41). Users can directly use the scMetabolism pipeline based on their Seurat object (22, 37).

To evaluate the performance of scMetabolism, we used four different data sets, including PBMCs, melanoma, head and neck squamous cell carcinoma, and testis (45–47). To test the accuracy of this pipeline, we compared VISION and AUCell with ssGSEA, which was previously reported to be appropriate for single-cell metabolism quantification (36). We found that VISION shows good correlation with ssGSEA, with the median Spearman ρ estimate over 0.8 in all data sets (Supplementary Fig. S4A). As for AUCell, its correlation with ssGSEA ranges from 0.5 to 0.8 (Supplementary Fig. S4A). Hence, ssGSEA and VISION show best consistency among all methods.

Recent pipeline from Locasale Lab by default imputes the missing gene-expression values. However, it still remains unclear how dropout affects the metabolic quantification and whether imputation can improve the performance. We hence simulated the dropout events by randomly select 5%, 10%, 20%, and 30% expressed genes and set their values to 0. In those four simulated dropout gene-expression matrix, we used ALRA (38) to impute the simulated missing data. Next, we used VISION, AUCell, and ssGSEA to quantify the metabolic activity and perform the correlation analysis with the raw data matrix. Interestingly, we found imputation may not improve the metabolic quantification performance (Supplementary Fig. S4B–S4E). In the 5% dropout-imputation data set, we observed a sharp decreased correlation value, whereas the correlation remains similar in the 30% dropout-imputation data set in all data sets. Among all three methods, results of VISION and ssGSEA are less affected by dropout. Among all methods, VISION is the fastest way for metabolism quantification. We hence set VISION as the default method and set “not imputation” as the default option. Finally, we developed a web server, scMetabolism online (http://www.cancerdiversity.asia/scMetabolism/), for experimental biologist to easily perform metabolism quantification.

We used shiny (V1.15.0) to construct the web server and runs on a Linux-based Apache Web server (http://www.apache.org) as previously described (101, 102). We tested this web server on different web browsers, including Chroma, Safari, and Internet Explorer. Users are not required to register or login to access features in the database.

This raw sequencing data generated during this study are available at The National Omics Data Encyclopedia (https://www.biosino.org/node/project/detail/OEP001756). The analyzed data are publicly available in the web portal (http://www.cancerdiversity.asia/scCRLM). The R package of scMetabolism is available on Github (https://www.github.com/wu-yc/scMetabolism) and can be directly installed. The online version of this package, named as scMetabolism online, is also publicly available (http://www.cancerdiversity.asia/scMetabolism).

No disclosures were reported.

Y. Wu: Software, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. S. Yang: Software, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. J. Ma: Software, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. Z. Chen: Software, formal analysis, validation, investigation, visualization, methodology, writing–original draft, writing–review and editing. G. Song: Investigation, writing–review and editing. D. Rao: Investigation, writing–review and editing. Y. Cheng: Investigation, writing–review and editing. S. Huang: Investigation, writing–review and editing. Y. Liu: Investigation, writing–review and editing. S. Jiang: Investigation, writing–review and editing. J. Liu: Investigation, writing–review and editing. X. Huang: Investigation, methodology, writing–review and editing. X. Wang: Investigation, methodology, writing–review and editing. S. Qiu: Investigation. J. Xu: Investigation, methodology, writing–review and editing. R. Xi: Investigation, methodology, writing–review and editing. F. Bai: Investigation, methodology, writing–review and editing. J. Zhou: Investigation, methodology, writing–review and editing. J. Fan: Conceptualization, investigation, methodology, writing–original draft, writing–review and editing. X. Zhang: Conceptualization, resources, supervision, funding acquisition, investigation, methodology, writing–original draft, writing–review and editing. Q. Gao: Conceptualization, resources, supervision, funding acquisition, investigation, visualization, methodology, writing–original draft, writing–review and editing.

This work was supported by the National Natural Science Foundation of China (no. 81961128025; to Q. Gao); National Natural Science Foundation of China (no. 91942313; to X. Huang); Program of Shanghai Academic Research Leader (no. 19XD1420700; to Q. Gao); Shanghai Municipal Key Clinical Specialty; The Strategic Priority Research Program (XDPB0303) and Frontier Science Key Research Project (QYZDB-SSW-SMC036) of Chinese Academy of Sciences (X. Zhang); and Shanghai Municipal Science and Technology Major Project (no. 2019SHZDZX02; to X. Zhang). We thank Prof. Jun Qin (Chinese Academy of Sciences), Dr. Lu Meng (Chinese Academy of Sciences), Prof. Shu Zhang (Fudan University), Dr. Liangqing Dong (Fudan University), Prof. Yihui Fan (Nantong University), and Dr. Jianming Zeng (University of Macau) for their support in this work.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.