Metagenomic analysis of rumen microbial population in dairy heifers fed a high grain diet supplemented with dicarboxylic acids or polyphenols.

BACKGROUND: The aim of this study was to investigate the effects of two feed supplements on rumen bacterial communities of heifers fed a high grain diet. Six Holstein-Friesian heifers received one of the following dietary treatments according to a Latin square design: no supplement (control, C), 60 g/day of fumarate-malate (organic acid, O) and 100 g/day of polyphenol-essential oil (P). Rumen fluid was analyzed to assess the microbial population using Illumina sequencing and quantitative real time PCR.
RESULTS: The P treatment had the highest number of observed species (P < 0.10), Chao1 index (P < 0.05), abundance based coverage estimated (ACE) (P < 0.05), and Fisher's alpha diversity (P < 0.10). The O treatment had intermediate values between C and P treatments with the exception of the Chao1 index. The PCoA with unweighted Unifrac distance showed a separation among dietary treatments (P = 0.09), above all between the C and P (P = 0.05). The O and P treatments showed a significant increase of the family Christenenellaceae and a decline of Prevotella brevis compared to C. Additionally, the P treatment enhanced the abundance of many taxa belonging to Bacteroidetes, Firmicutes and Tenericutes phyla due to a potential antimicrobial activity of flavonoids that increased competition among bacteria.
CONCLUSIONS: Organic acid and polyphenols significantly modified rumen bacterial populations during high-grain feeding in dairy heifers. In particular the polyphenol treatment increased the richness and diversity of rumen microbiota, which are usually high in conditions of physiological rumen pH and rumen function.

Background

The rumen is an anaerobic fermentation chamber housing diverse microorganisms such as bacteria, protozoa, fungi and viruses. Bacteria play a key role in rumen fermentation, which in turn greatly impact production and health of dairy cows. [1]. High-grain diets used to improve performance of dairy cows alter microbial communities in the rumen and the symbiosis between the host and these communities by causing the production of excessive amounts of organic acids (volatile organic acid (VFA) and lactate) [2-4]. High concentrations of VFA and lactate can exceed the buffering capacity of the rumen resulting in major changes in the rumen environment, such as a pH depression, and reduction in the populations of many beneficial bacteria [1, 5, 6]. Large individual variability of the rumen bacterial community, even in animals feeding on the same ration have been reported. [7, 8]. Other studies [1, 6, 7] showed that dietary changes, like an increased inclusion of cereals in the ration, lead to shifts in the microbial community even though the ruminal microbiota was able to maintain a stable core of bacterial taxa. The populations of starch-fermenting and lactic acid-utilizing bacteria increase when high starch diets are fed [9]. An excessive depression in rumen pH causes a reduction in the number of cellulolytic bacteria [2], major shifts in bacterial populations, and increases the lysis of Gram-negative bacteria resulting in a higher concentration of free lipopolysaccharide (LPS) endotoxins in rumen digesta [2, 10, 11]. When the translocation of free LPS from the rumen and/or lower gut into the interior circulation occurs, an immune response in the host is induced, causing metabolic alterations. Moreover, the increase in proportion of bacteria such as Enterobacteriacae during high starch feeding heightens the presence of virulence factors (fimbrial adhesins, heat-stable and heat-labile toxins and inflammatory peptides), which have the potential to cause inflammation [2, 10, 12]. Nutritional strategies to prevent the onset of SARA and translocation of LPS in dairy cows are based on feeding sufficient amounts of physically effective fiber (peNDF) and modulating the amount of easily degradable starch in the diet [13]. In addition the use of feed supplements to enhance the rumen microbial community and subsequently ruminal fermentation has also been suggested for this purpose. These supplements include the use of yeasts, like Saccharomyces cerevisiae, strain CNCM I-1077 [14], probiotic bacteria, like Megasphaera elsdenii, strain H6F32 [15], and dicarboxylic acids [16], flavonoids [17] and essential oils [18]. It has been suggested that the dicarboxylic acids malate and fumarate attenuate the ruminal pH drop during high grain feeding [16], and that these acids increase the activity of the succinate-propionate metabolic pathway in several rumen bacteria, resulting in increased lactic acid uptake and production of propionate [19, 20]. In in vitro rumen fermentative experiments, several essential oils (EO), or blends of EO have been demonstrated to enhance rumen fermentation [18, 21]. However, few in vivo studies have investigated the effects of EO on rumen fermentation and bacterial populations [22]. Furthermore, the effects of addition of polyphenolic compounds like flavonoids to diets of dairy cows may include prevention of the pH reduction and the decrease of the acetate-to-propionate ratio due to an increase of the numbers of lactate-consuming and propionate-producing bacteria [17].

