0.   引言

【研究意义】百香果（Passiflora edulis）学名西番莲，为西番莲科西番莲属多年生藤本植物。作为新兴树种，百香果已成为我国南方多个省区优先发展的特色果树，因其在栽培生产带动了全产业链发展并创造了极高经济价值，近年来栽培面积持续增加^[1]。镉（Cd）是一种对动/植物极具毒性的重金属元素^[2]。百香果植株对Cd的富集能力较强^[3]，Cd胁迫严重影响果实的品质与商品价值，探讨百香果响应Cd胁迫的分子机理，可为百香果的安全生产提供依据。【前人研究进展 】Cd极易被植物吸收，并在植物的根、茎、叶、果实中大量积累，从而对植物产生严重的毒性，威胁植物的生长发育^[4−7]，尤其是Cd还可通过食物链进入人体，对人体造成潜在损害^[8−9]。许文宝等^[2]在研究百香果对Cd的吸收累积特性时发现，百香果对Cd有较强的吸收能力，枝、叶对土壤镉的富集系数均大于1，且在植物体内的吸收量随着处理浓度的增加而增加，生产中需特别注意Cd的污染风险。【本研究切入点】目前有关百香果重金属污染方面的研究甚少，特别是Cd胁迫方面鲜见报道。【拟解决的关键问题】转录组作为研究生物有机体功能的重要手段，已广泛应用于植物的生长发育、对逆境胁迫适应等分子机制方面的研究。本研究基于转录组测序技术，并结合生物信息学分析挖掘百香果参与Cd胁迫响应的关键基因，从分子水平揭示百香果叶片响应Cd胁迫的机理，以期为研究百香果Cd胁迫的积累规律及其分子机理、百香果安全生产及后续的遗传分子改良提供理论依据。

1.   材料与方法

1.1   试验材料及处理

水培试验于 2023年12月在福建省农业科学院果树研究所培养室进行，供试材料为百香果钦密9号嫁接后60 d苗。环境温度 （24±2） ℃ ，相对湿度75%～90%，光源常规LED灯，昼夜循环，光照14 h/黑暗10 h。将根系基质洗脱后，选取长势良好、均匀一致的植株，一部分植株作为对照在 Hoagland 全营养液[其中KNO₃、Ca(NO₃)₂、(NH₄)H₂PO₄、MgSO₄分别为6、4、2、1 mmol·L⁻¹；H₃BO₃、ZnSO₄、CuSO₄、(NH₄)₆Mo₇O₂₄、FeSO₄、EDTA分别为10、2、2、0.5、0.065、20、20 μmol·L⁻¹，调整pH至 5.8]中培养（对照，处理编号为T0），一部分植株在 Hoagland 全营养液中加入50 μmol·L⁻¹ 的 CdCl₂ 培养，处理24、48、72 h后分别采集幼嫩叶片（处理编号分别为T24、T48和T72）。为保证重金属Cd²⁺浓度，每 24 h换营养液1次。每个处理设置3个生物学重复，每个重复5株苗。采集的样本，液氮速冻后放置于−80 ℃冰箱备用。

1.2   转录组测序

对于转录组测序，样本组织使用RNeasy Plant Mini Kit （Qiagen，上海玉博生物科技有限公司）提取总RNA反转录成cDNA。利用Illumina HiSeq 3000平台获得原始序列读取，并使用默认设置的Bowtie2 version 2.1.0映射到Ensemble MSU6转录集。严格过滤原始测序序列（Raw reads），去除杂质数据，获得干净序列（Clean reads）。对干净序列进行组装和拼接，再将组装的转录本进行聚类和去冗余处理后用于后续分析。

1.3   数据分析

FPKM值指每百万碱基对测序所映射的每千碱基转录本序列片段的预期数量，edgeR程序用于鉴定差异表达基因。错误发现率FDR （False discovery rate）＜ 0.05且变化大于1倍的基因被认为是差异表达基因。基因本体（Gene ontology, GO）富集分析和KEGG代谢通路进行富集分析使用R语言Bioconductor包，热图由软件TBtools绘制。

2.   结果与分析

2.1   转录组测序数据分析

