Integration of Non-Coding RNA and mRNA Profiles Reveals the Mechanisms of Rumen Development Induced by Different Types of Diet in Calves

Selecting suitable feed types and understanding the gastrointestinal digestive mechanism are helpful for the growth and health of calves in intensive dairy farming. However, the effects on rumen development of changing the molecular genetic basis and the regulatory mechanism by using different feed types are still unclear. Nine 7-day-old Holstein bull calves were randomly divided into GF (concentrate), GFF (alfalfa: oat grass = 3:2) and TMR (concentrate: alfalfa grass: oat grass: water = 0.30:0.12:0.08:0.50) diet experiment groups. Rumen tissue and serum samples were collected for physiological and transcriptomic analysis after 80 days. The results showed that serum α-amylase content and ceruloplasmin activity were significantly higher in the TMR group, and Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis ncRNAs and mRNAs were significantly enriched in the pathways of rumen epithelial development and stimulated rumen cell growth, including the Hippo signaling pathway, Wnt signaling pathway, thyroid hormone signaling pathway, ECM–receptor interaction and the absorption of protein and fat. The circRNAs/lncRNA-miRNAs-mRNA networks constructed, including novel_circ_0002471, novel_circ_0012104, TCONS_00946152, TCONS_00960915, bta-miR-11975, bta-miR-2890, PADI3 and CLEC6A, participated in metabolic pathways of lipid, immune system, oxidative stress and muscle development. In conclusion, the TMR diet could improve rumen digestive enzyme activities, stimulate rumen nutrient absorption and stimulate the DEGs related to energy homeostasis and microenvironment balance, and is thus better than the GF and GFF diets for promoting rumen growth and development.


Introduction
The dairy industry is an indispensable part of the national economy and is mainly based on natural grazing, semi-grazing and intensive farming. In the natural grazing and semi-grazing systems, due to the uncertainty of the ecosystem and environment, the DMI is low and ADF content is high in the diet, which reduces the efficiency of nutrient intake or absorption, and individual productivity cannot be fully developed [1]. Under foragebased livestock production systems, the quality and availability of feed determine the later growth and development of dairy cows [2]. Pasture is an important source of nutrition for calves. Selecting suitable feed types promotes the healthy growth of calves and reduces the feeding costs. Calves are fed solid feed after 7 days of life to increase productivity, fecundity and milk production, and to reduce the occurrence of disease [2,3]. The proportion and composition of solid feed can affect various aspects of development [4,5]. A high proportion of concentrated feed can improve ruminal starch digestibility, lower rumen pH and reduce the activity of the rumen papillae needed to absorb VFA [6][7][8]. The inclusion of dietary fiber through forage supplementation provides more physical stimulation, thereby enhancing the rumen volume, wall muscle content, motility and weight [9][10][11][12]. Therefore, with the rapid development of animal husbandry, determining whether different types of solid feed (including feed and concentrate) can provide certain nutrients for rumen growth and development is an urgent problem. Almost all studies to date on the effects of forage on rumen development have only been carried out at the level of nutrition, but the changes in the genetic basis and regulatory mechanism are still unclear and have never been studied.
Genomics takes all the nucleic acids of the organism as the research object. DNA, mRNA, lncRNA and other nucleic acid sequences are represented in the regulation of gene expression. The transcriptome is considered as the simplified research means for studying the genome, that is, the collection of all transcripts [13,14]. Therefore, RNA-seq is widely used to study the molecular mechanism of feed and to study the genes related to rumen growth and the development of ruminants, which play an important role in promoting rumen morphology and functional development [15,16]. miRNAs regulate the expression of target genes to maintain normal cellular homeostasis, digestion and nutrient absorption in the gastrointestinal tract [15,17]. lncRNA not only plays an important role in epigenetic inheritance, genomic imprinting and tissue cell formation, but also influences mRNA expression as an enhancer, promoter and transcriptional activator [18]. Studies have reported that being fed a high proportion of concentrate led to lncRNAs target genes activating PI3K-AKT, proteasome and HIF-1 signaling pathways, causing liver tissue inflammation and cell apoptosis, and aggravating SARA in dairy cows [15,19]. circRNA has a ring structure, is more stable than linear RNA and is rich in multiple binding sites of miRNA. It can regulate cell development, proliferation and differentiation by adsorbing miRNA and alleviating its inhibitory effect on mRNA in cells and tissues [20,21]. Zhang et al. showed that some candidate ceRNAs (cir_41, cir_115 and cir_171) might play key roles in immune and metabolic status in sheep [12].
Since the development of rumen in young animals is limited by their own physiological structure, the intake and type of solid feed are particularly important for nutrient absorption in rumen and the regulation of gene expression in epithelial cells to provide energy for themselves. In order to improve rumen growth and development, and the adult health and immunity of Holstein bulls before weaning, it is necessary to better understand the expression of ncRNAs and mRNAs in rumen tissues and the modifications involved in rumen growth-and development-related genes at the molecular level. At present, most studies on calves starter and abroad focus on nutrient levels and intestinal microorganisms, while there are few reports on the molecular genetic mechanism of rumen stimulation. Hence, the current study was designed to explore the serum immune factor content, digestive enzyme activities, rumen papilla size, ncRNAs and mRNAs regulatory mechanism, and their interactions under different feed types for the calf. The objective was to provide a comprehensive bioinformatics resource for selecting the most suitable feed type for rumen growth and development and the mechanism of ncRNA regulating rumen growth and development in calves.

