粗梗水蕨营养期和生殖期微RNA的测定与分析

孙 军 李贵生

（植物资源保护与利用湖南省高校重点实验室，吉首大学，生物资源与环境科学学院，吉首 416000）

拟南芥孢子体有营养期到生殖期的转变，而它们受到相互拮抗的miR156 和miR172 的控制。拟南芥有3个MIR156基因，它们产生的miR156能通过碱基互补配对的方式结合SQUAMOSA-PRO⁃MOTER BINDING PROTEIN-LIKE（SPL）基 因 的mRNA，并引起后者的断裂或者翻译抑制［1-2］。miR156 的浓度在幼期植株中很高，而在发育过程中逐渐下降［3-5］。反过来，SPL转录本的水平前期低、后期高［6］。MIR156的组成型表达延长植株幼期，而其功能失去的植株没有幼期的特征，譬如短的节间和光滑的叶下表面［7］。拟南芥SPL3/4/5的miR156结合位点被屏蔽后使得成熟和开花提前［8］，而加强SPL9的功能可以使幼期被跳过［4，9］。相应地，SPL9/15的双突变可打断从幼态到成熟的转变［4，9］。拟南芥的SPL2/10/11/13也在从幼态到开花的转变中发挥作用［3］。正因为如此，增加植株的糖分可以通过抑制MIR156的转录而加快生长时期的转变［10］。

拟南芥有5 个MIR172基因，它们编码的miR172 在营养期浓度较低、但在生殖期会积累到较高水平，而其APETALA2（AP2）-like 靶基因的转录本、特别是其蛋白产物的浓度会反方向变化［11-12］。特别是，MIR172D在顶端分生组织积累，而MIR172A/B积聚在叶的维管束，从而分别在短日照下通过年龄途径和在长日照下通过光周期途径去加速开花过程［13-14］。miR172水平过高会使拟南芥提前开花［11-12］。反过来，AP2-like靶基因的超表达抑制开花［11，15］，而它们的功能失去突变导致过早开花［11，16-17］。MIR172基因的转录能被SPL蛋白激活［4］，而miR172的积累水平在spl15的突变体里大大降低［18］。

蕨类植物的孢子体也有营养生长和生殖生长，而它们是否也由miR156 和miR172 调控不得而知。粗梗水蕨（Ceratopteris pteridoides）是一种水生的水龙骨蕨类［19-21］，在配子体世代首先发育成丝状体、然后发育成铲状雄性原叶体或者连指套状两性原叶体。粗梗水蕨孢子体的叶在营养期为通常的片状，但在生殖期叶为条状、其下表面长有孢子［22］。我们测定以及分析了粗梗水蕨营养期和生殖期的小RNA（sRNA），发现miR156 和miR172 以及它们的靶基因有和开花植物相似的积累模式。因此，miRNA 控制生长时期的转变可能在二者的共同祖先即原始真叶植物里就已出现。

1 材料与方法

1.1 材料收集

营养期的粗梗水蕨生长在无菌的MS 培养基中，28 ℃、16 h 光照。生殖期的植株位于培养钵中，25 ℃、连续光照、70%湿度。整株采集，营养期植株采集为受精3 周以后、因而无任何孢子发生；生殖期植株通过孢子囊的存在而确认，并在采集时去掉衰老的根和额外的叶。样品先用液氮速冻，然后保存于-80 ℃，运输时用干冰。

1.2 建库测序

总RNA用CTAB法提取［23］，18～30 nt RNA被切胶回收、接着先后连上5′和3′接头、然后逆转录，最后通过PCR扩增构建文库、并用Agilent 2100 Bioana⁃lyzer（Agilent Technologies，Waldbronn，Germany）和ABI StepOnePlus Real-Time PCR System（Applied Biosystems，CA，USA）进行质量检测［24］。测序平台为Hiseq 2000（Illumina，SanDiego，USA），测序公司为华大集团（中国，深圳）。测序原理为边合成边测序，预期测序读长为50 nt。测序数据编号为SRA022988（http：//www.ncbi.nlm.nih.gov/Traces/sra）。

1.3 数据分析

