基于高通量測序的極小種群野生植物長梗杜鵑轉(zhuǎn)錄組分析

李太強 劉雄芳 萬友名 李正紅 起國海 李鈺瑩 劉秀賢 和 銳 馬 艷 馬 宏

(中國林業(yè)科學(xué)研究院資源昆蟲研究所，昆明 650233)

基于高通量測序的極小種群野生植物長梗杜鵑轉(zhuǎn)錄組分析

李太強 劉雄芳 萬友名 李正紅 起國海 李鈺瑩 劉秀賢 和 銳 馬 艷 馬 宏*

(中國林業(yè)科學(xué)研究院資源昆蟲研究所，昆明 650233)

為加強特有瀕危植物長梗杜鵑資源的評價、保護和鑒定工作，及為今后遺傳育種和改善其農(nóng)藝性狀提供有益參考。本研究采用新一代高通量測序平臺Illumina HiSeq 4000對其轉(zhuǎn)錄組測序，得到的數(shù)據(jù)過濾后進行de novo組裝并聚類去冗余，獲得74 092個Unigenes，平均長度、N50、Q20、Q30以及GC含量分別為938 nt、1 616 nt、98.22%、95.20%和43.24%，其中1 Kb以上的Unigenes有23 879條。通過與七大功能數(shù)據(jù)庫比對，分別有39 876(NR:53.82%)、38 065(NT:51.38%)、27 384(Swissprot:36.96%)、16 099(COG:21.73%)、30 401(KEGG:41.03%)、17 518(GO:23.64%)以及29 676(Interpro:40.05%)條Unigenes獲得功能注釋。長梗杜鵑轉(zhuǎn)錄組中的Unigenes根據(jù)GO功能大致可分為生物學(xué)過程、細胞組分和分子功能3大類56亞類，其中執(zhí)行生物學(xué)過程的基因最多，占41.53%；與COG數(shù)據(jù)庫比對，根據(jù)其功能大致可分為25類；以KEGG數(shù)據(jù)庫為參考，可歸為6個代謝通路大類、32類代謝途徑，并發(fā)掘出176條與人類疾病相關(guān)的Unigenes，包括內(nèi)分泌及代謝疾病(167條)和抗藥性(9條)。根據(jù)注釋結(jié)果共檢測出39 418個CDS，未注釋上的Unigenes使用ESTScan預(yù)測后獲得3 194個CDS。同時,預(yù)出到1 488個編碼TF的Unigenes，以及檢測到57 927個SNP多態(tài)位點。該轉(zhuǎn)錄組分析為今后長梗杜鵑乃至杜鵑屬植物功能基因挖掘與利用、基因克隆、抗性機理分析、遺傳資源分類和進化、分子標(biāo)記開發(fā)以及分子輔助育種等研究提供了基礎(chǔ)數(shù)據(jù)和重要參考。

長梗杜鵑；轉(zhuǎn)錄組；Unigene；功能注釋；編碼序列；轉(zhuǎn)錄因子；單核苷酸多態(tài)性

我國具有世界上最豐富的杜鵑花野生資源，目前記錄有580余種(包括Flora of China出版后被描述的一些新種)。因此，歐美至今流傳著一句名言“沒有中國的杜鵑花,就沒有西方的園林”。以長江為界，長江以南種類最多，長江以北種類很少；云南最多，四川次之，西藏第三，離分布中心愈遠，種類愈少[1～3]。

杜鵑花是世界著名觀賞花卉，年產(chǎn)量和產(chǎn)值巨大。然而，相較于種類眾多的野生資源，我國在杜鵑花品種選育方面明顯落后，市場上常見的品種均為比利時、美國、英國、日本、德國等國家選育。近年來，我國的杜鵑花育種工作取得一定進展，培育了一些具有自主知識產(chǎn)權(quán)的新品種[4～5]。另一方面，由于發(fā)展旅游、商業(yè)開發(fā)等人為破壞的影響，不少野生杜鵑花種類的生存受到不同程度的威脅，其中不乏極度瀕危的種類。因此，對杜鵑花野生資源中瀕危且具有較高觀賞價值的種類優(yōu)先開展相關(guān)研究，已迫在眉睫。長梗杜鵑(Rhododendronlongipedicellatum)就是其中的典型代表。

長梗杜鵑系杜鵑屬、杜鵑亞屬(subg.Rhododendron)、越桔杜鵑組(sect.Vireya)、類越桔杜鵑亞組(subsect.Pseudovireya)常綠植物。其花朵質(zhì)地厚實、花色鮮黃，比同屬的滇越杜鵑(R.rushforthii)、羊躑躅(R.molle)、鮮黃杜鵑(R.xanthostephanum)、純黃杜鵑(R.chrysodoron)、硫磺杜鵑(R.sulfureum)等花色更為純正。國際上，杜鵑花的花色育種趨勢為純色花，特別是黃色和藍色更為珍貴[6]。更令人稱奇的是，與野生杜鵑花期多集中于3～6月不同，長梗杜鵑的自然花期為11月下旬至翌年的2月上旬。此外，本研究組近幾年通過多次專門的野外調(diào)查，僅發(fā)現(xiàn)位于云南省麻栗坡縣的4個野生居群，各居群植株數(shù)為80～350株，海拔介于1 183～1 316 m。因其生境狹窄、分布區(qū)單一、居群數(shù)量少，受人為干擾嚴重，加之種群內(nèi)實生苗較少，自然更新較弱，為典型的極小種群野生植物(Plant Species with Extremely Small Populations,簡稱PSESP)，野生資源狀況令人擔(dān)憂，亟待保護。

