甘蔗耐寒相关miRNA的生物信息学分析及其靶基因预测

朱鹏锦， 宋奇琦， 谭秦亮， 程 琴， 李佳慧， 庞新华， 周全光， 吕 平， 欧克纬， 卢业飞， 农泽梅， 唐桓伟， 龙盛风

(广西壮族自治区亚热带作物研究所， 南宁 530001 )

miRNA是一类含有20～24个碱基的小分子内源性非编码RNA，具有高度的保守性、时序性和组织特异性(Ambros, 2004)。miRNA主要在转录后水平上通过介导靶基因mRNA的切割或抑制翻译来调节基因的表达，在生物体代谢过程中起到多种调控作用，如参与调控植物器官的形态建成(Sunkar et al., 2012)、生长发育(Thiebaut, et al.， 2012)、激素分泌、信号转导以及对外界环境胁迫(Xiong & Zhu, 2003)的响应等过程。在拟南芥和水稻响应低温的研究中发现，miR-167、miR-169、miR-319和miR-171等miRNA家族在低温响应具有重要的生物学作用(Wang et al., 2010)，其中miR-169的靶基因为低温诱导的重要基因CBF (Sunkar et al., 2007)；还在水稻响应低温胁迫的研究中证明，miRNA-319和miRNA-171的靶基因属于MYB类的转录因子，两者的表达量互为消长关系，进而说明miRNA在水稻耐低温途径中所起的调控作用(Lü et al., 2010)。因此，准确高效分离和鉴定miRNA及其靶基因并分析它们的功能，更精确地掌握miRNA在植物抗逆胁迫过程中的调控机制，是目前植物miRNA研究领域的重要内容。

甘蔗是热带、亚热带地区重要的糖料作物，低温是限制其扩大种植区域和实现高产稳产的重要因素之一(吴棉国等, 2010)。出现的偶发性大范围冰冻雨雪或严重的霜冻天气会导致我国甘蔗主产区遭受巨大的经济损失(何燕等, 2009; 匡昭敏等, 2009)；广西地理环境特殊，冬季持续低温以及春季“倒春寒”形成的阴雨霜冻频发导致大面积甘蔗受到冷害，具体表现在甘蔗叶片枯萎、茎杆坏死以致蔗糖含量、甘蔗产量下降(李杨瑞等, 2011)。古丽等(2011)采集甘蔗主产市县40多年的气象数据，结合甘蔗种植面积、甘蔗产量和蔗糖产量等指标进行分析发现，低温霜冻灾害是影响甘蔗种植、甘蔗生产和蔗糖产量的主要环境因子。因此，了解甘蔗的耐低温调控机制是培育耐寒性强或适合我国热带北缘气候的甘蔗品种的前提。

本研究在观察不同甘蔗品种的田间农艺性状和模拟低温胁迫的生理生化研究基础上，筛选出抗寒能力较强的甘蔗品种作为材料，利用高通量测序技术及生物信息学方法，获得与甘蔗响应低温胁迫相关的miRNA，分析miRNA的差异表达情况，明确miRNA与靶基因的作用关系，并对预测所得的靶基因进行基因本体(gene ontology, GO)分析，挖掘低温胁迫应答基因，为选育耐寒性强的优良甘蔗新品种提供理论依据和技术支撑。

1 材料与方法

1.1 材料

以广西壮族自治区农业科学院甘蔗研究所选育的‘桂糖28号’(GT28)、广西蔗区主栽品种‘新台糖22号’(ROC22)和广西壮族自治区亚热带作物研究所新选育的‘桂热2号’(GR2)为材料。

1.2 方法

1.2.1 实验设计与管理 选择无病虫害、蔗茎大小均匀的种茎切成单芽段，用清水冲洗种茎、干布擦洗干净；用50%多菌灵可湿性粉剂1 000倍液浸12 h消毒；用蒸馏水浸洗种茎1 min后，再用蒸馏水浸洗拧干的棉布将种茎包好并做好标记；用橡皮筋扎好放进温度为25 ℃的恒温箱催芽。当种茎萌芽并长出幼根时，将其移植到装有营养土的塑料盆并做好标记。每盆1段种茎，盆高17.5 cm，盆宽16 cm。育苗期间，每株施用完全营养液2次，每次10 mL。当幼苗两叶一心时，选择长势、大小均一的幼苗进行低温处理。低温胁迫处理温度为4 ℃，光照强度为5 000 lx，处理24 h；对照(CK)温度为28 ℃，光照强度为5 000 lx。

