基于高通量測序的遼東櫟轉(zhuǎn)錄組學(xué)研究

劉玉林 李偉 張志翔

（1.北京林業(yè)大學(xué)生物科學(xué)與技術(shù)學(xué)院，北京 100083；2.北京林業(yè)大學(xué)自然保護(hù)區(qū)學(xué)院，北京 100083）

隨著煤、石油、天然氣等傳統(tǒng)能源的枯竭及燃燒所帶來的生態(tài)環(huán)境的惡化，燃料乙醇作為新興的生物質(zhì)能源得到越來越多的國家的重視［1，2]。在許多相對發(fā)達(dá)的國家，生產(chǎn)燃料乙醇主要是以玉米、甜高粱等作物為原材料［3，4]，而對于我國人均耕地面積處于世界平均耕地面積以下的現(xiàn)狀，利用糧食生產(chǎn)燃料乙醇明顯是不現(xiàn)實(shí)的［5]。因此，尋求廉價(jià)易得的“非糧”生物質(zhì)能源的原材料已成為國內(nèi)外學(xué)者競相研究的熱點(diǎn)。

櫟樹（橡樹或柞樹），為殼斗科櫟屬植物的統(tǒng)稱。全世界約有300余種，我國約有60余種［6]。櫟樹種子富含淀粉，含量可達(dá)到50%-70%，而我國約有櫟樹林1 800萬hm2，年產(chǎn)約1 000萬t櫟樹種子，可生產(chǎn)燃料乙醇約250萬t［7]，具有巨大的開發(fā)潛力。若利用櫟樹種子發(fā)展燃料乙醇，不僅在不影響國家糧食安全的情況下發(fā)展燃料乙醇產(chǎn)業(yè)，同時(shí)又能加大櫟樹種子的利用率，刺激林業(yè)經(jīng)濟(jì)的發(fā)展［8，9]。

利用淀粉植物生產(chǎn)燃料乙醇，淀粉含量是一重要指標(biāo)。目前，國內(nèi)外對櫟屬植物中淀粉的研究僅停留在物理與化學(xué)性質(zhì)的分析上［10，11]，還未對櫟屬植物中淀粉合成代謝的分子機(jī)理進(jìn)行深入的研究，進(jìn)而利用分子手段進(jìn)一步提高淀粉的含量和產(chǎn)量。雖然國外學(xué)者利用傳統(tǒng)的測序方法已對夏櫟（Quercus robur）和無?；担≦uercus petraea）進(jìn)行一定層次的研究［12-14]，但大部分研究和所獲得的基因序列來自于芽和葉等組織，還未包括主要積累淀粉的果實(shí)中的轉(zhuǎn)錄數(shù)據(jù)，因而對櫟樹中淀粉合成代謝的研究有一定的局限性。隨著高通量測序技術(shù)的發(fā)展，測序成本的降低，基于高通量測序技術(shù)的轉(zhuǎn)錄組分析越來越成為非模式植物中發(fā)掘功能基因的一種有效的手段。近年來，研究人員已對多種植物進(jìn)行了轉(zhuǎn)錄組測序和基因的發(fā)掘與研究［15]。

本研究利用Illumina Solexa Hiseq 2000高通量測序技術(shù)對廣泛分布于我國的陜西、寧夏、黑龍江、河北、河南、吉林、青海、遼寧、山西、四川、甘肅、內(nèi)蒙古、山東等地及朝鮮地區(qū)的溫帶、暖溫帶森林植被的主要建群種，且種子中淀粉含量高達(dá)62.88%左右［16]的遼東櫟的芽、花、葉和果實(shí)的混合樣品進(jìn)行轉(zhuǎn)錄組測序，旨在獲得遼東櫟更多的轉(zhuǎn)錄本和更為全面的轉(zhuǎn)錄組信息，發(fā)掘遼東櫟中參與淀粉合成代謝的基因和潛在的SSR（Simple sequence Repeat）標(biāo)記，為進(jìn)一步研究遼東櫟中淀粉的合成機(jī)制，選育優(yōu)良樹種及遼東櫟乃至櫟屬植物的進(jìn)化與多樣性分析奠定分子基礎(chǔ)。

1 材料與方法

1.1 材料

本研究中遼東櫟的芽、花、葉與果（5個(gè)階段：開花后20、40、60、80和100 d均采集）均采集于北京市門頭溝區(qū)小龍門國家森林公園（E11526’，N3959’），每份樣品至少取自10個(gè)單株。樣品采集后立即放入液氮中，后存放于-80℃冰箱備用。

