Transcriptome Prole of Key CircRNAs and MiRNAs in Oviduct that Affect Sheep Reproduction

Background: CircRNA and miRNA, as classes of non-coding RNA, have been found to play pivotal roles in sheep reproduction. There are many reports of circRNA and miRNA in ovary and uterus, but few in oviduct. In this study, RNA-Seq was performed to analyze the expression prole of circRNA and miRNA in oviduct between follicular phase and luteal phase in FecB BB (MM) and FecB ++ (WW) sheep. Results: The result showed that 15 circRNAs were found differentially expressed between the follicular phase and luteal phase of FecB BB genotype (MF VS. ML), no circRNAs were found differentially expressed between the follicular phase and luteal phase of FecB ++ genotype (WF VS. WL). 9 and 1 differentially expressed miRNAs were identied in MF VS. ML, and WF VS. WL, respectively. GO and KEGG analysis showed that the host genes of circRNAs in MF VS. ML were mainly involved in the Rap1 signaling pathway and PI3K-Akt signaling pathway. GO and KEGG analysis of miRNAs showed that the predicted target genes were mainly involved in the Insulin secretion and Rap1 signaling pathway in MF VS. ML. In WF VS. WL, the result showed that the predicted target genes were mainly involved in the TGF-β signaling pathway, Insulin secretion and Rap1 signaling pathway. In 3 comparison groups, GO and KEGG analysis indicated that a number of host genes and target genes (XPR1, LPAR3, SLC7A11, RAB3A, PLCB3, CREB3L4, LPAR1, LPAR2, FGF18, TACR3, BMP6, SMAD4, SKP1) were found to inuence sheep reproduction. Conclusion: The results of the study will help to elucidate the regulatory mechanisms of circRNAs and miRNAs in And, this study will enrich the sheep circRNA and miRNA databases, provide sheep


Background
Small Tail Han (STH) sheep is widely bred in China for their year-round estrus and polyembryony [1]. The lambing rate of primiparous ewes is about 200%, and in produced ewes is higher than 250% [2]. Studies have demonstrated the vital role of FecB (BMPR1B) gene in sheep fecundity [3][4][5][6]. Our previous research found that all three genotypes of FecB (FecB BB , FecB B+ and FecB ++ ) are distributed in STH sheep, and there is a signi cant correlation between three genotypes of FecB and the litter size of ewes [7]. Thus, we regard STH sheep as suitable animal model to study the molecular mechanism of FecB gene regulation of reproductive traits.
The process of reproduction of sheep is complicated, in which ovarian follicle development, ovulation, luteinization occurred. Most studies of fecundity have focused on ovaries and uterus [8,9], and little is known about oviduct, key tissue in sheep reproduction. In mammals, the oviduct is the rst maternal site that comes into contact with the embryo. This contact occurs to the rst four days after fertilization [10].
In order to ensure an optimal environment during implantation, the oviduct needs to exchange information with the embryo [11]. In this way, a molecular mechanism through which the oviduct interacting with the developing embryo is initiated, and the body can successfully perform epigenetic reprogramming and activate the embryonic genome [12,13]. This embryo-oviduct interaction may cause changes in transcript, which may have effects on offsprings and lead to changes in fecundity. In addition, mammalian oviduct epidermal cells can synthesize and secrete a series of proteins and affect embryonic development through a variety of signaling pathways, which highlight the role of fallopian tubes in sheep's reproductive process. Therefore, further understanding the molecular regulatory mechanisms and signaling pathways of oviduct-related functions is important for studying the reproductive characteristics of sheep. Circular RNAs are a new type of non-coding RNA with regulatory functions [14,15]. It has a closed loop structure and is abundant in eukaryotic transcriptome. Most circRNAs are composed of exon sequences, which are conserved in the transcriptome, and have tissue-speci c and expression speci city at different developmental stages. Studies have found that circRNA can act as a sponge for miRNA, which leads to the inhibition of miRNA activity and increased levels of target genes [16]. For example, circRNA-9119 binds to miR-26a by acting as a microRNA sponge, making it unable to participate in the mediation process, thereby affecting the endometrial receptivity of dairy goats [17]. The further development of high-throughput sequencing has led the research on circRNAs to the molecular level, and there have been more and more reports about sheep circRNAs, which are identi ed in sheep's pituitary [18], hypothalamus [19], uterus [20] and other tissues.
MicroRNAs are a class of non-coding single-stranded RNA molecules with a length of approximately 22 nucleotides, encoded by endogenous genes. They are involved in the regulation of post-transcriptional gene expression in mammals. MicroRNAs were initially pri-miRNAs, after processing, they became microRNAs precursor with a length of about 70 to 90 bases. Pre-miRNAs were processed by enzymes to become mature miRNAs. Some microRNAs are important regulators involved in the ovarian follicular and luteal development [21]. For example, miR-224, miR-378 and miR-383 were found to regulate aromatase expression during follicle development and miR-17-5p and let-7b were vital for development of luteum [22]. MicroRNAs also play a key role in growth and development, such as muscle growth [23] and neurodevelopment [24].
To explore the roles of circRNAs and miRNAs in the oviduct between STH sheep with FecB BB (MM sheep) and STH sheep with FecB ++ (WW sheep), RNA-Seq was performed and comparison of the expression pro les of circRNA and miRNA of MM sheep and WW sheep was conducted. In addition, GO and KEGG enrichment analysis was performed on the host genes of circRNA and predicted target genes of miRNA, which were common analytical methods [25][26][27]. Which may help us better understand the molecular mechanism of circRNAs and miRNAs in sheep reproduction.

