基于高速双倍自助法对蒙古马系统发育树的研究

任爱珍，白东义，赵一萍，侯 娜，芒 来*

(1. 内蒙古农业大学理学院，呼和浩特 010018； 2. 内蒙古农业大学动物科学学院内蒙古自治区蒙古马遗传资源保护及马产业工程实验室，呼和浩特 010018)

各种分子系统进化树的构建方法只给出了一个真实进化树的点估计，需要对其可靠性给出概率描述[1]。目前提出的常用的基于评价最大似然树的可靠性方法有:自助法、下平-长谷川法(shimodaira-hasegawa method，SH)、多尺度自助法(multiscale bootstrap method,MBP)、双倍自助法(double bootstrap method，DBP)、高速双倍自助法(speedy double boostrap method，SDBM)[2-5]，其中利用自助法计算出的P-值用自助P-值(bootstrapP-value，BP)表示。但是杨子恒[6-7]在其有关分子进化树的综述和著作中总结并指出，理论研究和计算机模拟比较都表明，BP与它所估计的“重建的树是真实树的概率”间存在着较大的差距，利用BP进行检验在统计上是保守的，即不易得出“估计的系统树受资料显著支持”的结论[8-9]。因此，提出了几种复杂的自助法以得到更高精度的研究结果，其中之一是多尺度自助法，Shimodaira[4]用渐近无偏(approximately unbiased，AU)表示其计算的P-值。这个方法具有3次精度，计算量是O(MB),其中M是重抽样的自助样本尺度的个数(重抽样的自助样本尺度是原样本容量n除以重抽样本容量m的比值)，B是重新采样的样本个数。这个方法已成功地在许多实际问题中得到应用[4]。然而，问题是它的计算量比较大。因此，研究者提出了高速双倍自助法[5]，用高速双倍自助P-值(speedy double boostrapP-value，SDBP)表示其P-值。这个方法克服了计算AU速度慢和因外插估计出的数值不稳定的缺陷，解决了为获得3次精度P-值所需计算量大的难题，且也是具有3次精度而计算量仅为O(B)。在哺乳动物系统树的应用研究表明，3次精度P-值中SDBP计算速度是最快的[5]。

国外关于马系统发育的研究较早，但是对于蒙古马系统全面的研究比较少[10-14]。国内有关马和蒙古马的进化与遗传多态性的研究比较多[15-21]。2006年李金莲[19]对内蒙古锡尼河马、巴尔虎马、乌审马、乌珠穆沁马4个类型的蒙古马及引入品种纯血马、培育品种三河马和地方品种广西矮马共309个样品的mtDNA D-loop序列利用非加权组平均法(unweighted pair-group method with arithmetic means, UPGMA)进行聚类分析得出，乌珠穆沁马和乌审马亲缘关系最近，锡尼河马和巴尔虎马亲缘关系较近，与二者次近的是三河马，它们与广西矮马和纯血马亲缘关系都较远。2010年王小斌等[21]对蒙古马与世界范围内的家马、古马、普氏野马的mtDNA D-loop序列进行系统发育分析，并根据邻接法(neighbor joining，NJ)分析结果，共分出A、B、C、D、E、F、G 7个分支，其中蒙古马分在A、B、C、D、E、F 6个分支中。2013年陈建兴等[20]针对亚欧马、蒙古马、普氏野马的ND4基因序列进行了遗传多样性分析，并根据NJ树分析结果得出，普氏野马和蒙古家马聚在一起，用数据证明了普氏野马和蒙古家马具有较近的亲缘关系。我国对于蒙古马系统发育分析研究有一定的成果，但是建树方法多采用UPGMA法和NJ法。杨子恒[6]指出，由于对绝大多数资料来讲分子钟的假定都不能成立，所以UPGMA法在使用中存在缺陷。另外，杨子恒[6]还指出，最大似然法由于计算量庞大，仅在少量研究中进行了比较，不过这些研究都表明其效率在现有方法中较高，比使用相同模型的NJ法或最小二乘距离矩阵法效率要高[22-24]。因此可以利用最大似然法估计蒙古马系统进化树，对估计出的系统进化树进行可靠性评价，因此，可靠性检验方法也需要进一步完善，推断出更精确的进化关系。

本研究以最大似然法估计蒙古马的系统发育树，然后对估计出的系统树利用高速双倍自助法计算其可靠度SDBP值，并给出了SDBP值最大的系统发育树。这将为蒙古马的适应性进化以及基因结构和功能的研究奠定基础。

1 材料与方法

1.1 试验材料

