Dynamic distribution of gut microbiota in meat rabbits at different growth stages and relationship with average daily gain (ADG).

BACKGROUND
The mammalian intestinal tract harbors diverse and dynamic microbial communities that play pivotal roles in host health, metabolism, immunity, and development. Average daily gain (ADG) is an important growth trait in meat rabbit industry. The effects of gut microbiota on ADG in meat rabbits are still unknown.


RESULTS
In this study, we investigated the dynamic distribution of gut microbiota in commercial Ira rabbits from weaning to finishing and uncover the relationship between the microbiota and average daily gain (ADG) via 16S rRNA gene sequencing. The results indicated that the richness and diversity of gut microbiota significantly increased with age. Gut microbial structure was less variable among finishing rabbits than among weaning rabbits. The relative abundances of the dominant phyla Firmicutes, Bacteroidetes, Verrucomicrobia and Cyanobacteria, and the 15 predominant genera significantly varied with age. Metagenomic prediction analysis showed that both KOs and KEGG pathways related to the metabolism of monosaccharides and vitamins were enriched in the weaning rabbits, while those related to the metabolism of amino acids and polysaccharides were more abundant in the finishing rabbits. We identified 34 OTUs, 125 KOs, and 25 KEGG pathways that were significantly associated with ADG. OTUs annotation suggested that butyrate producing bacteria belong to the family Ruminococcaceae and Bacteroidales_S24-7_group were positively associated with ADG. Conversely, Eubacterium_coprostanoligenes_group, Christensenellaceae_R-7_group, and opportunistic pathogens were negatively associated with ADG. Both KOs and KEGG pathways correlated with the metabolism of vitamins, basic amino acids, and short chain fatty acids (SCFAs) showed positive correlations with ADG, while those correlated with aromatic amino acids metabolism and immune response exhibited negative correlations with ADG. In addition, our results suggested that 10.42% of the variation in weaning weight could be explained by the gut microbiome.


CONCLUSIONS
Our findings give a glimpse into the dynamic shifts in gut microbiota of meat rabbits and provide a theoretical basis for gut microbiota modulation to improve ADG in the meat rabbit industry.

performance of food producing animals [4,5]. Recently, the role of gut microbiota in production performances has received more attention. For instance, Díaz-Sánchez et al. suggested using gut microbiota as biomarkers for prediction of production performance in the selective breeding process of chicken [6]. Warne et al. indicated that manipulation of gut microbiota during critical developmental windows could affect production performance in animals [7]. Hence, it is important to understand the dynamic distribution of gut microbiota in food producing animals at different growth stages and to recognize how microbial taxa and functions affect production performance.
Average daily gain (ADG) is an important index of production performance in commercial meat rabbit breeds as it is related to economic benefits for the meat rabbit industry. It has been reported that gut microbiota is intimately correlated with ADG in meat producing animals. He et al. revealed that the phylum Firmicutes showed a positive association with the ADG of meat duck, while the phylum Proteobacteria exhibited the inverse correlation [8]. Ma et al. indicated that broiler chickens with higher ADG harbored more abundant Christensenellaceae and Caulobacteraceae in the gut [9]. Jiao et al. found that Desulfopila was positively associated with ADG in beef cattle [10]. Ramayo-Caldas et al. suggested that increasing the proportion of Prevotella in the gut microbial community could improve porcine ADG [11]. However, the relationship between gut microbiota and ADG in meat rabbits remains unclear.
In this study, we investigated the dynamic distribution of gut microbiota in commercial Ira rabbits from weaning to finishing via high-throughput 16S rRNA gene sequencing. Additionally, we identified microbial taxa and potential functional capacities associated with ADG. Our results not only highlight the shifts and differences of gut microbial communities in meat rabbits at different growth stages, but also provide important information for improving ADG by modulating gut microbiota.

