Differential Gene and MicroRNA Expression between Etoposide Resistant and Etoposide Sensitive MCF7 Breast Cancer Cell Lines

In order to develop targeted strategies for combating drug resistance it is essential to understand it’s basic molecular mechanisms. In an exploratory study we have found several possible indicators of etoposide resistance operating in MCF7VP cells, including up-regulation of ABC transporter genes, modulation of miRNA, and alteration in copy numbers of genes.


Introduction
Major advances have been made in the treatment of breast cancer; however, it still remains the leading cause of cancer death among women worldwide [1]. In the United States breast cancer is the second leading cause of cancer death among women after lung cancer [2]. The chemotherapeutic drug etoposide is used as a salvage treatment for advanced stage breast cancer and causes DNA damage by stabilization of DNA topoisomerase II. Etoposide stabilizes the topoisomerase II -DNA covalent complex impairing the strand-rejoining activity of the enzyme causing double stranded DNA breaks to persist instead of being repaired [3]. Etoposide can delay progression of the cell cycle through the late S or early G2 phase but has no effect on tubulin assembly [3]. As a single agent, oral etoposide response rate for breast cancer was found to be around 35% [4] and the mean bioavailability of orally administered etoposide is approximately 50% [3]. In a more recent study, it was found that oral etoposide in combination with cisplatin was much more effective in patients with advanced breast cancer (pre-treated with anthracyclines) than paclitaxel [5]. Etoposide (in combination with other drugs) is also widely used for the treatment of small cell lung cancer, lymphoma, leukemia and testicular cancer among others [6]. Though etoposide is widely used as therapy for cancer patients, the fact remains that tumors often acquire resistance to the drug. The nature of drug resistance is inherently multi-factorial involving mechanisms which include alteration in drug targets, inactivation/detoxification of the drug, decreased drug up take, increased drug efflux and the disregulation of the apoptotic pathway [7]. Apart from these mechanisms, in recent years, it has been found that small, non-coding RNAs called microRNAs (miRNAs) may also have a role in regulating drug resistance. miRNAs are small (18-22 nucleotide) non-coding RNAs that are capable of silencing the expression of certain genes by binding to complementary sites in the 39UTR of these genes and causing either mRNA cleavage or translational repression [8]. miRNA abnormalities have become an emergent theme in cancer research. It was found that global miRNA expression patterns could classify human cancers according to developmental lineage and differentiation status much more accurately than mRNA expression profiling [9]. Gerard Wright first proposed the term 'resistome' to denote the genes responsible for antibiotic resistance and their precursors in bacteria [10,11]. Similar to the 'integrated network of antibiotic resistance elements' in bacteria that provide protection against chemical threats [11], some cancer cells may also possess an integrated network that protects them from chemotherapeutic agents. In our study we identified several indicators that may represent the resistance network specific to an etoposide resistant MCF7VP breast cancer cell line. This network may contain genes and miRNA related to multiple mechanistic pathways of drug resistance, however, the applicability of the study to a clinical setting remains to be determined.

Results
The MCF7VP Cell Line Shows a Higher Level of Resistance to Etoposide Compared to the Drug Sensitive Parental MCF7 Cell Line To determine the extent of etoposide resistance in both the drug sensitive and drug resistant cell lines we measured sensitivity using a cellular cytotoxicity assay. Briefly, 10,000 cells per well of a 96 well plate were plated and pre-incubated at 37 degrees C for 24 hrs after which 100 ml of various concentrations (0-100 mM) of etoposide was added. After 72 hr incubation 10 ml of CCK8 reagent was added. Two hours later the absorbance was measured at 450 nm. This data was fitted to a dose-response curve using GraphPad Prism and EC50's were calculated revealing that the MCF7VP cells were 12.5 fold more resistant to etoposide than the drug sensitive MCF7 cells (Fig. 1).
ABCC1 and ABCC6 Genes are Overexpressed in the MCF7VP Cell Line The ABC (ATP binding cassette) transporters are a diverse family that contain a number of energy dependent efflux pumps. Using TaqMan low density microfluidic arrays which contain the probes and primers (in triplicate) to detect the expression of all 48 human ABC transporters (and also endogenous controls such as beta actin) we found that ABCC1 and ABCC6 were highly expressed in the etoposide-resistant cell line compared to the sensitive cell line (Fig. 2). We validated the data using quantitative real time PCR (qRT-PCR) and found that ABCC1 was 20.5 fold over-expressed (Fig. 3A) and ABCC6 was 72.5 fold over-expressed in the MCF7VP cell line compared to the parental MCF7 cell line (Fig. 3B). We found that expression of other transporters was not appreciably increased (Fig. 2). We also validated this data using western blotting to show that both ABCC1 and ABCC6 proteins are expressed in MCF7VP cells (Fig. S1).

