Adv. Sci. 中国农业科学院-全基因组尺度代谢模型指导的大肠杆菌工程化实现菊酸的从头生物合成

菊酸是天然杀虫剂除虫菊酯中一种非常规的单萜结构单元，具有显著的抗昆虫活性，使其工业生物合成成为可持续农业的一条前景广阔的途径。本研究以基因组尺度代谢模型为指导，设计并构建了一株用于高效生产菊酸的大肠杆菌细胞工厂。在重建生物合成通路的基础上，模拟了外源模块引入所引起的代谢变化。通过该模型，识别出了由ispA催化的关键代谢分支节点，并利用合成小RNA抑制其表达，成功重定向了代谢流，使前体菊醇和菊酸的产量分别提高了162%和59%。GEM模型还预测了下游脱氢酶表达水平对菊酸产量的影响，通过进一步优化脱氢酶基因的拷贝数，实验中将菊酸产量显著提升了570%。通过将代谢支路阻断策略与拷贝数优化相结合，在生物反应器中实现了141.78 mg L⁻¹的菊酸产量，创下新高。这项工作将计算机模拟优化与湿实验操作无缝衔接，通过代谢网络工程显著提升了目标代谢物产量，为构建高效天然生物产物细胞工厂提供了新路径。

除虫菊酯是一种源自除虫菊（Tanacetum cinerariifolium）的天然杀虫剂，但由于植物提取成本高昂、化学合成工艺复杂，其发展面临限制。尽管在传统植物育种和农业优化方面付出了诸多努力，这些方法仍未能使除虫菊酯成为一种经济实惠的农用杀虫剂。菊酸作为除虫菊酯的酸组分之一，因其抗昆虫活性，同样被视为一种天然杀虫剂。此外，它还是合成丙烯菊酯、溴氰菊酯等拟除虫菊酯类杀虫剂的关键原料。然而，其复杂的立体化学结构使得化学合成极具挑战性，导致副产物众多且纯化步骤繁琐。

菊酸的生物合成途径已在除虫菊中得到阐明。其分子结构具有独特的单萜骨架，由两个二甲基烯丙基二磷酸（DMAPP）单元以头-中对位缩合方式形成。菊基二磷酸合酶利用质体甲基赤藓糖醇磷酸（MEP）途径产生的两个DMAPP分子，催化缩合生成菊基二磷酸。随后，包括Nudix1在内的磷酸酶将CDP水解为菊醇。接着，ADH2氧化菊醇生成菊醛，再由ALDH1进一步氧化最终生成菊酸（图1A）。

已有数项利用植物底盘进行菊酸生物工程的研究，证明了其异源生产的潜力。徐等人成功在富含异戊二烯的番茄果实中重建了该合成途径，为蔬菜害虫防治提出了新思路，其产量记录约为1.09微摩尔/克干重。同样，胡等人在菊花中过表达CDS，检测到挥发性菊醇的产量约为50皮摩尔/小时/克鲜重及1.1毫摩尔浓度，实现了对蚜虫的双重防护。尽管取得了这些进展，值得注意的是，目前尚无研究报道利用微生物底盘生产菊酸。

众多成功案例已证明大肠杆菌在生产天然代谢物（包括萜类化合物的生物合成）方面的有效性。Martin等人优化了青蒿素前体的异源合成，产量达到24微克/升；他们所构建的大肠杆菌MVA代谢途径为其他萜类化合物的合成奠定了基础。Ajikumar等人将紫杉二烯代谢途径在大肠杆菌中进行模块化设计，实现了约1克/升的产量。Park等人报道通过异源表达在大肠杆菌中生产了432毫克/升的虾青素。Sun等人通过基因工程手段，获得了高达3520毫克/升的番茄红素产量，成果令人瞩目。

