不同果肉厚度苦瓜转录组比较分析

裘波音，郑 旋，李大忠，林 珲，张前荣，刘建汀，朱海生，温庆放

(1.福建省蔬菜遗传育种重点实验室,福州 350013；2.福建省农业科学院作物研究所,福州 350013；3.福建省蔬菜工程技术研究中心,福州 350013；4.福建省种子总站(福建种子质量检验站)，福州 350003)

【研究意义】苦瓜(MomordicacharantiaL.)属于葫芦科苦瓜属一年生蔓性草本植物，起源于热带，喜光喜湿，对肥料要求较高，花、果期5—10月[1]。苦瓜因含有罗汉果苷和苦瓜素等成分而呈苦味，是一种重要的食药兼用园艺作物，具有重要经济价值。果肉厚度是瓜果类作物重要的品质性状，也是产量的决定性因素，即果肉越厚，果实的可食用部分就越大[2]。果肉发育是果实发育的重要阶段，也是关系果实整体品质的核心因子[3]。果实发育过程受一系列基因表达引起的生理生化变化调控，最终形成不同的颜色、质地、香气、营养成分等，导致产量及品质上的差异[4-5]。研究果肉发育一方面有利于深化植物发育相关机制的探讨，另一方面也对推进理论成果转化应用起到积极作用。【前人研究进展】高通量转录组测序作为后基因时代率先发展起来的技术，已经在生物学领域得到了广泛的应用。在蔬菜相关研究中，转录组测序技术主要用于分子标记开发，如刘峰等[6]基于5份辣椒的转录组数据筛选获得18 159个SNP，且选取的15对验证引物中有8对表现出多态性；林珲等[7]在24份花椰菜中挖掘出10 715个转录组SSR位点，分析获得具有潜在高多态性的引物6119对；刘建汀等[8]在印度南瓜11 871个Unigenes中发现其中5391个含有SSR位点，且含有至少2个SSR位点的有319个；遗传育种和进化分析，如Koenig等[9]利用转录组分析了番茄栽培品种和5个近缘野生品种的基因序列和表达差异，发现序列变异方面有50多个基因出现了正向选择，而在表达水平方面，有数千个基因出现了变化，其中很多包括果色、耐热、耐盐等在内的相关基因其表达差异是由选择压力造成的，表明人工驯化在栽培植物的形态、生理和生活史方面产生了重要影响；差异表达基因分析和功能注释，如王明霞等[10]对安徽乌菜与普通白菜共3份材料进行测序，其中200个KEGG注释基因中参与植物病原菌互作通路的占11%，位居第一，参与植物激素信号转导通路和类苯基丙烷生物合成通路的占9.5%，排第二；吴萍等[11]对来源于6个产地的不同种白芷进行分析共获得15 939个差异基因，其功能主要与细胞组分、催化活性、代谢过程等相关，其中在川白芷中上调最高的10个基因多与蛋白质加工、激素信号传导等相关，而下调最多的10个基因则与RNA运输、果糖代谢等相关；功能基因挖掘，如基于转录水平获得的BvM14-SAMS2基因在甜菜M14品系根中表达最高，并受盐胁迫诱导表达[12]；而WRKY家族各成员在甜菜不同组织间表达差异大，其中BvWRKY33在热胁迫和盐胁迫下均呈上调表达[13]，表明上述基因对甜菜逆境生理调控具有重要作用。【本研究切入点】近年来，苦瓜营养价值被不断开发和利用，促使其生产在全世界逐渐流行，这对苦瓜种质创新提出了挑战。了解不同果肉厚度苦瓜果实成熟过程中的基因表达情况，有利于探究果实发育调控机理，推动良种培育，满足人类对其持续消费的需求。【拟解决的关键问题】本研究立足果肉厚度，对不同苦瓜种质进行统计分析，筛选果肉厚度差异显著材料，再利用转录组测序技术对其进行比较分析，挖掘差异表达基因，预测转录因子，并对差异表达基因进行功能注释、代谢通路富集等分析，为进一步揭示苦瓜果肉发育分子调控机制提供理论依据。

1 材料与方法

1.1 试验材料

70份苦瓜品系由福建省农科院作物研究所蔬菜研究室提供，于2020—2021年春季种植于福州市试验基地(119.28°E, 26.08°N)。

