The Oviductal Extracellular Vesicles’ RNA Cargo Regulates the Bovine Embryonic Transcriptome

Oviductal extracellular vesicles (oEVs) are emerging as key players in the gamete/embryo–oviduct interactions that contribute to successful pregnancy. Various positive effects of oEVs on gametes and early embryos have been found in vitro. To determine whether these effects are associated with changes of embryonic gene expression, the transcriptomes of embryos supplemented with bovine fresh (FeEVs) or frozen (FoEVs) oEVs during in vitro culture compared to controls without oEVs were analyzed by low-input RNA sequencing. Analysis of RNA-seq data revealed 221 differentially expressed genes (DEGs) between FoEV treatment and control, 67 DEGs for FeEV and FoEV treatments, and minor differences between FeEV treatment and control (28 DEGs). An integrative analysis of mRNAs and miRNAs contained in oEVs obtained in a previous study with embryonic mRNA alterations pointed to direct effects of oEV cargo on embryos (1) by increasing the concentration of delivered transcripts; (2) by translating delivered mRNAs to proteins that regulate embryonic gene expression; and (3) by oEV-derived miRNAs which downregulate embryonic mRNAs or modify gene expression in other ways. Our study provided the first high-throughput analysis of the embryonic transcriptome regulated by oEVs, increasing our knowledge on the impact of oEVs on the embryo and revealing the oEV RNA components that potentially regulate embryonic development.


Introduction
Extracellular vesicles (EVs) are well recognized mediators of cell-to-cell communication [1], a function they carry out by transferring their bioactive molecular cargo (RNAs, proteins, lipids, metabolites, and genomic DNA) to recipient cells [2,3]. Although at least three different types of EVs have been described based on their biogenesis and physical characteristics-exosomes, microvesicles and apoptotic vesicles [4]-only the first two types have attracted much attention in recent years, due to their contribution to a wide range of physiological and pathological processes such as angiogenesis, cell survival, modulation of the immune response, inflammation, and cancer, as well as embryonic development [5,6]. In fact, EVs identified in the oviduct and in the uterus have emerged as key players in the embryo-maternal dialogue contributing to successful pregnancy [7][8][9].
In particular, the potential role of oviductal EVs (oEVs) has received growing attention in recent years, since the oviduct is the place that hosts and supports the first reproductive events [10,11], and oEVs could be key modulators of such events. To date, EVs have been identified in the oviduct of different species (bovine, mouse, porcine, avian, and turtle) and their functional effects have been studied in gametes and embryos (reviewed in Almiñana and Bauersachs [12]). For example, it has been shown that oviductal EVs (oEVs) support bovine embryonic development [13,14], canine oocyte maturation [15], modulate sperm capacitation and sperm fertilizing ability in the mouse and in the cat [11,16], and regulate polyspermy fertilization in the pig [17]. Regarding the effects of oEVs on embryonic development, our laboratory previously demonstrated that oEVs are taken up by the bovine embryo during in vitro culture, and that the supplementation of oEVs during in vitro embryo culture improved embryonic development and quality in terms of blastocyst rates, cell number, and hatching rates [13]. Moreover, we showed that frozen and fresh oEVs had different effects on embryonic development and quality [13]. Along the same lines, Lopera-Vásquez [14] reported that oEVs enhanced embryo cryosurvival. Furthermore, Lopera-Vásquez [14,18] showed that oEV supplementation during in vitro culture altered the expression of a few genes involved in embryonic development, metabolism, and epigenetic regulation, making the embryos more similar to their in vivo counterparts [14]. These two studies by Lopera-Vásquez et al. [14,18] provided a few hints about the potential role of oEVs in modulating embryonic gene expression by using a targeted RT-qPCR approach, and called for an in-depth analysis of the impact of oEVs on the embryonic transcriptome.
Given the wide range of oEV components recently identified in our laboratory (mRNAs, proteins, ncRNAs including miRNAs, snoRNAs, snRNAs, and metabolites) [19,20], it is difficult to select potential candidates as modulators of embryonic development. To date, only a few miRNAs and proteins have been proven to be responsible for oEVs' functional effects on spermatozoa [11,16], while the functional impact of the oEV cargo on embryos and the extent of those effects is not yet fully understood.
Therefore, in the present study we aimed to demonstrate that the RNA cargo in oEVs regulates early embryonic development by altering the embryonic transcriptome. We hypothesized that oEVs bring RNA components (mRNAs and miRNAs) and proteins into the embryo and thus alter the embryonic transcriptome. Moreover, we propose different modes of action by which the RNA cargo of oEVs could modify the embryonic transcriptome: (1) oEV-derived mRNAs could be incorporated into embryos via EVs and thereby increase the concentration of the delivered transcripts; (2) delivered mRNAs could be translated and the corresponding proteins could lead to regulation of embryonic gene expression; and (3) oEV-derived miRNAs and other ncRNAs could act by targeting embryonic mRNAs, and thus downregulate or modify embryonic gene expression in other ways (e.g., mRNA isoform expression, indirect effects on gene expression). In addition, given the differential effect of fresh and frozen oEVs on embryonic development, as demonstrated in our previous study [13], we also hypothesized that frozen and fresh oEVs affect the embryonic transcriptome differently. To this end, we used a low-input RNA-seq approach to profile the transcriptional responses of embryos cultured in vitro with fresh and frozen oEVs and controls without the addition of oEVs. Subsequently, to unveil the potential oEV RNA components capable of regulating the embryonic development, we performed an integrative analysis of mRNA and miRNA cargo identified in oEVs [19] and the embryonic transcriptome alterations induced by oEVs. The knowledge derived from our study will lead to a more meaningful understanding of the impact of oEVs on the embryo, while revealing the oEV RNA cargo potentially involved in the regulation of embryonic development.

Oviductal EV Supplementation during In Vitro Embryo Culture Altered the Embryonic Transcriptome
The oEVs used in this study were derived from cows in the postovulatory stage. In our previous studies, characterization of these oEVs revealed a population of small extracellular vesicles (30-100 nm) resembling exosomes (50%-60% of all vesicles) and a population of larger extracellular vesicles (>100 nm) resembling microvesicles (25%-30% 100-150 nm, 10% 150-200 nm) [13,19]. Typical EV marker proteins were detected by Western blotting, such as HSP70, ANXA1, MYH9, and HSPA8 [13,19]. In the present study, we focused on analysis of the transcriptome of embryos cultured in vitro with or without oEV supplementation. Hierarchical cluster analysis of the DEGs obtained for the three comparisons is shown in Figure  2. The main differences were identified between FoEV treatment and Co (221 DEGs; Figures 1 and  2A), the lowest for the comparison between FeEV and Co treatments (28 DEGs; Figures 1 and 2B), and 67 DEGs for the comparison between EV treatments (Figures 1 and 2C). Furthermore, when an additional cut-off was considered (minimum fold change of 1.5), 112, 19, and 46 DEGs were obtained for comparisons between FoEV treatment and Co, FeEV treatment and Co, and between EV treatments, respectively (Supplementary Data S1 (Tables S2-S4)). Hierarchical cluster analysis of the DEGs obtained for the three comparisons is shown in Figure 2. The main differences were identified between FoEV treatment and Co (221 DEGs; Figures 1 and  2A), the lowest for the comparison between FeEV and Co treatments (28 DEGs; Figures 1 and 2B), and 67 DEGs for the comparison between EV treatments (Figures 1 and 2C). Furthermore, when an additional cut-off was considered (minimum fold change of 1.5), 112, 19, and 46 DEGs were obtained for comparisons between FoEV treatment and Co, FeEV treatment and Co, and between EV treatments, respectively (Supplementary Data S1 (Tables S2-S4)). Color scale is from -2 (blue, lower than mean) to 2 (red, higher than mean).

