转录组及代谢组联合解析胡麻根部对盐胁迫的响应机制

赵 玮，赵 利，张建平，齐燕妮，王利民，谢亚萍，李闻娟，党 照，袁明璐，张艳萍

（1. 甘肃省农业科学院作物研究所, 甘肃 兰州 730070；2. 甘肃省农业科学院生物技术研究所, 甘肃 兰州 730070）

胡麻(Linure usitatissmum)即油用亚麻，属亚麻科亚麻属一年生草本植物[1]。胡麻油中含有45%～65%的α-亚麻酸，是ω-3 型不饱和脂肪酸的重要来源[2]。α-亚麻酸具有促进智力发育、强身健脑、降血脂、预防血栓等作用[3]，且α-亚麻酸及其衍生物二十碳五烯酸(eicosapntemacnioc acid, EPA)和二十二碳六烯酸(docose hexaenoie acid, DHA)具有预防肿瘤发生和抑制肿瘤细胞增殖的作用[4]。中国胡麻收获面积约35 万hm2，总产量近40 万t，均位居世界前列，我国胡麻加工和消费水平在国际上已处于绝对的优势地位[5]。近年来随着胡麻保健功能食品的飞速发展，其种植效益逐年提高，市场需求越来越大[5]。全世界盐渍土面积约10 亿hm2，在全球干旱和半干旱地区约有50%的灌溉土地受到盐碱化的影响，中国土地盐渍化面积大约占全球的1/10[6-7]。受到盐胁迫时，植物因吸收过量盐离子抑制其吸收部分营养元素而产生盐害，因此土壤盐渍化已成为现代农业发展的重要制约因素之一[8-9]。提高作物的耐盐性，培育适合在盐碱地生长的植物品种是有效利用盐渍化土地的主要方法之一[10]。我国胡麻主要分布在甘肃、内蒙古、宁夏、山西等干旱半干旱省区，灌溉水中一般都含有一定的易溶盐分，长期灌溉必然使土壤盐含量增加[11]。目前，胡麻响应盐胁迫的关键基因及分子机制尚不明确，制约了胡麻的耐盐高产遗传改良进程，威胁着胡麻产业的可持续发展[12]。

随着高通量测序技术快速发展，作物育种研究逐渐从常规育种转为多组学技术相结合的方式，通过对重要农艺性状、遗传基因的解析和注释，为作物育种工作提供了更有效的方案[13]。代谢组+转录组联合分析方法，能较为全面地实现从“因”到“果”深入探究作物的抗逆机理，相互验证作用更加明显[14]。本研究通过盐胁迫处理下的胡麻根部“广靶代谢组+转录组”差异表达基因和积累代谢物的共表达分析，锁定调控不同抗盐机制的关键代谢通路，找出关键调控因子及功能基因[15]。结合功能基因分析、代谢通路富集和分子互作，探明代谢物与耐盐调控机制之间的关联性，最终实现对胡麻耐盐代谢通路和调控机制的解析，筛选出核心代谢通路、代谢产物和功能基因，为提高胡麻耐盐性奠定基础。

1 材料与方法

1.1 试验材料

试验材料来自甘肃省农业科学院作物研究所筛选的耐盐材料R40 和敏感材料R24，水培20 d (温度为25 ℃、湿度75%、光照12 h、光照强度1 500 lx)，待生长出8～10 片叶后用150 mmol·L−1NaCl 溶液进行胁迫处理，分别于1、2、4、8、12 h 取材料的根部组成混合样，以0 h 取样为对照，采样设3 次重复，每个重复采样3 g。采样材料编码如表1 所列。

表1 试验材料编码Table 1 Code notations used for the experimental materials

1.2 样品检测方法

1.2.1 代谢组检测

