Identification of genes and pathways associated with MDR in MCF-7/MDR breast cancer cells by RNA-seq analysis.

Multidrug resistance (MDR) is a major problem in the treatment of breast cancer. In the present study, next-generation sequencing technology was employed to identify differentially expressed genes in MCF‑7/MDR cells and MCF‑7 cells, and aimed to investigate the underlying molecular mechanisms of MDR in breast cancer. Differentially expressed genes between MCF‑7/MDR and MCF‑7 cells were selected using software; a total of 2085 genes were screened as differentially expressed in MCF‑7/MDR cells. Furthermore, gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed using the DAVID database. Finally, a protein‑protein interaction network was constructed and the hub genes in the network were analyzed using the STRING database. GO annotation demonstrated that the differentially expressed genes were enriched in various biological processes, including 'regulation of cell differentiation', 'cell development', 'neuron development', 'movement of cell or subcellular component' and 'cell morphogenesis involved in neuron differentiation'. Cellular component analysis by GO revealed that differentially expressed genes were enriched in 'plasma membrane region' and 'extracellular matrix' terms. Furthermore, KEGG analysis demonstrated that the target genes were enriched in various pathways, including 'cell adhesion molecules (CAMs)', 'calcium signaling pathway', 'tight junction', 'Wnt signaling pathway' and 'pathways in cancer' terms. A protein‑protein interaction network demonstrated that certain hub genes, including cyclin D1, nitric oxide synthase 3 (NOS3), NOTCH3, brain‑derived neurotrophic factor (BDNF), paired box 6, neuropeptide Y, phospholipase C β (PLCB) 4, PLCB2 and actin α cardiac muscle 1, may be associated with MDR in breast cancer. Subsequently, RT‑qPCR confirmed that the expression of these 9 hub genes was higher in MCF‑7/MDR cells compared with MCF‑7 cells, consistent with the RNA‑sequencing analysis. Additionally, a Cell Counting Kit‑8 assay demonstrated that specific inhibitors of NOS3 and BDNF/neurotrophic receptor tyrosine kinase, type 2 signaling reduced the IC50 of MCF‑7/MDR cells in response to various anticancer drugs, including adriamycin, cisplatin and 5‑fluorouracil. The results of the present study provide novel insights into the mechanism underlying MDR in MCF‑7 cells and may identify novel targets for the treatment of breast cancer.

Introduction

Globally, breast cancer is the most commonly occurring type of cancer in women. Currently, despite substantial advances in treatment, multidrug resistance (MDR) remains a major clinical obstacle to effective treatment and results in a poor prognosis for these patients. MDR is a phenomenon that occurs when cancer cells exposed to anticancer drugs develop cross-resistance to a range of structurally and functionally dissimilar agents.

Previous studies have demonstrated that numerous mechanisms are involved in the regulation of MDR. The most recognized genes associated with MDR are the ATP-binding cassette (ABC) transporter family. ABC transporters are members of a large protein family that are involved in the import or export of an extensive spectrum of substrates, and increasing evidence has confirmed the role of ABC proteins, such as ABCB1, ABCC1 and ABCG2, in an MDR breast cancer cell phenotype (1–3). Altered drug metabolism is also a major mechanism that contributes to MDR. Glutathione S-transferases (GSTs) consist of a family of phase II metabolic isozymes that facilitate the addition of reduced glutathione to xenobiotic substrates for detoxification, and overexpression of GSTs has been implicated in breast cancer MDR (4). In addition, defects in apoptotic machinery are reported to lead to abnormal chemoresistance, abnormal expression of the Bcl-2 family has been implicated in breast cancer drug resistance and various anticancer drugs, such as imatinib, gefitinib, and bortezomib, are regarded as Bcl-2-like 11-targeting therapeutic agents (5,6). Furthermore, caspase-3, a member of the caspase family, is activated by various chemotherapeutic agents, including doxorubicin and etoposide (7,8), and has an important role in apoptotic pathways through its cleavage of cellular proteins. In addition, tumor suppressor p53 and Raf-1 proto-oncogene have also been reported to influence the response of breast cancer to doxorubicin treatment (9,10). Homeostasis disruption also influences chemoresistance in breast cancer (11), and interleukin-6, estrogen and nomegestrol acetate are reported to affect MDR in breast cancer (12–14).

miRNAs are a novel class of small, noncoding RNAs that regulate the expression of various target genes through RNA silencing and post-transcriptional gene regulation to regulate the activity of cells. Certain miRNAs have been reported to be associated with drug resistance in breast cancer; for example, miR-326 affected the chemoresistance of breast cancer through modulation of multidrug resistance-associated protein 1 (15), miR-221/222 inhibited tamoxifen resistance by targeting p27 (16), miR-125b suppressed Bcl-2 antagonist killer 1 expression and inhibited resistance in response to paclitaxel (17), and miR-34a was associated with docetaxel resistance in breast cancer (18). Increasing evidence has indicated that miRNAs have a key regulatory role in MDR through modulation of various mechanisms, thereby providing insights for the development of more effective individualized therapies for breast cancer treatment.

Numerous signaling pathways are also considered to be determinants of MDR in breast cancer, including phosphatidylinositol-3-kinase/Akt/mechanistic target of rapamycin kinase, Janus kinase/signal transducer and activator of transcription, AMP-activated protein kinase and epithelial-mesenchymal transition (EMT) (19–23). Although numerous genes and signaling pathways are reported to be associated with MDR in cancer, the exact mechanism of MDR remains unclear. RNA-sequencing (RNA-seq) is a recently developed method for transcriptome profiling that employs next-generation sequencing technologies and measures transcript levels with increased precision compared with other approaches (24). Recently, RNA-seq has been employed extensively to investigate mechanisms of drug resistance in cancer; differentially expressed genes screened by RNA-seq provide novel insight into the complex mechanisms of resistance to anticancer drugs (25,26).

To investigate the mechanism of MDR in breast cancer cells in vitro, drug resistant cell lines have been isolated by exposing cancer cell lines to increasing concentrations of various chemotherapeutic agents. MCF-7, a human breast carcinoma cell line that has been employed in studies longer than any other breast cancer cell model system, has been widely employed in the investigation of apoptosis, MDR and a malignant phenotype of breast cancer. MCF-7/MDR breast cancer cells have been reported to exhibit resistance to several chemotherapeutic agents for cancer, including doxorubicin, paclitaxel, salvianolic acid A, cisplatin, etoposide and melphalan (27–29). In the current study, RNA-seq was performed to analyze the differential expression of genes between MCF-7 and MCF-7/MDR human breast cancer cells to determine the underlying molecular mechanisms of chemoresistance in breast cancer.

Materials and methods

Cell culture and reagents

The MCF-7 human breast cancer cell line was obtained from the Cell Bank of Type Culture Collection of the Chinese Academy of Sciences (Shanghai, China). The multidrug-resistant human breast cancer cell MCF-7/MDR was purchased from the Cell Bank of Xiangya Medical College, Central South University (Changsha, China). Cells were cultured in Dulbecco's modified Eagle's medium (DMEM; Gibco; Thermo Fisher Scientific, Inc., Waltham, MA, USA) supplemented with 10% (v/v) fetal bovine serum (FBS; Thermo Fisher Scientific, Inc.), 100 U/ml penicillin (Gibco; Thermo Fisher Scientific, Inc.), 100 U/ml streptomycin (Gibco; Thermo Fisher Scientific, Inc.) and 10 mmol/l HEPES at 37°C in a humidified atmosphere containing 5% CO2. Adriamycin (ADM), 5-fluorouracil (5-FU), cisplatin (DDP) and L-N5-(1-Iminoethyl)-ornithine hydrochloride (LNIO) were purchased from Sigma-Aldrich (Merck KGaA, Darmstadt, Germany). ANA-12 was purchased from MedChemExpress (Monmouth Junction, NJ, USA).