1.2 方法

1.2.1 總RNA的提取和Solexa測序 總RNA的提取采用RNeasy Plant Mini Kits（Qiagen，Inc.，Valencia，CA，USA）試劑盒按照實(shí)驗(yàn)手冊操作步驟提取。來自芽、花、葉和果的等量RNA混合后，取10 μg用于cDNA文庫構(gòu)建，cDNA文庫的構(gòu)建參考文獻(xiàn)［17] 的方法。利用Illumina Solexa Hiseq 2000的paired-end測序方法進(jìn)行轉(zhuǎn)錄組測序。

1.2.2 序列的處理與拼接 鑒于Solexa數(shù)據(jù)錯(cuò)誤率對結(jié)果的影響，對原始數(shù)據(jù)進(jìn)行質(zhì)量預(yù)處理。先利用滑動(dòng)窗口法去除低質(zhì)量片段：質(zhì)量閾值20（錯(cuò)誤率=1%），窗口大小5 bp，長度閾值35 bp；再切除reads中含N部分序列：長度閾值35 bp；最后利用軟件Trinity對遼東櫟有效reads進(jìn)行從頭拼接［18]。

1.2.3 功能注釋 使用BLASTx程序?qū)⑵唇铀玫膗nigene與核酸、蛋白質(zhì)序列數(shù)據(jù)庫比對（E值<1e-5），并選取最佳注釋。蛋白質(zhì)數(shù)據(jù)庫包括Swiss-Prot、GenBank非冗余蛋白數(shù)據(jù)庫（Nr）、蛋白質(zhì)直系同源簇?cái)?shù)據(jù)庫（Clusters of Orthologous Groups of proteins，COGs）、京都基因和基因組百科全書（Kyoto Encyclopedia of Genes and Genomes，KEGG）及Gene Ontology（GO）。其中，unigene通過COG、GO和KEGG數(shù)據(jù)庫的分類的參考文獻(xiàn)［19] 的方法。

1.2.4 SSR位點(diǎn)的篩選 利用MISA軟件在所有unigene 中搜索SSR 位點(diǎn)，參數(shù)設(shè)置如下：二核苷酸至少重復(fù)次數(shù)為6，三核苷酸、四核苷酸、五核苷酸和六核苷酸至少重復(fù)次數(shù)均為5。

2 結(jié)果

2.1 Illumina Solexa 測序和序列拼接

采用 Illumina Solexa Hiseq 2000 測序技術(shù)對遼東櫟的芽、花、葉和果實(shí)的混合樣進(jìn)行轉(zhuǎn)錄組測序，共獲得46 400 862條原始序列，總長4.64 Gb。經(jīng)過預(yù)處理，最終得到了有效序列40 621 588條，數(shù)據(jù)量為3.8 Gb，平均長度為94.95 bp。使用軟件Trinity對有效reads進(jìn)行從頭拼接最終得到了151 339個(gè)長度大于200 bp的contig，總長度約為130.92 Mb，最大長度、平均長度以及N50分別為11 284、865和1 442 bp。取每個(gè)contig下最長的轉(zhuǎn)錄本作為unigene，得到了95 800個(gè)unigene，總長度約為73.57 Mb，平均長度與N50分別為768 bp和1 373 bp。其中，大于2 000 bp的序列共有7 975條，占unigene總數(shù)的8.32%（圖1-A）。

2.2 序列的比對、功能注釋及unigene的特征分析

將所獲得的unigene與公共數(shù)據(jù)Nr和Swiss-Prot進(jìn)行Blastx比對，通過gene的相似性進(jìn)行功能注釋，共有57 637條unigene獲得了基因注釋（hit），占總unigene總數(shù)的60.16%，而在其余未得到任何一個(gè)上述數(shù)據(jù)庫的注釋（no-hit）38 163條unigene（39.84%）中，有37 752條unigene（98.92%）小于或等于1 000 bp（圖1-B）。

圖1 contig和unigene的長度分布（A）及有注釋和無注釋的unigene的分布比例（B）

2.3 功能分類研究

為進(jìn)一步研究遼東櫟中unigene的功能分類，將所得到的95 800條unigene在COG和GO數(shù)據(jù)庫中進(jìn)行比對及功能注釋與分類。

在COG分類中，共有36 407 條unigene（占unigene總數(shù)的38%）被注釋到24個(gè)COG類別中。其中，“一般功能基因”是最大類別，包含8 330 條unigene，占被注釋到unigene總數(shù)的22.88%；其次是“蛋白質(zhì)翻譯后修飾與轉(zhuǎn)運(yùn)，分子伴侶”，包含4 128條unigene；而“核酸結(jié)構(gòu)”是最小的類別，僅包含11條unigene。此外，有2 491條unigene參與了碳水化合物運(yùn)輸與代謝（圖2）。

