Network-Based Approach to Identify the Antiproliferative Mechanisms of Bruceine D in Breast Cancer From the Cancer Genome Atlas

Bruceine D (BD) is a natural compound extracted from a Chinese herb Brucea javanica that has been used for anti-inflammatory and anti-cancer treatment. However, little is reported about BD's effects in breast cancer tumorigenesis. In this paper, we aimed to investigate the effect of BD in breast cancer and elucidate the potential mechanism of BD by integrated multiple databases. Our data suggested BD inhibited MCF-7 and MDA-MB-231 cells proliferation and promoted cells apoptosis. We integrated multiple bioinformatics analysis strategies to identify genes, hub modules and pathways associated with BD treatment. Three key pathways, including AMIT_SERUM_RESPONSE_40_MCF10A, BILD_HRAS_ONCOGENIC_SIGNATURE, and NAGASHIMA_NRG1_SIGNALING_UP were identified to be associated with therapeutic effects of BD in breast cancer. Additionally, we validated the key genes by using quantitative real-time PCR and western blot. In conclusion, these findings revealed potential molecular mechanisms of BD to treat breast cancer by affecting AMIT_SERUM_RESPONSE_40_MCF10A, BILD_HRAS_ONCOGENIC_SIGNATURE, and NAGASHIMA_NRG1_SIGNALING_UP pathways and regulating expression of ZFP36, EGR1, and FOS.


