不同硬度西瓜果皮的转录组测序及相关基因表达分析

张敬敬，李 冰，史宇凡，高秀瑞，潘秀清，宋 雪，武彦荣

(河北省农林科学院 经济作物研究所，河北 石家庄 050051)

西瓜在我国农业生产中占有重要地位，我国是西瓜的主要产区之一，每年西瓜产量约7 000万t[1]，占全球西瓜产量的67%[2]。河北省西瓜呈现区域化种植，优势产区明显，自“十二五”以来，衡水阜城和保定清苑西瓜播种面积在6 666.67 hm2以上[3]，西瓜生产区域化种植、上市时间集中，对品种的耐裂性、贮藏性、货架期品质提出了更高的要求。

西瓜果皮硬度在果实的耐裂性、贮藏性、货架期品质等方面起着重要作用[4]，生产上许多西瓜品种存在优质不耐裂、耐裂不优质现象，是造成西瓜生产过程中和采后产量损失的主要原因。通过对果皮硬度研究，在保证产量和品质的基础上改良果皮硬度，创制优质、果皮高硬度、耐裂的种质资源，在培育耐裂西瓜新品种上意义重大。

西瓜裂果是造成田间和产后减产的主要原因，果实开裂是一个复杂的性状，遗传和环境因素对西瓜裂果影响较大[5]，1994年赵尊练等[4]首次发现西瓜果皮硬度是果实耐贮性的主要内在因素，之后学者们分别从生理生化[6]、细胞水平[7-8]、分子水平[9]对西瓜果皮硬度进行研究，不同硬度西瓜的果皮含水量、纤维素、半纤维素、果胶含量、细胞排列方式、石细胞和外果皮层数相差较大，Liao 等[9]采用BSA-seq 技术和SNP 分析对西瓜果皮硬度进行研究，发现CIERF4基因与果皮硬度密切相关，该基因可能参与木质素的生物合成和细胞壁修饰，从而调控果皮硬度。而随着西瓜全基因组测序完成[10]，促进了转录组测序技术在西瓜果肉硬度[11-12]、果皮色泽[13]等方面的应用，但是，转录组测序技术在西瓜果皮硬度方面的研究鲜有报道。

本试验以河北省农林科学院经济作物研究所瓜类室创制的西瓜果皮高硬度的耐裂资源(901)4-1-1-M和低硬度易裂资源BSH为研究对象，通过对果皮硬度差异最大时期果皮的RNA-seq转录组测序，分析调节果皮硬度高的主要生物学过程和代谢通路，探讨果皮硬度差异大的基因表达调控机制，鉴定与果皮硬度相关的候选基因，从转录组水平解释西瓜果皮硬度差异的形成机制，旨在为耐贮运西瓜品种培育提供理论依据。

1.1 试验材料

供试西瓜为高硬度耐裂材料(901)4-1-1-M和低硬度易裂材料BSH，由河北省农林科学院经济作物研究所选育而成。西瓜开花后挂牌标记授粉日期，分别取授粉后第10，14，18，22，26，30天的西瓜，取样时选取果实大小一致的果实，使用质构仪测定西瓜阳面果皮硬度，3 次生物学重复；
选择检测部位两侧切割大小为3 cm×3 cm的果皮，厚度1 cm，每个生物学重复包含一个西瓜的2块西瓜果皮，3 次生物学重复，取样后用锡箔纸包裹后迅速放入液氮，保存于-80 ℃冰箱备用。选取果实成熟果皮硬度差异最大时期即授粉后第30天的西瓜果皮，用于后续的转录组测序和qRT-PCR验证。

1.2 果皮转录组分析

委托北京诺禾生物科技有限公司完成RNA提取、文库构建和转录组测序。以(901)4-1-1-M和BSH阳面果皮为转录组分析材料，利用Illumina HiSeqTM测序平台进行转录组测序分析，以西瓜97103基因组为参考基因组进行序列比对，对获得的差异基因进行GO富集分析、KEGG富集分析和差异基因功能分析。