Ethics Statement
This study was approved by and conducted in accordance with the ethical standards of the Institutional Animal Care and Use Committee of the College of Animal Science and Technology, Sichuan Agricultural University, Sichuan, 611130, China (20220236).

Animals and Feeding Strategy
The nine healthy Chinese Holstein bull calves (mean body weight, 41.58 ± 0.79 kg), 7 days old, were raised by the Mianyang Anzhou Hongfeng Dairy Cattle Breeding Co., Ltd. (Sichuan, China). They were randomly assigned to the GF, GFF and TMR groups, and each group was in a single pen. The GF, GFF and TMR groups were fed with "concentrate" (n = 3), "concentrate + hay" (n = 3), or a hay mixture (alfalfa: oat grass ratio = 3:2) and "concentrate: alfalfa grass: oat grass: water" = 0.30:0.12:0.08:0.50 (n = 3), respectively, with analysis after  (Table S1). Calves were given free access to both feed and water during the experimental time.

Detection of Rumen Papillae Height, Width and Number
All calves were slaughtered using electric-stunning slaughtering [22]. The back and abdomen tissue of rumen were collected. After fixation with 10% neutral formaldehyde, and after dehydration, pruning, embedding, sectioning, staining, sealing, etc., a Pannoramic 250 digital slice scanner produced by 3DHISTECH (Hungary) was used to collect images of the slices [23,24].

IgG Content and Digestive Enzyme Activity Detection
On the 80th day, before morning feeding, 5 mL of calf blood was collected by vacuuming the caudal vein with heparin sodium anticoagulant, and after standing for 30 min and centrifuging at 3000 r for 10 min to collect serum, it was stored at −80 • C in a refrigerator to determine serum indexes [15]. Serum samples were taken for IgG, α-amylase, ceruloplasmin, trypsin and lysozyme testing. An IgG ELISA kit was provided by Shanghai Zhuocai Biotechnology Co., LTD. The AMS Kit, CP Assay Kit, Trypsin Assay Kit and LZM Kit were all provided by Nanjing Jiancheng Institute of Bioengineering.

