应用16S rDNA分析青海喜马拉雅旱獭肠道菌群多样性

徐 鹏，张雪飞，马 英

青海是多种野生动物和珍稀保护动物的栖息地，生活着多种在中国乃至世界上都特有的物种，如喜马拉雅旱獭(Himalayan marmot)、高原鼠兔(Ochotona curzoniae)、藏羚羊(Tibetan antelope)等。1954年青海省从喜马拉雅旱獭体内分离出鼠疫菌，被判定为喜马拉雅旱獭鼠疫自然疫源地[1]。自此，喜马拉雅旱獭就成为鼠疫监测的重点对象，鼠防人员对其进行研究，发现喜马拉雅旱獭对鼠疫菌的留存、传播和流行发挥重要作用。近年，国内学者对青海部分地区的喜马拉雅旱獭所携带的病毒方面开展初步研究，相继在喜马拉雅旱獭粪便中发现了一些新的病毒[2-3]，说明我们对于该物种的研究还有很多的空缺需要去补充。

随着近年来国家加大对野生动物的保护力度，生态环境的日益改善，人们对于保护野生动物思想观念的提升，多种因素使得野生动物的数量迅速增多。另一方面，城市化进程的加剧，使人类活动和居住的范围越来越大，动物的栖息地却越来越小，野生动物不得不与人类一起共同生存，增加了人类与野生动物的接触，从而为病毒、细菌的跨物种传播提供了大量机会，通过野生动物的尿液、粪便或蜱类、螨类、蚤类等将贝纳柯克斯体、汉坦病毒、鼠疫菌等传播给人类[4]。现已确认的多种急性传染病中，来源于野生动物的占比高达43%[5]。于是，野生动物间、野生动物与人之间的传染病防治就成为了目前主要的防控要点。而肠道作为生物体内一个庞大、复杂的微生物寄生区域，蕴藏着丰富的细菌和病毒，对其进行细致深入的研究势必成为防控传染病的重要一环。

喜马拉雅旱獭作为一种大型地栖类啮齿动物，广泛分布在青藏高原以及与中国接壤的尼泊尔等国的青藏高原边缘山地。但目前的研究仅局限于传统形态分类、物种生态调查和动物鼠疫监测等方面，针对携带其它病原微生物的研究很少，肠道菌群微生物多样性方面的研究目前还未见报道。因此，我们拟对青海不同地理区域的喜马拉雅旱獭肠道菌群展开研究，为应对、预防和控制未来可能发生的动物源性新发传染病提供基础数据，同时也为鼠疫的防控工作拓展新思路。

1 材料与方法

1.1 实验材料

1.1.1 样品采集 采样点结合青海鼠疫流行重点样区、旱獭种群遗传特点和动物地理区划，选择青海省3个地理亚区的5个小区共采集45份青藏高原喜马拉雅旱獭粪便样本(使用一次性粪便采集器将新鲜粪便样本挑取放入5 mL冻存管，并快速置于液氮罐中冷冻保存)。按照不同的采样地区将45份样本分为14组，采样点详细信息见表1。

表1 采样点地理信息Tab.1 Geographical information for sampling points

1.1.2 样品处理 称取200 mg样品，放入灭菌的2 mL离心管中，加入1 mL 70%乙醇，震荡混匀，10 000 r/min室温离心3 min，弃置上层液体。加入1×PBS溶液，震荡混匀，10 000 r/min室温离心3 min，弃置上层液体。倒置2 mL管于吸水纸上1 min，直至没有液体流出。将样品管放入55 ℃烘箱10 min，使残留酒精完全挥发，保证后续实验操作。

1.2 数据处理

