干旱胁迫下偏关苜蓿保护酶及转录组差异性分析

卫 凯，朱慧森，岑慧芳，程映杰

（山西农业大学 草业学院，山西 太谷 030801）

紫花苜蓿（Medicago sativa）是黄土高原地区的当家牧草，营养物质丰富，粗蛋白、维生素和矿物质含量丰富，氨基酸的组成比较齐全，适口性好[1]。其中，偏关苜蓿是山西省农业科学院畜牧兽医研究所和偏关县畜牧局选育登记的地方品种，具有抗寒、抗旱性强以及营养价值高等特点[2]。山西处于山丘地带，属中纬度东亚季风气候，离海洋相对较远，旱灾时有发生[3]。这使得苜蓿产量通常很低，因而，了解紫花苜蓿抗旱的生理和分子遗传机制对于紫花苜蓿的可持续生产至关重要。

干旱胁迫下，植物细胞中的酶促防御系统可以通过维持活性氧等物质的产生和清除的平衡来缓解细胞受到的损伤[4]，植物保护酶系统主要包括超氧化物歧化酶（SOD）、过氧化氢酶（CAT）、过氧化物酶（POD）和抗坏血酸过氧化物酶（APX）等[5]。保护酶活性的主要研究方法在20世纪70年代得到发明和建立[6]。张新兰[7]测定了不同品种苜蓿叶片在干旱胁迫下抗氧化酶活性的动态变化，结果表明，抗旱性强的品种在干旱胁迫下草产量更高，叶片CAT、POD和SOD活性较高。张翠梅等[8]研究了PEG模拟干旱胁迫下紫花苜蓿叶片的膜脂过氧化程度，结果表明，抗旱性强的苜蓿品种叶片的膜脂过氧化程度更低。李红等[9]采用PEG对12种不同抗旱性的紫花苜蓿进行干旱胁迫，结果显示，不同抗旱性品种的叶片保护酶活性不同。因此，SOD、POD和CAT等保护酶的活性大小可以作为苜蓿抗旱性评价的指标。

植物在干旱条件下，能够诱发全部组织水平的应答，不仅涉及形态学水平的变化调整，还包括细胞水平的应答[10]。转录组学（Transcriptome）研究可以反映在特定时间和特定环境下植物与正常条件下的差异基因表达情况以及所在通路，有助于了解植物适应逆境胁迫的机制，从而有助于发现抵抗和缓解逆境胁迫的相关基因[11]。植物感受干旱信号后，通过信号转导激发一系列生理生化变化，在分子水平调控特定基因表达的过程尤为关键[12]，最早是1997年VELCULESCU等[13]提出在转录水平上进行的基因表达差异分析实际上就是进行转录组研究。MA等[14]对紫花苜蓿Dryland品种在干旱胁迫下的转录组学分析，筛选出1 690个上调基因和3 827个下调基因，富集分析发现，蔗糖合酶、胺氧化酶和脂酰辅酶A还原酶相关基因在紫花苜蓿响应干旱胁迫的过程中具有重要研究价值。刘佳月[15]在干旱胁迫下紫花苜蓿中草5号的研究发现，差异表达基因在苯丙烷生物合成和谷胱甘肽代谢通路显著富集，其数量在与代谢途径相关的方面注释最多，在次生代谢物的生物合成上也有相对较多的注释数量。对甘农三号品种抗旱性的分析发现，紫花苜蓿可以通过调控与脂质代谢、氨基酸代谢、碳水化合物代谢和次级代谢等通路相关的基因表达来响应干旱胁迫，增强抗旱能力[16]。

基于前期开展的偏关苜蓿种子萌发期的耐受性、显微结构的响应特征以及关键酶的活性等研究工作[17-20]，发现偏关苜蓿表现出的抗旱潜力还有待挖掘，尤其在分子水平调控特定基因表达的过程[4]。本试验以偏关苜蓿为研究材料，分析其在干旱胁迫和正常条件下保护酶及转录组的差异性，挖掘与抗旱相关的基因，为培育抗旱性强、性状优良的苜蓿新种质提供理论依据。

