Martín-Pardillos et al. BMC Cancer (2019) 19:666 https://doi.org/10.1186/s12885-019-5883-y RESEARCH ARTICLE Open Access The role of clonal communication and heterogeneity in breast cancer Ana Martín-Pardillos1,2* , Ángeles Valls Chiva1†, Gemma Bande Vargas1†, Pablo Hurtado Blanco3, Roberto Piñeiro Cid2,3, Pedro J. Guijarro1, Stefan Hümmer1,2, Eva Bejar Serrano1,2, Aitor Rodriguez-Casanova4, Ángel Diaz-Lagares2,4, Josep Castellvi5, Samuel Miravet-Verde6, Luis Serrano6,7,8, María Lluch-Senar6, Víctor Sebastian9,10, Ana Bribian11, Laura López-Mascaraque11, Rafael López-López4,12 and Santiago Ramón y Cajal1,2,5* Abstract Background: Cancer is a rapidly evolving, multifactorial disease that accumulates numerous genetic and epigenetic alterations. This results in molecular and phenotypic heterogeneity within the tumor, the complexity of which is further amplified through specific interactions between cancer cells. We aimed to dissect the molecular mechanisms underlying the cooperation between different clones. Methods: We produced clonal cell lines derived from the MDA-MB-231 breast cancer cell line, using the UbC-StarTrack system, which allowed tracking of multiple clones by color: GFP C3, mKO E10 and Sapphire D7. Characterization of these clones was performed by growth rate, cell metabolic activity, wound healing, invasion assays and genetic and epigenetic arrays. Tumorigenicity was tested by orthotopic and intravenous injections. Clonal cooperation was evaluated by medium complementation, co-culture and co-injection assays. Results: Characterization of these clones in vitro revealed clear genetic and epigenetic differences that affected growth rate, cell metabolic activity, morphology and cytokine expression among cell lines. In vivo, all clonal cell lines were able to form tumors; however, injection of an equal mix of the different clones led to tumors with very few mKO E10 cells. Additionally, the mKO E10 clonal cell line showed a significant inability to form lung metastases. These results confirm that even in stable cell lines heterogeneity is present. In vitro, the complementation of growth medium with medium or exosomes from parental or clonal cell lines increased the growth rate of the other clones. Complementation assays, co-growth and co-injection of mKO E10 and GFP C3 clonal cell lines increased the efficiency of invasion and migration. Conclusions: These findings support a model where interplay between clones confers aggressiveness, and which may allow identification of the factors involved in cellular communication that could play a role in clonal cooperation and thus represent new targets for preventing tumor progression. Keywords: Tumor, Breast, Cancer, Metastasis, Heterogeneity, Clone, Communication, Cooperation, MDA-MB-231 Background 20% [1, 2]. The prognosis is even worse in the case of Breast cancer is the most common cancer in women triple negative breast cancer (TNBC), which lacks the worldwide and the second leading cause of cancer-related major targets of approved therapies, limiting treatment death. In breast cancer, metastatic disease progression car- options to surgery, radiotherapy and/or systemic chemo- ries a poor prognosis, with 5-year survival rates of below therapy [3]. Breast cancer, like many other malignant tumors, shows a great molecular and phenotypic heterogeneity, both at * Correspondence: ana.martin@vhir.org; sramon@vhebron.net † Ángeles Valls Chiva and Gemma Bande Vargas contributed equally to this an inter- and intra-tumoral level [4–6]. This heterogeneity work. may be understood to be due to the accumulation of mo- 1 Translational Molecular Pathology Group, Vall d’Hebron Research Institute, lecular alterations from an initial clone that undergoes Barcelona, Spain Full list of author information is available at the end of the article Darwinian selection. During this evolution, certain core © The Author(s). 2019 Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. 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. Martín-Pardillos et al. BMC Cancer (2019) 19:666 Page 2 of 26 molecular alterations will be largely governed by cell- factors allow synergistic growth and ultimately contrib- extrinsic phenomena, such as environmental conditions, ute to tumor progression and metastasis. including the immune response [7–12]. This mechanism of tumor evolution has been observed in multiple tumor Methods types [13–15]. Tumor heterogeneity can result in an in- Cell culture and reagents complete or incorrect diagnosis when small tumor sam- The study was conducted using the MDA-MB-231 cell line. ples or biopsies are used; consequently, treatments may be This is an epithelial human breast cancer cell line, estab- directed against targets that are not expressed throughout lished from a pleural effusion of a 51-year-old Caucasian the entire tumor [16–18]. woman with metastatic breast adenocarcinoma. The MDA- Furthermore, the genetic-based type of tumor evolu- MB-231 breast cancer cell line was purchased from the tion cannot fully explain how all the functional and American Type Culture Collection (ATCC) and authenti- phenotypic alterations required to fulfill Weinberg’s cated by DNA profiling using short tandem repeat (STR) 10 hallmarks of cancer [19] can accumulate within a (GenePrint® 10 System, Promega) at Genomics Core Facil- single cell clone. For this reason, the theory of clonal ity, Instituto de Investigaciones Biomédicas “Alberto Sols” cooperation asserts that tumor clones have comple- CSIC-UAM. Viral, bacterial and parasitic pathogen analysis mentary genetic alterations that synergistically con- by RT-PCR/PCR was performed at Dynamimed Research tribute to tumor progression and metastasis [20–27]. Company: no genetic material was detected. Cells were This clonal cooperation, together with the influence maintained at 37 °C in a 5% CO2 humidified incubator of the tumor stroma and the immune system, is more (AutoFlow UN-5510, Nuaire), in Dulbecco’s modified Ea- likely to give rise to a functional consortium capable gle’s medium (DMEM) (Life Technologies) supplemented of altering all the biochemical pathways required for with 10% heat-inactivated fetal bovine serum (FBS) (Bio- tumor formation [20]. west) and antibiotics (penicillin, streptomycin) (Gibco®-Invi- We therefore reason that analysis of tumor heterogen- trogen). Cells were trypsinized/passaged every 2–3 days eity should not only be based on genetic alterations (mu- using TrypLE reagent (ThermoFisher Scientific). Cells were tations, amplifications, and translocations, among automatically counted using Countess cell counting cham- others) but should be complemented with functional ber slides (Invitrogen) and an Eve Automatic cell counter analysis based on protein expression and pathway ana- (NanoEntek), excluding dead cells by trypan blue (Invitro- lysis to reveal tumor heterogeneity at the phenotypic gen) staining. level [15, 28–31]. Phenotypic heterogeneity may also be caused by non-genetic alterations such as epigenetic changes or factors secreted from other cells within the Color-coding (transfection), clone isolation and tumor or tumor environment [14, 18, 27, 32, 33]. In fluorescence confirmation addition, conditions within the tumor, such as hypoxia, Ubc-StarTrack plasmids and the vector containing the oxidative stress or starvation are not reflected in genetic hyperactive transposase of the PiggyBac system (hyP- alterations, even though numerous adaptive changes Base) were kindly provided by Dr. López-Mascaraque. (e.g. metabolism) can be observed within the affected Ubc-StarTrack plasmids were generated as previously cells [28]. Since most of these changes are brought about described [35]. MDA-MB-231 cells were treated with by altered cell signaling pathways, the evaluation of ex- 20 μM chloroquine and transfected with 3.5 μg total pression levels and activity status (e.g. phosphorylation) UbC-StarTrack plasmid DNA plus 1.2 μg total HypBase of signaling factors is currently the best approach to (transposase) DNA, diluted in HBS (HEPES buffered sa- monitor functional tumor heterogeneity [14, 28, 32, 34]. line). CaCl2 was added to the medium to a final concen- In order to study intratumoral heterogeneity and tration of 150 mM. The medium was changed after six clonal cooperation on a functional level, we character- hours of incubation. ized the phenotypic features of individual clones isolated After transfection, the Ubc-StarTrack plasmids and from a breast carcinoma cell line, MDA-MB-231, and the transposase under the control of the ubiquitous compared these phenotypes with the parental cell popu- CMV promoter enter the nucleus. The Ubc-StarTrack lation. We observed differences in gene expression plasmids are flanked by two terminal repeat sequences among the different clones and were able to link many (ITRs) that are recognized and cut by the HypBase. The of these changes to alterations in cytokine-mediated released sequence is integrated into the genome of the intercellular signaling pathways. The importance of these transfected cell at TTAA repeat regions, allowing the alterations at the cellular level was demonstrated by cell to stably produce the fluorescent 161 protein (Add- clone-specific phenotypes in vitro and the metastatic po- itional file 1: Figure S1A). This system allows tracking by tential in vivo. Based on these results, we propose a color-coding of a complete population derived from the model in which clone-specific secretion and reception of transfected cell. Martín-Pardillos et al. BMC Cancer (2019) 19:666 Page 3 of 26 Forty-eight hours after transfection, cells were trypsi- Spheroid formation assay nized and isolated by fluorescence activated cell sorting 3D proliferation was evaluated by sphere growth measure- (High Speed Cell Sorter FacsAria (Becton Dickinson)). We ments at 1, 4 and 7 days after seeding. Ten thousand cells isolated cell clones expressing mT-Sapphire, EGFP, and per well were grown in ultra-low attachment 96 well monomeric Kusabira Orange (mKO) integrable plasmids. plates in the presence of 5% Matrigel® (v/v) after 10 min of Cells were individually seeded in a 96 well plate to obtain centrifugation at 2000 rpm. Pictures were taken with an clonal cell lines. Non-fluorescent cells were recovered to inverted microscope NIKON Eclipse TE2000–5 and proc- be used as a control (MDA-MB-231). The transfection ef- essed by ImageJ software to measure volume growth. ficiency was lower than 1%, showing that skewing of the parental cell line population was avoided. The clonal cell Caspase assay lines GFP C3, mKO E10, and Sapphire D7 and the MDA- Caspase activity was measured using the Caspase-Family MB-231 parental cell line were reauthenticated as de- Colorimetric Substrate Set (K132, Biovision). Five hun- scribed before. Fluorescence was confirmed on confocal dred thousand cells were seeded in a 100x17mm dish microscopy (Spectral Confocal Microscope FV1000 and trypsinized after 96 h in culture. The pellet was re- (Olympus) and FV10-ASW 4.2 software. For every fluor- suspended in Cell Lysis Buffer (Biovision) and protein escent protein the following excitation and emission wave- was quantified by Pierce™ BCA Protein Assay Kit (23, length settings were used: Ex.405-Em.520 (Sapphire), 225, Thermo Fisher). A positive control was included: Ex.488-Em.520 (GFP), and Ex.515-Em.527 (mKO). After MDA-MB-231 cell line treated with Camptothecin (208, every trypsinization the percentage of fluorescent cells 925, Merck) to a final concentration of 25 μg/mL for 24 was determinate by flow cytometry (LSR Fortessa (BD), h. 100 μg protein and pNA conjugated substrate to a software FCS Express DeNovo), ruling out the possibility final concentration of 0.5 mg/mL were added to each of fluorescent gene silencing. well in triplicate and incubated overnight at 37 °C. Ab- sorbance was measured at 405 nm with an Epoch Micro- Contrast phase image capture plate Spectrophotometer (BioTek). Forty-eight hours after seeding, phase contrast images were taken using an Olympus FSX100 microscope (amp- Spheroid invasion assay lification 40X) and FSX-BSW software. 3D invasion assay was performed by embedding 4-day- old spheres of each clone in a Matrigel-Collagen I mix- Cell count ture on a pre-coated 24-well plate with the same mix- Proliferation was evaluated by direct cell count. Fifty ture. Sphere invasion was analyzed by measuring the thousand cells were seeded in a six-well plate. Cells were invaded area at 24 and 48 h using the inverted micro- trypsinized and automatically counted at the following scope NIKON Eclipse TE2000–5 and ImageJ software. time-points: 0 (24 h after seeding), 24, 48, and 72 h. To calculate the cumulative population doublings (CPDs), 500,000 cells were seeded in a 100x17mm dish Migration assay and trypsinized after 7 days in culture. Cells were counted Five hundred thousand cells were seeded per well in 24- and seeded again, repeating the same process for 4 weeks well plates in triplicate. After 24 h, the medium was re- (28 days in culture). Cumulative population doublings placed by medium without FBS and maintained over- were calculated using the following formula: night. Then, a wound was made in the monolayer with a CPDs = [log(n° cells after 7 days) – log(n° seeded pipette tip, and the medium was replaced with complete cells)] / log102. medium. Pictures of the wounds were taken at 0 h and 8 h using an Olympus FSX100 microscope. Wound clos- Metabolic activity quantification (MTT assay) ure was measured using ImageJ software. Results repre- Metabolic activity was evaluated by the MTT (3-[4,5-di- sent the migrated distance between 0 and 8 h, expressed methylthiazol-25-yl]-2.5-diphenyl tetrazolium bromide; as a percentage relative to MDA-MB-231. Panreac-AppliChem) assay. 2500 cells in 150 μL were seeded in 96-well plates in triplicate and the MTT assay Transposon insertion event detection and was performed at time-point 0 (24 h after seeding), and characterization at 24, 48 and 72 h. MTT was added to the medium to a Paired-end sequencing raw reads were first processed to final concentration of 0.5 mg/mL and incubated for 3 h remove PCR duplicates using Fastuniq. Then, we trimmed at 37 °C. The medium was aspirated, and the formazan the specific Inverted Repeat (IR) associated with the trans- crystals were dissolved in 0.2 mL DMSO. Absorbance poson insertion process (GATTATCTTTCTAGGGTTAA, was measured at 595 nm using an Epoch Microplate length = 20) with Trimmomatic [36]. This step was re- Spectrophotometer (BioTek) and Gen 5 1.10 software. quired to later select those reads that were shorter than Martín-Pardillos et al. BMC Cancer (2019) 19:666 Page 4 of 26 the original read length (151), implying that they were All microarray data in this publication have been de- covering a transposition event. To map the reads, we posited in the NCBI Gene Expression Omnibus and are aligned the reads to the human reference genome g1k_v37 accessible through the GEO Series accession number [37] using Bowtie 2 and allowing 1 mismatch [38]. Later, GSE122008 (https://www.ncbi.nlm.nih.gov/geo/query/ we filtered the alignment with SAMtools [39] to select acc.cgi?acc=GSE122008). paired reads mapped unambiguously with a minimum Bioinformatic analysis was performed at the Statistics alignment quality of 30. The final step of the process relied and Bioinformatics Unit (UEB) of the Vall d’Hebron Re- on basic shell text processing tools (awk/grep/sed) to sub- search Institute (VHIR, Barcelona, Spain). A Robust set those pairs where one of the reads presented a length ≤ Multi-array Average (RMA) algorithm [41] was used for 131 (expected read length after removing the IR) and ex- pre-processing microarray data. Background adjustment, tract them from a specific position in the chromosome normalization and summarization of raw core probe ex- representing the first mapped genome base contiguous to pression values were defined so that the exon level the IR. This list of positions required an additional values were averaged to yield one expression value per prioritization step as the IR sequence was endogenously gene. Data were subjected to non-specific filtering to re- found in the human genome. To differentiate these arte- move low variability genes. Conservative thresholds were factual cases from genuine transposition events we took used to reduce possible false negative results. Selection advantage of one of the transposase properties: they do of differentially expressed genes was based on a linear not perform even cuts in both reverse/forward strands but model analysis with empirical Bayes modification for the a staggered cut that generates a duplication of 5 bases. A variance estimates [42]. To account for multiple testing, consequence of this is the representation of a specific in- P-values were adjusted to obtain stronger control over sertion event at two positions: n and n + 5. This criterion the false discovery rate (FDR), as described by the Benja- was used to depict actual transposition events and define mini and Hochberg method. The analysis of biological the set of insertion events that passed to the last step of significance was based on enrichment analysis against characterization. In this last process, we used BEDTools the Gene Ontology (GO; http://www.geneontology.org) map function and the annotation related to the human and KEGG (http://www.genome.ad.jp/kegg/) databases. genome reference selected to map the specific insertion positions to their genomic context [40]. DNA methylation analysis Microarray-based DNA methylation analysis was con- ducted with the Infinium MethylationEPIC BeadChip Analysis of the transcribed genome microarray (Illumina, San Diego, CA), that covers over 850, The microarray service was carried out at the High Tech- 000 CpG methylation sites (850 K). 1.5 × 106 cells were nology Unit (UAT) at Vall d’Hebron Research Institute seeded in a 100x17mm dish in triplicate. After 72 h cells (VHIR), Barcelona (Spain). Affymetrix GeneTitan micro- were trypsinized and centrifuged. Pellets were frozen until array platform and the Genechip Human Clariom D array DNA extraction. DNA quality checks, bisulfite modifica- cartridges were used for this experiment. This array ana- tion, hybridization, data normalization, statistical filtering, lyzes gene expression patterns on a whole-genome scale and beta (β) value calculations were performed as described on a single array with probes covering many exons on the elsewhere [43, 44]. The DNA concentration of the samples target genomes, and thus permits an accurate summary of was measured using the Quantifluor ONE dsDNA System gene expression. (Promega). A total of 500 ng of DNA samples were selected 1.5 × 106 cells were seeded in a 100x17mm dish in tripli- for bisulfite conversion with the EZ DNA Methylation™ Kit cate. After 72 h cells were trypsinized and centrifuged. RNA (Zymo Research). The Illumina Infinium HD methylation was extracted using the PureLink™ RNA Mini Kit (Thermo- protocol was followed for the hybridization to the Infinium Fisher). Quantification and assessment of RNA purity was MethylationEPIC BeadChips. Whole-genome amplification performed using a NanoDrop™ 2000/2000c Spectropho- and hybridization were performed on the BeadChips and tometer and confirmed according to the RIN (RNA Integ- followed by single-base extension and analysis on a HiScan rity Number) using RNA 6000 Nano Kit and 2100 (Illumina, San Diego, CA) to assess the cytosine methyla- Bioanalyzer (Agilent). Starting material was 200 ng of total tion states. The methylation score of each CpG was repre- RNA of each sample. Briefly, sense ssDNA suitable for la- sented as β value, and previously normalized for color bias belling was generated from total RNA with the GeneChip adjustment, background level adjustment and quantile WT Plus Reagent Kit from Affymetrix (Thermofisher-Affy- normalization across arrays. Probes and sample filtering in- metrix) according to the manufacturer’s instructions. Sense volved a two-step process for removing SNPs and unreli- ssDNA was fragmented, labelled and hybridized to the ar- able betas with a high detection P-value > 0.01. Sex rays with the GeneChip WT Plus Terminal Labeling and chromosome probes were also removed. After this filtering, Hybridization Kit from the same manufacturer. the remaining CpGs were considered valid for the study. Martín-Pardillos et al. BMC Cancer (2019) 19:666 Page 5 of 26 Cytokine quantification Protein extraction and immunoblotting 1.5 × 106 cells were seeded in a 100x17mm dish in tripli- Total protein extracts were generated using RIPA Lysis cate. After 72 h the medium was collected and filtered Buffer System (sc24948, SantaCruz Biotechnology) supple- through 0.22 μm filters (Merck Millipore) to remove float- mented with Protease Inhibitor Cocktail set III (539134) ing cells, apoptotic bodies and cell debris. The concentra- and Phosphatase Inhibitor Cocktail Set II (524625) from tions of the cytokines TNF alpha, IFN gamma, IL-1 beta, Calbiochem. Protein was quantified using a BCA Protein IL-2, IL-4, IL-5, IL-6, IL-8, IL-9, IL-10, IL-12p70, IL-13, IL- Assay kit (23,225, ThermoFisher). Protein extracts (10 μg 17A, GM-CSF, MCP1, MIP1a and MIP1b were determined per sample) were loaded onto SDS-PAGE gels and elec- using the FirePlex Human Key Cytokines - Immunoassay trophoretically transferred to PVDF membranes. The fol- Panel (ab229791) by Abcam FirePlex Service Lab. lowing primary antibodies were used: CD81 (sc-166,028, SantaCruz Biotechnology), TSG101 (ab83, Abcam) and b- Complementation assays actin (JLA20, Calbiochem). Goat anti-rabbit and anti- To test the presence of factors that could increase prolifer- mouse HRP secondary antibodies were from Pierce Ther- ation and metabolic activity, the media where the parental moScientific (31460) and Calbiochem (JA1200), respect- cell line or a clone were cultured for several days were re- ively. Immunodetection of proteins was performed using covered. The medium was tested by two independent as- Amersham ECL Western Blotting Detection Reagent (GE says (MTT and cell count). Exosomes were tested by Healthcare). MTT. “Donors cells” were cultured as follows: 17,600 cells/cm2 (medium complementation) or 10,400 cells/cm2 Transmission electron microscopy (exosome complementation) were seeded in 0.17 or 0.1 Exosomes were characterized by transmission electron mL complete medium/cm2, respectively. Twenty-four microscopy (TEM), using negative staining. Electron mi- hours before treatment the “receptor cells” were seeded: croscopy images were recorded on a T20-FEI Tecnai ther- 1250 cells in 150 μL were seeded in 96-well plates in tripli- moionic microscope operated at an acceleration voltage of cate (MTT assay) or 25,000 cells in 3 mL were seeded in 200 kV. Negative stained samples were prepared by drop- 6-well plates in duplicate (cell count). Seventy-two hours ping 20 μL of sample onto carbon coated copper grids after the “donor” seeding, the medium was recovered. The (200 mesh), which were then dried at room temperature medium was filtered through 0.22 μm filters (Merck Milli- and stained with phosphotungstic acid. Exosome size was pore) to remove possible floating cells and added to the measured using ImageJ software. “receptor cells” (complete medium, CM). The medium was also ultracentrifuged to extract exosomes, the super- Co-culture assay natant after two ultracentrifugations was recovered (sol- MDA-MB-231, GFP C3, mKO E10 and Sapphire D7 cell uble factors, SF) and the pellets were resuspended in PBS lines were cultured as individual cell lines and as combina- (exosomes). Exosome protein concentration was quanti- tions, with 2, 3 and 4 cell lines in co-culture. The co- fied by BCA Protein Assay kit (ThermoFisher). cultures started with equal percentages of every clone. Five To perform the direct cell count, donor medium was hundred thousand total cells were seeded in a 100x17mm added to new medium to a final concentration of 0, 10, dish and trypsinized after 7 days in culture. Cells were 25 and 50%. Cells were counted after 96 h’ treatment. counted and seeded again, repeating the same process for For the MTT assay, donor medium was added to new 4 weeks (28 days in culture). After every trypsinization the medium to a final concentration of 0, 25 and 50%. MTT percentage of every cell line per co-culture was determinate assay was performed at time-point 0 (24 h after seeding), by flow cytometry (LSR Fortessa (BD), software FCS Ex- 72 and 96 h. To study the effect of exosomes on cell press DeNovo). CPDs were calculated as previously de- metabolic activity, cells were treated with 50% complete scribed. Expected CPDs were calculated using data from medium and 50% soluble factors as a control. Exosomes the individual plates (single culture). Observed CPDs were were added to a final concentration of 5, 10, 20, 40, 80 calculated using data from every co-culture. and 160 μg protein/mL. MTT assay was performed at 72 h. Arsenite resistance assay Two hundred thousand total cells were seeded in a Exosome isolation 60x17mm dish. 24 h after seeding, cells were treated for Seventeen thousand six hundred cells/cm2 were seeded 90 min with 250 μM Arsenite (NaAsO2 in complete in 0.1 mL complete medium/cm2. After 72 h in culture, medium). After treatment, cells were washed twice with supernatant was recovered and sequentially centrifuged PBS and new complete medium was added. Cells were as follows: 500 rcf/10 min, 12,000 rcf/20 min and 100, counted after 72 h. 000 rcf/90 min (× 2). Finally, each pellet was resus- The co-culture experiment was performed in the pended in PBS. presence of arsenite. The co-cultures started with equal Martín-Pardillos et al. BMC Cancer (2019) 19:666 Page 6 of 26 percentages of every clone (100,000 cells per clone): and the experimental procedures were approved by the MDA-MB-231 + GFP C3 or MDA-MB-231 + mKO E10. Animal Experimentation Ethics Committee at the insti- 24 h after seeding, cells were treated as detailed before. tution (76/17 CEEA). Body weight and physical appear- The percentage of each cell line in the total population ance of the animals were monitored twice a week. The was detected by flow cytometry at seeding (day 0), 90 animals were euthanized by cervical dislocation follow- min (day 1) and 72 h (day 4) after treatment. ing the euthanasia standard operating procedure (SOP) of Laboratory Animals of the Vall d’Hebron University Invasion assay Hospital Animal Facility. All the in vivo studies were 8-μm pore inserts were covered with 65 μL Matrigel performed by the ICTS ‘NANBIOSIS’, more specifically, (356,231, Corning) at a final concentration of 1.5 mg/ by the CIBER-BBN in vivo Experimental Platform of the mL. Matrigel was incubated for 4 h at 37 °C. Once the Functional Validation & Preclinical Research (FVPR) Matrigel was polymerized over the Matrigel layer, 30,000 area (Barcelona, Spain). cells were seeded in 100 μL medium without FBS. The implanted cells were the MDA-MB-231 parental Complete medium (supplemented with FBS, which stim- cell line, three in vitro isolated clones (mKO E10, GFP C3 ulates the cells to cross the Matrigel layer) was added to and Sapphire D7), both individually and as a mix of the the well under the insert, covering the bottom of the in- four lines in a ¼ ratio. The parental cell line was included sert. Twenty-four hours after seeding, cells were fixed in the mix to provide factors probably not produced by with PFA 4% and counterstained with Hoechst. The top the clonal cell lines and which may be needed for tumor of the insert was cleaned with a cotton bud to remove growth and/or metastasis. All cell variants were confirmed the Matrigel and any cells that did not cross the layer. to be negative for viral, bacterial (including mycoplasma) Images were taken of the entire insert bottom where the and parasitic pathogens (VHIR Screening Humano Com- invading cells were located, using an Olympus FSX100 pleto). The test was performed by the external reference microscope (amplification 4.2X) and FSX-BSW software. laboratory Dynamimed Research Company. Prior to the Following the previous description, we performed two injection, the percentage of fluorescent cells was quanti- different invasion assays, with co-culture and with fied by flow cytometry, to confirm the injected cells medium complementation. Modifications to the previ- expressed the fluorophore and exclude the possibility of ous protocol are detailed as follows: fluorescent gene silencing. The tumor growth rate and invasive capacities were eval- (1) Co-culture: GFP C3 and mKO E10 were seeded as uated by implantation of 2.5 × 106 cells into the right ab- single cell lines and as an equal combination of dominal mammary fat pad (i.m.p.f.) of twelve animals per both cell lines. group (five groups). The tumor volume was measured by (2) Medium complementation: GFP C3 and mKO E10 caliper measurements twice a week and calculated accord- were seeded as single cell lines. Complemented ing to the formula D × d2/2, where D is the largest diameter complete medium was added to the well under the of the tumor and d the smallest one. All animals were eu- insert, covering the bottom of the insert. thanized 34 days after inoculation to compare the primary Complemented medium was obtained as detailed in tumor size and composition and the number and extent of “Complementation assays”. lung metastases between groups. The tumors and lungs were weighed, fixed with paraformaldehyde 4%, and later Animal experimentation processed for histopathological analyses (hematoxylin and Female athymic nude mice (Strain: Hsd:Athymic Nude- eosin staining). Foxn1nu) (ENVIGO, Spain) were kept in pathogen-free The metastasis growth rate of the MDA-MB-231 conditions and used at 7 weeks of age. Animals were and clonal cell lines was evaluated by intravenous randomly housed under Specific Pathogen Free (SPF) (IV) injection of 2.5 × 106 cells into the caudal tail conditions in autoventilated racks in groups of six. Two vein of 10 animals per group (five groups). All ani- enrichment elements were included, a cardboard tube mals were euthanized 36 days after inoculation. Ani- and a square nestlet for nesting and thermoregulation. mals underwent gross necropsy consisting of a The housing temperature ranged from 23 °C to 25 °C, macroscopic evaluation. Lungs were excised, weighed, relative humidity ranged between 47 and 55%. The cycle fixed and processed for histopathological analysis. gradually simulates twilight and sunset, giving 12 h of Immediately following dissection, the tumors and lungs light with an intensity of 300 Lux and 12 h of darkness were fixed for 24 h by immersion in paraformaldehyde 4%. for each 24 h period. Food and water were provided ad- After fixation, the tissue was dehydrated to enable embed- libitum. Animal care was handled in accordance with ding with paraffin. Five-micron-thick sections were cut the Guide for the Care and Use of Laboratory Animals from fixed, paraffin-embedded tissues and mounted of the Vall d’Hebron University Hospital Animal Facility, on poly-L-lysine-coated glass slides. Sections were Martín-Pardillos et al. BMC Cancer (2019) 19:666 Page 7 of 26 deparaffinized in xylene and rehydrated in graded al- (Tukey’s multiple comparisons test). Individual figure leg- cohol. The presence of clonal cells in the primary ends specify the test used in each case. Asterisks indicate tumor and metastases were detected by fluorescence significant differences when P-values are < 0.05 (*), < 0.01 on confocal microscopy (Spectral Confocal Micro- (**), and < 0.001 (***). scope FV1000 (Olympus)). To avoid spectral overlap- ping of the different fluorescent proteins, a Lambda scan was performed from 470 to 635 nm followed by Graphics spectral deconvolution using FV10-ASW 4.2 software. Graphics were created using Prism 5.0 software and Images were quantified using the program ImageJ. Adobe® Illustrator CC. In vivo zebrafish tumor xenograft assays Results Zebrafish (Danio rerio) embryos were generated by nat- Generation of color-coded clonal cell lines derived from ural mating of adult fish according to previously de- the MDA-MB-231 cell line scribed procedures (Westerfield, 2000). In order to In order to study tumor heterogeneity in a cell culture- remove surface pigmentation, embryos were incubated based model, we selected the triple negative human breast with PTU (N-Phenylthiourea, Sigma) at a 0.003% w/v con- cancer (TNBC) cell line MDA-MB-231. The presence of centration and maintained at 28 °C. For xenograft experi- genetically and phenotypically different clonal populations ments, 2 days postfertilization (dpf) zebrafish embryos has previously been described in this line [45]. To allow us were de-chorionized if necessary and anaesthetized with to study individual clonal populations and trace their fate 0.003% tricaine (Sigma) w/v in E3 water and positioned during tumor progression, cells were color-coded according on a 10 cm Petri dish coated with 1.5% agarose. Immedi- to the expression of fluorescent proteins [46, 47]. Genomic ately prior to injection, single cell suspensions of MDA- integration of fluorescent protein-encoding sequences was MB-231-GFP C3 cells and MDA-MB-231-mKO E10 cells achieved by the use of the UbC-StarTrack PiggyBac trans- were labelled with the lipophilic fluorescent tracking dyes poson system (Additional file 1: Figure S1a) [35, 46]. To CM-DiD and CM-DiI (Invitrogen) respectively, according generate the color-coded cell lines, we used different to manufacturer’s instructions. To remove unincorporated transposon-containing plasmids that encode for the fluores- dye, cells were centrifuged and rinsed twice with Dulbec- cent proteins mT-Sapphire, EGFP and mKO (Additional co’s phosphate-buffered saline (DPBS). Cells were kept on file 1: Figure S1b). From the clonal cell lines, we selected ice before implantation and implanted within 3 h. For GFP C3 (EGFP), mKO E10 (mKO) and Sapphire D7 (mT- microinjection, the individual cell populations (individual Sapphire) for further study, each of which expresses one injection) or the combination of both at equal numbers fluorescent protein (Additional file 1: Figure S1c). Addition- (co-injection) were loaded in borosilicate glass capillary ally, FACS-sorted non-fluorescent cells were pooled and needles (1 mm O.D. × 0.58 mm I.D, Harvard Apparatus) grown as a control population (named MDA-MB-231, par- and approximately 300 cells were injected into the duct of ental). The transfection efficiency was lower than 1%, show- Cuvier (DoC) using a micromanipulator and an IM 300 ing that skewing of the parental cell line composition was microinjector (Narishige) with an output pressure of 10 avoided. Finally, cell line authentication by DNA profiling psi and 0.03 ms injection time. After injection, embryos of short tandem repeats (STR) confirmed that all obtained were examined for the presence of a fluorescent cell mass clones were derived from the parental cell line MDA-MB- at the injection site in the DoC, and then transferred to 231 (Table 1). fresh PTU-containing E3 water and placed into an incuba- To confirm the presence of heterogeneity in the tor at 34 °C for up to 3 days post injection (dpi). Zebrafish MDA-MB-231 cell line and to validate our approach of embryos were photographed at 0 h post injection (hpi) studying tumor heterogeneity in a cell line-based model, and 72 hpi with a fluorescent microscope DMi8 (Leica) to we characterized different phenotypic features of the in- determine tumor cell dissemination. The number of dividual clonal cell lines. tumor cells disseminated in the tail vein of the fish was analyzed with the software ImageJ Fiji. For each condition, data are representative of at least three independent ex- Cell morphology periments, with at least 30 embryos per group. MDA-MB-231 is an epithelial cell line, characterized by elongated polygonal cells. Although all clonal cell lines Statistics maintained an epithelial morphology, the shape of the Prism 5.0 software (GraphPad Software Inc., La Jolla, CA) cells varied between the different clones, the clone mKO was used for statistics and data representation. Significant E10 being the most different to the others, with more differences were determined using unpaired t-test with polygonal and less elongated morphology (Fig. 1a, Add- Welch’s correction (equal SDs not assumed) and ANOVA itional file 2: Figure S2a). Martín-Pardillos et al. BMC Cancer (2019) 19:666 Page 8 of 26 Cell growth Cell growth was measured by cell counting after 72 h Cell growth was determined by counting the number of in culture. Cell growth is the net result of proliferation viable cells in an exponentially growing culture 72 h after and cell death. GFP C3 may have had increased apop- seeding. As shown in Fig. 1b, the growth rate of GFP C3 tosis that was compensated by higher proliferation; in was similar to that of the parental cell line. In contrast, the contrast, cell growth in low-growth-rate cell lines would clones mKO E10 and Sapphire D7 had a reduced growth be the result of decreased proliferation, rather than in- rate of up to 40 and 30%, respectively, compared to the creased apoptosis. parental cell line. These results were confirmed by meas- uring the metabolic activity of the cells via MTT assay Cell migration (Fig. 1c). To confirm these results were not a consequence The wound healing assay was used to measure cell mi- of 2D culture, we performed spheroid 3D culture. The re- gration, which is essential for many biological processes sults were confirmed, showing an even greater growth rate including tumor invasion and metastasis [48]. The ana- reduction in 3D culture than in 2D culture for mKO E10 lyzed clonal cell lines had similar migratory capacity to and Sapphire D7 (Fig. 1d-e). With the aim of determining the parental cell line (Additional file 2: Figure S2c-d). if this phenotype is heritable, we cultured the different cell lines for 28 days and calculated the cumulative popula- Invasion tions doublings (CPDs) after 7, 14, 21 and 28 days in cul- Spheres assembled in 3D cultures were embedded in a ture (Additional file 2: Figure S2b). Cell lines maintained Matrigel/collagen matrix and the expansion of the spheres stable proliferation rates over time. into the matrix was determined microscopically. Analysis of the obtained data revealed that, while the invasive cap- Apoptosis acity of the clones GFP and Sapphire was similar to the In order to assess if the reduction in cell growth was due parental cell line, cells from the clone mKO had a reduced to increased cell death, we determined the activation of invasion rate of up to 50% into the matrix (Fig. 1h-i). caspase 1 activity, known to induce pyroptosis, a lytic form of cell death. While the clone mKO E10 showed Tumorigenicity and metastasis of clonal cell lines some activation of caspase 1, a striking increase of cas- Tumorigenicity of clonal cell lines pase 1 activity was observed in the clone GFP-C3 (Fig. Having established that the isolated clonal cell lines dis- 1f ). The activity of caspases [2, 3, 6, 9] was measured as played phenotypic heterogeneity in vitro, we next analyzed a readout for the activation of pro-apoptotic signaling their tumorigenic potential in vivo. We orthotopically pathways. In comparison to the parental cell line MB- injected the clones individually (GFP C3, mKO E10, Sap- MDA231, cells of the clonal cell lines GFP C3, mKO phire D7 or the parental MDA-MB-231 cell line) and in E10 and Sapphire D7 exhibited significantly higher levels combination (mix of clones and parental cell line, 25% of these caspases (Fig. 1g). each) into the mammary fat pad of immunocompromised Table 1 Authentication profile for the selected clones vs. MDA-MB-231 parental cell line. Characterization of the clonal cell lines (GFP C3, mKO E10 and Sapphire D7) versus parental cell line (MDA-MB-231, also named ATCC® HTB-26) by STR profile evaluation (first column): number of shared alleles and percent match STR profile Expected: MDA-MB-231 GFP C3 mKO E10 Sapphire D7 ATCC® HTB-26 D5S818 12 12 12 12 12 D13S317 13 13 13 13 13 D7S820 8,9 8 8 8 8 D16S539 12 12 12 12 12 vWA 15,18 15 15 15,16 15 TH01 7,9.3 7,9.3 7,9.3 7,9.3 7,9.3 TPOX 8,9 8,9 8,9 8,9 8,9 CSF1PO 12,13 12,13 12,13 12,13 12,13 D21S11 30,33.2 30,33.2 30,33.2 30,33.2 30,33.2 Amelogenin x x x x x Number of shared alleles 14 14 14 14 Number of alleles in database 16 16 16 16 Percent match 87% 87% 87% 87% Martín-Pardillos et al. BMC Cancer (2019) 19:666 Page 9 of 26 Fig. 1 (See legend on next page.) Martín-Pardillos et al. BMC Cancer (2019) 19:666 Page 10 of 26 (See figure on previous page.) Fig. 1 Phenotypic characterization of MDA-MB-231, GFP C3, mKO E10, and Sapphire D7 cell lines. a Morphological evaluation by phase contrast images (scale bar = 200 μm). b Cell count after 72 h in culture, relative to MDA-MB-231 count. c Metabolic activity rate: quantification of MTT metabolization by metabolically active cells, absorbance relative to MDA-MB-231. d-e 3D growth rate: spheroid area at days 1, 4 and 7, relative to MDA-MB-231-day 0 quantification (d) and representative pictures (scale bar = 250 μm) (e). f-g Apoptotic analysis: absorbance relative to MDA- MB231. Expression of Caspase 1 (f) and Caspase 9, 3, 2 and 8 (g). h-i 3D invasion by spheroid cell migration: relative invasion (h) and representative picture at 0 and 48 h (scale bar = 250 μm) (i). Significant differences were determined using unpaired t-test with Welch’s correction. Asterisks indicate significant differences when P-values are < 0.05 (*), < 0.01 (**), and < 0.001 (***) nude mice. Surprisingly, despite the differences observed microscopy revealed that, although cells from all clonal for these clones in vitro, all of the clonal cell lines formed cell lines were able to form metastases, the metastasis primary tumors, and the tumor growth rates were similar count was significantly lower in the mKO E10 clonal cell for all individual clonal lines and the parental cell line line (Fig. 3a). (Fig. 2a-b). Similar results were observed when the cell To confirm these results, we determined the meta- clones and parental cell line were co-injected as a mixture static capacity of the individual cell lines in the lungs of (Fig. 2a-b). the mice after tail vein injection. Measurement of the Despite the similar tumor growth rates, we wondered metastatic areas within the lungs confirmed the reduced if the tumors would have morphological differences. capacity of mKO E10 to form metastases (Fig. 3b, Add- Pathological examination of the tumors revealed that the itional file 3: Figure S3). The lungs were weighed as an primary tumor cells from the different cells lines were of indirect measurement of the presence of metastases, and a similar shape, but with subtle differences. Cells within the results were confirmed by area measurement (Fig. the mKO E10-derived tumors displayed a more rhab- 3c). These results indicate that the inability of these cells doid, plasmacytoid shape, Sapphire D7 being the more to form metastases might be due to their homing or sur- compact tumor, formed by particularly small, round cells vival in the lungs. We think that the mKO E10 cell line (Fig. 2c). Thus, despite similar tumor growth rates, would be an appropriate subject of future research on tumor cell shape varied slightly, depending on the origin targeted therapies. of the clones. Hematoxylin/eosin stained lung sections confirmed The tumors formed after co-injection of the different inter-individual heterogeneity, in the percentage of lung clonal cell lines and the parental cell line were assessed area occupied by metastases formed by MDA-MB-231 on confocal microscopy. As shown in Fig. 2d-e, the (Mean ± SEM 25.14% ± 11.14), GFP C3, (19.42% ± 6.41), clonal composition varied depending on the section ana- mKO E10 (0.25% ± 0.20), Sapphire D7 (10.07% ± 2.95) lyzed, and the individual clones showed variable inter- and Mix (8.54% ± 3.78) (Fig. 3b). Photographs showing mixing (Fig. 2d). Figure 2e shows the number of the variability among individual animals are shown in fluorescent clonal cells and the parental cells (not ex- Additional file 3: Figure S3. pressing fluorescent proteins) in the analyzed sections. In addition to determining the metastatic capacity of Interestingly, the parental cells were predominant, while the individual clones, our experimental setup of the the clones GFP C3 and Sapphire D7 were detected in mixed clonal populations also allowed us to determine if various areas of the tumor and few mKO E10 cells were macrometastases were formed from individual or mul- observed (Fig. 2d-e). Finally, cell counting of different tiple cells. We used confocal microscopy to analyze the sections of the whole tumor revealed that the majority metastases induced by an equal mix of the parental cell of fluorescent clone cells within the tumor were either line and clonal cell lines (Fig. 3d). mKO E10 cells were GFP or Sapphire positive (46.96 and 49.74%, respect- not found in the resulting metastases. Fluorescence re- ively) and only 3.3% of the cells were derived from the vealed the metastases could be formed by one (Fig. 3d) clone mKO E10 (Fig. 2e). Considering that all the clones or more than one clone (Fig. 3e). However, even in me- had similar growth rates when injected individually into tastases formed by more than one clone, we found a the mice, this finding indicates that the interplay among main component formed from a single clone. the different clones determines their fate during tumor progression. Molecular characterization of clonal cell lines by gene Metastatic tumor formation capacity expression, methylomics and cytokine expression We next analyzed the ability of the clones to form metas- Gene expression profiling by microarray analysis tases derived from the orthotopic tumor. The lungs of the To uncover the molecular mechanisms of interplay be- sacrificed mice were stained with hematoxylin/eosin, and tween different clones in a heterogeneous population, Martín-Pardillos et al. BMC Cancer (2019) 19:666 Page 11 of 26 a b c d e Fig. 2 (See legend on next page.) Martín-Pardillos et al. BMC Cancer (2019) 19:666 Page 12 of 26 (See figure on previous page.) Fig. 2 Malignancy of individual clones and tumor heterogeneity. a-c Tumor growth over 34 days after orthotopic injection of MDA-MB-231, GFP C3, mKO E10, or Sapphire D7 clonal cell lines and an equal mix of the previous cell lines (mix): (a) volume relative to MDA-MB-231 tumors (n = 12 per group), (b) macroscopic images of the tumors (scale bar = 0.5 cm) and (c) hematoxylin/eosin stains of tumor sections (scale bar = 50 μm) (three different animals per cell line: ID1, ID2, ID3). d-e Tumor heterogeneity: (d) representative fluorescent images (scale bar = 50 μm) of three independent sections of tumor formed by MDA-MB-231 parental and GFP C3, mKO E10 and Sapphire clonal cell lines (mix) and (e) quantification of tumor composition formed by the clonal cell lines (n = 6). Significant differences were determined using unpaired t-test with Welch’s correction. Asterisks indicate significant differences when P-values are < 0.05 (*), < 0.01 (**), and < 0.001 (***) microarray analysis was performed to determine the detected between the different clonal cell lines and the relative alterations in gene expression profiles of the parental cell line. Statistically significant differences in clonal cell lines in comparison to the parental cell line gene expression compared to MDA-MB-231 cells were MDA-MB-231(Fig. 4). determined, separated into up- and downregulated genes Variation at the level of individual transcripts was quan- and displayed as a Venn diagram (Fig. 4b-c). The most tified by identifying which genes had a statistically signifi- striking differences were observed for the clone mKO E10: cant difference in expression level from average (p < 0.1) 1017 genes were downregulated, and 288 genes were up- and expressing these in a heatmap. As depicted in Fig. 4a, regulated, while the clonal cell lines GFP C3 and Sapphire significant differences in the gene expression profiles were D7 had fewer regulated genes (GFP C3, 262 down and 71 Fig. 3 Metastatic capacity of cell lines and an equal mix of all cell lines. a Absolute count of micrometastases derived from orthotopic primary tumors, 34 days after injection of tumor cells. b-c Metastasis growth induced 35 days after intravenous injection of tumor cells: (b) lung area occupied by metastases, and (c) lung weight, an indirect measurement of lung tumors. Significant differences were determined using unpaired t- test with Welch’s correction. Asterisks indicate significant differences when P-values are < 0.05 (*), < 0.01 (**), and < 0.001 (***). d-e Representative fluorescent images (scale bar = 50 μm) of four independent sections of lung metastasis formed by MDA-MB-231 parental and GFP C3, mKO E10 and Sapphire clonal cell lines (mix) Martín-Pardillos et al. BMC Cancer (2019) 19:666 Page 13 of 26 up; Sapphire D7, 106 down and 40 up), in comparison to multiple insertions of the same transposon into the gen- the parental cell line (Fig. 4b-c). omic DNA, but, excepting one insertion in the Sapphire To elucidate which alterations in gene expression D7 clone, these insertions were located in non-coding could potentially account for the observed phenotypic genomic regions (introns). Analysis of the microarray differences in vitro and in vivo, we performed bioinfor- data confirmed none of the insertions affected the corre- matic analysis. We narrowed down the number of po- sponding gene expression with respect to MDA-MB-231 tential candidates using gene ontology and KEGG parental cell line (Table 2). pathway analysis to identify the pathways and biological In summary, even though all clonal cell lines were de- processes that were altered in the different clones (Fig. rived from the same parental cell line, their mRNA ex- 4d-e). Gene ontology (GO) revealed both common and pression profiles varied widely, and this variation was not clone-specific changes in biological processes (Fig. 4d). caused by genomic mutations via transposon insertion. Most of the changes occurred in the expression levels of We can therefore conclude that heterogeneity within the genes involved in signaling cascades. Interestingly, many MDA-MB-231 is present at the level of gene expression. of the genes identified within the class of signaling fac- tors also belonged to the GO term “cytokine-mediated Epigenetic alterations of gene expression signaling pathways”, the third most abundant class of Epigenetic mechanisms, such as DNA methylation, pathway affected in the analysis. Therefore, intercellular are important regulators of gene expression. To estab- communication via cytokines and their respective down- lish if the observed changes in gene expression were stream signaling pathways displayed the greatest differ- due to DNA methylation modifications, we deter- ences between clonal cell lines and the parental MDA- mined the genomic methylation pattern of the indi- MB-231 cells. Following this line, we were able to iden- vidual clones and the parental cell line (Fig. 5). To tify cytokines and cytokine receptors that showed alter- extrapolate data from these measurements that are ations that were either common to the parental cell line probably responsible for the observed changes in gene or clone-specific. The detection of specific signaling expression, we focused our analysis on CpG islands pathways in the analysis, namely NOD, Toll-like recep- and shore regions of gene promoters. The most dif- tor and JAK-STAT signaling pathways, was also due to ferentially methylated CpGs in these regions are the cytokines and their receptors, acting upstream of the shown in the heatmap in Fig. 5a. The comparative respective pathways. KEGG pathway analysis confirmed analysis of the methylation signature revealed specific, these findings (Fig. 4e). characteristic methylation patterns for the individual To elucidate why the clone mKO E10 had a non- clones and for the parental cell line MDA-MB-231. metastatic phenotype, we matched our dataset with the To determine if these changes caused the observed Human Cancer Metastasis Database. This analysis re- alterations in gene expression, we analyzed the associ- vealed that of the 658 genes with altered mRNA expres- ation between gene expression levels and the methyla- sion in the clone mKO, 95 have been linked to tion status of their promoter regions. As represented metastasis (Additional file 4: Table S1). Among the genes in Fig. 5b-c, changes in the methylation patterns were with the biggest alterations compared to the parental cell linked to discrete changes in the expression of the re- line was TGF-beta receptor type-1 (TGFR-1, log2FC spective genes. Downregulated gene expression corre- -0.9). TGF-beta receptor signaling has been strongly lated to hypermethylation of the promoter in 40% linked to metastatic disease progression in breast cancer, (GFP), 19% (mKO) and 25% (Sapphire) of cases (Fig. so downregulation of TGFR-1 in the clone mKO might 5b). We also observed changes in gene expression have been the reason for the non-metastatic phenotype that were common to the different clonal cell lines. [49, 50]. Detailed analysis of these and other potential Four altered methylation patterns could be linked to factors will be required in the future to determine how similar changes in gene expression in two of the their altered expression can affect metastatic disease clonal cell lines, and two altered methylation patterns progression. were linked to gene expression changes in all three clones (Fig. 5d). In summary, the observed alterations in gene expres- Transposon sequencing sion patterns could indeed be partly linked to different In order to test if the observed changes in gene expres- epigenetic states of the cell. Notably, changes in the sion profiles reflected a natural heterogeneity within the methylation patterns often affect the expression of genes MDA-MB-231 cell line, or if these changes were induced that act upstream of signaling cascades. It is therefore by genomic alterations through the insertion of transpo- plausible to assume that changes in the methylation pat- sons, we performed genomic DNA sequencing, mapping tern of core regulatory genes may play a crucial role in the exact insertion sites of the transposons. We detected the observed alterations in gene expression (Fig. 5e-f ). Martín-Pardillos et al. BMC Cancer (2019) 19:666 Page 14 of 26 Fig. 4 (See legend on next page.) Martín-Pardillos et al. BMC Cancer (2019) 19:666 Page 15 of 26 (See figure on previous page.) Fig. 4 Transcriptional profile of MDA-MB-231 clonal cell lines. a Heatmap for P-value under 0.05 and absolute logFC> 1.5 of selected clonal cell lines. b-c Venn diagrams of downregulated (b) and upregulated (c) transcripts of selected clonal cell lines vs. MDA-MB-231. d Gene ontology (GO) with differentially expressed genes (logFC> 1 and FDR < 0.05) of each clonal cell line vs. MDA-MB-231 cell line. e KEGG with differentially expressed genes (logFC> 1 and FDR < 0.05) of each clonal cell line vs. MDA-MB-231 cell line Cytokine arrays single culture (Fig. 1b); however, the increase of CPDs To confirm these alterations of cytokine-mediated sig- for MDA-MB-231 in co-culture (Fig. 7a) induced a sig- naling pathways at the protein level, we used cytokine nificant decrease of the clone GFP C3 in the total popu- arrays to analyze the composition of secreted factors lation (Fig. 7b). After 28 days in culture, the percentage in the individual clones (Fig. 6a). Among the 17 key of MDA-MB-231 in the population was 28% higher than cytokines tested, quantifiable amounts were detected expected, while the observed percentage of GFP C3 for IL8, IL6, MCP1(CCL2), GM-CSF and IL13. (14.72%) was lower than expected (39.94%) (Fig. 7a). Among those, MCP1 was reduced more than 10-fold These results were confirmed in experiments in which in all clonal cell lines compared to the parental cell only the parental cell line and the clone GFP-C3 were line. Weaker but still significant alterations were de- combined (Additional file 5: Figure S4b-c). In all combi- tected for GM-CSF and IL6. Clone-specific differences nations involving GFP C3, the observed percentages of were detected for IL8, which was strongly reduced in MDA-MB-231 increased and GFP C3 decreased. the clones mKO E10 and GFP but unaltered in the clone Sapphire, and higher levels of IL13 were de- tected in the GFP C3 and mKO E10 clonal cell lines. Resistance to stress in co-culture experiments To complete this analysis, we merged two datasets We show in Additional file 5: Figure S4b how the pres- (cytokine expression and mRNA receptor expression) ence of clonal cell lines with slower proliferation in the to elucidate which cytokine-cytokine receptor interac- total population tends to decrease as time in co-culture tions had observable changes between the parental increases. In the absence of stressors, such cells would cell line and the different clonal cell lines (Fig. 6b). be outcompeted after several weeks in co-culture. We This comparison revealed that cytokine-mediated sig- aimed to test why these slower-proliferating cells were naling pathways might be affected by altered expres- still present when we sorted the parental cell line. We sion of cytokine receptors, amount of cytokines hypothesized that cells could persist in the population released, or by simultaneous changes in cytokines and was subject to a stressor that provided any benefit to the cytokine receptors. slower cell line or any disadvantage to the faster- proliferating cell lines. In Additional file 5: Figure S4d Interplay between the different clones we show that mKO E10 was the clonal cell line most re- Co-culture experiments sistant to arsenite treatment, followed by GFP C3 and Co-culture experiments were performed to test if the co- Sapphire, MDA-MB-231 being the least resistant. We operation between different clones might affect the wondered how the co-cultured population would evolve growth rate of the individual clones. MDA-MB-231, in the presence of this stressor. Additional file 5: Figure GFP C3, mKO E10 and Sapphire D7 cell lines were cul- S4e-f shows how the percentage of clones (GFP C3 or tured individually and in different combinations for 28 mKO E10) in co-culture was significantly higher (and days (Fig. 7, Additional file 5: Figure S4a-c). In line with consequently, the percentage of MDA-MB-231 cells was the results obtained in the in vivo mouse model (Fig. lower) in the presence of arsenite than in the absence of 2e), the mKO clonal cell line was progressively reduced arsenite. The clone population increases because they and after 28 days of co-culture was severely depleted, are more resistant than MDA-MB-231. We have in- representing just 3% of the total population (Fig. 7a, cluded this as an example of a stressor, but other Additional file 5: Figure S4a). The decrease of mKO cells stressors could also affect the resistance of the clones, when mixed with the parental cells may indicate that its for example, resistance to growth from a single cell (sub- growth is overtaken by that of other clones. This result, cloning). This allows us to explain why slower- and the result for Sapphire D7, matched the expected proliferating clones were present in the MDA-MB231 values determined by the calculated cumulative popula- cell line when we selected them, as they could be more tion doublings (CPDs) from culture of the individual cell resistant to a particular stressor than the whole popula- lines (Fig. 7b). In contrast, the percentage of cells from tion, and thus persist. Therefore, our results show that the clone GFP C3 was significantly lower than the ex- maintenance in the total population it is not only pected values (Fig. 7a). The proliferation rate of MDA- dependent on the proliferation rate; environmental con- MB-231 and GFP C3 was not significantly different in ditions can affect the presence of the clones. Martín-Pardillos et al. BMC Cancer Table 2 Location of transposons in the genome of MDA-MB-231 clonal cell lines. Chromosome location (chromosome number and position, region type, strand, gene name, intron/exon, number of isoforms with the insertion, analyzed or not in the microarray, microarray probe and Clone vs MDA-MB-231 logFC in microarray analysis Cell line Chromosome posL posR Type Strand Name Relative Isoforms N° isoforms Microarray Microarray probe Clone vs location with the analysis MDA-MB231 (2019) 19:666 insertion logFC MDA-MB-231 – – – – – – – – – – – – GFP C3 10 33,272,491 33,272,496 Gene Reverse ITGB1 Intron 2 0 Yes TC1000010265.hg.1 Non significant lncRNA Forward RP11-462 L8.1 Intron 2 0 No – mKO E10 10 69,810,770 69,810,775 Gene Reverse HERC4 Intron 8 2 Yes TC1000010831.hg.1 Non significant 10 71,001,553 71,001,558 Gene Forward HKDC1 Intron 2 0 Yes TC1000007890.hg.1 Non significant 10 41,117,156 41,117,161 Gene Reverse PPP1R14D Intron 2 0 Yes TC1500009127.hg.1 Non significant Sapphire D7 5 170,347,638 170,347,643 Gene Forward RANBP17 Intron 10 0 Yes TC0500009411.hg.1 Non significant 8 145,570,348 145,570,353 Intergenic – – – – – Yes – Non significant 9 113,552,211 113,552,216 Gene Forward MUSK Intron 3 0 Yes TC0900008415.hg.1 Non significant 9 78,576,637 78,576,642 Gene Forward PCSK5 Intron 3 0 Yes TC0900007618.hg.1 Non significant 10 112,724,387 112,724,392 Gene Forward SHOC2 Exon 5 3 Yes TC1000008908.hg.1 Non significant X 106,242,426 106,242,431 Gene Reverse MORC4 Intron 3 0 Yes TC0X00010459.hg.1 0.484 Page 16 of 26 Martín-Pardillos et al. BMC Cancer (2019) 19:666 Page 17 of 26 Fig. 5 (See legend on next page.) Martín-Pardillos et al. BMC Cancer (2019) 19:666 Page 18 of 26 (See figure on previous page.) Fig. 5 DNA methylation profile of MDA-MB-231 clonal cell lines. a Heatmap of the most differentially methylated CpGs (delta ≥0.20) localized at promoter, CpG island and shore regions (12,201 CpGs). b-c Identification of genes whose expression changes could be explained by methylation differences: (b) downregulated vs hypermethylated and (c) upregulated vs hypomethylated. d Venn diagrams of methylation gene changes of selected clonal cell lines vs. MDA-MB-231. e-f Gene Ontology (e) and KEGG (f) with differentially methylated genes (delta ≥0.20) at promoter and CpG island and shore Clonal cooperation on migration and invasion assays invasive capacity of the clone GFP C3 increased two- In addition to clonal cooperation by secreted factors, fold when seeded together with the clone mKO E10 (Fig. cells of a heterogenous tumor might use feature comple- 7d). mentation to achieve complex cellular processes like in- To determine if the increased invasion was the result vasion into the cellular matrix – for which not only cell of physical interactions among clonal cell lines or se- migration, but also digestion of the extracellular matrix creted stimulatory factors, we performed the invasion and chemotactic abilities are required. assay with GFP C3 and mKO E10 clonal cell lines To test if our cell line model could reproduce such a seeded individually, adding complemented medium complementary behavior between different clones, we under the insert. We compared the invasion efficiency in tested the invasive capacity of the phenotypically more the presence of their own complemented medium or distinct clones, mKO E10 and GFP C3, alone and in medium complemented by the other clonal line, and combination in the Matrigel invasion assay (Fig. 7c-d). found that invasion also increased in the presence of Comparison of the total number of invading cells re- medium complemented by the other cell line (Fig. 7e). vealed that cells in the clone mixture were significantly To confirm these data in a more physiological context, more effective at invading than the cells from the clone we performed an assay to measure cell dissemination in GFP C3 alone (Fig. 7c). To estimate the contribution of the tail of zebrafish. Cells of the different clones were each clone to this increased invasive capacity, we must injected alone (300 cells) or in combination (150 GFP take into consideration that although the total number C3 cells + 150 mKO E10 cells) into the duct of Cuvier of seeded cells was the same as when we seeded a single and the number of cells in the tail after 72 h was deter- clone (30,000 cells), the number of cells from each clone mined (Fig. 7f-i). We found that co-injection of both in the mix represents only 50% of the cells when seeded clonal cell lines resulted in significantly more cells in the individually (15,000 GFP C3 cells + 15,000 mKO E10 tail of the fish than injection of individual clonal cell cells). Normalization of the data revealed that the lines (Fig. 7f-g). In addition, a significant increase in the Fig. 6 Cytokine expression. a Heatmap with clustering of released cytokines from the previous clonal cell lines. b Released cytokine mRNA and protein levels vs receptor levels Martín-Pardillos et al. BMC Cancer (2019) 19:666 Page 19 of 26 a b c d e f g h i Fig. 7 (See legend on next page.) Martín-Pardillos et al. BMC Cancer (2019) 19:666 Page 20 of 26 (See figure on previous page.) Fig. 7 Clone interactions among parental and clonal cell lines in co-culture for 28 days (a-b). a Expected population percentage of total population calculated using cumulative population doublings (CPDs) per day vs. observed population percentage. b CPDs/day of single culture vs. CPDs/day of the same cell line in co-culture. c-e Cooperation of GFP C3 and mKO E10 in invasion assays: (c-d) Evaluation of cooperation on the invading capacity of single clones vs co-cultured clones: (c) Total number of invading cells and (d) Relative ratio of invading cells based on the cell number per clone cultured individually. (e) Invading capacity of cells conditioned with their own medium or medium from another clonal cell line. Relative ratio of invading cells based on the cell number per clone cultured with their own conditioned medium. (f-i) In vivo clonal cooperation. Cooperation of GFP C3 and mKO E10 in in vivo migration assays (Zebra fish model). Disseminated cells from the duct of Cuvier to the tail 72 h after injection: individual injection of GFP C3 and mKO E10 and co-injection of both clonal cell lines: (f) Total number of disseminated cells, (g) representative graph (scale bar = 70 μm), (h) relative ratio of invading cells based on the initial cell number per clone, (i) analysis of the coexistence of clones in the co-injected fish: percentage of co-injected fish with GFP C3 cells only, mKO E10 cells only, or both. Significant differences were determined using unpaired t-test with Welch’s correction. Asterisks indicate significant differences when P-values are < 0.05 (*), < 0.01 (**), and < 0.001 (***) disseminating capacity of both clones was observed faster in the presence of MDA-MB231, comparing the (four-fold GFP and two-fold mKO) when the number of percentage of observed populations (experimentally) injected cells was normalized for the respective condi- and expected population (calculated from proliferation tions (Fig. 7h). The efficiency of the injection did not rate when the clonal cell line is cultivated individually). show significant differences between individual and co- We observed this phenomenon when MDA-MB231 injected cells (Additional file 5: Figure S4g), excluding and GFP C3 were cultured together and in the presence injection efficiency as the cause of the differences in in- of a third and even a forth clonal cell line. As we have vasion. Determination of the number of fish that con- shown, GFP C3 proliferation does not increase with tained either both or only one of the clones after co- medium complementation from other cell lines, how- injection revealed that 40% of fish had a combination of ever all the other cell lines do (Fig. 8a-d). The lack of GFP C3 and mKO E10, 9.95% had GFP C3 only and response of GFP C3 to medium complementation can 10.29% had mKO E10 only (Fig. 7i), confirming that modify the expected evolution of the clonal compos- communication between the cells of the different clones ition: instead of being maintained at the same per- and physical interaction account for the positive effect in centage as MDA-MB-231 (expected), GFP C3 is dissemination. progressively reduced. Effect of secreted factors on cell growth Extracellular vesicle induced cell growth To identify We aimed to determine if cell growth could be modified factors responsible for the increased metabolic rates upon by secreted factors. treatment with conditioned medium, we isolated extracel- lular vesicles from the conditioned medium by ultracentri- Cell growth in the presence of conditioned medium fugation. We confirmed that all the clonal cell lines were To study the effect of secreted factors (including cyto- able to release factors that increase metabolic rate (Fig. kines and vesicles), conditioned medium (cell culture 7a-d) and that selected clonal cell lines were also able to supernatant obtained after 96 h of culture) was applied release exosomes (Additional file 6 Figure S5). We se- in different concentrations to the individual cell lines lected the parental cell line as the exosome donor based and the growth rates were determined by MTT. The on our previous results: (1) analysis of vesicle size and growth of the clonal cell lines mKO E10 and Sapphire protein composition did not reveal major differences in D7 as well as the parental cell line increased with the the vesicles isolated from the parental cell line and the in- conditioned medium from all other cell lines in a dose- dividual clones (Additional file 6: Figure S5a-e) and (2) cell dependent manner (Fig. 8a-d). These results therefore lines responded positively to the addition of the factors, suggest that paracrine signaling promotes cell growth independently of the donor (Fig. 8a-d). Isolated extracellu- (Fig. 8a-d). While this growth-supporting effect of con- lar vesicles from the parental cell line were then added to ditioned medium has long since been known, it was in- the parental cell line and the effect on cell growth was de- teresting that the clonal cell line GFP C3 did not termined. Addition of increasing amounts of vesicles to significantly respond to any of the applied media. We the cells modified the metabolic rate of the receiving cells could suppose that GFP C3, as the faster clonal cell (Fig. 8e). The maximal increase in metabolic activity was line, does not require additional supplements from observed with concentrations up to 1.25 μg exosome pro- other cell lines. The consequence of this lack of re- tein/mL. Higher concentrations of vesicles decreased the sponse is shown in co-culture with cells with a similar effect on metabolic activity. Thus, extracellular vesicle growth rate. In Fig. 7a and Additional file 4: Figure fractions contain factors that can modify the metabolic ac- S4b-c we see how the GFP C3 population decreases tivity of the receiving cells in a dose-dependent manner. Martín-Pardillos et al. BMC Cancer (2019) 19:666 Page 21 of 26 Fig. 8 Clone communication by released soluble factors. a-d Metabolic activity by MTT assay of clonal cell lines, 72 h after medium complementation with supernatant extracted from MDA-MB-231, GFP C3, mKO E10 and Sapphire D7 cell line. Each graph corresponds to a cell line donor medium, (a) MDA-MB-231, (b) GFP C3, (c) mKO E10 and (d) Sapphire D7. Color bars indicate the receiving cell line. e Metabolic activity measurement of MDA-MB-231 by MTT assay after complementation with non-complemented medium (C-), 50% complemented medium (CM) and exosomes. Significant differences were determined using unpaired t-test with Welch’s correction. Asterisks indicate significant differences when P-values are < 0.05 (*), < 0.01 (**), and < 0.001 (***) Martín-Pardillos et al. BMC Cancer (2019) 19:666 Page 22 of 26 Discussion epigenetic alterations accumulated. In such a situation, her- Breast cancer is the most frequently diagnosed neoplasm itability could be tested by studying the phenotype of sub- in women and a leading cause of cancer death [1, 2]. clones derived from the original clones, having cultured the Intertumor heterogeneity is one of the greatest issues in original clones for several months. oncology, but it is not the only factor that makes patient Moreover, in co-culture (Fig. 7, Additional file 5: treatment difficult: variation within a single tumor – Figure S4a-c) the parental cells grew faster than the intratumor heterogeneity – is also a key factor in deter- single clones, indicating that the parental cell line mining the best treatment [51]. harbors a sum of clones that may outperform single Cellular heterogeneity in established breast cell lines, clones through local or secreted factors, as yet not including the MDA-MB-231 cell line, has been described identified. These findings support the concept of the previously [45, 52]. In this study we corroborate clonal whole is greater than the sum of its parts and de- heterogeneity in the MDA-MB-231 cell line of breast scribe a model where the interplay of clones confers carcinoma, after the selection of clones marked with aggressiveness, and which may allow the identification fluorochromes, using the StarTrack system (Additional of factors involved in cellular communication and me- file 1: Figure S1c). Clonal cell lines were selected to tastasis. Thus, clonal heterogeneity allows the malig- study the relevance of clonal communication in tumoral nant cell line to acquire the greatest malignant heterogeneity. The results from the experiments in vivo potential. and in vitro demonstrated that despite all the observed Conventional models propose that each metastasis phenotypic differences, the parental cell line, composed originates from a single tumor cell [55–57]. However, of a mix of an undetermined high number of clones, was recent studies using mouse models of cancer have equal or superior in terms of malignancy (tumor com- demonstrated that multiple subclones undergo poly- position and metastasis) to the individual clonal cell clonal seeding and demonstrate interclonal cooper- lines in all tested conditions. Because the parental cell ation between multiple subclones [7, 58]. Our results line is a mix of different clones, these results led us to confirm (Fig. 3d) that metastasis could be formed ei- the hypothesis that the interplay between different clonal ther by a single (Fig. 3e) or several clonal cell lines populations within a heterogeneous population, by inter- (Fig. 3d). However, even in cases of metastasis formed cellular communication and/or feature complementa- by several clonal cell lines, each metastasis contained tion, provides advantages over isolated clonal cell lines. one predominant clone (Fig. 3d). Comparative studies Our results show that the mix of three clonal cell lines indicate monoclonal patterns of seeding, suggesting did not recover the features of a multiclonal population. that clones compete to metastasize. However, poly- For the study of tumor heterogeneity, the subcloning clonal seeding, in which multiple clones from the pri- method limits the number of clones and makes reconsti- mary tumor seed the same metastasis, is also tution of a highly multiclonal population difficult. There- observed, indicating subclones might cooperate as fore, the mix of several clones cannot be equated to a well as compete to metastasize [7, 59]. In our model multiclonal population; rather, a parental cell line must the cells were injected as a mix of single cells, there- be included to provide a true multiclonal population. fore the metastasis formed by more than one cell line However, subcloning may be a suitable method for the originated from several cells that reached the lung to- independent study of the single clones that constitute a gether, demonstrating that the cells physically interact cancer cell line or tumor. to form the metastasis. Selected clonal cell lines showed variations in cell morph- Several studies call into question the theory of clonal pro- ology (Fig. 1a, Additional file 2: Figure S2a), just as patholo- gression by the progressive accumulation of genetic alter- gists observe variation between different regions of the ations and selection of more aggressive clones, supporting same tumor [53, 54]. The clonal cell lines also showed dif- instead the proposed theory of clonal cooperation between ferences at the biological (Fig. 1), tumorigenic (Figs. 2, 3, tumor clones [20–23, 25–27]. Tumor multiclonality is also Additional file 3: Figure S3) and molecular levels, and in supported by the field cancerization theory [60, 61], which gene expression (Fig. 4), DNA methylation pattern (Fig. 5) states that there are many genetic alterations in the normal and cytokine expression (Fig. 6), with clear functional het- tissue surrounding tumors that can give rise to independent erogeneity. Parental cells showed a faster growth rate (Fig. clones. Similarly, supporting interpretations can be drawn 1b-e) and biological aggressiveness – a greater presence in from the stem cell hypothesis, as diverse clones can derive metastases – than the clones we studied (Fig. 3, Additional from more than one pluripotent stem cell [62, 63], and the file 3: Figure S3). Our results show that the phenotypes of Big Bang model of colorectal tumor growth where the the single derived clones were heritable in the conditions tumor grows predominantly as a single expansion popu- used in this study; however, culturing subclones for a longer lated by numerous intermixed subclones [64]. Clonal co- time could induce phenotypic changes if molecular or operation has recently been suggested in studies of single Martín-Pardillos et al. BMC Cancer (2019) 19:666 Page 23 of 26 cell sequencing [62, 65, 66]. The present study further sup- level, and factors involved in the communication between ports the idea that there are several clones that together tumor cells and microenvironmental or inflammatory cells. confer the properties of malignancy, thus strengthening the In summary, these data support that functional clonal concept of clonal cooperation, whereby clones synergistic- heterogeneity does not always reflect genetic heterogeneity. ally provide certain selective advantages for proliferation, Inhibition of clonal cooperation, by blocking cytokines or resistance to apoptosis, induction of angiogenesis, and other factors involved in the formation of clusters and the interaction with environmental factors and inflammatory interplay with environmental cells may represent a change cells [20, 21, 30, 63]. in the therapeutic paradigm to prevent the development of Our results show that cells are able to interact and the metastasis [21]. coexistence of clonal cell lines resulted in a positive ef- fect (Fig. 7c-h). We could conclude that physical inter- Conclusions action between clones (Fig. 7c-h) and secreted factors Interplay between clones confers malignancy, supporting (Fig. 7e, Fig. 8a-e) favored the tumorigenic capacities of the concept of clonal cooperation, whereby clones synergis- the cells. Cells are able to secrete factors (Fig. 6a) and tically provide certain selective advantages by physical inter- send messages to surrounding and distant cells (para- action, soluble factors and extracellular vesicles. crine signaling). Studies of cellular communication in breast cancer have demonstrated tumor cells’ abilities to Additional files secrete factors increasing breast cancer cell proliferation and metastasis [67, 68]. Additional file 1: Figure S1. Generation of fluorescent clonal cell lines. Recently, interesting discoveries have been made in the a. Mechanism to induce the expression of fluorescent proteins by field of paracrine signaling. Extracellular vesicles (EVs) have transposon integration. b. Fluorescent proteins codified by transposons. c. Representative images of the clonal cell lines (scale bar = 50 μm) been postulated as a highly efficient method to transform obtained by transposon integration and subcloning: phase contrast and cells. Several groups have demonstrated that EVs upregulate specific fluorescence for every clonal cell line. (PDF 994 kb) prometastatic and tumor angiogenesis pathways [69], even Additional file 2: Figure S2. Phenotypic characterization of MDA-MB- transforming normal cells into cancer cells [70, 71]. Our 231 GFP C3, mKO E10, and Sapphire D7 cell lines. a. Morphological evalu- ation by phase contrast images (scale bar = 200 μm). b. Cumulative popu- complementation assay confirmed intercellular communi- lation doublings per day (CPDs) (relative to MDA-MB-231 CPDs at 7 days cation via factors released from the cells, but we could not in culture) at 7, 14, 21 and 28 days in culture. c-d. Migration capability by distinguish between soluble factors (such as cytokines) or wound healing assay: measurements relative to MDA-MB-231 (c) and rep- resentative images (scale bar = 200 μm) (d). Significant differences were released EVs. Our results showed that exosomes were able determined using ANOVA (Tukey’s multiple comparisons test) (b) and to modify growth rate (Fig. 8e). However, the soluble factors Tukey’s unpaired t-test with Welch’s correction (c). Asterisks indicate sig- seemed to be more efficient than EVs, as the effect of exo- nificant differences when P-values are < 0.05 (*), < 0.01 (**), and < 0.001 (***). (PDF 1639 kb) somes was dose-dependent. Additional file 3: Figure S3. Homing capacity of parental, clonal cell Our results demonstrate how some clones, via soluble lines and an equal mix of all cell lines. Analysis of lung metastasis. factors, EVs and physical interactions, are able to induce in- Hematoxylin/eosin staining. Complete lung reconstruction (scale bar = creased aggressiveness in other clones. Therefore, we can 400 μm). (PDF 823 kb) conclude that the Darwinian clonal progression theory Additional file 4: Table S1. Analysis of mKO mRNA expression vs MDA- MB-231 with the Human Cancer Metastasis Database. (PDF 511 kb) should be complemented with a “cooperative model”, Additional file 5: Figure S4. Clone interactions among parental and where clones are able to interact and transfer “properties” clonal cell lines: Co-culture for 28 days. a. Percentage of total population among one another. represented by clonal cell lines and MDA-MB-231 at 0, 7, 14, 21 and 28 In our study, the biological differences observed between days in co-culture. b-c. Expected (E) percentage of total population calcu- lated using CPDs per day vs. observed (O) population percentage: (b) Co- the different clones have been corroborated by clear differ- culture of two clonal cell lines and (c) co-culture of three cell lines. d-f. Ef- ences in the RNA expression arrays (Fig. 4), which were fect of arsenite on co-culture: (d) Proliferation 72 h after arsenite treat- complemented with methylomic studies (Fig. 5). We have ment (250 μM-90 min). (e-f) Percentage of cell line populations after arsenite treatment: (e) MDA-MB-231 vs GFP C3, (f) MDA-MB-231 vs mKO also demonstrated clonal communication (Figs. 7, 8), and E10. The percentage of each cell line in the total population was de- the challenge now is to identify the factors released by the tected at seeding (day 0), 90 min (day 1) and 72 h (day 4) after treatment. parental cells that enhance survival and metastasis. The g. Co-injection in zebra fish model: percentage of fish with cells in the tail 72 h after injection with individual clones or co-injected with the mix. identification of factors and cytokines that modulate cellular Significant differences were determined using Tukey’s unpaired t-test communication and enhance malignant properties will be with Welch’s correction (b, c, d, e, f). Asterisks indicate significant differ- essential to understand the formation of clonal clusters and ences when P-values are < 0.05 (*), < 0.01 (**), and < 0.001 (***). (PDF 917 kb) prevent metastasis. We propose the term functional hetero- Additional file 6: Figure S5. Clone communication by exosomes. a-c. geneity to highlight the different expression between clones Characterization of exosomes: (a) area, (b) perimeter and (c) diameter. due to extracellular factors and epigenetic changes. This ap- Significant differences were determined using unpaired t-test with proach opens the way to new paradigms in the development Welch’s correction. Asterisks indicate significant differences when P-values are < 0.05 (*), < 0.01 (**), and < 0.001 (***). d. Representative transmission of metastases, including central factors at an intracellular Martín-Pardillos et al. BMC Cancer (2019) 19:666 Page 24 of 26 electron microscopy image of exosomes. e. Immunoblot of exosome Ethics approval and consent to participate markers (TSG101 and CD81) and housekeeping gene (β-actin). (PDF 868 kb) Animal care was handled in accordance with the Guide for the Care and Use of Laboratory Animals of the Vall d’Hebron University Hospital Animal Facility, and the experimental procedures were approved by the Animal Abbreviations Experimentation Ethics Committee at the institution (76/17 CEEA). ANOVA: ANalysis Of VAriance; ATCC: American Type Culture Collection; BCA: BicinChoninic Acid; CD81: Cluster of Differentiation 81; CM: Complete Consent for publication Medium; CPDs: Cumulative Population Doublings; DMBA: 2,4- Not applicable. DiMethoxyBenzAldehyde; DMEM: Dulbecco’s Modified Eagle’s Medium; DNA: DeoxyriboNucleic Acid; DoC: Duct Of Cuvier; DPBS: Dulbecco’s Phosphate-Buffered Saline; dpf: Days Post-Fertilization; eGFP: Enhanced Competing interests Green Fluorescent Protein; Em: Emission; EVs: Extracellular vesicles; The authors declare that they have no competing interests. Ex: Excitation; FACS: Fluorescence-Activated Cell Sorting; FBS: Fetal Bovine Serum; FC: Fold Change; FDR: False Discovery Rate; GEO: Gene Expression Author details 1 Omnibus; GFP: Green Fluorescent Protein; GM-CSF: Granulocyte Macrophage Translational Molecular Pathology Group, Vall d’Hebron Research Institute, Colony-Stimulating Factor; GO: Gene Ontology; HBS: HEPES Buffered Saline; Barcelona, Spain. 2CIBERONC (Centro de Investigación Biomédica en Red de hpi: Hours Post Injection; HRP: HorseRadish Peroxidase; i.m.p.f: Inoculation at Cáncer), Madrid, Spain. 3Cancer Modelling Lab, Roche-CHUS Joint Unit, Mammary Fat Pad; IL: Interleukin; IR: Inverted Repeat; ITRs: Terminal Repeat Santiago de Compostela, Spain. 4Cancer Epigenomics, Translational Medical Sequences; JAK-STAT: JAnus Kinase/Signal Transducers and Activators of Oncology (Oncomet), Health Research Institute of Santiago (IDIS), University Transcription; KEGG: Kyoto Encyclopedia of Genes and Genomes; Clinical Hospital of Santiago (CHUS), Santiago de Compostela, Spain. 5 MCP1: Monocyte Chemoattractant Protein 1; MIP1: Macrophage Hospital Vall d’Hebron, Anatomía Patológica, Barcelona, Spain. 6EMBL/CRG Inflammatory Protein; mKO: Monomeric Kusabira Orange; MTT: 3-[4,5- Systems Biology Research Unit, Centre for Genomic Regulation (CRG), The Dimethylthiazol-25-Yl]-2.5-Diphenyl Tetrazolium Bromide; NCBI: National Institute of Science and Technology, Barcelona, Spain. 7Universitat Pompeu Center for Biotechnology Information; NOD: Nucleotide-binding Fabra (UPF), Barcelona, Spain. 8Institució Catalana de Recerca i Estudis Oligomerization Domain; PBS: Phosphate-Buffered Saline; Avançats (ICREA), Barcelona, Spain. 9Department of Chemical Engineering, PFA: ParaFormAldehyde; PTU: N-Phenylthiourea; PVDF: PolyVinylidene Aragon Institute of Nanoscience (INA), University of Zaragoza, Zaragoza, DiFluoride; rcf: Relative Centrifugal Force; RNA: RiboNucleic Acid; RT- Spain. 10Networking Research Centre on Bioengineering, Biomaterials and PCR: Real-Time Polymerase Chain Reaction; SDS: Sodium Dodecyl Sulfate; Nanomedicine, CIBER-BBN, Madrid, Spain. 11Department of Molecular, Cellular SEM: Standard Error of the Mean; SF: Soluble Factors; STR: Short Tandem and Developmental Neurobiology, Instituto Cajal-CSIC, Madrid, Spain. 12 Repeat; TEM: Transmission Electron Microscopy; TGFR-1: TGF-beta receptor Roche-CHUS Joint Unit, University Clinical Hospital of Santiago (CHUS), type-1; TNBC: Triple Negative Breast Cancer; TPA: 12-O-TetradecanoylPhorbol- Santiago de Compostela, Spain. 13-Acetate; TSG101: Tumor Susceptibility Gene 101; w/v: Weight/volume; WHO: World Health Organization; WNT: Wingless-related iNtegration Site; Received: 28 May 2019 Accepted: 26 June 2019 WT: Wild Type Acknowledgements References Not applicable. 1. Siegel RL, Miller KD, Jemal A. Cancer statistics, 2017. CA Cancer J Clin. 2017; 67(1):7–30. Authors’ contributions 2. DeVita VT, Lawrence TS, Rosenberg SA, DePinho RA, Devita WRA. Hellman, AMP designed and supervised the experiments. AMP and SRYC interpreted and Rosenberg's cancer: principles & practice of oncology. 10th ed. all the results. AMP, SRYC and SH wrote the manuscript. AMP, AVC, GBV, PHB, Philadephia: Lippincott Williams and Wilkins; 2014. RPC, ARC, PJG and EBS carried out the experiments and interpreted the data. 3. Fan C, Oh DS, Wessels L, Weigelt B, Nuyten DS, Nobel AB, et al. AMP and SH interpreted the microarray data. ARC, ADL and RLL analyzed the Concordance among gene-expression-based predictors for breast cancer. N epigenomic experiments and interpreted the results. JC carried out and Engl J Med. 2006;355(6):560–9. interpreted the histopathology analysis. SMV, LS and MLS designed and 4. Gerlinger M, Horswell S, Larkin J, Rowan AJ, Salm MP, Varela I, et al. carried out the transposon sequencing analysis and interpretation of the Genomic architecture and evolution of clear cell renal cell carcinomas data. VS carried out the exosome microscopy preparation and photograph defined by multiregion sequencing. Nat Genet. 2014;46(3):225–33. capture. AB and LLM designed the transposon plasmids, including an 5. Gerlinger M, Rowan AJ, Horswell S, Math M, Larkin J, Endesfelder D, et al. exhaustive study of them, necessary for their use in this manuscript. All Intratumor heterogeneity and branched evolution revealed by multiregion authors read and approved the manuscript and agreed to be personally sequencing. N Engl J Med. 2012;366(10):883–92. accountable for the author’s own contributions and to ensure that questions 6. Nash I. Intratumor heterogeneity and branched evolution. N Engl J Med. related to the accuracy or integrity of any part of the work, even ones in 2012;366(22):2132–3 discussion 3. which the author was not personally involved, are appropriately investigated, 7. Turajlic S, Swanton C. Metastasis as an evolutionary process. Science. 2016; resolved, and the resolution documented in the literature. 352(6282):169–75. 8. Hiley CT, Swanton C. Spatial and temporal cancer evolution: causes and consequences of tumour diversity. Clin Med. 2014;14(Suppl 6):s33–7. Funding 9. Greaves M, Maley CC. Clonal evolution in cancer. Nature. 2012; SRYC acknowledges support from Fondo de Investigaciones Sanitarias (FIS; 481(7381):306–13. PI17/02247 and PI14/01320), Centro de Investigación Biomédica en Red de 10. Lloyd MC, Cunningham JJ, Bui MM, Gillies RJ, Brown JS, Gatenby RA. Cáncer (CIBERONC; CB16/12/00363) and Generalitat de Catalunya (AGAUR; Darwinian dynamics of Intratumoral heterogeneity: not solely random 2017 SGR 1799 and 2014 SGR 1131). AMP is funded by a FP7 Marie mutations but also variable environmental selection forces. Cancer Res. Sklodowska-Curie COFUND program under Grant Agreement n° 267128 2016;76(11):3136–44. (INCOMED program). ADL is funded by a “Juan Rodés” contract (JR17/00016) 11. Heppner GYK, Miller B, Miller F. Tumor heterogeneity in metastasis. Prog from ISCIII. Funding bodies have provided funding to to cover experimental Clin Biol Res. 1986;212:45–59. and publication expenses. 12. Heppner GH, Shekhar M. Tumor heterogeneity is fundamental to the tumor ecosystem. Oncology (Williston Park). 2014;28(9):780–1. Availability of data and materials 13. Gatenby RA, Cunningham JJ, Brown JS. Evolutionary triage governs fitness All microarray data in this publication have been deposited in the NCBI Gene in driver and passenger mutations and suggests targeting never mutations. Expression Omnibus and are accessible through the GEO Series accession Nat Commun. 2014;5:5499. number GSE122008 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc= 14. Alderton GK. Tumour evolution: epigenetic and genetic heterogeneity in GSE122008). metastasis. Nat Rev Cancer. 2017;17(3):141. Martín-Pardillos et al. BMC Cancer (2019) 19:666 Page 25 of 26 15. Almendro V, Kim HJ, Cheng YK, Gonen M, Itzkovitz S, Argani P, et al. Genetic 41. Irizarry RA, Hobbs B, Collin F, Beazer-Barclay YD, Antonellis KJ, Scherf U, et al. and phenotypic diversity in breast tumor metastases. Cancer Res. 2014;74(5): Exploration, normalization, and summaries of high density oligonucleotide 1338–48. array probe level data. Biostatistics. 2003;4(2):249–64. 16. Denisov EV, Litviakov NV, Zavyalova MV, Perelmuter VM, Vtorushin SV, 42. Smyth GK. Linear models and empirical bayes methods for assessing Tsyganov MM, et al. Intratumoral morphological heterogeneity of breast differential expression in microarray experiments. Stat Appl Genet Mol Biol. cancer: neoadjuvant chemotherapy efficiency and multidrug resistance 2004;3:Article3. gene expression. Sci Rep. 2014;4:4709. 43. Sandoval J, Heyn H, Moran S, Serra-Musach J, Pujana MA, Bibikova M, et al. 17. Russo M, Siravegna G, Blaszkowsky LS, Corti G, Crisafulli G, Ahronian LG, Validation of a DNA methylation microarray for 450,000 CpG sites in the et al. Tumor heterogeneity and lesion-specific response to targeted therapy human genome. Epigenetics. 2011;6(6):692–702. in colorectal Cancer. Cancer Discov. 2016;6(2):147–53. 44. Moran S, Arribas C, Esteller M. Validation of a DNA methylation microarray 18. Schneider G, Schmidt-Supprian M, Rad R, Saur D. Tissue-specific for 850,000 CpG sites of the human genome enriched in enhancer tumorigenesis: context matters. Nat Rev Cancer. 2017;17(4):239–53. sequences. Epigenomics. 2016;8(3):389–99. 19. Hanahan D, Weinberg RA. Hallmarks of cancer: the next generation. Cell. 45. Amaro A, Angelini G, Mirisola V, Esposito AI, Reverberi D, Matis S, et al. 2011;144(5):646–74. A highly invasive subpopulation of MDA-MB-231 breast cancer cells 20. Polyak K, Marusyk A. Cancer: clonal cooperation. Nature. 2014; shows accelerated growth, differential chemoresistance, features of 508(7494):52–3. apocrine tumors and reduced tumorigenicity in vivo. Oncotarget. 2016; 21. Ramon Y Cajal S, Capdevila C, Hernandez-Losa J, De Mattos-Arruda L, Ghosh 7(42):68803–20. A, Lorent J, et al. Cancer as an ecomolecular disease and a neoplastic 46. Garcia-Marques J, Lopez-Mascaraque L. Clonal identity determines astrocyte consortium. Biochim Biophys Acta. 2017;1868(2):484–99. cortical heterogeneity. Cereb Cortex. 2013;23(6):1463–72. 22. Cheung KJ, Ewald AJ. A collective route to metastasis: seeding by tumor cell 47. Bribian A, Figueres-Onate M, Martin-Lopez E, Lopez-Mascaraque L. clusters. Science. 2016;352(6282):167–9. Decoding astrocyte heterogeneity: new tools for clonal analysis. 23. Bertolaso M. Philosophy of cancer: a dynamic and relational view. New York: Neuroscience. 2016;323:10–9. Springer; 2016. 48. Grada A, Otero-Vinas M, Prieto-Castrillo F, Obagi Z, Falanga V. Research 24. Chapman A, Fernandez del Ama L, Ferguson J, Kamarashev J, Wellbrock C, techniques made simple: analysis of collective cell migration using the Hurlstone A. Heterogeneous tumor subpopulations cooperate to drive wound healing assay. J Invest Dermatol. 2017;137(2):e11–e6. invasion. Cell Rep. 2014;8(3):688–95. 49. Xie F, Ling L, van Dam H, Zhou F, Zhang L. TGF-beta signaling in cancer 25. Aceto N, Bardia A, Miyamoto DT, Donaldson MC, Wittner BS, Spencer JA, metastasis. Acta Biochim Biophys Sin. 2018;50(1):121–32. et al. Circulating tumor cell clusters are oligoclonal precursors of breast 50. Padua D, Massague J. Roles of TGFbeta in metastasis. Cell Res. 2009; cancer metastasis. Cell. 2014;158(5):1110–22. 19(1):89–102. 26. Zhao Q, Eichten A, Parveen A, Adler C, Huang Y, Wang W, et al. Single- 51. Hutchinson L. Breast cancer: challenges, controversies, breakthroughs. Nat cell transcriptome analyses reveal endothelial cell heterogeneity in Rev Clin Oncol. 2010;7(12):669–70. tumors and changes following antiangiogenic treatment. Cancer Res. 52. Khan GN, Kim EJ, Shin TS, Lee SH. Heterogeneous cell types in single- 2018;78(9):2370–82. cell-derived clones of MCF7 and MDA-MB-231 cells. Anticancer Res. 27. Shen H, Laird PW. Interplay between the cancer genome and epigenome. 2017;37(5):2343–54. Cell. 2013;153(1):38–55. 53. Fisher R, Pusztai L, Swanton C. Cancer heterogeneity: implications for 28. Ramon Y Cajal S, Castellvi J, Hummer S, Peg V, Pelletier J, Sonenberg N. targeted therapeutics. Br J Cancer. 2013;108(3):479–85. Beyond molecular tumor heterogeneity: protein synthesis takes control. 54. Aleskandarany MA, Vandenberghe ME, Marchio C, Ellis IO, Sapino A, Rakha Oncogene. 2018;37(19):2490–501. EA. Tumour heterogeneity of breast Cancer: from morphology to 29. Nguyen A, Yoshida M, Goodarzi H, Tavazoie SF. Highly variable cancer personalised medicine. Pathobiology. 2018;85(1–2):23–34. subpopulations that exhibit enhanced transcriptome variability and 55. Talmadge JE, Fidler IJ. Cancer metastasis is selective or random depending metastatic fitness. Nat Commun. 2016;7:11246. on the parent tumour population. Nature. 1982;297(5867):593–4. 30. Zhou H, Neelakantan D, Ford HL. Clonal cooperativity in heterogenous 56. Talmadge JE, Fidler IJ. AACR centennial series: the biology of cancer cancers. Semin Cell Dev Biol. 2017;64:79–89. metastasis: historical perspective. Cancer Res. 2010;70(14):5649–69. 31. Uhlen M, Zhang C, Lee S, Sjostedt E, Fagerberg L, Bidkhori G, et al. A 57. Poste G, Fidler IJ. The pathogenesis of cancer metastasis. Nature. 1980; pathology atlas of the human cancer transcriptome. Science. 2017; 283(5743):139–46. 357(6352). 58. Cheung KJ, Padmanaban V, Silvestri V, Schipper K, Cohen JD, Fairchild AN, 32. McDonald OG, Li X, Saunders T, Tryggvadottir R, Mentch SJ, Warmoes MO, et al. Polyclonal breast cancer metastases arise from collective dissemination et al. Epigenomic reprogramming during pancreatic cancer progression of keratin 14-expressing tumor cell clusters. Proc Natl Acad Sci U S A. 2016; links anabolic glucose metabolism to distant metastasis. Nat Genet. 2017; 113(7):E854–63. 49(3):367–76. 59. Massague J, Obenauf AC. Metastatic colonization by circulating tumour 33. Egeblad M, Nakasone ES, Werb Z. Tumors as organs: complex tissues that cells. Nature. 2016;529(7586):298–306. interface with the entire organism. Dev Cell. 2010;18(6):884–901. 60. Dotto GP. Multifocal epithelial tumors and field cancerization: stroma as a 34. Ramon Y Cajal S, De Mattos-Arruda L, Sonenberg N, Cortes J, Peg V. The primary determinant. J Clin Invest. 2014;124(4):1446–53. intra-tumor heterogeneity of cell signaling factors in breast cancer: p4E-BP1 61. Curtius K, Wright NA, Graham TA. An evolutionary perspective on field and peIF4E are diffusely expressed and are real potential targets. Clin Transl cancerization. Nat Rev Cancer. 2018;18(1):19–32. Oncol. 2014;16(11):937–41. 62. Michor F, Polyak K. The origins and implications of intratumor 35. Yusa K, Rad R, Takeda J, Bradley A. Generation of transgene-free induced heterogeneity. Cancer Prev Res. 2010;3(11):1361–4. pluripotent mouse stem cells by the piggyBac transposon. Nat Methods. 63. Parsons BL. Multiclonal tumor origin: evidence and implications. Mutat Res 2009;6(5):363–9. Rev Mutat Res. 2018;777:1–18. 36. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina 64. Sottoriva A, Kang H, Ma Z, Graham TA, Salomon MP, Zhao J, et al. A sequence data. Bioinformatics. 2014;30(15):2114–20. big bang model of human colorectal tumor growth. Nat Genet. 2015; 37. Genomes Project C, Auton A, Brooks LD, Durbin RM, Garrison EP, Kang HM, 47(3):209–16. et al. A global reference for human genetic variation. Nature. 2015; 65. Gao R, Davis A, McDonald TO, Sei E, Shi X, Wang Y, et al. Punctuated copy 526(7571):68–74 Epub 2015/10/04. number evolution and clonal stasis in triple-negative breast cancer. Nat 38. Langmead B, Salzberg SL. Fast gapped-read alignment with bowtie 2. Nat Genet. 2016;48(10):1119–30. Methods. 2012;9(4):357–9. 66. Kim C, Gao R, Sei E, Brandt R, Hartman J, Hatschek T, et al. Chemoresistance 39. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The evolution in triple-negative breast Cancer delineated by single-cell sequence alignment/map format and SAMtools. Bioinformatics. 2009; sequencing. Cell. 2018;173(4):879–93 e13. 25(16):2078–9. 67. Schlange T, Matsuda Y, Lienhard S, Huber A, Hynes NE. Autocrine WNT 40. Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing signaling contributes to breast cancer cell proliferation via the canonical genomic features. Bioinformatics. 2010;26(6):841–2. WNT pathway and EGFR transactivation. Breast Cancer Res. 2007;9(5):R63. Martín-Pardillos et al. BMC Cancer (2019) 19:666 Page 26 of 26 68. Lin C, Liao W, Jian Y, Peng Y, Zhang X, Ye L, et al. CGI-99 promotes breast cancer metastasis via autocrine interleukin-6 signaling. Oncogene. 2017; 36(26):3695–705. 69. Feng Q, Zhang C, Lum D, Druso JE, Blank B, Wilson KF, et al. A class of extracellular vesicles from breast cancer cells activates VEGF receptors and tumour angiogenesis. Nat Commun. 2017;8:14450. 70. Jia Y, Chen Y, Wang Q, Jayasinghe U, Luo X, Wei Q, et al. Exosome: emerging biomarker in breast cancer. Oncotarget. 2017;8(25):41717–33. 71. O'Brien K, Rani S, Corcoran C, Wallace R, Hughes L, Friel AM, et al. Exosomes from triple-negative breast cancer cells can transfer phenotypic traits representing their cells of origin to secondary cells. Eur J Cancer. 2013;49(8): 1845–59 Epub 2013/03/05. Publisher’s Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.