許薔薇，樓雄珍，楊 彬，林二培，童再康

（1.浙江農(nóng)林大學(xué) 林業(yè)與生物技術(shù)學(xué)院，浙江 杭州311300；2.浙江農(nóng)林大學(xué) 省部共建亞熱帶森林培育國(guó)家重點(diǎn)實(shí)驗(yàn)室，浙江 杭州311300）

隨著第2代測(cè)序技術(shù)的發(fā)展，轉(zhuǎn)錄組測(cè)序（RNA-Seq）已成為植物分子生物學(xué)研究的重要手段，廣泛應(yīng)用于功能基因挖掘、分子標(biāo)記開發(fā)、代謝通路和調(diào)控機(jī)制研究等方面［1-3］。該技術(shù)具有成本低、數(shù)據(jù)量大、效率高、準(zhǔn)確性高等優(yōu)點(diǎn)［4］，可對(duì)組織或者細(xì)胞中所有RNA進(jìn)行測(cè)序，并通過讀段（reads）的拼接和豐度統(tǒng)計(jì)獲得相應(yīng)的轉(zhuǎn)錄本序列信息及其表達(dá)水平［5］。云錦杜鵑Rhododendron fortunei為杜鵑花科Ericaceae杜鵑花屬Rhododendron植物，為中國(guó)特有種，主要分布于安徽、湖南、湖北、浙江等省［6］。云錦杜鵑葉形大，花形美麗具清香，具有較高的園藝觀賞價(jià)值；適應(yīng)性較強(qiáng)，易人工栽培，也常作為杜鵑屬種間雜交的親本；此外，枝葉等組織中含有槲皮素、山柰酚、楊梅素等有效活性成分，可開發(fā)入藥［7-8］。可見，云錦杜鵑作為一種兼具觀賞價(jià)值和藥用價(jià)值的優(yōu)良木本花卉，開發(fā)潛力巨大，具有良好的產(chǎn)業(yè)化前景。目前，云錦杜鵑的研究主要集中在野生資源調(diào)查［9］、群落結(jié)構(gòu)［10］、無性繁殖技術(shù)［11］、有效成分分析［8,12］、光合特性［13］、菌根共生等方面［14］。由于缺乏遺傳信息，云錦杜鵑的遺傳多樣性、分子輔助育種等的研究較為滯后。本研究通過RNA-seq高通量測(cè)序技術(shù)對(duì)云錦杜鵑轉(zhuǎn)錄組進(jìn)行測(cè)定，通過序列拼接、功能注釋和分析，獲取大量的序列信息，旨在為云錦杜鵑分子標(biāo)記開發(fā)，遺傳多樣性分析，以及功能基因挖掘、重要性狀形成分子機(jī)制等研究提供序列信息。

1 材料與方法

1.1 植物材料

前期以天臺(tái)華頂林場(chǎng)云錦杜鵑優(yōu)株莖段為外植體，通過組培獲得云錦杜鵑無性系ZL-01。以該無性系組培苗為材料，取其嫩葉、嫩莖和根，液氮冷凍后保存于-80℃用于RNA提取。

1.2 RNA提取與質(zhì)控

采用Trizol試劑（Invitrogen），按照試劑說明書分別提取組培苗根、莖、葉的總RNA。采用Agilent 2100 Bioanalyzer（Agilent），對(duì)總RNA樣品的純度、濃度和完整性進(jìn)行檢測(cè)評(píng)估，等量混合后用于測(cè)序文庫(kù)構(gòu)建。

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

取檢驗(yàn)合格的RNA，用帶有Oligo（dT）的磁珠富集mRNA，將得到的mRNA隨機(jī)打斷成短片段，用六堿基隨機(jī)引物合成一鏈cDNA，然后加入緩沖液、 dNTPs、DNA polymerase I和RNase H合成二鏈cDNA。雙鏈cDNA經(jīng)純化后進(jìn)行末端修復(fù)，加ploy（A）并連接測(cè)序接頭，通過PCR擴(kuò)增及純化后得到測(cè)序文庫(kù)。采用Illumina Hiseq 2500進(jìn)行轉(zhuǎn)錄組的測(cè)序，測(cè)序由北京諾禾致源生物公司完成。

1.4 轉(zhuǎn)錄組數(shù)據(jù)分析

