云錦杜鵑轉錄組分析

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

（1.浙江農林大學 林業與生物技術學院，浙江 杭州311300；2.浙江農林大學 省部共建亞熱帶森林培育國家重點實驗室，浙江 杭州311300）

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

1 材料與方法

1.1 植物材料

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

1.2 RNA提取與質控

采用Trizol試劑（Invitrogen），按照試劑說明書分別提取組培苗根、莖、葉的總RNA。采用Agilent 2100 Bioanalyzer（Agilent），對總RNA樣品的純度、濃度和完整性進行檢測評估，等量混合后用于測序文庫構建。

1.3 測序文庫構建

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

1.4 轉錄組數據分析

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

1.5 基因功能注釋

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

1.6 SSR分析

采用MISA軟件對云錦杜鵑單基因簇進行簡單序列重復（SSR）分析，鑒定其中6種類型的SSR；各SSR類型重復次數設定為：單核苷酸SSR，重復數≥10次；二核苷酸SSR，重復數≥6次；三至六核苷酸SSR，重復數≥5次。

2 結果與分析

2.1 測序結果與組裝

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

測序數據經處理后得到94252430個干凈測序片段（clean reads），包含了14.14 Gb核苷酸序列信息，GC含量為47.58%。測序質量評估結果顯示，堿基錯誤率為0.02%，Q20為96.70%，Q30為91.85%。這些表明該轉錄組的數據量和質量均較高，為后續的序列組裝提供了高質量的原始數據。

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

表1 云錦杜鵑轉錄組組裝結果Table 1 Summary of transcriptome assembly for R.fortunei

2.2 基因功能注釋

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

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

在Nr庫中，云錦杜鵑轉錄組注釋到其他物種的單基因簇基因序列共30700條。其中與葡萄Vitis vinifera基因序列相似的最多，所占比例為24.39%；其次為中粒咖啡Coffea canephora，所占比例為7.39%；第3為可可Theobroma cacao，占6.77%；其他相似性序列數量大于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數據庫的注釋，共有23215條單基因簇獲注釋信息，得到119389個GO功能注釋。由分類結果可知：生物學過程（biological process）最多55692條，占46.65%，其次是細胞組分（cellular component），35577條，占29.80%，最少的分子功能（molecular function），有28120條，占23.55%；這三大功能分類又可分為55個亞類，其中生物學過程23個亞類，細胞組分18個亞類，分子功能14個亞類（圖2）。生物學過程中，涉及細胞過程、代謝過程和單一有機體進程的單基因簇較多，分別有12698，12376和9581條；細胞組分中涉及較多的是細胞、細胞部分和大分子復合體，分別有7216，7214和4713條；分子功能中涉及較多的有結合功能和催化活性，分別有13302，10524條（圖2）。

圖1 云錦杜鵑單基因簇Nr數據庫比對相似性物種分布圖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數據庫對云錦杜鵑單基因簇進行注釋，結果顯示有11085條序列獲得12475個注釋信息，可分為26個功能分類（圖3）。從基因功能分類來看，涉及一般功能預測的序列最多，多達1970條；涉及翻譯后修飾、蛋白翻轉、分子伴侶功能的序列次之，有1527條；而涉及核結構、胞外結構和細胞運動的序列很少，僅有37，34和5條（圖3）。

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

2.5 KEGG通路注釋

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

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

2.6 MADS-box基因家族分析

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

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

2.7 SSR分析

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

3 結論與討論

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

表5 云錦杜鵑轉錄組SSR分析結果Table 5 Summary of SSR in R.fortunei transcriptome

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

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

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