基于转录组测序的慈竹笋木质素基因筛选及其表达分析

蒋海港 王树刚 刘 阳 沈彦合 李祥宁 曹 颖

(西南科技大学生命科学与工程学院，绵阳 621010)

木质素(lignin)是一种复杂的酚类高分子化合物，其生物合成主要是通过苯丙酮酸途径(phenylpropanoid pathway)和木质素特异合成途径[1]。苯丙酮酸途径上的苯丙氨酸氨裂解酶(Phenylalanine ammonia-lyase,PAL)、肉桂酸4-羟基化酶(Cinnamate 4-hydroxylase,C4H)、香豆酸3-羟化酶(C3H)、和4-香豆酸CoA连接酶(4-coumarate-CoA ligase,4CL)，它们表达活性的高低与木质素总量密切相关；而木质素特异合成途径上的咖啡酸O-甲基转移酶(Caffeic acid O-methyltransferase,COMT)、咖啡酰CoA-O-甲基转移酶(Caffeoyl-CoA-O-methyltransferase,CCoAOMT)和阿魏酸-5-羟基化酶(Ferulate 5-hydroxylase,F5H)，与木质素单体特异合成相关，它们的表达对木质素含量尤其木质素单体的特异合成影响较大，决定了各种单体在木质素总量中的比例。羟基肉桂酸脱氢酶(Cinnamoyl alcohol dehydrogenase,CAD)和肉桂酰-CoA还原酶(Cinnamoyl-CoA reductase,CCR)是位于木质素特异合成途径下游的还原酶，它们负责将各种羟基肉桂酰-辅酶A(CoA)酯还原成各种木质素单体，对木质素的合成与代谢也起关键作用[2]。另外，木质素单体的聚合过程还需要过氧化物酶(Peroxidase)和漆酶(Laccase,LAC)的催化，但在木质素聚合过程中的具体功能至今还不清楚[3]。

尽管木质素是植物生长发育所必须的，但由于木质素和纤维素高度胶联，木质素成为生物质能源利用如造纸和生物乙醇生产的主要障碍[4～6]。竹类植物纤维素含量较高，一般在40%～60%，竹浆的性能在针叶林和阔叶林之间，优于草浆类，是优良的造纸用纤维原料[7]。但是，竹类木质素含量比其他草本植物要高，研究竹类木质素生物合成途径，寻找其关键基因，遗传调控竹类植物的木质素合成，培育高纤维素、低木质素的竹材，对造纸行业及生物乙醇生产具有重大意义。近年来竹类植物的分子生物学研究取得了一定的进展。毛竹(Moso bamboo,Phyllostachysheterocycla)基因组草图的发布[8]，以及毛竹、麻竹(Ma bamboo,DendrocalamuslatiflorusMunro)、巨龙竹(Da bamboo,Dendrocalamussinicus)等竹种的转录组和(或)蛋白组研究，为从分子水平揭示竹类植物高速生长和开花机制，以及细胞壁发育与调控机制奠定了基础[9～12]。慈竹(Bambusaemeiensis)，是四川本土大型丛生竹之一，其秆径较小且秆壁薄，纤维含量和长度优良，是较好的造纸原料[13～14]。但与毛竹、麻竹等竹种相比，慈竹分子生物学研究相对滞后，其木质素生物合成途径和相关功能基因仍不清楚。

竹子为单子叶植物,没有次生生长。竹笋伸长主要是由节部居间分生组织的活动所决定的；笋高度的增加是节间伸长的结果。随着笋的伸长，笋各部分纤维素和木质素含量逐渐积累[10]。竹笋的这种独特的顺序发育机制，使其成为研究竹子高生长特性，以及纤维素和木质素生物合成机制的理想模型。近年来，竹类植物的基因组、转录组及蛋白质组等研究也多以竹笋为对象进行的[10～12]。我们前期利用Illumina HiSeqTM2000平台，对露出地面10、50、100和150 cm高的慈竹笋进行了转录组测序，获得了早期伸长阶段的慈竹笋转录组数据库，分析了其中纤维素合成关键基因[15]，但尚未对木质素合成途径和相关基因表达进行深入分析。

