# Maternal age affects equine day 8 embryo gene expression both in trophoblast and inner cell mass

Emilie Derisoud, Luc Jouneau, Cédric Dubois, Catherine Archilla, Yan Jaszczyszyn, Rachel Legendre, Nathalie N. Daniel, Nathalie N. Peynot, Michèle Dahirel, Juliette Auclair-Ronzaud, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Emilie Derisoud, Luc Jouneau, Cédric Dubois, Catherine Archilla, Yan Jaszczyszyn, et al. Maternal age affects equine day 8 embryo gene expression both in trophoblast and inner cell mass. BMC Genomics, 2022, 23 (1), pp.443. 10.1186/s12864-022-08593-7. hal-03695686v2

HAL Id: hal-03695686

https://hal.inrae.fr/hal-03695686v2

Submitted on 15 Jun 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>
RESEARCH

Open Access

<figure><img src="image_3.png" /></figure>

Check for
updates

# Maternal age affects equine day 8 embryo gene expression both in trophoblast and inner cell mass

Emilie Derisoud<sup>1,2*</sup>, Luc Jouneau<sup>1,2</sup>, Cédric Dubois<sup>3</sup>, Catherine Archilla<sup>1,2</sup>, Yan Jaszczyszyn<sup>4</sup>, Rachel Legendre<sup>5</sup>, Nathalie Daniel<sup>1,2</sup>, Nathalie Peynot<sup>1,2</sup>, Michèle Dahirel<sup>1,2</sup>, Juliette Auclair-Ronzaud<sup>3</sup>, Laurence Wimel<sup>3</sup>, Véronique Duranthon<sup>1,2</sup> and Pascale Chavatte-Palmer<sup>1,2*</sup>

**Abstract**

**Background:** Breeding a mare until she is not fertile or even until her death is common in equine industry but the fertility decreases as the mare age increases. Embryo loss due to reduced embryo quality is partly accountable for this observation. Here, the effect of mare's age on blastocysts' gene expression was explored. Day 8 post-ovulation embryos were collected from multiparous young (YM, 6-year-old, *N* = 5) and older (OM, > 10-year-old, *N* = 6) non-nursing Saddlebred mares, inseminated with the semen of one stallion. Pure or inner cell mass (ICM) enriched trophoblast, obtained by embryo bisection, were RNA sequenced. Deconvolution algorithm was used to discriminate gene expression in the ICM from that in the trophoblast. Differential expression was analyzed with embryo sex and diameter as cofactors. Functional annotation and classification of differentially expressed genes and gene set enrichment analysis were also performed.

**Results:** Maternal aging did not affect embryo recovery rate, embryo diameter nor total RNA quantity. In both compartments, the expression of genes involved in mitochondria and protein metabolism were disturbed by maternal age, although more genes were affected in the ICM. Mitosis, signaling and adhesion pathways and embryo development were decreased in the ICM of embryos from old mares. In trophoblast, ion movement pathways were affected.

**Conclusions:** This is the first study showing that maternal age affects gene expression in the equine blastocyst, demonstrating significant effects as early as 10 years of age. These perturbations may affect further embryo development and contribute to decreased fertility due to aging.

**Keywords:** Horse, Blastocyst, Deconvolution, RNA-sequencing, Mare

**Background**

Unlike other farm animals, the main purpose for the use of horses is athletic performance and leisure. In many cases, mares produce foals regularly to remain profitable after their sport career so that their reproductive career

may last until advanced age, even though mares' fertility declines with age (for review [1]). Indeed, in Thoroughbreds, pregnancy rates per estrous at Day 42 post-ovulation were shown to be reduced from 63.0% in young mares to 58.5% in mares older than 9 years of age, and to reach 43.8% in mares older than 18 years old [2]. It is difficult to determine an age limit defining when a mare is considered too old for breeding as aging is a progressive process during the life. One study on 535,746 fertility records in mixed breeds in France, however, suggested

<footnote>*Correspondence: emilie.derisoud@gmail.com; pascale.chavatte-palmer@inrae.fr</footnote>

<footnote n="2">² Ecole Nationale Vétérinaire d'Alfort, BREED, 94700 Maisons-Alfort, France
Full list of author information is available at the end of the article</footnote>
that if a 10% of fertility decline is chosen as cutoff value for retiring breeding mares, then mares older than 8 should not be bred [3].

Embryo loss seems to be the principal cause for the fertility decline in old mares. Using pregnancy data obtained from 898 cycles in Thoroughbred mares over period of 11 years, the cutoff for embryo mortality to reach values above average was calculated to be 10 years of age [4]. In practice, embryo loss frequency is generally observed between the first (14–22 days) and the second (42–45 days) pregnancy check in old than in young mares [2, 4, 5]. Embryo loss may, however, be unobserved in earlier stages. Indeed, although experimental oviductal flushes on Day 2 post-ovulation showed that the fertilization rate was similar between young and fertile mares vs old and subfertile mares (respectively, 10/14 vs 11/14 embryos recovered per group on Day 2), when the same mares were inseminated during another estrus period, the pregnancy rate was significantly higher at 14 days post ovulation in young and fertile mares compared to older ones (respectively, 12/15 vs 3/15 pregnant mares on Day 14 by ultrasonography) [6].

Although uterine degeneration, including large inflammatory cell infiltration, fibrosis and reduced gland density [7–9] is observed in old mares, oocyte quality seems to play a major role in the observed age-related fertility decrease: the transfer of Day 4 embryos collected from old and subfertile mares resulted in lower pregnancy rates at 14 days post-ovulation than when embryos collected from young and fertile mares were transferred. This remained true when only the subset of good quality embryos was considered [10].

Controversial effects of maternal age on oocyte morphology have been observed in several studies [11–15] but there is a clear effect of maternal age on oocyte quality and function. Old mares’ oocytes are less likely to reach metaphase II after in vitro maturation than oocytes from younger mares [11]. Alterations of the expression of several genes have been also observed during in vivo and in vitro oocyte maturation [16, 17]. Moreover, after in vitro maturation, higher spindle instability is observed in old versus young mares’ oocytes, resulting in increased risks of chromosome misalignment leading to embryo aneuploidy [18–20]. Controversial effects of mare age on mitochondria degeneration have been also reported [16, 21–23].

The effect of maternal age on embryo size is not so clear, as studies demonstrated a reduction of embryo diameter in old mares at several developmental stages [24–28] while others did not detect any effect [9, 13, 25, 29] of maternal aging on embryo size. In these studies, however, the developmental stage was not uniform and the method to measure embryo size could be different

(direct measurement or using ultrasonography). Furthermore, the size of the equine embryo is highly variable for the same gestational age [30]. The only study that directly measured more than 200 Day 7 to 10 embryos as part of an embryo transfer program observed that maternal aging induced a reduction in embryo size (in average, more than 740 µm in diameter for embryos produced by mares younger than 15 vs less than 620 µm in diameter for embryos produced by mares older than 16) [31].

Maternal age, however, seems to affect embryo morphology from the earliest stages of development, before the embryonic genome activation started (around the 8–16 cells stage) [32, 33]. Zygotes collected in vivo on Day 1.5 post ovulation, indeed, had the same number of cells, but a poorer morphology grade (uniformity of cell size and texture, imperfections and abnormalities) when they had been collected from mares older than 15 years of age compared to young ones [34]. In addition, at Day 3 post ovulation, zygotes collected from mares older than 20 years of age combined poorer morphology grade and reduced cell number [34]. In contrast, in other studies, no difference according to maternal age was observed in terms of number, size or texture of blastomeres of embryos collected in oviducts on Days 2 and 4 post ovulation [6, 10]. Finally, at blastocyst stage, a higher incidence of morphological abnormalities was observed in embryos that were collected from subfertile and/or 18 years old mares compared to 4 years old nulliparous mares [24].

Only a few studies have examined the effect of maternal age on equine embryo quality. As observed in oocytes, metabolic function and capacity were reduced in Day 2 embryos collected from old mares (≥20 years old) compared to those collected from young mares (≤14 years old) [22]. Moreover, a reduced quantity of mitochondrial DNA, but not genomic DNA, has been reported in Day 7 blastocysts from mares older than 16 years of age compared to younger counterparts, suggesting a reduction in mitochondria/cell ratio in old mares’ blastocysts [23]. Taken together, these alterations could lead to the alterations in embryo mobility and/or day of fixation that have been observed in some studies [9, 28].

Altogether, these data indicate that the quality of embryos is reduced in older mares. Considering the Developmental Origins of Health and Diseases (DOHaD) [35], dam age could thus also affect feto-placental and post-natal growth as well as long-term offspring health. To the authors’ knowledge, the effect of maternal age on gene expression in the equine embryo, however, has not yet been explored.

The aim of this study was to determine the effect of maternal age in the mare on embryo gene expression at the blastocyst stage. Young (<10 years) and older
(>10 years) multiparous mares were inseminated with the same stallion. Blastocysts were collected on Day 8 post-ovulation and bisected to separate the pure trophoblast (TE_part) from the inner cell mass enriched hemi-embryo (ICMandTE). Gene expression was analyzed by RNA-seq in each compartment. Deconvolution was used to isolate gene expression of inner cell mass from trophoblast cells in ICMandTE part.

## Results

### Embryo recovery rates, diameter, total RNA content and quality

Altogether, 36 uterine flushings were performed (12 in YM <10 years old and 24 in OM >10 years old) to obtain one embryo per mare. For 5 mares, after a negative embryo collection, a second attempt was performed. Characteristics of mares are detailed in Table 1. One YM and 3 OM had a double ovulation but no twin embryo was obtained. Post-breeding uterine fluid accumulation was checked at the time of ovulation confirmation and none of the mares presented sign of pathological fluid accumulation.

Embryo recovery rates (ERR) per mare and per ovulation were, respectively, 42 and 39% in YM and 25 and 22% in OM, with no significant difference between groups ($p$ = 0.45 for ERR/mare and $p$ = 0.46 for ERR/ovulation). All embryos were expanded blastocysts, grade I or II according to the embryo classification of McKinnon and Squires [36].

Embryo diameter, sex and RNA yield are presented in Additional Table 1. Embryo diameter ranged from 409 µm to 2643 µm, with no effect of group on embryo diameter ($p$ = 0.97). RNA yield per embryo ranged from 16 ng to 2915.5 ng and was not related to mares' age ($p$ = 0.92). There is a clear exponential relation between RNA yield/embryo and embryo size (Fig. 1a).

The median RNA Integrity Number (RIN) was 9.8 (8.7–10 range). Between 29.4 and 69.5 million reads

per sample were obtained after trimming. On average, 83.61% of reads were mapped on the modified EquCab 3.0 using STAR and 66.66% of reads were assigned to genes by featureCounts.

### Deconvolution of gene expression to discriminate ICM and TE gene expression on ICMandTE hemi-embryos

After keeping genes with more than 3 non-zero count values at least in one group (OM or YM) per hemi-embryo (ICMandTE or TE_part), 16,741 genes were conserved for deconvolution. In addition, twenty-seven genes were removed because their variance was null in the TE_part. However, for these genes, the mean counting in ICMandTE samples was above 10 counts, excepted for one which had a mean expression of 34 counts. On the other hand, this gene was only expressed in 5 ICMandTE over the 10 counts threshold.

During the analysis, 3 additional genes were excluded because the deconvolution quality for these genes was not sufficient. Therefore, at the end of deconvolution algorithm, 16,711 genes were available for differential analysis.