1.2.1 RNA提取和文库构建 以(901)4-1-1-M和BSH阳面果皮为转录组分析材料，采用天根DP441多糖多酚植物总RNA提取试剂盒提取西瓜果皮总RNA，使用Agilent 2100 bioanalyzer精确检测RNA完整性和总量，RNA质控合格后通过Oligo(dT)磁珠富集带有polyA尾的mRNA，再用二价阳离子将得到的mRNA随机打断，以片段化的mRNA为模板、随机寡核苷酸为引物，在M-MuLV逆转录酶体系中合成cDNA第1条链，随后用RNaseH降解RNA链，并在DNA polymerase Ⅰ体系下，以dNTPs为原料合成cDNA第2条链。纯化后的双链cDNA经过末端修复、加A尾并连接测序接头，用AMPure XP beads筛选250～300 bp的cDNA进行PCR扩增，并再次使用AMPure XP beads纯化PCR产物，最终获得文库。

1.2.2 转录组测序 文库构建完成，库检合格后，利用Illumina HiSeqTM测序平台进行转录组测序分析，在测序的flow cell中加入4种荧光标记的dNTP、DNA聚合酶以及接头引物进行扩增，在每一个测序簇延伸互补链时，每加入一个被荧光标记的dNTP就能释放出相对应的荧光，测序仪通过捕获荧光信号，并通过计算机软件将光信号转化为测序峰，从而获得待测片段的序列信息。

1.2.3 测序结果分析 对测序获得的原始数据去除接头、无法确定的碱基信息、低质量的信息，获得高质量的测序数据(Clean data)，以西瓜97103基因组为参考基因组进行有参转录组比对分析，采用Subread软件对基因表达水平进行定量，利用DESeq 2软件进行差异基因分析，以|log2(FoldChange)|≥1 & padj≤0.05为标准筛选差异基因。采用ClusterProfiler软件对差异基因集进行GO功能富集分析和KEGG通路富集分析。

1.3 实时荧光定量qRT-PCR验证

选择与果皮硬度发育相关的9个差异基因进行实时荧光定量PCR(qRT-PCR)验证，采用天根DP441多糖多酚植物总RNA提取试剂盒提取(901)4-1-1-M和BSH西瓜果皮的RNA，再用HiFiScript gDNARemoval cDNA Syntheris Kit反转录试剂盒进行反转录合成cDNA后作为模板，最后采用ChemoHS qPCR Mix(None ROX)试剂，以西瓜的Actin基因作为内参基因，在CFX connectTM(BIO-RAD Applied Biosystems，Foster City，CA，USA)实时荧光定量PCR仪上进行qRT-PCR，每个基因设置3次生物学重复，按以下程序进行：95 ℃，15 min；
40个PCR循环(95 ℃，10 s；
58 ℃，30 s；
72 ℃，30 s)。荧光定量数据采用2-ΔΔCt相对定量计算，各引物序列见表1。

表1 qRT-PCR引物序列Tab.1 The primer sequence of qRT-PCR

2.1 质构仪测定结果

不同时期的(901)4-1-1-M和BSH西瓜果皮硬度试验结果表明(图1)，(901)4-1-1-M果皮硬度在各个时期均极显著高于BSH(P<0.01)，果实授粉后第18天(901)4-1-1-M果皮硬度达最大值89.02 N，BSH果皮硬度为42.46 N；
授粉后30 d(901)4-1-1-M果皮硬度为73.77 N，比BSH果皮硬度高47.21 N，此时果皮硬度差异最大。

A、B为在P<0.01条件下差异显著。A and B means the significant difference under P<0.01.

2.2 转录组测序质量评估

