Effects of Gut Microbiome and Short-Chain Fatty Acids (SCFAs) on Finishing Weight of Meat Rabbits

Understanding how the gut microbiome and short-chain fatty acids (SCFAs) affect finishing weight is beneficial to improve meat production in the meat rabbit industry. In this study, we identified 15 OTUs and 23 microbial species associated with finishing weight using 16S rRNA gene and metagenomic sequencing analysis, respectively. Among these, butyrate-producing bacteria of the family Ruminococcaceae were positively associated with finishing weight, whereas the microbial taxa related to intestinal damage and inflammation showed opposite effects. Furthermore, interactions of these microbial taxa were firstly found to be associated with finishing weight. Gut microbial functional capacity analysis revealed that CAZymes, such as galactosidase, xylanase, and glucosidase, could significantly affect finishing weight, given their roles in regulating nutrient digestibility. GOs related to the metabolism of several carbohydrates and amino acids also showed important effects on finishing weight. Additionally, both KOs and KEGG pathways related to the membrane transportation system and involved in aminoacyl-tRNA biosynthesis and butanoate metabolism could act as key factors in modulating finishing weight. Importantly, gut microbiome explained nearly 11% of the variation in finishing weight, and our findings revealed that a subset of metagenomic species could act as predictors of finishing weight. SCFAs levels, especially butyrate level, had critical impacts on finishing weight, and several finishing weight-associated species were potentially contributed to the shift in butyrate level. Thus, our results should give deep insights into how gut microbiome and SCFAs influence finishing weight of meat rabbits and provide essential knowledge for improving finishing weight by manipulating gut microbiome.


INTRODUCTION
The gut microbiome comprises the collective genome of the trillions of microorganisms colonizing the intestinal tract, which is responsible for vital metabolic, immunological, and nutritional functions (Sidhu and van der Poorten, 2017). Short-chain fatty acids (SCFAs), mainly acetate, propionate, and butyrate, are recognized as important metabolites derived from the fermentation of indigestible dietary fibers by gut microbial species, which play fundamental roles in maintaining intestinal homeostasis and energy metabolism regulation (Priyadarshini et al., 2018).
Recently, the roles of gut microbiome and SCFAs in modulating growth performances of farm animals have been extensively studied. For instance,  indicated that certain bacteria of the family Ruminococcaceae, the genus Faecalibacterium, and the genus Prevotella are involved in metabolizing dietary fibers which subsequently affect SCFA production, which exerted important roles in feed intake and body weight modulation in pigs. Du et al. (2018) demonstrated that both butyrate-producing bacteria Butyricimonas and the digestive enzyme family metabolic pathway could potentially affect the body weight of beef calves. Additionally, research on the gut microbiome of broiler chickens suggested that increased abundance of Lactobacillus and Bifidobacterium, and facilitated generations of lactate and SCFAs, had a positive influence on feed efficiency and body weight gain (Yu et al., 2019).
Finishing weight is regarded as a crucial growth trait of meat rabbits, which directly affects broiler rabbits' meat production. However, only few studies have unraveled the relationship between the gut microbiome and finishing weight of rabbits. Zeng et al. reported that changes in the abundance of YS2, Bacteroides, Lactococcus, Lactobacillus, and Prevotella in the gut microbial communities of rabbits were associated with finishing weight (Zeng et al., 2015). Wang Q. et al. (2019) found that Coprococcus and Ruminococcaceae_UCG-004 may hinder the production of pro-inflammatory factors that exert beneficial effects on the finishing weight of rabbits. Moreover, North et al. (2019) indicated that rabbits with greater abundances of Eubacteriaceae, Natranaerobiaceae, Peptococcaceae, and Syntrophomonadaceae in the gut microbial communities tended to gain more weight at finishing. Nevertheless, these investigations are based on 16S rRNA gene sequencing analysis, which cannot determine species and functional capacities that affect finishing weight. Additionally, the potential roles of SCFAs in modulating the finishing weight of meat rabbits remain unknown. Hence, the present study aimed to assess the effects of gut microbiome and SCFAs on finishing weight using 16S rRNA gene sequencing, metagenomic sequencing, and fecal SCFA level data. We not only identified a number of microbial taxa and metabolic functions associated with finishing weight but also unraveled the correlations between SCFA levels and finishing weight. In addition, we evaluated the role of the gut microbiome in the phenotypic variation of finishing weight. Such findings could improve our understanding of how the gut microbiome and SCFAs modulate finishing weight and thus provide important information for the meat rabbit industry.

