瀕危樹(shù)種香木蓮轉(zhuǎn)錄組分析

苗藝明，石 松，楊 梅，劉世男

(1.廣西大學(xué)林學(xué)院，廣西 南寧 530004；2.廣西壯族自治區(qū)都安瑤族自治縣自然資源局，廣西 都安 530799)

香木蓮(Manglietiaaromatica)是木蘭科木蓮屬常綠喬木，屬?lài)?guó)家Ⅱ級(jí)重點(diǎn)保護(hù)植物.該樹(shù)種生長(zhǎng)較快，高可達(dá)35 m，胸徑可達(dá)1.9 m，主要分布于廣西、云南和貴州海拔400～1 600 m地區(qū).香木蓮具有較強(qiáng)的適應(yīng)性，在石灰?guī)r發(fā)育的鈣質(zhì)土，以及土壤貧瘠、生境惡劣的溶洞邊緣或石漠化嚴(yán)重的山坡能形成小群落或單株生長(zhǎng).香木蓮全株都具香氣且均可提取香油，用于調(diào)配名貴香料[1]；樹(shù)體形態(tài)優(yōu)美，花大香艷美麗[2]；木材較輕軟、紋理細(xì)、抗蟲(chóng)蛀、耐腐蝕、不裂、不變形[3].香木蓮是重要的香料，是園林綠化觀賞及優(yōu)良的用材樹(shù)種.長(zhǎng)期以來(lái)深受當(dāng)?shù)厝罕娤矏?ài)，其野生資源幾乎被砍伐殆盡，現(xiàn)存野生種群十分稀少，且多為老樹(shù)，幼樹(shù)、幼苗少，自然更新能力差.PAN等[4]揭示了香木蓮的瀕危原因，結(jié)果顯示，其花粉萌發(fā)能力低且傳粉條件較差，大孢子、雌配子退化和敗育，這些都大大降低了香木蓮的結(jié)實(shí)率.種苗繁育研究對(duì)瀕危物種擴(kuò)大種群數(shù)量具有重要意義.有性繁殖研究發(fā)現(xiàn)，香木蓮種皮含油質(zhì)，致使種子發(fā)芽率和保存率低，從而限制了種子苗來(lái)源[5]；無(wú)性繁殖主要是扦插繁育，ABT生根粉能提高扦插生根率[6].研究幼苗在低溫和干旱脅迫下的生理特性發(fā)現(xiàn)，香木蓮具有一定的抗寒能力，但耐干旱能力較弱[7-8].目前，有關(guān)香木蓮種苗繁育和瀕危原因的研究只有少量報(bào)道，而分子生物學(xué)方面的研究未見(jiàn)報(bào)道.本研究在高通量測(cè)序的基礎(chǔ)上，利用生物信息學(xué)方法分析香木蓮轉(zhuǎn)錄組序列信息，為后續(xù)的香木蓮遺傳多樣性、種質(zhì)資源保存利用和功能基因挖掘等提供基礎(chǔ)數(shù)據(jù).

1 材料與方法

1.1 試驗(yàn)材料

本次研究的試驗(yàn)材料采自廣西壯族自治區(qū)百色市都安縣三只羊鄉(xiāng)上遠(yuǎn)村.收集香木蓮嫩葉、嫩枝和花，用錫箔紙包裹迅速放入液氮中保存，帶回實(shí)驗(yàn)室，置于-80 ℃冰箱備用.

1.2 文庫(kù)構(gòu)建及測(cè)序

使用試劑盒提取嫩葉、嫩枝和花樣本的總RNA，質(zhì)量檢驗(yàn)合格后，等量混合不同組織的RNA樣品.香木蓮cDNA文庫(kù)構(gòu)建及測(cè)序由杭州科睿迪有限責(zé)任公司完成.

1.3 轉(zhuǎn)錄組序列組裝、注釋及分析

RNA-seq測(cè)序完成后，過(guò)濾帶有接頭、低質(zhì)量、冗余的Raw reads得到高質(zhì)量clean reads數(shù)據(jù)，統(tǒng)計(jì)數(shù)量、長(zhǎng)度、N%、Q20及GC%.利用Trinity software進(jìn)行denovo組裝，從長(zhǎng)到短排序，依次累加不小于總長(zhǎng)50%(N50)的拼接轉(zhuǎn)錄本長(zhǎng)度，統(tǒng)計(jì)各文庫(kù)subreads.將獲得的Unigenes分別與NR、GO、KOG、KEGG、SwissProt等數(shù)據(jù)庫(kù)進(jìn)行比對(duì)，分析相應(yīng)的功能注釋.