RNA-seq

Total RNA was isolated from MCF-7 and MCF-7/MDR cells using TRIzol reagent (Invitrogen; Thermo Fisher Scientific, Inc.), according to the manufacturer's protocol. First-strand reverse transcription (RT) primer: biotin, 5′-CAAGCAGAAGACGGCATACGAGTVN-3′ was conducted using a SuperScript III First-Strand Synthesis system (Thermo Fisher Scientific, Inc.) (30) was used to synthesize the first-strand cDNA according to the manufacturer's protocols, followed by purification of cDNA with a NucleoSpin gel and polymerase chain reaction (PCR) Clean-up kit (Macherey-Nagel GmbH & Co. KG, Düren, Germany) to remove free RT primer. Dynabeads® MyOne™ Streptavidin C1 (cat. no. 650.01; Thermo Fisher Scientific, Inc.) were used for further purification. The second-strand cDNA was synthesized using dNTP mix and Taq DNA polymerase with ThermoPol reaction Buffer (cat. no. M0237L; New England BioLabs, Inc., Ipswich, MA, USA) following the program: 25°C for 60 min, 68°C for 30 sec, and 75°C for 5 min. Subsequently, the second-strand cDNA was amplified by PCR using the following program: 10 min at 94°C for initial denaturation, 35 cycles of 30 sec at 94°C for denaturation, 30 sec at 58°C for annealing, and 30 sec at 72°C for elongation. Following purification, the pooled libraries were loaded in a lane of an 8-lane flow cell and sequenced with the HumanHT-12 v4 Expression BeadChip arrays (Illumina, Inc., Sand Diego, CA, USA) on a HiSeq 2500 Illumina sequencer (Illumina, Inc.). RNA-Seq libraries were generated and 12 libraries were pooled for multiplex (30,31).

Sequencing data analysis

The initial sequence reads were generated by the Illumina Genome Analyzer (Illumina, Inc.). The expression values were calculated using RSEM (32) for each gene (Ensembl, Homo.sapiens.GRCh38.87.gtf) in individual samples. A total of 1 million reads were randomly selected for quality control with Bowtie 2 (33) mapping to Genome Reference Consortium Human Build 38/hg38 (34). Differential expression was analyzed between different subgroups using Bioconductor edgeR (version 3.12.0) (35). A false discovery rate of 0.05 was applied to the resulting P-values to correct for multiple hypothesis testing. T-tests were employed to screen out differentially expressed genes between MCF-7/MDR and MCF-7 cells, with a threshold of P<0.05 and fold change >2.

Gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway mapping analysis

The Database for Annotation, Visualization and Integrated Discovery (DAVID) database (36) was used to analyze the differentially expressed genes between MCF-7/MDR and MCF-7 cells. The potential targets of differentially expressed genes were analyzed by the GO program and enriched pathways were analyzed by the KEGG program. The GO analysis included biological process and cellular component terms, and P<0.01 was regarded as statistically significant. In the KEGG pathway enrichment analysis, enriched pathways were identified according to P<0.05.

RT-quantitative PCR (RT-qPCR)

Total RNA was isolated from 5×106 cells using TRIzol reagent, according to the manufacturer's protocol. Total RNA concentration and purity were analyzed in duplicate samples using a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Inc.). cDNA was synthesized from the RNA using the TransScriptR Reverse Transcriptase kit (Beijing TransGen Biotech Co., Ltd., Beijing, China). The reaction program constituted: 65°C for 5 sec, 2 min on ice, 25°C for 5 min, 42°C for 30 min, and 85°C for 5 sec. qPCR assays were performed to evaluate gene expression profiles using an ABI PRISM 7300 Sequence Detection system. qPCR was performed with the TransStart® Top Green qPCR SuperMix kit with SYBR® Green I from Beijing Transgen Biotech Co., Ltd. with specified primers. The reaction program consisted of 95°C for 3 min followed by 40 cycles of 95°C for 30 sec, 55°C for 20 sec and 72°C for 15 sec. GAPDH served as an internal control and the relative mRNA levels were calculated by the 2−ΔΔCq method (37). Primers were designed and synthesized by Sangon Biotech Co., Ltd., (Shanghai, China), as listed in Table I.

Table I.

Primers for reverse transcription-quantitative polymerase chain reaction.

	Primer sequence (5′-3′)
Gene	Forward	Reverse
BDNF	CCAAGGCAGGTTCAAGAGG	TCCAGCAGAAAGAGAAGAGGA
GAPDH	ACCCACTCCTCCACCTTTG	CTCTTGTGCTCTTGCTGGG
PLCB2	AAGGTGAAGGCCTATCTGAGCCAA	CTTGGCAAACTTCCCAAAGCGAGT
PLCB4	GCACAGCACACAAAGGAATGGTCA	CGCATTTCCTTGCTTTCCCTGTCA
ACTC1	GCCCTGGATTTTGAGAATGA	ATGCCAGCAGATTCCATACC
NOS3	GTGGCTGGTACATGAGCACT	GTGGTCCACGATGGTGACTT
CCND1	GCTGCGAAGTGGAAACCATC	CCTCCTTCTGCACACATTTGAA
NPY	TGCTAGGTAACAAGCGACTG	CTGCATGCATTGGTAGGATG
NOTCH3	TCTCAGACTGGTCCGAATCCAC	CCAAGATCTAAGAACTGACGAGCG
PAX6	CGAGACTGGCTCCATCAGAC	CACTCCCGCTTATACTGGGC
TrkB	TGGTGCATTCCATTCACTGT	CGTGGTACTCCGTGTGATTG

BDNF, brained-derived neurotrophic factor; PLCB2, phospholipase C β2; PLCB4, phospholipase C β4; ACTC1, actin alpha cardiac muscle 1; NOS3, nitric oxide synthase 3; CCND1, cyclin D1; NPY, neuropeptide Y; PAX6, paired box 6; TrkB, neurotrophic receptor tyrosine kinase, type 2.

Western blot analysis

A total of 2×105 MCF-7/MDR cells/well were plated in 6-well culture plates to determine the optimal concentration of LNIO in subsequent experiments, cells were treated with 0, 5, 10, 20, 40 and 80 µM LNIO for 48 h at 37°C and western blotting was performed. Western blotting was also performed to analyze the difference between the protein expression of NOS3 in MCF-7 and MCF-7/MDR cells. Total protein was extracted from MCF-7 and MCF-7/MDR cells with 200 µl radioimmunoprecipitation assay lysis buffer (Beyotime Institute of Biotechnology, Haimen, China) containing 1 mM phenylmethylsulfonyl fluoride (Sigma-Aldrich; Merck KGaA), followed by centrifugation at 13,000 × g for 20 min at 4°C. Denatured protein (50 µg) was quantified via a Bicinchoninic Acid protein assay kit (Beyotime Institute of Biotechnology) was loaded onto 12% SDS polyacrylamide gels. Following electrophoresis, the proteins on the gel were transferred onto nitrocellulose membranes (EMD Millipore, Billerica, MA, USA). The membranes were blocked with 5% non-fat milk at 37°C for 1 h and incubated overnight at 4°C with the following primary antibodies: Monoclonal mouse β-actin antibody (1:1,000; cat. no. sc-130300; Santa Cruz Biotechnology, Inc., Dallas, TX, USA) and polyclonal rabbit NOS3 antibody (1:500; cat. no. wl01789; Wanleibo Co., Ltd., Shanghai, China). The membranes were subsequently incubated with anti-mouse IgG (1:2,000; cat. no. 7076; Cell Signaling Technology, Inc., Danvers, MA, USA) and anti-rabbit IgG (1:2,000; cat. no. 7074; Cell Signaling Technology, Inc.) horseradish peroxidase-conjugated secondary antibodies at room temperature for 1 h. Finally, the immunoreactive bands of NOS3 and GAPDH were visualized using a BeyoECL Star kit (Beyotime Institute of Biotechnology). Semi-quantitative analysis was conducted using Quantity-One software version 4.1 (Bio-Rad Laboratories, Inc., Hercules, CA, USA) to measure densitometric values for each band.

