Targeting DNA Methylation Depletes Uterine Leiomyoma Stem Cell–enriched Population by Stimulating Their Differentiation

Uterine leiomyoma (LM) is the most common tumor in women and can cause severe morbidity. Leiomyoma growth requires the maintenance and proliferation of a stem cell population. Dysregulated deoxyribonucleic acid (DNA) methylation has been reported in LM, but its role in LM stem cell regulation remains unclear. Here, we fluorescence-activated cell sorting (FACS)-sorted cells from human LM tissues into 3 populations: LM stem cell–like cells (LSC, 5%), LM intermediate cells (LIC, 7%), and differentiated LM cells (LDC, 88%), and we analyzed the transcriptome and epigenetic landscape of LM cells at different differentiation stages. Leiomyoma stem cell–like cells harbored a unique methylome, with 8862 differentially methylated regions compared to LIC and 9444 compared to LDC, most of which were hypermethylated. Consistent with global hypermethylation, transcript levels of TET1 and TET3 methylcytosine dioxygenases were lower in LSC. Integrative analyses revealed an inverse relationship between methylation and gene expression changes during LSC differentiation. In LSC, hypermethylation suppressed the genes important for myometrium- and LM-associated functions, including muscle contraction and hormone action, to maintain stemness. The hypomethylating drug, 5 ′ -Aza, stimulated LSC differentiation, depleting the stem cell population and inhibiting tumor initiation. Our data suggest that DNA methylation maintains the pool of LSC, which is critical for the regeneration of LM tumors.

Uterine leiomyoma (LM, or fibroid) is the most common tumor in women, affecting approximately 80% of women worldwide (1,2). It affects African American women disproportionately, beginning at an earlier age and with a higher prevalence (1,3). Although benign in origin, LM can cause significant morbidity, including excessive uterine bleeding, recurrent pregnancy loss, and pelvic pain; these symptoms may also mimic or mask malignant tumors (2). More than 200 000 surgeries are performed in the Unites States each year to remove LM, at an estimated cost of $5.9-34.4 billion (4,5). Despite their prevalence and impact on women's health, no new medical treatments have been approved in the United States for LM since the 1990s (6).
Stem cells are critical for normal tissue and disease development. Several studies have demonstrated the presence of stem cells in LM, which are indispensable for proliferation, regeneration, and tumor growth (7)(8)(9). Leiomyoma stem cells are deficient in estrogen and progesterone receptors, which may explain the frequent regrowth of LM that is observed after the discontinuation of hormone-based therapies (10). Using (10,11). Microarray-based gene expression profiling demonstrated that the LSC population has a distinct transcriptome from LIC and LDC but how LSC maintain their unique gene expression profile and stemness remains unclear (11).
Epigenetic alterations, such as deoxyribonucleic acid (DNA) methylation, play crucial roles in stem cell regulation and disease progression (12,13). Within cancer stem cells, the genes required for differentiation are suppressed by DNA methylation, leading to abnormal clonal expansion (14)(15)(16)(17). In several cancer types, researchers have found that DNA-hypomethylating agents can "reset" cancer stem-like cells towards a differentiation phenotype, thereby resensitizing them to chemotherapy (17). Previously, it was demonstrated that, compared with myometrium, LM contain a dysregulated DNA methylation landscape, with hypermethylation identified around the promoter regions of tumor suppressors such as KLF11, DLEC1, and KRT19 (18,19). It remains unclear how DNA methylation affects LSC differentiation and tumorigenesis.
In this study, we characterized the global transcriptomes and methylomes of each LM cell population and evaluated the role of DNA methylation during the LSC differentiation process. We hypothesized that DNA methylation "locks" genes critical for differentiation in a suppressed state, which maintains LSC function, and how removing this epigenetic suppression would drive LSC differentiation, resulting in the depletion of the LSC population.

Tissue collection
Northwestern University's Institutional Review Board approved the use of human tissue. Leiomyoma tissues were obtained from premenopausal women undergoing either myomectomy or hysterectomy (age 38 ± 9 years, range 29-47 years). Informed consent was obtained from all participants prior to tissue collection. Patients receiving hormone treatment within 6 months prior to surgery were excluded. Information of race/ethnicity and MED12 mutation status for patient samples used for next generation sequencing (RNA-Seq and MethyCap-Seq) are listed in Table 1. Tissues were dissociated and cells were isolated, as previously described (10).