MCF7VP Cells have a Distinct miRNA Profile when Compared to MCF7 Drug Sensitive Cells
Using miRNA microfluidic arrays from ABI, we determined the microRNA profile for MCF7VP cells compared to drug sensitive cells (Table 1, Data File S1). Multiple miRNAs showed differential expression among the cell lines including hsa-miR-382, hsa-miR-23b and hsa-miR-885-5p, which were up-regulated (.2-fold increase) in MCF7VP cells and hsa-mir-218, hsa-miR-758 and hsa-miR-548d-5p, which were down-regulated (.2-fold decrease) in MCF7VP cells (Table 1 and Data File S1), suggesting that etoposide resistant MCF7 cells have a miRNA profile that is distinct from the MCF7 drug sensitive cell line. Of particular interest was the 7.79-fold decreased expression of Hsa-miR-218 in MCF7VP cells.

Gene Expression Profiling and Enrichment Analysis Reveals Differences in Gene Expression between the Drug Resistant and Drug Sensitive MCF7 Cells
In order to obtain a complete picture of the differences in expression levels between the drug resistant and drug sensitive cell lines, gene expression profiling was carried out using Affymetrix HG-U133 Plus 2.0 arrays and the arrays were run in triplicate which is the industry standard for microarray data. The data was RMA normalized and analyzed with Partek Genomics Suite. The top 5000 differentially expressed genes ($2 fold change, P,0.05) are displayed as a heat map (Fig. 4A). Twenty-nine genes have greater than 10-fold higher expression in the MCF7VP cell line ( Table 2) including ABCC1 and ABCC6, confirming our earlier TaqMan assay results. In addition, the data revealed that topoisomerase 2 (TOP2A), which is a target of etoposide, was down-regulated in the MCF7VP cells (The entire data set is deposited in GEO: GEO accession no. GSE28415).

