宏基因組測序分析男性運動員腸道菌群物種組成及代謝通路特點

金澤宇, 李威, 孫寶林

(中國科學技術大學 先進技術研究院卓源健康人體微生物組研究聯合實驗室， 合肥 230088)

定植在人體內細菌的數量至少與人體體細胞數相同，而人體微生物編碼的基因數量則遠超人體基因組[1]。據估計有500～1000種細菌同時存在于人體內，其中的亞種數目更是超過這一數量級[2]。每株細菌的基因組包含上千種基因，提供了豐富的基因多樣性，從而使人體微生物組具備了更強的適應性[3-4]。腸道是人體菌群定植的主要部位，腸道菌群可以參與宿主代謝[5]、抑制病原微生物生長[6]等多種生理過程。不同的人群由于受其基因型[7]、飲食[8-9]和健康狀況[10]等不同因素的影響，腸道菌群的組成具有很大差別。目前雖有一些針對特定人群的腸道菌群研究，但多為按照年齡[11]、地區[12]或是疾病類型[13-14]來進行劃分，運動員群體具有生活飲食規律、運動訓練時間充足的特點，其腸道菌群的特征能夠在一定程度上反映規律的生活飲食和充分的運動訓練對腸道菌群的影響。本研究選取男性田徑運動員的腸道菌群樣本進行研究，其每日飲食標準如下：熱能3500～4500 kCal，蛋白質80～120 g，脂肪80～120 g，碳水化合物615～735 g，鈣600 mg，磷1000 mg，鐵12 mg，維生素A 1000 μg視黃醇當量，維生素B12 mg，維生素B22 mg，煙酸20 mg，維生素C 100 mg。采集其腸道微生物樣本后使用Illumina HiSeq平臺進行宏基因組測序，同時為排除性別對腸道菌群的影響，以相同地區同等年齡的普通男性大學生腸道菌群作為對照，在不同分類水平上對男性運動員的腸道菌群結構進行分析，并尋找其代謝通路相比對照組發生的改變，為男性運動員腸道菌群的物種組成及功能研究提供了參考。

1 材料與方法

1.1 志愿者招募及樣本采集

招募35名年齡為18～20周歲的健康男性運動員，身高體重指數(BMI)符合正常標準。作為對照，在相同地區招募35名同等年齡的普通健康男性大學生。采集糞便樣本，裝入50 mL無菌離心管，采用低溫運輸(0 ℃)送至北京微生物所進行測序。

1.2 序列分析軟件及流程

使用KneadData軟件(版本為0.7.2)、Trimmomatic軟件[15]、Bowtie2(版本為2.3)[16-17]去除低質量堿基以及人體基因組污染序列，隨后使用HUMAnN2流程進行分析(版本為0.11.2)[18]。首先使用MetaPhlAn2[19]下載chocophlan full數據庫用于物種信息的比對，獲得所有樣本的物種組成信息；其次使用DIAMOND(版本為0.8.36)[20]及UniRef90數據庫對序列進行比對，獲得所有樣本的基因家族豐度表，以及對應的代謝通路豐度表。

2 結果與分析

2.1 所有樣本的物種豐度分析

所有60個樣本中共鑒定出8個門、13個綱、18個目、44個科及109個屬，共326種已知細菌。圖1-A展示了全部60個樣本豐度前25的屬信息，其余屬歸為Other。從圖1中可以看到運動員樣本(J002-J059)與普通男性大學生樣本(S002-S034)之間未見明顯多樣性差別。進一步從兩組樣本的種水平維恩圖中可見運動員組(M)共鑒定出284種細菌，而學生組(S)共鑒定出286種細菌，其中244種細菌為兩組共有，占總菌種數的74.8%；運動員組包含40種特有菌種，占總菌種數的12.3%；而學生組包含42種特有菌種，占總菌種數的12.9%(圖1-B)。

