藜麦叶片响应低氮胁迫基因的可变剪接事件分析

杨利艳，杨小兰，朱满喜，张玉荣，杨雅舒，王创云

(1.山西师范大学生命科学学院，中国 临汾 041000；2.山西农业大学农学院，中国 太谷 030800)

藜麦(ChenopodiumquinoaWilld)为一年生双子叶植物，苋科藜亚科藜属，是一种起源于南美洲安第斯区的古老物种[1]。藜麦种子富含无谷蛋白以及人体无法合成的8种必需氨基酸，且铁、钙、磷含量丰富，具有良好的营养健康效益[2，3]。藜麦不仅有极高的营养价值，还具有很强的抗逆性，因此受到世界的关注并被广泛引种，是未来最具潜力的农作物之一[4]。

氮元素(N)是细胞结构不可或缺的一部分，担负着植物的生长发育及机体调控的重要责任[5]。氮素约占植物干重的0.3%～5%，缺氮或低氮是导致农作物生长不良、减产的重要原因[6]。氮胁迫条件下，植物会通过根的伸长来增强氮的吸收[7]，在甘蓝型油菜[8]、小麦[9]、水稻[10]中皆有报道。然而目前对藜麦低氮响应的研究鲜有报道。

可变剪接(alternative splicing，AS)是一种转录后调节机制。该机制可使蛋白的最终产物表现出不同的功能和结构特性，或者在相同的细胞中可能由于表达水平的不同而导致表型发生变化，因此对增强机体内转录组的可塑性、基因的表达以及蛋白质结构的多样性起重要作用[11,12]。目前，利用RNA-seq 技术对可变剪接的分析有不少报道，可变剪接基因(alternatively spliced genes，ASG)在一些生物中占较高比例，如人(95%)、线虫(71%)[13,14]，在禽类[15]、甲壳类动物[16]中也发现了新的剪接变体。作为模式植物，拟南芥AS事件的发生比例也较高，为61%[17]，约有42%的内含子参与拟南芥的生长发育及胁迫响应[18,19]。KHOKHAR等[20]发现在不同生态类型的拟南芥中，许多反式剪接数量性状位点(sqtl)在昼夜节律钟、开花和应激反应基因中富集，表明差异可变剪接在调节这些重要过程中发挥了作用。KANNO等[21]在拟南芥中通过对一个prp4ka突变体的定量磷酸化蛋白组学分析，检测到几个富含丝氨酸/精氨酸的蛋白的磷酸化变化，这些蛋白调节组成型和选择性剪接以及其他剪接相关的因素。因此，研究可变剪接有利于对基因的深入了解。本研究对不同氮处理的藜麦叶片进行转录组测序，统计并分析可变剪接基因，旨在挖掘响应低氮胁迫的ASG，从而探究AS对植物低氮响应基因的调控作用。

1 材料与方法

1.1 材料

藜麦ZK7种子来自于山西省农业科学院作物科学研究所。将种子浸泡在10% NaClO溶液中消毒15 min，用蒸馏水冲洗3～5次，置于放有湿润滤纸的培养皿中，在25 ℃下培养发芽。待芽长约0.5 cm时，将其移入干净的砂子中置于光照培养箱育苗，光暗周期16 h/8 h，温度18～22 ℃。待藜麦幼苗生长至2叶，采用霍格兰全营养液培养，每3天浇一次，长至6叶期进行以下处理。

1.2 实验设置

选取长势相同的幼苗进行两组处理，即对照(CK)：全营养液培养(N浓度为7.1 mmol·L-1)；低氮培养(LN)：氮元素浓度为全营养液的20% (N浓度1.42 mmol·L-1)。每个处理设置3个重复，处理5天后取相同叶位材料，液氮速冻后保存。提取总RNA，转录组数据委托生物公司测定。

1.3 总RNA的提取与质量检测

采用试剂盒RNAprep Pure Plant Kit(TIANGEN)的方法，分别提取藜麦叶片的总RNA。用Nanodrop ND-1000检测RNA的浓度，Caliper LabChip GX精确检测RNA的纯度和完整性(RIN)。

1.4 RNA-seq文库构建和测序

