Differential Expression of Circular RNAs in Polytocous and Monotocous Uterus during the Reproductive Cycle of Sheep

Simple Summary The uterus is an important reproductive organ that provides nutrition and place for embryonic development. In this study, we identified circular RNAs by deep sequencing and analyzed their expression in the uteri of polytocous and monotocous sheep (FecB++) during follicular and luteal phases. Gene Ontology (GO) and KEGG enrichment analyses revealed that the source genes of these differential circular RNAs (circRNAs) were mainly enriched in reproductive hormone- and energy metabolism-related pathways. These results provide information on the molecular mechanisms of sheep prolificacy. Abstract CircRNA plays important roles in cell proliferation, differentiation, autophagy and apoptosis during development. However, there are few reports on circRNAs related to livestock reproduction. In this study, we identified circRNAs by deep sequencing and analyzed their expression in the uteri of polytocous and monotocous sheep (FecB++) during follicular and luteal phases. There were 147 and 364 circRNAs with differential expression in the follicular and luteal phases, respectively. GO and KEGG enrichment analysis was performed for the host genes of the circRNAs to predict the functions of differentially expressed circRNAs. These source genes were mainly involved in the estrogen signaling pathway, TGFβ signaling pathway, GnRH signaling pathway, oxytocin signaling pathway, pentose phosphate pathway, and starch and sucrose metabolism related to reproduction and energy metabolism. CircRNA expression patterns were validated by RT-qPCR. Our findings provide a solid foundation for the identification and characterization of key important circRNAs involved in reproduction.


Introduction
Small Tail Han sheep are a Chinese indigenous breed that has attracted a lot of attention because of their high prolificacy and year-round estrus [1]. In Booroola and Small Tail Han ewes, the effect of Based on a TaqMan assay using the FecB mutation probe, pluriparous ewes with the FecB++ genotype (pluriparous ewes without the FecB mutation) were selected from nuclear herds of Small Tail Han sheep in the southwest region of Shandong Province, China. The estrus cycles of all experimental ewes were synchronized. A controlled internal drug release (CIDR) device (300 mg progesterone) was inserted into their vagina for 12 days. Six days after CIDR removal, the ovulation rate was detected by a laparoscopy procedure. The ovulation rate was equal to the total number of corpus lutea on both sides of a ewe's ovary. According to the litter size and ovulation rate, 12 Small Tail Han ewes were selected and equally divided into polytocous ewes (PG, litter size and ovulation number ≥ 2) and monotocous ewes (MG, litter size and ovulation number = 1). There was no significant difference in phenotypic indexes, including average age, average weight, body length, and chest circumference, between the polytocous and monotocous ewes. Estrus synchronization, as described above, was conducted again to collect uterine samples from six polytocous ewes and six monotocous ewes. Three polytocous ewes and three monotocous ewes were euthanized at 45-48 h after CIDR removal (follicular stage) and the other three polytocous ewes and three monotocous ewes were euthanized at 9 days after CIDR removal (luteal stage). Uterine tissues from polytocous group ewes at follicular and luteal stages were named PF (n = 3) and PL (n = 3) respectively, and the uterine tissues from monotocous group ewes at follicular and luteal stages were named MF (n = 3) and ML (n = 3), respectively. The whole fresh uterus from each ewe was collected and then washed with buffer. Uterine tissues were cut to small pieces, pulverized, and mixed with liquid nitrogen, and then stored at −80 • C until use.

RNA Extraction, Library Construction, and RNA-Seq
Total RNA was extracted from the uterine tissues of 12 ewes. Then, measurements of the RNA purity, concentration and integrity were conducted by 1% agarose electrophoresis, a Kaiao K5500 spectrophotometer (Beijing Kaiao Technology Development Co., Ltd., Beijing, China) and RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 System (Agilent Technologies, Santa Clara, CA, USA), respectively.
Library construction was performed using approximately 3 µg total RNA from each sample. Ribosomal RNAs (rRNAs) were removed from the total RNA using a Ribo-Zero™ Gold Kit (Epicentre, Madison, WI, USA). RNA libraries of 12 samples were generated using a NEB Next Ultra Directional RNA LibraryPrep Kit for Illumina (NEB, Ipswich, MA, USA), following the manufacturer's instructions. CircRNAs were randomly fragmented and reverse transcribed into cDNA with random primers. Second-strand cDNA was synthesized using DNA polymerase I, RNase H, dNTPs, and buffer, and the cDNA fragment was purified by QiaQuick PCR, repaired at the end, tagged was a poly(A), and ligated into Illumina sequencing adapters. The samples were amplified by PCR and sequenced in the PE150 sequencing mode (Illumina, San Diego, CA, USA) using the Illumina X-Ten platform. Raw data of the RNA-Seq have been deposited in the SRA public database (Accession number: SRP173986).

