茶树miR164a及其靶基因的鉴定与表达分析

孔雷，朱向向，王屹玮，谢小芳，江昌俊，李叶云

安徽农业大学茶树生物学与资源利用国家重点实验室，安徽 合肥，230036

MicroRNAs（miRNAs）是内源性单链非编码小 RNA，在动物和植物中主要是在促进靶向 mRNA的剪切或翻译抑制中起着重要的调节作用[1]。miR164家族是植物中广泛存在的一类保守microRNA，近年来，关于miR164家族成员参与植物抗逆胁迫的报道越来越多[2]。miR164的表达具有组织及发育阶段特异性，在植物细胞生长、组织发育、器官形成中起着重要的调控作用[3-4]。在拟南芥中，miR164在根部和花中的表达量比在叶片中的高[5]。在拟南芥和烟草花器官、分生组织以及原形成层中有miR164表达，且表达量较高[6]。在逆境胁迫条件下，植物miR164响应逆境胁迫，并且可能作为调控分子在防御系统中发挥着重要作用[3]。在非生物胁迫条件下，机械损伤会抑制毛果杨miR164的表达[7]，盐胁迫和干旱胁迫处理会抑制甜杨幼苗 miR164的表达[8]。

植物 miR164调控的靶基因属于植物NAC（NAM，ATAF1/2 和 CUC2）转录因子[3]。NAC转录因子是一个植物特有的转录因子家族[9-10]。总的来说，NAC蛋白的N端有由150个保守氨基酸组成的NAC结构域，C端有一个多样化的转录调控区域[11]。近年来，NAC蛋白因其在植物生长发育和对环境适应的过程中发挥着多重作用而被深入研究[12]。NAC转录因子可以由转录水平上的某些顺式作用元件和反式作用因子、转录后水平的 miRNA以及包括磷酸化、蛋白降解和二聚化的翻译后调节水平调控[13]。许多NAC蛋白在非生物胁迫反应中发挥作用[14]。拟南芥ANAC019、055和072响应干旱、盐和ABA等非生物胁迫[15]，拟南芥NAC转录因子LOV1响应冷胁迫[16]。过表达水稻SNAC1显著提高转基因水稻的抗旱性[17]。

研究表明，拟南芥 miR164家族（ath-miR164a、b和c）调控 5个 NAC转录因子基因（CUC1，CUC2，NACl，at5g07680，和at5g61430）的mRNA的降解是植物器官边界建成、侧根的产生、营养和花器官的形成以及与年龄相关的细胞死亡所需的[5,18-19]。小麦miR164和其靶基因 NAC（miR164-NAC）对干旱胁迫的响应呈负相关[20]。胡杨peu-miR164通过负调控PeNAC070来提高对干旱胁迫和盐胁迫的耐受力[21]。

近年来，在茶树抗逆分子机理研究中，茶树中2个NAC转录因子获得克隆，并响应温度胁迫[22]，茶树 miR164a也报道响应低温胁迫[23-24]。但是茶树miR164a与NAC转录因子家族中哪一个成员（miR164-NAC）存在调控关系以及如何响应逆境胁迫尚不清楚。本研究根据低温胁迫下茶树小RNA高通量测序结果获得的茶树miR164a，进行前体序列的克隆与分析；通过茶树降解组测序获得 miR164a靶基因，利用 RLM-5´RACE加以验证，并对靶基因进行克隆与序列分析。探讨茶树miR164a及其靶基因在茶树不同叶位、高温胁迫和低温胁迫中的表达模式，为揭示茶树miR164a-NAC在茶树生长发育、抗逆性中的作用提供一定的理论依据。

1 材料与方法

1.1 试验材料与处理

以盆栽1年生无性系茶树品种舒茶早（C.sinensiscv.Shuchazao）为试验材料，温度胁迫试验在人工气候室中进行[25]。选取芽、第1叶、第3叶、第5叶、第7叶，高温胁迫处理温度为白天38℃，晚上28℃，处理0、1、2、3、4 d（以处理0 d为对照）；低温胁迫处理温度为 4℃，处理 0、3、6、12、24 h（以处理0 h为对照），在各个时间点采集相同位置的叶片混样，并迅速放入液氮中，储存在-80℃冰箱备用[25]。

