植物血紅蛋白基因家族的分子進(jìn)化研究

齊小瓊+王艇

摘要： 植物血紅蛋白 （Hemoglobin） 是一類由珠蛋白 （Globin） 和血紅素 （Ferroheme） 組成的結(jié)合蛋白，廣泛存在于植物界中。根據(jù)序列特征、表達(dá)模式及配基結(jié)合性質(zhì)可將血紅蛋白分為共生血紅蛋白 （Symbiotic hemoglobin， sHb）、非共生血紅蛋白 （Nonsymbiotic hemoglobin， NsHb） 和截短的血紅蛋白 （Truncated hemoglobin， tHb） 3大類，并對(duì)植物血紅蛋白基因家族進(jìn)行了分子系統(tǒng)發(fā)育分析和進(jìn)化研究。結(jié)果表明，通過(guò)基因重復(fù)在被子植物中產(chǎn)生了class 1和class 2兩類血紅蛋白基因。該基因重復(fù)發(fā)生在單、雙子葉植物分化之前，但由于class 2類基因在單子葉植物中發(fā)生了一次基因丟失事件，致使class 1類基因廣泛分布于單、雙子葉植物中，而class 2類基因則僅見于雙子葉植物；class 2類基因的進(jìn)化速率高于class 1類基因，并分別在class 2、sHbs的祖先支上檢測(cè)到了正選擇的發(fā)生，表明基因重復(fù)發(fā)生后，class 2、sHbs通過(guò)插曲式進(jìn)化和適應(yīng)性進(jìn)化獲得了不同的功能；在class 1和class 2的功能性分歧位點(diǎn)檢測(cè)中獲得14個(gè)I類功能性分歧位點(diǎn)，但未能發(fā)現(xiàn)II類功能性分歧位點(diǎn)，這些功能性分歧位點(diǎn)是由正選擇作用或者放松的選擇限制產(chǎn)生的。該研究結(jié)果可以指導(dǎo)后續(xù)的試驗(yàn)工作，來(lái)驗(yàn)證那些具有重要功能的氨基酸置換位點(diǎn)和正選擇位點(diǎn)，從而揭示功能性分歧和適應(yīng)性進(jìn)化的分子基礎(chǔ)。

關(guān)鍵詞：血紅蛋白基因；基因家族；分子進(jìn)化分析；適應(yīng)性進(jìn)化；功能性分歧

中圖分類號(hào)：Q37 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：0439-8114（2014）18-4447-10

根據(jù)序列組成、表達(dá)模式以及配體結(jié)合性質(zhì)等特點(diǎn)，植物血紅蛋白可分為共生血紅蛋白（Symbiotic hemoglobin， sHb）、非共生血紅蛋白 （Nonsymbiotic hemoglobin， NsHb） 和截短的血紅蛋白（Truncated hemoglobin， tHb） 3大類。其中，根據(jù)與氧的親和力和序列的相似性，NsHb又可分為NsHb1和NsHb2兩大類。系統(tǒng)發(fā)育分析表明，植物血紅蛋白通過(guò)基因重復(fù)事件在被子植物中產(chǎn)生了class 1和 class 2兩大類基因[1-5]，其中class 1類由NsHb1組成，class 2類由NsHb2和sHb組成；而class 1和class 2中均不包含tHb，tHb是通過(guò)基因水平轉(zhuǎn)移從細(xì)菌中獲得[6]，與植物中的NsHb和sHb形成了不同的進(jìn)化譜系。

研究表明，class 1類基因廣泛存在于被子植物中，而class 2類僅見于雙子葉植物中，迄今為止沒(méi)有在單子葉植物中發(fā)現(xiàn)。因此，對(duì)于class 2類基因在植物中的分布范圍還沒(méi)有確定——該類基因在單子葉植物中是缺失的還是普遍存在或是僅存在于部分支系中由于數(shù)據(jù)較少而沒(méi)有檢測(cè)到。目前有兩種可能的解釋：①基因重復(fù)在單、雙子葉分化之前就已經(jīng)發(fā)生，但是class 2類基因在單子葉植物中丟失了；②基因重復(fù)只在雙子葉植物中發(fā)生，即發(fā)生了雙子葉植物特異性的重復(fù)事件，該過(guò)程不涉及單子葉植物[7]。顯然第二種假說(shuō)更簡(jiǎn)約，但要驗(yàn)證這兩個(gè)假說(shuō)，需要更多單子葉植物基因組的序列數(shù)據(jù)及對(duì)基因重復(fù)和單、雙子葉分歧發(fā)生時(shí)間的正確估算、綜合分析判定。針對(duì)基因家族成員在基因重復(fù)之后的功能分化，Ohno[8]提出了重復(fù)-保留-非/新功能化 （Duplication-Retention-Non/Neofunctionlization， DRNNF） 模型，Hughes[9]提出了蛋白質(zhì)亞功能化 （Subfunctionalization， SF）模型，F(xiàn)orce等[10]提出了DDC模型。對(duì)于植物血紅蛋白基因的兩個(gè)拷貝，均在進(jìn)化過(guò)程中保留了下來(lái)且具有了不同的功能，發(fā)生了功能性分歧。針對(duì)基因重復(fù)發(fā)生后，class 1、class 2類基因發(fā)生的功能性分歧，有兩個(gè)不同的假說(shuō)可以解釋：①基因重復(fù)之后，一個(gè)重復(fù)基因的功能被選擇性的保留，另一個(gè)基因逃脫了選擇限制而獲得了新的功能[8-12]；②基因重復(fù)發(fā)生之后，一個(gè)重復(fù)基因通過(guò)適應(yīng)性進(jìn)化獲取了新的功能[13]。研究表明，正選擇僅發(fā)生在特定時(shí)間、個(gè)別支系、蛋白質(zhì)的某些位點(diǎn)，因此特定的正選擇信號(hào)很容易被大量的負(fù)選擇及中性選擇信號(hào)所淹沒(méi)[14]，如果檢測(cè)方法不當(dāng)，很難檢測(cè)到。

