基于RNA-seq数据分析玉米抗虫响应基因的可变剪接事件

杨利艳， 张玉荣， 杨雅舒， 王美霞， 陈保国， 赵丽， 张丽光， 王创云

(1.山西师范大学生命科学学院,山西 临汾 041004;2.山西省农业科学院作物科学研究所,山西 太原 030031)

同一mRNA前体在产生成熟mRNA的过程中会在不同的排列方式下剪除内含子，这种机制最早由GILBERT[1]在1978年首先提出，称之为可变剪接(alternative splicing, AS) 。可变剪接被认为是导致蛋白质功能多样性的重要原因之一,它使一个基因产生结构和功能不同的mRNA和蛋白质[2-3]。选择性剪接事件在真核生物基因的转录过程中广泛发生，尽管最初有研究认为高等生物仅大约5%的基因发生可变剪接[4]，但随后在人和其他几种重要的模式生物研究中发现，可变剪接的发生概率远大于5%。在人类(Homosapiens)、线虫(Caenorhabditiselegans)、果蝇(Drosophilamelanogaster)、拟南芥(Arabidopsisthaliana)中发生AS事件的比例分别高达95%，71%，60.7%，61%[5-8]。基因组数据显示，由于选择性剪接的发生，拟南芥中约有42%的内含子基因参与了拟南芥的发育、生物及非生物胁迫反应[9-10]。可变剪接通过产生不同转录本，导致蛋白质的结构存在差异，并可能在蛋白质的功能和特异性方面进行修饰。可变剪接转录本编码的蛋白质中，约有超过一半与信号转导有关，且大部分可变剪接区域包含了蛋白质和其结合模体的相互作用位点，这些模体包括配体、DNA/RNA和其他的蛋白质[11]。通过可变剪接通过形成不同转录本通过调节信号途径中某些关键基因的表达，从而影响动植物生长发育、生理代谢以及抗逆能力。烟草中N抗性基因是植物抗性基因LRR的成员，通过选择性剪接可以形成2个转录本NS，NL。研究发现，转cDNA-NL和cDNA-NS烟草均不能对烟草花叶病毒(TMV)表现出完全抗性，只有cDNA-NS-内含子和含有N基因组序列的转基因植株，同时编码NS和NL转录本才表现出对TMV的完全抗性，即烟草N基因发生可变剪接所形成的2个转录本对于TMV的完全抗性是必需的[12]。拟南芥中蛋白磷酸激酶HAB1在ABA信号转导途径中起着关键作用，但若RBM25剪接因子丢失，其HAB1基因的表达量会下降，从而阻止ABA信号传导[13]。此外，可变剪接参与了拟南芥开花时间的控制，flowering time control locus a(FCA)前mRNA的选择性剪接参与了花期转换的调控[14-15]。由此可见，可变剪接是产生蛋白质结构和功能多样性的一个重要机制，鉴定基因可变剪接事件有利于全面研究基因的功能。近几年，已有多种方法用以鉴定AS事件[16]，包括基因芯片及RNA-seq 技术[17-20]。目前，利用RNA-Seq技术已分析了多个物种基因组的AS事件，包括人[5]、小鼠(Musmusculus)[21]、牛(Bovini)[22]、拟南芥[8]、中华猕猴桃(Actinidiachinensis)[23]等。本研究分别对黏虫取食前后玉米的叶片进行转录组测序，将两者发生可变剪接事件的基因进行比较分析，旨在挖掘响应抗虫取食、同时发生AS事件的基因，从而探究可变剪接对植物抗虫响应基因的调控作用。

1 材料与方法

1.1 材料

玉米(Zeamays)自交系B 73由中国农业科学院作物科学研究所提供。分别用2%次氯酸钠溶液消毒种子20 min，用蒸馏水清洗干净后，置于湿润滤纸上，于22～23 ℃培养箱，萌发后，置于25 ℃、相对湿度65%～80%、光暗周期16 h/8 h的光照培养箱，长出2 叶后，移至花盆，常规管理至5～6叶。