本研究从美国国家生物技术信息中心(National center of biotechnology information，NCBI)网站中获取了30匹中国蒙古家马与6匹普氏野马的mtDNA D-loop全序列(GenBank登录号见表1)，这些序列包含巴尔虎马、三河马、乌审马、乌珠穆沁马以及锡尼河马5个类群。从每个类群中选取了6个样本。此外，普氏野马(P)有6个样本，共计36匹马。

1.2 最大似然系统发育树的计算

设有一组生物种类数为S，每种生物DNA的长度为n的碱基序列见表2，为了一般性和简单明了,利用了虚拟的DNA碱基序列。在统计学中,记S种生物DNA碱基序列为矩阵Xs×n=(x1,…,xn)，其中xi(i=1,…,n)表示表2中DNA碱基序列的第i位点。S种(不包含外群的种数)生物可有K=(2S-3)!!棵候补无根系统树(本研究所建的树都是无根树)，每棵系统树记为Tk(k=1,…,K)，第k棵系统树的概率函数：

表1蒙古马D-loop区碱基序列的GenBank登录号

Table1GenBankaccessionnumbersofMongolianhorseD-loopsequences

品种Breed数量/匹QuantityGenBank登录号GenBank accession number巴尔虎马Baerhu 6JN790867～ JN790872三河马Sanhe6JN790879～JN790884乌审马Wushen 6JN790885～JN790890乌珠穆沁马Ujimqin6JN790897～JN790902锡尼河马Xinihe6JN790891～JN790896普氏野马Przewalski’s6DQ017347、DQ017373、GU565340、GU565529、JN398402、JN395403

fk(x;θk)

(1)

其中，k表示第k棵系统树，θk表示第k棵系统树的概率函数参向量，其包含碱基置换模型中的参数和树状拓扑中的若干个树枝长度参数，x为DNA碱基序列中的一个位点所表示的碱基列。

表2S种生物DNA碱基序列

Table2DNAbasesequencesofSspecies

生物的序号Number of speciesDNA碱基序列DNA base sequence12345678…n1TCACCTATG2TCACTTAT…G3TCATCTAT…G……………………………STCACCTAT…G

表中列的“1，2，3，…，S”表示生物的序号；行的“1，2，3，…，n”表示DNA碱基序列的位点序号

The "1, 2, 3,…， S" in the columns denote the serial number of species; the "1, 2, 3,…， n" in the rows denote the serial number of sites

根据以往的研究报道设定碱基置换模型中的参数，所以碱基置换模型是已知的，只有树枝长度是未知。本研究中的S种生物指的是S=5个类群的蒙古马，而K=(2·5-3)!!=105是5个类群蒙古马的所有候补树，5个类群蒙古马的DNA碱基序列矩阵为X5×309，其中位点总数n=309(以下计算步骤中都设成S=5,K=105,X5×309)。求解最大似然系统发育树有如下4个步骤。

步骤1：计算每棵系统发育树的最大似然值。在实际的计算中先利用Mega6对表2的DNA碱基序列和外群的DNA碱基序列一同进行多重比对，比对后去掉有空位的列，然后得到处理后的DNA碱基序列，建一个文件夹作为输入文件。S种(不包含外群的种数)生物的(2S-3)!!棵不同树形拓扑按Newick格式全部手写出来，建一个文件夹作为另一个输入文件。利用软件PAML4.9以最大似然估计法计算每棵系统树的对数似然函数：

(2)

(3)

步骤2：建立最大似然矩阵。利用软件PAML4.9计算每棵系统树对于每个位点的对数似然估计值，按行的形式摆放，得到最大对数似然矩阵：

简记为

步骤3：计算最大对数似然向量。利用软件CONSEL将公式(4)中每行数据相加，可得最大对数似然向量：

(5)

1.3 计算系统发育树的可靠度SDBP值

对BP以及AU的计算步骤这里省略掉。系统树可靠性评价首先要提出如下的假设，这个假设是属于多元假设检验问题。

…

第1棵系统树为真实系统树的假设检验：

原假设H1的概率值SDBP的计算步骤：

(6)

(2)重抽样

(7)

(8)

(9)

(10)

(3)计算SDBP值

(11)

其中d是公式(6)表示的量。

2 结 果

2.1 蒙古马系统进化树可靠度分析

将获取的36匹马的D-loop碱基全序列导入Mega6软件中，对所有马匹D-loop区片段序列进行多重序列比对，比对后把空位的位点全部删除，结果每匹马的碱基长度为309 bp。

2.1.1 PAML4.9软件输出最大似然进化树 将多重比对之后的碱基序列和所有的105棵候补树组成的两个文件夹输入PAML4.9软件中的baseml程序，并选择REV碱基替换模型分析。运行程序baseml后可得105棵树每个位点的对数似然值，即公式(4)。这些对数似然值可得出105棵树中的最大似然树的序号。

