基于KEGG通路分析紫花苜蓿的冷适应过程

徐洪雨， 李钰莹

(1.山西农业大学草业学院， 山西 太谷 030801； 2.农业农村部饲草高效生产模式创新重点实验室， 山西 太谷 030801)

紫花苜蓿(MedicagostaviaL.)是一种重要的多年生饲用豆科植物，具有营养价值高、适应性广等优点。它还是世界上栽培最广泛的作物之一，常作为种质资源用于饲料的育种[1]。我国苜蓿种植区分布在全国14个省，从低纬度到高纬度均有分布。然而，极端的环境条件影响苜蓿的生长和质量，并限制其生存范围，如低温、盐渍和干旱胁迫[2]。在美国中西部、加拿大、北欧及我国东北等地，苜蓿常常无法抵御严酷的冬季条件[3]。2012年，我国内蒙古通辽苜蓿的返青率在50%以下[4]；2012—2020年间，内蒙古阿鲁科尔沁旗优质苜蓿生产基地也多次发生返青率低的问题[5]。

当植物暴露在非冷冻温度下时，可通过调整自身获得抗寒性，从而适应后期的冰冻条件，并能够生存下来。从进化的角度来讲，此过程称之为冷适应或冷驯化[6]。植物冷适应过程，体内发生一系列生理生化变化以及转录因子调控冷诱导基因表达的变化[7]。其中，生理生化变化涉及生长形态[8]、膜结构和细胞骨架等[9]；积累可溶性糖(如蔗糖、棉子糖和海藻糖)和低分子量含氮化合物(脯氨酸、甘氨酸和甜菜碱)，改变膜的组成和性能[10-11]。有些耐寒植物能够产生特异蛋白，通过降低冰点来保护细胞结构[12-13]。在转录因子调控冷诱导基因变化方面，研究最广泛的耐冻途径为CBF(C-repeat Binding Factor)信号通路[14]。它包括三个基因CBF1，CBF2和CBF3，这些基因共享一个保守的AP2 DNA结合区，并通过启动子中特征良好的冷重复/干旱响应(CRT/DRE)调控元件激活冷响应基因(COR)的表达[15]。CBF1和CBF3基因受CBF2的负调控，并对诱导整个CBF调节子和冷适应的提升具有协同累加效应[16-18]。此外，CBF信号通路可提高植物的抗寒性，对拟南芥(Arabidopsisthaliana)[19]、大麦(Hordeumvulgare)[20]、白杨(Populussimonii)[21]和水稻(Oryzasativa)[22]等的研究均已证实此结果。MtCBF3(又名MtDREB1C)的克隆及在蒺藜苜蓿(Medicagotruncatula)中的过表达诱导了冷调节基因(COR)的表达，从而提高转基因植株的抗冻性[23]。冷适应过程非常复杂，几乎涉及到细胞的每一个过程，如何找出调控植物耐寒性的关键过程是一个巨大挑战[18]。

近年来，核糖核酸(RNA)测序分析已成为一种可以生成大量基因序列并构建表达谱的有效方法[24]。抗寒性因植物种类而异，在一定程度上也取决于特定组织或器官对冷冻胁迫的反应。近年来，在应对盐渍和干旱胁迫中，人们对紫花苜蓿进行了多项研究，包括单核苷酸多态性的发现、基因鉴定及转录组的分析[25-26]。在冷胁迫条件下，苜蓿转录组的变化也有学者进行了测定，但仍未能有效筛选出影响抗寒性的关键基因。植物物质代谢过程并不是由单个基因所决定，因此筛选单个基因并通过研究其功能来解析苜蓿的抗寒机制往往效率较低。2020年，苜蓿基因组的发布，对于转录组研究奠定了坚实基础[27]。借助转录组测序技术及KEGG数据库分析，筛选冷适应过程的关键生物通路，进一步分析参与的基因或基因模块，将为解析苜蓿冷适应调节机制提供新的思路。

