Integrated Chromatin Accessibility and Transcriptome Landscapes of Doxorubicin-Resistant Breast Cancer Cells

Introduction

Breast cancer has the highest morbidity and mortality rates of all cancers in the female population, and its incidence has gradually increased over the past few decades (Sung et al., 2021). Treatment options for breast cancer are dependent on the specific biological characteristics and genetic heterogeneity of the tumor. Despite significant progress has been made regarding treatment options for patients in recent years, drug resistance in cancer cells is one of the major obstacles for effective chemotherapy to treat cancer, which causes treatment failure in almost 90% of patients with metastatic cancer (Hammond et al., 2016). With regard to breast cancer treatment, although doxorubicin, known as an anthracycline antibiotic, is currently considered to be one of the most effective agents, the resistance to doxorubicin is still a major problem in clinical practice and leads to an unsuccessful outcome in many patients. Overcoming doxorubicin resistance would represent a major advance in effective breast cancer management.

Doxorubicin induces cell death and cell growth arrest primarily by inhibition of the topoisomerase II activity, DNA intercalation and production of free radicals, leading to apoptosis, autophagy, senescence and necrosis of the fast-growing cancer cells (Burgess et al., 2008; Meredith and Dass, 2016). In spite of the numerous genes, non-coding RNA and signaling pathways associated with doxorubicin resistance in breast cancer have been identified (Tormo et al., 2019; Ciocan-Cartita et al., 2020; Tanaka et al., 2020), the epigenetic regulatory mechanism, especially the role of chromatin-mediated processes in the regulation of gene expression, has not been investigated yet. Therefore, exploring the regulatory elements of doxorubicin-resistant genes and their corresponding transcription factors (TFs) in breast cancer is critically important for the elucidation of the disease process.

Eukaryotic genomes, as is well known, are extremely compacted in chromatin, a complex of DNA and proteins, and approximately 1% of the genome is composed of accessible elements, such as promoters, enhancers, and other regulatory sequences (Bickmore and van Steensel, 2013). Chromatin structure is dynamically shaped in a cell-specific manner, and the open chromatin regions allow TFs and other regulators of transcription to access cis-regulatory elements and activate gene expression, while closed chromatin regions impede the accessibility of promoters and enhancers to transcription machinery components leading to gene silencing (Klemm et al., 2019). Based on this principle, the genome-wide profiling of chromatin accessibility can be used to identify candidate regulatory elements associated with specific transcriptional and epigenetic signatures (Buenrostro et al., 2013). Recent advances in epigenomic profiling have enabled the development of ATAC-seq, which is a reliable tool for the generation of high-resolution chromatin accessibility landscape and direct recognition and sensitive detection of open chromatin regions (Buenrostro et al., 2013), both in bulk and at the single-cell level (Cao et al., 2018; Lareau et al., 2019). Depending on effective enzymatic cleavage and transposition of Tn5, being much faster and more sensitive, ATAC-seq has been widely applied in diverse development and disease-associated transitions to identify and localize epigenetic landscapes and regulators that drive disease pathogenesis (Corces et al., 2018; Wang et al., 2018; Xie et al., 2020).

In this study, we characterized the landscapes and relationship of chromatin accessibility and gene transcription in doxorubicin-resistant MCF7 cells, using ATAC-seq coupled with the analysis of RNA-seq data. The analysis of differentially accessible regions (DARs) allowed us to identify the regulatory DNA sequences, putative TFs, and landscape of TF binding events that may be responsible for these changes.

Materials and Methods

Cell Lines and Cell Culture

The MCF7 and doxorubicin-resistant MCF7 (MCF7-DR) cell lines originally obtained from the Cell Bank of the Type Culture Collection of the Chinese Academy of Sciences Committee (Shanghai, China), were kindly provided by the key laboratory of breast cancer in Shanghai, Fudan University. All cells were cultured in DMEM medium (KEL Biotech, China), supplemented with 2 mM L-glutamine (Gibco, United States), 10% fetal bovine serum (KEL Biotech, China), 0.1 M sodium pyruvate (Gibco, United States), 50 units/mL penicillin and 50 μg/mL streptomycin (Sangon Biotech, China). The MCF7 and MCF7-DR cells were grown in 10-cm culture dishes at 37°C in a humidified 5% CO₂ atmosphere, and the medium was renewed every 2–3 days. Whereas parental MCF7 cells were grown in a drug-free medium, resistant-derived cells were maintained in doxorubicin hydrochloride (MedChem Express, United States) at a concentration of 1 μM. Doxorubicin hydrochloride was stored at a concentration of 2 mM (2,000×) in DMSO (Sigma-Aldrich, United States) at −20°C.

Before any further experiments, the MCF7-DR cells were maintained in doxorubicin free medium for 3 days. Parental MCF7 cells were cultured in parallel with the MCF7-DR cells for comparison.

RNA Library Construction and Sequencing (RNA-seq)

Total RNA was isolated from MCF7 and MCF7-DR cells using TRIzol Reagent (Invitrogen, United States) and purified using RNeasy Mini Kit (QIAGEN, Germany). Agilent Bioanalyzer 2100 was used to assess RNA integrity. The RNA-seq library preparation from 1 μg of total RNA was performed with TruSeq RNA Sample Prep Kit (Illumina, United States) according to the manufacturer’s instructions. Finally, PCR products were purified to remove adaptor dimers with VAHTS DNA Clean Beads (Vazyme Biotech, China). Library quantity was assessed using a Qubit 2.0 Fluorometer (Invitrogen, United States). 150 bp pair-end sequencing was performed on the Illumina HiSeq X Ten platform (Novogene Biotech, China).

Assay for Transposase-Accessible Chromatin With High-Throughput Sequencing (ATAC-seq)

The ATAC-seq libraries were prepared as previously described with minor modifications (Buenrostro et al., 2013; Wu et al., 2016). Briefly, a total of 50,000 fresh cells were washed twice with 100 μl of cold PBS (KEL Biotech, China) and resuspended in 100 μl of lysis buffer [10 mM Tris–HCl (pH 7.4), 10 mM NaCl, 0.1% IGEPAL CA-630 (Sigma-Aldrich, United States) and 3 mM MgCl₂] for 10 min on ice to prepare the nuclei. Immediately after lysis, the suspension of nuclei was spun at 500 g at 4°C for 5 min to remove the supernatant. Nuclei were then incubated with the 50 μl of Tn5 transposition reaction mix (DIATRE Biotech, China) at 37°C for 30 min. The stop buffer was then added directly into the reaction to end the tagmentation. DNA was purified using the QIAquick PCR Purification Kit (QIAGEN, Germany). PCR was performed to amplify the library for 12 cycles using the following PCR conditions: 72°C for 3 min; 98°C for 30 s; and thermocycling at 98°C for 15 s, 60°C for 30 s and 72°C for 3 min; following by 72°C 5 min. Libraries were purified with VAHTS DNA Clean Beads (Vazyme Biotech, China) to remove contaminating primer dimmers. Library quantity was checked using Qubit 2.0 Fluorimeter (Invitrogen, United States). Last, all libraries were sequenced on the Illumina HiSeq X Ten with 150 bp paired-end reads (Novogene Biotech, China).

RNA-seq Data Processing

FASTQ files were evaluated for quality control using FastQC¹. The sequences were aligned to the human genome build hg19 using STAR (Dobin et al., 2013). The RNA expression of a gene was quantified by FPKM (fragments per kilobase million reads) based on RefSeq gene annotation using Cuffdiff (Trapnell et al., 2012). Differential expression genes (DEGs) were filtered by a q value < 0.05, the absolute value of fold change (FC) > 1.5 and the average FPKM > 10 in at least one of the two groups. The “ggplot2” (Wickham, 2009) in R and deepTools 2.0 (Ramirez et al., 2014) were used for data normalization and visualization. The aligned BAM files were converted to bigWig format and normalized by reads per kilobase million (RPKM) with the “bamCoverage” from deepTools 2.0 (Ramirez et al., 2014).

ATAC-seq Data Processing

