Small RNA-seq analysis of extracellular vesicles from porcine uterine flushing fluids during peri-implantation

Background: The extracellular vesicles (EVs) of uterine flushing fluids (UFs) mediate intrauterine communication between conceptus and uterus in pigs. The small RNAs of UFs-EVs are widely recognized as important factors that influence embryonic implantation in humans, mice, cattle and sheep. Therefore, it can be hypothesized that small RNAs of UFs-EVs are essential for implantation in pigs. Results: Cup-shaped EVs of porcine UFs were isolated and characterized, and then PKH67 stained extracellular vesicles were taken up by endometrial epithelial cells and porcine trophectoderm cells in a time-dependent manner. The expression of small RNAs in these EVs on days 10 (D10), 13 (D13) and 18 (D18) of pregnancy was comprehensively profiled through sequencing. From nine libraries, a total of 152 known microRNAs (miRNAs), 43 novel miRNAs and 6248 known Piwi-interacting RNAs (piRNAs) and 110 novel piRNAs were identified. Hierarchical cluster analysis and principal component analysis suggested that the expression profiles of small RNAs in D13 group were different from those in D10 and D18 groups. Bioinformatics analysis showed that the miRNAs differentially expressed in the comparisons between D10 and D13, D13 and D18, and D10 and D18 groups were involved in important processes and pathways related to immunization, proliferation, angiogenesis, apoptosis and development, which play important roles in embryonic implantation. Conclusions: Our results reveal that EVs from porcine UFs contain various small RNAs with potentially vital effects on implantation. This research also provides resources for studies of miRNAs and piRNAs in the cross-talk between embryo and endometrium.


Background
Previous studies have indicated that the majority of pre-natal losses occur during the peri-implantation period in pigs [1,2]. During this period, the porcine blastocyst expands from 0.5-1 mm in diameter at hatching to 2-mm spheres on day 10 of pregnancy (D10) and begins to rapidly elongate into tubular and filamentous shapes on day 11 or 12 [3,4]. The conceptus starts to attach to the surface of maternal endometrium at day 13 of pregnancy (D13) and the attachment is usually completed between days 18 (D18) and 24 of pregnancy [5,6]. The synchronous communication between conceptus and uterus exerts strong influence on embryo implantation [7]. Furthermore, the extracellular vesicles (EVs) of uterine flushing fluids (UFs) are capable of mediating embryo-endometrium cross-talk during implantation [8,9]. Therefore, in order to improve the rate of successful implantation in pigs, it is necessary to investigate the mechanism by which intrauterine EVs regulate the cross-talk.
Small non-coding RNAs mainly include three classes: microRNAs (miRNAs), Piwiinteracting RNAs (piRNAs) and small-interfering RNAs (siRNAs) [22]. Among these small RNAs, miRNAs are single-stranded and non-coding RNAs consisting of [18][19][20][21][22][23][24] nucleotides [23]. They can regulate gene expression at the post-transcriptional level and participate in many biological processes, including apoptosis, cell differentiation, cell proliferation, and tumorigenesis [24]. The piRNAs are small noncoding RNAs associated with PIWI family proteins that silence transposable elements in the animal germline [25]. The interference of siRNA on gene expression is considered as a naturally occurring biological strategy for silencing alleles in vertebrates [26]. Previous studies have indicated that small RNAs such as miR-21-5p [27], miR-30d [28], miR-26a [29] and let-7a [30] play vital roles in embryo implantation. Despite their acknowledged importance, the biological function and significance of small RNAs in implantation have not been fully elucidated.
In this study, UFs-EVs of gilts were isolated and characterized, and then the expression of small RNAs in these extracellular vesicles was profiled on days 10, 13 and 18 of pregnancy, aiming to provides a resource for studies of miRNAs and piRNAs in communication between embryo and endometrium.

