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

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.

1. 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 the FecB (BMPR1B) gene increased the ovulation rate and was partially dominant for litter size [2,3]. According to previously published literature, the FecB mutation has an additive effect on enhancement of the ovine ovulation number and litter size, such that one copy of the FecB mutation may increase the ovulation number by 1.5 and the litter size by 1, and two copies by 3 and 1.5, respectively [4,5,6]. However, some Small Tail Han sheep without the FecB mutation (FecB++) also had high fecundity and maintained stable heredity in production. Therefore, we speculated that, besides the FecB mutation effect, other genetic factors may also influence the litter size of Small Tail Han sheep.

Mammalian reproductive processes include ovarian follicle development, ovulation, luteinization, luteolysis, and remodeling of the endometrium [7,8]. Presently, studies of fecundity have mainly concentrated on the ovary [9,10], but the function of the uterus is also important for prolificacy. The uterus is an organ composed of endometrial and myometrial linings [11,12]. According to the steroid-induced reproductive cycle, the endometrium regenerates and differentiates during the follicular phase, and the metabolic requirements of the endometrium change during the luteal phase [13]. Embryos are supported by maternal secretions from the uterus during the preimplantation period. These secretions form a complex environment that contains proteins, amino acids, pyruvate, lactate, carbohydrates and fatty acids [14,15,16,17]. In addition, hormones [18], homeobox proteins [19], morphogens and other factors directly act on the endometrium [20] and regulated endometrial functions promote interactions between the endometrium and embryo [21]. These functions make it an important source of valuable prolificacy genes to improve sheep fecundity. Therefore, an in-depth understanding of the molecular mechanisms of uterine-related functions is important to study the reproductive features of female sheep.

Circular RNAs (circRNAs) are a unique class of non-coding RNAs originally discovered in viruses, which can be generated by direct ligation of the 5’ and 3’ ends of linear RNAs [22,23,24]. To date, RNAseq were widely used in animals and plants research [25,26,27], and thousands of circRNAs have been identified [28]. Although the functions of animal circRNAs are still being elucidated, some reports have shown that circRNAs act as a sponge for miRNAs [29], regulate gene transcription, and regulate mRNA stability [30]. For example, circRNA-9119 reduces the levels of miR-26a by acting as a microRNA sponge, thereby modulating endometrial receptivity in dairy goats [31]. With the breakthroughs of high-throughput deep sequencing technology, some circRNAs have been identified in the sheep pituitary gland [32], muscle [33,34], spleen [35] and other tissues. However, there are very few reported functions of circRNAs in the sheep reproductive system.

To date, our knowledge of the expression status of the sheep uterus is very limited. To investigate the relationship between circRNA and sheep uterine development and prolificacy, we performed total RNA-Seq for ribose depletion to reveal circRNA expression profiles of polytocous and monotocous Small Tail Han sheep. We also investigated the genomic features, length distribution, and other features of circRNA in the uterus. By comparing the circRNA expression profiles from the uteri of polytocous and monotocous Small Tail Han sheep, we identified differentially expressed circRNAs. We then performed Gene Ontology (GO) annotation and KEGG enrichment analysis, which were common analytical methods [36,37,38,39], for host genes of the differentially expressed circRNAs. CircRNA expression patterns were extensively validated by RT-qPCR. Our study provides the circRNA expression profiles of the polytocous and monotocous uteri, which may facilitate better understanding of the roles of circRNAs in the reproductive mechanisms of sheep.

2. Materials and Methods

2.1. Samples

All experiments were performed following the relevant guidelines and regulations set by the Ministry of Agriculture of the People’s Republic of China. Ethical approval was provided by the Animal Ethics Committee of IAS-CAAS (No. IASCAAS-AE-03, 12 December 2016).

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.

2.2. 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).

2.3. 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.

2.4. 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.

2.5. 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].

2.6. 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−∆∆ 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.

Table 1

Details of primer sequences used for qRT-PCR and expected product sizes of circRNAs.

Gene Name	Primer Sequence (5′–3′)	Product Size (bp)
oar_circ_0030419	F: CCCATCCCCGGCACGTCCAR: GTCCCGTCCTCACTGCACTCG	120
oar_circ_0012660	F: AAACCTGGCACGTACGCAGAR: AGGTTGACGGCATACATCAGC	199
oar_circ_0023919	F: CTTCCTTCAGACAGAATGCACR: CGAGCTCTCCAATATTGTCAC	169
oar_circ_0031129	F: GAGTGAATAAAGATGTCATCCGAR: ACCAGAATATTCAGAAATGTGCC	147
oar_circ_0014321	F: GAAACTGCACACAGATAACCAR: ACATGCTGGATCTAAAACCAC	160
oar_circ_0001853	F: AAGTCACTTAAAACGCACCTR: TAAATTACTGTTCTCCGCTTC	180
β-Actin	F:CCAACCGTGAGAAGATGACCR:CCCGAGGCGTACAGGGACAG	97

3. Results

3.1. 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).

Figure 1

General characteristics of circRNAs in the sheep uteri. (A) Percentages of the six types of circRNAs. (B) Exon length distribution of cricRNAs. (C) Number of exons per transcript of sheep circRNA.

3.2. 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).

Figure 2

Analysis of differentially expressed circRNAs. (A) Differentially expressed circRNAs in the follicular phase. (B) Differentially expressed circRNAs in the luteal phase. Red, green, and gray dots in the graph represent transcripts that were significantly up-regulated, down-regulated, or unchanged respectively, between polytocous and monotocous sheep.

3.3. 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).

Figure 3

GO and KEGG analyses of differentially expressed circRNA host genes. (A) GO functional analysis of differentially expressed circRNA host genes in the follicular phase. (B) KEGG analysis of differentially expressed circRNA host genes in the follicular phase. (C) GO functional analysis of differentially expressed circRNA host genes in the luteal phase. (D) KEGG analysis of differentially expressed circRNA host genes in the luteal phase. The longitudinal and horizontal axes represent the enrichment pathways and rich factor (amount of differentially expressed genes enriched in the pathway/amount of all genes in the background gene set) of these pathways, respectively. The spot size represents the number of differentially expressed genes enriched in each pathway, and the color of the spot represents the p-value of each pathway.

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).

3.4. 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).

Figure 4

Regulatory networks of circRNA-miRNA interactions in the follicular phase. Triangular nodes represent miRNA, and circular nodes represent circRNA.

Figure 5

Regulatory networks of circRNA-miRNA interactions in the luteal phase. Triangular nodes represent miRNA, and circular nodes represent circRNA.

3.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.

Figure 6

Validation of circRNA expression by qRT-PCR and Sanger sequencing. (A) qRT-PCR verification of differentially expressed circRNAs. (B) Back-spliced junctions of the six randomly selected circRNAs was confirmed by Sanger sequencing. The black vertical line represents the back-spliced junction sites. Green and red horizontal lines indicate the 5’ and 3’ ends of the circRNA sequence, respectively.

4. 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.

5. 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.