通过整合转录组与代谢组解析不同类型椰子的脂肪酸调控网络

邓 渊，赖 军，毛梦迪，张越冉，李 淳，杨 君，王守创，罗 杰

（海南大学 热带作物学院/海南省热带生物资源可持续利用重点实验室，海口 570228）

椰子(Cocos nuciferaL.)是许多热带国家最重要的经济作物[1]。椰子属于棕榈科椰子属，基因组大小约为 2.4 Gb[2-3]。根据植物形态和繁殖类型，椰子可以分为 2 种类型：typica(tall)和nana(dwarf)[4]。椰子根据进化起源可以分为2个亚群——印度洋种群和太平洋种群[5-6]。高种椰子异花授粉，基因型表现为一定程度的杂合，而矮种椰子自花授粉，基因型纯合。矮种椰子生长缓慢、植株矮小、树冠较小、果小而多、椰肉厚度薄、开花时间较早，通常在种植 4～6 a 后开花，株高 10～15 m[7-9]。高种椰子耐寒，果实大，能够适应更广泛的环境，株高 20～30 m，初次开花时间晚 (种植 8～10 a后开花)。在中国，椰子主要种植在海南和云南等地区，海南地区的高种椰子“Hainan Tall”有着 36 000 hm2的种植面积[10]，高种椰子耐盐碱和高温，但对低温敏感。目前普遍认为矮种椰子经高种椰子驯化而来[4]，且矮种椰子更易受不利环境的影响。椰油占椰肉重量的65%[8,11]，与其他植物油相比,椰油具有独特的理化性质和药用价值。椰子油富含饱和脂肪酸(93%)，主要成分是己酸(C6)、辛酸(C8)、癸酸(C10)、月桂酸(C12)、肉豆蔻酸(C14)、棕榈酸 (C16)、硬脂酸 (C18)和花生酸 (C20)，其中85%是中链脂肪酸 (medium-chain fatty acids,MCFA)。据报道,月桂酸在棕榈科中的主要来源是椰子和油棕，两者都仅限于在热带和亚热带气候中生长[11-14]。月桂酸能够提高氧化稳定性，具有低熔点，能产生稳定的乳液，这对食品和化学工业应用具有重要意义[11]。目前月桂酸在医药和食品工业等各个领域都具有广泛的应用，包括用于制造肥皂、洗涤剂、纺织品、油漆、清漆、化妆品、药品等。肉豆蔻酸在食品工业中常用作调味剂和食品添加剂[15]。肉豆蔻酸有一定的药理活性，相较于其他饱和脂肪酸，肉豆蔻酸有更强的抗氧化活性[16]，肉豆蔻酸还可以降低血糖[17]。除了饱和脂肪酸以外，椰子油还含有少量不饱和脂肪酸，如油酸(C18:1)、棕榈油酸(C18:1)和亚麻酸(C18:3)。目前有报道使用椰子油可以降低阿尔茨海默病的患病风险[13-14]，抑制肿瘤生长、抑制氧化应激和神经炎症，以及消除抗肿瘤药甲氨蝶呤引起的脑神经毒性[18-24]。植物脂肪酸代谢途径在其他物种中已有广泛报道[25-34]，但目前在椰子中罕有报道。不同类型椰子品种在脂肪酸积累和合成调控途径是否存在差异，仍有待阐明。椰肉作为脂肪酸积累的主要部位，是研究脂肪酸合成调控的理想场所。本研究利用非靶向代谢组学，比较了海南高种、黄矮椰子和绿矮椰子的脂肪酸积累水平，研究表明，高种椰子的脂肪酸含量最高，黄矮椰子其次，绿矮椰子最低。结合转录组学数据，重建了椰子脂肪酸代谢途径，研究发现，相较于矮种椰子而言，脂肪酸途径的基因大多在高种椰子中高表达。WGCNA结果表明，“Melightyellow”模块中的基因与差异积累的脂肪酸有着相似的变化趋势，在此基础上结合模块中结构基因启动子区的cis-element，预测并构建了完整的脂肪酸调控网络。本研究通过组学策略探究了不同类型椰子品种脂肪酸差异的分子机理，预测并构建了完整的椰子脂肪酸调控网络，为后续椰子脂肪酸合成调控的进一步解析提供了可用资源，为高油脂椰子品种的选育提供了理论依据。

1 材料与方法