鑒于此，本文以長梗杜鵑為研究材料，采用Illumina HiSeq 4000高通量測序技術(shù)先對其進行轉(zhuǎn)錄組測序、數(shù)據(jù)過濾與組裝，再用生物信息學(xué)的方法對得到的大量Unigenes進行功能注釋、代謝途徑、CDS預(yù)測及SNP檢測等分析，以期了解長梗杜鵑在一定時期的基因表達情況和功能總體特征，為轉(zhuǎn)錄組水平上的研究以及分子層面揭示其瀕危機制奠定基礎(chǔ)，亦為今后合理開發(fā)利用這一珍稀杜鵑花種類提供參考。

1 材料與方法

1.1 試驗材料

試驗材料取自云南省麻栗坡縣(104°93′E，23°16′N)，選擇5株距離15 m以上的長梗杜鵑樣株，各株采集質(zhì)量大體相同(5～6片)的當(dāng)年生幼嫩葉片，迅速放入液氮中，帶回實驗室后于-80℃冰箱中保存?zhèn)溆谩?/p>

1.2 長梗杜鵑葉片總RNA提取

各單株分別稱取0.1 g葉片等量混勻后，在液氮中迅速研磨至粉末狀，用RNeasy Plant Mini Kit(Qiagen,Vallencia,CA)試劑盒按照說明書的操作指南提取，提取出的總RNA使用DNaseⅠ將其DNA消化，樣品濃度、RNA完整性、28S/18S用Agilent 2100檢測，OD260/280用NanoDrop檢測，其RIN=8.3>8，28S/18S=1.84>1.8，OD260/280=1.97>1.8，所提取RNA質(zhì)量較高，完整性較好。

1.3 轉(zhuǎn)錄組測序

進行mRNA的富集、cDNA合成、加堿基“A”、連接接頭和PCR擴增等步驟，構(gòu)建cDNA文庫；構(gòu)建好的文庫用Agilent 2100 Bioanalyzer和ABI StepOnePlus Real-Time PCR System質(zhì)檢，合格后使用Illumina HiSeq 4000平臺進行測序。

1.4 數(shù)據(jù)過濾和de novo組裝

為保證結(jié)果的可靠性，先對測序的原始數(shù)據(jù)進行過濾，包括去除包含接頭的reads、未知堿基N含量大于5%的reads以及質(zhì)量值低于15的堿基占該reads總堿基數(shù)的比例大于20%的reads，過濾后的reads稱為“Clean Reads”。然后使用Trinity對clean reads進行de novo組裝，包含3個獨立模塊：Inchworm(構(gòu)建k-mer庫，K=25。以k-mer間overlap長度等于k-1對種子進行延伸，直到不能再延伸，形成線性Contig)；Chrysalis(把可能存在可變剪切及其它平行基因的Contigs聚類，用reads驗證，看每個Component的reads支持情況)以及Butterfly(合并在de Bruijn圖中有連續(xù)節(jié)點的線性路徑，以形成更長的序列)，Trinity的組裝結(jié)果稱為轉(zhuǎn)錄本。最后使用Tgicl進行聚類去冗余得到Unigenes。

1.5 Unigene注釋、CDS、TF編碼能力預(yù)測以及SNP檢測

使用Blast[7]對Unigenes進行NT、NR、COG、KEGG以及SwissProt注釋，使用Blast2GO[8]以及NR注釋結(jié)果進行GO注釋，使用InterProScan5[9]進行InterPro注釋。根據(jù)功能注釋結(jié)果，按照NR，SwissProt，KEGG，COG的數(shù)據(jù)庫優(yōu)先順序，挑選Unigenes的最佳比對片段作為該Unigenes的CDS。未能注釋上的Unigenes使用上一步預(yù)測的CDS作為模型進行建模，然后通過ESTScan[10]進行預(yù)測。使用getorf[11]檢測Unigenes的開放閱讀框(ORF)，在通過hmmsearch[12]將ORF比對到轉(zhuǎn)錄因子蛋白結(jié)構(gòu)域(數(shù)據(jù)來源于PlantfDB)，然后根據(jù)PlantfDB描述的轉(zhuǎn)錄因子家族特征對Unigenes進行能力鑒定。最后使用HISAT[13]把clean reads比對到Unigenes，借助GATK[14]檢測SNP。

2 結(jié)果與分析

2.1 長梗杜鵑轉(zhuǎn)錄組測序數(shù)據(jù)過濾