RNA-seq文库构建和测序由北京百迈客生物科技公司完成，将对照组(CK)和低氮胁迫组(LN)两个样品检测合格后构建文库，之后进行Illumina HiseqTM测序。测序去除含有接头和低质量的Reads，获得高质量Clean数据。本研究通过分析已获得的Clean数据鉴定在低氮胁迫条件下发生AS事件的差异表达基因，从中找寻CK和LN中ASG的共性和差异性，进而理解AS对低氮胁迫的调控作用。

1.5 可变剪接事件的鉴定

采用StringTie[22]对Hisat2的比对结果进行拼接，通过ASprofile[23]软件分析可变剪接类型，细分为12类，并获取CK和LN中相应的可变剪接类型及表达量[24]。

1.6 低氮胁迫相关基因GO分类及KEGG功能富集分析

为筛选响应低氮胁迫的ASG，对LN中特有ASG进行差异分析[25]。并以Fold Change(表达量差异倍数)≥2，FDR(错误发现率)<0.01且q<0.05(q为校正后的P值)为筛选标准，由于不同处理组基因存在差异，根据CK和LN组中基因的表达量值绘制各样本中的基因表达量热图，采用Heatmap软件作图。GO分类根据GO数据库(Gene Ontology Consortium)获得分子功能( molecular function )、生物学过程(biological process )和细胞组成(cellular component)的GO编号，功能富集采用cluster Profiler分析(q<0.05)。KEGG通路分析以Pathway为单位，应用超几何检验，统计并分析其与整个基因组背景相比显著富集的Pathway。

2 结果与分析

2.1 RNA的提取与检测

将提取的RNA(RIN值≥7.0)采用Caliper Lab Chip GX检测，28S/18S均在2.0～2.1，A260/A280均在2.1～2.2(图1)。

图1 CK与LN总RNA的GX结果Fig. 1 The GX results of total RNA in CK and LN

2.2 可变剪接事件的鉴定

对CK和LN两个样品的转录组数据分析后发现CK中26 094个基因发生了AS，共发生68 764个，LN中25 562个基因发生AS，共发生66 611个。

在各种AS类型中，CK组以TSS和TTS两种类型为主，分别占所有可变剪接事件的40.95%和40.87%；此外AE，IR及SKIP三种类型分别占所有剪接事件的7.06%，3.97%及3.98%；剩下7种类型共占3.17%。低氮处理组中，LN组的可变剪接事件仍以TSS和TTS两种类型为主，但较CK组有所增加，分别占41.45%和41.11%；AE，IR及SKIP三种类型分别占6.18%，4.58%及3.68%，其中IR剪接类型较CK组增加15.4%；剩下7种类型共占3.05%(表1)。

表1 CK和LN基因组中鉴定的可变剪接事件类型及数量

图2 CK和LN中可变剪接重合及特有的基因数量 Fig. 2 Number of overlapped and specific alternative splicing genes in CK and LN

对CK与LN中发生的ASG进行整理对比后发现，有24 520个AS发生重合。然而，有1 573个新基因仅在CK中发生AS。同时，LN组中有1 041个新基因发生AS，在CK中却没有发生(图2)。

2.3 LN特有可变剪接基因差异表达分析

对LN中1 041个ASG进行分析，差异表达显著的有147个基因，对其聚类分析见图3。这147个基因中AS类型以TSS和TTS为主，其中，上调基因有72个，下调基因有74个，还有一个基因(NW_018743732.1)在TSS事件中下调，在TTS事件中上调。表明藜麦在低氮条件下的基因表达与这两种可变剪接事件密切相关。同时，这147个可变剪接基因的功能主要有无机离子运输和传递，次生代谢物的合成(NW_018742204.1及NW_018742461.1)、DNA解旋酶及修复蛋白(NW_018742412.1)、翻译后修饰(NW_018742673.1)、信号转导(NW_018743014.1)、细胞骨架的形成(NW_018742814.1)及细胞周期调控(NW_018742966.1)。

图3 LN差异可变剪接基因表达聚类热图 Fig. 3 Cluster thermography of differentially expressed genes

2.4 LN中差异可变剪接基因GO注释及功能分析