利用GO對獲得遼東櫟unigene進(jìn)行功能分類，共有43 766條unigene被注釋到生物學(xué)過程、細(xì)胞組分和分子功能3 個(gè)大類別中。其中，36 372條unigene歸入生物學(xué)過程，18 876 條unigene歸入細(xì)胞組分以及40 358條unigene歸入分子功能。3個(gè)大的類別又被劃分為45個(gè)小的類別（圖3）。代謝過程是生物學(xué)過程中的最大類別，包含31 427條unigene；細(xì)胞是細(xì)胞組分中的最大類別，包含12 764條unigene；蛋白結(jié)合是分子功能中的最大類別，包含28 892條unigene。

圖2 遼東櫟unigene的COG分類

圖3 遼東櫟unigene的GO分類

2.4 代謝途徑分析和淀粉合成基因的篩選

將遼東櫟的unigene序列映射到KEGG數(shù)據(jù)庫的參考代謝通路中，共有11 468條unigene參與到185個(gè)代謝通路中。其中包含unigene最多的是代謝通路是剪接體（ko03010），共有1 161條unigene。其次是內(nèi)質(zhì)網(wǎng)上的蛋白加工（ko04141），包含630條unigene。而參與淀粉與蔗糖的代謝通路（ko00500）的unigene共有434條（表1）。

表1 遼東櫟中包含unigene數(shù)目最多的10個(gè)代謝通路

結(jié)合KEGG數(shù)據(jù)庫中參考pathway的關(guān)于淀粉與蔗糖代謝通路中發(fā)掘到的unigene及在其他公共數(shù)據(jù)庫中的注釋，共統(tǒng)計(jì)篩選出67條參與淀粉合成的unigene，編碼9個(gè)關(guān)鍵酶，其中5個(gè)unigene編碼α-糖苷酶；17個(gè)unigene編碼β-呋喃果糖苷酶；9個(gè)unigene編碼己糖激酶；10個(gè)unigene編碼果糖激酶；4個(gè)unigene編碼葡萄糖-6-磷酸異構(gòu)酶；4個(gè)unigene編碼葡萄糖磷酸變位酶；8個(gè)unigene編碼葡萄糖-1-磷酸腺苷酰基轉(zhuǎn)移酶；6個(gè)unigene編碼淀粉合成酶及4個(gè)unigene編碼1，4-α-葡聚糖分支酶（表2）。

表2 遼東櫟中參與淀粉合成的酶

2.5 SSR分析

利用MISA軟件在遼東櫟的13 380條unigene中共搜索到15 901個(gè)SSR位點(diǎn)，占unigene總序列數(shù)的13.97%，平均每1.28 kb出現(xiàn)1個(gè)SSR，其中包含有兩個(gè)及兩個(gè)以上SSR的unigene共有2 521條。二核苷酸和三核苷酸重復(fù)類型分別占SSR總數(shù)的60.07%和38.09%；四核苷酸重復(fù)類型占SSR總數(shù)的1.57%；而五核苷酸和六核苷酸重復(fù)類型在遼東櫟中轉(zhuǎn)錄組序列中含量較少，僅占SSR總數(shù)的0.19%和0.09%。除此之外，不同核苷酸的重復(fù)次數(shù)也有很大的變化（表3）。

3 討論

基于高通量測序技術(shù)的轉(zhuǎn)錄組學(xué)研究是一種非常高效、可靠的發(fā)掘功能基因手段，且在淀粉的合成及代謝中也已得到廣泛應(yīng)用［20，21]。目前，高通量測序技術(shù)主要是Roche公司的454測序技術(shù)、Illumina公司和ABI公司相繼推出的Solexa和SOLid測序技術(shù)。其中，Solexa測序技術(shù)相對于其他兩種技術(shù)在測序成本和數(shù)據(jù)量輸出方面更具優(yōu)勢［22]。前人的研究表明，不同組織的混合取樣，可在節(jié)約試驗(yàn)成本的基礎(chǔ)上發(fā)掘到更多的轉(zhuǎn)錄本［23]。因此，本研究采用遼東櫟的芽、花、葉和果實(shí)的混合樣品進(jìn)行轉(zhuǎn)錄組測序。櫟屬植物的基因組大小約為539-921 Mb［24]，本研究共獲得3.8 Gb的有效數(shù)據(jù)量，約覆蓋遼東櫟基因組的4.13-7.05倍。同時(shí)利用在對Illumina Solexa Hiseq 2000測序數(shù)據(jù)進(jìn)行拼接過程中表現(xiàn)非常優(yōu)異的Trinity軟件［25]對本研究的數(shù)據(jù)進(jìn)行拼接處理，共獲得95 800條unigene，其中有38 163條unigene未在Blastx同源性搜索中得到注釋，但大多數(shù)片段小于1 000 bp（37 752條，占98.92%），因此這些片段可能是由于較短而未與公共數(shù)據(jù)庫中的序列比對上，也可能是短的非編碼序列或者是新的基因。

