青花菜细胞质雄性不育相关miRNA的鉴定与表达特征分析

裴徐梨，荆赞革,2，唐 征，冯鹏宇

(1 昆明学院 农学与生命科学学院，昆明 650214；2 温州科技职业学院 农业与生物技术学院，浙江温州 325006)

miRNA是一类重要的内源性单链小分子RNA[1]。前人研究表明，miRNA及其靶基因在花粉发育和育性调控中具有重要的作用[2]。例如，miR156可通过调节SPL转录因子的表达来影响拟南芥雄配子的产生[3]；miR159a-1能够抑制MYB基因的转录水平来调控花芽分化和植株育性[4]；miR172与其靶基因(AP2类基因)共同调节拟南芥的开花时间和器官形成[5]。近年来，利用高通量测序技术很多与植物育性相关的miRNA已被成功鉴定出来[6-8]。

青花菜杂种优势非常明显。细胞质雄性不育系以其育性稳定、开花能力优良和杂交率高等优点已被广泛应用于青花菜杂交育种体系中。本研究利用Illumina测序技术构建了青花菜细胞质雄性不育系和保持系的sRNA文库，共鉴定到47个差异表达miRNA。同时，通过qRT-PCR技术检测了部分miRNA及其靶基因的表达特征，本研究结果为进一步探讨miRNA调控青花菜育性提供一定的理论基础。

1 材料和方法

1.1 试验材料

以回交8代以上的青花菜Ogu细胞质雄性不育系WN12-95A(MS)及其保持系WN12-95B(FS)为试验材料。于花期采集长度小于2 mm的不育系和保持系花蕾用于高通量测序。采集不育系和保持系3个不同时期花蕾(长度< 2 mm、2～4 mm、> 4 mm)以及保持系大于4 mm花蕾的不同部位(花萼、花瓣、雄蕊和雌蕊)用于qRT-PCR分析。

1.2 试验方法

1.2.1 青花菜sRNA文库构建和miRNA的鉴定从质量合格的总RNA中分离和纯化得到长度合适的small RNA进行文库构建，而后进行高通量测序。采用SOAP2将过滤得到高质量clean read比对到青花菜转录组参考序列[8]。已知miRNA的鉴定通过miRDeep2软件比对到miRBase数据库；新型miRNA的预测则通过碱基延伸未比对上的序列，再通过MIREAP软件进一步的筛选。

1.2.2 青花菜不育系和保持系差异miRNA的筛选将miRNA表达水平进行归一化处理，使用归一化的数据计算Fold Change(FC)和P值。差异表达miRNA的筛选以|log2(FC)| ≥ 1，P< 0.05 作为标准。

1.2.3 青花菜差异miRNA的靶基因预测及功能注释使用psRNA Target软件进行潜在靶基因的预测[9]，错配小于4的基因被选为候选靶基因。同时基于 sRNA 和 mRNA 测序数据，采用关联分析的方法鉴定雄性不育相关差异表达基因[10]与差异miRNA的关联网络。

1.2.4 实时荧光定量分析差异miRNA及靶基因的表达特性miRNA及总RNA的提取和反转录均使用TaKaRa公司试剂盒。实时荧光定量参照Xu等的方法[11]，使用SYBR Premix Ex Taq试剂盒(TaKaRa公司)在ABI 7500(Applied Biosystems公司)上运行。miRNA反应程序为：95 ℃ 10 s，95 ℃ 5 s，60 ℃ 35 s，40个循环。靶基因反应程序为：95 ℃ 30 s；95 ℃ 15 s、60 ℃ 1 min，40个循环。基因相对表达量采用2-ΔΔCt法来计算。以U6 和actin基因分别作为miRNA及其靶基因的内参基因。试验所用引物序列见表1。

表1 qRT-PCR引物信息

2 结果与分析

2.1 青花菜不育系和保持系sRNA数据分析

构建了青花菜雄性不育系和保持系花蕾的2个sRNA文库，共得到47.88 M原始数据(表2)。过滤掉经接头序列和低质量序列后，分别获得15.44 M(不育系)和24.26 M(保持系)的clean reads，分别占Reads总数的76.6%和87.5%。过滤后序列的总碱基对不育系和保持系分别为351 525 505 bp和535 664 137 bp。

2.2 青花菜不育系和保持系miRNA的鉴定

