‘仓方早生’桃及其早熟芽变不同发育时期果实的转录组分析

彭佳伟，张叶，寇单单，杨丽，刘晓飞，张学英，陈海江，田义

‘仓方早生’桃及其早熟芽变不同发育时期果实的转录组分析

1河北农业大学园艺学院，河北保定 071000；
2河北农业大学山区研究所/河北省山区农业技术创新中心/国家北方山区农业工程技术研究中心，河北保定 071001；
3河北省保定市农业农村局，河北保定 071000；
4河北省盐山县望树镇中学，河北沧州 061300

【目的】通过对桃品种‘仓方早生’及其早熟芽变不同发育时期的果实进行转录组分析，挖掘参与调控桃果实成熟的关键因子，为深入研究桃果实成熟调控机理提供理论依据。【方法】以桃品种‘仓方早生’及其早熟芽变为试材，每个品种分别选择长势一致的样品树5株，分别于花后30 d（对应‘仓方早生’c1、早熟芽变y1）、45 d（对应c2、y2）、59 d（对应c3、y3）、71 d（对应c4、y4）及89 d（对应c5）对不同发育时期的桃去皮果肉进行取样和转录组测序，并利用实时荧光定量PCR（qRT-PCR）对筛选的差异表达基因进行定量验证；
利用GO和KEGG对‘仓方早生’及其早熟芽变的差异表达基因进行分析；
基于差异表达基因构建加权基因共表达网络分析（weighted gene co- expression network analysis，WGCNA），从中鉴定出与果实成熟密切相关的枢纽模块和枢纽基因。【结果】将处于果实相同发育时期的转录组数据进行比较，得到y1与c1、y2与c2、y3与c4和y4与c5四组对比数据，共筛选出差异表达基因4 395个，其中上调表达基因2 212个，下调表达基因2 183个。其中包括10个乙烯、11个脱落酸和18个生长素合成及其信号转导途径基因，并构建了10个IAA蛋白与预测互作ARF蛋白间的相互作用网络。由GO分类统计结果可知，差异表达基因在生物过程板块主要集中于细胞过程、代谢过程和单体过程；
在细胞组分板块主要聚集于膜和细胞组分；
在分子功能板块主要富集于结合蛋白和催化活性等方面。‘仓方早生’及其早熟芽变果实的差异基因主要集中在y3与c4和y4与c5对比组中，这些差异基因大多被富集到分子功能中结合活力、氧化还原酶活性等方面。对差异表达基因进行KEGG通路分析表明，在果实生长发育成熟过程中伴随着多种次生代谢产物的变化，如倍半萜和三萜生物合成、类黄酮生物合成、类胡萝卜素的生物合成和-亚麻酸代谢等。同时，本研究发现生长素信号转导途径在不同时间节点均有富集，这意味着植物激素信号转导通路对果实成熟具有极为重要的作用。【结论】在‘仓方早生’及其早熟芽变不同发育时期果实的差异表达基因中，大量激素信号转导途径基因特别是生长素信号途径基因发生了富集，这些基因可能在调控果实发育中具有重要作用，可对这些候选的和功能及其如何通过相互作用调控果实成熟的机制进行进一步的解析。

桃；
果实成熟；
转录组分析；
加权基因共表达网络分析

【研究意义】桃是重要的落叶果树之一，因其果实口感甜美和气味芳香而深受消费者喜爱。但桃果实在贮藏和运输过程中极易损伤变质，生产中主要依靠种植不同成熟期的桃品种来延长市场供应期。桃的果实成熟期是一个十分重要的农艺经济性状，选育不同成熟期的桃品种，对丰富市场供给、提高种植者的经济收益具有重要意义。本研究通过对‘仓方早生’及其早熟芽变进行转录组分析，从转录组数据中筛选出与果实成熟有关的代谢通路，深入挖掘与桃果实成熟有关的候选基因，分析相关基因表达变化与差异，为挖掘调控桃果实成熟的关键基因及其分子机制奠定基础。【前人研究进展】果实成熟是一个非常复杂的过程，包括颜色、质地、香气和生物活性化合物等多个方面的变化[1-2]。大量研究表明，果实成熟涉及到多种植物激素的调控[3-4]。其中，乙烯（Ethylene，ETH）是呼吸跃变型果实成熟及其后续代谢变化的主要触发和协调因子，直接影响果实品质[5-7]。如HAYAMA等[8]所述，乙烯影响桃果实软化的起始和进程，其浓度是调节软化速度的重要因素，桃果实的成熟软化在转录水平上受乙烯的调控。乙烯的合成由ACC合酶（1-aminocyclopropane-1-carboxylate synthase，ACS）和ACC氧化酶（1-aminocyclopropane-1-carboxylic acid oxidase，ACO）两种酶催化[9-10]，是乙烯合成的关键基因，在桃果实中的转录随着成熟的进程而增加[11]。有研究表明，乙烯转录因子（ERFs）可以通过在不同程度上调控乙烯合成关键基因来影响果实成熟进程[12-13]。同时，脱落酸（abscisic acid，ABA）作为一种关键的植物激素，也参与了呼吸跃变型果实的成熟和衰老，被认为是包括桃在内的几种呼吸跃变型果实成熟的决定因素之一[14-15]。ABA可以启动桃果实成熟，而且在果实成熟后期与乙烯一起积累[16]。在桃果实成熟过程中，ABA浓度显著增加，加速桃果实的软化和成熟[17]。此外，ABA信号转导中的基因对果实成熟也有着相似的调控作用[18]。Zeng等[16]从桃中鉴定的10个2C型蛋白磷酸酶（protein phosphatase 2C，PP2C）基因中，和可以引发果实成熟，并在果实成熟的后期发挥重要作用，在番茄中也有类似报道[19]。WANG等[17]在桃基因组中检测到两个9-顺式-环氧类胡萝卜素双加氧酶（9-cis-epoxycarotenoid dioxygenase，NCED）基因的表达水平与ABA一致，表明两个基因可能是果实成熟过程中ABA积累的重要参与者。生长素（Auxin）也是果实发育成熟过程中的重要激素，有重要的调控作用[20-22]。通过TATSUKI等[23]的研究可以看出，在不同类型的桃果实成熟过程中，硬质桃果皮组织中的IAA浓度较低，且没有增长；
但溶质桃果皮组织中的IAA水平显著增加，导致果实软化。因此，桃成熟进程中IAA浓度可以调控果实的成熟软化。GUAN等[24]从桃基因组中鉴定到23个，表达模式分析得出和在完全成熟时显著上调表达，且上调倍数最大，为果实发育、果实成熟进程所必需；
此外，由基因表达变化可知，等均可能在果实发育和成熟过程中发挥功效。部分研究表明3种植物激素的调控因子协同调控桃果实成熟，如PAN等[25]发现，在处于低表达水平时，对硬质型（stony hard，SH）桃果实成熟期间生长素的积累产生影响，从而抑制乙烯生物合成基因的表达，导致果实不能软化；
WANG等[26]发现是调节桃成熟过程中生长素信号转导的重要因子，与一起在桃果实成熟过程中整合生长素和乙烯信号，并通过结合和激活和启动子促进ABA生物合成，来调节桃果实的成熟和软化。另有研究表明，通过抑制两个ABA生物合成基因（、）和一个细胞壁降解基因（）的表达来调节果实成熟[27]。【本研究切入点】目前，桃果实成熟的机制尚未完全阐明。本研究将结合转录组数据全面准确地对桃果实成熟的通路进行挖掘筛选与分析。【拟解决的关键问题】本研究以桃品种‘仓方早生’及其早熟芽变为研究试材，进行转录组测序，通过分析‘仓方早生’及其早熟芽变的转录组数据，挖掘出与参与调控果实成熟有关的激素信号途径中的关键基因。

