深度测序鉴定玉米雄穗花器官分化期响应干旱胁迫的miRNA和其靶基因

王业建，梁晓玲，阿布来提·阿布拉，韩登旭，杨 杰，郗浩江，刘 俊，李铭东

(新疆农业科学院粮食作物研究所， 乌鲁木齐 830091)

0 引 言

【研究意义】干旱能抑制植物的生长和发育，并导致植物的品质和产量下降。了解植物对非生物胁迫的反应对提高作物生产力有重要意义。【前人研究进展】在长期的进化过程中，植物已经开发出一系列调节机制，以应对不同层面的这些不利条件，包括细胞、生理、生物化学和分子过程[1-2]。已确定植物中激素介导的信号交叉参与对干旱和盐胁迫的反应，如脱落酸(ABA)、水杨酸和乙烯。基因表达调控是植物在转录后水平对抗干旱的关键策略[3-4]。基因表达调控最重要的参与者之一是microRNA(miRNA)。miRNA是一种长度约为21个核苷酸的非编码小RNA分子。研究表明，miRNA具有与其靶mRNA的完美或接近完美的互补核苷酸结合的特点，通过mRNA切割或通过翻译抑制负向调节靶基因的表达[5]。除了miRNA在发育和代谢中的作用外，miRNA还参与非生物和生物应激反应。例如，miR394是一种保守的miRNA，已在一系列植物物种中发现，如拟南芥、水稻、棉花、柳枝稷和甘蓝型油菜[6-7]。目前，有多达40个与非生物胁迫相关的植物miRNA家族，其中多种与干旱胁迫反应有关[8]。玉米(Zeamays)基因组中miRNA的表征和分析取得了重大进展[9-10]。虽然玉米是一种相对耐旱和耐盐的作物，但玉米高盐度和过度缺水可导致一系列代谢紊乱，包括渗透作用(脱水)，营养不平衡和离子毒性，都会产生极大的负面影响[11]。基于转录组和转基因分析，已显示大量基因在玉米或其他物种中响应盐度和干旱胁迫而显示出持久的表达。例如，转基因拟南芥中玉米CBL相互作用蛋白激酶基因(GhCIPK6)的过表达导致对盐，干旱和ABA胁迫的耐受性提高，GlyCIP6可能是对抗玉米中盐和干旱胁迫的正向调节因子[12]。转基因烟草表达玉米C群丝裂素活化蛋白激酶(MAPK)基因(GhMPK2)具有较低的水分损失率，表现出对盐和干旱的耐受性增强，GhMPK2可能正向调节烟草和玉米的盐和耐旱性[13]。基于微阵列的转录组分析揭示了玉米中某些盐/干旱介导的信号转导途径，其中许多候选基因表达差异，可能是耐受干旱胁迫的潜在标记，如WRKY，ERF，跨膜转运蛋白，丙酮酸脱羧酶和蔗糖合成酶等[14]。【本研究切入点】与亲本相比，近交亲本系的F1自交后代表现出优越的性能和胁迫耐受性[15-17]。这种效应被称为杂种优势，在植物育种中被广泛利用。研究运用现代分子生物学及高通量测序技术，测序鉴定玉米雄穗花器官分化期响应干旱胁迫的miRNA及其靶基因。【拟解决的关键问题】通过miRNA测序，筛选出耐旱自交系“PHBA6”和干旱敏感自交系“吉63”之间，以及PHBA6干旱胁迫和对照处理之间差异表达的miRNA，挖掘调控耐旱自交系“PHBA6”雄穗花器官分化期干旱适应性的miRNA。

1 材料与方法

1.1 材 料

材料为耐旱自交“PHBA6”(PHZ51×PHG47)和干旱敏感自交系“吉63”[(127-32×铁84)(W24×W20)辐)] , 由新疆农业科学院粮食作物研究所提供。

材料于2017年在新疆农科院海南三亚科技示范园基地大棚分期播种，前期进行正常的滴灌浇水处理，2个自交系的雄穗处于花器官分花期时进行干旱胁迫处理，耐旱种质“PHBA6”进行正常滴灌浇水处理作为对照，干旱处理15 d后，分别取干旱处理的“PHBA6”和“吉63”的雄穗(分别命名为PDST15和JDST15)，以及对照处理的“PHBA6”的雄穗(命名为PNTT15)，迅速包好放入液氮中, 于-80℃保存。

