紫花苜蓿细胞质雄性不育系相关microRNA的鉴定

郭 强，王英哲，陈晶晶，周 仂，徐 博*

(1. 吉林农业大学动物科学技术学院，吉林 长春 130118； 2. 吉林省农业科学院，吉林 长春 130033)

紫花苜蓿(MedicagosativaL.)是一种多年生豆科苜蓿属植物。1958年，加拿大学者发现紫花苜蓿的雄性不育现象，细胞质雄性不育(Cytoplasmic male sterility，CMS)指开花植物不能产生功能性花药、花粉或雄配子的一种遗传学现象[1]，这种现象主要由细胞质因子和细胞器基因控制[2]，其不育特性表现为雄器官的退化、畸形或功能缺失，导致花药、花粉母细胞和雄配子的生长发育产生异常[3-4]。常用于作物的杂交育种。

随着雄性不育系的发现，国外的众多学者开始采用分子技术、组织培养法研究遗传机理，并利用不育系培育杂交组合和转基因材料[5-7]；1978年，吴永敷等[8]在我国内蒙地区从‘草原1号’中发现了6株紫花苜蓿雄性不育系；2008年，于洪柱等[9]发现了紫花苜蓿细胞质雄性不育材料—MS-GN。随着紫花苜蓿细胞质雄性不育系及保持系的发现，对其机理的探索也在不断深入。同时，国内学者[10-14]对不育系展开了同工酶分析、细胞学观察、生殖生物学、生理生化机制等内容的研究；利用紫花苜蓿细胞质雄性不育系选育杂交组合并对其F1代进行生物学性状及抗逆性研究[15-18]；有学者[19-21]利用SSR标记以及BSR技术对紫花苜蓿细胞质雄性不育恢复基因进行了定位分析，并且对紫花苜蓿细胞质雄性不育系及其保持系的花药发育过程进行切片观察，测定其在不同发育阶段生理生化指标的变化。

植物microRNA是一类长度为20～26 nt的内源性非编码小RNA[22]。成熟的microRNA长度为20～24 nt，多为21 nt[23]，其中长度为21～24 nt的microRNA占主体地位。植物microRNA参与植物体内许多重要的生命过程，拟南芥(ArabidopsisthalianaL.)体内的miR160可通过负调控生长素响应因子进一步调控根冠细胞的发育[24]；Zea-miR390则可以调控玉米(ZeamayL.)茎尖细胞的发育及叶片的形态建成[25]；miR172，miR159，miR156/157则可以调控开花时期[26-27]。microRNA也参与调控植物的环境适应性，如miR398与植物的耐逆性密切相关[28]。另外，microRNA参与植物激素合成代谢与信号传导，如miR160，miR164，miR167，miR528等microRNA均与生长素的信号传递有关，其中miR167还可以参与脱落酸等其他的信号通路的传导[29-30]。

基于全基因组信息可知，在许多雄性不育植物的花粉发育过程中均有大量的microRNA参与[31-33]。在拟南芥中，miR165/166参与小孢子的发育，miR156和miR159负责调控植物绒毡层细胞的发育过程，miR159的靶基因属于GAMYB(Gibberellins acid v-myb avian myeloblastosis viral oncogene homolog)家族，主要调控植物花药细胞程序性死亡(Programmed Cell Death,PCD)过程中的赤霉素的传导，miR159的异常表达会使绒毡层无法降解，最终使小孢子不能正常发育而丧失育性[34-36]。拟南芥中若缺乏miR165/166会使花药只产生两个花粉囊，最终会导致拟南芥的育性降低[37]；水稻(OryzasativaL.)中的OsmiR528，OsmiR159，OsmiR172d；玉米中Zea-miR397；油菜(BrassicanapusL.)中的bra-miR167，bra-miR172等miRNA都可能与植物的小孢子发育有关[31,38]。另外，有研究发现，某些miRNA可通过参与植物激素合成和应答过程来调控花粉的发育，例如miR167缺失会引起茉莉酸的合成异常，从而导致花药不开裂以及散粉受阻；miR167的过表达同样会引起植物不育现象的产生[30,39-41]。

植物microRNA是揭示植物雄性不育机制的有效方式之一，基于之前的研究发现：紫花苜蓿细胞质雄性不育的败育可能与线粒体基因表达异常、淀粉蔗糖含量亏缺、绒毡层发育异常等现象有关[17,42]，因此，为深入探索不育机制，本试验利用small RNA-seq技术，找到雄性不育相关的miRNA，为更好地利用紫花苜蓿不育系进行杂交育种工作提供理论基础，以期选育出更具优良性状的种质资源。

