基于整合方法分析茶樹響應(yīng)病原真菌脅迫的共有模式

連玲麗, 陳強, 周穎, 付婧, 李婉瑩, 魏日鳳, 劉偉

基于整合方法分析茶樹響應(yīng)病原真菌脅迫的共有模式

連玲麗1a, 陳強1b, 周穎1a, 付婧1a, 李婉瑩1a, 魏日鳳1b, 劉偉2*

(1. 福建農(nóng)林大學(xué), a. 生命科學(xué)學(xué)院; b. 園藝學(xué)院, 福州 350002; 2. 寧德師范學(xué)院, 福建 寧德 352100)

為探討茶樹()對病菌脅迫的共有響應(yīng)模式和抗病機制，運用生物信息學(xué)方法對多組RNA-seq數(shù)據(jù)進(jìn)行提取、整合及功能富集，結(jié)合多種工具和數(shù)據(jù)庫資源對主要調(diào)控分子及蛋白互作模塊加以分析。結(jié)果表明，病原真菌脅迫下，茶樹有較多細(xì)胞色素P450家族成員表達(dá)顯著上調(diào)；類固醇和激素的代謝過程、苯丙烷合成途徑被激活，有絲分裂細(xì)胞周期調(diào)控、DNA甲基化等生物過程及光合作用途徑受到抑制；主要調(diào)控分子如轉(zhuǎn)錄因子和、激酶和等以上調(diào)為主。差異表達(dá)的蛋白互作模塊分析表明，有絲分裂周期調(diào)控、基于微管運動、淀粉和蔗糖代謝、細(xì)胞壁多糖合成、光合作用、類黃酮代謝模塊明顯下調(diào)，木質(zhì)素合成和萜類生物合成模塊上調(diào)；且模塊之間可能存在互作。病菌脅迫激活的木質(zhì)素和萜類合成途徑的關(guān)鍵基因包括阿魏酸-5-羥基化酶基因、過氧化物酶基因和萜類合成酶基因等。細(xì)胞色素基因可能在病菌脅迫中起關(guān)鍵作用，增強木質(zhì)素和萜類物質(zhì)的合成、削弱光合作用可能是茶樹響應(yīng)真菌脅迫的核心模式。

茶樹; 病原真菌; 病菌脅迫；整合分析

茶樹()作為一種重要經(jīng)濟植物,在我國東部和南部各大茶區(qū)廣泛種植，近年來種植范圍更是向北延伸至河北、向西進(jìn)入西藏[1]。然而, 隨著全球氣候變化加劇，惡劣天氣使病蟲害發(fā)生呈現(xiàn)加重趨勢，給茶樹種植業(yè)帶來嚴(yán)重?fù)p失。其中同樣由炭疽病菌屬(spp.)引起的茶葉炭疽病和云紋葉枯病作為茶園中重要的葉部病害，更是對茶葉的產(chǎn)量和品質(zhì)造成直接影響[2]。如何提高茶樹對病菌脅迫的抗性是生產(chǎn)中亟待解決的關(guān)鍵問題。因此，本研究基于整合分析方法探究茶樹響應(yīng)病菌脅迫的共有途徑或基因，以期為明確茶樹抗病機制、提高茶樹抗病性等深入的研究提供參考。

伴隨著茶樹中國種基因組草圖和精細(xì)測序數(shù)據(jù)的公布[3–4]，人們得以結(jié)合組學(xué)技術(shù)更深入地探討茶樹在病菌脅迫下的響應(yīng)機制，Wang等[5]報道茶樹抗感品種對炭疽病菌()的響應(yīng)差異體現(xiàn)在苯丙烷和類黃酮等物質(zhì)合成的差別；Lu等[6]的研究表明茶樹與炭疽病菌的互作主要與胼胝質(zhì)的沉積、多種不同激素的信號轉(zhuǎn)導(dǎo)有關(guān)；茶樹響應(yīng)葉斑病菌(、)侵染時，淀粉與蔗糖代謝、苯丙烷生物合成、植物激素信號、類黃酮生物合成等代謝途徑呈現(xiàn)顯著增強[7–8]。這從不同角度展現(xiàn)了茶樹響應(yīng)病菌的可能機制，但未能提供茶樹抵抗病菌的共性信息。另一方面，茶樹在生產(chǎn)中常出現(xiàn)不同病菌復(fù)合侵染的現(xiàn)象，因而突破以往僅針對單一病菌脅迫的研究范式、探討茶樹對多種病菌脅迫響應(yīng)的共性機制, 具有重要的實踐意義。

整合分析方法作為一類統(tǒng)計方法，對圍繞同一問題的多個獨立研究進(jìn)行綜合分析，既能充分利用多項研究的信息，也能使我們對目標(biāo)問題有更深的理解，因而被廣泛應(yīng)用于醫(yī)學(xué)、農(nóng)學(xué)等領(lǐng)域[9]，近年來也越來越多地被用于組學(xué)數(shù)據(jù)的分析。如Ashrafi-Dehkordi等[10]對番茄()在不同脅迫下的轉(zhuǎn)錄響應(yīng)進(jìn)行整合分析，檢測得到響應(yīng)不同脅迫的通用基因，也發(fā)現(xiàn)了54個新的差異表達(dá)基因；Panahi等[11]通過跨物種的整合分析確定了重要的鹽脅迫響應(yīng)途徑等。目前，利用整合分析方法挖掘茶樹轉(zhuǎn)錄組數(shù)據(jù)資源的相關(guān)研究仍未見報道，我們以茶樹生產(chǎn)中常見的不同病菌復(fù)合脅迫問題為切入點，結(jié)合生物信息學(xué)方法開展相關(guān)研究。通過整合茶樹響應(yīng)不同病菌脅迫的組學(xué)數(shù)據(jù), 克服單一研究可能存在的局限性，確定整合方法在茶樹相關(guān)研究領(lǐng)域應(yīng)用的可行性；在獲取、整合差異表達(dá)基因數(shù)據(jù)集的基礎(chǔ)上，結(jié)合功能富集和網(wǎng)絡(luò)模塊等分析方法明確共有核心響應(yīng)分子及途徑，為后續(xù)的關(guān)鍵分子功能研究及育種應(yīng)用奠定理論和實踐基礎(chǔ)。

