Transcriptome analysis of wax gourd grafted with pumpkin to improve resistance to Fusarium wilt

冬瓜品种选用海南主栽品种‘铁柱2号’，南瓜品种选用‘海砧1号’，由海南农业科学院蔬菜研究所提供。

育苗地点设在海南省农业科学院澄迈基地。2023年5月24日浸种冬瓜种子用于培养自接冬瓜的砧木，6月2日浸种冬瓜种子用于培养自接冬瓜和嫁接冬瓜的接穗，6月6日浸种南瓜种子用于培养自接南瓜和嫁接冬瓜的砧木，6月9日浸种南瓜种子用于培养自接南瓜的接穗，浸种时间均为12 h，砧木种子露白长度达到种子本身1/2播于60孔穴盘，接穗种子露白长度达到种子本身1/2播于托盘中。待接穗长至两片子叶平展，砧木第一片真叶成梭状时采用插接法嫁接与自接^[12]，得到嫁接冬瓜、自接冬瓜和自接南瓜3种幼苗，分别用GW、SW和SP表示。3种幼苗培养至两叶一心期备用。

实验菌种由广西农业科学院蔬菜研究所提供。该菌种为尖孢镰刀菌冬瓜专化型（Fusarium oxysporum Schl. F. sp. benincasae）生理小种1号（后用FOB表示），从冬瓜‘铁柱2号’中分离获得。

PDA与PDB培养基的制备参照^{[13 − 14]}。FOB于PDA培养基中培养7 d后，挑取菌丝转至PDB培养基，置25 ℃、150 r·min⁻¹转速的恒温振荡器培养3 d，脱脂棉过滤得到孢子悬浮液，测定得到孢子浓度为2.75×10⁷ CFU·g⁻¹，按1∶4的比例稀释孢子悬浮液得稀释液，备用。

土壤采自澄迈县瑞溪镇罗浮村，基本理化性质为：pH 5.94、有机质5.1 g·kg⁻¹、碱解氮28.1 mg·kg⁻¹、速效磷103.3 mg·kg⁻¹与速效钾129.6 mg·kg⁻¹。拣出实验土壤杂物，过0.18 mm筛子，灭菌3次，每次2 h^[15]，风干后得灭菌土壤，备用。

将孢子稀释液和水分别按土壤∶液体= 6 g∶1 ml混合拌匀，得接菌土壤。测土壤FOB浓度，为8.0×10⁵ CFU·g⁻¹和0 CFU·g⁻¹（用F和0表示）。装接菌土壤，每钵120 g，备用。

实验设3类幼苗，分别为GW、SW与SP，2种浓度接菌土壤，分别为F和0。选长势一致的3种幼苗，清洗根系后，定植在装有接菌土壤的方形塑料钵中，每钵1株。6种接菌处理的幼苗分别用GWF、GW0、SWF、SW0、SPF、SP0表示。GWF、SWF和SPF处理幼苗各25株，GW0、SW0和SP0处理幼苗各10株。

定植后将幼苗置光照培养箱培养，培养条件：光照周期12 h/12 h，温度26 ℃/21 ℃（昼/夜），光强为2000 Lux，培养前2 d无需浇水，后每2 d补充植株所蒸发的水分，实验过程中不施肥。

课题组前期通过不同浓度病菌的接菌实验发现，接菌土壤FOB浓度为2.8 × 10⁶ CFU·g⁻¹时，自接冬瓜幼苗第6 d开始出现枯萎的症状，接菌土壤FOB浓度为8.0 × 10⁵ CFU·g⁻¹时，第10 d开始出现枯萎的症状，而嫁接冬瓜均未出现枯萎的症状^[11]，故本研究选择病菌入侵阶段（第5 d）开展幼苗根系与叶片的转录组分析。

培养5 d后对GWF、GW0、SWF和SW0幼苗叶片（分别用GWFL、GWL、SWFL和SWL表示）和GWF、GW0、SPF和SP0幼苗根部（分别用GWFR、GWR、SPFR和SPR表示）进行采样，样本均设4个重复，每组4个生物学重复命名为 1、2、3、4。先用自来水冲洗第3～4片叶片的正反面（根部用自来水冲洗根部表面），再用去离子水再次冲洗、立即用滤纸吸干残留的水分，用灭菌后的剪刀剪下叶片或根部后置于锡箔纸中，去除叶片主脉后剪碎（根部首尾各1 cm，将中间根部剪碎），置于2 mL冻存管中，随后立即于液氮中速冻，后转到−80 ℃超低温冰箱中保存，随后立即将样品交由上海美吉生物医药科技有限公司完成RNA-seq测序工作。

