Impact of negative energy balance on transcriptomic profiles of three endometrial cell types isolated by laser capture microdissection in postpartum dairy cows


 Background: In postpartum dairy cows, the energy needs to satisfy high milk production induces a more or less pronounced Negative Energy Balance (NEB) status. NEB associated with fat mobilization impairs reproductive function. This study investigated the specific impact of NEB on gene expression in the three main types of endometrial cells at time planned for insemination and implantation. Endometrial cell types (stromal, glandular and luminal epithelial cells) were isolated by laser micro-dissection allowing the study of constitutive gene expression and their specific response to NEB. Methods: Nine Swedish Red cows receiving a control diet or a mild restricted diet to induce differences of energy balance were categorized into mild (MNEB, n = 5) and severe negative energy balance (SNEB, n = 4). The three endometrial cell types: luminal (LE), glandular (GE) epithelium and stroma (ST) were collected by laser microdissection from endometrial biopsies performed at 80 days postpartum. Results: Transcriptome profiles obtained by RNA sequencing revealed differences in constitutive gene expression between the three cells types and also differences in specific responses related to the severity of NEB. Number of differentially expressed genes between SNEB and MNEB cows was higher in ST than in LE and GE, respectively. SNEB was associated with differential expression of genes related to metabolic processes and embryo-maternal interactions in ST. Under-expression of genes related to cell structure was found in GE whereas genes related to pro-inflammatory pathways were over-expressed. Genes associated to adaptive immunity were under-expressed in LE. Conclusion: The three different main cells types of the endometrium, have very different patterns of gene expression. The severity of NEB after calving is associated with changes in gene expression at time of breeding. Specific alterations in GEs are associated with activation of pro-inflammatory mechanisms. Concomitantly, changes in the expression of genes related to cell to cell interactions and maternal recognition of pregnancy takes place in ST. The combination of these effects possibly altering the uterine environment and embryo maternal interactions may negatively influence the establishment of pregnancy.

play key roles for the establishment and maintenance of pregnancy through activation of the innate 1 immune system and secretion of chemokines [19] that support the recruitment and activation of 2 immune cells directed against pathogens. Moreover, LE and GE exhibit unique molecular signatures 3 having cooperative roles at time of establishment of pregnancy [16,20,21]. Their morphology [22] 4 and biochemical activity [23] differs at time of implantation. RNA-sequencing of the complete 5 transcriptome for the three cell types has been described for equine cells [24]. Laser capture 6 microdissection (LCM) has also been successfully used to retrieve two different uterine epithelial cell 7 types to define the transcriptome and proteomic analysis of the ovine and porcine endometrium, 8 respectively [25,26]. However, to our knowledge the transcriptomic profile of bovine endometrial 9 cells has not yet been documented. Previously published research, regarding the impact of NEB on 10 uterine function and endometrial transcriptome, suggests that NEB associated with elevated non-11 esterified fatty acids (NEFAs) concentrations induces infertility in postpartum cows through 12 dysregulation of immune pathways [12]. However, the understanding of molecular changes induced 13 by NEB from entire endometrial tissues is still unclear and difficult to interpret functionally as 14 responses may be affected by other cell types such as endothelial cells, smooth muscle cells and 15 leukocytes [27]. In vitro studies have clearly shown that NEFAs stimulate pro-inflammatory cytokine 16 production and lipid accumulation of endometrial cells [28] and oviductal epithelial cells [29] but the 17 results from these in vitro models need to be confirmed in vivo. 18 We hypothesized here that NEB may differentially influence the physiology of three endometrial cell 19 types. The objectives of the present study were i) to investigate transcriptomic profiles of luminal 20 epithelial cells, glandular epithelial cells and stromal cells which were harvested by LCM, and ii) 21 identify possible differences in the profiles between cows diagnosed with either mild or severe NEB 22 during the postpartum period. The collection of endometrial biopsies was performed at time of 23 planned AI and the observed changes in gene expression suggest the existence of long-term impacts 24 of NEB that are cell type-specific. given in sheets S8, S9 and S10 respectively of the additional file 12 (TableS4_TS8_TS9_TS10_DEG-SNEBvsMNEB.xlsx). In SNEB animals, a large proportion of 13 DEGs were identified as over-expressed in ST (72%) and GE (63%) whereas almost all DEGs were 14 under-expressed in LE (98%) ( Table 1). An overview of the differential patterns of gene expression in 15 ST, GE, and LE obtained by LCM between SNEB and MNEB cows are illustrated in volcano plots 16 ( Figure 5A to 5C). 17