Sequencing data and microbial diversity analysis
The 16S rRNA gene sequencing process generated a total of 30,080,332 paired-end reads in both weaning and finishing samples. After sequences filtering steps and chimera removal, a total of 29,742,985 high-quality reads were obtained (15,380,901 clean reads of weaning samples, 14,362,084 clean reads of finishing samples). Based on 97% sequences identity, 1460 and 1586 OTUs were obtained from samples at weaning and finishing, respectively. A total of 2072 OTUs were identified from all samples, with 974 of those defined as core OTUs due to their existence in both weaning and finishing samples (Additional file 1: Fig. S1a). Additionally, 486 OTUs were uniquely identified in weaning samples and 612 unique OTUs were found in finishing samples.
The observed species, Shannon, and Good's coverage index of alpha diversity were used to analyze the richness and diversity of the gut microbial community. As shown in Fig. 1a, the observed species index of the finishing samples was significantly higher than weaning samples (FDR adjusted p < 0.0001). Finishing samples also presented significantly higher Shannon index values than weaning samples did (Fig. 1b, FDR adjusted p < 0.0001). In addition, Good's coverage index values were 0.999 and 0.998 for weaning and finishing rabbits, respectively, all showing excellent coverage (Additional file 1: Fig. S1b).
PCoA was performed to uncover the changes in gut microbial community structures in weaning and finishing samples. Both the unweighted and weighted UniFrac distances indicated that the samples were clustered by age (Fig. 1c, Additional file 1: Fig. S1c). Besides, the unweighted and weighted UniFrac distance metric comparison analysis further revealed less variation in the gut microbiota among finishing samples than among weaning samples (Fig. 1d, Additional file 1: Fig. S1d, FDR adjusted p < 0.0001).

Gut microbiota composition at different taxonomical levels
Gut microbial community composition was analyzed at the phylum and genus level. A total of 10 phyla were shared by weaning and finishing samples, including Actinobacteria, Bacteroidetes, Cyanobacteria, Euryarchaeota, Firmicutes, Lentisphaerae, Proteobacteria, Saccharibacteria, Tenericutes, and Verrucomicrobia (Additional file 1: Fig. S2a). In addition, Synergistetes and Planctomycetes were uniquely found in weaning and finishing samples, respectively. Firmicutes, Bacteroidetes, Verrucomicrobia, Proteobacteria, Cyanobacteria, and Tenericutes were the top six phyla in all samples and they comprised more than 99% of the total sequences (Fig. 2a). Among these, four phyla showed significant differences in relative abundances between weaning and finishing samples (Fig. 2a, Additional file 2: Table S2). The relative abundance of Firmicutes significantly increased from weaning (55.38%) to finishing (75.85%). Conversely, a significant decrease in the relative abundance of Bacteroidetes was observed from weaning (33.39%) to finishing (17.65%). Moreover, weaning samples exhibited a significantly higher abundance of Verrucomicrobia than finishing samples did (8.62% vs. 2.22%, respectively), while finishing samples exhibited a significantly higher abundance of Cyanobacteria than weaning samples did (1.41% vs. 0.82%, respectively).
Akkermansia belongs to the phylum Verrucomicrobia, Alistipes, Bacteroides, and Parabacteroides are members of the phylum Bacteroidetes, and the other 14 genera derive from phylum Firmicutes. Among these genera, the relative abundances of 15 were significantly changed from weaning to finishing (Additional file 2: Table S2). Similar to  the dynamic changes in the relative abundances of the dominant phyla, we found that the relative abundances of Akkermansia, Alistipes, Bacteroides, and Parabacteroides significantly declined from weaning to finishing, while the relative abundances of most genera (e.g., Ruminococcaceae_ UCG-010, Ruminococcaceae_NK4A214_group, Christense-nellaceae_R-7_group, and Ruminococcaceae _V9D2013_ group) of the phylum Firmicutes significantly increased from weaning to finishing.

Potential functional profile of gut microbial community
To investigate the changes in gut microbial functional profiles of weaning and finishing samples, KOs and KEGG pathways were predicted based on 16S rRNA gene sequencing data using Tax4Fun. We identified 6367 common KOs in all samples, while 41 KOs were specific to weaning samples and 137 KOs were unique to finishing samples (Additional file 1: Fig. S3a). In total, 177 KEGG pathways were identified and shared by both Fig. 2 The dynamic distributions of gut microbiota at different taxonomic levels. a At phylum level. b At genus level. The IDs on the X-axis with the same number but different letters ("W" and "F") in the two groups represent the same rabbit at weaning and finishing stage, respectively  S3b). To identify different enrichment of functional capacities between weaning and finishing samples, we performed a LEfSe analysis using the relative abundances of common KOs and pathways. As shown in Fig. 3a Table S3). Similarly, KEGG pathways related to the metabolism of glycan (e.g., glycosaminoglycan and other glycan degradation), monosaccharides (e.g., fructose, mannose, and galactose metabolism), vitamins (e.g., lipoic acid and biotin metaboilsm) were more abundant in weaning samples, while KEGG pathways related to the metabolism of amino acids (e.g., cysteine, methionine, arginine and proline metabolism) and short-chain fatty acids (SCFAs) (e.g., butanoate mand propanoate metabolism) were enriched in finishing samples.

