Non-Darwinian dynamics in therapy-induced cancer drug resistance.

The development of drug resistance, the prime cause of failure in cancer therapy, is commonly explained by the selection of resistant mutant cancer cells. However, dynamic non-genetic heterogeneity of clonal cell populations continuously produces metastable phenotypic variants (persisters), some of which represent stem-like states that confer resistance. Even without genetic mutations, Darwinian selection can expand these resistant variants, which would explain the invariably rapid emergence of stem-like resistant cells. Here, by using quantitative measurements and modelling, we show that appearance of multidrug resistance in HL60 leukemic cells following treatment with vincristine is not explained by Darwinian selection but by Lamarckian induction. Single-cell longitudinal monitoring confirms the induction of multidrug resistance in individual cells. Associated transcriptome changes indicate a lasting stress response consistent with a drug-induced switch between high-dimensional cancer attractors. Resistance induction correlates with Wnt pathway upregulation and is suppressed by β-catenin knockdown, revealing a new opportunity for early therapeutic intervention against the development of drug resistance.

INTRODUCTION

The development of resistance to chemotherapy and target-selective drugs is the prime cause of therapy failure[1]. Mechanisms implicated range from the canonical “acquisition” of a resistant cell phenotype to whole-tumor level changes involving non-neoplastic stroma[2,3]. Cell-level resistance is still considered the major cause of loss of drug-sensitivity because it can readily be reproduced in cell cultures and explained by the nature of associated molecular changes, such as expression of detoxification proteins[4], structural alteration of the drug target protein[5], or activation of alternative growth and survival pathways [6-8]. One of the best characterised molecular mechanisms of cellular resistance is the expression of ABC-transporter family proteins, such as MDR1, which pumps a variety of drugs out of the cell[9], thereby conferring the multi-drug resistance (MDR) phenotype.