滤掉测序得分低的序列，并进一步去掉含有接头、少于18 nt、以及为polyA的序列，然后鉴定全长sRNA。用bowtie 程序（http：//bowtie-bio.source⁃forge.net/index.shtml）查询miRBase 数据库（http：//www.mirbase.org/ftp.shtml）以鉴定保守的miRNA，其中要求错配少于两个、没有间断，并最后在RNAcentral 平台（https：//rnacentral.org/）逐个比对以及肉眼检查。同样查询GenBank（tp：//ftp.ncbi.nlm.nih.gov/genbank/）和Rfam（http：//rfam.janelia.org/）以 鉴 定rRNA/scRNA/snoRNA/snRNA/tRNA，同时根据3′端2 个碱基的突出鉴定siRNA。最后用mireap（https：//sourceforge.net/projects/mireap/）对剩下的sRNA 进行新miRNA 预测。用psRobot（http：//omicslab.genetics.ac.cn/psRobot/）和Target⁃Finder（https：//github. com/carringtonlab/TargetFind⁃er）进行miRNA靶基因预测，并取其并集作为最终结果。用DEGseq（http：//www.r-project.org/）进行表达差异分析，以qvalue<0.05和|log2FoldChange|>1为阈值。

2 结果与分析

2.1 sRNA测序

营养期sRNA 测序产生了22 596 914 个原始序列，其中干净序列有21 930 111 个、占比97.78%。生殖期测序产生了21 612 191 个干净序列，占所有原始序列的94.20%。在营养期的干净序 列 中，miRNA 占9.72%，rRNA/tRNA/snRNA/snoRNA 占47.95%，而剩下的42.31%未能被注释（见图1A）。生殖期情况类似。21 nt 的干净序列最多，接着是20、24 nt的序列（见图1B）。

图1 干净序列的统计A.营养期样品干净序列对应的各种sRNAs的数目；B.2个样品中干净序列的长度分布Fig.1 Statistics of clean reads A.Clean reads corresponding to various sRNAs with sums in the vegetative sample；B.Length distribution of clean reads in the two samples

2.2 保守miRNA的鉴定

通过序列相似性分析鉴定了42 个保守的miRNA（见表1）。它们多数覆盖了最佳匹配的90%以上，而只有包括miR394 和miR845 在内的7个miRNA 覆盖了最佳匹配的76%～86%。通过对RNAcentral 所有动植物miRNA 的默认参数比对，有14 个miRNA 能在比蕨类植物更早的苔藓植物和石松植物里找到同源物，而剩下的绝大部分只能在种子植物里找到同源物。有12 个miRNA 只能找到1 个相似性匹配，而有着明确功能的miR⁃NA 能在多个物种里找到匹配。所有miRNA 在营养期和生殖期之间差异表达，其中25 个miRNA 具有2 倍以上的差异。9 个miRNA 在营养期表达上调，而更多的在生殖期表达上调（P<0.01）。

表1 两个样品中鉴定到的保守miRNATable 1 Conserved miRNAs identified in the two samples

miR319a-3p、 miR394a、 miR397、 miR3630-3p、miR4414a-5p 和miR7806 只在生殖期表达，而miR408 和miR845b 只在营养期表达。最后，miR156、miR529 和miR535 有着同样的靶基因（如下所述），而它们的差异表达也一样，只是miR157有点例外。同样，miR165和miR166有着近似的差异性表达，而miR166 的变体之间，以及miR396 的变体之间也是如此。

2.3 靶基因及其表达

我们也对这两个时期进行了转录组测序，这使得我们能够分析以上miRNA 的靶基因及其表达。参照种子植物里miRNA和靶基因之间的对应关系，我们发现粗梗水蕨的miR408 有3 个cupre⁃doxin靶基因和2个uclacyanin靶基因（见表2）。这些基因在生殖期表达上调、从而与对应的miRNA反向水平表达。4 个SPL基因可被预测为miR156a、miR157d、miR529d 和miR535d 的 靶 位点，且它们都在生殖期表达上调从而与所有相互作用的miRNA 反向水平表达。2 个NUCLEAR FACTOR Y（NF-YA）靶基因在生殖期表达上调，而它们的5′端不翻译区有miR169a-5p 的结合位点、并与之反向水平表达。一个LEAF CURLING RE⁃SPONSIVENESS（LCR）靶基因在生殖期的表达下调可能与此时miR394a 的表达上调有关。有3 个AP2-like 基因可能被miR172 靶向，只是其中2 个基因与miRNA 反向水平表达的同时、另一个与之正向水平表达。3 个LTR 逆转座子的引物结合位点能被miR845 识别，并且二者之间表达水平反向。5 个class Ⅲhomeodomain leucine zipper（HDZIPIII）基因可被预测为miR166a-3p 和miR166m、但不是miR165a-3p 的靶位点，不过二者都在生殖期表达上调。miR171a-5p 和它的一个SCARE⁃CROW-LIKE（SCL）靶基因同时在生殖期表达上调。miR160h 和它的3 个AUXIN RESPONSE FACTOR（ARF）靶基因，miR159c-5p 和它的2 个MYB靶基因也都同时在生殖期积累。miR319b和miR159有同样的MYB靶基因，二者之间也是正向水平表达。miR397 和它的1 个cupredoxin 靶基因之间也无反向水平的表达。最后，在粗梗水蕨里没有找到miR319 的TCP靶 基 因，miR167 的ARF靶 基 因，miR396 的GROWTNH-REGULATING FACTOR（GRF）靶基因，以及miR168 的Argonaute1（AGO1）靶基因。剩下的miRNA大部分缺乏靶基因方面的可信研究，故这里没有对它们进行探讨。

