Heterogeneous Mechanisms of Secondary Resistance and Clonal Selection in Sarcoma during Treatment with Nutlin

Nutlin inhibits TP53-MDM2 interaction and is under investigation in soft-tissue sarcomas (STS) and other malignancies. Molecular mechanisms of secondary resistance to nutlin in STS are unknown. We performed whole-transcriptome sequencing (RNA-seq) on three pretreatment and secondary resistant STS cell lines selected based on their high primary sensitivity to nutlin. Our data identified a subset of cancer gene mutations and ploidy variations that were positively selected following treatment, including TP53 mutations in 2 out of 3 resistant cell lines. Further, secondary resistance to nutlin was associated with deregulation of apoptosis-related genes and marked productive autophagy, the inhibition of which resulted in significant restoration of nutlin-induced cell death. Collectively, our findings argue that secondary resistance to nutlin in STS involved heterogeneous mechanisms resulting from clonal evolution and several biological pathways. Alternative dosing regimens and combination with other targeted agents are needed to achieve successful development of nutlin in the clinical setting.


Introduction
The tumor suppressor TP53 plays a crucial role in protection from malignant tumor development. It is a transcription factor which is activated following stress and regulates multiple downstream genes involved in cell cycle control, apoptosis, DNA repair, and senescence [1]. In non-stressed cells, the level of T53 is controlled tightly by MDM2 (murine double minute 2). MDM2 regulates p53 through a negative-feedback loop. When the nuclear TP53 level is elevated, it activates the transcription of the MDM2 gene. In turn, MDM2 binds to TP53 and blocks its transactivation domain. MDM2 also serves as a TP53 ubiquitin ligase that targets TP53 for ubiquitin-dependent degradation in the proteasome [2].
Treatment of cancer cells expressing wild-type TP53 with TP53-MDM2 interaction antagonists should result in the concurrent transcriptional activation of TP53 downstream genes, cell cycle arrest, and apoptosis.

Cell viability
Nutlin effects on cell viability were investigated using the MTT assay [3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyl tetrazolium bromide] (Sigma-Aldrich Chimie, Saint-Quentin-Fallavier, France) as an indicator of metabolically-active cells [16]. Known number (2000 or 3000) of STS cells was transferred into 96-well plates incubated for 24 h before addition of test compound. Cells were then exposed for 72 h at 37°C to an increasing concentration range of RG7388, for 48 h to doxorubicine or gemcitabine, MTT at the final concentration of 0.5mg/ml was added and, following incubation for 3h, Formazan crystals were dissolved in DMSO. Absorbance of the colored solution was measured on a microplate-photometre (Bio-Tek Instruments, Colmar, France) using a test wavelength of 570 nm and a reference wavelength of 630 nm. The concentration of substance required for 50% growth inhibition (IC50) was estimated with the Graphpad Prism software (GraphPad software INC, San Diego, USA).

Cell cycle analysis
Cell cycle distribution of sensitive and resistant cell lines was studied by DNA content by fluorescence-activated cells sorting and analyzed using the Cell Quest Pro software (BD Biosciences, San Jose, USA). After 48h of treatment with two different concentrations of RG7388, cells were centrifuged at 1500g for 5 min and washed twice with PBS. Cells were then fixed with 70% ethanol at 4°C overnight. Ethanol was removed and cells washed twice with PBS. 300μl of a propidium iodide, ribonuclease-containing solution were added to cells and then analyzed by FACS. Data were analyzed with FlowJo v.7.6.3. software and results expressed in terms of percentage of cells in a given phase of cycle.

Apoptosis
The mitochondrial membrane depolarization (ΔCm) was measured by tetramethyl-rhodamine ethyl staining (TMRME). 5000 cells were seeded in 96-well plates, incubated 24h before adding an increasing concentration of RG7388 over 48h. Cells were loaded with TMRM (200nM) for 30 min at 37°C and then trypsinized and harvested in saline extracellular solution. Cells were then analyzed on a FAC scan flow cytometer (BD biosciences, San Jose, USA).
For apoptosis assessment, 1.105 cells were seeded in 6-wells plates and after 24h cells were treated by Nutlin for 72hr and exposed to FITC-Annexin and propidium iodide (PI) according to the manufacturer's protocol (BD Biosciences, Erembodegem, Belgium). It allows distinguishing annexin V positive cells in early apoptosis, versus annexin V and PI positives cells in late apoptosis. Cells were analyzed in cytomerty using FL1 Annexin-V, whereas FL2 was used for Propidium iodide.

