A Single-Cell Atlas of Liver Metastases of Colorectal Cancer Reveals the Reprogramming of the Tumor Microenvironment in Response to Preoperative Chemotherapy

Metastasis is the primary cause of cancer-related mortality in colorectal cancer (CRC) patients. How to improve therapeutic options for patients with metastatic CRC (mCRC) is the core question for CRC treatment. However, the complexity and diversity of stromal context of the tumor microenvironment (TME) in liver metastases of CRC is not fully understood, and its inuence on response to chemotherapy is unclear. Here we provide an in-depth analysis of transcriptional landscape of primary CRC, matched liver metastases and blood at single-cell resolution, and perform a systematic examination of transcriptional changes and phenotypic alteration of the TME in response to preoperative chemotherapy (PC). Based on 111,292 single-cell transcriptomes, our study reveal that TME of treatment-naïve tumors is characterized by higher abundance of less-activated B cells and higher heterogeneity of tumor-associated macrophages (TAMs). By contrast, in tumors treated with preoperative chemotherapy, we nd reprogramming of B cell activation, lower diversity of TAMs with immature and less activated phenotype, lower abundance of both dysfunctional T cells and ECM-remodeling cancer-associated broblasts (CAFs), an accumulation of myobroblast. Our study may provide a foundation to investigate the cellular mechanisms underlying liver metastasis of CRC and response to preoperative chemotherapy, and open up new possibilities for predicting and guiding therapeutic responsiveness. performed on BD Inux (BD Biosciences) using 488nm (calcein AM, 530/40 lter), 638nm (TO-PRO-3, 670/30 lter) lasers. Singlets were captured and doublets were discard through forward scatter height and width parameters. Viable cells were recognized as Calcein AM (high) and TO-PRO-3 (low) cell cluster. For PBMC sample, only lymphocytes and monocytes clusters were gated for further sorting. Viable single cells were resuspended in HBSS with 0.04% BSA. Viability was conrmed to be >90% using trypan blue (ThermoFisher Scientic) exclusion prior to scRNA-seq process.


Introduction
Metastasis is the primary cause of cancer-related mortality in colorectal cancer (CRC) patients (Jemal et al., 2011;Vatandoust, 2015). Preoperative chemotherapy (PC) aims at reducing tumor load, which may reduce the risk of local relapse and promoting patients with initially unresectable mCRC to resectable liver metastases (Glimelius et al., 2013;Kapiteijn et al., 2001). Nevertheless, despite theoretical bene t and randomized trail demonstrations (Nordlinger et al., 2008), whether the patients undergoing chemotherapy and resection have long-term bene t is still questionable. How to provide optimal treatment of CRC patients with liver metastasis remains a pivotal issue.
Understanding the complex cellular and phenotypic diversity within the tumor microenvironment (TME) may pave the way for the development of effective treatment in cancer especially in metastatic disease.
1. single-cell analysis of tumor microenvironment of liver metastases of colorectal cancer.
To gain a better understanding of tumor microenvironment, and investigate how TME response to preoperative chemotherapy in liver metastases of CRC, we performed scRNA-seq of 15 samples from three sites (primary CRC, matched liver metastases and blood) of six CRC patients with liver metastatic disease (Table S1). While patients COL15, COL17 and COL18 have been treated with preoperative chemotherapy, the others are treatment naïve. All patients were classi ed as microsatellite-stable (MSS) with invasive adenocarcinomas and late-stage (IV) disease (Methods , Table S1).
Viable single cells were sorted and used for droplet-based single-cell RNA sequencing (scRNA-seq). After quality control (see Methods), we obtained transcriptome data for 111,292 single cells from primary CRC (n = 6), matched liver metastases (n = 6), and PBMCs (n = 3) for liver metastases of CRC (Table S2). Then we clustered single cells using shared nearest neighbor (SNN) clustering based on signi cant principal components (PCs), and visualized cell clusters using t-distributed stochastic neighbor embedding (t-SNE) (Fig. 1a, Methods). The major cell populations (including T cells /natural killer (NK) cells, B/plasma cells, cancer-associated broblasts (CAFs), endothelial and myeloid cells) were annotated with canonical marker genes (Fig. 1b), treatment state and tissue origin are mapped in Fig. 1c (also see Figure S1a), selected representative markers of each cell type is presented in Fig. 1d and Fig. 1e and Figure S1b.
To characterize the ecosystems of primary and metastatic tumor in CRC, and better understand how TME response to preoperative chemotherapy in liver metastases of CRC, we focus on major cell types of TME (T/ NK cells, B cells, Myeloid cells and CAFs), for each compartment, we re-centred, scaled, normalized and re-clustered the data. Ultimately, we obtained 28 myeloid clusters (6 dendritic cells, 18 tumorassociated macrophages (TAMs), 1 monocytes, 2 myeloid-derived suppressor-like cells (MDSCs-like) and 1 mast cells), 17 B cell clusters, 10 mesenchymal cell clusters (1 endothelial, 6 cancer-associated broblasts and 3 myo broblasts) and T/NK cell clusters. Each cluster was composed of cells from different patients, and for each cluster, the distribution of tissue and treatment status were different ( Figure S1c).

