基于高速雙倍自助法對蒙古馬系統發育樹的研究

任愛珍，白東義，趙一萍，侯 娜，芒 來*

(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法得到相似的結論。