In vitro drug sensitivity assay

In vitro drug cytotoxicity was measured by a Cell Counting Kit-8 (CCK-8; Dojindo Molecular Technologies, Inc., Kumamoto, Japan) assay. MCF-7/MDR cells were seeded into 96-well plates (3,000 cells/well) and treated for 48 h at 37°C in 100 µl DMEM medium containing 10% FBS and various concentrationsof ADM, DDP and 5-FU (ADM: 0, 5, 10, 20, 40 and 80 µM; DDP and 5-FU: 0, 0.5, 1, 2, 4 and 8 µg/ml) with or without pretreatment with 10 µM LNIO or 25 µM ANA-12 for 24 h at 37°C. The cells incubated without drugs were set at 100% survival and were used to calculate the concentration of each cytostatic drug that was lethal to 50% of the cells (IC50). CCK-8 reagent was subsequently added and the cells were incubated at 37°C for 1 h. The optical density (OD) of each well at 450 nm was recorded on a Varioskan Flash microplate reader (Thermo Fisher Scientific, Inc.). The cell viability (% of control) was expressed as the percentage of (ODtest-ODblank)/(ODcontrol-ODblank).

Statistical analysis

All experiments were performed in triplicate and data are presented as the mean ± standard deviation. Statistical analysis was performed with SPSS 19.0 software (SPSS, Inc., Chicago, IL, USA). Statistical significance was determined by a Student's t-test and one-way analysis of variance, followed by a Dunnett's test. P<0.05 was considered to indicate a statistically significant difference.

Results

Identification of differentially expressed genes between MCF-7 and MCF-7/MDR cells

edgeR was used to analyze the gene expression profiles and identify the differentially expressed genes between MCF-7/MDR and MCF-7 cells with the described criteria (38). A total of 56,233 genes were detected and 2,085 genes exhibited differential expression in MCF-7 compared with MCF-7/MDR cells. A fold change >2 and P<0.05 were required for differences to be considered as statistically significant. A total of 100 top differentially expressed genes, including 50 upregulated genes and 50 downregulated genes between MCF-7/MDR cells and MCF-7 cells were observed (Fig. 1).

Figure 1.

Heat map of the top 100 differentially expressed genes (50 upregulated genes and 50 downregulated genes) between MCF-7/MDR and MCF-7 cells. MDR, multidrug resistance; Red, upregulation; Blue, downregulation.

GO term enrichment analysis

Using the DAVID database, functional analysis of the differentially expressed genes was performed. In the GO analysis, biological process and cellular component terms were included. GO analysis of biological processes demonstrated that differentially expressed genes were enriched in ‘neurogenesis’ (60 genes), ‘neuron differentiation’ (53 genes) and ‘nervous system development’ (84 genes) biological process terms, with others enriched in ‘locomotion’ (57 genes), ‘cell motility’ (49 genes), ‘localization of cell’ (49 genes) and ‘movement of cell or subcellular component’ (67 genes). In addition, differentially expressed genes were also enriched in ‘cell-cell signaling’ (54 genes) and ‘regulation of extent of cell growth’ (10 genes) biological process terms, among others (Fig. 2 and Table II). Furthermore, GO analysis of cellular components demonstrated that the differentially expressed genes were primarily associated with ‘plasma membrane region’ (42 genes), ‘proteinaceous extracellular matrix’ (22 genes), ‘extracellular matrix’ (27 genes), ‘synapse’ (35 genes), ‘synapse part’ (30 genes) and ‘apical plasma membrane’ (18 genes; Fig. 3 and Table III).

Figure 2.

GO analysis of the enrichment of differentially expressed genes between MCF-7/MDR and MCF-7 cells in biological process terms. GO, gene ontology; MDR, multidrug resistance.

Table II.

GO biological process analysis of differentially expressed genes.