经转录组测序可知，6个样品共得到275 043 248条原始序列(Raw reads)，为保证信息分析质量，去除Raw reads里面含有带接头的、低质量的reads，得到滤后序列(Clean reads)272 077 872条，后续分析都基于Clean reads；
20(Q20)和30(Q30)分别表示碱基识别出错的概率为1%和0.1%，Q20和Q30值均达到90%以上；
GC百分率分布检查用于检测有无AT、GC 分离现象，本试验中西瓜果皮的GC百分率均在43%以上(表2)。

表2 6个样本转录组数据测序质量Tab.2 The sequencing quality of 6 sample transcriptome data

将Clean reads 与西瓜参考基因组比对，由比对结果可知(表3)，完全比对到参考基因组的reads占总reads数的71%以上，唯一比对位置序列比对率在70%以上，多个比对位置序列比对率在1.8%以下，外显子比对率在75%以上。据此可知，此次测序质量能够满足试验要求。

表3 与参考基因组比对结果Tab.3 Sequence comparison of samples with reference genome and genes

2.3 基因表达量分析

通过西瓜样本间相关性热图(图2)，对样本间基因表达量的生物学重复进行评估，样品间基因表达水平相关性是检验试验可靠性和样本选择是否合理的重要指标，相关系数R2越接近1，表明样品之间表达模式的相似度越高，由图3可知，6个西瓜样本的R2最小值为0.889，满足试验的生物学重复要求。

图2 6个西瓜样本间相关性分析Fig.2 The correlation analysis among 6 watermelon samples

2.4 基因表达数据分析

西瓜(901)4-1-1-M和BSH的共表达基因16 229个，(901)4-1-1-M特有的表达基因431个，BSH所特有的表达基因390个，采用DESeq 2软件对基因表达进行差异分析，以|log2(FoldChange)|≥1 & padj≤0.05为标准，筛选差异表达基因(DEG)。由图3可知，(901)4-1-1-M和BSH西瓜果皮中共筛选出1 085个差异表达基因，其中有555个上调表达基因，530个下调表达基因。

图3 差异基因火山图Fig.3 Volcanic map of differentially expressed genes

2.5 DEGs的GO富集分析

利用GO(Gene Ontology)数据库对获得的1 085个差异基因进行富集分析，1 085个差异基因被注释在生物学过程(Biological process)、细胞组分(Cellular component)及分子功能(Molecular function)三大功能方面，分别涉及254，41，213条功能分支。对富集程度前30的GO 条目作柱形图(图4)，差异基因显著富集在DNA复制(GO：0006260)(13个)、细胞壁(GO：0005618)(13个)、细胞边缘(GO：0071944)(16个)、外部封装结构(GO：0030312)(13个)、细胞外区(GO：0005576)(7个)、质外体(GO：0048046)(5个)、酶抑制剂活动(GO：0004857)(13个)、血红素结合(GO：0020037)(31个)、四吡咯结合(GO：0046906)(31个)、氧化还原酶活性，作用于配对供体(GO：0016705)(27个)、铁离子结合(GO：0005506)(26个)、转移酶活性，转移己糖基(GO：0016758)(30个)、果胶酯酶活性(GO：0030599)(8个)、蛋白质二聚化活性(GO：0046983)(27个)、酶调节剂活性(GO：0030234)(14个)和DNA结合(GO：0003700)(35个)方面。

2.6 DEGs的KEGG富集分析

利用KEGG数据库，1 085个差异基因被富集到81个代谢通路中，选择富集程度前20条代谢通路进行散点作图(图5)，差异基因主要在苯丙烷代谢(19个)，DNA复制(8个)，二萜生物代谢(5个)，角质、栓木和蜡的生物合成(5个)，戊糖、葡萄糖醛酸转换(8个)，α亚麻酸代谢(6个)，淀粉和蔗糖代谢(12个)，植物激素信号转导(17个)等通路富集，其中苯丙烷生物代谢(cmo00940)、DNA复制(cmo03030)、二萜生物合成(cmo00904)通路显著性富集(padj<0.05)。苯丙烷生物代谢途径是显著富集的代谢通路中注释到差异表达基因数目最多的通路，说明苯丙烷代谢物在西瓜果皮硬度差异方面存在重要作用。

