李静宇，许林兵，应帆，肖维强，陈谷

（1.华南理工大学食品科学与工程学院，广东广州 510640）（2.广东省农业科学院果树研究所，广东广州 510640）

香蕉枯萎病（也称为巴拿马病）是由真菌尖孢镰刀菌古巴专化型（Fusarium oxysporumF. sp.cubense，Foc）引起的土传性真菌流行病害，可造成香蕉枯萎黄化甚至死亡，严重威胁香蕉产业的发展[1,2]。它在20世纪中叶之前摧毁了大蜜哈类（Gros Michel，AAA）香蕉的出口贸易产业，并威胁到用来取代它的香蕉卡文迪许（Cavendish，AAA）品系的生存。卡文迪许品系，因果形大，产量高，货架期长，在产蕉区被广泛种植，卡文迪许品系产蕉量目前占所有香蕉产量的约45%。目前，卡文迪许品系仍是产蕉区主要种植品种。然而，香蕉枯萎病菌4号生理小种（Foc4）是侵染能力和毒性最强的尖孢镰刀菌生理小种，可以侵染所有卡文迪许品系[2]。因而详细研究Foc4具有重要科学意义和应用价值。

Micro RNA（miRNAs）是一类内源性的具有调控功能的非编码RNA，其大小长约20～24个核苷酸[3]。MiRNA在动物，植物，病毒和藻类中广泛存在，在植物中成熟体miRNA能通过切割互补配对的靶基因mRNA抑制其表达，而动植物中都存在成熟体miRNA结合到与其互补的mRNA的位点抑制靶基因的翻译[4]。真菌中miRNA产生机制与动植物中miRNA产生机制不同，动植物中miRNA成熟体依赖Dicer酶剪切作用产生，而在真菌中除了有依赖Dicer酶作用产生miRNA外，还发现有和Argonaute蛋白QDE-2、核酸外切酶QIP和包含RNase III结构域的蛋白MRPL3结合产生miRNA的机制，所以在真菌中它们被称为miRNA-like（milRNA）[5,6]。从2010年在链孢霉（Neurospora）中真菌milRNA首次被发现以来，在西瓜枯萎病菌（Fusarium oxysporumF. sp.niveum），红色毛藓菌（Trichophyton rubrum）、小麦锈菌（Puccinia triticina）、金龟子绿僵菌（Metarhizium anisopliae）、葡萄座腔菌（Botryosphaeria dothidea）、新月弯孢霉（Curvularia lunata）、核 盘 菌（Sclerotinia sclerotiorum）、马 尔 尼 菲 青 霉 菌（Penicillium marneffei）、番茄枯萎病菌（Fusarium oxysporumF. sp.lycopersici）、里氏木霉（Trichoderma reesei）和灰盖鬼伞菌（Coprinopsis cinerea）中相继报道了milRNA的鉴定和分析[4,5,7-16]。大部分真菌milRNA鉴定均先通过与miRbase中植物或动物前体和成熟体比对鉴定保守milRNA，再通过发夹结构原则鉴定新型milRNA；进而根据利用软件预测milRNA的靶标基因，预测其潜在功能。并发现真菌milRNA对真菌生长发育及致病机制具有重要调控作用。

迄今为止，国内外尚未见期刊文献报道香蕉枯萎病致病菌尖孢镰刀菌（Foc）的小分子RNA（milRNA），仅在孟春亮的硕士论文中提及对Foc4进行的小分子RNA测序鉴定[17]，但未对Foc4进行降解组测序寻找milRNA的靶基因并深入研究。相对于植物和动物，真菌滞后的miRNA研究阻碍了Foc4 milRNA的研究和应用。因而，本研究通过小分子RNA高通量测序鉴定Foc4菌丝中milRNA，并通过降解组测序预测milRNA的靶基因，阐述Foc4中milRNA可能的调控功能，为更好地研究Foc4和香蕉枯萎病奠定基础。

1 材料与方法

1.1 实验材料

本实验所用香蕉枯萎病菌4号生理小种购于广东省微生物研究所菌种保藏与应用重点实验室。购回后在本实验室进行活化培养，使用18S测序鉴定其种属。挑取Foc4菌丝于PDA培养基中28 ℃中恒温培养9 d，收集平板上真菌菌丝保存于负80 ℃冰箱中，用于后续测序。

1.2 试验方法

1.2.1 小RNA文库构建及测序