从幼苗叶片中分别提取总RNA，共构建了12个cDNA文库（T0、T24、T48和T72）。转录组测序结果显示，原始数据处理后共得到446937924 bp高质量碱基，各个样品的原始数据经去除接头和低质量序列后，比对基因组的平均比对率为94.34% ；检测到的有效基因数为18181个。各处理组重复之间的相关性分析结果显示（图1），除T48-2与处理组其他生物学重复的相关系数小于0.8（分别为0.55、0.56）外，其余各组样本间相关系数均大于0.8。因此，后续分析数据时将T48-2剔除，以确保测序结果的真实可靠性。

图  1  样品相关性分析

Figure  1.  Correlations on samples

下载: 全尺寸图片 幻灯片

与对照T0相比，本研究T24、T48和T72的差异表达基因数分别为3465、1262 和 2039个，其中上调基因数为1819、842和867个，下调基因数为1646、420和 1172 个（图2）。可以看出Cd处理后24 h，差异表达基因数目远大于48、72 h处理的试验组。因此，百香果叶片响应Cd胁迫差异基因主要是在胁迫后24 h内发生。

图  2  差异基因分析

Figure  2.  DEGs analysis

下载: 全尺寸图片 幻灯片

2.2   差异表达基因的GO功能富集分析

为明确 Cd 胁迫下百香果差异表达基因的生物学功能和属性，将 Cd 处理 24、48、72 h 与对照之间的差异表达基因进行 GO富集分析。GO注释结果显示，差异表达基因（Differentially expressed genes, DEGs）主要被注释到了 3 个功能组（图3）：第一个是生物学过程，对于 T24 和 T48 DEGs 富集最多的是光合作用（Photosynthesis）和光响应（Photosynthesis，light reaction），而对于 T72 DEGs 富集更多在于碳水化合物代谢过程（Carbohydrate metabolic processes）；第二个是细胞组分方面，对于 T24 和 T48 DEGs 富集最多的是细胞相关（Extracellular region, apoplast）及类囊体相关（Thylakoid, thylakoid membrane），而对于 T72 DEGs 富集在胞外及细胞壁等 5 个方面的数量均无太大变化；第三个是分子功能方面，对于 T24 和 T48 DEGs 富集最多是的氧化还原酶活性（Oxidoreductase activity），而对于 T72 DEGs 富集最多是催化活性（Cataltyic activity）。

图  3  差异表达基因的 GO 功能注释

A：T24 vs T0差异表达基因的GO功能注释散点图；B：T48 vs T0差异表达基因的GO功能注释散点图；C：T72 vs T0差异表达基因的GO功能注释散点图。富集因子代表富集到 GO功能注释 DEG 数量与 DEGs 总数的比值。

Figure  3.  Go functional annotation of DEGs

A: GO function annotation scatter plot of T24 vs. T0 DEGs enrichment; B: GO function annotation scatter plot of T48 vs. T0 DEGs enrichment；C: GO function annotation scatter plot of T72 vs. T0 DEGs enrichment. Rich factor indicated ratio of number of enriched DEGs in GO function annotations to total number of DEGs.

下载: 全尺寸图片 幻灯片

2.3   差异表达基因的KEGG代谢通路分析

为分析Cd胁迫下百香果差异表达基因的代谢调控通路，将Cd处理24、48、72 h分别与对照之间的差异表达基因进行KEGG代谢通路分析（图4）。结果表明，在所有处理中，DEGs主要富集在EC相关酶的数量（Enzymes with EC numbers）、转运蛋白（Transporters）以及激素信号转导（Plant hormone signal transduction）方面；对于T48和T72处理 DEGs在MAPK信号途径有富集，而前期T24处理未富集。表明百香果幼苗在响应Cd胁迫时，可能是通过激活激素来响应Cd胁迫，且在后期激活MAPK途径。

图  4  KEGG 富集散点图

A：T24 vs T0差异表达基因的KEGG富集散点图；B：T48 vs T0差异表达基因的KEGG富集散点图； C：T72 vs T0差异表达基因的KEGG富集散点图。富集因子表示富集到 KEGG 通路的 DEG 数量与 DEGs 总数的之比。

