Identification of potential therapeutic target of naringenin in breast cancer stem cells inhibition by bioinformatics and in vitro studies.

Cancer therapy is a strategic measure in inhibiting breast cancer stem cell (BCSC) pathways. Naringenin, a citrus flavonoid, was found to increase breast cancer cells' sensitivity to chemotherapeutic agents. Bioinformatics study and 3D tumorsphere in vitro modeling in breast cancer (mammosphere) were used in this study, which aims to explore the potential therapeutic targets of naringenin (PTTNs) in inhibiting BCSCs. Bioinformatic analyses identified direct target proteins (DTPs), indirect target proteins (ITPs), naringenin-mediated proteins (NMPs), BCSC regulatory genes, and PTTNs. The PTTNs were further analyzed for gene ontology, Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment, protein-protein interaction (PPI) networks, and hub protein selection. Mammospheres were cultured in serum-free media. The effects of naringenin were measured by MTT-based cytotoxicity, mammosphere forming potential (MFP), colony formation, scratch wound-healing assay, and flow cytometry-based cell cycle analyses and apoptosis assays. Gene expression analysis was performed using real-time quantitative polymerase chain reaction (q-RT PCR). Bioinformatics analysis revealed p53 and estrogen receptor alpha (ERα) as PTTNs, and KEGG pathway enrichment analysis revealed that TGF-ß and Wnt/ß-catenin pathways are regulated by PTTNs. Naringenin demonstrated cytotoxicity and inhibited mammosphere and colony formation, migration, and epithelial to mesenchymal transition in the mammosphere. The mRNA of tumor suppressors P53 and ERα were downregulated in the mammosphere, but were significantly upregulated upon naringenin treatment. By modulating the P53 and ERα mRNA, naringenin has the potential of inhibiting BCSCs. Further studies on the molecular mechanism and formulation of naringenin in BCSCs would be beneficial for its development as a BCSC-targeting drug.

Introduction

Cancer is a deadly disease triggered by impaired cell proliferation and apoptosis (Hanahan and Weinberg, 2011). According to the World Health Organization, breast cancer is the second leading cause of death in women worldwide, and its incidence is predicted to increase (WHO, 2015). Along with surgery and radiotherapy, chemotherapy is widely used as a neoadjuvant in breast cancer therapy (Alvarez et al., 2010). The emergence of acquired resistance and the progression of cancer cells and metastasis resulted in treatment failure, where only 50% of the patients have shown good response to chemotherapy (Liu, 2009).

The existence of resistance to the use of chemotherapy agents is partly due to the presence of a population of cancer stem cells (CSCs) that can differentiate, causing relapse and metastasis and death in patients (Peitzsch et al., 2017). Breast cancer stem cell (BCSC) emergence is characterized by the cell’s ability to inactivate drugs, alter target cells, and suppress the accumulation of drugs entering the cells through increased efflux pumping (Badve and Nakshatri, 2012, Moitra, 2015). Therefore, breast cancer therapy, with inhibition of BCSC pathways, is a strategic and targeted approach.

Naringenin is a citrus flavonoid compound potentially developed as a BCSC-targeted drug. Previously, naringenin (Fig. 1A), a citrus flavonoid, has been shown to have cytotoxic and antiproliferative effects on several cancer cell types, one of which is the MCF-7 cells (Kanno et al., 2005). Naringenin inhibits proliferation, migration, and invasion, triggers apoptosis of SGC7901 gastric cancer cells through downregulation (Bao et al., 2016). Naringenin also induces apoptosis, increases intracellular reactive oxygen species (ROS), and upregulates ASK-1, p53, p38, and JNK in SNU-213 pancreatic cancer cells (Park et al., 2017).

Fig. 1

(A). The structure of naringenin. (B). Interactions of naringenin and its DTPs. (C). A Venn diagram of PTT related to Naringenin and BCSCs. (C). PPI network of naringenin protein targets, analyzed using STRING. (D). Top 10 hub proteins based on degree score, analyzed by CytoScape. (E). GO enrichment of PTTNs, analyzed by WebGestalt. (F). Overview of changes in PRKCA, EGFR, ERBB4, AREG, ESR1, and STAT1 in genomics dataset from 16 studies of breast cancer. (G). Summary alterations of ESR1, AKT1, EGFR, WNT5A, MAPK8, JUN, IL6, and TP53 across breast cancer samples (based on a study by Levebvre et al., 2016). (H). Gene network and (I). Drug-gene network connected to ESR1, AKT1, EGFR, WNT5A, MAPK8, JUN, IL6, and TP53 in breast cancer samples (based on a study by Lefebvre et al., 2016).

ROS play a critical role in cancer progression, in which it induces angiogenesis, metastasis, and chemoresistance in low concentrations and promotes cell death in high concentrations (Aggarwal et al., 2019). On the one hand, naringenin exhibits anti-inflammatory activity by decreasing ROS and antioxidant content in CD45+ hematopoietic cells (Bussmann et al., 2019) and RAW 264.7 cell line (Teng et al., 2018); and on the other hand, it inhibits diethylnitrosamine/acetylaminofluorene-induced hepatocarcinogenesis by inducing apoptosis and modulating oxidative stress (Ahmed et al., 2019).