Functional Annotation Revealed Potential Regulatory Pathways Affected by Oviductal EVs
To obtain a more meaningful view of how oEVs can modulate the embryo transcriptome, gene ontology terms (GO) terms and pathways for DEGs from the three different comparisons (FoEV treatment vs. Co; FeEV treatment vs. Co, and FoEV vs. FeEv treatment) were analyzed using functional annotation databases in Metascape and with Ingenuity Pathway Analysis (IPA) tools.
In Figure 3, a Metascape heatmap plot represents the functional enrichment analysis for up-and downregulated genes in embryos treated with frozen or fresh EVs compared to control and between EV treatments, and shows clusters of enriched functional terms across DEGs from different comparisons. The heatmap illustrates that functional terms and pathways related to "embryo morphogenesis" were only enriched for genes downregulated in FoEV treatment vs. Co and FeEV treatment vs. Co sets. Furthermore, genes downregulated in FoEV treatment vs. Co were also enriched for biological processes related to "female gamete generation", "apoptosis", "response to external stimulus", "response to osmotic stress", "telomerase pathway", "glucose transmembrane transport", and "regulation of lipid metabolic process". For the upregulated genes in the FoEV treatment vs. Co comparison, we found high enrichment of "translation" and also of functional terms Each sample represents a pool of embryos. Mean-centered expression values (log2 counts per million of sample-mean of log2 CPM of all samples) for the samples of the three embryo treatments. Color scale is from -2 (blue, lower than mean) to 2 (red, higher than mean).

Functional Annotation Revealed Potential Regulatory Pathways Affected by Oviductal EVs
To obtain a more meaningful view of how oEVs can modulate the embryo transcriptome, gene ontology terms (GO) terms and pathways for DEGs from the three different comparisons (FoEV treatment vs. Co; FeEV treatment vs. Co, and FoEV vs. FeEv treatment) were analyzed using functional annotation databases in Metascape and with Ingenuity Pathway Analysis (IPA) tools.
In Figure 3, a Metascape heatmap plot represents the functional enrichment analysis for up-and downregulated genes in embryos treated with frozen or fresh EVs compared to control and between EV treatments, and shows clusters of enriched functional terms across DEGs from different comparisons. The heatmap illustrates that functional terms and pathways related to "embryo morphogenesis" were only enriched for genes downregulated in FoEV treatment vs. Co and FeEV treatment vs. Co sets. Furthermore, genes downregulated in FoEV treatment vs. Co were also enriched for biological processes related to "female gamete generation", "apoptosis", "response to external stimulus", "response to osmotic stress", "telomerase pathway", "glucose transmembrane transport", and "regulation of lipid metabolic process". For the upregulated genes in the FoEV treatment vs. Co comparison, we found high enrichment of "translation" and also of functional terms involved in "ribosome biogenesis", "membrane trafficking", "mitochondrial organization", and "protein methylation". Additionally, FoEV treatment vs. Co upregulated genes were involved in "negative regulation of actin filament polymerization". When both EV treatments were compared, the GO terms found to be enriched were "translation", "ribonucleoprotein complex biogenesis", and "ribosome biogenesis" for upregulated genes, and "membrane trafficking" for downregulated genes in the FoEV treatment vs. FeEV treatment comparison.
Int. J. Mol. Sci. 2020, 21, 1303 5 of 28 involved in "ribosome biogenesis", "membrane trafficking", "mitochondrial organization", and "protein methylation". Additionally, FoEV treatment vs. Co upregulated genes were involved in "negative regulation of actin filament polymerization". When both EV treatments were compared, the GO terms found to be enriched were "translation", "ribonucleoprotein complex biogenesis", and "ribosome biogenesis" for upregulated genes, and "membrane trafficking" for downregulated genes in the FoEV treatment vs. FeEV treatment comparison. More details for the obtained overrepresented functional terms and pathways for each DEG set obtained from embryos under different EV treatments compared to control, with highlighted genes related to interesting biological processes and pathways can be found summarized in Tables 1-6 and in Supplementary Data S2(Tables S1-S6).  More details for the obtained overrepresented functional terms and pathways for each DEG set obtained from embryos under different EV treatments compared to control, with highlighted genes related to interesting biological processes and pathways can be found summarized in Tables 1-6 and in Supplementary Data S2(Tables S1-S6).       GO:0000184 nuclear-transcribed mRNA catabolic process, nonsense-mediated decay RPL23A Table 5. Metascape functional annotation clusters for downregulated genes in embryos treated with FoEVs compared to FeEVs.  The analysis of the DEG lists performed with IPA software (gene expression core analysis) revealed basically similar results, i.e., similar overrepresented functional terms and pathways as were obtained with Metascape. Thus, the analysis was focused on genes related to embryonic development by searching for all enriched related functions and generating a network of the associated genes, including 17 miRNAs potentially targeting the genes in the network and previously identified in bovine oEVs [19] (Figure 4; Supplementary Data S2(Tables S9 and S10)). Overall, the majority of the assigned genes had lower mRNA levels in the embryos supplemented with frozen oEVs. The genes with the highest connectivity were Nanog homeobox (NANOG), MCL1 apoptosis regulator, BCL2 family member (MCL1), Fos proto-oncogene, AP-1 transcription factor subunit (FOS), activating transcription factor 2 (ATF2), apoptotic peptidase activating factor 1 (APAF1), cytochrome c, somatic (CYCS), and thrombospondin 1 (THBS1). The processes with the highest connectivity were "apoptosis of embryonic cell lines", "proliferation of embryonic cells", and "cell viability of embryonic cell lines". The genes targeted by the highest number of miRNAs were APAF1 (five miRNAs), FOS (three miRNAs), THBS1 (three miRNAs), and RING1 and YY1 binding protein (RYBP, three miRNAs).  MicroRNAs revealed by the software to be potential upstream regulators were added to the network (restricted to miRNAs known to be present in bovine oEVs) as well as further connected DEGs. Genes are colored by log2 fold change ratio, blue for downregulated and red for upregulated genes.

