基于高通量测序筛选‘紫娟’花青素合成相关的miRNA

陈林波，夏丽飞，刘悦，孙云南，蒋会兵，田易萍，陈亮

基于高通量测序筛选‘紫娟’花青素合成相关的miRNA

陈林波1,2，夏丽飞1，刘悦1，孙云南1，蒋会兵1，田易萍1，陈亮2*

1. 云南省农业科学院茶叶研究所/云南省茶树种质资源创新与配套栽培技术工程研究中心/云南省茶学重点实验室，云南 勐海 666201；2. 中国农业科学院茶叶研究所，浙江 杭州 310008

为筛选调控紫娟茶树花青素合成相关miRNA，以茶树品种紫娟（ZJ）、云抗10号（YK）和福鼎大白茶（FD）为材料，构建了miRNA文库。鉴定出46种已知的miRNAs和67种新的miRNAs，预测到具有注释的靶基因765个。通过差异表达分析，筛选出在ZJ与YK、ZJ与FD共有差异表达的miRNA 24个。通过对24个差异表达miRNA的靶基因分析，筛选出可能参与调控花青素合成的miRNA 4个，包括miR828a、miR845c、novel_14和novel_87，其预测的靶基因包括转录因子基因、、、、以及4-香豆酰辅酶A链接酶（）、二氢黄酮醇4-还原酶（）和UDP-葡萄糖黄酮3--葡糖基转移酶（）等基因。利用RT-PCR分析8个差异表达的miRNA，其结果与转录组分析一致。本研究结果为进一步开展茶树花青素生物合成的调控机制研究奠定基础。

茶树；花青素合成；高通量测序；miRNA；靶基因

MicroRNA又称miRNA，是广泛存在于生物体内的一种长度约21~25个核苷酸组成的内源非编码小分子RNA[1]。miRNA通过与靶mRNA互补配对来识别靶基因mRNA，在转录水平上降解靶mRNA或抑制靶mRNA的翻译，从而调节靶mRNA的丰度和功能[2]。研究表明，miRNA对于植物花青素的生物合成具有一定的调控作用。在拟南芥中，miR156通过靶向SPL转录因子对花青素的生物合成起着负调控作用[3-4]；miR828通过介导调控靶基因，在缺磷条件下对、和的表达进行调控，进而影响花青素的积累[5]。过量表达miR778的拟南芥在缺磷条件下会促进花青素的积累[6]。在番茄中，miR858可通过靶向R2R3 MYB转录因子参与花青素合成的调控[7]。

目前关于茶树花青素代谢的研究主要集中在关键酶和转录因子方面，例如MYB[8]、bHLH[9]、花青素合成酶[10-11]。然而，有关茶树中miRNA介导的茶树花青素生物合成的调节机制研究较少。紫娟是一种特异的茶树品种，其幼嫩新梢芽、叶、茎均为紫色，富含花青素。花青素是紫娟茶树叶片呈现紫色的特征成分，是区别于其他茶树的主要生化物质，也是决定紫娟茶健康功效的主要成分[12-13]。此外，紫娟茶还是提取天然花青素的重要来源物[14-15]。因此，我们利用Illumina HiSeq 2500测序平台和生物信息学技术对紫叶茶树和绿叶茶树中的花青素代谢相关的miRNA以及靶基因进行比较分析，为进一步研究miRNA对茶树花青素合成的调控机制提供理论基础。

1 材料与方法

1.1 试验材料

分别采摘云南省农业科学院茶叶研究所品种园内树龄均为10年的紫叶茶树品种紫娟（ZJ）、绿叶茶树品种云抗10号（YK）和福鼎大白茶（FD）的一芽二叶，其中紫娟和云抗10号为阿萨姆茶（var.），福鼎大白茶为茶种（var.），液氮速冻后置于–80℃冰箱保存备用。

1.2 试验方法

1.2.1 RNA的提取与质量评估

RNA的提取采用美国Invitrogen公司的Trizol试剂盒，分别提取ZJ、YK、FD的总RNA，采用1%的琼脂糖凝胶电泳分析RNA降解程度以及是否有污染，利用Nanodrop 2000微量分光光度计检测RNA的纯度，利用Qubit 4.0（美国）荧光定量仪对RNA浓度进行定量，利用Agilent 2100精确检测RNA的完整性。