The quality of sequenced ATAC-seq data was evaluated with FastQC (see text footnote 1). All sequencing data were mapped onto the hg19 human genome assembly using Bowtie2 (Langmead and Salzberg, 2012) with parameter –X 2000. SAMtools (Li et al., 2009) was used to remove duplicate reads, only non-duplicate reads kept in the BAM format were used for the subsequent analysis. Peak calling on nucleosome-free reads was performed by MACS2 (Zhang et al., 2008) using the module ‘‘callpeak’’ with the x following parameters --shift --100 --extsize 200. Picard Tools (v.2.2.4²) and ggplot2 (Wickham, 2009) were used to analyze and plot the distribution of paired-end sequencing fragment sizes. For data normalization and visualization, the BAM files were converted to the bigWig format using the “bamCoverage” (Ramirez et al., 2014) with RPKM. The heatmaps and average profiles were also generated using the “plotHeatmap” and “plotProfile” scripts in deepTools (Ramirez et al., 2014). To annotate the location of ATAC-Seq peaks in terms of important genomic features, we assigned their BED files to promoter-TSS (by default defined from −1 Kb to +100 bp of transcription start site), TTS (by default defined from −100 bp to +1 Kb of transcription termination site), intron, intergenic, exon, etc. using the annotatePeaks.pl script in HOMER (v4.10) (Heinz et al., 2010) with default parameters.

Identification of Differentially Accessible Regions

To assess the chromatin accessibility changes on a genome-wide scale, deepTools (Ramirez et al., 2014) was used to compute the average scores (RPKM normalized values) for transposase hypersensitive sites (THSs) in each sample. It was defined as hyper-accessible regions (hyper) associated with increased ATAC-seq signals if the regions showed the average fold change > 1.5 in MCF7-DR cells compared to parental MCF7 cells. On the contrary, the hypo-accessible regions (hypo) associated with decreased ATAC-seq signal was defined if the regions showed average fold change > 1.5 in parental MCF7 cells compared to MCF7-DR cells. The heatmaps and average profiles of regions expanded to ±1.5 Kb surrounding the DARs center were generated using the “plotHeatmap” and “plotProfile” functions in deepTools (Ramirez et al., 2014).

Correlation Analysis of DARs and DEGs

Annotation of DARs to genomic features was performed using the annotatePeaks.pl function in HOMER (Heinz et al., 2010). The DARs neighboring genes (distance to TSS < 100 Kb) were then overlapped with DEGs identified from RNA-seq. Fold changes of DARs and their nearest DEGs were used to calculate the Pearson’s correlation coefficient (PCC) and P value in R. Venn diagrams, boxplots and scatter plots were produced using “Venn Diagram” (Chen and Boutros, 2011) and “ggplot2” (Wickham, 2009) R packages.

Enrichment Analysis of GO Terms and KEGG Pathways

The clusterProfiler (Yu et al., 2012) was used to analyze GO terms and KEGG pathways, which is an R package for functional classification and enrichment of gene clusters using the hypergeometric distribution. The enriched terms and pathways with P value < 0.05 were considered significant. A Cytoscape network analysis³ was conducted to investigate links between significantly enriched top 20 KEGG pathways and key DEGs associated with MCF7-DR cells.

Transcription Factor Motif Enrichment and Occurrences Analysis

For the motif enrichment analysis of DARs, BEDTools (v.2.25.0) (Quinlan, 2014) was employed to obtain the peak summits overlapping with DARs, then we used findMotifsGenome (Heinz et al., 2010) to discover known motif surrounding ATAC-seq peak summits and to identify binding sites for a certain TF with the parameter –size −200,200. Motifs were only kept if the P value was < 0.01. The HOMER (v4.10) (Heinz et al., 2010) was used to identify the nearest genes of the top 50 TFs and analyze the occurrence probability of a certain TF using a 1 Kb window flanking the peak summits.

PPI Network Construction Using STRING

Based on RNA-seq and ATAC-seq data, differentially expressed transcription factors (DETFs) positively correlated with DARs were analyzed using the STRING database to identify groups of TFs that have known and predicted protein-protein associations (Szklarczyk et al., 2017). The network was subdivided into differentially colored nodes by the Markov cluster algorithm (MCL) with the default inflation parameter 3.0. The co-expression scores between proteins were calculated based on RNA expression patterns, and protein co-regulation provided by ProteomeHD (Kustatscher et al., 2019).

TCGA-BRCA Data Analysis

Based on the genomic and clinical data from the Cancer Genome Atlas Breast Invasive Carcinoma (TCGA-BRCA), the breast cancer samples were divided into two subgroups based on whether the Doxorubicin (a.k.a. Doxorubicin) was used alone, namely the doxorubicin alone group (Dox⁺; n = 300) and the other group (Dox^–; n = 762). The Kaplan–Meier survival curves and log-rank test were performed by utilizing the “survival” and “survminer” packages in R to analyze survival differences for critical DEGs. The mRNA expression levels of representative genes were compared between the Dox⁺ and Dox^– subgroups using the Wilcoxon signed rank test in R software (version 4.0.5). Pearson’s correlation coefficients and P values were also calculated to assess the associations between critical genes using the “Hmisc” R package. The P values less than 0.05 were regarded as statistically significant.

Statistical Analysis and Data Visualization

If not specified, R platform⁴ was used to compute statistics and generate plots throughout this manuscript. Student’s t-test, Pearson’s correlation coefficient, Wilcoxon signed rank test was used to assess the significance. Fold change, P values and false discovery rate (FDR) were calculated in analysis. Significant differences for all quantitative data were considered when ^∗P < 0.05, ^∗∗P < 0.01, ^∗∗∗P < 0.001, and ^****P < 0.0001. The genomic signal tracks were extracted and visualized using the WashU Epigenome Browser⁵.

Results

Characterization of Critical Genes, Functional and Pathway Enrichment Analyses Associated With Acquired Resistance to Doxorubicin

To investigate potential genes associated with the resistance to doxorubicin in MCF7 breast cancer cells, RNA-seq was performed with two biological replicates for each sample (Supplementary Figures 1A–C). A total of 3,963 differential expression genes (DEGs) were identified, among which 2,150 genes were up-regulated and 1,813 genes were down-regulated (Figure 1A, Supplementary Figures 1D,E, and Supplementary Tables 1, 2). Furthermore, to reveal crucial candidate genes associated with doxorubicin resistance, filtered by log₂(fold change) > 4, average FPKM > 150 in at least one group, we obtained 42 critical DEGs associated with doxorubicin resistance, including dramatically up-regulated metallothionein 1 family (MT1E, MT1M, and MT1L), GSTP1 (glutathione S-transferase pi 1), LDHB (lactate dehydrogenase B), and significantly down-regulated TFF1 (trefoil factor 1), UBB (ubiquitin B), DSCAM-AS1 (DSCAM Antisense RNA 1).

FIGURE 1

Figure 1. Identification of key genes and functional characteristics analysis. (A) Volcano plot showing the transcript levels of differentially expressed genes (DEGs) between the parental and doxorubicin-resistant MCF7 cells. Significantly up-regulated 2,150 genes are shown as red dots, whereas 1,813 down-regulated genes are shown as blue dots. Filtered by q value < 0.05, absolute log₂(FC) > 4 and average FPKM > 150, 42 significantly up- and down-regulated genes are highlighted in different colors according to log₂(FC). (B) Heatmap of 42 key genes labeled in the volcano plot (A), which are listed in a descending order based on log₂(FC). The numeric values are log₂(FPKM + 1). (C) Bubble plot showing GO enrichment analysis of DEGs. The top 20 GO terms of biological processes significantly enriched by up-regulated (left) and down-regulated (right) genes. (D) Bar plot showing top 20 KEGG pathways significantly enriched by up-regulated (left) and down-regulated (right) genes, respectively. The polygonal chain in black shows the gene count related to each GO term. (E) Venn diagram showing the overlap between the known histone-modifying enzymes and DEGs. (F) Bar plot indicating average expression levels (FPKM) of the 18 histone-modifying enzymes identified in panel (E), including histone acetyltransferase (purple), deacetylase (blue), methyltransferase (orange), demethylase (brown) and kinase (green). The P values were calculated using Student’s t test in R. *P < 0.05, **P < 0.01, and ***P < 0.001.