Experimental Rabbits and Sample Collection
One hundred and five Ira rabbits (28 ± 2 days, 53 males and 52 females) with similar weaning weight were randomly distributed to separated cages (one rabbit per cage) and fed with the fatten diet (details are shown in Supplementary Table S1) until finishing (72 ± 2 days). Finishing body weight was measured, and hard fecal samples were collected. Finishing weights were sorted, and the top 5 rabbits with the highest body weight (high group) and the bottom 5 with the lowest body weight (low group) were selected for metagenomic sequencing. All rabbits were healthy and had not received antibiotics, anti-coccidial drugs, probiotics, or prebiotics during the experimental period. All fecal samples were snap frozen in liquid nitrogen for transportation and stored at −80 • C until further processing.

16S rRNA Gene Sequencing
Microbial genomic DNA of fecal samples was extracted by the QIAamp Fast DNA Stool Mini Kit (QIAGEN, Germany) according to the manufacturer's instructions. The purity and integrity of total DNA was detected by using the NanoDrop ND-2000 spectrophotometer (Thermo Fisher Scientific, United States) and 1.5% agarose gel electrophoresis, respectively. Hypervariable regions V3-V4 of the 16S ribosomal RNA gene were then amplified by using the primer 341F:CCTACGGGNGGCWGCAG and 806R:GGACTACHVGGGTATCTAAT and sequenced by the HiSeq 2500 platform (Illumina, United States) according to the manufacturer's manuals. QIIME (v.1.9.1) was used for the quality control process of sequencing data including filtering out of the primers, barcodes, and low-quality sequences (quality score < 20) (Caporaso et al., 2010). High-quality paired-end reads were merged into tags by using FLASH (v.1.2.11) (Magoc and Salzberg, 2011). We rarefied the library size of microbial sequences to 40,000 tags per sample to normalize the sequencing depth (Fu et al., 2015). USEARCH (v.10.0) was used to cluster tags into operational taxonomic units (OTUs) at 97% sequence similarity (Edgar, 2010). OTUs which had relative abundance < 0.1% and were presented in less than 1% of the experimental rabbits were filtered out before further analysis. OTU taxonomic category assignments were performed by using the SILVA database (v.132) (Quast et al., 2013).

Two-Part Model Analysis
After sex and cage rearing effect corrections, residuals of finishing weight phenotypic values were used for association analysis between finishing weight and relative abundances of OTUs by using a two-part model method as described before (Fu et al., 2015). In brief, the two-part model analysis was consisted of binary, quantitative, and meta models. The binary model describes a binomial analysis that assesses for the effect of the presence/absence of the gut microbe on finishing weight. The quantitative model tests for the association between finishing weight and the abundance of each OTU, but only the samples where that microbe is present were included in the analysis. The meta model was used to integrate the effect of both binary and quantitative analysis. The final association p-value was assigned from the minimum of p-values from the binary analysis, quantitative analysis, and meta-analysis. Skewness correction was performed by 1000 × permutation tests. False discovery rate (FDR)-adjusted p < 0.05 was set as the significance threshold.

