蓖麻miR396基因家族及其靶基因GRF生物信息学分析及鉴定

邢海彦，唐杰松，于耸，，郑志民*

（1.东北林业大学生命科学学院，东北 盐碱植被恢复与重建教育部重点实验室，黑龙 江哈尔滨， 150040；2.东北林业大学林学院，林木 遗传育种国家重点实验室，黑龙 江哈尔滨， 150040）

MicroRNAs（miRNAs）是一类真核生物体内普 遍存在的内源性非编码小RNA，一般由21～22 个核苷酸组成。植物中的miR396高度保守，很多植物的miR396家族已经被鉴定，如拟南芥，水稻，番茄，大豆，杨树，烟草，小麦等[1～7]。miRNA 在转录后水平通过剪切或翻译抑制对其靶基因进行调控。植物中miR396主要靶基因是GRF转录因子家族，GRF是一类生长调节因子。GRF的N端由两个高度保守的结构域组成——QLQ 和WRC[8,9]。其中miR396 可识别GRF 的mRNA 上与之互补的序列，更准确的说是编码WRC 的最后一个氨基酸[10～12]。WRC 具有核定位信号，其特征是具有Trp-ArCys 基序和Cys 和His 残基的锌指状间隔（CX9CX10CX2H）[8～14]。另一个保守的QLQ 结构域含有芳香疏水性氨基酸结构，与酵母的SWI/SNF 的N端部分相似，介导了GRF和GIF（GRF-Interacting Factor）的结合[15,16]。番茄降解组分析发现，8个番茄GRF成员可被miR396剪切[17]。

先前的研究表明，miR396及其靶基因GRF参与调节细胞增殖[16]，叶片大小[18]，叶形态发生[19]，细胞分裂[20]，花器官发育[21]，种子大小[22]，根生长[23]等。据报道，在拟南芥中过表达miR396a不仅会导致花瓣和雄蕊数目均减少，表型严重时雌蕊由两心皮变为单心皮，而且会造成拟南芥角果变小[24]。在拟南芥中超表达miR396b，叶面积明显变小，细胞数目减少。当GRF上miR396的靶位点发生突变时（rGRF突变,通常设计为同义突变），miR396将不能识别rGRF，拟南芥叶面积明显增大，细胞数目増多，细胞面积基本不变[16]。大籽粒籼稻与小籽粒粳稻进行杂交后，Duan 等发现一个非常重要的QTL--GS2，编码转录因子OsGRF4，籼稻中的GS2在OsmiR396g/h/i的靶位点处突变了2个碱基,碱基由TC变成了AA，OsmiR396g/h/i对GS2AA的抑制作用解除。在粳稻ZH11超表达GS2AA后，OsGRF4表达上调，ZH11籽粒明显变大[25]。

miR396及其靶基因GRF不仅在植物生长发育过程中发挥重要作用，它还参与响应植物的生物胁迫和非生物胁迫[1,3,26～29]。干旱胁迫下，过表达miR396可提高拟南芥和烟草的存活率[1,3,30～31]。当受到紫外线UV-B 的照射时，miR396表达被诱导，而其靶基因AtGRF1、AtGRF2和AtGRF3表达下调，叶面积减小。过表达miR396和特定的grf突变可导致水稻[32]、烟草[3]和拟南芥对病原体更加敏感，而rGRF/GRF的过表达则增强了拟南芥和水稻对真菌和细菌病原体的抗性[31,33]。

蓖麻（Ricinus communisL.），大戟科蓖麻属单属种植物，双子叶，一年生草本或多年生木本植物，其种子含油量可达45%～54%[34]。蓖麻油为世界十大油料和四大不可食用的油料之一，酸值低、粘度大、比重高（0.958～0.968 g/cm3）、燃点高（322℃以上）、凝固点低（-18℃以下），是具有双键的十八碳羟基脂肪酸的唯一来源，享有“生物石油”的美誉。蓖麻具有耐旱、耐脊薄、耐盐碱等突出特点，因而适宜作为生物柴油原料[35]。但目前关于蓖麻中miR396 家族和其靶基因GRF家族的研究还未见报道。本研究中，我们在蓖麻基因组中鉴定了4 个miR396家族成员和11个GRF家族成员，并分析了GRF的系统发育关系，染色体分布和进化关系。同时利用RNA-seq数据，对GRF基因家族成员的表达模式进行了分析，为后续深入研究蓖麻中miR396-GRF 的生物学功能奠定了基础。

1 材料与方法

1.1 蓖麻miR396的预测

