Exploration of Global Gene Expression Changes During the Estrous Cycle in Equine Endometrium¹

Abstract

The equine endometrium exhibits characteristic morphological and functional changes during the estrous cycle controlled by the interplay of progesterone and estradiol. A microarray analysis of endometrial tissue samples derived from five time points of the estrous cycle (Day [D] 0, D3, D8, D12, and D16) was performed to study the dynamics of equine endometrial gene expression. Statistical analysis revealed 4996 genes differentially expressed during the estrous cycle. Clustering of similar expression profiles was performed to find groups of coregulated genes. This revealed eight major profiles: highest mRNA concentrations on D0, from D0 to D3, on D3, from D3 to D8, on D8, from D8 to D12, from D12 to D16, and on D16. Bioinformatics analysis revealed distinct molecular functions and biological processes for the individual expression profiles characterizing the different phases of the estrous cycle (e.g., extracellular matrix and inflammatory response during the estrus phase, cell division and cell cycle during early luteal phase, and endoplasmic reticulum, protein transport, and lipid metabolism in the luteal phase). A comparison to dynamic gene expression changes in bovine endometrium identified common and species-specific gene regulations in cyclic endometrium. Analysis of expression changes during the estrous cycle for genes previously found to be differentially expressed on D12 of pregnancy provided new evidence for possible regulation of these genes. This study provides new insights regarding global changes of equine endometrial gene expression as molecular reflections of physiological changes in the cyclic equine endometrium with regard to the crucial role of this tissue for successful reproduction.

Introduction

Throughout the estrous cycle, the equine endometrium exhibits characteristic morphological and functional changes, which are mainly induced by the ovarian hormones progesterone (P4) and estradiol (E2) [1]. The estrous cycle of the mare is approximately 22 days in length, with seasonal activity related to long days. According to the changes in steroid hormone concentrations and corresponding endometrial structural and functional changes, the equine estrous cycle can be divided into four different stages: estrous/ovulatory phase (Day 16 to Day 0), early luteal phase (Day 1 to Day 4), midluteal phase (around Day 8) and late luteal phase (Day 12 to Day 15) (for review, see [1]). During these stages, the endometrium undergoes pronounced changes. At estrus, endometrial edema and intense secretory activity (merocrine activity and exocytosis) can be observed [2]. The increased secretory activity during estrus, together with infiltration of the luminal epithelium by granulocytes and macrophages, is important for uterine clearance after breeding [2]. In the ovary, follicular maturation is controlled by the follicle-stimulating hormone and the luteinizing hormone (LH) from the pituitary, which are stimulated by gonadotropin-releasing hormone, a tropic peptide hormone, which in turn is synthesized and released from neurons within the hypothalamus [3]. As in ruminants, the high level of LH during this stage primarily stimulates the further development of the dominant follicle, producing high amounts of E2 responsible for the LH surge in a positive-feedback loop preceding ovulation. In contrast to other large domestic animals, no distinct periovulatory LH peak exists in the mare, and the period of elevated LH concentrations lasts for several days [4]. Following ovulation, P4 levels immediately start to increase in the course of the formation of the new corpus luteum (CL), whereas E2 decreases to basal levels [5]. The height of epithelial cells changes, the edema of lamina propria declines, and uterine glands widen at this time point of the equine estrous cycle [6].

In the mare, oxytocin has been shown to be produced by hypothalamic neurons and secreted from the posterior pituitary, but it is also produced in the endometrium [7]. Oxytocin and its endometrial receptor (OXTR) play an important role in stimulating the pulsatile release of prostaglandin (PG) F2α from the endometrium [8, 9]. This pulsatile uterine secretion of PGF2α is responsible for luteolysis starting around Day 14 after ovulation [10]. In contrast to ruminants, a systemic utero-ovarian pathway exists in the horse; that is, PGF2α acts via the peripheral circulation on the CL [11]. So far, expression of a number of genes involved in PG synthesis and signaling has been investigated during the estrous cycle and early pregnancy [12–14]. Altogether, increased capacity of PG production during the late luteal phase has been observed in cyclic mares (e.g., increased expression of PG-endoperoxide synthase 2 [PTGS2], also known as cyclooxygenase 2 [COX2], in uterine epithelial cells on Day 15 of the cycle), whereas increase of endometrial PTGS2 expression is inhibited in pregnant mares [13].

In addition to genes involved in PG metabolism and signaling, a number of other genes involved in regulation of cyclic changes in the endometrium have been investigated, such as estrogen and P4 receptor genes [15], genes related to endometrial vascular changes [16], and wingless-type (WNT) genes and their antagonists [17]. To get a holistic overview of dynamic gene expression changes during the estrous cycle in equine endometrium, a DNA microarray analysis of endometrial biopsy samples was performed, which were collected at five defined time points of the estrous cycle: Day 16 (beginning of follicular phase), Day 0 (day of ovulation), Day 3 (early luteal phase), Day 8 (midluteal phase), and Day 12 (late luteal phase). The obtained differentially expressed genes (DEGs) were analyzed using bioinformatics tools to describe biological processes and molecular pathways involved in the regulation of endometrial functions during the estrous cycle in the mare.

Materials and Methods

Collection of Endometrial Tissue Samples

The mares used for the present study were provided by and housed in the facilities of the Bavarian Principal and State Stud of Schwaiganger, Germany. Treatments of mares were approved and performed with permission of the local authorities (District Government of Upper Bavaria). Collection of biopsy and blood samples was performed in accordance with the International Guiding Principles for Biomedical Research Involving Animals, as proposed by the Society for the Study of Reproduction, and with the European Convention on Animal Experimentation and the German Animal Welfare Act.

Endometrial biopsy samples were collected from five normal-cycling mares (Bavarian Warmblood). Follicular development and ovulation were monitored routinely by transrectal palpation and ultrasound examination. Samples were obtained by transcervical biopsy and taken at five different time points during the estrous cycle (Days 0, 3, 8, 12, and 16), where Day 0 corresponds to the day of ovulation. Biopsy samples were collected during different cycles. When two biopsy samples were collected in the same cycle, the interval between biopsy collections was at least 9 days. For analysis of tissue composition, the biopsy samples were cut transversely into six equal and plane-parallel slices. Every second slice was transferred into embedding capsules with the right cut surface facing downward, covered with a foam sponge to avoid tissue sample distortion, and fixed by immersion in 4% buffered formaldehyde. The remaining pieces of the biopsy samples were immediately transferred into vials containing 3 ml of RNAlater (Life Technologies GmbH) for use in mRNA expression analysis. The vials were cooled on ice and incubated overnight at 4°C. Samples were stored at −80°C until further processing. Additionally, blood samples were collected in ethylenediaminetetra-acetic acid tubes from the jugular vein at the same time at biopsy collection for determination of peripheral plasma P4 concentrations. Blood samples were centrifuged at 2000 × g at 4°C for 10 min, and plasma was decanted and stored at −20°C until assay (mini VIDAS and VIDAS Progesterone kits; bioMérieux Deutschland GmbH).

Quantitative Stereological Analysis