1.2.2 sRNA文库的构建和测序

利用Small RNA Sample Pre Kit分别构建ZJ、YK、FD小RNA的cDNA文库。采用Qubit 4.0荧光定量仪进行初步定量，使文库浓度为1 ng·μL-1，再利用Agilent 2100对文库的质量和Insert size进行检测，测序与分析由北京诺禾致源生物有限公司完成。

1.2.3 sRNA测序数据分析

高通量测序得到的原始图像数据文件经碱基识别（Base calling）分析转化为Raw reads。再将Raw reads带接头、低质量的reads进行处理，得到Clean reads。筛选长度范围为18~30 nt的sRNA进行分析。用bowtie将长度为18~30 nt的sRNA定位到参考序列上，分析small RNA在参考序列上的分布情况。将mapped sRNA与miRBase中指定范围序列进行比对，得到各样品匹配上的sRNA详细情况。选取Rfam中的rRNA、tRNA、snRNA、snoRNA来注释测序所得的小RNA。通过检索miRBase21.0数据库（http://www.mirbase.org）来确定已知的miRNA，采用miRNA前体的标志性发夹结构来预测新的miRNA[16-17]。

1.2.4 茶树miRNA靶基因预测与功能分析

分析得到的已知miRNA和新的miRNA，利用psRobot软件[18]对前期获得的ZJ、YK和FD的转录组数据[19]进行靶基因预测，根据miRNA与其靶基因间的对应关系获得miRNA靶基因。靶基因通过在线网站（http://www.geneontology.org）进行Gene ontology分析和Kyoto Encyclopedia of Genes and Genomes数据库（https://www.genome.jp）进行Pathway显著性富集分析。

1.2.5 茶树叶色相关miRNA筛选

对3个样本中已知的和新的miRNA进行表达量的统计，采用TPM对miRNA表达水平的readcount数据进行标准化处理，再用DEGseq[20]进行差异分析。差异miRNA筛选条件为：qvalue<0.01 & |log2(foldchange)|>1。

1.2.6 实时荧光定量PCR的检测

表1 实时定量PCR所使用引物

2 结果与分析

2.1 茶树叶片总RNA的提取

采用Trizol试剂方法提取ZJ、YK和FD中的总RNA，经凝胶电泳（图1）和紫外分光光度计检测，其总RNA浓度分别为709、482、559 ng·μL-1，OD260/OD230值分别为2.05、2.00和1.90，OD260/OD280值分别为2.12、2.11、2.11，28 S和18 S条带亮度等均符合转录组分析的要求。

2.2 茶树叶片miRNA测序分析

利用高通量Illumina测序技术分别对ZJ、YK和FD的cDNA文库进行测序，通过去除接头、低质量reads、polyA/T/G/C以及N，分别获得12 361 695（98.16%）、13 930 750（98.20%）、13 528 588（98.37%）条clean reads。筛选长度在18~30 nt范围的sRNA分别有10 702 304（54.58%）、12 061 820（54.34%）、11 202 955（43.15%）条clean reads（表2）。

3个文库中sRNA主要分布在21~24 nt范围内，其中24 nt的最多（图2），这些sRNA分布情况与其他报道的结果一致[22-23]。再利用bowtie将筛选后的sRNA定位到参考序列上，分析small RNA在参考序列上的分布情况。ZJ、YK和FD的clean reads中能mapped到参考序列的分别有5 774 584、6 835 843、6 278 768条。将mapped的sRNA通过NCBI数据库（http://www.ncbi.nlm.nih.gov）和Rfam数据库进行注释。被注释到的非编码RNAs，包括已知miRNA、新miRNA、rRNAs、tRNAs、snRNAs、snoRNAs、TAS基因等（表3）。

注：M为DL5000 DNA Marker；1为紫娟；2为云抗10号；3为福鼎大白茶

表2 3个miRNA文库测序数据

表3 被注释到的各种miRNA数

注：ZJ：紫鹃，YK：云抗10号，FD：福鼎大白茶。下同

Note: ZJ: Zijuan. YK: Yunkang10. FD: Fuding dabaicha. The same as follow

