New insights in gene expression alteration as effect of doxorubicin drug resistance in triple negative breast cancer cells

Triple negative breast cancer (TNBC) is a heterogeneous disease with aggressive behavior and an unfavorable prognosis rate. Due to the lack of surface receptors, TNBC must be intensely investigated in order to establish a suitable treatment for patients with this pathology. Chemoresistance is an important reason for therapeutic failure in TNBC. The aim of this study was to investigate the effect of doxorubicin in TNBC cell lines and to highlight cellular and molecular alterations after a long exposure to doxorubicin. The results revealed that doxorubicin significantly increased the half maximal inhibitory concentration (IC50) values at P12 and P24 compared to parenteral cells P0. Modifications in gene expression were investigated through microarray technique, and for detection of mutational pattern was used Next Generation Sequencing (NGS). 196 upregulated and 115 downregulated genes were observed as effect of multiple dose exposure, and 15 overexpressed genes were found to be involved in drug resistance. Also, the presence of some additional mutations in both cell lines was observed. The outcomes of this research may provide novel biomarkers for drug resistance in TNBC. Also, this activity can highlight the potential mechanisms associated with drug resistance, as well as the potential therapies to counteract these mechanisms.


Background
Triple-negative breast cancer (TNBC) is immunohistochemically defined by the lack of three important receptors: estrogen (ER), progesterone (PR) and human epidermal growth factor receptor 2 (HER2) [1,2]. TNBC is a less frequent phenotype, being around 15-20% of all breast cancers. Despite this, TNBC is generally diagnosed in young patients, presenting a high metastatic capacity and an unfavorable prognostic rate [3]. TNBC is a relevant clinical challenge considering that this cancer subtype does not respond to endocrine therapy or other targeted agents [1,4]. Meanwhile, conventional chemotherapy and radiotherapy remain the main important alternative for these patients [5].
Chemotherapeutics is widely used as treatment strategy against TNBC tumors. However, the effectiveness of the treatment can be affected by the activation of the resistance related mechanisms [5,6]. A typical and common treatment used for breast cancer, is represented by doxorubicin (DNA damaging agent) in combination with paclitaxel (microtubule-stabilizing drug) or/and cyclophosphamide [7][8][9].
The mechanisms of primary or acquired chemoresistance to doxorubicin still remain to be deciphered. Thus, our study is focused on the regulatory pathways responsible for resistance to therapy and possible specific targets that could help optimize patient responses to this drug [10,11]. Chemoresistance is correlated with genetic alterations, that can activate pro-survival signaling, DNA damage repair, drug efflux or epithelial-mesenchymal transition [12][13][14]. The recognition of molecular features responsible for a particular response to chemotherapy, particularly transcriptomics and genetic alterations, proved to have a significant impact on cancer research. These features could be represented by biomarkers for resistance or sensitivity to a particular drug, or specific mechanistic alteration that can be a starting point for overcoming this resistance mechanism [15].
In this study, we developed two doxorubicin-resistant TNBC cell lines (Hs578T/Dox and MDA-MB-231/Dox) by multiple dose exposure of TNBC cells to doxorubicin, followed by evaluation of the alteration at genetic and transcriptomic levels. In this sense, the evaluation of gene expression patterns as effect of multiple dose exposure to doxorubicin was explored. The response of TNBC cell lines (Hs578T and MDA-MB-231) to doxorubicin was examined after exposure to 50 nM doxorubicin for 12, respectively 24 passages, followed by the evaluation of morphological alterations, along with genetic and transcriptomic patterns.

Assessment of doxorubicin sensitivity
The assessment of the half maximal inhibitory concentration (IC 50 ) of the parental (Hs578T and MDA-MB-231) and drug-resistant cells (Hs578T/Dox and MDA-MB-231/ Dox) was performed through MTT assays in order to evaluate the inhibitory effects on cell proliferation. In brief, at a seeding density of 1.2 × 10 4 cells/well, the cells were plated in 96-well plates and treated with stepwise concentrations of doxorubicin. Cellular viability and cytotoxicity were evaluated after 48 h by adding 1 mg/ml MTT solution and withdrawn after 2 h of incubation. As a final step, 100 μl of dimethyl sulfoxide was added in each well and the absorbance was detected at 570 nm using a Bio-Tek microplate reader.