GO analysis was applied to elucidate biological processes, which indicated that the up-regulated DEGs were mainly involved in the mitotic cell cycle process, cytoskeleton organization, regulation of protein kinase activity. Meanwhile, the down-regulated DEGs were significantly enriched in epithelial cell differentiation, oxidation-reduction process, lipid biosynthetic process and Ras protein signal transduction (Figure 1C). In addition, the KEGG pathway analysis revealed that the up-regulated DEGs were mainly enriched in focal adhesion, cell cycle, Hippo and HIF-1 signaling pathways. The down-regulated DEGs were enriched in endocrine and platinum resistance, breast cancer, FoxO and p53 signaling pathways (Figure 1D).

Based on RNA-seq data, we identified 18 differentially expressed known histone-modifying enzymes, which were involved in epigenetic control of critical cellular functions and post-translational modification of histone and non-histone substrates (Stillman, 2018; Morgan and Shilatifard, 2020; Figure 1E). The results revealed that several genes coding for histone acetyltransferases such as ELP3, EP300 and KAT7, were significantly up-regulated; Except for EZH2 and PRMT5, most of the differentially expressed histone methyltransferase and demethylase genes were down-regulated (Figure 1F).

Genome-Wide Changes in Chromatin Accessibility

Assay for transposase-accessible chromatin with high-throughput sequencing (ATAC-seq) is the most current method for probing accessible chromatin regions, based on the ability of hyperactive Tn5 transposase to insert Illumina sequencing adaptors into active regulatory regions (Buenrostro et al., 2013; Wu et al., 2016; Li et al., 2017). To generate ATAC-seq libraries, we employed parental and doxorubicin-resistant MCF7 breast cancer cells with two independent biological replicates for each group (Supplementary Table 3). The chromatin was fragmented into mono-, di- and tri-nucleosome patterns and the similar distribution of fragment sizes indicates that chromatin is accessible to transposase to the same degree in all samples (Figure 2A). Genome-wide analysis of THSs distribution over functional genomic elements revealed that the promoter-TSS, introns and intergenic regions were preferentially (more than 90%) accessible to Tn5 transposase (Figure 2B). Moreover, our results conform to the pattern reported in previous studies in multiple organisms by showing that ATAC-seq signal was significantly enriched at ±1 Kb region surrounding TSS (Figure 2C), which suggests that the majority of cis-regulatory regions were located in the vicinity of gene core promoters and that open chromatin at THSs is predictive of active transcription (Buenrostro et al., 2013; Doganli et al., 2017; Lu et al., 2017; Bozek et al., 2019).

FIGURE 2

Figure 2. Genome-wide landscape of chromatin accessibility in MCF7 andMCF7-DR cells. (A) Distribution of ATAC-seq fragment size ineach sample. Clear modulation of signal is evident for mono-, di- and tri-nucleosomes. (B) Annotation of ATAC-seq peaks to genomic features: exon, intergenic regions, introns, 3′ UTR, 5′ UTR, promoters-TSS, TES and no-coding regions. Peak summit located within upstream –1,000 bp to + 100 bp downstream of the TSS is determined to be “promoter-TSS.” (C) The average read density (top) and heatmaps (bottom) indicating ATAC-seq signal across a genomic window of upstream –5 Kb to +5 Kb downstream of the TES. (D) Distribution of hyper- (red) and hypo-accessible (blue) regions over chromosomes in MCF7-DR cells compared to MCF7.

In order to determine whether the alteration of chromatin accessibility was consistent with transcription changes, we assessed differential accessibility between the parental and doxorubicin-resistant MCF7 breast cancer cells (MCF7-DR) across the genome. In total, 11,686 accessible regions were > 1.5-fold more accessible in MCF7-DR cells (Figures 2D, 3A–C). Conversely, 6,542 regions were > 1.5-fold more accessible in the MCF7 cells, and thus we defined these chromatin states as “hyper-accessible” and “hypo-accessible.” In general, these DARs were largely intragenic and gene distal (Figures 3A,D and Supplementary Table 4), with relatively few promoter-proximal regions exhibiting differential accessibility, suggesting that DARs might represent distal regulatory elements such as enhancers, which were often located far away from the genes they control. Moreover, it is worth noting that there were some wide continuous areas with decreased chromatin accessibility on chr7 and chr13 (Figure 2D), and most of the DEGs located in these regions were significantly down-regulated (Supplementary Table 2), suggesting that changes in chromatin accessibility are closely associated with the concomitant differential gene expression.

FIGURE 3

Figure 3. Identification and analysis of differentially accessible regions (DARs). (A) Distance to the closest transcription start site (TSS) of all accessible regions and newly open and closed regions. (B) Heatmap representation of normalized ATAC-seq signal in parental and doxorubicin-resistant MCF7 cells over DARs. The top panel shows read signal over the 11,686 hyper-accessible regions in MCF7-DR cells, while the bottom panel shows read signal over the 6,542 hypo-accessible regions. Signals within 1.5 Kb surrounding the center of DARs are displayed in descending order. (C) Profiles of normalized tag density across a genomic window of ±1.5 Kb surrounding the center of hyper- (top) and hypo-accessible (bottom) regions. (D) Pie chart showing the proportion of hyper- (top) and hypo-accessible (bottom) sites within the indicated genomic regions: introns, exon, intergenic regions, 3′ UTR, 5′ UTR, promoters-TSS, TES and no-coding regions. Peak summits located up to 1 Kb upstream and 100 bp downstream of the TSS are determined as promoter-TSS region. (E) Venn diagram illustrating the overlap among DEGs and the nearest genes of differentially accessible regions. (F) Box plots for mRNA expression levels of DEGs relative to hyper- (left) and hypo-accessible (right) regions. Notches of the boxes indicate medians. The P values were calculated by Wilcoxon’s signed-rank test. ****P < 0.0001. (G) Correlation analysis between DARs and their nearest DEGs. Each blue-colored dot represents a gene that is significantly differentially expressed and associated with chromatin accessibility changes. The top- and bottom-ranked 10 DEGs are labeled and shown in different colors based on the log₂(FC) of the average FPKM. The Pearson’s correlation coefficient (r) and the corresponding P value are shown.

The DARs Were Positively Correlated With Nearest DEGs

To acquire the potential links between DARs and DEGs, we assigned differentially accessible sites to the nearest genes according to their locations in the hg19 genome. Combining RNA-seq data, a total of 1,062 and 785 nearby DEGs were identified to be associated with the hyper- and hypo-accessible regions, respectively (Figure 3E). Interestingly, our results showed that most of the DEGs associated with hyper-accessible regions were significantly up-regulated, while the DEGs associated with hypo-accessible regions displayed significantly down-regulated RNA levels in all biological replicates (Figure 3F), indicating that genes within the open chromatin regions tended to have higher expression levels. Upon further analysis, it was confirmed that the changes in gene expression were positively associated with DARs located within 100 Kb of the TSS (Pearson correlation coefficient r = 0.6), suggesting that the alterations in chromatin accessibility contribute to differential expression of the genes associated with doxorubicin in MCF7-DR cells, such as MT1E, GPX1, LDHB and TPM2 (Figure 3G).

Key Signaling Pathway Associated DEGs That Positively Correlated With DARs

To identify potential key DEGs and signal pathways associated with DARs in MCF7-DR cells, we performed KEGG pathway enrichment analysis for DEGs associated with hyper- and hypo-accessible regions respectively. It was revealed that the hyper-accessible regions associated up-regulated DEGs were significantly enriched in small cell lung cancer, cell cycle, MAPK and Rap1 signaling pathways, while the hypo-accessible regions associated down-regulated DEGs mainly involved in endocrine resistance, lysosome, breast cancer, TGF-beta and p53 signaling pathways (Supplementary Figures 2, 4A). By integrating ATAC-seq and RNA-seq data, we obtained a total of 71 key DEGs enriched in the top 20 KEGG pathways with an absolute value of log₂(FC) > 3 (Figures 4A,B), including 41 up-regulated and 30 down-regulated DEGs positively associated with hyper- and hypo-accessible regions respectively, such as GPX1, FN1, FOSL1, CDH1, NOTCH3, WNT7B, and MT1E (Figures 4C,D and Supplementary Figure 3). This study indicated that global changes in chromatin accessibility affected the gene expression of many components of signaling cascades, and played a critical role in acquired resistance to doxorubicin in MCF7 breast cancer cells.