2.3 茶树已知miRNA的鉴定

为了鉴定茶树已知miRNA，将3个文库中mapped到参考序列上的reads利用miRBase 21.0数据库中葡萄（）的miRNA数据进行比对，共鉴定出46条已知miRNA，分为26个家族，其中miR171、miR172和miR399分别有4个成员，miR167和miR319分别有3个成员，miR156、miR160、miR169、miR394、miR396和miR398分别有2个成员，其他均为1个成员。获得已知的miRNA的长度为19~23 nt，其中大多数miRNA的长度为21 nt（表4）。

图2 miRNA的长度分布

表4 已知的miRNA家族成员

续表4

2.4 茶树新miRNA的预测

由于茶树缺乏RNA背景，只能将独特的小RNA序列映射到茶树转录识别潜在的新miRNA序列。利用miRNA前体的标志性发夹结构预测新的miRNA[16-17]。通过BLASTn进行搜索和发夹结构预测，共发现了新的成熟miRNAs 67条。这些新miRNA序列长度为19~25 nt（表5），其中24 nt最多，占所有新预测miRNA的52.24%。这些没有被归类到已知miRNA家族的新miRNA，可能是茶树miRNA新家族成员。

2.5 茶树叶片miRNA靶基因预测

为了研究miRNAs对基因表达的调控作用，利用ZJ、YK和FD的转录组数据为参考[19]，采用psRobot软件[18]对分析得到的已知miRNA和新miRNA进行靶基因预测。113个miRNA共预测到2 200个基因，其中765个基因具有注释。这些被注释到的基因包括、、、、、和等转录因子家族。结果发现，单个miRNA可同时调节多个靶基因，而单靶基因也可由多个miRNA调节，说明茶树中miRNA具有多重靶向性。这些预测到的miRNA靶基因的准确性需要在未来的研究中进行验证。

2.6 差异miRNA的筛选

通过高通量测序产生的序列可通过在文库中的读取数来计算miRNA的丰度。ZJ、YK和FD中已知和新miRNA的读取数通过TPM进行标准化处理后，采用DEGseq进行差异分析[21]，差异miRNA筛选条件为|log2(foldchange)|>1和qvalue<0.01。筛选出在ZJ与YK、ZJ与FD之间共有的差异表达miRNA 24个（表6），在ZJ中上调表达的miRNA有13个，包括已知miRNA（miR845c、miR828a、miR482、miR399e、miR398a、miR398b、miR169a和miR162）和新miRNA（novel_93、novel_108、novel_36和novel_35）；在ZJ中下调表达的miRNA有9个，均为新miRNA（novel_87、novel_86、novel_84、novel_80、novel_46、novel_42、novel_37、novel_18和novel_14）。

2.7 差异miRNA靶基因分析

对获得的24个差异miRNA的靶基因进行分析，发现具有注释的靶基因227个，其中可能参与花青素合成调控相关的靶基因包括转录因子基因、、、和，参与花青素合成相关的结构基因包括4-香豆酰辅酶A链接酶5（）、UDP-葡萄糖黄酮3--葡糖基转移酶3（）和UDP-葡萄糖黄酮3--葡糖基转移酶6（）等（表7）。因此，推测在紫娟茶树中差异表达的miR828a、miR845c、novel_14和novel_87可能参与花青素生物合成的调控，其成熟miRNA见图3。这些预测到的参与花青素合成调控相关的miRNA需要在未来的研究中进行验证。

表5 新的miRNA

表6 差异表达的miRNA

表7 miRNA部分靶基因

2.8 miRNA的qRT-PCR验证

随机选取8个miRNA，利用qRT-PCR检测其在ZJ、YK、FD等茶树叶片中的表达情况（图4）。发现miR398a、miR398b、miR399e、miR828a在ZJ中的表达量高于YK和FD，novel_87、novel_84、novel_37和novel_14在ZJ中的表达量均低于YK和FD。8个差异表达miRNA在ZJ、YK和FD的相对表达量与测序分析的变化趋势一致，说明高通量测序技术分析结果可靠。

3 讨论