隨著研究的深入，越來(lái)越多物種的血紅蛋白序列被測(cè)定，尤其是對(duì)古老陸生植物中血紅蛋白基因的研究，推斷出血紅蛋白祖先基因的序列特征、表達(dá)模式、三級(jí)結(jié)構(gòu)、配體結(jié)合特征及進(jìn)化樣式等[15]，為該基因家族進(jìn)化研究奠定了基礎(chǔ)。但是對(duì)于被子植物血紅蛋白基因重復(fù)發(fā)生所涉及的物種范圍，尤其是class 2類基因在植物中的存在狀況及共生血紅蛋白基因的來(lái)源問(wèn)題等都沒(méi)有很好的解答，還需選取更多代表被子植物各分類群的血紅蛋白基因序列，通過(guò)分析血紅蛋白基因在這些物種中的存在狀況、序列組成、基因結(jié)構(gòu)等特征，闡明該基因家族在不同分類群中的進(jìn)化特征和進(jìn)化趨勢(shì)，并進(jìn)一步通過(guò)分子進(jìn)化分析的各種方法從時(shí)間尺度和空間尺度闡釋該家族基因在陸生植物中的起源和進(jìn)化樣式。

本研究從GenBank數(shù)據(jù)庫(kù)中選取盡可能多的被子植物NsHb和sHb基因序列，通過(guò)密碼子使用分析、相對(duì)進(jìn)化速率檢測(cè)、適應(yīng)性進(jìn)化分析、基因重復(fù)發(fā)生時(shí)間的估算、功能性分歧分析等，來(lái)初步確定基因重復(fù)涉及的物種范圍；闡明在該基因重復(fù)早期是放松的選擇限制還是正選擇作用起主導(dǎo)作用，同時(shí)進(jìn)一步鑒定發(fā)生選擇限制改變或者氨基酸性質(zhì)改變的功能性分歧位點(diǎn)；闡明植物血紅蛋白的系統(tǒng)發(fā)育關(guān)系和進(jìn)化樣式及基因重復(fù)發(fā)生之后的選擇壓力變化情況，并且進(jìn)一步鑒定出正選擇位點(diǎn)。

1 材料與方法

1.1 序列數(shù)據(jù)的收集和選取

用擬南芥（Arabidopsis thaliana）的AHB1和AHB2基因（登錄號(hào)：AAD26949， AAB82770） 的蛋白序列作為查詢序列，從GenBank數(shù)據(jù)庫(kù)中搜索盡可能多的被子植物NsHb和sHb基因序列。endprint

1.2 序列比對(duì)和系統(tǒng)發(fā)育分析

將血紅蛋白核苷酸序列翻譯成氨基酸序列，經(jīng)MEGA 4.0軟件[16，17]對(duì)位排列，再根據(jù)氨基酸對(duì)位排列結(jié)果進(jìn)行核苷酸序列對(duì)位排列。對(duì)位排列的結(jié)果根據(jù)需要進(jìn)行手工校正，將分歧度很大的區(qū)域或者空格等從比對(duì)中移除，以免產(chǎn)生比對(duì)的不確定性。采用ModelTest軟件[18]確定GTR+I+G為最適的核苷酸置換模型；利用PHYML[19]軟件，根據(jù)最大似然法（Maximum likelihood， ML）構(gòu)建系統(tǒng)發(fā)育樹，并進(jìn)行了1 000次的Bootstrap檢驗(yàn)。

1.3 同義密碼子使用偏好性分析

采用CodonW程序（http：//bioweb.pasteur.fr/seqanal/interfaces/codonw.html）計(jì)算有效密碼子數(shù) （ENC）[20]，最優(yōu)密碼子頻率（FOP）[21]，密碼子適應(yīng)指數(shù)（CAI）和同義密碼子第三位上的G+C含量 （GC3s），并且檢測(cè)同義密碼子使用偏好性在每個(gè)基因間的差異。

1.4 相對(duì)進(jìn)化速率檢測(cè)

采用RRTree軟件[22]計(jì)算基因重復(fù)形成的不同亞家族在進(jìn)化過(guò)程中的置換速率，并比較它們之間的差異，計(jì)算時(shí)以苔蘚植物小立碗蘚 （Physcomitrella patens）和角齒蘚（Ceratodon purpureus）中的血紅蛋白基因作為外類群。

1.5 基因重復(fù)的時(shí)間估計(jì)

根據(jù)系統(tǒng)發(fā)育構(gòu)建的ML樹將所研究的植物分為class 1和class 2兩個(gè)類群。采用BEAST軟件包中自帶的Uncorrelated Lognormal Distributed （UCLD） Relaxed Clock模型估算ML樹上每個(gè)節(jié)點(diǎn)的分歧時(shí)間[23]。主要分支的分歧時(shí)間根據(jù)最近共同祖先的時(shí)間 （tMRCA） 值來(lái)進(jìn)行校正。最適的核苷酸置換模型采用Modeltest確定的GTR+I+G模型，總共迭代計(jì)算60 000 000代，每1 000代保存1株樣本；最后選用Tracer v1.4.1軟件[24]檢測(cè)收斂程度。計(jì)算時(shí)把最初獲得的15 000株老化樣本去除，采用剩余樣本重建時(shí)間尺度下的發(fā)育樹。

1.6 適應(yīng)性進(jìn)化分析

將系統(tǒng)發(fā)育分析得到的ML樹分為3個(gè)分析組——1個(gè)大數(shù)據(jù)組 （全部58個(gè)序列組成） 和兩個(gè)小數(shù)據(jù)組 （分別僅含class 1或class 2中的全部序列），采用PAML 4分析方法中的3類模型：分支模型[25]、位點(diǎn)模型[24，26]以及分支-位點(diǎn)模型[14，27]進(jìn)行分析，序列位點(diǎn)處理采用cleandata=0進(jìn)行設(shè)置。

在分支模型中，分別進(jìn)行了單比率（One-ratio） 和二比率（Two-ratio）模型分析。在二比率模型設(shè)置中，為了檢測(cè)class 1和class 2是否存在差異，假定class 2有著不同的ω值，即在大數(shù)據(jù)組中設(shè)定：（Outgroup#ω0，NsHb1#ω0，NsHb2#ω1，sHb#ω1，#ω1，#ω0）；為了檢測(cè)共生的血紅蛋白與非共生血紅蛋白的差異，假定sHbs有著不同的ω值，即在僅含有class 2的小數(shù)據(jù)組中設(shè)定：（NsHb2#ω0， sHb#ω1）。分別將這些二比率模型同其相對(duì)應(yīng)的單比率模型比較，進(jìn)行似然比檢驗(yàn)（LRT），分析所得結(jié)果。