蓖麻全基因组数据由本实验室测序组装完成。从miRBase 数据库（http://www. mirbase. org/ftp. shtml/）下载所有miRNA 成熟序列，筛选出所有植物中miR396成熟序列。将蓖麻基因组与所有植物中miR396成熟序列进行blastp（https://blast.ncbi.nlm.nih. gov/Blast. cgi）比对，筛选出错配碱基数小于等于3 的相似序列[36,37]。将获取的相似序列剔除重复位置，通过TBtools 工具（https://github. com/CJChen/TBtools）获取相似序列上下游各200 bp的一段序列，作为候选的RcmiR396前体序列。根据miRNA二级结构特征，利用UNAFold（http://unafold.rna.albany. edu/? q=mfold/RNA-Folding-Form）在 线 工具，对候选RcmiR396前体序列进行二级结构的预测，去除无典型发卡结构的序列，从而获得蓖麻miR396家族成员。

1.2 蓖麻GRF基因家族成员的预测

从PlantTFDB 数据库（http://planttfdb. cbi. pku.edu.cn/）下载拟南芥的GRF转录因子蛋白序列。将蓖麻注释蛋白序列与拟南芥GRF 蛋白序列进行blastp（https://blast. ncbi. nlm. nih. gov/Blast. cgi）比对，获取E-value＜e-5相似序列[38]。使用Pfam 在线分析工具（http://pfam.xfam.org/）进一步验证，去除不完整GRF 结构域序列，从而获得蓖麻GRF 转录因子家族成员。使用ExPASy（https://web. expasy. org/compute_pi/）在线工具预测RcGRF 蛋白的等电点（pI）和相对分子质量（MW）。

1.3 RcmiR396与靶基因RcGRF的作用位点分析

将RcGRF的基因序列和RcmiR396 的成熟序列输入到在线靶位点分析软件psRNATarget（http://plantgrn.noble.org/psRNATarget/）中，分析两者的作用位点。

1.4 蓖麻GRF转录因子家族系统进化分析

为了进一步探究蓖麻GRF 转录因子家族成员在基因组内和在基因组间的进化关系，从PlantTFDB 数据库下载拟南芥、大豆、水稻、毛果杨和木薯的GRF 转录因子蛋白序列，构建系统发育树。用MEGA 10.0 自带的ClustalW 工具分别对包括蓖麻在内的84 个GRF 蛋白序列（拟南芥9 个，大豆21 个，木薯17 个，毛果杨14 个，水稻12 个）进行多序列比对，采用了基于bootstrap 抽样的邻接法（neighbor—joining，NJ），重复值设定为1000。并且通过MEGA 10.0 对蓖麻GRF 蛋白序列进行多序列比对，构建进化树，辅助分析基因结构和保守基序、顺式作用元件。

1.5 蓖麻GRF转录因子家族的结构与保守基序分析

从蓖麻基因注释信息GFF 文件中获取蓖麻GRF 的CDS 序列序列，使用TBtools 工具（https://github.com/CJ-Chen/TBtools）将其转换为蛋白序列，通过NCBI CDD-search 在线工具对GRF家族成员的保守结构域进行预测，获取其保守结构域和位置信息。同时将蓖麻GRF 蛋白序列上传至MEME-5.0.4 在线分析工具（http://meme-suite. org/），预测其保守基序，参数为默认值，最后将二者合并，通过TBtools工具对其进行可视化展示[39]。

1.6 蓖麻GRF转录因子家族启动子顺式作用元件分析

利用TBtools 工具，从蓖麻基因组注释文件中提取蓖麻GRF家族基因起始密码子ATG 上游2 kb 作为潜在的启动子区域，将上述启动子区域序列提交至PlantCARE 工具（http://bioinformatics. psb. ugent.be/webtools/plantcare/html/），预测该区域存在的顺式作用元件，筛选保留其中的生长调控、激素和非生物胁迫响应元件（包括：ACE、LTR、TCA-element、ABRE、ARE、AuxRR-core、RY-element、CGTCA-motif、TGACG-motif、O2-site、GCN4_motif、TATA-box、GARE-motif、P-box、AACA_motif、MBS、MBSI、MRE）进行分析和展示。

1.7 蓖麻GRF家族基因的表达模式分析

采用FPKM（fragments per kilobases per millionfragments）衡量基因的表达水平，筛选出RcGRF 基因，通过TBtools实现基因表达层次聚类。

2 结果与分析

2.1 蓖麻miR396家族的预测

通过多序列比对，获得了105 个非冗余的蓖麻miR396 成熟序列候选序列。通过对105 个候选序列的二级结构预测，得到了4 个RcmiR396 家族成员，分别命名为RcmiR396a、b、c、d，如图1 所示，4 个RcmiR396 均具有经典的发卡结构。如表1 所示，RcmiR396序列较为保守。

表1 RcmiR396基因家族的成熟序列及其在染色体上的位置Table 1 Mature sequence of RcmiR396 gene family and their location in chromosome

图1 RcmiR396家族成员的二级结构预测Fig.1 Secondary structure predictions of RcmiR396 family