随着对植物中miRNA研究的不断深入，它已被证实广泛作用于植物的生长发育[24-26]、次生代谢[27-29]、器官形态建成[30]等过程。近年来人们发现在转录后水平上miRNA对植物花青素合成具有不可替代的调控作用[3,5,31]。为了研究miRNA对茶树花青素合成的调控机制，本研究构建了紫娟（ZJ）、云抗10号（YK）、福鼎大白茶（FD）的miRNA文库并进行测序，分别获得12 361 695、13 930 750、13 528 588条clean reads。miRNA读数主要是21~24 nt，这与大多数被子植物中相应的读数相似；以24 nt的数量最多，与烟草[32]、芦笋[22]、杜仲[33]等作物一样。

图3 关键miRNA的前体茎环结构

图4 miRNA表达的qRT-PCR分析

在本研究获得的miRNA以及筛选在紫叶与绿叶茶树差异表达的miRNA中，新预测的miRNA数量均远高于已知miRNAs，表明本研究中新发现的miRNAs对紫娟茶树叶片呈色和花青素合成贡献较大。在拟南芥中miR828介导调控的靶基因与花青素合成有关，miR828能够诱发转录本的剪切，进而调控下游靶基因、和转录因子以调控花青素的合成[31]。本研究预测到miR828的靶向基因为、、。在拟南芥中超表达miR828能阻碍转录因子[31]，推测参与了花青素的合成。本研究中还发现了新的novel_87靶向UDP-葡萄糖黄酮3--葡糖基转移酶基因（）、4-香豆酰CoA连接酶基因（）和转录因子基因，novel_14靶向二氢黄酮醇4-还原酶基因（）。本课题组前期研究紫叶茶树和绿叶茶树叶片转录组差异，筛选花青素合成相关基因时，发现了4-香豆酰CoA连接酶（）、二氢黄酮醇4-还原酶（）、类黄酮糖基转移酶（）在紫叶茶树叶片中为上调表达[19]，这与本研究获得的miRNA存在相反的表达模式。推测新预测的novel_87和novel_14可能参与了紫娟茶树叶片呈色和花青素合成的调控。此外，预测到miR845c靶向和，MYB是重要的一类转录因子，可以直接或与bHLH和WD形成MYB-bHLH-WD复合物参与花青素合成[8-9]。在本研究中没有发现miR858[5]，这可能是由于不同植物具有不同调控机制所导致。因此，本研究获得的miRNA中一部分参与调控茶树花青素的合成，但对于那些功能未知的miRNA与茶树花青素合成的关系仍需进一步研究证实和确定。

[1] Sun F, Guo G, Du J, et al. Whole-genome discovery of miRNAs and their targets in wheat (L.) [J]. BMC Plant Biology, 2014, 14(1): 142. DOI: 10.1186/1471-2229-14-142.

[2] Baulcombe D. RNA silencing in plants [J]. Nature, 2004, 431: 356-363.

[3] Gou J Y, Felippes F F, Liu C J, et al. Negative regulation of anthocyanin biosynthesis inby a miR156-Targeted SPL transcription factor [J]. The Plant Cell, 2011, 23(4): 1512-1522.

[4] Ang G, Jun Y, Gu Y, et al. Effective small RNA destruction by the expression of a short tandem target mimic in[J]. The Plant Cell, 2012, 24(2): 415-427.

[5] Hsteh L C, LIN C I, Shih C C, et al. Uncovering small RNA-mediated responses to phosphate deficiency inby deep sequencing [J]. Plant Physiology, 2009, 151(4): 2120-2132.

[6] Wang L, Zeng H Q, Song J, et al. miRNA778 andare involved in phosphate homeostasis in[J]. Plant Science, 2015, 238: 273-285.

[7] Jia X Y, Shen J, Liu H, et al. Small tandem target mimic-mediated blockage of microRNA858 induces anthocyanin accumulation in tomato [J]. Planta, 2015, 242(1): 283-293.

[8] Jiang X L, Huang K Y, Zheng G S, et al. CsMYB5a and CsMYB5e fromdifferentially regulate anthocyanin and proanthocyanidin biosynthesis [J]. Plant Science, 2018, 270: 209-220.