样品提取流程：将真空冷冻干燥保存的生物样品利用研磨仪(MM 400, Retsch)研磨至粉末状；称取100 mg 的粉末溶解于1.0 mL 提取液中，4 ℃冰箱过夜，期间涡旋3 次；10 000 r·min−1离心10 min后，吸取上清液用0.22 μm 微孔滤膜过滤，后用于液相色谱串联质谱(LC-MS/MS)分析。

色谱质谱采集条件：仪器包括超高效液相色谱(UPLC)和串联质谱(MS/MS)。液相条件包括，1) 色谱柱：Waters, ACQUITY UPLC HSS, T3, C18, 1.8 μm,2.1 mm × 100 mm。2) 流动相：水相为超纯水(加入0.04%的乙酸)，有机相为乙腈(加入0.04%的乙酸)；洗脱梯度：0 min，水/乙腈 (V/V) 为95 ： 5；11.0 min，为 5 ： 95； 12.0 min， 为 5 ： 95； 12.1 min， 为 95 ： 5；15.0 min，为95 ： 5。3) 流速0.4 mL·min−1；柱温40 ℃；进样量22 μL。质谱条件：500 ℃电喷雾离子源(ESI)，5 500 V 质谱电压，帘气(CUR) 25 psi，碰撞诱导电离(CAD)参数设置为高。

试剂和标准品：试验主要试剂有乙腈(Merck)、乙醇(Merck)、 甲醇(Merck)和标准品(BioBioPha/Sigma-Aldrich)，标准品二甲基亚砜(DMSO)或甲醇作为溶剂溶解后，−20 ℃保存，质谱分析前用70%甲醇稀释成合适浓度。

用fold change 值和OPLS-DA 模型的VIP 值相结合的方法筛选代谢物。筛选标准：fold change ＞ 2和fold change ＜ 0.5 的代谢物为具有显著差异。

1.2.2 转录组检测

RNA 样品检测：琼脂糖凝胶电泳方法分析RNA的完整性及是否存在DNA 污染；分光光度计(Nano Photometer)检测RNA 纯度(OD260/280和OD260/230)；荧光计(Qubit 2.0)高精确度测量RNA 浓度；生物分析仪(Agilent 2100)精确检测RNA 完整性。

文库构建：通过磁珠Oligo (dT)与mRNA 的ployA 尾结合，富集真核生物的mRNA，然后将mRNA随机打断；以片段化的mRNA 为模版，使用6 碱基随机引物，在逆转录酶体系(M-MuLV)中合成cDNA第1 条链；随后用RNaseH 降解RNA 链，并在DNA polymeraseⅠ体系下，以dNTPs 为原料合成cDNA 第2条链。纯化后的双链cDNA 经过末端修复、加A 尾并连接测序接头；用AMPure XP beads 筛选200 bp左右的cDNA，进行PCR 扩增并再次使用AMPure XP beads 纯化PCR 产物，最终获得文库。

文库质检：使用荧光计定量(Qubit 2.0)把文库稀释至1.5 ng·μL−1后，使用生物分析仪(Agilent 2100)对文库的插入片段进行检测。将长度符合预期的插入片段用PCR 实时荧光定量(qRT-PCR)，对文库有效浓度进行准确定量(文库有效浓度高于2 nmol·L−1)。库检合格后用Illumina HiSeq 平台进行测序。差异基因的筛选条件为|log2Fold Change| ≥1，且错误发现率(false discovery rate, FDR) ＜ 0.05。

2 结果与分析

2.1 代谢物检测结果与分析

通过检测共筛选到具有显著差异的代谢物741个，其中下调409 个，上调332 个。R40 和R24 两个样本不处理条件下(CK1d vs CK2d)检测出117 个差异代谢物，其中下调66 个，上调51 个；处理后(STd vs SSd)检测出281 个差异代谢物，其中71 个下调，210 上调；R40 处理前后(CK1d vs STd)检测出246个差异代谢物，其中下调220 个，上调 26 个；R24 处理前后(CK2d vs SSd)检测出97 个差异代谢物，其中下调52 个，上调45 个。结果显示，盐胁迫处理后两份样本间筛选的差异代谢物最多，而且上调数是下调数的近3 倍，而敏感型材料R24 处理前后筛选的差异代谢物最少，耐盐型材料R40 处理前后筛选的代谢物数量是敏感型材料R24 的2.5 倍(表2)。