Confocal microscopy
Cells were seeded on coverslips and treated with RG-7388 for 72 hours. Slides were then washed twice with PBS, fixed in formaldehyde 4% and incubated with anti-LC3IIB monoclonal antibody (Sigma Aldrich, Saint Louis, USA) overnight, and then with a goat anti-rabbit Alexa fluor 488 antibody (Invitrogen, Paisley, United Kingdom). Slides were then counterstained by 4,6-diamidino-2-phenylindole (Hoechst).

RnaSeq protocols for Illumina HiSeq2000 NGS
RnaSeq was sequenced by the Centro nacional de análisis genómico (Barcelona, Spain) using standard Illumina (Illumina Inc., 9885 Towne Centre Drive, San Diego, CA, USA) protocols. Briefly, the library from total RNA was prepared using the TruSeq™ Stranded mRNA Sample Preparation kit (Illumina Inc.) according to manufacturer's protocol. 1μg of total RNA was used for poly-A based mRNA enrichment selection using oligo-dT magnetic beads followed by fragmentation by divalent cations at an elevated temperature resulting into fragments of 80-450nt, with the major peak at 160nt. First strand cDNA synthesis by random hexamers and reverse transcriptase was followed by the second strand cDNA synthesis, performed in the presence of dUTP instead of dTTP. Blunt-ended double stranded cDNA was 3´adenylated and the 3´-"T" nucleotide at the Illumina indexed adapters was used for the adapters' ligation. The ligation product was amplified with 15 cycles of PCR.
Each library was sequenced using TruSeq SBS Kit v3-HS, in paired end mode with the read length 2x101bp. We generated 135 million minimally-paired end reads for each sample run in one sequencing lane on HiSeq2000 (Illumina, Inc) following the manufacturer's protocol. Images analysis, base calling and quality scoring of the run were processed using the manufacturer's software Real Time Analysis (RTA 1.13.48) and followed by generation of FASTQ sequence files by CASAVA.