2.1.2 CONSEL、R3.0中SDBP包给出各系统树的可靠度 将所得105棵树的每个位点的对数似然值，利用CONSEL软件，计算各个系统树的AU值和BP值。最后利用R3.0中SDBP包计算各个树的SDBP值，可得SDBP值最大的蒙古马系统发育进化树。

表3展示了SDBP值大于等于0.05的60棵系统进化树的SDBP值，旁边的数值是AU值和BP值。表3中“排序”表示按每棵系统进化树的最大对数似然值降序生成的树的排序，“树形”表示PAML4.9软件中输入105棵系统进化树建立的树文件中系统树的序号。综合考察结果，按系统树的SDBP值分析系统树，序号为100的系统树的SDBP值最大为0.980。而按系统树的AU值和BP值分析进化树，序号为66和19的可靠度最大，分别为0.717和0.215。采用最大似然法结合SDBP值分析蒙古马的系统拓扑关系是较最大似然法结合渐近无偏(approximately unbiased，AU) 值以及结合自助(bootstrapP-value，BP)值分析蒙古马系统拓扑关系更有效。将t100号、t66号以及t19号系统树的拓扑结构画成二分系统树的树形，得到图1、图2、图3。

表330匹蒙古马的60棵候选树SDBP、AU以及BP值的计算结果

Table3TheresultsofSDBP,AUandBPvaluesof60candidatephylogenetictreesfor30Mongolianhorses

排序Rank树形TreeSDBP值SDBP valueAU值AU valueBP值BP value排序Rank树形TreeSDBP值SDBP valueAU值AU valueBP值BP value11000.9800.6970.0263140.4400.0120.0002660.9330.7170.09632790.3480.0930.0003190.7510.5630.21533700.3510.0970.0004890.7770.5700.161341020.3450.0980.0005300.7690.4990.10135410.3530.0950.0006780.7760.4970.1013690.3390.0810.0007430.7780.4980.10137530.3520.0960.00081030.7680.3660.0033850.3480.0590.0009370.7710.4960.10139630.3270.0000.10110240.7800.1600.00540420.3290.0000.00511490.7740.4960.10141350.3290.0000.10112990.6490.4460.21042520.3370.0000.21013730.6490.4470.20943500.3390.0000.20914820.5960.3740.14844580.3370.0000.14815770.6780.1540.00345930.3400.0000.00316160.5390.1440.0094660.3310.0000.00917550.5580.2340.05047210.3350.0000.05018320.5590.2380.05148830.3360.0000.05119740.5650.2370.05149170.3260.0000.05120760.5590.2370.05150220.3360.0000.05121130.5640.1990.00251200.3330.0000.00222440.5570.1280.00152690.3070.0000.00123940.5580.1340.00053600.3000.0000.00024250.5820.2260.00054230.2960.0000.00025260.5770.2210.00055870.3060.0330.00026920.5730.0020.00056980.3130.0310.00027340.5760.0020.00057720.3120.0320.00028950.5840.2320.00058100.3160.0320.0002970.5800.2310.00059140.3090.0320.00030620.5830.2350.00060120.3120.0320.000

图1 SDBP值最大系统树t100拓扑结构Fig.1 The topology of tree t100 with the highest SDBP (speedy double bootstrap P-value) value

图2 AU值最大系统树t66拓扑结构Fig.2 The topology of tree t66 with the highest AU (approximately unbiased) value

图3 BP值最大系统树t19拓扑结构Fig.3 The topology of tree t19 with the highest BP (bootstrap P-value) value

SDBP值最大的蒙古马系统树(图1)共有3个分支，巴尔虎马与乌审马合并为一个分支，具有较近的遗传关系；乌珠穆沁马与锡尼河马具有较近的遗传关系；三河马与其他4种蒙古马独立地分出来构成一个系统发育分支。SDBP值最大的蒙古马系统树与三河马和其他4类蒙古马具有较远的预测血缘关系是一致的，同时此结果的系统树也具有最大似然值。

2.2 蒙古马NJ树与UPGMA树分析

将多重比对之后的36个碱基序列输入Mega6软件中，利用Mega6软件工具栏中Data下的Setup/Select Taxa and Groups标签对36个碱基序列按表1分组，结果得到带有6组分组信息的36个碱基序列并保存。然后对带有分组信息的36个碱基序列，利用工具栏中Data下的Setup/Compute Between and Groups Means标签计算6组间的距离。其次对组间距离的结果利用工具栏中Phylogeny下的Construct Phylogeny标签采用NJ方法和UPGMA方法分别建立系统树。最后得到蒙古马的NJ系统树和UPGMA系统树，把它们画成二分系统树的树形，得到图4和图5。结合SDBP值分析蒙古马的系统拓扑关系较非加权组平均法(unweighted pair-group method with arithmetic means，UPGMA) 有效，而与邻接法(neighbor joining，NJ)得到了相近的结论。

