Transcriptomic intratumor heterogeneity of breast cancer patient-derived organoids may reflect the unique biological features of the tumor of origin
Breast Cancer Research volume 25, Article number: 21 (2023)
The intratumor heterogeneity (ITH) of cancer cells plays an important role in breast cancer resistance and recurrence. To develop better therapeutic strategies, it is necessary to understand the molecular mechanisms underlying ITH and their functional significance. Patient-derived organoids (PDOs) have recently been utilized in cancer research. They can also be used to study ITH as cancer cell diversity is thought to be maintained within the organoid line. However, no reports investigated intratumor transcriptomic heterogeneity in organoids derived from patients with breast cancer. This study aimed to investigate transcriptomic ITH in breast cancer PDOs.
We established PDO lines from ten patients with breast cancer and performed single-cell transcriptomic analysis. First, we clustered cancer cells for each PDO using the Seurat package. Then, we defined and compared the cluster-specific gene signature (ClustGS) corresponding to each cell cluster in each PDO.
Cancer cells were clustered into 3–6 cell populations with distinct cellular states in each PDO line. We identified 38 clusters with ClustGS in 10 PDO lines and used Jaccard similarity index to compare the similarity of these signatures. We found that 29 signatures could be categorized into 7 shared meta-ClustGSs, such as those related to the cell cycle or epithelial–mesenchymal transition, and 9 signatures were unique to single PDO lines. These unique cell populations appeared to represent the characteristics of the original tumors derived from patients.
We confirmed the existence of transcriptomic ITH in breast cancer PDOs. Some cellular states were commonly observed in multiple PDOs, whereas others were specific to single PDO lines. The combination of these shared and unique cellular states formed the ITH of each PDO.
There is mounting evidence that cancer cell heterogeneity within tumors plays an essential role in part of cancer resistance and recurrence [1, 2]. In breast cancer, stronger intratumor heterogeneity (ITH) in estrogen receptor (ER) expression, as assessed by immunostaining, is associated with a poorer prognosis . ITH in human epidermal growth factor receptor 2 (HER2) amplification and expression also serves as a predictive indicator of poor therapeutic response [4, 5]. Although these studies focused on single genes essential to breast cancer, they were still suggestive of a relationship between cancer cell diversity and its malignant potential. Therefore, we must elucidate the molecular mechanisms underlying ITH and their functional significance to improve therapeutic strategies. However, partly a consequence of the limited number of suitable experimental models for the analysis of ITH, the overall research landscape of ITH in human cancers is still largely unexplored.
In recent years, three-dimensional organoid culture methods have been actively utilized in cancer research [6, 7] because patient-derived organoids (PDOs) are easier to manipulate and can be analyzed more quickly than patient-derived xenografts, yet still retain some of the tissue architecture and biological characteristics of the patients’ tumor of origin [8, 9]. Therefore, the intratumor heterogeneity of cancer cells may be preserved. Although several reports analyze PDOs from the perspective of intratumor heterogeneity in some cancer types, such as lung and colorectal cancer [10, 11], no reports examine the detailed cellular transcriptomic ITH in breast cancer PDOs.
In the present study, we use single-cell RNA-seq (scRNA-seq) to show that breast cancer PDOs are composed of several cell clusters with distinct cellular states. Furthermore, we define a cluster-specific gene signature (ClustGS) for each cell cluster in PDOs, and successfully capture one aspect of ITH. We found that seven ClustGSs were shared among multiple PDOs and some ClustGSs were specific to single PDOs only. The former represent cellular states related to the cell cycle, EMT, and estrogen response, whereas the latter represent cellular states often not clearly defined. These unique gene signatures specific to a certain PDO may reflect the biological characteristics of that PDO, indicating the assessment of cellular transcriptomic ITH may be helpful for characterizing PDOs.
Establishing organoids from clinical specimens
Surgical specimens were obtained by the core needle biopsy of surgically removed tumors. The specimens were dissociated into single cells using a MACS Tumor Dissociation Kit and a gentle MACS dissociator (Miltenyi Biotec, North Rhine-Westphalia, Germany) in accordance with the manufacturer’s instruction. For pleural effusion specimens, drainage fluid was collected after trocar placement, centrifuged, and purified with Red Blood Cell Lysis Solution lysis buffer (Miltenyl Biotec). The organoids were established from the obtained cells as described previously . Briefly, cell pellets were resuspended in BME gel (R&D Systems, Minnesota, USA) at 4 °C and seeded to form domes in 24-well or 48-well plates (IWAKI, Tokyo, Japan), polymerized at 37 °C for 10–20 min, and the prepared medium was added at a volume of 500 mL per 50 µL BME dome. The medium was changed every 3–4 days and the PDO was passaged depending on growth. Cultrex Organoid Harvesting Solution (Thermo Fisher Scientific, Massachusetts, USA) and 0.25 w/v% Trypsin-1 mM EDTA 4Na Solution with Phenol Red (FUJIFILM, Tokyo, Japan) were used for passaging of organoids. Expression of ER and PgR in surgical specimens was assessed using the Allred score .
Single-cell RNA-seq experiment
Organoids were dissociated into a single-cell suspension using harvesting solution and trypsin, and the cells were cryopreserved for the scRNA-seq experiment. Moreover, BC specimens were dissociated into single cells, as described above. The scRNA-seq libraries were prepared using the BD Rhapsody Single-Cell Analysis system (BD, New Jersey, USA) in accordance with the manufacturer's instruction. Briefly, after thawing, cells were labeled with BD Single-Cell Multiplexing Kit. The labeled cells were washed, pooled, and loaded onto the BD Rhapsody microwell cartridge; then cDNA was synthesized using the BD Rhapsody Whole Transcriptome Analysis Amplification Kit. The resulting gene expression and cellular label libraries were sequenced on Illumina NextSeq 550 platform (Illumina, California, USA) with paired-end reads (read1, 75 bp; index1, 8 bp; read2, 75 bp). Sequencing data were processed using the BD Rhapsody Analysis pipeline on the Seven Bridges Genomics platform and converted to a gene expression count matrix.
We utilized the web-based analysis pipeline for BD Rhapsody and BD Precise ASSAYS (https://www.sevenbridges.com/bdgenomics/) to generate the count matrixes and Seurat  for downstream analysis. RSEC counts were used as the input count matrixes. Low-quality cells with over 40% of mitochondrial RNAs and < 400 or > 9000 features were filtered out. To avoid misassignment between cells and samples, we performed pre-clustering of each PDO using Seurat’s standard flow and manually removed clusters that exhibited distant embedding from major clusters and different distribution of each sample tag that were possible contaminants from other samples.
For high-quality scRNA-seq data of each PDO, transcripts count matrixes were normalized to the total number of counts for the cell and multiplied by a scaling factor of 10,000. The normalized values were subsequently natural-log transformed using Seurat’s “NormalizeData()” function and a linear transformation was applied by “ScaleData()”. The principal component analysis was performed by “RunPCA()” with top 2,000 highest variable features identified by “FindVariableFeatures()” with vst selection method. We excluded PCs that contained ribosomal protein-encoding genes in more than half of the top 20 features. Then we performed Seurat’s standard clustering procedures using “FindNeighbors()” and “FindClusters()” with the top 20 PCs and a resolution of 0.4. To visualize data, “RunUMAP()” was used with the same PCs to identify the clusters.
After clustering of each PDO, we identified differentially expressed features as cluster-specific gene expression signatures (ClustGS) by using the “FindAllMarkers()” function with “only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0”. We removed mitochondrial RNAs or ribosomal protein-encoding genes and selected the genes with an adjusted P-value of less than 0.05. The top 25 most significant genes in each cluster of each PDO were selected as ClustGS. Any clusters with fewer than 25 significant genes (adjusted P-value < 0.05) were not considered for further analysis. Meta-clusters were manually defined by the Jaccard similarity index. All genes included in at least two ClustGS were defined as meta-cluster specific genes (meta-ClustGS).
Gene enrichment analysis for each ClustGS was performed using the “enrichr()” function of ClusterProfiler v3.18.1  with MsigDB Hallmark genes obtained through msigdbr’s msigdbr(species = “Homo sapiens,” category = “H”) function. The false discovery rate (FDR) was used as the adjusted P-value.
Gene overlap analysis between ClustGS and cluster-specific expressed genes of corresponding primary tumors or different passage organoids was performed using “GeneOverlap” package’s newGeneOverlap() and testGeneOverlap() functions. To obtain cluster-specific expressed genes in each dataset, we selected the genes with < 0.05 adjusted P-values calculated by Seurat’s FindAllMarkers() with options of “only.pos = TRUE, min.pct = 0.25, logfc.threshold = 0”.
The extracellular cytokines in the culture medium were comprehensively semi-quantified using the Proteome Profiler Array Human XL Cytokine Array Kit (R&D Systems) in accordance with the manufacturer’s instruction. Briefly, culture medium from four organoids (PDOs 155, 165, 166, 180) was collected 72 h after the change of fresh medium and frozen at − 20 °C. Between 59.5 and 846.2 µL of the culture medium, normalized by RNA quantity from the PDOs in the same culture, was loaded into the cytokine array. The Odyssey XF Imaging System (LI-COR, Nebraska, USA) was used to detect fluorochromes; the absorbance was set to 685 nm and the exposure time was set to 10 min. The signal intensity corresponding to each factor in the array was quantified by ImageJ software, and the expression level was calculated relative to the positive control.
Tumor sections were stained with hematoxylin–eosin (HE) and immunohistochemically examined for CHI3L1 and CST3 using a Nichirei Histofine system (Nichirei Biosciences Inc.). Tissue sections were incubated with primary anti-CHI3L1 antibody (Cell Signaling, #47,066) and anti-CST3 antibody (Abcam, ab109508) at a 1:800 and 1:4000 dilution, respectively, and detected with HRP-labeled polymer-conjugated secondary antibody (Histofine Simple Stain MAX PO, multi, Nichirei). Color development was achieved with 3,3′-diaminobenzidine tetrahydrochloride. Organoids were fixed and paraffin-embedded using Epredia HistoGel Specimen Processing Gel (Thermo Scientific) and Tissue-Tec VIP (SAKURA, Japan). Immunohistochemical staining on the organoid sections was performed by incubating with anti-CHI3L1 antibody (Cell Signaling, #47,066) and anti-CST3 antibody (Abcam, ab109508) at a 1:400 and 1:500 dilution, respectively. Signals were detected using Dako EnVision + Dual Link System-HRP (Agilent).
Organoids established from breast cancer tissues contained morphologically heterogeneous cell populations
Ten PDO lines established during a certain period of time were subjected to the analysis in this study (Fig. 1A). The clinicopathological features of the 10 cases are presented in Table 1. Briefly, eight were from patients with primary breast cancer, one case (PDO210) was local recurrence in the chest wall after total mastectomy, and one case (PDO207P) was lung metastasis with pleural effusion. In addition, one of the eight cases of primary breast cancer (PDO155) developed metastasis soon after surgery. Tumors with high proliferative potential may have been more frequently established as organoids; there were five cases with nuclear grade 3 and Ki67 > 20%.
The speed of organoid growth differed for each PDO. PDO155, 166, and 207P had the fastest growth rate, with a doubling time of approximately 7 days, which continued even after ten passages. PDO 165 and 180 grew more slowly than the three PDO lines listed above, with a doubling time of approximately 14 days. The other lines grew more slowly.
Organoid morphology was also heterogeneous among PDOs (Fig. 1B, C). Spheroids were common in all PDOs, but some PDOs had mixed populations of organoids with different shapes. For example, PDO155 consisted of two different cell populations: spheroids and densely expanding structures. PDO154, 165, and 180 appeared to form sheet-like structures around the spheroids that dragged each spheroid along, altering their 3D position every few days. PDO207P contained grapevine cluster-like formations. PDO166, 195, 202, 203, and 210 formed homogenous spheroids.
Collectively, our results present the inter- and intratumor heterogeneity for organoid growth and morphology, suggesting that the PDOs reflect the biological properties of the original tumors.
Examining ITH of PDOs by single-cell transcriptome profiling
We performed single-cell transcriptome analysis of 10 PDOs by multiplexing scRNA-seq using the BD Rhapsody Single-Cell Analysis System. After assignment and quality filtering, we obtained the high-quality expression profiles of 9548 cells: an average of 6635 counts and 2203 genes per cell (Additional file 1: Figure S1; “Method” Section). For simple quality assessment, we used the Seurat package to cluster the cells and examine the expression of marker genes (Additional file 1: Figure S1). Almost all cells expressed epithelial markers such as EPCAM and keratins, and we considered that the PDO lines used in this study contained nearly no stromal cells. Therefore, we hereby note that the subsequent analysis assesses the ITH of cancer cells.
The cancer cells were clustered as per patient samples in scRNA-seq analysis (Additional file 1: Figure S1), which is consistent with previous reports [15,16,17,18]. To investigate the ITH in each PDO, we developed the analytical approach as follows: (1) cluster the cells in each PDO, (2) identify cluster-specific expressed genes, and (3) calculate their similarity and classify clusters (Fig. 2A). Through a standard analytical flow in Seurat, the cells were classified into 3–6 clusters per each PDO (3 clusters in PDO180, 195, and 202; 6 clusters in PDO210) (Fig. 2B, C). Thus, we identified a set of genes expressed specifically in each cluster within each PDO and defined the most significant 25 differentially expressed genes as the cluster-specific gene signature (ClustGS). Any clusters with fewer than 25 significant genes (adjusted P-value < 0.05) were not considered for further analysis. Each PDO had at least one cluster with ClustGS, and in total, 38 clusters with a ClustGS were identified (Fig. 2C).
For example, PDO165 had four clusters with a ClustGS (Fig. 2D). ClustGSs in PDO165#C1 contained luminal-lineage transcription factors (TFs) (FOXA1, GATA3); PDO165#C2 contained TGFB2 and DKK1; PDO165#C3 contained cell cycle-related genes (CDK1, MKI67); and PDO165#C4 contained basal keratins (KRT5, KRT14) (Fig. 2D; Additional file 2: Table S1). Consistently, the molecular signatures database (MSigDB) hallmark gene enrichment analysis suggested that PDO165#C1 was related to estrogen response; PDO165#C2 and PDO165#C4 were related to the epithelial–mesenchymal transition (EMT); and PDO165#C3 was related to the cell cycle (Additional file 1: Figure S2). These results suggest that each PDO is composed of multiple subpopulations with distinct cell states and lineages with the cell cycle, estrogen response, and EMT-like gene expression programs.
We next assessed whether the ITH observed in the PDOs reflected the ITH in the original tumor tissues and whether it was maintained after repeated PDO passages. We performed scRNA-seq experiments on two additional original breast cancer tissues (165 and 210) and two different passages of PDO lines (180 and 210). Subsequently, we performed the same analysis as above, and confirmed the existence of ITH in each sample (Additional file 1: Figure S3; Additional file 2: Table S2). Comparison of the characteristics of the identified cell populations by looking at the overlap of cluster-specific expressed genes revealed that they were composed of somewhat similar cell populations between original tissues and PDOs as well as between different passages of the same PDO lines (Additional file 1: Figure S4; Additional file 2: Table S3), although not completely one-to-one correspondence. Thus, we decided to analyze these PDOs further, considering that they retain some extent of ITH of the original tissues even after successive passages.
Identifying cell populations defined by common expression programs across PDOs
To compare the degree of ITH between PDOs, we calculated the Jaccard similarity index across 38 ClustGS from each cluster. Based on the similarity of ClustGS, seven meta-clusters were identified with the hierarchical clustering guide, representing a set of cell clusters characterized by the same expression programs across PDOs (Fig. 3A). To annotate each of the meta-clusters, we selected the overlapping genes across ClustGS in at least two PDOs as a meta-cluster (“meta-ClustGS”; Additional file 2: Table S4). Gene enrichment analysis for meta-ClustGSs using the MsigDB Hallmark gene set suggested that meta-ClustGSs 1–2 were associated with cell cycle, notably the G2/M phase for meta-ClustGS 1 and G1/S phase for meta-ClustGS 2; meta-ClustGSs 3, 4, and 6 were associated with the EMT; and meta-ClustGSs 5 and 7 were associated with estrogen response (Fig. 3B; Additional file 1: Figure S5). Although meta-ClustGSs 3, 4, and 6 all contained EMT-associated genes, each gene signature was different: basal cytokeratin (such as KRT5, KRT14) in meta-ClustGS 3, mesenchymal marker ACTA2 in meta-ClustGS 4, and the metastasis-associated gene ANXA1 in meta-ClustGS 6. There were also two meta-ClustGSs associated with estrogen response: luminal-lineage specific TFs FOXA1 and GATA3 in meta-ClustGS 5, CCND1 in meta-ClustGS 7. These results suggested that each PDO contained several cell populations that were regulated by the defined gene expression programs that were common across PDOs.
Common gene expression programs highlighted ITH across PDOs
We subsequently characterized each PDO using meta-ClustGSs. The meta-ClustGS 1 (Cell cycle—G2/M) was present in 7 of 10 PDO lines (Figs. 3A, 4A, B; Additional file 1: Figure S6), indicating that the proliferating cell populations were common in PDOs. Interestingly, both meta-ClustGS 3 (EMT1) containing basal markers (KRT5 and KRT14) and meta-ClustGS 5 (estrogen response 1) containing luminal-lineage TFs (FOXA1 and GATA3) were observed in four PDO lines (PDO155, 165, 166, and 180) derived from patients with ER+ /HER2– or ER–/HER2+ tumors, suggesting that these PDOs had distinct ITH of co-existing of basal-like and luminal cells (Fig. 4C–F; Additional file 1: Figure S6). Clusters of cells with meta-ClustGS 5 expression were observed in PDO155 derived from ER–/HER2+ tumors, suggesting that a small number of cells expressed ER target genes that were not detected by the bulk assessment, i.e., immunohistochemistry, but were observable at a single-cell resolution (Fig. 4E; Additional file 1: Figure S6). Most cells in three PDOs (PDO180, 195, and 202) derived from ER+ tumors expressed the meta-ClustGS 3 (EMT1, basal-like), implying that the property of ER+ cancer cells was altered in organoid culture or that part of organoids was derived from juxta-tumoral basal epithelial cells in these PDOs (Fig. 4C, D).
Another estrogen response-related meta-cluster, meta-ClustGS 7, was observed in PDO154, 207P, and 210. Interestingly, they had distinct clinical features such as ER + /HER2 + (PDO154), pleural effusion sample (PDO207P), and early recurrence (PDO210). These suggested that meta-ClustGS 7 may reflect a different cellular state to normal ER-positive breast cancer cells with a good prognosis, even though it is annotated with the same estrogen response by MsigDB analysis. For example, meta-ClustGS 7 contained CCND1, which has been implicated in endocrine resistance and poor outcome [19, 20], suggesting that this gene signature might reflect the peculiarity of a cell population (Additional file 1: Figure S6).
These results indicate that meta-ClustGSs can describe each PDO with a qualitative and quantitative assessment of ITH.
PDO-specific cell clusters reflected distinct properties for PDOs
In the process of identifying meta-clusters, we also identified nine PDO-specific clusters (Fig. 5A; Additional file 2: Table S5). These ClustGSs corresponded to four PDO lines, each derived from cases with characteristic clinical courses or histopathology: PDO155 from inflammatory breast cancer (IBC), PDO203 from invasive micropapillary carcinoma, PDO207P from metastatic breast cancer, and PDO210 from recurrent breast cancer. Thus, ClustGSs not shared with any other PDOs possibly reflected the clinical or molecular features for each PDO.
For example, PDO207P#C4 ClustGS included the mesenchymal marker VIM and CYP4Z1, which are associated with poor prognosis and high-grade tumors in breast cancer . This ClustGS was highly expressed in a certain population of PDO207P (Fig. 5B, C). As PDO207P was derived from the metastatic pleural effusion, this cluster probably reflected the distinct features of mesenchymal and aggressive disease in PDO207P.
Organoids from a patient with inflammatory breast cancer
We further investigated PDO155 which was derived from a patient with IBC. The patient developed liver metastases early after surgery and had a poor clinical outcome [22, 23]. To our knowledge, no reports are available on PDO lines being established from patients with IBC and then precisely analyzed. We therefore considered that it might be possible to find the biological features of IBC in the PDO155. We focused on a specific cell population, PDO155P#C1 (Fig. 2A), with a ClustGS that was not similar to other ClustGSs. In the gene signature, genes associated with RNA transport and translation initiation functions were highly enriched (EIF3H, PABPC1, EIF3E, and EIF2S3; Additional file 2: Table S5). In this cell population, translation was expected to occur preferentially to cell cycle progression (Fig. 4C), suggesting that cell size may increase [24, 25]. This may be related to the morphological heterogeneity observed specifically in PDO155, in which some of the organoids formed gigantic spheroids (Fig. 1B, C).
The gene signature also included several genes encoding cytokines, such as CST3 and CHI3L1. We thought that the specific cytokines released from cancer cells might be related to the pathophysiological characteristics of IBC. Thus, we examined cytokines in the culture medium in PDO155 and other PDOs using a cytokine array (Fig. 6A). CST3 and CHI3L1 were highly detected in PDO155; in particular, the secretion of CST3 was specifically observed in PDO155. CST3 is known to have various functions on cancer progression [26, 27], which may be related to the feature that part of PDO155 forms huge clumps. Interestingly, these genes were not uniformly expressed in all cells; they were heterogeneously and non-randomly expressed (Fig. 6B). By performing immunohistochemistry of PDO155 and pathology specimens of this patient, we confirmed that CST3 and CHI3L1 proteins were heterogeneously expressed in PDO155 and even in the original tumor tissue (Additional file 1: Figure S7). This indicates that only a subset of the cell population may only exhibit the specific characteristics observed in the bulk assay; conversely, a unique feature visible in the ITH may display the representative feature of that PDO.
This study aimed to investigate the intratumor transcriptomic heterogeneity in breast cancer PDOs. We performed scRNA-seq analysis and found that PDOs were composed of 3–6 cell populations with distinct cellular states. Furthermore, we characterized the cellular states of these cell populations based on the top 25 genes expressed specifically in each population, which we defined as the cluster-specific gene signature (ClustGS). We then compared the cellular states within and between PDOs by examining the similarity of the ClustGSs. Seven meta-ClustGSs, reflecting cellular states, were common among multiple PDOs, and nine ClustGSs were found to be unique to individual organoids. The common meta-ClustGSs included cell states related to the cell cycle, EMT, and estrogen response. In contrast, the cellular states unique to the individual PDOs appeared to display unique characteristics for that specific PDOs. These results indicated that we may be able to estimate the unique characteristics of a particular PDO by examining the ITH within the PDO; a comparison with other PDOs may not be necessary.
Using single-cell RNA-seq technology, we confirmed the existence of ITH in breast cancer PDOs. Several studies have reported the ITH of cancer cells and its characteristics in cancer tissues and cell lines [15, 16, 18, 28]. Kinker et al.  reported 12 recurrent heterogeneous programs, including cell cycle, EMT, stress, and senescence. They also demonstrated that many transcriptomic ITH patterns of cancer cells reflected the intrinsic plasticity of cancer cells, even outside of the native microenvironment. Our data show ITH related to cell cycle and EMT within breast cancer PDOs, which is consistent with the above study. In addition, several PDOs also displayed ITH related to estrogen response and basal cells, which may be specific to the lineage plasticity of breast cancer.
We also observed patient-specific ITH patterns in several PDOs. Previous studies tended to focus on the recurrent patterns of ITH and ignore unique ITH [28, 29]. However, we hypothesized that the patient-specific ITH might represent the unique characteristics of the PDO and contain useful information. As the ITH pattern in cancer cells is thought to reflect the cell-of-origin and the activity of the lineage specific transcription factors , the ITH pattern should provide some information on the origin or pathophysiology of the cancer cell population.
From this perspective, we focused on PDO155 derived from a patient with aggressive inflammatory breast cancer (IBC). IBC is a cancer subtype with a poor prognosis characterized by extensive lymphatic embolism and cutaneous inflammation, and its etiology and molecular mechanisms are poorly understood [23, 30]. PDO155 was distinct from the other nine organoids in terms of morphology. It grew quickly, and some of the organoids became huge (Fig. 1B, C). These unique characteristics of morphology may possibly reflect the nature of IBCs, which tend to form emboli in the lymphatic vessels. The scRNA-seq analysis also identified cell populations in a particular state that were not observed in other organoids, which may hold some clues. The signature (PDO155P#C1 ClustGS) included several genes related to translation. PDO155 had MYC amplification (data not shown) and may have enhanced ribosome biogenesis . The abnormal increase in translation may be related to the increase in cell size and the giant organoids [24, 32]. Although further functional verification is needed, these data indicate that the specific ITH pattern may predict and explain the characteristics of the original IBC.
We described the landscape of transcriptomic ITH within breast cancer PDOs. We observed cell populations with cellular states common to various PDOs and cell states specific to individual PDOs. The latter could reflect the unique biology of the tumor of origin.
Availability of data and materials
R code for performing the analyses in available at (https://github.com/KoheiKumegawa/BC_10PDOs_scRNA). Processed scRNA-seq data are available from the lead contact upon request.
Actin alpha 2, smooth muscle
Cyclin-dependent kinase 1
Chitinase 3 like 1
Cluster-specific gene signature
Cytochrome P450 family 4 subfamily Z member 1
Eukaryotic translation initiation factor 2 subunit gamma
Eukaryotic translation initiation factor 3 subunit E
Eukaryotic translation initiation factor 3 subunit H
Epithelial cell adhesion molecule
Forkhead box protein A1
GATA binding protein 3
Human epidermal growth factor receptor 2
Inflammatory breast cancer
Marker of proliferation ki-67
Molecular signatures database
MYC proto-oncogene, bHLH transcription factor
PolyA binding protein cytoplasmic 1
McGranahan N, Swanton C. Clonal heterogeneity and tumor evolution: past, present, and the future. Cell. 2017;168(4):613–28.
Li Z, Seehawer M, Polyak K. Untangling the web of intratumour heterogeneity. Nat Cell Biol. 2022;24(8):1192–201.
Lindström LS, Yau C, Czene K, Thompson CK, Hoadley KA, Van’t Veer LJ, Balassanian R, Bishop JW, Carpenter PM, Chen YY, et al. Intratumor heterogeneity of the estrogen receptor and the long-term risk of fatal breast cancer. J Natl Cancer Inst. 2018;110(7):726–33.
Horii R, Nitta H, Nojima M, Maruyama R, Ueno T, Ito Y, Ohno S, Banks P, Kanda H, Akiyama F. Predictive significance of HER2 intratumoral heterogeneity, determined by simultaneous gene and protein analysis, for resistance to trastuzumab-based treatments for HER2-positive breast cancer. Virchows Arch. 2021;479(1):13–21.
Filho OM, Viale G, Stein S, Trippa L, Yardley DA, Mayer IA, Abramson VG, Arteaga CL, Spring LM, Waks AG, et al. Impact of HER2 heterogeneity on treatment response of early-stage HER2-positive breast cancer: phase II neoadjuvant clinical trial of T-DM1 combined with pertuzumab. Cancer Discov. 2021;11(10):2474–87.
Kim J, Koo BK, Knoblich JA. Human organoids: model systems for human biology and medicine. Nat Rev Mol Cell Biol. 2020;21(10):571–84.
Tuveson D, Clevers H. Cancer modeling meets human organoid technology. Science. 2019;364(6444):952–5.
Sachs N, de Ligt J, Kopper O, Gogola E, Bounova G, Weeber F, Balgobind AV, Wind K, Gracanin A, Begthel H, et al. A living biobank of breast cancer organoids captures disease heterogeneity. Cell. 2018;172(1–2):373-386e310.
Yan HHN, Siu HC, Law S, Ho SL, Yue SSK, Tsui WY, Chan D, Chan AS, Ma S, Lam KO, et al. A comprehensive human gastric cancer organoid biobank captures tumor subtype heterogeneity and enables therapeutic screening. Cell Stem Cell. 2018;23(6):882–897811.
Monberg ME, Geiger H, Lee JJ, Sharma R, Semaan A, Bernard V, Wong J, Wang F, Liang S, Swartzlander DB, et al. Occult polyclonality of preclinical pancreatic cancer models drives in vitro evolution. Nat Commun. 2022;13(1):3652.
Okamoto T, duVerle D, Yaginuma K, Natsume Y, Yamanaka H, Kusama D, Fukuda M, Yamamoto M, Perraudeau F, Srivastava U, et al. Comparative analysis of patient-matched PDOs revealed a reduction in OLFM4-associated clusters in metastatic lesions in colorectal cancer. Stem Cell Rep. 2021;16(4):954–67.
Harvey JM, Clark GM, Osborne CK, Allred DC. Estrogen receptor status by immunohistochemistry is superior to the ligand-binding assay for predicting response to adjuvant endocrine therapy in breast cancer. J Clin Oncol. 1999;17(5):1474–81.
Hao Y, Hao S, Andersen-Nissen E, Mauck WM 3rd, Zheng S, Butler A, Lee MJ, Wilk AJ, Darby C, Zager M, et al. Integrated analysis of multimodal single-cell data. Cell. 2021;184(13):3573-3587e3529.
Yu G, Wang LG, Han Y, He QY. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–7.
Tirosh I, Venteicher AS, Hebert C, Escalante LE, Patel AP, Yizhak K, Fisher JM, Rodman C, Mount C, Filbin MG, et al. Single-cell RNA-seq supports a developmental hierarchy in human oligodendroglioma. Nature. 2016;539(7628):309–13.
Venteicher AS, Tirosh I, Hebert C, Yizhak K, Neftel C, Filbin MG, Hovestadt V, Escalante LE, Shaw ML, Rodman C, et al. Decoupling genetics, lineages, and microenvironment in IDH-mutant gliomas by single-cell RNA-seq. Science. 2017;355(6332):eaai8478.
Filbin MG, Tirosh I, Hovestadt V, Shaw ML, Escalante LE, Mathewson ND, Neftel C, Frank N, Pelton K, Hebert CM, et al. Developmental and oncogenic programs in H3K27M gliomas dissected by single-cell RNA-seq. Science. 2018;360(6386):331–5.
Tirosh I, Izar B, Prakadan SM, Wadsworth MH 2nd, Treacy D, Trombetta JJ, Rotem A, Rodman C, Lian C, Murphy G, et al. Dissecting the multicellular ecosystem of metastatic melanoma by single-cell RNA-seq. Science. 2016;352(6282):189–96.
Roy PG, Pratt N, Purdie CA, Baker L, Ashfield A, Quinlan P, Thompson AM. High CCND1 amplification identifies a group of poor prognosis women with estrogen receptor positive breast cancer. Int J Cancer. 2010;127(2):355–60.
Bostner J, Ahnstrom Waltersson M, Fornander T, Skoog L, Nordenskjold B, Stal O. Amplification of CCND1 and PAK1 as predictors of recurrence and tamoxifen resistance in postmenopausal breast cancer. Oncogene. 2007;26(49):6997–7005.
Wu W, Warner M, Wang L, He WW, Zhao R, Guan X, Botero C, Huang B, Ion C, Coombes C, et al. Drivers and suppressors of triple-negative breast cancer. Proc Natl Acad Sci U S A. 2021;118(33):e2104162118.
Faldoni FLC, Villacis RAR, Canto LM, Fonseca-Alves CE, Cury SS, Larsen SJ, Aagaard MM, Souza CP, Scapulatempo-Neto C, Osorio C, et al. Inflammatory breast cancer: clinical implications of genomic alterations and mutational profiling. Cancers (Basel). 2020;12(10):2816.
Lim B, Woodward WA, Wang X, Reuben JM, Ueno NT. Inflammatory breast cancer biology: the tumour microenvironment is key. Nat Rev Cancer. 2018;18(8):485–99.
Thomas G. An encore for ribosome biogenesis in the control of cell proliferation. Nat Cell Biol. 2000;2(5):E71-72.
Pelletier J, Thomas G, Volarevic S. Ribosome biogenesis in cancer: new players and therapeutic avenues. Nat Rev Cancer. 2018;18(1):51–63.
Breznik B, Mitrovic A, Lah TT, Kos J. Cystatins in cancer progression: more than just cathepsin inhibitors. Biochimie. 2019;166:233–50.
Cox TR. The matrix in cancer. Nat Rev Cancer. 2021;21(4):217–38.
Kinker GS, Greenwald AC, Tal R, Orlova Z, Cuoco MS, McFarland JM, Warren A, Rodman C, Roth JA, Bender SA, et al. Pan-cancer single-cell RNA-seq identifies recurring programs of cellular heterogeneity. Nat Genet. 2020;52(11):1208–18.
Suva ML, Tirosh I. Single-cell RNA Sequencing in cancer: lessons learned and emerging challenges. Mol Cell. 2019;75(1):7–12.
Chakraborty P, George JT, Woodward WA, Levine H, Jolly MK. Gene expression profiles of inflammatory breast cancer reveal high heterogeneity across the epithelial-hybrid-mesenchymal spectrum. Transl Oncol. 2021;14(4):101026.
van Riggelen J, Yetil A, Felsher DW. MYC as a regulator of ribosome biogenesis and protein synthesis. Nat Rev Cancer. 2010;10(4):301–9.
Lempiainen H, Shore D. Growth control and ribosome biogenesis. Curr Opin Cell Biol. 2009;21(6):855–63.
The authors are grateful to Ms. Satoko Baba (Division of Pathology, Cancer Institute, Japanese Foundation for Cancer Research) for helping our pathological assessment. The authors would like to thank Enago for the English language review.
This work was supported in part by JSPS KAKENHI grant numbers JP20K07708 [to R.M.], JP17K18336 [to R.M.], the Research Program on Hepatitis from Japan Agency for Medical Research and Development, AMED (JP21fk0210054, JP21fk0210092) [to R.M.], the Vehicle Racing Commemorative Foundation [to R.M., T.U., L.Y.], and the Novartis Research Grant [to R.M.].
Ethics approval and consent to participate
The participant gave written informed consent before the collection of specimens. The protocol was approved by the institutional ethical committee of Cancer Institute Hospital, Japanese Foundation for Cancer Research (No. 2018-1168).
Consent for publication
The authors declare no potential competing interest.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
. Figure S1: A and B. QC metrics showing unique feature counts and percentage of mitochondrial RNA per cell. C. UMAP visualization of scRNA-seq profiles from 10 PDOs. Each dot represents a single-cell colored by its corresponding PDOs. D and E. UMAP overlay representing expression levels of EPCAM and KRT8. Figure S2: Hallmark terms and adjusted p-values detected by MSigDB hallmark gene enrichment analysis for each ClustGS in PDO165. Figure S1: A and C. UMAP visualization of 1,822 and 6,362 cells annotated as EPCAM-positive cells within scRNA-seq profiles of primary breast cancer tissues of patients 165 (A) and 210 (C). B and D. Heatmaps showing the number (top) and significance (bottom) of overlapping genes between each ClustGS identified in indicated PDOs and each cluster-specific expressed gene list identified in its original tumor tissues. Figure S4: A and B. UMAP visualization of scRNA-seq profiles of PDO180 P6 (A) and PDO210 P8 (B). C and D. Heatmaps showing the number (top) and significance (bottom) of overlapping genes between each ClustGS identified in indicated PDOs and each cluster-specific expressed gene list identified in PDOs. Figure S5: Hallmark terms and adjusted p-values detected by MSigDB hallmark gene enrichment analysis for 7 meta-ClustGSs. Figure S6: Heatmap showing which meta-ClustGS are present or absent in each PDO. Figure S7: HE staining and immunohistochemistry (CHI3L1 and CST3) of PDO155 and its original tumor tissue.
.Table S1: 38 ClustGSs across 10 PDOs.Table S2: Cluster-specific expressed gene lists in P165 and P210.Table S3: Cluster-specific expressed gene lists in PDO180 P6 and PDO210 P8.Table S4: Meta-ClustGSs.Table S5 PDO-specific ClustGSs.
About this article
Cite this article
Saeki, S., Kumegawa, K., Takahashi, Y. et al. Transcriptomic intratumor heterogeneity of breast cancer patient-derived organoids may reflect the unique biological features of the tumor of origin. Breast Cancer Res 25, 21 (2023). https://doi.org/10.1186/s13058-023-01617-4
- Intratumor heterogeneity
- Breast cancer
- Patient-derived organoids
- Cancer cell diversity
- Inflammatory breast cancer