对LN基因组中147个可变剪接基因进行GO注释和功能分析，发现39个GO注释的基因(图4) 。其中，显著富集(q<0. 05)的GO项在生物学过程中有81个，在细胞构成中有26个，分子功能中有60个。在生物学过程中，发生可变剪接的差异基因主要集中于代谢过程(GO:0008152)和刺激反应(GO:0050896)；在细胞构成中，发生可变剪接的差异基因主要集中于核(GO:0005634)、细胞质(GO:0005737)、细胞壁(GO:0005618)以及质外体(GO:0048046)上；在分子功能中，发生可变剪接的差异基因主要富集于转移酶活性(GO:0016740)和水解酶活性(GO:0016787)相关GO项。

2.5 差异可变剪接基因KEGG通路分析

对LN中147个差异可变剪接基因进行KEGG通路分析(表2)，发现显著富集到DNA复制、异源末端链接、抗坏血酸和醛酸代谢、错配修复系统等通路上。DNA复制、异源末端链接参与新细胞的合成；抗坏血酸和醛酸代谢参与植物逆境响应、细胞生长；错配修复系统参与DNA双链上错配碱基对的矫正，泛素介导的蛋白水解调控细胞周期，囊泡运输的相互作用参与分子的转运，谷胱甘肽代谢参与抗氧化和解毒作用。

注:数字代表在该功能下注释的基因数；A：生物学过程；B：细胞构成；C：分子功能。图4 LN中特有可变剪接差异基因的GO注释Fig. 4 GO annotations for specific differential expressed alternative splicing genes in LN

3 讨论

选择性剪接在广大真核生物产生潜在的基因组信息方面发挥着重要作用，它是一个重要的转录后调节机制，可以增加蛋白质多样性，影响mRNA的稳定性[26]。然而，关于AS如何参与藜麦对低氮胁迫的响应，人们知之甚少。在本实验中，笔者在LN组中共识别到66 611个AS事件，包括TSS，TTS，AE，IR，SKIP等12种。其中TSS及TTS事件最为丰富，其次是AE，IR和SKIP。与CK组相比，鉴别出低氮处理组(LN)基因组中有1 041个新基因发生AS，这些基因的生物学过程富集在代谢过程和刺激反应等过程中。LEVIATAN等[27]发现拟南芥在低温胁迫下，其表达水平保持相对不变，但其剪接方式在异构体的相对丰度上有明显变化。大多数低温调控的选择性剪接将提前终止密码子(PTC)引入到转录本中，通过mRNA衰变(NMD)过程产生潜在的降解目标，这可能导致新蛋白质功能改变或产生负面影响。FANG等[28]发现水稻中11个OsAAT选择性剪接基因，其中1个基因剪接变异的表达受氮(N)的调控，5个基因在各种形态氮浓度较高时表达上调，9个基因在硝酸盐浓度较高时表达上调。剪接变异体的mRNA水平在水稻叶和穗发育过程中也受到调控，结果表明水稻OsAAT家族的可变剪接变异体在不同氮素条件下的表达率均有变化。本实验发现，藜麦响应低氮的差异可变剪接基因富集于与转移酶活性和水解酶活性相关的通路。GIARETTA等[29]研究表明由GGT1和GGT2编码的Apoplastic-谷氨酰胺转移酶活性对植物营养和生殖发育具有重要意义，因此笔者推测转移酶对藜麦缓解低氮胁迫起一定作用。对LN中差异可变剪接基因进行KEGG分析，发现显著富集到DNA复制、异源末端链接、抗坏血酸和醛酸代谢、错配修复系统等通路上，这些信号通路与藜麦氮吸收密切相关，表明在藜麦生长过程中可能通过差异可变剪接基因调节氮吸收来响应低氮胁迫。

4 结论

本试验利用转录组样品的高通量测序数据鉴定了藜麦响应低氮的特异可变剪接的差异表达基因，藜麦叶片主要通过TSS和TTS两种剪接方式形成差异基因，其中NW_018742556.1,NW_018742616.1,NW_018744348.1,NW_018742616.1,NW_018742733.1等基因通过DNA复制和异源末端链接等通路来响应低氮。