INTRODUCTION
Breast cancer is one of the most common malignancies with high rates of morbidity and mortality in women (1)(2)(3)(4). Over the past 20 years, the incidence rate of breast cancer is increasing in China (5)(6)(7). Although, several therapeutic strategies, including radiation therapy, chemotherapy, surgery and molecule-targeted therapy, improve survival rates greatly, many patients still face recurrence, drug resistance and serious side effects (8)(9)(10)(11). Therefore, it is urgent to identify a more effective therapy.
Recently, natural compounds have been considered as an important source of anti-tumor drugs (12,13). Bruceine D (BD) is a natural compound derived from a herb of Brucea javanica. It has been reported that BD has good pharmacological activities including cytotoxic effects on cancer cell, anti-inflammatory and hypoglycemic activities (2,(14)(15)(16)(17). However, little is reported about BD's effects in breast cancer. In our present study, we used breast cancer cells MCF-7 and MDA-MB-231 to test the effects of BD. Results illustrated that the intervention with BD significantly inhibits the proliferation and induces apoptosis of MCF-7 and MDA-MB-231 cells. In recent years, integrated bioinformatics analysis has been applied to understand molecular mechanisms of cancer and drug action mechanism (18,19). A large amount of genomic data was generated which has aided cancer studies considerably. Thus, we investigate the underlying pharmacological mechanism of BD in breast cancer by integrated bioinformatics analysis. We identified three pathways, including AMIT_SERUM_RESPONSE_40_MCF10A BILD_HRAS_ONCOGENIC_SIGNATURE and NAGASHIMA_NRG1_SIGNALING_UP pathways and three genes, including ZFP36, EGR1, and FOS. Finally, we adopt quantitative real-time PCR and western blot to validate the effects of BD on ZFP36, EGR1 and FOS expression in MCF-7 and MDA-MB-231 cells.
In this work, we firstly reported that BD significantly inhibited MCF-7 and MDA-MB-231 breast cancer cells proliferation and induced S phase cell cycle arrest and apoptosis in MCF-7 and MDA-MB-231 cells. Integrated computational analysis results indicated that BD mainly affects AMIT_SERUM_RESPONSE_40_MCF10A, BILD_HRAS_ONCOGENIC_SIGNATURE, and NAGASHIMA_NRG1_SIGNALING_UP pathways and regulates expression of ZFP36, EGR1, and FOS.

Cell Culture
Human breast cells MCF-7 and MDA-MB-231 were cultured in DMEM supplemented with 10% Fetal Bovine Serum (FBS), 100 µg/ml of streptomycin, and 100 U/ml of penicillin in a humidified incubator at 37 • C where the carbon dioxide content is 5%. The cells were grown to 70% confluence in monolayer and treated with indicated concentrations of BD for experiment, dissolved in dimethyl sulfoxide and subsequently diluted in cell culture medium.

Cell Cycle Distribution Assay and Apoptosis Assay
MCF-7 and MDA-MB-231 cells were treated with 0 µM (control), 0.1, 0.5, and 1.5 µM BD for 72 h. According to the instructions of the cell cycle detection kit, the cell cycle distributions were identified by PI staining and followed analyzed by flow cytometric. The BD's effect on cell apoptosis was evaluated using the Annexin-V/PI double-labeled flow cytometry kit. Cells were analyzed by flow cytometry with a FACS Calibur flow cytometer (BD Biosciences, San Jose, CA, USA).

Data Collection and Preprocessing
The mRNA raw count profiles of the TCGA-BRCA project were downloaded from the GDC data portal (https:// portal.gdc.cancer.gov). Total of 1109 BRCA patient samples and 113 control samples were available in TCGA. Simultaneously, the GSE85871 datasets were obtained from the Gene Expression Omnibus database (https://www.ncbi.nlm.nih.gov/geo). Raw CEL files of MCF-7 cells treated with BD (GSM2286384 and GSM2286384) and the respective controls (GSM2286316 and GSM2286317) were downloaded from the GSE85871 dataset which contained the gene expression profiles of 102 TCM ingredients used to treat MCF-7 cells. Raw CEL files were normalized by Robust Multiarray Average using affy package.

Differential Expression Analysis
Using DESeq2 package, we performed differential expression analysis at the cutoff of |log2 fold change| > 2 and padj < 0.05 (p-value adjusted using Benjamini-Hochberg method) in TCGA dataset. For GEO datasets, the limma package was used to screen the differential expression genes at the cutoff of |logFC| > 1 and padj < 0.05 (p-value adjusted using Benjamini-Hochberg method). In addition, a pre-ranked gene list of all genes ordered by the absolute value of fold change were constructed in GEO dataset for further analysis.

Establishment of Gene Co-expression Network
First, we normalized the counts data from TCGA by the TMM algorithm implemented in the edgeR package. Next, the weighted correlation network analysis (WGCNA) was conducted to construct the co-expression network in TCGA dataset using WGCNA package. The detailed procedure for WGCNA construction could be found in our previous study (20). Briefly, according to an appropriate β, we built a weighted adjacency matrix and transformed the adjacency into a topological overlap matrix (TOM). Finally, we performed an average linkage hierarchical clustering according to the TOM-based dissimilarity measure (20). In this study, we chose a minimum module size of 30 for the gene dendrogram and a cut-line (0.25) for the module dendrogram (20).

Determination of Significant Modules Associated With Breast Cancer
According to gene expression similarities in samples, gene modules were discerned. Simultaneously, using two methods, the correlation was calculated between the gene modules and sample information to find significant modules associated with breast cancer. First, we calculated the module eigengene (ME) of a module to represent the overall expression level of the module. Subsequently, we also calculated gene significance (GS) between gene expression and a sample trait, which was defined as the log10 transformation of the P-value (GS = lgP) in the linear regression (20). Finally, the gene module most correlated with sample traits were identified as a module of interest.

Differential Co-expression Analysis and Gene Set Enrichment Analysis
The gene set net correlations analysis (GSNCA) was conducted to find pathways that were differentially co-expressed between BRCA samples and control samples. There was an assumption that the co-expression network of a pathway did not alter between two conditions in GSNCA, which was implemented in GSAR package. Using the most correlated genes, it builds a core of coexpression network in each condition and finds a hub gene with the highest correlations with the other genes in a pathway. It hinted regulatory changes of the pathway when hub genes in a pathway were different between two conditions. The colors of nodes represent weight factor (w)'s value distributed to all genes for suggesting the mean correlation with any other gene in the gene set co-expression (21). In this paper, it is worth noting that pathways only with 10 ≤ n ≤500 where n represents the number of genes were kept. Then, using structure of minimum spanning tree-2 (MST2), the most highly correlated pathways were distinguished (21). In this study, the pathway gene sets were used in C2 collection (curated gene sets) from Molecular Signatures Database (MSigDB). In addition, a pre-ranked gene list of all genes ordered by fold change were constructed in BD samples vs. model samples. Then, we performed gene set enrichment analysis (GSEA) for C2 collection (curated gene sets) of MSigDB using clusterProfiler package. Simultaneously, other functional enrichment analysis also was did using clusterProfiler package at a cutoff of P < 0.05.

Quantitative Real-Time PCR Analysis
Human breast cancer cells MCF-7 and MDA-MB-231 were treated with BD at indicated doses for 72 h. The total RNA was extracted from MCF-7 and MDA-MB-231 cells using a commercial kit (TianGen Biotech Co., Ltd, China) according to the manufacturer's instructions, and 500 ng of RNA were reverse transcribed to cDNA using high capacity cDNA reverse transcription kit (TaKaRa Biotechnology Co., Ltd, Japan). The expression of ZFP36, EGR1 and FOS genes were determined by real-time PCR using the ABI Prism 7900 sequence detection system (Applied Biosystems Co., Ltd, USA) and SYBR PrimerScript RT-PCR Kit (TaKaRa Biotechnology Co., Ltd, Japan). The primers were: ZFP36 (Gene bank: 7538) sense, CTGTCACCCTCTGCCTTCTC; anti-sense, TCCCAGGGACT GTACAGAGG, EGR1 (Gene bank:1958) sense, TGACCGC AGAGTCTTTTCCT; anti-sense, TGGGTTGGTCATGCTCA CTA and FOS (Gene bank: 2535) sense, AGAATCCGAA GGGAAAGGAA; anti-sense, CTTCTCCTTCAGCAGGTTGG. To control variation in mRNA concentration, all results were normalized to the housekeeping gene, GAPDH.

Western Blot Analysis
Human breast cancer cells MCF-7 and MDA-MB-231 were grown in 6-wells plates up to 60% confluence, and treated with BD at indicated doses for 72 h. The cells were washed 2 times with phosphate buffered saline (PBS), lysed in 40 µL of cold RIPA lysis buffer containing 1X Phenylmethylsulfonyl Fluoride (PMSF) and then used for western blot assay. Western blot assay was performed as described previously (22). Membranes were probed using rabbit antibodies specific for ZFP36, EGR1, FOS, and GAPDH (1:1000), followed by incubation with HRPconjugated secondary antibody (1:2000).

Statistical Analysis
All statistical tests were did using Student's t-test or one-way analysis of variance (ANOVA) when appropriate. GraphPad and R software were used for statistical analysis of the original data and the data were expressed as the means ± standard deviation (SD). For all statistic al analysis, a P < 0.05 was considered statistically significant.

The Effect of BD on MCF-7 and MDA-MB-231 Breast Cancer Cells Proliferation
After BD treatment of MCF-7 and MDA-MB-231 (triple negative breast cancer) cells at different concentrations for 24, 48, and 72 h, we used CCK8 method to determine the inhibition rates of cell proliferation. As shown in Figure 1

Identification of Modules by Weighted Co-expression Network Analysis
Firstly, the genes for co-expression network analysis were selected by differential expression analysis in TCGA dataset. A total of 1844 DEGs (1141 up-regulated and 703 downregulated) were identified using TCGA dataset ( Figure 3A). Meanwhile, 1106 DEGs (707 up-regulated and 399 downregulated) were identified in BD samples vs. model samples in GSE85871 ( Figure 3B). Then, the 1844 DEGs were analyzed for co-expression modules by the average linking clustering implemented in WGCNA package. In this work, the power of β = 6 (scale free R 2 = 0.9) was chosen to ensure a scale-free network (Figures 4A-D). Then, the hierarchical clustering tree was determined by conducting hierarchical clustering for dissTOM and 6 modules were determined ( Figure 5A). Subsequently, two methods were used to test the relations between each module and disease condition. Modules with a larger MS were considered to have more relations with disease condition. We found that the ME in brown module displayed a highest correlation with disease condition (Figure 5B). In addition, the ME of the brown module also showed the highest gene significance ( Figure 5C). Finally, the brown module was selected for further analysis.

Identification of Key Pathways and Key Genes
According to the above analysis results, the hub module was selected to identify the key pathways and key genes. Finally, the AMIT_SERUM_RESPONSE_40_MCF10A, BILD_HRAS_ONCOGENIC_SIGNATURE, and NAGASHIMA_NRG1_SIGNALING_UP pathways were identified and ZFP36, EGR1, and FOS genes were identified as key genes by MST2 plot in BRCA samples (Figures 6A-C). Simultaneously, we also performed GSEA to identify the potential therapeutic mechanism of BD in GSE85871 dataset. Interesting, the pathways were also significantly enriched (Figure 7). In addition, the w-values of three key genes above are increased in BRCA samples, indicating their regulatory roles in breast cancer samples and their loss in the control samples (21). Although the w-values of the three genes was decreased in control samples, these genes were still close to each other in structure of MST2 in control samples, this suggests that correlation of them with other genes was not completely lost (21). Additionally, to further find the intervention biological processes of BD, we also performed functional enrichment analysis. The GO terms of biological process were most significantly enriched in regulation of lipid metabolic process, protein-DNA complex assembly, intracellular receptor signaling pathway, intrinsic apoptotic signaling pathway, regulation of interferon-beta production and DNA conformation change (Figure 8). Detail parameters of pathways and GO terms are shown in Tables 1-3. These results suggest the BD may affect these pathways and biological process and thus inhibited breast cancer cells proliferation and induced their apoptosis.

Validation of Key Genes by qRT-PCR and Western Blot Assays
To validate the effects of BD on key genes, we initially quantified its effect on the level of ZFP36, EGR1 and FOS genes mRNA. As shown in Figure 9A, at three concentration

DISCUSSION
Bruceine D, a tetracyclic triterpene quassinoid, is one of the active compounds of the Chinese herb Brucea javanica. Several studies have shown that BD exhibits good pharmacological activity including anti-inflammatory (23), anti-parasitic (24), hypoglycemic activities (25), and anti-tumor activity such as hepatocellular carcinoma (2), human chronic myeloid leukemia (26), and pancreatic cancer (27). However, by literature search, we found that therapeutic effects and mechanisms of BD have been rarely reported in breast cancer. In this study, we firstly demonstrated that BD significantly inhibited MCF-7 and MDA-MB-231 proliferation and induced S phase cell cycle arrest and apoptosis, which indicates BD could exhibit anti-breast cancer pharmacological activity. Then, we utilized computational analysis to investigate the mechanism of BD intervention in breast cancer. We built a weighted co-expression network and identified six modules from BRCA. And the brown module showed the best correlation with BRCA by correlation analyses. In addition, GSNCA and GSEA was conducted to identify key pathways and key genes to elucidate molecular mechanisms of BD. The results showed BD in breast cancer treatment may take actions through affecting AMIT_SERUM_RESPONSE_40_MCF10A, BILD_HRAS_ONCOGENIC_SIGNATURE, and NAGASHIMA_NRG1_SIGNALING_UP pathways and regulation expression of ZFP36, EGR1, and FOS(c-FOS). More importantly, the key pathways are closely associated with cancer. For example, it was reported that a number of genes within AMIT_SERUM_RESPONSE_40_MCF10A pathway are downregulated in diverse tumor types and have a correlation with clinical outcomes. BILD_HRAS_ONCOGENIC_SIGNATURE pathway can discriminate cells expressing activated HRAS from control cells expressing GFP and distinguish between     Additionally, some studies have revealed three key genes play an important role in cancer development. ZFP36 is a tandem zinc finger protein in TTP family, which can directly bind AUrich elements (ARE) and destabilize the host transcript. Several studies showed that ZFP36 binds to and degrades the GM-CSF, TNFα, and VEGF mRNAs, in turn affecting the function of transcription factors such as NF-κB, which has a vital impact on cell viability (28,29). In breast cancer, ZFP36 gene defects will contribute to cancer progression and development. Our results showed that BD could significantly increase ZFP36 mRNA and protein expression in a dose-dependent manner in MCF-7 and MDA-MB-231 breast cancer cells. EGR1 is a tumor suppressor and an oncogene (30)(31)(32). As a transcription factor, EGFR can active or repress gene transcription (33,34). In addition, the transcription of EGR1-mediated brings together a multitude of signaling cascades, which have a vital impact on cell growth, apoptosis and differentiation (33,35). Recent studies have demonstrated the in breast cancer EGR1 defects are associated with poor prognosis, which is independent of subtype (36). In accordance with this context, we illustrated that EGR1 may be a BD's target for breast cancer treatment. FOS (c-FOS) is one of FOS family protein, which contains transcription factors such as Fra-1, Fra-2, FosB, FosB2, and FosB (37). FOS plays a relevant role in the regulation of normal cell growth, differentiation, and cellular transformation processes (38,39). The c-Jun family members can interact with nuclear protein, which was encoded by FOS, to form the heterodimeric activating protein-1 (AP-1) transcription factor complex (40). This complex binds to DNA at its specific sites in promoters and enhancer regions and causes changes in gene expression patterns, which might be critical for cell cycle progression. Normally, FOS protein plays a major role in several cellular functions, and its overexpression has been observed in many neoplasias (41). In the present study, we demonstrated that BD significantly decreased the level of FOS mRNA and protein expression in MCF-7 and MDA-MB-231 breast cancer cells. Hence, these key genes are likely to be potential targets of BD. However, there were some limitations in our study. Above all, in our paper, we just chose differentially expressed genes for weighted correlation network analysis. This did not include all the genes that may be associated with breast cancer. Then, there was only two biological replicates in GSE85871 dataset, which may affect the accuracy of the results. Finally, experimental studies need to be carried out to verify the three potential targets of BD in future work.
In summary, by integrated multiple computation methods, we found that the BD may influence some key pathways that are closely associated with cancer, such as AMIT_SERUM_RESPONSE_40_MCF10A, BILD_HRAS_ONCOGENIC_SIGNATURE, and NAGASHIMA_NRG1_SIGNALING_UP pathways. Simultaneously, the key genes, including ZFP36, EGR1, and FOS may be potentially effective targets of BD and also are crucial for breast cancer development. Additionally, this is the first work to report therapeutic effects and possible molecular mechanisms of BD in breast cancer which would future improve our understanding of BD. The effects of BD on FOS, ZFP36, and EGR1 protein expression in MCF-7 and MDA-MB-231 breast cancer cells were detected by western blot. *p < 0.05 to control group; **p < 0.01 to control group; ***p < 0.001 to control group (0 µM).

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study, these can be found in The Cancer Genome Atlas (https://portal. gdc.cancer.gov/) (TCGA-BRCA); the NCBI Gene Expression Omnibus (GSE85871).

AUTHOR CONTRIBUTIONS
ST and RJ collected and analyzed the data, drafted, and revised the manuscript. WZ designed the study and revised the manuscript. All authors contributed to the article and approved the submitted version.