PEG和NaCI胁迫下毛竹萌发种子的MicroRNAS表达谱分析

王晓静 王涛 杨凯 李潞滨

关键词：毛竹；种子萌发；PEG；NaCI；miRNA

毛竹（PhyHostachys edulis （Garr.）H.deLehaie）属于禾本科（Gramineae）竹亚科（Bambusoideae）刚竹属（PhyHostachys Sieb.Et Zucc.），是一种多年生木本散生竹，广泛分布在460 N以南地区，是我国种植面积最广的竹种之一，具有重要的经济、生态和文化价值。毛竹开花后产生典型颖果类果实，薄而干的膜质果皮与种皮紧密贴合，在种子生物学和苗木培育上毛竹颖果通常被称为种子。毛竹种子萌发是实生苗造林的基础，通过毛竹种子萌发进行实生苗造林能降低运输成本，便于移栽，有利于提高造林整齐度、成活率和遗传多样性，延长竹林高产时间，且实生苗竹林萌笋多、竹鞭生长较快、竹材利用率更高，因此，开展毛竹种子萌发研究对于竹类资源的培育和应用等具有重要意义。

随着气候变化的加剧，干旱和由土壤初级或次级盐碱化造成的盐胁迫已经成为全球性问题，并对种子萌发、植物生长等造成严重影响。干旱和盐胁迫抑制了毛竹种子的萌发，进而影响了毛竹实生苗的生长，限制了通过毛竹种子萌发进行实生苗培育和应用。目前，研究者逐步开展了干旱或盐胁迫对毛竹种子萌发率、种子活力、种子寿命和萌发后幼苗生长等生理方面的研究，但相关的分子调控机制鲜有报道。

MicroRNA（miRNA）是一类长度约为21～24nt的小分子非编码RNA，能够通过互补配对抑制或降解靶基因进行转录后调控。miRNA广泛参与了植物的生长发育，在种子萌发和非生物胁迫响应中具有重要的调控作用。在竹类植物中，研究者已经在组学水平开展了关于毛竹叶片、花发育、雄蕊发育、茎秆快速生长和竹笋发育等的miRNA研究，但毛竹种子中miRNA的表达及其在干旱和盐胁迫下萌发时的调控仍是未知的，因此，本研究通过small RNA测序对毛竹种子露白阶段的miRNA进行系统分析，并探究其在不同PEG或NaCI胁迫下萌发时的表达模式，挖掘毛竹种子萌发阶段调控干旱或盐胁迫抗性潜在的关键miRNA，以期为后续的研究提供基础和参考。

1材料与方法

1.1试验材料

试验所用毛竹种子于2020年9月收集自广西壮族自治区桂林市灵川县，剥除外稃后挑选大小均一、颗粒饱满、色泽明亮的完整种子备用（图1a）。种子的千粒质量为21.85g±0.04g，含水量为12.05%±0.35%。种子表面消毒后在培养皿中使用纸床法进行萌发。聚乙二醇（PEG6000）是一种亲水性大分子物质，是模拟植物干旱胁迫的常用试剂之一，本研究使用PEG6000模拟干旱胁迫环境，同时使用氯化钠（NaCI）模拟盐胁迫。使用完全随机实验设计，将种子分为5组（A～E）：A组为对照组，B组、C组模拟干旱胁迫，D组、E组模拟盐胁迫。参考前期预实验的结果，分别在培养皿中加入H20（A）、10% PEG（B）、15%PEG（C）、50mmoI.L-1 NaCI （D）和100mmoI.L-1NaCI（E）．每组萌发200粒种子，并重复3次。每天观察并记录种子萌发状态，在第4天选择生长状态一致且均到达种皮破裂阶段的种子（图1b）。对于每个处理，在3个重复中各取5粒种子并将15粒种子混合取样后用于建库。取样后使用液氮速冻并置于-80℃保存。

1.2建库、测序与数据过滤

使用植物多糖多酚RNA提取试剂盒（DP441，TIANGEN）提取样品总RNA。在建库时首先根据small RNA的结构特点将总RNA中small RNA加接头并进行纯化，随后进行反转录、PCR扩增和聚丙烯酰胺凝胶电泳，切胶回收目的条带即得到small RNA的cDNA文库。文库有效浓度合格后进行测序，使用lllumina高通量测序平台NovaSeq6000进行单端测序，测序长度为50bp。对获得的原始数据首先进行过滤，去除接头序列、低质量序列、含N比例>10%的序列、5接头污染的序列、没有3接头序列和插入片段的序列及含连续的/T/G/C的序列；以及长度异常的序列。