A previous study described the effect of organic acid (O) and polyphenols (P) on reticular pH drop and acute phase response in dairy heifers fed a high grain diet [23]. In this manuscript we report the effects of these supplements on rumen bacterial populations in dairy heifers fed a high-grain diet.

Results

Illumina sequencing produced 317,369 sequences across treatments. The average numbers of sequences were 16,413, 18,814 and 28,247 for the C, O and P treatment, respectively. The number of the generated sequences was not affected by treatment (P = 0.60).

The P treatment led to the highest richness of the rumen microbial population (Table 1), characterized by the highest number of observed species (P < 0.10) and the highest value of ACE (P < 0.05). The O treatment resuted in intermediate values between P and C for both the number of species observed and the index ACE.

Table 1

Number of observed species, richness (Chao1 and ACE) and diversity estimators (Shannon, Simpson and Fisher’alfa) in control (C), organic acid (O) and polyphenols (P) dietary treatments

	Observed species	Richness		Diversity
		Chao1	ACE	Shannon	Simpson1	Fisher’s alpha
Control	2,116β	3,420b	3,731b	5.9	0.99	645β
Organic acid	2,255α, β	3,834b	4,214a, b	6.1	0.99	691α, β
Polyphenols	2,742α	7,164a	6,927a	6.1	0.98	890α
SEM	188.7	653.1	711.8	0.22	(0.99–0.98)	78.2
P-value	0.08	0.03	0.06	0.90	0.93	0.07

1Statistical analysis was conducted on natural logarithm (ln) transformed data that are presented as ln back transformed and 95 %-confidence interval in brackets

a, bWithin column indicate statistical differences (P < 0.05); α, βWithin column indicate statistical differences (P < 0.10)

The Chao1 index was higher during the P treatment (P < 0.05) compared to the C and O treatments. Moreover, Fisher’s alpha diversity tended (P < 0.10) to be higher during the P treatment compared to the C treatment, whereas the O resulted in intermediate values. The Shannon and Simpson indices were not affected by treatments.

The PERMANOVA analysis showed a significant (P < 0.10) effect of treatments on Bray-Curtis and unweighted Unifrac distances (Table 2). In addition, differences were observed in bacterial population among periods (P = 0.02 and P = 0.06, respectively). However the interaction of treatment and period (D × P) on these measures was not significant pointing out the lack of a carry-over effects of treatments. The multiple-comparisons based on Bray-Curtis distance showed only a tendency (P = 0.11) towards a significant difference between dietary treatments (C vs P), while the multivariate analysis based on unweighted Unifrac distance indicated a significant (P = 0.05) separation between the C and P treatments (Table 2). The comparison of the bacterial communities by principal coordinate analysis (PCoA) with unweighted Unifrac distance (Fig. 1) confirmed a separation between the C and P treatments.

Table 2

PERMANOVA analysis of the effect of dietary treatments, periods and heifers on rumen bacteria dissimilarities based on Bray-Curtis and unweighted Unifrac

Source	Bray-curtis		Unweighted unifrac
	pseudo-F	P-value	pseudo-F	P-value
