Small RNA-seq analysis of exosomes from porcine uterine luminal fluid during peri-implantation

Background The exosomes of uterine luminal fluid (ULF) mediate intrauterine communication between conceptus and uterus in pigs. The small RNAs of ULF exosomes 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 ULF exosomes are essential for implantation in pigs. Results Cup-shaped exosomes (diameters: 40–160 nm) of porcine ULF were isolated and characterized, and then PKH67 stained exosomes were taken up by endometrial epithelial cells and porcine trophectoderm cells. The expression of small RNAs in these exosomes 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 145 novel Piwi-interacting RNAs (piRNAs) were identified. Among these small RNAs, ssc-miR-21-5p, ssc-let-7a, ssc-let-7i-5p and ssc-let-7g were differentially and highly expressed during the three stages. 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 molecular functions and pathways such as “G-protein coupled receptor activity” and “Insulin signaling pathway”. Conclusions Our results reveal that porcine ULF exosomes 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. and enriched in profile 5 or profile 6 in STEM analysis. These results and previous research reveal that ssc-miR-21-5p, ssc-let-7a, ssc-let-7i-5p and ssc-let-7g may act on embryo implantation.

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 molecular functions and pathways such as "G-protein coupled receptor activity" and "Insulin signaling pathway". Conclusions Our results reveal that porcine ULF exosomes 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-3 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 exosomes of uterine luminal fluid (ULF) are capable of mediating embryo-endometrium cross-talk during implantation [8]. Therefore, in order to improve the rate of successful implantation in pigs, it is necessary to investigate the mechanism by which ULF exosomes regulate the cross-talk.
Small non-coding RNAs mainly include three classes: microRNAs (miRNAs), Piwi-interacting RNAs (piRNAs) and small-interfering RNAs (siRNAs) [20]. Among these small RNAs, miRNAs are single-stranded and non-coding RNAs consisting of 18-24 nucleotides, and regulate gene expression at the post-transcriptional level and participate in many biological processes, including apoptosis, cell differentiation, cell proliferation, and tumorigenesis [21,22]. The piRNAs are small non-coding RNAs associated with PIWI family proteins that silence transposable elements in the animal germline [23]. The interference of siRNA on gene expression is considered as a naturally occurring biological strategy for silencing 4 alleles [24]. Previous studies have indicated that small RNAs such as miR-181 [25], miR-199a [26], miR-200a [27] and SINE/B1 piRNAs [28] 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, ULF exosomes of gilts were first isolated and characterized, and then the expression of small RNAs in these exosomes was profiled on days 10, 13 and 18 of pregnancy, aiming to gain new insights into the regulatory mechanisms of these small RNAs during peri-implantation period and a better understanding of the communication between conceptus and uterus.

Characterization of exosomes from porcine ULF
Exosomes were obtained from porcine ULF by filtration and ultracentrifugation. Western Blot (WB) analysis revealed the presence of exosome-related protein markers (CD9, HSP70 and GAPDH) in these exosomes ( Figure 1A). The cup-shaped ULF exosomes were observed by Transmission electron microscope (TEM) ( Figure 1B). In addition, separated ULF exosomes were subjected to nanoparticle tracking analysis (NTA), revealing that their sizes ranged from 40 to 160 nm ( Figure 1C).
Then, we tested whether porcine ULF exosomes can be taken up by endometrial epithelial cells (EECs) and porcine trophectoderm porcine trophectoderm cells (PTr2). These porcine ULF exosomes were labeled with the fluorescent dye PKH67 and added into the culture medium of EECs or PTr2 cells. After 24 h, green fluorescent dots were observed in both EECs and PTr2 cells, indicating that ULF exosomes were taken up by these cells ( Figure   1D, E).

