羽衣甘蓝红色系叶片形成相关microRNAs的鉴定及分析

张新鹏 徐宗大 刘林馨 夏甜甜 李晓艳*

(1.山东建筑大学 风景园林科学研究中心，济南 250101； 2.山东农业大学 林学院，山东 泰安 271018； 3.山东建筑大学 建筑城规学院，济南 250101)

羽衣甘蓝(

Brassica

oleracea

var.

acephala

)为十字花科芸薹属食用甘蓝的园艺变种，为二年生草本观叶植物。羽衣甘蓝叶片，特别是其心叶的叶色表型艳丽多彩，株形如牡丹盛开的花朵，业界对其有“叶牡丹”的美誉”。同时羽衣甘蓝耐寒性极强，观赏期较长，对我国北方地区景观单调和少花的冬季而言，该植物对园林景观绿化和美化的意义重大。羽衣甘蓝心叶具体有纯白、乳白、淡黄、浅粉、粉红、玫红和紫红等类型，其中当属红色系叶片观赏色彩最为突出，而且其花青苷含量相对最高，对人具有很好的保健功能。因此，探索和揭示红色系叶片形成的调控机理，对于培育色彩突出且营养价值含量高的羽衣甘蓝新品种具有重要意义。在许多园林植物中，叶色是其非常重要的观赏性状之一。对色彩鲜艳的叶片(如红色和紫色等)而言，在生理层面，花青苷往往是呈色的关键花色素。先前许多研究表明花青苷合成结构基因的差异表达通常是影响植物花青苷的差异积累最直接的因素，而这些结构基因的表达模式和强度往往会受到其上游转录因子在不同时空上的调控。另外，在转录后水平上，microRNA(miRNA)也是植物中一类非常重要的基因表达调控因子，它是一类长度约21～24 nt的非编码单链小分子RNA。近年来大量研究发现，许多miRNAs(如miR156、miR828和miR858等)可通过靶向某些与花青苷生物合成相关的转录因子(如

MYB

和

SPL

等)并抑制其表达，进而影响下游花青苷合成关键结构基因的表达水平，并最终调控植物组织中花青苷的合成和积累。

近年来，在许多植物中高通量测序技术已被广泛应用于鉴定花青苷生物合成相关的重要miRNAs，如美人蕉、芍药、蝴蝶兰、牡丹和茶树等。然而，对于羽衣甘蓝而言，目前尚无有关叶色形成相关miRNA的报道。为了探索羽衣甘蓝红色系叶片中在转录后水平miRNAs对花青苷的调控机制，本研究以白色系叶片作为对比试材，以红色系的羽衣甘蓝叶片作为主要研究试材，基于高通量测序分析，以期筛选到红色系叶片与花青苷合成相关的关键miRNAs，为今后在miRNA-mRNA层面系统探索羽衣甘蓝叶色的分子调控机制提供一定的理论基础。

1 材料与方法

1.1 试验材料

红色系和白色系的羽衣甘蓝品种‘紫鸽’‘白鸽’材料于2020年种植于山东建筑大学风景园林实训基地室外试验田的同一区域，日常培养和外界环境条件均保持一致。在经冬季低温叶片转色后的观赏期，采集这2个品种的心部(紫红色和白色)叶片，并去掉略带绿色的边缘区域，分别用编号Bo-P表示紫红色叶片，Bo-W表示白色叶片(图1)。采集的上述组织材料，经液氮速冻后保存于-80 ℃超低温冰箱。

图1 研究用紫红色(Bo-P)和白色(Bo-W)叶片Fig.1 Purple-red (Bo-P) and white (Bo-W) leaves used in this study

1.2 RNA提取和小RNA测序

将紫红色和白色叶片组织样品在液氮冷冻下进行研磨后，使用Trizol法提取总RNA，每种颜色叶片进行3次生物学重复。检测合格的RNA样品再进行进一步的文库构建和小RNA测序。本研究构建了6个小RNA(sRNA)文库，紫红色叶片的文库分别命名为Bo-P1、Bo-P2和Bo-P3，白色叶片的分别命名为Bo-W1、Bo-W2和Bo-W3。文库构建和高通量测序均由生工生物工程(上海)股份有限公司完成。

1.3 miRNA的鉴定和靶基因预测

