Transcriptomic Profiling of Circular RNAs in the Goat Rumen During Fetal and Prepubertal Period

Circular RNAs (circRNAs) are key regulatory factors with vital functions in various biological activities. However, little has been reported concerning the genetic regulation of circRNAs during rumen development in goats. The aim of this study was to identify the genome-wide expression profiles of circRNAs in the rumen of goats during fetal development and before and after weaning. Histological morphology showed that from the fetal period (days 60 and 135 of gestation) to the prepuberal period (days 60 and 150 of age) the rumen papilla developed gradually, and the thickness of the rumen muscular layer increased. A total of 11,149 circRNAs were identified in the four development stages by RNA-sequencing. From this, 1,518 were differentially expressed circRNAs (DECs). Fifty-eight DECs were up-regulated from 60 to 135 days of gestation, and 93 from day 135 of pregnancy to 30 days after birth. A large proportion (598) of DECs were down-regulated from day 135 of gestation to 30 days after birth. The expression levels of six randomly selected circRNAs were validated by qPCR, and their back-splicing junction (BSJ) sites were also confirmed. Ontology and pathway analyses revealed that the parental genes of DECs were mainly involved in the signaling pathways related to cell proliferation and apoptosis. The interaction network of circRNAs with their target miRNAs showed its involvement in cell proliferation and apoptosis signaling pathways. In conclusion, we identified the genome-wide expression profiles of circRNAs in the rumen of goats during fetal development and before and after weaning. These results provide a basis for further study on the regulatory effect of circRNAs on the development of rumen tissues.

Introduction

The rumen is the main organ for digestion, absorption, and metabolism in ruminants. The development of the rumen includes anatomic changes (increase in rumen absolute mass and relative size in relation to the other stomach cavities, and growth of the rumen papillae), functional achievement (acquisition of the fermentation capacity and enzyme activity), and microbial colonization (development of the bacteria, fungi, protozoa, and archaea populations) (Wang et al., 2016). From fetal period to adulthood, the morphology and functionality of the rumen epithelium are in continuous development. During the fetal life, rumen function is not required as there is no income from the digestive system. During this period, maternal glucose is the primary energy source for the fetus, being even more critical during the last third of gestation when its greatest growth occurs (Trotta and Swanson, 2021). At the end of gestation, the rumen papillae are entirely developed, and its’ surface is keratinized (Garcia et al., 2012). During the first hours after birth, the ingestion of colostrum is the main energy source for the neonate (Nowak and Poindron, 2006), beginning the income of nutrients through the digestive system. The rumen of the neonate is covered by a smooth, stratified, and squamous epithelium without any prominent papillae, and its physiological functionality is incomplete as the microbiome is not established (Wang et al., 2016). During the first 4 weeks of age, milk is the main energy source for the offspring (Degen and Benjamin, 2003, 2005). Solid food consumption increases gradually with age, being the unique food source for the kids after weaning (Weary et al., 2008). Overall, during ruminant development, the changes in age, diet, and environment are accompanied by rumen growth, and metabolic and physiological adaptations, which are regulated by specific genes.

The studies on rumen development have been mainly focused on changes in its morphology, internal environment, establishment and maintenance of microflora, absorption and transport of nutrients (Lesmeister and Heinrichs, 2004; Jiao et al., 2015; Ragionieri et al., 2016; Wang et al., 2016; Do et al., 2019). However, there is limited information on the candidate genes/ncRNAs related to rumen development. Circular RNAs (circRNAs), a large category of recently known endogenous non-coding RNAs, are characterized as a closed-loop structure covalently linked by the 5′- and 3′- termini (Qu et al., 2015). CircRNAs are generated by the back-splicing of precursor messenger RNAs during the post-transcriptional processes, being more stable than linear RNAs (Jeck et al., 2013). Emerging evidences indicate that circRNAs modulate miRNA activity by competing with endogenous RNA or acting as microRNA sponges, regulating a variety of biological and disease processes (Hansen et al., 2013; Chen et al., 2017; Kristensen et al., 2018; Zhong et al., 2018). Numerous circRNAs have been identified in various tissues of farm animals using the high-throughput sequencing, including yak adipose tissue, goat endometrium, and pig longissimus dorsi muscle (Le Bras, 2019; Song et al., 2019; Wang et al., 2019; Zhang Y. et al., 2020). CircRNAs are ubiquitously expressed and abundant in animals, and they are differentially expressed in diverse types of cells, specific tissues, and developmental stages (Memczak et al., 2013; Ji et al., 2019; Li et al., 2021). Recent studies have identified candidates of miRNAs differently expressed, which are mainly involved in cell proliferation, growth, and apoptosis during the development of goat rumen (Zhong et al., 2017, 2020). Therefore, CircRNAs could participate in those processes, modulating miRNA activity.