2 結(jié)果與分析

2.1 香木蓮RNA-seq與組裝

通過(guò)轉(zhuǎn)錄組測(cè)序，共得到36 737 304條原始序列，過(guò)濾處理后獲得35 321 846條有效序列，占總數(shù)的96.15%，Q20、Q30序列分別占總數(shù)的98.29%和95.07%，GC含量占總數(shù)的48.27%，堿基錯(cuò)誤率為0.02%.以上結(jié)果說(shuō)明，通過(guò)高通量測(cè)序平臺(tái)獲得的香木蓮序列數(shù)量和質(zhì)量較高，可以用于后續(xù)的相關(guān)生物信息學(xué)分析.

利用Trinity軟件組裝處理后的片段，共獲得48 123條Unigene，全部堿基數(shù)達(dá)46 188 480 bp.組裝后的Unigene長(zhǎng)度分布見(jiàn)圖1.結(jié)果可知：Unigene平均長(zhǎng)度為960 nt，N50為1 331 nt.序列長(zhǎng)度為200～500 nt的有15 932條，占33.1%；500～1 000 nt的有15 089條，占31.4%；1 000～1 500 nt的有8 098條，占16.85%；1 500～2 000 nt的有4 553條，占9.5%；大于等于2 000 nt的有4 451條，占9.2%.

2.2 香木蓮轉(zhuǎn)錄組Unigene功能注釋

利用Blast軟件將香木蓮48 123條Unigene與各數(shù)據(jù)庫(kù)進(jìn)行比對(duì)，結(jié)果顯示：比對(duì)到NR、GO、SwissProt、KEGG和KOG數(shù)據(jù)庫(kù)的Unigene分別有37 877、32 125、27 988、30 143和37 199條，占比依次為78.7%、66.8%、58.2%、62.6%和77.3%.

2.2.1 香木蓮轉(zhuǎn)錄組Unigene的NR功能注釋

利用Blast軟件將香木蓮全部Unigene與NR數(shù)據(jù)庫(kù)進(jìn)行比對(duì)，結(jié)果見(jiàn)圖2.由圖2可見(jiàn)：在匹配的近緣物種中，香木蓮與蓮花(Nelumbonucifera)同源序列最多，為10 427條，占總數(shù)的27.5%；其次為博落回(Macleayacordata)，有6 084條，占總數(shù)的16.1%；葡萄(Vitisvinifera) 2 287條、棕櫚(Elaeisguineensis) 1 993條、海棗(Phoenixdactylifera)1 622條、洛磯山耬斗菜(Aquilegiacoerulea)1 465條、無(wú)油樟(Amborellatrichopoda)854條、核桃(Juglansregia) 553條、栓皮櫧(Quercussuber)521條、小果野芭蕉(Musaacuminatasubsp.malaccensis)516條、蘆筍(Asparagusofficinalis) 501條、菠蘿(Ananascomosus) 435條、橡膠樹(shù)(Heveabrsiliensis)427條、甜橙(Citrussinensis) 346條，分別占總數(shù)的6.0%、5.3%、4.3%、3.9%、2.3%、1.5%、1.4%、1.4%、1.3%、1.2%、1.1%和0.9%；其余9 846條分布于其他物種中，占總數(shù)的26%.

圖1香木蓮轉(zhuǎn)錄組組裝序列長(zhǎng)度分布Fig.1Length distribution of transcriptome Unigenes for Manglietia aromatica

圖2香木蓮轉(zhuǎn)錄組Unigene注釋匹配的物種分布Fig.2Species distribution of Manglietia aromatic with Unigenes annotation

2.2.2 香木蓮轉(zhuǎn)錄組Unigene的GO功能注釋