Results
Overview of circRNA Pro les of Small Tail Han Sheep Oviduct GO and KEGG Pathway Enrichment Analyses of miRNAs In MF VS. ML group, GO analyses of predicted target genes revealed signi cantly enriched terms in biological process, molecular function, and cellular components (Supplemental Table S6). We found that some GO terms that were related to lipid, enzymatic activity, junction, and binging are signi cantly enriched (Fig. 4a), including lipid transport, protein binding, lipid binding, serine hydrolase activity, adherens junction and so on. In MF VS. ML group, the KEGG pathway analysis showed that ribosome and nucleotide excision repair were the most enriched pathway (Fig. 4b, Supplemental Table S6). We also found some signaling pathways related to reproduction, including Insulin secretion, cAMP signaling pathway, cGMP-PKG signaling pathway, Rap1 signaling pathway. Interestingly, some related pathway involved protein synthesis and transport, such as biosynthesis of amino acids, ABC transporters and ribosome.
In WF VS. WL group, the GO analyses showed that some term were related to lipid, enzymatic activity, junction, and binging ( Fig. 5a, Supplemental Table S6). The KEGG pathway analysis showed that ribosome and nucleotide excision repair was the most enriched pathway (Fig. 5b). Also, we found some signaling pathways related to reproduction, including TGF-β signaling pathway, Insulin secretion, Protein processing in endoplasmic reticulum, cGMP-PKG signaling pathway, Rap1 signaling pathway. Similarly, biosynthesis of amino acids, ABC transporters and ribosome which related to protein synthesis and transport were included.

Validation of circRNA Expression
RT-qPCR was conducted to con rm the sequencing data of circRNAs. Our results indicated that the four selected circRNAs revealed expression trends similar to the sequencing data, suggesting the reliability of our sequencing results (Fig. 6).

Validation of miRNA Expression
RT-qPCR was conducted to con rm the sequencing data of miRNAs. Our results indicated that the ve selected miRNAs revealed expression trends similar to the sequencing data, suggesting the reliability of our sequencing results (Fig. 7).