Considering all this information, we hypothesized that the universally expressed circRNAs act as potential regulators of the rumen development in goats. The aim of this study was to identify the genome-wide expression profiles of circRNAs in the rumen of goats during fetal development and before and after weaning.

Materials and Methods

Ethics Approval and Consent to Participate

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.

Location, Animals, and General Management

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.

Experimental Groups and Sampling Procedures

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.

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

Rumen RNA Isolation, Library Preparation, and Sequencing

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.

Sequencing Data Analysis, circRNAs’ Identification, and Expression

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.

GO Enrichment, KEGG Pathway, and CircRNA-miRNA Network

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

Validation of CircRNAs by qRT-PCR, Sequencing, and RNase R Resistance

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.

Results

Histological Observation of the Goat Rumens

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

FIGURE 1

Figure 1. General view of goat rumen tissue during fetal and prepubertal periods. Photomicrograph of rumen wall of the fetal goat at F60 (A) and at F135 (B) days of gestation, and at 30 (BW30) and 150 days (AW150) of age (C,D), respectively. Bar, 200 μm.

Expression Profiles of CircRNAs in Goat Rumens

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

TABLE 1

Table 1. Summary of circRNAs sequencing results in goat rumen tissues at four developmental stages.

FIGURE 2

Figure 2. Properties of circRNAs in goat rumen tissue samples. (A) Venn diagram showing the number of circRNAs either shared between or uniquely expressed in goat rumen tissues during fetal and prepubertal periods. (B) Pie chart showing the prevalence of different circRNA types. circRNAs from introns are represented by the orange section, while the blue region indicates circRNAs from exons. (C) Distribution of circRNAs on the chromosomes of the goat reference genome. (D) Lengths of the circRNAs.

Differentially Expressed CircRNAs in Goat Rumens

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

FIGURE 3

Figure 3. (A) Hierarchical cluster analysis (heat map) represents all differentially expressed circRNAs (DECs) between F60 vs. F135, BW30 vs. AW150 groups. Red and green denoted high and low expressions, respectively. (B) Venn diagram showing the DECs at the four developmental stages. (C) Numbers of up-regulated or down-regulated circRNAs in goat rumens during fetal and prepubertal periods.

Biological Functions of the Parental Genes of Differentially Expressed CircRNAs

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

FIGURE 4

Figure 4. (A) The significantly enriched GO terms (p < 0.01) in the biological process, molecular function, and cellular component categories. (B) Top 20 KEGG enrichment pathways.

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

Construction of the CircRNA-miRNA Network

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.

FIGURE 5

Figure 5. The circRNA-miRNA interaction network in the goat rumen during fetal and prepubertal periods, Including 27 circRNAs and 71 miRNAs.

Verification of the Obtained CircRNAs

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.

FIGURE 6

Figure 6. Validation of DECs in goat rumens in the fetal and prepubertal periods. (A) Expression profiles of the six randomly selected circRNAs based on RNA-Seq and qPCR. (B) Schematic for back-splicing junction (BSJ) validation of circRNAs, represented by circ25300-5. The (BSJ) sequence was confirmed by Sanger sequencing. (C) RNase R assay for the circularity of the selected circRNAs by resistance to exonuclease digestion. OligdT shows the template cDNA for RT-PCR synthesized with the OligdT primer, while others were generated using random primers. +, indicates the samples were treated with exonuclease.

Discussion

