絨毛白蠟（Fraxinus velutina Torr）韌皮部響應白蠟窄吉丁（Agrilus planipennis Fairmaire）危害的轉錄組變化

封潤霞，趙 婕，張蘇芳，王建軍，魏建榮，劉建鳳＊

(1.河北大學生命科學學院，河北 保定 071000；2.中國林業科學研究院森林生態環境與保護研究所，北京 100091；3.遼寧省林業科學研究院，遼寧沈陽 110032)

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

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

1 材料與方法

1.1 試驗材料

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

1.2 RNA 抽提

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

1.3 文庫構建及RNA-Seq 測序

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

1.4 數據的處理與篩選

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

2 結果與分析

2.1 幼蟲受害后白蠟韌皮部轉錄組數據質量分析

分別取健康（Non-infested 1-3）和受害（ASF 1-3）的白蠟樹韌皮部，對樣品文庫進行了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)以及低質量(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測序結果比例分別為98.11%、98.34%、98.3%、98.29%、98.3%和98.25%；堿基G 和C 的數量總和占總的堿基數量的百分比分別為44.75%、43.58%、43.3%、43.57%、43.4%和43.37% (表1)。結果表明測序質量合格，可進行下一步的生物學分析。

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

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

表1 轉錄組優化組裝結果評估Table 1 Evaluation of transcriptome optimized assembly results

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

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

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

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

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

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

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

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

為獲取健康與受害白蠟樹韌皮部差異表達基因的主要參與代謝通路，作者對集中的差異基因進行KEGG 功能富集分析，發現健康與受害白蠟樹韌皮部差異基因主要在122 條通路(pathway)中顯著富集。圖4 表示q-value 值最小（0.0～0.2）的前20 條pathway，差異表達基因數量以點的大小來表示，不同的q-value 范圍以點的顏色區分。由圖可知，內質網蛋白加工（Protein processing in endoplasmic reticulum）的差異表達基因最多，達到78 個；其次是內吞作用 (Endocytosis)，達到75個（表4）。另一方面，植物激素信號轉導（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 值最小，說明這些通路富集程度最高。結果表明，蟲害脅迫時以上通路中相關基因的差異表達最為活躍。

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

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

2.6 韌皮部響應害蟲侵害后轉錄因子家族分析

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

圖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 健康與受害白蠟樹韌皮部轉錄因子家族統計Fig.5 Family statistics of transcription factors in phloem of Fraxinus velutina

表5 轉錄因子家族基因表達量差異統計Table 5 Statistical table of gene expression difference of transcription factor family

3 討論

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

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

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

4 結論

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