根據獲得的物種組成豐度表，分析了兩組樣本不同水平的優勢細菌。從圖1-C中可以看到，兩組樣本中平均豐度最高的3個門類為厚壁菌門(Firmicutes)、放線菌門(Actinobacteria)以及擬桿菌門(Bacteroidetes)。其中厚壁菌門的平均豐度最高，大部分運動員組中相對豐度在60%～85%，而學生組中的相對豐度為78%～90%。相對豐度次之的門為放線菌門，運動員大部分相對豐度為8%～33%，而學生組多數為6%～20%。再次為擬桿菌門，運動員組大部分為4%～12%，高于學生組。總體來說，兩組樣本中均遵循某些物種的含量較高，而大多數物種的豐度較低的分布規律(如圖1-D和圖1-E中呈現出的長尾效應)。

A：所有樣本豐度前25的屬相對豐度條形圖;B：運動員組(M)和學生組(S)種水平維恩圖;C：門水平豐度前3的優勢物種線箱圖;D：屬水平豐度前15的優勢物種線箱圖;E：種水平豐度前20的優勢物種線箱圖

圖1所有樣本的物種組成情況

Figure 1 Taxonomic composition of all samples

2.2 運動員和學生樣本物種組成的主成分分析

獲得所有樣本的物種組成相對豐度后，對所有樣本的物種組成進行了主成分分析(principal components analysis，PCA)。從圖2中可以看到，3條主成分坐標軸分別為PC1(39.2%)、PC2(16.5%)、PC3(13.1%)，即每條坐標軸分別能解釋所有樣本中39.2%、16.5%及13.1%的方差。PC1和PC2構成的圖中，運動員組(M)和學生組(S)的樣本間大部分重合，說明兩組樣本在區分度最大的平面內(54.7%)的差異較小，物種組成相似度較高。而在PC1和PC3構成平面內，雖然兩組間在右側還有部分樣本重疊，但運動員組的樣本大部分分布在中下部，而學生組的樣本大部分分布在中上部，已經可以逐漸彼此區分，這種分離的趨勢同樣呈現在PC2和PC3構成的二維平面內。總體來說，兩組樣本在主成分貢獻率較高的平面組成相似，而在貢獻率較低的平面內逐漸分離，說明兩組樣本在總體的物種組成上大部分相似，但仍有一些物種不同。

圖2 所有樣本物種組成的PCA圖

2.3 STAMP分析運動員和學生差異菌屬

STAMP可用于提供不同樣本間物種和功能組成分析的統計學假設檢驗以及作圖[21]。運動員組與學生組的差異菌屬見圖3-A，共篩選出10種組間分布差異顯著的菌屬，運動員組與學生組的豐度依次為：雙歧桿菌屬(Bifidobacterium)為12.1045%和6.1745%、Dorea屬為4.2856%和8.2763%、擬桿菌屬(Bacteroides)為3.3391%和1.2562%、鏈球菌屬(Streptococcus)為0.7833%和2.1674%、Bilophila屬為0.1437%和0.0171%、Burkholderiales_noname屬為0.0486%和0.0036%、Parasutterella屬為0.0493%和0.0048%、孿生菌屬(Gemella)為0.0061%和0.0483%、Holdemania屬為0.0121%和0.0437%、丙酸桿菌屬(Propionibacterium)為0.1278%和0.0096%。其中前4種菌屬所有樣本的平均相對豐度在1%以上，而后6種菌屬的平均相對豐度較低。