表2 已知miRNAs靶基因的差异表达Table 2 Differential expression of target genes of known miRNAs

2.4 tasiR-ARF的产生

如前所述，miR390 在粗梗水蕨的生殖期有更高水平的积累，我们也从相应的转录组里找到了它的初级转录本（见图2A）。粗梗水蕨有2 个TAS3基因，一个是含有2 个串联tasiR-ARF 的长的CpTAS3a，一个是只有1 个tasiR-ARF 的短的Cp⁃TAS3b。有24 225 个干净序列及其对应的288 个sRNA 可以对应到CpTAS3a，而其中24 014 个干净序列及226 个对应的sRNA 位于它的2 个miR390互补位点之间（见图2B）。这些sRNA 的长度在20～25 nt，而其中173 个（～60%）为21 nt。进一步，4 个21 nt sRNA 可对应33 个干净序列、并且和3′端miR390 切割位点相位吻合；特别是其中1 个对应着10 个干净序列、并刚好位于5′D1（+）区域。最后，4 个sRNAs 与5′端miR390 切割位点相位吻合，并且其中一个刚好从切割位点起始。但是，没有任何sRNA 对应着CpTAS3a的2 个预测tasiRARF。不过，CpTAS3a的第1 个tasiR-ARF 与它的5′端miR390 切割位点相位吻合。在CpTAS3b中，有238 个sRNAs 位于2 个miR390 互补位点之间，8个sRNAs 与3′端的miR390 切割位点相位吻合、5个与5′端的miR390 切割位点相位吻合，并且有1个sRNA 位于5′D1（+）区域。最后，一个相位吻合的sRNA 在营养期有158 个对应的干净序列、在生殖期有154 个对应的干净序列，并且准确落在了CpTAS3b预测的tasiR-ARF位点。

图2 tasiR-ARF的鉴定A.MIR390 序列及miR390 和miR390*位点；B.TAS3 序列及对应sRNAs；Cp.Ceratopteris pteridoides；Cg.Cyrtomium guizhouense；Pn.Parathelyp⁃teris nipponica；和5′端miR390 切割位点相位吻合的sRNA 用前向箭头表示；和3′端切割位点相位吻合的sRNA 用后向箭头表示；CpTAS3b中对应于tasiR-ARF的sRNA用红色实线箭头表示；CpTAS3a中可能的tasiR-ARF的sRNA用红色虚线箭头表示Fig.2 Characterization of tasiR-ARF A.MIR390 sequences with sites of miR390 and miR390*；B.TAS3 sequences with corresponding sRNAs；Cp.Ceratopteris pteridoides；Cg.Cyrtomium guizhouense；Pn.Parathelypteris nipponica；sRNAs in phase with the 5′cleavage site of miR390 were indicated by forward arrows and those with the 3′cleavage site by backward arrows；One sRNA corresponding to tasiR-ARF was indicated by a filled red arrow in CpTAS3b，and a putative sRNA corresponding tasiR-ARF was shown by a punctuated red arrow in CpTAS3a

3 讨论

3.1 粗梗水蕨的miRNA

Illumina HiSeq 仪器能检测每个细胞中低至一个的分子，并且能同时检测miRNA、siRNA、piR⁃NA、rRNA、tRNA 以 及mRNA 的 降 解 产 物［25］。sRNA 被界定为小于200 nt，而其中miRNA 为20～24 nt。植物中经典的miRNA 为21 nt，它们被AGO1 蛋白复合体招募；长的miRNA 有24 nt，被AGO4 蛋 白 复 合 体 招 募［26-27］。siRNA 长 为24 nt，piRNA长为30 nt，而sRNA的峰值长度在动物中为20 nt。我们的测序参数与这些经验吻合（见图1），因此粗梗水蕨营养期到生殖期转变相关的miRNA可能就在所测的sRNA当中。

