基于WGCNA鑒定茶樹響應(yīng)草甘膦相關(guān)的基因共表達(dá)模塊

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

基于WGCNA鑒定茶樹響應(yīng)草甘膦相關(guān)的基因共表達(dá)模塊

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

1福建農(nóng)林大學(xué)園藝學(xué)院/茶學(xué)福建省高校重點(diǎn)實(shí)驗(yàn)室，福州 350002；2福建中醫(yī)藥大學(xué)藥學(xué)院，福州 350122

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

茶樹；草甘膦；莽草酸；轉(zhuǎn)錄組學(xué)；WGCNA

0 引言

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

1 材料與方法

1.1 材料處理

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

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；進(jìn)樣量：10 μL。

1.3 轉(zhuǎn)錄組測序的獲取與分析

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

1.4 差異表達(dá)基因和富集分析

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

1.5 共表達(dá)模塊構(gòu)建與分析

利用RESM軟件對基因表達(dá)數(shù)據(jù)進(jìn)行背景校正和標(biāo)準(zhǔn)化，過濾表達(dá)量過低和變異系數(shù)較小的基因[14]。利用R軟件（R version 3.4.4）和WGCNA（R version 1.6.6）包構(gòu)建基因共表達(dá)網(wǎng)絡(luò)[19]，根據(jù)過濾數(shù)據(jù)（平均表達(dá)水平1，變異系數(shù)0.1）識別與代謝物高度相關(guān)的基因模塊。過濾后的18 976個(gè)基因和3個(gè)代謝物的豐度通過計(jì)算Pearsons相關(guān)性構(gòu)建共表達(dá)網(wǎng)絡(luò)。使用Cytoscape 3.4.0對核心的共表達(dá)模塊進(jìn)行可視化[20]。

1.6 轉(zhuǎn)錄組數(shù)據(jù)驗(yàn)證

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

2 結(jié)果

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

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

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

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

2.2 茶樹草甘膦處理的轉(zhuǎn)錄組分析

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

基于所有樣本中的TPM分布情況進(jìn)行皮爾森相關(guān)性分析，結(jié)果表明生物學(xué)重復(fù)樣本之間具有很強(qiáng)的相關(guān)性，后續(xù)分析的結(jié)果也更可信（圖3-A）。將草甘膦處理不同時(shí)間點(diǎn)的數(shù)據(jù)與未處理的數(shù)據(jù)比對，共篩選出12 568個(gè)DEGs（圖3-B）。0.25 d與0 d比對篩選到5 542個(gè)DEGs（上調(diào)3 319個(gè)，下調(diào)2 223個(gè)）；1 d與0 d比對篩選到2 201個(gè)DEGs（上調(diào)1 308個(gè)，下調(diào)893個(gè)）；3 d與0 d比對篩選到2 481個(gè)DEGs（上調(diào)1 491個(gè)，下調(diào)990個(gè)）；7 d與0 d比對篩選到2 344個(gè)DEGs（上調(diào)1 267個(gè)，下調(diào)1 077個(gè)）。其中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個(gè)基因是特有的DEGs；有242個(gè)基因在4個(gè)比較組中都呈現(xiàn)差異表達(dá)，這些DEGs可在草甘膦處理7 d時(shí)持續(xù)發(fā)揮作用。

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

表1 參試樣品的轉(zhuǎn)錄組數(shù)據(jù)的質(zhì)量

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

縱軸表示多個(gè)基因集中富集到的KEGG pathway，橫軸表示Rich factor（Rich factor越大，表示富集的程度越大），點(diǎn)的大小表示此pathway中的基因個(gè)數(shù)，點(diǎn)的顏色對應(yīng)于不同的P值范圍

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

通路圖中的紅色邊框表示各通路上注釋到的DEGs；熱圖表示所有與通路相關(guān)DEGs的表達(dá)水平，由3個(gè)重復(fù)樣本計(jì)算得出的平均值進(jìn)行l(wèi)og2轉(zhuǎn)換生成。下同

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

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

2.3.1 基因共表達(dá)模塊構(gòu)建 經(jīng)過濾TPM＜1的基因，18 976個(gè)差異基因用于構(gòu)建加權(quán)基因共表達(dá)網(wǎng)絡(luò)。通過計(jì)算15個(gè)樣本表達(dá)水平的相關(guān)系數(shù)聚類分析，樣本間聚類良好，未出現(xiàn)離群樣本（圖7-A）。依據(jù)無尺度容適曲線處于平滑處，設(shè)定冪指數(shù)加權(quán)值即軟閾值為14（圖7-B）。

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

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

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

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

2.4 qRT-PCR驗(yàn)證

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

A：DEGs富集到的植物激素信號通路圖；B：熱圖表示與植物激素信號轉(zhuǎn)導(dǎo)相關(guān)DEGs的表達(dá)水平

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

3 討論

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

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