构建微生物细胞工厂的关键在于设计可控且高效的生物合成系统，这需要在平衡物质与能量代谢的前提下，对细胞模块间的协同作用以及模块对微生物宿主的适应性进行迭代优化。基因组尺度代谢模型已成为系统水平代谢工程设计中的重要工具，能对代谢网络进行全面描述。该模型网络整合了生物化学与基因组学知识，包括代谢反应的精确化学计量学、代谢物的化学式与电荷，以及基因、蛋白质与反应之间的关联。GEM能够在系统层面预测细胞表型，从而指导遗传操作以促进生长并最大化目标化合物的生产效率。

在本研究中，通过一种开源、经济高效的策略，提升了大肠杆菌中菊醇和菊酸的产量。这包括引入DMAPP供给模块，并整合菊醇与菊酸合成模块。借助转录组学与代谢组学数据，构建了一个基因组尺度代谢模型来模拟改造菌株。利用机器学习工具UniKP，筛选了DMAPP相关的支路反应，并结合代谢网络建模，鉴定出ispA为关键分支节点，随后通过实验验证了这一点。在模型指导下，确定了菊酸合成中的限速步骤，并通过优化关键基因拷贝数显著提高了产量。经过包括模块整合在内的迭代优化，在摇瓶培养中实现了菊醇产量超过117.0 mg L⁻¹、菊酸产量达15.3 mg L⁻¹。放大至2升发酵罐后，菊酸产量进一步提升至141.78 mg L⁻¹，这为未来基于除虫菊酯的天然杀虫剂大规模发酵生产展现了可行性。

为构建稳健的菊酸生物合成底盘，设计了三个模块（图1B）。首先，为确保菊酸关键前体DMAPP的持续供应，在大肠杆菌BL21(DE3)中通过共表达预先验证的pMevT和pMBI质粒，引入了外源甲羟戊酸途径，构建了菌株CS3。随后，通过引入负责菊醇生物合成的第二模块——包含来自除虫菊的支架形成基因CDS和Nudix1基因——构建了菌株CS5-1。最终，为了生产菊酸，在CS5-1中通过引入来自除虫菊的两种脱氢酶ADH2和ALDH1整合了第三模块，从而构建了菌株CS6。

为进一步提升菊酸产量，在对工程菌株进行代谢建模的基础上，执行了一系列综合性的优化策略。该模型整合了多组学数据，能够精确定位关键代谢通量节点，并促进代谢流的理性重定向。具体优化措施包括：利用合成sRNA抑制ispA基因（菌株CS5-3和CS11），以及精细调控ADH和ALDH基因的拷贝数（菌株CS7-CS10）。这些优化策略的详细阐述将在后续章节中提供（图1B）。

通过共表达来自酿酒酵母和大肠杆菌的7个相关基因，重建了甲羟戊酸途径（图2A）。经GC-MS分析，与野生型菌株相比，工程菌株CS3表现出高浓度的多种单萜和倍半萜（如香叶醇和法尼醇）积累（图2B）。该结果虽证实了萜类前体的充足供应，但这些挥发性萜类的大量积累凸显了对该底盘进行代谢网络优化的迫切需求。

为模拟异源甲羟戊酸途径引起的代谢扰动下的代谢通量分布，收集并分析了转录组与代谢组数据。从五个菌株（野生型、CS3、CS4、CS5-1和CS6）在三个时间点（诱导前、诱导后16小时、诱导后48小时）采集了包含45个样本的综合数据集，每个菌株设置三个生物学重复，进行RNA测序。获得了不同菌株间的基因表达谱。对这些表达谱的主成分分析显示，在诱导中期，野生型与实验菌株之间存在明显分离（支持信息图S3）。此外，差异基因表达模式以及基因本体和KEGG富集分析表明，与野生型相比，实验菌株的碳代谢与次级代谢显著上调。

同时，通过GC-MS测定了相应菌株的挥发性代谢组。如预期所示，实验菌株中的大多数代谢物与对照组相似，变化主要局限于部分萜类代谢途径，对初级代谢影响甚微。代谢组数据的主成分分析进一步支持了这一观察。差异代谢物分析显示，与对照组相比，外源萜类代谢途径的引入显著影响了实验菌株中少数化合物的生物合成。同时，正交偏最小二乘判别分析表明，通过引入的萜类途径合成的化合物（如香叶醇和法尼醇）显著影响了OPLS-DA图中菌株的分离。

