田 郎 譯

(中國熱帶農(nóng)業(yè)科學院橡膠研究所,海南 儋州 571737)

巴西橡膠樹的基因組序列框架圖

田 郎 譯

(中國熱帶農(nóng)業(yè)科學院橡膠研究所,海南 儋州 571737)

巴西橡膠樹(Heveabrasiliensis)為大戟科多年生木本植物,也是迄今為止全世界商品天然橡膠(NR)的主要來源。NR是一種具有高彈性、高可塑性及高回彈力的膠乳聚合物，自1876年以來, NR在世界經(jīng)濟發(fā)展中一直都起著至關重要的作用。本文報道巴西橡膠樹的基因組序列框架圖。已組裝的基因組序列大約1.1 Gb,約占整個單倍體基因組估計長度(2.15 Gb)的51%。該基因組中大約有78%為重復DNA。通過基因預測共鑒定出68955個基因模型，其中12.7%為橡膠樹所特有。此外，大多數(shù)與橡膠生物合成、橡膠木材形成、病害抗性以及致敏性有關的關鍵基因也被成功鑒定。從該基因組獲得的序列信息將有助于人們進一步培育高產(chǎn)無性系，以充分滿足全世界對天然橡膠不斷增長的需求。

巴西橡膠樹；大戟科；天然橡膠；基因組

1 背景

橡膠是生產(chǎn)和制造5 000多種產(chǎn)品不可或缺的重要材料[1]。全世界大約有2 500種植物均能合成橡膠[2],但其中只有巴西橡膠樹(Heveabrasiliensis)是商品天然橡膠(NR)的主要來源。該大戟科物種起源于南美亞馬孫盆地，但直到19世紀才開始大規(guī)模用于商業(yè)開發(fā)，不過其人工馴化栽培并非始于其南美起源地。如今，人工栽培的橡膠樹主要分布于亞洲、非洲及拉丁美洲的熱帶地區(qū)。橡膠樹大約種植7 a后可開始收獲膠乳，其經(jīng)濟壽命約25～30 a。根據(jù)國際橡膠研究組織(www.rubberstudy.com)統(tǒng)計, 2011年全球天然橡膠的產(chǎn)量將近1100萬t,其中亞洲約占93%。在20世紀，全世界橡膠(包括天然和合成橡膠)的需求量一直呈現(xiàn)穩(wěn)步上升的態(tài)勢，預計今后還會繼續(xù)增長。

NR是一種具有彈性強、可塑性高、回彈力強、 抗沖擊力強以及散熱性良好等諸多性能的膠乳聚合物[2]。這些優(yōu)良性能使得NR在許多產(chǎn)品的制造上難以被合成橡膠所取代，如醫(yī)用手套及飛機和卡車的重型輪胎。NR由94%的順式-1,4-聚異戊二烯和6%的蛋白質(zhì)及脂肪酸組成[3]。順式-1,4-聚異戊二烯這種生物高聚物由含5個碳的異戊烯焦磷酸(IPP)單體構(gòu)成，并在橡膠粒子的表面通過縮合反應而生成。橡膠鏈的延伸由順式-異戊烯基轉(zhuǎn)移酶(CPTs)，也稱橡膠多聚酶催化[4]，所產(chǎn)生的聚合物的分子量是決定橡膠質(zhì)量的重要因素。僅有為數(shù)不多的植物能夠大量合成高質(zhì)量NR(分子量>100萬Da)，其中包括巴西橡膠樹以及潛在的橡膠樹替代作物銀膠菊(Partheniumargentatum) 和俗稱橡膠草的俄羅斯蒲公英(Taraxacumkoksaghyz)[5]。

除了NR之外，當橡膠樹進入膠乳生產(chǎn)末期而失去生產(chǎn)價值時，它們還可作為木材加以利用。如今，橡膠木已成為東南亞的一種主要出口木材[1]。橡膠木材的天然色澤及優(yōu)良物理性質(zhì)使其非常適于制作實木地板及家具。有鑒于此，人們業(yè)已培育出了多個膠木兼優(yōu)的橡膠無性系。

橡膠樹的病害及橡膠制品的致敏性也是橡膠業(yè)關注和擔憂的重要問題。真菌病害，例如由病原菌Microcyclusulei引起的南美葉疫病(SALB)，以及分別由病原菌Colletotrichum、Oidium和Corynespora引起的落葉病是橡膠生產(chǎn)的主要威脅[1]。在20世紀30年代中期，SALB摧毀了整個巴西的植膠業(yè)。亞洲的橡膠園雖然尚未受到該病害的侵襲，但此病一旦在該地區(qū)爆發(fā)將會產(chǎn)生毀滅性的后果。因直接接觸某些天然橡膠制品(如橡膠手套)而引發(fā)的人體過敏反應一直是醫(yī)學上關注的重要問題之一。引起這類過敏癥的過敏源主要為橡膠樹的膠乳中存在的某些蛋白質(zhì)。近年來，人們已開始嘗試將銀膠菊作為低致敏性膠乳的生產(chǎn)來源[2]。

常規(guī)育種方法的局限性及基因組信息的匱乏和不足已嚴重地阻礙和制約巴西橡膠樹的遺傳改良。分子標記輔助選擇通過對目標基因型的直接選擇可有效地提高育種效率。標記間的遺傳連鎖分析以及理想表型的遺傳定位將進一步提高選擇的準確度。近年來，高通量測序技術的迅猛發(fā)展及應用[6-12]在一定程度上拓寬并豐富了巴西橡膠樹的基因資源。然而迄今為止，該物種的全基因組信息依然十分欠缺。之前人們的研究大都集中于轉(zhuǎn)錄組分析,但事實上,基因組的非編碼區(qū)域?qū)α私饪刂苹虮磉_的調(diào)控元件以及建立一整套更加全面的分子標記必不可少。本文首次報道巴西橡膠樹的基因組序列框架圖,以期為該重要經(jīng)濟作物的進一步改良提供一個良好的平臺。

2 結(jié)果與討論

2.1 基因組測序及注釋

本研究對馬來西亞橡膠研究院培育的巴西橡膠樹高產(chǎn)無性系RRIM600進行了全基因組測序，該無性系的親本組合為Tjir1×PB86。巴西橡膠樹的基因組序列分布于18對染色體上[13]，通過孚耳根微顯影[14]其單倍體基因組估計約2.15 Gb。我們采用全基因組鳥槍法(whole-genome shotgun，WGS)并利用Roche/454, Illumina, 以及SOLiD測序平臺獲得大約43×覆蓋率的序列數(shù)據(jù)(表1)。由于大多數(shù)測序數(shù)據(jù)來自Roche/454平臺, reads(讀取片段)相對較長，尤其是鳥槍片段(shotgun)為單端測序,故我們選用Newbler軟件[15]進行最后的組裝[16]。對于初步組裝的序列我們還對其進行了重復基序的鑒定,并去除與重復區(qū)域匹配的原始測序reads，同時整個拼裝過程均采用嚴格的組裝參數(shù)。在濾除重復匹配reads之后，我們僅根據(jù)27.86Gb即大約13×覆蓋率的測序數(shù)據(jù)(表1)最終組裝獲得總長為1119Mb的序列支架(scaffolds),其N50 scaffold長度為2972bp(表2)。之后,根據(jù)154個微衛(wèi)星標記[17],我們將143個scaffolds以及1325個關聯(lián)基因錨定到了巴西橡膠樹的18個連鎖群上。在已定位的scaffolds中，還鑒定出另外74個已報道過的標記。這些標記大部分被定位在基因間隔區(qū)。

本研究基于全基因組鳥槍法(WGS)策略,并完全采用新一代測序技術成功構(gòu)建巴西橡膠樹的全基因組框架圖,而其它植物中僅有少數(shù)幾種采用類似方法獲得全基因組序列圖譜。與我們的研究方法一樣,Shulaev V等人在草莓中也是結(jié)合采用Roche/454，Illumina以及SOLiD平臺的reads(39×覆蓋率)構(gòu)建其全基因組圖譜，不過草莓基因組(240Mb)幾乎比橡膠樹小9倍，而且重復DNA的比例要小得多(22%)，在此情況下所得到的contigs/scaffolds相對要大得多[18]。然而,本研究中我們所拼裝得到的最大scaffold(531.5 kb)與大麻基因組利用Illumina和Roche/454測序數(shù)據(jù)組裝獲得的最大scaffold的長度(565.9 kb)幾乎相當[19]。橡膠基因組組裝的主要挑戰(zhàn)在于其大量重復序列的存在,這也是KFX Mayer等人在大麥基因組(5.1Gb，84%重復DNA)組裝過程中所遇到的一個難題，他們在全基因組鳥槍法測序中對Illumina平臺產(chǎn)生的短reads進行拼接僅獲得長度相對較短的contigs(N50=1425 bp)[20]。不過，當將其與BAC物理圖譜及高分辨率遺傳圖譜相結(jié)合時，依然可實現(xiàn)染色體水平上的序列重建。因此,對橡膠樹而言,下一步要做的就是結(jié)合物理圖譜或其它能夠提供更多序列毗連信息的方法以進一步改進其基因組拼裝的質(zhì)量。

