絨毛白蠟（Fraxinus velutina Torr）韌皮部響應(yīng)白蠟窄吉?。ˋgrilus planipennis Fairmaire）危害的轉(zhuǎn)錄組變化

封潤霞，趙 婕，張?zhí)K芳，王建軍，魏建榮，劉建鳳＊

(1.河北大學(xué)生命科學(xué)學(xué)院，河北 保定 071000；2.中國林業(yè)科學(xué)研究院森林生態(tài)環(huán)境與保護(hù)研究所，北京 100091；3.遼寧省林業(yè)科學(xué)研究院，遼寧沈陽 110032)

絨毛白蠟為木犀科（Oleaceae）落葉喬木，原產(chǎn)于北美，現(xiàn)在我國華北、內(nèi)蒙古和長江中下游地區(qū)均有栽培，是重要的木材和觀賞物種。白蠟窄吉?。ˋgrilus planipennisFairmaire）屬鞘翅目吉丁總科吉丁科（Coleoptera:Buprestidae），是木犀科梣屬（Fraxinus）樹木的一種重要蛀干害蟲[1]，其幼蟲在樹木的韌皮部、形成層和木質(zhì)部淺層蛀食為害，嚴(yán)重危害我國北方地區(qū)引種栽植的各個白蠟樹種，如洋白蠟（Fraxinus pennsylvanicaMarsh.）和絨毛白蠟（Fraxinus velutinaTorr）[2]。白蠟窄吉丁幼蟲的危害極具隱蔽性，其在樹皮內(nèi)蛀食為害時不將碎木和蟲糞推出坑道，從樹皮外很難發(fā)現(xiàn)樹木受害。成蟲的飛翔能力較弱，新羽化的成蟲常常在同一棵樹上或附近樹上繼續(xù)產(chǎn)卵為害[3]。因此，成蟲一旦侵染樹干成功，危害會逐年加重，通常1～3 年即可導(dǎo)致樹木死亡[4]。由于白蠟窄吉丁幼蟲生活的隱蔽性，常規(guī)的防治方法很難奏效，所以目前對于白蠟窄吉丁的防治主要集中在成蟲期，一般采用誘捕器、粘蟲板、阻蟲網(wǎng)以及藥劑噴霧等。目前國內(nèi)針對白蠟吉丁幼蟲天敵的研究較多，主要是為利用天敵昆蟲的生物防治服務(wù)[5-6]。

植物在長期進(jìn)化過程中，形成了通過物理結(jié)構(gòu)和有毒次生代謝產(chǎn)物抵御植食性昆蟲危害的機(jī)制[7]。次級代謝產(chǎn)物的應(yīng)答策略可以通過分子水平上的變化來解析其產(chǎn)生的分子機(jī)制。目前，國內(nèi)外對于白蠟樹抗白蠟窄吉丁的相關(guān)研究相對較少，尤其是針對白蠟樹響應(yīng)害蟲侵害時的韌皮部組織中轉(zhuǎn)錄組水平變化的研究更是鮮見報道。鑒于此，作者對絨毛白蠟與白蠟窄吉丁幼蟲互作過程中韌皮部的分子表達(dá)譜進(jìn)行了研究。通過對健康與受害白蠟樹的韌皮部材料進(jìn)行轉(zhuǎn)錄組水平測序、組裝和注釋，獲得了寄主樹與幼蟲互作過程中被激活的各類轉(zhuǎn)錄因子家族，并分別對差異表達(dá)基因進(jìn)行GO (Gene Ontology)、KEGG(Kyoto Encyclopedia of Genes and Genomes)注釋和富集分析。研究結(jié)果將為絨毛白蠟樹抗蟲分子機(jī)制研究奠定理論基礎(chǔ)，同時也為豐富木犀科植物抗蟲相關(guān)基因庫，促進(jìn)白蠟抗逆分子育種和優(yōu)良品系培育提供重要科學(xué)依據(jù)。