随后，获取的多组学数据被用于构建代谢模型。添加到原始iEC1356_Bl21DE3模型中的代谢物和反应的详细信息在补充材料中提供。使用标准化基因组尺度代谢模型质量检查工具MEMOTE，对模型BL21_WT、BL21_CS3和BL21_CS5-1的质量进行了严格评估。生成的报告证实了这些模型的化学计量一致性和质量平衡。接着，采用通量变异性分析²³和通量采样方法，研究了野生型菌株与引入甲羟戊酸途径的CS3菌株之间代谢通路分布的差异（图2C，D）。

引入甲羟戊酸途径不仅影响了萜类生物合成通路，还改变了大肠杆菌中超过一半的代谢网络。随后，通量采样结果预测，CS3菌株中的甲羟戊酸途径可能将本应流向菊醇合成路径的通量分流至副产物（C10、C15萜类）的生产，包括香叶醇和法尼醇（图2E）。这一预测与离子色谱检测到的单萜和倍半萜含量增加的结果一致（图2B）。这些建模结果共同表明，阻断竞争性途径的分流通量可将副产物的生物合成重定向至菊醇生产，从而提高菊醇产量。

菊醇的生物合成以DMAPP为底物，由CDS和Nudix1催化其转化（图3A）。在CS3菌株（菌株CS4）中仅表达CDS，经48小时培养后可产生33.31 mg L⁻¹菊醇（图3B）。菊醇的提取离子流色谱图在质荷比123.0处生成。虽然分子离子峰（m/z 154.2）存在，但碎片离子峰m/z 123.0的响应强度显著更高，这提高了检测灵敏度，使化合物检测更为可靠（图3B）。该菌株中菊基二磷酸水解为菊醇的过程，可能是由于野生型大肠杆菌中磷酸水解酶的广谱性所致。为提升菊醇产量，我们在CS4菌株中进一步引入了除虫菊来源的Nudix1磷酸水解酶（菌株CS5-1）。CS5-1经48小时发酵后，菊醇产量达到44.69 mg L⁻¹，是CS4菌株产量的134.16%。

为确定导致DMAPP通量偏离菊醇生产的关键反应，模拟了DMAPP代谢分支点的通量分布。首先，利用转录数据并结合基于机器学习的周转率预测算法UniKP²¹工具，计算了前体代谢物DMAPP在其合成途径中参与各分支反应的最大反应速率。结果显示，在所有已识别的分支反应中，ispA催化的DMAPP利用速率最高（图3C）。随后，在不同底物浓度梯度下，计算了与底物DMAPP相互作用的各酶的通量比例。计算分析表明，在不同底物浓度水平下，ispA对DMAPP的分流比例均显著高于其他参与旁路生物合成途径的酶。该假说认为，在大肠杆菌中催化DMAPP和IPP合成GPP与FPP的法尼基二磷酸合酶（ispA），显著分流了菌株CS5-1中的关键前体代谢物DMAPP。

因此，提出了通过抑制ispA表达的分流阻断策略来提高菊醇产量。然而，与先前文献报道一致，对ispA敲除的模型模拟显示生物量为零，这表明ispA是大肠杆菌生长所必需的，直接敲除会导致菌株死亡。利用通量平衡分析在BL21_CS5-1模型中模拟了ispA基因抑制的效果（图3D，E）。模拟前，对CS5-1模型设置了特定约束：葡萄糖被指定为唯一碳源，摄取速率限制为20 mmol/g·DW/h。此外，由于ispA对细胞生长至关重要，MEP和MVA途径的总通量，以及GPP和FPP通量，最初均根据与细胞生长和菊醇产量相关的实验数据进行了约束。

