Characterization of the endometrial transcriptome in early diestrus influencing pregnancy status in dairy cattle after transfer of in vitro-produced embryos

Characterization of the endometrial transcriptome in early diestrus inﬂuencing pregnancy status in dairy cattle after transfer of in vitro-produced embryos. Physiol Genomics of the endometrial transcriptome at day 7 of the estrus cycle are crucial to maintain gestation after transfer of in vitro-produced (IVP) embryos, although these changes are still largely unknown. The aim of this study was to identify genes, and their related biological mechanisms, important for pregnancy establishment based on the endometrial transcriptome of recipient lactating dairy cows that become pregnant in the subsequent estrus cycle, upon transfer of IVP embryos. Endometrial biopsies were taken from Holstein Friesian cows on day 6–8 of the estrus cycle followed by embryo transfer in the following cycle. Animals were classiﬁed retrospectively as pregnant (PR, n (cid:2) 8) or nonpregnant (non-PR, n (cid:2) 11) cows, according to pregnancy status at 26–47 days. Extracted mRNAs from endometrial samples were sequenced with an Illumina platform to determine differentially expressed genes (DEG) between the endometrial transcriptome from PR and non-PR cows. There were 111 DEG (false discovery rate (cid:3) 0.05), which were mainly related to extracellular matrix interaction, histotroph metabolic composition, prostaglandin synthesis, transforming growth factor- (cid:4) signaling as well as inﬂammation and leukocyte activation. Comparison of these DEG with DEG identiﬁed in two public external data sets conﬁrmed the more fertile endometrial molecular proﬁle of PR cows. In conclusion, this study provides insights into the key early endometrial mechanisms for pregnancy establishment, after IVP embryo transfer in dairy cows. found at the top or bottom of the list are those that contribute the most to the enrichment results; these genes deﬁne the core enrichment. The enrichment score is computed using a weighted Kolmogorov–Smirnov-like statistic. The normalized enrichment score is computed to account for sizes of the gene sets. We used the gene


INTRODUCTION
The fertilization rate after artificial insemination (AI) of dairy cattle is estimated to be 80 -90%, whereas the average calving rate is 55% (14,62). This high pregnancy loss occurs during the first part of the embryonic period (78), mainly between days 7 and 16 of gestation (5), due mostly to the inability of the uterus to support embryo implantation and growth (67).
Successful establishment and maintenance of a pregnancy are dependent on both embryonic viability and the uterine environment sustaining embryonic development. With respect to embryonic viability, transfer of in vitro-produced (IVP) embryos generally results in lower pregnancy rates (30 -40%) as compared with embryos produced in vivo after superovulation and AI of the donor cow (55-60%) (53,55). IVP embryo quality is dependent on oocyte competence, sperm quality, and in vitro handling, and despite many efforts to improve media and methodologies, IVP and in vivo-produced embryos still deviate in several ways (10,18,35), due to which optimization of the IVP process is subject to intense investigation yet [reviewed by Lonergan (39)]. With respect to the uterine environment, estradiol produced by the follicles around estrus primes the endometrium (41), and progesterone produced by the corpus luteum (CL) in diestrus induces endometrial receptivity (66), which is known as the period when the uterine environment is physiologically ideal for establishment of pregnancy. During diestrus, the endometrial stroma is thickened and uterine glands are enlarged, resulting in an increased production of histotroph (49,64) that is absorbed by the trophoblast of the potential embryo (16). The thickness of the endometrium depends on the concentration of circulating estradiol and progesterone, and it has been shown to vary between natural or induced estrus (72).
Endometrial changes occur in every reproductive cycle, independently of the presence of a conceptus (21). Interestingly, the endometrial transcriptome underlying these changes may convey higher endometrial receptivity in some animals than in others. Hence, the endometrial transcriptome at day 7 differs between heifers that maintain the pregnancy and those that do not (54). On top of these cyclic changes in the endometrial transcriptome, signals from an embryo further modify endometrial gene expression as early as day 7 in the estrus cycle (68). Together, the endometrial transcriptome varies over the estrous cycle, between receptive and less receptive individuals and with respect to whether an embryo is present or not.
Regarding endometrial receptivity, experiments performed in Nellore cows and Simmental heifers (6,54,60) after estrus synchronization have unveiled important differences in the endometrial transcriptome at day 6 or 7 of animals that become pregnant after AI or embryo transfer (ET), respectively, compared with the nonpregnant ones. However, disparities in embryo survival between heifers and lactating cows (62), the different effects of progesterone on the fertility of dairy and beef cows (69), and the variations in hormonal concentration and uterine thickness between induced or natural estrus (72) indicate the need of further studies for a better comprehension of the mechanisms controlling endometrial receptivity.
The aim of this study was to determine biological processes and pathways that convey receptivity in the endometrium for establishment of pregnancy based on the identification of differentially expressed genes (DEG) in the endometrium on day 6 -8 in diestrus of recipient lactating dairy cows that become pregnant in the subsequent estrus cycle after transfer of IVP embryos compared with those that do not establish pregnancy. Thus, cows were classified retrospectively as pregnant (PR) or nonpregnant (non-PR) according to their status.
Furthermore, results obtained from the present study were compared with the endometrial transcriptome differentially expressed at day 7 between cows with high or low merit for fertility (47) or cows with high or lower levels of progesterone (19). The goal was to employ a bioinformatic approach to test the hypothesis that the observed differences in the transcriptome are attributed more to the inherent endometrial molecular signature of the animal than to hormonal -and other-variations between cycles.