To overcome chemotherapy resistance, studies on naringenin as a combination agent with daunomycin, mitoxantrone, doxorubicin, and vincristine in MCF-7/ADR breast cancer cells, NCI-H460 lung cancer cells, and A549 lung cancer cells were carried out (Mitsunaga et al., 2000, Chung et al., 2005). The combination of naringenin with doxorubicin shows an increase in the effectiveness of doxorubicin treatment in monolayer MCF-7 cells (Fitriasari et al., 2010). Moreover, combining naringenin and TNF-related apoptosis-inducing ligand (TRAIL) is an appropriate and safe strategy in treating apoptosis-resistant non-small cell lung cancer (Jin et al., 2011). The combination of naringenin and ABT-737 (Bcl-2 and Bcl-xL dual inhibitor) increases the effect of ABT-737 on gastric cancer cells (Zhang et al., 2016). Although naringenin promotes apoptosis in cancer cells, it can prevent apoptosis in normal cells stimulated by radiation by means of p53-, Bax-, and Bcl-2-mediated expansion modulation (Kumar and Tiku, 2016). However, research on the effects of naringenin as well as the molecular mechanisms of its action in BCSCs has not been conducted.

This study aimed to examine the potential of naringenin in BCSC inhibition using a 3D tumorsphere from MCF-7 breast cancer cells. Moreover, an integrated bioinformatics analysis was also used in predicting the potential therapeutic target of naringenin (PTTN) in BCSCs, therefore highlighting the association between naringenin and BCSC inhibition. The results of this study are available as reference for furture studies in the development of BCSC-targeted drugs capable of overcoming chemoresistance.

Methods

Bioinformatics analysis