1 材料与方法

1.1 试验材料

本试验所选用的材料为紫花苜蓿细胞质雄性不育系MSGN1A及其对应的同型保持系MSGN1B，MSGN1A和MSGN1B是以具有雄性不育性的‘公农3号’为供体母本株系，通过保持系筛选，回交4代后选育出来的。所有植物材料均由吉林省农业科学院提供。MSGN1A和MSGN1B的个体植株于2016年4月在吉林省长春市吉林农业大学草业科学试验田种植。选取现蕾期24 h的花药部分，用于提取总RNA，备用。试验设置3个生物学重复。

1.2 试验方法

1.2.1总RNA的提取及质量控制 参照Garg等[43]的方法提取MSGN1A和MSGN1B花药总RNA，用紫外分光光度计分别在230，260和280 nm处测定吸光值(OD230，OD260和OD280)，OD260/OD280和OD260/OD230均须≥1.8。采用Agilent 2100生物分析仪检测总RNA的完整性，总RNA完整性值在8.0～10.0的样品可用于后续试验。

1.2.2small RNA文库构建 将提取到的总RNA送北京百迈克生物技术公司，完成small RNA文库的构建及测序工作。提取两个样品总RNA，使用聚丙烯酰胺凝胶电泳法根据片段的长度分离RNA并回收长度在18～30 nt之间的条带，构建2个small RNA文库。

1.2.3miRNA的鉴定与分析 对已知miRNA的鉴定，采用miRBase(v21)软件将比对上参考基因组的reads序列与known miRNA数据库中的成熟miRNA序列进行序列比对。若序列与known miRNA完全一致，则这个reads将被认为是鉴定到的known miRNA。

考虑到miRNA的生物学特征，对于没有被鉴定到known miRNA的序列，使用miRDeep 2软件对novel miRNA做预测。利用miRDeep 2软件包，将reads比对到基因组上的位置信息，得到可能的前体序列，根据reads在前体序列上的分布信息及前体结构能量信息，使用贝叶斯模型经打分最终实现novel miRNA的预测。

1.2.4差异表达miRNA筛选及其靶基因功能注释 以|log2(FC)|≥1；FDR≤0.01作为标准筛选差异表达的microRNA。差异倍数(Fold Change，FC)表示两样品间表达量的比值。

依据known miRNA和novel miRNA与对应物种的基因序列信息，使用Target Finder软件[40]对miRNA靶基因进行预测。

