Single-cell genomic and transcriptomic landscapes of primary and metastatic colorectal cancer tumors

Abstract

Background: Colorectal cancer (CRC) ranks as the second-leading cause of cancer-related death worldwide with metastases being the main cause of cancer-related death. Here, we investigated the genomic and transcriptomic alterations in matching adjacent normal tissues, primary tumors, and metastatic tumors of CRC patients.

Methods: We performed whole genome sequencing (WGS), multi-region whole exome sequencing (WES), simultaneous single-cell RNA-Seq, and single-cell targeted cDNA Sanger sequencing on matching adjacent normal tissues, primary tumors, and metastatic tumors from 12 metastatic colorectal cancer patients (n=84 for genomes, n=81 for exomes, n=9120 for single cells). Patient-derived tumor organoids were used to estimate the anti-tumor effects of a PPAR inhibitor, and self-renewal and differentiation ability of stem cell-like tumor cells.

Results: We found that the PPAR signaling pathway was prevalently and aberrantly activated in CRC tumors. Blocking of PPAR pathway both suppressed the growth and promoted the apoptosis of CRC organoids in vitro, indicating that aberrant activation of the PPAR signaling pathway plays a critical role in CRC tumorigenesis. Using matched samples from the same patient, distinct origins of the metastasized tumors between lymph node and liver were revealed, which was further verified by both copy number variation and mitochondrial mutation profiles at single-cell resolution. By combining single-cell RNA-Seq and single-cell point mutation identification by targeted cDNA Sanger sequencing, we revealed important phenotypic differences between cancer cells with and without critical point mutations (KRAS and TP53) in the same patient in vivo at single-cell resolution.

Conclusions: Our data provides deep insights into how driver mutations interfere with the transcriptomic state of cancer cells in vivo at a single-cell resolution. Our findings offer novel knowledge on metastatic mechanisms as well as potential markers and therapeutic targets for CRC diagnosis and therapy. The high-precision single-cell RNA-seq dataset of matched adjacent normal tissues, primary tumors, and metastases from CRCs may serve as a rich resource for further studies.

Keywords: Genotype-phenotype relationship; Lineage tracing; Metastatic colorectal cancer; Mitochondrial mutations; PPAR signaling pathway; Patient-derived organoids; Single-cell transcriptome profiling.

Fig. 1

Single-cell transcriptome analysis of colorectal cancer. A The workflow illustrates the strategy for cell collection from matching adjacent normal tissues and primary and metastatic colorectal tumors for single-cell RNA-Seq, single-cell cDNA Sanger sequencing, and bulk level whole genome sequencing and whole exome sequencing. B Bar plot showing the number of cells collected from each patient. Patients are ordered according to the total number of cells. Color represents cell origin. Low-quality cells are removed with strict criteria. Only cells that expressed at least 1000 genes and showed second maximum pairwise Pearson correlations greater than 0.6 were retained for further analysis. Finally, 8085 cells (88.7%) were retained for subsequent analyses. The details can be found in the “Methods” section. PT: primary tumors. LN: adjacent normal tissue collected from liver. N: adjacent normal tissues. MT: metastatic tumors, including lymph node metastasis, liver metastasis, and omentum metastasis. Removed: cells that have not passed the quality control and not used for subsequent analysis. C UMAP plot of cell clusters. Cell types were identified base on the regulon activity matrix and then visualized by UMAP. Cells were colored according to annotated cell types. According to their expression of known marker genes shown in Supplementary Figure 1B, we annotated these clusters as epithelial cells, endothelial cells, fibroblasts, T-cells, B-cells, pre-B-cells, macrophages, and mast cells. Associated with Supplementary Figure 1A and Supplementary Figure 1C-D. Most of the immune cells, fibroblast, and endothelial cells come from the tumor area. D Immunofluorescence staining of the shared endothelial and fibroblast marker SPARC and T-cell marker CD3D in adjacent normal tissues, primary tumor, and liver metastasis. The boxplot shows the SPARC- or CD3D-positive cell ratio in different regions (N: adjacent normal tissue; PT: primary tumor; LM: liver metastasis). Scale bar, 100μm. E The dot plot shows the ratio of cells that highly expressed enterocyte marker (CA2) and intestinal stem cell marker (SOX9) in normal and tumor regions for each patient. F The expression levels of enterocyte marker (CA2) and intestinal stem cell marker (SOX9) were projected on epithelial cells tSNE maps. Colors from yellow to red represent expression levels from low to high

