Profiling of Circadian Genes Expressed in the Uterus Endometrial Stromal Cells of Pregnant Rats as Revealed by DNA Microarray Coupled with RNA Interference

The peripheral circadian oscillator plays an essential role in synchronizing local physiology to operate in a circadian manner via regulation of the expression of clock-controlled genes. The present study aimed to evaluate the circadian rhythms of clock genes and clock-controlled genes expressed in the rat uterus endometrial stromal cells (UESCs) during the stage of implantation by a DNA microarray. Of 12,252 genes showing significantly expression, 7,235 genes displayed significant alterations. As revealed by the biological pathway analysis using the database for annotation, visualization, and integrated discovery online annotation software, genes were involved in cell cycle, glutathione metabolism, MAPK signaling pathway, fatty acid metabolism, ubiquitin mediated proteolysis, focal adhesion, and PPAR signaling pathway. The clustering of clock genes were mainly divided into four groups: the first group was Rorα, Timeless, Npas2, Bmal1, Id2, and Cry2; the second group Per1, Per2, Per3, Dec1, Tef, and Dbp; the third group Bmal2, Cry1, E4bp4, Rorβ, and Clock; the fourth group Rev-erbα. Eleven implantation-related genes and 24 placenta formation-related genes displayed significant alterations, suggesting that these genes involved in implantation and placenta formation are controlled under circadian clock. Some candidates as clock-controlled genes were evaluated by using RNA interference to Bmal1 mRNA. Down-regulation of Igf1 gene expression was observed by Bmal1 silencing, whereas the expression of Inhβa was significantly increased. During active oscillation of circadian clock, the apoptosis-related genes Fas and Caspase3 remained no significant changes, but they were significantly increased by knockdown of Bmal1 mRNA. These results indicate that clock-controlled genes are up- or down-regulated in rat UESCs during the stage of decidualization. DNA microarray analysis coupled with RNA interference will be helpful to understand the physiological roles of some oscillating genes in blastocyst implantation and placenta formation.