FIGURE 4

Figure 4. Key KEGG pathways associated DEGs related to DARs. (A) Top 20 KEGG pathways associated DEGs related to hyper- (left) and hypo-accessible (right) regions. DEGs are shown in different colors according to log₂(FC). (B) Bar plot showing average expression levels (FPKM) of DEGs related to hyper- (red shaded area) and hypo-accessible (blue shaded area) regions which enriched in top 20 KEGG pathways. DEGs are listed in a descending order based on log₂(FC). The P values were calculated using Student’s t test in R. *P < 0.05, **P < 0.01, ***P < 0.001, and ****P < 0.0001. (C) The WashU Epigenome Browser tracks show ATAC-seq (red) and mRNA-seq (blue) signals of representative up-regulated GPX1 and FN1 genes. The hyper-accessible regions are shaded with red. (D) Genomic snapshots of ATAC-seq (red) and mRNA-seq (blue) signal of representative down-regulated CDH1 and NOTCH3 genes. The hypo-accessible regions are shaded with blue.

TF Families Related to Doxorubicin-Induced Chromatin Accessibility Changes

TFs are molecules involved in regulating gene expression through binding to cis-regulatory specific sequences, known as motifs, in the promoters or enhancers of their target genes (Tsankov et al., 2015; Lambert et al., 2018). In most cases, chromatin accessibility to DNA is a prerequisite for TF binding (Huertas et al., 2020), which plays an important role in development and diseases (Whyte et al., 2013; Denny et al., 2016; Lambert et al., 2018; Tome-Garcia et al., 2018; Wang et al., 2018). To examine potential TFs located in the DARs in MCF7-DR cells, we performed TF enrichment analysis using motif discovery software HOMER (v4.10) (Heinz et al., 2010). As TFs binding to their cognate motifs, often obligates nucleosome eviction and creation of accessible chromatin regions, the integration of known TF motifs and chromatin accessibility data from ATAC-seq can predict TF occupancy on the chromatin across the whole genome (Buenrostro et al., 2013; Bao et al., 2015; Qu et al., 2017; Corces et al., 2018). The ± 200 bp flanking sequences around ATAC-seq peak summits in DARs were scanned for TF motif occurrences and TF binding sites (TFBSs). By default, 268 and 266 potential TFs were identified (P < 0.01) in hyper- and hypo-accessible regions, respectively (Supplementary Table 5). The activator protein-1 (AP-1) family members, including Jun (c-Jun, JunB, JunD), Fos (c-Fos, FosB, FosB2, Fra-1, and Fra-2), ATF (ATFa, ATF2, ATF3, BATF), and Maf (Maf1, MafA, MafK) proteins, were significantly enriched in both hyper- and hypo-accessible regions (Figures 5A,B), which always bind to the DNA sequence TGACTCA and play a critical role in various cellular events including proliferation, differentiation and apoptosis (Eferl and Wagner, 2003; Papoudou-Bai et al., 2017; Lambert et al., 2018). In addition, it was revealed that the hyper-accessible regions were also enriched with transcriptional enhanced associate domain (TEAD) family TFs, while the hypo-accessible regions were enriched for forkhead box (FOX) family members (Figures 5A,B and Supplementary Table 5).

FIGURE 5

Figure 5. Transcription factor (TF) motif enrichment analysis of DARs. (A,B) Top 20 enriched known TF motifs of hyper- (A) and hypo-accessible (B) sites, with P values estimated from HOMER v4.10. The AP-1, TEAD and FOX family members were colored red, green and blue, respectively. The polygonal chain in black shows the percentages of target sequences with TF motifs. (C,D) Putative TF occupancy of representative up- (C) and down-regulated (D) DEGs positively correlated to chromatin accessibility changes.

To obtain potential regulators of DEGs associated with doxorubicin resistance in breast cancer cells, we performed the analysis of the binding sites of top 50 TFs (Supplementary Table 5) enriched in DARs using HOMER (v4.10) (Heinz et al., 2010). The data results showed that the hyper-accessible sites relative to TSS of key up-regulated DEGs, including MT1E, LDHB, ANXA1, and TGFBl, were significantly enriched by TFBSs of KLF, ETV, RUNX and TEAD family members (Figure 5C). Meanwhile, the FOX family members and GRHL2 were mainly bound to the hypo-accessible sites associated with down-regulated DEGs, such as CDH1, GADD45G, LYPD3, and UBB (Figure 5D).

Differentially Expressed Transcription Factors (DETFs) and Their Target DEGs Related to DARs

In order to identify the DETFs that might be involved in the regulation of chromatin accessibility and gene transcription, the enriched motif cognate TFs (P < 10^–10) were intersected with the DEGs (Figures 6A,B). A total of 13 up-regulated TFs enriched in the hyper-accessible regions and 16 down-regulated TFs enriched in the hypo-accessible regions. For instance, besides AP-1 family members, TEAD1 and RUNX1 were significantly up-regulated and enriched in hyper-accessible regions, while down-regulation of GRHL2 and FOXA1 might be responsible for global decreases in chromatin accessibility. As expected, further analysis demonstrated that TEAD1 and RUNX1 showed higher binding probability around peak summits in hyper-accessible regions, while GRHL2 and FOXA1 were generally distributed more frequently in hypo-accessible regions (Figure 6C). Moreover, the differential expression patterns of these TFs were positively correlated with the alterations in chromatin accessibility (Figure 6D). The majority of predicted target genes of TEAD1 and RUNX1 were significantly up-regulated, including SFTA1P and LY6K, while most of GRHL2 and FOXA1 targets were down-regulated, such as SLC44A2 and GPNMB (Figure 6E).

FIGURE 6

Figure 6. Differentially expressed TFs (DETFs) and their target DEGs. (A) Venn diagram (left) showing overlap of the corresponding TFs for each motif in hyper-accessible regions and up-regulated DEGs. The dot plot (right) showing the identified up-regulated TFs. The color of each dot represents the P value of enrichment (calculated using the cumulative binomial distribution in HOMER v4.10) for each motif and cell type, and the size of each dot represents the gene expression (FPKM) of the enriched motif cognate TFs. Only TFs motif enrichment P value < 10^{– 10} were included. (B) Venn diagram (left) showing overlap of the corresponding TFs for each motif in hypo-accessible regions and down-regulated DEGs. The dot plot (right) showing the identified down-regulated TFs. (C) Distribution probability of TEAD1, RUNX1, GRHL2 and FOXA1 binding motifs around ATAC-seq peak summits within DARs. (D) The WashU Epigenome Browser views show the hyper-accessible chromatin regions (red shaded area) near two up-regulated TFs (TEAD1 and RUNX1), and the hypo-accessible chromatin regions (blue shaded area) near two down-regulated TFs (GRHL2 and FOXA1). (E) Bar charts showing the expression levels of potential TF (TEAD1, RUNX1, GRHL2 and FOXA1) target DEGs in DARs predicted using HOMER. Only DEGs with TFBSs in the promoter-TSS region were included. Genes are listed in a descending order based on log₂(FC). TEAD1 and RUNX1 common targets were colored in red, while GRHL2 and FOXA1 common targets were colored in blue. The P values were calculated using Student’s t test in R. *P < 0.05, **P < 0.01, ***P < 0.001, and ****P < 0.0001.

Loss of the Pioneer Factor FOXA1 Expression Might Be Responsible for Doxorubicin Resistance in Breast Cancer