Dietary treatment	1.208	0.10	1.155	0.09
Period	1.328	0.02	1.174	0.06
Heifer	1.175	0.11	1.135	0.09
Pair-wise tests
Control vs Organic acid		0.14		0.22
Organic acid vs Polyphenols		0.30		0.15
Control vs Polyphenols		0.11		0.05

P-values were calculated on 999 possible permutations

Fig. 1

Principal coordinate analysis (PCoA) using unweighted Unifrac to explore dissimilarities in microbial composition among dietary treatments (Control, C; Organic acids, O and Polyphenols, P). The axes (PC1 = 12.4 % and PC2 = 11.0 %) account for 23.4 % of the total variation of the model

The abundance of the 20 phyla detected by using Illumina was not affected by dietary treatment. Bacteroidetes, Firmicutes, Tenericutes and Euryarchaeota bacterial phyla were abundant (>1 %) and accounted for 94.8 % of the total bacterial community. Spirochaetes, Proteobacteria, Actinobacteria, Fibrobacteres, Chloroflexi, Cyanobacteria and Verrucomicrobia were in low abundance (0.1–1 %), and unclassified bacteria accounted for 2.3 % of the total bacterial community.

Illumina sequencing detected 32 classes, 55 orders, 87 families and 88 genera. The predominant sequences (>1 %) within Bacteroidetes belong to order Bacteroidales (22.7 %), genus Prevotella (21.2 %), genus Paludibacter (2.24 %), family BS11 (2.04 %), family S24-7 (1.77 %), and genus CF231 (1.06 %). The predominant sequences in Firmicutes were genus Butyrivibrio (4.30 %), and family Christensenellaceae (2.43 %), and in phylum Tenericutes they were order RF39 (1.74 %). In particular, the cellulolytic bacteria Ruminococcus, Treponema and Fibrobacter represented 8.15, 0.58 and 0.19 %, respectively on average of total bacteria. The relative abundance of the starch-fermenting bacteria were low and included: Bifidobacterium (0.43 %), Eubacterium (0.03 %), Selenomonas (0.02 %), Succinivibrio (0.02 %), and Ruminobacter (<0.01 %). The relative abundance of Butyrivibrio (4.30 %) was high, while Succimonas, Lactobacillus and Streptococcus were not detected. Among lactic acid-utilizing bacteria, Anaerovibrio (0.01 %) was detected, whilst Fusobacterium was not detected. The statistical analyses of these data (OTUs–Illumina), refering to the single taxa, showed a total of 19 significant pair-wise comparisons (O vs C, P vs C, and P vs O) as reported in Fig. 2. The comparison between O and C showed that only the family Christensenellaceae was significantly (P < 0.05) more abundant in the O treatment (Fig. 2a). The greatest number of significant differences between OTUs (11) were observed in the comparison between the P and C treatments. In particular, the P dietary treatment showed a significant higher number of OTUs in the case of phylum Bacteroidetes, order Bacteroidales (phylum Bacteroidetes), family BS11 (phylum Bacteroidetes), genus Paludibacter (phylum Bacteroidetes), genus YRC22 (phylum Bacteroidetes), genus CF231 (phylum Bacteroidetes), genus Butyrivibrio (phylum Firmicutes), family Christensenellaceae (phylum Firmicutes), family Mycoplasmataceae (phylum Tenericutes), and genus RFN20 (phylum Tenericutes). The opposite occurred for the family WCHB1-25 (phylum Verrucomicrobia) (Fig. 2b). Moreover, during the P treatment the abundance of genera YRC22 (phylum Bacteroidetes), CF231 (phylum Bacteroidetes), Anaeroplasma (phylum Tenericutes), and RFN20 (phylum Tenericutes), was increased compared to the O treatment; the contrary for the families Pseudomonadaceae (phylum Proteobacteria), WCHB1-25 (phylum Verrucomicrobia) and the order Actinomycetales (phylum Actinobacteria) (Fig. 2c).

Fig. 2

