峨眉山與黃山藏酋猴腸道菌群組成的比較

翟子豪， 宋飏， 王俊茵， 張可俊， 孫丙華， 李靜*

(1. 四川大學生命科學學院，生物資源與生態環境教育部重點實驗室，四川省瀕危野生動物保護生物學重點實驗室，成都610065； 2. 峨眉山景區管委會，峨眉山生物多樣性保護研究所，四川峨眉山614200；3.安徽大學資源與環境工程學院，安徽省黃山生物多樣性與短尾猴行為生態學國際聯合研究中心，合肥230039)

哺乳動物通過母體產道出生時，首次獲得腸道微生物(Vaishampayanetal.，2010)。隨著時間的推移，在食物、社會接觸或周圍環境的影響下，腸道微生物多樣性逐漸增加(Davidetal.，2014)。這些微生物之間及微生物與宿主之間形成了相互依存、相互作用的不可分割的整體。微生物菌群和宿主相互交換能量物質、傳遞信息，在宿主的營養、免疫、代謝中起重要作用(Egertetal.，2006)。隨著食物、生境等的改變，腸道微生物的組成也會發生相應的變化。研究表明，在不同季節(Sunetal.，2016)、圈養與野生(Wangetal.，2016)、不同海拔(Zhaoetal.，2018)的條件下，同種動物的腸道微生物組成均有明顯的差異。Clayton等(2016)發現，野生白臀葉猴Pygathrixnemaeus和鬃毛吼猴Alouattapalliata具有明顯不同的腸道菌群，而在人工圈養條件下卻表現出了較大的相似性，在人體中占主導地位的普氏菌屬Prevotella和擬桿菌屬Bacteroides物種均增加。Zhao等(2018)發現高海拔地區的恒河猴Macacamulatta較低海拔地區的腸道微生物多樣性更高，且有多種獨有的可操作分類單元(operational taxonomic units，OTUs)。

藏酋猴Macacathibetana，又稱毛面短尾猴、大青猴，隸屬于靈長目Primates猴科Cercopithecidae獼猴屬Macaca，是我國特有的靈長類，國家Ⅱ級重點保護野生動物，世界自然保護聯盟(IUCN)瀕危物種紅色名錄近危(NT)物種(Long & Richardson，2008)，《中國脊椎動物紅色名錄》易危(VU)物種(蔣志剛等，2016)。由于適宜棲息地不斷縮小，藏酋猴現主要分布在四川、安徽、貴州、福建等地，分為4個亞種：指名亞種M.t.thibetana、貴州亞種M.t.guizhouensis、黃山亞種M.t.huangshanensis和福建亞種M.t.pullus(蔣學龍等，1996)。

峨眉山和黃山分別位于我國西南部和中東部地區，不同的地理位置塑造了兩地獨有的自然環境。黃山植被主要為次生常綠闊葉林，峨眉山植被以亞熱帶常綠闊葉林為主。同時，兩地都以其得天獨厚的藏酋猴資源開展了靈長類生態旅游，并且兩地管理部門都對游客向藏酋猴的投喂行為做出了限制(黃山不允許游客投喂，由管理人員定時、定點、定量投食；峨眉山允許游客購買特定的猴糧投喂)。目前，關于黃山、峨眉山藏酋猴的研究主要集中在棲息活動范圍(Zhao，1994)、形態和社群行為(Zhao，1997)等方面，發現兩地藏酋猴的猴群大小差異較大(孫丙華等，2010)；Sun等(2016)報道了不同季節的黃山藏酋猴腸道微生物在組成和多樣性上均有顯著差異。但兩地猴群腸道微生物組成上的差異迄今尚無相關報道。

本文以16SrRNA基因的V3-V4高變區為分子標記，使用高通量測序技術研究了峨眉山藏酋猴群體的腸道微生物，并與黃山藏酋猴的數據進行比較分析，為揭示藏酋猴腸道微生物的特征提供資料，也為藏酋猴的保護與管理提供科學依據。

1 材料與方法

1.1 研究地點和樣本信息