It has been recently reported that FOXA1/GRHL2 collaboration could establish a targetable pathway in endocrine therapy-resistant breast cancer (Cocce et al., 2019), and facilitate depositing H3K4me1 at potential enhancer elements by recruiting the chromatin modifier MLL3 (Jozwik et al., 2016). GATA3 and FOXA1 are marker genes that define an ESR1-positive breast cancer (Cancer Genome Atlas, 2012). GATA3 acts upstream of FOXA1, loss of which results in changes in FOXA1 binding and alteration of gene expression (Theodorou et al., 2013). Our data demonstrated that a large proportion of the predicted target DEGs were potentially co-regulated by both FOXA1 and GRHL2 in MCF7 cells (Figure 6E), such as SLC44A2, GPNMB, BATF, and GOLT1A, which suggested that loss of FOXA1/GRHL2 collaboration might be responsible for the down-regulation of critical genes in doxorubicin-resistant breast cancer cells. To identify DETFs (P < 10^–2) that might act coordinately to control their gene expression, we constructed a protein-protein interaction (PPI) network using the STRING database (Szklarczyk et al., 2017), which serves as an entry point for a better interpretation of relationships between different TFs on a genome-wide scale. It revealed that clustered SMAD4/KLF5/HIF1A and FOXA1/GRHL2/ELF5 potentially worked together to regulate gene expression in hyper- and hypo-accessible regions, respectively (Figure 7A). Moreover, there were higher observed co-expression scores for FOXM1/E2F1, FOXA1/GRHL2, FOXA1/GATA3 (Figure 7B), suggesting that absence of FOXA1 might contribute to doxorubicin resistance in breast cancer cells through collaboration with the down-regulated GRHL2 and/or GATA3.

FIGURE 7

Figure 7. Protein-protein interaction (PPI) and co-expression network analysis of DETFs. (A) PPI network showing DETFs identified by HOMER enriched in hyper- (left) and hypo-accessible (right) regions. Network nodes indicate proteins, and the colored edges represent evidence of different types of protein-protein interactions. Strong relationships and networks were clustered using Markov cluster algorithm by default inflation parameters and indicated by solid lines. (B) Observed co-expression in Homo sapiens identified by STRING for DETFs associated with DARs.

To further testify the expression patterns and the clinical significance of critical DEGs related to chromatin accessibility changes in MCF7-DR cells as well as predict their role in the development of acquired resistance to doxorubicin in breast cancer, we retrieved the genomic and clinical data from the Cancer Genome Atlas Breast Invasive Carcinoma (TCGA-BRCA) and analyzed expression levels of the DEGs and their impact on the prognosis of breast cancer patients treated with (Dox⁺; n = 300) and without (Dox^–; n = 762) doxorubicin alone. Remarkably, the results showed that lower expression of several significantly down-regulated genes in MCF7-DR cells was verified to be significantly (P < 0.05) associated with the poor prognosis of the Dox⁺ patients when compared to the Dox^– patients, including FOXA1, GATA3, RARA, and SLC27A2 (Figure 8A and Supplementary Figure 4). On the contrary, higher expression of up-regulated GSTP1 and LDHB was related with the poor prognosis of Dox⁺ patients (Figure 8A). The expression levels of FOXA1 and GATA3 were significantly (P < 0.05) lower in the Dox⁺ patients than Dox^–, while GSTP1 and LDHB expression levels were higher in the Dox⁺ patients (Figure 8B). Similarly, the expression patterns of RARA and SLC27A2 were generally consistent with the previously described transcriptomic changes in MCF7-DR cells (Figure 1B, Supplementary Figure 5, and Supplementary Table 2). Moreover, in the Dox⁺ patients, there was a strong positive correlation between the mRNA expression levels of FOXA1 and GATA3 (r = 0.827), GSTP1 and LDHB (r = 0.614), FOXA1 and SPDEF (r = 0.857), while FOXA1 expression was negatively correlated with GSTP1 (r = −0.653) and LDHB (r = −0.7) as expected (Figure 8C and Supplementary Figure 6). These results revealed that some of the above-mentioned DEGs associated with chromatin accessibility changes, such as FOXA1, GATA3 and GSTP1, played a critical role in the development of doxorubicin resistance in breast cancer cells, and might be important indicators for judging the prognosis or a therapeutic target in breast cancer patients treated with doxorubicin.

FIGURE 8

Figure 8. Integrative analysis of clinical and genomic data from TCGA-BRCA datasets. (A) Kaplan–Meier survival analyses for FOXA1, GATA3, GSTP1 and LDHB of breast invasive carcinoma (BRCA) patients treated with (Dox⁺; top panel) and without (Dox^–; bottom panel) doxorubicin (adriamycin). Expression values were sorted according to high and low expression (FPKM) from the median of the corresponding genes. The shaded area indicates the 95% confidence interval. All the survival duration was censored to 10 years. Statistical analysis was performed by the log-rank test, and P < 0.05 was statistically significant. (B) Box plots showing the mRNA expression levels (FPKM) of FOXA1, GATA3, GSTP1 and LDHB in Dox⁺ and Dox^– breast cancer subgroups derived from TCGA datasets. The P values were calculated by Wilcoxon’s signed-rank test. *P < 0.05, **P < 0.01, and ***P < 0.001. (C) Correlation analysis of mRNA expression levels among FOXA1, GATA3, GSTP1 and LDHB in the Dox⁺ breast cancer subgroup. The Pearson’s correlation coefficients (r) and P values were calculated using R packages. The shaded gray area around the red line represents the 95% confidence interval. The P value less than 0.05 was considered to be statistically significant.

Discussion

Despite a few molecular mechanisms involved in the process of doxorubicin resistance in cancer have been reported (Huang et al., 2018; Lovitt et al., 2018; Stefanski et al., 2019; Jin et al., 2020; Wang et al., 2020; Mirzaei et al., 2021; Yang et al., 2021), drug resistance is almost inevitable in advanced breast cancer patients, which remains a major unresolved issue in the treatment of cancer patients. An increasing number of studies have revealed that the alteration in chromatin accessibility plays a critical role in acquired resistance against cancer drugs (Jing et al., 2018; Zhang et al., 2020), but we still lack precise knowledge of how accessible chromatin state define transcriptional regulatory networks and which are the regulatory elements involved in doxorubicin resistance in breast cancer cells.

In the present study, we integrated chromatin accessibility and transcription analysis in parental and doxorubicin-resistant MCF7 cells via ATAC-Seq in combination with RNA-seq, which are mutually authenticating. Compared to existing FAIRE-seq, DNase-seq and MNase-seq technologies, ATAC-seq offers substantial advantages and becomes the most widely used method to analyze the genome-wide profile of chromatin accessibility due to its low input cell number requirement, simplicity, speediness and fine accuracy (Buenrostro et al., 2013).

Based on the analysis of differential expression level of mRNA, we identified several dramatically up-regulated genes associated with doxorubicin resistance, such as glutathione S-transferases P1(GSTP1), metallothionein gene family (MT1E, MT1M, and MT1L) and glutathione peroxidase 1 (GPX1). It has been reported that overexpression of metallothionein confers resistance to anti-cancer drugs (Kelley et al., 1988); GPX1 promotes cisplatin resistance via ROS-Induced activation of PI3K/AKT pathway in non-small cell lung cancer (Chen et al., 2019). Up-regulation of GSTP1 maintains resistance to doxorubicin through inhibition of PI3K/AKT/mTOR activity to promote autophagy in breast cancer cells (Dong et al., 2019); Simultaneously, there were a lot of significantly down-regulated genes in MCF7-DR cells, including trefoil factor 1 (TFF1), ubiquitin B (UBB), and SAM pointed domain containing ETS TF (SPDEF). Furthermore, from a perspective of epigenetic regulation, we analyzed the changes in the expression of known histone-modifying enzymes especially. The results suggested that the up-regulation of histone deacetylase 2 (HDAC2) and histone methyltransferase EZH2 was potentially involved in the down-regulation of certain doxorubicin-resistant genes.

In general, our RNA-seq results are broadly consistent with previous studies regarding the identification of DEGs and pathway enrichment analysis. However, as far as we know, most of the identified DEGs have not yet been studied in sufficient detail to identify their genetic functions involved in doxorubicin-resistant breast cancer.

Genome-wide analysis of chromatin accessibility revealed that ATAC-Seq signal was significantly enriched in promoter-TSS and distal regulatory elements such as enhancers. Combined with RNA-seq data, a significant positive correlation (r = 0.6) was found between DARs and DEGs, validating our expectation that increased chromatin accessibility results in higher expression of critical genes associated with doxorubicin resistance, such as MT1E, GPX1, and TPM2. In the same way, condensed chromatin, known as closed chromatin, restricts the binding of TFs and transcriptional regulators to the promoter-TSS and/or enhancer, which leads to gene silencing, including TFF1, UBB and CDH1. Furthermore, we screened the DEGs that were significantly positively correlated with chromatin accessibility changes, and then, Cytoscape (version 3.7.1) software was used to display their enriched functional pathways.