利用基因組重復序列分析軟件RepeatModeler及RepeatMasker,在我們組裝的序列中有72.01%被鑒定為重復DNA(不包括低復雜度區(qū)域及RNA基因),由此估計重復序列約占橡膠樹基因組的78%,這與玉米(85%)[21]和大麥(84%)[20]相當。分析還顯示,長末端重復序列(LTR)類反轉(zhuǎn)錄轉(zhuǎn)座子為最主要的一類轉(zhuǎn)座元件(占總重復序列的46.15%),其中又以Gypsy類(38.20%)和Copia類(7.38%)為最多?？偟闹貜驮袃H有不到2%為DNA轉(zhuǎn)座子。重復元件中的大部分(50.24%)與任何已知的家族均無關。

注：* 來自Illumina測序平臺的 reads 通過去除與重復區(qū)域匹配的reads進行過濾，來自SOLiD測序平臺的 reads通過只保留長度至少為50 pb的雙末端測序reads進行過濾。

表2 巴西橡膠樹基因組組裝及注釋統(tǒng)計資料

注：*指將所有scaffolds按長度從長到短排序并累加，當累加長度達到所有scaffolds總長度50%時最后1個scaffold的長度；**指將所有scaffolds按長度從長到短排序并累加，當累加長度達到所有scaffolds總長度50%時scaffolds的總數(shù)。

結(jié)合來自多個從頭基因預測程序以及轉(zhuǎn)錄組和蛋白質(zhì)比對的證據(jù)和信息,我們采用EVidenceModeler(EVM)軟件[22]從已屏蔽的組裝序列中預測出68 955個基因模型，基因、外顯子及內(nèi)含子的平均長度分別為1 332 bp、238 bp及332 bp(表2)。在已公布的137 540個橡膠樹表達序列標簽(ESTs)及組裝的轉(zhuǎn)錄本中，有95.4%存在于該基因組中。為了給基因模型預測及確認提供進一步的支持，我們還進行了葉片轉(zhuǎn)錄組測序,并利用Roche/454及Illumina測序平臺分別獲得1 085 Mb和4.89 Gb的轉(zhuǎn)錄組序列。之后這些序列被從頭組裝進73 060個contigs(重疊群)。分析顯示，這些contigs中99%以上均與組裝的基因組序列相匹配，同時，來自Roche/454測序平臺的轉(zhuǎn)錄組讀取片段中有81%也與之相匹配。這些結(jié)果表明,我們組裝的基因組框架序列中已包含了該植物的大部分基因。

本研究借助不同數(shù)據(jù)庫對來自預測基因的蛋白質(zhì)序列作了進一步的注釋,這些數(shù)據(jù)庫包括NCBI非冗余數(shù)據(jù)庫，SwissPro蛋白序列數(shù)據(jù)庫[23], InterPro 蛋白質(zhì)家族數(shù)據(jù)庫[24], 以及KEGG生物學通路數(shù)據(jù)庫[25]。真核直系同源蛋白簇(KOG)[26]分析顯示,蛋白質(zhì)的數(shù)量在“信號轉(zhuǎn)導機制”( 5 216個),“翻譯后修飾, 蛋白質(zhì)周轉(zhuǎn), 分子伴侶”(2 886個),以及“碳水化合物轉(zhuǎn)運與代謝”(1 665個)等功能范疇上明顯占更高比例。此外, 富含亮氨酸重復(LRR)是該基因組中最為豐富的 Pfam[27]結(jié)構(gòu)域。在預測的基因中，有6.7%含有信號肽，其中大多數(shù)被定位在質(zhì)體及細胞外結(jié)構(gòu)中。

除了蛋白質(zhì)編碼基因之外，我們還鑒定出了729個tRNA基因，其中包括12個抑制(Sup)tRNA，32個假基因，以及4個功能未定的tRNA基因。令人感興趣的是，我們還觀察到了tRNA基因成簇分布的現(xiàn)象，例如所有Sup tRNA基因僅出現(xiàn)在2個scaffolds之中(9個在scaffold 134 351中，3個在scaffold 134 362中)。另外，我們在組裝的序列中還鑒定出了5S(113個拷貝)、5.8S(18個拷貝)、18S(11個拷貝)以及28S (21個拷貝) rRNA基因。

2.2 系統(tǒng)發(fā)生分析及種系特異性基因

利用來自17個已測序植物基因組的114個單拷貝直系同源基因簇進行系統(tǒng)發(fā)育分析，結(jié)果顯示，巴西橡膠樹與木薯(Manihotesculenta)共有一個最近的祖先(圖1)，這與該植物葉綠體基因組的鑒定分析結(jié)果也相一致[28]。除了大戟科之外，與之關系最為密切的已測序基因組為毛果楊(Populustrichocarpa)。與Shulaev V等人基于154個核基因所作被子植物系統(tǒng)發(fā)育分析的結(jié)果一致[18]，我們的分析結(jié)果也顯示，金虎尾目(包括楊柳科的楊屬及大戟科)與錦葵類植物的其它成員擁有一個共同的祖先。

對13個有代表性的植物基因組(分為大戟科、雙子葉植物、單子葉植物及低等植物4組)進行比較分析，結(jié)果發(fā)現(xiàn)，該4組基因組共有一套包括7140個基因家族的核心基因，而有9516個基因家族則為大戟科組所特有(圖2a)。比較4個已測序的大戟科植物基因組(即麻風樹、蓖麻、木薯和巴西橡膠樹),其結(jié)果顯示，共包含8748個基因的2708個基因家族為橡膠樹所特有(圖2b)。這些橡膠樹特有基因可被歸入526個GO[29]功能范疇，266個InterPro 結(jié)構(gòu)域，以及267個Pfam結(jié)構(gòu)域。最為豐富的InterPro 及Pfam結(jié)構(gòu)域是富含亮氨酸重復序列(LRR)及蛋白激酶。KOG分析顯示，這些基因中的大多數(shù)與信號轉(zhuǎn)導，細胞骨架，以及翻譯后修飾有關。

采用PhyML軟件并利用分布在17個已測序物種的144個單拷貝直系同源基因簇構(gòu)建基于最大似然法的系統(tǒng)發(fā)生樹。最優(yōu)樹的尋找選用SPR(子樹修剪與嫁接)和NNI(最近鄰居互換)算法中最佳的一個。分析顯示巴西橡膠樹在分類學上屬于錦葵類植物,并與木薯(M.esculenta)的親緣關系更近。 采用自舉法(bootstraping) 檢驗系統(tǒng)樹的可信度(隨機重復次數(shù)=100,隨機種子數(shù)目=7),自舉檢驗值(bootstraping值)示于圖中每個分支的節(jié)點處。

a.利用OrthoMCL軟件鑒定出13個植物物種特有及共有的基因家族；b. 利用OrthoMCL軟件鑒定出4個已測序大戟科物種特有及共有的基因家族

2.3 橡膠生物合成

橡膠生物合成涉及碳固定、同化物裝載運輸、橡膠前體物質(zhì)合成與代謝以及乳管中聚異戊二烯的儲存等諸多過程。蔗糖既為橡膠生物合成提供碳骨架又為合成提供能量,而乳管則是其強大的吸收匯[30]。本研究中,我們重構(gòu)了巴西橡膠樹橡膠生物合成的整個代謝途徑。橡膠的生物合成過程包括12個主要步驟(圖3)并涉及417個基因。我們還確證了這些來自葉片和(或)膠乳cDNA庫的基因的表達，并且在大多數(shù)基因家族中至少檢測到了1個異構(gòu)體。此外，我們將巴西橡膠樹的基因組與得自銀膠菊產(chǎn)膠樹皮組織的ESTs[31]進行了比較，結(jié)果發(fā)現(xiàn)，有360個橡膠樹的基因與銀膠菊ESTs相匹配，其中有205個基因與最佳匹配物的序列一致性高于70%。

注:白色小框中數(shù)字為每個步驟所涉及酶與蛋白的數(shù)量,黑色小框中數(shù)字為橡膠樹中直系同源體的數(shù)量。圖3 天然橡膠生物合成的12個主要途徑示意圖

蔗糖轉(zhuǎn)運蛋白及其表達模式與割膠和橡膠生產(chǎn)有直接關系[32]。蔗糖及單糖借助蔗糖和單糖轉(zhuǎn)運蛋白(SUT和MT)被輸送進乳管細胞溶質(zhì)，該2種蛋白質(zhì)分別由7個和30個基因編碼。β-呋喃果糖苷酶和果聚糖β-果糖甙酶能將蔗糖轉(zhuǎn)變成單糖,并且在巴西橡膠樹中發(fā)現(xiàn)大量此類基因(31個),這充分表明該功能在橡膠生物合成中的重要性。過多的蔗糖則以果聚糖和淀粉形式被儲存，并能夠用作以后橡膠生物合成的碳源。果聚糖代謝包括9個酶促步驟(由54個基因編碼)，而淀粉代謝則涉及11個酶促反應(由38個基因編碼)。碳通過醣酵解(由52個基因編碼)、可供選擇的磷酸戊糖途徑(由14個基因編碼)以及乙酰輔酶A生物合成途徑(由120個基因編碼)定向生成中間底物，用于橡膠前體物的生物合成。

