Identification of Circular RNAs in the Ovaries of Hu Sheep and Local Fat-tailed Sheep during the Luteal Phase

The ovary is a vital gonadal organ which has the functions of synthesis and secretion of sex hormones, endocrine and reproduction. The maturation and ovulation of primary follicles are regulated by hormones and signaling pathways (Johnson, 2015). During the post-ovulation, the antral follicle and surrounding tissues are reorganized to form a gland with an internal analysis function under the action of gonadotropin, called the corpus luteum (Wang et al., 2020a). The function of the corpus luteum is not only regulated by gonadotropins, but also closely related to cell proliferation and transformation, angiogenesis and metabolite transport in the corpus luteum. Guo et al., (2015) demonstrated that ERK1/2 and AMPK metabolic pathways are associated with luteinizing hormone. Moreover, genes, such as StAR, P450 (CYP11A) and 3β-HSD have appeared many times in many studies and they are the hallmark genes of ovarian corpus luteum development and follicular maturation and ovulation (Song et al., 2021).

In recent years, with the development of high-throughput sequencing, novel bioinformatics methods and related experimental validation, several classes of non-coding RNAs, such as circRNA, lncRNA and miRNA, have emerged and appear to play an important role in various biological processes (Chang et al., 2018; He et al., 2019; Chen et al., 2020; Zhou et al., 2020; Mei et al., 2021; Wu et al., 2021). CircRNA is different from linear RNA without 5' cap and 3' end poly-A tail. It is mainly produced by RNA alternative splicing or reverse splicing and has a closed loop structure with valence (Qu et al., 2015). Due to the special structure, compared with long non-coding RNA, circRNA could resist the degradation of RNase R, which was a particular path to identify circRNA (Suzuki and Tsukahara, 2014). Especially in terms of tissue specificity, spatiality and stability, circRNAs are abundant in eukaryotic cells (Chen and Yang, 2015). Jeck and Sharpless (2014) proved that circRNA was widely present in the cytoplasm. Xu et al. (2019) identified that circ-TCP1 act as competing endogenous RNAs (ceRNAs) that bind miR-183 that regulates porcine ovarian in biological processes related to litter size. In the recent studies, circRNAs have become markers for disease diagnosis (Wang et al., 2019; Sheng et al., 2020). However, studies on sheep luteal ovaries were very few.

In this study, high-throughput transcriptomic circRNA sequencing was performed to reveal changes in the expression profile of circRNAs in luteal phase ovaries in two different breeds of sheep. An in-depth analysis to explore the potential biological functions of specific circRNAs on sheep ovary development was also undertaken. This study adds new comprehensive knowledge on circRNA expression in ovarian tissue, complements key pathways of follicular development in luteal phase ovaries and provides candidate circRNAs for further functional studies and biomarker evaluation.

All experimental procedures involving sheep followed the policies and guidelines on the Henan University of Science and Technology Animal Welfare and Use Committee.
 
Experimental animal preparation
 
The experiment was conducted in 2019-9 and 2019-10 at Henan University of Science and Technology pastures. The geographical location belongs to the western Henan Province of China at 34°32'-34°45'N and 112°16'-112°37'E. Hu sheep and fat-tailed sheep were purchased from Luoning Nongben Biotechnology Co., Limited located in Luoning of Luoyang Province of China. Eight ewes with the same feeding and other management conditions, aged 3.5±0.40 years and weighing 50.5±2.0 kg, each ewes with three consecutive births and divided into OAL group (Hu sheep, n=4, litter size per ewe=3) and OZL group (local Yuxi fat-tailed sheep, n=4, litter size per ewe=1). Ovarian tissues were collected on day 5 of the luteal phase, washed with diethyl pyrocarbonate water (DEPC) and placed in lyophilized tubes, numbered, rapidly frozen in liquid nitrogen and stored in a -80°C refrigerator until follow-up experiments.
 
Extraction of RNA, construction of cDNA library and sequencing
 
Total RNA was extracted from the ovarian tissues using TRlzol reagent (Invitrogen). An amount of 5 μg of RNA from each sample was used to prepare the input material and the circRNA library was prepared according to the steps. 1) Pretreatment of 5μg of total RNA, enrichment of circRNA was digested with 3U of RNase R enzyme reagent (Gisai, Guangzhou, China) per μg of RNA. 2) the RNA pretreament using NEBNext UltraTM Directional RNA Library Prep Kit (NEB, USA). 3) Use the BioAnalyzer 2100 instrument for library quality control and quantification. According to Illunima sequencing instructions, the library was denatured into single-stranded DNA molecules, the second strand cDNA synthesis was subsequently performed using DNA Polymerase I and RNase H, captured on Illunima flowcel, amplified in situ into clusters and selected fragments at 150 to 200 bp in length on an Illumian HiSeq sequencer. The sequencing work was performed by NovoHozyme Biotech.

Identification of differentially expressed circRNA
 