1 材料和方法

1.1 试验材料

供试材料为籽粒饱满、外形完整、大小均匀的偏关苜蓿种子，由山西农业大学草业学院牧草种质资源库提供。

1.2 试验方法

试验在山西农业大学草业学院实验室中进行。将偏关苜蓿种子播种于营养钵，置于14 h光照、10 h黑暗的培养箱中，昼夜温度分别为25、20˚C，用Hoagland营养液每日定时浇灌。30 d后挑选生长情况一致的偏关苜蓿，在20%PEG溶液（用Hoagland营养液配制）中进行72 h的胁迫处理。未进行干旱处理的材料为对照组，胁迫72 h的为干旱组，采集偏关苜蓿叶片，用于SOD、POD、CAT、APX活性的测定及RNA提取。

1.3 保护酶活性的测定

SOD活性测定参照文献[21]进行，采用NBT光还原法测定其在560 nm波长处的吸光值，酶活性以抑制3 mL NBT反应液光化学反应的50%所需酶量为1个活性单位；POD活性测定采用愈创木酚显色法[21]，在分光光度计上测量470 nm下5 min时的吸光值，以每分钟每克样品鲜质量吸光值的变化为1.0表示1个活性单位；CAT活性测定参照高俊凤[22]中的H2O2分解法，以1 min内A240减少0.1的酶量记为1个酶活性单位，共计时4 min；APX活性采用高俊凤[22]的方法测定，在20℃下测定10～30 s内A290的变化，计算单位时间内抗坏血酸减少量及酶活性。

1.4 RNA提取及高通量测序

试验采用TRIzol法提取RNA。检测合格的样品采用Illumina Hiseq TM 2500进行高通量测序，将得到的原始数据进行过滤，再利用TRINITY[23]软件进行组装，采用CD-HIT[24]方法选择相似类中最长的作为唯一基因，得到基因的序列集合（Unigene）。

1.5 基因功能注释

将Unigene序列与NR（NCBI non-redundant protein sequences）、NT（NCBInucleotidesequences）、COG（Clusters of Orthologous Groups of proteins）、Swiss-Prot、KEGG（Kyoto Encyclopedia of Genes and Genomes）、GO（Gene Ontology）六大数据库进行基因功能注释[25]。

1.6 差异表达基因分析

基因表达量的计算使用RPKM法[26]（Reads per kilobase million）。差异基因的筛选条件为：FDR＜0.05和|log2FC≥1|。GO分 析 的 方 法 是GOseq[27]。利用KEGG数据库进行通路（Pathway）分析。

1.7 数据分析

试验数据采用Microsoft Excel 2010和SAS 9.0统计分析软件进行处理。

2 结果与分析

2.1 干旱胁迫下偏关苜蓿SOD、POD、CAT和APX活性变化

从表1可以看出，偏关苜蓿叶片的干旱组SOD、POD、CAT和APX活性相对于对照组均显著上升（P＜0.05），其中，干旱组的POD活性是对照组的4.87倍。

表1 对照组与干旱组叶片中保护酶活性比较Tab.1 Comparison of protective enzyme activities in leaves of the control and drought groups

2.2 测序数据统计

测序结果如表2所示，对照组和干旱组的Q20百分率分别达到98.09%和98.05%，GC百分率平均值分别为51.56%和51.08%，这些数据说明测序质量偏好。

表2 测序数据输出质量情况Tab.2 The output quality of sequencing data

2.3 转录组的拼装

过滤后的测序数据经拼装后最终得到79 794条Unigene，总长度为81201 250 bp，平均长度为1 018bp，N50为1 599 bp（表3）。采用组装序列的长度分布来衡量转录组组装质量，可知序列满足分析的要求。

表3 Unigene的长度及数量统计Tab.3 Unigene length and quantity statistics

2.4 基因功能注释

由表4可知，成功注释到56 911条Unigenes，注释到6个数据库的数量分别是53 963（94.82%）、50 745（89.17%）、36 588（64.29%）、35 212（61.87%）、23 866（41.94%）、40 252个（70.73%），NR数据库的注释百分比最高。

