Dynamics of rumen bacterial composition of yak (Bos grunniens) in response to dietary supplements during the cold season

Ethics statement and site description

All procedures involved in animal care and their use were in strict agreement with the guidelines approved by the Animal Administration and Ethics Committee of Lanzhou Institute of Husbandry and Pharmaceutical Sciences (LIHPS) of CAAS, China (SYXK-2018-0011). After collection, the samples were processed strictly in agreement with the rules of LIHPS. The yaks were delivered for experiments by the yak owner and after the completion of the experiment, all yaks keep on living for other experimental work.

The feeding trial was conducted for 90 days from December 2018 to February 2019 at Hezuo Yak Breeding Cooperatives (35°08′N, 102°99′E, 2,960 m above sea level) of the Tibetan Autonomous Prefecture of Gannan, Gansu province, China. The pastures in the area was characterized by alpine meadows and inhabited by vegetation species of Kobresia graminifolia, Kobresia humilis, Elymus nutan, Kobresia pygmaea, Anaphalis lacteal, Polyginum viviparum, Potentilla fruticose, Cortaderia jubata, and Sibiraea angustata. The average temperature and precipitation during trial were −3.9 °C and 4.8 mm, respectively.

Animal management, diets, and experimental design

Twenty-one 3-year old female Gannan yaks having an initial average body weight of 178 ± 2.71 kg from the same herd were randomly placed into three dietary treatment groups (n = 7). The animals were lodged individually in tie-stalls and had 9 m² space with bedding in the night for regular activities. The experiment started after 15 days of adaptation period to familiarize the animals with diet, facilities, and staff. All animals remained together in the same flock and three groups grazed freely in the natural alpine pastures from 7:00 to 18:00 with free access to water. To check pasture availability, the height of dominant species from pasture was randomly measured with a ruler. About 10 samples were collected by hand as representative of forage selected by yak, weighted, and dried to calculate dry matter.

The first group served as a control group and grazed the natural pasture without any supplement (GG). The grazing yaks in the second group were supplemented with 1 kg/d oats hay (OG), while grazing yaks in the third group were fed with 0.9 kg/d concentrate pellets (CG), which was based on the amount of feed voluntarily consumed by yaks during adaptation period to ensure the complete absence of orts. The concentrate supplement consisted of corn 500 g/kg, wheat bran 220 g/kg, rapeseed cake 240 g/kg, urea 16 g/kg, salt 17 g/kg, additives 1.5 g/kg, calcium hydrogen phosphate 5.5 g/kg, calcium 5.6 g/kg and phosphorous 2.6 g/kg on DM basis. The nutritional composition of pasture herbage, oats hay, and concentrate pellets is listed in Table 1. Before this study, all yaks only grazed natural pasture and did not receive any supplement and during the experimental period, about 2–3 yaks grazed per hectare of pasture. The yaks were offered oats hay and concentrate pellets twice a day at 07:00 and 18:00.

Forage DMI was measured with lignin as an internal marker and chromic oxide as an external marker as described previously (Sampaio et al., 2011). For this purpose, forage and fecal samples were collected (Ding et al., 2010). Each yak was given 20 g of chromic oxide in a gelatin capsule orally once daily for 9-days preliminary dosing period and 10-days fecal sampling period. The fecal samples were collected daily during the period of administration of chromic oxide. The samples were dried at 60 °C for 48 h, ground to pass through 1 mm sieve, and concentrations of internal and external markers were measured (Myer et al., 2015; Van Soest, Robertson & Lewis, 1991). Forage DMI was estimated from the fecal output of the internal marker adjusted for the supplement contribution. Finally, the total dry matter intake of supplemented groups was calculated by adding forage and supplement intake (Kartchner, 1980).

Sample collection

The mixed forage samples were taken from the pasture grazed by yaks. The quadrat of 0.5 × 0.5 m² size was randomly placed across the pasture to mimic the yak selection of forage to graze. The forage samples were clipped in triplicates and stored for later nutritional composition analysis.

Rumen sample (liquid contents, 100 mL/yak) was taken at the end of the feeding trial via an oral stomach tube in the morning before feeding to explore rumen fermentation and bacterial composition. The tube was properly washed with fresh water every time before taking a new sample and the first 15–20 ml of the sample was always discarded to prevent saliva contamination. The rumen contents were immediately stored in liquid nitrogen and then preserved at −80 °C for future analysis.

Nutrient analysis