Metagenomic Sequencing
A paired-end (PE) DNA library was constructed for each sample according to the manufacturer's instruction (Illumina, United States), and sequencing was performed on an Illumina HiSeq 4000 platform. Quality control, adapter trimming, and low-quality read filtering of raw data were performed by using fastp (v.0.19.4) (Chen et al., 2018) and obtained an average of 64,042,801 high-quality reads for each sample. We assembled the high-quality reads into contigs using the MEGAHIT (v.1.1.3) . The contigs with more than 200 bp in length were used for open reading frame (ORF) prediction with MetaGeneMark (v.2.10) (Zhu et al., 2010). A non-redundant gene catalog was established by removing the redundant genes from the predicted ORFs using Cd-hit (v.4.6.1) (Fu et al., 2012). The gene abundance profile was generated by mapping the high-quality reads against the non-redundant gene catalog using MOCAT (v2.0) (Kultima et al., 2016). DIAMOND (v.0.9.24) (Buchfink et al., 2015) was used for taxonomic category and GO term assignments of the predicted genes by aligning against the non-redundant (NR) database and Gene Ontology (GO) database, respectively. Carbohydrate-Active enZYmes (CAZymes) were annotated by hmmscan program in HMMER (v.3.1) (Eddy, 2011). GhostKOALA was used to extract KEGG Orthology (KO) and KEGG pathway annotation information from the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (Kanehisa et al., 2016).

Co-abundance Group (CAG) Analysis
In order to construct CAG, we first calculated the correlation coefficients among the finishing weight-associated microbial taxa using the Sparse Correlations for Compositional data (SparCC) algorithm (Friedman and Alm, 2012). Then, CAGs were defined by Ward's linkage hierarchical clustering method and PERMANOVA (999 permutations, p < 0.05) based on the SparCC correlation coefficient matrix through made4 and vegan R packages, respectively (Zhang et al., 2016). The CAG interaction networks were established based on FDR-adjusted p < 0.05 and | r| > 0.35 by CYTOSCAPE (v.3.7.0) (Shannon et al., 2003). Spearman's correlation analysis with FDR correction was performed to compute the association between CAGs and finishing weight.

Estimating Phenotypic Variance Explained by the Gut Microbiome
To evaluate the contribution of the gut microbiome to the variation of finishing weight, we performed a 100× cross-validation (Fu et al., 2015). The data set was randomly divided into an 80% discovery data set and a 20% validation data set. In the discovery data set, we performed two-part model analysis to identify a number of (n) OTUs that were significantly associated with finishing weight at a certain p-value and assessed the effect sizes of binary and quantitative features (β 1 and β 2 ) of each OTUs. In the validation data set, the effect of the gut microbiome on finishing weight (r m ) was estimated by an additive model: r m = n j=1 (β 1 + b j + β 2 q j ), b j and q j corresponding to the binary and quantitative feature of j OTU, respectively. We calculated the squared correlation coefficient (R 2 ) between r m and the phenotypic value, which represents d or the phenotypic variance explained by the gut microbiome. To ensure validity and stability of the estimation, we repeated the cross-validation for 100 times and calculated the average value of the explained variations.

Determinations of SCFA in Fecal Samples
The concentrations of SCFAs were measured in the fecal samples using gas chromatography on an Agilent 7890A GC system with a flame ionization detector (Agilent Technologies, United States) according to the previously described method with some modifications (Reddivari et al., 2017). Hundred and fifty milligram fecal samples were homogenized in 1 ml of water -0.5% phosphoric acid by Vortex-Genie 2 (Scientific Industries, United States). The mixture was centrifuged at 15,000 rpm for 10 min at 4 • C. The supernatant (600 µl) which was mixed with an equal volume of ethyl acetate was used for filtering by a 0.35 mm filter membrane and then transferred to a glass gas chromatography vial prior to injection into a GC instrument. The internal standard of 2-methyl butyrate at a final concentration of 10 mM was added to the organic phase to correct injection variability among the samples. 10-100 mM acetate, 1-10 mM propionate, and 1-10 mM butyrate were set as external standards, and five sample intervals were prepared to calculate retention times and create standard curves. The concentrations of SCFAs in fecal samples were determined by comparing to standard curves. The correlations between finishing weight-associated species and SCFA levels were analyzed by using the Spearman method.