INTRODUCTION
Circadian rhythms are primarily synchronized with environmental time by the 24-h period of light-dark cycle. In mammals, the master clock in the suprachiasmatic nucleus coordinates the subsidiary oscillators in the majority of peripheral tissues. However, autonomic circadian oscillators are also functional in peripheral tissues. The peripheral circadian oscillators play critical roles in synchronizing local physiology and metabolism to operate in a circadian manner via regulation of the expression of clock-controlled genes (1). These physiological processes mainly include hormonal secretion, gluconeogenesis, lipogenesis, and bile acid homeostasis (2)(3)(4). In addition, cellular differentiation may cause the suspension of the cyclic expression of clock genes (5,6). Recent studies suggest that the circadian system is not only required for proper growth control, but also involved in the circadian regulation of cell proliferation and apoptosis. It has been reported that 2-10% of all mammalian genes are controlled under the circadian clock (1,7,8). Most of these genes are involved in organ functions and show tissue-specific expression. Only a small set of the clock-controlled genes are expressed in multiple organs. Among them are genes that encode key regulators of cell cycle progression (9,10).
Several recent studies have demonstrated that circadian clock genes are rhythmically expressed in the uterus (2,(11)(12)(13)(14)(15). In the uterus composed of heterogeneous cell types, ovarian steroids regulate the proliferation and differentiation of uterus endometrial stromal cells (UESCs). In rodents and humans, the UESCs undergo proliferation and differentiation into decidual cells in response to ovarian steroids and blastocyst implantation at the early stage of pregnancy (16)(17)(18). This process ultimately results in the formation of the placenta.
Mice lacking the Clock gene display abnormal estrus cycles and are infertile (19,20). Furthermore, implantation fails in Bmal1 deficient mice, due to impaired steroidogenesis (21), and mutations of Per1 and Per2 in mice display reproductive deficits in the middle-aged mutant females (22). More recently, we proved the Per2 expression is down-regulated in the UESCs during decidualization that influences the expression of vascular endothelial growth factor A (Vegfa) gene (23). Deregulation of the circadian clock may attenuate or disrupt expression of the clock-controlled genes and can have a profound influence on organ functions. Studies have demonstrated that the circadian clock function is very important for cell cycle, DNA damage response, and tumor suppression in vivo (24)(25)(26).
In the present study, to search the clock-controlled genes expressed during the period of implantation, we analyzed the expression of the clock genes and clock-controlled genes in cultured UESCs prepared from pregnant rats at the stage of implantation using DNA microarray technology. We used transgenic rats constructed with mouse Per2 promoter-destabilized luciferase (Per2-dLuc) reporter gene (27) to precisely adjust the time of gene expression. In addition, several genes of significantly expressed genes including growth factor genes and apoptosisrelated genes were analyzed using RNA interference (siRNA) to Bmal1 mRNA to determine whether these were controlled under circadian clockwork.

ANIMALS
Mouse Per2 promoter region, assembly by NCBI and the Mouse Genome Sequencing Consortium, was fused to a dLuc reporter gene (27). Per2-dLuc transgenic rats were generated in accordance with the method described in the patent publication number WO/2002/081682 (Y.S. New Technology Institute, Utsunomiya, Japan). Adult females were mated with fertile males, and 12:00 p.m.

PREPARATION AND CULTURE OF UESCs
The UESCs were isolated from Per2-dLuc transgenic rats on day 4.50 of gestation as reported previously (6, 28, 29). The harvested cells were washed thrice with fresh DMEM/F12, and seeded onto 35 mm collagen-coated dishes at the density of 2 × 10 5 cells/dish with 2 mL of culture medium (phenol red-free DMEM/F12 supplemented with 10% charcoal-treated FBS and 1× PS). The culture medium was replaced at 15 min after cell seeding to remove epithelial cells. Cells were cultured in a humidified atmosphere of 95% air and 5% CO 2 at 37°C for 2 days. Then, cells were cultured in serumfree medium supplemented with 1× antibiotic-antimycotic (AA; Nacalai Tesque, Kyoto), 1× Insulin-Transferrin-Selenium (ITS, Life Technologies, Grand Island, NY, USA), 0.1% bovine serum albumin (BSA, Sigma Chemicals), and 100 nM progesterone (P 4 , Sigma Chemicals) for additional 2 days prior to other treatments.

REAL-TIME MONITORING OF Per2-dLuc OSCILLATIONS
The cultured UESCs were synchronized with 100 nM dexamethasone for 2 h in the serum-free medium containing 1× AA. Then, cells were given the serum-free medium DMEM/F12 supplemented with 15 mM HEPES, 0.1 mM luciferin (Wako, Tokyo), 0.1% BSA, 1× AA, and 1× ITS, and subjected to luminescence determination. Luciferase activity was chronologically monitored at 37°C with a Kronos Dio AB-2550 luminometer (ATTO, Tokyo) interfaced to a computer for continuous data acquisition, as described previously (6, 13,14). The amplitude and period of Per2-dLuc oscillations were documented by the single Cosinor method using Timing Series Single 6.3 (Expert Soft Tech., Richelieu, France).

MICROARRAY ANALYSIS
RNA samples isolated from cultured UESCs at 30, 36, 42, and 48 h after dexamethasone synchronization were used for microarray analysis using the Whole Rat Genome Microarray 4 × 44 K Ver3.0 (Agilent Technologies, Santa Clara, CA, USA) representing 30,367 probe sets. The preparation of the samples, microarray hybridizations, and bioinformatics analysis were performed by the Cell Innovator at the Kyushu University (Fukuoka, Japan). Bioinformatics analysis was performed using Agilent Future Extraction software (Agilent Technologies). The data were filtered for signal intensity values (p ≤ 0.05, detectable), which allowed removing very low signal values. After this filtering, 12,252 probe sets remained. Probe sets passing the quality control were then analyzed by ANOVA, and p-values of <0.05 were considered significant. The ratio from signal intensity values of four time points was calculated. For the pathway analysis, probe sets were then used for gene ontology analysis using Database for Annotation, Visualization, and Integrated Discovery (DAVID) Bioinformatics Resources (http://david.abcc.ncifcf.gov/) (30). Pathways regulated at p < 0.05 were considered significant.

BMAL1 siRNA TRANSFECTION
Three sequences targeting the Bmal1 mRNA and no silencing RNA for rat were purchased from BOVAC Co. (Kurume, Japan). The sequences of RNA oligos used are listed in Table 1. The scrambled RNA for rat was used as no silencing RNA (BOVAC Co.). Both the Bmal1 siRNA and no silencing RNA were used at final concentrations of 25 nM. The cells were maintained with transfection medium for duration of 12 h (31). Then the medium was replaced with DMEM/F12 supplemented with 1× AA, 1× ITS, 0.1% BSA, and 100 nM P 4 .

RNA EXTRACTION AND RT-qPCR
Cultured cells were harvested at the indicated times, and total RNA was isolated using RNeasy Mini kit (Qiagen), according to the www.frontiersin.org manufacturer's protocol. RNA samples were treated with RNasefree DNase I (Qiagen). The cDNAs were generated by reverse transcription with random primers using a High Capacity Reverse Transcription kit (Applied Biosystems). One microgram of total RNA were reverse transcribed in 20 µL of mixture using MMLV High Performance Reverse Transcriptase (Epicentre Biotechnologies) and Oligo-dT primer according to the manufacture's protocol. Primer sets used for real-time PCR were listed in Table 2.
PCR was performed with a 1:15 dilution of cDNA samples in Master SYBR Green I mixture (Roche Diagnostics) with specific primers (0.25 µM final of each primer) using Mx3000P Real-time QPCR System (Stratagene). Relative quantification of the mRNA levels was performed using the comparative cycle threshold (∆Ct) method. The ∆Ct for each sample was normalized to Gapdh and expressed as relative to the control values (23).

STATISTICAL ANALYSES
All data are presented as means ± SEM of at least three separate experiments, each performed with triplicate samples. The amplitude of Per2-dLuc was documented by Cosinor analysis using Timing Series Single 6.3 (Expert Soft Tech., Richelieu, France). The statistical differences of examined values of target genes in cultured UESCs were determined by Student's t -test using SigmaPlot software (Ver. 11.2; Systat Software, San Jose, CA, USA). Differences were considered significant at p < 0.05 or less.

DNA MICROARRAY ANALYSIS OF OSCILLATING GENES EXPRESSED IN THE RAT UESCs AFTER DEXAMETHASONE SYNCHRONIZATION
During monitoring of Per2-dLuc activity oscillation, RNA samples were prepared at 30, 36, 42, and 48 h after dexamethasone synchronization and global gene expression patterns were determined using DNA microarray technology (Figure 1). The analysis revealed 7,235 significantly altered genes in the UESCs. The increased expression (357 genes) and decreased expression (202 genes) of genes showing with fold change were obtained from 7,235 significantly altered genes at four time points during oscillation ( Table 3). The majority of fold-changed genes were identified in both RNA samples 30 versus 48 h and 36 versus 48 h, including growth factors, transcription factors, receptors, channels, and enzymes ( Table 4). For example, the expression of growth factors such as Inhβb, Gdf10, and Gdf15 were increased during the interval, whereas the expression of Igf1 was decreased. To place the microarray results in the cellular context, we performed the biological pathway analysis using the DAVID online annotation software. This analysis revealed that genes were involved in cell cycle, glutathione metabolism, MAPK signaling pathway, fatty acid metabolism, ubiquitin mediated proteolysis, focal adhesion, PPAR signaling pathway, and so on ( Table 5). It was identified that the expression of 38 cell cycle related genes were significantly altered during the second to third phase of Per2-dLuc activity oscillation (p < 0.00001).

FUNCTIONAL MOLECULAR MACHINERY OF A CIRCADIAN CLOCK EXISTS IN THE RAT UESCs AFTER DEXAMETHASONE SYNCHRONIZATION
In order to determine whether 18 clock genes displayed diurnal rhythmic expression, we performed gene clustering on the microarray results. The gene expression profiles were mainly divided into four groups: the first group was Rorα, Timeless, Npas2, Bmal1 (Arntl), Id2, and Cry2; the second group Per1, Per2, Per3, Dec1 (Bhlhe40), Tef, and Dbp; the third group Bmal2 (Arntl2),   Figure S1 in Supplementary Material). The Per2 intensity profile was consistent with Per2 circadian oscillation (Figure 2). Profiles of Bmal1 and Rev-erbα expression were obviously antiphase in the UESCs, because Bmal1 promotes Rev-erbα transcription with Clock and in turn Rev-erbα inhibits Bmal1 transcription through its binding to the RORE in the promoter region. The Per2 and Dbp expression profiles were also out-of-phase with the Bmal1 rhythm. The protein kinase genes such as Csnk1ε and Csnk1δ also displayed the weak rhythmic expression (Figure 2). The expression profiles of core clock genes such as Bmal1, Clock, Per1, Per2, and Rev-erbα were further confirmed by RT-qPCR. As the results are shown in Figure 3, their profiles were very similar to the microarray results. The Cosinor analysis method was used to determine the rhythmic expression of examined genes. Bmal1, Per1, Per2, and Rev-erbα displayed rhythmic expression (p < 0.05).

SIGNIFICANT ALTERATIONS OF GENES RELATED TO EITHER IMPLANTATION OR PLACENTA FORMATION
In order to determine whether genes involved in implantation or placenta formation displayed significant alteration, we performed gene clustering on the microarray results. As the gene clustering is shown in Figure S2 in Supplementary Material, 11 genes of 33 implantation-related genes displayed significant alterations (p < 0.05). Of 11 significantly altered genes, three genes such as Bsg, Pcsk5, and Klf9 were up-regulated, whereas Tgfβr2, Pparδ, Ptgs2, Grn, and Fkbp4 were significantly down-regulated. The results of gene clustering for placenta formation are shown in Figure S3 in Supplementary Material. Twenty-four genes of 107 placenta formation-related genes displayed significant alterations (p < 0.05). Of 24 significantly altered genes, nine genes Txnrd1, Pgf, Pkd2, Hif1a, Ccnf, Plac9, Prdx3, Adm, and Pbrm1 were upregulated, whereas Plac8, Prdm1, Rspo3, Tmed2, Dcn, Epas1, and Stk3 were significantly down-regulated. These results indicate that these genes involved in implantation and placenta formation may be controlled under circadian clock.

UP-REGULATION AND DOWN-REGULATION OF GENES RELATED TO CELL GROWTH AND APOPTOSIS
Of fold-changed genes ( Table 4), we focused on Inhβb and Igf1 that were regulated with opposite directionality. Both Inhβa and Inhβb were high in signal intensity, whereas Inha was low, indicating that activin, but not inhibin, is highly produced in the UESCs. The Fst gene coding the activin-binding protein follistatin was down-regulated. In contrast to Inhβa and Inhβb, Igf1 was down-regulated, whereas its binding proteins Igfbp3 and Ifgbp4 were rather increased (Figure 4). Consequently, the activities of activin and IGF1 may be differentially regulated during circadian rhythms. During active oscillation of circadian clock, however, the apoptosis-related genes Fas and Caspase3 remained no significant changes (Figure 4).

EFFECT OF Bmal1 KNOCKDOWN ON THE EXPRESSION OF GENES RELATED TO CELL GROWTH AND APOPTOSIS IN THE RAT UESCs
In order to suppress the cellular clock in the UESCs, Per2-dLuc oscillations were investigated using Bmal1 siRNA (siRNA) or no silencing RNA (CONT). The UESCs transfected with either RNA displayed several Per2-dLuc oscillations. A decline of Per2-dLuc bioluminescence oscillation and significantly decreased amplitude (p < 0.01, versus CONT) were observed in Bmal1 siRNA treated cells (Figure 5A). There was not a significant effect of Bmal1 siRNA treatment on the peak time of the first phase and the cycle time (data not shown). The transcript levels of core clock genes were estimated at the peak time of the first phase. The results are shown in Figure 5B. Bmal1 mRNA expression was inhibited after transfection with Bmal1 siRNA, the Bmal1 transcript level being reduced by 60% compared with that in the CONT group. Rev-erbα and Dbp are regarded as core circadian oscillator genes that fall under the regulation of Bmal1-Clock heterodimer www.frontiersin.org through E-box elements located in their promoter regions. Thus Rev-erbα and Dbp were significantly down-regulated in Bmal1 siRNA cells.
In addition, we observed that the mRNA expression of several uterus genes (Igf1, Inhβa, Fas, and Caspase3) was differentially affected by the Bmal1 siRNA treatment. The expression of Igf1 gene was significantly down-regulated by Bmal1 silencing (p < 0.05), whereas the Inhβa was up-regulated approximately sevenfold (p < 0.001) (Figure 5C). Although Fas, and Caspase3 did not exhibit cellular circadian rhythms as described above, these genes were significantly up-regulated by Bmal1 silencing (p < 0.01).

DISCUSSION
The present analysis of circadian genes regulated by clockwork system significantly contributes to our understanding of the relationship of cellular circadian clock with implantation and decidualization. Our present study on cultured UESCs of pregnant rats is the first to analyze the global gene expression profiles using DNA microarray technology. We found the specific expression profiles of 18 clock genes and two protein kinase genes. Bmal1 (Arntl) and Rev-erbα (Nr1d1) displayed typically anti-phase expression patterns, showing that the proper range of clock gene oscillation was successfully covered. We identified 7,235 significantly altered genes in the UESCs and 559-fold-changed genes with up or down directionality. The biological pathway analysis using the DAVID online annotation software revealed that many genes involved in cell cycle (38 genes), glutathione metabolism (19 genes), MAPK signaling pathway (55 genes), and so on were significantly expressed with circadian regulation. In addition, some genes involved in implantation or placenta formation displayed significant alteration with up or down directionality: 11 genes involved in implantation and 24 genes involved in placenta formation.

FIGURE 3 | Expression profiles of core clock genes in the uterus endometrial stromal cells after dexamethasone synchroni zation as revealed by RT-qPCR.
According to the Per2-dLuc oscillation as illustrated in Figure 1, total RNA samples were collected at the indicated times after synchronization. qRT-PCR analyses of transcript levels were performed using their specific primers. Gapdh was used as an internal control. Each value represents the mean ± SEM (n = 3). in the decidual stromal cells and trophoblasts (34). GPR4, a proton-sensing G-protein-coupled receptor, is also expressed in human umbilical vein endothelial cells (35). GPR4 signaling regulates endothelial cell adhesion through the cAMP pathway. It is known that the oxidized low density lipoprotein receptor 1 (Olr1) regulates growth of a variety of cells and is important in inflammation, oxidative stress, and tissue remodeling (36). However, the expression and regulation of Olr1 remain unclear in the uterus. Aquaporin-1 (Aqp1) is proposed as a mediator of estrogen-induced angiogenesis in breast cancer and endometrial cancer (37). The expression of Aqp1 was reported in the pregnant rat uterus (38, 39). Although more than 10 Aqp1 have been identified up to date, only Aqp1 was expressed in the cultured UESCs.
In contrast, the down-regulated representative genes with fold change were also listed, including Fst, Igf1, Zmiz1, Fbxo9, Rgs2, Rgs5, Mt1a, and Mt2A. IGF family genes including Igf1 are differentially expressed throughout gestation, and especially the IGF system contributes phenotypic and functional changes of myometrium smooth muscle cells (40). The overexpression of Zmiz1, recently identified as a candidate oncogene, was reported in human breast, ovarian, and colon cancers (41).
Zmiz1 was also reported to function in regulating the activity of the androgen receptor as a co-activator (42). However, the expression and function of Zmiz1 are unclear in the uterus. Fbxo9 was reported to contribute survival of myeloma cells (43). Rgs2 and Rgs5, regulators of G-protein signaling proteins, are abundantly expressed in pregnant human myometrium (44). Rgs2 was reported to up-regulate in response to stress (45). The expression of Mt1a and Mt2A, typical metal response parameters, was reported in the proliferative phase of endometrium (46).
Of up-regulated and down-regulated genes, we focused on the expression of Inhβb, Inhβa, and Igf1 that displayed fold change with up or down directionality. Both the Inhβa and Inhβb genes showed high signal intensities of DNA microarray, whereas Inhα was low. This result indicates the expression of activin. The binding protein of activin, follistatin (Fst ), was downregulated, indicating the relative action of activin is increased. Expression of activin A was increased with reduced stromal cell mitosis, tissue growth, and mitogenic signaling in the decidual endometrium (47). We analyzed using Bmal1 siRNA whether the expression of Inhβa and Igf1 was controlled under circadian clock. Transfection of Bmal1 siRNA into the UESCs induced Frontiers in Endocrinology | Systems and Translational Endocrinology a significant decline of Per2-dLuc bioluminescence oscillation. However, the transfection did not alter other parameters in the oscillation. The core clock genes Rev-erbα and Dbp as well as Bmal1 were down-regulated. The expression of Igf1 was significantly decreased after Bmal1 silencing, suggesting that Igf1 is positively regulated by circadian clock. Conversely, the expression of Inhβa was enhanced approximately sevenfold after Bmal1 silencing, suggesting that Inhβa is negatively regulated by circadian clock. Although both Inhβa and Igf1 are clock-controlled genes, thus, the two gene are differentially regulated. Rev-erbα is a critical component of circadian clock (48). It acts as a domain clock-controlled gene under the regulation of Bmal1-Clock heterodimer; in turn it has been shown to suppress Bmal1 transcription through binding to RORα response elements (RORE) presenting in Bmal1 promoters (49). Recent studies have shown that Rev-erbα is a transcriptional silencer (49,50). However, Reverbα increases the transcription of Star by binding to its agonist heme (51).
In the uterus endometrial cells of pregnant rats, the early proliferative phase is characterized by tissue remodeling, angiogenesis, and modulation of inflammation; the mid-proliferative phase is characterized not only by proliferation but also marks the onset of expression of genes required for endometrial receptivity. Thus, cell growth and apoptosis support the process of tissue remodeling and implantation. The death receptor Fas/ligand system is a key regulator of apoptotic cell death and irregularity of this signaling pathway has been shown to participate in immune-mediated βcell apoptosis (52,53). It is also known that the expression level of the Fas gene is down-regulated in a variety of malignancies including malignant melanoma, adenocarcinoma, and squamous cell carcinoma (54,55). In the cultured UESCs, however, the expression of Fas and Caspase3 displayed no or small changes. Interestingly, these genes were significantly up-regulated after Bmal1 silencing, similar to Inhβa, suggesting that these apoptosisrelated genes are suppressed during active circadian rhythms of clock genes. It is possible that dysfunction of the circadian clockwork during the stage of decidualization is necessary to increase or decrease expression of clock-controlled genes for formation of the placenta.

ACKNOWLEDGMENTS
This work was supported by a Grant-in-Aid for Scientific Research (B) from the Japan Society for the Promotion of Sciences (JSPS; 22380152) (to Masa-aki Hattori).