本研究选择我国种植广泛且抗寒性差异较大的苜蓿品种‘肇东’和‘WL440HQ’作为试验材料，借助高通量测序技术对冷适应处理前后苜蓿的基因表达进行测定。根据转录组测定结果，筛选冷适应处理前后差异表达基因，进一步从代谢的角度分析与苜蓿冷适应过程抗寒性变化相关的代谢过程，以及冷适应后2品种苜蓿抗寒性差异的原因。

1 材料与方法

1.1 试验材料

本研究选择我国北方常用的苜蓿品种‘WL440HQ’(WL，秋眠级6)和‘肇东’(ZD，秋眠级2)作为试验材料。2020年5月，将苜蓿种子播种在沙质试验田内，幼苗期每2日浇水一次，保持土壤湿润。此外，及时控制杂草生长及病虫害的发生。

播种一个月选择生长一致的苜蓿幼苗移栽至PVC管，而后在培养箱里进行试验处理，移栽定植每管4株。PVC管的规格为直径10 cm，高15 cm。移栽时，管内的填充物是由土壤、珍珠岩和蛭石按重量比100∶30∶55混合而成，其田间持水量为151%([湿重-干重]/干重)。移栽后，先将材料放在培养箱(中国广东泰宏器械)内进行2周的恢复生长，期间培养箱内的温度设置为24℃/20℃(昼/夜)、光周期为14 h/10 h(昼/夜)和光照为300～400 μmol·m-2·s-1。恢复生长后，进行2个阶段的温度控制试验，包括适应性生长和冷适应处理。本研究试验处理设置3次生物学重复，共需苜蓿苗48株，其中包括2品种×2阶段×3生物学重复×每管4株。

1.2 试验处理

恢复生长后进行适应性生长和冷适应处理(详见图1)，期间每日补水控制PVC管内的基质含水量为田间持水量的75%，且每日随机变换PVC管的位置。阶段1为适应性生长，18℃(昼/夜)环境下培养1周后在S1点取样。所得样品采用“品种+试验阶段_生物学重复”的方式进行命名，S1点所得样品分别命名为WL1_1，WL1_2，WL1_3，ZD1_1，ZD1_2和ZD1_3。阶段2为冷适应处理，先在2℃(昼/夜)下培养1周，而后将温度降低到-2℃(昼/夜)再培养1周，在S2点再次取样，所得样品分别命名为WL2_1，WL2_2，WL2_3，ZD2_1，ZD2_2和ZD2_3。在2个试验阶段，光周期设置为12 h/12 h(昼/夜)，光强100～150 μmol·m-2·s-1。

图1 试验处理过程Fig.1 The experimental process of two phases

1.3 取样

取样时，每个品种随机选择3管，每PVC管里的4株材料上摘取的叶片作为1个生物学重复的样品。取样完成后，叶片样品分别用做低温半致死温度(LT50)和转录组的测定。其中，新鲜样品通过测定细胞液相对渗透率，计算LT50以评价样品的抗寒性；另一份样品液氮速冻后放于-80℃备用，用于RNA测序。

1.4 低温半致死温度

本研究根据LT50评价样品的抗寒性，即细胞内离子相对渗透性达到50%时所对应的处理温度，其测定方法参考Xu[28]。离体冷冻测试前，将新鲜的叶片放于2 mL的离心管内，4℃下放置2 h，然后放在冰上做短暂的保存。冷冻测试时，在ZX-5C恒温循环器(智信，中国上海)内通过一系列逐渐降低的温度对样品进行离体冷冻。对于取样点S1的样品，冷冻温度分别为-4℃，-6℃，-8℃，-10℃，-12℃，-14℃，-16℃，-18℃和-20℃；对于取样点S2的样品，冷冻温度设置为-5℃，-8℃，-11℃，-14℃，-17℃，-20℃，-23℃，-26℃和-29℃。在每一温度下冷冻1 h后，取出一份装有相应样品的离心管并放置在冰里过夜保存。次日，将装有样品的离心管先在4℃下缓冲2 h。而后，将样品转移到15 mL的离心管内，加入5 mL的去离子水后在摇床HZQ-A(恒瑞，江苏)上摇12 h。随后，使用电导率仪FE38(Mettler Toledo，上海)测定样品的电导率并记为EL1。测定完EL1后，将样品在120℃下蒸煮30 min，测定电导率记为EL2。细胞液相对渗透率大小及样品LT50大小分别根据公式(1)和公式(2)进行计算。在公式(1)中，EL为去离子水的电导率，在公式(2)中，x代表冷冻温度，y代表相对电导率，A和B是常数。在抗寒性比较分析中，采用最小显著差数法进行差异显著性判断，差异显著性水平P< 0.05。

