江西铅山红芽芋和青秆芋的转录组比较分析

尹明华，曹 晴，陈 红，邓思宇，邓燕梅

(上饶师范学院 生命科学学院，江西 上饶 334001)

江西铅山红芽芋种植历史较为悠久，1 300多年前东晋时期江西铅山县就已开始种植红芽芋。至明代，万历《铅书》记载“斋菜有苋、茄、莴苣、芋、莲、姜、萝卜，皆能散秽去积”。更有《铅山县志》记录河口郊区农民于清咸丰五年(1855年)从事红芽芋生产，远销长江流域等地[1]。江西铅山红芽芋因芽红肉白而名，药食兼优，为江西省名优农产品，2013年核准为国家地理标志农产品[2]。食用，富含淀粉、蛋白质和脂肪，鲜美细腻，粉而不黏，味甘可口，不愧为芋类佳品[3]；药用，性平味甘，益脾补胃，调气解乏，含多种微量元素，可增强免疫、防癌抑瘤，含氟较高，可洁齿防龃护齿[4]。研究表明，红芽芋芋秆长而肥大，约占芋头重量30%，对人体有用的营养物质较为丰富，且含有优良的食用纤维，但其刺激性涩味较为强烈，容易造成刺舌、刺喉感，芋秆中一水合草酸钙针晶的特殊晶形及含量与芋秆涩味相关，因此，红芽芋芋头采收后芋秆大多作为废料处理，一般不用作食用[5]。但江西铅山县有种本土种植的青秆芋，地下芋头不大无收获价值，而地面部分繁殖快，割去也可快速长出，且芋秆肥大粗壮，食用无涩味，一般用作餐桌的美味佳肴。目前关于江西铅山红芽芋和江西铅山青秆芋的转录组分析，尚无相关研究报道。

随着第2代测序技术的发展，转录组测序(RNA-Seq)已逐渐成为植物分子生物学研究的重要手段，广泛应用于功能基因挖掘、分子标记开发、代谢通路和调控机制研究等研究[6-7]。转录组测序能够获取植物特定组织或器官在特定时期的基因表达信息，有利于挖掘品种之间的差异基因[8]。转录组分析成本低、数据量大、效率高、准确性高，可对组织或者细胞中所有RNA进行测序，并通过读段(reads)的拼接和丰度统计获得相应的转录本序列信息及其表达水平[9]，这些结果将为园艺植物分子标记开发提供重要序列信息，也为园艺植物种质资源遗传多样性分析、功能基因挖掘以及分子辅助育种等工作提供了重要基础[10]。江西铅山红芽芋芋秆刺激性涩味较为强烈，而青秆芋芋秆无涩味可作餐桌的美味佳肴，为了解其原因，本研究将以江西铅山红芽芋为对象，以江西铅山青秆芋为对照，利用高通量测序技术进行转录组基因测序，以期获得叶片中差异基因的种类、表达量及差异表达信息，探究江西铅山红芽芋芋秆具刺激性涩味而青秆芋芋秆食用无涩味的相关基因差异表达情况，为江西铅山红芽芋和江西铅山青秆芋的品种鉴定和分子育种奠定基础。

1 材料与方法

1.1 材料

江西铅山红芽芋试管苗(HYY组)和江西铅山青秆芋试管苗(QGY组)。

1.2 方法

1.2.1 RNA提取

参照说明书用Trizol试剂提取江西铅山红芽芋(HYY)和江西铅山青秆芋(QGY)试管苗(全株)总RNA。所提RNA的质量由 Nanodrop2000和琼脂糖凝胶电泳检测，Agilent2100测定RIN值。单次建库若RNA总量5 μg，浓度≥200 ng·μL-1，D260/D280介于1.8～2.2即为可用。

1.2.2 Illumina Hiseq测序

RNA测序送由上海美吉生物科技股份有限公司完成。利用带有Oligo(dT)的磁珠将3′端带有Poly A尾结构的mRNA从总RNA分离出来。然后利用Truseq TM RNA sample prep kit试剂盒构建cDNA文库。利用Illumina hiseq 4000 进行测序和分析。

