CXCL1 Clone Evolution Induced by the HDAC Inhibitor Belinostat Might Be a Favorable Prognostic Indicator in Triple-Negative Breast Cancer

Triple-negative breast cancer (TNBC) is the most lethal subtype of breast cancer due to its lack of treatment options. Patients with TNBC frequently develop resistance to chemotherapy. As epigenetic-based antineoplastic drugs, histone deacetylase inhibitors (HDACis) have achieved particular efficacy in lymphoma but are less efficacious in solid tumors, and the resistance mechanism remains poorly understood. In this study, the GSE129944 microarray dataset from the Gene Expression Omnibus database was downloaded, and fold changes at the transcriptome level of a TNBC line (MDA-MB-231) after treatment with belinostat were identified. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were used to identify the critical biological processes. Construction and analysis of the protein-protein interaction (PPI) network were performed to screen candidate genes related to cancer prognosis. A total of 465 DEGs were identified, including 240 downregulated and 225 upregulated genes. The cytokine-cytokine receptor pathway was identified as being significantly changed. Furthermore, the expression of CXCL1 was implicated as a favorable factor in the overall survival of breast cancer patients. With in vitro approaches, we also showed that belinostat could induce the expression of CXCL1 in another 2 TNBC cell lines (BT-549 and HCC-1937). We speculate that belinostat-induced CXCL1 expression could be one of the results of the stress clone evolution of cells after HDACi treatment. These findings provide new insights into clone evolution during HDACi treatment, which might guide us to a novel perspective that various mutation-targeted treatments should be implemented during the whole treatment cycle.

As the first successful application of epigenetic-based cancer therapy, HDAC inhibitors have been discovered to have specific anticancer activities in preclinical studies and clinical treatments [5]. In 2006, suberoylanilide hydroxamic acid (SAHA, vorinostat) was approved by the US FDA to treat cutaneous T-cell lymphoma (CTCL). After that, several HDACis have been approved by the FDA for the treatment of other cancers, including peripheral T-cell lymphoma (PTCL) and multiple myeloma [6,7]. Until now, epigenetic therapy has achieved specific effects in hematologic neoplasia, which has stimulated growing interest and enthusiasm for further developing epigenetic therapies for other malignancies.
In addition to some progression in T-cell lymphomas and other hematological malignancies, several HDAC inhibitors are efficacious in solid tumors. The results of most clinical trials were in favor of using HDAC inhibitors either before the initiation of chemotherapy or in combination with other treatments [8]. However, their overall effects on chromatin, which can be viewed as positively modulating the expression of many genes, are also likely to generate clinical toxicities, limiting their clinical use. It has been reported that limited penetration and extensive tissue distribution result in clinical ineffectiveness and off-target effects, such as myelosuppression, fatigue, and cardiac toxicity [9]. Hence, the molecular mechanism underlying the inadequate response of HDAC inhibition in solid tumors still needs to be elucidated.
Breast cancer can be divided into 3 subtypes based on the presence or absence of different proteins in breast cancer cells, which have different prognostic and therapeutic implications. Hormone receptor-positive breast cancer accounts for approximately 70% of breast cancer cases and has either estrogen receptor (ER) or progesterone receptor (PR) protein in the cancer cells; human epidermal growth factor receptor 2 (HER2, also known as ERBB2) breast cancer makes up 15% to 20% of breast cancer cases; TNBC is more heterogeneous and lacks ER, PR, and HER2 protein expression, accounting for approximately 15% of all breast cancer cases [10]. Patients with TNBC show poor prognosis and frequently develop resistance to chemotherapy. Since they lack ER, PR, and HER2 receptors, they are not eligible for hormone or anti-HER2 therapy. TNBC patients harbor high levels of somatic mutations, frequent mutations in TP53 (83%), and complex aneuploid rearrangements (80%) that result in extensive intratumor heterogeneity (ITH) [11,12]. To date, few HDACis have been approved by the US FDA for breast cancer treatment. Notably, chidamide, which selectively targets class I HDACs (subtypes 1, 2, and 3) and class IIb HDACs (subtype 10), has been officially approved by the Chinese National Medical Products Administration (NMPA) for use in combination with aromatase inhibitors in the treatment of locally advanced or metastatic hormone receptorpositive breast cancer that overexpresses HER2 [13]. The benefits from the breakthroughs in HDAC inhibitor research have greatly encouraged the study of their mechanisms in solid tumors. Table 1 lists several HDAC inhibitors currently being evaluated in a few phase II/III clinical trials. Belinostat is an HDAC inhibitor that has been broadly used to treat PTCL [14]. As a hydroxamic acid-derived pan-HDAC inhibitor that has a high affinity for class I HDACs 1, 2, and 3, class II HDACs 6, 9, and 10, and class IV HDAC 11 [15,16], belinostat was also evaluated for the treatment of solid cancers, such as lung squamous cell carcinoma [17], renal cancer [18], and hepatocellular carcinoma [19]. Several clinical trials are based on targeting advanced breast cancer, including NCT04315233, NCT03432741, and NCT00413075. The column number of Table 1 should be corrected in sequence.
In recent years, the development of bioinformatics has extensively promoted progress in the field of life sciences. Integrating and reanalyzing the RNA-seq or microarray data using bioinformatics methods may help identify gene regulatory pathways, essential genes, and their associated networks in different diseases, providing information on the possible molecular mechanisms of diseases, potential drug research, and development directions [20,21]. Some studies have precisely revealed details on tumor microenvironment crosstalk [22,23] and the effect of epigenetic modification on tumor development [24]. In this study, we found that some cytokines, especially CXCL1, showed a significant difference in expression after treatment with belinostat. Previous literature has reported that CXCL1 is a critical factor in inflammatory diseases and tumor progression. It increases the expression of matrix metalloproteinase (MMP) 2/9 through the ERK1/2 pathway as well as breast cancer metastasis and invasion [25,26]. This study is aimed at better understanding the potential molecular alterations underlying HDAC inhibitors in breast cancer via bioinformatics methods and in vitro experiments, providing a rationale to explore the clinical value of HDACi in solid cancer.

