基于QIIME2平臺擴增子測序分析探究肉牛瘤胃和直腸中厭氧真菌區(qū)系結(jié)構(gòu)的差異

張亞偉 陶薪燕 張月嬌 王月紅 張元慶

(山西農(nóng)業(yè)大學(xué)動物科學(xué)學(xué)院，晉中 030801)

新毛孢菌門(Neocallimastigomycota)真菌廣泛存在于草食動物尤其是大型哺乳動物的腸道微生態(tài)系統(tǒng)中，是一類嚴(yán)格厭氧型真菌。它們憑借自身強大的纖維降解酶體系[1]，可以分別提高7%～9%的腸道纖維消化率和高達(dá)40%的飼糧采食量[2]，因此在反芻動物營養(yǎng)研究中受到普遍關(guān)注，近年來也在生物技術(shù)領(lǐng)域得到了廣泛的研究[3]。截至目前，共有20個屬的厭氧真菌得以鑒定和命名[4-6]，而非培養(yǎng)學(xué)基礎(chǔ)的環(huán)境分子微生物學(xué)研究表明，尚存在更多屬水平的厭氧真菌有待分離、鑒定和描述[7]。隨著研究的不斷深入，許多研究發(fā)現(xiàn)，腸道厭氧真菌在分類學(xué)或功能水平上的變化與甲烷生成或剩余采食量密切相關(guān)[8-9]，因此準(zhǔn)確評估腸道厭氧真菌的區(qū)系組成和多樣性對評價其與宿主動物生產(chǎn)性能之間的聯(lián)系至關(guān)重要。

迄今，最常用的探究腸道乃至環(huán)境中微生物區(qū)系組成和多樣性的方法是標(biāo)記基因基礎(chǔ)的擴增子測序，即宏分類組學(xué)分析[10]。內(nèi)轉(zhuǎn)錄間隔區(qū)1(ITS1)是環(huán)境中厭氧真菌分類學(xué)和多樣性研究常用的標(biāo)記基因區(qū)域，然而厭氧真菌中該基因區(qū)長度的同源異質(zhì)性和單核苷酸多態(tài)性的存在，可能導(dǎo)致經(jīng)典的以操作分類單元(operational taxonomic unit，OTU)為基礎(chǔ)的分類學(xué)分析方法產(chǎn)生很大的誤差，同時也會受其他高通量數(shù)據(jù)處理過程的影響[6,11]。QIIME2是一個免費、開源且可擴展的微生物組信息學(xué)分析平臺，它可以完成從原始測序數(shù)據(jù)至統(tǒng)計分析的全部微生物組高通量數(shù)據(jù)處理過程，同時可以實現(xiàn)單體代表性序列基礎(chǔ)的分類學(xué)注釋和分析[12]，后者在一定程度上可以減少序列相似性聚類產(chǎn)生OTU過程導(dǎo)致的誤差。然而，目前尚沒有針對腸道厭氧真菌群落的QIIME2平臺擴增子測序分析流程的研究報道。

另外，許多研究表明，碳源和宿主品種對厭氧真菌的區(qū)系結(jié)構(gòu)具有顯著的影響，即厭氧真菌群落產(chǎn)生了生態(tài)位分化[13-14]。也就是說，在不同的生態(tài)位中存在著不同的優(yōu)勢厭氧真菌，具有不同的生物學(xué)特性，同時發(fā)揮著不同的生態(tài)學(xué)功能和作用。然而，迄今我們對其生物學(xué)特性的理解多來自于瘤胃厭氧真菌，因此了解厭氧真菌在草食動物前后腸道內(nèi)區(qū)系結(jié)構(gòu)的差異，將有助于加深對厭氧真菌生物學(xué)特性及其在腸道微生態(tài)營養(yǎng)中作用的理解，例如不同種類的厭氧真菌菌體或孢子(囊)被宿主動物消化的程度是否相同等。為此，本研究擬介紹一種基于QIIME2平臺的厭氧真菌擴增子測序分析流程，并探究肉牛瘤胃和直腸中厭氧真菌區(qū)系結(jié)構(gòu)的差異。

