Velasco-Galilea et al. Animal Microbiome https://doi.org/10.1186/s42523-020-00059-z (2020) 2:40 Animal Microbiome RESEARCH ARTICLE Open Access Breeding farm, level of feeding and presence of antibiotics in the feed influence rabbit cecal microbiota María Velasco-Galilea1* , Miriam Guivernau2, Miriam Piles1, Marc Viñas2, Oriol Rafel1, Armand Sánchez3,4, Yuliaxis Ramayo-Caldas1, Olga González-Rodríguez1 and Juan P. Sánchez1 Abstract Background: The effect of the production environment and different management practices in rabbit cecal microbiota remains poorly understood. While previous studies have proved the impact of the age or the feed composition, research in the breeding farm and other animal management aspects, such as the presence of antibiotics in the feed or the level of feeding, is still needed. Characterization of microbial diversity and composition of growing rabbits raised under different conditions could help better understand the role these practices play in cecal microbial communities and how it may result in different animal performance. Results: Four hundred twenty-five meat rabbits raised in two different facilities, fed under two feeding regimes (ad libitum or restricted) with feed supplemented or free of antibiotics, were selected for this study. A 16S rRNA gene- based assessment through the MiSeq Illumina sequencing platform was performed on cecal samples collected from these individuals at slaughter. Different univariate and multivariate approaches were conducted to unravel the influence of the different factors on microbial alpha diversity and composition at phylum, genus and OTU taxonomic levels. The animals raised in the facility harboring the most stable environmental conditions had greater, and less variable, microbial richness and diversity. Bootstrap univariate analyses of variance and sparse partial least squares-discriminant analyses endorsed that farm conditions exerted an important influence on rabbit microbiota since the relative abundances of many taxa were found differentially represented between both facilities at all taxonomic levels characterized. Furthermore, only five OTUs were needed to achieve a perfect classification of samples according to the facility where animals were raised. The level of feeding and the presence of antibiotics did not modify the global alpha diversity but had an impact on some bacteria relative abundances, albeit in a small number of taxa compared with farm, which is consistent with the lower sample classification power according to these factors achieved using microbial information. Conclusions: This study reveals that factors associated with the farm effect and other management factors, such as the presence of antibiotics in the diet or the feeding level, modify cecal microbial communities. It highlights the importance of offering a controlled breeding environment that reduces differences in microbial cecal composition that could be responsible for different animal performance. Keywords: Cecal microbiota, Meat rabbit, Breeding farm, Feed restriction, Antibiotics, 16S MiSeq Illumina sequencing, Analysis of variance, Multivariate approach * Correspondence: maria.velasco@irta.es 1 Institute of Agrifood Research and Technology (IRTA) – Animal Breeding and Genetics, E08140 Caldes de Montbui, Barcelona, Spain Full list of author information is available at the end of the article © The Author(s). 2020 Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. Velasco-Galilea et al. Animal Microbiome (2020) 2:40 Page 2 of 16 Background antibiotics. Quantitative feed restriction is a widely Microbial communities that inhabit the gastrointes- applied commercial practice which consists of redu- tinal tract (GIT) of animals constitute a complex eco- cing the amount of feed the animal would consume system whose members constantly interact between by a certain percentage when the food is provided themselves and with their host [1]. These interactions ad libitum. Gidenne et al. (2009) [11] demonstrated ensure homeostatic balance maintenance since GIT that feed restriction, despite penalizing animal ecosystem components are involved in many physio- growth, improves feed efficiency and reduces mortal- logical and immunological processes [2]. In the case ity due to enteric disorders. It is hypothesized that of the domestic meat rabbit (Oryctolagus cuniculus), a these positive effects could be partially explained by small herbivorous mammalian belonging to the family changes in gut microbial composition or activity Leporidae, cecum is the main organ for microbial fer- originated by the application of feed restriction. mentation. Thus, it is not surprising that the rabbit However, techniques used so far to study this pos- cecum hosts the richest and the most diverse micro- sible association have found no evidence of it [11]. bial community of its GIT [3]. For this reason, the This study, which comprises a large number of ani- cecum has been the organ preferably chosen in previ- mals in an experimental design involving different man- ous rabbit gut microbiota assessments [4–7]. agement and environmental factors, is intended to Thanks to the development of next generation se- unravel changes in diversity and composition of rabbit quencing (NGS) technologies, and their rapidly de- cecal microbial communities associated with these fac- creasing costs, it is currently possible to characterize tors. It will allow for a better understanding of how the the gut microbiota of a large number of animals. farm where the animal was raised, the presence of anti- This characterization allows a deeper comprehension biotics in the feed, and feed restriction shape the cecal of the differences between animals concerning their microbiota of growing rabbits. microbial composition and diversity. It is hypothe- sized that the production environment could par- Results tially mediate these differences. Our general aim is Sequence processing to provide further evidence of the effect of different After the removal of doubletons and samples with management and environmental factors on cecal mi- low sequence counts, 425 rabbit cecal samples (Add- crobial composition and diversity. In relation to this itional file 1) were represented on 14,928,203 se- topic, there is a certain amount of information quence counts clustered into 963 different OTUs. already published. A growing number of studies have Each sample had on average 35,125 final sequences revealed changes in rabbit cecal microbial communi- (range: 10,157-678,798) and 677 OTUs (range: 197– ties exerted by age [8] or the type of feed provided 841) (Additional files 2 and 3). Figure 1 shows two to the kits after weaning [6, 7]. Another factor that histograms representing the sample richness and the causes variation is the administration of antibiotics proportion of OTUs present across samples. Most of in the feed. Different molecules have been widely ad- the samples had more than 700 different OTUs ministered in rabbit meat production, especially after (mode = 748) and nearly 140 OTUs were present in weaning, to curb mortality peaks (sometimes over all the samples. 20%) as a result of the onset of gastrointestinal Taxonomic assignment of representative OTUs against symptoms [9]. Multiple studies have shown alter- the Greengenes reference database gg_13_5_otus (Add- ations caused in gut microbiota by the administra- itional file 4) revealed the presence of 8 different known tion of antibiotics in the feed [5, 10]. Despite the phyla with an average of 8 phyla per sample (range: 7–8) European Union having banned the use of antibiotics (Additional file 5), and 28 different known genera with in animal feeds as growth promoters since 2006 (EC an average of 24 genera per sample (range: 17–28) 1831/2003), at the time this experiment was con- (Additional file 6). ducted, the administration of a mix of up to four an- tibiotics was permitted to prevent or treat the emergence of potential infectious diseases on farms. Animal management and farm environment shaping Nowadays, the administration of only one antibiotic cecal microbial alpha diversity molecule is allowed and substantial efforts are being The study of alpha diversity was performed after rar- made towards searching for efficient alternatives efying the prefiltered and unnormalized OTU table which allow for a complete withdrawal of antibiotics to 10,000 sequences per sample. Rarefaction gener- in animal feeds. In this context, the application of ated a table which contained the sequence counts of feed restriction during the growing period was pro- 963 different OTUs for 425 samples. The average posed as an interesting alternative to the use of (standard deviation) number of observed OTUs Velasco-Galilea et al. Animal Microbiome (2020) 2:40 Page 3 of 16 Fig. 1 Sample richness and presence of CSS-normalized OTUs across samples within animal was 560.52 (75.03) and the average Furthermore, larger variability in both indexes was Shannon index within animal was 5.09 (0.26). The observed in farm A than in farm B. No significant comparison of alpha diversities revealed that the differences for the two alpha diversity indexes were group of animals raised in farm B had greater alpha found between feeding regimes within both farms diversity than the group of animals raised in farm A (Fig. 2, PFDR > 0.05), nor between the presence and (estimated differences of 40.20 (9.83) observed OTUs the absence of antibiotics in the feed within farm B and 0.17 (0.03) Shannon indexes; PFDR < 0.001). (Fig. 2, PFDR > 0.05). Fig. 2 Microbial richness and diversity between samples grouped according to management that animals received. The cecal microbial richness and diversity were estimated by the observed number of different OTUs and the Shannon indexes, respectively Velasco-Galilea et al. Animal Microbiome (2020) 2:40 Page 4 of 16 Animal management and farm environment shaping Odoribacter, Methanobrevibacter, Blautia, Butyrivibrio, cecal microbial composition Coprococcus, Dehalobacterium, Dorea, Oscillospira, rc4–4 According to the taxonomic assignment of representa- and Oxalabacter were overrepresented in cecal samples of tive sequences (Additional file 4) performed with the rabbits raised in farm B. Interestingly, genera Epulopis- UCLUST consensus taxonomy assigner on the Green- cium, p-75-a5, Phascolarctobacterium, Campylobacter genes reference database gg_13_5_otus, Firmicutes and Desulfovibrio were only found in samples collected (76.74%), Tenericutes (7.22%) and Bacteroidetes (6.26%) from farm A (Table 2). were the predominant phyla, accounting for more than The analyses on the CSS-normalized OTUs revealed 90% of the microbial diversity, in the rabbit cecal sam- that 648 out of the 946 OTUs showed signatures signifi- ples studied (Fig. 3). cantly different between farms. Out of these, 276 were overrepresented in farm A, while 372 were overrepre- Differential growth and cecal microbial composition across sented in farm B. Table S1 shows the estimated differ- farms ence between farms for these OTUs, their sequences The facility where the animals were raised affected their and their assignment at the lowest taxonomic level. Only growth performance. Animals raised in farm B exhibited 9 of them could be assigned at species level and 129 a faster growth (47.11 g/day) than those raised in farm A were assigned to known genera. These results showed (44.19 g/day). The estimated average daily gain differ- remarkable coincidences with those obtained from the ence between farm B and farm A was 2.92 ± 0.94 g per analyses directly performed on the relative abundance of day (P < 0.005). Cecal samples of rabbits raised in farm taxa at phylum and genera levels. An example that illus- A showed an overrepresentation of phyla Bacteroidetes, trates this match is the overrepresentation of genus Proteobacteria and Verrucomicrobia while phyla Eur- Akkermansia in farm A. This genus is encompassed by yarchaeota, Cyanobacteria and Firmicutes were found to phylum Verrucomicrobia that was also overrepresented be overrepresented in cecal samples of rabbits raised in in rabbits raised in farm A, as well as 6 out of the 7 farm B (Table 1). OTUs assigned to this phylum. Genera Ruminococcus (4.32%), Blautia (2.96%) and Oscillospira (2.37%) dominated the meat rabbit cecal Differential growth and cecal microbial composition across microbiota. Most of the relative abundance differences at feeding regimes genus level were found differentially represented between The feeding regime affected the rabbits’ growth perform- animals raised in the different farms: genera Bacteroides, ance in both facilities. Animals fed AL had a higher Parabacteroides, Rikenella, Anaerofustis, Anaerostipes, growth (48.74 and 55.77 g/day in farms A and B, re- Clostridium, Coprobacillus, Anaeroplasma and Akkerman- spectively) than those fed R (38.95 and 38.65 g/day in sia were overrepresented in cecal samples of rabbits raised farms A and B, respectively). The estimated average daily in farm A while genera Adlercreutzia, Butyricimonas, gain difference between AL and R groups was 9.79 ± Fig. 3 Phyla relative abundances of samples grouped according to farm, level of feeding and presence of antibiotics in the feed Velasco-Galilea et al. Animal Microbiome (2020) 2:40 Page 5 of 16 Table 1 Microbial composition at phylum level in cecal samples of rabbits grouped by farm Phylum Mean relative abundance in farm A Mean relative abundance in farm B Estimated difference farm A - farm PFDR (%) (SD) (%) (SD) B ± SE Actinobacteria 1.62 (0.67) 1.84 (0.33) − 0.14 ± 0.08 0.09 Bacteroidetes 6.84 (1.81) 4.03 (0.70) 2.74 ± 0.22 0.00 Cyanobacteria 0.77 (0.40) 1.05 (0.36) −0.39 ± 0.05 0.00 Euryarchaeota 0.13 (0.19) 0.44 (0.17) −0.28 ± 0.02 0.00 Firmicutes 75.83 (3.34) 79.66 (1.53) −3.78 ± 0.41 0.00 Proteobacteria 1.83 (0.62) 0.66 (0.12) 1.14 ± 0.07 0.00 Tenericutes 7.21 (1.47) 7.25 (0.93) 0.00 ± 0.18 0.99 Verrucomicrobia 1.62 (0.45) 0.91 (0.24) 0.68 ± 0.05 0.00 0.58 and 17.12 ± 1.08 g per day in farms A and B, re- per day; P = 0.78). Cecal samples of rabbits that received spectively (P < 0.001). An overrepresentation of phyla feed free of antibiotics showed an overrepresentation of Cyanobacteria (estimated difference R - AL = 0.11 ± 0.04; phyla Cyanobacteria compared to those that received PFDR = 0.04) and Verrucomicrobia (estimated difference feed supplemented with antibiotics (estimated difference R - AL = 0.11 ± 0.05; PFDR = 0.04) was found in cecal without antibiotics - with antibiotics = 0.49 ± 0.09; PFDR < samples of rabbits fed R and raised in farm A. On the 0.001). In addition, the analyses on the CSS-normalized other hand, phylum Euryarchaeota was overrepresented OTUs revealed an overrepresentation of 15 and 29 in animals fed R and raised in farm B (estimated differ- OTUs in cecal samples of rabbits that received a feed ence R - AL = 0.14 ± 0.04; PFDR < 0.001). At genus level, supplemented or free of antibiotics; respectively. Table the only significant contrast was observed for rc4–4 S3 shows the estimated difference between the presence which resulted overrepresented in samples from animals and the absence of antibiotics in the feed for the OTUs fed AL in farm A (estimated difference R - AL = − in which the differences reached the significance thresh- 0.03 ± 0.01; PFDR < 0.001) while in farm B none of the old. The OTU sequences as well as their assignment at genera resulted differentially represented (PFDR > 0.05) the lowest taxonomic level are also shown in Table S3. between feeding regimes. The contrasts based on the Only 1 of these OTUs could be assigned at species level CSS-normalized OTUs revealed 51 and 9 OTUs differ- (Bacteroides fragilis) and 2 OTUs at genus level (Oscil- entially represented between feeding regimes within lospira and Coprococcus). farms A and B, respectively. Within farm A, 32 OTUs were overrepresented in cecal samples of rabbits that Microbial information as a classifier of cecal samples were fed AL and 19 OTUs in the samples from rabbits according to farm environment and animal management fed R. Within farm B, 7 OTUs were overrepresented in Sparse partial least squares-discriminant analyses (sPLS- cecal samples of rabbits that were fed AL and 2 OTUs DA) on the CSS-normalized OTUs were conducted to were overrepresented in rabbits that were fed R. Table discriminate samples according to the factors considered S2 shows the estimated difference between feeding re- in this study (i.e., the farm where the animal was raised, gime within farm of these OTUs, their sequences and the presence or the absence of antibiotics in the feed their assignment at the lowest taxonomic level. The ana- and the feeding regime). The tuning process of the lyses based on the CSS-normalized OTUs within farm A sPLS-DA conducted to discriminate samples according were in full accordance with the analyses performed at to the farm where the rabbits were raised selected 5 genus level given that all OTUs assigned to genus rc4–4 OTUs for component 1 and 1 OTU for component 2 (phylum Firmicutes) were overrepresented in cecal sam- (Fig. 4). Component 1 explained 7.00% of the total vari- ples of rabbits fed AL. ance while component 2 explained 0.67%. The classifica- tion performance of this sPLS-DA could be said to be Effect of the presence of antibiotics in the feed perfect since its overall and balanced error rate (BER) The effect of the presence of antibiotics in the feed per class across 1000 replicates of 5-folds cross- could only be assessed within farm B given that all rab- validation runs was 0.00 (0.00). Furthermore, two OTUs bits raised in farm A received feed supplemented with of component 1 had a stability higher than 0.9. antibiotics. Animals that received antibiotics had a The sPLS-DA performed to discriminate samples slightly higher growth (47.29 g/day) than those that did across feeding regimes within farm A selected 70 OTUs not (46.59 g/day). The estimated average daily gain dif- for component 1 and 65 OTUs for component 2 (Fig. 5). ference between groups was not significant (0.69 ± 2.43 g Component 1 explained 2.34% of the total variance while Velasco-Galilea et al. Animal Microbiome (2020) 2:40 Page 6 of 16 Table 2 Relative abundances of genera, grouped by phylum, differentially represented between farms (PFDR < 0.05) Genus Mean relative abundance Mean relative abundance Estimated difference in farm A (%) (SD) in farm B (%) (SD) farm A - farm B ± SE Actinobacteria Adlercreutzia 0.89 (0.47) 1.14 (0.23) −0.19 ± 0.06 Bacteroidetes Bacteroides 1.88 (0.67) 0.80 (0.35) 1.10 ± 0.08 Butyricimonas 0.16 (0.19) 0.35 (0.17) −0.19 ± 0.02 Odoribacter 0.23 (0.21) 0.44 (0.20) −0.21 ± 0.03 Parabacteroides 0.25 (0.18) 0.07 (0.07) 0.18 ± 0.02 Rikenella 0.39 (0.24) 0.18 (0.13) 0.25 ± 0.03 Euryarchaeota Methanobrevibacter 0.13 (0.19) 0.44 (0.17) −0.28 ± 0.02 Firmicutes Anaerofustis 0.12 (0.08) 0.08 (0.04) 0.03 ± 0.01 Anaerostipes 0.17 (0.08) 0.12 (0.04) 0.06 ± 0.01 Blautia 2.86 (0.67) 3.22 (0.46) −0.36 ± 0.08 Butyrivibrio 0.10 (0.07) 0.13 (0.06) −0.03 ± 0.01 Clostridium 1.09 (0.26) 0.87 (0.13) 0.21 ± 0.03 Coprobacillus 0.20 (0.27) 0.14 (0.08) 0.08 ± 0.03 Coprococcus 1.96 (0.42) 2.26 (0.29) −0.28 ± 0.05 Dehalobacterium 0.05 (0.08) 0.18 (0.03) −0.13 ± 0.01 Dorea 0.46 (0.12) 0.51 (0.09) −0.05 ± 0.02 Epulopiscium 0.14 (0.11) 0.00 (0.00) 0.15 ± 0.01 Oscillospira 2.11 (0.53) 2.85 (0.31) −0.79 ± 0.07 p-75-a5 0.13 (0.06) 0.00 (0.00) 0.13 ± 0.01 Phascolarctobacterium 0.27 (0.24) 0.00 (0.00) 0.26 ± 0.03 rc4–4 0.13 (0.06) 0.23 (0.03) −0.10 ± 0.01 Proteobacteria Campylobacter 0.08 (0.08) 0.00 (0.00) 0.08 ± 0.01 Desulfovibrio 0.58 (0.22) 0.00 (0.00) 0.57 ± 0.03 Oxalabacter 0.10 (0.06) 0.13 (0.03) −0.03 ± 0.01 Tenericutes Anaeroplasma 0.23 (0.18) 0.10 (0.09) 0.12 ± 0.02 Verrucomicrobia Akkermansia 1.62 (0.45) 0.91 (0.23) 0.68 ± 0.05 component 2 explained 5.58%. The cross-validation as- Component 1 explained 3.05% of total variance and de- sessment of the classification performance of this sPLS- fined the discrimination between samples from animals DA showed an overall and BER per class of 0.27 (0.02). fed with antibiotics and those fed without antibiotics. The stability of 18 and 5 OTUs selected in components On the other hand, component 2 explained 3.05% of 1 and 2, respectively, across the different cross- total variance and defined the discrimination between validation folds was higher than 0.9. samples from animals fed R and those belonging to ani- Finally, the sPLS-DA conducted to discriminate sam- mals fed AL. The cross-validation assessment of the ples of animals raised within farm B according to the classification performance of this sPLS-DA showed an combination of the presence or not of antibiotics in the overall BER of 0.32 (0.15). The BER per class was 0.34 feed and the feeding regime selected 9 OTUs for compo- (0.12) for samples fed R without antibiotics, 0.46 (0.14) nent 1 and 70 OTUs for component 2 (Fig. 6). for samples fed AL without antibiotics, 0.29 (0.11) for Velasco-Galilea et al. Animal Microbiome (2020) 2:40 Page 7 of 16 Fig. 4 Sparse partial least squares discriminant analysis representing cecal samples of rabbits raised in farm A (blue) and in farm B (orange) samples fed R with antibiotics, and 0.20 (0.07) for sam- with previous studies that have characterized the rabbit ples fed AL with antibiotics. The stability of 3 and 11 cecal microbiota [5, 7, 12] and reported Firmicutes as the OTUs selected in components 1 and 2, respectively, predominant phylum. However, there are discrepancies across the different cross-validation folds was higher between studies in establishing which other phyla are also than 0.9. prevalent in this ecosystem. Whereas we found phyla Tenericutes and Bacteroidetes representing 7.22 and 5.93% Discussion of the cecal microbial composition, respectively, Chen The influences of farm environment and common com- et al. 2019 [7] and Zou et al. 2016 [5] reported Bacteroi- mercial practices of animal management on their gut detes as the second predominant phylum representing 18 microbiota are not yet well known in many livestock spe- and 20% of New Zealand White and Rex rabbit cecal mi- cies. In this study, we have aimed to disentangle potential crobial composition, respectively. Conversely, other stud- changes in microbial diversity and composition of meat ies that have previously characterized meat rabbit fecal rabbit cecal communities as a result of being raised in dif- microbiota identified higher relative abundances of phyla ferent farms and subjected to different handling during Proteobacteria and Verrucomicrobia [10, 13]. Velasco- their growing period. To shed light on this matter, we Galilea et al. 2018 [12] reported Firmicutes (76.42%), conducted a microbiota comparison of a large number of Tenericutes (7.83%) and Bacteroidetes (7.42%) as the pre- rabbits raised in different farms, feeding regimes, and fed dominant phyla of meat rabbit fecal and cecal microbial with feed supplemented or free of antibiotics. communities. These discrepancies found across studies could be attributed to technical issues (e.g., pair of 16S rRNA gene-based characterization of meat rabbit primers, sequencing platform, bioinformatic pipeline cecal microbiota employed to process raw sequences or reference database The Illumina MiSeq sequence processing of samples col- used for the taxonomic assignment of the representative lected from these animals revealed that phyla Firmicutes, sequences) or to purely biological reasons (e.g., breed, age Tenericutes and Bacteroidetes dominate the growing meat or section of the GIT sampled). Nonetheless, Kylie et al. rabbit cecal ecosystem representing more than 90% of its (2018) [13] depicted that the relative increase in less bene- entire microbial composition. This fact is in accordance ficial phyla, such as Proteobacteria, could be related to Velasco-Galilea et al. Animal Microbiome (2020) 2:40 Page 8 of 16 Fig. 5 Sparse partial least squares discriminant analysis representing cecal samples of rabbits raised in farm A and fed R (blue) or AL (orange) seasonal climate changes that directly impact rabbits’ indexes observed in this farm [13]. Despite not having health. This impact affects the susceptibility to enteritis observed significant differences between the presence and possibly feed conversion efficiency. In any case, this or not of antibiotic in the feed, nor between feeding phylum was more prevalent in farm A where the animals regimes, it is noteworthy to mention that samples col- were more exposed to changes in climate conditions. lected from animals fed AL in both farms had a greater, although not significant, richness than those Farm environment modify alpha diversity fed R. This fact is consistent with previous studies in Regarding the alpha diversity assessment, Shannon mice that observed a lower alpha diversity in animals and the observed number of OTUs indexes revealed with a restricted level of feeding [15–17]. Surprisingly, the existence of significant differences between the but in agreement with our results, studies performed experimental farm where the rabbits were raised. in pigs [18], chicken [19] and Rex rabbits [5] also did Cecal samples collected from rabbits raised in farm B not show clear significant differences on alpha diver- had greater richness and diversity than those belong- sity indexes between animals fed on diets with antibi- ing to animals raised in farm A. This could be ex- otics with respect to those on diets free of antibiotics. plained by more stable environmental conditions in Nevertheless, these studies were able to detect differ- farm B (i.e., facility better insulated) than in farm A. ences in the relative abundances of some specific spe- It has been already shown that intestinal health is cies between diets. For example, Kumar et al. 2018 positively associated with microbial diversity [14]. In [19] found that the inclusion of bacitracin in the feed our case, this better health could be said to be did not affect the chicken bacterial phyla. However, granted by the more stable environmental conditions they observed differences between the control and the offered by farm B. The most exposed environmental bacitracin-fed group in the ileal and cecal bacterial conditions of farm A, combined with the fact that populations at lower taxonomic levels. It is worth samples of animals raised in this facility were col- noting that the antibiotic withdrawal at the beginning lected from rabbits produced in 4 different batches, of the last week of the rabbits’ lives equalized the di- could also explain the larger variability in both ets of both groups and possibly their microbial Velasco-Galilea et al. Animal Microbiome (2020) 2:40 Page 9 of 16 Fig. 6 Sparse partial least squares discriminant analysis representing cecal samples of rabbits raised in farm B and fed R without antibiotics (blue), fed AL without antibiotics (orange), fed R with antibiotics (gray) and fed AL with antibiotics (green) populations, which may explain some lack of differ- them in meat rabbit cecum. Species of this genus are an- ences between them. aerobic Gram-negative members of the family Bacteroi- daceae that play an important role in the degradation of Farm environment has a large impact on rabbit cecal vegetal polysaccharides and amino acid fermentation in microbiota the mammal GIT [20, 21]. Moreover, this genus is in- Despite the lack of differences in microbial diversity and volved in propionic acid and lactate formation depend- richness across management factors; univariate studies ing on nitrogen organic availability. Nonetheless, some revealed differential microbial composition across the authors showed that great amounts of Bacteroides could studied factors. In addition, the performed multivariate predict obesity tendency. Parabacteroides is also an an- analysis evidenced a certain classification power of the aerobic Gram-negative bacterium (family Porphyromo- samples on the different levels of management and en- nadaceae) involved in amino acid transport and vironment factors based on the microbial composition metabolism, energy production and conversion, lipid of the samples. transport and metabolism, recombination and repair, cell As it might be expected, analyses of variance con- cycle control, cell division, and cell motility in the intes- firmed that the breeding farm strongly impacts meat tinal microbiota of the growing rabbit [22]. This genus rabbit cecal microbial composition. Our results revealed was specifically found in the cecal microbiota of mice that the relative abundances of 6 out of 8 phyla are dif- raised in conventional conditions and absent in those ferentially represented between both farms. At genus raised in pathogen-free facilities in a study performed level, we detected significant differences in the relative under different housing conditions [23]. abundances of almost all of them. Genera Bacteroides, Within the phylum Firmicutes, genus Clostridium Parabacteroides, Rikenella, Anaerofustis, Anaerostipes, (family Clostridiaceae) is an anaerobic Gram-positive Clostridium, Coprobacillus, Anaeroplasma and Akker- bacterium that inhabits the GIT of many mammals mansia were enriched in cecal samples of rabbits housed where it acts by degrading cellulose. However, some in farm A. The first three belong to phylum Bacteroi- Clostridium species (e.g., C. perfringens and C. difficile) detes and genus Bacteroides is the most abundant of are pathogenic, and an enrichment of this genus has Velasco-Galilea et al. Animal Microbiome (2020) 2:40 Page 10 of 16 previously been described in rabbits affected by epizootic Bacteroidales (phylum Bacteroidetes), and genera Desulfo- rabbit enteropathy [24]. This genus, together with genus vibrio (phylum Proteobacteria) and Bacteroides (phylum Bacteroides, was found enriched in the cecal microbiota Bacteroidetes). It is worth mentioning that these 5 OTUs of mice housed in open cages compared with those kept were also declared as differentially represented between in individual ventilated cages [25]. Both genera have farms by the univariate analyses. been associated with an exacerbation of the intestinal in- flammatory response in mammals [26]. Genus Anaero- fustis (family Eubacteriaceae) has been found enriched Administration of antibiotics impact on some taxa in cecal samples of rabbits affected by paratuberculosis relative abundances infection (Mycobacterium avium) [27]. Within farm B, the effect of the presence of antibiotics Within the phylum Verrucomicrobia, genus Akker- in the feed was assessed by comparing the microbial mansia is an anaerobic Gram-negative bacterium that cecal composition of rabbits fed with antibiotics with encompasses mucin degrader species [28]. In the cecum, that of some animals that received feed without antibi- a proper enrichment of this genus could maintain a suit- otics. As stated above, we did not detect significant dif- able mucosal turn-over, thus exerting a protective effect ferences in alpha diversity, nor in genera relative that could help the animal to deal with inflammatory abundances, between both groups. However, some sig- processes. nificant differences were observed at phylum and OTU It is worth mentioning that we have detected genera levels. An overrepresentation of phylum Cyanobacteria Epulopiscium, p-75-a5, Phascolarctobacterium, Cam- was found in rabbits fed without antibiotics. The detec- pylobacter and Desulfovibrio only in the cecal samples of tion of this bacterial phylotype, commonly assigned to rabbits housed in farm A. The first three are encom- photosynthetic activity, in the rabbit cecum could sug- passed within the phylum Firmicutes. Genus Epulopis- gest contamination during the GIT sampling. However, cium is a large size Gram-positive bacterium that has a Zeng et al. 2015 [31] previously reported its presence in nutritional symbiotic relationship with surgeonfish that rabbit feces. In the present study, all OTUs taxonomic- eats algae and detritus. This bacterium is physically simi- ally assigned to phylum Cyanobacteria are as well lar to the phylogenetically related Metabacterium poly- encompassed in the order YS2. Interestingly, it was dem- spora which is an endospore-producing bacterium onstrated that this order does not really have photosyn- isolated from the cecum of guinea pigs [29]. On the thetic capacity and it is currently classified within the other hand, genera Campylobacter and Desulfovibrio are candidate phylum Melainabacteria [32]. The non- Gram-negative bacteria that belong to phylum Proteo- photosynthetic cyanobacteria YS2, now named Gastra- bacteria. Some species of these genera are pathogens re- naerophilales, is a fermenter gut-associated order sponsible for infections and diarrheas in mammals. The present in humans and other animals such as squirrels, exclusive presence of these genera in farm A could indi- where its exact role is unknown but it has the capacity cate the existence of a potential dysbiosis of the animals to produce hydrogen, fix nitrogen and synthesize vita- raised in that facility that could affect their sanitary sta- mins B and K [32–34]. Our results, in accordance with tus and growth. While farm A was a semi-open-air facil- Kylie et al. 2018 [13], revealed that rabbits fed without ity, farm B was artificially ventilated and offered more antibiotics exhibited higher abundances of OTUs controlled environmental conditions that favor animal assigned to phylum Bacteroidetes than those fed with an- growth. Moreover, the presence of sulfate-reducing bac- tibiotics. In addition, samples of rabbits that received an- teria (SRB) such as Desulfovibrio could be enhanced by tibiotics had a significant increase of an OTU sulfate-secreting bacteria (SSB) such as Rikenella in farm taxonomically assigned to genus Coprococcus. Interest- A where this genus is significantly more predominant. It ingly, a study that evaluated the differences in bacterial is noteworthy to mention that SRB could also obtain sul- communities of Rex rabbits fed with different antibiotics fate via cross-feeding mediated by Bacteroides-encoded also found an overrepresentation of this bacterium in sulfatases [30], and interestingly, this phylum is more animals treated with zinc bacitracin [5]. Coprococcus is prevalent in farm A. an anaerobic bacterium that may protect against colon Regarding sample classification based on the sPLS-DA cancer in humans by producing butyric acid [35]. We study, given the important differences in gut microbial hypothesized that the administration of antibiotics could composition found between farms, a perfect classification modulate the abundance of some Coprococcus species to of the samples can be achieved with only 5 OTUs. One of provide intestinal protection on meat rabbits. However, these 5 OTUs was overrepresented in farm B and it is important to recognize that the reduced sample size belonged to family S24–7 (phylum Bacteroidetes). The of the group of rabbits fed without antibiotics may have remaining 4 were overrepresented in farm A and belonged limited the statistical power to detect microbial compos- to family Barnesiellaceae (phylum Bacteroidetes), order ition differences associated with this factor. Velasco-Galilea et al. Animal Microbiome (2020) 2:40 Page 11 of 16 Feed restriction modify Euryarchaeota and some bacteria Firmicutes and it is known as an obesity-associated bac- relative abundances terium [41] and as a pathogenic candidate identified in Within farm B, the effect of the feeding regime in micro- mice with multiple sclerosis [42]. A potential pro- bial composition was also assessed by comparing sam- inflammatory role has been proposed for this genus [42] ples of animals fed R with those fed AL. The main what could be related to a reduced incidence of enteric difference found was for phylum Euryarchaeota which disorders when feed restriction is applied. It is worth was overrepresented in animals fed R in farm B. All Eur- mentioning that family Peptococcaceae, which encom- yarchaeota species found in the rabbit cecum belong to passes genus rc4–4, is strongly related to total rabbit genus Methanobrevibacter that encompasses different weight gain from weaning to 12-week old [43]. Although hydrogenotrophic methane-producing species. Previous in our study this genus was prevalent in animals fed AL, studies in humans [36] and cattle [37, 38] found an over- its association with weight gain is not clear since the representation of Methanobrevibacter species in individ- greater growth exhibited by these animals was conse- uals submitted to feed restriction and a negative quence of higher feed intake. correlation between the abundance of this bacterium and body mass index. A prevalence of Methanobrevibac- Rabbit cecal microbiota is shaped by farm environment ter species could be a positive indicator of a healthy and animal management microbiota since restricted animals showed an overrep- Different approaches have been applied in this study to resentation of this genus. The main purpose of applying evaluate the effect of different environments and man- feed restriction is to improve intestinal health, reducing agement practices, commonly used in rabbit production, weaning mortality. The growth of Methanobrevibacter is in their cecal microbial composition and diversity. Those supported by fermenters such as Gastranaerophilales animals raised in the best insulated facility (farm B) ap- and butyrate-producing bacteria such as Anaereostipes pear to have a microbiota characteristic of healthier ani- via interspecies formate/hydrogen transfer [39]. A study mals than those raised in the open-air facility (farm A). in mice determined that Methanobrevibacter smithii fa- It is worth mentioning that the rabbits were housed in cilitates Bacteroides thetaiotaomicron capacity to digest cages interspersed with feeding regime. This fact could glycans resulting in increased production of short-chain make possible the exchange of microorganisms between fatty acids [40]. The same study defined M. smithii as a animals of different feeding regimes and therefore have “power broker” that regulates polysaccharide fermenta- reduced the differences observed between regimes. How- tion efficiency that influences the fat stores. The lower ever, the joint consideration of 70 OTUs in the sPLS-DA prevalence of methanogenic archaea in farm A could be made possible a certain discrimination power of samples explained by the high presence of SRB that outcompete according to the level of feeding received by each animal with methanogens for hydrogen consumption. This fact raised in farm A. It implies the existence of cecal micro- could favor hydrogen sulfide production and comprom- biota content patterns characteristic of each regime ise the rabbits’ health. which could be revealed thanks to the univariate ana- Regarding the sample classification based on the sPLS- lyses conducted at different taxonomic levels. Similarly, DA study conducted within farm B, component 1 and the sPLS-DA performed within farm B also involved the component 2 discriminated between animals that did or consideration of 70 OTUs to discriminate samples ac- did not received antibiotics in the feed and between cording to the amount of feed consumed. Within this feeding regimes, respectively. It is worth mentioning that farm, the classification of samples regarding the presence 8 out of 9 OTUs selected in component 1 were also de- or the absence of antibiotics in the feed needed a smaller clared as differentially represented between the presence number of OTUs than the feeding regime but greater or the absence of antibiotics in the feed by the univariate than the farm. This suggests that the effect of the pres- analyses. Within farm A, an sPLS-DA was also per- ence of antibiotic in feed is stronger than the feeding formed to classify samples according to the feeding re- level. The lack of a group of samples collected from ani- gime using microbial information. Although a large mals that did not receive antibiotics precluded the evalu- number of OTUs were selected as classifier variables in ation of the magnitude of importance of this effect over the tuning process of this sPLS-DA, the classification the feeding level on the cecal microbiota of animals error rate was high. It implied a poor discrimination raised farm A. It might have been possible that the mag- capacity of samples according to the feeding regime the nitude of the effect of the presence of antibiotics in the animal received. Nevertheless, bootstrap univariate ana- feed was larger in farm A than that observed in farm B. lyses of variance detected some significant differences at The experimental design of this study prevented the all taxonomic levels analyzed between feeding regimes comparison of the effect of antibiotic treatments across within farm A. At genus level, rc4–4 was overrepre- farms on rabbits’ microbial communities. The implica- sented in animals fed AL. This genus belongs to phylum tion of the discussed microbial composition and Velasco-Galilea et al. Animal Microbiome (2020) 2:40 Page 12 of 16 diversity differences originated by the studied manage- those in farm B received oxytetracycline, valnemulin and ment and environmental factors on the animals’ per- neomycin. At the time this experiment was conducted, it formance still needs to be investigated. In future studies was possible to use up to four types of molecules to pre- the role of specific groups of bacteria in rabbit growth vent or treat the emergence of potential infectious dis- and feed efficiency will be analyzed. eases on farms. However, nowadays, only one antibiotic molecule is allowed. During the last fattening week all Conclusions the animals received an antibiotic free diet. Feed was The analysis of a large number of animals from a paternal supplied once per day in a feeder with three places for rabbit line has allowed a deeper comprehension of the role the 4–5 weeks that the fattening lasted. Water was pro- played by different management and environmental fac- vided ad libitum during the whole fattening period. The tors shaping the composition and diversity of cecal micro- animals were under two different feeding regimes: (1) ad bial communities. It reveals that the farm environment libitum (AL) or (2) restricted (R) to 75% of the AL feed offered to the rabbits during their growth play a key role intake. The amount of feed supplied to the animals that can result in different microbial alpha diversity and under R feeding regime in a given week for each batch composition of almost all species that inhabit the rabbit was computed as 0.75 times the average feed intake of GIT. This highlights the importance that a stable and con- kits on AL from the same batch during the previous trolled environment could have in the intestinal health week, plus 10% to account for a feed intake increase as and, consequently, in animal performance. It seems clear the animal grows. Kits were randomly assigned to one of that the better insulated conditions of farm B favored the these two feeding regimes after weaning (32 days of age). presence of a gut microbiota characteristic of healthier an- They were categorized into two groups according to imals. Although the level of feeding and the presence of their size at weaning (big if their body weight was antibiotics in the feed did not modify the global diversity greater than 700 g or small otherwise) aiming to obtain of cecal microbial communities, these factors can increase homogenous groups regarding animal size within feed- or decrease the prevalence of specific bacteria which could ing regime. A maximum of two kits of the same litter lead to a microbial composition potentially beneficial for were assigned to the same cage in order to remove the the animal or, at the other extreme, to an origin of future possible association between cage and maternal effects intestinal dysbiosis. on animal growth during the fattening period. The dis- tribution of these animals across the different levels of Methods management factors is shown in Table 3. The body Animals and experimental design weight of each animal was weekly recorded. The individ- All biological samples used in the study were collected ual average daily gain was computed as the slope of the from animals of an experiment conducted at the Insti- within animal regression of all body weight measure- tute of Agrifood Research and Technology (IRTA) in dif- ments recorded during the growing period. ferent periods and involving two different farms. The objective of that experiment was to estimate the effect of the interaction between the genotype and the feeding re- Sample processing, DNA extraction and sequencing gime (i.e., the amount of feed provided during fattening) Animals were slaughtered (at 66 and 60 days of age in on growth, feed efficiency, carcass characteristics, and farm A and farm B, respectively) and cecal sample of health status of the animals [44]. For this particular each rabbit were collected in a sterile tube, kept cold in study, 425 meat rabbits from Caldes line [45] of that ex- the laboratory (4 °C) and stored at − 80 °C. DNA extrac- periment were randomly selected. Most of them (336) tion, amplification, Illumina library preparation and se- were raised in 4 different batches in a semi-open-air fa- quencing followed methods described previously [12]. cility (farm A). The remaining animals (89) were pro- Whole genomic DNA was extracted from 250 mg of duced in a single batch in another facility under better each cecal samples using ZR Soil Microbe DNA Mini- controlled environmental conditions (farm B). Rabbits Prep™ kit (ZymoResearch, Freiburg, Germany) according raised in farm A were housed in collective cages con- to manufacturer’s instructions with the following modifi- taining 8 kits each one while those raised in farm B were cation: cecal samples were mechanically lysed in a housed in cages with 6 kits each one. All animals were FastPrep-24™ Homogenizer (MP Biomedicals, LLC, raised under the same management conditions and re- Santa Ana, CA, United States) at a speed of 1 × 6 m/s for ceived the same standard pelleted diet. Twenty-three 60 s facilitating an efficient lysis of archaea and bacteria rabbits raised in farm B received a diet free of antibiotics species. Integrity and purity of DNA extracts were mea- and the remaining sixty-six received the same diet but sured with Nanodrop ND-1000 spectrophotometer supplemented with antibiotics. Those raised in farm A equipment (NanoDrop products; Wilmington, DE, received oxytetracycline, valnemulin, and colistin while United States) according to Desjardins and Conklin’s Velasco-Galilea et al. Animal Microbiome (2020) 2:40 Page 13 of 16 Table 3 Distribution of rabbits in groups according to different management factors Farm Batch Feed Feeding regime Number of rabbits A 1 With antibiotics Ad libitum 27 A 1 With antibiotics Restricted 30 A 2 With antibiotics Ad libitum 35 A 2 With antibiotics Restricted 41 A 3 With antibiotics Ad libitum 61 A 3 With antibiotics Restricted 53 A 4 With antibiotics Ad libitum 57 A 4 With antibiotics Restricted 32 B 5 With antibiotics Ad libitum 32 B 5 With antibiotics Restricted 34 B 5 Without antibiotics Ad libitum 12 B 5 Without antibiotics Restricted 11 protocol [46]. All DNA extracts had adequate integrity Service (SGB) of the Autonomous University of Barce- and purity (absorbance ratio 260 nm/280 nm > 1.6) to lona (UAB). avoid PCR inhibition issues. A fragment of the 16S rRNA gene including the V4- Bioinformatic pipeline for OTU calling V5 hypervariable regions was amplified with F515Y/ Sequence processing was performed using QIIME soft- R926 primer combination (5′-GTGYCAGCMGCC ware (version 1.9.0) [48]. In a first step, the resulting GCGGTAA-3′, 5′-CCGYCAATTYMTTTRAGTTT-3′) paired-ended V4-V5 16S rRNA gene reads were assem- [47] and then re-amplified in a limited-cycle PCR reac- bled into contigs with the python script multiple_join_ tion to add sequencing adaptors and 8 nucleotide dual- paired_ends.py. Then the contigs were curated using the indexed barcodes of multiplex Nextera® XT kit (Illumina, script split_libraries.py with default parameters in order Inc., San Diego CA, United States) following manufac- to assign them to samples and to discard those with a turer’s instructions. The initial PCR reactions were per- low-quality (Q19 was the minimum acceptable quality formed for each sample using 12.5 μl 2x KAPA HiFi score). Chimeric sequences generated during the process HotStart Ready Mix, 5 μl forward primer, 5 μl reverse of DNA amplification were detected with a UCHIME al- primer and 2.5 μl template DNA (5 ng/ μl). The initial gorithm [49] and removed. The totality of filtered con- PCR conditions were as follows: initial denaturation for tigs were clustered into operational taxonomic units 3 min at 95 °C, 25 cycles of 30 s at 95 °C, 30 s at 55 °C (OTUs) with a 97% similarity threshold using the script and 30 s at 72 °C; and final extension for 2 min at 72 °C. pick_open_reference_otus.py with default parameters [50] The addition of indexes and sequencing adaptors to both that grouped, through a UCLUST algorithm [51], the se- ends of the amplified regions took place in a second quences against Greengenes reference database (version PCR by using 25 μl 2x KAPA HiFi HotStart Ready Mix, gg_13_5_otus) and also made a de novo clustering of 5 μl index i7, 5 μl index i5, 10 μl PCR Grade water and those that did not match the database. The generated 5 μl concentrated amplicons of initial PCR. The second OTU table was filtered at: (1) sample level: by discarding PCR conditions were as follows: initial denaturation for samples with less than 5000 final sequence counts and 3 min at 95 °C, 8 cycles of 30 s at 95 °C, 30 s at 55 °C and at (2) OTU level: by removing the doubleton ones. The 30 s at 72 °C; and final extension for 5 min at 72 °C. Final filtered OTU table contained the sequence counts of 963 libraries were cleaned up with AMPure XP beads, vali- OTUs for 425 samples. Taxonomic assignment of repre- dated by running 1 μl of a 1:50 dilution on a Bioanalyzer sentative sequences of each OTU defined (963) was con- DNA 1000 chip (Agilent Technologies, Inc., Santa Clara, ducted by mapping them to the Greengenes reference CA, United States) to verify their size, quantified by database gg_13_5_otus with the UCLUST consensus tax- fluorometry with PicoGreen dsDNA quantification kit onomy assigner (QIIME default parameters). The raw (Invitrogen, Life Technologies, Carlsbad, CA, United sequence data were deposited in the sequence read arch- States), pooled at equimolar concentrations and paired- ive of NCBI under the BioProject accession number end sequenced in 5 parallel plates in a Illumina MiSeq PRJNA524130. Metadata, the prefiltered and normalized 2 × 250 platform at the Genomics and Bioinformatics OTU tables, and corresponding taxonomic classifications Velasco-Galilea et al. Animal Microbiome (2020) 2:40 Page 14 of 16 are also included as Additional files 1, 2, 3 and 4, or greater than that obtained with the original dataset. respectively. P-values were corrected defining a false discovery rate (FDR) of 0.05 [55]. This bootstrap analysis of variance Models and statistical methods approach was also implemented in order to study the ef- In order to study differences in diversity and richness fect of the management factors on the relative abun- between rabbits grouped according to farm environment dance of bacteria at phylum and genus levels. and management that they received, two alpha diversity The value of the microbial information to classify sam- indexes (Shannon and the observed number of OTUs) ples into the three factors considered in our study was were computed from the OTU table rarified to 10,000 explored using multivariate techniques. In particular, sequences per sample with “phyloseq” R package [52]. sparse partial least squares-discriminant analysis (sPLS- The statistical method chosen to assess alpha diversity DA) [56] was used to find the combination of OTUs that differences between these groups of animals was an ana- allowed the best classification of cecal samples according lysis of variance that included a factor resulting from the to: (1) the farm where the animals were raised, (2) the combination of four factors (the farm where the animal feeding regime within farm A and (3) the combination was raised, the batch, the presence or the absence of an- of feeding regime and the presence or absence of antibi- tibiotics in the feed and the feeding regime). The signifi- otics in the feed for the animals raised in farm B. This cance threshold was set at 0.05 for type I error. approach was implemented through the R package Different approaches were considered to assess the in- “mixOmics” [57]. In a first step, the function fluence of the environments and management factors on tune.splsda() was used to select the optimal sparsity pa- microbial composition. A bootstrap analysis of variance rameters of the sPLS-DA model: the number of compo- was individually implemented for each OTU to test nents and the number of variables (OTUs) per whether it was differentially represented between the dif- component. For the tuning process, a 5-fold cross- ferent categories of the factors studied. This univariate validation repeated 10 times was performed one compo- analysis was conducted by normalizing the OTU table nent at a time, with a maximum of 4 components, on an with the cumulative sum scaling (CSS) method [53] and input grid of values that indicate the number of variables only for those OTUs which were detected in at least 5% to select on each component. The sparsity parameters of the samples and had a sum of its counts resulting in a were defined, based on the BER and centroids distance, frequency greater than 0.01% of the total sum of all and then included in the final sPLS-DA model. Samples OTUs counts across all samples. It was implemented by were represented on the first two components and col- fitting a model defined by the combination of the four ored according to their class (e.g., R or AL in the case of aforementioned factors by using lm() function in R [54]. the feeding regime) in a sample plot with the function Then, the differences between the CSS-normalized plotIndiv(). The performance of the sPLS-DA model was OTUs counts in the different levels of the studied factors assessed with a 5-fold cross-validation repeated 1000 were tested. The significance between the levels of the times that randomly split the data in training and valid- main factors: farm, presence of antibiotics in the feed ation sets. In this data partition, it was ensured that 20% and feeding regime was assessed using an F statistic. of the samples within each level of the discriminant fac- When the involved interaction terms were significant, tor were assigned to the validation set. Five different par- the contrasts of interest were studied nested within the titions were performed for each replicate to guarantee a levels of other interacting factors, i.e. feeding regime different sample distribution in each validation set. The were studied within farm levels. When the interaction sPLS-DA model with the sparsity parameters previously terms were not significant, the effects of the different defined was adjusted in the training set and its classifica- levels were averaged, i.e. the effects of the levels of the tion performance was assessed in the validation set using batches within farm A were averaged to present the ef- the overall and BER per class as criteria. The stability of fect associated with this farm. In the performed F tests, the OTUs selected on each component was also assessed instead of relying on the theoretical distribution of the in the cross-validation by computing the selection fre- statistic under the null hypothesis to define the p-values, quency of each variable across the replicates. they were empirically computed using bootstrap after 1000 permutations of the dependent variable with re- spect to the design matrix of factors in the model. The Supplementary Information use of bootstrapping enabled the hypothesis test to be The online version contains supplementary material available at https://doi. done without the necessity of assuming that data are org/10.1186/s42523-020-00059-z. normally distributed, which is an assumption that fails for OTUs counts. P-value was defined as the proportion Additional file 1. Metadata associated with the 425 rabbit cecal samples analyzed in this study. of bootstrap rounds having an F statistic value equal to Velasco-Galilea et al. Animal Microbiome (2020) 2:40 Page 15 of 16 Additional file 2. Prefiltered and unnormalized OTU table used for Consent for publication statistical analyses in this study. Not applicable. Additional file 3. Filtered and CSS-normalized OTU table used for statis- tical analyses in this study. Competing interests The authors declare that they have no competing interests. Additional file 4. Taxonomic assignments for all OTUs in Additional file 2. Author details Additional file 5. Relative abundances phyla table built from the 1 Institute of Agrifood Research and Technology (IRTA) – Animal Breeding collapse of the filtered and CSS-normalized OTU table at phylum level. and Genetics, E08140 Caldes de Montbui, Barcelona, Spain. 2Institute of Additional file 6. Relative abundances genera table built from the Agrifood Research and Technology (IRTA) - Integral Management of Organic collapse of the filtered and CSS-normalized OTU table at genus level. Waste, E08140 Caldes de Montbui, Barcelona, Spain. 3Animal Genomics Group, Centre for Research in Agricultural Genomics (CRAG) Additional file 7: Table S1. OTUs differentially represented between CSIC-IRTA-UAB-UB, Campus UAB, Catalonia, Spain. 4Unit of Animal Science, farms. Department of Animal and Food Science, Autonomous University of Additional file 8: Table S2. OTUs differentially represented between Barcelona, Barcelona, Spain. feeding regimes within farms. Additional file 9: Table S3. OTUs differentially represented between Received: 3 July 2020 Accepted: 16 October 2020 the presence and the absence of antibiotics in the feed within farm B. References Abbreviations 1. Gaskins HR. Immunological aspects of host/microbiota interactions at the AL: Ad libitum; BER: Balanced error rate; CSS: Cumulative sum scaling; intestinal epithelium. Gastrointestinal Microbiol. 1997;2:537–87. FDR: False discovery rate; GIT: Gastrointestinal tract; NGS: Next generation 2. Belkaid Y, Hand TW. Role of the microbiota in immunity and inflammation. sequencing; OTU: Operational taxonomic unit; PCR: Polymerase chain Cell. 2014;157(1):121–41. reaction; R: Restricted; sPLS-DA: Sparse partial least squares-discriminant ana- 3. Gouet PH, Fonty G. Changes in the digestive microflora of holoxenic∗ lysis; SRB: Sulfate-reducing bacteria (SRB); SSB: Sulfate-secreting bacteria (SRB) rabbits from birth until adulthood. Ann Biol Anim Biochim Biophys. 1979;19: 553–66. Acknowledgements 4. Abecia L, Fondevila M, Balcells J, Lobley GE, McEwan NR. The effect of We would like to thank Oscar Perucho, Josep Ramon and Carmen Requena medicated diets and level of feeding on caecal microbiota of lactating (staff of Unitat de Cunicultura, IRTA) for their contribution to data recording rabbit does. J Appl Microbiol. 2007;103(4):787–93. and animal care during the experiment. We also acknowledge Nicolas 5. Zou F, Zeng D, Wen B, Sun H, Zhou Y, Yang M, et al. Illumina Miseq Boulanger and Joana Ribes (Genomics and NGS Unit, CRAG) for their platform analysis caecum bacterial communities of rex rabbits fed with assistance in massive libraries preparation. The English revision of the different antibiotics. AMB Express. 2016;6(1):100. manuscript conducted by Mr. Roderick Cantlay-Hollis is also acknowledged. 6. Zhu Y, Sun Y, Wang C, Li F. Impact of dietary fibre: starch ratio in shaping caecal archaea revealed in rabbits. J Anim Physiol Anim Nutr. 2017;101(4): Authors’ contributions 635–40. JS, MP and OR conceived the experimental design. JS, OR, MP and MVG 7. Chen SY, Deng F, Jia X, Liu H, Zhang GW, Lai SJ. Gut microbiota profiling collected biological samples. MVG, OGR, MP, MG and AS processed the with differential tolerance against the reduced dietary fibre level in rabbit. samples in the laboratory. MVG processed and analyzed the sequencing Sci Rep. 2019;9(1):1–9. data, interpreted data, prepared figures and tables, and wrote the 8. Combes S, Michelland RJ, Monteils V, Cauquil L, Soulié V, Tran NU, et al. manuscript. JS and YRC helped analyzing the sequencing data. JS, MG, MP, Postnatal development of the rabbit caecal microbiota composition and MV and YRC helped interpreting the data, and wrote and revised the activity. FEMS Microbiol Ecol. 2011;77(3):680–9. manuscript. All authors read and approved the final manuscript. 9. Gidenne T, García J, Lebas F, Licois D. Nutrition and feeding strategy: interactions with pathology. Nutr Rabbit. 2010;20:179. Funding 10. Eshar D, Weese JS. Molecular analysis of the microbiota in hard feces from The experimental design of this work was conducted thanks to funding from healthy rabbits (Oryctolagus cuniculus) medicated with long term oral INIA project RTA2011–00064–00-00. This study was part of the Feed-a-Gene meloxicam. BMC Vet Res. 2014;10(1):62. project that received funding from the European Union’s H2020 program 11. Gidenne T, Combes S, Feugier A, Jehl N, Arveux P, Boisot P, et al. Feed under grant agreement no. 633531, and the Spanish project RTI2018- restriction strategy in the growing rabbit. 2. Impact on digestive health, 097610R-I00. MVG is a recipient of a “Formación de Personal Investigador growth and carcass characteristics. Animal. 2009;3(4):509–15. (FPI)” pre-doctoral fellowship from INIA, associated with the research project 12. Velasco-Galilea M, Piles M, Viñas M, Rafel O, González-Rodríguez O, RTA2014–00015-C2–01. YRC was funded by Marie Skłodowska-Curie grant (P- Guivernau M, et al. Rabbit microbiota changes throughout the intestinal Sphere) agreement no. 6655919 (EU). tract. Front Microbiol. 2018;9:2144. 13. Kylie J, Weese JS, Turner PV. Comparison of the fecal microbiota of Availability of data and materials domestic commercial meat, laboratory, companion, and shelter rabbits The raw sequence data were deposited in the sequence read archive of (Oryctolagus cuniculi). BMC Vet Res. 2018;14(1):143. NCBI under the accession number SRP186982 (BioProject PRJNA524130). 14. Larsen OFA, Claassen E. The mechanistic link between health and gut Metadata, the prefiltered and unnormalized OTU table, the filtered and CSS- microbiota diversity. Sci Rep. 2018;8(1):1–5. normalized OTU table and corresponding taxonomic assignments have all 15. O’Neil DS, Stewart CJ, Chu DM, Goodspeed DM, Gonzalez-Rodriguez PJ, been included as Additional files 1, 2, 3 and 4, respectively. Relative abun- Shope CD, et al. Conditional postnatal deletion of the neonatal murine dances phyla and genera table have also been included as Additional files 5 hepatic circadian gene, Npas2, alters the gut microbiome following and 6, respectively. OTUs differentially represented between the studied fac- restricted feeding. Am J Obstet Gynecol. 2017;217(2):218–e1. tors, their sequences and their assignment at the lowest taxonomic level 16. Chen J, Toyomasu Y, Hayashi Y, Linden DR, Szurszewski JH, Nelson H, et al. have been included as Additional files 7, 8 and 9. Altered gut microbiota in female mice with persistent low body weights following removal of post-weaning chronic dietary restriction. Genome Ethics approval and consent to participate Med. 2016;8(1):103. This study was carried out in accordance with the recommendations of the 17. Zarrinpar A, Chaix A, Yooseph S, Panda S. Diet and feeding pattern animal care and use committee of the Institute for Food and Agriculture affect the diurnal dynamics of the gut microbiome. Cell Metab. 2014; Research and Technology (IRTA). The protocol was approved by the 20(6):1006–17. committee of the Institute for Food and Agriculture Research and 18. Soler C, Goossens T, Bermejo A, Migura-García L, Cusco A, Francino O, et al. Technology (IRTA). Digestive microbiota is different in pigs receiving antimicrobials or a feed Velasco-Galilea et al. Animal Microbiome (2020) 2:40 Page 16 of 16 additive during the nursery period. PLoS One. 2018;13(5):e0197353. https:// 39. Bui TPN, Schols HA, Jonathan M, Stams AJ, de Vos WM, Plugge CM. Mutual doi.org/10.1371/journal.pone.0197353. metabolic interactions in co-cultures of the intestinal anaerostipes 19. Kumar S, Chen C, Indugu N, Werlang GO, Singh M, Kim WK, et al. Effect of rhamnosivorans with an acetogen, methanogen, or pectin-degrader antibiotic withdrawal in feed on chicken gut microbial dynamics, immunity, affecting butyrate production. Front Microbiol. 2019;10:2449. growth performance and prevalence of foodborne pathogens. PLoS One. 40. Samuel BS, Gordon JI. A humanized gnotobiotic mouse model of host– 2018;13(2):e0192450. https://doi.org/10.1371/journal.pone.0192450. archaeal–bacterial mutualism. Proc Natl Acad Sci. 2006;103(26):10011–6. 20. Fang S, Xiong X, Su Y, Huang L, Chen C. 16S rRNA gene-based association 41. Ziętak M, Kovatcheva-Datchary P, Markiewicz LH, Ståhlman M, Kozak LP, study identified microbial taxa associated with pork intramuscular fat Bäckhed F. Altered microbiota contributes to reduced diet-induced obesity content in feces and cecum lumen. BMC Microbiol. 2017;17(1):162. upon cold exposure. Cell Metab. 2016;23(6):1216–23. 21. Dai ZL, Wu G, Zhu WY. Amino acid metabolism in intestinal bacteria: links 42. Gandy KAO, Zhang J, Nagarkatti P, Nagarkatti M. The role of gut microbiota between gut ecology and host health. Front Biosci. 2011;16(1):1768–86. in shaping the relapse-remitting and chronic-progressive forms of multiple 22. Sun X, Shen J, Liu C, Li S, Peng Y, Chen C, et al. L-Arginine and N- sclerosis in mouse models. Sci Rep. 2019;9(1):1–17. carbamoylglutamic acid supplementation enhance young rabbit growth 43. North MK, Dalle Zotte A, Hoffman LC. Composition of rabbit caecal and immunity by regulating intestinal microbial community. Asian Australas microbiota and the effects of dietary quercetin supplementation and sex J Anim Sci. 2020;33(1):166. thereupon. World Rabbit Sci. 2019;27(4):185–98. 23. Müller VM, Zietek T, Rohm F, Fiamoncini J, Lagkouvardos I, Haller D, et al. 44. Piles M, Sánchez JP. Use of group records of feed intake to select for feed Gut barrier impairment by high-fat diet in mice depends on housing efficiency in rabbit. J Anim Breed Genet. 2019;136(6):474–83. conditions. Mol Nutr Food Res. 2016;60(4):897–908. 45. Gómez EA, Rafel O, Ramon J. The caldes strain (Spain). Options 24. Bäuerl C, Collado MC, Zúñiga M, Blas E, Pérez Martínez. Changes in cecal méditerranéennes. Série B: Etudes et recherches (CIHEAM); 2002. microbiota and mucosal gene expression revealed new aspects of epizootic 46. Desjardins P, Conklin D. NanoDrop microvolume quantitation of nucleic rabbit enteropathy. PLoS One. 2014;9(8):e105707. https://doi.org/10.1371/ acids. J Vis Exp. 2010;45:2565. Published 2010 Nov 22. https://doi.org/10. journal.pone.0105707. 3791/2565. 25. Thoene-Reineke C, Fischer A, Friese C, Briesemeister D, Göbel UB, 47. Parada AE, Needham DM, Fuhrman JA. Every base matters: assessing small Kammertoens T, et al. Composition of intestinal microbiota in immune- subunit rRNA primers for marine microbiomes with mock communities, deficient mice kept in three different housing conditions. PLoS One. 2014; time series and global field samples. Environ Microbiol. 2016;18(5):1403–14. 9(11):e113406. https://doi.org/10.1371/journal.pone.0113406. 48. Caporaso JG, Kuczynski J, Stombaugh J, Bittinger K, Bushman FD, Costello 26. Terán-Ventura E, Roca M, Martin MT, Abarca ML, Martinez V, Vergara P. EK, et al. QIIME allows analysis of high-throughput community sequencing Characterization of housing-related spontaneous variations of gut data. Nat Methods. 2010;7:335–6. microbiota and expression of toll-like receptors 2 and 4 in rats. Microb Ecol. 49. Edgar RC, Haas BJ, Clemente JC, Quince C, Knight R. UCHIME improves 2010;60(3):691–702. sensitivity and speed of chimera detection. Bioinformatics. 2011;27:2194– 27. Arrazuria R, Elguezabal N, Juste RA, Derakhshani H, Khafipour E. 200. Mycobacterium avium subspecies paratuberculosis infection modifies gut 50. Rideout JR, He Y, Navas-Molina JA, Walters WA, Ursell LK, Gibbons SM, et al. microbiota under different dietary conditions in a rabbit model. Front Subsampled open-reference clustering creates consistent, comprehensive Microbiol. 2016;7:446. OTU definitions and scales to billions of sequences. PeerJ. 2014;2:e545. 28. Belzer C, De Vos WM. Microbes inside—from diversity to function: the case 51. Edgar RC. Search and clustering orders of magnitude faster than BLAST. of Akkermansia. ISME J. 2012;6(8):1449–58. Bioinformatics. 2010;26:2460–1. 29. Angert ER, Brooks AE, Pace NR. Phylogenetic analysis of Metabacterium 52. McMurdie PJ, Holmes S. phyloseq: an R package for reproducible interactive polyspora: clues to the evolutionary origin of daughter cell production in analysis and graphics of microbiome census data. PLoS One. 2013;8(4): Epulopiscium species, the largest bacteria. J Bacteriol. 1996;178(5):1451–6. e61217. 30. Rey FE, Gonzalez MD, Cheng J, Wu M, Ahern PP, Gordon JI. Metabolic niche 53. Paulson JN, Stine OC, Bravo HC, Pop M. Differential abundance analysis for of a prominent sulfate-reducing human gut bacterium. Proc Natl Acad Sci. microbial marker-gene surveys. Nat Methods. 2013;10(12):1200. 2013;110(33):13582–7. 54. R Development Core Team. R: A language and environment for statistical computing. 2010 [http://cran.r-project.org]. 31. Zeng B, Han S, Wang P, Wen B, Jian W, Guo W, et al. The bacterial 55. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical communities associated with fecal types and body weight of rex rabbits. Sci and powerful approach to multiple testing. J R Stat Soc Ser B Methodol. Rep. 2015;5:9342–50. 1995;57(1):289–300. 32. Di Rienzi SC, Sharon I, Wrighton KC, Koren O, Hug LA, Thomas BC, et al. The 56. Lê Cao KA, Rossouw D, Robert-Granié C, Besse P. A sparse PLS for variable human gut and groundwater harbor non-photosynthetic bacteria selection when integrating omics data. Stat Appl Genet Mol Biol. 2018;7(1). belonging to a new candidate phylum sibling to cyanobacteria. Elife. 2013; 57. Rohart F, Gautier B, Singh A, Le Cao KA. mixOmics: An R package for 'omics 2:e01102. feature selection and multiple data integration. PLoS Comput Biol. 2017; 33. Monchamp ME, Spaak P, Pomati F. Long term diversity and distribution of 13(11):e1005752. non-photosynthetic cyanobacteria in peri-alpine lakes. Front Microbiol. 2019; 9:3344. 34. Liu PY, Cheng AC, Huang SW, Lu HP, Oshida T, Liu W, et al. Body-size Publisher’s Note scaling is related to gut microbial diversity, metabolism and dietary niche of Springer Nature remains neutral with regard to jurisdictional claims in arboreal folivorous flying squirrels. Sci Rep. 2020;10(1):1–12. published maps and institutional affiliations. 35. Ai D, Pan H, Li X, Gao Y, Liu G, Xia LC. Identifying gut microbiota associated with colorectal cancer using a zero-inflated lognormal model. Front Microbiol. 2019;10:826. 36. Shen Q, Maitin V. Chapter 11 - Obesity-associated gut microbiota: characterization and dietary modulation. In: Diet-Microbe Interactions in the Gut. San Marcos: Academic Press; 2015. p. 149–71. https://doi.org/10.1016/ B978-0-12-407825-3.00011-3. 37. McCabe MS, Cormican P, Keogh K, O’Connor A, O’Hara E, Palladino RA, et al. Illumina MiSeq phylogenetic amplicon sequencing shows a large reduction of an uncharacterised Succinivibrionaceae and an increase of the Methanobrevibacter gottschalkii clade in feed restricted cattle. PLoS One. 2015;10(7):e0133234. https:// doi.org/10.1371/journal.pone.0133234. 38. McGovern E, McCabe MS, Cormican P, Popova M, Keogh K, Kelly AK, et al. Plane of nutrition affects the phylogenetic diversity and relative abundance of transcriptionally active methanogens in the bovine rumen. Sci Rep. 2017; 7(1):1–10.