黏虫(Mythimnaseparata)来自中国农业科学院植保所培养的2～3龄幼虫。

1.2 方法

待幼苗长至六叶期，选长势一致、健康的幼苗，每株幼苗即为1个重复，每组设置3个重复。每株材料上接种4头黏虫，使其均匀分布在植株上，12 h后取材，以不接虫材料为对照(M0)。各处理组(M12)植株置于用木架(长50 cm、宽50 cm、高80 cm)中进行取食试验，外罩一层薄纱以防黏虫逃逸。12 h后取相同叶位材料并用液氮和干冰保存，提取总RNA，由生物公司进行转录组测序。

1.3 RNA的提取

采用Trizol法分别提取玉米试验组M12、对照组M0中叶片总RNA。将所提取的RNA以1% 的琼脂糖凝胶进行检测，并利用超微量紫外分光光度计对RNA进行定量分析。

1.4 RNA-seq文库构建和测序

RNA-seq文库构建和测序均由上海生工生物工程有限公司完成，分别取3个重复样品的总RNA等量混合后构建2个文库，即M12(试验组)、M0(对照组)，进行Illumina HiseqTM测序。测序得到的原始数据，里面含有带接头的、低质量的序列。为了保证信息分析质量，使用Trimmomatic进行对原始数据进行过滤，去除带N碱基的序列(N表示无法确定碱基信息)、 去除reads中的接头序列、去除低质碱基(Q<20)、使用滑窗法去除reads尾部质量值在20以下的碱基(窗口大小为5 bp)，得到Clean数据。本研究仅使用得到的RNA-seq测序数据鉴定响应黏虫胁迫的发生可变剪接事件的差异表达基因，分析M12与M0发生可变剪接事件的共性和特异性，研究可变剪接对抗虫的调控作用。

1.5 可变剪接事件的鉴定

使用HISAT2将质控后的测序序列与参考基因组进行比对，再运用 Cufflinks 软件[24-25]预测基因模型，使用ASprofile分析[26]，并根据每个样本预测出来的基因模型，对可变剪切事件进行分类。ASprofile分析存在12种可变剪接事件类型(图1)，即单外显子跳跃(skipped exon，SKIP)、单外显子跳跃(模糊边界)(approximate SKIP，XSKIP)、多外显子跳跃(multi-exon SKIP，MSKIP)、多外显子跳跃(模糊边界)(approximate MSKIP，XMSKIP)、单内含子滞留(intron retention，IR)、单内含子滞留(模糊边界)(approximate IR，XIR)、多内含子滞留(multi-IR，MIR)、多内含子滞留(模糊边界)(approximate MIR， XMIR)、可变5′ 或3′ 端剪切(alternative exon ends，AE)、可变 5′ 或3′ 端剪切(模糊边界)(approximate AE，XAE)、第1个外显子可变剪切(alternative 5′first exon，TSS)、最后1个外显子可变剪切(alternative 3′ last exon，TTS)。

A．单外显子跳跃和多外显子跳跃(灰色和红色均代表外显子)；B. 单内含子滞留和多内含子滞留(灰色代表外显子、红色代表内含子)；C. 可变5′ 或3′ 端剪切(灰色和红色均代表外显子)；D. 第1个外显子可变剪切(灰色和红色均代表外显子)；E. 最后1个外显子可变剪切(灰色和红色均代表外显子)。

A. Skipped exon (SKIP) and multi-exon SKIP(MSKIP)(Both gray and red represent exons); B. Intron retention(IR) and multi-IR (MIR)(Gray for exons and red for introns); C. Alternative exon ends (AE)(Both gray and red represent exons); D. Alternative 5′ first exon(TSS)(Both gray and red represent exons); E. Alternative 3′last exon (TTS)(Both gray and red represent exons).

图1 可变剪接事件的分类

Fig.1 Classification of alternative splicing events

1.6 抗虫相关基因分析及GO功能富集分析