1 材料與方法

1.1 試驗動物管理和樣品采集

選擇6頭14月齡、裝有永久性瘤胃瘺管的晉南牛[體重(350±20)kg]作為試驗牛，飼養(yǎng)于山西農(nóng)業(yè)大學(xué)動物科學(xué)學(xué)院畜牧獸醫(yī)研究所試驗基地。所有試驗用瘺管牛均飼喂相同的飼糧，該飼糧由75%的全株玉米青貯和25%的精料補充料混合配制而成(干物質(zhì)基礎(chǔ))。精料補充料購自山西易大飼料有限公司，由玉米、菜籽粕、糖蜜、石粉、磷酸氫鈣、氯化膽堿、蛋氨酸、復(fù)合維生素和復(fù)合微量元素等原料配制而成。全混合日糧的主要營養(yǎng)成分含量為粗蛋白質(zhì)10.9%、粗脂肪3.1%、粗纖維24.9%、無氮浸出物52.3%和粗灰分8.7%(干物質(zhì)基礎(chǔ))。試驗動物每日08：00和17：00各飼喂1次，自由飲水。試驗動物的飼養(yǎng)管理及后續(xù)的樣品采集方案均得到了山西農(nóng)業(yè)大學(xué)實驗動物倫理委員會的批準(zhǔn)。

預(yù)飼15 d后，于第16天晨飼前分別經(jīng)瘤胃瘺管和直腸直接采集6頭試驗牛的腸道內(nèi)容物樣品，所采集的瘤胃內(nèi)容物和直腸糞便均為原始樣品，不進行固液分離。所有樣品均直接置于已滅菌離心管中，隨即液氮冷凍，-80 ℃保存?zhèn)溆谩?/p>

1.2 DNA的提取

瘤胃內(nèi)容物和直腸糞便樣品中全基因組DNA的抽提參照Yu等[15]的方法，使用糞便微生物DNA提取試劑盒(Omega Bio-Tek，Norcross，GA，美國)進行。DNA的純度和濃度使用NanoDrop2000(NanoDrop Technologies，Wilmington，DE，美國)進行檢測，完整性使用10 g/L瓊脂糖凝膠電泳進行分析。純度和完整性合格的DNA樣品方可用于后續(xù)擴增子文庫構(gòu)建和高通量測序。

1.3 擴增子文庫構(gòu)建和高通量測序

為了特異性地評估腸道厭氧真菌的群落結(jié)構(gòu)，本研究選用厭氧真菌特異性引物對[16-17]MN100F(5′-TCCTACCCTTTGTGAATTTG-3′)和MNGM2R(5′-CTGCGTTCTTCATCGTTGCG-3′)對真菌ITS1區(qū)域進行PCR擴增。擴增過程在PCR循環(huán)溫控儀(ABI GeneAmp?9700，Applied Biosystems，新加坡)中進行，具體程序如下：94 ℃預(yù)變性2 min，然后進行35個溫控循環(huán)(94 ℃變性30 s，52 ℃退火30 s，72 ℃延伸30 s)，最后72 ℃穩(wěn)定延伸7 min。PCR反應(yīng)體系為：2×Pro Taq預(yù)混液10 μL，上、下游引物溶液(5 μmol/L)各0.8 μL，模板DNA 10 ng，最后用ddH2O補足至20 μL。每個樣本3個重復(fù)。

將同一樣本的3個PCR重復(fù)產(chǎn)物混合并使用20 g/L瓊脂糖凝膠電泳進行富集，隨后切取目的條帶使用DNA凝膠提取試劑盒(Axygen Biosciences,Union City,CA,美國)進行回收純化。純化后的擴增子樣品濃度使用熒光定量儀(QuantusTM，Promega Corporation，美國)進行檢測，大小和完整性使用20 g/L瓊脂糖凝膠電泳進行檢測。條帶大小正確、濃度合適的擴增子樣品方可進行后續(xù)的建庫和測序操作。