Microbial taxa and potential functional capacities associated with ADG
To identify microbial taxa associated with ADG, we performed the two-part model analysis using the ADG phenotypic values adjusted for sex and cage effects, and the relative abundances of OTUs in finishing samples. A total of 34 OTUs were identified that exhibited significant associations with ADG. Among these OTUs, 18 were positively associated with ADG and 16 had negative associations with ADG. We annotated these OTUs to microbial taxa using the SILVA database ( Fig. 4, Additional file 2: Table S4), the phylogenetic relationships and abundances of these OTUs were shown in Additional file 1: Fig. S4 and Fig. S5. Among the positive ADG-associated OTUs, five were annotated to family level, including two OTUs to each of Bacteroidales_S24-7_group and Ruminococcaceae, and one OTU to Erysipelotrichaceae. At the genus level, four OTUs were annotated to Ruminococcaceae_UCG-014, two OTUs to Ruminococcaceae_UCG-010, and one OTU to each of Ruminiclostridium_5, Erysipelatoclostridium, Lactonifactor, Rikenella, Coprococcus_1, and Lachnospiraceae_NK4A136_group. One OTU was annotated to species Clostridium sp. ID11. Among the negative ADG-associated OTUs, three were annotated to each of the families Clostridiales_ vadinBB60_group and Lachnospiraceae. Four OTUs were annotated to the genus Eubacterium_coprostanoligenes_group, two to the genus Christensenellaceae_ R-7_group, and one to each of the genera Ruminococcaceae_UCG-009, Ruminococca-ceae_UCG-005, and Parasutterella. At the species level, one OTU was annotated to Sphingomonas paucimobilis.
To identify potential functional capacities correlated with ADG, Spearman correlation analysis was performed for the relative abundances of KEGG items and ADG phenotypic values. As shown in Fig. 5a (and Additional file 2: Table S5) Table S5). Among these, pantothenate and CoA biosynthesis, biotin metabolism, lysine degradation, arginine and proline metabolism, butanoate metabolism, propanoate metabolism, and glycine, serine and threonine metabolism were positively correlated with ADG, while NOD-like receptor signaling pathway, D-Alanine metabolism, antigen processing and presentation, cyanoamino acid metabolism, PPAR signaling pathway, phenylalanine metabolism, glycerolipid metabolism, and D-Glutamine and D-glutamate metabolism were negatively associated with ADG.

Phenotypic variation of ADG explained by gut microbiome
To estimate what proportion of variation in ADG could be explained by the microbiome, we performed 100 times cross-validation analyses at different p value thresholds (ranging from 10 − 4 to 0.1). As shown in Fig. 6, we found that the OTUs identified at p = 1.0 × 10 − 4 could explain 5.83% of the variations in ADG. At p = 0.1, the gut microbiome explained 10.42% of the variations in ADG, given that more OTUs were included in the analysis as the threshold increased.

