畜牧兽医学报  2018, Vol. 49 Issue (9): 1861-1869. DOI: 10.11843/j.issn.0366-6964.2018.09.007    PDF    
基于高速双倍自助法对蒙古马系统发育树的研究
任爱珍1, 白东义2, 赵一萍2, 侯娜1, 芒来2     
1. 内蒙古农业大学理学院, 呼和浩特 010018;
2. 内蒙古农业大学动物科学学院 内蒙古自治区蒙古马遗传资源保护及马产业工程实验室, 呼和浩特 010018
摘要:旨在利用高速双倍自助法研究蒙古马系统发育树的问题。本研究选用30匹蒙古马的mtDNA D-loop区碱基序列,其中包括5个蒙古马类群(巴尔虎马、三河马、乌审马、乌珠穆沁马、锡尼河马)各6匹,对其所有105个可能系统树(普氏野马作为外群,6匹)进行分析,利用生物信息学软件Mega6和PAMAL4.9,以最大似然法估计蒙古马的最大似然系统树。最后利用生物信息软件CONSEL、R3.0中的SDBP包等计算105个候补系统树的基于高速双倍自助法的可靠度(speedy double bootstrap P-value,SDBP)。结果表明,SDBP值最大的蒙古马系统树共有3个分支,巴尔虎马与乌审马合并为一个分支,具有较近的遗传关系;乌珠穆沁马与锡尼河马具有较近的遗传关系;三河马独立于其他4种蒙古马,构成一个系统发育分支。SDBP值最大的蒙古马系统树与三河马和其他4类蒙古马具有较远的预测血缘关系是一致的,同时此结果的系统树也具有最大似然值。本研究结果还表明,采用最大似然法结合SDBP值分析蒙古马的系统拓扑关系是较最大似然法结合渐近无偏(approximately unbiased,AU)值以及结合自助(bootstrap P-value,BP)值分析蒙古马系统拓扑关系更有效。同时也比非加权组平均法(unweighted pair-group method with arithmetic means,UPGMA)有效,而与邻接法(neighbor joining,NJ)得到了相近的结论。上述结果提示,我们的研究结果使蒙古马的进化理论奠定在更坚实的遗传学基础上,可推动蒙古马的进化理论进一步发展与完善。
关键词蒙古马    普氏野马    高速双倍自助法    系统进化树    
Phylogenetic Tree Analysis of Mongolian Horses Using the Speedy Double Bootstrap Method
REN Ai-zhen1, BAI Dong-yi2, ZHAO Yi-ping2, HOU Na1, DUGARJAVIIN Mang-lai2     
1. College of Science, Inner Mongolia Agricultural University, Hohhot 010018, China;
2. Inner Mongolia Autonomous Region Mongolian Horse Genetic Resources Protection and Horse Industry Engineering Laboratory, College of Animal Science, Inner Mongolia Agricultural University, Hohhot 010018, China
Abstract: The aim of this study was to analyze the phylogenetic tree of Mongolian horses using the speedy double bootstrap method. The mitochondrial DNA D-loop base sequences of 30 modern horses were used, which included 5 Mongolian horse subspecies (Baerhu, Sanhe, Wushen, Ujimqin and Xinihe horses) 6 for each subspecies and 6 Przewalski's horse as the outgroup, and all 105 possible phylogenetic trees were analyzed. Mega6 and PAMAL4.9 were used, and maximum likelihood phylogenetic tree of Mongolian horse were estimated by maximum likelihood method. Finally, the 3rd order accurate P-value (speedy double bootstrap P-value, SDBP) of the 105 phylogenetic trees was calculated with the CONSEL, SDBP packages in R3.0. The phylogenetic tree of Mongolian horse with the largest SDBP value had 3 branches; Baerhu and Wushen horses were clustered together, Ujimqin and Xinihe horses were clustered together, and Sanhe horse formed an independent branch which was sister to all other Mongolian subspecies. The Mongolian horse phylogenetic tree with the largest SDBP value was consistent with the predicted relationship between Sanhe and the other 4 Mongolian horse subspecies, and this topological tree also had the highest likelihood value. The results of this study also showed that it was more effective to analyze the topological relationships of Mongolian horses using the maximum likelihood method combined with the SDBP value than using the maximum likelihood method combined with the approximately unbiased (AU) value and the bootstrap P-value (BP) value. Additionally, this approach was more effective than the unweighted pair-group method with arithmetic means (UPGMA), and had similar results to the neighbor joining (NJ) method. The results of this study provided robust genetic data which would produce a more complete picture of Mongolian horse evolution.
Key words: Mongolian horses     Przewalski's horse     speedy double bootstrap method     phylogenetic tree    

