基于高通量測序的遼東櫟轉錄組學研究

劉玉林 李偉 張志翔

（1.北京林業大學生物科學與技術學院，北京 100083；2.北京林業大學自然保護區學院，北京 100083）

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

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

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

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

1 材料與方法

1.1 材料

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

1.2 方法

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

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

1.2.3 功能注釋 使用BLASTx程序將拼接所得的unigene與核酸、蛋白質序列數據庫比對（E值<1e-5），并選取最佳注釋。蛋白質數據庫包括Swiss-Prot、GenBank非冗余蛋白數據庫（Nr）、蛋白質直系同源簇數據庫（Clusters of Orthologous Groups of proteins，COGs）、京都基因和基因組百科全書（Kyoto Encyclopedia of Genes and Genomes，KEGG）及Gene Ontology（GO）。其中，unigene通過COG、GO和KEGG數據庫的分類的參考文獻［19] 的方法。

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

2 結果

2.1 Illumina Solexa 測序和序列拼接

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

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

將所獲得的unigene與公共數據Nr和Swiss-Prot進行Blastx比對，通過gene的相似性進行功能注釋，共有57 637條unigene獲得了基因注釋（hit），占總unigene總數的60.16%，而在其余未得到任何一個上述數據庫的注釋（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 功能分類研究

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

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

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

圖2 遼東櫟unigene的COG分類

圖3 遼東櫟unigene的GO分類

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

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

表1 遼東櫟中包含unigene數目最多的10個代謝通路

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

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

2.5 SSR分析

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

3 討論

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

表3 遼東櫟中的SSR位點的數量與分布

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

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

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

4 結論

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

［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] 謝光輝, 郭興強, 王鑫, 等.能源作物資源現狀與發展前景［J] .資源科學, 2007, 29（5）：74-80.

［4] 吳創之, 周肇秋, 陰秀麗, 等.我國生物質能源發展現狀與思考［J] .農業機械學報, 2009, 40（1）：91-99.

［5] 李碧芳.發展生物質能源對能源安全和糧食安全的影響［J] .生態經濟, 2010, 3：41-42.

［6] 方炎明, 張聰穎, 虞木奎, 等.櫟樹繁殖生物技術進展［J] .生物技術通報, 2011（4）：60-65.

［7] 陳婧, 馬履一, 段劼, 等.中國林業生物質能源發展現狀［J] .林業實用技術, 2012, 11：16-19.

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

［9] 李安平, 田玉峰, 謝碧霞, 等.橡實淀粉生料酒精發酵與傳統酒精發酵的能耗和成分組成比較［J] .江西農業大學學報, 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] 劉紅亮, 鄭麗明, 劉青青, 等.非模式生物轉錄組研究［J] .遺傳, 2013, 35（8）：955-970.

［16] 羅偉祥, 郝懷曉, 薛安平.橡樹資源-優質林木生物質能源發展戰略研究［J] .生物質化學工程, 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.