測序共獲得58.30 Mb的Raw Reads，過濾后得到44.85 Mb的Clean Reads，占原始讀長的76.93%，包含6.73 Gb Clean Bases。其中，Q20為98.22%，即質(zhì)量值大于20的堿基數(shù)目有6.61 Gb；Q30為95.20%。從Clean reads的堿基含量分布圖可以看出(圖1:A)，除RNA-Seq測序?qū)е虑? bp出現(xiàn)正常波動外，其余reads每個位置的堿基含量分布穩(wěn)定，無AT或GC分離現(xiàn)象。其次，在Clean reads的堿基質(zhì)量分布圖中(圖1:B)，顏色深的點質(zhì)量值集中于30～40，說明達到該質(zhì)量值范圍的堿基數(shù)目較多，而低質(zhì)量(質(zhì)量值小于20)的堿基比例較低，表明測序質(zhì)量較好，測序所得reads準(zhǔn)確性較高。

圖1 Clean reads的堿基含量和質(zhì)量分布圖Fig.1 Distribution of base composition and quality on clean reads

2.2 長梗杜鵑轉(zhuǎn)錄組de novo組裝

通過組裝得到94 906個轉(zhuǎn)錄本，組裝的連續(xù)性N50較大，為1 470 nt，說明組裝效果較好；平均長度、GC含量分別為864 nt和43.31%。進一步聚類去冗余得到N50、平均長度以及GC含量分別為1 616 nt、938 nt和43.24%的Unigenes 74 092條，序列信息達69 505 225 nt，其中clusters(同一個clusters里面有若干條相似度大于70%的Unigenes)為51 505條，singletons(單獨的Unigenes)為22 587條。從組裝的轉(zhuǎn)錄本長度分布來看(圖2:A)，隨著轉(zhuǎn)錄本長度的增加，對應(yīng)的數(shù)目呈遞減的趨勢；且大多數(shù)集中于200～2 000 nt，占轉(zhuǎn)錄本總數(shù)的89.85%，其中≥1 000 nt的占28.51%，僅3.36%的轉(zhuǎn)錄本≥3 000 nt。對聚類后Unigenes的長度分布進行統(tǒng)計(圖2:B)，有32.23%的序列長度大于1 000 nt，3 000 nt及以上的序列僅占4.12%。表明聚類去冗余的效果較好，所得Unigenes平均長度、N50以及1 000 nt以上的序列數(shù)都有所提高，對Trinity組裝結(jié)果起到了很好的優(yōu)化作用。

圖2 長梗杜鵑轉(zhuǎn)錄本和Unigene的長度分布圖Fig.2 Transcripts and Unigenes length distribution for R.longipedicellatum

2.3 長梗杜鵑轉(zhuǎn)錄組Unigene的NR功能注釋

將去冗余后的74 092條Unigenes比對到NR等七大功能數(shù)據(jù)庫，最終被任一數(shù)據(jù)庫注釋上的序列共有45 538條，注釋率為61.46%。其中，在NR庫中(E≤1e-5)獲得注釋的Unigenes數(shù)量最多，有39 876條(占Unigenes總數(shù)的53.82%)，且與其它物種已知基因序列具有不同程度的匹配性。如注釋物種分布圖(圖3)所示：注釋為葡萄(Vitisvinifera)相關(guān)基因的序列最多，達10 261條(占NR庫中被注釋Unigenes的25.73%)，其次，序列同源性大于4%的近緣物種還有芝麻(Sesamumindicum)、羅布斯塔咖啡樹(Coffeacanephora)、可可(Theobromacacao)以及蓮(Nelumbonucifera)，其余近3/5的被注釋Unigenes分布于其它511個物種中。由于缺乏長梗杜鵑基因組和轉(zhuǎn)錄組信息，尚存一定數(shù)量的Unigenes在NR庫中仍未能獲得注釋。

圖3 長梗杜鵑轉(zhuǎn)錄組Unigene的NR注釋物種分布圖Fig.3 NR annotated species distribution of Unigenes of transcriptome for R.longipedicellatum

2.4 長梗杜鵑轉(zhuǎn)錄組Unigene的COG注釋及其分類

將獲得的長梗杜鵑Unigenes比對到COG蛋白質(zhì)直系同源數(shù)據(jù)庫(E≤1e-5)，根據(jù)注釋結(jié)果預(yù)測Unigenes可能的功能，并對其進行分類統(tǒng)計(圖4)。

結(jié)果表明，有16 099條Unigenes(占總數(shù)的21.73%)被注釋到25個COG的功能分類中，其各個分類的基因表達豐度各不相同；共包含26 605個功能注釋信息，基本涵蓋了長梗杜鵑大部分的生命活動。其中，一般功能基因(General function prediction only)為最大類別，占COG庫總功能注釋信息的17.08%；其次是轉(zhuǎn)錄(Transcription)功能，復(fù)制、重組和修復(fù)(Replication,recombination and repair)，翻譯后修飾、蛋白質(zhì)折疊和分子伴侶(Posttranslational modification,protein turnover and chaperones)以及信號傳遞機制(Signal transduction mechanisms)，分別占總功能注釋信息的10.22%、8.75%、7.39%和6.76%。而核結(jié)構(gòu)類基因(Nuclear structure)和胞外結(jié)構(gòu)類基因(Extracellular structures)較少，僅有4(0.02%)和9(0.03%)個。另外有1 251個功能注釋信息未明確其生物學(xué)功能，占4.7%。

圖4 長梗杜鵑轉(zhuǎn)錄組Unigene的COG功能注釋分布統(tǒng)計圖Fig.4 COG functional annotation distribution of Unigenes of transcriptome for R.longipedicellatum