1.2 果肉厚度差异苦瓜筛选

于成熟期采样，每个品系采长势相似的4个瓜，在离瓜蒂2/3处横切，用游标卡尺测量横径和瓤腔宽度，计算果肉厚度，公式如下：果肉厚度=(横径-瓤腔宽度)×0.5。根据果肉厚度的显著性差异筛选厚肉和薄肉材料各1份进行后续试验。

1.3 RNA提取与检测

于样品成熟期对厚肉型苦瓜和薄肉型苦瓜各1份进行取样，每份材料取9个瓜，设3个生物学重复，即每个重复为3个瓜混样。采摘后立即置于液氮中保存，用于总RNA的提取。利用OmniPlant RNA(DNase I)试剂盒(康为世纪，北京)进行RNA提取；1%琼脂糖凝胶电泳分析样品RNA完整性及是否存在DNA污染；利用NanoPhotometer spectrophotometer检测RNA纯度(OD260/280及OD260/230比值)；利用Agilent 2100 bioanalyzer精确检测RNA完整性。样品送福建华诺通达公司进行转录组测序。

1.4 文库构建与质检

使用Illumina的NEBNext®UltraTM RNA Library Prep Kit对总RNA进行建库。文库构建完成后，先使用Qubit2.0 Fluorometer进行初步定量，随后使用Agilent 2100 bioanalyzer对文库的insert size进行检测，qRT-PCR对文库有效浓度进行准确定量(文库有效浓度高于2 nmol/L)，以保证文库质量。

1.5 上机测序及功能注释

采用Illumina HiSeq2500进行测序，对测序所获得的原始数据，按照华诺通达公司的数据处理方法进行过滤，去除带接头(adapter)、含N(N表示无法确定碱基信息)和低质量reads(Qphred≤20的碱基数占整个read长度的50%以上的reads)，再对clean data进行Q20，Q30和GC含量计算。后续所有分析均是基于clean data进行的高质量分析。