NGS RnaSeq sequence alignment and quality control pipeline
Raw RnaSeq sequences where quality controlled using a set of published tools in order to produce curated reads. Briefly, the package sickle [18] was used to trim raw reads of the 5' and 3' low quality bases (phred cut off 20, max trim size 30nc, package SeqPrep (St John GitHub, https://github.com/jstjohn/SeqPrep) was used to remove sequencing adaptors from the raw reads. The application of these standard packages revealed a proportion of RNA fragments of smaller size, whose R1 and R2 paired end reads were overlapping. To keep exploiting these fragments without biasing the count of reads in reads overlapping regions, we developed a home-made script (awk language) that merged the overlapping R1 and R2 and split it into new non-overlapping read sequences.
Curated read sequences were then aligned using TopHat2 [19] on both genome and transcriptome with fragment size 175 +/-70. Following, we applied strict quality control of aligned reads, by filtering the reads whose alignment score was lower than 20 (samtools) and removing the reads with identical align starting positions that were possibly caused by PCR amplification during the creation of the RnaSeq libraries (picard MarkDuplicates at http://broadinstitute. github.io/picard/).
The statistics of the number of raw, curated and aligned sequences and curated aligned sequences was that created with a home-made tool (bash, awk)

NGS allelic variant detection and annotation pipeline
Positions in the transcriptome showing at least one alternative allele were detected via samtools mpileup [20] taking care to filter out bases with less than 20 as phred quality score. Variant Calling format information was created via bcftools view (v 0.1.19) of the samtools suite.
Positions detected as variant in several samples (sensitive and resistance cell lines) were merged using a homemade script (bash awk).
Variant positions were annotated using the Annovar tool( [21]. Briefly, Annovar determines the type of genomic location of the variant for example exonic, intergenic, and its class (e.g. missense, nonsense), its presence in population genetics databases dbSNP [22] and 1000g database [23] and somatic variants (COSMIC, http://cancer.sanger.ac.uk/cancergenome/projects/ cosmic). For variants in coding positions the impact of the aminoacid change in the protein domains is reported via the use of Sift [24] and PolyPhen2 [25] databases.
We obtained on average 135M and 120M Paired End reads respectively after sequencing and post-alignment quality controls. As a starting point of our analysis we produced a table, based on post-alignment controlled reads, collecting all genetic variants (either polymorphic or not, single nucleotides or Indels) found in sequences in at least one of the sensitive or resistant samples, of any of the three cell lines. The table resumes the variants that present at least 1 alternative allele counts in any sample, each base having sequencing quality (phred score) of at least 20, without additional filters. The table, containing 285,999 variants (herein called Merged Variant Table) allows checking the alternative and reference allele counts at a position, in any of the sequenced samples.
We extracted from the Merged Variant table the variants detected in sensitive and resistant cells. We then included in the analysis only variants whose genotype was estimated as a nonhomozygous reference in both samples, with more than 10 reads at the variant position, localized in coding regions and of non-synonymous type. We included all variants whose positions were reported in the COSMIC database (http://cancer.sanger.ac.uk/cancergenome/projects/cosmic) and whose alternative allele frequency (AF) in the Caucasian population (CEU), as reported in 1000g, was lower than 5%. If a variant was neither reported in 1000g CEU nor dbSNP, it was included only if either SIFT or Polyphen2 HDIV functional interpretations were reported as "deleterious" for the protein. If a variant was not reported in 1000g CEU but reported in dbSNP, it was included only if both SIFT or Polyphen2 HDIV functional interpretations were reported as "deleterious" for the protein. Also, if a variation was located in a region of the genome that is annotated as Segmental Duplication [26], then it was not reported. In order to distinguish variants that were shared between sensitive and resistant cells, we tested whether the Alternative Allele Ratios (AAR) were statistically different between the two conditions (Fisher test P< = 10-3, and difference between AAR < = 0.1). Variants whose differences in AAR were not significant and superior to 0.30 in both conditions were associated to the ancestral primary tumor clone, while the other variants identified the differences between the clones sensitive and secondary resistant to RG7388.

Detection of genomic ploidy and Differential Gene Expression/IPA
RnaSeq was used to estimate the ploidy of the tumor samples, where the term ploidy defines the total number of parental alleles in large genomic regions (see S1 Methods). The RnaSeq Differential Gene Expression between two NGS samples was calculated using the methodology described in S1 Methods and implemented in R language. The selection criteria for detection of Differentially Expressed Genes were based on a Fold Change greater or equal to 2.5 and an adjusted p-value of less than or equal to 1% (Benjamini Hochberg adjustment). Data were analyzed through the use of QIAGEN's Ingenuity 1 Pathway Analysis (IPA 1 , QIAGEN Redwood City,www.qiagen.com/ingenuity) in order to predict the interactions between the genes and the potential effects on cellular pathways and networks.

Results
RG7388 activates the TP53 pathway, induces significant proliferation inhibition, cell-cycle arrest and apoptosis in MDM2-amplified and nonamplified STS As predicted by the mechanistic model of TP53 regulation, the nutlin compound RG7388 inhibited significantly the proliferation of 5 out of 11 STS cell lines with no TP53 mutations as assessed by Sanger sequencing (IC50: 2-50 nM), but not of the 6 out of 11 cell lines with TP53 mutations (Fig 1A). The most sensitive cell lines were the DDLPS cell lines IB111 and IB115 characterized by an amplification of the MDM2 gene and the extraskeletal osteosarcoma cell line IB128 characterized with no alteration of the MDM2 gene copy numbers. In agreement with the mechanism of action of nutlins, treatment of wild-type TP53 (Sanger sequencing) STS cell lines with RG7388 showed an accumulation of the TP53 protein and its targets, P21 and MDM2, as revealed by Western blotting (Fig 1B and 1C).
One of the main cellular functions of activated TP53 is blocking cell cycle progression in the G1 and G2 phase. Treatment of exponentially proliferating STS cell lines with RG7388 for 48 hours led to a dose-dependent cell cycle block in G1 and G2/M phase and a depletion of the S phase compartment (Fig 2A and 2B). One of the other main functions of activated TP53 is induction of apoptosis. Exposure of exponentially proliferating STS cell lines to RG7388 RO5503781 for 72 hours led to the induction of apoptosis in a dose-dependent manner as revealed by an increase in the percentage of TMRM-staining cells ( Fig 2C). The 3 STS cell lines that elicited the most significant apoptotic responses were IB111, IB115 and IB128.
Repeated exposures of STS cell lines sensitive to RG7388 lead to the emergence of sub-cell lines strongly resistant to the RG7388 anti-tumor activity Among the entire panel of STS cell lines tested, RG7388 displayed the highest cell growth inhibition effects in IB111, IB115 and IB128. In order to investigate mechanisms of secondary resistance to RG7388, IB111, IB115 and IB128 cells were treated with RG7388 for 3 days. The cells were then rinsed to remove RG7388, and the remaining cells were expanded in normal medium (minus RG7388). We repeated the process four times, and we obtained populations that survived 1-4 rounds of RG7388 treatment (P1 for one round-P4 for four rounds). We compared the extent to which IB111, IB115 and IB128 underwent apoptosis and cell cycle inhibition when treated for 3 days with RG7388. Results indicated that the selected populations became progressively more resistant to cell growth inhibition ( Fig 3A) to apoptosis ( Fig 3B) and cell cycle inhibition (data not shown). For instance, whereas 1 μM of RG7388 induced significant apoptosis in parental IB115 cells (80.7% apoptosis) after 3 days of treatment with RG7388, this effect was significantly reduced in IB115P2 (36.7% apoptosis) and IB111P4 (8.8% apoptosis) cells. Beside, we confirmed resistance to apoptosis at 1μM of RG7388 in resistant cell lines IB115P4 (69.4% apoptosis) and IB111P4 (28.4% apoptosis) versus parental cell lines IB115 (98.7% apoptosis) and IB111 (76.7% apoptosis) by FITC annexin-V and propidium iodide assay ( Fig 3C).  STS sub-lines with secondary resistance to RG7388 display sensitivity to cytotoxic drugs similar to the parental cell lines In order to assess specificity of secondary resistance to RG7388, we compare the sensitivity profiles of IB111, IB115, IB128 and their RG7388-resistant counterparts IB111P4, IB115P4 and IB128P4 to doxorubicin and gemcitabine, two cytotoxic drugs commonly used to treat STS patients with advanced disease. Doxorubicin and gemcitabine significantly inhibited the proliferation with the same extent in the parental cell lines and the RG7388-resistant sublines (S1 Fig). This indicated that secondary resistance to RG7388 in IB111P4, IB115P4, IB128P4 did not result from a multidrug resistance phenotype and may involve specific mechanisms related to the mechanisms of action of RG7388.

Deep sequencing identifies ploidy variations involved in secondary resistance of STS cells to RG7388
Several studies have shown correlations between gains and losses of chromosomes and cancerspecific drug resistance [27][28][29]. We wanted to investigate whether regional ploidy variations were associated with secondary resistance to RG7388 in STS cells. Results are summarized in Table 1 and S2 Fig. IB111 had the greater number of allelic imbalances, followed by IB115 and IB128. Some of them were shared by both the parental and secondary resistant cell lines. However, all the secondary resistant cell lines were characterized by at least one specific ploidy variation not found in the parental ones. For instance, the secondary resistant IB111 cell line was also characterized by an unbalanced tetraploidy of the 12q13.2-qter region in which MDM2 is located. The secondary resistant IB128 cell line was also characterized by an unbalanced pentaploidy of the short arm of chromosome 17 including the TP53 gene.

Secondary resistance of STS cells to RG7388 is characterized by clonal selection of cancer gene mutations
To identify changes in the mutation profiles of secondary resistant cells, we compared the abundance of somatic mutations found in the parental and secondary resistant cell lines. For each cell line, we examined a conservative list of mutations and used Fisher's exact test to identify those specifically associated with secondary resistance (see material and methods). In all Table 1. Ploidy variations between secondary resistant and parental soft-tissue sarcoma cell lines.
only sensitive clone only resistant clone -: absence of allelic imbalance cell lines the majority of variants were common to the parental and resistant cells (Fig 4). The number of non-synonymous mutations with significant changes in normalized abundance between parental and resistant cell lines ranged from 10 to 23 for each case. These include mutations in well-known cancer genes, genes linked to drug resistance and drug metabolism, and genes not previously associated with carcinogenesis or therapy resistance ( Gene expression analysis identifies defects in apoptosome activity and autophagy induction as a mechanism of secondary resistance to RG7388 We also analyzed the sequencing data to identify differences in gene expression between secondary resistant and parental cell lines. We found 196, 105 and 370 genes with 2.5X or more fold changes (adjusted p-value < 0.01) between the parental and secondary resistant cell lines in IB111, IB115 and IB128 respectively (S2 Table). We then applied these genes to Ingenuity Pathway Analysis software (IPA 1 ,QIAGEN Redwood City, www.qiagen.com/ingenuity). The results showed that these genes were mainly enriched in proliferative, growth and movement networks ( Table 3). We focus particularly on the genes differentially expressed between the resistant and parental IB111 cells in order to identify some genetic and transcriptomic alterations that may compensate for the lack of TP53 mutations observed in the IB111 resistant cells. Strikingly, in the resistant IB111 cells, we found a significant alteration of genes involved in the regulation of apoptosis and more particularly, a down-regulation of pro-apoptotic activators such as BMF, BIM and PUMA (S3 Fig). Among them, the most strongly down-regulated was BMF with a fold change of 17. BMF also has a role as a regulator of autophagy. Since autophagy has been associated with resistance to radiation, chemotherapy, and targeted agents [30][31][32][33][34][35][36][37][38][39], we decided to investigate whether secondary resistance to RG7388 was associated with autophagy induction in STS cell lines. When autophagy is not activated, LC3 is localized homogeneously in the cytoplasm, while upon initiation of autophagy, it associates with the membrane of autophagosomes. Since an increase in LC3-II levels or GFP-LC3 vesicles can occur not only due to increased autophagosome synthesis but also due to impaired autophagosome-lysosome fusion, we assessed LC3-II levels also in the presence of chloroquine, a blocker of LC3-II degradation. Analysis of LC3-II levels by western blotting and by fluorescence microscopy allowed us to detect autophagy induction in the resistant IB111 cells but not in the IB115 and the IB128 resistant cells (Fig 5A and 5B). We therefore wondered whether blocking the induction of autophagy with chloroquine could restore sensitivity to RG7388 in IB111 resistant cells. The RG7388 IC50's were significantly reduced (divided by 6) in the secondary resistant IB111 cells when they were treated with chloroquine whereas RG7388 IC50 was almost unchanged in parental IB111 cells treated with chloroquine in the same conditions. (Fig 5C).

Discussion
The new paradigm of targeted therapies has dramatically impacted oncology practice with the discovery and development of 'personalized' anti-cancer drug medicines that produce remarkable clinical responses in a subset of patients with advanced disease. While there is rapidly growing enthusiasm for this new paradigm, there is also increasing realization that such targeted therapies suffer from the same major limitation as traditional chemotherapy-clinical benefits are of limited duration due to secondary drug resistance.
Establishing specific molecular mechanisms of resistance to classical cytotoxic drugs has been challenging due to the nonspecific nature of their antitumor mechanisms of action. On the other hand, identification of mechanisms of acquired resistance to targeted therapies is crucial since the discovery of such mechanisms can prompt the development of strategies specifically designed to overcome them. The design of new drugs to treat lung cancer patients who become resistant to EGFR-directed therapy due to the emergence of the T790M secondary mutations represents a perfect example of the importance of research efforts in the field of secondary resistance [40].
In Phase I studies, RG7388 has been associated with disease control in a significant subset of sarcoma patients. However all of them had disease progression within a median of 6 months    after treatment onset suggesting the occurrence of mechanisms of secondary resistance involved in tumor progression [14]. We report here the first comprehensive analysis of genetic mechanisms involved in secondary resistance to nutlin using next-generation sequencing of isogenic pairs of nutlin-sensitive and nutlin-resistant STS cell lines. The aberrations involved in secondary resistance we described here can be considered specific to RG7388 given that no cross-resistance was noted with other drugs commonly used for the management of STS patients such as doxorubicin or gemcitabine.
By performing deep sequencing of secondary resistant and parental sensitive cell lines, we showed here that secondary resistance to nutlin was associated with several mutations whose Among them, we identified TP53 mutations in the secondary resistant cells that were not present in the parental counterpart. Previous reports investigating mechanisms of secondary resistance to nutlin in neuroblastoma and osteosarcoma cell lines suggested that exposure to nutlin induces the emergence of TP53-mutated clones [41,42]. Whether these mutations appeared de novo in the secondary resistant cells or were already present in a minority of clones was not known. Using deep sequencing, our results suggest that nutlin favors the occurrence of TP53 mutations in initially TP53 wild-type cells. Besides TP53 mutations, our data identified a subset of other gene mutations that were positively selected in secondary resistant cells, many of which have been previously involved in the control of cell growth or apoptosis. Some of them may represent 'passenger' mutations or false-positives, but some are likely to contribute to resistance to RG7388. We highlight two examples here. In IB128 resistant cells, we observed a significant increase in the abundance of an activating mutation in N-RAS in the IB128 resistant cells. Activation of the MAPK/ERK kinase (MEK)1/2 kinase pathway has been shown to impair TP53-dependent apoptosis in U20S cells treated with nutlin [43]. Recent data showed that combining the MDM2-p53 protein-protein interaction inhibitor SAR405838 and the MEK inhibitor pimasertib resulted in synergic anti-tumor activity in several K-RAS, N-RAS and B-RAF mutant tumor models [44]. The possibility that combining a MEK inhibitor with nutlin could overcome secondary resistance deserves further investigation. In IB111 resistant cells, we observed an increase in abundance of a missense mutation of the ECT2 gene, the expression of which has been shown to be normally down-regulated by Nutlin in a TP53-dependent manner [45]. Inactivation of ECT2 was shown to be sufficient to prevent cell death induced by ionizing radiation underlying its potential important role in resistance to cancer therapy resistance [46].
We have also found that secondary resistance to RG7388 was associated with few acquired ploidy variations. This observation fits with the results of a previous study showing that instable aneuploidy is associated with acquired drug resistance and that nutlin treatment of U2OS osteosarcoma cells and HCT116 colon cancer cells can result in the emergence of tetraploid cells [47]. Here we have extended these previous observations, using deep sequencing, demonstrating that the majority of allelic imbalances were shared by the parental and resistant cells. Altered intrinsic tumor sensitivity to RG7388 was associated with only a few more aberrations. Some of them targeted genomic regions including genes that are directly linked to nutlin activity such as MDM2 (12q13-15) for IB111 resistant cells or TP53 (17p) for the IB128 resistant cells. This implies that some ancestral cells that carried this complement of ploidy aberrations emerged under RG7388 selection pressure. These additional ploidy aberrations demarcate the split between fully clonal versus subclonal ploidy variations.
It was recently observed that tumor cells vary in terms of their apoptotic activity potential, and that alteration of the apoptotic machinery is a cause of chemoresistance and oncogenic transformation [48][49] By comparing gene expression of parental and secondary resistant cells, we observed in IB111 resistant cells a strong down-regulation of genes encoding for BH3-only factors such as BIM, BMF and PUMA. Interestingly, down-regulation of these genes has been associated with resistance to targeted therapies such as imatinib in bcr-abl+ leukemia cells [50] or vemurafenib in human melanoma cells [51]. Besides being a pro-apoptotic molecule, BMF is also an important regulator of autophagy [52]. Since autophagy has been associated with resistance to radiation, chemotherapy, and targeted agents [30][31][32][33][34][35][36][37][38][39], we decided to investigate whether secondary resistance to RG7388 was associated with autophagy induction in STS cell lines. Interestingly, we found that IB111 resistant cells showed pronounced autophagy whereas no autophagy was observed in IB115 and IB128 resistant cells. Moreover, we have demonstrated that autophagy inhibition significantly enhanced RG7388-induced cell death in IB111 resistant cells. We report here that down-regulation of genes encoding pro-apoptotic molecules and autophagy induction are important mechanisms involved in secondary resistance to nutlin. These observations indicate that in the presence of nutlin autophagy can promote cell adaptation and survival.
In conclusion, our results demonstrate that secondary resistance of STS cells to nutlin involves multiple and complex mechanisms and that the interplay between these mechanisms warrants extensive investigation. As such events do not appear to impact sensitivity to cytotoxic agents commonly given to STS patients, therapeutic strategies preventing their emergence should be designed. One of them could be administering nutlin for a restricted period, alternating with other drugs and then restarting, irrespective of disease progression status in the period between treatments.  Table. (DOCX) S2 Table. (DOCX)