表4 Unigene在各大数据库中的功能注释情况Tab.4 Statistics of Unigene functional annotation accor ding to the major databases

利用NR数据库进行同源物种比对，结果发现（图1），与偏关苜蓿同源基因相似度最高的是高粱，相似率达38.60%；其次为玉米（Zea mays），相似率为26.10%；日本水稻（Oryza sativa）的相似率为14.30%；二穗短柄草（Brachypodium distachyon）的相似率为6.90%。

对所有的Unigenes进行SSR分析，共识别出7 868个SSR。其中，二碱基重复2 570个，三碱基重复4 931个，四碱基重复85个，五碱基重复169个，六碱基重复113个，各碱基的分布频率如表5所示。

表5 SSR分布频率Tab.5 SSR distr ibution fr equency

2.5 基因差异性表达分析

进行差异基因筛选，获得5 867个差异表达基因，上调基因占总数的46.28%，下调基因占总数的53.72%。对基因的表达量制作火山图进行分析，如图2所示。

2.6 差异表达基因的GO分析

对处理组和干旱组之间的差异表达基因进行GO功能分析,结果如图3所示，图中从左往右依次为生物附着、生物调节、细胞组分组织或生物合成、细胞过程、发育过程、定位系统的建立、生长、免疫系统过程、定位、代谢过程、多元有机体过程、多细胞生物过程、生物过程负调控、生物过程正调控、生物过程的调节、繁殖、繁殖过程、应急响应、节律过程、信号传导、单生物过程、细胞、细胞连接、细胞部分、细胞外基质、细胞外基质部分、胞外器、胞外区部分、高分子复合物、膜、膜部分、膜包围腔、拟核、细胞器、细胞器部分、共质体、病毒体、病毒体部分、抗氧化活性、结合、催化活性、电子载体活性、酶调节活性、金属伴侣活性、分子转导活性、核酸结合转录因子、营养贮藏、蛋白质结合转录因子、受体活性、结构分子、转运活性，注释到51个类别，3个部分中数目最多的类别分别是代谢过程（2 145个）、细胞（2 785个）、催化活性（2 035个）；进行GO富集分析，选择富集度高的前30个如图4所示，3个部分中富集最多的类别分别是作用于糖基键水解酶活性、染色体和DNA包装。

对照组与干旱组蛋白激酶、O-糖基化合物水解酶、过氧化物酶、作用于糖基键的水解酶、木葡聚糖内糖基转移酶、抗氧化酶和蛋白磷酸酶活性等相关基因功能有明显差异。基因功能注释表明，代谢过程、催化活性、金属伴侣活性、抗氧化活性和酶调节活性等与干旱胁迫联系紧密的生物学过程和分子功能具有明显的响应变化。

2.7 差异表达基因的KEGG富集分析

偏关苜蓿共有35 215条Unigenes可在KEGG中得到注释，占Unigene总数的61.88%，注释到代谢途径的数目最多，为9 106个。KEGG富集分析表明，有3 471个差异表达基因可注释到122条代谢途径。

统计分析可知，122条代谢途径中最显著的4个是次生代谢产物的生物合成、植物病原体相互作用、苯丙烷的生物合成和其他聚糖降解，最显著的前20个通路相关的差异表达基因占总量的59.69%。由表6和图5可知，信号转导、物质合成与降解、代谢途径和次生代谢等通路相关基因的数量最多。差异表达基因在植物激素信号转导和次生代谢物的生物合成富集程度相对较高。

表6 差异表达基因富集程度排名前20的Pathway条目Tab.6 Top 20 enrichment pathway in differentially expressed genes

从次生代谢物的生物合成、苯丙烷生物合成和植物信号转导等代谢通路中挑选出部分与MYC2转录因子、苯丙氨酸酶、类胡萝卜素结合蛋白、光合作用蛋白、淀粉酶、蔗糖合酶有关的基因作为偏关苜蓿响应干旱胁迫的候选基因，主要有Unigene12274_All、Unigene10105_All、Unigene35365_All、CL7158.Contig1_All、CL 11234.Contig1_All、CL 12445.Conti g1_All、Unigene12274_All、CL4041.Contig1_All。