使用fastp v 0.19.5软件（https://github.com/OpenGene/fastp）对原始数据进行过滤，过滤掉带接头、首尾端低质量、碱基N含量过高以及舍弃去adapter及质量修剪后长度小于30 bp的序列，从而得到高质量的测序数据（clean data）以保证后续分析的顺利进行。

使用HiSat2v2.1.0软件（http://ccb.jhu.edu/software/hisat2/index.shtml）将质控后的clean data （reads），与参考基因组比对获得用于后续转录本组装、表达量计算等的mapped data （reads），叶片参考基因来源: Benincasa_hispida （https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_009727055.1/），根系参考基因来源：Cucurbita_moschata （https://www.ncbi.nlm.nih.gov/datasets/genome/GCF_002738365.1/），同时对该次转录组测序的比对结果进行质量评估，主要包括测序饱和度、基因覆盖度、Reads在参考基因组不同区域分布以及Reads在不同染色体分布分析。

获得基因的Read Counts数后，使用DESeq2 v 1.10.1软件（http://biocon-ductor.org/pac-kages/stats/bioc/DESeq2）进行样本间基因的表达差异分析，将|log₂FC| ≥ 1且FDR （False Discovery Rate）<0.05的基因定义为差异表达基因（differentially expressed genes, DEGs）。

利用GO数据库，将基因按照生物学过程（Biological Process, BP）、细胞组分（Cellular Component, CC）和分子功能（Molecular Function, MF）进行分类，对差异表达基因的主要功能进行GO注释。

利用KEGG数据库，将基因按照pathway通路或行使的功能进行分类，对差异表达基因进行KEGG注释。

利用植物抗性基因数据库（PRGdb;http://prgdb.org/prgdb4/），将差异表达基因的蛋白序列提交到DRAGO 3工具中预测R基因的类别及结构域，进一步筛选抗病基因。

此外，GO功能富集和KEGG功能富集分析分别通过软件Goatools v 0.6.5 （https://git-hub.com/tanghaibao/GOatool）和Python scipy （https://scipy.org/install/）进行，使用 4种多重检验方法 （BH, BY, holm 和 bonferroni）对p值进行了校正，当p值≤0.05时，认为GO与KEGG功能存在显著富集情况。

选取叶片4个、根部2个差异表达基因进行实时荧光定量PCR，验证转录组结果的准确性。

以F-box protein基因作为叶片的内参基因^[16]，以Actin基因作为根部的内参基因^[17]进行qRT-PCR反应。实验设3个重复，使用2^−∆∆ct法计算样本基因相对表达量，使用Origin 2024软件进行绘图。利用NCBI中的Pick Primers设计引物（表1），使用Blast进行引物特异性检验，将设计好的引物交于擎科生物科技有限公司海南分公司合成。

16个叶片样品测序共获得原始序列700 019 546条，原始数据经过滤后获得高质量的reads数共649 983 334条，共获得103.28 Gb Clean Data，各样品Clean Data均达到6.04 Gb以上。Q20碱基百分比含量的平均值为97.73％，Q30碱基百分比含量的平均值为95.88％，GC的含量的平均值为45.89％（表2）。16个根部样品测序共获得原始序列796 045 970条，原始数据经过滤后获得高质量的reads数共788 820 340条，共获得118.44 Gb Clean Data，各样品Clean Data均达到6.29 Gb以上。Q20碱基百分比含量的平均值为97.90％，Q30碱基百分比含量的平均值为93.75％，GC的含量的平均值为45.93％（表3），说明文库构建质量较高，可满足下一步的数据分析。

将叶片组织样品比对到冬瓜参考基因组后，利用HiSat2软件分别将各样品的Clean Reads与指定的参考基因组进行序列比对，比对率从96.65%到97.64%不等。该比对率较高，可以保证后续分析的准确性。与参考基因组相比，最终鉴定得到检测到表达基因共22 808个，其中已知基因有21 121个，新基因1 687个。将根部组织样品比对到南瓜参考基因组后，利用HiSat2软件分别将各样品的Clean Reads与指定的参考基因组进行序列比对，比对率从72.43%到86.66%不等。该比对率较高，可以保证后续分析的准确性。与参考基因组相比，最终鉴定得到检测到表达基因共29 436个，其中已知基因有28 566个，新基因870个。