1.1 材料处理

试材‘仓方早生’和早熟芽变于2020年取自保定市顺平西闫庄试验园，两品种盛花期相近，均在4月6—7日，花后59 d早熟芽变的果实开始着色，花后71 d达到成熟；
‘仓方早生’花后89 d达到成熟，两品种果实均属于溶质型。

每个品种分别选择长势一致的样本树5株。根据果实生长动态的4个发育时期（图1）共采样5次：花后30 d（第1次膨大期，对应c1与y1）、花后45 d（硬核期，对应c2与y2）、花后59 d（早熟芽变第2次膨大期，对应c3与y3）、花后71 d（早熟芽变成熟期，‘仓方早生’第2次膨大期，对应c4与y4）及花后89 d（‘仓方早生’成熟期，对应c5）。每次分别随机采取‘仓方早生’和早熟芽变果实10个，取5 g去皮果肉，切碎，迅速用液氮处理后置于-80℃保存。

1.2 方法

1.2.1 总RNA提取及检测、cDNA文库构建和转录组测序 采用多糖多酚植物总RNA提取试剂盒（TIANGEN，China）分别对‘仓方早生’及其早熟芽变在不同采样时间的果实进行总RNA的提取，所有操作按照说明书进行。取2 μL RNA样品进行琼脂糖凝胶电泳，检测RNA样品的质量。用NanoDrop One（Thermo Fisher Scientific，USA）对RNA质量进行初步检测。利用HiFi Script gDNA Removal RT Master Mix（Beijing ComWin Biotech Co.,Ltd.）试剂盒合成第一链cDNA。采用Illumina HiSeqTM 2000进行转录组表达谱测序，由北京华诺时代科技有限公司完成。

c：仓方早生Kurakato；
y：早熟芽变Early-ripening mutant。下同 The same as below

1.2.2 转录组测序质量评估 测序reads经fastp（version 0.18.0）[28]过滤后得到clean reads，通常以Q20、Q30大小衡量测序准确度，通常≥85%[29]。将得到的clean reads与桃参考基因组（http://www.rosaceae.org/species/prunus_persica/genome_v1.0）比对，进行基因结构注释、基因表达分析和基因功能预测等。

1.2.3 转录组差异基因的实时荧光定量PCR验证 使用Primer Premier Software设计qRT-PCR引物，利用天根反转录试剂盒（TIANGEN，China）合成cDNA。以为内参基因，按照Ultra SYBR One Step RT-qPCR Kit （Beijing ComWin Biotech Co.,Ltd.）提供的说明书，在Light Cycler_96进行qRT-PCR分析，反应条件：94℃预变性60 s；
95℃ 10 s，58℃ 10 s，72℃ 10 s，45个循环；
95℃ 10 s，65℃ 60 s，97℃ 1 s。采用2-ΔΔCT法[30]计算基因的相对表达量。

1.2.4 样品间差异表达基因及富集分析 通过EdgeR包（http://www.bioconductor.org/packages/2.12/bioc/ html/edgeR.html）计算筛选差异表达基因（differentially expressed genes，DEG），筛选条件是|log2Fold Change|＞1，-value＜0.05，log2（FC）为样品1与样品2的FPKM差异倍数的对数值[31]。采用GO和KEGG两个数据库[32]，对筛选出的差异基因进行功能注释以及通路显著性富集分析。

1.2.5 植物激素相关代谢途径中差异表达基因的筛选与热图分析 根据样品间差异表达基因的GO功能分析与KEGG代谢途径富集分析结果，从转录组数据中找出与桃果实成熟有关的代谢途径。本试验重点分析植物激素相关途径并从中筛选出差异表达基因，进一步明确基因的功能注释，并利用Tbtools软件制作热图对这些基因进行分析。

1.2.6 加权基因共表达网络的构建与分析 首先对前期得到的数据进行预处理和质量控制。利用GSG（Good Samples Genes）函数来检测数据集中是否有离群值和缺失值。使用R中的WGCNA（v1.47）包构建共表达网络[31]。筛选样本间表达量差异最大的前10%的基因作为后续构建加权共表达网络的数据，直至筛选到枢纽基因。将枢纽基因的蛋白序列上传至KEGG官网（http://www.KEGG.jp/）得到基因的功能注释。对富集到生长素信号转导通路上的AUX/IAA进行蛋白互作预测，得到若干个ARF蛋白，将IAA与ARF的相互作用关系通过Cytoscape进行可视化分析。