2.8 抗旱相关转录因子的预测

对材料中与转录调控相关的差异基因进行了功能注释，共得到136个转录因子（Transcription factors，TFs）基因，被分为18个基因家族，其中，数目较多的是bHLH（29个）、WRKY（22个）、ERFs（18个）和MYB（11个）。筛选出表达差异倍数2以上的转录因子基因（表7），这些转录因子在植物响应干旱胁迫中可能发挥重要的作用。

表7 可能与抗旱相关的转录因子Tab.7 Possibly dr ought-r esistance-related transcr iption factors

3 讨论

3.1 干旱胁迫对偏关苜蓿保护酶的影响

干旱胁迫下植物体内产生过量的活性氧，活性氧的氧化作用会对植物造成损伤，植物清除活性氧的能力强弱是抗旱性评价的重要生理指标[28]。SOD是一种广泛存在于植物内的酶，可以清除活性氧自由基从而使植物免受伤害。OLGA等[29]在干旱胁迫下研究不同抗旱性玉米品种的SOD活性差异，结果表明，正常条件下不同抗性品种的SOD活性无明显差异，但在干旱胁迫下抗旱性强的品种SOD活性显著高于抗旱性差的品种，可能的原因是抗旱性强的品种保护酶活性强，SOD相关基因表达量更高。本研究结果表明，偏关苜蓿叶片SOD活性显著高于未胁迫时，与韩刚等[30]研究结果一致，说明干旱胁迫下SOD活性提高有利于偏关苜蓿清除植物体内产生的过量活性氧，增强偏关苜蓿的干旱胁迫耐受性。

POD是酶促反应系统的保护酶，通过防御活性氧来缓解干旱胁迫对植物的伤害。植物感受干旱信号后会调控相关基因表达，多数POD的合成属于诱导表达型，与植物信号转导有关[31]。干旱处理后小麦幼苗根和叶中的POD活性显著增强[32]，与本研究一致，偏关苜蓿在干旱胁迫下POD活性显著升高，发挥其清除活性氧的生理功能。

CAT对H2O2具有较高的亲和力，可以清除植物体内过多的H2O2，将H2O2分解为对细胞无害的水和分子氧。研究表明，干旱胁迫下植物体内CAT活性会发生改变，高粱CAT活性在干旱胁迫下显著增加，证明CAT是高粱重要的保护性酶[33]。也有研究发现，小麦抗旱性与干旱胁迫下CAT基因表达量密切相关，CAT等保护酶活性较高的品种抗旱性更强[34]，与本试验结果相一致，干旱胁迫下偏关苜蓿叶片CAT活性显著高于未胁迫时，说明偏关苜蓿可以增强CAT活性来响应干旱胁迫。

APX是植物体内清除H2O2的关键酶。研究表明，甘薯在干旱胁迫下APX活性增强，H2O2含量减少，从而削弱其对光合的抑制，保护叶绿体[35]。转APX枣树在干旱条件下，APX活性高于非转基因植株，表明转基因枣树耐旱能力强于非转基因植物[36]。本研究结果表明，偏关苜蓿叶片APX活性显著高于未胁迫时，APX通过有效清除活性氧在增强偏关苜蓿的耐旱能力中起一定作用。

3.2 干旱胁迫下偏关苜蓿转录组的差异性

