Metagenome sequencing to analyze the impacts of thiamine supplementation on ruminal fungi in dairy cows fed high-concentrate diets

Animal care and procedures were in accordance with the Chinese guidelines for animal welfare and approved by Animal Care and Use Committee of the Chinese Academy of Agricultural Sciences. Twelve Chinese Holstein dairy cows (627 ± 19.9 kg BW; 180 ± 8 DIM) in second-parity fitted with 10 cm ruminal cannulas (Bar Diamond, Parma, ID) were divided into three treatments. Treatments included a control diet (CON; 20% starch, dry matter basis), a high-concentrate diet (HC; 33.2% starch, dry matter basis) and high-concentrate diet supplemented with 180 mg thiamine/kg dry matter intake (HCT). The diets were formulated according to NRC (2001) to meet or exceed the energy requirements of Holstein dairy cows yielding 20 kg of milk/d with 3.5% milk fat and 3.0% true protein. Details of ingredient analysis and chemical composition of diets were shown in Additional file: Table S1. Cows were fed twice equally at 06:00 and 18:00 h each day, and thiamine (thiamine hydrochloride, purity ≥ 99%; Wanrong Science and Technology Development Co., Ltd., Wuhan, China) was add via the rumen cannula twice daily after diets were offered for a 21 day period. Throughout the experimental periods, the cows were housed in individual stalls and had free access to fresh water. 1

During the experimental period, automatic feeding equipments (made by Institute of Animal Science Chinese Academy of Agricultural Sciences, Beijing, China and NanShang Husbandry Science and Technology Ltd. Henan, China) were used to record dry matter intake. Milking facilities of Afimilk were applied to record milk production of each cow. On day 21, rumen contents were sampled from cranial, caudal, dorsal, and ventral sites of rumen at 3 h after the morning feeding. Collected samples were strained through four layers of cheesecloth to obtain rumen fluid. Rumen fluid was then divided into two parts. One part was processed to analyze the pH value, thiamine content and rumen volatile fatty acid (VFA) content. The other part was put into the liquid nitrogen immediately after adding stabilizer and then stored at − 80 °C for DNA extraction. Ruminal pH value was measured using a portable type pH meter (Testo 205, Testo AG, Lenzkirch, Germany) immediately after rumen fluid sample collection. Thiamine concentration in rumen fluid was detected by high performance liquid chromatography (HPLC) according to Analytical Methods Committee (2000). Individual and total VFA (TVFA) of ruminal fluid were determined by gas chromatograph (GC-2010, Shimadzu, Kyoto, Japan).

DNA for metagenomics was extracted from the rumen fluid samples by using the QIAamp DNA Stool Mini Kit (Qiagen, Hilden, Germany) according to manufacturer’s protocols. The DNA concentration and purity were quantified with TBS-380 and NanoDrop2000, respectively. DNA quality was examined with the 1% agarose gels electrophoresis system.