Pasture forage and feed samples were dried at 65 ^∘C for 72 h in an oven, pulverized, and passed through a 1 mm sieve. Dry matter (DM) was measured by drying the samples at 135 °C for 3 h (AOAC # 930.15). Crude protein (CP, 984.13) was measured by AOAC methods (Association of Official Analytical, 1990). The neutral detergent fiber (NDF) and acid detergent fiber (ADF) were measured according to the method reported by Van Soest (Van Soest, Robertson & Lewis, 1991). The sodium sulfite was used in the NDF method and the values were corrected for ash content.

Rumen fermentation parameters analysis

Rumen pH was measured immediately after sampling with a portable pH meter (Kadiya, 6010, Shenzhen, China). The ammonia nitrogen (NH₃-N) concentration was determined by a colorimetric method (Weatherburn, 1967). The VFA concentration was determined as previously described in (Ahmad et al., 2020). Briefly, the frozen rumen liquid sample was thawed, mixed by vortexing, and five mL of rumen liquid was centrifuged at 3000 × g for 10 min. In a centrifuge tube, 1 mL of supernatant was taken and 0.2 mL of a metaphosphoric acid solution was added to it. The solution was mixed and held for 30 min in an ice bucket and then centrifuged at 1,0000× g at 4 °C. The supernatant was taken into a new centrifuge tube and preserved for testing at 4 °C. The gas chromatography (7890A GC system Agilent Technologies Inc, Santa Clara, CA, USA) with FID detector was used to measure VFA concentration. The gas chromatographic conditions and subsequent test procedures were conducted as described previously (Isac et al., 1994).

DNA extraction and PCR amplification

The metagenomic DNA of 21 rumen fluid samples was extracted by the cetyltrimethylammonium bromide method. DNA was dissolved in the 200 µL of elution buffer and then preserved at −20 °C. DNA integrity and concentration were determined by 1.5% agarose gel electrophoresis and NanoPhotometer^® spectrophotometer (Implen, Westlake Village, CA, USA), respectively. For rumen bacterial analysis, universal primer pairs (515F-806R) with barcodes were used to PCR amplify the V4 region of the 16S rRNA gene (Li et al., 2016). PCR amplification was done by using Phusion^® High-Fidelity PCR Master Mix with GC Buffer from New England BioLabs. The amplified products were detected by 2% agarose gel electrophoresis and purified by QIAquick PCR Purification Kit to build libraries. The libraries were built according to the manufacturer’s protocol using TruSeq^® DNA PCR-Free sample preparation kit. The libraries were then quantified by using Qubit dsDNA High Sensitivity Assay kit by Invitrogen. The Illumina HiSeq2500 PE250 system was used for paired-end sequencing using the standard protocol.

Sequencing and data processing

Raw sequences were prepared and QIIME (Quantitative Insights Into Microbial Ecology, Version 1.7.0) software was used for their analysis (Caporaso et al., 2010). After sequencing, primer sequences, barcodes, and low-quality sequences were truncated (Li et al., 2016). Chimeric sequences were identified and discarded from the dataset using UCHIME Algorithm with reference to Gold database to finally get effective tags (Edgar et al., 2011). UCLUST software version 7.1 (http://drive5.com/uparsel/) was used to cluster effective tags into OTUs with 97% identity (Edgar, 2013) and bacterial taxa were identified from Greengenes database using representative sequences (Wang et al., 2007).

Alpha diversity indices such as Chao1, ACE, Simpson, Shannon, and Goods-coverage were calculated using QIIME software. STAMP (Statistical analysis of taxonomic and functional profiles, Version 2.1.3) was used to construct a principal component analysis (PCA) plot to illustrate the significant difference between samples at the OTU level. A correlation heat map was generated in GraphPad Prism version 8.00 for Windows (http://www.graphpad.com/).

Statistical analysis

Intake data, rumen fermentation parameters and relative abundances of bacterial taxa (Phylum and genus levels) were analyzed using the mixed procedure in SPSS software (Version 20.0) (IBM, Armonk, NY, United States). The dietary group was considered as a fixed effect and the individual yak was taken as a random effect. Residuals were checked for normality using Shapiro–Wilk test. When the residuals did not follow a normal distribution the data was log-transformed. Bonferroni correction method was used for multiple comparisons of groups (Haynes, 2013). Analysis of similarities (Anosim) was performed in R studio by using anosim function of the vegan package to examine the grouping variations of each group separately in the PCA plot. Spearman’s correlation coefficient of rumen fermentation parameters with the 10 more abundant ruminal bacterial genera and genera significantly affected by the treatments were calculated using GraphPad Prism 8.0.2 and a heat map was also generated. Significance was declared at P <0.05 and P values were modified using a false discovery rate to exclude false-positive results.

Data availability

The sequencing data for this study was stored in the European Nucleotide Archive (ENA) at EMBL-EBI under accession number PRJEB36248.

Rumen fermentation parameters

Total dry matter intake (TDMI) of the CG group was significantly higher compared to GG and OG groups (Table 2). The results of rumen fermentation parameters are displayed in Table 2. In comparison to GG and OG groups, the CG group showed significantly (P < 0.01) higher concentrations of NH₃-N and TVFA in rumen fluid. OG group also showed higher concentrations of NH3-N and TVFA in rumen fluid than GG group. GG group showed a higher concentration of acetate (P = 0.03) than CG groups, while both groups displayed no significant differences with the OG group. GG and OG groups showed significantly lower concentrations of propionate (P = 0.01), butyrate (P = 0.01), isobutyrate (P < 0.01), and isovalerate (P < 0.01), and AP ratio (P < 0.01) than CG group. However, no significant difference was observed between GG and OG groups. Rumen pH was significantly (P < 0.001) lower in the CG group followed by OG and GG groups.

Bacterial diversity analysis

After filtering, quality control, and chimera removal, a total of 1,097,941 high-quality bacterial sequences were produced and the Good’s coverage was more than 95.2% with an average sequence read length of 253 bp. The rarefaction curve reached the plateau, indicating that the OTUs number did not increase with the increasing size of data and the sequencing data volume was sufficient. The clean reads were clustered into 6,340 OTUs on a 97% sequence similarity basis. The analysis of OTU recognized 3,378 common OTUs in three groups, while 651, 177, and 870 unique OTUs were identified in the GG, OG, and CG groups, respectively (Fig. 1).

Simpson index did not show a significant difference among groups (P = 0.21, Table 3). Shannon, chao1, and ACE indices were significantly lower in the OG group compared to GG and CG groups, while no significant variations were observed between GG and CG groups.

The STAMP software was used to further analyze the difference in the bacterial community among groups by the PCA plot (Fig. 2). PCA plot was not very informative in studying the differences among groups as all the samples assembled close to each other. However, Anosim analysis displayed a significant difference in the composition of rumen bacteria at OTU level between GG and CG (R = 0.44, P < 0.001) groups, between GG and OG (R = 0.21, P = 0.005) groups and between CG and OG (R = 0.18, P = 0.038) groups.

Rumen bacterial composition

The taxonomic analysis showed that the yak rumen bacterial community comprised of 23 phyla and 487 genera. At the phylum level, Bacteroidetes (42.0% ± 9.1%), Proteobacteria (14.6% ± 11.0%) and Firmicutes (15.9% ± 3.34%) were found to be the most dominant phyla in three groups (Fig. S1). Verrucomicrobia (6.9% ± 2.7%), Fibrobacteres (4.8% ± 2.7%), and Tenericutes (4.2% ± 5.2%) were other abundant phyla present in three groups. The low abundance phyla among three groups included Cyanobacteria, Actinobacteria, Acidobacteria, and Lentisphaerae.

Changes at the phylum level in the three groups are shown in Table 4. The abundance of Tenericutes in the OG group was significantly lower (P = 0.02) compared to CG and GG groups, while no significant variations were recorded between GG and CG groups. The phylum Cyanobacteria was significantly higher (P = 0.04) in the CG group compared to the OG group, while the GG group showed no significant difference with OG and CG groups.

To further explore the microbial abundance, classification was performed at the genus level. Prevotella (16.2% ± 6.4%) was the dominant genus followed by Fibrobacter (4.84% ± 2.7%), Ruminobacter (2.18% ± 2.4%) and CF231 (1.6% ± 0.03%), accounting for 24.8% of all genera (Fig. S2). Other dominant genera included RFN20 (1.25% ± 0.04%), Anaeroplasma (1.0% ± 0.03%), BF311 (1.6% ± 0.96%), Halomonas (1.6% ± 0.83%), Bacteroides (0.29% ± 0.05%) and Lactobacillus (0.22% ± 0.07%).

The variations in the relative abundance of the top 20 genera in three groups were examined (Table 5). The abundance of Ruminobacter was significantly (P < 0.01) lower in GG compared to OG and CG groups. However, no significant difference was displayed between OG and CG groups. The relative abundance of Oscillospira was significantly (P = 0.014) higher in the GG group compared to the CG group, while no significance was recorded between GG and OG and between OG and CG groups. The relative abundance of Ruminococcus (P = 0.04) was significantly higher in the GG group compared to the OG group, while no significant difference was observed between GG and CG groups and between OG and CG groups. The abundance of Clostridium (P = 0.007) was significantly higher in the GG group compared to OG and CG groups. However, no significant difference was observed between OG and CG groups. The genera Selenomonas (P = 0.014), and Corynebacterium (P = 0.043) were significantly higher in the CG group compared to GG and OG groups, while no difference was observed between GG and OG groups.

Correlation analysis between bacterial genera and rumen fermentation parameters

Spearman’s correlation was implemented to identify the association of rumen fermentation parameters with the abundance of top 10 and significantly varied bacterial genera (Fig. 3). Genus Ruminococcus was positively correlated with acetate (r = 0.44, P = 0.04) ruminal proportions and negatively associated with ruminal TVFA (r = −0.43, P = 0.016) concentrations and isobutyrate (r = −0.71, P < 0.001 0.02), and isovalerate (r = −0.58, P < 0.01) proportions. Genus Clostridium was negatively correlated with ruminal TVFA (r = −0.44, P = 0.04) and NH₃-N (r = −0.44, P = 0.04) concentrations, and isobutyrate (r = −0.56, P < 0.01) proportions. Finally, genus RFN20 was negatively correlated with ruminal isobutyrate (r = −0.48, P = 0.02) proportions. The rest of the genera did not show any significant relationship with any parameters.