Gap-free genome assembly and comparative analysis reveal the evolution and anthocyanin accumulation mechanism of Rhodomyrtus tomentosa
无缺口基因组组装及比较分析揭示了桃金娘的进化和花青素积累机制
摘要
桃金娘（Rhodomyrtus tomentosa）是桃金娘科的重要肉质果树和著名的药用植物，广泛种植于世界热带和亚热带地区。然而，由于缺乏参考基因组，桃金娘的进化和基因组育种研究受到阻碍。在此，我们使用PacBio和ONT长读长测序技术，提出了桃金娘的染色体级无缺口T2T基因组组装。我们组装了大小为470.35 Mb的基因组，contig N50约为43.80 Mb，共有11条拟染色体。在该基因组中注释了33,382个基因和239.31 Mb的重复序列。系统发生分析阐明了桃金娘从14.37百万年前（MYA）开始的独立进化，并与其他桃金娘科物种共享一个最近的全基因组复制（WGD）事件。我们鉴定了桃金娘中四种主要的花青素化合物及其合成途径。比较基因组和基因表达分析表明，桃金娘的着色和高花青素积累倾向于由花青素合成途径的激活决定。MYB转录因子的正向选择和上调是这一过程中隐含的因素。在桃金娘中还检测到与花青素运输相关的OMT和GST基因的拷贝数增加。表达分析和途径鉴定丰富了在桃金娘科肉质果实发育过程中淀粉降解、对刺激的响应、激素效应和细胞壁代谢的重要性。我们的基因组组装为研究桃金娘科物种的起源和分化提供了基础，并加速了桃金娘的遗传改良。
介绍
桃金娘（Rhodomyrtus tomentosa (Ait.) Hassk），中文名为“桃金娘”，也被称为玫瑰桃金娘，是桃金娘科的一种开花植物，原产于南亚和东南亚，包括中国、印度、菲律宾、马来西亚和苏拉威西岛[1]。在传统中医药（TCM）中，桃金娘被视为一种具有滋补血液系统、抗风湿、治疗咯血、腹泻和子宫出血的药用草药[2]。现代药理学研究已鉴定出从桃金娘中提取的代谢物具有多种药理作用，如抗菌、抗肿瘤、抗炎和抗氧化[3]。
果实着色和花青素积累是桃金娘在桃金娘科中的显著特征。基于其独特的结构特性，花青素具有多种生物活性，包括抗氧化、抗癌、抗炎、抗动脉粥样硬化、抗高血压和抗菌活性[2]。在过去的几十年里，花青素的生物机制研究引起了全球学术界的兴趣。先前的研究初步确定了桃金娘中花青素的合成类型及其药理作用[3]。然而，关于桃金娘中花青素和其他生物活性物质的遗传基础的研究很少。
分子机制分析对于了解花青素和其他生物活性化合物的形成至关重要。许多研究表明，花青素合成是黄酮类化合物合成的一部分[4]。该途径受各种转录因子调控，主要是骨髓瘤（MYB）转录因子[5]。花青素的糖基化和运输在果实着色中也起着重要作用[6, 7]。然而，不同桃金娘科物种之间的变异导致桃金娘独特的紫色果实特征仍然未知。
根据果实类型的传统基础，分类学研究者将桃金娘科物种分为两个亚科：桃金娘亚科（肉质果实）和白千层亚科（蒴果）[8]。具有肉质果实的桃金娘被归类为桃金娘亚科。对另一个桃金娘亚科物种番石榴（Psidium guajava）的先前研究表明，其果实软化由淀粉降解和细胞壁活性促进[9]。这种趋势是否在桃金娘亚科物种中普遍存在，显然需要更多的研究来解决。
此外，一般观点认为桃金娘的食用和药用部位主要集中在果实上。近年来，桃金娘叶和茎中的有效药用成分，如桃金娘酮B和桃金娘酮C，逐渐引起关注[10, 11]。然而，由于缺乏基因组信息，桃金娘中药用类黄酮的全面和系统研究仍然缺乏。
无缺口基因组提供了尽可能完整地识别基因组信息的机会，被视为基因组组装的最终目标。先前的研究已经在包括水稻[12, 13]、拟南芥[14]和西瓜[15]在内的多种物种中发布了无缺口基因组，但对于桃金娘科物种，尚未报告无缺口基因组。
在这项研究中，我们报告了优质桃金娘自交系“LFSTJN-1”的端粒到端粒（T2T）无缺口基因组，结合了PacBio HIFI测序、Oxford Nanopore技术（ONT测序）以及Hi-C技术。我们确定了桃金娘中花青素的主要化合物及其合成途径。比较基因组提供了花青素积累的可能调控机制。表达趋势分析和途径鉴定进一步丰富了淀粉降解、信号传导和细胞壁活性在桃金娘科肉质果实形成中的重要性。我们的基因组组装为研究桃金娘科的起源和分化提供了基础，并加速了桃金娘的遗传改良。
结果
桃金娘的T2T无缺口参考基因组
选择广泛种植于中国南部的桃金娘品种LFSTJN-1进行T2T无缺口参考基因组组装（图1a）。基因组调查表明，桃金娘的基因组大小约为450.77 Mb，杂合率为0.29%（图S1，见在线补充材料）。使用多种测序平台开发高质量的桃金娘基因组组装。我们使用PacBio sequel II平台生成了约33.40 Gb（约74.08×覆盖率）的HiFi读数，并使用Oxford Nanopore技术（ONT）生成了约11.02 Gb（约24.44×覆盖率）的超长读数（表S1，见在线补充材料）。HiFi读数和ONT读数的N50长度分别大于16 kb和100 kb（图S2，见在线补充材料）。初步组装采用Hifiasm对HiFi读数进行处理，生成了N50长度为29.87 Mb的482.47 Mb contigs（表1）。通过NextDenovo使用ONT数据生成了包含19个contigs、大小为471.22 Mb、N50长度为39.17 Mb的组装，并通过NextPolish使用HiFi读数及Illumina平台的短读数进行校准。这两个组装使用3D-DNA与Hi-C数据处理以纠正组装错误并删除冗余序列。使用基于HiFi的组装填补了ONT组装的缺口。在填补所有剩余缺口后，生成了470.35 Mb的无缺口桃金娘基因组，包含11条染色体，contig N50长度为43.80 Mb。染色体数目与先前的记录相对应（https://ccdb.tau.ac.il/）。最后，使用七碱基端粒重复序列（‘CCCTAAA’）作为序列查询，我们鉴定了所有22个端粒，并构建了桃金娘基因组的11个T2T拟染色体（图1C；表S2，见在线补充材料）。
展示桃金娘的图像、基因组组装及基因组特征。A 桃金娘的图像，包括花和叶。B 基于组装的HI-C相互作用矩阵。C 桃金娘基因组中的基因密度、端粒和着丝粒的分布。D 桃金娘的基因组特征，从圈的内侧到外侧依次为基因组共线性位点（基因组序列同源块）的分布、GC含量、长末端重复（LTR）Copia的密度、LTR Gypsy的密度、总LTR的密度以及总转座子（TE）的密度。
Table 1
Statistics of genome assembly among Myrtaceae species
| Statistics without reference | R. tomentosa (HIFI-contig) | R. tomentosa (ONT-contig) | R. tomentosa (Gap-free) | P. guajava | E. grandis |
|---|---|---|---|---|---|
| Number of contigs | 269 | 19 | 11 | 44 | 4943 |
| Number of contigs (≥ 50 000 bp) | 113 | 19 | 11 | 35 | 288 |
| Largest contig | 43 631 835 | 49 958 455 | 50 111 023 | 50 577 630 | 83 952 244 |
| Total length | 482 474 435 | 471 221 806 | 470 350 250 | 443 755 635 | 691 346 258 |
| Total length (≥ 50 000 bp) | 477 256 938 | 471 221 806 | 470 350 250 | 443 448 379 | 651 047 186 |
| N50 | 29 878 240 | 39 171 969 | 43 802 537 | 40 370 300 | 57 472 304 |
| L50 | 7 | 6 | 5 | 5 | 5 |
| GC (%) | 40.63 | 40.56 | 40.56 | 39.48 | 39.29 |
| # N’s (Mismatches) | 0 | 0 | 0 | 2900 | 50 906 790 |
| # N’s per 100 kbp (Mismatches) | 0 | 0 | 0 | 0.65 | 7363.43 |
进一步使用TRF分析，在桃金娘基因组中鉴定出两种主要候选着丝粒串联重复序列，单体长度分别为422 bp和153 bp（图S3，见在线补充材料）。在11条拟染色体中的每一条上都鉴定出一个假定的着丝粒，长度范围从0.35 Mb到3.49 Mb。桃金娘中的基因密度从着丝粒向染色体末端增加，这与多种物种中的分布一致（图1C；表S3，见在线补充材料）。
桃金娘基因组组装的质量评估
为了确认桃金娘基因组组装的遗传背景，实施了多种数据和方法。从序列的角度来看，从生成的HI-C短读序列库中得到的相互作用矩阵表明，11条染色体已完全且合理地组装（图1B），其数量与先前的记录一致。同时，调查和组装中使用的短读和HIFI读数数据的比对显示了99.96%和99.93%的映射率。从遗传角度来看，基于eudicots_odb10和embryophyta_odb10的BUSCO评估表明，分别在组装中完全找到了97.7%（2284个中的2067个）和99.0%的核心保守植物基因（表S4）。此外，我们成功地将不同组织和发育阶段生成的RNA-Seq数据集的92.4%至96.5%映射到组装中。长末端重复（LTR）的完整性测试显示，组装的LTR组装指数（LAI）值为16.16，这与无缺口组装如猕猴桃（Actinidia chinensis cv. 'Hongyang'）相似，高于番石榴（~9.38）[9, 16]。这些数据表明，桃金娘基因组组装质量很高。
基因组注释和重复序列识别
基于de novo和同源证据的组合表明，重复序列占组装基因组的50.92%（239.51 Mb），其中190.23 Mb（79.42%）属于长末端重复逆转录转座子（图1d；表S5，见在线补充材料）。利用重复掩盖的基因组，通过结合de novo、同源和转录组预测结果来预测基因结构。最终，在我们组装的基因组中预测了33,382个编码蛋白的基因，平均编码序列大小为1483.1 bp，平均每个基因有4.92个外显子。通过BUSCO评估注释的完整性。在桃金娘的基因集中检测到约93.2%的完整基因元件，这表明大多数保守基因被精确预测，反映了预测结果的高可靠性（表S6，见在线补充材料）。
基因组进化
系统发生和进化事件分析对于澄清桃金娘在桃金娘科中的地位是必要的。基于617个单拷贝直系同源基因的时间尺度系统发生树表明，桃金娘与石榴的分化时间约为1437万年前（MYA），而在桃金娘中，953个和714个基因家族显示出扩张和收缩（图2A和B）。三个桃金娘科物种以及石榴（Punica granatum）的共线性分析显示了完整且连续的共线性。进一步的染色体变异鉴定揭示，石榴染色体的重组在桃金娘科物种的形成中起了重要作用。同时，与桃金娘的单独比较中，在大桉树（Eucalyptus grandis）中检测到比番石榴（P. guajava）更多的染色体倒位，这可能推动了桃金娘科内部的分化（图2C）。
桃金娘及其他物种的系统发育和共线性分析。A 桃金娘及其他九种植物的系统基因组分析和基因家族分布。节点上的值表示从现在起的分化时间（单位：百万年前，MYA）；红色三角形显示了桃金娘科中的WGD事件。B 十种物种中所有基因家族的基因分布。C 桃金娘科三种物种和石榴的基因组共线性比较。D 物种和石榴WGD事件的同义置换位点（ks）图。
桃金娘科物种包括桃金娘在内的同源基因的同义置换位点（Ks）分布在Ks ≈ 1.25处出现一个峰值，这表明桃金娘与其他桃金娘科物种共享了一个最近的全基因组复制事件（WGD）（图2d）。共线性分析也表明桃金娘保留了与其他桃金娘科物种（包括番石榴）相似的祖先多倍性特征（图S4，见在线补充材料）。
基于Ks的分布，我们进一步确定了桃金娘目物种中共有WGD事件的年龄。基于在石榴和桃金娘目物种之间的直系同源基因检测到的Ks ≈ 1.4的峰值，我们进一步推断出桃金娘目的平均Ks/年速率约为6.55 × 10^−9至9.39 × 10^−9（95%的置信区间）。这些数据表明，桃金娘科的WGD事件发生在约66.58至95.50百万年前（图2a）。
桃金娘器官间的基因表达模式和与果实软化相关的代谢过程
植物的生长和生命活动离不开基因表达。为了探索桃金娘器官中的基因表达模式，使用来自不同器官和发育阶段的十种样品中表达的25 038个基因进行了加权基因相关网络分析（WGCNA）。选择了阈值功率为八，这是适当符合无尺度拓扑指数的最低功率（图S5B，见在线补充材料）。在合并动态分析后，揭示了包括12 042个基因在内的共25个共表达模块，其余的基因被划分为其他模块（灰色）（表S7，图S5c，见在线补充材料）。
先前对番石榴的研究表明，果实软化和成熟过程与细胞壁活性和淀粉降解有关。桃金娘是进一步探索桃金娘科肉质果实发育机制的良好材料。相关性分析检测到器官和发育阶段样品的高度相关的共表达模块（r > 0.8）（图3a）。进一步的基因富集分析表明，与细胞壁活性相关的途径在模块16中得到富集，该模块与F2密切相关（r = 0.96）。激素反应的基因本体（GO）术语在与F3相关的模块（模块24；r = 0.88；图3B）中得到富集（图S6，见在线补充材料）。在果实发育的所有阶段都检测到外部响应和压力的基因富集。激素调控细胞壁降解在包括番石榴在内的多种果实成熟研究中已有报道，该物种属于同一亚科并已有基因组发布[9]。同时，外部响应对肉质果实发育的巨大影响得到了广泛确认[17]。
不同组织中的基因表达模式以及桃金娘的淀粉降解途径。A 热图显示了桃金娘各模块与不同组织和发育阶段之间的相关性。给定模块与样本之间的相关系数由细胞的颜色和细胞内部的文本表示（上方数字是Pearson相关系数的值（r），下方数字是P值）。红色和蓝色分别表示正相关和负相关。B 模块24中基因的GO富集；显示了调整后的P < 0.05的本体论，而本文中提到的本体论以文本形式显示：详细信息见表S12（在线补充材料）。C 桃金娘果实发育期间淀粉含量的变化，x轴表示开花后的时间和相应的果实发育时期（F1–F4）。D 桃金娘中的淀粉降解途径及相应的基因表达。
此外，类似于番石榴，桃金娘的淀粉含量随着果实发育而减少（图3c）。基于京都基因与基因组百科全书（KEGG）注释，我们在桃金娘中鉴定了属于7个家族的20个淀粉降解相关基因。表达分析表明，与淀粉代谢过程相关的基因（GWD/ISA）主要在F1阶段表达（图3d；图S7，见在线补充材料）。在F3和F4阶段鉴定出多个高度特异性表达的基因副本，这些基因涉及单糖或多糖的代谢，包括RmAGL-1、RmAGL-2、RmAMY-1和RmBAM-5。相应地，模块25中与糖苷代谢相关的基因在F4阶段富集（r = 0.83；图3A）（图S6，见在线补充材料）。
此外，参与淀粉降解的那些基因在叶片或茎的衰老过程中也高度表达。这与多种物种中的器官老化过程一致。
与桃金娘果实发育中的着色和花青素合成相关的代谢物和基因表达模式
果实成熟过程中的着色是桃金娘的一个显著特征（图4a）。形态观察和总花青素含量测定显示，桃金娘果实中的花青素含量在F3到F4期间急剧增加，果实颜色也变为紫色（图4a和b）。相应地，模块24中的基因在与黄酮合成和代谢，尤其是花青素合成和代谢相关的途径中高度富集（图3b）。
代谢物丰度聚类和类黄酮（主要是花青素分支）化合物合成途径。A 桃金娘果实在不同发育阶段的图像。B 桃金娘果实发育期间总花青素含量的趋势。C 基于K-means的代谢物丰度聚类。D 桃金娘果实发育期间花青素含量的变化趋势。E 类黄酮合成途径、代谢物含量分布及相应的基因表达、位点和桃金娘与番石榴的同源基因共线性。
为了研究桃金娘果实着色与类黄酮积累之间的机制，进行了类黄酮的代谢组学研究。通过UP-MS检测了来自桃金娘三个器官六种样品中的189种类黄酮。这些代谢物被分为九个丰度聚类（图4c；表S8，图S8，见在线补充材料）。我们重点关注了在成熟果实中积累的化合物（cluster8），包括20种化合物，如花青素糖苷、表儿茶素和槲皮素。成熟桃金娘果实中的花青素（cluster8）主要以四种糖苷的形式存在，包括花青素-3-O-葡萄糖苷、芍药苷-3-O-葡萄糖苷和芍药苷-3-O-葡萄糖苷。天竺葵中的花青素在果实中含量太低，无法检测到。
为了研究共表达模块与代谢组学数据之间的关系，我们分别将果实中代谢物的含量和积累率与不同模块中的基因表达进行了相关性分析。模块25与cluster8中的代谢物含量之间检测到最高相关性（图S9，见在线补充材料）。功能分析表明这些基因主要富集在19个GO术语中，包括三个与液泡相关的术语（图S6b，见在线补充材料）。先前的研究表明，果实颜色变化与代谢物向液泡的运输密切相关[18]。F4中花青素的积累是否与该模块中的基因相关值得进一步探讨。更重要的是，模块24中的基因表达与cluster8中代谢物的积累率密切相关（平均R = 0.873）。这与F3相关模块（模块24）的GO富集结果相对应，表明模块24中的基因在桃金娘果实中花青素积累和着色中起重要作用。
我们进一步基于KEGG中记录的花青素生物合成途径研究了基因表达模式和代谢物波动。根据记录，共有八种类型的基因被认为编码花青素合成中的核心酶。在果实着色过程中，尽管模块25中的基因表达趋势与20种代谢物和总花青素的含量趋势差异密切相关，但在花青素合成途径中仅检测到一个基因Rm4CL-2。相比之下，在花青素合成途径中至少有六种类型的酶（4CL、CHS、CHI、F3H、ANS、DFR）在模块24中聚集，并在F3中上调。在果实发育过程中，显示在F3上调的两种酶（PAL、C4H）也检测到相应的基因副本（RmPAL-1、RmC4H-1）。在花青素合成的核心过程中（从柚皮苷到花青素），包括RmCHI-1、RmCHI-2、RmF3H、RmDFR和RmANS在内的所有已知核心酶在F3阶段高度表达。这一趋势与桃金娘果实从F3到F4（开花后75至90天）的颜色和花青素含量的急剧变化和增加相对应（图4c和d）。
桃金娘科中的果实着色和花青素积累
桃金娘在桃金娘科中独特果皮颜色背后的自然选择进化轨迹尚不清楚。番石榴（P. guajava）是另一种具有绿色果皮和果肉的桃金娘科肉质果实物种[9]，是进行比较分析的最合适物种。与桃金娘不同，转录组分析检测到参与花青素核心合成的基因在番石榴果实中没有上调，而是下调（与桃金娘中的同源基因呈负相关；表S9，见在线补充材料）。这种表达波动的差异表明花青素合成活性在番石榴果实中没有增强，这与番石榴果实的淡绿色一致（图4D和E）[9]。
此外，比较研究两种物种之间的花青素合成下游途径报告了OMT基因的拷贝数变异（CNV），这些基因需要进行花青素的糖基化。在染色体1末端定位的串联重复扩展结果中有两个基因副本（RmOMT-4/RmOMT-5）（图5A）。这些基因副本与番石榴中的同源基因（Pgu20700）具有高度的序列同源性和结构相似性。表达分析显示，RmOMT-4的表达在果实中特异性增加，而RmOMT-5在所有果实成熟阶段都有表达，但最高表达出现在根部（图4d）。同时，花青素糖苷的运输和积累极大地影响了植物的颜色。先前报道，GST超级基因家族参与了花青素的运输过程[18]。基于KEGG注释，我们在桃金娘中检测到更多的GST基因（76 > 51）相比番石榴（图5B）。基于系统发育关系的亚家族分析表明，桃金娘中GSTU亚家族基因比番石榴中更多（56 > 32；图5B）。这个亚家族被证明参与花青素的运输[19]。CNVs与花青素糖基化以及运输的潜在差异之间的关系，进一步影响了桃金娘和番石榴果实的着色，值得进一步研究。
基因拷贝数、MYB和GST基因家族系统发育、MYB差异表达基因以及正选择基因的分析。A 与番石榴中的同源基因相比，桃金娘在Chr1上多出一个OMT基因拷贝。B 基因家族分析表明，与番石榴相比，桃金娘中的GST家族基因拷贝数增加（76 > 51）。C 桃金娘模块24中MYB基因及其在番石榴中的同源基因的表达。D 高度相关于RmPAP-1和RmPAP-2的模块24中的相关网络。花青素合成途径中的基因节点用绿色表示，边缘的透明度表示两个基因之间的权重值。节点的大小代表连接基因的数量。E 正选择基因及其在桃金娘中的分布。
正选择的MYB基因对桃金娘花青素合成的潜在影响
本研究进一步分析MYB基因家族，检测到桃金娘和番石榴中的MYB基因数量相似（251/249）（图S10，见在线补充材料）。其中，4个MYB基因的表达在WGCNA中聚集于模块24。相比之下，这四个基因的同源基因在番石榴果实中未检测到类似的特定高表达时期（图5C）。系统发育分析表明，这四个基因中的两个（RmPAP1和RmPAP2）与拟南芥中的PAP1（MYB75）、PAP2（MYB90）以及AtMYB113聚集在一起（图S10，见在线补充材料），这些基因的高表达已被证明能促进花青素的合成。这两个MYB转录因子的共表达网络包含花青素合成途径中的九个基因（图5D）。
同时，在桃金娘和其他桃金娘科物种中检测到的正选择基因表明，桃金娘中有2432个基因受到正选择。这些基因涵盖了32个MYB转录因子，包括RmPAP-2和RmMYB113。物种间表达趋势的差异和花青素正调控的MYB基因在进化中的正选择倾向导致了桃金娘果实的着色（图5E；表S9，见在线补充材料）。其余的正选择MYB转录因子倾向于参与植物的多个生物过程，包括花粉发育（MYB80）、种子萌发（MYBC1）和表皮发育（ETC2）（表S10，见在线补充材料）。
桃金娘叶和茎中的类黄酮药用活性
先前的研究表明，从相应药用植物的叶或茎中提取了具有药用价值的多种类黄酮代谢物。在本研究中，作为花青素前体的二氢杨梅素在植物叶片中的含量显著升高，而其波动与下游花青素不同。二氢杨梅素还被报道为另一种药用类黄酮杨梅素的直接前体，其含量在叶片中也有所增加。相关性分析表明，二氢杨梅素形成和杨梅素转化所需的RmF3' 5' H和RmFLS的表达分别与其代谢产物的含量波动密切相关（表S11，见在线补充材料；图4D）。这表明叶片中的更多二氢杨梅素倾向于参与杨梅素的积累。进一步的研究检测到在具有药用价值的多种类黄酮代谢物中存在类似机制，如在茎中含量较高的木犀草素和在叶中含量较高的槲皮素。在具有巨大利用潜力的草本茶中也发现了这种药用代谢物的类似分布模式[20]。需要更多实验证据和研究来验证和利用桃金娘中的这种模式（图4D）。
此外，在桃金娘中检测到四种B型异构体（B1-B4）的原花青素，这是一种类似花青素的抗氧化活性物质。有趣的是，两种异构体（B2/B4）的含量仅在F4的果实中升高，而其余两种仅在叶片中高度检测到。这些异构体在不同器官中的含量差异背后的隐藏功能变异值得进一步研究（图4D）。
讨论
根据先前的统计数据，桃金娘科中有超过3000个物种主要分布在世界的热带和亚热带地区[21]。桃金娘是该科中最具代表性的物种之一，是常见于东亚和东南亚的常绿灌木[3]。具有高质量和代表性的基因组资源可以促进桃金娘和桃金娘科的分子育种和进化研究。在本研究中，通过结合PacBio HiFi和ONT数据以及Hi-C数据，发布了无缺口端粒到端粒的桃金娘基因组。完整的基因组组装无疑将成为该物种和桃金娘科进一步研究的基准。本研究发布的基因组信息、代谢物丰度和基因表达模式也将成为分子育种和改进收获后储存策略的宝贵遗传资源。
水果和蔬菜中的花青素和其他类黄酮可以预防多种疾病，尤其是心血管疾病和某些类型的癌症。同时，花青素被证明在给人类和动物服用时显著改善视力[22]。先前关于重要药用和营养特性的研究表明，类黄酮化合物，尤其是花青素，是桃金娘果实中的重要生物活性物质[3]。这与相应途径合成基因的高表达一致，也促进了桃金娘果实的着色。代谢研究揭示了这一过程归因于四种花青素的积累。由于含量极低，未检测到天竺葵中的花青素积累，这与先前的定量研究一致[2]。
花青素的高合成和积累通常是由花青素合成途径中编码核心酶（CHI、F3H、F3' H、DFR、ANS）的基因高表达引起的[23]。基因表达和代谢物波动分析表明，该途径在桃金娘的F3中被激活。先前的研究证明了这一机制在各种由MYB家族基因代表的转录因子调控的广泛存在[5]。比较桃金娘和番石榴之间的途径和表达趋势表明，这两种物种中与花青素合成相关的核心基因表达趋势存在显著差异，这与它们果实的颜色差异一致。基因家族和基因表达的进一步分析表明，这两种桃金娘科物种之间的分歧是由正选择以及包括PAP1、PAP2、MYB113和MYB114在内的多个MYB转录因子的各种表达导致的，这些基因在拟南芥中已被证明能促进花青素合成[24, 25]。类似的调控机制在胡萝卜和萝卜等多个物种中被检测到[26, 27]。需要更多的比较基因组和生理学实验来证明这是否是导致物种间花青素积累差异的普遍且重要的途径。
此外，糖基化和运输对于在植物器官中积累花青素和颜色变化是必要的。OMT和GST家族的基因是这个过程的关键因素[6, 7]。先前的研究报道，GST家族基因的丧失会严重影响花和桃皮中花青素的积累，导致白色的花和果皮[28]。在本研究中，检测到桃金娘中RmOMT基因拷贝在果实成熟过程中表达增加。与番石榴相比，桃金娘中OMT和GST家族的基因拷贝数增加。CNV影响的形态特征，如人类的疾病抗性或水稻的粒大小[29, 30]，已被频繁报道。在本研究中发现的CNV是否对桃金娘中的花青素代谢产生影响值得进一步研究。
桃金娘科肉质果实成熟伴随着淀粉降解和细胞壁裂解[9]。种子传播的优化取决于果实成熟，这由复杂的转录因子和遗传调控网络控制。这个调控过程与植物激素的作用密切相关。RNA-Seq分析表明，这一机制也有助于桃金娘的果实发育。这与桃金娘科另一种肉质果实物种番石榴的果实软化过程相对应。类似的机制在先前的水果发育研究中被检测到，如芒果（漆树科）、草莓（蔷薇科）和梨（蔷薇科）[9]。
结论
在这里，我们展示了桃金娘科家族中的第一个无缺口T2T基因组——桃金娘基因组及其详细的基因组信息。我们确定了桃金娘中的主要花青素化合物及其合成途径。基因表达模式分析和途径鉴定进一步为桃金娘科的肉质果实发育提供了新的见解。比较基因组和基因表达分析提供了花青素积累和果实着色的可能机制。我们的基因组组装为研究桃金娘科肉质果实的起源以及加速桃金娘的遗传改良奠定了基础。
材料和方法
桃金娘样品收集和测序
从中国广东省惠州市博罗区（23.238691°N，114.044166°E）收集了‘LFSTJN-1’的根、茎、叶、花和果实组织样品（图1A）。使用QIAGEN Genomic Kit提取高质量基因组DNA，通过0.75%琼脂糖凝胶电泳检查DNA的完整性。分别使用Nanodrop（OD260/280 = 1.8–2.0，OD260/230 = 2.0–2.2）和Qubit检测DNA的纯度和数量。
转录组和代谢组样品涵盖了五种组织的十个和六个时期（表S1，见在线补充材料）。
Hi-C实验
使用NEBNext Ultra II DNA library preparation kit和DpnII酶（Ipswich, MA, USA）构建新鲜嫩叶的Hi-C文库。使用Illumina NovaSeq 6000平台对Hi-C文库进行测序。此过程中生成了59.45 GB的清洁数据（表S1，见在线补充材料）。
基因组调查、染色体组装构建和缺口填补
通过软件jellyfish和Genomescope评估桃金娘的基因组大小和杂合性[31]。对于PacBio组装，使用CCS模式（HiFi reads）生成的PacBio reads通过软件Hifiasm进行组装，使用默认参数[32]。NextDenovo（https://github.com/Nextomics/NextDenovo）用于组装ONT reads，参数为：genome_size = 460 M，read_cutoff = 50 000，seed_cutoff = 55 959，seed_depth = 45。使用Hi-C数据分离嵌合体，通过3D-DNA处理并通过juicerbox可视化[33, 34]。使用HiFi-based组装填补ONT组装参考的缺口。
重复元素识别、基因预测和功能注释
通过Pipeline EDTA进行重复元素识别和屏蔽。屏蔽的基因组序列用于进一步的基因预测[35]。
通过整合de novo预测、同源蛋白证据和EST证据，使用pipeline Maker进行基因预测[36]。de novo基因预测使用软件Augustus，而hisat2、stringties和PASA用于生成186.71Gb RNA-seq数据的EST证据[37]。从桉树和番石榴中导入同源蛋白序列到exonerate模块以生成同源证据。PASA还用于拼接位点调整和UTR预测。
功能注释通过使用NCBI BLAST与公共数据库比较预测蛋白质并使用1e-5的截止e值实现。最佳匹配的BLAST结果被认为是基因功能。GO和KEGG的注释通过使用EggNOG和KAAS获得[38]。端粒检测和着丝粒定位通过结合TRF和CD-HIT基于先前研究进行处理[16, 39, 40]。
基因组组装质量评估
采用多种方法评估组装基因组的准确性和完整性。首先，我们计算了以N50和LAI评分表示的染色体长度和基因组组装完整性的基本指标[41]。此外，将配对末端短读数和HIFI读数映射到基因组，使用BWA-MEM，同时使用hisat2对RNA-seq数据进行比对以评估其完整性。从遗传角度，使用数据库eudicots_odb10和embryophyta_odb10的基准单拷贝直系同源基因（BUSCO）评估基因组的完整性[42]。
系统发育分析、分化时间估计和基因家族分析
通过OrthoFinder2确定10种物种的基因直系同源物，包括单拷贝直系同源基因[43]。使用MAFFT完成单拷贝直系同源物的编码DNA序列（CDS）的比对，选项为-maxiterate 1000[44]。通过软件RaxML在GTRGAMMA模型和100次bootstrap下用最大似然（ML）方法推导物种间的系统发育树[57]。使用PAML软件的MCMCtree模块确定每个树节点的分化时间[45]。系统发育过程的校准使用了各种化石记录。通过在分裂节点上使用来自timeree数据库的记录（TimeTree :: The Timescale of Life）放置软边界来估计分子分化时间。使用默认参数的程序CAFE识别每个分化节点和物种的基因家族扩展和收缩[46]。
共线性、正选择基因和WGD检测
使用python版本的MCScanX（JCVI）在桃金娘科物种和石榴之间识别共线性块[47] [9, 48, 49]。通过使用Orthofinder2（e值≤1 × 10^–5）执行多序列比对以提取此数据的保守同源物蛋白序列[43]。使用这一过程生成的3372个单拷贝直系同源物识别物种之间的系统发育关系以检测桃金娘中的正选择基因。使用PAML的codeml模块计算桃金娘分支的选择压力，采用分支位点模型[45]。随后，通过卡方检验（P值<1e-5）验证正选择基因的显著性。
对于WGD的检测，使用WGDI识别和提取共线性间隔和相应的基因对。通过NG86算法的方法计算同义替换每同义位点（Ks）以预测WGD事件[50]。基于公式：分化时间 = Ks/（2 × r）估计WGD的年龄[58, 59]，通过直系同源基因的Ks分布推断桃金娘目的r（植物平均Ks/年速率）。
转录组数据处理
转录组数据的比对通过pipeline Hisat2处理，数据排序和文件格式转换使用Samtools。使用软件FeatureCounts生成基因表达统计数据[51]。R模块DESeq2用于数据归一化和鉴定差异表达基因（DEGs）[52]。番石榴器官和发育阶段的转录组数据从NCBI获取，Bioproject PRJNA631442 [9]。
WGCNA和相关性分析
WGCNA用于识别与桃金娘器官相关的基因模块[53]。首先，我们基于基因表达水平对样本进行聚类。然后，使用软件包的TOMsimilarity模块和PickSoftThreshold函数生成基因之间的共表达相似系数和加权过程。这种相关性进一步转化为基因之间的关联。我们使用绝对值>0.8和P值<0.001来识别与器官显著相关的基因模块。其余细节基于先前研究中的协议[54]。
R包Hmisc（CRAN - Package Hmisc）用于评估基因、基因模块和代谢物聚类之间的相关性，使用皮尔逊相关系数。特别是，果实中代谢物积累率的计算如下：Vt, n = (Mt + 1, n - Mt, n) / (Tt + 1 - Tt). 上述公式中的符号：n：代谢物类型；t：采样时间点（本研究中的值范围为1 ~ 3）；M：化合物含量；T：每个采样点对应的开花后时间。
MYB和GST基因家族分析
使用基于PF00249的hmmsearch从Pfam（InterPro）中识别桃金娘和番石榴中的MYB基因和GST基因，并使用eggNOG进行功能注释。为了进一步分析这两个家族基因的系统发育关系和识别亚家族，通过iQtree构建系统发育树[60]。基于功能注释信息，将拟南芥中的两个家族基因提供用于序列比对和进一步的系统发育树构建作为亚家族分类的基础。
类黄酮相关化合物的提取
生物样品在真空冷冻干燥机（Scientz-100F）中冻干。使用混合研磨机（MM 400，Retsch）和锆珠以30 Hz研磨1.5分钟将冻干样品粉碎。50毫克的冷冻干燥粉末溶解于1.2毫升70％甲醇溶液中，每30分钟涡旋30秒，共涡旋六次。以12000 rpm离心3分钟后，过滤提取液（SCAA-104，孔径0.22μm，ANPEL，上海，中国，上海安谱实验科技股份有限公司官网），然后进行UPLC-MS分析。
代谢物的定量、表征和变异分析
基于先前研究中的协议进行桃金娘中类黄酮相关代谢物的定量和表征过程[55]。使用k-means算法进行代谢物含量的聚类。R包vegan中的cascadeKM函数用于识别最合适的聚类数[61]。
淀粉和总花青素的测定
收集不同发育阶段的果实，称重，并在4°C下过夜用1毫升100%冷却的甲醇进行恒定振荡提取。以1300 tpm离心10分钟后，吸取200μL上清液，并在A650和A666波长下记录数据。此外，吸取200μL上清液用于花青素测定，加入1%体积的浓盐酸，充分混合。分别在A530和A657波长下记录数据，并使用公式((A530–0.25*A 657)/鲜重)计算花青素值[56]。淀粉提取和定量检测由Solarbio Life science BC0700试剂盒（北京，中国）执行。