1 材料與方法

1.1 試驗材料

2019 年9 月中旬在遼寧省凌海市（40°48′～441°26′ N，120°42′～121°45′ E）白蠟樹種植園內(nèi)，選取長勢一致樹齡10 年，直徑6～8 cm 的健康（Non-infested）與受害絨毛白蠟樹（infested，ASF）為實驗材料。選區(qū)部位：截取距離樹干基部1.5～2.0 m處的韌皮部，此區(qū)域約有白蠟窄吉丁幼蟲30 頭，在此部位周徑上分別采集3 塊10 cm×5 cm(長×寬)的韌皮部材料，并混合為一個重復(fù)，直接放入液氮中保存?zhèn)溆?。該試驗設(shè)置3 個重復(fù)。

1.2 RNA 抽提

采用TRIzol (Invitrogen) 法提取蜈蚣組織中的總RNA，并使用DNase I (TaKara)去除基因組DNA。分別采用2 100 Bioanalyser (Agilent)、ND-2000 (NanoDrop Technologies)方法檢測RNA 樣品的質(zhì)量，以保證使用合格的樣品（OD260/280=1.8～2.2,OD260/230 ≥ 2.0,RIN ≥ 6.5,28S:18S ≥1.0,>2μg）進(jìn)行轉(zhuǎn)錄組測序。

1.3 文庫構(gòu)建及RNA-Seq 測序

RNA 文庫的建立采用TruSeqTM RNA sample preparation Kit(Illumina,San Diego,CA)試劑盒。首先利用帶有Oligo(dT)的磁珠從5 ug 總RNA 中富集有poly-A 尾的mRNA。再加入fragmentation buffer，將mRNA 隨機(jī)斷裂成200 bp 左右的小片段。接著采用SuperScript double-stranded cDNA synthesis kit (Invitrogen,CA)試劑盒、加入六堿基隨機(jī)引物(Illumina)，以mRNA 為模板反轉(zhuǎn)合成一鏈cDNA，隨后進(jìn)行二鏈合成，形成穩(wěn)定的雙鏈結(jié)構(gòu)。雙鏈的cDNA 結(jié)構(gòu)為粘性末端，加入End Repair Mix 將其補(bǔ)成平末端，隨后在3′末端加上一個A 堿基，用于連接Y 字形的接頭，具體步驟參見說明書。cDNA 經(jīng)過PCR 富集后，利用2%瓊脂糖膠回收200～300 bp 的條帶。經(jīng)TBS380（Picogreen）定量后，文庫使用Illumina HiSeq Xten/NovaSeq 6000 測序平臺進(jìn)行高通量測序，測序讀長為PE 150。

1.4 數(shù)據(jù)的處理與篩選

序列數(shù)據(jù)(raw data 或raw reads)由Illumina HiSeq Xten/NovaSeq 6 000 測序得到的原始圖像數(shù)據(jù)經(jīng)base calling 轉(zhuǎn)化形成FASTQ 格式的文件，該序列數(shù)據(jù)包含reads 的序列以及堿基的測序質(zhì)量，經(jīng)過去除和過濾含adapter、N 過多或低質(zhì)量堿基的dirty raw reads 后獲得高質(zhì)量的clean reads，用于后續(xù)的信息分析。用FPKM (Fragments Per Kilobase of transcript per Million mapped reads)值衡量計算表達(dá)量。根據(jù)FDR <0.05 且|log 2FC| >1 篩選出差異表達(dá)基因(differentially expressed genes,DEGs)。

2 結(jié)果與分析

2.1 幼蟲受害后白蠟韌皮部轉(zhuǎn)錄組數(shù)據(jù)質(zhì)量分析

