基于高通量測序的慈竹筍轉(zhuǎn)錄組分析與基因功能注釋

陳宇鵬，曹穎，胡尚連，黃艷，盧學琴，徐剛，龍治堅

西南科技大學 生命科學與工程學院 四川省生物質(zhì)資源利用與改性工程技術(shù)研究中心，四川 綿陽 621010

基于高通量測序的慈竹筍轉(zhuǎn)錄組分析與基因功能注釋

陳宇鵬，曹穎，胡尚連，黃艷，盧學琴，徐剛，龍治堅

西南科技大學 生命科學與工程學院 四川省生物質(zhì)資源利用與改性工程技術(shù)研究中心，四川 綿陽 621010

陳宇鵬, 曹穎, 胡尚連, 等. 基于高通量測序的慈竹筍轉(zhuǎn)錄組分析與基因功能注釋. 生物工程學報, 2016, 32(11): 1610-1623.

Chen YP, Cao Y, Hu SL, et al. Transcriptome analysis and gene function annotation of Bambusa emeiensis shoots based on high-throughput sequencing technology. Chin J Biotech, 2016, 32(11): 1610-1623.

慈竹是我國四川當?shù)氐膬?yōu)勢叢生竹種之一，其纖維長度和質(zhì)量較優(yōu)異，是造紙、紡織等工業(yè)的良好原料。本文利用Illumina HiSeqTM2000平臺，對10、50、100和150 cm高的慈竹筍進行轉(zhuǎn)錄組分析，共得到69.28 M條讀長 (Reads)，經(jīng)從頭拼接、組裝和聚類后得到111 137條非重復(fù)序列基因Unigene，其中共有63 094條注釋到COG、GO、KEGG、Swiss-Prot和Nr數(shù)據(jù)庫中。這些Unigene不僅具有一般的功能，如轉(zhuǎn)錄和信號轉(zhuǎn)導等，還涉及到蔗糖轉(zhuǎn)運與代謝、次級代謝產(chǎn)物及細胞壁的生物合成等方面。不同高度慈竹筍的纖維素合成酶基因存在差異表達，發(fā)現(xiàn)了可能調(diào)控慈竹生長發(fā)育以及纖維素和木質(zhì)素生物合成的相關(guān)基因，為慈竹品種改良提供一定的理論基礎(chǔ)。

慈竹，轉(zhuǎn)錄組，高通量測序，功能注釋

竹子是重要的非木材可再生資源，以生長速度快、產(chǎn)量高、易加工、用途廣泛、經(jīng)濟價值高等特點著稱。近年來隨著分子生物學研究手段的發(fā)展，尤其是基因組和轉(zhuǎn)錄組測序技術(shù)的廣泛應(yīng)用，毛竹、麻竹等竹種的基因組和 (或)轉(zhuǎn)錄組數(shù)據(jù)已經(jīng)被公布，一些與竹筍快速生長以及竹材纖維形成的相關(guān)基因被報道。如Peng等[1-2]通過基因組和轉(zhuǎn)錄組測序發(fā)現(xiàn)了19個毛竹纖維素合酶基因，和大量與植物激素、細胞周期調(diào)控、細胞壁合成、細胞形態(tài)等相關(guān)的候選基因，為毛竹竹材特性的改良提供了依據(jù)。Liu等[3]分析了麻竹的轉(zhuǎn)錄組，篩選出105個木質(zhì)素生物合成途徑相關(guān)關(guān)鍵酶基因，并與毛竹的相關(guān)基因進行了比對，確定了潛在的與生長和發(fā)育相關(guān)的候選基因。

慈竹Bambusa emeiensis，是四川本土大型叢生竹之一，其稈徑較小，稈壁較薄，纖維長度和質(zhì)量均優(yōu)于其他竹種，是造紙、紡織等行業(yè)的優(yōu)良原料之一，具有比較重要的應(yīng)用前景[4]。目前，對于慈竹的研究主要集中于生長狀況[5]、基因克隆[6]、遺傳多樣性[7]、生理生化指標[8]等方面。與毛竹、麻竹等竹種相比，慈竹的分子生物學研究相對滯后。有研究表明慈竹可能有70條染色體，為典型的多倍體復(fù)合體[9]，其基因組測序的難度較大，而且目前慈竹轉(zhuǎn)錄組分析以及纖維素、木質(zhì)素生物合成相關(guān)的分子機理等也未見報道。

本文以不同高度的慈竹筍為材料，比較了不同高度慈竹筍的纖維素含量，并采用Illumina HiSeqTM2000平臺對其進行了轉(zhuǎn)錄組測序。將得到的Unigene進行Cluster of Orthologous Groups of proteins (COG)、Gene Ontology (GO)、Kyoto Encyclopedia of Genes and Genomes (KEGG)、Swiss-Prot和Nonredundant protein (Nr)注釋，并與其他相似物種進行比對，旨在發(fā)現(xiàn)調(diào)控慈竹生長發(fā)育以及與纖維素、木質(zhì)素生物合成相關(guān)基因，挖掘慈竹中的珍貴基因資源。同時，本文也進一步分析了4個樣品中與纖維素合成相關(guān)的差異表達基因，為慈竹的分子遺傳育種以及品種改良等提供一定的理論基礎(chǔ)。

