The protocol for collecting ruminal biopsies was approved by the Institutional Animal Care and Use Committee of the College of Animal Science and Technology of Sichuan Agricultural University, Sichuan, China, under the permit of No. DKY-2019302081.

The animals used in this study were raised at the Jianyang Dageda husbandry Co., Ltd., Sichuan, China (34°4′ N; 104°5′ E). Jianzhou Da’er Goat is a meat breed with a large body size, fast growth, high resistance to rough feeding, and desired reproductive performance (China National Commission of Animal Genetic Resources, 2011). All the goats were kept indoors in a naturally illuminated and ventilated shed with playgrounds. Estrous cycles of all the does were synchronized, and after estrus detection, does were inseminated (Day 0). Pregnancy status and fetal number were determined on Day 30 with transrectal ultrasound. Twelve pregnant multiparous goat (2–3 years old), homogenous in body condition score (between 3.0 and 3.5 on a scale of 1–5), were included in this study. The female fetuses were identified by the PCR amplification using the SE47/SE48 sex identification primers (Ennis and Gallagher, 1994). The pregnant goats were fed with 65% forage and 35% concentrate twice daily at 08:00–9:00 h and 17:00–18:00 h. The kids were breast-fed by their mothers and received concentrate supplement ad libitum. The kids were weaned at 60 days after birth. The feed nutrition level was designed according to the local feeding standard of Jianzhou Da’er Goat (DB51/T 1750-2014). The weaned goat kids received 70% forage and 30% concentrate twice daily.

The ages for sampling were selected according to the level of rumen development. Considering that the rumen develops gradually during fetal life, and at the end of gestation, the rumen papillae are entirely developed, the fetuses on days 60 and 135 of gestation were selected. After birth, considering that the diet affects rumen development, female kids were selected before and after weaning (30 and 150-days-old) for further study. Rumen tissues during the fetal period were collected at 60 and 135 days of gestation (groups F60 and F135, respectively; n = 3 in each group), and after birth, during the prepubertal period, before and after weaning (at 30 and 150 days of age; groups BW30 and AW150, respectively; n = 3 in each group). The pregnant does from groups F60 and F135 were anesthetized (intravenous injection of Zoletil 50^®, 4 mg/kg, Carros, France) according to Zhong et al. (2017), followed by the cesarean section to separate fetus. The rumen of fetuses was sampled by a veterinary surgery immediately after a cesarean. The kids from groups BW30 and AW150 were euthanized with anesthesia followed by immediate bleeding. Thereafter, the rumen tissues were dissected, and each sample was divided into two similar pieces. One piece was immediately frozen in liquid nitrogen and then stored at −80^°C for the following RNA-Seq analysis. The other half of the rumen sample was rinsed in phosphate buffer saline (PBS) solution and fixed with 4% formaldehyde for hematoxylin and eosin staining.

The fixed rumen samples were dehydrated in increasingly concentrated ethanol (80°, 90°, and 95°), dehydrated in xylene, and embedded in paraffin wax. Sections (5-μm) were cut using a microtome (LeicaRM2016, Leica Microsystems, Wetzlar, Germany) perpendicular to the longitudinal axis of the rumen and placed on glass slides. The hematoxylin-eosin stained sections were captured using an Olympus BX51 microscope (Olympus, Tokyo, Japan) equipped with a Nikon DS-Fi1 camera (Nikon, Tokyo, Japan).

Total RNA was extracted from rumen tissues using the TRIzol reagent (Invitrogen, CA, United States). The quality, quantity, and integrity of RNA were detected by the NanoDrop spectrophotometer (Thermo Fisher Scientific, DE, United States) and Agilent 2100 Bioanalyzer (California, United States). Approximately 5 μg of total RNA was used to remove ribosomal RNAs by the Epicenter Ribo-zero™ rRNA Removal Kit (Epicenter, Madison, WI, United States) and then randomly interrupt rRNA-depleted RNA by adding the Fragmentation Buffer. Sequencing libraries were generated with the NEB Next^® Ultra™ Directional RNA Library Prep Kit (New England Biolabs, Beverly, MA, United States). Lastly, the libraries were qualified with the Agilent 2100 Bioanalyzer. The clustering of the index-coded samples was performed on a cBot Cluster Generation System using the HiSeq PE Cluster Kit v4 cBot (Illumina, San Diego, CA, United States). After cluster generation, the library preparations were sequenced on the Illumina HiSeq 4000 platform, and 150 bp paired-end reads were generated.

