基于WGCNA鑒定茶樹響應草甘膦相關的基因共表達模塊

郭永春，王鵬杰，金珊，侯炳豪，王淑燕，趙峰，葉乃興

基于WGCNA鑒定茶樹響應草甘膦相關的基因共表達模塊

郭永春1，王鵬杰1，金珊1，侯炳豪1，王淑燕1，趙峰2，葉乃興1

1福建農林大學園藝學院/茶學福建省高校重點實驗室，福州 350002；2福建中醫藥大學藥學院，福州 350122

【】分析茶樹響應草甘膦相關的基因表達規律和調控途徑，在轉錄水平上探究草甘膦對茶樹的作用，確定茶樹響應草甘膦的關鍵基因。以茶樹品種‘金觀音’為試驗材料，將推薦使用濃度的草甘膦施于茶樹土壤基質表面，經0、0.25、1、3和7 d后，取葉片進行轉錄組測序，并測定莽草酸含量。利用WGCNA方法聯合分析轉錄組和莽草酸含量數據，鑒定與草甘膦響應相關的共表達基因模塊，篩選關鍵調控基因。茶樹葉片中的莽草酸含量在草甘膦處理后0.25、1和3 d降低，而在7 d時大量積累（為未處理的6.99倍）。從表達譜數據中篩選得到12 568個差異表達基因（DEGs），草甘膦處理不同時間點與未處理數據比對的DEGs均顯著富集在苯丙烷、類黃酮生物合成及植物激素信號轉導途徑；此外，草甘膦處理分別誘導茶樹莽草酸代謝及其下游苯丙烷、類黃酮生物合成和激素信號轉導途徑相關的24、52、31和69個基因差異表達。通過加權基因共表達網絡（WGCNA）方法鑒定得到19個基因模塊，將轉錄組與莽草酸含量數據相關聯，篩選到兩個與草甘膦響應高度相關的關鍵基因模塊，分別包含2 024和2 305個基因。選取關鍵模塊中連通度最高的前50個基因進行共表達分析，獲得6個關鍵調控基因，包括2個抗性基因（和）、1個耐藥性基因（）、1個離子轉運基因（）、1個膜轉運基因（）和1個轉錄因子（）。草甘膦通過干擾茶樹葉片中莽草酸代謝，影響其下游代謝途徑苯丙烷、類黃酮生物合成及激素信號轉導途徑的基因轉錄。此外，研究還鑒定到2個草甘膦響應密切相關的共表達模塊，發現、、和等多個潛在候選基因和轉錄因子在茶樹抵御草甘膦逆境中發揮重要作用。

茶樹；草甘膦；莽草酸；轉錄組學；WGCNA

0 引言

【研究意義】茶樹[(L.) O. Kuntze]是我國重要的經濟作物，其芽葉可制成風味獨特的成品茶，具有很大的經濟、健康和文化價值[1]。草甘膦是一種最常見的非選擇性除草劑，在殺滅雜草的同時，會被茶樹根部吸收并轉運至全株[2]。草甘膦在茶樹體內的殘留時間較長，并可能造成茶葉產品中草甘膦的殘留量超標，影響茶葉生產安全[3-4]。現階段，我國福建、江西、安徽的部分茶園已禁用草甘膦，但因其高效廉價，目前在實際生產中完全停止使用草甘膦較為困難。中國農藥信息網數據庫顯示，草甘膦仍是我國茶園主要登記使用的除草劑。因此，深入探究草甘膦對茶樹的影響，對保證茶葉質量安全性有重要意義。【前人研究進展】草甘膦主要通過抑制植物莽草酸代謝途徑中的關鍵酶活性，阻止莽草酸向芳香氨基酸（色氨酸、酪氨酸和苯丙氨酸）轉化，導致植物體內代謝紊亂[5]。研究證實草甘膦的使用對植物次生代謝物的生物合成存在不利作用，例如，JIANG等[6]研究表明草甘膦除草劑通過抑制大豆莽草酸代謝途徑，抑制了大豆頂芽中色氨酸和類黃酮生物合成相關基因的轉錄；RAINIO等[7]通過對馬鈴薯的土壤基質進行草甘膦處理，發現植株內的糖苷生物堿濃度下降；MALALGODA等[8]發現使用草甘膦除草劑干擾了小麥體內正常碳代謝，從而對小麥的氨基酸代謝產生負面影響。已有研究表明，草甘膦經茶樹根部吸收并富集至葉片，是茶葉中草甘膦殘留的主要來源[4,9-11]。此外，筆者課題組前期研究發現，草甘膦的使用不易使茶樹出現藥害表征，但可顯著改變葉片中游離氨基酸、兒茶素和生物堿類化合物的含量[12]。【本研究切入點】目前，草甘膦處理茶樹不同時期對葉片基因表達及相關代謝通路的影響有待研究，尤其是基于加權基因共表達網絡（weighted gene co-expression network analysis，WGCNA）的茶樹響應草甘膦除草劑的基因共表達模塊尚未被鑒定。【擬解決的關鍵問題】本研究對草甘膦處理5個時間梯度（第0、0.25、1、3和7天）的茶樹葉片進行轉錄組測序，分析其基因表達水平和調控途徑，并利用WGCNA構建共表達基因模塊，與莽草酸數據關聯分析，挖掘出與草甘膦誘導相關的關鍵模塊，進而確定模塊中與抵御草甘膦逆境相關的核心基因。

