PAHs污染胁迫下翅碱蓬基因表达响应机制研究

多环芳烃(Polycyclic Armomatic Hydrocarbons，PAHs)是一种三致性环境污染物，倾向于在环境中积累，一旦被生物吸收，会对细胞造成损害，抑制生物正常生长发育[1-2]。

由于沿海城镇高速发展、企业生产工业化和石油开采等原因，PAHs对海岸带造成的污染日益严重[3]。盘锦红海滩湿地位于辽河入海口，对阻拦并净化辽河尾端污染物、保持辽东湾环海地区及近海健康生态系统平衡具有重要作用。近年来，由于环境污染加剧，红海滩翅碱蓬(Suaeda heteroptera)生态系统出现生态退化现象[4]。

碱蓬是一种耐盐碱的一年生滩涂先锋植物，研究表明其对污染物等具有良好的去除能力[5]。辽河油田作为典型化石能源开采区，其周边环境介质中普遍存在以菲和苯并[a]芘为代表的多环芳烃复合污染物。二者的赋存具有明确的石油源特征，主要源于油田勘探、开采、运输及加工过程中的原油及其副产物的泄漏与排放。大量实证研究证实，在油田区的土壤、底泥及水体中，菲与苯并[a]芘呈现出显著富集的态势，其空间分布与工业生产活动的强度呈现高度相关性[6-7]。PAHs作为滨海湿地的重要污染物，翅碱蓬响应该胁迫的分子机制目前尚不清楚。

组学测序分析是研究植物非生物胁迫抗逆机制的有效手段。转录组测序可以得到植物在胁迫条件下富集的转录序列，对于挖掘植物抗逆机制的关键基因有重要意义。Wang等[8]通过转录组分析进一步揭示生长素处理协同增强了植物激素信号转导途径与ABC转运蛋白通路，并特异性调控SAUR生长素早期响应基因及ABC转运蛋白家族基因的表达，从基因表达层面揭示了生长素调控PAHs吸收的分子网络。陆程张[9]研究Na2CO3胁迫下玉米苗根部转录组基因表达时发现，植物MAPK 信号通路、植物激素信号转导通路的基因被显著富集。

Zeng等[10]从分子水平揭示了紫花苜蓿响应盐胁迫的主要因子。通过对紫花苜蓿的转录组进行分析，筛选出了紫花苜蓿中与盐胁迫相关的差异表达基因[11]。植物在胁迫条件下可以通过调节自身的基因表达、信号传导、能量代谢及蛋白质合成等过程，抵御胁危害并提高耐胁迫能力[12-13]。

本研究中通过高通量测序技术对PAHs胁迫的翅碱蓬植株进行转录组测序，分析不同PAHs浓度下翅碱蓬基因表达的差异，从转录水平揭示翅碱蓬适应PAHs环境的分子机制，以期为培育翅碱蓬等盐生植物消减PAHs、改良盐碱地提供参考资料，丰富耐盐植物响应持久性有机污染物(Persistent Organic Pollutants，POPs)胁迫机制数据库。

1 材料与方法

1.1 材料

1.1.1 仪器与药品主要使用的仪器为电子天平、超声波清洗器、高效液相色谱仪、加速溶剂提取仪和色谱柱等。主要试剂为丙酮(分析纯、色谱纯)、二氯甲烷(分析纯、色谱纯)、正己烷(分析纯、色谱纯)和PAHs混合标准品等。

1.1.2 供试植物试验所用翅碱蓬为前期土培种植的幼苗。翅碱蓬种子购于辽宁省盘锦市，土壤为市面售卖的花土。

1.2 方法

1.2.1 前期翅碱蓬育苗准备将翅碱蓬种子种植在花盆中，每盆播种翅碱蓬种子1.0 g(约200粒)。将种植好的翅碱蓬置于室外有光照且通风处，每2天用自来水浇灌一次，保持土壤湿润。60 d后选取长势相近的植株幼苗，株高约为10 cm，用清水清洗植物的根部及叶片，防止其对试验造成干扰。然后将植物置于营养液中培养3 d，使其适应新环境，之后进行试验。

1.2.2 PAHs胁迫方法本研究设置的PAHs胁迫浓度(1、10、100 μg/L)基于盘锦红海滩湿地及周边辽河油田区域实际监测数据。据文献报道，辽河口湿地0～20 cm土壤16种PAHs总量变化范围为268.7～2 853.8 ng/g，平均值为1 241.9 ng/g，且表层海水中PAHs均属于轻度-中度污染，污染浓度范围为50～1 000 ng/L[14-15]。为模拟从背景值到污染热点区的浓度梯度，并覆盖植物可能遭遇的实际污染范围，本试验设置了低(1 μg/L)、中(10 μg/L)、高(100 μg/L)3个胁迫浓度，以期系统揭示翅碱蓬在不同污染水平下的转录响应模式。

用去离子水将生长时间相同、生长条件一致的翅碱蓬(约10 cm)清洗干净，放入盛有霍格兰(Hoagland)营养液的烧杯中，采用水培的试验方法并在光照培养箱中稳定2 d，在第3天时，使用不同浓度的菲和苯并[a]芘混合溶液对翅碱蓬进行处理，每个浓度约放入25根翅碱蓬，设置3组平行，每3天更换一次营养液，保持PAHs浓度不变，试验周期为12 d。于第12天每组各取3株翅碱蓬(约1.2 g)，取样时使用去离子水将翅碱蓬苗冲净，吸干翅碱蓬样品表面的液体，置于干净的离心管中，用记号笔标记好后，立即置于液氮速冻。样品经液氮冷冻后立即送华大基因用于翅碱蓬转录组分析。