将实验样本Foc4菌丝送至杭州联川生物有限公司进行小RNA及降解组测序。采用天根DP441试剂盒提取其总RNA，然后采用TruSeq Small RNA Sample Prep Kits（Illumina，San Diego，USA）试剂盒按照标准流程进行小RNA文库制备。文库制备工作完成后，对构建好的文库使用Illumina Hiseq2000/2500进行测序，测序读长为单端1×50 bp。

1.2.2 小RNA测序数据处理

首先对下机的原始数据进行测序质量评估，评估之后再对数据进行处理。本实验使用联川生物miRNA数据分析软件ACGT101-miR（LC Sciences，Houston，Texas，USA）进行分析，该软件的分析流程如下：（1）去除3’接头和垃圾序列：获取clean data；（2）长度筛选：筛选出clean data中碱基长度在18～25 nt的数据；（3）各种RNA数据库比对分析：将18～25 nt序列与不包含miRNA的mRNA、RFam[18]和Repbase数据库进行比对，若比对上这些数据库，则过滤这些非miRNA数据，剩余序列为有效数据。

1.2.3 miRNA鉴定

将有效数据与miRbase[19]数据库中植物miRNA和文献中报道的真菌[4,7,9,10,12]的前体及成熟体序列比对，详见附表3，其中使用的Foc4参考基因组为第三代测序组装所得[20]，并且根据其比对情况对其进行分组，gp1a是测序物种本身已报道的，known miRNAs；gp1b，2a，2b，3是保守的conservative miRNAs，gp2有基因组位置支撑，gp3没有基因组位置支撑，gp2a有前体基因支撑，gp2b没有前体基因支撑；gp4是新的novel miRNAs；具体各组分类标准如下所示：

（1）gp1a为可以比对到miRBase已知的首选物种前体（pre-miRNAs），并且前体可以进一步比对到基因组。

（2）gp1b Reads可以比对到miRBase已知的选择物种的前体（pre-miRNAs），并且前体可以进一步比对基因因组。

（3）gp2a Reads可以比对到miRBase选择物种的前体（pre-miRNAs），前体不能进一步比对到该物种的基因组，但是reads可以比对到基因组上。延伸的基因组序列可以形成满足11条原则发夹结构。

（4）gp2b Reads可以比对到miRBase选择物种的前体（pre-miRNAs），前体不能进一步比对到该物种的基因组，但是reads可以比对到基因组上。延伸的基因组序列不可以形成满足11条原则发夹结构。

（5）gp3a Reads可以比对到miRBase选择物种的前体（pre-miRNAs），前体不能进一步比对到基因组，reads也不可以比对到基因组上。cluster中比对上的Reads数目拷贝数>1，代表性的Reads分值>260，没有错配。

（6）gp4a Reads不能比对到miRBase选择物种的前体（pre-miRNAs），reads可以比对到基因组上。延伸的基因组序列可以形成满足11条发夹结构原则。

1.2.4 降解组文库构建及测序

降解组测序的原理：在植物体内绝大多数的miRNA是利用剪切作用调控靶基因的表达，且剪切常发生在miRNA与mRNA互补区域的第十位核苷酸上。靶基因经剪切产生两个片段，5’剪切片段和3’剪切片段。其中3’剪切片段，包含有自由的5’单磷酸和3’polyA尾巴，可被RNA连接酶，连接产物可用于下游高通量测序；而含有5’帽子结构的完整基因，含有帽子结构的5’剪切片段或是其他缺少5’单磷酸基团的RNA是无法被RNA酶连接，因而无法进入下游的测序实验；对测序数据进行深入地比对分析，可以直观地发现在mRNA序列的某个位点会出现一个波峰，而该处正是候选的miRNA剪切位点。

降解组建库包含以下步骤：1.通过磁珠捕获mRNA，3.5’adaptor连接；2.Biotinylated Random Primers和mRNA的混合反转录；3.PCR扩增，完成整个文库制备工作后，构建好的文库用Illumina Hiseq2000/2500进行测序，测序读长为单端1×50 bp。

1.2.5 降解组测序数据处理及分析将下机后的降解组数据进行如下分析：

（1）测序获得的原始数据通过一系列数据处理得到可用于后续分析的可比对测序序列。

（2）将可比对序列与测序物种的cDNA数据库序列比对生成降解组密度文件（degradome density file）。

（3）通过剪切位点预测软件（GSTAr）预测出与测序物种小RNA序列配对的靶基因mRNA序列。

