宏基因组测序分析男性运动员肠道菌群物种组成及代谢通路特点

金泽宇, 李威, 孙宝林

(中国科学技术大学 先进技术研究院卓源健康人体微生物组研究联合实验室， 合肥 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 研究意义

整体实验结果系统地分析了运动员肠道菌群物种组成在各个物种水平上与对照学生组的差异，详细地描述了运动员肠道菌群结构的特点。同时在物种组成研究结果的基础上，补充了肠道菌群功能的研究，获得了运动员肠道菌群不同于对照学生组的代谢功能，并将运动员肠道菌群的物种组成和功能联系在了一起，全面地阐释了运动员肠道菌群的特点，为特定人群肠道菌群宏基因组的研究提供了参考。