表2 代谢物检测结果Table 2 Metabolite testing results

筛选出差异代谢物排名前十的为有机酸及衍生物115 个，黄酮93 个，氨基酸及衍生物和核苷酸及衍生物均为80 个，脂质75 个，苯丙素74 个，生物碱29 个，黄酮类25 个，维生素及衍生物24 个，黄酮醇19个。除了氨基酸及衍生物和核苷酸及衍生物在R40处理前后(CK1d vs STd)差异数最高外，其他代谢物均在两份材料处理后(STd vs SSd)差异数最多。黄酮数量最多(STd vs SSd)与最少(CK2d vs SSd)的差距为9.5 倍，脂质最多(STd vs SSd)与最少(CK1d vs CK2d)的差距为8.3 倍，黄酮醇最多(STd vs SSd)与最少(CK2d vs SSd)的差距为7 倍(表3)。

表3 排名前十的代谢物结果分析Table 3 Analysis of metabolite results in the top 10

2.2 转录组检测与分析

2.2.1 差异表达基因筛选与注释

本研究共筛选到差异基因30 233 个，其中上调15 827 个，下调14 406 个。敏感型材料R24 处理前后(CK2z vs SSz)筛选的差异基因最多，达到13 171个，下调6 088 个，上调7 083 个，其次是耐盐型材料R40 处理前后(CK1z vs STz)，差异基因达到12 948个，下调6 334 个，上调6 614 个；两份对照材料(CK1z vs CK2z)差异基因数2 720 个，下调1 378 个，上调1 342 个；两份材料处理后(STz vs SSz)差异基因数1 394 个，下调606 个，上调788 个(表4)。

表4 差异基因数量统计Table 4 Differentially expressed genes in different samples

本研究共筛选到有注释结果的转录因子有13 015个，其中细胞生成过程和信号功能5 007 个、信息存储与处理功能2 567 个以及新陈代谢功能5 441 个。根据功能注释结果，差异基因数前十依次为：翻译修复以及蛋白质转换1 779 个，信号转导机制1 545个，碳水化合物的运输和代谢1 222 个，次生代谢物的生物合成、运输和分解代谢1 121 个，转录877个，核糖体结构和生物发生837 个，氨基酸的运输和代谢793 个，能源生产与转换685 个，脂质运输和代谢685 个，无机离子的转运与代谢582 个。差异基因最少的功能分别是细胞外结构和细胞核心结构，分别有85 和36 个(表5)。

表5 差异基因数及功能注释Table 5 Function annotation and gene dosage of the differentially expression genes

续表5Table 5 (Continued)

2.2.2 差异表达基因KEGG 通路

本研究通过差异表达基因KEGG 通路，筛选出差异基因KO 节点24 个，其中上调基因的KO 节点有9 个，下调差异基因的KO 节点有5 个，同时包含上下调基因的KO 节点有10 个。结果显示，调控抗盐代谢途径的重要基因主要以参与调控糖类代谢为主，包括蔗糖、葡萄糖、葡聚糖、葡糖苷、纤维二糖、纤维糊精、腺苷二磷酸葡萄糖、尿苷二磷酸葡萄糖、海藻糖、麦芽糖糊精和麦芽糖等。

2.3 代谢组和转录组联合分析

2.3.1 代谢物KEGG 联合通路分析