本文根据慈竹笋RNA高通量测序结果，用拟南芥(Arabidopsisthaliana)、水稻(Oryzasativa)等植物的已知木质素生物合成基因作为查询序列，通过BLASTp分析，结合系统进化分析，在笋转录组数据库中比对鉴定出了一系列与木质素生物合成相关的基因。同时，对这些木质素合成基因在不同发育阶段笋中的差异表达进行了分析，探讨了慈竹笋早期伸长过程中木质素含量变化及其与分子表达的关系，预测到一些可能参与发育性木质素积累的候选基因，为慈竹的分子遗传育种以及品种改良等提供一定的理论基础。

1 材料和方法

1.1 植物材料

慈竹材料采自西南科技大学(四川绵阳)资源圃(年平均气温17.2℃，年平均降雨量793.5 mm)。分别取立地条件一致的当年生的10，50，100和150 cm高度的笋，每个高度取3株独立的笋(测量误差小于2 cm)，取基部第二个节间组织切碎、混匀，液氮速冻后，储存于-80℃冰箱中备用。

1.2 组织解剖和木质素含量分析

笋组织样品(0.5 cm3)取自不同高度慈竹笋基部第二个节间中部，经FAA溶液固定后保存，常规的石蜡切片法进行组织切片，番红固绿染色[16]。Leica LDM 2500显微照相，用内置的标尺进行显微测量，数字图像用用Image ProPlus software(IPP6.0,Media Cybernetics,Inc.)分析。木质素含量的测定采用乙酰溴法[17]，每个样品取3次生物学重复，并重复测定3次。

1.3 转录组数据分析

采用Plant RNA Kit(OMEGA BIO-TEK)分别提取露出地面10、50、100和150 cm高的慈竹笋总RNA，构建4个样品的cDNA测序文库。利用Illumina HiSeqTM2000平台进行转录组测序。测序得到的原始数据经过去除杂质和去冗余处理，以及序列拼接和聚类分析后，组装获得Unigene序列。以拟南芥、水稻等已知木质素基因为查询序列，通过BLASTp在慈竹笋Unigene数据库中筛选假定的木质素合成相关序列。根据其与拟南芥、水稻同源基因注释命名基因。将获得的序列经DAMBE检测饱和度[18]后导入MEGA7.0，用最大似然法(Maximum Likehood,ML)，选择WAG+G+F模型构建系统进化树。对生成的系统进化树进行Bootstrap校正(参数设为500 replicates)，生成最终的系统进化树。

1.4 基因差异表达和定量qRT-PCR分析

基于RNA高通量测序结果，统计考察基因的相对表达丰度(reads per kilo bases per million mapped reads,RPKM)。同时，将各样品的RPKM值经均一化处理后，以10cm笋的RPKM值为参照，以|log2 FC(Fold Change)|≥2(FDR<0.01)为标准，并运用MeV软件做层次聚类分析，筛选差异表达基因。选择一些差异表达基因，用SyBR Green PCR Master Mix(天根生物)和基因特异引物(表2)在iQ5荧光定量PCR系统(BIO-RAD)进行qRT-PCR分析。每个样品被扩增3次。以慈竹Tublin基因作为内参，采用2-ΔΔCt法[19]计算基因相对表达量。

2 结果与分析

2.1 早期伸长阶段慈竹笋的木质素含量和组织解剖结构

慈竹在发笋初期出现木质化现象。早期生长阶段慈竹笋(10～150 cm高度)木质素含量分析表明，随着笋伸长，慈竹笋基部的木质素含量逐渐增加(图1A)。解剖结构分析表明，幼嫩的10 cm笋基部节间出现可辨识的导管组织，维管束结构基本形成。随着笋的伸长，维管束中的纤维细胞的数量随笋龄的增加而增加，发生木质化。50 cm笋中初步出现了木质化的纤维细胞，随着笋高度的增加，木质化程度加深(图1B)。这与木质素含量分析结果相一致。