1 材料與方法

1.1 材料處理

2019年12月，將一年生‘金觀音’茶樹（Jin-guanyin）植于體積為2 L的盆缽進行適應性培養，栽培基質為泥炭土（10—30 mm）。于2020年8月選取株高約30 cm、樹幅15—20 cm的植株，將5 g·L-1（推薦施用濃度）的草甘膦異丙胺鹽在土壤基質表面均勻噴施0.3 L，分別在第0（噴施前）、0.25、1、3和7天采集茶樹嫩梢第二葉，液氮速凍，-80℃保存。每個時期取3個生物學重復。

1.2 莽草酸含量測定方法

將參試樣品分別在液氮中研磨成粉，稱取0.5000 g置于4 mL棕色容量瓶中，加入2 mL去離子水，振蕩混勻300 s。60℃水浴超聲提取4 h，4 000 r/min離心10 min，取上清液1 mL過0.45 μm濾膜待測。采用高效液相色譜儀（waters E2695）測定樣品中的莽草酸含量，色譜條件：C18色譜柱（250 mm×4.6 mm，0.5 μm）；柱溫：25℃；波長：213 nm；流動相：甲醇﹕1%磷酸水溶液=3﹕97；流速：0.80 mL·min-1；進樣量：10 μL。

1.3 轉錄組測序的獲取與分析

采用Plant RNA Purification Reagent（Invitrogen）提取各樣本的總RNA，分別使用2100 Bioanalyser（Agilent）和ND-2000（NanoDrop Technologies）檢測RNA質量（RIN≥6.5，OD260/280=1.8—2.2，OD260/230≥2.0），瓊脂糖凝膠電泳檢測RNA完整性。按照Illumina TruSeqTMRNA sample preparation Kit方法構建轉錄組文庫，通過上海美吉生物醫藥科技有限公司在Illumina Novaseq 6000測序平臺上對文庫片段進行雙末端測序。通過SeqPrep（https://github.com/jstjohn/ SeqPrep）和Sickle（https://github.com/najoshi/sickle）對原始數據進行質控，從而得到高質量的質控數據。利用TopHat 2.1.1軟件將質控數據比對到‘黃棪’茶樹基因組[13]，基因組從中國核酸數據庫網站（https:// bigd.big.ac.cn/gsa/index.jsp）下載，GSA：CRA003208。通過Cufflinks軟件組裝并得到本次與原有注釋的差異信息[14]，使用RESM軟件根據每百萬讀轉錄量（TPM）法計算基因表達水平[15]。

1.4 差異表達基因和富集分析

通過EdgeR包（http://www.bioconductor.org/ packages/2.12/bioc/html/edgeR.html）計算篩選差異表達基因（differentially expressed genes，DEGs）[16]，差異基因界定為：|log2FoldChange|＞1，-value＜0.05。差異倍數（FoldChange）表示兩樣品組間表達量的比值。采用Goatools和KOBAS軟件對差異表達基因進行GO和KEGG功能富集分析，當經過校正的值（-value）＜0.05時，認為此GO功能和KEGG pathway功能存在顯著富集情況[17]。使用Tbtools軟件制作熱圖對基因表達量進行可視化[18]。

1.5 共表達模塊構建與分析

利用RESM軟件對基因表達數據進行背景校正和標準化，過濾表達量過低和變異系數較小的基因[14]。利用R軟件（R version 3.4.4）和WGCNA（R version 1.6.6）包構建基因共表達網絡[19]，根據過濾數據（平均表達水平1，變異系數0.1）識別與代謝物高度相關的基因模塊。過濾后的18 976個基因和3個代謝物的豐度通過計算Pearsons相關性構建共表達網絡。使用Cytoscape 3.4.0對核心的共表達模塊進行可視化[20]。

1.6 轉錄組數據驗證

采用Primer3plus在線網站（http://www.bioinformatics. nl/cgi-bin/primer3plus/primer3plus.cgi）設計引物（附表1），然后利用全式金Easyscript One-step gDNA Removal and cDNA synthesis superMix試劑盒（北京全式金生物技術有限公司）將RNA合成cDNA。參照王鵬杰等[21]的方法，選用茶樹（登錄號GE651107）作為內參基因[22]，按照Transstart?Tip Green qPCR superMix試劑盒（北京全式金生物技術有限公司）的使用說明在CFX96 Touch熒光定量PCR儀（伯樂生命醫學產品（上海）有限公司）上進行qRT-PCR分析，反應程序：94℃ 30 s；94℃ 5 s，60℃ 30 s，40個循環，每個時間點設3次生物學重復。用2-ΔΔCT法計算基因相對表達水平[23]，GraphPad Prism 5軟件制作柱狀圖，并通過SPSS 17.0進行差異顯著性分析（＜0.05）。

2 結果

2.1 草甘膦處理后茶樹葉片表型及莽草酸積累量