But how are these molecular changes acquired? Current paradigm holds that treatment selects for cancer cells carrying a random genetic mutation that preexists before treatment and happens to confer a survival advantage in the presence of the drug, resulting in the “emergence” of drug-resistant clones[10,11]. Point mutations in drug target proteins that alter the drug docking site [5] or genomic rearrangements in the MDR1 regulatory region[12,13] that cause high expression of the MDR1 protein support the theory of a somatic evolution of cancer cells that follows the Darwinian scheme of mutation and selection. However, the increasing realization that cancer cells exhibit a rich intra-clonal dynamics manifested as “non-genetic heterogeneity” complicates the picture[14]. Such phenotype heterogeneity in the absence of genetic variation is the combined consequence of multi-stability in gene expression dynamics (the coexistence of multiple stable steady-states, or “attractors”, in gene regulatory networks[15-17] and of gene expression noise [18]. Thus, one genome can produce multiple enduring (stable) or transient (metastable) phenotypic states. This departure from the simple one-to-one mapping between genotype and phenotype necessitates a reexamination of the standard scheme of somatic evolution driven by random genetic mutations[19,20].

Non-genetic cell state dynamics (that manifest in the behavior of tumor cells) have recently received renewed interest and are best understood in terms of stochastic multi-attractor dynamics: Tumor cells within a clonal population spontaneously switch between several (meta)stable attractor states which represent different developmental states, including mesenchymal, epithelial as well as cancer stem cell-like states[21-26]. Cells in the latter state are naturally endowed with increased xenobiotics resistance [22,27]. In an unperturbed cell population, multi-stability is manifest as broad quasi-continuous or as multi-modal distribution of a phenotypic marker across the entire population[15]. Switching between attractor states can occur in two ways: in a spontaneous and stochastic fashion[2] due to noise-induced attractor state transition[28,29] and in a directed way following a perturbation by external signals that alter gene expression. Both have consequences for resistance development[30].

First, stochastic non-genetic phenotype switching can act as a source of random variability – the substrate for Darwinian selection[31,32]: Cells that by chance occupy states that are more resilient to cytotoxic stress, including therapy-induced cytotoxicity, can be transiently selected for during treatment. Because these cell states are non-genetically inherited over many cell generations they can in principle promote evolution according to the Darwinian scheme – in the absence of mutations[29,33]. The transient selection of cells in the resistant state allows a subpopulation of temporarily “fitter” cells to expand, thereby increasing the probability for adaptive genetic mutations. Such “mutation-less pre-selection” could accelerate classical Darwinian evolution of drug resistance[31-33] as observed for antibiotic resistance in microorganisms (“persisters”)[34-36].

Second, induction of attractor state switching by external signals opens the possibility for a Lamarckian scheme of evolution[37]: A perturbation by cytotoxic agents may “instruct” the cell to enter an attractor state that confers the stem-like, more stress-resistant phenotype –perhaps recapitulating a generic, physiological stress-response– which can be passed on to subsequent cell generations. The non-genetic inheritance of an acquired adaptive trait at the cellular, not organismal level, does not violate the neo-Darwinian dogma [37-39]. In fact, chemotherapy and irradiation appear to “cause” the emergence of resistant, stem-like or mesenchymal cancer cells [26,40-48].

Promoter analyses have shown that chemotherapy leads to changes in DNA methylation and histone modification in the MDR1 gene locus [49]. Such findings are compatible with both schemes, Darwinian selection and Lamarckian instruction. However, a distinction is rarely explicitly articulated. While clinicians often take the apparent “induction” of resistance markers following treatment for granted, even viewing it as a form of “active” adaptation because of their rapid and nearly inevitable occurrence, biological orthodoxy assumes by default a Darwinian selection[11,21]. The rapid appearance of stemness markers or MDR1 expression following treatment has typically been assessed at the level of tumor tissues or whole cell populations[43]. Thus, it remains open to what extent the increase of MDR1 expression after treatment is the result of very rapid selection of cells already residing in a state with an active MDR1 locus (Darwinian scheme) or of cell-autonomous gene induction in individual cells (Lamarckian scheme).

To understand the elementary dynamics of resistance development here we determine the relative contribution of these two (non-genetic) schemes of emergence of the MDR phenotype in HL60 acute myeloid leukemia cells, which have long served as a model for MDR1-dependent drug resistance[50]. We show by quantitative measurement and modeling that appearance of MDR1-positive cells 1-2 days after treatment with vincristine is predominantly mediated by cell-individual induction of MDR1 expression and not by selection of MDR1 expressing cells. We confirm this by single-cell longitudinal monitoring. Drug-induced resistance and MDR1 expression correlated with up-regulation of Wnt-signaling pathway and could be suppressed by knock-down of β-catenin. Following transient low-dose chemotherapy, surviving cells exhibited a persistent transcriptome change indicative of a lasting stress-response state, consistent with switching between high-dimensional cancer attractors [51]. Acknowledging that resistance can be promoted by a non-genetic Lamarckian mechanism opens a new window for pharmacological interference with resistance.

RESULTS

Spontaneous non-genetic drug-resistant state in tumor cells

We observed that within a clonally derived population of cultured leukemia cells (HL60 cell line) a small subpopulation (~1-2%) of cells consistently expresses high levels of multidrug resistance protein 1 (MDR1High) on its cell surface and exhibits the MDR-phenotype as measured by the fluorescent dye efflux assay (=effluxHigh) in the absence of drugs (). The MDR1Low and MDR1High subpopulations also differed greatly in their sensitivity to killing by the chemotherapeutic agent vincristine (VINC) [52] (). Both sub-populations correspond to metastable epigenetic states [15,53] (Supplementary ) because they both re-established the initial population distribution after isolation by fluorescence activated cell sorting (FACS) (). The MDR1Low cells accomplished the repopulation of the original distribution within <1 day while the MDR1High cells took ~18 days to do so (Supplementary ). Similar relaxation dynamics was observed when the functional MDR phenotype was measured using the efflux of fluorescent dyes.

When cells were treated with low-dose VINC (10nM), within 48hrs the MDR1High subpopulation increased from <2% to >25% in a dose-dependent manner (Supplementary ). Appearance of effluxHigh (MDR) cells was even more rapid and pronounced, typically reaching 30-40% within 24-48h (. A similar MDR response could be induced with another drug, doxorubicin (Supplementary ). This rapid response, in line with prior biochemical analyses[54], raises the question whether this population shift was driven by a cell-individual induction of the MDR phenotype or by (non-genetic) selection of the pre-existing “epigenetic” MDR1High cells in the stationary populations because of their survival advantage in the presence of VINC (). In the absence of the drug, the sorted cells in the MDR1Low/effluxLow state had a ~5-fold net growth advantage over cells in the MDR1High/effluxHigh state (Supplementary ) and accordingly, the DNA content of live cells measurement revealed a smaller fraction of the MDR1High/effluxHigh cells in the S/G2 state of the cell cycle (Supplementary ). This difference in cell cycle status disappeared when the two subpopulations re-equilibrated to the same MDR1 distribution (Supplementary ). In the presence of VINC (10nM) the relative fitness was reversed: the effluxLow cells exhibited reduced growth and were growth-arrested after three days while the effluxHigh cells survived, displaying slow net population growth (Supplementary ).

Quantitative model of cell state interconversion

The current paradigm in tumor biology assumes a predominantly selection-based mechanism to explain population-level shifts of phenotypes, even for acquisition of epigenetic states[21,55]. This process requires that the spontaneous presence of 1-2% cells in the MDR1High/effluxHigh state () exploit their growth advantage in the presence of VINC (Supplementary ). To quantify the theoretical contribution of selection vs. instruction (), we analyzed the observed relative growth and induction kinetics in two mathematical models. First, we used a simple kinetic state transition Markov model to examine the contribution of growth rate and transition rate in a population at steady state as often employed for modeling tumor cell dynamics[21] (Methods). With the observed numbers for the ratio of the observed effective growth rate constants (g, g) (Table 1), we found that a stationary state with a proportion of 1.5% effluxHigh cells, as observed, does not exist unless one accepts substantial state transitions between the effluxHigh and effluxLow states (). Consequently, already under basal conditions (absence of drugs) the measured difference in the effective growth rates of the two subpopulations in isolation is far from being able to ensure the preservation of the observed steady proportions of these two fractions.

Table 1

Growth rate constants.

	Effective Growth Rates, g
Control, no drug	gL = 0.50/d	gH = 0.10/d
+ Vincristine	g’L = 0.25/d	g’H = 0.37/d

Measured net growth rate constants for the effluxLow and effluxHigh subpopulations in the presence (g, g) and absence (g’, g’) of vincristine 10 nM.

In a second model, we examined the response to treatment, where stationarity/equilibrium is not reached within the observed time period. We used ordinary differential equations (ODEs) to model the joint effect of differential growth and transitions between the two states on the cell population dynamics (Methods): where x and x denote the population fraction of effluxLow (L) and effluxHigh (H) cells and g and k are the respective rate constants for effective growth and transition, and are separately measured in the absence (g, k) and presence (g’, k’) of the drug (table in Fig. 1C). Because the observed change after 24h is far from the new equilibrium state and we are only interested in the relative contribution of state transition to the departure from the existing stationary population proportions, we assumed a simple first-order transition, kx. This equation unites the Darwinian and Lamarckian principles because the effect of selection will come from the difference in g and induction is captured by the difference in k for the two states, in the presence vs. absence of the drug ().

Figure 1

Dynamical heterogeneity of MDR1 expression within a clonal population of HL60 cells

( A distinct subpopulation of 1-2% of the cells of a clonally derived HL60 cell population consistently expresses high levels of MDR1 on the cell surface in the absence of drugs exposure. The MDR1High (red) and MDR1Low sub-populations (blue) differ in sensitivity to vincristine after 48h ( Measurements of population dynamics and effective growth were obtained in three different laboratories using different culture of HL60 cells and representative results are shown. Error bar, standard deviation of one representative experiment with (n=2) biological replicates Scheme of the state transition model for distinguishing between drug-induced shifts in state transition rates (cell-individual switch to the MDR phenotype) versus drug-induced growth rate differences (selection of the MDR phenotype) x, population fraction of cells in the respective state indicated by the index: H, MDRHigh (=effluxHigh) and L, MDRLow (=effluxLow). k, kinetic rate constant for the first-order state transition represented by the arrows. P state transition probability used in the Markov model. ( Results of the steady-state Markov model. The state transition and “self-renewal” probabilities required to reach the steady state, shown as heat map with colors indicating the steady-state ratio x (color bar) as a function of the ratios of the Markov model probabilities P (see METHODS). Change in ratio of transition probabilities P (vertical axis) visibly affects x while change in the ratio P/P does not result in significant change of x. Undefined regions are marked by *. ( Results of the non-equilibrium ODE model. Color map represents the parameter space indicating which combination of the two sets of parameters, the ratio of the relative growth rate constants, g/g (horizontal axis), and the ratio of the state transition rate constants, k/k (vertical axis), causes which population fraction x/x (color map) 24h after addition of the VINC.

Instead of fitting the unknown parameters k and k we computed, based on the measured numbers of differential growth rates, how a change of the values for k and g due to presence of the drug would account for the observed ratio of the two subpopulations, effluxLow and effluxHigh, at 24h after administration of VINC, . . The parameter plane () displays the non-steady state population ratio, r(24h) = x/x after 24h treatment with VINC, as a function (color of the map) of the ratios of the growth and state transition rate constants. The almost horizontal course of the color contour lines, parallel to the x-axis that represents variation of the growth parameters g/g, indicates that a shift of x/x (color) is minimally affected by change of the relative growth rates change but instead is predominantly defined by a change in the relative state transition rates. Clearly, to achieve the observed appearance of a fraction of 30-40% effluxHigh cells after 24h (), corresponding to a ratio x/x ≈ 0.5-0.7 (= green zone in parameter space in ), the measured growth advantage of the effluxHigh cells in the presence of VINC, at g’/g’ ≈ 0.25/0.37 = 0.67, is far from sufficient (dotted vertical line in the parameter space of ). If there were no cell-individual state transitions, then, with the observed growth differential g’/g’ (), selection alone could account for only an increase of MDRHigh cells to x/x = 0.04 after one day (corresponding to a population fraction of MDRHigh of ~ 4%) instead of the observed x/x = 0.67 (=40% MDRHigh).

The rapid appearance of hnRNA for MDR1 following a 24h-pulse of VINC by targeting the RT-PCR to the first intron-exon junction, with a >20-fold induction of MDR1 pre-mRNA at the whole-population level within 30 min of VINC treatment (, followed by detectable expression of mature mRNA followed within 24h (Supplementary ) supports an induction by a molecular change. However this finding does not prove induction because it could in principle reflect an extreme selection of “fitter cells” that display an intrinsic high constitutive synthesis of the MDR1 transcript.

Validation of cell-individual induction of resistant state

Unequivocal demonstration of cell-individual induction (“instruction”) of the MDR phenotype requires the direct observation of the actual induction event in the very same cell before and after addition of the drug to the medium by real-time longitudinal monitoring of the cell culture during treatment. The drug-treated cells preloaded with fluorescent dye (=effluxLow) displayed a visible reduction of fluorescence starting 12h after addition of the drug. In contrast, no change in fluorescence was detectable in the untreated cells. We also observed onset of apoptosis as indicated by DNA condensation in the VINC treated sample after >24h ( Supplementary ). Counting after a typical 48h longitudinal monitoring revealed 63% of the live cells treated with VINC exhibited elimination of the dye, representing the switch to the effluxHigh phenotype compared to 16% of untreated cells (n=80 cells) ().

To demonstrate that selection per se plays no significant role in the emergence of cells with the MDR phenotype, we decoupled MDR1's functional activity from its expression by blocking MDR1 mediated drug efflux with verapamil [56,57]. By separating fitness function of a trait from expression of that trait we can expose the role of instructive (non-selective) factors of phenotype change. We first confirmed that baseline or VINC-induced efflux of the fluorescent dye calceinAM was reduced almost completely by saturating doses of verapamil (Supplementary ). However, this inhibition of efflux did not measurably suppress the early emergence of the subpopulation of MDR1 expressing cell after VINC treatment (). The observation that cells increased MDR1 expression independently of the increase of the length of pre-incubation with verapamil over the time window in which verapamil gradually unfolds its inhibitory function further suggests that MDR1 induction was independent of the pump function of this protein (Supplementary ). This corroborates the role of a selection-independent, instructive mechanism, at least for the rapid appearance of this new phenotype after chemotherapy.

This persistent effluxHigh state induced by a transient drug treatment was also associated with improved survival when the same cells were re-exposed to VINC after washout of the drug (). Thus, drug-induced resistance is non-genetically inherited across cell generations independently of the presence of the drug – at least for a limited period of time (Supplementary ). This dynamics represents a Lamarckian scheme for acquisition of a new phenotype.

Wnt-signaling mediates state transition into a stress state

Since the presence of discrete effluxHigh and effluxLow subpopulations may reflect transitions between distinct stable cellular states (attractor states) we next measured their transcriptomes after a 24h pulse of 10 nM VINC (), when the cell population exhibited a stable bimodal distribution. Even after this short time, globally distinct gene expression pattern were seen when the effluxHigh and effluxLow subpopulations were compared (). Comparison revealed 974 significantly differentially expressed genes in these two subpopulations, indicating a fractional, globally diverse response within a clonal cell population. Gene-Ontology (GO) analysis unveiled the enrichment in this set for genes involved in cell cycle, translation, ribosome and rRNA synthesis as well as response to DNA damage, metal-binding, oxidative phosphorylation and mitochondrial function (), suggesting that these two transcriptomes represented biologically distinct, high-dimensional attractor states [17].

To determine if a stem-cell like state has actually been induced by VINC treatment [26,40-48]., we next performed pairwise comparisons for all sorted subpopulations, now including untreated cells: cells treated and sorted for effluxHigh, cells treated and sorted for effluxLow, as well as untreated mock-sorted and treated mock-sorted cells. The set of 2096 genes that were significantly expressed above background (p<0.05, Methods), and whose expression level differed for each comparison pair by more than 4-fold, were first manually examined for the relevant functional annotations using the NCBI Gene database and the stem-cell signatures reported in Brandenberger et al. [58]. Of note are the alterations in the expression of genes that belong to the Wnt and Polycomb pathways, consistent with the role of a stemness signature in drug-resistant tumor cells. These differentially expressed genes were also subjected to unbiased gene set enrichment analysis (GSEA) () which also extracted the Wnt signaling [59] gene set (Supplementary ), in line with previous studies [60,61] showing an apparent “induction” of Wnt during therapy. Many of the gene expression changes induced by VINC treatment were not detected when we simply compared treated and untreated whole (not sorted) cell cultures, highlighting the cell population heterogeneity and the importance of cell sorting to isolate relevant cell subpopulations for biochemical cell analysis.

A feature of a high-dimensional attractor state is the memory of the perturbation, i.e. a lasting change of a large set of responding genes that persists after removal of the perturbation. To determine if the transcriptome-wide adaptive response exhibited such memory, VINC was washed out after the 24h treatment. Whereas the effluxHigh phenotype persisted for a week after transient exposure to VINC (), the VINC-treated cells remained globally altered beyond 17 days, long after the population had relaxed to the native distribution with baseline efflux ( ). Crucially, if the reappearance of the effluxLow state, which has a growth advantage in the absence of drug (), is the result of selection of naïve effluxLow cells (which either have never responded or fully reverted back) one would expect to see reappearance of the transcriptome of untreated cells. However, in these post-treatment cells, the global changes in gene expression triggered by VINC persisted for > 17 days, suggesting that VINC induced an adaptive, slowly reversible response with respect to the efflux phenotype as a non-specific stress-response, but a stable (apparently irreversible) attractor transition with respect to other state space dimensions, which are orthogonal to those conferring the efflux phenotype.

Since many Wnt-pathway components () were highly induced in the effluxHigh cells we next knocked-down β-catenin in HL60 cells to determine whether it plays an active role in the induction of the MDR phenotype (). When β-catenin was suppressed, 60h of VINC treatment reduced both efflux capacity () and MDR1 expression by half (). This inhibition had functional consequences because the viability of β-catenin-knockdown cells in just 1 nM VINC was reduced to half compared to the control (). Rescue by ectopic overexpression of an RNAi-resistant β-catenin construct completely restored VINC induced efflux and expression of MDR1 in the knock-down cells, confirming specificity of the knockdown (Supplementary ).

DISCUSSION

Our analysis of the rapid “appearance” of the MDR phenotype and of MDR1 expression following chemotherapy provides evidence that this early drug resistance phenotype can be induced by a Lamarckian instruction, independent of selection. Our conclusion is supported by a series of findings: First, the experimental measurement of a rapid drug-triggered induction of a MDR1High (effluxHigh) subpopulation comprising 40% of the cells within one day and associated mathematical cell population dynamics modeling show that the observed moderate fitness advantage of the effluxHigh cells cannot account for this response kinetics. Second, single-cell longitudinal monitoring directly demonstrates “true”, cell-individual adaptation. Third, the decoupling of the expression of a trait from its function (which conveys selective advantage) excludes a role for Darwinian selection. Finally, transcriptome analysis reveals a genome-wide, distinct and enduring stress-induced state that is unlikely orchestrated by random mutations.

The results of course do not preclude a role of canonical Darwinian somatic evolution driven by selection of random genetic mutants at later stages, as amply supported by the observed genomic alterations whose nature readily offer a mechanistic rationale for selective advantage. However, the argumentation for this evolutionary scheme as the sole mechanism requires the assumption of a substantial amount of preexisting mutations[11]. This requirement is alleviated by admitting non-genetic processes as catalyzer[31,32] which is mediated by cells that either transiently, and by chance, occupy a stem-like attractor states and hence survive the treatment, or are induced by the cytotoxic stress to enter such a protective state. Here we show that the latter dominates.

Although the apparent “activation” of resistance mechanisms and alternative survival pathways or of stem-like states after drug or radiation therapy is frequently observed[26,40-48] the distinction between Lamarckian induction and Darwinian selection is rarely explicitly articulated. Existing thinking in cancer biology tacitly implies the latter but often communication is blurred by the use of metaphoric shorthand expressions that are common in evolution biology, such as the “the tumor adapts to the therapy”, which suggests the former. Awareness of this dualism and specifically, of non-genetic dynamics may help to explain several non-intuitive tumor behaviors, such as: why does treatment not only cause drug resistance but inseparably also increases malignancy in recurrent tumors; why can drug resistance not be suppressed by just blocking MDR1[57]; why do early tumor cell clones disappear and reappear[62] or why do recurrent tumors, after developing resistance to target-selective drugs, become sensitive to the same therapy again[63,64]. Considering non-genetic and drug-triggered cell state dynamics may open new opportunities for the management of resistant tumors, such as molecular pathways that can be targeted prior to conventional treatment to prevent therapy-induced tumor progression.

METHODS

Cell culture

Acute leukemic cell line HL60 was obtained from ATCC and independently re-cloned twice from individual cells and cultured in three independent laboratories (see author affiliations). HL60 cells were cultured in Iscove's Modified Dulbecco's Medium (IMDM, Invitrogen) supplemented with 20% fetal bovine serum (FBS, Sigma), 1% L-glutamine, penicillin (100 U/mL, Invitrogen) and streptomycin (100 mg/mL, Invitrogen). Cell number was monitored daily and culture was maintained at a density of 2×105 – 2×106 cells/mL.

Viable Cell Count

To determine the viability and number of cells, 0.4% trypan blue solution was used. Cell suspensions were diluted 1:5 with trypan blue and viable cells, that exclude the dye, were scored on a hemocytometer under a light microscope.

Flow cytometry

HL60 cells were labeled with MDR1/P-glycoprotein APC-conjugated mouse anti-human monoclonal antibody (e-Biosciences, 2.5ng/μL). Flow cytometry analysis was performed on a BD FACSCalibur cell cytometer or a Guava cell cytometer. To determine the percentage of labeled cells, a quadrant gate using an equal concentration of a relevant mouse isotype control was placed. For each analysis 10,000 of viable events that exclude PI were saved. For cell cycle studies, cells were labeled with 1 mM DRAQ5 (Axxora, San Diego, CA) according to manufacturer's instructions.

MDR Functional Assay

HL60 cells (2×105) were washed in HBSS/ 5% FBS (washing buffer) and then incubated with 1nM of calcein AM (CaAM, Invitrogen) for 15 minutes at 37°C. Cells were then washed in cold buffer and re-suspend in PI staining buffer (PI in washing buffer, 1:200 dilution) and samples were kept on ice until analysis. For each analysis in the cytometer, 10,000 of viable events that exclude PI were saved. As controls, cells that had not been loaded with CaAM were used. In addition, dead cells were gated out also using scatter characteristics.

FACS

Fluorescence-activated cell sorting was conducted on a BD Biosciences Aria II at the Hematologic Neoplasia Flow Cytometry Facility of the Dana Farber Cancer Institute. For studies of the dynamics of sorted subpopulations, antibodies (2.5ng/μL) were removed following cell sorting using brief incubation in low-pH buffer [29].

RNA isolation

~7-10×106 cells were collected for RNA isolation from each condition or subpopulation. RNA was extracted following standard protocol from RNeasy Mini Kit (Qiagen).

qRT-PCR

200ng of total cellular RNA was used to prepare cDNA following the standard protocol from iScript cDNA synthesis kit (BioRad). qPCR was performed according to the manufacturers protocol using iTaq SYBR Green Supermix with Rox (BioRad). 300nM each of forward and reverse MDR1 primers were used. MDR1 primer sequences were as follows: For the hnRNA(Ex1-Intr1): 5’-CTCACTTCAGGAAGCAACCA-3’ (forward) and 5’-TGATTGCAAACTTCTAGTCAAGACA-3’ (reverse); for the hnRNA(Intr1-Ex2): 5’-TGGAGAGGTCGGAGTTTTTG-3’ (forward) and 5’-GGTTGAATTTCCAGGAGGAATG -3’ (reverse); for the coding region: 5’-TACAGTGGAATTGGTGCTGGG-3’ (forward) and 5’CCCAGTGAAAAAATGTTGCCA-3’ (reverse). qRT-PCR was carried out with a BioRad CFX96 real time system C1000 and the iQ5 thermal cycler.

Single cell imaging

HL60 cells were labeled with CaAM using the Vybrant Multidrug Resistance Assay or 10μg/mL of Rhodamine123 (Rh123, Invitrogen). Single cells were seeded into 96-well optical plates with a BD FACS Aria (Dana Farber Hematologic Neoplasia Core facility) based on size/forward scatter, regardless of CaAM fluorescence. Following plating, a baseline fluorescence image was obtained for each individual cell. 10nM vincristine or an equivalent volume of drug-free control media was added to each well and the dye pumping function of individual cells was measured by imaging CaAM fluorescence 24h and 48h after drug administration. Phase contrast imaging was conducted in parallel to assist in locating the cell within the well.

Live cell microscopy

HL60 cells at a density of 1×106 cells/mL were seeded for 30min in a glass-bottom dish (Iwaki) and imaged with a Zeiss Exciter (Plan-Neofluar 40×NA 1.3 Oil DIC). Cells were stained with Rhodamine123 (Invitrogen, 1:20,000) or Hoechst 33342 (Invitrogen, 1:200,000). For both dyes, expression level was monitored until it reached a plateau (uptake) before initiation of experiment. Image analysis was conducted with LSM Aim Software.

β-catenin constructs

β-catenin small hairpin construct: β-catenin targeting shRNA was kindly provided by Jürgen Deka. It was cloned using a published β-catenin target sequence [65]. The oligonucleotide was subcloned (Xho1 and EcoR1) into the pInducer11 backbone (kindly provided by Stephen Elledge [66]) to result in pInd-β-catenin.

shRNA resistant β-catenin: a shRNA resistant β-catenin cDNA was kindly provided by Jürgen Deka and Frédérique Baruthio. It is based on a β-catenin D164A plasmid [67] which was kindly provided by Konrad Basler. The D164A mutation was reversed by PCR. The pInd-β-catenin construct targets the following sequence (underlined) in β-catenin: (1580bp)GTCTGCCAATTGTGCGTACTGTCCTTCGGGCT (1630bp). Wobble base mutations were introduced under consideration of highest human codon usage into the target sequence (underlined), resulting in a rescue cDNA fragment: GTGGGGCATGAGCC. The β-catenin rescue fragment was amplified by PCR using AttB containing primers and subsequently subcloned into pInducer20 (kindly provided by Stephen Elledge [66]) resulting in pInd-β-cat-rescue.

Western blot analysis

HL60 parental cells, sh-β-Catenin knockdown cells and rescued (sh-resistant cDNA) knockdown cells were cultured in regular growing medium supplemented with 2ug/ml Doxycycline for 96 hours. 107 cells were suspended in CelLytic MT cell lysis reagent (Sigma-Aldrich) containing protease inhibitors (Complete mini cocktail, Roche). The lysis solution was sonicated and centrifuged and the total protein extraction contained in the supernatant was quantified based on the Bradford colorimetric assay (Bio-Rad). Samples were resolved by electrophoresis on a 8% SDS-PAGE gel, transferred to a PVDF membrane and probed with primary antibodies against β-Catenin (1:2000 in TBST +5%BSA, BD Biosciences) and α-Tubulin (loading control, 1:3000 in TBST, Calbiochem). A horseradish peroxidase-conjugated goat anti-mouse antibody was added, and secondary antibodies were detected through chemiluminescence (ECL, Amersham) on a Fusion FX7 CCD image acquisition system.

Statistical analysis

Two-way ANOVA followed by Dunnett's multiple comparisons test was performed using GraphPad Prism version 6.00 for MacOS, GraphPad Software, La Jolla California USA, www.graphpad.com.

Mathematical modelling

For the Markov state transition model we make two assumptions. First, we discretize the continuous distribution of cell population as two discrete states, MDRLow and MDRHigh. This is warranted given the bimodal distribution (that is better visible in the drug induced conditions). Second, we assume that the probabilities at each time-step for staying in the same state (“self-renewal”) and for transition are linear (first-order kinetics). The governing equations of the Markov model of the state transition model are (Eq.1): where , denote the initial relative proportions of cells in the two respective states MDRLow (L) and MDRLow (L) and MDRHigh (H) and , represent the relative proportions of cells in the MDRLow and MDRHigh states in the n-the time step (x + x = 1). The P represents the probabilities of transitioning from state i to j (the basis of instruction by the drug) and P represents the “self-renewal” (probability of staying in the same state), as previously employed [21] although the Markov mode is not a growth model. Because of its Markov property, the transition matrix P satisfies: It also has a maximum eigenvalue λ = 1 and a corresponding eigenvector, which determines the final steady state of cell population. With this equation, we can study the influence of the transition rates P,P and the self-renewal rates P,P on the steady-state. From Eq.(3) we obtain for the steady-state ratio r* = x*/x* of the two populations, effluxLow and effluxHigh: With this model we determine the parameter space which maps the relative values of growth and transition rates to a given steady-state composition, r*. Using the numbers for the observed effective growth rate constants (g, g) () in Eq.(4) for the relative magnitudes of P and P, we display r* as a function of the ratios of the transition rates P and P and of the self-renewal rates P and P on the steady-state ratio r* (). One cannot measure the effective transition rates accurately because the observed apparent reconstitution of the effluxHigh population from the sorted effluxLow cells is a composite result of effective transition and effective growth and the observed state transition is slower than cell proliferation. By contrast, the initial effective growth rates of the sorted subpopulations can be reliably estimated because in this initial period the transition to and from the other state is negligible for both subpopulations due to the small probabilities P or the small population size n.

For the Linear ordinary differential equation model as with the Markov model, we discretize the continuous distribution of cell population as two distinct states, MDRLow and MDRHigh, and assume again that the rate of both cell effective growth (death and birth combined) and transitions follow first-order kinetics. The governing equations of the linear ODEs model can be written as follows (Eq.5): where x and x denote the population fraction of effluxLow (L) and effluxHigh (H) cells as in the last Markov model. g and k are the respective rate constants for effective growth and transition and are separately measured for absence (g, k) and presence of the drug (g’, k’). The mathematical solutions of Eq. (5) depend on the eigenvalue of following matrix:

With the quadratic equation for the eigenvalue λ and the solution λ1 and λ2, the solution of ODEs is: We numerically estimated the ratios of growth rate constants constants g/g and the transition rate constants k/k that produce a given ratio of the subpopulations r(24h)=x(24h)/x(24h) after 24h of drug treatment as displayed in using the values shown in the .

Microarray analysis

Microarray studies were performed by the Molecular Genetics Core Facility at Children's Hospital Boston supported by NIH-P50-NS40828 and NIH-P30-HD18655. RNA samples were hybridized to Illumina Human Ref-3 expression BeadChips and initial processing and normalization was performed with BeadStudio software. BeadChip internal p-values (technical bead replicates) were used for filtering out genes significantly expressed above background noise. To filter out genes with signals not significant above background noise, p-value from Illumina BeadStudio (see above) of 0.05 was used as cutoff and only genes with a p-value smaller than 0.05 in all four samples passed the filter. From the original set of 18,401 probes, 2096 genes met this criterion. Significant genes were identified with Significance Analysis of Microarrays (SAM) using a FDR of <1% on all genes, resulting in 974 genes. Self-organized maps of significant gene lists were prepared using the Gene Expression Dynamics Inspector software (GEDI, http://www.childrenshospital.org/research/ingber/GEDI/gedihome.htm)[68]. Hierarchical clustering based on Euclidean distance was performed on time course gene expression data normalized by z-score. Gene expression data have been deposited in the GEO databank under this publication ID. For Gene Set Enrichment analysis (GSEA). the dataset analyzed corresponds to the four samples shown in (MDR- High subpopulation after vincristine (sample 1); MDR-Low subpopulation after vincristine (sample 2); no treatment/mock sort, (sample 3); and vincristine treatment/mock sort (sample 4). The gene set with 2096 genes (see above) was normalized using the quantilenorm function in MatLab®. The aim was to identify genes differentially expressed between the untreated population (sample 3) and the other three samples. To do that, we calculated the pair-wise ratios between sample 3 and samples 1, 2 and 4, respectively. Genes with a log fold change (log2(ratio)) higher or equal than 2 or smaller or equal than −2 are listed in Supplementary Data 1 The genes in Supplementary Data 1, columns 2 (sample 3 vs. sample 2) and column 3 (sample 3 vs. sample 4) were manually checked for their biological function using NCBI – Gene database. For the genes listed in Table 1 column 1 (sample 3 vs. sample 1) we manually mapped it to the stemness pathways identified by Branderberg et al [58].

The set of 2096 genes was subjected to GSEA (Gene Set Enrichment Analysis) [69] to identify gene sets in MSigDB enriched for genes of KEGG pathways in all the samples when using untreated mock-sorted samples as reference (Supplementray ). A ranked list is presented as a heat map plot in Supplementary showing the most significantly differentially expressed genes between mock sorted and the other 3 samples (treated/mock-sorted, treated/sorted low efflux, treated/sorted high efflux). The GSEA analysis expands the manual approach and confirms the stem-cell signature of Wnt-pathway, whose genes are differentially expressed after treatment.

For Gene Ontology analysis. the same gene set used for GSEA was tested in DAVID (Database for Annotation, Visualization and Integrated Discovery) against the Homo sapiens gene reference set [70]. Sixteen functional categories were found to be significantly enriched in the dataset with a p-value less than 0.001. These GO terms, their significance and the number of associated genes in the dataset are summarized in Supplementary Table S1.