徐 鵬，張雪飛，馬 英

青海是多種野生動物和珍稀保護(hù)動物的棲息地，生活著多種在中國乃至世界上都特有的物種，如喜馬拉雅旱獺(Himalayan marmot)、高原鼠兔(Ochotona curzoniae)、藏羚羊(Tibetan antelope)等。1954年青海省從喜馬拉雅旱獺體內(nèi)分離出鼠疫菌，被判定為喜馬拉雅旱獺鼠疫自然疫源地[1]。自此，喜馬拉雅旱獺就成為鼠疫監(jiān)測的重點對象，鼠防人員對其進(jìn)行研究，發(fā)現(xiàn)喜馬拉雅旱獺對鼠疫菌的留存、傳播和流行發(fā)揮重要作用。近年，國內(nèi)學(xué)者對青海部分地區(qū)的喜馬拉雅旱獺所攜帶的病毒方面開展初步研究，相繼在喜馬拉雅旱獺糞便中發(fā)現(xiàn)了一些新的病毒[2-3]，說明我們對于該物種的研究還有很多的空缺需要去補充。

隨著近年來國家加大對野生動物的保護(hù)力度，生態(tài)環(huán)境的日益改善，人們對于保護(hù)野生動物思想觀念的提升，多種因素使得野生動物的數(shù)量迅速增多。另一方面，城市化進(jìn)程的加劇，使人類活動和居住的范圍越來越大，動物的棲息地卻越來越小，野生動物不得不與人類一起共同生存，增加了人類與野生動物的接觸，從而為病毒、細(xì)菌的跨物種傳播提供了大量機會，通過野生動物的尿液、糞便或蜱類、螨類、蚤類等將貝納柯克斯體、漢坦病毒、鼠疫菌等傳播給人類[4]?，F(xiàn)已確認(rèn)的多種急性傳染病中，來源于野生動物的占比高達(dá)43%[5]。于是，野生動物間、野生動物與人之間的傳染病防治就成為了目前主要的防控要點。而腸道作為生物體內(nèi)一個龐大、復(fù)雜的微生物寄生區(qū)域，蘊藏著豐富的細(xì)菌和病毒，對其進(jìn)行細(xì)致深入的研究勢必成為防控傳染病的重要一環(huán)。

喜馬拉雅旱獺作為一種大型地棲類嚙齒動物，廣泛分布在青藏高原以及與中國接壤的尼泊爾等國的青藏高原邊緣山地。但目前的研究僅局限于傳統(tǒng)形態(tài)分類、物種生態(tài)調(diào)查和動物鼠疫監(jiān)測等方面，針對攜帶其它病原微生物的研究很少，腸道菌群微生物多樣性方面的研究目前還未見報道。因此，我們擬對青海不同地理區(qū)域的喜馬拉雅旱獺腸道菌群展開研究，為應(yīng)對、預(yù)防和控制未來可能發(fā)生的動物源性新發(fā)傳染病提供基礎(chǔ)數(shù)據(jù)，同時也為鼠疫的防控工作拓展新思路。

1 材料與方法

1.1 實驗材料

1.1.1 樣品采集 采樣點結(jié)合青海鼠疫流行重點樣區(qū)、旱獺種群遺傳特點和動物地理區(qū)劃，選擇青海省3個地理亞區(qū)的5個小區(qū)共采集45份青藏高原喜馬拉雅旱獺糞便樣本(使用一次性糞便采集器將新鮮糞便樣本挑取放入5 mL凍存管，并快速置于液氮罐中冷凍保存)。按照不同的采樣地區(qū)將45份樣本分為14組，采樣點詳細(xì)信息見表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，使殘留酒精完全揮發(fā)，保證后續(xù)實驗操作。

1.2 數(shù)據(jù)處理

1.2.1 數(shù)據(jù)質(zhì)控 測序后得到原始測序數(shù)據(jù)(Raw PE)，為了得到能夠用于后續(xù)分析的數(shù)據(jù)，首先去除干擾數(shù)據(jù)，通過PEAR(V0.9.6)將序列進(jìn)行拼接[9]，去除Barcode和引物序列得到Raw Tags[8]；依據(jù)Qiime(V1.7.0)對Raw Tags的質(zhì)控過程[10]，分別進(jìn)行Tags截取和Tags長度過濾等操作后得到Clean Tags；Clean Tags經(jīng)過UCHIME Algorithm(V4.2.40)去除嵌合體等步驟[11]，得到最終的有效數(shù)據(jù)(Effective Tags)[12]。