a, b, c. Relative changes (log2 fold) in rumen bacterial taxa (based on the classification of 16s RNA gene) in the comparison of a organic acid vs control, b polyphenols vs control, and c polyphenols vs organic acid determinated with Illumina. †: P < 0.10; *: P < 0.05. Bar = standard error

The comparison of single species based on qPCR and submitted to the linear MIXED model showed a significant effect of the treatments only in the case of Prevotella brevis. This species was more abundant (P = 0.03) during the C treatment compared to the other treatments (Fig. 3).

Fig. 3

Relative changes (log2 fold) in predominant rumen microorganisms in the comparison of a organic acid vs control, b polyphenols vs control, and c polyphenols vs organic acid determined with quantitative real-time PCR. *: P < 0.05. Bar = standard error

Discussion

This study was conducted to investigate the effects of the two feed supplements, i.e. organic dicarboxylic acids (O) and polyphenols (P), on the rumen microbiota in dairy heifers fed a high-grain diet. The effects of these additives on reticular pH, VFA, lactate and immune response was published earlier by De Nardi et al. [23]. In this companion study both additives, O and P, attenuated the reticular pH drop due to high-grain feeding wich resulted in 16, 18 and 199 min/d spent below pH 5.6, for O, P and control, respectively. In particular, the P treatment was more effective in reducing the inflammatory response that resulted from the high grain feeding without interfering with dry matter intake (DMI), whereas the O treatment caused the highest decline of the acetate to propionate ratio.

Both supplements enhanced the richness and α-diversity of the rumen microbiota compared to control. This is in agreement with other studies [1, 5, 22], which report a higher richness and diversity for cows with physiological rumen pH and rumen function compared to cows with impaired rumen parameters. This finding could have been related to the different ruminal pH conditions found in heifers fed P and O treatments, which spent daily less time below pH 5.6 [23] compared to heifers on the control treatment. However, the richness and diversity were higher in the P than in the O treatment, despite the similar rumen pH values of these treatments, suggesting that these changes are also based on other factors than ruminal pH. Other authors [7] in fact found that fermentation products, like VFA, have a large influence on the composition of the bacterial community in the rumen. Differences in the composition of the rumen microbiota between the O and P treatments could be explained by the different mechanisms of action of the supplements. However because of the complexity of the interactions among taxa it is difficult to clearly understand these biochemical processes in detail. The O treatment may have increased electron sinks for H2 allowing an increase of lactate utilization by bacteria which use the succinate-propionate pathway to synthesize succinate and (or) propionate [16, 19]. This hypothesis was confirmed by the lowest acetate to propionate ratio observed during the O treatment as reported by De Nardi et al. [23]. The blend supplemented during the P treatment contained flavonoids, e.g. polyphenols. Thanks to antimicrobial activity of theser products [24], the P treatment was reported to promote bacteria that metabolize lactate to VFA avoiding lactate accumulation in the rumen [17]. The higher richness and α-diversity found during the P treatment, which represent the diversity within sample [25], could be the consequence of a restraining effect towards some of the dominant taxa present in the rumen during high grain feeding, thereby leaving nutrient substrates available for other microbial groups. Despite differences in richness and diversity among treatments, the multivariate analysis based on PCoA explained a small percentage of the variance. Based on the microbial analysis results, the PCoA discriminant model was able to separate the control treatment from the P and O treatments. This discrimination demonstrated an association between the rumen pH drop resulting from high grain feeding and a decrease of Prevotella brevis as reported by Fernando [1], and an increase of abundance of the family Christensenellaceae. To the best of our knowledge this is the first contribute that reported a relation between rumen pH changes and the abundance of family Christensenellaceae. The change in the abundance of this family could be used as a microbial marker for a balanced ruminal pH condition. This relationship needs to be confirmed by further studies.

Results from Illumina sequencing showed that treatments affected the abundance of various phyla, orders, families, and genera. Among the species investigated with qPCR only Prevotella brevis was reduced by the P and O treatments. Matsui et al. and Khafipour et al. [5, 26] reported that Prevotella spp. are involved in starch degradation and utilization and thus grow better at low pH conditions.