圖5 長梗杜鵑轉(zhuǎn)錄組Unigene的GO功能分類統(tǒng)計圖Fig.5 GO functional classification of Unigenes of transcriptome for R.longipedicellatum

2.5 長梗杜鵑轉(zhuǎn)錄組Unigene的GO注釋及其分類

根據(jù)NR注釋結(jié)果進行GO基因功能分類(參數(shù)為Blast2GO軟件默認值)，綜合描述長梗杜鵑葉片中相關(guān)基因的生物學(xué)特征，從宏觀上了解其表達基因的功能分布情況，對獲得相應(yīng)GO條目的Unigenes進行統(tǒng)計分析(圖5)?？芍?7 518(23.64%)條Unigenes得到92 861個GO注釋，平均每條5.3個；近3/4 Unigenes在GO庫中無同源基因與之匹配，進一步說明當(dāng)前長梗杜鵑含有的基因信息較為缺乏。

獲得的GO注釋可劃分為生物過程、細胞組分和分子功能三大類，其中執(zhí)行生物學(xué)過程功能的注釋最多(占41.53%)，涉及分子功能的最少(占20.73%)。3個大的類別又被劃分為更加詳細的56個小類別。依據(jù)其參與的生物學(xué)過程，38 567個GO條目被劃分為24個亞類，其中，以細胞過程(cellular process)最多，占該亞類的24.22%；其次是代謝過程(metabolic process，23.90%)和單一有機體過程(single-organism process，16.65%)；而以細胞殺傷(cell killing)最低，僅占生物過程全部的0.002 6%。根據(jù)其在細胞中所處的位置，35 048個GO條目被劃分為18個亞類，其中，以細胞(cell)、細胞組分(cell part)以及膜(membrane)較多，Unigenes分別占該類型的20.06%、19.90%和17.38%，而以細胞外基質(zhì)組分(extracellular matrix component)和細胞外基質(zhì)(extracellular matrix)較少，僅占0.003%和0.014%。分子功能包含的14個亞類中，以催化活性(catalytic activity)和結(jié)合(binding)居多，分別占該類別的45.71%和40.86%，而以蛋白標(biāo)簽(protein tag)和金屬伴侶蛋白活性(metallochaperone activity)較少，分別占0.016%和0.021%。從這些結(jié)果可以看出，長梗杜鵑在細胞活動、代謝活動的基因表達豐度較高，表明長梗杜鵑自身具有較強的代謝能力。

2.6 長梗杜鵑轉(zhuǎn)錄組Unigene的KEGG代謝通路分析

將長梗杜鵑的Unigenes序列映射到KEGG代謝庫(E≤1e-5)，根據(jù)注釋信息對其可能參與或涉及的代謝途徑進行統(tǒng)計分析。結(jié)果表明，30401條Unigenes獲得注釋(41.03%)；其參與的代謝通路可歸為6大類別、21個亞類(圖6)。

圖6 長梗杜鵑轉(zhuǎn)錄組Unigene的KEGG功能注釋分布統(tǒng)計圖Fig.6 KEGG functional annotation distribution of Unigenes of transcriptome for R.longipedicellatum

由圖6可知，6個代謝通路大類中，代謝(Metabolism)和遺傳信息處理(Genetic Information Processing)相關(guān)的通路所占比例最高，分別為60.79%和23.25%，而以人類疾病(Human Diseases)相關(guān)的通路最少，僅占0.57%，環(huán)境信息處理(Environmental Information Processing)、生物系統(tǒng)(Organismal Systems)以及細胞過程(Cellular Processes)相關(guān)的通路分別占4.84%、5.04%和5.52%。進一步細分為21類代謝途徑，其代謝包含的11個亞類中，全局整體映射(Global and overview maps)所占比例最高，為38.36%，其次是碳水化合物代謝(Carbohydrate metabolism，14.65%)和氨基酸代謝(Amino acid metabolism，8.12%)，而萜類化合物和聚酮化合物代謝(Metabolism of terpenoids and polyketides)所占比例最低，為4.10%。細胞過程和生物系統(tǒng)僅分別包括運輸和代謝(Transport and catabolism)以及環(huán)境適應(yīng)(Environmental adaptation)1個亞類，其余遺傳信息處理、人類疾病和環(huán)境信息處理相關(guān)的通路分別包括4、2、2個亞類，這些亞類中，翻譯，折疊、分類和降解(Folding,sorting and degradation)，所占比例較高，分別占總注釋量的9.38%和6.78%；而內(nèi)分泌及代謝疾病(Endocrine and metabolic diseases，0.537%)以及抗藥性(Antimicrobial resistance，0.029%)較少?？偟腒EGG注釋結(jié)果表明，該時期長梗杜鵑代謝活動和遺傳信息處理能力較強。

2.7 長梗杜鵑轉(zhuǎn)錄組Unigene的CDS預(yù)測

