Introduction {#s1}
============

Uterine leiomyomas are the most common uterine tumors in reproductive-age women with a prevalence of about 25% [@pone.0066632-Stewart1]. Uterine leiomyomas frequently cause serious gynecological problems such as pelvic pain, menorrhagia, dysmenorrhea, infertility and recurrent pregnancy loss [@pone.0066632-Stewart1], [@pone.0066632-Bajekal1]. In addition, uterine leiomyomas are the most common indication for hysterectomy. Risk factors for uterine leiomyomas include African descent, high body mass index, meat consumption, early menarche, hypertension and a history of pelvic inflammatory disease. On the other hand, factors that lower the risk include use of hormonal contraception, smoking, giving birth and consumption of green vegetables [@pone.0066632-Chiaffarino1], [@pone.0066632-Faerstein1], [@pone.0066632-Faerstein2]. These findings suggest that both genetic and environmental factors are involved in the development of uterine leiomyomas. In addition, uterine leiomyomas often show multifocal tumorigenesis with various sizes from the corresponding myometrium. These findings suggest that smooth muscle cells of normal myometrium in the uterus with leiomyomas already acquire the potential in molecular levels to develop into leiomyomas in future.

DNA methylation is one of the well characterized epigenetic marks and plays a crucial role in the regulation of gene expression. DNA methylation is specific to each cell type and has been used to characterize abnormal cells [@pone.0066632-Shiota1], [@pone.0066632-Lieb1], [@pone.0066632-Shiota2]. Maintaining the specific DNA methylation profile of the cell is necessary for cellular integrity, and alterations in DNA methylation may have serious health consequences. Environmental factors can be shown to affect DNA methylation [@pone.0066632-Iwatani1], [@pone.0066632-Wu1]. We previously demonstrated that uterine leiomyomas have an aberrant DNA methylation profile using genome-wide DNA methylation analysis methods [@pone.0066632-Yamagata1], [@pone.0066632-Maekawa1]. Navarro et al. also analyzed the DNA methylation patterns in uterine leiomyomas using genome-wide analysis method in African American women [@pone.0066632-Navarro1]. They found several aberrantly methylated and concomitantly expressed genes in leiomyomas, and suggested that DNA methylation plays a key role in the pathogenesis of uterine leiomyomas by altering the normal myometrial mRNA expression profile. The study by Navarro et al. examined one or two CpG sites in the promoter region in each gene. Recently, a new type of genome-wide DNA methylation analysis has been developed, which covers about 16 CpG sites from the promoter to the gene body region. Therefore, the new analysis will provide more detailed information on genome-wide DNA methylation profiles in uterine leiomyomas.

Although uterine leiomyomas develop only after menarche and their growth depends on estrogen, the mechanism is not fully clarified. The biological effect of estrogen is mediated by estrogen receptor (ER), which is a transcription factor [@pone.0066632-Stewart1]. DNA methylation of the gene promoter region directly affects the response to transcription factors [@pone.0066632-Hitchler1]. Thus, aberrant DNA methylation of the promoter regions of genes targeted by ER alpha directly may affect the response to estrogen. Recently, exposure to environmental estrogen during reproductive tract development has been shown to cause genes targeted by ER alpha to become hyper-responsive to estrogen in the adult myometrium [@pone.0066632-Greathouse1]. These findings raise the possibility that aberrant DNA methylation of the promoter of ER alpha-target genes may cause the aberrant responses to estrogen exposure in uterine leiomyomas, which may be involved in the development of uterine leiomyomas after menarche.

We previously reported that aberrant DNA hypomethylation is enriched on the X chromosome in uterine leiomyomas compared with the adjacent normal myometrium [@pone.0066632-Maekawa1]. Interestingly, the X chromosome was found to have more hypomethylated genes than other chromosomes [@pone.0066632-Maekawa1]. It has been reported that breast cancers, ovarian cancers, and cervical cancers have aberrant DNA hypomethylation on the X chromosome such as loss of inactive X chromosome or aberrant replication of active X chromosome [@pone.0066632-Kawakami1], [@pone.0066632-Sirchia1]. Therefore, an analysis of the genotype of the X chromosome might provide clues to the cause of uterine leiomyomas.

The main objective of this study is to determine whether DNA methylation patterns are associated with the pathogenesis and development of uterine leiomyomas. To this end, we 1) compared the DNA methylation profiles and mRNA expression profiles among uterine leiomyomas and myometrium with or without leiomyomas using an advanced genome-wide DNA methylation analysis method, 2) attempted to identify a subset of genes whose differential DNA methylation in the promoter is correlated with differential mRNA expression, 3) identified ER alpha-target genes that have aberrant DNA methylation in the promoter, and 4) examined the genotype of the X chromosome and the chromosome distribution of aberrant DNA methylation in uterine leiomyomas.