图5 1 085个DEGs的前20条KEGG的功能分类Fig.5 The 20 tops of Enriched KEGG analysis for 1 085 DEGs

2.7 苯丙烷生物代谢中DEGs分析

对1 085个差异基因进行分析，(901)4-1-1-M和BSH在苯丙烷生物合成中存在19个差异基因(表4)。通过分析苯丙烷代谢通路发现(图6)，下游代谢物质木质素、对羟基苯酚木质素、愈创木酚木质素和5羟基愈创木酚木质素信号通路中富集到7个差异基因上调表达(Cla97C01G025380、Cla97C03G055890、Cla97C04G070210、Cla97C05G097050、Cla97C10G203590、Cla97C02G045070、Cla97C02G049940)，8个差异基因下调表达(Cla97C09G165820、Cla97C02G044950、Cla97C04G075800、Cla97C09G165630、Cla97C10G195660、Cla97C07G139710、Cla97C07G144540、Cla97C02G032910)，其中，下游基因Cla97C01G025380在(901)4-1-1-M果皮中基因表达量是BSH的60.24倍，而该基因在BSH中基因表达量最低，平均含量为0.056(表4)。

图6 苯丙烷生物合成途径中的差异基因表达Fig.6 DEGs in phenylpropanoid biosynthesis pathways

表4 苯丙烷生物合成路径差异基因表达量Tab.4 Expression analysis of genes related to phenylpropanoid biosynthesis

2.8 差异基因功能分析

从苯丙烷代谢路径中筛选出的19个差异基因，其中，6个过氧化氢酶基因(Cla97C02G049940、Cla97C05G097050、Cla97C03G055890、Cla97C02G045070、la97C04G070210)，4个转移酶基因(Cla97C02G033110、Cla97C09G165630、Cla97C07G144540、Cla97C02G032910)，3个糖基水解酶基因(Cla97C11G210660、Cla97C06G123590、Cla97C11G224410)，2个辅酶A连接酶基因(Cla97C09G165820、Cla97C10G195660)，脱氨酶、羟化酶、磷脂酶和桥小檗碱酶各1个，依次为Cla97C04G075800、Cla97C02G044950、Cla97C07G139710、Cla97C01G025380。

2.9 实时荧光验证

为了验证转录组数据的准确性，从富集程度最高的苯丙烷代谢路径中选择与果皮硬度相关的差异表达基因进行实时荧光验证(qRT-PCR)，对转录组测序数据(RNA-seq)和qRT-PCR的结果进行相关性分析，相关系数R2为0.791，证明RNA-seq结果试验数据可靠(图7)。

图7 RNA-seq和qRT-PCR相关性分析Fig.7 The correlation analysis results between RNA-seq and qRT-PCR

验证基因中包括与木质素代谢物质形成过程相关的9个差异基因，苯基丙氨酸在Cla97C04G075800、Cla97C02G044950、Cla97C09G165820、Cla97C10G195660基因作用下形成香豆酰辅酶A，再在Cla97C09G165630、Cla97C07G139710、Cla97C02G032910、Cla97C07G144540、Cla97C02G045070基因作用下最终形成不同种类的木质素；
qRT-PCR结果显示，下游基因Cla97C02G045070在(901)4-1-1-M中表达量较BSH提高7.13倍，这与RNA-seq结果中富集到苯丙烷通路的基因上调表达结果一致。

