應用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)?？梢钥闯鱿柴R拉雅旱獺在腸道菌群的組成上與上述研究的哺乳動物腸道菌群存在一定的共性，但是也存在一定的區別，說明不同種類的動物，腸道菌群的組成比例完全不同，同一種細菌可能在一種動物腸道中占比較高，而在另外一種動物腸道中占比極少。

本文中喜馬拉雅旱獺腸道菌群分布以變形菌門的含量為最高，其次是厚壁菌門、擬桿菌門。變形菌門不僅是喜馬拉雅旱獺，也是大熊貓腸道中的主要菌群[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]，這無疑嚴重影響了喜馬拉雅旱獺的生存繁衍。但正因這些惡劣的環境給我們提供了獨特的研究機會，研究高原野生動物腸道微生物菌群多樣性組成，可為微生物學、生態學和獸醫學提供信息數據，進而對公共衛生事業的進步具有一定的促進作用。