對編碼序列(Coding Sequence,CDS)進行預(yù)測分析能為長梗杜鵑的基因功能發(fā)掘利用和遺傳連鎖圖譜構(gòu)建提供重要參考。本研究按NR、Swissprot、KEGG和COG注釋庫的順序共檢測出39 418個Unigenes的最佳比對片段作為CDS，其總長度為35 487 000 nt，平均長度、N50以及GC含量分別為900 nt、1 350 nt和46.10%；未注釋上的Unigenes使用ESTScan預(yù)測后獲得3 194個CDS，其總長度為1 046 334 nt，平均長度、N50以及GC含量分別為327 nt、321 nt和49.31%，其長度分布如圖7所示。從圖7A中可看出，根據(jù)注釋結(jié)果檢測出的CDS序列較長，1 000 nt以上的占33.34%，其中1 000～2 000 nt的占總數(shù)24.61%，2 000～3 000 nt的占6.13%，3 000 nt以上的僅占2.6%。而通過ESTScan預(yù)測后獲得的CDS序列較短，大部分集中在200～800 nt(占96.87%)，且大于1 600 nt的序列僅6條(圖7：B)。

圖7 長梗杜鵑轉(zhuǎn)錄組Unigene的CDS長度分布圖Fig.7 CDS length distribution of Unigenes of transcriptome for R.longipedicellatum

2.8 長梗杜鵑轉(zhuǎn)錄組Unigene的TF編碼能力及SNP檢測

本研究共預(yù)測出60個基因家族共1 488個編碼轉(zhuǎn)錄因子(Transcription Factor,TF)的Unigenes(圖8)；其中MYB、MYB-related、AP2-EREBP、bHLH以及mTERF類轉(zhuǎn)錄因子為多基因家族，分別占總預(yù)測量的16%、12.03%、7.19%、6.72%和6.05%，而ULT、LFY、DBP所占比例較少，均為0.07%。同時，還預(yù)測出80條與非生物脅迫有關(guān)的TF，包括NAC(67條)和HSF(13條)；72條與生物脅迫的抗性密切相關(guān)的TF，包括WRKY轉(zhuǎn)錄因子超家族(58條)和bZIP類轉(zhuǎn)錄因子(14條)以及29條可能參與調(diào)節(jié)花發(fā)育的MADS框基因。這些TF的發(fā)現(xiàn)表明長梗杜鵑轉(zhuǎn)錄調(diào)控的復(fù)雜性，也為后續(xù)長梗杜鵑逆境相關(guān)基因的表達調(diào)控以及性狀的利用和改良等研究提供了新的思路。

此外，還檢測出57 927個單核苷酸多態(tài)性(Single Nucleotide Polymorphsims,SNP)位點(圖9)，且在Unigenes序列上呈現(xiàn)不均勻分布，轉(zhuǎn)換類型高于顛換。其中，轉(zhuǎn)換類型(A-G和C-T)有36 171個，占總變異數(shù)的62.44%；顛換類型(A-C、A-T、C-G和G-T)有21 756個，占37.56%。在所有類型中，C-T變異的SNP數(shù)目最多(50.02%)，與杜瑋南等[15]所得結(jié)論一致，SNP多發(fā)生在C和T之間。長梗杜鵑SNP分析，結(jié)合Unigenes的表達量及表型，可在mRNA層面的基因型差異區(qū)別于其它物種，對了解其親屬關(guān)系和進化關(guān)系具有重要的生物學(xué)意義，也為長梗杜鵑SNP標(biāo)記的后續(xù)分子實驗驗證及其在自然種群中的應(yīng)用奠定了基礎(chǔ)。

圖8 長梗杜鵑轉(zhuǎn)錄組Unigene的轉(zhuǎn)錄因子家族分類圖Fig.8 Transcription Factor family classification of Unigenes of transcriptome for R.longipedicellatum

圖9 長梗杜鵑轉(zhuǎn)錄組Unigene的SNP變異類型分布圖Fig.9 SNP variant type distribution of Unigenes of transcriptome for R.longipedicellatum

3 討論

隨著高通量測序技術(shù)的快速發(fā)展以及測序平臺的升級改進，通過轉(zhuǎn)錄組測序獲取大量EST、發(fā)掘功能基因成為一種非常高效、可靠且高輸出的重要手段，既降低了成本，又保證了測序的通量和所得序列的長度與質(zhì)量，適合缺乏基因組信息的非模式物種開展生物信息學(xué)相關(guān)研究。雖然杜鵑屬植物種類超過1 000種，但基因和基因組學(xué)研究相對滯后，遺傳信息較為缺乏。與云錦杜鵑(R.fortunei)根系(4.66 Gb)轉(zhuǎn)錄組相比[16]，本研究獲得的信息量較飽和，共產(chǎn)生6.73 Gb的數(shù)據(jù)量，包含74 092條Unigenes，為后續(xù)Unigenes的功能注釋和分類、CDS預(yù)測、有利基因克隆及功能研究等方面奠定了堅實的基礎(chǔ)。De novo組裝后獲得的Unigenes平均長度、N50值均較大，分別為938 nt和1 616 nt，遠高于其它杜鵑屬植物如酒紅杜鵑(R.catawbiense)[17]、沼澤杜鵑(R.viscosum)[18～19]、云錦杜鵑[16]等，說明本研究所得序列中長片段較多，組裝效果較好；堿基Q20和Q30分別為98.22%、95.20%。當(dāng)Q30值在80%以上就認為測序質(zhì)量非?？煽縖20]，而堿基Q20與堿基識別的錯誤率呈對數(shù)相關(guān)，其表示每100個序列堿基中僅有一個出錯的概率[21]。長梗杜鵑轉(zhuǎn)錄組Unigenes在1 kB以上的有23 879條，且clusters有51 505條，singletons有22 587條，所得clusters數(shù)目為單獨Unigenes的2倍多，進一步說明所得序列的高效性。