Characterization of extracellular vesicles from porcine UFs
Extracellular vesicles were obtained from porcine UFs by filtration and ultracentrifugation. Western Blot (WB) analysis revealed that the UFs-EVs fraction was positive for exosome-related protein markers (CD63, CD9 and HSP70) and negative for cytochrome C and calnexin ( Figure 1A). The cup-shaped UFs-EVs were observed by transmission electron microscope (TEM) ( Figure 1B). In addition, separated UFs-EVs were subjected to nanoparticle tracking analysis (NTA), revealing that their sizes ranged from 40 to 160 nm ( Figure 1C). Results from Agilent 2100 Bioanalyzer show that the RNA content and RIN/RQN value of UFs-EVs are 2.63±0.78 6 ng/μL and 2.6±0.2, respectively. RNase treatment result in an absence of miRNA, meaning that EVs protected miRNAs from degradation ( Figure 1D).
Then, we tested whether porcine UFs-EVs can be taken up by endometrial epithelial cells (EECs) and porcine trophectoderm porcine trophectoderm cells (PTr2). These porcine UFs-EVs (20 μg/mL) were labeled with the fluorescent dye PKH67 and added into the culture medium of EECs or PTr2 cells for 12 h or 24 h. PKH67 labeled UFs-EVs were successfully internalized by EECs and PTr2 cells in a time-dependent manner. Relative green fluorescenc e emitted by the UFs-EVs was used to assess the amount of EV uptake. The results suggested that there was significant increase in the relative fluorescence indicating an increased uptake over a period of time in EECs group (

Small RNA sequencing data
Nine small RNA-seq libraries of porcine UFs-EVs were generated from nine gilts on days 10, 13 and 18 of pregnancy. The length distribution of reads showed that the highest peak was at 22 nucleotides (nt), which was consistent with the expected length of miRNAs. It is worth noting that the second peak appeared in the length range of 28-33 nt, which might correspond to piRNAs ( Figure 3). Altogether, we obtained 297,715,036 raw reads and 209,120,304 clean reads (Additional file 1: Table S1). Among the clean reads, known miRNAs (17.67%) and piRNAs (48.69%) were identified, while 0.84% and 1.31% of clean reads were aligned to novel miRNAs and piRNAs, respectively (Additional file 2: Figure S1).
Then, quantitative PCR (qPCR) was performed on 16 chosen differently expressed small RNAs, including five known miRNAs, five novel miRNAs, three known piRNAs and three novel piRNAs, to verify the reliability of the sequencing data. All Pearson correlation coefficients between the qPCR data and sequencing data were greater than 0.7, indicating that the expression patterns of miRNAs and piRNAs determined by qPCR were generally consistent with those obtained by sequencing ( Table 1).

Prediction of miRNA target genes
Piano and miRNADeeps2 were used to predict novel piRNAs and miRNAs, respectively. We normalized the relative expression levels of miRNAs to TPM. After filtering (TPM >28 in at least one library), we identified 152 known miRNAs, 43 novel miRNAs, 6248 known piRNAs and 110 novel piRNAs (Additional file 3: Table   S2). The target genes of the miRNAs are shown in Table S4, which were derived from the intersection of miRnada, RNAhybird and Targetscan predictions (Addition file 4: Table S3).

Clustering analysis
To explore the differences in small RNA expression pattern of the porcine UFs-EVs during different phases of peri-implantation, hierarchical cluster analysis (HCA) was preformed to cluster these miRNAs and their expression profiles in the three stages with 43 novel miRNAs and 152 known miRNAs. The results revealed that the expression profiles of small RNAs in D13 group were different from those in D10 and D18 groups ( Figure 4). The results also showed that these miRNAs were grouped into 3 categories: category A (e.g., ssc-let-7a), category B (e.g., ssc-let-7c, ssc-miR-30d, ssc-miR-21-5p) and category C ( Figure 4). In addition, the expression patterns of the small RNAs at three stages were investigated by principal component analysis (PCA). The results indicated that samples of the same stage were clustered together (Additional file 5: Figure S2).

Functional annotation of target genes of highly expressed known miRNAs
To determine whether the expression of miRNAs in porcine UFs-EVs is related to that of miRNAs in endometrium and embryos during peri-implantation period, the expression of top 20 highly expressed known miRNAs in UFs-EVs was compared with that in embryos at D10, endometrium at D10 and endometrium at D12 [31][32][33]. The results indicated that most of these miRNAs were expressed in both the embryo and endometrium (Additional file 6: Table S4). Importantly, ssc-miR-21-5p was the most highly differently expressed miRNA in porcine embryos, endometrium and UFs-EVs, and ssc-let-7a, ssc-let-7g, ssc-let-7f-5p, ssc-let-7i-5p, ssc-miR-26a and ssc-miR-429 were also highly expressed ( Figure 5). The predicted target genes of these miRNAs were used for functional analysis. Gene Ontology (GO) analysis (ClueGo network of GO terms) indicated that the target genes were mainly involved in "G proteincoupled receptor signaling pathway", "cell differentiation", "cell-substrate adhesion" and "cell surface receptor signaling pathway" (Figure 6). Kyoto Encyclopedia of Genes and Genomes database (KEGG) analysis showed that they were mainly enriched in "Cell cycle", "mTOR signaling pathway", "Insulin signaling pathway", and "Progesterone-mediated oocyte maturation" (Additional file 7: Table   S5).

Identification of differentially expressed small RNAs
We compared the expression of small RNAs on D10, D13 and D18 of pregnancy with each other (D10 vs D13, D13 vs D18 and D10 vs D18). Small RNAs (TPM >28 in at least one library) with P-values <0.05 and foldchange ≥2 as differentially expressed (DE) small RNAs (Additional file 8: Table S6). As a result, 172 DE miRNAs and 5602 DE piRNAs were obtained. It is noteworthy that 55 DE (e.g., ssc-let-7a) miRNAs and 2407 DE piRNAs were commonly found in all three comparisons ( Figure 7A, B).

Time-series expression analysis
The short time-series expression miner (STEM) was used to investigate the pattern of DE miRNA expression during the three stages (Additional file 9: Table S7). The typical expression pattern (P<0.05) was profile 6 ( Figure 8B, Additional file 10: Table S8) and profile 5 ( Figure 8C, Additional file 11: Table S9). Expression of 51 miRNAs (e.g., ssc-let-7a, ssc-let-7f-5p, ssc-let-7i-5p, ssc-let-7g) assigned to profile 6 which first increased from D10 to D13 and then stay stable until D18 ( Figure 8B). In addition, 36 miRNAs were enriched in profile 5 (e.g., ssc-miR-21-5p and ssc-miR-26a), and the expression levels of these miRNAs reached the peak on D13 and significantly decreased thereafter ( Figure 8C). Notably, Profile 6 and profile 5 have the same in that miRNA expression level elevated from D10 to D13.

Functional annotation of predicted target genes of DE miRNAs
To analyze the functions of DE miRNAs in porcine UFs-EVs during peri-implantation, GO and KEGG analyses were performed on the predicted target genes of these miRNAs in three comparisons (D10 vs D13, D13 vs D18 and D10 vs D18).

Discussion
The intrauterine EVs could mediate embryo-endometrium cross talk during embryo implantation [8]. In the present study, porcine UFs-EVs were isolated and identified through NTA, TEM and WB. A comprehensive profiling of the expression of small RNAs in these extracellular vesicles on D10, D13 and D18 was carried out via sequencing. The small RNAs (TPM >28 in at least one library) included a total of 152 known miRNAs, 43 novel miRNAs, 6248 known piRNAs and 110 novel piRNAs, among which 130, 42, 5492 and 110 were differentially expressed during the three stages.
Importantly, Functional annotation of predicted target genes of highly expressed known miRNAs and DE miRNAs have been done to define the potential role of UFs-EVs miRNAs during embryo implantation.

Characterization of UFs-EVs
The UFs-EVs were identified as cup-shaped vesicles with a diameter of 40-160 nm and measurement of exosomal protein markers CD63, CD9 and HSP70 [34]. In order to test whether UFs-EVs can be transferred to EECs or PTr2 cells, we established an in vitro model that cell internalize labelled EVs. In our study, PKH67 labelled UFs-EVs were successfully internalized by EECs and PTr2 cells in a time-dependent manner. We concluded that UFs-EVs could be taken up by both endometrium and embryo trophoblast.

D10 Group Directly Cluster with D18 Group
In the clustering analyses, the PCA results suggested that the expression profiles of small RNAs in porcine UFs-EVs vary among different stages. HCA of miRNAs showed that D10 group was directly clustered with D18 group but not with D13 group. These 11 results imply that the expression profiles of small RNAs in porcine UFs-EVs vary significantly from D10 to D13, possibly because of the constant elongation and subsequent attachment of blastocysts to the surface of endometrium during this period.

Development
Previous studies have shown that highly expressed miRNAs may play vital roles in the regulation of embryo implantation [35,36]. The expression of top 20 most highly expressed known miRNAs in UFs-EVs was compared with that in endometrium at D10 or D12 and embryo at D10 [31][32][33]. Among these miRNAs, ssc-miR-21-5p was the most abundant in all libraries. MiR-21-5p was also found in ovine mammary gland [37] and milk [38] during early pregnancy. MiR-21 (miR-21-5p) in extracellular vesicles is conducive to embryo development and prevention of high embryo death rate in mice [39]. Moreover, ssc-miR-21-5p, ssc-let-7a, ssc-let-7i-5p and ssc-let-7g were highly expressed in the endometrium during per-implantation in Meishan and Yorkshire pigs [36]. A previous study also indicated that let-7a alters miRNA expression to influence the implantation competency of the activated blastocysts by regulating Dicer expression [30]. Other studies indicated that miR-26a influence endometrial receptivity by targeting PTEN in dairy goats [29], and miR-429 suppresses epithelial-mesenchymal transition (EMT) in vivo and in vitro via targeting protocadherin 8 during embryo implantation in mice [40]. The ClueGO analysis of highly expressed known miRNAs indicated that the predicted target genes were enriched in biological processes and molecular functions associated with blastocyst activation (e.g., "G protein-coupled receptor signaling pathway"), embryo development(e.g., "cell differentiation") and endometrial receptivity (e.g., "cell-substrate adhesion" and "cell surface receptor signaling pathway") [41,42].
Notably, a previous study has suggested that G-protein coupled receptor may play a vital role in embryo-maternal crosstalk during the period of maternal recognition and establishment of pregnancy [43]. Both previous research and our results reveal that these most highly abundant miRNAs in porcine UFs-EVs may act on regulating embryo development and endometrial receptivity at the time of implantation.

DE miRNAs Related to Proliferation, Immunization and Angiogenesis
Of 172 DE miRNAs in three comparisons, 18 DE miRNAs (e.g., ssc-miR-30d, ssc-miR-125b, ssc-miR-200b) were observed for D10 vs D13 and for D13 vs D18 comparisons, and 50 DE miRNA (e.g., ssc-miR-451, ssc-miR-181c, ssc-let-7i-5p) were observed for D10 versus D13 and for D10 versus D18 comparisons. In these miRNAs, a previous study indicated that miR-451 and its target Ankrd46 have effect on embryo implantation in mice [44]. It was also shown previously that downregulation of miR-451 in mouse and human oocytes influenced pre-implantation embryogenesis via inhibiting the Wnt signaling pathway [45]. Krawczynski and coworkers [46] found that ssc-miR-125b was present in porcine UFs-EVs and could regulate expression of vital genes for establishment of pregnancy involved in immune system function, endometrial receptivity or embryo/trophoblast development and implantation.
Interestingly, The expression levels of miR-451 and miR-125b in the endometrium of pregnancy woman with high progesterone levels were lower than them with normal progesterone levels [47]. In addition, the evidence has shown that miR-200b was able to induced EMT via targeting ZEB1 and SIP1 during tumour metastasis [48].
The successful pregnancy during peri-implantation is associated with a series of important biological processes [52]. Function annotation analysis can provide insight into of biological processes and molecular functions of target genes. The GO analysis showed that the main processes enriched in the three comparisons were associated with immune regulation (e.g., immune effect process), signal transduction (e.g., Ras protein signal transduction) and proliferation (e.g., cell cycle). In KEGG analyses, "TNF signaling pathway", "apoptosis", "Axon guidance", "Phosphatidylinositol signaling system" and "Insulin signaling pathway" were common signaling pathways in the three comparisons. Among these pathways, axon guidance is a means of directing the axon process of differentiated nerve cells to its target during embryonic development [53]. TNF signaling pathway was associated with the regulation of apoptosis in murine blastocysts in vitro and in vivo [54]. A previous study has suggested that DE mRNAs between patient with repeated 14 implantation failure and women who conceived after one IVF/ICSI attempt during window of implantation participated in Phosphatidylinositol signaling system [55]. In addition, activation of placental insulin signaling pathway contributes to increases in fetal growth and stimulation of placental amino acids in obese mice [56]. In another study, research data indicated that insulin signaling pathway was related to cell proliferation [57]. As to "AGE-RAGE signaling pathway in diabetic complications", which is significantly enriched in D10 vs D13 comparisons. A previous study indicated that AGE / RAGE signaling can mediate the decline of VEGFR-2 protein and mRNA expression, leading to functional angiogenesis disorders in bone marrow-derived endothelial progenitor cells [58]. Collectively, these results suggest that DE miRNAs in porcine UFs-EVs may have effects on proliferation, immunization, angiogenesis and embryo development during embryo implantation.
Although the DE miRNAs and highly expressed known miRNAs were found and their biological functions were predicted by using bioinformatics analysis in this study, further studies are needed to definitely validate their roles in vitro. Besides small RNAs, other biological information, including mRNAs, proteins, long non-coding RNAs and circle RNAs, is also transported by extracellular vesicles of UFs and involved in the cross-talk between conceptus and endometrium. More sequencing analyses will facilitate a more comprehensive understanding of the regulatory mechanism of porcine UFs-EVs. Furthermore, The piRNAs have been found in spent culture media of embryo [59], our results shown that lots of piRNAs were present in porcine UFs-EVs. However, there are currently not enough studies on porcine piRNAs or effective methods that accurately predict their target genes in pigs, so this analysis of porcine UFs piRNAs is limited. 15

Conclusions
In conclusion, porcine UFs-EVs were isolated and characterized, and proved to be taken up by EECs and PTr2 cells in a time-dependent manner. The small RNA expression in these extracellular vesicles was profiled on D10, D13 and D18 of pregnancy. The target gene functional annotation indicated that these DE miRNAs may influence embryonic implantation by regulating pathways or processes related to proliferation, immunization, angiogenesis and embryo development. This research also provides a resource for studies of miRNAs and piRNAs in communication between embryo and endometrium.

Sample preparation
The animal experiment has been described in Wang et al [60]. shower to relax the pigs, then they were stunned with low-voltage electric shock to reduce the pain and were exsanguinated by puncturing carotid artery to death.
Nine purebred Yorkshire gilts with similar weights (120 ± 10 kg), ages (8 months) and genetic background were selected in this study. Gilts were artificially inseminated using extended semen from one boar (the breeding pig farm of Huazhong Agricultural University) at the onset of estrus (day 0) and again 12 h after. Uteri were obtained from pigs slaughtered on days 10, 13  At the same time, these uteri were used for the separation of endometrial epithelial cells.

Transmission electron microscope (TEM)
The morphology of isolated extracellular vesicles was visualized with high-resolution transmission electron microscope (Hitachi HT7700) based on a previous method [61]. In brief, the resuspended extracellular vesicles were placed on carbon-coated copper grid, and then subjected to 2% phosphotungstic acid (catalog number: DZ0035, Leagene) staining for 2 min. The excess liquid was blotted off by filter paper, and grids were allowed to dry overnight.

Nanoparticle tracking analysis (NTA)
The size distribution and concentration of isolated extracellular vesicles were analyzed by NTA with Zeta View (Particle Metrix). Isolated extracellular vesicles were diluted with particle-free PBS at the ratio of 1:200 and added into the chamber. The results of size distribution and vesicle concentrations were assessed with the software Zeta View.

RNA extraction and RNase treatment
The RNAs of EVs were extracted by miRNeasy Serum/plasma Kit (catalog number: 18 217184, Qiagen, German) according to manufacturer's protocol. The quantity and quality of RNA were assessed using a Small RNA kit and a Small RNA assay on an Agilent 2100 Bioanalyzer (Agilent Technologies, USA). Extracted RNA was treated with 10ng/ul of RNase A (Sigma-Aldrich) at room temperature for 30 min. and then treated RNA was assessed as above.

Cell culture
The uterus of non-pregnant gilts was opened longitudinally on a sterile clean bench.
The endometrium was separated and shredded with a sterile scissor. After washing twice with PBS, the tissue pieces were incubated at 37°C for 2.5 h with collagenase I (catalog number: 17100017, Gibco) and shaken vigorously every half hour.
Undigested tissue pieces were removed by screen filtration. Then, the filtrate was centrifuged at 500 g for 10 min to remove supernatant (fraction rich of endometrial stromal cells). The pellet (epithelial-rich fraction) was resuspended twice in PBS and

Small RNA sequencing
Nine small RNA libraries were generated under standard procedures at Beijing Genomics Institute (BGI, China) according to a previous study [64]. Sequencing of the libraries was carried out with the BGISEQ-500 platform (BGI, China; http://www.seq500.com/en/).

Prediction of novel miRNAs, novel piRNAs and miRNA target genes
Before data analysis, clean reads were obtained by eliminating low-quality reads, adaptors, and other contaminants. Clean reads were then mapped to reference genome and other sRNA databases by Bowtie2 [65]. Because each unique small RNA was supposed to be mapped to only one category, we followed a precedence rule of MiRbase > piRBase > snoRNA (human/plant) > Rfam > other sRNA. Piano [66] and miRNADeeps2 [67] were used to predict novel piRNAs and miRNAs, respectively. (RNAhybrid) and "the 2-8nt of 5 'end of small RNA as seed sequences and 3′-UTR region of porcine transcript were used for prediction" (Targetscan).

Differentially expressed (DE) small RNAs
The expression levels of small RNAs were calculated by using Transcript Per Kilobase Million (TPM) as described in a previous report [69]. Small RNAs with expression value >28 TPM in at least in one of the nine libraries were selected for subsequent data analysis [68]. Raw P-value were convert to adjusted P-value using the Benjamin-Hochberg false discovery rate [70]. To improve the accuracy of the identification of DE small RNAs, we defined a miRNA as a DE miRNA when reads number foldchange ≥2 and P-value ≤ 0.05.

Real time quantitative PCR (RT-qPCR)
First-strand cDNAs of small RNAs were synthesized by Mir-X™ miRNA First Strand Synthesis kit (TaKaRa, Dalian, China) according to the manufacturer's instructions.
The qPCR was carried out by using SYBR Green PCR Master Mix (TaKaRa, Dalian, China) and LightCycler 480II Real-Time PCR System. The forward primers of small RNAs were designed based on the mature sequences of small RNAs (Additional file 12: Table S10), and the reverse primer was a universal primer supplied by the above kit. The expression levels of small RNAs were normalized with U6 to obtain the relative expression by comparative CT method.

Hierarchical cluster analysis (HCA)
The expression values of the miRNAs (>28 TPM in at least one library) were normalized by using the Z-score. HCA was conducted using the R package.

22
The expression values of the small RNAs (> 28 TPM in at least one library) were normalized by using the lg (TPM). PCA was performed using the imageGP (http://www.ehbio.com/ImageGP/index.php/Home/Index/PCAplot.html).

Functional annotation of the predicted target genes
Functional analyses of the predicted target genes were conducted on the basis of Kyoto Encyclopedia of Genes and Genomes database (KEGG) and Gene Ontology slim database (GO-slim) by using KOBAS (http://kobas.cbi. pku.edu.cn/) and PANTHER (http://www.pantherdb.org/), respectively. An enriched functional class with a P-value <0.05 was defined as a significantly enriched target gene candidate.
The network of GO terms based on the GO database was constituted by ClueGO, which is a plugin in Cytoscape (http://www.cytoscape.org/).

Time-series expression analysis
The profiles of DE miRNAs over time were visualized by Short Time-series

Statistical analysis
Calculate the significance of relative fluorescence was assessed by two sets of comparisons using the Bonferroni t-test in SAS 9.1, and shown as * p < 0.05 and ** p < 0.01. Correlation tests were performed to calculate the correlation between the results of RNA-seq and qPCR.

Consent for publication
Not applicable.

24
The datasets supporting the conclusions of this article are deposited in the NCBI SRA repository (https://www.ncbi.nlm.nih.gov/sra/), the accession number is PRJNA530644.

Competing interests
We declared this manuscript have no competing interests. Top 20 most highly expressed miRNAs. D10: day 10 of pregnancy; D13: day 13 of pregnancy; Figure 6 ClueGO network of GO terms. Each node represents a GO term. The class of GO term is reflec