擴增子測序文庫使用NEXTflexTM快速DNA測序試劑盒(Bioo Scientific，美國)按照廠家使用說明進行構(gòu)建，構(gòu)建好的測序文庫使用Illumina MiSeq PE300平臺(Illumina Inc.，San Diego，CA，美國)進行高通量測序，產(chǎn)生長度為300 bp的雙端序列文件用于后續(xù)分析。測序過程由上海美吉生物醫(yī)藥科技有限公司完成。

1.4 測序數(shù)據(jù)處理

擴增子測序數(shù)據(jù)使用QIIME2平臺(版本2021.08)進行處理，具體流程參見圖1。首先，使用q2-tools插件中的import方法將下機數(shù)據(jù)集導(dǎo)入至QIIME2平臺中；隨后，使用q2-cutadapt插件中的trim-paired方法在默認(rèn)條件下將barcode和primer堿基序列從雙端序列中移除[18]；最后，使用q2-dada2插件中集成的denoise-paired方法對雙端序列進行降噪、合并、去嵌合體和去重復(fù)等處理[19]，該過程設(shè)置參數(shù)為移除雙端序列中第203位后的所有堿基序列(由于測序質(zhì)量下降)，并獲得代表性序列(amplicon sequence variants，ASVs)及其頻率分布表。

圖1 厭氧真菌擴增子測序數(shù)據(jù)處理流程

本研究采用機器學(xué)習(xí)的方法進行ASVs的分類學(xué)注釋[20]，參比數(shù)據(jù)庫及分類學(xué)注釋信息選用UNITE真菌ITS數(shù)據(jù)庫[21](版本8.3，2021-10-05發(fā)布)。具體過程如下：首先使用q2-tools插件中的import方法將99% OTU參比數(shù)據(jù)集和相應(yīng)的分類學(xué)注釋信息文檔導(dǎo)入至QIIME2平臺，隨后使用q2-feature-classifier插件中集成的fit-classifier-naive-bayes方法將二者訓(xùn)練成樸素貝葉斯分類注釋文件[22]。然后，基于該注釋文件，根據(jù)機器學(xué)習(xí)的原理，使用q2-feature-classifier插件中的classify-sklearn方法對ASVs進行分類學(xué)注釋[22]。ASVs及其頻率分布和分類學(xué)注釋信息用于下游多樣性和統(tǒng)計分析。

1.5 菌群多樣性分析

對于菌群多樣性分析，首先基于1.4得到的ASVs利用q2-phylogeny插件中集成的align-to-tree-mafft-fasttree處理流程生成系統(tǒng)發(fā)育樹[23]；隨后利用q2-diversity插件中集成的alpha-rarefaction方法分析alpha多樣性指標(biāo)隨測序深度變化的趨勢(即稀疏度曲線分析)，并設(shè)置最大測序深度為33 627，以確保所有樣品的測序深度足以用于下游分析菌群多樣性。隨后，以系統(tǒng)發(fā)育樹和1.4得到的ASVs頻率分布表為基礎(chǔ)，利用q2-diversity插件中集成的core-metric-phylogenetic方法流程計算香農(nóng)指數(shù)(Shannon index)、費氏系統(tǒng)發(fā)育學(xué)多樣性指數(shù)(Faith PD)、可觀測特征序列(observed futures)、辛普森指數(shù)(Simpson index)和皮氏均勻度指數(shù)(Pielou evenness)等alpha多樣性指標(biāo)和Unweighted UniFrac和Weighted UniFrac距離等beta多樣性指標(biāo)，該過程設(shè)置所有樣品的測序深度均為33 627。

1.6 統(tǒng)計分析

