The AFWS used in this study were raised in the AFWS Stud Farm of Inner Mongolia Autonomous Region and fed according to the farm’s feeding plan. Twelve healthy AFWS ewes of similar age (3–5 years old), body weight (55–60 kg), and body size were selected. Estrus of the 12 ewes was synchronized, and artificial insemination was performed during September. The ewes and lambs were anesthetized with sodium pentobarbital at a dose of 25 mg/kg by intravenous injection. After samples collection, the ewes and born lambs were released, whereas the fetuses from E90d and E120d were placed, still under anesthesia, inside a closed chamber, which was filled with 20% carbon dioxide per minute. When gas concentration had reached 80%, the fetuses died. The anesthesia procedure was performed following published protocols [67, 68].
The 2-cm-diameter skin tissue samples (about 0.5–1.0 g per fetus/lamb) were collected from the shoulder area at the three developmental stages (E90d, E120d, and Birth), three individuals for each stage, nine in total. The collected samples were placed into clean RNAase-free Eppendorf tubes and stored under liquid nitrogen, pending total RNA extraction. Skin samples were also fixed in 4% formaldehyde, and paraffin sections were prepared and stained with H&E for histological observations.
RNA isolation and quality assessment
To extract total RNA from the nine samples, TRIzol reagent (Life Technologies, CA, USA) was used. RNase-free DNase (Tiangen, Beijing, China) was used to remove DNA contamination from the extracted RNA. RNA degradation and contamination were monitored by 1% agarose gel electrophoresis and RNA purity was measured at an OD260/280, using a NanoDrop ND-2000 instrument (Thermo Fisher Scientific, MA, USA). We also assessed RNA integrity by testing the RIN of the samples.
High-throughput whole transcriptome sequencing and subsequent bioinformatics analyses were performed by Annoroad Technologies (Beijing, China) as follows: A total of 3 μg RNA per sample were used for circRNA sample preparation. The Ribo-ZeroTM Gold Kit was used to remove rRNA from the samples, and different index tags were selected to build the library according to the specifications of NEB Next Ultra Directional RNA Library Prep Kit for Illumina (NEB, Ispawich, USA). The specific steps of library construction were as follows: Ribosomal RNA was removed using a kit, RNase R was added to remove linear RNA. Fragmentation Buffer was added to the reaction system to fragment the RNA, and then this fragmented RNA was used as a template for first strand cDNA synthesis, using random primers (Random Hexamers). Second strand cDNA was synthesized by adding buffer, dNTPs, RNase H, and DNA Polymerase I. After purification by QiaQuick PCR kit and elution with EB buffer, the following steps were performed: repair ending, adenine addition, sequencing linker addition, and target size fragments’ recovery (approximately 350 bp) by agarose gel electrophoresis. Uracil N-glycosylase (UNG) was then added to digest the DNA strand prior to PCR amplification. Finally, agarose gel electrophoresis was used to recover the DNA fragments of the target size. The constructed library was sequenced using Illumina X Ten and PE150 sequencing strategy.
Sequencing analysis of circRNA
Sheep genome oar_v4.0 was selected as reference genome for comparison with the RNA-Seq data. Reads were mapped to the reference genome using the BWA-MEM method, which is fast and efficient in aligning reads, and allows mapping fragment reads to genomes as well. The raw reads generated by Illumina sequencing were processed to create clean reads by several processes, including de-junction contamination and removal of rRNA. For mapping, the BWA-MEM algorithm was first used for sequence splitting and alignment. The resulting Sam files were scanned in search of PCC (paid Chinese clipping) and PEM (paid end mapping) sites, as well as GT-AG splicing signals. Finally, sequences with junction sites were re-aligned with dynamic programming algorithm to ensure the reliability of circRNA identification. CIRI , an efficient and rapid tool for circRNA recognition, was also used. All subsequent analyses were based on the clean reads. The process of analyzing the circRNAs sequencing information in this study was divided into seven parts: (1) sequencing data quality control, (2) data alignment analysis, (3) circRNAs identification and classification, (4) circRNAs characteristics analysis, (5) circRNAs differential analysis, (6) differentially-expressed circRNAs source genes functions, and (7) miRNA molecular sponge analysis.
Identification of differentially-expressed circRNAs
We used SRPBM as a normalization method to quantify the expression of circRNA. The DEseq2  software was used to analyze the differentially-expressed circRNAs. The three fetuses/lambs at each stage were used as biological replicates. Differentially-expressed circRNAs were detected by comparing one stage with another. CircRNAs with P < 0.05 and absolute fold-change values of > 1.5 in any of the pairwise comparisons were considered to be significantly differentially-expressed. Upregulated and downregulated circRNA numbers were eventually obtained. The calculation formula of SRPBM is: , where SR is the number of spliced reads, and N is the total number of mapped reads in the sample.
Gene Ontology and Kyoto Encyclopedia of Genes and Genomes pathway enrichment analyses
Gene Ontology and KEGG pathway analyses were used to annotate the source genes of differentially-expressed circRNAs. The Blast2GO method  was used for GO functional analysis, while KOBAS software was used to test the statistical enrichment of differential gene expression in the KEGG pathway analysis . Enrichment was considered significant in the GO term and KEGG pathway analyses when P < 0.05.
Prediction of miRNAs targeted by circRNA
To explore the functions of circRNAs, predict the targeting relationship, and thus predict which of the circRNAs function as miRNAs sponges, we used miRanda V.3.3a (http://www.microrna.org/microrna/home.do) . In view of known reports and extractability of the sequences, we selected only CLASSIC and ANTISENSE circRNA types for the prediction of the miRNA targeting relationship.
Experimental validation of circRNAs
Quantitative real-time PCR (qRT-PCR) was used to validate circRNAs expression. We randomly selected seven circRNAs for validation. The expression levels of the selected circRNAs were normalized against the expression of a housekeeping gene, GAPDH. Primers were designed and synthesized by Sangon Biotech Co., Ltd. (Shanghai, China). Total RNA was converted into cDNA using random hexamers with Transcriptor First Strand cDNA Synthesis Kit (Roche, Australia). The qRT-PCR analysis was carried out in triplicate with iTaqTMUniversal [email protected] Green Supermix (Bio-Rad, CA, USA) on a Bio-Rad CFX96 instrument (Bio-Rad, CA, USA). The total 20 μL reaction mixture contained 10 μL 2×iTaqTMUniversal SYBR@ Green Supermix, 1 μL cDNA, 8 μL ddH2O, and 0.5 µL each of forward and reverse primers. The following program was used: 95 °C for 10 min; 45 cycles of 95 °C for 10 s, 60 °C for 10 s, and 72 °C for 10 s; 72 °C for 6 min. The 2-ΔΔCt method was used to analyze the relative expression levels of the selected circRNAs.
To determine the resistance of the selected seven circRNAs to RNase R digestion, total RNA and RNase R (Geneseed Biotech, Guangzhou, China) were mixed together. The mix was incubated at 37°C for 15 minutes, cDNA was then synthesized, and expression level of circRNAs was finally detected by qRT-PCR.