基因功能注释由北京百迈客生物技术公司进行，差异表达miRNA靶基因的功能注释所使用的数据库为GO (Gene Ontology，http://www.geneontology.org/)和KEGG (Kyoto Encyclopedia of Genes and Genomes，http://www.genome.jp/kegg/)。

2 结果与分析

2.1 测序数据分析

对测序得到的数据进行质量控制，每个样品的Clean Data均大于14.91 M。细胞质雄性不育系MSGN1A及其同型保持系MSGN1B的Raw reads分别为27.83 M和28.28 M，Clean reads分别为20.24 M和14.91 M，Q30(质量值大于30)分别为98.87%和99.03%。

表1 测序数据产出统计表

2.2 miRNA的鉴定

为了进一步鉴定已知的miRNA，将比对上参考基因组的reads序列与已知miRNA数据库miRBase(v21)中的成熟miRNA序列进行比对。序列与已知miRNA完全相同的reads被认为是鉴定到的已知miRNA。针对miRNA的生物特征，对于未鉴定到known miRNA的序列，利用miRDeep 2软件进行novel miRNA的预测。通过鉴定，本试验共检测到1 478个miRNA，其中已知miRNA 1 351个，novel miRNA 127个(表2)。

表2 miRNA统计结果

试验检测到1 351个known miRNA，分属于136个miRNA家族，部分结果如图1所示。每个家族中所含miRNA的数量和表达水平不尽相同，其中包括属于miR166家族的247个miRNA，miR167家族的159个miRNA，miR160的家族124个的miRNA，miR172家族的105个miRNA，mi156家族的98个miRNA，miR319家族的88个miRNA等。成熟miRNA长度主要集中在20 nt到24 nt的范围内，其中植物的miRNA以21 nt或24 nt为主，试验鉴定出的novel miRNA的长度主要为21 nt，其次是24 nt和22 nt。

图1 部分已知miRNA在miRNA家族中所含数量

2.3 差异表达miRNA的筛选与鉴定

在差异表达miRNA检测过程中，使用|log2(FC)|≥1；FDR≤0.01作为筛选标准。试验共计筛选出差异表达miRNA 130个，包括90个上调表达的miRNA和40个下调表达的miRNA。试验共鉴定到40个差异表达的已知miRNA，分属于15个miRNA家族，与保持系MSGN1B相比，紫花苜蓿细胞质雄性不育系MSGN1A中有24个上调表达miRNA，有16个下调表达的差异miRNA。

表3 差异表达的已知miRNA

2.4 miRNA靶基因预测

根据known miRNA和novel miRNA与其对应物种的基因序列信息，使用Target Finder软件进行靶基因的预测，1 351个known miRNA中，有1 347个miRNA预测到靶基因，靶基因数目为9 848个；127个novel miRNA中，有126个miRNA预测到8 252个靶基因。注释到的靶基因数目为16 458个。

表4 miRNA靶基因数目预测统计

2.5 差异表达miRNA靶基因注释

试验共得到16 458个靶基因，其中10 637个靶基因获得功能注释信息，紫花苜蓿细胞质雄性不育系MSGN1A与其同型保持系MSGN1B间注释的差异表达miRNA靶基因在各注释库中的数目统计见,表5，共有10 637个miRNA靶基因被各个注释库注释到，其中COG注释库中有3 104个miRNA靶基因被注释到，包括815个差异表达miRNA靶基因；GO注释库中有5 147个miRNA靶基因被注释到，包括1 264个差异表达miRNA靶基因；KEGG注释库中有3 479个miRNA靶基因被注释到，包括861个差异表达miRNA靶基因；KOG注释库中有5 620个miRNA靶基因被注释到，包括1 286个差异表达miRNA靶基因；Pfam注释库中有8 309个miRNA靶基因被注释到，包括1 871个差异表达miRNA靶基因；Swissprot注释库中有6 535个miRNA靶基因被注释到，包括1 662个差异表达miRNA靶基因；eggNOG注释库中有9 133个miRNA靶基因被注释到，包括2 261个差异表达miRNA靶基因；Nr注释库中有6 535个miRNA靶基因被注释到，包括2 612个差异表达miRNA靶基因。

表5 miRNA靶基因注释统计表

2.5.1GO注释分类 对差异表达miRNA靶基因进行GO分类，共计1 264个miRNA靶基因得到注释。按照GO的三个大类进行下一次级分类，共分成41个功能群，其中包括18个生物学过程，12个细胞成分，11个分子功能(图2)。

图2 差异表达miRNA靶基因GO注释分类

在生物学过程中，“代谢过程”所占比例最大，其次是“细胞过程”、“信号生物过程”、“生物调控”、“刺激应答”、“定位”、“细胞成分组织或生物合成”等生物学过程，它们的占比均大于10%；“发育过程”、“多细胞生物过程”、“繁殖过程”、“信号”、“多生物过程”、“生长”、“繁殖”、“免疫系统过程”等生物学过程的占比均在1%～10%之间；“生物附着”、“生物相”、“节律过程”的占比在0.1%～1.0%之间。

在细胞组分类别中，“细胞部分”和“细胞”占比最大，紧接着是“细胞器”、“膜”、“细胞器部分”、“高分子配合物”，其占比均大于10%；“胞外区”和“细胞连接”的占比在1%～10%之间；“膜封闭腔”、“病毒”以及“病毒部分”的占比在0.1%～1.0%之间。

在分子功能类别中，“合成”和“催化活性”是最大的功能类别；其次是“转运蛋白活性”、“核酸合成转录因子活性”、“电子传递体活性”、“结构分子活性”、“酶调节活性”以及“分子传感器活性”，占比均在占比在1%～10%之间；最后是“抗氧化剂活性”、“受体活性”、“蛋白合成转录因子活性”和“鸟苷酸交换因子活性”，占比均不到1%。

图3 差异表达miRNA靶基因KEGG通路注释

2.5.2KEGG通路分析 对差异表达miRNA靶基因进行KEGG通路分析，差异表达miRNA靶基因所涉及的通路的类别主要包括：“细胞进程”、“环境信息过程”、“基因信息过程”、“代谢”、“生物系统”等五个大类，在此基础上进一步分为82个类别，其中有24个靶基因注释到“代谢”中的“氨基酸生物合成”通路中，占比为5.61%；18个靶基因注释为“内吞作用”通路中，占比4.21%；16个靶基因注释到“核糖体”通路中，占比为3.74%；“RNA转运”和“碳代谢”通路中分别有15个靶基因被注释到，占比均为3.50%；“内质网中的蛋白质过程”、“泛素介导的蛋白水解作用”、“植物激素信号传导”通路中分别有14个靶基因被注释到，占比均为3.27%；有13个靶基因被注释到“剪接体”通路中，占比为3.04%；“氨基糖和核糖代谢”和“氧化磷酸化”通路中有12个靶基因被注释到，占比均为2.8%；“淀粉蔗糖代谢”通路中有8个靶基因被注释到，占比为1.87%。

2.6与花药发育相关的miRNA

在差异表达的miRNA中，mtr-miR5256的靶基因的功能注释为“glucan ando-1,3-beta-glucosidase 4”,该基因参与淀粉蔗糖代谢。miR159靶基因的功能注释为“MYB transcription factor”，“MYB转录因子”是花分生组织识别基因的启动子。miR172靶基因功能注释为“AP2 Domain transcription factor”，AP2类似转录因子主要调控花器分化和花期。miR396靶基因功能注释为“生长调控因子HSP70蛋白”，miR169靶基因功能注释为“CBF HAP2 like protein”，其主要作用是调控花期，这些差异表达的miRNA均有可能导致花粉败育。

3 讨论

MiR159靶基因的功能注释为“MYB transcription factor”，“MYB转录因子”是花分生组织识别基因的启动子。MYB是谷物糊粉蛋白细胞中赤霉素(Gibberellic acid，GA)信号的组成部分，在花的发育过程中起着重要作用[44]。陶园[45]认为，拟南芥转录因子ICE1(Inducer of CBF Expression 1)通过调控GAMYB基因来调节赤霉素信号传导，从而影响花粉育性并调控植物生长。miR159的过量表达会引起与花发育过程相关的转录因子MYB的表达异常，从而进一步影响到花粉绒毡层的正常代谢，包括绒毡层的形态和结构异常等[46]。在本研究中miR159b和miRNA159c均下调表达，miRNA159d上调表达，这些miRNA的表达异常让我们有理由认为miR159家族与紫花苜蓿细胞质雄性不育有着密切联系。

试验所得结果中miR172，miR172a，miR172d相较于保持系来说，在不育系中均上调表达，miR172靶基因功能注释为“AP2 Domain transcription factor”，AP2-like转录因子主要调控花器分化和花期。AP2-like转录因子对花器官的生成是至关重要的，如果花器官发育出现异常情况(如花丝和花粉囊发育异常)就很可能会使花粉粒无法正常发育，从而导致花粉败育[47]。Chen等[48]研究认为，miR172过表达会导致植株的花期提前以及花器官发育畸形，其表型与ap2突变体表型类似，最终导致配子异常发育。因此，miR172在不育系中的表达量降低也可能导致花粉败育，具体的调控模式还需进一步研究。

miR396靶基因功能注释为“生长调控因子HSP70蛋白”。试验结果显示miR396家族的miR396，miR396a，miR396b-5p，miR396c-5p，miR396d，miR396e等miRNA在不育系中上调表达，而miR396在初级和成熟的miRNA中都具有高度保守性，是“拟南芥生长调节因子”家族成员[49-50]。Yang等[51]认为miRNA396介导的NtGRF-like基因表达调控通路与烟草叶片和花的发育过程密不可分，这些miRNA对应的靶基因调控花的发育过程，也就是说，miR396的异常表达可能影响花的发育。

由结果可知，与保持系相比，miR169c在不育系中上调表达，miR169靶基因功能注释为“CBF HAP2 like protein”，其主要作用是调控花期，拟南芥miR169家族成员是发育过程和刺激通路中普遍存在的调控因子[52]。因此，miR169的表达量降低会可能导致花期的提前或延迟，影响花药的发育等过程。

4 结论

差异表达microRNA靶基因功能注释结果显示，与紫花苜蓿细胞质雄性不育相关的差异表达microRNA家族主要有4个，包括miR159，miR172，miR396，miR169；其靶基因的功能主要是调控花药的发育，其次是调控花期以及信号传导；micro RNA靶基因主要涉及的通路有：“代谢过程”、“信号生物过程”、“植物激素信号传导”、“氧化磷酸化”和“淀粉蔗糖代谢”等，而差异表达microRNA靶向的转录因子中多为控制花期、花药发育和信号传导等过程的因子，对不育系的败育起到一定的影响。本试验结果可初步揭示紫花苜蓿不育系的分子机理，为利用紫花苜蓿不育系进行杂交育种提供理论基础，从而选育出具有优良性状的紫花苜蓿种质资源。