分別取健康（Non-infested 1-3）和受害（ASF 1-3）的白蠟樹韌皮部，對樣品文庫進(jìn)行了RNA-Seq測序，分別獲得Raw Reads 49 287 036、44 918 932、50 780 898、49 722 100、45 082 112 和 45 802 214條，去除有接頭(Adapter related)、N 含量超過10%(Containing N)以及低質(zhì)量(Low quality)的reads，最終得到(Clean reads)分別為48 778 052、44 592 026、50 393 962、49 349 330、44 708 080 和45 432 400 條；去除Phred 值小于20 的堿基片段，樣本高于Q20測序結(jié)果比例分別為98.11%、98.34%、98.3%、98.29%、98.3%和98.25%；堿基G 和C 的數(shù)量總和占總的堿基數(shù)量的百分比分別為44.75%、43.58%、43.3%、43.57%、43.4%和43.37% (表1)。結(jié)果表明測序質(zhì)量合格，可進(jìn)行下一步的生物學(xué)分析。

2.2 健康與受害的韌皮部差異表達(dá)基因DEGs 的篩選

對健康與受害白蠟樹韌皮部進(jìn)行樣品間相關(guān)性及基因表達(dá)差異分析，結(jié)果顯示，健康與受害白蠟樹韌皮部之間相關(guān)性較低，均分布在0.7～0.8 之間，而各自生物學(xué)重復(fù)間的相關(guān)性相對較高，均高于0.9，說明我們前期分別采集健康與受害白蠟樹韌皮部并設(shè)置3 組重復(fù)的實驗設(shè)計較為合理（圖1A）。進(jìn)一步對健康與受害白蠟樹韌皮部的表達(dá)差異基因通過火山圖進(jìn)行分析，顯著上調(diào)的基因以紅色表示，顯著下調(diào)的基因以綠色表示，無顯著性差異的基因以灰色表示(圖1B)。共得到了3 388 個差異表達(dá)基因(DEGs)，其中受害相對于健康白蠟樹韌皮部表達(dá)下調(diào)的DEGs 有1 247 個，表達(dá)上調(diào)的DEGs 有2 141 個，上調(diào)的基因數(shù)明顯高于下調(diào)基因數(shù)。結(jié)果表明，蟲害脅迫下白蠟樹韌皮部的大部分基因受到激活，少數(shù)基因受到抑制。

表1 轉(zhuǎn)錄組優(yōu)化組裝結(jié)果評估Table 1 Evaluation of transcriptome optimized assembly results

圖1 健康與受害樣本間相關(guān)性熱圖（A）與差異基因火山圖（B）Fig.1 Heat map of correlation between health and hazard samples（A）and differential gene volcano map（B）

2.3 健康與受害的韌皮部差異表達(dá)基因GO 分析

為明確健康與受害絨毛白蠟樹韌皮部差異表達(dá)基因的生物學(xué)功能，對這些基因進(jìn)行GO 功能注釋分析，將差異表達(dá)基因分為生物學(xué)過程(biological process)、細(xì)胞組成(cellular component)和分子功能(molecular function)三大類。差異表達(dá)基因可主要歸納于9 個生物學(xué)過程、7 個細(xì)胞組分以及4 個分子功能(圖2)。在生物學(xué)過程中以細(xì)胞過程（cellular process）以及代謝進(jìn)程(metabolic process)差異表達(dá)基因數(shù)量居多；在分子功能中差異表達(dá)基因集中于催化活性(catalytic activity)和結(jié)合元件（binding）兩個方面，另外，轉(zhuǎn)運活性（transporter activity）以及核酸結(jié)合轉(zhuǎn)錄因子活性(nucleic acid binding transcription factoractivity)也出現(xiàn)了基因差異表達(dá)（表2）。表明白蠟樹韌皮部在白蠟窄吉丁脅迫下主要通過以上途徑的差異表達(dá)基因來進(jìn)行相關(guān)應(yīng)答及調(diào)控。

圖2 差異表達(dá)基因GO 富集分析Fig.2 Enrichment analysis of differentially expressed Gene Ontology