为了明确基因与代谢物之间的关系，将盐胁迫下胡麻根部筛选的差异代谢物与转录组差异基因，同时映射到13 条KEGG 通路图上并进行联合分析。主要代谢物包括：生长素(auxin, IAA)、脱落酸(abscisic acid, ABA)、乙烯(ethylene)、菜籽类固醇(brassinosteroid, BR)、茉莉酸(jasmonic acid, JA)、细胞分裂素(cytokinin, CTK)和异亮氨酸(l-isoleucine,JA-Ile)等。主要代谢通路包括：色氨酸-生长素通路、玉米素生物合成-细胞分裂素通路、二萜生物合成-赤霉素通路、类胡萝卜素-脱落酸通路、半胱氨酸和蛋氨酸的代谢-乙烯通路、菜籽类固醇生物合成-菜籽类固醇通路、α-亚麻酸-茉莉酸通路、苯丙氨酸生物合成-水杨酸通路。其中菜籽类固醇和茉莉酸代谢通路上调，生长素、脱落酸、乙烯等代谢通路下调。在色氨酸-生长素通路中发现GH3基因，其在植物生长素信号途径、光信号途径以及植物抗逆防卫反应中起重要作用。

2.3.2 差异代谢物和差异基因KEGG 富集分析

根据本研究的差异代谢物分析结果，结合转录组差异基因分析结果，将相同分组的差异基因及差异代谢物同时映射到KEGG 通路图上，差异基因和差异代谢物的富集P值，用-logP来表示，结果显示：差异基因P＜ 0.05 的代谢通路一共有46 条，其中对比组CK2d vs SSd 有20 条，CK1d vs CK2d 有4 条，STd vs SSd 有4 条，CK1d vs STd 有18 条；差异基因P＜ 0.01 的代谢通路24 条，其中对比组CK2d vs SSd 有8 条，CK1d vs CK2d 有3 条，STd vs SSd 有2 条，CK1d vs STd 有11 条。差异代谢物P＜ 0.05 的代谢通路一共有27 条，其中对比组CK2d vs SSd 有22 条，STd vs SSd 有1 条，CK1d vs STd 有4 条；差异代谢物P＜ 0.01 的代谢通路11 条，其中对比组CK2d vs SSd 有9 条，CK1d vs STd 有2 条。基因富集数量较多的代谢途径分别有代谢途径相关、抗生素生物合成、次生代谢产物的合成、苯丙素的生物合成、淀粉和蔗糖代谢等通路；代谢物富集数量较多的代谢途径分别有代谢途径相关、次生代谢产物的合成、抗生素生物合成、苯丙素的生物合成、ABC转运蛋白、类黄酮生物合成等通路(表6)。

表6 P ＜ 0.05 水平的代谢物及其关联基因Table 6 Metabolites and their correlated genes at the 0.05 level

续表6 (1)Table 6 (Continued)

续表6 (2)Table 6 (Continued)

续表6 (3)Table 6 (Continued)

图2 CCA 图Figure 2 CCA maps

2.3.3 差异基因和代谢物相关性网络图

对每组差异分组检测到的基因和代谢物进行相关性分析，统计皮尔森相关系数大于0.8 的基因和代谢物的差异倍数情况，通过网络图来表示基因和代谢物之间的相关性(图1)。

结果显示基因与代谢物差异表达模式呈正相关的有207 408 对，负相关的有430 949 对。从代谢物和基因相关性网络图可以看出，Lus_GLEAN_10025781、Lus_GLEAN_10045022、Lus_GLEAN_10038907、Lus_GLEAN_10028484、Lus_GLEAN_10017173、Lus_GLE AN_10022204、Lus_GLEAN_10022072、Lus_GLEAN_10011859共8 个基因，与N-乙酰-L-谷氨酸、烟酸核糖核苷、L-(-)-酪氨酸、7-甲氧基香豆素、L-色氨酸5 个代谢物具有明显的相关性，每一个代谢物与一个或多个基因相关联(图1)。

图1 代谢物与转录组相关性网络图Figure 1 The network between the metabolites and the transcriptome

