Spatial gene analysis
To explore the spatial gene expression, we performed the ST sequencing of two primary colorectal cancer tissues and two matched liver metastatic tissues after simultaneous colorectal and liver resection. Cryosections of the four tissues were mounted on spatially barcoded ST microarray slides. Before sequencing, H&E staining and brightfield imaging were performed, and distinct histological features were annotated by professional pathologists. The samples were processed for ST analysis, including complementary DNA synthesis, amplification by in vitro transcription, library construction, and sequencing (Figure 1).
Transcriptome heterogeneity for matched primary tumor and liver metastasis
To explore the ITH of primary CRC tissues and their matched metastatic tissues, we drew a boundary to distinguish the malignant tissues from the nonmalignant tissues (Figures 2a, 2b). Subsequently, we performed ST sequencing across the four tissues. The ST data consisted of 2700-4300 spots per tissue, and each spot contained approximately 8000-15000 UMIs and 3200-4600 genes (Figures 2e, 2g). After clustering and uniform manifold approximation and projection analysis, we identified a total of 19 clusters from the four tissues (Figure 2c). ST lacks cellular resolution, and we used tissue-specific gene expression analysis to further determine the annotation of spots. MUC1 and PIGR encode two epithelial-related plasma membrane proteins rarely expressed in intestinal tissue. FGA and APOB are two liver-specific markers, while OLFM4 is a cancer stem cell marker and MYC is a classic proto-oncogene. The expression of PTPRC is selective in the cytoplasm of lymphoid tissue and immune cells (Figures 2f, 2h).
In the cancerous regions, the ST data showed apparent transcriptome heterogeneity in both primary and metastatic tumors, with three clusters of spots in primary tumors and five clusters in metastatic tumors (Figure 2c). Next, we performed differential gene expression analysis in all clusters to obtain a detailed view of the heterogeneity (Figure 2d). Cluster 3 variably expressed COL12A1, MMP1, and KRT17. Cluster 4 showed the upregulation of CD24, TNNC2, and RUBCNL, and cluster 16 differentially expressed RPL36A, RPL22A, and ZFAS1. Uniform manifold approximation and projection analysis showed that clusters 2, 6, 11, 12, and 13 were colorectal liver metastatic tissues; however, clusters 2 and 12 were consistent with the observed cancer cells based on the H&E staining, and the remaining clusters were stromal cells. Cluster 2 demonstrated the upregulation of DEFA5, REG1A, and SPP1, and cluster 12 variably expressed DAPL and AC078993. Interestingly, the stromal region was heterogeneous in the metastatic tissues. Cluster 6 was mainly associated with immune-related genes such as IGHG2, IGLC3, and IGHG1. Cluster 11 was enriched with MGP, AEBP1, and IGFBP7, whereas cluster 13 was differentially expressed CHIT1, APOE, and GPNMB. Combined with the spatial information, we observed that clusters 2 and 12 were surrounded by cluster 6, where immune-related genes were enriched. Considering this phenomenon, metastatic tumors may be infiltrated with more immune cells, and the crosstalk between immune cells and metastatic cancer cells may facilitate metastasis. Taken together, the ST maps provided a global view of ITH in the primary tumor and its paired liver metastases.
Spatial expression patterns of cancer stem cell markers revealed the enrichment of stem cells in metastatic tumors
The development of cancer is a multi-step process and involves unlimited replicative potential and acquired invasive and metastatic abilities, among others. We have identified ITH in primary tumor and liver metastases, and each subgroup of heterogeneous cells can have distinct involvements in tumor progression, drug resistance, and metastasis. To further analyze the unique biological functions of each cluster, we evaluated the expression of proliferation-related genes in each cluster. The results showed that clusters 6, 9, 10, 11, 13, 14, 15, 17, and 18 had relatively higher G1 phase scores and, thus, the least proliferation ability (Figure 3a). Epithelial-mesenchymal transition is another feature of cancer progression, during which epithelial cells acquire features of mesenchymal cells for metastasis. We initially evaluated the expression of sets of genes related to epithelial-mesenchymal transition to identify the clusters that have the highest possibility of metastasis (Figure 3a). However, considering that the boundaries of spots have no natural correspondence to the boundaries of cells, the gene expression profiles of spots can encompass multiple cells and/or portions of cells rather than truly resolving individual cells. Clusters with the highest p-EMT scores may contain a higher proportion of mesenchymal cells than those with the highest metastatic ability. Thus, we chose clusters 2, 4, 12, and 16 for further analysis, as these four clusters had the highest proliferation ability and abundance of epithelial proportions, which indicated that these clusters had the highest proportion of tumor cells.
In primary tumors, 508 genes (e.g., PIGR, IGHA1, IGHA2, IGKC, and IGHG4) were upregulated in cluster 4 while 282 genes were upregulated in cluster 16 (e.g., SLC25A5, TXN, HSPE1, SNHG6, and RPS15A) (Figure 3b). In the metastatic tumors, 783 genes were upregulated (e.g., BGN, MT2A, VIM, APOE, and APOC1) in cluster 2, while 151 genes (e.g. AGR2, SLC6A2, EPCAM, VAV3, and FOXD1) were upregulated in cluster 12 (Figure 3b). When the details of differentially expressed genes were analyzed, several cancer stem cell (CSC) markers in cluster 12, including EPCAM (19), LGR5 (20), OLFM4 (21) and FOXD1 (22), were upregulated, which suggested that CSCs were enriched in cluster 12 (Figure 3b).
CSCs belong to a subpopulation of tumor cells capable of self-renewal, multilineage differentiation, and survival in adverse microenvironments and are key to tumor initiation and tumor relapse (23, 24). CSCs were identified based on the expression of specific markers. To further explore the global distribution of CSCs in cancerous tissues using ST, we investigated a set of stem cell-related genes in all clusters. BMI1, HOPX, LRIG1, and LGR5 are recognized for their ability to self-renew and differentiate into normal intestinal epithelium (25). The violin plots indicated that cancerous clusters (clusters 2, 3, 4, 6, 11, 12, 13, and 16) had higher expression of intestinal stem cell markers than other regions of the four tissues (Figure 3c). A group of CSC markers (ALCAM, CD44, EPCAM, EPHB2, OLFM4, and PROM1) (19, 26) were analyzed, and the results demonstrated that CSC markers were significantly enriched in cancerous clusters. Surprisingly, metastatic clusters (clusters 3, 4, and 16) also showed a higher expression of these markers than the primary clusters (Figure 3c). DPP4 (27) and PROX1 (28) have been reported as metastatic CSC markers, and our results showed an increase in the expression of these two markers in metastatic clusters (Figure 3c). We identified FOXD1 as another metastatic CSC marker because of its similar expression as DPP4 and PROX1 (Figure 3c).
Trajectory and pseudo-time analysis revealed a dedifferentiation-differentiation process during metastasis
Cellular plasticity allows cancer cells to shift between differentiated and undifferentiated states, and this contributes to cancer progression (29, 30). To gain insights into the development of liver metastasis, we performed trajectory analysis of eight cancerous clusters, and the results revealed two distinct trajectories (Figure 4). To quantify the progression through differentiation, we calculated the pseudo-time for each branch and observed three states of differentiation (Figure 4). Integrating the pseudo-time and spatial information, most spots of cluster 12 were collected at the top of the evolution tree; this indicated more dedifferentiation (Figure 4), which was consistent with the high abundance of CSCs in cluster 12. The bifurcation spots from metastatic and primary tumors were clustered at each branch (Figure 4), suggesting that the tumor cells in primary tumors undergo dedifferentiation to metastasize before differentiating into seeds in the liver. Figure 4b shows that the expression of genes significantly varied over pseudo-time, and the heatmap demonstrated six clusters based on the expression changes. Next, we analyzed the expression patterns of several stem cell markers over pseudo-time, including ALCAM, FOXD1, BMI1, HOPX1, CD44, LGR5, EPCAM, and OLFM4 (Figure 4c). Apart from the low expression of ALCAM over the pseudo-time, we observed three different expression patterns. FOXD1, BMI1, and LGR5 were highly expressed during early differentiation; HOPX and CD44 were upregulated midway through pseudo-time, while EPCAM and OLFM4 were downregulated midway through pseudo-time. Moreover, the expression patterns of FOXD1 and LGR5 were almost identical, implying crosstalk between FOXD1 and LGR5. Therefore, we performed a correlation analysis of FOXD1 and LGR5, and the results revealed that the expression of FOXD1 was positively correlated with LRG5 (Figure 4d). We observed the upregulation of FOXD1 and LGR5 in cluster 12. Therefore, we conducted a functional analysis of cluster 12, and the Hippo signaling pathway was significantly enriched (Figure 4e). Compared with other clusters, the Hippo signaling pathway was hyper-activated in clusters 2 and 12 (Figure 4f).
Dominant crosstalk of CD74-MIF between clusters in metastatic tumors.
Cell-cell communication is fundamental to various processes in CRC, such as cell fate decisions, proliferation, migration, and homeostasis (31). For deeper insights into the direct interaction networks of cancer cells, we performed CellPhoneDB analysis of the cancerous clusters in primary and metastatic tissues. Overall, we observed a significant increase in the cell-cell interactions in metastatic tissues relative to the primary tissues (Figure 5a), suggesting a more complex cell-cell communication in the development of metastasis. In the primary tumor, the dominant interaction was between C5AR1 and RPS19 in clusters 3 and 4. RPS19 is known to be associated with a high incidence of cancer in humans and upregulated in colon cancer (32). RPS19 interacts with C5aR1 expressed on tumor-infiltrating myeloid-derived suppressor cells (MDSCs), facilitating MDSC recruitment in tumors to promote tumor growth (33). In metastatic tissues, communication between CD74 and MIF was the most dominant. MIF secreted by the cells of cluster 12 interacted with CD74 in clusters 2, 3, 6, 11, and 13. MIF is a multifunctional cytokine and a ligand of the CD74 receptor (34). Activated by the secreted MIF from cluster 12, the MIF/CD74 axis in clusters 2, 3, 6, 11, and 13 could be initiated to promote cell survival by modulating glucose uptake and metabolism through the activation of AKT and AMPK proteins (35), which was verified by KEGG pathway analysis. However, the reserve communication between MIF and CD74 among the clusters was not observed.
FOXD1 was associated with cancer cell stemness, and the protein-protein interaction network showed that FOXD1 interacted directly with WNT9B, which functions in the canonical Wnt/β-catenin signaling pathway. To further examine the clinical implications of FOXD1, we explored the TCGA dataset and found that patients with lower expression of FOXD1 demonstrated better overall survival than those with higher expression of FOXD1. Multivariate analysis verified the independent protective role of FOXD1 in surviving patients (AUC = 0.6725) and its usefulness for predicting clinical survival with satisfactory accuracy.