, 黄艳, 卢学琴, 徐刚, 龙治坚西南科技大学 生命科学与工程学院 四川省生物质资源利用与改性工程技术研究中心，四川 绵阳 621010
基金项目：国家自然科学基金(Nos. 31400257, 31400333)，四川省“十三五”育种公关项目，四川省生物质资源利用与改性工程技术研究中心基金(Nos. 12zxsk07, 13zxsk01, 14tdgc05)，四川省教育厅科研项目(No. 16ZA0145)，西南科技大学研究生创新基金(No. 15ycx089)资助
摘要: 慈竹是我国四川当地的优势丛生竹种之一，其纤维长度和质量较优异，是造纸、纺织等工业的良好原料。本文利用Illumina HiSeqTM 2000平台，对10、50、100和150 cm高的慈竹笋进行转录组分析，共得到69.28 M条读长(Reads)，经从头拼接、组装和聚类后得到111 137条非重复序列基因Unigene，其中共有63 094条注释到COG、GO、KEGG、Swiss-Prot和Nr数据库中。这些Unigene不仅具有一般的功能，如转录和信号转导等，还涉及到蔗糖转运与代谢、次级代谢产物及细胞壁的生物合成等方面。不同高度慈竹笋的纤维素合成酶基因存在差异表达，发现了可能调控慈竹生长发育以及纤维素和木质素生物合成的相关基因，为慈竹品种改良提供一定的理论基础。
关键词: 慈竹 转录组 高通量测序 功能注释
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, 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
Received: March 9, 2016; Accepted: July 20, 2016
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)
Corresponding authors:Shanglian Hu. Tel/Fax: +86-138-81194095; E-mail: hushanglian@126.com
Abstract: 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 HiSeqTM 2000 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.
Key words: Bambusa emeiensis transcriptome high-throughput sequencing function annotation
竹子是重要的非木材可再生资源，以生长速度快、产量高、易加工、用途广泛、经济价值高等特点著称。近年来随着分子生物学研究手段的发展，尤其是基因组和转录组测序技术的广泛应用，毛竹、麻竹等竹种的基因组和(或)转录组数据已经被公布，一些与竹笋快速生长以及竹材纤维形成的相关基因被报道。如Peng等[1-2]通过基因组和转录组测序发现了19个毛竹纤维素合酶基因，和大量与植物激素、细胞周期调控、细胞壁合成、细胞形态等相关的候选基因，为毛竹竹材特性的改良提供了依据。Liu等[3]分析了麻竹的转录组，筛选出105个木质素生物合成途径相关关键酶基因，并与毛竹的相关基因进行了比对，确定了潜在的与生长和发育相关的候选基因。
慈竹Bambusa emeiensis，是四川本土大型丛生竹之一，其秆径较小，秆壁较薄，纤维长度和质量均优于其他竹种，是造纸、纺织等行业的优良原料之一，具有比较重要的应用前景[4]。目前，对于慈竹的研究主要集中于生长状况[5]、基因克隆[6]、遗传多样性[7]、生理生化指标[8]等方面。与毛竹、麻竹等竹种相比，慈竹的分子生物学研究相对滞后。有研究表明慈竹可能有70条染色体，为典型的多倍体复合体[9]，其基因组测序的难度较大，而且目前慈竹转录组分析以及纤维素、木质素生物合成相关的分子机理等也未见报道。
本文以不同高度的慈竹笋为材料，比较了不同高度慈竹笋的纤维素含量，并采用Illumina HiSeqTM 2000平台对其进行了转录组测序。将得到的Unigene进行Cluster of Orthologous Groups of proteins (COG)、Gene Ontology (GO)、Kyoto Encyclopedia of Genes and Genomes (KEGG)、Swiss-Prot和Nonredundant protein (Nr)注释，并与其他相似物种进行比对，旨在发现调控慈竹生长发育以及与纤维素、木质素生物合成相关基因，挖掘慈竹中的珍贵基因资源。同时，本文也进一步分析了4个样品中与纤维素合成相关的差异表达基因，为慈竹的分子遗传育种以及品种改良等提供一定的理论基础。
1 材料与方法1.1 材料供试材料采自西南科技大学生命科学与工程学院资源圃(年平均气温17.2 ℃，年平均降雨量793.5 mm)。采样时，以露出地面10、50、100和150 cm的笋为对象，每个高度取3株独立的笋(测量误差小于2 cm)，并取其基部第2个节间的部分，分别切碎、混匀，其中一部分储存于-80 ℃冰箱中备用，另一部分置于烘箱中烘干用于纤维素含量测定。
1.2 不同高度慈竹笋纤维素含量测定纤维素含量的测定参照FOSS公司FibertecTM M6 1020/1021型纤维素测定仪的操作手册进行。通过酸性洗涤纤维(ADF)[10]法，得到纤维素含量，每个样品生物重复3次，每组样品运行1个空白对照。
1.3 转录组测序1.3.1 总RNA琼脂糖凝胶电泳检测及测序RNA提取采用OMEGA BIO-TEK公司Plant RNA Kit试剂盒中的试剂及提取方法。总RNA样品采用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)构建上机文库。制备好的文库用1.8%琼脂糖凝胶电泳进行检测文库插入片段大小，然后用Library Quantification Kit-Illumina GA Universal (Kapa, KK4824)进行QPCR定量，再在cDNA片段两端加上接头。检测合格的文库在illumina cbot上进行簇的生成，最后用Illumina HiSeqTM 2000进行测序。
1.3.2 数据分析转录组测序得到的原始数据经过去除杂质和冗余处理后，利用Trinity软件[11]对经过过滤后的高质量数据进行从头拼接，获得重叠群Contig。根据其结果利用双末端信息将来自同一转录本的不同Contig连在一起，做进一步的序列拼接，得到转录本Transcripts。在Transcripts聚类单元中选取最主要的转录本作为Unigene序列，并对Unigene数据进行聚类分析和进一步去冗余处理，最终得到非冗余Unigene库。对于多样品的组装，由于后续的Open reading frame (ORF)预测、Simple sequence repeats (SSR)分析、表达丰度分析等都建立在同一套参考基因的基础上，因此对各样品得到的Unigene作进一步的聚类分析，最后得到慈竹的Unigene数据库。
1.3.3 功能注释使用Basic Local Alignment Search Tool (Blast)将组装后的Unigene序列与COG、GO、KEGG、Swiss-Prot和Nr数据库比对，得到最高序列相似性的Unigene。
1.3.4 基因的差异表达分析在NR、Swiss-Rort注释结果中进行关键字检索，记录“GeneID”以及相对应的注释结果。根据“GeneID”记录相应基因的相对表达丰度值，并将值导入MeV 4.9.0分析软件中，制作热图并进行聚类分析。
1.4 差异表达基因qRT-PCR验证利用primer premier 5.0对已测序正确的序列和慈竹内参基因Tublin进行qRT-PCR引物的设计。在iQ Multicolor Real-Time PCR自动扩增仪进行反应，每个样品重复3次，反应体系20 μL。以内参基因Tublin为对照，利用公式2-△△CT计算其相对表达量。
2 结果与分析2.1 测序数据结果统计利用Illumina HiSeqTM 2000平台，对供试材料进行转录组测序，共获得13.98 Gb的原始数据，得到69.28 M条的读长(reads) (表 1)。其中10、50和150 cm慈竹笋总碱基数均大于4 Gb，reads数大于20 M，GC含量分别为52.78%、53.12%和52.21%；而100 cm慈竹笋因测序深度较低，仅获得7 829 716条reads，GC含量为49.49%。同时，供试材料的不确定碱基数(N%)所占比例均≤0.03%，测序质量值≥30的碱基(Q 30%)所占比例均 > 80%，这表明测序质量可靠。
表 1 测序数据统计Table 1 Statistics of sequencing data
| Samples | Total reads | Total nucleotides (bp) | GC (%) | N (%) | Q20 (%) | Q30 (%) |
| 10 cm | 20 125 331 | 4 062 680 957 | 52.78 | 0.01 | 89.12 | 80.27 |
| 50 cm | 20 823 597 | 4 200 352 710 | 53.12 | 0.02 | 89.08 | 80.35 |
| 100 cm | 7 829 716 | 1 579 745 023 | 49.49 | 0.03 | 92.01 | 87.12 |
| 150 cm | 20 504 408 | 4 137 822 934 | 52.21 | 0.01 | 89.12 | 80.35 |
表选项
通过进一步分析后得到的非冗余Unigene库，含有111 137条慈竹笋的Unigene。其中200 bp-300 bp之间的Unigene，为41 814条，占37.62%；大于1 000 bp的Unigene为21 518条，占19.35% (图 1A，1B)。
 |