[9] Cui X, Wang Y X, Liu Z W, et al. Transcriptome-wide identification and expression profile analysis of thefamily genes in[J].Functional & Integrative Genomics, 2018, 18(5): 489-503.

[10] Punyasiri PAN, Abeysinghe ISB, Kumar V, et al. Flavonoid biosynthesis in the tea plant: properties of enzymes of the prominent epicatechin and catechin pathways [J]. Arch. Biochem. Biophys, 2004. 431(1): 22-30.

[11] Singh K, Rani A, Kuma S, et al. An early gene of flavonoid pathway, flavanone 3-hydroxylase, exhibits a positive relationship with the concentration of catechins in tea () [J]. Tree Physiol, 2008, 28(9): 1349-1356.

[12] Lv H P, Dai W D, Tan J F, et al. Identification of the anthocyanins from the purple leaf coloured tea cultivar Zijuan (var.) and characterization of their antioxidant activities [J]. Journal of Functional Foods, 2015, 17: 449-458.

[13] Shen J Z, Zou Z W, Zhang X Z, et al. Metabolic analyses reveal different mechanisms of leaf color change in two purple-leaf tea plant(L.) cultivars [J]. Horticulture Research, 2018, 5(1): 7. DOI: 10.1038/s41438-017-0010-1.

[14] 费旭元, 林智, 梁名志, 等. 响应面法优化“紫娟”茶中花青素提取工艺的研究[J]. 茶叶科学, 2012, 32(3): 197-202.

[15] 吕海鹏, 梁名志, 张悦, 等. 特异茶树品种“紫娟”不同茶产品主要化学成分及其抗氧化活性分析[J]. 食品科学, 2016, 37(12): 122-127.

[16] Wen M, Shen Y, Shi S H, et al. miREvo: An Integrative microRNA evolutionary analysis platform for next-generation sequencing experiments [J]. BMC Bioinformatics, 2012, 13(1): 140. DOI: 10.1186/1471-2105-13-140.

[17] Friedlander M R, Mackowiak S D, Li N, et al. miRDeep2 accurately identifies known and hundreds of novel microRNA genes in seven animal clades [J]. Nucleic Acids Research, 2012, 40(1): 37-52.

[18] Wu H J, Ma Y K, Chen T, et al. PsRobot: a web-based plant small RNA meta-analysis toolbox [J]. Nucleic Acids Research, 2012, 40(W1): W22-W28. DOI: 10.1093/nar/gks554.

[19] 蒋会兵, 夏丽飞, 田易萍, 等. 基于转录组测序的紫芽茶树花青素合成相关基因分析[J]. 植物遗传资源学报, 2018, 19(5): 967-978.

[20] Wang L, Feng Z, Wang X, et al. DEGseq: an R package for identifying deferentially expressed genes from RNA-seq data [J]. Bioinformatics, 2010, 26(1): 136-138.

[21] 谢小芳, 添先凤, 江昌俊, 等. 茶树低温胁迫下microRNA实时定量PCR内参基因的筛选[J]. 茶叶科学, 2015, 35(6): 596-604.

[22] Zhang Y, Zhu X J, Chen X, et al. Identification and characterization of cold-responsive microRNAs in tea plant () and their targets using high-throughput sequencing and degradome analysis [J]. BMC Plant Biology, 2014, 14: 271. DOI:10.1186/s12870-014-0271-x.

[23] Chen J L, Zheng Y, Qin L, et al. Identification of miRNAs and their targets through high-through put sequencing and degradome analysis in male and female[J]. BMC Plant Biology, 2016, 16(1): 80. DOI: 10.1186/s12870-016-0770-z.

[24] Mecchia M A, Debernardi J M, Rodriguez R E, et al. MicroRNA miR396 andsynergistically regulate leaf development [J]. Mechanisms of Development, 2013, 130(1): 2-13.

[25] Zhang W, Xie Y, Xu L, et al. Identification of microRNAs and their target genes explores miRNA-mediated regulatory network of cytoplasmic male sterility occurrence during anther development in radish (L.) [J]. Frontiers in Plant Science, 2016, 7: 1054. DOI: 10.3389/fpls.2016.01054.