为了检验样本的合理性，分别对叶片与根部各样品间所有基因表达量进行相关性分析。如图1所示，叶片生物学重复样品间Pearson相关系数的平方值（R²）均在0.8以上（图1−A），根部生物学重复样品间Pearson相关系数的平方值（R²）均在0.9以上（图1−B），说明二者样品之间生物学重复性良好，可以满足后续分析工作。

Figure 1.  Sample correlation heat map of leaf（A）and root（B）

将|log₂FC| ≥ 1且FDR < 0.05的基因定义为差异表达基因。在叶片当中，GWFL vs SWFL组共鉴定出4 310个DEGs，其中包含1 956个上调DEGs和2 354个下调DEGs；GWFL vs GWL组共鉴定出6 808个DEGs，其中包含3 130个上调DEGs和3 678个下调DEGs；SWFL vs SWL组共鉴定出3 759个DEGs，其中包含2 290个上调DEGs和1 469个下调DEGs。在根部当中，GWFR vs SPFR组共鉴定出56个DEGs，其中包含49个上调DEGs和7个下调DEGs；GWFR vs GWR组共鉴定出803个DEGs，其中包含403个上调DEGs和401个下调DEGs；SPFR vs SPR组共鉴定出242个DEGs，其中包含115个上调DEGs和127个下调DEGs。

为进一步筛选关键抗枯萎病基因，本研究以接菌后在嫁接冬瓜中上调，且在自接冬瓜或南瓜中下调的差异表达基因定义为潜在抗病基因^[18]，并对其进行进一步分析。对叶片GWFL vs GWL组的3 130个上调DEGs、GWFL vs SWFL组的1 956个上调DEGs与SWFL vs SWL组的1469个下调DEGs进行Venn分析，筛选出284个叶片潜在抗病基因（图2−A）。对根部GWFR vs SPFR组的49个上调DEGs、GWFR vs GWR组的402个上调DEGs与SPFR vs SPR组的127个下调DEGs进行Venn分析，筛选出2个根部潜在抗病基因（图2−B）。

Figure 2.  Venn map of DEGs in leaf（A）and root（B）

对所筛选出来的284叶片关键抗病基因进行GO功能分析，在生物过程中，叶片的关键抗病基因主要富集在细胞过程、代谢过程与生物调节；在细胞组分中，叶片的关键抗病基因主要富集在细胞组分、膜组分和细胞器；在分子功能中，叶片的关键抗病基因主要富集在催化活性与结合（图3−A）。

Figure 3.  GO functional classification of potential disease resistance genes in leaf（A）and root（B）

同时对所筛选出来的2个根部关键抗病基因进行GO功能分析，在生物过程中，根部的关键抗病基因主要富集在代谢过程、细胞过程和应激反应；在细胞组分中，根部的关键抗病基因主要富集在含蛋白质复合物、膜组分、细胞器组分、细胞膜和细胞组分；在分子功能中，根部的关键抗病基因主要富集在催化活性与结合（图3−B）。

叶片关键抗病基因在丙氨酸、天门冬氨酸和谷氨酸代谢（Alanine, aspartate and glutamate metabolism）、苯丙烷生物合成（Phenylpropanoid biosynthesis）、角质、亚蜡和蜡的生物合成（Cutin, suberine and wax biosynthesis）与氮代谢（Nitrogen metabolism）等通路显著富集（图4−A）。

Figure 4.  KEGG enrichment analysis of potential disease resistance genes in leaf（A）and root（B）

根部的关键抗病基因在芪类、二芳基庚烷类和姜酚生物合成（Stilbenoid, diarylheptanoid and gingerol biosynthesis）、光合作用天线蛋白（photosynthesis-antenna proteins）、类黄酮生物合成（Flavonoid biosynthesis）与苯丙烷生物合成（Phenylpropanoid biosynthesis）等通路显著富集（图4−B）。

叶片中丙氨酸、天门冬氨酸和谷氨酸代谢（Alanine, aspartate and glutamate metabolism）通路涉及到3个基因，分别是LOC120076313、LOC120089431与LOC120083554。苯丙烷生物合成（Phenylpropanoid biosynthesis）涉及到5个基因，分别是LOC120067870、LOC120073918、LOC120081818、LOC120085310与LOC120091171。角质、亚蜡和蜡生物合成（Cutin, suberine and wax biosynthesis）涉及到2个基因，分别是LOC120086723与LOC120071639。氮代谢（Nitrogen metabolism）也涉及到2个基因，分别为LOC120074685与LOC120075868（表4）。