| 图 1 慈竹功能基因长度分布和所占比例 Figure 1 Distribution and percentage of Unigene length from Bambusa emeiensis. |
| 图选项 |
2.2 基因功能注释2.2.1 基因功能注释通过数据库比对，得到最高序列相似性的Unigene，其中有16 151条Unigene被匹配到COG数据库中；52 872条Unigene被匹配到GO数据库中；11 032条Unigene被匹配到KEGG数据库中；44 552条Unigene被匹配到Swiss-Port数据库中；而62 795条Unigene被匹配到Nr数据库中(表 2)。
表 2 慈竹转录组功能注释Table 2 Transcriptome function annotation of Bambusa emeiensis
| Database | Annotated number | 3 00≤Length < 1 000 | Length≥1 000 |
| COG_Annotation | 16 151 | 5 991 | 8 640 |
| GO_Annotation | 52 872 | 23 297 | 19 361 |
| KEGG_Annotation | 11 032 | 4 315 | 5 066 |
| Swissprot_Annotation | 44 552 | 19 210 | 18 551 |
| Nr_Annotation | 62 795 | 28 348 | 20 587 |
| All_Annotated | 63 094 | 28 491 | 20 599 |
表选项
2.2.2 功能基因的COG注释COG的功能注释与基因产物有直接的联系[12]。在慈竹笋Unigene数据库的COG功能注释中，有16 151条Unigene具有蛋白质功能定义。在COG功能分类中包含复制、重组、修复、信号转导、新陈代谢、分子过程等25个分类。其中一般功能基因代表最大的一类，有4 221条Unigene；其次是复制、重组和修复，有2 701条。此外，还涉及到多个植物生长发育相关以及纤维素、木质素生物合成相关的生理生化过程，如碳水化合物转运与代谢、氨基酸转运与代谢、辅酶转运与代谢、次级代谢产物生物合成、运输和分解代谢、细胞壁与细胞膜的生物合成、信号转导机制等过程(图 2)。
 |