1.2 植物总RNA的提取

试验材料经液氮研磨后，使用 miRcute miRNA提取分离试剂盒（TIANGEN，北京）提取茶树总RNA。用核酸定量仪测定总RNA样品的浓度、纯度，用 1.2%琼脂糖凝胶电泳检测提取的总RNA样品完整性，储存在-80℃冰箱备用[25]。总 RNA反转录为 cDNA利用PrimeScriptTM 1st Strand cDNA Synthesis Kit（TaKaRa，大连）试剂盒。

1.3 miR164a前体序列的克隆及其二级结构预测

使用植物MiniBEST Plant Genomic DNA Extraction Kit提取分离试剂盒（TaKaRa，大连）提取茶树基因组 DNA，根据从本实验室茶树低温miRNA高通量测序数据中得到的茶树miR164a的前体基因序列：TCTAGATTCGT CGAACAGTAGGTGGGCAGCTCCTTGTTGG AGAAGCAGGGCACGTGCAAAAAGCTAAT GCTAGTAAGTTTTGTGTACGTGCTCCCCT TCTCCAACACGGGTTCCCTCACTCCCTCG ATTGCGAGCTC，用Primer 5软件设计其特异性引物（表1）进行PCR扩增。利用RNAfold web server在线软件（http://rna. tbi. univie.ac.at/cgi-bin/RNAWeb Suite/RNAfold.cgi）对 miR164a的前体序列进行二级结构的预测[26]。

1.4 利用降解组测序数据预测靶基因

降解组测序和预测工作由联川生物公司（LC Bio Tech，杭州）完成，步骤如下：首先将测序获得的原始数据通过一系列数据处理得到可用于后续分析的可比对测序序列，然后将可比对序列与测序物种茶树的 cDNA数据库序列比对生成降解组密度文件（Degradome density file），通过Targetfinder靶基因预测软件预测出与测序物种茶树小RNA序列配对的靶基因mRNA序列，最后将预测的miRNA对应的靶基因和生成降解组密度文件中的 mRNA进行结合运算，找出共同具有的mRNA，该mRNA即为miRNA的靶基因，并给出降解组的峰值分类和分值，并对产生的预测结果进行作图（t-plots）[2]。

1.5 RLM-5´RACE验证miR164a靶基因裂解位点

从降解组测序数据中获知 gi319875685（NAC21/22）为miR164a的靶基因，根据从茶树转录组中获得的 gi319875685（NAC21/22）的基因序列片段，用 EditSeq和 MegAlign软件查询 ORF框，根据其 ORF两端的序列设计RT-PCR引物（表1），测序验证该序列的准确性。将获得的cDNA序列在NCBI中进行BLASTn比对分析，结果显示与公开的茶树CsNAC1基因序列（GenBank:JN857066.1）相一致。

依据CsNAC1的 cDNA序列用 Primer Premier 5.0软件设计用于RLM-5´RACE的引物（表1），RLM-5´RACE 试验参照 5´-Full RACE试剂盒（TaKaRa，大连）说明书进行，将扩增出的特异产物进行克隆测序[2]。

1.6 靶基因CsNAC1序列生物信息学分析

用 EditSeq分析氨基酸序列。从 NCBI（https://blast.ncbi.nlm. nih.gov/Blast.cgi）得到保守结构域，在 NCBI数据库中分别用BLASTn和 BLASTx进行核酸序列及氨基酸序列比对分析，用 MAFFT version 7（ https://mafft.cbrc.jp/alignment/server/index.html）比对蛋白质的同源性序列，拟南芥、水稻 NAC家族转录因子基因序列的获得以及NAC结构域和 miR164结合区域的确认参照文献[21,27]，将序列用ClustalX比对后，利用MEGA 6.06软件并通过邻接（Neighborjoining）算法构建系统进化树[25]。

1.7 茶树miR164a及其靶基因在茶树中的表达

