腹泻柴达木马亚成体粪便微生物多样性分析及其生物标记物的筛选

王小琪 郝文静 张志超 韩晶,4 王儒敬 段子渊*

（1中国科学院遗传与发育生物学研究所，北京 100101）（2中国科学院合肥物质研究院，合肥 230031）（3中国科学技术大学，合肥 230026）（4宁夏大学农学院，银川 750021）

腹泻是马(Equus caballus)常见的一种消化道疾病，严重影响马的健康和正常活动，尤其对未成年马的健康状况影响极大。但动物腹泻的病因复杂，包括结肠炎、寄生虫感染、病毒或细菌感染等 (Oliver-Espinosa,2018;Shaw and Stämpfli,2018)，临床症状不一致。肠道微生物影响宿主营养代谢和能量供应，还可能介导宿主疾病的发生发展(翟齐啸等，2013；杨利娜等，2014)。尽管引起马腹泻的原因并非全部都是微生物，但已证明大肠杆菌(Escherichia coli)、沙门氏菌(Salmonella)、嗜水气单胞菌(Aeromona hydrophyla)和新立克次氏体(Neorickettsia risticii)等均与马腹泻相关(Shaw and Stämpfl,2018;Costa and Weese,2018)。马肠道微生物群落的平衡稳态对胃肠疾病等非常敏感，发生腹泻的马肠道菌群种类和丰度会发生改变，而肠道菌群的紊乱也会引起腹泻，甚至导致死亡，两者密切相关(Al Jassim and Andrews,2009;Garberet al.,2020)，因此有必要从肠道菌群多样性方面分析马在腹泻期间肠道菌群的结构和组成。

柴达木马主要分布在青海省都兰、乌兰、格尔木和德令哈地区，目前相关研究极少，关于柴达木马肠道微生物群落多样性的研究未见报道，非常不利于柴达木马的遗传资源保护、开发与利用。采用常规的培养方法只能分离肠道或粪便的部分微生物，而16S rRNA高通量测序可以获得粪便中大部分微生物的信息，进而可以有效分析粪便微生物群落多样性的改变。本研究对发生腹泻和健康的柴达木马亚成体的粪便微生物群落进行了16S rRNA测序数据分析，对两者微生物群落组成及结构差异进行了比较，并利用机器学习的相关算法(随机森林法，random forest)对影响柴达木马亚成体肠道菌群的特征菌属进行了筛选。在青藏高原独特的自然条件下形成了有较好脂肪囤积能力的柴达木马，主要生活在我国柴达木盆地及周边地区，是对荒漠、沼泽盆地适应性极强的优良地方品种，也是当地牧民重要的生产对象(青海省质量技术监督局，2011)，因此，柴达木马发生腹泻会对当地居民的生产生活产生严重影响。虽然对于柴达木马腹泻的患病率并没有确切的报道，但根据当地的调研情况，腹泻一直是马医学科中最常见的疾病，并且动物腹泻的临床病例占比高达27%(Xieet al.,1997)，由此，通过与健康的柴达木马的粪便微生物进行对比分析，了解发生腹泻的柴达木马肠道细菌群落情况，为柴达木马的健康养殖提供参考。

1 研究方法

1.1 样品采集

柴达木马粪便样品采集自青海省乌兰县(海拔3 300 m左右)同一个群体的雄性亚成体马(1～2岁)。采样时人跟随马身后等待，一旦排泄，立即用无菌手套将未沾土的粪便装入无菌收集袋，并立即置于保温箱干冰之中，样品采集结束当天运输至北京实验室，-80°C保存。在该群体共采集到同时发生腹泻的3匹柴达木马粪便样本(粪便不成形)，设为腹泻组 (Diarrhea,Dia)，编号Dia1，Dia2和Dia3；同群体13匹健康柴达木马粪便样本，设为健康组(Naive control,NC)，编号NC1～NC13。

1.2 基因组DNA提取及16S rRNA测序