位點(diǎn)模型分析時(shí)將class 1和class 2單獨(dú)組成的兩個(gè)小數(shù)據(jù)組分別進(jìn)行，采用了三對(duì)模型M1a （近中性模型）和M2a（正選擇模型），M0（單比率模型）和M3（離散模型），M7（beta 模型）和M8（beta & ω模型），并進(jìn)行似然比檢驗(yàn)。模型M8a（beta &ω=1）同M8類似，只是將ω值固定為1，比較模型M8和M8a可以大大降低假陽(yáng)性。

分支-位點(diǎn)模型可以檢測(cè)出在特定支系上的特定位點(diǎn)上存在的正選擇。本研究采用該模型分別檢測(cè)在class 2和sHbs的祖先分支上的正選擇位點(diǎn)，因此運(yùn)用大數(shù)據(jù)組將class 2的祖先支作為前景支，運(yùn)用小數(shù)據(jù)組class 2將sHbs的祖先支作為前景支。model A設(shè)為：model=2，NSsite=2，fix_omega=0；model A1 設(shè)為：model=2，NSsite=2，fix_omega=1，omega =1。將model A和model A1進(jìn)行LRT。

1.7 功能性分歧分析

為了研究基因重復(fù)發(fā)生之后各重復(fù)基因的功能性分歧及結(jié)構(gòu)差異，采用Gu[28，29]提出的方法計(jì)算class 1和class 2之間的兩類功能性分歧系數(shù)（θI和θII），通過(guò)DIVERGE v2.0軟件包[30]實(shí)現(xiàn)。第一種類型（type I）的功能性分歧氨基酸位點(diǎn)具有如下特點(diǎn)：它們?cè)谝粋€(gè)亞組中非常保守，而在另一個(gè)亞組中高度可變，反之亦然[28，31-33]?？赡馨凳具@些殘基經(jīng)歷了改變的功能限制而導(dǎo)致了位點(diǎn)特異的速率差異；而第二種類型（type II） 的功能性分歧位點(diǎn)分別在兩個(gè)亞組中都非常保守，但是它們的生物化學(xué)性質(zhì)差別很大，因此導(dǎo)致了成簇特異的氨基酸性質(zhì)的轉(zhuǎn)換（Shift of cluster-specific amino acid property），例如正電比負(fù)電，暗示著這些殘基可能有助于使不同亞組的功能實(shí)現(xiàn)特化[29]。如果通過(guò)計(jì)算得到θI>0或者θII>0，并且具有統(tǒng)計(jì)顯著性，那么就說(shuō)明基因重復(fù)之后發(fā)生了位點(diǎn)特異的選擇限制的改變或者氨基酸性質(zhì)的巨大轉(zhuǎn)變。

2 結(jié)果與分析

2.1 序列相似性搜索和多序列比對(duì)結(jié)果

通過(guò)序列相似性搜索，從被子植物中獲得了48個(gè)物種58條血紅蛋白序列。其中，有56條是完整的編碼區(qū)序列，還有2條是部分編碼區(qū)序列。從表1中可以看出，大部分蛋白質(zhì)序列都是由150個(gè)左右的氨基酸所組成；但蓖麻（Ricinus communis） 的NsHb比其他序列要長(zhǎng)很多，含有537個(gè)氨基酸，其也包含保守的血紅素結(jié)合位點(diǎn)（Heme-binding site），但是在5端含有一段序列，目前對(duì)于這段延長(zhǎng)序列的作用還有待研究。endprint

2.2 被子植物血紅蛋白基因家族系統(tǒng)發(fā)生關(guān)系和進(jìn)化分析結(jié)果

獲得的58條植物血紅蛋白序列，采用ML法構(gòu)建了系統(tǒng)發(fā)育樹（圖1）。采用PHYML v2.4.4程序[19]構(gòu)建時(shí)，選取的模型及參數(shù)如下：核酸置換模型GTR+I+G、核酸位點(diǎn)分為6類、符合離散的伽馬分布、伽馬形狀參數(shù)為1.341 4、不變位點(diǎn)的比例為0.045 2。用苔蘚植物中的兩條非共生血紅蛋白序列作為外類群，從圖1中可以看出，被子植物的血紅蛋白基因分為兩個(gè)亞家族，定義為class 1和class 2；兩個(gè)分支之間的Bootstrap支持率高達(dá)100%。class 1中既包含雙子葉植物也包含單子葉植物，但在class 2中卻只含有雙子葉植物的基因序列。截至目前，在class 2中未檢測(cè)到有單子葉植物的序列。結(jié)果表明，通過(guò)基因重復(fù)在被子植物中產(chǎn)生了class 1和class 2兩類血紅蛋白基因，該基因重復(fù)可能發(fā)生在單、雙子葉植物分開之前，但后來(lái)class 2類基因在單子葉植物發(fā)生了丟失；也有可能基因重復(fù)僅僅發(fā)生在雙子葉植物中，即發(fā)生了雙子葉植物特異的基因重復(fù)。

2.3 同義密碼子特征

