氮胁迫下小麦叶片转录组分析

张文云, 张建诚, 姚景珍*

(1.临汾市农业种子管理站, 山西 临汾 041000; 2.山西农业大学棉花研究所, 山西 运城 044000)

我国是世界小麦种植面积较大、总产量较高的国家之一，近年来种植面积已达2 427万hm2，总产量约为1.2亿t，分别约占世界的13%和19%[1]。氮肥作为提高作物产量和改善品质的主要栽培手段，在农业生产中起着重要作用[2]。在小麦生产中，为了追求高产往往过量施用氮肥，造成氮肥利用率降低，不仅导致资源浪费，而且给生态环境带来巨大压力。因此，深入研究作物体内的氮素循环和提高作物氮素利用效率的技术已成为当前国际的研究热点。

转录组(transcriptome)是某个物种或者特定细胞类型产生的所有转录产物的集合，包括信使RNAs、非编码RNAs和小RNA[11]。转录组学(tanscriptomics)是功能基因组学的重要组成部分，从整体水平上研究细胞中基因转录情况及转录调控规律，对特定条件下基因表达水平改变做定量分析[12]。转录组学获得的大量数据经过专业的生物信息学分析，既可以还原细胞内基因表达的特征，也可以获得多种基因表达调控的重要信息。目前，转录组学分析已应用于小麦在不同胁迫条件的响应[13-14]，孟晨[15]对耐盐碱小麦山融4号进行转录组分析，发现在碱胁迫条件下其根中有大量被诱导表达的基因，部分碱胁迫诱导表达基因提高了转基因植物的耐盐碱性。Li等[16]利用转录组测序分析田间小麦响应冷胁迫的分子机制，发现CBF和ABA通路对小麦的冷适应非常重要，且在小麦6A和6D染色体上与冷胁迫有关的差异表达基因分布最多。目前通过转录组测序研究小麦适应低氮胁迫的分子机制的相关报道较少。因此，本文对氮胁迫下的小麦叶片转录组进行了研究和分析，旨在阐明小麦在氮胁迫下的基因转录和调控规律。

1 材料与方法

1.1 材料培养及处理

供试小麦品种为晋麦47，其对低氮条件具有较好适应性[17]，由山西省临汾市农业种子管理站提供。小麦种子经75%乙醇表面消毒，并用蒸馏水冲洗4～5次，置于湿润滤纸在24 ℃黑暗条件下培养至萌发，然后转至人工气候箱光照培养(16 h光照/8 h黑暗)，生长至二叶一心时选取整齐一致的幼苗，分别转移至限氮营养液和氮充分营养液中生长。

氮充分营养液：0.2 mmol·L-1KH2PO4、1.0 mmol·L-1MgSO4·7H2O、0.75 mmol·L-1K2SO4、2.5 mmol·L-1CaCl2、0.000 05 mmol·L-1(NH4)6Mo7O24·4H2O、0.001 mmol·L-1H3BO3、0.000 5 mmol·L-1CuSO4·5H2O、0.001 mmol·L-1ZnSO4·7H2O、0.001 mmol·L-1MnSO4·H2O、0.1 mmol·L-1Fe-EDTA、1.0 mmol·L-1NH4NO3，pH 5.5。限氮营养液把NH4NO3浓度调整为0.05 mmol·L-1，其他成分与氮充分营养液相同。为保持小麦生长条件恒定，每3 d换一次营养液。

1.2 转录组测序

待小麦幼苗在限氮营养液中和氮充分营养液生长2周后取样，应用天为时代公司RNA提取试剂盒提取总RNA。由北京百迈客生物科技有限公司采用Illumina-Hiseq测序平台进行转录组测序。

1.3 转录组数据的组装

为保证后续组装质量，提高生物信息分析结果的科学性，应用SeqPrep和Sickle软件去除原始测序数据中的测序接头、引物序列和低质量值的读段，处理后得到高质量测序数据。

然后，对原始数据过滤得到的RNA-seq测序数据进行从头组装。使用软件为Trinity、Inchworm搜索方法建立K-mer graph (K=25)中全部可能路径；采用Chrysalis对上述得到可能路径进行分组后，对Chrysalis生成路径进行修剪和整合，保存主要路径。

1.4 差异表达基因分析

将测序得到的Reads与Unigene库进行比对，根据比对结果结合RSEM进行表达量水平估计。利用RPKM (reads per kilobase per million reads)表示对应Unigene的表达丰度[18]。

在差异表达分析中采用校正后的P值，即FDR(false discovery rate)<0.05,且差异倍数FC(fold change)≥2作为差异表达基因(DEGs)的筛选标准。

