鸭茅响应水淹胁迫的miRNA 差异表达分析

杨兴云，乔丹丹，张雅洁，王少青，任俊才，李明阳，屈明好，尚盼盼，杨成，黄琳凯，曾兵，*

（1. 西南大学动物科学技术学院，重庆402460；2. 草食动物科学重庆市重点实验室，重庆 400715；3. 福建傲农生物科技集团股份有限公司，福建 厦门 361111；4. 四川农业大学草业科技学院，四川 成都 611130）

鸭茅（Dactylis glomerata）又名果园草（orchard grass）、鸡脚草（cocksfoot），属禾本科（Poaceae）早熟禾亚科（Festucoideae）鸭茅属（Dactylis）［1］。鸭茅是一种高度异花授粉的冷季型多年生禾本科牧草，自然分布于欧洲及我国，其作为饲料作物有极高的营养价值以及较强的适应性和耐阴性［2−3］，但对水淹胁迫较敏感。在国内鸭茅主产区，由于洪涝灾害频繁，严重影响鸭茅的生长发育、品质和产量。但目前国内外关于鸭茅抗逆性研究主要集中在干旱胁迫［4］、耐热胁迫［5］、锈病胁迫［6］等方面，因此探究鸭茅在水淹胁迫下的抗逆分子机制有重要意义。

涝害作为植物正常生长发育的非生物胁迫因子，对植物产生的危害主要是阻碍植物根系的供氧，抑制根系呼吸，导致植物重要代谢过程的根细胞能量供应不足［7］，最终植物枯萎甚至坏死。水淹胁迫对植物生理生化的影响主要包括光合作用、根部抗氧化系统和根部激素合成三方面［8］。长期水淹后，光合作用产生的光合产物不能通过韧皮部向下运输，导致植物的代谢循环受阻，叶片净光合速率、光合色素、气孔导度降低［9−10］，从而抑制植物正常生长。此外，根系抗氧化系统也会发生相应变化，水淹胁迫后植株体内的抗坏血酸过氧化物酶（ascorbic acid peroxidase，APX）、谷胱甘肽还原酶（glutathione reductase，GR）、超氧化物歧化酶（superoxide dismutase，SOD）及过氧化物酶（peroxidase，POD）等含量急剧变化［11］。同时，水淹胁迫还影响植物体内激素合成与运输动态平衡，生 长 素（indole-3-acetic acid，IAA）、细 胞 分 裂 素（cytokinin，CTK）、脱 落 酸（abscisic acid，ABA）、赤 霉 素（gibberellin，GA）、乙烯（ethylene，ETH）等激素在不同部位增加。其中，乙烯作为衡量植物耐涝性重要的指标之一，在有氧部位积累，造成叶片老化脱落，抑制正常植物生长［12］。

小RNA（microRNAs，miRNA）是一类内源性的大小为21～22 个核苷酸的小分子单链RNA，位于基因组非编码区［13］。通过靶向mRNA 的切割或抑制翻译过程参与基因的表达调控［14］。miRNA 作为一种重要的调控因子，通过在转录后水平与其靶基因mRNA 的完美配对参与植物对环境胁迫的反应［15］。研究表明，miRNA 在植物抗逆胁迫反应中发挥着重要作用［16］，例如盐胁迫［17］、温度胁迫［18］、营养胁迫［19］和水分胁迫［20］等。植物遭遇水淹时，气体在水中的扩散速度降低，植物组织会出现严重的缺氧表现，已有相关研究揭示低氧环境下miRNA 如何参与相关基因表达调控［21−23］。Licausi 等［24］将拟南芥（Arabidopsis thaliana）根进行不同程度的缺氧处理，通过对1900个转录因子（transcription factors，TF）和180 个miRNA 差异表达分析后发现，低氧诱导基因的调控是通过各种TF 组合的同时相互作用来控制的，miRNA 在植物低氧基因调控的早期阶段起次要作用［24］。刘智捷［25］对3 种不同耐涝表型的玉米（Zea mays）幼苗进行短期渍水胁迫处理，分析各处理条件下根系miRNA 与其靶基因的表达动态，结果表明miR159、miR164、miR167、miR393、miR408和miR528基因主要参与根系发育和水淹胁迫响应。进一步分析发现，约有1/2 目的基因与其miRNA 表达量呈负相关，其中大部分是参与低氧胁迫反应的转录因子。在其他物种中，Jin 等［26］对水淹胁迫下的莲花（Nelumbo nucifera）进行高通量测序研究耐涝相关的miRNA，发现NNU_Far-miR159、NNU_GMA-miR393h 和NNU_Aly-miR319c-3p 在莲花对淹水的响应中起着重要的调节作用。