Cytoskeletal evaluation
Morphological traits were evaluated through confocal microscopy using specific dyes for actin-filaments (Phalloidin, green dye) and nucleus (DAPI, blue dye). Modifications that occur post-therapy were evaluated at passage P0, P12 and P24 in order to increase the effects of doxorubicin on both TNBC cell lines. For this evaluation, treated cells were fixed and permeabilized with 4% paraformaldehyde, respectively 0.5% Triton X. Moreover, treated cells were stained with 200 nM Phalloidin dye for 30 min followed by an additional 1 min incubation with DAPI dye. The coverslips were mounted with 90% glycerol. The fluorescence microscopically images were captured using UPLSAPO40X2 (NA:0.95, Olympus Japan).
DNA fragmentation using comet assay DNA fragmentation using Comet assay as effect of serial exposure to doxorubine (P0, P12 and P24) on Hs578T and MDA-MB-231 was evaluated using the alkaline single cell gel electrophoresis assay or Comet assay, by Tice' protocol as described previously [16,17].
Genetic alteration evaluation using next-generation sequencing panel on ion torrent DNA was extracted using the Purelink Genomic DNA minikit following the manufacturer instruction. We used two triple negative breast cancer cell lines (HS578T and MDA-MB-231) samples that were treated with doxorubicin. DNA was extracted from samples at passage 0 (P0), 12 (P12) and 24(P24) after treatment. The DNA concentration was quantified using NanoDrop and were obtained concentrations between 54.19-115.4 ng/μl. 20 ng of DNA were used for sequencing using the Ion AmpliSeq Cancer Hotspot Panel v2 (ThermoFisher Scientific) and the Ion AmpliSeq Library 2.0 kit (Ther-moFisher Scientific). The Ion AmpliSeq Cancer Hotspot Panel v2 consists of primers for hotspot evaluation in the following genes : ABL1, AKT1,ALK, APC, ATM,  BRAF, CDH1, CDK2A, CSF1R, CTNNB1, EGFR, ERBB2,  ERBB4, EZH2, FBXW7, FGFR1, FGFR2, FGFR3, FLT3,  GNA11, GNAQ, GNAS, HNF1A, IDH1, IDH2, JAK2 After library preparation, the samples were purified using the AMpure XP Beads (Bechman Coulter). The purified libraries were quantified using the fluorometer Qubit 2.0 and the Qubit HS DNA kit. For template synthesis, libraries were diluted to 100pM and multiplex for libraries on Ion 316 Chip (ThermoFisher Scientific). The sequencing process was performed on the Ion Torrent PGM Machine (ThermoFisher Scientific) using the Ion PGM HI-Q Sequencing 200 kit. The data obtained after sequencing were analyzed using the Torrent Suit 5.6 and Ion Reporter 5.6 software for data trimming, alignment and variant calling. The obtained variants were filtered using the following conditions: p value ≤0.05, coverage ≥500.

Gene expression microarray evaluation
Total RNA extraction, from TNBC treated and untreated cells, was performed using TriReagent (Invitrogen) and purified using RNeasy miniprep kit (Qiagen) according to the manufacturer's instruction. The RNA concentration and quality were evaluated using Nanodrop-1000 spectrophotometer (Thermo Scientific) and Bioanalyzer (RIN ≥ 7).
The alteration of gene expression pattern was done using Agilent microarray technology using SurePrint G3 Gene Expression Microarrays (8x60k), covering 26,083 genes and 30,606 lncRNA transcripts starting from 200 ng of total RNA following the Agilent standard protocol. After hybridization step, 17 h at 65°C at 10 rpm, the arrays were washed and scanned with the Agilent scanner. Probe features were extracted from the microarray scan data using Feature Extraction software (Agilent Technologies).