1.2.2 OTU聚類與物種注釋 利用Uparse軟件[13]對所有樣品的全部Effective Tags進(jìn)行聚類，以97%的一致性歸為一類，即一個可操作分類單元(Operational Taxonomic Units，OTUs)，其中無法被聚類的Tags序列和Tags序列頻數(shù)為1的序列將被舍棄，Tags 序列中出現(xiàn)頻數(shù)最高的將被作為 OTUs 的代表序列。使用Mothur(V1.30.1)[10]與SILVA的SSUrRNA數(shù)據(jù)庫對其進(jìn)行物種注釋，得到物種在域(domain)、門(phylum)、綱(class)、目(order)、科(family)、屬(genus)上的構(gòu)成。

1.2.3 微生物多樣性分析 基于各樣本數(shù)據(jù)的均一化處理，探究各個樣品的復(fù)雜度，采用Alpha多樣性分析。運用Miseq SOP[14]繪制等級豐富度圖(Rank Abundance)展現(xiàn)樣本的物種豐富度與均勻度、稀釋曲線反映測序深度是否合適；計算Chao 1與ACE指數(shù)[15]評估菌群豐富度、計算Shannon和Simpson指數(shù)[16]評估菌群多樣性。

1.2.4 物種分類與豐度聚類 為了得到每個OTU對應(yīng)的物種分類信息，需要對OTU進(jìn)行物種分類，我們采用RDP classifier方式。RDP classifier(V2.12)[6]基于Bergey’s taxonomy，采用Naive Bayesian assignment算法對每條序列在不同層級水平上計算其分配到此rank中的概率值。Bergey’s taxonomy分為6層，它們依次為domain、phylum、class、order、family、genus，根據(jù)其制作物種分類表。物種分類數(shù)據(jù)庫采用NCBI 16S(http://ncbi.nlm.nih.gov/)。根據(jù)分類學(xué)分析結(jié)果，統(tǒng)計在各個分類層級水平上的每個樣品的群落組成。

基于物種相對豐度的樣本聚類樹圖可以通過樹枝結(jié)構(gòu)直觀的反應(yīng)出多個樣品間的相似性和差異關(guān)系。首先使用R(V3.2)根據(jù)各樣本物種相對豐度計算beta多樣性距離矩陣，使用Bray Curtis法獲得樣本間距離后進(jìn)行層次聚類(Hierarchical clustering)分析，再使用非加權(quán)組平均法UPGMA(Unweighted pair group method with arithmetic mean)算法構(gòu)建樹狀結(jié)構(gòu)，得到樹狀關(guān)系形式用于可視化分析。

1.2.5 腸道菌群與環(huán)境因子的關(guān)系 根據(jù)不同的采樣地區(qū)將45份樣本分為14組，在R上先使用species-sample數(shù)據(jù)(97%相似性的樣品OTU表)做DCA(Detrended Correspondence Analysis)分析，根據(jù)分析結(jié)果中Lengths of gradient的第一軸的大小選擇冗余分析(Redundancy Analysis，RDA)或典范對應(yīng)分析(Canonical Correspondence Analysis，CCA)。RDA或者CCA是基于對應(yīng)分析發(fā)展而來的一種排序方法，將對應(yīng)分析與多元回歸分析相結(jié)合，每一步計算均與環(huán)境因子回歸，主要是用來反映每一組樣本的腸道菌群與環(huán)境因子(溫度、海拔、經(jīng)度、緯度)之間的關(guān)系。

為探究環(huán)境因子具體對腸道哪一種菌群產(chǎn)生影響，將腸道菌群在屬水平含量上從高到低排序，通過使用R的pearson correlation算法計算相關(guān)性系數(shù)(coefficient)來分析前20個含量較高的細(xì)菌與海拔、溫度、經(jīng)度、緯度的相關(guān)性，最終獲得顯著影響某些物種的環(huán)境因子。其中P值越小相關(guān)性越大，表明該細(xì)菌屬與該環(huán)境因子的關(guān)系較為密切。

1.2.6 菌群相對豐度差異及功能預(yù)測分析 基于物種分類結(jié)果，計算在不同水平上各rank的豐度，比較樣本間相對豐度差異，找出由于生存環(huán)境不同而存在顯著差異的物種分類，關(guān)注樣本中存在明顯差異的菌屬。篩選條件為P≤0.05，比較對象為樣本間兩兩比較，采用Fisher’s精確檢驗(Fisher exact test)，最后將檢驗得到的P值采用FDR(False Discovery Rate)做Multiple test correction得到q值。

利用PICRUSt軟件(V1.0.0)[17]對喜馬拉雅旱獺腸道菌群的功能進(jìn)行預(yù)測分析，基于COG基因預(yù)測結(jié)果對所有樣品進(jìn)行聚類與柱狀圖組合分析，樣本聚類樹圖可以通過樹枝結(jié)構(gòu)直觀的反應(yīng)出多個樣品間的相似性和差異關(guān)系。

2 結(jié) 果

2.1 Alpha多樣性分析 通過數(shù)據(jù)拼接、質(zhì)控和嵌合體過濾等步驟，將45份喜馬拉雅旱獺糞便樣本進(jìn)行16S rDNA基因測序分析，所有樣本菌群豐富度與菌群多樣性均符合實驗要求，覆蓋率除XG3-18為96.55%、GD1-18為97.78%、WT1-18為97.79%以外，其余均高于98%。Rank Abundance曲線較寬且平坦，說明物種組成豐富且均勻度高；稀疏曲線趨于平穩(wěn)且其中42份樣本覆蓋率均高于98%，說明樣本測序深度足夠[18]。

2.2 腸道菌群物種分類 從門分類水平上分析喜馬拉雅旱獺的腸道菌群，排名前6位的分別是：變形菌門(Proteobacteria)、厚壁菌門(Firmicutes)、擬桿菌門(Bacteroidetes)、放線菌門(Actinobacteria)、疣微菌門(Verrucomicrobia)、軟壁菌門(Tenericutes)。其中變形菌門占比最高，為38.77%；軟壁菌門占比最低，為0.15%。

在屬水平上排名前10位的包括：不動桿菌屬(Acinetobacter)、節(jié)桿菌屬(Arthrobacter)、鞘氨醇桿菌屬(Sphingobacterium)、埃希氏菌屬(Escherichia)、另枝菌屬(Alistipes)、黃桿菌屬(Flavobacterium)、瘤胃球菌屬(Ruminococcus)、肉桿菌屬(Carnobacterium)、梭菌屬 IV(ClostridiumIV)、叢毛單胞菌屬(Comamonas)，其中不動桿菌屬占比最高，為20.36%，是樣本W(wǎng)T2-18、CD11-18、GH11-18、WK8-18中的優(yōu)勢菌屬；叢毛單胞菌屬占比最低，為1.94%，見表2。

表2 樣品分類系統(tǒng)組成表Tab.2 Composition of all sample classification systems at the genus level

2.3 物種相對豐度聚類 根據(jù)各樣本物種相對豐度計算beta多樣性距離矩陣，使用Bray Curtis法獲得樣本間距離后進(jìn)行層次聚類得到3個聚類，第一聚類又分為兩個類別，第一類別不動桿菌屬(Acinetobacter)較高，不動桿菌屬為革蘭氏陰性專性需氧菌，第二類別鞘氨醇桿菌屬(Sphingobacterium)較高，其屬于好氧或兼性厭氧非發(fā)酵革蘭氏陰性桿菌；第二聚類中未分類菌(unclassified)含量較高，并且是在屬分類狀況下，表明其中含有大量新奇且未被充分研究的新微生物；第三聚類由于樣本太少(2份)，故不做分析，結(jié)果如圖1所示。

圖分為左、中、右3個部分，圖左：基于Bray-Curtis的樣本聚類樹圖；圖中：根據(jù)聚類順序排序的物種豐度柱狀圖；圖右：物種的圖示。圖1 基于物種相對豐度的所有樣品聚類樹圖與柱狀圖組合分析圖Fig.1 Combination analysis graph of clustering dendrogram and histogram for all samples based on relative abundance of species

2.4 樣本與環(huán)境因素的RDA分析 根據(jù)DCA的分析結(jié)果，Lengths of gradient的第一軸<3，于是對收集到的樣本環(huán)境因子信息與樣本采用RDA分析，可得：溫度(T)≈海拔(ASL)>經(jīng)度(E)>緯度(N)，如圖2所示。T對1組與2組影響較大。ASL對8組、3組、7組、9組、10組、6組影響較大。E對12組、4組、1組和2組影響較大，N對各組的影響較小。海拔與溫度對腸道菌群的影響最大，海拔、溫度不同，腸道菌群的組成也不同[19-22]。通過細(xì)菌屬與環(huán)境因子的相關(guān)性分析后發(fā)現(xiàn)不動桿菌屬(Acinetobacter)與海拔、溫度、經(jīng)度關(guān)系密切；瘤胃球菌屬(Ruminococcus)與溫度、經(jīng)度相關(guān)；馬賽菌屬(Massilia)與海拔、緯度相關(guān)；肉桿菌屬(Carnobacterium)與乳桿菌屬(Lactobacillus)僅與溫度相關(guān)；假單胞菌屬(Pseudomonas)僅與海拔相關(guān)。

圖2 14組樣本的RDA圖Fig.2 RDA chart of 14 sets of samples

2.5 樣本間菌群相對豐度差異及功能預(yù)測分析 探究單個樣本由于生存環(huán)境不同，腸道菌群中哪些細(xì)菌屬存在差異，使用Fisher exact test檢驗后發(fā)現(xiàn)45個樣本中有14個屬的細(xì)菌存在統(tǒng)計學(xué)差異(表3)。

表3 45個樣本中存在明顯差異的細(xì)菌Tab.3 Bacteria with obvious differences in 45 samples

為探究樣本細(xì)菌功能，利用PICRUSt軟件，基于COG的功能結(jié)構(gòu)對喜馬拉雅旱獺腸道菌群的功能進(jìn)行預(yù)測分析，結(jié)果顯示喜馬拉雅旱獺腸道菌群功能預(yù)測大多集中在遺傳信息處理與代謝通路上，如基因一般功能預(yù)測，未知功能轉(zhuǎn)錄、復(fù)制、重組與修復(fù)，核糖體結(jié)構(gòu)翻譯與生物發(fā)生，氨基酸運輸與代謝，能源生產(chǎn)與轉(zhuǎn)換，細(xì)胞壁與細(xì)胞膜的生物發(fā)生。表明同一物種，不同生存環(huán)境，其腸道微生物群落的功能存在一定的差異。其中有6個樣本(XH16-18、XH12-18、XH15-18、GH11-18、GM6-18、TR89-18)整體功能豐度較低且結(jié)構(gòu)較為相似，且對照采樣點信息，這6個樣本都不處于喜馬拉雅旱獺的最適生存環(huán)境[23]；整體功能豐度最高的3組樣本(XG3-18、DH3-18、XB3-18)的生存環(huán)境均處于青藏高原喜馬拉雅旱獺的最適生存環(huán)境，結(jié)果如圖3所示。

圖分為左、中、右三個部分，圖左：基于Bray-Curtis的樣本聚類樹圖；圖中：根據(jù)聚類順序排序的功能豐度柱狀圖；圖右：物種的圖示。 圖3 基于COG的所有樣品聚類樹與柱狀圖組合分析圖Fig.3 Combination analysis graph of the clustering tree and histogram based on COG for all samples

3 討 論

喜馬拉雅旱獺樣本中分類操作單元(OTU)數(shù)目越多，表示細(xì)菌分類越多樣，反之則分類越單一。本研究中第一聚類的樣本平均OTU要遠(yuǎn)小于第二聚類，提示第一聚類的樣本所面臨的環(huán)境復(fù)雜程度要小于第二聚類，回顧采樣點的生境條件，第二聚類樣本的環(huán)境復(fù)雜度確實要高于第一聚類的樣本，但之前沒有人進(jìn)行過相關(guān)的實驗研究，所以我們目前還只是推測，具體的驗證工作有待于開展下一步的實驗。

Ley等[24]利用宏基因組技術(shù)對來自13個生物分類的60種哺乳動物的腸道菌群分析后發(fā)現(xiàn)相對豐度排序前6的門為：厚壁菌門(Firmicutes)、擬桿菌門(Bacteroidetes)、變形菌門(Proteobacteria)、放線菌門(Actinobacteria)、疣微菌門(Verrucomicrobia)、梭桿菌門(Fusobacteria)；馮天舒等[25]對高原鼢鼠的盲腸內(nèi)容物進(jìn)行16SrRNA基因測序后獲得該物種的腸道菌群組成，其主要的腸道微生物有：擬桿菌門(Bacteroidetes)、厚壁菌門(Firmicutes)、變形菌門(Proteobacteria)、螺旋體門(Spirochaetes)、放線菌門(Actinobacteria)；譚春桃等[26]對青藏高原夏季野生鼠兔的腸道菌群進(jìn)行研究后發(fā)現(xiàn)其腸道菌群主要集中在擬桿菌門(Bacteroidetes)、變形菌門(Proteobacteria)、厚壁菌門(Firmicutes)、放線菌門(Actinobacteria)、藍(lán)藻菌門(Cyanobacteria)、螺旋體門(Spirochetes)?？梢钥闯鱿柴R拉雅旱獺在腸道菌群的組成上與上述研究的哺乳動物腸道菌群存在一定的共性，但是也存在一定的區(qū)別，說明不同種類的動物，腸道菌群的組成比例完全不同，同一種細(xì)菌可能在一種動物腸道中占比較高，而在另外一種動物腸道中占比極少。

本文中喜馬拉雅旱獺腸道菌群分布以變形菌門的含量為最高，其次是厚壁菌門、擬桿菌門。變形菌門不僅是喜馬拉雅旱獺，也是大熊貓腸道中的主要菌群[27]。喜馬拉雅旱獺作為植食性動物，禾木科和莎草科植物的綠色部分為主要飼料，間或食種子及花序，變形菌門細(xì)菌通過降解動物食物中的木質(zhì)素，破壞植物細(xì)胞壁，使動物快速汲取營養(yǎng)物質(zhì)；厚壁菌門作為纖維素分解細(xì)菌，可將植物纖維降解為揮發(fā)性脂肪酸以供宿主利用；擬桿菌門則主要幫助宿主降解碳水化合物、蛋白質(zhì)，提高宿主營養(yǎng)利用率，維持腸道的微生物生態(tài)平衡，促進(jìn)宿主的免疫系統(tǒng)發(fā)育，增強宿主免疫水平[28-30]。各種菌屬協(xié)調(diào)共生，相輔相成，共同提高宿主生命健康水平。

有研究表明青藏高原喜馬拉雅旱獺的最適生長溫度為14.3～24 ℃，最適生長高度為3 000～4 000 m[23]，基于Bray-Curtis建立的樣本聚類樹圖中，第一聚類第一分類中不動桿菌屬含量較高，這是一種專性需氧菌，高寒低氧環(huán)境可使其活力下降。第二分類中鞘氨醇桿菌屬含量較高，該種細(xì)菌屬含有大量的鞘磷脂[31]，是細(xì)胞膜的主要組成成分。綜合采樣信息，第一聚類中有大量樣本不處于喜馬拉雅旱獺最適生存環(huán)境，為了適應(yīng)高海拔給機體帶來的影響，不動桿菌屬和鞘氨醇桿菌屬以相對較高的含量來應(yīng)對高海拔所造成的低氧環(huán)境；第二聚類的喜馬拉雅旱獺腸道菌群中未分類菌含量較多，究其原因是之前很少全面研究高寒低氧環(huán)境下野生動物的腸道微生物菌群，我們猜測在這一聚類中遇到了某些較為新奇、尚未被充分研究的微生物，這些微生物可能在宿主處于極端環(huán)境下發(fā)揮作用，但實際情況是否如此，有待于后續(xù)進(jìn)一步的深入研究。

對樣本豐度存在顯著差異的物種進(jìn)行t檢驗，得到14種有統(tǒng)計學(xué)意義的細(xì)菌屬別，我們猜測這些差異可能與樣本的腸道菌群的功能分布有關(guān)。因為同一物種處于不同環(huán)境之下，為了促進(jìn)宿主適應(yīng)環(huán)境，腸道菌群也隨之發(fā)生了不同程度的結(jié)構(gòu)、組成改變，從而腸道菌群的功能分布也發(fā)生變化，以達(dá)到一種自然狀態(tài)下的保護(hù)性機制[32]。并且考慮采樣時間，該階段喜馬拉雅旱獺處于交配期，存在交配競爭、食物競爭、領(lǐng)地競爭等現(xiàn)象。為了在競爭中獲得有利地位，必然會通過多種途徑提升自己的體重，攝入多元充足的食物，使身體各項機能指標(biāo)達(dá)到最佳狀態(tài)[33]，從而導(dǎo)致其相應(yīng)腸道菌群的功能豐度上升。

PICRUSt分析發(fā)現(xiàn)喜馬拉雅旱獺的腸道菌群功能預(yù)測多集中在遺傳信息處理與代謝通路上，說明低氧可能對宿主腸道微生物的基因造成了損傷，這與馬燕等[34]的研究結(jié)果一致。高寒缺氧帶給喜馬拉雅旱獺的不僅僅是生存壓力，其繁殖壓力也隨之而來[35-36]。本研究有64.45%的樣本含有不動桿菌，該種細(xì)菌為機會致病菌，在機體免疫力低下時導(dǎo)致機體疾病的發(fā)生，加之暴露于高寒低氧條件本身就會導(dǎo)致動物疾病的發(fā)生機會增加[37-38]，這無疑嚴(yán)重影響了喜馬拉雅旱獺的生存繁衍。但正因這些惡劣的環(huán)境給我們提供了獨特的研究機會，研究高原野生動物腸道微生物菌群多樣性組成，可為微生物學(xué)、生態(tài)學(xué)和獸醫(yī)學(xué)提供信息數(shù)據(jù)，進(jìn)而對公共衛(wèi)生事業(yè)的進(jìn)步具有一定的促進(jìn)作用。