对小RNA测序得到的raw reads进行数据过滤，得到高质量的clean reads。对得到的clean reads与甘蓝基因组(https:∥www.ncbi.nlm.nih.gov/genome/?term=Brassica+oleracea)、miRbase(https:∥www.mirbase.org/)和Rfam(http:∥rfam.xfam.org/)等已知的sRNA数据库进行比对，对sRNA分类注释。其中，通过与miRbase 数据库(版本号22.1)中登录的已知植物的miRNA成熟体序列进行比对，来鉴定羽衣甘蓝叶片中的保守miRNAs(known miRNAs)。同时，基于miRNA的前体能够形成发夹二级结构，用miRdeep2软件再对未注释的clean reads进行miRNA的预测，鉴定新miRNAs(novel miRNAs)。

以甘蓝基因组序列作为参考序列，用TargetFinder软件对羽衣甘蓝叶片中鉴定出的保守和新的miRNA进行靶基因的预测。

1.4 miRNA的差异表达分析

用R语言包DESeq2进行miRNA的差异表达分析。本研究以白色叶片样本作为对照，采用

P

value≤0.05和|log(Fold change)|≥1这2个阈值作为筛选紫红色和白色叶片这2组样本中差异表达miRNA的条件。同时，对差异表达miRNA预测的靶基因进行GO和KEGG pathway富集分析。

1.5 候选miRNA的荧光定量(qRT-PCR)分析

利用TaKaRa公司的Mir-XmiRNA First-Strand Synthesis Kit反转录试剂盒对miRNA进行反转录，得到的cDNA用于候选miRNAs的进一步定量。再利用TaKaRa 公司的TB GreenPremix Ex

Taq

II试剂盒，通过qRT-PCR法对候选miRNAs(miR828和miR858-3)在紫红色和白色叶片中的差异表达情况进行验证分析。其中，内参基因为

U

6，每个样品设置3次重复。以白色叶片的基因表达量为对照，通过2相对定量法计算miRNA的相对表达水平。引物合成由生工生物工程(上海)股份有限公司完成，引物的具体序列信息见表1。

表1 候选miRNA荧光定量PCR所用的引物序列
Table 1 Primer sequences of the candidate miRNAs used for qRT-PCR

miRNA名称miRNA name上游引物(5'-3')Forward primer (5'-3')下游引物(5'-3')Reverse primer (5'-3')miR828ACACTCCAGCTGGGTCTTGCTTAAAT-GAGTGGTGTCGTGGAGTCGmiR858-3ACACTCCAGCTGGGTTCGTTGTCTGTTCTGGTGTCGTGGAGTCGU6TTGGAACGATACAGAGAAGATTAGCATTGGACCATTTCTCGATTTGTG

2 结果与分析

2.1 测序数据质量评估

本研究利用高通量测序法对紫红色和白色叶片的6个小RNA文库(Bo-P1、Bo-P2、Bo-P3、Bo-W1、Bo-W2和Bo-W3)进行了测序，从紫红色和白色叶片文库获得的原始数据量分别为61 388 537和52 363 761 条raw reads。通过去除包含poly-A、带接头和低质量等reads，从紫红色叶片文库中最终获得了35 183 047条Q20为97.42%的clean reads，从白色叶片文库中获得了33 023 141条Q20为97.49% 的clean reads(表2)。

对于sRNA的长度分布，在紫红色和白色叶片这2个文库中，大多数sRNA的长度分布区间均为21～24 nt。在Bo-P和Bo-W文库中，长度为24 nt的sRNA所占比例均为最高，分别为16.4%和17.6%，其他主要长度所占比例的排序依次为21、22和23 nt(图2)。

2.2 羽衣甘蓝叶片中miRNA的鉴定及靶基因的预测

小RNA(sRNA)的组成种类比较多，通过将紫红色和白色叶片这2个文库得到的clean reads分别与甘蓝基因组和不同的sRNA数据库进行比对，可以将miRNA、rRNA、tRNA、snRNA和snoRNA等分类注释出来。结果显示，在Bo-P(图3(a))和Bo-W(图3(b))这2个文库中分别有2.48%和1.99% 的clean reads比对到miRNA上。经进一步数据分析，在这6个sRNA文库中共鉴定出1 027个保守miRNAs(known miRNAs)。在Bo-P和Bo-W文库中分别鉴定出799和835个known miRNAs，其中这2个文库中的大多数known miRNAs(607个)都是相同的，另外Bo-P和Bo-W文库分别各有192和228个不同的known miRNAs。初步推测紫红色或白色叶片中这些不同的known miRNAs可能与羽衣甘蓝叶色的形成有一定的关联。