qRT-PCR data validation
Validation of the microarray results was done using RT-PCR technique on both TNBC cells. In this regard, genes involved in drug resistance mechanisms were selected (IL-6, CLU, JUNB and TNSF10). GAPDH and B2M were used as reference genes. In brief, 1000 ng of total RNA was reversed transcribed into cDNA using High Capacity cDNA Reverse Transcription Kit (Applied Biosystems) and amplified using SYBR Select Master Mix (Applied Biosystems) on ViiA™7 System (10 μl reaction volume in 384-well plate). Relative quantification was done using the 2 -ΔΔCT method.

Statistical analysis
Resulted data were expressed as mean ± SD (standard deviation). The difference between experimental conditions and controls were analyzed using t test (statistically significant was considered p < 0.05). Statistical analyses were carried out using GraphPad Prism version 6, Panther, Venny and String 8.0 free version software.

Cellular viability and doxorubicin assessment on TNBC cell lines
The effect of doxorubicin on TNBC cell lines was investigated in control cells (P0), cells with multiple dose exposure for 12 serial passages (P12), respectively 24 passages (P24), with doxorubicin between 0 and 100 μM.
The MTT values are presented as % of control in relation with the log (concentration, μM) (Fig. 1a). IC 50 concentrations were calculated at each time point using GraphPad Prism software. According to IC 50 concentrations, important modifications in antiproliferative activity of treated TNBC cells were observed. The increased values of IC 50 corresponding for passage P12 respectively P24 are related to the activation of drug resistance mechanisms.
Further, morphological traits were investigated as a result of a long-term exposure to doxorubicin therapy. In Fig. 1b the significant alterations occurred in cellular morphology are highlighted through the presence of modifications in actin-filaments and nucleus structure. In the control group, both cell lines exhibited normal features, meanwhile treated TNBC cells presented modifications in their structure. After doxorubicin exposure both cells lines started to exhibit alterations in their structure suggesting the activation of apoptotic processes, as compared to the control cells (P0). Thereby, the presence of micronucleus (magenta arrows), cytoskeleton damage (blue arrow) and polynuclear cells (yellow arrow) were showed at P12 and 24 in Hs578T and MDA-MB-231 cells, respectively after a prolonged exposure with doxorubicin. No such alterations were observed in P0 cells. The micronucleus alteration is detailed in Fig. 1c, by nuclear DAPI staining, on blue channel. Also, MDA-MB-231 cells treated at passage P12 and P24 exhibited giant cells suggesting that doxorubicin induces cell death through mitotic catastrophe where cells become multinucleated and enlarged, as well as loss of plasma membrane integrity. The alterations observed in Hs578T P24 were associated with stress fibers that are related to the EMT process and to the mesenchymal traits, proving that doxorubicin treatment activates more aggressive cells.
DNA fragmentation using comet assay DNA breaks (single or double strand breaks) leads to a relaxation of the DNA, migration under electric field creating comet tails. DNA cross-links do not permit DNA unwinding. Thus, when they are produced, the DNA does not migrate and comet tails are shorter. The treatment with doxorubicin led to DNA lesions which were observed under the microscope as comet tails of different grades, being described as lesion scores (LS) and tail factor (TF), as presented in Fig. 2. DNA breaks form comet tails with different lengths depending on the severity of the lesions. In our experiment the exposure to doxorubicin did not produce significant difference between the LS and TF, compared to control cells (P0) in both cell lines, except for the Hs578T cells at P12, where the LS and TF values were lower than the control cells, probable as effect of DNA repair mechanism activation. This phenomenon was not retrieved at P24, meaning that DNA breaks were similar as in control P0 cells.