Before deconvolution, 341 genes were differentially expressed (FDR < 0.05) between the ICMandTE and the TE_part (Fig. 2a). After deconvolution, comparison between DeMixT_ICM1 and DeMixT_TEpart yielded 7569 differentially expressed genes while the comparison DeMixT_ICMpart vs TE_part yielded 6085 differentially expressed genes, with 5925 in common (77% of the total number of unique DEG in all comparisons). Moreover, except one of the initially 341 differentially expressed genes before deconvolution, were also identified as differentially expressed in both post-deconvolution analyses. On the PCA graphic of individuals, ICMandTE and TE_part were not really separated before deconvolution (Fig. 2b). DeMixT_TEpart and TE_part were partly superposed, suggesting that datasets before and after deconvolution have a

<table><caption>**Table 1** Mares' characteristics at embryo collection time</caption><thead><tr><th rowspan="2">Characteristics</th><th colspan="2">Young mares (YM) (n = 11)</th><th colspan="2">Old mares (OM) (n = 20)</th></tr><tr><th>Total</th><th>With embryo (n = 5)</th><th>Total</th><th>With embryo (n = 6)</th></tr></thead><tbody><tr><td>Breed</td><td>AA n = 9; SF n = 1; SB n = 1</td><td>AA</td><td>AA n = 13; SF n = 7</td><td>AA n = 4; SF n = 2</td></tr><tr><td>Age (in years)</td><td>6.00 ± 0.00</td><td>6.00 ± 0.00</td><td>11.95 ± 2.11</td><td>12.00 ± 2.90</td></tr><tr><td>Parity (number of foalings)</td><td>1.00 ± 0.00</td><td>1.00 ± 0.00</td><td>3.20 ± 1.15</td><td>3.17 ± 1.47</td></tr><tr><td>Weight (in kg)</td><td>536.62 ± 44.12</td><td>524.02 ± 61.39</td><td>589.35 ± 57.34</td><td>598.19 ± 76.83</td></tr><tr><td>BCS (scale 1–5)</td><td>2.16 ± 0.28</td><td>2.3 ± 0.33</td><td>2.61 ± 0.32</td><td>2.5 ± 0.59</td></tr><tr><td>Withers' height (in cm)</td><td>157.95 ± 3.74</td><td>157.40 ± 5.41</td><td>161.21 ± 4.23</td><td>162.00 ± 6.16</td></tr></tbody></table>

All mares were multiparous

AA Anglo Arab or Anglo-Arabian type, SF Selle Français, SB Saddlebred, BCS Body Condition Score. Age, parity, weight, BCS and height are presented as mean ± SD
<figure><img src="image_4.png" /></figure>

similar global gene expression; whereas the DeMixT_ICMCells group is clearly separated from both on Axis 2 (14.3% of variance), indicating that the deconvolution effectively enabled the separation of gene expression in the two cell types.

Six of the 12 genes previously identified as more expressed in ICM [37] were also more expressed in

ICMandTE vs TE_part comparison (Table 2). After deconvolution (comparison DeMixT_ICMCells vs TE_part), 3 more genes previously identified were also differentially expressed. However, for two of them (POU Class 5 Homeobox 1, *POU5F1* and Undifferentiated Embryonic Cell Transcription Factor 1, *UTF1*), their expression was increased in TE_part
after deconvolution. In TE, any genes previously identified were differentially expressed in the comparison ICMandTE vs TE_part, i.e., before deconvolution. After deconvolution, the expression of 2 genes was increased in TE_part compared to DeMixT_ICM CELLS. For further analysis, DeMixT_ICM ``cells and TE_part gene expression datasets are presented. Nevertheless, the comparison of maternal age on ICMandTE dataset (graphics, differential gene expression analysis and GSEA) was also performed and this analysis is available in Additional Figure 1 and Additional Tables 2 and 3.

Differential gene expression in deconvoluted ICM cells

After filtering out of genes with an average expression
<10 counts/maternal age group/hemi-embryo, 13,735
genes were considered as expressed in the ICM cells
from OM or YM embryos including 255 differentially
expressed genes (227 downregulated and 28 upregu-
lated in OM, with YM as the reference group) (Fig. 3a
and Additional Table 4). Only 218 and 12 genes out of
the down- and upregulated genes in OM compared to
YM, respectively, were associated to a known protein in
human. These downregulated genes in ICM of embryos
from OM compared to YM mainly belonged to cellular
process, metabolic process and biological regulation in
GO biological process gene sets (Fig. 3a). From down-
regulated genes in OM ICM compared to YM ICM,
84 pathways were statistically overrepresented in GO
Biological process leading to 16 most specific subclass
pathways according to hierarchy sorting by PANTHER
(Fig. 3b). Overrepresented pathways in downregulated
genes in OM ICM compared to YM ICM were linked
to cell division, embryo morphogenesis, histone meth-
ylation, protein modification processes and cytoskeleton
organization.

Upregulated differentially expressed genes (DEG) in
OM ICM compared to YM ICM were mainly part of bio-
logical regulation, response to stimulus, cellular process,
signaling and metabolic process in GO biological process
gene sets (Fig. 3a). No pathway was statistically overrep-
resented in upregulated genes in ICM of embryos from
old mares compared to ICM of embryos from young
mares using PANTHER overrepresentation testing.

Differential gene expression in the TE part

In the TE, 13,178 genes were considered as expressed in
OM or YM. Fourteen were differentially expressed (Addi-
tional Table 5) with 11 genes being downregulated and 3
being up regulated in OM, YM being the reference group
(Fig. 3a). Only 12 genes were associated to a known pro-
tein in human. Downregulated genes in OM compared to
YM mainly belonged to both cellular and metabolic pro-
cess in GO biological process gene sets.

Analysis of the variability of differentially expressed genes

As the embryo size was highly variable within each
group, a PCA analysis was made on DEG in ICM and
TE to verify that observed differences were effec-
tively due to maternal age more than to embryo size
(Fig. 1b). Respectively 3 YM and 3 OM were classified
in small (<1100 µm in diameter) while the remain-
ing 2 YM and 3 OM were classified in large group
(> 1100 µm in diameter).

In the ICM, the first axis of the PCA explained 58.3% of the variance of the data and clearly separated OM from YM, indicating than more than half of the variance in DEG is due to maternal aging. The second axis slightly separated genes according to embryo size but explained only 11.9% of the variance, indicating that in the ICM, embryo size explains a minor part of the variance. These observations were confirmed by the Partial redundancy analysis. Indeed, in ICM, 43% of the variance was explained by maternal age while only 14% was explained by embryo size. Moreover, the 20 variables that contrib- uted the most to the separation are mostly explaining dif- ferences related to maternal age.

In the TE part, both the first and the second axis
(explaining respectively, 67.3 and 16.5% of data vari-
ability) separated both maternal age and embryo size
groups. The partial redundancy analysis also showed
a less clear impact of maternal age in DEG from TE.
Indeed, mares' age explained 31% of the variance of
DEGs whereas 22% was explained by embryo size.
Moreover, although most of genes showed a clear con-
tribution in the construction of the first axis, some
such as cystathionine beta-synthase (CBS) and CD69
molecule (CD69) were clearly involved in the construc-
tion of the second axis.

(See figure on next page.)

**Fig. 2** Gene expression in ICM and TE before and after deconvolution using DeMixT. **A** Venn diagrams of the differential gene expression in ICMandTE vs TE part (before deconvolution), DeMixT_ICM ``cells vs DeMixT_TE ``cells (after deconvolution) and DeMixT_ICM ``cells vs TE ``part (gene expression of ICM after deconvolution vs gene expression in TE_part without deconvolution); **B** Principal Component Analysis of gene expression of DeMixT_ICM ``cells, DeMixT_TE ``cells, ICMandTE and TE part datasets. Deconvolution was used to isolate gene expression of ICM and TE cells in ICMandTE hemi-embryos. ICMandTE: inner cell mass + trophoblast; TE part: pure trophoblast. Here trophoblast represents trophectoderm + endoderm
<figure><img src="image_5.png" /><figcaption>**Fig. 2** (See legend on previous page.)</figcaption></figure>
<table><caption>**Table 2** Comparison of selected genes expression before and after deconvolution</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Gene name</th><th rowspan="2">Ensembl ID</th><th colspan="2">ICMandTE vs TE_part</th><th colspan="2">DeMixT_ICMpart</th></tr><tr><th>log2FC from DeSeq2</th><th>padj</th><th>log2FC from DeSeq2</th><th>Padj</th></tr></thead><tbody><tr><td rowspan="10">ICM</td><td>SOX2</td><td>ENSECAG00000010653</td><td>4.49</td><td>0.02</td><td>6.49</td><td>1.7E-06</td></tr><tr><td>NANOG</td><td>ENSECAG00000012614</td><td>4.24</td><td>1.0E-07</td><td>5.98</td><td>2.4E-15</td></tr><tr><td>SPP1</td><td>ENSECAG00000017191</td><td>4.02</td><td>8.4E-09</td><td>5.74</td><td>5.1E-16</td></tr><tr><td>LIN28B</td><td>ENSECAG00000020994</td><td>3.23</td><td>1.3E-10</td><td>5.00</td><td>4.7E-22</td></tr><tr><td>SMARCA2</td><td>ENSECAG00000024187</td><td>0.62</td><td>0.97</td><td>1.25</td><td>0.07</td></tr><tr><td>POU5F1 (OCT4)</td><td>ENSECAG00000008967</td><td>0.60</td><td>0.009</td><td>-2.04</td><td>4.6E-07</td></tr><tr><td>ID2</td><td>ENSECAG00000008738</td><td>0.44</td><td>0.29</td><td>1.24</td><td>9.0E-05</td></tr><tr><td>DNMT3B</td><td>ENSECAG00000012102</td><td>0.35</td><td>0.24</td><td>0.8</td><td>0.06</td></tr><tr><td>DPPA4</td><td>ENSECAG00000013271</td><td>0.35</td><td>0.04</td><td>1.14</td><td>1.8E-05</td></tr><tr><td>SALL4</td><td>ENSECAG00000018533</td><td>0.10</td><td>1</td><td>0.60</td><td>0.01</td></tr><tr><td>KLF4</td><td>ENSECAG00000010613</td><td>-0.04</td><td>1</td><td>-0.81</td><td>0.39</td></tr><tr><td>UTF1</td><td>ENSECAG00000039888</td><td>-0.13</td><td>1</td><td>-2.21</td><td>7.3E-05</td></tr><tr><td rowspan="7">TE</td><td>TFAP2A</td><td>ENSECAG00000017468</td><td>-0.17</td><td>0.46</td><td>-0.30</td><td>0.004</td></tr><tr><td>CDX2</td><td>ENSECAG00000027754</td><td>-0.14</td><td>0.77</td><td>-0.23</td><td>0.07</td></tr><tr><td>ELF3</td><td>ENSECAG00000014608</td><td>-0.13</td><td>1</td><td>-0.19</td><td>0.27</td></tr><tr><td>GATA2</td><td>ENSECAG00000016768</td><td>-0.13</td><td>1</td><td>-0.29</td><td>0.048</td></tr><tr><td>GATA3</td><td>ENSECAG00000024574</td><td>-0.10</td><td>0.87</td><td>-0.05</td><td>0.64</td></tr><tr><td>TEAD4</td><td>ENSECAG00000011303</td><td>-0.10</td><td>0.96</td><td>0.00</td><td>1</td></tr><tr><td>FREM2</td><td>ENSECAG00000020410</td><td>0.07</td><td>1</td><td>0.19</td><td>0.44</td></tr></tbody></table>