Among the OTUs detected by Illumina sequencing, the family Pseudomonadaceae, order Actinomycetales, and Verrucomicrobia were reduced by the P treatment compared to the O treatment. These reductions were probably due to a potential antimicrobial activity of flavonoids as reported earlier [17, 24]. Decreasing the presence of these microbials was associated with the increase of several genera belonging to Bacteroidetes and Tenericutes. Similar trends were observed when comparing the P treatment with the control treatment, but in this case there were also increases in the abundances of some Firmicutes, such as genus Butyrivibrio. The high grain diet fed to heifers in this study resulted in a rumen microbiota typical of high concentrate feeding [6], although especially the P treatment affected this microbiota. In fact, the higher pH and the antimicrobial activity towards some microbial families, due to the P treatment, allowed an increase of bacteria which are usually favoured by low starch diets and a higher reticular pH. These include bacteria belonging to the family BS11 as reported by Zened et al. [22] and Paludibacter [27]. The latter was found to decrease when steers were transitioned from a hay-based to cereal grain-based diet. Even the genus Butyrivibrio showed the same response, i.e. increasing due to a higher rumen pH, since it is usually involved in fibre degradation [28].

The use of supplements over short times altered microbial populations. This shows that the rumen is a highly dynamic environment, and that the microbiota are capable of adapting to dietary changes. Polyphenols, in particular, were found to increase the richness and diversity of rumen microbiota in heifers feeding high starch diets.

Conclusions

Both additives (organic acid and polyphenols) changed microbial communities in the rumen, including the decline of Prevotella brevis and the increase of Christensenellaceae. These modifications were partially explained by the attenuation of the low reticulo-rumen pH caused by both additives. Polyphenols increased the microbial biodiversity in terms of richness and α-diversity more than organic acids, suggesting a benefit in using this additive in cattle fed high-grain diets. Furthermore, principal coordinate analysis demonstrated the discrimination between control and supplemented diets. The differences in the composition of the microbiota in the rumen between organic acid and polyphenol supplemented cows were likely due to an antimicrobial activity of the latter that widened the competition among different taxa.

Methods

Animals, experimental design and dietary treatments

The experimental protocol and all the procedures used in this study were approved by the Animal Ethics Committee of the University of Padova, Italy (CEASA, approval number 73/2012).

The experimental design and dietary treatments are described in De Nardi et al. [23]. Briefly, the study was carried out using six purebred Italian Holstein-Friesian non- pregnant heifers with an average body weight (BW) of 556 ± 33 kg (mean ± SD) according to a 3 × 3 Latin square design, with 3 experimental periods and 3 dietary treatments.

During each experimental period that lasted 22 days, a low starch (LS) diet was given ad libitum to the heifers for 14 d, followed by feeding a high starch (HS) diet for 8 d (from d 15 to d 22). Diets (Table 3) were provided as total mixed ration (TMR) once daily at 0800 h. From the d 18 to the d 22, barley meal was top dressed and its dose was gradually increased from 0.5 to 1.5 kg (250 g/day) to cause a rumen pH drop. Diet composition and the protocol to reduce rumen pH were similar to those used in an earlier study [29] and were adopted to prevent acute ruminal acidosis. During the HS diet feeding, the heifers were offered one of the following three dietary treatments: i) no supplement, C treatment; ii) a daily dose of 60 g of fumarate-malate mixture (RumenStabiliser®, DR. Eckel, Niederzissen, Germany), O treatment; iii) a daily dose of 100 g polyphenol-essential oil mixture (Anta®Phyt RU, DR. Eckel, Niederzissen, Germany), P treatment. The O supplement is a buffer and an organic acid mixture made of magnesium fumarate, sodium acetate, malic acid and sodium bicarbonate, whereas P is a mixture made of plant extracts, characterised by a high content of phenolic compounds comprising mostly flavonoids (1.88 mg/g). Feed samples were dried at 60 °C for 48 h, ground to 1-mm and then analysed for DM, crude protein (CP), crude fat and crude ash according to AOAC [30]. The NDF and ADF were analyzed according to Van Soest et al. [31], whilst the starch content was determined using high performance liquid chromatography [32].