近年来，随着全球自然环境的持续发展恶化以及频繁的极端天气，导致某些地区暴雨频发，洪水灾害越来越严重［27−28］。中国也是遭受洪涝灾害较为严重的国家，特别是长江中下游地区，受灾面积占全国总受灾面积的3/4［29］。2020年8月，四川、重庆等地严重水灾造成了巨大社会经济损失。此外，持续的气候变暖预示着世界所有区域的洪涝事件的频率和强度还会增加［30］，这将导致更大的经济损失。鸭茅作为优质牧草，深度解析植物耐涝的分子响应机制对加快作物耐涝品种选育具有重要意义。虽有人对水淹胁迫后“斯巴达”鸭茅进行了转录组测序研究［31］，但对于鸭茅在水淹胁迫下miRNA 表达特征的研究尚未见报道。前期用人工模拟涝害胁迫，通过形态学、生理学根部显微结构以及组学研究［32−33］，筛选到了耐涝性较强的“滇北”鸭茅。

鉴于此，本研究以耐涝性较强的“滇北”鸭茅为试验材料，利用高通量测序与生物信息学技术，筛选与鸭茅水淹抗性相关的miRNA，分析它们在水淹胁迫中的差异表达，旨在为今后研究鸭茅抗涝作用机制及抗涝遗传改良工作提供参考，同时为耐涝草本植物新品种的选育提供借鉴。

1 材料与方法

1.1 试验材料与处理

以鸭茅品种“滇北”（简称“DB”）为试验材料，种子由四川农业大学草业科技学院提供。试验于2018年9月开始，筛选均匀饱满的鸭茅种子播种于育苗盘中，挑选长势大体一致的幼苗移栽至花盆（口径16 cm，高14 cm），再放入人工培养箱，温度设置为22 ℃/15 ℃（昼/夜），光周期14 h/10 h（昼/夜），光照10000 lx，土壤相对含水量保持在50%左右。根据植株生长需要每周浇2 次1/2 Hoagland 溶液，待植株长至5 cm 及叶片覆盖率超过80%时进行试验。

试验采用盆栽水淹胁迫，将3 盆长势良好且一致的鸭茅放入水箱（长80 cm×宽57 cm×高50 cm）中，然后加水淹没至植株叶片最高处，并始终保持水位淹没鸭茅顶部。另外3 盆不水淹，做对照处理。于水淹胁迫处理0 h取样对照组（CK），水淹胁迫处理8 和24 h 取样处理组（D），取样3 次，每次取样3 个重复，共9 个样品。采用混样法剪取鸭茅叶片，液氮迅速冷冻，后置于−80 ℃低温冰箱保存备用。样品送至武汉菲沙基因公司进行测序。

1.2 小RNA 文库的构建与高通量测序

利用Trizol 试剂盒（Invitrogen，USA）分别提取9 个鸭茅叶片样品中的总RNA。总RNA 质量检测：1）通过进行琼脂糖凝胶电泳分析RNA 降解和污染的程度；2）采用Nanodrop 测试RNA 浓度（OD260/280）；3）利用Qubit 对RNA 浓度进行精确定量；4）采用Agilent 2100 检测RNA 完整性。最后仅选择质量检测合格的总RNA 样品进行小RNA（small RNA，sRNA）文库构建及测序工作（Illumina Hiseq TM2100）。

1.3 测序原始数据预处理

对高通量测序获得的sRNA 进行生物信息学分析前，需先对原始数据进行杂质处理。数据过滤的步骤如下：1）去除低质量的reads，当测序reads 中含有的低质量（≤5）碱基数超过该条reads 长度比例的50%时，过滤掉该reads；2）去除含N（任一个碱基）比例较高的reads，当测序reads 中含有的N 的含量超过该条reads 长度比例的10%时，过滤掉该reads；3）去除有5′接头污染的reads；4）去除没有3′接头序列和插入片段的reads；5）去除3′接头序列；6）去除含polyA/T/G/C 的reads；7）去除最终长度小于17 bp 的reads；经过上面的步骤过滤数据之后，得到了大量的sRNA 序列。

