腹瀉柴達木馬亞成體糞便微生物多樣性分析及其生物標記物的篩選

王小琪 郝文靜 張志超 韓晶,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認為兩組間具有顯著性差異趨向?；贠TU表格，使用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)?？梢姼篂a通過影響消化道內纖維分解相關菌群的豐度，間接削弱了植食性宿主對多糖的利用。本研究顯示，豐度發生顯著變化的還有丁酸弧菌屬，由于該菌屬與脂肪酸合成相關，是潛在的益生菌屬，在腹瀉組中豐度增加是比較異常的現象，在旅行性腹瀉(traveler’s diarrhea)的人群中也發生了類似的情況(Stampset al.,2020)。彎曲桿菌和沙門氏桿菌都屬于變形桿菌門，都是導致人類腹瀉的主要食源性病原體，在人和禽類中極為常見，但與我們的檢測結果不同?？赡芤驗轳R對彎曲桿菌反應不敏感，彎曲桿菌與馬腹瀉之間沒有顯著關聯的報道(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)。

綜上所述，腹瀉和健康狀態的柴達木馬亞成體糞便細菌群落結構差異極大，目前尚無法確定腹瀉和腸道菌群多樣性的因果關系，但是通過隨機森林模型篩選出的生物標記物，具有成為治療柴達木馬腹瀉靶點的潛力。本研究更傾向于對高原地區發生腹瀉的亞成體柴達木馬糞便微生物群落結構評估，以利于對青藏高原地區家畜尤其是馬屬動物的飼養。

致謝：在此特別感謝中國科學院微生物研究所東秀珠研究員對本研究提供的實驗指導。