相对电导率(%)=(EL1-EL)/(EL2-EL)×100

(1)

y(%)=A/(1+B×e-kx)×100

(2)

1.5 RNA的提取、文库的建立及转录组测序

本研究使用Trizol试剂(Invitrogen，USA)并根据说明书提取叶片的总RNA。使用AMPure XP系统(Beckman Coulter，Austrilia)对RNA进行纯化，并用Agilent Bioanalyzer 2100系统(Agilent Technologies，USA)对RNA质量进行评估。RNA测序工作用生物标记技术完成(北京，中国)，使用NEBNext®UltraTMRNA Library Prep Kit for Illumina®(NEB，USA)构建RNA-Seq建库。在总RNA中使用带有磁珠的poly-T oligo对mRNA进行纯化。用逆转录酶将打碎后的RNA片段合成第一链cDNA，然后使用DNA聚合酶I和RNase H合成第二条cDNA链。将这些片段连接到测序适配器上，在Illumina HiSeq 2500平台上对文库制备物进行测序，生成平均长度为150 bp的原始数据，本研究原始序列存入NCBI-SRA数据库(PRJNA814482)。

1.6 序列比对及功能注释

测序完成后，去除接头和低质量序列(含有poly-N及碱基质量值小于Q20的Reads)，对原始数据进行过滤，从而获得Clean Data。本研究最终采用碱基质量值≥Q30的Reads。而后，使用HISAT2软件将Clean Reads参考紫花苜蓿基因组进行比对[27]。比对分析完成后利用String Tie软件将比对上的Reads进行重新组装和定量。根据重新组装结果，通过以下数据库对基因功能做进一步注释分析，其中包括Nr(NCBI冗余蛋白质序列数据库)、Nt(NCBI冗余核苷酸序列数据库)、Pfam(蛋白质家族数据库)、KOG/COG (同源蛋白质集群数据库)、Swiss-Prot(手工注释和审查蛋白质序列数据库)、KO(全基因组及代谢途径数据库)和GO(基因本体数据库)。

1.7 差异表达基因

在转录组分析中，本研究从阶段间比较的角度横向分析苜蓿对冷适应过程的响应，从阶段内比较的角度纵向分析2品种苜蓿响应冷适应低温的差异性。其中，阶段间比较包括ZD2 vs ZD1和WL2 vs WL1对比组，阶段内比较包括ZD1 vs WL1和ZD2 vs WL2对比组。其中，各样品对比组差异表达基因的筛选方法如下：使用Bowite 2软件将Clean Reads与组装序列进行比对，将比对上的Reads进行计数，并使用RSEM方法估计丰度，以获得FPKM、TPM和预期计数；基于负二项分布模型，用DESeq 2软件鉴别样本间的差异表达基因。为控制错误发现率，对得到的P值用Benjamini和Hochberg方法进行调整。当调整后的P<0.01时，DESeq 2软件便将此基因选定为差异表达基因。随后，在对差异表达基因进行KEGG通路分析时，根据某生物通路富集差异基因的数目进行排序，以此评价各生物通路在试验中的活跃程度。

2 结果与分析

2.1 抗寒性

图2为冷适应前后2品种苜蓿LT50的大小，LT50越低，抗寒性越强。与S1点相比，经过冷适应处理后，在S2点WL和ZD的抗寒性均显著提高(P<0.05)，其中ZD抗寒性提高的程度较大。冷适应结束后，WL的LT50达到-15.72℃，而ZD的LT50达到-24.70℃。冷适应结束后，ZD的抗寒性显著高于WL(P<0.05)。

2.2 转录组测序

