基于高通量测序的低温胁迫下赤桉miRNAs的挖掘与分析

张紫阳，刘 艳，魏瑞研，林元震

(华南农业大学 林学与风景园林学院/广东省森林植物种质创新与利用重点实验室，广东 广州 510642)

赤桉Eucalyptuscamaldulensis是桃金娘科Myrtaceae桉属Eucalyptus比较广泛栽培的树种，原产于澳大利亚，与其他桉树相比，赤桉具有高抗寒性，是研究桉树抗寒分子机理的理想材料[1]。

miRNA是由microRNA基因所编码、II型RNA 聚合酶 (Pol II)和 DCL1 酶 (Dicer-like 1)等加工形成的长度约18～24 nt的单链非编码RNA[2]。研究证实，miRNA通过对靶基因的剪切或抑制来实现转录后水平上的负调控，在植物的生长发育、形态建成和信号转导以及逆境胁迫响应中都发挥着重要的调控作用[3-5]。

低温是一种常见的环境因子，它限制植物只能在特定的地理位置和季节生长。植物为适应低温胁迫，在生理生化和分子水平上形成了一系列的抵御机制，包括瞬时Ca2+和脱落酸(ABA)含量的提高、植物体液组成的改变、抗氧化物含量的增加，以及一些渗透调节物质的集聚[6]。最早在拟南芥Arabidopsisthaliana中发现低温胁迫下miRNA的表达丰度会发生变化[7]，之后许多研究报道了木本植物低温胁迫下差异表达的miRNA[4,8]。华南农业大学桉树育种课题组前期对赤桉耐寒关键性转录因子ICE1(Inducer of CBF expression 1)[9]、SIZ1(SAP and Miz)苏素化[10]和 HOS1(High expression of osmotically responsive gene 1)泛素化[11]通路进行了一些研究，但在赤桉ICE1通路及其上下游中是否有miRNA参与调节，仍然未知。本研究拟通过高通量测序和生物信息学软件相结合的方式，挖掘、分析涉及低温胁迫应答的赤桉miRNA，以进一步弄清赤桉低温胁迫应答的完整分子调控网络，并为后续的桉树抗寒分子育种提供一些理论依据和技术支撑。

1 材料与方法

1.1 植物材料与处理

赤桉无性系CV103无菌组培苗取自广东省森林植物种质创新与利用重点实验室。在MS培养基 (pH=5.8)附加 0.5 mg·L−16-BA、0.1 mg·L−1NAA、30 g·L−1蔗糖和 6 g·L−1琼脂，并接种赤桉无菌苗，25 ℃ 条件下 16 h 光照、8 h 黑暗光周期培养3～4个月。低温处理组于4 ℃低温条件下处理24 h，对照组(CK)不进行处理。分别选取长势健壮无菌苗的茎尖组织，液氮速冻后磨碎，用Trizol(上海生工生物工程技术服务有限公司)提取总RNA。

1.2 Small RNA文库构建与测序

用 Qubit 2.0 (Thermo, 美国)对总 RNA 进行定量，使用 T4 RNA 连接酶 (NEB, 美国)将 RNA 3′端和5′端与接头连接，随后用 M-MuLV Reverse Transcriptase(NEB, 美国)合成 cDNA 单链，经PCR扩增后通过12 g·L−1聚丙烯酰胺凝胶电泳回收140～150 bp左右的PCR产物条带，通过质检后得到满足Illumina平台测序的文库，并委托上海生工生物工程技术服务有限公司进行高通量测序。

1.3 测序数据处理

测序下机的原始序列(Raw reads)经过去接头，去低质量、冗余序列等处理后得到纯净序列(Clean reads)，使用 Bowtie(参数为-v 1)将纯净序列匹配到巨桉E.grandis基因组上，得到的序列再与Rfam14.1数据库对比，过滤其他非编码RNA。

1.4 miRNA的预测

以miRBase21.0库中植物和巨桉[12]的miRNA序列为参考源，使用blastn软件预测已知miRNA。具体参数为：-task blastn-short-word_size 16-evalue 0.01。将错配数≤1、e＜0.01的序列作为潜在miRNA。同时，将这些潜在的miRNA 用Bowtie(参数为-v 2)匹配回巨桉基因组后，通过miREAP和RNAfold软件预测所有已知miRNA及其二级结构。未匹配上的序列，继续使用miRDeep2(v2.0.0.8)[13]软件预测新miRNA，采用默认参数。