Raw reads were filtered to acquire clean reads by removing adapter-containing reads, poly-N-containing reads (over 5%), and low-quality reads using the FASTQC (Cock et al., 2010). The clean reads were aligned to the goat reference genome (ARS1, GCF_001704415.1) using the TopHat v2.0.9. To predict circRNAs, the sequence alignment map format file was generated by aligning to the ARS1 using the BWA-MEM algorithm (Houtgast et al., 2018), which was performed in the CIRI software (Gao et al., 2015). Prediction detected reverse concatenated reads (at least two unique back-spliced reads whereas obey the GU AG rule) using default parameters.

The relative expression of each circRNA was estimated using the mapped back-splicing junction (BSJ) reads and normalized by transcripts per million reads (TPM) (Arcinas et al., 2019). The DESeq2 package was applied to identify differentially expressed circRNAs (DECs) in each comparison between two different developmental stages. The p-value was adjusted by the Benjamini-Hochberg approach for controlling the false discovery rate (FDR < 0.05). The circRNAs with adjusted FDR < 0.05 and | log₂(Fold Change)| > 1 were considered as DECs.

The GO and KEGG analyses were performed using the DAVID tool (Huang et al., 2009), and KOBAS software (Xie et al., 2011) was used to explore the functions of the parental genes of the DECs identified in the rumen tissues. The corrected p-value < 0.05 was used as the threshold to define the significant enriched GO terms or KEGG pathways.

CircRNAs are proposed to regulate gene expression through competing endogenous RNAs acting as miRNA sponges (Hansen et al., 2013). The target miRNAs, which could be bound to the circRNAs, were predicted using the miRnada and regRNA 2.0 database (Chang et al., 2013) with a set of default parameters. Then, the goat rumen miRNAs from our previous studies (Zhong et al., 2017), which were identical to the above predicted miRNAs, were used to construct the circRNA-miRNA network by the Cytoscape software (Shannon et al., 2003).

The circRNAs sequencing data was validated by randomly selecting six DECs to assess the reliability of the obtained results. Firstly, divergent primers (Supplementary Table 1) were designed and used to detect the expression by the RT-qPCR, and the BSJ sites of each circRNA were confirmed by Sanger sequencing. The expression of the selected circRNA was determined by the 2^–ΔΔCt method with the internal GAPDH gene. Amplified fragments of PCR were visualized using 2% agarose gel electrophoresis and then performed cloning sequencing on an ABI 3730XL DNA analyzer (Applied Biosystems, CA, United States). Secondly, the resistance of selected circRNAs to exonuclease digestion was assessed by RNase R assay. Total RNA was treated with the RNase R reagent according to the manufacturer’s instructions (Geneseed, Guangzhou, China). cDNA was synthesized using the random primer and OligdT primer, respectively. Then, PCR was performed to compare the exonuclease resistance of circular and linear transcripts. Data were displayed as mean ± SEM and were considered statistically significant when the P-value was less than 0.05.

Evident stratification of the rumen wall was observed in the rumen tissues from the F60 group, and the rumen papilla was present (Figure 1A). The papilla of the rumen became longer, and the layers of the rumen wall were thicker in the F135 group than those in theF60 group (Figure 1B). The length and width of the rumen papilla were longer and wider, and the thickness of the muscular layer was wider in the BW60 than in the AW150 group (Figures 1C,D).

A total of 1,492,819,832 clean reads were obtained from the rumen tissues collected from the groups of F60, F135, BW30, and AW150, to capture the expression landscape of circRNAs during rumen development. Each library generated more than 11 GB clean bases after removing the adapter sequences and low-quality reads (Table 1). Over 91% of the valid reads were successfully mapped to the goat reference genome ARS1. A total of 11,149 circRNAs were identified, some of them are listed (Supplementary Table 2), including 8,388 circRNAs expressed in the fetal rumens (F60 and F135), 7,157 circRNAs in the prepubertal rumens (BW30 and AW150), and 1,949 common circRNAs across the four developmental stages (Figure 2A). These circRNAs identified in the rumen tissues were mainly derived from exons (up to 84.6%, Figure 2B) and displayed a broad distribution on each chromosome with the highest number on chromosome 10 (Chr10, 663) but the lowest number on Chr27 (145, Figure 2C). We further analyzed the length distribution of the circRNAs, ranging from 400 to 1,400 nt among 8,966 circRNAs (80.42%, Figure 2D).