針對(duì)測(cè)序得到的原始數(shù)據(jù)（raw data）進(jìn)行接頭和低質(zhì)量測(cè)序片段（reads）去除等處理，獲得高質(zhì)量的干凈數(shù)據(jù)（clean data）；計(jì)算干凈數(shù)據(jù)的Q20和Q30，GC含量和堿基錯(cuò)誤率，評(píng)估測(cè)序質(zhì)量；利用序列的重疊，通過Trinity軟件對(duì)干凈數(shù)據(jù)數(shù)據(jù)進(jìn)行序列組裝，將短測(cè)序片段延伸成較長(zhǎng)的片段，并得到片段集合， 利用 De Bruijin方法［15］得到轉(zhuǎn)錄本（transcripts）和單基因簇（unigene）序列。

1.5 基因功能注釋

利用BLAST和HMMER軟件將云錦杜鵑單基因簇序列與公共數(shù)據(jù)庫(kù)進(jìn)行比對(duì)，根據(jù)基因的相似性進(jìn)行功能注釋，得到與給定單基因簇具有最高序列相似性的蛋白，從而得到該單基因簇的蛋白功能注釋信息，其中BLAST參數(shù)E≤1e-5，HMMER參數(shù)E≤1e-10作為篩選標(biāo)準(zhǔn)。比對(duì)采用的公共數(shù)據(jù)庫(kù)包括美國(guó)生物信息中心（NCBI）非冗余蛋白數(shù)據(jù)庫(kù)（Non-redundant protein database，Nr），NCBI核酸序列數(shù)據(jù)庫(kù)（Non-redundant nucleotide sequences， Nt）， 蛋白質(zhì)序列數(shù)據(jù)庫(kù)（Swiss Prot protein database， SwissProt），真核直系同源基因數(shù)據(jù)庫(kù)（Eukaryotic orthologous groups，KOG），蛋白質(zhì)家族數(shù)據(jù)庫(kù)（Protein families database，Pfam），基因本體論數(shù)據(jù)庫(kù)（Gene ontology，GO）及京都基因與基金組百科全書（Kyoto encyclopedia of genes and genomes， KEGG）。

1.6 SSR分析

采用MISA軟件對(duì)云錦杜鵑單基因簇進(jìn)行簡(jiǎn)單序列重復(fù)（SSR）分析，鑒定其中6種類型的SSR；各SSR類型重復(fù)次數(shù)設(shè)定為：?jiǎn)魏塑账酳SR，重復(fù)數(shù)≥10次；二核苷酸SSR，重復(fù)數(shù)≥6次；三至六核苷酸SSR，重復(fù)數(shù)≥5次。

2 結(jié)果與分析

2.1 測(cè)序結(jié)果與組裝

分別提取無性系ZL-01組培苗根、莖和葉的RNA。檢測(cè)結(jié)果表明：根RNA的質(zhì)量濃度為362 mg·L-1，28S/18S為1.6，RNA完整性計(jì)數(shù)（RNA integrity number，RIN）為9.6；莖RNA的質(zhì)量濃度為515 mg·L-1，28S/18S為 1.8，RIN值為 9.8；葉 RNA的質(zhì)量濃度為 485 mg·L-1，28S/18S為 1.8，RIN值為9.7。這些指標(biāo)均符合建庫(kù)測(cè)序要求，因此，將這3個(gè)組織的RNA等量混合后進(jìn)行測(cè)序。

測(cè)序數(shù)據(jù)經(jīng)處理后得到94252430個(gè)干凈測(cè)序片段（clean reads），包含了14.14 Gb核苷酸序列信息，GC含量為47.58%。測(cè)序質(zhì)量評(píng)估結(jié)果顯示，堿基錯(cuò)誤率為0.02%，Q20為96.70%，Q30為91.85%。這些表明該轉(zhuǎn)錄組的數(shù)據(jù)量和質(zhì)量均較高，為后續(xù)的序列組裝提供了高質(zhì)量的原始數(shù)據(jù)。

通過Trintiy軟件組裝得到112777個(gè)轉(zhuǎn)錄本，總長(zhǎng)度為91621682 bp，平均長(zhǎng)度為812.4 bp，N50為1410，其中長(zhǎng)度在1 kb以上的有29225條，占25.92%；2 kb以上的10418條，占9.24%（表1）。對(duì)轉(zhuǎn)錄本進(jìn)行聚類和組裝得到84633個(gè)單基因簇，總長(zhǎng)度為58517298 bp，平均長(zhǎng)度為691.4 bp，N50為1177 bp；其中超過1 kb的有16631條，占19.65%；超過2 kb的5760條， 占6.81%（表1）。

表1 云錦杜鵑轉(zhuǎn)錄組組裝結(jié)果Table 1 Summary of transcriptome assembly for R.fortunei

2.2 基因功能注釋