1.2.3 RNA提取、文库构建和测序委托华大基因对受不同浓度PAHs胁迫后的翅碱蓬进行转录组测序工作，转录组测序方法均来自该公司。为方便后续分析，将对照组和1、10、100 μg/L的PAHs处理组记为A、B、C、D组，由于每组设置3个平行，故将12个样品分别记为A1、A2、A3、B1、B2、B3、C1、C2、C3、D1、D2、D3；将A、B、C、D组两两对比分为6个对比组，记为A-vs-B、A-vs-C、A-vs-D、B-vs-C、B-vs-D和C-vs-D。

采用 Trizol 试剂(Invitrogen 公司)提取翅碱蓬总RNA。使用Nanodrop检测RNA的纯度(OD260 nm/OD280 nm)；利用 Qubit 对 RNA 浓度进行精确定量；参照Agilent 2100电泳法精确检测RNA的完整性。

RNA质检合格后，对mRNA进行富集。随后用打断buffer把获得的RNA片段化。用六碱基随机引物进行反转录，再合成cDNA二链形成双链DNA。将双链 cDNA进行纯化、末端修复、加Poly(A)尾并连接测序接头。连接产物通过特异的引物进行PCR扩增。PCR产物经热变性解链为单链DNA，用引物将单链DNA环化得到单链环状DNA文库[16-17]。之后通过DNBSEQ平台(华大智造)上机测序。

1.2.4 转录组的组装、注释和差异表达基因(DEGs)分析测序的原始数据包含低质量、接头污染及未知碱基N含量大于5%的reads，数据分析之前需要去除这些reads以保证结果的可靠性。采用过滤软件SOAPnuke进行过滤。对过滤后得到的Clean Data的质量和数据量进行统计，包括质量值大于30的碱基数比例(Q30)、数据量和碱基含量等。

采用 Trinity 软件对clean reads进行组装，然后使用Tgicl对转录本进行聚类去冗余得到Unigene。对组装的转录本进行质量评估，通过与保守基因进行比较，评价转录组组装的完整性。为获得全面的基因信息，将组装得到的Unigene进行7大功能数据库注释(KEGG 、GO 、NR、NT、SwissProt、Pfam和KOG)。然后使用RSEM计算各个样品的基因表达水平[18]。采用DEseq2方法对差异表达基因进行分析[19]。对DEGs(Q≤0.05 且|log2 ratio|≥1)进行通路富集分析。

1.2.5 实时荧光定量PCR(qRT-PCR)分析为验证RNA-seq结果的可靠性，本研究以β-actin(CL7574.Contig3_All)为内参基因，从转录组6个对比组中分别挑选出11、7、9、8、8、13条差异基因。根据转录组文库中的Unigene 序列，用GenScript设计特异性引物(表1～3)。RNA提取和质量检测方法与“第1.2.3节”相同。RNA反转录采用TIANGEN公司FastKing gDNA Dispelling RT SuperMix 试剂盒，合成cDNA。反应结束后，cDNA模板保存于-20 ℃。使用实时荧光定量PCR仪(Eppendorf，德国)和TB Green Premix Ex TaqTM Ⅱ试剂进行检测。每个样品设置3组重复。试验结果采用2-ΔΔCt法[20]计算基因的相对表达量。

2 结果与分析

2.1 测序结果及组装

将A、B、C、D组进行测序，每组3个样品共12个样品，每个样品的原始读取数据(raw reads)数量均为4 382万条，将一些质量低、受接头污染的reads过滤掉后得到的高质量clean bases数量为4 208～4 242万条，12组clean bases数量为76.11 Gb，且Q20均达到97%以上，Q30均达到93%以上，Clean Reads Ratio为95.95%～96.79%。确认转录组组装的完整性可通过单拷贝直系同源数据库BUSCO完成，12个样品及All-Unigene在BUSCO的组装评价可知，All-Unigene有131个能完全比对上BUSCO数据库的单拷贝序列，162个完全比对上BUSCO数据库的多拷贝序列；A组中A1和A3的Complete(C) and Single-copy(S)最多，但A3的Complete(C) and duplicated(D)略大于A1；B组Complete(C) and Single-copy(S)数量排序为B3

2.2 基因功能注释

组装后共有80 222条Unigenes在7大公共数据库(KEGG 、GO 、NR、NT、SwissProt、Pfam和KOG)中完成基因功能注释。以上7大数据库中注释的Unigenes数量从37 083条(GO占46.23%)～61 825条(Nr占77.07%)，其中能成功在7大公共数据库完成注释的Unigenes，占总Unigene数比的25.91%(共20 786条)；共64 526条Unigenes能成功在7大公共数据库中至少一个数据库获得注释，占总Unigene数比的80.43%。图1为选取组装后的Unigenes在5大数据库(NR、NT、GO、Pfam和KOG)注释结果韦恩图，共有23 930条Unigenes在5大数据库(NR、NT、GO、Pfam和KOG)均有注释，占总Unigene数比的29.82%。

2.2.1 GO分类基因在细胞、分子和组织水平上的功能通常利用GO来描述，将所有比对上NR数据库的Unigene结果使用GO 数据库进行了注释，并针对GO的3个分类— 生物学过程(biological process)、细胞组分(cellular component)及分子功能(molecular function)，对基因产物的功能进行描述。从图2可见，三大分类共获得45个子节点(term)，生物学过程分类的子节点最多，数量为26，其中，细胞进程(cellular progress)的基因数量最多，为18 521条，其次为代谢过程(metabolic process)有15 173条基因；再次是分子功能分类有16个子节点，在该分类中基因数量超过一万的为binding和catalytic activity；细胞组分分类的子节点最少，数量为3，其中cellular anatomical entity所占基因数目最多。