Identification of mutation signatures in drug-resistant TNBC cells
We analyzed the mutation patterns for both TNBC cell lines, Hs578T and MDA-MB-231, at passage P0, P12 respectively P24. Figure 3a and b present the mutations identified in each gene for each passage. In the case of Hs578T cell line, we observed that all three passages present the same mutations signature, with two exceptions; mutation c.215C > G in TP53 gene presented in both passages, P12 respectively P24, as well as the presence of mutation c.4732_4734delGTG in NOTCH1 gene presented in passage P24. Also, the mutation presented in TP53 gene is associated with drug response in clinical database ClinVar (Fig. 3). Meanwhile, the mutation observed in NOTCH1 gene exhibits unknown clinical implication (based on ClinVar or FATHMM data base) but is already described in the public databases dbSNP and COSMIC, the clinical significance of this mutation remains to be demonstrated. For MDA-MB-231 cell line, the mutation signatures are similar for passage P12 and P24. For passage P0 (used as control) we found only the presence of three mutated genes, BRAF, KRAS, TP53. The mutated genes in TNBC cell lines have both intronic and exonic localization. TP53 has a very low activity in the analyzed cell lines. As can be observed in the IntoGene software, the main driver genes in breast cancer are TP53 and PIK3CA. Also, by reevaluating the CliVar database we observed that mutations of ERBB4 (c.421 + 58 A > G), PIK3CA (c.352 + 40 A > G) and KDR (c.3849-24C > A) have unknown significance (Fig. 3). Also, the TP53 c.469C > T was observed in some studies on breast cancer and classified as likely pathogenic or pathogenic [18,19]. The TP53 c.839G > A mutation was also observed in early onset familial prostate cancer and classified as likely pathogenic [20].

Identification of altered genes and lncRNAs expression profiles in TNBC cell lines as effect of doxorubicin therapy
In the present investigation of expression profile was used Agilent microarray technology (8x60k slides) to identify the altered transcriptomics pattern as response to multiple dose exposure to doxorubicin therapy. Thus, it was possible the identification of the most relevant altered genes in TNBC cell lines, as well as in passage P12, respectively P24, versus control cells (P0). A cut-off value for FC of ±2 and p ≤ 0.05 was selected to determine the modifications occurred in transcriptomic patterns.  Tables 1  and 2, where top 20 altered transcripts on each experimental setting are presented. The heatmap for this data is presented in Fig. 4a and b. We overlapped the genes and lncRNAs profiling data obtained from the datasets mentioned above and the Venn diagram represented in Fig. 4c-j. The common genes and lncRNAs were altered in both TNBC cell lines in all three passages (P0, P12 respectively P24). In Fig. 5 are presented the common genes/lncRNAs in Hs578T and MBA-MB-231 cells. In this regard, we can observe the presence of a high number of genes/lncRNAs identical in both cell lines, transcripts involved mainly in modulation of different biological processes.

Pathway analysis in TNBC cell lines
The genes with an altered expression level were input into specialized software in order to obtain the important biological processes activated after the exposure to doxorubicin. Panther software highlights the main altered biological functions based on the genes signature from the TNBC cell lines at passage P12, respectively P24. Therefore, 10 relevant biological processes were identified for P12 and P24 as can be observed in Fig. 6. Moreover, the biological processes between genes with down-or upregulated expression profile in P12 and P24 conditions are similar. The altered genes are mainly    involved in cell proliferation and reproduction, and immune system processes. Using Venn diagram, the main common up-and downregulated genes between P12 and P24 in same cell lines were highlighted (Fig. 7a). The altered up-and downregulated genes were used to generate graphical representation of molecular networks using String software, which displays the specific interactions between transcripts (Fig. 7b). The altered genes critically regulated progression and cell fate by activation of tremendous processes such as TNF signaling pathway or cytokinecytokine receptor interaction.
Altered genes can regulate essential processes responsible for cancer development. Therefore, an additional gene enrichment analysis using GOrilla software presents the involvement of this alterations in different molecular functions proving that doxorubicin has intricate features on nucleic acid binding activity and particular DNA binding. This gene enrichment analysis highlighted the fact that Hs578T cell line (Fig. 8a) exhibits more altered processes compared to MDA-MB-231 cell line (Fig. 8b).