將長梗杜鵑Unigenes序列比對到NR、GO、COG、KEGG等7大功能數(shù)據(jù)庫，共有45538條獲得注釋，占總體的61.46%，且所得Unigenes的表達量FPKM值集中在0.2～28.98，表明本研究檢測到長梗杜鵑低表達水平基因的比例很高。在NR數(shù)據(jù)庫成功注釋的Unigenes中，25.73%的Unigenes注釋為與葡萄同源的基因，比例遠遠高于其它被注釋的515個物種,其可能是由于葡萄具有參考基因組且與長梗杜鵑在生活史和進化關(guān)系上較為接近的緣故。而與其比對上的17個同屬物種(包括栽培品種)僅120條Unigenes，其中映山紅(R.simsii)、臺北杜鵑(R.kanehirai)、錦繡杜鵑(R.pulchrum)、酒紅杜鵑以及臺灣杜鵑(R.formosanum)相對較多，分別有54(0.14%)、14(0.04%)、13(0.03%)、10(0.03%)和7(0.02%)條，其余12個種與之比對上的基因數(shù)所占比例均≤0.01%，這可能是現(xiàn)有杜鵑屬植物基因組信息較為缺乏或是所得Unigenes多為長梗杜鵑特有的新功能基因而未能獲得相應(yīng)匹配?？偟奈幢蛔⑨孶nigenes有28 554條，占總體的38.54%，這些片段可能是因為較短或是本身為非編碼序列而未能在7大功能數(shù)據(jù)庫中獲得同源性比對[22]。在28 554個無匹配序列中，長度在300～1 000 nt的Unigenes有27 398個，占96%，且隨著序列長度增加，無同源序列的Unigenes比例明顯降低，與賈新平[20]所得結(jié)論一致。轉(zhuǎn)錄組中Unigenes的序列越長，獲得注釋信息的可能性就越高。

利用COG數(shù)據(jù)庫對所有Unigenes序列進行比對，獲得了大量長梗杜鵑某一時期生長發(fā)育過程中表達的基因，16 099條序列共26 605個功能注釋信息被劃分為25類，一般功能、轉(zhuǎn)錄和翻譯的基因表達豐度較高。通過GO數(shù)據(jù)庫比對，根據(jù)其參與的生物學(xué)過程、構(gòu)成的細胞組分和實現(xiàn)的分子功能分為56個小類別，以細胞活動、代謝活動和催化活性較為豐富，表明長梗杜鵑具有較強的代謝能力。根據(jù)KEGG數(shù)據(jù)庫分析其可能參與或涉及的代謝通路，共30 401個Unigenes注釋到6個代謝通路大類、21類代謝途徑，其中定位到代謝相關(guān)通路的基因數(shù)最多，約占總注釋量的3/5，進一步證實長梗杜鵑這一時期具有較強的代謝活動。更為值得一提的是，KEGG的代謝通路分析發(fā)掘出176條與人類疾病相關(guān)的Unigenes，包括內(nèi)分泌及代謝疾病(167條)和抗藥性(9條)。長梗杜鵑KEGG代謝途徑的分析為進一步探索其抗性機理、挖掘花期調(diào)控相關(guān)功能基因，加快在引種馴化和雜交育種等方面的研究開發(fā)提供了分子依據(jù)。

轉(zhuǎn)錄因子在轉(zhuǎn)錄水平上與DNA發(fā)生相互作用的多樣性決定了其生物學(xué)功能的多樣性，本研究根據(jù)組裝結(jié)果預(yù)測出60個基因家族共1 488個編碼轉(zhuǎn)錄因子的Unigenes，這些家族在長梗杜鵑生長發(fā)育過程中發(fā)揮著不同的作用。同擬南芥、棉花、草莓和番茄等大多數(shù)植物一樣[23]，MYB基因家族在長梗杜鵑中普遍存在，是其中數(shù)量最大、功能最多樣的轉(zhuǎn)錄家族之一，可能與其代謝調(diào)控、激素的信號傳導(dǎo)、細胞周期調(diào)控有關(guān)。同時，107個植物特有轉(zhuǎn)錄因子AP2-EREBP和100個植物第二大類轉(zhuǎn)錄因子bHLH的發(fā)現(xiàn)，也為長梗杜鵑生長發(fā)育和應(yīng)答外界環(huán)境信號過程、進化和基因間的相互作用以及抗逆調(diào)控網(wǎng)絡(luò)的研究奠定了基礎(chǔ)。EST文庫是篩選SNP的重要來源，本研究通過轉(zhuǎn)錄組測序分析，共檢測出57 927個SNP多態(tài)位點，包括36 171個轉(zhuǎn)換類型，21 756個顛換類型，將為今后開展長梗杜鵑遺傳圖譜構(gòu)建、重要性狀基因定位、遺傳多樣性分析、群體選擇分析、品種鑒定以及分子標(biāo)記輔助選擇育種等研究提供了重要參考。

4 結(jié)論