Fig. 2

SOX9/MKI67-positive cells may have potential of self-renewal and differentiation. A Immunofluorescence staining of cell proliferation marker (MKI67) and intestinal stem/progenitor cell marker (SOX9) for adjacent normal tissue and tumor tissue in vivo as well as patient-derived cancer organoids that were cultured for 2 weeks and 2 months. CO: cancer organoid. 2W: 2 weeks. 2M: 2 months. Scale bar, 100 μm. B Dot plot showing the expression level of MKI67 and SOX9 of patient #1. Colors represent cells that collected from different regions. N, normal region. PT, primary tumor. LM, liver metastasis. PT organoid, in vitro cultured organoid that are derived from the primary tumor. The proportion of double positive cells is written in the upper-right corner of each dot plot. C Bar plot shows the ratio of SOX9/MKI67 double positive cells of patient #1 from different regions. N, normal tissue; PT, primary tumor; LM, liver metastasis; PT organoid, primary tumor-derived organoid. D The bar plot shows the ratio of cells that expressed different level of enterocyte marker CA2. D0 represents thawed tumor organoid, and 1mon represents tumor organoid that is cultured for 1 month. P-value was calculated through t-test. The color of the bar represents the expression levels (log₂(TPM+1)) of CA2. E Immunofluorescence staining of enterocyte marker CA2 on long-term cultured tumor organoid (O#H). Scale bar, 50μm. F Heatmap shows cell-type specific genes expression in one of five tumor organoids of one patient (O#S5). Colors from blue to red represents expression level (log₂(TPM+1)) from low to high. The black boxes highlight cells that do not express or low express stem/pluripotency markers (SOX9, OLFM4, LGR5, ALCAM, LRIG1, and MKI67) but highly express differentiated marker CA2

Fig. 3