Construction of a gene expression network involved in drug resistance
Despite important progress in cancer treatment, acquired resistance to chemotherapeutic drugs still remain a major obstacle in patient treatment and overall outcome. Anticancer drug resistance is activated via numerous mechanisms, and microarrays technique offer a new approach to analyze cellular pathways involved in drug resistance mechanisms being used to predict the unexpected side effects.
Using Venn diagram, we overlapped and highlighted the common genes between both TNBC cell lines at passage P12 and P24, and drug resistance genes list (list downloaded from NCBI). This overlapping analysis identified 15 genes with a significant involvement in drug resistance mechanisms. Between the 15 genes we identified also ABCC3 and ABCC6, members of the superfamily ABC transporters. Through String software, the interaction network between common genes involved in drug resistance highlighted the JUNB, CLU, IL-6, TNFSF10 genes as important regulators involved in processes including invasion/metastasis, apoptosis and resistance (Fig. 9).

Validation of genes involved in drug resistance by RT-PCR technique
The calculation of gene expression FC used the ΔΔC T method and B2M and GAPDH as the housekeeping gene (Fig. 10). In this section, genes including IL-6, CLU, JUNB and TNFSF10 were analysed in both TNBC cell lines, specifically for each passage. In Hs578T cell line, the expression levels for IL-6 (p = 0.0382) and TNSFS10 (p = 0.0368) are statistically significant in P12 group compared to P0 group. Regarding the CLU gene, the relative expression level is slightly overexpressed, meanwhile JUNB gene exhibits no alteration level compared to control group. In the case of P24, the relative expression levels are statistically increased compared to P0 group for IL-6 (p = 0.0453), CLU (p = 0.0181), JUNB (p = 0.0083) and TNFSF10 (p = 0.0005). In MDA-MB-231 cell line, we observed that the gene expression profile for the selected genes is statistically overexpressed in P12, respectively P24 group compared to P0 group IL-6 (p = 0.0065, p = 0.0020), CLU (p = 0.0235, p = 0.0017), JUNB (p = 0.0051) and TNFSF10 (p = 0.0012, p < 0.0001). Regarding JUNB gene, the relative gene expression level is slightly increased but not statistically significant in P12 group compared to control group (P0).
CXCL1, IL-6 and TNF-α protein quantification using ELISA As additional validation step was asses the evaluation of CXCL1, IL-6 and TNF-α released for cell culture after