干旱胁迫下，植物细胞通过信号转导来调控抗旱相关基因表达，使自身结构、生理生化以及物质代谢发生一系列变化。利用高通量测序技术可以研究偏关苜蓿干旱响应基因，对测序结果的统计和评估可知本次数据质量较好，本研究共得到79 794条Unigenes，平均长度1 018 bp，所得到的Unigenes数量低于黄花苜蓿的114 347条[37]，而高于中苜一号紫花苜蓿的41 734条[15]，造成这种差异的原因可能是品种和试验环境不同。本研究所有的Unigenes在六大功能注释数据库中进行了注释，其中在NR数据库的注释率最高，与刘佳月[15]对紫花苜蓿的研究结果相同。Unigenes在NR数据库中的注释率高达99.1%，结果发现与偏关苜蓿同源率最高的是高粱，其次是玉米和日本水稻。SSR标记有多态性高、稳定性好和不受环境影响等优点，是遗传变异研究最好的标记之一，目前已在苜蓿种质鉴定、亲缘关系、DNA遗传图谱构建、分子标记育种和遗传多样性等领域广泛应用[38]。本研究识别的7 868个SSR标记可为偏关苜蓿遗传图谱构建和分子育种提供依据。

本研究中GO功能分析结果可知，分子功能显著差异在作用于糖基键的水解酶活性、蛋白质异二聚化活性、O-糖基化合物水解酶活性等；细胞组分显著差异在细胞质囊、微管、DNA蛋白复合物等；生物学过程显著差异在DNA包装、DNA-蛋白质复合体亚基结构、DNA构象变化等。说明偏关苜蓿通过复杂的基因调控网络来响应干旱胁迫。白雪梅[39]对紫花苜蓿的研究结果表明，蛋白激酶基因是紫花苜蓿在响应干旱胁迫中的关键信号传递体，参与磷酸化和去磷酸化，在抗旱中起关键作用。本研究发现，对照组与干旱组在过氧化物酶活性、抗氧化酶活性和蛋白磷酸酶活性相关基因功能有显著差异。其中，抗氧化酶和过氧化物酶的基因差异性表达可以反映偏关苜蓿对干旱的分子响应，抵抗活性氧诱导的氧化应激损伤。另外，有研究表明，蛋白磷酸酶活性相关基因参与了草地早熟禾对干旱胁迫的响应[40]。

KEGG分析发现，本研究中偏关苜蓿差异表达基因显著富集在植物激素信号转导与淀粉和蔗糖代谢通路，这表明植物感受干旱信号后会将信号传递给相关细胞器来引起抗旱基因表达，使淀粉和蔗糖代谢相关途径向有利于渗透调节的方向激活或抑制。一般而言,Ca2+信号传导、植物激素信号、蛋白激酶和转录因子是植物体内重要的信号传递体[41]，在本研究得到的差异表达基因中均有注释。淀粉和蔗糖代谢能够影响植物的生长及其胁迫响应，这表明偏关苜蓿受到干旱胁迫后通过形成一些糖类小分子化合物来进行渗透调节[42]，其原理是干旱胁迫下淀粉和蔗糖可以降解为易溶于水的葡萄糖，同等条件下单糖溶液渗透压大于双糖和多糖溶液，从而避免植物细胞损失更多的水分[43]。

能够显著响应干旱胁迫的基因可以作为苜蓿处在逆境时分子水平的应答特征，同时具有提高苜蓿耐旱和抗旱的潜质。本研究中干旱胁迫下偏关苜蓿蔗糖合成酶基因的下调以及淀粉水解酶基因的上调表达，一方面可以为植物响应干旱胁迫提供更直接的能量，另一方面可以防止细胞内渗透压下降造成的失水。JACOBSEN等[44]研究也发现，随着叶片水势下降，大麦幼苗叶片淀粉酶活性增加。WANG等[45]对茶树在干旱胁迫下转录组的表达进行分析，发现差异表达基因在黄酮类生物合成和苯丙烷类生物合成途径大量富集。本研究中差异表达基因在苯丙烷生物合成中显著富集，其中苯丙氨酸酶基因上调表达，表明了与苯丙烷生物合成途径相关的关键基因在偏关苜蓿响应干旱胁迫中起重要作用。类胡萝卜素相关基因表达的差异也是苜蓿品种抗旱差异的重要原因之一[16]。相关上述研究的差异表达基因有助于提高耐旱性，可作为与偏关苜蓿干旱胁迫相关的潜在抗旱基因。