表2 差異基因GO 功能注釋分類統(tǒng)計Table 2 Statistical table of functional annotation classification of differential gene GO

2.4 健康與受害的韌皮部KEGG 代謝途徑分析

為了系統(tǒng)分析健康與受害白蠟樹韌皮部差異表達(dá)基因的基因功能、聯(lián)系基因組信息和功能信息，作者將差異表達(dá)基因按照參與的pathway 通路或行使的功能進(jìn)行分類，共建立了20 條代謝通路(pathway)。圖3 表示健康與受害白蠟樹韌皮部差異基因KEGG 代謝途徑，由圖可知，這些差異基因參與六大類代謝，分別為代謝（Metabolism）、遺傳信息處理（Genetic Information Processing）、環(huán)境信息處理（Environmental Information Processing）、細(xì)胞過程（Cellular Processes）、生物體系統(tǒng)（Organismal Systems）和人類疾?。℉uman Diseases）。這些差異基因主要參與代謝通路較多，其中碳水化合物代謝（Carbohydrate metabolism）的差異表達(dá)基因最多，達(dá)到6 397 個，其次是氨基酸代謝 (Amino acid metabolism)，達(dá)到3 632 個，另外這些差異基因參與的其他代謝通路有能量代謝、脂質(zhì)代謝以及其他次生代謝物的生物合成（表3）。結(jié)果表明，這些基因在多個方面較多的參與受害白蠟樹響應(yīng)白蠟窄吉丁脅迫的生理過程，而次級代謝產(chǎn)物極有可能在此過程中發(fā)揮重要作用。

2.5 健康與受害的韌皮部差異基因KEGG 功能富集分析

為獲取健康與受害白蠟樹韌皮部差異表達(dá)基因的主要參與代謝通路，作者對集中的差異基因進(jìn)行KEGG 功能富集分析，發(fā)現(xiàn)健康與受害白蠟樹韌皮部差異基因主要在122 條通路(pathway)中顯著富集。圖4 表示q-value 值最?。?.0～0.2）的前20 條pathway，差異表達(dá)基因數(shù)量以點的大小來表示，不同的q-value 范圍以點的顏色區(qū)分。由圖可知，內(nèi)質(zhì)網(wǎng)蛋白加工（Protein processing in endoplasmic reticulum）的差異表達(dá)基因最多，達(dá)到78 個；其次是內(nèi)吞作用 (Endocytosis)，達(dá)到75個（表4）。另一方面，植物激素信號轉(zhuǎn)導(dǎo)（Plant hormone signal transduction）、植物-病原互作（Plantpathogen interaction）、MAPK信號通路-植物（MAPK signaling pathway-plant）、RNA 運輸（RNA transport）、氨基糖和核苷酸糖代謝（Amino sugar and nucleotide sugar metabolism）、淀粉和蔗糖代謝（Starch and sucrose metabolism）、嘌呤代謝（Purine metabolism）的q-value 值最小，說明這些通路富集程度最高。結(jié)果表明，蟲害脅迫時以上通路中相關(guān)基因的差異表達(dá)最為活躍。

圖3 差異基因KEGG 代謝途徑分類統(tǒng)計柱狀圖Fig.3 Histogram of metabolic pathway classification of differential gene KEGG

表3 差異基因KEGG 代謝途徑分類統(tǒng)計Table 3 Classification and statistics of metabolic pathway of differential gene KEGG

2.6 韌皮部響應(yīng)害蟲侵害后轉(zhuǎn)錄因子家族分析