峨眉山位于四川省峨眉山市境內，地理位置103°10′～103°37′E，29°16′～29°43′N，最高海拔3 099 m；黃山位于安徽省南部黃山市境內，地理位置118°01′～118°17′E，30°01′～30°18′N，最高海拔1 864.8 m，兩地相距超1 460 km。峨眉山藏酋猴(EM)31個糞便樣本采集于生態猴區和雷洞坪，采樣時間為5月。在猴群每天出現的游步道守候，等待猴群進食完畢，采集新鮮糞便，回到住所后立即置于液氮中保存。為確保所有樣本均來自不同個體，樣本間距需大于1.5 m。此外，通過本實驗室篩選的微衛星位點進一步確認，共有25只不同個體的樣品。黃山藏酋猴(HS)樣本下載于NCBI(GenBank 登錄號：SRP073002)，共25個樣品，為黃山野生猴谷的魚鱗坑YA1群體，采樣時間為春季，測序平臺為Illumina MiSeq PE250，共產生290 542條高質量序列，平均長度為415.28 bp±9.81 bp(Sunetal.，2016)。

1.2 DNA提取、PCR擴增

取糞便內部部分，使用MoBio PowerSoil?DNA Isolation Kit(12888)提取DNA。16SrRNAV3-V4區的PCR擴增引物為338F(5’-ACTCCTACGGGAGGCAGCA-3’)和806R(5’-GGACTACHVGGGTWTC-TAAT-3’)。反應體系25 μL：Q5 high-fidelity DNA polymerase 0.25 μL，5×Reaction Buffer 5 μL，5×High GC Buffer 5 μL，dNTP(1×104μmol)0.5 μL，模板DNA 1 μL，正、反向引物(10 μmol)各1 μL， ddH2O 11.25 μL。PCR程序：98 ℃ 2 min；98 ℃ 15 s，55 ℃ 30 s，72 ℃ 30 s，27個循環；72 ℃ 5 min。擴增產物使用2%瓊脂糖凝膠電泳檢測，使用Axygen凝膠回收試劑盒回收目的片段。采用Illumina MiSeq PE250平臺對糞便微生物群落DNA片段進行雙向測序，原始測序數據與下載的黃山數據合并后進行后續分析。

1.3 數據處理

使用Perl語言編寫腳本，采用滑動窗口法對FASTQ格式的雙端序列作質量控制：窗口大小為10 bp，步長為1 bp，從5’端第一個堿基位置開始移動，要求窗口中堿基平均質量≥Q20，從第一個平均質量值低于Q20的窗口處截斷序列，且截斷后的序列長度≥150 bp，不允許存在模糊堿基N。使用FLASH(Magoc & Salzberg，2011)對通過質量初篩的雙端序列進行拼接，之后根據每個樣本所對應的Barcode序列，使用QIIME 1.9.1(Caporasoetal.，2010)將拼接后的序列按樣本分開。擴增子測序結果中序列重復度高，并且大量出現1次或幾次的序列，統計學和功能上意義不大。因此，使用Usearch(Edgar，2010)去冗余并過濾低豐度序列，得到非冗余序列后，調用silva.gold數據庫(Quastetal.，2013)去除嵌合體序列，之后選用Uparse算法(Edgar，2013)按97%的序列相似度進行聚類和OTUs劃分，并生成代表性序列。隨后，與Greengenes數據庫(Release 13.8)(McDonaldetal.，2012)比對，獲取每個OTU所對應的分類學信息，并生成具有注釋信息的OTUs表。

1.4 統計分析

去除豐度值低于全體樣本測序總量0.01%的OTUs(Bokulichetal.，2013)，通過mothur(Schlossetal.，2009)繪制稀疏曲線。使用Clustal Omega(Sieversetal.，2011)對代表序列進行多序列比對，數據過濾后，使用FastTree 2.1.9(Priceetal.，2009)構建系統發育樹。對OTUs進行拉平處理，之后計算每個樣本的α多樣性和β多樣性。α多樣性由反映群落豐富度的ACE指數、Chao1指數，反映群落多樣性的Simpson指數、Shannon指數和反映樣品覆蓋度的Good’s coverage指數進行評估，并通過箱線圖進行可視化。β多樣性選擇Unweighted UniFrac距離(Lozuponeetal.，2011)對樣本進行比較，通過主坐標分析(principal co-ordinates analysis，PCoA)、主成分分析(principal component analysis，PCA)和樣本的UPGMA聚類分析進行可視化表達。使用LEfSe(Segataetal.，2011)進行組間豐度差異的可視化，首先采用非參數因子Kruskal-Wallis秩和檢驗檢測組間豐度差異顯著的物種，之后對這些物種進行成對Wilcoxon秩和檢驗，最后通過線性判別分析(LDA)對數據進行降維和評估差異顯著的物種的影響力。采用PICRUSt(Langilleetal.，2013)對序列進行基因功能預測分析，并使用LEfSe進行組間差異可視化。本研究中的顯著性水平檢驗均選用單因素方差分析(One-Way ANOVA，α=0.05)。

