辣椒不同发育时期果实的转录组分析

王红娟，唐荣莉，蒋晓英，官 玲，雷开荣，黄启中，林 清

(1. 重庆市农业科学院生物技术研究所，重庆 401329；2. 重庆市农业科学院蔬菜花卉研究所, 重庆 401329)

【研究意义】辣椒(Capsicumsp.)是一类重要的经济作物，它既能够满足全世界人们对辛辣味的喜爱，同时还具有非常重要的营养价值和医学价值。辣椒素类物质的合成受遗传控制，并受栽培环境的影响。辣椒素类物质积累从果实开始发育的第10～20天开始，并在绿熟期或转色期达到最高[1-2]。不同栽培品种的辣椒素类物质含量变异极大，其史高维尔指数(scoville heat unit, SHU)低至0，高达2.5×106[3]。研究辣椒素生物合成途径，揭示其调控因子，可为辣椒生物技术育种工作提供分子基础。【前人研究进展】辣椒的辣味有无由位于2号染色体上的Pun1基因控制[4]。Qin等[5]推测Pun1基因(AT3)的两个拷贝可能协同作用，共同合成并积累含量不同的辣椒素，决定不同品种辣椒的辣度。通过转录水平分析和基因沉默技术，PAL、Ca4H、COMT、KAS、pAMT、AT3、ACL、DHAD、DT和FAT基因被证明参与辣椒素类物质的生物合成途径[6-12]。Mazourek等[13]建立了辣椒素类物质生物合成模型(CapCyc Model)，共包括52个候选酶类相关基因。近年，随着新一代测序技术的迅速发展和测序成本的降低，辣椒组学研究有了重大进展[5, 11, 14-17]。研究者通过转录组和全基因组分析获得了CapCyc模型中部分或者全部候选基因[5, 11, 15]，通过对转录组数据的比较分析发现，参与辣椒素合成途径的基因呈现组织/发育阶段特异性表达，但是仅有少数基因在辣椒果实合成辣椒素阶段显著表达[5]。这些结果进一步为辣椒素的生物合成提供了重要的候选基因。然而，目前的研究主要关注辣椒素合成途径中的结构基因，对调控因子的研究相对较少，目前报道的只有MYB和ERF 2个转录因子家族[18-24]。其中最近的研究表明，一个R2R3-MYB转录因子CaMYB31可能是辣椒素类物质合成途径的主要调控因子，该转录因子和合成途径中的Ca4H、COMT、KAS、pAMT和AT3等结构基因的表达呈正相关，且和果实发育过程中辣椒素和二氢辣椒素的含量也呈正相关[21-22, 25]。另一个MYB家族转录因子CaMYB108受到茉莉酸诱导，参与调控辣椒素类物质的合成[24]。此外，Keyhaninejad等[18]研究显示2个ERF转录因子ERF和JERF的表达与不同辣椒的辣椒素含量正相关，从而推测其对辣椒素类物质有调节功能。毛艳辉[19]通过表达分析和基因沉默判定一个ERF转录因子(CaERF1)参与辣椒素合成的调控。酵母双杂试验显示JERF与Pun1之间存在相互作用关系[23]。【本研究切入点】本研究以高辣和低辣度的2个一年生辣椒品种为研究对象，选取品种开始积累辣椒素和高含辣椒素的2个时期，利用新一代测序技术进行RNA-seq测序。【拟解决的关键问题】比较不同辣度的2个品种及其不同发育时期的转录组数据，筛选并鉴定差异表达基因，重点研究差异表达基因的功能分类，差异表达基因中的转录因子家族和CapCyc模型中的候选基因，研究辣椒素合成的调控机制，以期为辣椒的生物技术育种工作提供分子基础。

1 材料与方法

1.1 试验材料

本研究使用高辣度的辣椒材料754-3-1-1-1和低辣度的渝椒5号，种子来源于重庆市农业科学院蔬菜花卉研究所。754-3-1-1-1为朝天椒优良自交系，高效液相色谱测得辣椒素和二氢辣椒素含量之和为4.61 mg/g，为高辣材料[3]。渝椒5号是由83-1-1为母本、76-3-1-1为父本配制成的适宜春秋栽培的一代杂种，该品种果实长牛角形，微辣。实验材料露天种植，分别采集花后16 d左右(定义为KS快速生长期)和花后35 d左右LS绿熟期的果实。每个单株为1个生物学重复，754-3-1-1-1取3个生物学重复，渝椒5号取1个生物学重复。样品采集后放入液氮速冻，然后保存于-80 ℃超低温冰箱中，用于后续试验。

