绒毛白蜡（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 个转录因子家族，表达量均达到显著差异。研究结果为揭示白蜡树应对虫害胁迫反应的分子机制提供分子与理论依据。