Skip to main content

Transcriptomic intratumor heterogeneity of breast cancer patient-derived organoids may reflect the unique biological features of the tumor of origin

Abstract

Background

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.

Methods

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.

Results

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.

Conclusions

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.

Background

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 [3]. 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.

Methods

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 [8]. 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 [12].

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.

Data analysis

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 [13] 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 [14] 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”.

Cytokine measurement

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.

Immunohistochemistry

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).

Results

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%.

Fig. 1
figure 1

Process for establishing PDOs. A Schematics for collecting clinical specimens and establishing organoids. B Representative images showing the changes in the morphology of organoids over time. Arrow: densely expanding structures. Arrowhead: sheet-like formations around the spheroids. C Morphology of the 10 PDOs. The number after the P represents the number of passages of indicated PDO

Table 1 Clinicopathological characteristics of 10 patients

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).

Fig. 2
figure 2

Single-cell transcriptome analysis of 10 PDOs. A Schema of analytical flow. B Uniform manifold approximation and projection (UMAP) after clustering of scRNA-seq profiles from 10 PDOs. C Summary of the number of cells and clusters of each PDO. D Heatmap of scaled expression of cluster-specific gene expression signatures (ClustGSs)

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.

Fig. 3
figure 3

Meta-clusters designated by common expression signatures. A Heatmap of the Jaccard similarity index for 38 ClustGSs. Right annotations represent PDO samples and receptor status of the original tumors. B Annotation for meta-clusters and profiling of meta-ClustGSs

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).

Fig. 4
figure 4

Expression status of meta-ClustGSs in each PDO. A, C, E UMAP overlay representing the expression levels of meta-ClustGS 1 (cell cycle, G2/M phase) in (A), meta-ClustGS 3 (EMT1) in C, and meta-ClustGS 5 (estrogen response 1) in (E). B, D, F Violin plot representing expression levels of meta-ClustGS 1 (cell cycle, G2/M phase) in B, meta-ClustGS 3 (EMT1) in D, and meta-ClustGS 5 (estrogen response 1) in (F)

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.

Fig. 5
figure 5

PDO-specific ClustGSs. A Summary for nine PDO-specific ClustGSs. B Ridge plot of expression level of PDO207P#C4 ClustGS. Black arrows show a cell population with high expression of the ClustGS. C UMAP overlay representing expression levels of PDO207P#C4 ClustGS

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 [21]. 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.

Fig. 6
figure 6

Measurement of cytokines in PDO culture medium. A Cytokine arrays measured 102 cytokines in culture supernatant from the indicated PDOs. Positive controls at the upper left, upper right, and lower left edges and negative controls at the lower right edge. B Heatmap of the quantitative results of detected cytokines relative to the positive control. C UMAP overlay representing expression levels of the indicated genes

Discussion

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. [28] 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 [2], 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 [31]. 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.

Conclusions

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.

Abbreviations

ACTA2:

Actin alpha 2, smooth muscle

ANXA1:

Annexin A1

CCND1:

Cyclin D1

CDK1:

Cyclin-dependent kinase 1

CHI3L1:

Chitinase 3 like 1

ClustGS:

Cluster-specific gene signature

CST3:

Cystatin C

CYP4Z1:

Cytochrome P450 family 4 subfamily Z member 1

EIF2S3:

Eukaryotic translation initiation factor 2 subunit gamma

EIF3E:

Eukaryotic translation initiation factor 3 subunit E

EIF3H:

Eukaryotic translation initiation factor 3 subunit H

EMT:

Epithelial–mesenchymal transition

EPCAM:

Epithelial cell adhesion molecule

ER:

Estrogen receptor

FOXD1:

Forkhead box protein A1

GATA3:

GATA binding protein 3

HER2:

Human epidermal growth factor receptor 2

IBC:

Inflammatory breast cancer

ITH:

Intratumor heterogeneity

KRT5:

Keratin 5

KRT14:

Keratin 14

MKI67:

Marker of proliferation ki-67

MSigDB:

Molecular signatures database

MYC:

MYC proto-oncogene, bHLH transcription factor

PABPC1:

PolyA binding protein cytoplasmic 1

PDO:

Patient-derived organoid

PgR:

Progesterone receptor

scRNA-seq:

Single-cell RNA-seq

TF:

Transcription factors

VIM:

Vimentin

References

  1. McGranahan N, Swanton C. Clonal heterogeneity and tumor evolution: past, present, and the future. Cell. 2017;168(4):613–28.

    Article  CAS  PubMed  Google Scholar 

  2. Li Z, Seehawer M, Polyak K. Untangling the web of intratumour heterogeneity. Nat Cell Biol. 2022;24(8):1192–201.

    Article  CAS  PubMed  Google Scholar 

  3. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  4. 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.

    Article  CAS  PubMed  Google Scholar 

  5. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  6. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  7. Tuveson D, Clevers H. Cancer modeling meets human organoid technology. Science. 2019;364(6444):952–5.

    Article  CAS  PubMed  Google Scholar 

  8. 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.

    Article  CAS  PubMed  Google Scholar 

  9. 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.

    Article  CAS  PubMed  Google Scholar 

  10. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. 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.

    Article  CAS  Google Scholar 

  12. 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.

    Article  CAS  PubMed  Google Scholar 

  13. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  16. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  17. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  19. 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.

    CAS  PubMed  Google Scholar 

  20. 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.

    Article  CAS  PubMed  Google Scholar 

  21. 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.

    Article  CAS  PubMed Central  Google Scholar 

  22. 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.

    Article  CAS  PubMed  Google Scholar 

  23. 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.

    Article  CAS  PubMed  Google Scholar 

  24. Thomas G. An encore for ribosome biogenesis in the control of cell proliferation. Nat Cell Biol. 2000;2(5):E71-72.

    Article  CAS  PubMed  Google Scholar 

  25. Pelletier J, Thomas G, Volarevic S. Ribosome biogenesis in cancer: new players and therapeutic avenues. Nat Rev Cancer. 2018;18(1):51–63.

    Article  CAS  PubMed  Google Scholar 

  26. Breznik B, Mitrovic A, Lah TT, Kos J. Cystatins in cancer progression: more than just cathepsin inhibitors. Biochimie. 2019;166:233–50.

    Article  CAS  PubMed  Google Scholar 

  27. Cox TR. The matrix in cancer. Nat Rev Cancer. 2021;21(4):217–38.

    Article  CAS  PubMed  Google Scholar 

  28. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Suva ML, Tirosh I. Single-cell RNA Sequencing in cancer: lessons learned and emerging challenges. Mol Cell. 2019;75(1):7–12.

    Article  CAS  PubMed  Google Scholar 

  30. 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.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  31. 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.

    Article  PubMed  Google Scholar 

  32. Lempiainen H, Shore D. Growth control and ribosome biogenesis. Curr Opin Cell Biol. 2009;21(6):855–63.

    Article  PubMed  Google Scholar 

Download references

Acknowledgements

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.

Funding

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.].

Author information

Authors and Affiliations

Authors

Contributions

SS, YT, KO, and YO obtained clinical specimens. SS, KO, KY, and YS established organoids. SS and LY performed single-cell RNA-seq experiments. TO and MY performed immunohistochemical analysis. SS, KK, KO, and RM performed data analysis. KM, JS, MY, and RY supported organoid culture and experiments. TT, TS, TN, SO, RY, TU, and RM supervised projects. SS, KK, and RM wrote the manuscript. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Reo Maruyama.

Ethics declarations

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

Not applicable.

Competing interests

The authors declare no potential competing interest.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1

. 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.

Additional file 2

.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.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

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

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13058-023-01617-4

Keywords