3.2 草甘膦影響茶樹苯丙烷、類黃酮生物合成及激素信號轉(zhuǎn)導(dǎo)途徑相關(guān)基因的轉(zhuǎn)錄

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

橫坐標(biāo)代表不同表型，縱坐標(biāo)代表不同模塊。圖中左側(cè)一列數(shù)字表示該模塊的基因數(shù)目，右側(cè)每組數(shù)據(jù)表示模塊與表型的相關(guān)性系數(shù)r值及顯著性P值（括號內(nèi)）。紅色代表模塊與表型的相關(guān)性較大，藍(lán)色代表模塊與表型的相關(guān)性較小

3.3 基于WGCNA鑒定茶樹響應(yīng)草甘膦相關(guān)的共表達(dá)基因

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

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

三角形節(jié)點(diǎn)為連通性排名前3的基因，黑色邊框節(jié)點(diǎn)為轉(zhuǎn)錄因子基因

熱圖分別表示RNA-Seq和qRT-PCR的基因表達(dá)水平；兩個(gè)熱圖之間的值代表每個(gè)基因qRT-PCR和RNA-seq值的相關(guān)系數(shù)

4 結(jié)論

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

[1] 王鵬杰, 岳川, 陳笛, 鄭玉成, 鄭知臨, 林浥, 楊江帆, 葉乃興. 茶樹CsWRKY6、CsWRKY31和CsWRKY48基因的分離及表達(dá)分析. 浙江大學(xué)學(xué)報(bào)(農(nóng)業(yè)與生命科學(xué)版), 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] 唐杏燕, 邵增瑯, 楊路成, 裴少芬, 岳鵬翔, 王曉霞. 茶園中草甘膦在靶標(biāo)雜草和非靶標(biāo)茶樹中的吸收、轉(zhuǎn)運(yùn)、分布和代謝. 食品安全質(zhì)量檢測學(xué)報(bào), 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] 楊亞琴, 馮書惠, 胡永建, 李圓圓, 王會鋒, 劉進(jìn)璽, 鐘紅艦. 氣相色譜-質(zhì)譜法測定綠茶中草甘膦和氨甲基膦酸殘留量. 茶葉科學(xué), 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] 郭永春, 陳金發(fā), 趙峰, 王淑燕, 王鵬杰, 周鵬, 歐陽立群, 金珊, 葉乃興. 草甘膦及其代謝物氨甲基膦酸在茶樹體中的分布研究. 茶葉科學(xué), 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] 于惠林, 賈芳, 全宗華, 崔海蘭, 李香菊. 施用草甘膦對轉(zhuǎn)基因抗除草劑大豆田雜草防除、大豆安全性及雜草發(fā)生的影響. 中國農(nóng)業(yè)科學(xué), 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] 高萬君, 張永志, 童蒙蒙, 馬慧勤, 錢珊珊, 王天雨, 李葉云, 吳慧平, 侯如燕. 茶園常用除草劑田間藥效試驗(yàn)與殘留動態(tài). 茶葉科學(xué), 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] 高萬君, 李葉云, 侯如燕. 茶葉中草甘膦殘留現(xiàn)狀與對策. 中國茶葉, 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] 郭永春, 王淑燕, 王鵬杰, 陳金發(fā), 周鵬, 歐陽立群, 趙峰, 葉乃興. 草甘膦對茶樹葉片主要生化成分的影響. 天然產(chǎn)物研究與開發(fā), 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] 王鵬杰, 曹紅利, 陳丹, 陳笛, 陳桂信, 楊江帆, 葉乃興. 茶樹脂肪酸去飽和酶家族基因的克隆與表達(dá)分析. 園藝學(xué)報(bào), 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-磷酸合成酶基因的克隆及表達(dá). 應(yīng)用與環(huán)境生物學(xué)報(bào), 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基因的表達(dá)調(diào)控. 應(yīng)用與環(huán)境生物學(xué)報(bào), 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] 陳延儒. 采后硝普鈉處理對蘋果果實(shí)品質(zhì)、莽草酸和苯丙烷代謝途徑的影響[D]. 錦州: 渤海大學(xué), 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的馬鈴薯根系抗旱相關(guān)共表達(dá)模塊鑒定和核心基因發(fā)掘. 作物學(xué)報(bào), 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 轉(zhuǎn)錄組數(shù)據(jù)RT-qPCR 驗(yàn)證引物序列

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協(xié)同創(chuàng)新中心”中國烏龍茶產(chǎn)業(yè)協(xié)同創(chuàng)新中心專項(xiàng)（閩教科[2015]75號）、福建農(nóng)林大學(xué)科技創(chuàng)新專項(xiàng)基金（CXZX2017181）、福建農(nóng)林大學(xué)茶產(chǎn)業(yè)鏈科技創(chuàng)新與服務(wù)體系建設(shè)項(xiàng)目（2020A01）、福建農(nóng)林大學(xué)園藝學(xué)院優(yōu)秀碩士學(xué)位論文資助基金（2019S01）

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

（責(zé)任編輯 趙伶俐）