1 材料與方法

1.1 材料

供試材料采自西南科技大學生命科學與工程學院資源圃 (年平均氣溫17.2 ℃，年平均降雨量793.5 mm)。采樣時，以露出地面10、50、100和150 cm的筍為對象，每個高度取3株獨立的筍 (測量誤差小于2 cm)，并取其基部第2個節(jié)間的部分，分別切碎、混勻，其中一部分儲存于-80 ℃冰箱中備用，另一部分置于烘箱中烘干用于纖維素含量測定。

1.2 不同高度慈竹筍纖維素含量測定

纖維素含量的測定參照FOSS公司FibertecTMM6 1020/1021型纖維素測定儀的操作手冊進行。通過酸性洗滌纖維 (ADF)[10]法，得到纖維素含量，每個樣品生物重復(fù)3次，每組樣品運行1個空白對照。

1.3 轉(zhuǎn)錄組測序

1.3.1 總RNA瓊脂糖凝膠電泳檢測及測序

RNA提取采用OMEGA BIO-TEK公司Plant RNA Kit試劑盒中的試劑及提取方法?？俁NA樣品采用Nanodrop檢測和Agilent Technologies 2100 Bioanalyzer檢測。檢測合格后用NEBNext Poly (A) mRNA Magnetic Isolation Module (NEB, E7490) 富集mRNA，將mRNA片段化處理，以mRNA片段為模板，采用隨機引物法，用NEBNext mRNA Library Prep Master Mix Set for Illumina (NEB, E6110) 和NEBNext Multiplex Oligos for Illumina (NEB, E7500) 構(gòu)建上機文庫。制備好的文庫用1.8%瓊脂糖凝膠電泳進行檢測文庫插入片段大小，然后用Library Quantification Kit-Illumina GA Universal (Kapa, KK4824) 進行QPCR定量，再在cDNA片段兩端加上接頭。檢測合格的文庫在Illumina cbot上進行簇的生成，最后用Illumina HiSeqTM2000進行測序。

1.3.2 數(shù)據(jù)分析

轉(zhuǎn)錄組測序得到的原始數(shù)據(jù)經(jīng)過去除雜質(zhì)和冗余處理后，利用Trinity軟件[11]對經(jīng)過過濾后的高質(zhì)量數(shù)據(jù)進行從頭拼接，獲得重疊群Contig。根據(jù)其結(jié)果利用雙末端信息將來自同一轉(zhuǎn)錄本的不同Contig連在一起，做進一步的序列拼接，得到轉(zhuǎn)錄本Transcripts。在Transcripts聚類單元中選取最主要的轉(zhuǎn)錄本作為Unigene序列，并對Unigene數(shù)據(jù)進行聚類分析和進一步去冗余處理，最終得到非冗余Unigene庫。對于多樣品的組裝，由于后續(xù)的Open reading frame (ORF) 預(yù)測、Simple sequence repeats (SSR) 分析、表達豐度分析等都建立在同一套參考基因的基礎(chǔ)上，因此對各樣品得到的Unigene作進一步的聚類分析，最后得到慈竹的Unigene數(shù)據(jù)庫。

1.3.3 功能注釋

使用Basic Local Alignment Search Tool (Blast) 將組裝后的Unigene序列與COG、GO、KEGG、Swiss-Prot和Nr數(shù)據(jù)庫比對，得到最高序列相似性的Unigene。

1.3.4 基因的差異表達分析

在NR、Swiss-Rort注釋結(jié)果中進行檢索，記錄“GeneID”以及相對應(yīng)的注釋結(jié)果。根據(jù)“GeneID”記錄相應(yīng)基因的相對表達豐度值，并將值導入MeV 4.9.0分析軟件中，制作熱圖并進行聚類分析。

1.4 差異表達基因qRT-PCR驗證

利用primer premier 5.0對已測序正確的序列和慈竹內(nèi)參基因Tublin進行qRT-PCR引物的設(shè)計。在iQ Multicolor Real-Time PCR自動擴增儀進行反應(yīng)，每個樣品重復(fù)3次，反應(yīng)體系20 μL。以內(nèi)參基因Tublin為對照，利用公式2-CT△△計算其相對表達量。

表1 測序數(shù)據(jù)統(tǒng)計Table 1 Statistics of sequencing data

圖1 慈竹功能基因長度分布和所占比例Fig. 1 Distribution and percentage of Unigene length from Bambusa emeiensis.

2 結(jié)果與分析

2.1 測序數(shù)據(jù)結(jié)果統(tǒng)計