| 图 2 Unigene的COG功能分类 Figure 2 Unigene of COG function categories. |
| 图选项 |
2.2.3 功能基因的GO分类GO数据库可对基因和蛋白质进行分类和注释，其对很多物种都是适用的[13]。GO数据库包含3大类，即细胞组分、分子功能和生物学过程。在慈竹笋Unigene数据库中，有52 872条基因注释到GO数据库中。其中第一大类中的细胞部分、细胞和细胞器所占比例较高，分别达到85.71%、84.85%和80.39%；第二大类中的结合和催化活性所占比例较高，分别达到68.90%和50.82%；而第三大类中的细胞过程和代谢过程所占比例较高，分别达到75.86%和73.31% (图 3)。
 |
| 图 3 Unigene的GO功能分类 Figure 3 Unigene of GO function categories. |
| 图选项 |
2.2.4 功能基因的KEGG注释KEGG数据库是有关基因产物代谢通路的主要数据库，基于代谢通路的分析对于进一步分析解读基因的功能非常有帮助[14]。在慈竹笋Unigene数据库中，有11 032条Unigene注释到KEGG数据库中。其中注释序列较多的3条代谢通路为RNA转运、植物激素信号转导与核糖体。而与植物纤维素生物合成相关的代谢通路，如淀粉和蔗糖代谢、光合生物的固碳作用和光合作用分别有257条、189条和105条注释到KEGG数据库中；与木质素生物合成相关的代谢途径，如苯丙烷类的生物合成和苯丙氨酸代谢分别有118条和110条在KEGG代谢通路中得到注释(表 3)。
表 3 Unigene数目最多的5个代谢通路及与纤维素、木质素生物合成相关的代谢通路Table 3 Top five metabolic pathways and the metabolic pathways of cellulose and lignin biosynthesis involving Bambusa emeiensis Unigene
| Name of metabolic pathway | Number of Unigene | Percent (%) |
| RNA transport | 408 | 5.47 |
| Plant hormone signal transduction | 407 | 5.46 |
| Ribosome | 404 | 5.42 |
| Spliceosome | 371 | 4.98 |
| Protein processing in endoplasmic reticulum | 345 | 4.63 |
| Starch and sucrose metabolism | 257 | 3.45 |
| Carbon fixation in photosynthetic organisms | 189 | 2.54 |
| Phenylpropanoid biosynthesis | 118 | 1.58 |
| Phenylalanine metabolism | 110 | 1.48 |
| Photosynthesis | 105 | 1.41 |
表选项
2.3 慈竹Unigene与其他物种的比较分析首先将慈竹的Unigene与水稻Oryza sativa、毛竹Phyllostachys edulis、短柄草Brachypodium、玉米Zea mays和高粱Sorghum的Unigene进行功能比对。所得结果如表 4所示，慈竹Unigene与毛竹匹配度最高，达到52.46%，可能是因为慈竹与毛竹同属竹类植物的原因。其次是与水稻相匹配达到50.41%，可能是因为水稻为模式植物，其基因序列信息较其他几种禾本科更全面。
表 4 慈竹Unigene与其他物种的比较分析Table 4 Analysis on Unigene between Bambusa emeiensis and other species
| Number of Unigene | Percent (%) | |
| With hits to Oryza sativa | 28 666 | 50.41 |
| With hits to Phyllostachys edulis | 29 831 | 52.46 |
| With hits to Brachypodium | 27 106 | 47.67 |
| With hits to Zea mays | 25 662 | 45.13 |
| With hits to Sorghum | 27 046 | 47.56 |
| With hits to all above | 21 524 | 37.85 |
| With no hits to all above | 23 062 | 40.56 |
表选项
根据表 4的结果，再将慈竹的Unigene与水稻转录本序列和毛竹编码序列数据库进行Blast比对，选择E-value < 1e-10。63 094个慈竹Unigene中34 127个(54.09%) Unigene与毛竹转录组相匹配，33 685个(53.39%) Unigene与水稻相匹配。在这些配对的序列中29 591个(46.90%) Unigene与毛竹和水稻都相匹配，4 536个(7.19%) Unigene只与水稻相匹配，4 094个(6.49%) Unigene只与毛竹相匹配。还有许多预测的Unigene与水稻和毛竹都不匹配，它们编码的蛋白质主要与生理过程、结合活性和催化功能有关，其中有大量预测的Unigene可能是慈竹特有的(图 4)。
 |