Term	Count	Genes	FDR
GO:0007399~nervous system development	84	DLC1, CTHRC1, HMX2, TUBB2B, MYO7A, USH1G, UCHL1, JAG2, PAX6, ANKRD1, FGF12, RTN1, ARHGAP4, PCDH1, GSTM3, BDNF, SPINT2, PAX8, SEMA3E, SNPH, SEMA3D, POU4F1, NRG1, INA, KIF17, CLMN, EFNB1, EFNB2, DLL3, SIX3, GPER1, HMGA2, PCLO, HOXD9, EPB41L3, PRKCQ, ADM, SIX1, FOXG1, MYRF, NRGN, EXT1, GAP43, NEURL1, DPF3, FGFR1, WNT5B, NDN, TH, BEX1, MYEF2, OXTR, CLDN11, FKBP1B, MDK, ITM2C, CDH4, EPHB1, ITM2A, ALCAM, BCL11A, GBX2, APBA2, NKX2.1, NEFH, ADRA2C, FN1, FZD8, TRNP1, PTPRZ1, GABRA5, NTNG1, AXL, NTN1, NKX6.1, NOTCH3, SALL2, HOXB2, NPY, PLK2, H2AFY2, GFRA1, SCN8A, APBB1	0.0000335
GO:0048699~generation of neurons	58	CTHRC1, TUBB2B, MYO7A, USH1G, UCHL1, JAG2, PAX6, ANKRD1, RTN1, ARHGAP4, BDNF, SNPH, SEMA3E, SEMA3D, POU4F1, NRG1, CLMN, EFNB1, EFNB2, SIX3, DLL3, GPER1, HMGA2, HOXD9, EPB41L3, PRKCQ, ADM, SIX1, EXT1, GAP43, NEURL1, FGFR1, WNT5B, NDN, TH, MYEF2, FKBP1B, ITM2C, CDH4, EPHB1, ALCAM, BCL11A, GBX2, NKX2-1, NEFH, ADRA2C, FN1, FZD8, PTPRZ1, GABRA5, AXL, NTNG1, NTN1, NKX6-1, NOTCH3, NPY, PLK2, APBB1	0.0004872
GO:0022008~neurogenesis	60	CTHRC1, TUBB2B, MYO7A, USH1G, UCHL1, JAG2, PAX6, ANKRD1, RTN1, ARHGAP4, BDNF, SNPH, SEMA3E, SEMA3D, POU4F1, NRG1, CLMN, KIF17, EFNB1, EFNB2, DLL3, SIX3, GPER1, HMGA2, HOXD9, EPB41L3, PRKCQ, ADM, SIX1, MYRF, EXT1, GAP43, NEURL1, FGFR1, WNT5B, NDN, TH, MYEF2, FKBP1B, ITM2C, CDH4, EPHB1, ALCAM, BCL11A, GBX2, NKX2-1, NEFH, ADRA2C, FN1, FZD8, PTPRZ1, GABRA5, AXL, NTNG1, NTN1, NKX6-1, NOTCH3, NPY, PLK2, APBB1	0.0007804
GO:0007409~axonogenesis	27	NDN, UCHL1, PAX6, CDH4, EPHB1, ARHGAP4, ALCAM, BDNF, BCL11A, SEMA3E, GBX2, SEMA3D, NEFH, NKX2-1, POU4F1, NRG1, FN1, PTPRZ1, EFNB1, EFNB2, NTNG1, NTN1, NKX6-1, PRKCQ, EXT1, APBB1, GAP43	0.0011163
GO:0030182~neuron differentiation	53	CTHRC1, MYO7A, USH1G, UCHL1, PAX6, JAG2, ANKRD1, RTN1, ARHGAP4, BDNF, SNPH, SEMA3E, SEMA3D, POU4F1, NRG1, CLMN, EFNB1, EFNB2, SIX3, HOXD9, PRKCQ, EPB41L3, ADM, SIX1, EXT1, GAP43, NEURL1, FGFR1, WNT5B, NDN, TH, MYEF2, ITM2C, CDH4, FKBP1B, EPHB1, ALCAM, BCL11A, GBX2, NKX2-1, NEFH, ADRA2C, FN1, FZD8, PTPRZ1, GABRA5, NTNG1, NTN1, NKX6-1, NOTCH3, NPY, PLK2, APBB1	0.0013166
GO:0061564~axon development	28	NDN, UCHL1, PAX6, FKBP1B, CDH4, EPHB1, ARHGAP4, ALCAM, BDNF, BCL11A, SEMA3E, GBX2, SEMA3D, NEFH, NKX2-1, POU4F1, NRG1, FN1, PTPRZ1, EFNB1, EFNB2, NTNG1, NTN1, NKX6-1, PRKCQ, EXT1, APBB1, GAP43	0.0013888
GO:0046903~secretion development	49	TACR2, ANKRD1, KISS1, PAX8, SNPH, CHRNA7, NOS2, NRG1, PTPRN2, ACTN1, SERPING1, GPER1, CD40, PCDH7, HMGA2, PTPRN, PCLO, INHBB, CCND1, ADM, STXBP6, RAB15, NEURL1, XDH, FGFR1, IL1R2, RAB3B, NKD2, WASH1, HLA-DRB1, LMF1, OXTR, FKBP1B, AHSG, C1QTNF5, HLA-DRB5, ADRA2C, OLFM2, FN1, LPL, AXL, IGF2, FOXP3, NKX6-1, ABCB4, FBLN5, APBB1, SYTL1, SLC17A9	0.0016024
GO:0048468~cell development	73	CTHRC1, MAEA, MYO7A, USH1G, UCHL1, PAX6, ANKRD1, YBX2, ARHGAP4, GSTM3, BDNF, FRMD6, SPINT2, PAX8, SEMA3E, SEMA3D, POU4F1, NRG1, RAMP2, CLMN, EFNB1, EFNB2, DLL3, SIX3, SIP A1L3, ACTN1, GPER1, HMGA2, FLNC, HOXD9, INHBB, EPB41L3, PRKCQ, REC8, ADM, SIX1, EPB41L5, MYRF, RYR1, EXT1, LRRK1, GAP43, NEURL1, FGFR1, WNT5B, NDN, TH, FHL2, FKBP1B, ITM2C, CDH4, EPHB1, ALCAM, BCL11A, GBX2, NKX3.2, NKX2.1, NEFH, ADRA2C, FN1, ACTC1, PTPRZ1, GABRA5, NTNG1, AXL, NTN1, NKX6.1, NOTCH3, NPY, PLK2, KLF2, APBB1, RFLNB	0.0021369
GO:0051241~negative regulation of multicellular organismal process	46	XDH, IL1R2, FGFR1, HLA-DRB1, TACR2, PAX6, OXTR, FKBP1B, ITM2C, AHSG, ARHGAP4, BCL11A, PAX8, SEMA3E, NKX3-2, SEMA3D, HLA-DRB5, NKX2-1, CHRNA7, NOS3, ADRA2C, ADAMTS12, NRG1, ARHGDIB, FN1, F12, DLL3, AXL, SIX3, UBE2L6, SERPING1, FOXP3, HMGA2, NTN1, NKX6-1, NOTCH3, INHBB, ADRB2, TNNT1, PLK2, ADM, SIX1, SERPINB2, KLF2, TMEM176B, RFLNB	0.0026439
GO:0006928~movement of cell or subcellular component	67	DLC1, CTHRC1, TUBB2B, MLPH, TNNC1, MYO7A, UCHL1, JAG2, PAX6, FGF12, AMOTL1, MMP1, ARHGAP4, BDNF, FRMD6, SPINT2, SEMA3E, SEMA3D, NOS3, CHRNA7, POU4F1, NRG1, ARC, KIF17, PDPN, EFNB1, ARHGEF5, EFNB2, SIX3, ACTN1, GPER1, SIRPA, PRKCQ, TNNT1, SIX1, EPB41L5, EXT1, CTSH, GAP43, NEURL1, FGFR1, WNT5B, APC2, WASH1, NDN, ABI3, PLVAP, MDK, CDH4, EPHB1, ALCAM, IFT46, GBX2, NKX2-1, NEFH, ADAMTS12, FN1, ARHGDIB, ACTC1, FLT4, AXL, NTN1, NKX6-1, LAMA1, NPY, PTP4A3, HOXB9	0.0056175
GO:0000904~cell morphogenesis involved in differentiation	37	NDN, MYO7A, UCHL1, PAX6, CDH4, EPHB1, ARHGAP4, ALCAM, BDNF, FRMD6, SPINT2, BCL11A, PAX8, SEMA3E, GBX2, SEMA3D, NKX2-1, NEFH, POU4F1, NRG1, FN1, PTPRZ1, EFNB1, EFNB2, NTNG1, AXL, SIPA1L3, ACTN1, HMGA2, NTN1, NKX6-1, PRKCQ, EPB41L5, EXT1, APBB1, GAP43, RFLNB	0.0059248
GO:0048666~neuron development	43	FGFR1, CTHRC1, NDN, MYO7A, USH1G, UCHL1, TH, PAX6, ANKRD1, FKBP1B, ITM2C, CDH4, EPHB1, ARHGAP4, ALCAM, BDNF, BCL11A, SEMA3E, GBX2, SEMA3D, NKX2-1, NEFH, POU4F1, NRG1, FN1, PTPRZ1, CLMN, EFNB1, GABRA5, EFNB2, NTNG1, NTN1, NKX6-1, HOXD9, EPB41L3, PRKCQ, NPY, PLK2, ADM, EXT1, APBB1, GAP43, NEURL1	0.0093174
GO:0032940~secretion by cell	43	IL1R2, FGFR1, RAB3B, NKD2, HLA-DRB1, WASH1, TACR2, LMF1, OXTR, ANKRD1, FKBP1B, AHSG, KISS1, C1QTNF5, PAX8, SNPH, HLA-DRB5, CHRNA7, ADRA2C, NOS2, OLFM2, FN1, LPL, PTPRN2, AXL, ACTN1, SERPING1, IGF2, PCDH7, CD40, GPER1, PTPRN, FOXP3, HMGA2, PCLO, NKX6-1, INHBB, ADM, STXBP6, RAB15, SYTL1, APBB1, SLC17A9	0.0094679
GO:0048667~cell morphogenesis involved in neuron differentiation	28	NDN, MYO7A, UCHL1, PAX6, CDH4, EPHB1, ARHGAP4, ALCAM, BDNF, BCL11A, SEMA3E, GBX2, SEMA3D, NEFH, NKX2-1, POU4F1, NRG1, FN1, PTPRZ1, EFNB1, EFNB2, NTNG1, NTN1, NKX6-1, PRKCQ, EXT1, APBB1, GAP43	0.0107608
GO:0031175~neuron projection development	37	FGFR1, NDN, UCHL1, PAX6, ANKRD1, FKBP1B, ITM2C, CDH4, EPHB1, ARHGAP4, ALCAM, BDNF, BCL11A, SEMA3E, GBX2, SEMA3D, NEFH, NKX2-1, POU4F1, NRG1, FN1, PTPRZ1, CLMN, EFNB1, EFNB2, NTNG1, NTN1, NKX6-1, EPB41L3, PRKCQ, NPY, PLK2, ADM, EXT1, APBB1, GAP43, NEURL1	0.0355076
GO:0040011~locomotion	57	DLC1, CTHRC1, TUBB2B, JAG2, PAX6, AMOTL1, MMP1, ARHGAP4, BDNF, SPINT2, CCR10, SEMA3E, SEMA3D, CHRNA7, POU4F1, NOS3, NRG1, FOSL1, ARC, PDPN, EFNB1, ARHGEF5, EFNB2, SIX3, GPER1, SIRPA, PRKCQ, SIX1, EPB41L5, EXT1, CTSH, GAP43, NEURL1, FGFR1, WNT5B, WASH1, NDN, APC2, ABI3, PLVAP, MDK, CDH4, EPHB1, ALCAM, DOCK2, GBX2, NKX2-1, ADAMTS12, ARHGDIB, FN1, FLT4, AXL, NTN1, NKX6-1, LAMA1, PTP4A3, HOXB9	0.0414058
GO:0048812~neuron projection morphogenesis	28	NDN, UCHL1, PAX6, CDH4, EPHB1, ARHGAP4, ALCAM, BDNF, BCL11A, SEMA3E, GBX2, SEMA3D, NEFH, NKX2-1, POU4F1, NRG1, FN1, PTPRZ1, EFNB1, EFNB2, NTNG1, NTN1, NKX6-1, EPB41L3, PRKCQ, EXT1, APBB1, GAP43	0.0077004
GO:0009887~organ morphogenesis	40	DLC1, FGFR1, FGD1, CTHRC1, HMX2, TNNC1, USH1G, MYO7A, TH, JAG2, PAX6, FHL2, ANKRD1, EPHB1, GPC3, PAX8, GBX2, NKX3-2, NKX2-1, POU4F1, NRG1, ACTC1, EFNB2, SIX3, TLE1, HMGA2, NTN1, ZNF22, NKX6-1, HOXD9, HOXB4, HOXD8, HOXB2, ADM, SIX1, RYR1, MFAP2, CTSH, TMEM176B, RFLNB	0.0264522
GO:0007267~cell-cell signaling	54	CTHRC1, TACR2, FGF12, SYNGR1, AMOTL1, KISS1, RSPO4, BDNF, PCDH1, GPC3, PAX8, SNPH, CHRNA7, NOS2, PLCB2, ARC, PTPRN2, EFNB1, EFNB2, SIX3, TLE1, GPER1, HMGA2, PTPRN, PCLO, INHBB, CCND1, ADRB2, ADM, NRGN, LRRK1, NEURL1, FGFR1, RAB3B, NKD2, WNT5B, HLA-DRB1, APC2, TH, OXTR, FKBP1B, AMPH, APBA2, HLA-DRB5, ADRA2C, NOVA1, FZD8, GABRA5, NKX6-1, DKK3, NPY, PLK2, HOXB9, SYTL1	0.0264522
GO:0051674~localization of cell	49	DLC1, CTHRC1, TUBB2B, PAX6, JAG2, AMOTL1, MMP1, ARHGAP4, SPINT2, SEMA3E, SEMA3D, CHRNA7, POU4F1, NOS3, NRG1, ARC, PDPN, EFNB1, EFNB2, ARHGEF5, SIX3, GPER1, SIRPA, PRKCQ, EPB41L5, SIX1, CTSH, NEURL1, FGFR1, WNT5B, WASH1, NDN, APC2, ABI3, PLVAP, MDK, EPHB1, GBX2, NKX2-1, ADAMTS12, ARHGDIB, FN1, FLT4, AXL, NTN1, NKX6-1, LAMA1, PTP4A3, HOXB9	0.0294439
GO:0048870~cell motility	49	DLC1, CTHRC1, TUBB2B, PAX6, JAG2, AMOTL1, MMP1, ARHGAP4, SPINT2, SEMA3E, SEMA3D, CHRNA7, POU4F1, NOS3, NRG1, ARC, PDPN, EFNB1, EFNB2, ARHGEF5, SIX3, GPER1, SIRPA, PRKCQ, EPB41L5, SIX1, CTSH, NEURL1, FGFR1, WNT5B, WASH1, NDN, APC2, ABI3, PLVAP, MDK, EPHB1, GBX2, NKX2-1, ADAMTS12, ARHGDIB, FN1, FLT4, AXL, NTN1, NKX6-1, LAMA1, PTP4A3, HOXB9	0.0294439
GO:2000026~regulation of multicellular organismal development	59	CTHRC1, PAX6, ANKRD1, ARHGAP4, BDNF, GPC3, PAX8, SEMA3E, SEMA3D, POU4F1, CHRNA7, NOS3, NRG1, RAMP2, EFNB2, DLL3, SIX3, GPER1, CD40, HMGA2, CCND1, ADRB2, ADM, SIX1, EPB41L5, MYRF, CTSH, NEURL1, XDH, FGFR1, OXTR, FKBP1B, ITM2C, CDH4, EPHB1, AHSG, FAM213A, BCL11A, NKX3-2, NKX2-1, ADRA2C, ADAMTS12, FN1, ARHGDIB, TESC, PTPRZ1, IL1RL2, AXL, FOXP3, NTN1, NKX6-1, NOTCH3, LAMA1, PLK2, RASSF2, H2AFY2, APBB1, TMEM176B, RFLNB	0.0335908
GO:0048598~embryonic morphogenesis	27	DLC1, FGFR1, CTHRC1, HMX2, USH1G, MYO7A, TH, PAX6, JAG2, GPC3, SPINT2, PAX8, GBX2, NKX3-2, FN1, SIX3, HMGA2, NTN1, DUSP5, HOXD9, HOXB4, HOXB2, ADM, SIX1, EPB41L5, MFAP2, EXT1	0.0349494
GO:0000902~cell morphogenesis	46	DLC1, FGD1, NDN, MYO7A, UCHL1, PAX6, CDH4, EPHB1, ALCAM, ARHGAP4, BDNF, FRMD6, IFT46, SPINT2, BCL11A, PAX8, SEMA3E, GBX2, SEMA3D, NKX2-1, NEFH, POU4F1, NRG1, FN1, CDC42EP5, COCH, PTPRZ1, PDPN, EFNB1, EFNB2, NTNG1, AXL, SIPA1L3, ACTN1, HMGA2, NTN1, NKX6-1, EPB41L3, PRKCQ, EPB41L5, KLF2, EXT1, APBB1, GAP43, RFLNB, NEURL1	0.0371107
GO:0045595~regulation of cell differentiation	54	CTHRC1, PAX6, DLK2, ANKRD1, ARHGAP4, BDNF, PAX8, SEMA3E, SEMA3D, POU4F1, NRG1, RAMP2, PRAME, EFNB2, SIX3, DLL3, GPER1, HMGA2, CCND1, ADM, SIX1, EPB41L5, NEURL1, XDH, FGFR1, WNT5B, APC2, FKBP1B, ITM2C, CDH4, FAM213A, BCL11A, NKX3-2, NKX2-1, ADRA2C, ADAMTS12, OLFM2, FN1, LPL, TESC, PTPRZ1, IL1RL2, AXL, FOXP3, NTN1, NKX6-1, NOTCH3, HOXB4, PLK2, RASSF2, H2AFY2, APBB1, RFLNB, TMEM176B	0.0381061
GO:0009790~embryo development	38	DLC1, FGFR1, CTHRC1, HMX2, USH1G, MYO7A, TH, JAG2, PAX6, GPC3, SPINT2, PAX8, GBX2, NKX3-2, APBA2, NOS3, ETNK2, NRG1, FOSL1, FN1, EFNB1, DLL3, SIX3, HMGA2, NTN1, HOXD9, DUSP5, HOXB4, SALL2, HOXD8, HOXB2, ADM, SIX1, EPB41L5, MFAP2, HOXB9, KLF2, EXT1	0.0373992
GO:0032989~cellular component morphogenesis	48	DLC1, FGD1, NDN, MYO7A, UCHL1, PAX6, ANKRD1, CDH4, EPHB1, ALCAM, ARHGAP4, BDNF, FRMD6, IFT46, SPINT2, BCL11A, PAX8, SEMA3E, GBX2, SEMA3D, NKX2-1, NEFH, POU4F1, NRG1, FN1, CDC42EP5, COCH, ACTC1, PTPRZ1, PDPN, EFNB1, EFNB2, NTNG1, AXL, SIPA1L3, ACTN1, HMGA2, NTN1, NKX6-1, EPB41L3, PRKCQ, EPB41L5, KLF2, EXT1, APBB1, RFLNB, GAP43, NEURL1	0.0375119
GO:0009914~hormone transport	18	FGFR1, HLA-DRB1, TACR2, PTPRN2, GPER1, HMGA2, PTPRN, PCLO, FKBP1B, NKX6-1, INHBB, KISS1, SLC16A2, ADM, PAX8, HLA-DRB5, ADRA2C, NOS2	0.0415878
GO:0007411~axon guidance	15	EFNB1, EFNB2, PAX6, CDH4, NTN1, EPHB1, ALCAM, PRKCQ, BDNF, SEMA3E, GBX2, SEMA3D, NKX2-1, EXT1, GAP43	0.0410559
GO:0097485~neuron projection guidance	15	EFNB1, EFNB2, PAX6, CDH4, NTN1, EPHB1, ALCAM, PRKCQ, BDNF, SEMA3E, GBX2, SEMA3D, NKX2-1, EXT1, GAP43	0.0414269
GO:0061387~regulation of extent of cell growth	10	ARHGAP4, BDNF, BCL11A, SEMA3E, SEMA3D, NRG1, NTN1, CDH4, NKX6-1, FN1	0.0486247

