不同发育时期蓝果忍冬果实的转录组分析

孙燕 王金刚 臧丹丹 赵恒田 刘淑华

（1.东北农业大学园艺园林学院，哈尔滨 150030；2.中国科学院东北地理与农业生态研究所黑土区农业生态重点实验室，哈尔滨 150081）

蓝果忍冬（Lonicera caeruleaL.）为忍冬科忍冬属植物，原产于北半球温带地区，现广泛分布于我国各地，是我国重要的野生果树资源［1-3］。蓝果忍冬果实出汁率高，有香气，营养价值丰富，具有调节血压、减缓衰老等作用［4-5］，因其富含花青素、氨基酸、维生素、类黄酮等多种活性物质，被誉为“饮料之王”［6］。‘中科蓝1 号’是我国第一个具有自主知识产权和国家备案的蓝果忍冬新品种，植株健壮，丰产性强，果实为长椭圆形，口感酸甜无苦味，是优良的鲜果品种，适宜东北寒区广泛种植［7］。

花青素是植物体内一类重要的次生代谢物质，属于类黄酮类化合物，在赋予植物色彩时起主要作用［8-9］。在果实生长发育过程中，果皮中色素会发生一系列变化，随着果实的成熟，叶绿素在水解酶的作用下逐渐分解消失，类胡萝卜素略有增加；花青素在细胞质中合成，以苯丙氨酸为直接前体，由一系列结构基因编码的合成酶催化形成，并经过一系列的运输后贮存于液泡中；浆果的大小、色泽、硬度均会随着果实的成熟发生显著变化［10-11］。近年来，关于蓝果忍冬果实花青素的研究受到广泛关注，Zhou 等［12］对蓝果忍冬品种‘蓓蕾’果实中的花青素提取并纯化，发现提取的花青素ABL-2 的主要成分是花青素-3-葡萄糖苷、花青素-3，5-二葡萄糖苷、花青素-3-芦丁糖苷和芍药-3-葡萄糖苷，且能够显著抑制细胞生长，诱导DNA 损伤，杀死肿瘤细胞，改善小鼠的生存状态；Paereckait 等［13］利用分光光度法对11 种蓝果忍冬果实中的酚类化合物和花青素含量等进行提取研究并比较不同蓝果忍冬品种的抗菌活性。然而，目前有关蓝果忍冬基因信息尚不明确，对蓝果忍冬果实转录组的相关研究少有报道，严重限制了蓝果忍冬的开发和利用。

本研究以蓝果忍冬‘中科蓝1 号’果实为试验材料，使用高通量测序平台对绿熟期果实和成熟期蓝果忍冬果实进行转录组测序分析，通过组装建立蓝果忍冬转录组数据库，探究蓝果忍冬果实不同发育时期的功能基因和调控网络，旨在为后续了解蓝果忍冬果实发育提供理论参考。

1 材料与方法

1.1 材料

试验材料为蓝果忍冬‘中科蓝1 号’品种，由中国科学院东北地理与农业生态研究所培育，分别采取绿熟期（花后28 d，果实成长到体积不再发生变化，果皮颜色为绿色）和成熟期（花后38 d，果皮颜色全部变为蓝紫色，果实未软化）两个时期的果实（果皮和果肉），各取3 个生物学重复，用液氮速冻后保存在-80℃冰箱中用于后续转录组测序和花色苷含量的测定（图1）。

图1 不同发育时期的‘中科蓝1 号’果实Fig.1 Fruits of ‘Zhongkelan No.1’ at different developmental stages

1.2 方法

1.2.1 不同时期的花色苷总量测定 参照高效液相色谱法（high performance liquid chromatography，HPLC）［14］，利用安捷伦液相色谱系统（Agilent Technologies 1260 Infinity Ⅱ Prime，US）检测花青素种类及含量，流动相A 为含3% 甲酸的双蒸水，流动相B 为100% 乙腈。洗脱程序为0 min 90%A，10% B；16 min 86% A，14% B；20 min 90% A，10% B；流速为1 mL/min，柱温为30℃，检测波长为520 nm，每个样品设3 个重复；分别以矢车菊素-3，5-双葡萄糖苷、飞燕草素-3-O-葡萄糖苷、飞燕草素-3-O-芸香糖苷、矢车菊素-3-O-葡萄糖苷、矢车菊素-3-O-芸香糖苷、天竺葵素-3-O-葡萄糖苷、芍药素-3-O-葡萄糖苷、芍药素-3-O-芸香糖苷标准品绘制标准曲线，含量由峰面积计算方法得出，单位为mg/100 g FW。