2.1 数据质量评估

对‘仓方早生’及其早熟芽变果实不同生长发育时期的27个样品进行RNA提取以及转录组测序分析。经过滤去除低质量的reads后与桃参考基因组进行比对，比对率均在90%以上。

2.2 转录组数据的实时荧光定量PCR验证

从转录组中选取变化趋势明显的8个DEG进行验证，包括5个果实成熟相关的DEG（、、、、）和3个随机挑选的DEG（、、），结果如图2所示。转录组测得的FPKM值与这8个基因相对表达量的变化趋势基本一致，且与基因已知正常表达量一致，表明RNA-seq数据可靠，可以进行后续的数据分析和关键基因挖掘。

2.3 桃‘仓方早生’与早熟芽变果实差异表达基因

把同一发育时期的‘仓方早生’和早熟芽变果实转录组数据进行两两对比，共获得差异表达基因4 395个，其中上调基因2 212个，下调基因2 183个。如图3-A所示，y1与c1、y2与c2、y3与c4和y4与c5中分别有282、152、1 409和2 552个差异表达基因。结果发现，果实发育后期两个阶段（y3与c4、y4与c5）的差异基因数量较发育前期明显增多。

通过维恩图探究果实不同成熟期各样品间转录组的不同，如图3-B所示，鉴定到‘仓方早生’及其早熟芽变在不同成熟时期涉及的差异表达基因存在较大差异，各个发育时期均存在差异表达的基因有13个。对上述13个基因进行KEGG通路分析，发现其中的光合作用途径（ko00195: Photosynthesis）与果实发育密不可分，铁氧还蛋白-NADP+还原酶（ferredoxin-NADP+reductase，FNR）作为光合电子传递链的末端氧化酶直接参与光合作用；
细胞色素P450途径（ko00199: Cytochrome P450）中注释到的，具有强大的催化作用，能通过植物新陈代谢来调节植物生长发育。

图2 差异基因的qRT-PCR验证

图3 y1与c1、y2与c2、y3与c4、y4与c5四组差异基因数（A）及样品间差异表达基因维恩图（B）

2.4 差异基因的GO富集分析

GO富集分析可以反映桃果实发育过程中DEG的主要生物学功能。本研究将处于同一发育时期的数据进行对比发现，果实中差异基因在生物过程（biological process）富集数量最多（图4）。其中，y1与c1组间比较发现282个DEG，被富集到24个GO通路；
y2与c2组间比较发现152个DEG，被富集到25个GO通路；
y3与c4组间比较发现1 409个DEG，被富集到32个GO通路；
y4与c5组间比较发现2 552个DEG，被富集到35个GO通路（图5）。

图4 ‘仓方早生’及其早熟芽变果实不同发育时期GO富集差异基因数目

2.5 差异表达基因的KEGG富集分析

将4个对比组中-value≤0.05的代谢途径列为显著富集的代谢途径，通过KEGG分析显示，y1与c1、y2与c2、y3与c4和y4与c5中分别有53、36、342和522个差异基因又被分别注释到多条代谢途径上，分别是66、50、189、179条。分析发现，在y1与c1组间主要富集在倍半萜和三萜生物合成、类黄酮生物合成、类胡萝卜素的生物合成途径；
在y2与c2组间中GABAergic突触途径富集程度最高，然后是FoxO信号通路和自噬调节途径；
在y3与c4组间，主要富集的代谢途径是光合作用、DNA复制、碳代谢；
y4与c5组间显著富集在减数分裂-酵母、细胞周期等方面（图6）。但是，在处于果实发育后期的y3与c4及y4与c5对比组中，植物信号转导途径中包含差异基因数均较多，分别为27与40个DEG。

2.6 植物激素相关代谢途径中差异表达基因的筛选与热图分析

本研究重点关注与植物激素有关的通路，从转录组中筛选到10个乙烯、11个脱落酸与18个生长素生物合成及其信号转导基因，并对此进行分析。

图5 ‘仓方早生’及其早熟芽变果实不同发育时期差异表达基因的GO功能类别

图6 ‘仓方早生’及其早熟芽变果实不同发育时期差异表达基因的KEGG富集

筛选出的10个与乙烯有关的差异表达基因中，包括2个乙烯生物合成基因和8个乙烯信号转导基因。如图7-A所示，的表达在不同时间点均具有显著差异，其中表达量高的（ppa001786m）在果实二次膨大期差异显著，早熟芽变较‘仓方早生’表达量上调。

在ABA代谢与信号转导通路中（图7-B），合成ABA的关键基因（ppa002804m）和分解ABA的关键基因P450单加氧酶（ppa005059m、MSTRG.14953.1）在‘仓方早生’及早熟芽变果实生长发育中表达差异显著。在果实发育及成熟过程中，表达丰度良好，两品种均于花后45 d达到最高值，成熟时，‘仓方早生’中表达量高于早熟芽变；
2个表达模式相近，呈下降趋势，且在果实成熟过程中，早熟芽变中的表达量较‘仓方早生’显著下调，通过调控ABA合成与分解途径中相关基因表达水平最终影响果实不同发育阶段的ABA积累。在ABA信号转导途径中，‘仓方早生’及早熟芽变从59 d之后到果实成熟，调控ABA受体基因及参与下游信号转导的基因表达差异显著，在PYR/PYL-PP2C-SnRK2 ABA信号转导系统中，与表达都可能对果实成熟发挥重要作用。PP2C基因（ppa008381m、ppa005286m、ppa006696m）等表达量较‘仓方早生’下调；
表达低且变化幅度小。