用于橡膠生物合成的類異戊二烯前體由非自主型的細胞質(zhì)甲羥戊酸途徑(MVA)以異戊烯基焦磷酸(IPP)的形式提供[33]。有人認為自主型的質(zhì)體甲羥戊酸途徑(MEP)也能為橡膠生物合成提供IPP[34]。最近，橡膠樹實生苗的13C-標記研究顯示，MEP途徑為類胡羅卜素的生物合成提供IPP，但并不為橡膠生物合成提供IPP[35]。不過，MVA和MEP途徑基因的表達分析顯示，在并不產(chǎn)生大量類胡羅卜素的成齡橡膠樹或無性系中，MEP途徑很可能是提供IPP的另一條替代途徑[8]。在巴西橡膠樹基因組中，我們鑒定出了18個和29個分別參與MVA及MEP途徑的酶基因?；罨南┍沽姿?法尼基焦磷酸, 香葉基香葉基焦磷酸，或者十一異戊二烯焦磷酸)為引發(fā)橡膠生物合成所必需[36]。這些活化物的產(chǎn)生涉及5個酶促步驟并與組裝序列中的21個基因有關。

參與類異戊二烯聚合的橡膠聚合酶屬順式-異戊烯基轉(zhuǎn)移酶(CPTs)家族[37]。我們從橡膠樹的基因組中鑒定出了8個CPTs(命名為CPT 1-8)，并采用最大似然法和MEGA5.05軟件[38]推斷其進化歷史，根據(jù)進化關系將其分成3個類群(圖4)。我們發(fā)現(xiàn)，第2和第3類群中的巴西橡膠樹CPTs與其它植物的CPTs(十一異戊二烯焦磷酸合成酶及脫氫多萜醇焦磷酸合成酶)具有同源性，在這些植物中，CPTs 負責短鏈C5-異戊二烯(C55-C120)的延伸。由CPT 1-3組成的第1類群為巴西橡膠樹所特有,該組成員的作用被證明是催化長鏈C5-異戊二烯的形成[4]。這些CPTs中,僅第3類群中的CPT 4含有內(nèi)含子,而且該CPT與其它CPT的同源性最小。小橡膠粒子(SRPP)及橡膠延伸因子(REF)是參與橡膠生物合成的另外2種關鍵蛋白質(zhì)[5]，在我們組裝的序列中分別鑒定出其10個和12個相關基因。

2.4 橡膠木材

成年橡膠樹當進入膠乳生產(chǎn)末期后，還可作為一種木材來源用于家具及其它產(chǎn)品的制造。木材質(zhì)量與幾個木質(zhì)纖維素基因密切相關[39]。本研究中，我們在巴西橡膠樹中鑒定出了127個相關基因。與擬南芥和楊樹相比,我們在橡膠樹中鑒定出了36個纖維素合酶(CesA)編碼基因,而前2者中僅分別鑒定出10個[40]和18個[41]這樣的基因。本研究在巴西橡膠樹中也鑒定出與半纖維素生物合成有關的基因，而且α-L-巖藻糖苷酶的數(shù)量要多于α-L-巖藻糖基轉(zhuǎn)移酶，這與楊樹相似[42]。木質(zhì)素是一種由木質(zhì)素單體組成的雜聚物，它決定著木材的質(zhì)地和硬度。分析顯示，橡膠樹基因組中參與木質(zhì)素單體生物合成的一整套基因與楊樹基因組存在最大相似性。與楊樹基因組相比，橡膠樹中咖啡酸O-甲基轉(zhuǎn)移酶(COMT)和肉桂醇脫氫酶(CAD)的數(shù)量明顯要少(橡膠樹和楊樹中的COMT分別為10個和41個，CAD分別為5個和24個)。這可能和楊樹木材的硬度[39]與橡膠木材不同有關。本研究還鑒定出了一系列參與木質(zhì)素單體運輸、貯藏、調(diào)動以及其最后聚合形成木質(zhì)素的基因。

采取最大似然法并利用MEGA5.05軟件推斷其進化歷史。去除所有包含比對空格及缺少數(shù)據(jù)的位點。CPT，順式-異戊烯基轉(zhuǎn)移酶；UPPS，十一異戊烯焦磷酸合酶；DHDDS，脫氫多萜醇焦磷酸合酶；自舉檢驗值(重復100次)示于圖中分支節(jié)點處。

2.5 病害抗性

橡膠樹高度感染真菌病害，因此抗病基因的鑒定歷來是橡膠樹的研究重點之一。超敏反應(HR)是一種能引起植物組織壞死或細胞死亡從而達到限制病原菌生長的早期防御機制。據(jù)報道，水楊酸和茉莉酸作為植物信號分子在引發(fā)系統(tǒng)獲得性抗性(SAR)及誘導某些病程相關蛋白(PR)中起著關鍵作用[43]。核苷酸結(jié)合位點(NBS)編碼抗性基因家族在植物中是最大的一類抗性基因[44]。在巴西橡膠樹中，我們鑒定出了618個該家族成員，與水稻(Oryza sativa)中鑒定出的數(shù)量相當。這些基因被分成6個亞類，即toll-白細胞介素類似受體(TIR)-NBS、卷曲螺旋(CC)-NBS、NBS、TIR-NBS-LRR、CC-NBS-LRR以及NBS-LRR，其中多數(shù)不含LRR結(jié)構(gòu)域，而與此相比，其它植物中含有LRR的類型通常更為豐富。我們在組裝的基因組中還鑒定出147個病程相關蛋白(PR)和96個早期防御(SAR及HR)相關基因。所有這些病害抗性基因分布在665個序列支架中，并且NBS編碼基因通常被發(fā)現(xiàn)于這群序列中(如有9個NBS-LRR基因存在于骨架序列409956中)。此外，我們還重構(gòu)了巴西橡膠樹的植物系統(tǒng)獲得性抗病性(SAR)及超敏反應(HR)信號途徑。所有這些信息均有可能用于對植物生物脅迫的管控。

2.6 膠乳過敏源

天然橡膠膠乳(NRL)過敏癥是全球關注的一個重要醫(yī)學問題。國際公認的NRL過敏源有14種，分別被稱為Hevb 1至Hevb 14(www.allergen.org)。在橡膠樹中，這些過敏源被100個基因編碼。膠乳中的大多數(shù)過敏源是含量極為豐富且與脅迫和防御相關的蛋白[45]。導致對天然橡膠膠乳(NRL)過敏的主要過敏源中，Hevb 6 (橡膠素)由16個基因編碼，而Hevb 5則是一個單拷貝基因。與橡膠粒子有關的Hevb 1 (REF,橡膠延伸因子)和Hevb 3 (SRPP,小橡膠粒子)分別由12個和10個基因編碼。由5個基因編碼的Hevb 4 ( 卵磷脂酶同系物)和由9個基因編碼的Hevb 13 (酯酶)也被稱為含糖過敏源。編碼交叉反應蛋白Hevb 8 (抑制蛋白)、Hevb 9( 烯醇酶)以及Hevb 10 (錳超氧化物歧化酶)的基因分別有6個、4個及2個。防御相關過敏源Hevb 2(β-1,3-葡聚糖酶)及Hevb 11 (幾丁質(zhì)酶)各自有11個基因。Hevb 11的結(jié)構(gòu)域分析顯示，該結(jié)構(gòu)域存在1個長18～23個氨基酸的信號肽，此信號肽能引發(fā)植物的系統(tǒng)傷害反應[46]。除了上述的膠乳過敏源之外，我們還在橡膠樹中鑒定出了4類非膠乳過敏源，即花粉過敏源、α-膨脹素、β-膨脹素以及異黃酮還原酶。

2.7 轉(zhuǎn)錄因子

巴西橡膠樹的基因組大約含有6000個轉(zhuǎn)錄因子,它們分布在50個主要家族之中。轉(zhuǎn)錄因子約占巴西橡膠樹基因模型的8.5%,而其中bHLH、MYB、C3H、G2-like以及WRKY等家族的轉(zhuǎn)錄因子所占比例又大大超過其它家族。大多數(shù)植物中數(shù)量最多的bHLH轉(zhuǎn)錄因子家族在橡膠樹中多達752個成員。MYB是一個多樣化的轉(zhuǎn)錄因子家族，該家族的轉(zhuǎn)錄因子與bHLH家族轉(zhuǎn)錄因子互作從而調(diào)節(jié)次生代謝[47]及生物和非生物脅迫，在橡膠樹中該家族成員的數(shù)量也多達570個。C3H家族轉(zhuǎn)錄因子涉及花發(fā)育、胚胎發(fā)生、越冬以及葉片衰老[48]，其成員數(shù)為470個。之后為G2-like和WRKY家族的轉(zhuǎn)錄因子，它們在橡膠樹中的成員數(shù)分別為461和445個。據(jù)報道該2個轉(zhuǎn)錄因子家族分別參與了植物的光合調(diào)控[49]及免疫應答[50]。此外，橡膠樹中成員數(shù)較多的轉(zhuǎn)錄因子家族還有MYB-related(397個)、NAC(336個)以及ERF(246個)等。MADS-box基因是一類編碼與花器官發(fā)育有關的轉(zhuǎn)錄因子的同源異型基因。MADS-box基因被分為5類，即Mα、Mβ、Mγ、 Mδ (即 MIKC*) 和MIKCc[51]，它們在橡膠樹中共有112個成員。在這112個成員中有79個屬于II型MADS-box基因(Mδ和MIKCc組)，而該類型基因在擬南芥、楊樹及水稻中的數(shù)量為54～67個。與此相比, I型的 MADS-box基因 (包括Mα、Mβ及Mγ組)在橡膠樹中的數(shù)量只有33個，而其它3個物種中則存在29～94個。此外, 巴西橡膠樹中的MADS-box基因僅有12.5%(112個中的14個)以成簇形式存在,而與此相比，番木瓜的MADS-box基因中則有47% 在基因組中成簇存在[52]。