To the best of our knowledge, this is the first study that identified the genome-wide expression profiles of rumen circRNAs during fetal development and before and after weaning in ruminants. The development of the rumen changes with age and is determined by genetic regulation (Zhong et al., 2020), diet composition (Xue et al., 2019), and microbial colonization (Lin et al., 2019). In the early stage of rumen development at the second week of newborn kids, the volatile fatty acids generated by the microbiota can interact with the host transcriptome and mRNA to influence the metabolism of rumen tissue and epithelial development (Lin et al., 2019). In terms of transcriptome analysis, the involvement of miRNA and mRNA in the regulation of rumen development has been reported (Zhong et al., 2017, 2020; Do et al., 2019), but the regulatory mechanism of circRNAs in rumen development has rarely been reported (Do et al., 2019). In the present study, we identified a total of 11,149 circRNAs during the fetal (at 60 and 135 days of gestation) and prepubertal (at 30 and 150 days of age) periods in goats, among which 1,518 circRNAs were differentially expressed. These results suggest that these circRNAs may play an active role in rumen development. We found that most of the DECs were up-regulated during the development of the rumen in the fetal period (F60 and F135), and in the pre- and post-weaning period (BW30 and AW150). In fetal and prepubertal periods, the number of down-regulated circRNAs in the rumen tissues was more than that of up-regulated circRNAs. We proposed that those circRNAs may be involved in the adaptation of rumen function, together with the changes in age and diet.

Previous studies have shown that rumen development is related to specific biological pathways, including pteridine-containing compounds, calcium signaling, and transcription factor PPAR-α. In this study, we found that the parental genes of DECs were significantly enriched in the signal pathways related to cell proliferation and apoptosis, such as MAPK (Han et al., 2020), Ras (Takino et al., 2019), FoxO (Xie et al., 2019), and Rap1 signaling pathways (Zhou et al., 2020). In addition, 83 circRNAs were only expressed in AW150 when we compared BW30 vs. AW150 groups (Supplementary Table 3), suggesting that the functions of circRNAs change with weaning and/or with age. The parental genes in the rumen tissue from the AW150 group were also significantly enriched in MAPK, RAS, and RAP1 signaling pathways (Supplementary Figure 2), indicating that such changes were related to rumen development. In our previous study, the target genes of differently expressed miRNAs identified in goat rumens were mainly enriched with similar signal pathways, such as MAPK, Ras, and Jak-STAT signaling pathways (Zhong et al., 2020). In view of the KEGG analysis, we hypothesized that these related pathways might play an essential role in the development of rumen epithelial cells.

CircRNAs could act as the sponges of miRNAs by competitively binding to miRNAs and modulating the activities of their target genes (Li et al., 2018). In fact, one circ9917-5, derived from the MAP3K5 gene, had 13 miRNA targets (Figure 5). We proposed that this circRNA could play its function as a miRNA sponge. In gastric cancer cells, circDUSP16 acted as a sponge of miR-145-5p and affected cell proliferation (Zhang Z. et al., 2020). Interestingly, three circRNAs (circ5813-1, circ5813-2, and circ5813-3) were derived by cyclization of the DUSP16 transcripts, which could bind with chi-miR-145-5p in this study. Thus, we suggest that these three circRNAs may play biological roles in the proliferation of rumen epithelial cells. In addition, circ5813-1, circ5813-2, and circ5813-3 were found to have a binding site of miR-148a-3p (Supplementary Table 6). In the previous study, we found that chi-miR-148a-3p could regulate the proliferation of gastric mucosal epithelial cells (Zhong et al., 2020). Therefore, it was speculated that these circRNAs could act as molecular sponges of miR-148a-3p in the rumen tissues. Several miRNAs have been identified as the key regulators of early rumen development in cattle, including miR-29b, miR-145, and miR-493 (Do et al., 2019). In the present study, we found that chi-miR-29a-5p, chi-miR-145-5p, chi-miR-145-3p, and chi-miR-493-3p could bind to more than one circRNAs (Supplementary Table 6). We speculated that those circRNAs could bind with miRNAs during rumen development and play potential regulatory roles during rumen development in goats.