Total RNA Extraction from Calf Tissue
The collected rumen tissues were divided into 2 mL cryopreservation tubes and stored at −80 • C for further RNA extraction. The centrifuge was opened in advance and prechilled at 4 • C; after sterilization at high temperature, the mortar was wiped again with 75% alcohol, liquid nitrogen was prepared, the rumen tissue sample was ground into powder, an appropriate sample was put into a 1.5 mL EP tube, and 1 mL Trizol was added for total RNA extraction. The following steps were then carried out: (1) After adding 200 µL chloroform, it was shaken up and down for 2 min until the chloroform fully reacted with the rumen tissue lysate; then, after standing on ice for 5 min, it was centrifuged at 4 • C for 15 min at 12,000 r; (2) The supernatant of about 600 µL was sucked into a new 1.5 mL EP tube (the supernatant was sucked several times with a 200 µL pipette gun); (3) After adding 600 µL isopropyl alcohol at a ratio of 1:1, it was shaken up and down until the reaction was complete; then, after standing on ice for 10 min, it was centrifuged at 12,000 r and 4 • C for 10 min; (4) After quickly pouring out the supernatant and adding 1 mL of 75% ethanol prepared in advance, it was shaken upside down, and then centrifuged at 7500 r and 4 • C for 5 min; (5) After pouring out the supernatant, a 200 µL pipette gun was used to suck out the remaining supernatant at the bottom; (6) After instantaneous centrifugation, the remaining supernatant at the bottom was sucked out with a 100 µL pipette gun, and white precipitates could be seen; (7) According to the amount of white precipitation, DEPC water was added to dissolve the RNA, and a 20 µL pipette gun was used to blow and mix. RNA integrity was determined by 1.5% agarose gel electrophoresis, and RNA quality was determined using a nucleic acid protein analyzer and an Agilent 2100 bioanalyzer. When the bands of the glue running results were between 28S and 18S without dragging, this indicated that the integrity of the proposed RNA was good. When the concentration of the measured RNA was 1.8 ≥ OD ≥ 2.0, this indicated that the purity was good, which met the requirements of the subsequent test. The concentration and purity of RNA (OD260/280) were detected using a Nanodrop instrument. Part of the extracted RNA was used for full transcriptomic sequencing and part was used for differential gene expression RT-qPCR verification.

Transcriptomic Sequencing
The method of removing ribosomal RNA (circRNA library building to increase the process of removing linear RNA) was used to construct strain-specific libraries using Novogene [22]. Removing ribosomal RNA from total RNA, two strands were synthesized using RNA as a template. USER enzyme was used to degrade the second strand of cDNA containing U, and finally PCR amplification was performed to obtain the library [25]. First, Qubit was used for preliminary quantification, and the library was diluted to 1 ng/µL. Then, the Agilent 2100 BioAnalyzer was used to detect the insert size of the library [26]. The insert size was about 250-300 bp, in line with the expected Insert size. The effective concentration of the library was accurately quantified by qPCR. The effective concentration of the library was >2 nm and sequencing was performed using an Illumina PE150 [27].
2.6. Bioinformatics Analysis 2.6.1. Quality Control Trimmomatic Software was used to filter the Sequencing data to obtain clean readings [28]. The clean reads obtained by sequencing were compared to the reference genome using the comparison software, and the mapping rate was further analyzed [29].

ncRNA Identity
Clean reads obtained by sequencing were compared with reference genomes using HISAT2 software to identify lncRNA [30]. The CuffCompare software was used for comparison with the known database, and the known transcripts in the database were filtered out. Finally, the coding potential of the screened new transcripts was predicted, and Novel_lncRNAs and mRNAs were obtained [31]. For the identity of miRNAs, clean reads of each sample were screened, and sRNA with lengths ranging from 18 to 35 nt were selected for subsequent analysis [32]. Bowtie was used to locate sRNA after length screening to the reference sequence and analyze the distribution of small RNA on the reference sequence [33]. The mapped reads were aligned on the reference sequence with the specified range sequence in miRBase to get the sRNA details for each sample match. Then, we integrated miRNAs prediction software miREvo and miR-deep2 to analyze new miRNAs [34][35][36]. For circRNA identity, we first analyzed the numbers in the CIGAR file. We scanned the numbers in the Sam file and paired chiastic clipping signals. Then, the Junction reads were filtered based on PEM and GT-AG signals, and finally, the DM algorithm was used to detect the Junction reads, which further filters and prevents the false positives caused by homologous gene similarity and repeated sequences [37].

Quantification of Gene Expression Levels
StringTie was used to splice and quantify transcripts and genes using a network flow algorithm, and then the expression levels of all the known transcripts after alignment, splicing and screening, and the Novel_ncRNAs and mRNAs transcripts predicted, were quantified. Gene expression levels of mRNA and ncRNA were expressed as FPKM. FPKM (Fragments Per Kilobase of Transcript sequence Per Millions base pairs sequenced) is derived from a gene/transcript per kilobase length per million fragments count, which simultaneously corrects for the effect of sequencing depth and gene length on fragment counts, and is currently the most commonly used method for estimating gene expression levels [38].