将参与生长素信号转导的18个差异表达基因进行表达量分析（图7-C）。在花后45 d之前，相关基因（ppa002065m）在‘仓方早生’及其早熟芽变中表达无显著差异；
硬核期之后，早熟芽变果实中ppa002065m表达量逐渐下降，‘仓方早生’中该基因表达量先缓慢增加后降低。在花后45 d，早熟芽变果实中属于家族的相关基因（ppa013846m、ppa013844m、ppa013765m、ppa013847m）表达量较‘仓方早生’上调，随着果实逐渐成熟，其表达量逐渐增加。在果实生长发育过程中，AUX/IAA家族差异显著基因共6个，其中除ppa010871m、ppa010501m在前期不表达外，其余4个基因表达量呈上升趋势。在花后59、71 d，早熟芽变中AUX/IAA相关基因ppa010342m、ppa010871m、ppa009416m表达量较‘仓方早生’上调，其中ppa010871m表达量上调幅度最大。随着果实发育成熟，ppa007483m、ppa010342m、ppa011935m、ppa011570m表达量逐渐增高。编码的ppa003134m表达丰度较高且在‘仓方早生’及其早熟芽变果实发育后期表达量差异显著。生长素合成限速酶基因（ppa007054m）的表达量在早熟芽变果实中自硬核期便迅速增加，而在‘仓方早生’中的表达趋势则是先缓慢增加，后期才迅速增加。

A：乙烯；
B：脱落酸；
C：生长素 A: Ethylene; B: Abscisic acid; C: Auxin

2.7 加权基因共表达网络分析

为了鉴定与果实发育相关的基因，对转录组数据中测得的34 375个基因进行加权基因共表达网络分析。首先筛选出样本间表达量差异最大的前10%的基因（3 483个）作为网络构建中的样本-基因数据集。然后计算网络权重，当power=10时，相关性达到0.85，此时平均连接度趋近于0，为最优软阈值数。根据power值获得临近矩阵和拓扑矩阵，得到的拓扑矩阵使用相异度对基因进行聚类，然后使用动态剪切法将树剪切成不同的模块，如图8所示。本研究鉴定到14个模块，各模块所包含的基因数目见图9。另外，各个模块均具有相对独立性，有利于后续模块分析的公平性和客观性。

经分析鉴定，模块turquoise中基因的差异表达可能是导致两品种果实成熟期出现差异的主要原因，以MM值和GS值筛选出隶属于该模块中的枢纽基因643个。将上述643个基因进行KEGG通路分析，共有356个基因被注释到了各种通路。经统计，生长素信号转导途径注释到的基因数目最多，包含了10个AUX/IAA基因，其表达量变化如图10所示。AUX/IAA基因在‘仓方早生’和早熟芽变果实中表达量总体呈上升趋势，但表达丰度呈现明显不同。当‘仓方早生’及早熟芽变果实达到成熟时，表达量最高，最小，且8自花后45 d表达量迅速增加。在‘仓方早生’及其早熟芽变的果实中，在各采样时间点表达丰度良好，而在果实发育早期不表达，在成熟时达到峰值。在表达丰度中等的AUX/IAA基因中，的表达量随果实发育成熟持续上升，当果实进入成熟时，上升趋势最明显。将上述IAA蛋白序列输入到String中，搜索与之互作的ARF蛋白，并绘制IAA与ARF之间的蛋白互作网络图，结果通过Cytoscape进行可视化（图11）。以Degree绘制点的大小和颜色衡量该节点在网络图中关键程度，以combined-score绘制线的粗细和颜色衡量互作关系的强度，点越大、越绿则节点越重要（ARF5＞ARF1＞AUX/IAA6＞AUX/IAA27＞ARF4），线越粗、越绿则互作更强（ARF5-AUX/IAA8＞ARF5-AUX/IAA20D＞ARF9-AUX/IAA20D＞ARF5-ARF4）。

图8 基于diss TOM的动态混合算法构建系统聚类树

图9 各个模块中基因数量

从转录组中查找预测到的8个ARF蛋白，除ARF9以外，其余均被鉴定到。因此，对余下7个ARF基因进行表达量分析，发现在‘仓方早生’及早熟芽变果实中，随着果实的发育进程，不同ARF基因间的表达趋势大为不同。在果实成熟终阶段，表达水平达到峰值。多数ARF基因在‘仓方早生’及其早熟芽变中表达趋势相似，而在两个品种中表达趋势相异（图12）。在果实发育成熟过程中的表达丰度保持良好，且蛋白互作预测表明（图11），与筛选得到的枢纽AUX/IAA基因均有互作关系。

图10 ‘仓方早生’和早熟芽变中IAA的表达量

图11 IAA和ARF的蛋白互作网络图

图12 ‘仓方早生’和早熟芽变中ARF基因表达量

3.1 ‘仓方早生’及其早熟芽变间的差异基因表达分析

本研究对‘仓方早生’及其早熟芽变不同发育时期果实进行转录组测序分析，发现早熟芽变较‘仓方早生’果实多数基因下调表达。花后30 d，‘仓方早生’及其早熟芽变之间的DEG数量较低，而花后71 d时‘仓方早生’及其早熟芽变之间的差异表达基因数量最高，且各组间除y4 vs c5上调基因数目多于下调基因外，其他3个时期的下调基因数目均多于上调基因，因此推测发育后期可能是影响早熟芽变果实成熟差异的关键。GO富集显示，4个对比组中注释到生物过程类的基因数量最多，结合KEGG通路分析发现，在果实生长发育成熟过程中，伴随着多种次生代谢途径基因的差异表达，猜测果实在成熟过程中次生代谢产物大量合成，这为更深入研究果实内在品质形成的分子代谢机理奠定了基础。

CHEN等[33]分析了两种肉质桃果实在3个不同发育阶段的RNA-seq数据，鉴定了包含52个在内的120个差异基因可能与果实成熟软化相关。ONIK等[34]对3个不同阶段的苹果进行RNA-Seq分析，结果显示，众多差异表达基因参与了乙烯、脱落酸、生长素、赤霉素等激素生物合成途径，MYB、NAC、WRKY和HSF等转录因子在成熟前和成熟后的果实间也存在差异表达；
除乙烯外，ABA等激素在调节苹果果实成熟中也发挥着关键作用，并可能与乙烯信号相互作用。殷亚蕊[35]通过对‘京红’桃与其晚熟芽变不同发育阶段果实进行RNA-Seq分析，从中筛选出23个生长素信号转导途径的差异基因，和表达量随着果实发育呈下调趋势，且果实发育后期在晚熟突变体中的表达量显著低于野生型；
此外，随着果实发育成熟进程，多个IAA基因在野生型果实发育后期的表达呈明显上调。本研究也发现植物激素信号转导途径在不同时间节点均有富集。‘仓方早生’及其早熟芽变果实发育进程差异主要表现在硬核期及以后，早熟芽变果实比‘仓方早生’的硬核期短，推测可能是硬核期基因表达差异的影响引起后续植物中生长素、乙烯、脱落酸相关基因表达差异，激素相关基因表达差异又进一步加速果实成熟相关基因的表达。