利用Illumina HiSeqTM2000平臺，對供試材料進行轉(zhuǎn)錄組測序，共獲得13.98 Gb的原始數(shù)據(jù)，得到69.28 M條的讀長 (reads) (表1)。其中10、50和150 cm慈竹筍總堿基數(shù)均大于4 Gb，reads數(shù)大于20 M，GC含量分別為52.78%、53.12%和52.21%；而100 cm慈竹筍因測序深度較低，僅獲得7 829 716條reads，GC含量為49.49%。同時，供試材料的不確定堿基數(shù)(N%) 所占比例均≤0.03%，測序質(zhì)量值≥30的堿基(Q 30%) 所占比例均>80%，這表明測序質(zhì)量可靠。

通過進一步分析后得到的非冗余Unigene庫，含有111 137條慈竹筍的Unigene。其中200 bp-300 bp之間的Unigene，為41 814條，占37.62%；大于1 000 bp的Unigene為21 518條，占19.35% (圖1A，1B)。

2.2 基因功能注釋

2.2.1 基因功能注釋

通過數(shù)據(jù)庫比對，得到最高序列相似性的Unigene，其中有16 151條Unigene被匹配到COG數(shù)據(jù)庫中；52 872條Unigene被匹配到GO數(shù)據(jù)庫中；11 032條Unigene被匹配到KEGG數(shù)據(jù)庫中；44 552條Unigene被匹配到Swiss-Port數(shù)據(jù)庫中；而62 795條Unigene被匹配到Nr數(shù)據(jù)庫中 (表2)。

表2 慈竹轉(zhuǎn)錄組功能注釋Table 2 Transcriptome function annotation of Bambusa emeiensis

圖2 Unigene的COG功能分類Fig. 2 Unigene of COG function categories.

2.2.2 功能基因的COG注釋

COG的功能注釋與基因產(chǎn)物有直接的聯(lián)系[12]。在慈竹筍Unigene數(shù)據(jù)庫的COG功能注釋中，有16 151條Unigene具有蛋白質(zhì)功能定義。在COG功能分類中包含復(fù)制、重組、修復(fù)、信號轉(zhuǎn)導、新陳代謝、分子過程等25個分類。其中一般功能基因代表最大的一類，有4 221條Unigene；其次是復(fù)制、重組和修復(fù)，有2 701條。此外，還涉及到多個植物生長發(fā)育相關(guān)以及纖維素、木質(zhì)素生物合成相關(guān)的生理生化過程，如碳水化合物轉(zhuǎn)運與代謝、氨基酸轉(zhuǎn)運與代謝、輔酶轉(zhuǎn)運與代謝、次級代謝產(chǎn)物生物合成、運輸和分解代謝、細胞壁與細胞膜的生物合成、信號轉(zhuǎn)導機制等過程 (圖2)。

2.2.3 功能基因的GO分類

GO數(shù)據(jù)庫可對基因和蛋白質(zhì)進行分類和注釋，其對很多物種都是適用的[13]。GO數(shù)據(jù)庫包含3大類，即細胞組分、分子功能和生物學過程。在慈竹筍Unigene數(shù)據(jù)庫中，有52 872條基因注釋到GO數(shù)據(jù)庫中。其中第一大類中的細胞部分、細胞和細胞器所占比例較高，分別達到85.71%、84.85%和80.39%；第二大類中的結(jié)合和催化活性所占比例較高，分別達到68.90%和50.82%；而第三大類中的細胞過程和代謝過程所占比例較高，分別達到75.86%和73.31% (圖3)。

2.2.4 功能基因的KEGG注釋

KEGG數(shù)據(jù)庫是有關(guān)基因產(chǎn)物代謝通路的主要數(shù)據(jù)庫，基于代謝通路的分析對于進一步分析解讀基因的功能非常有幫助[14]。在慈竹筍Unigene數(shù)據(jù)庫中，有11 032條Unigene注釋到KEGG數(shù)據(jù)庫中。其中注釋序列較多的3條代謝通路為RNA轉(zhuǎn)運、植物激素信號轉(zhuǎn)導與核糖體。而與植物纖維素生物合成相關(guān)的代謝通路，如淀粉和蔗糖代謝、光合生物的固碳作用和光合作用分別有257條、189條和105條注釋到KEGG數(shù)據(jù)庫中；與木質(zhì)素生物合成相關(guān)的代謝途徑，如苯丙烷類的生物合成和苯丙氨酸代謝分別有118條和110條在KEGG代謝通路中得到注釋 (表3)。

圖3 Unigene的GO功能分類Fig. 3 Unigene of GO function categories.

表3 Unigene數(shù)目最多的5個代謝通路及與纖維素、木質(zhì)素生物合成相關(guān)的代謝通路Table 3 Top five metabolic pathways and the metabolic pathways of cellulose and lignin biosynthesis involving Bambusa emeiensis Unigene

2.3 慈竹Unigene與其他物種的比較分析