GO, gene ontology; FDR, false discovery rate.

Figure 3.

GO analysis of the enrichment of differentially expressed genes between MCF-7/MDR and MCF-7 cells in cellular component terms. GO, gene ontology; MDR, multidrug resistance.

Table III.

Gene Ontology cellular component analysis of differentially expressed genes.

Term	Count	Genes	FDR
GO:0098590~plasma membrane region	42	DLC1, NKD2, APC2, CLDN6, MYO7A, OXTR, PLVAP, AMOTL1, AMN, GPR143, KISS1, C1QTNF5, SH2D5, SNPH, CHRNA7, NOS3, IQSEC3, SLC39A4, NRG1, FOSL1, OLFM2, GABRP, FN1, KCNMA1, GABRG3, TESC, ARC, KIF17, PDPN, GABRA5, SIPA1L3, ABCB1, AJAP1, GPER1, ABCB4, EPB41L3, ADRB2, EPB41L5, CLDN1, APBB1, GAP43, NEURL1	0.0117429
GO:0005578~proteinaceous extracellular matrix	22	COCH, CTHRC1, WNT5B, PTPRZ1, EFEMP2, CCDC80, EMILIN2, NTN1, MMP1, LAMA1, ADAMTS6, COL9A2, CILP2, GPC3, MMP23B, FBLN5, PI3, MFAP2, ADAMTS12, TFPI2, MFAP5, FN1	0.0147915
GO:0045202~synapse	35	RAB3B, MYO7A, TH, SYNGR1, AMPH, PLCB4, SH2D5, SNPH, NEFH, APBA2, CHRNA7, SV2A, IQSEC3, NRG1, FOSL1, OLFM2, SAMD4A, GABRP, NMNAT2, GABRG3, ARC, PTPRN2, EFNB1, GABRA5, ATP1A3, GPER1, PTPRN, PCLO, EPB41L3, NRGN, APBB1, SYTL1, ADD2, GAP43, NEURL1	0.0119474
GO:0044456~synapse part	30	RAB3B, TH, SYNGR1, AMPH, PLCB4, SH2D5, SNPH, APBA2, NEFH, CHRNA7, SV2A, IQSEC3, FOSL1, OLFM2, GABRP, ARC, GABRG3, PTPRN2, GABRA5, ATP1A3, GPER1, PTPRN, PCLO, EPB41L3, NRGN, APBB1, SYTL1, GAP43, ADD2, NEURL1	0.0119478
GO:0031012~extracellular matrix	27	CTHRC1, AEBP1, WNT5B, MMP1, AHSG, COL9A2, CILP2, GPC3, SBSN, ADAMTS12, TFPI2, FN1, COCH, LPL, PTPRZ1, EFEMP2, CCDC80, EMILIN2, NTN1, SOD3, LAMA1, ADAMTS6, MMP23B, FBLN5, PI3, MFAP2, MFAP5	0.0139638
GO:0016324~apical plasma membrane	18	KCNMA1, MYO7A, OXTR, SIPA1L3, ABCB1, AJAP1, GPR143, AMN, AMOTL1, ABCB4, KISS1, C1QTNF5, ADRB2, CLDN1, CHRNA7, NRG1, SLC39A4, FN1	0.0254454