1.3miRNA鉴定

首先将过滤后的数据比对到miRbase数据库（V22）鉴定已知miRNA；将未比对至miRBase的数据比对至Rfam数据库中的ncRNA序列（rRNA、tRNA、snRNA、snoRNA）和参考基因组重复序列；随后对于上述均未比对中的cleanreads，使用mirEvo和miRdeep2进行新miRNA预测。对于已知和新的miRNA序列，将其前体序列比对至miRBase数据库进行microRNA家族的鉴定。

1.4表达水平分析

为了探究毛竹种子萌发露白阶段miRNA的表达，首先统计各样本中已知和新miRNA的readscount．进行TPM（Transcripts per million）归一化转化（TPM=readscount×107／library size）。随后使用edgeR软件进行差异表达分析，以鉴定露白阶段毛竹种子响应干旱或盐胁迫的miRNA，差异表达miRNA筛选条件为llog2foldchangel≥1且p<0.05。

1.5靶基因预测及功能富集分析

为了探究毛竹萌发种子中miRNA的潜在功能，对本研究中已知和新miRNA的靶基因进行预测，并进行靶基因的KEGG和GO功能富集；同时通过对差异表达miRNA靶基因的功能富集，进一步分析毛竹萌发种子中miRNA响应干旱或盐胁迫的潜在通路或途径。使用软件TargetFinder和psRNAtarget进行miRNA靶基因的预测，并对预测结果取交集。使用OmicShare在线平台（（https：／／www．omicshare．com／tooIs））进行靶基因KEGG和GO功能富集，富集显著性阈值为p<0.05。

1.6差异表达miRNA验证

总RNA提取方法同1.2，反转录使用miRcute增强型miRNA cDNA第一链合成试剂盒（KR211，TlANGEN），荧光定量检测使用增强型荧光定量试剂盒（FP411，TIANGEN）。荧光定量正向引物见表1，反向引物使用试剂盒通用引物，内参引物为U6。荧光定量检测仪器为qTOWER3.0（AnaIytik jene），每个基因进行3次生物学重复和3次技术重复，采用2方法计算相对表达水平。

2结果与分析

2.1small RNA测序结果

small RNA测序结果表明：每个文库产出数据约0.9～1.8GB，过滤后每个样本中分别有12984751条（A， H20）、14982423条（B，10%PEG）、13401719条（C，15% PEG）、13995450条（D，50mmoI.L-1NaCI）和11452796条（E，100mmoI.L-1NaCI）Clean Reads（表2）.

CleanReads的Q20值约为98%，Q30值均≥93%，与毛竹参考基因组的比对率≥88.42%，表明测序数据质量较好。具体数据产出及比对情况见表2。

2.2miRNA鉴定与miRNA家族分析

对small RNA数据库中miRNA进行比对和预测，结果显示：每个样品中有0.08%～0.11%的序列被注释为已知miRNA，0.10%～0.11%的序列被预测为novel miRNA。对miRNA长度进行统计发现本研究成熟miRNA长度主要分布在21nt和24nt（图2）。本研究共鉴定miRNA成熟体序列508条，前体序列845条，其中，已知miRNA成熟体序列数量为246，前体序列数量为574，鉴定出novel miRNA成熟体数量共有262条，前体序列数量为271。

对miRNA前体序列进行家族分析，共鉴定到包含514条前体序列的45个家族，包括MIR159、MIR166、MIR156、MIR408、MIR399、MIR530等，其中，有10条新miRNA前体序列也被鉴定属于已知家族，如novel 57属于MIR408，novel 14和novel 66属于MIR169 2。统计每个家族包含的前体成员数目表明，本研究毛竹萌发种子中最大的miRNA家族为MIR159，包含54个成员，其次为MIR166 （48）、MIR156 （46）、MIR167 1（37）、MIR396（32）等（详细数据未列出）。

2.3miRNA表达水平分析

毛竹萌发种子中表达量前5的已知miRNA分别为phe-miR166a-3p、phe-miR159a.1、phe-miR319a-3p.2-3p、phe-miR6478、phe-miR156a-5p，表达量前5的新miRNA分别为novel 1、novel 199、novel 311、novel 241、novel 2，其中，novel 1在5个样本所有miRNA中表达量均为最高。表达量前10的新miRNA折叠最小自由能值为-187.0～-13.8

KJ·mol-1（表3、4）。