1.2.3 数据分析

使用Trinity软件(https://github.com/trinityrnaseq/trinityrnaseq/wiki)对clean reads进行denovo拼接，形成一个不能在两端扩展的Unigene序列。利用Trinity软件组装得到转录本与参考蛋白质数据库NR进行BLAST比对。使用软件Trans Decode对组装结果进行 CDS 预测。再将序列注释到Swiss-Protein、KEGG和GO数据库中获得基因功能信息。利用RPKM(每百万个reads映射到外显子区的每千个碱基上的reads数)对Unigene的read count 值进行标准化。根据错误发现率(false discovery rate，FDR≤0.001)和表达量倍数变化(|foldchange|>1.5)的阈值筛选出两个样本间差异表达的基因(DEGs)。对显著性差异表达基因进行KEGG pathway显著性富集分析。

2 结果与分析

2.1 测序数据产出统计

经过测序质量控制，HYY组Raw reads为41 952 830，Raw bases为6 334 877 330，Clean reads为41 298 676，Clean bases为6 147 303 730，Q20和Q30分别为98.79%和96.13%，GC含量为53.09%；QGY组Raw reads为46 252 016，Raw bases为6984054416，Clean reads为45 551 860，Clean bases为6 745 871 513，误差率为0.0228%，Q20和Q30分别为998.92%和96.45%，GC含量为51.17%。

2.2 转录组从头组装结果评估

HYY组和QGY组组装得到的Unigene和Transcript的序列条数分别为51 941和91 532；组装得到的所有Unigenes和Transcripts的碱基个数分别为44 181 196和99 976 031；组装得到的最长Unigene和Transcript长度均为8 734 bp；组装得到的最短Unigene和Transcript长度均为201 bp；组装得到的所有Unigene和Transcript平均长度分别为850.60和1 092.25 bp。

2.3 测序数据与组装结果比对

HYY组Clean reads为41 298 676，Mapped reads为31 959 832，Mapped ratio为77.39%。QGY组Clean reads为45 551 860，Mapped reads为30 904 580，Mapped ratio为67.84%。

2.4 表达量统计

图1展示各样本表达Unigene/Transcript的表达量分布整体情况和表达量密度分布情况。HYY组和QGY组表达量密度在0～0.7。

图1 HYY组和QGY组表达量分布盒形图(A)和密度图(B)Fig.1 Box diagram (A) and density diagram (B) of expression quantity distribution in HYY group and QGY group

2.5 表达量样本间Venn分析

从图2可知，HYY组和QGY组表达的共有基因数为14 038，HYY组单独表达的基因数为17 843，QGY组单独表达的基因数为18 634。

图2 HYY组和QGY组表达量样本间Venn图Fig.2 Venn map between HYY group and QGY group

2.6 表达量样本间相关性分析

由图3可知，HYY组和QGY组表达量的相关系数为0.035，表明基因在HYY组和QGY组2个样本间的表达量相似度越低，即样本间相关性不够。

图3 HYY组和QGY组表达量样本间相关性热图Fig.3 Thermogram of correlation between HYY group and QGY group

2.7 差异表达分析

HYY组和QGY组共产生差异表达基因(differentially expressed gene，DEG) 32 555个，与QGY组比较，HYY组上调基因15 887个，下调基因16 668个(图4)。从表1可知，HYY组和QGY组前10个差异表达基因为假设蛋白GLYMA_13G018000、类丝氨酸/苏氨酸蛋白激酶D6PK、植物/T3A5-60蛋白质、非特征蛋白LOC110018203亚型X3、过硫双加氧酶ETHE1同系物(线粒体样)、肉桂酸-4-羟化酶、生长素应答蛋白IAA17、磷脂酰肌醇4-激酶γ7、非特征蛋白LOC105044622亚型X2。

图4 HYY组和QGY组表达量差异火山图(A)和表达量差异散点图(B)Fig.4 Volcano map (A) and scatter map (B) of expression difference between HYY group and QGY group

表1 HYY组和QGY组表达量差异统计表(前10)Table 1 Statistical table of expression difference between hyy group and QGY group (top10)

2.8 差异表达基因GO富集分析

HYY组和QGY组差异表达基因GO富集分析结果如图5和表2所示。从图5和表2可知，通过富集，多糖代谢过程、基因沉默、果胶分解代谢过程、碳水化合物代谢过程、生长素激活的信号通路、细胞外区、膜的固有成分、膜的组成成分、膜部件、凋亡细胞、作用于糖基键水解酶活性、水解O-糖基化合物的水解酶活性、蛋白质二聚活性、DNA结合、氧化还原酶活性对二酚及其相关物质的作用等GO term富集显著。

表2 部分差异表达基因GO富集结果Table 2 GO enrichment results of some differentially expressed genes

图5 HYY组和QGY组差异表达基因GO富集分析柱形图Fig.5 Column chart of GO enrichment analysis of differentially expressed genes in HYY group and QGY group

2.9 差异表达基因KEGG注释

从图6和表3可知，HYY组和QGY组差异表达基因富集的KEGG通路主要有植物激素信号转导、植物病原相互作用、半胱氨酸与蛋氨酸代谢、苯丙酸生物合成、皮质醇合成与分泌、过氧化物酶体、植物MAPK(丝裂原活化蛋白激酶)信号途径、不匹配修复、淀粉和蔗糖代谢、硫代谢、戊糖和葡萄糖醛酸相互转化、C型凝集素受体信号通路、甘油酯代谢、不饱和脂肪酸的生物合成、α-亚麻酸代谢、二萜生物合成、类黄酮生物合成、脂肪酸延长、光合作用、吲哚生物碱生物合成等。其中植物激素信号转导、植物病原相互作用、植物MAPK信号途径、淀粉和蔗糖代谢、苯丙酸生物合成、半胱氨酸与蛋氨酸代谢、过氧化物酶体注释的基因数均超过80,这可能是导致参试品种差异的主要原因。

图6 差异表达基因KEGG富集分析柱形图(A)和气泡图(B)Fig.6 Column chart (A) and bubble chart (B) of enrichment analysis of differential expression gene KEGG

表3 差异表达基因KEGG富集部分分析Table 3 Analysis of enrichment of differential expression gene KEGG

2.10 与芋秆涩味相关差异表达基因的筛选

芋秆中所含的草酸钙针晶为芋秆主要刺激性涩味成分之一。由表4可知，与芋秆涩味相关差异表达基因如线粒体样钙单转运蛋白2(calcium uniporter protein 2，mitochondrial-like)基因上调，而L-抗坏血酸氧化酶(L-ascorbate oxidase)基因下调，导致形成草酸钙针晶可能性下降。

表4 与芋秆涩味相关差异表达基因Table 4 Differentially expressed genes related to astringency of taro

3 讨论

转录组一般指所有表达基因的身份以及其转录水平，广义上指细胞内在某一生理条件下的所有转录产物集合(mRNA、核糖体RNA、转运RNA及非编码RNA)；狭义上特指所有mRNA的集合[11]。转录组测序有助于新基因的发现、基因功能注释，基因差异表达和分子标记筛选[12]。通过第二代测序——Illumina测序技术，能够全面快速地获得植物组织在某一状态下几乎所有转录本序列信息，已经广泛用于模式和非模式植物的研究[13]。

KEGG数据库是全面研究基因产物在细胞中的代谢途径和研究基因产物的功能的数据库，基于KEGG的代谢途径的分析有助于更进一步了解基因的生物学功能[14]。紫花苜蓿差异基因的GO功能分析和KEGG途径分析表明，核糖体、剪接体、内质网中的蛋白质加工和碳代谢途径可能导致2个品种紫花苜蓿的蛋白质含量差异，而植物激素信号转导代谢途径可能是导致2个紫花苜蓿品种中性洗涤纤维和酸性洗涤纤维含量差异的重要原因之一[15]。黄果、满天星和紫果3个品种百香果的KEGG富集分析表明，黄果与满天星的细胞胞吞作用有显著差异，黄果与紫果的玉米素生物合成和糖酵解途径差异基因富集显著，满天星与紫果的缬氨酸、亮氨酸和异亮氨酸生物合成中存在较大差异[16]。次生代谢物生物合成和萜类和聚酮化合物代谢途径可能与4个茶树品种(薮北、湘妃翠、牛皮茶和铁香茗)的苯甲醇合成有关，编码茶树苯甲醇合成途径可能有5条关键酶β-葡萄糖苷酶基因[17]。在本实验中，通过对江西铅山红芽芋和江西铅山青秆芋的转录组分析，结果表明，HYY组和QGY组差异表达基因富集的KEGG通路主要有植物激素信号转导、植物病原相互作用、半胱氨酸与蛋氨酸代谢、苯丙酸生物合成、皮质醇合成与分泌、过氧化物酶体、植物MAPK信号途径、不匹配修复、淀粉和蔗糖代谢、硫代谢、戊糖和葡萄糖醛酸相互转化、C型凝集素受体信号通路、甘油酯代谢、不饱和脂肪酸的生物合成、α-亚麻酸代谢、二萜生物合成、类黄酮生物合成、脂肪酸延长、光合作用、吲哚生物碱生物合成等。其中植物激素信号转导(175个基因)、植物病原相互作用(151个基因)、植物MAPK(丝裂原活化蛋白激酶)信号途径(118个基因)、淀粉和蔗糖代谢(118个基因)、苯丙酸生物合成(95个基因)、半胱氨酸与蛋氨酸代谢(84个基因)、过氧化物酶体(81个基因)注释的基因数较多，这可能是导致参试品种差异的主要原因。

植物激素实际上是一类信号小分子，植物的生长发育需要通过植物激素信号转导系统进行协调。一些因素包括内因或外因可能会诱导一系列植物激素基因的表达，继而作用于激素受体或组件，最终显现出不同性状[18]。江西铅山红芽芋和江西铅山青秆芋可能由于其植物激素(生长素、细胞分裂素、赤霉素、脱落酸、乙烯、油菜素内酯、茉莉酸和水杨酸等)在这些激素信号转导系统中有的受体或关键组分因互作(interact)或交互应答(crosstalk)会产生协同或拮抗的作用而使信号途径网络化，从而对江西铅山红芽芋和江西铅山青秆芋的生长发育、防御、环境适应等产生作用[19]。植物病原真菌是造成植物病害的主要病原菌之一，也是影响园艺植物产量和品质最主要因素之一。研究表明，植物细胞信号调控的普遍机制之一是蛋白磷酸化修饰，在植物—病原微生物互作过程中，关键调控蛋白的磷酸化状态直接影响到免疫信号的激活[20]。在江西铅山红芽芋和江西铅山青秆芋的植物—病原微生物互作过程中，可能依赖于蛋白激酶和蛋白磷酸酶来调节多种蛋白的磷酸化状态最终调控江西铅山红芽芋和江西铅山青秆芋的免疫途径，并且病原微生物效应蛋白可能还可改变植物代谢以满足侵染病原微生物的营养需要，这可能也是影响江西铅山红芽芋和江西铅山青秆芋品质差异的一个原因。MAPK信号转导通路由高度保守的丝/苏氨酸蛋白激酶MAPKKK、MAPKK和MAPK组成，广泛存在于真核生物体内[21]。近年来研究表明，在植物中，该信号转导通路通过其3个组分将胞外信号经细胞膜上的受体传递至下游应答因子，在个体生长发育及防御反应中发挥重要作用[22]。江西铅山红芽芋和江西铅山青秆芋的生长发育需要各细胞、组织和器官间的精准协调，蛋白质磷酸化可能参与了所有细胞信号转导机制，参与了细胞间以及细胞与外界信号间的交流，MAPK信号转导通路在江西铅山红芽芋和江西铅山青秆芋生长发育及器官发生中可能起到了核心作用。碳水化合物是植物光合作用的最终产物，碳水化合物含量高低体现了光合作用效率高低[23]。蔗糖和淀粉是光合作用同化产物的主要分配形式，蔗糖是光合产物的主要运输和贮存形式，为植物的生长发育提供能量，而淀粉是一种暂贮形式，供植物夜晚呼吸消耗[24]。淀粉和蔗糖的积累和转换不仅与物质基础的形成有关，而且也与产量和品质密切相关[25]。江西铅山红芽芋和江西铅山青秆芋可能通过增加体内蔗糖—淀粉代谢关键酶基因的表达，进而形成两者的品种特性和独特品质。苯丙酸类代谢是抗菌途径之一，通过苯丙酸类代谢能合成多种次生酚类物质，多属于植保素类物质(如酚类化合物、香豆酸酯类、类黄酮、木质素等)能抑制病原菌的生长，对植物抗病具有重要作用[26]。江西铅山红芽芋和江西铅山青秆芋的差异可能也在于苯丙氨酸解氨酶活性的大小以及木质素和绿原酸含量的高低。研究表明，植物蛋白普遍存在氨基酸不均衡的现象，有些植物缺乏蛋氨酸和半胱氨酸等含硫氨基酸[27]。蛋氨酸和半胱氨酸参与合成含硫维生素(如维生素 B1、生物素和硫辛酸等)、含硫辅酶(如辅酶A)、形成谷胱甘肽以及许多含硫防御性有机物(如金属硫蛋白、植物鳌合肽、蒜素和芥子油苷等)和含硫信号分子(如结瘤因子和膨胀素)等在新陈代谢上重要的有机化合物[28]。江西铅山红芽芋和江西铅山青秆芋可能通过蛋氨酸和半胱氨酸代谢广泛参与种芋萌发、幼苗发育、 胁迫响应和组织分化衰老等过程。过氧化物酶体是高度动态、代谢活跃的细胞器，主要参与脂肪酸等脂质的代谢及产生和清除不同的活性氧(ROS)[29]。因此，过氧化物酶体的动态平衡由ROS、过氧化物酶体蛋白酶及自噬过程调节，对于维持细胞的氧化还原平衡至关重要[30]。近年来，越来越多的证据表明，过氧化物酶体是ROS的重要来源，包含一套复杂的抗氧化防御机制来调控ROS的积累，ROS作为植物激素信号的第二信使，可以调节植物发育过程和应答反应[31]。过氧化物酶体、线粒体和叶绿体与ROS依赖的信号网络之间的相互作用是充分发挥江西铅山红芽芋和江西铅山青秆芋过氧化物酶体功能的关键。随着生物信息学发展与更新，有望开发出对江西铅山红芽芋和江西铅山青秆芋过氧化物酶体蛋白功能的预测软件，及通过信息技术与分子技术的结合可以多层面对ROS介导的途径如何与江西铅山红芽芋和江西铅山青秆芋生长发育、自噬过程及胁迫应答信号网络的整合有更深认识。

芋秆中所含的草酸钙针晶为芋秆主要刺激性涩味成分之一[5]。与芋秆涩味相关差异表达基因主要有线粒体样钙单转运蛋白2(calcium uniporter protein 2，mitochondrial-like)基因和L-抗坏血酸氧化酶(L-ascorbate oxidase)基因[32]。江西铅山青秆芋通过上调线粒体样钙单转运蛋白2(calcium uniporter protein 2，mitochondrial-like)基因活性，造成胞质钙离子浓度下降，并再通过下调L-抗坏血酸氧化酶(L-ascorbate oxidase)基因活性，这样以抗坏血酸作为前体物质通过氧化降解后生成草酸的可能也下降，最终导致江西铅山青秆芋形成草酸钙针晶可能性下降。这也是江西铅山红芽芋芋秆具刺激性涩味而青秆芋芋秆食用无涩味的根本原因。