For qualitative histological and quantitative stereological analyses, formalin-fixed slices of each biopsy sample were routinely processed and embedded in paraffin with the right cut surface facing downward. Histological sections were cut at a nominal thickness of 3 μm with a rotary microtome, transferred onto glass slides, and stained with hematoxylin and eosin. Quantitative stereological analyses were carried out with the newCAST software (Visiopharm A/S). Slides were displayed on a monitor at a 400× final magnification via a camera (Universal Camera DP72; Olympus Deutschland GmbH) coupled to a microscope (BX41 Laboratory Microscope; Olympus Deutschland GmbH), and images were superimposed by an adjustable point counting grid. Approximately 2000 points were evaluated per biopsy sample to determine volume densities of surface epithelium, glandular epithelium, blood vessels, and remaining tissue. The volume densities (Vv) of the different tissue compartments were obtained by dividing the number of points hitting a compartment (P(compartment), such as points hitting blood vessels, or P(blood vessels)) by the total number of points hitting the biopsy sample (P(sample)) considering the number of points in the counting grid used: Vv(compartment/sample) = P(compartment)/P(sample) × points used in grid.

P4 Assay

Progesterone concentrations in peripheral blood plasma were measured with mini VIDAS and VIDAS Progesterone kits, a system based on the enzyme-linked fluorescent assay technique. A detection limit of 0.25 ng/ml and a correlation coefficient of 0.89 toward radioimmunoassay are certified for the assay by the manufacturer.

Microarray Analysis

Total mRNA was isolated from the 25 endometrial biopsy samples using TRIzol Reagent (Life Technologies GmbH) according to the manufacturer's instructions. Quantity and purity of RNA were measured with a NanoDrop 1000 (PEQLAB Biotechnologie GmbH). Quality of total RNA was determined electrophoretically with an 2100 Bioanalyzer (Agilent Technologies). RNA integrity values ranged from 7.1 to 9.2. Microarray analysis was performed using 4×44k Horse Gene Expression Microarray (Agilent Biotechnologies). Cy3-labeled cRNA was produced with the Quick Amp Labeling Kit, One-Color (Agilent Technologies) and hybridized to the microarray according the manufacturer's instructions. Hybridized and washed slides were scanned at a resolution of 3 μm with a G2505C DNA Microarray Scanner (Agilent Technologies). Image processing was performed with Feature Extraction Software (Version 10.7.3.1; Agilent Technologies). Processed signals were filtered based on “Is well above background” flags (detection in four of five samples in at least one of the five experimental groups) and subsequently normalized with the BioConductor package vsn. A heatmap based on pairwise distances (BioConductor package geneplotter) was generated for quality control of 25 samples. Significance analysis was performed using the significance analysis of microarrays function of the BioConductor siggenes package (multiclass, false-discovery rate [FDR] = 1%) The data discussed in this publication have been deposited in the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO; http://www.ncbi.nlm.nih.gov/geo/) and are accessible through GEO Series accession number GSE39043.

Functional Analysis of Array Data