The molecular changes observed in the rumen tissue from goats during fetal and prepubertal periods were accompanied by morphological modifications. In this study, it was observed that the rumen papilla was already present in the fetal period. During the prepubertal period, the rumen papilla and the thickness of the muscle layer increased rapidly, probably due to the stimulation of the diet and the increase in age. In this sense, the rumen papilla can increase the surface area of the rumen epithelium, thereby prompting the absorption of short-chain fatty acids and minerals in the rumen (Lavker and Matoltsy, 1970). Previous studies have found that volatile fatty acids in the diet can enhance rumen development (Gorka et al., 2011). These results, as a whole, indicate that circRNAs could regulate the rumen development in goats.

In summary, a total of 11,149 circRNAs and 1,518 DECs were found in four developmental stages of goat rumens. Among them, the number of DECs between F135 and BW30 was the largest, suggesting that rumen development before and after birth was regulated by circRNA. GO and KEGG analyses showed that the parental genes of circRNAs were involved in cell division, proliferation, and apoptosis. Meanwhile, we predicted some of the circRNA target miRNAs and constructed a circRNA-miRNA network to expand our understanding. Circ5813-1 was predicted to combine with 12 binding miRNAs, so it could be a key circRNA involved in the process of rumen development in goats. Overall, our results can support further studies on the involvement of circRNA in the regulatory mechanisms of rumen development of goats.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi.nlm. nih.gov/bioproject/, PRJNA720177.

Ethics Statement

The animal study was reviewed and approved by the Institutional Animal Care and Use Committee of the College of Animal Science and Technology of Sichuan Agricultural University.

Author Contributions

TZ and LN conceived, designed the experiments, and contributed to reagents and materials. CW and MW performed the experiments. CW, SZ, and JC analyzed the data. LW, DD, JG, LL, and HZ participated in sample collection. CW and TZ wrote the manuscript. AF-M contributed to the interpretation and made substantial contributions to the manuscript. All authors have read and approved the final manuscript.

Funding

This work was supported by the Dual project of discipline construction of Sichuan Agricultural University.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Acknowledgments

We are grateful to Jianlin Han from International Livestock Research Institute (ILRI) for processing data and proofreading the article.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fphys.2022.858991/full#supplementary-material

Supplementary Figure 1 | Sequence alignment of goat circRNA and human circRNA.

Supplementary Figure 2 | Pathway analysis of the parental genes of differentially expressed circRNAs (DECs). (A) The top 20 pathways of the parental genes of the stage-specific DECs in F135. (B) The top 20 pathways of the parental genes of the stage-specific DECs in AW150.

References

Arcinas, C., Tan, W., Fang, W., Desai, T., Teh, D., Degirmenci, U., et al. (2019). Adipose circular RNAs exhibit dynamic regulation in obesity and functional role in adipogenesis. Nat. Metab. 1, 688–703. doi: 10.1038/s42255-019-0078-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Chang, T. H., Huang, H. Y., Hsu, J. B. K., Weng, S. L., Horng, J. T., and Huang, H. D. (2013). An enhanced computational platform for investigating the roles of regulatory RNA and for identifying functional RNA motifs. BMC Bioinformatics 14:S4. doi: 10.1186/1471-2105-14-S2-S4

PubMed Abstract | CrossRef Full Text | Google Scholar

Chen, Y. G., Satpathy, A. T., and Chang, H. Y. (2017). Gene regulation in the immune system by long noncoding RNAs. Nat. Immunol. 18, 962–972.

Google Scholar

China National Commission of Animal Genetic Resources (2011). Animal Genetic Resources in China: Sheep and Goats. Beijing: China Agricultural Press.

Google Scholar

Cock, P. J. A., Fields, C. J., Goto, N., Heuer, M. L., and Rice, P. M. (2010). The Sanger FASTQ file format for sequences with quality scores, and the Solexa/Illumina FASTQ variants. Nucleic Acids Res. 38, 1767–1771. doi: 10.1093/nar/gkp1137

PubMed Abstract | CrossRef Full Text | Google Scholar