Gene expressions were obtained from RNA of 11 equine embryos bisected in two hemi-embryos: one part is composed only of trophoblast (TE), TE_part, while the other part is composed of TE and inner cell mass (ICM), ICMandTE. As it is impossible to estimate the proportion of each cells in ICMandTE, deconvolution algorithm (package DeMixT) was used to estimate gene expression of these different kind of cells. DeMixT_ICM ``cells'' dataset is the deconvoluted gene expression of ICM cells from ICMandTE. Log2 fold change (log2FC) and padj (adjusted *p*-value with Benjamini-Hochberg correction) were obtained with Deseq2 package. TE_part is the reference group in both analyses: when log2 fold changes (log2FC) > 0, gene is more expressed in ICMandTE or DeMixT_ICM ``cells'' while when log2FC < 0, gene is more expressed in TE_part. Green is used to represent gene differentially expressed in the present study. Orange is used to represent gens that have been previously identified as predominant in the ICM [37] but which are identified here as predominant in the TE

**Gene set enrichment analysis in deconvoluted ICM cells**

After Entrez Gene ID conversion, 12,344 genes were considered expressed in ICM. In the GO Biological Process database, 106 identified gene sets differed according

to maternal age in the ICM (11 enriched in OM and 95 enriched in YM) (Additional Table 6). In the KEGG database, 14 gene sets were perturbed (Additional Table 6).
<figure><img src="image_6.png" /><figcaption>**Fig. 3** Analysis of differentially expressed genes (DEG) in embryos according to maternal age. A representation of down (blue) and upregulated (red) DEG in ICM (from DeMixT_ICMDown (blue) and upregulated (red) DEG in ICM (from DeMixT_ICM ``cells data obtained after deconvolution of ICMandTE using DeMixT R package [38, 39]) and TE (from TE_part dataset) of embryos from OM. *Nota bene*: YM is the reference group. Functional classifications (for down regulated genes in ICM and TE and for up-regulated genes in ICM) using GO Biological Process obtained with PANTHER online software are presented as pie charts. For upregulated DEG in TE, no functional classification was made. B Bar chart presenting p-adjusted of most specific subclass pathways provided by PANTHER online software after statistical overrepresentation test (Fisher's Exact test and Benjamini-Hochberg correction) with the Human GO Biological Process annotation on down-regulated DEGs in ICM of OM. *Nota bene*: YM is the reference group. DEG: Differentially Expressed Genes (FDR < 0.05); TE: Trophoblast; ICM: Inner Cell Mass; OM: Old mares; YM: young mares</figcaption></figure>
<figure><img src="image_7.png" /></figure>

Using SUMER analysis (Fig. 4), gene sets enriched in OM ICM compared to YM ICM were mainly grouped under terms “Respiratory electron transport chain” and “Translational elongation”. The first term represents

mitochondrial function while the second represents translation function. Moreover, terms relative to “Systemic lupus erythematosus” are conserved by SUMER analysis and enriched in OM ICM compared to YM
ICM. Genes involved in these enrichments were also present in dysregulated metabolic pathways.

Gene sets enriched in YM ICM compared to OM ICM were involved in different biological processes such as “Mitotic nuclear division” which is relative to mitosis. Under this term, “Spindle organization” and “Negative regulation of supramolecular fiber organization” terms represent the assembly of spindle during mitosis. “Chromosome segregation” is also found in this grouping. Furthermore, terms under “Golgi vesicle transport” seem to be involved in cytokinesis and vesicle trafficking. The catabolism of protein and mRNA is represented by “Ubiquitin mediated proteolysis”. Translation function is represented by several terms grouped under “covalent chromatin modification” in the mitosis part.

“Leukocyte transendothelial migration” term is also enriched in YM ICM compared to OM ICM and seems to represent cell signaling and adhesion. Indeed, this term regroups several pathways such as “Focal adhesion”, “Adherens junction” and “WNT signaling pathway”.

Finally, further terms related to embryo development could be found in different groups such as “Dendrite development”, “Cell differentiation in spinal cord” or “endothelial cell development”.

## Gene set enrichment analysis in the TE part

After Entrez Gene ID conversion, 11,874 genes were considered expressed in TE from OM or YM embryos. In GO BP database, 38 gene sets were perturbed by mares' age in the TE of embryos (11 enriched in YM and 27 enriched in OM) (Additional Table 7). In the KEGG database, 5 pathways were enriched in the TE of OM compared to TE from YM embryos (Additional Table 7).

As in ICM part, using SUMER analysis, translation (represented by both “Translational termination” and “Ribosome” terms) and mitochondrial function (represented by “Oxidative phosphorylation” associated terms) were enriched in OM embryos compared to YM embryos. “Proteasome” term is enriched in OM embryos, compared to YM representing protein catabolism.

Moreover, transcription (represented by “Peptidyl lysine trimethylation” term) and embryo development (represented by “spinal cord development”) were, as in ICM, also, enriched in YM embryos compared to OM. Finally, in TE, “Regulation of action potential” term and associated terms were enriched in YM compared to OM. These terms could be related to ion movement.

## Discussion

Maternal aging did not affect embryo recovery rate (ERR), embryo size nor total RNA quantity. It is important to notice that, as the experimental protocol required

only 5 embryos per group, a limited number of embryo collection attempts were performed, resulting in a limited statistical power for these analyses. ERR, especially in young mares, was not as good as in commercial embryo transfer programs [31, 42]. Moreover, total RNA content in equine embryo was exponentially correlated to embryo diameter and was not affected by maternal age.

As ovulation check was only performed 48 hours after induction of ovulation, there is a bias in the estimation of the embryo age in the present study. This certainly partly explains the large variation in embryo size observed here. Different studies, however, have shown that embryo size is highly variable [43–45], even when ovulation determination is accurate to the nearest 6 hours [46]. Accurately determining fertilization time in vivo is not possible and the size of 142 equine embryos recovered at 8 days post-ovulation was reported to range from 120 to 3980 µm with an average of 1132 µm [30]. In the present project, all embryos remain within these limits.

Although embryo size could be considered as normal for Day 8 embryos, embryo sizes were highly variable and not equally distributed between groups. Embryo size was considered as a cofactor in the differential analysis. In addition, partial redundancy analysis showed that embryo size explained less than 25% of the variability of DEGs in both ICM and TE. As effects of maternal age on equine embryo growth are not so clear [24–28], an interaction between both effects on embryo gene expression can not be excluded. Maternal age, however, explained 43 and 31% of the total variability of DEGs in both ICM and TE, respectively, suggesting that identified DEGs were mainly due to maternal age effects. It is important to note that in ICM, the proportion of variability explained by maternal age was 3 times greater than that of embryo size, whereas in the TE, this difference was smaller. This could be explained by the fewer number of DEGs identified in TE compared to ICM but also because trophoblast expansion was possibly more correlated to embryo size than the ICM.

As in tumors, today, the only way to analyze separate expression of ICM and TE cells in mammalian embryos is to use micromanipulation but micromanipulation is expensive in time and requires adapted skills and material. In this study, micromanipulation to properly separate ICM and TE was not possible. Therefore, deconvolution appeared to be a good option to analyze the effects of maternal ICM gene expression in mixed ICM and TE samples. Deconvolution seemed to have been effective on embryo to separate ICM and TE gene expression. Indeed, the 341 genes differentially expressed in ICMandTE vs TE_part were also identified by comparing TE with post-deconvolution ICM cells gene expression dataset. Moreover, after deconvolution, more genes
were differentially expressed between TE and ICM data-
sets and among them, more genes identified in litera-
ture were differentially expressed. Finally, PCA showed
that deconvolution separate DeMixT_ICM CELLS from
other datasets while before deconvolution ICMandTE
and TE_part were partly superposed. Only POU5F1 and
UTF1 differed in terms of allocation between previously
published data [37] and the present study. In data before
deconvolution, expression of POU5F1 and UTF1 were
extremely variable according to embryo diameter in both
compartments. Deconvolution reduced the variability
of both gene expression in ICM but not in TE. Studies
showed that POU5F1 is expressed in both compartments
until Day 10 post ovulation in equine in vivo produced
embryos [47] but its expression decreases progressively
from Day 7 in TE cells. Thus, the statistically increased
expression of POU5F1 in TE on Day 8 is compatible with
previously published data.

Therefore, it was chosen to only present the analysis
of the effects of age in DeMixT_ICM CELLS dataset. In all
cases, deconvolution did not change the principal results
of maternal aging on embryo gene expression because the
comparison of maternal age in ICMandTE and DeMixT_
ICM CELLS gave comparable results.

Genes with higher expression in both the ICM and the TE of OM embryos were overrepresented for functional terms related to mitochondrial function and translation with ribosome biogenesis. Those with higher expression in YM embryos were involved in transcription and chromatin modification (Fig. 5). Mitosis, signaling and adhesion pathways were additionally altered and the expression of genes involved in embryo development was particularly reduced in the ICM hemi-embryos from old mares. Finally, in the TE, ion movement related genes were affected.

Several studies on mature oocytes in women and mice, using microarrays and RNAseq, had shown a differential expression of genes involved in mitochondrial, cell cycle, cytoskeleton functions and related pathways according to maternal age [50–52], as observed in the present study in equine blastocysts, especially in the ICM. In equine follicles and oocytes, the expression of genes essential for oocyte maturation, as studied by RT-qPCR, quantitatively and temporally differed with maternal age, suggesting a desynchronization [16]. Moreover, a study on the expression of 48 genes that have been previously identified as related to aneuploidy and maternal aging in humans showed differential expression in cumulus-oocytes complexes from young (< 12 years of age) and old (> 18 years old) mares [17]. Genes related to mitochondrial replication and function were also down-regulated in oocytes of old mares [23]. Even effects of maternal aging on gene expression in equine oocytes, using high

throughput techniques, are not yet available, results on
equine Day 8 blastocysts suggest that mitochondrial per-
turbations present in the oocyte persist until the blasto-
cyst stage. In humans, single-embryo RNA-seq analysis
showed that maternal aging alters the transcriptome of
blastocysts derived from ICSI. Genes altered by maternal
age were related to segregation of chromosomes, embry-
onic growth and cell signaling [53, 54], as observed in our
study on equine blastocysts.

Unlike the human embryo, the equine embryo remains spherical, grows very fast and is surrounded by its glycoprotein capsule up to 22 days post ovulation [43]. At 8 days post ovulation, equine embryos are therefore growing and the endodermal layer entirely delineates the future yolk sac (Fig. 5) [48]. Moreover, whereas in Day 7 embryos, the epiblast is lined by a continuous layer of trophoblast named the polar trophoblast. As the blastocyst grows, epiblast cells phagocyte polar trophoblast cells that become discontinuous. During this time, epiblast and polar trophoblast cells divided and are bound by tight junctions [49]. The turnover in focal adhesion is essential for cell movement and phagocytosis [55].

Cell signaling and adhesion function were altered in OM embryos compared to YM. Indeed, KEGG pathways “focal adhesion” and “adherens junction” regrouped in “Leukocyte transendothelial migration” on Fig. 3 are disturbed in old mares’ embryos. Focal adhesions are subcellular structures observed in adherent cells, that are clusters of different structural and signaling molecules and integrins that form a physical link between the actin cytoskeleton and the extracellular matrix [56, 57]. Integrins are the major receptors for cell adherence to extracellular matrix proteins in focal adhesions by playing an important role in transmembrane connections to the cytoskeleton and activation of several signaling pathways, as growth factor receptors [58]. Integrin Subunit Alpha V (ITGAV) expression was reduced in the ICM of embryos collected from OM compared to the ones from YM. In different human breast cancer cell lines, cell proliferation, invasion and renewal are restricted when ITGAV is inhibited [59].

Moreover, the activity and localization of small GTPases within the Rho family must be precisely regulated as the communication between them is indispensable to achieve the turnover of actin filaments and cell adhesions necessary for cell migration, that are essential processes for the proper development of the embryo [60]. The Rho GTPase Activating Protein 9 and 35 (ARHGAP9 and ARHGAP35) were downregulated in OM ICM compared to YM ICM. ARHGAP35 codes for P190 Rho GTPase activating protein that mediates signalization through integrin-dependent adhesion in cells.
<figure><img src="image_8.png" /><figcaption>**Fig. 5** Schematic representation of the observed effects of maternal aging on ICM and TE in equine blastocysts. Normal blastocyst development from day 7 to day 9 is represented inside black square using literature [48, 49]. Light blue boxes represent biological processes that are enriched in YM embryos and light red boxes represent biological processes that are enriched in OM embryos</figcaption></figure>
Aberrant neural morphogenesis, related to an accumulation of polymerized actin in embryoneural tube cells, is reported in mice lacking P190 Rho GTPase activating proteins [61]. Furthermore, P21 activated kinases (PAK) are effector kinases for other small Rho GTPases such as Rac that regulate the polymerization of actin and can affect the microtubule organization [62]. Once activated, PAKs promote the turnover of focal adhesions by regulating transcription [63]. PAK3 was downregulated in the OM compared to YM ICM. However, in TE, SLIT-ROBO Rho GTPase Activating Protein 3, *SRGAP3* was one of the upregulated genes in OM compared to YM. This gene regulates cytoskeletal reorganization and plays an important role in the formation of cell-cell adhesion [64]. Therefore, maternal aging seems to alter focal adhesion in both ICM and TE of embryos but, while in ICM, it seems to reduce cohesion and signalization, in TE it seems to reinforce them.

Several signaling pathways required for development were altered and numerous genes related to focal adhesions and signaling pathways were down-regulated in OM ICM compared to YM (*ARHGAP9* and *ARH-GAP35*; B-Raf Proto-Oncogene, Serine/Threonine Kinase, *BRAF*; Filamin C, *FLNC*; FYN Proto-Oncogene, Src Family Tyrosine Kinase, *FYN*; *PAK3*; Protein Phosphatase 1 Regulatory Subunit 12A, *PPP1R12A*; Talin 2, *TLN2*; X-Linked Inhibitor Of Apoptosis, *XIAP*). Among all signaling pathways, transforming growth factor beta (TGF-β), Wg et Int (Wnt), mitogen-activated protein kinases (MAPK) and apoptosis signaling were altered by maternal age in equine embryos. Indeed, transforming growth factor beta receptor 1 (*TGFBR1*) was downregulated in the ICM of embryos from old mares compared to the ones from young mares. *TGFBR1* is one of the two dimers for the receptor for transforming growth factors beta that is essential in the generation of axes and cell fate commitment during mammal embryogenesis [65]. Moreover, Casein Kinase 1 Alpha 1, (*CSNK1A1*) and *Protein Phosphatase 2 Regulatory Subunit B*'Epsilon (*PPP2R5E*), involved in Wnt signaling pathway, were downregulated in OM embryos compared to YM. Wnt signaling pathway is also involved in nervous development system and may act on cell determination in mice embryos [66]. Moreover, Doublesex And Mab-3 Related Transcription Factor 3 (*DMRT3*) was downregulated in the OM ICM compared to YM. *DMRT3* is required for the proper development of the brain in prenatal chicken and mice [67].

Cell cycle related pathways and genes were also disturbed by maternal aging. Indeed, the expression of Pleckstrin Homology domain-Interacting Protein (PHIP), also known as replication initiation determinant, was reduced in the OM ICM compared to YM. This gene is required for the initiation of DNA replication. Fibroblasts

originated from mice embryos with depleted *PHIP* expression presented less replication-initiation as well as abnormal replication fork progression events [68]. Moreover, several cyclins (CCN) were less expressed in the ICM of embryos from old mares compared to the ones from young mares (*CCNB3*, *CCNI* and *CCNT2*). Cyclins are key cell cycle regulators and, for instance cyclin B3 degradation allows the passage from metaphase (chromatids aligned on metaphase plate) to anaphase (segregation of sister chromatids) [69]. During metaphase and anaphase, spindle assembly and chromosome segregation could also be compromised. Indeed, genes involved in the formation of microtubules (Tubulin Gamma 2, *TUBG2*) [70], the generation of the mitotic spindle (Never in Mitosis A related kinase 9, *NEK9*) [71], the production of cohesion and condensing complexes (Structural Maintenance Of Chromosomes 4, *SMC4*) [72], the control of centriole length, important for centrosomes and microtubules organization (Centriolar Coiled-Coil Protein 110, *CCP110*) [73], the amplification of intra-spindle microtubules allowing the formation of a sufficient quantity of kinetochore (HAUS Augmin Like Complex Subunit 6, *HAUS6*) [74], the stabilization and protection of the cohesin complex association with chromatin (Cell Division Cycle Associated 5, *CDC45*; Shugosin 2, *SGO2*) [75, 76] and the production of microtubule-based motor proteins (Kinesin Family Member 14 and 23, *KIF14* and *KIF23*) [77] were downregulated in OM ICM compared to YM. Furthermore, Centromere protein F (CENPF), an important regulator of chromosome segregation during mitosis, was reduced in OM ICM compared to YM. CENPF is indispensable for murine embryo development from the first stages of development until at least the blastocyst stage [78]. Impaired CENPF expression could induce chromosome mis-segregation in the ICM of old mares' embryos compared to the ones from young mares.

Cytokinesis seems to also be altered by maternal age. CD2-associated protein (CD2AP) and Vacuolar Protein Sorting 4 Homolog A (VPS4A) were downregulated in OM ICM compared to YM. Both are involved in the late phase of cytokinesis which correspond to the abscission of the cell, producing two daughter cells [79, 80]. Globally, results suggest that cell division is impaired or indicate that there is less division in the ICM of old mares' embryos compared to the ones from young mares.

Transcription is regulated by epigenetic marks. Several pathways linked to methylation and acetylation were enriched in YM ICM and TE in relation to numerous upregulated genes in YM ICM compared to OM (i.e., Histone Deacetylase 1 and 2, *HDAC1* and *HDAC2*; Helicase, Lymphoid Specific, *HELLS*; DNA Methyltransferase 3 Alpha, *DNMT3A*; Lysine Acetyltransferase 6B, *KAT6B*; Lysine Demethylase 4A, *KDM4A*; Lysine
Methyltransferase 2A, *KMT2A*; Nuclear Receptor Binding SET Domain Protein 3, *NSD3*; SET Domain Containing 5, *SETD5*; Tet Methylcytosine Dioxygenase 1, *TET1*). At blastocyst stage, de novo methylation is more important in the inner cell mass than in trophoblast cells and is essential for the proper development of mammalian embryos [81]. *DNMT3A* is essential for de novo methylation and required for embryo viability [82]. As it has been shown that maternal aging increases the acetylation of lysine 12 of histone H4 in mice oocytes, affecting fertilization and embryo development [83], it is likely that the methylation and acetylation dysregulation observed in the embryos of old mares could be inherited from the oocyte rather than directly due to embryo environment.

Both in OM ICM and TE compared to YM, several pathways related to mitochondrial function were enriched. Mitochondria are long known to be the “powerhouse” of cells, producing most of the necessary ATP for cellular function [84]. ATP production is possible through oxidative phosphorylation enabled by the electron transport chain in the inner membrane of mitochondria. This process begins with the oxidation of NADH or FADH that generates a proton gradient across the mitochondrial membrane, driving ATP synthesis through the phosphorylation of ADP via ATPase [85]. Enriched pathways in OM compared to YM embryos demonstrate that the NADH dehydrogenase (complex I), the oxidoreductase acting on a heme group of donors (complex III) with the global electron carrier, and the proton-transporting ATP synthase (1st part of the complex V) activities were particularly enriched while the second part of complex V is depleted. Since all complexes are required to produce ATP, the impaired complex generation may lead to mitochondrial stress and uncoupled respiratory chain, resulting in poor ATP production. In human IVF embryos, reduced mitochondrial respiration efficiency, with a negative effect on embryo development, is observed when oocytes are collected from old women [86]. In mammalian expanded blastocysts, energy requirement is high, probably because of the need for ion transport enabling quick blastocoel cavity expansion. Defects in metabolic requirements during blastocyst development are associated with impaired implantation [87]. Equine embryos are not an exception [88].

Mitochondria being maternally inherited, uncoupling of the mitochondrial respiratory chain may originate from mitochondrial defects already present in oocytes [89]. Studies have shown that, in human, poor oocyte quality is responsible for the decline in fertility associated with maternal aging. The reduced oocyte developmental competence in older women is related to mitochondrial dysfunction [90]. In horses' oocytes, before in vitro maturation, no difference between young and old mares is

observed while after, reduced mitochondrial DNA copy numbers, increased swelling and reduced mitochondrial cristae are reported in oocytes collected from mares >12 years. This suggests that oocyte mitochondria from old mares present deficiencies that become problematic when energy needs are high, as during in vitro maturation or subsequent embryo development [21]. Because mitochondrial replication in equine embryos does not begin until the blastocyst stage [91], oocyte mitochondrial numerical or functional deficiencies may compromise early embryogenesis.

In TE, the GO biological process “Regulation of action potential” and “Transmission nerve impulse”, related to ion movement, were enriched in YM compared to OM embryos. These processes are associated with genes involved in the action potential dependent ion channels. These ion channels such as NA⁺, K⁺, ATPase are essential for the accumulation of blastocoel and subsequently yolk sac fluid in equine embryos [92]. Their inhibition is detrimental for embryo development [93]. Thus, both reduced ATP production and disturbed ion channels' function may impact embryo growth, although effects of mares' age on embryo growth are controversial [24–28].

## Conclusions

This is the first study showing an effect of maternal age on gene expression in the equine blastocyst at Day 8 post ovulation. Maternal age, even for mares as young as 10 years old, disturbs several processes related to chromatin segregation, cytokinesis, apposition of epigenetic marks and mitochondrial function and that the ICM appears to be more affected than the TE. These perturbations on Day 8 post ovulation may affect the further development of the embryo and as a result, may contribute to decreased fertility due to aging. So far, maternal age does not appear to have an impact on foal morphology at birth but, as for placental function, more subtle effects may exist. Here, as embryos were collected long before implantation, it is impossible to know if these embryos would have implanted and if they would have produced a live foal or not. Studies on long-term growth and health of foals born to the same mares are currently being performed.

## Methods

### Ethics

The experiment was performed at the experimental farm of IFCE (research agreement C1903602 valid until March 22, 2023). The protocol was approved by the local animal care and use committee (“Comité des Utilisateurs de la Station Expérimentale de Chamberet”) and by the regional ethical committee (“Comité Régional d’Ethique pour l’Expérimentation Animale du Limousin”, approved
under N° C2EA - 33 in the National Registry of French Ethical Committees for animal experimentation) under protocol number APAFIS#14963-2018050316037888 v2. All experiments were performed in accordance with the European Union Directive 2010/63EU. The authors complied with the ARRIVE guidelines.

## Embryo collection

Thirty-one non-nursing multiparous Saddlebred mares (French Anglo-Arabian and Selle Français breeds) were included in this study. Multiparous mares were defined as dams that had already foaled at least once. Mares were bred and raised in the "Institut Français du Cheval et de l'Equitation" experimental farm (Chamberet, France, 45° 34'55.17"N, 1°43'16.29"E, 442 m). During the experimental protocol, mares were managed in herds in natural pastures 24 h/day with free access to water. The experiments took place from April 1st to May 23rd, 2019. All mares remained healthy during this period.

Mares were allocated to 2 groups according to their age: young mares (YM, < 10 years, n = 11) and older mares (OM, ≥ 10 years, n = 20). As several studies had shown an increase embryo loss and/or reduced pregnancy rates as early as 10 years old, the age threshold was defined in accordance with the literature [3, 4]. Before experimentation, mare's withers' height was measured. During the experimentation, mares were weighted and mares' body condition score (BCS, scale 1 to 5) was determined by one trained operator [94].

The mares' estrous period was monitored routinely by ultrasound with a 5 MHz trans-rectal transducer. During estrus, ovulation was induced with a single injection of human chorionic gonadotropin (i.v.; 750 - 1500 IU; Chorulon® 5000; MSD Santé animale, France) as soon as one ovarian follicle > 35 mm in diameter was observed, together with marked uterine edema. Ovulation usually takes place within 48 h, with > 80% occurring 25 to 48 h after injection [95, 96]. At the same time, mares were inseminated once with fresh or fresh overnight cooled semen containing at least 1 billion motile spermatozoa from a single fertile stallion. Ovulation was confirmed within the next 48 hours by ultrasonography.

Embryos were collected by non-surgical uterine lavage using prewarmed (37°C) lactated Ringer's solution (B.Braun, France) and EZ-Way Filter (IMV Technologies, France) 10 days after insemination, i.e. approximately 8 days post ovulation. All embryo collections were performed by the same operator, properly trained but not performing embryo collection routinely. Just after embryo collection, mares were treated with luprotiol an analogue of prostaglandin F2α (i.m; 7.5 mg; Prosolvin, Virbac, France).

The aim of the embryo collection was to obtain 5 embryos/group with each embryo coming from a different mare. Therefore, some mares that failed to produce an embryo at their first attempt had been bred again.

## Embryo bisection and RNA extraction

Collected embryos were immediately observed using a binocular magnifier and photographed with a size standard. Embryo diameter was subsequently determined on photographs using ImageJ® software (version 1.52a; National Institutes of Health, Bethesda, MD, USA). Embryos were washed 4 times in commercially available Embryo holding medium (IMV Technologies, France) at 34°C and then bisected with a microscalpel under binocular magnifying glass to separate trophoblast part, at this stage composed of trophectoderm and endoderm (TE_part) from inner cell mass (composed of epiblast cells) enriched hemi-embryos (ICMandTE) (Fig. 6). Each hemi-embryo was transferred to an Eppendorf DNA LoBind tube (Eppendorf, Montesson, France) with extraction buffer (PicoPure RNA isolation kit, Applied Biosystems, France) and incubated for 30 min at 42°C prior to storage at -80°C.

RNA was extracted from each hemi-embryo following the manufacturer's instructions using PicoPure RNA isolation kit (PicoPure RNA isolation kit, Applied Biosystems, France), which included a DNAse treatment. RNA quality and quantity were assessed with the 2100 Bioanalyzer system using RNA 6000 Pico or nano kit (Agilent Technologies, France) according to the manufacturer's instructions (Additional Figure 2).

## RNA sequencing

Five nanograms of total RNA were mixed (Additional Figure 2) with ERCC spike-in mix (Thermofisher Scientific, France) according to manufacturer recommendations. Messenger RNAs were reverse transcribed and amplified using the SMART-Seq V4 ultra low input RNA kit (Clontech, France) according to the manufacturer recommendations. Nine PCR cycles were performed for each hemi-embryo. cDNA quality was assessed on an Agilent Bioanalyzer 2100, using an Agilent High Sensitivity DNA Kit (Agilent Technologies, France). Libraries were prepared from 0.15 ng cDNA using the Nextera XT Illumina library preparation kit (Illumina, France). Libraries were pooled in equimolar proportions and sequenced (Paired end 50–34 pb) on NextSeq500 instrument, using a Next-Seq 500 High Output 75 cycles kit (Illumina, France) (Additional Figure 2). Demultiplexing was performed with bcl2fastq2 version 2.2.18.12 (Illumina, France) and
<figure><img src="image_9.png" /></figure>

adapters were trimmed with Cutadapt version 1.15 [97]. Only reads longer than 10pb were kept.

## Embryo sexing

Embryo sex was determined for each embryo in order to take in consideration embryo sex in the statistical analysis as sexual dimorphism was observed in blastocysts in other species [98].

The expression of *X Inactive Specific Transcript* (*XIST*) was analyzed. *XIST* is a long noncoding RNA from chromosome X responsible for X inactivation when two X chromosomes are present, meaning that only female embryos express *XIST*. *XIST* was shown to be expressed in equine female blastocysts on Day 8 post ovulation [37, 44]. *XIST* is not annotated in the available equine Ensembl 3.0.99 or NCBI annotation. Thus, a de novo transcript assembly was made using StringTie version 2.1.2 [99] following manufacturer instructions to identify *XIST* in equine embryo samples (Additional Figure 2). The genome annotation file output was aligned on Ensembl 99 EquCab3.0 assembly using IGV version 2.7.2 [100]. As *XIST* position and length on X chromosome has been previously determined [37] and as *XIST* is annotated in humans, the most similar contig was selected. The obtained *XIST* sequence was 22,600 bp long and composed of 8 exons and 7 introns. *XIST* sequence was added to

FASTA files from Ensembl 99 EquCab3.0 assembly and to the genome annotation file from Ensembl 99 EquCab3.0 assembly as a unique fictitious chromosome with one gene composed of one transcript.

*XIST* expression differed between individuals but was relatively close in quantification between related hemi-embryos. Therefore, seven embryos were determined as females (4 in the YM group and 3 in the OM group) while 4 were considered as males (1 in the YM group, and 3 in the OM group).

## RNA mapping and counting

Alignment was performed using STAR version 2.6 [101] (2-pass mapping with the genome annotation file using the following parameters: --alignSoftClipAtReferenceEnds No --alignEndsType EndToEnd --alignIntronMax 86,545 --alignMatesGapMax 86,545) on previously modified Ensembl 99 EquCab3.0 assembly and annotation. Genes were then counted with featureCounts [102] from Subreads package version 1.6.1 (parameters: -p -B -C -M) (Additional Figure 2).

## Data analysis

All statistical analyses were performed by comparing OM to YM (set as reference group) using R version 4.0.2 [103] on Rstudio software version 1.3.1056 [104].
## Embryo recovery rate, embryo diameter and total RNA content analysis

Embryo recovery rate (ERR) per mare was calculated as number of attempts with at least one embryo collected/total number of attempts. ERR per ovulation was determined as number of collected embryos/total number of observed ovulations. They were analyzed using the Exact Fisher test to determine if age influences ERR.

For total RNA content analysis, as embryos were bisected without strict equality for each hemi-embryo, a separate analysis of ICMandTE and TE_part RNA quantities would not be meaningful. Thus, ICMandTE and TE_part RNA quantities were initially summed up and analyzed using a mixed linear model from nlme package version 3.1-148 [105], followed by 1000 permutations using PermTest function from pgirmess package version 1.6.9 [106]. Maternal age and embryo sex were considered as random effects.

Embryo diameter was analyzed with a fixed effects linear model of nlme package version 3.1-148 [105] including maternal age and embryo sex, followed by 1000 permutations using PermTest function from pgirmess package version 1.6.9 [106]. Variables were kept in models when statistically significant differences were observed. Differences were considered as significant for *p* < 0.05.

Moreover, the relation between embryo size and total RNA quantity/embryo was studied and represented using the nonlinear regression analysis (fitness to exponential growth equation) in GraphPad Prism software 8.0.1 for Windows (Graphpad Software, San Diego, California USA, www.graphpad.com).

## Deconvolution of gene expression in ICMandTE using DeMixT

A filter on counts obtained with featureCounts was applied: all genes with more than 3 non null count values in at least one group (OM or YM) per hemi-embryo (ICMandTE or TE_part) were conserved. The ICMandTE hemi-embryo was composed of both trophoblast and inner cell mass in unknown relative proportions. In order to estimate the relative gene expression of both cell types, the DeMixT R package version 1.4.0 was used [38, 39]. Starting from a dataset obtained from one cell type, DeMixT estimates the proportion of cells in the mixed samples and performs a deconvolution algorithm on gene expression. TE_part (reference cells) and ICMandTE (mixed cells type) datasets were thus analyzed using the following options: niter = 10; nbin = 50; if.filter = TRUE; ngene. selected.for.pi = 1500; ngene. Profileselected = 1500; nspikein = 0. Since the deconvolution algorithm

cannot be performed in the presence of null values, the value "1" was added to all count values in ICMandTE and TE_part data tables. Moreover, variance in reference cell dataset (TE_part) must not be null and gene with null variance had been removed prior to the deconvolution. Output datasets were DeMixT_ICM1 cells and DeMixT_TE1 cells, corresponding to the deconvoluted gene expression in ICM cells and TE cells of ICMandTE, respectively.

## Quality checks for deconvolution

At the end of deconvolution, a quality check is automatically performed by the analysis: for every gene, if the difference of average of deconvoluted expression of reference cells in mixed samples (here, gene expression of TE cells in ICMandTE dataset, i.e., DeMixT_TE_part cells) and average of expression of reference cells (here, gene expression of TE cells in TE_part) is >4, then, the deconvolution for concerned genes is not reliable. Therefore, genes in this case are filtered out.

Moreover, outputs of DeMixT_ICM ``cells vs DeMixT_TE ``cells, DeMixT_ICM ``cells vs TE ``part and ICMandTE vs TE ``part were compared with Deseq2 version 1.28.1 [107] to confirm that the deconvolution was effective at separating gene expression. Moreover, several genes proper to ICM and TE cells in equine embryos were identified using literature search [37]. Their raw and deconvoluted expressions were analyzed to check the reliability of deconvolution. Results of these analyses were represented through manually drawn Venn diagrams as well as principal component analysis (PCA) graphics of individuals, using FactoMineR version 2.4 [108], ggplot2 version 3.3.3 [109] and factoextra version 1.0.7 [110].

## Maternal age comparison for gene expression

For the comparison according to maternal age, ICMandTE, DeMixT_ICM ``cells and TE ``part datasets were used and a filter was applied on datasets: all genes with an average expression equal or above 10 counts in at least one group (OM or YM) per hemi-embryo (ICM or TE) were conserved. Differential analyses were performed with Deseq2 version 1.28.1 [107] with YM group as reference, without independent filtering. Genes were considered differentially expressed (DEG) for FDR < 0.05 after Benjamini-Hochberg correction (also known as false discovery rate, FDR). As ovulation was checked only every 48h and because embryos growth is exponential in the uterus, embryo diameter was considered as a cofactor in the model as well as embryo sex. Moreover, PCAs using factominer version 2.4 [108] on R software followed with a Partial
redundancy analysis using rda function from vegan package version 2.5-7 [111] were performed. For these analyses, embryo size was considered as a factor of two modalities: embryos < 1100 µm were considered as small and > 1100 µm were considered as large. Results were represented with a biplot graphic using the package factoextra version 1.0.7 [110]. The combined analyses were used to represent and determine the specific contribution of each of the two factors in the total variability observed in DEGs. The objective was to verify that the age of the mare was the main factor explaining the difference between the DEGs and not the differences in embryo size. On biplot graphics, only the 20 variables, i.e., genes, that contributed the most in axis construction, were represented.

Equine Ensembl IDs were converted into Human Ensembl IDs and Entrez Gene names using gorth function in gprofiler2 package version 0.1.9 [112]. Genes without Entrez Gene names using gprofiler2 were manually converted when Entrez Gene names were available, using Ensembl web search function [113].

Moreover, DEG were analyzed for classification, over- and under- representation (FDR < 0.05) using Fisher's exact test corrected with FDR and using the database Gene Ontology (GO) Biological Process complete on PANTHER version 16.0 web-based software [114]. Only the most specific subclass pathways, provided by PANTHER's automated hierarchy sorting, were conserved and represented.

## Gene set enrichment analyses (GSEA)

The count values were log transformed using RLOG function of DESeq2 version 1.28.1. Gene set enrichment analyses (GSEA) were performed on expressed genes using GSEA software version 4.0.3 (Broad Institute, Inc., Massachusetts Institute of Technology, and Regents of the University of California) [115, 116] to identify biological gene sets disturbed by age. GSEA was performed using the following parameters: 1000 gene set permutations, weighted enrichment analysis, gene set size between 15 and 500, signal-to-noise metrics. Molecular Signatures Database [117] version 7.1 (C2: KEGG: Kyoto Encyclopedia of Genes and Genomes, C5: BP: GO biological process) were used to identify most perturbed pathways. Pathways were considered significantly enriched for FDR < 0.05. YM group was considered as the reference group, which indicates that when the normalized enrichment score (NES) was positive, the gene set was enriched in the OM group while when NES was negative, the gene set was enriched in the YM group.

Enriched terms from GO BP and KEGG database were represented using SUMER analysis from SUMER R

package version 1.1.5 and using FDR q-values [40]. First, SUMER analysis reduces redundancy using weighted set cover algorithm and then clusters the results using affinity propagation algorithm. The weighted set cover algorithm performs a selection of the fewest number of gene sets that include all genes associated with the enriched sets, giving priority for the most significant terms. Then, the affinity propagation algorithm groups similar gene sets and identifies one representative gene set. This algorithm will repeat this clustering to obtain the most representative gene sets. Results were represented with graph modified using Cytoscape version 3.8.2 [41]. Gene sets were represented by nodes and the gene set size was represented by the size of the node. Node shape represented the gene set database (GO BP or KEGG). Blue nodes represented gene sets enriched in YM (NES < 0) while red nodes represented gene sets enriched in OM (NES > 0). Edge width represents the level of connection between representative gene sets (thinner edges represent the first clustering while thicker edges represent the second clustering of the affinity propagation algorithm).

## Data availability statement

The RNA sequencing data supporting the conclusions of this article are available in the repository [accession: GSE162893; https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE162893] from the GEO SuperSeries [accession: GSE193676; https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE193676].

## Supplementary Information

The online version contains supplementary material available at https://doi.org/10.1186/s12864-022-08593-7.

**Additional file 1: Supplementary Figure 1.** Analysis of ICM enriched and TE part before gene expression deconvolution using DeMixT. On first page, the analysis of differential expressed genes is represented. The second page represent the results of the gene set enrichment analysis of ICMandTE gene expression. The third page presents the clustering of gene sets altered by maternal age from SUMER analysis of the GO BP and KEGG terms of mixed ICMandTE and TE_part.

**Additional file 2: Supplementary Figure 2.** Pipeline of biostatistical analysis. Once extracted, total RNA qualification and quantification were obtained using a Bioanalyzer. For sequencing, 5 ng/sample of total RNA were used for paired end RNA sequencing (I2BC platform) using Illumina NextSeq 500 high technology. Sequences were trimmed using Cutadapt. To determine XIST sequence, a de novo assembly was made using a mapping with STAR and StringTie. The GTF observed was aligned in IGV and XIST sequence was found at the known position. XIST sequence was included in EquCab 3.0 and a new mapping was performed using STAR. Counting was performed using featureCounts. The differential analysis was performed including embryo sex and diameter using Deseq2. A functional enrichment analysis of DEGs was performed using PANTHER web software. In another time, counts of all genes were normalized using Deseq2 and a gene set analysis (GSEA) was performed using the Gene Ontology (GO) biological process (BP) and the Kyoto Encyclopedia of Genes and Genomes (KEGG) databases with GSEA software from the
Broad Institute. To reduce redundancy between terms, SUMER analysis was performed on enriched gene sets.

**Additional file 3: Supplementary Table 1.** Sex, diameter, RNA yield and mare age group of each equine embryo at Day 8 post-ovulation. Embryo ID, mare age group, embryo sex and diameter and RNA yield for each embryo are represented in columns while each line represents one individual embryo. OM: old mares; YM: young mares.

**Additional file 4: Supplementary Table 2.** Differential gene analysis using DeSeq2 in ICMandTE of equine embryo at Day 8 post-ovulation. Equine ensemble ID, orthologue human Ensembl ID, Orthologue human Entrez Gene ID, gene description, normalized counts for each sample of ICM and parameters obtained after Deseq2 analysis (log2FoldChange, pvalue and padj (after FDR correction)) of genes expressed in ICM enriched part (before gene expression deconvolution using DeMixT) of OM and YM embryos. Young group was the reference group. ICM: Inner cell mass; OM: old mares; YM: young mares.

**Additional file 5: Supplementary Table 3.** Gene set enrichment analysis results on gene expression comparing ICMandTE of embryos from old mares to young mares. Gene Set Enrichment Analysis results (pathway name, GO accession when possible and size, Normalized Enrichment Score, *p*-value and FDR corrected *q*-value) for GO biological process and KEGG databases on ICMandTE gene expression table. Young group was the reference group.

**Additional file 6: Supplementary Table 4.** Differential gene analysis using DeSeq2 in DeMixT_ICMells of equine embryo at Day 8 post-ovu-lation. Equine ensemble ID, orthologue human Ensembl ID, Orthologue human Entrez Gene ID, gene description, normalized counts for each sample of ICM and parameters obtained after Deseq2 analysis (log2Fold-Change, pvalue and padj (after FDR correction)) of genes expressed in ICM (after gene expression deconvolution of ICMandTE using DeMixT) of OM and YM embryos. Young group was the reference group. ICM: Inner cell mass; OM: old mares; YM: young mares.

**Additional file 7: Supplementary Table 5.** Differential gene analysis using DeSeq2 in TE_part of equine embryo at Day 8 post-ovulation. Equine ensemble ID, orthologue human Ensembl ID, Orthologue human Entrez Gene ID, gene description, normalized counts for each sample of ICM and parameters obtained after Deseq2 analysis (log2FoldChange, pvalue and padj (after FDR correction)) of genes expressed in TE of OM and YM embryos. Young group was the reference group. TE: Trophoblast; OM: old mares; YM: young mares.

**Additional file 8: Supplementary Table 6.** Gene set enrichment analysis results on gene expression of DeMixT_ICM ``cells of embryos from young and old mares. Gene Set Enrichment Analysis results (pathway name, GO accession when possible and size, Normalized Enrichment Score, *p*-value and FDR corrected *q*-value) for GO biological process and KEGG databases on DeMixT_ICM ``cells gene expression table (after gene expression deconvolution on ICMandTE using DeMixT). Young group was the reference group.

**Additional file 9: Supplementary Table 7.** Gene set enrichment analysis results on gene expression of TE_part of embryos from young and old mares. Gene Set Enrichment Analysis results (pathway name, GO accession when possible and size, Normalized Enrichment Score, *p*-value and FDR corrected *q*-value) for GO biological process and KEGG databases on TE_part gene expression table. Young group was the reference group.

## Acknowledgments

The authors are grateful to the staff of the Institut Français du Cheval et de l'Equitation (IFCE) experimental farm (Domaine de la Valade, Chamberet, France) for care and management of animals. We acknowledge the High-throughput sequencing facility of I2BC for its sequencing and bioinformatics expertise. The bioinformatics analyses were performed thanks to Core Cluster of the Institut Français de Bioinformatique (IFB) (ANR-11-INBS-0013). Many thanks to Matthias Zytnicki and Christophe Köpp for their advice on RNA-seq de novo analysis. Many thanks to Pablo Ross who kindly provided the coordinates for the XIST gene. The authors thank Shavahn Loux for her information

on PLAC8A. The authors are grateful to Alice Jouveau and Sophie Calderari who kindly revised this article.

## Authors' contributions

PCP obtained the funding. PCP and VD conceived the project. LW, PCP, and VD supervised the study. ED, CD, CA, ND, NP, LW, VD and PCP adapted the methodology for the project. ED, CD, CA, YJ, JAR, and LW performed the experiments. CD, CA, ND, NP, MD and LW provided the resources. ED, LJ, YJ and RL performed data curation. ED and LJ analysed the data. ED wrote the original draft. All authors read, revised, and approved the submitted manuscript.

## Funding

This work was supported by the "Institut Français du Cheval et de l'Equitation" (grant numbers CS_2018_23, 2018). The National Research Institute for Agriculture, Food and Environment (INRAE) department Animal Physiology and Breeding Systems also supported this research.

## Availability of data and materials

The RNA sequencing data supporting the conclusions of this article are available in the repository [accession: GSE162893; https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE162893] from the GEO SuperSeries [accession: GSE193676; https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE193676].

## Declarations

### Ethics approval and consent to participate

The experiment was performed at the experimental farm of IFCE (research agreement C1903602 valid until March 22, 2023). The protocol was approved by the local animal care and use committee ("Comité des Utilisateurs de la Station Expérimentale de Chamberet") and by the regional ethical committee ("Comité Régional d'Ethique pour l'Expérimentation Animale du Limousin", approved under N° C2EA - 33 in the National Registry of French Ethical Committees for animal experimentation) under protocol number APAFIS#14963-2018050316037888 v2. All experiments were performed in accordance with the European Union Directive 2010/63EU.

### Consent for publication

Not applicable.

### Competing interests

The authors declare that they have no competing interests.

### Author details

Received: 7 February 2022 Accepted: 22 April 2022

Published online: 15 June 2022

## References

1. Scoggin CF. Not just a number: effect of age on fertility, pregnancy and offspring vigour in thoroughbred brood-mares. Reprod Fertil Dev. 2015;27(6):872.

2. Allen WR, Brown L, Wright M, Wilsher S. Reproductive efficiency of Fla-trace and National Hunt Thoroughbred mares and stallions in England. Equine Veterinary J. 2007;39(5):438-45.

3. Langlois B, Blouin C. Statistical analysis of some factors affecting the number of horse births in France. Reprod Nutr Dev. 2004;44(6):583-95.

4. Malheiros de Souza JR, Dias Gonçalves PB, Bertolin K, Ferreira R, Sardinha Ribeiro AS, Ribeiro DB, et al. Age-dependent effect of foal heat breeding on pregnancy and embryo mortality rates in Thoroughbred mares. J Equine Veterinary Sci. 2020;5:102982.

<footnote n="1">¹Université Paris-Saclay, UVSQ, INRAE, BREED, 78350 Jouy-en-Josas, France.</footnote>

<footnote n="2">²Ecole Nationale Vétérinaire d'Alfort, BREED, 94700 Maisons-Alfort, France.</footnote>

<footnote n="3">³IFCE, Plateau technique de Chamberet, 19370 Chamberet, France. ⁴Institute for Integrative Biology of the Cell (I2BC), UMR 9198 CNRS, CEA, Paris-Sud University F-91198, Gif-sur-Yvette, France. ⁵Institut Pasteur—Bioinformatics and Biostatistics Hub—Department of Computational Biology, Paris, France.</footnote>
5. Chevalier-Clément F. Pregnancy loss in the mare. Anim Reprod Sci. 1989;20(3):231-44.

6. Ball BA, Little TV, Hillman RB, Woods GL. Pregnancy rates at days 2 and 14 and estimated embryonic loss rates prior to day 14 in normal and subfertile mares. Theriogenology. 1986;26(5):611-9.

7. Doig PA, McKnight JD, Miller RB. The use of endometrial biopsy in the infertile Mare. Can Vet J. 1981;22:72-6.

8. Siemieniuch MJ, Gajos K, Kozdrowski R, Nowak M. Advanced age in mares affects endometrial secretion of arachidonic acid metabolites during equine subclinical endometritis. Theriogenology. 2017;103:191-6.

9. Carnevale EM, Ginther OJ. Relationships of age to uterine function and reproductive efficiency in mares. Theriogenology. 1992;37(5):1101-15.

10. Ball BA, Little TV, Weber JA, Woods GL. Survival of Day-4 embryos from young, normal mares and aged, subfertile mares after transfer to normal recipient mares. Reproduction. 1989;85(1):187-94.

11. Brinsko SP, Ball BA, Ellington JE. In vitro maturation of equine oocytes obtained from different age groups of sexually mature mares. Theriogenology. 1995;44(4):461-9.

12. Carnevale E, Uson M, Bozzola J, King S, Schmitt S, Gates H. Comparison of oocytes from young and old mares with light and electron microscopy. Theriogenology. 1999;51(1):299.

13. Frank BL, Doddman CD, Stokes JE, Carnevale EM. Association of equine oocyte and cleavage stage embryo morphology with maternal age and pregnancy after intracytoplasmic sperm injection. Reprod Fertil Dev. 2019;31(12):1812.

14. Carnevale EM, Coutinho da Silva MA, Panzani D, Stokes JE, Squires EL. Factors affecting the success of oocyte transfer in a clinical program for subfertile mares. Theriogenology. 2005;64(3):519-27.

15. Altermatt JL, Suh TK, Stokes JE, Carnevale EM. Effects of age and equine follicle-stimulating hormone (eFSH) on collection and viability of equine oocytes assessed by morphology and developmental competency after intracytoplasmic sperm injection (ICSI). Reprod Fertil Dev. 2009;21(4):615.

16. Campos-Chillon F, Farmerie TA, Bouma GJ, Clay CM, Carnevale EM. Effects of aging on gene expression and mitochondrial DNA in the equine oocyte and follicle cells. Reprod Fertil Dev. 2015;27(6):925.

17. Cox L, Vanderwall DK, Parkinson KC, Sweat A, Isom SC. Expression profiles of select genes in cumulus-oocyte complexes from young and aged mares. Reprod Fertil Dev. 2015;27(6):914.

18. Rizzo M, Stout TAE, Cristarella S, Quartuccio M, Kops GJPL, De Ruijter-Villani M. Compromised MPS1 Activity Induces Multipolar Spindle Formation in Oocytes From Aged Mares: Establishing the Horse as a Natural Animal Model to Study Age-Induced Oocyte Meiotic Spindle Instability. Front Cell Dev Biol. 2021;9:657366.

19. Rizzo M, du Preez N, Ducheyne KD, Deelen C, Beitsma MM, Stout TAE, et al. The horse as a natural model to study reproductive aging-induced aneuploidy and weakened centromeric cohesion in oocytes. Aging. 2020; [cité 9 nov 2020]; Disponible sur: http://www.aging-us.com/article/104159/text.

20. Rizzo M, Ducheyne KD, Deelen C, Beitsma M, Cristarella S, Quartuccio M, et al. Advanced mare age impairs the ability of in vitro-matured oocytes to correctly align chromosomes on the metaphase plate. Equine Vet J. 2019;51(2):252-7.

21. Rambags BPB, van Boxtel DCJ, Tharasanit T, Lenstra JA, Colenbrander B, Stout TAE. Advancing maternal age predisposes to mitochondrial damage and loss during maturation of equine oocytes in vitro. Theriogenology. 2014;81(7):959-65.

22. Catandi GD, Obeidat YM, Broeckling CD, Chen TW, Chicco AJ, Carnevale EM. Equine maternal aging affects oocyte lipid content, metabolic function and developmental potential. Reproduction. 2021;161(4):399-409.

23. Hendriks WK, Colleoni S, Galli C, Paris DBBP, Colenbrander B, Roelen BAJ, et al. Maternal age and in vitro culture affect mitochondrial number and function in equine oocytes and embryos. Reprod Fertil Dev. 2015;27(6):957.

24. Woods GL, Baker CB, Hillman RB, Schlafer DH. Recent studies relating to early embryonic death in the mare. Equine Vet J. 1985;17(S3):104-7.

25. Carnevale EM, Griffin PG, Ginther OJ. Age-associated subfertility before entry of embryos into the uterus in mares. Equine Vet J. 1993;25(S15):31-5.

26. Brinsko SP, Ball BA, Miller PG, Thomas PGA, Ellington JE. In vitro development of day 2 embryos obtained from young, fertile mares and aged, subfertile mares. Reproduction. 1994;102(2):371-8.

27. Willmann C, Schuler G, Hoffmann B, Parvizi N, Aurich C. Effects of age and altrenogest treatment on conceptus development and secretion of LH, progesterone and eCG in early-pregnant mares. Theriogenology. 2011;75(3):421-8.

28. Camargo Ferreira J, Linhares Boakari Y, Sousa Rocha N, Saules Ignácio F, Barbosa da Costa G, de Meira C. Luteal vascularity and embryo dynamics in mares during early gestation: effect of age and endometrial degeneration. Reprod Dom Anim. 2019;54(3):571-9.

29. Morel MCGD, Newcombe JR, Swindlehurst JC. The effect of age on multiple ovulation rates, multiple pregnancy rates and embryonic vesicle diameter in the mare. Theriogenology. 2005;63(9):2482-93.

30. Vanderwall DK. Early embryonic development and evaluation of equine embryo viability. Vet Clin N Am Equine Pract. 1996;12(1):61-83.

31. Panzani D, Rota A, Marmorini P, Vannozzi I, Camillo F. Retrospective study of factors affecting multiple ovulations, embryo recovery, quality, and diameter in a commercial equine embryo transfer program. Theriogenology. 2014;82(6):807-14.

32. Goszczynski DE, Tinetti PS, Choi YH, Hinrichs K, Ross PJ. Genome activation in equine in vitro-produced embryos. Biol Reprod. 2022;106(1):66-82.

33. Grondahl C, Hyttel P. Nucleogenesis and ribonucleic acid synthesis in preimplantation equine embryos. Biol Reprod. 1996;55(4):769-74.

34. Carnevale EM, Bergfelt DR, Ginther OJ. Aging effects on follicular activity and concentrations of FSH, LH, and progesterone in mares. Anim Reprod Sci. 1993;31(3-4):287-99.

35. Chavatte-Palmer P, Velazquez MA, Jammes H, Duranthon V. Review: epigenetics, developmental programming and nutrition in herbivores. Animal. 2018;12(s2):s363-71.

36. McKinnon AO, Squires EL. Morphologic assessment of the equine embryo. J Am Vet Med Assoc. 1988;192(3):401-6.

37. Iqbal K, Chitwood JL, Meyers-Brown GA, Roser JF, Ross PJ. RNA-Seq Transcriptome Profiling of Equine Inner Cell Mass and Trophectoderm. Biol Reprod. 2014;90(3) [cité 2 août 2019] Disponible sur: https://academicroup.com/biolreprod/article-lookup/doi/10.1095/biolreprod.113.113928.

38. Cao S, Wang JR, Ji S, Yang P, Chen J, Shen JP, et al. Differing total mRNA expression shapes the molecular and clinical phenotype of cancer. bioRxiv. 2020;57. https://doi.org/10.21203/rs.3.rs-600171/v1.

39. Wang Z, Cao S, Morris JS, Ahn J, Liu R, Tyekucheva S, et al. Transcriptome Deconvolution of Heterogeneous Tumor Samples with Immune Infiltration. iScience. 2018;9:451-60.

40. Savage SR, Shi Z, Liao Y, Zhang B. Graph algorithms for condensing and consolidating gene set Analysis results. Mol Cell Proteomics. 2019;18(8):S141-52.

41. Shannon P. Cytoscape: A Software Environment for Integrated Models of Biomolecular Interaction Networks. Genome Res. 2003;13(11):2498-504.

42. Marinone AI, Losinno L, Fumuso E, Rodríguez EM, Redolatti C, Cantatore S, et al. The effect of mare's age on multiple ovulation rate, embryo recovery, post-transfer pregnancy rate, and interovulatory interval in a commercial embryo transfer program in Argentina. Anim Reprod Sci. 2015;158:53-9.

43. Betteridge KJ, Eaglesome MD, Mitchellt D, Flood PF, Beriault R. Development of horse embryos up to twenty two days after ovulation: observations on fresh specimens. J Anat. 1982;135(Pt 1):191-209.

44. Beckelmann J, Budik S, Bartel C, Aurich C. Evaluation of Xist expression in preattachment equine embryos. Theriogenology. 2012;78(7):1429-36.

45. Colchen S, Battut I, Fiéni F, Tainturier D, Siliart B, Bruyas JF. Quantitative histological analysis of equine embryos at exactly 156 and 168 h after ovulation. J Reprod Fertil Suppl. 2000;56:527-37.

46. Leisinger CA, Medina V, Markle ML, Paccamonti DL, Pinto CRF. Morphological evaluation of day 8 embryos developed during induced aluteal cycles in the mare. Theriogenology. 2018;105:178-83.

47. Hisey E, Ross PJ, Meyers SA. A review of OCT4 functions and applications to equine embryos. J Equine Veterinary Sci. 2021;98:103364.
48. Enders AC, Schlafke S, Lantz KC, Liu IKM. Endoderm cells of the equine yolk sac from day 7 until formation of the definitive yolk sac placenta. Equine Vet J. 1993;25(S15):3-9.

49. Enders AC, Lantz KC, Liu IKM, Schlafke S. Loss of polar trophoblast during differentiation of the blastocyst of the horse. J Reprod Fertil. 1988;83(1):447-60.

50. Grøndahl ML, Andersen CY, Bogstad J, Nielsen FC, Meinertz H, Borup R. Gene expression profiles of single human mature oocytes in relation to age. Hum Reprod. 2010;25(4):957-68.

51. Reyes JM, Silva E, Chitwood JL, Schoolcraft WB, Krisher RL, Ross PJ. Differing molecular response of young and advanced maternal age human oocytes to IVM. Human Reprod. 2017;32(11):2199-208.

52. Hamatani T, Falco G, Carter MG, Akutsu H, Stagg CA, Sharov AA, et al. Age-associated alteration of gene expression patterns in mouse oocytes. Human Mol Genet. 2004;13(19):2263-78.

53. Kawai K, Harada T, Ishikawa T, Sugiyama R, Kawamura T, Yoshida A, et al. Parental age and gene expression profiles in individual human blastocysts. Sci Rep. 2018;8(1):2380.

54. McCallie BR, Parks JC, Trahan GD, Jones KL, Coate BD, Griffin DK, et al. Compromised global embryonic transcriptome associated with advanced maternal age. J Assist Reprod Genet. 2019;36(5):915-24.

55. Hauck CR, Hsia DA. Schlaepfer DD. The Focal Adhesion Kinase--A Regulator of Cell Migration and Invasion. IUBMB Life. 2002;53(2):115-9.

56. Provenzano PP, Keely PJ. Mechanical signaling through the cytoskeleton regulates cell proliferation by coordinated focal adhesion and Rho GTPase signaling. J Cell Sci. 2011;124(8):1195-205.

57. Zhao X, Guan J-L. Focal adhesion kinase and its signaling pathways in cell migration and angiogenesis. Adv Drug Delivery Rev. 2011;63(8):610-5.

58. Hynes RO. Integrins: bidirectional, allosteric signaling machines. Cells. 2002;110(6):673-87.

59. Cheuk IW-Y, Siu MT, Ho JC-W, Chen J, Shin VY. ITGAV targeting as a therapeutic approach for treatment of metastatic breast cancer. Am J Cancer Res. 2020;10(1):211-23.

60. Riento K, Ridley AJ. ROCKs: multifunctional kinases in cell behaviour. Nat Rev Mol Cell Biol. 2003;4(6):446-56.

61. Brouns MR. p190 RhoGAP is required for neural development. Development. 2000;127:4891-903.

62. Rane CK, Minden A. P21 activated kinases: structure, regulation, and functions. Small GTPases. 2014;5(1):e28003.

63. Chan PM, Lim L, Manser E. PAK Is Regulated by PI3K, PIX, CDC42, and PP2Cα and Mediates Focal Adhesion Turnover in the Hyperosmotic Stress-induced p38 Pathway. J Biol Chem. 2008;283(36):24949-61.

64. Bacon C, Endris V, Rappold GA. The cellular function of srGAP3 and its role in neuronal morphogenesis. Mechanisms Dev. 2013;130(6-8):391-5.

65. Wu MY, Hill CS. TGF-β superfamily signaling in embryonic development and homeostasis. Dev Cell. 2009;16(3):329-43.

66. Kemp C, Willems E, Abdo S, Lambiv L, Leyns L. Expression of all Wnt genes and their secreted antagonists during mouse blastocyst and postimplantation development. Dev Dyn. 2005;233(3):1064-75.

67. Smith CA, Hurley TM, McClive PJ, Sinclair AH. Restricted expression of DMRT3 in chicken and mouse embryos. Mech Dev. 2002;119:573-6.

68. Zhang Y, Huang L, Fu H, Smith OK, Lin CM, Utani K, et al. A replicator-specific binding protein essential for site-specific initiation of DNA replication in mammalian cells. Nat Commun. 2016;7(1):11748.

69. Yuan K, O'Farrell PH. Cyclin B3 is a mitotic cyclin that promotes the metaphase-anaphase transition. Curr Biology. 2015;25(6):811-6.

70. Vinopal S, Černohorská M, Sulimenko V, Sulimenko T, Vosecká V, Flemgr M, et al. γ-Tubulin 2 Nucleates Microtubules and Is Downregulated in Mouse Early Embryogenesis. PLoS ONE. 2012;7(1):e29919.

71. Fry AM, O'Regan L, Sabir SR, Bayliss R. Cell cycle regulation by the NEK family of protein kinases. J Cell Sci. 2012;125(19):4423-33.

72. Losada A, Hirano T. Shaping the metaphase chromosome: coordination of cohesion and condensation. Bioessays. 2001;23(10):924-35.

73. Schmidt Tl, Kleylein-Sohn J, Westendorf J, Le Clech M, Lavoie SB, Stierhof Y-D, et al. Control of centriole length by CPAP and CP110. Curr Biol. 2009;19(12):1005-11.

74. Watanabe S, Shioi G, Furuta Y, Goshima G. Intra-spindle microtubule assembly regulates clustering of microtubule-organizing centers during early mouse development. Cell Reports. 2016;15(1):54-60.

75. Ladurner R, Kreidl E, Ivanov MP, Ekker H, Idarraga-Amado MH, Busslinger GA, et al. Sororin actively maintains sister chromatid cohesion. EMBO J. 2016;35(6):635-53.

76. Rivera T, Losada A. Shugoshin and PP2A, shared duties at the centromere. Bioessays. 2006;28(8):775-9.

77. Welburn JPI. The molecular basis for kinesin functional specificity during mitosis. Cytoskeleton. 2013;70(9):476-93.

78. Zhou C-J, Wang X-Y, Han Z, Wang D-H, Ma Y-Z, Liang C-G. Loss of CENPF leads to developmental failure in mouse embryos. Cell Cycle. 2019;18(20):2784-99.

79. Monzo P, Gauthier NC, Keslair F, Loubat A, Field CM, Marchand-Brustel YL, et al. Clues to CD2-associated protein involvement in cytokinesis. *Trends in Mol Biol. 2005;16:12.

80. Babst M, Davies BA. Regulation of Vps4 during MVB sorting and cytokinesis. Traffic. 2011;12:1298-305.

81. Hales BF, Grenier L, Lalancette C, Robaire B. Epigenetic programming: from gametes to blastocyst. Birth Defects Res Part A Clin Mol Teratol. 2011;91(8):652-65.

82. Marcho C, Cui W, Mager J. Epigenetic dynamics during preimplantation development. Reproduction. 2015;150(3):R109-20.

83. Suo L, Meng Q-G, Pei Y, Yan C-L, Fu X-W, Bunch TD, et al. Changes in acetylation on lysine 12 of histone H4 (acH4K12) of murine oocytes during maternal aging may affect fertilization and subsequent embryo development. Fertil Steril. 2010;93(3):945-51.

84. Dumollard R, Duchen M, Carroll J. The role of mitochondrial function in the oocyte and embryo. In: The mitochondrion in the germline and early development: Elsevier; 2007. p. 21-49. [cited 18 sept 2020] (Current topics in developmental biology; vol. 77). Disponible sur: https://linkinghub.elsevier.com/retrieve/pii/S0070215306770028.

85. Balaban RS, Nemoto S, Finkel T. Mitochondria, oxidants, and aging. Cell. 2005;120(4):483-95.

86. Wilding M, Dale B, Marino M, di Matteo L, Alviggi C, Pisaturo ML, et al. Mitochondrial aggregation patterns and activity in human oocytes and preimplantation embryos. Human Reproduction. 2001;16(5):909-17.

87. Gardner DK, Harvey AJ. Blastocyst metabolism. Reprod Fertil Dev. 2015;27(4):638.

88. Lane M, O'Donovan MK, Squires EL, Seidel GE, Gardner DK. Assessment of metabolism of equine morulae and blastocysts. Mol Reprod Dev. 2001;59(1):33-7.

89. Hutchison CA, Newbold JE, Potter SS, Edgell MH. Material inheritance of mammalian mitochondrial DNA. Nature. 1974;251(5475):536-8.

90. Krisher RL. Maternal age affects oocyte developmental potential at both ends of the age spectrum. Reprod Fertil Dev. 2019;31(1):1.

91. Hendriks WK, Colleoni S, Galli C, Paris DBBP, Colenbrander B, Stout TAE. Mitochondrial DNA replication is initiated at blastocyst formation in equine embryos. Reprod Fertil Dev. 2019;31(3):570.

92. Waelchli RO, MacPhee DJ, Kidder GM, Betteridge KJ. Evidence for the presence of sodium- and potassium-dependent adenosine Triphosphatase alpha 1 and beta 1 subunit isoforms and their probable role in blastocyst expansion in the Preattachment horse conceptus. Biol Reprod. 1997;57:630-40.

93. do Nascimento AD, JCC M, ARR C, Batista AM, Kastelic JP, Câmara DR. Inhibition of Na+, K+ - ATPase with ouabain is detrimental to equine blastocysts. Anim Reprod. 2020;17(1):e20190079.

94. Arnaud MG, Dubroeucq MH. Rivot MD. Notation de l'état corporel des chevaux de selle et de sport: guide pratique. 1ère ed. Paris: Institut de l'élevage; 1997. p. 40.

95. Duchamp G, Bour B, Combarnous Y, Palmer E. Alternative solutions to hCG induction of ovulation in the mare. J Reprod Fertil Suppl. 1987;35:221-8.

96. Bucca S, Carli A. Efficacy of human chorionic gonadotropin to induce ovulation in the mare, when associated with a single dose of dexamethasone administered at breeding time: efficacy of human chorionic gonadotropin to induce ovulation when associated with dexamethasone. Equine Vet J. 2011;43:32-4.

97. Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBNet J. 2011;17(1):10-2.

98. Gutiérrez-Adán A, Perez-Crespo M, Fernandez-Gonzalez R, Ramirez M, Moreira P, Pintado B, et al. Developmental consequences of sexual
dimorphism during pre-implantation embryonic development. Reprod Domest Anim. 2006;41(s2):54-62.

99. Pertea M, Pertea GM, Antonescu CM, Chang T-C, Mendell JT, Salzberg SL. StringTie enables improved reconstruction of a transcriptome from RNA-seq reads. Nat Biotechnol. 2015;33(3):290-5.

100. Robinson JT, Thorvaldsdóttir H, Winckler W, Guttman M, Lander ES, Getz G, et al. Integrative genomics viewer. Nat Biotechnol. 2011;29(1):24-6.

101. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29(1):15-21.

102. Liao Y, Smyth GK, Shi W, et al. Bioinformatics. 2014;30(7):923-30.

103. R Core Team. A Language and Environment for Statistical Computing. Vienna: R foundation for Statistical Computing; 2020. Disponible sur: https://www.R-project.org/

104. RStudio Team. RStudio: integrated development for R. RStudio. Boston: PBC; 2020. Disponible sur: http://www.rstudio.com/

105. Pinheiro J, Bates D, Debroy S, Sarkar D, R Core Team. nlme: Linear and Nonlinear Mixed Effects Models. 2020. Disponible sur: https://CRAN.R-project.org/package=nlme

106. Giraudoux P. pgirmess: Spatial Analysis and Data Mining for Field Ecologists. 2018. Disponible sur: https://CRAN.R-project.org/package=pgirmess

107. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550.

108. Lê S, Josse J, Husson F. FactoMineR: An R package for multivariate analysis. J Stat Soft. 2008;25(1):18.

109. Gómez-Rubio V. ggplot2 - Elegant Graphics for Data Analysis (2nd Edition). J Stat Soft. 2017;77 [cité 9 août 2020] (Book Review 2). Disponible sur: http://www.jstatsoft.org/v77/b02/.

110. Kassambara A, Mundt F. factoextra: Extract and Visualize the Results of Multivariate Data Analyses. 2020. Disponible sur: https://CRAN.R-project.org/package=factoextra

111. Oksanen J, Blanchet FG, Friendly M, Kindt R, Legendre P, McGlinn D, et al. vegan: Community Ecology Package. 2020. Disponible sur: https://CRAN.R-project.org/package=vegan

112. Kolberg L, Raudvere U. gprofiler2: Interface to the "g:Profiler" Toolset. 2020. Disponible sur: https://CRAN.R-project.org/package=gprofiler2

113. Yates AD, Achutan P, Akanni W, Allen J, Allen J, Alvarez-Jarreta J, et al. Ensembl 2020.

114. Mi H, Ebert D, Muruganujan A, Mills C, Albou L-P, Mushayamaha T, et al. PANTHER version 16: a revised family classification, tree-based classification tool, enhancer regions and extensive API. Nucleic Acids Res. 2021;49(D1):D394-403.

115. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci. 2005;102(43):15545-50.

116. Mootha VK, Lindgren CM, Eriksson K-F, Subramanian A, Sihag S, Lehar J, et al. PGC-1α-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes. Nat Genet. 2003;34(3):267-73.

117. Liberzon A, Birger C, Thorvaldsdóttir H, Ghandi M, Mesirov JP, Tamayo P. The molecular signatures database Hallmark gene set collection. Cell Systems. 2015;1(6):417-25.

## Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Ready to submit your research? Choose BMC and benefit from:

* fast, convenient online submission

* thorough peer review by experienced researchers in your field

* rapid publication on acceptance

* support for research data, including large and complex data types

* gold Open Access which fosters wider collaboration and increased citations

* maximum visibility for your research: over 100M website views per year

At BMC, research is always in progress.

Learn more biomedcentral.com/submissions

<figure><img src="image_10.png" /></figure>