1.1 材料本试验材料由中国热带农业科学院椰子研究所提供，包括2份高种椰子(Hainan Tall,HT)，2份黄矮椰子 (Yellow Dwarf,YD)和 2份绿矮椰子 (Green Dwarf, GD)。每份材料随机选取3株长势一致的植株，每株植株采摘3个大小和发育时期一致的椰果，取椰肉混合后作为1个生物学重复，并将其一分为二，一份作为转录样，一份作为代谢样。将样品立即置于液氮冷冻，-80 ℃超低温冰箱保存备用。

1.2 代谢样品处理椰肉代谢样品经冷冻干燥后，采用BM500球磨仪研磨成粉末(频率30 Hz，持续时间 1.5 min)。称取 100 mg 研磨好的干粉样品，加入 1.2 mL 含有 0.1 mg·L-1利多卡因内标的φ= 70% 甲醇提取液，涡旋 30 s，低温静置 10 min，重复上述操作3次后，放入4℃冰箱继续萃取代谢产物 12 h，10 300 r·min-14 ℃ 离心 10 min，取上清液，用0.22 μm的滤膜过滤后，放入进样瓶中即可用于LC-MS进行代谢物检测分析。

1.3 转录数据处理对于原始的下机数据使用Fastp软件的默认参数去除低质量序列和接头序列，用HISAT2软件默认参数将其比对到矮种椰子的参考基因组，Samtools排序后使用featureCounts进行定量，用 TPM (transcripts per million)对原始定量结果进行标准化。R软件中的FactoMineR 和factoextra用于主成分分析。DEseq2用于样本间差异表达基因的鉴定，同时满足|log2FoldChange| ＞1 和 FDR (False Discovery Rate) ＜ 0.05 的基因被认为是2个样本比较间的差异表达基因。

2 结果与分析

2.1 基于非靶向代谢组学的椰子代谢谱比较为了探究不同类型椰肉的代谢组差异，选择了6份具有代表性的椰子种质资源，包括2份高种椰子(HT)、2份黄矮椰子(YD)、2份绿矮椰子(GD)，进行代谢物检测。通过 Q Exactive-Orbitrap UHPLC-ESI-MS/MS 在正离子模式下进行非靶向扫描，共检测到12 579种质谱信号。为了研究椰子果肉的代谢组构成，利用多种策略对所检测的质谱信号进行代谢物结构解析。对于能获取商业标准品的代谢物，通过比对保留时间(retention time, RT)，母离子(Q1)及质谱碎片信息等因素来鉴定代谢物。如 RT为 9.15 min的代谢信号(Cnm04191)，其母离子精确的m/z为 279.230 28，并且有149.0、219.1、57.1和105.0等一系列的特征碎片。与标准品比对，代谢信号Cnm04191由于相同的质谱碎片信息被鉴定为α-Linolenate(图1-A，1-B)。对于不能获得商业标准品的代谢物，将试验获取的质谱信息与已经公开发表的文献或者代谢数据库进行比对来初步解析代谢物结构。RT为 8.39 min的质谱信号 (Cnm09776)的MS1谱图中检测到准确的m/z为518.322 20，二级质谱碎片中观察到[Y0]+离子184.1，这是由于形成了磷酸胆碱这一特征碎片，因此，推测质谱信号Cnm09776 为 LysoPC (18:3/0:0) (图1-C)。质谱信号 Cnm06436被注释为 Glyceryl monolinoleate(图1-D)。通过 Compound Discoverer 3.1 软件自动批量匹配在线代谢数据库，包括mzCloud、ChemSpider、MassBank等进行代谢物结构鉴定。最终，超过 564 种代谢物被鉴定，主要为氨基酸、黄酮类化合物、有机酸、脂质 、苯丙烷 、糖类 、多酚、生物碱、核苷酸及其衍生物、植物激素。

图1 利用 UHPLC-HRMS 对次生代谢物进行检测和鉴定