| 图 4 慈竹Unigene与水稻和毛竹的比较分析 Figure 4 Analysis on Unigene between Bambusa emeiensis, Oryza sativa and Phyllostachys edulis. |
| 图选项 |
2.4 纤维素、木质素生物合成相关基因与其他物种比对分析2.4.1 纤维素生物合成相关酶基因与其他物种比对分析将鉴定出的慈竹纤维素生物合成途径中的关键酶基因与水稻、短柄草、玉米、高粱和毛竹进行比对分析，毛竹的相关数据引用自Peng等[1]。由表 5可知，在Swiss-Prot数据库中，慈竹纤维素合酶Cellulose synthase (CesA)的数量与其他几个物种相比差异较大，CesA仅有9条。纤维素合酶相似蛋白Cellulose synthase-like protein (Csl)则是毛竹的最多有76条，慈竹仅有11条，与毛竹相比差异比较明显。而短柄草中最少，只在Swiss-Prot数据库中鉴定到7条Csl基因。蔗糖合酶Sucrose synthase (SuSy)数量最多的是水稻(21)，其次是短柄草(11)，慈竹只在Swiss-Prot数据库中鉴定出7条。
表 5 慈竹纤维素合成的功能基因与水稻、短柄草、玉米、高粱和毛竹的比较Table 5 Bambusa emeiensis functional Unigene of cellulose biosynthesis comparison with Oryza sativa, Brachypodium, Zea mays, Sorghum and Phyllostachys edulis
| Gene | Oryza sativa | Brachypodium | Zea mays | Sorghum | Bambusa emeiensis | Phyllostachys edulisa |
| CesA | 21 | 10 | 21 | 17 | 9 | 19 |
| SuSy | 21 | 11 | 7 | 8 | 7 | - |
| Csl | 18 | 7 | 10 | 13 | 11 | 76 |
| aThe results were cited from Peng et al. [1]. | ||||||
表选项
2.4.2 木质素生物合成相关酶基因与其他物种比对分析经Swiss-Prot数据库功能注释，筛选出188个木质素生物合成途径相关关键酶基因。与其他禾本科植物相比，慈竹苯丙氨酸解氨酶Phenylalanine ammonia-lyase (PAL)、4-香豆酸辅酶A连接酶4-coumarate-CoA ligase (4CL)、肉桂酰辅酶A还原酶Cinnamoyl-CoA reductase (CCR)、肉桂醇脱氢酶Cinnamoyl alcohol dehydrogenase (CAD)等7个关键酶基因均略少于水稻(表 6)。在Swiss-Prot数据库中，慈竹、短柄草、玉米和高粱的漆酶Laccase (LAC)基因是数量最丰富的，分别有17、32、20和20条，毛竹则未检测到。水稻的CCR，咖啡酰辅酶A-O-甲基转移酶Caffeoyl-CoA O-methyltransferase (CCoAOMT)和咖啡酸/5-羟基阿魏酸-O-甲基转移酶Caffeic acid-O-methyltransferase (COMT)基因数量较其他物种丰富，分别为45条、10条和10条。在慈竹中数量较少的是COMT，为2条。在毛竹中PAL最多，为8条，其余则较少，与慈竹相比差异较大。毛竹的相关数据引用自Peng等[1]。
表 6 慈竹木质素合成的功能基因与水稻、短柄草、玉米、高粱和毛竹的比较Table 6 Bambusa emeiensis functional Unigene of lignin biosynthesis comparison with Oryza sativa, Brachypodium, Zea mays and Sorghum and Phyllostachys edulis
| Gene | Oryza sativa | Brachypodium | Zea mays | Sorghum | Bambusa emeiensis | Phyllostachys edulisa |
| PAL | 7 | 4 | 5 | 3 | 5 | 8 |
| 4CL | 26 | 0 | 0 | 0 | 13 | 6 |
| CCR | 45 | 24 | 0 | 0 | 2 | 3 |
| CAD | 18 | 10 | 10 | 21 | 10 | 1 |
| LAC | 23 | 32 | 20 | 20 | 17 | - |
| CCoAOMT | 10 | 0 | 0 | 0 | 3 | 2 |
| COMT | 10 | 0 | 1 | 4 | 2 | 1 |
| aThe results were cited from Peng et al. [1]. | ||||||
表选项
2.4.3 转录因子相关基因与其他物种比对分析转录因子是由核基因编码的一类蛋白质，能选择性地激活或抑制某些基因的转录与表达[15]。从Swiss-Prot数据库中得到的转录因子如表 7所示。慈竹MYB和bHLH转录因子数量较水稻、短柄草、玉米和高粱多，且转录因子整体数量上与毛竹相比差异较明显。
表 7 慈竹转录因子相关基因与水稻、短柄草、玉米、高粱和毛竹的比较Table 7 Transcription factor of Bambusa emeiensis comparison with Oryza sativa, Brachypodium, Zea mays and Sorghum and Phyllostachys edulis
| Gene | Oryza sativa | Brachypodium | Zea mays | Sorghum | Bambusa emeiensis | Phyllostachys edulisb |
| MYB | 4 | 1 | 0 | 0 | 24 | 116 |
| WRKY | 133 | 82 | 119 | 100 | 47 | 116 |
| NAC | 120 | 71 | 101 | 87 | 23 | 125 |
| bZIP | 6 | 7 | 6 | 8 | 3 | 107 |
| bHLH | 1 | 1 | 1 | 1 | 54 | 168 |
| HD-ZIP | 0 | 0 | 0 | 0 | 0 | 51 |
| MADS-box | 97 | 87 | 79 | 78 | 23 | 16 |
| bThe results were cited from plant transcription factor database. | ||||||
表选项
2.5 差异表达基因分析2.5.1 差异表达基因的筛选与统计为了解4个不同高度慈竹笋的差异表达基因，采用DESeq软件，设定FDR < 0.01、Fold Change (FC)≥2作为筛选差异表达基因的关键指标。将10、50、100和150 cm四个样品进行两两比对后，得到以下结果，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条差异表达基因。将同一组合在不同的数据库间进行比较，提取并统计筛选得到的差异表达基因的注释信息(表 8)。
表 8 注释的差异表达基因数目统计Table 8 Summary of differential expression gene annotated
| Combination types | COG | GO | KEGG | Swiss-port | Nr |
| 10 vs 50 cm | 311 | 1 171 | 192 | 111 | 1 347 |
| 10 vs 100 cm | 271 | 1 089 | 253 | 972 | 1 254 |
| 10 vs 150 cm | 587 | 1 881 | 457 | 1 779 | 2 086 |
| 50 vs 100 cm | 222 | 807 | 213 | 725 | 961 |
| 50 vs 150 cm | 400 | 1 188 | 348 | 1 112 | 1 332 |
| 100 vs 150 cm | 268 | 799 | 255 | 720 | 919 |
表选项
2.5.2 纤维素生物合成相关差异表达基因分析慈竹纤维素含量直接影响着造纸效率。通过酸性洗涤纤维(ADF)法分析了4个不同高度的慈竹笋基部的纤维素含量。结果表明，在笋的早期生长阶段，随着笋的伸长，慈竹笋基部的纤维素含量逐渐增加(图 5)。通过进一步的转录组数据分析，可以从分子机理上了解慈竹笋在发育过程中纤维素含量的变化。
 |