Functional Analysis of the Microarray Data Predicted Several Possible Pathways of Etoposide Resistance
Gene Ontology (GO) analysis of the microarray data revealed a number of gene clusters that could be associated with drug resistance. DAVID, an annotation and integration tool, was used to explore the differentially expressed genes (top 219 genes, p,0.05) in etoposide-sensitive and resistant cell lines. We identified possible indicators in pathways that might be associated with drug resistance including extra-cellular matrix (ECM) organization, JAK-STAT signaling pathway, MAP kinase signaling pathway, and ATP binding ( We also found that extra-cellular matrix (ECM) genes are upregulated in etoposide resistant cells. The extra-cellular matrix primarily provides support to animal cells and tissues. Analysis showed that over 66% of extra-cellular region genes are upregulated in the MCF7VP cells (Fig. S2). These genes include: collagens (particularly collagen type 3A1 and collagen type 12A1), fibronectin and matrix-gla protein, which are involved in remodeling of the ECM. Some studies have shown that the ECM may mediate drug resistance by either disrupting integrin signaling or by the formation of an actual physical barrier by remodeling of the ECM [12] however, very little is known about this mode of drug resistance. IPA analysis validated that integrin signaling was one of the pathways upregulated in the microarray dataset (Fig. S3). In addition, we validated the fact that TOP2A gene was down-regulated using qRT-PCR, which showed that the gene expression of TOP2A was decreased by 2.1 fold in the MCF7VP cell line (Fig. 5D). This may be indicative of a possible modulation of drug resistance through the downregulation of TOP2A where the MCF7VPcells may circumvent drug sensitivity by down-regulating the drug target -in this case Topoisomerase 2A. Copy Number Analysis Identified Differences in Copy No. of Genes between the MCF7 and MCF7VP Cell Lines As expected, the number of copy number changes detected with the Affymetrix 6.0 chip varied greatly between the 2 cell lines. Numerous amplifications and deletions were detected on chromosomes in the MCF7VP cells when compared to the baseline of the parental MCF7 cells ( Table 2). The copy number analysis was further verified by analyzing loss of heterozygosity (LOH) of the SNPs on different chromosomes. LOH is a form of allelic imbalance that can result from the complete loss of an allele or from an increase in copy number of one allele relative to the other. We compared the copy number analysis with the gene expression data and found the both ABCC6 and ABCC1 were amplified in the chromosomal DNA and also at the RNA level. The GHR gene or growth hormone receptor involved in various signal transduction pathways (eg;cellular growth) was also up-regulated at both DNA and RNA levels ( Table 2). The MAFF oncogene which is thought to be associated with cellular stress was upregulated and showed an increase in copy number. A Circos plot was generated incorporating differential gene expression and loss of heterozygosity (Fig. 6).

Discussion
In this study we have described gene and miRNA indicative of resistance in the MCF7VP etoposide-resistant breast cancer cells. Gene expression data and copy number analysis (Table 2) revealed that the transporter genes ABCC1 and ABCC6 are significantly upregulated in MCF7VP cells. It has been previously found that the transfection of MRP1/ABCC1 expression vectors into HeLa cells resulted in these cells becoming moderately cross-resistant to several drugs including etoposide [13] and that ABCC6 transfected Chinese hamster ovary (CHO) cells show a low level of resistance to etoposide [14]. Our findings are consistent with these studies. However, we did not observe significant up-regulation of other ABC transporters including ABCB1 or ABCC3 [15] that are thought to be associated with etoposide resistance. It was recently shown that ABCC6 was five-fold up-regulated in gemcitabine resistant A549 non-small cell lung cancer cells [16] suggesting that it also may be a factor in influencing gemcitabine resistance. Heat map depicting the upregulation of ABCC6 and ABCC1 genes in MCF7VP cells. It can be seen that both ABCC6 (red) and ABCC1 (red) genes are upregulated in the MCF7VP cell lines. cDNA from the MCF7 and MCF7VP cells was used to run the Taqman Low Density ABC transporter arrays. The data was analyzed using Microsoft excel and R-Bioconductor to generate heat maps of the ABC array data. doi:10.1371/journal.pone.0045268.g002 It has been previously reported that Hsa-miR-326 can target ABCC1 39UTR [17] in MCF7VP cells and accordingly we found that this miRNA was down-regulated in the etoposide resistant cells when compared to the sensitive cells (Data File S1). Other transporters such as ABCB1 (P-gp) are known to be regulated by microRNA. It was previously found that antagomirs of hsa-miR-451 and 27a can target MDR1 in A2780DX5 cells [18]. Taken together, this suggests that microRNA's and their antagomirs could be used to modulate the expression of ABC transporters.
We also found a set of extra-cellular matrix genes were upregulated in the etoposide resistant cells. It has been proposed that the extra-cellular matrix may be an important factor in contributing to drug resistance either by the disruption of integrin signaling or by the formation of an actual physical barrier achieved by remodeling of the ECM, [12] however, very little is known about this particular mode of drug resistance. Currently there is scant evidence to support the role of ECM in drug resistance, but our findings suggest that this may be a valid factor in etoposide resistance.
Finally, we found that the copy numbers of both ABCC1 and ABCC6 were increased significantly in the MCF7VP cell line and this data correlated well with our gene expression data ( Table 2).
As expected we found differences in copy number in a several genes between the drug resistant and drug sensitive cell lines which correlated for the most part with the gene expression data (Fig. 6). This would suggest that, along with upregulation of gene expression, alteration in gene copy number may be a mechanism of activating or down-regulating genes in chemoresistance. This was also suggested in a previous study by Yasui et al 2004 [19].
In the current study we identify a number of molecular indictors of etoposide resistance in the MC7VP etoposide resistant cell line. Our work may serve as a model on which to base future investigations of etoposide resistance in diverse cancer cell lines and also in clinical patient samples.

Cell Lines and Culture Conditions
The MCF7VP (etoposide resistant) cell line was obtained from Rob Robey -originally sourced from Dr. Kenneth Cowan, University of Nebraska Medical Center, NE, USA [20] and was cultured for at least 8 passages in 4 mM etoposide in IMEM with  10% FBS and 1% Penicillin-Streptomycin. The paired MCF7 parental cell line was also cultured in IMEM with with 10% FBS and 1% Pen-Strep. MCF7 was obtained from the DTP (Developmental Therapeutics Program) NCI-Frederick that was also originally sourced from Dr. Cowan [20].

Cellular Cytotoxicity Assay to Determine Sensitivity to Etoposide
The cellular cytotoxicity kit from Dojindo (Rockville, MD) was used for this assay as per manufacturers instructions.
The cellular cytotoxicity assay allows sensitive colorimetric assays for the determination of cell viability. A highly water-soluble tetrazolium salt, WST-8*, is reduced by dehydrogenase activities in cells to give a yellow-color formazan dye, that is soluble in the tissue culture media. The amount of the formazan dye, generated by the activities of dehydrogenases in cells, is directly proportional to the number of living cells hence the toxicity of etoposide can be determined by measuring absorbance of the formazan dye. Briefly 10,000 cells per well of a 96 well plate were plated and preincubated at 37 degrees C for 24 hrs. after which 100 ml of various concentrations (0-100 mM) of etoposide was added. After 72 hrs incubation 10 ml of CCK8 reagent was added. Two hours later the absorbance was measured at 450 nm. Each experiment was carried out twice and the data pooled. Calculations were carried out with MSExcel and GraphPad prism and R.

Isolation of Total RNA, Including the miRNA Fraction
RNA was isolated from MCF7 and MCF7VP cells using the mirVana kit from ABI as per manufacturers instructions (Applied Biosystems (ABI), Carlsbad, CA) and analyzed using an Agilent bioanalyzer to determine the quality and quantity of the RNA.

Reverse Transcription and Running ABC Transporter TLDA (Taqman Low Density Arrays)
The RNA was reverse transcribed into cDNA using the high capacity reverse transcriptase kit (ABI). The cDNA was used to run TLDA-ABC (ATP binding cassette) transporter arrays or individual reactions with the relevant primers on a 7900HT real time PCR machine. The data was analyzed using Microsoft excel and R-Bioconductor to generate heat maps of the ABC array data.

Conversion of RNA to cDNA and Running miRNA Array
The isolated RNA was converted to cDNA using human miRNA primer pools A & B respectively (ABI) and the cDNA was used to run the real-time reaction in two human miRNA TLDA (Taqman low density array ABI cards version 2.0) on a 7900HT

Reverse Transcription and qRTPCR
The RNA was reverse transcribed into cDNA using the high capacity reverse transcriptase kit from ABI. Briefly, we made a reaction mixture of the following components-106 rev transcription buffer (10 ml) 256 dNTP's, (4 ml), 106 random primers (10 ml), multiscribe reverse transcriptase (5 ml), Nuclease Free Water (21 ml) and 30 ng of RNA. The PCR reaction was carried out on a thermal cycler with the following parameters: 25uC for 10 min, 37uC for 120 min, 85uC for 5 min to synthesize the cDNA. The real time reaction was carried out using taqman primers and probes specific for a particular gene. We added 2X Taqman Universal PCR master mix (10 ml) RNAase free water (8 ml), 20X target primers and probes (1 ml), cDNA sample 30 ng in a microfuge tube. This was centrifuged for 1 min at 1000 rpm. The PCR cycle was carried out in a 7900 real time instrument (ABI) under the following conditions: 50uC (2 mins) hold, 95uC (10 mins) hold. Then 40 cycles, 95uC for 15 secs and 60uC for 1 minute. The qRTPCR data was analyzed using RQ manager (ABI).

Affymetrix U133 plus 2 Arrays
Affymetrix HG-U133 Plus 2.0 arrays were run in triplicates with RNA isolated from MCF7 and MCF7VP cells at the Laboratory of Molecular Technology, Frederick, MD. Briefly, for each technical replicate (3 replicates), 100 ng of total RNA were amplified and labeled using the message amp II-biotin enhanced reagents (Ambion) according to the protocol provided by the supplier. Arrays were hybridized with 11 mg of labeled cRNA, washed, stained, and scanned according to the protocol described in Affymetrix GeneChip Expression Analysis Manual (Fluidics protocol FS450_0001). Arrays were scanned using an Affymetrix GeneChip scanner 7G, the .CEL files for each array was imported into Partek Genomics Suite and normalized using the RMA (Robust Multichip Averaging) algorithm ANOVA was used to determine significantly differentiated probe sets between samples Data was deposited in GEO -accession no. GSE28415, http:// www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc = GSE28415. Gene ontology (GO) analysis was carried out using DAVID v6.7 [21,22], gene set enrichment analysis was done with GSEA [23,24] and pathway maps were constructed with IPA (Ingenuity Pathway Analysis) software. All microarray experiments were done in compliance with MIAME guidelines.

Determination of DNA Copy Numbers Using the SNP 6.0 Array
We followed the manufacturer's instructions for the Affymetrix Genome-wide Human SNP array 6.0. A SNP array can be used to generate a virtual karyotype using software to determine the copy number of each SNP on the array and then align the SNPs in chromosomal order. To detect copy number variations genomic DNA from the MCF7 and MCF7VP samples was labeled, fragmented and hybridized to Affymetrix SNP6.0 arrays according to Affymetrix protocols.

Data Analysis
Analysis of miRNA data. The miRNA data obtained was analyzed using the StatMiner software suite (from Integromics) to determine differential expression of miRNA. Briefly, the data was imported into StatMiner, visually inspected for any anomalies, filtered and unexpressed detectors were flagged and the data was normalized with the help of the appropriate endogenous control. We then performed hierarchical clustering to generate heatmaps and also bar charts of the resulting analyzed data to determine miRNA expression profiles both in StatMiner and in Rbioconductor.
Analysis of microarray data. Using the gene expression workflow in Partek (PartekH software, version 6.4 Copyright ß 2010,Partek Inc., St. Louis, MO, USA,) we converted the .CEL files generated by the Affymetrix data file into.FMT files (which can be later converted to .XLS files) and normalized the files using RMA normalization procedure. We identified differentially expressed genes using ANOVA (analysis of variance) analysis, filter (,2 fold change, p value .0.05) and displayed the data as a heat map and exported to an excel file to create a gene list of differentially expressed genes between the drug sensitive and drug resistant cells.
Quantitative Real Time PCR data analysis. Relative quantification of the expression level of each transcript in each sample was calculated using the delta delta CT method in RQ manager (ABI-Prism). Human beta actin was used as the endogenous reference gene. DAVID (v6.7 from NIAID, NIH). We imported the filtered microarray gene list using affymetrix gene identifiers into the DAVID workflow and conducted a functional annotation clustering analysis. This yielded functional annotation clustering scores for groups of genes to discover enriched functionally related gene groups. This was then compared between the drug sensitive and drug resistant cell lines. Using DAVID gene ontology analysis we also investigated the function of these genes and redirected to relevant literature.
GSEA (Broad Institute). We used GSEA analysis to also calculate gene enrichment scores and also graphically represented these scores. Briefly, we imported the filtered microarray data and ran the analysis using a fixed number of permutations.
Copy number data analysis. Paired copy number analysis was performed on DNA from both cell lines with the Partek Genomics Suite, using SNP 6.0 intensity data (Affymetrix '.cel' files) with the DNA of MCF7 used as a reference. Differences were analyzed between the paired samples based on the gene involved and the base pair start and end of the region of change. Changes were described as either amplification (if the average intensity value was over the diploid value of 2) or deletion (if under 2). The genomic segmentation algorithm was used with default parameters to identify regions of amplification or deletion and a t-test was used  Fig. 5A. Gene set enrichment analysis of microrray data depicting the enrichment of genes in the JAK-STAT signaling pathway. The GSEA software was used to calculate the enrichment levels. Fig. 5B. Gene set enrichment analysis of microrray data depicting the enrichment of genes in the MAP kinase signaling pathway. The GSEA software was used to calculate the enrichment levels. Figure 5C. Gene set enrichment analysis depicting enrichment of ECM (extra-cellular matrix) genes in MCF7VP cells. The GSEA algorithm was used to calculate the enrichment levels. Figure 5D. Down-regulation of TOPO2A (the drug target of etoposide) in MCF7VP cells. The microarray data was validated by qRTPCR as depicted in the bar chart which shows differences in fold change. doi:10.1371/journal.pone.0045268.g005 to assess significance of the regions. The copy number and gene expression data were displayed as a circos plot [25].  . Circos plot incorporating differential gene expression and LOH. Chromosome numbers and bands are identified in the outermost ring. Other tracks from outer to inner represent: amplifications (red) and deletions (blue) in MCF7VP compare to MCF7 using Affy 6.0 SNP data (Partek); differential expression (p,0.05) between MCF7 and MCF7VP (AffyU133), ln(MCF7/MCFVP) values plotted; LOH in MCF7 derived from Affy 6.0 SNP data (green); LOH in MCF7VP derived from Affy 6.0 SNP data (black).