为了模拟ispA表达下调对细胞生长和菊醇生产的可能影响，部分限制了ispA通量。预测结果表明，当ispA表达降至CS5-1的1/16同时维持细胞活力时，菊醇总通量将增加2.64倍（图3D，E）。因此，该模型模拟表明，抑制ispA表达可显著提高菊醇产量并避免致死表型，从而证实了通过“分流阻断”策略将通量重定向至目标产物的有效性。

为通过实验验证模型的预测，在菌株CS5-1中利用大肠杆菌异源RNA干扰系统应用合成小RNA抑制了ispA表达，构建了菌株CS5-3。定量PCR验证了该抑制的成功，显示CS5-3中ispA表达水平比CS5-1对照组降低了87.77%（图3F）。同时，对两菌株中的菊醇含量进行了定量分析。培养48小时后，CS5-3表现出显著更高的菊醇浓度（117.03 mg L⁻¹），相比菌株CS5-1的53.24 mg L⁻¹增加了约119.83%（图3F）。该结果与GEM的预测完全吻合，表明菊醇产量随ispA表达降低而增加。此外，与菌株CS5-1相比，菌株CS5-3的香叶醇（降低60.94%）和法尼醇（降低97.32%）浓度均显著下降（图3G）。这些发现最终证实，抑制ispA表达是显著提升菊醇产量的一种可行且关键的策略。

为高效探索多种优化策略并缩短总体实验周期，在菌株CS5-1中开展了一系列平行试验，以优化菊醇生产条件（包括培养时间、温度、诱导剂浓度和培养基），同时进行代谢通量优化。通过间歇取样监测了菌株长达102小时的生长与产物曲线。结果显示，细胞密度在26小时达到峰值，但菊醇在培养液中仍持续积累。进一步试验表明，在2×YT培养基中，最佳IPTG浓度和温度分别为1.0 mM和20°C。结合上述优化条件，CS5-1菌株在54小时发酵后，于特定条件下实现了80.6 mg L⁻¹的菊醇产量。进一步分析显示菊醇在胞内与胞外的分布情况，约95%存在于培养基中。这一结果可能归因于菊醇的亲脂性和低分子量特性，使其易于自由扩散穿过细胞膜。

在CS5-1的基础上，尝试使用大肠杆菌偏好密码子优化CDS和Nudix1的cDNA序列（菌株CS5-2）。遗憾的是，优化后的CS5-2菌株获得的菊醇水平反而更低。这一结果与先前研究一致，即改变密码子使用模式并非总能提高目标产物产量，并突显了基因优化的复杂性。优化基因的过表达可能导致包涵体形成，这通常与错误折叠或聚集的蛋白质有关。

随后的工作集中于菊酸生物合成模块（模块III），该模块包含从菊醇出发的两个氧化步骤（图4A）。通过在菌株CS5-1中引入脱氢酶ADH2和ALDH1（通过pRSFDuet-ADH2-ALDH质粒）进行初步工程改造，使得菌株CS6能够检测到由菊醇转化生成的0.63 mg L⁻¹菊酸（图4B，F）。然而，仍有高达49.30 mg L⁻¹的菊醇未被转化，凸显了优化ADH2和ALDH1表达以提高底物转化率和产量的必要性（图4F）。

为确定菊酸生物合成中的瓶颈反应，我们利用CS6的基因组尺度代谢模型和米氏方程比较了ADH2与ALDH1的动态反应速率。首先建立了代谢网络建模工作流程图（图4C）。第一步，通过整合文献报道的酶动力学参数与实验测得的基因表达水平，计算了两种酶的初始最大反应速率，并以此约束其在代谢网络模型中的最大反应速率比。第二步，为识别瓶颈反应，模拟了反应表达增强策略（如结合基因拷贝数增加n倍）以预测产物产量。计算结果显示，ALDH1的最大反应速率低于ADH2（图4D）。随后，模型模拟了在不同过表达策略下的理论最大通量，提供了关键见解：单独提高ALDH1通量能显著提升菊酸产量至约3.3 mmol/g·DW/h（绿色曲线），而单独增加ADH2通量则无明显影响（橙色曲线）。关键在于，模型预测同时优化ALDH1和ADH2会产生协同放大效应，将产量提升至约7 mmol/g·DW/h（蓝色曲线）（图4E）。这强调了协同调控这两种酶以克服已识别的瓶颈并大幅提高菊酸产量的必要性。