CircRNAs were detected and identified using find-circ (Version 1.1) software (Memczak et al., 2013) and CIRI2 (Version 1.2) (Gao et al., 2018). Circos (Version 0.62-1) software was used to map and display in co-expressed networks. The DESeq2 (Version 1.10.1) R software package was performed to analyze the data for differential expression. DESeq2 provides a model based on the negative binomial distribution to determine the regular statistics of differential expression in digital gene expression databases. In brief, we identified differentially expressed genes with |fold change|≥2.0 and P<0.05 between the two groups as differentially expressed circRNAs.
 
RNase R treatment experiments
 
The RNA samples extracted from the ovarian tissue of Hu sheep were divided into two parts, one was treated with RNase R enzyme (RNase R+) and the other was used as a control (RNase R-). For the treated group of RNA (2.5 ng/ μL), 10 units of RNase R enzyme was added while for the control group, similar quantity of buffer was added and the control group added the same amount of buffer, incubated at 37°C for 30 min and then inactivated the enzyme at 70°C for 10 min and then proceeded to reverse transcription reaction directly.
 
qRT-PCR analysis and sanger sequencing
 
To verify the accuracy of the sequencing results, the ovis glyceraldehyde 3 phosphate dehydrogenase enzyme gene (GAPDH) was selected as the internal reference gene and the expression levels of eight randomly selected differentially expressed circRNAs were verified by applying real-time fluorescence quantitative PCR. Primers were designed using Oligo 7 software and the primer sequences of the eight circRNAs and the GAPDH internal reference gene are shown in Table 1. Then after, the RNA samples were extracted and the first strand of cDNA was synthesized was used as a template for real-time fluorescence quantitative PCR to detect the circRNAs expression level. The qPCR reaction mix for circRNA (20 μl) consisted of 10 μl of 2 x ChamQ Universal STBR qPCR Master Mix, 0.4 μl each for upstream and downstream primers, 1 μl of cDNA and 8.2 μl ddH2O. The qPCR reaction conditions were as follows: 95°C for 15 min, followed by 40 cycles of 95°C for 10 s and 60°C for 30 s. Each sample contains three technical replicates. The qPCR products were sent to Shanghai Biotechnology Co., Ltd for sequencing. According to the 2-DDCt method was used to calculate the relative expression (Rao et al., 2013). Dependent sample t-test was performed to calculate relative expression and SPSS13.0 was used to analyze the significant difference.


 
Co-expression network of circRNA-miRNA construction
 