基于检测到的97种脂质信号进行了主成分分析 (principal component analysis, PCA)，结果表明，每个品种的3个生物学重复都能较好地聚在一起，黄矮椰子和绿矮椰子聚集较近，但两者与高种椰子相聚较远(图2-A)。结果表明，同一样本不同重复间的代谢数据有较高的一致性，且不同类型的椰子品种在脂质积累水平有较大的差异，高矮种间的差异更大。进一步分析结果表明，高种椰子的脂质水平极显著(P＜ 2.22E-16)高于矮种椰子(图2-C)，且绿矮椰子的脂质积累水平最低(图2-D)。热图结果表明，检测到的97种脂质在3种不同类型的椰子果肉中具有不同的积累模式(图2-B)。为进一步探究三者间脂质积累存在差异的代谢产物，利用OPLS-DA对样本两两间差异代谢物进行鉴定 (YD vs HT, GD vs HT, GD vs YD)。在两两比较间脂质积累存在差异的物质共有50种，其中5种脂质在3组比较间都表现为积累差异(图2-E)。进一步对积累存在差异的代谢物进行注释，脂肪酰、甘油脂、鞘脂是主要的差异物质，其中脂肪酰占比最高达74%，鞘脂和甘油酯分别各占8% (图2-F)。糖脂主要与光合作用及磷响应相关[35]，鞘脂主要与植物生长、植物育性和胁迫应答等相关[36-38]。以上结果表明，这些物质在不同类型的椰子中迥异的积累模式，可能与它们响应生物和非生物胁迫以及适应不同地区的环境条件有关。

图2 不同类型椰子品种脂质代谢组的比较

2.2 不同类型椰子的转录组比较为进一步研究不同类型椰子品种在脂质通路的转录水平是否存在差异。对代表性品种的椰肉组织进行转录组测序。为了解不同类型椰子的转录动态，基于转录组数据进行PCA (图3-A)。结果表明，绝大多数的生物学重复彼此聚集，表明了转录数据的高可靠性。此外，这些转录本根据样本来源聚为3类，其中PC1捕获了43.5%的数据方差，代表不同类型的椰子样本之间的差异，这些结果均表明三者在基因的转录水平上存在较大的差异。随后对不同类型椰子品种的转录数据进行两两比较来鉴定不同样本间的差异表达基因(differently expressed genes, DEGs)。笔者在 HT vs YD，HT vs GD，GD vs YD 间分别鉴定到了 2 139、4 729 和 2 599 个DEGs，在两两比较间存在差异的基因共计5 841个，在3个比较中都存在差异表达的基因有384个(图3-B)。为了进一步探究差异基因的分子功能及其参与的生物学途径，对5 841个差异基因进行Kyoto Encyclopedia of Genes and Genomes (KEGG)和 Gene Ontology (GO)富集分析，KEGG 结果表明这些差异基因主要与植物激素信号转导、植物与病原菌互作、氨基酸代谢相关(图3-C)。GO富集分析结果也表明，不同椰子品种的表达差异主要体现在赤霉素信号通路、脂质储存和油菜素内酯响应等途径(图3-D)。

图3 不同品种椰子差异表达基因的鉴定与功能富集

2.3 WGCNA揭示脂肪酸合成调控网络采用加权共表达网络分析 (weighted gene co-expression network analysis, WGCNA)，将高连接度的基因定义为基因簇，同一簇内的基因具有较高的相关性。这些表达模式相近的基因在功能上具有相关性。根据表达模式的不同可以将5 841个差异基因划分为20个模块，其中模块“Melightyellow”中基因的表达模式与差异脂肪酸的积累模式高度相关(图4-A)。随后对“Melightyellow”中的基因进行蛋白功能注释和功能富集，其中参与“Fatty acid metabolism”，“Glycerolipid metabolism”，“Biosynthesis of unsaturated fatty acids”等脂肪酸相关通路的基因作为节点构建调控网络。基于保守功能结构域、蛋白功能注释和同源比对，笔者初步预测了13个可能参与脂肪酸代谢的基因，其中包括long chain acyl-CoA synthetase，glycerol-3-phosphate dehydrogenase，1-acyl-sn-glycerol-3-phosphate acyltransferase，lysophospholipid acyltransferase等。随后基于基因-基因、基因-代谢物、代谢物-代谢物间的皮尔森相关系数，结合脂肪酸代谢途径基因启动子区cis-elements的潜在结合能力，在“Melightyellow”模块中鉴定到了10个转录因子，包括ERF，MYB 和bHLH等。这些转录因子与“Melightyellow”模块中鉴定的13个候选基因，差异积累的脂肪酸都表现出较高的相关性(图4-B)。其中ADH (alcohol dehydrogenase,AZ13G0249870)基因参与 20-Hydroxyarachidonic acid 到 20-Oxoarachidonic acid (Cnm05349)的 合 成 ，而 LPC 18:1在LPCAT (lysophosphatidylcholine acyltransferase,AZ11G0210280)基因的参与下与Oleic acid结合生 成 Phosphatidylcholine (18:1/18:1, Cnm11965)(图4-C)。研究发现 20-Oxoarachidonic acid和Phosphatidylcholine (18:1/18:1)在高种椰子中的含量显著高于矮种椰子(图4-D)，转录组的数据也表明，ADH和LPCAT在高种椰子中有着更高的转录水平(图4-E)。为了进一步验证转录数据的可靠性，笔者通过qRT-PCR证实，ADH和LPCAT在高种椰子中的表达水平显著高于黄矮椰子和绿矮椰子(图4-F)。这些基因在高种椰子中的高水平表达可能是导致高种椰子脂质含量显著高于黄矮椰子和绿矮椰子的重要原因。