雙歧桿菌屬可發酵腸道內的碳水化合物，尤其是寡糖類分子[22]。運動員組的平均含量幾乎是學生組的一倍。擬桿菌屬是哺乳動物腸道菌群的主要部分，長期的富含蛋白和動物脂肪的飲食可引起擬桿菌屬含量的上升[23]。擬桿菌屬還可阻止潛在的病原體定植在腸道，但擬桿菌屬同時可能成為腸道菌群內的抗性庫[24]。肉食為主的飲食會使Bilophila屬和擬桿菌屬的含量有所上升，同時腸道內代謝植物多糖的厚壁菌門細菌數量會下降[25]。綜合來看，運動員組的擬桿菌屬和Bilophila屬含量均高于學生組，同時厚壁菌門的含量低于學生組。丙酸桿菌屬是革蘭陽性厭氧菌，具有獨特的羧基轉移酶用以合成丙酸，被用于益生菌產業，也有一些種與人體皮膚感染相關[26]。所有差異菌屬中，運動員組的3種益生菌雙歧桿菌屬、擬桿菌屬以及丙酸桿菌屬的含量都高于學生組。鏈球菌屬是革蘭陽性球菌，總體來說鏈球菌屬對人體的健康狀態影響較為復雜；Dorea屬的功能目前未見報道，運動員組中這兩種菌屬的含量都低于學生組。Burkholderiales_noname屬、Parasutterella屬、孿生菌屬和Holdemania屬等幾種菌屬在運動員組和學生組間中平均豐度較低(小于0.1%)，差異較小。

在得到差異菌屬后，利用差異菌屬的分布繪制了所有樣本的相對豐度熱圖(圖3-B)。從圖3中可以看到，大部分的運動員樣本(藍色)均聚類在左側，證明其差異菌屬的分布模式較為相似，而中間部分多數為學生樣本(橙色)的聚類，最后末尾有部分運動員和學生的樣本相互夾雜。說明僅憑差異菌屬的分布，能區分大部分的運動員樣本和學生樣本，但也有一些樣本的差異菌屬分布模式并不規律。運動員和普通大學生雖然在飲食及生活習慣上存在一定的差異，但均屬于健康人群。鑒于腸道菌群的復雜性和多樣性，無法進行清晰地分組也屬于正常情況。

A：運動員組(M)和學生組(S)差異菌屬。左側為每種差異菌屬在兩組樣本間的平均豐度(Mean proportion)，中間為兩組平均豐度差值(Difference in mean proportions)及95%置信區間(confidence intervals)，右側為每種差異菌屬的P值(P-value)。B：運動員組(M)和學生組(S)差異菌株差異菌屬的豐度熱圖。其中橫軸為所有樣本，縱軸為差異菌屬，均使用Ward算法進行聚類，顏色深度代表每個樣本對應菌屬的相對豐度(relative abundance)，顏色越深相對豐度越高

圖3運動員和學生差異菌屬STAMP分析結果

Figure 3 STAMP analysis of taxonomic composition (genus level)

2.4 LEfSe分析運動員和學生差異菌種

LEfSe可判斷最有可能解釋類別不同的因素，能夠尋找具有統計學差異的生物標識(Biomarker)。圖4-A為從門到種所有水平下運動員組(M)和學生組(S)的線性判別分析結果，圖4-B為根據差異的物種所構建的進化樹(門到屬水平)。結合所有結果可以看到，運動員組含量較高的細菌主要集中在擬桿菌門、變形菌門的β-變形菌綱(Betaproteobacteria)和δ-變形菌綱(Deltaproteobacteria)下；而學生組含量較高的細菌則主要集中在厚壁菌門下。

2.5 所有樣本代謝通路分析

HUMAnN2除了可以采用標記基因檢索已知物種外，還可以翻譯未知物種的序列，與已知的蛋白數據庫進行比對，從而計算基因家族和代謝通路的豐度(包括群體和物種水平)。經過與UniRef90數據庫進行比對后，一共獲得了1 091 690個基因家族的信息。進一步通過MetaCyc數據庫計算出代謝通路后，在所有樣本中共發現了409條代謝通路(圖5-A)，其中382條代謝通路為兩組樣本共有，占通路總數的93.4%；運動員組(M)包含10條特有代謝通路，占通路總數的2.4%；而學生組包含17條特有代謝通路，占通路總數的4.2%。