对不同长度miRNA的首位碱基进行统计并绘制直方图，分析赤桉miRNA碱基的偏好性。本文除了miRNA二级结构采用RNAfold软件生成，其他所有图形均使用R语言绘制[14]。

1.5 差异表达miRNA的挖掘与分析

低温处理下赤桉差异表达的miRNA采用Bioconductor中的DEGSeq[15]包分析。首先将2组数据使用TPM进行归一化处理，归一化表达量=(Readcount×1 000 000)/libsize，其中：Readcount为miRNA序列数目，libsize为样品所有miRNA 序列数目之和。设置DEGSeq的算法为MARS，进行差异分析。定义表达量变化2倍以上且P<0.001的为显著差异表达miRNA。

1.6 miRNA靶基因预测与富集分析

使用psRNATarget[16]软件对miRNA进行靶基因预测，采用默认参数。同时对差异表达miRNA的靶基因分别进行 Gene ontology[17](http://www.geneontology.org/) 和 KEGG[18](http://www.genome.jp/kegg/)富集分析。GO注释使用在线网站 agriGO v2.0[19](http://systemsbiology.cau.edu.cn/agriGOv2/index.php)，并定义P<0.05为显著差异的类别。KEGG使用KOBAS3.0[20]软件，该软件以KEGG通路为单位，应用超几何检验进行分析：

式中，p为KEGG通路富集概率，i为求和中的特定项，N为所有基因中具有通路注释的基因数目，n为N中候选靶基因的数目，M为所有基因中注释为某特定通路的基因数目，m为注释为某特定通路的候选靶基因数目。同时，以拟南芥的低温胁迫响应网络和华南农业大学桉树育种课题组前期研究为参考，进一步挖掘涉及ICE1-CBFs-COR通路的miRNA。

2 结果与分析

2.1 测序数据分析

对照组和低温处理组的小RNA高通量测序分别得到 29 725 524 和 31 104 120 条原始序列，过滤掉低质量、两端接头和重复序列后，最后分别得到15 813 414 和 14 456 077 条纯净序列。质控数据显示，低温处理组的Q30为96.47%，对照组的Q30为96.46%，结果可靠性高。此外，纯净序列的长度统计结果(图1)显示：2组小RNA长度的分布特征基本一致，以21 nt分布最多，其次是20和24 nt。

图 1 赤桉小RNA纯净序列的长度分布图Fig.1 Length distribution of clean reads of small RNA in Eucalyptus camaldulensis

2.2 赤桉已知miRNA和新miRNA的预测

以RFam数据库对纯净序列进行注释，结果如表1所示，除了其他类别外，核糖体RNA(rRNA)占比最高，在对照和低温处理组中分别占比30.27%和26.08%，而核仁小RNA (snoRNA)占比最少，分别占比0.21%和0.17%。过滤掉rRNA、snRNA、snoRNA和tRNA序列后，以巨桉和miRBase21中所有植物的成熟miRNA为参考组，进行miRNA预测，结果如图2所示：2组数据一共预测到392个已知miRNAs(隶属于 54个家族)，其中对照组282个，低温处理组329个，它们共同的有219个；预测到97个新miRNA，其中对照组65个，低温处理组51个，它们共同的有19个。在预测到的所有miRNA中，长度均以 21 nt为主，其次，已知miRNA为20 和 22 nt，新 miRNA为24 nt(图 3)。图4展示了部分预测miRNA前体的二级结构。

表 1 赤桉小RNA分类统计Table 1 Classification statistics of small RNAs in Eucalyptus camaldulensis

图 2 赤桉miRNA维恩图Fig.2 Venn diagram of miRNAs in Eucalyptus camaldulensis

图 3 赤桉预测miRNA长度分布Fig.3 The length distribution of the predicted miRNAs in Eucalyptus camaldulensis

图 4 赤桉部分预测miRNA前体二级结构Fig.4 The secondary structure of some predicted miRNAs in Eucalyptus camaldulensis

miRNA是由DCL1酶剪切得到的，DCL1酶在剪切pre-miRNA时，对其5′端第1个碱基具有强烈的偏好性，尤其体现在对尿嘧啶(U)的偏好上。miRNA首位碱基偏好性结果表明，对照和低温处理组均呈现出相似规律，具体为：长度为18 nt的miRNA其主要为腺嘌呤(A)，19～22 nt的miRNA其主要为尿嘧啶(U)，而23～26 nt的miRNA其以腺嘌呤(A)或鸟嘌呤(G)为主(图5)。

图 5 赤桉不同长度miRNA首位碱基偏好性Fig.5 Bias of the first base of different length miRNAs in Eucalyptus camaldulensis

2.3 差异表达miRNA分析

以TPM标准化后的序列数为源数据，使用DEGSeq包分析对照组和低温处理组2个样品间差异表达的miRNA。共筛选出80个显著差异表达的miRNA(图6)，包括55个显著上调(含10个新miRNA)和25个显著下调的miRNA(含16个新miRNA)。上调的miRNA归属于miR395、miR167、miR399、miR8282、miR477、miR156、miR171 和 miR482(eca-miR482a-3p, eca-miR482d-5p)等 8 个家族，而下调的miRNA归属于miR862、miR482(ecamiR482b-5p)、miR6300、miR397和 miR10506 等5个家族。

图 6 赤桉差异表达miRNA火山图Fig.6 Volcanic map of differentially expressed miRNAs in Eucalyptus camaldulensis

2.4 靶基因预测和基因富集分析

默认参数下，psRNATarget软件共预测到17 614个靶基因，其中低温胁迫下差异表达miRNA的靶基因有 5 983个。这些差异表达miRNA的靶基因GO注释结果(图7)显示，定位到87个功能群中，包括分子功能 (Molecular function, MF)49 个、生物学过程 (Biological process, BP)31 个和细胞组成(Cell component)7个。这些靶基因中最多的一类是分子功能，主要集中在结合作用；在生物学过程中，主要集中在单一生物过程、生物调节和信号转导；在细胞成分中，靶基因主要与细胞膜有关。KEGG富集分析结果表明，这些靶基因参与的通路有23条，其中代谢通路(Metabolic pathway)和次生代谢物的生物合成 (Biosynthesis of secondary metabolite)最多，分别占到了36.7%和17.8%(图 8)。

图 7 赤桉差异表达miRNA靶基因的GO注释Fig.7 GO annotation of the target genes of differentially expressed miRNAs in Eucalyptus camaldulensis

图 8 赤桉差异表达miRNA靶基因KEGG富集分析Fig.8 KEGG pathway analysis of the target genes of differentially expressed miRNAs in Eucalyptus camaldulensis

2.5 赤桉ICE1-CBFs-COR通路相关miRNA的筛选

ICE1-CBFs-COR通路是植物低温胁迫应答的重要信号转导途径[21]，本试验进一步筛选到25个与该通路相关的miRNA(表2)。这些靶基因主要编码蛋白激酶、E3泛素化连接酶和转录因子，具体包含：ICE1，直接调控ICE1的MPK6(Mitogen-activated protein kinase 6)、OST1(Open stomata 1)、HOS1、SIZ1，与调控ICE1下游的CBF1～CBF3有关的BTF3(Basic transcription factor 3)、CAMTA3(Calmodulin-binding transcription activator 3)、BZR1(Brassinazole-resistant 1)、EIN3(Ethylene in sensitive 3)、PIF3(Phytochrome-interacting factor 3)和SOC1(Suppressor of constans overexpression 1)等。其中，仅eca-miR-n60有2个靶基因，其余均有1个靶基因；多个miRNA靶向同一个基因，比如eca-miR5780b-5p、eca-miR390b-5p、eca-miR-n51 和eca-miR-n60均靶向OST1。这些miRNA长度主要为20～21 nt，包括 8 个保守型 miRNA(miR156、miR159、miR164、miR171 和 miR482 等家族)、10个非保守型miRNA以及7个新miRNA。此外，有3个 miRNA(eca-miR-n51，eca-miR-n41和 ecamiR862a-5p)经过4 ℃低温处理24 h后显著差异表达，其中，eca-miR-n51靶向OST1，而eca-miR-n41和eca-miR862a-5p均靶向EBF1，其余的miRNA则差异不显著。

表 2 赤桉ICE1-CBFs-COR通路相关miRNATable 2 The miRNAs associated with ICE1-CBFs-COR pathway in Eucalyptus camaldulensis

续表 2 Continued table 2

3 讨论与结论

本试验通过高通量测序共预测到隶属于54个家族的392个已知miRNA和97个新miRNA，并通过软件预测到17 614个靶基因。与对照组相比，低温处理组在4 ℃处理24 h后，预测到80个显著差异表达的miRNA，包括55个上调miRNA和25个下调miRNA。GO基因功能注释表明，低温胁迫下miRNA的靶基因主要与结合作用、单一生物过程、生物调节和信号转导以及细胞膜有关；KEGG富集分析表明，差异表达miRNA靶基因可能在代谢通路和次生代谢物的生物合成中起关键作用。同时，本试验还分析了可能参与到赤桉ICE1-CBFs-COR通路有关的miRNA。

在miRNA预测过程中，发现Lin等[12]的巨桉成熟miRNA序列存在一些序列无法匹配到巨桉基因组的现象，因此我们将Lin等[12]的miRNA序列与其他植物的miRNA整合后，作为miRNA参考组，再与测序下机后的纯净序列匹配，得到潜在的已知miRNA。随后，将这些潜在已知miRNA匹配回巨桉基因组，再使用miREAP和RNAFold软件寻找前体序列。但是，由于赤桉和巨桉属于近缘种，基于巨桉基因组预测的miRNA有可能存在假阳性，后续仍需试验进一步验证。

在植物中，行使重要功能的miRNA，往往是保守的。本试验中，miR156、miR167、miR395和miR399家族的多个成员在低温处理后均显著差异表达。Zhou等[22]指出，miR156家族广泛响应植物低温胁迫，比如低温胁迫会诱导杨树miR156表达量下调[23-24]，而甘蔗抗寒品种‘FN39’在低温胁迫下miR156表达量上调[25]。研究发现，miR167和ARF(Auxin response factors)有关，并且 miR167 通过调控生长素水平抵御冷胁迫[26]。薄维平等[27]报道，低温胁迫会诱导木薯miR395的上调表达，其中miR395abcd的上调有利于减轻木薯的低温伤害。Gao等[28]发现在番茄中过表达拟南芥athmiR399d可以提高番茄的耐寒性。除了这些保守型miRNA，本试验中其他显著差异表达miRNA在其他植物低温胁迫中也有过类似报道[29]。

ICE1作为关键性抗寒转录因子，其主导的ICE1-CBFs-COR途径是植物响应低温胁迫的重要通路。在该通路中，ICE1通过与MYB15[30]或HOS1[31]的互作来负调控CBFs及其下游基因，但其通过与SIZ1[32]或OST1[33]的互作来实现正调控。最近，Li等[34]报道MPK3/MPK6可与ICE1相互作用并将其磷酸化，降低ICE1的稳定性及其转录活性，从而负调控拟南芥CBF表达及其抗寒力。本试验进一步筛选到25个与该通路相关的miRNA，其中，eca-miR-n33靶向ICE1，可正调控ICE1-CBFs-COR通路的miRNA有靶向OST1的ecamiR5780b-5p、eca-miR390b-5p、eca-miR-n51 和 ecamiR-n60以及靶向SIZ1的eca-miR482f-3p，而负调控该通路的miRNA有靶向HOS1的eca-miR23a-5p和eca-miR-n2以及靶向MPK6的eca-miR171g-3p和eca-miR171j-3p。此外，Li等[35]报道BZR1可以正调控CBF1/2，增强植物抵御低温胁迫的能力。Jiang等[36]报道与光形态建成有关的转录因子PIF3可与CBF启动子直接结合，实现对低温胁迫的负调控；而 EBF1(EIN3-BINDING F-BOX 1)可直接靶向PIF3进行26S蛋白酶体介导的降解，提高植物的抗寒力。综上可知，上述的这些miRNA可为后续深入研究其参与调控赤桉低温胁迫应答的分子网络奠定基础。