Ethical Statement
This study was carried out in strict accordance with the Danish law on animal experimental work (law no. 474 of 15 May 2014). Accord-ing to this law, an application to perform the live animal work in this study was submitted and approved by the ethical committee "Danish Animal Experiments Inspectorate" before starting the experiment (license no. 2015-15-0201-00636).

Animals
The animals were Holstein Friesian cows in the middle of their first, second, or third lactation, with no clinical disease, in good body condition, and with no history of reproductive problems. All estruses were spontaneous during the experimental period, so no hormonal synchronizations were used. Estrus observation was performed visually twice daily and by automatic monitoring of cow activity in 6 h intervals. Before being included in this experiment, the cows had at least three regular estrus cycles and passed a thorough gynecological examination. Endometrial biopsies were taken from cows on day 6 -8 in diestrus (day 0 ϭ last day of standing estrus) followed by ET on day 6 -8 in the following cycle. The pregnancy status was determined when all cows were slaughtered on day 26 -47. Embryos were produced in two or three replicates per week to ensure availability of fresh embryos for transfer, according to the recipient's natural cycle, after a strict selection of the best-quality embryos for transfer as well as attempting synchrony between embryo developmental stage and recipient (3,500 oocytes fertilized with semen from one bull in 35 IVP rounds). For transfer, 28 blastocysts were selected, derived from 16 IVP replicates (1,600 oocytes fertilized) with an average blastocyst rate of 43% (min-max 26 -56%).
The experiment was continued until 12 pregnancies were achieved to obtain enough replicates to get reliable results with DESeq2 (40). All recipients received exclusively a single embryo, and none were, hence, reused for later embryo transfer. The overall pregnancy rate obtained was 43% from a total of 28 cows for ET. From the 16 non-PR cows sampled, 12 were selected for further analysis. Table 1 summarizes the main descriptive characteristics of the cows used in this study. The whole experiment took around 6 mo.

Endometrial Biopsies
On day 6 -8 of the cycle preceding embryo transfer, the cows had an epidural analgesia (5 ml Procamidor, Salfarm Danmark A/S, Cows were classified as pregnant (PR) or nonpregnant (non-PR) after transfer of in vitro-produced embryos. *Values are shown as means Ϯ standard deviation in parentheses. There were no statistical differences for these values between PR and non-PR cows (P Ͼ 0.05 after t test and Fisher exact test for continuous and categorical variables, respectively). ET, embryo transfer. Denmark), and an endometrial biopsy of 50 -60 mg was taken from the middle to cranial part of the uterine horn ipsilateral to the CL using a biopsy instrument (141700, Kruuse A/S, Denmark). The biopsies were transferred into Eppendorf Safe-Lock Microtube 2.0 ml with 1.5 ml RNAlater (76154, Qiagen), and kept at room temperature for a minimum of 2 h. Each biopsy was divided into two replicates (each of 25-30 mg tissue) and transferred into 0.5 ml RNAse-free Microfuge tubes (Ambion) and stored at Ϫ80°C.
The biopsies were taken from live cows. Therefore, they could not be taken precisely in relation to the caruncle sites. The biopsies were deep enough to include the full endometrium (luminal epithelium, glandular epithelium, stroma, blood vessels, nerves, and invading leukocytes). Exclusively endometrial biopsies of the expected weight, i.e., (50 mg) were included in the full bioinformatic analyses to achieve a comparable uniform material. Thus, two samples were excluded from the RNA-Seq analysis. These samples belonged to cows that were classified later as PR.  (29). On day 6 -8, each recipient cow had one fresh embryo transferred. On the transfer day (day 7 of embryo development, IVF ϭ day 0), the embryos were evaluated and selected based on developmental stage and morphology. Morphological good and excellent blastocysts/expanded blastocysts were used for transfer. For synchronization between embryo development stage and recipient, day 6 recipients had a blastocyst transferred while day 7 and 8 recipients had a blastocyst or an expanded blastocyst transferred. Before transfer, the embryos were washed in holding medium (Syngro Holding, Bioniche), loaded into 0.25 ml straws, and transported to the stable at 38°C. The cows had an epidural analgesia (5 ml Procamidor) before the one embryo was transferred to the middle to cranial part of the uterine horn ipsilateral to the CL.