密碼子偏好性也就是同義密碼子間不相等的使用頻率，可能是為了防止蛋白翻譯成本上升而降低氨基酸錯(cuò)誤摻入的可能性[34]。本研究通過(guò)計(jì)算NsHb1、NsHb2和sHb中同義密碼子的偏好性指數(shù)，并比較這些指數(shù)在不同組間的變化，試圖探討核苷酸組成與密碼子使用頻率之間的關(guān)系。其中，F(xiàn)OP和CAI測(cè)度了一個(gè)基因中使用的密碼子與該基因的最適密碼子的偏離程度；ENC則測(cè)度了一個(gè)基因中使用的密碼子與同義密碼子無(wú)偏好使用之間的偏離程度。從表2中可以看出，NsHb1的ENC最低，說(shuō)明其具有更高的同義密碼子使用偏好性；對(duì)比GC3s和ENC之間，發(fā)現(xiàn)兩者呈負(fù)相關(guān)關(guān)系，即ENC越低GC3s越大。GC3s和ENC關(guān)系見圖2，圖2中包含了一條基準(zhǔn)線 （GCref），這條基準(zhǔn)線代表沒(méi)有同義密碼子使用偏性時(shí)的情況[20]。絕大多數(shù)血紅蛋白基因ENC基本沿著基準(zhǔn)線的走向排布，說(shuō)明這些基因的GC3s是ENC的主要決定因素，但是也有個(gè)別基因的ENC偏離了基準(zhǔn)線，表明可能受到了一定的選擇作用。FOP與ENC反映出相同的密碼子使用偏好性，即NsHb1的密碼子偏好性最強(qiáng)；CAI與FOP、ENC分析結(jié)果并不一致，CAI結(jié)果表明NsHb2的密碼子偏好性最強(qiáng)?？赡苁且?yàn)榕c計(jì)算方式有關(guān)，F(xiàn)OP是基于一組物種特異的最適密碼子的使用頻率，ENC是基于每個(gè)氨基酸所使用的密碼子的觀察數(shù)目[35]，而CAI則是基于高表達(dá)基因中的最優(yōu)密碼子表。因此，CAI與基因的表達(dá)情況正相關(guān)，CAI越高，表達(dá)量越高；CAI越低，表達(dá)量越低。

2.4 相對(duì)進(jìn)化速率分析結(jié)果

用苔蘚植物（角齒蘚和小立碗蘚）兩個(gè)血紅蛋白基因作為外類群，計(jì)算并比較不同基因亞家族之間的相對(duì)進(jìn)化速率，結(jié)果見表3。由表3可以看出，NsHb2和sHb亞家族相對(duì)于NsHb1均表現(xiàn)出了非同義置換速率明顯加快的跡象，NsHb2和sHb之間的非同義置換速率沒(méi)有明顯的差異。結(jié)果表明，基因重復(fù)之后不同拷貝之間的進(jìn)化速率發(fā)生了明顯的變化，為后來(lái)的功能性分歧奠定了基礎(chǔ)。

2.5 基因重復(fù)的時(shí)間估計(jì)

使用玉米/水稻 （大約72 Ma）、玉米/小麥 （大約62 Ma） 的分歧時(shí)間[36]作為時(shí)間校對(duì)點(diǎn)，粗略估計(jì)基因重復(fù)發(fā)生的時(shí)間為439.96 Ma。

2.6 功能性分歧分析

本研究用DIVERGE[28]考察了class 1和class 2之間的功能性分歧系數(shù)θ，得到θI=0.592 511 （P<0.01），表明在基因重復(fù)形成的兩個(gè)亞家族中存在置換速率明顯轉(zhuǎn)換的位點(diǎn)；大部分氨基酸位點(diǎn)的后驗(yàn)概率值大于0.5，它們可能在功能性分歧中起主要作用。將后驗(yàn)概率的閾值設(shè)為0.87和0.95，分別在這兩個(gè)水平上檢測(cè)到14個(gè)和7個(gè)改變功能限制的氨基酸位點(diǎn)，但是未能發(fā)現(xiàn)II類功能性分歧位點(diǎn) （圖3）。這些功能性分歧位點(diǎn)是由正選擇作用或放松的選擇限制產(chǎn)生的。

2.7 選擇壓力的檢測(cè)

2.7.1 分支模型 以大數(shù)據(jù)組用ML法構(gòu)建的基因樹為PAML分析的樹文件，進(jìn)行分支模型檢驗(yàn)。單比率模型下估算的ω值為0.198 1（表4）。二比率模型估算的ω1為0.270 3，背景支上的ω0為0.156 2；通過(guò)比較二比率模型和單比率模型，LRT檢驗(yàn)得到2Δl=41.82，P<0.001，df=1，說(shuō)明二比率模型更適合分析的數(shù)據(jù)組，class 2亞家族經(jīng)歷了選擇限制的放松。

以class 2小數(shù)據(jù)組構(gòu)建的基因樹為樹文件進(jìn)行分支模型檢驗(yàn)，單比率模型估算的ω值為0.261 6 （表4）。二比率模型估算的ω1為0.328 7，背景支上的ω0為0.175 2，在sHbs上的ω值有明顯的提升；通過(guò)比較二比率模型和單比率模型，LRT檢驗(yàn)得到2Δ =20.09，P<0.001，df=1，說(shuō)明二比率模型更適合分析的數(shù)據(jù)組，作用于sHb上的選擇壓力明顯不同于其他的背景分支。同時(shí)將sHbs中的每個(gè)成員分別同擬南芥NsHb2成對(duì)比較dS和dN（圖1），發(fā)現(xiàn)sHbs每個(gè)成員的dS和dN同NsHb2相比都有不同程度的提升，說(shuō)明sHb亞家族經(jīng)歷了更多的選擇限制放松。

2.7.2 位點(diǎn)模型 蛋白質(zhì)的不同功能區(qū)域可能會(huì)經(jīng)受不同的選擇壓力，這種位點(diǎn)可變的選擇壓力可以通過(guò)包含不同位點(diǎn)分類的模型來(lái)檢測(cè)[35]。在檢測(cè)兩個(gè)小數(shù)據(jù)組的時(shí)候都分別采用了三對(duì)嵌套模型來(lái)成對(duì)檢驗(yàn)位點(diǎn)上的選擇壓力。對(duì)于只包含class 1的小數(shù)據(jù)組（表5）位點(diǎn)模型的比較中，離散模型M3要顯著優(yōu)于單比率模型M0，說(shuō)明位點(diǎn)間承受的選擇壓力具有異質(zhì)性；選擇模型M2a并不優(yōu)于近中性模型M1a（2Δl=0，df=2），而M8模型也不優(yōu)于M7模型（2Δl<0，df=2），因此從兩對(duì)模型的比較中無(wú)法發(fā)現(xiàn)正選擇（即使發(fā)現(xiàn)正選擇位點(diǎn)也不可信）。模型M3中3類密碼子位點(diǎn)p0（0.557 1）、p1（0.353 2） 和 p2（0.089 7）分別對(duì)應(yīng)的ω值為0.047 2，0.249 0， 0.858 2，說(shuō)明有55.7%左右的位點(diǎn)處于較強(qiáng)的凈化選擇約束；通過(guò)這3對(duì)位點(diǎn)模型比較結(jié)果，說(shuō)明class 1類基因整體上經(jīng)受凈化選擇或者中性選擇作用。對(duì)于僅含class 2的小數(shù)據(jù)組 （表5），得出同class 1數(shù)據(jù)組相類似的結(jié)論：在3對(duì)成對(duì)比較的模型中，M3顯著優(yōu)于M0，模型M3中3類密碼子位點(diǎn) p0（0.345 5）、p1（0.550 6） 和 p2（0.103 9） 分別對(duì)應(yīng)的ω值為0.052 0，0.330 7，0.916 4，說(shuō)明有33.6%左右的位點(diǎn)處于較強(qiáng)的凈化選擇約束。endprint

