The chromosome-scale assembly of the Notopterygium incisum genome provides insight into the structural diversity of coumarins

羌活（Notopterygium incisum）基因组的染色体级别组装为香豆素的结构多样性提供了新的见解

摘要

香豆素是由苯丙素途径衍生而来的，是被子植物中的主要代谢物之一。四氢吡喃（THP）和四氢呋喃（THF）环与内酯结构的对齐形成至少四种复杂的香豆素。然而，香豆素结构多样性的机制尚不清楚。在本研究中，我们报告了羌活（Notopterygium incisum）基因组的染色体级别组装，总长度为1.64 Gb，contig N50值为22.7 Mb，注释到60,021个蛋白质编码基因。此外，我们鉴定了塑造香豆素结构多样性的关键酶，包括两个对简单香豆素基本骨架结构至关重要的对羟基肉桂酰辅酶A 2'-羟化酶，两个负责角型或线型香豆素生物合成的UbiA异戊烯转移酶，以及五个参与THP和THF环形成的CYP736环化酶。值得注意的是，我们首次鉴定了两个能够同时催化7-去甲基软木花椒素和王草酚的双功能酶。进化分析表明，CYP736亚家族的串联和异位重复，特别是在伞形科中出现，促进了羌活中香豆素结构多样性的形成。总结，本研究提出了一种不同被子植物中复杂香豆素生物合成途径的平行进化情景，并为异源工业生产香豆素提供了基本的合成生物学元素。

图形摘要

我们组装了高质量的羌活基因组，并验证了香豆素生物合成中涉及的关键基因。进化分析表明，CYP736的重复有助于香豆素结构的多样性。

关键词

羌活，平行进化，生物合成，吡喃香豆素，呋喃香豆素

引言

香豆素是高等植物中广泛存在的主要次级代谢产物之一。2H-1-苯并吡喃-2-酮核心结构及其四氢吡喃（THP）和四氢呋喃（THF）环是天然和合成生物活性分子中的常见结构单元。特别是，将THP和THF环与内酯核心结构对齐可形成至少四种基本THP和THF香豆素骨架的亚群。这导致了数千种香豆素化合物的形成，包括吡喃香豆素和呋喃香豆素。研究表明，香豆素表现出多种生理和医疗生物活性，包括抗真菌感染、防御素或诱导素、抗癌、抗炎、抗氧化和钙通道阻断等。例如，呋喃香豆素、氧化桔皮素和佛手柑内酯是美国食品药品监督管理局（FDA）批准用于治疗皮肤病的常用临床药物。这些多样的生物活性主要归因于其结构的多样性。然而，关于这些香豆素结构多样性的生物合成机制知之甚少。

香豆素衍生自苯丙素途径，在植物王国中具有保守且广泛分布的初始上游合成过程。该途径可通过苯丙氨酸解氨酶（PAL）、肉桂酸4-羟化酶（C4H）和4-香豆酸：辅酶A连接酶（4CL）生产关键前体，如肉桂酸、对香豆酸和对香豆酰-CoA。不同的酶参与分支对香豆酰-CoA形成不同的代谢物，如木质素、类黄酮和苯乙烯。在这些分支酶中，两种2-氧代戊二酸依赖的双加氧酶（2-OGD）家族蛋白，即对香豆酰CoA 2′-羟化酶（C2′H）和阿魏酰CoA 6′-羟化酶（F6′H），被认为参与了简单香豆素骨架（鬼臼素或东莨菪内酯）的形成。此外，从不同谱系的被子植物中克隆了四种C2′Hs/F6′Hs，包括十字花科、旋花科、芸香科和伞形科。随后，通过伞形素6/8-异戊烯转移酶（PT）在6或8位取代伞形素上的异戊烯基团，分别形成7-去甲基软木花椒素（DMS）或王草酚。这一步是复杂香豆素合成的入口点，决定了角型或线型香豆素的形成。然而，已鉴定的PTs仅限于桑科、芸香科和伞形科成员。这一现象与复杂香豆素主要存在于桑科、芸香科和伞形科，而简单香豆素广泛分布于至少75个科的事实一致。最后，7-去甲基软木花椒素或王草酚的环化形成相应的THP和THF环，这是塑造复杂香豆素基本骨架的关键步骤。先前的研究表明，这种环化是由CYP450蛋白催化的，其中桑科的特定CYP76F112蛋白参与线性呋喃香豆素环化。最近，我们发现两个CYP736家族蛋白也参与了线性/角型THP和THF环的形成。然而，植物生成香豆素结构多样性的过程以及不同被子植物使用不同酶元件合成相同化合物的原因仍不清楚。随后，羟基化、甲氧基化和糖基化等结构修饰步骤发生，完成复杂香豆素的后修饰。然而，到目前为止，仅鉴定出一种类型的羟基化或甲氧基化酶。因此，香豆素的生物合成机制以及扩展结构多样性的进化机制尚待阐明。

图1. 提议的香豆素生物合成途径。提议的香豆素生物合成途径包括三个阶段：简单香豆素，从简单香豆素到复杂香豆素，以及复杂香豆素如吡喃香豆素和呋喃香豆素。每个阶段都用不同的颜色标记。缩写：PAL，苯丙氨酸解氨酶；C4H，肉桂酸4-羟化酶；4CL，4-香豆酸：辅酶A连接酶；C2′H，对香豆酰CoA 2′-羟化酶；F6′H，阿魏酰CoA 6′-羟化酶；6 PT，伞形素6-异戊烯转移酶；8 PT，伞形素8-异戊烯转移酶；DC，7-去甲基软木花椒素环化酶；OC，王草酚环化酶。多个箭头表示多个已知或未知的步骤。