[26] Yang X, Zhao Y, Xie D, et al. Identification and functional analysis of microRNAs involved in the anther development in cotton genic male sterile line Yu98-8A [J]. International Journal of Molecular Sciences, 2016, 17(10): 1677. DOI: 10.3390/ijms17101677.

[27] Liu J, Yuan Y, Wang Y L, et al. Regulation of fatty acid and flavonoid biosynthesis by miRNAs in[J]. Royal Society of Chemistry, 2017, 7: 35426-35437.

[28] Sun Y, Qiu Y, Duan M, et al. Identification of anthocyanin biosynthesis related microRNAs in a distinctive Chinese radish (L.) by high-throughput sequencing [J]. Molecular Genetics & Genomics, 2017, 292(1): 215-229.

[29] Shen E M, Singh S K, Ghosh J S, et al. The miRNAome of: identification, expression analysis, and potential roles of microRNAs in regulation of terpenoid indole alkaloid biosynthesis [J]. Scientific Reports, 2017, 7: 43027. DOI: 10.1038/srep43027.

[30] Liu N, Tu L, Wang L, et al. MicroRNA 157-targetedgenes regulate floral organ size and ovule production in cotton [J]. BMC Plant Biology, 2017, 17: 7. DOI: 10.1186/s12870-016-0969-z.

[31] Yang F X, Cai J, Yang Y, et al. Overexpression of microRNA828 reduces anthocyanin accumulation in[J]. Plant Cell Tiss Organ Cult, 2013, 115(2): 159-167.

[32] Baksa I, Nagy T, Barta E, et al. Identification ofmicroRNAs and their targets using high through put sequencing and degradome analysis [J]. BMC Genomics, 2015, 16(1): 1025. DOI: 10.1186/s12864-015-2209-6.

[33] Wang L, Du H Y, Ta-na W Y. Genome-wide identification of MicroRNAs and their targets in the leaves and fruits ofusing high-through put sequencing [J]. Frontiers in Plant Science, 2016, 7: 1632. DOI: 10.3389/fpls.2016.01632.

Screening of miRNA Related to Anthocyanin Synthesis in Tea Cultivar ‘Zijuan’ Based on High Throughput Sequencing

CHEN Linbo1,2, XIA Lifei1, LIU Yue1, SUN Yunnan1, JIANG Huibing1,TIAN Yiping1, CHEN Liang2*

1. Tea Research Institute, Yunnan Academy of Agricultural Sciences/Yunnan Engineering Research Centerof Tea Germplasm Innovation and Matching Cultivation /Yunnan Provincial Key Laboratory of Tea Science, Menghai 666201, China; 2. Tea Research Institute, Chinese Academy of Agricultural Sciences, Hangzhou 310008, China

In order to screen the miRNAs related to anthocyanin synthesis in 'Zijuan', the miRNA library was constructed by using the tea cultivars 'Zijuan' (ZJ), 'Yunkang 10' (YK) and 'Fuding Dabaicha' (FD). In this study, 46 known miRNAs and 67 unknown miRNAs were identified, and 765 annotated target genes were predicted. A total of 24 miRNAs were screened out, which were differentially expressed between ZJ and YK, and between ZJ and FD. Four target miRNAs involved in the regulation of anthocyanin synthesis, namely miR828a, miR845c, novel_14 and novel_87, were identified by target gene analysis of 24 differentially expressed miRNAs. The predicted target genes include transcription factor genes,,,,and(4-coincyl-CoA-linked enzyme),(dihydroflavonol 4-reductase) and(UDP-glucoside flavonoid 3--glucosyltransferase). Eight differentially expressed miRNAs were analyzed by RT-PCR and the results were consistent with transcriptome analysis. Finally, a theoretical basis for further study of regulation mechanisms related to anthocyanin biosynthesis in tea plant was provided.

, anthocyanin synthesis, high throughput sequencing, microRNAs, target genes

S571.1；Q52

A

1000-369X(2019)06-681-11

2019-03-04

2019-05-25

国家自然科学基金项目（No.31560220）、国家重点研发项目（2016YFD0200903）、国家茶叶额产业体系（CARS-19）、云南省人才培养计划项目（2015HB105）

陈林波，男，研究员，主要从事茶树资源育种研究。

liangchen@mail.tricaas.com