Antibody and primer information
All antibodies and primers used in this manuscript are listed in Table 2.
Ribonucleic acid isolation, real-time quantitative polymerase chain reaction, and RNA-Seq

MethylCap-Seq
Genomic DNA was extracted from sorted LM cells using the Qiagen Allprep RNA/DNA mini kit (80204, Qiagen, Hilden, Germany) and was fragmented to 300-500 bp using Covaris M220 (Covaris, Woburn, Massachusetts). Methylated DNA fragments were captured using the MethylCap Kit (C02020010, Diagenode, Sparta, New Jersey) following the manufacturer's protocol. Briefly, 500 ng of fragmented genomic DNA was incubated with H6-GST-MBD fusion proteins that bind methylated cytosines. Previous studies have shown that a single, fully methylated CpG is sufficient for the interaction between the H6-GST-MBD fusion protein and methylated DNA (22). The protein-DNA complex was then precipitated with antibody-conjugated beads specific to the protein tag.
The immunoprecipitated DNA fragments were purified and subjected to library construction and next-generation sequencing, as described below. Deep sequencing provides greater genome coverage, representing the majority of MBD-bound methylated DNA.

Native ChIP-Seq for histone modification
We followed the protocol published by Peter van Galen et al for native ChIP-Seq on histone modifications using low cell numbers (23). All sequencing barcodes used in this study were provided in the original publication (23).

Colony formation assay
Primary LM cells were treated with vehicle (DMSO; Sigma-Aldrich, St. Louis, Missouri), 5′-Aza, or RU486 for 6 days, as described under primary cell and tissue explant culture. After treatment, 500 viable cells were seeded in each well of 12-well plates and cultured for 21 days without any treatment. Plates were then washed with 2 mL PBS, fixed with 3 mL 10% formalin (252549, Sigma) for 15 minutes, and stained with crystal violet for 5 minutes (0.025% w/v; HT90132, Sigma-Aldrich, St. Louis, Missouri). The number of colonies was counted in each well, excluding small colonies (<50 cells) (32).

Animal studies
Northwestern University's Animal Care and Use Committee approved all procedures involving animals. For each xenograft, 10 6 live cells were grafted underneath the kidney capsules of ovariectomized nonobese diabetic-scid IL2Rgnull mouse hosts (NSG, Jackson Laboratory, Bar Harbor, Maine), as previously described (20). To test whether 5′-Aza affects the tumor-initiation capacity of primary LM cells, cells were treated with vehicle (DMSO), 5′-Aza (100 nM), or RU486 (10 -6 M) for 6 days prior to xenograft. Mice were euthanized 4 weeks postsurgery. Images of regenerated tumors on the kidney surface were taken from the x-, y-, and z-axes using a dissecting microscope connected to a computer with Leica Application Suite, version 3.8 software (Leica Microsystems Inc, Wetzlar, Germany). Tumors were measured by 2 individuals who were blinded to the treatment group allocation; the average measurement taken by the 2 individuals was used for data plotting. Tumor volume was quantified using the following formula: volume (mm 3 ) = 0.52 (derived from π/6) * length * width * height (mm 3 ) (33).

Statistical analysis
All statistical analyses were performed using GraphPad Prism 8 (GraphPad Inc., San Diego, California) and R (3.6.0); detailed statistical test descriptions are reported in each figure legend. No sample was excluded during the analysis. Values were considered statistically significant when P < 0.05. All experiments were repeated with samples from at least 3 patients, with the patient number (n) noted in the figure legends. Data points in the bar plots represent biological replicates from different patients or mice, and the error bars represent Standard error of the mean.