厭氧真菌生物學(xué)分類信息在屬水平上進行分析，并定義相對豐度在至少1個樣本中大于0.01%的菌屬是可被鑒定的，而平均相對豐度大于0.1%且在每個組半數(shù)以上的樣本中存在的菌屬是可被檢測的，只有可被檢測的菌屬才被用于組間差異分析。菌群alpha多樣性參數(shù)和屬水平相對豐度的組間差異采用Kruskal-Wallis檢驗進行分析，顯著性水平為P<0.05，當(dāng)0.05≤P<0.10時認(rèn)為有顯著性變化趨勢，該過程使用R(版本4.1.2)中的rcompanion軟件包(版本2.4.13)完成。此外，各樣本beta多樣性的組間差異使用主坐標(biāo)分析(PCoA)進行分析和展示，顯著性采用置換多元方差分析(PERMANOVA)進行檢驗[24]，該過程使用qiime2R軟件包(版本0.99.20)完成。

2 結(jié)果與分析

2.1 測序情況介紹

擴增子測序共產(chǎn)生813 343對雙端序列，平均每個樣本(67 779±10 392)對雙端序列。經(jīng)序列降噪、合并和去嵌合體后，共產(chǎn)生676 759條、平均(56 397±9 867)條有效ITS1序列，用于后續(xù)分析。有效序列經(jīng)去除重復(fù)序列后，共產(chǎn)生256條代表性ITS1序列，即ASVs，序列平均長度(233±16)個堿基。所有樣本的覆蓋度指數(shù)(Good’s coverage)均為1.000，表明本研究的測序深度足以覆蓋全部厭氧真菌群落。該結(jié)果也可由alpha多樣性稀疏度曲線分析(圖2)進一步證明：當(dāng)測序深度大于15 000條序列時，所有樣本的香農(nóng)指數(shù)和費氏系統(tǒng)發(fā)育學(xué)多樣性指數(shù)的稀疏度曲線均趨于水平，即厭氧真菌物種豐度和系統(tǒng)發(fā)育學(xué)多樣性均不再增加，而該值遠(yuǎn)小于單樣本最小有效ITS1序列數(shù)(33 627)，表明所有樣本的測序深度均足以覆蓋全部厭氧真菌群落。

F2007、F2022、F2023、F2201、F2202和F2203代表直腸糞便樣本，R2007、R2022、R2023、R2201、R2202和R2203代表瘤胃內(nèi)容物樣本。

2.2 腸道內(nèi)容物采樣位點對厭氧真菌群落多樣性的影響

表1列出了瘤胃和直腸內(nèi)容物中厭氧真菌群落的alpha多樣性指標(biāo)。由表可知，除了費氏系統(tǒng)發(fā)育學(xué)多樣性指數(shù)外，其他厭氧真菌群落alpha多樣性指標(biāo)在瘤胃和直腸之間均沒有統(tǒng)計學(xué)顯著差異(P>0.05)；直腸糞便中厭氧真菌群落費氏系統(tǒng)發(fā)育學(xué)多樣性指數(shù)顯著大于瘤胃內(nèi)容物(P<0.05)。

表1 不同腸道位點厭氧真菌群落的alpha多樣性指標(biāo)

主坐標(biāo)分析顯示，基于Unweighted UniFrac(圖3-A)和Weighted UniFrac(圖3-B)距離的主坐標(biāo)1(PCo1)和主坐標(biāo)2(PCo2)共分別能解釋總變異的70.20%和69.92%；同時，基于2種距離得到的瘤胃內(nèi)容物和直腸糞便樣品真菌區(qū)系結(jié)構(gòu)明顯分離，該現(xiàn)象也由PERMANOVA分析結(jié)果(P<0.05)進一步證實。值得注意的是，瘤胃內(nèi)容物中厭氧真菌群落beta多樣性置信橢圓均小于直腸糞便樣品，即厭氧真菌群落結(jié)構(gòu)在瘤胃中的變異程度更低。

PERMANOVA即置換多元方差分析；藍(lán)色和紅色橢圓分別代表糞便和瘤胃內(nèi)容物厭氧真菌群落beta多樣性的95%置信區(qū)間。

2.3 腸道內(nèi)容物采樣位點對屬水平厭氧真菌相對豐度的影響

