Multi-omics analyses reveal that the gut microbiome and its metabolites promote milk fat synthesis in Zhongdian yak cows

Housing and feeding systems

Zhongdian yak cows from Shangri-La in the Yunnan province (26−34°N, 94−102°E; 3600 m) were selected and gut samples were collected using rectal palpation as the main research objects. To minimize environmental variation, before sampling commenced, all yak cows foraged freely in the same pasture for at least 14 days.

Collection of milk and gut samples

All samples from Zhongdian yak cows were collected according to the Animal Welfare Committee of Shandong Agricultural University (permit number SDAUA-2018-022). A total of 36, 8-year-old Zhongdian yak cows from Shangri-La, Yunnan province were selected and milked once daily during October 2021 (mean monthly temperature: 15 °C; altitude: 3,600 m). Milk samples were collected from each cow at milking, and analyzed for fat using mid-infrared instruments (FOSS, Denmark). The five Zhongdian yak cows with the highest milk fat percentage (>7%, 8.70 ± 1.89%, H group) and the five with the lowest milk fat percentage (<5%, 4.12 ± 0.43%, L group) were selected for gut analysis. Gut samples were obtained by rectal palpation and stored in liquid nitrogen until analyzed.

DNA extraction, metagenomic sequencing, and metagenomics data processing

DNA was extracted using the PowerSoil^® DNA Isolation kit (MO BIO Laboratories, Inc., Carlsbad, CA). The quality and quantity of DNA were measured using Qubit V.2.0 fluorometer quantitation (Thermo Fisher Scientific, Inc., Waltham, MA, USA). The genomic DNA was randomly sheared into fragments and the library was constructed using TrueLib DNA Library Rapid Prep Kit for Illumina (ExCell Bio, Shanghai, China). Metagenomic library sequencing was performed on an Illumina Novaseq PE150 (150 bp paired-end sequencing, 500 pb inserts) at Biomarker Technologies Co. Ltd. (Beijing, China).