獲得所有樣本的代謝通路相對豐度后，對所有樣本的代謝通路組成進行了主成分分析。從圖5-B中可以看到，3條主成分坐標軸分別為PC1(23.1%)、PC2(16.0%)和PC3(12.4%)，即每條坐標軸分別能解釋所有樣本中23.1%、16.0%及12.4%的方差。與物種組成PCA圖相比，代謝通路的主成分軸所能解釋的樣本方差更少。

2.6 STAMP分析運動員和學生差異代謝通路

使用STAMP獲得11條運動員組(M)與學生組(S)的差異代謝通路(圖6-A)，每條通路對應的功能簡介見表1。

A：運動員組(M)和學生組(S)各物種水平下LDA結果。取LDA score=2為閾值篩選差異結果，按LDA score由大到小排列。LDA score>0說明在學生組中含量較高，而LDA score<0說明在運動員組中含量較高。B：根據LDA結果構建的進化樹(門到屬水平)。其中綠色節點代表該水平的細菌在學生組(S)中含量較高，而紅色節點代表該水平的細菌在運動員組(M)中含量較高，黃色節點代表該水平的細菌在兩組之間的分布無顯著性差異

圖4運動員和學生差異菌種LEfSe分析結果

Figure 4 LEfSe analysis of taxonomic composition

A：運動員組(M)和學生組(S)代謝通路維恩圖。B：所有樣本代謝通路組成的PCA圖

圖5 所有樣本的代謝通路分析結果

在獲得差異代謝通路后，利用差異代謝通路的分布繪制了所有樣本的相對豐度熱圖(圖6-B)。從圖中可以看到，大部分的運動員樣本(藍色)差異代謝通路的分布較為相似，而中間部分多數樣本為學生樣本(橙色)的聚類，最后左側部分同樣出現了運動員和學生樣本相互夾雜的情況。與物種組成熱圖類似，僅憑差異的代謝通路無法完整的區分兩組樣本。

A：運動員組(M)和學生組(S)差異代謝通路。左側為每條差異代謝通路在兩組樣本間的平均豐度(Mean proportion)，中間為兩組平均豐度差值(Difference in mean proportions)及95%置信區間(Confidence intervals)，右側為每條差異代謝通路的P值(P-value)。B：運動員組(M)和學生組(S)差異菌株差異代謝通路的豐度熱圖。其中橫軸為所有樣本，縱軸為差異代謝通路，均使用Ward算法進行聚類，顏色深度代表每個樣本對應菌屬的相對豐度(relative abundance)，顏色越深相對豐度越高

圖6運動員和學生差異代謝通路STAMP分析結果

Figure 6 STAMP analysis of metabolic pathways

2.7 LEfSe分析運動員和學生差異代謝通路

除了STAMP外，還使用了LEfSe對兩組樣本間的差異代謝通路進行了分析，結果(圖7-A)表明，共有13條代謝通路的LDA score絕對值在2以上，其中運動員組(M)有5條，而學生組(S)有8條，每條通路對應的功能簡介見表1。在獲得LEfSe分析的差異代謝通路后，還總結了這13條代謝通路在所有樣本中的分布情況(圖7-B)。

2.8 兩組間差異代謝通路的物種貢獻組成

結合STAMP和LEfSe的分析結果，選取了在兩種方法下都具有顯著差異的8條通路進行進一步分析。總體來說，8條差異代謝通路中，2條與B族維生素及吡哆醇的合成相關(PWY0-845、PWY-7282)，而這兩條通路在運動員組內的相對豐度要高于學生組，預示運動員腸道內B族維生素的合成能力可能強于普通學生；還有2條通路與IPP合成相關(PWY-6270、PWY-7560)，其中PWY-7560在運動員組內的相對豐度也要高于學生組，預示運動員腸道內合成IPP的能力可能也強于普通學生，但運動員組的PWY-6270豐度有所下降，而PWY-6270的作用是將IPP最終轉化成異戊二烯，此通路的豐度下降可能會造成IPP一定程度的積累，但腸道內IPP的實際含量仍需進一步檢測驗證。IPP可最終合成各種萜類化合物，包括所有類固醇分子的前體羊毛甾醇、維生素A(視黃醇及視黃醛)、胡蘿卜素等重要分子。剩余的GALACTUROCAT-PWY和PWY-6507都涉及碳水化合物的代謝和3-磷酸-D-甘油醛的生成，而且兩者的分布模式呈現很高的一致性，物種貢獻模式也很相近，均由柔嫩梭菌貢獻絕大部分豐度，這兩條通路在運動員組中的豐度也要高于學生組。這兩條通路貢獻的3-磷酸甘油醛可以作為很多下游代謝通路的原料，包括IPP的合成。而剩余的NAD合成通路(NAD-BIOSYNTHESIS-II)和L-丙氨酸合成通路(PWY0-1061)則是在學生組樣本中的豐度更高。