Table 3

Ingredients and proximate composition of diets fed to heifers

	Dieta
Item	LS	HS
Ingredients, % of DM
Corn meal (0.5 mm)	23.0	34.0
Hay	24.0	22.0
Dehydrated alfalfa hay	18.0	16.0
Extruded soybean hull	7.0	2.0
Barley meal	7.0	9.0
Straw	6.5	5.5
Molasses	2.0	2.0
Sugar beet dry pulps	4.0	1.0
Corn gluten meal	6.0	2.0
Sunflower	-	2.0
Soybean meal	-	2.0
Vitamin and mineral mix	1.5	1.5
Extruded de-hulled soybean	1.0	1.0
Proximate composition
DM, %	88.9	88.8
Crude protein, % of DM	12.5	12.5
Crude fat, % of DM	3.5	3.3
NDF, % of DM	39.8	33.6
Crude ash, % of DM	7.8	7.4
ADF, % of DM	21.4	19.2
NFC,b % of DM	36.4	43.2
Starch, % of DM	24.0	30.0

aDiets: LS low starch, HS high starch

bNFC = 100 – (NDF + crude protein + crude fat + crude ash)

Sampling and analyses of ruminal fluid were also reported by De Nardi et al. [23]. Rumen fluid was transferred to 50-ml sterile tubes and immediately frozen in liquid nitrogen and then stored at −80 °C until analysis for quantitative real-time PCR (qPCR) and Illumina sequencing.

DNA extraction

Frozen rumen fluid samples were thawed at room temperature quickly and then kept on ice. One millilitre of rumen fluid was centrifuged at 15,000 × g followed by removing the supernatant. The DNA was extracted from pellets (200 mg of each sample) using ZR fecal DNA MiniPrepTM kit (Zymo Research Corp., Irvine, CA), which included a bead-beating step for the mechanical lysis of the microbial cells.

At the last step of the procedure, DNA was eluted from the column with elution buffer, and the concentration and purity of DNA were subsequently determined using a A260/280 spectrophotometer (NanoDrop 2000, Thermo Scientific, DE, USA). The DNA in all samples was diluted with the same elution buffer to a final nominal concentration of 20 ng/μl and quality was checked by PCR amplification of the 16S rRNA gene using universal primers 27 F (5’-GAAGAGTTTGATCATGGCTCAG-3′) and 342R (5′-CTGCTGCCTCCCGTAG-3′), as described by Khafipour et al. [5]. Amplicons were verified by agarose gel electrophoresis [5].

The DNA samples for qPCR analysis were diluted 2 ng/μl and aliquoted into 10 μL/vials, which is sufficient for testing one set of primers, in order to avoid repeated freeze-thaw cycles. All DNA samples were stored at −80 °C until analysis.

Sequencing

Library Construction and Illumina Sequencing.

Throughout the use of modified F515/R806 primers [33], the V4 region of 16S rRNA gene was adopted for PCR amplification. The reverse PCR primer was indexed with 12-base Golay barcodes and samples were analysed in the PCR reaction. ZR-96 DNA Clean-up Kit™ (ZYMO Research, CA, USA) was used to purify the PCR products. The V4 library was quantified using Picogreen dsDNA (Invitrogen, NY, USA). Multiple dilution steps with pre-chilled hybridization buffer (HT1) (Illumina, CA, USA) were used to obtain a 5 pM concentration of the pooled amplicons, which was determined using Qubit® 2.0 Fluorometer (Life technologies, ON, Canada). In order to reduce unbalanced and biased base compositions, 15 % of PhiX control library was spiked into the amplicon pool. Customized sequencing primers for read1 (5′-TATGGTAATTGTGTGCCAGCMGCCGCGGTAA-3′), read2 (5′-AGTCAGTCAGCCGGACTACHVGGGTWTCTAAT-3′) and index read (5′-ATTAGAWACCCBDGTAGTCCGGCTGACTGACT-3′) were synthesized and purified using polyacrylamide gel electrophoresis (Integrated DNA Technologies, IA, USA). A MiSeq platform (Illumina, CA, USA) was used to perform the 150 paired-end sequencing reaction.