（4）将预测的miRNA对应的靶基因和生成降解组密度文件中的mRNA进行结合运算，找出共同具有的mRNA，该mRNA即为miRNA的靶基因。并给出降解组的峰值分类和分值，并对产生的预测结果进行作图（t-plots）。降解组峰值分类Category 0表示原始数据片段在该位置，丰度等于该转录RNA上的丰度最大值，并且只有1个最大值；Category 1表示原始数据片段在该位置，丰度等于该转录RNA上的丰度最大值，并且不只1个丰度最大值；Category 2表示原始数据片段在该位置，丰度小于该转录RNA上的丰度最大值，但却大于该转录RNA上的丰度中间值；Category 3表示原始数据片段在该位置，丰度小于或等于该转录RNA上的丰度中间值；Category 4表示原始数据片段在该位置，只有1个与该转录RNA比对上。降解组等级代表milRNA剪切靶基因得到片段的丰度情况，降解组等级越低代表预测的靶基因越可靠。

1.2.6 milRNA功能分析

对通过降解组测序获得的milRNA靶基因进行GO和KEGG功能富集分析，并结合KEGG数据库中的代谢通路描述、基因功能预测和基因注释信息，查阅相关文献中基因功能，对milRNA靶基因功能进行解析。

2 结果与讨论

2.1 小RNA测序数据的统计与分析

如表1所示，小RNA测序共获得11783990条原始数据（raw reads），原始数据去除重复数据后得到unique reads为2326111条。首先去除原始数据中接头和垃圾序列等低质量片段，然后对清洁读数进行长度筛选，保留18～25 nt长度的序列。图1为18～25 nt长度序列的统计分布。将18～25 nt序列与不包含miRNA的mRNA、RFam和Repbase数据库进行比对，若比对上这些数据库，则过滤这些非miRNA数据，剩余序列为有效数据（valid reads）。比对数据库后，获得有效数据条数3351578条，占原始数据的比例为28.44%。有效数据去除重复后得到unique reads为744349条，用于鉴定Foc4中milRNA。

表1 小RNA文库组成Table 1 Composition of the small RNA library

图1 小RNA长度分布Fig.1 Size distribution of small RNAs

2.2 milRNA的鉴定分析

将测序得到的valid reads与miRbase中植物和文献中收集得到的真菌miRNA成熟体及前体序列（附表3）比对，如表2所示，鉴定获得Foc4中7个保守的（conservative）milRNA，包括4个与真菌比对上的milRNA（Foc-milR1，Foc-milR2 Foc-milR3和Foc-milR4）和3个与植物比对上的milRNA（Foc4-milR5，Foc4-milR6和Foc4-milR7），还有3个新型的（novel）milRNA（Foc4-novel milR1，Foc4-novel milR2和Foc4-novel milR3），详细信息见表1。其中4个与真菌比对上的milRNA分组为gp1b，即比对到miRBase已知的选择物种的前体（pre-miRNAs），并且前体可以进一步比对到基因组。例如，Foc4中预测的Foc4-milR4与番茄枯萎病菌（Fon）中fon-miR-1-m0017比对上，R-1代表Foc4-milR4右端少掉1个碱基，且其前体可以进一步比对到Foc4基因组。三个与植物比对上的milRNA分组为gp2a，表示可与植物前体比对上，但该植物前体并不能比对上Foc4基因组，该成熟体milRNA可以比对到Foc4基因组上，且在成熟体milRNA周围延伸的基因组序列可以形成满足11条发夹结构原则的前体。Foc4-milR5，Foc4-milR6和Foc4-milR7分别与多个miRbase已知miRNA比对上被注释为：mtr-MIR5229a-p3_2ss16AG18TG，ppe-MIR399k-p5_2ss13GC20GA和zma-MIR162-p3_2ss4GA19AT。Foc4-milR7的注释为zma-MIR162-p3_2ss4GA19AT表示可与miRbase中玉米（Zea mays）前体MIR162比对上且位于前体的3’端，且2ss4GA19AT表示在第4个碱基G由A替换，在第19个碱基A由T替换，共计2个替换，且其分组为gp2b，表示玉米前体MIR162不能比对上Foc4基因组。可以发现，三个与植物比对上的Foc4-milRNA成熟体中都存在两个碱基的替换，且比对上的前体都不能比对上Foc4基因组，表明miRNA虽然在植物和真菌间有保守的成熟体，但其前体在植物和真菌间变化差异较大。

3个novel milRNA其分组为gp4a，表示该序列不能比对上miRbase中已选择物种的前体，但可以比对上Foc4基因组上，且延伸基因组序列可以形成满足11条发夹结构原则的前体。

图2 milRNAs前体发夹结构Fig.2 The hairpin structures of milRNAs precursors