Materials and Methods
2.1. Acquiring RNA Sequencing Profiles and Analysis of Differentially Expressed Genes. Expression profiles from high-throughput sequencing were obtained using the Gene Expression Omnibus (GEO) database. We downloaded the mRNA expression microarray dataset GSE129944, which contained data about MDA-MB-231 cells after treatment with 17-AAG, belinostat, and the combination of 17-AAG with belinostat. Each treatment group was performed in duplicate. The mRNA expression profiles were determined using the Illumina HiSeq 2500 mRNA sequencing platform (Illumina Inc., USA) and normalized as fragments per kilobase of exons per million mapped reads (FPKM) data [27]. Among those sets, we chose the belinostat treatment group to probe the molecular alterations responding to HDAC inhibitors. Gene expression alterations in the treated group were normalized to the corresponding control treated with vehicle. We used the ggplot2 package in R version 4.0 to explore the differentially expressed genes (DEGs) and adopted a threshold cutoff of p < 0:05 with absolute log 2fold change ð|log2 FC | Þ ≥ 1:5. To understand the functional roles of the differentially expressed mRNAs, DEGs were input into the DAVID (https://david.ncifcrf.gov) [28] and subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses [29,30]. To better understand the differential expression, upregulated and downregulated genes were analyzed separately. Enrichment analysis was carried out to measure the function's significance; the higher the value of enrichment, the more specific the corresponding function, by which the GO term of the associated biological process was identified. KEGG was used to understand the high-level functions and utilities of biological systems as previously mentioned [31]. We used the "topGO" and "pathview" R packages to implement the enrichment analyses. GO terms and KEGG pathways with a corrected p < 0:05 were considered significantly enriched.

Protein-Protein Interaction
Network. Based on the GO and KEGG pathway analysis, the STRING database (http:// string-db.org) was used to construct the protein-protein interaction network of differentially expressed genes [32]. Furthermore, Cytoscape v3.7 was used to identify the core motifs [33]. Through the MCODE function of Cytoscape, modules with the highest score were filtered out. Besides, using the cytoHubba plug-in in Cytoscape v3.7 software, the top 20 nodes were ranked by degree [34].

Hub Gene Screening and Survival Analysis.
After acquiring the hub genes, the prognosis was assessed via GEPIA (http://gepia.cancer-pku.cn/) [35]. As the hub genes involved in the GO analysis and pathway analysis were associated with the tumor's biological characteristics, we assessed their prognostic significance. Also, we used the online tool TIMER (https://cistrome.shinyapps.io/timer/) to investigate CXCL1 levels in different cancers [36], and the UALCAN website (http://ualcan.path.uab.edu/) was used to acquire the expression of CXCL1 in BRCA based on breast cancer subtypes [37]. 2.6. Real-Time Quantitative PCR. According to the manufacturer's protocol, total RNA was isolated from cells at the logarithmic phase using the TRIzol reagent (Sigma, USA). Firststrand cDNA was synthesized using the GoScript Reverse Transcription System Kit (Promega, USA). Real-time PCR was performed with GoTaq qPCR Master Mix (Promega, USA) using a C1000 Thermal Cycler apparatus (Bio-Rad) in a 20 μl reaction volume to the manufacturer's protocols. The procedure was as follows: 95°C for 3 min, (95°C for 15 s, 60°C for 60 s, and 72°C for 30 s), and 95°C for 10 s, followed by a melt curve analysis (60°C to 95°C, increments of 0.5°C for 20 s) to confirm the specificity of the PCR primers. Ct values for mRNA were normalized to GAPDH. The primers for CXCL1 and GAPDH were as follows:  overnight. According to the manufacturer's instructions, the proteins were detected by an enhanced chemiluminescence system (Pierce, USA). Three independent experiments were carried out, and the CXCL1 protein data were normalized to α-tubulin.
2.8. Statistical Analysis. Fold change and Student's t-test were employed to evaluate the statistical significance of the results. Differences with p < 0:05 between the two groups were considered significant. The p value was corrected by calculating the false discovery rate. GO analysis was performed via the DAVID with its statistical tool, the Bonferroni correction method, the Benjamini-Hochberg false discovery rate, and the bootstrap method. KEGG analysis was carried out by Fisher's exact probability test and the gene enrichment analysis included in the DAVID. Pearson correlation coefficients were used to construct the PPI network with Cytoscape software. Then, the gene expression data were clustered hierarchically, and the connection methods used were average linkage and median standardization. We profiled the general survival information and the transcripts using the univariate Cox proportional hazards regression model. The difference in OS between the high mRNA expression and low mRNA expression groups was assessed by Kaplan-Meier survival curves and the log-rank test. p values by two sides less than 0.05 were considered statistically significant. Real-time quantitative PCR and Western blot analysis were performed using one-way ANOVA.

Results
3.1. Data Summary of RNA Sequencing Identification. To evaluate the effect of the HDAC inhibitor belinostat on TNBC cells, the GSE129944 dataset and expression profiles from high-throughput sequencing were collected from the GEO database. We used absolute log2 ðfold changeÞ ≥ 1:5 and p value < 0.05 as cutoff values value. As expected, belinostat treatment led to alterations in the global transcription profile, and 465 DEGs were obtained from pairwise comparisons of samples, including 240 downregulated and 225 upregulated genes. As shown in Figure 1, the DEGs were visualized in the form of a volcano plot.

Belinostat Treatment Results in the Expressional
Reprogramming of MDA-MB-231 Cells. To further evaluate the transcriptome's functional distribution, the DAVID was used to perform GO and KEGG pathway enrichment analyses of the DEGs [38]. The GO enrichment results showed that the upregulated DEGs' biological processes were mainly involved the immune response, female pregnancy, cell-cell signaling, response to lipopolysaccharide, and negative regulation of endopeptidase activity. The downregulated DEGs were associated with cell adhesion, negative regulation of transcription from the RNA polymerase II promoter, extracellular matrix organization, protein ubiquitination involved in ubiquitin-dependent protein catabolic processes, and negative chemotaxis. Phosphorylation of the RNA polymerase II large subunit is necessary for initiation and elongation of transcription [39]. After the HDACi treatment, the negative regulation of transcription from the RNA polymerase II promoter may lead to cell death in preclinical models of TNBC. Regarding molecular functions, the upregulated DEGs were mainly related to cytokine activity, receptor binding, iron ion binding, growth factor activity, and cysteine-type endopeptidase inhibitor activity. The downregulated DEGs were mainly involved in transcription factor activity, sequencespecific DNA binding, chemorepellent activity, neuropilin binding, semaphorin receptor binding, and extracellular matrix binding. Additionally, the cell component analysis results suggested that the DEGs might be involved in the extracellular region and integral component of the membrane. GO analysis results have some similarities with the previous RNA-seq of TNBC [40], which indicated that the identified genes were associated with the regulation of various biological processes (as shown in Figures 2(a) and 2(b)).
Some studies have shown that the cytokine-cytokine receptor interaction pathway plays a critical role in generating an immune-suppressive microenvironment and participates in metastasis and proliferation [34]. Besides, the changes in pathways related to some immune diseases were statistically significant. These annotations provide a valuable resource for investigating biological pathways and gene functions.

Survival Curve and Bioinformatics Prediction.
In this study, the hub genes associated with the tumor's biological     Table 1). The results showed that the expression of CXCL1 was near related to superior OS in breast cancer patients. The survival curve suggested that OS was significantly shortened in breast cancer patients with low expression of CXCL1 ( Figure 5(a)). Furthermore, the UALCAN dataset, which is based on publicly available cancer omics data (TCGA and MET500), was used to estimate the expression of CXCL1 in BRCA based on different subclasses (Figure 5(b) and Supplementary Table 2). The statistical analysis showed that normal vs. TNBC and luminal vs. TNBC had significant differences (p < 0:05), but HER2-positive vs. TNBC did not have a noticeable difference. Through the analysis of immune infiltration across diverse cancer types, we found that compared to normal tissue, CXCL1 has low expression in breast cancer, while in cholangiocarcinoma,   BioMed Research International colon adenocarcinoma, esophageal carcinoma, head and neck squamous cell carcinoma, and hepatocellular liver carcinoma, the situation is the opposite (Figure 6). Commonly, different tumor types or different periods with the same type of tumor may be related to the unique immune microenvironment. CXCL1 has been reported to be associated with metastasis, angiogenesis, and chemoresistance [40,41]. The objective function of CXCL1 in tumor development and the alteration mechanisms and properties after the HDACi treatment still need to be elaborated.    (Figure 7(b)). The changes were statistically significant.

Belinostat Could
In summary, the induction of CXCL1 expression after belinostat treatment might be an indicator of superior prognosis for TNBC patients, which needs further attention for precision therapy in the setting of this life-threatening disorder.

Discussion
Although HDAC inhibitors have achieved some success in treating the hematological system's malignant tumors, the results in solid cancers in terms of benefits are less clear. Cur-rently, some clinical trials are underway and tend to combine HDAC inhibitors with chemotherapy or other targeted therapies to enhance clinical efficacy [42,43]. Generally, HDAC inhibitors lead to the inhibition of tumor growth and apoptosis of cancer cells. A previous study built an acetylation model to predict and verify the cellular metabolic state's impact on sensitivity to drugs that disrupt acetylation and demonstrated the interconnection between metabolism and acetylation [44]. Besides, class II HDAC inhibitors can selectively reprogram monocytes and macrophages in the tumor; reprogramming activates a robust antitumor immune response, mainly mediated by macrophages, CD8+ T-cells, and IFNγ, and reduces both primary and metastatic tumor burdens [45]. For example, a selective class IIa histone deacetylase inhibitor, TMP195, influenced human monocyte responses to the colony-stimulating factors CSF-1 and CSF-2 in vitro [46]. It alters the tumor microenvironment and reduces tumor burden and metastasis by modulating 9 BioMed Research International macrophage phenotypes, including the recruitment and differentiation of highly phagocytic and stimulatory macrophages within tumors [47]. Belinostat broadly inhibits class I HDACs 1, 2, and 3, class II HDACs 6, 9 and 10, and class IV HDAC 11. Indeed, in some human tumors, the overexpression of HDAC6 is associated with more advanced tumor stages and higher tumor invasiveness, so the survival rate is low in cholangiocarcinoma, ovarian cancer, and acute myeloid leukemia (AML) [48][49][50]. Besides, there might be a connection between the overexpression of HDAC6 in ERpositive cells and ineffective endocrine therapy and poor prognosis [51]. A previous study demonstrated that MDA-MB-231 cells show increased invasiveness and migration compared with MCF-7 cells because they overexpress HDAC6 and matrix metalloproteinase (MMP) 9 [52].
In this study, the cytokine-cytokine receptor pathway, as the top-ranked pathway, stimulated our interest. The DEGs involved in the interaction were transcriptionally modulated, suggesting that the reprogramming of the network in cancer cells was triggered by HDAC inhibition. According to previous reports, CXC chemokines are critical to malignant initiation and cancer progression, in addition to their role in inflammation. Some CXC chemokine family members act as promoters of angiogenesis, including CXCL1, CXCL2, CXCL3, CXCL5, CXCL6, CXCL7, and CXCL8 [53]. Based on the analysis of immune infiltrates across diverse cancer types, abnormal expression of CXCL1 has been found in numerous types of malignancies and has been associated with metastasis, angiogenesis, and chemoresistance [54,55]. Zeng et al. discovered that the mRNA level alteration of cytokine pathways caused by HDAC led to the downstream response via the LIFR-JAK1-STAT3 signaling-centered feedback loop, which restrained the efficacy of HDAC inhibitors in breast cancer [56].
Notably, the lack of response could be due to druginduced compensatory alterations arising in both malignant cells and the tumor microenvironment. A previous study reported the evolution of genetic mechanisms of resistance to palbociclib plus fulvestrant in ER-positive breast cancer. It showed that clonal evolution is frequent in response to therapy since acquired driver mutations in growth factor receptors and signal transduction pathways are frequently detected [57]. According to large-scale sequencing results, the level of CXCL1 in normal tissue is higher than that in breast cancer. In our study, CXCL1 was upregulated after treatment with HDACis and predicted a better prognosis, which is consistent with previous results. The CXCL1 mutation induced by belinostat treatment stress might result in DNA damage in tumor cells, which causes the microenvironment to become inadaptive to tumor cell growth. Drugmediated stress-induced mutagenesis probably acts as a double-edged sword, causing neoplasm necrosis in the context of belinostat intervention [58].
The selective pressure of treatment may lead to intrinsic tumor resistance and the acquisition of an adaptive response. It is not clear whether CXCL1 plays a positive or negative role. In this scenario, defining the various antitumor activities' molecular events is vital for selecting the appropriate HDACi therapy in solid tumors. There are also some limita-tions in this work. This article was primarily an in silico one. Limited wet lab studies were performed to have the first line of validation of the bioinformatics outputs, but larger sequencing groups of TNBCs need to be analyzed and detailed studies are required for further validation. In the future, we hope to conduct high-throughput sequencing and explore the potential mechanism between HDACis and elevated cytokine levels. In summary, the present findings revealed HDACi treatment-related RNA sequencing alterations in TNBC cells via bioinformatics analyses from the GEO database and in vitro experiments. Furthermore, the interaction network of differentially expressed genes revealed that changes in the cytokine-cytokine receptor pathway and CXCL1 secretion might have potential relevance for pharmacogenetic resistance to histone deacetylase inhibitors.

Conclusions
We downloaded the GSE129944 microarray dataset from the Gene Expression Omnibus database to identify changes at the transcriptome level in MDA-MB-231 breast cancer cells after treatment with belinostat. Unexpectedly, cells that underwent belinostat treatment gained a drug sensitivity advantage to epigenetic medication via CXCL1. That is, breast cancer patients with high expression of CXCL1 have a better prognosis, indicating that CXCL1 could be a novel favorable prognostic indicator. The induction of CXCL1 expression by belinostat treatment was also found in other TNBC cell lines. Detecting the emergence of clonal evolution under the pressure of tumor treatment selection will indicate prognosis and help us understand the disease's development.

Data Availability
The datasets used and analyzed during the present study are available from the GSE129944 (https://www.ncbi.nlm.nih .gov/geo/query/acc.cgi?acc=GSE129944).

Conflicts of Interest
We declare that we have no financial or personal relationships with other people or organizations that can inappropriately influence our work. There are no professional or other personal interests of any nature or kind in any product, service, and/or company that could be construed as influencing the position presented in this manuscript.

Authors' Contributions
Xin-le Han and Jun Du contributed equally to this work.