| 图 5 不同高度慈竹笋纤维素含量 Figure 5 The content of cellulose in different height of Bambusa emeiensis. |
| 图选项 |
CesA在纤维素生物合成过程中起着关键的调节作用。转录组数据分析表明，CesA基因的RPKM值在10和50 cm的样品中聚在一簇；而100和150 cm的样品聚为另一簇，表明CesA基因随慈竹笋高度的增长呈现不同的表达模式(图 6A，6B)。其中，CesA1、CesA5、CesA6和CesA8同源的Unigene有着较高的相对表达量，但却随着笋的伸长，呈现出降低的趋势(图 6A，6B)。在4个不同生长高度的慈竹笋中，与CesA4、CesA7、CesA9同源的10条Unigene，虽然呈现着较低的表达，却随着笋的伸长，呈现出增高的趋势(图 6A，6B)。这些Unigene在GO数据库的功能注释中都涉及多个GO功能条目，包括细胞壁、初生细胞壁生物合成、次生细胞壁生物合成、纤维素合酶活性、纤维素生物合成过程、木质部发育、细胞壁增厚等，而在KEGG注释中，10条Unigene都参与的KEGG通路是K10999 (纤维素合酶) (表 9)。
 |
| 图 6 不同生长高度的慈竹笋中CesA的相对表达水平(RPKM值)和差异表达水平(Log2值) Figure 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
| Unigene | KEGG pathway | Annotation |
| CL10Contig1 | K10999 | Cellulose synthase A catalytic subunit 9 |
| T2Unigene23059 | K10999 | Cellulose synthase A catalytic subunit 7 |
| CL1669Contig1 | K10999 | Cellulose synthase A catalytic subunit 7 |
| T2Unigene28744 | K10999 | Cellulose synthase A catalytic subunit 9 |
| T1Unigene23998 | K10999 | Cellulose synthase A catalytic subunit 7 |
| T1Unigene22002 | K10999 | Cellulose synthase A catalytic subunit 9 |
| T4Unigene34723 | K10999 | Cellulose synthase A catalytic subunit 9 |
| CL1525Contig2 | K10999 | Cellulose synthase A catalytic subunit 4 |
| T2Unigene23060 | K10999 | Cellulose synthase A catalytic subunit 7 |
| T3Unigene29112 | K10999 | Cellulose synthase A catalytic subunit 7 |
表选项
2.6 差异表达基因qRT-PCR验证根据纤维素生物合成相关差异表达基因热图分析，筛选出CesA9 (T1_Unigene_BMK.22002)进行qRT-PCR验证。如图 7所示，在100 cm的慈竹笋中，CesA9的相对表达量大约是10 cm笋的9.5倍，表明其表达水平可能直接影响着笋的纤维素含量。
 |