2.7.3 分支-位點(diǎn)模型 對(duì)于基因重復(fù)的研究表明，在基因重復(fù)完成后的不同時(shí)期，重復(fù)基因所承受的選擇的方式有所不同。自然選擇幫助一部分突變固定下來(lái)使生物產(chǎn)生多樣化的功能分歧，一旦新的或增強(qiáng)的功能建立之后，正選擇效應(yīng)終止，此時(shí)凈化選擇開始維持這些新功能[37]。本研究通過(guò)分支-位點(diǎn)模型來(lái)檢測(cè)在基因復(fù)制形成的class 2的祖先分支上可能發(fā)生的正選擇作用（表6）。Model A的參數(shù)評(píng)估結(jié)果表明，大約70%的氨基酸位點(diǎn)在所有的支系中是高度保守的（ω0=0.174 5），16%的位點(diǎn)是近中性的（ω1=1），14%的前景支上的位點(diǎn)處于正選擇壓力之下；進(jìn)一步鑒定出這些正選擇位點(diǎn)，其中在95%水平上有兩個(gè)（108P，126G），在99%水平上檢測(cè)到兩個(gè)（33L，132E）。為了檢測(cè)class 2中sHb的出現(xiàn)是否經(jīng)歷過(guò)適應(yīng)性進(jìn)化，采用僅含class 2的小數(shù)據(jù)組，將sHbs的祖先分支作為前景支進(jìn)行分支-位點(diǎn)的檢測(cè)。在該檢驗(yàn)中只檢測(cè)到1個(gè)位點(diǎn) （151K） 在統(tǒng)計(jì)上具有顯著性 （>0.95）。在以上這兩組檢測(cè)中，LRT成對(duì)檢驗(yàn)時(shí)，備則假說(shuō)模型都優(yōu)于零假說(shuō)模型，因此每個(gè)備則假說(shuō)模型檢測(cè)出的正選擇位點(diǎn)都是可信的。

3 結(jié)論

本研究通過(guò)分子進(jìn)化研究，闡明了被子植物血紅蛋白基因的系統(tǒng)發(fā)育關(guān)系及基因復(fù)制產(chǎn)生的class 1 和class 2兩類基因的分布格局、進(jìn)化機(jī)制以及在進(jìn)化中發(fā)生的功能性分歧和適應(yīng)性進(jìn)化，這些認(rèn)識(shí)為該基因家族進(jìn)一步的進(jìn)化研究奠定了基礎(chǔ)，同時(shí)也可以指導(dǎo)后續(xù)的試驗(yàn)工作，來(lái)驗(yàn)證那些具有重要功能的氨基酸置換位點(diǎn)和正選擇位點(diǎn)，揭示功能性分歧和適應(yīng)性進(jìn)化的分子基礎(chǔ)。

4 討論

4.1 血紅蛋白基因家族的起源和進(jìn)化

血紅蛋白基因在植物中廣泛存在，通過(guò)基因復(fù)制在被子植物中產(chǎn)生了class 1和class 2兩個(gè)亞家族，這兩個(gè)亞家族的基因具有不同的功能和表達(dá)模式[5，19，38]。本研究深入分析了被子植物中這兩類基因的系統(tǒng)發(fā)育及進(jìn)化特點(diǎn)，結(jié)果表明class 1和 class 2類血紅蛋白基因的姊妹關(guān)系。并充分利用數(shù)據(jù)庫(kù)現(xiàn)有的已完成測(cè)序或者正在組裝加工的單子葉植物 （水稻、玉米、小麥、大麥和高粱）序列信息，以擬南芥中非共生血紅蛋白AHB2作為搜索序列，在上述單子葉植物數(shù)據(jù)庫(kù)中做本地BLAST比對(duì)，沒(méi)有發(fā)現(xiàn)相似性比較高的序列，這些結(jié)果初步表明class 2類基因在單子葉植物中是缺失的。

本研究也估算了class 1和class 2的分歧時(shí)間為439.96 Ma，遠(yuǎn)早于單、雙子葉植物的分歧時(shí)間147 Ma[39]。因此，說(shuō)明class 1和class 2類血紅蛋白基因的出現(xiàn)是通過(guò)祖先基因發(fā)生的一次復(fù)制事件，該基因重復(fù)發(fā)生在單、雙子葉植物分化之前?？墒?，class 2類基因在單子葉植物中發(fā)生了丟失，致使形成了現(xiàn)在的分布格局，即class 1類基因廣泛分布于單、雙子葉植物中，而class 2類基因則僅見于雙子葉植物。

4.2 插曲式進(jìn)化和適應(yīng)性進(jìn)化

很多研究都表明很大比例的基因都是通過(guò)基因重復(fù)產(chǎn)生的?；蛑貜?fù)最主要的貢獻(xiàn)是為新的變異、漂變和選擇提供了豐富的基因材料[40]。正選擇可能發(fā)生在基因重復(fù)發(fā)生的早期而且可以使某些基因拷貝產(chǎn)生新的功能[41]。本研究中，在class 2的祖先支上檢測(cè)到正選擇的發(fā)生，這表明在基因重復(fù)發(fā)生之后，可能在其中一個(gè)拷貝中（class 2）的一些中性或者近中性的置換位點(diǎn)創(chuàng)造了新的但是較弱活性的功能，后來(lái)正選擇作用加速了這些有利突變的固定，并且提升了這些新功能的活性[42]。