在十个鉴定milRNA中，Foc4-milR3和Foc4-novel milR3的表达量较高，表明香蕉枯萎病菌Foc4的milRNA可能具有重要作用。图2所示为gp2a分组和gp4a分组milRNA的前体发夹结构。所鉴定的milRNA中最短的Foc4-milR7的前体长度为56 nt，而最长的Foc4-novel milR1的前体长度为233 nt，可以发现真菌milRNA的前体长度变化非常大，有研究也发现在链孢霉和灰盖鬼伞菌中milRNA的前体长度也非常的灵活[5,16]。

表2 Known and predicted milRNATable 2 Summary of known and predicted milRNA in Foc4

2.3 milRNA靶基因的分析鉴定

进而对Foc4菌丝进行降解组测序，产出原始数据27984500条，去除重复后得到unique raw reads为6695437条，除去3’接头后剪切位点标签较短（<15 nt）的序列后，能比对上Foc4转录本数据库的序列（Mapped Reads）为19773216条，占raw reads的比率为70.66%，unique Mapped Reads为4848061条，占Unique Raw Reads的比率为72.44%（表3），依此生成降解组密度文件。将Targetfinder软件根据碱基互补配对预测的靶基因和降解组测序所收集到的剪切片段结合运算，找出共同拥有的mRNA，即为Foc4 milRNA的靶基因。

表3 降解组数据总览Table 3 Data summary of degradome sequencing

降解组测序总共预测出53对milRNA-mRNA，有7个milRNA通过降解组测序发现其潜在靶基因，分别 为Foc4-milR2，Foc4-milR4，Foc4-milR5，Foc4-milR6，Foc4-milR7，Foc4-novel milR1和Foc4-novel milR3。其中，Foc4-novel milR1所调控靶基因数目最多为27个，其靶基因注释包括染色质结构重塑复合亚基snf21（Chromatin structure-remodeling complex subunit snf21），腺苷酸脱氢酶（AMP deaminase），细胞壁合成蛋白psu1（Cell wall synthesis protein psu1），线粒体ATP合成酶亚基α（ATP synthase subunit alpha，mitochondrial），苹果酸脱氢酶（malate dehydrogenase），ABC转运器CDR4（ABC transporter CDR4），假定螺旋酶mug81（Putative helicase mug81），线粒体细胞色素b-c1复合亚基Rieske（cytochrome b-c1 complex subunit Rieske，mitochondrial），过氧化物生物发生因子19（Peroxisomal biogenesis factor 19）和孢子壁成熟蛋白DIT1（Spore wall maturation protein DIT1）等。此外，Foc4-milR5调控14个靶基因，其中靶基因编码F型转运氢离子的ATPase亚基C（F-type H+-transporting ATPase subunit C），碱性神经酰胺酶（alkaline ceramidase），羟异羟乙酸水解酶（hydroxyisourate hydrolase）和CAMKK蛋白激酶（CAMKK protein kinase）等。Foc4-novel milR3是表达水平最高的milRNA，其靶基因编码乌头水合酶（Aconitate hydratase），3-植酸酶A（3-phytase A）和真核翻译起始因子3亚基B（eukaryotic translation initiation factor 3 subunit B）等。其他milRNA靶基因详细信息见表2。图3列出了几个milRNA与其靶标基因的T-plot图，例如图3a中，Foc4-milR2靶向基因Foc4_580068650，降解组等级（category）为2，milRNA剪切位点在靶基因1880处（T-plot图中红色点标记处）。

图3 T-plot图Fig.3 T-plot images

从表2中发现，Targetfinder软件通过碱基互补配对对milRNA靶基因进行评分（AllenScore），milRNA所有靶基因均不是完全匹配的，且匹配率最高的milRNA-mRNA的AllenScore为3.5分。在链孢霉中绝大多数milRNA也没有发现完美匹配的靶基因，推断真菌milRNA可能靶向不完全互补的序列[5]。Foc4中milRNA也可能是通过靶向不完全互补的基因发挥调控作用。

2.4 靶基因功能分析

图4 milRNA靶基因GO分类Fig.4 Gene ontology categories of the target genes of Foc4 milRNAs

对通过降解组测序预测的milRNA的靶基因进行GO富集分析，从图4可以看出，在分子功能（molecular function）分类中，归类为分子功能的ATP结合（ATP binding）条目的基因较多。在细胞定位（cellular component）分类中，细胞核（nucleus）条目和细胞质基质（cytosol）条目较多，在生物学过程（biological process）分 类 中，蛋 白 质 磷 酸 化（protein phosphorylation）和致病机制（pathogenesis）条目中基因数较多。表明Foc4中milRNA可能调控多个生物学过程，具有重要意义。