Differential Expression Analysis
Differential expression analysis of the three groups was performed using the DESeq R package (1.10.1) [39]. Genes with an adjusted p-value found by DESeq were assigned as differentially expressed. Significantly differential expression genes were screened based on the following criteria: corrected p-value < 0.05 and |(fold change)| > 2 [40].

Function Annotations of RNA
GO enrichment analysis for host genes of differentially expressed ncRNAs was implemented by the GO-seq R package, in which gene length biases were corrected [41,42]. GO terms with corrected p-value < 0.05 were considered significantly enriched by differential gene expression. We used KOBAS software to test the statistical enrichment of differential expression genes or circRNA host genes in KEGG pathways [42,43]. The sequences of the DEGs were blast (blastx) to the genome of a related species (the protein interaction of which exists in the String data base: http://stringdb.org/ (accessed on 12 August 2022)) to get the predicted PPI of these DEGs. Then, the PPI of the DEGs was visualized in Cytoscape [15]. The interaction network was built and visually displayed using Cytoscape software based on the screening of cirRNA/lncRNA-miRNA-mRNA pairs [44].

Validation of RNA-Seq Results Analyzed by Quantitative PCR (qPCR)
To validate the repeatability of the RNA-seq analysis results, 18 candidate genes were randomly selected and evaluated by using qRT-PCR (for the primers, see Table S2). We performed qRT-PCR in a CFX96 Real-Time PCR Detection system (Bio-RadCo., Hercules, CA, USA) and detected RNA expression using SYBR Green Real-Time PCR Master Mix (Takara Co., Dalian, China). The relative expression levels were calculated using the 2 −∆∆Ct method and normalized to those of the reference genes ACTB and U6 [45,46].

The Data Statistical Analysis
IBM SPSS (27.0) software was used to analyze the experimental data by ANOVA. For quantitative variables, the data of the three groups were expressed as mean ± SD and compared using a t-test for statistical analysis. A value of p < 0.05 indicated a statistically significant difference [47].

Development of Rumen Dorsal and Ventral Papillae
The results showed that the three different diets greatly affected the growth and development of rumen dorsal and ventral papillae in calves (Table 1). In particular, there were significant differences in the width and number of papillae in calves in the posterior abdominal sac between the three groups (p < 0.01). In the GF group, the height and width of the middle dorsal and middle ventral papillae of calves were significantly greater than those in the GFF and TMR groups ((p < 0.01). Note: Values with superscript letters a , b and ab are significantly different across lines (p < 0.05). First, all the averages were arranged in descending order, and then the largest marked with the letter a . This average was compared with the other two sets of averages and marked with the letter a if the difference was not significant, and significant differences were marked with the letter b . The mean labeled b was taken as the standard and compared to the mean of the other group, without any letters if the difference was not significant, or with ab if it was in the same line as b .

Effects of IgG Content and Digestive Enzyme Activity
The serum content of IgG, α-amylase, ceruloplasmin, trypsin and lysozyme activity were measured ( Table 2). The results showed that the α-amylase content and ceruloplasmin activity in the TMR group were significantly higher than those in the GF and GFF groups (p < 0.01). There were no significant differences in trypsin activity, lysozyme activity and IgG in the GF, GFF and TMR groups (p < 0.05). Note: Values with superscript letters a and b are significantly different across lines (p < 0.05). First, all the averages were arranged in descending order, and then the largest marked with the letter a . This average was compared with the other two sets of averages and marked with the letter a if the difference was not significant, and significant differences were marked with the letter b .

Histopathological Examination (HE) for Pathological Change
The tissue section method was used to determine if there was pathological change in the rumen of calves fed on the different types of fodder. In this study, HE revealed in the middle dorsum of rumen neutrophils and a few lymphocytes in the GF group. In the GFF group, it was more severe and the local epithelial necrosis penetrated into the lamina propria, the cells in the necrotic area disintegrated, and the nuclei were fragmented, accompanied by a large number of infiltrated neutrophils and a small amount of bleeding. In the TMR group, the keratinized layer of the mucosal layer was obvious, the morphology of cells in the middle layer and basal layer was intact, and the lamina propria was mildly edemic (Figure 1). The results showed that lesions were more severe in the GFF group, followed by the GF group, and less severe in the TMR group.
Genes 2023, 14, x FOR PEER REVIEW 7 of 23 edemic ( Figure 1). The results showed that lesions were more severe in the GFF group, followed by the GF group, and less severe in the TMR group.

Overview of Identified ncRNA and mRNA
The Illumina HiSeq2500 sequencing platform was used to sequence 27 sample libraries, which generated a total of 72.18 GB of original data. After removing the reads with connectors and low quality in the original sequence, a total of 71.48 GB of clean reads were obtained, as shown in Table 3. The average base mismatch rate of the 27 sequenced libraries was 0.03%, the Q20 was higher than 97%, the Q30 was higher than 93%, and the GC

Overview of Identified ncRNA and mRNA
The Illumina HiSeq2500 sequencing platform was used to sequence 27 sample libraries, which generated a total of 72.18 GB of original data. After removing the reads with connectors and low quality in the original sequence, a total of 71.48 GB of clean reads were obtained, as shown in Table 3. The average base mismatch rate of the 27 sequenced libraries was 0.03%, the Q20 was higher than 97%, the Q30 was higher than 93%, and the GC content was 51.71%. The above sequencing quality control data showed that the experiment had high sequencing accuracy, a low base mismatch rate, reasonable distribution of base composition and good randomness in the sequencing process, which met the requirements of the later analysis. Then, we mapped the clean reads to the UMD 3.1 reference genome, and the alignment efficiency of ncRNA exceeded 90.63%. To identify ncRNA in calf rumen, transcripts were filtered by following a rigorous set of criteria. Finally, 50 miRNA, 1241 lncRNA and 1447 circRNA were obtained (Figure 2). To assess the distribution of the putative RNA, we analyzed the RNA expression ratio across the chromosomes. The putative RNAs were not equally distributed across the rumen chromosomes in Figure 3. Most mRNA were found on chromosome 1 and 7, miRNA were found on chromosome 5, whereas lncRNA were found on chromosome 7. In addition, most circRNA were found on chromosome 20 and X. To identify ncRNA in calf rumen, transcripts were filtered by following a rigorous set of criteria. Finally, 50 miRNA, 1241 lncRNA and 1447 circRNA were obtained ( Figure  2). To assess the distribution of the putative RNA, we analyzed the RNA expression ratio across the chromosomes. The putative RNAs were not equally distributed across the rumen chromosomes in Figure 3. Most mRNA were found on chromosome 1 and 7, miRNA were found on chromosome 5, whereas lncRNA were found on chromosome 7. In addition, most circRNA were found on chromosome 20 and X.

Expression Patterns of mRNAs and ncRNAs
To determine whether ncRNAs are involved in rumen growth and development, the DE ncRNAs and mRNAs from the different feed type groups were visualized using a histogram ( Figure 4

Expression Patterns of mRNAs and ncRNAs
To determine whether ncRNAs are involved in rumen growth and development, the DE ncRNAs and mRNAs from the different feed type groups were visualized using a histogram (Figure 4

Differentially Expressed ncRNAs and mRNA
The study identified the overlap between lncRNA target genes and DE mRNAs (Figure 5A-C). It was found that the relationships between XLOC_588071 and SLC51A, ARNTL2, CFAP299, XLOC_1002614 and C9, HNF4A, XLOC_609407, XLOC_994154 and

Differentially Expressed ncRNAs and mRNA
The study identified the overlap between lncRNA target genes and DE mRNAs ( Figure 5A-C). It was found that the relationships between XLOC_588071 and SLC51A, ARNTL2, CFAP299, XLOC_1002614 and C9, HNF4A, XLOC_609407, XLOC_994154 and ANXA13, XLOC_589694 and MYH11, LDB3 were co-expressed, which indicated that those lncRNAs are potential regulators.  Co-expression network analysis was performed on the selected differentially expressed circRNAs related to rumen growth and development and their targeted miRNAs. In the interaction network, there were four specific differentially expressed cirRNAs and their interacting miRNAs. According to the degree value, novel_circ_0019095, novel_circ_0011652, novel_circ_0016549 and bta-miR-211, novel_440, novel_circ_0012141, Co-expression network analysis was performed on the selected differentially expressed circRNAs related to rumen growth and development and their targeted miRNAs. In the interaction network, there were four specific differentially expressed cirRNAs and their interacting miRNAs. According to the degree value, novel_circ_0019095, novel_circ_0011652, novel_circ_0016549 and bta-miR-211, novel_440, novel_circ_0012141, novel_circ_0012117 and bta-miR-2300b-3p, novel_circ_0008114 and bta-miR-1, etc., played an important regulatory role in the network ( Figure 6 and Table 4). novel_circ_0012117 and bta-miR-2300b-3p, novel_circ_0008114 and bta-miR-1, etc., played an important regulatory role in the network ( Figure 6 and Table 4).

Differentially Expressed ncRNAs and mRNA Functional Enrichment Analysis
We further explored the potential regulatory roles of ncRNAs in rumen growth, which predicted the functions of ncRNAs using GO. The study shows that enrichment analysis based on the co-expression of lncRNAs ( Figure 7 and Table S3). The most significantly enriched in the three groups were immune response, immune system process, lymphocyte activation, extracellular region and extracellular space, hormone activity and Ras GTPase binding biological processes. The study revealed the most enriched pathways for mRNAs and circRNA/lncRNA-miRNA-mRNAs using p < 0.05. Interestingly, the ncRNA target genes exhibited KEGG enrichment patterns similar to those of the DE mRNAs ( Figure 8 and Table S4). The most enriched pathways were nutrient digestion and absorption (pathways in fat digestion and absorption, protein digestion and absorption, glycerophospholipid metabolism, α-linolenic acid metabolism, ABC transporter pathways and dorso-ventral axis formation) and signaling molecule and interaction pathways (NF-kappa B signaling pathway, intestinal immune network for IgA production, Wnt signaling pathway, Hippo signaling pathway, extracellular matrix (ECM)-receptor interaction and thyroid hormone signaling pathway) in the GF and TMR groups. In the GFF group, the most significant pathways were pathways in steroid hormone biosynthesis, cytokine-cytokine receptor interactions, systemic lupus erythematosus, alcoholism, DNA replication, base excision repair and The study revealed the most enriched pathways for mRNAs and circRNA/lncRNA-miRNA-mRNAs using p < 0.05. Interestingly, the ncRNA target genes exhibited KEGG enrichment patterns similar to those of the DE mRNAs ( Figure 8 and Table S4). The most enriched pathways were nutrient digestion and absorption (pathways in fat digestion and absorption, protein digestion and absorption, glycerophospholipid metabolism, αlinolenic acid metabolism, ABC transporter pathways and dorso-ventral axis formation) and signaling molecule and interaction pathways (NF-kappa B signaling pathway, intestinal immune network for IgA production, Wnt signaling pathway, Hippo signaling pathway, extracellular matrix (ECM)-receptor interaction and thyroid hormone signaling pathway) in the GF and TMR groups. In the GFF group, the most significant pathways were pathways in steroid hormone biosynthesis, cytokine-cytokine receptor interactions, systemic lupus erythematosus, alcoholism, DNA replication, base excision repair and arrhythmogenic right ventricular cardiomyopathy (ARVC). These genes are related to energy production and disease. arrhythmogenic right ventricular cardiomyopathy (ARVC). These genes are related to energy production and disease.

Construction of a Potential lncRNA/circRNA-miRNA-mRNA Regulatory Network
The results indicate that differentially expressed growth-and development-related genes in each group may participate in multiple different ceRNA networks [48]. More importantly, this study found that miRNAs may participate in different ceRNA networks in different rumen development nodes. For example, in the groups of GF vs. GFF, bta-miR-11973 involves a total of 3 mRNAs, 33 lncRNAs and 3 circRNAs in the network. However, in the GF vs. TMR groups, bta-miR-11973 may also be regulated by 14 lncRNAs and 2 circRNAs (Figure 9 and Table S5). The results present the regulatory relationships between ncRNAs and mRNAs in the processes of rumen growth and development.

Construction of a Potential lncRNA/circRNA-miRNA-mRNA Regulatory Network
The results indicate that differentially expressed growth-and development-related genes in each group may participate in multiple different ceRNA networks [48]. More importantly, this study found that miRNAs may participate in different ceRNA networks in different rumen development nodes. For example, in the groups of GF vs. GFF, bta-miR-11973 involves a total of 3 mRNAs, 33 lncRNAs and 3 circRNAs in the network. However, in the GF vs. TMR groups, bta-miR-11973 may also be regulated by 14 lncRNAs and 2 circRNAs (Figure 9 and Table S5). The results present the regulatory relationships between ncRNAs and mRNAs in the processes of rumen growth and development.

Validation of DE Genes by qRT-PCR
To ensure the accuracy and validity of the sequencing results, 18 genes (ncRNAs and mRNAs) were randomly selected and verified by qRT-PCR. Compared with the RNA-seq data, similar expression trends for qRT-PCR were observed, showing the strong consistency between the qRT-PCR and RNA-seq data ( Figure 10).

Validation of DE Genes by qRT-PCR
To ensure the accuracy and validity of the sequencing results, 18 genes (ncRNAs and mRNAs) were randomly selected and verified by qRT-PCR. Compared with the RNA-seq data, similar expression trends for qRT-PCR were observed, showing the strong consistency between the qRT-PCR and RNA-seq data ( Figure 10).

Discussion
The growth of the rumen papillae was the main manifestation of the physiological development of the rumen [49]. Most obviously, the height and width of the ruminal papillae of calves in the GF group were significantly greater than those in other groups. High concentrate diets contain a large amount of starch that could be converted into butyrate VFA by microorganisms to stimulate the growth of the rumen papillae, increasing the area of the papillae [50,51]. In the current study, the effective surface area of the rumen papilla fed high-concentrate diet was enlarged, which promoted nutrient absorption by the calves. Digestive enzyme activity plays an important role in the digestion of rumen nutrients, the regulation of pancreatic secretion and the promotion of pancreatic enzyme synthesis, and it is a good indicator for judging feed efficiency and the growth and development of calves [52]. In this experiment, we found that the TMR group diet increased the rumen α-amylase content and ceruloplasmin activity of calves, which were beneficial to the intake, digestion and utilization of solid feed, thus laying a good foundation for the intake of solid feed for calves [53]. Appropriate intake of concentrate could stimulate ruminal papillae growth and development. Likewise, a certain amount of roughage was also important for maintaining the degree of keratinization of the rumen epithelium [52]. Current results indicated that the TMR-type diet was more conducive to growth and development in calves.
The quality control and alignment efficiency results indicated that a successful library was constructed with high sequencing quality. In the study, 2935 mRNA, 50 miRNA, 1241 lncRNA and 3508 circRNA were obtained from the clean data. The results indicated that ncRNA play an important role in rumen growth and development at the transcriptional level. Furthermore, RT-qPCR was performed to confirm the reliability of sequencing data [54]. In order to understand the potential functions of the ncRNAs identified, functional enrichment analyses on the target genes of DE ncRNAs were performed. It was obvious that the functional GO enrichment of these DEGs was mainly related to the immune response. Fueling immune response is an energetically expensive process [55]. The transfer of nutrients from growth-to immune-related processes may increase maintenance requirements in animals during the immune response [12,56]. Therefore, the results of GO enrichment of these DEGs led to speculation that feeding calves with different feed types might result in less energy being consumed for combating systemic inflammation, and hence more efficient utilization of nutrients for growth and protein accretion.
The results of KEGG pathway analysis showed that the Hippo signaling pathway, Wnt signaling pathway, thyroid hormone signaling pathway and ECM-receptor interaction were significantly enriched in the TMR group. Those signaling molecule pathways are vital to cell differentiation, cell proliferation, regulated development and homeostasis, and rumen epithelial development [57][58][59][60]. In the present study, the TMR group might induce the differentiation of rumen cells and promote the maturation of rumen cells by signaling molecules and interaction. Systemic lupus erythematosus, DNA replication and alcoholism metabolic pathways were highly expressed in the GFF group in particular. Subsequently, we screened DEGs related to these pathways, such as HIST1H2BJ, HIST1H2AG, HIST1H2BB, HIST1H2BD and H4. Importantly, H4 binds to DNA to act in DNA replication and DNA repair by activating B-cell receptor signaling and autoantibodies. At the same time, the abnormal expression of DEGs was involved in the occurrence of cervical cancer [61]. It is worth mentioning that the Ras signaling pathway and Fc γ R-mediated phagocytosis single pathway were highly expressed in the GF and GFF groups. The literature reports that supplementation with high-concentrate diet, which activates the immune system, requires considerable energy, and the immune response might be an important factor contributing to the poorer feed efficiency and decreased performance of finished cattle [42,55,62,63]. Consistent with the results of this study, the literature indicates that feeding high-concentration concentrate to calves could trigger a series of inflammatory reactions, disorder of the immune system and glucose and lipid metabolism, which is not conducive to rumen growth and development. The above results indicated that the TMR group diet is more beneficial to rumen growth and development in calves.
Finally, to further understand how ncRNAs regulate rumen development, we constructed circRNAs/lncRNAs-miRNAs-mRNAs interaction networks based on competing ceRNAs to reveal the main functions of these ncRNAs. In the interaction network of the TMR group, some circRNAs containing several target sites of miRNAs (novel_circ_0002471, novel_circ_0012104 and bta-miR-11975) were found. Simultaneously, nine targeted key genes (PADI3, CLEC6A, FEV, MGLL, CACNG2, PLA2G2E, PLEKHS1, PSD and FHL3) were identified from five miRNAs (bta-miR-11975, bta-miR-1, bta-miR-2890, bta-miR-11974 bta-miR-211 and bta-miR-486), which participate in lipid, immune system, oxidative stress and muscle development and in the immune system. In addition, the highly expressed TCONS_00946152, TCONS_00960915 could act as ceRNAs to adsorb bta-miR-11975, and mainly participate in energy homeostasis and microenvironment balance by participating in the phosphatidylinositol signaling system, protein digestion and absorption, and arachidonic acid metabolism. The results might be related to the development of rumen, but the specific mechanism of action needs to be further investigated.
The TMR-type diet for calves in this experiment promoted enhanced rumen growth and development (as determined by α-amylase content and ceruloplasmin activity). In the rumen tissue from calves fed with the different type diets, the target genes of differentially expressed ncRNAs were significantly enriched in pathways closely related to rumen epithelial development and stimulated rumen cell growth, including the Hippo signaling pathway, Wnt signaling pathway, thyroid hormone signaling pathway and ECMreceptor interaction. Some target genes, including novel_circ_0002471, novel_circ_0012104, TCONS_00946152, TCONS_00960915, bta-miR-11975, bta-miR-2890, PADI3 and CLEC6A, play important roles in rumen development in calves. These target genes and pathways might be involved in the absorption and transport of protein and fat, further improving the growth performance of calves. The pathways and affected genes might serve as markers for further investigation into the mechanisms regulating growth and impacting rumen development, and provide favorable conditions for calf growth.

Conclusions
In this study, the GF group diet significantly increased the area of rumen papillae, the miRNAs and lncRNAs target genes were highly expressed in inflammatory-responserelated pathways, and single concentrate could cause a series of inflammatory reactions, which were not conducive to rumen growth and development. lncRNAs target genes in the GFF group were abnormally expressed during energy production and consumption, which reduced the feed utilization rate of calves, and led to immune system and lipid metabolism disorders in calves. The TMR diet was beneficial to maintain the keratinization degree of rumen epithelium, provide abundant conjugated linoleic acid, and induce rumen cell differentiation and promote rumen cell maturation through signaling molecule interaction. Therefore, in order to improve the rumen growth and development of calves and the longterm stable development of animal husbandry, the TMR feeding strategy is recommended. There are some limitations to this study. Considering animal welfare and the 3Rs principle, the number of calf samples used in this study was limited, which made the biological repetition of this experiment lower, but it still provides a comprehensive information resource for future research on related aspects. Therefore, in vitro experiments are needed to further improve relevant studies.

Conflicts of Interest:
The authors declare no conflict of interest.