2.8 植物激素生物合成及信號轉(zhuǎn)導

分析顯示, 巴西橡膠樹中存在為數(shù)眾多的植物激素生物合成及信號轉(zhuǎn)導相關基因。被子植物基因組通常都含有較大比例的生長素信號通路相關基因,已鑒定出的基因家族達12個。不過，與其它植物相比，橡膠樹中有些生長素基因家族成員明顯減少，尤其是生長素上調(diào)小RNA(SAUR)及生長素(IAA)阻遏蛋白家族成員的數(shù)量顯著少于楊樹、擬南芥及水稻。赤霉素-20-氧化酶是赤霉素生物合成中的一種關鍵調(diào)節(jié)酶，在巴西橡膠樹中共鑒定出5個直系同源基因，而蓖麻及水稻中僅存在1個這樣的基因。乙烯響應元件結(jié)合因子(ERF)在橡膠樹中共有246個直系同源基因，該數(shù)量大大多于其它植物。ERF轉(zhuǎn)錄因子數(shù)量的增加可能與橡膠樹中特定的乙烯依賴途徑有關。12-氧-植物二烯酸還原酶(OPR7)是茉莉酸生物合成途徑中的一個關鍵酶，在橡膠樹中存在13個編碼基因。一氧化氮合成酶參與一氧化氮的生物合成，而一氧化氮在很多植物抗病防御機制中起著特別重要的作用。在擬南芥、蓖麻、楊樹以及水稻中，一氧化氮合成酶基因是一個高度保守的單拷貝基因，但在巴西橡膠樹種卻有4個拷貝。

2.9 光信號轉(zhuǎn)導及生物鐘相關基因

光信號轉(zhuǎn)導通路與生物鐘相互聯(lián)系并對植物的生理活動具有十分重要的影響。光是最重要的環(huán)境信號之一，來自環(huán)境的光信號作用于植物內(nèi)在的計時機制,即生物鐘使植物的生理節(jié)奏與環(huán)境的晝夜及季節(jié)性變化相一致[53]。與楊樹和擬南芥相比,巴西橡膠樹中涉及光感受及晝夜節(jié)律的基因數(shù)明顯增加,數(shù)量多達154個,而楊樹和擬南芥中的數(shù)量分別僅有77和66個。這些結(jié)果表明，環(huán)境信號強烈參與橡膠樹生理活動及機能的調(diào)控。

2.10 F-box蛋白質(zhì)

F-box蛋白質(zhì)是SCF(Skp1p-cullin-F-box)蛋白質(zhì)復合體的組成部分，在植物中SCF復合物廣泛參與泛素/26S蛋白酶體途徑所介導的蛋白質(zhì)選擇性降解[54]。F-box蛋白質(zhì)在其N-端存在一個保守F-box結(jié)構(gòu)域(40-50個氨基酸)。據(jù)報道,這類蛋白質(zhì)參與植物眾多生長發(fā)育過程的調(diào)節(jié),如葉片衰老、開花、分枝、光敏素與植物激素信號轉(zhuǎn)導以及自交不親和等[55]。巴西橡膠樹基因組中含有655 個F-box基因,而該家族基因在釀酒葡萄(Vitisvinifera)、番木瓜 (Caricapapaya)、毛果楊、擬南芥 (Arabidopsisthaliana)以及水稻中的數(shù)量分別為315、198、425、897以及971個[56]。這一結(jié)果頗令人感興趣，因為與人們的預期正好相反,一年草本生植物的F-box基因家族其基因數(shù)量較之木本的多年生植物反而出現(xiàn)了明顯的擴張[55]。

2.11 類胡蘿卜素

類胡蘿卜素在光捕獲、光保護、光形態(tài)建成、脂質(zhì)過氧化以及植物的很多發(fā)育過程中都起著至關重要的作用[57]。類胡蘿卜素在包括橡膠膠乳中的弗萊-威士林粒子(即FW粒子)在內(nèi)的幾乎所有類型質(zhì)體中都有所發(fā)現(xiàn)，而正是FW粒子的存在使得有些橡膠無性系的膠乳呈現(xiàn)出黃色。盡管膠乳中類胡蘿卜素的作用尚不十分確定和清楚，但它可能是乳管中異戊烯焦磷酸(IPP)的一個競爭匯。有人認為，在膠乳類胡蘿卜素含量低的橡膠無性系中，來自2C-甲基-D-赤蘚醇-4-磷酸途徑(MEP)的IPP被用于順式-聚異戊二烯的合成[8]。分析顯示，橡膠樹中涉及類胡蘿卜素生物合成途徑的基因數(shù)量(48個)較之擬南芥基因組中相關基因的數(shù)量(28個)明顯增多。八氫番茄紅素合成酶和八氫番茄紅素脫氫酶是催化類胡蘿卜素生物合成中第一個關鍵步驟的酶，它們在橡膠樹中的數(shù)量高度增加，基因數(shù)分別達5個和9個，而在擬南芥中均為單拷貝基因。這些觀察結(jié)果表明，橡膠樹具有更加有效的類胡蘿卜素生物合成機制，而類胡蘿卜素在該植物中的合成可能有著多方面的作用和功能。

3 結(jié)論

鑒于天然橡膠生產(chǎn)及其可持續(xù)性在世界工業(yè)及經(jīng)濟發(fā)展中的重要作用和地位,橡膠樹基因組序列框架圖的構(gòu)建無疑為大戟科植物研究增添了一份寶貴的資料，而這也將有力促進和加速橡膠樹分子育種的開展及基因資源的開發(fā)利用。我們的觀察顯示，巴西橡膠樹的基因組中存在較大比例的重復序列(大約78%)，這有可能緣于該植物非同源重組率及外顯子重排事件的增加[58]，并由此也降低了后代遺傳純度的一致性。大量重復序列的存在及染色體水平上的信息的缺乏是巴西橡膠樹基因組組裝的主要障礙。來自全基因組的信息以及與目標基因連鎖或相關聯(lián)的分子標記的鑒定將有助于通過分子育種實現(xiàn)一些重要性狀的遺傳改良，從而有力地促進天然橡膠生產(chǎn)的發(fā)展。與木本植物楊樹、桉樹及草本模式植物擬南芥相比，全基因組信息的開發(fā)利用將特別有助于橡膠樹在膠乳產(chǎn)量,木材開發(fā),抗病性,以及降低膠乳致敏性等關鍵領域取得突破性進展。

4 方法

4.1 基因組測序及組裝

從巴西橡膠無性系RRIM600的嫩葉中提取高質(zhì)量染色體DNA。按試劑盒使用說明書制備鳥槍(Shotgun)及雙末端(Paired-end,PE)法測序文庫。之后,分別采用Illumina、Roche/454以及SOLiD測序儀產(chǎn)生高質(zhì)量的讀取片段(reads)。Illumina 測序文庫包括200 bp及500 bp PE，Roche/454 測序文庫包括shotgun、8 kb PE以及20 kb PE, SOLiD 測序文庫為2 kb PE。

采用2個為新一代測序數(shù)據(jù)的從頭組裝而設計的軟件[59],即CLC Workbench (CLC bio, Denmark)和Newbler (version 2.3)進行序列初步組裝,該2個組裝軟件各自采用不同的輸入數(shù)據(jù)內(nèi)容及組裝策略。借助基于de Bruijn圖的短序列拼接算法，從經(jīng)過質(zhì)量修剪的Illumina 200bp reads生成CLC組裝的basic contigs,所有經(jīng)過修剪的Illumina reads、Roche/454 reads以及SOLiD reads則用于contigs(重疊群)的延伸和連接。Newbler 組裝的組裝參數(shù)設定為:大的或復雜的基因組，reads限定到某一個contig中，最小重疊長度50bp，最小重疊一致性95%。每個組裝軟件初步組裝的序列中，長度≥200bp的Contigs被留下用于進一步分析。