Small RNA sequencing data
Nine small RNA-seq libraries of porcine ULF exosomes were generated from nine gilts on 5 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 2). Altogether, we obtained 297,715,036 raw reads and 209,120,304 clean reads (Additional file 2: Table S2). Among the clean reads, 17.67% were identified as known miRNAs, while 0.84% and 1.69% were aligned to novel miRNAs and piRNAs, respectively ( Figure 3).
Then, quantitative PCR (qPCR) was performed on 16 randomly chosen small RNAs, including five known miRNAs, five novel miRNAs and six 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 and 145 novel piRNAs (Additional file 3: Table S3). The target genes of the miRNAs are shown in Table   S4, which were derived from the intersection of miRnada and RNAhybird predictions (Addition file 4: Table S4).

Clustering analysis
To explore the differences in small RNA expression pattern of the porcine ULF exosomes 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 6 of small RNAs in D13 group were different from those in D10 and D18 groups ( Figure 4). Importantly, "ssc-let-7d-5p", "ssc-let-7i-3p", "ssc-let-7i-5p", "ssc-let-7e", "ssc-let-7f-5p", "ssc-let-7a", and "ssc-let-7g" were clustered into a class ( 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 S1)

Functional analysis of target genes of highly expressed known miRNAs
To determine whether the expression of miRNAs in porcine ULF exosomes is related to that of miRNAs in endometrium and embryos during peri-implantation period, the expression of top 20 highly expressed known miRNAs in ULF exosomes was compared with that in embryos and endometrium at D10 [29]. The results indicated that most of these miRNAs were expressed in both the embryo and endometrium (Additional file 6: Table S5).
Kyoto Encyclopedia of Genes and Genomes database (KEGG) analysis showed that they were mainly enriched in "Hippo signaling pathway", "Ras signaling pathway", "Insulin signaling pathway", and "Calcium signaling pathway" (Additional file 8: Table S7).

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 with P-values <0.05 and TPM >28 in at least one library were defined as differentially expressed (DE) small RNAs 7 (Additional file 9: Table S8). As a result, 172 DE miRNAs and 144 DE piRNAs were obtained. It is noteworthy that 55 DE miRNAs and 65 DE piRNAs were commonly found in all three comparisons (Figure 6 A, 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 10: Table S9). The most typical expression pattern was profile 6 (Additional file 11: Table S10), which first increased and then stayed stable. In addition, 36 miRNAs were enriched in profile 5 (Additional file 12: Table S11), and the expression levels of these miRNAs reached the peak on D13 and significantly decreased thereafter (Figure 7 A, B, C).

Functional analysis of predicted target genes of DE miRNAs
To analyze the functions of DE miRNAs in porcine ULF exosomes 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). Significantly enrichment GO terms mainly included "plasma membrane", "transferase activity", "regulation of transcription by RNA polymerase II", "G-protein coupled receptor activity" and "enzyme linked receptor protein signaling pathway" in all three comparisons ( Figure 8). The KEGG analysis showed that the significantly enriched pathways in three comparisons were "Metabolic pathways", "Pathways in cancer", "Axon guidance", "Insulin resistance pathway" and "Insulin signaling pathway" (Figure 9).

Discussion
In the present study, porcine ULF exosomes were isolated and identified through NTA, TEM, WB and high-content screening (HCS). A comprehensive profiling of the expression of small RNAs in these exosomes 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 8 novel miRNAs and 145 novel piRNAs, among which 130, 42 and 144 were differentially expressed during the three stages. Remarkably, bioinformatics analysis showed that the target genes of both DE miRNAs and highly expressed miRNAs were enriched in "Insulin signaling pathway" and "G-protein coupled receptor activity".
In the clustering analyses, the PCA results suggested that the expression profiles of small RNAs in porcine ULF exosomes vary among different stages. HCA of miRNAs showed that D10 group was directly clustered with D18 group but not with D13 group. These results imply that the expression profiles of small RNAs in porcine ULF exosomes may 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.
A previous study has shown that highly expressed miRNAs may play vital roles in the regulation of embryo implantation [30]. The expression of top 20 most highly expressed known miRNAs in ULF exosomes was compared with that in endometrium [31] and embryo [29] collected on D10. Among these miRNAs, ssc-miR-21-5p was the most abundant in all libraries. MiR-21 (miR-21-5p) in extracellular vesicles is conducive to embryo development and prevention of high embryo death rate [32]. MiR-21-5p was also found in ovine mammary gland [33] and milk [34] during early pregnancy. Moreover, let-7a, let-7c, let-7g and let-7i are involved the Lin28B/let-7 axis, which is the key regulatory mechanism of the immune inflammatory response during early pregnancy in cattle [35]. Let-7a alters miRNA expression to influence the implantation competency of the activated blastocysts by regulating Dicer expression [36]. 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 [30]. These RNAs were DE miRNAs clustered in a class and enriched in profile 5 or profile 6 in STEM analysis. These results and previous research reveal that ssc-miR-21-5p, ssc-let-7a, ssc-let-7i-5p and ssc-let-7g may act on embryo implantation.

