观赏海棠‘火焰’果皮中可变剪接基因分析

段诗瑶，田 佶，张 杰

(北京农学院植物科学技术学院/农业应用新技术北京市重点实验室，北京 102206)

苹果(Malusdomestica)是中国栽培面积最广的果树[1]。苹果果实的品质包括果实的大小、形状、色泽、香气、糖分和有机酸的比例等，其直接影响苹果的经济价值。在产量大幅度提升的同时，苹果品质的下降严重阻碍中国苹果产业的发展[2-3]。提高苹果果实品质是维持其产业发展的重要发展方向，而果皮色泽是果实品质的重要判断标准。苹果果实发育从开花到成熟一般需要150 d左右，是一个漫长的过程，有大量基因参与其成熟的调控。通过基因组与表观组甲基化测序分析，发现SPL13、ACS8等基因在果实发育早期通过控制细胞数量的变化而可能对成熟期果实大小产生影响[4]。MdMYB9和MdMYB11则通过与bHLH3、bHLH33相互作用来促进苹果果皮花色素苷的合成[5-6]。利用单核苷酸多态性(single nucleotide polymorphism，SNP)标记，可以了解到苹果基因组单核苷酸的变异情况。苹果MdCoLBD1/2基因序列通过SNP分析，显示出LBD基因家族各基因的特异性，进而导致各基因功能的特异性[7]。SNP4299和SNP4432与苹果抗炭疽菌叶枯病相关的基因位点紧密连锁[8]。

可变剪接在植物生长发育过程中起着非常重要的调控作用，其通过不同的剪接方式使同一个mRNA前体产生两个或更多成熟mRNA，进而产生多种不同异构体来增加蛋白的多样性[9-10]。拟南芥中开花阻遏因子FLOWERINGLOCUSM通过可变剪接产生一种被称为FLM-β的剪接体，其受温度的调节进而影响拟南芥的开花时间[11]。番茄MADS-box基因Lemads1在转录时存在4种可变剪接体用于调节萼片及果实的形成[12]。可变剪接不仅可以通过产生蛋白异构体来调控植物生长发育，还可以通过改变转录本的可读框使终止密码子提前产生，进入降解途径。目前针对苹果果皮发育进程中的可变剪切研究较少[13]。

观赏海棠作为砧木在苹果属植物的生产及育种中发挥重要作用。其叶片和果皮的多色性为研究果皮发育过程提供良好的研究材料。该研究通过转录组测序(RNA Sequencing，RNA-seq)技术对观赏海棠‘火焰’5个发育阶段果皮进行转录组测序分析，比较‘火焰’在不同发育阶段与果皮成熟相关的共有及特有可变剪接基因，为探索苹果果皮发育的遗传调控机理提供新的思路，进而为苹果属植物的育种奠定理论基础。

1 材料与方法

1.1 测序数据

测序数据为课题组前期测定的苹果属观赏海棠‘火焰’5个果实发育时期测序数据(NCBI登录号：PRJNA546083)。观赏海棠‘火焰’果实选自北京农学院观赏海棠种质资源圃(40.l°N，116.6°E)长势大小一致的东南方向着生的5个不同发育时期的绿果品种‘火焰’(Maluscv. ‘Flame’)。在花后35、60、95、120和150 d进行果实的采摘，分别命名为S1、S2、S3、S4和S5。

1.2 建库及分析流程

分别利用Qubit©2.0荧光仪(Life Technologies，CA)和Nano Photometer©2000分光光度计(IMPLEN,CA)对观赏海棠果实样品中提取的总RNA进行浓度和纯度的测定。质检合格后利用富含Oligo dT的磁珠富集mRNA，经过mRNA片段化后，将mRNA反转成cDNA，连接adaptor上机测序(Illumina Hiseq x Ten Miseq测序仪)。下机后Raw data通过Trimmomatic[14]软件过滤得到Clean data。