Sequence Mapping and Circular RNAs (circRNA) Prediction
Clean reads were mapped to the reference genome (Oar_v3.1) by the HiSAT2 alignment method. CIRI is an efficient and fast tool to identify circRNAs [40]. To ensure the reliability of other circRNAs, the BWA-MEM algorithm was used to perform a sequence splitting comparison, and then the SAM file was scanned to find PCC (paired chiastic clipping) and PEM (paired-end mapping) sites, and GT-AG splicing signals [41]. Finally, the sequence with the junction site was re-aligned with the dynamic programming algorithm. CircRNAs were blasted against the circBase for annotation. The circRNAs that could not be annotated were defined as novel circRNAs. Statistical analysis of the type, chromosome distribution and length distribution of the identified circRNA was conducted.

Differential Expression Analysis of circRNAs
We used SRPBM (spliced reads per billion mapping) to estimate the expression of circular RNA [42]. To identify differentially expressed circRNAs across groups, the DEseq2 package was used [43]. Since the population we selected are the same sheep breed and comes from the same farm, the growth traits are the same. In addition, all individuals are FecB++ genotype, but the litter size has difference. So, we identified circRNAs with a fold change of >1.5 and a p-value of <0.05 between two groups as differentially expressed circRNAs.

Bioinformatics Analysis
KEGG pathway annotations of the host genes for differentially expressed circRNAs were determined with the KEGG database (http://www.genome.jp). GO enrichment analysis of the host genes for differentially expressed circRNAs was performed with the GOseq R package [44]. GO and KEGG pathways with p < 0.05 were considered as significant, p < 0.05 was considered to be significant for functional enrichment. Interactions between circRNAs and miRNAs were predicted based on sequence data from miRanda (3.3a) [45].

Validation of the Expression of circRNAs
We randomly selected six DE-circRNAs (fold change > 1.5, p < 0.05), for subsequent validation. We designed primers encompassing circRNA-specific, back-splice junctions for each candidate circRNA. Details of the primer sequences are summarized in Table 1. qRT-PCR was performed using SYBR Green Real-time PCR Master Mix (TOYOBOCO, LTD, Osaka, Japan) in a Roche LightCycler 480II (Roche, Basel, Sweden), according to the manufacturer's instructions. Real-time PCR was performed at 95 • C for 10 min, followed by 45 cycles of 95 • C for 15 s, 60 • C for 60 s, and 72 • C for 30 s. β-Actin was used as an internal reference to normalize target gene expression. The 2 −∆∆Ct method was used to calculate relative changes in gene expression between control and experimental groups. Finally, PCR products were gel extracted and subjected to Sanger sequencing.

Overview of circRNA Profiles in Polytocous and Monotocous Sheep Uteri
To identify differentially expressed circRNAs during follicular and lutein phases in the uteri of polytocous and monotocous Small Tail Han sheep, we used deep sequencing to analyze circRNAs. We obtained 371,999,010 and 362,973,390 clean reads from polytocous and monotocous sheep during the follicular stage, and 380,416,106 and 373,756,026 clean reads from polytocous and monotocous sheep during the luteal phase, respectively. A total of 32,687 candidate circRNAs were identified by at least one read spanning a head-to-tail splice junction, according to the splice site information of the circular RNA and the relative position of the gene structure. These circRNAs were divided into six types, alter_exon (8.0%), classic (75.3%), intron (0.6%), overlap_exon (11.6%), intergenic (4.0%) and antisense (0.4%) ( Figure 1A). When the circular RNA consisted of one exon, the length of the exon was significantly longer than that of the circular RNA consisting of multiple exons ( Figure 1B). Many circRNAs had only three or four exons ( Figure 1C).

Differential Expression Analysis of circRNAs
Differentially expressed circRNAs in the uterine tissues of polytocous and monotocous sheep were identified by DEseq2, based on a fold change of ≥1.5 and a p-value < 0.05. In the follicular phase, 147 DE-circRNAs were found in PF and MF groups, including 80 up-regulated and 67 downregulated circRNAs ( Figure 2A, Table S1). In the lutein phase, 364 DE-circRNAs were found in PL and ML groups, including 178 up-regulated and 186 down-regulated circRNAs ( Figure 2B, Table S1).

Differential Expression Analysis of circRNAs
Differentially expressed circRNAs in the uterine tissues of polytocous and monotocous sheep were identified by DEseq2, based on a fold change of ≥1.5 and a p-value < 0.05. In the follicular phase, 147 DE-circRNAs were found in PF and MF groups, including 80 up-regulated and 67 down-regulated circRNAs ( Figure 2A, Table S1). In the lutein phase, 364 DE-circRNAs were found in PL and ML groups, including 178 up-regulated and 186 down-regulated circRNAs ( Figure 2B, Table S1).

Differential Expression Analysis of circRNAs
Differentially expressed circRNAs in the uterine tissues of polytocous and monotocous sheep were identified by DEseq2, based on a fold change of ≥1.5 and a p-value < 0.05. In the follicular phase, 147 DE-circRNAs were found in PF and MF groups, including 80 up-regulated and 67 downregulated circRNAs ( Figure 2A, Table S1). In the lutein phase, 364 DE-circRNAs were found in PL and ML groups, including 178 up-regulated and 186 down-regulated circRNAs ( Figure 2B, Table S1).

Gene Ontology (GO) and KEGG Pathway Enrichment Analyses
In the follicular phase, GO terms analyses of the 147 DE-circRNA host genes revealed significantly enriched terms (p < 0.05) in the categories of biological process, molecular function, and cellular components (Table S2). Figure 3A shows the significantly enriched biological processes such as reproduction, regulation of a cell cycle, cell proliferation and growth hormone receptor signaling pathway. A total of 154 terms were enriched in the KEGG pathway analysis in which the estrogen signaling pathway, biotin metabolism, prolactin signaling pathway, TGF-beta signaling pathway and oxytocin signaling pathway were related to reproduction ( Figure 3B, Table S3).
In the follicular phase, GO terms analyses of the 147 DE-circRNA host genes revealed significantly enriched terms (p < 0.05) in the categories of biological process, molecular function, and cellular components (Table S2). Figure 3A shows the significantly enriched biological processes such as reproduction, regulation of a cell cycle, cell proliferation and growth hormone receptor signaling pathway. A total of 154 terms were enriched in the KEGG pathway analysis in which the estrogen signaling pathway, biotin metabolism, prolactin signaling pathway, TGF-beta signaling pathway and oxytocin signaling pathway were related to reproduction ( Figure 3B, Table S3).
In the luteal phase, GO terms analyses of the 364 DE-circRNA host genes revealed significantly enriched terms (p < 0.05) in the categories of biological process, molecular function, and cellular components (Table S2). Figure 3C shows the significantly enriched biological processes such as the developmental process, regulation of cell development, carbohydrate derivative catabolic process, placenta blood vessel development, and glucan catabolic process. A total of 228 terms were enriched in the KEGG pathway analysis in which the GnRH signaling pathway, VEGF signaling pathway, estrogen signaling pathway, oxytocin signaling pathway, pentose phosphate pathway, and starch and sucrose metabolism were related to reproduction and energy metabolism ( Figure 3D, Table S3).  In the luteal phase, GO terms analyses of the 364 DE-circRNA host genes revealed significantly enriched terms (p < 0.05) in the categories of biological process, molecular function, and cellular components (Table S2). Figure 3C shows the significantly enriched biological processes such as the developmental process, regulation of cell development, carbohydrate derivative catabolic process, placenta blood vessel development, and glucan catabolic process. A total of 228 terms were enriched in the KEGG pathway analysis in which the GnRH signaling pathway, VEGF signaling pathway, estrogen signaling pathway, oxytocin signaling pathway, pentose phosphate pathway, and starch and sucrose metabolism were related to reproduction and energy metabolism ( Figure 3D, Table S3).

Target miRNAs of Differentially Expressed circRNAs in PG and MG
A total of 100 differently expressed circRNAs targeting 26 miRNAs were found. Subsequently, we analyzed the interaction between circRNAs and predicted target miRNAs. The resulting circRNA-miRNA association network provided nodes and linkages between circRNAs and their target miRNAs (Table S4). In the follicular phase, a comparison of the circRNA-miRNA networks of PG and MG groups revealed eight miRNAs interacting with 43 differentially expressed circRNAs associated with uterine functions, and oar-miR-370-3p targeted up to 12 cirRNAs (Figure 4). In the luteal phase, a comparison of the circRNA-miRNA networks of PG and MG groups revealed 10 miRNAs interacting with 64 differentially expressed circRNAs associated with uterine functions, and oar-miR-665-3p and oar-miR-370-3p targeted up to 10 and 18 cirRNAs, respectively ( Figure 5).

Target miRNAs of Differentially Expressed circRNAs in PG and MG
A total of 100 differently expressed circRNAs targeting 26 miRNAs were found. Subsequently, we analyzed the interaction between circRNAs and predicted target miRNAs. The resulting circRNA-miRNA association network provided nodes and linkages between circRNAs and their target miRNAs (Table S4). In the follicular phase, a comparison of the circRNA-miRNA networks of PG and MG groups revealed eight miRNAs interacting with 43 differentially expressed circRNAs associated with uterine functions, and oar-miR-370-3p targeted up to 12 cirRNAs (Figure 4). In the luteal phase, a comparison of the circRNA-miRNA networks of PG and MG groups revealed 10 miRNAs interacting with 64 differentially expressed circRNAs associated with uterine functions, and oar-miR-665-3p and oar-miR-370-3p targeted up to 10 and 18 cirRNAs, respectively ( Figure 5).

Target miRNAs of Differentially Expressed circRNAs in PG and MG
A total of 100 differently expressed circRNAs targeting 26 miRNAs were found. Subsequently, we analyzed the interaction between circRNAs and predicted target miRNAs. The resulting circRNA-miRNA association network provided nodes and linkages between circRNAs and their target miRNAs (Table S4). In the follicular phase, a comparison of the circRNA-miRNA networks of PG and MG groups revealed eight miRNAs interacting with 43 differentially expressed circRNAs associated with uterine functions, and oar-miR-370-3p targeted up to 12 cirRNAs (Figure 4). In the luteal phase, a comparison of the circRNA-miRNA networks of PG and MG groups revealed 10 miRNAs interacting with 64 differentially expressed circRNAs associated with uterine functions, and oar-miR-665-3p and oar-miR-370-3p targeted up to 10 and 18 cirRNAs, respectively ( Figure 5).

Validation of circRNA Expression
We selected six differentially expressed circRNAs for quantitative PCR (qPCR) to confirm the sequencing data ( Figure 6A) and verify the specificity by Sanger sequencing (Figure 6B). The expression levels of these circRNAs detected by real-time PCR were consistent with those obtained by sequencing, confirming the reliability of our sequencing results.

Validation of circRNA Expression
We selected six differentially expressed circRNAs for quantitative PCR (qPCR) to confirm the sequencing data ( Figure 6A) and verify the specificity by Sanger sequencing (Figure 6B). The expression levels of these circRNAs detected by real-time PCR were consistent with those obtained by sequencing, confirming the reliability of our sequencing results.

Discussion
In previous studies, circRNAs had been considered as RNA splicing errors [46]. Recent studies have demonstrated that circRNA plays an important role in the biology and developmental processes of sheep [32,35]. However, uterine circRNAs in sheep have received little attention compared with other tissue non-coding RNAs [33,34]. Here, we identified 32,687 circRNAs in the uterus ofpolytocous and monotocous Small Tail Han sheep. Many circRNAs had only three or four exons, which is similar to previous findings in sheep skeletal muscle [34]. These observations in the uterus of Small Tail Han sheep combined with previous findings suggest that the exon number is a universal feature of the uterus.
Research regarding the functional roles of circRNAs in this field is still at its infancy [47]. Presently, there are reports about the expression and potential biological functions of circRNAs in reproductive organs in xenopus tropicalis [48], mice [49,50] and humans [51][52][53]. These early studies laid the foundation for further research. In mice, 2891 circRNAs have been found in oocytes and early embryos, of which about 91% consist of multiple exons. In addition, host genes producing these circRNAs revealed that circRNAs in mouse early embryos may be responsible for cell division and DNA repair [48]. Eleven circRNAs were significantly downregulated, and 46 circRNAs are significantly upregulated in advanced age female granulosa cells [54]. In this study, 147 and 364 circRNAs were differentially expressed in polytocous and monotocous groups in the follicular phase and luteal phases, respectively. Therefore, we predicted that circRNAs might play important roles in the prolificacy of sheep. The apparent regulation of circRNA appears to be a common phenomenon.
The response to prolificacy in animals is a complicated process, involving several genes and metabolic networks. Hormone synthesis and nutrients are important factors [55][56][57]. In our study, four host genes, including AKT3, ITPR1, ITPR2, and PRKACA, are involved in reproduction-related pathways such as GnRH, estrogen, oxytocin and TGF-beta signaling pathways. Studies have shown that GnRH and GnRHR are expressed in human endometrial epithelial cells and stromal cells [58,59]. During the follicular phase, the GnRH/GnRHR system induces a rapid increase in the plasma E2 concentration, leading to the LH preovulatory surge [60]. During the luteal phase, FSH gradually increases with decreasing LH [61]. The GnRH/GnRHR system might play a key role in promoting trophoblast invasion of the maternal endometrium during embryo implantation [62,63]. Estrogen is the major regulator of placental growth and uterine functions in sheep [64]. In the luteal phase, estrogen plays a regulatory role in vascular functions of the uterus and placenta, such as regulation of blood flow, vascular tone, promoting angiogenesis, and vascular remodeling [65][66][67]. One of the most common effects of oxytocin (OT) is the induction of uterine contractions and childbirth [68]. OT stimulates uterine activity of estrus sows and promotes uterine contractions and muscle hyperplasia [69,70]. It is well known that the TGF-β signaling pathway is essential for the normal follicular development and functions of ovarian cells [71][72][73].
In addition, host genes PGM2 and UGP2 expressed in the luteal phase are involved in energy metabolism-related pathways such as the pentose phosphate pathway, galactose metabolism, and starch and sucrose metabolism. As the animal transitions from ovulation to the seventh day of the estrous cycle, the uterine metabolome also changes. Amino acids, benzoic acid, lipid molecules, carbohydrates, purines, pyrimidines, vitamins, and other intermediate and secondary metabolites reach maximum intensities on either day five or seven relative to ovulation [74,75]. In the present study, the results of the KEGG pathway analysis showed that circRNAs were associated with hormone synthesis and energy metabolism, suggesting their potential effects on reproduction by regulating uterine hormone synthesis and energy metabolism, although the specific mechanism remains to be elucidated.
CircRNA can bind to miRNAs and compete with endogenous RNA [76]. Therefore, some circRNAs may inhibit or relieve repression of miRNA for translation [77]. In our study, 100 differently expressed circRNAs targeting 26 miRNAs were found. We deduced that various circRNAs containing common miRNA-binding sites might act as miRNA sponges to regulate the response of the sheep uterus to prolificacy.

Conclusions
We examined the abundance, genomic characteristics, and length distribution of uterine circRNAs, which lay the foundation for the study of uterine circRNAs in Small Tail Han sheep. Meanwhile, we find that several circRNAs take part in the process of development and uterine growth by KEGG pathway analysis. Function prediction indicates that circRNAs play an important role in the reproduction. Our study provides a valuable resource for circRNAs biology, as well as contributes to understanding circRNAs function in growth and development of uterine.

Conflicts of Interest:
All authors declare no conflicts of interest.