牦牛是位于青藏高原及其邻近地区的大型反刍动物。目前，我国牦牛约1 600万头，占世界牦牛总数的95%[1-2]。青海省牦牛数量占全国牦牛总数的40%。牦牛可以提供牛奶、肉以及可制成织物的毛发[3-4]。反刍动物具有独特的消化特性，其肠道菌群在物种生物学中的地位比较突出[5]。目前，反刍动物肠道微生物已被广泛研究，如羚牛[6]、绵羊[7]、林麝[8]等。影响动物肠道微生物的组成和多样性变化的因素很多[9-10]，如季节变化会引起林麝的肠道菌群中厚壁菌门和拟杆菌门的比例波动，且冬春季节的肠道菌群多样性高于夏秋季节[11]。本研究利用16S rRNA测序技术初步测定和分析冷暖季节放牧牦牛粪便菌群的组成、多样性及功能预测的变化，为放牧条件下牦牛补饲调节肠道菌群和后续研究牦牛肠道微生物变化与饲养管理提供参考。1材料与方法1.1试验设计在青海省玉树州歇武镇牧业村的同一牧场中，选取体重相近的1.5岁和3.5岁的雄性牦牛各5头。牦牛自由放牧，采食牧草。该区域位于东经96°53'、北纬34°56'，平均海拔4 500 m，年平均温度-1.6 ℃。冷暖季节不同年龄放牧牦牛体重见表1。10.13557/j.cnki.issn1002-2813.2021.15.017.T001表1冷暖季节不同年龄放牧牦牛体重年龄暖季冷季1.5岁119.40±5.82117.70±9.283.5岁189.40±13.69176.00±16.43kg分别于2019年9月30日（暖季）和2019年12月30日收集粪样（冷季）。暖季采集的1.5岁牦牛粪样为BOC组，冷季采集的1.5岁牦牛粪样为BTC组，暖季采集的3.5岁牦牛粪样为DOC组，冷季采集的3.5岁牦牛粪样为DTC组。观察牦牛采食牧草的情况，分别于2019年9月30日（暖季）和2019年12月30日（冷季）采集天然牧草样品，在青海省畜牧兽医科学院青海省牦牛工程技术研究中心实验室进行常规营养成分测定。牧草营养水平见表2。10.13557/j.cnki.issn1002-2813.2021.15.017.T002表2牧草营养水平（风干基础）项目暖季冷季灰分（Ash）5.325.45粗蛋白质（CP）11.282.81粗脂肪（EE）3.662.58中性洗涤纤维（NDF）42.7048.60酸性洗涤纤维（ADF）31.7340.64钙（Ca）0.310.32磷（P）0.170.07%1.2样品采集禁食约8 h，牦牛固定于限制槽中，戴上无菌橡胶手套，从牦牛直肠中取粪样。将每头牦牛的粪样收集到两个5 mL的无菌冻存管中，液氮下保存，用于提取DNA。1.3试验方法1.3.1DNA提取与PCR扩增将冷冻样品在冰上解冻，在涡流混合器中均匀混合。使用CTAB/SDS法提取样品总基因组DNA。在1%琼脂糖凝胶上监测DNA浓度和纯度。根据浓度，使用无菌水稀释DNA至1 mg/L。利用通用引物V338F（5'-ACTCCTACGGG AGGCAGCAG-3'）和V806R（5'-GGACTACHVGGG TWTCTAAT-3'）对微生物基因组总DNA的16S rRNA V3-V4区进行PCR扩增。PCR反应体系（20 μL）为上下引物（10 μmol/L）各0.8 μL、1 μL 10 ng/μL模板DNA、2 μL 2.5 mmol/L dNTPs、0.4 μL Fastpfu聚合酶和4 μL 5×Fastpfu缓冲液，加ddH2O至20 μL。扩增程序为95 ℃预变性3 min；95 ℃变性30 s、55 ℃退火30 s、72 ℃延伸45 s，共30个循环；将相同体积的1×缓冲液（含SYBR Green Ⅰ ）与PCR产物混合，在2%琼脂糖凝胶上进行电泳检测。选择条带亮度在400~450 bp的样品进行试验。PCR产物以等密度比混合。使用德国桥根凝胶提取试剂盒对PCR产物进行纯化。在诺禾致源生物技术公司的Illumina MiSeq平台上对细菌扩增子进行测序。1.3.2序列及多样性分析基因序列文库使用离子加片段库工具包48 rxns（Thermo Scientific）生成，建库质量采用Qubit@2.0荧光计（Thermo Scientific）评估。在离子S5TM XL平台上对该文库进行测序，得到400 bp/600 bp的单端读数。根据样本的唯一条形码将单端读取分配给样本，并通过切断条形码和底漆序列来截断。根据Cutadapt[12]质量控制过程，在特定的过滤条件下对原始读数进行质量过滤，以获得高质量的干净读数。使用UCHIME算法（UCHIME Algorithm）[13]将读取的数据与参考数据库（Silva database）[14]进行比较，以检测嵌合体序列，去除嵌合体序列[15]，得读数。采用Uparse软件（Uparse v7.0.1001）进行序列分析[13]。相似性≥97%的序列被分配给相同的OTU（操作分类单元）。为更好地注释每个OTU的代表序列，先对每个OTU的代表序列进行筛选。对于每个有代表性的序列，Silva数据库[14]基于Mothur算法用于注释鉴定信息。为研究不同OTU的系统发育关系以及不同样本（群）中优势种的差异，利用MUSCLE软件（3.8.31版）[16]进行多序列比对。OTUs丰度信息使用与序列最少的样本对应的序列号标准进行归一化。α多样性和β多样性分析均基于这一输出标准化数据进行。应用α多样性方法，通过Observed-species、Chao 1、Shannon、Simpson、ACE、Good-coverage等6个指标，分析样本物种多样性的复杂性。这些指标均使用QIIME（1.7.0版）计算，并用R软件（2.15.3版）显示。使用β多样性分析评价物种复杂度的差异，使用QIIME软件（1.7.0版）计算加权和未加权Unifrac的β多样性。在聚类分析之前先进行主成分分析（PCA），在R软件（版本2.15.3）中使用FactoMineR包和ggplot2包降低原始变量的维数。通过主坐标分析（PCoA）从复杂的多维数据中获取主坐标并可视化。将之前得到的样本间加权或未加权均匀度的距离矩阵变换为1组新的正交轴，通过这些正交轴，第一主坐标表示最大变异系数，第二主坐标表示第二最大变异系数，依此类推。PCoA分析由R软件（2.15.3版）中的WGCNA包、stat包和ggplot2包显示。采用算术平均聚类法（UPGMA）对距离矩阵进行层次聚类，并用QIIME软件（1.7.0版）进行聚类分析。1.3.3细菌代谢途径与功能预测Tax4Fun功能预测是基于最小16S rRNA序列相似性的最近邻法，通过提取KEGG（京都基因与基因组百科全书）数据库原核全基因组16S rRNA基因序列，并使用BLASTN算法（BLAST Bitscore1500）将其与SILVA SSU Ref NR数据库对齐来实现的建立相关矩阵，将UProC和PAUDA注释的KEGG数据库的原核全基因组功能信息映射到SILVA数据库，实现SILVA数据库功能注释。利用SILVA数据库序列作为参考序列，将已测序的样本从OTU中聚类出来，获得功能注释信息。1.4数据统计和分析数据采用Excel进行初步整理，采用SPSS 20.0软件进行数据独立性、正态性和方差齐性检验，进行ANOVA单因素方差分析，采用Duncan氏法进行多重比较。数据结果以“平均值±标准差”表示，P0.05表示差异显著。2结果与分析2.1冷暖季节放牧牦牛粪便细菌Mothur与α多样性分析结果（见表3、图1）在3%的距离上，从20个不同的样本总共获得1 999个OTU，每个样本的平均OTU为1 079.95个。试验共分析20个不同样本的942 695个有效序列。由表3可知，在BTC组牦牛（冷季）瘤胃细菌的OUT数目显著高于BOC组（暖季）（P0.05），DTC组牦牛（冷季）在的OUT数目极显著高于在DOC组（暖季）（P0.01）。冷季放牧牦牛粪便菌群的ACE指数、Chao 1指数、Shannon指数和Simpson指数均显著或极显著高于暖季（P0.05或P0.01）。10.13557/j.cnki.issn1002-2813.2021.15.017.T003表3冷暖季节放牧牦牛粪便细菌的OTU数目和α多样性分析项目1.5岁P值3.5岁P值BOCBTCDOCDTCOTUs1 017.00±31.88b1 123.40±30.63a0.0301049.20±37.60b1 130.20±28.40a0.005测序深度指数1.00±01.00±0—1.00±01.00±01.000ACE指数1 004.06±28.65b1 118.27±24.80a0.0001 053.21±40.53b1 133.06±14.20a0.003Chao 1指数998.01±27.76b1 108.68±23.56a0.0001 064.86±35.85b1 139.15±37.32a0.004Shannon指数6.99±0.15b7.95±0.13a0.0007.18±0.11b7.82±0.22a0.000Simpson指数0.97±0.05b0.99±0.02a0.0000.98±0.04b0.99±0.04a0.011注：同年龄同行数据肩标不同字母表示差异显著（P0.05），相同字母或无字母表示差异不显著（P0.05）。具有97%相似性的稀疏曲线表明，足够的测序覆盖率可以解释每个样本中的大多数细菌多样性。由图1可知，覆盖指数大于99%，表明微生物群落反映准确。10.13557/j.cnki.issn1002-2813.2021.15.017.F001图1冷暖季节放牧牦牛粪便细菌的稀疏曲线图2.2冷暖季节放牧牦牛粪便细菌在门水平的优势菌群分析（见图2、图3）采用鉴定学分析方法，共检测到20个菌门。由图2可知，各组牦牛粪便细菌中，最丰富的门包括厚壁菌门（Firmicutes）、拟杆菌门（Bacteroidetes）、螺旋体门（Spirochaetes）、疣微菌门（Verrucomicrobia）、变形菌门（Proteobacteria）等。BOC组和BTC组以厚壁菌门（79.96%和69.74%）、拟杆菌门（18.99%、25.35%）、螺旋体门（0.01%、0.41%）、疣微菌门（0.11%、1.09%）和变形菌门（0.41%、0.72%）最常见。DOO组和DTC组中厚壁菌门（75.07%、59.97%）、拟杆菌门（23.22%、33.96%）、螺旋体门（0.01%、2.71%）、疣微菌门（0.07%、0.47%）和变形菌门（1.15%、0.66%）为主。10.13557/j.cnki.issn1002-2813.2021.15.017.F002图2冷暖季节放牧牦牛粪便细菌门水平相对丰度图由图3可知，与暖季组相比，冷季组牦牛粪便中拟杆菌门的相对丰度显著增加，厚壁菌门的相对丰度则显著降低。在1.5岁牦牛粪便中，在冷季中厚壁菌门的相对丰度极显著低于暖季（P0.01），而螺旋体门、疣微菌门、变形菌门和拟杆菌门在冷季中显著或极显著高于暖季（P0.05或P0.01）。在3.5岁牦牛粪便中厚壁菌门在冷季中显著低于暖季（P0.05），螺旋体门、疣微菌门显著或极显著增高（P0.05或P0.01），拟杆菌门和变形菌门无显著变化（P0.05）。10.13557/j.cnki.issn1002-2813.2021.15.017.F003图3冷暖季节放牧牦牛粪便细菌门水平热图2.3冷暖季节放牧牦牛粪便细菌在属水平的优势菌群分析（见图4、图5）在属水平中，共检测到145个细菌属。由图4可知，在BOC和BTC组中，最丰富的菌属为拟杆菌属（Bacteroides）（3.98%、4.07%）、阿里叶柄属（Alistipes）（2.02%、2.23%）、未分类-瘤胃球菌属（unidentified_Ruminococcaceae）（1.21%、1.47%）、阿克曼菌属（Akkermansia）（0.11%、1.19%）、Romboutsia（0.36%、2.38%）和霍氏真杆菌属（Faecalibacterium）（1.25%、0.01%）。在DOC和DTC组中，最常见的为拟杆菌属（9.08%、4.44%）、阿里叶柄属（1.60%、2.72%）、未分类-瘤胃球菌属（1.37%、1.47%）、Romboutsia（0.31%、2.69%）、霍氏真杆菌属（2.16%、0.01%）和阿克曼菌属（0.07%、0.46%）。10.13557/j.cnki.issn1002-2813.2021.15.017.F004图4冷暖季节放牧牦牛粪便细菌属水平相对丰度图在门水平上，厚壁菌门的相对丰度在冷季极显著低于暖季（P0.01）；在属水平上，阿克曼菌属、Romboutsia等菌属在冷季显著或极显著高于暖季（P0.05或P0.01）。由图5可知，在1.5岁牦牛粪便中，冷季时阿克曼菌属、Romboutsia、未分类-瘤胃球菌属等菌属极显著高于暖季（P0.01），霍氏真杆菌属显著低于暖季（P0.05）。在3.5岁牦牛粪便中，在冷季阿克曼菌属、Romboutsia等显著或极显著高于暖季（P0.05或P0.01），而拟杆菌属和霍氏真杆菌属显著低于暖季（P0.05）。10.13557/j.cnki.issn1002-2813.2021.15.017.F005图5冷暖季节放牧牦牛粪便细菌属水平热图2.4冷暖季节放牧牦牛粪便细菌的聚类差异分析（见图6）由图6可知，暖季与冷季牦牛粪便细菌群落具有很大差异。冷季组与暖季组的粪便细菌群落有明显的PC1差异，占总变异的57.17%，PC2差异占总变异的21.23%。10.13557/j.cnki.issn1002-2813.2021.15.017.F006图6冷暖季节放牧牦牛粪便细菌的主坐标分析（PCoA）图2.5冷暖季节放牧牦牛粪便细菌代谢途径和功能预测比较分析（见图7~图10）由图7可知，KEGG水平1鉴定代谢功能显示出最高的相对丰度，在每组中超过总阅读量的44%（BOC：44.15%；BTC：44.83%；DOC：44.27%；DTC：45.20%）。10.13557/j.cnki.issn1002-2813.2021.15.017.F007图7冷暖季节放牧牦牛粪便细菌KEGG水平1功能预测由图8可知，KEGG水平2鉴定，大多数基因属于膜转运、碳水化合物代谢（carbohydrate metabolism）、氨基酸代谢（amino acid metabolism）、复制修复（replication and repair）、翻译（translation）、能量代谢（energy metabolism）和核苷酸代谢（nucleotide metabolism）。10.13557/j.cnki.issn1002-2813.2021.15.017.F008图8冷暖季节放牧牦牛粪便细菌KEGG水平2功能预测由图9可知，在BOC和BTC组的比较中，碳水化合物代谢、氨基酸代谢、辅酶和维生素代谢（metabolism of cofactors and vitamins）、聚糖生物合成与代谢（glycan biosynthesis and metabolism）等功能在冷季的表达量显著或极显著高于暖季（P0.05或P0.01），而复制修复、信号转导（signal transduction）、细胞运动（cell motility）等则显著极低于暖季（P0.01）。10.13557/j.cnki.issn1002-2813.2021.15.017.F009图9使用PICRUSt测定的BOC和BTC组的KEGG路径比较由图10可知，在DOC和DTC组中，碳水化合物代谢、辅酶和维生素代谢、聚糖生物合成与代谢和运输和分解代谢（transport and catabolism）在冷季的表达量显著或极显著高于暖季（P0.05或P0.01），复制修复、翻译、信号转导等则显著或极显著低于暖季（P0.05或P0.01）。10.13557/j.cnki.issn1002-2813.2021.15.017.F010图10使用PICRUSt测定的DOC和DTC组的KEGG路径比较3讨论反刍动物后肠约占总消化道的17%，其活动可以为生物体提供5%~30%的能量，说明后肠菌群对营养物质的消化吸收起着重要作用，能够在一定程度上反映宿主的能量利用效率[17]。有试验表明，小尾寒羊后肠聚集了大量产生短链脂肪酸（SCFA）的细菌，是产生SCFA的重要区域[18]，类似的细菌也常见于牦牛和其他反刍动物肠道内[19]。这些菌群在不同反刍动物肠道中的功能作用并不相同。一些反刍动物可以为宿主转化更多的SCFA，而另一些反刍动物可以产生甲烷[20]。本试验结果表明，厚壁菌门和拟杆菌门放牧牦牛粪便细菌区系结构中的优势菌门，与Nie等[21]测定的牦牛肠道菌群结构的结果相似，表明利用粪便样本研究牦牛后肠菌群的结构和功能是合理的。大量研究表明，厚壁菌门和拟杆菌门在反刍动物和其他动物的胃肠道中的主要优势菌[22-23]。厚壁菌门是促进宿主胃肠道微生物降解纤维素最重要的菌门[24]，拟杆菌门是促进机体降解低聚糖最重要的菌门[25]。本试验中，冷季放牧牦牛粪便中的厚壁菌门含量明显增多而拟杆菌门含量明显减少，可能因为冷季天然牧草数量及质量严重不足。在属水平研究中发现，冷季牦牛粪便中嗜黏蛋白阿克曼氏菌属、Romboutsia和霍氏真杆菌属的相对丰度明显上升。有研究表明，Romboutsia具有益生菌特性，可以参与肠道代谢反应等[26-27]。阿克曼菌属可以利用黏蛋白作为能量来源，通过竞争防止病原体侵害肠道[28]。阿克曼菌属可能参与恢复肠道黏蛋白储备的机制[29]。阿克曼菌属降解黏蛋白时，会释放出各种副产物，包括乙酸等短链脂肪酸。肠道中阿克曼菌属的数量可能与乙酸盐含量有关[30]。由此可知，放牧牦牛在冷季时，为应对外界寒冷的环境及牧草的不足，会主动调节自身肠道免疫功能及能量供应机制。拟杆菌属和阿里叶柄属在暖季至冷季呈上升趋势。有研究表明，拟杆菌属在维持肠道菌群结构平衡方面发挥重要作用[31-32]，阿里叶柄属在提高机体免疫功能方面发挥重要作用[33]。青藏高原寒冷季节气温极低，牦牛特别是幼龄牦牛易患肠道疾病。拟杆菌属数量的增加，表明放牧牦牛在冷季能够主动调节并维持肠道生态平衡。本试验中，碳水化合物代谢、氨基酸代谢、辅酶和维生素代谢、聚糖生物合成与代谢等功能基因的数量在冷季显著高于暖季。放牧牦牛肠道菌群结构和健康会受到冷暖季节的变化的影响。本试验中，通过Alpha多样性分析发现，冷季中放牧牦牛粪便细菌中ACE、Chao 1、Shannon和Simpson指数均显著高于暖季。根据冷季和暖季放牧牦牛肠道菌群的PCoA分析结果，冷季和暖季牦牛肠道菌群群落有明显差异，但同一年龄和季节牦牛肠道菌群群落变化不大。冷季牦牛肠道细菌的丰富度和多样性显著增加，说明肠道菌群结构随冷季而变化，有利于牦牛适应高原冷季的寒冷环境。4结论冷暖季节的变化是影响放牧牦牛肠道菌群结构和健康的重要因素。从暖季到冷季，牦牛粪便菌群的丰富度和多样性均显著增加。放牧牦牛在冷季为了应对外界寒冷的环境及牧草的不足，会主动调节自身肠道免疫功能及能量供应机制。