直接提取基因的蛋白序列，用eggNOG-Mapper[15]进行功能注释后，构建数据库。

利用HISAT2[16]软件构建参考基因组的索引，将质量控制后的高质量数据序列与参考苹果基因组序列(http://www.rosaceae.org或http://www.ncbi.nlm.nih.gov)进行比对，通过StringTie将比对上的reads进行组装和定量。

1.3 可变剪接事件鉴定

基于各样品reads与参考基因组序列的Hisat2比对结果，使用GATK[17]软件识别测序样品与参考基因组间的单碱基错配，识别潜在的SNP位点，并分析这些SNP位点是否影响基因的表达水平或者蛋白产物的种类。变异注释(SNP、InDel)和预测变异影响是利用SnpEff[18]软件。根据变异位点在参考基因组上的位置以及参考基因组上的基因位置信息，获得变异位点在基因组发生的区域(基因间区、基因区或CDS区等)，以及变异产生的影响(同义突变或非同义突变等)。

根据SNP位点碱基替换方式的不同，可以将SNP位点分为转换(Transition)和颠换(Transversion)两种类型[19]。根据SNP位点的等位基因(Allele)数目，可以将SNP位点分为纯合型SNP位点(只有一个等位基因)和杂合型SNP位点(两个或多个等位基因)[20]。采用StringTie[21]对Hisat2的比对结果进行拼接，通过ASprofile[22]软件获取每个样品存在的可变剪接类型及相应表达量并将可变剪接类型分为可变5′端或3′端剪接(Alternative exon ends-5′，3′，or both，AE)，内含子滞留(Intron retention，IR)和外显子跳跃(Skipped exon，SKIP)，第一个外显子可变剪接(Alternative 5′ first exon-transcription start site，TSS)和最后一个外显子可变剪接(Alternative 3′ last exon-transcription terminal site，TTS)4种类型。

1.4 发生可变剪接基因的注释分析

KEGG富集采用KOBAS[23]软件进行分析。将5个时期(每个时期3次生物学重复)的共同发生可变剪接基因注释到KEGG数据库，分析其代谢通路。

2 结果与分析

2.1 测序质量分析

对‘火焰’5个发育阶段(花后的35、60、95、120、150 d，分别对应S1、S2、S3、S4和S5时期)的果皮进行测序分析。Clean reads的总长度10 430 201～13 583 562，有82.96%到88.64%的reads被映射到参考基因组上，见表1。GC含量均维持在47.10%～48.15%范围内。每个时期的3个生物学重复高度一致，且所有样本测序所得的纯净数据Q30值均在90%以上(表1)。RNA-seq数据质量较好，满足后续分析需要。

表1 测序数据质量统计Tab.1 Statistical analysis of RNA-seq data

2.2 观赏海棠‘火焰’果皮可变剪接事件类型鉴定

分析可变剪接事件(图1)发现在‘火焰’果皮发育的5个不同时期均存在可变5′端或3′端剪接，内含子滞留和外显子跳跃，第一个外显子可变剪接和最后一个外显子可变剪接这4种可变剪接类型。

其中，第一个外显子可变剪接和最后一个外显子可变剪接占总可变剪接事件的绝大多数，且二者比例非常接近。在S1期，这两种类型共占比88.44%；在之后的4个时期中，共占比均高达92%以上。可变5′端或3′端剪接，内含子滞留和外显子跳跃发生较少。

在整个发育过程中，发生可变剪接的基因数量(图1)趋势较平稳，随着发育阶段的不同略有减少。S1期发生可变剪接的基因数量最多为46 952个，S5期最少为43 669个，减少的数量占比7.19%。

2.3 观赏海棠‘火焰’果皮可变剪接基因数量分析

在‘火焰’果皮发育过程中,每个时期发生可变剪接的基因总数不同，各类型在不同时期的基因数

也不同。通过构建维恩图(图2)，共有15 506个基因是5个发育时期共同拥有的，占总可变剪接基因数(23 062)的67.24%。每两个相邻时期(S1和S2、S2和S3、S3和S4、S4和S5)共有的可变剪接基因数分别为18 298、18 217、18 326及16 946个。每个发育阶段特有的基因数依次是432、261、373、270和775个。

2.4 观赏海棠‘火焰’果皮发育中共有可变剪接基因的KEGG功能富集

将‘火焰’果皮整个发育过程中共有的15 506个可变剪接基因进行KEGG功能富集，选择显著富集基因最多的20条代谢通路(图3)。显著富集基因最多的是嘌呤代谢，其次是氧化磷酸化、糖酵解/糖异生和嘧啶代谢，同时发现柠檬酸循环、磷酸戊糖途径、光合作用等通路富集在果皮发育过程中。在果皮发育过程中，可变剪接过程可能主要与NADH、ATP的形成有关，参与呼吸作用、光合作用等多种生理过程。

3 讨 论

可变剪接是广泛存在于生物体内基因表达过程中的一种重要的调控机制，且在植物的各个组织器官及发育过程中有重要作用[24]。同一个pre-mRNA经过可变剪接可以形成不同mRNA剪接异构体，从而增加蛋白的多样性[25]。冯雅岚等[26]借助RNA-seq证明可变剪接主要在转录后水平对植物发育和逆境胁迫响应进行调控，进而在发育和非生物胁迫响应中发挥重要作用。孙鸿等[27]发现小麦中的4个ATG18基因均具有2种可变剪接方式，在抗白粉菌侵染和抗高盐、干旱、低温、黑暗和缺氮等逆境中起到重要作用。张盼娃等[28]对干旱胁迫下对玉米蛋白磷酸酶2C基因ZmPP2C26两个可变剪接体的功能进行分析，证明两个剪接体均可增加植物对干旱胁迫的敏感性。

可变剪接在植物生长发育过程中起着十分重要的作用，其在增加植物抗逆性方面被广泛验证[5-6,29-33]。可变剪接基因参与苹果果皮发育调控鲜有研究。大多对果皮发育的研究是集中在转录因子的调控作用上，MYB-bHLH-WD40复合物、NAC转录因子家族、WRKY转录因子家族以及部分锌酯蛋白及乙烯合成相关基因能够广泛参与苹果果皮的发育调控[5-6,34-37]。芦笋两性花发育过程中的可变剪切的分析证明可变剪切不仅能够在增加植物抗性中起重要作用，而且其对植物的整个发育过程非常关键[38]。

该研究中利用RNA-seq测序技术分析观赏海棠‘火焰’果皮在5个不同发育时期中发生可变剪接的基因，在不同的发育时期，主要发生可变剪切类型是第一个外显子可变剪接和最后一个外显子可变剪接，且二者比例非常接近。发生可变剪接的基因数量随着果皮的发育有所减少，但幅度不大，说明可变剪接在观赏海棠果皮的整个发育过程中普遍存在，且可能主要在果实发育的前期起十分重要的作用。

为了进一步明确这些可变剪接基因的作用，通过KEGG对‘火焰’果皮发育过程中5个阶段共有的可变剪接基因进行功能富集，发现显著富集的代谢通路主要是嘌呤代谢，氧化磷酸化，糖酵解/糖异生和嘧啶代谢。氧化磷酸化与多种生物化学代谢过程有着密切关系，尤其是糖酵解、柠檬酸循环、β氧化等。另外，果糖和甘露糖代谢、半乳糖代谢、磷酸戊糖途径等代谢通路也被发现。糖类主要负责提供生物体所需的能量，也可作为信号分子影响糖运输方向、寄生抗病性等多样的细胞反应[37]。这些证明可变剪接确实在植物免疫及增加抗性中存在关联[36-39]，并且可能借助信号转导在其他生理过程中发挥重要作用。