Comparative and Integrative Analysis of Embryonic mRNAs Altered under oEV Treatment and oEV-Derived mRNAs.
To unveil the oEV mRNA cargo responsible for the embryonic transcriptomic changes observed in the present study, we compared the mRNAs identified in oEVs in our previous study [19] at the post-ovulatory stage (S1) and the embryonic mRNAs identified under oEV treatments (from the postovulatory stage in FoEV-treated embryos) from the present study. The Venn diagram in Figure 5A shows a comparison of all detectable genes in the embryos (Supplementary Data S1 (Table S1)), the DEGs for the comparison FoEV treatment vs. Co (Supplementary Data S1 (Table S2)), all detectable control DEGs was generated using ingenuity pathway analysis software. MicroRNAs revealed by the software to be potential upstream regulators were added to the network (restricted to miRNAs known to be present in bovine oEVs) as well as further connected DEGs. Genes are colored by log2 fold change ratio, blue for downregulated and red for upregulated genes.

Comparative and Integrative Analysis of Embryonic mRNAs Altered under oEV Treatment and oEV-Derived mRNAs
To unveil the oEV mRNA cargo responsible for the embryonic transcriptomic changes observed in the present study, we compared the mRNAs identified in oEVs in our previous study [19] at the post-ovulatory stage (S1) and the embryonic mRNAs identified under oEV treatments (from the post-ovulatory stage in FoEV-treated embryos) from the present study. The Venn diagram in Figure 5A shows a comparison of all detectable genes in the embryos (Supplementary Data S1 (Table S1)), the DEGs for the comparison FoEV treatment vs. Co (Supplementary Data S1 (Table S2)), all detectable genes in the postovulatory-stage oEVs (Supplementary Data S1 (Table S7)), and the top 500 genes with highest expression in postovulatory-stage oEVs (based on transcripts per million, TPM; Supplementary Data S1 (Table S8)). The comparison of all detectable genes revealed an overlap of 9404 genes. Considering that the transcripts with low concentrations in oEVs might not be very likely to have an effect on the embryo, we selected the top 500 most abundant transcripts that showed a frequency of more than 25 TPM [19], resulting in 453 embryonic mRNAs in total and 14 of the DEGs of the FoEV treatment vs. Co comparison (Supplementary Data S1 (Table S9)) as being in common with the top 500 mRNAs in oEVs.
To further test the first of our hypotheses, that a part of the alterations of the embryonic transcriptome could be directly due to RNAs derived from oEVs that simply increase the number of these transcripts in the embryo due to oEV delivery, we compared upregulated DEGs in frozen-EV-derived embryos compared to control (125 genes, Supplementary Data S1 (Table S5)) to all mRNAs in oEVs and the top 500 most abundant mRNAs in oEVs ( Figure 5B). This comparison revealed 86 upregulated genes in common with all oEV mRNAs (Supplementary Data S1 (Table S9)) and 14 upregulated genes in common with the top 500 most abundant mRNAs in oEVs (Supplementary Data S1 (Table S9)). In contrast, the genes downregulated in frozen-EV-supplemented embryos (96 genes, Supplementary Data S1 (Table S6)) did not overlap with the top 500 most abundant mRNAs in oEVs ( Figure 5C), whereas a list of 87 genes in common with all mRNAs in oEVs was obtained (Supplementary Data S1 (Table S9)).
Furthermore, we used gene set enrichment analysis (GSEA) [21] to identify significantly enriched genes in embryos which might be derived from oEVs by comparison to (1) the top 500 genes with the highest expression in oEVs and (2) non-coding RNAs contained in oEVs (Supplementary Data S1 (Table S10)) [19]. The GSEA plot in Figure 5D shows a substantial enrichment of the top 500 most abundant mRNAs identified in oEVs towards the upregulated genes in oEV-supplemented embryos. The list of the 19 genes with the highest enrichment scores obtained from this analysis is also shown in Figure 5D. Moreover, a strong enrichment was also found for a group of ncRNAs identified in oEVs containing mainly small nucleolar RNAs (snoRNAs) and spliceosomal RNAs, as illustrated in Figure 5E. The complete lists of GSEA results can be found in the Supplementary Data S1(Tables S11 and S12). Comparison of Venn diagram and GSEA results revealed a set of 14 genes common to both analyses as potential mRNAs contained in oEVs contributing to upregulation of embryonic gene expression, as well as a set of small non-coding RNAs, including snoRNAs and spliceosomal RNAs.
Moreover, transcripts delivered by oEVs to the embryo could be translated and regulate gene expression in the embryo. To identify potential factors upregulating embryonic gene expression, the upregulated genes in FoEV-treated embryos were subjected to transcription factor (TF) analysis using ChIP-X Enrichment Analysis Version 3 (ChEA3) [22]. This analysis revealed five factors, HMGN3, JUND, NME2, PIN1, and YBX1, of which the mRNAs were contained in the top 500 most highly abundant mRNAs in oEVs, potentially regulating 53 of the upregulated genes (Supplementary Data S1 (Table S13)). upregulated genes in common with all oEV mRNAs (Supplementary Data S1 (Table S9)) and 14 upregulated genes in common with the top 500 most abundant mRNAs in oEVs (Supplementary Data S1 (Table S9)). In contrast, the genes downregulated in frozen-EV-supplemented embryos (96 genes, Supplementary Data S1 (Table S6)) did not overlap with the top 500 most abundant mRNAs in oEVs ( Figure 5C), whereas a list of 87 genes in common with all mRNAs in oEVs was obtained (Supplementary Data S1 (Table S9)).

Comparative Analysis of Embryonic mRNAs Altered under oEV Treatment and Potential Genes Targeted by oEV-Derived miRNAs
To test the second of our hypotheses, that the alterations of the embryonic transcriptome could also be due to miRNAs contained in oEVs that downregulate mRNAs in the embryo, we used three different approaches.
The first approach was based on identifying potential miRNAs that could target the identified downregulated genes in embryos treated with oEVs and comparing them to identified miRNAs in oEVs. Figure 6A represents the comparison of potential miRNAs targeting identified downregulated genes in embryos using miTarBase and TargetScan databases to the miRNAs identified in oEVs (62 miRNAs, Supplementary Data S1 (Table S14)). Based on the overlapping miRNAs, potential target genes derived from miTarBase and TargetScan datasets were compared to identified downregulated genes in embryos ( Figure 6B, Supplementary Data S1 (Table S14)). This first approach provided a list of 75 predicted target mRNAs of miRNAs contained in oEVs common to genes downregulated in embryos supplemented with FoEVs.  The second approach was based on identifying potential target genes of the miRNAs identified in oEVs and comparing them to the observed downregulated genes in embryos treated with oEVs. The Venn diagram illustrated in Figure 6C shows the results of this comparison, representing potential target genes derived from miTarBase and TargetScan databases and the downregulated genes identified in embryos (Supplementary Data S1 (Table S14)). This second approach provided a list of 63 predicted targets of miRNAs in oEVs and common to downregulated genes in embryos supplemented with FoEVs. A comparison of the lists of target genes derived from the two approaches provided a list of 57 commonly identified genes, of which expression is probably downregulated in frozen-oEV-supplemented embryos due to miRNAs derived from oEVs ( Figure 6D, Supplementary Data S1 (Table S14)).
Since we found genes downregulated in frozen-oEV-supplemented embryos to be mainly associated with functions such as embryonic development, embryo death, embryonic pluripotent cell lines etc. in the IPA analysis (Figure 4), we looked for potential activated upstream regulators and selected significant miRNAs that were then added to the network. Seventeen miRNAs identified in oEVs were connected to potential targets in the network according to the IPA knowledge base. Fifteen of these miRNAs were also contained in the 62 miRNAs potentially targeting identified downregulated genes in embryos using the miTarBase and TargetScan databases ( Figure 6A).