2 結果

2.1 測序數據

經質量控制后，50個樣品共得到高質量序列714 259條[(14 285條±3 274條)/樣品]，其中，EM為 423 717條，HS為290 542條，每條序列平均長度為413.87 bp±12.34 bp。序列按97%相似度聚類后得到500個OTUs，其中，EM獨有119個，HS獨有34個，兩地共有347個(分別占峨眉山和黃山的74%和91%)。按樣品最低序列數(4 203條)進行拉平處理。Good’s coverage指數平均值為99.15%±0.14%(98.76%～99.46%)，表明99%的樣品均得到了鑒定。50個樣品的稀疏曲線趨于平緩(圖1)，表明測序數據量已覆蓋樣本中的絕大多數物種。

圖1 峨眉山藏酋猴和黃山藏酋猴各樣品稀疏曲線
Fig. 1 The rarefaction curves of each sample ofMacacathibetanabetween Mount Emei and Mount Huangshan

2.2 EM與HS腸道微生物組成

Greengenes數據庫(Release 13.8)匹配到的OTUs可劃分為9門46科79屬。其中，EM腸道中，相對豐度大于1%的門有8個：厚壁菌門Firmicutes(69.04%±11.81%)、擬桿菌門Bacteroidetes(21.59%±10.05%)、放線菌門Actinobacteria(2.73%±2.17%)、螺旋體門Spirochaetes(1.69%±2.74%)、藍藻菌門Cyanobacteria(1.49%±2.60%)、變形菌門Proteobacteria(1.27%±1.75%)、軟壁菌門Tenericutes(1.16%±1.75%)、疣微菌門Verrucomicrobia(1.00%±1.72%)，占門水平總豐度的99.97%；HS腸道中，相對豐度大于1%的門有3個：厚壁菌門(46.34%±8.15%)、擬桿菌門(36.75%±6.38%)、變形菌門(14.91%±8.06%)，占門水平總豐度的97.97%。EM腸道中，厚壁菌門與擬桿菌門相對豐度的比值(3.20)顯著高于HS(1.26，P<0.05)。在屬級水平上，EM腸道中相對豐度最高的為顫螺菌屬Oscillospira(23.49%±16.63%)，其次為普氏菌屬(19.94%±14.48%)、柔嫩梭菌屬Faecalibacterium(18.40%±12.41%)；HS腸道中相對豐度最高的為普氏菌屬(36.35%±9.15%)，其次為琥珀酸弧菌屬Succinivibrio(13.94%±11.47%)、顫螺菌屬(13.76%±10.59%)(圖2)。

為了比較兩地藏酋猴腸道微生物的差異，選用LEfSe發掘在豐度上有顯著差異的類群。在EM中，厚壁菌門、梭菌綱Clostridia、梭菌目Clostrisiales、瘤胃菌科Ruminococcaceae和毛螺菌科Lachnospiraceae均顯著富集，主要表現在布勞特氏菌屬Blautia、糞球菌屬Coprococcus、顫螺菌屬、瘤胃球菌屬Ruminococcus等的增加，同時螺旋體門的螺旋體綱Spirochaetia、螺旋體目Spirochaetales、螺旋體科Spirochaetaceae、密螺旋體屬Treponema也在EM中顯著增加(|LDA score|>4且P<0.05)(圖3)。而在HS中，擬桿菌門的擬桿菌綱Bacteroidia、擬桿菌目Bacteroidales、普雷沃氏菌科Prevotellaceae、普氏菌屬，變形菌門的胃螺桿菌屬Flexispira和琥珀酸弧菌屬均顯著富集。