采用RepeatModeler(version 1.0.4) 軟件包[60],并按默認參數(shù)對2個組裝軟件初步拼裝的序列進行重復序列的從頭識別。我們分別從CLC及Newbler軟件的初步組裝序列中鑒定出2323個和1520個重復模塊。之后，通過BLASTX相似性搜索(采用E-值界限10-5)，從來自該2個初步組裝序列的重復文庫中篩選與NR(non-redundant,非冗余)、KEGG以及TrEMBL[61]蛋白質(zhì)數(shù)據(jù)庫中未分類重復序列相關的潛在基因家族,并將其合并為一個含有3771個重復序列的巴西橡膠樹特異性重復文庫。以重復序列的出現(xiàn)頻率為標準，對該巴西橡膠樹重復文庫進行篩選，凡是每一個初裝序列中出現(xiàn)100次以上的重復序列被保留下來，并將其組合成一個巴西橡膠樹高頻重復文庫，用于軟件RepeatMasker (version 3.2.9)[62]的輸入。之后將其用于鑒定和屏蔽Newbler初步組裝序列中的重復區(qū)域，但不包括序列中的低復雜度區(qū)域和RNAs。將已屏蔽重復的Newbler初步組裝序列作為模板用以篩選通過Illumina平臺產(chǎn)生的與重復序列相關的測序reads。在進行篩選之前，Illumina reads已經(jīng)過質(zhì)量控制且所有測序堿基位的質(zhì)量值≥25的reads被保留下來,其中200 bp測序文庫的read長度為100 bp，500 bp的雙末端測序文庫其正向和反向read長度分別為85 bp和75 bp。采用BOWTIE (version 0.12.7)軟件將每一個經(jīng)過質(zhì)量篩選的高質(zhì)量Illumina read的前50 bp與Newbler初步組裝的序列進行比對，允許的錯配堿基不超過3個。與重復區(qū)域匹配的成對或非成對reads均從read數(shù)據(jù)集中排除。就通過SOLiD平臺產(chǎn)生的測序reads而言,其質(zhì)量控制首先是將SOLiD reads與已經(jīng)SOAP軟件[63]校正的Illumina 200 bp reads用較早版本CLC軟件初步拼裝所得到的序列進行比對，任何種類的匹配錯誤(包括彩色空隙、單核苷酸差異，或者插入缺失)允許出現(xiàn)2個。兩條reads都匹配的SOLiD 成對reads其對應的參考序列被分離出來，并將其作為一個read對。將所有兩端讀長至少達到50 bp的read對保留于read數(shù)據(jù)集中。

利用所有Roche/454 reads、部分挑選出的Illumina reads(包括來自200 bp測序文庫的成對或非成對reads和來自500 bp測序文庫的成對reads)以及SOLiD 成對reads，并采用Newbler組裝軟件獲得最終的基因組組裝序列(表1)。從TIGR植物重復序列數(shù)據(jù)庫及TIGR玉米重復序列數(shù)據(jù)庫獲得其它植物的重復序列文庫，并且核糖體DNA序列已從這些數(shù)據(jù)庫中被去除。將巴西橡膠樹特異性重復文庫與TIGR植物重復序列數(shù)據(jù)庫合并為RepeatMasker的輸入重復文庫，以鑒定和屏蔽最終拼裝所得contigs中的重復區(qū)域，但不包括其所含低復雜度區(qū)域及RNAs基因。采用 BOWTIE程序(version 0.12.7)[64]將每一個基因組測序reads的前50 bp比對到最終的基因組拼裝序列上,同時利用TopHat (version 1.1.4)軟件[65]將通過Roche/454平臺產(chǎn)生的轉(zhuǎn)錄組測序reads與最終的基因組拼裝序列進行比對，以評估編碼區(qū)域的完整性和覆蓋度。

為了鑒定帶有細胞器來源的contigs，采用BLASTN搜索程序?qū)⑵囱b序列與巴西橡膠樹葉綠體基因組序列進行相似性比對，同時，還通過BLASTX程序?qū)⑵囱b序列與來自豆目植物細胞器基因組，巴西橡膠樹、麻風樹 (Jatropha curcas)和木薯的葉綠體基因組，以及蓖麻(Ricinuscommunis)、西瓜(Citrulluslanatus)和西葫蘆(Cucurbitapepo)線粒體基因組的蛋白質(zhì)進行比對分析。源于細菌污染的contigs則通過搜索GenBank數(shù)據(jù)庫進行鑒定并將其從最終的拼裝序列中去除。

4.2 轉(zhuǎn)錄組測序

從巴西橡膠樹的葉片中分離總RNA。按試劑盒使用說明書步驟，將分離得到的mRNA隨機打斷后，通過逆轉(zhuǎn)錄合成雙鏈cDNA片段并用其構(gòu)建測序文庫，之后采用Illumina和Roche/454平臺測序。利用CLC Workbench組裝軟件并通過拼接Illumina reads獲得最初的轉(zhuǎn)錄組序列。將來自Illumina轉(zhuǎn)錄組組裝的Contigs切割成最大1999 bp的短片段,并與Roche/454轉(zhuǎn)錄組測序reads相結(jié)合，用于Newbler組裝程序的輸入用以優(yōu)化EST數(shù)據(jù)。最后,通過對非冗余蛋白質(zhì)數(shù)據(jù)庫進行BLASTX搜索(E-值界限10-5)對轉(zhuǎn)錄組組裝序列的Contigs進行注釋,并以此檢測轉(zhuǎn)錄組的完整性及多樣性。

4.3 基因組注釋

利用EVidenceModeler(EVM；version r03062010)軟件，通過整合來自轉(zhuǎn)錄組比對、蛋白質(zhì)比對以及基因從頭預測的證據(jù)最終完成橡膠樹已屏蔽基因組的基因模型預測。利用PASA (version r09162010)[66]和Exonerate(version 2.2.0)[67]軟件將來自橡膠樹轉(zhuǎn)錄組組裝的Contigs與其基因組序列進行比對。采用GMAP (version 20100727)程序[69]將從PlantGDB數(shù)據(jù)庫[68]獲得的植物唯一轉(zhuǎn)錄本(PUTs)與基因組進行比對。利用AAT(version 1.52)[70]及BLAT(version 34)[71]軟件將從PlantGDB數(shù)據(jù)庫所獲基因組測序計劃的植物蛋白質(zhì)序列與組裝的橡膠樹基因組進行比對。將來自橡膠樹轉(zhuǎn)錄組組裝的Contigs作為訓練集，用以訓練從頭預測軟件Fgenesh[72]。之后，將基于轉(zhuǎn)錄組比對的PASA預測結(jié)果進一步作為訓練集并用以訓練從頭開始的基因預測軟件Augustus (version 2.5)[73]、 GlimmerHMM (version 3.0.1)[74]以及SNAP(version 2010-07-28)[75]。此外，在結(jié)構(gòu)注釋過程中,我們還一并使用了用擬南芥訓練的從頭預測軟件GeneMarkHMME (version 3)[76]及用甜瓜(Cucumis melo)訓練的從頭預測軟件Geneid (version 1.4.4)[77]進行橡膠樹的基因模型預測。最后,利用EVM軟件整合來自不同方法及軟件的預測結(jié)果以獲得基因的一致性序列。綜合考慮各種預測方法的可信度及有效性等因素，我們以橡膠樹轉(zhuǎn)錄組比對的PASA預測結(jié)果作為最高標準，將各種來源證據(jù)的權(quán)重人工設置為，橡膠樹轉(zhuǎn)錄組組裝：PASA 1，Exonerate 0.5；植物PUTs:GMAP 0.2；植物蛋白質(zhì):AAT 0.2,BLAT 0.2；基因從頭預測：Fgenesh 0.6,Augustus 0.5,SNAP 0.3,GlimmerHMM 0.3,Geneid 0.2,GeneMarkHMME 0.2。

為了確保質(zhì)量和精確的注釋,除了常用的功能注釋程序外,我們還設置了人工校對及若干驗證標準。采用BLASTP程序?qū)碜灶A測基因的蛋白質(zhì)序列與數(shù)據(jù)庫Swiss-Prot、TrEMBL、PlantGDB、UniRef100、NCBI非冗余數(shù)據(jù)庫、STRING (version 8.3)[78]以及KEGG GENES進行相似性比對(E-值界限10-5)，從而得到相應的功能注釋?；陂_放閱讀框(ORFs)的功能篩查其長度覆蓋率及相似性均至少達到70%。對未經(jīng)完整ORFs得到注釋的序列，則進一步通過掃描InterPro、PANTHER[79]、PRINTS[80]、PROSITE[81]patterns、Pfam以及SMART[82]等數(shù)據(jù)庫進行結(jié)構(gòu)域檢測，隨后與已經(jīng)得到明確注釋的序列模板(如擬南芥和蓖麻)進行比對得到相關功能信息。之后,基于相互最佳匹配法并借助Pathway Studio (Ariadne Genomics Inc.)軟件從嚴格審核的植物特定數(shù)據(jù)庫中尋找直系同源物，并以此對功能注釋作進一步的分類。利用Pathway Studio功能類別及KEGG直系同源評估獲得酶的EC分類號。根據(jù)STRING及GENES數(shù)據(jù)庫的BLASTP匹配片段獲得其KO分類號。從InterPro數(shù)據(jù)庫的搜索結(jié)果獲得蛋白質(zhì)對應的GO分類號。通過蛋白質(zhì)與PlantRefSeq、KEGG、Swiss-Prot以及InterPro數(shù)據(jù)庫的比較進行人工審核校對。超過10 000條蛋白質(zhì)的相應功能得到人工注釋。鑒定百分率、序列比對可靠性得分(bit-score)以及長度覆蓋率的比較連同其它分析通過計算機的綜合處理分析最終使一個特定ORF的推定功能得到可靠認定。