3.2 乙烯、脱落酸和生长素相关基因在果实发育过程中的差异表达分析

前期研究表明，在桃果实乙烯生物合成中发挥关键作用，是控制桃果实成熟软化的主要因子[36-39]。本研究发现，随着果实成熟，表达量增加，当果实成熟时，早熟芽变中表达量显著高于‘仓方早生’。番茄、柿和甜樱桃的表达量也表现出在成熟转变期达到高峰，后熟期迅速下降的趋势[40-41]。本研究发现在果实成熟过程中，相较于‘仓方早生’，早熟芽变中2个表达水平呈下降趋势，且ABA合成关键酶显著下调，多数PP2C基因表达量下调，暗示脱落酸积累的多少对果实发育成熟具有十分重要的作用。此外，鉴定到了多个生长素信号响应因子、和，这些基因的启动子区均包含保守生长素响应元件，生长素响应因子（auxin response factor，ARF）通过与这些元件特异性结合来调控生长素的响应[42-44]，其中对于生长素极为敏感[45-47]。在果实发育前期，早熟芽变与‘仓方早生’相比，表达量上调，推测当生长素达到一定浓度时，首先响应生长素信号，通过对生长素的合成及运输调控果实发育。另有研究表明IAA基因可以直接或间接促进乙烯的生物合成，最终调控果实成熟进程[42,48]。本研究中，与的表达量随果实发育逐渐升高，并在果实成熟时表达量达到最高（图10）。根据基因表达量的变化与早熟芽变比‘仓方早生’的上调倍数，推测、、等参与果实的成熟。此外，本研究还鉴定到在果实达到成熟时表达丰度良好，但在果实发育后期，在两个品种中出现显著的差异表达，因此推测可能参与果实成熟差异调控过程。

RNA-seq获得的数据量大且条目繁多，WGCNA能够有效从大量数据中提取到关键信息来解释生物学现象[49]，因此，利用WGCNA解析具体的农艺性状，分析模块的生物功能，挖掘模块中的目标基因，将为解析复杂的农艺性状提供重要的参考依据。本研究运用WGCNA算法从转录组数据中挖掘出了关于桃果实成熟的相关基因，并应用该算法对‘仓方早生’及其早熟芽变果实不同发育时期的转录组数据通过构建无尺度的基因共表达网络，根据权重参数将数据进行聚类形成不同的模块，再以相关系数值的大小对基因模块的关联性划分，得到相关显著的模块turquoise，对模块中具有特异性的基因深度分析。然后对上述基因进行KEGG功能富集分析，进一步筛选到与果实成熟相关的10个AUX/IAA基因，为阐明调控桃果实成熟的分子网络奠定了基础。

本研究从桃‘仓方早生’及其早熟芽变果实不同发育时期的RNA-seq结果中筛选出与桃果实成熟相关的差异表达基因4 395个，其中上调基因2 212个，下调基因2 183个。通过WGCNA将差异基因划分为14个共表达模块，其中模块turquoise是与桃果实成熟度显著相关的特异性模块，结合GO和KEGG分析，共筛选到39个植物激素相关通路的差异表达基因，涉及10个乙烯、11个脱落酸和18个生长素生物合成及其信号转导基因。利用枢纽基因挖掘到10个作为候选基因，在此基础上预测到了7个ARF基因，并用其构建了基因互作调控网络。

[1] GIOVANNONI J, NGUYEN C, AMPOFO B, ZHONG S L, FEI Z J. The epigenome and transcriptional dynamics of fruit ripening. Annual Review of Plant Biology, 2017, 68: 61-84.

[2] SHEN J Y, TIEMAN D, JONES J B, TAYLOR M G, SCHMELZ E, HUFFAKER A, BIES D, CHEN K S, KLEE H J. A 13-lipoxygenase, TomloxC, is essential for synthesis of C5 flavour volatiles in tomato. Journal of Experimental Botany, 2014, 65(2): 419-428.

[3] KLIE S, OSORIO S, TOHGE T, DRINCOVICH M F, FAIT A, GIOVANNONI J J, FERNIE A R, NIKOLOSKI Z. Conserved changes in the dynamics of metabolic processes during fruit development and ripening across species. Plant Physiology, 2014, 164(1): 55-68.

[4] KUMAR R, KHURANA A, SHARMA A K. Role of plant hormones and their interplay in development and ripening of fleshy fruits. Journal of Experimental Botany, 2014, 65(16): 4561-4575.

[5] LERSLERWONG L, THIPPAN S, RUGKONG A, IMSABAI W. Expression pattern of ethylene-related genes in response to preharvest chemical treatments during development and ripening of mangosteen fruit (L.). Walailak Journal of Science and Technology, 2021, 18(3): 6663-6672.

[6] YANG Y Y, SHAN W, KUANG J F, CHEN J Y, LU W J. Four HD-ZIPs are involved in banana fruit ripening by activating the transcription of ethylene biosynthetic and cell wall-modifying genes. Plant Cell Reports, 2020, 39(3): 351-362.

[7] LI S, ZHU B Z, PIRRELLO J, XU C J, ZHANG B, BOUZAYEN M, CHEN K S, GRIERSON D. Roles of RIN and ethylene in tomato fruit ripening and ripening-associated traits. The New Phytologist, 2020, 226(2): 460-475.

[8] HAYAMA H, SHIMADA T, FUJII H, ITO A, KASHIMURA Y. Ethylene-regulation of fruit softening and softening-related genes in peach. Journal of Experimental Botany, 2006, 57(15): 4071-4077.