2.2.2 KOG分类通过KOG数据库对Unigene序列进行注释以获得基因同源物的分类信息(真核生物)，将所得结果绘制成图(图3)。在KOG功能注释中有46 736条Unigenes得到功能注释，根据功能进一步细分成25组，将各已知功能组别按照注释基因数量从多到少排序，前4名依次为仅一般功能预测(general function prediction only)，信号转导机制(signal transduction mechanisms)，翻译后修饰(posttranslational modification)、蛋白质更新(protein turnover)、分子伴侣(chaperones)，以及转录(transcription)，基因个数分别为10 368、5 409、3 993和2 616条。

2.2.3 KEGG分类从图4可见，统计在KEGG代谢通路获得注释Unigenes并将其分类，获得注释上KEGG数据库level 1层级共5组，进一步细分注释到level 2层级有19类。5组level 1层级当中有2组organismal systems、cellular processes均只有1类level 2层级，注释分别为环境适应(environmental adaptation)(1 751条Unigene)、运输与分解代谢(transport and catabolism)(1 970条Unigene)；environmental information processing中信号转导(signal transduction)的注释最多(2 448条Unigene)；genetic information processing中蛋白质合成(translation)和折叠、分选与降解(folding，sorting and degradation)数量较多(3 462条Unigene、2 912条Unigene)；代谢(metabolism)组有11类level 2层级，其中全局与概览图(global and overview maps)数量最多(10 443条Unigene)。

2.3 基因差异表达结果

翅碱蓬在不同PAHs浓度下各组差异基因数量统计如图5所示，A-vs-C、A-vs-B对比组产生的差异基因数量较多，分别为641、638条，其中上调差异基因分别为202、408条，下调差异基因分别为439、230条；A-vs-D对比组产生差异性基因较少(378条)，其下调差异基因数量要大于上调差异基因，由此可推断，随着PAHs浓度的增加，翅碱蓬响应PAHs的差异基因数先增加后减少，下调差异基因数大于上调差异基因数。

翅碱蓬差异基因数量的“先增后减”模式，本质上反映了植物从主动防御适应转向被动毒性抑制的过程。这一现象提示在利用植物修复PAHs污染时，存在一个“最佳胁迫浓度窗口”，过低则诱导不足，过高则抑制植物活性，数据表明，PAHs胁迫存在一个毒性阈值[21-22]。低于该阈值，植物启动主动基因调控网络以应对胁迫；高于该阈值，则进入被动损伤与抑制状态。简言之，差异基因数量的“先增后减”是植物从积极适应转向耐受极限的分子体现，反映了胁迫强度与生物响应能力之间的非线性关系。

图6根据A-vs-D、A-vs-C和A-vs-B 3个对比组差异基因绘制的差异基因韦恩图，3个对比组的差异重叠基因有6个，在A-vs-D和A-vs-C对比组中有45个差异重叠基因；A-vs-D、A-vs-B对比组中差异重叠基因有37个；A-vs-B、A-vs-C对比组的差异重叠基因数量(75个)大于A-vs-D和A-vs-C对比组、A-vs-D和A-vs-B对比组。

2.4 差异基因分析

2.4.1 差异基因的GO富集分析将PAHs浓度胁迫的翅碱蓬基因表达同对照组基因表达对比，共获得378条差异基因，针对差异基因检测结果进行GO功能分类(图7(a))，并注释到GO的3个分支(BP、CC和MF)中，其中BP和MF分支注释差异基因较多，分别占差异基因总数的55.82%和59.26%；差异基因注释到分支CC的数量略少，约占差异基因总数的46.30%，这表明大多数差异基因都涉及某些生物学过程和相关分子功能。BP分支有14个GO term，其中细胞过程差异表达基因数量最多，为86条；CC分支有3个GO term，细胞解剖实体(cellular anatomical entity)的差异表达基因为106条，数量最多；MF分支有6个GO term，差异基因参与比例较高的为binding和catalytic activity，分别为101、102条。

把A-vs-D对比组注释到GO的差异基因，根据功能分类后进行富集计算，Q value≤0.05的功能视为显著富集。按Q value值从小到大对GO Term进行排序(Q value≤0.05，且值越接近0富集越显著)，取前20个绘制GO富集结果图，即图7(b)。从富集结果发现，富集在CC分支的内质网膜(endoplasmic reticulum membrane)与富集在MF分支的磷酸吡哆醛结合(pyridoxal phosphate binding)部分的差异基因富集效果显著且数量较多，分别为6、5个。磷酸吡哆醛是维生素B6的主要辅酶形式，不仅可催化氨基酸代谢[23]，也对磷脂的合成有影响，在植物生长过程中，磷酸吡哆醛结合可促进植物根的发育，通过调节离子泵增强根系对养分离子的吸收，提升植物抗逆性[24]，以便植物适应生存环境；内质网是细胞质膜系统，用于合成膜蛋白和分泌蛋白，张敏等[25]指出，后沟树在盐胁迫后内质网膜各蛋白组分表达量会受到影响，且内质网在根与叶的蛋白质组成不同，有研究表明，内质网膜可以调节细胞内Ca2+的信号传导，与植物抗逆性有密切关系。由此推测，PAHs胁迫激活了翅碱蓬磷酸吡哆醛结合和内质网膜功能，以抵御不良的生存环境。

从图8(a)可见，对A-vs-C对比组产生的差异基因进行分类，注释到GO 3个分支(MF、BP和CC)的差异基因分别有408、385、295条，分别占差异基因总数的63.65%、60.06%和46.02%。MF分支有8个GO term，其中binding和catalytic activity发生差异基因较多，分别有172、182条；BP分支有17个GO term，其中cellular process差异表达基因有137条，数量最多；CC分支有3个GO term，其中cellular anatomical entity差异基因数量最多，为166个。