KEGG pathway enrichment analysis

The enrichment of dysfunctional signaling pathways was screened by the KEGG pathway analysis, and the results demonstrated that the pathways that differentially expressed genes were enriched in included ‘pathways in cancer’ (16 genes), ‘calcium signaling pathway’ (12 genes), ‘cell adhesion molecules (CAMs)’ (11 genes), ‘tight junction’ (9 genes), ‘amoebiasis’ (9 genes), ‘Wnt signaling pathway’ (8 genes), ‘retrograde endocannabinoid signaling’ (7 genes), ‘type I diabetes mellitus’ (5 genes), ‘intestinal immune network for IgA production’ (5 genes), ‘asthma’ (4 genes), ‘viral myocarditis’ (5 genes), ‘allograft rejection’ (4 genes; Fig. 4 and Table IV).

Figure 4.

KEGG analysis of the enrichment of differentially expressed genes between MCF-7/MDR and MCF-7 cells in specific pathways. KEGG, Kyoto Encyclopedia of Genes and Genomes; MDR, multidrug resistance.

Table IV.

Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis of differentially expressed genes.

Term	Count	P-value	Genes	FDR
hsa04514: Cell adhesion molecules (CAMs)	11	0.0010925	ALCAM, HLA-DQB1, HLA-DRB1, CLDN6, ICAM2, CLDN1, NTNG1, HLA-DRB5, CLDN11, CD40, CDH4	0.0131104
hsa04020: Calcium signaling pathway	12	0.0010925	PTGER1, GNA15, ADRB2, PLCB4, TACR2, TNNC1, RYR1, OXTR, CHRNA7, NOS3, NOS2, PLCB2	0.0113229
hsa05146: Amoebiasis	9	0.0010925	IL1R2, LAMA1, GNA15, PLCB4, SERPINB2, ACTN1, NOS2, PLCB2, FN1	0.0088997
hsa04530: Tight junction	9	0.0104316	EPB41L3, PRKCQ, RAB3B, CLDN6, CLDN1, ACTN1, CLDN11, AMOTL1, PPP2R2C	0.0312949
hsa04940: Type I diabetes mellitus	5	0.0131689	HLA-DQB1, HLA-DRB1, PTPRN2, HLA-DRB5, PTPRN	0.0316055
hsa04672: Intestinal immune network for IgA production	5	0.0192920	HLA-DQB1, HLA-DRB1, CCR10, HLA-DRB5, CD40	0.0385841
hsa04723: Retrograde endocannabinoid signaling	7	0.0192920	GABRG3, GNGT2, PLCB4, GABRA5, GNG4, PLCB2, GABRP	0.0402837
hsa05200: Pathways in cancer	16	0.0192920	LAMA1, FZD8, FGFR1, PTGER1, CCND1, GNGT2, PLCB4, WNT5B, APC2, PAX8, FGF12, NOS2, GNG4, PLCB2, MMP1, FN1	0.0384532
hsa05310: Asthma	4	0.0274460	HLA-DQB1, HLA-DRB1, HLA-DRB5, CD40	0.0365947
hsa04310: Wnt signaling pathway	8	0.0321127	FZD8, CCND1, WNT5B, NKD2, PLCB4, APC2, PLCB2, FOSL1	0.0385353
hsa05416: Viral myocarditis	5	0.0321127	HLA-DQB1, CCND1, HLA-DRB1, HLA-DRB5, CD40	0.0394291
hsa05330: Allograft rejection	4	0.0321127	HLA-DQB1, HLA-DRB1, HLA-DRB5, CD40	0.0470282