本研究采用Illumina HiSeq 4000最新高通量測序平臺建立了新物種長梗杜鵑轉(zhuǎn)錄組數(shù)據(jù)庫，獲得了74 092條與該物種生長發(fā)育相關(guān)的Unigenes，通過生物信息學(xué)方法對其進行功能注釋、代謝相關(guān)通路、CDS預(yù)測、TF編碼能力預(yù)測、SNP檢測以及表達量計算等分析。綜合COG、GO和KEGG三大功能數(shù)據(jù)庫的注釋結(jié)果，表明長梗杜鵑在該時期轉(zhuǎn)錄、復(fù)制和翻譯等基因表達豐度較高，具有較強的代謝活動和遺傳信息處理能力。所得轉(zhuǎn)錄組信息作為杜鵑屬植物基因組序列的重要補充，進一步豐富了該屬植物的基因表達信息平臺，為長梗杜鵑乃至同屬植物基因克隆及功能分析、SSR、SNP分子標(biāo)記的規(guī)模化開發(fā)奠定了數(shù)據(jù)基礎(chǔ)。同時，隨著新一代高通量測序技術(shù)的廣泛應(yīng)用，將長梗杜鵑Unigenes與其它同屬植物或近緣物種進行比較分析，可為進一步研究其進化和有利基因的篩查提供更為重要的信息。

1.方瑞征,閔天祿.杜鵑屬植物區(qū)系的研究[J].云南植物研究,1995,17(4):359-379.

Fang R Z,Min T L.The floristic study on the genusRhododendron[J].Acta Botanica Yunnanica,1995,17(4):359-379.

2.中國科學(xué)院中國植物志編輯委員會.中國植物志[M].北京:科學(xué)出版社,1999:1-3.

Delectis Florae Reipublicae Popularis Sinicae,Agendae Academiae Sinicae Edita.Flora reipublicae popularis sinicae[M].Beijing:Science Press,1999:1-3.

3.張長芹,高連明,薛潤光,等.中國杜鵑花的保育現(xiàn)狀和展望[J].廣西科學(xué),2004,11(4):354-359,362.

Zhang C Q,Gao L M,Xue R G,et al.A general review of the research and conservation statue of ChineseRhododendron[J].Guangxi Sciences,2004,11(4):354-359,362.

4.張長芹,羅吉鳳.杜鵑花新品種‘金躑躅’和‘紫艷’[J].園藝學(xué)報,2002,29(5):502.

Zhang C Q,Luo J F.NewRhododendronvarieties-‘Jinzhizhu’ and ‘Ziyan’[J].Acta Horticulturae Sinica,2002,29(5):502.

5.張長芹,羅吉鳳,馮寶均.杜鵑花新品種‘朝暉’和‘紅暈’[J].園藝學(xué)報,2002,29(3):296.

Zhang C Q,Luo J F,Feng B J.NewRhododendronhybrid-‘Zhaohui’ and ‘Hongyun’[J].Acta Horticulturae Sinica,2002,29(3):296.

6.蘭熙,張樂華,張金政,等.杜鵑花屬植物育種研究進展[J].園藝學(xué)報,2012,39(9):1829-1838.

Lan X,Zhang L H,Zhang J Z,et al.Research progress ofRhododendronbreeding[J].Acta Horticulturae Sinica,2012,39(9):1829-1838.

7.Altschul S F,Gish W,Miller W,et al.Basic local alignment search tool[J].Journal of Molecular Biology,1990,215(3):403-410.

8.Conesa A,G?tz S,García-Gómez J M,et al.Blast2GO:a universal tool for annotation,visualization and analysis in functional genomics research[J].Bioinformatics,2005,21(18):3674-3676.

9.Quevillon E,Silventoinen V,Pillai S,et al.InterProScan:protein domains identifier[J].Nucleic Acids Research,2005,33(S2):W116-W120.

10.Iseli C,Jongeneel C V,Bucher P.ESTScan:a program for detecting,evaluating,and reconstructing potential coding regions in EST sequences[C].//Proceedings of the Seventh International Conference on Intelligent Systems for Molecular Biology.Menlo Park:AAAI Press,1999:138-148.

11.Rice P,Longden I,Bleasby A.EMBOSS:the European molecular biology open software suite[J].Trends in Genetics,2000,16(6):276-277.

12.Mistry J,Finn R D,Eddy S R,et al.Challenges in homology search:HMMER3 and convergent evolution of coiled-coil regions[J].Nucleic Acids Research,2013,41(12):e121.

13.Kim D,Langmead B,Salzberg S L.HISAT:a fast spliced aligner with low memory requirements[J].Nature Methods,2015,12(4):357-360.

14.Mckenna A,Hanna M,Banks E,et al.The genome analysis toolkit:a MapReduce framework for analyzing next-generation DNA sequencing data[J].Genome Research,2010,20(9):1297-1303.

15.杜瑋南,孫紅霞,方福德.單核苷酸多態(tài)性的研究進展[J].中國醫(yī)學(xué)科學(xué)院學(xué)報,2000,22(4):392-394.

Du W N,Sun H X,Fang F D.The research development of single nucleotide polymorphism[J].Acta Academiae Medicinae Sinicae,2000,22(4):392-394.