Samples at Slaughter
From day of transfer until day 26 -47, the recipient cows had not shown any signs of estrus and were then slaughtered to collect their reproductive organs. The cows were slaughtered across different days to describe gradual changes in pregnancies based on IVP embryos, which were destined for another experiment. After slaughter, the uterine horns were opened, and based on the findings of embryos (26 -42 days of age), fetuses (above 42 days of age), and fetal membranes, the cows were registered as PR or non-PR. By this approach, we expected to find signs of an eventual early loss of an early pregnancy, because of the size of an early embryo with membranes, and because the uterus had not been cleaned by an estrus period.

RNA Extraction and Sequencing
Endometrial biopsies were homogenized with Tissue Lyser LT. Total RNA was isolated and purified with RNAeasy mini Kit. Quantity and quality of the RNA samples were assessed by Agilent 2100 Bioanalyzer. One sample was not sequenced because its RNA integrity number (RIN), a measure of the RNA quality in the range from 1 to 10, was 1.5, resulting in 11 samples from non-PR cows. The average RIN for the remaining samples was 7.6 (minimum ϭ 4.2, maximum ϭ 10).
The libraries for the sequencing were generated with an Illumina TruSeq-stranded total RNA protocol.

Availability of Data
All RNA-Seq data have been deposited in the National Center for Biotechnology Information's (NCBI) Gene Expression Omnibus (GEO) and are accessible without restriction through GEO accession number GSE115756 (https://www.ncbi.nlm.nih.gov/geo/query/acc. cgi?accϭGSE115756).
RNA samples were paired-end sequenced with the Illumina HiSeq 2500 platform with a read length of 100 nt, in two lanes pooling all the samples together.
Bioinformatic Analysis RNA processing. The quality of the reads was checked with FastQC (v. 0.11.2) (2). Adapters and quality trimmer and filtering were performed using Trimmomatic (v. 0.36) (7) to compute (1) trimming of the 3=-end and 5=-end with a quality threshold of 20 (2), window-based filtering using a window of four bases and an average quality threshold value of 15 (3), filtering reads with a final length of Ͻ25 nucleotides. The remaining 41,799,192 read pairs were mapped to the bovine reference genome (Bos taurus UMD3.1) (80) with STAR aligner (v. 2.5.2) (15), including the gene Bos taurus release 87 annotation in the genome index. We allowed for a maximum of five mismatches, while setting the other parameters to STAR default values. Read counts were estimated at gene level using HTSeq-count (v. 0.6.0) (1), setting the model of intersection as "intersectionnonempty" using the same annotation file.
RNA-Seq preprocessing. The samples 20 and 21, from PR cows, showed different profiles during library fragmentation and were excluded from the rest of the analysis (Supplemental Fig. S2). Therefore, the analyzed sample size for each group was PR ϭ 8 and non-PR ϭ 11. Further quality control on count data was performed with Noiseq (v.2.14.0) (73) and with the R package cqn (27). The quality control did not reveal any GC or length bias between samples (Supplemental Fig. S3).
Differential gene expression analysis. We filtered out genes with fewer than one count per million in more than eight samples (smaller class) (56), retaining 14,231 genes. The gene expression counts were normalized by library size with DESeq2 methods.
The differential gene expression analysis was performed by fitting a logistic regression model to the gene counts (modeled by a negative binomial distribution) in DESeq2 as: where yi is the gene normalized gene expression count for gene i, ␤0 is a fixed intercept term; parity i(z,ƒZ-1) is a set of (Z Ϫ 1) dummy variables representing the parity for the i th animal; and calving, ET representing the days between the last calving and the embryo transfer was fitted as covariate. 〉parity,z and ␤c_ET are the corresponding solutions. The differential expression analysis was based on the embryo transfer outcome that was fitted as a dummy variable (0 ϭ non-PR or 1 ϭ PR). We applied a Wald test statistic (40) to test for model significance. P values were adjusted with Benjamini-Hochberg method. DEG were called at false discovery rate (FDR) Ͻ0.05.
Functional enrichment analysis. The Gene Set Enrichment Analysis (GSEA) was performed using the GseaPreranked tool from GSEA v2.2.2 (Broad Institute) (71). This tool tests a functionally related set of genes that are defined a priori. The DEG overlapping with the gene sets that are found at the top or bottom of the list are those that contribute the most to the enrichment results; these genes define the core enrichment. The enrichment score is computed using a weighted Kolmogorov-Smirnov-like statistic. The normalized enrichment score is computed to account for sizes of the gene sets. We used the gene sets from Molecular Signatures Database (MSigDB Collections) v5.1. GSEA was able to map 12,391 genes based on the gene symbol.
The functional analysis with Ingenuity Pathway Analysis (IPA, QIAGEN) (https://www.qiagen.com/ingenuity) was performed using all the DEG, employing the Entrez ID as identifier. IPA mapped successfully 99 out of the 111 DEG. IPA defines upstream regulators as any type of molecules that could affect the expression of a specific set of genes. The upstream regulators, which are significantly overrepresented for the DEG list, are selected. If the directionality of the gene expression (upregulation, downregulation) of the DEGs is consistent with the a priori information of that upstream regulator, additional information about its predicted status (activated or inhibited) is provided.
The proportions of the DEG were not correlated to the length of the genes (Supplemental Fig. S4). Thus, the GO term enrichment analysis was performed with NETwork-based human gene enrichment (NET-GE) (8), using the symbol name for Biological Process (BP), Cellular Component (CC), and Molecular Function (MF) GO terms. NET-GE is network-based method for gene enrichment analysis. For each functional gene set NET-GE relies on an interactome to build a module of functionally related genes. The modules are then tested via Fisher's exact test.
Results of the functional analysis were considered statistically significant with an FDR Ͻ 5%.

Comparison with Other Studies
The experimental model used in this study is similar to the one used by Salilew-Wondim et al. (60) and Ponsuksili et al. (54) in the sense that endometrial tissue was collected in the cycle previous to the cycle where the embryo was transferred, to avoid interfering with establishment of pregnancy. A flaw in this model is that changes in estradiol and progesterone concentrations between cycles could be modifying the endometrial transcriptome, and thus the observed differences would be attributed to variations in the level of steroid hormones rather than to the inherent endometrial molecular signature of the animal. Unfortunately, hormonal concentrations in both cycles were not measured in the present experiment for reasons of logistics.
To rule out this concept, two publicly available gene expression data sets compared with our results. These data sets were downloaded from the NCBI GEO repository, with the accession numbers GSE52438 (47) and GSE33030 (19). In the first study (GSE52438), endometrial biopsy samples from nonpregnant lactating dairy cows on day 7 of the estrus cycle with similar genetic merit for milk production traits but with very good genetic merit for fertility [high fertility (HF), n ϭ 7] and cows with very poor genetic merit for fertility [low fertility (LF), n ϭ 7] were analyzed by RNA-Seq technology. The matrix of gene counts was downloaded and reanalyzed with the DESeq2 package for R (40). Genes with low expression counts (Ͻ1 CPM) in seven or more samples were filtered out before normalization. Normalization method by library size was performed according to DESeq2 methods. DEG were determined using a generalized linear model where the counts for each gene on each sample are modeled using a negative binomial distribution with fitted mean and a gene-specific dispersion parameter. The Wald test statistic is applied to test for significance of the generalized linear model coefficients. In the second study (GSE33030), samples were obtained from pregnant or cyclic cows at different time points in the cycle (19 -21) but only samples collected at day 7 from cyclic heifers were considered here. Animals were treated with a progesterone-releasing intravaginal device (PRID) from day 3 of the estrus cycle (HP group, n ϭ 5) or did not receive a device (NP group, n ϭ 5). Extracted RNAs were hybridized to the Affymetrix Bovine Genome Array platform. The raw data were obtained and processed with the gcRMA package (79), which was used to import the raw data into R, perform background correction, and then transform and normalize the data using the quantile normal-ization method. DEG were determined with moderated t statistics, a variation of the classical t test (65).
For the sake of comparisons, significant DEGs in these and our study were defined as those with a P value Ͻ 0.05 and a fold change Ͼ1.4. Up-or downregulated genes in our study were compared with the up-and downregulated genes in the two external studies. The ENSENMBL IDs were used for comparisons. The test of independence (Pearson's 2 test) was employed to determine the relatedness of up-and downregulated DEGs observed in PR vs. non-PR with upand downregulated DEGs in HF vs. LF or HP vs. NP. Similarities between PR and non-PR samples according to the expression of the overlapping genes were assessed through a multidimensional scaling plot (MDS), performed with the Glimma package (70).
Furthermore, we employed an additional approach to confirm the similarities between the endometrial profiles of HF and PR cows. The expression values of DEG (FDR Ͻ 0.05) for the PR vs. non-PR were employed as gene signatures to predict the groups in the HF/LF or HP/LP animals by logistic regression, as applied in other studies (63,74). The reverse approach was also performed, using the expression values of the DEG (FDR Ͻ0.05) between HF vs. LF cows as gene signature. These procedures were executed with the caret package (34). Add-on normalization across the data sets was performed with the bapred package (31).

RNA-Seq Data Preprocessing and Alignment
The sequencing generated, on average, 44,252,601 read pairs per sample, and 41,799,192 of the read pairs (94% of the read pairs) were retained after the preprocessing step. Among them, on average, 91% were uniquely mapped to the bovine reference genome (43% of the read pairs mapped to exonic regions, 27% of the read pairs mapped to intronic regions, and 30% of the read pairs mapped to intergenic regions; Supplemental Fig. S1).

Differential Expression Analysis
We identified 111 DEG (FDR Ͻ 0.05) in PR vs. non-PR cows. Among them, 60 were up-and 51 were downregulated DEG in PR cows compared with non-PR cows (Supplemental Table S1).

Functional Analysis
GSEA. GSEA was applied to identify significantly enriched pathways considering the entire expression pattern change between PR and non-PR cows. The functional enrichment with GSEA determined 11 pathways enriched of genes upregulated in PR cows and four pathways for genes downregulated in PR cows (Supplemental Table S2).
Gene ontology term analysis and functional classification. The NET-GE gene ontology (GO) term enrichment analysis performed on all DEG (up-and downregulated) identified a significant enrichment for two CC GO terms (cell surface and extracellular matrix part), two MF GO terms [leukotriene-C4 synthase activity, type I transforming growth factor (TGF)-␤ receptor binding], and 41 BP GO terms (Supplemental Table  S3). The BP GO terms were mainly related to response to growth factor, negative regulation of TGF-␤ signaling, extracellular matrix organization, and cell adhesion proliferation and migration, as well as to arachidonic acids and eicosanoid metabolic process and immune response and leukocyte proliferation. Other interesting GO terms were response to proges-terone and estrogen. The functional classification provided by IPA revealed involvement of cellular processes (cellular, function, movement, morphology, development, growth, and proliferation free radical scavenging) as well as the importance of cellular metabolism of lipids, carbohydrate, vitamins, and amino acids (Supplemental Table S4).
Upstream regulator analysis with IPA. We selected upstream regulators that were predicted as significantly activated after bias correction by IPA. We identified 11 as being activated and two as being inhibited in PR cows ( Table 2). The top 25 identified upstream regulators are listed in Supplemental Table S5).

Comparison with Other Studies
Results from this analysis are shown in Fig. 2A. There was a significant association between overlapping genes going in the same direction between PR and cows with good genetic merit for fertility, or HF cows (when compared, respectively, to non-PR and LF cows). On the contrary, some genes overlapped in opposite direction between PR and cows with high levels of progesterone (HP) (versus their respective controls), but there was no significant association. The MDS plot (Fig.  2B) depicts the distribution of PR and non-PR samples based on the expression of the overlapping genes. Samples were distributed apart when the expression of the overlapping genes with HF were considered, but not with the expression of the overlapping genes with HP. Using the expression values of DEG between PR and non-PR or between HF and LF cows, animals in the HF/LF groups, or PR/non-PR groups, were correctly classified with an accuracy of 78.5 or 78.9%, respectively (P Ͻ 0.05 for the one-sided test, implying that the accuracy is better than the "no information rate"). Sensitivity was 85.7 and 87.5%, respectively. However, the prediction accuracy was 66.6 and 44.4% (P Ͼ 0.05), respectively, and sensitivity 50%, for cows in the HP/NP groups.
These results reinforce the concept that the observed differences in the endometrial transcriptome between PR and non-PR cows is given due the inherent molecular signature for each cow, rather than to variations in hormonal concentrations.

DISCUSSION
Understanding the physiology of the bovine endometrium and its role in the establishment of pregnancy has been an overarching theme in cattle reproduction. Transcriptomic analyses have become a powerful tool for obtaining in depth molecular knowledge about endometrial function. Exploration of the transcriptome provides an unbiased system-wide analysis of physiologically relevant gene expression levels, determining endometrial receptivity [reviewed by Forde and Lonergan (22)]. It has been demonstrated that the endometrium is already modified at day 5-6 after ovulation, independently of the presence of an embryo, to create a microenvironment to receive the embryo from the oviduct (22). Interestingly, results from this and other studies indicate that the endometrial transcriptome at day 7 of animals that become pregnant differs from those that do not (6,54,60). Therefore, this period is ideal to determine the receptivity of a cow, which is particularly relevant in ET procedures when highly valuable embryos are transferred.
Previous Both studies implemented a similar model used in the present investigation, in the sense that in vivo-produced or IVP embryos, respectively, were transferred in the subsequent cycle Results from these studies contribute greatly to the knowledge of the endometrial transcriptomic profile of the period when an embryo is transferred into a recipient. However, they also demonstrate that the endometrial transcriptome is affected by many factors, including parity (heifers or cows), breeds (dual-purpose, beef or dairy breed), and induction of estrus (72).
Establishment of pregnancy is dependent on a relevant transcriptome conveying endometrial receptivity, but also on embryonic viability. As alluded to earlier, a multitude of factors influence embryo viability, a matter of particular importance with IVP embryos. Here, we aimed at minimizing these sources of variation as much as possible by e.g., using a single bull, a consistent IVP system, strict selection of the best-quality embryos for transfer, fresh transfers, careful matching the day of transfer and stage of embryonic development, as well as a single person performing the ET.
In the present study, lactating dairy cows were observed and monitored for spontaneous estrus, and hence, no hormonal synchronization method was applied. Cows were classified retrospectively according to pregnancy status of IVP embryos transferred on the next cycle after the endometrial biopsies were obtained, which were done with an instrument to include the full endometrium (luminal epithelium, glandular epithelium, stroma, blood vessels, nerves and invading leukocytes). Also, biopsies were taken consistently from the cranial to middle part of the uterine horn ipsilateral to the CL, to investigate specifically the part that receives the embryo at ET Rounded rectangle in gray ϭ secondary pathways and upstream regulators: Hedgehog signaling (i), ribosome (ii), angiogenesis (iii), other upstream regulators (iv). Gray filled rectangles ϭ genes upregulated in PR cows. Rectangles ϭ upstream regulators and biological processes predicted to be upregulated-regulated (orange) or downregulated (blue) in PR cows or with no specific direction (black). Ovals ϭ other molecules from literature. The discontinuous lines point to the biological functions, genes, and upstream regulators belonging to the five main mechanism. Red arrows ϭ positive effect exerted on PR endometrium (embryo development and endometrial receptivity). The black lines points to the three type of cell adhesion molecules (CAD) colored in purple (N-cad, N-cadherins; Ig, immunoglobulins; It, integrins) that are encoded by the upregulated genes. The black arrows show the mechanisms of TNF and TGFb, which bind to the correspondent receptors (pink molecules) in the endometrial surface. and to avoid regional differences in the uterine transcriptome (68). Therefore, our findings could be of enormous help to identify key genes that characterize a receptive endometrium. Results from this and other studies could be used to predict endometrial receptivity of a cow, as it has been developed for humans (12), ideally through the development of measurable biomarkers from accessible samples, such as the blood.
Results from differential gene expression between PR and non-PR cows revealed 111 candidate genes, with 60 upregulated and 51 downregulated in PR cows. Furthermore, the functional analysis of these DEG identified five main biological mechanisms associated with endometrial receptivity (depicted in Fig. 1): 1) ECM interaction, 2) control of histotroph metabolic composition, 3) PG synthesis), 4) TGF-␤ signaling, and 5) inflammation and leukocyte activation. Each of these mechanisms is briefly discussed in the following paragraphs.