Integrative Analysis of mRNAs and miRNAs Contained in oEVs and Embryonic Transcriptome Alterations Induced by oEVs
Datasets of oEV mRNAs (top 500) and predicted targets of oEV miRNA cargo as well as DEGs in embryos in response to FoEVs were plotted in a circular layout (Figure 7 and Supplementary Data S1 (Table S15)), providing an integrative view of embryonic gene expression data and RNAs derived from oEVs and their potential effects on embryonic gene expression. Figure 7 shows the overlaps among gene lists at the gene level. This meta-enrichment summary analysis illustrated the two hypotheses postulated in the present study: (1) mRNAs derived from oEVs increase concentrations of a proportion of mRNAs in embryos treated with oEVs, which was supported by the fact that a considerable proportion of DEGs that were upregulated in embryos supplemented with frozen oEVs were present in the top 500 mRNAs in oEVs, and no overlap was found with genes downregulated in supplemented embryos; and (2) miRNAs derived from oEVs downregulate mRNAs in embryos treated with oEVs, which was supported by the overlap of the genes downregulated in embryos supplemented with frozen oEVs (DW_FoEVsEmb) and the genes predicted as targets of oEV miRNAs (PG_miRNA_EVs_T1 and PG_miRNA_EVs_T2).
Further functional annotation analysis of this integrative data was represented in a heatmap of enriched terms across input gene lists ( Figure 8A, Supplementary Data S2 (Tables S7 and S8)). Interestingly, functional terms related to "mitochondria organization", "ribosome biogenesis", "ribonucleoprotein complex biogenesis", "cellular response to external stimuli", "SRP-dependent co-translational protein targeting to membrane", "protein methylation", and "vesicle-mediated transport" were highly enriched for the upregulated genes in embryos supplemented with frozen oEVs and the top 500 oEV mRNA datasets, while functional terms involved in "response to inorganic substances", "cellular protein catabolic process", "response to endoplasmic reticulum stress", "regulation of intrinsic apoptosis signaling pathways", and "positive regulation of apoptotic process" were highly enriched in downregulated genes (frozen oEV embryo treatment), predicted target genes derived for miRNAs in oEVs, and the top 500 oEV mRNA datasets. Interestingly, functional terms related to "reproductive structure development" were highly enriched in downregulated genes (frozen oEV embryo treatment), predicted target genes derived for miRNAs in oEVs, and the top 500 oEV mRNA datasets. Figure 8B provides a Metascape network of these enriched terms and the protein-protein interaction enrichment analysis. of a proportion of mRNAs in embryos treated with oEVs, which was supported by the fact that a considerable proportion of DEGs that were upregulated in embryos supplemented with frozen oEVs were present in the top 500 mRNAs in oEVs, and no overlap was found with genes downregulated in supplemented embryos; and (2) miRNAs derived from oEVs downregulate mRNAs in embryos treated with oEVs, which was supported by the overlap of the genes downregulated in embryos supplemented with frozen oEVs (DW_FoEVsEmb) and the genes predicted as targets of oEV miRNAs (PG_miRNA_EVs_T1 and PG_miRNA_EVs_T2).  (Tables S7 and S8)). Interestingly, functional terms related to "mitochondria organization", "ribosome biogenesis", "ribonucleoprotein complex biogenesis", "cellular response to external stimuli", "SRP-dependent cotranslational protein targeting to membrane", "protein methylation", and "vesicle-mediated transport" were highly enriched for the upregulated genes in embryos supplemented with frozen oEVs and the top 500 oEV mRNA datasets, while functional terms involved in "response to inorganic substances", "cellular protein catabolic process", "response to endoplasmic reticulum stress", "regulation of intrinsic apoptosis signaling pathways", and "positive regulation of apoptotic process" were highly enriched in downregulated genes (frozen oEV embryo treatment), predicted target genes derived for miRNAs in oEVs, and the top 500 oEV mRNA datasets. Interestingly, functional terms related to "reproductive structure development" were highly enriched in downregulated genes (frozen oEV embryo treatment), predicted target genes derived for miRNAs in oEVs, and the top 500 oEV mRNA datasets. Figure 8B provides a Metascape network of these enriched terms and the protein-protein interaction enrichment analysis.

Discussion
This study demonstrated that the supplementation of IVC media with oEVs alters the embryonic transcriptome. Moreover, differences in embryonic gene expression were found between embryos supplemented with oEVs that were frozen before adding to IVC (FoEVs) and with fresh oEVs (FeEVs). In the following sections, we discuss these different effects of FoEVs and FeEVs on the embryo transcriptome, which supported the differential functional changes observed in embryonic development when FeEVs and FoEVs were used in a previous study [13]. Moreover, we discuss different modes of action through which the oEV RNA cargo could regulate the embryo transcriptome and embryonic development by providing different modes of integrative analysis of data derived from our studies.