根据样本间的相关性热图3判断，WL1，WL2，ZD1和ZD2四组样品内的生物学重复的重复性均较好。在WL2组内的生物学重复中，只有WL2_1与另外2样品的重复略差，但相关性系数仍达到0.70。本研究完成了12个样品的转录组测序，共获得38.2 Gb Clean Data，其中各样品Clean Data中Q30碱基百分比≥93%。各组样品Clean Reads与参考基因组进行序列比对，比对效率为90.75%～94.08%。基于比对结果，进行基因表达量分析，最终根据基因在不同样品中的表达量识别筛选差异表达基因，并进一步进行KEGG富集分析。

图2 冷适应处理前后苜蓿的低温半致死温度Fig.2 Semilethal temperature of alfalfa before and after cold acclimation注：“*”表示同一时间不同品种间差异显著(P<0.05)，“ns”表示品种间无显著差异Note:“*” means significant difference between two alfalfa at the same time point at the 0.05 level,and “ns” means no significant difference between two alfalfa

2.3 同一苜蓿品种不同时间点间转录表达的对比分析

根据本研究基因差异表达热图结果判断，ZD2 vs ZD1和WL2 vs WL1的基因表达可以被明显分开。经过冷适应处理后，ZD2 vs ZD1中的差异基因数共有26 986个，包含14 699个上调基因和12 287个下调基因；WL2 vs WL1中的差异基因数共有25 018个，包含13 123个上调基因和11 895个下调基因。在冷适应过程中，ZD和WL共同上调的基因有9 416个，共同下调的基因有7 890个(图4)。此外，在ZD中为上调基因而在WL中却为下调的基因有12个；在ZD中为下调基因而在WL中为上调的基因共有13个。

在ZD2 vs ZD1和WL2 vs WL1对比组中，借助KEGG数据库分析差异基因参与的生物通路。图5展示了富集差异基因数目最多的前30条生物通路，即冷适应过程中2品种苜蓿体内较活跃的通路。经对比，2品种苜蓿体内活跃度排名前30的生物通路名称及排序无明显差别，且这些生物通路均涉及到“细胞的信号转导”和“物质代谢(包括糖、蛋白质、脂、氨基酸、脂肪酸)”等过程。这些生物通路与苜蓿响应冷适应过程中的低温环境有关，如植物激素信号转导、内质网中的蛋白质加工、淀粉和蔗糖代谢、植物MAPK信号通路、糖酵解和糖异生、甘油磷脂代谢、甘油脂代谢和泛素介导的蛋白质水解作用等过程。

图3 试验样本之间的相关性热图Fig.3 Heat map of correlation between experimental samples

图4 阶段间比较中差异基因韦恩图Fig.4 Venn diagrams of differentially expressed genes in comparisons within phases

2.4 同一时间点不同苜蓿品种间转录表达的对比分析

根据本研究基因差异表达热图结果判断，ZD1 vs WL1和ZD2 vs WL2的基因表达可以被明显分开。据图6分析，冷适应处理前，ZD1 vs WL1中共筛选差异基因1 490个，包含 635个上调基因和855个下调基因。冷适应处理后，ZD2 vs WL2中共筛选差异基因1 968个，包含946个上调基因和1 022个下调基因。其中，ZD2 vs WL2与ZD1 vs WL1共有99个相同的上调基因和154个相同的下调基因。此外，在ZD2 vs WL2样品对比组中为下调基因，而在ZD1 vs WL1中为上调基因的数量只有1个(MS.gene26754，GO:0016758 transferase activity,transferring hexosyl groups)。

本研究在同一试验阶段内进行了差异基因的筛选，进一步又对这些差异基因进行了KEGG富集分析。根据KEGG分析后所富集的差异基因数目，对富集的生物通路进行了排序。图7展示了与ZD1 vs WL1相比，在ZD2 vs WL2中活跃度排序位次有所提高的生物通路。冷适应处理前后，在ZD1 vs WL1和ZD2 vs WL2中差异表达基因的变化，从而导致图7中统计的生物通路在2品种苜蓿中的活跃度有所改变，这些活跃度提高的生物通路，可能是导致冷适应后2品种苜蓿抗寒性产生差异的原因。

图5 阶段间比较中差异基因数最多的30条生物通路Fig.5 The 30 metabolic pathways with the largest number of differentially expressed genes in comparison between phases