Basic Statistical Analysis
The Wilcoxon rank-sum test with FDR correction was used to identify the metagenomic species and function capacities of the gut microbiome having significantly different enrichments between high and low finishing rabbits. To identify important species that could predict finishing weight variation, random forest analysis was performed by using species relative abundance data, and a predictor species was determined based on the mean decrease in accuracy of discrimination > 2 (Breiman, 2001). Two-sided unpaired Student's t-test with FDR correction was performed to detect the differences in levels of SCFAs in fecal samples between high and low finishing weight rabbits. Except the differential metagenomic species which were visualized online by Interactive Tree Of Life (iTOL) (Letunic and Bork, 2019), the other plots were accomplished in R software.

Gut Microbiota Associated With Finishing Weight
The phylogenetic analysis of gut microbiota in experimental rabbits has been previously described (Fang et al., , 2020. Here, we focused on identifying the gut microbiota associated with finishing weight. The phenotypic distribution of finishing weight in experimental rabbits is shown in Supplementary Figure S1A. The phenotypic value of rabbits in the high finishing weight group was significantly different than that in the low finishing weight group (FDR adjusted P < 0.05, Supplementary Figure S1B).
To identify microbial taxa associated with finishing weight, we performed a two-part model analysis using the relative abundances of OTUs and phenotypic values in the experimental population. We identified 15 OTUs associated with finishing weight (Figure 1 and Supplementary Table S2, FDR adjusted P < 0.05), of which 8 OTUs exhibited positive associations and the remaining showed negative associations. Among the OTUs positively associated with finishing weight, one was annotated FIGURE 1 | The OTUs significantly associated with finishing weight (FDR adjusted P < 0.05) are shown as Z scores.
to the order Mollicutes_RF9. Two were annotated to the family Ruminococcaceae, and one to the family Lachnospiraceae. Three OTUs were annotated to the genus Ruminococcaceae_UCG-014, and one to the genus Ruminococcus_1. As for the OTUs negatively associated with finishing weight, three OTUs were annotated to the family level (two belonging to the family Coriobacteriaceae and one to the family Alcaligenaceae) and four to the genus level (Ruminococcaceae_UCG-005, Tyzzerella_3, Ruminiclostridium, and Anaerotruncus one OTU each).

Interactions of Microbial Taxa Correlated With Finishing Weight
Previous studies have demonstrated that the interactions of gut microbes play vital roles in host physiology, development, and growth (Ramayo-Caldas et al., 2016;Ke et al., 2019). To investigate the potential effect of interactions among the identified species on finishing weight, we constructed a co-abundance network based on SparCC correlation coefficients. The co-abundance network was consisted of two co-abundance groups (CAGs; Figure 3A and Supplementary Figure S2). The species that were positively associated with finishing weight were clustered into CAG 1 through intra-positive interactions, whereas CAG 2 was formed by intra-positive interactions among the species that were negatively associated with finishing weight. In addition, CAG 1 had positive associations with finishing weight, but CAG 2 was negatively correlated with finishing weight (Figures 3B,C). Notably, when we analyzed the entire experimental population, the co-abundance network could also be established using finishing weight associated OTUs, and similar correlations FIGURE 2 | The phylogenetic relationships and abundances of metagenomic species enriched in high and low finishing weight rabbits. The phylogenetic tree of species is shown as the innermost layers. Labels with red and blue color represent high and low finishing weight groups, respectively. The different color strips in the third layer correspond to different families as indicated by the color code on the left. The outermost layers depict the relative abundances of differential species between high and low finishing weight individuals.
were observed between the CAGs and finishing weight (Supplementary Figure S3).

Functionalities of the Gut Microbiome Related to Finishing Weight
Functionalities of the gut microbiome related to finishing weight were detected by comparing the abundances of CAZymes, GOs, and KEGG items between rabbits with high and low finishing weights.
A total of 30 CAZymes with significantly different abundances in rabbits with high and low finishing weights were identified (Figure 4A and Supplementary Table S4, FDR adjusted P < 0.05). Among these, 16 CAZymes especially galactosidase and xylanase (e.g., GH27, GH95, GH4, GH51, GH120, and GH11) were significantly enriched in the gut microbiome of rabbits with high finishing weight. 14 CAZymes, mainly glucosidase (e.g., GH30, GH1, and GH3) had significantly higher abundances in the gut microbiome of low finishing weight rabbits.
On the other hand, 25 differential enriched KEGG pathways were identified in the gut microbiome of rabbits with distinct finishing weights (Figure 5B and Supplementary Table S4, FDR adjusted P < 0.05). Similarly, we found that aminoacyl-tRNA biosynthesis, glycolysis/gluconeogenesis, butanoate metabolism, and cysteine and methionine metabolism were more active in the gut microbial communities of high finishing weight rabbits. Meanwhile, in rabbits with low finishing weight, gut microbiota were more capable of operating ABC transporters and two-component system; additionally, the metabolism of tyrosine, fructose, and mannose was more active in these animals.

Gut Microbiome as a Predictor of Finishing Weight Variation
To assess whether the gut microbiome was able to predict finishing weight, we investigated how much degree of phenotypic variance of finishing weight was explained by the gut microbiome by performing 100 times cross-validation analyses at different p-value thresholds (ranging from 10 −5 to 0.1) using the OTU data. The OTUs identified at p = 1 × 10 −5 could explain 8.01% of the variations in finishing weight ( Figure 6A); at p = 0.1, the explained variation increased to 10.85%, given that more OTUs were included in the analysis as the threshold increased.
Then, to evaluate whether a subset of metagenomic species within the gut microbial community could predict the finishing weights of rabbits, we performed random forest analysis in high and low finishing weight rabbits. Thirty-one species were identified as being able to predict finishing weight and showed a substantial overlap with the species above identified as being associated with finishing weight (Figure 6B).

Changes in Fecal SCFAs Linked to Finishing Weight
As SCFAs produced by gut microbiota exert important effects on host energy metabolism and intestinal health status regulation in animals (den Besten et al., 2013), we analyzed the levels of acetic acid, propionic acid, and butyric acid in the feces of rabbits with high and low finishing weights ( Figure 7A). Our results showed  that the level of butyric acid was significantly greater in the feces of rabbits with high finishing weight in comparison to that in rabbits with low finishing weight (FDR adjusted P < 0.05), but the levels of propionic acid tended to decrease. No significant change in the levels of acetic acid was observed between high and low finishing weight rabbits. In addition, correlation analysis indicated that the level of butyric acid had a significantly positive association with finishing weight, while the levels of acetic acid and propionic acid showed the tendency of positive and negative association with finishing weight, respectively (Figures 7B-D).
We also analyzed the correlations between finishing weight-associated species and SCFAs levels. As shown in Figure 8, nine out of ten species enriched in the high finishing weight individuals were positively associated with the level of butyric acid, and nine out of thirteen species augmented in the low finishing weight individuals were negatively correlated with the level of butyric acid. Both Anaerotruncus sp. CAG:390 and Bacteroides thetaiotaomicron were more abundant in low finishing weight rabbits and had positive associations with the level of propionic acid. However, no significant correlations were observed between the level of acetic acid and finishing weight-associated species.

DISCUSSION
Accumulating evidence indicates that the gut microbiome plays a pivotal role in nutrient metabolism, energy utilization, and health maintenance in domestic animals (O'Callaghan et al., 2016). Accordingly, the gut microbiome is considered to be an essential factor that affects growth of livestock animals. Thus, identifying gut microbiota as key for production performances could be a game changer in the livestock production industry (Maltecca et al., 2020). SCFAs are important metabolites generated by gut microbial fermentation of complex carbohydrates, which have emerged as key regulators in intestinal health and energy homeostasis regulation. However, the effects of gut microbiome and SCFAs on finishing weight of meat rabbits have received far less attention. Here, we systematically and comprehensively evaluated how the gut microbiome and SCFAs affect the finishing weight of meat rabbits.
We identified 15 OTUs that were significantly associated with finishing weight, of which eight OTUs were annotated to members of the family Ruminococcaceae (Figure 1). Among these, Ruminococcaceae_UCG-014 and Ruminococcus_1 were positively associated with finishing weight, while Ruminococcaceae_UCG-005 and Ruminiclostridium showed the opposite correlations. Ruminococcaceae_UCG-014 belongs to butyrate-producing bacteria capable of degrading cellulose and hemicellulose in feeds and has a potential role in maintenance of intestinal health . Ruminococcus_1 is another important player in butyrate production by fermenting complex non-digestible polysaccharides, which is considered to be related to intestinal anti-inflammatory responses (Xie et al., 2019). In contrast, Ruminococcaceae_UCG-005 and Ruminiclostridium were positively correlated with diarrhea incidence and intestinal villi damage, respectively (Hung et al., 2019;Xing et al., 2019).
Family Ruminococcaceae encompasses many species. It is important to know which are associated with a good intestinal health status. Metagenomic sequencing analysis sheds light on this. Our results showed that different species from family Ruminococcaceae exert contrasting effects on finishing weight (Figure 2). For example, Ruminococcus albus and Ruminococcus flavefaciens were enriched in rabbits with high finishing weight, but Ruminococcus gauvreauii and Ruminococcus bromii were more abundant in low finishing weight individuals. Both R. albus and R. flavefaciens are well-known cellulolytic bacteria involved in the butyrate metabolic pathway and abundant in ruminants with high body weight gain (Carberry et al., 2012;Mao et al., 2017;Izuddin et al., 2019). However, both R. gauvreauii and R. bromii are mucolytic bacteria which can induce chronic intestinal inflammation (Lyra et al., 2009;FIGURE 7 | Fecal short-chain fatty acid levels in high and low finishing weight rabbits and their correlations with finishing weight. (A) Comparisons of fecal short-chain fatty acid levels between high and low finishing weight rabbits. (B-D) Correlations between fecal short-chain fatty acid levels and finishing weight. Crost et al., 2018;Fernandez et al., 2018). Moreover, R. bromii is negatively associated with body weight in pigs (Tran et al., 2018).
We identified several other species that could significantly affect finishing weight. For instance, Faecalibacterium prausnitzii is a major butyrate producer in the intestine, which has an important role in providing energy sources for enterocytes and fulfilling anti-inflammatory actions (Foditsch et al., 2014). Lactobacillus ruminis is a member of lactic acid bacteria that has multiple probiotic properties, including inhibition of intestinal pathogens, fortification of epithelial barrier functions, and modulation of immune responses (Yu et al., 2017). In agreement with previous studies (Oikonomou et al., 2013;Gao et al., 2017), we found that these two species were more abundant in the gut microbiome of individuals with a higher body weight. In contrast, although both Akkermansia muciniphila and Bacteroides fragilis were considered to be promising candidates as probiotics (Gilbert et al., 2013;Zhang T. et al., 2019), the overgrowth of these two species in the gut of rabbits could result in the decreased butyrate yield and lead to the incidence of epizootic rabbit enteropathy (ERE, a severe gastrointestinal syndrome disease) (Jin et al., 2018). Hence, these two species augmenting in the gut microbiome of meat rabbits could be a potential trigger for reduction of finishing weight.
In order to make a microbial community-based inference, it is important not only to identify microbial taxa associated with finishing weight but also to understand the effects of the interactions among such taxa on finishing weight. Thus, co-abundance networks were established at both metagenomic species and OTU level. We found that co-occurrence and co-exclusion relationships exist among finishing weight-associated microbial taxa (Figure 3 and Supplementary Figure S3). Importantly, their interactions which could affect the finishing weight have intimate relations with their distinct roles in host metabolic and intestinal health regulation. Similarly, earlier studies have revealed that interactions of gut microbes affect the body weight in farm animals. Yang et al. demonstrated that the strong co-occurrence relationship between butyrate-producing bacteria (e.g., Ruminococcaceae) and lactic acid bacteria (e.g., Lactobacillus) may affect porcine body weight gain by regulating host appetite and feeding behavior (Yang et al., 2018). Zhang et al. (2018) suggested that the significant co-exclusion relationship between Enterobacteriaceae and Prevotellaceae was associated with the metabolism of SCFAs and energy, which provides an important direction for manipulating gut microbiota to improve the body weight in pigs. Additionally, Gao et al. (2017) and Ma et al. (2018) highlighted the important role of gut microbial interactions in improving growth performances of broiler chickens.
Distinct metagenomic functional capacities were also found to be associated with finishing weight. The CAZyme comparison analysis showed that galactosidase and xylanase were more abundant in rabbits with high finishing weight, while glucosidase was more abundant in low finishing weight rabbits ( Figure 4A). Adding extra galactosidase and xylanase into feeds has been reported to enhance nutrient (e.g., fibers and amino acids) digestibility and improve the body weight of farm animals Jasek et al., 2018). This could be used to explain why rabbits with a greater abundance of these two enzymes showed greater finishing weights. However, the effects of gut microbial glucosidase on body weight are conflicting. De Cesare et al. (2017) demonstrated that a higher abundance of glucosidase in the gut microbiome contributed to increased body weight gain of broiler chickens, whereas Shokryazdan et al. (2017) suggested that the reduced glucosidase enzyme activity in the intestine led to an increase in body weight in broiler chickens.
Our results also showed that GOs related to xylan, galactose, and arabinose metabolism were more abundant in rabbits with high finishing weight, while those associated with glucose metabolism showed higher abundances in low finishing weight rabbits ( Figure 4B). Additionally, we found that the gut microbiota of rabbits with high finishing weight preferred to be involved in the biosynthesis of asparagine, cysteine, and arginine, but gut microbial communities in low finishing weight rabbits were more active in metabolizing tryptophan, tyrosine, and phenylalanine. Although asparagine is a non-essential amino acid, emerging evidence has demonstrated that it plays a key role in attenuating intestinal injury and improving the energy status of enterocytes, all of which are good for the health and growth of animals Patra et al., 2019). Cysteine is able to regulate the antioxidant status and expression of anti-inflammatory genes in intestinal cells, which exert beneficial effects on body weight gain in pigs and chickens (Dilger and Baker, 2007;Yang and Liao, 2019). Arginine is a hub precursor for the synthesis of various important metabolic molecules, including NO and polyamines, exerting regulatory roles in the host's metabolic processes (Flynn et al., 2002). Thus, it has been used to improve meat production in meat-producing animals (Hu et al., 2017;Liu et al., 2019). Conversely, gut microbiota is involved in aromatic amino acids (tryptophan, tyrosine, and phenylalanine) metabolism has been linked to chronic low-grade inflammation and metabolic disorders, which may hinder development and growth of animals (Wang et al., 2011;Agus et al., 2018).
On the other hand, the aminoacyl-tRNA biosynthesis pathway and the metabolic pathways of butanoate, cysteine, and methionine were more active in high finishing weight rabbits, whereas ABC transporters, and fructose, mannose, and tyrosine metabolic pathways were overrepresented in rabbits with low finishing weight ( Figure 5). The aminoacyl-tRNA biosynthesis pathway is an essential metabolic function of microorganisms and is intimately correlated with maturation of gastrointestinal microbiota (Zhang K. et al., 2019). More importantly, microbial mature status is fundamental for host development and growth . Several studies have demonstrated that the gut microbial butanoate metabolic pathway exerts beneficial effects on body weight of livestock due to its potential impact on increasing energy intake, improving intestinal histomorphometric characteristics (e.g., villi height and crypt depth), and modulating the immune system (Ndou et al., 2018;Che et al., 2019;Jacquier et al., 2019). Similar to the abovementioned cysteine, methionine is crucial in maintaining the functions of intestinal cells by regulating the redox status (Azad et al., 2019). Thus, enhanced methionine metabolism in the gut microbial communities contributed to bovine body weight gain (Jiao et al., 2017). On the other hand, ABC transporters involved in transporting a large variety of sterols, lipids, drugs, and primary and secondary metabolites (Hou et al., 2017) have recently been related to fat deposition, drug resistance, and colonic inflammation (Andersen et al., 2015;Fang et al., 2017), which may do harm for growth performances in animals. Furthermore, enhanced fructose, mannose, and tyrosine metabolic pathways in the gut microbiome have been associated with host intestinal permeability and metabolic endotoxemia that implies their negative effects on growth of animals (Wang et al., 2011;Do et al., 2018).
To further emphasize the predictive role of the gut microbiome in finishing weight variation, we estimated that the gut microbiome could explain 7.23-10.85% of the variation in finishing weight (Figure 6A), which is similar to host genetics on finishing weight (6.3-12.1%) (Piles et al., 2007). Moreover, we identified 31 species that could act as predictors of finishing weight, most of which are strongly associated with finishing weight (Figure 6B). Consistent with prior studies, our findings suggested that the gut microbiome should be regarded as a key variable of body weight prediction model of farm animals Wen et al., 2019).
Our results also revealed that the level of fecal butyric acid was significantly higher in rabbits with high finishing weight in comparison to rabbits with low finishing weight ( Figure 7A). Importantly, the level of butyric acid showed positive associations with finishing weight (Figure 7B). Butyrate is not only recognized as an essential energy source but also acts as a signal transduction molecule of G-protein-coupled receptors (FFAR3, GPR109A) and as epigenetic regulators of gene expression by inhibiting histone deacetylase (HDAC) (Kasubuchi et al., 2015). Due to these important properties, butyrate exhibits wide energy and metabolic regulation abilities and strong anti-inflammatory effects that positively affect the body weight gain of animals (Bedford and Gong, 2018). Consistently, we found that several microbial species positively associated with finishing weight (such as, Faecalibacterium prausnitzii, Roseburia sp. CAG:303, and butyrate-producing bacterium SS3/4) were well-known butyrate producers (Qin et al., 2012) that had positive associations with butyric acid level (Figure 8). Interestingly, we also observed that both Lactobacillus ruminis and Bifidobacterium saeculare were positively associated with butyric acid level. This is in accordance with previous findings which revealed that high levels of intestinal Lactobacillus sp. and Bifidobacterium sp. contributed to increased butyrate levels (Le Roy et al., 2015;Kim et al., 2020). In addition, our results showed the changes in the levels of acetic acid and propionic acid in the feces of rabbits with high and low finishing weight and their potential correlations with finishing weight (Figures 7A,C,D). These findings were in agreement with previous studies suggested that both acetate and propionate were important players in modulating growth performances in animals Min et al., 2019).
Our study, although limited by a small rabbit population, provides important information that supports the potential of gut microbiota manipulation in improving finishing weight of meat rabbits. However, from potential to action, the next steps are to validate the causative roles of some gut bacterial species in modulating finishing weight through specific pathogenfree (SPF) rabbits' intervention experiment and to gain more mechanistic insights into the cross talk between gut microbiome and host and its effects on growth performances of meat rabbits by multi-omics studies.