如圖1所示，本試驗處理下茶樹嫩梢葉片無明顯藥害特征。圖2顯示，葉片中莽草酸含量在草甘膦處理0.25 d時顯著下降，而后升逐漸升高，草甘膦處理7 d時莽草酸含量最高，達到處理前的6.99倍。

圖1 草甘膦處理后茶樹葉片的表型

不同小寫字母表示差異顯著（P＜0.05）

2.2 茶樹草甘膦處理的轉錄組分析

2.2.1 測序質量及KEGG富集分析 轉錄組測序共產生713 187 436個原始序列，在去除測序接頭污染和低質量的原始序列后，共獲得高質量的706 923 928個過濾序列，共產生104 531 971 833個過濾堿基，過濾序列質量值大于30的堿基所占的百分比（Q30%）在93.69%以上，將各樣品的過濾序列與茶樹參考基因組比對，比對結果顯示各樣本的比對率為87.69%—89.79%，序列主要分布在參考基因組外顯子區域。轉錄組測序結果與茶樹參考基因組的比對率高，測序質量好，可用于后續分析（表1）。

基于所有樣本中的TPM分布情況進行皮爾森相關性分析，結果表明生物學重復樣本之間具有很強的相關性，后續分析的結果也更可信（圖3-A）。將草甘膦處理不同時間點的數據與未處理的數據比對，共篩選出12 568個DEGs（圖3-B）。0.25 d與0 d比對篩選到5 542個DEGs（上調3 319個，下調2 223個）；1 d與0 d比對篩選到2 201個DEGs（上調1 308個，下調893個）；3 d與0 d比對篩選到2 481個DEGs（上調1 491個，下調990個）；7 d與0 d比對篩選到2 344個DEGs（上調1 267個，下調1 077個）。其中0 d vs 0.25 d 、0 d vs 1 d、0 d vs 3 d和0 d vs 7 d比較組中分別有2 971、409、796和879個基因是特有的DEGs；有242個基因在4個比較組中都呈現差異表達，這些DEGs可在草甘膦處理7 d時持續發揮作用。

將4個比較組的DEGs進行生物學代謝途徑富集分析，結果僅顯示富集程度前10的代謝途徑（圖4）。與碳水化合物代謝相關的乙醛酸和二羧酸代謝（map00630）僅在“0 d vs 0.25 d”顯著富集；與能量代謝相關的硫代謝（map00920）和氮代謝（map00910）僅在“0 d vs 1 d”顯著富集；與脂質代謝相關的脂肪酸伸長率（map00062）和-亞麻酸代謝（map00592）僅在“0 d vs 3 d”富集；單萜類生物合成及與氨基酸代謝相關的谷胱甘肽代謝（map00480）僅在“0 d vs 7 d”顯著富集。值得關注的是，4個比較組的DEGs均顯著富集在苯丙烷和類黃酮等次生代謝物合成（map00940和map00941）及植物激素信號轉導（map04075）途徑。

表1 參試樣品的轉錄組數據的質量

A：相關性。B：左邊水平柱狀圖表示各集合的DEGs。中間矩陣中，單個點表示某個集合特有的元素，點和點之間的連線表示不同集合特有的交集，豎直柱狀圖中則分別表示對應交集的DEGs

縱軸表示多個基因集中富集到的KEGG pathway，橫軸表示Rich factor（Rich factor越大，表示富集的程度越大），點的大小表示此pathway中的基因個數，點的顏色對應于不同的P值范圍

2.2.2 苯丙烷、類黃酮生物合成及莽草酸代謝相關的DEGs分析 茶樹在響應草甘膦處理過程中，與苯丙烷、類黃酮生物合成相關的基因在0.25、1、3和7 d均受誘導而差異表達。圖5展示了苯丙烷、類黃酮生物合成及莽草酸代謝的通路，以及分別顯著富集在各通路上DEGs（24、52和31個）的轉錄水平。、、、、、和等莽草酸代謝相關的10個基因在處理過程中較高表達（13.85＜TPM＜179.64）。其中，（）、（）、（）和（）在處理各時間點均上調表達，其余6個基因呈現先下降后升高的趨勢。、、、、、、、和等苯丙烷生物合成相關的14個基因較高表達（5.28＜TPM＜353.77），這14個基因在0.25 d下調表達，隨后上調表達。、、、、、、和等類黃酮生物合成途徑中11個DEGs在處理過程較高表達（4.72＜TPM＜365.41），并呈現先下調后上調的趨勢。

通路圖中的紅色邊框表示各通路上注釋到的DEGs；熱圖表示所有與通路相關DEGs的表達水平，由3個重復樣本計算得出的平均值進行log2轉換生成。下同

2.2.3 植物激素信號轉導途徑相關的DEGs分析 草甘膦處理過程中，植物激素信號轉導途徑中生長素（Auxin）、細胞分裂素（Cytokinine）、脫落酸（Abscisic acid）、赤霉素（Gibberellin）、乙烯（Ethylene）、油菜素內酯（Brassinosteroids）、茉莉酸（Jasmonic acid）和水楊酸（Salicylic acid）等重要防御相關的植物激素信號通路均受到影響（圖6-A）。69個基因表達具有差異，并顯著富集于植物激素信號轉導途徑，14個DEGs（，2個；，1個；，1個；，2個；，1個；，1個；，1個；，1個；，1個；，1個；，1個；，1個）在處理過程中較高水平表達（圖6-B）。其中，乙烯信號途徑下游基因（）在處理不同時間點均上調表達（最高上調3.08倍），且表達水平較高（87.88＜TPM＜271.03）；水楊酸信號途徑下游基因（）在0.25 d顯著下調（200倍），而后上調表達（最高上調4.22倍）；脫落酸信號途徑下游基因（）在0.25 d下調表達（5.19倍），而后上調表達（最高上調1.71倍）。

