郭由兵，程廷才，李开和，卢燕回，王卫峰，夏庆友*

（1.家蚕基因组生物学国家重点实验室，西南大学，重庆 400715；2.中国烟草总公司广西壮族自治区公司，南宁 530022）

打顶期烟草组织基因表达分析及打顶对腋芽基因表达的影响

郭由兵1，程廷才1，李开和2，卢燕回2，王卫峰2，夏庆友1*

（1.家蚕基因组生物学国家重点实验室，西南大学，重庆 400715；2.中国烟草总公司广西壮族自治区公司，南宁 530022）

打顶是烟叶生产重要的技术环节之一。为了研究烟草（Nicotiana tabacum）打顶期基因表达模式，选取广西大宁地区晒黄烟现蕾期10个打顶或非打顶的组织样本，构建链特异总RNA测序文库，利用Illumina HiSeq 2000完成测序。通过质量控制，获得9千多万条测序数据，分析获得50 706个基因的表达信息。筛选组织特异和打顶差异表达基因，共富集687个组织特异表达基因和444个腋芽打顶差异表达基因，为烟草组织基因表达模式和打顶分子响应机制研究提供数据支撑，也为调控烟草营养和生殖生长提供分子靶标。

烟草；打顶；转录组；生长调控

烟草（Nicotiana tabacum）是以收获烟叶为主的经济作物，在其营养生长后期，茎顶形成生殖中心，大量营养物质集中到顶端，影响烟叶产量和品质[1]。在烟叶生产过程中，打顶抹杈是抑制生殖生长和促进营养生长的有效措施，将引起烟草植株一系列生理生化响应，如增加营养元素的吸收利用、代谢物质的积累[1-2]。研究人员利用miRNA测序技术分析烟草根部对打顶响应的分子机制，发现打顶差异基因主要涉及激素代谢、转录、应激反应等[3-4]。杨洁等[5]报道打顶可影响腋芽内激素的含量和腋芽的生长，而杨惠娟等[6]报道了打顶差异蛋白能促进光合作用、增强植株抗性。随着测序技术的发展，3个烟草品种[Basma Xanthi (BX, Oriental)，K326 (Flue-cured)，TN90 (Burley)]的基因组测序完成，共鉴定到9万多个基因[7]。但迄今为止，打顶期全基因组范围的基因表达模式却鲜有报道。本研究采用RNA-seq技术，对烟草根、茎、叶、腋芽、花蕾、顶芽等组织进行转录组分析，比较打顶和非打顶样品的基因表达差异，为揭示烟草打顶响应分子机制和调控营养和生殖生长提供数据支撑。

1 材料与方法

1.1 试验材料

本试验选用广西大宁地区晒黄烟作为试验材料。收集现蕾期打顶3 d的根、茎、叶、腋芽等4份样品，同期未打顶根、茎、叶、腋芽、顶芽和花蕾等6份样品，共10份组织样本。所有样本用液氮保存，用于提取总RNA。

1.2 Illumina测序文库构建