基于序列相似性在粗梗水蕨中发现了多个保守miRNA，但是其中一些并没有显示很高的整体相似度，这可能是数据库中没有与之最近缘的。粗梗水蕨的miRNA 多数也见于种子植物，这可能与蕨类植物是它的姊妹类群、并且二者同属于真叶植物有关［28-29］。粗梗水蕨的一些miRNA 只能在一个物种中找到匹配，所以它们可能是假阳性、从而没有被收入到某些数据库中［30-31］。其他蕨类植物研究没有发现miR397，这可能是因为它只存在于叶以外的组织；因此，它和miR845、miR894、miR394、miR172、miR168 以及smiR169 可能在陆地植物的进化中首次出现于蕨类植物［32-33］。

3.2 粗梗水蕨miRNA的功能

miR156a、miR157d、miR529d 和miR535d 在粗梗水蕨中都靶向4个不同的SPL基因，并且前者几乎都在生殖期表达下调而后者都在此时表达上调；这与它们在种子植物中的情况非常相似，因此它们可能也参与蕨类植物从营养期到生殖期的转变。鉴于生殖期表达下调的miRNA 显著较少（见表1），这一点更令人信服。在小立碗藓（Physcomi⁃trella patens）配子体的发育过程中，miR156 会逐渐积累而SPL靶基因受其影响同时表达下调，这表明了该相互作用在配子体中有相反的方向［34］。并且，在更早期的地钱（Marchantia polymorpha）中，没有miR156、而SPL受到物种特异性的miRNA 的控制［35］。miR172 不存在于苔藓植物和石松植物，并且它对AP2-like 基因的调控可能首先出现在蕨类植物［32，36-37］。粗梗水蕨中miR172 在生殖期表达上调，与此同时它的一些靶基因表达下调，因此可能与miR156-SPL 一样参与蕨类植物生长时期转变的调控。值得指出的是，AP2-like 靶基因对miR172 的反应也可能是翻译水平上的［11-12］。综上所述，蕨类植物的孢子体可能在进化中第1次同时使用了miR156 和miR172 去控制营养期到生殖期的转变。miR171 和miR159 都在粗梗水蕨生殖期表达上调，这和它们在种子植物中是miR156—SPL 相互作用的拮抗者和修饰者相一致［38-39］。如此看来，蕨类植物的营养期到生殖期的调控除了没有导致开花之外，与被子植物基本相同。

miR160 在粗梗水蕨的生殖期积累，而它在番茄中超表达会造成向下的内卷叶［40］，这与粗梗水蕨的可育叶非常相似。并且，它的一个靶基因在生殖期表达下调，而在拟南芥中它的靶基因的超表达会使叶向上卷［41］。粗梗水蕨中miR166 和miR319 在生殖期的表达上调也可能影响可育叶的卷曲［42］。miR394 在粗梗水蕨生殖期的表达上调可能调控叶的卷曲方向、也可能与胁迫应答有关［43］。这是因为，我们的营养期样品生长在无菌培养基中，而生殖期样品生长在开放的培养钵中。同样，miR408、miR169、miR397 和miR845 的差异性表达可能源自试验中铜、氮和水供应的不一致［44-48］。

3.3 tasiR-ARF的起源

粗梗水蕨的一个21 nt 的sRNA 受到多个干净序列的支持，并且与TAS3b的3′端miR390 切割位点的相位吻合（见图2）。这个sRNA的序列能与粗梗水蕨ARF2-like 基因高度互补，并且这些基因确实在二者的结合位点受到了切割［29］。其他蕨类也有MIR390、TAS3，以及含有该结合位点的ARF2-like基因。因此，蕨类植物确实存在类似于种子植物的tasiR-ARF。ARF2-like 基因在苔藓植物中受到其他分子的切割，而石松植物没有MIR390、其ARF2-like 基因也没有tasiR-ARF 结合位点［29］。因此，tasiR-ARF 可能只局限于真叶植物。TAS3a的一个预期切割物极其类似tasiR-ARF，只是可能在测序中没有被检测到（见图2）。果真如此，则ta⁃siR-ARF 来源于不同长度版本的TAS3基因在真叶植物的祖先里就已出现。