菌群注釋分析顯示，本試驗共從肉牛腸道內(nèi)容物樣品中鑒定出7個屬水平腸道真菌群落，均處于可檢測水平。其中，厭氧真菌屬(Anaeromyces)、瘤胃壺菌屬(Piromyces)、盲腸真菌屬(Caecomyces)、腸道真菌屬(Cyllamyces)、奧賓氏真菌屬(Orpinomyces)和未歸類新毛孢菌科(unclassified Neocallimastigaceae)等6個屬水平菌群為已確定厭氧真菌，分別占真菌群落總量的23.3%、16.8%、27.6%、25.6%、5.4%和0.6%。另外，盲腸真菌屬、腸道真菌屬、厭氧真菌屬和瘤胃壺菌屬為優(yōu)勢菌屬，共約占全部厭氧真菌群落的93.2%。差異豐度分析顯示，直腸糞便樣品中盲腸真菌屬和未歸類真菌(unclassified fungi)的平均相對豐度分別為33.7%和1.7%，較瘤胃內(nèi)容物樣品中的21.4%和0具有升高的趨勢(0.05≤P<0.10)；更為嚴(yán)格的線性判別效應(yīng)尺寸分析(LEfSe)顯示，上述2個菌群的線性判別分析(LDA)得分均大于4.0，但僅未歸類真菌的相對豐度差異達(dá)到了顯著性水平(P<0.05)。

3 討 論

本研究介紹了一種特異性的分析腸道厭氧真菌區(qū)系結(jié)構(gòu)的擴增子分析流程，該流程基于QIIME2平臺，利用DADA2對高通量數(shù)據(jù)進行降噪以獲得所有單體代表性序列集，并利用該代表性序列集基于機器學(xué)習(xí)算法進行分類學(xué)注釋。與經(jīng)典的OTU聚類模式分析流程相比，本研究所用流程既不需要獨立進行序列過濾、雙端序列合并、去嵌合體和去重復(fù)序列等數(shù)據(jù)處理過程，也無需先對單體代表性序列進行聚類以獲得OTU后再進行分類學(xué)注釋，操作相對較為簡便，因此由所選數(shù)據(jù)處理方法不同導(dǎo)致的誤差也相對較少，有利于不同研究間的橫向?qū)Ρ?。更為重要的是，OTU是基于序列相似度進行聚類而產(chǎn)生的，并非生物學(xué)意義上的分類水平，在聚類的過程中不可避免的會導(dǎo)致某些原本非同一物種的單體代表性序列被人為的歸為同一分類簇，造成系統(tǒng)發(fā)育信息丟失[19]。因此，對于同源異質(zhì)性和單核苷酸多態(tài)性普遍存在的厭氧真菌ITS1區(qū)域[6,16]，OTU聚類過程造成的系統(tǒng)發(fā)育分類信息損失理論上會更為嚴(yán)重。本研究所用數(shù)據(jù)處理流程是在單體代表性序列的基礎(chǔ)上進行分類學(xué)注釋分析的，因此可以在一定程度上減少由于OTU聚類導(dǎo)致的系統(tǒng)發(fā)育學(xué)分類信息的損失，理論上更適用于腸道厭氧真菌群落結(jié)構(gòu)的分析。值得一提的是，目前有關(guān)腸道真菌的研究多使用通用引物將包括需氧真菌在內(nèi)的所有真菌生物標(biāo)記物進行全克隆測序[25-27]，僅有少數(shù)研究專注于厭氧真菌群落[28-29]，考慮到腸道內(nèi)環(huán)境特點以及厭氧真菌和需氧真菌各自的生活習(xí)性和生物學(xué)特性的較大差異，我們認(rèn)為特異性地構(gòu)建厭氧真菌測序文庫對探究厭氧真菌在草食動物腸道營養(yǎng)中的生物學(xué)功能更具意義。

*代表該菌屬在瘤胃和直腸內(nèi)容物中的相對豐度具有差異變化的趨勢(0.05≤P<0.10)。