根部中芪类、二芳基庚烷类和姜酚生物合成（Stilbenoid, diarylheptanoid and gingerol biosynthesis）、光合作用天线蛋白（photosynthesis-antenna proteins）、类黄酮生物合成（Flavonoid biosynthesis）与苯丙烷生物合成（Phenylpropanoid biosynthesis）等通路均涉及到1个基因LOC111431282，光合作用天线蛋白（photosynthesis-antenna proteins）涉及到1个基因为LOC111434154（表5）。因此推测LOC111431282基因可能在抵御尖孢镰刀菌侵染发挥重要的作用。

从KEGG富集分析可知，苯丙烷生物合成是冬瓜叶片与南瓜根部共有的富集代谢通路。因此对苯丙烷生物合成进行代谢通路分析（ko00904）。该通路以苯丙氨酸化学分子为底物，经过苯丙氨酸解氨酶（Phenylalanine Ammonia-Lyase, PAL）催化作用变成肉桂酸（Cinnamic acid），而后经过3条不同的代谢通路形成对-羟基苯基木质素（p-hydroxyphenyl lignin, H木质素）、愈创木基木质素（Guaiacyl lignin, G木质素）和紫丁香基木质素（Syringyl lignin, S木质素）3种不同的木质素（图5）。

Figure 5.  Map of phenylpropanoid biosynthesis pathway in leaf（A）and root（B）

分别对嫁接冬瓜叶片284个及根部2个潜在抗病基因进行蛋白序列提取，利用PRGdb网站的DRAGO 3工具对蛋白结构域进行预测。在叶片284个潜在抗病基因中，共预测到16与R基因结构域相关的基因，根据其不同的结构域主要分为7类，分别为CK、CL、KIN、LECRK、N、RLP、TRAM（表6），而根部2个潜在抗病基因没有预测到与R基因相关的结构域。

对叶片所有转录组序列进行TFs分析，共鉴定出2 828个差异表达的TFs，这些TFs属于48个基因家族，主要包括MYB_related、bHLH、MYB、HB-other、ERF、bZIP、MIKC、B3、WRKY、NAC与M_type等（图6−A）。对根部所有转录组序列进行TFs分析，共鉴定出5 426个差异表达的TFs，这些TFs也属于48个基因家族，主要包括MYB_related、bHLH、MYB、ERF、NAC、bZIP、HB-other、ARF、MIKC、B3、WRKY、M_type、FAR1、HD-ZIP、C3H、GRAS与Dof等（图6−B）。

Figure 6.  TFs in response to Fusarium oxysporum infection in leaf（A）and root（B）

叶片与根部共有的基因家族有MYB_related、bHLH、MYB、ERF、NAC、bZIP、HB-other、MIKC、B3、WRKY与M_type。这些TFs调控不同靶基因的表达，可能参与抵抗病原菌入侵反应，在抵御冬瓜枯萎病菌侵染中起着重要作用。

分别选取叶片4个DEGs与根部2个DEGs进行qRT-PCR以验证转录组结果的准确性。结果表明，除冬瓜叶片LOC120086723基因的相对表达量在某些处理间的变化趋势与RNA-seq结果有出入，其他的冬瓜叶片3个DEGs与南瓜根部2个DEGs的相对表达量表现与其对应的RNA-seq结果一致（图7），说明本研究转录组测序的结果可靠。

Figure 7.  qRT-PCR verification of 4 genes from wax gourd in leaves（A-D）and in roots（E-F）

本研究将在嫁接冬瓜中上调，且在自接冬瓜中下调的差异表达基因定义为潜在抗病基因，分别在叶片与根部中筛选出了284个和2个潜在抗病基因。对叶片与根部的潜在抗病基因进行GO与KEGG分析，发现苯丙烷生物合成是嫁接冬瓜叶片与根部共有的代谢通路。苯丙烷代谢的限速酶是PAL，是连接植物新陈代谢中初生代谢和次生代谢的桥梁，在植物木质素、类黄酮、香豆酸酯类等次生代谢物的形成过程以及植物的抗逆调控中发挥重要作用^[19]。因此，推测该代谢通路可能在防御尖孢镰刀菌侵染上扮演着重要角色。