为分析A-vs-C对比组的差异基因功能分类，将GO Term以Q value值从小到大排序选前20个绘制的差异基因GO富集气泡图结果(图8(b))，显著富集(Q value值≤0.05)在CC分支的细胞解剖实体分类中的chloroplast的差异基因有19个，数量最多。发生光合作用的场所以叶绿体为主，叶绿体对植物的生长发育有至关重要的影响，说明随着PAHs浓度增加，光合作用进一步受到影响。

A-vs-B对比组产生638条差异基因，其中有11个GO term属于分支MF，3个GO term属于分支CC，19个GO term属于分支BP，参与到3个分支的差异性基因分别有399、263、355条，分别占差异基因总数62.54%、41.22%和55.64%。差异基因参与细胞成分(CC)功能中cellular anatomical entity比例最高，为158条；MF分支中差异基因所占比例最高的是binding和catalytic activity，分别为182、181条；BP分支cellular process差异表达基因的数量最多，为124条(图9(a))。

A-vs-B对比组差异基因GO富集结果如图9(b)所示，GO富集气泡结果显示，显著富集(Q value≤0.05)MF分支催化活性功能的核酮糖-二磷酸羧化酶(ribulose-bisphosphate carboxylase activity)与BP分支代谢功能的光呼吸(photorespiration)部分差异基因最多，数量均有4条。核酮糖-二磷酸羧化酶是与光合作用中CO2的固定有关的酶，由此可推断核酮糖-二磷酸羧化酶的显著富集与光合作用有关，目前已有研究指出，核酮糖-二磷酸羧化酶可维持植物光合作用，植物叶片固定CO2的能力与该酶有关[26]；光合作用里消耗能量副反应之一为光呼吸，光呼吸强度与光合作用和净光合速率成反比[27]，由此推测在低PAHs胁迫下，翅碱蓬光合作用受到影响。

2.4.2 差异基因的KEGG富集分析基于KEGG注释结果和7类KEGG代谢通路，将A-vs-B、A-vs-C和A-vs-D 3个对比组的差异基因分别进行通路分类，富集分析利用phyper函数计算，FDR用于p value校正，并分别对各对比组功能按照Q value值，从小到大进行排序，取各对比组取top20用于绘制差异基因KEGG代谢通路富集气泡图(图10)。

A-vs-B对比组显著富集的通路是乙醛酸和二羧酸代谢(glyoxylate and dicarboxylate metabolism)通路，注释到该通路的差异基因数量为17条(9条上调)；A-vs-C对比组代谢通路富集程度较为显著的有果糖和甘露糖代谢(fructose and mannosemetabolism)，光合生物中的碳固定(carbon fixation in photosynthetic organisms)、核黄素代谢(riboflavin metabolism)、碳代谢(carbon metabolism glyoxylate)和乙醛酸和二羧酸代谢(glyoxylate and dicarboxylate metabolism)，差异基因注释到以上通路的数量分别为8、12、7、26、13条(上调基因数分别是0、2、6、7、6条)；A-vs-D对比组有2条代谢通路富集显著，为磷酸戊糖通路、碳代谢通路，注释到这两个通路的差异基因为9、20条(上调基因分别是5、9条)。

2.5 实时荧光定量PCR验证

为验证本次转录组数据的可靠性，本试验根据转录组测序结果对部分基因进行验证，在A-vs-B、A-vs-C、A-vs-D、B-vs-C、B-vs-D和C-vs-D对比组分别选取了11、7、9、8、8、13条差异基因进行qRT-PCR检测。各组qRT-PCR检测结果与转录组数据表达趋势一致，说明本试验转录组数据可靠。

3 讨论

本研究中，差异基因主要KEGG富集通路见表1。通过富集程度对比发现翅碱蓬突变体主要通过乙醛酸和二羧酸代谢、碳代谢和磷酸戊糖通路来响应PAHs胁迫。

3.1 乙醛酸和二羧酸代谢通路

乙醛酸和二羧酸代谢通路是植物在遇到外界逆境干扰时，平衡体内局部紊乱，输送能量以提高抗逆性的代谢通路[28]。根据KEGG注释结果，A-vs-D对比组中有8条差异基因注释到乙醛酸和二羧酸代谢所在通路。从图4～图13可见，红框为上调基因，绿框为下调基因，其中6条为下调基因，包括编码甘氨酸脱羧酶复合物(GDC)3条、丝氨酸羟甲基转移酶(E2.1.2.1)1条。

GDC是一种由4个亚基组成的多酶复合物，在植物叶片线粒体基质中浓度极高。该复合物催化光呼吸过程中从过氧化物酶体溢出的甘氨酸分子的分解[29]。陈帅等[30]研究了以甘氨酸脱羧酶P-蛋白亚基参与嵌合结构在水稻中的特异性表达情况，发现β-葡萄糖苷酸酶活性(GUS)不仅受水稻发育程度影响，在不同组织中活性也有差异。另有研究表明，上述嵌合结构在转基因烟草中亦可特异性表达，而在不同组织GUS活性与水稻不同，这可能与各类植物组织细胞起源与分化不同有关[31]。E2.1.2.1不仅参与到多种催化反应，也是丝氨酸与甘氨酸相互转化的关键酶，且植物细胞中该酶的同工酶种类要多于动物细胞，常见于线粒体或细胞质中[32]。E2.1.2.1可在嘌呤和蛋白质合成过程中源源不断的提供甘氨酸，亦为Cl库提供N5，为合成DNA等多种产物需用到的专一性Cl化合物打下基础[33]。有研究指出，在植物正常生长发育过程中，丝氨酸经一系列的催化反应最终可使甘氨酸完成氧化，并提供大量的Cl化合物，这与由E2.1.2.1与GDC协同催化密切相关，植物的生长与光合作用会受到GDC含量与E2.1.2.1活性的影响[34]。随着PAHs浓度升高，乙醛酸和二羧酸代谢通路富集下调差异基因比重增加，由此可推断在PAHs胁迫下该途径过程受到了抑制，导致编码GDC与E2.1.2.1的相关基因下调，抑制翅碱蓬的生长与光合作用。