DNA was fragmented to an average size of 300 bp using Covaris M220 (Gene Company Limited, China) for paired-end library construction. Paired-end library was prepared by using TruSeqTM DNA Sample Prep Kit (Illumina, San Diego, CA, USA). Adapters containing the full complement of sequencing primer hybridization sites were ligated to the Blunt-end fragments. Paired-end sequencing was performed on Illumina HiSeq 4000 platform (Illumina Inc., San Diego, CA, USA) using HiSeq 3000/4000 PE Cluster Kit and HiSeq 3000/4000 SBS Kits according to the manufacturer’s instructions (http://www.illumina.com↗).

The quality control methods used the following criteria: 3′ and 5′ ends were stripped using SeqPrep (https://github.com/jstjohn/SeqPrep↗) (Aronesty 2013). Low-quality reads (length < 50 bp or with a quality value < 20 or having N bases) were removed by Sickle (https://github.com/najoshi/sickle↗). Sequences that lost their mated reads were considered as single reads and were used in the assembly procedure. Reads were aligned to the cow genome by NCBI (https://www.ncbi.nlm.nih.gov↗) and any hit associated with the reads and their mated reads were removed. Filtered reads were considered for the next step of the analysis.

Considering that K-mers, varying from 1/3 to 2/3 of reads length were used in assembly by SOAP denovo (http://soap.genomics.org.cn↗, Version 1.06), which is based on De Brujin graph construction (Li et al. 2010). Scaffolds with a length over 500 bp were retained for statistical tests; we evaluated the quality and quantity of scaffolds generated by each assembly and finally chose the best K-mer which yielded the minimum scaffold number and the maximum value of N50 and N90. Then, scaffolds with a length over 500 bp were extracted and broken into contigs without gaps. Contigs were used for further gene prediction and annotation.

Open reading frames (ORFs) from each sample were predicted using MetaGene (http://metagene.cb.k.u-tokyo.ac.jp/↗) (Bao et al. 2014). The predicted ORFs with length being or over 100 bp were retrieved and translated to amino acid sequences using the NCBI translation table (http://www.ncbi.nlm.nih.gov/Taxonomy/taxonomyhome.html/index.cgi?chapter=tgencodes#SG1↗). All sequences from gene sets with a 95% sequence identity (90% coverage) were clustered as the non-redundant gene catalog by the CD-HIT (http://www.bioinformatics.org/cd-hit/↗). Reads after quality control were mapped to the representative genes with 95% identity using SOAP aligner (http://soap.genomics.org.cn/↗) and gene abundance in each sample was evaluated (Qin et al. 2014). DIAMOND (Huson et al. 2016), a new alignment tool for aligning short DNA sequencing reads to a protein reference database was employed for taxonomic annotations by aligning non-redundant gene catalogs against NCBI NR database with e-value cutoff of 1e⁻⁵ and Score > 60. Based on NCBI Microbial Taxonomy Information Database, species annotation information of genes was obtained and relative abundance of species was calculated at the level of Kingdom, Phylum, Class, Order, Family, Genus and Species.

Ruminal pH and thiamine content were analyzed using one-way ANOVA of SAS 9.2. P value < 0.05 was considered to be significance and 0.05 ≤ P < 0.10 was considered as a tendency. Spearman correlations between fungi communities and ruminal fermentation variables or thiamine content were assessed using the PROC CORR procedure of SAS 9.2. Relative abundance of all phyla, and the top 30 species of fungi were chosen to conduct the correlation analysis. A correlation matrix was created and visualized in a heatmap format using R package version 3.3.1. The abundance of fungi communities and ruminal variables were considered to be correlated with each other when the absolute value of correlation coefficients (r) were above 0.55 and P-value below 0.05. Barplot, principal coordinate analysis (PCoA), hierarchical clustering analysis (HCA) and heat map for different rumen fungi species were conducted using R package version 3.3.1.

All the raw sequences were submitted to the NCBI Sequence Read Archive (SRA; http://www.ncbi.nlm.nih.gov/Traces/sra/↗), under Accession Number SRP144478.

As shown in Table 1, HC feeding significantly decreased dry matter intake, milk production, ruminal pH value and the concentrations of thiamine and acetate in rumen fluid, while significantly increased the concentration of propionate and isovalerate compared with CON diet (P < 0.05). However, these changes caused by HC diet were inversed by thiamine supplementation (P < 0.05). Butyrate, valerate, isobutyrate and total VFA concentrations were not affected by dietary treatments (P > 0.05).

In the present study, 12 metagenomic libraries were constructed. Based on the quality control methods, all pollution data were removed and approximately 45,000,000 reads of per sample were acquired. The GC content of each sample was between 40–60% except HCT4, which had a lower GC content of 36.36%. The assembly results showed that a minimum of 130,000 contigs per sample were obtained. The N50 value was between 760 and 890 bp for the samples (shown in Table 2). The aligned reads were annotated for gene prediction. The results indicated that numbers of prediction genes were mostly between 200,000 and 400,000. However sample HC3 was detected to contain a significant higher prediction gene numbers while sample HCT4 contained a significant lower prediction gene numbers.

Principal coordinates analysis (PCoA) based on unweighted UniFrac distance metrics were conducted to compare the three treatments. As shown in Fig. 1, PCoA axes 1 and 2 accounted for 94.24% and 1.19% of the total variation, respectively. The cows offered HC separated from those offered CON and HCT by PCo1 effectively. Similarly, cows fed CON separated from those fed with HCT by PCo2.

In the present study, approximately 1050 species of rumen fungi were identified across all samples in seven characterized phyla and one uncharacterized phylum. All taxonomy results are shown in Additional file 2: Excel. As shown in Fig. 2, Ascomycota, Basidiomycota and Chytridiomycota were the three most phyla for all three treatments. Relative abundance of Ascomycota was significantly higher in HC treatment than the other two treatments. Neocallimastigomycota, characterized as a phylum of strictly anaerobic fungi, was identified including 3 genera (g__Neocallimastix, g__Orpinomyces, and g__Piromyces) and 11 species (s__Neocallimastix_frontalis, s__Neocallimastix_patriciarum, s__Orpinomyces_joyonii, s__Orpinomyces_sp._C1A, s__Orpinomyces_sp._OUS1, s__Orpinomyces_sp._PC-2, s__Orpinomyces_sp._ukk1, s__Piromyces_equi, s__Piromyces_sp, and s__Piromyces_sp._E2) (shown in Additional file 3). All identified fungi and their relative abundances in the level of phyla are displayed in Table 3. HC feeding significantly decreased the abundance of all phyla of fungi (P < 0.05) when compared with CON treatment. In contrast, thiamine supplementation significantly increased the abundance of all phyla of fungi (P < 0.05) compared with HC treatment.

Hierarchical clustering analysis (HCA) and heat map analysis was conducted to further understand the effects of HC feeding and thiamine supplementation on ruminal fungi profile. Because little information of fungi phyla was detected, taxonomy level of classes and species were used for HCA and heat map analysis. As shown in Fig. 3, samples of HCT treatment gathered into a cluster which was significantly separated from the other two treatments no matter at the level of class or species. Samples of HC treatment were separated from those of CON at the level of class but not at the level of species. Across all three treatments, classes of Ascomycota, Basidiomycota, Chytridiomycota, and Blastocladiomycota gathered into a big cluster. At the species level, species of Mucorales gathered into a big cluster and Saccharomyces_cerevisiae which is often used in dairy feed, gathered into a single cluster.

As shown in Fig. 4, the relative abundance of all phyla was negatively correlated with the concentration of TVFA, acetate, propionate and isovalerate (r < − 0.55, P < 0.05). In particular, Chytridiomycota was positively correlated with butyrate concentration (r > 0.55, P < 0.05). There was a trend that the abundances of all phyla of fungi were positively correlated with ruminal pH and thiamine content (r < 0.55, 0.10 > P > 0.05). Similarly, all the top 20 species of fungi had significantly negative correlated to TVFA, acetate, propionate and isovalerate (r < − 0.55, P < 0.05). Whilst, Tetrapisispora_blattae and Allomyces_macrogynus both had positively relationships with butyrate content (r > 0.55, P < 0.05). In addition, a trend was also found that the abundances of all the top 20 species of fungi were positively correlated with ruminal pH and thiamine content (r < 0.55, 0.10 > P > 0.05).

In the present study 7 phyla and almost 1050 species were identified, the sequencing depth and data size were improved when compared with previous studies (Kumar et al. 2015; Mercado-Garcia 2017; O’Malley et al. 2014), and many species of the rumen fungi were first reported such as Blastocladiomycota, Chytridiomycota, Microsporidia and Glomeromycota result. The proportion of Ascomycota in the whole 7 phyla of fungi in HC treatment significantly increased compared with that in the other two treatments, which indicated that phylum Ascomycota had better resistibility to the changed condition caused by HC diet. As the most representative species of Ascomycota, Saccharomyces cerevisiae was always used as the additive to alleviate subacute ruminal acidosis in lactating dairy cows (Alzahal et al. 2014), and the result further elaborated the function of Saccharomyces cerevisiae in alleviating SARA (Alzahal et al. 2014). Neocallimastigomycota, although took up a small proportion, attracted more attention in the rumen microbiome (Hibbett et al. 2007) because of its strictly anaerobic characterization and the ability of producing wide array of cell-bound and cell-free cellulolytic, hemicellulolytic, glycolytic, and proteolytic enzymes (Liggenstoffer et al. 2010). Sequencing and functional researches of Neocallimastigomycota were proceeded in recent years in order to explore the mechanism of its functions in rumen such as the high-throughput expressed sequence tag (EST) analysis of Neocallimastix frontalis (Mi et al. 2009) and the analysis of the Orpinomyces genome (Youssef et al. 2013). The results of the present study that the proportion of Neocallimastigomycota was increased after thiamine supplementation in the whole rumen micro-biomass may provide useful information in investigating the mechanism of Neocallimastigomycota.

In the present study, HC significantly decreased the abundance of fungi (P < 0.05) compared with CON treatment. This finding was in line with several previous studies in which a greater fungal diversity and higher fungal abundance were detected in animals fed high-fiber diets than those offered high-grain diets (Belanche et al. 2012; Kumar et al. 2015). The main function of rumen fungi is to degrade diet fiber and product absorbable compounds for the host animal (Li 2003; Theodorou et al. 1996). Because of the higher gain and lower fiber content in HC, less substrates were provided for fungi growth, therefore leading to the decreased abundance of rumen fungi in HC treatment.

However, the fungal abundance increased significantly after thiamine supplementation in HC treatment which indicated that thiamine supplementation promoted the proliferation and metabolism of rumen fungi. Thiamine is an essential nutrient for dairy cows and plays important roles in carbohydrate metabolism (Kamke et al. 2010; Subramanya et al. 2010). The promotion of carbohydrate metabolism and its interaction with ruminal microbiome might have contribute to the increased fungal abundance.

In ruminal anaerobic environment, ruminal fungi energy generation organelles are represented by hydrogenosomes. These organelles, under anoxic conditions, decarboxylate malate into acetate, CO2 and H2 with concomitant production of energy in the form of ATP (Van 2010; Yarlett and Hackstein 2005). Rumen fungi exhibit tremendous diversity of CAZymes, which work together to degrade carbohydrates and cellulose into monosaccharides and then the monosaccharides are further hydrolyzed into pyruvate (Hall and Mertens 2017). As the cofactor of pyruvate formate-lyase (PFL) and pyruvate:ferredoxin oxidoreductase (PFO), thiamine supplementation promoted the reaction of pyruvate generating into acetyl-CoA (Xue et al. 2018), enhanced the utilization of pyruvate and promoted the degradation efficiency of carbohydrates which lead to more ATP formed to support the metabolism and proliferation of ruminal fungi in the current study. Besides, thiamine supplementation significantly decreased rumen VFAs content and increased ruminal pH (Pan et al. 2016), which promoted the activity of CAZymes (Fliegerova et al. 2015). Therefore, fungal abundance significantly increased in the current study.

Interaction with other microorganisms may also contribute to the change of ruminal fungi. Ruminal microorganism works in harmony to keep the health and productivity of ruminants (Hobson and Stewart 1997; Mickdam et al. 2016). In particular, many positive associations were observed between fungi and bacteria, such as Caecomyces1—Ruminococcus1 and Piromyces3—Succiniclasticum (Smidt et al. 2010; Tapio et al. 2017). Thiamine was detected to induce a significant increase of fungal growth in fungi and bacteria co-culture system (Shimizu et al. 2016) which indicated that thiamine promoted the mutualism of fungi and bacteria. In addition, previous study of our team indicated that HC significantly decreased the abundance of Bacteroides and Ruminococcus 1 (Pan et al. 2017). As stated above, the decreased Ruminococcus 1 may partly contribute to the decreased rumen fungi in HC treatment. Thiamine supplementation stimulated the growth of the lactate utilizing bacteria Ruminococcus 1 (Pan et al. 2017), and thus increased the abundance of rumen fungi indirectly.

In conclusion, HC significantly decreased while thiamine supplementation significantly increased the abundance of rumen fungi. These results indicated that thiamine supplementation may effectively attenuate rumen metabolic disorder caused by HC diet through buffering the ruminal pH, shifting the rumen fermentation pattern and increasing the abundance of ruminal fungi. The findings in this study could therefore contribute to the further understanding of the mechanism of thiamine’s function in dairy cows and provide effective strategies to improve dairy cows’ health under high-concentrate feeding regime.

We thank the Beijing Key Laboratory for Dairy Cow Nutrition, Beijing University of Agriculture, Beijing, China for supporting the experimental equipment. We sincerely thank Dr. Yiguang Zhao for the contribution to English editing. This study was funded by the Project of National Nature Science Foundation of China (Grant No. 31572435) and the National Key Research and Development Plan (2016YFD0700205, 2016YFD0700201).

All authors declare that they have no competing interests.

All the raw sequences were submitted to the NCBI Sequence Read Archive (SRA; http://www.ncbi.nlm.nih.gov/Traces/sra/↗), under Accession Number SRP144478.

All authors carefully read and agree to be accountable for all aspects of the work.

This article does not contain any studies with human participants or animals performed by any of the authors.

The study was funded by the Project of National Nature Science Foundation of China (Grant No. 31572435) and the National Key Research and Development Plan (2016YFD0700205, 2016YFD0700201).

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

All the raw sequences were submitted to the NCBI Sequence Read Archive (SRA; http://www.ncbi.nlm.nih.gov/Traces/sra/↗), under Accession Number SRP144478.