首先將慈竹的Unigene與水稻Oryza sativa、毛竹Phyllostachys edulis、短柄草Brachypodium、玉米Zea mays和高粱Sorghum的Unigene進行功能比對。所得結(jié)果如表4所示，慈竹Unigene與毛竹匹配度最高，達到52.46%，可能是因為慈竹與毛竹同屬竹類植物的原因。其次是與水稻相匹配達到50.41%，可能是因為水稻為模式植物，其基因序列信息較其他幾種禾本科更全面。

表4 慈竹Unigene與其他物種的比較分析Table 4 Analysis on Unigene between Bambusa emeiensis and other species

根據(jù)表4的結(jié)果，再將慈竹的Unigene與水稻轉(zhuǎn)錄本序列和毛竹編碼序列數(shù)據(jù)庫進行Blast比對，選擇E-value＜1e-10。63 094個慈竹Unigene中34 127個 (54.09%) Unigene與毛竹轉(zhuǎn)錄組相匹配，33 685個 (53.39%) Unigene與水稻相匹配。在這些配對的序列中29 591個(46.90%) Unigene與毛竹和水稻都相匹配，4 536個 (7.19%) Unigene只與水稻相匹配，4 094個(6.49%) Unigene只與毛竹相匹配。還有許多預(yù)測的Unigene與水稻和毛竹都不匹配，它們編碼的蛋白質(zhì)主要與生理過程、結(jié)合活性和催化功能有關(guān)，其中有大量預(yù)測的Unigene可能是慈竹特有的 (圖4)。

圖4 慈竹Unigene與水稻和毛竹的比較分析Fig. 4 Analysis on Unigene between Bambusa emeiensis, Oryza sativa and Phyllostachys edulis.

表5 慈竹纖維素合成的功能基因與水稻、短柄草、玉米、高粱和毛竹的比較Table 5 Bambusa emeiensis functional Unigene of cellulose biosynthesis comparison with Oryza sativa, Brachypodium, Zea mays, Sorghum and Phyllostachys edulis

表6 慈竹木質(zhì)素合成的功能基因與水稻、短柄草、玉米、高粱和毛竹的比較Table 6 Bambusa emeiensis functional Unigene of lignin biosynthesis comparison with Oryza sativa, Brachypodium, Zea mays and Sorghum and Phyllostachys edulis

2.4 纖維素、木質(zhì)素生物合成相關(guān)基因與其他物種比對分析

2.4.1 纖維素生物合成相關(guān)酶基因與其他物種比對分析

將鑒定出的慈竹纖維素生物合成途徑中的關(guān)鍵酶基因與水稻、短柄草、玉米、高粱和毛竹進行比對分析，毛竹的相關(guān)數(shù)據(jù)引用自Peng等[1]。由表5可知，在Swiss-Prot數(shù)據(jù)庫中，慈竹纖維素合酶Cellulose synthase (CesA) 的數(shù)量與其他幾個物種相比差異較大，CesA僅有9條。纖維素合酶相似蛋白Cellulose synthase-like protein (Csl) 則是毛竹的最多有76條，慈竹僅有11條，與毛竹相比差異比較明顯。而短柄草中最少，只在Swiss-Prot數(shù)據(jù)庫中鑒定到7條Csl基因。蔗糖合酶Sucrose synthase (SuSy) 數(shù)量最多的是水稻 (21)，其次是短柄草 (11)，慈竹只在Swiss-Prot數(shù)據(jù)庫中鑒定出7條。

2.4.2 木質(zhì)素生物合成相關(guān)酶基因與其他物種比對分析

經(jīng)Swiss-Prot數(shù)據(jù)庫功能注釋，篩選出188個木質(zhì)素生物合成途徑相關(guān)關(guān)鍵酶基因。與其他禾本科植物相比，慈竹苯丙氨酸解氨酶Phenylalanine ammonia-lyase (PAL)、4-香豆酸輔酶A連接酶4-coumarate-CoA ligase (4CL)、肉桂酰輔酶A還原酶Cinnamoyl-CoA reductase (CCR)、肉桂醇脫氫酶Cinnamoyl alcohol dehydrogenase (CAD) 等7個關(guān)鍵酶基因均略少于水稻 (表6)。在Swiss-Prot數(shù)據(jù)庫中，慈竹、短柄草、玉米和高粱的漆酶Laccase (LAC) 基因是數(shù)量最豐富的，分別有17、32、20和20條，毛竹則未檢測到。水稻的CCR，咖啡酰輔酶A-O-甲基轉(zhuǎn)移酶Caffeoyl-CoA O-methyltransferase (CCoAOMT) 和咖啡酸/5-羥基阿魏酸-O-甲基轉(zhuǎn)移酶Caffeic acid-O-methyltransferase (COMT) 基因數(shù)量較其他物種豐富，分別為45條、10條和10條。在慈竹中數(shù)量較少的是COMT，為2條。在毛竹中PAL最多，為8條，其余則較少，與慈竹相比差異較大。毛竹的相關(guān)數(shù)據(jù)引用自Peng等[1]。