Quantitative PCR analysis

Rumen fluid samples were analyzed by qPCR using primers for Eubacteria, Prevotella brevis, Prevotella bryantii, Fibrobacter succinogenes, Megasphaera elsdenii, Ruminococcus albus, Ruminococcus flavefaciens, Butyrivibrio fibrisolvens, Lactobacillus spp., and Streptococcus bovis as described in Wang et al. [34], Ozutsumi et al. [35], Denman and McSweeney [36], Khafipour et al. [5], Fernando et al. [1]. Real-time PCR was carried out using a CFX connect Real Time system (Bio-Rad Laboratories, Inc., USA). Each reaction mixture was run in triplicate in a volume of 15 μl in optical reaction plates (Thermo Fisher Scientific, U.K.) sealed with optical adhesive film (Applied Biosystems, Foster City, CA, USA). Amplification reactions were carried out with 7.5 μl Power SYBR green PCR master mix (Applied Biosystems, Foster City, CA, USA) mixed with the selected primer set at a concentration of 0.5 μM for each primer and 10 ng of genomic DNA. Amplification consisted of one cycle of 95 °C (10 min), 40 cycles of denaturation at 95 °C (15 s), and annealing/extension at 60 °C (1 min). Final melting analysis was conducted by slow heating from 65 °C to 95 °C in order to assess the specificity of the primer set. Data were normalized for Eubacteria using the universal bacteria 16S RNA gene primer sets, which detect all bacterial strains.

Bioinformatic analysis

The PANDAseq assembler [37] was used to merge overlapping paired-end Illumina fastq files. All the sequences with mismatches or ambiguous calls in the overlapping region were removed. The QIIME software [38] was used to analyze the output fastq file and the reads obtained were filtered and selected to discard the unsuitable ones. The sequences were assigned to Operational Taxonomic Units (OTUs) after being processed by UCHIME [39] and UCLUST [40]. The RDP classifier [41] was used to find correspondence between each OTU and known taxa. Chao1, abundance based coverage estimated (ACE), Shannon and Simpson indices were calculated using the Phyloseq package (version 1.9.2) in R (version 3.0.2).

Statistical analyses

After verifying the normality of residuals (PROC UNIVARIATE) of the observed species, richness, diversity and qPCR microbial variables were analyzed using the MIXED procedure with a compound symmetry structure using a linear model that included the fixed effects dietary treatment, period and their interaction, and the random effect of animal. If significant treatment effects were detected (P < 0.10), the LSmeans were compared using the probability of differences (PDIFF) option and the Tukey adjustment test. In order to obtain a normal distribution and homogeneous residual error, Simpson (diversity) and qPCR microbial variables were log transformed (SAS, release 9.3, 2010). Using the same linear model, the β-diversity indices (Bray-Curtis and unweighted Unifrac) were submitted to distance-based permutational multivariate analyses of variance (PERMANOVA) according to Anderson [42].

Principal coordinate analysis (PCoA) was conducted to evaluate differences in community structure among dietary treatments (β-diversity). PCoA was generated using unweighted Unifrac distance [43] (R package vegan version 2.0–10, 2013).

At the bacterial taxa level, pair-wise differential expression analyses were performed in the nbinomWaldTest method of DESeq2 [44] according to a statistical model consisting in the effects of dietary treatment, period and heifer; these were considered significant at P < 0.05 and trends were discussed at 0.05 < P < 0.10 (DESeq2, R package version 1.5.34, 2012).