通過對健康與受害白蠟樹韌皮部進(jìn)行轉(zhuǎn)錄因子家族分析，共發(fā)現(xiàn)了20 個轉(zhuǎn)錄因子家族，分別為MYB 超級家族、bZIP、C2H2、C2C2、AP2/ERF、bHLH、C3H、NAC、WRKY、LBD、B3 超級家族、GRAS、HSF、FAR1、MADS、LOB、TCP、SBP、NF-Y 和CAMTA 轉(zhuǎn)錄因子家族（圖5）。其中，C3H 轉(zhuǎn)錄因子家族差異顯著基因最多為73 個，上調(diào)48 個，下調(diào)25 個；其次為BHLH 轉(zhuǎn)錄因子家族，共58 個基因差異顯著，上調(diào)40 個，下調(diào)18個；另外，NAC、MYB、B3、GRAS、SBP 等轉(zhuǎn)錄因子家族基因表達(dá)量達(dá)到顯著差異（表5）。CAMTA 轉(zhuǎn)錄因子家族沒有差異顯著基因，推測此轉(zhuǎn)錄因子家族對白蠟窄吉丁脅迫無明顯響應(yīng)作用。綜上所述，以上20 個轉(zhuǎn)錄因子家族中除CAMTA轉(zhuǎn)錄因子家族外均可能參與白蠟樹響應(yīng)白蠟窄吉丁脅迫的過程，其中C3H、BHLH、NAC、MYB 和B3 轉(zhuǎn)錄因子家族極有可能在白蠟樹抵御白蠟窄吉丁受害的過程中發(fā)揮重要作用，該結(jié)果為本試驗后續(xù)研究奠定理論基礎(chǔ)。

圖4 差異基因KEGG 富集分析Fig.4 Enrichment analysis of differential gene KEGG

表4 絨毛白蠟受白蠟窄吉丁危害下差異基因富集程度排名前20 的pathway 條目Table 4 The top 20 pathway entries of differential gene enrichment under pest stress of ash tree

圖5 健康與受害白蠟樹韌皮部轉(zhuǎn)錄因子家族統(tǒng)計Fig.5 Family statistics of transcription factors in phloem of Fraxinus velutina

表5 轉(zhuǎn)錄因子家族基因表達(dá)量差異統(tǒng)計Table 5 Statistical table of gene expression difference of transcription factor family

3 討論

白蠟樹是我國重要的木材與觀賞性植物，具有重要的經(jīng)濟(jì)及生態(tài)價值[8]，是我國城市園林綠化的重要樹種。白蠟樹易受到蛀干害蟲和食葉害蟲的危害，尤其蛀干害蟲對白蠟樹造成的危害極其顯著，其中以白蠟窄吉丁危害較為突出[9-10]，其幼蟲生活在白蠟樹韌皮部并以其為食，形成S 型蟲道，從而切斷營養(yǎng)及水分的運輸，最終導(dǎo)致白蠟樹的死亡[11]。目前國內(nèi)外對于白蠟窄吉丁習(xí)性及其入侵機(jī)制的研究較為深入，但生產(chǎn)上對白蠟窄吉丁的防治主要依靠化學(xué)農(nóng)藥，該措施雖然可以降低蟲口密度，但同時造成農(nóng)藥殘留和環(huán)境污染等問題。因此，如何采用綠色防控手段以高效、環(huán)保、可持續(xù)的方式控制白蠟窄吉丁為害，是白蠟樹產(chǎn)業(yè)亟待解決的問題，其中利用寄主樹本身的抗性來抵御害蟲的為害是當(dāng)前害蟲防控研究的熱點之一。本研究以健康與受害絨毛白蠟樹韌皮部為試驗材料，利用Illumina 測序技術(shù)進(jìn)行轉(zhuǎn)錄組測序，但由于白蠟樹屬于木犀科，目前還尚未有相關(guān)同屬或同科林木基因組或轉(zhuǎn)錄組的報道，故而采用無參轉(zhuǎn)錄組分析。研究結(jié)果為后續(xù)開展白蠟樹與白蠟窄吉丁互作的分子機(jī)制及挖掘關(guān)鍵抗蟲基因提供科學(xué)依據(jù)。

