Effects of particle size of ground alfalfa hay on caecal bacteria and archaea populations of rabbits

Animal experiment and sample collection

A total of 120 New Zealand rabbits (half male and female) with average body weight (950.3  ± 8.82 g), weaned at 35 days of age were selected and allocated into 4 treatments, with five replicates in each treatment and six rabbits in each replicate, and the rabbits were kept in three pairs with two rabbits per cage. Rabbits were housed in the same building (Sichuan Agricultural University, China) in flat-deck cages (600 × 250 × 330 mm) for the 49 d experiment. The temperature inside the house was maintained at 15–25 °C and rabbits had ad libitum access to water and diets during the whole experimental period. Neither feed nor drinking water was medicated with antibiotics, but a coccidiostatic (robenidine) was provided in the feed. In this study, the diets (Table S1) were formulated according to NRC (1977) Growing Rabbit Feeding Standards, and the alfalfa meal incorporated into the four diets were 2,500, 1,000, 100 and 10 µm, respectively. Firstly, alfalfa meal with particle sizes of 2,500, 1,000, 100 and 10 µm were produced, and then the rest of the ingredients were milled through a 2.5-mm grinder screen. After all the ingredients were ready, they were mixed and granulated (diameter was 3 mm). After the finish of the growth experiment, 24 animals were slaughtered (six rabbits per diet) by cervical dislocation 1 h before dark (19:00) to avoid soft feces excretion. Once slaughtered, 50 g of cecal content was collected in sterile conditions. The samples were immediately frozen at −80 °C until DNA extraction.

DNA extraction

Caecal samples were thawed at room temperature and kept on ice during the extraction process. According to the manufacturer’s instructions, total genomic DNA was extracted from caecal samples using a DNeasy PowerSoil Kit (Qiagen, Valencia, CA, USA) that included a step of the microbiological cells to be mechanically broken. The quality and concentration of the extracted DNA was detected, respectively, using 0.8% agarose gel electrophoresis and a NanoDrop ND-1,000 spectrophotometer (Nyxor Pharmacia, Paris, France). All extracted DNA samples were diluted to 10 ng/µL using sterile ultrapure water and stored at −20 °C until used for real-time PCR and Illumina sequence analyses.

PCR amplification, Illumina library generation and sequencing