1.2.1 数据质控 测序后得到原始测序数据(Raw PE)，为了得到能够用于后续分析的数据，首先去除干扰数据，通过PEAR(V0.9.6)将序列进行拼接[9]，去除Barcode和引物序列得到Raw Tags[8]；依据Qiime(V1.7.0)对Raw Tags的质控过程[10]，分别进行Tags截取和Tags长度过滤等操作后得到Clean Tags；Clean Tags经过UCHIME Algorithm(V4.2.40)去除嵌合体等步骤[11]，得到最终的有效数据(Effective Tags)[12]。

1.2.2 OTU聚类与物种注释 利用Uparse软件[13]对所有样品的全部Effective Tags进行聚类，以97%的一致性归为一类，即一个可操作分类单元(Operational Taxonomic Units，OTUs)，其中无法被聚类的Tags序列和Tags序列频数为1的序列将被舍弃，Tags 序列中出现频数最高的将被作为 OTUs 的代表序列。使用Mothur(V1.30.1)[10]与SILVA的SSUrRNA数据库对其进行物种注释，得到物种在域(domain)、门(phylum)、纲(class)、目(order)、科(family)、属(genus)上的构成。

1.2.3 微生物多样性分析 基于各样本数据的均一化处理，探究各个样品的复杂度，采用Alpha多样性分析。运用Miseq SOP[14]绘制等级丰富度图(Rank Abundance)展现样本的物种丰富度与均匀度、稀释曲线反映测序深度是否合适；计算Chao 1与ACE指数[15]评估菌群丰富度、计算Shannon和Simpson指数[16]评估菌群多样性。