2.3 茶樹草甘膦處理的WGCNA分析

2.3.1 基因共表達模塊構建 經過濾TPM＜1的基因，18 976個差異基因用于構建加權基因共表達網絡。通過計算15個樣本表達水平的相關系數聚類分析，樣本間聚類良好，未出現離群樣本（圖7-A）。依據無尺度容適曲線處于平滑處，設定冪指數加權值即軟閾值為14（圖7-B）。

根據差異基因的TPM值進行相關度分析并聚類，相關度較高的基因被分配到同一個模塊中，圖7-C中聚類樹的不同分支代表不同的基因共表達模塊，不同顏色代表不同的模塊，共劃分為19個模塊，顏色為灰色代表沒有劃分到其他模塊的基因。

2.3.2 基因共表達模塊篩選 圖8展示了19個模塊各包含的基因個數及其與莽草酸含量的相關性。不同模塊間包含的基因個數差異較大，其中，Blue模塊的基因個數最多，為3 437個；Royalblue模塊的基因個數最少，為60個。相關性系數（絕對值）較大且顯著性值較小的模塊與表型高度相關，鑒定到與莽草酸極顯著正相關的為Green模塊（= 0.864，=0.0000329），其次為Brown模塊（=0.771，=0.00296）；極顯著負相關的為Tan模塊（=-0.875，=0.0000195），其次為Blue模塊（=-0.754，= 0.00117）。

為了解各個共表達模塊中基因的生物學功能，對上述各個基因模塊進行KEGG富集分析（＜0.05），發現Green和Brown模塊顯著富集到酪氨酸代謝（map00350）、苯丙氨酸代謝（map00360）等芳香氨基酸代謝和花青素的生物合成（map00942）、類黃酮生物合成（map00941）、苯丙烷生物合成（map00940）等次生代謝物合成通路，根據功能關聯原則可知這2個基因模塊與草甘膦響應相關（圖9）。

2.3.3 共表達網絡可視化分析 選取Green和Brown模塊內連通度最高的前50個基因作為模塊的核心基因，并對核心基因的互作網絡利用Cytoscape軟件進行可視化（圖10）。在Green模塊中的基因（7 d上調）連通度最高的3個依次是和，是一個植物抗性基因，與614個基因相連；是一種抗病基因，與612個基因相連；是一個是離子轉運基因，與608個基因相連。在該網絡中發現6個轉錄因子（，3個；和，1個）與植物防御反應相關。在Brown模塊中連通度最高的3個基因依次為、和，是一種多效性耐藥性基因，與500個基因相連；是乙烯響應因子，與494個基因相連；是一種膜轉運蛋白基因，與480個基因相連。轉錄因子與植物防御反應相關。

2.4 qRT-PCR驗證

為了驗證RNA-Seq數據的可靠性，采用qRT-PCR方法檢測15個DEGs的表達水平，包括6個莽草酸代謝相關的DEGs、4個共表達模塊篩選到的關鍵基因和5個隨機挑選的DEGs（圖11）。結果表明，qRT-PCR和RNA-Seq的基因表達變化趨勢基本一致，二者數據結果呈正相關，表明本研究中基于轉錄組數據得到的DEGs是可信的。

A：DEGs富集到的植物激素信號通路圖；B：熱圖表示與植物激素信號轉導相關DEGs的表達水平

A：樣本聚類樹；B：無尺度容適曲線及平均連通度曲線；C：基因聚類樹及模塊切割，基因聚類樹的每一個分支對應一個模塊

3 討論

3.1 草甘膦能夠干擾茶樹葉片莽草酸代謝

莽草酸作為草甘膦抑制莽草酸代謝的積累產物，其含量是最能反映草甘膦毒性的指標之一[24]。在草甘膦處理3 d內，茶樹葉片中的莽草酸含量低于處理前，而在7 d時莽草酸含量大量積累，表明茶樹莽草酸代謝在試驗處理3 d后明顯受到了草甘膦的抑制。前人研究表明，莽草酸途徑上的關鍵酶基因和在抵御茶樹逆境脅迫中具有重要作用[25-26]，此外，作為草甘膦脅迫的作用靶點酶基因，已被廣泛用于轉基因抗草甘膦作物的開發[27]。本研究中（）和（）在處理不同時間點的表達量上調表達，說明它們可能參與茶樹對草甘膦脅迫的抵御反應，在草甘膦處理前期保護茶樹以免莽草酸積累中毒。

3.2 草甘膦影響茶樹苯丙烷、類黃酮生物合成及激素信號轉導途徑相關基因的轉錄