植物抵御生物脅迫是一個極其復(fù)雜的過程，涉及到許多相關(guān)基因的調(diào)控[12]。本研究從健康與受害韌皮部的轉(zhuǎn)錄組數(shù)據(jù)中篩選到應(yīng)答脅迫的關(guān)鍵差異轉(zhuǎn)錄因子家族，同時從DEGs 的GO 功能富集方面來看，生物學(xué)過程中以細(xì)胞過程和代謝過程差異表達(dá)基因居多，且已達(dá)到顯著水平；分子功能中以催化活性和結(jié)合元件差異表達(dá)基因居多，證明其在絨毛白蠟樹蟲害脅迫生理調(diào)控中有重要作用。這與Nalam[13]研究的9-脂氧合酶（LOXs）的活性影響擬南芥與病原菌和昆蟲相互作用的結(jié)果一致。KEGG代謝途徑分析結(jié)果顯示代謝途徑差異表達(dá)基因最多，而Melvin[14]等人研究植物可通過Hsp 蛋白誘導(dǎo)甲基乙二醛來調(diào)節(jié)糖酵解及其他代謝途徑從而達(dá)到抵御生物和非生物脅迫的作用，因此作者推測白蠟樹主要在代謝途徑抵御蟲害脅迫。另外，Misra[15]等人在煙草中表達(dá)了擬南芥的轉(zhuǎn)錄因子AtMYB12，使轉(zhuǎn)基因煙草增加了蘆丁的積累而對害蟲產(chǎn)生抗性，基于此，本研究通過對健康與受害白蠟樹韌皮部轉(zhuǎn)錄因子家族差異表達(dá)基因的統(tǒng)計挖掘，初步鎖定并推測C3H、BHLH、NAC、MYB、B3 和GRAS等轉(zhuǎn)錄因子家族參與了絨毛白蠟樹響應(yīng)白蠟窄吉丁脅迫的生理過程。綜上，這些代謝通路可能參與了絨毛白蠟韌皮部響應(yīng)白蠟窄吉丁危害的過程，通路上的一些關(guān)鍵基因或轉(zhuǎn)錄因子基因的表達(dá)可能是植物抵御蟲害脅迫的關(guān)鍵。由此，作者推測一些差異表達(dá)的基因是影響絨毛白蠟樹抵御蟲害脅迫的關(guān)鍵因子或候選基因，這將是本課題組下一步將要開展的研究內(nèi)容。

基于以上結(jié)果，本實驗室已初步獲得了白蠟樹抵御白蠟窄吉丁脅迫的分子證據(jù)，這將為我們從分子水平進(jìn)一步挖掘影響白蠟樹抗蟲的關(guān)鍵基因，并對其進(jìn)行生物學(xué)分析及功能驗證，對關(guān)鍵基因作用的分子機(jī)制及相關(guān)代謝通路進(jìn)行研究，進(jìn)而從形態(tài)學(xué)、轉(zhuǎn)錄組學(xué)和分子生物學(xué)水平解析白蠟樹抵御白蠟窄吉丁脅迫的應(yīng)答策略，為白蠟樹種與窄吉丁昆蟲互作的分子機(jī)制研究提供科學(xué)依據(jù)。

4 結(jié)論

本研究通過對健康與受害絨毛白蠟樹韌皮部的轉(zhuǎn)錄組分析，共鑒定出DEGs 3 388 個，其中上調(diào)2 141個，下調(diào)1 247 個。將差異基因劃分為20 個GO 功能類別，參與20 個KEGG 代謝途徑，另外差異表達(dá)基因分別在122 條通路中均有富集，包括植物-病原體互作等。健康與受害絨毛白蠟樹韌皮部共有20 個轉(zhuǎn)錄因子家族，表達(dá)量均達(dá)到顯著差異。研究結(jié)果為揭示白蠟樹應(yīng)對蟲害脅迫反應(yīng)的分子機(jī)制提供分子與理論依據(jù)。