The trimmomatic (version 0.33) was used to control the quality of each dataset by trimming the 3′-end and the 5′-end of the reads. After cutting low-quality bases (quality scores <20) and removing short reads (<50 bp), the reads were aligned to the bovine genome (bosTau8 3.7, DOI: https://doi.org/10.18129/B9.bioc.BSgenome.Btaurus.UCSC.bosTau8) using bowtie2 (version 2.2.4). MEGAHIT (version 2.2.4) (Li et al., 2015) and QUAST (Gurevich et al., 2013) were used to de novo assemble and assess the filtered reads for each sample. MetaGeneMark (version 3.26) (Zhu, Lomsadze & Borodovsky, 2010) was used to predict the coding region from the assembled contigs with length >500 bp. MMseqs (version 4.6.6) (Steinegger & Soding, 2017) was used to pool the assembled contigs and construct non-redundancies based on identical contigs (>95% identity). DIAMOND (version 0.9.24) was used to estimate the abundances after mapping the original sequences to the predicted genes.

Taxonomic and functional annotation from gut metagenomes

The taxonomic assessment of gut microbiota was performed using DIAMOND (version 0.9.24) against the non-redundant protein sequence database (Nr database) (Yu & Zhang, 2013). Taxonomic profiles were conducted at the phylum, class, family, genus, and species levels, with relative abundances calculated at each level. Microbial taxa with a relative abundance >0.1% in at least 50% of the Zhongdian yak cows in each group were used for downstream analysis. Contigs were annotated using DIAMOND against the EggNOG (Evolutionary Genealogy of Genes: Non-supervised Orthologous Groups, v4.0) (Powell et al., 2014) and KEGG databases (Kyoto Encyclopedia of Genes and Genomes, 2017-03) (Kanehisa et al., 2021) with an E-value <1e−5. The CAZymes annotation was performed using hmmer (version 3.0) against the CAZy database (carbohydrate-active enzymes database, version 6.0) (Helbert et al., 2019). In a downstream analysis, the abundances of EggNOG categories, KEGG Orthology (KO) pathways, and CAZymes were normalized into counts per million reads. The EggNOG categories, KEGG pathways, and CAZymes with cpm >5 in at least 50% of the animals in each group were used for the downstream analysis. The microbial network was constructed using the top 80 species, and associations with R > 0.5 and P < 0.05 were used for the downstream analysis.

Analysis of metabolites in gut samples

All Zhongdian yak cow gut samples were prepared according to the manufacturer’s instructions of the Waters Acquity I-Class PLUS ultra-high performance liquid tandem Waters Xevo G2-XS QT high resolution mass spectrometer (Macioszek et al., 2021). After freezing in liquid nitrogen, the frozen gut samples were extracted with one mL of extracting solution, and then the samples were sonicated with KQ-500DE ultrasonic cleaner at 0 °C for 15 min, and centrifuged at 4 °C and 12,000 g for 15 min. The supernatants were collected into new tubes, and the extracts were dried in a vacuum concentrator, and then filtered through 0.22 µm nylon membrane filters.

The liquid chromatography was equipped with a Waters Acquity UPLC HSS T3 column (1.8 um 2.1*100 mm). Solvent A was composed of water and 0.1% formic acid, solvent B was composed of acetonitrile and 0.1% formic acid, and the mobile phase was composed of solvent A and solvent B. The gradient profile is outlined as follows: 0–2 min, 98% A; 2–11 min, 2–98% B; 11–13 min, 98% B; 13–15 min, 2–98% A. The flow rate was 200 µL/min, and the injection volume was 10 µL. The sample temperature was maintained at 10 °C, and the column temperature was maintained at 35 °C.

The Waters Xevo G2-XS QTOF high resolution mass spectrometer was used to collect the primary and secondary mass spectrometry data in MSe in the accompanying software (MassLynx V4.2, Waters) and to test metabolites eluted from the column in both positive and negative ion modes. In each data acquisition cycle, dual-channel data acquisition was performed on low collision energy (2V) and high collision energy (10−40 V), with a 0.2 s scan time. The capillary voltages were as follows: 2,000 V in positive ion mode and −1,500 V in negative ion mode. The sampling cone voltage was 30 V. The ion source temperature was set at 150 °C, the desolvent gas temperature was 500 °C, the backflush gas flow rate was 50 L/h, and the desolventizing gas flow rate was 800 L/h.

MassLynx V4.2 was used to collect the raw data. The peak extraction, peak alignment, and other data processing operations were then processed by the Progenesis QI software based on the METLIN database and Biomark’s self-built library for identification. The theoretical fragment identification and mass deviation were all within 100 ppm.

Statistical significance testing

The milk fat percentage of Zhongdian yak cows was statistically analyzed using the R-package (R v3.02; P < 0.05). Gut microbial domains, phyla, genera, and species were compared using a Wilcoxon rank-sum test, with an FDR adjusted P value <0.05 considered significantly different (Lin et al., 2021). The abundances of microbial EggNOG categories, KEGG Orthology pathways, and CAZymes were compared between the two groups using a linear discriminate analysis effect size (LEfSe), with an LDA score >2 and P value <0.05 indicating significant differences. A T test was performed between the two groups, with an FDR adjusted P value <0.05 and VIP (Variable Importance in the Projection) >1 considered significantly different metabolites.

Milk fat percentage and the gut metagenome of the selected yak cows

Five Zhongdian yak cows with high milk fat percentages (8.7 ± 1.89%) and five Zhongdian Yak cows with low milk fat percentages (4.12 ± 0.43%) were selected for gut metagenome analyses (Table 1). After quality control and the removal of host genes, a total of 828,374,098 reads were generated using metagenome sequencing, with 82,837,409.80 ± 6,903,448.80 per sample (10.6084/m9.figshare.19467455; DOI: 10.6084/m9.figshare.19434545). A total of 3,1451,262 contigs were generated (the N50 length of 622 ± 46 bp) after de novo assembly, with 3,145,126 ± 223,620 per sample. The gut metagenome of Zhongdian yak cows consisted of 76.81% bacteria, 1.20% archaea, 0.43% viruses, and 0.15% fungi, with 21.36% of the gut metagenome unassigned and 0.51% unclassified (Fig. 1A). The downstream comparison of the gut microbial taxa between the high milk fat (H) and low milk fat (L) Zhongdian yak cows focused on bacteria.

The gut bacterial composition and differences between the H and L Zhongdian yak cows

Moving forward, we analyzed the composition and differences in the gut bacteria of H and L Zhongdian yak cows (Figs. 1B–1F). In the analysis of the abundance of bacteria, at the phyla level, the dominant bacteria were Firmicutes (47.40 ± 1.59%), Bacteroidetas (16.82 ± 2.27%), and Proteobacteria (1.79 ± 0.19%); at genus level, the dominant bacteria were Bacteroides (5.12 ± 0.75%), Clostridium (3.56 ± 0.22%), and Alistipes (2.45 ± 0.40%); at species level, the dominant bacteria were Firmicutes_bacterium_CAG:110 (3.85 ± 0.47%), Clostridiales_bacterium (2.40 ± 0.27%), and Ruminococcaceae_bacterium (1.78 ± 0.32%). We then compared the differences in the gut microbial domains between H and L Zhongdian yak cows using a Wilcoxon rank-sum test, and found significant differences at the genus and species levels (P < 0.05). At the genus level, 91 genera, including Halolamina, Aphanocapsa, and Limnobacter, were significantly more abundant in the gut of H Zhongdian yak cows, while 53 genera, including Epicoccum, Haloarcobacter, and Palleronia, showed significant enrichment in the gut of L Zhongdian yak cows. At the species level, 42 species were significantly more abundant in the gut of H Zhongdian yak cows, with the three most significantly elevated being Ruminococcus_sp._CAG:724, Ruminococcus_sp._CAG:382, and Acetobacter_sp._CAG:267. There were 47 species that were significantly less abundant in the gut microbiome of the H group compared to the L Zhongdian yak cows, including Clostridia_bacterium_UC5.1-1D1, Anaerorhabdus_furcosa, and Ruminococcus_sp._TF12-2.

Functional analysis of the gut microbiome and differential functions between H and L Zhongdian yak cows

To analyze the different functions of the gut microbes between the H and L Zhongdian yak cows, we identified the functions of the gut microbes using EggNOG, KEGG profiles, and CAZy. The EggNOG results identified three main categories including “Metabolism” (22.51 ± 1.83%), “Information storage and processing” (14.87 ± 3.00%), and “Cellular processes and signaling” (10.96 ± 1.12%). A total of 25 EggNOG classes were obtained (Fig. 2A), including “Translation, ribosomal structure and biogenesis” (6.42 ± 0.32%), “Carbohydrate transport and metabolism” (5.62 ± 0.22%), “Energy production and conversion” (4.16 ± 0.21%), and “Lipid transport and metabolism” (1.29 ± 0.05%). There were no differences in the EggNOG classes identified between the H and L Zhongdian yak cows.

The KEGG analysis identified 175 endogenous third-level pathways belonging to four first-level categories (Fig. 2B) including “Metabolism” (77.14%), “Genetic information processing” (12%), “Cellular processes” (6.86%), and “Environment information processing” (4%). At the second level, 22 categories were obtained, with “Xenobiotics biodegradation and metabolism” (9.71%), “Lipid metabolism” (8.57%), “Carbohydrate metabolism” (8.57%), “Amino acid metabolism” (8.00%), and “Metabolism of terpenoids and polyketides” (8.00%) being the most abundant. We then compared the identified KEGG pathways (Fig. 2C), and found five pathways that were significantly enriched in the H gut microbiome, including “Phosphonate and phosphinate metabolism” (P = 0.006, df = 8), “Carbapenem biosynthesis” (P = 0.023, df = 8), “Histidine metabolism” (P = 0.014, df = 8), “Biosynthesis of unsaturated fatty acids” (P = 0.029, df = 8), and “Biosynthesis of ansamycins” (P = 0.008, df = 8); while only one pathway, “Bacterial chemotaxis” (P = 0.028, df = 8), was significantly enriched in the gut of L Zhongdian yak cows.

The CAZyme analysis identified 125 glycoside hydrolases (GHs), 85 glycosyltransferases (GTs), 79 carbohydrate-binding modules (CBMs), 16 carbohydrate esterases (CEs), 23 polysaccharide lyases (PLs), and 11 auxiliary activities (Poon et al., 2021) (Fig. 2D). The five most dominant CAZymes were GT2 (0.21 ± 0.01%), GH13 (0.17 ± 0.01%), GH109 (0.17 ± 0.01%), GT4 (0.12 ±0.00%), and CE1 (0.11 ± 0.00%). Among all the CAZymes identified, 22 were significantly enriched in the gut of H Zhongdian yak cows (P < 0.05), including GH123 (P = 0.007, df = 8), PL12 (P = 0.031, df = 8), and CE11 (P = 0.012, df = 8), and 20 CAZymes were significantly enriched in the gut of L Zhongdian yak cows (P < 0.05), including PL1 (P = 0.028, df = 8), PL9 (P = 0.016, df = 8), and CBM77 (P = 0.019, df = 8). There were significantly more CAZymes identified in H Zhongdian yak cows than in L yak cows (Fig. 2E).

Differences analysis of association networks of the gut microbiome in H and L Zhongdian yak cows

Microbial network for both H and L yak cows were constructed at the species level based on the abundance and correlation of gut microbes (Fig. 3A). Both positive and negative associations were identified in both networks. Notably, in both H and L Zhongdian yak cows, the gut microbial network contained a greater number of positive associations than negative associations. Interestingly, the number of negative associations in H Zhongdian yak cows was significantly less than the number of positive associations in L Zhongdian Yak cows, and the majority of associations found within both networks were unique to that network. We hypothesized that this dramatic difference may reflect changes in the ecosystem associated with milk fat synthesis. We then analyzed the taxonomic relationships between the species within both networks. The microbial networks in the gut of both H and L Zhongdian yak cows exhibited a large abundance of taxonomically distant networks, and many of the bacterial associations identified were unique to that network. Notably, the negative associations found in H Zhongdian yak cows occurred between species that were taxonomically similar.

Our analysis of microbial networks found that Ruminococcaceae_bacterium and Clostridiales_bacterium were positively associated with other species in H Zhongdian yak cows, and negatively correlated with other species in L Zhongdian yak cows; Firmicutes_bacterium exhibited a negative association with other species in H Zhongdian yak cows and a positive association with other species in L Zhongdian yak cows. In H Zhongdian yak cows, there were unique positive correlations between Anaerotruncus_sp_CAG:390 and Clostridium_sp_CAG:413, and between Sarcina_sp_DSM_11001 and other species, and a unique negative correlation between Firmicutes_bacterium_CAG:110_56_8 and Firmicutes_bacterium_CAG:631. All species that exhibited unique correlations with others in H Zhongdian yak cows belonged to the phylum Firmicutes. This may indicate that these are influential bacterial species within the networks and that they may play important roles in milk fat synthesis in Zhongdian yak cows.

Gut metabolome analysis

A total of 2,400 compounds were identified in the gut metabolome. Using the T-test and Variable Importance in the Projection (VIP) to filter the relative concentrations of gut metabolites, 57 metabolites were significantly higher and 49 metabolites were significantly lower in the gut of H Zhongdian yak cows (P < 0.05, VIP >1, Fig. 3B). The KEGG annotation based on these 106 significantly different gut metabolites revealed the enrichment of nine first-level categories in H Zhongdian yak cows (Fig. 3C), including “Lipid metabolism” (23.61%), “Amino acid metabolism” (18.06%), “Global and overview maps” (12.5%), “Nervous system” (11.11%), “Metabolism of cofactors and vitamins” (11.11%), “Digestive system” (9.72%), “Signal transduction” (6.94%), “Membrane transport” (4.17%), and “Cancer: overview” (2.78%). Of these nine pathways, only “Metabolism of cofactors and vitamins” and “Signal transduction” were also enriched in L Zhongdian yak cows. To further identify the potential gut microbiome-metabolome interactions, a Spearman’s rank correlation and constrained correspondence analysis between the gut metabolome and microbiome were performed, with the results revealing 27 significant correlations (—CC—>0.80, CCP <0.05). Among these 27 correlations, 12 positive correlations were identified between: myristic acid and six genera, choline and four genera, and 5-aminopentanal and two genera (Fig. 3D). Firmicutes and Proteobacteria were the dominant bacterial phyla in the gut microbiome of H Zhongdian yak cows. These results indicate that the regulation of Firmicutes (Desulfocucumis, Anaerotignum, Dolosiccus)-myristic acid and Proteobacteria (Catenovulum, Comamonas, Rubrivivax, Marivita, Succinimouas)-choline interactions may be the main contributors to milk fat synthesis in Zhongdian yak cows.