Discussion
Recently, investigations into mammalian gut microbiota have revealed its vital role in host metabolism, physiology, and immunology. However, few reports have been published on the dynamic distributions of gut microbiota at different growth stages of commercial meat rabbits. ADG is an important growth trait, which is inevitably related to growth stage, but the relation between the gut microbiota and ADG in commercial meat rabbits is still elusive. In the present study, we analyzed the gut microbial diversities and abundances of Ira rabbits at weaning and finishing age. The potential functional profiles of gut microbial communities were then predicted and compared between weaning and finishing rabbits. In addition, we identified microbial taxa and functional capacities associated with ADG.
We found that the richness and diversity of gut microbial communities significantly increased with age ( Fig. 1a-b). This result is in agreement with previous findings for animals and humans. Niu et al. reported that gut microbial richness and diversity significantly increased with age in pigs, and suggested that growth stages and conditions were important factors affecting gut microbiota [12]. He et al. conclued that diet and physiological changes may have resulted in the observed increase of the richness and diversity of gut microbiota with increasing age in camels [13]. Additionally, maternal enteric microbes exposure level and breast-milk feeding rate significantly affected gut microbial diversity and richness in the early life was in human infants [14]. Furthermore, we investigated the variation in microbial community structure from weaning to finishing. PCoA analysis suggested that samples were clustered according to age. Unifrac distances further indicated that the gut microbiota of finishing samples showed significantly greater similarity than those of weaning samples (Fig. 1c-d, Additional file 1: Fig. S1c-d).
These results could be explained by the fact that all the weaned rabbits fed with the same fatten until finishing, which could minimize the effects of maternal environment and genetics on gut microbiota [15,16].
As with previous studies [17,18], we found that Firmicutes, Bacteroidetes, Verrucomicrobia, Proteobacteria, Cyanobacteria, and Tenericutes were the dominant phyla in both weaning and finishing samples. Moreover, the significant increases in the relative abundances of Firmicutes and Cyanobacteria were observed as the rabbits aged, along with a significant decline in the abundances of Bacteroidetes and Verrucomicrobia (Fig. 2a, Additional file 2: Table S2). Zhu et al. also noted that the percentage of Firmicutes increased as rabbits aged and demonstrated that Firmicutes played a vital role in dietary fiber and cellulose degradation [19]. During the fattening period, the rabbits were fed a diet with a high fiber content which should further stimulate the growth of Firmicutes (Additional file 2: Table S1). The intestinal environment became more anaerobic as the host gradually matured, and the increase in the abundances of Cyanobacteria could facilitate obligate anaerobic fermentation and synthesis of vitamins [20]. Bacteroidetes can break down polysaccharides and proteins in breast milk and diet, and facilitate the development of intestinal immune system [21,22]. Hence, the high prevalence of Bacteroidetes present in the gut microbial community of Fig. 5 Heatmap of predicted KEGG Orthologs (a) and pathways (b) significantly associated with ADG (FDR adjusted p < 0.05, |r| > 0.4). The correlation coefficient was used to plot weaning rabbits fits with feeding pattern (breast milk and solid feed mix-feeding) and immune physiology function development stage. Verrucomicrobia is a relatively newly-defined phylum with largely unknown functions. Interestingly, Akkermansia muciniphila, a wellknown mucin-degrading bacterium belongs to this phylum, which could improve nutrients extraction during cecotrophy in weaning rabbits. However, its overgrowth in finishing rabbits is closely related to the incidence of epizootic rabbit enteropathy [23,24]. This could be used to explain the decrease in relative abundances of Verrucomicrobia in healthy finishing rabbits.
As shown in Fig. 2b (and Additional file 2: Table S2), from weaning to finishing, we observed a significant decrease in the relative abundances of genera Alistipes, Bacteroides, and Parabacteroides, which belong to the phylum Bacteroidetes. According to the published reports, Alistipes is able to degrade dietary polysaccharides and flourishes in the intestine when host innate immunity at immaturity [25,26]. Bacteroides constitutes essential components of the mammalian intestinal microbiota that is capable of degrading breast milk polysaccharides and stimulating the formation of intestinal mucosa during infancy [27,28]. Parabacteroides is another gut bacteria that participates in breast milk polysaccharides metabolism and its abundance significantly declines with formula milk feeding [29,30]. Additionally, the abundance of Akkermansia (a genus of the phylum Verrucomicrobia) significantly declinesd as rabbits aged, which is in accordance with the results observed for gut microbiota development study in foals [31]. By contrast, genera from the phylum Firmicutes, such as Ruminococcaceae_UCG-010, Rumino-coccaceae_NK4A214_group, Christensenellaceae_R-7_group, Ruminococcaceae_V9D2013_group, and Ruminococcaceae_ UCG-014 exhibited higher abundances in finishing samples. These genera are suggested to exert key roles in dietary cellulose, hemicellulose, and lignocellulose fermentation and SCFAs production [32][33][34].
Comparison analysis of gut microbial functional capacities indicated that both KOs and KEGG pathways related to the metabolism of monosaccharides and vitamins were enriched in the weaning samples, while those related to the metabolism of amino acids and polysaccharides were more abundant in the finishing samples (Fig. 3, Additional file 2: Table S3). These alternations in the functional profiles should be correlated with the dynamic shifts in gut microbiota at different taxonomic levels. For instance, mannose is an important monosaccharide for protein glycosylation in mammals, and mannose metabolism associated with a higher percentage of Bacteroides has been reported in previous studies [35,36]. Galactose is an essential component of milk oligosaccharides, and a higher relative abundance of Bacteroides correlated with the utilization of galactose was observed in the gut microbiome of nursing piglets [37]. Both Bacteroides and Akkermansia are involved in the metabolism of vitamins, such as biotin and retinol [38,39]. Previous studies demonstrated that specific bacteria of phylum Firmicutes involved in amino acids (e.g., cysteine and arginine) metabolism could enhance intestinal mucosa immunity and reduce intestinal oxidative stress during the post-weaning period [40,41]. Besides, both Meale et al. and Ke et al. suggested that dietary polysaccharides metabolism is enhanced with an increase in abundance of Firmicutes [42,43].
Thirty-four OTUs were significantly associated with ADG (Fig. 4, Additional file 1: Fig. S4 and Fig. S5, Additional file 2: Table S4). Among these, OTUs positively associated with ADG were mostly annotated to members of family Ruminococcaceae (e.g., Ruminococcaceae_ UCG-014, Ruminiclostridium_5, and Ruminococcaceae_ UCG-010) and Bacteroidales_S24-7_group. These bacteria are able to produce butyrate by degrading indigestible fibers and polysaccharides [44][45][46]. Butyrate is not only an energy source for gut microbial growth, but has also been linked to intestinal epithelial cell proliferation and heat shock protein 70 (Hsp70) production [47]. Hsp70 plays an important role in maintaining the functional and structural properties of intestinal epithelial cells in response to pathogens challenge and oxidative stress during weaning transition [48,49]. These actions of butyrate that promote the intestinal development and health of productive animals are crucial for improving  [50]. Importantly, butyrate can maintain metabolic homeostasis and modulate immune and inflammatory responses via binding to G protein-coupled receptors FFAR3 and GPR109A, respectively [51,52], and the optimized metabolic and immune status is beneficial for the growth of farm animals [53]. In addition, butyrate produced by gut microbiota has been found to modulate the growth of animals by inducing secretion of intestinal satiety hormones [e.g., peptide tyrosine tyrosine (PYY) and glucagon-like peptide 1 (GLP1)] and growth hormones [e.g., growth factor insulin like growth factor 1 (IGF-1)] [54,55]. In contrast, OTUs annotated to Eubacterium_coprostanoligenes_group, Christense-nellaceae_R-7_group, Parasutterella, and Sphingomonas paucimobilis showed negative associations with ADG. Oral administration of Eubacterium_coprostanoligenes_ group bacteria in mice reduced body weight by affecting cholesterol metabolism [56,57]. A Study of the relationship between gut microbiota and ADG in meat ducks indicated that a decline in the abundance of Eubacterium_ coprostanoligenes_group was accompanied by an increase in ADG [58]. Christensenellaceae has long been known to affect host body weight [59]. This may be due to its negative correlation with the ratio of goblet cell to villus height, which affects host's inherent immunity and nutrient absorption [60,61]. Additionally, Parasutterella and Sphingomonas paucimobilis are opportunistic pathogens related to intestinal barrier dysfunction and inflammation, which are detrimental to host growth [62,63].
We also identified 125 KOs and 25 KEGG pathways showed potential correlations with ADG (Fig. 5, Additional file 2: Table S5). Our results suggested that both KOs and KEGG pathways related to the metabolism of vitamins, basic amino acids, and SCFAs were positively associated with ADG. Vitamins produced by gut microbiota could prevent oxidative stress and anti-inflammation, which maintains intestinal homeostasis and has been linked to promote food intake and body weight gain in mice [64,65]. In addition, the metabolism of basic amino acids and SCFAs mediated by gut microbiota was shown to affect body weight in pigs. Li et al. found that the metabolism of arginine, butanoate, and propanoate was more active in the gut microbiota of pigs with higher body weight gain [66]. Cheng et al. suggested that gut microbiota modulate increases in the concentrations of acetate, propionate, and butyrate, which contribute to a higher ADG in pigs [67]. In contrast, gut microbiota that participated in aromatic amino acids metabolism and immune response exhibited negative associations with ADG. Previous studies revealed that gut microbiota involved in aromatic amino acids metabolism related to colitis, and producing specific metabolites could reduce host body weight gain [68][69][70]. Immune response pathways mediated by gut microbiota included the NOD-like receptor signaling pathway and PPAR signaling pathway. These pathways are related to pro-inflammatory cytokines production and immune cells proliferation which exert negative effects on host health and growth [71,72].
Additionally, we found that the gut microbiome could explain 5.83-10.42% (Fig. 6) of the variation in ADG that effects similar to host genetics on ADG (5.2-9.6%) [73]. This result implies that the effects of gut microbiota should not be underestimated in attempts to improve growth performances of meat rabbits.