Differential Effect of Fresh and Frozen oEVs on Embryonic Transcriptome
The differential functional effects of fresh and frozen oEVs on in vitro embryonic development, as demonstrated in on our previous study [13], were reflected in the observed changes in the embryonic transcriptome in the present study. It is worth mentioning here that both oEV treatments were obtained from the same oviducts and split into two samples for FoEV and FeEV treatments. Thus, the only difference in the oEVs used for IVC supplementation was the freeze-thaw step for the FoEV samples. Studies on the effect of storage of EVs at −80 • C have shown controversial results. While some studies have indicated that storing EV samples at −80 • C does not alter EV morphology or size [23,24], others have demonstrated that EV integrity can be disrupted by freezing-thawing [25,26]. Boch and colleagues [25] evaluated the impact of one and four freeze-thaw steps on EV damage and showed that no measurable differences in the particle size distribution and concentration measured by nanoparticle tracking analysis (NTA) were observed after the first freeze-thaw step. However, four freeze-thaw cycles induced a slight increase in the particle size distribution and particle number for EVs stored in PBS in contrast to EVs stored in Trehalose, a protein stabilizer and cryoprotectant widely used in the food and pharmaceutical industries. We hypothesized that freezing might induce membrane alterations in FoEVs and faster leaking or release of EV content compared to FeEVs, thus resulting in a better transfer of the cargo to target cells. The study by Teng et al. (2015) [27] supported our hypothesis, showing a significant reduction in the bi-layer membrane of frozen vesicles (−80 • C) when compared to fresh EVs, even when only imperceptible changes of size and concentration were observed in EVs after one step of freezing-thawing [25]. Moreover, Maroto et al. (2016) [26] showed that in exosome preparations stored at −80 • C, proteins appeared in the supernatant fraction, suggesting that distinct protein groups leak from exosomes even at a −80 • C storage temperature. On the other hand, studies on EV cargo have shown that freezing seems to almost completely preserve the EV-associated proteins [28] and does not impair their functionality [23,29], as we observed in our study. It is worth noting that most EV functional studies have been based on frozen EVs for practical reasons, or due to the impossibility of performing functional studies immediately after EV isolation and characterization. However, not many studies have compared the functional effects of frozen and fresh EVs on cells to date. Further studies are required to elucidate the impact of the freezing process on EVs' integrity, cargo, and functional effects, and the use of protective substances such as Trehalose. Nevertheless, the results from our laboratory provide clear evidence of the positive effect of fresh and frozen oEVs on embryonic development [13], and the extent of such effects to inducing modifications in embryonic development.

Oviductal EVs Regulate Early Embryonic Development by Altering the Embryonic Transcriptome
To date, our study represents the first high-throughput analysis of the effects of oEVs on the embryonic transcriptome. Two previous studies using a targeted approach by RT-qPCR showed that in vitro oEV supplementation induced few changes in gene expression related to epigenetic alterations, metabolism, and embryonic development [14,18]. Five genes, PAG1, AQP3, LDLR, DNMT3A, and SNRPN, were found to be upregulated in embryos when oEVs or isthmus-derived oEVs were used in in vitro embryo culture compared to control (with serum without EV depletion), while only two genes, interferon tau (IFNT) and PLAC8, were found to be downregulated in embryos compared to control (in absence of serum). Interestingly, a few genes have been studied in embryos cultured in media different to the common SOF media used in bovine embryo culture [30] (TCM199 versus DMEM), showing that the effect of oEVs on embryonic gene expression of CX43/GJA1, GAPDH, and G6PD9 also depends on the medium used. However, all of the 10 genes altered due to oEV treatments in Lopera-Vasquez et al. [14,18] were found among the DEG sets in our study, which could be explained by different reasons: different oEV sources (in vitro versus in vivo origin), oviducts related to the side of ovulation and stage of the estrous cycle (ipsi-or/and contralateral oviducts; mid-luteal phase versus post-ovulatory), anatomical region of the oviduct (isthmus and ampulla versus complete oviductal fluid), EV isolation method, concentration of oEV vesicles used, the medium used for embryo culture in which the oEVs were diluted and the use of serum (without EV depletion vs. with depletion) as a supplement, and the methodology used for analyzing gene expression (RT-qPCR vs. low-input RNA-seq). The current knowledge about the effect of oEVs on embryonic development and the differences found among studies emphasizes the need for further studies to examine different variables that could affect the embryo oEV treatment (e.g., medium of embryo culture, supplements, duration of co-incubation, EV source, isolation method, and concentration) in order to establish reliable protocols that can be used to optimize in vitro embryo production in different species.
On the other hand, da Silveira et al. [31] analyzed the effect of EVs from follicular fluid during in vitro maturation and in vitro embryo culture on embryonic development and embryonic gene expression, and showed that genes involved in embryonic development (ACSL6, CDH1, REST) or methylation (DNMT3A) were altered due to follicular EV treatment during IVC. Although CDH1, REST, and DNMT3A have been identified in oEVs (not ACSL6), expression of any of these genes was altered in our study. It is worth noting that da Silveira et al. [31] found differences in embryonic gene expression between embryos supplemented with EVs from pre-ovulatory follicles and control for ACSL6, between treatment with EVs from follicles 3-6 mm in diameter and control for CDH1, and between both EV treatments and control for REST, indicating an important effect of the source of the EVs. The data derived from this study indicated the importance of the EV source for follicular or oviductal EVs.
In our study, we found that 25 of the 125 mRNAs upregulated in FoEV embryos compared to controls were annotated as different types of small non-coding RNAs (sncRNAs), such as small nucleolar RNA, C/D (known as SNORDs), small nucleolar RNA, H/ACA (known as SNORA), and six small nuclear RNAs, which play important roles in RNA splicing. To date, with the exception of the studies addressing the regulatory function of miRNAs in embryonic development [32,33], the regulatory function of other classes of sncRNAs in the embryo remains largely unknown. Small nucleolar RNAs form a specific class of small (60-170 nucleotides) non-coding RNAs that is best known for guiding the post-transcriptional modification of other non-coding RNAs such as ribosomal RNAs (rRNAs) and small nuclear RNAs (snRNAs) [34,35]. Additionally, snoRNAs have key regulatory functions in various other cellular processes, and their altered expression has been associated with a wide range of disorders; they have been suggested as useful diagnostic tools and biomarkers in endometrial and lung cancers [36,37]. El Hajj et al. [38] suggested that decreased SNORDs promoter region methylation in ICSI children may modulate cancer susceptibility, based on the identification of SNORDs such as SNORD11 as being differentially methylated in ICSI versus control umbilical cord blood samples, showing a relationship to processes during early embryonic development. Interestingly, several SNORDs have been identified in EVs secreted by embryos during in vitro culture [39]. From a total of 32 snoRNAs identified, two (SNORD110, SNORD81) were exclusively present in media from viable embryos, whereas media from non-viable embryos had three exclusive snoRNAs (SCARNA24, SNORD97, SNORD48), indicating them as potential biomarkers of embryo viability. However, none of these snoRNAs were found in oEVs or found to be altered in the embryonic transcriptome due to oEV treatment. By contrast, the small nucleolar RNA host gene 3 (Snhg3), identified in oEVs, has been shown to be essential for self-renewal and pluripotency maintenance of murine embryonic stem cells (mESCs), and knockdown of Snhg3 disrupted mouse early embryonic development [40]. Furthermore, this study showed that Snhg3 formed a positive feedback network with Nanog and Oct4 and interacted with 126 proteins in mESCs. Interestingly, the bovine SNHG3 gene hosts the snoRNA genes LOC112443630, LOC112443631, and LOC112443632, belonging to the snoRNA SNORA73 family, of which LOC112443631 was found to be upregulated in FoEV-treated embryos. Another very recent study analyzed sncRNA expression in bovine IVF embryos during the maternal-to-zygotic transition (MZT) period and found a marked increase of sncRNAs, including snoRNAs, during the time of embryonic genome activation [41]. Based on the obtained results, the authors suggested a possible regulatory role of snoRNAs during the MZT in mammalian embryogenesis [41]. Although the role of snoRNAs in embryonic development is not very clear, the results of our and other studies call for further investigation.
To conclude this section, we would like to mention that the observed changes in embryonic gene expression as a result of embryo treatment with oEVs were analyzed after 8 days of IVC with oEVs. We hypothesize that greater transcriptome changes in embryos treated with oEVs could be found if embryos were analyzed after a shorter period of IVC, for example after a few hours. Bland et al. [42] observed changes in the T-cell transcriptome after as little as 0.5 h of exosome treatment, while other transcriptome changes were observed after 8 h of treatment. This study suggested that exosome treatment elicits a dynamic transcriptomic signature in cytotoxic T cells that becomes apparent for some clusters of genes at 0.5 h, while others needed a longer treatment period. Therefore, it is likely that a dynamic transcriptome response to oEV treatment might also be observed in embryos at different culture time points, and also depending on the embryo stage.