The Agilent Horse Microarray was reannotated based on Ensembl, Entrez Gene, NCBI BLAST, and University of California, Santa Cruz BLAT analyses to obtain equine and human (putative orthologous genes) Entrez Gene identifiers and the corresponding gene information. To classify the DEGs according to their gene expression profile during the estrous cycle, clustering with the self-organizing tree algorithm (SOTA; MeV 4.7.1, TM4 software suite; http://www.tm4.org/mev/) was performed. Functional classification of the DEGs was done for each expression cluster with the functional annotation clustering tool of the Database for Annotation, Visualization, and Integrated Discovery (DAVID) [18]. This analysis was performed on the basis of Entrez Gene IDs of the putative human orthologous genes.

Quantitative Real-Time RT-PCR

The same RNA samples used for microarray analysis were also used for quantitative real-time RT-PCR (qPCR). First-strand cDNA was synthesized starting from 1 μg of total RNA with the Sprint RT Complete-Double PrePrimed kit (Takara Bio Europe/Clontech). The qPCR experiments were performed in accordance with the StepOne Software (Version 2.1; Applied Biosystems) protocol with SYBR Green reagent. Primer sequences, annealing temperatures, and melting points are shown in Supplemental Table S1 (all Supplemental Data are available online at www.biolreprod.org). The threshold cycle (Ct), which inversely correlates with the target mRNA level, was measured as the cycle number at which the SYBR Green fluorescence signal appeared above the background threshold by StepOne Software. Ct-values of the target genes were normalized against the reference gene ubiquitin B (UBB). Additional reference genes were tested but did not show stable expression over all five time points. All PCRs were performed in duplicates. Amplified PCR fragments were sequenced to verify the resulting PCR product.

Results

Quantitative Stereological Analysis

The volume percentage of luminal epithelium, blood vessels, glandular epithelium, and the remaining tissue was estimated in the endometrial biopsy samples using quantitative stereological techniques. The estimated volume fractions ranged during the estrous cycle from 0.5% to 1.1% for luminal epithelium, from 1.4% to 2.4% for blood vessels, from 30.1% to 34.4% for glandular epithelium, and from 63.5% to 68.1% for the remaining tissue (Supplemental Fig. S1). For all analyzed compartments, no statistically significant difference was found during the estrous cycle, but a trend (P = 0.084) for higher volume percentage on Days 12 and 16 was observed for glandular epithelium (Supplemental Table S2). Overall, tissue composition was quite consistent within the biopsy samples for the percentage of glandular epithelium and remaining tissue. Relatively high deviations were obtained for the estimated proportions of luminal epithelium and blood vessels.

Peripheral Plasma P4 Concentrations

Blood samples for measurement of plasma P4 concentrations were collected in different cycles (when the corresponding biopsy sample was taken), except for samples from Day 0. P4 concentrations were similar between all five mares at each time point of the estrous cycle (Supplemental Fig. S2). P4 concentrations were low on Day 16 (1.2–3.7 ng/ml), basal on Day 0 (0.8–1.4 ng/ml), rising on Day 3 (4.1–11.2 ng/ml), highest on Day 8 (14.8–30.0 ng/ml), and on average, second highest on Day 12 (11.0–26.8 ng/ml). The SD of P4 concentrations was highest on Days 8 and 12 (5.5 and 6.5, respectively).

Microarray Analysis

After data processing and normalization, the microarray data set was initially analyzed with a heatmap based on pairwise correlations (Fig. 1). Samples collected on Day 8 as well as samples from Day 12 clustered together in separate groups clearly distant from the samples collected on Days 16, 0, and 3. Also, the samples from Day 0 and Day 3 formed relatively homogeneous groups. Samples from Day 16 were the most heterogeneous group, where four samples clustered together and one sample (D16M1) was more similar to samples from Day 0.

Fig. 1

Heatmap of pairwise correlations of microarray data sets. Microarray data of the 25 samples were filtered for detectable probes and normalized with the BioConductor package vsn. Normalized data were used for calculation of pairwise distances and drawing of a heatmap by use of the BioConductor package geneplotter. Each column represents one sample and shows the correlation to all samples (including itself), with red for correlation = 1 and blue for the lowest observed correlation.

Statistical analysis of microarray data revealed almost 10 000 probes with significant differences in signal intensity during the estrous cycle (BioConductor package siggenes, multiclass; FDR = 1%). These microarray probes represented 4996 different genes/transcripts (Supplemental Table S3). To classify these genes according to their gene expression during the estrous cycle, a clustering based on similarity in expression profiles (SOTA) was performed. This analysis revealed eight major expression profiles (Fig. 2 and Supplemental Table S3). The genes in cluster 1 showed, on average, highest expression levels at Day 0 and low levels during the luteal phase. Expression of some of these genes had already increased on Day 16. Cluster 2 contained mRNA expression profiles with highest levels from Day 0 to Day 3 and low levels from Day 8 to Day 16. Expression of the majority of genes in cluster 3 peaked on Day 3. For the genes in cluster 4, highest mRNA levels were found from Day 3 to Day 8. Genes of cluster 5 showed profiles similar to those of cluster 4 but with later increase and higher levels up to Day 12. Cluster 6 displayed higher mRNA levels from Day 8 to Day 12 and lowest levels on the remaining days of the estrous cycle. The highest expression levels in cluster 7 were found on Day 12, with a slight decrease to Day 16 and lower levels on Days 0, 3, and 8. Cluster 8 comprised genes with highest mRNA levels on Day 16 and lowest levels on Days 3 and 8.

Fig. 2

Clustering based on similarity in gene expression profiles during the estrous cycle. Relative gene expression levels during the estrous cycle were obtained for all DEGs by calculation of mean-centered expression values (vsn-normalized value of one time point minus mean of vsn-normalized values of all five time points; vsn-normalized values are in log2 scale). SOTA of MeV software (Version 4.7.1) was used to obtain groups of genes with similar expression profiles. Vertical axis is in log2 scale and shows the difference to the mean expression during the estrous cycle. Horizontal axis shows the day of the estrous cycle. The SOTA dendrogram (bottom right) shows the similarity between clusters and mean expression changes during the cycle.

Microarray results were confirmed with analysis of selected genes by qPCR. Genes were selected as representatives from the obtained gene expression clusters (SOTA cluster). Overall, expression differences found by microarray analysis were clearly confirmed (Fig. 3 and Supplemental Table S4).

Fig. 3

Confirmation of microarray data. Quantitative real-time RT-PCR was used to confirm microarray data for selected genes. To illustrate correlation of microarray and qPCR data, relative expression levels were calculated for microarray data (M; vsn-normalized value of one time point minus mean vsn-normalized values of all five time points) and qPCR data (P; Ct-value of one time point minus mean Ct-values of all five time points) for each gene and are shown as a heatmap (MeV software, Version 4.7.3). Color scale is from 4-fold lower (log2-fold = −2) than the mean (blue) to 4-fold higher (log2-fold = 2) than the mean (red). D, day of the estrous cycle.

Bioinformatics Analysis of Microarray Data

To characterize molecular functions and assigned biological processes overrepresented during different phases of the estrous cycle, functional annotation clustering (DAVID) was performed for the individual clusters of genes with similar expression profiles obtained by SOTA clustering (Fig. 2). This analysis resulted in a relatively large number of significant annotation clusters of related functional terms that represented distinct overrepresented functions/processes at different stages of the estrous cycle (Supplemental Table S5). Table 1 shows the most representative functional terms from enriched DAVID annotation clusters for the eight major expression profiles obtained by SOTA clustering. For graphical illustration of the processes overrepresented during distinct phases of the estrous cycle, a heatmap based on the frequency patterns (numbers of assigned genes at different cycle stages) during the cycle was plotted (Fig. 4). The selected functional categories were clustered according to their frequency patterns during the estrous cycle. The bottom half of Figure 4 shows categories with highest numbers of genes in clusters 1 and 8 (highest levels on Days 16 and 0). These functional groups represented processes related to extracellular matrix (ECM) synthesis, focal adhesion, vasculature development, and immune response. ECM-related categories had more genes on Day 0 and immune-related categories contained more genes on Day 16. For the category translational elongation and the category ribosome, highest numbers of genes were found for clusters 1 and 2 (genes with highest expression levels from Day 0 to Day 3). Specifically for cluster 3 (highest expression levels on Day 3), categories related to cell division showed highest numbers of genes. The category protein disulfide isomerase activity likewise showed a distinctive pattern and was represented only in SOTA clusters 4 and 5 (n = 3 and 4 genes, respectively). Functional categories with highest numbers of genes during the luteal phase (SOTA clusters 4, 5, and 6) were related to a variety of functions and processes, such as metabolic processes (lipids, fatty acid, and glucose), secretory processes (vesicles, endoplasmic reticulum, and Golgi apparatus), and cell junctions. Some functional categories showed a kind of biphasic pattern of their gene frequencies during the estrous cycle, such as angiogenesis and blood circulation, with highest numbers of genes on Day 0 and on Day 8 or 12, respectively (Fig. 4).

Fig. 4

Representation of selected functional categories/terms during the estrous cycle. Functional terms overrepresented during distinct phases of the estrous cycle were selected from the results of DAVID functional annotation clustering analysis. Numbers of genes assigned to these functional terms (in rows) are shown for each of the eight SOTA clusters (in columns) as a percentage of the number of all DEGs obtained for a term. Color scale is from 0% (white) to 60% (brown).

Table 1

Representative functional terms from enriched DAVID annotation clusters.

Expression cluster	Representative enriched functional termsa	Enrichment scoreb
Cluster 1: Highest level at Day 0; 636 genes; 621 human IDs	Extracellular region part (88, 2.6); extracellular matrix (53, 4.4)	17.05
	Collagen (24, 8.0); focal adhesion (24, 2.9); ECM-receptor interaction (20, 5.8)	11.38
	Glycoprotein (204, 1.5); signal peptide (177, 1.7); extracellular region (125, 1.8)	10.82
	Translational elongation (30, 8.7); ribosome (29, 3.8)	8.81
	Cytokine-cytokine receptor interaction (17, 1.6); chemotaxis (16, 2.9); cytokine activity (11, 1.7)	2.37
	Vasculature development (21, 2.4); blood vessel morphogenesis (13, 1.8); angiogenesis (11, 2.2)	2.32
Cluster 2: Highest level from Day 0 to Day 3; 583 genes; 563 human IDs	Ribonucleoprotein complex (52, 3.2); ribosome (33, 4.9); translational elongation (28, 9.0)	11.80
	Nuclear lumen (77, 1.7); nucleoplasm (52, 1.9)	5.43
	RNA processing (34, 2.0); RNA splicing (21, 2.4); spliceosome (15, 3.5)	3.10
	RNA transport (11, 3.7); RNA export from nucleus (5; 3.9)	2.84
	Cell cycle (49, 2.0); cell division (19, 2.1); mitosis (16, 2.3)	2.79
Cluster 3: Highest level at Day 3; 721 genes; 713 human IDs	Cell cycle (134, 4.1); mitosis (70, 7.6); cell division (67, 5.5)	42.47
	Nuclear lumen (168, 2.8)	38.08
	Chromosome (86, 4.5); kinetochore (32, 13.8)	27.09
	Nucleotide binding (169, 1.9); ATP binding (114, 2.0)	11.92
	Ribonucleoprotein complex biogenesis (39, 5.2); noncoding RNA processing (33, 4.2); ribosome biogenesis (29, 5.7)	10.97
	Cytoskeleton (81, 1.4); microtubule cytoskeleton (67, 2.9); spindle (38, 6.2)	9.53
	RNA processing (74, 3.3); RNA splicing (36, 3.0); spliceosome (19, 3.5)	8.81
	Mitochondrion (89, 2.0); mitochondrial matrix (24, 2.5)	5.74
	Lipid biosynthetic process (30, 2.2); steroid biosynthetic process (14, 4.0); isoprenoid biosynthetic process (5, 6.0)	3.36
Cluster 4: Highest level from Day 3 to Day 8; 501 genes; 495 human IDs	Protein transport (44, 2.2); intracellular transport (36, 2.1)	4.27
	Monosaccharide metabolic process (19, 3.3); glucose metabolic process (12, 3.0)	3.68
	Mitochondrion (53, 1.6); mitochondrial matrix (18, 2.6)	3.22
	Endoplasmic reticulum lumen (9, 3.8)	2.88
	Endoplasmic reticulum part (26, 2.5); endoplasmic reticulum membrane (17, 2.1)	2.80
	Nucleotide binding (87, 1.4); ATP binding (53, 1.3)	2.26
	Cytoplasmic vesicle (34, 1.8); membrane-bounded vesicle (27, 1.6)	2.20
Cluster 5: Highest level at Day 8; 701 genes; 692 human IDs	Electron carrier activity (15, 1.8); flavoprotein (13, 3.5)	2.79
	Protein transport (50, 1.8); intracellular transport (39, 1.6)	2.72
	Nucleotide binding (107, 1.3); ATP binding (74, 1.4)	2.46
	Endomembrane system (43, 1.4); endoplasmic reticulum part (28, 2.1)	2.40
	Endoplasmic reticulum lumen (10, 3.2); thioredoxin domain (7, 6.9); protein disulfide isomerase activity (4, 12.0)	2.36
	Lipid metabolism (14, 2.6); fatty acid metabolic process (13, 1.8)	2.21
	Monosaccharide metabolic process (16, 2.0); glucose metabolic process (12, 2.1)	1.97
	Kinase (45, 1.9); protein serine/threonine kinase activity (25, 1.6); h_mapkPathway:MAPKinase Signaling Pathway (12, 3.2)	1.74
Cluster 6: Highest level from Day 8 to Day 12; 709 genes; 683 human IDs	Cytoplasmic vesicle (49, 2.0)	5.40
	Golgi apparatus (58, 1.8)	4.58
	Vacuole (19, 2.0); lysosome (16, 2.0)	2.32
	GTPase regulator activity (21, 1.5); regulation of Ras GTPase activity (9, 2.4); Rab GTPase activator activity (8, 4.7)	1.85
	Regulation of apoptosis (39, 1.4); negative regulation of apoptosis (23, 1.8)	1.76
Cluster 7: Highest level from Day 12 to Day 16; 501 genes; 464 human IDs	Vacuole (21, 3.3); lysosome (19, 3.5)	5.00
	Glycoprotein (141, 1.4); signal peptide (111, 1.4)	4.31
	Response to organic substance (38, 2.1); response to hormone stimulus (20, 2.2); response to insulin stimulus (10, 4.1)	3.03
	Defense response (26, 1.7); inflammatory response (19, 2.4)	2.86
	Cell junction (21, 1.6); adherens junction (11, 2.8); focal adhesion (9, 3.5)	2.41
	Cytoplasmic vesicle (29, 1.8); membrane-bounded vesicle (23. 1.6)	1.95
Cluster 8: Highest level at Day 16; 644 genes; 619 human IDs	Glycoprotein (216, 1.6); signal peptide (175, 1.7)	10.09
	Defense response (46, 2.2); inflammatory response (35, 3.1)	6.81
	Cell adhesion (47, 1.9); cell-cell adhesion (22, 2.3)	4.59
	Cell-cell adhesion (22, 2.3); calcium-dependent cell-cell adhesion (8, 10.1); synaptogenesis (6, 5.4)	3.39
	Carbohydrate binding (23, 1.9); pattern binding (15, 2.9)	2.98
	Ion homeostasis (30, 2.1); metal ion homeostasis (18, 2.5); calcium ion homeostasis (14, 2.2)	2.61
	Cell motion (34, 2.1); cell migration (19, 2.0)	2.54
	Regulation of growth (21, 1.8); regulation of cell growth (17, 2.5); negative regulation of cell growth (9, 2.8)	1.99

^a

Values within parentheses indicate the number of genes and fold-enrichment, respectively, of the functional term.

^b

Geometric mean of member's P-values of the corresponding annotation cluster (in −log10 scale).

Because to our knowledge the present study is the first systematic analysis of gene expression changes during the estrous cycle in equine endometrium, a comparison of the obtained expression profiles during the cycle to a similar study in bovine endometrium was performed. The data set of Mitko et al. [19] consisted of 268 known genes. Comparison was performed on the basis of assigned putative human orthologous genes. This revealed 133 genes found in both species as differentially expressed during the estrous cycle (Supplemental Table S6). Genes with high correlation and genes with inverse correlation of expression profiles during the cycle, respectively, are shown in Figure 5. Functional annotation clustering of genes with good correlation revealed overrepresentation of genes related to ECM, cell adhesion, and vasculature development. Genes with poor or even inverse correlation were predominantly assigned to functional categories such as endoplasmic reticulum, defense response, and extracellular region.

Fig. 5

Comparison of gene expression changes during the estrous cycle in equine and bovine endometrium. Gene expression changes during the estrous cycle in equine endometrium were compared to data obtained from bovine endometrium in a previous study [19]. Correlation of expression changes was calculated in Microsoft Excel. Relative gene expression levels (vsn-normalized value of one time point minus mean of vsn-normalized values of all five time points) are shown for the 30 genes with highest correlation (correlation > 0.83; left) and for the 30 genes with lowest correlation (correlation < −0.41; right). Genes were clustered according to their expression profiles (MeV software, Version 4.7.3). Color scale is from 4-fold lower (log2-fold = −2) than the mean (blue) to 4-fold higher (log2-fold = 2) than the mean (red). D, day of the estrous cycle.

Genes involved in PG metabolism and signaling and in P4 and estrogen signaling play an important regulatory role during the estrous cycle. Related genes were selected from the DEGs and used for a cluster analysis (Fig. 6). These genes showed a variety of distinct expression profiles during the cycle. Estrogen receptor 1 mRNA (ESR1) showed a profile very similar to that of G protein-coupled estrogen receptor 1 mRNA (GPER), with highest levels on Day 3 and lowest on Day 12, whereas P4 receptor mRNA (PGR) had highest expression on Day 0, decreasing to Day 12. In contrast, estrogen-related receptor alpha (ESRRA), estrogen-related receptor gamma (ESRRG), and P4 receptor membrane component 1 (PGRMC1) mRNAs showed completely different expression changes during the estrous cycle compared to the classical nuclear receptors. Messenger RNA levels of PGE synthases (PTGES, PTGES2, and PTGES3) were, on average, highest during the luteal phase. PGE receptor mRNAs (PTGER2, PTGER3, and PTGER4) showed completely different expression profiles. For most of the genes involved in arachidonic acid metabolism, highest mRNA levels were found during estrus. Two genes involved in PGD synthesis, PGD2 synthase 21kDa (PGDS) and hydroxyprostaglandin dehydrogenase 15-(NAD) (HPGD) were also regulated during the estrous cycle with highest levels on Day 3 and Day 0, respectively.

Fig. 6

Expression profiles of genes related to PG metabolism and signaling and to steroid hormone signaling. For PG metabolism and signaling, related genes were obtained from corresponding Gene Ontology categories (http://www.geneontology.org), a search of the Entrez Gene database (http://www.ncbi.nlm.nih.gov/gene), and the literature. For steroid hormone signaling, only selected genes were used (e.g., classical and nonclassical hormone receptors and related receptors and some genes encoding steroid hormone metabolic enzymes). Relative expression values of DEGs matching this candidate list were subjected to hierarchical clustering for gene expression during the estrous cycle in equine endometrium (MeV software, Version 4.7.3). Color scale is from 2-fold lower (log2-fold = −1) than the mean (blue) to 2-fold higher (log2-fold = 1) than the mean (red). D, day of the estrous cycle.

Finally, the expression of genes previously found as differentially expressed on Day 12 of pregnancy (compared to Day 12 of the estrous cycle) in equine endometrium [14] was analyzed during the estrous cycle. Supplemental Table S7 shows relative expression levels of these genes during the estrous cycle and corresponding q-values obtained from statistical analysis and expression changes for Day 12 of pregnancy versus Day 12 of the estrous cycle. Genes were sorted based on their expression profiles using SOTA analysis, resulting in six major expression profiles as shown in Supplemental Figure S3 (details can be found in Supplemental Table S7). Table 2 summarizes the number of genes assigned to the expression clusters of the SOTA analysis limited to the DEGs for Day 12 of pregnancy (Supplemental Fig. S3) and the corresponding expression clusters of the first SOTA analysis (Fig. 2 and Supplemental Table S2). Basically, almost all genes down-regulated on Day 12 of pregnancy showed their highest expression during the luteal phase (i.e., on Day 8 of the estrous cycle). In contrast, the majority of genes with higher expression in pregnant compared to Day 12 nonpregnant endometrium showed their highest expression on Day 0 or Day 16 of the estrous cycle (Table 2). A considerable number of genes with higher expression in pregnant endometrium showed no change during the cycle (Table 2, last column). Results of a DAVID analysis are shown in Table 3 (details can be found in Supplemental Table S8).

Table 2

Summary of the analysis of gene expression profiles during the estrous cycle for genes found as differentially expressed on Day 12 of pregnancy compared to Day 12 of the estrous cycle in equine endometrium.

		SOTA clustera
SOTA2 clusterb	Genesc	1 (Day 0)	2 (Days 0–3)	3 (Day 3)	4 (Days 3–8)	5 (Day 8)	6 (Days 8–12)	7 (Days 12–16)	8 (Day 16)	NSd
Genes with lower levels on Day 12 of pregnancy
1 (Day 0)	1	0	0	0	0	0	0	0	0	1
2 (Days 3–8)	23	0	0	0	5	16	0	0	0	2
3 (Days 8–12)	17	0	0	0	0	9	6	0	0	2
4 (Day 12)	1	0	0	0	0	0	0	0	0	1
5 (Days 12–16)	0	0	0	0	0	0	0	0	0	0
6 (Day 16)	0	0	0	0	0	0	0	0	0	0
Genes with higher levels on Day 12 of pregnancy
1 (Day 0)	134	46	8	2	0	0	0	0	34	44
2 (Days 3–8)	11	0	0	0	2	0	0	0	0	9
3 (Days 8–12)	13	0	0	0	0	1	3	0	0	9
4 (Day 12)	26	0	0	0	0	0	2	11	0	13
5 (Days 12–16)	32	0	0	0	0	0	0	25	0	7
6 (Day 16)	114	0	0	1	0	0	0	4	74	35

^a

Number of the expression cluster obtained from the analysis of all genes found as differentially expressed during the estrous cycle (see Fig. 2).

^b

Number of the expression cluster obtained from the analysis of genes found as differentially expressed on Day 12 of pregnancy (see Supplemental Fig. S3).

^c

Number of genes assigned to an expression cluster.

^d

NS, no significant change (q-value > 0.01).

Table 3

Results of a DAVID analysis of the genes differentially expressed on Day 12 of pregnancy compared to Day 12 of the estrous cycle in equine endometrium.

Group of genes	Functional term	Counta	P-valueb	FEc
Down-regulated on Day 12 of pregnancy, highest expression during luteal phase (SOTA2 clusters 2 and 3)	Transferase	7	0.033	2.8
	GO:0005739∼mitochondrion	6	0.065	2.6
	Transport	6	0.168	2.0
Up-regulated on Day 12 of pregnancy, highest expression during estrus (SOTA2 cluster 1)	GO:0005576∼extracellular region	34	0.000	2.9
	GO:0001568∼blood vessel development	10	0.000	7.9
	GO:0007155∼cell adhesion	14	0.000	3.9
	GO:0001558∼regulation of cell growth	6	0.003	6.0
	GO:0030246∼carbohydrate binding	7	0.008	3.9
Up-regulated on Day 12 of pregnancy, highest expression during late luteal phase (mainly genes from SOTA2 clusters 4 and 5)	Signal peptide	17	0.003	2.1
	GO:0008277∼regulation of G -protein coupled receptor protein signaling pathway	3	0.007	23.4
	GO:0005576∼extracellular region	12	0.010	2.2
	GO:0055085∼transmembrane transport	5	0.048	3.5
Up-regulated on Day 12 of pregnancy, highest expression at early estrus (mainly genes from SOTA2 cluster 6)	IPR006210:EGF-like	8	0.000	8.0
	GO:0005576∼extracellular region	26	0.000	2.2
	GO:0001568∼blood vessel development	8	0.000	6.0
	GO:0031410∼cytoplasmic vesicle	11	0.003	3.0
	GO:0007155∼cell adhesion	9	0.034	2.4
Up-regulated on Day 12 of pregnancy, no significant change during the estrous cycle (fold-change from highest to lowest expression < 2, q-value > 0.01)	GO:0005576∼extracellular region	19	0.001	2.2
	IPR013783:Immunoglobulin-like fold	8	0.009	3.3
	IPR000742:EGF-like, type 3	5	0.009	6.0
	GO:0007155∼cell adhesion	8	0.025	2.7
	GO:0048545∼response to steroid hormone stimulus	4	0.045	4.9
	GO:0060429∼epithelium development	4	0.068	4.2
	GO:0001568∼blood vessel development	4	0.081	3.9

^a

Number of genes assigned to a term.

^b

Significance of overrepresentation (enrichment) of a functional term.

^c

FE, fold-enrichment (overrepresentation) of a functional term.

Discussion

To our knowledge, the present study provides the first insights regarding global dynamic changes in mRNA expression profiles during the estrous cycle in equine endometrium. To analyze these changes, endometrial biopsy samples were collected from five mares at five different time points of the estrous cycle. Two time points were chosen for the estrous phase (low P4 concentrations): Day 16 for beginning of estrus and Day 0 as the day of ovulation. For the luteal phase, three time points were selected for sample collection: Day 3 for the early luteal phase, with P4 concentrations on the rise; Day 8 for the midluteal phase; and Day 12 for the late luteal phase, with P4 concentrations falling. Overall, the analysis of sample correlation (pairwise correlation of microarray data sets) showed grouping of samples according to the different phases of the cycle, with the greatest difference between luteal and ovulatory phase but also between the time points of these two stages. This result is similar to previous findings in bovine endometrium during the estrous cycle [19].

Statistical analysis of microarray data revealed almost 5000 differentially expressed genes/transcripts. A classification of these genes according to their gene expression during the estrous cycle revealed eight clusters of similar expression profiles, which contained relatively equal numbers of genes (501–721 genes). This result is different to the findings of the bovine study [19], in which most of the genes had highest expression levels either at estrus (Day 0) or diestrus (Day 12). This difference can be explained by the limited number of genes that were analyzed in the bovine study using a specialized cDNA array [20], which mainly contained genes found as differentially expressed between estrus and diestrus in bovine endometrium and oviduct and on Day 18 of pregnancy in bovine endometrium [21–25]. The analysis of all bovine genes during the estrous cycle would probably give results similar to those obtained for equine endometrium. However, the results of the present study agree well with cyclic changes of ovarian steroid hormones and findings from histological studies. The bioinformatics analysis of the clusters of genes with similar expression profiles during the estrous cycle revealed specific gene functions and biological processes for different phases of the cycle.

Genes with Highest mRNA Concentrations During Estrus

Similar to cattle [23], many of the genes with highest mRNA concentrations on Day 0 were related to the ECM (e.g., 18 genes coding for collagens) and the process of focal adhesion (i.e., adhesion of cells to the ECM). The process of protein translation was also highly overrepresented during late estrus (after ovulation) and the early luteal phase. The changes during the ovulatory phase are mainly controlled by E2. The endometrial mRNA concentration for ESR1 also showed significant differences during the estrous cycle, started to increase on Day 16, and was highest on Day 3. This agrees with the results of Hartt et al. [15], who found highest mRNA levels on Days 0, 17, and 20 of the estrous cycle but did not investigate Day 3. Fibroblast growth factor (FGF) 9 has been described as one of the mediators of estrogen signaling in the endometrium [26]. FGF9 mRNA concentration increased on Day 16 and was highest on Day 0 following serum E2 levels. The mRNAs of two estrogen-related receptors, ESRRA and ESRRG, and of GPER were also found as differentially expressed during the estrous cycle. GPER showed an expression profile similar to that of ESR1, with highest levels on Day 3. In human endometrium, GPER expression is also regulated during the cycle, with highest mRNA levels during the proliferative phase [27]. Based on the present finding that a cluster of genes showing a peak of their expression levels on Day 3 was highly enriched for functions related to cell proliferation, ribosome biogenesis, and mitochondrion, GPER expression seems to be associated with increased cell proliferation in mare endometrium as well. Furthermore, increased cell proliferation has been found in a previous study of equine endometrium during early diestrus [28]. During this time, E2 levels decrease, and P4 levels are rising. In contrast to ESR1 and GPER, ESRRG showed lowest expression levels on Day 3 and highest levels on Day 0, and ESRRA had lowest levels on Day 0 and highest on Day 8, with higher levels on Day 16. Estrogen-related receptors are a subfamily of orphan nuclear receptors closely related to the classical estrogen receptors and have been shown to share target genes, coregulatory proteins, and so on with classical ESRs and can influence the estrogenic response [29]. Although ESRRs lack identified natural ligands (orphan nuclear receptor), binding and activation by phytoestrogens have been shown [30]. Altogether, the expression patterns of these selected genes involved in estrogen signaling indicate a complex regulation of E2 effects on the endometrium, and further detailed studies of dynamic protein expression and localization of expression in different endometrial compartments need to be performed to understand the interaction of these genes and their products.

In addition to genes related to ECM constituents, several functional categories for immune system functions were found to be overrepresented for the genes with highest mRNA concentrations during the estrous phase. These categories comprised more than 100 genes, including 25 genes coding for cytokines, 16 genes involved in innate immune response, 9 genes involved in adaptive immune response, 8 genes involved in complement activation, and 8 cytokine receptor genes. The cytokine genes represented the chemokine (C-C motif) ligand family (CCL2, CCL3L3, CCL5, CCL8, CCL13, CCL15, and CCL21), chemokine (C-X-C motif) ligand family (CXCL1, CXCL5, CXCL6, CXCL9, CXCL10, CXCL16, and CXCL17), and the interleukin family (IL8, IL15, IL32, and IL34), which are attractors for almost all types of immune cells, such as monocytes, macrophages, dendritic cells, neutrophils, B cells, T cells, and natural killer cells. This finding and the increased mRNA concentrations of genes coding for specific surface antigens indicate infiltration with immune cells and activation of the endometrium during estrus. So far, this process has been investigated only in the context of mating or insemination [31]. The present results clearly show that the endometrium prepares for the uterine clearance needed during and after mating in the estrous phase, which is probably further triggered by spermatozoa and components of the seminal plasma [31], leading, for example, to influx of polymorphonuclear neutrophils into the uterine lumen [32].

Genes with Highest mRNA Concentrations During the Luteal Phase

Similar to cattle and swine, expression of PGR in equine endometrium decreases during the luteal phase, particularly in luminal and glandular epithelium [15]. Effects of P4 during the time of highest P4 levels are mainly mediated via stromal cells that still express PGR [33] and via other types of nonclassical P4 receptors. In the present study, expression of PGRMC1 mRNA was highest at Day 8, in contrast to PGR expression, which was highest at Day 0. The associated “late P4 response” genes (highest levels on Days 8 and 12) were predominantly related to protein secretion, as indicated by functional categories, such as protein transport, endoplasmic reticulum, and Golgi apparatus, to signaling processes and to metabolic processes, such as monosaccharide metabolic process and lipid metabolism. Several overrepresented functional categories related to lipid and fatty acid metabolism were obtained. A more detailed analysis for the categories lipid metabolism, lipid biosynthetic process, and fatty acid metabolic process revealed a complex regulation of genes belonging to different related pathways. Specifically, genes involved in fatty acid metabolism, adipocytokine signaling pathway, peroxisome proliferator-activated receptor signaling pathway, and glycerophospholipid metabolism were found in the expression clusters with highest mRNA levels on Days 8 and 12 of the cycle. For the adipocytokine signaling pathway, both adiponectin receptors 1 and 2 (ADIPOR1 and ADIPOR2) were identified. Expression of these receptors also has been described in human and murine endometrium during the window of implantation time and in preimplantation endometrium, respectively [34, 35]. Highest mRNA concentrations for genes coding for enzymes involved in fatty acid metabolism, such as acyl-coenzyme A (CoA) synthetases; acyl-CoA dehydrogenase (ACADL); carnitine palmitoyltransferases, acetyl-CoA acyltransferase 2 (ACAA2); acetyl-CoA carboxylase beta (ACACB); enoyl-CoA delta isomerase 2 (ECI2); hydroxyacyl-CoA dehydrogenase/3-ketoacyl-CoA thiolase/enoyl-CoA hydratase (trifunctional protein), beta subunit (HADHB); and ELOVL fatty acid elongase 5 (ELOVL5), were observed on Days 8 and 12 of the estrous cycle, indicating increased fatty acid synthesis, which is needed for provisioning of the preattachment equine conceptus (in case of pregnancy) with essential fatty acids and retinoids. Suire et al. [36] suggested that equine uterocalin (official name and symbol: P19 lipocalin, P19) binds fatty acids and lipids and transports them to the conceptus. Messenger RNA expression profile of P19, coding for one of the major P4-dependent proteins secreted into the uterine lumen of the mare [37], was similar to those of genes involved in fatty acid metabolism, with highest expression on Days 8 and 12 and lowest on Day 0 of the estrous cycle.

The functions of the genes with highest levels on Days 8 and 12 of the estrous cycle (protein secretion, signaling, and metabolic processes) indicate that the equine endometrium prepares for support of embryo growth and development. Studies of asynchronous embryo transfer showed a relatively wide window for establishment of pregnancy in the mare [38]. Also, conceptus growth and development was retarded when Day 10 embryos were placed in negatively asynchronous uteri (up to −7 days) [38]. This agrees well with the highly significant overrepresentation of molecular functions and biological processes related to cell proliferation for genes with highest mRNA levels on Day 3 and highest expression levels of “supportive” genes with functions related to protein secretion and nutrient supply on Days 8 and 12 of the estrous cycle. In this regard, the situation in equids seems to be different compared to that in ruminants, where an adequate timing of P4-regulated changes in the endometrium is the prerequisite for normal embryo development and conceptus elongation [39, 40].

Comparison to Dynamic Gene Expression Changes During the Estrous Cycle in Bovine Endometrium

The comparison to published data from bovine endometrium when comparable phases of the estrous cycle were analyzed [19] revealed a number of genes with very similar mRNA expression profiles. Some of those genes, such as collagen genes, tenascin C (TNC), and claudin 4 (CLDN4), showed similar expression changes in murine and human endometrium as well [41]. The genes with low or even inverse correlation of expression profiles during the estrous cycle were predominantly related to immune system functions, protein secretion, and proteins secreted to the extracellular region. Whereas highest mRNA levels for immune-related genes have been found during the mid and late luteal phase in bovine endometrium [19], immune-related genes were mainly up-regulated during estrus in equine endometrium. This agrees well with the reported role of the immune system in uterine clearance after mating [32]. Although the mares used in the present study were not inseminated, the corresponding immune-related genes (inflammatory response, cytokine-cytokine receptor interaction, complement and coagulation cascades, and innate immune response) are up-regulated during estrus probably by effects of estrogen in preparation for the shock of intrauterine insemination in the mare. The different expression profiles of genes related to protein secretion—higher levels during diestrus in the mare and during estrus in cattle [19]—could be due to specific requirements of the developing conceptus in cattle and horses. The major protein secreted to the uterine lumen in the mare is encoded by P19 [42], which had highest mRNA concentrations on Days 8 and 16 of the estrous cycle (18-fold higher on Day 8 vs. Day 0). P19 is thought to function as a transporter to provide the equine conceptus with lipids essential for growth, morphogenesis, and pattern formation [36]. The genes coding for extracellular proteins also contained, in addition to immune-related genes, some genes related to ECM remodeling (transglutaminase 2, TGM2 [43]) and to epithelial barrier formation (dystonin, DST [44]; cystatin E/M, CST6 [45]).

Overall, the identification of genes that show completely different expression profiles in different species, such as cattle and horses, indicates an adaptation of endometrial gene expression to specific differences in embryonic development and differences in the mode of maternal recognition of pregnancy (MRP). In the horse, a spherical conceptus is enveloped by a tough and elastic glycoprotein capsule [46]; in cattle, the conceptus becomes filamentous and secretes a type I interferon to achieve MRP [47]. Although, the equine conceptus secretes a number of different potential factors for MRP, the nature of MRP in the mare remains unclear [48].

Genes Involved in PG Metabolism and Signaling

Most of the genes involved in arachidonic acid metabolism (i.e., PG metabolism) and PG signaling showed highest expression levels during the late luteal phase and the ovulatory phase (Days 12, 16, 0, and 3) of the estrous cycle. The expression profiles identified in the present study show a higher time resolution compared to data from previous publications, because fewer or less-well-defined time points have been analyzed in other studies [12, 13]. Overall, the present results were in accordance with those of other studies, except for the PGF synthase gene (PGFS), which showed expression changes during the cycle that did not reach statistical significance (q-value = 0.03; threshold q-value = 0.01). In our recent microarray study of equine endometrium [14], we found increased expression of PTGER3 and PTGER4; phospholipase A2, group IVA (PLA2G4A); and solute carrier organic anion transporter family, member 2A1 (SLCO2A1; PG transporter) on Day 12 of pregnancy compared to nonpregnant controls. For all four genes, highest mRNA concentrations were found on Days 16 and 0 during the estrous cycle, indicating an earlier rise of mRNA levels during pregnancy compared to the cycle. Altogether, the complex regulation of genes involved in lipid, fatty acid, and PG metabolism and signaling suggests an important role of the corresponding genes in the regulation of the estrous cycle. This is further supported by the results of a recent study in which the intrauterine administration of plant oils was shown to postpone luteolysis [49]. The authors of that study concluded that fatty acids contained in plant oils could modulate the synthesis or the release of PGs from the endometrium [49].

Genes Related to Angiogenesis, Vasculature Development, and Blood Circulation

The expression levels of genes assigned to the categories angiogenesis, vasculature development, and blood circulation was lowest during the early luteal phase. Most of the genes assigned to the Gene Ontology category angiogenesis showed their highest mRNA concentrations either at estrus or on Day 8. The genes assigned to the category vasculature development with highest levels during estrus were almost identical to the genes assigned to angiogenesis. During estrus, mRNA levels were highest for members of the angiopoietin system (angiopoietin 2, ANGPT2; angiopoietin-like 1, 2, and 4, ANGPTL1, ANGPTL2, and ANGPTL4, respectively; TEK tyrosine kinase, TEK), the VEGF system (a coreceptor for VEGF, neuropilin 2, NRP2; vascular endothelial growth factor A, VEGFA; kinase insert domain receptor, KDR, also known as fetal liver kinase-1 [FLK1] or vascular endothelial growth factor receptor 2 [VEGFR2], genes related to endothelial cell regulation (endothelins 1 and 2, EDN1 and EDN2; the receptor for EDN1, endothelin receptor type A, EDNRA; endothelial cell surface expressed chemotaxis and apoptosis regulator, ECSCR), angiotensin I-converting enzyme (peptidyl-dipeptidase A) 2 (ACE2), renin (REN), and FGF2. Although a study of vascular density in the equine endometrium did not reveal significant differences between the follicular and the luteal phase [50], the up-regulation of these genes are signs of vascular changes, which could be related to the observation that “inflammatory response” genes are overrepresented during estrus, because inflammatory responses and vascular responses are tightly connected [51]. Furthermore, up-regulation of a number of angiogenic factors by E2 has been shown in sheep [52]. During diestrus, the genes for nitric oxide synthase 3 (NOS3), angiotensin I-converting enzyme (peptidyl-dipeptidase A) 1 (ACE), angiotensinogen (AGT), angiomotin (AMOT), angiopoietin-like 6 (ANGPTL6), and endothelial PAS domain protein 1 (EPAS1) showed their highest mRNA levels. For NOS3 protein, differential expression has been described in human endometrial vessels, with highest expression in the early secretory phase [53]. Angiotensinogen is cleaved by renin to angiotensin and subsequently to peptides with either vasodilator or vasoconstrictor activities by angiotensin-converting enzymes (see KEGG pathway ko04614 Renin-angiotensin system) [54]. Interestingly, mRNA coding for ACE2 that produces vasodilators from angiotensin I and II was highest at estrus, whereas mRNA coding for ACE that produces the potent vasopressor angiotensin II was highest on Day 8 of the estrous cycle. In correlation with the observation that the genes assigned to the functional categories blood circulation and angiogenesis were predominantly found in the expression clusters with highest mRNA concentrations during estrus and in expression clusters with highest mRNA concentrations on Days 8 or 12, studies of uterine blood flow in cyclic mares found highest blood flow on Days 0 and 10/11 [55, 56]. The comparison to studies of gene expression during the time around MRP [14, 16] shows that expression of similar genes is induced during estrus and on Days 12 and 14 of pregnancy, indicating a specific induction of those genes by the presence of the conceptus, probably by conceptus estrogens.

Genes Found as Differentially Expressed on Day 12 of Pregnancy Compared to Day 12 of the Estrous Cycle

To draw additional conclusions for the genes previously identified as differentially expressed on Day 12 of pregnancy [14], expression data for these genes during the estrous cycle were analyzed separately. The main findings—that genes down-regulated on Day 12 of pregnancy show their highest expression during the luteal phase and that the majority of genes up-regulated on Day 12 of pregnancy show their highest expression on Day 16 or Day 0 of the cycle—suggest regulation of these genes by conceptus estrogens. For example, one of the genes with highest up-regulation in Day 12 and Day 13.5 pregnant endometrium, solute carrier family 36 (proton/amino acid symporter), member 2 (SLC36A2), showed highest expression on Day 16 of the cycle. The observation that the probes covering the 3′-untranslated region of the SLC36A2 transcript show this strong regulation but not the single probe located in the open reading frame revealed only 2-fold higher expression in pregnant endometrium, and no change during the cycle could indicate transcript isoforms for this gene. Vice versa, the genes with strongest down-regulation in Day 12 pregnant endometrium, FXYD domain containing ion transport regulator 4 (FXYD4) and keratin 4 (KRT4), showed 50-fold higher levels on Day 8 compared to Day 0 of the estrous cycle. KRT4, a marker of mucosal epithelia [57], has been shown to be down-regulated by E2 [58]. Furthermore, a number of the genes up-regulated in Day 12 pregnant endometrium showed their highest expression during the luteal phase. The observation that almost all of these genes had their highest expression levels on Day 12 of the estrous cycle and only a few of them during the early (Day 3) and mid (Day 8) luteal phase indicates that expression of these genes correlates with down-regulation of PGR in luminal and glandular epithelium [15]. Further up-regulation of these genes in pregnant endometrium could also be due to conceptus estrogen. For example, a similar expression pattern has been found for the gene coding for uterine serpins (SERPINA14) in bovine endometrium, with lowest levels on Day 3.5 of the estrous cycle and highest at early estrus [59]. In vitro, a further increase of SERPINA14 mRNA by treatment with E2 was shown [59]. In addition, up-regulation by E2 has been shown for the genes transient receptor potential cation channel, subfamily V, members 5 and 6 (TRPV5 and TRPV6) [60, 61]. However, because E2 has been shown not to be the equine pregnancy recognition signal [49], all these potential E2-regulated genes may only be of secondary interest. No genes were observed that show a clear increase of mRNA expression from Day 8 to Day 12 of the estrous cycle and decreased expression in Day 12 pregnant endometrium, which could be expected from genes involved in the initiation of luteolysis. The last group of genes (88 genes) with increased expression in samples from Day 12 pregnant endometrium showed no significant change during the estrous cycle (q-value > 0.01 and fold-change from lowest to highest signal < 2). No change of expression during the analyzed days of the cycle indicates that these genes are not regulated by E2 and/or P4. Two of these genes, nuclear receptor subfamily 2, group F, member 2 (NR2F2) and tumor necrosis factor (ligand) superfamily, member 10 (TNFSF10), were previously suggested to have a conserved role in different mammalian species for establishment and maintenance of pregnancy [14, 62]. DAVID functional annotation revealed a variety of different functional terms related to immune functions; development of hematopoietic cells, endothelial cells, and epithelial cells; and cell adhesion processes. In particular, the genes related to developmental processes could have a role in endometrial differentiation in the context of establishment of pregnancy.

In conclusion, the present study revealed new insights regarding the complex dynamic gene expression changes during the estrous cycle in equine endometrium. The comparison to data from bovine endometrium revealed similarities but also specific differences, pointing to species-specific regulations. Data analysis with a focus on specific functional categories and biological processes in relation to dynamic gene expression profiles provided the molecular reflections of important physiological changes in the cyclic equine endometrium with regard to the crucial role of this tissue for successful reproduction.

Acknowledgment

We thank Dr. E. Senckenberg from the Bavarian Principal and State Stud of Schwaiganger, Germany, for providing mares for the present study. Furthermore, we thank Karin Groβ and Andrea Klanner for excellent technical assistance.

Author notes

¹

Supported by the German Ministry for Education and Research (BMBF, FUGATO-plus, COMPENDIUM) and the German Research Foundation (DFG, Research Training Unit 1029 Functional Genome Research in Veterinary Medicine). The data discussed in this publication have been deposited in the NCBI Gene Expression Omnibus (GEO, ) and are accessible through GEO Series accession number GSE39043.