1.2.2 RNA 提取 采用CTAB 法提取蓝果忍冬不同发育时期的总RNA，经PrimeScript One Step RT-PCR Kit Ver.2（TaKaRa）反转录为cDNA 并将cDNA 稀释10 倍用于后续RT-qPCR 实验。样品总RNA 经检测合格后送至北京百迈客生物技术有限公司运用高通量平台进行RNA-Seq 测序。转录组样品编码为G：“绿熟期”，B：“成熟期”。

1.2.3 转录组测序与分析 为了获得高质量的转录组数据，采用illumina novaseq 6000 平台进行双末端测序，截除掉原始序列接头和引物序列，对低质量的序列进行过滤处理，并利用Trinity 软件对生成的高质量测序数据进行组装，得到转录本序列Unigene；利用BLAST 软件在NR、GO、Swiss-Prot、COG、KEGG 和eggNOG4.5 数据库中进行比较分析并预测Unigene 的氨基酸序列，使用HMMER软件将Unigene 序列在Pfam 数据库比对从而获得Unigene 的功能信息，并使用 TBtools 软件绘制热图。

1.2.4 差异表达基因的筛选和富集分析 用DESeq2软件对蓝果忍冬绿熟期和成熟期果实的转录组数据进行组间差异表达分析，筛选标准为FDR（False Discovery Rate）＜0.01 且FC（Fold Change）≥2，并利用Goatools 软件对筛选出的差异表达基因进行COG 注释分析、GO 功能富集分析和KEGG 富集分析。

1.2.5 差异表达基因的RT-qPCR 验证 根据蓝果忍冬转录组数据随机筛选出10 个差异基因，在NCBI数据库中比对并设计定量引物（表1）。选用LcACT1（MT344113）、LcTUB1（MT344114）作为内参基因，利用Agilent AriaMx 实时荧光定量PCR 仪（Agilent Stratagene，USA）对其进行实时荧光定量PCR 反应，程序为94℃预变性30 s，94℃变性15 s，57℃退火30 s，72℃延伸35 s，共44 个循环结束，65-95℃熔解曲线读取，用2-ΔΔCT算法计算处理。

表1 转录组数据RT-qPCR 验证引物序列Table 1 Primer sequences used in RT-qPCR validation of transcriptome data

2 结果

2.1 蓝果忍冬‘中科’两种时期的花色苷含量变化

利用高效液相色谱（high performance liquid chromatography，HPLC）检测技术和紫外检测器对520 nm 波长下进行扫描，检测到的峰图如图2所示，可以看出，蓝果忍冬‘中科’绿熟期果实未检测出花色苷成分（图2-a），成熟期果实中检测到的花青素种类共有8 种（A1-A8）（图2-b），总花色苷含量为598.55 mg/100 g FW。其中，矢车菊素-3-O-葡萄糖苷含量最高，占总花色苷的86.28%，是蓝果忍冬‘中科’果实中主要的花色苷（表2），说明不同发育阶段蓝果忍冬‘中科’果实中花色苷类化合物组成和含量有较大差异。

表2 ‘中科蓝1 号’成熟期果实中花色苷含量Table 2 Anthocyanins contents in the ripening stage fruits of ‘Zhongkelan No.1’

图2 ‘中科蓝1 号’不同发育时期果实高效液相色谱Fig.2 HPLC of ‘Zhongkelan No.1’ fruit at different developmental stages

2.2 转录组测序数据统计

选择蓝果忍冬绿熟期果实和成熟期果实各3 次生物学重复，共6 个样品，转录组测序共得到40.32 Gb Clean Data，每个样品Clean Data 均大于6.07 Gb，GC 含量为46.32%-49.38%，Q30 均值在93.89%及以上，表明测序数据较好，可用于后续分析（表3）。

表3 测序数据统计结果Table 3 Statistics of sequencing data

6 个蓝果忍冬果实的转录组数据经过质量分析处理后共得到236 033 条转录本序列，平均长度为1 371.99 bp，N50 长2 186 bp；组装后获得81 777个单基因序列，且长度大于1 kb 的Unigene 共有17 339 个（表4），对所获得的单基因序列与NR、Swiss-Prot、KEGG、COG、KOG、GO 和Pfam 数据库进行功能注释比对，注释成功的Unigene 共有28 992 条，占全部Unigene 的35.45%（表5）。