按照粪便基因组DNA提取试剂盒(DP328-02，天根，北京)说明书，从16个粪便样本中分别提取微生物基因组DNA。1%琼脂糖凝胶电泳，酶标仪(Epoch，BioTek，美国)检测提取DNA浓度和纯度，确保提取DNA的OD 260/280值为1.6～2.0，且浓度大于100 ng/μL。使用通用引物组(341F和805R)对细菌16S rRNA基因的V3～V4区域进行扩增(Herlemannet al.,2011)，扩增产物送至北京健云康生物科技有限公司建库和双端测序(Illumina HiSeq2500测序平台)。

1.3 实时荧光定量PCR检测菌属含量

普雷沃菌属(Prevotella)、拟杆菌属(Bacteroides)、甲烷短杆菌属(Methanobrevibacter)和双歧杆菌属(Bifidobacterium)的特异性PCR引物由华大基因(BIG)合成(表1)，并依据引物对应的参考文献中的方法进行PCR扩增(PTC-200，Bio Rad，美国)。2.5%琼脂糖凝胶电泳检测PCR产物，胶回收试剂盒(DC301-01，诺唯赞生物科技公司，南京)对其回收纯化。之后使用pMD19-T Vector(6013，Takara，美国)和TreliefTM5α感受态细胞(TSC01，擎科生物科技有限公司，北京)进行连接转化，得到阳性克隆。阳性质粒也由BIG进行测序，利用GenBank的BLAST功能对测序结果进行同源性分析。鉴定成功的菌液使用普通质粒小提试剂盒(RTP2102-03，中科瑞泰生物科技有限公司，北京)提取质粒DNA，酶标仪检测浓度。通过以下公式计算拷贝数：

表1 柴达木马粪便细菌16S rRNA基因PCR扩增引物Table 1 PCR amplification primers of bacterial 16S rRNA gene in Qaidam horses’feces

依照引物对应的参考文献中反应体系及方法，使用SYBR qPCR预混液(Q711，诺唯赞生物科技公司，南京)和实时荧光定量PCR仪(CFX Connect，Bio Rad，美国)对计算出拷贝数的质粒标准品进行5倍连续梯度稀释后，以每梯度3次重复建立标准曲线；对16个样品的普雷沃菌属、拟杆菌属、甲烷短杆菌属和双歧杆菌属进行荧光定量检测。

1.4 生物信息与统计分析

根据Wang等(2021)的分析方法，通过QIIME1.9.1对数据进行质量过滤等处理，用usearch61去除嵌合体序列，使用UCLUST将获得的上述有效序列进行聚类分析(以相似度97%为阈值)，生成可操作分类单元(Operational taxonomic units,OTU)表格，根据RDP分类器(Ribosomal database project,RDP)和Greengenes13.8数据库对OTU进行物种注释。QIIME1.9.1计算Alpha多样性 (Chao1和Shannon指数)。Alpha多样性指样本内的物种多样性，是反映物种丰富度和均匀度的综合指标。Shannon是菌群多样性指数，数值与群落多样性成正比；Chao1是菌群丰富度指数，指数越大，菌群丰富度越高(夏瑜等，2018)。Beta多样性可以反映组间的微生物群落差异(陈圣宾等，2010)，本研究基于Unweighted Unifrac距离进行了Beta多样性分析。通过R4.0.3的pheatmap包和stats包对数据进行中心化和标准化处理，complete聚类法绘制热图；利用random forest、pROC、vegan及ggplot2软件包筛选细菌标志物，进行ROC曲线分析、ADONIS分析和绘制科学图表。采用Student’st检验对Alpha多样性指标及细菌门、属的丰度进行组间差异比较，P＜0.05为差异显著，0.05＜P＜0.10认为两组间具有显著性差异趋向。基于OTU表格，使用BugBase对细菌菌群进行表型分类(Wardet al.,2017)。在OTU水平上计算非加权Unifrac(Unweighted unifrac)距离矩阵并进行主坐标分析(Principal co-ordinates analysis,PCoA)。代表每个样本的点与点之间在PCoA图上的距离越远，表明组间微生物群落结构组成差异越大。ADONIS通常用于分析不同分组因素对样品差异的解释度，并利用F检验进行显著性分析。