1 材料和方法

1.1 數(shù)據(jù)獲取

以“([Title]) AND infected[Text Word]”為關(guān)鍵詞在NCBI-SRA數(shù)據(jù)庫檢索，再通過過濾詞“Source：RNA”、“Library Layout：paired”和“Platform：Illumina”進(jìn)行過濾，獲取含有對照和病菌脅迫的茶樹RNA-seq數(shù)據(jù)；同時，以“(Tea Plant [Title]) OR ([Title]) AND (resistance [Abstract]) AND (transcriptome[Abstract])”為檢索詞搜索NCBI PubMed數(shù)據(jù)庫，手工選擇搜索結(jié)果中與病菌侵染茶樹有關(guān)的文獻(xiàn)，并從文獻(xiàn)中查找數(shù)據(jù)源。

通過數(shù)據(jù)庫檢索確定了PRJNA396805、PRJNA 528172、PRJNA637492和PRJNA564655等4項同時包含對照和病菌處理的茶樹轉(zhuǎn)錄組研究[5,7–8,12]，經(jīng)樣本間PCA分析，去除樣本分組效果不佳的部分研究,最終保留了3項研究(表1)共計24份表達(dá)數(shù)據(jù)用于整合分析。

表1 用于整合分析的轉(zhuǎn)錄組數(shù)據(jù)集

CK: 對照; TR: 接種病原。

CK: Control; TR: Inoculation of pathogen.

1.2 數(shù)據(jù)處理

采用FastQC v 0.11.8對每個樣本進(jìn)行質(zhì)量檢測,再由Trimmomatic 0.39去除接頭和低質(zhì)量區(qū)域，獲得清理后的讀段；繼而采用HISAT2將清理后的讀段匹配至茶樹“舒茶早”基因組，采用feature Counts對匹配結(jié)果進(jìn)行計數(shù)；最后使用DESeq2 (v 1.32.0)對各項研究進(jìn)行差異表達(dá)分析，視Log2Fold Change(Log2FC)絕對值>1、校正后值<0.05的基因為具有統(tǒng)計意義的差異表達(dá)基因(DEGs)。

1.3 整合分析鑒定差異基因

基于值合并法，通過metaRNASeq[13]的Fisher結(jié)合概率法對獨立研究的差異表達(dá)基因(individual- DEGs)進(jìn)行整合，其中采用Benjamini-Hochberg錯誤發(fā)現(xiàn)率對原始值進(jìn)行校正，取校正后值<0.05且表達(dá)趨向一致的基因作為整合分析后的差異表達(dá)基因(meta-DEGs)；繼而構(gòu)建individual-DEGs和meta-DEGs的韋恩圖，并提取單項研究與整合分析的共有差異表達(dá)基因(common-DEGs)用于后續(xù)功能分析。

1.4 功能預(yù)測與分析