16.Wei X Y,Chen J J,Zhang C Y,et al.Differential gene expression inRhododendronfortuneiroots colonized by an ericoid mycorrhizal fungus and increased nitrogen absorption and plant growth[J].Frontiers in Plant Science,2016,7:1594.

17.Wei H,Dhanaraj A L,Rowland L J,et al.Comparative analysis of expressed sequence tags from cold-acclimated and non-acclimated leaves ofRhododendroncatawbienseMichx[J].Planta,2005,221(3):406-416.

18.Susko A Q,Rinehart T,Bradeen J,et al.EST-SSR development and validation in deciduous azalea breeding populations from transcriptome sequencing[C].//Plant and Animal Genomic XXIII Conference.San Diego,CA:University of Minnesota,2015.

19.Susko A Q.Phenotypic and genetic variation for rhizosphere acidification,a candidate trait for pH adaptability,in deciduous azalea(Rhododendronsect.Pentanthera)[D].Minnesota:University of Minnesota,2016.

20.賈新平,孫曉波,鄧衍明,等.鳥巢蕨轉(zhuǎn)錄組高通量測序及分析[J].園藝學(xué)報,2014,41(11):2329-2341.

Jia X P,Sun X B,Deng Y M,et al.Sequencing and analysis of the transcriptome ofAspleniumnidus[J].Acta Horticulturae Sinica,2014,41(11):2329-2341.

21.Ewing B,Hillier L,Wendl M C,et al.Base-calling of automated sequencer traces usingPhred.Ⅰ.accuracy assessment[J].Genome Research,1998,8(3):175-185.

22.蔡年輝,鄧麗麗,許玉蘭,等.基于高通量測序的云南松轉(zhuǎn)錄組分析[J].植物研究,2016,36(1):75-83.

Cai N H,Deng L L,Xu Y L,et al.Transcriptome analysis forPinusyunnanensisbased on high throughput sequencing[J].Bulletin of Botanical Research,2016,36(1):75-83.

23.牛義嶺,姜秀明,許向陽.植物轉(zhuǎn)錄因子MYB基因家族的研究進展[J].分子植物育種,2016,14(8):2050-2059.

Niu Y L,Jiang X M,Xu X Y.Research advances on transcription factor MYB gene family in plant[J].Molecular Plant Breeding,2016,14(8):2050-2059.

The Cultivation Project of Technology Innovation Talent of Yunnan Province

introduction：LI Tai-Qiang(1993—),male,graduate student,mainly engaged in the study of conservation biology inRhododendron.

date:2017-03-29

TranscriptomeAnalysisforRhododendronlongipedicellatum(PlantSpecieswithExtremelySmallPopulations)BasedonHighThroughputSequencing

LI Tai-Qiang LIU Xiong-Fang WAN You-Ming LI Zheng-Hong QI Guo-HaiLI Yu-Ying LIU Xiu-Xian HE Rui MA Yan MA Hong*

(The Research Institute of Resource Insects,Chinese Academy of Forestry,Kunming 650233)

To strengthen the research of resources evaluation, protection and identification of endemic and endangered species ofRhododendronlongipedicellatum, and to provide a helpful reference for genetic breeding and improvement of its agronomic traits, the transcriptome was sequenced by using Illumina Hiseq 4 000, in total, 74 092 Unigenes with an average length of 938 nt, N50of 1616 nt, Q20of 98.22%, Q30of 95.20% and GC content of 43.24% were obtained by de novo assembly and cluster with filtered data, and there were 23 879 Unigenes with more than 1 kB. Then, the Unigenes were annotated by 7 functional databases, and finally, 39 876(NR: 53.82%)，38 065(NT: 51.38%), 27 384(Swissprot: 36.96%), 16 099(COG: 21.73%), 30 401(KEGG: 41.03%), 17 518(GO: 23.64%), and 29 676(Interpro: 40.05%) Unigenes were annotated. The Unigenes were roughly divided into three functional categories(i.e. biological processes, cellular components and molecular function) and 56 sub-categories according to GO function. Most of the genes performed biological processes. KEGG functional annotation analysis showed that Unigenes could be grouped into 6 categories, 32 metabolic pathways. 176 Unigenes relating to human diseases were detected, including endocrine and metabolic diseases(167) and antimicrobial resistance(9). The 39 418 CDS were detected with functional annotation results, and after 3 194 CDS also were predicted by ESTScan with the remaining Unigenes. The 1 488 Transcription Factor(TF) coding Unigenes were predicted and 57 927 SNP polymorphic loci were detected. The analysis of the transcriptome could lay a foundation for further study of functional gene discovery and utilization, resistance mechanism analysis, classification and evolution of genetic resources, molecular marker development and molecular assisted breeding ofR.longipedicellatumand other congeneric species.

Rhododendronlongipedicellatum;transcriptome;Unigene;gene annotation;CDS;Transcription Factor;SNP

“云南省技術(shù)創(chuàng)新人才”培養(yǎng)對象項目(2016HB007)

李太強(1993—)，男，碩士研究生，主要從事杜鵑屬植物保護生物學(xué)研究。

* 通信作者：E-mail:hortscience@163.com

2017-03-29

* Corresponding author：E-mail:hortscience@163.com

S685.21

A

10.7525/j.issn.1673-5102.2017.06.004