Preoperative chemotherapy promotes the activation of B cells in the primary CRC
First, we investigate how preoperative chemotherapy in uenced the ability of B cells. Interestingly, B cells were prevalent in primary CRC (account for 0.1% of all stromal cells in primary CRC), but depleted signi cantly in liver metastases (account for 0.01% of all stromal cells in liver metastases) (Fig. 2a). Subclustering of B cells revealed 17 subpopulations (Fig. 2b). Among them, 15 clusters were mature B cells (with eight subsets from tumor lesion and ve clusters from peripheral blood, Fig. 2b), which are characterized by highly expression of CD20 (MS4A1). While cluster 3 represents plasma cells which is characterized by highly abundant immunoglobulin (IGHG1, IGHG2, IGHG3, IGHG4 and IGHA2), cluster 16 represent plasmablast, characterized by upregulation of immunoglobulin and cell proliferation markers (e.g. MKI67, CDC20, CDKN3 and CCNB2) (Table S3).
Since few B cells in ltrated in the liver metastases, we mainly focused on B cells in the primary CRC. In primary tumor, B cells were separated into different compartments, the B cell populations from tumor treated with preoperative chemotherapy were quite distinct from that in treatment-naïve tumors ( Fig. 2c and 2d). We found that cluster 0, 2, 6 and 11 were enriched speci c in treatment-naïve tumors, whereas cluster 1, 4, 9, 10 and 12 were almost exclusively present in treated tumors (Fig. 2c, 2d).
On closer examination of the two groups, we nd that they had distinct phenotypes (Table S3). Untreated tumors-derived B cells exhibited a naïve and in ammatory phenotype, with cluster 6 expressing IgD (IGHD), cluster 0 expressing immature B cell marker VPREB3 and cluster 2 expressing in ammatory transcription factor NF-kB (NFKBIA), lipid molecules (e.g. APOE and APOC1) and cytokines (e.g. AREG).Whereas B cells derived from treated tumors (cluster 1, 4, 9, 8 and 12) have an activated immune activation phenotype with upregulation of immunocostimulatory molecules (e.g.,CD82, CD83 and CLECL) and MHCs molecules (e.g.,HLA-DRA, HLA-DRB5) (Table S3). GO analysis also con rms that upregulated genes associated with treated tumors were enriched in antigen processing and presentation ( Figure S2a). , higher expression of immunoglobulin, such as IGLC3,  JCHAIN, IGHG1, IGHG3, IGHG4 and IGHA1, were observed in treated tumors (Fig. 2e, 2f, and Figure S2b), implying class switch recombination (CSR, also known as isotype switching) may occur after preoperative chemotherapy. Moreover, the expression of MHCs molecules (such as HLA-A, HLA-DQB1, HLA-DQA2 and HLA-DRB1) were also elevated in treated tumors (Fig. 2g). In line with these results, GO analysis shows that B cells in tumors treated with PC enriched in immune activation signaling pathway, including immune response-activating cell surface receptor, antigen processing and presentation signal pathway, and T cell co-stimulation (Fig. 2h) Figure S3a. While TAMs and DCs were enriched in both primary tumor and liver metastases, monocytes and MDSCs-like cells were prevalent in blood and mast cells were mainly enriched in primary CRC ( Fig. 3a and 3b). , S100A family genes are highly expressed, including S100A12, S100A9, S100A8, together with FCN1 and VCAN, but MHC I and MHC II molecules tend to be lowly expressed (Table S3). The rst subset, including cluster 0, 12 and 13, were classi ed as MHC low TAMs. MHC low TAMs exhibited weak capacity of antigen presentation and immune activation, with lowly expressed MHC I and MHC II genes (e.g. HLA-A, HLA-C, HLA-DRA and HLA-DRB1), and lowly expressed immune-costimulatory genes (e.g. CD80 and CD86). Cluster 12, classi ed as IL1B + MHC low TAMs, characterized by upregulation of a large number of in ammatory and chemokine genes (e.g., IL1B, IL6, S100A9, S100A8, CXCL8, CXCL3 and CXCL1) which involved in recruiting and regulating of immune cells. Cluster 12 were presented in the primary CRC, in patients who are treatment-naïve (Fig. 3c, right panel, Fig. 3d). Cluster 0 (THBS1 + MHC low TAMs) were prevalent in liver (Fig. 3c, right panel), with highly expression of THBS1, MARCO and genes promoting proliferation of epithelia cells and angiogenesis, such as EREG, AREG and VEGFA, which could stimulate tumor growth and progression. Cluster 13 (MARCO + MHC low TAMs), mainly derived from liver resident Kupffer cells (KCs; the tissue resident macrophage of liver), may exert a tolerogenic or immune inhibitory functions in the liver metastases through MARCO, VSIG4 and CD163.