在GO(gene ontology)数据库对差异表达基因进行功能注释和富集分析。利用KEGG(kyoto encyclopedia of genes and genomes)对差异表达基因进行通路注释分析。应用STRING软件对上调和下调的差异表达基因进行了蛋白质相互作用分析。

2 结果与分析

2.1 测序数据分析

利用Illumina HiseqTM2500高通量测序，从氮充分营养液生长的小麦叶片(CKL)和限氮营养液生长的小麦(NTL)样本中分别获得53 920 024、53 817 412条原始序列，平均读长100 bp，原始碱基数分别为5 445 922 424、5 435 558 612 bp，碱基质量分值达到20分(测序错误率小于1%)的碱基比例均超过95%。

将低质量数据过滤后，从CKL和NTL样本中获得高质量序列分别为52 383 726和52 192 061条，高质量碱基数分别为5 108 134 511、5 091 803 373 bp，碱基质量分值达到20分(测序错误率小于1%)的碱基比例均超过99%。

2.2 转录组NR库注释

将获得的转录组在NR数据库中进行比对，发现有21%的isogene注释到粗山羊草中，15.1%注释到大麦中，13.6%注释到乌拉尔图小麦中，其余分别为二穗短柄草5.4%、水稻2.5%、小麦2.0%、高粱0.77%和玉米0.61%。15%的E值分布在10-20～10-5之间，6%的E值位于10-30～10-20之间，小于10-30的E值占78%。根据Blastx的结果可知，98%转录组序列的相似度大于60%，84%转录组序列的相似度大于80%。

2.3 转录组COG注释

如图1所示，基于序列同源性，27 194个isogene得到了COG注释。所参与的代谢过程共25个分类，主要包括：一般功能预测，复制、重组和修饰，转录，信号传递机制，翻译后修饰、蛋白翻转、伴侣，翻译、核糖体结构和生物发生，碳水化合物转运和代谢，氨基酸转运和代谢等。

2.4 低氮胁迫下差异表达基因的筛选

对转录组的表达量进行计算，并筛选差异表达基因(DEGs)后发现，低氮胁迫下叶片差异表达基因为1 267个，其中179个基因上调表达，1 088个基因下调表达(图2)。

图2 低氮胁迫下叶片差异表达基因散点图Fig.2 Scatter plot of differentially expressed genes under low nitrogen stress

2.5 叶片差异表达基因GO富集分析

对氮胁迫条件下小麦差异表达基因进一步进行GO富集分析，确定差异表达基因主要行使的生物学功能。叶片中差异表达基因在GO注释中共分为3个功能组(表1)：生物过程、细胞组分和分子功能。生物过程中，注释分类到与植物光合作用相关过程的差异表达基因较多；细胞组分中，差异表达基因主要分布在质体、叶绿体类囊体膜(腔)等部位。以上结果表明，低氮胁迫会导致小麦叶片光合作用等过程发生显著变化。

表1 低氮胁迫下叶片差异表达基因的GO功能注释Table 1 GO functional annotation of differentially expressed genes under low nitrogen stress

2.6 叶片差异表达基因KEGG通路显著性富集分析

叶片差异表达基因被注释到 KEGG 数据库6个一级通路的38个二级通路中，更进一步可细分到178个三级代谢或信号转导途径中(表2)。其中二级通路中氨基酸代谢(amino acid metabolism)、碳水化合物代谢(carbohydrate metabolism)、脂类代谢(lipid metabolism)、信号传导(signal transduction)、免疫系统(immune system)被注释的三级通路数目较多。氨基酸代谢、碳水化合物代谢和脂类代谢是植物生物合成的主要成员，其三级通路数目均为10左右，表明氮胁迫对这些通路影响较大，导致小麦的营养生长受到抑制。信号转导的三级通路数目为15，则说明氮代谢在三级通路水平上开始进行级联的响应，通过信号转导影响更多的代谢过程。

表2 叶片差异表达基因KEGG通路显著性富集分析Table 2 Enriched KEGG pathways of differentially expressed genes in leaf