根据从本实验室茶树低温miRNA高通量测序数据中获得的茶树 miR164a成熟序列设计具有茎环结构的反转录引物，应用PrimerSelect软件设计 miRNA上游引物[28]。利用 PrimeScriptTMRT reagent kit（TaKaRa，大连）进行反转录反应。荧光定量实验参照SYBR Premix Ex TaqTMII试剂盒（TaKaRa，大连）说明书进行。pc-3p-222和CsEF-1α分别作为茶树 miR164a qRT-PCR和靶基因CsNAC1qRT-PCR的内参基因[29-30]。采用CFX96TMReal-Time PCR Detection System扩增仪检测，每个样品3次生物学重复。验证目标基因和内参基因的扩增效率和扩增范围，采用法进行相对表达量的计算[25]。

表1 引物序列Table 1 Primer sequences

2 结果与分析

2.1 Csn-miR164a前体序列的克隆及其二级结构的预测

从本实验室茶树低温 miRNA高通量测序数据中得到茶树miR164a的前体序列，以茶树总DNA为模板进行PCR扩增验证，结果表明两者序列完全一致，其长度为 126 bp（图1）。Csn-miR164a前体序列的GC含量是53.97%，最小折叠自由能（Minimal folding free energies，MFE）为-65.60 kcal·mol-1，最小折叠自由能指数（Minimal folding free energies index，MEFi）为 0.96，MEFi大于 0.85，符合miRNA前体的特点[31]。RNAfold分析结果显示，Csn-miR164a前体序列可折叠成较为典型的二级结构，其成熟片段位于前体结构的5´端茎结构上（图2），其成熟序列为：5´-UGGAGAAGCAGGGCACGUGCA-3´，21 bp，与miRbase 21数据库中拟南芥、水稻和烟草等22个物种的miR164a成熟序列相一致。

2.2 茶树Csn-miR164a靶基因的预测

将预测的Csn-miRNA164a对应的靶基因和生成降解组密度文件中的 mRNA进行结合运算，找出共同具有的mRNA（Gi319875685）。对产生的预测结果进行作图（Target plot，t-plot），Alignment Score（Target finder的打分）为2、Alignment Range（配对作用位点）为631～651、Cleaveage Site（剪切位点）为 642、Category（分型）为 0、P-value（P值）为0.003 95、Rep_Norm_reads（匹配多位点迚行位点归一化后的reads）为18.667。

在t-plot图中，miRNA剪切靶基因mRNA的位点用箭头标出，降解组测序预测的miR164a剪切位点位于靶基因mRNA的第642位核苷酸上（图3）。

2.3 RLM-5´RACE验证茶树miR164a靶基因

利用RLM-5´RACE对降解组测序预测的茶树 miR164a靶基因进行验证，扩增出单一目的条带（图4），大小约为 231 bp，与预期的一致。为了准确验证靶基因的裂解位点，挑取10个阳性克隆单菌落测序，测序结果全部显示茶树miR164a对靶基因的mRNA进行单一剪切，与 miR164a结合的区域为5´-AGCAAGUGCCCUGCUUCUCCA-3´ ， 剪切位点位于miRNA/target mRNA结合区域的第10位核苷酸上，第10位和第11位碱基之间，与降解组测序结果一致（图5）。

图1 Csn-miR164a前体序列的克隆电泳图Fig. 1 Electrophoresis results of Csn-MIR164a cloning

图2 Csn-miR164a前体二级结构预测Fig. 2 Predicted secondary stem-loop structure of Csn-MIR164a

图3 利用降解组测序技术获得茶树miR164a靶基因的t-plot图Fig. 3 Target plots (t-plots) of miR164a target confirmed by degradome sequencing in Camellia sinensis

2.4 Csn-miR164a靶基因全长ORF序列的验证及其序列分析