羌活（Notopterygium incisum Ting ex H. T. Chang）是一种重要的传统中药和藏药植物，具有悠久的历史背景（图2A及支持信息图S2）。它主要分布在中国青海、西藏、陕西、四川、甘肃等高原地区，海拔2000–4000米（支持信息图S3）。考虑到羌活丰富的香豆素含量和结构多样性，我们组装了该植物的高质量基因组，以探索香豆素合成途径的生物合成和进化机制。通过结合基因组学、代谢组学和转录组学分析，我们鉴定了所有参与简单和复杂香豆素骨架生物合成的关键基因。因此，本研究为复杂香豆素的异源工业生产提供了基础。此外，它提供了该物种的详细遗传信息，并探讨了羌活和被子植物中香豆素结构多样性的潜在因素。

图2. 羌活的基因组组装、系统发育树重建和全基因组复制分析。 (A) 羌活的形态。 (B) 羌活的基因组特征。圆形图从外到内的圆圈分别表示染色体级别的伪染色体（Chr1-Chr11）(i)、基因密度(ii)、重复序列密度(iii)、GC含量(iv)，圆形图中心的每条连线表示一对同源基因(v)。 (C) 羌活与其他33种植物之间的系统发育关系和分化时间。物种树是使用1730个低拷贝同源基因基于最大似然法生成的，支持率很高。 (D) 羌活中锚定的旁系同源基因对和全基因组重复（WGD）事件的Ks分布。羌活的Ks分布显示了两个峰值，一个在约0.30（α-WGD），另一个在约0.83（β-WGD）。虚线表示在校正进化速率后基因组间的Ks分布峰值。所有WGD和物种分化事件的Ks峰值都被标注出来。

材料和方法

2.1. 实验材料

羌活、川芎和栽培的羌活样本采自青海省贵德（北纬36°21′34″，东经101°32′6″），并通过DNA分子鉴定确认。一个凭证标本（编号CPUZYC2021013）存放在中国药科大学植物园，南京，中国。基因组DNA从羌活的幼叶中提取用于测序。羌活的叶、果、根和茎，以及川芎和栽培羌活的根用于转录组测序。用于提取和生化分析的溶剂为色谱级。用于酶反应的化合物和用于微生物培养的抗生素购自Herbest（宝鸡，中国）和Sigma-Aldrich（上海，中国）。用于基因片段扩增的酶、试剂盒和体外表达载体购自TransGen Biotech（北京，中国）。

2.2. 基因组大小估计、基因组测序、组装和注释

基因组大小基于72.5 Gb Illumina双端短读序列的k-mer分布分析估计。使用Jellyfish 2.2.6计算k-mer频率，k-mer大小为19。提取的羌活DNA在PacBio Sequel II平台上测序，生成38.93 Gb HiFi长读序列数据，平均长度为18.42 Kb。经过SMRTLink v8.0的初步质量控制后，使用hifiasm（v0.16.1）首先获得了羌活核基因组的de novo组装的contig级别组装。采用HiC技术辅助羌活contig的搭建，总共使用了387.7 Gb HiC数据，经过FastQC（v0.11.9）过滤。最后，使用AllHiC（v0.9.8）管道将抛光和去除haplotigs的contigs大多锚定到11个伪分子上。

羌活基因组中的蛋白质编码基因通过转录组、同源性和ab initio预测方法相结合进行注释。首先，重复序列注释结合基于Repbase库的同源预测和基于自序列比对的de novo预测。在这两种方法中，使用RepeatModeler（版本：open-1.0.11）和RepeatMasker（版本：open-4.0.9）分别构建de novo重复序列数据库并从基因组和Repbase库中搜索重复序列。在掩蔽重复序列后，使用GeneMarkS-T（v5.1）通过从RNA-Seq读序列组装的转录本预测基因模型。RNA-Seq读序列通过HISAT2（v2.0.4）映射到掩蔽的基因组上，并使用StringTie（v1.2.3）组装成转录本。然后，同源性方法使用GeMoMa（v1.7）结合五个物种的参考基因模型，包括胡萝卜、香菜、芹菜、三七和人参。之后，使用Augustus（v2.4）和SNAP（2006-07-28）预测de novo基因模型。最后，使用EVidenceModeler（v1.1.1）整合所有三种方法的基因模型，并通过PASA（v2.0.2）更新。

通过使用Diamond（v0.9.24）针对NCBI NR（202009）、SwissProt/TrEMBL（202005）和EggNOG（v5.0）数据库进行blast进行蛋白质编码基因的功能注释。使用HMMER（v3.1）在Pfam（v33.1）数据库下搜索蛋白质结构域。使用LTRharvest和LTR_finder鉴定末端重复反转录转座子。串联重复序列通过TRF（v4.09.1）和MISA（v2.1）注释。

2.3. 全基因组复制（WGD）和进化速率校正

WGD事件的鉴定依赖于使用wgd（v1.1）为羌活、葡萄和其他三个伞形目物种的所有旁系同源物（paranome）构建基于Ks的年龄分布。在WGD管道中，编码的核苷酸序列首先翻译成肽序列，然后根据要求输入BLASTP（v2.13.0）进行全对全blast，E值设置为1×10^-10。随后，调用MCL（v14-137）和MUSCLE（v5.1）获得每个旁系同源基因家族的蛋白质水平的多序列比对（MSA）。在将此MSA反向翻译为密码子比对后，使用PAML（v4.9）中的CODEML包在默认控制文件下获取Ks值的最大似然估计，并最终通过结合每个锚对的权重（无异常值）构建基于Ks的年龄分布。对于直系同源基因对（orthologous Ks-based age distributions），MCL聚类过程被替换为相互最佳匹配（RBH）搜索以识别直系同源基因对，并仅推断一对一的直系同源基因对。此外，通过使用Ksrates（v1.1.1）校正不同物种之间的进化速率差异。

2.4. 系统发育树构建与分子年代测定