1.2.2 RNA 提取及 Illumina 测序 采集低温处理和对照样本的叶片，及时用TRNzol Reagent试剂盒(天根生化科技有限公司, 北京)提取并检测RNA。将低温胁迫处理组与对照样本组提取获得的高质量RNA送往华大基因科技有限公司(深圳)进行文库构建，每处理3个重复。采用Illumina HiSeqTM 2000平台对质量合格的RNA文库进行测序(Ali et al., 2008)。通过去接头、去低质量、去污染等过程完成数据处理获得干净序列后，使用Bowtie 2软件按照MiRbase > pirnabank > snoRNA(human/plant) > Rfam > other sRNA的优先级顺序将小RNA(sRNA)遍历注释获得未注释RNA片段。

1.2.3已知miRNA的鉴定及新miRNA预测 用miRDeep 2软件对所得的未注释sRNA序列与参考基因组 [割手密(Saccharumspontaneum)基因组]进行序列比对分析 (Maurits et al., 2015)，鉴定出已知的甘蔗耐寒相关miRNA；将过滤获得的未比对序列比对参考序列，通过碱基数目延伸、miRNA结构预测方法获得新的miRNA。

1.2.4 差异表达miRNA的鉴定 对耐寒型甘蔗样品的已知miRNA的读数进行分析，判断低温胁迫前后不同耐寒型甘蔗样品中miRNA的差异表达，以|log 2(FoldChange)|≥1，P<0.05为筛选标准，获得低温胁迫前后差异表达 miRNA。

1.2.5 miRNA靶基因预测 用psRNATarget(Wu et al., 2012)、TargetFinder(Fahlgren & Carrington, 2010)和Tapirhybrid(Peer，2010)软件进行差异表达miRNA的靶基因预测，并在GO数据库对预测靶基因进行同源性搜索，确定其参与的信号传导及生物代谢途径。

1.2.6 实时荧光定量PCR验证 选取14个差异表达miRNA及其靶基因，先用Premier 5.0设计其成熟miRNA特异正向引物、通用反向引物(表4)，再用LightCycler®480 Instrument II进行 RT-PCR扩增，设阴性对照(不添加cDNA模板)以监控可能的污染。以甘蔗GAPDH为内参基因, 每个样品均设置3次重复, 采用2-ΔΔCt法计算基因相对表达。

2 结果与分析

2.1 高通量测序数据分析

选择3个耐寒性不同甘蔗品种的幼苗进行低温胁迫，每个品种的对照(CK)和低温处理(T)3个重复，分别取叶片共18个样本进行高通量测序，构建低温胁迫前后sRNA文库；对原始测序数据进行3′端去接头、Trim低质量、片段大小选择等质控处理，获取高质量clean read序列(表1)。从表1可以看出，在GR2的CK和T叶片中分别挖掘出24 310 558、23 925 673条测序数据，最终处理分别得到21 343 194、21 576 594条，高质量clean read序列分别占总序列的87.83%、90.25%，且所构建的文库与参考基因组进行比对，其CK和T分别有86.71%与84.32%的序列与参考基因组匹配；GT28的CK和T叶片中分别挖掘出23 915 244、23 528 103条测序数据，最终处理分别得到22 045 186、21 481 129条，高质量clean read序列分别占总序列的92.18%与91.28%，所构建的文库与参考基因组进行比对，其CK和T分别有71.75%与84.10%的序列与参考基因组匹配；ROCC22的CK和T叶片中分别挖掘出27 069 867、23 631 208条测序数据，最终处理分别得到21 007 371、20 640 357条，高质量clean read序列占总序列的79.50%与87.35%，所构建的文库与参考基因组进行比对，其CK和T分别有92.62%与85.67%的序列与参考基因组匹配(表1)。

表 1 样品测序数据统计Table 1 Sequencing data statistics of samples

sRNA长度与不同功能有关，其中21～22 nt的主要与mRNA切割和转录后基因沉默相关，24 nt的主要与RNA导向的DNA甲基化和转录基因沉默相关。对total clean reads进行统计发现，sRNA主要集中在21～24 nt之间，但不同抗寒能力的材料之间有一定差异，且不同长度的测序频率不同(图1)。

图 1 sRNA长度分布图Fig. 1 Graphs sRNA length distribution

根据基因表达量进一步对每个样品进行Pearson相关系数计算，并将这些系数以热图的形式反映出来(图2)。从Pearson相关系数可以看出，同一品种相同处理的重复样品间相关性较高(0.667～0.990)，表明所有样品的表达量基本保持一致，说明其符合样本重复性实验标准，可以满足后续的差异表达分析。