1.4 响应水淹胁迫miRNA 的鉴定

将排除核糖体RNA（ribosomalRNA，rRNA）等的sRNA 序列在miRBase［34］数据库中进行比对分析，综合成熟体相似度和前体茎环结构两方面比对情况得到鸭茅中miRNA 的表达情况。具体策略条件如下：1）成熟体相似度方面考虑到鸭茅与其他物种的差异，比对中允许两个碱基错配；2）对条件1 中成熟体miRNA 进行pre-miRNA 的预测鉴定，若出现基因组上无pre-miRNA 位点或不能形成二级茎环结构的状况，则排除其为鸭茅保守miRNA，获得更为准确和高可行性的验证。

miRNA 前体的标志性发夹结构能够用来预测新的miRNA。基本原理是通过截取一定长度sRNA 比对上的参考序列，通过探寻其二级结构及Dicer 酶切位点信息、能量等特征进行分析，预测样品中新miRNA，并进行各样本中匹配上的sRNA 的序列、长度、出现的次数等信息，以及不同长度miRNA 的首位点碱基分布和所有miRNA的各位点碱基分布情况的统计。对于各个数据库中均未注释但可比对到基因组序列的sRNA，将作为候选的新miRNA。采用mirdeep2［35］软件进行新miRNA 的分析，综合考虑到sRNA 测序序列表达丰度，前体miRNA 最低自由能和复杂的二级结构，预测得到的候选新miRNA 前体的茎环结构图，完成新miRNA 的预测。

1.5 miRNA 差异表达分析

基因差异表达分析的输入数据为miRNA 表达水平分析中得到的miRNA 对应的read counts 数据。用R 语言包DESeq2 进行差异分析，筛选阈值为|log2 fold change|＞1&P＜0.05。对于之前得到的已知miRNA 和新miRNA，做miRNA 靶基因预测分析。

由于单个软件预测的miRNA 往往会得到特别多的靶基因，因此本研究选用TargetFinder 和psRobot 两种软件结果的交集作为最后的靶基因预测分析［36−37］，分别获取各组差异表达的miRNA，根据miRNA 与其靶基因的关系，对每组差异表达miRNA 的靶基因集合分别进行GO (gene ontology)功能分析和KEGG (kyoto encyclopedia of genes and genomes)富集分析。

2 结果与分析

2.1 鸭茅小RNA 测序数据及质量控制

应用高通量测序技术，对水淹胁迫0、8 和24 h 后的“滇北”鸭茅叶片进行小RNA 测序。测序结果由表1 可知，水淹胁迫后“滇北”叶片中获得109825656 原始测序数据（raw reads）；通过质量控制后，样品中未发现质量值较低序列，经过滤获得高质量序列（clean reads）数量是68789053；可从表中看出所有样本的GC 含量为51.2%～52.9%（GC 含量：碱基G 和C 的数量总和占总的碱基数量的百分比）；Q20≥99.2%，Q30≥98.1%［Q20、Q30：碱基质量值Qphred=−10 log10（error rate）数值大于20、30 的碱基占总体碱基的百分比］。这说明该测序结果可信且质量较高，可用于进一步的分析。

表1 不同鸭茅样品测序数据统计Table 1 Statistical results of sequencing data of different D.glomerata samples

2.2 鸭茅小RNA 分类注释

采用bowtie2［38］将小RNA 比对到mirBase（v21）数据库指定范围的miRNA 序列上，得到小RNA 与miRNA 比对情况。在水淹胁迫处理0、8 和24 h 的“滇北”鸭茅中鉴定到能比对上的miRNA 有3259、3309 和3188 个，分别占能比对上miRNA 的小RNA 总数（total sRNA）的4.22%、3.13%、2.59%。

根据与Rfam 数据库的比对作为非编码RNA（non-coding RNA，ncRNA）的注释结果，“滇北”鸭茅非编码RNA 序列可以注释为核仁小RNA（small nucleolar RNA，snoRNA）、核糖体RNA（ribosomalRNA，rRNA）、小核RNA（small nuclear RNA，snRNA）、转运RNA（transfer ribonucleic acid，tRNA）4 种。结合基因组上外 显子（exon），内含子（intron）和重复序列的信息，以及与Rfam，mirBase（v21）数据库的比对结果，比对上的小RNA 做分类统计。由表2 可知，比对到的miRNA 序列从11.00%上升到15.80%，再下降到14.00%，说明水淹胁迫处理8 h 时，鸭茅可能积极响应水淹胁迫。

