Targeting Class I Histone Deacetylases in Human Uterine Leiomyosarcoma

Uterine leiomyosarcoma (uLMS) is the most frequent subtype of uterine sarcoma that presents a poor prognosis, high rates of recurrence, and metastasis. Currently, the molecular mechanism of the origin and development of uLMS is unknown. Class I histone deacetylases (including HDAC1, 2, 3, and 8) are one of the major classes of the HDAC family and catalyze the removal of acetyl groups from lysine residues in histones and cellular proteins. Class I HDACs exhibit distinct cellular and subcellular expression patterns and are involved in many biological processes and diseases through diverse signaling pathways. However, the link between class I HDACs and uLMS is still being determined. In this study, we assessed the expression panel of Class I HDACs in uLMS and characterized the role and mechanism of class I HDACs in the pathogenesis of uLMS. Immunohistochemistry analysis revealed that HDAC1, 2, and 3 are aberrantly upregulated in uLMS tissues compared to adjacent myometrium. Immunoblot analysis demonstrated that the expression levels of HDAC 1, 2, and 3 exhibited a graded increase from normal and benign to malignant uterine tumor cells. Furthermore, inhibition of HDACs with Class I HDACs inhibitor (Tucidinostat) decreased the uLMS proliferation in a dose-dependent manner. Notably, gene set enrichment analysis of differentially expressed genes (DEGs) revealed that inhibition of HDACs with Tucidinostat altered several critical pathways. Moreover, multiple epigenetic analyses suggested that Tucidinostat may alter the transcriptome via reprogramming the oncogenic epigenome and inducing the changes in microRNA-target interaction in uLMS cells. In the parallel study, we also determined the effect of DL-sulforaphane on the uLMS. Our study demonstrated the relevance of class I HDACs proteins in the pathogenesis of malignant uLMS. Further understanding the role and mechanism of HDACs in uLMS may provide a promising and novel strategy for treating patients with this aggressive uterine cancer.


Introduction
Uterine leiomyosarcoma (uLMS) is a rare and aggressive uterine cancer, representing 1-2% of all uterine malignancies [1]. The annual incidence of uLMS is approximately 0.8 per 100,000 women [2]. The five years survival for all patients is between 25 and 76%, with survival for women with metastatic disease at the initial diagnosis approaching only 10-15% [3]. Although irrespective of treatment, the uLMS is characterized by poor prognosis [4], the present treatment for uLMS patients exhibits resistance to currently

Proliferation Assay
A trypan blue exclusion assay was performed for Cell proliferation measurement. Cells were seeded into 12-well tissue culture plates and treated with the Tucidinostat and DL-sulforaphane at a dose range of 1-25 µM for 48 hr. An equal amount of DMSO was used as vehicle control. After treatment, the cells were trypsinized and collected by centrifuge. The cells were resuspended in a serum-free medium. Equal volume of 0.4% trypan blue and cell suspension was mixed and applied to a hemacytometer for cell counting. Viable cells were unstained. This assay was performed three times in triplicate.

Protein Extraction and Western Blot
Protein extraction and specific protein bands visualization were performed as described previously [27]. The information about usage of primary antibodies is listed in Table 1. The antigen-antibody complex was detected with Trident Femto Western HRP substrate (GeneTex, Irvine, CA, USA). Specific protein bands were visualized using ChemiDoc maging system (Bio-Rad, Hercules, CA, USA).

RNA-Sequencing
The uLMS cell line (SK-UT-1) was treated with 5 µM Tucidinostat or 5 µM DL-Sulforaphane for 48 hr. Cells were subjected to RNA isolated using Trizol. RNA and library quality and quantity were assessed as described previously [27]. An Illumina NovaSEQ6000 was used for library sequencing.

Transcriptome Profiles Analysis
We designed the bioinformatics analyses on the basis of the flow diagram, as shown in Figure 1.