2.4.3 轉(zhuǎn)錄因子相關(guān)基因與其他物種比對分析

轉(zhuǎn)錄因子是由核基因編碼的一類蛋白質(zhì)，能選擇性地激活或抑制某些基因的轉(zhuǎn)錄與表達[15]。從Swiss-Prot數(shù)據(jù)庫中得到的轉(zhuǎn)錄因子如表7所示。慈竹MYB和bHLH轉(zhuǎn)錄因子數(shù)量較水稻、短柄草、玉米和高粱多，且轉(zhuǎn)錄因子整體數(shù)量上與毛竹相比差異較明顯。

2.5 差異表達基因分析

2.5.1 差異表達基因的篩選與統(tǒng)計

為了解4個不同高度慈竹筍的差異表達基因，采用DESeq軟件，設(shè)定FDR＜0.01、Fold Change (FC) ≥2作為篩選差異表達基因的關(guān)鍵指標。將10、50、100和150 cm四個樣品進行兩兩比對后，得到以下結(jié)果，10 vs 50 cm組有1 828條差異表達基因；10 vs 100 cm組有1 542條差異表達基因；10 vs 150 cm組有2 601條差異表達基因；50 vs 100 cm組有1 224條差異表達基因；50 vs 150 cm組有1 745條差異表達基因；100 vs 150 cm組有1 231條差異表達基因。將同一組合在不同的數(shù)據(jù)庫間進行比較，提取并統(tǒng)計篩選得到的差異表達基因的注釋信息(表8)。

表7 慈竹轉(zhuǎn)錄因子相關(guān)基因與水稻、短柄草、玉米、高粱和毛竹的比較Table 7 Transcription factor of Bambusa emeiensis comparison with Oryza sativa, Brachypodium, Zea mays and Sorghum and Phyllostachys edulis

表8 注釋的差異表達基因數(shù)目統(tǒng)計Table 8 Summary of differential expression gene annotated

2.5.2 纖維素生物合成相關(guān)差異表達基因分析

慈竹纖維素含量直接影響著造紙效率。通過酸性洗滌纖維 (ADF) 法分析了4個不同高度的慈竹筍基部的纖維素含量。結(jié)果表明，在筍的早期生長階段，隨著筍的伸長，慈竹筍基部的纖維素含量逐漸增加 (圖5)。通過進一步的轉(zhuǎn)錄組數(shù)據(jù)分析，可以從分子機理上了解慈竹筍在發(fā)育過程中纖維素含量的變化。

圖5 不同高度慈竹筍纖維素含量Fig. 5 The content of cellulose in different height of Bambusa emeiensis.

CesA在纖維素生物合成過程中起著關(guān)鍵的調(diào)節(jié)作用。轉(zhuǎn)錄組數(shù)據(jù)分析表明，CesA基因的RPKM值在10和50 cm的樣品中聚在一簇；而100和150 cm的樣品聚為另一簇，表明CesA基因隨慈竹筍高度的增長呈現(xiàn)不同的表達模式(圖6A，6B)。其中，CesA1、CesA5、CesA6和CesA8同源的Unigene有著較高的相對表達量，但卻隨著筍的伸長，呈現(xiàn)出降低的趨勢(圖6A，6B)。在4個不同生長高度的慈竹筍中，與CesA4、CesA7、CesA9同源的10條Unigene，雖然呈現(xiàn)著較低的表達，卻隨著筍的伸長，呈現(xiàn)出增高的趨勢 (圖6A，6B)。這些Unigene在GO數(shù)據(jù)庫的功能注釋中都涉及多個GO功能條目，包括細胞壁、初生細胞壁生物合成、次生細胞壁生物合成、纖維素合酶活性、纖維素生物合成過程、木質(zhì)部發(fā)育、細胞壁增厚等，而在KEGG注釋中，10條Unigene都參與的KEGG通路是K10999 (纖維素合酶) (表9)。

圖6 不同生長高度的慈竹筍中CesA的相對表達水平 (RPKM值) 和差異表達水平 (Log2值)Fig. 6 The relative expression level (RPKM value) and the differential expression level (Log2 value) of CesA in the shoots of different height from Bambusa emeiensis.

表9 纖維素合酶基因的功能注釋Table 9 Annotation of cellulose synthase genes

2.6 差異表達基因qRT-PCR驗證

根據(jù)纖維素生物合成相關(guān)差異表達基因熱圖分析，篩選出CesA9 (T1_Unigene_BMK.22002)進行qRT-PCR驗證。如圖7所示，在100 cm的慈竹筍中，CesA9的相對表達量大約是10 cm筍的9.5倍，表明其表達水平可能直接影響著筍的纖維素含量。

圖7 CesA9 (T1_Unigene_BMK.22002) 基因qRT-PCR驗證Fig. 7 Quantitative RT-PCR verification of CesA9 (T1_Unigene_BMK.22002) gene.

3 討論與結(jié)論

