李太強(qiáng) 劉雄芳 萬(wàn)友名 李正紅 起國(guó)海 李鈺瑩 劉秀賢 和 銳 馬 艷 馬 宏

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

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

李太強(qiáng) 劉雄芳 萬(wàn)友名 李正紅 起國(guó)海 李鈺瑩 劉秀賢 和 銳 馬 艷 馬 宏*

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

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

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

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

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

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

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

1 材料與方法

1.1 試驗(yàn)材料

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

1.2 長(zhǎng)梗杜鵑葉片總RNA提取

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

1.3 轉(zhuǎn)錄組測(cè)序

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

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

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

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

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

2 結(jié)果與分析

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

測(cè)序共獲得58.30 Mb的Raw Reads，過濾后得到44.85 Mb的Clean Reads，占原始讀長(zhǎng)的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測(cè)序?qū)е虑? bp出現(xiàn)正常波動(dòng)外，其余reads每個(gè)位置的堿基含量分布穩(wěn)定，無AT或GC分離現(xiàn)象。其次，在Clean reads的堿基質(zhì)量分布圖中(圖1:B)，顏色深的點(diǎn)質(zhì)量值集中于30～40，說明達(dá)到該質(zhì)量值范圍的堿基數(shù)目較多，而低質(zhì)量(質(zhì)量值小于20)的堿基比例較低，表明測(cè)序質(zhì)量較好，測(cè)序所得reads準(zhǔn)確性較高。

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

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

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

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

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

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

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

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

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

結(jié)果表明，有16 099條Unigenes(占總數(shù)的21.73%)被注釋到25個(gè)COG的功能分類中，其各個(gè)分類的基因表達(dá)豐度各不相同；共包含26 605個(gè)功能注釋信息，基本涵蓋了長(zhǎng)梗杜鵑大部分的生命活動(dòng)。其中，一般功能基因(General function prediction only)為最大類別，占COG庫(kù)總功能注釋信息的17.08%；其次是轉(zhuǎn)錄(Transcription)功能，復(fù)制、重組和修復(fù)(Replication,recombination and repair)，翻譯后修飾、蛋白質(zhì)折疊和分子伴侶(Posttranslational modification,protein turnover and chaperones)以及信號(hào)傳遞機(jī)制(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%)個(gè)。另外有1 251個(gè)功能注釋信息未明確其生物學(xué)功能，占4.7%。

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

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

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

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

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

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

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

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

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

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

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

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

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

本研究共預(yù)測(cè)出60個(gè)基因家族共1 488個(gè)編碼轉(zhuǎn)錄因子(Transcription Factor,TF)的Unigenes(圖8)；其中MYB、MYB-related、AP2-EREBP、bHLH以及mTERF類轉(zhuǎn)錄因子為多基因家族，分別占總預(yù)測(cè)量的16%、12.03%、7.19%、6.72%和6.05%，而ULT、LFY、DBP所占比例較少，均為0.07%。同時(shí)，還預(yù)測(cè)出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)表明長(zhǎng)梗杜鵑轉(zhuǎn)錄調(diào)控的復(fù)雜性，也為后續(xù)長(zhǎng)梗杜鵑逆境相關(guān)基因的表達(dá)調(diào)控以及性狀的利用和改良等研究提供了新的思路。

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

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

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

3 討論

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

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

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

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

4 結(jié)論

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

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.中國(guó)科學(xué)院中國(guó)植物志編輯委員會(huì).中國(guó)植物志[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.張長(zhǎng)芹,高連明,薛潤(rùn)光,等.中國(guó)杜鵑花的保育現(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.張長(zhǎng)芹,羅吉鳳.杜鵑花新品種‘金躑躅’和‘紫艷’[J].園藝學(xué)報(bào),2002,29(5):502.

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

5.張長(zhǎng)芹,羅吉鳳,馮寶均.杜鵑花新品種‘朝暉’和‘紅暈’[J].園藝學(xué)報(bào),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ìn)展[J].園藝學(xué)報(bào),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ìn)展[J].中國(guó)醫(yī)學(xué)科學(xué)院學(xué)報(bào),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)錄組高通量測(cè)序及分析[J].園藝學(xué)報(bào),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.蔡年輝,鄧麗麗,許玉蘭,等.基于高通量測(cè)序的云南松轉(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.牛義嶺,姜秀明,許向陽(yáng).植物轉(zhuǎn)錄因子MYB基因家族的研究進(jìn)展[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)對(duì)象項(xiàng)目(2016HB007)

李太強(qiáng)(1993—)，男，碩士研究生，主要從事杜鵑屬植物保護(hù)生物學(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