本研究采用有参转录组测序(Momordica charantia ASM199503v1，https://ftp.ncbi.nlm.nih.gov/genomes/all/GCF/001/995/035/GCF_001995035.1_ASM199503v1/)，使用HISAT2 v2.0.5构建参考基因组的索引，并将配对末端clean reads与参照基因组比对。通过BlastX软件，与NCBI的non-redundant(NR)、UniProt/Swiss-Prot、Gene Ontology(GO)、Kyoto Encyclopedia of Genes and Gnomes(KEGG)等数据库进行比对，以获得其功能注释信息。

1.6 差异基因鉴定及富集分析

通过FPKM法来计算基因的表达量，即每百万fragments中来自某一基因每千碱基长度的fragments数目。差异基因筛选基于“有生物学重复”类型，利用DESeq2软件和负二项分布的P-value计算模型，以“|log2(Fold Change)| > 0，padj <0.05”为标准进行。用BLAST2GO软件进行GO注释，用KOBAST软件进行KEGG注释。并采用clusterProfiler软件对差异基因集进行GO功能富集分析和KEGG通路富集分析。

1.7 数据分析

数据采用SPSS 21.0进行统计分析，采用WPS Office进行作图。

2 结果与分析

2.1 不同果肉厚度苦瓜筛选

如表1所示，测定的70个品系材料果肉厚度分布范围较广，最大值为12.85 mm，最小值为8.65 mm，平均厚度为9.84 mm，整体呈正态分布，其中大部分厚度在9～11 mm，小于9 mm的有7个，大于11 mm的有3个。基于一般线性模型的单因素分析进行多重比较，方法选用Tamhane’s T2(P≤0.01)，发现LX1-3(12.85 mm)和ZK54(8.83 mm)果肉厚度具有极显著差异，前者是后者的1.46倍，因此，将LX1-3作为厚肉型材料，ZK54作为薄肉型材料进行后续分析。

表1 苦瓜不同果肉厚度统计Table 1 Summary of flesh thickness of different bitter gourds

2.2 转录组测序统计

经转录组测序分析统计，LX1-3共计得到241 306 304条初始序列，过滤后得到234 825 400条干净序列，17.61 G干净碱基对，平均Q30值为92.74%；比对到参考基因组上的序列平均为34 455 982条，占总序列数的88.07%，其中，单一比对序列为34 027 357条，占比86.98%，多重比对序列为428 625条，占比1.10%，GC含量为46.54%；ZK54共计得到250 983 052条初始序列，过滤后得到241 318 020条干净序列，18.1 G干净碱基对，平均Q30值为92.26%；比对到参考基因组上的序列平均为35 079 431条，占总序列数的87.23%，其中，单一比对序列为34 607 181条，占比86.06%，多重比对序列为472 250条，占比1.17%，GC含量为46.75%(表2)。

表2 样品转录组数据统计Table 2 Summary of transcriptome sequencing results of different samples

对序列长度统计分析(图1)显示，长度小于200 bp的有54条，占0.28%；长度为200～500 bp的有684条，占3.49%；长度为501～1000 bp的有3672条，占18.73%；长度为1001～2000 bp的合计共有8578条，占43.76%；长度为2001～3000 bp的有4169条，占21.27%；长度为3001～4000的有1555条，占7.93%；长度超过4000 bp的有890条，占4.54%。总体呈“两头低中间高”的趋势，数量最多的序列长度集中在1001～2000 bp。

图1 苦瓜转录组序列的长度分布Fig.1 The length distribution of unigenes of transcriptome in bitter gourd

2.3 差异基因鉴定及转录因子预测

经筛选共得到1492个差异基因，其中796个基因上调表达，696个基因下调表达，分别占差异基因数量的53.35%和46.65%，从图2中可以看岀,2组材料之间的基因表达水平呈显著差异。本次测序共预测到598个转录因子，分属于216类，│log2(Fold Change)│值在0.62～11.24。其中，Pkinase数目最多(45个，7.53%)，其次是p450(29个，4.85%)，接着为Pkinase_Tyr(21个，3.51%)，2OG-FeII_Oxy(17个，2.84%)，UDPGT(16个，2.68% )，Myb_DNA-binding(14个，2.34%)，AP2(12个，2.01%)和HLH(10个，1.67%)，其他家族共有基因434个，占72.58%(图3)。

图2 不同果肉厚度苦瓜差异基因Fig.2 Differential gene of bitter gourd with different sarcocarp thickness

图3 苦瓜序列的转录因子预测 Fig.3 Transcription factor prediction of Unigenes in bitter gourd

2.4 GO功能富集分析

共有1452个差异基因被注释到了GO数据库中，得到616条功能注释，共分为三大功能，分别是生物过程、细胞组分和分子功能，依次得到308、93和215条功能注释，占比分别为50.00%、15.10%、34.90%。每项功能中重要注释(P<0.05)如表3所示。生物过程含DEGs 556个，共获得308条注释。最主要的11条功能注释共包含151个DEGs，占27.16%。从差异基因数量上可看出，“碳水化合物代谢过程”条目居于第一位，含47个DEGs，其次为“细胞碳水化合物代谢过程”和“多糖代谢过程”，分别含18和13个DEGs，条目“细胞葡聚糖代谢过程”，“葡聚糖代谢过程”，“细胞多糖代谢过程”和“抗氧化反应”并居第四，分别含12个DEGs。从差异基因的上下调变化上得出，11个条目中上调基因数目多于下调基因数的有7条，且其中6条位于前六位。细胞组分含DEGs 180个，功能注释93条。最主要的前9条功能注释(P<0.05)共含75个DEGs，占41.67%。其中，“细胞周边”条目含DEGs数目最多，为16个，其次为“细胞壁”和“外部封装结构”，均含有11个。9个条目共计含有上调基因55个，下调基因20个，上调基因数目是下调基因的2.275倍。分子功能含DEGs 716个，功能注释215条。前27条主要注释条目含有587个DEGs，占81.98%。“血红素结合”和“四吡咯结合”最多，均含有43个DEGs，其次为“转移酶的活性/转移糖基” “DNA结合转录因子活性” “转录因子活性”，各含有41个DEGs。27个条目共计含有上调基因328个，下调基因259个，上调基因数目多于下调基因。从上述结果看出，生物过程，细胞成分和分子功能三方面的上调基因数目均大于下调基因数，说明果肉厚度的发育受上调基因的表达影响较大。

表3 不同果肉厚度苦瓜差异表达基因的GO富集分析Table 3 GO enrichment analysis of differentially expressed genes of bitter gourd with different sarcocarp thickness

续表3 Continued table 3

2.5 KEGG通路富集分析

共有308个差异基因富集到了KEGG通路上，包括188个上调基因和120个下调基因，结果如图4所示。上调基因有178个富集到20条主要的代谢通路。富集因子最高的为“卟啉和叶绿素代谢”通路，其次为“光合作用生物的碳固定作用”和“光合作用”通路。而富集数目最多的通路为“核糖体”，其次为“碳代谢”，两者所占DEGs比列分别为12.77%和9.04%。下调基因有111个富集到20条主要的代谢通路。其中，富集因子最高的通路为“苯丙素生物合成”，其次为“二萜生物合成”通路。富集下调差异基因最多的通路为“苯丙素生物合成”，占差异基因比例为10.00%，其次为“淀粉和蔗糖代谢”“氨基酸的生物合成”和“碳代谢”，三者数量相当，所占比列为7.5%。从富集DEGs的数量角度也可以看出，厚肉苦瓜材料的代谢途径中上调基因数量高于下调基因数量，是下调基因数量的1.56倍，说明果肉厚度发育可能主要受上调基因表达调控。

圆点的位置代表代谢途径，圆点的大小代表基因的数量The position of the dot represents the enriched pathway，the size of the dot indicates the number of differential genes图4 差异表达基因的KEGG分析Fig.4 The KEGG analysis of the differentially expressed genes

3 讨 论

苦瓜是重要的蔬菜作物，其果实生长发育受细胞分裂、激素水平、信号传导、糖碳基础代谢、环境条件等多种因素的影响。果肉厚度作为一个重要的性状，对其产量和质量起决定性作用。立足果肉厚度的早期研究报道不多，何晓明[14]、徐强[15]等对黄瓜产量和果实性状进行调查，得出果肉厚度对黄瓜产量影响重大，基于相关分析也得出果实的鲜重与果肉厚度的遗传呈正相关，杨涛等[16]也发现辣椒果肉厚度与果实横径、单果重也极显著正相关。本研究通过对70份苦瓜的果肉厚度测定分析，发现不同材料的果肉厚度差异较大，最大值与最小值差值占比平均果肉厚度超四成，说明果肉厚度受基因型调控表现出种系间差异，进而可能对产量及品质造成影响。LX1-3和ZK54两者果肉厚度具有极显著差异，前者是后者的1.46倍，但其它生长性状基本相似，因此将两者作为后续转录组测序的研究材料对保证结果真实性和准确性具有有利的一面，对深入研究其次生代谢、产量形成、品质进化等具有重要意义。

转录组测序是目前生物学研究中的主流技术，通过分析不同材料在不同时期或环境下的转录组差异，有助于了解它们的差异表达在不同生物学过程中的作用[17]。早期研究发现，果实品质在成熟过程中会发生一系列生理生化变化，这些变化与各种基因的表达调控有关，如黄瓜基因Csa2GO58670与两亲本(厚果肉和薄果肉)果肉厚度均显著相关，果实发育不同时间该基因在两个品种中相对表达量存在差异[18]。本研究基于不同果肉厚度苦瓜进行转录组分析，发现多个差异表达基因，其中上调基因数目明显多于下调基因数，表明其果肉发育可能受上调基因影响较大。作为转录过程的重要元件，转录因子能与5’端上游特定序列专一性结合，从而保证目的基因以特定的强度在特定的时间与空间表达，在真核生物中，转录因子能与RNA聚合酶Ⅱ形成转录起始复合体，共同协助其完成复杂的转录过程。 本研究预测得到598个转录因子，但分属于多个家族，说明参与苦瓜果肉发育调控的转录因子种类繁多，再次证实植物发育是一个极其复杂的过程，最终果实的质量和品质是多种因素综合调控的结果。

GO和KEGG分析是基因注释的重要手段，GO的功能分类有利于了解基因功能，而KEGG富集分析除了可以推测基因功能，还可以研究基因在不同代谢通路中的位置及作用[19]。本研究GO注释表明，生物过程在三大功能中占主要地位，深入挖掘发现糖类代谢相关注释所含差异基因最多且上调基因数均多余下调基因数，同时，KEGG富集分析也发现，上调基因与叶绿素、碳固定以及光合作用相关途径富集因子最高，而核糖体及碳代谢途径富集的差异基因数量最多，表明在维持苦瓜果肉正常发育中，基础糖碳代谢和光合作用以及源库关系在整个生物过程中至关重要。光合作用是植物进行物质生产的基本生理过程，植物绝大部分干物质是由光合作用提供的，光合作用的重要产物及其运输的主要形式是蔗糖[20]，光合组织中蔗糖不仅影响着光合同化产物的分配运输，其含量的多少直接影响籽实中糖分含量，对糖类整体代谢、果实的膨大与生长速率均造成干扰，进而影响作物产量与品质[21-23]。同时，各类酶在保障叶绿体内光合碳代谢、蔗糖运输、多糖积累平衡等生物过程中起着重要的调节作用，并且与多种调控作物生长和产量的指标密切相关[24-25]。本研究基于转录组测序发现参与苦瓜果肉发育生物过程的重要调控因子与早期研究结果类似，再次证明糖碳代谢是保证作物生长发育的关键。早期研究表明，源库器官的数量及其机能的相互协调，是作物生产力的重要决定因素。在本次试验中，厚肉型苦瓜相较薄肉型苦瓜在糖代谢和光合作用相关基因表达水平更高，表明其具有较高的“源”性能，为多糖形成提供基础，而“淀粉和蔗糖代谢”基因下调，表明“库”性能在化合物分配时将较多的基底物质倾向于葡聚糖等多糖的形成，因此推测，改善光合作用性能、调节源库平衡可能是提高苦瓜产量及品质的重要途径。

众所周知，生长素、赤霉素、细胞分裂素等植物内源激素对果实的发育和成熟也发挥着重要作用，不同的激素在不同品种以及同一品种的不同发育阶段作用不一。IAA含量的提高能促进发育早期的番茄子房细胞的分裂和膨大以及诱导果实发育[26]；GA种类众多，在幼果中含量较高，能刺激果实和种子发育[27]；茄子子房(果实)发育前期的IAA、ZR含量明显高于发育后期，其中单性结实品系的IAA、ABA含量和ABA/(IAA+GA4+ZR)均高于非单性结实品系[28]；黄瓜Q30(CsSUN/CsLNG1)的BR/ABA与GA3/ABA在开花前6 d显著低于其近等基因系材料[27]。本研究对下调基因进行KEGG分析发现，“苯丙素生物合成”是富集因子最高且基因数最多的通路，这个通路主要涉及一些过氧化物酶，其次为“二萜生物合成”，以赤霉素加双氧酶为主。植物体内的过氧化物酶对生长素有破环作用，对调控生长素水平，维持生长素库源的平衡非常关键。而赤霉素加双氧酶参与赤霉素代谢，对维持植物体内活性赤霉素的水平及含量具有重要作用。上述结果表明植物激素参与苦瓜果肉发育，再次证实激素平衡在果实发育中的重要地位。

细胞壁组分和代谢酶对植物生长发育具有重要影响，如基于甜瓜和黄瓜材料在成熟过程的转录组测序数据，发现了果实膨大和发育与细胞壁相关基因包括扩张蛋白基因(Csa6G014540)、木葡聚糖内转葡糖基酶基因(MelO3C014468，Csa1G422480)等关系密切[29-30]；而随着果实的不断成熟，砂梨中β-半乳糖苷酶基因PpGal1和PpGal4表达量呈逐渐上升趋势[31]，鳄梨中PaGal2表达量增高但PaGal4则不表达[32]。本次研究发现细胞组分条目中涉及细胞周边、细胞壁以及外部封装结构的基因如果胶酯酶基因、木葡聚糖内糖基转移酶/水解酶基因等数目占比较大，表明细胞壁相关结构在苦瓜果肉发育中也起到了较大的作用。

4 结 论

通过研究获得的果肉厚度差异苦瓜种质为培育高产优质苦瓜提供潜力材料；高通量转录组测序挖掘的差异基因以及相关功能注释和富集分析结果，为苦瓜果实发育分子调控机制的探讨提供理论基础。