Results
The LSC population has a unique DNA methylation landscape For our unbiased genome-wide evaluation of the DNA methylome and transcriptome, we performed highthroughput sequencing analyses, including MethylCap-Seq and RNA-Seq, on the 3 LM cell populations freshly sorted from patient LM tissues (Fig. 1A). MethylCap-Seq characterized the global DNA methylation landscape of each LM cell population. Leiomyoma stem cell-like cells showed a unique DNA methylome compared to LIC and LDC (Fig. 1B, n = 3 patients). On average, LSC had the highest DNA methylation level, whereas LIC had the lowest DNA methylation level (Fig. 1C). In LSC, we also identified 7020 regions with higher DNA methylation and 1842 with lower DNA methylation compared with LIC, and 8364 regions with higher DNA methylation and 1080 with lower DNA methylation compared with LDC ( Fig. 1D). Leiomyoma intermediate cells and LDC had comparable DNA methylomes, with only 288 differentially methylated regions (DMRs) detected, indicating that DNA methylation plays an important role during the initial differentiation of LSC. We evaluated the genomic distribution of DMRs between LSC and the nonstem cell populations. Compared with randomly picked genomic regions (0.7% at CpG islands and 3.5% at shores), DMRs in LSC were highly enriched in CpG islands and shores (Fig. 1E). This feature was more significant in regions with lower DNA methylation in LSC: 19.9% CpG islands and 35.7% shores were found in the LIC and LSC comparison; and 29.3% CpG islands and 47.9% shores were found in the LDC and LSC comparison (Fig. 1E). These findings were consistent with the DNA methylation pattern changes observed during intestinal stem cell differentiation using whole-genome shotgun bisulfite sequencing (13).

Transcript levels of TET1 and TET3 methylcytosine dioxygenases are the lowest in LSC
We then performed RNA-seq to reveal possible underlying mechanisms for the distinct DNA methylation pattern in LSC. The principal component and hierarchical clustering analyses of RNA-Seq data demonstrated distinct clustering of the 12 samples (n = 4 patients) into the 3 cell populations and a matching pattern between the transcriptomes and methylomes ( Fig. 2A and 2B). Consistent with previous observations based on microarray data, heatmap clustering ( Fig. 2A) indicated a significant transcriptome difference between LSC and the nonstem cell populations (11). Compared to both LIC and LDC, differential gene expression analysis detected 2401 genes with higher expression and 1956 genes with lower expression in LSC (Fig. 2C). Pathway enrichment analysis demonstrated that genes highly expressed in LSC were enriched in pathways essential for stem cell maintenance and function, including immune response and G-protein signaling (34,35), whereas genes highly expressed in LIC or LDC were enriched in pathways related to the committed cell function, including smooth muscle contraction, ESR1 nuclear pathway, and TGF-β signaling, which are known to be important for LM tumorigenesis (33,36) (Fig. 2D).
Supporting the distinct methylation landscape identified by MethylCap-Seq, genes involved in DNA methylation and demethylation processes (Gene Ontology Term [GO:0044728]) maintained a unique expression pattern in LSC compared with the other 2 cell populations (Fig. 2E). Among the DNA methylation machineries, real-time qPCR analysis demonstrated that DNMT1 and MBD1 showed slightly higher expression in LSC, but these differences did not reach statistical significance (data not shown). Importantly, the mRNA levels of TET1 and TET3, 2 DNA demethylases, were expressed at the lowest levels in LSC (Fig. 2F), providing a possible molecular explanation for the global hypermethylation observed in the LSC population.