表3 遼東櫟中的SSR位點(diǎn)的數(shù)量與分布

COG和GO的功能分類對初步了解基因的功能起著重要作用，而KEGG數(shù)據(jù)庫中的參考pathway不僅可以推測基因的功能，而且可以研究基因在不同代謝通路中所在位置及作用，三者相輔相成，成為新物種中發(fā)掘功能基因的重要手段［19]。本研究通過KEGG 數(shù)據(jù)庫中的pathway分析篩選與淀粉合成相關(guān)的基因，同時(shí)結(jié)合所篩選到的unigene在本研究中其他數(shù)據(jù)庫中的功能注釋，從而進(jìn)一步確保了所獲得的基因的可靠性。

鑒于分布廣泛和多態(tài)性較高的特性，SSR標(biāo)記已被廣泛應(yīng)用于分子標(biāo)記輔助選擇育種（Molecular marker-assisted selection，MAS）和利用分子標(biāo)記通過關(guān)聯(lián)分析（Association mapping）發(fā)掘與好的農(nóng)藝性狀連鎖緊密的相關(guān)基因［26]，尤其是對于多年生的植物，可大大縮短育種年限［27]。本研究在遼東櫟中發(fā)掘到15 901個(gè)SSR位點(diǎn)，其中二核甘酸和三核苷酸的重復(fù)占到總數(shù)的98.16%，為了保證SSR位點(diǎn)的潛在多態(tài)性，我們在篩選過程中對于四、五和六核苷酸的最小重復(fù)次數(shù)同樣設(shè)置為5，一定程度上影響到了這3類核苷酸重復(fù)在總SSR位點(diǎn)中所占比例。

本研究通過高通量的測序，獲得了大量的遼東櫟轉(zhuǎn)錄組序列，不僅豐富了遼東櫟的基因庫，而且為遼東櫟及其他櫟類淀粉合成基因的克隆與功能研究奠定了基礎(chǔ)。發(fā)掘到的SSR位點(diǎn)可通過進(jìn)一步的開發(fā)為遼東櫟的進(jìn)化與多樣性分析及遼東櫟的遺傳圖譜的構(gòu)建和QTL（Quantitative Trait Loci）的定位提供了數(shù)據(jù)支持。

4 結(jié)論

本研究通過Solexa Hiseq 2000高通量測序，獲得3.8 Gb的遼東櫟轉(zhuǎn)錄組序列，拼接獲得95 800條unigene，發(fā)掘出67條參與淀粉合成的unigene以及15 901個(gè)SSR位點(diǎn)。

［1] Balat M, Balat H. Recent trends in global production and utilization of bio-ethanol fuel ［J] . Applied Energy, 2009, 86：2273-2282.

［2] Hahn-H?gerdal B, Galbe M, Gorwa-Grauslund MF, et al. Bioethanol-the fuel of tomorrow from the residues of today［J] . Trends in Biotechnology, 2006, 24（12）：549-556.

［3] 謝光輝, 郭興強(qiáng), 王鑫, 等.能源作物資源現(xiàn)狀與發(fā)展前景［J] .資源科學(xué), 2007, 29（5）：74-80.

［4] 吳創(chuàng)之, 周肇秋, 陰秀麗, 等.我國生物質(zhì)能源發(fā)展現(xiàn)狀與思考［J] .農(nóng)業(yè)機(jī)械學(xué)報(bào), 2009, 40（1）：91-99.

［5] 李碧芳.發(fā)展生物質(zhì)能源對能源安全和糧食安全的影響［J] .生態(tài)經(jīng)濟(jì), 2010, 3：41-42.

［6] 方炎明, 張聰穎, 虞木奎, 等.櫟樹繁殖生物技術(shù)進(jìn)展［J] .生物技術(shù)通報(bào), 2011（4）：60-65.