為確定茶樹響應(yīng)病菌脅迫的關(guān)鍵生物過程或途徑，利用eggNOG在線工具(http://eggnog5.embl.de/)重新注釋茶樹基因組蛋白序列的GO和KEGG信息，再由clusterProfiler對meta-DEGs進(jìn)行功能富集分析，保留其中校正值<0.05的結(jié)果。借助iTAK v 1.6 (http://itak.feilab.net/cgi-bin/itak/online_itak. cgi)鑒定meta-DEGs中的轉(zhuǎn)錄因子和蛋白激酶，取值<1E-5的結(jié)果；同時基于miRNA的數(shù)據(jù)庫資源(plant microRNA database, http://bioinformatics.cau. edu.cn/PMRD/)，通過TAPIR靶標(biāo)預(yù)測工具(http://bio informatics.psb.ugent.be/webtools/tapir/)確定可能受miRNA調(diào)控的差異表達(dá)基因。借助STRING數(shù)據(jù)庫(http://string-db.org)基于高可信度互作記錄繪制meta-DEGs的蛋白互作網(wǎng)絡(luò)，以確定茶樹響應(yīng)病菌脅迫的關(guān)鍵模塊。

2 結(jié)果和分析

2.1 整合分析結(jié)果

將各項研究的基因表達(dá)信息分別匹配到茶樹“舒茶早”基因組的33 932個編碼基因上，并經(jīng)差異分析和整合處理，獲得3個獨立差異基因集(individual- DEGs)和1個整合差異基因集(meta-DEGs)。其中, 獨立差異基因集的基因數(shù)量差別較大，尤其是同樣與葉斑病脅迫相關(guān)的study2和study3，分別涉及2 515和10 699個差異基因。通過值合并法整合得到的meta-DEGs包括3 335個差異基因，涉及2 093個上調(diào)表達(dá)基因和1 242個下調(diào)表達(dá)基因。從individual- DEGs與meta-DEGs結(jié)果的交疊情況來看(圖1)，獨立差異基因集內(nèi)有23.76%～40.83%的基因仍保留于整合結(jié)果；而整合結(jié)果也有近40%的差異基因與至少2個獨立研究的表達(dá)趨勢一致，表明整合結(jié)果在較大程度上包含了多項研究的信息。

進(jìn)一步了解4個基因集的共有差異基因和整合分析的特有差異基因，可知共有差異基因集(common- DEGs)包含的305個基因中有約30%的基因編碼未知蛋白，40%的基因編碼酶類物質(zhì)，包括激酶、甲基化轉(zhuǎn)移酶和維生素C氧化酶等；整合分析獲得192個新的差異基因，除了多數(shù)編碼未知蛋白的基因之外，還有部分編碼受體樣蛋白、激酶及抗病蛋白的基因，總體上表達(dá)量的變化倍數(shù)為1～2。

圖1 獨立研究與整合分析的差異表達(dá)基因集的關(guān)系

以差異倍數(shù)的絕對值為主要排序依據(jù)，取meta-DEGs中表達(dá)量變化倍數(shù)最大且FDR值較小的前30個DEGs列于表2。編碼細(xì)胞色素P450的基因出現(xiàn)次數(shù)最多(8/30)且均上調(diào)表達(dá)，編碼漆酶、病程相關(guān)蛋白、激酶和幾丁質(zhì)酶等基因的差異表達(dá)水平也較高, 同樣呈現(xiàn)上調(diào)表達(dá)；這些基因中有不少是獨立研究及整合分析的共有基因，如細(xì)胞色素P450基因(TEA010837)、漆酶基因(TEA020596)等，表明它們更可能是茶樹響應(yīng)不同病菌脅迫的關(guān)鍵基因。

表2 整合分析鑒定的前30個差異表達(dá)基因

*: 共有DEGs。下同

*: Common-DEGs. The same below

2.2 整合結(jié)果的功能富集分析

對整合分析獲取的上調(diào)和下調(diào)差異表達(dá)基因分別進(jìn)行GO生物過程和KEGG代謝途徑的富集分析(表3)。上調(diào)表達(dá)基因的富集結(jié)果中，GO富集分析確定了病菌脅迫下茶樹的類固醇代謝、激素代謝、有毒物質(zhì)響應(yīng)、苯丙烷代謝等生物過程增強, KEGG分析同樣富集到油菜素類固醇合成、細(xì)胞分裂素(玉米素)合成、苯丙烷生物合成等相似的途徑; 而下調(diào)表達(dá)基因的富集結(jié)果則表明，病菌脅迫下茶樹的細(xì)胞壁組織生成、有絲分裂細(xì)胞周期調(diào)節(jié)、DNA復(fù)制、微管解聚等生物過程減弱，黃酮類和萜類等物質(zhì)的合成、淀粉與蔗糖代謝、光合作用等相關(guān)代謝途徑受到抑制。

表3 整合分析后差異表達(dá)基因的功能富集

以上調(diào)為主的苯丙烷合成途徑和顯著下調(diào)的光合作用途徑也分別被common-DEGs的上調(diào)和下調(diào)子集所富集，其中富集到苯丙烷合成途徑的關(guān)鍵基因有編碼阿魏酸-5-羥基化酶的基因(TEA032005和TEA000057)、過氧化物酶基因(TEA001789、TEA028696)、糖苷水解酶基因(TEA004253)和-葡萄糖苷酶基因(TEA001710)；富集到光合作用途徑的基因有編碼光系統(tǒng)I的和基因(TEA017139、TEA001351、TEA002548和TEA033829),表明這2個代謝途徑在茶樹響應(yīng)病菌脅迫中起重要作用。

2.3 響應(yīng)病菌脅迫的轉(zhuǎn)錄因子、激酶和miRNA的鑒定

轉(zhuǎn)錄因子(TFs)和激酶(PKs)在植物響應(yīng)生物脅迫中發(fā)揮重要作用。通過將meta-DEGs比對到轉(zhuǎn)錄因子數(shù)據(jù)庫中，獲得207個轉(zhuǎn)錄因子，歸屬于38個家族，其中、/、、、和等家族均有10多個成員發(fā)生顯著表達(dá)變化。差異表達(dá)基因中，上調(diào)表達(dá)的數(shù)量明顯少于下調(diào)表達(dá)的(63 up/174 down)，表明病菌脅迫下轉(zhuǎn)錄因子的表達(dá)普遍受到抑制。但也有部分轉(zhuǎn)錄因子的表達(dá)以激活為主，如家族(13up/1down)和家族(9 up/7 down)；還有部分轉(zhuǎn)錄因子的表達(dá)雖以抑制為主，但仍有一定數(shù)量的成員被激活，如家族(5 up/20 down)和(4 up/21 down), 表明這些轉(zhuǎn)錄因子相較于其他轉(zhuǎn)錄因子在茶樹響應(yīng)病菌脅迫中發(fā)揮更為重要的作用。

圖2 差異表達(dá)基因中轉(zhuǎn)錄因子家族的分布

利用iTAK工具在meta-DEGs中鑒定得到308個蛋白激酶，涉及AGC、CAMK、CMGC、Group-Pl、RLK、STE和TKL等多種激酶家族，其中受體蛋白激酶RLK的成員數(shù)量最多(269個)、其次是鈣調(diào)激酶CAMK(6個)。這2個家族的成員均以上調(diào)表達(dá)為主，分別有160和4個成員表達(dá)量顯著上升; 其他家族成員則多數(shù)是下調(diào)表達(dá)。在差異表達(dá)的家族中，亞家族、、、和的絕大多數(shù)成員呈現(xiàn)上調(diào)表達(dá)；相反地，亞家族、等的成員則多為下調(diào)表達(dá)，表明的亞家族可能在茶樹響應(yīng)病菌脅迫中扮演不同角色。

miRNA參與基因表達(dá)的調(diào)控。通過miRNA-靶標(biāo)基因關(guān)系的預(yù)測，可知病菌脅迫下78個下調(diào)表達(dá)的茶樹基因可能受到53個miRNA (歸屬于23個miRNA家族)的調(diào)控。從表4可見，調(diào)控的基因數(shù)量最多(34個)，這些基因參與有絲分裂細(xì)胞周期過程、植物類細(xì)胞壁組織發(fā)生等生物過程，其次和分別抑制漆酶和生長調(diào)節(jié)因子的表達(dá)，表明這些miRNA可能是重要的脅迫相關(guān)調(diào)控因子。此外，和等保守家族也參與茶樹的脅迫響應(yīng)，調(diào)控的基因負(fù)責(zé)編碼轉(zhuǎn)錄因子、K+外排逆向蛋白和同源盒亮氨酸拉鏈蛋白。

圖3 差異表達(dá)基因中蛋白激酶RLK亞家族的分布

2.4 蛋白互作網(wǎng)絡(luò)分析

將meta-DEGs的基因映射到擬南芥蛋白互作關(guān)系中，預(yù)測病菌脅迫下茶樹差異表達(dá)基因之間的可能互作模式。由構(gòu)建的蛋白互作網(wǎng)絡(luò)中提取核心區(qū)域作圖，并分析其中的主要功能模塊(圖4)，獲得的8個功能模塊分別是有絲分裂周期調(diào)控(M1)、基于微管運動(M2)、淀粉與蔗糖代謝(M3)、細(xì)胞壁多糖的合成(M4)、光合作用(M5)、類黃酮和類固醇合成代謝(M6)、木質(zhì)素合成(M7)、萜類合成(M8)。其中模塊M1～M5的節(jié)點對應(yīng)基因幾乎都呈現(xiàn)下調(diào)表達(dá), 模塊M6中與黃酮醇的合成、花青素合成有關(guān)的節(jié)點對應(yīng)基因(、、、)下調(diào)表達(dá)，與油菜素甾醇合成有關(guān)的節(jié)點(細(xì)胞色素P450)對應(yīng)基因(TEA015379、TEA015397)上調(diào)表達(dá)；M7和M8的節(jié)點均以上調(diào)表達(dá)為主，分別與木質(zhì)素的合成與代謝、萜類骨架合成有關(guān)。

從圖4可見，模塊之間相對獨立，但彼此又存在不同程度的交互作用，其中模塊M1和M2、M6和M7間的關(guān)系緊密，有較多跨模塊的互作，如模塊M6的黃酮醇合酶FLS (TEA010328)、二氫黃酮醇還原酶DFR (TEA024762)均與模塊M7的莽草酸羥基肉桂酸酰轉(zhuǎn)移酶HCT (TEA022314、TEA011691)互作；其他模塊間則通過單對蛋白節(jié)點的互作建立聯(lián)系，如模塊M5的早期光誘導(dǎo)蛋白ELIP1 (TEA009875)與模塊M6的DFR互作、模塊M5的PSII色素結(jié)合蛋白NPQ4 (TEA019348)與模塊M1的蛋白激酶STN8 (TEA020777)互作、模塊M5的核酮糖-1,5-二磷酸羧化酶/加氧酶Rubisco (TEA004730)與模塊3的26S蛋白酶體非ATP酶調(diào)節(jié)亞基RPN (TEA000905)互作等。這些參與模塊間互作的橋梁分子可能在模塊的協(xié)同響應(yīng)病菌脅迫中發(fā)揮作用。

表4 預(yù)測的miRNA及其靶向的下調(diào)基因

*: 受2個miRNA家族調(diào)控; **: 功能未知, 共計21個。

*: Regulated by two miRNA families; **: Unknown function, 21 in total.

圖4 茶樹響應(yīng)病菌脅迫的核心網(wǎng)絡(luò)的模塊分布。三角形: 表達(dá)上調(diào); 圓形: 表達(dá)下調(diào); 方形: 表達(dá)上下調(diào)兼有。M1:有絲分裂周期調(diào)控; M2: 基于微管運動; M3: 淀粉與蔗糖代謝; M4: 細(xì)胞壁多糖合成; M5: 光合作用; M6: 類黃酮和類固醇合成; M7: 木質(zhì)素合成; M8: 萜類合成。

進(jìn)一步提取以上調(diào)表達(dá)為主的模塊M7和M8的互作關(guān)系，并將其映射到代謝途徑中。從圖5可見，模塊M7中木質(zhì)素合成相關(guān)酶類的基因(、、、、和)及萜類合成相關(guān)的催化酶基因(、、、和) 幾乎都是上調(diào)表達(dá)，其中(TEA000057和TEA 032005)、(TEA028696)和(TEA021767)不僅在不同病菌脅迫下均呈現(xiàn)上調(diào)表達(dá)，且上調(diào)程度最為顯著，表明它們是這些代謝途徑中更為重要的分子。從代謝途徑的蛋白互作情況來看，由l-苯丙氨酸到木質(zhì)素的代謝流內(nèi)存在的互作關(guān)系主要發(fā)生在相鄰酶分子之間，也有個別酶分子與非相鄰酶分子互作，如PAL與COMT、CcoAOMT、CAD等的互作，表明這類酶可能在整個代謝調(diào)控中占據(jù)重要地位。

圖5 木質(zhì)素合成模塊(A)和萜類的甲羥戊酸合成模塊(B)及蛋白互作。紅色: 差異表達(dá)的酶基因; 黃色: 上調(diào)表達(dá)基因; 綠色: 下調(diào)表達(dá)基因; *: 共有基因; 虛線: 互作關(guān)系。

3 結(jié)論和討論

測序技術(shù)的發(fā)展和轉(zhuǎn)錄組研究數(shù)量的增多，使得我們有可能基于多份組學(xué)數(shù)據(jù)獲取共有的表達(dá)信息，以減小單次研究可能存在的分析偏頗。本研究首次嘗試對多份茶樹響應(yīng)病菌脅迫的組學(xué)數(shù)據(jù)進(jìn)行整合分析，以挖掘共有信息。在整合得到的3 335個差異表達(dá)基因(DEGs)中，顯著下調(diào)的DEGs居多(62.76%)，但在表達(dá)倍數(shù)變化最大的前30個DEGs中，顯著上調(diào)表達(dá)的更多(76.67%)，表明茶樹在遭受病菌侵染時可能從整體上降低代謝水平, 并將有限的能量用于激活特定的生物分子。這些被激活的分子中，存在多個細(xì)胞色素P450家族成員, 還包括碳酸酐酶(TEA010186)、病程相關(guān)蛋白(TEA017523)和漆酶(TEA020596)等抗病相關(guān)基因, 它們的功能在以往研究中雖有報道[14–16]，但在茶樹響應(yīng)病菌脅迫中的作用仍有待明確。

為了解茶樹在病菌脅迫下基因調(diào)控的有關(guān)信息，對差異表達(dá)的轉(zhuǎn)錄因子、激酶和miRNA進(jìn)行分析。結(jié)果表明，病菌脅迫下來自38個家族的207個轉(zhuǎn)錄因子的表達(dá)發(fā)生明顯變化，其中和轉(zhuǎn)錄因子的成員數(shù)目最多，以下調(diào)表達(dá)為主;和轉(zhuǎn)錄因子成員數(shù)目次之，以上調(diào)表達(dá)為主，表明這些家族可能參與茶樹響應(yīng)真菌脅迫的負(fù)向或正向調(diào)控。在這4類轉(zhuǎn)錄因子中，茶樹家族的功能研究較多，不僅通過qRT-PCR技術(shù)驗證了多個茶樹基因在植株接種炭疽病菌后呈現(xiàn)表達(dá)量的顯著上調(diào)[17]、還利用基因沉默試驗證明缺失基因的植株對病原真菌更為敏感[18];相比而言，其他3類茶樹轉(zhuǎn)錄因子的功能研究較少。盡管如此，其他物種的相關(guān)研究證實了它們同樣在植物響應(yīng)病菌脅迫中發(fā)揮重要功能，轉(zhuǎn)錄因子對水稻()的稻瘟病抗性、小麥()的白粉病抗性、番茄的灰霉病抗性及馬鈴薯()的晚疫病抗性等均具有正向調(diào)控作用[19]；煙草()的基因沉默使植株對晚疫病菌更加敏感[20]，辣椒()的基因的過表達(dá)增強植株對病原真菌的抗性[21]；棉花()的轉(zhuǎn)錄因子正向調(diào)控植株對黃萎病菌的抗性[22]等。因此，這些轉(zhuǎn)錄因子很可能也在茶樹抗病響應(yīng)的基因調(diào)控中扮演重要角色。對茶樹在病菌脅迫下差異表達(dá)的蛋白激酶(PKs)進(jìn)行分析，結(jié)果檢測到308個PKs，其中絕大部分來自家族，其次是家族。在其他植物全基因組水平的PKs分析中，這2個家族的數(shù)量也較多, 兩者均參與植物的逆境脅迫響應(yīng)，葡萄()的CAMK家族成員在植株響應(yīng)非生物脅迫中發(fā)揮作用[23]；豇豆()接種病毒后，亞家族和家族上調(diào)表達(dá)的成員數(shù)量明顯多于其他家族[24]。在茶樹中，家族成員參與植株對寒害脅迫的響應(yīng)[25]，至于亞家族、家族在茶樹響應(yīng)病菌脅迫中的作用則未見報道。因此，開展相關(guān)的后續(xù)研究有助于解析兩者在茶樹對抗真菌病害中的具體功能。對茶樹在病菌脅迫下差異表達(dá)的miRNA進(jìn)行預(yù)測，得到來自23個家族的53個miRNA，其中靶標(biāo)基因數(shù)量或成員數(shù)量較多的miRNA來自、、、等家族，它們的靶向基因主要是轉(zhuǎn)錄因子。Jeyaraj等[26]預(yù)測、、和在茶樹受炭疽菌侵染后發(fā)生顯著的表達(dá)變化，且可能靶向調(diào)控轉(zhuǎn)錄因子(如、)和抗氧化酶基因(如、)的表達(dá)；更多的研究表明，和家族成員參與植物對不同病原真菌的防御響應(yīng)，如黃瓜()對霜霉病菌侵染[27]、水稻對稻瘟菌侵染的響應(yīng)[28]等。因此，這2個miRNA家族可能在茶樹抗病響應(yīng)中發(fā)揮相似的作用。

為進(jìn)一步明確茶樹在病菌脅迫下發(fā)生變化的生物過程及代謝模塊，對整合差異表達(dá)基因(meta- DEGs)進(jìn)行功能富集分析。結(jié)果表明，病菌脅迫下類固醇代謝、激素代謝和苯丙烷代謝等顯著增強,這與單項研究的結(jié)果基本一致。結(jié)合KEGG富集分析結(jié)果，相關(guān)代謝被進(jìn)一步細(xì)化至苯丙烷、油菜素內(nèi)酯和玉米素的合成途徑，其中苯丙烷合成途徑也在油茶響應(yīng)炭疽病菌脅迫時被激活[29]，表明該合成途徑可能是茶樹響應(yīng)病菌的核心途徑。結(jié)合蛋白互作網(wǎng)絡(luò)模塊來看，與苯丙烷代謝密切相關(guān)的2個分支途徑中，木質(zhì)素合成模塊內(nèi)合成相關(guān)酶對應(yīng)的基因明顯上調(diào)表達(dá)，類黃酮合成模塊內(nèi)與花青素、黃酮醇合成相關(guān)的酶類基本呈現(xiàn)下調(diào)表達(dá)，其中前者的代謝與植物抗性相關(guān)，后者的代謝則與植物品質(zhì)有關(guān)[30]，表明茶樹在響應(yīng)病菌脅迫時可能以損失品質(zhì)為代價來加強自身的抗病性。類似地，甲羥戊酸(油菜素內(nèi)酯及玉米素的前體物質(zhì))的合成途徑也在病菌脅迫時被激活。在月季()抗灰霉病菌的轉(zhuǎn)錄組學(xué)分析中，與油菜素內(nèi)酯和乙烯相關(guān)的基因表達(dá)水平都明顯高于水楊酸和茉莉酸的，且油菜素內(nèi)酯在脅迫響應(yīng)中發(fā)揮主導(dǎo)作用[31]。由此可見，木質(zhì)素合成途徑、油菜素內(nèi)酯合成途徑可能是茶樹防御病菌脅迫的重要代謝途徑。

除上調(diào)表達(dá)的代謝途徑外，整合差異表達(dá)基因的功能富集分析還表明，茶樹在病菌脅迫下受到抑制的代謝途徑，包括基于微管的過程、有絲分裂細(xì)胞周期及其調(diào)節(jié)、淀粉和蔗糖代謝、光合作用等途徑的顯著下調(diào)。其中，光合作用不僅在其他植物響應(yīng)病原真菌脅迫中發(fā)生明顯的下調(diào)變化[32]，也在植物響應(yīng)寒冷、高溫、病原細(xì)菌等多種脅迫中受到抑制[33]，表明光合作用削弱可能是植物逆境響應(yīng)的一般模式。編碼與光合作用相關(guān)的光系統(tǒng)I組成蛋白psaA和psaB、光系統(tǒng)II的PSBS蛋白和PSBO蛋白、ATP合酶亞基等的基因均呈現(xiàn)下調(diào)表達(dá)，盡管這些光系統(tǒng)組分在植物抗真菌病害中的作用仍鮮見報道，但已有研究[34]表明，病原細(xì)菌效應(yīng)蛋白及病毒蛋白均能靶向光系統(tǒng)組分蛋白以破壞光系統(tǒng), 也有研究表明[35–36]增強光系統(tǒng)組分的表達(dá)或活性能提高植株對病原細(xì)菌或花葉病毒病的抗性，這些研究表明光合作用系統(tǒng)在植物-病原互作中占據(jù)重要位置。除此之外，結(jié)合蛋白互作網(wǎng)絡(luò)模塊分析的結(jié)果，光合作用模塊與有絲分裂周期調(diào)控模塊、淀粉與蔗糖代謝模塊及次生代謝物模塊之間可能存在模塊間的互作，這與“光合作用可能是調(diào)節(jié)及平衡基礎(chǔ)代謝與防御相關(guān)代謝的重要樞紐[34]”的觀點一致。因此，明確茶樹光合系統(tǒng)在調(diào)控防御響應(yīng)中發(fā)揮的作用將有助于加速抗病機制、品種選育等研究的進(jìn)程。

結(jié)合蛋白互作預(yù)測方法了解茶樹在病菌脅迫下可能存在的蛋白互作信息，著重對以上調(diào)表達(dá)為主的木質(zhì)素合成模塊和萜類合成模塊的互作進(jìn)行分析。結(jié)果表明，病菌脅迫下茶樹的這2個合成途徑中相關(guān)合成酶之間存在較多互作。在木質(zhì)素合成途徑中，預(yù)測的PAL與4CL、COMT、CcoAOMT的互作，HCT與4CL、CSE的互作，以及CAD與CcoAOMT的互作均在毛果楊()木質(zhì)素生物合成酶的互作試驗中得到驗證，且多數(shù)為瞬時互作[37]，其中HCT與4CL能形成蛋白復(fù)合體促進(jìn)木質(zhì)素的合成[38]；合成途徑中POD與CAD、UGT72E、F5H的互作則僅出現(xiàn)于其他相似的預(yù)測研究[39]中，尚未有試驗驗證的報道。在萜類合成途徑中，合成相關(guān)酶AACT、HMGS、HMGR、MPD和IDI之間可能存在互作，其中AACT與HMGS的互作在古細(xì)菌和真細(xì)菌中得到證實[40]，而HMGR和HMGS、MPD的互作則出現(xiàn)于葡萄的HMGR蛋白互作預(yù)測結(jié)果[41]中。這些研究暗示了合成途徑中酶的互作既有保證正常代謝的穩(wěn)定互作，也有應(yīng)對環(huán)境改變而產(chǎn)生的瞬時互作，同時酶的互作可能不同程度地調(diào)控代謝終產(chǎn)物的合成量，進(jìn)而影響植物表型。因而，探明茶樹中與病菌防御相關(guān)代謝途徑的完整互作網(wǎng)絡(luò)及逆境脅迫下的動態(tài)互作網(wǎng)絡(luò), 將有助于從分子層面揭示茶樹的抗病分子機制。

綜上，茶樹在響應(yīng)不同病原真菌脅迫時均增加木質(zhì)素合成、降低光合作用，且更傾向于通過甲羥戊酸途徑合成萜類物質(zhì)或激素。、、和等轉(zhuǎn)錄因子為參與脅迫響應(yīng)的重要分子,這為后續(xù)的關(guān)鍵途徑與分子的功能研究、茶樹的抗病品種選育提供參考信息。

[1] TANG J X, WANG P J, E Y H, et al. Climatic suitability zoning of tea planting in China’s mainland [J]. J Appl Meteor Sci, 2021, 32(4): 397– 407. [唐俊賢, 王培娟, 俄有浩, 等. 中國大陸茶樹種植氣候適宜性區(qū)劃 [J]. 應(yīng)用氣象學(xué)報, 2021, 32(4): 397–407.]

[2] CAI X M, LUO Z X, BIAN L, et al. Tea pest prevention and control progress during the 13th five-year plan period and development direction in the 14th five-year plan period [J]. China Tea, 2021, 43(9): 66–73. [蔡曉明, 羅宗秀, 邊磊, 等. 茶園綠色防控“十三五”進(jìn)展及“十四五”發(fā)展方向 [J]. 中國茶葉, 2021, 43(9): 66–73. doi: 10.3969/ j.issn.1000-3150.2021.09.010.]

[3] WEI C L, YANG H, WANG S B, et al. Draft genome sequence ofvar.provides insights into the evolution of the tea genome and tea quality [J]. Proc Natl Acad Sci USA, 2018, 115(18): E4151–E4158. doi: 10.1073/pnas.1719622115.

[4] CHEN J D, ZHENG C, MA J Q, et al. The chromosome-scale genome reveals the evolution and diversification after the recent tetraploidi- zation event in tea plant [J]. Hort Res, 2020, 7: 63. doi: 10.1038/s414 38-020-0288-2.

[5] WANG Y C, HAO X Y, LU Q H, et al. Transcriptional analysis and histochemistry reveal that hypersensitive cell death and H2O2have crucial roles in the resistance of tea plant [(L.) O. Kuntze] to anthracnose [J]. Hort Res, 2018, 5: 18. doi: 10.1038/s414 38-018-0025-2.

[6] LU Q H, WANG Y C, XIONG F, et al. Integrated transcriptomic and metabolomic analyses reveal the effects of callose deposition and multihormone signal transduction pathways on the tea plant-interaction [J]. Sci Rep, 2020, 10(1): 12858. doi: 10. 1038/s41598-020-69729-x..

[7] YANG H, LUO P G. Changes in photosynthesis could provide impor- tant insight into the interaction between wheat and fungal pathogens [J]. Int J Mol Sci, 2021, 22(16): 8865. doi: 10.3390/ijms 22168865.

[8] JIANG S L, YIN Q X, LI D X, et al. Integrated mRNA and small RNA sequencing for analyzing tea leaf spot pathogen, underconditions and the course of infection [J]. Phytopathology, 2021, 111(5): 882–885. doi: 10.1094/PHYTO-07-20- 0297-A.

[9] XU R F, YANG K, DING J, et al. Effect of green tea supplementation on blood pressure: A systematic review and meta-analysis of rando- mized controlled trials [J]. Medicine, 2020, 99(6): e19047. doi: 10. 1097/MD.0000000000019047.

[10] ASHRAFI-DEHKORDI E, ALEMZADEH A, TANAKA N, et al. Meta-analysis of transcriptomic responses to biotic and abiotic stress in tomato [J]. PeerJ, 2018, 6: e4631. doi: 10.7717/peerj.4631.

[11] PANAHI B, FRAHADIAN M, DUMS J T, et al. Integration of cross species RNA-seq meta-analysis and machine-learning models identifies the most important salt stress-responsive pathways in[J]. Front Genet, 2019, 10: 752. doi: 10.3389/fgene.2019.00752.

[12] WANG S S, LIU L, MI X Z, et al. Multi-omics analysis to visualize the dynamic roles of defense genes in the response of tea plants to gray blight [J]. Plant J, 2021, 106(3): 862–875. doi: 10.1111/tpj.15203.

[13] RAU A, MAROT G, JAFFRéZIC F. Differential meta-analysis of RNA-seq data from multiple studies [J]. BMC Bioinform, 2014, 15: 91. doi: 10.1186/1471-2105-15-91.

[14] PANDIAN B A, SATHISHRAJ R, DJANAGUIRAMAN M, et al. Role of cytochrome P450 enzymes in plant stress response [J]. Antioxidants, 2020, 9(5): 454. doi: 10.3390/antiox9050454.

[15] FLORYSZAK-WIECZOREK J, ARASIMOWICZ-JELONEK M. The multifunctional face of plant carbonic anhydrase [J]. Plant Physiol Biochem, 2017, 112: 362–368. doi: 10.1016/j.plaphy.2017.01.007.

[16] BARROS J, SERK H, GRANLUND I, et al. The cell biology of lignification in higher plants [J]. Ann Bot, 2015, 115(7): 1053–1074. doi: 10.1093/aob/mcv046.

[17] WEI R F, LAI J D, PENG C B, et al. cDNA-AFLP reveals differential gene expression profiles of tea plant (cv. Maoxie) induced bysp.1 infection [J]. J Tea Sci, 2020, 40(1): 26–38. [魏日鳳, 賴建東, 彭成彬, 等. 利用cDNA-AFLP技術(shù)分析炭疽菌危害誘導(dǎo)茶樹的差異表達(dá)基因 [J]. 茶葉科學(xué), 2020, 40(1): 26–38. doi: 10.3969/j.issn.1000-369X.2020.01.003.]

[18] LIU S Y, ZHANG Q Q, GUAN C F, et al. Transcription factor WRKY14 mediates resistance of tea plants [(L.) O. Kuntze] to blister blight [J]. Physiol Mol Plant Pathol, 2021, 155: 101667. doi: 10.1016/j.pmpp.2021.101667.

[19] BIAN Z Y, GAO H H, WANG C Y. NAC transcription factors as positive or negative regulators during ongoing battle between patho- gens and our food crops [J]. Int J Mol Sci, 2021, 22(1): 81. doi: 10. 3390/ijms22010081.

[20] IMANO S, FUSHIMI M, CAMAGNA M, et al. AP2/ERF transcription factor NbERF-IX-33 is involved in the regulation of phytoalexin production for the resistance ofto[J]. Front Plant Sci, 2022, 12: 821574. doi: 10.3389/ fpls.2021.821574.

[21] HUH S U. Functional analysis of hot pepper ethylene responsive factor 1A in plant defense [J]. Plant Signal Behav, 2022: 2027137. doi: 10. 1080/15592324.2022.2027137.

[22] ZHU Y T, HU X Q, WANG P, et al., an R2R3-type MYB transcription factor, positively regulates cotton resistance tothe lignin biosynthesis and jasmonic acid signaling pathway [J]. Int J Biol Macromol, 2022, 201: 580–591. doi: 10.1016/j. ijbiomac.2022.01.120.

[23] ZHU K K, WANG X L, LIU J Y, et al. The grapevine kinome: Annotation, classification and expression patterns in developmental processes and stress responses [J]. Hort Res, 2018, 5: 19. doi: 10.1038/ s41438-018-0027-0.

[24] FERREIRA-NETO J R C, DA COSTA B A N, DA SILVA M D, et al. The cowpea kinome: Genomic and transcriptomic analysis under biotic and abiotic stresses [J]. Front Plant Sci, 2021, 12: 667013. doi: 10. 3389/fpls.2021.667013.

[25] DING C Q, LEI L, YAO L N, et al. The involvements of calcium- dependent protein kinases and catechins in tea plant [(L.) O. Kuntze] cold responses [J]. Plant Physiol Biochem, 2019, 143: 190–202. doi: 10.1016/j.plaphy.2019.09.005.

[26] JEYARAJ A, WANG X W, WANG S S, et al. Identification of regulatory networks of microRNAs and their targets in response toin tea plant (L.) [J]. Front Plant Sci, 2019, 10: 1096. doi: 10.3389/fpls.2019.01096.

[27] JIN W B, WU F L. Identification and characterization of cucumber microRNAs in response toinfection [J]. Gene, 2015, 569(2): 225–232. doi: 10.1016/j.gene.2015.05.064.

[28] LI Y, LU Y G, SHI Y, et al. Multiple rice microRNAs are involved in immunity against the blast fungus[J]. Plant Physiol,2014, 164(2): 1077–1092. doi: 10.1104/pp.113.230052.

[29] YANG C C, WU P F, YAO X H, et al. Integrated transcriptome and metabolome analysis reveals key metabolites involved indefense against anthracnose [J]. Int J Mol Sci, 2022, 23(1): 536. doi: 10.3390/ijms23010536.

[30] JIA X L, WANG G L, XIONG F, et al.assembly, trans- criptome characterization, lignin accumulation and anatomic character- ristics: Novel insights into lignin biosynthesis during celery leaf development [J]. Sci Rep, 2015, 5: 8259. doi: 10.1038/srep08259.

[31] LIU X T, CAO X Q, SHI S C, et al. Comparative RNA-Seq analysis reveals a critical role for brassinosteroids in rose () petal defense againstinfection [J]. BMC Genet, 2018, 19(1): 62. doi: 10.1186/s12863-018-0668-x.

[32] YANG R, JIANG S L, LI D X, et al. Integrated mRNA and small RNA sequencing for analyzing leaf spot pathogenand its host, tea (), during infection [J]. Mol Plant Microbe Interact, 2021, 34(1): 127–130. doi: 10.1094/MPMI-07-20- 0207-A.

[33] ZHANG Y, ZHANG A H, LI X M, et al. The role of chloroplast gene expression in plant responses to environmental stress [J]. Int J Mol Sci, 2020, 21(17): 6082. doi: 10.3390/ijms21176082.

[34] KANGASJ?RVI S, TIKKANEN M, DURIAN G, et al. Photosynthetic light reactions: An adjustable hub in basic production and plant immunity signaling [J]. Plant Physiol Biochem, 2014, 81: 128–134. doi: 10.1016/j.plaphy.2013.12.004.

[35] SHASMITA, SAMAL P, MOHAPATRA P K, et al. Improved photo- system II and defense enzymes activity in rice () by biopriming againstpv.[J]. Funct Plant Biol, 2021, 48(3): 298–311. doi: 10.1071/FP20221.

[36] BWALYA J, ALAZEM M, KIM K H. Photosynthesis-related genes induce resistance against soybean mosaic virus: Evidence for involve- ment of the RNA silencing pathway [J]. Mol Plant Pathol, 2022, 23(4): 543–560. doi: 10.1111/mpp.13177.

[37]LIU B G. Transcription regulation of lignin biosynthesis and interaction of monolignol biosynthesis enzymes in[D]. Harbin: Northeast Forestry University, 2020. [劉寶光. 毛果楊木質(zhì)素生物合成的轉(zhuǎn)錄調(diào)控及單體生物合成酶互作研究 [D]. 哈爾濱: 東北林業(yè)大學(xué), 2020.]

[38] LIN CY, SUN Y, SONG J N, et al. Enzyme complexes of Ptr4CL and PtrHCT modulate co-enzyme a ligation of hydroxycinnamic acids for monolignol biosynthesis in[J]. Front Plant Sci, 2021, 12: 727932. doi: 10.3389/fpls.2021.727932.

[39] LI X, ZHANG Y Y, ZHAO S J, et al. Omics analyses indicate the routes of lignin related metabolites regulated by trypsin during storage of pitaya () [J]. Genomics, 2021, 113(6): 3681– 3695. doi: 10.1016/j.ygeno.2021.08.005.

[40] V?GELI B, ENGILBERGE S, GIRARD E, et al. Archaeal aceto- acetyl-CoAthiolase/HMG-CoA synthase complex channels the inter- mediate via a fused CoA-binding site [J]. Proc Natl Acad Sci USA, 2018, 115(13): 3380–3385. doi: 10.1073/pnas.1718649115.

[41] ZHENG T, GUAN L B, YU K, et al. Expressional diversity of grape- vine 3-hydroxy-3-methylglutaryl-CoAreductase () in different grapes genotypes [J]. BMC Plant Biol, 2021, 21(1): 279. doi: 10.1186/s 12870-021-03073-8.

Common Pattern in Response to Pathogenic Fungal Stress of Tea Plants Based on Meta-analysis

LIAN Lingli1a, CHEN Qiang1b, ZHOU Ying1a, FU Jing1a, LI Wanying1a, WEI Rifeng1b, LIU Wei2*

(1a. College of Life Sciences; 1b. College of Horticulture, Fujian Agriculture and Forestry University,Fuzhou 350002, China; 2. Ningde Normal University, Ningde 352100, Fujian, China)

To explore the common response mode and disease resistance mechanism of tea plants () to pathogenic stress, bioinformatics methods were used to extract, integrate and function enrich of multiple sets of RNA-seq data, and the main regulatory molecules and protein interaction modules were analyzed by combining various tools and database resources. This results showed that the expression of cytochrome P450 family members in tea plant was significantly up-regulated under the fungal pathogen stress. The metabolic processes of steroid and hormone, and phenylpropanoid synthesis pathway were activated, and the biological processes, such as mitotic cell cycle regulation, DNA methylation and photosynthesis pathway were inhibited. The major regulatory molecules, such asandtranscription factors, theandfamily of kinases were mainly up-regulated. The differentially expressed protein interaction modules showed that the modules involved in mitotic cycle regulation, microtubule motion-based, starch and sucrose metabolism, cell wall polysaccharide synthesis, photosynthesis, flavonoid metabolism were down-regulated, while lignin synthesis and terpenoid biosynthesis were up-regulated. There may be interactions between modules. The key genes in lignin and terpenoid synthesis pathways activated by pathogen stress included ferulic acid-5-hydroxylase gene (), peroxidase gene () and terpenoid synthase gene. Cytochromegene might play a key role in fungus stress of tea plants. Enhancing the synthesis of lignin and terpenoids, and weakening photosynthesis might be the core modes of tea plants responding to fungus stress.

Tea plant; Fungal pathogen; Pathogenic fungi stress; Meta-analysis

10.11926/jtsb.4628

2022-02-24

2022-03-27

福建農(nóng)林大學(xué)科技創(chuàng)新專項基金項目(KFA20047A, KFA20143A)；福建農(nóng)林大學(xué)大學(xué)生創(chuàng)新項目(202110389107)資助

This work was supported by the Project for Science and Technology Innovation in Fujian Agricultre and Forestry University (Grant No. KFA20047A, KFA20143A), and the Project for Undergraduate Innovation in Fujian Agricultre and Forestry University (Grant No. 202110389107).

連玲麗(1979年生)，女，博士，副教授，研究方向為逆境生物信息學(xué)。E-mail: lianll2002@163.com

E-mail: liuwei0593@126.com