为了全面阐明伞形目 (Apiales) 内的系统发育关系，我们选择了34个已报道的被子植物基因组和转录组（参见支持信息表 S1），并据此重建了系统发育树。我们使用 Orthofinder (v2.5.5) 对这些34个物种的肽序列数据集进行直系同源群的分类。筛选出1730个低拷贝的直系同源基因家族（拷贝数 <3，覆盖率 >90%），并从每个物种中选择这些家族中最长的拷贝作为精炼的1730个低拷贝直系同源基因，用于系统发育分析。通过去冗余脚本去除直系同源群中的冗余拷贝，并将其输入 Gblock (v0.91b) 以获得保守序列。随后，将所有保守序列串联成一个超级基因。在构建系统发育树之前，使用FastTree (v2.1.11) 对数据矩阵进行评估并选择最佳模型。接着，使用 JTT 模型进行数据集的系统发育分析。MCMCtree 使用 PAML (v4.9) 包对基于1730个串联基因和15个基于化石的年代校准 (参见支持信息表 S2) 所得的树形结构进行物种分化时间的估计。GTR 模型被用作肽替换模型。节点年龄的后验分布通过马尔科夫链蒙特卡罗采样估计，每10步抽样一次，共100,000步，随后进行400,000步的抽样。

2.5. 代谢组学分析