The V4 variable region of 16S rRNA gene was amplified by using caecal total DNA as template. Selecting 515F (5′-GTGCCACMCCGCGGTAA-3′) and 806R (5′-GGACTACHVGGGTWTCTAAT-3′) as the primer pair for the bacteria (Caporaso et al., 2011), while selecting A516F (5′-TGYCAGCCGCCGCGGTAAHACCVGC-3′) and U806R (5′-GGACTACHVGGGTWTCTAAT-3′) as the primer pair for the archaea (Kuroda et al., 2015). 16S rRNA genes were amplified using the specific primer with 12 nt unique barcode, and used the same reaction system for PCR amplification for bacteria and archaea. The total PCR mixture (25 µL) contained 1 × PCR buffer, 1.5 mM MgCl₂, each deoxynucleoside triphosphate at 0.4 µM, each primer at 1.0 µM, 0.5 U of KOD-Plus-Neo (Toyobo, Tokyo, Japan) and 10 ng template DNA. The PCR amplification program consists of initial denaturation at 94 °C for 1 min, followed by 30 cycles of denaturation at 94 °C for 20 s, annealing at 54 °C for 30 s, elongation at 72 °C for 30 s, and a final extension at 72 °C for 5 min. Three replicates of PCR reactions for each sample were combined together, PCR products mixed with 1/6 volume of 6× loading buffer were loaded on 2% agarose gel for detection. We used a 2 step process for library preparation. First, 16S amplicons were generated using PCR to generate 200–400 bp amplicons. Second, we added unique barcodes to samples using emulsion PCR (EmPCR) to prevent chimera formation (Williams et al., 2006), and the PCR was stopped at linear stage. This methodological basic principle is dilution and compartmentalization of template molecules in water droplets in a water-in-oil emulsion. Ideally, the dilution is to a degree where each droplet contains a single template molecule and functions as a micro-PCR reactor (Kanagal-Shamanna, 2016). This method reduces the formation of artifactual molecules, as often seen in conventional PCR, thus preserving library complexity (Williams et al., 2006; Qiu et al., 2001). The electrophoresis band was purified using OMEGA Gel Extraction Kit (Omega Bio-Tek, USA), the gel purified barcoded amplicons were pooled with equal molar amount and quantified on a Qubit@ 2.0 Fluorometer (Thermo Scientific). Finally, 15% of PhiX control library was spiked into the amplicon pool to improve the unbalanced and biased base composition. In brief, sequencing libraries were generated using TruSeq DNA PCR-Free Sample Prep Kit following manufacturer’s recommendations and index codes were added, and using ZymoBIOMICS Microbial Community Standard (Cat#D6300) as the positive control. The library quality was assessed on the Qubit@ 2.0 Fluorometer (Thermo Scientific) and Agilent Bioanalyzer 2100 system. At last, the library was applied to paired-end sequencing (2 × 250 bp) with the Hiseq Illumina Sequencing Platform (Rhonin Biosciences Co., Ltd., Chengdu, China), according to the protocols described by Caporaso et al. (2012).

Bioinformatics analysis

The sequences were analyzed according to Usearch (version 7.1 http://drive5.com/uparse/) and QIIME (Caporaso et al., 2010a) pipeline. Paired-end reads from the original DNA fragments were merged using FLASH (Magoč & Salzberg, 2011). Then sequences were assigned to each sample according to the unique barcode. The first step was to filter low quality reads (length < 200 bp, more than two ambiguous base ‘N’, or average base quality score < 30) and truncated sequences where quality scores decay (score < 11) using Trimmomatic (Bolger, Lohse & Usadel, 2014) and Usearch. After finding duplicated sequences, discarded all the singletons, which may be sequencing errors (http://www.drive5.com/usearch/manual/singletons.html) and lead to false positive results for overestimation of diversity. Sequences were clustered into operational taxonomic units (OTUs) at 97% identity threshold using UPARSE algorithms (Edgar, 2013), and picked representative sequences and removed potential chimeras using Uchime algorithm (Edgar et al., 2011). Taxonomy were assigned using the Silva (https://www.arb-silva.de) database (Quast et al., 2013) and uclust classifier in QIIME (version 1.8.0). Representative sequences were aligned using PyNAST (Caporaso et al., 2010b) embedded in QIIME. After quality checking, phylogenetic trees were reconstructed based on maximum likelihood–approximation method using the generalised time-reversible (GTR) model in FastTree (Price, Dehal & Arkin, 2010).

Then the filtered OTU is removed from the OTU table in the process of evolutionary tree reconstruction, and the OTU table is resampled so that each sample has the same sequence number. Phylogenetic diversity (Faith’s PD18) was calculated using Picante (Kembel et al., 2010). Weighted and Unweighted Unifrac distances were calculated in GUniFrac version 1.1 (Chen et al., 2012). The analysis of alpha and beta diversity metrics was conducted with Vegan (version 2.0-2. R CRAN package). Rarefaction curves were generated based on these three metrics. Principal component analysis (PCA) was applied to reduce the dimensions of original community data. Hierarchical cluster analysis was done using R (https://www.R-project.org/) function hclust. To identify if there were significant differences among different groups, permutational multivariate analysis of variance was performed based on the dissimilarity matrix.

Statistical analysis

The effects of alfalfa meal particle size on cecal microbial communities were analyzed using the MIXED procedure of SAS v9.4 (SAS Institute Inc., Cary, NC, USA). The model used for the analysis was ${\displaystyle Y_{ij}=\mu +T_i+e_{ij}}$ where Y_ij is an observation on the dependent variable ij, μ is the overall population mean, T_i isthe fixed effect of alfalfa meal particle size, and e_ij is the random error associated with the observation ij. Tukey-Kramer multiple comparison tests were performed after differences were detected. Differences between means with P < 0.05 were accepted as statistically significant differences.

Sequencing depth, coverage and alpha diversity

A rarefaction test was performed at the OTU level and the results are presented (Fig. S1). As can be seen from the rarefaction curve that all of the samples tended to reach a plateau, indicating that the sequencing quantity of both bacteria and archaea covered most microorganisms. Measures of alpha diversity were shown in Table 1. As for the bacteria, the OTUs, Chao1 index, Shannon-Wiener index and PD value (phylogenetic diversity) were all numerically higher in group 10 µm than that of the other three groups (P > 0.05). As for the archaea, the Shannon-Wiener index of group 1,000, 100 and 10 µm were significantly decreased compared with group 2,500 µm (P = 0.044). However, no significant difference was obtained among 1,000, 100 and 10 µm groups. The above results suggested that rabbit caecum alpha diversity of bacteria and archaea experienced different alterations when the alfalfa particle size decreased.

A total of 745,946 bacterial sequences were generated after quality control with a mean of 31,081  ± 13,901 (mean  ± standard deviation [SD], n = 24) per sample, and the mean length of tag N50 is 304  ± 6 base pairs (bp). Based on the principle that the similarity is greater than 97%, the obtained effective sequences were clustered into a total of 9,953 OTUs, the average OTU was 1,896  ± 251 (mean  ± SD, n = 24) per sample. A total of 2107 OTUs were shared in all of four treatments, while 1,390, 587, 632 and 1,791 OTUs were exclusively found in group 2,500, 1,000, 100 and 10 µm, respectively (Fig. 1A). Illumina HiSeq sequence of the archaeal 16S rRNA yielded 539,227 sequences for the 24 caecal samples, with a mean of 22,468  ± 2,443 sequences per sample, and the mean length of tag N50 is 292  ± 4. A total of 2246 OTUs were assigned, and each sample contained 336  ± 21 OTUs. Furthermore, 351 OTUs were shared by the four groups, while 404, 391, 378 and 398 OTUs were exclusively found in group 2,500, 1,000, 100 and 10 µm, respectively (Fig. 1B).

Analysis of microbial community similarity in caecum of rabbits

Community OTU comparisons were visualised by clustering analysis (OTU ≥ 97% identity, species level similarity) using weighted unifrac clustering in Fig. 2. The closer samples and the shorter branches are indicating the more similar the species composition of the two samples. The relationship of bacterial community between different treatments showed that the samples were clustered into three branches (Fig. 2A). All samples in 10 µm group clustered into the same clade, however, sample 2,500-2, 2,500-4 and 100-3 were clustered into this clade. The samples of group 2,500, 1,000 and 100 µm were blended together. This means that the bacterial composition changed when the particle size of alfalfa meal decreased from 100 to 10 µm. Similarly, a difference was also observed in the composition of archaea between the four treatments (Fig. 2B). The 2,500 and 1,000 µm groups were clustered thoroughly into one clade first and then mingled together, and the distance was larger than 0.2 between the two groups. The 100 and 10 µm groups mixed together and could not be distinguished, and the distance was less than 0.1. This means archaea is kept stable when the particle size of alfalfa meal larger than 1,000 µm, and it is altered when the particle size meal is lowered less than 100 µm. In addition, the archaea in the caecum of rabbit is changed earlier than bacteria as alfalfa meal particle size decreased.

A PCoA plot of overall diversity based on weighted UniFrac metric was generated (Fig. 3). The closer the distance between points means the more similar the community structure of the two samples. The PCoA plot demonstrated that the bacterial community between 10 µm group and other three treatments could be clearly distinguished (P < 0.01). However, the bacterial community among 2,500, 1,000 and 100 µm groups could not be obvious isolated (P > 0.05), and PCo1 accounted for 32.2% (Fig. 3A; Table S2). In addition, the six archaeal samples of group 100 and 10 µm were mixed together (P = 0.995) but could be separated thoroughly from group 2,500 and 1,000 µm (P = 0.025), and the PCo1 accounted for 98.7% (Fig. 3B; Table S2). And the same result was obtained in our unweighted UniFrac metric analysis (Fig. S3). This again, proved that the archaeal composition structure changed before the bacterial structure, and the uniformity of archaea is better with the decrease of alfalfa meal particle size.

Taxonomy of rabbit caecum microbial composition

The sequences in the present experiment were finally annotated as bacteria and archaea, a total of 26 bacterial phyla and 3 archaeal phyla are identified. As for the bacteria, among which 23 phyla were detected in the 2,500 and 100 µm groups, and 21 and 24 phyla were obtained in 1,000 and 10 µm group, respectively. Figure 4 presents the relative abundance of microbial composition at the phylum level of different treatments. There were only 3 bacterial phyla with relative abundance larger than 1% in 2,500, 1,000 and 100 µm groups (Fig. 4A), and the two most abundant phyla were Firmicutes (58.69 ± 0.114%, 67.77 ± 0.105% and 67.82 ± 0.097%) and Bacteroidetes (36.05 ± 0.118%, 26.24 ± 0.114% and 26.20 ± 0.105%), followed by Proteobacteria (2.69 ± 0.003%, 3.37 ± 0.005% and 3.40 ± 0.004%). At the meanwhile, the 10 µm group had 4 bacterial phyla with relative abundance above 1%. In addition to the three dominating phyla, Firmicutes (68.50 ± 0.052%), Bacteroidete (23.96 ± 0.056%) and Proteobacteria (3.34 ± 0.002%), the relative abundance of Tenericutes (2.08 ± 0.012%) was also greater than 1%. The archaeal abundance analysis showed that the highest relative abundance was Euryarchaeota (over 99.9%) in all of the four treatment groups (Fig. S2). The remaining archaeal abundance (Thaumarchaeota and Miscellaneous Crenarchaeotic Group) were less than 0.1%.

At the genus level, the bacteria kingdom is composed of 465 genera, whereas only 26 genera with relative abundance larger than 1%, among which 10 of them were unclassified. The number of genera with relative abundance greater than 1% in the four treatment groups were 19 (2,500 µm), 23 (1,000 µm), 22 (100 µm) and 20 (10 µm), respectively. The most abundant genera (relative abundance more than 4%) were Unclassified Bacteroidales S24-7 group (20.39 ± 0.112%), Ruminococcaceae UCG-014 (13.25 ± 0.101%), Lachnospiraceae NK4A136 group (10.73 ± 0.065%) and Rikenellaceae RC9 gut group (5.00 ± 0.066%) in the 2,500 µm group. Unclassified Bacteroidales S24-7 group (14.48 ± 0.087% and 14.56 ± 0.088%), Ruminococcaceae UCG-014 (9.34 ± 0.040% and 11.55 ± 0.097%), Lachnospiraceae NK4A136 group (8.80 ± 0.046% and 9.61 ± 0.035%), unclassified Clostridiales vadinBB60 group (5.72 ± 0.060% and 4.86 ± 0.015%) and Ruminococcaceae NK4A214 group (4.71 ± 0.024% and 4.39 ± 0.027%) were the most two abundant genera in 1,000 and 100 µm groups. In the 10 µm group, Ruminococcaceae UCG-014 (33.82 ± 9.752%) and Unclassified Bacteroidales S24-7 (11.81 ± 2.876%) group were the two most abundant genera. The relative abundance of the 10 unclassified genera accounted for 29.45%, 20.55%, 30.34% and 24.86% in 2,500, 1,000, 100 and 10 µm groups, respectively (Fig. 5A).

Table 2 shows the significantly affected bacteria by the decrease of alfalfa particle size. Ruminococcaceae UCG-014 (P < 0.001) and Lactobacillus (P = 0.043) were increased while the alfalfa particle size decreased, while the relative abundance of Lachnospiraceae NK4A136 group (P = 0.016), Ruminococcaceae NK4A214 (P = 0.044), Ruminococcaceae UCG-005 (P = 0.012), Christensenellaceae R-7 group (P = 0.019), Lachnospiraceae other (Family) (P = 0.011) and Ruminococcaceae UCG-013 (P = 0.021) were decreased .

Meanwhile, a total of 10 genera were assigned from the sequences of archaea, however only 4 genera were classified. Archaeal classification at the genus level demonstrated that the dominating genus in caecum of rabbits was Methanobrevibacter, and its relative abundance was significantly increased from 62.48% to 90.40% as the particle size of alfalfa decreased from 2,500 to 10 µm (P < 0.001). The Methanosphaera was the second largest genus in the caecum of rabbits, and its relative abundance decreased from 35.47% to 8.26% with the decrease of alfalfa particle size (P < 0.001) (Fig. 5B; Table 2).

Nucleotide sequence accession numbers

Sequences of this project have been deposited into the Sequence Read Archive (SRA) of the NCBI nucleotide database under accession number PRJNA542420.