香木蓮轉(zhuǎn)錄組Unigene的GO功能注釋見(jiàn)圖3.結(jié)果表明：32 125條Unigene共獲得246 974個(gè)GO功能注釋?zhuān)殖缮飳W(xué)過(guò)程、分子功能以及細(xì)胞組分3大類(lèi).其中，生物學(xué)過(guò)程獲得了最多的注釋?zhuān)?07 579個(gè)，占43.6%；其次是細(xì)胞組分，獲得98 271個(gè)注釋?zhuān)?9.8%；第3為分子功能，獲得41 124個(gè)注釋?zhuān)?6.6%.3大功能進(jìn)一步又分成58個(gè)亞類(lèi)，其中，生物學(xué)過(guò)程包含28個(gè)亞類(lèi)，細(xì)胞過(guò)程、代謝過(guò)程亞類(lèi)獲得的注釋偏多，分別占該類(lèi)型的20.0%和17.0%，碳利用率、細(xì)胞失活、行為以及硫利用率亞類(lèi)注釋最少，均占該類(lèi)型的0.01%及以下；細(xì)胞過(guò)程包含18個(gè)亞類(lèi)，其中，細(xì)胞、細(xì)胞組分亞類(lèi)獲得的注釋偏多，均占該類(lèi)型的22.2%；分子功能包含12個(gè)亞類(lèi)，其中，結(jié)合、催化活性亞類(lèi)獲得的注釋偏多，分別占該類(lèi)型的47.1%和39.3%，翻譯調(diào)節(jié)器活性和蛋白標(biāo)亞類(lèi)獲得的注釋較少，分別占該類(lèi)型的0.02%和0.01%.

圖3香木蓮轉(zhuǎn)錄組Unigene的GO功能分類(lèi)Fig.3GO functional annotation of transcriptome Unigenes for Manglietia aromatica

2.2.3 香木蓮轉(zhuǎn)錄組Unigene的KOG功能注釋

香木蓮轉(zhuǎn)錄組Unigene序列KOG蛋白數(shù)據(jù)庫(kù)分類(lèi)注釋見(jiàn)圖4.結(jié)果顯示：有37 199條Unigenes能夠匹配在KOG數(shù)據(jù)庫(kù)中，共獲得25 525條注釋?zhuān)煞譃?5個(gè)功能大類(lèi).其中，一般功能預(yù)測(cè)基因最多，有4 542條，占總數(shù)的17.8%；其次是信號(hào)傳導(dǎo)機(jī)制，有2 999條，占總數(shù)的11.8%；再次是翻譯后修飾、蛋白質(zhì)轉(zhuǎn)換、伴侶有2 492條，占9.8%；轉(zhuǎn)錄功能有1 412條，占5.5%；碳水化合物運(yùn)輸和代謝有1 314條，占5.2%；細(xì)胞內(nèi)、分泌、囊泡運(yùn)輸有1 293條，占5.1%；胞外結(jié)構(gòu)和細(xì)胞運(yùn)動(dòng)所占比例最少，分別占0.3%和0.05%.

圖4香木蓮轉(zhuǎn)錄組Unigene的KOG注釋分類(lèi)Fig.4KOG functional classification of transcriptome Unigenes for Manglietia aromatica

2.2.4 香木蓮轉(zhuǎn)錄組Unigene的KEGG代謝通路分析

香木蓮轉(zhuǎn)錄組Unigene的KEGG分類(lèi)見(jiàn)圖5.結(jié)果顯示：共有30 143條Unigenes獲得注釋并涉及142個(gè)代謝通路.進(jìn)一步將這142個(gè)代謝通路劃分為5大類(lèi)，包括代謝、遺傳信息處理、環(huán)境信息處理、細(xì)胞過(guò)程以及生物系統(tǒng)相關(guān)通路.該5大類(lèi)又分為19個(gè)亞類(lèi)，其中，代謝相關(guān)通路中11個(gè)亞類(lèi)，以全局和概述地圖居多，占總數(shù)的33.3%，第2是碳水化合物代謝相關(guān)通路，占該通路的17.5%，萜類(lèi)和聚酮代謝、多糖合成和代謝以及核苷酸代謝相關(guān)通路相對(duì)較少，占比均在3.0%以下；遺傳信息處理相關(guān)通路有4個(gè)亞類(lèi)，其中，翻譯相關(guān)通路所占比例最大，占總數(shù)的38.7%，折疊、分選和降解處理相關(guān)通路占32.3%，復(fù)制和修復(fù)相關(guān)通路最少，占比僅為13.0%；環(huán)境處理相關(guān)代謝通路中包含2個(gè)亞類(lèi)，信號(hào)傳導(dǎo)通路明顯居多，高達(dá)82.3%；細(xì)胞過(guò)程和生物系統(tǒng)相關(guān)通路都僅包括1個(gè)亞類(lèi).

圖5香木蓮轉(zhuǎn)錄組Unigene的KEGG分類(lèi)Fig.5KEGG classification of transcriptome Unigenes for Manglietia aromatica

3 小結(jié)與討論