图1 早期伸长阶段慈竹笋木质素含量和解剖结构 A.不同高度笋木质素含量；B.不同高度笋基部节间解剖分析 F.纤维细胞；LF.纤维素细胞木质化；MX.后生木质部；P.薄壁组织；PP.原生韧皮部；PX.原生木质部Fig.1 Lignin content and anatomical analysis of B.emeiensis shoots during the early elongation stages A. Lignin content in developing bamboo shoots; B. Paraffin transverse and longitudinal sections of culms at the basal 2nd internode of different height shoots F. Fiber cells; LF. Lignification of fiber cells; MX. Metaxylem; P. Parenchyma cells; PP. Protophloem; PX. Protoxylem

2.2 慈竹笋中木质素生物合成相关基因

BLASTp分析从慈竹笋转录组数据库比对鉴定出351个与木质素生物合成相关Unigenes(E-value<10-10)：其中与木质素单体聚合相关的Peroxidase和LAC的转录子最为丰富，分别为155条和51条(表1)；与木质素单体合成相关Unigenes，135个，其中尤以编码4CL和CCR的序列最为丰富，分别为37条和34条，其次是PAL和CAD转录子，分别为26和25条，而参与甲基化过程的COMT和CCoAOMT和参与羟基化C3H、C4H和F5H蛋白，仅有少数转录子被检测到(表1)。这与已报道的麻竹、毛竹、巨龙竹等其他的竹类植物是相似的[9,11～12]。

为避免注释的错误，我们从BLASTp筛选到的135个木质素单体合成相关Unigenes(序列)和51条LAC序列中，选取了59条可信度较高的Unigene序列(E value=0)，进行了系统进化分析。系统进化分析结果证明了BLAST比对和注释结果的准确性。在59条Unigene序列中，仅有4个序列(2个CCR、1个LAC和1个F5H)没有聚类到与注释结果对应的进化簇分支中，这可能与这4个Unigene序列长度不足有关(<1 000 bp)(图2)。

2.3 木质素单体合成相关基因的差异表达

转录组数据分析表明，在4个不同高度慈竹笋中，控制木质素生物合成碳流入口的苯丙酮酸途径中的一些关键酶基因如PAL和4CL，以及一些参与特异木质素单体合成的CCR和CAD基因呈现较高的相对表达丰度，其相对表达水平(RPKM值)在162～698范围内(图3A)。在50、100和150 cm笋中，大部分基因的表达水平较为接近，但与较幼嫩10cm笋距离较远(图3A)。

基因的差异表达分析表明，随着笋的伸长，木质素合成途径中4种关键酶的编码基因，如PAL、4CL、CCR、CAD呈现明显的上调表达。其中2条PAL1(CL56Contig1和CL56Contig4)，2条4CL3(CL1213Contig2和T3Unigene34303)和1个4CL4(CL1213Contig1)、5条CCR1(T3Unigene19278、T3Unigene31394、CL21423Contig1、T4Unigene31965和T1Unigene17952)和2条CADH2(CL25770Contig、CL17216Contig1)基因，其表达明显上调(│log2 FC(Fold Change)│≥2，图3B)，表明这些基因的表达可能与发育性木质化进程相关。

表1慈竹、巨龙竹、麻竹和毛竹等竹类植物的木质素合成基因数目

Table1NumberofligningenesfoundinB.emeiensisandotherbamboospeciesincludingDabamboo,MabambooandMosobamboo

酶Enzymes慈竹Bambusa emeiensisa巨龙竹Da bamboo[12]a麻竹Ma bamboo[9]a毛竹Moso bamboobPAL26151713C4H28411C3H36034CL3726356HCT122704CCR3447103CCoAOMT3942CAD251521F5H151313COMT2272DAHPS10672LAC51273416Peroxidase155174163143A1dOMT0002