2 结果

2.1 粪便16S rRNA测序数据预处理统计

对原始数据进行去除低质量片段、拼接、去嵌合体等预处理后，得到用于后续分析的有效序列575 354条。腹泻和健康两个组的样本测序覆盖度均超过0.97，说明测序深度可以反映样品中细菌的物种情况(表2)。对有效序列以97%相似度聚类，鉴定出6 122个OTU，属于165个细菌的属。16S rRNA测序的原始数据已提交至国家基因组科学数据中心(https://ngdc.cncb.ac.cn/)的GSA数据库，编号CRA004854。

表2 柴达木马粪便细菌16S rRNA测序数据Table 2 16S rRNA sequencing data of feces from Qaidam horses

2.2 微生物群落Alpha多样性和Beta多样性分析

由图1A可知，健康组Shannon指数显著高于腹泻组(P＜0.05)，说明健康组微生物物种多样性高于腹泻组。此外，健康组Chao1指数虽略高于腹泻组，但无显著差异(P＞0.05)，说明两组物种丰富度相似。由此可见，两个组微生物物种多样性的差异主要来自于菌群均匀度的差异。

图1 柴达木马粪便微生物物种多样性.A：粪便微生物群落Alpha多样性(均值±标准差)；B：基于非加权Unifrac矩阵的主坐标分析图.*P＜0.05Fig.1 The fecal microbial diversity of Qaidam horses.A:The Alpha diversity of fecal bacteria(mean±SD);B:PCoA plot based on Unweighted Unifrac matrix.*P＜0.05

PCoA分析显示健康组和腹泻组的样本明显分开(图1B)，表明健康组和腹泻组的微生物群落结构组成存在差异。通过ADONIS分析，两组间的解释度较高(R2=0.14,P=0.007)，表明两组间的粪便微生物群落结构具有统计学差异。第一主坐标和第二主坐标的贡献率分别是25.81%和16.44%，较高的贡献率也增加了结果的可靠度。在PCoA1上两组样本可以完全分开，说明PCoA1是造成腹泻组和健康组微生物群落差异的主维度。

2.3 微生物群落结构组成分析

对OTU进行物种分类学归类，健康组和腹泻组样本粪便细菌和古菌相对丰度排名前十的菌门见图2A，两组菌门相同。厚壁菌门(Firmicutes)、拟杆菌门(Bacteroidetes)和疣微菌门(Verrucomicrobia)占主导地位，其相对丰度超过86.50%。相对于健康组，腹泻组亚成体马粪便中拟杆菌门、变形杆菌门(Proteobacteria)相对丰度更高(P＜0.05)，而厚壁菌门、疣微菌门、螺旋体菌门(Spirochaetes)、纤杆菌门(Fibrobacteres)的相对丰度较低(P＜0.05)。此外，两组间低丰度菌门，即浮霉菌门 (Planctomycetes)、迷踪菌门 (Elusimicrobia)、蓝细菌门(Cyanobacteria)、放线菌门(Actinobacteria)和广古菌门(Euryarchaeota)差异显著(图2B)。

图2 腹泻组和健康组柴达木马粪便微生物门水平的相对丰度比较.A：粪便细菌相对丰度组成(门水平)；B：两组间有显著性差异的菌门.*P＜0.05Fig.2 Comparison of relative abundance of fecal microbial in diarrhea(Dia)group and naive control(NC)group of Qaidam horses at phylum level.A:The histograms of microbiota abundances at phylum level;B:The phyla with significant difference between two groups.*P＜0.05

在健康组粪便中共检测到163个属的细菌，其中27个属相对丰度大于0.5%；腹泻组检测到149个属的细菌，相对丰度大于0.5%的属有25个。在相对丰度排名前15%的菌属中，鉴定出9个属，分别是厚壁菌门的瘤胃球菌属(Ruminococcus)、梭菌属(Clostridium)、颤螺菌属(Oscillospira)、土源芽孢杆菌属(Solibacillus)，拟杆菌门的普雷沃菌属(Prevotella)，螺旋体菌门的密螺旋体属(Treponema)，纤杆菌门的纤杆菌属(Fibrobacter)，变形杆菌门的反刍杆菌属(Ruminobacter)和疣微菌门的阿克曼氏菌属(Akkermansia)(图3A)。为了更清晰地展示腹泻组和健康组柴达木马亚成体粪便细菌在属水平上的差异，进一步比较了两组样本的全部菌属，并对其中有显著性差异(P＜0.05)的菌属进行热图分析(图3B)。结果显示，两组间有显著性差异的菌属的相对丰度为0.01%～5.00%，属于低丰度属。其中，健康组纤杆菌属、瘤胃球菌属、梭菌属、拟杆菌属和Paludibacter的相对丰度 (2.58%、2.90%、0.68%、0.23%和0.51%)远高于腹泻组(0.15%、1.90%、0.39%、0.12%和0.06%)；腹泻组甲烷短杆菌属(Methanobrevibacter)、Sphaerochaeta、丁酸弧菌属(Butyrivibrio)、迷踪菌属(Elusimicrobium)和月形单胞菌属 (Selenomonas)的相对丰度 (0.24%、0.12%、0.05%、0.13%和0.06%)远高于健康组(0.05%、0.03%、0.02%、0.01%和0.01%)。

图3 腹泻组和健康组柴达木马粪便微生物属水平的相对丰度比较.A：粪便细菌相对丰度组成(属水平)；B：热图展示的两组间有显著性差异的菌属Fig.3 Comparison of relative abundance of fecal microbial in diarrhea(Dia)group and naive control(NC)group of Qaidam horses at genus level.A:The histograms of microbiota abundance at genus level;B:Clustering heatmap of genera with significant difference between two groups

2.4 实时荧光定量PCR定量检测结果

通过构建相应质粒，利用荧光定量PCR对粪便DNA中的普雷沃菌属、拟杆菌属、甲烷短杆菌属和双歧杆菌属进行绝对定量检测(图4)。其中拟杆菌属和甲烷短杆菌属在腹泻组和健康组间差异显著，两组的拟杆菌和甲烷短杆菌拷贝数的比较与16S rRNA测序的结果一致。尽管双歧杆菌在肠道菌群中含量低，却是代表性的益生菌，由于本次研究样本中双歧杆菌含量较低，16S rRNA测序未检测出双歧杆菌，而实时荧光定量PCR结果揭示在腹泻柴达木马的肠道内双歧杆菌含量虽低，却高于健康组，但两组差异不显著。普雷沃菌是与宿主饮食、健康密切相关的菌属，根据PCR定量结果，与纤维消耗相关的普雷沃菌和拟杆菌的比值在健康组中略高，但两组差异不显著。

图4 腹泻组和健康组柴达木马粪便主要菌属16S rRNA基因拷贝数差异分析Fig.4 Compared average copy number of 16S rRNA gene from the main fecal bacteria between diarrhea group and naive control group of Qaidam horses

2.5 腹泻相关属筛选

通过随机森林分析旨在确定一个小的特征菌属集合，可以最大程度上解释腹泻组和健康组柴达木马间的差异，进而用于腹泻的诊断和病程检测。图5A显示利用随机森林筛选的健康组和腹泻组柴达木马亚成体粪便微生物差异形成中的重要特征菌属(生物标记物)。根据十倍交叉验证选择重要性和可靠度排名最高的12个菌属，已被鉴定的有甲烷短杆菌属、纤杆菌属、Paludibacter、肉食杆菌属(Carnobacterium)和迷踪菌属。较少的样本量对机器学习的准确性会有一定的影响，在此使用ROC曲线检验随机森林预测的准确性(图5B)。在筛选出的生物标记物中，超过9个标记物的AUC值大于0.85，说明其具有较好的分类效果。这些菌属作为潜在微生物标记物对健康组和腹泻组柴达木马亚成体粪便细菌间的差异具有显著性影响。

图5 基于随机森林的健康和腹泻柴达木马粪便微生物差异分析(A)及ROC曲线分析(B)Fig.5 Random Forest bar chart of variable importance on genus level(biomarkers)(A)and ROC curve analysis(B)

2.6 物种的功能分析

为进一步了解腹泻组和健康组间菌群表型的差异情况，使用BugBase对两组菌群进行表型分类比较。腹泻组和健康组粪便的厌氧细菌、好氧细菌和兼性厌氧细菌丰度无显著差异。腹泻组的革兰氏阴性细菌(P=0.057)和潜在致病性细菌(P=0.025)的丰度高于健康组，而革兰氏阳性细菌(P=0.057)和氧化胁迫耐受性细菌(P=0.039)的丰度则低于健康组(图6)。

图6 腹泻组和健康组柴达木马粪便微生物表型相对丰度Fig.6 The microbiota phenotype relative abundance in Qaidam horses of diarrhea(Dia)group and naive control(NC)group

3 讨论

本研究采集的粪便样本来自于青藏高原地区同一群体1～2岁龄的亚成体柴达木马，通过16S rRNA高通量测序技术，对其粪便样本的微生物群落结构多样性进行分析比较。其中，每组样本数量≥3，与多项利用16S rRNA测序技术探讨实验动物或野生动物粪便菌群多样性的研究相同(Donget al.,2020；李栋梁等，2021)，可有效且合理地支撑实验结果。通过Alpha多样性分析发现，相比健康组，腹泻组的粪便菌群多样性显著下降。相似的研究也表明，腹泻与比利时成年马、加拿大安大略省饲养场马驹的肠道菌群多样性、丰富度和均匀度下降有关(Rodriguezet al.,2015;Schosteret al.,2017)。通过Beta多样性分析可见，与大多数单胃动物腹泻报道相同(Zhuet al.,2018;Wanget al.,2020)，本研究涉及的腹泻组和健康组亚成体马的粪便样本间菌群结构差异极大，远高于组内差异。Wang等(2020)采用与本研究相同的测序技术研究了健康和腹泻的野猪(Sus scrofa)粪便微生物，发现与健康样本相比，腹泻野猪粪便微生物中拟杆菌门、厚壁菌门的丰度都有所下降，而变形杆菌门丰度上升。Costa等(2012)基于16S rRNA基因的V3～V5区序列分析，发现来自加拿大安大略省的健康成年马粪便中优势菌门是厚壁菌门，其次是拟杆菌门和变形杆菌门；腹泻马的粪便中，拟杆菌门丰度最高，其次是厚壁菌门和变形杆菌门。Rodriguez等(2015)基于16S rRNA基因V1～V3区序列数据，统计比利时成年马的粪便细菌丰度后发现，不论是健康还是腹泻的样本，菌门丰度排名第三的是疣微菌门。虽然因地域、饮食等差异，目前尚没有肠道微生物多样性水平与健康的具体标准，但高丰富度和均匀度的菌群结构更有利于肠道微生态紊乱后的恢复(Fassarellaet al.,2021)。

Costa和Weese(2018)的研究表明，变形杆菌门的大肠杆菌、沙门氏菌、嗜水气单胞菌和新立克次氏体等导致马腹泻，以上致病菌都属于革兰氏阴性细菌。也有研究表明，变形杆菌门丰度的增加是宿主代谢紊乱和免疫失调的标志(Shinet al.,2015)。由于高丰度的厚壁菌门和梭菌属对于维持马的健康很重要(Costaet al.,2012)，因此，厚壁菌门和梭菌属可能是马肠道微生物菌群的关键核心组成部分，其中只有少数梭菌被确定是肠道病原体。本研究显示，腹泻样本中的变形杆菌门、革兰氏阴性细菌和潜在致病性细菌丰度发生了显著性或趋向显著性上升，而厚壁菌门丰度下降。由此我们推测，厚壁菌门丰度下降、变形杆菌门丰度上升及属于革兰氏阴性细菌的潜在致病性细菌丰度的上升，都与腹泻相关。通常纤杆菌门的细菌可以分解纤维素以产生短链脂肪酸(Metzler and Mosenthin,2008)，腹泻致使纤杆菌门丰度下降可能直接影响宿主对多糖的利用。低丰度菌门在腹泻中的研究相对较少，有研究报道腹泻患者的浮霉菌门和迷踪菌门丰度显著高于健康个体(Drancourtet al.,2014;Xiet al.,2021)；蓝藻菌门、放线菌门丰度在腹泻患者中显著低于健康人群(Xiet al.,2021;Meiet al.,2021)。以上结果均与本研究结果一致。

广古菌门属于古菌，其中已知能够定植于动物体的有甲烷短杆菌属和甲烷球形菌属(Methanosphaera)。相比于健康组，腹泻组中丰度显著降低的大部分菌属在动物体内参与纤维素的代谢，如参与纤维素分解、果胶分解和半纤维素分解的有瘤胃球菌属、纤杆菌属、拟杆菌属、Paludibacter和梭菌属(Metzler and Mosenthin,2008;Panditet al.,2018)。可见腹泻通过影响消化道内纤维分解相关菌群的丰度，间接削弱了植食性宿主对多糖的利用。本研究显示，丰度发生显著变化的还有丁酸弧菌属，由于该菌属与脂肪酸合成相关，是潜在的益生菌属，在腹泻组中丰度增加是比较异常的现象，在旅行性腹泻(traveler’s diarrhea)的人群中也发生了类似的情况(Stampset al.,2020)。弯曲杆菌和沙门氏杆菌都属于变形杆菌门，都是导致人类腹泻的主要食源性病原体，在人和禽类中极为常见，但与我们的检测结果不同。可能因为马对弯曲杆菌反应不敏感，弯曲杆菌与马腹泻之间没有显著关联的报道(Netherwoodet al.,1996)也间接支持了我们的结果。

实时荧光定量PCR检测结果表明，各种菌属的差异情况与16S rRNA测序结果基本相符。据报道，普雷沃菌和拟杆菌的比值与宿主的纤维消耗、葡萄糖代谢呈正相关关系(Kovatcheva-Datcharyet al.,2015)，腹泻组此比值的下降虽未达到显著水平，但也可以提示腹泻通过影响消化道内纤维分解相关菌群，间接影响宿主对多糖的利用。本研究中，相比于健康组，腹泻组粪便中普雷沃菌含量低而双歧杆菌含量高，两者的含量可能具有负相关的趋势。虽然尚未有相关的报道，但对人类肠道微生物的研究发现，普雷沃菌可以抑制双歧杆菌效用(bifidogenic effect)(Chunget al.,2020)。

随着机器学习算法的不断优化，其在菌群数据分析中越来越多地发挥独特优势。其中，随机森林法被认为是高维数据领域中表现最好的分类器之一，近几年不断被用于肥胖、疾病预测等研究领域(Liuet al.,2019;Zenget al.,2019)。本研究利用随机森林分类方法筛选出对健康和腹泻柴达木马亚成体粪便微生物差异具有较大影响的12个潜在生物标记物，其中部分已被证实与肠道疾病相关，如甲烷短杆菌属中的M.smithii与肠易激综合征(IBS)相关(Stern and Brenner,2018)，迷踪菌属与人腹泻有关(Xiet al.,2021)；其他筛选出的如纤杆菌属、Paludibacter和肉食杆菌属与腹泻的关联性还有待进一步验证。由于生物标记物在药物开发、临床试验和治疗评估策略中有广泛的应用场景，本研究筛选出的生物标记物可用于诊断及疾病程度分类或疾病预后的标志(BDW Group,2001)。

综上所述，腹泻和健康状态的柴达木马亚成体粪便细菌群落结构差异极大，目前尚无法确定腹泻和肠道菌群多样性的因果关系，但是通过随机森林模型筛选出的生物标记物，具有成为治疗柴达木马腹泻靶点的潜力。本研究更倾向于对高原地区发生腹泻的亚成体柴达木马粪便微生物群落结构评估，以利于对青藏高原地区家畜尤其是马属动物的饲养。

致谢：在此特别感谢中国科学院微生物研究所东秀珠研究员对本研究提供的实验指导。