劉玉林 李偉 張志翔

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

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

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

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

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

1 材料與方法

1.1 材料

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

1.2 方法

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

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

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

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

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

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

圖1 contig和unigene的長(zhǎng)度分布（A）及有注釋和無(wú)注釋的unigene的分布比例（B）

2.3 功能分類(lèi)研究

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

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

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

圖2 遼東櫟unigene的COG分類(lèi)

圖3 遼東櫟unigene的GO分類(lèi)

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

將遼東櫟的unigene序列映射到KEGG數(shù)據(jù)庫(kù)的參考代謝通路中，共有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ù)庫(kù)中參考pathway的關(guān)于淀粉與蔗糖代謝通路中發(fā)掘到的unigene及在其他公共數(shù)據(jù)庫(kù)中的注釋，共統(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-磷酸腺苷?；D(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ù)類(lèi)型分別占SSR總數(shù)的60.07%和38.09%；四核苷酸重復(fù)類(lèi)型占SSR總數(shù)的1.57%；而五核苷酸和六核苷酸重復(fù)類(lèi)型在遼東櫟中轉(zhuǎn)錄組序列中含量較少，僅占SSR總數(shù)的0.19%和0.09%。除此之外，不同核苷酸的重復(fù)次數(shù)也有很大的變化（表3）。

3 討論

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

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

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

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

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

4 結(jié)論

本研究通過(guò)Solexa Hiseq 2000高通量測(cè)序，獲得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)之, 周肇秋, 陰秀麗, 等.我國(guó)生物質(zhì)能源發(fā)展現(xiàn)狀與思考［J] .農(nóng)業(yè)機(jī)械學(xué)報(bào), 2009, 40（1）：91-99.

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

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

［7] 陳婧, 馬履一, 段劼, 等.中國(guó)林業(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] 羅偉祥, 郝懷曉, 薛安平.橡樹(shù)資源-優(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.