依据计算预测，通过实验评估了ADH2基因拷贝数对菊酸产量的影响。对构建的含有一、二、三个ADH2拷贝的工程菌株进行了表征。双拷贝菌株（CS7）的ADH2表达量约为单拷贝菌株的5.4倍。虽然这导致了菊醇的显著减少，表明其消耗增加，但菊酸的总产量并未大幅提升。同时，检测到了中间产物菊醛的存在，这突显了后续转化步骤中持续存在的瓶颈（图4B，F）。在更高ADH2拷贝数下菊醛产量的下降，可能源于ADH2表达量的降低，正如菌株CS8与CS7相比所示。推测这是由于转录-翻译偶联所致，即过量的基因拷贝导致mRNA上的"交通拥堵"，从而降低了整体蛋白质合成。这一实验观察与计算模型的预测相符——单独增加ADH2通量影响有限，从而强调了协同优化ALDH1和ADH2的必要性。

基于这些综合模拟与初步实验结果，ALDH1被确定为菊酸合成中的主要限速酶。为实验验证该假设并进一步提高菊醇向菊酸的转化，随后在CS7菌株的骨架基础上优化了ALDH1的拷贝数。对具有不同ALDH1拷贝数的菌株（CS9为两个拷贝，CS10为三个拷贝，对照CS7为一个拷贝）的评估表明，CS10的ALDH1表达水平显著更高（比CS7增加120%）。我们的qPCR结果显示，虽然在16小时时间点菌株CS9的ALDH1表达低于CS7，但两者ALDH1的时间表达谱存在差异。在诱导后4、8、12和16小时的测量显示，CS9中ALDH1的表达峰值出现时间早于CS7。这种更早的表达可能是CS9最终菊酸产量更高的关键因素，尽管其最终时间点的表达水平低于CS7。同时，在CS10中引入三个ALDH1基因拷贝显著改善了菊醇转化，在16小时内实现了6.48 mg L⁻¹的菊酸产量，是单拷贝ALDH1菌株的5.7倍（图4B，F）。该实验结果精准验证了模型预测，明确证实了ALDH1是该途径的关键瓶颈。

为进一步提升菊酸产量，实施了整合上游分流阻断与下游强化转化的协同策略。通过将ispA沉默系统（先前在CS5-1中优化）引入拷贝数优化的CS10菌株，构建了菌株CS11（图5A）。在已建立的菊酸生产最适培养条件下，菌株CS11在48小时内实现了15.30 mg L⁻¹的产量，而CS10的产量为9.59 mg L⁻¹（图5B）。

为评估可放大性，随后在2升生物反应器中进行了初步实验，采用正十四烷进行产物收集。通过联动控制搅拌转速（400-1000 rpm范围）将溶解氧维持在30%。通过自动添加氨水溶液将pH维持在7.0。采用补料分批策略，使用0.5 g mL⁻¹葡萄糖溶液，在补料阶段将葡萄糖浓度维持在3-5 g L⁻¹。该工艺最终在168小时实现了141.78 mg L⁻¹的菊酸产量（图5C），这是目前达到的最高生产水平。

随着现代农业产出的增长，农药需求也日益上升。随着合成生物学的兴起，对环境保护和可持续发展的需求日益增长。随着化学合成的拟除虫菊酯类农药（如苯锈啶和溴氰菊酯）的广泛使用，昆虫对其抗性迅速增强，而天然除虫菊酯依然保持高效杀虫活性。