| 图 7 CesA9 (T1_Unigene_BMK.22002)基因qRT-PCR Figure 7 Quantitative RT-PCR verification of CesA9 (T1_Unigene_BMK.22002) gene. |
| 图选项 |
3 讨论与结论测序技术日趋成熟，水稻[16]和小麦[17]等禾本科比较有代表性植物的基因组测序已经完成，并且已经转向实际应用方面的研究。在竹子中毛竹[18]基因组测序也已经完成，为竹类植物的分子生物学研究奠定了良好基础。转录组测序已经不再仅限于拟南芥等模式植物，越来越多非模式植物的转录组也开始得到研究，其对应的生长发育状态与机理也逐渐被人们所了解。本文利用转录组测序技术，对处于早期生长的4个不同高度的慈竹笋进行转录组从头测序，得到了非冗余的111 137条慈竹笋的Unigene。COG和KEGG功能注释表明，除了基本的功能外，这些Unigene涉及到蔗糖转运与代谢、次级代谢产物(如与木质素相关的苯丙烷类的生物合成)及细胞壁的生物合成等方面(图 2，表 3)。
值得一提的是竹类植物MYB类转录因子数量明显多于另外4个物种。有研究指出MYB转录因子家族能通过调控参与苯丙烷代谢途径的相关基因表达，进而影响木质部细胞的形成[19]。此外，MYB类转录因子也会调控部分纤维素生物合成相关的基因表达，如MYB46可直接调控与次生壁相关的部分纤维素合酶基因，从而影响纤维素的合成[20]。因此，进一步探讨MYB类转录因子在慈竹中的功能，对慈竹分子遗传改良具有重要的意义。
慈竹纤维素含量丰富，是较好的造纸原料，其与硬头黄[21]等竹种在四川已形成了以其为龙头的竹浆造纸工业产业链。竹材中的纤维素含量直接影响竹浆造纸的效率。10、50、100和150 cm高度的笋中，纤维素含量具有明显差异。因此本文进一步深入分析了4个样品中与纤维素合成相关的差异表达基因。纤维素的生物合成是一个高度复杂的生物过程，有多个酶参与其生物合成，其中包括CesA、Csl、SuSy等[22]。CesA是目前研究最多的纤维素生物合成基因，它编码的纤维素合酶不仅是纤维素生物合成的场所，而且在纤维素生物合成过程中起着关键的调节作用。纤维素合酶本身具有多基因现象，在纤维素生物合成过程中，多个CesA基因彼此相关，协同作用[23]。在慈竹Unigene数据库中，分别与Swiss-port和Nr数据库比对时，与CesA基因同源的Unigene分别有49条和33条；与Csl基因同源的Unigene分别有36条和27条；与SuSy基因同源的Unigene分别有14条和7条。这些基因序列与毛竹、水稻、玉米等物种的CesA、Csl和SuSy基因具有很高的同源性。进一步分析表明，不同高度慈竹笋的纤维素合酶CesA4、CesA7和CesA9基因存在差异表达且随着笋的伸长，总体呈上升的趋势(图 6)。有研究表明，竹类植物笋期的高速生长与笋的居间分生组织中细胞的快速生长与伸长有关，初生壁的发育确保了细胞的正常生长，而次生壁的发育所起到的支撑作用也确保了笋的快速生长[24-25]。另有研究表明，不同的CesA基因分别在细胞初生壁和次生壁合成时期起不同的作用[23]。水稻的OsCesA4、OsCesA7和OsCesA9是细胞次生壁合成所必需的，在木质部的细胞中高度表达。OsCesA4、OsCesA7或OsCesA9突变后会导致水稻茎杆变脆且纤维素含量降低，进而影响水稻的正常生长[26]。因此推测慈竹的纤维素合酶基因CesA4、CesA7和CesA9与细胞次生壁的合成有关且可能影响慈竹笋的生长发育，这有待进一步研究。
参考文献
| [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.DOI: 10.1038/ng.2569 |
| [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.DOI: 10.1371/journal.pone.0078944 |
| [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.DOI: 10.1371/journal.pone.0046766 |
| [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). 齐锦秋, 池冰, 谢九龙, 等. 慈竹纤维形态及组织比量的研究.中国造纸学报, 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). 刘庆, 何海, 沈昭萍. 成都地区慈竹生长状况及其与环境因子关系的初步分析.四川环境, 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基因的克隆及其生物信息学分析.西北农林科技大学学报:自然科学版, 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). 陈其兵, 蒋瑶, 卢学琴, 等. 四川不同地区慈竹的遗传多样性研究.西北农林科技大学学报:自然科学版, 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种丛生竹理化特性及纤维形态研究.植物研究, 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). 李秀兰, 林汝顺, 冯学琳, 等. 中国部分丛生竹类染色体数目报道.植物分类学报, 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.DOI: 10.1038/nbt.1883 |
| [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.DOI: 10.1093/nar/28.1.33 |
| [13] | Ashburner M, Ball CA, Blake JA, et al. Gene Ontology:tool for the unification of biology.Nat Genet,2000, 25(1): 25–29.DOI: 10.1038/75556 |
| [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). 张椿雨, 龙艳, 冯吉, 等. 植物基因在转录水平上的调控及其生物学意义.遗传, 2007, 29(7): 793-799.DOI:10.1360/yc-007-0793 |
| [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.DOI: 10.1016/j.plaphy.2013.02.007 |
| [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.DOI: 10.1016/j.ygeno.2015.04.003 |
| [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). 赵广枝, 孙化雨, 赵韩生, 等. 毛竹基因组测序及数据应用研究现状.世界竹藤通讯, 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.DOI: 10.1104/pp.105.060202 |
| [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.DOI: 10.1111/j.1365-313x.2012.05124.x |
| [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.(in Chinese). 陈其兵, 高素萍, 刘丽. 四川省优良纸浆竹种选择与竹纸产业化发展.竹子研究汇刊, 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.DOI: 10.1146/annurev-arplant-050213-040240 |
| [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). 甘小洪, 丁雨龙. 毛竹茎秆纤维发育过程中细胞壁的变化规律研究.林业科学研究, 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). 崔凯, 何彩云, 张建国, 等. 毛竹茎秆组织速生的时空发育特征.林业科学研究, 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.DOI: 10.1104/pp.103.022442 |