为筛选出响应抗虫胁迫的可变剪接基因，对M12中特有可变剪接基因用DEseq进行差异分析。为了得到表达显著差异的基因，将筛选条件为qvalue<0.05(q为校正后的P值<0.05)且差异倍数|FoldChange|> 2(表达量差异倍数>2)。根据差异基因的表达量FPKM(fragments per kilobase of exon per million reads mapped)值，采用Heatmap软件绘制它们在各样本中的表达量热图。

为确定其基因功能，将显著差异表达的可变剪接基因ID比对到DAVID数据库，获得根据GO分类的分子功能(molecular function)、生物学过程(biological process)和细胞组成(cellular component)的GO编号，再使用clusterProfiler进行功能富集分析(qvalue<0.05，认为该功能存在显著富集)。

2 结果与分析

2.1 RNA的提取和检测

将所提取的RNA以1% 的琼脂糖凝胶进行检测(图2)，并使用超微量紫外分光光度计对RNA进行定量测定，其中A260/A280均为1.9～2.1。

2.2 可变剪接事件的鉴定

将M0与M12 2种材料的转录组数据分别与玉米基因组进行比对后，M0中15 701个基因发生了可剪接事件，共发生79 363个， M12中11 791个基因发生了AS， 共发生了39 385个。

在不同的可变剪接事件类型中，M0以TSS，TTS，IR，AE 4种类型为主，数量分别为21 148，18 805，10 252，8 283，占所有可变剪接事件的26.65%，23.69%，12.92%，10.44%，其他的8种类型占比为1.2% ～5.7%。在经黏虫取食后，M12主要发生的可变剪接事件类型没有变化，仍为TSS，TTS，IR，AE等4种，但可变剪接事件发生数量有所减少，分别为14 421，13 625，2 735，2 706，其比例为36.62%，34.59%，6.94%，6.87%，余下8种AS事件类型比例为0.2% ～ 4.3%(表1，图3)。

M：DNA mark；1, 2, 3: 样品RNA 。

表1 M0与M12基因组中鉴定的可变剪接事件

图3 M0与M12中各类型可变剪接事件所占比例

对M0与M12中发生可变剪接的基因进行统计和比较发现，M0与M12之间有10 670个可变剪接事件是重合的。然而，在M0中存在5 031个新基因发生可变剪接，这些基因在M12中没有发生AS事件。在黏虫取食后的M12基因组中存在相同现象，有1 121个新基因发生可变剪接，这些基因没有在M0中发生可变剪接事件(图4)。

2.3 M12特有AS事件基因差异表达分析

对M12 中1 121个可变剪接基因进行差异表达分析，得到177个显著差异表达基因。差异基因的聚类分析见图5。其中，上调基因有142个，下调基因有35个。在这177个差异表达基因发生的可变剪接事件类型中，主要为第1个外显子可变剪切(TSS)和最后1个外显子可变剪切(TTS)，表明玉米响应抗虫基因的表达与这2种剪切方式之间存在密切的联系。

值得关注的是，M12中的177个发生AS事件的差异表达基因有多个基因为抗虫相关基因。这些基因分别是Zm00001d048736,Zm00001d004555。Zm00001d022453，主要参与生物碱的合成；Zm00001d039394，主要参与萜类物质的合成；Zm00001d032079，主要参与植物-病原互作信号通路。还有一些与RNA剪接、RNA转运、RNA降解相关的基因。

2.4 可变剪接基因GO注释及功能富集分析

为了解可变剪接对玉米抗虫能力的调控，对经黏虫取食的M12基因组中这差异表达的177个发生AS事件的基因进行了GO注释和功能分析，统计P≤1×10-10的GO项，发现有132个基因存在GO注释(图6)。其中，在生物学过程、分子功能和细胞构成中显著(qvalue<0.05)富集的GO项分别有84，93，25个。在生物学过程中，发生可变剪接的差异基因主要集中于氧化还原过程(GO:0055114)和代谢过程(GO:0008152)；在分子功能中，较多可变剪接基因被富集于金属离子的结合(GO:0046872)、氧化还原酶的活性(GO:0016491)相关GO项中；在细胞构成中，发生可变剪接的基因主要集中于膜的组成部分(GO:0016021)。