转录因子在植物抗病起着重要的调控作用。逆境胁迫下转录因子可与结构基因上游启动子序列中的顺式作用元件特异性结合，激活或抑制特定抗病基因的表达，从而影响植物对病原体的抗性^[20]。本研究发现，在病原菌侵染后，叶片与根部大量转录因子集中在MYB、b HLH、ERF与NAC等转录因子家族上，MYB、NAC、bHLH、ERF、bZIP等转录因子家族在植物抗病过程中发挥重要作用^{[21 − 23]}。MYB转录因子几乎参与植株生理的各个活动过程，在植株抗逆，尤其是抗病中发挥重要的作用^{[24 − 25]}，NAC转录因子也被发现在参与植物抗病^{[26 − 27]}，且在次生壁的合成上发挥着重要的调控作用^[28]。研究发现转录因子MYB与NAC通过结合PAL1、C4H、4CL1、C3'H1、HCT、CCoAOMT、CCR1、F5H1与CAD6等关键酶基因启动子中的AC元件，激活其表达并调控木质素合成^[29]。本课题组前期研究发现，随着处理时间的延长，病原菌侵染下嫁接冬瓜的木质素含量呈现逐渐上升的趋势^[11]，结合本研究推测，可能是由于大量的转录因子MYB与NAC结合了木质素关键酶基因的启动子中的AC原件，从而促进了木质素的持续合成。木质素的合成沉积有利于细胞壁形成更加坚硬与牢固的物理屏障，进一步减少病原菌的侵染扩散^[11]。

bHLH基因家族是植物中第二大类转录因子家族，广泛参与植物生长发育、信号传导、代谢及逆境胁迫等过程^{[30 − 31]}。bHLH转录因子在抵御各种细菌性与真菌性病害侵染上扮演重要的角色^[32]。南瓜bHLH转录因子CmbHLH87在烟草植株中的异源表达诱导细胞发生了程序性死亡，促进了H₂O₂积累与病程相关蛋白（PRs）的表达，提高了对南瓜白粉病的抗性^[33]。本研究发现，bHLH转录因子数量在叶片与根部都是排名第二，其可能在植株抵御尖孢镰刀菌侵染上发挥重要的作用。

ERF俗称乙烯应答因子（ethylene responsive factor），可以与乙烯诱导顺式作用元件GCC-box结合，在植物抵御植物逆境胁迫中发挥重要的作用^{[34 − 35]}。研究发现，过表达EFR1基因的拟南芥植株对土传尖孢镰刀菌十字花科专化型菌株（Fusarium oxysporum f. sp. Conglutinans）、番茄尖孢镰刀菌番茄专化型菌株（Fusarium oxysporum f. sp. Lycopersici） 引起的枯萎病的抗性增强^[36]。本研究发现，叶片与根部的ERF转录因子数量分别是146与315个，可能也在植株抵御病原菌侵染过程中发挥了重要的防御作用。

综上，MYB、b HLH、ERF与NAC等转录因子在嫁接冬瓜叶片与根部抵御病原菌侵染上发挥重要的作用，但这些转录因子是如何参与尖孢镰刀菌侵染下嫁接冬瓜叶片与根部的调节机制与功能尚未清楚，还需进一步深入研究。

植物抗病基因（R基因）是植物在与病原体的长期斗争中逐渐形成的，它们为植物提供了特定的受体和防御机制^[37]。这些基因共有一些保守的结构域，根据其特定结构域被分成不同的类别，主要分为CNL、TNL, RLK、RLP、LYK、LYP和LECRK^[38]。富含亮氨酸重复序列（LRRs）是大多数R基因蛋白的主要结构域，对识别特异性和抗性至关重要^[39]。例如，多聚半乳糖醛酸酶抑制蛋白（polygalacturonase-inhibiting proteins，PGIPs）包含一个LRR结构域，它限制植物病原体的生长并诱导植物的防御反应^[40]。有些受体蛋白也存在除LRR外的结构域，如LECRK，其含有一个凝集素样基序（lectin-like motifs，LECM）而不是LRR结构域，水稻Pid2蛋白就属于这一类别，其能诱导水稻对稻瘟病菌的抗性^[41]。本研究通过DRAGO 3在叶片的差异表达基因中进一步预测出16个与抗病相关的基因，其中4个具有典型R基因结构域，分别是RLP类：LOC120072345与LOC120085476，LECRK类：LOC120090305与LOC120090267，因此推测嫁接冬瓜可能通过胞外模式识别受体识别病原菌以触发植物的免疫反应，但还需要对预测出的抗病基因进行进一步的分析与功能验证。