測序技術(shù)日趨成熟，水稻[16]和小麥[17]等禾本科比較有代表性植物的基因組測序已經(jīng)完成，并且已經(jīng)轉(zhuǎn)向?qū)嶋H應(yīng)用方面的研究。在竹子中毛竹[18]基因組測序也已經(jīng)完成，為竹類植物的分子生物學研究奠定了良好基礎(chǔ)。轉(zhuǎn)錄組測序已經(jīng)不再僅限于擬南芥等模式植物，越來越多非模式植物的轉(zhuǎn)錄組也開始得到研究，其對應(yīng)的生長發(fā)育狀態(tài)與機理也逐漸被人們所了解。本文利用轉(zhuǎn)錄組測序技術(shù)，對處于早期生長的4個不同高度的慈竹筍進行轉(zhuǎn)錄組從頭測序，得到了非冗余的111 137條慈竹筍的Unigene。COG和KEGG功能注釋表明，除了基本的功能外，這些Unigene涉及到蔗糖轉(zhuǎn)運與代謝、次級代謝產(chǎn)物 (如與木質(zhì)素相關(guān)的苯丙烷類的生物合成) 及細胞壁的生物合成等方面(圖2，表3)。

值得一提的是竹類植物MYB類轉(zhuǎn)錄因子數(shù)量明顯多于另外4個物種。有研究指出MYB轉(zhuǎn)錄因子家族能通過調(diào)控參與苯丙烷代謝途徑的相關(guān)基因表達，進而影響木質(zhì)部細胞的形成[19]。此外，MYB類轉(zhuǎn)錄因子也會調(diào)控部分纖維素生物合成相關(guān)的基因表達，如MYB46可直接調(diào)控與次生壁相關(guān)的部分纖維素合酶基因，從而影響纖維素的合成[20]。因此，進一步探討MYB類轉(zhuǎn)錄因子在慈竹中的功能，對慈竹分子遺傳改良具有重要的意義。

慈竹纖維素含量豐富，是較好的造紙原料，其與硬頭黃[21]等竹種在四川已形成了以其為龍頭的竹漿造紙工業(yè)產(chǎn)業(yè)鏈。竹材中的纖維素含量直接影響竹漿造紙的效率。10、50、100和150 cm高度的筍中，纖維素含量具有明顯差異。因此本文進一步深入分析了4個樣品中與纖維素合成相關(guān)的差異表達基因。纖維素的生物合成是一個高度復(fù)雜的生物過程，有多個酶參與其生物合成，其中包括CesA、Csl、SuSy等[22]。CesA是目前研究最多的纖維素生物合成基因，它編碼的纖維素合酶不僅是纖維素生物合成的場所，而且在纖維素生物合成過程中起著關(guān)鍵的調(diào)節(jié)作用。纖維素合酶本身具有多基因現(xiàn)象，在纖維素生物合成過程中，多個CesA基因彼此相關(guān)，協(xié)同作用[23]。在慈竹Unigene數(shù)據(jù)庫中，分別與Swiss-port和Nr數(shù)據(jù)庫比對時，與CesA基因同源的Unigene分別有49條和33條；與Csl基因同源的Unigene分別有36條和27條；與SuSy基因同源的Unigene分別有14條和7條。這些基因序列與毛竹、水稻、玉米等物種的CesA、Csl和SuSy基因具有很高的同源性。進一步分析表明，不同高度慈竹筍的纖維素合酶CesA4、CesA7和CesA9基因存在差異表達且隨著筍的伸長，總體呈上升的趨勢 (圖6)。有研究表明，竹類植物筍期的高速生長與筍的居間分生組織中細胞的快速生長與伸長有關(guān)，初生壁的發(fā)育確保了細胞的正常生長，而次生壁的發(fā)育所起到的支撐作用也確保了筍的快速生長[24-25]。另有研究表明，不同的CesA基因分別在細胞初生壁和次生壁合成時期起不同的作用[23]。水稻的OsCesA4、OsCesA7和OsCesA9是細胞次生壁合成所必需的，在木質(zhì)部的細胞中高度表達。OsCesA4、OsCesA7或OsCesA9突變后會導致水稻莖桿變脆且纖維素含量降低，進而影響水稻的正常生長[26]。因此推測慈竹的纖維素合酶基因CesA4、CesA7和CesA9與細胞次生壁的合成有關(guān)且可能影響慈竹筍的生長發(fā)育，這有待進一步研究。

REFERENCES

[1] Peng ZH, Lu Y, Li LB, et al. The draft genome of the fast-growing non-timber forest species moso bamboo (Phyllostachys heterocycla). Nat Genet, 2013, 45(4): 456-463.

[2] Peng ZH, Zhang CL, Zhang Y, et al. Transcriptome sequencing and analysis of the fast growing shoots of moso bamboo (Phyllostachys edulis). PLoS ONE, 2013, 8(11): e78944.

[3] Liu MY, Qiao GR, Jiang J, et al. Transcriptome sequencing and De novo analysis for ma bamboo (Dendrocalamus latiflorus Munro) using the illumina platform. PLoS ONE, 2012, 7(10): e46766.