[9] 阚娟, 刘俊, 金昌海. 桃果实成熟软化与细胞壁降解相关糖苷酶及乙烯生物合成的关系. 中国农业科学, 2012, 45(14): 2931-2938. doi: 10.3864/j.issn.0578-1752.2012.14.016.

KAN J, LIU J, JIN C H. Study on the relationship between peach fruit softening, cell wall degradation related glycosidase and ethlylene biosynthesis. Scientia Agricultura Sinica, 2012, 45(14): 2931-2938. doi: 10.3864/j.issn.0578-1752.2012.14.016. (in Chinese)

[10] 刘进平. 乙烯生物合成关键酶基因研究进展. 热带农业科学, 2013, 33(1): 51-57.

LIU J P. Advances in research on key enzyme genes of ethylene biosynthesis. Chinese Journal of Tropical Agriculture, 2013, 33(1): 51-57. (in Chinese)

[11] TRAINOTTI L, TADIELLO A, CASADORO G. The involvement of auxin in the ripening of climacteric fruits comes of age: The hormone plays a role of its own and has an intense interplay with ethylene in ripening peaches. Journal of Experimental Botany, 2007, 58(12): 3299-3308.

[12] HAN Y C, KUANG J F, CHEN J Y, LIU X C, XIAO Y Y, FU C C, WANG J N, WU K Q, LU W J. Banana transcription factor MaERF11 recruits histone deacetylase MaHDA1 and represses the expression of MaACO1 and expansins during fruit ripening. Plant Physiology, 2016, 171(2): 1070-1084.

[13] LI T, XU Y X, ZHANG L C, JI Y L, TAN D M, YUAN H, WANG A D. The jasmonate-activated transcription factor mdmyc2regulates ethylene response factor and ethylene biosynthetic genes to promote ethylene biosynthesis during apple fruit ripening. The Plant Cell, 2017, 29(6): 1316-1334.

[14] NICOLAS P, LECOURIEUX D, KAPPEL C, CLUZET S, CRAMER G, DELROT S, LECOURIEUX F. The basic leucine zipper transcription factor abscisic acid response element-binding factor2 is an important transcriptional regulator of abscisic acid-dependent grape berry ripening processes. Plant Physiology, 2014, 164(1): 365-383.

[15] LIU Y F, ZHAO Y N, CHAI L P, ZHOU J Q, YANG S, KOU X H, XUE Z H. Transcriptome profiling of abscisic acid-related pathways in-silenced tomato fruits. Transactions of Tianjin University, 2021, 27(6): 473-486.

[16] ZENG W F, TAN B, DENG L, WANG Y, ASLAM M, WANG X B, QIAO C K, WANG Z Q, FENG J C. Identification and expression analysis of abscisic acid signal transduction genes during peach fruit ripening. Scientia Horticulturae, 2020, 270: 109402.

[17] WANG X B, ZENG W F, DING Y F, WANG Y, NIU L, YAO J L, PAN L, LU Z H, CUI G C, LI G H, WANG Z Q.PpERFpositively regulates ABA biosynthesis by activatingtranscription during fruit ripening in peach. Horticulture Research, 2019, 6: 19.

[18] KLINGLER J P, BATELLI G, ZHU J K. ABA receptors: The start of a new paradigm in phytohormone signalling. Journal of Experimental Botany, 2010, 61(12): 3199-3210.

[19] SUN L, WANG Y P, CHEN P, REN J, JI K, LI Q, LI P, DAI S J, LENG P. Transcriptional regulation of SlPYL, SlPP2C, and SlSnRK2 gene families encoding ABA signal core components during tomato fruit development and drought stress. Journal of Experimental Botany, 2011, 62(15): 5659-5669.

[20] LEYSER O. Auxin signaling. Plant Physiology, 2018, 176(1): 465-479.

[21] WANG Y C, WANG N, XU H F, JIANG S H, FANG H C, SU M Y, ZHANG Z Y, ZHANG T L, CHEN X S. Auxin regulates anthocyanin biosynthesis through the Aux/IAA-ARF signaling pathway in apple. Horticulture Research, 2018, 5: 59.

[22] 余佳, 李阳, 龚硕, 关伟, 刘悦萍. 桃果实生长素结合蛋白ABP1的组织定位及蛋白表达分析. 中国农业科学, 2015, 48(5): 921-930. doi: 10.3864/j.issn.0578-1752.2015.05.10.

YU J, LI Y, GONG S, GUAN W, LIU Y P. Tissue location and protein expression analysis of auxin binding protein ABP1 in peach fruit (L.). Scientia Agricultura Sinica, 2015, 48(5): 921-930. doi: 10.3864/j.issn.0578-1752.2015.05.10. (in Chinese)

[23] TATSUKI M, NAKAJIMA N, FUJII H, SHIMADA T, NAKANO M, HAYASHI K I, HAYAMA H, YOSHIOKA H, NAKAMURA Y. Increased levels of IAA are required for system 2 ethylene synthesis causing fruit softening in peach (L. Batsch). Journal of Experimental Botany, 2013, 64(4): 1049-1059.

[24] GUAN D, HU X, DIAO D H, WANG F, LIU Y P. Genome-wide analysis and identification of the Aux/IAA gene family in peach. International Journal of Molecular Sciences, 2019, 20(19):4703-4721.

[25] PAN L Z W F, NIU L, LU Z H, LIU H, CUI G C, ZHU Y Q, CHU J F, LI WP, FANG W C, CAI Z G, LI G H, WANG Z Q., a strong candidate gene for the stony hard phenotype in peach (L. Batsch), participates in IAA biosynthesis during fruit ripening. Journal of Experimental Botany, 2015, 66(22): 7031-7044.

[26] WANG X B, PAN L, WANG Y, MENG J R, DENG L, NIU L, LIU H, DING Y F, YAO J L, NIEUWENHUIZEN N J, AMPOMAH- DWAMENA C, LU Z H, CUI G C, WANG Z Q, ZENG W F.andPpERFform a positive feedback loop to regulate peach fruit ripening by integrating auxin and ethylene signals. Plant Science: an International Journal of Experimental Plant Biology, 2021, 313: 111084.