利用TRNAscan-SE v.1.23軟件包且EufindRNA檢測軟件按非嚴格設置(Int Cutoff=-32.1)鑒定所組裝的基因組序列中的tRNA 基因[83]。采用BLASTN 2.2.24軟件，通過與擬南芥及水稻所含5S、5.8S、18S及28S rRNA的比對鑒定rRNA基因(長度覆蓋率至少為80%，序列一致性至少達到50%[84])。利用SignalP 4.0軟件預測和鑒定基因組中的信號肽[85]。

4.4 巴西橡膠樹基因家族的鑒定和注釋

利用CLC軟件并采用合適的模板序列鑒定與橡膠及木質(zhì)纖維素生物合成、系統(tǒng)獲得性抗性、超敏反應、病程相關蛋白、過敏源、轉(zhuǎn)錄因子、植物激素代謝及信號轉(zhuǎn)導、生物鐘及光信號轉(zhuǎn)導、F-box以及類胡蘿卜素生物合成等相關的基因。通過在PlantGDB、UniProtKB/TrEMBL以及Plant_refseq蛋白質(zhì)數(shù)據(jù)庫進行BLASTX搜索(E-值<10-5)對鑒定的基因作進一步注釋。

4.5 NBS-LRR基因家族的鑒定和注釋

采用BLASTP程序(E-值<10-5)從PlantGDB、NCBI以及TAIR數(shù)據(jù)庫中搜索與擬南芥NBS-LRR(核苷酸結(jié)合位點—富含亮氨酸重復序列)蛋白質(zhì)相似且覆蓋率大于90%的巴西橡膠樹蛋白質(zhì)，并通過與NCBI保守結(jié)構(gòu)域數(shù)據(jù)庫進行結(jié)構(gòu)域匹配加以進一步確認。之后，利用數(shù)據(jù)庫Pfam、InterPro及HMMPanther對這些來自橡膠樹基因組的已注釋基因模型進行結(jié)構(gòu)域掃描和搜索，并分別獲得相應的模體及ID號，其中包括：TIR[PF01582、IPR000157]，NBS[PF0931、IPR002182]，TIR-NBS [PF01582、PF00931、IPR000157、IPR002182]，NBS-LRR[PF0931、PF00560、PF07723、PF07725及PTHR23155:SF236]，TIR-NBS-LRR[PF01582、PF0931、PF00560、PF07723、PF07725、IPR000157、IPR002182、IPR001611、IPR011713及PTHR23155:SF300]，CC-NBS-LRR[PTHR23155: SF231]，以及三種類型的LRR，如LRR_1[PF00560、IPR001611]、LRR_2[PF07723]以及LRR_3[PF07725、IPR011713]。利用COILS 程序預測和鑒定[86]卷曲螺旋(CC)結(jié)構(gòu)域。通過合并分析、人工驗證及片段匹配度的檢驗，我們從巴西橡膠樹中鑒定出了618個NBS-LRR基因。

4.6 橡膠生物合成相關基因與銀膠菊ESTs的比較

將來自巴西橡膠樹與橡膠生物合成相關的基因翻譯成蛋白質(zhì)序列，并用作查詢序列與銀膠菊ESTs[GenBank:GW775573-GW 787311]作TBLASTN比對分析。所得結(jié)果采用E-值界限10-5進行甄別和過濾。

4.7 通路重建

利用Resnet-Plant 3.0數(shù)據(jù)庫及代謝通路數(shù)據(jù)庫(Metabolic Pathway Databases，MPW)[87]并借助Pathway Studio software軟件(Ariadne Genomics Inc.)重建巴西橡膠樹的代謝通路。Resnet-Plant 3.0數(shù)據(jù)庫集成了從擬南芥AraCyc 4.0數(shù)據(jù)庫引入的303條代謝通路信息，每條代謝通路體現(xiàn)為一系列不同功能類群的蛋白質(zhì)(酶)及一組相應的化學反應。數(shù)據(jù)庫中的每一個功能類群可包含數(shù)量不限并具有相應酶活性編碼序列的蛋白質(zhì)成員。同一功能類群中通常包含若干行使酶活功能所必需的催化及調(diào)節(jié)亞基的旁系同源物。

利用Pathway Studio初步重建的代謝通路表征為一系列不同功能類群的人工審核蛋白質(zhì)。該重建過程實際上也相當于一個在代謝網(wǎng)絡中填補空缺的過程。帶有橡膠基因組標識的蛋白質(zhì)在Resnet-Plant 3.0數(shù)據(jù)庫中得到注釋并去除非橡膠蛋白質(zhì)之后，有311個功能類群不存在任何與其對應的蛋白質(zhì)成員。利用根據(jù)BLASTP結(jié)果采用雙向最佳匹配法鑒定出的直系同源蛋白質(zhì),并通過對已組裝橡膠基因組DNA序列的TBLASTN分析人工鑒定自動注釋未能發(fā)現(xiàn)的蛋白質(zhì)。填補橡膠代謝網(wǎng)絡中的空缺首先是從GenBank及UniProt數(shù)據(jù)庫下載能執(zhí)行所缺失酶蛋白功能的蛋白質(zhì)序列，然后將其作為查詢序列用于TBLASTN分析。這些用于查詢的蛋白質(zhì)序列或者來自擬南芥基因組，或者是來自其它植物或細菌基因組。通過與擬南芥AraCyc數(shù)據(jù)庫的比較，鑒定僅僅存在于水稻RiceCyc或楊樹PoplarCyc數(shù)據(jù)庫中的代謝通路。擬南芥AraCyc數(shù)據(jù)庫中所缺少的代謝通路采用手動方式添加到Pathway Studio的巴西橡膠數(shù)據(jù)庫中。MPW數(shù)據(jù)庫也被用于橡膠生物合成途徑的重建。

4.8 將序列支架錨定到連鎖群

根據(jù)已公布的連鎖圖譜[17]，將組裝的scaffolds錨定和定位到橡膠樹的18個連鎖群。所用154個微衛(wèi)星標記序列均得自公共序列數(shù)據(jù)庫。通過標記序列與所有scaffolds之間的 BLAST相似性分析對各個scaffolds進行鑒定，而所在scaffold上的基因模型也因此得以鑒定并被錨定在連鎖群的相應位置。當某個scaffold上存在一個以上標記時，基因則有可能被錨定在正確的位置上，也有可能在其它位置，其定位不確定。位于scaffolds上的其它額外標記則通過全部scaffolds對GenBank數(shù)據(jù)庫的BLAST相似性分析進行鑒定。

4.9 特有及共有基因簇分析

采用OrthoMCL軟件[88]鑒定和估算大戟科植物以及不同植物類群間旁系和直系同源基因簇的數(shù)量。采用標準設置(BLASTP, E-值<10-5)計算序列的多對多相似性。

4.10 系統(tǒng)發(fā)生分析

采用17個已測序的植物基因組,即萊茵衣藻(Chlamydomonasreinhardtii)、江南卷柏(Selaginellamoellendorfii)、玉米(Zeamays)、水稻、二穗短柄草(Brachypodiumdistachyon)、馬鈴薯(Solanumtuberosum)、葡萄、番木瓜、擬南芥、毛果楊、麻瘋樹、蓖麻、巴西橡膠樹、木薯、野草莓(Fragariavesca)、大豆(Glycinemax)以及黃瓜(Cucumissativus)構(gòu)建系統(tǒng)發(fā)生樹。采用BLAST算法及E-值界限10-5對蛋白質(zhì)序列進行多對多的相似性比較，并從該BLAST分析結(jié)果計算鑒定百分率。利用OrthoMCL軟件鑒定內(nèi)旁系同源基因、直系同源基因以及共直系同源基因。通過尋找某一物種內(nèi)的所有成對蛋白質(zhì)確定潛在的內(nèi)旁系同源基因?qū)?，這些成對蛋白質(zhì)之間的序列匹配度高于或等于它們中的任一蛋白質(zhì)與其它物種中所有蛋白質(zhì)間的匹配度。通過尋找存在于2個物種間的所有成對蛋白質(zhì)確定所有潛在的直系同源基因?qū)?，這些成對蛋白質(zhì)之間的序列匹配度等于或高于二者中的任一蛋白質(zhì)與另一物種中任何其它蛋白質(zhì)間的匹配度,也即該成對的2個蛋白質(zhì)互為最佳匹配蛋白質(zhì)。潛在的共直系同源基因?qū)σ餐ㄟ^尋找兩個物種間的所有成對蛋白質(zhì)加以確定，某一物種中的內(nèi)旁系同源基因通過與另一物種中的內(nèi)旁系同源基因間的直系同源關系而結(jié)為共直系同源基因?qū)?。對采用MCL程序?qū)ν椿蜻M行聚類分析，所有存在內(nèi)旁系及直系同源關系的蛋白質(zhì)被歸入同一類，共產(chǎn)生57 250個同源基因簇。去除不包括所有17個物種的同源基因簇后，共得到1 364個基因簇。通過選留至少在該17種植物的14個成員中均存在單拷貝的同源基因簇對其進一步過濾，最后獲得144個單拷貝直系同源基因簇。利用ClustalX軟件進行序列比(采用Gonnet系列矩陣，空位開放罰分=10，空位延伸罰分=0.1)，之后采用Gblocks軟件從比對結(jié)果中提取保守序列區(qū)。最后，按最大似然法并利用PhyML軟件[89]構(gòu)建系統(tǒng)發(fā)生樹。采用SPR(子樹修剪與嫁接)和NNI(最近鄰居互換)算法中最佳的一種尋找最優(yōu)樹。通過自舉分析(Bootstrapping)檢驗系統(tǒng)樹的可信度(隨機重復100次,隨機種子數(shù)目為7)。