Regulation of the Embryonic Transcriptome by oEV-Derived mRNAs: mRNAs in oEVs Upregulate mRNAs in oEV-Treated Embryos
In our study, 125 genes were upregulated in embryos upon FoEV treatment during IVC. We hypothesize that this upregulation could be in part due to the incorporation of transcripts delivered by the oEVs to the embryo after uptake. An example of this oEV control of gene expression could be the higher gene expression of SPINT2 found in FoEV-treated embryos compared to controls, and its presence in the top 500 most abundant transcripts identified in oEVs. Moreover, the functional analysis revealed that SPINT2 is associated with "epithelial cell morphogenesis involved in placental branching" and "epithelial cell differentiation involved in embryonic placenta development". Moreover, it has been shown that in mice, Spint2 contributes to the appropriate development of the embryo, as indicated by Spint2 knockout embryos showing clefting of the embryonic ectoderm, neural tube defects, and defective placental branching morphogenesis [43].
Furthermore, once the transcripts have been incorporated into the embryos, they could be translated into proteins that could regulate embryonic gene expression of the embryos, such as transcription factors. Interestingly, for one of the predicted regulators of a proportion of the FoEV-treated embryo upregulated genes, Y-box binding protein 1 (YBX1), the mRNA and the protein were both contained in the oEVs. YBX1 is a multifunctional protein, regulating cellular processes as a TF, involved in regulation of apoptosis, translation, cell proliferation, mRNA splicing, repair, differentiation, and stress response, and is also found in extracellular vesicles [44,45]. Recently, YBX1 has been identified as being involved in small ncRNA and mRNA sorting into exosomes [46,47]. Regarding a potential function in embryonic development, Ybx1 has been shown to play essential roles in maternal mRNA stability and early embryogenesis of zebrafish during the maternal-to-zygotic transition [48]. Additionally, PRDX2 and NEAT1, genes of which the expression was also upregulated in embryos under FoEV treatment and which are targeted by transcription factor JUND, also identified in oEVs, play an important role during embryonic development. Higher expression of PRDX2 has been associated with a greater viability of oocytes and embryos [49,50]. Moreover, a significantly lower expression of PRDX2 was found in the first-trimester villous cytotrophoblasts of patients with recurrent miscarriage compared to healthy controls. The knockdown of PRDX2 inhibited proliferation and increased apoptosis of trophoblast cells [51], while the lncRNA Neat1 has been found to play a key role in corpus luteum formation and the establishment of pregnancy in a subpopulation of mice [52]. On the other hand, PSMD4 (also known as Rpn10) was also upregulated in embryos under FoEV treatment and potentially regulated by PIN1; it has been shown to be essential for mouse development, since Rpn10 knockout resulted in early embryonic lethality [53]. Moreover, PSMD4 protein seems to be involved in sperm-oocyte binding ability in pigs [54].

Regulation of the Embryonic Transcriptome by oEV-Derived miRNAs: miRNAs in oEVs Downregulate mRNAs in the Embryo
Our data analysis predicted that 62 miRNAs present in oEVs could target 57 out of the 96 obtained downregulated genes in embryos under FoEV treatment compared to control. Interestingly, some of these genes are involved in reproductive structure development processes. For example, Sp3 transcription factor (SP3) is involved in "trophectodermal cell differentiation" and "embryonic process involved in female pregnancy", and can be targeted by four miRNAs identified in oEVs: miR-27a-3p, miR-484, miR-1260b, and miR-218-5p. Sp3 transcription factor is ubiquitously expressed in early embryos and, as an Sp1-like transcription factor, is a regulator of embryonic development in vertebrates [55]. Moreover, Sp3-/mutant mice showed growth retardation and died at birth [56]. Another example is the gene NANOG, which is associated with "functional embryonic morphogenesis", "embryonic morphogenesis", and "gastrulation" and can be targeted by seven different miRNAs present in oEVs: miR-34a-5p, miR-34c-5p, miR-34b-3p, miR-335-5p, miR-128-3p, miR-150-5p, and miR-125b-2-3p. These two examples, among others discussed in the next section, revealed how oEV-derived miRNAs can modulate embryonic gene expression by downregulating mRNAs in the embryo.
Among the 96 genes downregulated in embryos treated with FoEVs, we found different molecules involved in toll-like receptor (TLR) pathways (TLR 2-10 pathways), as shown by the Metascape functional analysis. The TLRs are a family of innate immune system receptors which recognize various molecular patterns of microbial pathogens, inducing antimicrobial immune responses [57], and also have an important role at the maternal-fetal interface [58]. It has been proposed that exosomal miRNAs such as let-7, miR-21, and miR-29a, identified in oEVs, can act as an unconventional mode to activate TLR7 in mice and cause neurodegeneration and tumor growth and metastasis [59,60]. In the same line, milk exosomal miRNAs have been shown to decrease LPS-induced TLR4/NF-κB signaling pathway activation, reducing LPS-induced inflammation through the NF-κB pathway and inhibiting LPS-induced apoptosis via the p53 pathway [61]. Altogether, these data suggest that miRNAs present in oEVs could regulate components of TLR pathways in the embryo, modulating the maternal immune system at the maternal-fetal interface.