图4 5个蒙古马类群的NJ树拓扑结构Fig.4 The topology of NJ trees of 5 Mongolian horse subspecies

这次的试验对进化时间没有做分析，只是对拓扑结构做的分析。所有105棵树都是以普氏野马为外群的无根树，所以图1～图5中的进化树是没有时间标尺和时间方向的进化系统树。

图5 5个蒙古马类群的UPGMA树拓扑结构Fig.5 The topology of UPGMA tree of 5 Mongolian horse subspecies

3 讨 论

从统计学假设检验理论来判断，每棵系统发育树的SDBP值相当于原假设Hk的P-值，k=1,…,105。所以如果以显著性水平检验各原假设的话，各原假设的SDBP值大于或等于0.05的60个系统树都不被拒绝。如何从这60个系统树中找出真实的系统树是分子进化学中的一个重要问题。一般采用生物学等其他信息在选择过的范围里再选择。在SDBP值大于或等于0.05的60个系统树都不被拒绝的情况下，我们认为结合一个或几个种群马与其他种群马的遗传关系的已知信息，从60棵系统树中判断出一个符合这个已知信息的系统树作为蒙古马的系统进化树在逻辑与统计学方法上都是正确的。所以根据已有研究报道，三河马的血缘主要是在后贝加尔马和其杂种马的基础上，混入呼伦贝尔草原蒙古马的血液，经几十年精心培育而成[28]。根据三河马这样的血缘关系，预测三河马应该从系统进化发育树首先分离出来独立作为一个分支，而其他4种 蒙古马种群作为另一个分支。从t100、t66以及t19 3个系统树看符合这一条件的只有t100，而t100是最大似然树，其SDBP值也是最大的。所以，笔者认为t100比较真实地反映了蒙古马的系统进化关系。对比3次精度AU值最大的系统树t66与SDBP 值最大的t100，它们都有以乌珠穆沁马与锡尼河马为叶子组成的一个独立小分支，而在BP值最大的t19里没有这样的独立小分支。

对比2006年李金莲[19]用UPGMA的结果与本研究中3次精度SDBP值最大的t100，它们没有共同的由两个叶子组成的独立小分支，在拓扑结构上也不一样。2006年李金莲[19]用UPGMA的结果与3次精度的AU值最大的系统树t66以及BP值最大的t19相比，无论在拓扑结构还是独立分支上也没有相似之处。由于2010年王小斌等[21]是把蒙古马与世界范围内的家马、古马、普氏野马一起进行系统发育分析，即没有进行分类分析，所以本研究结果与王小斌等[21]的NJ树结果没有可比性，陈建兴等[20]的结果也存在同样的问题，因而也无可比性。对比本研究用NJ方法与UPGMA分析出的结果(图4与图5)可知，蒙古马的NJ树与SDBP值最大的t100树，它们的分支在系统树中的顺序完全一样，只是系统树拓扑结构上有些差异。而蒙古马的UPGMA树拓扑关系与BP值最大的t19树非常接近。这说明在每匹马序列长度比较少的情况下NJ树的结果要比UPGMA树的结果更接近SDBP值最大的系统树。综上所述，本研究通过利用30匹蒙古马线粒体DNA中的D-loop高变区序列对5类蒙古马的105个候补系统发育树进行分析，结果表明，采用最大似然法结合SDBP值分析生物的系统拓扑关系是较最大似然法结合AU值以及结合BP值分析生物的系统拓扑关系更有效的方法。同时也比UPGMA法有效，而与NJ法得到的结论相似。

4 结 论

本研究利用30匹蒙古马的线粒体DNA中的D-loop高变区序列对蒙古马的105个候补系统发育树进行了分析，采用最大似然法，并利用生物信息软件Mega6、PAMAL4.9等估计出蒙古马的最大似然树。最后利用CONSEL和R3.0中的SDBP软件包等计算了105个候补系统发育树的SDBP值。结合三河马的特殊血缘关系分析出：除去普氏野马共有3个分支，巴尔虎马与乌审马合并为一个分支，具有较近的遗传关系;乌珠穆沁马与锡尼河马具有较近的遗传关系;三河马与其他4种蒙古马独立地分出来构成一个系统发育分支。通过对蒙古马的分析表明，采用最大似然法结合SDBP值分析生物的系统拓扑关系是较最大似然法结合AU值以及结合BP值分析生物的系统拓扑关系更有效的方法。同时也比UPGMA法有效，而与NJ法得到相似的结论。