DNA methylation changes correlate with gene expression changes globally during LSC differentiation
To investigate the potential influence of DNA methylation on gene expression changes during LSC differentiation, we integrated the RNA-Seq and MethylCap-Seq data. We compared differentially methylated and expressed genes between LSC and LIC/LDC. As shown in Fig. 3A, among 8862 DMRs between LSC and LIC, 3230 regions were associated with 1772 genes differentially expressed between these 2 cell populations. Similarly, between LSC and LDC, 3446 out of 9444 DMRs were annotated to 1741 genes differentially expressed between LSC and LDC (Fig. 3A). The majority of these overlapped genes were hypermethylated, with lower expression in LSC (LSC vs LIC: 58.3% DMRs; LSC vs LDC: 67.2% DMRs; Fig. 3B, blue dots), suggesting an inverse relationship between DNA methylation and gene expression changes. Pearson correlation analysis further confirmed that DNA methylation loss during differentiation was significantly correlated with expression increase of the nearest gene ( Fig. 3B; LSC vs LIC: r = -0.38, P < 2.2e-16; LSC vs LIC: r = -0.30, P < 2.2e-16; e = exponential). We found that the correlation between DNA methylation status and the expression of the nearest gene varied across the genome. For example, higher inverse correlations were found when the differentially methylated regions were near transcription start sites, including promoter, first exon, and first intron (Fig. 3C). To further investigate the possible role of DNA methylation in mediating the transcriptional activities of DNA regulatory elements, we performed ChIP-Seq for the active histone marks, H3K4me3 and H3K27Ac, in the 3 sorted cell populations (n = 3 patients). As expected, we found that during differentiation from LSC to LIC and LDC, loss of DNA methylation was strongly associated with the enrichment of H3K4me3 and H3K27Ac active histone marks (Fig. 3D and 3E).
We next performed enrichment analysis to identify pathways and networks enriched in the genes found to be both differentially methylated and expressed during LSC differentiation. We identified several processes that play crucial roles in LM tumorigenesis, including the WNT/β-catenin pathway, muscle development, BMP/ TGFβ signaling, and steroid hormone-related signaling (Fig. 3F) (7,33,36). Notably, genes hypermethylated and downregulated in LSC (blue bars) were more enriched in pathways regulating myometrium-and LM-associated functions (eg, muscle contraction, muscle development, hormone action, and ESR1 signaling), while genes hypomethylated and upregulated (yellow bars) were involved in stem cell-related pathways.

Genes critical for stem cell differentiation are hypermethylated and suppressed in LSC
Previous studies have shown that genes governing the differentiation process are repressed by Polycomb in human embryonic stem cells (16). We found significantly higher DNA methylation and lower mRNA levels at these Polycomb targets in LSC versus LIC and LDC (16), suggesting that DNA methylation might "lock" differentiation-associated genes in a repressed state in LSC ( Fig. 4A and 4B). To explore the potential role of DNA methylation during LSC differentiation, we investigated the DNA methylationmediated regulation of several Polycomb target genes with established roles in stem cell differentiation or LM pathogenesis. Among those Polycomb targets adjacent to hypermethylated regions and inhibited in LSC, we focused on ESR1, TIMP3, ROR2, and MYH11. The Estrogen/ESR1 signaling pathway plays important roles in LM tumorigenesis (37,38). TIMP3, tissue inhibitor of metalloproteinase 3, regulates extracellular matrix at the tumor-stromal interface (39) and is critical in driving differentiation and decreasing tumorigenesis in the mammary gland and muscle tissues (40,41). The ROR2 pathway plays complex roles in stem cell function and carcinogenesis (42,43), and MYH11 is a marker of fully differentiated smooth muscle cells (44).
MethylCap-Seq showed that the ESR1, TIMP3, ROR2, and MYH11 gene loci were hypermethylated at several intronic regions in LSC; the ESR1 gene was also hypermethylated at the promoter region in LSC (Fig. 4C). Opposite from the DNA methylation status, mRNA levels of ESR1, TIMP3, ROR2, and MYH11 were the lowest in LSC (Fig. 4D). To assess the effect of DNA methylation on the transcriptional activities of these genes, we treated individual cell populations with DNA methylation inhibitor 5′-Aza (100 nM) for 96 hours. 5′-Aza treatment significantly increased the mRNA levels of these genes in LSC, suggesting that the transcriptional activity of genes significant for the differentiation process were inhibited by DNA methylation in LSC (Fig. 4E).