表2 6个小RNA文库测序数据
Table 2 Sequencing data of small RNAs (sRNAs) in the six libraries

样品名称Sample name过滤后的readsClean readsQ20值/%Q20 valueGC含量/%GC content平均长度/bpAverage lengthBo-P111 610 46097.5544.2222.30Bo-P211 980 71797.1943.9922.63Bo-P311 591 87097.5144.0822.64Bo-W111 443 05897.4343.5723.00Bo-W210 014 61697.4843.6122.97Bo-W311 565 46797.5743.7822.95Bo-P总计 Bo-P summary35 183 04797.4244.1022.52Bo-W总计 Bo-W summary33 023 14197.4943.6522.97

图2 sRNA在紫红色叶片(a)和白色叶片(b)中的长度分布Fig.2 Length distribution of sRNA in Bo-P leaf (a) and Bo-W leaf (b)

对于那些与已知sRNA数据库未比对上的clean reads，对其进行了新miRNA(novel miRNA)的预测。结果发现，在这6个sRNA文库中共鉴定出100个novel miRNAs。其中，在Bo-P和Bo-W文库中分别鉴定出100和99个novel miRNAs。另外，对于在紫红色和白色叶片的6个sRNA文库中鉴定出的known miRNAs和novel miRNAs，通过分析发现它们分别有892和86个miRNAs存在靶基因，其相应的靶基因个数为8 428和1 923个。

括号内的数字和百分号数字分别代表该项的reads数和其所占clean reads的比例。 The numbers and percentages in brackets represent the number of reads of this item and its proportion in clean reads, respectively.图3 sRNA在紫红色叶片(a)和白色叶片(b)中的分类注释Fig.3 Classification of sRNA annotations in Bo-P leaf (a) and Bo-W leaf (b)

2.3 不同颜色叶片差异表达的miRNA分析

为了解影响羽衣甘蓝红色系叶片形成的关键miRNA，对紫红色和白色叶片中鉴定出的miRNAs进行了差异表达分析。基于阈值

P

value≤0.05和|log(Fold change)|≥1，在这2种颜色叶片中共鉴定出110个(属于33个家族)差异表达的miRNAs(DEMs)，其中32个在紫红色叶片中表达上调，78个在紫红色叶片中表达下调(图4)。

2.4 差异表达miRNA靶基因的功能富集分析

为了更好地明确紫红色和白色叶片中这些靶基因行使的主要生物学功能，对其进行了GO和KEGG功能富集分析。GO功能显著富集(矫正后的

P

value≤0.05)分析结果(图5)显示，对于生物学过程类别而言，富集主要集中于转录(Transcription)和转录调控(Regulation of transcription)这2个子类别。对于细胞成分类别而言，富集主要集中于CCAAT结合位点复合体(CCAAT-binding factor complex)和细胞核(Nucleus)这2个子类别。在分子功能这1大类中，富集主要集中于转录因子活性及序列特异性DNA结合(Transcription factor activity, sequence-specific DNA binding)和DNA结合(DNA binding)这2个子类别。

对差异表达miRNA的靶基因也进行了KEGG pathway富集分析，结果显示这些差异表达miRNA的靶基因主要在植物MAPK信号通路(MAPK signaling pathway-plant)、植物激素信号转导(Plant hormone signal transduction)和类黄酮合成(Flavonoid biosynthesis)等通路富集(图6)。

图4 差异表达miRNAs的火山图Fig.4 Volcanic map of differentially expressed miRNAs

图5 差异表达miRNAs的靶基因的GO功能富集Fig.5 GO functional enrichment of target genes of differentially expressed miRNAs

图6 差异表达miRNA的靶基因的KEGG pathway富集Fig.6 KEGG pathway enrichment of target genes of differentially expressed miRNAs (DEMs)

2.5 羽衣甘蓝红色系叶片形成相关的差异表达miRNA

在这110个差异表达miRNAs中，基于靶基因的功能注释和前人花青苷合成相关miRNA的研究，筛选出了2个可能与羽衣甘蓝红色系叶片形成相关的重要miRNAs(miR828和miR858-3)。其中，miR828的靶基因为XM_013735358.1 [R2R3-MYB: Transcription factor MYB114-like]。miR858-3的靶基因包括XM_013758669.1 [Transcription factor MYB12-like]和XM_013763976.1 [R2R3-MYB: Transcription factor MYB86-like]。小RNA测序结果显示，与白色叶片相比，miR828和miR858-3在紫红色叶片中的表达均明显上调(表3)。