圖2 峨眉山與黃山藏酋猴腸道菌群中門(上)和屬(下)的相對豐度Fig. 2 Relative abundance of gut microbiome taxa on phylum (upper) and genus (bottom) category in Macaca thibetana from Mount Emei and Mount Huangshan

圖3 LEfSe分析峨眉山和黃山藏酋猴腸道微生物組成的差異Fig. 3 Analysis of gut microbial composition in Macaca thibetana from Mount Emei and Mount Huangshan based on LEfSe

A. LDA值分布柱狀圖， B. 進化分支圖： 由內到外輻射的圓圈代表由門至屬的分類級別， 圓圈的大小代表相對豐度的大小， 節點代表起重要作用的微生物類群， 黃色代表差異無統計學意義的類群

A. LDA value distribution histogram， B. evolutionary branching diagram： the circles from the inside to the outside represent the classification level from the phylum to the genus， the size of the circle represents the relative abundance level， the node represents the important microbial group， and the yellow node represents microbiomes of no significant differences

2.3 EM與HS腸道微生物多樣性

EM的ACE指數(P<0.05)、Shannon指數(P<0.01)均顯著高于HS，而Chao1指數(P>0.05)和Simpson指數(P>0.1)雖然在EM中較高，但在組間的差異無統計學意義(圖4)。PCA、PCoA分析結果均表明，EM與HS在第一軸均明顯分開(圖5)。

2.4 EM與HS腸道微生物基因功能分析

PICRUSt對EM和HS微生物群落在KEGG代謝途徑的差異分析結果顯示，共匹配到41個KEGG二級代謝通路，其中差異有統計學意義的代謝通路15個(|LDA score|>2，P<0.05)。EM顯著富集的通路5個：脂代謝、外源化學物的生物降解與代謝、信號轉導、轉錄、膜運輸；HS顯著富集的通路10個：多糖生物合成與代謝、核苷酸代謝、輔因子及維生素代謝、能量代謝、細胞通路和信號、細胞生長與死亡、萜類化合物及多肽代謝、氨基酸代謝、運輸和分解代謝、信號分子與互作。

圖4 峨眉山藏酋猴和黃山藏酋猴腸道微生物α多樣性指數Fig. 4 α diversity index of gut microbiome in Macaca thibetana from Mount Emei and Mount Huangshan

圖5 基于PCA和PCoA分析峨眉山和黃山藏酋猴腸道菌群的差異Fig. 5 Difference in the gut microbiome of Macaca thibetana between Mount Emei and Mount Huangshan based on PCA and PCoA

3 討論

宿主的生活環境、飲食組成決定了其腸道微生物的組成，而腸道微生物組成的變化也反過來影響宿主的健康、食物消化、營養獲取。對EM和HS腸道菌群的比較研究有利于了解動物對不同生境的適應性進化，也有利于評估不同生態旅游管理模式對藏酋猴的潛在影響。本研究發現，EM和HS的腸道微生物組成以厚壁菌門、擬桿菌門、變形菌門、放線菌門為主，符合脊椎動物腸道微生物組成的共同特征(Deng & Swanson，2015)，也與其他非人靈長類動物相似，如恒河猴(Yasudaetal.，2015)。EM和HS共有的OTUs數量分別占其總數量的74%和91%，表明藏酋猴腸道微生物群落結構相對穩定。

EM和HS腸道微生物在組成、豐度和多樣性方面存在明顯的差異。峨眉山和黃山地理位置的差異造就了兩地不同的植被分布，而這也影響著藏酋猴食物的組成。藏酋猴食性復雜，通常以枝葉、果實、樹根等為主，也取食某些昆蟲(Bermanetal.，2007)。尤碩愚(2013)報道HS取食共計26科50種植物，主要包括禾本科Poaceae、殼斗科Fagaceae、樟科Lauraceae、杜鵑花科Ericaceae、金縷梅科Hamamelidaceae、山茶科Theaceae、豆科Leguminosae等，其中以殼斗科、樟科的取食量尤為突出，占全年總取食量的51.26%～59.75%。春季時，它們還偏愛取食竹筍，在食物中占據一定比例(熊成培，1984)。而EM取食的植物則高達177種，主要包括五加科Araliaceae、殼斗科、禾本科、木通科Lardizabalaceae、蕁麻科Urticaceae等，此外，EM還被觀察到取食地面小型無脊椎動物(Zhaoetal.，1991)。食物組成上的差異可能是兩地藏酋猴腸道菌群組成和多樣性產生差異的主要原因。無論是評價豐富度的ACE指數還是評價多樣性的Shannon指數，EM均顯著高于HS，這可能與峨眉山豐富的食物種類密切相關。而HS腸道中具有較高豐度的普氏菌屬，它們在分解利用水果、谷類和嫩葉中的半纖維素、果膠、淀粉、單糖等時發揮著重要作用(Russell & Baldwin，1979)，這可能與HS攝食更多的谷物、竹筍等有關。而EM腸道中富集的顫螺菌屬和瘤胃球菌屬則可能與食物中的高比例纖維含量相關(Jindouetal.，2006)。