图6 阶段内比较中差异基因韦恩图Fig.6 Venn diagrams of differentially expressed genes in comparison within phases注：位次变化=在ZD2 vs WL2比较中的位次-在ZD1 vs WL1比较中的位次。此图只统计了位次有所提高的生物通路，数字表示位次上升的级数Note:Rank change = Rank in ZD2 vs WL2-Rank in ZD1 vs WL1. This figure only counts the biological pathways that the rank increased,and the numbers indicate the progression of rank

3 讨论

本研究在冷适应前后，分别对比了‘肇东’和‘WL440HQ’2品种苜蓿的抗寒性，结果表明经过冷适应处理后，2品种苜蓿抗寒性均有显著性提高，但提高的程度有所不同，冷适应处理后品种间抗寒性存在显著性差异。为了解冷适应过程温度变化对苜蓿代谢的影响，本研究分别对同一品种苜蓿冷适应前后进行了转录表达分析，筛选2品种苜蓿在冷适应前后的差异表达基因，并进一步进行KEGG分析。结果如图5所示，对于‘肇东’和‘WL440HQ’而言，冷适应处理前后差异表达基因所参与的通路并无明显差异，这些通路涉及“细胞的信号转导”和“物质代谢(包括糖、蛋白质、脂、氨基酸、脂肪酸)”等过程。此外，配合试验处理过程分析，这些通路不仅是苜蓿对冷适应低温环境的响应，还可能与苜蓿抗寒性的提高有关，具体包括植物激素信号转导、内质网中的蛋白质加工、淀粉和蔗糖代谢、植物MAPK信号通路、糖酵解和糖质新生、甘油磷脂代谢、甘油脂代谢和泛素介导的蛋白质水解作用等。

图7 与ZD1 vs WL1相比，在ZD2 vs WL2中活跃度排序位次提高的生物通路Fig.7 Metabolic pathways with increased activity rank in ZD2 vs WL2 compared to ZD1 vs WL1

在本研究中，植物激素信号转导参与了苜蓿对冷适应低温环境的响应。有研究表明，植物激素信号介导非生物胁迫响应存在多种机制，包括ABA、水杨酸、生长素、乙烯和茉莉酸[29]。激素作为信号分子，在调节基因表达中发挥着关键作用。如ABA是一种苜蓿体内重要的应激激素，它能是苜蓿对寒冷及冰冻胁迫做出响应，可能与ABA依赖通路中PP2C基因的表达有关[30]。

本研究中，2品种苜蓿冷适应前后部分差异表达基因均显著富集于泛素介导的蛋白质水解系统(见图5)，是苜蓿对低温响应的代谢通路之一。泛素介导的蛋白质水解作用参与许多生物学过程，包括植物的生长以及对新环境条件的适应[31]。3种酶参与靶蛋白的泛素化，可解决环境胁迫导致的蛋白质异常积累，包括E1泛素激活酶、E2泛素缀合酶和E3泛素连接酶[32]。E3连接酶是决定靶蛋白特异性的关键酶，据其基序可分为不同的家族，与E6-AP羧基端、Ring/U-box和后期酶促复合体、Skp1-Cullin-F-box复合体同源。F-box蛋白是Skp1-Cullin-F-box复合体的亚基，参与逆境响应，可由不同的环境胁迫条件而触发，如低温、干旱或盐碱胁迫[33-35]。近年来，在植物中发现了许多涉及激素信号转导、转录调控和细胞周期转换的F-box蛋白[36-38]，通过过表达F-box蛋白可以提高转基因植物的抗氧化能力和抗旱能力[39]。在苜蓿对低温胁迫响应的研究中，学者们也发现了冷冻处理可诱导F-box蛋白的表达，在低温胁迫中发挥重要作用[30]。据前人的研究基础及本研究结果判断，泛素介导的蛋白水解可能是调控苜蓿抗寒性的重要途径之一。