Conclusions
Our study characterized the gut microbiota profiles of weaning and finishing Ira rabbits. Gut microbial richness and diversity increased with age. Significant differences in gut microbial structure were observed between weaning and finishing rabbits. Dynamic shifts in microbial taxa at the phylum and genus level were uncovered between weaning and finishing rabbits. The metagenomic predicted KOs and KEGG pathways exhibited differential enrichment in weaning and finishing rabbits. Our results emphasized the importance of both butyrate producing bacteria and gut microbiota that modulate the metabolism of vitamins, basic amino acids, and SCFAs in promoting the ADG of meat rabbits. In addition, we found that gut microbiome had a similar effect size on ADG as host genetics. Taken together, our results improve our comprehensive understanding of the dynamic distributions of gut microbial communities in meat rabbits, and offer a direction for gut microbiota modulation to improve ADG in the meat rabbit industry.

Animals and sample collection
ADG is a complex quantitative trait [74], therefore, a total of 105 Ira rabbits (53 males and 52 females) were used in the present study derived from Laidewang Animal Husbandry Co., Ltd., Sanming, China. Six to eight pup rabbits per cage were raised with their dam under natural light and room temperature in the same commercial farm. A commercial pellet diet (details are shown in Additional file 2: Table S1) was provided to lactating dams twice a day and pup rabbits had free access to the feed. Pup rabbits were weaned at 28 ± 2 days, at which point one or two rabbits were randomly selected from 90 cages to measure weaning body weight. Hard fecal samples were collected by stimulating the anus. To reduce the effect of differences in initial body weight, 105 rabbits had similar weaning weight (0.9 ± 0.06 kg) were selected and randomly assigned to separate cages (one rabbit per cage) and fed with a fattening diet (details are shown in Table S1) until finishing (72 ± 2 days). Finishing body weight was measured to calculate ADG and hard fecal samples were collected. All rabbits were healthy and had not received antibiotics, anticoccidial drugs, probiotics or prebiotics during experimental period. All fecal samples were dipped into liquid nitrogen for transportation and stored at − 80°C in the laboratory. At the end of experiments, all rabbits were transported to the local slaughterhouse, stunned with electronarcosis (80 V for 10 s) and quickly bled by cutting the jugular veins and carotid arteries. To avoid the effect of artificial bias, the authors were not involved in the processes of rabbits' selection, grouping, body weight measurement, and further 16S rRNA gene sequencing.