Protein-protein interaction network

STRING, a database that comprises established and predicted protein interactions, was employed in the present study to predict protein interactions among the differentially expressed genes. A total of 9 hub nodes were screened based on their association with other proteins. These hub genes included cyclin D1 (CCND1), nitric oxide synthase 3 (NOS3), NOTCH3, brain-derived neurotrophic factor (BDNF), paired box 6 (PAX6), neuropeptide Y (NPY), phospholipase C β (PLCB) 4, PLCB2 and actin α cardiac muscle 1 (ACTC1) (Table V). Among these genes, CCND1 exhibited the highest node degree of 22. A protein-protein interaction network was constructed as presented in Fig. 5. RT-qPCR was subsequently performed to verify the expression of these 9 genes in MCF-7 and MCF-7/MDR cells; the expression of all genes was higher in MCF-7/MDR cells compared with MCF-7 cells (Fig. 6), which was consistent with the RNA-seq data.

Table V.

Key hub genes in the protein-protein interaction network.

Gene	Degree
CCND1	22
NOS3	19
NOTCH3	17
BDNF	15
PAX6	14
NPY	13
PLCB4	11
PLCB2	11
ACTC1	10

CCND1, cyclin D1; NOS3, nitric oxide synthase 3; BDNF, brain-derived neurotrophic factor; PAX6, paired box 6; NPY, neuropeptide Y; PLCB, phospholipase C β; ACTC1, actin α cardiac muscle 1.

Figure 5.

Protein-protein interaction network of the differentially expressed genes between MCF-7/MDR and MCF-7 cells. Circles represent proteins, lines represent strong association between proteins. MDR, multidrug resistance.

Figure 6.

Reverse transcription-quantitative polymerase chain reaction in MCF-7/MDR and MCF-7 cells verified the mRNA expression of the 9 hub genes identified in the protein-protein interaction network. *P<0.05 and **P<0.01 vs. MCF-7 cells. MDR, multidrug resistance; BDNF, brain-derived neurotrophic factor; PLCB, phospholipase C β; ACTC1, actin α cardiac muscle 1; NOS3, nitric oxide synthase 3; CCND1, cyclin D1; NPY, neuropeptide Y; PAX6, paired box 6.

Inhibition of NOS3 and BDNF decreases chemoresistance properties of MCF-7/MDR cells

To verify the accuracy and practicality of RNA-seq in investigating the mechanism involved in breast cancer chemoresistance, specific inhibitors of NOS3 and BDNF were employed. As the expression of NOS3 was higher at mRNA and protein levels in MCF-7/MDR cells compared with MCF-7 cells (Figs. 6 and 7A), different concentrations of LNIO, which is a specific inhibitor of NOS3, was added to determine the optimal inhibitory concentration for subsequent drug sensitivity assays. The results of western blotting demonstrated that LNIO decreased the expression of NOS3 even at a low concentration (5 µM) and 10 µM was selected as the optimal concentration (Fig. 7B). LNIO (10 µM) was added to MCF-7/MDR cells prior to treatment with anticancer drugs and a CCK-8 assay was performed to measure cell viability. IC50 values were subsequently calculated in response to different anticancer drugs. The results demonstrated that 10 µM LNIO decreased the IC50 of ADM, DDP and 5-FU in MCF-7/MDR cells, as demonstrated in Table VI. Furthermore, 25 µM ANA-12, a specific inhibitor of BDNF/neurotrophic receptor tyrosine kinase, type 2 signaling, decreased the mRNA expression of BDNF and TrkB (Fig. 8), and decreased the IC50 of ADM and 5-FU in MCF-7/MDR cells (Table VII).

Figure 7.

LNIO inhibited the protein expression of NOS3 in MCF-7/MDR cells. (A) NOS3 expression in MCF-7/MDR and MCF-7 cells was detected by western blot analysis. (B) NOS3 protein expression under different concentrations of LNIO was detected by western blot analysis in MCF-7/MDR cells. *P<0.05 and **P<0.01 as indicated. LNIO, L-N5-(1-Iminoethyl)-ornithine hydrochloride; NOS3, nitric oxide synthase 3; MDR, multidrug resistance.

Table VI.

Inhibition of nitric oxide synthase 3 decreased the IC50 of MCF-7/MDR cells in response to various anticancer drugs.

	IC50 value
Anticancer drug	MCF-7/MDR	MCF-7/MDR+LNIO	P-value
ADM, µM	44.24±10.37	26.86±2.05	0.0160001
DDP, µg/ml	4.63±1.45	4.01±1.62	0.0276806
5-FU, µg/ml	5.335±1.23	4.14±1.81	0.0443251

MDR, multidrug resistance; LNIO, L-N5-(1-Iminoethyl)-ornithine hydrochloride; ADM, adriamycin; DDP, cisplatin; 5-FU, 5-fluorouracil.

Figure 8.

Reverse transcription-quantitative polymerase chain reaction results demonstrated that various concentrations of ANA-12 inhibited the mRNA expression of BDNF and TrkB in MCF-7/MDR cells. *P<0.05 and **P<0.01, as indicated. BDNF, brain-derived neurotrophic factor; TrkB, neurotrophic receptor tyrosine kinase, type 2; MDR, multidrug resistance.

Table VII.

Inhibition of brain-derived neurotrophic factor decreased the IC50 of MCF-7/MDR cells in response to various anticancer drugs.

	IC50 value
