干旱胁迫下花生转录组与miRNA测序及相关基因的表达

徐 媛，陈锦玲，陈玉梅，李璐璐，李惠敏，秦新民

(广西师范大学 生命科学学院，广西 桂林 541004)

0 引言

【研究意义】花生(ArachishypogaeaL.)是我国重要的油料作物，也是植物油和蛋白质的来源之一。2014年我国花生种植面积达460.39万hm2，总产量为1 648.2万t，成为产量最大的油料作物[1]。我国花生产区主要分布在干旱、半干旱地区。干旱是花生生产上影响最大的限制因素，据统计，由于干旱引起的花生减产占全国总产的20%以上[2]。干旱除引起减产外，还使花生品质下降、黄曲霉污染及病虫害加重等[3]。因此，开展花生响应干旱的分子机制研究，在提高花生抗逆性、增加产量和提高品质方面均极为重要。【前人研究进展】近年来，转录组测序技术在花生干旱研究中得到广泛应用，为进一步认识花生抗旱分子机制提供了重要依据。赵小波等[4]利用转录组测序技术对抗旱花生品种J11在干旱胁迫下转录因子的表达变化分析，发现236个差异表达转录因子。孙爱清等[5]利用 Solexa 高通量测序技术对干旱胁迫下的花生叶片 cDNA文库进行差异基因表达谱分析，发现9个类黄酮代谢相关基因在干旱胁迫下显著上调表达，推测类黄酮在干旱胁迫下可能发挥重要作用。万丽云等[6]通过对10个不同的花生材料苗期进行干旱-复水试验，筛选抗感材料进行转录组分析，发现抗感材料的差异表达基因主要富集在氧化磷酸化、光合作用和植物代谢途径。胡述浩[7]研究发现花生60条保守miRNA序列，其靶基因为392个，这些靶点包括编码转录因子和其他功能蛋白，如NAC转录因子、乙酰辅酶A羧化酶、色氨酸合成酶α链、叶绿体核糖核酸蛋白等，同时发现随着高渗溶液浓度升高，miR2111、miR482的表达量明显降低，而分别对应的靶基因表达量则显著升高，植株的抗旱能力也随之增加。【研究切入点】花生在干旱胁迫过程中，差异基因和miRNA的表达变化以及差异表达的miRNA所调控基因种类和功能尚需深入研究。【拟解决的关键问题】近年来，PEG6000常用于植物模拟干旱研究[8-11]，通过对干旱胁迫下花生根系进行转录组测序和miRNA测序，结合转录组测序数据中基因的表达水平和miRNA测序中靶基因的表达水平进行统计分析，旨在挖掘出与抗旱相关的关键基因和miRNA，为花生的抗旱分子机理研究与抗旱育种提供分子依据。

1 材料与方法

1.1 材料

花生品种为桂花37，广西农业科学院经济作物研究所选育的珍珠豆型高油酸花生新品种。

1.2 方法

1.2.1 花生种子干旱处理 花生种子经75%酒精消毒10 min，在28℃培养箱黑暗条件下蒸馏水浸泡10 h后种在灭过菌的蛭石中，浇灌蒸馏水，待其发芽后浇灌霍格兰营养液；用营养液配置15%浓度的PEG6000，待花生长至3叶一心时，选择生长健壮、长势一致的花生苗小心冲洗根部，将其进行15% PEG6000胁迫处理0 d、1 d、2 d和 3 d。收集处理后的花生根系，贮存在-80℃的冰箱，作为转录组和miRNA测序的材料。0 d、1 d、2 d、3 d处理后的材料分别以HS15-0、HS15-1、HS15-2和HS15-3命名。

1.2.2 花生总RNA的提取、文库构建及测序 样品总RNA的提取按照trizol试剂盒说明书进行。通过Agilent 2100 Bioanalyzer分析仪和紫外分光光度计对总RNA的质量、浓度和完整性进行检测。4个RNA文库(HS15-0、HS15-1、HS15-2、HS15-3)构建及测序由华大基因科技有限公司完成。