Transcription factors play an important role in cancer drug resistance (Fan et al., 2017; Lambert et al., 2018; Huh et al., 2019; Bi et al., 2020). By analyzing the accessibility of chromatin, we cannot only identify the regulatory regions but also infer TF motifs. Therefore, we performed known motif analysis of the DARs to identify potential TFs and their target DEGs using the motif discovery software HOMER. We found that, besides the AP-1 family, members of TEAD and FOX families were significantly enriched in hyper- and hypo-accessible regions in MCF7-DR cells. The AP-1 is another family of bZIP TFs and comprises multiple proteins that bind as heterodimers (Lambert et al., 2018), which is also a common TF family in motif analysis of ATAC-Seq data. In this study, we would like to pay more attention to other TFs with important regulatory functions. Finally, we identified several significantly up- and down-regulated TFs (P < 10^–10) associated with the hyper- and hypo-accessible regions, respectively, such as TEAD1, RUNX1, GRHL2 and FOXA1. It has been previously reported that the pioneers TFs FOXA1 and GATA3 mediated ESR1 binding by shaping enhancer accessibility (Theodorou et al., 2013); FOXA1 collaborates with GRHL2 and MLL3 to establish a targetable collateral pathway in endocrine therapy-resistant breast cancer (Werner et al., 2013; Cocce et al., 2019). The PPI network analysis showed that GATA3/FOXA1/GRHL2 exhibited strong correlations, suggesting loss of GATA3/FOXA1/GRHL2 expression might play an important role in the down-regulation of key genes in MCF7-DR cells, such as SLC44A2, BATF, and PRR15L. Moreover, the analysis of genomic and clinical data from TCGA-BRCA further revealed that low expression levels of FOXA1 and GATA3 in patients treated with doxorubicin alone experienced significantly worse clinical outcomes, while GSTP1 and LDHB had opposite prognostic effects. As was expected, all of them exhibited similar gene expression patterns in both DOX⁺ breast cancer patients and MCF7-DR cells.

Most of the previous research on drug resistance has focused on the up-regulated genes but paid relatively little attention to the down-regulated genes. Our findings suggest that down-regulation or loss of certain regulators may also play an important role in promoting the development of cancer drug resistance. According to previous studies and our data analysis, we concluded that chromatin accessibility changes could promote or inhibit the expression of critical genes associated with doxorubicin resistance by affecting the DNA-binding ability of TFs and transcriptional regulators in promoter and/or enhancer. It is probably that, in doxorubicin-resistant breast cancer cells, loss of the pioneer TF FOXA1 leads to a significant decrease in chromatin openness, which impairs the ability of co-regulators GATA3 and/or GRHL2 to occupy their specific binding sites in the promoter and/or enhancer region of specific genes, such as BATF, SLC44A2, and GPNMB. Besides, the opposite changes in gene expression of FOXA1 and GSTP1 might be potential therapeutic targets for overcoming or reversing doxorubicin resistance in breast cancer.

Conclusion

The genome-wide profiling of open chromatin by ATAC-seq represents a novel strategy to study accessibility dynamics and transcriptional regulation in cancer drug resistance. The application of this technique allows us to identify the accessible regulatory elements and TFs that are probably responsible for the acquired resistance to doxorubicin in breast cancer cells. Our ATAC-seq data contribute to close the loop for the association between TF binding, transcriptional and chromatin states, providing new insights into non-genetic mechanisms of doxorubicin resistance and putative targets for epigenetic therapy.

Data Availability Statement

All relevant datasets supporting the key findings of this study are available within the article and its Supplementary Material or from the corresponding author on reasonable request. All sequencing data have been deposited in the GEO repository: GSE174152.

Author Contributions

GW and XW conceived the idea. XW prepared samples and generated the data. XW analyzed the data with assistance from GW and JY. XW wrote the manuscript. GW and BS supervised the study and revised the manuscript. JY provided helpful comments on this study. All authors reviewed and approved the final manuscript.

Funding

This research was funded by National Key Research and Development Plan (2016YFA0100703 and 2015CB964800) and the National Natural Science Foundation of China, grant numbers 31771431 and 81871906.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We would like to thank the experimental support of the Uli Schwarz public laboratory platform in SINH, CAS. We would also like to thank and appreciate Yuanyuan Zhao, Zhao Ye, Yujie Chen, and other colleagues for their technical assistance and valuable suggestions for this study.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fcell.2021.708066/full#supplementary-material

Footnotes

1. ^ http://www.bioinformatics.babraham.ac.uk/projects/fastqc/
2. ^ https://broadinstitute.github.io/picard/
3. ^ http://www.cytoscape.org/
4. ^ https://www.R-project.org/
5. ^ https://epgg-test.wustl.edu/browser/

References

Bao, X., Rubin, A. J., Qu, K., Zhang, J., Giresi, P. G., Chang, H. Y., et al. (2015). A novel ATAC-seq approach reveals lineage-specific reinforcement of the open chromatin landscape via cooperation between BAF and p63. Genome Biol. 16:284.

Google Scholar

Bi, M., Zhang, Z., Jiang, Y. Z., Xue, P., Wang, H., Lai, Z., et al. (2020). Enhancer reprogramming driven by high-order assemblies of transcription factors promotes phenotypic plasticity and breast cancer endocrine resistance. Nat. Cell Biol. 22, 701–715. doi: 10.1038/s41556-020-0514-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Bickmore, W. A., and van Steensel, B. (2013). Genome architecture: domain organization of interphase chromosomes. Cell 152, 1270–1284. doi: 10.1016/j.cell.2013.02.001

PubMed Abstract | CrossRef Full Text | Google Scholar

Bozek, M., Cortini, R., Storti, A. E., Unnerstall, U., Gaul, U., and Gompel, N. (2019). ATAC-seq reveals regional differences in enhancer accessibility during the establishment of spatial coordinates in the Drosophila blastoderm. Genome Res. 29, 771–783. doi: 10.1101/gr.242362.118

PubMed Abstract | CrossRef Full Text | Google Scholar

Buenrostro, J. D., Giresi, P. G., Zaba, L. C., Chang, H. Y., and Greenleaf, W. J. (2013). Transposition of native chromatin for fast and sensitive epigenomic profiling of open chromatin. DNA-binding proteins and nucleosome position. Nat. Methods. 10, 1213–1218. doi: 10.1038/nmeth.2688

PubMed Abstract | CrossRef Full Text | Google Scholar

Burgess, D. J., Doles, J., Zender, L., Xue, W., Ma, B., McCombie, W. R., et al. (2008). Topoisomerase levels determine chemotherapy response in vitro and in vivo. Proc. Natl. Acad. Sci. U.S.A. 105, 9053–9058. doi: 10.1073/pnas.0803513105

PubMed Abstract | CrossRef Full Text | Google Scholar

Cancer Genome Atlas, N. (2012). Comprehensive molecular portraits of human breast tumours. Nature 490, 61–70. doi: 10.1038/nature11412

PubMed Abstract | CrossRef Full Text | Google Scholar

Cao, J. Y., Cusanovich, D. A., Ramani, V., Aghamirzaie, D., Pliner, H. A., Hill, A. J., et al. (2018). Joint profiling of chromatin accessibility and gene expression in thousands of single cells. Science 361, 1380–1385. doi: 10.1126/science.aau0730

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, B., Shen, Z., Wu, D., Xie, X., Xu, X., Lv, L., et al. (2019). Glutathione peroxidase 1 promotes NSCLC resistance to cisplatin via ROS-induced activation of PI3K/AKT pathway. Biomed Res. Int. 2019:7640547.

Google Scholar

Chen, H., and Boutros, P. C. (2011). VennDiagram: a package for the generation of highly-customizable Venn and Euler diagrams in R. BMC Bioinformatics 12:35. doi: 10.1186/1471-2105-12-35

PubMed Abstract | CrossRef Full Text | Google Scholar