A total of 1,518 circRNAs were found to be differentially expressed in the rumen tissues between the four development stages. 76, 798, and 125 DECs were identified between the comparisons of F60 and F135, F135 and BW30, BW30 and AW150 groups, respectively (Supplementary Table 3), and their expression profiles are shown using the hierarchical cluster (Figure 3A). The result showed that the expression of DECs was clearly distinguished and clustered into the fetal (F60 and F135) and prepubertal periods (BW30 and AW150), represented by a well-consistent pattern among the three samples at each developmental stage. A Venn diagram was constructed to explore the interactions between the DECs in each comparison, represented by 1,518 rumen DECs (Figure 3B). It is noticed that 12, 58, 2, and 50 DECs were the stage-specific DECs in the relative period, and 517 DECs were jointly identified in the four developmental stages. Totally, 125 circRNAs were significantly altered in the AW150 rumens compared with the BW30 rumens. Of which, 93 circRNAs were significantly up-regulated, and 32 circRNAs were remarkably down-regulated in AW150 rumens (Figure 3C). Notably, we observed more numbers of up-regulated DECs in the comparisons between F60 vs. F135, BW30 vs. AW150 groups, whereas a large number of down-regulated DECs were identified in the comparison between F135 and BW30 groups (Figure 3C). In addition, we selected six circRNAs and compared them with human circRNAs according to their parental genes and sequences, and found that these six circRNAs had high sequence similarity with human circRNAs (Supplementary Figure 1).

GO term and KEGG pathway analyses were performed to explore the function of the parental genes of the DECs identified in the three comparisons between the following groups: F60 vs. F135, F135 vs. BW30, and BW30 vs. AW150. 80, 158, 144 GO terms were identified in these three comparisons (Supplementary Table 4). The biological progress was the largest component of the GO terms (Figure 4A). When comparing F60 with F135 groups, and BW30 with AW150 groups, the parental genes of DECs were involved in cellular protein metabolic process (GO:0044267, p < 0.01), regulation of cell division (GO:0051246, p < 0.01), and regulation of cell division (GO:0051302, p < 0.01). When comparing F135 with BW30 groups, the most enriched pathway was the cytoskeletal part (GO:0044430, p < 0.01).

KEGG pathway analysis can reveal molecular interactions and the pathways associated with genes. The MAPK signaling pathway related to cell proliferation and apoptosis was identified across the three comparisons (Figure 4B). In contrast, the Ras and Rap1 signaling pathways were enriched in F135 compared with BW30, and in BW30 relative to AW150 (Figure 4B). In addition, the Adherens junction was enriched in F60 compared with F135, and in F135 vs. BW30, with the highest enrich factor in this latter comparison. Enrichment information of all the KEGG pathways in three comparisons is shown in the Supplementary Table 5. Subsequently, the pathways of the stage-specific DECs in either F135 or AW150 were represented in the Supplementary Figure 1. Their parental genes were also significantly enriched in the MAPK, Rap1, and Ras signaling pathways in F135 and AW150 groups. Furthermore, we selected several circRNAs with high junction ratio and performed KEGG analysis on their parental genes. Interestingly, we found that some of them were involved in RNA transport and basic transcription factors (Supplementary Table 6).

In order to explore the role of the parental genes of circRNAs involved in the MAPK signaling pathway, we used RegRNA2.0 to predict binding miRNAs with a score of 155 and minimum free energy of 15 kcal/mol as the criteria. Twenty-five circRNAs were found to be combined with 71 binding miRNAs (Figure 5 and Supplementary Table 7). Among them, 19 circRNAs had multiple miRNA targets, and the circRNAs derived from MAP3K5 (circ9917-5) or DUSP16 (circ5813-1, circ5813-2, and circ5813-3) had the most number of target miRNAs.

Six DECs (circRNA25300-5, circRNA25471-2, circRNA16964-3, circRNA2068-2, circRNA1831-1, and circRNA2158-2) were randomly selected to validate the reliability of RNA-Seq data. the qPCR results of all the six circRNAs showed transcription tendencies highly similar to those in the RNA-Seq data (Figure 6A). The BSJ sequences of the six circRNAs were also verified by Sanger sequencing using divergent primers, represented by circRNA25300-5 (Figure 6B). Subsequently, we performed RT-PCR to compare the exonuclease resistance of circular and linear transcripts. All the six circRNAs showed a resistance to the RNase R digestion (Figure 6C). In addition, the PCR reaction, using the cDNA synthesized by OligdT primer, did not amplify any band because the circRNAs lacked 3′ poly-A tail. These verifications indicated that the circRNAs sequencing data were reliable.