通過整合轉錄組與代謝組解析不同類型椰子的脂肪酸調控網絡

鄧 淵，賴 軍，毛夢迪，張越冉，李 淳，楊 君，王守創，羅 杰

（海南大學 熱帶作物學院/海南省熱帶生物資源可持續利用重點實驗室，海口 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家族轉錄因子進行了進一步鑒定，預測并構建了椰子脂肪酸途徑的完整調控網絡，為后續脂肪酸途徑的進一步解析提供了有力資源。本研究不僅闡明了高矮種椰子間脂肪酸差異及相關分子機制，而且可為椰子油脂的遺傳改良和高油脂椰子品種的選育提供一定的依據。