[4] Qi JQ, Chi B, Xie JL, et al. Study on variations of fiber morphology and tissue proportion of Neosinocalamus affinis culm. Trans China Pulp Paper, 2013, 28(3): 1-4 (in Chinese).齊錦秋, 池冰, 謝九龍, 等. 慈竹纖維形態(tài)及組織比量的研究. 中國造紙學報, 2013, 28(3): 1-4.

[5] Liu Q, He H, Shen ZP. Study on the growth status and correlative environment factors of bamboo, Neosinocalamus affinis from Chengdu region. Sichuan Environ, 2001, 20(4): 43-46 (in Chinese).劉慶, 何海, 沈昭萍. 成都地區(qū)慈竹生長狀況及其與環(huán)境因子關(guān)系的初步分析. 四川環(huán)境, 2001, 20(4): 43-46.

[6] Hu SL, Cao Y, Huang SX, et al. Cloning and bioinformation analysis of 4CL gene in Neosinocalamus affinis. J Northwest A&F Univ: Nat Sci Ed, 2009, 37(8): 204-210 (in Chinese).胡尚連, 曹穎, 黃勝雄, 等. 慈竹4CL基因的克隆及其生物信息學分析. 西北農(nóng)林科技大學學報: 自然科學版, 2009, 37(8): 204-210.

[7] Chen QB, Jiang Y, Lu XQ, et al. Assessment of genetic diversity in Neosinocalamus affinis from different regions of Sichuan province. J Northwest A&F Univ: Nat Sci Ed, 2009, 37(6): 187-193 (in Chinese).陳其兵, 蔣瑤, 盧學琴, 等. 四川不同地區(qū)慈竹的遺傳多樣性研究. 西北農(nóng)林科技大學學報: 自然科學版, 2009, 37(6): 187-193.

[8] Hu SL, Jiang Y, Chen QB, et al. Physical and chemical properties of 2 species in bamboos from the different regions in Sichuan province. Bull Botan Res, 2010, 30(6): 708-712 (in Chinese).胡尚連, 蔣瑤, 陳其兵, 等. 四川2種叢生竹理化特性及纖維形態(tài)研究. 植物研究, 2010, 30(6): 708-712.

[9] Li XL, Lin RS, Fung HL, et al. Chromosome numbers of some caespitose bamboos native in or introduced to China. Acta Phytotaxonom Sin, 2001, 39(5): 433-442 (in Chinese).李秀蘭, 林汝順, 馮學琳, 等. 中國部分叢生竹類染色體數(shù)目報道. 植物分類學報, 2001, 39(5): 433-442.

[10] Tang WW, Wu XL, Zhou J. Comparison of neutral detergent fiber, acid detergent fiber and acid detergent lignin contents of Dongxiang wild rice and cultivated rice. Agric Sci Technol, 2010, 11(3): 11-14.

[11] Grabherr MG, Haas BJ, Yassour M, et al. Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat Biotechnol, 2011, 29(7): 644-652.

[12] Tatusov RL, Galperin MY, Natale DA, et al. The COG database: a tool for genome-scale analysis of protein functions and evolution. Nucleic Acids Res, 2000, 28(1): 33-36.

[13] Ashburner M, Ball CA, Blake JA, et al. Gene Ontology: tool for the unification of biology. Nat Genet, 2000, 25(1): 25-29.

[14] Kanehisa M, Goto S, Kawashima S, et al. The KEGG resource for deciphering the genome. Nucleic Acids Res, 2004, 32(Suppl 1): D277-D280.

[15] Zhang CY, Long Y, Feng J, et al. Transcriptional regulation of plant genes and it’s significance in biology. Hereditas (Beijing), 2007, 29(7): 793-799 (in Chinese).張椿雨, 龍艷, 馮吉, 等. 植物基因在轉(zhuǎn)錄水平上的調(diào)控及其生物學意義. 遺傳, 2007, 29(7): 793-799.

[16] Lin CY, Trinh NN, Lin CW, et al. Transcriptome analysis of phytohormone, transporters and signaling pathways in response to vanadium stress in rice roots. Plant Physiol Biochem, 2013, 66: 98-104.

[17] Zhang J, Liu W, Han H, et al. De novo transcriptome sequencing of Agropyron cristatum to identify available gene resources for the enhancement of wheat. Genomics, 2015, 106(2): 129-136.

[18] Zhao GZ, Sun HY, Zhao HS, et al. A study of genome sequencing of Phyllostachys edulis and its data applications. World Bamb Ratt, 2015, 13(3): 8-13 (in Chinese).趙廣枝, 孫化雨, 趙韓生, 等. 毛竹基因組測序及數(shù)據(jù)應(yīng)用研究現(xiàn)狀. 世界竹藤通訊, 2015, 13(3): 8-13.