轉(zhuǎn)錄組測(cè)序技術(shù)在植物分子標(biāo)記開(kāi)發(fā)、功能基因挖掘及鑒定研究中起著重要作用[9-10]，已在木蘭科樹(shù)種中廣泛應(yīng)用.已有研究利用RNA-Seq技術(shù)對(duì)多種植物進(jìn)行了測(cè)序，結(jié)果顯示：紅花玉蘭(Magnoliawufengesis)共獲得94 805條Unigene，平均長(zhǎng)度為695 nt，N50為1 038 nt[11]；景寧玉蘭(M.sinostellata)獲得52 441條Unigene，平均長(zhǎng)度為648 nt，N50為1 126 nt[12]；樂(lè)東擬單性木蘭(Parakmerialotugensis)獲得273 252條Unigene，平均長(zhǎng)度為590 nt，N50為752 nt[13]；鵝掌楸(LiriodendronchinenseSarg)獲得162 092條Unigene，平均長(zhǎng)度為547 nt，N50為719 nt[14].本研究利用IIIumina HiSeqTM4000對(duì)香木蓮轉(zhuǎn)錄組進(jìn)行了測(cè)序，共獲得48 123條Unigene，平均長(zhǎng)度為960 nt，N50為1 331 nt.香木蓮的N50和平均長(zhǎng)度都高于上述木蘭科樹(shù)種.N50是評(píng)價(jià)組裝序列完整性的重要指標(biāo)，長(zhǎng)度越長(zhǎng)，代表組裝的完整性越好[15].研究結(jié)果表明，香木蓮測(cè)序獲得的序列質(zhì)量高且拼接完整性較好，有利于后續(xù)開(kāi)展基因組方面的研究.

通過(guò)幾種木蘭科樹(shù)種轉(zhuǎn)錄組Unigene與NR數(shù)據(jù)庫(kù)比對(duì)發(fā)現(xiàn)，紅花玉蘭Unigene獲得注釋數(shù)較高的前3個(gè)物種為蓮花、葡萄和可可樹(shù)(Thenobromacacao)[11]；華木蓮獲得注釋較高的前3個(gè)物種為蓮花、油棕和海棗[16]；樂(lè)東擬單性木蘭和景寧木蘭獲得注釋數(shù)較高的前3個(gè)物種均為葡萄、海棗和可可樹(shù)[12-13].而香木蓮獲得匹配率較高的前3個(gè)物種為蓮花、博落回和葡萄，且博落回是以上樹(shù)種都沒(méi)有匹配到的.由此表明，香木蓮與同科的樹(shù)種存在相似功能基因，也可能存在特異功能基因.由香木蓮GO功能注釋分析結(jié)果可知，被注釋的32 125條Unigene在功能上劃分為3大類(lèi)58個(gè)亞類(lèi)，其中，注釋到生物學(xué)過(guò)程大類(lèi)的Unigene數(shù)量最多，主要是細(xì)胞過(guò)程和代謝過(guò)程；通過(guò)KEGG pathway分析可知，有30 143條Unigene獲得注釋?zhuān)婕?大類(lèi)19個(gè)亞類(lèi)共142個(gè)代謝通路，其中，以代謝相關(guān)通路和碳水化合物代謝相關(guān)通路為主.這一結(jié)果與華木蓮、紅花玉蘭的Unigene GO功能注釋和KEGG pathway分析結(jié)果類(lèi)似[11，16]，都以生物學(xué)過(guò)程中的細(xì)胞過(guò)程和代謝過(guò)程為主.可見(jiàn)，香木蓮在細(xì)胞和代謝活動(dòng)的基因表達(dá)量相對(duì)較高，具有較強(qiáng)的代謝能力和豐富的生物過(guò)程.此外，分析玉蘭的KEGG pathway發(fā)現(xiàn)，通路“類(lèi)黃酮合成”富集44個(gè)Unigene，通路“花青素合成”富集9個(gè)Unigene，通路“黃酮和黃酮醇”富集6個(gè)Unigene[11]，這些通路可能參與花青素苷合成，進(jìn)一步影響紅花玉蘭花色.本研究中的KEGG pathway分析表明，通路“類(lèi)黃酮合成” “花青素合成” “黃酮和黃酮醇”分別富集140、4和19個(gè)Unigene.前期的野外調(diào)查發(fā)現(xiàn)，不同居群香木蓮的花被片顏色有淡紅色和白色兩種，因此，通過(guò)測(cè)序分析香木蓮轉(zhuǎn)錄組，有利于挖掘花青素苷合成的相關(guān)通路及關(guān)鍵基因，可為闡明不同花色形成機(jī)理提供重要的研究基礎(chǔ).