峨眉山和黃山在藏酋猴生態旅游管理方式上的顯著不同也是導致兩地藏酋猴腸道微生物組成差異的重要原因之一。黃山為吸引藏酋猴到固定地方供游客觀賞，管理人員會每天定點定量喂食6 kg玉米(Mathesonetal.，2006)，人工食物已成為藏酋猴的固定食物組成；峨眉山則允許游客購買指定的猴糧(主要為花生和紅薯干)自行投喂，因此，藏酋猴獲取人工食物的量會隨著季節和人流量而波動，一般夏秋季居多，冬季最少(本實驗調查)。長期固定的人工飲食模式可能導致動物腸道菌群α多樣性降低(Uenishietal.，2007)，這可能也是HS腸道菌群多樣性較EM低的原因。此外，EM具有較高豐度的厚壁菌門與較低豐度的擬桿菌門。厚壁菌門微生物攜帶許多編碼能量代謝相關酶的基因，并且可以產生許多消化酶來分解各種物質，從而幫助宿主消化和吸收養分(Nuriel-Ohayonetal.，2016)。擬桿菌門則在蛋白質、碳水化合物(特別是多糖)等其他能量物質吸收方面發揮著重要作用，從而增加宿主中養分的利用率(Backhedetal.，2004)。研究表明，在經常食用高脂飲食的西方人中，會出現厚壁菌門豐度增加、擬桿菌門豐度減少的現象(Filippoetal.，2010)。而高厚壁菌門、低擬桿菌門也是肥胖人群腸道菌群的重要表現(Clarkeetal.，2012)。腸道菌群KEGG功能分析也發現，EM在與脂類代謝有關的功能通路上有更多的富集，這可能是由于峨眉山游客投喂的猴糧為花生等脂類含量較高的食物，脂類攝入量的增加導致EM腸道中厚壁菌門類群增多，并加強了其對脂類食物的代謝。此外，野外觀察也發現部分體型較胖的個體，但這是否與游客投喂食物有關尚需進一步研究。

值得注意的是，在EM中發現了一定豐度的螺旋體門和密螺旋體屬，其通常被認為是潛在的病原體(Schwan，1996)，可能導致各種慢性傳染性皮膚病，如雅司和品他病(Karlssonetal.，2013)。作為靈長類生態旅游的重要組成部分，近距離觀賞藏酋猴越來越受到追捧。在管理上，黃山實施人猴隔離，禁止人猴接觸及游客投喂；而EM則是逗留在游步道旁，與游客直接接觸，給疾病傳播提供了條件，可能給靈長類動物和游客都帶來潛在威脅(Pedersen & Davies，2009)。KEGG功能分析也發現，EM腸道菌群在外源化學物的代謝與降解途徑有顯著的功能富集。外源化學物主要包括一些人工合成物或由外界環境攝入，而非機體代謝產生的內源物質。而游客與藏酋猴的直接接觸或投喂不適當的食物都可能是導致它們在該代謝通路活動增強的原因。

通過比較和分析峨眉山和黃山藏酋猴的腸道菌群組成，發現兩地的藏酋猴在厚壁菌門、擬桿菌門等腸道微生物組成、豐度，以及多樣性上存在顯著差異，這可能與兩地藏酋猴在食物組成、生態旅游的管理模式上的差異密切相關。本研究為藏酋猴的保護和合理開發生態旅游提供了科學參考。

致謝：感謝峨眉山景區管委會郝國歉在采樣過程中給予的支持與幫助，謹此致謝。