ECM Interaction
Involvement of many DEG (upregulated in PR cows) in the ECM was confirmed by the enrichment for MF GO terms (glycosaminoglycan binding) and for the BP GO term (regulation of adhesion, ECM organization, extracellular structure organization). The GO terms belonging to the CC domain revealed that a huge component of the proteins encoded by the DEG are secreted or membrane bound and that many of them are part of the ECM. Among them, 16 proteins showed to function as anchoring junctions important for the endometrial organization. In human, three superfamilies of proteins are involved in ECM adhesion and embryo implantation (43), which were part of the DEG: three immunoglobulins, two cadherins, and two integrins. These findings reflect that the endometrial tissue is rich in ECM molecules for interaction and adhesion. A balanced ECM integrity has been associated with pregnancy maintenance, and the combined action of metalloproteinase inhibitor (TIMP2) and metalloproteinases (75) has been suggested to be responsible for the regulation of ECM integrity (45,75). In addition, genes related with the ECM have been also deregulated in several of the studies mentioned above (6,47,60).

Histotroph Metabolic Composition
The histotroph is the mixture of molecules secreted into the uterine lumen, and, thus, its composition cannot be determined through an endometrial biopsy. However, the endometrial transcriptome reflects pathways that are involved in the cellular production of components of the histotroph. Our transcriptomic results demonstrated alterations related to the histotroph synthesis that could be involved in increased receptivity in the PR cows. A set of genes downregulated in PR cows was enriched for the KEGG pathway valine, leucine, and isoleucine degradation. Amino acids are important components of the histotroph (4) and are essential for conceptus elongation (9,23,25). Bovine embryos utilize amino acids (23,32), and the con- centrations of L-leucine, L-valine, and L-isoleucine increase in pregnant animals (23,26). We speculate that the early inhibition of endometrial genes involved in degradation of these amino acids would lead to a higher concentration of them during the implantation window, improving the endometrial receptivity. In addition, some of the DEG were involved in molecular transport and in the metabolic pathways of carbohydrates, vitamins, minerals, and amino acids, according to the IPA classification. Again, these are important components of the histotroph that show how PR and non-PR cows are different regarding the quality of this essential nourishment for the embryo (45,57).