表4 组装结果统计Table 4 Assembly result statistics

表5 单基因功能注释统计Table 5 Single-gene function annotation statistics

2.3 不同发育时期差异表达基因筛选结果

根据评估标准FDR（False Discovery Rate）＜0.01 且FC（Fold Change）≥2，共筛选出3 247 个差异表达基因，其中有1 642 个基因在上调表达，1 605 个基因在下调表达，从图3可以清楚地看出，在两个生长时期蓝果忍冬果实中基因表达水平存在较大差异。

图3 差异表达基因的火山图Fig.3 Volcano map of differentially expressed genes

2.4 不同发育时期差异表达基因COG功能分类

利用COG 数据库对蓝果忍冬果实两个发育时期的差异表达基因进行对比分析（图4）发现，注释到COG 通路中的差异表达基因有1 081 个，分布在23 条途径中，占全部差异基因的37.88%，其中，碳水化合物运输和代谢（carbohydrate transport and metabolism）是最大通路，所占数量最多，其次是一般功能预测（general function prediction only），差异基因在脂质运输和代谢（lipid transport and metabolism）通路和次级代谢物生物合成、运输和分解代谢（secondary metabolites biosynthesis，transport and catabolism）通路占比也较高，说明在果实成熟发育过程中代谢物质变化较大。

图4 不同发育时期差异表达基因COG 注释分类统计图Fig.4 COG annotation classification of differentially expressed genes at different developmental stages

2.5 不同发育时期差异表达基因GO功能富集结果

对蓝果忍冬果实两个发育时期的差异表达基因进行对比分析，发现注释到GO 通路中的差异表达基因有1 908 个，分别映射到51 个代谢通路中（图5），在生物过程（biological process）中，主要富集在代谢过程（metabolic process）和细胞过程（cellular process），其次是单个有机体过程（single-organism process）；在细胞组分（cellula component）通路中，膜部分（membrane）占比最多、其次是细胞（cell）和细胞部分（cell part）；在分子功能（molecular function）通路中，催化活性（catalytic activity）和结合功能（binding）较其他通路显著富集。花青素是植物次生代谢产物，糖苷衍生物。推测可能归类于代谢过程（metabolic process）和催化活性（catalytic activity）的Unigenes 中，可能有参与蓝果忍冬花青素生物合成积累差异相关的良好候选基因。

图5 不同发育时期差异表达基因GO 功能富集结果Fig.5 Gene ontology annotation of differentially expressed genes at different developmental stages

2.6 不同发育时期差异表达基因KEGG功能富集结果

以q-value＜0.05 作为筛选差异表达基因在KEGG数据库通路的标准，由图6可知，2 854 个差异基因中共有1 034 个DEGs 注释到KEGG 通路中，其中植物激素信号转导（plant hormone signal transduction）途径是最大通路，有55 个，占比10.07%；其次是氨基酸生物合成（biosynthesis of amino acid）和碳代谢（carbon metabolism）途径，均为46 个，占比为8.42%；此外，淀粉和蔗糖代谢（starch and sucrose metabolism）通路与果实的发育也显著富集，有37 个，占比6.78%。

图6 不同发育时期差异表达基因KEGG 分类图Fig.6 KEGG classification of differentially expressed gene at different developmental stages

花青素是类黄酮化合物，且与苯丙氨酸代谢途径有关。对KEGG 代谢通路进行富集显著性筛选，得出富集最显著的20 条代谢通路（图7），其中，苯丙氨酸，酪氨酸和色氨酸的生物合成（phenylalanine，tyrosine and tryptophan biosynthesis）、类黄酮生物合成（flavonoid biosynthesis）、半胱氨酸和蛋氨酸代谢（cysteine and methionine metabolism）这3 条代谢途径被显著富集，说明可能存在大量与蓝果忍冬花青素生物合成积累相关的差异基因。

图7 不同发育时期差异表达基因KEGG 功能富集结果Fig.7 KEGG analysis of differentially expressed genes at different developmental stages

2.7 蓝果忍冬转录组测序数据的RT-qPCR验证

为了验证RNA-seq 的结果，对测序结果随机挑选了10 个差异表达基因进行RT-qPCR 验证，分析RT-qPCR 与RNA-seq 结果之间的相关性，每个基因进行3 次生物学重复。结果显示（图8），随机挑选的10 个差异表达基因的表达量变化与转录组基因的变化完全一致，说明转录组的测序数据和RT-qPCR的数据之间具有良好的一致性，表明本研究RNAseq 数据具有较高的可靠性。