Discussion
So far, studies have shown that circRNA and miRNA have regulatory potency, especially in sheep's growth and reproduction [18][19][20][21][22][23]. However, little research has been done on circRNA and miRNA in sheep oviduct. Here, we performed RNA-Seq to analyze the circRNA and miRNA of sheep oviduct tissue between follicular phase and luteal phase in MM sheep and WW sheep, and host genes of DE circRNAs and predicted target genes of miRNAs associated with fecundity were identi ed. We also analyzed the distribution of circRNAs in the genome regions and length distribution of circRNAs. Most of the circRNAs genome composition in sheep uterus was intron [28], while the majority of oviduct was composed of exons. The circRNAs with a length of less than 20,000 nt account for the majority, which was consistent with circRNAs of sheep mammary gland. The length distribution of miRNAs in pituitary and ovary of sheep was similar to our results [29,30], with 22nt miRNAs accounting for the majority. Thus, circRNAs may be tissue-speci c and miRNAs are conservative in different tissues.
We identi ed 15 DE circRNAs in MF VS. ML, but no DE circRNAs were found in WF VS.WL. This result indicated a difference in the expression pro les of MM and WW sheep. Among the DE circRNAs, the top two with the highest expression levels were novel_circ_0009938 and novel_circ_0007904, whose host gene were PAWR and SMC6. Studies have found pro-apoptotic WT1 regulator (PAWR) conducted cell apoptosis, which inhibited the growth of prostate cancer cells [31]. Also, PAWR regulates apoptosis in rat follicles in the ovary but suppressed by FSH by activating PKCζ-dependent anti-apoptotic pathway [32].
But, the expression of PAWR was observed up-regulated in granulosa cells, indicating the increased susceptibility of GCs to undergo apoptosis [33]. In this study, the expression level of novel_circ_0009938 was low at follicular phase and increased at luteal phase, which may be explained by these studies, implying the ovary and oviduct coordinate with each other and stay in sync during the estrous cycle [34,35]. Studies have shown structural maintenance of chromosomes 6 (SMC6) is essential for DNA repair and maintenance of genomic integrity [36]. And, SMC6 plays a key role in spermatogenesis and oocyte meiosis [37,38], which indicated that SMC6 may maintain the genomic integrity of sperm and embryo to ensure fertility.
GO analysis showed that most of the host genes of circRNAs were related to response to stimulus, positive/negative regulation, binding and protein activity. The top enriched GO terms in MF VS.ML (at BP level), whose host gene were both XPR1, were response to external stimulus, MAPK cascade, regulation of response to food. Xenotropic and polytropic retrovirus receptor 1 (XPR1) is a gene encoding cellular inorganic phosphate export protein, and its mutation can cause primary familial brain calci cation. The normal development of the fetus is inseparable from phosphorus. This nutrient is mainly transported from the maternal blood to the fetus via the placenta. Xu et al. [39] found that XPR1 was expressed at a high level in the murine placenta, but the placenta of the murine that knocks out this gene was severely calci ed. Another top enriched host gene in GO analysis (at MF level) is SLC7A11. Soluble carrier family 7 member 11 (SLC7A11) gene is a target of p53-mediated transcriptional repression, and p53 can inhibit the uptake of cystine by repressing the expression of SLC7A11. Studies in mutant mice have found that p53 plays an important role in embryonic development [40]. In addition, SLC7A11 exists in the sperm of stallions and regulates the oxidation-reduction status of sperm by exchanging extracellular cystine (Cyss) for intracellular glutamate [41]. Moreover, host gene lysophosphatidic acid receptor 3 (LPAR3) is involved in Rap1 signaling pathway, PI3K-Akt signaling pathway and neuroactive ligand-preceptor interaction. Rap1 combined with GTP activates the PI3K-Ark signaling pathway, and the PI3K-Ark signaling pathway is widely involved in various important processes of mammalian ovarian development [42], and is related to the survival and activation of primitive follicles [43], hormone secretion, and so on. In addition, neuroactive ligand-preceptor interaction is related to the effect of GnRH and GnRHR. LPAR3 was found to be expressed in mouse oviduct, placenta, and uterus, and its essential role in the female reproductive system was reported [44]. Studies have found that progesterone is likely to have a direct effect on LPAR3, and progesterone treatment can increase the expression of LPAR3 mRNA in endometrium [45]. Besides, dynamic changes that occur in the organization of luminal and glandular epithelia in endometrium during the estrous cycle are necessary to modulate the appropriate environment for the developing embryo and to allow implantation of the conceptus [46]. The differentiation of oviductal epithelial cells is also affected by progesterone. Given the key role of progesterone, we suppose that LPAR3 may play crucial roles in sheep reproduction. Therefore, three novel circRNAs-novel_circ_0012086, novel_circ_0001794 and novel_circ_0014274, whose host genes are XPR1, LPAR3 and SLC7A11, respectively, may have key roles in reproduction. However, it remains to be validated.
GO analysis of predicted target gene of miRNAs showed that most terms were related to lipid, enzymatic activity, junction, and binging in MF VS.ML, including lipid transport, protein binding, lipid binding, serine hydrolase activity, adherenss junction and so on. In addition, KEGG enrichment analysis implied that the target genes were mainly involved in the Insulin secretion, cAMP signaling pathway, cGMP-PKG signaling pathway, Rap1 signaling pathway and Calcium signaling pathway. And the DE miRNAs may affect sheep reproduction by modulating target genes associated with the above signaling pathways and biological processes. In MM sheep, these pathways were enriched for miR-370-3p and its target genes (RAB3A, PLCB3, CREB3L4, LPAR1, LPAR2, FGF18, TACR3). RAB3A, with evolutionarily conserved proteins, is essential regulator of membrane tra cking [47]. Researches have found that RAB3A is indispensable in human sperm acrosome reactions and may regulate the quality of oocytes [48,49]. Phospholipase C beta 3 (PLCB3) was considered to be a candidate gene for litter size in Finnsheep, and playing key roles in the development of folliculogenesis and LH signaling [50]. PLCB3 was also differentially expressed in endometrial of heifers with high and low fertility [51]. CAMP responsive element binding protein (CREB) was found to be expressed in follicular granulosa cells [52]. CREB protein concentration increases during sexual maturation and ovarian follicular development. CREB3L4 was detected in different stages of embryogenesis [53]. Same as LPAR3, LPAR1 and LPAR2 play crucial roles in female reproductive system. Studies have found that LPA medium can improve the survival and development potential of follicles, and can stimulate the cell function and E2 synthesis of mouse ovarian tissue [54]. In addition, oviduct is an important part where gamete transport and fertilization happened. LPA was found to be involved in gamete transport, fertilization, and cell signal transmission between oviduct tissue and cumulus oocyte complex [55,56]. And LPAR2 was found to be abundantly expressed in the oviduct of cattle, suggesting that the oviduct is an important target of LPA [57]. Fibroblast growth factor 18 (FGF18) inhibits the secretion of estradiol and progesterone, and is a candidate factor that regulated the steroidogenesis during ovarian development [58]. Moreover, FGF18 is likely to cause granulosa cell apoptosis, thereby affecting follicular atresia [59,60]. Tachykinin receptor 3 (TACR3) plays a key role in regulating gonadotropin secretion and sex hormone feedback regulation of the reproductive axis [61]. Tacr3 may also be related to the regulation of granulosa cell function and changes in ovarian function [62]. Besides, the expression of TACR3/TAC3 can promote the secretion of GnRH [63], which may affect sheep reproduction. These genes are likely to participate in the reproductive process of MM sheep, but speci c molecular mechanism needs further veri cation. GO analysis of predicted target gene of miRNAs in WF VS. WL showed that most terms were related to lipid, enzymatic activity, junction, and binging, which is similar to MM sheep. KEGG enrichment analysis indicated that the target genes were mainly involved in the TGF-β signaling pathway, Insulin secretion, Protein processing in endoplasmic reticulum, cGMP-PKG signaling pathway, Rap1 signaling pathway. In WW sheep, in addition to target genes above (RAB3A, PLCB3, CREB3L4, LPAR1, LPAR2, FGF18, TACR3), miR-370-3p and its target genes (BMP6, SMAD4, SKP1) were enriched in TGF-β signaling pathway as well. Bone morphogenetic proteins 6 (BMP6) is a member of the transforming growth factor-β (TGF-β) superfamily and was found to be highly expressed in mammalian oocytes and granulosa cells [64,65]. Studies have found that BMP6 is involved in primary/secondary follicle transition, dominant follicle selection, ovarian steroid production, follicular atresia, prevention of luteinization, and luteolysis [66][67][68].
Besides, mice genetically missing BMP6 was characterized by reduced ovulation rate, impaired oocyte quality and impaired embryo implantation, resulting in reduced litter size [69]. Seven TGF-1 receptor subtypes and ve type 2 receptor subtypes associated with signal transduction ligand of the TGF-β superfamily have been found in mammals [70]. BMP ligand activates cellular activity by binding to two types of receptors and phosphorylates the responsive SMAD proteins SMAD1/5/8 (R-SMADs) [71]. Then, R-SMADs bind to SMAD, SMAD4, translocating into the nucleus, to regulate the transcription of target genes by combining with other transcription factors [72]. Mothers against decapentaplegic homolog 4 (SMAD4) is a key signal transduction molecule in the TGFβ/SMAD signaling pathway, which plays an important role in the development of mammalian follicles and the proliferation and differentiation of granulosa cells. [73]. Studies have found that speci cally knocking out the SMAD4 gene in ovaries led to premature failure of mouse follicles, premature luteinization of granulosa cells, and decreased fertility [74]. In addition, mice died in the embryonic stage after knocking out SMAD4 [75]. S-phase kinase association protein 1 (SKP1), as a downstream regulator of the TGF-β/SMAD signaling pathway, regulates follicle formation and ovulation in mammals [76]. SKP1 is a key skeleton protein in Skp1-Cull-Fbox protein (SCF), which mediates the ubiquitination and degradation of different cyclins [77], thereby promoting the cell cycle [78]. SCF has also been found to be crucial for oocyte division and maturation [79], as well as the process of fertilization and implantation [80]. Here, the results imply that these target genes are probably related to sheep reproduction process, but the molecular mechanism through which they affect litter size is still unclear. Further experiments are needed to verify these target genes.