1.2.4 物种分类与丰度聚类 为了得到每个OTU对应的物种分类信息，需要对OTU进行物种分类，我们采用RDP classifier方式。RDP classifier(V2.12)[6]基于Bergey’s taxonomy，采用Naive Bayesian assignment算法对每条序列在不同层级水平上计算其分配到此rank中的概率值。Bergey’s taxonomy分为6层，它们依次为domain、phylum、class、order、family、genus，根据其制作物种分类表。物种分类数据库采用NCBI 16S(http://ncbi.nlm.nih.gov/)。根据分类学分析结果，统计在各个分类层级水平上的每个样品的群落组成。

基于物种相对丰度的样本聚类树图可以通过树枝结构直观的反应出多个样品间的相似性和差异关系。首先使用R(V3.2)根据各样本物种相对丰度计算beta多样性距离矩阵，使用Bray Curtis法获得样本间距离后进行层次聚类(Hierarchical clustering)分析，再使用非加权组平均法UPGMA(Unweighted pair group method with arithmetic mean)算法构建树状结构，得到树状关系形式用于可视化分析。

1.2.5 肠道菌群与环境因子的关系 根据不同的采样地区将45份样本分为14组，在R上先使用species-sample数据(97%相似性的样品OTU表)做DCA(Detrended Correspondence Analysis)分析，根据分析结果中Lengths of gradient的第一轴的大小选择冗余分析(Redundancy Analysis，RDA)或典范对应分析(Canonical Correspondence Analysis，CCA)。RDA或者CCA是基于对应分析发展而来的一种排序方法，将对应分析与多元回归分析相结合，每一步计算均与环境因子回归，主要是用来反映每一组样本的肠道菌群与环境因子(温度、海拔、经度、纬度)之间的关系。

为探究环境因子具体对肠道哪一种菌群产生影响，将肠道菌群在属水平含量上从高到低排序，通过使用R的pearson correlation算法计算相关性系数(coefficient)来分析前20个含量较高的细菌与海拔、温度、经度、纬度的相关性，最终获得显著影响某些物种的环境因子。其中P值越小相关性越大，表明该细菌属与该环境因子的关系较为密切。

1.2.6 菌群相对丰度差异及功能预测分析 基于物种分类结果，计算在不同水平上各rank的丰度，比较样本间相对丰度差异，找出由于生存环境不同而存在显著差异的物种分类，关注样本中存在明显差异的菌属。筛选条件为P≤0.05，比较对象为样本间两两比较，采用Fisher’s精确检验(Fisher exact test)，最后将检验得到的P值采用FDR(False Discovery Rate)做Multiple test correction得到q值。

利用PICRUSt软件(V1.0.0)[17]对喜马拉雅旱獭肠道菌群的功能进行预测分析，基于COG基因预测结果对所有样品进行聚类与柱状图组合分析，样本聚类树图可以通过树枝结构直观的反应出多个样品间的相似性和差异关系。

2 结 果

2.1 Alpha多样性分析 通过数据拼接、质控和嵌合体过滤等步骤，将45份喜马拉雅旱獭粪便样本进行16S rDNA基因测序分析，所有样本菌群丰富度与菌群多样性均符合实验要求，覆盖率除XG3-18为96.55%、GD1-18为97.78%、WT1-18为97.79%以外，其余均高于98%。Rank Abundance曲线较宽且平坦，说明物种组成丰富且均匀度高；稀疏曲线趋于平稳且其中42份样本覆盖率均高于98%，说明样本测序深度足够[18]。

2.2 肠道菌群物种分类 从门分类水平上分析喜马拉雅旱獭的肠道菌群，排名前6位的分别是：变形菌门(Proteobacteria)、厚壁菌门(Firmicutes)、拟杆菌门(Bacteroidetes)、放线菌门(Actinobacteria)、疣微菌门(Verrucomicrobia)、软壁菌门(Tenericutes)。其中变形菌门占比最高，为38.77%；软壁菌门占比最低，为0.15%。

在属水平上排名前10位的包括：不动杆菌属(Acinetobacter)、节杆菌属(Arthrobacter)、鞘氨醇杆菌属(Sphingobacterium)、埃希氏菌属(Escherichia)、另枝菌属(Alistipes)、黄杆菌属(Flavobacterium)、瘤胃球菌属(Ruminococcus)、肉杆菌属(Carnobacterium)、梭菌属 IV(ClostridiumIV)、丛毛单胞菌属(Comamonas)，其中不动杆菌属占比最高，为20.36%，是样本WT2-18、CD11-18、GH11-18、WK8-18中的优势菌属；丛毛单胞菌属占比最低，为1.94%，见表2。

表2 样品分类系统组成表Tab.2 Composition of all sample classification systems at the genus level

2.3 物种相对丰度聚类 根据各样本物种相对丰度计算beta多样性距离矩阵，使用Bray Curtis法获得样本间距离后进行层次聚类得到3个聚类，第一聚类又分为两个类别，第一类别不动杆菌属(Acinetobacter)较高，不动杆菌属为革兰氏阴性专性需氧菌，第二类别鞘氨醇杆菌属(Sphingobacterium)较高，其属于好氧或兼性厌氧非发酵革兰氏阴性杆菌；第二聚类中未分类菌(unclassified)含量较高，并且是在属分类状况下，表明其中含有大量新奇且未被充分研究的新微生物；第三聚类由于样本太少(2份)，故不做分析，结果如图1所示。

图分为左、中、右3个部分，图左：基于Bray-Curtis的样本聚类树图；图中：根据聚类顺序排序的物种丰度柱状图；图右：物种的图示。图1 基于物种相对丰度的所有样品聚类树图与柱状图组合分析图Fig.1 Combination analysis graph of clustering dendrogram and histogram for all samples based on relative abundance of species

2.4 样本与环境因素的RDA分析 根据DCA的分析结果，Lengths of gradient的第一轴<3，于是对收集到的样本环境因子信息与样本采用RDA分析，可得：温度(T)≈海拔(ASL)>经度(E)>纬度(N)，如图2所示。T对1组与2组影响较大。ASL对8组、3组、7组、9组、10组、6组影响较大。E对12组、4组、1组和2组影响较大，N对各组的影响较小。海拔与温度对肠道菌群的影响最大，海拔、温度不同，肠道菌群的组成也不同[19-22]。通过细菌属与环境因子的相关性分析后发现不动杆菌属(Acinetobacter)与海拔、温度、经度关系密切；瘤胃球菌属(Ruminococcus)与温度、经度相关；马赛菌属(Massilia)与海拔、纬度相关；肉杆菌属(Carnobacterium)与乳杆菌属(Lactobacillus)仅与温度相关；假单胞菌属(Pseudomonas)仅与海拔相关。

图2 14组样本的RDA图Fig.2 RDA chart of 14 sets of samples

2.5 样本间菌群相对丰度差异及功能预测分析 探究单个样本由于生存环境不同，肠道菌群中哪些细菌属存在差异，使用Fisher exact test检验后发现45个样本中有14个属的细菌存在统计学差异(表3)。

表3 45个样本中存在明显差异的细菌Tab.3 Bacteria with obvious differences in 45 samples

为探究样本细菌功能，利用PICRUSt软件，基于COG的功能结构对喜马拉雅旱獭肠道菌群的功能进行预测分析，结果显示喜马拉雅旱獭肠道菌群功能预测大多集中在遗传信息处理与代谢通路上，如基因一般功能预测，未知功能转录、复制、重组与修复，核糖体结构翻译与生物发生，氨基酸运输与代谢，能源生产与转换，细胞壁与细胞膜的生物发生。表明同一物种，不同生存环境，其肠道微生物群落的功能存在一定的差异。其中有6个样本(XH16-18、XH12-18、XH15-18、GH11-18、GM6-18、TR89-18)整体功能丰度较低且结构较为相似，且对照采样点信息，这6个样本都不处于喜马拉雅旱獭的最适生存环境[23]；整体功能丰度最高的3组样本(XG3-18、DH3-18、XB3-18)的生存环境均处于青藏高原喜马拉雅旱獭的最适生存环境，结果如图3所示。

图分为左、中、右三个部分，图左：基于Bray-Curtis的样本聚类树图；图中：根据聚类顺序排序的功能丰度柱状图；图右：物种的图示。 图3 基于COG的所有样品聚类树与柱状图组合分析图Fig.3 Combination analysis graph of the clustering tree and histogram based on COG for all samples

3 讨 论

喜马拉雅旱獭样本中分类操作单元(OTU)数目越多，表示细菌分类越多样，反之则分类越单一。本研究中第一聚类的样本平均OTU要远小于第二聚类，提示第一聚类的样本所面临的环境复杂程度要小于第二聚类，回顾采样点的生境条件，第二聚类样本的环境复杂度确实要高于第一聚类的样本，但之前没有人进行过相关的实验研究，所以我们目前还只是推测，具体的验证工作有待于开展下一步的实验。

Ley等[24]利用宏基因组技术对来自13个生物分类的60种哺乳动物的肠道菌群分析后发现相对丰度排序前6的门为：厚壁菌门(Firmicutes)、拟杆菌门(Bacteroidetes)、变形菌门(Proteobacteria)、放线菌门(Actinobacteria)、疣微菌门(Verrucomicrobia)、梭杆菌门(Fusobacteria)；冯天舒等[25]对高原鼢鼠的盲肠内容物进行16SrRNA基因测序后获得该物种的肠道菌群组成，其主要的肠道微生物有：拟杆菌门(Bacteroidetes)、厚壁菌门(Firmicutes)、变形菌门(Proteobacteria)、螺旋体门(Spirochaetes)、放线菌门(Actinobacteria)；谭春桃等[26]对青藏高原夏季野生鼠兔的肠道菌群进行研究后发现其肠道菌群主要集中在拟杆菌门(Bacteroidetes)、变形菌门(Proteobacteria)、厚壁菌门(Firmicutes)、放线菌门(Actinobacteria)、蓝藻菌门(Cyanobacteria)、螺旋体门(Spirochetes)。可以看出喜马拉雅旱獭在肠道菌群的组成上与上述研究的哺乳动物肠道菌群存在一定的共性，但是也存在一定的区别，说明不同种类的动物，肠道菌群的组成比例完全不同，同一种细菌可能在一种动物肠道中占比较高，而在另外一种动物肠道中占比极少。

本文中喜马拉雅旱獭肠道菌群分布以变形菌门的含量为最高，其次是厚壁菌门、拟杆菌门。变形菌门不仅是喜马拉雅旱獭，也是大熊猫肠道中的主要菌群[27]。喜马拉雅旱獭作为植食性动物，禾木科和莎草科植物的绿色部分为主要饲料，间或食种子及花序，变形菌门细菌通过降解动物食物中的木质素，破坏植物细胞壁，使动物快速汲取营养物质；厚壁菌门作为纤维素分解细菌，可将植物纤维降解为挥发性脂肪酸以供宿主利用；拟杆菌门则主要帮助宿主降解碳水化合物、蛋白质，提高宿主营养利用率，维持肠道的微生物生态平衡，促进宿主的免疫系统发育，增强宿主免疫水平[28-30]。各种菌属协调共生，相辅相成，共同提高宿主生命健康水平。

有研究表明青藏高原喜马拉雅旱獭的最适生长温度为14.3～24 ℃，最适生长高度为3 000～4 000 m[23]，基于Bray-Curtis建立的样本聚类树图中，第一聚类第一分类中不动杆菌属含量较高，这是一种专性需氧菌，高寒低氧环境可使其活力下降。第二分类中鞘氨醇杆菌属含量较高，该种细菌属含有大量的鞘磷脂[31]，是细胞膜的主要组成成分。综合采样信息，第一聚类中有大量样本不处于喜马拉雅旱獭最适生存环境，为了适应高海拔给机体带来的影响，不动杆菌属和鞘氨醇杆菌属以相对较高的含量来应对高海拔所造成的低氧环境；第二聚类的喜马拉雅旱獭肠道菌群中未分类菌含量较多，究其原因是之前很少全面研究高寒低氧环境下野生动物的肠道微生物菌群，我们猜测在这一聚类中遇到了某些较为新奇、尚未被充分研究的微生物，这些微生物可能在宿主处于极端环境下发挥作用，但实际情况是否如此，有待于后续进一步的深入研究。

对样本丰度存在显著差异的物种进行t检验，得到14种有统计学意义的细菌属别，我们猜测这些差异可能与样本的肠道菌群的功能分布有关。因为同一物种处于不同环境之下，为了促进宿主适应环境，肠道菌群也随之发生了不同程度的结构、组成改变，从而肠道菌群的功能分布也发生变化，以达到一种自然状态下的保护性机制[32]。并且考虑采样时间，该阶段喜马拉雅旱獭处于交配期，存在交配竞争、食物竞争、领地竞争等现象。为了在竞争中获得有利地位，必然会通过多种途径提升自己的体重，摄入多元充足的食物，使身体各项机能指标达到最佳状态[33]，从而导致其相应肠道菌群的功能丰度上升。

PICRUSt分析发现喜马拉雅旱獭的肠道菌群功能预测多集中在遗传信息处理与代谢通路上，说明低氧可能对宿主肠道微生物的基因造成了损伤，这与马燕等[34]的研究结果一致。高寒缺氧带给喜马拉雅旱獭的不仅仅是生存压力，其繁殖压力也随之而来[35-36]。本研究有64.45%的样本含有不动杆菌，该种细菌为机会致病菌，在机体免疫力低下时导致机体疾病的发生，加之暴露于高寒低氧条件本身就会导致动物疾病的发生机会增加[37-38]，这无疑严重影响了喜马拉雅旱獭的生存繁衍。但正因这些恶劣的环境给我们提供了独特的研究机会，研究高原野生动物肠道微生物菌群多样性组成，可为微生物学、生态学和兽医学提供信息数据，进而对公共卫生事业的进步具有一定的促进作用。