本研究也揭示了在共生血紅蛋白基因出現(xiàn)時(shí)發(fā)生了插曲式的進(jìn)化速率加速和適應(yīng)性進(jìn)化。在圖1 中發(fā)現(xiàn)共生血紅蛋白的祖先支相對(duì)于其他非共生蛋白NsHb2的枝長(zhǎng)要長(zhǎng)，這表明在物種特異性的基因重復(fù)之后，共生血紅蛋白以較快的速率進(jìn)化，以插曲式的進(jìn)化模式產(chǎn)生了物種特異的適應(yīng)性；通過(guò)將 sHbs中的每個(gè)成員分別同擬南芥NsHb2成對(duì)比較dS和dN，發(fā)現(xiàn)在共生血紅蛋白中的非同義的置換速率和同義的置換速率都明顯加快（圖1），結(jié)果也進(jìn)一步證實(shí)了上述插曲式進(jìn)化的結(jié)論。在sHbs中置換速率的加快可能是祖先豆科物種中發(fā)生適應(yīng)性進(jìn)化的前提。

4.3 class 1和class 2的功能性分歧分析

現(xiàn)在普遍認(rèn)為基因重復(fù)可以導(dǎo)致基因功能的創(chuàng)新，這些包括假基因的形成[43，44]、亞功能化[45]、新功能化[46]和亞新功能化[11]，同時(shí)基因重復(fù)也可以為旁系同源基因之間產(chǎn)生多樣的功能性分歧創(chuàng)造機(jī)會(huì)。通過(guò)分析功能性分歧結(jié)果表明，class 1和class 2亞家族由于在某些氨基酸位點(diǎn)上具有顯著差異的進(jìn)化速率，導(dǎo)致了這兩個(gè)亞家族之間明顯的功能性分歧。通過(guò)分析這些功能性分歧位點(diǎn)在兩類基因中的組成特點(diǎn)發(fā)現(xiàn)，它們?cè)赾lass 1中非常保守而在class 2中高度可變。說(shuō)明通過(guò)基因復(fù)制產(chǎn)生的兩個(gè)重復(fù)基因中，class 1維持了原有功能，class 2經(jīng)歷了選擇限制的放松獲得了新的功能，即發(fā)生了新功能化。通過(guò)設(shè)置0.87和0.95兩個(gè)閾值，一些關(guān)鍵的功能性分歧位點(diǎn)被鑒定出來(lái)。在0.87水平下，預(yù)測(cè)了14個(gè)I類候選功能性分歧位點(diǎn)；當(dāng)水平提高到0.95時(shí)，只剩下7個(gè)候選 I 類功能性分歧位點(diǎn)。位點(diǎn)特異的圖譜分析（圖3）表明在14個(gè) I 類功能性分歧位點(diǎn)中只有1個(gè)位點(diǎn)經(jīng)受正選擇作用，說(shuō)明在這兩組旁系同源基因的某些位點(diǎn)中有著很明顯的置換速率轉(zhuǎn)換，這些位點(diǎn)的功能性分歧的產(chǎn)生或者通過(guò)放松的選擇限制或者通過(guò)正選擇作用。通過(guò)上述結(jié)果也說(shuō)明放松的功能限制是class 1和class 2中的功能性分歧位點(diǎn)形成的主要驅(qū)動(dòng)力。將這14個(gè)位點(diǎn)定位到血紅蛋白的三維結(jié)構(gòu)上發(fā)現(xiàn)，大部分位點(diǎn)分布在α-螺旋區(qū)，特別是在螺旋B、E和H （分別有5個(gè)、3個(gè)和3個(gè)位點(diǎn)），其他的均分布在CD-loop、螺旋G、螺旋E和F的連接區(qū) （每個(gè)區(qū)域均有1個(gè)位點(diǎn)分布）。這些功能性分歧位點(diǎn)定位的結(jié)構(gòu)組件中都包含著與血紅素結(jié)合的關(guān)鍵位點(diǎn)，它們定位在這里可能會(huì)引起血紅素口袋（Hemo pockete） 中的一些關(guān)鍵位點(diǎn)的位置改變，從而導(dǎo)致其在跟配基發(fā)生作用時(shí)的結(jié)合能力和結(jié)合模式發(fā)生變化，最終影響class 1和class 2兩類基因?qū)2的不同的結(jié)合特性。endprint

參考文獻(xiàn)：

[1] HUNT P W， WATTS R A， TREVASKIS B， et al. Expression and evolution of functionally distinct haemoglobin genes in plants [J]. Plant Mol Biol， 2001， 47（5）：677-692.

[2] ANDERSON C R， JENSEN E O， LLEWELLYN D J， et al. A new hemoglobin gene from soybean： A role for hemoglobin in all plants [J]. Proc Natl Acad Sci U S A， 1996， 93（12）：5682-5687.

[3] ARREDONDO-PETER R， HARGROVE M S， MORAN J F， et al. Plant hemoglobins [J]. Plant Physiol， 1998， 118（4）：1121-1125.

[4] STROZYCKI P M， KARLOWSKI W M， DESSAUX Y， et al. Lupine leghemoglobin I： Expression in transgenic Lotus and tobacco tissues [J]. Mol Gen Genet， 2000， 263（2）：173-182.

[5] TREVASKIS B， WATTS R A， ANDERSSON C R， et al. Two hemoglobin genes in Arabidopsis thaliana： The evolutionary origins of leghemoglobins [J]. Proc Natl Acad Sci U S A， 1997， 94（22）：12230-12234.

[6] VINOGRADOV S， HOOGEWIJS D， BAILLY X， et al. A phylogenomic profile of globins[J]. BMC Evol Biol， 2006， 6（1）：31-48.

[7] GULDNER E， DESMARAIS E， GALTIER N， et al. Molecular evolution of plant haemoglobin： Two haemoglobin genes in Nymphaeaceae Euryale ferox[J]. J Evol Biol， 2004， 17（1）：48-54.

[8] OHNO S. Evolution by Gene Duplication[M]. Berlin： Springer-Verlag， 1970.

[9] HUGHES A L. The evolution of functionally novel proteins after gene duplication[J]. Proc Biol Sci， 1994， 256（1346）：119-124.

[10] FORCE A， LYNCH M， PICKETT F B， et al. Preservation of duplicate genes by complementary， degenerative mutations [J]. Genetics， 1999， 151（4）：1531-1545.