2.2 蓖麻miR396靶基因GRF家族成员的鉴定和命名

通过blastp 序列比对与Pfam 保守结构域验证，共获得11 个蓖麻GRF家族基因（表2），根据他们在染色体上的位置，分别命名为RcGRF1—RcGRF11。这11个GFR基因均包含完整的WRC和QLQ保守结构域，其中RcGRF5上存在两个WRC 结构域，他们的蛋白质长度在318～619 aa 之间，理论等电点在5.54～9.34 之间，分子量最小为35 179.57 Da，最大值则为67 639.97 Da，由此可见，蓖麻GRF在基因结构相对保守的情况下，理化性质又各不相同，这可能暗示着蓖麻GRF不尽相同的功能。

表2 鉴定的11个RcGRF家族成员Table 2 List of 11 RcGRF genes identified

2.3 蓖麻miR396 与其靶基因RcGRF 作用位点的分析

根据在线工具预测，结果如表3 所示，RcmiR396a,b 靶向相同的6 个RcGRF成员（RcGRF2、3、4、6、8、11），RcmiR396d 则可以靶向其余5 个GRF成员，RcmiR396c 仅靶向RcGRF3。同时，蓖麻GRF 家族成员的靶点大多位于外显子区和3′UTR区，仅Rc⁃GRF3的靶点在内含子区，这也与在水稻中的研究一致。

表3 RcmiR396的靶基因Table 3 Target genes of RcmiR396

2.4 蓖麻GRF 转录因子家族的系统发育分析和分类

通过构建来自五种植物GRF 蛋白序列的系统进化树，分析GRF 基因家族间进化关系。结果如图2 所示，84 个基因聚类为5 个亚族，蓖麻GRF 在五个亚族中均有分布，拟南芥中的GRF 分布在第II、III、IV、V 亚族中，水稻中的GRF 分布在第I、II、III、V 类中，并且IV 亚族中GRF 成员均来自于双子叶植物，表明亚族II、III、V 亚族在单子叶双子叶分化之前已经存在，在单子叶分化之后产生了IV 亚族，而Rc-GRF 在IV 亚族分布最多，这也说明了蓖麻在进化上与同属双子叶植物的大豆、木薯、毛果杨进化关系更为紧密。如图3 所示，蓖麻GRF 转录因子每个亚家族成员间WRC 和QLQ 结构域氨基酸序列相对保守。不仅如此，RcGRF2 和3 在C 端包含FFD、TQL和GGPL 结构域，RcGRF6、8、10 则分别包含了部分的结构域。

图2 RcGRF基因家族的系统发育分析Fig.2 Phylogenetic analysis of RcGRF gene family

图3 RcGRF家族成员保守结构域WRC和QLQ的氨基酸序列比对Fig.3 Alignment of multiple RcGRF WRC and QLQ domain amino acid sequences

2.5 蓖麻GRF转录因子家族的基因结构分析

为进一步探究RcGRF基因结构特征，本研究中阐明了该家族基因的内含子数量、外显子个数和保守基序。如图4 所示，11 个RcGRF基因中，均含有内含子，RcGRF1，5，7，9均包含2 个内含子，其余RcGRF成员包含3 个内含子，RcGRF成员外显子个数也在4～5 个之间。不同基因之间开放阅读框（ORF）中内含子-外显子结构差异不大。本研究中使用MEME 在线分析工具共鉴定出RcGRF 转录因子中的10 个保守基序，其中，Motif1 是WRC 结结构域，Motif2是QLQ结构域。

图4 蓖麻GRF家族的进化树，保守基序和基因结构Fig.4 Phylogenetic tree,motif distribution and gene structure of RcGRF gene family

2.6 RcGRF家族启动子顺式作用元件的预测

对RcGRFs启动子序列的分析结果显示，多个RcGRF启动子区富集了GA 响应、ABA 响应、IAA 响应、茉莉酸甲酯响应、水杨酸响应、干旱胁迫、厌氧诱导响应、胚乳表达调控、低温响应、分生组织表达、光响应等顺式作用元件（图5）。RcGRF启动子序列共包含了4种激素响应元件，5种非生物胁迫响应元件和3 种生长发育相关调控元件。所有的Rc⁃GRF启动子序列均包含厌氧诱导响应元件和GA 响应元件，某些基因可能含有多个某一顺式作用元件，RcGRF4和10均含有8 个响应茉莉酸甲酯的元件，RcGRF2和5均含有4 个响应ABA 的元件。这表明RcGRF的功能可能涉及调控生长发育过程、激素信号转导与环境胁迫响应。

图5 蓖麻GRF家族成员启动子顺式作用元件的预测Fig.5 Prediction of cis-regulatory elements of the GRFs in Ricinus communis