青花菜不育系和保持系中共鉴定到69个家族中的181个已知miRNA，包括138个保守miRNA和43个非保守miRNA。其中miR156家族成员数最多，为10个，而大部分保守miRNA家族仅发现1个成员。miR403和miR167表达丰度分别为351 989和 278 824，在2个样品中呈现高表达；而miR1520、miR529和miR535的表达丰度均较低。

青花菜不育系和保持系中共鉴定到8个新型miRNA(表3)，其中1个新型miRNA在保持系中特异性表达。对这些新型miRNA的基本信息进行分析，结果显示它们的最小自由能(the minimum free energy，MFE)为-40.70～-24.10 kcal/mol，长度介于21到25 nt。

表3 青花菜不育系和保持系中鉴定的新型miRNA

2.3 青花菜雄性不育相关miRNA的鉴定

通过比较2个文库中的miRNA的表达量，结果(图1)发现, 43个已知miRNA和4个新型miRNA在青花菜不育系和保持系中差异表达。其中，24个已知miRNAs和4个新型miRNAs表达上调，其余miRNA表达下调；6个差异表达miRNA的差异倍数大于8；2个差异miRNA只在不育系中表达，7个差异miRNA为保持系所特有。

A.已知miRNA； B.新型miRNA图1 青花菜不育系(MS)和保持系(FS)差异miRNA表达量比较分析Fig.1 Comparative relative expressions of differentially expressed miRNAs in male sterile line (MS) and maintainer line (FS) of broccoli

2.4 青花菜差异表达miRNA的靶基因预测

利用psRNA Target软件共预测到188条靶序列。通过进一步的靶基因功能分析，发现其中一些靶基因是雄性不育发生和花粉发育过程中的转录因子。例如AP2-like转录因子、锌指结构域蛋白基因、激素响应因子(auxin response factors家族，ARF)和SPL(squamosa promoter binding protein-like)结合蛋白基因等。同时，一些靶基因编码的磷酸甘油酸激酶、三角状五肽重复蛋白质、油质蛋白、谷氨酸脱羧酶和γ-微管蛋白复合体等都是该过程中重要的生物蛋白。

图2显示，差异表达miRNA的靶基因经GO注释为17个生物学过程类目，10个细胞成分类目和11个分子功能类目，表明这些靶基因参与了多种代谢和发育过程。KEGG注释结果(图3)将靶基因分为23个代谢途径，其中“植物病原菌的相互作用”(ko04626)、“N-糖链的合成”(ko00510)和“线粒体的脂肪酸延伸”(ko00062)为靶基因参与最多的途径。COG注释结果(图4)可将靶基因分为19个不同的功能类。其中属于“一般功能预测”类群的靶基因数最多(图4)。

图2 青花菜差异表达miRNA靶基因的GO注释Fig.2 Gene ontology analysis for target genes of differentially expressed miRNAs in broccoli

图3 青花菜差异表达miRNA靶基因的KEGG注释Fig.3 Kyoto Encyclopedia of Genes and Genomes analysis for target genes of differentially expressed miRNAs in broccoli

图4 青花菜差异表达miRNA靶基因的COG注释Fig.4 Clusters of orthologous groups of proteins analysis for target genes of differentially expressed miRNAs in broccoli

2.5 青花菜雄性不育相关miRNAs与mRNAs的关联分析

将本实验室的青花菜sRNA 和 mRNA数据进行关联分析，鉴定到2个关联的miRNA-mRNA对(miR859-1/T2-Unigene_BMK.21608和miR5174e-1/T2_Unigene_BMK.31772)(图5，A)。分析发现2个miRNA均负向调控相应的靶基因(图5，B)。对miRNA的靶基因进行注释发现，miR859-1靶向的花粉外壁蛋白基因(T2/Unigene_BMK.21608)和miR5174e-1调控的氨基酸通透酶基因(T2_Uni-gene_BMK.31772)均可参与雄性不育的发生过程。

图5 青花菜雄性不育相关miRNAs与mRNAs的关联分析(A)及表达特性(B)Fig.5 The association analysis (A) and expression patterns (B) of CMS-related miRNAs and mRNAs in broccoli

2.6 青花菜育性相关miRNA表达特性分析