[27] WANG X B, ZENG W F, DING Y F, WANG Y, NIU L, YAO J L, PAN L, LU Z H, CUI G C, LI G H, WANG Z Q. Peach ethylene response factor PpeERFrepresses the expression of ABA biosynthesis and cell wall degradation genes during fruit ripening. Plant Science: an International Journal of Experimental Plant Biology, 2019, 283: 116-126.

[28] CHEN S F, ZHOU Y Q, CHEN Y R, GU J. Fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics, 2018, 34(17): i884-i890.

[29] 赵湾湾. 基于RNA-seq的红肉番木瓜果实成熟相关基因的挖掘[D]. 福州: 福建农林大学, 2019.

ZHAO W W. Mining related genes of red meat papaya fruit ripening based on RNA-seq [D]. Fuzhou: Fujian Agriculture and Forestry University, 2019. (in Chinese)

[30] 徐献斌, 耿晓月, 李慧, 孙丽娟, 郑焕, 陶建敏. 基于转录组分析ABA促进葡萄花青苷积累相关基因. 中国农业科学, 2022, 55(1): 134-151. doi: 10.3864/j.issn.0578-1752.2022.01.012.

XU X B, GENG X Y, LI H, SUN L J, ZHENG H, TAO J M. Transcriptome analysisof genes involved inABA-induced anthocyanin accumulation in grape. Scientia Agricultura Sinica, 2022, 55(1): 134-151. doi: 10.3864/j.issn.0578-1752.2022.01.012. (in Chinese)

[31] 郭永春, 王鹏杰, 金珊, 侯炳豪, 王淑燕, 赵峰, 叶乃兴. 基于WGCNA鉴定茶树响应草甘膦相关的基因共表达模块. 中国农业科学, 2022, 55(1): 152-166. doi: 10.3864/j.issn.0578-1752.2022. 01.013.

GUO Y C, WANG P J, JIN S, HOU B H, WANG S Y, ZHAO F, YE N X. Identification of Co-expression gene related to tea plant response to glyphosate based on WGCNA. Scientia Agricultura Sinica, 2022, 55(1): 152-166. doi: 10.3864/j.issn.0578-1752.2022.01. 013. (in Chinese)

[32] 何平, 李林光, 王海波, 常源升, 李慧峰. 遮光性套袋对桃果实转录组的影响. 中国农业科学, 2017, 50(6): 1088-1097. doi: 10.3864/j.issn.0578-1752.2017.06.010.

HE P, LI L G, WANG H B, CHANG Y S, LI H F. Effects of shading fruit with opaque paper bag on transcriptome in peach. Scientia Agricultura Sinica, 2017, 50(6): 1088-1097. doi: 10.3864/j.issn.0578- 1752.2017.06.010. (in Chinese)

[33] CHEN C W, GUO J, CAO K, ZHU G R, FANG W C, WANG X W, LI Y, WU J L, XU Q, WANG L R. Identification of candidate genes associated with slow-melting flesh trait in peach using bulked segregant analysis and RNA-seq. Scientia Horticulturae, 2021, 286: 110208.

[34] ONIK J C, HU X J, LIN Q, WANG Z D. Comparative transcriptomic profiling to understand pre- and post-ripening hormonal regulations and anthocyanin biosynthesis in early ripening apple fruit. Molecules, 2018, 23(8): 1908-1927.

[35] 殷亚蕊. ‘京红’桃及其晚熟芽变果实发育特征和转录组测序分析[D]. 秦皇岛: 河北科技师范学院, 2020.

YIN Y R. Development characteristics and transcriptome sequencing analysis of ‘Jinghong’ peach and its late-maturing bud mutation fruit [D]. Qinhuangdao: Hebei Normal University of Science & Technology, 2020. (in Chinese)

[36] 王小贝. 桃和协同调控桃成熟软化的分子机制研究[D]. 武汉: 华中农业大学, 2019.

WANG X B. Study on the molecular mechanism of synergistic regulation of peach ripening and softening byand[D]. Wuhan: Huazhong Agricultural University, 2019. (in Chinese)

[37] YUAN H, YUE P T, BU H D, HAN D G, WANG A D. Genome-wide analysis ofandgenes in pear ().Cellular & Developmental Biology-Plant, 2020, 56(2): 193-199.

[38] LI T, TAN D M, YANG X Y, WANG A D. Exploring the apple genome reveals six ACC synthase genes expressed during fruit ripening. Scientia Horticulturae, 2013, 157: 119-123.

[39] 王宏, 杨王莉, 蔺经, 杨青松, 李晓刚, 盛宝龙, 常有宏. 早熟砂梨‘苏翠1号’与其亲本‘翠冠’‘华酥’成熟果实差异代谢产物及差异基因比较分析. 园艺学报, 2022, 49(3): 493-508.

WANG H, YANG W L, LIN J, YANG Q S, LI X G, SHENG B L, CHANG Y H. Comparative metabolic and transcriptomic analysis of ripening fruit in pear cultivars of ‘Sucui l’ ‘Cuiguan’ and ‘Huasu’. Acta Horticulturae Sinica, 2022, 49(3): 493-508. (in Chinese)

[40] JI K, KAI W B, ZHAO B, SUN Y F, YUAN B, DAI S J, LI Q, CHEN P, WANG Y, PEI Y L, WANG H Q, GUO Y D, LENG P.and: Key genes involved in ABA metabolism during tomato fruit ripening. Journal of Experimental Botany, 2014, 65(18): 5243-5255.

[41] ZHAO S L, QI J X, DUAN C R, SUN L, SUN Y F, WANG Y P, JI K, CHEN P, DAI S J, LENG P. Expression analysis of theandgenes that regulate homeostasis of abscisic acid during the maturation of persimmon fruit. The Journal of Horticultural Science and Biotechnology, 2015, 87(2): 165-171.

[42] ITO J, FUKAKI H, ONODA M, LI L, LI C Y, TASAKA M, FURUTANI M. Auxin-dependent compositional change in Mediator in ARF7- and ARF19-mediated transcription. Proceedings of the National Academy of Sciences of the United States of America, 2016, 113(23): 6562-6567.