表2 小RNA 分类统计Table 2 Small RNA classification statistics（%）

2.3 鸭茅保守miRNA 鉴定及其家族分析

植物中成熟的miRNA 序列高度保守，生物信息学方法用于比较基因组学研究从而探索不同物种中miRNA 基 因 家 族的 保 守 性［39］。基 于miRNA 比 对 分析结果，对miRNA 在鸭茅水淹胁迫叶片中的表达进行统计。通过表3 可知，推测8 h 是植物遭受水淹胁迫的代谢活跃期。Ren 等［40］对水淹胁迫处理8 h 后的毛白杨（Populus tomentosa）进行小RNA 测序分析，结果表明7 个保守的miRNA 家族和5 个新的miRNAs 在水淹胁迫下的表达发生了显著变化，对差异表达miRNAs 潜在靶点的注释表明，许多编码转录因子、酶和信号转导成分的基因参与了非生物胁迫反应。与水淹胁迫24 h 相比，8 h“滇北”鸭茅有72 个miRNA 发生差异表达，有38 个上调，34 个下调。由此可推断，水淹胁迫处理8 h 时，鸭茅通过调节生理代谢活动响应水淹胁迫。对检测到的已知miRNA 和新miRNA进行家族分析，探索其所属的miRNA 家族在其他物种中的存在情况。如果来源于同一个家族的miRNA，它们在2～7 这些位点（这几个位点是作用于靶基因的种子位点）表现为一致的保守型，就会被定义为一个miRNA 家族［41］，其中序列相似性比对，是最为简单直接的方法。对鉴定到的miRNA 进行DESeq2 差异分析，“滇北”鸭茅中总共鉴定到228 个差异表达的miRNA。可以发现鸭茅在水淹胁迫下发生差异表达的miRNA 主要属于miR166（78 个）、miR159（35 个）、miR167（29 个）、miR396（17个）、miR156（21 个）这5 个miRNA 家族（表4），表明这几个家族的miRNA 在植物响应水淹胁迫过程中发挥了重要作用。

表3 差异基因统计Table 3 Statistical of differentially expressed genes（DEGs）

表4 鸭茅中响应水淹胁迫的差异表达miRNA 家族Table 4 Differentially expressed miRNA families in D.glomerata in response to flooding stress

2.4 鸭茅miRNA 的靶基因预测

miRNA 主要通过与靶mRNA 的结合，或促使mRNA 降解，或阻碍其翻译，从而抑制目的基因的表达。用生物信息学方法准确快速地预测miRNA 的靶基因，可以为研究miRNA 功能提供线索。对不同水淹胁迫处理下的鸭茅材料中得到的已知miRNA 和新miRNA 做miRNA 靶基因预测分析，获得了miRNA 与靶基因的对应关系。利用TargetFinder 和psRobot，进一步预测了上述miRNA 的靶基因，共获得9136 个靶基因，5 个主要的miRNA 家族中分别鉴定到miR166（1817 个）、miR167（320 个）、miR159（835 个）、miR396（924 个）、miR156（5204 个）。再选取5 个基因家族中差异显著的miRNA 及其对应的靶基因进行分析，通过表5 可知，在研究中发现某些miRNA 有多个靶基因与之对应，表明miRNA 在响应水淹胁迫过程中并不是单基因调控，而是一个miRNA 可以参与调控多个生物代谢过程，这也证实了鸭茅在水淹胁迫中响应过程复杂。

表5 不同的靶基因对应同一个miRNATable 5 Different target genes correspond to the same miRNA

2.5 差异miRNA 靶基因GO 富集分析

GO（gene ontology）是基因功能国际标准分类体系。在预测水淹胁迫下响应miRNA 的靶基因后，研究这些靶基因在GO 中所映射到terms 的情况并对其进行生物学富集分析，可以对水淹响应miRNA 靶基因进行功能分类集合分析和具有生物学意义的基因功能描述。GO 为进一步解析miRNA 的功能和相关靶基因研究提供了思路，该数据库主要分为三大类：生物过程（biological process，BP）、细胞组成（cellular component，CC）、分子功能（molecular function，MF）［42−43］。