4.11 基因的實驗確認

利用RNeasy Mini Kit (Qiagen)試劑盒并按照使用說明書從巴西橡膠無性系RRIM 600的嫩葉和膠乳中分離總RNA。之后，采用試劑盒SuperScript VILO cDNA Synthesis Kit (Invitrogen)合成cDNA第一鏈。以cDNA為模板并采用基因特異性引物對基因進行PCR擴增。將純化的PCR產(chǎn)物克隆到pCR4Blunt-TOPO (Invitrogen)或pJET1.2/blunt (Fermentas)載體中并測序。與橡膠生物合成、木質(zhì)素生物合成、病害抗性、過敏源、轉(zhuǎn)錄因子以及植物激素生物合成相關的基因均按以上步驟擴增、克隆和測序。

4.12 數(shù)據(jù)的可得性

該全基因組鳥槍法測序計劃已提交到DDBJ/EMBL/GenBank數(shù)據(jù)庫，其GenBank登錄號：AJJZ00000000。本文所描述的版本為第一版。組裝的轉(zhuǎn)錄本也已提交到NCBI轉(zhuǎn)錄組鳥槍組裝數(shù)據(jù)庫，其GenBank登錄號：JT914190-JT981478。該橡膠基因組游覽器可供下載和使用[90]。

本項目受馬來西亞高等教育部頂尖大學計劃(APEX)資助,由馬來西亞理科大學、中國南開大學泰達生物技術研究院、美國夏威夷大學、美國Ariadne Genomics公司、丹麥CLC bio公司、南非德班理工大學、美國德克薩斯A&M大學等機構(gòu)共同完成。Ahmad Yamin Abdul Radman等所有署名作者均已閱讀并同意本論文的最終稿。

[1] Prabhakaran Nair KP. The Agronomy and Economy of Important Tree Crops of the Developing World. Burlington: Elsevier， 2010.

[2] Mooibroek H et al.. Alternative sources of natural rubber. Appl Microbiol Biotechnol 2000, 53:355-365.

[3] Sakdapipanich JT. Structural characterization of natural rubber based on recent evidence from selective enzymatic treatments. J Biosci Bioeng 2007, 103:287-292.

[4] Asawatreratanakul K et al. Molecular cloning, expression and characterization of cDNA encoding cis-prenyltransferases from Hevea brasiliensis: a key factor participating in natural rubber biosynthesis. Eur J Biochem 2003, 270:4671-4680.

[5] Gronover CS et al.. Natural rubber biosynthesis and physic-chemical studies on plant derived latex. In Biotechnology of Biopolymers. Edited by Elnashar M. Croatia: Intech Open Acess Publisher; 2011:75-88.

[6] Triwitayakorn K et al.. Transcriptome sequencing of Hevea brasiliensis for development of microsatellite markers and construction of a genetic linkage map. DNA Res 2011, 18:471-482.

[7] Xia Z et al.. RNA-Seq analysis and de novo transcriptome assembly of Hevea brasiliensis. Plant Mol Biol 2011, 77: 299-308.

[8] Chow KS et al.. Metabolic routes affecting rubber biosynthesis in Hevea brasiliensis latex. J Exp Bot 2012, 63:1863-1871.

[9] Gebelin V et al.. Identification of novel microRNAs in Hevea brasiliensis and computational prediction of their targets. BMC Plant Biol 2012, 12:18.

[10] Li D et al.. De novo assembly and characterization of bark transcriptome using Illumina sequencing and development of EST-SSR markers in rubber tree (Hevea brasiliensis Muell. Arg.). BMC Genomics 2012, 13:192.

[11] Lertpanyasampatha M et al.. Genome-wide analysis of microRNAs in rubber tree (Hevea brasiliensis L.) using high-throughput sequencing. Planta 2012, 236:437-445.

[12] Pootakham W et al.. Development of genomic-derived simple sequence repeat markers in Hevea brasiliensis from 454 genome shotgun sequences. Plant Breeding 2012, 131:555-562.

[13] Leitch AR et al.. Molecular cytogenetic studies in rubber, Hevea brasiliensis Muell. Arg. (Euphorbiaceae). Genome 1998, 41: 464-467.

[14] Bennett MD et al.. Nuclear DNA amounts in angiosperms-583 new estimates. Ann Bot 1997, 80:169-196.

[15] Margulies M et al.. Genome sequencing in microfabricated high-density picolitre reactors. Nature 2005, 437:376-380

[16] Barthelson R et al.. Plantagora: modeling whole genome sequencing and assembly of plant genomes. PLoS One 2011, 6: e28436.

[17] Le Guen V et al.. A rubber tree’s durable resistance to Microcyclus ulei is conferred by a qualitative gene and a major quantitative resistance factor. Tree Genet Genomes 2011, 7: 877-889.

[18] Shulaev V et al.. The genome of woodland strawberry (Fragaria vesca). Nat Genet 2011, 43:109-116.

[19] van Bakel H et al.. The draft genome and transcriptome of Cannabis sativa. Genome Biol 2011, 12:R102.

[20] The International Barley Genome Sequencing Consortium. A physical, genetic and functional sequence assembly of the barley genome. Nature 2012, 491:711-716.

[21] Schnable PS et al.. The B73 maize genome: complexity, diversity, and dynamics. Science 2009, 326:1112-1115.

[22] Haas BJ et al.. Automated eukaryotic gene structure annotation using EVidenceModeler and the program to assemble spliced alignments. Genome Biol 2008, 9:R7.

[23] Boutet E et al.. UniProtKB/Swiss-Prot. Methods Mol Biol. 2007,406:89-112.

[24] Zdobnov EM et al.. InterProScan - an integration platform for the signature-recognition methods in InterPro. Bioinformatics 2001, 17:847-848.

[25] Kanehisa M et al.. The KEGG resource for deciphering the genome. Nucleic Acids Res 2004, 32:D277-D280.

[26] Tatusov RL et al.. The COG database: an updated version includes eukaryotes. BMC Bioinformatics 2003, 4:41.

[27] Finn RD et al.. Pfam: clans, web tools and services. Nucleic Acids Res 2006, 34:D247-D251.

[28] Tangphatsornruang S et al.. Characterization of the complete chloroplast genome of Hevea brasiliensis reveals genome rearrangement, RNA editing sites and phylogenetic relationships. Gene 2011, 475:104-112.

[29] Ashburner M et al.. Gene Ontology: tool for the unification of biology. Nat Genet 2000, 25:25-29.

[30] Silpi U et al.. Carbohydrate reserves as a competing sink: evidence from tapping rubber trees. Tree Physiol 2007, 27: 881-889.

[31] Ponciano G et al.. Transcriptome and gene expression analysis in cold-acclimated guayule (Parthenium argentatum) rubber-producing tissue. Phytochemistry 2012, 79:57-66.

[32] Tang C et al.. The sucrose transporter HbSUT3 plays an active role in sucrose loading to laticifer and rubber productivity in exploited trees of Hevea brasiliensis (para rubber tree). Plant Cell Environ 2010, 33:1708-1720.

[33] Kekwick RGO: The formation of isoprenoids in Hevea latex. In Physiology of Rubber Tree Latex. Edited by d’Auzac J, Jacob JL, Chrestin L. Boca Raton: CRC Press; 1989:145-164.

[34] Ko JH et al.. Transcriptome analysis reveals novel features of the molecular events occurring in the laticifers of Hevea brasiliensis (para rubber tree). Plant Mol Biol 2003, 53: 479 -492.

[35] Sando T et al.. Cloning and characterization of the 2-C-methyl-D-erythritol 4-phosphate (MEP) pathway genes of a natural-rubber producing plant. Hevea brasiliensis. Biosci Biotechnol Biochem 2008,72: 2903-2917.

[36] Rattanapittayaporn A et al.. Significant role of bacterial undecaprenyl diphosphate (C55-UPP) for rubber synthesis by Hevea latex enzymes. Macromol Biosci 2004, 4:1039-1052.

[37] Kharel Y et al.. Molecular analysis of cis-prenyl chain elongating enzymes. Nat Prod Rep 2003, 20:111-118.

[38] Tamura K et al.. MEGA5: molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol Biol Evol 2011, 28:2731-2739.