The Regrna 2.0 database (http://regrna2.mbc.nctu.edu.tw/) was performed to predict the regulation of miRNA by circRNA. The miRNA target sites in the exons of the circRNA locus were identified using miRanda (Version 3.3a) software. According to the predicted binding site of miRNA to target gene of circRNA, the circRNA–miRNA interaction network was constructed using Cytoscape (Version 3.6.4) software for visual data analysis.

Identification features and characteristics of circRNAs
 
The RNA sequencing was undertaken to detect the identities and abundances of circRNAs in luteal ovaries obtained from Hu sheep and fat-tailed sheep. As shown in Table 2, after excluding low-quality raw reads, a total of 92,420,022 and 92,846,510 clean reads from two groups (OAL and OZL), which mapping to the sheep reference genome (Oar-v4.0) of OAL and OZL were 80,858,124 (87.49%) and 80,298,201 (86.48%), the GC contents were 54.8% and 56.44%, respectively. The distribution of reads in known types of genes was shown in (Fig 1A). Totally, 7 711 circRNAs were identified together in the two groups, with 2,875 annotated circRNAs could be compared to the reference genome. In addition, 278 novel circRNAs were obtained in this experiment. The statistics of the density of each chromosome-circRNA show that circRNAs were mainly distributed on chromosomes 1 to 9 and X (Fig 1B,1C). Results of the present study revealed that the length of the sheared circRNAs were between 200 to 500 bp, which was consistent with the content of the previous study (Yang et al., 2020). These circular RNAs are mainly composed of three types (Fig 1D), introns (8%~3%), intergenic (3.4%~8.4%) and exons (both 88.6%). The sequencing results in the ovaries of previous studies, such as pig (Wang et al., 2021), goat (La et al., 2019) and sheep (Liu et al., 2021) animals, were similar to the results obtained in the present study.




 
Differentially expressed analysis of circRNAs
 
Using RNA-seq transcriptome data to analysis the transcriptional levels of circRNA in the ovaries of two sheep breeds. Based on the standard (|fold change|≥2.0 and P<0.05), there were we identified 2,875 differentially expressed circRNAs (DE-circRNAs) by DEseq2. A volcano map of the circRNA expression profiles among the OAL and OZL showed that the 1,173 DE-circRNAs were up- and 1,702 DE-circRNAs were down-regulated (Fig 2).

​
 
Enrichment analysis of differentially expressed host gene
 
In the GO analysis, as shown in Fig 3A, the top three biological processes consisted of the metabolic process (GO:0008152), organic substance metabolic process (GO:0071704) and cellular metabolic process (GO: 0044237). In terms of molecular function, OAL and OZL were mainly enriched in catalytic activity (GO:0003824), binding (GO:0005488), protein binding (GO:0005515), organic cyclic compound binding (GO:0097159) and heterocyclic compound binding (GO:1901363). The cellular components identified were significantly enriched in terms of the intracellular part (GO:0044424) and intracellular (GO:0005622). As shown in the Fig 3B, the three GO domains (cell composition, biological process and molecular function) were represented by a separate fundamental body term. All terms in a domain could trace their parent source to a root term. There might be many different paths through the intermediate terms to the ontology root. KEGG analysis showed that the host genes (STAT1 and UBE2L6) of the DE-circRNAs were mainly enriched in the ubiquitin mediated pathways, such as ubiquitin-mediated proteolysis and protein processing in endoplasmic reticulum; hormone metabolism-related pathways like thyroid hormone signaling pathway and propanoate metabolism (Fig 3C). The host gene, such as MAP9 (Microtubule-associated protein 9) and UBE2L6 also enrichment in GO terms. MAP9, located on chromosome 4q32.1, was a gene responsible for spindle assembly and cytokinesis (Zhang et al., 2020). For example, UBE2L6 could adsorb miR-146a-5p and inhibit the apoptosis of Mtb infected macrophages, which can be used as a potential biomarker of TB. The formation and dissipation of the luteal body provide energy reserves for follicular development, follicular recruitment, ovulation and other biological processes of Hu sheep (Wang et al., 2020b). Therefore, the results of GO and KEGG analyses provide a valuable resource for studying high yield in luteal phase ewes.

​
 
Co-expression network of circRNA-miRNA and function prediction
 
The miRNA binding sites of the selected eight differentially expressed circular RNAs were predicted by miRNA software and the relationship network of circRNA-miRNA was obtained. As shown in the Fig 4, a total of 35 nodes (including 27 miRNAs and 8 circRNAs) and 39 edges constitute the circRNA-miRNA interaction network. These differentially expressed circRNAs can regulate the expression of downstream target genes by antagonizing 28 miRNAs and its molecular mechanism needs to be further verified. Importantly, the predicted results indicate that new-circ-0013661 has an FC value of -2.5 in the OAL/OZL group and targets miR-23a and miR-23b, in addition miR-23b can act on its host gene STAT1 (ENSORARG0 0000013903). These predicted circRNAs are the focus of in-depth exploration.

​
 
Verification of candidate circRNAs using RT-qPCR and Sanger sequence
 
To verify the accuracy of the RNA-seq results, we randomly selected eight DE circRNAs and designed back-to-back divergent primers. The RNA incubated with RNase R was reverse transcribed to cDNA and the tissue expression levels of these eight circRNAs were determined by quantitative real-time PCR (qRT-PCR). Also, the results showed similar expression trends to our RNA-seq data (Fig 5A). The dorsal spore node was verified by Sanger sequencing of the products amplified by the PCR procedure (Fig 5B). It indicated that the loop RNA in this experiment had a specific loop structure. And the tissue expression levels of these eight circRNAs were determined by quantitative real-time PCR (qRT-PCR). Meanwhile, the results showed a similar expression trend with our RNA-seq data. The reverse shear point was verified by Sanger sequencing of the products amplified by the PCR procedure. This indicates that the loop RNA in this experiment has a specific loop structure.

​

Hansen et al. (2013) identified circRNA that contained binding sites of miRNA due to the miRNA response element (MRE), which acted as a competitive endogenous RNA (ceRNA) to modulate gene expression. The ceRNA hypothesis (Wu et al., 2021) states that circRNA can act as sponges for miRNAs and restrict the expression of miRNA target genes. In present study, miRanda software was to the identified miRNA target sites in sheep circRNAs. The ceRNA network had been constructed to explore further the role of the small RNAs that regulated in the ovary. Interestingly, the interaction network displays that oar-miR-16b, oar-miR-23a, oar-miR-23b, oar-miR-544-3p, oar-miR-376b-3p and oar-miR-432 are associated with two circRNAs. The miR-16 inhibited BCL-2 and the NF-KBI/MMp-9 signaling pathway, a tumor suppressor gene in glioma and invasiveness (Yang et al., 2014). The previous study detected that miR-544 inhibits cell proliferation, migration and invasion through down-regulating the 3' UTR of PARK7 (Jin et al., 2016). The study showed that miR-544 regulated dairy goat male germline stem cell self-renewal by targeting PLZF (Song et al., 2015).

In this study, the results of transcriptome sequencing demonstrated the specific expression of circRNA in the ovaries of Hu sheep and local fat-tailed sheep. Functional enrichment analysis of significantly different circRNA host genes revealed that STAT1 (novel-circ-0004071) and UBE2L6 (novel-circ-0013661) were predicted to influence the molecular mechanism of ovarian development in sheep through the ubiquitination pathway. These results contribute to an improved understanding of circRNA and the biological mechanisms of ovarian reproduction in sheep.

This research was supported by the China Agriculture Research System of MOF and MAR, CARS-38 and The National Natural Science Foundation of China (31472095).

The authors declared that they have no conflict of interest.