1.2 RNA-seq测序

754-3-1-1-1的3个生物学重复分别编号为L754KS1、L754KS2、L754KS3和L754LS1、L754LS2、L754LS3，渝椒5号2个时期分别编号为YJ5KS和YJ5LS。采用Invitrogen(美国)TRIzol Reagent试剂，参照说明书操作规程提取果实(包括胎座)总RNA。质检合格后送上海美吉生物医药科技有限公司构建RNA-seq文库，并在Illumina HiSeq 2000平台上测序。

1.3 转录组数据分析

转录组数据通过美吉生物云平台进行分析。为保证数据质量，在分析前对测序得到的原始数据进行过滤以减少数据噪音。过滤后保留下来的clean reads与辣椒(C.annuum)参考基因组比对，获得的Mapped Data用于下一步分析。为鉴定754-3-1-1-1不同生物学重复数据间的关系，对不同生物学重复数据进行主成分分析(PCA)。根据测序得到的基因表达水平的RPKM值[26]，通过Cufflinks软件进行基因差异表达分析[27]。差异显著基因筛选标准为|log2(Fold Change)|≥1且P<0.05。然后利用Goatools[28]对差异表达基因进行GO富集分析。

1.4 荧光定量RT-qPCR分析

为了验证RNA-seq数据的可靠性，随机选择7个差异表达基因进行实时荧光定量RT-qPCR表达量验证。利用在线引物设计工具Primer-BLAST设计引物(表1)，内参基因选用β-actin。RT-qPCR扩增每个样品重复3次。采用2-ΔΔCt方法分别计算差异表达基因在2个辣椒品种的2个果实发育时期的相对表达量[29]。

表1 RT-qPCR分析使用的引物序列

2 结果与分析

2.1 转录组测序质量和Reads分析

所有样品均测序得到超过4 Mb的原始数据。对原始数据进行过滤后，各样品均有超过99%的clean reads，低质量的reads占比不到1%，说明测序质量好。将clean reads与辣椒基因组进行比对，除L754KS2和L754LS3 2个样品外，其余样品均有超过84%的unique reads进入后续分析(表2)。

图1 754-3-1-1-1生物学重复主成分分析Fig.1 Principal component analysis for 754-3-1-1-1 samples

表2 样本Reads数分布情况

主成分分析可以降低数据的复杂性，深入挖掘样品之间的关系和变异大小。为鉴定754-3-1-1-1果实发育的2个时期各3个生物学重复转录组数据的关系，判断是否存在离群样品，本研究对其进行主成分分析。结果发现754-3-1-1-1LS期的3个重复聚在一起，说明这3个样品的数据相似性较高；KS期的2个重复(L754KS2和L754KS3)聚在一起，而样品L754KS1偏离另外2个重复，属离群样品(图1)，故后续分析不包括样品L754KS1。

2.2 荧光定量RT-qPCR验证基因表达

RT-qPCR结果显示，除了1个基因(ID: LOC107868336)在L754LS_vs_YJ5LS对比中表达情况不符外，RT-qPCR和RNA-seq测序得到的基因差异表达情况一致(图2)。表明RNA-seq数据可靠，可进行后续的差异表达基因分析。

图2 差异表达基因RT-qPCR验证Fig.2 Validation of DEG by RT-qPCR

2.3 基因差异性表达分析

高辣材料754-3-1-1-1和低辣品种渝椒5号的KS期和LS期分别设置4组对比进行基因差异表达分析。结果显示，L754LS_vs_L754KS有139个基因差异表达，其中17.27%的基因在LS期上调表达，82.73%的基因下调表达；YJ5LS_vs_YJ5KS发现1903个基因差异表达，其中61.59%的基因在LS期上调表达，38.41%下调表达。L754KS_vs_YJ5KS有5550个基因差异表达，其中51.17%在753-3-1-1-1中上调表达，48.83%的基因下调表达；L754LS_vs_YJ5LS有3220个基因差异表达，其中55.59%的基因在754-3-1-1-1中上调表达，44.41%的基因下调表达(表3)。除了L754LS_vs_L754KS，其余3组对比表达量上调的基因数目均多于下调基因数目。高辣材料754-3-1-1-1的KS和LS期差异表达基因最少，仅有139个。相反，2个品种在果实KS期和LS期都有较多的基因差异表达，尤其是L754KS_vs_YJ5KS，有5550之多。这些结果表明高辣材料754-3-1-1-1果实发育前后期基因表达模式变化较少，而754-3-1-1-1和渝椒5号两者间从果实发育早期就有较多的基因表达模式不同。