图4 加权共表达分析鉴定与脂肪酸合成相关的基因

2.4 椰子脂质合成途径的重建脂质的代谢和转录数据表明高种椰子在脂质积累水平和脂质合成基因的转录水平都显著高于矮种椰子。为了进一步探究脂肪酸合成通路在高矮种椰子间是否存在差异，构建了亚油酸相关的合成代谢通路。该通路由丙酮酸起始，在一系列酶的催化下生成亚油酸，亚油酸在酶的作用下最终生成茉莉酸甲酯(图5)。基于蛋白的同源性关系，对椰子脂肪酸合成通路的基因进行注释，包括ACCase、FAD3、LOX2S、AOS、AOC、OPR、OPCL1、JMT等基因。结合转录组数据的定量结果和代谢组的相对物质含量，笔者发现在脂肪酸上游途径，参与亚油酸、亚麻酸合成的的基因，如ACCase在高种椰子中的表达量高于矮种，而代谢物含量变化也与转录数据相一致。亚麻酸至9,10-Epoxy-10,12Z-octadecadienoate途径和亚油酸至8,9-DiHETrE途径及亚油酸至10-Oxo-11,15-phytodienoic acid和Heptadecatrienal途径，相关的代谢物在高种椰子中也始终表现为高积累。在代谢通路下游，茉莉酸甲酯合成途径中的关键基因如AOC、OPR、JMT等均在高种椰子中有着更高的表达水平，相关代谢物也在高种椰子中表现为高积累。以上结果均表明，高种椰子椰肉中的脂肪酸含量高于矮种椰子，这一现象产生的原因应该是脂肪酸合成通路的基因大多在高种椰子中特异性高表达。

图5 高矮种椰子脂肪酸合成途径基因表达水平及代谢产物积累水平的比较

3 讨 论

本研究不仅成功构建了不同类型椰子椰肉组织的代谢数据库，而且发现高种椰子的脂质积累水平远高于矮种椰子，不同类型的椰子品种在脂质的积累上具有特异性。此外通过整合转录组学和代谢组学的相关数据，发现“Melightyellow”模块与脂肪酸高度相关，并从中鉴定到了大量参与脂肪酸合成的结构基因。笔者通过分析脂肪酸合成途径基因启动子区的cis-element，发现这些基因的启动子区存在大量的激素信号响应元件和MYB家族转录因子的结合位点，这表明MYB家族在脂肪酸合成调控中发挥重要作用，而且激素作为一种信号分子可能通过响应外界环境条件的变化来参与到脂肪酸的合成途径中。蜡质主要由脂肪酸和相关衍生物构成，目前大量的研究结果表明，蜡质的合成受到多种环境因素的诱导，如干旱胁迫等[26]，ABA等激素可以正向调控蜡质合成基因 的表 达[26-27]。 熊 程 等[28]的研 究结 果表 明，MYB31通过调控CER6的表达，协同调控番茄表皮蜡质的合成，并且CER6的表达进一步受到ABA的诱导。鞠延仑等[29]在葡萄中发现外源施加茉莉酸甲酯(MeJA)可以调控脂肪酸代谢途径的脂氧合酶，进而对葡萄的风味造成影响。此外，研究者鉴定到大量的参与脂质途径的MYB转录因子，如MYB30和MYB41能够调控超长链脂肪酸的合成来影响棉花纤维的伸长[30-31]，MYB96和MYB94能够激活蜡质合成基因的表达[32-33]，烟草中的MYB12a能够促进脂肪酸的降解等[34]，这些均与本研究的结果相契合。本研究通过结合共表达分析和分析脂质途径结构基因启动子区的ciselement，对模块中可能参与脂质途径的MYB家族转录因子进行了进一步鉴定，预测并构建了椰子脂肪酸途径的完整调控网络，为后续脂肪酸途径的进一步解析提供了有力资源。本研究不仅阐明了高矮种椰子间脂肪酸差异及相关分子机制，而且可为椰子油脂的遗传改良和高油脂椰子品种的选育提供一定的依据。