表3 2个可能参与羽衣甘蓝红色系叶片形成的候选差异表达miRNAs的表达水平
Table 3 Expression levels of two candidate differentially expressed miRNAs possibly involved in formation of red-color leaf in var.

miRNA名称miRNA namelog2(Bo-P/Bo-W)靶基因Target gene注释AnnotationmiR8286.06XM_013735358.1Transcription factor MYB114-likemiR858-31.85XM_013758669.1Transcription factor MYB12-likeXM_013763976.1Transcription factor MYB86-like

本研究进一步采用qRT-PCR的方法对这2个候选重要miRNAs的表达模式进行了验证分析。结果显示，它们均与小RNA测序的结果呈现出一致的趋势，其中高通量测序和qRT-PCR这2次结果显示miR828在紫红色叶片中的表达量分别为白色叶片中的66.7倍和44.4倍(图7)。

不同字母代表在0.05水平上存在显著性差异。 Different lowercase letters indicate significant differences at 0.05 level. 图7 2个候选差异表达miRNAs的定量分析Fig.7 qRT-PCR analysis of 2 candidate differentially expressed miRNAs

3 讨论与结论

羽衣甘蓝红色系内叶的色彩表型极为突出，该色系的叶片富含花青素且可食用，这使其在一定程度上兼具观赏、营养和商业价值。因此，探索其红色系叶片形成的机理具有重要意义。对于羽衣甘蓝红色系在生理层面的形成机理，前人研究发现花青苷，特别是矢车菊素类糖苷，是红色系叶片呈色的关键花色素。羽衣甘蓝红色系叶片在分子层面的形成机理，目前多数研究主要集中在mRNA(结构基因和转录因子)层面。对花青苷合成相关结构基因而言，研究发现

BoDFR

应该是羽衣甘蓝红色系叶片形成的关键结构基因。任杰通过对

DFR

启动子区域序列分析发现，该基因上存在花青苷合成相关转录因子MYB的结合位点。对花青苷合成相关转录因子而言，研究发现MYB类转录因子(如BoMYB2、BoMYB1、BoPAP1和BoMYBL2.1)可能对羽衣甘蓝红色系叶片中的花青苷积累起着重要的调控作用。但在转录后水平上，特别是miRNA层面目前尚未见报道，本研究首次报道了羽衣甘蓝红色系叶片中与叶色调控相关的miRNAs。

近年来在植物花青苷生物合成领域(包括叶色、花色和果色等)，miRNAs作为在转录后层面重要的基因调控者，引起了众多研究者的关注。本研究分别构建了‘紫鸽’和‘白鸽’羽衣甘蓝品种红色和白色系叶片的小RNA文库，并进行了高通量测序和相关分析。在这2种颜色叶片的文库中，均显示21～24 nt长度的sRNAs占多数，而24 nt的sRNA比例最高。这一结果与枸杞、甘薯和茶树等的研究结果一致，而与大丽花、芍药和牡丹等的报道不同，这表明在不同植物中sRNA的生物合成途径可能存在一些差异。

本研究在紫红色和白色叶片文库中共鉴定了110个显著差异表达的miRNAs，同时对其进行了靶基因的预测及功能分类。基于前人花青苷合成相关miRNAs的研究和对相应靶基因的进一步分析，本研究在Bo-P和Bo-W 文库中共筛选到2个与花青苷合成相关的DEMs(miR828和miR858-3)，并发现它们的靶基因中均有R2R3-MYB转录因子(XM_013735358.1和XM_013763976.1)。对MYB类转录因子家族而言，许多研究已表明其参与调控花青苷生物合成的亚族多为R2R3-MYB转录因子。

对于miR828和miR858这2个家族的一些miRNAs，前人许多研究也已表明它们均可以通过调控其靶基因R2R3-MYB转录因子的表达，进而最终影响花青苷的生物合成和积累。本研究的qRT-PCR分析结果显示，miR828和miR858-3的表达趋势与高通量测序结果基本一致。同时这2个结果均显示，miR828在紫红色和白色叶片中的差异高达44倍以上。

因此，本研究认为miR828和miR858-3(特别是miR828)应该是转录后水平上调控羽衣甘蓝红色系叶片形成的关键miRNAs。该结果将有助于人们进一步理解miRNAs在羽衣甘蓝叶色形成中潜在的功能角色，以期为今后深入揭示羽衣甘蓝红色系叶片形成的分子机理提供一定程度的理论参考。