参照TRIZOL Reagent (Cat#15596-018, Life technologies, Carlsbad, CA, US)提取总RNA，采用RNase-Free DNase Set (Cat#79254, QIAGEN, GmBH, Germany)纯化。使用Agilent 2100检查RNA完整性和纯度。采用TruSeq® Stranded Total RNA kit (Cat#RS-122-2201, Illumina, CA, US)构建Paired-end测序文库，使用Illumina HiSeq 2000完成测序。

1.3 转录组数据分析

使用软件fqtrim (v0.9.4, http://ccb.jhu.edu/ software/fqtrim/index.shtml)去除ployA，Bowtie（version 2.2.5）[8]去除rRNA序列，Trimmamotic（Version 0.36）[9]去除接头、低质量碱基和长度较短的读段(Reads)，fastqc (v0.11.2, http://www. bioinformatics.babraham.ac.uk/projects/fastqc/)检测数据质量。处理后的数据使用RSEM[10](v1.2.22)比对到烟草K326（Nicotiana tabacum K326）的参考基因集上，对样本中的基因表达定量。烟草K326的参考基因集数据下载自烟草基因组数据库(http://solgenomics.net)[11]。通过脚本程序abundance_estimates_to_matrix.pl将基因表达结果整合成表达矩阵，并采用TMM（M值的截尾均值，Trimmed Mean of M-values）方法均一化。采用软件edgeR[12]鉴定样本间差异表达的基因，脚本程序get_tissue_enriched_DE_one_vs_all.pl分析各样本富集的高表达基因。

1.4 GO和KEGG注释与富集分析

基因功能注释使用Trinonate流程(v2.0.2, http:// trinotate.github.io/)完成，采用KAAS（KEGG自动注释系统，KEGG Automatic Annotation Server）[13]对基因进行KEGG（京都基因和基因组百科全书，Kyoto Encyclopedia of Genes and Genomes）注释。通过WEGO (在线基因本体注释图示，Web Gene Ontology Annotation Plot)[14]比较打顶差异表达基因的功能。GO（基因本体，Gene Ontology）和KEGG富集通过程序GO_seq[15]软件完成。以上perl脚本来自Trinity[16]软件，相关图表由软件R和Excel生成。

2 结 果

2.1 转录组数据统计

经过质量控制的测序数据比对到参考基因集上，共获得50 706个基因的表达信息。由于参考序列仅为基因集，不能代表全基因组，因此各个样本的比对率较低（表1）。样本间比对率存在差异，如 打顶样本中，叶和根的比对率为23.38%和31.94%，其余样本的比对率都在35%以上。以FPKM（每百万测序片段覆盖每千碱基片段数，Fragments Per Kilobase per Million）值大于等于2为限，各样本表达的基因占总基因数的2/3（66.84%）以上，表达的转录本占总转录本数的1/3（36.63%）以上。为避免高表达和零表达基因的干扰，对FPKM矩阵进行log10(FPKM+1)转化（图1），结果显示基因在各样本中的表达相对均衡。上述结果表明，尽管各样本测序数据和比对率存在差异，但基因表达相对均衡，其结果可用于下游分析。

表1 RNA-seq数据统计Table 1 Summary of RNA-seq data

图1 各样本中表达基因分布Fig. 1 Distribution of expressed genes in each sample

2.2 组织特异基因表达模式分析

通过差异表达基因分析，获得烟草组织特异表达的基因，即在某一组织高表达且与其他所有组织都有显著差异的基因。利用组织特异表达基因的FPKM值绘制热图（图2），结果显示根部含有数目较多的特异表达基因，暗示地上部分和地下部分在基因表达上的差异，在转录水平上揭示了烟草根部组织的特异性。叶片中存在一类高表达的基因，且在花蕾、腋芽、芽和茎中也有表达，功能注释显示这些基因可能与光合作用及色素合成代谢相关。

2.3 组织特异表达基因功能注释

图2 组织差异基因表达模式Fig. 2 Expression patterns of tissue-specific genes

通过软件GO_seq对特异表达基因进行GO富集。相较于其他组织，根部有较多的特异表达基因，主要涉及响应刺激物（GO:0050896 response to stimulus）、物质运输（GO:0006810 transport）、氧化还原（GO:0055114 oxidation- reduction process）、天然免疫（GO:0045087 innate immune response）和细胞程序性死亡（GO: 0012501 programmed cell death）等。其中烟草生物碱合成过程中的腐胺N甲基转移酶全部在根部特异表达，和生物碱合成相关的腐胺、聚胺等物质合成相关基因也在根部特异表达，表明烟草根系是生物碱合成的重要场所。氧化还原过程中的血红素结合、氧化还原酶、谷胱甘肽转移酶等相关酶类，在根部都特异表达，暗示根部存在特殊的氧化还原保护机制。这与研究人员报道的生物碱和抗氧化基因在根部上调表达是一致的[3,7]。叶片特异表达的基因主要和光合作用相关，组成包括叶绿体、类囊体膜以及光合体系等组分，涉及光合作用、光呼吸、丙酮酸代谢等生物学过程，包括果糖二磷酸醛缩酶、醛裂解酶、果糖1,6-二磷酸1-磷酸酶等功能基因。这些基因的特异表达是叶片进行光合作用的分子基础。花蕾特异表达的基因主要涉及苯酸盐羟化降解、天冬氨酸代谢等生物学过程，相应的分子功能为酵母赖氨酸脱氢酶（GO:0004753 saccharopine dehydrogenase activity）、丙氨酸乙醛酸转氨酶（GO:0008453 alanine- glyoxylate transaminase activity）、alpha-L-阿拉伯呋喃糖苷酶（GO:0046556 alpha-L- arabinofuranosidase activity）和天冬氨酰合酶（GO: 0004066 asparagine synthase activity）等。这些基因的表达，暗示花蕾中存在复杂的物质代谢过程。

腋芽位于烟草茎上叶的腋部，具有萌发成分枝的潜力。腋芽组织特异表达的基因主要是与环境响应、转录和几丁质酶活性等功能相关，主要包括水份或压力应激相应，DNA结合转录调控因子等功能，可能与基因表达调控相关。参与几丁质酶降解过程的几丁质酶也在腋芽中富集表达，几丁质是仅次于纤维素的细胞内最丰富的有机化合物，其降解在细胞营养循环中发挥着重要作用。这些结果暗示打顶期的腋芽正加快物质循化和功能分化。

利用KAAS对基因集进行KEGG功能注释，共注释11 318个基因。组织特异表达基因的KEGG富集分析显示，叶中富集到的ALDO、rbcS、GAPA、FBP等基因与光合作用相关（表2），GGAT参与碳固定和氨基酸代谢[17]，DOX1参与长链脂肪酸a氧化的起始步骤，与植物抗菌或其他病原反应相关[18]。

表2 组织特异表达基因KEGG注释富集Table 2 Enrichment of KEGG annotation of tissue-specific genes

2.4 打顶差异基因分析

利用edgeR软件，以样本间基因表达差异倍数大于等于4，且基因在样本间差异显著的P值小于等于0.001（即P-vlaue≤1e-3和|log2(FC)|≥2）为标准，计算样本间的基因表达差异。比较叶、茎、腋芽、根四组打顶与未打顶样本的基因表达差异，发现腋芽中差异基因最多，茎、根中存在少量差异基因，暗示打顶对腋芽的影响远大于其他组织。同时，打顶腋芽和顶芽仅有16个差异表达基因，暗示打顶后腋芽的顶端优势抑制解除，替代顶芽成为新的生长点。以未打顶的腋芽为参照样本（图3），打顶腋芽和未打顶腋芽相比，上调152个基因，其中105个在芽中同样也上调，约占69.08%，下调292个基因，其中245个在芽中也同样上调，约占83.90%，暗示这些基因在顶芽萌发和打顶后腋芽发育中起着关键作用。

图3 打顶处理的腋芽和顶芽差异表达基因数量比较Fig. 3 Comparison of differentially expressed genes in apical buds and axillary buds upon topping

2.5 打顶差异基因的功能注释

比较腋芽打顶差异表达基因的GO功能注释（图4），结果显示打顶上调的基因主要分布在GO注释的细胞组分类目下，涉及细胞壁、细胞质部分、蛋白复合物、转运载体和形态学构成等。打顶下调的基因主要分布在GO注释的分子功能和生物学过程类目下，涉及的生物学过程有DNA结合，蛋白质结合，RNA合成和转录，应激和耐受性，激素响应和抗毒素的合成等。RNA合成和转录受抑制，延缓细胞功能的进一步分化；激素和抗毒素的合成下降，减缓顶端优势，促进腋芽生长。腋芽中的细胞组分类蛋白上调则是腋芽向营养生长转变的开始，暗示打顶后腋芽由先前的功能分化向营养生长转变。

图4 腋芽打顶差异表达基因GO功能注释Fig. 4 GO functional output of differentially expressed genes after topping in axillary buds

表3 腋芽打顶差异表达基因KEGG注释富集Table 3 Enrichment of KEGG annotation differentially expressed genes in axillary buds after topping

依据KAAS的注释结果，对腋芽打顶差异表达基因进行KEGG富集分析（表3）。与未打顶腋芽相比，打顶腋芽中的基因表达模式发生了显著变化。上调的组蛋白H2A、H4为核小体组装因子，可导致核内染色质结构重塑[19-20]；RPS10为线粒体中核糖体的组成部分之一，上调表达导致线粒体大小亚基比例失调，影响植株的休眠与复苏[21-22]；TIP和PIP为水通道蛋白，可维系细胞内水分平衡，其上调表达可增进细胞对水分的吸收[23-24]；LHCA3为光合作用触角蛋白，以较慢的速率吸收红光[25]，暗示打顶后腋芽中光合作用也受到影响。打顶腋芽下调的TPS和sacA参与海藻糖和蔗糖的代谢，有研究显示TPS与冬小麦耐寒特性相关[26]；WRKY33和HSFF为转录因子，WRKY33参与植物与病原菌的相互作用[27]；PRR5抑制植物的避荫反应，影响植物周期节律[28]。

3 讨 论

本研究采集打顶期10个组织样本，基于RNA-Seq技术，获得50 706个基因的表达信息。由于公共数据库中烟草基因组的注释信息并不完整，因此，将测序数据比对到基因集上，分析基因表达模式。分析中测序数据比对到基因集上的比对率较低，可能是由于参考基因集仅来自根和叶[7]，且根、叶与环境接触较多，取样过程中难免存在污染，同时基因表达具有的时期和条件的特异性，共同导致分析的组织中根和叶的比对率不高。即使如此，以FPKM值大于等于2为限，超过总基因和转录本的2/3（66.84%）和1/3（36.63%）是表达的。各样本表达基因数目和FPKM分布相对均衡，数据可以用于比较分析。

分析组织特异表达基因发现，根和叶中含有较多的特异基因（351和208个），这与根和叶具有活跃代谢活动的组织特性是相符的[7]。与其他组织相比，根部有较多的特异基因表达，在转录水平上揭示了地上组织与地下组织的区别。对特异表达基因功能注释发现，根部特异表达的基因多参与环境响应、氧化还原以及物质运输，从转录水平反映了根与环境的紧密接触和在植株摄取营养中的关键作用。另外，生物碱合成相关的基因在根部富集表达，暗示根部是生物碱合成的重要场所，这与先前的报道是一致的[7]。作为光合作用的主要场所，叶片中富集表达了一系列的酶、叶绿体组分、暗反应等光合作用相关的基因。花蕾作为生殖器官，富集表达的基因中较多地涉及苯酸盐羟化降解、天冬氨酸类氨基酸代谢等相关过程，暗示花蕾中存在较为复杂的生殖代谢过程[29]。与此对应，腋芽中特异表达的基因除与环境相关，还包括转录因子，能调控相关基因的表达，暗示腋芽有进一步分化或增殖的潜力[2,30]。

打顶抹杈是烟草生产的重要技术环节，可以打破生殖与营养生长的平衡，使较多的营养富集到叶片中，从而提高烟叶产量和品质。比较打顶对各组织中基因表达的影响，发现腋芽中含有较多的差异基因，根、茎及叶等组织中基因表达差异不明显，打顶后腋芽中的基因表达模式更趋向于顶芽。这些结果暗示打顶后如不及时抹杈，腋芽可能取代顶芽成为新的生长中心，根、叶等组织不再产生明显变化，无法实现提高叶片产量的目标，从转录水平证实了打顶后抹杈的必要[31-34]。

4 结 论

本研究采用转录组测序技术，对10个烟草组织样本进行转录组测序，获得超过5万个基因的表达信息，揭示各组织样本的基因表达特征以及植株对打顶的响应机制。结果表明根部主要富集生物碱等物质合成代谢相关基因，叶片富集光合作用相关基因。比较打顶差异表达基因，显示腋芽具有440余差异表达基因，明显高于其他组织样本，暗示打顶主要影响烟草腋芽生长发育。

[1] 林桂华，周冀衡，范启福，等. 打顶技术对烤烟产质量和生物碱组成的影响[J]. 中国烟草科学，2002，23（4）：9-13.

[2] 刘进平. 植物腋芽生长与顶端优势[J]. 植物生理学通讯，2007（3）：575-582.

[3] Qi Y, Guo H, Li K, et al. Comprehensive analysis of differential genes and miRNA profiles for discovery of topping-responsive genes in flue-cured tobacco roots[J]. FEBS Journal, 2012, 279(6): 1054-1070.

[4] Tang S, Wang Y, Li Z, et al. Identification of wounding and topping responsive small RNAs in tobacco (Nicotiana tabacum)[J]. BMC Plant Biol, 2012, 12: 28.

[5] 杨洁，胡日生，童建华，等. 打顶对烟草腋芽生长及植物激素含量的影响[J]. 烟草科技，2013（10）：72-75.

[6] 杨惠娟，王景，王红丽，等. 打顶前后烤烟叶片microRNAs表达差异的研究[J]. 中国烟草学报，2014，20（5）：110-115.

[7] Sierro N, Battey J N, Ouadi S, et al. The tobacco genome sequence and its comparison with those of tomato and potato[J]. Nat Commun, 2014, 5: 3833.

[8] Langmead B, Trapnell C, Pop, M, et al. Ultrafast and memory-efficient alignment of short DNA sequences to the human genome[J]. Genome Biology, 2009, 10(3): R25-R25.

[9] Bolger A M, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data[J]. Bioinformatics, 2014, 30(15): 2114-2120.

[10] Li B, Dewey C N. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome[J]. BMC Bioinformatics, 2011, 12: 323.

[11] Fernandez-Pozo N, Menda N, Edwards J D, et al. The Sol Genomics Network (SGN)--from genotype to phenotype to breeding[J]. Nucleic Acids Res, 2015, 43(Database issue): D1036-1041.

[12] Robinson M D, McCarthy D J, Smyth G K. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data[J]. Bioinformatics, 2010, 26(1): 139-140.

[13] Moriya Y, Itoh M, Okuda S, et al. KAAS: an automatic genome annotation and pathway reconstruction server[J]. Nucleic Acids Res, 2007, 35(Web Server issue): W182-185.

[14] Ye J, Fang L, Zheng H, et al. WEGO: a web tool for plotting GO annotations[J]. Nucleic Acids Res, 2006, 34(Web Server issue): W293-297.

[15] Young, M D, Wakefield M J, Smyth G K, et al. Gene ontology analysis for RNA-seq: accounting for selection bias[J]. Genome Biol, 2010, 11(2): R14.

[16] Haas B J, Papanicolaou A, Yassour M, et al. De novo transcript sequence reconstruction from RNA-seq using the Trinity platform for reference generation and analysis[J]. Nat Protocols, 2013, 8(8): 1494-1512.

[17] Liepman A H, Olsen, L J Alanine aminotransferase homologs catalyze the glutamate: glyoxylate aminotransferase reaction in peroxisomes of Arabidopsis[J]. Plant Physiol, 2003, 131(1): 215-227.

[18] Shimada T L, Takano Y, Shimada T, et al. Leaf oil body functions as a subcellular factory for the production of a phytoalexin in Arabidopsis[J]. Plant Physiol, 2014, 164(1): 105-118.

[19] Tolstorukov M Y, Goldman, J A, Gilbert C, et al. Histone variant H2A. Bbd is associated with active transcription and mRNA processing in human cells[J]. Mol Cell, 2012, 47(4): 596-607.

[20] Miotto B, Struhl K. HBO1 histone acetylase activity is essential for DNA replication licensing and inhibited by Geminin[J]. Mol Cell, 2010, 37(1): 57-66.

[21] Kwasniak M, Majewski P, Skibior R, et al. Silencing of the nuclear RPS10 gene encoding mitochondrial ribosomal protein alters translation in Arabidopsis mitochondria[J]. Plant Cell, 2013, 25(5): 1855-1867.

[22] Majewski P, Woloszynska M, Janska H. Developmentally early and late onset of Rps10 silencing in Arabidopsis thaliana: genetic and environmental regulation[J]. J Exp Bot, 2009, 60(4): 1163-1178.

[23] Venkatesh J, Yu J W, Park S W. Genome-wide analysis and expression profiling of the Solanum tuberosum aquaporins[J]. Plant Physiol Biochem, 2013, 73: 392-404.

[24] Reuscher S, Akiyama M, Mori C, et al. Genome-Wide Identification and Expression Analysis of Aquaporins in Tomato[J]. PLoS ONE, 2013, 8(11): e79052.

[25] Wientjes E, van Stokkum I H, van Amerongen H, et al. The role of the individual Lhcas in photosystem I excitation energy trapping[J]. Biophys J, 2011, 101(3): 745-754.

[26] Xie D W, Wang X N, Fu L S, et al. Identification of the trehalose-6-phosphate synthase gene family in winter wheat and expression analysis under conditions of freezing stress[J]. J Genet, 2015, 94(1): 55-65.

[27] Eulgem T, Rushton P J, Robatzek S, et al. The WRKY superfamily of plant transcription factors[J]. Trends Plant Sci, 2000, 5(5): 199-206.

[28] Takase M, Mizoguchi T, Kozuka T, et al. The unique function of the Arabidopsis circadian clock gene PRR5 in the regulation of shade avoidance response[J]. Plant Signal Behav, 2013, 8(4): e23534.

[29] Bokvaj P, Hafidh S, Honys D. Transcriptome profiling of male gametophyte development in Nicotiana tabacum[J]. Genom Data, 2015, 3: 106-111.

[30] 杨洁. 植物激素与烟草腋芽生长的关系及其调控研究[D]. 长沙：湖南农业大学，2013.

[31] 杨惠娟，许俐，王景，等. 打顶前后烤烟叶片蛋白质表达差异的研究[J]. 中国烟草学报，2012，18（5）：73-78.

[32] 陈常伟，姚忠达，张新华，等. 打顶时期及仿生型信号分子施用时间对抑制烟草生物碱效应的影响[J]. 中国烟草学报，2011，21（6）：32-36.

[33] 王茂. 打顶前后烟草miRNA表达谱的生物信息学分析及靶基因的电子克隆[D]. 郑州：河南农业大学，2010.

[34] 林中麟，石健林，周益. 烟草打顶研究进展[J]. 江西农业学报，2009（6）：32-36.

Analysis of Expressed Genes at the Topping Stage and Effects of Topping to Genes Expression in Axillary Buds in Nicotiana tabacum

GUO Youbing1, CHENG Tingcai1, LI Kaihe2, LU Yanhui2, WANG Weifeng2, XIA Qingyou1*

(1. State Key Laboratory of Silkworm Genome Biology, Southwest University, Chongqing 400715, China; 2. China Tobacco Guangxi Zhuang Autonomous Region Corporation, Nanning 530022, China)

Topping is one of the important technical procedures of Nicotiana tabacum production. In order to study the gene expression patterns at the topping stage in Nicotiana tabacum, we selected the Daning tobacco growing area in Guangxi and collected ten tissue samples at the topping stage. The samples were used to extract total RNAs and to construct stranded specific total RNA sequencing library and to sequence using the Illumina HiSeq 2000 sequencer. Through quality control, more than 91 million clean reads were obtained and a total of 50706 genes were detected to be expressed from ten tissues at the topping stage. After normalization, screening for tissue-specific and genes differentially expressed among tested tissues, a total of 687 genes were identified to be tissue-specifical and 444 genes showed differential expression before and after topping. Above data and results will not only provide some support for exploring the molecular responses of topping, but also provide molecular targets for the regulatiion of Nicotiana tabacum vegetative and reproductive growth.

Nicotiana tabacum ; topping; transcriptome; growth regulation

S572.03

1007-5119（2016）06-0014-07

10.13496/j.issn.1007-5119.2016.06.003

广西壮族自治区烟草公司科技创新项目“广西贺州大宁晒烟基因组测序、表达谱芯片和代谢组学研究”(GYK201226)

郭由兵（1989-），男，硕士研究生，研究方向为生物信息学。E-mail：gyoubing@yeah.net。*通信作者，E-mail：xiaqy@swu.edu.cn

2016-05-06

2016-09-30