草甘膦干擾植物莽草酸代謝的同時，也會影響以芳香族氨基酸為前體物質的次生代謝物（如植物激素、類黃酮）合成[28]。ZHAI等[29]研究表明草甘膦誘導了水稻谷胱甘肽代謝、介導激素信號通路相關基因的差異表達。MESNAGE等[30]認為草甘膦能夠干擾土壤真菌的蛋白合成、次生代謝及應激反應，導致土壤質量下降。課題組前期發現，草甘膦處理顯著改變了茶樹葉片中游離氨基酸、兒茶素和生物堿總量及組分構成[12]。本研究中草甘膦處理不同時期（0.25、1、3和7 d）與處理前（0 d）比對的差異基因均顯著富集于苯丙烷、類黃酮生物合成途徑及植物激素信號轉導途徑。說明草甘膦干擾茶樹莽草酸代謝的同時，也會誘導下游代謝途徑苯丙烷、類黃酮生物合成及激素信號轉導相關基因的差異表達。此外，苯丙烷代謝途徑是參與植物防御反應的主要次級代謝通路之一，苯丙烷途徑相關基因大量表達有利于加速下游分支代謝產物（如類黃酮、花青素）的生成[31-32]。本研究中參與茶樹苯丙烷（、、和）及類黃酮（、、、和）生物合成途徑的多個關鍵基因在草甘膦處理0.25 d表達量下調，說明草甘膦在處理能夠迅速干擾茶樹苯丙烷、類黃酮化合物的生物合成，在前期對茶樹次生代謝產物的合成產生不利影響。

橫坐標代表不同表型，縱坐標代表不同模塊。圖中左側一列數字表示該模塊的基因數目，右側每組數據表示模塊與表型的相關性系數r值及顯著性P值（括號內）。紅色代表模塊與表型的相關性較大，藍色代表模塊與表型的相關性較小

3.3 基于WGCNA鑒定茶樹響應草甘膦相關的共表達基因

利用WGCNA分析法，可特異地篩選出與目標性狀相關的基因，并進行模塊化分類，得到具有高度生物學意義的共表達模塊[33]。莽草酸含量可直接反映植物對草甘膦的響應情況[24]。本研究利用WGCNA法將轉錄組與莽草酸含量數據相關聯，并對模塊進行KEGG功能富集分析，最終篩選到2個與草甘膦響應高度相關的關鍵模塊，分別是green模塊（2 024個基因）和brown模塊（2 305個基因）。WGCNA共表達分析中連通度高的基因通常起重要調控作用[33]。ZHAO等[34]認為草甘膦誘導了蜜蜂許多農藥解毒和抗性基因上調表達，以防御草甘膦對蜜蜂生長發育產生不利影響。本研究中、和在草甘膦處理7 d明顯上調表達（圖10），說明、和可能在草甘膦處理后期降低草甘膦對茶樹的不利作用。DOGRAMACI等[35]研究表明雜草通過細胞轉運、激素信號傳導和解毒機制等相互作用，以防御草甘膦誘導的莽草酸積累中毒。JIANG等[6]發現大豆轉運蛋白基因能夠降低草甘膦除草劑毒性作用。由此可知，、和等可能通過參與離子運輸、膜運輸和激素信號轉導等相互作用，在茶樹抵御草甘膦逆境過程中發揮重要作用。此外，（乙烯響應因子）是模塊中篩選出響應草甘膦的關鍵基因之一，本研究中乙烯信號途徑下游基因（）在處理不同時間點表達水平較高（87.88＜TPM＜271.03），并上調表達。因此，在茶樹響應草甘膦脅迫的激素信號轉導中，乙烯可能起重要作用。

縱軸表示pathway名稱，橫軸表示Rich factor；點的大小表示此pathway中基因個數多少，點的顏色對應于不同的P值范圍。圖中僅展示P-value＜0.05的前提下富集程度在前20的KEGG Pathway

三角形節點為連通性排名前3的基因，黑色邊框節點為轉錄因子基因

熱圖分別表示RNA-Seq和qRT-PCR的基因表達水平；兩個熱圖之間的值代表每個基因qRT-PCR和RNA-seq值的相關系數

4 結論

本研究利用草甘膦處理后茶樹葉片的RNA-seq及莽草酸表型數據，發現草甘膦能夠干擾茶樹葉片中莽草酸代謝，并誘導其下游苯丙烷、類黃酮生物合成及植物激素信號轉導途徑相關基因的差異表達，而乙烯可能在茶樹響應草甘膦脅迫的激素信號轉導中起重要作用。通過WGCNA方法發現2個與草甘膦響應高度相關的關鍵基因模塊，共表達分析發現6個關鍵調控基因（、、和）在茶樹抵御草甘膦逆境過程中發揮重要作用。

[1] 王鵬杰, 岳川, 陳笛, 鄭玉成, 鄭知臨, 林浥, 楊江帆, 葉乃興. 茶樹CsWRKY6、CsWRKY31和CsWRKY48基因的分離及表達分析. 浙江大學學報(農業與生命科學版), 2019, 45(1): 30-38.

WANG P J, YUE C, CHEN D, ZHENG Y C, ZHENG Z L, LIN Y, YANG J F, YE N X. Isolation and expression analysis of CsWRKY6, CsWRKY31 and CsWRKY48 genes in tea plant. Journal of Zhejiang University (Agriculture and Life Sciences), 2019, 45(1): 30-38. (in Chinese)