菊酸作为除虫菊酯的关键组分，具有天然的驱虫特性。尽管除虫菊酯的代谢途径尚未完全明晰，但构建菊酸的微生物细胞工厂可作为生产除虫菊酯的基础，或作为独立的杀虫剂，具有广阔的生态效益前景。本研究致力于利用已阐明的除虫菊酯合成相关基因进行微生物工程改造，旨在降低除虫菊酯的生产成本。我们已成功在大肠杆菌中合成了除虫菊酯的酸基组分——菊酸，其产量高于菊花和番茄中的菊酸产量。

最优化的代谢工程需要对设计的代谢途径和调控机制有系统性的理解。在本研究中，通过引入菊醇合成基因，实现了菊醇在大肠杆菌中的生产。通过整合代谢组学与转录组学分析，得出结论：DMAPP的分流限制了该宿主中除虫菊酯的产量。为克服这一限制，结合转录组学和代谢组学数据构建了高分辨率代谢模型，并预测了DMAPP的分支途径。通过比较使用合成小RNA抑制ispA表达的大肠杆菌与未抑制ispA表达的大肠杆菌中菊醇和菊酸的产量，结果显示合成的菊醇和菊酸分别增加了162%和59%。

利用基因组尺度代谢模型，预测了菊酸合成的关键限速步骤。模型揭示了一个意外现象：与单独上调相比，协同增强关键限速酶与非限速酶能显著提高产量上限。实验验证证实了这一预测的协同效应，凸显了该模型在指导高效代谢途径工程方面的潜力。具体而言，增加ADH和ALDH的拷贝数，结合代谢模型预测与基因沉默技术，使得48小时后菊酸积累量提高了23.29倍。

整合了基因组、生物化学及代谢网络数据的代谢模型，被证明在预测基因敲除影响和识别过表达靶点方面具有重要价值，即便在数据有限的情况下也展现了其实用性。然而，这些模型存在局限性，包括其静态特性——这忽略了动态调控信息，未能考虑转录调控的动态性和酶活性变化（如变构效应），可能导致预测偏差，尤其对于瞬时代谢反应。另一局限性在于某些尚未充分研究的酶缺乏全面的动力学参数。在未来的工作中，计划通过整合多组学数据及人工智能模型生成的基因调控模式来更好地约束模型，从而提高其准确性。此外，还将考虑整合热力学约束以进一步完善模型精度。

总之，研究利用了模型指导的代谢工程方法及其提供的高精度细胞表型预测，显著降低了开发高产微生物菌株相关的实验成本和时间。该方法有助于设计更合理、更高效的代谢工程策略。除这些技术外，还采用了传统代谢工程方法，如培养基优化、温度控制和诱导强度优化，以提高菊酸产量。实现了141.78 mg L⁻¹的菊酸产量，这是目前报道的最高产量。在生物反应器发酵过程中，细菌培养物的光密度低于预期，且在诱导后观察到细胞密度缓慢下降。推测这可能由三个因素导致。首先，菌株CS11中的五个质粒及其相应的抗生素抗性标记带来了显著的代谢负担，减缓了生长速率。其次，用于重定向代谢流的ispA基因抑制策略已知会在诱导剂激活后抑制细胞生长。最后，发酵后期菊酸的积累可能产生了细胞毒性，同样导致了细胞密度下降。解决这些问题将在未来进一步提高菊酸产量。本工作为天然产物的生物合成提供了宝贵的参考与启示。

本研究通过引入除虫菊植物的催化基因，构建了用于菊酸生物合成的大肠杆菌底盘。为提高产率，基因组尺度代谢模型根据时序转录组与代谢组数据模拟了代谢扰动。如模型预测所示，识别并部分阻断了关键旁路途径，使得产物菊醇和菊酸的产量分别提高了162%和59%。为提升下游转化效率，模型指导识别了限速步骤，通过优化关键基因拷贝数使产量提高了470%。在生物反应器中，菊酸产量达到创纪录的141.78 mg L⁻¹。限制大肠杆菌中菊醇与菊酸产量的两个关键因素是支路代谢的分流和下游转化速率。