在本研究中，2品种苜蓿在冷适应处理前后，部分差异表达基因均可显著富集于淀粉和蔗糖代谢通路。可溶性糖可提高细胞渗透压、减少胞间水分含量、降低胞间冰形成的数量[40]，还可降低胞液的冰点[41]，是植物应对低温胁迫的重要响应。即使在结冰脱水的情况下，可溶性糖仍可以保护细胞的结构[42]和蛋白质的结构[40]。有学者认为蔗糖是苜蓿应对低温胁迫的重要响应，甚至将其含量用来评估品种的抗寒性[43]；但也有学者认为棉子糖和水苏糖才是保护苜蓿不受冻害最有效的糖类物质[44]。也有研究表明，一分子的蔗糖与一分子的肌醇半乳糖苷结合，在棉子糖合成酶的作用下可生成一分子的棉子糖。棉子糖在水苏糖合成酶的作用下，再结合一分子的肌醇半乳糖苷有可生成一分子的水苏糖。因此，蔗糖分子可称为棉子糖和水苏糖合成的原料[45]，而蔗糖又可通过淀粉水解得到[46]。由此可见，淀粉和蔗糖的代谢通路是苜蓿对低温环境的重要响应。

感知低温信号，调控冷诱导基因的表达，是植物应对低温胁迫的重要路径，MAP级联反应就是其中一种重要的感知低温信号和调控基因表达的途径。低温胁迫引起细胞Ca2+浓度变化，从而激活钙调蛋白受体激酶CRLK1/2，引起MAPK级联反应，从而诱导冷调节基因COR表达量的变化[47]。冷调节基因的表达，进一步调控细胞生理生化过程以应对胁迫环境。在植物冷适应过程中，细胞膜组分的变化，可提高细胞膜流动性从而提高植物的抗寒能力[48]。SFR2是一种对拟南芥耐冻性至关重要的基因，它可编码叶绿体外膜的半乳糖脂质重塑酶，通过一系列加工转移过程改变膜的组分，最终起到稳定细胞膜的作用[49]。在本研究中，冷适应前后差异表达基因涉及甘油磷脂代谢和甘油脂代谢，二者为细胞膜重要组分，2代谢过程可能与冷适应过程温度环境的降低以及苜蓿抗寒性的提高有关。

为进一步了解2品种苜蓿对冷适应低温环境响应的差异，本研究在冷适应前后的2时间点分别进行了品种之间转录表达的对比分析，筛选出2苜蓿品种之间在冷适应前后的差异表达基因，并进一步进行生物通路分析。结果发现，与冷适应处理前相比，冷适应处理后共有66条生物通路的活跃程度有所提高(见图7)。其中，这些通路主要涉及代谢，包括脂肪酸、氨基酸、脂质、聚糖、萜类和多酮类化合物等物质的代谢和能量代谢。在以往的研究中，已表明脂肪酸[48]、氨基酸[44]、脂质[50]、聚糖[51]、萜类和多酮类化合物[52]等物质对苜蓿抗寒性的积极影响，同样也有研究表明了能量代谢与植物抗寒性的关系[53]。因此，冷适应处理后2品种苜蓿抗寒性的差异受多种代谢通路共同影响，特别是活跃度变化比较明显的氨基酸、脂肪酸类物质的代谢。对于这些代谢可进一步开展更深入的研究，如对代谢通路所涉及的基因模块及基因的功能的研究。

4 结论

经过冷适应处理，‘肇东’和‘WL440HQ’苜蓿的抗寒性均有所增强，‘肇东’苜蓿提高的幅度较大。“细胞的信号转导”和“物质代谢(包括糖、蛋白质、脂、氨基酸和脂肪酸)”等相关生物通路，不仅是苜蓿对冷适应降温环境的响应，还可能与此过程中苜蓿抗寒性的提高有关，具体包括植物激素信号转导、内质网中的蛋白质加工、淀粉和蔗糖代谢、植物MAPK信号通路、糖酵解和糖质新生、甘油磷脂代谢、甘油脂代谢和泛素介导的蛋白质水解作用等。造成冷适应处理后2品种苜蓿抗寒性出现差异是多种生物通路共同作用的结果，主要涉及“物质代谢(包括脂肪酸、氨基酸、脂质、聚糖、萜类和多酮类化合物等)”和“能量代谢”等相关生物通路，特别是脂肪酸和氨基酸类物质的代谢。