对每个样本中表达水平前10的已知miRNA和新miRNA的进行TPM值统计，结果表明：排名前10的已知miRNA在每个样本总miRNA中的读序丰度（total miRNA reads）占45%～51.8%，排名前10的新miRNA在每个样本总miRNA中的读序丰度中占23.1%～29.5%。PEG和NaCI胁迫下这些miRNA在毛竹露白阶段种子表达水平均较高，推测其在毛竹种子萌发的调控中可能具有保守的重要作用。

2.4差异表达分析

通过差异表达分析，本研究共鉴定181个差异miRNA，其中，有84个上调表达，97个下调表达（图3a）。与对照组（A）相比，10%PEG（B）、15%PEG（C）、50mmoI.L-1NaCI（D）、100mmoI.L-1NaGI（E）

4种胁迫下分别有20、26、41和24个差异表达miRNA（DEmiRNA），此外在B-C和D-E比较组分别有16和54个DEmiRNA（图3b）。

依据miRNA的差异表达情况和表达水平，对本研究具有高丰度且显著差异表达的miRNA进行聚焦。与对照组相比，10% PEG、15% PEG、50mmoI.L-1 NaCI、100mmoI.L-1 NaCI中表达水平最高的DEmiRNA分别为novel 14、novel 311、novel 14、phe-miR159a.1，表达水平最高的已知DEmiRNA分别为phe-miR171e-5p、phe-miR3630-3p、phe-miR171e-5p和phe-miR159a.1，这些miRNA在毛竹种子萌发露白阶段中大量积累，同时在PEG或NaCI胁迫下差异表达，可能参与了种子萌发露白阶段miRNA对干旱或盐胁迫的调控。

2.5miRNA靶基因预测

通过Targetfinder预测，共获得505个miRNA的31729个靶基因，psRNAtarget预测结果中501个miRNA可以靶向22215个基因。对2个软件的预测结果取交集后共有17666对miRNA-靶基因关系，包含489个miRNA和8812个靶基因。预测的miRNA-靶基因关系中，novel 153的靶基因数目最多，有188个靶基因，其次是phe-miR396b（167个）、phe-miR396e-5p（161个）、phe-miR396h （141个）、phe-miR164b（138个）.对同- miRNA家族成员的靶基因数目进行统计，靶基因数目前十的家族中MIR396家族靶基因数目最多。此外，根据预测结果8812个靶基因中每个基因能够受到1～28个miRNA靶向，表明在毛竹种子萌发过程存在复杂的miRNA调控网络。具有miRNA-gene关系最多的基因中，排名前20的基因有7个属于GAMYB（GibberellinM、YB）家族，飞2个属于SPL（Squamosapromoter-binding-like）家族，推测GAM、YB和SPL基因家族在毛竹萌发种子miRNA调控网络中可能具有重要意义。

2.6差异表达miRNA靶基因功能富集分析

根据靶基因预测结果，A-B、A-C、A-D、A-E、B-C、D-E 6个比较组差异表达miRNA分别有609、548、1011、812、460和1866个靶基因。GO富集结果显示：6个组合差异表达miRNA的靶基因分别显著富集在520（A-B）、516（A-C）、620（A-D）、446（A-E）、483（B-G）、610（D-E）个条目中，总计1989个涉及细胞成分、生物学过程和分子功能的条目，推测miRNA可以通过这些靶基因广泛地参与不同的生物学过程或分子功能等调控。表5为不同比较组GO富集前10的生物学过程条目，与对照组相比，10% PEG（B）、15% PEG（C）．50mmoI.L-1 NaCI（D）、100mmoI.L-1 NaCI（E）胁迫下DEmiRNA靶基因富集最显著的生物学过程分别为跨膜转运（GO：0055085）、DNA复制正向调控（G0：0045740）、发育过程（G0：0032502）、海藻糖生物合成（GO：0005992），推测这些生物学过程对于miRNA调控PEG和NaCI胁迫下的种子萌发具有重要作用。

KEGG富集结果（表6）表明：苯丙烷生物合成途径在A-B、A-C、A-E 3个组合中被显著富集，推测苯丙烷生物合成途径对于毛竹种子萌发露白阶段miRNA响应PEG和NaCI胁迫可能具有重要意义，但毛竹种子萌发过程miRNA与苯丙烷合成途径关键基因的靶向关系及其是否共同调控干旱或盐胁迫下的种子萌发需要进一步的研究和验证。此外，果糖和甘露糖代谢、玉米素生物合成途径仅在盐胁迫下（A-D和A-E）显著富集，油菜素内酯生物合成和脂肪酸降解途径仅在A-B中显著富集，表明在毛竹种子露白阶段miRNA响应10%PEG、15% PEG、50mmoI.L-1NaCI、100mmoI.L-1 NaCI胁迫的调控通路也存在差异。

