《Microbiology Spectrum》：Metagenomic and metabolomic analyses of rumen fiber digestion in Mongolian cattle fed fresh grass versus hay
ABSTRACT
蒙古牛因其瘤胃微生物组具有强大的纤维降解能力而表现出卓越的粗饲料耐受性，使其能够在蒙古高原季节性波动的条件下有效利用低质量饲草。本研究比较了饲喂青草（FG）与干草的牛只的瘤胃微生物组成、碳水化合物活性酶（CAZymes）谱、发酵参数和代谢途径，以阐明纤维消化背后的微生物-代谢物相互作用。将30头非妊娠雌性蒙古牛（460 ± 35 kg，3-4岁）随机分为两组（每组n = 15）：一组放牧采食FG，另一组舍饲并饲喂秋季收获的干草（HG）。每组随机选取6头牛进行瘤胃液采集和多组学分析（每组n = 6，总计n = 12）。与FG组相比，HG组在瘤胃发酵参数中显示出乙酸摩尔比例和乙酸/丙酸比值的增加，同时丙酸和丁酸的摩尔比例降低。宏基因组分析显示HG组中拟杆菌目（Bacteroidales）细菌和厌氧真菌（包括Neocallimastixsp. JGI-2020a和Piromycessp. E2）的丰度更高。功能注释进一步表明HG组中碳水化合物代谢途径富集，同时CAZymes多样性更高，尤其是那些参与半纤维素和果胶降解的酶。代谢组学鉴定出13种差异丰度的碳水化合物代谢物，其中葡萄糖酸内酯在HG组中上调。此外，代谢组中鉴定出的碳水化合物代谢途径证实了宏基因组功能注释的可靠性。相关性网络分析揭示了拟杆菌科（Bacteroidaceae）细菌、Neocallimastixsp. JGI-2020a和Piromycessp. E2与乙酸、降解半纤维素的糖苷水解酶（GHs）以及碳水化合物代谢途径呈正相关。总之，干草饲喂通过增加拟杆菌目和厌氧真菌、多样化CAZymes（特别是半纤维素酶/果胶酶）以及上调碳水化合物代谢，增强了蒙古牛的瘤胃纤维降解，反映了微生物对低质量饲草的适应性。
INTRODUCTION
锡林郭勒草原是中国北方典型的温带大陆性草原，对区域畜牧业生产和生态系统稳定至关重要。该生态系统由严酷的气候塑造，年均气温约3°C，年降水量低（约290毫米），寒冷季节长，导致饲草可用性和质量存在明显的季节性波动。在生长季节（6月至9月），牧草表现出高粗蛋白（CP）含量（10-20%），足以满足牲畜营养需求。相反，寒冷季节（冬春季）CP含量急剧下降（5-8%），中性洗涤纤维（NDF）含量增加（60-70%），显著降低了饲草质量并限制了动物性能。因此，及时在秋季收获和加工天然牧草对于最大化其营养价值并确保牲畜在缺料的冬春季节生存至关重要。随着饲草营养价值随植物成熟度动态变化，新鲜牧草提供高营养效益，适合即时放牧需求，而干草则提供稳定的饲料以支持牲畜在饲草稀缺时期。因此，了解牛和绵羊如何利用这两种饲草形式对于优化该草原生态系统的畜牧业生产至关重要。
蒙古牛是锡林郭勒草原的本地品种，对当地牧民不可或缺，提供奶、肉和役力，同时是区域生计和经济稳定的基石。在传统放牧实践下，天然牧草是这些牛的唯一饲料来源。为了应对饲草可用性和营养成分的季节性波动，蒙古牛进化出了具有高效纤维降解能力的独特瘤胃微生物群落，使其能够在极端寒冷、漫长旱季和营养缺乏的情况下持续生长和生产。瘤胃通常被称为天然的“发酵罐”，蕴藏着包含细菌、古菌、真菌和原生动物在内的多样化微生物群落，它们协同将植物纤维素降解为可溶性糖。这些糖随后被发酵成挥发性脂肪酸（VFAs）、微生物蛋白和其他对宿主能量和生长至关重要的营养物质。瘤胃微生物群落不直接降解纤维素，而是产生多种多样的CAZymes来驱动纤维降解。Hagen等人和He等人的研究表明，细菌主要参与半纤维素降解，而厌氧真菌专门负责木质纤维素降解，其中GH48和GH5家族酶最为普遍。此外，瘤胃微生物组成和酶活性受反刍动物物种、放牧季节、日粮成分和饲料添加剂的影响。对其他反刍动物（如牦牛和青年母牛）的研究表明，富含纤维的日粮促进了纤维降解细菌的富集，特别是在拟杆菌门（Bacteroidota）和厚壁菌门（Firmicutes）中，从而增强了纤维降解和VFA生产。
尽管认识到瘤胃微生物在纤维消化中的作用，但对蒙古牛纤维降解的具体机制研究仍然有限——特别是核心细菌和真菌类群、关键CAZymes以及将顽固植物生物质转化为VFAs的代谢途径。高通量多组学技术（如宏基因组学和代谢组学）的出现，现在能够系统地探索这些复杂的微生物-宿主相互作用。本研究假设：（i）干草饲喂的蒙古牛表现出纤维降解细菌丰度增加，增强了用于纤维素降解的CAZymes分泌；（ii）它们激活了额外的碳水化合物代谢途径，增加了VFA合成以满足能量需求。为了验证这些假设，我们比较了饲喂青草（FG）与干草的蒙古牛的瘤胃微生物组成、CAZymes谱、发酵参数和代谢途径，以阐明驱动瘤胃纤维消化的微生物-代谢物相互作用。
MATERIALS AND METHODS
动物管理与样品采集
本研究在内蒙古自治区锡林郭勒盟乌拉盖的一个蒙古牛育种保护基地进行。2023年7月，从锡林郭勒草原的一个放牧牛群中随机选择30头雌性蒙古牛（460 ± 35 kg，3-4岁，非妊娠）。将牛随机分配到两组（每组n = 15）。一组（FG）自然放牧以针茅属（Stipaspp.）、榆属（Ulmusspp.）和羊茅属（Festucaspp.）为主的FG，另一组（HG）饲养在同一围栏内并饲喂前一年秋季收获的干草。未向任何一组提供补充饲料，所有动物均自由采食和饮水。试验包括15天的适应期和随后45天的试验期。
在试验期的第43-45天，使用四点取样法每天从两组收集饲草样品。将3天的样品彻底混合并合并以测定日粮营养成分（表1）。在第45天，从每组随机选择6头牛（每组n = 6，总计n = 12）。在早晨饲喂前（约晚上饲喂后12小时）使用胃管采样器（GCYQ-1-A, 武汉科利博设备有限公司）从每头动物收集约200 mL瘤胃液。选择这个时间点是为了与昼夜周期的低谷对齐，尽量减少餐后发酵高峰带来的变异性。最初的100 mL瘤胃液被丢弃以减少唾液污染。剩余的液体通过四层无菌纱布过滤，分装到冷冻管中，在液氮中速冻，并储存在-80°C直至分析。
饲草营养分析与瘤胃发酵谱
将200克饲草样品在65°C下干燥至恒重以测定干物质（DM）含量。将干燥样品使用切割磨（ZM200, Retsch GmbH, 北京, 中国）研磨过1 mm筛用于营养分析。使用凯氏定氮法测定CP含量。使用AOAC方法测定NDF和酸性洗涤纤维（ADF）含量。
对于VFA分析，将瘤胃液样品在4°C解冻，向1 mL瘤胃液中加入0.25 mL 25%（重量/体积）的偏磷酸。混合物在4°C下以20,000 × g离心15分钟。上清液使用气相色谱仪（TP-2060T, 北京北分天普仪器技术有限公司, 北京, 中国）分析VFA浓度，如Xu等人所述。分析条件如下：2 m × 6 mm玻璃柱，填充涂有10% FFAP的Chromosorb W AW-DMCS；高纯N2作为载气，流速30 mL/min；温度程序：100°C（保持2分钟），以8°C/分钟升至180°C（保持5分钟）；进样口温度220°C；进样体积1 μL；FID检测器温度250°C。燃料气体为H2，流速30 mL/min，氧化气体为空气，流速300 mL/min。使用外标法和峰面积积分进行定量。
瘤胃微生物的宏基因组谱分析
使用E.Z.N.A. Soil DNA Kit（Omega Bio-tek, Norcross, GA, USA）按照制造商方案从瘤胃液样品中提取总基因组DNA。分别使用TBS-380荧光计和NanoDrop 2000分光光度计评估DNA浓度和纯度。通过1%琼脂糖凝胶电泳评估DNA完整性。高质量DNA使用Covaris M220超声破碎仪（Gene Company Limited, China）剪切至约400 bp用于双末端文库制备。使用NEXTFLEX Rapid DNA-Seq Kit（Bioo Scientific, Austin, TX, USA）按照制造商说明构建测序文库。双末端测序（2 × 150 bp）在Illumina NovaSeq 6000平台（Illumina Inc., San Diego, CA, USA）上使用NovaSeq 6000 S4 Reagent Kits在美格生物技术有限公司（上海, 中国）进行。
对原始测序读数进行质量过滤和修剪，以去除低质量碱基和接头。使用BWA（v0.7.9a）将过滤后的读数与牛（Bos taurus）参考基因组比对，以去除宿主来源的序列。使用MEGAHIT（v1.1.2）将高质量非宿主读数组装成重叠群。使用MetaGene预测长度≥300 bp的重叠群的开放阅读框。使用CD-HIT（v4.8.1）在95%同一性下构建非冗余基因目录。使用SOAPaligner（v2.21）以≥95%同一性将clean reads映射到目录，以生成基因丰度谱（以TPM表示）。使用DIAMOND（v0.8.35）针对NCBI NR数据库（访问日期2020年6月4日）对代表性序列进行物种注释。针对KEGG数据库（v94.2）进行功能注释，并使用hmmscan（HMMER）针对CAZy数据库进行CAZymes注释。使用自定义R脚本分析物种水平微生物类群的功能贡献。每个类群对功能途径的相对贡献计算为其基因丰度除以与该途径相关的所有类群的总丰度。将值标准化为每个样本总和为1。
瘤胃液的非靶向代谢组学分析
收集瘤胃液样品（100 mL），将100 μL等分试样转移到1.5 mL离心管中。每个等分试样与400 μL含有0.02 mg/mL L-2-氯苯丙氨酸（作为内标）的提取溶液（乙腈-甲醇，1:1，体积/体积）混合。样品涡旋30秒，随后在5°C和40 kHz下超声处理30分钟。通过在-20°C孵育样品30分钟诱导蛋白质沉淀。然后在4°C下以13,000 × g离心样品15分钟。收集所得上清液并在氮气流下蒸发至干。将干燥的样品重新溶解在100 μL乙腈-水（1:1，体积/体积）中，在5°C（40 kHz）下超声处理5分钟，并在4°C下以13,000 × g离心10分钟。将最终上清液转移到液相色谱-质谱（LC-MS）进样瓶中进行分析。在Vanquish UHPLC系统 coupled with a Q Exactive HF Orbitrap质谱仪（均来自Thermo Fisher Scientific, Waltham, MA, USA）上进行代谢组学分析。12个样品符合非靶向代谢组学分析的质量标准，并被纳入下游分析。
使用Progenesis QI（Waters Corporation, Milford, MA, USA）处理原始LC-MS数据，并通过美格生物云平台进行分析。保留在任何组中80%样品中检测到的代谢特征。对于低于定量限的特征，使用最小值插补值，所有数据通过总和进行标准化。在质量控制样品中相对标准偏差>30%的变量被排除。剩余数据经过log10转换以生成用于后续统计分析的最后数据矩阵。
对处理后的矩阵进行多变量分析。使用ropls R包进行主成分分析（PCA）和正交偏最小二乘判别分析（OPLS-DA）。通过七折交叉验证评估模型性能，并使用200次置换检验进一步验证。所有置换的Q2值均为负或接近零（R2 = 0.806, Q2 = -0.144），确认不存在过拟合。使用OPLS-DA模型中的变量重要性投影（VIP > 1）和Student's t-检验（P< 0.05）的组合鉴定差异代谢物。代谢组学t-检验未应用多重检验校正，因为VIP过滤作为特征选择步骤。
使用KEGG数据库将差异代谢物映射到生化途径。进行途径分类和富集分析以确定差异代谢物的生物学功能。使用Python的“scipy.stats”模块进行富集分析，识别与实验处理显著相关的途径。
统计分析
统计分析旨在评估两种日粮处理对瘤胃VFAs和瘤胃微生物组的影响。对于满足正态性和方差齐性假设的瘤胃VFA数据，使用SPSS（版本27.0, IBM Corp., Armonk, NY, USA）进行独立样本t-检验以比较组间均值差异。
使用美格生物云平台分析宏基因组微生物组和功能谱数据。基于Bray-Curtis相异度进行主坐标分析（PCoA）和PCA，以可视化样本间群落结构的差异。然后应用ADONIS（999次置换）检验组间整体群落组成差异的统计显著性。在个体特征水平，使用Wilcoxon秩和检验鉴定分类学和功能特征的差异丰度。采用线性判别分析效应大小（LEfSe）检测显著富集的KEGG途径（LDA得分 > 2.0, P< 0.05）。基础的Kruskal-Wallis检验经过FDR校正，仅报告P< 0.05的特征为显著。
RESULTS
瘤胃发酵参数
如表2所示，FG组中丙酸和丁酸的浓度及摩尔比例显著高于HG组（P< 0.05）。相反，HG组中乙酸的摩尔比例和乙酸/丙酸比值显著高于FG组（P< 0.05）。相比之下，总VFA浓度或乙酸浓度未观察到显著差异（P> 0.05）。
瘤胃微生物群落的变化
瘤胃宏基因组测序产生了1,069,072,464条原始读数，平均每个样本89,089,372 ± 5,799,900条读数（平均值 ± 标准差）。经过质量过滤和宿主基因组去除后，保留了1,050,475,550条高质量读数，平均每个样本87,539,629 ± 5,897,312条读数。重叠群组装产生了19,817,220个重叠群，平均每个样本1,651,435 ± 238,969个重叠群。
蒙古牛的瘤胃宏基因组主要由细菌（91.03%）组成，其次是古菌（5.72%）、真核生物（2.21%）、病毒（1.03%）和未分类微生物（0.02%）。Wilcoxon秩和检验显示，HG组中细菌和真核生物的相对丰度显著高于FG组，而古菌丰度显著较低（P< 0.05）。为了阐明瘤胃纤维降解的机制，我们分析了细菌和真菌群落的组成和组间差异。基于Bray-Curtis距离的β-多样性分析并通过PCA可视化显示，HG组和FG组之间细菌和真菌群落存在明显的聚类和显著分离，表明不同饲喂制度下的群落结构存在差异。在门水平，拟杆菌门（Bacteroidota）（57.91%）和芽孢杆菌门（Bacillota）（34.49）占主导地位，其中Prevotellasp.（21.08%）和拟杆菌目细菌（Bacteroidales bacterium）（19.87%）为最丰富的物种。优势真菌门为壶菌门（Chytridiomycota）（68.58%）、子囊菌门（Ascomycota）（17.10%）和毛霉门（Mucoromycota）（13.87%），其中Neocallimastixsp. JGI-2020a（22.92%）、Piromycessp. E2（13.37%）、Neocallimastix californiae（12.24%）和Anaeromyces robustus（12.09%）为主要物种。在前10种细菌物种中，HG组中Bacteroidales bacterium和Bacilli bacterium的丰度显著高于FG组，而FG组中Muribaculaceaebacterium的丰度显著更高（P< 0.05）。在前10种真菌物种中，Neocallimastixsp. JGI-2020a、Piromycessp. E2、Neocallimastix californiae、Anaeromyces robustus和Piromyces finnis在HG组中丰度显著更高，而FG组中Rhizopus arrhizus和Fusarium avenaceum的丰度显著更高（P< 0.05）。
瘤胃微生物群落的CAZyme谱
基于Bray-Curtis距离并通过PCoA可视化的β-多样性分析显示，FG组和HG组之间的CAZyme谱存在显著差异。在所有鉴定出的CAZyme家族中，糖苷水解酶（GH）家族在两组中最为丰富，并且在HG组中的丰度显著高于FG组（P< 0.05）。功能贡献分析表明，Prevotellasp.和拟杆菌目细菌（Bacteroidales bacterium）是GHs的主要细菌贡献者，而Neocallimastixsp. JGI-2020a、Piromycessp. E2和Neocallimastix californiae是主要的真菌贡献者。比较分析鉴定出37个在组间丰度存在显著差异的GHs（P< 0.05）。特别是，HG组表现出与半纤维素和果胶降解相关的GH家族丰度升高。具体而言，GH43、GH10、GH130、GH8、GH54和GH113与半纤维素降解相关，而GH105和GH67与果胶降解相关。
瘤胃微生物群落的功能特征
基于Bray-Curtis距离并通过PCA可视化的β-多样性分析显示，FG组和HG组之间瘤胃微生物功能谱存在显著差异。根据KEGG一级功能分类，代谢是主要类别，突出了其在瘤胃微生物活动中的核心作用。为了进一步阐明代谢差异，将LEfSe应用于代谢类别内的KEGG三级功能，鉴定出38个显著富集的途径（LDA得分 > 2.0, P< 0.05）。这些富集的途径主要与碳水化合物和氨基酸代谢相关。鉴于本研究关注瘤胃纤维降解机制，特别强调了碳水化合物代谢途径。七个KEGG三级途径在组间显示出显著差异，均在HG组中富集（P< 0.05）。这些途径包括氨基糖和核苷酸糖代谢、糖酵解/糖异生、半乳糖代谢、磷酸戊糖途径、戊糖和葡萄糖醛酸相互转化、果糖和甘露糖代谢以及乙醛酸和二羧酸代谢。
瘤胃代谢组学分析
为了进一步阐明瘤胃代谢物的特征及其与纤维降解细菌的关联，对12个瘤胃液样品进行了非靶向LC-MS代谢组学分析。PCA显示HG组和FG组之间瘤胃代谢物谱存在明显分离，表明它们的代谢组存在显著差异。OPLS-DA进一步证实了这种分离，并具有稳健的预测性能（R2Y = 0.748, Q2Y = 0.997）。火山图分析鉴定出组间612种差异代谢物，其中与FG相比，HG中有403种上调，209种下调。这612种差异代谢物使用KEGG数据库进行注释，映射到36条KEGG途径。包含超过10种化合物的途径包括与癌症、消化系统、氨基酸代谢、脂质代谢、其他氨基酸代谢、辅因子和维生素代谢、碳水化合物代谢和膜转运相关的途径。为了研究瘤胃纤维降解的机制，我们重点关注了与碳水化合物代谢途径相关的13种代谢物。六种代谢物在HG组中富集：石胆酸3-O-葡萄糖醛酸、3α,6β,7α,12α-四羟基-5β-胆烷酸、蔗糖、海藻糖、葡萄糖酸内酯和棕榈酰葡萄糖醛酸。相反，七种代谢物在FG组中富集：α-d-葡萄糖、富马酸、d-葡萄糖醛酸内酯、2-氧代丁酸、琥珀酸、γ-氨基丁酸（GABA）和l-谷氨酸。
对碳水化合物代谢中13种差异代谢物的KEGG富集分析鉴定出14条富集途径：氨基糖和核苷酸糖代谢、磷酸戊糖途径、糖酵解/糖异生、柠檬酸循环（TCA循环）、淀粉和蔗糖代谢、丙酮酸代谢、C5-支链二元酸代谢、丙酸代谢、戊糖和葡萄糖醛酸相互转化、半乳糖代谢、乙醛酸和二羧酸代谢、果糖和甘露糖代谢、抗坏血酸和醛糖酸代谢以及丁酸代谢。其中，七条途径在代谢组学和宏基因组学分析中重叠：氨基糖和核苷酸糖代谢、糖酵解/糖异生、半乳糖代谢、磷酸戊糖途径、戊糖和葡萄糖醛酸相互转化、果糖和甘露糖代谢以及乙醛酸和二羧酸代谢。
组学与表型关联的整合分析
进行Spearman等级相关分析以评估前10种差异微生物物种、GH家族、微生物功能途径和VFAs之间的关系。仅保留|r| > 0.7且P< 0.05的相关性用于解释。拟杆菌科细菌（Bacteroidaceae bacterium）与乙酸以及与