Degen, A. A., and Benjamin, R. W. (2003). Milk intake and growth rate of Awassi lambs sucking ewes grazing on natural pasture in the semi-arid Negev. Anim. Sci. 76, 455–460. doi: 10.1017/S1357729800058677

CrossRef Full Text | Google Scholar

Degen, A. A., and Benjamin, R. W. (2005). Milk and herbage intakes and growth rate of lambs from 32 to 130 days of age raised on natural pasture in the semi-arid Negev. Small Ruminant Res. 58, 39–45.

Google Scholar

Do, D. N., Dudemaine, P. L., Fomenky, B. E., and Ibeagha-Awemu, E. M. (2019). Integration of miRNA weighted gene co-expression network and miRNA-mRNA co-expression analyses reveals potential regulatory functions of miRNAs in calf rumen development. Genomics 111, 849–859. doi: 10.1016/j.ygeno.2018.05.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Ennis, S., and Gallagher, T. F. (1994). A PCR-based sex-determination assay in cattle based on the bovine amelogenin locus. Anim. Genet. 25, 425–427. doi: 10.1111/j.1365-2052.1994.tb00533.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Gao, Y., Wang, J., and Zhao, F. (2015). CIRI: an efficient and unbiased algorithm for de novo circular RNA identification. Genome Biol. 16:4. doi: 10.1186/s13059-014-0571-3

PubMed Abstract | CrossRef Full Text | Google Scholar

Garcia, A., Masot, J., Franco, A., Gazquez, A., and Redondo, E. (2012). Histomorphometric and immunohistochemical study of the goat rumen during prenatal development. Anat. Rec. (Hoboken) 295, 776–785. doi: 10.1002/ar.22431

PubMed Abstract | CrossRef Full Text | Google Scholar

Gorka, P., Kowalski, Z. M., Pietrzak, P., Kotunia, A., Jagusiak, W., Holst, J. J., et al. (2011). Effect of method of delivery of sodium butyrate on rumen development in newborn calves. J. Dairy Sci. 94, 5578–5588. doi: 10.3168/jds.2011-4166

PubMed Abstract | CrossRef Full Text | Google Scholar

Han, F., Jiang, H., Qu, W., and Rui, Y. J. (2020). KLF11 protects chondrocytes via inhibiting p38 MAPK signaling pathway. Eur. Rev. Med. Pharmacol. Sci. 24, 6505–6516. doi: 10.26355/eurrev_202006_21634

CrossRef Full Text | Google Scholar

Hansen, T. B., Jensen, T. I., Clausen, B. H., Bramsen, J. B., Finsen, B., Damgaard, C. K., et al. (2013). Natural RNA circles function as efficient microRNA sponges. Nature 495, 384–388. doi: 10.1038/nature11993

PubMed Abstract | CrossRef Full Text | Google Scholar

Houtgast, E. J., Sima, V. M., Bertels, K., and Al-Ars, Z. (2018). Hardware acceleration of BWA-MEM genomic short read mapping for longer read lengths. Comput. Biol. Chem. 75, 54–64. doi: 10.1016/j.compbiolchem.2018.03.024

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, D. W., Sherman, B. T., and Lempicki, R. A. (2009). Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat. Protoc. 4, 44–57. doi: 10.1038/nprot.2008.211

PubMed Abstract | CrossRef Full Text | Google Scholar

Jeck, W. R., Sorrentino, J. A., Wang, K., Slevin, M. K., Burd, C. E., Liu, J., et al. (2013). Circular RNAs are abundant, conserved, and associated with ALU repeats. RNA 19, 141–157. doi: 10.1261/rna.035667.112

PubMed Abstract | CrossRef Full Text | Google Scholar

Ji, P., Wu, W., Chen, S., Zheng, Y., Zhou, L., Zhang, J., et al. (2019). Expanded expression landscape and prioritization of circular RNAs in mammals. Cell Rep. 26, 3444–3460 e5. doi: 10.1016/j.celrep.2019.02.078

PubMed Abstract | CrossRef Full Text | Google Scholar

Jiao, J., Huang, J., Zhou, C., and Tan, Z. (2015). Taxonomic identification of ruminal epithelial bacterial diversity during rumen development in goats. Appl. Environ. Microbiol. 81, 3502–3509. doi: 10.1128/AEM.00203-15