1.3 数据处理和分析

转录组测序和miRNA测序数据的处理分别参照文献[12-13]进行。转录组数据使用过滤软件SOAPnuke对BGISEQ-500平台测序获得的原始数据进行筛选,去除包含接头的reads、未知碱基N含量大于5%的reads,以及去除低质量的reads。得到clean reads之后,以花生基因组 (ftp://ftp.ncbi.nlm.nih.gov/genomes/Arachis_duranensis)作为参考基因组，使用HISAT (http://www.ccb.jhu.edu/software/hisat)将 clean reads比对到参考基因组序列,使用Cufflinks(http://cole-trapnell-lab.github.io/cufflinks)将比对结果组装构建转录本。miRNA测序将获得的sRNA经过blast或者bowtie比对到miRBase数据库，鉴定得到已知的miRNA。新miRNA采用华大自主开发的Mireap进行预测。

2 结果与分析

2.1 转录组序列的拼接与组装

转录组HS15-0、HS15-1、HS15-2和HS15-3 4个文库分别获得43 625 768、44 792 456、43 100 476和36 638 768条Clean reads，过滤后的Reads数比对到参考基因组(花生：NCBI_Aradu1.1 ftp://ftp.ncbi.nlm.nih. gov/genomes/ Arachis_duranensis/)，比对上参考基因组的占比分别是77.19%、74.76%、69.83%和51.78%。同时，测序数据表明具有蛋白编码潜力的转录本有13 387条，其中属于已知基因的具有编码潜力的转录本为11 590条，占86.57%，而属于新基因的具有编码潜力的转录本为1 797条，占13.42%。

2.2 转录组差异表达基因

按照FPKM法计算差异基因的表达量，当倍数变化>0时为上调表达，反之，则为下调表达。同时以log(FPKM)≥2、FDR≤0.001为条件，筛选显著差异表达基因。在HS15-0与HS15-1比较组中，差异基因总数为14 018个，其中上调和下调差异表达的基因分别为3 836个和10 182个；HS15-1与HS15-2比较组中，差异基因总数是8 811个，其中上调和下调差异表达的基因分别有2 068个和6 743个；HS15-2与HS15-3比较组中，基因总数是5 936个，其中上调和下调差异表达的基因分别为2 264个和3 672个。

2.2.1 差异基因的GO功能 将3个比较组的差异基因进行GO功能分类，在HS15-0与HS15-1组的差异表达基因中，共有44个功能小类：其中分子功能有12类，共28 479个差异表达基因；细胞组分有14类，共20 730个差异表达基因；生物过程有18类，共24 065个差异表达基因。HS15-1与HS15-2组的差异表达基因中，共有为44小类：11类分子功能、14类细胞组分及19类生物过程，差异表达基因分别为19 032个、13 898个和16 107个。HS15-2与HS15-3组的差异表达基因中，共有42小类：10类分子功能、14类细胞组分以及18类生物过程，差异表达基因分别为12 065个、8 841个和10 631个。

2.2.2 差异基因的GO功能富集 对各个比较组中的差异表达基因富集到的GO Term进行统计分析发现，HS15-0与HS15-1组的差异表达基因被富集到4 160个GO Term中，其中生物过程2 377个、细胞组分545个、分子功能1 238个；HS15-1与HS15-2组的差异表达基因被富集到3 867个GO Term中，其中生物过程2 210个、细胞组分500个、分子功能1 157个；HS15-2与HS15-3组中的差异表达基因被富集到3 398个GO Term中，其中生物过程1 946个、细胞组分472个、分子功能980个。从3个功能(生物过程、细胞组分和分子功能)被富集到的条目数发现，3个比较组中参与生物过程差异表达基因的数目最多，而细胞组分的差异基因数目最少。

2.2.3 差异基因的KEGG代谢通路 对HS15-0与HS15-1、HS15-1与HS15-2、HS15-2与HS15-3 3个比较组中的差异基因进行KEGG Pathway注释分类，共分为6个分支：细胞过程(Cellular Processes)、环境信息处理(Environmental Information Processing)、人类疾病(Human Disease)(仅限动物)、遗传信息处理(Genetic Information Processing)、有机系统(Organismal Systems)和代谢(Metabolism)。3个比较组差异基因在代谢分支中的注释的数目最多。

根据 KEGG Pathway 注释分类结果，采用函数进行富集分析得出，3个比较组中的差异表达基因共同富集到133个代谢通路中。在HS15-0与HS15-1比较组中，植物激素信号转导(Plant hormone signal transduction)、MAPK信号通路-植物(MAPK signaling pathway - plant)、核糖体(Ribosome)、植物病原体互作(Plant-pathogen interaction)、碳代谢(Carbon metabolism)通路中差异基因数目比较多，分别为368个、357个、341个、339个和319个。在HS15-1与HS15-2比较组中，核糖体(Ribosome)氨基酸的生物合成(Biosynthesis of amino acids)、碳代谢(Carbon metabolism)植物激素信号转导(Plant hormone signal transduction)和RNA转运(RNA transport)通路中的差异基因数目比较多，分别为251个、241个、207个和202个；在HS15-2与HS15-3比较组中，核糖体(Ribosome)、氨基酸的生物合成(Biosynthesis of amino acids)、碳代谢(Carbon metabolism)、植物激素信号转导(Plant hormone signal transduction)、RNA转运遗传信息处理(RNA transport Genetic Information Processing)通路中的差异基因数目比较多，分别为251个、241个、207个和202个；核糖体(Ribosome)、氨基酸的生物合成(Biosynthesis of amino acids)、碳代谢(Carbon metabolism)、RNA转运、植物病原菌互作(Plant-pathogen interaction)通路中的基因数目最多，分别为251个、172个、172个、144个和143个。

2.2.4 持续上调和下调表达基因 在3个比较组中共筛选出13个持续上调的差异基因，其中7个为HSP分子伴侣家族基因、2个RING E3泛素连接酶、1个为谷胱甘肽-S-转移酶、1个柱头特异性STIG1样蛋白1、1个免疫球蛋白结合蛋白和1个 DnaJ同源B亚家族成员13。同时筛选到13个持续下调的差异基因。其表达量及注释见表1。

表1 持续上调和下调基因的表达量及注释Table 1 Expression and annotation of continuously up- and down-regulated gene

2.3 样品miRNA的测序

对HS15-0、HS15-1、HS15-2 和HS15-3 4个样品进行的miRNA测序，分别获得30 182 740条、29 044 095条、29 702 900条和28 200 685条 Raw reads，去除低质量、没有3’接头序列、5’接头污染、小于18nt和包含poly A的序列，分别获得28 005 836条、26 862 969条、26 668 939条和25 970 037 条clean reads。将4个样本中所有的 clean reads 和已知的miRNA数据库(miRBase、Rfam、siRNA、piRNA、snoRNA等)进行比对，能比对上的 reads 分别有21 804 361条、21 797 056条、21 192 973条和15 074 860条，分别占 clean reads 的 77.86%、81.14%、79.47%和 58.05%。

2.3.1 已知 miRNAs的鉴定 将HS15-0、HS15-1、HS15-2和 HS15-3 4个样本文库中比对上基因组的sRNA比对到miRbase数据库中，分别鉴定到 28个、27个、27个和 26个已知miRNA，其中4个文库中共有的已知 miRNA有25个。

2.3.2 新miRNA的鉴定 在样本HS15-0、HS15-1、HS15-2和 HS15-3中分别鉴定出116个、113个、110个和106个新的miRNA，其中有105个新 miRNAs在4个样本中为共表达。

2.4 差异表达miRNA

2.4.1 差异表达miRNA的筛选 在HS15-0与HS15-1比较组中，差异表达的miRNA有23个，其中有7个为上调差异miRNA，16个为下调差异miRNA；在HS15-1与HS15-2比较组中，差异表达的miRNA有36个，其中11个为上调差异miRNA，25个为下调差异miRNA；在HS15-2与HS15-3比较组中，差异表达的miRNA有57个，其中7个为上调差异miRNA，50个为下调差异miRNA。在3个比较组中，下调差异表达的miRNA数目多于上调差异表达的miRNA。

2.4.2 差异表达miRNA的靶基因 利用Target Finder等软件对差异表达miRNA进行靶基因预测。HS15-0与HS15-1比较组中15个miRNAs为差异表达，共预测到 392个靶基因，其中5个已知 miRNA调控 330个靶基因，10个新 miRNA调控62个靶基因；HS15-1与HS15-2比较组中有25个差异表达 miRNAs，预测到的靶基因数为913个，其中6个已知miRNA调控440个靶基因，19个新miRNA调控473个靶基因；HS15-2与HS15-3比较组中的 44个miRNAs呈现出差异表达，预测的靶基因为1 718个，其中8个已知miRNA调控的的靶基因数为1 007个，而36个新 miRNA仅调控 711个靶基因。

2.4.3 差异表达miRNA靶基因的GO和KEGG富集 将3个比较组的差异表达 miRNA的靶基因进行GO功能分类。在HS15-0与HS15-1比较组差异表达miRNA靶基因中，生物过程、细胞组分和分子功能三大类所占的靶基因数目分别为495个、429个和246个；HS15-1与HS15-2比较组的靶基因数目分别是1 089个、860个和504个；HS15-2与HS15-3比较组的靶基因数目分别是1 479个、1 324个和766个。数据分析发现，HS15-2与HS15-3比较组中差异表达miRNA的靶基因数最多，HS15-0与HS15-1比较组最少。在GO分类中，3个比较组的功能分类模式相似。在生物过程分类中，3个比较组注释到差异表达miRNA的靶基因数主要分布在细胞过程(cellular process)和代谢过程(metabolic process)类别中；在细胞组分分类中，3个比较组中的靶基因主要集中在细胞(cell)、膜(membrane)、膜部分(membrane part)和细胞器(organelle)4个小类别中；在分子功能类别中，3个比较组中的差异表达miRNA靶基因主要分布在结合(binding)和催化活性(catalytic activity)类型中。该分类结果与转录组测序得到的差异基因注释分类结果一致。在HS15-0与HS15-1比较组中，共富集到76条代谢通路，其中代谢途径(Metabolic pathways)、次生代谢产物的生物合成(Biosynthesis of secondary metabolites)、昼夜节律-植物(Circadian rhythm - plant)和类黄酮生物合成(Flavonoid biosynthesis)通路中差异miRNA靶基因的数目比较多，分别为98个、80个、27个和21个。在HS15-1与HS15-2比较组中，共富集到91条代谢通路，其中代谢途径(Metabolic pathways)、次生代谢产物的生物合成(Biosynthesis of secondary metabolites)、内质网中的蛋白质加工(Protein processing in endoplasmic reticulum)和剪接体(Spliceosome)通路中差异miRNA靶基因数目比较多，分别为123个、53个、45个和44个；在HS15-2与HS15-3比较组中，共富集到105条代谢通路，其中代谢途径(Metabolic pathways)、次生代谢产物的生物合成(Biosynthesis of secondary metabolites)、植物激素信号转导(Plant hormone signal transduction)和胞吞作用(Endocytosis)通路中差异miRNA靶基因数目比较多，分别为205个、104个、69个和65个。数据分析发现，在花生干旱胁迫下，有较多的靶基因参与调控代谢和次生代谢产物的生物合成途径。

2.5 与花生干旱相关的miRNA及靶基因预测

通过miRNA的靶基因功能注释分析发现，有13条miRNA的靶基因注释为与干旱相关的转录因子或基因，其中9条为持续下调miRNA，涉及35个基因，源异型域-亮氨酸拉链蛋白(homeobox-leucine zipper protein, HD-Zip)编码基因1个、MYB类LHY基因2个，WRKY转录因子编码基因3个、谷胱甘肽S-转移酶(glutathione S-transferase，GST)基因2个、 SPL转录因子(Squamosa promoter binding protein-like)基因20个、 CCCH型锌指蛋白(zinc finger CCCH domain protein)基因7个(表2)。ahy-miR156a分别参与WRKY转录因子、谷胱甘肽S-转移酶和SPL转录因子基因的调控。

表2 花生干旱相关的靶基因功能注释及miRNATable 2 Functional annotation and miRNA of target gene related to peanut under drought stress

续表2

3 讨论

3.1 抗旱功能相关基因

植物热激蛋白(HSPs)在植物减轻逆境胁迫引起的伤害方面发挥作用[14]。王利彬等[15]发现转HSP70烟草的耐旱性和耐热性得到显著提高。WANG等[16]在研究大麦干旱胁迫差异基因时发现HSP90蛋白在不同处理上出现差异表达。西瓜在水分胁迫时热激蛋白也出现上调表达[17]。此外，在小麦研究中也发现干旱胁迫会诱导抗旱小麦材料中HSP60和HSP90表达上调[18]。最近有研究表明，泛素连接酶E3在植物的抗旱功能中发挥至关重要的作用[19-20]，持续上调的RING finger E3泛素连接酶可通过调控脱落酸(ABA)信号通路参与干旱胁迫响应[21]。XERICO是一种RING-H2E3泛素连接酶，ZENG等[22]发现，过表达XERICO能增强ABA的生物合成，从而提高拟南芥的抗旱性。谷胱甘肽-S-转移酶(GSTs)是一类重要的多功能超家族酶，GSTs可通过亲电取代反应、降低毒害作用和过氧化酶清除功能，增加植物适应各种逆境胁迫的能力[23-24]。如过表达PjGSTU1基因可提高烟草的抗旱能力[25]。近年来，DnaJ蛋白在植物干旱胁迫的研究方面取得了较大进展[26]。SCHAFLITNER等[27]鉴定出6 个受干旱胁迫诱导上调表达的DnaJ蛋白基因。RODRIGUES等[28]从甘蔗中鉴定出 3个DnaJ蛋白基因，并发现其中有 1个在耐旱品种中下调表达，而在不耐旱品种中则上调表达。THUMMA等[29]鉴定出1个桉树DnaJ蛋白，发现在干旱胁迫条件下，其3个等位基因的表达差异极其显著。研究结果表明，在HS15-0与HS15-1、HS15-1与HS15-2和HS15-2与HS15-3 3个比较组中共筛选出13个持续上调表达基因，分别为7个HSP分子伴侣家族基因、2个RING E3泛素连接酶、1个GSTs，1个柱头特异性STIG1样蛋白1、1个免疫球结合蛋白和1个 DnaJ同源B亚家族成员13。其中HSPs分子伴侣家族基因、 RING E3泛素连接酶、 GSTs和DnaJ基因很可能与花生应激干旱胁迫相关。而柱头特异性STIG1样蛋白1和免疫球结合蛋白基因也为持续上调表达，但目前尚未见这2个基因与植物干旱胁迫相关的研究报道。此外，13个持续下调表达的基因其功能注释与花生干旱胁迫的关系也不明确，有待深入研究。

3.2 花生根系抗旱功能相关miRNA

HD-Zip是植物特有的转录因子，其可通过与相关激素基因及其下游基因互作增强植物耐逆性[30]。研究中miRNA中novel_mir70的靶基因为同源异型域-亮氨酸拉链蛋白基因。随着胁迫时间延长，novel_mir70表达持续减少，通过负调控机制使得靶基因表达增强，有利于植物抗旱能力的提高。MYB-related transcription factor LHY属于生物钟节律基因，其主要功能不仅调控生物昼夜节律，还可以通过识别胁迫信号通路基因上游的调控元件，调控植物对非生物胁迫作出反应[31]。PATADE等[32]在甘蔗的PEG胁迫中发现MYB和ahy-miR159的表达水平彼此相反。

ahy-miRNA159可以通过控制ABI3的表达来调控 MYB33转录因子的含量，从而参与ABA信号通路响应干旱[33]。研究中靶向MYB-related transcription factor LHY的2条miRNA ahy-miR156b-3p和ahy-miR159均呈持续下降的表达模式，可能通过负调控来增强花生的抗旱性。WKRY基因家族是高等植物中最大的转录因子家族之一，其蛋白结构基本含有 1～2个 WRKY结构域，参与非生物逆境应答如干旱、涝渍、高盐、高温、低温、损伤等，如BhWRKY1 与BhGolS1(半乳糖醇合成酶)启动子的 W-box 结合，转录激活BhGolS1的表达，可以提高拟南芥的耐旱性[34]。研究发现ahy-miR156a、novel_mir103和 novel_mir115靶向4个WKRY转录因子基因，其中ahy-miR156a和novel_mir103在胁迫的1～3 d过程中表达量也表现为逐步下降模式，使得靶基因WKRY转录因子基因高表达，从而提高花生的抗旱能力。谷胱甘肽S-转移酶(GSTs)为细胞内多功能蛋白，在植物体内参与初生代谢、次生代谢、解毒以及非生物胁迫过程[35-36]。研究发现2个miRNA(novel_mir83和ahy-miR156a)的靶基因为谷胱甘肽S-转移酶基因，其表达量在胁迫过程中呈持续下降趋势。 该结果与转录组中GSTs基因表达持续上调一致，推测novel_mir83和ahy-miR156a通过负调控作用，增加GSTs基因的表达，增强花生的抗旱能力。

鳞状启动子结合蛋白样(Squamosa promoter binding protein-like，SPL)是植物中特有的转录因子。已有研究证明SPL可以增强植物在干旱胁迫下的耐受性[37]。在植物生长发育过程中miR156与SPL的表达水平通常呈负相关，miR156可以通过mRNA剪切或翻译抑制来调控SPL的表达[38]。研究共筛选到2个下调表达的miR156家族蛋白(ahy-miR156a和ahy-miR156b-5p)，所靶向的20个靶基因功能均注释为SPL转录因子。由此推断，干旱胁迫下ahy-miR156a、ahy-miR156b-5p通过下调，使其靶基因表达水平增强来提高自身的抗逆作用。CCCH型锌指蛋白(zinc finger CCCH domain protein)是一类重要的调控因子，在动植物生长发育以及对外界逆境的反应过程中发挥着重要的作用。徐倩倩[39]研究发现，PEG可以诱导玉米中的ZmC3H54的表达，干旱处理下过表达的 Zm C3H54 转基因水稻植株抗旱性增强。本研究中发现1个下调的ahy-miR167-5p，其7个靶基因功能注释为CCCH型锌指蛋白，由此推测干旱胁迫下ahy-miR167-5p可能负调控CCCH型锌指蛋白，从而减少逆境对植物的伤害。

以上几个家族的miRNA可能通过作用于各自的靶基因参与到花生的干旱胁迫中，但这些miRNAs参与花生干旱胁迫应答机制尚待进一步研究。

4 结论

研究采用15% PEG6000对花生幼苗进行了干旱胁迫处理，通过对材料的转录组和miNRA测序，在转录组测序数据中发现了13个显著持续上调的差异基因和13个显著下调的差异基因；在miRNA测序数据中筛选到8个持续下调表达的miRNA，分别为ahy-mir156-3p、 ahy-miR159、ahy-miR167-5p、ahy-miR156a、ahy-miR156b-5p、novel_mir103、novel_mir83和novel_mir70。这些差异表达基因和miRNA参与了花生的干旱胁迫过程。