Unveiling the Potential of oEV RNA Components to Regulate Embryonic Development
To unveil the potential of oEV RNA components to regulate embryonic development, we performed an integrative analysis of mRNA and miRNA cargo identified in oEVs [19] and the embryonic transcriptome alterations induced by oEVs observed in the present study. Meta-enrichment functional analysis of these integrative data highlighted interesting functional terms highly enriched for the upregulated genes in embryos supplemented with frozen oEVs and the top 500 oEV mRNAs dataset, such as "protein methylation" and "part of non-histone protein methylation", which is a prevalent post-translational modification and important regulator of cellular signaling and function [62]. Functional terms enriched for downregulated genes in frozen oEV embryo treatment, predicted target genes of the miRNAs in oEVs, and the top 500 oEV mRNAs dataset were related to "reproductive structure development", among others. The parent functional term "reproductive structure development" involves child terms such as "embryonic process involved in female pregnancy", "in utero embryonic development", "blastocyst development", "stem cell population maintenance", "in utero embryonic development", "embryo development ending in birth or egg hatching", "trophectodermal cell differentiation", "decidualization", "maternal placenta development", and "epithelial cell differentiation involved in embryonic placenta development". Moreover, our data showed that functional terms related to both gamete formation processes such as "oocyte maturation", "polar body extrusion after meiotic divisions", "spindle assembly involved in meiosis", "spermatogenesis", "binding of sperm to zona pellucida", "sperm-egg recognition regulation", and "fertilization" were also enriched for downregulated genes in FoEV-treated embryos, predicted target genes of miRNAs contained in oEVs, and the top 500 oEV mRNAs. Altogether, the biological processes and pathways mentioned above bring up potential functional roles of oEVs in modulating embryonic development and contributing to successful pregnancy.
Therefore, based on the proprietary knowledge database of the IPA software, we focused the pathway analysis on the identification of genes in the list of FoEV treatment vs. Co DEGs associated with functions and processes related to embryonic development. The majority of the assigned genes were downregulated in embryos supplemented with frozen oEVs and were mainly related to apoptosis, proliferation, and viability of embryonic cells and cell lines in the IPA analysis. Many of these downregulated genes were assigned to enriched functional terms related to apoptosis, cell growth, proliferation, and embryonic morphogenesis in the Metascape analysis. One of the highly connected downregulated genes was NANOG, well-known as a pluripotency-sustaining factor in embryonic stem cells [63]. In rabbits, lower NANOG expression has been found in in vivo embryos at the 16 cell, morula, and blastocyst stages compared to the same stages in in vitro produced embryos [64]. Since NANOG expression is restricted to the embryoblast (EB) and is repressed in trophectoderm (TE) cells [65], this lower expression in the total blastocyst could derive from a different ratio of the number of EB and TE cells at a given stage in in-vitro-produced and in vivo embryos. It was shown in a recent study [14] that culture of bovine embryos in the presence of BOEC-derived EVs increased the number of TE cells. Altogether, this could be an indication that supplementation of embryo culture with oEVs leads to a more in-vivo-like gene expression of Nanog. Sirtuin 1 (SIRT1) has been shown to regulate apoptosis and Nanog expression in embryonic stem cells in the mouse [66]. In our study, SIRT1 was downregulated in FoEV-treated embryos, which could have been due to miR-29b-3p, which is found in bovine oEVs [19]. In mouse embryonic stem cells, Sirt1 has been identified as a direct target of miR-29b [67]. In addition to Sirt1, many more genes in the network were connected to the process "apoptosis of embryonic cell lines". Some of the downregulated genes were typical apoptosis-promoting genes, such as APAF1, or genes that may have proapoptotic functions depending on the context, like TRIM32 [68], MAP3K1 [69], and ATF2 [70]. One of the downregulated and highly connected genes, MCL1, is actually described as an anti-apoptotic factor, but is expressed in different isoforms, whereas the short version of the MCL1 protein has a proapoptotic function [71,72]. For miR-10a, targeting APAF1 in the embryo-development-related network, anti-apoptotic effects have been shown after delivery by exosomes to follicular cells [73]. Interestingly, miR-101 has been shown recently to improve the early development of bovine somatic cell nuclear transfer embryos when overexpressed in the donor cells by increasing proliferation and vitality, and improving the early embryonic development [74].
Overall, the results revealed a network of genes and their regulatory miRNAs that suggested a mode of action of transferred oEV cargo inducing changes in embryonic gene expression which lead to a decrease in apoptosis of embryonic cells and improved embryo viability. For a variety of the identified miRNAs, a positive effect on embryonic development has been shown, for example miR-23a-3p with increased expression in outgrowth embryos in the mouse [75], and miR-21-5p with positive effects on the development of murine embryos [76]. In-vitro-produced bovine embryos cultured in the presence or absence of oEVs for 8 days were used for transcriptomic analysis by RNA sequencing in this study. The oEVs were obtained from the oviductal fluid of bovine oviducts collected from cows at a slaughterhouse, as described in detail in our previous studies [13,19]. According to the status of the ovary, the animals were determined to be in the postovulatory stage. Pools of 8-14 embryos cultured under three different treatments during in vitro embryo culture (IVC) were used for RNA isolation: (1) embryo samples supplemented with fresh oEVs (seven replicates, a total of 74 embryos in seven embryo pools); (2) embryo samples supplemented with frozen oEVs (eight replicates, a total of 81 embryos in eight embryo pools); and (3) embryo samples without supplementation (control) (six replicates, a total of 68 embryos in six embryo pools). Total RNA from these 21 embryo pool samples was isolated using the RNeasy Micro kit (QIAGEN AG, Hombrechtikon, Switzerland) according to the manufacturer's instructions. RNA quality and concentration were analyzed using the Agilent Bioanalyzer 2100 (Agilent Technologies (Schweiz) AG, Basel, Switzerland), NanoDrop (Thermo Fisher Scientific (Schweiz) AG, Basel, Switzerland) and Quantus Quantiflour ® RNA system (Promega AG, Dübendorf, Switzerland). Samples with best RNA quality and concentration were selected for the generation of RNA-seq libraries (15 libraries, five replicates/embryo treatments). For all samples, the RNA integrity number (RIN) was between 9 and 10.