1.2 方 法

1.2.1 总RNA提取及质量控制

分别合并来自6个对照和PEG处理的玉米雄穗，使用miRNA Mini Kit(Qiagen)分离小RNA。 使用Agilent 2100生物分析仪(Agilent)测试RNA样品的质量和浓度。 将小RNA与3'和5'适配子连接，并使用TruSeq TM Small RNAkit(Illumina)根据说明书进行RT-PCR。 从凝胶中分离预期的最终PCR产物，纯化并使用Hiseq2000测序仪(Illumina)测序。sRNA测序和基本生物信息学分析分别在MAGrogen(韩国)和LC Sciences(Texas, USA)上进行。

1.2.2 文库制备和miRNA测序

用SolexaQA软件包(Solexa)计算质量统计数据，采用Illumina HiSeq 2500仪器进行深度测序得到50 bp 的单端测序读数(reads)。来自4个库中每个库的生成数据文件用于后续分析。从miRNA读数中去除适应序列后，建立小RNA的总读数和独特读数。 从读数中，映射到rRNA，tRNA，snoRNA和snRNA的小RNA被丢弃，剩余的读数被映射到玉米基因组和转录本中。这些基因组合转录本数据可使用Bowtie软件在美国国家医学图书馆(NLM)的Genbank(http://ftp.ncbi.nlm. nih.gov)和法国的膜生物学国家重点实验室获得(http://www.biomemb. cnrs.fr/ date_ contigs.doc)。为了使用Bowtie进行对齐，需要去除缺乏3’的序列、5’污染的序列和小于18个核苷酸的序列。

1.2.3 保守和新miRNA的预测

为了鉴定保守的miRNA，使用Bowtie v 0.12.7将sRNA与NCBI中发表在miRBase 20.0(http：/ / www.miRBase.org)植物成熟miRNA比对。具有完整匹配序列的miRNA被认为是保守序列。具有1～3个错配的miRNA被认为是miRNA变体，并且与miiRNA数据库中没有相似性的miRNA被认为是新的miRNA。 在该研究中鉴定的miRNA前体序列的长度范围为57～264 bp。 使用UAEsmall RNA workbench软件验证了二级发夹结构的形成和推定的新miRNA序列在相应前体上的比对。

1.2.4 miRNA表达分析比较和功能

按照先前描述的方法，用miRNA Read/Clear×Total Read Count来进行表达量归一化处理，即计算每个miRNA中的序列计数的标准化。文库大小参数被认为是特定样本和伪参考样本的比率之间的中值。在去除了文库中未显示读数的miRNA序列后，在文库中技术miRNA独特序列。使用数字基因表达数据软件(degseq包)进行差异表达分析。将变异的生物系数(BCV)值设定为0.2。如果多个测试校正的P值并且错误发现率(FDR)<0.05，则预测样品池间miRNA表达是差异显著的。使用PermutMatrix软件和KOBAS软件进行层次聚类分析和KEGG 通路的富集分析。

1.2.5 miRNA靶基因的计算预测

使用psRNAtarget软件在RNA-seq重叠群上搜索miRNA来完成对从4个文库测序的保守的，变体和新的miRNA的靶基因的预测。参数设置为：最大期望值(ME)为3，完成度评分(hspsize)为20，目标位点周围的侧翼长度(上游17 bp、下游13 bp和中心错配区9～11 nt的翻译抑制为靶基因可分析范围。

注释潜在的mRNA靶标，并使用blast2GO v2.3.5软件注释功能，基于推断的蛋白质相似性，使用BLASTx在Swiss-Prot / Uniprot蛋白质数据库中获得的其他序列和保守结构域。每个序列被注释一个基因本体(GO)或更的细胞组分，分子功能和生物过程。

2 结果与分析

2.1 测序数据对比

研究表明，分别从JDST15、PDST15和PNTT15的miRNA文库中获得12133269、14669833和1009941 reads。去除低质量的数据后，共剩下9632540、550036和6868703条有效读取(valid reads)被用miRNA比对鉴定与预测分析。样本中测序序列的分布概况，其中一些低比重占有率的大片段序列，例如rRNA或mRNA，样本的总RNA质量较高未被降解。3个文库都显示出与其他RNA家族相似的分布，包括rRNA(Uniq～1.34%)，sn RNA(Uniq～0.02%)，snoRNA(Uniq～0.01%)和tRNA(Uniq～0.13%)。表1

选取Rfam数据库来注释测序得到的小RNA序列，尽可能地发现并去除其中可能的rRNA、 scRNA、snoRNA、snRNA、tRNA等非miRNA序列。分别对总的测序数据(总数)以及Unique数据(种数)进行了统计。选取Rfam数据库来注释测序得到的小RNA序列，尽可能地发现并去除其中可能的rRNA、 scRNA、 snoRNA、 snRNA、tRNA等非miRNA序列。分别对总的测序数据(总数)以及Unique数据(种数)进行了统计。PNTT15组特异性rRNA、tRNA、snoRNA和snRNA所占比例分别为：54.50%、26.74%、6.39%和6.62%；PDST15组分别为为42.09%、26.09%、18.43%和10.97；JDST15组分别为53.12%、24.49%、15.58%和7.50%。图1

图1 Rfam数据比对Fig.1 Rfam data comparison results

2.2 Valid小RNA长度分布

研究表明，统计过滤后的Valid数据的总数及种数长度，大部分数据分布在20～24 nt，符合Dicer 酶切割的典型特征。图2

图2 Valid数据长度分布统计Fig.2 Valid data length distribution statistics

2.3 PDST15 vs PNTT15和JDST15 vs PDST15差异miRNA的检出和发现

研究表明，在玉米中已经报道的前体miRNA检出176个。其中在JDST15、PDST15和PNTT15组中分别检出160、151和145个。成熟体miRNA检出199个在3组中分别检出156、159和146个。在JDST15、PDST15和PNTT15组分别新发现前体miRNA 83、81和50个，成熟体miRNA7、76和42个。

根据miRBase(Release 20)(Kozomara和Griffiths-Jones，2011)中所有已知植物miRNA的3个文库的所有清洁读数的比对发现，miRNA在物种间高度保守，共有709个已知植物miRNA家族。这些miRNA家族占总独特读数序列的约0.47%和平均总冗余阅读序列的23.59%。在这些miRNA家族中，71个miRNA在2个不同物种中响应干旱有特异型。47个miRNA家族仅在PHBA6自交系对照处理中发现。例如，miR1868和miR2099分别仅在干旱处理的样品中表达。干旱处理的文库共享65个miRNA家族，这些家族没有出现在对照组文库中。在3个文库中共鉴定出709个miRNA家族中的357个，表明在维持正常生物活性中的关键作用，例如miR156 / 157、miR159、miR168和miR172。所有3个文库都存在相似的最常见miRNA家族，包括miR156，miR157、miR166、miR167和miR3954。 Pearsonχ2检验显示，709个(79.69%)miRNA家族中有565个在3个文库中差异表达(P值≤0.05)。共有443个(61.37%)miRNA家族在三种处理的成对比较中表现出显着差异，包括miR157，miR159，miR2948和miR3694。 miR1854和miR1148在对照与干旱和盐与干旱之间具有最大的倍数变化，高达≥10倍，干旱胁迫强烈抑制其表达。表1～3

表1 3个文库中总miRNA标签分类Table 1 The classification of total small RNA tags in 3 library

2.4 miRNA靶基因预测和富集性

研究表明，共预测到4 658个基因，利用Fisher 精确假设检验对差异miRNA进行分析，得出miRNA和基因富集的KEGG通路和GO信号通路。在PDST15 vs PNTT15组中这些基因主要参与生物学过程中的蛋白磷酸化，转录本的调节，ARF蛋白信号转导的调控；细胞组成中的维持细胞膜和核完整性；分子过程中的蛋白结合，ATP结合。在JDST15 vs PDST15组中，差异基因参与的生物学过程主要为：转录的调节，脂肪酸生物合成过程和蛋白磷酸化；细胞组成中的维持细胞膜和核完整性；分子过程中的蛋白结合，ATP结合以及DNA结合。

差异表达miRNAs靶向基因富集的通路主要为：碳水化合物形成，脂质形成，氨基酸代谢产物，辅助因子和维生素生成，其他次生代谢的生物合成，聚糖生物合成和环境适应等。图3～5

注：PHBA6干旱胁迫和对照前50种最丰富的保守miRNA(左)，不同品种玉米干旱胁迫后前50种最丰富的保守miRNA的热图。红色，上调; 绿色，下调

图4 通过差异表达的miRNAs预测的靶基因注释的部分基因GOFig.4 GO analysis of partial genes annotated by target genes predicted by differentially expressed miRNAs

图5 预测差异表达miRNAs靶向的前10个途径Fig.5 Top 10 pathways for predicting differentially expressed miRNAs

3 讨 论

在干旱胁迫下，许多保守的和新的miRNA表达存在差异；一些miRNA甚至在干旱处理中特异性表达。在3个文库中获得了大量的miRNA读数，同时在3个文库中观察到冗余读数，独特读数和匹配的独特读数的类似大小分布，其中长度为24个核苷酸的读数占大多数。匹配的冗余读数在21个核苷酸类中具有最多的读数，其次是24个核苷酸。玉米中的小RNA丰度和大小与拟南芥[18]和水稻[19]中报道的结果大致一致。

在所有3种处理中，miR156/157家族是最丰富的小RNA，占总小RNA含量的60%。与对照相比，miR156和miR157在干旱处理中分别下调0.44和1.23倍。据报道，玉米根和叶中的miR156和miR157均在高浓度盐(> 2.5%)下表达量均下降，并且随着PEG的增加以剂量依赖性方式下调。目前的在干旱胁迫下雄穗中的结果也与前人报道基本一致。研究表明，miR156 / 157在植物中负向地靶向SPL转录因子，并且miR156 / 157过表达导致拟南芥中成年性状和开花的延迟[20]。miR156(Corngrass1，Cg1)的过量产生导致玉米幼年营养期的延长。目前关于miR156 / 157的研究主要集中在其在形态变化和开花调控中的作用[21]。研究中，小RNA测序提供了干旱胁迫干扰miR156 / 157表达的证据，miR156 / 157在干旱胁迫下的新作用。

大豆中NF-Y复合物的NF-YA(GmNFYA3)可被ABA和非生物胁迫诱导，包括干旱，NaCl和冷刺激[22]。拟南芥中GmNFYA3的过量表达导致叶片水分流失减少和干旱胁迫增加，并提高其对高盐度和外源ABA的敏感性。烟草中的体内试验显示miR169指导GmNFYA3 mRNA切割[23]。在玉米中，预测NF-YA3(contig16841和contig22907)是ghr-miR169的靶标。此外，与对照相比，干旱和处理中的ghr-miR169a / b / c表达显着下调0.04至0.93倍。相反，ghr-miR169d / e / f / g在干旱处理中显着上调0.04～0.85倍。Ghr-miR169i / j表达在干旱处理中受到抑制。推测至少ghr-miR169a / b / c可能通过作用于玉米中的NF-YA3而在抗旱中发挥积极作用。

植物水通道蛋白是大型主要内在蛋白家族的一类，众所周知，通过不同的小分子(包括水和其他小营养素)的生物膜发挥作用[24]。除了在水和养分的吸收和运输中的作用外，还发现水通道蛋白涉及一系列非生物胁迫，例如干旱和冷胁迫。例如，转基因烟草中质膜水溶蛋白的过表达在有利的生长条件下提高了植物活力，但在干旱或盐胁迫下没有，因为通过质膜水通道蛋白的共转运水在干旱或盐胁迫期间具有有害作用[25]。已发现从玫瑰中分离出的TIP型水通道蛋白基因Rh-TIP1受到乙烯和水缺乏处理的抑制[26]。确定20个miRNA可能靶向玉米中的22个水通道蛋白，其中两个miRNA-tar-get对(ghr-miR4371和contig780;以及ghr-miR4371和BK007054.1)也通过降解组序列分析检测到。水通道蛋白可能通过作为miRNA靶标(例如ghr-miR4371)参与对玉米中的干旱胁迫的响应。

过量表达miR393的转基因水稻植物对盐和碱处理更敏感，miR393是通过靶向非生物相关基因对盐和碱胁迫响应的负调节因子[27]。在该研究中，获得了相反的结果，其中ghr-miR393a / b / c / d / e表达通过干旱和盐度处理上调。然而，玉米中预测的ghr-miR393靶标也是应激相关基因，以及几种激素应答基因，包括NADPH：细胞色素P450还原酶(CPR1)，III类过氧化物酶(POX4)，蛋白质AUXIN SIGNALING F-BOX 3，和TIR1(运输抑制器响应1)。如果玉米中的miR393也是水稻中的负压胁迫调节因子，一种可能的解释是，玉米幼苗中miR393s在干旱胁迫下的上调有助于扩大玉米的胁迫信号，触发通过信号转导串扰来对抗干旱胁迫的更有效或更强大的途径，可能是通过生长素相关途径。

使用基于CitationRank的算法，在植物的干旱和盐度胁迫的背景下编码相关的编码基因。基于GO分析允许确定基因属于哪些GO术语(生物过程，分子功能和细胞成分)[28]。基于GO的分析可以提供更多关于理解miRNA功能的想法。研究中总共274个miRNA(256个保守的miRNA和18个新的miRNA)和其1验364个靶标被分类为542个分子功能，679个生物学过程和138个细胞组分。至少137种miRNA及其与应激反应相关。229种靶标能够分为1.2种分子功能，118种生物过程和37种细胞成分。35个miRNA-靶标对属于响应于干燥的生物学过程(GO：0009269)和对水剥夺的反应(GO：0009414)，例如ghr-miR159b、ghr-miR166h、ghr-miR399i、 ghr-n26和ghr-miR399f。许多分类的生物过程与信号转导有关，如生长素代谢(GO：0009850)和生物合成论(GO：0009851)，乙烯介导的信号通路(GO：0009873)，对生物刺激的反应(GO) ：0009607)，细胞分裂素代谢(GO：0009690)和生物合成(GO：0009691)，以及茉莉酸代谢(GO：0009694)和生物合成。其中，MYB3R被评为干旱环境中最重要的基因。最近，越来越多的研究报道MYB与植物的耐旱性有关。例如，烟草中的R2R3型MYB基因NbPhan被认为在解决干旱胁迫方面发挥了作用，其烟草沉默导致烟叶严重萎缩和增加水分流失率[29]。MYB对纤维发育也至关重要。玉米中miRNA和MYB转录因子之间的相关调节机制有望用于双重目的，即纤维发育和耐旱性。对于盐度反应，包括颗粒蛋白重复半胱氨酸蛋白酶，醇脱氢酶1，富含甘氨酸的RNA结合蛋白3和赤霉素2-β-双加氧酶7的基因被确定为在盐度响应的重要性方面排名较高。有许多报道显示实际上颗粒蛋白重复半胱氨酸蛋白酶与干旱反应无关，而与植物中的程序性细胞死亡无关[30-31]。干旱胁迫可能更容易通过半胱氨酸蛋白酶触发植物中的程序性细胞死亡。

4 结 论

鉴定了337前体miRNA，其中包含289种已知的miRNA和48种新的miRNA。在3个文库，两个分组中共有155种差异表达miRNA。目标预测，基于GO功能分类和基于遗传和基因组(KEGG)的功能富集显示，miRNA可能通过靶向一系列与胁迫相关的基因而在干旱胁迫中发挥作用。至少55个预测的靶基因进一步被60个miRNA调控。NAC、MYB和MAPK基因家族在干旱胁迫下评分最高，在植物抗旱中重要作用。根据目标基因预测，一系列玉米miRNA与这些排名靠前的基因相关，包括miR164、miR172、miR1520、miR6158、ghr-n24和ghr-n56等。miRNAs可能在玉米雄穗花器官分化期耐旱中发挥重要作用，miRNAs的筛选为分子辅助育种和转基因育种将提供新的靶标。