从茶树转录组中获取个i319875685 mRNA序列，NCBI BLASTx比对显示它具有完整的开放阅读框（ORF）。设计引物扩增验证其全长 ORF，PCR结果显示扩增出的条带大小与预期目的片段大小基本一致，约1 000 bp（图6）。测序后，BLASTn分析表明，该序列与NCBI登录的一个茶树NAC1基因（登录号：JN857066.1）的同源性达到99%。而与文献报道的两个茶树NAC蛋白的同源性性分别只有30.26%和29.23%。该序列ORF长度为912 bp,编码 303个氨基酸。保守结构域预测如图7所示，CsNAC1中的第 33位至第 109位氨基酸序列为 NAM 超家族的保守结构域，这是NAC主要的特征功能域。

在NCBI数据库中对CsNAC1氨基酸序列进行 BLASTx比对，结果显示它与木豆（XP_020223961.1）的NAC21/22氨基酸序列相似性最高，为 73%，与毛果杨（XP_002310688.1）、野生大豆（KHN26800.1）的NAC相似性分别为72%、71%。

用MAFFT version 7对拟南芥、水稻和茶树3种植物的14条NAC家族转录因子进行蛋白质序列分析，如图8所示，这14条NAC家族转录因子的N端均有NAC结构域的5个亚结构域A—E，C端含有miR164的结合区域[21,27]。

选取20种不同植物中的NAC21/22构建进化树分析，结果显示，茶树CsNAC1与藜麦（Chenopodium quinoa）NAC21/22的亲缘关系最近（图9）。

图4 RLM-5´RACE扩增产物Fig. 4 PCR Production of RLM-5´RACE

图5 miR164a靶基因裂解位点的验证Fig. 5 Identification of cleavage sites of miR164a target

图6 茶树CsNAC1 ORF扩增产物电泳图Fig. 6 The PCR product of CsNAC1 ORF in tea plant

图7 CsNAC1的保守结构域Fig. 7 Conserved domain of CsNAC1 protein

图8 茶树、拟南芥和水稻miR164靶基因NAC的多序列比对Fig. 8 Multi-sequence alignments of the miR164 NAC target genes in tea plant, rice and Arabidopsis thaliana

图9 茶树CsNAC1的进化树分析Fig. 9 Phylogenetic tree analysis of CsNAC1

2.4 茶树miR164a和其靶基因CsNAC1表达模式分析

2.4.1 茶树Csn-miR164a和CsNAC1在不同叶位中的表达分析

通过荧光定量PCR，检测了Csn-miR164a及其靶基因CsNAC1在茶树不同叶位的表达情况（图10）。结果显示，Csn-miR164a和CsNAC1在茶树不同叶位的表达存在差异，Csn-miR164a在芽中的表达量最高，而靶基因CsNAC1在芽中的表达量最低。随着叶片的发育，Csn-miR164a随着叶位的降低表达量依次降低，在第7叶位（老叶）的表达量最低，仅为芽的47%。CsNAC1的表达量随叶位的降低表达量逐渐增高，在第7叶位表达量最高，是芽的8倍。可见，Csn-miR164a随着叶片的衰老表达量下降，而CsNAC1与Csn-miR164a呈现负相关，随着叶片的衰老表达量增加。

2.4.2Csn-miR164a和CsNAC1对高温胁迫的响应

从图11可以看出，高温诱导茶树miR164a及其靶基因CsNAC1的表达。在胁迫处理的第1天高温诱导miR164a的表达显著下调，表达量为对照（0 d）的 58%，靶基因CsNAC1的表达量显著上调，相对表达量是对照的3.5倍。miR164a表达量都在高温胁迫处理的第4天达到最低，仅为对照的 46%，靶基因CsNAC1在第4天的表达量最高，是对照的6.6倍。在高温胁迫处理下茶树 miR164a的表达量随胁迫处理时间的增加呈下调趋势，而靶基因CsNAC1的表达量呈上调趋势。茶树在高温胁迫处理下，靶基因CsNAC1的表达量与Csn-miR164a的表达量呈负相关。

2.4.3Csn-miR64a和CsNAC1对低温胁迫的响应