3.2 碳代谢通路

植物的基础代谢中最重要的就是碳代谢，王诗雅[35]指出，碳代谢对植物生长发育和产量的影响巨大，与植物体内核酸与蛋白质的合成有也有密切联系。整体来看，随着PAHs浓度增加，碳代谢通路富集的差异基因数呈先增加后减少趋势，上调差异基因数占比先降低后增加，由此可推测低浓度PAHs胁迫对翅碱蓬碳代谢起促进作用，高PAHs浓度抑制翅碱蓬碳代谢。在A-vs-D对比组碳代谢通路中检测到20条差异基因，其中编码Glucose-6P的转化过程与核糖5-磷酸(Ribose-5P)与核酮糖5-磷酸(Ribulose-5P)相互转化过程均有基因上调。

Glucose-6P不仅是胞质内葡萄糖的存在形式，其转化也是碳代谢的关键步骤，进入植物体内的葡萄糖会由磷酸葡糖激酶生成Glucose-6P[36]，有研究指出，Glucose-6P可转化成参与糖酵解的关键化合物果糖-6磷酸，继而参与能量代谢过程；也可转化为磷酸戊糖途径所需要的NADPH和五碳糖，进而参与糖合成过程[37]。Ribose-5P与Ribulose-5P相互转化是在D-ribose-5-phosphate aldose-ketose-isomerase(E5.3.1.6)参与下进行的，且该反应在碳代谢等基础代谢中有重要地位，并在生物生长初期阶段格外显著[38]，这与尚瑞沙等[39]研究结果一致。已有研究表明[40]，Ribose-5P是合成核苷酸的原料之一，其转化对植物生长发育具有重要作用，缺乏Ribose-5P与Ribulose-5P相互转化，将导致戊糖磷酸代谢途径不能正常进行。本研究中编码Glucose-6P的转化与Ribose-5P与Ribulose-5P相互转化过程有基因上调，由此推断翅碱蓬受到PAHs胁迫后，通过调节上述两种转化过程来促进碳代谢途径，以抵御PAHs的胁迫。

3.3 磷酸戊糖代谢通路

磷酸戊糖代谢可以为细胞核酸代谢、脂肪酸与氨基酸的合成提供大量的NADPH，是植物糖代谢的重要途径[41]。翅碱蓬受到PAHs胁迫后，富集在戊糖磷酸途径的下调差异基因占比有所增加，可推测翅碱蓬受PAHs胁迫抑制了其磷酸戊糖代谢，减损细胞合成代谢所需诸多原料。通过KEGG注释，发现在A-vs-D对比组中磷酸戊糖通路中检测到4条基因下调，其中2条为D-ribulose-5-phosphate 3-epimerase(E5.1.3.1)，2条为果糖-1，6-二磷酸醛缩酶(E4.1.2.13)。

E5.1.3.1广泛存在于真核生物、真菌和细菌内，隶属ribulose phosphate binding家族[42]，是催化戊糖磷酸的非氧化分支d-核酮糖-5-磷酸和d-木酮糖-5-磷酸的相互转化的关键酶。关于E5.1.3.1的动力学研究较少，但有研究表明，该酶在某些生物体中可作为一种金属酶参与反应，并且已显示使用二价阳离子作为辅助因子或活化剂。作为一种金属酶，大肠杆菌中的E5.1.3.1被认为是 H2O2 产生的氧化应激的主要目标之一[43]，除此之外，D-ribulose-5-phosphate 3-epimerase还被认为具有调控木薯光合效率的功能，对其光合作用及产量产生影响[44]。E4.1.2.13属于裂解酶，是戊糖磷酸途径中裂解果糖1，6-二磷酸可逆反应的关键酶，以为生物体提供物质合成能量代谢所需的底物与能量的方式，影响生物体的生长发育[45]。E4.1.2.13在生物体内种类不同，在真菌、绿藻内多为Ⅱ型，同源二聚体构成，催化反应时不产生酶-底物中间物；植物或动物体内的多为Ⅰ型，同源四聚体构成，产生中间代谢物，EDTA与硼氢化物分别对上述两种类型的酶活性有抑制作用[46]。除此之外，Pontremoli等[47]研究发现，在正常进食与禁食的兔子肝脏中提取到的E4.1.2.13活性有较大差异，但该酶浓度差距不大，同时也指出在禁食期间该酶的修饰受到较大影响，导致该酶催化活性损失，进而造成免疫反应性损失。

此外，不同代谢通路在各浓度下的富集强度亦呈现动态变化。如乙醛酸与二羧酸代谢通路在低浓度组(A-vs-B)富集最显著(Q=0.006 9)，随浓度升高富集程度减弱，提示该通路可能在低浓度胁迫初期即被激活，参与能量补给与抗氧化防御；而碳代谢和磷酸戊糖通路在高浓度组(A-vs-D)仍显著富集，说明碳代谢与磷酸戊糖途径在高浓度PAHs胁迫下仍持续参与能量供应与还原力平衡，但其调控方向(上调/下调)可能已发生转变，反映出植物在不同胁迫强度下的代谢策略调整。

3.4 GO功能与代谢通路联系

本研究中，GO富集结果显示，不同浓度PAHs胁迫下翅碱蓬的关键细胞结构与分子功能发生特异性响应，这些响应与KEGG富集的代谢通路密切相关：

内质网膜(endoplasmic reticulum membrane) 在A-vs-D组显著富集，提示高浓度PAHs可能引起内质网应激，影响蛋白质合成、折叠与分泌过程，进而可能影响碳代谢通路中关键酶的合成与功能。