果实硬度变化是自然界常见的生理生化现象[14]，近年来，对果实硬度形成机理研究较多如苹果[15-16]、李子[17]、蓝莓[18]、西瓜[19-20]等，但对果皮硬度研究较少，Liao 等[9]采用BSA-seq 技术对西瓜果皮硬度进行研究，将果皮硬度基因精细定位到9 kb 的CIERF4(Cla97C10G187120)基因，该基因可能参与木质素的生物合成和细胞壁修饰从而调控果皮硬度。秦改花等[21]通过对石榴籽粒硬度与木质素、纤维素含量的相关性研究发现，籽粒木质素、纤维素含量与籽粒硬度显著相关，其可作为评价石榴籽粒硬度的指标。Qin等[22]通过对软硬籽石榴不同时期的种皮转录组、基因组和代谢组联合分析，筛选出与硬籽石榴中木质素单体芥子醇大量积累的相关候选基因POD(Pgr011171.1)和木质素单体松柏醇积累的相关候选基因F5H(Pgr013634.1)、CAD(Pgr022328.1，Pgr006310.1，Pgr006318.1和Pgr006319.1)。刘恋等[23]通过采用质构仪、生理生化指标测定和转录组分析发现，高硬度金桔果皮中木质素含量显著高于低硬度金桔，转录组测序结果表明，苯丙烷生物合成代谢引起的木质素含量差异导致了3个品种金桔果皮韧性存在差异。

本试验采用RNA-seq 技术分析了果皮硬度差异显著的2份资源(901)4-1-1-M和BSH的果皮转录本差异，与果皮低硬度易裂资源BSH相比，(901)4-1-1-M中555个基因上调表达，530个基因下调表达；
GO富集分析发现，1 085个差异表达基因显著富集在细胞壁、细胞外围、胞外区、四吡咯骨架、氧化还原酶活性、果胶酯酶活性等功能类别；
KEGG富集分析发现，19个差异表达基因显著富集在苯丙烷代谢通路，其中下游基因Cla97C01G025380在(901)4-1-1-M中显著上调表达，导致木质素单体芥子醇、松柏醇、5羟基松柏醇形成，进而在Cla97C03G055890、Cla97C04G070210、Cla97C05G097050、Cla97C10G203590、Cla97C02G045070、Cla97C02G049940基因作用下形成木质素、5羟基愈创木酚木质素、愈创木酚木质素和对羟基苯酚木质素代谢物。研究结果与Liao 等[9]发现的西瓜果皮硬度与木质素的生物合成路径相关的结论一致；
与Qin等[22]发现的硬籽石榴中木质素单体芥子醇大量积累和木质素单体松柏醇积累的研究结果一致；
与刘恋等[23]发现的金桔果皮硬度与木质素正相关，果皮韧性差异富集在苯丙烷代谢通路的结论一致。

因此，推测苯丙烷代谢通路中与木质素代谢物形成相关的差异表达基因可能是造成西瓜果皮硬度差异的主要原因，但与Liao 等[9]发现的CIERF4(Cla97C10G187120)基因参与木质素代谢结果不一致，主要原因可能是由于西瓜材料不同，但本试验中发现，下游基因Cla97C01G025380在BSH中基因表达量低，影响芥子醇、松柏醇、5羟基松柏醇代谢物形成，这与Qin等[22]研究结果一致，本试验中发现的西瓜果皮硬度富集到苯丙烷代谢路径中木质素与刘恋等[23]的研究结果一致。

本研究通过对西瓜果皮硬度差异显著的2种资源(901)4-1-1-M和BSH成熟期的果皮进行转录组分析，发现了一些与西瓜果皮硬度密切相关的代谢路径，尤其是苯丙烷代谢途径最终形成的不同种类的木质素在西瓜果皮硬度形成中起重要作用，其中，19个差异基因被富集到该通路，包括Cla97C04G075800、Cla97C02G044950、Cla97C10G195660、Cla97C09G165820、Cla97C09G165630、Cla97C02G032910、Cla97C07G144540、Cla97C02G045070等基因。该研究从转录水平解析了(901)4-1-1-M和BSH西瓜果皮硬度存在差异的分子机制，为挖掘西瓜果皮硬度相关的关键基因提供了理论基础。