X、Y轴均代表每个样品。颜色代表相关性系数，颜色越蓝代表相关性越高，颜色越浅代表相关性越低。X and Y axes represent each sample. Color represents the correlation coefficient, the bluer the color, the higher the correlation, and the lighter the color, the lower the correlation.图 2 样品间相关性分析热图Fig. 2 Correlation heatmap of different samples

2.2 已知miRNA分析及新miRNA预测

用AASRA软件将核苷酸序列比对到参考基因组及miRBase数据库，鉴定得到分属于84个已知miRNA家族的322个miRNA，其中对照组的297个miRNA分属于69个家族，低温处理组的305个miRNA分属于74个家族，这些家族中拥有最多家族成员数量的为miR169(32个)，其次分别为miR166、miR171、miR167、miR156、miR396 (图3)。

图 3 miRNA家族数量分布前18位Fig. 3 Top 18 miRNA family members

根据成熟植物miRNA序列的高度保守型及miRNA前体拥有标志性发夹结构的特性，利用同源搜寻比对的方法在割手密miRBase数据库中进行搜索，对具有同源性的miRNA序列进行提取，选取其中表达量最高的作为甘蔗响应低温的保守miRNA的候选者，最终得到甘蔗响应低温的110个新miRNA。

2.3 差异表达miRNA分析

通过比较对照组与低温处理组miRNA的表达量变化情况，可以判断低温胁迫下抗寒能力不同的甘蔗品种miRNA的差异表达情况。在差异表达miRNA 检测过程中，以|log 2(FC)|≥2，P<0.05作为筛选标准，低温胁迫时差异表达miRNA 如表2所示，包括 100个已知miRNA (61个上调，39个下调)，37个新miRNA (15个上调，22个下调)。呈下调趋势的miRNA家族包括miR8175、miR5564、miR444、miR166在内的21个家族，呈上调趋势的miRNA家族包括miR156、 miR169、miR172、miR393、miR397、miR408等10个家族。这些具有差异性表达的miRNA在甘蔗响应低温胁迫情况下可能发挥特定的功能。

表 2 甘蔗响应低温差异表达miRNATable 2 Differential expressions of miRNAs in sugarcane responding to low temperature

续表2

续表2

2.4 miRNA作用靶基因预测与分析

根据筛选差异表达的miRNA与对应物种的基因序列信息，使用psRNATarget、TargetFinder和Tapurhybrid软件进行靶基因预测，结果见图4。3种预测方法共预测了1 844个相同的潜在靶基因，其中1 696个属于已知miRNA靶基因，148个属于新miRNA靶基因。为进一步分析预测所得靶基因的潜在生物学功能，对预测所得的靶基因进行GO分析，共鉴定出甘蔗响应低温胁迫miRNA靶基因：在生物过程(biological process)有13个功能亚类，在细胞组分(cellular component)有11个功能亚类，在分子功能(molecular function)有7个功能亚类(图4)。在生物过程类别中，预测靶基因的主要生物学功能富集于细胞过程和代谢过程；在细胞组分类别中，预测靶基因生物学功能主要富集于细胞、细胞器和生物膜；在分子功能类别中，预测靶基因生物学功能主要富集于结合和催化活性及转运活性。大多靶基因功能均与这些结合功能及其他相近的结合功能相关。大多miRNA通过直接或间接介导靶基因的表达调控相关代谢途径响应低温胁迫，这些miRNA所调控的靶基因对甘蔗的耐寒性起关键的调控作用(表 3)。

1. 细胞过程; 2. 代谢过程; 3.生物调节; 4. 发育过程; 5. 定位; 6. 响应刺激; 7. 多细胞有机体过程; 8. 繁殖; 9. 繁殖过程; 10. 信号传导; 11. 细胞成分组织或生物发生; 12. 多有机体过程; 13. 生长; 14. 细胞; 15. 细胞器; 16. 细胞膜; 17. 细胞膜部分; 18. 蛋白复合体; 19. 细胞器部分; 20. 细胞连接; 21. 共质体; 22. 膜包管腔; 23. 超分子复合体; 24. 胞外区; 25. 结合; 26. 转录因子活性; 27. 催化活性; 28. 转运活性; 29. 结构分子活性; 30. 分子功能因子; 31. 抗氧化活性.1. Cellular process; 2. Metabolic process; 3. Biological regulation; 4. Developmental process; 5. Localization; 6. Response to stimulus; 7. Multicellular organismal process; 8. Reproduction; 9. Reproductive process; 10. Signaling; 11. Cellular component organization or biogenesis; 12. Multi-organism process; 13. Growth; 14. Cell; 15. Organelle; 16. Membrane; 17. Membrane part; 18. Protein-containing complex; 19. Organelle part; 20. Cell junction; 21. Symplast; 22. Membrane-enclosed lumen; 23. Supramolecular complex; 24. Extracellular region; 25. Binding; 26. Transcription regulator activity; 27. Catalytic activity; 28. Transporter activity; 29. Structural molecule activity; 30. Molecular function regulator; 31. Antioxidant activity.图 4 差异表达miRNA靶基因的GO功能分析Fig. 4 GO function analysis of miRNAs target genes for differentially expressed