將單基因簇序列與Nr，Nt，SwissProt，GO，KOG，KEGG和Pfam 7個(gè)數(shù)據(jù)庫(kù)進(jìn)行比對(duì)，共有35526條單基因簇獲得成功注釋，占單基因簇的41.97%（表2）。其中，30700條單基因簇獲得Nr數(shù)據(jù)庫(kù)注釋，占36.27%； 23215條單基因簇獲得GO數(shù)據(jù)庫(kù)注釋，占27.43%； 22882條單基因簇獲得Pfam數(shù)據(jù)庫(kù)注釋，占27.03%；22202條單基因簇獲得Swiss-prot數(shù)據(jù)庫(kù)注釋，占26.23%；18046條單基因簇獲得Nt數(shù)據(jù)庫(kù)注釋，占21.32%；11085條單基因簇獲得KOG數(shù)據(jù)庫(kù)注釋，占13.09%；9887條單基因簇獲得KEGG數(shù)據(jù)庫(kù)注釋占11.68%（表2）。

表2 云錦杜鵑單基因簇功能注釋Table 2 Annotation of unigenes in transcriptome of R.fortunei

在Nr庫(kù)中，云錦杜鵑轉(zhuǎn)錄組注釋到其他物種的單基因簇基因序列共30700條。其中與葡萄Vitis vinifera基因序列相似的最多，所占比例為24.39%；其次為中?？Х菴offea canephora，所占比例為7.39%；第3為可可Theobroma cacao，占6.77%；其他相似性序列數(shù)量大于3%的物種有煙草Nicotiana tomentosiformis（6.19%）， 毛果楊Populus trichocarpa（4.09%）， 甜橙Citrus sinensis（3.82%）， 麻瘋樹Jatropha curcas（3.76%）和梅花Prunus mume（3.16%）， 其他物種占 40.43%（圖 1）。

2.3 單基因簇的GO功能分類

通過GO數(shù)據(jù)庫(kù)的注釋，共有23215條單基因簇獲注釋信息，得到119389個(gè)GO功能注釋。由分類結(jié)果可知：生物學(xué)過程（biological process）最多55692條，占46.65%，其次是細(xì)胞組分（cellular component），35577條，占29.80%，最少的分子功能（molecular function），有28120條，占23.55%；這三大功能分類又可分為55個(gè)亞類，其中生物學(xué)過程23個(gè)亞類，細(xì)胞組分18個(gè)亞類，分子功能14個(gè)亞類（圖2）。生物學(xué)過程中，涉及細(xì)胞過程、代謝過程和單一有機(jī)體進(jìn)程的單基因簇較多，分別有12698，12376和9581條；細(xì)胞組分中涉及較多的是細(xì)胞、細(xì)胞部分和大分子復(fù)合體，分別有7216，7214和4713條；分子功能中涉及較多的有結(jié)合功能和催化活性，分別有13302，10524條（圖2）。

圖1 云錦杜鵑單基因簇Nr數(shù)據(jù)庫(kù)比對(duì)相似性物種分布圖Figure 1 Similarity of unigenes of R.fortunei with those of other species in Nr database

圖2 云錦杜鵑單基因簇的GO功能分類Figure 2 GO functional categories of R.fortunei unigenes

2.4 單基因簇的KOG功能注釋

通過KOG數(shù)據(jù)庫(kù)對(duì)云錦杜鵑單基因簇進(jìn)行注釋，結(jié)果顯示有11085條序列獲得12475個(gè)注釋信息，可分為26個(gè)功能分類（圖3）。從基因功能分類來看，涉及一般功能預(yù)測(cè)的序列最多，多達(dá)1970條；涉及翻譯后修飾、蛋白翻轉(zhuǎn)、分子伴侶功能的序列次之，有1527條；而涉及核結(jié)構(gòu)、胞外結(jié)構(gòu)和細(xì)胞運(yùn)動(dòng)的序列很少，僅有37，34和5條（圖3）。

圖3 云錦杜鵑單基因簇的KOG功能分類Figure 3 KOG functional categories of R.fortunei unigenes

2.5 KEGG通路注釋

