分布于我国华南地区的野生蕉——阿宽蕉（Musa itinerans）是一类观叶、观果的芭蕉属植物,广泛应用于园林绿化中。福建野生蕉资源非常丰富,许多野生蕉因富含花青素[1]而具有紫红色的果实,观赏性很强,为芭蕉属植物观赏性状改良提供了优异的基因资源。花青素作为一种水溶性天然可食用色素,不仅在植物显色及生长发育中起作用,还具有医疗保健功效,是植物观赏性状和品质改良的一个重要研究内容。随着分子生物学的快速发展,植物花青素的生物合成途径已经比较明确,但花青素生物合成调控的分子机制非常复杂,除了受结构基因和转录因子的共同调控外,近年来的许多研究也表明miRNAs可以通过靶向调控花青素合成的MYB、bHLH和WD40等多种转录因子,实现对花青素的转录后水平调控。如在拟南芥中,miR156的靶基因SPL9通过降低MBW复合物的稳定性来抑制花青素的合成和积累[2];miR827通过调控靶基因NLA参与花青素的生物合成[3];miR858通过靶基因MYB12调控黄酮类物质的生成途径来影响花青素的生物合成[4,5];miR828通过对TAS4剪切产生许多siRNA,共同靶向花青素合成相关的PAP1/MYB75、PAP2/MYB90、MYB113,从而影响花青素的合成[6,7]。此后,有越来越多的研究表明miRNA参与植物花青素的转录后调控[8,9]。
1 材料与方法
1.1 材料
系列试验所用材料均为种植于福建农林大学园艺植物生物工程研究所试验地的闽侯野生蕉[11]。该野生蕉为通信作者于2011年在闽侯三叠井采集后移栽至试验基地。经过几年的观察,果皮紫色性状稳定。试验取同一果树上不同颜色（绿色和紫色）的果皮组织,用编号GP表示绿色果皮,PP表示紫色果皮,见图1。
1.2 方法
1.2.1 样品采集、RNA的提取、文库构建与测序 取树上刚采下的野生蕉果梳,用削皮刀剥下表面带颜色的果皮,在液氮下迅速研磨成粉末,用北京天恩泽公司生产的柱式植物RNAOUT 2.0试剂盒（Column Plant RNAOUT 2.0 Kit）提取总RNA后,交由北京百迈客生物科技有限公司进行文库构建,并在Illumina Hiseq 2500平台进行测序。
1.2.2 数据处理 （1）将测序得到的原始序列去除接头序列和短于18 nt或长于30 nt的序列等低质量序列,获得高质量序列（即clean reads）用于后续分析。（2）Bowtie（http://bowtie-bio. Sourceforge.net/）,用于序列比对,获得包含miRNA的Unannotated reads。（3）将Unannotated reads与植物miRNA数据库（miRbase 20.0,http://www. MiRbase.org/）中的已知miRNA进行比对,获得保守miRNA的数目及家族分布。同时使用Musa_Itinerans_v1作为参考基因组（http:// banana- genome-hub.southgreen.fr/organism/ Musa/Itinerans）进行新的miRNA预测、序列比对和后续分析。（4）用TPM算法归一化处理miRNA的表达量,对2个样本中miRNA表达量进行统计。差异表达miRNA的筛选以校正后的p-value即错误发现率FDR（false discovery rate）≤0.05,2个样本间表达量差异倍数FC在2倍以上即|log2(FC)|≥1作为条件。（5）对得到的miRNA（已知的miRNA和新预测的miRNA）用Target Finder软件进行靶基因预测。（6）用DAVID 软件（http://david.abcc. Ncifcrf. gov/）进行KEGG通路分析和GO分析,用富集因子（enrichment factor）分析差异表达miRNA靶基因的KEGG通路富集程度。
1.2.3 cDNA及引物的合成 对提取的总RNA用北京全式金生物技术有限公司生产的miRNA反转录试剂盒（TransScript miRNA RT Enzyme Mix）进行cDNA的合成。用ABI公司的Primer Express Software v2.0设计候选基因用于荧光定量PCR的引物,委托深圳华大基因公司合成引物,引物序列信息见表1。
表1 荧光定量引物及内参Tab. 1 The primers of miRNA and reference gene |
| miRNA | 上游引物序列 Primer sequences (5°-3°) | miRNA序列 miRNA sequences (5°-3°) |
|---|---|---|
| mit-miR156k | TTGACAGAAGAGAGTGAGC | TTGACAGAAGAGAGTGAGCAC |
| mit-miR319c | CTTGGACTGAAGGGAGCT | CTTGGACTGAAGGGAGCTCCC |
| mit-miR167b-5p | GAAGCTGCCAGCATGATC | TGAAGCTGCCAGCATGATCTGA |
| mit-miR166a | TCGGACCAGGCTTCATTC | UCGGACCAGGCUUCAUUCC |
| mit-miR528-5p | TTGGAAGGGGCATGCAGA | TGGAAGGGGCATGCAGAGGAG |
| mit-miR858b | TTCGTTGTCTGTTCGACCT | TTCGTTGTCTGTTCGACCTTG |
| 内参U6_F | ACAGAGAAGATTAGCATGGCC | |
| 内参U6_R | GACCAATTCTCGATTTGTGCG | |
| qPCR下游通用引物 | TGCGTGTCGTGGAGTCG | |
| miRNA反转录引物 | GTCGTATCCAGTGCGTGTCGTGGAGTCGGCAATTGCACTGGATACGA | |
1.2.4 Real-time qPCR分析 用Real-time qPCR检测野生蕉不同颜色果皮中差异表达miRNA的相对表达量。反应体系包括5 mL Nuclease-Free Water、8 μL 2×SYBGEEN PCR mix（QIAGEN）、1 μL上游引物、1 μL下游引物、1 μL样品cDNA,构成总体积为20 μL反应体系。反应条件为95 ℃变性2 min,94 ℃ 10 s,60 ℃ 10 s,72 ℃ 40 s,50个PCR循环,采用默认设置自动生成Ct值。RT-qPCR相对定量的内参选择U6。每个样品设置3个重复。以绿色果皮的基因表达量为对照,用相对定量2-ΔΔCt法计算分析基因的表达水平。
2 结果与分析
2.1 测序数据质量评估
通过高通量测序,从GP和PP 2个小RNA文库中分别获得了14 505 778和19 523 133条原始reads,原始数据已上传至NCBI SRA数据库,登录号分别为SRR9945688和SRR9945687。去除低质量读数、未知碱基N含量≥10%读数、接头序列、poly(A)读数以及短于18 nt或长于30 nt的读段后,分别获得11 311 826和17 420 104条clean reads,占原始reads的77%以上（表2）。
表2 闽侯野生蕉2个小RNA文库的测序数据统计Tab. 2 Statistics of sequencing reads from two small RNA libraries in Minhou wild banana |
| 类型 Type | 绿色果皮 Green peel | 紫色果皮 Purple peel |
|---|---|---|
| 原始读段Raw reads | 14 590 578 | 19 523 133 |
| N含量超过10%的读段Containing ‘N’ reads | 2 760 (0.02) | 11 511 (0.06) |
| 长度小于18 nt的读段Length<18 nt | 875 514 (6.04) | 959 572 (4.92) |
| 长度大于30 nt的读段Length>30 nt | 2 400 478 (16.55) | 1 131 946 (5.80) |
| 有效读段Clean reads | 11 311 826 (77.53) | 17 420 104 (89.23) |
| 碱基质量值Q30/% | 94.39 | 95.09 |
Note: The percentages in the parentheses indicate that different types of data accounts for the raw reads. | |
注：括号内的百分数表示不同类型数据占原始数据的百分比。 |
利用Bowtie软件将clean reads与Silva、GtRNAdb、Rfam和Repbase数据库进行比对,去除了rRNA、tRNA、snRNA、snoRNA和重复序列后,获得了包括miRNA在内的未经注释reads,这部分数据分别占2个文库总小RNA原始reads的67.05%和71.52%（图2）。此外,2个数据集中rRNA的比例小于31%,表明2个sRNA文库的质量很高。未注释reads中,2个文库（GP和PP）与Musa itinerans基因组相匹配的reads分别占54.54%和55.52%。
许多研究表明,不同长度的小RNA具有不同的功能。例如,21 nt的小RNA通常调控转录后的基因沉默,而24 nt的小RNA通过RNA依赖的DNA甲基化和异染色质维持介导基因的沉默[12,13]。因此,我们调查了2个文库小RNA序列的长度分布,结果表明,2个文库小RNA序列的长度分布相似,主要在18 nt到30 nt长度之间变化（图3）。长度在21~24 nt的小RNA序列,占全部数据的60%以上。2个文库中24 nt小RNA均占优势,分别占总数据的28.08%（GP）和28.93%（PP）。典型的21 nt长度的miRNA,在2个小RNA文库中的数量位居第二,分别占14.51%（GP）和18.99%（PP）。这些结果与在拟南芥[14]、水稻[15]、普通小麦[16]、大豆[17]、苹果[18]、茄科茄属[19]和萝卜[20]等成熟植物观察到的小RNA长度分布一致。另外,在紫色果皮中,长度小于或等于24 nt的小RNA比绿色果皮中的多,特别是21 nt小RNA的丰度差异尤为显著;长度大于24 nt的小RNA相对丰度则相反,在绿色果皮中要比在紫色果皮中的多。
2.2 野生蕉果皮miRNA的鉴定
使用Ssearch36软件将与基因组匹配的读本与miRBase（V21.0）数据库中已知的植物成熟miRNA进行比对,获得闽侯野生蕉保守的miRNA序列,结果共鉴定出132个已知的miRNA序列和122个新的miRNA。对获得的254个miRNA序列进行miRNA家族分析,研究这些miRNA在进化中的保守性。结果表明,这些基因来自58个miRNA家族（图4）。这些保守的miRNA家族成员数目显著不同。例如,miR156家族成员数量最多,有21名成员,其次是miR396和miR171-1家族,分别有20名和18名成员。
表3 野生蕉果皮中表达量前10位的miRNATab. 3 The most abundantly expressed miRNA in wild banana peel (Top 10) |
| miRNA名称 miRNA name | 表达量TPM | 成熟序列 Mature sequences (5°-3°) | |
|---|---|---|---|
| GP | PP | ||
| mit-miR319a | 561 988.90 | 594 844.20 | TTGGACTGAAGGGAGCTCCCT |
| mit-miR159a | 142 135.00 | 113 958.90 | TTTGGATTGAAGGGAGCTCTA |
| mit-miR166a-3p | 93 678.99 | 92 398.41 | TCGGACCAGGCTTCATTCCCC |
| mit-miR319a-3p | 36 221.95 | 27 648.01 | TTGGACTGAAGGGAGCTCCC |
| mit-miR319g | 25 399.34 | 34 598.86 | TTGGACTGAAGGGAGCTCCCA |
| mit-miR162 | 21 767.16 | 21 216.95 | TCGATAAACCTCTGCATCCGG |
| mit-miR166m | 15 158.68 | 14 156.14 | TCGGACCAGGCTTCATTCCCT |
| mit-miR319e | 12 287.68 | 9 653.40 | TTGGACTGAAGGGAGCTCCT |
| mit-miR168a-5p | 7 264.34 | 7 565.15 | TCGCTTGGTGCAGGTCGGGAA |
| mit-miR396b-5p | 8 638.88 | 4 540.91 | TTCCACAGCTTTCTTGAACTT |
由于miRNA长度影响其调控作用,我们统计了254个miRNA的长度分布范围。结果表明,所鉴定的miRNA长度在18~25 nt范围内分布,其中长度为21 nt的miRNA所占比例最高,为48.03%,另有19.29%的miRNA长度在24 nt,而18 nt的miRNA所占的比例最低,为0.39%;已知的132个miRNA长度分布相对较集中,范围在19~22 nt之间,其中20 nt和21 nt的miRNA占了85.61%;122个新的miRNA长度分布范围同总的miRNA分布范围一致,18~25 nt的miRNA都有,其中,24 nt长度的miRNA数量最多,占40.16%,其次是20 nt的miRNA,占33.61%（图6）。
2.3 不同颜色果皮差异表达的miRNA分析
为了了解果皮2种颜色差异表达的miRNA,用绿色果皮为对照,计算Log2FC,以校正后的P-value(FDR)≤0.05,|log2(FC)|≥1为标准,筛选出2种颜色果皮差异表达的miRNA,结果从新鉴定的254个野生蕉果皮miRNA中,筛选出36个差异表达的miRNA,它们来源于16个已知的miRNA家族,结果如图7所示。这36个差异表达的miRNA中,有10个在紫色果皮中上调表达,26个在绿色果皮中上调表达。
为了更直观地了解野生蕉果皮差异表达miRNA的表达情况,对筛选出的差异表达miRNA做层次聚类分析,将具有相同或相似表达行为的miRNA进行聚类,结果见图8。从图8可以看出,这36个miRNA在不同颜色的果皮中具有差异表达模式,根据表达量及差异显著性可以划分为4类：（1）平均表达量TPM<20的miRNA有9个,包括mit-miR157d-3p、mit-miR172a-3p-2、mit- miR160-5p、mit-miR482-2、mit-miR156f-5p、mit- miR156b-5p、mit-miR156、mit-miR399和mit- miRn4;（2）平均表达量20<TPM<100的miRNA有11个,包括mit-miR172a-3p-1、mit-miRn1、mit-miR156-1、mit-miR156-2、mit-miR171、mit- miR408-1、mit-miR408-2、mit-miR159、mit-miR156- 3、mit-miR7536和mit-miRn3;（3）平均表达量100<TPM<1100的miRNA有10个,包括mit- miR8275、mit-miR5272、mit-miR535、mit-miR396、mit-miR167d-5p、mit-miR482-1、mit-miR535d、mit- miR3627、mit-miR156j和mit-miR156k;（4）平均表达量TPM>1100的miRNA有6个,包括mit- miR167b-5p、mit-miR167、mit-miR396b-3p、mit- miR528-5p、mit-miR535和mit-miRn2。
2.4 差异表达miRNA靶基因的预测与注释
为了更好地了解不同颜色野生蕉果皮差异表达miRNA的功能,用Target Finder软件对在不同颜色果皮中差异表达的36个miRNA进行了靶基因预测,考察了它们的功能注释和富集情况。结果表明,36个差异表达的miRNA共预测到1164个靶基因,其中有1085个靶基因获得公共数据库的注释。其中,经过与GO数据库比对,共有482个差异表达miRNA的靶基因获得注释。这些靶基因分布在40个GO功能二级分类中,结果见图9。差异表达miRNA的靶基因功能主要涉及生物学过程分类中的代谢过程、细胞过程和单一有机体过程,细胞组分分类中的细胞部分、细胞和生物体,分子功能分类中的结合和催化活性等GO功能条目。
通过差异miRNA靶基因的KEGG注释结果统计,共有371个靶基因注释到KEGG的酶,其中有198个靶基因关联到59个KEGG代谢通路上,如图10所示。其中,涉及植物-病原互作、植物激素信号传导、核糖体和氨基酸生物合成等4个具体途径的靶基因数较多。而注释到次生代谢物合成的靶基因,主要集中在苯丙烷生物合成（Phenylpropanoid biosynthesis）途径中,包含预测出的7个靶基因。
在苯丙烷生物合成通路中（图11）,多数目标基因集中在苯丙烷生物合成的下游,其中4个基因在紫色果皮中下调表达,3个基因上调表达,且这些差异表达miRNA的靶基因主要受到miR156和miR482家族的调控。可以看到,一个靶基因受到多个miRNA的共同调控,而同一个家族的miRNA可以调控多个靶基因。
2.5 差异表达miRNA靶基因GO功能和KEGG通路富集分析
为了更好地描述不同颜色果皮中检测到的差异表达miRNA靶基因,将预测出的靶基因进行了GO功能和KEGG途径富集分析。结果表明,差异表达miRNA靶基因主要在DNA结合（DNA binding）、离子跨膜运输活动（anion transmembrane transporter activity）、高渗的响应（hyperosmotic response）、DNA代谢过程（DNA metabolic process）、蛋白激酶活性（protein kinase activity）、有机离子转运（organic anion transport）、细胞分化（cell differentiation）和蛋白磷酸化（protein phosphorylation）等8个生物学过程显著富集（表4）。
表4 差异表达miRNA靶基因GO功能富集分析（KS<0.01）Tab. 4 The GO enrichment result of differentially expressed miRNA (KS<0.01) |
| GO术语编号 GO ID | 功能术语 Term | GO 功能描述 GO Term Description | 注释到该功能的基因数 Annotated | 差异表达miRNA注释到该功能的基因数 Significance | 显著性KS |
|---|---|---|---|---|---|
| GO: 0003677 | Molecular function | DNA binding | 107 | 36 | 0.00075 |
| GO: 0008509 | Molecular function | anion transmembrane transporter activity | 22 | 15 | 0.00107 |
| GO: 0006972 | Biological process | hyperosmotic response | 10 | 8 | 0.0011 |
| GO: 0006259 | Biological process | DNA metabolic process | 46 | 13 | 0.0022 |
| GO: 0004672 | Molecular function | protein kinase activity | 129 | 47 | 0.00287 |
| GO: 0015711 | Biological process | organic anion transport | 14 | 10 | 0.0033 |
| GO: 0030154 | Biological process | cell differentiation | 43 | 20 | 0.0048 |
| GO: 0006468 | Biological process | protein phosphorylation | 99 | 40 | 0.0068 |
KEGG通路富集分析表明,差异表达miRNA靶基因主要在植物-病原互作通路富集,富集因子达1.93（图12）。
2.6 差异表达miRNA靶向的转录因子分析
鉴于植物miRNA具有倾向于靶定转录因子的特点,本文考察了差异表达miRNA的靶基因中转录因子的情况,结果找到了160个靶基因与转录因子有关,其中有65个注释为MYB/MYB- related基因,19个注释为bHLH,3个注释为WD40,12个注释为WRKY,17个注释为AP2基因,35个注释为SPL（图13）。将野生蕉果皮中靶定转录因子的miRNA进行统计,结果表明,注释为MYB转录因子的靶基因主要是由miR159、miR858和miR482等家族调控,注释为转录因子AP2的靶基因主要是由miR172和miR482家族调控,注释为SPL转录因子的靶基因主要是由miR156和miR535家族调控,注释为bZIP的靶基因主要是由miR396和miR482家族调控（表5）。
图13 差异表达miRNA靶基因中的转录因子类型Fig. 13 Transcription factor type in differentially expressed miRNA target gene |
表5 部分调控转录因子的miRNATab. 5 miRNAs that regulate transcription factors |
| miRNA 家族miRNA family | 差异表达miRNA Differently expressed miRNA | 靶基因 Target gene |
|---|---|---|
| miR482 | mit-miR482-1, mit-miR482-2 | Myb4/TT2/C1, SAP, bZIP |
| miR482 | mit-miR482-2 | MYB101, AP2, TAF |
| miR858 | mit-miR858b | MYB5, MYB61 |
| Unknown | mit-miRn4 | MYB5, WRKY, AP2 |
| Unknown | mit-miRn1 | MYB61 |
| miR159/319 | mit-miR159 | MYB101, GAMYB |
| miR156 | mit-miR156, mit-miR156f-5p, mit-miR157d-3p | bHLH |
| miR156 | mit-miR156f-5p, mit-miR157d-3p | WRKY |
| miR156 | mit-miR156b-5p, mit-miR156j, mit-miR156k, mit-miR156, mit-miR156f-5p | SPL |
| miR535 | mit-miR535, mit-miR535d | |
| miR172 | mit-miR172a-3p, mit-miR172a-3p | AP2 |
| miR535 | mit-miR535 | TAF |
| miR396 | mit-miR396b-5p | bZIP |
| miR156 | mit-miR157d-3p | SPT6 |
| miR156 | mit-miR156 | NFYA |
| miR528 | mit-miR528-5p | HSF/CHI |
2.7 差异表达miRNA定量PCR验证
为了验证RNA-seq获得的miRNA表达模式,选取了6个在野生蕉不同颜色果皮中差异表达的miRNA（mit-miR156k、mit-miR319c、mit-miR167b-5p、mit-miR166a、mit-miR528-5p、mit-miR858b）进行RT-qPCR,将2种分析结果（RNA-seq和RT-qPCR）进行比较,进行定量验证（图14）。结果表明,供试的6个基因中,mit-miR156k、mit-miR166a、mit-miR528-5p在绿色果皮中上调表达,而mit-miR319c、mit-miR167b-5p、mit-miR858b 3个miRNA在紫色果皮中上调表达,与高通量测序的结果一致。
3 讨论
3.1 野生蕉果皮miRNA具有品种特异性
Ghag等[21]在2个栽培香蕉品种Grand Naine（AAA subgroup Cavendish）和Rasthali（AAB subgroup Silk）的研究中,发现了miRNA具有明显的品种特异性,2个品种分别具有44%~53%的特异性miRNA。宋顺等[22]在3个植物保守miRNAs响应香蕉不同品种（巴西蕉、大蕉、皇帝蕉）的抗枯萎病特性表达分析中也发现miRNA具有明显的品种特异性。本研究利用RNA-seq测序技术,共从野生蕉果皮中鉴定了254个miRNA,其中132个miRNA具有保守miRNA序列,比对上已知的其他物种miRNA,而有122个miRNA为野生蕉果皮特异的,占总数的48.03%。此外,本研究获得的miRNA中,保守的已知miRNA以21 nt长度为主,野生蕉特异的新的miRNA以24 nt为主,21 nt次之,这与2个栽培香蕉品种Grand Naine (AAA subgroup Cavendish)和Rasthali (AAB subgroup Silk)的miRNA长度也有所不同,它们24nt miRNA数量远大于其他长度的miRNA数量[21],但却与香蕉Musa acuminata Colla（AAA Group）cv. Berangan 根中的miRNA长度分布一致[23]。这些都证实了芭蕉属植物miRNA具有明显的品种特异性。
3.2 可能参与野生蕉果皮花青素合成的miRNAs
通过对差异表达miRNA的靶基因GO功能和KEGG途径富集分析显示,靶基因主要在DNA结合和离子跨膜运输活动等功能和植物病原互作等途径富集,也有一些靶基因富集在苯丙烷代谢途径中,有些靶基因可能是参与花青素合成的MYB、bHLH、WD40、WRKY、SPL、AP2等转录因子,暗示野生蕉果皮颜色调控的分子机制复杂而多样。
在不同颜色果皮中我们共获得36个差异表达的miRNA,它们多数（26个miRNA）在紫色果皮中呈下调表达,而上调表达的基因主要来自miR160、miR167-1、miR171-1、miR399、miR482和miR5272家族的10个成员。已有研究报道Sly-miR167在番茄果实破色期的丰度增强[24],在我们的研究中预测到的miR167-1（mit-miR167b-5p、mit-miR167）在紫色果皮中表达量TPM大于3000,这2个miRNA可能与Sly-miR167一样,在果皮花青素的合成中起调控作用,但这些都还需进一步通过试验验证。
另外,通过差异表达miRNA靶基因的预测也发现,野生蕉果皮中表达的miRNA159/miRNA319和miRNA858b靶定的基因,同陆地棉等植物一样[25],也多数是编码MYB转录因子的基因,可能参与植物花青素的调控。许多证据证明,miRNA858调控与花青素相关的MYB基因的表达,如番茄中miR858能在转录后水平通过对靶基因Sly-MYB12和Sly-myb-like的作用负调控花青素的生物合成。通过与转录组数据[26]的关联分析也发现野生蕉果皮靶向MYB的miRNA也有miRNA858b,而且它在紫色果皮中上调表达,但差异未能达到显著（log2FC=0.96）,它是否在野生蕉果皮颜色调控中起作用还需更进一步的研究。此外,研究中还发现2个属于miR482家族的miRNA（mit-miR482-1、mit-miR482-2）靶向调控MYB基因,但目前研究表明,miR482家族主要靶向NBS-LRR病害抗性基因,而对MYB基因的调控尚未见报道。总而言之,以上这些预测均需进一步的研究才能证实。