Figure  4.  KEGG enrichment scatter plot

A: KEGG pathway scatter plot of T24 vs. T0 DEGs enrichment; B: KEGG pathway scatter plot of T48 vs. T0 DEGs enrichment；C: KEGG pathway scatter plot of T72 vs. T0 DEGs enrichment. Rich factor indicated ratio of number of enriched DEGs in KEGG pathways to total number of DEGs.

下载: 全尺寸图片 幻灯片

2.4   差异表达基因分析及候选基因的筛选

通过不同胁迫时间差异表达基因之间的交集，可获得其连续性差异表达基因，百香果叶片响应Cd胁迫的差异基因主要富集在代谢途径和生物合成相关途径（图3、图4）。对T0、T24、T48和T72处理的差异表达基因进行维恩分析（图5A），结果表明， Cd处理下T24、T48、T72与T0之间3个时期连续性表达的DEGs有416个，而不连续性表达的DEGs有4882个。

图  5  不同时间点 Cd 处理下的连续性差异表达基因及部分转录因子的差异表达基因

A：不同时间点Cd处理下的差异表达基因韦恩图；B和C ：部分差异表达基因的表达热图。

Figure  5.  Continuous and partial DEGs under Cd stress at different times

A: Venn diagram of DEGs treated with Cd at different time point; B, C: heat map of partial DEGs.

下载: 全尺寸图片 幻灯片

本研究发现63个转录因子包括14个WRKY（WRKY）、22个NAC（NAM、ATAF1/2、CUC1/2）、13个MYB（v-myb avian myeloblastosis viral oncogene homolog）、7个bHLH（basic helix-loop-helix）、5个ERF（Ethylene-responsive factor）及2个bzip[Basic （region）leucine zipper]等在Cd处理下差异表达；同时还发现转运Cd²⁺的Nramp（Natural Resistance-Associated Macrophage Protein）和HMA（Heavy metal transporting ATPase）家族中各一个成员被诱导表达。在这些转录因子中，NAC（Pe3g02665, NAC domain-containing protein 86）和EFR（Pe4g00254, Ethylene-responsive transcription factor ERF039-like）在Cd胁迫下被诱导表达超过100倍；而WRKY中的两个基因Pe9g01955（WRKY transcription factor 35 family protein）和Pe4g02129（WRKY family transcription factor）被抑制表达倍数最高。

2.5   差异表达基因的验证

为了进一步验证转录组数据可靠性，随机挑选WRKY、NAC、bHLH、MYB以及HAM家族各1个共5个差异表达基因，以前人发表的S60作为内参^[10]，用Primer 5软件设计特异性引物（表1），利用实时荧光定量PCR（RT-qPCR）检测基因表达水平的变化，PCR反应体系为：2× Master Mix 5 μL，上游引物（F）0.2 μL，下游引物（R）0.2 μL，cDNA 1 μL，ddH2O补充至10 μL。反应程序为：50 ℃ 2 min；95 ℃ 5 min；95 ℃15 s，60 ℃ 1 min，40个循环；95 ℃ 15 s；60 ℃ 15 s；95 ℃ 15 s。用2−ΔΔCT计算基因相对表达量。结果如图6所示，百香果幼苗Cd胁迫处理组的基因表达量变化趋势与转录组测序结果基本一致，表明转录组测序结果有较高的可靠性。

表  1  实时荧光定量PCR的引物序列

Table  1.  Sequences of qRT-PCR primers applied

基因代码 Gene ID	正向序列 Forward sequences (5′–3′)	反向序列 Reverse sequences (5′–3′)
60S	AGGTGGGTAACAGGATTATC	TGGCTGTCTTTTGGTGCTG
WRKY31	CAGGTTCAACGTTGTGCAGA	AGCAAGGAAGGATGGCTCTT
NAC100	AAAGCCTCTGATCAACCCCA	TGCAGCTTCTCCATGACAGA
MYB306	TGCAGGCTTAGATGGACCAA	GTCTTGTGTCAGAGGGTCCA
bHLH87	GCCGAACGTTCATCCAAAGA	TTCTGGGTCAGCTGGTTCTT
HMA3	AGTGGTAGGAACAATCGCCA	CCGGTTTCTGCTATGACTGC