DNA extraction and 16S rRNA gene sequencing
Microbial genomic DNA was extracted from feces using the QIAamp Fast DNA Stool Mini Kit (QIAGEN, Hilden, Germany) following the manufacturer's instructions. Before PCR amplification and sequencing, we assessed the purity and integrity of total DNA by using the Nanodrop ND-2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) and 1.5% agarose gel electrophoresis, respectively.
The V3-V4 hypervariable region of the 16S rRNA gene was amplified by the barcoded fusion primers 341F (5′-CCTACGGGNGGCWGCAG-3′) and 806R (5′-GGACTA CHVGGGTATCTAAT-3′). The PCR conditions were as follows: initial denaturation step at 95°C for 3 min, followed by 28 cycles of 95°C for 30 s, 55°C for 30 s and 72°C for 30 s, and a final extension step at 72°C for 10 min. The PCR products purification was performed using the Agencourt AMPure XP system (Beckman Coulter, Brea, CA, USA). The final DNA libraries were sequenced on a Hiseq-2500 platform (Illumina, San Diego, CA, USA) according to the manufacturer's instructions.

16S rRNA gene sequencing data analysis
Quality control of raw data was performed using QIIME (v.1.9.1), including the removal of the primers, barcodes, and low quality sequences [75]. FLASH (v.1.2.11) was used to merge high-quality paired-end reads into tags [76]. To normalize the sequencing depth, we rarefied the library size of microbial sequences to 40,000 tags per sample before further analysis [77]. Tags were clustered into operational taxonomic units (OTUs) at 97% sequence identity using USEARCH (v.10.0) [78]. We filtered out those OTUs with relative abundance < 0.1% and those that were presented in less than 3% of the experimental rabbits from further analysis. The SILVA database (v.132) was used to assign taxonomic category to OTUs [79]. The alpha and beta diversity indices were calculated using Mothur (v.1.41.1) and QIIME (v.1.9.1), respectively [75,80]. Potential functional capacities of gut microbiota were predicted using Tax4Fun [81].