Comparing with B cells in untreated tumors
The second subset, including cluster 6, 9, 14, 18 and 19, were classi ed as lipid-associated TAMs (LAMs), which were characterized by upregulation of genes involved in lipid metabolism (e.g., APOC1, APOE and LPL), extracellular matrix (ECM) degradation (e.g., MMP7, MMP9 and SPRAC) and complement activation (e.g., C1QA, C1QB and C2). They also highly expressed TREM2 (encoding lipid receptor) and LGALS3 (associated with immune suppression), which were recently found to be associated with metabolic diseases (Jaitin et al., 2019). Among the subset, cluster 6 and cluster 9 showed higher expression levels of matrix metalloproteinases (MMPs) compared to other clusters. GO analysis of upregulated genes in each cluster demonstrated that TAMs were enriched in neutrophil activation and degranulation ( Figure S3c). In the primary CRC, lipid-associated TAMs (LAMs) is mainly presented in treatment-naïve tumors, but is shared in treated and untreated in the liver metastases (Fig. 3c, right panel, Fig. 3d).
The third subset, including cluster 4 and 11, were enriched for genes involved in regulation of myeloid leukocyte differentiation, were prevalent in tumors treated with PC (Fig. 3c, right panel, Fig. 3d). With highly expression of FOS, cluster 4, mainly from liver metastases, may be monocyte-derived, which could differentiate into macrophage (Azizi et al., 2018). In addition to the canonical myeloid marker gene CD33, cluster 11 was also associated with highly expression of CD4, which is typically expressed by monocytes and involved in triggering cytokine expression and the differentiation of monocytes into functional . Moreover, MHC genes tend to be expressed at a higher level in cluster 11 cells than cluster 4, suggesting cluster 11 cells are more mature and activated than the latter. Thus, we classify this subset as monocyte-derived immature TAMs.
The fourth subset, including cluster 20 and 23, are immune regulatory TAMs characterized by upregulation of immune suppressive genes, such as CD274, CCL2, IL10 and TGFB2, and prevalent in primary CRC in treatment-naïve tumors (Fig. 3c, right panel). Strkingly, both of the two clusters are also associated with highly expression of MHC and co-stimulating genes, a signature of immune activation and anti-tumor activities. Together, we uncovered "double-agent" immune regulatory TAMs with coexpression of both immune activated and immune suppressive genes, which indicate complex interactions between anti-tumor and pro-tumor activities.
In addition to the four clusters described above, there are three more clusters, named as M01, M03 and M26, that can not be classi ed into the above-mentioned ve subsets. Genes highly expressed in these three clusters are associated with different bioligical processes and cellular functions. M26 (MKI67 + TAMs) is marked with highly expression of genes in cell proliferation (e.g. MKI67, Table S3). In addition to EGF and MACRO, heat-shock genes are also highly expressed in M01. M03 (CXCL10 + TAMs) highly expressed genes involved in response to interferon-gamma (GBP1, STAT1, IFITM3 and PARP14) and cellular response to zinc ion (MT2A, MT1X and MT1F).
Preoperative chemotherapy reprogramming TAMs from highly heterogeneity to immature and less activated phenotypes TAMs from treatment-naïve tumors and chemotherapy-treated tumors exhibited distinct phenotypes in both primary tumor and liver metastases. In the primary CRC, TAMs in untreated tumors showed higher heterogeneity, they present distinct phenotypes compared with TAMs in tumors treated with preoperative chemotherapy (Fig. 3d). Focusing on the treatment state, we observed some phenotypes shared between treated and untreated tumors, such as MKI67 + TAMs (M26). However, immature TAMs (M04, M11) and MHC low TAMs (M00, THBS1 + MHC low TAMs) were largely speci c to treated tumors ( TAMs from tumors treated with PC are distinct from those in treatment-naïve tumors ( Figure S3d) at the transcriptome level. In the niche of treatment-naïve tumors, TAMs from the primary CRC were enriched for processes of neutrophil activation, responding to IFN-γ and regulation of broblast proliferation ( Figure  S3d), indicating its proin ammatory phenotype . TAMs from the metastases enriched for antigen processing and presentation, neutrophil activation and responding to IFN-γ. In the ecosystem of treated tumors, TAMs were characterized by unregulated genes involved in protein targeting to endoplasmic reticulum (ER) and RNA catabolic progress in the primary tumors, while in the metastatic sites, TAMs were enriched for myeloid leukocyte chemotaxis and migration processes, which could be related to monocyte-like TAMs aggregated in this lesion.
In general, our results showed that preoperative chemotherapy suppressed the diversity of TAMs, and only the immature TAMs and THBS1 + MHC low TAMs in ltrated in the primary tumors, while more TAMs aggregated in liver metastases and tended to be immature and less activated. Thus, chemotherapy reprogramming TAMs from highly heterogeneity to immature and less activated phenotypes 4. Preoperative chemotherapy decreases the abundance of ECM-remodeling CAFs and dysfunctional T cells, but promotes the accumulation of myo broblast In our study, 1383 CAFs were detected and classi ed into 9 clusters. Notably, CAFs were signi cantly more abundant in primary CRC than liver metastatic tumor (Fig. 4a). Using gene expression pro le, we classi ed CAFs into three major subsets (Fig. 4b, c), including secretory CAFs (cluster 0, 1, 6 and 7), ECMremodeling CAFs (cluster 2 and 8) and contractile CAFs (cluster 3, 4 and 5) (Fig. 4b). The subset of secretory CAFs highly expressed secretory genes such as various growth factor (e.g. IGF1, PDGFD, FGF7 and VEGFB) which mediate angiogenesis and cancer cell proliferation, some signal molecules (e.g. BMP4 and WNT2B) that were able to maintain cancer stem cell niche, complements (e.g. C1S and C3) and chemokines (e.g. CCL2, CXCL12 and CXCL14) that regulated tumor immunity and in ammation. The ECM-remodeling CAFs highly expressed abundant extracellular matrix (ECM) proteins (such as ECM collagens and bronectin), strongly associated with a brotic matrix (Fig. 4c). They also expressed a large number of ECM proteases, which altered ECM and assisted angiogenesis and metastasis (Koliaraki et al., 2017). The contractile CAFs enriched for genes involved in regulation of cell contraction (Fig. 4c), suggesting some distinct phenotypes. Since CAFs have numerous potential cellular sources, we suggested the cluster 4 tended to be myo broblast nature for their upregulation of myo broblast markers (e.g. ACTA2 and TAGLN) and genes involved in myogenesis (e.g. MYH11, PLN and CNN1) (Elyada et al., 2019; Koliaraki et al., 2017). Cluster 5, highly expressing pericyte-associated markers (e.g. RGS5 and CSPG4), were largely pericytes source (Koliaraki et al., 2017). Cluster 3 upregulate genes involved in stress response (e.g., JUN, BAG3 and HSPA2) and only present in liver metastases, suggesting they may be triggered as adaptation to tumor microenvironments.
CAFs derived from treated and untreated tumor exhibited distinct phenotypes (Fig. 4b, right panel). ECMremodeling CAFs were prevalent in the primary CRC in treatment-naïve tumors (Fig. 4d), whereas contractile CAFs were more prevalent in tumors treated with preoperative chemotherapy both in primary (cluster 4) and liver metastases (cluster 3 and cluster 5). Secretory CAFs were observed in the primary CRC, and were mainly enriched in treated tumors (except cluster 0) (Fig. 4d). Comparison the CAFs from treatment-naive tumors and treated tumors, CAFs from treatment-naive tumors were strongly enriched for processes of ECM organization and collagen metabolism (Fig. 4e), whereas CAFs in tumors treated with PC were signi cantly enriched pathways involved in regulating muscle cell differentiation, immune system (T cell activation) and epithelial cell proliferation (Fig. 4e). As we know, ECM remodeling is one important feature of CAFs common to progressive tumor and promotes metastasis (Bonnans et al., 2014). The observations here indicated that preoperative chemotherapy suppressed ECM remodeling in CAFs, but promotes accumulation of myo broblast and diverse of secretory CAFs in metastases of CRC.

Preoperative chemotherapy suppress the CD8 + dysfunctional T cells.
Here, we also identi ed different phenotypes of T cells (Fig. 5a) In addition, we also identi ed two MKI67 + CD8 + T cell populations (T36 and T24) in the T cell populations (Fig. 5b). Closer examination of those two clusters revealed cluster 36 were prevalent in the liver metastases (Fig. 5b, 96% in liver and 4% in the primary CRC), whereas cluster 24 were enriched in both the primary CRC and the liver (53% in CRC and 45% in liver), indicating that primary CRC and liver metastases share cluster 24, and the liver metastases has their speci c T cells characterized by high proliferation property (T36). To identify differences of those two subpopulations of T cells, we detected differentially expressed genes (DEGs) between T24 and T36 with p-value less than 1%, log2-fold change more than l and more than 10% cells expressing the gene. The results revealed heat-shock protein genes such as HSP90AA1, HSPA6, HSPA1A, HSPA1B and DNAJA4, and some molecular chaperones (e.g. BAG3 and HSPB1) were greatly upregulated in cluster 36 (Fig. 5b). Gene ontology (GO) analysis showed that cluster 36 was mainly enriched for protein folding, and response to stress (response to temperature stimulus, response to heat, response to topologically incorrect protein), suggesting they may be activated as adaptation to tumor microenvironments.
Many studies revealed that chemotherapy could impact the compartment of T cells. For example, chemotherapy could increase TIL (tumor-in ltrating lymphocyte) in ltration and decrease Tregs accumulation and proliferation in CRC patients (Maeda et al., 2011;Terme et al., 2013). Comparison of the cellular diversity of T cells in treatment-naïve tumors and tumors treated with PC, we found most subsets of CD4 + T cell were shared in chemotherapy-treated and untreated tumors (Fig. 5c, 5d), but the phenotypes of CD8 + T cells were signi cantly different (Fig. 5c, 5d).
In the treatment-naïve tumors, different phenotypes of CD8 + T cell were in ltrated in the primary tumor, including effector T cells and exhausted T cells, but in the metastatic sites, for the CD8 + T cells, only dysfunctional or exhausted T cells were accumulated in the liver (Fig. 5c). Most signi cantly, preoperative chemotherapy suppressed accumulation of dysfunctional T cell signi cantly both in the niches of primary CRC and liver metastases (Fig. 5e), which were validated by immuno uorescence ( Figure S5). Consistent with previous studies (Maeda et al., 2011;Terme et al., 2013), chemotherapy decrease the accumulation of Tregs in the primary CRC, but in liver metastases, the abundance of Tregs were comparable between treated and untreated tumors (Fig. 5f). The differentiation trajectories of CD8 + and CD4 + T cells (Fig. 2e, 2f, Figure S4e, S4f) also con rm the CD8 + dysfunctional T cells were most prevalent in treatment-naïve tumors, while the Tregs is shared in treated and untreated tumors.
6. Cell-cell cross-talks within the tumor microenvironments in liver metastasis of CRC Tumor microenvironment are a complex ecosystems, the crosstalks between each cell determine tumor biology and response to therapies (Tirosh and Suvà, 2019). In order to systematically map cellular interactions especially different immune cells and mesenchymal cells in the TME of the primary CRC and metastases, and investigate the potential cellular communication which contribute to the cancer progression, metastases and immune evasion, we used CellPhoneDB v2.0.6 to study the cross-talks between stromal cells in tumor microenvironment. To visualize the cross-talks between different cells types, a chord diagram was built using the circlize package (Gu et al., 2014) in R.
Firstly, we provided a landscape of cross-talks within major stromal cell types both in the niche of the primary CRC and the liver metastases. Then based on each sub cell types annotated above, we focused on sub-cell types, investigated intercellular communications in tumors treated with preoperative chemotherapy and untreated tumors. Selected LR paired were summarized in Fig. 6. The full list of results that were unique to different tumor microenvironment were available in Supplementary Tables (Table S4).
The cross-talks within major stromal cells reveals that TAMs had the most broad cross-talks with other cells, both in the niche of primary CRC and liver metastases ( Figure S6a, S6b). They displayed a rich ligands and receptors (LR) pro le, broadly communicated with mesenchymal compartment (CAFs and

endothelial cells), and immune compartment (including T cells, NK cells, mast cells, DCs and TAMs) in
both the primary and liver metastases ( Figure S6a). Comparison of the interactions in the primary CRC and liver metastases, we found that TAMs in the primary CRC communicated more frequently with CAFs than that in liver metastases, whereas the cross-talk between TAMs and DCs was greatly increased in the metastatic niche when compared with that in primary CRC ( Figure S6a, S6b).
Then on the basis of the annotation above, we focused on cell subpopulations, to identify key interactions between each cell subpopulations.The cross-talks between different phenotypes of TAMs were most abundant in the TME, compared with other cell types. Subsequently, the cross-talks between sub populations of DCs, CAFs, TAMs, dysfunctional T cells and endothelial cells were relatively more frequent. Whereas immune cell such as B cell, plasma cells, other T cell subtypes (including T EM , T RM , T EMRA and Tregs), NK cells and mast cells communicated less in TME ( Figure S6c).
Interestingly, among the TAMs compartments, monocyte-derived immature TAMs, MKI67 + TAMs and M03-CXCL10 + TAMs were more activated, with more interactions with other cell types ( Figure S6c), hinting that they mediate signals between diverse cell types. Moreover, compared with other T cell subtypes, we observed MKI67 + CD8 + T cells and T EX displayed a rich LR pro le ( Figure S6c Cell-Cell communications in the niches of tumors treated with and without preoperative chemotherapy Whether preoperative chemotherapy can reprogram the interactions within stromal cells is still unclear. To investigate the effects of chemotherapy on tumor microenvironment of liver metastases of CRC, we further compared cell interaction networks between tumors treated with PC and without PC both in primary tumor and metastatic sites.
Our data reveal the phenotypes of TAMs were highly heterogeneous, and the dysfunctional T cells and ECM-remodeling CAFs were expanded in the primary CRC of treatment-naïve tumors. While in the microenvironment of tumors treated with preoperative chemotherapy, in the primary CRC, chemotherapy promoted the activation of B cells, increased the abundance of immature and less activated TAMs, and expanded the myo broblasts.
Moreover, preoperative chemotherapy promoted the activation of B cells. In the niche of primary tumors treated with PC, B cells, MKI67 + T cells (T24) and dysfunctional T cells expressed HLA-F, and/ or HLA-DPB1, interacted with PC-enriched TAMs through HLA-F -LILRB2 and HLA-DPB1 -NRG1 (Fig. 6a, right).
In liver metastases, chemotherapy promoted the abundance of DCs and myo broblasts, but decreased the dysfunctional T cells.
In tumors treated with PC, in line with the observations in the primary tumors, Notch signaling was also identi ed in the liver metastases after tumors treated with chemotherapy (Fig. 6b, left). Importantly, TNF signaling was unique present in tumors treated with PC. Comparison of tumors after treatment, in treatment-naive tumors, cross-presenting DCs expressed DPP4, interacted with TAMs through CXCL2 -DPP4, CXCL10 -DPP4 and CCL3L3-DPP4 (Fig. 6b, left). DPP4 (also known as CD26) has been showed to be positive correlated with distance metastasis in colorectal cancer (Lam et al., 2014), and the high expression of DPP4 had signi cantly worse overall survival in CRC (Lam et al., 2014). In addition, Cross-presenting DCs also interacted with TAMs through LGALS9-HAVCR2, implying their roles in immunosuppression.
In summary, our LR interaction map highlights ACKR1, Notch signaling and molecules mediated immune regulation may contribute to the reprogramming of TME after tumors treated with preoperative chemotherapy.

Discussion
How to improve therapeutic options for patients with metastatic CRC is a core question for CRC treatment. In this study, we performed a single cell atlas of primary CRC and their matched liver metastases with 111,292 cells, provides a fundamental and comprehensive understanding of cellular composition in TME of liver metastases of CRC. More importantly, for the rst time, we provided a dynamic and comprehensive stromal cell mapping to illustrate how preoperative chemotherapy reprogramed the TME of both primary CRC and matched liver metastases in CRC patients.
We nd that B cells mainly existed in primary CRC. Importantly, our results indicated preoperative  (Fig. 2i, Fig. 2j). Thus, our results suggested that the in ltration of the switched, activated B cell may play an important role in anti-tumor response, and they could be a potential predictor of effective chemotherapy and good prognosis of CRC.
In this study, we identi ed 18 clusters of TAMs, and classi ed them into four major heterogeneous subclasses based on gene expression pro les (Fig. 3c) The subset of in ammatory TAMs had a strong in ammatory phenotype and may contribute to tumor progression seriously. Taken together, TAMs in TME of treatment naïve tumors may be tumor-promoting. Whereas, the vision of TAMs in of tumors treated with preoperative chemotherapy were very different. We found immature TAMs and less activated TAMs were more enriched. Specially, gene signature of immature TAMs (M11) is associated with better prognosis in MSS cohort. According to recent studies, high TAMs density were closely correlated with poor survival in many cancers (Gabrilovich, 2017;Ruffell and Coussens, 2015). However, TAMs targeted therapy has not been effective yet, which may be hampered by TAMs heterogeneity and elusive molecular phenotype. Our study provided a full-scale illustration of TAMs composition and molecular characteristics, which were an important basic and resource information for TAMs targeted therapy research to inhibit and cure metastatic CRC. Our results also found preoperative chemotherapy suppressed the abundance of dysfunctional T cells and ECM-remodeling CAFs, induced the generation of myo broblasts. The accumulation of myo broblasts indicated tissue injury and brosis related to chemotherapy.
Therefore, based on the phenotypic alteration we observed in the primary CRC treated with and without preoperative chemotherapy, we suggested that chemotherapy destroyed tumor cells, which released tumor antigens and activated immune microenvironment (Fig. 7), including B cell maturation and antibody generation. Moreover, chemotherapy inhibited the diversity of TAMs and remolded the character of TAMs, which transformed into an immature and less activated phenotype. While the TME of liver metastasis displayed clearly different scenario, the number of B cells in the liver metastasis had a signi cant reduction, which might be a cause for the feasibility of liver metastasis of CRC. Similar to the situation of primary CRC, myo broblast were aggregated in response to chemotherapy, which led to the brosis of liver metastasis.
Taken together, this atlas provides a fundamental reference for the complex cellular and phenotypic diversity within TME of both primary CRC and liver metastasis. Our systematic investigation of transcriptional changes and phenotypic alteration in TME at single-cell level can provide valuable insight into our understanding of therapeutic outcome. This may open up new possibilities to develop or improve therapeutic strategy to CRC.

Tumor specimens and Patient Clinical Characteristics
Primary CRC, matched liver metastases and blood sample were collected from six CRC patients with metastatic disease. All patients were classi ed as microsatellite-stable (MSS) with invasive adenocarcinomas and late-stage (IV) disease (Table S1). For each lesion, we collected the tissue in the core of the tumor after the surgery. Single cells were isolated from fresh tumor tissues without surface marker pre-selection. All patients underwent curative intent surgery of synchronous colectomy with liver resection. In addition, all patients who provided specimens signed an informed consent form and agreed to the specimens being used for scienti c research, and more pathological and clinical information of patients see in Table S1. Among the patients, patients COL15, COL 17 and COL18 treated with preoperative chemotherapy, and the others were treatment naïve. More details could see in Table S1.
Tumor disaggregation and single cell collection Venous blood was collected before surgery in EDTA anticoagulant tubes and used to isolate PBMC immediately. Fresh biopsies of the primary colorectal cancer and the matched liver metastases were collected at the time of surgical resection. Once the specimens were separated from body, they were processed for single cell RNA sequencing (scRNA-seq) immediately.

PBMC Isolation
Peripheral blood mononuclear cells (PBMCs) were isolated using Ficoll (TBD) solution according to manufacturer's instructions. In brief, 5ml fresh peripheral blood was layered onto equal Ficoll, following by centrifugation at 450x g for 25min. After centrifugation, lymphocytes layer remained at the plasma-Ficoll interface and were carefully transferred to a new tube and washed twice with Phosphate-Buffered Saline (PBS, ThermoFisher Scienti c). Lymphocyte pellets were re-suspended with sorting buffer (Hank's Balanced Salt Solution (HBSS, ThermoFisher Scienti c) with 0.04% bovine serum albumin (BSA, MRC)) for ow cytometry analyses.

Tumor Dissociation
The primary CRC and metastatic tumor tissue were dissociated using MACS® Tumor Dissociation Kit (low) cell cluster. For PBMC sample, only lymphocytes and monocytes clusters were gated for further sorting. Viable single cells were resuspended in HBSS with 0.04% BSA. Viability was con rmed to be >90% using trypan blue (ThermoFisher Scienti c) exclusion prior to scRNA-seq process.

Droplet-based sing-cell RNA sequencing and library preparation
The scRNA-seq libraries were constructed by using the Chromium™ Single Cell 3' Reagent Kits v2 (10x Genomics) according to manufacturer instruction. Brie y, cells were suspended in HBSS with 0.04% BSA at a concentration approximately 1000 cells/μl and appropriate suspension loading volume were determined by calculating for a target capture 8000 cells. Cell suspension of corresponding volume were loaded onto the 10X Chromium Single Cell Platform (10X Genomics). Generation of gel beads in emulsion (GEMs), barcoding, GEM-RT clean-up, complementary DNA ampli cation and library construction were all performed according to manufacturer's protocol. Sequencing library quality was checked with Bioanalyzer (Agilent Bioanalyzer 2100). Library quanti cation was measured using Qubit before pooling. The nal library pool was sequenced on the Illumina X Ten instrument using 150-basepair paired-end reads. Sequencing data of individual sample were summarized in Table S2.
Preprocessing of single-cell RNA-seq data analysis The raw base call (BCL) les were demultiplexed into FASTQ le by bsl2fastq. Droplet-based sequencing data were quali ed by FastQC software. Then reads were aligned against GRCh38 human reference genome provided by Cell Ranger (version 2.0, 10x Genomics), unique molecular identi er (UMI) counts were summarized for each cell of each gene. The raw UMI count matrices were ltered to (1) remove cells with a low number of unique detected genes (<500); (2) for each batch, remove cells for which the total number of UMI (after log 10 transformation) was not within the three standard deviations of the mean; (3) for each batch, remove cells that showed an unusually high or low number of genes; (4) discard cells in which the proportion of the UMI count attributable to mitochondrial genes was greater than 15%. 731 cells were ltered in step1, while step 2 to step 4 remove only a small number of cells (0.1%). After exclude low-quality cells, 25,121 protein-coding genes across 111,292 single cells remained for downstream processing.

Identi cation of cell types and subtypes by dimensional reduction
After quality control, raw UMI counts were lognormalized using the scale of 10,000. The genes with normalized expression between 0.0125 and 3, and dispersion >0.5 were selected as highly variable genes. 1511 highly variable genes were identi ed based on dispersion and mean. "var.to.regress" option UMI's and percent mitochondrial content were used to regress out unwanted sources of variation. The resultants were rst summarized by principle component analysis (PCA). We used the function FindClusters on 50 PCs with resolution 1.0 to perform the rst-round cluster and annotation. The annotation of each cell cluster was con rmed by the expression of many canonical marker genes. As shown in Figure S1b

Pathway enrichment analysis
We used the R package limma (Ritchie et al., 2015) to identi ed differentially expressed genes between cells from treatment-naïve tumors and chemotherapy-treated tumors. The Benjamini-Hochberg multiple testing correction was applied to estimate the FDR.
Biological-process Gene ontology (GO) enrichment (P<0.01) were performed using clusterPro ler packages (version 3.9.2) with a Benjamini-Hochberg multiple testing adjustment. Gene sets with a FDRcorrected P<0.01 were considered to be signi cantly enriched.

Putative Interactions between Cell Types
We used CellPhoneDB v2.0.6 (www.cellphonedb.org) to study the crosstalk between stromal cells in tumor microenvironment. CellPhoneDB is a repository of curated receptor, ligand and their interactions to predict communicating pairs. They integrated multiple databases and account for subunit architecture

Multicolor immunohistochemistry (IHC)
Multicolor IHC staining of formalin-xed, para n-embedded (FFPE) tissues section were used to con rm the presence of novel subpopulations. Multicolor IHC staining was obtained using PANO 7-plex IHC kit (0004100100, Panovue). We detected the expression of CD3, CD8, PD-1. Each formalin-xed, para nembedded (FFPE) tissues section (4μm) was put through four sequential rounds of staining. Brie y, FFPE tissue section (4μm) were melted at 60℃ for 1 hour followed by depara nizing and rehydrating. Heat mediated antigen retrieval was performed in EDTA buffer pH 9 using microwave incubation. The section were blocked with blocking buffer (hydrogen peroxide ) for 10 minutes. In the rst round of staining, the section were stained with an anti-PD-1 mouse monoclonal antibody (ZM-0381, Zsbio). Then, they were incubated with an anti-rabbit and mouse horseradish peroxidase (HRP)-conjugated secondary antibody (0004100100, Panovue). The stained signal was further ampli ed using PPD 520 tyramide signal ampli cation (TSA) reagents through incubated with TSA diluent (0004100100, Panovue). After removing the bound antibody by a method similar to the antigen retrieval process using EDTA buffer pH 9, the same process was repeated in sequence for then following antibodies and uorescent dyes: anti-CD8 (ZM-0508, Zsbio) /PPD560, anti-CD3 (85061, CST) /PPD690. Nuclei were stained with DAPI after all the staining had done.

Multispectral imaging
The obtain multispectral images, the stained slides were scanned using the Mantra system (Perkin Elmer), which captures the uorescent spectra at 20-nmwavelength intervals from 420 to 720 nm with identical exposure time. The scans were combined to build a single stack image. Images of unstained and single-stained sections were used to extract the spectrum of auto uorescence of tissues and each uorescein, respectively. The extracted images were further used to establish a spectral library required for multispectral unmixing by InForm image analysis software (PerkinElmer). Using this spectral library, we obtained reconstructed images of sections with the auto uorescence removed.          Wilcoxon rank-sum test was used for statistical analysis. e, GO analyses for genes that were differentially expressed between CAFs from tumor treated with PC versus treatment-naive tumors. Benjamini-Hochberg corrected P < 0.01.

Figure 7
Diagram illustrating the reprogram of TME to response to preoperative chemotherapy.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.