植物处于干旱胁迫时转录因子起不可或缺的作用，本研究中bHLH、WRKY、ERFs和MYB家族基因在偏关苜蓿响应干旱胁迫过程意义重大。通常认为，WRKY、bHLH、NAC等转录因子家族与干旱胁迫有关[15]。MYB转录因子对干旱胁迫有负调控作用[46]，干旱胁迫诱导NAC转录因子上调[47]。LIU等[48]采用PEG胁迫小麦1、6 h后取叶片进行转录组测序，差异基因主要富集到FAR1、NAC、bZIP、bHLH、AP2/ERF、WRKY、Myb-related和Myb转录因子家族。对杜梨干旱胁迫的转录组分析结果显示，数目最大的转录因子家族是MYB，有72个MYB基因差异表达，其次是bHLH，有57个bHLH基因响应干旱胁迫，45个WRKY基因在干旱的诱导下差异表达，39个NAC基因、15个HSF基因在干旱条件下差异表达[49]。对不抗旱香蕉品种和抗旱香蕉品种进行干旱胁迫后取样测序，结果表明，与香蕉干旱胁迫相关的转录因子大多是MYB、NAC、WRKY、bHLH、bZIP、ERF和HSFs转录因子家族的成员[50]。

植物复杂的信号网络可以使其迅速地感知外界环境变化，调控基因的表达，从而形成应对干旱胁迫的多种机制。WANG等[51]研究了过表达葡萄V vbHLH 1拟南芥与非转基因拟南芥在干旱胁迫下的差异性，结果表明，VvbHLH 1的过表达增加了植株中SOD和POD的活性，H2O2和丙二醛含量显著下降；与脯氨酸生物合成、黄酮生物合成和抗氧化酶相关的基因上调。这说明VvbHLH 1通过提高抗氧化酶活性和增加黄酮类化合物的积累增强了拟南芥对干旱胁迫的适应能力。在本研究中也发现，黄酮和黄酮醇的生物合成相关基因和抗氧化相关基因显著表达，并且挖掘出bHLH转录因子，这表明在偏关苜蓿中可能也存在类似的抗旱过程。有研究发现，在干旱胁迫下过表达甘薯IbWRKY 2的拟南芥中SOD活性较高，并且MDA和H2O2的含量较低，与脱落酸信号转导途径、抗氧化酶和脯氨酸生物合成有关的基因差异显著[52]。这表明WRKY转录因子具有提高转基因拟南芥耐旱性的作用，本研究中偏关苜蓿的WRKY转录因子具有提高苜蓿抗旱性的潜力。有研究表明，苯丙烷类代谢主要受MYB转录因子家族的调控[53]；在本研究中也发现，苯丙烷生物合成以及黄酮和黄酮醇的生物合成途径相关基因显著表达，所挖掘的偏关苜蓿bHLH及MYB转录因子可能通过调控保护酶活性相关基因和注释到上述代谢通路的基因来提高抗旱性。基于偏关苜蓿响应干旱胁迫的差异表达基因，根据功能注释结果可知，偏关苜蓿可能通过调控抗氧化活性和酶调节活性相关的基因以及bHLH等转录因子，使SOD、POD、CAT及APX酶活性上升，减少干旱胁迫造成的损伤。本研究中受干旱胁迫诱导而表达的基因，可在今后苜蓿耐旱育种中作进一步研究。

4 结论

本研究以偏关苜蓿为材料，以正常条件和干旱胁迫下的偏关苜蓿叶片进行保护酶活性测定和转录组分析，结果表明，偏关苜蓿可以通过调节抗氧化酶活性来响应干旱胁迫，转录水平上总共获得56 911条基因序列，识别出7 868个SSR位点。通过差异表达基因的分析得出，Unigene12274_All、Unigene10105_All、Unigene35365_All、CL7158.Contig1_All、CL 11234.Contig1_All、CL 12445.Con tig1_All、Unigene12274_All、CL 4041.Contig1_All以及bHLH 35、ERF25、PCF7和RF2b可作为偏关苜蓿潜在抗旱基因，为后续苜蓿抗旱分子机制研究提供理论依据。