Statistical analysis
Wilcoxon test with false discovery rate (FDR) correction was used to determine differences in observed species, Shannon, Good's coverage, weighted and unweighted UniFrac distance metric between weaning and finishing rabbits. Principal coordinate analysis (PCoA) was performed using both unweighted and weighted UniFrac distances. The dynamic changes in the relative abundances of microbiota at the phylum and genus level were presented as alluvial diagrams. The numbers of shared OTUs, phyla, genera, KOs, and KEGG pathways were showed as the Venn diagrams. Linear discriminant analysis Effect Size (LEfSe) was used to analyze the differential enrichment of KOs and KEGG pathways in weaning and finishing samples.
After sex and cage effects correction, the residuals of ADG were used for further association analysis between ADG phenotypic values and the relative abundances of OTUs. To overcome the problem of non-normal distribution of OTUs, a two-part model analysis was performed to identify microbial taxa associated with ADG as described previously [82]. Briefly, the two-part model association analysis consisted of binary, quantitative, and meta models. The binary model describes a binomial analysis that tests for associations between ADG and detection of a microbe. The quantitative model tests for associations between ADG and the abundance of microbes, but only the samples where the microbe was present were included in the analysis. The meta model was used to combine the effects 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 metaanalysis. Skewness correction was performed by 1000 × permutation tests. FDR < 0.05 was set as the significance threshold. The phylogenetic relationships of the identified microbes were analyzed using neighbor-joining algorithm and heatmap of their abundances were generated using the SEED2 program [83]. Spearman's correlation analysis with FDR correction was performed to uncover ADGassociated KOs and KEGG pathways.
To investigate the contribution of the gut microbiome to the variation in ADG, a 100 × cross-validation was performed as described by Fu et al. [84]. We randomly divided the dataset into an 80% discovery dataset and a 20% validation dataset. In the discovery dataset, the two-part model association analysis was performed to identify a number of (n) OTUs that were significantly associated with phenotype at a certain p value and assessed the effect sizes of binary and quantitative features (β 1 and β 2 ) of each OTUs. In the validation dataset, the effect of gut microbiome on ADG (r m ) for each individual was estimated by an additive model: r m = P n j¼1 ðβ 1 þ b j þ β 2 q j Þ , where b j and q j represents the binary and quantitative feature of j OTU, respectively. We calculated the squared correlation coefficient (R 2 ) between r m and the phenotypic value (corrected for sex and cage), which represents the phenotypic variance explained by the gut microbiome. We repeated the crossvalidation for 100 times and calculated the average value of the explained variations to ensure validity and stability of the estimation.
Additional file 1: Figure S1. Comparison of OTUs and other diversity index between weaning and finishing samples. (a) Venn diagram to describe the common and unique OTUs between the two groups. (b) Good coverage index (c) PCoA analysis based on Weighted Unifrac distance. (d) Weighted Unifrac distance metric. Figure S2. The Venn diagram representation of the share phyla (a) and genera (b) between weaning and finishing samples. Figure S3. The Venn diagram representation of the share KOs (a) and KEGG pathways (b) between weaning and finishing samples. Figure S4. The phylogenetic relationships of ADG-associated microbial taxa. The coral and blue tree labels represent for positive and negative ADG associated OTUs. Bootstrap values are shown on the branches. Figure S5. The heatmap of abundances of ADG-associated microbial taxa. The coral and blue strips correspond to positive and negative ADG associated OTUs.
Additional file 2: Table S1. Composition of pellet diet at different growth stages. Table S2. Differences in the relative abundances of dominant phyla and predominant genera. Table S3. Different enriched KOs and KEGG pathways between weaning and finishing samples. Table  S4.
OTUs showing significant associations with ADG and taxonomic assignment using SILVA database. Table S5. KOs and KEGG pathways significantly associated with ADG.