14
Functional classification using DAVID identifies also a first cluster of nine genes encoding proteins 1 including mainly G-protein coupled receptors (GO: 0005887; integral component of plasma 2 membrane), which were under-expressed in ST from SNEB. Six genes encoding membrane proteins 3 with immunoglobulin-like domains and related to cytokine are part of a second cluster and a last 4 group includes 28 genes coding for component of membrane. 5

6
Differential expression in GE (Table 4, Table 5 and additional TableS7_david_GE-7 overexpressed.pdf) 8 Only seven known genes are under-expressed in GE cells from SNEB cows when compared to 9    6 Differential expression in LE (Table 6 and additional TableS8_david_LE-underexpressed.pdf) 7 In LE samples, only B4GALT5 is over-expressed in SNEB. No significant enriched GO terms is 8 related to the 55 under-expressed DEGs at FDR p value <0.05. By raising the FDR p value at 0. 25,9 over-represented DEGs corresponds to biological processes associated with complement activation, B 10 cell mediated immunity, defense response to bacterium, cell differentiation and cellular component 11 link with plasma membrane and organelle.

KEGG pathway analysis of the DEGs. 4
Significantly enriched KEGG pathways from DAVID database were found in GE and ST, whereas no 5 significant KEGG pathway was detected in LE. In ST cells, DEGs between SNEB and MNEB cows 6 were significantly enriched in four different KEGG pathways. 25 KEGG pathways were recognized 7 by David with the overexpressed genes. Two were found significantly enriched. They are related to 8 calcium signalling pathway (KEGG map04020, fold enrichment = 3.4; 17 DEGs) and tight junctions 9 (KEGG map04530; fold enrichment = 4.8; 11DEGs. With under-expressed DEGs, two KEGG 10 pathways associated with viral infectious diseases (KEGG "measles" map05162 and KEGG 11 "Influenza A" map05164; fold enrichment respectively = 5.0 and 4.1; 11 DEGs) are overrepresented 12 (Table 7). The names of these two KEGG pathways do not make sense with endometrial physiology. 13 The genes of these pathways are known to be important partners of interferon signalling that is a 14 critical mechanism for establishment of pregnancy (reactome pathways: BTA-913531, BTA-877312 The corresponding STRING-generated interaction network obtained from DEGs belonging to the 5 1 KEGG pathways associated to ST and GE cells revealed strong interactions (PPI enrichment value < 2 1.0E-16) between these sets of DEGs that are related to the JAK/STAT signalling ( Figure 6). 3 DISCUSSION 4 During negative energy balance (NEB), lipolysis in adipose tissue is increased resulting in decreased 5 BCS and increased NEFAs in blood [31]. Changes in BCS and NEFA concentrations were correlated 6 with NEB nadir and plasma NEFA concentrations in SNEB cows were greater than in MNEB cows in 7 the prepartum and early post-partum. Both observations are consistent with earlier findings [32] and 8 shows that the two groups were in a different metabolic status before and during the two first weeks [40]. This may result from differences between species but could also reveal the changes induced by 24 the conceptus on the endometrial transcriptome previously reported from full tissue [37,41] and 25 epithelial cells [26]. 26 Using a cut-off of 10 TPM, different numbers of genes were expressed in the three endometrial cell 1 types. ST expressed 5% and 25% more genes than GE and LE, respectively. However, as reported 2 before from a large variety of tissues [42], and the three laser-dissected cell types of porcine 3 endometrium [43], our results confirm that a high number of genes are expressed in common in 4 different endometrial cell types. In the present study, 70 to 85% of genes were expressed in all cells 5 suggesting either "house-keeping" functions or genes encoding proteins with functions common to the 6 endometrium while lower proportions (5%, 7% and 15% for LE, GE and ST, respectively) were 7 restricted to each cell type indicating that they code for proteins supporting the functional specialized 8 signature of each cell type. When compared to porcine endometrium [43], the number of genes 9 showing cell-specific expression is in the same order of magnitude for GE and LE, but appears 10 different for ST cells where this number is ten times higher. These differences in specific expression 11 between cell types, especially the large number of functions enriched in ST are well reflected by the 12 REVIGO analysis (Figure 4). 13 Regardless of the cut off chosen and related limitations, these studies illustrate huge differences in the 14 gene expression patterns between cell types corresponding to specialized functions. This confirms that 15 separating cell types is more appropriate and possibly less biased to decipher the impacts of any factor 16 on a given tissue than former approaches based on full tissue. The clear clustering obtained when 17 analysing the full transcriptome, indicates that luminal and glandular epithelial cells are closely 18 related. These similarities may reflect common functional properties and/or may be related to the 19 common epithelial nature of these cells. The genes associated with GE and LE, which distinguish 20 these two epithelial cells from stromal cells, are all related to GO terms typical of epithelia (GO: 21 0030855, epithelial cell differentiation; GO: 0060429, epithelium development; GO: 0045216, cell-22 cell junction organization). Examples are given below for critical genes previously cited as key 23 regulators of endometrial epithelial cells. CDH1 is involved in organization of epithelium in mouse 24 and its ablation causes the absence of endometrial glands. Occludin is an important protein for tight 25 junction assembly which preserves the epithelial barrier function. The REVIGO analysis showed that 26 in both epithelial cell types, genes encoding proteins related to metabolism were under-represented. 27 On the contrary, genes related to cilium function are enriched in GE, whereas those involved in 1 binding/ receptor function and adhesion are over-represented in LE. In addition, LE cells differentiate 2 from GE by the expression of genes like CYP26A1, that encodes a key enzyme of trans retinoic acid 3 inactivation, already shown as strongly expressed in luminal epithelial cells of rat endometrium and 4 playing a role in embryo implantation [44]. Endometrial expression of HCTR1 has been reported to 5 be, with its main ligand orexin-A, an important local regulator of endometrial functions in porcine 6 uterus [45,46]. 7 As mentioned above for GE and LE, the REVIGO analysis showed that genes involved in metabolism 8 were also under-represented in ST. In contrast, a very large number of functions including but not 9 limited to, cell structure, angiogenesis, extra cellular matrix and immunity are enriched in ST whereas 10 a lack of strong expression of these genes is observed in GE and LE. As awaited, among the genes 11 most discriminating stromal cells, those involved in the production of extracellular matrix and 12 PRLEP gene expression has been reported to be regulated by the embryo in the bovine oviduct [56]. 24 Contrary to the porcine endometrium where its expression was located in epithelial cells, NTRK2 was 25 mainly expressed here in stromal cells [57]. The expression of the NTRK2 gene, which encodes the 26 receptor of brain derived neurotrophic factor, is conserved in mammalian uterus but its signalling function is not yet understood in the female reproductive system [58]. Genes known to be key 1 regulators of uterine receptivity in different species such as, HOXA10 and HOXA11 belong also to the 2 Finally, among these first 50 genes that best separate ST from epithelial cells, numerous ones have not 9 been described so far in the mammalian endometrium. For instance, we could not find any

Impact of NEB on the three endometrial cell types 15
Overall, our results show that NEB impacts mainly ST whereas GE and LE cells are less affected. 16 More than 10% (13%) of the total number of genes expressed in ST were impaired by NEB status 17 while less than 1% were affected in GE and LE (0.3% and 0.7% respectively). When considering the 18 sub groups of genes showing a specific expression related to cell type, NEB did not affect any of differences between studies, however due to the lack of effect on GE, these discrepancies may result 2 also from differences in time postpartum and severity of NEB. The glandular epithelium plays a major 3 role in the activation of the innate immune system as reviewed by [87]. In our study, most of the 4 DEGs in GE related to chemokines, immune response processes, TLRs and TNF signalling pathways, 5 such as CCL2, CCL3, CCL4, CCL11, FOS, JUNB, and SOCS3 were strongly over-expressed in SNEB 6 cows. Some of those genes belonged to the C-C motif chemokine ligands (CCLs) family and play an 7 important role in monocyte recruitment in the endometrium [88]. Increased expression of CCL2 8 mRNA was found associated with lipid accumulation induced uterine inflammation in obese rats [89]. 9 The present results are similar with previous studies performed with full endometrial tissue, showing 10 the up-regulation of inflammatory response genes in SNEB cows [38]. This is also consistent with 11 several studies in mammals showing that metabolic imbalance, increased lipolysis and most 12 particularly NEFAs, play essential functions in the activation of TNF and TLRs signalling to promote 13 the release of pro-inflammatory molecules [90,91]. Taken together, these studies and our present 14 findings suggest that SNEB and NEFAs activate pro-inflammatory pathways in the glandular 15 epithelium and stromal cells. On the contrary, in luminal epithelium, the adaptive immune response 16 (B cell-mediated immunity) and innate immunity, was represented by under-expressed genes such as 17 tracheal antimicrobial peptide (TAP), a beta-defensin gene, which was associated to the NF-κB 18 pathway [92], and by genes coding for immunoglobulin heavy variable chains that participates in the 19 antigen recognition. These observations need further confirmation. Our results indicate that SNEB 20 induces changes in immune responses, which are different in the three endometrial cell types. They 21 show also that these changes are still present, long after NEB has disappeared suggesting long term 22 effects of metabolic imbalance and NEFAs on the pro-inflammatory status of the glandular epithelium 23 and the stroma. 24

Effect of NEB on genes related to maternal-conceptus recognition. A large set of IFN-inducible 25
genes such as MX1, MX2, STAT1, JAK1, IFIH1, IFNGR2, ISG15, LY6G6C, OAS1Y, OAS1Z and IRF7 26 were under-expressed in ST of SNEB cows. A weaker expression of those genes that encode proteins 27 involved in IFN-T signalling could account for the decreased endometrium-related fertility in SNEB 1 cows. In pregnant ruminants, IFN-T is the main pregnancy recognition signal [93], that allows the 2 persistence of the corpus luteum and maintaining elevated progesterone concentrations by blocking 3 oxytocin signalling and PGF2 secretion [94]. Oxytocin signalling has been associated with the AREG was over-expressed in GE of SNEB cows. AREG gene is known as an epidermal growth factor 2 receptor and is involved in cell growth, proliferation, differentiation and migration. It is highly 3 expressed in luminal and glandular epithelium during the secretory phase of menstrual cycle and early 4 pregnancy in human and primate [102]. As for SOCS3, it could be speculated that the over-expression 5 of AREG mRNAs in GE may be part of a compensatory mechanisms in response to the increased 6 expression of cytokines in these cells. It would be interesting to compare the amplitude of over-7 expression of SOCS3 and AREG in the present situation (luteal phase under cyclic conditions) and in 8 pregnancy to evaluate possible impacts of NEB on implantation. 9

10
The present study provides novel and specific information about gene expression in three endometrial 11 cell types from postpartum dairy cows and illustrates specific signatures in ST, LE and GE cells. We

Ethics approval and consent to participate 3
All experimental protocols were approved by the Uppsala Animal Experiment Ethics Board 4 (application C329/12, PROLIFIC)(Uppsala University, Sweden) and were carried out in accordance 5 with the terms of the Swedish Animal Welfare Act. After the study was completed, all cows were kept 6 alive under normal husbandry conditions. 7

Consent for publication 8
« Not applicable » 9

Availability of data and materials 10
The data will be deposited pending acceptance of publication. The datasets generated and/or analyzed 11 during the current study will be available in the [NCBI/Gene expression omnibus)] repository, 12 [https://www.ncbi.nlm.nih.gov/geo/info/seq.html] 13

Competing interests 14
The authors declare that they have no competing interests. 15

Funding 16
This work was funded by from European Union [FP7-KBBE-2012-6, Prolific (Pluridisciplinary study  17 for a RObust and sustainanLe Improvement of the Fertility In Cows,) Grant agreement number 18 311776) for animals, reagents, sequencing and travel. WC was supported by the Rajamangala 19 University of Technology Srivijaya (RMUTSV), Thailand. The funders had no role in study design, 20 data collection, and interpretation, or the decision to submit the work for publication. 21

Authors' contributions 22
W.C., P.H. and G.C. contributed to the conception and design of the study. S.L., M.R., C.R., C.B. and 23 T.N. contributed to sample collection and preparation. G.C., D.M., Y.G., W.C. and P.H. performed 24 bioinformatics analysis and integration of data. W.C. performed the experiment, sample collection 25 and preparation, data analyses and W.C., G.C., and P.H. drafted the manuscript. All authors provided 26 critical feedback and helped shape research, analyses and manuscript. GC and PH are both senior co-27 authorship. 28

Acknowledgements 29
We acknowledge Pierrette Reinaud (INRAE, Jouy en Josas, France) and Olivier Dubois for consulting 30 during the LCM process and RNA analysis. The authors would like to thank the staff of Swedish 31 Research Center, Lövsta, Uppsala, Sweden and Biology of Reproduction, Epigenetic, Environment 32 and Development (BREED), INRA, Jouy en Josas, France for their help and support. 33

Authors' information (optional) 34
Methods 1 Animals and experimental design. This study was approved by the Uppsala Animal Experiment Ethics 2 Board (application C329/12, PROLIFIC). After the study was conducted all cows have been kept in 3 usual farm living conditions. The animals used in this study were second lactation cows of the 4 Swedish Red breed (SRB; n = 12) fed two different diets i.e. i) high-energy diet (control, n=6) 5 targeting 35 kg energy-corrected milk (ECM) and ii) low-energy diet targeting (n=6) 25 kg energy-6 corrected milk (ECM) which was achieved by giving to these cows 50% concentrate. All cows were 7 conducted at the Swedish Livestock Research Centre in Lӧvsta, Uppsala, Sweden. For each cow, the 8 differential diets were given between 30 days prepartum and 120 days postpartum. The animals were 9 kept in a loose housing barn with a voluntary milking system (VMS, DeLaval, Tumba, Sweden), and 10 had free access to drinking water. The dietary details and management conditions were previously 11 described [32]. During the experiment, consumption of concentrate was individually adjusted with an 12 automatic feeding machine while forage was fed ad libitum. At day 60 after calving, estrous was 13 synchronized using an intra-vaginal progesterone device (CIDR, Zoetis, Parsippany, NJ, USA) for a 14 week followed by i.m. injection of 500 µg of prostaglandin analog (Estrumate ® , MSD animal health, 15 Madison, NJ, USA) intramuscular as described [103]. Fifteen days after visual oestrus detection, 16 endometrial tissue biopsies were collected under epidural anesthesia with 0.5 mg/kg of 1% lidocaine 17 hydrochloride (1% Xylocaine ® , Astra Zeneca, Cambridge, UK). Timeline for samplings and analysis 18 of phenotypic responses are presented in supplemental Figure.S1. 19 Energy balance (EB) calculation and classification. The energy balance (EB) (residual feed intake 20 (RFI) expressed in MJ/day) was calculated as the difference between energy consumed and energy 21 used for milk production, body maintenance, growth and pregnancy for each individual cow. 22 Calculations were performed once per week from first week after calving to day 120 as described in 23 [104]. All data used were routinely recorded in the university herd and energy balance calculation was 24 performed with NorFor used as the reference system in the Nordic countries. Based on most 25 differentiated EB profiles, nine out of twelve cows were classified into two NEB groups with either a 26 mild negative energy balance (MNEB) group (n = 5) or a severe negative energy balance (SNEB) 27 group (n = 4). Residual feed intake values in the first week postpartum of these nine cows ranged 1 from -52.77 to 21.26 MJ/day and means (± s.e.m.) of 1.30 ± 6.35 and -29.48 ± 7.10 MJ/day were 2 observed in the MNEB and SNEB groups, respectively. 3 Body condition score (BCS) and plasma NEFA measurements. Body condition score (BCS) was 4 evaluated and recorded by the same person every two weeks, from 30 days prepartum until 120 days 5 postpartum. BCS was used on a 5 point scale with 0.5 point increments, 1 = very lean to 5 = fat [105]. 6 Blood samples were taken every two weeks from the coccygeal vein in EDTA containing tubes (BD 7 Vacutainer, Kremsmünster, Austria) from 30 days prepartum to 56 days postpartum and then 8 centrifuged at 4000 g for 10 min at 4°C. Following centrifugation, plasma samples were distributed 9 into 0.5 mL aliquots and stored in -20° C until NEFA analyses were performed. NEFA concentrations 10 were measured in duplicate by using a non-esterified fatty assay kit (Bio Scientific Corporation, 11 Austin, TX, USA) with detection range 0 -4 mM. The intra-and inter-assay variability was 4.19 ± 12 3.99% and 2.63 ± 1.08%, respectively. UK) as previously published [32]. The progesterone concentration profile was used to determine the 18 estrous cycle stage at the time of biopsy sampling. All cows selected were in the luteal phase at time 19 of endometrial biopsy as shown by their mean (± s.e.m.) progesterone concentration (8.78 ± 2.12 20 ng/mL; range from 6.66 to 10.90 ng/ml). 21 Collection of endometrial biopsies. Endometrial biopsies were collected from the uterine horn 22 ipsilateral to the corpus luteum by using Kevorkian-Younge uterine biopsy forceps (Alcyon, Paris, 23 France). Biopsies were cut into 3 pieces (sizes ≈ 4×4 mm 2 ). One of them was snap frozen in cold 24 isopentane (2-Methylbutane, Sigma Aldrich, Saint Louis, MO, USA) previously placed in liquid 25 nitrogen for 5 min, and immediately embedded in ≈1 cm 3 optimal cutting temperature (OCT) 26 compound (VWR, Radnor, PA, USA). OCT conditioned biopsies were then put into dry ice and kept 1 at -80°C until sectioning. Tissue blocks were 8 µm sectioned with a cryostat (Leica CM1860 Cryostat, 2 Wetzlar, Germany) at -20°C under RNA-free conditions. Tissue section slices were mounted on Super 3 Frost slides RNA-free which were chilled on ice, following immersion in ice-cold 75% RNA-free 4 ethanol and stored at -80°C until staining [106]. 5 Laser capture microdissection (LCM) and RNA isolation. All procedures used were those previously 6 published [106]. Tissue sections were mounted on RNAse-free glass slides which were chilled on ice, 7 following immersion in cold 75% RNA-free ethanol at -20°C in the cryostat and then transferred into 8 75% ethanol at RT (30 sec), stained with 1% cresyl violet in ethanol (15 sec), rinsed successively with 9 75% ethanol (30 sec), 95% ethanol (2 x 1 min), and 100% ethanol (2 x 1 min) (anhydrous Ethanol 10 absolute). Finally, the slides were completely dehydrated by immersion in pure xylene (M-xylene, 11 Sigma-Aldrich, Saint-Quentin-Fallavier, France) for 2 × 5 min. Stained tissue sections were then 12 immediately air dried. The LCM process was performed by using an ArcturusXT™ Laser Capture 13 Microdissection System and software (Applied Biosystems ® , Arcturus, ThermoFisher Scietific,  Table S9). 2 RNA sequencing and data analysis. RNA sequencing libraries prepared from 24 samples (number of 3 samples in each NEB group and endometrial cell types presented in Table S9) were prepared and 4 sequenced on GenomEast Platform (IGBMC, Cedex, France; http://genomeast.igbmc.fr/). Libraries 5 were built using the Clontech SMART-Seq v4 Ultra Low Input RNA kit for Sequencing. Full length 6 cDNA were generated from 4 ng of total RNA using Clontech SMART-Seq v4 Ultra Low Input RNA 7 kit for Sequencing (Takara Bio Europe, Ozyme, Montigny-Le-Bretonneux, France) according to 8 manufacturer's instructions, with 10 cycles of PCR for cDNA amplification by Seq-Amp polymerase. 9 Then, 600 pg of pre-amplified cDNA were then used as input for Tn5 transposon tagmentation using 10 the Nextera XT DNA Library Preparation Kit (Illumina, San Diego, CA) followed by 12 cycles of 11 library amplification. Following purification with Agencourt AMPure XP beads (Beckman-Coulter, 12 Roissy, France), the size and concentration of libraries were assessed by capillary 13 electrophoresis. Sequencing was performed on an Illumina HiSeq 4000 with 50 bp paired-end reads. 14 Image analysis and base calling were performed using RTA 2.7.3 and bcl2fastq 2.17.1.14. Gene level 15 exploratory analysis and differential expression were performed using the RNAseq workflow 16  20). DEGs of specific-endometrial cell samples were identified in comparison between 1 SNEB and MNEB group with an adjusted p-value of 0.05. Volcano plot was applied to gene lists of 2 each endometrial cell type considering the log2 fold change between SNEB and MNEB on the x axis 3 and the negative log10 of the adjusted p-value on the y axis. 4 Gene ontology and KEGG Pathway Analysis. Lists of genes expressed by the three types of 5 endometrial cells as well as sets of over-or under-expressed DEGs between SNEB and MNEB were 6 annotated into three categories of Gene Ontology (GO) pathways such as biological process (BP), 7 cellular component (CC) and molecular function (MP) using PANTHER classification system 8 (Protein Analysis THrough Evolutionary Relationships version 14.0, http://pantherdb.org). 9 PANTHER overrepresentation tests were performed using all genes from the whole Bos taurus 10 genome or from specified list. Lists of GO terms were summarized and visualized in semantic space 11 by REVIGO (http://revigo.irb.hr/) [111]. The SimRel semantic similarity score was used and the 12 threshold was set at 0.15. Moreover, the analysis of enriched Kyoto Encyclopedia of Genes and 13 Genomes (KEGG) pathways was performed using Database for Annotation, Visualization and 14 Integrated Discovery software (DAVID version 6.8, https://david.ncifcrf.gov/summary.jsp). If a 15 KEGG pathway was determined to be significantly enriched (Benjamini-adjusted p-value < 0.05), 16 this significant process/pathway was reported. By using DEGs which are involved in significant 17 KEGG pathways, a molecular interaction network analysis was generated by using STRING database 18

Statistical analysis 22
The statistical analyzes for phenotype parameters (BCS, NEFA concentrations) were performed using 23 the Statistical Analysis System Software (SAS ® version 9.4, SAS Institute Inc., Cary, NC, USA) and 24 analyzed by mixed models with repeated measurement (Proc MIXED). All variables were checked for 25 normality and data were log10-transformed if needed. The effect of the cow was considered as 26 random when running the models. The model included NEB group, diet group and time of sampling 1 defined as fixed effects and their second order interactions. Non-significant effects were progressively 2 removed from models. Scheffe's post hoc test was used for multiple comparisons and also the 3 "estimate" and "contrast" statements under Proc MIXED were used for pairwise comparisons. 4 Individual BCS loss from start of experiment until a nadir (the lowest postpartum value) and a nadir 5 of feed residual intake (RFI) value after calving were recorded. Pearson correlation coefficients 6 between the different variables were calculated using the Proc CORR function. The results of BCS, 7 NEFA's concentration, and milk progesterone concentration are presented as LSmeans ± S.E.M. 8 Differences with associated p-value < 0.05 were considered to be significant. In the statistical analysis 9 of transcriptome profiles, generalized linear model was fitted and Wald test were performed to 10 determine which of the observed fold changes were significantly different between severe and mild 11 negative energy balance groups. p-values < 0.05 were considered to identify DEGs according to 12   ST and GE endometrial cell types selected from significant KEGG pathways (Table 8)