2.7 蓖麻GRF的组织表达差异分析

GRF 转录因子家族成员广泛参与调控植物对根、茎、叶、花等组织的发育。为探究该家族基因与蓖麻生长发育的关系，本研究根据实验室已有的转录组测序数据，系统分析了RcGRF基因在根、茎、叶、花中的表达模式。结果如图6 所示，除RcGRF4外，其余蓖麻GRF 成员在蓖麻的各个组织中均有表达，RcGRF4在蓖麻各个组织中未检测到表达，可能是具有时空表达特异性。蓖麻叶片和茎、花中表达量最高的为RcGRF9，种子中表达量最高的为Rc⁃GRF6，根中表达量最高的为RcGRF2。聚类分析表明，RcGRF成员在花和茎中的表达模式相似且表达量相对较高，叶和根、种子中的表达模式更为一致表达量更低，这可能暗示着蓖麻GRF 在调控株高和发育方面的作用。同时，RcGRF5、RcGRF7和Rc⁃GRF9在进化上与OsGRF1、3、4、6、10更为紧密，而水稻这些GRF成员已被报道对开花时间和花器官发育及茎伸长有明显影响，蓖麻GRF成员在花和茎中的表达明显高于其他组织中，与在水稻中研究一致，这可能暗示了蓖麻GRF在花和茎发育中起到的作用。

图6 RcGRF基因家族表达模式图Fig.6 Expression pattern of RcGRF gene family

3 讨论

蓖麻种子含油量可达45%～54%，我国蓖麻产量位居世界前列，然而对蓖麻分子生物学的研究较少。主要原因是受限于蓖麻基因组信息的不完整及缺乏高效稳定的遗传转化体系，实验室在前期的研究中组装了高质量染色体级别的蓖麻基因组，为蓖麻生物学研究和功能基因的挖掘奠定了基础。在水稻和拟南芥中对miR396-GRF的研究已十分清晰，miR396-GRF 参与调控植物的生长发育过程包括叶形态发生、花器官发育、株高、籽粒大小，所以对于蓖麻中miR396及其靶基因GRF家族成员的研究是十分重要且又有现实意义的。

本研究利用实验室测序组装的蓖麻基因组数据鉴定了4 个miR396和11 个GRF家族成员，并对miR396的二级结构和GRF家族成员特征进行了全面分析。结果发现，4 个miR396 前体序列均可以形成发卡结构，且与已发现的miR396 成熟序列相比，错配碱基数在1～3 个，RcmiR396a、b也仅有一个碱基的差别，这表明了miR396的保守性。RcmiR396与RcGRFs 作用位点的分析显示RcmiR396a、b、d可以靶向不同的的RcGRF家族成员，这暗示了不同的miR396家族成员在蓖麻中可能发挥不同的作用。从组织表达分析结果来看，RcmiR396a、b、d靶向了多个RcGRF成员，而这些RcGRF在不同组织中的表达并不一致。如RcmiR396d的靶基因中组织表达模式可分为两类，RcGRF1、5、7、9在花和茎中均表达量最高，而RcGRF10仅在种子中有较高表达，这说明了RcGRF的表达多样性。拟南芥中进化分析表明，6个物种共84个GRF 聚类为5个亚族，蓖麻GRF在五个亚族中均有分布。在RcGRF 分布最多的IV亚族中，大豆、木薯、毛果杨GRF 成员数目占比明显增多，这也说明了蓖麻在进化上与同属双子叶植物的大豆、木薯、毛果杨进化关系更为紧密。在部分RcGRF的C 端包含TQL、FFD 和GGPL 结构域，这在拟南芥[40]、水稻[8]、卷心菜[41]和玉米[42]中都有发现。

RcmiR396通过调控RcGRF的表达量来影响植物的生长发育，因此对蓖麻GRF 家族成员表达模式的分析有助于研究其功能。在水稻中，大部分GRF成员在多个组织中表达，并参与水稻生长过程中多种生理过程。对RcGRF启动子顺式作用元件分析表明，RcGRF对激素诱导、胁迫响应、生长发育均有响应，这也与前人在水稻和拟南芥中的研究一致。系统发育分析表明，RcGRF1、7在进化上与OsGRF6较近，水稻中的研究证明OsGRF6在调控株高和茎伸长中发挥着重要作用[43]。在RcGRF1、6、7启动子区存在分生组织表达调控元件，组织表达分析表明这些RcGRF在茎组织中的表达量较高，这可能暗示了它们在茎生长中起着重要调节作用。在蓖麻中，GRF4在各个组织中均无表达，这暗示GRF的表达可能有时空特异性。除RcGRF4、6、10外，其余GRF成员在花中均表达量很高，这暗示了RcGRF成员可能调控了蓖麻从营养生长到生殖生长的转变。