磷酸吡哆醛结合(pyridoxal phosphate binding) 同样在A-vs-D组富集，磷酸吡哆醛作为多种酶的辅因子，参与氨基酸代谢、磷脂合成等过程，其相关基因的上调可能与glyoxylate and dicarboxylate metabolism中丝氨酸/甘氨酸代谢的增强有关，以维持一碳单位供应与抗氧化代谢。

叶绿体(chloroplast) 在A-vs-C组富集显著，结合KEGG中carbon fixation in photosynthetic organisms通路的富集，说明中浓度PAHs胁迫下光合作用相关结构与功能受到显著影响，这可能直接关联到碳代谢通路的调控，尤其是与光合同化产物的分配与利用有关。

综上所述，GO富集结果从细胞结构、分子功能层面揭示了翅碱蓬在PAHs胁迫下的响应重点，这些功能与KEGG富集的三大代谢通路(乙醛酸代谢、碳代谢、磷酸戊糖途径)共同构成一个协同应对胁迫的调控网络：内质网参与蛋白代谢，叶绿体关联光合碳固定，磷酸吡哆醛则作为关键辅因子支持多条代谢途径的顺利进行。

4 结论

本研究中通过不同浓度PAHs对翅碱蓬进胁迫试验，利用高通量测序技术对翅碱蓬植株进行转录组测序分析，PAHs胁迫不利于翅碱蓬的生长代谢，同时植物也通过自身机制应对胁迫，主要得到以下结论：

1)高通量测序技术可以对PAHs胁迫下的翅碱蓬植株进行转录组测序分析，且测序质量好，数据分析可靠。

2)PAHs胁迫下，差异基因表达主要与内质网膜、磷酸吡哆醛结合、叶绿体、核酮糖-二磷酸羧化酶活性和光呼吸等相关。

3)不同浓度PAHs胁迫下，翅碱蓬主要通过乙醛酸和二羧酸代谢、碳代谢、磷酸戊糖3个通路产生响应PAHs的胁迫。翅碱乙醛酸和二羧酸代谢通路中编码甘氨酸脱羧酶复合物(GDC)与丝氨酸羟甲基转移酶(E2.1.2.1)的基因下调。碳代谢通路中编码Glucose-6P转化过程、核糖5-磷酸(Ribose-5P)与核酮糖5-磷酸(Ribulose-5P)相互转化过程相关基因上调。磷酸戊糖通路中 D-ribulose-5-phosphate 3-epimerase(E5.1.3.1)和果糖-1，6-二磷酸醛缩酶(E4.1.2.13)基因下调。

[1] HUANG W X,WANG Z Y,YAN W.Distribution and sources of polycyclic aromatic hydrocarbons (PAHs) in sediments from Zhanjiang Bay and Leizhou Bay,South China[J].Marine Pollution Bulletin,2012,64(9):1962-1969.

[2] COOK J W,HEWETT C L,HIEGER I.106.The isolation of a cancer-producing hydrocarbon from coal tar.Parts Ⅰ,Ⅱ,and Ⅲ[J].Journal of the Chemical Society (Resumed),1933:395.

[3] BLUMER M.Benzpyrenes in soil[J].Science,1961,134(3477):474-475.

[4] 张明亮.滨海盐沼湿地退化机制及生态修复技术研究进展[J].大连海洋大学学报,2022,37(4):539-549.ZHANG M L.Research advancement on degradation mechanism and ecological restoration technology of coastal salt-marsh:a review[J].Journal of Dalian Ocean University,2022,37(4):539-549.(in Chinese)

[5] LIU Q,YI T X,LI Q Y,et al.Bioaccumulation of heavy metals by Suaeda salsa in the tidal flat of the Liaohe Estuary[J].Separations,2022,9(11):374.

[6] 华正韬.土壤和石油烃性质对溶剂萃取法修复石油污染土壤的影响[D].天津:天津大学,2013.HUA Z T.The influence of soil and petroleum hydrocarbon properties on the efficiency of petroleum contaminated soil remediation by solvent extraction[D].Tianjin:Tianjin University,2013.(in Chinese)

[7] 廖书林.辽河口湿地土壤中多环芳烃的分布特征及来源解析[D].青岛:中国海洋大学,2011.LIAO S L.Distribution and sources apportionment of PAHs from Liaohe estuarine wetland soils[D].Qingdao:Ocean University of China,2011.(in Chinese)

[8] WANG D R,FENG Q R,ZHU S L,et al.Auxin-regulated PAHs absorption in wheat roots:insights from kinetics,enzymatic activity,and transcriptomic profiling[J].Journal of Agricultural and Food Chemistry,2025,73(36):22257-22271.

[9] 陆程张.Na2CO3胁迫下玉米苗期根部转录组表达谱分析[D].延吉:延边大学,2021.LU C Z.Analysis of transcriptome expression profile of maize seedling roots under Na2CO3 stress[D].Yanji:Yanbian University,2021.(in Chinese)

[10] ZENG N B,YANG Z J,ZHANG Z F,et al.Comparative transcriptome combined with proteome analyses revealed key factors involved in alfalfa (Medicago sativa) response to waterlogging stress[J].International Journal of Molecular Sciences,2019,20(6):1359.

[11] PARVAIZ A,SATYAWATI S.Salt stress and Phyto-biochemical responses of plants—a review[J].Plant,Soil and Environment,2008,54(3):89-99.

[12] ZHANG J L,LI J Y,WANG Y C,et al.Identification and analysis of differential genes for salt stress response in southern-type alfalfa mutant leaves[J].Journal of Agricultural Biotechnology,2017,25(10):1588-1599.

[13] YANG Y Q,GUO Y.Unraveling salt stress signaling in plants[J].Journal of Integrative Plant Biology,2018,60(9):796-804.