图5 milRNAs靶基因KEGG分类Fig.5 KEGG pathways categories of the target genes of Foc4 milRNAs

MilRNA的靶基因能归入KEGG通路中的基因有20个。如图5所示，靶基因归入氧化磷酸化通路（Oxidative phosphorylation）和嘌呤代谢通路（Purine metabolism）的 基 因 数 最 多，脂 代 谢 通 路（Glycerophospholipid metabolism），丙酮酸盐代谢（Pyruvate metabolism）和核糖体通路（Ribosome）次之。从KEGG通路注释分类图，可以发现Foc4 milRNA的靶基因参与多个代谢过程。其中表达量较高的milRNA Foc4-milR4的靶基因Foc4_580004970编码假定蛋白参与调控甘油磷脂代谢通路（Glycerophospholipid metabolism），Foc4-novel milR3靶向的Foc4_580043960基因编码3-植酸酶A（3-phytase A）可参与调节硫胺素代谢（Thiamine metabolism）和核黄素代谢（Riboflavin metabolism）。表明Foc4 milRNA可能对代谢过程有重要影响。香蕉枯萎病菌能量代谢过程对其生长发育和致病性有重要作用，本实验中发现Foc4-novel milR3的靶基因通过编码乌头水合酶（Aconitate hydratase）参与调控三羧酸循环通路，Foc4-milR5靶基因为Foc4_580004190编码F型转运H+的ATPase亚基C和Foc4-novel milR1靶基因Foc4_580041730编码线粒体ATP合成酶亚基α，两者都是调控与ATP合成相关的酶参与氧化磷酸化通路。在葡萄座腔菌响应病毒胁迫时也发现差异表达milRNA调控代谢过程[10]。

真菌milRNA不仅可以调控代谢过程影响其生长发育，也会影响其致病性。例如，在西瓜枯萎病菌中，发现Fon-miR7696a-3p和Fon-miR6108a分别靶向单端孢霉毒素（trichothecene）和乙烯诱导肽1（NEP1）调节毒素合成[4]；在金龟子绿僵菌和新月弯孢酶中发现milRNA对菌丝体生长和孢子的形成有重要调控作用，可能影响致病性[9,11]。在本研究中也发现了一些milRNA调控致病性相关基因表达，例如，Foc4-novel milR1靶向基因Foc4_580051390编码ABC转运器CDR4。ABC（ATP binding cassette）转运器由赋予特异性的跨膜结构域和含有高度保守的氨基酸基序的结构保守的核苷酸结合结构域组成，它们不仅充当转运蛋白，也是ATP水解产生的能量的受体或通道[21]。例如CDR4是链孢菌唑耐药的主要原因[22,23]，而在葡萄座腔菌中Bd-milRNA1147靶向ABC转运器可能会抑制其侵染梨树[10,24]。在植物致病菌禾谷镰刀菌中敲除编码ABC转运器基因FgABC1和FgABCC9，发现其对小麦根系和小麦冠的毒力减弱[25,26]。推测Foc4-novel milR1可能通过调控ABC转运器CDR4影响其毒力。此外，Foc4-novel milR1靶基因Foc4_580122700编码孢子壁成熟蛋白DIT1（Spore wall maturation protein DIT1），可能与Foc4孢子形成相关。在新月弯孢酶中也发现milRNA参与孢子形成[11]。

Foc4在宿主香蕉内的生长发育和致病性都是影响香蕉枯萎病的重要因素，本研究将为香蕉枯萎病的防治提供新的思路和一定的参考。

3 结论

本研究利用小分子RNA测序对香蕉枯萎病致病菌Foc4进行测序，共鉴定7个保守的milRNA和3个新型的milRNA。利用降解组测序分析鉴定出53个受milRNA调控的潜在靶标基因，它们参与调控多个代谢通路，包括嘌呤代谢、甘油磷脂代谢、硫胺素代谢和三羧酸循环通路，对真菌生长发育有重要作用；同时鉴定到milRNA靶向ABC转运器CDR4和孢子壁成熟蛋白DIT1，可能对Foc4导致香蕉枯萎病发病具有重要影响。本研究首次在香蕉枯萎病致病菌Foc4中利用降解组测序预测milRNA靶基因，对研究真菌milRNA调控功能提供了重要参考，也为香蕉枯萎病的防治提供了新的思路。