图4 M0和M12共有和特有发生可变剪接基因数量

数字代表注释到该功能下的可变剪接基因数目.A：生物学过程；B：细胞成分；C：分子功能。The numbers represent the AS gene number annotated to the corresponding function.A: Biological process; B: Cellular components; C: Molecular functions.

3 结论与讨论

可变剪接在植物的生长发育、病虫害防御中扮演着重要的角色。KANG等[27]发现，拟南芥miR400基因在位于内含子306 bp所发生的1个可变剪接事件是由热应激引起的。玉米ZmPP2C26基因转录后加工过程中会发生内含子保留型可变剪接，产生2个转录本，第1个转录本内含子尽管没有改变其编码蛋白的PP2C功能，但有可能使其作用由细胞核和细胞膜延伸至叶绿体[28]。

本试验中，与M0相比，黏虫取食后的M12基因组中有1 121个新基因发生可变剪接。对差异表达的可变剪接基因进行GO分析表明，分子功能结果主要与氧化还原酶活性、转移酶活性、水解酶活性、DNA结合、ATP结合相关，而这些功能均与植株的抗性相关。如研究发现，过表达2-半胱氨酸氧化还原酶基因(2-Cysperoxiredoxins,Zj2-CP)的拟南芥与野生型拟南芥相比，过表达Zj2-CP基因植株增加了对干旱和盐胁迫的敏感性[29]，依赖于2OG-Fe(II)的氧化还原酶可以特异性调控种子的萌发过程[30]。对拟南芥精氨酸酶(AtArgAH1和AtArgAH2)基因与植物耐盐性之间的关系分析发现，2个精氨酸酶基因(AtArgAH1和AtArgAH2)表达与盐胁迫存在一定的应答关系[31]。生物过程结果主要与氧化还原过程、代谢过程、以DNA为模板的转录及调控相关。这些功能与次生代谢物的产生、植株抗性基因表达调控相关，从而增加植物的抵抗力。植物能够产生类萜、生物碱、多炔、类黄酮、精油、独特的氨基酸和糖类等多种性质各异的次生代谢物, 这些代谢物在作物保护上有着广泛用途[32]。

本研究发现的Zm00001d048736基因发生2个AE可变剪接事件，TTS和TSS事件各1个，主要参与生物碱的合成。大多数生物碱都具有毒性作用，是许多植物杀虫剂中最常见的主要成分[33]。生物碱对昆虫的影响是多方面的，可影响昆虫的神经系统，刺激或抑制昆虫取食，影响昆虫产卵等[34]。Zm00001d039394基因发生TTS和TSS事件各1个，主要参与萜类物质的合成。利马豆叶暴露于蜘蛛螨取食植物后所诱导的任何萜类挥发物质，可诱导利马豆叶PR蛋白以及LOX，PAL和FPS[1]等基因的表达[35-36]，将植物暴露于外源茉莉酸中，这些基因的表达水平也有所提高，表明这些基因会对植物的损伤胁迫作出响应[37-38]，即萜类物质可以诱导植物防御基因的表达。还有一些与RNA剪接、RNA转运、RNA降解相关的基因也发生可变剪接，这些基因通过产生不同转录本对其功能产生影响，通过影响基因的转录翻译水平而发挥功能。以上表明，可变剪接对于植物抗虫的调控不仅发生于植物的诱导抗性，还对参与转录的基因有所调控。

本试验利用转录组样品的高通量测序数据鉴定了玉米虫食胁迫后所特异发生可变剪接的差异表达基因，通过寻找可变剪接事件的发生规律和可变剪接参与植物基因表达调控的途径，为探究可变剪接对植物抗虫能力的调控，挖掘植物对虫害的适应机制提供了基础研究数据。