Discussions
Doxorubicin remains one of the most active and used chemotherapeutic agent in the treatment of early and advanced breast cancer. Tumor resistance has limited the effectiveness of this therapeutic agent in single drug treatment regiments [15,21]. It is known that doxorubicin can induce drug resistance, but the most important aspect remains the exact mechanisms behind the resistance that is still poorly understood [22]. As an outcome, some biological processes are modified, including cellular state or adaptive response and chemotherapeutic tolerance, which are reflected through an increase of IC 50 values at passage P12 and P24 for both cell lines. Also was observed the presence of modifications in the morphological traits and the presence of multinucleated giant cancer cells. This aspect was also found in literature, the frequency of giant cells being correlated with significant anticancer treatment alteration [23]. Other important alteration observed is related to the presence of stress fibers formation (Fig. 1), a mechanism necessary for EMT [24].
NGS is a very valuable tool used in disease characterization, as well as in cancer. Alterations in key genes can affect the response to therapy. In our study, we identified types of mutations involved in several cellular pathways. According to our NGS data, mutations occurred in TP53 gene were identified. TP53 is a gene that can be used as an independent prognostic factor and associated with a worse prognosis, but further investigation might be needed in order to predict the response to specific therapeutic agents [25]. In Hs578T cell line, TP53 c.215 C > G mutation is presented only in P12 and P24, meanwhile delGTG NOTCH1 c.4732_4734, a variant of unknown significance mutation was observed in P24. MDA-MB-231 cell Further, microarray data for coding and noncoding genes show a differential transcriptomics pattern between both cell lines as well as both passages (P12 and P24). Using bioinformatics tools, we found that most of the enriched GO terms were mainly common among the P12 and P24 for both cell lines treated with doxorubicin. GO analysis revealed that the differentially expressed genes were involved in cellular response, which can be an adaptive response mechanism to doxorubicin, sustained also by the comet assay data.
According to GOrilla software, an activation of the DNA binding activity was observed. This is a mechanism that is still under debate in the case of doxorubicin, but an important aspect is related to the chemoresistance process and correlated with TP53-deficient cell context [26]. The DNA damage might induce stabilization of tumor suppressor gene, TP53, and might affect cell cycle progression [26] or apoptosis related mechanisms via TNFα signaling (Fig. 7). Moreover, TNF-α was initially found to induce apoptosis in different cancer types by tumor-promoting activities including transformation, proliferation, invasion and metastasis of cancer cells [27]. A high level of TNFα is characteristic to breast cancer and has frequently been associated with a poor prognosis and an aggressive behavior [28]. TNFα has a particular role in enhancing migration and invasion in tumor cells, but the underlying mechanisms are still elusive [29,30]. In MDA-MB-231 cells was demonstrated that TNF-α increased the expression profile and activity of MMP-9 by inducing JUNB DNA binding activity, thus strengthening the concept of a pro-tumorigenic effect of TNFα in breast cancer [31]. Another process that involves TNF-α activity is the EMT, responsible for the loss of cell adhesion, disruption of cell-cell junctions, extensive actin cytoskeleton reorganization, resistance to cell death, increasing the mobility and invasiveness [32]. Doxorubicin acts as an intercalating agent having the ability to block DNA synthesis and transcription, as well as inhibiting the activity of topoisomerase type II enzyme. This process cause breaks in the genomic DNA that ultimately can lead to apoptosis [33]. These alterations induced in genomic DNA can lead to modifications responsible for the activation of critical biological processes, mainly involved in drug resistance [7]. In Fig. 8 are highlighted the main common genes involved in drug resistance mechanisms. The presence of ATPbinding cassette (ABC) family transports suggest the implication and activation of the mechanisms involved in drug response. ABCC3 and ABCC6 members, which participate directly in the active transport of drugs into subcellular organelles or influence drug distribution indirectly, were observed in our study. Balaji et al. showed that overexpression of ABCC3 is correlated with decreased drug retention; meanwhile knockdown of ABCC3 increased drug retention and cell death [34]. Overexpression of ABCC6 is able to confer low levels of resistance to several anticancer agents including doxorubicin, etoposide, daunorubicin and actinomycin [35]. Other gene found in drug resistance network is JUNB which has been associated with invasion/metastasis in breast cancer and represents an important target in diseases, associated with EMT. Also, JUNB has been involved in the earliest events of resistance development in breast cancer [36]. In addition, CLU (Clusterin) is involved in anti-apoptotic processes, development of therapy resistance, induction of EMT, all associated with cancer metastasis. Moreover, CLU has the capability to protect metastatic cells from cell death, leading to cell survival in different environment [37]. Also, CLU is a stress-activated cytoprotective chaperone targeted and upregulated by a wide range of anticancer therapies that confer treatment resistance [38]. In breast cancer, overexpression of CLU was also associated with resistance to neoadjuvant chemotherapy [39]. The significant altered gene expression signature identified can reflect biological regulation and metabolic response in the presence of doxorubicin. We found that TNBC cells developed an adaptive response against the intracellular stress induced by doxorubicin therapy during the stimulation of the resistance process. We identified a 15 gene common signature (TNF, VEGFA, IL-6, TNFSF10, CLU, ABCC6, EGR1, SNAI1, ABCC3, EPHX1, FASN, CXCL1, IL24, JUNB, TP53I11) correlated with drug resistance.

Conclusion
By understanding the molecular bases of chemoresistance in TNBC and functional pathways that are involved in drug resistance mechanisms, we tried to reveal critical aspects to improve prognosis rate and to counteract drug resistance mechanisms. Our study provided important information related to doxorubicin response mechanism correlating data of genetic and transcriptomic alteration with gene networks associated with drug resistance and cellular response. This information could be a valuable starting point for follow-up experiments that might test novel drug targets for anti-cancer treatment in order to prevent activation of drug resistance mechanisms.