如图12所示，低温诱导茶树miR164a及其靶基因CsNAC1的表达，在低温处理3 h后，Csn-miR164a的表达量显著下调，仅为对照（0 h）的42%，CsNAC1上调表达，相对表达量是对照的5.7倍。在低温处理24 h后，miR164a的表达量达到最低，相对表达量为对照的41%，而CsNAC1的表达量达到最高，相对表达量是对照的 20.3倍。在低温胁迫处理下，靶基因CsNAC1的表达量与Csn-miR164a的表达量呈负相关。

3 讨论

研究表明许多高度保守的miRNA响应逆境胁迫并有助于植物适应环境，miRNA靶基因大部分是转录因子，如AP2、NF-Y、MYB和NAC。植物miR164及其靶基因NAC基因在植物中是保守的，NAC蛋白代表一个植物特有的广泛存在的转录因子家族，其调控植物的发育和应激反应[21]。在拟南芥中，miR164家族靶向7个NAC基因的转录物；在水稻中，miR164家族与6个NAC基因互补；在玉米中，发现7个NAC基因是预测的miR164的靶基因[21]。在茶树中，初步发现NAC转录因子家族有45个成员[32]，然而本试验仅发现与验证一个 NAC转录因子为茶树 miR164a的靶基因，是否还存在其它家族成员是 miR164a的靶基因，需要进一步研究。

图10 Csn-miR164a和其靶基因CsNAC1在茶树不同叶位中的表达分析Fig. 10 Expression analysis of Csn-miR164a and target gene CsNAC1 in different leaf position of tea plant

图11 茶树Csn-miR164a及其靶基因CsNAC1在高温胁迫下的表达分析Fig. 11 Expression analysis of Csn-miR164a and target gene CsNAC1 in tea plant under heat stress

图12 茶树Csn-miR164a及其靶基因CsNAC1在低温胁迫下的表达分析Fig. 12 Expression analysis of Csn-miR164a and target gene CsNAC1 in tea plant under cold stress

在水稻中，RLM-RACE分析显示miR164定向剪切NAC转录因子家族成员OMTN1和OMTN2的剪切位点在第 10和 11位碱基之间[27]。在胡杨中，5´-RACE分析显示miR164定向剪切PeNAC070、PeNAC012和PeNAC028mRNA，剪切位点在miR164第10和11位碱基之间[21]。在光皮桦中，运用5´-RACE验证了NAC1为miR164靶基因，剪切位点也在第10和11位碱基之间[33]。在本研究中，我们通过降解组测序数据和RLM-5´-RACE验证，获知茶树miR164a定向剪切CsNAC1，剪切位点同样也在 miR164a第 10和 11位碱基之间。这表明miR164剪切位点可能是非常保守，没有物种差异。

miR164表达的时序性和组织特异性可能决定了植物组织和细胞的功能特异性，对细胞生长和组织发育的调节起重要作用[3]。在拟南芥中，miR164随着拟南芥叶片的衰老表达量下降[19]，甜菊新叶中miR164的表达量高于老叶[34]。在本研究中，随着茶树叶片的衰老，miR164表达量也呈下降趋势，与其他植物一致，且Csn-miR164a与CsNAC1呈负相关关系，推测Csn-miR164a参与茶树叶片生长发育过程，但是起何种作用尚不清楚。

植物miR164通过调控靶基因NAC转录因子参与植物抗逆胁迫[2]。对毛竹进行高温（42℃）和低温（4℃）处理，发现 miR164b的表达均呈现不同程度的下调，其靶基因PeNAC1上调表达[35]。对大豆进行低温（4℃）处理，miR164也是下调表达[36]。与毛竹、大豆一致，茶树Csn-miR164a低温温度下调表达，靶基因CsNAC1上调表达，呈负相关关系。但是 0℃低温处理下甜杨幼苗、4℃低温处理下小麦中 miR164a的表达均上调[37-38]，这表明不同物种miR164对低温胁迫响应模式存在差异。

本研究初步揭示了Csn-miR164a-CsNAC1参与茶树叶片的生长发育过程和响应温度胁迫。但是，茶树miR164a功能与作用仍需通过过表达转基因等手段进一步的研究，miR164a-NAC1通过何种调控途径来响应温度胁迫，提高茶树抗逆性的分子机制也有待深入研究。