Conclusions
In this study, we established the rst circRNA and miRNA expression pro le in sheep oviduct. We also identi ed several key circRNAs and miRNAs (novel_circ_0012086, novel_circ_0001794, novel_circ_0014274 and miR-370-3p). Furthermore, GO and KEGG analysis indicated that a number of host genes and target genes (XPR1, LPAR3, SLC7A11, RAB3A, PLCB3, CREB3L4, LPAR1, LPAR2, FGF18, TACR3, BMP6, SMAD4, SKP1) were found to in uence sheep reproduction. The results of the study will help to elucidate the regulatory mechanisms of circRNAs and miRNAs in sheep reproduction. And, this study will enrich the sheep circRNA and miRNA databases, which provide a basis for further research on sheep reproduction.

Animal Processing and Tissues Acquirement
All animals in this experiment were approved by the Science Research Department (in charge of animal welfare issues) of IAS-CAAS and ethical approval was given by the Animal Ethics Committee of the IAS-CAAS (No. IAS 2019-49). All ewes were from Yuncheng Breeding Sheep Farm (Yuncheng County, China), where they all obtained similar feeding conditions.
First, the TaqMan MGB probe with FecB mutation [81] was used to identify the genotype of Small Tail Han Sheep. And, six MM sheep and six WW sheep were selected according to average age, average weight, body length and chest circumference.
Then, all of the selected STH sheep were treated with CIDR (controlled internal drug releasing; Zoetis Australia Pty., Ltd., NSW, Australia; progesterone 300 mg) for 12 days. Three MM sheep and three WW sheep were euthanized (Intravenous pentobarbital (100 mg/kg)) within 45-48 h of CIDR removal (follicular phase), and the remaining three MM sheep and three WW sheep were euthanized (Intravenous pentobarbital (100 mg/kg)) seven days after CIDR removal (luteal phase), and oviduct tissues were collected. The oviduct tissues of MM sheep at follicular phase and luteal phase were named MF (n = 3) and ML (n = 3), and the oviduct tissues of the WW sheep at the follicular phase and luteal phase were named WF (n = 3) and WL (n = 3), respectively. All tissues were snap-frozen in liquid nitrogen and then stored at -80 °C for RNA extraction.