2.7差异表达miRNA qPCR验证

随机选取10个在不同样本间存在差异表达的miRNA进行qPCR验证，miRNA的表达量在不同处理组之间被上调或下调，如phe-miR6478在10% PEG胁迫下比对照组上调1.45倍，phe-miR3630-3p在100mmoI.L-1 NaCI胁迫下上调1.33倍，phe-miR171e-5p在10% PEG、15% PEG、50 mmoI.L-1 NaCI、100mmoI.L-1 NaCI胁迫下分别下调0.23、0.50、0.56和0.61倍（图4）。miRNA荧光定量结果和small RNA测序数据在表达趋势上整体一致，表明测序数据的可靠性。

3讨论

在本研究中，miR166、miR159、miR319、miR156、miR396、miR167、miR168、miR894、miR160等在对照组和不同浓度的PEG和NaCI胁迫条件下均具有较高表达水平，推测这些miRNA可能在毛竹种子萌发中有重要作用，其中，miR159、miR156、miR396、miR160已经被证明参与了种子萌发调控。在大麦（Hordeum vulgare L）萌发种子中，miR156、miR166、miR167、miR168等也具有高表达水平，同时miR5071被大量积累，而在本研究中未检测到miR5071的表达，表明不同物种间miRNA对种子萌发的调控具有特异性。在已有研究中，miR402、miR163和miR417被证明能够参与种子萌发期的干旱或盐胁迫响应，而在本研究对照和处理组中均未检测到miR402、miR163和miR417的表达，推测这些miRNA在本研究中可能不是响应干旱或盐胁迫的关键miRNA，其分布和表达可能与物种有关。

与对照组相比，本研究中miRNA主要在单一胁迫下显著差异表达，仅有2个miRNA能够同时响应2种胁迫，表明miRNA在10% PEG、15%PEG、50mmoI.L-1 NaCI、100 mmoI.L-1 NaCI4种胁迫下的调控存在差异。本研究丰度较高且显著差异表达的miRNA中，phe-miR171e-5p能够响应10% PEG和50mmoI.L-1 NaCI胁迫，phe-miR3630-3p在15%PEG胁迫下显著上调表达，phe-miR159a.1在100mmoI.L-1 NaCI胁迫下差异表达，推测这3个miRNA对于毛竹种子萌发响应PEG或NaGI胁迫可能具有重要意义。

前人的研究中，PEG胁迫下甘蓝型油菜（Brassica napus L）萌发种子中的miR171比对照组显著下调表达；此外，在拟南芥（Arabidopsisthaliana （L.） Heynh.）中过表达桑树（Morus albaL_） mno-miR171可提高转基因植株在NaCI和甘露醇胁迫下的种子发芽率，这与本研究phe-miR171e-5p能够响应PEG和NaCI胁迫是一致的。miR171主要通过靶向GRAS家族基因影响GA和ABA信号通路，参与胚胎发生潜能的维持、花药发育、芽分枝和复叶形态的调控、顶端优势的调控等，其是否通过GRAS家族参与毛竹种子萌发过程中对干旱和盐胁迫的调控需要进一步研究。在苜蓿（Medicago sativaL．）中，miR3630被鉴定为响应干旱胁迫的miRNA，本研究中phe-miR3630-3p在15% PEG胁迫下显著上调表达，但目前关于miR3630的研究较少，本研究为miR3630响应干旱胁迫提供了新的证据。此外，phe-miR159a.1在A-E比较组显著下调。拟南芥中miR159通过介导GAMYB家族MYB101和MYB33的转录本切割，参与种子萌发过程糊粉层细胞的细胞程序性死亡，推测miR159对于100mmoI.L-1 NaCI胁迫下毛竹种子的萌发也可能具有重要的调控作用。

4结论

本研究首次系统鉴定了毛竹萌发种子中的miRNA，并对其在PEG和NaCI胁迫下毛竹萌发露白阶段种子中的表达模式进行研究，探究了响应不同PEG和NaCI胁迫的差异表达miRNA。下一步将对本研究中涉及的重要miRNA进行靶基因的验证，并对其调控机制和调控功能进行深入研究和探讨。