利用KEGG注釋系統(tǒng)對(duì)云錦杜鵑單基因簇涉及的代謝途徑進(jìn)行分析，結(jié)果顯示：9887條單基因簇得到15455個(gè)注釋，歸屬于272條通路。按獲得注釋的基因數(shù)量進(jìn)行排序，取前20個(gè)途徑，發(fā)現(xiàn)含有200條單基因簇以上的通路有10個(gè)，涉及碳代謝的單基因簇最多，有377條；其次是與核糖體相關(guān)的單基因簇，有364條；第3位是與氨基酸合成相關(guān)的單基因簇，有348條；涉及其他途徑有內(nèi)質(zhì)網(wǎng)蛋白加工（319）、剪切體（256）、植物激素信號(hào)轉(zhuǎn)導(dǎo)（244）、淀粉和蔗糖代謝（240）、RNA轉(zhuǎn)運(yùn)（227）、氧化磷酸化（224）和植物病原物相互作用（210）；其他通路的單基因簇?cái)?shù)量均在200以下（表3）。

表3 云錦杜鵑單基因簇的KEGG分析Table 3 Summary KEGG pathway of R.fortunei transcriptome

2.6 MADS-box基因家族分析

MADS-box基因家族在植物花分生組織形成、花器官發(fā)育等過程中發(fā)揮關(guān)鍵作用［16］。鑒定轉(zhuǎn)錄組中的MADS-box基因有助于進(jìn)一步研究或調(diào)控云錦杜鵑成花過程。通過與Nr，Nt和Swiss-prot三大數(shù)據(jù)庫(kù)比對(duì)，共找出編碼MADS-box基因的單基因簇序列24條，分別屬于10個(gè)不同的亞家族（表4）。其中AGL17亞家族成員最多，包含c41091_g4，c39737_g1，c36841_g1，c38538_g1和c40334_g4；SQUA和SVP亞家族則各有 4個(gè)成員，分別為 c50592_g1，c10093_g1，c9572_g1，c2583_g1與 c37773_g1，c31351_g1，c30064_g1，c33061_g2；TM3亞家族有3個(gè)成員，其他亞家族僅有2個(gè)或1個(gè)成員（表4）。這些單基因簇的同源基因分別與花分生組織發(fā)育、花期調(diào)控、花器官發(fā)育、果實(shí)發(fā)育等重要生物學(xué)過程相關(guān)。

表4 云錦杜鵑成花相關(guān)MADS-box基因鑒定Table 4 Identification of Floral related MADS-box genes of R.fortunei

2.7 SSR分析

利用MISA軟件對(duì)云錦杜鵑轉(zhuǎn)錄組序列進(jìn)行分析，共發(fā)現(xiàn)21900個(gè)SSR位點(diǎn)，分布在17414條單基因簇中，其中有含有1個(gè)以上SSR位點(diǎn)的單基因簇有3606條（表5）。在所有SSR位點(diǎn)中，雙堿基重復(fù)SSR最多，有12294個(gè)，占總數(shù)的56.14%；其次為單堿基重復(fù)SSR，有6448個(gè)，占總數(shù)的29.44%；三堿基重復(fù)SSR有2970個(gè)，占13.56%；四堿基重復(fù)SSR有140個(gè)，占0.63%；五堿基和六堿基重復(fù)SSR分別僅有25個(gè)和23個(gè)（表5）。進(jìn)一步分析這些SSR重復(fù)基序，可以發(fā)現(xiàn)，在單堿基SSR中，A/T發(fā)生頻率最高；雙堿基SSR中，發(fā)生頻率最高的是AG/CT，其次是AC/GT；三堿基重復(fù)中發(fā)生頻率最高是AAG/CTT， 其次是AGG/CCT（圖4）。

3 結(jié)論與討論

中國(guó)是世界上杜鵑花屬植物資源最為豐富的國(guó)家，為世界杜鵑花育種做出了巨大貢獻(xiàn)。但中國(guó)杜鵑花育種尤其常綠杜鵑育種水平較歐美、日本等國(guó)仍有較大差距，資源開發(fā)利用水平低，優(yōu)良品種少。云錦杜鵑作為中國(guó)特有常綠杜鵑，觀賞價(jià)值高、抗性好，野生資源也較為豐富，具有良好的開發(fā)潛力。但由于缺乏相關(guān)的遺傳背景信息，云錦杜鵑的遺傳多樣性、雜交子代鑒定和優(yōu)異基因型挖掘等遺傳育種研究一直受到制約。近年來，隨著高通量測(cè)序技術(shù)的發(fā)展，轉(zhuǎn)錄組測(cè)序已成為非模式生物遺傳背景解析的重要手段，在標(biāo)記開發(fā)、表達(dá)分析、功能基因挖掘等方面得到廣泛應(yīng)用［17-21］。本研究利用Illumina測(cè)序技術(shù)對(duì)云錦杜鵑組培苗的轉(zhuǎn)錄組進(jìn)行測(cè)序和分析，以獲得其轉(zhuǎn)錄組序列信息。云錦杜鵑轉(zhuǎn)錄組的測(cè)序數(shù)據(jù)分析結(jié)果表明：數(shù)據(jù)的Q30值為91.85%，拼接后共獲得84633條單基因簇，平均長(zhǎng)度為691.4 bp，N50值為1177 bp。一般認(rèn)為Q30在80%以上就認(rèn)為測(cè)序質(zhì)量可靠；N50值越大就表示長(zhǎng)片段越多，且不小于800 bp就說明組裝得到序列完整性較好［22］。上述結(jié)果表明本研究測(cè)序數(shù)據(jù)的質(zhì)量和組裝長(zhǎng)度達(dá)到了轉(zhuǎn)錄組分析的基本要求，為進(jìn)一步分析利用奠定了基礎(chǔ)。