本研究对水淹胁迫下“滇北”鸭茅0 h vs 8 h 的70 个差异表达基因进行GO 功能富集分析，共注释到431 个GO terms，其中生物过程（biological process）分类富集到256 个GO terms，细胞组成分类富集到75 个GO terms，分子功能分类富集到105 个GO terms。并挑选了富集程度最靠前的31 个GO terms（图1），其中，显著富集的GO terms 有10 个（Q-value＜0.05），生物过程类别中显著富集的分别是：细胞过程（cellular process，GO：0009987）、代谢 过 程（metabolic process，GO：0008152）、生 物 调 控（biological regulation，GO：0065007）、生 物 过 程 调 控（regulation of biological process，GO：0050789）、发育过程（developmental process，GO：0032502）；细胞组成类别中显著富集的分别是：细胞（cell，GO：0005623）、细胞部分（cell part，GO：0044464）、细胞器（organelle，GO：0043226）；分子功能类别中显著富集的分别是：黏合物（binding，GO：0005488）、核酸结合转录因子活性（nucleic acid binding transcription factor activity，GO：0001071）。

图1 “滇北”淹水0 h vs 8 h 差异miRNA 靶基因GO 分类Fig.1 Histogram of GO classification of differential miRNA target genes for 0 h vs 8 h flooding in“Dianbei”1：细胞过程Cellular process；2：代谢过程Metabolic process；3：生物调控Metabolic process；4：生物过程调控Regulation of biological process；5：单一有机体过程Single-organism process；6：发育过程Developmental process；7：对刺激的反应Response to stimulus；8：定位Localization；9：生物过程的正向调节Positive regulation of biological process；10：信号Signaling；11：多细胞生物的过程Multicellular organismal process；12：细胞组成组织或生物发生Cellular component organization or biogenesis；13：生物过程的负调控Negative regulation of biological process；14：繁殖Reproduction；15：生殖过程Reproductive process；16：细胞Cell；17：细胞部分Cell part；18：细胞器Organelle；19：膜Membrane；20：膜部分Membrane part；21：细胞器部分Organelle part；22：大分子复合体Macromolecular complex；23：细胞外区域Extracellular region；24：超分子络合物Supramolecular complex；25：细胞连接Cell junction；26：合胞体Symplast；27：黏合物Binding；28：催化活性Catalytic activity；29：核酸结合转录因子活性Nucleic acid binding transcription factor activity；30：结构分子活性Structural molecule activity；31：转运活力Transporter activity.

对水淹胁迫下“滇北”鸭茅0 h vs 24 h 的66 个差异表达基因进行GO 功能富集分析，共注释到404 个GO terms 上，其中生物过程分类富集到242 个GO terms，细胞组成分类富集到73 个GO terms，分子功能分类富集到94 个GO terms。并挑选了富集程度最靠前的29 个GO terms（图2），其中，显著富集的GO terms 有11 个（Qvalue ＜0.05），其中生物过程类别中最显著富集的分别是：细胞过程（cellular process，GO：0009987）、代谢过程（metabolic process，GO：0008152）、生物调节（biological regulation，GO：0065007）、生物过程调控（regulation of biological process，GO：0050789）、发育过程（developmental process，GO：0032502）；细胞组成类别中最显著富集的分别是：细胞（cell，GO：0005623）、细胞部分（cell part，GO：0044464）、细胞器（organelle，GO：0043226）、超分子络合物（supramolecular complex，GO：0099080）；分子功能类别中最显著富集的分别是：黏合物（binding，GO：0005488）、核酸结合转录因子活性（nucleic acid binding transcription factor activity，GO：0001071）。