表 3 响应低温的部分甘蔗miRNA及其靶基因Table 3 Some miRNAs and their target genes in sugarcane responding to low temperature

2.6 qRT-PCR验证

利用qRT-PCR技术对测序所得到的miRNA及靶基因的差异丰度进行验证，miRNA及其靶基因的引物序列如表4所示。在对照组与低温处理组中共筛选出14个差异表达的miRNA及其靶基因，并对它们进行qRT-PCR验证，包括miR156、miR160g_1、miR167d、miR169 a-3p_3、miR171b-3p_3、miR172d-5p_4、miR393-3p_1、miR396b、 miR397a_3、miR398b、miR399k_1、miR408d、novel_mir36、novel_mir89(图5)。根据表4和图5的结果，除novel-miR36外，选择的其余13个miRNA在qRT-PCR实验中的表达模式与高通量测序检测到的表达模式一致(图5：A)，这表明大部分高通量测序结果均能被荧光定量 PCR 技术所验证，此次研究数据为真实可信的。所筛选的14个miRNA，除novel-miR36外，其余13个miRNA均与其靶基因呈负调控关系(图5：B)，证明miRNA通常以负调控或沉默靶基因的方式来调控植物的生长、发育及对环境的应激反应。

表 4 miRNA和靶基因qRT-PCR引物设计Table 4 Design of qRT-PCR primer for miRNAs and target genes

A. miRNA qRT-PCR 检验; B. 靶基因qRT-PCR 检验。A. Results of miRNA qRT-PCR; B. Results of target genes.图 5 miRNA及其靶基因qRT-PCR分析Fig. 5 qRT-PCR analysis of miRNA and target genes

3 讨论与结论

甘蔗热带种原产热带地区，喜温，现代甘蔗栽培品种以热带种的种质为主体。作为广西传统的主导产业，而低温不仅是限制其扩大种植区域和实现高产稳产的重要因素之一，还影响蔗农收益和糖业稳定发展(苏永秀等, 2006)。为了解甘蔗响应低温的内在分子机制，挖掘其与耐寒相关的miRNA及相关靶基因，本研究对不同基因型甘蔗进行低温胁迫处理，通过高通量测序技术及生物信息学方法，系统分析不同抗寒能力的甘蔗对低温胁迫的响应。