注：木质素合成相关基因的数据分别来自于转录组数据库a和基因组数据库(Bamboo Genome Database，http://www.bamboogdb.org/)b

Note：Number of lignin biosynthesis genes found in the bamboo transcriptomeaand bamboo genome( http://www.bamboogdb.org/)b,respectively.

图2 慈竹木质素生物合成基因系统进化分析Fig.2 The phylogenetic analysis of lignin biosynthesis genes identified in shoots of B.emeiensis

图3 不同生长高度的慈竹笋木质素基因的相对表达(RPKM值)和差异表达(log2值)Fig.3 The relative expression(RPKM value,A) and the differential expression(log2 value,B) of lignin biosynthesis genes in the shoots of B.emeiensis with different height

2.4 参与木质素单体聚合的漆酶(LAC)

木质素单体在细胞质中形成后，需要运输到细胞壁进行聚合反应生成木质素。LAC被报道参与植物木质素的聚合反应[3]。在慈竹笋转录组数据库中，有51条假定编码LAC的序列被发现，经Swiss-Prot数据库功能注释，分别与17个水稻LAC基因如OsLAC3、5、22、25等同源。由图4知，在4个不同高度慈竹笋中，大部分LAC基因随着笋的伸长而表达上调，上调表达程度比下调表达程度高。差异表达分析显示10 vs 100和10 vs 150聚为一簇，表明100和150 cm生长高度慈竹笋中LAC同源基因表达情况相近，而与50 cm样品差异较大。随着笋伸长，LAC5、LAC10、LAC11和LAC22同源基因的表达明显上调；而LAC17、LAC19、LAC24、LAC25呈现下调趋势。将这些差异表达基因与已知水稻LAC同源基因进行了聚类分析，结果表明上调的LAC基因与下调的LAC基因分别聚在a、b两个亚簇上，呈现明显的系统进化差异(图4)。

2.5 差异表达基因的qRT-PCR分析

根据木质素生物合成相关差异表达基因热图分析，选择4个差异表达基因LAC5(CL33Contig7)、CADH2(CL17216Contig1)、4CL3(T3Unigene34303)和PAL1(CL56Contig4)，进行qRT-PCR验证。

如图5所示，在基部第二个节间，LAC5、CADH2、4CL3和PAL的表达均随着笋的伸长而明显增加。尤其170 cm笋中这4个基因的相对表达量是30cm和50cm笋的数十到数百倍，表明这4个差异表达基因的表达水平可能直接影响着笋的木质素含量。

表2 RT-PCR引物

图4 不同生长高度的慈竹笋中LAC的差异表达热图(log2值，A)和系统进化分析(B)Fig.4 The heat map of(log2 value,B) of LAC differential expressions in the shoots of B.emeiensis with different height and their phylogenetic analysis

图5 不同高度笋中部分木质素基因的表达Fig.5 The expression of lignin genes in shoots with different height(P<0.01)

3 讨论

对竹秆的解剖研究表明，竹子维管组织形成较早，未伸出地面的毛竹笋的基部节间已经有可见的木质部导管和基本的维管束组织出现[20]。毛竹笋的纤维组织的木质化，从竹秆6 m左右时开始[10]。与毛竹相比，慈竹笋纤维细胞的木质化发生较早。在50 cm高的慈竹笋基部就出现了木质化的纤维细胞，随着笋高度的增加，木质化程度加深(图1B)。对慈竹早期发育过程中木质素生物合成基因的表达水平分析也表明，大部分木质素基因的表达水平在50、100和150 cm笋中较为接近，而与较幼嫩10 cm笋距离较远。因此早期伸长阶段的慈竹笋是研究生木质素合成和积累的理想模型。

竹类植物含有较多的木质素含量，这可能是由木质素合成关键酶的数量和表达水平上的差异引起的。在所有已报道的巨龙竹、麻竹和毛竹等竹类植物转录组中，除了Peroxidase外，LAC和4CL转录子是最显著的，另外也有较多数量的假定编码PAL、CAD和CCR的序列被检测到。在慈竹中，我们检测到了37个4CL，26个PAL和25个CAD转录子，其数量均明显高于其他已报道的竹类植物(表1)。另外，慈竹笋中可检测到34个CCR转录子，虽然少于巨龙竹(47个)，但明显高于麻竹(10个)和毛竹(19个)。另外，在慈竹笋中没有检测到5-羟基松柏醛O-甲基转移酶相关(5-hydroxyconiferyl aldehyde O-methyltransferase,AldOMT)的转录子，这与在麻竹、毛竹、巨龙竹中的研究结果一致(表1)，可能是由于竹类植物中其它表达的甲基转移酶取代了AldOMT的活性，这些特征可能是竹类植物独有的，其中内在的调控机制需要进一步的研究加以阐明。

木质素主要沉积在次生壁上，对于细胞乃至整个植株都起到支撑的作用，与慈竹笋的快速生长关系紧密。转录组数据和RT-PCR分析表明，在4个不同高度慈竹笋中，控制木质素生物合成碳流入口的苯丙酮酸途径中的一些关键酶基因如PAL和4CL，以及一些参与特异木质素单体合成的CCR和CAD基因呈现较高的相对表达丰度，其RPKM值在162～698。PAL、4CL、和CAD基因的上调表达，表明在笋的早期生长过程中木质素就已经开始合成，该现象与笋的快速生长需要机械强度增加[21]的结果相一致。LAC(漆酶)介导的木质素单体的聚合是木质化过程中必要的一个步骤。而且，在植株发育的过程中漆酶能在一定程度上调控维管束的木质化[22]，LAC为调控目标的木质素工程是改良能源植物的一个新的可行途径[3]。与其他已报道的竹类植物类似，慈竹笋中检测到51条假定编码LAC的序列。随着笋的伸长，a亚簇上LAC同源的基因表达上调，而b亚簇上LAC同源的基因呈现下调趋势(图5)。这些结果表明，这些木质素合成基因的表达存在时空差异，木质素合成和积累可能存在多基因的协调作用。

1.Boerjan W,Ralph J,Baucher M.Lignin biosynthesis[J].Annual Review of Plant Biology,2003,54:519-546.

2.Carocha V,Soler M,Hefer C,et al.Genome-wide analysis of the lignin toolbox ofEucalyptusgrandis[J].New Phytologist,2015,206(4):1297-1313,doi:10.1111/nph.13313.

3.Wang J H,Feng J J,Jia W T,et al.Lignin engineering through laccase modification:a promising field for energy plant improvement[J].Biotechnology for Biofuels,2015,8:145.

4.Himmel M E,Ding S Y,Johnson D K,et al.Biomass recalcitrance:engineering plants and enzymes for biofuels production[J].Science,2007,315(5813):804-807.

5.Simmons B A,Loqué D,Ralph J.Advances in modifying lignin for enhanced biofuel production[J].Current Opinion in Plant Biology,2010,13(3):313-319.

6.Harrington M J,Mutwil M,Barrière Y,et al.Molecular biology of lignification in grasses[J].Advances in Botanical Research,2012,61:77-112.

7.陈其兵,高素萍,刘丽.四川省优良纸浆竹种选择与竹纸产业化发展[J].竹子研究汇刊,2002,21(4):47-51.

Chen Q B,Gao S P,Liu L.the selection of paper-pulp bamboo species and the development of bamboo paper sector in Sichuan[J].Journal of Bamboo Research,2002,21(4):47-51.

8.Peng Z H,Lu Y,Li L B,et al.The draft genome of the fast-growing non-timber forest species Moso bamboo(Phyllostachysheterocycla)[J].Nature Genetics,2013,45(4):456-461.

9.Liu M Y,Qiao G R,Jiang J,et al.Transcriptome sequencing andDenovoanalysisfor ma bamboo(DendrocalamuslatiflorusMunro) using the illumina platform[J].PLoS One,2012,7(10):e46766.

10.Cui K,He C Y,Zhang J G,et al.Temporal and spatial profiling of Internode elongation-associated protein expression in rapidly growing culms of bamboo[J].Journal of Proteome Research,2012,11(4):2492-2507.

11.Peng Z H,Zhang C L,Zhang Y,et al.Transcriptome sequencing and analysis of the fast growing shoots of Moso bamboo(Phyllostachysedulis)[J].PLoS One,2013,8(11):e78944.

12.Cui K,Wang H Y,Liao S X,et al.Transcriptome sequencing and analysis for culm elongation of the world’s largest bamboo(Dendrocalamussinicus)[J].PLoS One,2016,11(6):e0157362.

13.胡尚连,蒋瑶,陈其兵,等.四川2种丛生竹理化特性及纤维形态研究[J].植物研究,2010,30(6):708-712.

Hu S L,Jiang Y,Chen Q B,CAO Y,SUN X ,LU X Q.Physical and chemical properties of 2 species in Bamboos from the different regions in Sichuan Province[J].Bulletin of Botanical Research,2010,30(6):708-712

14.齐锦秋,池冰,谢九龙,等.慈竹纤维形态及组织比量的研究[J].中国造纸学报,2013,28(3):1-4.

Qi J Q,Chi B,Xie J L,et al.Study on Variations of Fiber Morphology and Tissue Proportion ofNeosinocalamusaffinisCulm[J].Transactions of China Pulp and Paper,2013,28(3):1-4.

15.陈宇鹏,曹颖,胡尚连,等.基于高通量测序的慈竹笋转录组分析与基因功能注释[J].生物工程学报,2016,32(11):1610-1623.

Chen Y P,Cao Y,Hu S L,et al.Transcriptome analysis and gene function annotation ofBambusaemeiensisshoots based on high-throughput sequencing technology[J].Chinese Journal of Biotechnology,2016,32(11):1610-1623.

16.崔凯,何彩云,张建国,等.毛竹茎秆组织速生的时空发育特征[J].林业科学研究,2012,25(4):425-431.

Cui K,He C Y,Zhang J G,et al.Characteristics of temporal and spatial tissue development during the rapidly growing stage of moso bamboo culms[J].Forest Research,2012,25(4):425-431.

17.李靖,程舟,杨晓伶,等.紫外分光光度法测定微量人参木质素的含量[J].中药材,2006,29(3):239-241.

Li J,Cheng Z,Yang X L,et al.Determination of lignin content in TinyPanaxginsengby UV spectrophotometry[J].Journal of Chinese Medicinal Materials,2006,29(3):239-241.

18.Xia X H,Xie Z,Salemi M,et al.An index of substitution saturation and its application[J].Molecular Phylogenetics and Evolution,2003,26(1):1-7.

19.Livak K J,Schmittgen T D.Analysis of relative gene expression data using real-time quantitative PCR and the 2-ΔΔCTmethod[J].Methods,2001,25(4):402-408.

20.Wei Q,Jiao C,Guo L,et al.Exploring key cellular processes and candidate genes regulating the primary thickening growth of Moso underground shoots[J].New Phytologist,2017,214(1):81-96,doi:10.1111/nph.14284.

21.徐有明,郝培应,费本华.竹笋形态发育构建过程中解剖结构及化学成分的动态变化[J].东北林业大学学报,2008,36(4):8-11,16.

Xu Y M,Hao P Y,Fei B H.Dynamic change of anatomical structure and chemical constitutes of bamboo shoots forPhyllostachyspubescenceat developmental stage[J].Journal of Northeast Forestry University,2008,36(4):8-11,16.

22.Zhao Q,Nakashima J,Chen F,et al.LACCASEis necessary and nonredundant withPEROXIDASEfor lignin polymerization during vascular development inArabidopsis[J].The Plant Cell,2013,25(10):3976-3987.