表3 不同对比组别的差异表达基因上调下调数目

借助Venn图对2个品种分别2个时期的差异表达基因进行展示(图3)，4组对比共有的差异表达基因有12个，高辣和低辣品种共有的差异表达基因有13个。L754LS_vs_L754KS有7个(5%)特有的差异表达基因，YJ5LS_vs_YJ5KS有429个(23%)特有的差异表达基因。2个品种在KS期有2667个(48%)特有的差异表达基因，在LS期有871个(27%)特有的差异表达基因。这些结果表明辣度不同的2个品种在辣椒素开始合成的时期有较多的基因表达模式不同，而随着果实的发育，基因表达差异变小。

图3 差异表达基因的Venn图Fig.3 Venn diagram of DEG

2.4 差异表达基因的GO富集分析

GO功能富集分析结果表明，754-3-1-1-1 2个生长时期的差异表达基因主要富集在转录调控区序列特异性DNA结合和儿茶酚氧化酶活性(图4-A)；渝椒5号2个生长时期的差异表达基因主要富集在多糖代谢过程、葡聚糖代谢过程等各类代谢过程(图4-B)；2个品种KS期差异表达基因主要富集在代谢过程、防御应答以及胁迫和刺激应答(图4-C)，LS期主要富集在光合作用、光捕获和蛋白质—生色团连锁(图4-D)。

***表示FDR<0.001，**表示FDR<0.01,*表示FDR<0.05*** indicates FDR<0.001, ** indicates FDR<0.01, * indicates FDR<0.05图4 差异表达基因的GO功能富集分析Fig.4 GO enrichment analysis of DEG identified from 4 pair comparisons

2.5 CapCyc模型基因和转录因子鉴定

本研究从4组对比数据的差异表达基因中共鉴定到18个CapCyc模型基因，其中PAL等5个基因分别有2～3个拷贝，因此共有24条基因序列差异表达(表4)。其中L754KS_vs_YJ5KS和 L754LS_vs_YJ5LS 2组对比中分别有13和17个CapCyc模型基因差异表达，YJ5LS_vs_YJ5KS对比发现3个，L754LS_vs_L754KS未发现。4个基因的表达差异倍数大于3倍，表达差异较大。

表4 差异表达基因中的CapCyc模型基因及其表达差异情况

续表4 Continued table 4

此外，本研究从L754LS_vs_L754KS、YJ5LS_vs_YJ5KS、L754KS_vs_YJ5KS和L754LS_vs_YJ5LS 4组对比数据的差异表达基因中鉴定到MYB家族转录因子36个，ERF家族转录因子4个(表5)。L754KS_vs_YJ5KS差异表达的转录因子最多，有25个(13个下调，12个上调)；L754LS_vs_YJ5LS次之，有19个(7个下调，12个上调)。其中11个为2组共有，表达量差异方向一致的有9个(表5)。36个MYB转录因子中有9个表达差异较大，在至少一组对比中表达差异超过5倍，其中LOC107852848在L754LS_vs_YJ5LS对比中、LOC107874571在L754KS_vs_YJ5KS对比中差异表达超过8倍，LOC107852202在YJ5LS_vs_YJ5KS对比中差异表达高达10.62倍。另有4个MYB转录因子(LOC107852599、LOC107879588、LOC107867180和LOC107846862)在除L754LS_vs_L754KS之外的3组对比中都差异显著。

表5 差异表达基因中的MYB和ERF转录因子及其表达差异情况

续表5 Continued table 5

3 讨 论