Transcriptome Data Analysis
The classical alignment-based mapper STAR, version 2.6.1d (GitHub, Inc., San Francisco, CA, USA) (23) was used to map sequencing reads to a human reference transcriptome. The results of STAR mapping were quantified by Salmon, version 1.4.0. Then, Bioconductor (https://bioconductor.org/packages/release/bioc/html/tximport.html, accessed on 27 October 2021) was used to read Salmon outputs into the R environment. Downstream analyses were performed as described previously [27].

Differential Gene Expression Analysis
To identify the differentially expressed genes (DEGs) between treatment and control groups, three count-based algorithms were implemented in R packages DESeq2 [29], edgeR [30], and Limma + voom [31]. For each of these three methods, we used a cutoff −1.5 > fold-change > 1.5 and a p-value of 0.05. In addition, Benjamini and Hochberg's (BH) method was performed to control the false discovery rate of all the genes with adjusted P-value less than 0.05.

Gene List Enrichment Analysis
Comprehensive gene set enrichment analysis for regulation machinery was carried out using the enrichR (version 3.1) [32] package in R (https://maayanlab.cloud/Enrichr/ (accessed on 27 October 2021). We used ENCODE Histone Modifications 2015 for histone modification enrichment, ENCODE and ChEA Consensus TFs for transcription factor enrichment and TargetScan microRNA for microRNAs enrichment in EnrichR to determine the mechanisms underlying the regulation of DEGs.

Drug Similarity Analysis
The L1000CDS2 is a pharmacogenetic search engine which enables users to find consensus L1000 small molecule signatures that match user input signatures. On the other hand, L1000CDS2 provides prioritization of thousands of small-molecule signatures, and

Transcriptome Data Analysis
The classical alignment-based mapper STAR, version 2.6.1d (GitHub, Inc., San Francisco, CA, USA) (23) was used to map sequencing reads to a human reference transcriptome. The results of STAR mapping were quantified by Salmon, version 1.4.0. Then, Bioconductor (https://bioconductor.org/packages/release/bioc/html/tximport.html, accessed on 27 October 2021) was used to read Salmon outputs into the R environment. Downstream analyses were performed as described previously [27].

Differential Gene Expression Analysis
To identify the differentially expressed genes (DEGs) between treatment and control groups, three count-based algorithms were implemented in R packages DESeq2 [29], edgeR [30], and Limma + voom [31]. For each of these three methods, we used a cutoff −1.5 > fold-change > 1.5 and a p-value of 0.05. In addition, Benjamini and Hochberg's (BH) method was performed to control the false discovery rate of all the genes with adjusted P-value less than 0.05.

Gene List Enrichment Analysis
Comprehensive gene set enrichment analysis for regulation machinery was carried out using the enrichR (version 3.1) [32] package in R (https://maayanlab.cloud/Enrichr/ (accessed on 27 October 2021). We used ENCODE Histone Modifications 2015 for histone modification enrichment, ENCODE and ChEA Consensus TFs for transcription factor enrichment and TargetScan microRNA for microRNAs enrichment in EnrichR to determine the mechanisms underlying the regulation of DEGs.

Drug Similarity Analysis
The L1000CDS2 is a pharmacogenetic search engine which enables users to find consensus L1000 small molecule signatures that match user input signatures. On the other hand, L1000CDS2 provides prioritization of thousands of small-molecule signatures, and their pairwise combinations. We performed drug similarity analysis using L1000CDS2 to investigate whether Tucidinostat and DL-sulforaphane induce the similar transcriptional changes with well-known HDAC inhibitors or not.

Visualization of Aggregates and Intersection DEGs on UpSet Plot
The UpSet plot was constructed to explore the interactive gene sets among different drug-induced expression profiles using R package UpSetR (version 1.4.0) [33]. For upregulated and down-regulated genes in all treatment, an eigengenes was calculated using the function module Eigengenes from WGCNA R package. Eigengenes is the vector that best describes the expression behavior of all genes within the module in the samples included in the analysis. To deep dive into pathways and relevant mechanisms of drugs, we integrated the present drugs (Tucidinostat and DL-sulforaphane) with our previous drug (TP472) to compare drug-induced expression profile.

Epithelial-Mesenchymal Transition Score Calculation
We used 76GS [34], KS [35], and ssGSEA [36] methods to calculate the Epithelialmesenchymal transition (EMT) scores for samples in the data set to investigate whether the drugs induced EMT or not.

Co-Expression Network Analysis
The WGCNA R package was applied to construct co-expression network [37]. The top 10% most variable genes were selected for co-expression analysis. Cluster tree sampling was calculated using the flashCust function in R to find out and exclude outlier samples in which Z.K values were under −2.5. Scale-free topology criterion was used to choose the power value in which degree of independence fits 0.85. WGCNA has detected modules with high correlation, which a minimum number of genes in each module determined as 30, a cut height of 0.25, and a deep split level of 2.
Module membership (MM) and module-clinical trait relationship was calculated using a correlation between module eigengene (ME) -the best summary of module expression based on the first principal component-and the special phenotype (EMT score, etc.) that can lead to discover key biological functions (BP) and recognize associated key biomarkers. MM calculation was used to choose the module for further analysis.
The EnrichR web server (version 3.1) was used to perform pathway enrichment analyses (https://maayanlab.cloud/Enrichr/ (accessed on 27 October 2021)). First, STRING database (https://string-db.org/ (accessed on 27 October 2021)), which is the online search tool to protein-protein network (PPI) construction, was used to reconstruct the modules network (A combined score ≥ 0.4 of PPI pairs was considered significant), then the Cytoscape software [38] was employed to analyze and visualize the network. The Drug Gene Interaction Database (DGIdb, v4.2.0) was used to identify potential drugs for EMT inhibition [39]. This database contains drug-gene interaction information from 22 source databases. To identify drug-gene interactions, only the approved interactions were considered.

Module-Based Drug Prediction
As an alternative approach, we used L1000CDS2 to identify candidate drugs for inhibition of EMT module. With the L1000CDS2 tool, we prioritized small molecules that can potentially reverse gene expression [40]. In this study, the L1000CDS2 were applied to prioritize small molecules that are predicted to reverse the expression profile of the EMT module.

Statistical Analysis
A comparison of the two and multiple groups were carried out as described previously [27]. Data were presented as mean ± standard error (SE), and the significant difference was defined as p < 0.05.

The Expression Levels of Class I HDACs Members Are Upregulated in uLMS Tissues
Compared to Adjacent Myometrium from Women with uLMS To determine the differential expression levels of Class I HDAC proteins between uLMS (n = 9) and MM +LMS (n = 7), IHC staining for HDAC1, 2, and 3 was performed. We examined the IHC images of three HDAC proteins in uLMS tissues vs. myometrium (MM). Figures 2 and 3 showed that the HDAC1 and HDAC2 positive cells were significantly higher in uLMS compared to MM +LMS . The H-score of HDAC1 and HDAC2 was also significantly increased in uLMS (n = 9) compared to MM (n = 7). Although HDAC3-positive cells showed no significant difference between uLMS vs. MM +LMS , the H-score of HDAC3 was significantly increased in uLMS compared to MM (Figures 2 and 3), indicating the critical role of class I HDACs in the pathogenesis and progression of uLMS. Figure 2 (right column) revealed an increase in expression density of HDAC1, 2, and 3 in uLMS compared to MM +LMS .

Module-Based Drug Prediction
As an alternative approach, we used L1000CDS2 to identify candidate drugs for inhibition of EMT module. With the L1000CDS2 tool, we prioritized small molecules that can potentially reverse gene expression [40]. In this study, the L1000CDS2 were applied to prioritize small molecules that are predicted to reverse the expression profile of the EMT module.

Statistical Analysis
A comparison of the two and multiple groups were carried out as described previously [27]. Data were presented as mean ± standard error (SE), and the significant difference was defined as p < 0.05.

The Expression Levels of Class I HDACs Members Are Upregulated in uLMS Tissues Compared to Adjacent Myometrium from Women with uLMS
To determine the differential expression levels of Class I HDAC proteins between uLMS (n = 9) and MM +LMS (n = 7), IHC staining for HDAC1, 2, and 3 was performed. We examined the IHC images of three HDAC proteins in uLMS tissues vs. myometrium (MM). Figures 2 and 3 showed that the HDAC1 and HDAC2 positive cells were significantly higher in uLMS compared to MM +LMS . The H-score of HDAC1 and HDAC2 was also significantly increased in uLMS (n = 9) compared to MM (n = 7). Although HDAC3-positive cells showed no significant difference between uLMS vs. MM +LMS , the H-score of HDAC3 was significantly increased in uLMS compared to MM (  IHC staining for HDAC1, 2, and 3 is presented with three representative cases. The right column Figure 2. IHC staining of HDAC1, 2, and 3 in human uLMS tissues and adjacent myometrium. IHC staining for HDAC1, 2, and 3 is presented with three representative cases. The right column showed the density map of HDAC1, 2, and 3 for the same representative case. Blue color: negative; Yellow color: low expression; brown color: moderate expression; red color: strong expression. Scale bars in black color: 100 µm; Scale bars in red color: 1mm. showed the density map of HDAC1, 2, and 3 for the same representative case. Blue color: negative; Yellow color: low expression; brown color: moderate expression; red color: strong expression. Scale bars in black color: 100 µm; Scale bars in red color: 1mm

Class I HDAC Protein Levels Are Upregulated in uLMS Cell Lines
The constitutive (basal) expression levels of Class I HDAC components in UTSM, HuLM, and uLMS cell lines were evaluated by immunoblot analysis. We demonstrated that the protein levels of HDAC1, 2, and 3 exhibited a graded increase from normal and benign UF tumor cells to malignant uLMS cells (p < 0.05) ( Figure 4A-C). In addition, although the protein levels of HDAC8 were not increased in uLMS compared to HuLM cells, a significant upregulation of HDAC8 was observed in two uLMS cell lines compared to the UTSM cell line (p < 0.05) ( Figure 4D).

Class I HDAC Protein Levels Are Upregulated in uLMS Cell Lines
The constitutive (basal) expression levels of Class I HDAC components in UTSM, HuLM, and uLMS cell lines were evaluated by immunoblot analysis. We demonstrated that the protein levels of HDAC1, 2, and 3 exhibited a graded increase from normal and benign UF tumor cells to malignant uLMS cells (p < 0.05) ( Figure 4A-C). In addition, although the protein levels of HDAC8 were not increased in uLMS compared to HuLM cells, a significant upregulation of HDAC8 was observed in two uLMS cell lines compared to the UTSM cell line (p < 0.05) ( Figure 4D).

Inhibition of HDACs Decreased the Cell Proliferation in uLMS Cells
Abnormal cell proliferation is common in many types of cancer. We detected the proliferation in UTSM, HuLM, and SK-UT-1 cell lines by Western blot using the antibody against PCNA, the cell proliferation marker. As shown in Figure 5A, the levels of PCNA in SK-UT-1 were highest among the three cell lines, suggesting that the uLMS cell line SK-UT-1 grew fastest. In addition, the PCNA levels in HuLM cell line are higher than in UTSM cell lines. In addition, we detected the Ki67-positive cells in uLMS and adjacent myometrium tissues. As shown in Figure 5B, a significant increase in Ki67-positive cells was observed in uLMS tissues compared to myometrium tissues.
Tucidinostat (chidamide), as a Class I HDAC inhibitor, has been shown to inhibit a variety of cancer growth [41][42][43][44][45][46]. Therefore, we selected Tucidinostat in our in vitro cell model to assess its effect on uLMS cell growth. The trypan blue exclusion assay was performed in SK-UT-1 and MES-SA cell lines treated with dose ranges from 1-25 µM.

Inhibition of HDACs Decreased the Cell Proliferation in uLMS Cells
Abnormal cell proliferation is common in many types of cancer. We detected the proliferation in UTSM, HuLM, and SK-UT-1 cell lines by Western blot using the antibody against PCNA, the cell proliferation marker. As shown in Figure 5A, the levels of PCNA in SK-UT-1 were highest among the three cell lines, suggesting that the uLMS cell line SK-UT-1 grew fastest. In addition, the PCNA levels in HuLM cell line are higher than in UTSM cell lines. In addition, we detected the Ki67-positive cells in uLMS and adjacent myometrium tissues. As shown in Figure 5B, a significant increase in Ki67-positive cells was observed in uLMS tissues compared to myometrium tissues.
Sulforaphane is an isothiocyanate present naturally in widely consumed vegetables and has been shown to have an inhibitory effect on various cancers [47][48][49][50][51][52][53]. In addition, DL-sulforaphane occurs naturally as L-isomer in edible cruciferous vegetables such as broccoli [54,55]. This study tested the effect of DL-sulforaphane on uLMS cell growth. Similarly, treatment of uLMS cell lines with DL-sulforaphane exhibited a dose-dependent cell growth inhibition for 48 h ( Figure 5E,F). Therefore, both Tucidinostat and DLsulforaphane elicit the anti-proliferation effect on uLMS cells.
Sulforaphane is an isothiocyanate present naturally in widely consumed vegetables and has been shown to have an inhibitory effect on various cancers [47][48][49][50][51][52][53]. In addition, DL-sulforaphane occurs naturally as L-isomer in edible cruciferous vegetables such as broccoli [54,55]. This study tested the effect of DL-sulforaphane on uLMS cell growth. Similarly, treatment of uLMS cell lines with DL-sulforaphane exhibited a dose-dependent cell growth inhibition for 48 h ( Figure 5E,F). Therefore, both Tucidinostat and DL-sulforaphane elicit the anti-proliferation effect on uLMS cells.

The Expression of Cell Cycle-and Apoptosis-Related Genes Is Altered upon Tucidinostat and DL-Sulforaphane Treatment
To determine genes regulating the cell cycle and apoptosis in uLMS cells in response to Tucidinostat and DL-sulforaphane treatments, we checked the expression of CDKN1A (P21), BAK1, CDK1, CDK3, CDK10, and HDAC6. As shown in Figure 9, both Tucidinostat and DL-sulforaphane increased the expression of CDKN1A and BAK1 and reduced the expression of CDK3 and CDK10. In addition, Tucidinostat but not DL-sulforaphane reduced the expression of CDK1. Moreover, both Tucidinostat and DL-sulforaphane decreased the expression of HDAC6 in uLMS cells. BAK localizes to mitochondria, and functions to induce apoptosis. HDAC6 is involved in cell cycle and apoptosis [56], and CDKN1A and CDK members were critical in cell cycle progression [57][58][59]. Therefore, our results suggested that Tucidinostat and DL-sulforaphane suppressed the uLMS proliferation via cell cycle arrest and apoptosis.
HDACs have been implicated as a target for sulforaphane. A number of studies have demonstrated that sulforaphane impacts the cellular acetylome and decreases HDAC activity in several models and diseases, including prostate epithelial cells, mouse colon tissues, satellite cells, and human peripheral blood mononuclear cells [60][61][62],

The Expression of Cell Cycle-and Apoptosis-Related Genes Is Altered upon Tucidinostat and DL-Sulforaphane Treatment
To determine genes regulating the cell cycle and apoptosis in uLMS cells in response to Tucidinostat and DL-sulforaphane treatments, we checked the expression of CDKN1A (P21), BAK1, CDK1, CDK3, CDK10, and HDAC6. As shown in Figure 9, both Tucidinostat and DL-sulforaphane increased the expression of CDKN1A and BAK1 and reduced the expression of CDK3 and CDK10. In addition, Tucidinostat but not DL-sulforaphane reduced the expression of CDK1. Moreover, both Tucidinostat and DL-sulforaphane decreased the expression of HDAC6 in uLMS cells. BAK localizes to mitochondria, and functions to induce apoptosis. HDAC6 is involved in cell cycle and apoptosis [56], and CDKN1A and CDK members were critical in cell cycle progression [57][58][59]. Therefore, our results suggested that Tucidinostat and DL-sulforaphane suppressed the uLMS proliferation via cell cycle arrest and apoptosis.
using L1000CDS2 and demonstrated that Tucidinostat-induced transcriptional pattern has similarity with signatures of multiple HDAC inhibitors, including BRD-K11663430, HDAC6 inhibitor ISOX, mocetinostat, vorinostat, entinostat, among others, which were shown in Table S1. However, the transcriptional signature induced by DL-sulforaphane did not show similarity with any other HDACi-induced transcriptional pattern (Table S2), indicating that DL-sulforaphane may not target HDAC directly in uLMS cells.

Tucidinostat and DL-Sulforaphane Altered the Gene Expression Correlating to Histone Modifications
To understand epigenetic-mediated transcriptional changes in response to the Tucidinostat and DL-sulforaphane treatment, we performed enrichment analysis of epigenetic histone markers using the Enrichr web server. As shown in Figure 10, several histone modifications, including H3K27me3 and H3K9me3, associated with upregulated DEGs in response to Tucidinostat were identified ( Figure 10A). In addition, downregulated DEGs associated with histone modifications, including H3K27ac, H3K4me3, HDACs have been implicated as a target for sulforaphane. A number of studies have demonstrated that sulforaphane impacts the cellular acetylome and decreases HDAC activity in several models and diseases, including prostate epithelial cells, mouse colon tissues, satellite cells, and human peripheral blood mononuclear cells [60][61][62], concomitantly with an increase in the levels of acetylated histones and p21. In addition, separate studies demonstrated that other types of cells, such as HCT116 is at least partially resistant to the nuclear HDAC inhibitory effect of Sulforaphane [63]. Therefore, sulforaphane exerted a distinct action on transcriptome changes in a cell-dependent manner.
To assess if Tucidinostat and DL-sulforaphane exhibited similar transcriptional pattern in uLMS with other HDAC inhibitors, we performed drug similarity analysis using L1000CDS2 and demonstrated that Tucidinostat-induced transcriptional pattern has similarity with signatures of multiple HDAC inhibitors, including BRD-K11663430, HDAC6 inhibitor ISOX, mocetinostat, vorinostat, entinostat, among others, which were shown in Table S1. However, the transcriptional signature induced by DL-sulforaphane did not show similarity with any other HDACi-induced transcriptional pattern (Table S2), indicating that DL-sulforaphane may not target HDAC directly in uLMS cells.

Tucidinostat and DL-Sulforaphane Altered the Gene Expression Associated with Transcriptional Factors
Transcriptional factors play an important in many biological processes and their control is disrupted in cancer cells [64]. The dysregulation of these core transcription factors forms interconnected transcriptional loops to establish and reinforce the abnormal gene-expression program in cancer cells [65]. Our studies demonstrated that Tucidinostat and DL-sulforaphane induced-DEGs are putative targets of multiple transcriptional factors, as shown in Figure 11A-D.

Tucidinostat and DL-Sulforaphane Altered the Gene Expression Associated with Transcriptional Factors
Transcriptional factors play an important in many biological processes and their control is disrupted in cancer cells [64]. The dysregulation of these core transcription factors forms interconnected transcriptional loops to establish and reinforce the abnormal gene-expression program in cancer cells [65]. Our studies demonstrated that Tucidinostat and DL-sulforaphane induced-DEGs are putative targets of multiple transcriptional factors, as shown in Figure 11A-D.  Cells 2022, 11, x FOR PEER REVIEW 15 of 28

Tucidinostat and DL-Sulforaphane Altered the Gene Expression Correlating to the miRNA Regulation
We used TargetScan microRNA analysis in Enrichr web server to determine the mechanism underlying the regulation of DEGs associated with miRNAs in response to HDAC inhibitors treatment. As shown in Figure 12A-D, Tucidinostat and DLsulforaphane induced-DEGs are putative targets of multiple miRNAs.

Tucidinostat and DL-Sulforaphane Altered the Gene Expression Correlating to the miRNA Regulation
We used TargetScan microRNA analysis in Enrichr web server to determine the mechanism underlying the regulation of DEGs associated with miRNAs in response to HDAC inhibitors treatment. As shown in Figure 12A-D, Tucidinostat and DL-sulforaphane induced-DEGs are putative targets of multiple miRNAs.

UpSet Plot Visualization
To further visualize the intersections of three drugs-induced DEGs (Tucidinostat up, Tucidinostat down, DL-sulforaphane up, DL-sulforaphane down, TP472 up, TP472 down), the upset plot was utilized to present the distribution characteristics of DEGs upon treatments. As shown in Figure 13A, the common up and down DEGs with Tucidinostat and DLsulforaphane contain 71 and 140 genes, respectively. The TP472-down/DL-sulforaphanedown groups had 513 genes, the group with the largest number of genes in all groups with genes involved in two types of gene regulations. In addition, the Tucidinostat-down/DLsulforaphane-down/TP472-down groups contained 358 genes, the group with the largest number of genes in all groups involved in three types of regulations. Figure 13B,C exhibited the expression pattern in the red intersection (Tucidinostat-down, DL-sulforaphane-down, and TP472-down), the blue intersection (Tucidinostat-up, DL-sulforaphane-up, and TP472up) presented in Figure 13A.

UpSet Plot Visualization
To further visualize the intersections of three drugs-induced DEGs (Tucidinostat up, Tucidinostat down, DL-sulforaphane up, DL-sulforaphane down, TP472 up, TP472 down), the upset plot was utilized to present the distribution characteristics of DEGs upon treatments. As shown in Figure 13A, the common up and down DEGs with Tucidinostat and DL-sulforaphane contain 71 and 140 genes, respectively. The TP472-down/DLsulforaphane-down groups had 513 genes, the group with the largest number of genes in all groups with genes involved in two types of gene regulations. In addition, the Tucidinostat-down/DL-sulforaphane-down/TP472-down groups contained 358 genes, the group with the largest number of genes in all groups involved in three types of regulations. Figure 13B,C exhibited the expression pattern in the red intersection (Tucidinostat-down, DL-sulforaphane-down, and TP472-down), the blue intersection (Tucidinostat-up, DL-sulforaphane-up, and TP472-up) presented in Figure 13A.
Functional enrichment analysis revealed that the genes from the blue and red intersections in Figure 13A exhibited multiple enriched pathways, including inflammatory response, p53 pathway, KRAS signaling, apoptosis, estrogen response, UV response, complement, angiogenesis, IL-2/STAT5 signaling, EMT, among others ( Figures  13D,E). Functional enrichment analysis revealed that the genes from the blue and red intersections in Figure 13A exhibited multiple enriched pathways, including inflammatory response, p53 pathway, KRAS signaling, apoptosis, estrogen response, UV response, complement, angiogenesis, IL-2/STAT5 signaling, EMT, among others ( Figure 13D,E).
Since the number of shared DEGs between Tucidinostat and DL-sulforaphane are much less compared to the number of DEGs from each individual group, one may consider that Tucidinostat and DL-sulforaphane promoted the inhibitory effect on uLMS via different network mechanisms.

Apoptosis and EMT
To further determine the apoptosis pathway affected by the treatments, we used the single-sample gene set enrichment analysis (ssGSEA) to measure the apoptosis scoring in response to the Tucidinostat, DL-sulforaphane, and TP472. We demonstrated that treatments with above three agents significantly increased the apoptosis scoring compared to the control (DMSO) group. Moreover, Tucidinostat exhibited the most effective of activating apoptosis among the three drugs ( Figure 14A).
Since the number of shared DEGs between Tucidinostat and DL-sulforaphane are much less compared to the number of DEGs from each individual group, one may consider that Tucidinostat and DL-sulforaphane promoted the inhibitory effect on uLMS via different network mechanisms.  For epithelial-mesenchymal transition (EMT) scoring, we used three different methods to quantify EMT; 76GS, KS, and ssGSEA. The KS method has a predefined scale for EMT scores [−1, 1], with higher scores indicating mesenchymal signatures. On the other hand, there is no pre-defined range of scores calculated by 76GS, and the higher the 76GS score, the more epithelial the sample is. The score derived from ssGSEA reflects the degree to which the input gene signature is coordinately up or downregulated within samples ( Figure 14A,B).

Co-Expression Network Analysis
The weighted gene co-expression network was constructed with 1356 most variable genes (the top 10% of genes). (Figure 15A) The scale-free topology R2 did not reach the soft threshold of 0.85, so the recommended power value of 12 was chosen ( Figure 14B). The WGCNA revealed two gene co-expression modules by average linkage hierarchical clustering. Modules with more than 0.25 expression profiles similarity were merged. Each module was shown in a unique color ( Figure 15C). The Module-Blue was found to have the highest significant association with EMT (correlation = 0.99 and p-value < 9 × 10 −13 ). (Figure 15D,E). The 255 genes were found in the Module-Blue, which used the STRING tool to reconstruct the gene-gene interactions (https://string-db.org/ (accessed on 27 October 2021). Cytoscape software identified the hub genes associated with the Module-Blue. Based on intramodular connectivity, top hub genes in the co-expression network are shown in Figure 16.

Potential Drugs Prediction
This study used DGIdb as the first approach for identifying the possible drugs with the effects on reversing the increased expression of EMT hub genes. Using DGIdb, we found 186 candidate drugs that target the top 10 percent of EMT modules' genes (based on the network's connectivity). These potential drugs are shown in Table S3.
Our research used L1000CDS2 as the second approach to identify drug candidates for EMT inhibition. The list of the EMT modules DEGs and their related fold changes between the Tucidinostat group (as the most significant treatments group) vs. Control was entered to the L1000CDS2 web tool to search for molecular compounds that can reverse the expression changes. The top 50 drugs are identified as an output and are shown in Table S4. This study proposes future research to use a combination treatment strategy (HDACi in combination with EMT inhibitors) to achieve a better outcome in treating human uLMS.
( Figure 15D,E). The 255 genes were found in the Module-Blue, which used the STRING tool to reconstruct the gene-gene interactions (https://string-db.org/ (accessed on 27 October 2021). Cytoscape software identified the hub genes associated with the Module-Blue. Based on intramodular connectivity, top hub genes in the co-expression network are shown in Figure 16.   ( Figure 15D,E). The 255 genes were found in the Module-Blue, which used the STRING tool to reconstruct the gene-gene interactions (https://string-db.org/ (accessed on 27 October 2021). Cytoscape software identified the hub genes associated with the Module-Blue. Based on intramodular connectivity, top hub genes in the co-expression network are shown in Figure 16.

Discussion
ULMS is a highly aggressive tumor with high tumor recurrence rates, progression, and metastasis [5]. The malignant potential of uterine fibroids is extremely low [66,67]. The origin and molecular mechanism underlying driving its clinical and biological behavior remain unclear [6].
This study demonstrated that Class I HDACs are abnormally upregulated in uLMS compared to myometrial tissues, which may contribute to the uLMS pathogenesis.
We compared the expression levels of main class I HDAC members HDAC1, 2, and 3 in uLMS and adjacent myometrium, as well as the expression of full members of Class I HDACs in UTSM, HuLM, and uLMS cell lines. The protein levels of HDAC1, 2, and 3 are significantly upregulated in two uLMS cell lines compared to UTSM and HuLM cell lines. In addition, the expression levels of HDAC 1,2 and 3 are upregulated in uLMS tumors compared to myometrium tissues suggesting that Class I HDACs may contribute to the pathogenesis of uLMS.
Abnormal cell proliferation via decreasing the cell cycle arrest and apoptosis is common in many cancers [68][69][70][71]. We revealed that uLMS cells grow faster than myometrial cells in vitro and in vivo. To determine if targeting HDACs can suppress the uLMS phenotype, we determined the uLMS cell growth in response to HDAC inhibitor treatment. Our study demonstrated that HDAC inhibition suppressed the uLMS cell proliferation. The anti-tumor effect of HDAC inhibition observed in our model was consistent with literature from other types of cancer. Targeting class I HDACs with Tucidinostat (benzamide HDAC inhibitor) has been studied in different types of cancers showing beneficial effects. For instance, in pancreatic cancer, Tucidinostat treatment synergistically enhances gemcitabine cytotoxicity in pancreatic cancer cells [72]. In multiple myeloma (NN), Tucidinostat inhibited the proliferation and invasion of MM cells. In addition, Tucidinostat in combination with lenalidomide and a low dose of borezomid exhibited a synergistic effect in MM [73]. In hepatic cancer, Tucidinostat showed an anti-tumor activity [74]. Moreover, Tucidinostat targeted stem and progenitor cells of acute myeloid leukemia [75]. All those studies demonstrated the critical role of Class I HDACs in cancer development and targeting class I HDACs showed beneficial effects in several types of neoplasms.
To further determine the mechanisms associated with inhibition, we performed a genome-wide RNA-sequencing experiment comparing the DMSO-treated with Tucidinostat -treated uLMS cells. In addition, DL-sulforaphane-induced transcriptional changes were also compared. The transcriptome analysis revealed that targeting HDACs with Tucidinostat altered several critical biological pathways that may contribute to uLMS pathogenesis. In addition, we determined the transcriptional changes in response to DL-sulforaphane treatment. Although DL-sulforaphane showed the activity of inhibiting cell proliferation, the specificity of the drug targeting class I HDACs has not been reported in uLMS. Our analysis demonstrated that 21.5% and 27.5% of down DEGs shared common genes between Tucidinostat and DL-sulforaphane treatment groups. Similarly, 18.1% and 40.8% of up DEGs showed common genes between Tucidinostat and DL-sulforaphane treatment groups. This overlapped analysis is consistent with the hallmark analysis demonstrating that Tucidinostat and DL-sulforaphane impaired common pathways. However, the alteration of some other biological pathways differed in response to Tucidinostat and DL-Sulforaphane, respectively.
In contrast to Tucidinostat, the inhibitory effect of Sulforaphane on HDAC activity is controversial. For example, sulforaphane treatment does not reduce nuclear HDAC activity, but decreases the levels of HDAC1-4 and 6 in keratinocytes [63]. However, other studies demonstrated the inhibitory effect of sulforaphane on HDAC activity. For example, sulforaphane repressed the HDAC activity by 40%, 30%, and 40% in BPH-1, LnCaP, and PC-3 prostate epithelial cells, respectively. The HDAC inhibition was accompanied by a 50-100% increase in acetylated histones in all three prostate cell lines [61]. A similar study showed that sulforaphane dramatically reduced HDAC activity in porcine satellite cells with an increase in global acetylated H3 and H4 levels [62]. The in vivo studies demonstrated that sulforaphane could decrease HDAC activity by~65%, concomitantly with an increase in acetylated histones globally, as well as locally on the promoters of genes such as P21 and BAX [60]. Therefore, sulforaphane may have a distinct impact on different cell types. In this study, our drug similarity analysis demonstrated that Tucidinostat-induced transcriptional signature had a similar pattern with several other known HDAC inhibitors-induced patterns. Still, DL-sulforaphane did not, indicating that DL-sulforaphane may affect uLMS via different mechanisms, rather than HDAC inhibition in uLMS. Notably, we demonstrated that the more inhibitory efficiency of Tucidinostat over DL-sulforaphane is consistent with the expression levels of cell cycle-related genes, including P21, CDK1, and CDK3, as shown in Figure 9. For the apoptosis-related gene, both Tucidinostat and DL-sulforaphane increased the expression levels of BAK1, which localizes to mitochondria, and functions to induce apoptosis. Moreover, Tucidinostat probably promoted apoptosis more strongly than DL-sulforaphane identified by GSEA Score. Therefore, targeting class I HDACs with Tucidinostat proves superior to DL-sulforaphane by inhibiting the uLMS proliferation via inducing apoptosis and cell cycle arrest in uLMS cells.
EMT is a cell biological process crucial for tumor aggressiveness, including cancer metastasis and drug response [76,77]. Therefore, the EMT status of cancer cells can be proved to be a critical estimate of patient prognosis. In this regard, we used three distinct metrics that score EMT on a continuum, based on the transcription signature of Tucidinostat and DL-sulforaphane and control groups. However, our results demonstrated that both drugs increase the EMT levels compared to the DMSO control ( Figure 13B). It has been reported that dysregulation of apoptosis and EMT are linked with various pathological progress, including tumor formation and progression [78]. Notably, TGF-β, as a potent pleiotropic molecule, induced apoptosis and simultaneously induced the EMT of AML-12 cells. The question is if targeting the HDACs induces these concurrent but distinct events in uLMS cells. A study by Yang et al. demonstrated that TGF-β1-induced apoptosis and EMT were closely related to the cell cycle stage, and TGF β 1-induced concurrent apoptosis and EMT are independent of each other [79]. Therefore, deep diving into the mechanism underlying the Tucidinostat and DL-sulforaphane-induced changes in apoptosis, EMT, and other pathways is worthwhile.
It has been reported that cross-talk between different epigenetic mechanisms regulates gene expression [80][81][82][83]. For instance, HDAC inhibitors can elicit transgenerational effects via altered DNA and histone methylation [84]. HDAC4 and HDAC9 can differentially influence H3K27 acetylation, which can explain the pleiotropic actions of MEF2 transcription factors in uLMS [85]. In this study, Tucidinostat altered the up-DEGs correlated with multiple histone modifications in uLMS, including H3K27me3, H3K79me2, 3, and H3K9ac, among others. The down-DEGs correlated with H3K27ac, H3K4me1, 2, 3, and H3K27me3, among others. Our studies indicated that HDAC inhibition might alter the transcription by reprogramming the oncogenic transcriptome in uLMS. This observation was further confirmed by transcriptional factor regulation analysis inferred from integrating genome-wide ChIP-X (ChIP-chip, ChIP-seq, ChIP-PET, and DamID) studies. By combining transcriptome data with ChIP-X transcriptional factor, and histone modification studies, we identified multiple putative networks linking to enriched pathways, that can help target specific transcription factor activity in uLMS cells using combination drug treatment strategy.
Notably, the interplay between miRNAs and HDACs has been widely reported [86][87][88]. We used Targetscan miRNA as a prediction tool to assess which miRNAs may interact with their targets. A number of miRNAs are shown to putatively interact with Tucidinostatinduced DEGs. In addition, DL-sulforaphane treatment also altered the interaction between miRNAs and target mRNA. Subsequent experiments for functional validation of miRNAtarget interactions in response to Tucidinostat and Dl-sulforaphane treatments are needed.
HDACs are divided into five groups based on sequence homology to the original yeast enzymes and domain organization. Among them, HDAC1, 2, 3, and 8 belong to Class I. Class IIa contains HDAC4, 5, 7, and 9. A previous study demonstrated that Tucidinostat targets Class I HDACs and sulforaphane attenuates class ILa HDACs and HDAC2 enzyme activities [89]. This study revealed several differences between Tucidinostat and DL-sulforaphane treatments in uLMS cells. (1) Tucidinostat showed a more potent inhibitory effect on cell proliferation than DL-sulforaphane. (2) Venn diagrams analysis demonstrated a distinct transcriptome pattern between Tucidinostat and DL-sulforaphane.
(3) The expression levels of key genes such as p21 differed between Tucidinostat and DL-sulforaphane treatments. In addition, the enriched pathways and reprogramming differed between the two treatments. Tucidinostat showed similarity with other HDAC inhibitors in the transcriptome pattern, but DL-sulforaphane did not. From upset plot analysis, the groups of Tucidinostat and DL-sulforaphane contained a smaller number of common genes than either Tucidinostat/TP472 groups or DL-sulforaphane/TP472 groups, respectively. Since TP472 mainly targets BRD9, one may consider that Tucidinostat and DL-sulforaphane trigger the distinct expression pattern in uLMS cells markedly.
Per our studies, we proposed a mechanism model for targeting HDACs in uLMS based on our novel findings that (1) Class I HDACs expression is dysregulated in uLMS tumors and cells, (2) targeting HDACs with Tucidinostat alters the uLMS phenotype with a decrease in cell proliferation and modulation of cell-cycle related genes, and increases the apoptosis process; (3) Tucidinostat and DL-sulforaphane reversed the phenotype of uLMS via different mechanism; (4) Class I HDACs constitutes a distinguished vulnerability in malignant uLMS, and HDAC inhibitors, such as Tucidinostat, alter key pathways and reprogram the oncogenic epigenome and miRNA network to suppress the uLMS phenotype ( Figure 17). The concept that preoperative therapy leads to the improvement of oncological results is widely approved for long-term survival in potentially curative cases and in those with metastatic diseases [90][91][92]. Unfortunately, no single preoperative test can reliably differentiate benign from malignant uterine disease [93]. Therefore, identifying biomarkers to differentiate malignant uLMS from benign UFs will help to initiate preoperative therapy to achieve a better outcome for patients with uLMS. The concept that preoperative therapy leads to the improvement of oncological results is widely approved for long-term survival in potentially curative cases and in those with metastatic diseases [90][91][92]. Unfortunately, no single preoperative test can reliably differentiate benign from malignant uterine disease [93]. Therefore, identifying biomarkers to differentiate malignant uLMS from benign UFs will help to initiate preoperative therapy to achieve a better outcome for patients with uLMS.

Conclusions
Our study demonstrated that uLMS tumors and cells exhibited an aberrant upregulation of class I HDAC proteins. Targeting HDACs in uLMS may impart beneficial effects in uLMS and provide a promising and novel strategy for treating patients with this aggressive uterine cancer. Furthermore, the present study provided a list of potential medications as inhibitors for EMT that can be used in combination with antineoplastic drugs to increase the sensitivity of tumor cells to improve responses to therapy in uLMS cells.