5′-Aza drives LSC differentiation to reduce stemness
We demonstrated that DNA methylation contributes to the expression changes of critical genes during LSC  differentiation. We then tested the ability of 5′-Aza to regulate LSC function and compared its effect with that of RU486, a progesterone antagonist shown to inhibit LM growth (33). We treated LM tissue explants with vehicle (DMSO), RU486 (1 µM), or 5′-Aza (100 nM) for 48 hours and analyzed the proportions of each LM cell population.
As shown in Fig. 5A and 5B, 5′-Aza treatment decreased around 40% of the LSC population (5.93 ± 1.38% vs 3.58 ± 1.01%). The treatment also decreased the LIC population and increased the LDC population compared to the vehicle-treated cells, whereas RU486 did not significantly change the LM cell composition. We also tested the effect of RU486 or 5′-Aza on the clonogenic activity of passage zero (unpassaged) primary LM cells, a marker of tumor stem cells (45). Cells were treated with vehicle (DMSO), RU486 (1 µM), or 5′-Aza (25 nM, 50 nM, or 100 nM) for 6 days, and 500 viable cells from each treatment group were plated in each well of a 12-well plate and cultured for 21 days without further treatment. We found that pretreatment with 5′-Aza markedly decreased colony formation in primary LM cells even at a very low dose (25 nM), whereas RU486 did not have a significant effect ( Fig. 5C and 5D). In addition, we compared the tumor initiation capacity of passage zero primary LM cells (1 x 10 6 viable cells) pretreated with vehicle, 5′-Aza, or RU486 for 6 days. Although the alteration of cell surface marker gene expression during in vitro culture hindered us from distinguishing cellular components of primary LM cells after culture, our previous studies and the current colony formation assay indicate the presence of LSC in cultured primary LM cells (7,46). We found that primary LM cells pretreated with 5′-Aza regenerated significantly smaller tumors (36.30 ± 3.57% of vehicle size) compared with RU486-pretreated (76.31 ± 1.86% of vehicle size) or vehicle-pretreated primary LM cells (Fig. 5E).
We further explored the molecular mechanisms underlying 5′-Aza-mediated LSC population decrease and LM growth inhibition. We performed RNA-Seq on vehicle-(DMSO) or 5′-Aza-(100 nM, 6 days) treated primary LM cells. The global transcriptome changed markedly with 5′-Aza treatment with clear clustering based on treatment status (n = 3 patients, Fig. 6A). We identified 3046 genes that were upregulated and 5120 genes that were downregulated in the 5′-Aza group compared with the vehicle group. When comparing these genes to the genes differentially expressed in LSC, we found that 5′-Aza regulated 41.5% differentially expressed genes in LSC versus LIC and 40.1% differentially expressed genes in LSC versus LDC, suggesting that DNA methylation plays a critical role in regulating transcription activities during stem cell differentiation in LM (Fig. 6B). Pathway enrichment analysis revealed that genes downregulated by 5′-Aza were related to cell cycle regulation, whereas genes upregulated by 5′-Aza were enriched in immune response activation, stem cell differentiation, and extracellular matrix remodeling pathways (Fig. 6C). Gene set enrichment analysis showed similar results, with 5′-Aza treatment significantly downregulating proliferation-related genes and upregulating differentiation-and endocrine-related genes (Fig. 6D), consistent with our in vivo findings. Together, these findings suggest that inhibition of DNA methylation decreases LM tumor regeneration ability, potentially through inducing LSC differentiation and exhaustion.