Materials and Methods {#s2}
=====================

Ethics Statement {#s2a}
----------------

This study was reviewed and approved by the Institutional Review Board of Yamaguchi University Graduate School of Medicine. Written informed consent was obtained from the participants before the collection of any samples, and the specimens were irreversibly de-identified. All experiments handling human tissues were performed in accordance with Tenets of the Declaration of Helsinki.

Tissue Preparations {#s2b}
-------------------

Paired specimens of leiomyoma and adjacent normal myometrium were obtained from 10 Japanese women who had a single leiomyoma nodule to limit biological heterogeneity ([Table 1](#pone-0066632-t001){ref-type="table"}). The women underwent hysterectomy, and their ages were 44.1+/−4.5 years old (mean +/− SD). The size of leiomyoma ranged from 40 to 160 mm (mean +/− SD; 85.5+/−34.0 mm). All patients showed intramural leiomyoma. As control myometrium tissues, we obtained normal myometrium tissue from 3 women who underwent hysterectomy under the diagnosis of early stage of cervical cancer (Control 1, Control and Control 3 in [Table 1](#pone-0066632-t001){ref-type="table"}). The ages of control women were 38.3+/−5.0 years old (mean +/− SD), and did not significantly differ from women with uterine leiomyomas. None of the women had received previous treatment with sex steroid hormones or gonadotropin releasing hormone analogs. Specimens were dissected immediately from the full thickness uterine wall without the endometrium and serosa after removal of the uterus, immersed in liquid nitrogen and stored at -80°C until DNA/RNA extraction. Leiomyomas (L1, L2 and L3) and adjacent myometrium (M1, M2 and M3) obtained from 3 cases (Case 1, Case 2 and Case 3 in [Table 1](#pone-0066632-t001){ref-type="table"}) were examined by following genome-wide DNA methylation analysis and mRNA expression analysis. The myometrium from 3 control women were also examined (C1, C2, and C3). The other 7 paired samples of leiomyoma and adjacent myometrium were used in bisulfite restriction mapping and RT-PCR.

10.1371/journal.pone.0066632.t001

###### Patients with uterine leiomyomas and characteristics of their leiomyomas.

![](pone.0066632.t001){#pone-0066632-t001-1}

               Age (y.o.)   Menstrual phase      Type      Diameter (mm)
  ----------- ------------ ----------------- ------------ ---------------
  Control 1        43         Follicular          --      
  Control 2        33         Follicular          --      
  Control 3        39           Luteal            --      
  Case 1           49           Luteal        intramural        160
  Case 2           49         Follicular      intramural        90
  Case 3           40           Luteal        intramural        80
  Case 4           43         Follicular      intramural        90
  Case 5           39         Follicular      intramural        70
  Case 6           40         Follicular      intramural        85
  Case 7           38           Luteal        intramural        120
  Case 8           48         Follicular      intramural        40
  Case 9           47         Follicular      intramural        60
  Case 10          48           Luteal        intramural        60

Controls 1--3 and Cases 1--3 were used for analysis with Illumina Infinium HumanMethylation450 BeadChip and Human Gene 1.0 ST array. Cases 1--10 were used for validation of the genes with differential DNA methylation and mRNA expression.

Illumina Infinium HumanMethylation450 BeadChip Assay {#s2c}
----------------------------------------------------

Genomic DNA was isolated from 20 mg frozen tissues using the Qiagen Genomic DNA (Qiagen, Valencia, CA, USA). The DNA methylation analysis was performed using the Illumina infinium assay with the HumanMethylation450 BeadChip (Illumina, San Diego, CA, USA), which interrogates a total of 482,421 CpG sites spread across the distal promoter regions of transcription start sites to 3′-UTR of consensus coding sequences. Methylated and unmethylated signals were used to compute beta values, which are quantitative scores of the DNA methylation levels, ranging from "0,", thus indicating completely unmethylated, to "1," indicating completely methylated. The BeadChip was scanned on a BeadArray Reader (Illumina) according to the manufacturer's instructions. CpG sites with "detection p values" \>0.05 (computed from the background based on negative controls) and CpG sites on Y chromosome were eliminated from further analysis, leaving 482,005 CpGs valid for use with the nine samples tested. To perform hierarchical clustering analysis and principal component analysis, we extracted the CpG sites in which the difference of maximum beta value and minimum beta value in 9 tissues is greater than 0.3. To identify the CpGs that showed methylation changes between leiomyoma and adjacent myometrium, the difference of beta value of CpGs were calculated in each case. Then, CpGs expressing greater than 0.3 were extracted. The microarray data of DNA methylation is available at the Gene Expression Omnibus Web site (<http://www.ncbi.nlm.nih.gov/geo/>) under accession No. GSE45187.

Transcriptome Analysis {#s2d}
----------------------

Gene expression was analyzed using a GeneChip® Human Gene 1.0ST Array (Affymetrix, Santa Clara, CA, USA) containing 764,885 probes supporting 28,869 genes. Target cDNA was prepared from 200 ng of total RNA with the Ambion® WT Expression kit (Ambion, Austin, TX, USA) and the Affymetrix® GeneChip WT Terminal Labeling kit (Affymetrix). Hybridization to the microarrays, washing, staining and scanning were performed using the GeneChip® system (Affymetrix) composed of the Scanner 3000 7G Workstation Fluidics 450 and the Hybridization Oven 645. The scanned image data were processed using a Gene Expression Analysis with the Partek® Genomics Suite 6.5 software program (Partech, Munster, Germany). To perform hierarchical clustering analysis and principal component analysis, we extracted the gene in which the ratio of maximum value to minimum value in 9 tissues is greater than 1.5-fold. To identify the genes that show the difference in mRNA expression between leiomyoma and adjacent myometrium, a fold-change greater than 1.5-fold or lesser than 0.67-fold was recognized as a significant difference. The microarray data of mRNA expression is available at the Gene Expression Omnibus Web site (<http://www.ncbi.nlm.nih.gov/geo/>) under accession No. GSE45188.

Bisulfite Restriction Mapping {#s2e}
-----------------------------

Bisulfite reactions were performed as previously reported [@pone.0066632-Asada1]. The bisulfite-converted DNA was amplified by PCR using the following sets of primers: IRS1 F: 5′-TGTTTTATTTTAAATTTTTAGTGGAGAGTAG-3′, R: 5′-AATATCCCCCACCCAAACTATC-3′; COL4A1 F: 5′-GGAGGGAGGTAGAGGATTTATGT-3′, R: 5′-AAAAACTCCTAAATAAATCTAAAAAAAACT-3′; GSTM5 F: 5′-GGGAGGGGGTTTATTGATTTTAGTT-3′, R: 5′-ACCAATCCTACAACTACTCCACACATC-3′. The thermocycling program was an initial cycle of 94°C for 10 min, then 43 cycles of 94°C for 30 sec, 59°C for 30 sec and 72°C for 1 min and finally 10 min of final extension at 72°C. For restriction mapping, 1/2 each of the PCR products was treated with or without restriction enzyme, and the resulting DNA fragments were assessed by agarose gel electrophoresis.

Real-time Quantitative RT-PCR {#s2f}
-----------------------------

Total RNAs were isolated from tissues using TRIzol reagent (Invitrogen, Carlsbad, CA, USA) and reverse-transcribed using a QuantiTect Reverse Transcription Kit (Qiagen) according to the manufacturer's protocol. The synthesized DNA was subjected to PCR reactions using the following sets of primers: IRS1 F: 5′- CAAGACCATCAGCTTCGTGA -3′, R: 5′- AGAGTCATCCACCTGCATCC -3′; COL4A1 F: 5′- TGGTGACAAAGGACAAGCAG -3′, R: 5′- GGTTCACCCTTTGGACCTG -3′; GSTM5 F: 5′- AATGCCATCCTGCGCTAC -3′, R: 5′- TTCTCCAAAATGTCCACACG -3′. All PCRs were performed using SYBR Premix Ex Taq (TAKARA) and a LightCycler (Roche Applied Science, Basel, Switzerland). A primer set for glyceraldehyde-3-phosphate dehydrogenase (GAPDH) was used as an internal control. All samples were run in duplicate. Melting curves of the products were obtained after cycling by a stepwise increase of temperature from 55 to 95°C. At the end of 40 cycles, reaction products were separated electrophoretically on an agarose gel and stained with ethidium bromide for visual conformation of the PCR products.

X chromosome Genotyping {#s2g}
-----------------------

The analysis was done using a panel of primers specific for the following X-linked short tandem repeats markers, obtained from Genome Database (<http://www.gdb.org/>) and National Center for Biotechnology Information (<http://www.ncbi.nlm.nih.gov/genomes/static/euk_g.html>): DXS7132, DXS6789, DXS8377 and DXS6807 spanning p and q arms. After PCR using one fluorescent-labeled primer, the products were separated electrophoretically on the Fluorescent Capillary System ABI PRISM 310 and analyzed with GeneScan software (Applied Biosystems). The allelic status in each case was examined by comparing the allelic statuses of leiomyoma and adjacent myometrium.

Bioinformatics {#s2h}
--------------

MultiExperiment Viewer (MeV in the TM4 microarray Software Suite; <http://www.tm4org/mevhtml>) was used for principal component analysis (PCA) [@pone.0066632-Saeed1]. The analysis for ER alpha-target genes analysis was performed using the MAPPER database (<http://mapper.chip.org/>) [@pone.0066632-Marinescu1]. Chromosome distribution analysis was performed using DAVID 6.7 (<http://david.abcc.ncifcrf.gov/home.jsp>) [@pone.0066632-Huangda1]. Integrative Genomics Viewer (<http://www.broadinstitute.org/igv/>) was used for chromosomal distribution map generation [@pone.0066632-Robinson1]. QUMA (<http://quma.cdb.riken.jp/>) was used to analyze the bisulfite sequencing data [@pone.0066632-Kumaki1].

Results {#s3}
=======

Analysis of DNA Methylation Profiles of Leiomyomas and Myometrium with and without Leiomyomas {#s3a}
---------------------------------------------------------------------------------------------

To compare the genome-wide DNA methylation patterns among leiomyomas (L1, L2, and L3), the corresponding myometrium with leiomyomas (M1, M2, and M3), and the myometrium without leiomyomas (C1, C2, and C3), the hierarchical clustering was analyzed according to the DNA methylation status. Myometrium without leiomyomas (C1, C2, and C3) and myometrium with leiomyomas (M1, M2, and M3) showed a similar DNA methylation profiles, and these 6 tissues were classified into the same cluster group ([Figure 1A](#pone-0066632-g001){ref-type="fig"}). On the other hand, leiomyomas (L1, L2, and L3) were clustered separately from myometrium with or without leiomyomas ([Figure 1A](#pone-0066632-g001){ref-type="fig"}). The DNA methylation profiles of these 9 tissues were further compared with a principal component analysis. In [Figure 1B](#pone-0066632-g001){ref-type="fig"}, the component-1 and component-2 axes clearly distinguished leiomyomas (L1, L2, and L3) from myometrium with (M1, M2, and M3) or without (C1, C2, and C3) leiomyomas, while myometrium with leiomyomas (M1, M2, and M3) was very close to myometrium without leiomyomas (C1, C2, and C3). These results indicate that myometrium with and without leiomyomas have quite similar DNA methylation profiles whereas leiomyomas are clearly distinguished from the myometrium in terms of DNA methylation profiles.

![DNA methylation profiling and mRNA expression profiling of leiomyomas and myometrium with and without leiomyomas.\
DNA methylation profiles and mRNA expression profiles of leiomyomas (L1, L2 and L3), myometrium with leiomyomas (M1, M2 and M3) and myometrium without leiomyomas (C1, C2 and C3) were compared using hierarchical clustering analyses and principal component analyses. A: Hierarchical clustering analyses according to DNA methylation profiles. The Heat map in hierarchical clustering analysis indicates DNA methylation levels from unmethylated (blue) to completely methylated (yellow). Distances of DNA methylation pattern (Euclidean Distances) were calculated by MultiExperiment Viewer. B: Principal component analyses according to DNA methylation profiles. Vertical axis and horizontal axis show principal component numbers, respectively. The principal component analyses were performed using MultiExperiment Viewer. C: Hierarchical clustering analyses according to mRNA expression profiles. The Heat map in hierarchical clustering analysis indicates mRNA expression levels from low (blue) to high (yellow). Distances of mRNA expression pattern (Euclidean Distances) were shown on the left side. D: Principal component analyses according to mRNA expression profiles. Vertical axis and horizontal axis show principal component numbers, respectively.](pone.0066632.g001){#pone-0066632-g001}

We also analyzed the hierarchical clustering and principal components of the mRNA expression profiles. In the hierarchical clustering analysis, myometrium with and without leiomyomas were classified into the same cluster group, while leiomyomas were clearly separated from myometrium with or without leiomyomas ([Figure 1C](#pone-0066632-g001){ref-type="fig"}). In the principal component analysis, the component-1 axis clearly distinguished leiomyomas (L1, L2, and L3) from myometrium with (M1, M2, and M3) or without (C1, C2, and C3) leiomyomas, whereas myometrium with leiomyomas (M1, M2, and M3) and myometrium without leiomyomas (C1, C2, and C3) were not separated into the same group by the component-2 axis ([Figure 1D](#pone-0066632-g001){ref-type="fig"}). The data in [Figures 1A and 1B](#pone-0066632-g001){ref-type="fig"} suggest that profiling by DNA methylation is more useful than profiling by mRNA expression in defining cell identity.

Analysis of DNA Methylation and mRNA Expression in Leiomyoma and Matched Adjacent Myometrium {#s3b}
--------------------------------------------------------------------------------------------

We first compared the DNA methylation profiles between leiomyomas and adjacent myometrium in each of the three cases examined in this study. Case1, Case2 and Case3 had 2,386, 1,327 and 3,078 genes, respectively, that were less methylated in the leiomyomas than in the myometrium ([Figure 2A](#pone-0066632-g002){ref-type="fig"}, hypomethylated genes). Of these hypomethylated genes, 478 genes were shared by the three cases. Case1, Case2 and Case3 had 2,390, 1,702 and 3,567 genes, respectively, that were more methylated in the leiomyomas than in the myometrium ([Figure 2A](#pone-0066632-g002){ref-type="fig"}, hypermethylated genes). Of these hypermethylated genes, 1,014 genes were shared by the three cases.

![Aberrantly methylated or expressed genes in Leiomyomas.\
A: Venn diagrams show the number of aberrantly hypomethylated or hypermethylated genes in 3 paired samples. B: Venn diagrams show the number of aberrantly upregulated or downregulated genes in 3 paired samples. C: Integration of DNA methylation and mRNA expression data in uterine leiomyoma and adjacent myometrium. The numbers of genes which have both aberrant methylation and aberrant mRNA expression are shown. D: Functional analysis of 24 genes that have aberrant hypomethylation and upregulation of mRNA. The functional categories are shown with the number of included genes and p-values. Functional analysis was performed using IPA. E: Functional analysis of 65 genes that have aberrant hypermethylation and downregulation of mRNA. The functional categories are shown with the number of included genes and p-values. Functional analysis was performed using IPA.](pone.0066632.g002){#pone-0066632-g002}

Messenger RNA expression profiles were also analyzed among the three cases. Case1, Case2 and Case3 had 1,139, 1,000 and 1,401 genes, respectively, that were upregulated in the leiomyomas compared with the myometrium ([Figure 2B](#pone-0066632-g002){ref-type="fig"}, upregulated genes). Of these upregulated genes, 349 genes were shared by the three cases. Case1, Case2 and Case3 had 2,173, 1,058 and 1,197 genes, respectively, that were downregulated in the leiomyomas compared with the myometrium ([Figure 2B](#pone-0066632-g002){ref-type="fig"}, downregulated genes). Of these downregulated genes, 576 genes were shared by the three cases.

We then selected 120 genes that had both aberrant DNA methylation and aberrant mRNA expression. Of these genes, 24 genes were less methylated and transcriptionally upregulated in the leiomyomas compared with the myometrium, and 65 genes were more methylated and transcriptionally downregulated in the leiomyomas compared with the myometrium. According to the Ingenuity Pathways Analysis (IPA) [@pone.0066632-Feng1], in the 24 hypomethylated and upregulated genes, the most specific pathway was "cancer process" in the section of diseases and disorders ([Figure 2D](#pone-0066632-g002){ref-type="fig"}). The genes included in the category "cancer process" are shown in [Table 2](#pone-0066632-t002){ref-type="table"}. Collagen-related genes such as COL4A1 and COL4A2 have been previously reported to be upregulated in leiomyomas [@pone.0066632-Vanharanta1], [@pone.0066632-Gilden1], [@pone.0066632-Weston1]. In the 65 hypermethylated and downregulated genes, the most specific pathway was also "cancer process" ([Figure 2E](#pone-0066632-g002){ref-type="fig"}). The genes included in the category "cancer process" are shown in [Table 2](#pone-0066632-t002){ref-type="table"}. Glutathione S-transferase mu 5 (GSTM5) is a member of glutathione S-transferase family and protects cells as an anti-oxidant enzyme against reactive oxygen species. GSTM5 is also known to be involved in cancer development [@pone.0066632-Peng1].

10.1371/journal.pone.0066632.t002

###### Aberrantly methylated and expressed genes contained in \"cancer process\".

![](pone.0066632.t002){#pone-0066632-t002-2}

  Gene symbol                                                                                  Gene                               fold change
  --------------------------------------------------------------- -------------------------------------------------------------- -------------
  **Hypomethylated and transcriptionally upregulated genes**                                                                     
  CACNA1C                                                          calcium channel, voltage-dependent, L type, alpha 1C subunit      1.974
  COL4A1                                                                            collagen, type IV, alpha 1                       1.822
  COL4A2                                                                            collagen, type IV, alpha 2                       2.022
  COL6A3                                                                            collagen, type VI, alpha 3                       1.950
  CYP1B1                                                              cytochrome P450, family 1, subfamily B, polypeptide 1          1.844
  IRS1                                                                             insulin receptor substrate 1                      2.395
  KIAA1199                                                                                   KIAA1199                                8.507
  POPDC2                                                                            popeye domain containing 2                       4.952
  PRL                                                                                       prolactin                                8.828
  RAD51L1(RAD51B)                                                                  RAD51-like 1 (S. cerevisiae)                      3.712
  TACR2                                                                               tachykinin receptor 2                          2.862
  UNC5D                                                                            unc-5 homolog D (C. elegans)                      3.851
  VCAN                                                                                       versican                                2.677
  WBSCR17                                                                 Williams-Beuren syndrome chromosome region 17              2.177
  **Hypermethylated and transcriptionally downregulated genes**                                                                  
  AIM1                                                                                 absent in melanoma 1                         −1.917
  ANXA3                                                                                     annexin A3                              −4.190
  CD74                                                               CD74 molecule, major histocompatibility complex, class I       −1.620
  CFB                                                                                  complement factor B                          −2.170
  DAPK1                                                                         death-associated protein kinase 1                   −2.007
  DPP4                                                                                dipeptidyl-peptidase 4                        −2.825
  DUSP6                                                                           dual specificity phosphatase 6                    −2.332
  EFEMP1                                                             EGF-containing fibulin-like extracellular matrix protein       −9.696
  EPAS1                                                                          endothelial PAS domain protein 1                   −2.369
  EZR                                                                                         ezrin                                 −2.227
  FOXP1                                                                                  forkhead box P1                            −0.626
  GATA2                                                                               GATA binding protein 2                        −2.460
  GSTM5                                                                           glutathione S-transferase mu 5                    −9.998
  HTATIP2                                                                    HIV-1 Tat interactive protein 2, 30 kDa                −3.253
  IGF2BP3                                                              insulin-like growth factor 2 mRNA binding protein 3          −1.684
  KAT2B                                                                           K(lysine) acetyltransferase 2B                    −2.715
  KDR                                                                             kinase insert domain receptor                     −2.565
  LDB2                                                                                 LIM domain binding 2                         −2.034
  LIMA1                                                                           LIM domain and actin binding 1                    −1.759
  MEOX2                                                                               mesenchyme homeobox 2                         −1.850
  MYEF2                                                                             myelin expression factor 2                      −2.304
  NR3C1                                                                  nuclear receptor subfamily 3, group C, member 1            −2.132
  NR4A2                                                                  nuclear receptor subfamily 4, group A, member 2            −1.931
  NR4A3                                                                  nuclear receptor subfamily 4, group A, member 3            −2.136
  NR5A2                                                                  nuclear receptor subfamily 5, group A, member 2            −2.317
  NRN1                                                                                      neuritin 1                              −3.641
  NTRK2                                                                   neurotrophic tyrosine kinase, receptor, type 2            −2.175
  NUAK1                                                                          NUAK family, SNF1-like kinase, 1                   −2.568
  PLCE1                                                                             phospholipase C, epsilon 1                      −4.468
  SDPR                                                                              serum deprivation response                      −2.895
  SORBS2                                                                        sorbin and SH3 domain containing 2                  −2.291
  SPTBN1                                                                        spectrin, beta, non-erythrocytic 1                  −1.957
  TIAM1                                                                     T-cell lymphoma invasion and metastasis 1               −1.730

The fold change was calculated as values of leiomyoma relative to adjacent myometrium in each case in the mRNA expression microarray.

Validation of the Genes with Differential DNA Methylation and mRNA Expression {#s3c}
-----------------------------------------------------------------------------

Among the hypomethylated and upregulated genes, we focused on insulin receptor substrate 1 (IRS1) and COL4A1 which have aberrant DNA hypomethylation in the gene body closed to the promoter region. IRS1 is initially characterized as a cytosolic adaptor protein involved in insulin receptor (INSR) and insulin-like growth factor 1 receptor (IGF1R) signaling. More recently, it has been shown to be involved in proliferation and transformation of cancer cells [@pone.0066632-Morelli1], [@pone.0066632-Esposito1]. First, we examined the DNA methylation status of the IRS1 and COL4A1 using combined bisulfite restriction analysis in Case1, Case2, Case3 and additional 7 paired samples of leiomyoma and matched myometrium. As shown in [Figure 3A](#pone-0066632-g003){ref-type="fig"}, the methylation levels of the IRS1 were lower than those in myometrium in all the 10 samples. In addition, mRNA levels of IRS1 in leiomyomas were higher than those in myometrium in all the 10 samples ([Figure 3D](#pone-0066632-g003){ref-type="fig"}). In COL4A1, the methylation levels in leiomyomas were lower than those in myometrium in all samples except Case7 ([Figure 3B](#pone-0066632-g003){ref-type="fig"}), and mRNA levels in leiomyomas were higher than those in myometrium in all samples except Case7 ([Figure 3E](#pone-0066632-g003){ref-type="fig"}).

![Analyses of DNA methylation levels and mRNA expression levels of IRS1, COL4A1 and GSTM5.\
DNA methylation levels and mRNA expression levels of IRS1, COL4A1 and GSTM5 were examined on 10 paired samples with myometrium (Myomet) and leiomyoma (Leio). DNA methylation levels of IRS1 (A), COL4A1 (B) and GSTM5 (C) were evaluated by bisulfite restriction mapping, quantified by densitometry, and expressed as %methylated DNA. The levels of mRNA expression of IRS1 (D), COL4A1 (E) and GSTM5 (F) were quantified by real-time quantitative RT-PCR and were normalized to GAPDH. The mRNA level in the myometrium of Case1 was expressed as 1.](pone.0066632.g003){#pone-0066632-g003}

To examine the possibility that IRS1 mRNA expression is under the regulation of DNA methylation, immortalized human uterine smooth muscle cells (1×10^6^ cells/25 cm^2^ tissue culture flask) were treated by an inhibitor of DNA methylation (5-aza-dC, 1 uM) for 96 h. IRS1 mRNA levels were significantly (p\<0.05) increased by 5-aza-dC, suggesting that DNA methylation is involved in the regulation of IRS1 mRNA expression (data not shown). The treatment with 5-aza-dC caused no significant changes in cell morphology and cell proliferation.

Among the hypermethylated and downregulated genes, we focused on GSTM5 which has aberrant DNA hypomethylation in the promoter region. As shown in [Figure 3C and 3F](#pone-0066632-g003){ref-type="fig"}, the methylation levels of the GSTM5 promoter in leiomyomas were higher than those in myometrium while mRNA levels of GSTM5 in leiomyomas were remarkably lower than those in myometrium in all the 10 samples.

ER Alpha- target Genes with Aberrant DNA Methylation and mRNA Expression {#s3d}
------------------------------------------------------------------------

In the aberrantly methylated and expressed 120 genes, 22 genes which have the consensus sequences of ER response element (ERE) in the promoter region (from −1,000 bp to the transcription start site) were extracted using the MAPPER database [@pone.0066632-Marinescu1] ([Table 3](#pone-0066632-t003){ref-type="table"}). In addition to COL4A1, COL6A3 and GSTM5, DAPK1 (death-associated protein kinase 1) and NUAK1 (novel kinase family 1), which were hypermethylated and transcriptionally downregulated, were included. DAPK1, which is a tumor suppressor gene, has been shown to be associated with apoptosis, and the downregulation and DNA hypermethylation were reported in several tumors [@pone.0066632-Martoriati1], [@pone.0066632-Raval1], [@pone.0066632-Claus1], [@pone.0066632-Missaoui1]. NUAK1 is known to possess tumor suppressive properties through the control of cellular senescence [@pone.0066632-Bernard1], [@pone.0066632-Hou1].

10.1371/journal.pone.0066632.t003

###### Estrogen receptor alpha-target genes with aberrant DNA methylation in the promoter and aberrant mRNA expression.

![](pone.0066632.t003){#pone-0066632-t003-3}

  Gene                                                                                                   mRNA expression[\#](#nt103){ref-type="table-fn"}
  -------------------------------------- -------------------------------------------------------------- --------------------------------------------------
  **Hypomethylated ERa target genes**                                                                   
  COL4A1                                                   collagen, type IV, alpha 1                                           ↑
  COL6A3                                                   collagen, type VI, alpha 3                                           ↑
  RPL39                                                      ribosomal protein L39                                              ↑
  ZMAT3                                                    zinc finger, matrin-type 3                                           ↑
  OXTR                                                         oxytocin receptor                                                ↓
  **Hypermethylated ERa target genes**                                                                  
  BCAN                                                              brevican                                                    ↑
  KIF5C                                                     kinesin family member 5C                                            ↑
  NPTX2                                                      neuronal pentraxin II                                              ↑
  TFAP2C                                                transcription factor AP-2 gamma                                         ↑
  SCIN                                                             scinderin                                                    ↑
  THBS2                                                         thrombospondin 2                                                ↑
  OCIAD2                                                    OCIA domain containing 2                                            ↓
  CRHBP                                         corticotropin releasing hormone binding protein                                 ↓
  NUAK1                                                 NUAK family, SNF1-like kinase, 1                                        ↓
  CCDC68                                                coiled-coil domain containing 68                                        ↓
  NR5A2                                         nuclear receptor subfamily 5, group A, member 2                                 ↓
  ELTD1                                   EGF, latrophilin and seven transmembrane domain containing 1                          ↓
  FAM162B                                        family with sequence similarity 162, member B                                  ↓
  IGF2BP3                                     insulin-like growth factor 2 mRNA binding protein 3                               ↓
  DAPK1                                                death-associated protein kinase 1                                        ↓
  GSTM5                                                  glutathione S-transferase mu 5                                         ↓
  EZR                                                                ezrin                                                      ↓

Microarray mRNA expression statuses of leiomyoma relative to adjacent myometrium are shown.

Analysis of X chromosome Genotype {#s3e}
---------------------------------

Recent reports have shown that in females, some cancer cells lose an inactive X chromosome and replicate an active X chromosome. It is called "loss-of-inactive X chromosome and gain-of-active X chromosome" [@pone.0066632-Kawakami1], [@pone.0066632-Sirchia1]. The cancer cells affected by this event have two active X chromosomes, resulting in the DNA hypomethylation status in X chromosome. Therefore, we decided to investigate whether uterine leiomyomas suffer this event and result in the hypomethylation status of X chromosome. Since the cells affected by this event should have two genetically same X chromosomes, we investigated the heterozygosity of X chromosomes in leiomyomas using following X-linked short tandem repeats markers: DXS7132, DXS6789, DXS8377 and DXS6807 spanning p and q arms. The result showed that all three cases showed heterozygosity in both myometrium and leiomyomas in all markers ([Figure 4](#pone-0066632-g004){ref-type="fig"}), indicating that both the leiomyoma and myometrium have a normal X chromosome genotype.

![X chromosomal microsatellite analysis of leiomyoma and adjacent myometrium.\
Chromatographs using four tandem repeat markers (DXS7132, DXS6789, DXS8377 and DXS6807) are indicated. Chromatographs of DNA from leiomyoma (L) and adjacent myometrium (M) obtained from Case1 (Top), Case2 (Middle) and Case3 (Bottom) are shown. The PCR products of DXS7132, DXS6789, DXS8377 and DXS6807 were separated electrophoretically on a Fluorescent Capillary System ABI PRISM 310 and analyzed with GeneScan software (Applied Biosystems). The allelic status in each case was examined by comparing the allelic statuses of the leiomyoma and adjacent myometrium.](pone.0066632.g004){#pone-0066632-g004}

Chromosome Distribution of Aberrantly Methylated CpG Sites {#s3f}
----------------------------------------------------------

Aberrant hypomethylation and hypermethylation occurred on all chromosomes in all three cases ([Figure 5](#pone-0066632-g005){ref-type="fig"}). Hypermethylated CpG sites (red) were enriched compared with hypomethylated CpG sites (green) in the autosomes in all three cases ([Figure 5](#pone-0066632-g005){ref-type="fig"}). On the other hand, X chromosome was preferentially enriched with hypomethylated CpG sites (green) in all three cases with remarkably low p values (5.31E-95 in Case 1, 1.34E-55 in Case 2, 1.63E-14 in Case 3), and the distribution of the hypomethylated CpG sites was observed throughout the X chromosome ([Figure 5](#pone-0066632-g005){ref-type="fig"}).

![The chromosomal distribution of aberrantly methylated CpGs in uterine leiomyomas.\
The locations of CpGs, in which the difference of beta value between leiomyoma and adjacent myometrium are greater than 0.3, are shown with green bar (hypomethylated CpG in leiomyoma) or red bar (hypermethylated CpG in leiomyoma). The chromosomal distribution of aberrantly methylated CpGs was shown in Case1 (A), Case2 (B) and Case3 (C). Autosomal chromosomal number and sex chromosome are shown on the bottom. D: All chromosomes are aligned side by side in a row (Top), and X chromosome is shown in magnified scale (Bottom). Green bar (hypomethylated CpG in leiomyoma) and red bar (hypermethylated CpG in leiomyoma) are shown in the same line.](pone.0066632.g005){#pone-0066632-g005}

Discussion {#s4}
==========

We performed a genome-wide DNA methylation analysis in uterine leiomyomas using a method that interrogates over 450,000 CpG sites of DNA methylation statuses of consensus coding sequences. Our results suggest that uterine leiomyomas are associated with genome-wide alterations in DNA methylation at multiple gene promoter regions.

Does the Adjacent Myometrium with Leiomyoma have the Potential to Develop into Uterine Leiomyomas? {#s4a}
--------------------------------------------------------------------------------------------------

Uterine leiomyomas often show multifocal tumorigenesis from the myometrium. We speculated that the adjacent myometrium, which looks normal, has the potential to develop into leiomyomas, i.e., already has aberrant DNA methylation. In this study, we selected paired tissue samples of leiomyoma and adjacent myometrium from women who had a single leiomyoma nodule. Interestingly, the DNA methylation status of the myometrium adjacent to the leiomyoma was quite similar to that in the normal uterine myometrium. This result suggests that adjacent myometrium with a leiomyoma nodule does not acquire the potential in DNA methylation levels to develop into leiomyomas. Another possibility is that cells with aberrant DNA methylation are present in the adjacent myometrium, but are too few to detect. In other words, aberrant DNA methylation may occur only in a limited number of cells. Recently, Ono et al. reported that side population cells in the myometrium and leiomyoma tissues act as tissue stem cells and have the potential to differentiate and proliferate [@pone.0066632-Ono1], [@pone.0066632-Ono2]. It has also been shown that each leiomyoma nodule shows monoclonal cell feature [@pone.0066632-Zhang1], indicating that each leiomyoma nodule is derived from just one affected cell as the origin of the tumor. Taken together, these findings suggest that tumor stem cells with aberrant DNA methylation are present in the myometrium and develop into leiomyoma nodules.

Identification of Genes with Aberrant DNA Methylation that are Associated with Tumorigenesis of Uterine Leiomyomas {#s4b}
------------------------------------------------------------------------------------------------------------------

Of 120 genes that had both aberrant DNA methylation and aberrant mRNA expression in leiomyomas, 24 genes were less methylated and transcriptionally upregulated in the leiomyomas compared with the myometrium, and 65 genes were more methylated and transcriptionally downregulated in the leiomyomas compared with the myometrium. The biological relevance of these aberrantly methylated and concomitantly expressed genes was "cancer process" in the section of diseases and disorders according to IPA. The genes associated with "cancer process" contained potentially novel or relevant candidate genes for leiomyoma formation such as IRS1, COL4A1, COL4A2, COL6A3, and GSTM5.

Gene expression of IRS1 was found to be increased in uterine leiomyomas in this study. IRS1 has been shown to be involved in cancer progression and metastasis by mediating proliferative and anti-apoptotic functions of the INSR and IGF1R signaling [@pone.0066632-Esposito1], [@pone.0066632-White1]. Recently, IRS1 has been found to have a novel function in addition to its role in intracellular signal transduction. Interestingly, nuclear presence of IRS1 contributes to malignant transformation [@pone.0066632-Reiss1]. Overexpression or ectopic expression of IRS1 causes cell transformation, including development of the ability to form colonies and tumors [@pone.0066632-DeAngelis1], [@pone.0066632-Valentinis1]. In contrast, when IRS1 expression is experimentally decreased, cancer cells lose their transformed phenotype [@pone.0066632-Cesarone1], [@pone.0066632-delRincon1], [@pone.0066632-DAmbrosio1].

In the present study, the gene expressions of COL4A1, COL4A2 and COL6A3 were found to be increased in uterine leiomyomas. Extracellular matrix (ECM) plays an important role in the pathophysiology of uterine leiomyomas. Collagens are most predominant components of ECM of uterine leiomyomas [@pone.0066632-Stewart1]. Increased production of collagens contributes to the volume expansion and fibroid formation [@pone.0066632-Stewart1]. IRS1 is also reported to be involved in the upregulation of collagen genes, including COL4A1, COL4A2 and COL6A3, suggesting that IRS1 possibly contributes to the leiomyoma nodule formation by upregulating the gene expression of COL4A1, COL4A2, and COL6A3 [@pone.0066632-Tseng1].

On the other hand, gene expression of GSTM5 was decreased in uterine leiomyomas. GST family has been shown to play a crucial role in preventing DNA damage [@pone.0066632-Strange1], [@pone.0066632-Townsend1], and overexpression of GST family inhibits oxidative stress-induced apoptosis and results in cell proliferation [@pone.0066632-Strange1], [@pone.0066632-Townsend1]. Interestingly, GSTM5, one of mu-classes of GST family, may have a different function in tumorigenesis. Peng et al. [@pone.0066632-Missaoui1] reported that promoter DNA hypermethylation and low mRNA expression of GSTM5 was observed in Barrett's adenocarcinoma, and suggested that DNA hypermethylation of GSTM5 is an early event in carcinogenesis [@pone.0066632-Peng1]. Aberrant DNA hypermethylation and downregulation of GSTM5 have also been reported in brain tumors, salivary gland cancers and leukemia [@pone.0066632-Etcheverry1], [@pone.0066632-Bell1], [@pone.0066632-Kearns1]. However, the detailed role of GSTM5 in cancer development remains to be explored.

It is interesting to note that uterine leiomyomas develop only after menarche, indicating that growth of uterine leiomyomas depends on estrogen. We hypothesize that aberrant DNA methylation of the promoter regions of genes targeted by ER alpha causes an abnormal response to estrogen in uterine leiomyomas. Thus, tumor stem cells of leiomyomas in the myometrium may have become hyper-responsive or hypo-responsive to estrogen. In fact, exposure to environmental estrogen during reproductive tract development caused genes targeted by ER alpha to become hyper-responsive to estrogen in the adult myometrium [@pone.0066632-Greathouse1]. The present study detected 22 ER alpha-target genes with both aberrant DNA methylation and aberrant mRNA expression, including COL4A1, COL6A3, DAPK1, and NUAK1. For example, if the promoter is aberrantly hypomethylated in ER alpha-target genes, ER alpha can bind to the promoter, resulting in a hyper-response to estrogen. In the case of the COL4A1 and COL6A3 genes, estrogen exposure causes overproduction of collagens after menarche. On the other hand, if the promoter is aberrantly hypermethylated in the ER alpha-target genes, ER alpha cannot bind to the promoter, resulting in a poor response to estrogen. In the case of DAPK1 and NUAK1, which are apoptosis-related genes that are induced by estrogen, apoptosis is not caused by estrogen, resulting in cell proliferation after menarche. Taken together, these findings suggest that aberrant responses of ER alpha-target genes to estrogen caused by aberrant DNA methylation in their promoter are involved in the development of uterine leiomyomas after menarche.

One may have a question whether the genes with aberrant expression in "cancer process" identified in uterine leiomyomas are also involved in the development of uterine leiomyosarcomas. Skubitz et al. previously reported aberrant overexpression of 16 genes in leiomyosarcomas [@pone.0066632-Skubitz1], but those genes were different from the genes with aberrant expression in this study. Also, uterine leiomyosarcoma development is generally independent on estrogen exposure.

It is reported that tumor stem cells of uterine leiomyomas do not respond to estrogen because of low ER expression [@pone.0066632-Ono2]. It is unclear how tumor stem cells of leiomyomas acquire the responsiveness to estrogen. Interestingly, recent evidence has shown that IRS1 regulates the transcription of ER alpha-target genes by forming a complex with ER alpha [@pone.0066632-Morelli2]. The IRS1/ER-alpha complex translocates into the nucleus and interacts with the promoter region containing ER-responsive elements, which in turn influences transcriptional activities via decreasing proteolytic turnover of ER alpha [@pone.0066632-Morelli2]. IRS1 also enhances ER alpha activities via the PI-3K/Akt pathway [@pone.0066632-Sun1], [@pone.0066632-Campbell1]. From these findings, we speculate that tumor stem cells with overexpression of IRS1 respond to estrogen by acquiring ER alpha expression and stimulating ER alpha activities, resulting in the development of leiomyomas.

Recently, Navarro et al. reported 10 hypomethylated and transcriptionally upregulated genes, and 36 hypermethylated and transcriptionally downregulated genes in uterine leiomyomas compared with the corresponding myometrium using genome-wide approach [@pone.0066632-Navarro1]. However, the aberrant genes reported by them were not consistent with the genes identified in this study. There seem to be a number of differences between two studies. First, the method used for DNA methylation analysis differs. The genome-wide approach used by Navarro et al. covered a few CpG sites in the promoter region in each gene (1--2 CpG sites per gene) while our method can cover much more CpG sites (about 16 CpG sites per gene) broadly from the promoter to the gene body region. Second, the gene extraction algorithm is different between the two studies. We compared the difference in DNA methylation and mRNA expression in a paired sample in each case, and extracted the aberrantly methylated and differentially expressed gene, in which all the three cases showed common aberrations. On the other hand, Navarro et al. compared the difference in DNA methylation and mRNA expression between the uterine leiomyoma group and the myometrium group with the limma package which is based on t-statistic, but did not do in a paired sample in each case [@pone.0066632-Navarro1], [@pone.0066632-Feng1]. Furthermore, we analyzed Japanese women who had a single leiomyoma nodule while Navarro et al. analyzed uterine leiomyomas with multiple nodules from African American women. More importantly, there must be inevitable differences in DNA methylation among individuals, which may be due to the difference in various local and environmental factors among individuals. In fact, we found that the DNA methylation status of the promoter region of ER alpha gene and GS20656 gene (differently methylated genomic loci in uterine leiomyomas) varies in even normal myometrium among individuals [@pone.0066632-Yamagata1], [@pone.0066632-Asada1], which may represent a within-physiological change in a certain cell type such as smooth muscle cells in myometrium [@pone.0066632-Yamagata1], [@pone.0066632-Asada1]. These differences between the two studies and the individual-related difference should be taken into consideration with some caution, and further studies are needed to make a more definitive conclusion.

X Chromosome Genotyping and Chromosome Distribution of Aberrant DNA Methylation {#s4c}
-------------------------------------------------------------------------------

In the present study, we used an advanced genome-wide DNA methylation analysis to confirm that hypomethylated-CpG sites were remarkably enriched on the X chromosome in uterine leiomyomas compared with adjacent myometrium. Loss of the inactive X chromosome and aberrant replication of the active X chromosome has been shown in several female-related cancer cells [@pone.0066632-Kawakami1], [@pone.0066632-Sirchia1]. Our analysis of the X chromosome genotype demonstrated that these events do not occur in uterine leiomyomas. The mechanism of DNA hypomethylation in the X chromosome in uterine leiomyomas is likely to be different from that in other cancer cells [@pone.0066632-Kawakami1], [@pone.0066632-Sirchia1]. We have suggested that X-inactivation machineries are disturbed in uterine leiomyomas [@pone.0066632-Maekawa1]. However, the detailed mechanism of the enriched DNA hypomethylation in X chromosome in uterine leiomyomas, and the relationship between DNA hypomethylation status of X chromosome and the pathogenesis of uterine leiomyomas remain to be elucidated. In the course of this study, we identified 11 genes in the X chromosome of uterine leiomyomas that were hypomethylated in all three cases (data not shown). We are now examining these genes to see if they have roles in the pathogenesis of uterine leiomyomas.

Conclusions {#s4d}
-----------

We propose the following hypothesis for pathogenesis of uterine leiomyomas: tissue stem cells in the myometrium transform to tumor stem cells by suffering genome-wide aberrant DNA methylation by unknown factors such as unfavorable environmental exposures. Aberrant DNA methylation and its related transcriptional aberration in cancer-related genes such as IRS1 may represent a critical initial mechanism that triggers transformation of a single tissue stem cell to a tumor stem cell, which will eventually develop into a monoclonal leiomyoma tumor. In addition, the aberrant DNA methylation of the promoter of ER alpha-target genes, e.g. COL4A1, COL6A3, DAPK1 and NUAK1, is responsible for the aberrant response to estrogen. The tumor stem cells aberrantly respond to estrogen after menarche, and gradually proliferate and differentiate to form leiomyoma nodules.

Our results suggest that aberrant DNA methylation contributes to the pathogenesis of uterine leiomyomas and may lead to the development of new strategies for treatment.

[^1]: **Competing Interests:**The authors have declared that no competing interests exist.

[^2]: Conceived and designed the experiments: RM NS. Performed the experiments: RM. Analyzed the data: RM SS YY ET AN NS. Contributed reagents/materials/analysis tools: SS HA IT LL MO HT. Wrote the paper: RM NS.