通过高通量测序技术及生物信息学分析发现，sRNA主要集中在21～24 nt之间，不同抗寒能力的材料之间有一定差异且不同长度的测序频率不同。sRNA在不同物种间的长度分布会有所区别，如拟南芥(Pasquinelli et al., 2000)、小麦(Meng et al., 2013)、棉花(Sripathi et al., 2014)的sRNA长度最多分布在24 nt，杨树(李明娜等, 2014)、 大豆 (Turner et al., 2012)和番茄 (Pilcheret al., 2007)的sRNA长度最多分布在21 nt，这些结果与本研究结果高度相似。本研究通过生物信息学手段共挖掘137个miRNA在低温胁迫前后进行差异性表达，其中100个已知miRNA (61个上调，39个下调)，37个新miRNA(15个上调, 22个下调)。在低温胁迫下植物通过调节miRNA的表达水平，进而调节对应靶基因的表达，从而引起相关代谢与信号转导途径的变化来实现对逆境的响应，其信号转导途径主要包括胞外信号途径、胞内第二信使、转录因子以及功能基因等。Wu和Poethig (2006)的研究发现，低温胁迫导致miR156下调而其靶基因表达量增加， 进而调控拟南芥营和养期延长、生长代谢变缓，以应对不良环境，这与本研究低温胁迫后miR156下调的结果一致。王丽丽等(2017)在拟南芥、小麦和水稻研究中发现miR160响应低温胁迫。而本研究中，在低温胁迫后miR160g_1表达下调，靶基因表达受抑制且作用于生长素信号通路从而在抵抗低温中发挥作用。Pourcel等( 2005)研究表明miR397的靶基因与漆酶有关，与细胞壁木质素的合成、抗病、对环境的适应过程密切相关。而本研究发现，低温胁迫后甘蔗叶片miR397a_3呈上调表达，其靶基因与抗坏血酸氧化酶(L-ascorbate oxidase)有关，可能参与调控抗坏血酸氧化酶的表达而增强低温响应能力。Li等(2014)研究旱芹miRNA对低温胁迫的响应发现， miR160、 miR164、miR394、 miR395、miR408具有表达差异。Sunkar和Zhu (2004)在研究分析拟南芥响应胁迫miRNA中，发现一些miRNA能被多种胁迫因素诱导，如miR393受到低温、干旱、高盐的诱导表达。Gupta等(2014)研究小麦miRNA对低温胁迫、盐胁迫、渗透胁迫的响应发现，在盐胁迫和低温胁迫下，miR168、miR397均表达下调，而miR172表达上调；miR393在渗透和盐胁迫下表达量上升，在低温胁迫下表达量下降。Sun等(2015) 研究发现，在低温胁迫下对葡萄的miR169的表达量上调，但在拟南芥、扁桃中研究发现miR169的表达量下调。结合前人的研究结果，本研究认为特定miRNA对低温胁迫的响应可能因植物种类、同一植物的不同基因型、不同的组织类型、胁迫的时间等而有所差异。

本研究利用psRNATarget、TargetFinder 及Tapirhybrid三种分析方式预测靶基因，对预测所得的靶基因进行GO分析发现，在生物过程类别中预测靶基因的主要生物学功能富集于细胞过程和代谢过程；在细胞组分类别中预测靶基因生物学功能主要富集于细胞、细胞器和生物膜；在分子功能类别中预测靶基因生物学功能主要富集于结合和催化活性及转运活性。在逆境胁迫条件下植物通过调节miRNA的表达水平，进而调节对应靶基因的表达，从而引起相关代谢与信号转导途径的变化来实现对逆境的响应。因此，进一步对参与调控植物激素信号传导、光合色素合成、抗氧化酶系统、泛素介导蛋白水解、淀粉与蔗糖代谢等通路的靶基因进行qRT-PCR验证。本研究发现，miR156靶向调控SBP转录因子，低温胁迫前后，miR156下调表达，负向调控SBP转录因子，使甘蔗生长代谢变缓而耐寒能力增强。梅琳(2007)研究冬小麦发现，低温胁迫下miR160的靶基因ARF17在过表达miR160f拟南芥植株中表达量明显下降。在拟南芥、水稻、玉米研究中证明miR160、miR167的靶基因是ARF(AuxinResponseFactors), 主要通过调控生长素信号通路应对逆境。而本研究也发现甘蔗在低温胁迫后miR160g_1上调表达且其作用于转录因子ARF。张译云等(2012)研究发现毛白杨在低温胁迫下miR169ac表达下调，对其靶基因转录因子NAC负向调控。党春艳(2013)在研究高山离子芥低温胁迫响应中发现，miR169a在冷胁迫下无显著表达，而miR169在冷胁迫下表达下调，表明miR169家族之间响应低温存在差异。本研究中，miR169 a-3p_3在低温胁迫后表达上调，作用于细胞核转录因子使其表达下调，推测miR169 a-3p_3与调控甘蔗低温耐受能力密切相关。

此外，本研究通过生物信息学挖掘甘蔗响应低温的miRNA，不仅参与抗氧化酶系统、植物激素信号传导、遗传信息等，而且参与类胡萝卜素代谢、卟啉与叶绿素代谢、淀粉与蔗糖代谢途径相关基因的表达。有研究表明叶绿体内的类胡萝卜素既可作为吸收光能的辅助色素将能量传递给叶绿素a，又可在持续胁迫条件下将过剩光能安全耗散而保护光合机构(Liu et al., 2004; Holt et al., 2005)。我们前期的研究也发现，低温胁迫时甘蔗能够通过增加光合机构过剩激发能的耗散和调整其叶片光合色素含量与构成进行有效的光能利用和分配。本研究对部分差异表达miRNA及其靶基因进行qRT-PCR检验发现，参与类胡萝卜素代谢的miR167d对靶基因具有负向调控作用。因此，后续我们将对低温条件下参与光合生理过程的miRNA分子调控机制进行深入研究，为抗寒育种提供重要理论依据和实验基础。