下载: 导出CSV 

| 显示表格

图  6  部分差异表达基因的 qRT-PCR 验证

*：P ＜ 0.05 水平差异显著； **： P ＜ 0.01 水平差异极显著。

Figure  6.  Verification of DEGs by qRT-PCR

*: significant difference at P＜0.05; **: extremely significant difference at P＜0.01.

下载: 全尺寸图片 幻灯片

3.   讨论与结论

Cd是植物生长发育的非必需元素，但它能通过与膜蛋白竞争性结合进入植（动）物细胞，从而破坏正常的生理过程，对细胞造成损害^[11]。近年来已有学者通过转录组测序技术分析多种作物在Cd胁迫处理下的基因表达模式，并揭示了菜用甘薯、艾纳香及高粱等作物对Cd胁迫的响应机制 ^[12−13]。本文首次应用转录组测序技术分析了百香果幼苗叶片响应中度Cd胁迫的基因表达模式，最终获得了63个与Cd胁迫相关的转录因子及蛋白，研究结果丰富了百香果对Cd胁迫响应的调控网络，为阐明百香果响应Cd胁迫的分子机制提供依据。

前人研究表明，在植物体内，Cd可影响叶绿素合成，降低叶片光合速率，进而影响植株干物质积累^[14]。本研究发现，在百香果幼苗响应中度Cd胁迫的过程中，差异表达基因主要富集于光合作用、碳水化合物代谢的相关通路中（图4），表明Cd胁迫对百香果可对叶片的光合作用造成干扰，通过影响光合产物生成而影响植物的正常生长。结合KEGG代谢通路分析还发现大量差异基因富集在相关酶和各种胞膜上，当大量Cd²⁺进入细胞时，可能导致细胞内膜系统发生不可逆损伤，引起细胞器的功能紊乱进而影响植物的正常代谢^[15]。

植物激素不仅在植物生长发育调节中发挥着重要作用，还深度参与植物对多种逆境胁迫的响应^[16]。如番茄响应Cd胁迫时，茉莉酸可经机体渗透调节和抗氧化酶系统功能提升增强植株抗性^[17]；乙烯和油菜素内酯能够通过抗氧化剂表达水平的提升来降低Cd的毒害作用^[18−19]；在生菜中，赤霉素则可通过增强Cd²⁺转运基因的表达抑制Cd²⁺吸收从而提高植株对Cd毒害的抗性^[20]。百香果叶片响应Cd胁迫过程的差异表达基因中，前期与植物激素信号转导相关的基因不在少数，因此百香果响应Cd胁迫处理过程中也可能是通过激活某些激素而后激活MAPK途径完成的。

在植物响应非生物胁迫中，转录因子也扮演了极为重要的角色。如AtHMA3是影响拟南芥叶片Cd含量的关键基因^[21]；Zm HMA3a是控制玉米种子Cd积累的重要基因^[22]；TdHMA3不同等位基因的活性与硬粒小麦（Triticum turgidum var. durum）籽粒Cd的积累量有一定的关系 ^[23]；而OsNramp5中第337位点的谷氨酰胺突变会影响水稻对Cd的吸收 ^[24]。本研究发现百香果幼苗Cd胁迫处理下有大量转录因子差异表达，同时发现转运Cd²⁺的Nramp家族和HMA家族成员也被诱导表达，这也与前人的研究结果一致。

综上，本研究通过对百香果幼苗中度Cd胁迫的转录组进行分析，发现在不同时间处理下共有63个特异性表达的相关转录因子及蛋白，差异基因主要富集在细胞结构体、催化活性和转录调节活性三个功能组；对GO分析其主要富集在光合作用和碳水化合物代谢通路，表明Cd可能通过干扰百香果的光合作用，影响光合产物生成，从而影响百香果的正常生长。本试验仅在百香果幼苗期进行Cd处理，百香果植株生育后期对Cd的胁迫响应以及各种转录因子对Cd胁迫的响应机制有待于进一步研究。