Transcriptomic differences between adjacent normal tissues and tumor tissues. A Heatmap showing the differentially expressed genes (DEGs) between adjacent normal and tumor tissues in each patient. Colors from blue to orange represent the expression level from low to high. The enriched gene ontology (GO) terms and p-value were showed in the right side of the heatmap. The full list of DEGs is summarized in supplementary table 2. B Immunofluorescence staining of CEACAM6 and LY6E in adjacent normal tissue and primary tumor tissue for another two patients. Both CEACAM6 and LY6E are highly expressed by tumor epithelial cells. Scale bar, 100 μm. C Bright-field images of patient-derived normal and tumor organoids (left panel). Representative bright-field images of PPAR inhibitor (FH535)-treated and untreated organoids (right panel). D Cell viability assay of PPAR signaling pathway inhibitor (FH535)-treated and untreated organoids. ** represents p-value < 0.05. E Cell apoptosis assay of PPAR inhibitor (FH535)-treated and control organoids. The results of fluorescence-activated cell sorting (FACS) are shown. Q2 represents late apoptotic cells, Q3 represents early apoptotic cells, and Q4 represents living cells. F Bar plot shows the ratio of apoptotic cells, light gray represents early apoptotic cells, and dark gray represents late apoptotic cells. P-values are calculated by t-test. ** represents p-value < 0.01. G Dot plot showing the proliferation ratio of cells. Different colors represent different drug treatment. 5-FU: 30 μM; GW9662: 30 μM; FH535: 15 μM. P-values are calculated by t-test. ** represents p-value < 0.01. NS represents not significant. There are 3 replications for each treatment. H The line plot showing the cell viability of tumor organoid (O#L) under different concentrations of drugs. The x-axis represents the concentration of PPAR inhibitor GW9662, and the y-axis represents the cell viability of tumor organoid with 5 days drug treatment. Colors from blue to red represent the concentration of WNT inhibitor XAV939 from low to high

Fig. 4

Tumor clonal architectures and tumor lineage inferences. A CNV pattern and sampling strategy of patient #9. Different subclones with distinct CNV patterns are found in lymph node and liver metastases. N: adjacent normal tissues. PT: primary tumors and R1-R3 represent different regions of primary tumor. LyM: lymph node metastasis. LM: liver metastasis and LM1 and LM2 represent two separated liver metastatic tumors. The pie charts reflect the proportion of different subclones in each region. Copy number gain and copy number loss are indicated with red and blue respectively. Different tissues are divided by bold black solid lines, and different clones within the same tissue block are divided by black dashed lines. The squares above the heatmap represent different chromosomes, black squares represent odd-numbered chromosomes and chromosome X, and light gray squares represent even-numbered chromosomes and chromosome Y. The squares on the left of the heatmap represent different subclones, and the same color represents the same clone. The most obvious copy number differences between different clones are highlighted by different colored boxes. Red squares in chromosome 2 show the additional CNVs of metastatic tumor (A2 and B2 subclones) compared to the primary tumor (A1 and B1 subclones). Pink square in chromosome 8 shows D subclone-specific CNV. Black squares in chromosome 18 show the A (A1 and A2) and C (C1 and C2) clones specific CNVs. Green squares in chromosomes 20 and 21 show C (C1 and C2) clone-specific CNVs. The numbers next to the heatmap show the number of cells of each subclone in different tissue blocks. B The diagram showing the tumor metastasis path of patient #9. The color of box represents different areas; blue, orange, and red boxes represent the primary tumor, lymphatic and liver metastasis tumor respectively. The circle and triangle represent the mutation state at the position 2897- and 1350-point site of mitochondria. Their colors represent the mutation status: blue, orange, and red represent wild-type, heterozygous, and homozygous mutations respectively. PT: primary tumors and R1-R3 represent different regions of primary tumor. C Heatmap showing selected mitochondrial mutations of patient #9. Orange represent heterozygous mutations and red represents homozygous mutations. Blue represent wild-type and gray represents read depth lower than 9. Cells from PTR1 and liver metastasis have a chrM:2897 heterozygous mutation. PTR2 has region-specific mutations on chrM:11380 and chrM:8534. The bars above the heatmap shows the CNV subclones and tissue origin of the cells. The full list of mitochondrial mutations can be found in supplementary table 5. PT: primary tumors and R1-R3 represent different regions of primary tumor. LyM: lymph node metastasis. LM: liver metastasis and LM1 and LM2 represent two separated liver metastatic tumors. D Heatmap showing the reginal distribution of somatic mutations in all samples from patient #9. In total, 136 somatic mutations that were shared by at least two samples were identified. Red represent mutant state and gray represent wild-type state (left panel). Phylogenetic tree of lesions of patient #9 based on somatic mutations calling by WES data and the pie charts reflect the proportion of different subclones in each region based on CNVs that inferred by scRNA-Seq data

Fig. 5

Possible presence of tumor precursor cells. A CNV pattern and sampling strategy of patient #8. N: adjacent normal tissues. PT: primary tumors and R1-R4 represent different regions of primary tumor. LyM: lymph node metastasis LyM1 and LyM2 represent two separated lymph node metastatic tumors. LM: liver metastasis. R1-R2 represents different regions of the liver metastatic tumors. The pie charts reflect the proportion of different subclones in each region. Copy number gain and copy number loss were indicated with red and blue respectively. Different tissues are divided by bold black solid lines, and different clones within the same tissue block are divided by black dashed lines. The squares above the heatmap represent different chromosomes, black squares represent odd-numbered chromosomes and chromosome X, and light gray squares represent even-numbered chromosomes and chromosome Y. The squares on the left of the heatmap represent different subclones, and the same color represents the same clone. The most obvious copy number differences between different clones were highlighted by different colored boxes. B The diagram showing the tumor metastasis path of patient #8. CNVs, somatic mutations and mitochondrial mutations are also shown in the diagram. The color of box represents different areas; blue, orange, and red box represent the primary tumor, lymphatic, and liver metastasis tumor respectively. Most of primary tumor cells have homozygous mutation at MT-ND5:13,536G>A mutation; thus, we speculated that there may be a group of tumor progenitor cells that we did not capture, which have MT-ND5(13,536G>A) mutation. With the development of the tumor, this group of tumor progenitor cells have additional mutations at MT-ND6(14,504 T>TA) and MT-ND3(10,396T>C) respectively, thus forming two clones (PTR1 A clones and PTR2-PTR4 B&C clones). Then the clones with MT-ND3(10,396T>C) mutation were further metastasized to the liver. C Heatmap showing selected mitochondrial mutations of patient #8. Almost all tumor cells have a chrM:13,526 homozygous mutation. Cells from the PTR1 (primary tumor region 1) have a chrM:14,504 mutation, and cells from PTR2-PTR4 and lymph node metastasis have a chrM:10,396 mutation. Blue represents wild-type and gray represents read depth lower than 9. Orange represent heterozygous mutation and red represent homozygous mutation. The bars above the heatmap shows the CNV subclones and tissue origin of the cells. D CNV pattern and sampling strategy of patient #1. PT: primary tumors. LM: liver metastasis. Copy number gain and copy number loss were indicated with red and blue respectively. Different tissues were divided by bold black solid lines, and different clones within the same tissue block were divided by black dashed lines. The squares above the heatmap represent different chromosomes, black squares represent odd-numbered chromosomes and chromosome X, and light gray squares represent even-numbered chromosomes and chromosome Y. The squares on the left of the heatmap represent different subclones, and the same color represents the same clone. The most obvious copy number differences between different clones are highlighted by different colored boxes. The black square in chromosome 1 shows A1–A3 subclone-specific CNV. The black squares in chromosome 2 and chromosome 6 show additional CNVs of LM (E1–E3 clone-specific clones) compared with PT. E The diagram showing the tumor metastasis path of patient #1. CNVs, somatic mutations and mitochondrial mutations are also shown in the diagram. F Heatmap showing selected mitochondrial mutations of patient #1. All tumor cells have homozygous mutation at 13,368 point site. N: adjacent normal tissue. PT: primary tumor. LM: liver metastasis

Fig. 6

Integrated analyses of the associations between mitochondrial mutations and mitochondrial gene expression profiles. A CNV pattern and sampling strategy of patient #6. Adenomas have no obvious CNVs and different subclones with distinct CNV patterns are found in primary tumors. Lymph node metastasis has different CNV patterns with primary tumors. A: adenoma. PT: primary tumors. LyM: lymph node metastasis; LyM1, LyM2, and LyM3 represent three separated lymph node metastatic tumors. The pie charts reflect the proportion of different subclones in each region. Copy number gain and copy number loss were indicated with red and blue respectively. Different tissues were divided by bold black solid lines, and different clones within the same tissue block were divided by black dashed lines. The squares above the heatmap represent different chromosomes, black squares represent odd-numbered chromosomes and chromosome X, and light gray squares represent even-numbered chromosomes and chromosome Y. The squares on the left of the heatmap represent different subclones, and the same color represents the same clone. The most obvious copy number differences between different clones were highlighted by different colored boxes. B Heatmap showing selected mitochondrial mutations of patient #6. Blue, orange, red, and gray represent wild-type, heterozygous, and homozygous mutations and undefined (read depth < 9X) respectively. The triangle, circle and star represent the mutation state at the position 1,670-, 927-, and 8277-point site of mitochondria. C The diagram showing the tumor metastasis path of patient #6. The triangle, circle, and star represent the mutation state at the position 1,670-, 927-, and 8277-point site of mitochondria. Their colors represent the mutation status: blue, orange, red, and gray represent wild-type, heterozygous, and homozygous mutations respectively. D The boxplot showing the relative expression levels of MT-TV for patient #6 and patient #1. Other patients without MT-TV mutations have a similar pattern as patient #1 (data not shown). N: adjacent normal tissue. A: adenoma. PT: primary tumor. LyM: lymph node metastasis and LyM1-LyM3 represent three separated lymph node metastatic tumors. E The mutation state of mitochondrial gene MT-TV (site: 1,670) that is defined by whole exome sequencing data. Gray rectangle represents the wild-type, and the red rectangle represents the mutant type. The bar plot shows the allele frequency of mutations at this site. N: adjacent normal tissue. A: adenoma. PT: primary tumor. LyM: lymph node metastasis. LM: liver metastasis. For patients #10–#12, cells in the adjacent normal tissues were enriched by MACs or FACs; N represents EPCAM-positive cells and N_Neg represents EPCAM-negative cells. 1 cm, 3 cm, 10 cm, and 20 cm represent the distance to the edges of the tumor

Fig. 7

Integrated analyses of the associations between TP53 mutation and global gene expression profiles. A The tSNE map of patient #4, which has a TP53 frame-shift mutation. Cells were colored according to tissue origin, mitochondrial mutation, merge CNV types, and average expression levels (log₂(TPM+1)) of 306 primary tumor OLFM4⁺SOX9⁺ group specific genes, which are shown in B. N: adjacent normal tissue. PT: primary tumor. LyM: lymph node metastasis. LM: liver metastasis. Cells from primary tumor were clustered into two groups, PT-a and PT-b. Associated with Fig. S10A. B Heatmap shows the expression pattern of DEGs between OLFM4⁺ SOX9⁺, and OLFM4⁺ in all four cell groups of patient #4 C. The expression levels of IDH1, TPRKB (TP53RK-binding protein), and TP53RK (TP53 regulation kinase) were projected on tSNE maps. Colors from yellow to red represent expression levels from low to high. D PCA clustering of patient #9 epithelial cells. Cells were colored by TP53 Sanger results, CNV type, cell origin, and mitochondrial mutation (chrM:2897). N: adjacent normal tissue. PT: primary tumor. LyM: lymph node metastasis. LM: liver metastasis. E The expression levels of TP53, TPRKB, TP53PK, and IDH1 were projected to the PCA maps of patient #9 epithelial cells. TP53^High group means tumor epithelial cells that expressed TP53 gene according to the single-cell RNA sequencing data, while TP53^Low group represent tumor epithelial cells that did not expressed TP53. F The ratio of tumor epithelial cells that expressed TPRKB genes for each region. Metastatic tumor regions had more cells expressing TPRKB. ** represents p-value < 0.05. N: adjacent normal tissue. PT: primary tumor. LyM: lymph node metastasis. LM: liver metastasis. Neg represent TPRKB negative cells and Pos represent TPRKB positive cells. G The three dimensions plot shows the relationship between IDH1, TPRKB, TP53RK expression with TP53 mutations. The XY plane shows the PCA result of Patient #9. The three layers of IDH1, TPRKB, and TP53RK show the expression levels of these three genes on PCA map. The summary layer shows cells with TP53 mutations (according to Sanger sequencing). Red represents cells with TP53 mutation and expressed all three genes (IDH1, TPRKB and TP53RK) simultaneously, while yellow represent cells with TP53 mutation and expressed two of these genes. The TP53 Mutant layer shows the tissue origin of cells with TP53 mutations. Red represents cells collected from metastatic tumor and blue represents cells collected from primary tumor