CONCLUSION
In conclusion, the current study identified key microbial taxa associated with finishing weight. We also determined the gut microbial CAZymes, GOs, KOs, and KEGG pathways associated with finishing weight. We found that the gut microbiome could act as a key predictor of finishing weight variation. In addition, we emphasized the important effects of SCFAs on finishing weight, especially butyrate. Hence, our findings provide essential insights into how the gut microbiome and SCFAS affect finishing weight and imply that manipulating the gut microbial community could be an efficient strategy to improve finishing weight in the meat rabbit industry.

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 in the article/ Supplementary Material.

ETHICS STATEMENT
The animal study was reviewed and approved by the Animal Care and Use Committee (ACUC) in Fujian Agriculture and Forestry University.

AUTHOR CONTRIBUTIONS
QG designed the experiments, analyzed the data, and wrote and revised the manuscript. SF and XC performed the experiments, analyzed the data, and wrote the manuscript. LZ, XY, and SX performed the experiments. All authors read and approved the final manuscript.

ACKNOWLEDGMENTS
We are grateful to the manager Mr. Zhoulin Chen of the commercial Ira rabbit farm who provides the opportunity to collect samples and phenotype measurement and the workers who are making contributions to the management of the experimental rabbit population.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb. 2020.01835/full#supplementary-material FIGURE S1 | Finishing weight phenotypic values of all rabbits (A) and high and low individuals (B). FIGURE S2 | The finishing weight associated species formed two clusters using Ward clustering algorithm.