[14] 廖书林,郎印海,王延松.辽河口湿地土壤多环芳烃的分布与生态风险评价[J].环境化学,2011,30(2):423-429.LIAO S L,LANG Y H,WANG Y S.Distribution and ecological risk assessment of pahs in soils from Liaohe estuarine wetland[J].Environmental Chemistry,2011,30(2):423-429.(in Chinese)

[15] 张玉凤,吴金浩,宋永刚,等.辽东湾海水中PAHs分布与来源特征及风险评估[J].环境科学研究,2017,30(6):892-901.ZHANG Y F,WU J H,SONG Y G,et al.Distribution,sources and ecological risk assessment of polycyclic aromatic hydrocarbons in surface seawater in Liaodong Bay,China[J].Research of Environmental Sciences,2017,30(6):892-901.(in Chinese)

[16] MIZUTANI S,TEMIN H M.An RNA-dependent DNA polymerase in virions of Rous sarcoma virus[J].Cold Spring Harbor Symposia on Quantitative Biology,1970,35:847-849.

[17] GUBLER U,HOFFMAN B J.A simple and very efficient method for generating cDNA libraries[J].Gene,1983,25(2/3):263-269.

[18] LI B,DEWEY C N.RSEM:accurate transcript quantification from RNA-Seq data with or without a reference genome[J].BMC Bioinformatics,2011,12(1):323.

[19] LOVE M I,HUBER W,ANDERS S.Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2[J].Genome Biology,2014,15(12):550.

[20] ANDERS S,HUBER W.Differential expression analysis for sequence count data[J].Genome Biology,2010,11(10):R106.

[21] 苟娜娜,刘泽梁,吴蔓莉,等.总石油烃中当量烷烃和多环芳烃组分的毒性风险及联合毒性效应[J].生态毒理学报,2025,20(1):156-169.GOU N N,LIU Z L,WU M L,et al.Toxicological risks and combined toxic effects of equivalent alkanes and polycyclic aromatic hydrocarbon components in total petroleum hydro-carbons[J].Asian Journal of Ecotoxicology,2025,20(1):156-169.(in Chinese)

[22] 冯爽.石油烃污染对蚯蚓和玉米的毒性效应研究[D].西安:西安建筑科技大学,2024.FENG S.Study on toxic effects of petroleum hydrocarbon pollution on earthworm and corn[D].Xi’an:Xi’an University of Architecture and Technology,2024.(in Chinese)

[23] 刘菁.烟草维生素B6代谢酶磷酸吡哆醛磷酸酶基因的克隆与分析[D].合肥:安徽农业大学,2018.LIU J.Cloning and analysis of vitamin B6 metabolizing enzyme pyridoxal phosphate phosphatase gene in tobacco[D].Hefei:Anhui Agricultural University,2018.(in Chinese)

[24] CHEN H,XIONG L M.Pyridoxine is required for post-embryonic root development and tolerance to osmotic and oxidative stresses[J].The Plant Journal,2005,44(3):396-408.

[25] 张敏,黄利斌,蒋泽平,等.NaCl胁迫对构树幼苗内质网膜脂肪酸及蛋白质组成的影响[J].植物资源与环境学报,2009,18(2):9-14.ZHANG M,HUANG L B,JIANG Z P,et al.Effects of NaCl stress on fatty acid composition and protein composition in endoplasmic reticulum membrane of Broussonetia papyrifera[J].Journal of Plant Resources and Environment,2009,18(2):9-14.(in Chinese)

[26] 赵相山,陈国仓,张承烈.不同生境芦苇1,5-二磷酸核酮糖羧化酶活性与二氧化碳固定关系[J].兰州大学学报,1993,29(2):153-154.ZHAO X S,CHEN G C,ZHANG C L.Relationship between 1,5- rubp carboxylase activity and carbon dioxide fixation in Phragmites communis in different habitats[J].Journal of Lanzhou University,1993,29(2):153-154.(in Chinese)

[27] 周天骄,丁晓辉,王君晖.植物光呼吸途径的调控和优化策略[J].浙江大学学报(农业与生命科学版),2020,46(3):271-279.ZHOU T J,DING X H,WANG J H.Strategies for modulation and optimization of the photorespiration pathway in plants[J].Journal of Zhejiang University (Agriculture and Life Sciences Edition),2020,46(3):271-279.(in Chinese)

[28] DOUCE R,BOURGUIGNON J,NEUBURGER M,et al.The Glycine decarboxylase system:a fascinating complex[J].Trends in Plant Science,2001,6(4):167-176.

[29] 周鸿铭,雷娜,鲁亚平.甘氨酸神经递质研究进展[J].生物学杂志,2011,28(1):79-81.ZHOU H M,LEI N,LU Y P.Research advances on Glycine neurotransmitters[J].Journal of Biology,2011,28(1):79-81.(in Chinese)

[30] 陈帅,瞿南,曹守云,等.C3-C4中间型植物Flaveria anomala中甘氨酸脱羧酶P-蛋白亚基上游调控序列在转基因水稻中的表达[J].科学通报,2001,46(11):939-942,970.CHEN S,QU N,CAO S Y,et al.Expression of upstream regulatory sequence of Glycine decarboxylase P- protein subunit in C3-C4 intermediate plant Flaveria anomala in transgenic rice[J].Chinese Science Bulletin,2001,46(11):939-942,970.(in Chinese)

[31] 张椿雨.芸薹属作物与Moricandia nitens间甘氨酸脱羧酶P亚基基因的比较分析[D].武汉:华中农业大学,2005.ZHANG C Y.Comparative analysis of Glycine decarboxylase P subunit gene between Brassica crops and Moricandia nitens[D].Wuhan:Huazhong Agricultural University,2005.(in Chinese)