[2] 唐杏燕, 邵增瑯, 楊路成, 裴少芬, 岳鵬翔, 王曉霞. 茶園中草甘膦在靶標雜草和非靶標茶樹中的吸收、轉運、分布和代謝. 食品安全質量檢測學報, 2018, 9(18): 140-145.

Tang X Y, Shao Z L, Yang L C, PEI S F, YUE P X, WANG X X. Uptake, translocation, distribution and metabolism of glyphosate in target weeds and non-target tea trees in tea garden. Journal of Food Safety Quality Inspection, 2018, 9(18): 140-145. (in Chinese)

[3] 楊亞琴, 馮書惠, 胡永建, 李圓圓, 王會鋒, 劉進璽, 鐘紅艦. 氣相色譜-質譜法測定綠茶中草甘膦和氨甲基膦酸殘留量. 茶葉科學, 2020, 40(1): 125-132.

Yang Y Q, Feng S H, Hu Y J, LI Y Y, WANG H F, LIU J X, ZHONG H J. Determination of glyphosate and aminomethylphosphonic acid residue in green tea by gas chromatography-mass spectrometry. Journal of Tea Science, 2020, 40(1): 125-132. (in Chinese)

[4] 郭永春, 陳金發, 趙峰, 王淑燕, 王鵬杰, 周鵬, 歐陽立群, 金珊, 葉乃興. 草甘膦及其代謝物氨甲基膦酸在茶樹體中的分布研究. 茶葉科學, 2020, 40(4): 510-518.

GUO Y C, CHEN J F, ZHAO F, WANG S Y, WANG P J, ZHOU P, OUYANG L Q, JIN S, YE N X. Study on the distribution of glyphosate and its metabolite aminomethylphosphonic acid inJournal of Tea Science, 2020, 40(4): 510-518. (in Chinese)

[5] 于惠林, 賈芳, 全宗華, 崔海蘭, 李香菊. 施用草甘膦對轉基因抗除草劑大豆田雜草防除、大豆安全性及雜草發生的影響. 中國農業科學, 2020, 53(6): 1166-1177.

YU H L, JIA F, QUAN Z H, CUI H L, LI X J. Effects of glyphosate on weed control, soybean safety and weed occurrence in transgenic herbicide-resistant soybean. Scientia Agricultura Sinica, 2020, 53(6): 1166-1177. (in Chinese)

[6] JIANG L X, JIN L G, GUO Y, TAO B, QIU L J. Glyphosate effects on the gene expression of the apical bud in soybean (). Biochemical and Biophysical Research Communications, 2013, 437(4): 544-549.

[7] RAINIO M J, MARGUS A, VIRTANEN V, LINDSTR?M L, SALMINEN J P, SAIKKONEN K, HELANDER M. Glyphosate- based herbicide has soil-mediated effects on potato glycoalkaloids and oxidative status of a potato pest. Chemosphere, 2020, 258: 127254.

[8] MALAGODA M, OHM J, HOWATT K A, GREEN A, SIMSEK S. Effects of pre-harvest glyphosate use on protein composition and shikimic acid accumulation in spring wheat. Food Chemistry, 2020, 332: 127422.

[9] Tong M M, Gao W J, Jiao W, Hou R Y. Uptake, translocation, metabolism, and distribution of glyphosate in nontarget tea plant (L.). Journal of Agricultural and Food Chemistry, 2017, 65(35): 7638-7646.

[10] 高萬君, 張永志, 童蒙蒙, 馬慧勤, 錢珊珊, 王天雨, 李葉云, 吳慧平, 侯如燕. 茶園常用除草劑田間藥效試驗與殘留動態. 茶葉科學, 2019, 39(5): 587-594.

Gao W J, Zhang Y Z, Tong M M, Ma H Q, Qian S S, Wang T Y, Li Y Y, Wu H P, Hou R Y. Weeds control effect and residues of several herbicides in tea gardens. Journal of Tea Science, 2019, 39(5): 587-594. (in Chinese)

[11] 高萬君, 李葉云, 侯如燕. 茶葉中草甘膦殘留現狀與對策. 中國茶葉, 2021, 43(4): 20-24.

Gao W J, Li Y Y, Hou R Y. Status and countermeasures of glyphosate residues in tea. China Tea, 2021, 43(4): 20-24. (in Chinese)

[12] 郭永春, 王淑燕, 王鵬杰, 陳金發, 周鵬, 歐陽立群, 趙峰, 葉乃興. 草甘膦對茶樹葉片主要生化成分的影響. 天然產物研究與開發, 2021, 33(3): 394-401.

GUO Y C, WANG S Y, WANG P J, CHEN J F, ZHOU P, OUYANG L Q, ZHAO F, YE N X. The effect of glyphosate on the main biochemical components of tea leaves. Natural Product Research and Development, 2021, 33(3): 394-401. (in Chinese)

[13] WANG P J, YU J X, JIN S, CHEN S, YUE C, WANG W L, GAO S L, CAO H L, ZHENG Y C, GU M Y, CHEN X J, SUN Y, GUO Y Q, YANG J F, ZHANG X T, YE N X. Genetic basis of high aroma and stress tolerance in the oolong tea cultivar genome. Horticulture Research, 2021, 8(1): 107.