9
The successful pregnancy during peri-implantation is associated with a series of important biological processes [37]. The GO analysis showed that "G-protein coupled receptor activity" was enriched in the three comparisons (D10 vs D13, D13 vs D18 and D10 vs D18) as well as in the group of highly expressed known miRNAs. G-protein coupled receptor may play a vital role in embryo-maternal crosstalk during the period of maternal recognition and establishment of pregnancy [38]. In KEGG analyses, insulin signaling pathway was the common pathway for the target genes of the above groups of miRNAs.
Insulin signaling pathway is associated with the target genes of exosomal miRNAs in EEC culture medium [39]. Activation of placental insulin and mTOR signaling pathway contributes to increases in fetal growth and stimulation of placental amino acids in obese mice [40]. These results suggest that porcine exosomes may have certain effects on embryo development, maternal recognition and 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 exosomes of ULF 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 ULF exosomes. Furthermore, there are currently not enough studies on piRNAs or effective methods that predict their target genes in pigs, so this analysis of porcine ULF piRNAs is limited. On the other hand, the sequencing data of piRNAs is beneficial for further research on porcine piRNAs in the future.

Conclusions
In conclusion, porcine ULF exosomes were isolated and characterized, and proved to be taken up by EECs and PTr2 cells. The small RNA expression in these exosomes was profiled during three stages of implantation. Notably, ssc-miR-21-5p, ssc-let-7a, ssc-let-7i-5p and ssc-let-7g were highly and differentially expressed in porcine ULF exosomes, and thus may play crucial roles in embryo implantation. Bioinformatics analysis showed that DE miRNAs are involved in important molecular functions and pathways such as "G-protein coupled receptor activity" and "Insulin signaling pathway". 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 [41]. All experiments involving animals were carried out in accordance with regulations (No. 5 proclamation of the Standing Committee of Hubei People's Congress) approved by the Standing Committee of the Hubei People's Congress in P. R. China. All animal sample collection procedures were approved by the Ethics Committee of Huazhong Agricultural University. The animals were housed at Jingpin farm of Huazhong Agricultural University (Wuhan, China) and were slaughtered at the slaughterhouse of Huazhong Agricultural University. Before every gilt was slaughtered, a warm shower to relax the pigs, then they were stunned with lowvoltage 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. After estrus, gilts were artificially inseminated twice using extended semen from one boar (the breeding pig farm of Huazhong Agricultural University). Uteri were obtained from pigs slaughtered on days 10, 13 and 18 of pregnancy, and ULF was collected by flushing with 30 mL phosphate buffer 11 saline (PBS; pH 7.4). The presence of filamentous conceptuses in the flush of the uterine horns was used to confirm pregnancy. For characterization of exosomes from porcine ULF, three additional nonpregnant gilts with similar weights, ages and genetic background (see above) were slaughtered and the uteri were collected. Exosomes were obtained by using the above method. At the same time, these uteri were used for the separation of endometrial epithelial cells.

ULF exosome isolation
The ULF was clarified by centrifugation (2000 g for 30 min at 4°C), filtered through a 0.22 μm filter (MILLEX-GP, USA) to remove impurity. Then, these samples were ultracentrifuged at 130,000 g for 2 h at 4°C (Beckman Optima XE-90, SW32 Ti rotor, Beckman Coulter. USA). Exosomes were re-suspended in 200 μL PBS after washing with PBS (130,000 g for 2 h). Exosomes were stored at -80°C until the RNA was isolated.

Western blot (WB)
The exosomal proteins were extracted by DNA/RNA/protein Isolation Kit (catalog number: R6734-02). The sample was denatured by heating, separated by SDS-PAGE and transferred to a PVDF (polyvinylidene fluoride) membrane. Next, the membranes were blocked with 5% skimmed milk powder and separately probed with rabbit anti CD9 (catalog number:

Transmission electron microscope (TEM)
The morphology of isolated exosomes was visualized with high-resolution transmission electron microscope (Hitachi HT7700 ) based on a previous method [42]. In brief, the resuspended exosomes were placed on carbon-coated copper grid, and then subjected to standard uranyl acetate staining and dried overnight.

Nanoparticle tracking analysis (NTA)
The size distribution and concentration of isolated exosomes were analyzed by NTA with Zeta View (Particle Metrix). Isolated exosomes 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.

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 pellets mainly containing endometrial epithelial cells (EECs) were suspended in Dulbecco's modified Eagle's medium/F-12 (DMEM/F12; 1:1) medium (Gibco, NY, USA) supplemented with 10% fetal bovine serum (Gibco) and 1% penicillin-streptomycin (Gibco) and cultured in 37°C and 5% CO 2 incubator. Endometrial stromal cells were further removed by 0.25% trypsin without edetic acid disodium salt (EDTA) after 2-4 d. Isolation and culture of porcine primary EECs referred to pervious research [41,43].

PKH67 staining and high-content screening (HCS)
Isolated exosomes were stained with PKH67 green fluorescent mini kit (catalog number: MINI67-1KT, Sigma-Aldrich) following the manufacturer's instructions. Briefly, exosomes and PBS (negative control) were respectively added to 1 mL Diluent C, and then incubated with 4 ul PKH67 for 4 min. At last, the mixture was closed with 2 mL 0.5% bovine serum albumin (BSA), washed once with PBS, and then ultracentrifuged at 130,000 g for 2 h at 4°C.
After incubation with PKH67-labeled exosomes or negative control for 24 h, EECs and PTr2 cells were imaged by a high-content screening confocal microscope (PerkinElmer, USA).
The number of fluorescent spots was counted with Opera Phenix TM High Content System as described in a previous study [45].

RNA extraction and small RNA sequencing
The exosomal RNAs were extracted by miRNeasy Serum/plasma Kit (catalog number: 217184, Qiagen, German) according to manufacturer's protocol. The quantity and quality of RNA were assessed with the Agilent 2100 Bioanalyzer (Agilent Technologies, USA). Nine small RNA libraries were generated under standard procedures at Beijing Genomics Institute (BGI, China) according to a previous study [46]. 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 14 and other sRNA databases by Bowtie2 [47]. Because each unique small RNA was supposed to be mapped to only one category, we followed a precedence rule of MiRbase > pirnabank > snoRNA (human/plant) > Rfam > other sRNA. Piano [48] and miRNADeeps2 [49] were used to predict novel piRNAs and miRNAs, respectively. Prediction of the target genes of miRNAs was carried out by RNAhybrid and miRnada.

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 [50]. Small RNAs with expression value >28 TPM in at least in one of the nine libraries were selected for subsequent data analysis [51]. The P-values calculated for every small RNA were adjusted to Q-values for multiple testing corrections by two alternative strategies. To improve the accuracy of the identification of DE small RNAs, the miRNAs with expression value >28 TPM in at least one of the libraries and P-values < 0.05 between two groups were selected as DE small RNAs.

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 1: Table S1), 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.

Principal component analysis (PCA)
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 analyses 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.

Time-series expression analysis
The profiles of DE miRNAs over time were visualized by Short Time-series Expression Miner v 1.3.8 (STEM, http://www.cs.cmu.edu/~jernst/stem/).

Statistical analysis
Statistical significance was assessed by two sets of comparisons using the Bonferroni ttest 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.

Availability of data and materials
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.