[32] 马莉,陈丽梅.植物丝氨酸羟甲基转移酶基因研究进展[J].生物技术通报,2008,24(2):15-19.MA L,CHEN L M.The research advances on serine hydroxymethyltransferase gene in plants[J].Biotechnology Bulletin,2008,24(2):15-19.(in Chinese)

[33] 林颖辉,王文磊,徐燕,等.坛紫菜丝氨酸羟甲基转移酶基因的克隆及表达特征[J].渔业科学进展,2018,39(5):122-129.LIN Y H,WANG W L,XU Y,et al.Cloning and expression analysis of serine hydroxyl methyltransferase (SHMT) genes from Pyropia haitanensis[J].Progress in Fishery Sciences,2018,39(5):122-129.(in Chinese)

[34] 马莉,陈丽梅,刘迪秋,等.植物丝氨酸羟甲基转移酶及其生理作用研究进展[J].安徽农业科学,2008,36(4):1357-1359,1404.MA L,CHEN L M,LIU D Q,et al.Research progresses on the molecular properties and physiological functions of plant serine hydroxymethyltransferase[J].Journal of Anhui Agricultural Sciences,2008,36(4):1357-1359,1404.(in Chinese)

[35] 王诗雅.初花期淹水胁迫下烯效唑对大豆碳代谢和产量的缓解效应[D].大庆:黑龙江八一农垦大学,2021.WANG S Y.Alleviation effect of uniconazole on carbon metabolism and yield of soybean under waterlogging stress at initial flowering stage[D].Daqing:Heilongjiang Bayi Agricultural University,2021.(in Chinese)

[36] 张冬冬.浑浊红球菌脂质积累机制的碳代谢流分析及关键基因改造[D].无锡:江南大学,2016.ZHANG D D.Carbon metabolic flux analysis of lipid accumulation mechanism and key genetic modification in Rhodococcus opacus[D].Wuxi:Jiangnan University,2016.(in Chinese)

[37] 李梦娇.灵芝碳代谢途径中UDP葡萄糖焦磷酸化酶和磷酸葡萄糖变位酶基因功能研究[D].南京:南京农业大学,2015.LI M J.Study on the function of UDP glucose pyrophosphorylase and phosphoglucomutase gene in carbon metabolism pathway of Ganoderma lucidum[D].Nanjing:Nanjing Agricultural University,2015.(in Chinese)

[38] WANG J,YANG W T.Concerted proton transfer mechanism of Clostridium thermocellum ribose-5-phosphate isomerase[J].The Journal of Physical Chemistry B,2013,117(32):9354-9361.

[39] 尚瑞沙,齐静茹,陈红丽,等.家蚕微孢子虫核糖-5-磷酸异构酶A基因的克隆及表达特征分析[J].蚕业科学,2019,45(1):61-66.SHANG R S,QI J R,CHEN H L,et al.Cloning and expression characteristics of ribose-5-phosphate isomerase a gene of Nosema bombycis[J].Acta Sericologica Sinica,2019,45(1):61-66.(in Chinese)

[40] WAMELINK M M C,GRÜNING N M,JANSEN E E W,et al.The difference between rare and exceptionally rare:molecular characterization of ribose 5-phosphate isomerase deficiency[J].Journal of Molecular Medicine,2010,88(9):931-939.

[41] 侯夫云.水稻戊糖磷酸途径两个关键酶基因的克隆与功能分析[D].南京:南京农业大学,2005.HOU F Y.Cloning and functional analysis of two key enzyme genes in pentose phosphate pathway from rice (Oryza sativa L.)[D].Nanjing:Nanjing Agricultural University,2005.(in Chinese)

[42] AKANA J,FEDOROV A A,FEDOROV E,et al.D-ribulose 5-phosphate 3-epimerase:functional and structural relationships to members of the ribulose-phosphate binding (β/α)8-barrel superfamily[J].Biochemistry,2006,45(8):2493-2503.

[43] GONZALEZ S N,VALSECCHI W M,MAUGERI D,et al.Structure,kinetic characterization and subcellular localization of the two ribulose 5-phosphate epimerase isoenzymes from Trypanosoma cruzi[J].PLoS One,2017,12(2):e0172405.

[44] 宋雁超,姚惠,吕亚,等.花叶木薯变种和木薯栽培种ZM-Seaside叶片光合参数及蛋白组学分析[J].植物遗传资源学报,2016,17(5):935-941.SONG Y C,YAO H,LÜ Y,et al.The analysis of photosynthetic parameters and proteomics of leaves from cassava (Manihot esculenta Crantz) mosaic-leaf mutation and cultivar ZM-seaside[J].Journal of Plant Genetic Resources,2016,17(5):935-941.(in Chinese)

[45] 唐功利,杨春松,鲍建绍,等.丙糖磷酸异构酶、果糖-1,6-二磷酸醛缩酶及果糖-1,6-二磷酸酶的共表达[J].生物化学与生物物理学报,2001(1):131-136.TANG G L,YANG C S,BAO J S,et al.Co-expression of triosephosphate isomerase,fructose-1,6-bisphosphate aldolase and fructose-1,6-bisphosphatase in E.coli[J].Acta Biochimica et Biophysica Sinica,2001(1):131-136.(in Chinese)

[46] 路玮.拟南芥果糖1,6-二磷酸醛缩酶家族分析[D].泰安:山东农业大学,2011.LU W.Genome-wide analysis of the fructose bisphosphate aldolases in Arabidopsis[D].Tai’an:Shandong Agricultural University,2011.(in Chinese)

[47] PONTREMOLI S,MELLONI E,SALAMINO F,et al.Changes in activity of fructose-1,6-bisphosphate aldolase in livers of fasted rabbits and accumulation of crossreacting immune material[J].Proceedings of the National Academy of Sciences of the United States of America,1979,76(12):6323-6325.