PG Synthesis
The GO term analysis revealed a set of genes involved in arachidonic acid metabolic processes and leukotriene C4 synthases. We identified PTGES as the upstream regulator of four DEG (CYBB, EZR, PTGS2, and VIM). PTGES is an enzyme that generates PGE from PGH2 (52), which is derived from an intermediate step performed by PTGS2. Interestingly, PTGS2 was differentially upregulated in PR cows. In sheep, PGs generated by the PTGS2 in the conceptus were responsible for regulation of genes involved in conceptus elongation (9).

TGF-␤ Signaling
One of the TGF-␤ signaling pathways is the MAPK cascade that was one of the biological functions enriched with the DEG list (both up-and downregulated in PR cows). TGF-␤1 and its related molecules were identified as upstream regulators of the DEG, and their encoded proteins were predicted to be activated in PR cows. Furthermore, the TGF-␤I and TGFBR2 were upregulated in PR cows. The TGF-␤s control many cellular processes, and the mechanism is exerted through interaction with the ECM. This signaling pathway stimulates ECM formation and prevents its degradation (24,58). Genes belonging to the TGF-␤ pathway regulate preimplantation development in both embryo and endometrium (36). Furthermore, TGF-␤1 and TGF-␤2 stimulate proliferation of bovine trophoblast and fetal endometrial epithelial cell in vitro (37,48). Therefore, activation of the TGF signaling at diestrus could lead to a better embryo receptivity, supporting embryo development and elongation and at the same time stimulating ECM formation.