PubMed Abstract | CrossRef Full Text | Google Scholar

Kristensen, L. S., Hansen, T. B., Veno, M. T., and Kjems, J. (2018). Circular RNAs in cancer: opportunities and challenges in the field. Oncogene 37, 555–565. doi: 10.1038/onc.2017.361

PubMed Abstract | CrossRef Full Text | Google Scholar

Lavker, R. M., and Matoltsy, A. G. (1970). Formation of horny cells: the fate of cell organelles and differentiation products in ruminal epithelium. J. Cell Biol. 44, 501–512. doi: 10.1083/jcb.44.3.501

PubMed Abstract | CrossRef Full Text | Google Scholar

Le Bras, A. (2019). Repertoires of circRNAs in mammals. Lab Anim. (NY) 48:138. doi: 10.1038/s41684-019-0300-3

CrossRef Full Text | Google Scholar

Lesmeister, K. E., and Heinrichs, A. J. (2004). Effects of corn processing on growth characteristics, rumen development, and rumen parameters in neonatal dairy calves. J. Dairy Sci. 87, 3439–3450. doi: 10.3168/jds.S0022-0302(04)73479-7

CrossRef Full Text | Google Scholar

Li, T. T., Luo, R. R., Wang, X., Wang, H. H., Zhao, X. X., Guo, Y. X., et al. (2021). Unraveling stage-dependent expression patterns of circular RNAs and their related ceRNA modulation in ovine postnatal testis development. Front. Cell. Dev. Biol. 9:627439. doi: 10.3389/fcell.2021.627439

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, X., Yang, L., and Chen, L. L. (2018). The biogenesis, functions, and challenges of circular RNAs. Mol. Cell 71, 428–442. doi: 10.1016/j.molcel.2018.06.034

PubMed Abstract | CrossRef Full Text | Google Scholar

Lin, L., Xie, F., Sun, D., Liu, J., Zhu, W., and Mao, S. (2019). Ruminal microbiome-host crosstalk stimulates the development of the ruminal epithelium in a lamb model. Microbiome 7:83. doi: 10.1186/s40168-019-0701-y

PubMed Abstract | CrossRef Full Text | Google Scholar

Memczak, S., Jens, M., Elefsinioti, A., Torti, F., Krueger, J., Rybak, A., et al. (2013). Circular RNAs are a large class of animal RNAs with regulatory potency. Nature 495, 333–338. doi: 10.1038/nature11928

PubMed Abstract | CrossRef Full Text | Google Scholar

Nowak, R., and Poindron, P. (2006). From birth to colostrum: early steps leading to lamb survival. Reprod. Nutr. Dev. 46, 431–446. doi: 10.1051/rnd:2006023

CrossRef Full Text | Google Scholar

Qu, S., Yang, X., Li, X., Wang, J., Gao, Y., Shang, R., et al. (2015). Circular RNA: a new star of noncoding RNAs. Cancer Lett. 365, 141–148.

Google Scholar

Ragionieri, L., Cacchioli, A., Ravanetti, F., Botti, M., Ivanovska, A., Panu, R., et al. (2016). Effect of the supplementation with a blend containing short and medium chain fatty acid monoglycerides in milk replacer on rumen papillae development in weaning calves. Ann. Anat. 207, 97–108. doi: 10.1016/j.aanat.2016.04.035

PubMed Abstract | CrossRef Full Text | Google Scholar

Shannon, P., Markiel, A., Ozier, O., Baliga, N., Wang, J., Ramage, D., et al. (2003). Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 13, 2498–2504. doi: 10.1101/gr.1239303

PubMed Abstract | CrossRef Full Text | Google Scholar

Song, Y., Zhang, L., Liu, X., Niu, M., Cui, J., Che, S., et al. (2019). Analyses of circRNA profiling during the development from pre-receptive to receptive phases in the goat endometrium. J. Anim. Sci. Biotechnol. 10:34. doi: 10.1186/s40104-019-0339-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Takino, J. I., Sato, T., Nagamine, K., and Hori, T. (2019). The inhibition of Bax activation-induced apoptosis by RasGRP2 via R-Ras-PI3K-Akt signaling pathway in the endothelial cells. Sci. Rep. 9:16717. doi: 10.1038/s41598-019-53419-4