图2 “滇北”淹水0 h vs 24 h 差异miRNA 靶基因GO 分类Fig.2 Histogram of GO classification of differential miRNA target genes for 0 h vs 24 h flooding in“Dianbei”1：细胞过程Cellular process；2：代谢过程Metabolic process；3：生物调节Biological regulation；4：生物过程调控Regulation of biological process；5：单一有机体过程Single-organism process；6：发育过程Developmental process；7：对刺激的反应Response to stimulus；8：定位Localization；9：信号Signaling；10：生物过程的正向调节Positive regulation of biological process；11：多细胞生物的过程Multicellular organismal process；12：细胞组成组织或生物发生Cellular component organization or biogenesis；13：生物过程的负调控Negative regulation of biological process；14：繁殖Reproduction；15：生殖过程Reproductive process；16：细胞Cell；17：细胞部分Cell part；18：细胞器Organelle；19：膜Membrane；20：膜部分Membrane part；21：细胞器部分Organelle part；22：高分子复合物Macromolecular complex；23：超分子络合物Supramolecular complex；24：细胞外区域Extracellular region；25：黏合物Binding；26：催化活性Catalytic activity；27：核酸结合转录因子活性Nucleic acid binding transcription factor activity；28：结构分子活性Structural molecule activity；29：转运活力Transporter activity.

综上可知，“滇北”鸭茅水淹0 h vs 8 h 和0 h vs 24 h 中，GO 3 个组分中富集到较多靶基因的条目基本相同。但注释到的GO term 数目以及富集在各条目上的靶基因，8 h 多于24 h，这表明“滇北”在水淹胁迫8 h 响应比24 h更复杂；推测这些鸭茅差异miRNA 可能通过调控参与以上过程的靶基因的表达来影响鸭茅叶片发育及响应环境。

2.6 差异miRNA 靶基因KEGG 富集分析

在生物体内，不同基因相互协调行使其生物学功能，通过Pathway 显著性富集能确定差异表达基因参与的最主要生化代谢途径和信号转导途径。KEGG（kyoto encyclopedia of genes and genomes）是有关Pathway 的主要公共数据库［44］，该数据库提供的整合代谢途径（pathway）有助于进一步解读基因的功能，以Q-value＜0.05 的pathway 为差异表达显著富集的pathway。在富集散点图中Rich factor 指该pathway 中富集到的差异基因个数与注释基因个数的比值；Rich factor 越大，表示富集的程度越大；Q-value 是进行多重假设检验校正之后的P-value。

“滇北”在水淹胁迫0 h vs 8 h 有73 个差异表达基因注释到19 条pathway（图3），发现差异miRNA 的靶基因主要参与IL-17 信号通路（IL-17 signaling pathway）、Th17 细胞分化（Th17 cell differentiation）、雌激素信号通路（estrogen signaling pathway）、NOD-样 受 体 信 号 通 路（NOD-like receptor signaling pathway）、坏 死 性 凋 亡（necroptosis）、PI3K-Akt 信号 通 路（PI3K-Akt signaling pathway）、植物−病原体相 互作用（plant−pathogen interaction）、抗原处理和呈递（antigen processing and presentation）、内质网中的蛋白质加工（protein processing in endoplasmic reticulum）。

图3 “滇北”水淹0 h vs 8 h 差异miRNA 靶基因KEGG 富集Fig.3 KEGG enrichment scatter plot of differential miRNA gene KEGG after 0 h vs 8 h flooding in“Dianbei”

“滇北”在水淹胁迫0 h vs 24 h 有69 个差异表达基因注释到48 条pathway，前20 条如图4 所示，发现差异miRNA 的靶基因主要参与Th17 细胞分化（Th17 cell differentiation）、IL-17 信号通路（IL-17 signaling pathway）、雌激素信号通路（estrogen signaling pathway）、NOD-样受体信号通路（NOD-like receptor signaling pathway）、坏死性凋亡（necroptosis）、PI3K-Akt 信号通路（PI3K-Akt signaling pathway）、植物−病原体相互作用（plant−pathogen interaction）、抗原处理和呈递（antigen processing and presentation）、内质网中的蛋白质加工（protein processing in endoplasmic reticulum）。

图4 “滇北”水淹0 h vs 24 h 差异miRNA 靶基因KEGG 富集Fig.4 KEGG enrichment scatter plot of differential miRNA gene KEGG after 0 h vs 24 h flooding in“Dianbei”

3 讨论