我们共选择了42个样品进行代谢组学分析，包括7个来自 羌活 根、茎、叶和果实的生物重复样本，以及 宽叶羌活 根和栽培的 羌活 根样本。这些样品用预冷的80%甲醇进行提取，并将提取物进行超高效液相色谱 - 四极杆飞行时间质谱 (UPLC-Q-TOF-MS) 分析。化合物的色谱分离使用了 ACQUITY UPLC 系统（沃特世公司，马萨诸塞州，USA）。采用 C18 反相柱（50 mm × 2.1 mm，1.5 μm，赛默飞世尔科技，美国）进行UPLC分析，柱温为30°C。梯度包括0.1%甲酸（v/v，A）和乙腈（B），其梯度为：0分钟，10%B；3分钟，15%B；8分钟，38%B；12分钟，50%B；16分钟，95%B。流速为0.4 mL/min。质谱分析使用 Synapt G2-Si Q-TOF（沃特世质谱技术，曼彻斯特，UK）。数据采集使用 MassLynx V4.2 软件进行，后续数据处理包括去除背景噪音、使用参考样本进行标准化、保留时间校正和峰对齐等操作均使用 Progenesis QI V2.0（沃特世公司，马萨诸塞州，USA）进行。代谢物的注释使用了公共数据库如 METLIN ( https://metlin.scripps.edu/index.php)，Lipidmaps 数据库 (http://www.lipid_maps.org/)，KEGG 数据库 ( KEGG: Kyoto Encyclopedia of Genes and Genomes) 和 HMDB 数据库 (Human Metabolome Database)。SIMCA 14.1 软件（Umetrics，乌梅奥，瑞典）用于主成分分析 (PCA)。代谢物的鉴定和定量按我们之前的报道进行。

2.6. 香豆素生物合成相关基因的筛选

为了筛选与简单香豆素骨架形成相关的候选基因，我们筛选了基因组的功能注释信息，并选择了已识别的伞形科植物基因的蛋白质序列。这些序列然后进行 BLASTPv2.11.0 分析 (E-value < 1e‒5)。芳香化合物的PT属于 UbiA 家族，并与负责维生素合成的 VTE2-1 进行聚类分析。我们使用 HMMER (版本 3.1b2) 扫描所有预测的 羌活 中的 PT 基因，保留那些含有 "UbiA 异戊烯转移酶家族" (PF01040) 域的基因。候选 PTs 通过系统发育分析进行筛选，与拟南芥中的 PT 基因和其他植物中已验证功能的蛋白质序列进行对比。羌活 的基因组、转录组和代谢组数据被用于构建转录组-代谢组共表达矩阵。基于皮尔森相关系数算法开发了一个相关网络。基因共表达理论和加权相关网络分析 (WGCNA) 被用来分析候选的环化酶基因，这些基因可能与 PT 或 C2'H 有相似的表达模式。

2.7. 在大肠杆菌中的异源表达及活性测定

从 羌活 的cDNA中克隆了候选的C2'H 和 PT 基因，然后分别插入到 pET28a（在 BamH I 和 EcoR I 限制性位点之间）和 pETDuet-1（在 BamH I 和 Sac I 限制性位点之间）载体中。成功测序的重组质粒转化到 E. coli BL21 (DE3) 受体细胞中。细胞在含有卡那霉素 (50 mg/L) 和氨苄青霉素 (100 mg/L) 的 Luria-Bertani (LB) 培养基中培养。用 IPTG (0.5 mmol/L) 诱导后，收获细胞并通过超声裂解用于活性测试或蛋白质纯化。由于PT在大肠杆菌中的膜蛋白特性及其不溶性，直接使用粗酶液进行酶活性测试。相比之下，C2'H以可溶形式表达，允许使用 Ni-NTA 树脂进行成功纯化，并使用 SDL-030-F2（苏州盛德尔仪器有限公司）蛋白纯化系统进一步纯化，获得纯化的重组蛋白用于酶活性测试。所有的体外酶活性测定均在30°C的振荡培养箱 (220 rpm) 上进行2小时。对于PT活性测试，反应系统包含50 μL粗蛋白，200 μmol/L伞形酮，100 μmol/L DMAPP，和200 μmol/L MgCl2 在200 μL的100 mmol/L Tris-HCl (pH 7.5) 中。对于C2'H活性测试，酶反应（200 μL）包含10 μg纯化蛋白，100 mmol/L Tris-HCl (pH 6.5)，0.5 mmol/L FeSO4，5 mmol/L抗坏血酸钠，5 mmol/L 2-氧代戊二酸 (2OG)，和1 mmol/L p-香豆酰辅酶A或阿魏酰辅酶A。反应样品用乙酸乙酯提取并溶解在甲醇中用于HPLC和LC-MS分析。

2.8. NiPTs 在烟草中的功能表达

通过瞬时表达在烟草 (Nicotiana benthamiana) 中进一步功能表征 NiPTs。将 NiPTs 的开放阅读框 (ORF) 通过 AgeI 和 XhoI 限制性酶切位点引入到 pEAQ 二元质粒中。正确测序的重组质粒通过冻融法转化到农杆菌 A. tumefaciens 株 GV3101 中。阳性转化子在含选择性 LB 琼脂平板 (50 μg/mL 卡那霉素和50 μg/mL 利福平) 中于28°C培养。阳性转化子接种到10 mL液体LB培养基中，在28°C下摇动1天。随后通过6,000 rpm离心5分钟收集细胞。然后将细胞沉淀重悬在农杆菌诱导培养基（10 mmol/L MES，10 mmol/L MgCl2，100 μmol/L 乙酰丁香酮，pH 5.8）中，在25°C下孵育1小时。通过测量其光密度 (OD600) 确定细胞悬浮液的浓度，使用无针注射器将这些细菌悬浮液注入到4-5周龄烟草叶片的背面。三天后，在感染区域注射100 μL的伞形酮 (100 μmol/L PBS 溶液)。一天后，收集已注射底物的叶片，用甲醇提取后进行 LC-MS 分析。用含有空载体的 A. tumefaciens 侵染的 N. benthamiana 用作阴性对照。每个实验包括至少三株植物作为平行实验。

2.9. 在酿酒酵母中的异源表达及活性测定

我们使用整合了拟南芥 CYP450 还原酶 (AtCPR) 的酿酒酵母 WAT11 筛选候选的 CYP450 蛋白。携带 CYP450 基因的 pYES2-Ura 载体转化到酿酒酵母株 WAT11 中，并在 30°C 的 SD 培养基中培养至 OD600 达到 0.8–1.0。然后收集细胞并用无菌水洗涤。将细胞洗涤三次后用含半乳糖的 SG 培养基诱导蛋白表达。将诱导的酵母细胞收集并重悬于 Tris-HCl 缓冲液中进行全细胞催化。所有反应均在30°C 的振荡培养箱 (220 rpm) 中进行2小时。反应系统包含500 μmol/L NADPH 和 1 mmol/L 底物在200 μL的100 mmol/L Tris-HCl (pH 7.5) 中。反应样品用乙酸乙酯提取并溶解在甲醇中用于HPLC和LC-MS分析。

2.10. HPLC 和 LC-MS 分析

使用岛津LC-2010C系统（岛津，日本）进行HPLC分析。样品在 Hedera ODS-2 C18 柱（4.6 mm × 250 mm，5 μm）上分离。流动相为含0.1%甲酸（v/v，A）的水和甲醇（B）。流速为0.5 mL/min，柱温为30°C。对于所有酶促产物，检测波长为340 nm。检测C2'H活性的梯度洗脱程序如下：0分钟，10%B；5分钟，15%B；15分钟，60%B；22分钟，60%B。检测PT活性的梯度洗脱程序如下：0分钟，40%B；5分钟，65%B；15分钟，75%B；22分钟，95%B。检测环化酶活性的梯度洗脱程序如下：0分钟，30%B；5分钟，65%B；12分钟，95%B；14分钟，95%B；16分钟，30%；25分钟，30%。质谱分析使用配备加热ESI源的安捷伦6545 LC/Q-TOF（安捷伦技术，美国）进行。参数如下：MS扫描范围100-2000 Da；鞘气温度350°C；离子喷射电压3500 V；碰撞能量44 V。MassHunter Qualitative Analysis 软件用于观察和处理 LC-MS 数据。

3. 结果

3.1. 基因组测序、组装与注释

为了阐明 羌活 中复杂香豆素的生物合成途径，我们首先组装了 羌活 的核基因组，结合了约 72.1 Gb 的 Illumina 短读序列、约 38.93 Gb 的 PacBio HiFi 读序列和约 387.71 Gb 的 HiC 数据（参见支持信息表 S3）。在对两个单倍型进行去冗余处理后，最终组装的基因组总长度为 1.64 Gb，略小于通过 k-mer 分析估算的 1.75 Gb（参见支持信息图 S4 和表 S4）。总体上，97.94% 的重叠群被锚定到一个单倍型的 11 个拟染色体上（图 2B，参见支持信息图 S5 和表 S5、S6），与染色体计数数据库 (CCDB, https://ccdb.tau.ac.il/) 记录的染色体数目（2n = 2x = 22）相对应。

为了评估组装的完整性，我们使用了基准通用单拷贝直系同源基因 (BUSCO) 对基因组进行了比对，评估结果显示，组装的胚胎植物_odb10 基因的完整性为 98.7%（参见支持信息表 S7 和 S8）。通过 Ab initio、基于同源性和转录组辅助注释的方法，我们注释了 60,021 个蛋白质编码基因。这些基因与七个数据库中的蛋白质序列进行了比较，包括 KEGG、Pathway、Nr、UniProt、GO、Pfam 和 InterPro，其中 93.03% 的基因被分配了推测的功能注释（参见支持信息表 S9）。此外，注释还表明 羌活 含有高达 75.79% 的重复序列。其中，长末端重复逆转座子是最常见的转座元件类型，占基因组的 44.86%，包括 12.35% 的 Gypsy 和 32.38% 的 Copia 逆转座子（参见支持信息表 S10）。

3.2. 羌活 的系统发育位置与多倍化事件

为了精确确定 羌活 与伞形科其他物种之间的进化关系，我们使用 羌活 和其他 33 种被子植物的 1730 个精炼的低拷贝直系同源基因（参见表 S1）基于最大似然法（ML）构建了系统发育关系（图 2C）。羌活 被解析为伞形亚科内的一个早期分支，继 Chamaesium paradoxum 和 Bupleurum scorzonerifolium 分化后出现。此外，羌活 被观察到与剩余的伞形亚科分类群是姐妹关系。利用这一解析的系统发育树和 15 个化石校准点（参见表 S2），我们推测伞形科的冠群起源于约 4900 万年前 (MYA)，而 羌活 大约在 2200 万年前分化（参见支持信息表 S11），可能代表了 Physospermopsis 类群的出现时间。

全基因组重复事件被认为是基因组进化的关键驱动因素。为了识别发生在 羌活 基因组中的全基因组重复 (WGD) 事件，我们使用其他伞形目物种和葡萄 (V. vinifera) 作为参考进行了全面分析。通过计算 羌活 和伞形目其他八个已报告基因组中共线基因组块锚定对中的平行基因的同义替换率 (Ks)，我们在 羌活 中观察到了明显的 Ks 密度分布，显示出两个显著的峰值，分别在 0.37 和 0.83（图 2D）。值得注意的是，这些峰值在所有其他伞形科分类群中也一致出现（参见支持信息图 S6‒S8）。尽管由于进化速率的不同，伞形科不同植物表现出的这两个峰值位置略有不同，但调整后的物种分化峰值（0.17）证实了伞形科内的物种分化发生在这两个 WGD 事件之后。结果支持伞形科植物在被子植物共享的全基因组三倍化事件 (γ-WGT) 之后，经历了两轮额外的共享全基因组重复，分别称为 α-WGD 和 β-WGD。然而，对于这两个 WGD 事件的发生时间和涉及的分类群，尤其是它们是否发生在伞形科和五加科的共同祖先中，还是特定于伞形科，存在着不同的观点。我们的分析表明，五加科植物 Aralia elata 中代表 WGD 的 Ks 峰值（0.37）远小于对应伞形科和五加科分化的峰值（0.83），表明伞形科和五加科的分化早于在这两个家族中检测到的所有 WGD，并且这两个家族自 γ-WGT 以来没有共享的 WGD。因此，在 羌活 中检测到的两个 WGD 事件可能特定于伞形科。

3.3. 羌活 主要化合物的代谢物谱

最初进行了非靶向代谢组学分析，以表征 羌活 和 宽叶羌活 中香豆素的代谢物谱，这两种植物是中国药典的主要来源。样品包括 羌活 的叶子、果实、根和茎，以及 宽叶羌活 的根和栽培的 羌活 根（参见支持信息图 S9）。对 羌活 样品的主成分分析 (PCA) 表明，地下部分（根）与地上部分（叶、果实和茎）明显不同（参见支持信息图 S10）。此外，根部在总离子色谱中显示出更多的峰值和更高的丰度，暗示根部可能特异性地积累了更多的香豆素（参见支持信息图 S11）。

随后进行了靶向代谢组学分析，以鉴定 羌活 中参与香豆素生物合成的候选代谢物。通过与参考标准（tR 和 MS 数据）的比较或与理论数据或商业库的匹配，我们从 羌活 的所有组织中鉴定出39种化合物（图 3A，参见支持信息图 S12 和表 S12），其中20种为香豆素。因此，我们选择了 羌活 中的九种代表性香豆素，即当归素、当归醇、藁本内酯、异欧前胡素、反式双烯香豆素、隐黄素、柠檬素、佛手柑内酯和补骨脂素，作为标志化合物，以研究它们在不同组别中的丰度（参见支持信息图 S13）。如图 3B 所示，大多数标志化合物在根部的丰度较高，这与传统中药使用 羌活 根部作为药用成分的事实一致。此外，野生 羌活 中标志化合物的总体丰度高于栽培的 羌活（图 3C 和参见支持信息图 S14）。这可能解释了 羌活 能够成功栽培的原因。然而，人们仍然追求野生 羌活，因为其质量较高。我们还发现，一些化合物如当归素和异欧前胡素在 宽叶羌活 中含量较高（图 3C），这可能解释了为什么它们被选择作为 羌活 的替代品。然而，在栽培的 羌活 中检测到的所有化合物的含量均低于野生 羌活，这表明在习惯化栽培过程中，模仿野生策略可能是一种获得高质量药材的好方法（图 3B）。

图 3. 羌活 中主要化合物的代谢物谱。(A) 羌活 中鉴定出的化合物的化学结构。(B) 羌活 不同组织中九种主要香豆素的含量。R：根，S：茎，L：叶，F：果实。(C) 羌活 (Ni)、宽叶羌活 (Nf) 和栽培的 羌活 (Nif) 中九种主要香豆素的含量。误差条表示三个生物独立样本 (n = 3) 的 ±SEM 平均值。

3.4. 简单香豆素生物合成相关基因的鉴定

简单香豆素的生物合成机制相对明确 ，因此我们鉴定了54个可能编码参与简单香豆素生物合成酶的基因，包括5个PAL基因、1个C4H基因、25个4CL基因、12个C2′H基因和11个PT基因（参见支持信息图S15）。鉴于2-OGD和PT酶在推动从苯丙素途径到香豆素生物合成途径中的关键作用，我们主要集中在识别 羌活 中这两类基因。基于 羌活 中香豆素的代谢物谱，我们推测在根部高表达的NincChr3G00107490.1和NincChr11G00498940.1基因最有可能参与香豆素的生物合成（图3和图S15）。此外，代谢组-转录组相关性分析表明，这两个基因与 羌活 的主要香豆素产物异欧前胡素显著相关（图4A）。因此，我们选择了这两个基因进行功能评估。在重组蛋白表达、纯化（参见支持信息图S16）和活性测试后，发现只有NincChr3G00107490.1能够将对香豆酰辅酶A羟基化，形成简单的香豆素骨架——伞形酮（图4B和支持信息图S17A）。然而，当我们用2-OGD家族蛋白的另一种潜在底物（阿魏酰辅酶A）测试其活性时，NincChr11G00498940.1显示出相应的活性并形成了东莨菪素（图4C和图S17B）。因此，我们分别将NincChr3G00107490.1和NincChr11G00498940.1命名为NiC2′H和NiF6′H。NiC2′H和NiF6′H与所有已鉴定的2-OGD基因的氨基酸序列比对表明，这两个基因含有高度保守的Fe(II)结合基序His-X-Asp-XnHis（His233, Asp235和His291）和2-氧代戊二酸C5羧基结合基序Arg-X-Ser（Arg301和Ser303）（参见支持信息图S18）。系统发育分析表明，NiC2′H和NiF6′H属于先前描述的DOXC30亚类（参见支持信息图S19）。这两个基因与我们之前表征的PpC2′H聚类，显示出高度的序列相似性，暗示C2′H在伞形科中具有直系同源关系和功能一致性。

图 4. 羌活 中 2-OGD 和 PT 基因的功能验证。(A) 羌活 的代谢组-转录组相关性分析。列出了简单香豆素生物合成中的候选基因。(B, C) 羌活 中 2-OGD 基因的功能验证。(D, E) 羌活 中 PT 基因的功能验证。使用液相色谱 (LC, 340 nm) 和正离子化模式下的 LC‒MS 检测产物。图中显示了标准伞形酮、东莨菪素、DMS（去甲基藁本内酯）的 LC 图谱及其对应的化学结构。煮沸的酶用作对照。红色框标示了分子离子峰。

近一半的 PT 基因 (6/13) 与 羌活 中九种主要香豆素中的八种显著相关 (图 4A)，这表明 PT 基因在香豆素生物合成中的重要性。由于所有已鉴定的 PT 基因的编码序列 (CDS) 限制在约 1200 bp，因此选择了三个高度相关的 PT 基因 (3/6, NincChr3G00107340.1, NincChr3G00107500.1 和 NincChr3G00107320.1) 进行功能验证 (参见图 S15)。这些基因被克隆到原核表达载体中，按照我们之前描述的方法测试其活性。正如图 4D 所示，只有 NincChr3G00107340.1 和 NincChr3G00107500.1 在与标准去甲基藁本内酯（C6-伞形酮的异戊烯化产物）相同的保留时间上显示出显著的峰值。产物与标准去甲基藁本内酯的质子化分子离子 m/z 231.10 相同，与 [M+H]+ 的分子量一致 (图 4E)。因此，我们分别将 NincChr3G00107340.1 和 NincChr3G00107500.1 命名为 Ni6PT1 和 Ni6PT2。此外，我们在 N. benthamiana 中验证了 Ni6PT1 和 Ni6PT2 的活性。当伞形酮用作底物时，生成了一个新峰，并与标准样品直接比较，表明产生了 DMS (参见支持信息图 S20)。我们观察到，所有已鉴定的香豆素特异性 PT 基因的活性较低，这可能归因于它们的膜蛋白特性，尤其是在 PT 的 CDS 中夹杂至少六个跨膜区 (参见支持信息图 S21)。Ni6PT1 和 Ni6PT2 均属于伞形科特异的 UbiA PT 分支，并且与桑科分支 (FcPT1a 和 FcPT1b) 中的 PT 基因不同，这表明尽管两者都属于 UbiA 家族，但这两个 PT 组来自不同的祖先 (参见支持信息图 S22)。

3.5. 复杂香豆素生物合成中关键基因的鉴定

PT 产物（去甲基藁本内酯和隐黄素）的环化对于复杂香豆素的生物合成至关重要，这最终导致两类复杂香豆素的形成：呋喃香豆素和吡喃香豆素。复杂香豆素可根据其线性和角状结构分为四类：线性呋喃香豆素、线性吡喃香豆素、角状呋喃香豆素、角状吡喃香豆素。然而，目前仅鉴定出少数负责该最终步骤的环化酶。其中一种是来自桑科的 CYP76 家族蛋白，另一种是来自伞形科的 CYP736 家族蛋白。鉴于 羌活 也属于伞形科，我们推测 羌活 中的环化酶也可能属于 CYP736 蛋白。因此，我们鉴定了 羌活 中的 44 种 CYP736 蛋白作为潜在候选者。我们首先排除了一些具有较短 CDS (<1000 bp) 的蛋白质。基于 羌活 中主要化合物的代谢物谱（图 3），选择了七个在根中高表达的基因进行功能验证（图 5A）。经过活性测试后，我们发现四个基因（NincChr3G00106470.1、NincChr3G00121510.1、NincChr6G00302030.1 和 NincChr6G00302100.1）对去甲基藁本内酯表现出酶促活性，其保留时间与脱氢香豆素和马尔梅辛相同（图 5B 和支持信息图 S23）。此外，我们使用隐黄素测试了它们的活性，两个基因 NincChr6G00302030.1 和 NincChr6G00302100.1 显示出相应的活性（图 5C）。因此，NincChr3G00106470.1 和 NincChr3G00121510.1 被鉴定为单功能酶，以去甲基藁本内酯为底物，而 NincChr6G00302030.1 和 NincChr6G00302100.1 则是能够催化去甲基藁本内酯和隐黄素的双功能酶。为了确保包括其他潜在的环化酶，克隆了如 NincChr6G00302130.1、NincChr6G00302140.1、NincChr6G00302160.1 和 NincChr7G00344610.1 等低表达基因，并测试了它们的活性。有趣的是，NincChr7G00344610.1 对隐黄素表现出单功能酶促活性，并且只能在酵母中产生哥伦比亚素（图 5C 和支持信息图 S24）。因此，NincChr3G00106470.1 和 NincChr3G00121510.1 被命名为去甲基藁本内酯环化酶（NiDC1 和 NiDC2）；NincChr7G00344610.1 被命名为隐黄素环化酶（NiOC1）；而 NincChr6G00302030.1 和 NincChr6G00302100.1 分别被命名为隐黄素/去甲基藁本内酯环化酶（NiOD1 和 NiOD2）。值得注意的是，NiODs 是迄今为止首次报道的具有双重功能的酶。这些酶填补了复杂香豆素骨架生物合成中的缺失步骤，并为复杂香豆素的合成生物学提供了重要的酶元素。

图 5. 羌活 环化酶的功能验证。(A) 热图显示了候选基因在 羌活 不同组织中的相对表达水平。本文中验证活性的基因用红色三角形标记。R：根，S：茎，L：叶，F：果实。(B) 在 340 nm 下的去甲基藁本内酯环化酶活性的 LC 分析。(C) 在 340 nm 下的隐黄素环化酶活性的 LC 分析。在此工作中鉴定的 NiODs 以玫红色显示，NiDCs 以绿色显示，NiOC 以蓝色显示，对照（煮沸的酶）以黑色显示。标准马尔梅辛、脱氢香豆素、哥伦比亚素的 LC 图谱及其对应的化学结构显示在补充图中。

3.6. 羌活 中 CYP736 家族的进化及环化酶活性的出现

趋同或平行进化被认为可以解释不同植物独立获得呋喃香豆素的现象。来自桑科的环化酶是一种 CYP76 蛋白，而本研究中从伞形科鉴定的环化酶被分类为 CYP736 蛋白（图 5），进一步支持了这一假设。我们通过结合系统发育和比较基因组学的方法探讨了这些环化酶的进化历史。此外，考虑到 羌活 的环化酶表现出多样的活性（NiDC、NiOD 和 NiOC），我们还研究了这些环化酶功能分化的机制。

通过构建包含 18 个物种的所有 CYP450 蛋白的系统发育树（参见支持信息图 S25 和表 S13），我们研究了 羌活 环化酶的进化位置及其与其他 CYP450 蛋白的关系。系统发育分析将 羌活 的环化酶置于一个单系分支中，该分支包含至少 15 个 羌活 的 CYP736 蛋白（图 6A）。对该分支的深入研究表明，所有基因均属于伞形科植物，表明伞形科特有的基因复制导致了这五个环化酶基因的出现。在这15个CYP736蛋白中，绝大多数来自于3号和6号染色体，而常见的拓扑结构是3号染色体上的基因总是与 羌活 的6号染色体上的基因形成姐妹群。羌活 6号染色体上CYP736基因与其他物种的基因也形成了群组，这表明这些基因在伞形科物种之间具有祖先同源性，而 羌活 3号染色体上的CYP736基因仅与6号染色体上的CYP736基因聚类。这种拓扑结构表明，羌活 6号染色体上的CYP736基因可能比3号染色体上的平行基因出现得更早，而3号染色体上的CYP736基因可能是6号染色体上基因的重复拷贝。

图 6. 羌活 中 CYP736 基因家族的进化。(A) 伞形科中 CYP736 基因家族的系统发育树。羌活 的十五个基因分别来自 3 号染色体（用蓝色标记的五个基因）、6 号染色体（用粉色标记的九个基因）和 7 号染色体（用绿色标记的一个基因）。蓝色扇区表明 羌活 3 号和 6 号染色体上的 CYP736 基因是姐妹关系的平行基因。同时，粉色扇区显示 羌活 6 号染色体上的 CYP736 基因未与 3 号染色体上的基因群组，而是与伞形科其他物种的基因聚类。(B) 羌活 中 CYP736 基因家族的推测进化过程。虚线、橙色箭头和深蓝色箭头分别表示串联复制、异位复制和 WGD（片段）复制事件。数字 1–4 显示了复制事件期间或之后的氨基酸突变，这些突变可能与 羌活 中不同环化酶（OC、DC 和 OD）的功能分化有关。

这一复制事件可能是覆盖多个基因组块的异位复制，如系统发育分析和两个位点之间缺乏共线性所示（参见支持信息图 6B）。此外，推测在祖先染色体（6 号染色体）出现之前至少发生了一次串联复制事件，因为仅观察到两个包含 3 号和 6 号染色体基因的单系群。在异位复制之后，串联复制可能发生在 3 号和 6 号染色体上。这一系列复杂的复制事件最终导致在这两条染色体上发现了四个已鉴定的环化酶基因。位于 7 号染色体上的第五个环化酶可能是由于 WGD/染色体片段复制导致的，因为该基因与 6 号染色体上的环化酶基因具有良好的共线性（参见图 S26）。推测这一复制可能发生在 3 号和 6 号染色体之间的异位复制之前，如在 CYP736 系统发育树中的基底位置所示（图 6A）。

推测了五个环化酶的进化历史，NiDCs 通过 WGD/染色体片段复制产生，NiOCs 通过异位复制产生，而 NiODs 通过 6 号染色体上的串联复制产生（图 6B）。然而，导致酶活性分化的机制，如序列突变，仍需进一步研究。通过比较 15 个 羌活 CYP736 的氨基酸序列（参见支持信息图 S27），我们使用四个系统发育上早期分化的 羌活 CYP736 作为参考（NincChr3G00105000.1、NincChr3G00105010.1、NincChr3G00105660.1 和 NincChr6G00302180.1），鉴定了保守和可变位点。结合进化关系，我们进一步推测了氨基酸变异的详细过程（图 6B 和图 S27），这可能有助于解释酶活性分化。总的来说，识别出了 18 个变异位点，可能对 15 个 羌活 CYP736 的新功能化和亚功能化至关重要。在这 18 个氨基酸位点中，有 13 个（L47、H113、K210、V233、A250、K294、V308、S345、A348、Q399、V405、R416 和 F417）可能与 DC 活性相关，11 个（H36、I47、K113、H210、V246、H294、M308、V417、Y426、D442 和 G513）与 OC 活性相关，而 12 个（P36、F47、D113、A233、F250、R345、T348、H399、M405、N416、F426 和 Q513）与 OD 双功能活性相关。

总之，所有的复制事件和序列突变均特异于伞形科，其中一些可能仅限于 羌活。这些复制事件的有限系统性等级支持了复杂香豆素生物合成途径在被子植物中多次独立起源的推测，表明植物生物合成产物的多样化进化。

4. 讨论

香豆素广泛分布于植物界，至少有 75 个科的数千种植物含有香豆素。然而，这通常是指简单香豆素。复杂香豆素，如呋喃香豆素和吡喃香豆素，主要限于伞形科、芸香科、桑科和豆科。简单香豆素与复杂香豆素的主要区别在于异戊烯化和环化步骤（图 1），这些步骤确立了复杂香豆素的基本骨架。此前鉴定的所有 PTs 和环化酶均限于伞形科、芸香科和桑科。参与简单香豆素生物合成的基因已在十字花科、旋花科、芸香科和伞形科中被鉴定。因此，香豆素途径中 PT 和环化酶的功能进化对复杂香豆素的形成至关重要。

2-OGD 家族蛋白负责通过从苯丙素途径向香豆素生物合成途径引导代谢通量来形成简单香豆素。该家族主要包括两种功能不同的酶，C2′H 和 F6′H。C2′H 具有广泛的底物识别能力。除了产生伞形酮外，它还可以催化阿魏酰辅酶 A 形成东莨菪素，而 F6′H 只能接受阿魏酰辅酶 A 作为底物。由于以下原因：(1) 鉴定出的 C2′H 主要来自伞形科和芸香科，而 F6′H 主要来自十字花科和旋花科；(2) 由 C2′H 产生的伞形酮是复杂香豆素生物合成的关键中间体；(3) 植物中的简单香豆素主要以东莨菪素及其衍生物为代表；(4) 在富含复杂香豆素的植物中，简单香豆素的形成由 C2′H 催化，而由 F6′H 催化的植物只能积累简单香豆素，因此我们推测，C2′H 基因的有无也是生产简单香豆素或复杂香豆素的一个优先指标。总之，只拥有 F6′H 的物种只能产生简单香豆素，而拥有 C2′H 的物种则具有生产复杂香豆素的能力。然而，由于目前仅从有限的物种中鉴定出了少量 C2′H 和 F6′H 基因，因此还需要进一步的证据。此外，C2′H 和 F6′H 在蛋白质同一性上没有显著差异，其功能通常可以互换（参见图 S18）。

另一个值得注意的现象是香豆素途径在被子植物中多次起源。不同物种中 PT 的鉴定，以及参与复杂香豆素羟化的基因，表明平行或趋同进化可能有助于香豆素途径的多次起源。此外，尽管在伞形科和桑科中已鉴定出许多 PT 和羟化酶，但在芸香科或豆科中尚未鉴定出 C-PTs，尽管这些科富含复杂香豆素。这些结果表明，芸香科和豆科中的 PTs 可能具有不同于伞形科和桑科的祖先。在本研究中，我们发现来自 羌活（伞形科）的环化酶属于 CYP736 家族，而来自桑科的环化酶属于 CYP76 家族，这支持了它们独立起源的观点（参见图 S25）。此外，参与代表性呋喃香豆素植物蛇床子素羟化步骤的基因在伞形科中为 CYP71AZ1/6，在芸香科中为 CYP82D64，在十字花科中为 CYP82C2/4。因此，我们推测芸香科和豆科中的环化酶可能不属于 CYP76 或 CYP736 家族。我们的推测可能对试图探索芸香科和豆科中复杂香豆素生物合成途径的研究人员提供有限的帮助。正如在 PT 的平行进化中所观察到的那样，环化酶在关系较远的类群中独立出现并不令人惊讶。

香豆素的代谢谱在被子植物中显著不同，不仅在丰度上不同，而且在结构特征上也不同。这些结构差异包括线性和角状配置以及呋喃型和吡喃型结构的变化。一个流行的假设认为，角状香豆素起源于线性香豆素。该假设主要基于以下观察：(a) 没有植物物种单独存在角状香豆素，而是总是与线性香豆素共存；(b) 在线性香豆素中发现的家族更多，包括伞形科、芸香科、桑科和豆科，而角状香豆素仅限于伞形科的某些物种；(c) 角状香豆素可以与线性香豆素协同作用，从而进一步提高香豆素对昆虫的毒性。然而，迄今为止，尚无直接证据支持该假设。PT 的功能进化可能为该假设提供支持，但在伞形科中可用于研究这一主题的案例（PsPT2 和 PpPT2）太少。此外，某些线性 PT 表现出极弱的角状活性，这与 (a) 一致。关于呋喃香豆素和吡喃香豆素，由于极少数环化酶被鉴定出来，因此其结构多样性的机制知之甚少。尽管根据我们的分析，不同的环化酶可能源于不同的基因复制方式（串联、异位和 WGD/片段），最终产生不同的香豆素结构（图 6），控制底物和/或产物特异性的关键氨基酸位点突变可能对于结构多样性更为重要。根据我们最近的研究，反应溶液的酸碱性显著影响呋喃香豆素和吡喃香豆素的形成。然而，确切的机制需要在蛋白质水平上进一步研究。