[14] TRAPNELL C, ROBERTS A, GOFF L, PERTEA G, KIM D, KELLEY D R, PIMENTEL H, SALZBERG S L, RINN J L, PACHTER L. Differential gene and transcript expression analysis of RNA-seq experiments with TopHat and Cufflinks. Nature Protocols, 2012, 7(3): 562-578.

[15] LI B, DEWEY C N. RSEM: Accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics, 2011, 12: 323.

[16] ROBINSON M D, MCCARTHY D J, SMYTH G K. edgeR: a bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics, 2010, 26(1): 139-140.

[17] XIE C, MAO X Z, HUANG J J, DING Y, WU J M, DONG S, KONG L, GAO G, LI C Y, WEI L P. KOBAS 2.0: a web server for annotation and identification of enriched pathways and diseases. Nucleic Acids Research, 2011, 39: W316-W322.

[18] CHEN C, XIA R, CHEN H, XIA T. TBtools, a toolkit for biologists integrating various HTS-data handling tools with a user-friendly interface. BioRxiv, 2018.

[19] LANGFELDER P, HORVATH S. WGCNA: An R package for weighted correlation network analysis. BMC Bioinformatics, 2008, 9: 559.

[20] KOHL M, WIESE S, WARSCHEID B. Cytoscape: Software for visualization and analysis of biological networks. Methods in Molecular Biology, 2011, 696: 291-303.

[21] 王鵬杰, 曹紅利, 陳丹, 陳笛, 陳桂信, 楊江帆, 葉乃興. 茶樹脂肪酸去飽和酶家族基因的克隆與表達分析. 園藝學報, 2020, 47(6): 1141-1152.

WANG P J, CAO H L, CHEN D, CHEN D, CHEN G X, YANG J F, YE N X. Cloning and expression analysis of fatty acid desaturase family genes in. Acta Horticulturae Sinica, 2020, 47(6): 1141-1152. (in Chinese)

[22] WU Z J, TIAN C, JIANG Q A, LI X H, ZHUANG J. Selection of suitable reference genes for qRT-PCR normalization during leaf development and hormonal stimuli in tea plant (). Scientific Reports, 2016, 6: 19748.

[23] LIVAK K J, SCHMITTGEN T D. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) method. Methods, 2001, 25: 402-408.

[24] ZELAYA I A, ANDERSON J A H, OWEN M D K, LANDES R D. Evaluation of spectrophotometric and HPLC methods for shikimic acid determination in plants: Models in glyphosate-resistant and -susceptible crops. Journal of Agricultural and Food Chemistry, 2011, 59(6): 2202-2212.

[25] 郭永春, 王鵬杰, 谷夢雅, 王淑燕, 趙峰, 葉乃興. 茶樹5-烯醇式丙酮酰莽草酸-3-磷酸合成酶基因的克隆及表達. 應用與環境生物學報, 2020. https://doi.org/10.19675/j.cnki.1006-687X.2020.10031.

GUO Y C, WANG P J, GU M Y, WANG S Y, ZHAO F, YE N X. Cloning and expression of 5-enolpyruvylshikimate-3-phosphate synthase gene from tea plants. Chinese Journal of Applied and Environmental Biology, 2020. https://doi.org/10.19675/j.cnki.1006- 687X.2020.10031. (in Chinese)

[26] 孫平, 章國營, 向萍, 林金科, 賴鐘雄. 茶樹中莽草酸途徑DHD/SDH基因的表達調控. 應用與環境生物學報, 2018, 24(2): 322-327.

SUN P, ZHANG G Y, XIANG P, LIN J K, LAI Z X. Expression and regulation of the shikimic acid pathway gene DHD/SDH in tea plant (). Chinese Journal of Applied and Environmental Biology, 2018, 24(2): 322-327. (in Chinese)

[27] ACHARY V M M, SHERI V, MANNA M, PANDITI V, BORPHUKAN B, RAM B, AGARWAL A, FARTYAL D, TEOTIA D, MASAKAPALLI S K, AGRAWAL P K, REDDY M K. Overexpression of improved EPSPS gene results in field level glyphosate tolerance and higher grain yield in rice. Plant Biotechnology Journal, 2020, 18(12): 2504-2519.

[28] 陳延儒. 采后硝普鈉處理對蘋果果實品質、莽草酸和苯丙烷代謝途徑的影響[D]. 錦州: 渤海大學, 2019.

CHEN Y R. Effect of postharvest sodium nitroprusside treatment on the quality, shikimate and phenylpropanoid pathway of apple fruit [D]. Jinzhou: Bohai University, 2019. (in Chinese)

[29] ZHAI R R, YE S H, ZHU G F, LU Y T, YE J, YU F M, CHU Q R, ZHANG X M. Identification and integrated analysis of glyphosate stress-responsive microRNAs, lncRNAs, and mRNAs in rice using genome-wide high-throughput sequencing. BMC Genomics, 2020, 21(1): 238.