2.3.4 差异基因和代谢物CCA 分析

典型相关性分析(canonical correlation analysis,CCA)，是利用综合变量对之间的相关关系来反映两组指标之间的整体相关性的多元统计分析方法。对相关性网络图中的差异基因及差异代谢物进行CCA 分析(图2)，结果表明，CCA 图中以十字区分出4 个区域，在同一个区域内，距原点越远关联性越高。其中代谢物L-谷氨酸与基因Lus_GLEAN_10025781、代谢物2-异丙基苹果酸与基因Lus_GLEAN_10022072、代谢物L-(-)-酪氨酸与基因Lus_GLEAN_10017173存在较高的关联。

3 讨论

近年来，转录组学和代谢组学联合分析方法在植物逆境胁迫研究中得到了广泛的应用[16-17]，如对盐胁迫处理的葡萄(Vitis vinifera)叶进行转录组分析发现，差异表达基因主要富集在碳水化合物代谢、氨基酸代谢、次生代谢产物的合成和脂类代谢等代谢通路[18]。本研究利用代谢组和转录组对NaCL 胁迫下胡麻耐盐型材料R40 和敏感型材料R24 根部代谢物含量和基因表达进行分析，结果表明，胡麻根部的差异代谢物和与之对应的功能转录因子富集性高度一致。富集前十的代谢物几乎都有相对应的功能转录因子富集，如黄酮和黄酮醇的生物合成、苯丙素的生物合成、氨基糖和核苷酸糖的代谢、花青素生物合成、类黄酮生物合成等均有对应的转录因子参与其中。本研究发现黄酮代谢物的数量最多(STd_vs_SSd)与最少(CK2d_vs_SSd)的差距达到9.5 倍，说明盐胁迫下胡麻根部黄酮类代谢物的富集较高，而敏感型材料在盐胁迫下黄酮的代谢积累最少，其主要原因是盐胁迫导致大量有害的活性氧产生，类黄酮主要功能是增强植物对活性氧的清除能力从而提高其耐盐性[19]，所以敏感型材料黄酮积累少也是其耐盐差的原因之一。有研究也同样证实类黄酮生物合成途径在高粱(Sorghum bicolor)抗盐胁迫中起到了重要作用[20]。另外氨基酸等有机物富集量较大，主要作用是为了增强细胞的渗透压，从而有效降低盐碱环境对植物的伤害。差异基因前十的功能注释结果，从宏观上认识了盐胁迫下胡麻的基因功能分布特征，即转录因子功能在盐胁迫下主要以新陈代谢功能和细胞生成为主，信息存储和处理的功能其次；差异基因功能注释显示，参与翻译修复及蛋白转换、信号转导机制以及碳水化合物、次生代谢物的生物合成运输和代谢的基因较多，而其他功能的基因数量也相对平衡；差异基因功能注释最少的两个分别是细胞外结构和细胞核心结构，可见盐胁迫下基因参与调控的主要是代谢物的修复、运输以及信号传导等功能，对细胞结构的调控相对较少。将差异基因及差异代谢物同时映射到13 条KEGG 通路图上，结果显示，基因与代谢物差异表达模式负相关数量是正相关的2 倍，主要原因可能是胡麻根部盐胁迫下代谢物的敏感性较基因更强，这是因为盐胁迫首先造成植株代谢通路损害，导致大量的代谢产物不能富集，而细胞的失活过程相对较慢。