表5 云錦杜鵑轉(zhuǎn)錄組SSR分析結(jié)果Table 5 Summary of SSR in R.fortunei transcriptome

圖4 云錦杜鵑SSR類型數(shù)量分布Figure 4 Distribution of SSR motif number of R.fortunei

基因功能注釋是轉(zhuǎn)錄組分析的重要內(nèi)容，是進(jìn)行重要功能基因挖掘的前提。因此，本研究利用Nr和Swiss-prot等七大數(shù)據(jù)庫(kù)對(duì)云錦杜鵑轉(zhuǎn)錄組序列進(jìn)行功能注釋，結(jié)果表明：共有35526條單基因簇獲得注釋信息，仍有約5萬條序列沒有獲得注釋。這與薏苡Coix lachryma-jobi［20］和巖穴蕨Monachosorum maximowiczii［23］的情況類似，可能是由于云錦杜鵑是未測(cè)序物種，在相關(guān)數(shù)據(jù)庫(kù)中缺乏對(duì)應(yīng)的功能注釋信息，也可能是部分云錦杜鵑單基因簇序列本身太短造成的。GO和KOG注釋功能分類的結(jié)果顯示，云錦杜鵑單基因簇的功能涉及了各類生命活動(dòng)；KEGG通路注釋到9887條單基因簇，涉及到272條代謝通路。這些結(jié)果表明：對(duì)于云錦杜鵑等這一非模式植物，轉(zhuǎn)錄組測(cè)序可以有效地解析遺傳背景，獲得大量序列信息?；蚬δ茏⑨屢彩峭诰蚺c特定途徑或功能相關(guān)基因的有效手段。如在紫色黃秋葵Abelmoschus esculentus中，通過轉(zhuǎn)錄組的KEGG注釋，獲得與花色素苷、黃酮、類黃酮、二萜類和萜類骨架等生物合成相關(guān)的單基因簇［24］。本研究通過功能注釋，也鑒定獲得24個(gè)編碼MADS-box基因的單基因簇，屬于10個(gè)不同的亞家族，它們可能與花分生組織發(fā)育、花期調(diào)控、花器官發(fā)育等重要成花過程相關(guān)。

簡(jiǎn)單序列重復(fù)（SSR）又稱微衛(wèi)星序列，具有共顯性、密度大、信息量豐富等優(yōu)勢(shì)，廣泛應(yīng)用于遺傳圖譜構(gòu)建、遺傳多樣性分析、基因定位、分子標(biāo)記輔助育種等方面［25］。利用轉(zhuǎn)錄組序列開發(fā)SSR標(biāo)記具有通量高，成本低的優(yōu)勢(shì)，已在多種植物中獲得成功［19,26-27］。在大王杜鵑R.rex轉(zhuǎn)錄組序列中鑒定獲15314個(gè)SSR位點(diǎn)，占比最高的為雙堿基重復(fù)SSR，其次為單堿基重復(fù)SSR和三堿基重復(fù)SSR，且利用這些SSR位點(diǎn)開發(fā)了相應(yīng)引物對(duì)20份大王杜鵑種質(zhì)進(jìn)行了遺傳多樣性評(píng)價(jià)［3］。本研究也在云錦杜鵑單基因簇序列中鑒定獲得21900個(gè)SSR位點(diǎn)，發(fā)現(xiàn)其中雙核苷酸重復(fù)SSR最多，達(dá)到12294個(gè)；其次為單核苷酸重復(fù)和三核苷酸重復(fù)SSR，這與大王杜鵑中發(fā)現(xiàn)的規(guī)律類似。這些結(jié)果將為云錦杜鵑SSR標(biāo)記開發(fā)提供重要序列信息，也為云錦杜鵑種質(zhì)資源遺傳多樣性分析、功能基因挖掘以及分子輔助育種等工作提供了重要基礎(chǔ)。