［7] 陳婧, 馬履一, 段劼, 等.中國林業(yè)生物質(zhì)能源發(fā)展現(xiàn)狀［J] .林業(yè)實(shí)用技術(shù), 2012, 11：16-19.

［8] 田玉峰, 李安平, 謝碧霞, 等.橡實(shí)淀粉生物乙醇化橡實(shí)品種和菌種的篩選［J] .食品科學(xué), 2011, 32（7）：207-210.

［9] 李安平, 田玉峰, 謝碧霞, 等.橡實(shí)淀粉生料酒精發(fā)酵與傳統(tǒng)酒精發(fā)酵的能耗和成分組成比較［J] .江西農(nóng)業(yè)大學(xué)學(xué)報(bào), 2012,34（5）：1032-1038.

［10] Soni PL, Sharma H, Dun D, et al. Physicochemical properties ofQuercus leucotrichophora（Oak）starch ［J] . Starch/St?rke, 1993,45（4）：127-130.

［11] Stevenson DG, Jane JL, Inglett GE. Physicochemical properties of pin oak（Quercus palustrisMuenchh.）acorn starch ［J] . Starch/St?rke, 2006, 58（11）：553-560.

［12] Derory J, Léger P, Garcia V, et al. Transcriptome analysis of bud burst in sessile oak（Quercus petraea）［J] . New Phytologist,2006, 170（4）：723-738.

［13] Lesur I, Durand J, Sebastiani F, et al. A sample view of the pedunculate oak（Quercus robur）genome from the sequencing of hypomethylated and random genomic libraries ［J] . Tree Genetics& Genomes, 2011, 7（6）：1277-1285.

［14] Ueno S, Le Provost G, Léger V, et al. Bioinformatic analysis of ESTs collected by Sanger and pyrosequencing methods for a keystone forest tree species：oak ［J] . BMC Genomics, 2010, 11：650.

［15] 劉紅亮, 鄭麗明, 劉青青, 等.非模式生物轉(zhuǎn)錄組研究［J] .遺傳, 2013, 35（8）：955-970.

［16] 羅偉祥, 郝懷曉, 薛安平.橡樹資源-優(yōu)質(zhì)林木生物質(zhì)能源發(fā)展戰(zhàn)略研究［J] .生物質(zhì)化學(xué)工程, 2006, 40（B12）：147-152.

［17] Huang LL, Yang X, Sun P, et al. The first Illumina-based de novo transcriptome sequencing and analysis of safflower flowers ［J] .PloS One, 2012, 7（6）：e38653.

［18] Haas BJ, Papanicolaou A, Yassour M, et al.De novotranscript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis ［J] . Nature Protocols, 2013,8（8）：1494-1512.

［19] Wang R, Xu S, et al.De novosequence assembly and characterization ofLycoris aureatranscriptome using GS FLX Titanium platform of 454 pyrosequencing ［J] . PloS One, 2013, 8（4）：e60449.

［20] Tao X, Fang Y, Xiao Y, et al. Comparative transcriptome analysis to investigate the high starch accumulation of duckweed（Landoltia punctata）under nutrient starvation ［J] . Biotechnology for Biofuels, 2013, 6（1）：72.

［21] Kaminski KP, Petersen AH, S?nderk?r M, et al. Transcriptome analysis suggests that starch synthesis may proceed via multiple metabolic routes in high yielding potato cultivars ［J] . PloS One,2012, 7（12）：e51248.

［22] Shu S, Chen B, et al.De novosequencing and transcriptome analysis ofWolfiporia cocosto reveal genes related to biosynthesis of triterpenoids ［J] . PloS One, 2013, 8（8）：e71350.

［23] Zhou Y, Gao F, Liu R, et al.De novosequencing and analysis of root transcriptome using 454 pyrosequencing to discover putative genes associated with drought tolerance inAmmopiptanthus mongolicus［J] . BMC Genomics, 2012, 13：266.

［24] Kole, Chittaranjan. Fagaceae Trees. //In genome mapping and molecular breeding in plants［M] . Springer, 2007, 7：161-184.

［25] Li SW, Yang H, Liu YF, et al. Transcriptome and gene expression analysis of the rice leaf folder,Cnaphalocrosis medinalis［J] . PloS One, 2012, 7（11）：e47401.

［26] O’Malley DM, McKeand SE. Marker assisted selection for breeding value in forest trees ［J] . Forest Genetics, 1994, 1：207-218.

［27] Neale DB, Kremer A. Forest tree genomics：growing resources and applications ［J] . Nature Reviews Genetics, 2011, 12：111-122.