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

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

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

1 材料與方法

1.1 試驗材料

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

1.2 文庫構建及測序

使用試劑盒提取嫩葉、嫩枝和花樣本的總RNA，質量檢驗合格后，等量混合不同組織的RNA樣品.香木蓮cDNA文庫構建及測序由杭州科睿迪有限責任公司完成.

1.3 轉錄組序列組裝、注釋及分析

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

2 結果與分析

2.1 香木蓮RNA-seq與組裝

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

利用Trinity軟件組裝處理后的片段，共獲得48 123條Unigene，全部堿基數達46 188 480 bp.組裝后的Unigene長度分布見圖1.結果可知：Unigene平均長度為960 nt，N50為1 331 nt.序列長度為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 香木蓮轉錄組Unigene功能注釋

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

2.2.1 香木蓮轉錄組Unigene的NR功能注釋

利用Blast軟件將香木蓮全部Unigene與NR數據庫進行比對，結果見圖2.由圖2可見：在匹配的近緣物種中，香木蓮與蓮花(Nelumbonucifera)同源序列最多，為10 427條，占總數的27.5%；其次為博落回(Macleayacordata)，有6 084條，占總數的16.1%；葡萄(Vitisvinifera) 2 287條、棕櫚(Elaeisguineensis) 1 993條、海棗(Phoenixdactylifera)1 622條、洛磯山耬斗菜(Aquilegiacoerulea)1 465條、無油樟(Amborellatrichopoda)854條、核桃(Juglansregia) 553條、栓皮櫧(Quercussuber)521條、小果野芭蕉(Musaacuminatasubsp.malaccensis)516條、蘆筍(Asparagusofficinalis) 501條、菠蘿(Ananascomosus) 435條、橡膠樹(Heveabrsiliensis)427條、甜橙(Citrussinensis) 346條，分別占總數的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條分布于其他物種中，占總數的26%.

圖1香木蓮轉錄組組裝序列長度分布Fig.1Length distribution of transcriptome Unigenes for Manglietia aromatica

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

2.2.2 香木蓮轉錄組Unigene的GO功能注釋

香木蓮轉錄組Unigene的GO功能注釋見圖3.結果表明：32 125條Unigene共獲得246 974個GO功能注釋，分成生物學過程、分子功能以及細胞組分3大類.其中，生物學過程獲得了最多的注釋，有107 579個，占43.6%；其次是細胞組分，獲得98 271個注釋，占39.8%；第3為分子功能，獲得41 124個注釋，占16.6%.3大功能進一步又分成58個亞類，其中，生物學過程包含28個亞類，細胞過程、代謝過程亞類獲得的注釋偏多，分別占該類型的20.0%和17.0%，碳利用率、細胞失活、行為以及硫利用率亞類注釋最少，均占該類型的0.01%及以下；細胞過程包含18個亞類，其中，細胞、細胞組分亞類獲得的注釋偏多，均占該類型的22.2%；分子功能包含12個亞類，其中，結合、催化活性亞類獲得的注釋偏多，分別占該類型的47.1%和39.3%，翻譯調節(jié)器活性和蛋白標亞類獲得的注釋較少，分別占該類型的0.02%和0.01%.

圖3香木蓮轉錄組Unigene的GO功能分類Fig.3GO functional annotation of transcriptome Unigenes for Manglietia aromatica

2.2.3 香木蓮轉錄組Unigene的KOG功能注釋

香木蓮轉錄組Unigene序列KOG蛋白數據庫分類注釋見圖4.結果顯示：有37 199條Unigenes能夠匹配在KOG數據庫中，共獲得25 525條注釋，可分為25個功能大類.其中，一般功能預測基因最多，有4 542條，占總數的17.8%；其次是信號傳導機制，有2 999條，占總數的11.8%；再次是翻譯后修飾、蛋白質轉換、伴侶有2 492條，占9.8%；轉錄功能有1 412條，占5.5%；碳水化合物運輸和代謝有1 314條，占5.2%；細胞內、分泌、囊泡運輸有1 293條，占5.1%；胞外結構和細胞運動所占比例最少，分別占0.3%和0.05%.

圖4香木蓮轉錄組Unigene的KOG注釋分類Fig.4KOG functional classification of transcriptome Unigenes for Manglietia aromatica

2.2.4 香木蓮轉錄組Unigene的KEGG代謝通路分析

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

圖5香木蓮轉錄組Unigene的KEGG分類Fig.5KEGG classification of transcriptome Unigenes for Manglietia aromatica

3 小結與討論

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

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