陳宇鵬，曹穎，胡尚連，黃艷，盧學(xué)琴，徐剛，龍治堅(jiān)

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

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

陳宇鵬，曹穎，胡尚連，黃艷，盧學(xué)琴，徐剛，龍治堅(jiān)

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

陳宇鵬, 曹穎, 胡尚連, 等. 基于高通量測(cè)序的慈竹筍轉(zhuǎn)錄組分析與基因功能注釋. 生物工程學(xué)報(bào), 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.

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

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

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

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

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

1 材料與方法

1.1 材料

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

1.2 不同高度慈竹筍纖維素含量測(cè)定

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

1.3 轉(zhuǎn)錄組測(cè)序

1.3.1 總RNA瓊脂糖凝膠電泳檢測(cè)及測(cè)序

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

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

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

1.3.3 功能注釋

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

1.3.4 基因的差異表達(dá)分析

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

1.4 差異表達(dá)基因qRT-PCR驗(yàn)證

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

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

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

2 結(jié)果與分析

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

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

通過(guò)進(jìn)一步分析后得到的非冗余Unigene庫(kù)，含有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 基因功能注釋

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

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

圖2 Unigene的COG功能分類(lèi)Fig. 2 Unigene of COG function categories.

2.2.2 功能基因的COG注釋

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

2.2.3 功能基因的GO分類(lèi)

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

2.2.4 功能基因的KEGG注釋

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

圖3 Unigene的GO功能分類(lèi)Fig. 3 Unigene of GO function categories.

表3 Unigene數(shù)目最多的5個(gè)代謝通路及與纖維素、木質(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進(jìn)行功能比對(duì)。所得結(jié)果如表4所示，慈竹Unigene與毛竹匹配度最高，達(dá)到52.46%，可能是因?yàn)榇戎衽c毛竹同屬竹類(lèi)植物的原因。其次是與水稻相匹配達(dá)到50.41%，可能是因?yàn)樗緸槟Ｊ街参?，其基因序列信息較其他幾種禾本科更全面。

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

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

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

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

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

經(jīng)Swiss-Prot數(shù)據(jù)庫(kù)功能注釋?zhuān)Y選出188個(gè)木質(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個(gè)關(guān)鍵酶基因均略少于水稻 (表6)。在Swiss-Prot數(shù)據(jù)庫(kù)中，慈竹、短柄草、玉米和高粱的漆酶Laccase (LAC) 基因是數(shù)量最豐富的，分別有17、32、20和20條，毛竹則未檢測(cè)到。水稻的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)基因與其他物種比對(duì)分析

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

2.5 差異表達(dá)基因分析

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

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

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

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

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

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

圖6 不同生長(zhǎng)高度的慈竹筍中CesA的相對(duì)表達(dá)水平 (RPKM值) 和差異表達(dá)水平 (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 差異表達(dá)基因qRT-PCR驗(yàn)證

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

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

3 討論與結(jié)論

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

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

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

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)及組織比量的研究. 中國(guó)造紙學(xué)報(bào), 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ū)慈竹生長(zhǎng)狀況及其與環(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基因的克隆及其生物信息學(xué)分析. 西北農(nóng)林科技大學(xué)學(xué)報(bào): 自然科學(xué)版, 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).陳其兵, 蔣瑤, 盧學(xué)琴, 等. 四川不同地區(qū)慈竹的遺傳多樣性研究. 西北農(nóng)林科技大學(xué)學(xué)報(bào): 自然科學(xué)版, 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).李秀蘭, 林汝順, 馮學(xué)琳, 等. 中國(guó)部分叢生竹類(lèi)染色體數(shù)目報(bào)道. 植物分類(lèi)學(xué)報(bào), 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)控及其生物學(xué)意義. 遺傳, 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).趙廣枝, 孫化雨, 趙韓生, 等. 毛竹基因組測(cè)序及數(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).程曦, 郝懷慶, 彭勵(lì). 植物細(xì)胞壁中纖維素合成的研究進(jìn)展. 熱帶亞熱帶植物學(xué)報(bào), 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ò)程中細(xì)胞壁的變化規(guī)律研究. 林業(yè)科學(xué)研究, 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).崔凱, 何彩云, 張建國(guó), 等. 毛竹莖稈組織速生的時(shí)空發(fā)育特征. 林業(yè)科學(xué)研究, 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.

(本文責(zé)編 陳宏宇)

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).

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