Inflammation and Leukocyte Activation
The DEG enriched for inflammation response and for leukocyte chemotaxis, adhesion, and proliferation. In cattle, leukocyte presence seems not to be a requirement for placentation (17,51) as in humans (13). However, the innate immune system is activated during early pregnancy in the bovine uterus (77), and several chemokines known to affect immune tolerance are upregulated in pregnant animals. These actions are required for establishment of the right balance between proand anti-inflammatory molecules during initial placentation (77). Immune-related genes were found also to be deregulated in the endometrial samples from day 6 -7 of the estrus cycle in most of the studies mentioned at the beginning of this section (33,47,54,60).
In addition, secondary pathways and upstream regulators previously associated with endometrial receptivity were identified, such as: 1) the Hedgehog pathway (Fig. 1i), constituted by genes upregulated in PR cows, involved in implantation and early embryo development (28,30); 2) the ribosome KEGG pathways (Fig. 1ii), sustained by genes downregulated in PR cows, is likewise downregulated in the endometrium during diestrus (3) corresponding to the phase of initial embryonic development and adhesion; and 3) the angiogenesis process (Fig. 1iii), including the upstream regulator predicted to be activated in PR cows (AGT), which contributes positively to the establishment of endometrium receptivity (46,54).
The upstream regulators of the DEG were explored in the present study as well. Progesterone was identified as the top upstream regulator of the DEG, being associated with the profile of 15 genes. Moreover, response to progesterone and estrogen were enriched processes among the DEG. Lower progesterone concentration has a negative effect on expression of genes that regulate histotroph composition and conceptus elongation (21,22). The upstream regulator analysis did not predict any significant accumulation or depletion of progesterone from the analysis of gene expression changes of its target molecules. However, 11 upstream regulators were predicted to be activated and two to be inhibited in PR cows. Among them, IGF1 (33,59,76), AGT (42,46), TNF (11,50), AKT and PI3K (38,61), FGF2 (68), TGF-␤1 (previously discussed) (activated in PR cows), and NFE2L2 (54) (inhibited in PR cows) have been previously associated with endometrial processes and receptivity.
Given the central role of progesterone and estrogen in defining endometrial receptivity, an additional analysis was performed to characterize the molecular profile of PR and non-PR cows. Two external data sets, obtained from public domains, were reanalyzed to compare the corresponding lists of DEG with the ones from the present study. One of the data sets, introduced earlier in this section, corresponded to endometrial samples obtained at day 7 from lactating dairy cows with good (HF) or poor (LF) genetic merit for fertility. Cows classified as HF presented greater body condition score and circulating IGF1 concentrations throughout lactation, more favorable metabolic status during the transition period and faster uterine recovery, and earlier resumption of cyclicity after parturition (47). The other data set corresponded to a large study with several experiments where samples were obtained from pregnant or cyclic beef cows at different time points in the cycle (19 -21). For the present study, only samples collected at day 7 from cyclic heifers with high (HP) or normal levels of progesterone (NP) were considered. Animals in the HP group were treated with a progesterone-releasing intravaginal device (PRID) from day 3 of the estrus cycle (19). Thus, that experiment aimed to identify the effects of elevated progesterone levels on the endometrial transcriptome.
Comparison of the up-and downregulated DEG on the present study with the DEG from these two studies showed a significant association of genes going in the same direction for PR and HF cows, while no association was observed with DEG in HP cows (Fig. 2). In addition, the expression values of DEG for the PR/non-PR cows were able to accurately classify HF/LF cows and vice versa, but not HP/NP animals. These results suggest that the endometrial transcriptome of PR cows is more comparable with the one from HF cows than with the one influenced by high progesterone levels, although intrinsic variations between the endometria of dairy and beef cows could be affecting the results as well. Nevertheless, inherent characteristics of PR and non-PR cows are probably defining the differences in the endometrial transcriptome, rather than hormonal variations between these two groups, since other authors have reported similar circulating progesterone levels between fertile and infertile cows (33,44). Therefore, we can speculate that the endometrial transcriptome of PR cows elicits a different response than non-PR cows to the same progesterone concentration. In other words, progesterone, as an upstream regulator, could induce the expression of genes favoring pregnancy in PR cows, but such gene expression is not elicited at the same favorable level in the non-PR animals.

Conclusion
Findings from this study provide knowledge about the endometrial molecular signature required for pregnancy establishment, around the time of embryo transfer, in lactating dairy cows. Five main relevant mechanisms were identified: ECM interaction, histotroph metabolic composition, prostaglandin synthesis, TGF-␤ signaling and inflammation, and leukocyte activation; all of which appear to be regulated by progesterone and estrogen concentrations. Moreover, we identified candidate upstream regulators that were predicted to be activated (IGF1, AGT, TNF, AKT, PI3K, FGF2, TGF-␤1) or inhibited (NFE2L2, alpha catenin) in PR cows. In addition to our own experimental data sets, we analyzed two public gene expression data sets from similar experiments and compared up-and downregulated DEG from our study with the DEG from these two studies. Results confirmed significant association of genes going in the same direction for highly receptive cows. Hence, these candidate gene regulators and enriched pathways could be used for identification of potential biomarkers for predicting pregnancy success in IVP-ET in artificial reproduction technology. However, further experimental validation in large samples is recommended.