图8 转录组测序数据的RT-qPCR 验证Fig.8 Validation of transcriptome sequencing data by RT-qPCR

3 讨论

近年来，随着分子生物学技术的发展，新一代高通量测序技术已被广泛应用于生物分析和医学研究等领域，因其具有成本低，测序速度快，能够估计不同发育时期和不同组织中的整体基因的表达量等特点，是当前无参考基因组的非模式生物的经济有效的选择［15-18］。目前，已有研究表明在蓝果忍冬果实或种子中分离和鉴定出酚类物质、抗氧化成分和挥发性物质等多种生物活性成分，包括黄酮类化合物、花青素、脂肪酸、甘油三酯和总甾醇等化合物，并证实花青素-3-葡萄糖苷为其主要的抗氧化成分［19-22］。本研究基于RNA-seq 测序技术，对‘中科蓝1 号’蓝果忍冬两个发育时期的果实进行转录组分析，序列经过组装之后，共获得81 777 条平均长度为763.35 nt 的Unigene，与前人报道的‘蓓蕾’品种果实发育阶段的基因数目相近［23］，比马尾松、枇杷、荔枝等木本植物组装的Unigene 多［24-26］，结果表明，组装的质量以及片段长度可以满足转录组数据分析的需要。有28 992 条Unigene 得到注释，占比35.45%，虽然比大叶女贞（50.20%）和蓝莓（73.60%）的注释率低，但与文冠果（36.15%）和红掌（38.87%）相近［27-30］，造成注释率低的原因可能是缺乏近缘物种的参考基因组信息，未得到注释的基因也可能是蓝果忍冬果实特异性新基因，有待进一步分析研究。

对比两个测序样本检测到的差异表达基因发现，表达量上调的基因数目多于下调基因数目，这表明成熟期果实相对于绿熟期果实的样本有更多的基因表达量。对这些基因进行KEGG 分析后发现有大量的色素形成相关途径被显著富集，包括与花青素代谢直接相关的苯丙氨酸、酪氨酸和色氨酸的生物合成途径、类黄酮生物合成途径和氨基酸生物合成途径，表明蓝果忍冬果实色泽调控的复杂性。GO 富集分析发现在果实绿熟期到转色期果实成熟过程中，参与细胞过程和代谢的基因较多，同时涉及到的功能最多的基因为结合和催化功能，从分子水平说明了这一过程中代谢旺盛。蓝果忍冬不同发育时期的差异表达基因在植物激素信号转导、碳代谢、氨基酸生物合成、脂肪代谢等多条代谢途径中富集，与类黄酮合成紧密相关，具体受哪些基因控制还需要做进一步信息学分析。

园艺产品品质形成的调控包括多个方面，在植物激素相关研究中，表明GA 对果实成熟软化［31］、有机酸及糖分积累［32］、果实着色［33］、果实大小和形状［34］等方面都有一定影响。比如，赵荣华等［35］研究表明GA 水平的调控可同时影响葡萄果实的形态和品质，另有研究发现生长素和赤霉素的相互作用在辣椒的生长发育过程中能平衡细胞分裂和细胞扩张［36］。在转录因子相关研究中，Medina-Puche等［37］研究表明YB9/FaMYB11、FabHLH3 和FaTTG1 在调控草莓花青苷合成过程中起关键作用，调控果皮花青苷的生物合成［38］。本研究发现蓝果忍冬‘中科蓝1 号’成熟期果实含有8 种花色苷成分，其中含量最高的花色苷为矢车菊素-3-O-葡萄糖苷，占总花色苷含量的86.28%，而绿熟期果实未检测出花色苷成分，推测在果实成熟过程中存在调控花色苷含量的转录因子，还需进一步研究。

4 结论

本研究对‘中科蓝1 号’蓝果忍冬绿熟期和成熟期果实进行了转录组测序，共有40.32 Gb Clean Data，每个样品Clean Data 均大于6.07 Gb，Q30 碱基百分比均在93.89%及以上，GC 含量为46.32%-49.38%；共筛选出3 247 个差异表达基因，包含1 642 个上调基因和1 605 个下调基因，为研究蓝果忍冬果实转色过程的的分子调控机制提供了参考。