RNA Extraction, Library Construction, and RNA-Seq
Total RNA was extracted from the oviduct tissues of 12 ewes. TRIzol reagent (Invitrogen, Carlsbad, California, USA) was used for total RNA extraction, according to the manufacturer's instructions. To obtain high-quality RNA, 1% agarose electrophoresis and an Agilent 2100 Bioanalyzer (Agilent Technologies, Santa Clara, CA, USA) were used to examine the integrity and concentration of the extracted RNA. The purity of isolated RNA was also ensured using an Agilent RNA 6000 Nano Kit (Agilent Technologies).
The circRNA library was constructed with 3 µg of total RNA using the NEBNext ® Ultra™ RNA Library Prep Kit for Illumina ® (#E7530L, NEB, USA) according to the manufacturer's recommendations. Ribosomal RNAs (rRNAs) were removed from the total RNA using a Ribo-Zero™ Gold Kit (Epicentre, Madison, WI, USA). CircRNAs were randomly fragmented and reverse transcribed into cDNA with random primers. Then, a fragmentation buffer was added to break the RNA into fragments, which were used as templates to synthesize the rst strand of complementary DNA. The obtained double-stranded cDNA was processed with end-repair, the addition of base A and sequencing adaptors, and Uracil-N-Glycosylase (UNG) enzyme digestion. Subsequently, the polymerase chain reaction was performed to construct a circRNA library.
The fragments with lengths of 18-30 nt, which were obtained from total RNA through the gel separation technique, were used as templates to synthesize the rst strand of complementary DNA (cDNA). The second strand of cDNA was also synthesized in the presence of deoxynucleoside triphosphates (dNTPs), ribonuclease H, and DNA polymerase I. Then the obtained double-stranded cDNA was processed with end-repair, the addition of base A and sequencing adaptors, and uracil-Nglycosylase (UNG) enzyme digestion. Finally, polymerase chain reaction was conducted to build the miRNA library. Raw data of the RNA-Seq have been deposited in the SRA database (Accession number: PRJNA658731).
In addition, all the sequencing works were conducted in Novogene Bioinformatics Technology Co. Ltd. (Beijing, China).

Analysis of circRNA sequencing data
Clean reads of circRNAs were mapped to the reference genome (Oar_v3.1) by the HiSAT2 v. 2.0.4 (https://ccb.jhu.edu/software/hisat2/index.shtml) alignment method, both the sheep reference genome and genome annotation le were downloaded from ENSEMBL (http://www.ensembl.org/index.html). CIRI is an e cient and fast tool to identify circRNAs [82]. The BWA-MEM algorithm was used to conduct a sequence splitting comparison to ensure the reliability of other circRNAs, and then the SAM le was scanned to nd PCC (paired chiastic clipping) and PEM (paired-end mapping) sites, and GT-AG splicing signals [83]. Moreover, we used dynamic programming algorithm to re-align the sequence with the junction site. CircRNAs were annotated against the circBase. Statistical analysis was performed on the identi ed circRNA types, chromosome distribution, and length distribution. Relative expression of circRNAs were analyzed by TPM (transcripts per million reads) [84]. The DEseq2 package which was based on negative binomial distribution was used to identify DE circRNAs across groups [85]. For the purpose of screening key circRNAs, the thresholds of fold change > 2 and q-value < 0.05 (corrected by FDR) were set to identify DE circRNAs. In addition, miRanda v3.3a was used to predict the miRNA binding site of circRNAs (Supplemental Table S7) [86].
Analysis of miRNA sequencing data Several criteria were conducted to obtain clean miRNA reads, including removing reads without a 3′ adapter, reads without insert fragment, reads with lengths beyond the normal range, raw reads containing too much A/T, and some low-quality reads using in-house scripts. Then, the clean data of miRNA were matched against the sheep reference genome (Oar v3.1) by Bowtie v1.1.2 [87]. Relative expression of miRNAs were analyzed by TPM. The DEseq2 package which was based on negative binomial distribution was used to identify DE miRNAs across groups. Also, the thresholds of fold change > 2 and q-value < 0.05 were set to identify DE circRNAs. In addition, miRanda v3.3a, PITA and RNAhybrid v2.1.2 was used to predict the target genes of miRNAs (Supplemental Table S7).

Bioinformatics Analysis
Gene ontology analysis was performed on host genes of circRNA and predicted target genes of miRNA by using GOseq R package [88]. Kyoto Encyclopedia of Genes and Genomes annotations on host genes of circRNA and predicted target genes of miRNA were also conducted based on KEGG database (http://www.genome.jp). The hypergeometric test method was applied to assess signi cantly enriched GO terms and KEGG pathways, and those with p < 0.05, were thought to be signi cantly enriched.

Validation of the Expression of circRNAs and miRNAs
To further con rm the circRNA and miRNA sequencing data, ve DE circRNAs and ve DE miRNAs were selected randomly. We designed forward and reverse primers encompassing circRNA-speci c, back-splice junctions for each candidate circRNA (Supplemental Table S8). The reverse transcription of circRNA was performed using PrimeScript™ RT reagent kit (TaKaRa, Dalian, China). Then, RT-qPCR 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.
And, the reverse primers of miRNAs were designed using tailing reaction, which increase the accuracy and speci city of detection (Supplemental Table S8). The forward primer is included in miRcute Plus miRNA qPCR Kit (SYBR Green). Afterwards, the reverse transcription of miRNA was performed using miRcute Plus miRNA First-Strand cDNA Kit (TIANGEN, Beijing, China), followed by the use of miRcute Plus miRNA qPCR Kit (SYBR Green) (TIANGEN, Beijing, China) to conduct RT-qPCR through the Roche Light Cycler®480II.
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. The β-Actin and U6 small nuclear RNA were used as internal control to normalize target gene expression, respectively (Supplemental Table S8). The results obtained from RT-PCR were calculated using the 2 −∆∆Ct method [89] and then processed by SPSS22.0. Finally, PCR products were gel extracted and subjected to Sanger sequencing.

Statistical Analyses
All data are presented as the mean ± SEM. Student's t tests were performed to compare means, and p < 0.05 was considered statistically signi cant.     Validation results of four selected circRNAs in MF and ML by RT-qPCR.