Anticancer drug	MCF-7/MDR	MCF-7/MDR+ANA-12	P-value
ADM, µM	43.98±1.52	16.80±2.07	0.0000521
DDP, µg/ml	6.49±0.17	6.78±0.19	0.1182792
5-FU, µg/ml	7.06±0.41	4.74±0.67	0.0007725

MDR, multidrug resistance; ADM, adriamycin; DDP, cisplatin; 5-FU, 5-fluorouracil.

Discussion

MDR is a major obstacle in the treatment of breast cancer. While the mechanism of MDR is complex, it is particularly difficult to identify the key genes and signaling pathways involved in regulating MDR. Large-scale investigation of gene expression levels is required in order to screen and identify the important factors involved in this phenomenon. As RNA-seq is considered to be a promising technology for measuring mRNA expression and may be used to identify differentially expressed genes in a manner that is superior to certain existing array-based approaches (39), the present study applied RNA-seq technology to investigate the mechanisms involved in MDR in human breast cancer cells. While a total of 2085 differentially expressed genes were screened out, GO analysis was performed to analyze the differentially expressed genes in terms of their enrichment in biological process and cellular component terms. GO analysis revealed that MDR in breast cancer MCF-7 cells was associated with ‘neuron differentiation’ and ‘development of nervous system’ biological process terms, which was consistent with the results of a previous study that demonstrated the involvement of the central nervous system in patients with advanced breast cancer receiving effective drugs (40), though the mechanism involved remains unclear. In addition, MDR was also associated with ‘cell motility’, ‘locomotion’ and ‘movement of cell or subcellular component’ biological process terms. In various cancer types, such as lung, atypical teratoid and ovarian cancers, cell motility has been reported to be closely correlated with drug resistance (41,42). In the GO analysis of cellular component terms, ‘plasma membrane region’ and ‘extracellular matrix’ were among the terms that differentially expressed genes were demonstrated to be enriched in; previously, cell surface markers CD9 and CD44 were reported to be associated with chemoresistance (43,44) and, in small cell lung cancer, it was also reported that extracellular matrix (ECM) is an important regulator of drug resistance (45).

KEGG analysis was performed in the present study to determine the enrichment of signaling pathways associated with MDR in MCF-7 cells. The results revealed that 12 different pathways may be involved. A number of these identified pathways were consistent with previous reports, for example, CAM signaling was previously reported to mediate chemoresistance in renal cell carcinoma and pancreatic adenocarcinoma (46), the calcium signaling pathway has been reported to affect MDR in cancer chemotherapy (47) and the Wnt signaling has been implicated in the regulation of chemoresistance via various mechanisms, including long non-coding RNA (48) and the canonical Wnt/β-catenin signaling pathway (49). Furthermore, activation of Wnt signaling was important in the response to certain anticancer drugs, particularly 5-fluorouracil and platinum (50,51).

The results of GO biological process analysis of differentially expressed genes also demonstrated that cell mobility may be a key regulator of MDR, and mobility is the basis of tumor cell metastasis. Matrix metalloproteinases (MMPs), a family of enzymes that degrade ECM components, also mediated MDR in various cancers, including MMP-2, MMP-7 and MMP-9 enzymes (52–54). Claudins (CLDNs), which are an important component of tight junctions, have also been reported to be factors involved in MDR (55–58), as 3 of 27 CLDN family members were observed to be differentially expressed in MCF-7/MDR cells, including CLDN-1, CLDN-6 and CLDN-11. In our previous study, the results demonstrated that CLDN-6 had an important role in breast cancer (59), verifying that tight junctions, which was a term screened in the KEGG enrichment analysis of the present study, may be involved in the MDR of breast cancer.

Through protein-protein interaction network construction, a series of hub genes were demonstrated to be associated with MDR. In total, 9 hub genes were identified in the present study, including CCND1, NOS3, NOTCH3, BDNF, PAX6, NPY, PLCB4, PLCB2 and ACTC1. These genes were previously reported to be associated with drug resistance in various cancers. For example, CCND1, which has an established role in the cell cycle, was also reported to be associated with chemoresistance in various cancers, including testicular germ cell tumors and cell lymphoma (60). NOS3, also termed endothelial NOS, contributed to DDP resistance through various signaling pathways in human ovarian cancer cells (61). NOTCH3 has an important role in calcium signaling and its overexpression enhanced chemoresistance in response to carboplatin, gemcitabine and paclitaxel via the Notch signaling pathway (62–64). Furthermore, a recent study demonstrated that NOTCH3 may be regarded as a predictor of distant relapse-free survival in estrogen receptor-positive breast cancers and its expression was associated with clinical chemoresistance in patients with breast cancer (65). BDNF was the gene with the largest difference in expression between MCF-7 and MCF-7/MDR cells when verified by RT-qPCR, and BDNF/TrkB signaling has been reported to be responsible for chemoresistance in neuroblastoma cells and head and neck squamous cell carcinoma (66,67). Additionally, PAX6 was reported to mediate the temozolomide sensitivity of glioblastoma stem cells (68), NPY signaling was activated in chemotherapy-treated neuroblastoma tumors, which indicated that NPY may be associated with chemoresistance in neuroblastoma (69), and Taxotere® (docetaxel) treatment influenced the expression of ACTC1 in breast cancer (70). Numerous genes involved in biological process, cellular component and signal transduction pathways may be involved in regulating MDR. As 9 genes were screened in the protein-protein interaction network, taking NOS3 as an example, NOS3 was demonstrated to be associated with ‘negative regulation of multicellular organismal process’, ‘movement of cell or subcellular component’, ‘locomotion’, ‘localization of cell’, ‘cell motility’, ‘regulation of multicellular organismal development’ and ‘embryo development’ terms in the GO biological process analysis. NOS3 has also been demonstrated to be involved in embryo development (71). In addition, other factors such as EMT, an important process, have also been associated with embryo development (72); it has been verified that EMT promoted chemoresistance in various cancers (73). In the GO cellular component analysis of the present study, NOS3 belonged to the ‘plasma membrane region’ term. Water-soluble drugs, including DDP, nucleoside analogues and antifolates, are only able to traverse the plasma membrane via membrane transporters or through hydrophilic channels within the membrane (74). Therefore, the plasma membrane region has an essential role in chemoresistance. In the KEGG pathway enrichment analysis of the present study, NOS3 was associated with the ‘calcium signaling pathway’ term, and it has been reported that the regulation of RhoA/calcium signaling may regulate chemoresistance in head and neck squamous cell carcinoma (75). Furthermore, the current study confirmed that LNIO, a specific inhibitor of NOS3, significantly decreased the resistance of MCF-7/MDR cells to various anticancer drugs, including ADM, 5-FU and DDP, which supported a role for NOS3 in chemoresistance. Additionally, a specific inhibitor of BDNF/TrkB signaling, ANA-12, was also employed to verify the role of BDNF in breast cancer chemoresistance; ANA-12 significantly decreased resistance to ADM and 5-FU, but not DDP, in MCF-7/MDR cells. Thus, we hypothesized that those hub genes of NOS3 and BDNF, may contribute to drug resistance through regulating biological processes and signaling pathways in breast cancer cells.

In conclusion, the results of the current study provide a comprehensive analysis of differentially expressed genes in MCF-7/MDR cells and the associated signaling pathways, which may be involved in the development of chemoresistance in breast cancer. As the hub genes screened, including NOS3 and BDNF, may provide insights into the mechanism of drug resistance, genome-wide mRNA profiling may be useful in identifying novel targets for the treatment of breast cancer that exhibits MDR.