诸多研究已经系统地阐述了辣椒素类物质的生物合成，并对其结构基因的功能有了一定的认识。但是对于不同辣度品种间辣椒素含量差异巨大的原因知之甚少，尤其是对调控辣椒素生物合成的转录因子研究较少[30]。本研究基于RNA-seq测序技术，对不同辣度的2个辣椒品种754-3-1-1-1和渝椒5号的2个不同发育时期的果实进行转录组分析。通过L754LS_vs_L754KS、YJ5LS_vs_YJ5KS、L754KS_vs_YJ5KS和L754LS_vs_YJ5LS 4组数据对比，鉴定了一批差异表达基因。结果显示，除了L754LS_vs_L754KS，其他3组比较表达量上调的基因数目均多于下调基因数目。这表明辣度更高的样本相对于辣度更低的样本有更多的基因表达量升高。高辣材料754-3-1-1-1的KS和LS期差异表达基因仅有139个，且Venn图显示仅有7个(5%)是这2个时期特有，表明高辣材料754-3-1-1-1果实发育前后期基因表达模式变化较少。相反，2个品种在果实KS期和LS期都有较多的基因差异表达，尤其是L754KS_vs_YJ5KS，高达5550个。Venn图显示在KS期有2667个(48%)特有的差异表达基因，在LS期有871个(27%)特有的差异表达基因。这些结果显示754-3-1-1-1和渝椒5号两者间在辣椒素开始合成的时期有较多的基因表达模式不同，而随着果实的发育，基因表达差异变小。Zhang等[31]关于“鬼椒王”研究结果也显示同样的规律，10 vs 20、20 vs 30、30 vs 40 DPA(花后10、20、30、40 d对比)差异表达基因数量递减。我们推测辣椒素含量受果实发育早期基因表达差异的影响更大。

Mazourek等[13]建立了辣椒素类物质生物合成模型(CapCyc Model)，共包括52个结构基因。本研究中我们鉴定到18个CapCyc模型基因在不同辣度和果实不同发育时期中差异表达，其中PAL、CCoAOMT、CCR3个基因在至少一组对比中表达量差异超过3倍。以往的研究表明辣椒素生物合成途径中结构基因的表达水平与辣椒素的含量相关，并受到协同调节[31-32]。并且，QTL定位和GWAS分析表明结构基因能够控制辣椒素的含量[33]。本研究中，PAL、CCoAOMT和CCR等结构基因差异表达较大，表明这些基因很可能具有调控辣椒素合成的作用。

对辣椒素生物合成途径调控因子的研究较少，目前只有2个R2R3-MYB家族转录因子CaMYB31、CaMYB108和2个ERF家族转录因子被证明是辣椒素类物质合成途径的主要调控因子[18, 21, 24]。本研究从4组对比的差异表达基因中鉴定到36个MYB家族转录因子和4个ERF家族转录因子，其中包括前人研究中的CaMYB31(LOC107877282)和ERF(LOC107864060)基因。在本研究中，CaMYB31在样品L754KS中表达水平较高，在L754LS_vs_L754KS和L754KS_vs_YJ5KS 2组对比间表达差异显著，且差异表达倍数均超过4倍。之前的研究表明CaMYB31与辣椒素合成途径中的Ca4H、COMT、KAS、pAMT、AT3等结构基因的表达呈正相关[21]。本研究发现pAMT和AT3-D1基因在L754KS_vs_YJ5KS对比中表达差异显著，且差异表达的方向和CaMYB31一致，表明这2个基因可能受CaMYB31调控。此外，我们还发现9个转录因子(7个MYB家族和2个ERF家族)在L754KS_vs_YJ5KS和L754LS_vs_YJ5LS 2组对比中同时差异表达，且差异方向一致，说明这些转录因子在果实积累辣椒素的不同时期都在发挥作用。

4 结 论

本研究对2种不同辣度的辣椒材料754-3-1-1-1和渝椒5号果实不同发育时期的转录组进行了测序和分析。通过对比不同发育时期果实的基因表达情况，我们推测果实发育早期基因的表达水平对辣椒素的含量有更大的影响。PAL、CCoAOMT、CCR等辣椒素合成途径中的结构基因在不同样品中差异表达较大，表明这些基因可能具有调控辣椒素合成的作用。此外，本研究还鉴定出一批可能参与调控辣椒素合成的MYB和ERF家族转录因子，为后续研究提供参考，也为辣椒的生物技术育种工作提供分子基础。