[19] Zhao CS, Craig JC, Petzold HE, et al. The xylem and phloem transcriptomes from secondary tissues of the Arabidopsis root-hypocotyl. Plant Physiol, 2005, 138(2): 803-818.

[20] Kim WC, Ko JH, Kim JY, et al. MYB46 directly regulates the gene expression of secondary wall associated cellulose synthases in Arabidopsis. Plant J, 2013, 73(1): 26-36.

[21] Chen QB, Gao SP, Liu L. The selection of paper-pulp bamboo species and the development of bamboo paper sector in Sichuan. J Bamb Res, 2002, 21(4): 47-51 (inChinese).陳其兵, 高素萍, 劉麗. 四川省優(yōu)良紙漿竹種選擇與竹紙產(chǎn)業(yè)化發(fā)展. 竹子研究匯刊, 2002, 21(4): 47-51.

[22] Cheng X, Hao HQ, Peng L. Recent progresses on cellulose synthesis in cell wall of plants. J Trop Subtr Bot, 2011, 19(3): 283-290 (in Chinese).程曦, 郝懷慶, 彭勵. 植物細胞壁中纖維素合成的研究進展. 熱帶亞熱帶植物學報, 2011, 19(3): 283-290.

[23] McFarlane HE, D?ring A, Persson S. The cell biology of cellulose synthesis. Annu Rev Plant Biol, 2014, 65: 69-94.

[24] Gan XH, Ding YL. Investigation on the variation of fiber wall in Phyllostachys edulis culms. For Res, 2006, 19(4): 457-462 (in Chinese).甘小洪, 丁雨龍. 毛竹莖稈纖維發(fā)育過程中細胞壁的變化規(guī)律研究. 林業(yè)科學研究, 2006, 19(4): 457-462.

[25] Cui K, He CY, Zhang JG, et al. Characteristics of temporal and spatial tissue development during the rapidly growing stage of moso bamboo culms. For Res, 2012, 25(4): 425-431 (in Chinese).崔凱, 何彩云, 張建國, 等. 毛竹莖稈組織速生的時空發(fā)育特征. 林業(yè)科學研究, 2012, 25(4): 425-431.

[26] Tanaka K, Murata K, Yamazaki M, et al. Three distinct rice cellulose synthase catalytic subunit genes required for cellulose synthesis in the secondary wall. Plant Physiol, 2003, 133(1): 73-83.

(本文責編 陳宏宇)

March 9, 2016; Accepted: July 20, 2016

Shanglian Hu. Tel/Fax: +86-138-81194095; E-mail: hushanglian@126.com

Transcriptome analysis and gene function annotation of Bambusa emeiensis shoots based on high-throughput
sequencing technology

Yupeng Chen, Ying Cao, Shanglian Hu, Yan Huang, Xueqin Lu, Gang Xu, and Zhijian Long
Engineering Research Center for Biomass Resource Utilization and Modification of Sichuan Province, School of Life Science and Engineering, Southwest University of Science and Technology, Mianyang 621010, Sichuan, China

Bambusa emeiensis is one of the preponderant species of sympodial bamboos in Sichuan province of China, and has excellent fiber length and quality as raw materials for papermaking, textile and other industries. In this study, with the application of Illumina HiSeqTM2000 platform, we analyzed transcriptome in B. emeiensis with different heights of 10, 50, 100 and 150 cm. A total of 69.28 M reads were obtained, and a sum up of 111 137 bands of Unigenes were acquired following de novo stitching, assembly and clustering, among which there were 63 094 bands that had been integrated in the COG, GO, KEGG, Swiss-Prot and Nr databases using annotated methods. These Unigenes not only had general functions, such as transcription and signal transduction, but were also involved in sucrose transport and metabolism, secondary metabolites and cell wall biosynthesis. There was significant difference regarding the expression of cellulose synthase gene in B. emeiensis at different heights, relevant genes were found that might be responsible for the regulation of the growth and development of B. emeiensis as well as the biosynthesis of cellulose and lignin. Our findings could provide some elementary theories for breed improvement of B. emeiensis.

Bambusa emeiensis, transcriptome, high-throughput sequencing, function annotation

Supported by: National Natural Science Foundation of China (Nos. 31400257, 31400333), Breeding Program Fund Project by the 13th Five-Year Plan of Sichuan Province, Fund of Engineering Research Center for Biomass Resource Utilizaiton and Modification of Sichuan Province (Nos. 12zxsk07, 13zxsk01, 14tdgc05), Major Project of Education Department in Sichuan Province (No. 16ZA0145), Postgraduate Innovation Fund Project by Southwest University of Science and Technology (No. 15ycx089).

國家自然科學基金 (Nos. 31400257, 31400333)，四川省“十三五”育種公關(guān)項目，四川省生物質(zhì)資源利用與改性工程技術(shù)研究中心基金 (Nos. 12zxsk07, 13zxsk01, 14tdgc05)，四川省教育廳科研項目 (No. 16ZA0145)，西南科技大學研究生創(chuàng)新基金 (No. 15ycx089) 資助。