Transcriptomic
For library preparation, the Ovation SOLO RNA-Seq System Kit (NuGEN Technologies, Inc. For Europe (Leek, The Netherlands) was used. Library preparation was done following the manufacturer's instructions. In brief, aliquots of 1ng of total RNA were prepared for each embryo sample as starting material for RNA-seq library preparation. First, samples were subjected to DNase treatment and primer annealing, followed by cDNA processing and second-strand synthesis. After end repair, adapters were ligated and the first step of library amplification was performed by qPCR. To remove fragments derived from ribosomal RNAs, NuGEN's insert-dependent adaptor cleavage (InDA-C) technology was applied in the next step. At the same time, strand selection was performed. After this step, the second library amplification and purification were performed for each sample using a universal primer and a set of barcode primers for sample multiplexing. Once RNA-seq libraries were prepared, quantitative and qualitative analyses were performed for each of the libraries using Agilent Bioanalyzer 2100 DNA High Sensitivity assays and Quantus Quantiflour ® ONE dsDNA system (Promega). Sequencing of the libraries was done on an Illumina HiSeq 2500 instrument (Illumina Inc., San Diego, CA, USA) at the Functional Genomics Center Zurich (FGCZ). Pooled barcoded libraries were run on two lanes of a single-end flow cell, generating between 4 and 11 million single-end reads (125 bp) per sample.

RNA-Sequencing Data Analysis
The obtained sequence reads (Fastq files) were analyzed with an established analysis pipeline integrated in a local Galaxy installation [77] at the Animal Physiology group, ETH, Zurich. Processing, quality control, mapping, and quantification of the obtained sequences was performed as previously described [78]. The bovine genome assembly ARS-UCD1.2 (bosTau9) was used along with the corresponding GFF3 annotation file from NCBI. Based on mapping information for the reads (BAM files), a read count table for all annotated bovine genes was generated using QuasR Qcount. This count table was filtered to remove sequences with negligible read counts by using the counts per million (CPM) per sample filtering tool [79]. The mean library size and potential CPM cutoff (Counttable statistics, custom Galaxy tool) were calculated and the cutoff set to 4.21 CPM (corresponding to an average of 20 reads per library) for at least 3 out of 20 libraries. This count table was the basis for the subsequent statistical analysis.
The analysis of differential gene expression was performed using BioConductor package EdgeR [80]. Data normalization was performed based on library size (TMM normalization) [81] and with the GLM robust (estimateGLMRobustDisp) [82] function. For comparison of the experimental groups, the following contrasts were set: Frozen versus Fresh, Frozen vs. Control, and Fresh vs. Control. An adjusted p-value (false discovery rate, FDR) of 0.05% was used as the threshold for significance of differentially expressed genes for the Frozen vs. Control comparison. Because of the much lower number of differentially expressed genes (DEGs) obtained for the other two comparisons, the likelihood ratio (LR) of 10.81 corresponding to FDR 0.05% in the Frozen vs. Fresh comparison was used as a threshold for the other two comparisons to achieve a comparable stringency and sensitivity of the significance analysis [83]. RNA-seq data were deposited in NCBI's Gene Expression Omnibus and are accessible through GEO Series (GSE143596, https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE143596).

Data Mining and Bioinformatics Analysis
Gene symbols and Entrez Gene IDs (bovine and putative human orthologs) were mapped for all transcripts using bioinformatics custom tools integrated in a local Galaxy installation. Custom database tools (NCBI annotation mapper, Mammalian Ortholog and Annotation database, MOADb [84]) were used to assign known or putative human orthologous genes. Human gene identifiers or symbols were used for subsequent functional annotation.
To compare mRNAs altered in the presence of oEVs to oEV-derived mRNAs and potential genes targeted by miRNA-derived oEVs, Jvenn, an integrative tool for comparing lists of genes with Venn diagrams, was used [85]. Furthermore, statistical comparison of genes altered due to oEV treatment and mRNA and ncRNAs contained in oEVs was performed using GSEA software [86]. For GSEA, all identified genes in embryos were ranked based on log2-fold change and FDR (log2(fold change + 2) × −log10(FDR)) [87]. The resulting preranked gene list containing the most significantly upregulated genes at the top of the list and the most significantly downregulated genes at the bottom was compared with mRNAs and ncRNAs derived from the oEV datasets [19]. To obtain information about overrepresented biological functions and pathways for the gene sets altered due to the oEV treatments and for further comprehensive comparison with oEV-derived RNAs, the Metascape tool was used [88]. QIAGEN Ingenuity Pathway Analysis software, Winter 2019 release 2019.4 (QIAGEN Inc., https://www.qiagenbioinformatics.com/products/ingenuitypathway-analysis, [89]) was also used for functional analysis, integration, and to understand embryonic gene expression data. MIENTURNET, an interactive web tool for microRNA target enrichment analysis based on miRTarBase (an up-to-date tool for validated interactions) and TargetScan (an up-to-date tool for sequence-based miRNA target predictions) was used [90]. To identify potential transcription factors regulating embryonic gene expression, ChEA3 was used [22].

Conclusions
This study revealed a broad impact of oEVs on the embryo by providing the first high-throughput analysis of the embryonic transcriptome regulated by oEVs. Our results showed a complex embryonic molecular signature modulated by oEVs, wherein the effects of oEVs are in fact the sum of multiple effects induced by the wide range of RNA components of oEVs (mRNAs, miRNAS, SNORDs, mRNAs encoding transcription factors). By integrating data from these different oEV components and their potential effects on the embryonic transcriptome, we proposed different modes of action of oEVs on the embryo. Our study provides the basis for in-depth functional investigations of the role of specific oEV RNA cargoes (mRNAs and miRNAs) controlling early embryonic development, which could impact embryo-maternal interactions and thus have key implications for reproductive success.
Supplementary Materials: Supplementary materials can be found at http://www.mdpi.com/1422-0067/21/4/1303/s1. Supplementary Data S1. Transcriptomic data of embryos under different oEV treatments versus controls and data derived from integrative analysis with oEVs. Contains: Table S1. All identified transcripts in all embryos examined under different in vitro culture treatments. Table S2. Differentially expressed genes (DEGs) in embryos for the FoEV treatment vs. Co comparison. Table S3. Differentially expressed genes (DEGs) in embryos for the FeEV treatment vs. Co comparison. Table S4. Differentially expressed genes (DEGs) in embryos for the FoEV vs. FeEV treatment comparison. Table S5. Upregulated genes in FoEV-treated vs. Co embryos. Table S6. Downregulated genes in FoEV-treated vs. Co embryos. Table S7. All mRNAs identified in oEVs. Table S8. Top 500 most abundant mRNAs in oEVs. Table S9. Data derived from Venn diagram comparisons in Figure 5. Table S10. Data used for GSEA analysis in Figure 5. Table S11. Data derived from GSEA analysis for the Top 500 mRNAs in oEVs in Figure 5. Table S12. Data derived from GSEA analysis for the ncRNAs in oEVs in Figure 5. Table S13. Transcription factor analysis by ChEA3. Table S14. Data derived from Venn diagram comparisons in Figure 6. Table S15. Data used for Metascape circular comparisons. Supplementary Data S2. Functional analysis (FA) of transcripts up-or downregulated in embryos under different oEV treatments compared to control. Contains: Table S1. FA of transcripts downregulated in FoEV-treated embryos versus control. Table S2. FA of transcripts upregulated in FoEV-treated embryos versus control. Table S3. FA of transcripts downregulated in FeEV-treated embryos versus control. Table S4. FA of transcripts upregulated in FeEV-treated embryos versus control. Table S5. FA of transcripts downregulated in FoEV-treated versus FeEV-treated embryos. Table S6. FA of transcripts upregulated in FoEV-treated versus FeEV-treated embryos. Table S7 and Table S8. FA of transcripts included in Metascape circular analysis. Table S9 and Table S10. Molecules and pathways derived from IPA analysis and interaction pathways. Funding: This research was funded by INRA-CI_PHASE founds and the SNSF grant 31003A_173171 and EU in the framework of H2020-SFS-2015-2 under grant agreement n • 677353-2 IMAGE. This work was also supported by COST-Action EPICONCEPT-FA1201 and COST-Action SLAAM-BM1308.