注：A—RNA 加工与修饰；B—染色质的结构和动力学；C—能量产生和转化；D—细胞周期控制、细胞分裂和染色体分区；E—氨基酸的运输和代谢；F—核苷酸的运输和代谢；G—碳水化合物的运输和代谢；H—辅酶的运输和代谢；I—脂质运输和代谢；J—翻译，核糖体结构和生物合成；K—转录；L—复制，重组和修复；M—细胞壁/膜/包膜的生物发生；N—细胞运动；O—翻译后修饰，蛋白质更新，分子伴侣；P—无机离子的运输和代谢；Q—次生代谢产物的生物合成，转运和分解代谢；R—一般功能预测；S—功能未知；T—信号转导机制；U—细胞内运输，分泌和囊泡运输；V—防御机制；W—细胞外结构；Y—核结构；Z—细胞骨架。Note: A—RNA processing and modification; B—Chromatin structure and dynamics; C—Energy production and conversion; D—Cell cycle control, cell division, chromosome partitioning; E—Amino acid transport and metabolism; F—Nucleotide transport and metabolism; G—Carbohydrate transport and metabolism; H—Coenzyme transport and metabolism; I—Lipid transport and metabolism; J—Translation, ribosomal structure and biogenesis; K—Transcription; L—Replication, recombination and repair; M—Cell wall/membrane/envelope biogenesis; N—Cell motility; O—Posttranslational modification, protein turnover, chaperones; P—Inorganic ion transport and metabolism; Q—Secondary metabolites biosynthesis, transport and catabolism; R—General function prediction only; S—Function unknown; T—Signal transduction mechanisms; U—Intracellular trafficking, secretion, and vesicular transport; V—Defense mechanisms; W—Extracellular structures; Y—Nuclear structure; Z—Cytoskeleton.图1 氮胁迫条件下小麦叶片转录组COG注释Fig.1 Function classification in clusters of orthologous groups (COG) of proteins

2.7 转录因子

大量研究表明，转录因子在植物生长或发育的各个阶段以及适应各种逆境胁迫过程中都发挥着重要作用。低氮胁迫条件下，鉴定出转录因子类型有WRKY转录因子、MYB转录因子、乙烯反应转录因子和NAC等转录因子。经分析后发现，发生差异表达的共有23个，其中4个转录因子上调表达，19个下调表达(表3)。这些转录因子通过影响下游氮转运相关基因的表达，参与植物氮缺乏胁迫的应答调控[19-21]。

表3 转录因子差异表达基因Table 3 Differentially expressed genes of transcription factor

2.8 差异表达基因蛋白互作网络分析

蛋白质相互作用在生命活动中起核心作用，由蛋白质相互作用构成的PPI网络的拓扑特性分析是后基因组时代最重要内容。对上调和下调的差异表达基因进行了蛋白质相互作用分析，在叶片的上调差异表达基因中未找到蛋白质相互作用。如图3所示，下调差异表达基因的PPI网络分析中，BRADI2G40600.1与BRADI1G77050.1两个结点分别与14个、12个蛋白相结合。进一步分析表明BRADI2G40600.1结合点参与6个GO功能，分别为DNA修复(GO: 0006260)、核苷酸绑结合(GO: 0000166)、DNA结合(GO: 0003677)、ATP结合(GO: 0005524)、ATP酶活性(GO: 0016887)和核苷酸三磷酸酶活性(GO: 0017111)。

图3 差异表达基因的蛋白质互作网络Fig.3 Protein interaction network of differentially expressed genes

3 讨论

比较不同生长环境、组织、器官和发育阶段样本间的基因表达差异是揭示特定分子机制的有效方法。本文对氮胁迫下转录组的表达量进行统计并标准化处理、筛选后，发现低氮胁迫下叶片差异表达基因为1 267个，其中叶片中179个基因上调表达，1 088个基因下调表达。练兴明等[22]利用cDNA芯片技术分析了水稻幼苗中10 422个基因在受到低氮胁迫后的反应，3个时间点上共有471个(4.5%)基因表达量发生了改变，其中上调表达的有158个，下调表达的有355个。蔡红梅等[23]在研究水稻低氮胁迫时发现，所检测的32 341个基因中，有3 518个(10.88%)发生了差异表达，其中根中有2 214个，茎中有1 766个基因，根和茎共同基因有462个。本研究所检测的27 194个基因中1 267个(4.66%)发生差异表达，差异表达基因比例与上述研究结果近似，低氮胁迫的共同反应相同，但数目低于后者研究结果，这些基因可能是导致晋麦47氮高效利用的主要原因。

PPI网络分析表明BRADI2G40600.1 和 BRADI1G77050.1是氮胁迫下小麦叶片下调差异表达基因的结点。进一步分析表明BRADI1G77050.1与钙离子结合相关。钙离子是第二信使，参与真核生物对生物和非生物胁迫信号传导途径[28]。钙参与植物胁迫反应中各种生理过程调节，如水和溶质的运动，细胞分裂，细胞壁的合成，呼吸和置换[29]。研究结果表明，BRADI1G77050.1通过钙离子参与小麦氮胁迫反应。另外，BRADI2G40600.1在PPI网络中与14个蛋白相关联并且参与DNA复制，其在细胞循环和凋亡、植物生长发育和代谢过程中发挥重要作用[11,30]。