选取青花菜育性相关的8个miRNA(5个已知miRNA和3个新型miRNA)进行荧光定量PCR分析，在花蕾不同时期的表达模式分析结果(图6，Ⅰ)显示，miR159a-1、miR164a-1、miR319a-1和miRn11在不育系(A1)中的表达量高于保持系(B1)。其中，miR164a-1、miR319a-1和miRn11的表达水平随着不育系花蕾的发育逐渐降低，但在保持系花蕾中的表达量随着花蕾发育呈现上升趋势；miR159a-1的表达量在不育系和保持系花蕾发育初期均最高，以后表达量逐渐下降。此外，miR397a-2、miR168a-4、miRn4和miRn16在保持系中上调表达，但表达模式各有差异。miRn4在不育系花蕾中的表达量逐渐增加，在保持系中表达趋势相反；其余3个miRNA在不育系和保持系中的表达趋势一致。同时，进一步比较靶基因和miRNA的表达特性后，结果发现miRNA与其对应的靶基因基本呈现相反的表达趋势。

A. 不育系;B. 保持系;1、2、3分别代表 < 2 mm、2～4 mm和 > 4 mm的花蕾图6 青花菜miRNA及其对应靶基因在不同发育时期花蕾和不同部位中的表达特征A. Male sterile line; B. Maintainer line; 1, 2 and 3. were buds with length of < 2 mm, 2-4 mm, and > 4 mm, respectivelyFig.6 Expression characteristics of miRNAs and their corresponding targets in different developmental stages and different parts of flower buds in broccoli

由图6，Ⅱ可以看出，miR168a-4、miR159a-1、miR164a-1和miRn4在青花菜不同部位均有表达,且表达量差异不大；miR319a-1在雄蕊中的表达量最高，在花萼、花瓣和柱头中的表达量降低；miRn11在花萼中的表达量明显高于其他组织，而与之相反，miRn16在花萼中的表达量极低；miR397a-2在雄蕊中表达量最高，但在花瓣和柱头中几乎不表达，推测其具有组织特异性。

3 讨 论

miRNA的基本特征分析是了解其调控作用的基础。本试验中sRNA的长度分析结果符合植物miRNA的特征。然而，青花菜的23 nt miRNAs的数量高于21 nt miRNAs，表明这些miRNA调控育性的机制可能具有其特点。近年来，大量研究表明已知miRNA通常具有较高的表达水平，且其功能在植物进化中较保守。而新型miRNA的表达量明显低于保守miRNA。在本研究中鉴定到的青花菜miRNA也存在此现象。

3.1 青花菜表达上调miRNA的调控机制探讨

植物雄性不育的发生是一系列相关miRNA参与的结果。本研究共鉴定到29个上调表达的miRNA，其中miR397a-2靶向laccase-like转录因子。该转录因子能够调控植物发育过程中油菜素内酯的信号转导。油菜素内酯作为雄性不育发生的重要因子，其含量升高可降低拟南芥花粉的释放效率[12]。此外，上调表达水稻OsLAC(laccase-like gene)基因可导致部分雄性不育[13]。本试验中miR397a-2在青花菜不育系中的表达量显著低于保持系，因此推测miR397a-2可能是通过调控LAC基因的表达来参与青花菜的雄性不育发生。

激素响应因子(ARF类基因)在大豆中可以通过结合激素响应启动子调控花粉发育。miR395的靶基因中就有ARF3(auxinresponsefactor3)，且该miRNA在保持系中的表达量高出不育系4.7倍。因此，推测miR395能够通过调节激素介导的信号转导来影响青花菜的育性。

3.2 青花菜表达下调miRNA的调控机制探讨

miRNA与其靶基因相互作用来调控花粉发育[14]。miR156通过介导SBP-box基因来调控花粉发育。在拟南芥中，过表达SPL可以保护花粉育性[3]。青花菜miR156在不育系中高表达，推测miR156可能通过调控SBP-box基因的表达来影响其育性。三角状五肽重复蛋白(pentatricopeptide repeat gene,PPR)在RNA编辑中具有重要功能。miR158、miR400、miR859和miR5067均靶向PPR基因，且在青花菜保持系中的表达量低于不育系。现有研究表明PPR基因可使萝卜育性恢复[15-18]。因此，我们推测miRNA抑制雄性不育相关的PPR基因的表达是青花菜雄性不育发生的可能原因之一。

Ca2+运输功能障碍可影响花粉发育。钙/钙调依赖性蛋白激酶是钙离子流入细胞的调节器[19]。预测靶向钙/钙调依赖性蛋白激酶基因的miR399在青花菜保持系中下调表达，推测其在不育系中的高表达引起花粉发育过程中钙稳态的失调和Ca2+介导的信号通路的减弱，从而导致雄性不育的发生。