盐胁迫促进胡麻根部菜籽类固醇生物合成-菜籽类固醇通路和α-亚麻酸-茉莉酸通路上调。油菜素类固醇(brassinosteroids, BR)，在盐碱等逆境下能够通过增强作物根系吸水性从而增强植物的抗逆性[21]；茉莉酸及其衍生物的含量显著增加，参与根的内部组织层诱导耐盐机制[22]，主要诱导一系列与抗逆有关的基因表达，如蛋白酶抑制剂、硫蛋白和苯丙氨酸转氨酶(phenylalanine aminotransferase, PAL)等，提高酯氧合酶活性，从而增强植物的抗性[23]。本研究中盐胁迫条件下色氨酸-生长素通路、类胡萝卜素-脱落酸通路、半胱氨酸和蛋氨酸的代谢-乙烯等代谢通路下调。关于脱落酸对根系发育的调控机制了解相对较少[24]，有研究认为盐胁迫下脱落酸通过抑制侧根的起始和分生组织的活性来抑制其生长和发育，但是在逆境条件下参与侧根发育调控的信号网络还没有建立起来[25]。脱落酸调控植物耐盐性的关键在于SnRK2s 蛋白激酶对多种下游靶标基因的激活及转录调控，从而调控植物耐盐性[26]。生长素是侧根发育的主要调节物质，脱落酸可以影响生长素的积累和分布[27]。色氨酸是生长素生物合成的重要前体代谢物，而生长素的生物合成在盐胁迫下选择哪条途径，会因为器官不同而有所不同，但这些途径可以保证胡麻在盐胁迫环境下维持体内生长素的正常水平[28]。另外玉米素能促进植物细胞分裂，阻止蛋白质和叶绿素降解，让呼吸作用减慢，增强细胞活力，延缓衰老，从而达到抗逆的效果[29]。胡麻根部在盐胁迫下产生丰富的糖类代谢通路，原因是糖类、脂类、氨基酸三大营养素的代谢通路，又是其代谢联络与转化的枢纽，即糖类和蛋白质在体内是可以相互转化的，几乎所有组成蛋白质的天然氨基酸都可以转变成糖类，另外糖类代谢的中间产物可以转化成脂肪，而脂肪分解产生的甘油、脂肪酸又可以转化成糖类而被植物吸收[22]。细胞分裂素和脱落酸在逆境响应过程中可能存在拮抗关系，但是其调控的分子机制并不清楚，同样乙烯信号在盐胁迫反应中的作用还需要进一步研究[24]。

代谢组和转录组联合分析共筛选出8 个候选基因与5 个代谢物具有明显的相关性，每一个代谢物与一个或多个基因相关联。在色氨酸-生长素通路发现了GH3基因(IAA 酰胺合成酶)，主要作用是调控植物激素的平衡以及响应激素变化，在生长素信号途径和光信号途径以及植物的抗逆应急中发挥重要功能[30]。有关代谢物如何诱导相关基因富集并如何响应盐胁迫的机制等问题还需要进一步验证。

4 结论

本研究筛选出具有显著差异的代谢物741 个，其中下调409 个，上调332 个，筛选出差异基因共30 233 个，其中上调15 827 个，下调14 406 个，有注释结果的转录因子有13 015 个。将差异基因及差异代谢物同时映射到13 条相关KEGG 通路上，功能注释结果表明：胡麻功能转录因子主要富集在黄酮和黄酮醇的生物合成、苯丙素的生物合成、氨基糖和核苷酸糖的代谢、次生代谢产物的合成、花青素生物合成、类黄酮生物合成等通路，调控抗盐代谢途径的重要基因主要是以参与调控糖类代谢为主，增强细胞的渗透压，从而有效降低盐碱环境的伤害。本研究筛选出8 个候选基因与5 个代谢物具有明显的相关性，每一个代谢物与一个或多个基因相关联，其中发现3 对具有较高的关联度，分别是L-谷氨酸与Lus_GLEAN_10025781, 2-异丙基苹果酸与Lus_GLEAN_10022072, L-(-)-酪氨酸与Lus_GLEAN_10017173。另外色氨酸-生长素通路发现抗逆相关基因GH3，是一种典型的与植物生长素原初反应相关的基因，它在植物的抗逆防卫反应中起着重要作用，这也为后续胡麻差异表达基因的分析研究提供了思路。