水淹胁迫会导致植物根系缺氧，并引起代谢从有氧呼吸向厌氧发酵转变，从而导致基因转录和蛋白质合成的快速变化［39］。前期研究结果表明，在鸭茅对淹水胁迫的响应过程中，多个差异基因和生物学过程参与了水淹的响应，特别是活性氧（reactive oxygen species，ROS）的清除过程［33］，在本研究中，鸭茅水淹处理0、8 和24 h 共鉴定到208 个 差 异 表 达的miRNA，差 异 表 达的miRNA 主 要属 于miR166、miR159、miR167、miR396 和miR156 这5 个miRNA 家族，这些结果表明“滇北”鸭茅中广泛存在各类保守miRNA，可能对“滇北”适应水淹胁迫后的缺氧环境起着重要的调节作用。在植物对外界环境响应过程中，miRNA 通过对靶基因的调控来影响植物的生理功能［45］。基于对差异miRNA 进行靶基因预测及靶基因的GO 和KEGG 功能分析，发现这些靶基因主要参与细胞生理过程（cellular process）、代谢过程（metabolic process）、IL-17 信号通路（IL-17 signaling pathway）、Th17 细胞分化（Th17 cell differentiation）等植物逆境响应过程，表明“滇北”鸭茅差异表达miRNA 通过调节对环境适应相关靶基因的表达，影响不同时间点叶片的抗逆性。相关研究发现水淹胁迫能够影响植物miRNA 的形成，部分响应水淹胁迫的miRNA 在植物中被鉴定出来，这些miRNA 在水淹胁迫中的调控功能也得到证实［40］。

miR166 通过调控靶基因HD-ZipⅢ参与调控植物的生长发育过程，也可以作用于HD-ZipⅢmRNA 参与调节拟南芥逆境胁迫响应［46−47］。在玉米中，水淹胁迫下zam-miR166 的积累，影响根分生组织细胞的分化从而导致维管系统的形成，木质部组织扩张［48］。有研究表明，苜蓿（Medicago truncatula）中过表达的miR166 会导致侧根数量形成减少和根维管束的变化，并且miR166 是少数在低氧胁迫下降低的家族之一［49］。拟南芥在遭遇涝渍胁迫后，miR166 依赖于线粒体功能，在基因表达调控和低氧环境反应中发挥作用［50］。也有关于miR166 在其他逆境胁迫中的报道，如经高温处理后的不同品种蕹菜（Ipomoea aquatica），表现出截然相反的表达模式，是由于其在长时间高温后茎顶端分生组织受损造成［51］。水淹胁迫后“滇北”鸭茅差异表达的miRNA 中，miR166 数目最多，推测miR166 可能在植物遭受水淹胁迫中整合关键的缺氧信号途径，以此抵御该非生物胁迫。

miR167 可通过调控不同的靶基因参与到植物发育和逆境胁迫适应的不同阶段中。miR167 的靶基因是生长激素反应因子ARF基因家族中的ARF6和ARF8，生物学功能是增强生长素反应活性以及侧根生长，miR167 与其靶基因ARF6和ARF8之间存在反馈调控，ARF6激活miR167 的表达，ARF8抑制miR167 的表达，它们之间的相互激活或抑制调节着不定根的形成［52］。白羽扇豆（Lupinus albus）在低磷胁迫条件下，miR167 靶向调节转录因子增强生长素反应活性，促使侧根生长上调表达［53］。张丹凤［54］对水淹胁迫下玉米幼苗的适应机制进行研究发现，水淹胁迫影响玉米幼苗根系细胞中生长素（IAA）的动态平衡。IAA 含量的变化会诱导zma-miR167 的表达，zmamiR167 通过与AGO1相互作用，负调控靶基因ARF17的表达。ARF17抑制基因GH3的表达，进一步诱导下游基因ra2的表达，最终导致不定根的产生。“滇北”鸭茅在不同的处理时间段下，miR167 差异表达基因的数量为8 h少于24 h，这表明miR167 受水淹胁迫后，不定根的形成与鸭茅水淹胁迫密切相关。

miR156，miR159 和miR396 这3 个miRNA 是植物发育阶段过渡的关键调控因子［55］，在非生物胁迫中也有重要作用。miR156 是在植物中发现的第一个miRNA，其调节发育阶段转换的功能最初在拟南芥中得到证实，它通过靶向表达SPL基因来实现［56］。SPL是一种植物特异性转录因子，具有保守的SBP 结构域，可识别并特异性结合下游靶基因的启动子区域，从而调节植物的生长发育［57］。在拟南芥中，miR156b 的过表达会产生更多的叶片并且迫使开花延迟，同时伴随着多个SPL基因（SPL2/3/4/5/6/9/10/11/13/15）表达的下调［56］。miR156 在植物暴露于各种环境胁迫时异常表达，研究表明miR156 是植物热应激记忆所必需的，并且植物表现出对热应激的耐受性增强［58］。在水淹胁迫方面，CSA-miR156G-3p 在渍水的黄瓜（Cucumis sativus）下胚轴中下调，参与植物的非生物胁迫响应［59］。本研究中，8 h 差异表达的miR156 多于24 h，说明8 h 时鸭茅为降低水淹胁迫对植物生长发育的影响，通过表达更多相关基因来适应该非生物胁迫。