[39] Dillon SK et al.. Allelic variation in cell wall candidate genes affecting solid wood properties in natural populations and land races of Pinus radiata. Genetics 2010, 185:1477-1487.

[40] Richmond TA et al.. The cellulose synthase superfamily. Plant Physiol 2000, 124:495-498.

[41] Djerbi S et al.. The genome sequence of black cottonwood (Populus trichocarpa) reveals 18 conserved cellulose synthase (CesA) genes. Planta 2005, 221: 739-746.

[42] Tuskan GA et al.. The genome of black cottonwood, Populus trichocarpa (Torr. & Gray). Science 2006,313:1596-1604.

[43] Durrant WE, Dong X. Systemic acquired resistance. Annu Rev Phytopathol 2004, 42:185-209.

[44] Mun JH et al.. Genome-wide identification of NBS-encoding resistance genes in Brassica rapa.Mol Genet Genomics 2009, 282:617-631.

[45] Yeang HY et al.. Allergenic protein of natural rubber latex. Methods 2002, 27:32-45.

[46] Ryan CA et al.. Systemic wound signaling in plants: a new perception. Proc Natl Acad Sci USA 2002, 99:6519-6520.

[47] Feller A et al.. Evolutionary and comparative analysis of MYB and bHLH plant transcription factors. Plant J 2011, 6:94-116.

[48] Wang D et al.. Genome-wide analysis of CCCH zinc finger family in Arabidopsis and rice.BMC Genomics 2008, 9:44.

[49] Riano-Pachon DM et al.. Green transcription factors: a Chlamydomonas overview. Genetics 2008, 179:31-39.

[50] Rushton PJ et al.. WRKY transcription factors. Trends Plant Sci 2010, 15:247-258.

[51] Arora R et al.. MADS-box gene family in rice: genome-wide identification, organization and expression profiling during reproductive development and stress. BMC Genomics 2007, 8:242.

[52] Ming R et al.. The draft genome of the transgenic tropical fruit tree papaya (Carica papaya Linnaeus). Nature 2008, 452:991 -996.

[53] Mas P: Circadian clock signaling in Arabidopsis thaliana: from gene expression to physiology and development. Int J Dev Biol 2005, 49:491-500.

[54] Jain M et al.. F-box proteins in rice. Genome-wide analysis, classification, temporal and spatial gene expression during panicle and seed development, and regulation by light and abiotic stress. Plant Physiol 2007, 143:1467-1483.

[55] Yang X et al.. The F-box gene family is expanded in herbaceous annual plants relative to woody perennial plants. Plant Physiol 2008, 148:1189-1200.

[56] Hua Z et al.. Phylogenetic comparison of F-box (FBX) gene superfamily within the plant kingdom reveals divergent evolutionary histories indicative of genomic drift. PLoS One 2011, 6:e16219.

[57] Chaudhary N et al.. Carotenoid biosynthesis genes in rice: structural analysis, genome-wide expression profiling and phylogenetic analysis. Mol Genet Genomics 2010, 283:13-33.

[58] Yang S et al.. Repetitive element-mediated recombination as a mechanism for new gene origination in Drosophila. PLoS Genet 2008, 4:e3.

[59] Miller JR et al.. Assembly algorithms for next-generation sequencing data.Genomics 2010, 95:315-327.

[60] RepeatModeler. http://www.repeatmasker.org/RepeatModeler.html webcite

[61] Boeckmann B et al.. The SWISS-PROT protein knowledgebase and its supplement TrEMBL in 2003. Nucleic Acids Res 2003, 31: 365-370.

[62] RepeatMasker. http://www.repeatmasker.org.

[63] Li R et al.. De novo assembly of human genomes with massively parallel short read sequencing. Genome Res 2010, 20:265-272.

[64] Langmead B et al.. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome. Genome Biol 2009, 10:R25.

[65] Trapnell C et al.. TopHat: discovering splice junctions with RNA-Seq. Bioinformatics 2009, 25:1105-1111.

[66] Haas BJ et al.. Improving the Arabidopsis genome annotation using maximal transcript alignment assemblies. Nucleic Acids Res 2003, 31: 5654-5666.

[67] Slater GS et al.. Automated generation of heuristics for biological sequence comparison. BMC Bioinformatics 2005, 6:31.

[68] Dong Q et al.. Comparative plant genomics resources at PlantGDB. Plant Physiol 2005, 139:610-618.

[69] Wu TD et al.. GMAP: a genomic mapping and alignment program for mRNA and EST sequences. Bioinformatics 2005, 21: 1859-1875

[70] Huang X et al.. A tool for analyzing and annotating genomic sequences. Genomics 1997, 46:37-45.

[71] Kent WJ: BLAT-the BLAST-like alignment tool. Genome Res 2002,12:656-664.

[72] Salamov AA et al.. Ab initio gene finding in Drosophila genomic DNA. Genome Res 2000, 10:516-522.

[73] Stanke M et al.. AUGUSTUS: a web server for gene prediction in eukaryotes that allows user-defined constraints. Nucleic Acids Res 2005, 33:W465-W467.

[74] Majoros WH et al.. TigrScan and GlimmerHMM: two open source ab initio eukaryotic gene-finders. Bioinformatics 2004, 20: 2878-2879

[75] Korf I: Gene finding in novel genomes. BMC Bioinformatics 2004,5:59.

[76] Ter-Hovhannisyan V et al.. Gene prediction in novel fungal genomes using an ab initio algorithm with unsupervised training. Genome Res 2008, 18:1979-1990.

[77] Blanco E et al.. Using geneid to identify genes. In Current Protocols in Bioinformatics. Edited by Baxevanis AD, Davison DB. New York: John Wiley and Sons Inc; 2002.Unit 4.3.

[78] Szklarczyk D et al.. The STRING database in 2011: functional interaction networks of proteins, globally integrated and scored. Nucleic Acids Res 2011, 39:D561-D568.

[79] Thomas PD et al.. PANTHER: a browsable database of gene products organized by biological function, using curated protein family and subfamily classification. Nucleic Acids Res 2003, 31:334-341.

[80] Attwood TK et al.. The PRINTS protein fingerprint database: functional and evolutionary applications. In Encyclopaedia of Genetics, Genomics, Proteomics and Bioinformatics. Edited by Dunn M, Jorde L, Little P, Subramaniam A. New Jersey: John Wiley and Sons Ltd; 2006.

[81] Sigrist CJA et al.. PROSITE, a protein domain database for functional characterization and annotation. Nucleic Acids Res 2010, 38:D161-D166.

[82] Letunic I et al.. SMART 7. recent updates to the protein domain annotation resource. Nucleic Acids Res 2012, 40:D303-D305.

[83] Lowe TM et al.. tRNAscan-SE: a program for improved detection of transfer RNA genes in genomic sequence. Nucleic Acids Res 1997, 25:955-964.

[84] Altschul SF et al.. Gapped BLAST and PSI-BLAST: a new generation of protein database search programs. Nucleic Acids Res 1997, 25: 3389-3402.

[85] Petersen TN et al.. SignalP 4.0: discriminating signal peptides from transmembrane regions. Nat Methods 2011, 8:785-786.

[86] Lupas A et al.. Predicting coiled coils from protein sequences. Science 1991, 252:1162-1164.

[87] Selkov E Jr et al.. MPW: the metabolic pathways database. Nucleic Acids Res 1998, 26:43-45

[88] Li L et al.. OrthoMCL. identification of ortholog groups for eukaryotic genomes. Genome Res 2003, 13: 2178-2189.

[89] Guindon S et al.. New algorithms and methods to estimate maximum-likelihood phylogenies: assessing the performance of PhyML 3.0. Syst Biol 2010, 59:307-321.

[90] Rubber Genome Browser. http://bioinfo.ccbusm.edu.my/cgi-bin/gb2/gbrowse/Rubber/

(全文譯自BMC Genomics 2013，14:75)

Draft genome sequence of the rubber treeHeveabrasiliensis

Translated by TIAN Lang

RubberResearchInstitute,ChineseAcademyofTropicalAgriculturalSciences,Danzhou,Hainan571737

Heveabrasiliensis,a member of the Euphorbiaceae family,is the major commercial source of natural rubber (NR).NR is a latex polymer with high elasticity,flexibility,and resilience that has played a critical role in the world economy since 1876.Here,we report the draft genome sequence of H.brasiliensis.The assembly spans～1.1 Gb of the estimated 2.15 Gb haploid genome.Overall,78% of the genome was identified as repetitive DNA.Gene prediction shows 68,955 gene models, of which 12.7% are unique to Hevea.Most of the key genes associated with rubber biosynthesis,rubberwood formation,disease resistance,and allergenicity have been identified.The knowledge gained from this genome sequence will aid in the future development of high-yielding clones to keep up with the ever increasing need for natural rubber.

Heveabrasiliensis; Euphorbiaceae; Natural rubber; Genome

2016-04-16

田 郎(1961-)，男，侗族，湖南新晃人，碩士，副研究員，現(xiàn)從事植物組織培養(yǎng)及分子生物學研究工作。

S

A

1001-2117(2017)02-0101-15