[11] HE X，ZHANG J. Rapid subfunctionalization accompanied by prolonged and substantial neofunctionalization in duplicate gene evolution [J]. Genetics， 2005， 169（2）：1157-1164.

[12] DYKHUIZEN D， HARTL D L. Selective neutrality of 6PGD allozymes in E. coli and the effects of genetic background [J]. Genetics， 1980， 96（4）：801-817.

[13] YANG Z. Maximum likelihood estimation on large phylogenies and analysis of adaptive evolution in human influenza virus A [J]. J Mol Evol， 2000， 51（5）：423-432.

[14] ZHANG J， NIELSEN R， YANG Z. Evaluation of an improved branch-site likelihood method for detecting positive selection at the molecular level [J]. Mol Biol Evol， 2005， 22（12）：2472-2479.

[15] ARREDONDO-PETER R， HARGROVE M S， SARATH G， et al. Rice hemoglobins. Gene cloning， analysis， and O2-binding kinetics of a recombinant protein synthesized in Escherichia coli [J]. Plant Physiol， 1997， 115（3）：1259-1266.endprint

[16] KUMAR S， NEI M， DUDLEY J， et al. MEGA： A biologist-centric software for evolutionary analysis of DNA and protein sequences [J]. Brief Bioinform， 2008， 9（4）：299-306.

[17] TAMURA K， DUDLEY J， NEI M， et al. MEGA4： Molecular evolutionary genetics analysis （MEGA） software version 4.0 [J]. Mol Biol Evol， 2007， 24（8）：1596-1599.

[18] POSADA D，CRANDALL K A. MODELTEST： Testing the model of DNA substitution [J]. Bioinformatics， 1998， 14（9）：817-818.

[19] GUINDON S， LETHIEC F， DUROUX P， et al. PHYML Online--a web server for fast maximum likelihood-based phylogenetic inference [J]. Nucleic Acids Res， 2005， 33（Web Server issue）：W557-W559.

[20] WRIGHT F. The 'effective number of codons' used in a gene [J]. Gene， 1990， 87（1）：23-29.

[21] IKEMURA T. Correlation between the abundance of Escherichia coli transfer RNAs and the occurrence of the respective codons in its protein genes： A proposal for a synonymous codon choice that is optimal for the E. coli translational system [J]. J Mol Biol， 1981， 151（3）：389-409.

[22] ROBINSON-RECHAVI M，HUCHON D. RRTree： Relative-rate tests between groups of sequences on a phylogenetic tree [J]. Bioinformatics， 2000， 16（3）：296-297.

[23] DRUMMOND A，RAMBAUT A. BEAST： Bayesian evolutionary analysis by sampling trees [J]. BMC Evol Biol， 2007， 7（1）：214-322.

[24] NIELSEN R， YANG Z. Likelihood models for detecting positively selected amino acid sites and applications to the HIV-1 envelope gene [J]. Genetics， 1998， 148（3）：929-936.

[25] YANG Z. Likelihood ratio tests for detecting positive selection and application to primate lysozyme evolution [J]. Mol Biol Evol， 1998， 15（5）：568-573.

[26] YANG Z. Maximum likelihood estimation on large phylogenies and analysis of adaptive evolution in human influenza virus A [J]. J Mol Evol， 2000， 51（5）：423-432.

[27] YANG Z， WONG W S， NIELSEN R. Bayes empirical bayes inference of amino acid sites under positive selection [J]. Mol Biol Evol， 2005， 22（4）：1107-1118.

[28] GU X. Statistical methods for testing functional divergence after gene duplication [J]. Mol Biol Evol， 1999， 16（12）：1664-1674.

[29] GU X. Maximum-likelihood approach for gene family evolution under functional divergence [J]. Mol Biol Evol，2001， 18（4）：453-464.

[30] GU X， VANDER V K. DIVERGE： Phylogeny-based analysis for functional-structural divergence of a protein family [J]. Bioinformatics， 2002， 18（3）：500-501.

[31] GAUCHER E A， GU X， MIYAMOTO M M， et al. Predicting functional divergence in protein evolution by site-specific rate shifts [J]. Trends Biochem Sci， 2002， 27（6）：315-321.endprint

[32] KNUDSEN B， MIYAMOTO M M. A likelihood ratio test for evolutionary rate shifts and functional divergence among proteins [J]. Proc Natl Acad Sci U S A， 2001，98（25）：14512-14517.

[33] LANDGRAF R， XENARIOS I， EISENBERG D. Three-dimensional cluster analysis identifies interfaces and functional residue clusters in proteins[J]. J Mol Biol， 2001， 307（5）：1487-1502.

[34] CARLINI D B， STEPHAN W. In vivo introduction of unpreferred synonymous codons into the Drosophila Adh gene results in reduced levels of ADH protein[J]. Genetics， 2003， 163（1）：239-243.

[35] WANG M， ZHANG X， ZHAO H， et al. FoxO gene family evolution in vertebrates[J]. BMC Evol Biol， 2009， 9（1）：222.

[36] ZHONG B， YONEZAWA T， ZHONG Y， et al. Episodic evolution and adaptation of chloroplast genomes in ancestral grasses [J]. PLoS ONE， 2009， 4（4）：e5297.

[37] BIELAWSKI J P，YANG Z. Maximum likelihood methods for detecting adaptive evolution after gene duplication[J]. J Struct Funct Genomics， 2003， 3（1-4）：201-212.

[38] JACOBSEN-LYON K， JENSEN E O， JORGENSEN J E， et al. Symbiotic and nonsymbiotic hemoglobin genes of Casuarina glauca [J]. Plant Cell， 1995， 7（2）：213-223.

[39] HEDGES S B，KUMAR S. The Timetree of Life[M]. New York： Oxford University Press， 2009. 161-165.

[40] ZHANG J. Evolution by gene duplication： An update [J]. Trends Ecol Evol， 2003， 18（6）：292-298.

[41] MOORE R C，PURUGGANAN M D. The evolutionary dynamics of plant duplicate genes[J]. Curr Opin Plant Biol， 2005， 8（2）：122-128.

[42] ZHANG J， ROSENBERG H F， NEI M. Positive Darwinian selection after gene duplication in primate ribonuclease genes [J]. Proc Natl Acad Sci U S A， 1998， 95（7）：3708-3713.