[30] MESNAGE R, OESTREICHER N, POIRIER F, NICOLAS V, BOURSIER C, VéLOT C. Transcriptome profiling of the fungusexposed to a commercial glyphosate-based herbicide under conditions of apparent herbicide tolerance. Environmental Research, 2020, 182: 109116.

[31] MENG J, WANG B, HE G, WANG Y, TANG X F, WANG S M, MA Y B, FU C X, CHAI G H, ZHOU G K. Metabolomics integrated with transcriptomics reveals redirection of the phenylpropanoids metabolic flux in. Journal of Agricultural and Food Chemistry, 2019, 67(11): 3284-3291.

[32] DONG N Q, LIN H X. Contribution of phenylpropanoid metabolism to plant development and plant-environment interactions. Journal of Integrative Plant Biology, 2021, 63(1): 180-209.

[33] 秦天元, 孫超, 畢真真, 梁文君, 李鵬程, 張俊蓮, 白江平. 基于WGCNA的馬鈴薯根系抗旱相關共表達模塊鑒定和核心基因發掘. 作物學報, 2020, 46(7): 1033-1051.

QIN T Y, SUN C, BI Z Z, LIANG W J, LI P C, ZHANG J L, BAI J P. Identification of drought-related co-expression modules and hub genes in potato roots based on WGCNA. Acta Agronomica Sinica, 2020, 46(7): 1033-1051. (in Chinese)

[34] ZHAO H, LI G L, GUO D Z, WANG Y, LIU Q X, GAO Z, WANG H F, LIU Z G, GUO X Q, XU B H. Transcriptomic and metabolomic landscape of the molecular effects of glyphosate commercial formulation onmellifera ligustica andcerana. The Science of the Total Environment, 2020, 744: 140819.

[35] DOGRAMACI M, ANDERSON J V, CHAO W S, HORVATH D P, HERNANDEZ A G, MIKEL M A, FOLEY M E. Foliar glyphosate treatment alters transcript and hormone profiles in crown buds of leafy spurge and induces dwarfed and bushy phenotypes throughout its perennial lifecycle. The Plant Genome, 2017, 10(3): 98.

附表1 轉錄組數據RT-qPCR 驗證引物序列

Table S1 Transcriptome data RT-qPCR validation primer sequences

Identification of Co-Expression Gene Related to Tea Plant Response to Glyphosate Based on WGCNA

1College of Horticulture, Fujian Agriculture and Forestry University/Key Laboratory of Tea Science in Universities of Fujian Province, Fuzhou 350002;2School of Pharmacy, Fujian University of Chinese Medicine, Fuzhou 350122

【】This study aimed at analyzing both expression patterns and regulatory pathways of tea plants in response to glyphosate stressing, which could revealed the effect of glyphosate herbicides on tea plants at transcriptional level and identify key genes of tea plants. 【】cv Jin-guanyin was applied as material plant. A recommended concentration of glyphosate was irrigated to test plants. The leave samples were collected at different time intervals (0, 0.25, 1, 3 and 7 d). The samples were sequenced by transcriptome, the content of shikimic acid was also quantified.The WGCNA method was used to jointly analyze transcriptome and shikimic acid content data, to identify co-expressed gene modules related to glyphosate response, and to screen out key regulatory genes. 【】The content of shikimic acid in tea leaves reduced gradually during first 3 days. However, it suddenly reached a peak on the 7th day (6.99 times compared with no glyphosate treated sample). A total of 12 568 differential expression genes (DEGs) were also identified, which mainly enriched in phenylpropane, flavonoid biosynthesis and plant hormone signal transduction pathways. In addition, the glyphosate treatment induced 24, 52, 31 and 69 genes respectively which related to shikimic acid metabolism, phenylpropane, flavonoid biosynthesis and hormone signal transduction pathways. A total of 19 modules were screened out by WGCNA method. The correlation analysis of transcriptome and shikimic acid content indicated two key modules, including 2 024 and 2 305 genes, respectively. The top 50 genes with the highest connectivity in the key modules were selected for co-expression analysis, and 6 key response genes were obtained, including 2 resistance genes (and), 1 drug resistance gene (), 1 ion transport gene (), 1 membrane transport gene (), and 1 transcription factor gene ().【】Glyphosate could affect downstream genes transcription of phenylpropane, flavonoid biosynthesis and hormone signal transduction pathways by interfering shikimic acid metabolism of tea plants. In addition, this study also identified two co-expression modules closely related to glyphosate response, and found that multiple potential candidate genes and transcription factors could resist glyphosate stress, such as,,,,and

; glyphosate; shikimic acid; transcriptome; WGCNA

10.3864/j.issn.0578-1752.2022.01.013

2021-03-12；

2021-05-10

福建省“2011協同創新中心”中國烏龍茶產業協同創新中心專項（閩教科[2015]75號）、福建農林大學科技創新專項基金（CXZX2017181）、福建農林大學茶產業鏈科技創新與服務體系建設項目（2020A01）、福建農林大學園藝學院優秀碩士學位論文資助基金（2019S01）

郭永春，E-mail：1062941682@qq.com。通信作者葉乃興，E-mail：ynxtea@126.com。通信作者趙峰，E-mail：zhaofeng0591@fjtcm.edu.cn

（責任編輯 趙伶俐）