Discussion
We demonstrated that LSC harbor a unique transcriptome and DNA methylome. Globally, hypermethylation in LSC suppresses key genes essential for stem cell differentiation, including ESR1, TIMP3, ROR2, and MYH11. The DNA methylation inhibitor 5′-Aza upregulated expression of genes necessary for differentiation, decreased the stem cell population by stimulating LSC differentiation, and reduced LM regeneration.
DNA methylation is a common mechanism for cell programming in both normal tissue growth and tumor development (47). Integrated analyses of genome-wide DNA methylation maps and expression profiles have revealed a hierarchy of cellular differentiation in hematopoietic tissue and the breast (48,49). Transcriptomic analyses of LM populations have also strongly suggested a hierarchical differentiation order: LSC→LIC→LDC, while MethylCap-seq analysis showed a similar pattern of alterations in DNA methylation, particularly a transition from LSC to LIC or LDC. We found, however, very little difference between LIC and LDC with respect to DNA methylation, suggesting that a change in DNA methylation status is a critical epigenetic modification during the initial differentiation of LSC. We uncovered regions that had lost or gained methylation during LSC differentiation and that were overrepresented in CpG islands and shores, consistent with the DNA methylation pattern changes seen during intestinal stem cell differentiation observed using whole-genome shotgun bisulfite sequencing (Fig. 1E) (13). The functional significance of this regional preference in stem cell differentiation warrants further investigation. DNA methylation, the most common and well-studied epigenetic mark, has a major influence on transcription activity. In our study, we found an inverse relationship between DNA methylation and gene expression changes when comparing LSC with LIC and LDC. The correlation was stronger when the methylation changes took place near the transcription start sites. It is well-known that hypermethylation at the promoter, especially a CpG dense promoter, is an indicator of suppressed transcription activity, but the role of methylation within the gene body regions remains ambiguous (50). Previously, gene body methylation was reported to be positively correlated with gene expression, but recent studies have demonstrated that methylation at the first intron and exon possesses a negative influence on the transcription activity, which coincides with our observation (51)(52)(53). Importantly, hypermethylated regions in LSC were functionally associated with transcriptional suppression of genes enriched in pathways regulating differentiated cell function, such as muscle contractions and hormone action, which is consistent with other tumor stem cell systems (Fig. 3F) (14,15).
Consistent with previous reports showing that DNA hypomethylating agents can remove the epigenetic "brake" and allow stem cells to undergo differentiation (17), we found that the DNA methylation inhibitor 5′-Aza not only upregulated expression of genes imperative for differentiation but also functionally depleted the LSC population and inhibited its colony formation (20,54). 5′-Aza-induced LSC differentiation was supported by the observation that 5′-Aza treatment significantly increased the expression of markers of fully differentiated smooth muscle cells, such as MYH11 (44). Notably, we found that RU486 pretreatment did not significantly affect the LM cells' colony-forming activity and was less efficient in inhibiting LM regeneration compared with 5′-Aza treatment, suggesting that LSC deficient in steroid hormone receptors may remain quiescent, constituting a small population of cells resistant to RU486 treatment. Upon treatment cessation, the persistent stem cells potentially expand and differentiate, eventually leading to LM regrowth, which may partially explain the clinically observed frequent LM regrowth after hormone therapy, a phenomenon similar to endocrine therapy resistance in breast cancer (55). Our findings that 5′-Aza treatment-activated genes related to endocrine processes, including ESR1 (Figs. 4E and 6C), indicated that it may sensitize LSC to hormone treatment. Here, we used 5′-Aza to inhibit DNA methylation as a proof-of-concept study that DNA methylation affects LSC differentiation and stemness. As a demethylation reagent, 5′-Aza is toxic (56). Thus, more studies are needed to understand the differentiation processes of LSC to identify novel interventions, which can potentiate hormone treatments to prevent LM regrowth.
In this study, we utilized the MethylCap-Seq method to evaluate the genome-wide methylation landscape of each LM cell population. MethylCap-Seq uses MBD protein to immunoprecipitate methylated DNA fragments. Besides whole-genome bisulfite sequencing, which is the gold standard technique for methylation profiling, MethylCap-Seq has the largest coverage of the genome (57). Using this method, we identified approximately 20 000 differentially methylated regions among the 3 LM cell populations, which were widely distributed across all genomic features. However, compared with methylation microarrays (EPIC array) and reduced representation bisulfite sequencing, MethylCap-Seq has the limitation of precisely quantifying methylation levels at the singlenucleotide resolution (57). In a recent study, De Mayer et al compared results obtained by using MethylCap-seq and methylation microarrays and confirmed the consistency of the 2 methods (58).
A strength of our study is the performance of genomewide studies in freshly sorted populations of LM cells to explore the mechanisms underlying LSC differentiation. The findings derived from patient samples undergoing minimal experimental manipulation better represent in vivo conditions, producing insights with high translational relevance. Unfortunately, technical challenges with long-term in vitro culture and propagation of LSC prohibited us from performing detailed downstream mechanistic studies specifically on stem cells. Further work is needed to optimize and establish a LSC culture system that will help us better understand the functional properties of these cells.
In summary, we found that hypermethylation maintains LSC stemness through suppression of genes critical for differentiation. This finding not only sheds mechanistic light on the molecular mechanisms that drive LSC differentiation but also suggests a treatment strategy targeting the methylation status of the LSC population to prevent the formation of new LM tumors.