[43] LYNCH M， CONERY J S. The evolutionary fate and consequences of duplicate genes[J]. Science， 2000， 290（5494）：1151-1155.

[44] WAGNER A. The fate of duplicated genes： Loss or new function[J] Bioessays， 1998， 20（10）：785-788.

[45] LYNCH M，F(xiàn)ORCE A. The probability of duplicate gene preservation by subfunctionalization [J]. Genetics， 2000， 154（1）：459-473.

[46] LYNCH M， O'HELY M， WALSH B， et al. The probability of preservation of a newly arisen gene duplicate[J]. Genetics， 2001，159（4）：1789-1804.endprint

[32] KNUDSEN B， MIYAMOTO M M. A likelihood ratio test for evolutionary rate shifts and functional divergence among proteins [J]. Proc Natl Acad Sci U S A， 2001，98（25）：14512-14517.

[33] LANDGRAF R， XENARIOS I， EISENBERG D. Three-dimensional cluster analysis identifies interfaces and functional residue clusters in proteins[J]. J Mol Biol， 2001， 307（5）：1487-1502.

[34] CARLINI D B， STEPHAN W. In vivo introduction of unpreferred synonymous codons into the Drosophila Adh gene results in reduced levels of ADH protein[J]. Genetics， 2003， 163（1）：239-243.

[35] WANG M， ZHANG X， ZHAO H， et al. FoxO gene family evolution in vertebrates[J]. BMC Evol Biol， 2009， 9（1）：222.

[36] ZHONG B， YONEZAWA T， ZHONG Y， et al. Episodic evolution and adaptation of chloroplast genomes in ancestral grasses [J]. PLoS ONE， 2009， 4（4）：e5297.

[37] BIELAWSKI J P，YANG Z. Maximum likelihood methods for detecting adaptive evolution after gene duplication[J]. J Struct Funct Genomics， 2003， 3（1-4）：201-212.

[38] JACOBSEN-LYON K， JENSEN E O， JORGENSEN J E， et al. Symbiotic and nonsymbiotic hemoglobin genes of Casuarina glauca [J]. Plant Cell， 1995， 7（2）：213-223.

[39] HEDGES S B，KUMAR S. The Timetree of Life[M]. New York： Oxford University Press， 2009. 161-165.

[40] ZHANG J. Evolution by gene duplication： An update [J]. Trends Ecol Evol， 2003， 18（6）：292-298.

[41] MOORE R C，PURUGGANAN M D. The evolutionary dynamics of plant duplicate genes[J]. Curr Opin Plant Biol， 2005， 8（2）：122-128.

[42] ZHANG J， ROSENBERG H F， NEI M. Positive Darwinian selection after gene duplication in primate ribonuclease genes [J]. Proc Natl Acad Sci U S A， 1998， 95（7）：3708-3713.

[43] LYNCH M， CONERY J S. The evolutionary fate and consequences of duplicate genes[J]. Science， 2000， 290（5494）：1151-1155.

[44] WAGNER A. The fate of duplicated genes： Loss or new function[J] Bioessays， 1998， 20（10）：785-788.

[45] LYNCH M，F(xiàn)ORCE A. The probability of duplicate gene preservation by subfunctionalization [J]. Genetics， 2000， 154（1）：459-473.

[46] LYNCH M， O'HELY M， WALSH B， et al. The probability of preservation of a newly arisen gene duplicate[J]. Genetics， 2001，159（4）：1789-1804.endprint

[32] KNUDSEN B， MIYAMOTO M M. A likelihood ratio test for evolutionary rate shifts and functional divergence among proteins [J]. Proc Natl Acad Sci U S A， 2001，98（25）：14512-14517.

[33] LANDGRAF R， XENARIOS I， EISENBERG D. Three-dimensional cluster analysis identifies interfaces and functional residue clusters in proteins[J]. J Mol Biol， 2001， 307（5）：1487-1502.

[34] CARLINI D B， STEPHAN W. In vivo introduction of unpreferred synonymous codons into the Drosophila Adh gene results in reduced levels of ADH protein[J]. Genetics， 2003， 163（1）：239-243.

[35] WANG M， ZHANG X， ZHAO H， et al. FoxO gene family evolution in vertebrates[J]. BMC Evol Biol， 2009， 9（1）：222.

[36] ZHONG B， YONEZAWA T， ZHONG Y， et al. Episodic evolution and adaptation of chloroplast genomes in ancestral grasses [J]. PLoS ONE， 2009， 4（4）：e5297.

[37] BIELAWSKI J P，YANG Z. Maximum likelihood methods for detecting adaptive evolution after gene duplication[J]. J Struct Funct Genomics， 2003， 3（1-4）：201-212.

[38] JACOBSEN-LYON K， JENSEN E O， JORGENSEN J E， et al. Symbiotic and nonsymbiotic hemoglobin genes of Casuarina glauca [J]. Plant Cell， 1995， 7（2）：213-223.

[39] HEDGES S B，KUMAR S. The Timetree of Life[M]. New York： Oxford University Press， 2009. 161-165.

[40] ZHANG J. Evolution by gene duplication： An update [J]. Trends Ecol Evol， 2003， 18（6）：292-298.

[41] MOORE R C，PURUGGANAN M D. The evolutionary dynamics of plant duplicate genes[J]. Curr Opin Plant Biol， 2005， 8（2）：122-128.

[42] ZHANG J， ROSENBERG H F， NEI M. Positive Darwinian selection after gene duplication in primate ribonuclease genes [J]. Proc Natl Acad Sci U S A， 1998， 95（7）：3708-3713.

[43] LYNCH M， CONERY J S. The evolutionary fate and consequences of duplicate genes[J]. Science， 2000， 290（5494）：1151-1155.

[44] WAGNER A. The fate of duplicated genes： Loss or new function[J] Bioessays， 1998， 20（10）：785-788.

[45] LYNCH M，F(xiàn)ORCE A. The probability of duplicate gene preservation by subfunctionalization [J]. Genetics， 2000， 154（1）：459-473.

[46] LYNCH M， O'HELY M， WALSH B， et al. The probability of preservation of a newly arisen gene duplicate[J]. Genetics， 2001，159（4）：1789-1804.endprint