[43] 欧春青, 姜淑苓, 王斐, 赵亚楠. 梨全基因组生长素反应因子(ARF)基因家族鉴定及表达分析. 中国农业科学, 2018, 51(2): 327-340. doi: 10.3864/j.issn.0578-1752.2018.02.012.

OU C Q, JIANG S L, WANG F, ZHAO Y N. Genome-wide identification and expression analysis of auxin response factor () gene family in pear. Scientia Agricultura Sinica, 2018, 51(2): 327-340. doi: 10.3864/j.issn.0578-1752.2018.02.012. (in Chinese)

[44] TAO S B, ESTELLE M. Mutational studies of the Aux/IAA proteins inreveal novel insights into their function. The New Phytologist, 2018, 218(4): 1534-1542.

[45] HOU K, WU W, GAN S S., a small auxin up RNA gene, is involved in the promotion of leaf senescence in. Plant Physiology, 2013, 161(2): 1002-1009.

[46] 朱宇斌, 孔莹莹, 王君晖. 植物生长素响应基因SAUR的研究进展. 生命科学, 2014, 26(4): 407-413.

ZHU Y B, KONG Y Y, WANG J H. Research advances in auxin-responsive SA UR genes. Chinese Bulletin of Life Sciences, 2014, 26(4): 407-413. (in Chinese)

[47] KONG Y Y, ZHU Y B, GAO C, SHE W J, LIN W Q, CHEN Y, HAN N, BIAN H W, ZHU M Y, WANG J H. Tissue-specific expression of SMALL AUXIN UP RNA41 differentially regulates cell expansion and root meristem patterning in. Plant and Cell Physiology, 2013, 54(4): 609-621.

[48] LIU M C, CHEN Y, CHEN Y, SHIN J H, MILA I, AUDRAN C, ZOUINE M, PIRRELLO J, BOUZAYEN M. The tomato ethylene response factor Sl-ERF.B3 integrates ethylene and auxin signaling via direct regulation of Sl-Au/IAA27. The New Phytologist, 2018, 219(2): 631-640.

[49] 秦天元, 孙超, 毕真真, 梁文君, 李鹏程, 张俊莲, 白江平. 基于WGCNA的马铃薯根系抗旱相关共表达模块鉴定和核心基因发掘. 作物学报, 2020, 46(7): 1033-1051.

QIN T Y, SUN C, BI Z Z, LIANG W J, LI P C, ZHANG J L, BAI J P. Identification of drought-related co-expression modules and hub genes in potato roots based on WGCNA. Acta Agronomica Sinica, 2020, 46(7): 1033-1051. (in Chinese)

Transcriptome Analysis of Peach Fruits at Different Developmental Stages in Peach Kurakato Wase and Early-Ripening Mutant

1Horticultural Department, Agricultural University of Hebei, Baoding 071000, Hebei;2Mountainous Areas Research Institute, Hebei Agricultural University/Technology Innovation Center for Agriculture in Mountainous Areas of Hebei Province/National Engineering Research Center for Agriculture in Northern Mountainous Areas, Baoding 071001, Hebei;3Baoding Municipal Bureau of Agriculture and Rural Affairs, Baoding 071000, Hebei;4Wangshu Town Middle School of Yanshan County, Hebei Province, Cangzhou 061300, Hebei

【Objective】 In this study, transcriptome analyses were carried out on the fruits of Kurakato Wase peach and its early-ripening mutant at different developmental stages. The key factors involved in fruit ripening regulation were explored, so as to provide a theoretical basis for further study on the regulation mechanism of fruit ripening. 【Method】 The flesh of Kurakato Wase peach and its early-ripening mutant was sampled at 30 d, 45 d, 59 d, 71 d, and 89 d after anthesis, and transcriptome analyses were performed on the above samples. The candidate differentially expressed genes (DEGs) were verified by quantitative real-time PCR (qRT-PCR). The biological function of DEGs were analyzed through GO function and KEGG pathway. The weighted gene co-expression network analysis (WGCNA) was constructed to identify the hub modules and hub genes closely related to fruit ripening. 【Result】 Four comparison groups including y1 vs c1, y2 vs c2, y3 vs c4 and y4 vs c5 were obtained based on fruit development stages. A tatal of 4 395 DEGs were identified with 2 212 up- and 2 183 down-regulated genes.There were 10, 11 and 18 candidate genes involved in ethylene, abscisic acid and auxin synthesis and signal transduction, respectively. The interaction networks between 10 IAA proteins and their predictive interacting proteins ARF were constructed. GO function revealed that the DEGs were mainly enriched in cellular processes, metabolic processes and monomeric processes in the biological process category; in cell component category, DEGs were mainly enriched in membranes and cellular components; in molecular function category, DEGs were mainly enriched in binding protein and catalytic activity. There were more DEGs in comparison groups y3 vs c4 and y4 vs c5, and these DEGs mainly enriched in molecular functions, such as binding and catalytic activity. The KEGG pathway analysis showed that a variety of secondary metabolites changed during fruit development and ripening, such as sesquiterpene and triterpenoid biosynthesis, flavonoid biosynthesis, carotenoid biosynthesis, and-linolenic acid metabolism. In addition, auxin signal transduction pathway was found to be enriched at different time nodes. 【Conclusion】 Among DEGs, a large number of hormone signal transduction pathway genes, especially auxin signal pathway genes, were enriched, and these genes might play an extremely important role during fruit ripening. The functions of candidate genesandand the molecular regulation of fruit ripening would be further elucidated in the future studies.

peach; fruit ripening; RNA-seq; WGCNA

10.3864/j.issn.0578-1752.2023.05.012

2022-04-27；

2022-07-21

财政部和农业农村部：国家现代农业产业技术体系（CARS-30-2-03）、热杂果现代种业科技创新团队（21326310D）

彭佳伟，E-mail：1358345303@qq.com。张叶，E-mail：20528036@qq.com。彭佳伟和张叶为同等贡献作者。通信作者张学英，E-mail：yyzxy@hebau.edu.cn。通信作者陈海江，E-mail：chenhaijiang2001@163.com。通信作者田义，E-mail：tianyi@hebau.edu.cn

（责任编辑 赵伶俐）