在該分析流程的基礎(chǔ)上，本研究利用擴增子測序進一步比較了晉南牛瘤胃和直腸中厭氧真菌區(qū)系結(jié)構(gòu)的差異。alpha多樣性的組間差異分析表明，厭氧真菌代表性序列的豐富度和均勻度在瘤胃和直腸中均沒有統(tǒng)計學(xué)差異，說明絕大多數(shù)的瘤胃真菌均可能存在特殊的休眠結(jié)構(gòu)[30]，該結(jié)構(gòu)可以抵抗腸道消化過程而到達(dá)后腸道并再次萌發(fā)和定植。同時，直腸中真菌群落的系統(tǒng)發(fā)育學(xué)多樣性要顯著高于瘤胃中，即直腸中真菌菌群間的親緣關(guān)系相對較遠(yuǎn)；此外，beta多樣性的組間差異分析證實，瘤胃和直腸內(nèi)容物中真菌區(qū)系結(jié)構(gòu)具有統(tǒng)計學(xué)差異。鑒于瘤胃和直腸內(nèi)環(huán)境的差異，本研究結(jié)果進一步說明，除了可利用碳源[13]和宿主動物類型[14,29]外，即便是同一動物個體的前后腸道中，真菌群落也表現(xiàn)出了明顯的生態(tài)位分化。與本研究結(jié)果相似，Mura等[31]發(fā)現(xiàn)，厭氧真菌群落結(jié)構(gòu)在馬消化道的不同區(qū)域間也具有明顯的差異。類似的生態(tài)位分化現(xiàn)象在細(xì)菌等其他類型腸道微生物中也有報道[32]。

由菌群注釋分析得到的屬水平厭氧真菌中，瘤胃壺菌屬和盲腸真菌屬為單中心厭氧真菌，共分別占瘤胃和直腸內(nèi)容物真菌群落的37.8%和50.9%；厭氧真菌屬、腸道真菌屬和奧賓氏真菌屬為多中心厭氧真菌，共分別占瘤胃和直腸內(nèi)容物真菌群落的61.2%和47.3%，表明多中心型厭氧真菌在瘤胃內(nèi)環(huán)境中數(shù)量相對更多，而單中心型厭氧真菌則更適應(yīng)于后腸道內(nèi)環(huán)境。結(jié)合直腸中盲腸真菌屬的相對豐度較瘤胃中有升高的趨勢，同時考慮到反芻動物后腸道內(nèi)容物中較高比例的難降解飼料組分，提示單中心厭氧真菌尤其是盲腸真菌屬可能有用作肉牛直飼型微生物或厭氧發(fā)酵飼料(如青貯飼料)接種劑進而改善飼料消化的潛力。與本研究結(jié)果相似，Guo等[29]在牦牛瘤胃中也鑒定出了本研究中所有注釋出的屬水平厭氧真菌，同時還鑒定出了新毛孢菌屬(Neocallimastix)和駱駝?wù)婢鷮?Oontomyces)的存在；彭全輝等[27]在西雜牛和犏牛雜交牛的瘤胃中也鑒定到了瘤胃壺菌屬、奧賓氏真菌屬和新毛孢菌屬的存在。然而與本研究不同的是，這2項研究中優(yōu)勢菌屬分別為瘤胃壺菌屬和新毛孢菌屬，產(chǎn)生這種差異的原因可能是測序文庫構(gòu)建和數(shù)據(jù)處理方法的差異，以及宿主動物類型和所采食飼糧的不同所致。

4 結(jié) 論

① 本研究介紹了一種適用于腸道厭氧真菌擴增子測序數(shù)據(jù)的處理流程，該流程簡便且可操作性強，便于研究結(jié)果間的橫向比較，同時理論上可以減少經(jīng)典的OTU聚類方法導(dǎo)致的系統(tǒng)發(fā)育學(xué)信息丟失。

② 厭氧真菌群落在肉牛瘤胃和直腸中具有一定的生態(tài)位分化，該現(xiàn)象主要是由盲腸真菌屬為主的單中心厭氧真菌所引起的。