miR159 的靶基因编码MYB，Ma 等［55］研究发现miR159 的上调可以抑制MYB29 的表达，降低脂肪族脂肪硫代葡萄糖苷（giucosinolates，GSLS）水平，从而调节茎向根茎转化过程中的生长，以响应地上茎到根茎的转变。在玉米的根中，miR159 被水涝诱导上调，并负责沉默编码GAMYBs、MYB33 和MYB101 同系物的两个mRNA［22］。有研究表明，miR159 对各种环境胁迫有反应，与野生型相比，过量表达miR159 的转基因水稻（Oryza sativa）对热胁迫更敏感，这表明miR159 的下调可能有助于增强水稻对热胁迫耐受性［60］。Zhang 等［48］发现Zam-miR159 在淹水初期被抑制，24 h 后被诱导。miRNA 的预测涉及碳水化合物和能量代谢，包括淀粉合成酶、转化酶、苹果酸酶和ATP 酶，表明淹水响应miRNA 在转录后水平上参与调节玉米根系的代谢、生理和形态适应。水淹胁迫处理后miR159 明显变化且显著上调，说明其与植物耐涝性的发挥密切相关。

miR396 通过调控其靶基因生长调节因子（growth regulating factor，GRF）在植物不同发育阶段发挥多种作用。Zhang 等［61］通过对发育中的玉米籽粒进行RNA 测序，分析miR396 和GRF基因在玉米有效灌浆过程中的动态变化，发现在该过程中随着玉米灌浆进程的推进，miR396 高表达并逐渐下降，与其靶基因GRF呈显著负相关，且大多持续增加。miR396 是一种已知的应激反应miRNA，在盐和碱胁迫处理期间，结构性过表达miR396 的转基因拟南芥和水稻植株显著阻碍了根的生长，并影响植株的生长和发育［62］。在本研究中，miR156、miR159 和miR396 在水淹胁迫处理后高表达，这可能是由于水淹胁迫限制了鸭茅正常生长发育的关键调控因子，意味着“滇北”鸭茅需要更多更强的抗性基因来维持植物的存活。

植物在遭遇水淹胁迫时会导致基因的表达变化，可通过生物信息学分析对DEGs 进行功能预测，探究植物耐涝机理和分子调控网络。Xu 等［59］通过研究黄瓜不定根的形成对内涝胁迫的影响，GO 富集分析表明，在细胞氧化还原稳态、细胞骨架、光合作用和细胞生长等方面功能的靶基因被过度表达，说明细胞相关生物过程在响应水淹胁迫中发挥重要作用。林延慧等［63］通过转录组和蛋白组的关联分析，在大豆（Glycine max）中筛选出3 个与水淹胁迫丝裂原活化蛋白激酶（mitogen-activated protein kinase，MAPK）信号通路相关的差异表达基因，分别是Glyma. 07G127100、Glyma. 05G124000 和Glyma. 05G123700，发现水淹胁迫引发了植物类似于病原菌侵染的基本防御反应，传递信号通过FLS2 蛋白进一步激活了MAPK 信号通路来抵御涝害胁迫。

4 结论

本研究利用小RNA 测序方法，对水淹胁迫处理后“滇北”鸭茅miRNA 在不同时间点的表达模式及其靶基因进行了分析。初步发现，miRNA 在水淹胁迫不同时间点表现出一定的动态变化。同时，根据miRNA 序列及靶标功能预测，推测鸭茅miR166、miR159、miR167、miR396 和miR156 很可能参与了鸭茅的水淹胁迫应答过程。这些结果为鸭茅中与水淹胁迫相关miRNA 的进一步研究以及探究miRNA 在鸭茅抗涝过程中的潜在功能和作用机理提供了参考。但有关响应水淹胁迫的miRNA 如何调控其靶基因，预测的耐涝相关靶基因又如何发挥其作用等问题还需要进一步探究。