Ciocan-Cartita, C. A., Jurj, A., Zanoaga, O., Cojocneanu, R., Pop, L. A., Moldovan, A., et al. (2020). New insights in gene expression alteration as effect of doxorubicin drug resistance in triple negative breast cancer cells. J. Exp. Clin. Cancer Res. 39:241

Google Scholar

Cocce, K. J., Jasper, J. S., Desautels, T. K., Everett, L., Wardell, S., Westerling, T., et al. (2019). The lineage determining factor GRHL2 collaborates with FOXA1 to establish a targetable pathway in endocrine therapy-resistant breast cancer. Cell Rep. 29, 889–903.e10.

Google Scholar

Corces, M. R., Granja, J. M., Shams, S., Louie, B. H., Seoane, J. A., Zhou, W., et al. (2018). The chromatin accessibility landscape of primary human cancers. Science 362:eaav1898.

Google Scholar

Denny, S. K., Yang, D., Chuang, C. H., Brady, J. J., Lim, J. S., Gruner, B. M., et al. (2016). Nfib promotes metastasis through a widespread increase in chromatin accessibility. Cell 166, 328–342. doi: 10.1016/j.cell.2016.05.052

PubMed Abstract | CrossRef Full Text | Google Scholar

Dobin, A., Davis, C. A., Schlesinger, F., Drenkow, J., Zaleski, C., Jha, S., et al. (2013). STAR: ultrafast universal RNA-seq aligner. Bioinformatics 29, 15–21. doi: 10.1093/bioinformatics/bts635

PubMed Abstract | CrossRef Full Text | Google Scholar

Doganli, C., Sandoval, M., Thomas, S., and Hart, D. (2017). Assay for transposase-accessible chromatin with high-throughput sequencing (ATAC-Seq) protocol for zebrafish embryos. Methods Mol. Biol. 1507, 59–66. doi: 10.1007/978-1-4939-6518-2_5

CrossRef Full Text | Google Scholar

Dong, X., Yang, Y., Zhou, Y., Bi, X., Zhao, N., Zhang, Z., et al. (2019). Glutathione S-transferases P1 protects breast cancer cell from adriamycin-induced cell death through promoting autophagy. Cell Death Differ. 26, 2086–2099. doi: 10.1038/s41418-019-0276-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Eferl, R., and Wagner, E. F. (2003). AP-1: a double-edged sword in tumorigenesis. Nat. Rev. Cancer 3, 859–868. doi: 10.1038/nrc1209

PubMed Abstract | CrossRef Full Text | Google Scholar

Fan, F., Bashari, M. H., Morelli, E., Tonon, G., Malvestiti, S., Vallet, S., et al. (2017). The AP-1 transcription factor JunB is essential for multiple myeloma cell proliferation and drug resistance in the bone marrow microenvironment. Leukemia 31, 1570–1581. doi: 10.1038/leu.2016.358

PubMed Abstract | CrossRef Full Text | Google Scholar

Hammond, W. A., Swaika, A., and Mody, K. (2016). Pharmacologic resistance in colorectal cancer: a review. Ther. Adv. Med. Oncol. 8, 57–84. doi: 10.1177/1758834015614530

PubMed Abstract | CrossRef Full Text | Google Scholar

Heinz, S., Benner, C., Spann, N., Bertolino, E., Lin, Y. C., Laslo, P., et al. (2010). Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. Mol. Cell 38, 576–589. doi: 10.1016/j.molcel.2010.05.004

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, J. F., Wen, C. J., Zhao, G. Z., Dai, Y., Li, Y., Wu, L. X., et al. (2018). Overexpression of ABCB4 contributes to acquired doxorubicin resistance in breast cancer cells in vitro. Cancer Chemother. Pharm. 82, 199–210. doi: 10.1007/s00280-018-3603-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Huertas, J., MacCarthy, C. M., Scholer, H. R., and Cojocaru, V. (2020). Nucleosomal DNA dynamics mediate oct4 pioneer factor binding. Biophys. J. 118, 2280–2296. doi: 10.1016/j.bpj.2019.12.038

PubMed Abstract | CrossRef Full Text | Google Scholar

Huh, H. D., Kim, D. H., Jeong, H. S., and Park, H. W. (2019). Regulation of TEAD transcription factors in cancer biology. Cells 8:600. doi: 10.3390/cells8060600

PubMed Abstract | CrossRef Full Text | Google Scholar

Jin, Y., Wang, H., Zhu, Y., Feng, H., Wang, G., and Wang, S. (2020). miR-199a-5p is involved in doxorubicin resistance of non-small cell lung cancer (NSCLC) cells. Eur. J. Pharmacol. 878:173105. doi: 10.1016/j.ejphar.2020.173105

PubMed Abstract | CrossRef Full Text | Google Scholar

Jing, D., Huang, Y., Liu, X., Sia, K. C. S., Zhang, J. C., Tai, X., et al. (2018). Lymphocyte-specific chromatin accessibility pre-determines glucocorticoid resistance in acute lymphoblastic leukemia. Cancer Cell 34, 906–921.e8.

Google Scholar

Jozwik, K. M., Chernukhin, I., Serandour, A. A., Nagarajan, S., and Carroll, J. S. (2016). FOXA1 directs H3K4 monomethylation at enhancers via recruitment of the methyltransferase MLL3. Cell Rep. 17, 2715–2723. doi: 10.1016/j.celrep.2016.11.028

PubMed Abstract | CrossRef Full Text | Google Scholar

Kelley, S. L., Basu, A., Teicher, B. A., Hacker, M. P., Hamer, D. H., and Lazo, J. S. (1988). Overexpression of metallothionein confers resistance to anticancer drugs. Science 241, 1813–1815. doi: 10.1126/science.3175622

PubMed Abstract | CrossRef Full Text | Google Scholar

Klemm, S. L., Shipony, Z., and Greenleaf, W. J. (2019). Chromatin accessibility and the regulatory epigenome. Nat. Rev. Genet. 20, 207–220. doi: 10.1038/s41576-018-0089-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Kustatscher, G., Grabowski, P., Schrader, T. A., Passmore, J. B., Schrader, M., and Rappsilber, J. (2019). Co-regulation map of the human proteome enables identification of protein functions. Nat. Biotechnol. 37, 1361–1371. doi: 10.1038/s41587-019-0298-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Lambert, S. A., Jolma, A., Campitelli, L. F., Das, P. K., Yin, Y., Albu, M., et al. (2018). The human transcription factors. Cell 175, 598–599.

Google Scholar

Langmead, B., and Salzberg, S. L. (2012). Fast gapped-read alignment with Bowtie 2. Nat. Methods 9, 357–359. doi: 10.1038/nmeth.1923

PubMed Abstract | CrossRef Full Text | Google Scholar

Lareau, C. A., Duarte, F. M., Chew, J. G., Kartha, V. K., Burkett, Z. D., Kohlway, A. S., et al. (2019). Droplet-based combinatorial indexing for massive-scale single-cell chromatin accessibility. Nat. Biotechnol. 37, 916–924. doi: 10.1038/s41587-019-0147-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, D., Liu, J., Yang, X., Zhou, C., Guo, J., Wu, C., et al. (2017). Chromatin accessibility dynamics during iPSC reprogramming. Cell Stem Cell. 21, 819–833.e6.

Google Scholar

Li, H., Handsaker, B., Wysoker, A., Fennell, T., Ruan, J., Homer, N., et al. (2009). The sequence alignment/map format and SAMtools. Bioinformatics 25, 2078–2079. doi: 10.1093/bioinformatics/btp352

PubMed Abstract | CrossRef Full Text | Google Scholar

Lovitt, C. J., Shelper, T. B., and Avery, V. M. (2018). Doxorubicin resistance in breast cancer cells is mediated by extracellular matrix proteins. Bmc Cancer 18:41. doi: 10.1186/s12885-017-3953-6

PubMed Abstract | CrossRef Full Text | Google Scholar

Lu, Z., Hofmeister, B. T., Vollmers, C., DuBois, R. M., and Schmitz, R. J. (2017). Combining ATAC-seq with nuclei sorting for discovery of cis-regulatory regions in plant genomes. Nucleic Acids Res. 45:e41. doi: 10.1093/nar/gkw1179