PubMed Abstract | CrossRef Full Text | Google Scholar

Trotta, R. J., and Swanson, K. C. (2021). Prenatal and postnatal nutrition influence pancreatic and intestinal carbohydrase activities of ruminants. Animals (Basel) 11:171. doi: 10.3390/ani11010171

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, J., Ren, Q., Hua, L., Chen, J., Zhang, J., Bai, H., et al. (2019). Comprehensive analysis of differentially expressed mRNA, lncRNA and circRNA and their ceRNA networks in the longissimus dorsi muscle of two different pig breeds. Int. J. Mol. Sci. 20:1107. doi: 10.3390/ijms20051107

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, W., Li, C., Li, F., Wang, X., Zhang, X., Liu, T., et al. (2016). Effects of early feeding on the host rumen transcriptome and bacterial diversity in lambs. Sci. Rep. 6:32479. doi: 10.1038/srep32479

PubMed Abstract | CrossRef Full Text | Google Scholar

Weary, D. M., Jasper, J., and Hotzel, M. J. (2008). Understanding weaning distress. Appl. Anim. Behav. Sci. 110, 24–41. doi: 10.1016/j.applanim.2007.03.025

CrossRef Full Text | Google Scholar

Xie, C., Mao, X. Z., Huang, J. J., Ding, Y., Wu, J. M., Dong, S., et al. (2011). KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Res. 39, W316–W322. doi: 10.1093/nar/gkr483

PubMed Abstract | CrossRef Full Text | Google Scholar

Xie, Y., Li, X., and Ge, J. (2019). Cyclophilin A-FoxO1 signaling pathway in endothelial cell apoptosis. Cell Signal. 61, 57–65. doi: 10.1016/j.cellsig.2019.04.014

PubMed Abstract | CrossRef Full Text | Google Scholar

Xue, M., Wang, K., Wang, A., Li, R., Wang, Y., Sun, S., et al. (2019). MicroRNA sequencing reveals the effect of different levels of non-fibrous carbohydrate/neutral detergent fiber on rumen development in calves. Animals (Basel) 9:496. doi: 10.3390/ani9080496

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Y., Guo, X., Pei, J., Chu, M., Ding, X., Wu, X., et al. (2020). CircRNA expression profile during yak adipocyte differentiation and screen potential circRNAs for adipocyte differentiation. Genes (Basel) 11:414. doi: 10.3390/genes11040414

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Z., Wang, C., Zhang, Y., Yu, S., Zhao, G., and Xu, J. (2020). CircDUSP16 promotes the tumorigenesis and invasion of gastric cancer by sponging miR-145-5p. Gastric Cancer 23, 437–448. doi: 10.1007/s10120-019-01018-7

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhong, T., Hu, J., Xiao, P., Zhan, S., Wang, L., Guo, J., et al. (2017). Identification and characterization of micrornas in the goat (Capra hircus) rumen during embryonic development. Front. Genet. 8:163. doi: 10.3389/fgene.2017.00163

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhong, T., Wang, C., Hu, J., Chen, X., Niu, L., Zhan, S., et al. (2020). Comparison of MicroRNA transcriptomes reveals the association between MiR-148a-3p expression and rumen development in goats. Animals (Basel) 10:1951. doi: 10.3390/ani10111951

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhong, Y., Du, Y., Yang, X., Mo, Y., Fan, C., Xiong, F., et al. (2018). Circular RNAs function as ceRNAs to regulate and control human cancer progression. Mol. Cancer 17:79. doi: 10.1186/s12943-018-0827-8

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhou, S., Liang, Y., Zhang, X., Liao, L., Yang, Y., Ouyang, W., et al. (2020). SHARPIN promotes melanoma progression via Rap1 signaling pathway. J. Invest. Dermatol. 140, 395–403 e6. doi: 10.1016/j.jid.2019.07.696

PubMed Abstract | CrossRef Full Text | Google Scholar