各种分子系统进化树的构建方法只给出了一个真实进化树的点估计,需要对其可靠性给出概率描述[1]。目前提出的常用的基于评价最大似然树的可靠性方法有:自助法、下平-长谷川法(shimodaira-hasegawa method,SH)、多尺度自助法(multiscale bootstrap method, MBP)、双倍自助法(double bootstrap method,DBP)、高速双倍自助法(speedy double boostrap method,SDBM)[2-5],其中利用自助法计算出的P-值用自助P-值(bootstrap P-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 boostrap P-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 蒙古马D-loop区碱基序列的GenBank登录号 Table 1 GenBank accession numbers of Mongolian horse D-loop sequences
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棵系统树的概率函数:

${f_k}\left( {x, {\theta _k}} \right) $ (1)

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

表 2 S种生物DNA碱基序列 Table 2 DNA base sequences of S species

根据以往的研究报道设定碱基置换模型中的参数,所以碱基置换模型是已知的,只有树枝长度是未知。本研究中的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以最大似然估计法计算每棵系统树的对数似然函数:

$ {l_k} = {l_k}\left( {{\theta _k};{X_{S \times n}}} \right) = \sum\limits_{i = 1}^n {\log {f_k}\left( {{x_i}, {\theta _k}} \right)} $ (2)

其中,k表示第k棵系统树,θk表示树枝长度向量,xi表示DNA碱基序列中的第i位碱基列,Xs×n表示S种生物DNA碱基序列矩阵,fk(x; θk)是第k棵系统树的概率函数。利用软件PAML4.9计算出的对数似然函数可估计出树枝长度向量θk的最大似然估计${\hat \theta _k}$

$ {{\hat \theta }_k} = \mathop {\arg \;\max }\limits_{{\theta _k} \in {{\mathit{\Theta}} _k}} \sum\limits_{i = 1}^n {\log {f_k}\left( {{x_i}, {\theta _k}} \right)}, k = 1, \ldots, K $ (3)

其中,kθkxi所表示的含义与公式(2)相同,Θk表示第k棵树的树枝长度向量θk所在的向量空间,K表示所有可能的候补树。把${\hat \theta _k}$代入lk=lk(θk; Xs×n)=$\sum\limits_{i = 1}^n {\log {f_k}\left({{x_i}, {\theta _k}} \right)} $中得到每棵系统树的对数似然估计值${{\hat l}_k} = {{\hat l}_k}\left({{{\hat \theta }_k}; {X_{s \times n}}} \right)$,简记为${{\hat l}_k}$

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

$ \left[{\begin{array}{*{20}{c}} {\log {f_1}\left( {{x_1};{{\hat \theta }_1}} \right)}& \cdots &{\log {f_1}\left( {{x_n};{{\hat \theta }_1}} \right)}\\ \vdots &{}& \vdots \\ {\log {f_K}\left( {{x_1};{{\hat \theta }_K}} \right)}& \cdots &{\log {f_K}\left( {{x_n};{{\hat \theta }_K}} \right)} \end{array}} \right] $ (4)

简记为

$ \left[{\begin{array}{*{20}{c}} {{{\hat l}_{11}}}& \cdots &{{{\hat l}_{1n}}}\\ \vdots &{}& \vdots \\ {{{\hat l}_{K1}}}& \cdots &{{{\hat l}_{Kn}}} \end{array}} \right] $

其中,${{\hat l}_{kj}} = {\rm{log}}{f_k}({x_j}; {{\hat \theta }_k}), k = 1, \ldots, K; j = 1, \ldots, n$是第k棵树的第j位点序列的最大对数似然估计值。

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

$\hat l = \left( \begin{array}{l} {{\hat l}_1}\\ \vdots \\ {{\hat l}_K} \end{array} \right) $ (5)

其中,分量${{\hat l}_k}$, k=1, …, K是公式(3)中的${{\hat \theta }_k}$代入lk=lk(θk; Xs×n)=$\sum\limits_{i = 1}^n {\log {f_k}\left({{x_i}, {\theta _k}} \right)} $中得到的${{\hat l}_k}$。然后记${{\hat l}_k} = {\rm{max}}\{ {{\hat l}_1}, {{\hat l}_2}, \ldots, {{\hat l}_K}\} $为所有系统树最大对数似然值${{\hat l}_k}$, k=1, …, K中的最大数值。

步骤4:估计出最大似然树。由${{\hat l}_k} = {\rm{max}}\{ {{\hat l}_1}, {{\hat l}_2}, \ldots, {{\hat l}_K}\} $得到最大似然树编号$\hat k = \arg \; \max \{ {{\hat l}_1}, {{\hat l}_2}, \ldots, {{\hat l}_K}\} $

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

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

${H_1}:{\mu _1} = \mathop {\max }\limits_{i = 1, \cdots, K} {\mu _i}$H1表示第1棵系统树为真实系统树;

${H_K}:{\mu _K} = \mathop {\max }\limits_{i = 1, \cdots, K} {\mu _i}$HK表示第K棵系统树为真实系统树。

其中,${\mu _{_K}} = {E_q}\left[{{l_k}\left({{{\hat \theta }_k}; X{\prime _{S \times n}}} \right)} \right]$k=1, …, K表示第k棵系统树对数似然值的期望值,这里XS×n设为一组未知的DNA碱基序列;${\mu _{_K}} = {E_q}\left[{{l_k}\left({{{\hat \theta }_k}; X{\prime _{S \times n}}} \right)} \right]$中的q表示q(x),q(x)是每个位点碱基列x的真实分布,${{l_k}\left({{{\hat \theta }_k}; X{\prime _{S \times n}}} \right)}$表示第k棵系统树的对数似然值,它与利用公式(3)估计出的${{\hat l}_k} = {l_k}\left({{{\hat \theta }_k}; {X_{S \times n}}} \right)$是不同的未知对数似然值,${{\hat \theta }_k}$是公式(3)中的已知向量。由于碱基位点序列x的真实分布q(x)是未知的,所以${{l_k}\left({{{\hat \theta }_k}; X{\prime _{S \times n}}} \right)}$${\mu _{_K}} = {E_q}\left[{{l_k}\left({{{\hat \theta }_k}; X{\prime _{S \times n}}} \right)} \right]$, k=1,…, K也都是未知的。以H1为例说明计算SDBP的值,计算其他假设的SDBP值的步骤类似。

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

${H_1}:{\mu _1} = \mathop {\max }\limits_{i = 1, \cdots, K} {\mu _i}$为原假设,即认为第1棵系统树为真实系统树,否则为备择假设。

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

(1) 计算每棵系统树的最大对数似然值${{\hat l}_k}$, k=2, …,K与第1棵系统树${{\hat l}_1}$的差值,然后取最大值,即符号距离d

$ d = \mathop {\max }\limits_{j = 2, \cdots, K} {{\hat l}_j}{\rm{-}}{{\hat l}_1} $ (6)

(2) 重抽样

重抽样的样本取自K元正态总体NK($\hat \mu, \Sigma $),即l*~NK($\hat \mu, \Sigma $)。其中$\hat \mu = \left({\hat \mu, \cdots {{\hat \mu }_K}} \right)$表示公式(5)中的最大对数似然向量${\hat l}$在假设H1所表示的领域G的边界上的射影,$\hat \mu = \left({{{\hat \mu }_1}, \cdots {{\hat \mu }_K}} \right)$的计算可用PAVA (pool adjacent violators algorithm method)[25-26]来计算:

$ {{\hat \mu }_1} = \frac{{{\Sigma _{j \in W}}{{\hat l}_j}}}{{\# W}}, {{\hat \mu }_j} = \min \left( {{{\hat \mu }_1}, {{\hat \mu }_j}} \right), j = 2, \cdots, K $ (7)

其中,集合W是集合{1, …,K}的子集,并且包含元素1,符号#W表示集合W所含元素的个数。另一方面,最大对数似然向量${\hat l}$的协方差矩阵$\Sigma $=(σij)的元素σij,可按如下的公式计算[27]

$ \begin{array}{l} V\left( {{{\hat l}_i}} \right) = n \times \frac{1}{{n- 1}}\sum\limits_{h = 1}^n {[\log {f_i}\left( {{x_h};{{\hat \theta }_i}} \right)}-\\ \;\;\;\;\;\;\;\;\;\;\;\frac{1}{n}\sum\limits_{h = 1}^n {\log {f_i}\left( {{x_h};{{\hat \theta }_i}} \right)} {]^2} \end{array} $ (8)
$ \begin{array}{l} Cov\left( {{{\hat l}_i}, {{\hat l}_j}} \right) = n \times \frac{1}{{n- 1}}\sum\limits_{h = 1}^n {\{ [\log {f_i}\left( {{x_h};{{\hat \theta }_i}} \right)}-\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{1}{n}\sum\limits_{h = 1}^n {\log {f_i}\left( {{x_h};{{\hat \theta }_i}} \right)}] \times \\ \left[{\log {f_j}\left( {{x_h};{{\hat \theta }_j}} \right)-\frac{1}{n}\sum\limits_{h = 1}^n {\log {f_j}\left( {{x_h};{{\hat \theta }_i}} \right)} } \right]\} \end{array} $ (9)

从总体NK($\hat \mu, \Sigma $)产生B(一般设为1 000或10 000)个自助样本向量:

$ l_1^* = \left( \begin{array}{l} l_{11}^*\\ \vdots \\ l_{K1}^* \end{array} \right), \cdots, l_B^* = \left( \begin{array}{l} l_{1B}^*\\ \vdots \\ l_{KB}^* \end{array} \right) $ (10)

(3) 计算SDBP值

对每个自助样本lb*, b=1,…,B按公式(6)计算符号距离db*, b=1, …,B,最终得到SDBP计算公式:

$ SDBP = \frac{{\# \left\{ {d_b^* > d} \right\}}}{B}, b = 1, \cdots, B $ (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)值以及结合自助(bootstrap P-value,BP)值分析蒙古马系统拓扑关系更有效。将t100号、t66号以及t19号系统树的拓扑结构画成二分系统树的树形,得到图 1图 2图 3

表 3 30匹蒙古马的60棵候选树SDBP、AU以及BP值的计算结果 Table 3 The results of SDBP, AU and BP values of 60 candidate phylogenetic trees for 30 Mongolian horses
图 1 SDBP值最大系统树t100拓扑结构 Figure 1 The topology of tree t100 with the highest SDBP (speedy double bootstrap P-value) value
图 2 AU值最大系统树t66拓扑结构 Figure 2 The topology of tree t66 with the highest AU (approximately unbiased) value
图 3 BP值最大系统树t19拓扑结构 Figure 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树拓扑结构 Figure 4 The topology of NJ trees of 5 Mongolian horse subspecies
图 5 5个蒙古马类群的UPGMA树拓扑结构 Figure 5 The topology of UPGMA tree of 5 Mongolian horse subspecies

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

3 讨论

从统计学假设检验理论来判断,每棵系统发育树的SDBP值相当于原假设HkP-值,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法得到的结论相似。

另外,当种与种之间分类数或种内部的分类数超过6类(包括外群),例如是7类时,利用最大似然法估计最大似然树时,所有候补树的个数为945个。由于树个数的增加,导致最大似然法估计的计算量也增加。所以当分类数超过6类时基于最大似然法的AU和SDBP方法的计算量也很大,这也是AU和SDBP方法的限制所在。而BP方法除了在最大似然法的应用外,也可以在NJ法以及UPGMA法中使用。另一方面,SDBP方法中的$\hat \mu = \left({{{\hat \mu }_1}, \cdots {{\hat \mu }_K}} \right)$是最大对数似然向量${\hat l}$在假设H1所表示的领域G的边界上的射影,$\hat \mu = \left({{{\hat \mu }_1}, \cdots {{\hat \mu }_K}} \right)$的计算在其他系统树构树方法如NJ法与UPGMA中很难估计。在这些情况下, 更适合用AU或BP计算系统树的可靠度。

4 结论

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

参考文献
[1] 吕宝忠. 分子进化树的构建[J]. 动物学研究, 1993, 14(2): 186–193.
LÜ B Z. The construction of molecular evolutionary trees[J]. Zoological Research, 1993, 14(2): 186–193. (in Chinese)
[2] FELSENSTEIN J. Confidence limits on phylogenies:an approach using the bootstrap[J]. Evolution, 1985, 39(4): 783–791. DOI: 10.1111/j.1558-5646.1985.tb00420.x
[3] SHIMODAIRA H, HASEGAWA M. Multiple comparisons of log-likelihoods with applications to phylogenetic inference[J]. Mol Biol Evol, 1999, 16(8): 1114–1116. DOI: 10.1093/oxfordjournals.molbev.a026201
[4] SHIMODAIRA H. An approximately unbiased test of phylogenetic tree selection[J]. Syst Biol, 2002, 51(3): 492–508. DOI: 10.1080/10635150290069913
[5] REN A Z, ISHIDA T, AKIYAMA Y. Assessing statistical reliability of phylogenetic trees via a speedy double bootstrap method[J]. Mol Phylogenet Evol, 2013, 67(2): 429–435. DOI: 10.1016/j.ympev.2013.02.011
[6] 杨子恒. 分子进化树的统计推断[J]. 遗传, 1995, 17(S1): 92–96.
YANG Z H. Statistical estimation of molecular evolutionary trees[J]. Hereditas (Beijing), 1995, 17(S1): 92–96. (in Chinese)
[7] 杨子恒.计算分子进化[M].钟扬, 译.上海: 复旦大学出版社, 2008: 207-215.
YANG Z H.Computational molecular evolution[M].ZHONG Y, trans.Shanghai: Fudan University Press, 2008: 207-215.(in Chinese)
[8] ZHARKIKH A, LI W H. Statistical properties of bootstrap estimation of phylogenetic variability from nucleotide sequences.Ⅱ.Four taxa without a molecular clock[J]. J Mol Evol, 1992, 35(4): 356–366. DOI: 10.1007/BF00161173
[9] HILLIS D M, BULL J J. An empirical test of bootstrapping as a method for assessing confidence in phylogenetic analysis[J]. Syst Biol, 1993, 42(2): 182–192. DOI: 10.1093/sysbio/42.2.182
[10] XU X F, ÁRNASON Ú. The complete mitochondrial DNA sequence of the horse, Equus caballus:extensive heteroplasmy of the control region[J]. Gene, 1994, 148(2): 357–362. DOI: 10.1016/0378-1119(94)90713-7
[11] ISHIDA N, OYUNSUREN T, MASHIMA S, et al. Mitochondrial DNA sequences of various species of the genus Equus with special reference to the phylogenetic relationship between Przewalskii's wild horse and domestic horse[J]. J Mol Evol, 1995, 41(2): 180–188.
[12] KIM K I, YANG Y H, LEE S S, et al. Phylogenetic relationships of Cheju horses to other horse breeds as determined by mtDNA D-loop sequence polymorphism[J]. Anim Genet, 1999, 30(2): 102–108. DOI: 10.1046/j.1365-2052.1999.00419.x
[13] PETERSEN J L, MICKELSON J R, RENDAHL A K, et al. Genome-wide analysis reveals selection for important traits in domestic horse breeds[J]. PLoS Genet, 2013, 9(1): e1003211. DOI: 10.1371/journal.pgen.1003211
[14] SCHUBERT M, JÓNSSON H, CHANG D, et al. Prehistoric genomes reveal the genetic foundation and cost of horse domestication[J]. Proc Natl Acad Sci U S A, 2014, 111(52): E5661–E5669. DOI: 10.1073/pnas.1416991111
[15] 李金莲, 芒来, 石有斐. 利用微卫星标记对蒙古马和纯血马遗传多样性的研究[J]. 畜牧兽医学报, 2005, 36(1): 6–9.
LI J L, MANG L, SHI Y F. Evaluation of genetic diversity of Mongolian horse and thoroughbred horse using microsatellite markers[J]. Acta Veterinaria et Zootechnica Sinica, 2005, 36(1): 6–9. DOI: 10.3321/j.issn:0366-6964.2005.01.002 (in Chinese)
[16] 任秀娟, 赵一萍, 萨如拉, 等. 马属动物进化特征概述[J]. 畜牧兽医学报, 2017, 48(3): 385–392.
REN X J, ZHAO Y P, SARULA, et al. The review for the characteristics of the equus evolution[J]. Acta Veterinaria et Zootechnica Sinica, 2017, 48(3): 385–392. (in Chinese)
[17] 赵启南, 芒来, 白东义, 等. 蒙古马高负荷运动训练前后转录组差异表达分析[J]. 畜牧兽医学报, 2017, 48(6): 1007–1016.
ZHAO Q N, MANG L, BAI D Y, et al. Mongolian horses transcriptome differential expression analysis before and after a high load exercise training[J]. Acta Veterinaria et Zootechnica Sinica, 2017, 48(6): 1007–1016. (in Chinese)
[18] 芒来, 杨虹. 蒙古马遗传多样性研究进展[J]. 遗传, 2008, 30(3): 269–276.
MANG L, YANG H. Progress in the study of genetic diversity of Mongolian horse[J]. Hereditas (Beijing), 2008, 30(3): 269–276. DOI: 10.3321/j.issn:0253-9772.2008.03.004 (in Chinese)
[19] 李金莲.中国蒙古马遗传多样性和分子系统进化研究[D].呼和浩特: 内蒙古农业大学, 2006.
LI J L.Study on genetic diversity and molecular phylogeny evolution in Chinese Mongolian horse[D].Hohhot: Inner Mongolia Agricultural University, 2006.(in Chinese) http://cdmd.cnki.com.cn/Article/CDMD-10129-2006101295.htm
[20] 陈建兴, 孙玉江, 乌尼尔夫, 等.基于ND4基因分析家马的母系起源[J].黑龙江畜牧兽医, 2013(3): 8-10, 14.
CHEN J X, SUN Y J, WUNIERFU, 等.Analysis of the maternal origins of domestic horses based on ND4 gene[J].Heilongjiang Animal Science and Veterinary Medicine, 2013(3): 8-10, 14.(in Chinese) http://kns.cnki.net/KCMS/detail/detail.aspx?filename=HLJX201303004&dbname=CJFD&dbcode=CJFQ
[21] 王小斌, 秦芳, 张云生, 等. 蒙古马mtDNA D-loop区遗传多样性与多重母系起源[J]. 西北农林科技大学学报:自然科学版, 2010, 38(3): 47–51.
WANG X B, QIN F, ZHANG Y S, et al. Mitochondrial DNA D-loop genetic diversity and multiple maternal origins in Mongolian horses[J]. Journal of Northwest A & F University:Natural Science Edition, 2010, 38(3): 47–51. (in Chinese)
[22] FUKAM-KOBAYASHI K, TATENO Y. Robustness of maximum likelihood tree estimation against different patterns of base substitutions[J]. J Mol Evol, 1991, 32(1): 79–91. DOI: 10.1007/BF02099932
[23] HASEGAWA M, KISHINO H, SAITOU N. On the maximum likelihood method in molecular phylogenetics[J]. J Mol Evol, 1991, 32(5): 443–445. DOI: 10.1007/BF02101285
[24] YANG Z H. JSTOR:systematic biology[J]. Syst Biol, 1994, 43(3): 329–342. DOI: 10.1093/sysbio/43.3.329
[25] AYER M, BRUNK H D, EWING G M, et al. An empirical distribution function for sampling with incomplete information[J]. Ann Math Stat, 1955, 26(4): 641–647. DOI: 10.1214/aoms/1177728423
[26] ZHAO H B. Comparing several treatments with a control[J]. J Stat Plan Inference, 2007, 137(9): 2996–3006. DOI: 10.1016/j.jspi.2006.12.005
[27] KISHINO H, MIYATA T, HASEGAWA M. Maximum likelihood inference of protein phylogeny and the origin of chloroplasts[J]. J Mol Evol, 1990, 31(2): 151–160. DOI: 10.1007/BF02109483
[28] 陈建兴.中国蒙古马的遗传多样性与系统发育及起源研究[D].呼和浩特: 内蒙古农业大学, 2012.
CHEN J X.Investigation of the genetic diversity, phylogeny and origins of Chinese Mongolian horses[D]. Hohhot: Inner Mongolia Agricultural University, 2012.(in Chinese) http://cdmd.cnki.com.cn/Article/CDMD-10129-1013153640.htm