PubMed Abstract | CrossRef Full Text | Google Scholar

Meredith, A. M., and Dass, C. R. (2016). Increasing role of the cancer chemotherapeutic doxorubicin in cellular metabolism. J. Pharm. Pharmacol. 68, 729–741. doi: 10.1111/jphp.12539

PubMed Abstract | CrossRef Full Text | Google Scholar

Mirzaei, S., Zarrabi, A., Hashemi, F., Zabolian, A., Saleki, H., Azami, N., et al. (2021). Nrf2 signaling pathway in chemoprotection and doxorubicin resistance: potential application in drug discovery. Antioxidants (Basel) 10:349. doi: 10.3390/antiox10030349

PubMed Abstract | CrossRef Full Text | Google Scholar

Morgan, M. A. J., and Shilatifard, A. (2020). Reevaluating the roles of histone-modifying enzymes and their associated chromatin modifications in transcriptional regulation. Nat. Genet. 52, 1271–1281. doi: 10.1038/s41588-020-00736-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Papoudou-Bai, A., Hatzimichael, E., Barbouti, A., and Kanavaros, P. (2017). Expression patterns of the activator protein-1 (AP-1) family members in lymphoid neoplasms. Clin. Exp. Med. 17, 291–304. doi: 10.1007/s10238-016-0436-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Qu, K., Zaba, L. C., Satpathy, A. T., Giresi, P. G., Li, R., Jin, Y., et al. (2017). Chromatin accessibility landscape of cutaneous T cell lymphoma and dynamic response to HDAC inhibitors. Cancer Cell 32, 27–41.e4.

Google Scholar

Quinlan, A. R. (2014). BEDTools: the swiss-army tool for genome feature analysis. Curr. Protoc. Bioinformatics. 47, 11.12.1–34.

Google Scholar

Ramirez, F., Dundar, F., Diehl, S., Gruning, B. A., and Manke, T. (2014). deepTools: a flexible platform for exploring deep-sequencing data. Nucleic Acids Res. 42, W187–W191.

Google Scholar

Stefanski, C. D., Keffler, K., McClintock, S., Milac, L., and Prosperi, J. R. (2019). APC loss affects DNA damage repair causing doxorubicin resistance in breast cancer cells. Neoplasia 21, 1143–1150. doi: 10.1016/j.neo.2019.09.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Stillman, B. (2018). Histone modifications: insights into their influence on gene expression. Cell 175, 6–9. doi: 10.1016/j.cell.2018.08.032

PubMed Abstract | CrossRef Full Text | Google Scholar

Sung, H., Ferlay, J., Siegel, R. L., Laversanne, M., Soerjomataram, I., Jemal, A., et al. (2021). Global cancer statistics 2020: GLOBOCAN estimates of incidence and mortality worldwide for 36 cancers in 185 countries. CA Cancer J. Clin. 71, 209–249. doi: 10.3322/caac.21660

PubMed Abstract | CrossRef Full Text | Google Scholar

Szklarczyk, D., Morris, J. H., Cook, H., Kuhn, M., Wyder, S., Simonovic, M., et al. (2017). The STRING database in 2017: quality-controlled protein-protein association networks, made broadly accessible. Nucleic Acids Res. 45, D362–D368.

Google Scholar

Tanaka, I., Chakraborty, A., Saulnier, O., Benoit-Pilven, C., Vacher, S., Labiod, D., et al. (2020). ZRANB2 and SYF2-mediated splicing programs converging on ECT2 are involved in breast cancer cell resistance to doxorubicin. Nucleic Acids Res. 48, 2676–2693. doi: 10.1093/nar/gkz1213

PubMed Abstract | CrossRef Full Text | Google Scholar

Theodorou, V., Stark, R., Menon, S., and Carroll, J. S. (2013). GATA3 acts upstream of FOXA1 in mediating ESR1 binding by shaping enhancer accessibility. Genome Res. 23, 12–22. doi: 10.1101/gr.139469.112

PubMed Abstract | CrossRef Full Text | Google Scholar

Tome-Garcia, J., Erfani, P., Nudelman, G., Tsankov, A. M., Katsyv, I., Tejero, R., et al. (2018). Analysis of chromatin accessibility uncovers TEAD1 as a regulator of migration in human glioblastoma. Nat. Commun. 9:4020.

Google Scholar

Tormo, E., Ballester, S., Adam-Artigues, A., Burgues, O., Alonso, E., Bermejo, B., et al. (2019). The miRNA-449 family mediates doxorubicin resistance in triple-negative breast cancer by regulating cell cycle factors. Sci. Rep. 9:5316.

Google Scholar

Trapnell, C., Roberts, A., Goff, L. A., Pertea, G., Kim, D., Kelley, D. R., et al. (2012). Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nat. Protoc. 7, 562–578. doi: 10.1038/nprot.2012.016

PubMed Abstract | CrossRef Full Text | Google Scholar

Tsankov, A. M., Gu, H., Akopian, V., Ziller, M. J., Donaghey, J., Amit, I., et al. (2015). Transcription factor binding dynamics during human ES cell differentiation. Nature 518, 344–349. doi: 10.1038/nature14233

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, D., Diao, H., Getzler, A. J., Rogal, W., Frederick, M. A., Milner, J., et al. (2018). The transcription factor Runx3 establishes chromatin accessibility of cis-regulatory landscapes that drive memory cytotoxic T lymphocyte formation. Immunity 48, 659–674.e6.

Google Scholar

Wang, F., Xiang, Z., Huang, T., Zhang, M., and Zhou, W. B. (2020). ANLN directly interacts with RhoA to promote doxorubicin resistance in breast cancer cells. Cancer Manag. Res. 12, 9725–9734. doi: 10.2147/cmar.s261828

PubMed Abstract | CrossRef Full Text | Google Scholar

Werner, S., Frey, S., Riethdorf, S., Schulze, C., Alawi, M., Kling, L., et al. (2013). Dual roles of the transcription factor grainyhead-like 2 (GRHL2) in breast cancer. J. Biol. Chem. 288, 22993–23008. doi: 10.1074/jbc.m113.456293

PubMed Abstract | CrossRef Full Text | Google Scholar

Whyte, W. A., Orlando, D. A., Hnisz, D., Abraham, B. J., Lin, C. Y., Kagey, M. H., et al. (2013). Master transcription factors and mediator establish super-enhancers at key cell identity genes. Cell 153, 307–319. doi: 10.1016/j.cell.2013.03.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Wickham, H. (2009). Ggplot2: Elegant Graphics for Data Analysis. New York, NY: Springer Publishing Company.

Google Scholar

Wu, J., Huang, B., Chen, H., Yin, Q., Liu, Y., Xiang, Y., et al. (2016). The landscape of accessible chromatin in mammalian preimplantation embryos. Nature 534, 652–657. doi: 10.1038/nature18606

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, H., Zhang, W., Zhang, M., Akhtar, T., Li, Y., Yi, W., et al. (2020). Chromatin accessibility analysis reveals regulatory dynamics of developing human retina and hiPSC-derived retinal organoids. Sci Adv. 6:eaay5247. doi: 10.1126/sciadv.aay5247

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, B., Ke, W., Wan, Y., and Li, T. (2021). Targeting RNF8 effectively reverses cisplatin and doxorubicin resistance in endometrial cancer. Biochem. Biophys. Res. Commun. 545, 89–97. doi: 10.1016/j.bbrc.2021.01.046

PubMed Abstract | CrossRef Full Text | Google Scholar

Yu, G., Wang, L. G., Han, Y., and He, Q. Y. (2012). clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 16, 284–287. doi: 10.1089/omi.2011.0118

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, L. L., Lu, J., Liu, R. Q., Hu, M. J., Zhao, Y. M., Tan, S., et al. (2020). Chromatin accessibility analysis reveals that TFAP2A promotes angiogenesis in acquired resistance to anlotinib in lung cancer cells. Acta Pharmacol. Sin. 41, 1357–1365. doi: 10.1038/s41401-020-0421-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Liu, T., Meyer, C. A., Eeckhoute, J., Johnson, D. S., Bernstein, B. E., et al. (2008). Model-based analysis of ChIP-Seq (MACS). Genome Biol. 9:R137.

Google Scholar