A：運動員組(M)和學生組(S)各物種水平下LDA結果。取LDA score=2為閾值篩選差異結果，按LDA score由大到小排列，LDA score>0說明在學生組中含量較高，而LDA score<0說明在運動員組中含量較高。B：LEfSe差異代謝通路各樣本分布圖。每條差異代謝通路在所有樣本中的相對豐度(Relative abundance)，左側為運動員組(M)，右側為學生組(S)，黑色實線代表該組中所有樣本相對豐度的平均值，黑色虛線代表相對豐度的中位數

圖7運動員和學生差異代謝通路LEfSe分析結果

Figure 7 LEfSe analysis of metabolic pathways

3 討論與結論

3.1 運動員差異菌群

采用宏基因組學策略，分析了32名男性運動員以及28名同齡男性大學生的腸道菌群，發現運動員腸道內的厚壁菌門含量較低，而擬桿菌屬和Bilophila屬含量較高。同時運動員腸道的擬桿菌屬、雙歧桿菌屬以及丙酸桿菌屬等益生菌的含量較高，而鏈球菌屬、孿生菌屬、Dorea屬和Holdemania屬等細菌的含量較低。

3.2 運動員差異代謝通路

發現運動員腸道菌群的B族維生素合成(主要為吡哆醇相關通路，PWY0-845和PWY-7282)的能力更強，而吡哆醇作為輔酶可以參與體內多種氨基酸、葡萄糖和脂質的代謝反應。此外，運動員體內代謝碳水化合物生成3-磷酸甘油醛的通路(GALACTUROCAT-PWY和PWY-6507)豐度也有所增加，可能會給下游通路提供更多的代謝底物3-磷酸甘油醛；而下游以3-磷酸甘油醛為起始底物的IPP合成通路(PWY-7560)也有所增強，但IPP轉化異戊二烯的通路(PWY-6270)有所下降，最終可能導致運動員組IPP的含量上升。同時還鑒定出運動員組的NAD合成及L-丙氨酸合成通路有一定程度的下降。以上所有的代謝通路水平變化均為測序數據預測結果，代謝物的體內真實水平仍需后續實驗驗證。

縱軸為每條代謝通路的相對豐度，單位為cpm(copy per million，每百萬堿基拷貝數)；橫軸為每個樣本，左側為運動員組(M)，右側為學生組(S)；每種菌種所貢獻的通路豐度用不同顏色標出，其中灰色代表未知菌種(Unclassified)。取貢獻程度最大的前10種細菌進行標注，其余菌種記為其他(Other)

圖8差異代謝通路的物種貢獻比例

Figure 8 Stratified results of significantly differentially abundant pathways

3.3 研究意義

整體實驗結果系統地分析了運動員腸道菌群物種組成在各個物種水平上與對照學生組的差異，詳細地描述了運動員腸道菌群結構的特點。同時在物種組成研究結果的基礎上，補充了腸道菌群功能的研究，獲得了運動員腸道菌群不同于對照學生組的代謝功能，并將運動員腸道菌群的物種組成和功能聯系在了一起，全面地闡釋了運動員腸道菌群的特點，為特定人群腸道菌群宏基因組的研究提供了參考。