In cells, a compound can target or interact with proteins. Proteins that interact or become targets of compounds are called direct target proteins (DTPs). Naringenin DTPs were searched from STITCH (http://stitch.embl.de) (Kuhn et al., 2007) and DrugBank (https://www.drugbank.ca) (Wishart et al., 2018), and the results were used for the following analysis. Proteins that interact with DTPs are defined as indirect target proteins (ITPs). The ITPs of each DTP were obtained from the STRING database (https://string-db.org) (Szklarczyk et al., 2015), with a minimum interaction score of 0.4 and a maximum number of interactors of 20. After removing repetitive proteins, the ITPs of all DTPs were generated. All DTPs and ITPs are naringenin-mediated proteins (NMPs). Breast cancer stem cell regulatory genes were retrieved from the PubMed gene database with the keywords “breast cancer stem cells.“ The PTTN against BCSCs are proteins that potentially target BCSCs and poses as NMPs. Analysis of protein–protein interaction (PPI) network of the PTTN was constructed with STRING-DB v11.0 (Szklarczyk et al., 2015), with confidence scores of > 0.7 and visualized by Cytoscape software (Shannon et al., 2003). As analyzed by CytoHubba plugin, proteins with a degree score of > 10 were selected as hub proteins (Chin et al., 2014). Analysis of Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment of the PTTN was conducted using the Database for Annotation, Visualization and Integrated Discovery (DAVID) v6.7 (Huang da et al., 2009), with p < 0.05 selected as the cutoff value.). The genetic alterations of selected genes were analyzed using cBioPortal (http://www.cbioportal.org) (Cerami et al., 2012, Gao et al., 2013). In the present study, protein genes such as EGFR, AKT1, WNT5A, MAPK8, JUN, TP53, IL6, and ESR1 were screened for genetic alterations in all breast cancer studies available in the cBioportal database, and the study with the highest number of genetic alterations was chosen for further connectivity analysis.

Cell culture and mammosphere generation

MCF-7 cells (2D cell culture) were cultured in Dulbecco's Modified Eagle's medium (DMEM) (Gibco) with a 10% fetal bovine serum (FBS) (Sigma-Aldrich), 1.5% penicillin–streptomycin (Sigma-Aldrich), and 0.5% fungizone (Sigma-Aldrich) supplementation. To generate the mammosphere (3D cell culture), MCF-7 cells were cultured in 5% polyHEMA (Sigma-Aldrich)-coated plate using DMEM (Gibco) with the addition of 10 ng/mL of epidermal growth factor (EGF) (Sinobiological), 10 ng/mL of B27 supplement (Gibco), 5 µg/mL of insulin (Gibco), 1.5% of penicillin–streptomycin (Sigma-Aldrich), and 0.5% of fungizone (Thermo Fischer Scientific), as previousy described (Grimshaw et al., 2008, Pickl and Ries, 2009, Oak et al., 2012). Cells were incubated at 37 °C with 5% CO2 and were routinely tested and confirmed as mycoplasma free.

Drugs

Naringenin or NAR (N5893, purity ≥ 95%), metformin (MET, D150959, purity 97%), and doxorubicin (DXR, 44583, purity 98.0%–102.0%) were purchased from Sigma-Aldrich. In previous studies, MET and DXR were used as a positive control for the mammosphere forming experiment and for the clonogenic assay, respectively (Lou et al., 2006, Hirsch et al., 2009).

CD44+/CD24− phenotype enrichment assay

The CD44+/CD24− population phenotype enrichment assay was performed according to the previous study with modifications (Seo et al., 2015). The mammosphere was separated into a single cell by trypsinization. The cell suspension was incubated with phycoerythrin (PE)-conjugated anti-human CD44 antibody and fluorescein isothiocyanate (FITC)-conjugated anti-human CD24 mouse antibody (Molecular Probes) on ice and dark conditions for 20 min. Cells were washed after incubation, and were measured by flow cytometer (BD FACS Calibur) at a wavelength of 561 and 488 nm, respectively.

Cytotoxic assay

Cytotoxicity testing was performed based on a previous study (Mosmann, 1983). Cells were cultured on 96-well plates for the 2D model and polyHEMA-coated 96-well plate for the 3D model. These cells were treated with samples in a series concentration for 72 h, followed by MTT reagent addition, and were incubated for 2–4 h. In each well, SDS was added after incubation to stop the reaction. Using a microplate reader set at a 570-nm wavelength, the absorbance was measured, and the cell viability percentage was calculated.

Mammosphere forming potential

With reference to Oak et al., mammosphere forming potential (MFP) was calculated (Oak et al., 2012). Briefly, cells were treated with samples for 72 h then cultured on polyHEMA-coated 96-well plates for 96 h. The number of mammospheres formed was calculated manually.

Clonogenic survival assay

Cells were cultured on 6-well plates and treated with samples for 72 h. The media was replaced with fresh media. Cells were incubated for the following 14 days, and were fixed with methanol and stained with gentian violet solution after incubation. The Colony Area plugin of ImageJ was used to analyze the remaining colony formation, as previously described (Guzman et al., 2014).

Cell cycle modulation assay

Cells were cultured on 6-well plates (as both 2D and 3D models) and were treated with samples for 72 h. These cells were harvested and separated into cell suspensions, which were stained using propidium iodide (PI) reagent, transferred into a cytometric tube, and analyzed using a flow cytometer (BD FACS Calibur, BD Bioscience).

Apoptosis assay

As both 2D and 3D models, cells were cultured on 6-well plates and were treated with samples for 96 h. These cells were harvested, separated into a cell suspension, and stained using an FTIC Annexin V (Roche) staining kit according to the manufacturer’s instructions and were analyzed using a flow cytometer (BD FACS Calibur, BD Bioscience).

Scratch wound-healing assay

Mammosphere-derived MCF-7 cells (70,000 cells/well) were seeded in a 24-well plate and incubated for 24 h and were starved with the serum-free medium for another 24 h. The cells were scratched after starvation using a sterile pipette tip and were treated with naringenin. Images of the cells were captured at 0, 18, 24, 42, and 48 h after treatment by a digital camera (Nikon, Japan). ImageJ software was used to analyze the results, presented as percentage closure.

Quantitative real-time polymerase chain reaction (q-RT PCR)

Total mRNA was isolated using GeneJet RNA Purification Kit (Thermoscientific). A RevertAid First Strand cDNA Synthesis Kit (Thermoscientific) was utilized to transcribe RNA to cDNA. To quantify real-time polymerase chain reaction (RT PCR), SsoFast Evagreen Supermix (Bio-Rad) was used. Stem cell markers SOX2, CD133, and ALDH1 were examined, whereas GAPDH was used as a housekeeping gene. The results were analyzed using a comparative threshold cycle (ΔΔCT). GAPDH was used as an internal control. RT PCR was performed using the selected primers (Table S1).

Statistical analysis

From at least three independent experiments, data were presented as mean ± SD and were analyzed using GraphPad Prism 5.0 software. To analyze the CD44+/CD24− population enrichment assay and cytotoxic effect of NAR, two-way ANOVA with post-hoc Bonferroni’s multiple comparisons test was used, whereas one-way ANOVA with post-hoc Bonferroni’s multiple comparisons test was used to analyze the effect of NAR on mammosphere and colony formation. Student’s t-test was used to analyze the anti-migratory effect of NAR and quantitative RT PCR (q-RT PCR) results.

Results

To obtain DTPs, we searched the DTPs of naringenin in STITCH and DrugBank databases and obtained 17 DTPs such as ESR1, AKR1C1, and CYP1B1 (Table 1). Fig. 1B shows the analysis of the interactions between naringenin and DTPs. A total of 249 NMPs containing 17 DTPs and 232 ITPs were obtained (Table S2), whereas 1041 BCSC regulatory genes from PubMed were retrieved (Table S3). A Venn diagram showed 65 PTTNs in BCSCs (Fig. 1C; Table S4).

Table 1

Direct target protein of naringenin, from DrugBank and STITCH.

No	Protein symbol	Protein name	Database
1	ESR1	Estrogen receptor alpha	DrugBank
2	AKR1C1	Aldo-keto reductase family 1 member C1	DrugBank
3	CYP1B1	Cytochrome P450 1B1	DrugBank
4	KANSL3	KAT8 regulatory NSL complex subunit 3	DrugBank
5	SHBG	Sex hormone-binding globulin	DrugBank
6	CYP19A1	Cytochrome P450 19A1	DrugBank
7	ESR2	Estrogen receptor beta	DrugBank
8	RAPGEF1	Rap guanine nucleotide exchange factor 1	STITCH
9	ABCB1	ATP binding cassette subfamily B1	STITCH
10	CYP1A2	Cytochrome P450 1A2	STITCH
11	CYP1B1	Cytochrome P450 1B1	STITCH
12	LDLR	Low-density lipoprotein receptor	STITCH
13	APOB	Apolipoprotein B-100	STITCH
14	PPARA	Peroxisome proliferator-activated receptor alpha	STITCH
15	CCL2	C-C motif chemokine 2	STITCH
16	HMOX1	Heme oxygenase 1	STITCH
17	BDNF	Brain-derived neurotrophic factor	STITCH

A PPI network of PTTNs was generated using a STRING database consisting 65 nodes, 553 edges, average node degree of 17, high confidence interaction of 0.7, average local clustering coefficient of 0.576, and PPI enrichment value < 1.10e-16 (Fig. 1D). Furthermore, hub proteins were selected from the PPI network based on a degree score (Fig. 1E). ESR1 was the only DTP among the hub proteins with a degree score of 26, indicating strong correlation with the biological effect of naringenin (Table 2).

Table 2

The top 20 hub protein by degree score.

No	Protein symbol	Degree score
1	STAT3	41
2	IL6	40
3	AKT1	36
4	JUN	34
5	EGFR	34
6	INS	34
7	MAPK8	32
8	TNF	31
9	TP53	28
10	IGF1	27
11	MYC	26
12	ESR1	26
13	IL10	26
14	CXCL8	25
15	PIK3R1	24
16	CCND1	24
17	IL1B	24
18	MAPK14	24
19	IL4	23
20	ICAM1	22

GO analysis of PTTNs was classified into three groups: biological process, cellular component, and molecular function. Among the PTTNs (1F) were ones taking part in the biological process of apoptosis and intracellular signaling cascades. The PTTNs were located in the extracellular space and cytosol. Moreover, the PTTNs played a molecular function in growth factor and transcription factor activity. Pathway enrichment analysis using KEGG of the NMPs (Table S5) showed the regulation of approximately 53 pathways. Among these are pathways in cancer, Wnt/ß-catenin signaling, and TGF-beta signaling.

In exploring genomic alterations in breast cancer, a total of 8 PTTNs were analyzed using cBioportal. Regulators of cancer and the Wnt pathway are PTTNs EGFR, AKT1, WNT5A, MAPK8, JUN, TP53, and IL6, whereas ESR1 is the only hub protein from DTP based on highest degree score. Due to high alterations (60%), Lefebvre et al. (2016) was selected for further analysis among 12 breast cancer studies (Fig. 1G; (Lefebvre et al., 2016). Moreover, OncoPrint analysis showed that genetic alterations of the 8 PTTNs occurred in 2.8% to 42% of patients (Fig. 1H), and mutation is the most common gene alteration. From a study by INSERM 2016, additional mutual exclusivity analysis showed that gene pairs MAPK8-JUN, ESR1-TP53, EGFR-IL6, WNT5A-MAPK8, EGFR-MAPK8, WNT5A-JUN, and EGFR-JUN exhibited significant (p < 0.05) co-occurrence in breast cancer samples (Table 3). Further analysis of the relationship between the 7 PTTN and the neighboring genes can be found in the INSERM 2016 study (Fig. 1I). Moreover, TP53, ESR1, and EGFR were the main target of most cancer drugs, indicating the potential of TP53, ESR1, and EGFR as naringenin targets in BCSCs (Fig. 1J)

Table 3

Mutual exclusivity analysis of selected genes in metastatic breast cancer study.

A	B	Log2 Odds Ratio	p-Value	Tendency
MAPK8	JUN	>3	<0.001	Co-occurrence
ESR1	TP53	−2.583	<0.001	Mutual exclusivity
EGFR	IL6	>3	<0.001	Co-occurrence
WNT5A	MAPK8	>3	<0.001	Co-occurrence
EGFR	MAPK8	>3	0.001	Co-occurrence
WNT5A	JUN	>3	0.001	Co-occurrence
EGFR	JUN	2.67	0.029	Co-occurrence
IL6	TP53	2.85	0.049	Co-occurrence
AKT1	IL6	2.807	0.065	Co-occurrence
AKT1	TP53	1.292	0.074	Co-occurrence
AKT1	EGFR	1.731	0.103	Co-occurrence
AKT1	MAPK8	2.208	0.111	Co-occurrence
WNT5A	IL6	>3	0.157	Co-occurrence
ESR1	AKT1	1.095	0.181	Co-occurrence
ESR1	EGFR	1.095	0.181	Co-occurrence
ESR1	IL6	1.642	0.207	Co-occurrence
JUN	IL6	2.159	0.25	Co-occurrence
MAPK8	TP53	−1.166	0.268	Mutual exclusivity
JUN	TP53	−0.799	0.326	Mutual exclusivity
EGFR	TP53	0.495	0.342	Co-occurrence
AKT1	WNT5A	1.379	0.373	Co-occurrence
EGFR	WNT5A	1.379	0.373	Co-occurrence
ESR1	JUN	0.61	0.434	Co-occurrence
WNT5A	TP53	0.471	0.5	Co-occurrence
AKT1	JUN	0.501	0.545	Co-occurrence
ESR1	WNT5A	0.263	0.61	Co-occurrence
ESR1	MAPK8	−0.239	0.678	Mutual exclusivity
MAPK8	IL6	<−3	0.795	Mutual exclusivity

Enrichment of BCSCs

As described in Section 2, BCSCs enrichment was conducted by generating mammospheres using different methods from previous studies with modification. Cell preparation was initiated by harvesting previously cultured MCF-7 cells in DMEM (Fig. 2A). Cells (40,000 cells/well) were seeded in 50 mg/mL polyHEMA-coated 24-well plate, in SFM supplemented with EGF, insulin, and B27. Mammosphere was stable until day 7 after seeding.

Fig. 2

(A). Generation of mammospheres from MCF-7 cells. Formation of mammospheres from MCF-7 cells in serum-free media, as described in the methods section. Cells were seeded in a poly HEMA-coated plate (B). Identification of CD44+/CD24− population in MCF-7 cells using flow cytometry. (C). Quantification of CD44+/CD24− population in MCF-7 cells using flow cytometry. Cells in Q1-UL corresponded to CD44+/CD24− cells. An isotype control of MCF-7 cells (2D) was shown. (D). Gene expression of apoptosis, cell cycle, EMT, and stemness markers in 2D and 3D cells. Gene expression was determined by q-RT PCR. GAPDH was used as an internal control. The results were analyzed using the comparative threshold cycle (ΔΔCT) and are presented as fold change to the 2D cells. (E). Cytotoxicity of naringenin in MCF-7 cell monolayer (2D) and mammospheres (3D). A total of 10,000 cells per well were seeded in 96-well plates, as described in the method section. Cells were treated with naringenin and incubated for 72 h. MTT assay is used to determine cell viability. (F). Naringenin inhibits mammosphere formation from MCF-7 cells. Cells were seeded in 6-well plates and then incubated to 80% confluence. Subsequently, cells were treated with samples for 72 h. After that, the media was replaced with new media, and the cells were incubated for a further 24 h. Furthermore, 10,000 cells per well were seeded in 96-well plates and incubated for five days. At the end of incubation, the number of mammospheres formed is then calculated manually and analyzed as MFP. (G). Naringenin inhibits colony formation in MCF-7 cells. Naringenin treatment was given at 500 cells/well for 72 h, then the culture medium was replaced, and cells were grown for 14 days and stained at the end of the incubation period using gentian violet. The number of colonies was quantified with ColonyArea. Results represent the average of three independent experiments (mean ± SD). Two-way ANOVA with post-hoc Bonferroni’s multiple comparisons test was used to analyze the CD44+/CD24− population enrichment assay and cytotoxic effect of naringenin, while one-way ANOVA with post-hoc Bonferroni’s multiple comparisons test was used to analyze the effect of naringenin on mammosphere and colony formation. *** or **** indicates p < 0.001 or p < 0.0001, respectively. In Fig. 2C **** indicates significant to MCF-7 2D (p < 0.0001), while ++ in indicates significant to day 3 (p < 0.01).

CD44+/CD24− was selected as a molecular marker to identify the population of BCSCs and examine the stem cell enrichment characteristics in the mammosphere. Q1-UL area in Fig. 2B showed cell population in CD44+/CD24−. Quantification of the cell population showed a significant increase of CD44+/CD24− population after 3 days, a further rise after 5 days, and was stable after 7 days (Fig. 2C, indicating that CD44+/CD24− population at day 5 was significant than that on day 3 but was not significant on day 7. In addition, Gene expression of BCL2 and ALDH1 were significantly upregulated in 3D cells (Fig. 2D), however, p53 and ESR1 were significantly downregulated. Taken together, BCSCs were enriched in mammospheres from MCF-7 cells. Therefore, the authors decided to use mammospheres for the subsequent naringenin treatment after day 5.

Naringenin effects on cell viability and colony and mammosphere formation

After successful formation and characterization, mammospheres were treated with naringenin. Measured using MTT, with a parameter of cell viability, naringenin showed cytotoxicity in MCF-7 3D cells (Fig. 2E). Based on the number of mammospheres per cell seeded or the MFP, the effect of naringenin on mammosphere formation inhibition was analyzed by the number of mammospheres formed after treatment. Results revealed that naringenin inhibited mammosphere formation from MCF-7 cells (Fig. 2F). Moreover, the MFP value of naringenin is more or less similar to metformin, which targets cancer stem cell.

The drugs’ effect on colony formation was examined using the colony formation assay. Preliminary experiment reveal that naringenin concentration for colony and mammosphere formation was respectively 200 and 100 μM, which produced>80% cell viability (Fig. 2E) to ensure that cell growth inhibition does not influence the colony's inhibitory effect and mammosphere formation. The effect of naringenin on cell proliferation or cancer colonization was observed for 14 days; and after>24 h of observation, naringenin inhibited colony formation similar to that of doxorubicin (Fig. 2G).

Effects of naringenin on cell cycle arrest, apoptosis, and migration

The effect of naringenin in modulating the cell cycle was tested using flow cytometry with propidium iodide (PI) staining on DNA. In G1 phase, an increase in 3D cells was observed compare with 2D cells (Fig. 3A); naringenin increased G1 cell cycle arrest in 2D cells. The cell cycle profile of naringenin-treated 3D cells is not significantly different from the untreated cells. Flow cytometry with annexin V and PI reagents was used to carry out confirmation of cell death upon naringenin treatment. Naringenin was shown to induce apoptosis in both 2D and 3D cell culture (Fig. 3B). Moreover, the effect of apoptosis induction in 2D cells was higher than in 3D cells. One indication of BCSC is migration, and a scratch wound-healing assay was performed to measure the potency of naringenin to inhibit the migration of mammosphere-derived MCF-7 cells. Fig. 2E shows a naringenin concentration of 100 μM produced>80% cell viability, ensuring that cell growth inhibition does not influence the inhibitory effect of migration. The results in Fig. 3C showed that naringenin inhibited cell migration, which agrees with the present study that naringenin inhibits proliferation and migration and induces apoptosis in BCSCs.

Fig. 3

The effect of naringenin on the (A). Cell cycle profile. Cells were harvested after naringenin treatment for 72 h, stained with propidium iodide reagents, before analysis of DNA content using flow cytometry. (B). Apoptosis induction. Cells were harvested after naringenin treatment for 120 h, stained with annexin V and propidium iodide reagents, incubated, and analyzed using flow cytometry. The total percentage of cells consists of living cells, and cells undergoing early apoptosis, late apoptosis, and necrosis. (C). Naringenin inhibits migration. Mammosphere-derived MCF-7 cells were seeded and incubated for 24 h and starved with serum-free medium for another 24 h. After starvation, the cells were scratched using a sterile pipette tip and treated with naringenin. Images of the cells were captured at 0, 18, 24, 42, and 48 h after treatment. The results were analyzed using ImageJ and presented as percentage closure (n = 6). (D). Gene expression of EMT and stemness regulators in 2D cells. (E). Gene expression of apoptosis and cell cycle regulators upon naringenin treatment in 2D cells. (F). Gene expression of EMT and stemness regulators in 3D cells. (G). Gene expression of apoptosis and cell cycle regulators upon naringenin treatment in 3D cells. The effect of naringenin on p53 and ESR1 gene expression in 2D (H) and 3D cells (I). Gene expression was determined by q-RT PCR. GAPDH was used as an internal control. The results were analyzed using a comparative threshold cycle (ΔΔCT) and presented as fold change to the untreated control. Results represent the average of three independent experiments (mean ± SD). Statistical analyses were conducted using Student's t-test. *, **, and **** indicate p < 0.05, p < 0.01, and p < 0.0001, respectively.

Alterations of gene expression upon naringenin treatment

q-RT PCR was performed to examine the molecular mechanism of naringenin on BCSCs. The expression of stemness regulatory gene ß-catenin and ALDH1 was affected by naringenin treatment, as well as the epithelial-mesenchymal transition (EMT) regulators E-cadherin and Vimentin and migration regulator MMP9. The significant downregulation of ß-catenin and MMP9 and the reduced the expression of ALDH1 in 2D cells shown in Fig. 3D was due to naringenin treatment. Moreover, naringenin also increased Bax expression shown in Fig. 3E. In 3D cells, naringenin significantly reduced the expression of Vimentin and ALDH1 and BCL2 (Fig. 3F,G). In addition, naringenin reduced the expression of survivin and p21, but not statistically significant (Fig. 3G), and significantly increased p53 and ESR1 gene expression in 2D and 3D cells (Fig. 3I,J).

Discussion

This study aimed to explore the effect of naringenin in BCSC inhibition using bioinformatics approach and in vitro experiments. Integrated bioinformatics analysis retrieved DTPs, ITPs, and BCSCs-related genes, which resulted in the PTTNs in BCSCs. ERα is the only DTP having the highest degree score among hub proteins, indicating the pivotal role of ERα in the effect of naringenin on BCSCs. A similar bioinformatics work was published, where the association between gene-prostate cancer and naringenin by functional network analysis was identified (Fu et al., 2019).

The present study uses a similar approach in identifying key genes in naringenin effects against BCSCs. Hub proteins include EGFR, MAPK8, PIK3R1, Cyclin D1, and p53. This result is supported by previous studies, which showed that MAPK is one of the molecular targets of naringenin in MCF-7 cells (Eanes and Patel, 2016). A further study demonstrated that naringenin induces apoptosis in PC3 and LNCaP prostate cancer cells through inhibition of the PI3K/AKT and MAPK pathways (Lim et al., 2017).

A GO enrichment analysis revealed that PTTNs participated in the biological processes of apoptosis and intracellular signaling cascades. The PTTNs are located in the extracellular space and cytosol. Moreover, the PTTNs play a molecular function in growth factor and transcription factor activity. KEGG’s pathway enrichment of the PTTNs showed the regulation of approximately 53 pathways, including cancer, Wnt, and TGF-beta signaling pathways. Wnt signaling plays an important role in the proliferation, differentiation, invasion, migration, metastasis, and metabolism of cancer cells (Asem et al., 2016).

CSC formation can be achieved using 3D cell culture (Zhang et al., 2019), and to enrich BCSC properties, mammospheres were utilized. A serum-free medium in a polyHEMA-coated plate, with a diameter of around 100 µm, was used in culturing to generate the mammospheres. Based on our study, mammospheres were stable until day 7 and examination using CD44+/CD24- revealed that mammospheres BCSC properties were maintained. Mammospheres after day 5 were used for subsequent naringenin treatment because BCSC population at day 5 is not significantly different from the population at day 7. Results from this study, supported by a previous study, demonstrated that culturing spheroids to over 5–7 days will result to a growth of 100–500 µm in diameter (Liu et al., 2016). The effects of naringenin to the drug resistance phenotype of BCSCs, including drug uptake and MDR gene expression, needs further studies.

Naringenin exerts cytotoxicity in 2D and 3D cells, but its effects are stronger in 2D than in 3D cells. Previous studies showed that the IC50 value of naringenin is high. Naringenin showed cytotoxicity in MCF-7 (Fitriasari et al., 2010) and T47D (Junedi et al., 2010) breast cancer cells with IC50 values of 400 and 500 μM, respectively. The IC50 values of naringenin on MCF-7, HT-29, PC-12, and L-929 ranges from 780 to 880 µM (Kocyigit et al., 2016). According to Hakimuddin et al. (2004), Naringenin showed selectivity toward cancer cells (Hakimuddin et al., 2004). Naringenin treatment on CRL1554 human fibroblast cells, with a concentration of 4 mM for 24 h, induced cytotoxicity by about 30% (Abaza et al., 2015). In addition, Salehi et al. discussed clinical trials of naringenin on bioavailability and its effect on cardiovascular systems with a dosage between 600 and 800 μM/day (Salehi et al., 2019). Generally, the concentration of 100 and 200 µM is not toxic toward normal cells. A concentration of 100–500 µM was chosen for this study to determine cytotoxicity and a maximum dose of 200 µM was chosen for the other experiments on BCSCs (i.e., mammosphere forming analysis, clonogenic assay, cell cycle, apoptosis, anti-migratory assay, and q-RT PCR experiments).

Naringenin inhibited mammosphere formation based on MFP value. MFP is an accessible and suitable assay in estimating BSCS behavior (Lombardo et al., 2015) and it represents the population of CSCs responsible for chemoresistance and disease recurrence (Ji et al., 2016b). In this study, naringenin inhibited colony formation similar to doxorubicin (Fig. 2G). The colony formation or clonogenic assay is an in vitro technique that measures the ability of a single cell to grow into a larger colony, and it is a sensitive indicator of an undifferentiated CSC (Rajendran and Jain, 2018). To validate the results of this study, a study using methylcellulose based clonogenic assay is needed to replate the surviving colonies in the primary clonogenic assay and check the secondary colony formation potential.

The importance of p53 and ERα in the effect of naringenin on BCSCs were emphasized int the OncoPrint analysis, and additional mutual exclusivity analysis revealed the co-occurrence of p53 and ERα. Gene network analysis showed the importance of TP53 among neighboring genes and that p53 and ERα are the most druggable targets. Aside from naringenin, drugs that target p53 are PRIMA-1-Me, RITA, PRIMA-1, and those that target ERα are tamoxifen, raloxifene, and fulvestrant. Naringenin is clinically proven safe to use as a daily dietary supplement in healthy subjects (Rebello et al., 2020), and it can prevent osteosarcoma relapse in post-surgery osteosarcoma patients (Zhang et al., 2018), which is advantageous compared with other drugs.

In vitro results showed that mRNA of p53 and ERα were downregulated in 3D cells (Fig. 2D). Previous studies supported these findings, which demonstrated that ERα is downregulated in luminal A breast cancer that is resistant to chemotherapy and displays estrogen-independent signaling (Morimoto et al., 2009), and that loss of p53 promotes the maintenance of BCSCs (Chiche et al., 2017) and promotes Myc-activation and stem cell reprogramming of progenitors (Santoro et al., 2019). In this study, naringenin treatment significantly increased p53 and ERα mRNA expression in 3D cells (Fig. 3I). Bioinformatics finding was supported by the in vitro results, which emphasizes the pivotal role of p53 and ERα as PTTNs on BCSCs. A study using monolayer cell culture by Kanno et al. (2005) revealed no correlation between naringenin cytotoxicity and p53 status in cancer cells (Kanno et al., 2005). Nevertheless, the correlation between naringenin cytotoxicity and p53 status in BCSCs and the protein level regulation of p53 and ERα upon naringenin treatment remains unclear and needs to be confirmed in future studies.

The KEGG pathway enrichment analysis showed that cancer, β-catenin, and TGF-β signaling pathways may be regulated by naringenin in BCSC. Naringenin inhibition of Wnt/ß-catenin pathway in 2D cells are indicated in this study. These findings were supported by Li et al., (2013) showing that through the downregulation of survivin, naringin, the glycoside form of naringenin, inhibited the Wnt/ß-catenin pathway in triple-negative breast cancer cells through (Li et al., 2013). Activation of Wnt/β-catenin signaling plays a critical role in mesenchymal stem cell aging, which involves interactions with the p21 pathway (Zhang et al., 2011). Previous studies demonstrated the crosstalk between p53, ERα, and Wnt/ß-catenin signaling. Sadot et al., (2001) showed that activation of p53 leads to the downregulation of ß-catenin signaling (Sadot et al., 2001). In addition, the loss of wild-type p53 expression is driven by the excessive accumulation of ß-catenin (Sadot et al., 2001). Gupta et al., (2011) demonstrated that ß-catenin is a positive regulator of ERα function in breast cancer cells (Gupta et al., 2011). Moreover, the protein level of ERα is correlated with DKK1, a target gene of Wnt/ß-catenin signaling, highlighting the crosstalk between ERα and ß-catenin signaling in endometrial cancer (Kasoha et al., 2019). Generally, naringenin’s effect on the Wnt/ß-catenin pathway, p53, and ERα in inhibiting BCSCs is highlighted as an interesting topic for future research.

In the present study, TGF-β signaling is also inhibited by naringenin in 3D cells. Through increased expression of Vimentin, a marker of EMT, TGF-β activation triggers an EMT in breast cancer cells (Yoshida et al., 2013). EMT is important for the acquisition and maintenance of stem cell properties and CSCs. CSCs often exhibit EMT properties (Liu and Fan, 2015). In this study, Vimentin was upregulated in 3D cells compared with 2D cells, and this was also supported by previous studies that reported vimentin upregulation in BCSCs (Makki et al., 2015). The expression of Vimentin in 3D cells decreases with naringenin treatment. Results presented in this study demonstrated that naringenin inhibits cell migration. MMP9 was downregulated upon naringenin treatment in 2D cells; however, MMP9 expression was not downregulated upon naringenin treatment in 3D cells. This may be due to MMP9 activity inhibition rather than MMP9 expression. Previous studies supported this hypothesis which showed inhibition of MMP9 activity upon naringenin treatment in A549 lung cancer cells via regulation of the NFkB pathway (Chang et al., 2017) and in GBM 8901 glioblastoma cells through modulation of the ERK and p38 pathways (Chen et al., 2019). A study by Lou et al., (2012) demonstrated that by inhibiting TGF-β-induced EMT, naringenin reduced migration and invasion in pancreatic cancer cells (Lou et al., 2012).

TGF-ß signaling is also regulated by p53 and ERα. The transcriptional activity of TGF-ß target genes is suppressed by a p53 mutant gene (Elston and Inman, 2012); one of which is the p21 gene, which requires p53 for its transcriptional regulation (Cordenonsi et al., 2003). Moreover, p53 function inactivation modulates TGF-ß signaling (Kawarada et al., 2016). Stope et al. (2010) showed that ERα signaling activation inhibits TGF-ß signaling in breast cancer cells (Stope et al., 2010). In general, the mechanism of naringenin in BCSCs by regulating TGF-ß signaling, p53, and ERα is interesting topic to be evaluated.

The results of this study indicated that naringenin induces G1 cell cycle arrest in 2D cells, but cell cycle does not change compared with controls. In addition, it demonstrated that naringenin inhibits BCSCs by downregulating mRNA levels of ALDH1 and p21 in 3D cells. In correlation with Meng et al. (2014), it showed the importance of ALDH1 and p21 interaction in maintaining the CSC properties of ovarian cancer, in which downregulation of ALDH1 regulates cell cycle by decreasing p21 expression (Meng et al., 2014). Therefore, the interaction between ALDH1 and p21 and treatment with naringenin in mammospheres are interesting topics for further study.

The present study was done using the MCF-7 breast cancer cell line, which is a luminal subtype and has wild-type p53. In previous studies, naringenin exhibited cytotoxicity on different subtypes of breast cancer, including luminal, HER2+, and triple-negative breast cancer and it showed anticancer activity in triple-negative breast cancer cells 4 T1 (Qin et al., 2011) and MDA-MB-231 (Filho et al., 2014, Wang et al., 2019). Naringenin also exhibited cytotoxicity in MCF-7 breast cancer cells (Eanes and Patel, 2016, Xu et al., 2018) and HER2 positive breast cancer cells (Chandrika et al., 2016). In Kanno et al. (2005), cytotoxicity in various subtypes of monolayer breast cancer cells did not show any difference in naringenin treatment (Kanno et al., 2005). Breast cancer is a complex and heterogeneous disease; therefore, further study using the mammosphere from luminal breast cancer with mutant p53 or other types of breast cancer cell types is needed.

Several studies have shown the development of databases for predicting in vitro dose extrapolation to in vivo or clinical trials (Yoshida et al., 2017, Fabian et al., 2019). In addition, databases like ADMETlab (Dong et al., 2018) and SwissADME (Daina et al., 2017) can also be used for dose adjustments related to absorption, distribution, metabolism, excretion of a drug. A pharmacokinetic study of naringenin revealed the pharmacokinetic profile of naringenin in humans and rats (Bai et al., 2020). Taken together, the databases and results of the pharmacokinetics study can be used as basis in determining the clinical trial dose of naringenin in breast cancer patients.

One limitation of naringenin development is poor bioavailability. A previous study showed that the bioavailability of naringenin is 4% in oral administration (Hsiu et al., 2002). Naringenin undergoes further metabolism into several metabolites, including mainly sulfates and glucoronides (Zeng et al., 2019). Pharmacokinetics studies on rats, dogs, and humans revealed that apigenin, eriodictyol, hesperetin, naringenin-o-glucoronides, naringenin-o-glucosides, and 5,7-dihydroxychromone are naringenin metabolites (Bai et al., 2020). Another previous study showed that metabolites such as apigenin inhibited CSC properties in breast (Kim et al., 2016, Li et al., 2018) and prostate cancer (Erdogan et al., 2016) and eriodictyol in glioblastoma cells (Li et al., 2020). Therefore, the effect of naringenin on BCSCs inhibition in vivo might be due to its metabolites. However, this hypothesis need further investigation.

Studies were carried out to improve naringenin bioavailability. The naringenin nanoparticles formulation was shown to increase the physicochemical properties and hepatoprotective effects in vivo (Yen et al., 2009) and the cytotoxicity of naringenin in HeLa cervical cancer cells (Krishnakumar et al., 2011). Solid lipid nanoparticle naringenin was shown to increase the bioavailability of naringenin in the pulmonary delivery system in vitro and in vivo (Ji et al., 2016a). Therefore, further studies of the effects of naringenin nanoparticles on BCSCs are warranted.

In summary, this study highlighted the potential therapeutic use of naringenin for BCSC inhibition. The critical role of p53 and ERα in the effects of naringenin on BCSCs was emphasized through bioinformatics analysis and in vitro results. KEGG pathway enrichment analysis revealed TGF-ß and Wnt/ß-catenin pathways regulation by PTTN, p53 and ESR1 mRNA downregulation in mammospheres, and significant increased p53 and ESR1 mRNA levels with naringenin treatment. However, these results need to be confirmed at protein expression level. Naringenin has a more significant effect on mammosphere and colony formation. Validation of inhibiting those key genes could reverse the anti-proliferative, anti-migratory, or proapoptotic action of naringenin in BCSCs is needed. This study has a limitation as it is challenging to have a physiological significance because in vivo data is not available. Data presented in the manuscript is derived from a single cell line that needs to complement another cell line. Molecular mechanism of naringenin, including p53 and ER, in the inhibition of BCSCs requires further investigation for its development as a BCSC-targeted drug.

Conclusions

The results in this study showed that naringenin exhibits cytotoxicity in mammospheres and inhibits EMT and migration. More importantly, PTTNs that play a pivotal role in the effects of naringenin on BCSCs are p53 and ERα. Naringenin has the potential to eliminate BCSCs. Further in vitro and in vivo studies on naringenin’s molecular mechanism in BCSCs, as well as formulation for increasing bioavailability of naringenin, and clinical trials would be beneficial for its development as a BCSC-targeted drug.

Compliance with ethical standards

This study does not contain any studies with human participants or animals performed by any of the authors.

Availability of data and material

All data and material could be available from the corresponding author upon request.

Competing interests

The authors declare that they have no competing interests.

Funding

This work was supported by the Penelitian Unggulan Perguruan Tinggi (PUPT) 2017 and 2018, Contract No. 2398/UN1.P.III/DIT-LIT/LT/2017) and No. 189/UN1/DITLIT/DIT-LIT/LT/2018.

Authors' contributions

AH—conception and design of the study, acquisition, analysis and interpretation of data, drafting, and revising the article, MI, RIJ, AK, IPN, HP, SMA, and HAM—acquisition and analysis of data. All authors had final approval of the submitted manuscript.

Acknowledgements

The authors thank Badan Penerbit dan Publikasi Universitas Gadjah Mada for their assistance in writing.