摘要
主要结论：
一种基于单倍型解析的、染色体级别的AAC Mountainview紫云英基因组序列，结合了细胞器组装、甲基化组以及组织/胁迫特异性的转录组数据，为后续的紫云英育种工作提供了坚实的基础。
紫云英（Onobrychis viciifolia Scop.）是一种异交繁殖的多年生牧草作物，因其豆科特性、营养价值、适口性以及能够产生缩合单宁而受到重视。这些缩合单宁可以降低反刍动物发生牧场胀气的概率，并提高饲料效率。然而，尽管具有这些优势，针对该物种的研究和育种工作仍然相对较少。在本研究中，我们报告了利用PacBio HiFi、Oxford Nanopore、Illumina短读长序列和Hi-C数据生成了一种四倍体紫云英品种AAC Mountainview的单倍型解析染色体级参考基因组序列。该基因组序列包含了所有28条假染色体（对应于7条基染色体上的4种单倍型），具有较高的连续性和基因完整性，并且所有单倍型的长末端重复序列（LTR）组装指数（LAI）得分均达到参考标准。基于Nanopore技术的甲基化分析显示了典型的基因体CG甲基化现象以及高水平的转座元件甲基化。我们注释了117,890个高置信度的蛋白质编码基因，并识别出位于非锚定但与染色体相关的rDNA串联阵列。此外，还成功组装了该品种的线粒体和叶绿体基因组。组织特异性和胁迫特异性的转录组分析揭示了不同组织和胁迫类型下的共有和独特的基因表达模式。等位基因特异性表达分析表明，在胁迫条件下单倍型的活性基本保持平衡，但存在细微但一致的表达优势变化。本研究提供的数据为这种有前景的牧草作物的后续育种工作提供了宝贵的资源。
引言
紫云英（Onobrychis viciifolia Scop.; 豆科）是一种多年生豆科牧草作物，以其营养价值、适口性和能够产生缩合单宁而闻名（Sheppard等人，2019年）。它分为两种类型：一种是“普通型”（O. sativa var. communis Ahlefed；原产于中欧），这种类型在首次收割后不再开花；另一种是“巨型”（O. sativa var. bifera Hort.；原产于中东），这种类型在再生过程中会再次开花（Poudel等人，2023年）。虽然紫云英在欧洲、亚洲和北美已有数百年的栽培历史（例如Mora-Ortiz和Smith，2018年），但由于集约化农业方法的采用以及高产牧草品种（如紫花苜蓿Medicago sativa L.）的广泛种植，其产量在过去一个世纪里大幅下降（Hayot Carbonero等人，2011年）。然而，近年来，作为一种可持续的半干旱环境牧草作物，紫云英再次受到关注，特别是在加拿大草原地区（例如Sheppard等人，2019年；Poudel等人，2023年）。
这种关注部分源于紫云英在多种组织中合成缩合单宁（CTs；也称为原花青素），这与仅在其种皮中产生大量CTs的紫花苜蓿不同（Goplen等人，1980年）。在反刍动物体内，这些CTs会在酸性的瘤胃中与蛋白质结合，随后在pH值较高的皱胃和小肠中释放出来，从而促进消化和吸收（Mueller-Harvey等人，2019年）。这有助于减少牧场胀气的发生（Min等人，2006年），提高饲料中的氮利用率，并降低牲畜的氨排放量（Mueller-Harvey，2006年）。此外，研究表明，牧草中的CTs还具有驱虫效果（Desrues等人，2016年），并可能减少食用它们的反刍动物的甲烷排放（Hatew等人，2016年）。尽管这些特性非常有益，但紫云英的产量通常低于紫花苜蓿，且在收割或放牧后的再生能力较差（例如Lauriault和Marsalis，2024年）。虽然紫云英由于其深根系能够在一定程度上耐受水分不足（Poudel等人，2023年），但其抗旱性似乎不如紫花苜蓿（例如Biligetu等人，2014年、2021年），在这方面仍有很大的改进空间。此外，紫云英对水涝/洪水特别敏感（Heinrichs，1970年），这大大限制了其种植面积。尽管目前的育种工作已经在某些方面有所改进（Acharya，2015年；Monirifar，2022年；Poudel和Acharya，2022年），但与其他牧草豆科植物（如三叶草属Trifolium spp.）相比，针对该物种的研究仍然非常有限（Poudel等人，2023年），用于指导育种工作的分子信息也相对匮乏。这种基因组资源的缺乏限制了对CT生物合成机制等性状的深入研究，以及提高其耐放牧性、再生能力和非生物胁迫耐受性的努力。
尽管最近发布了一种未具体说明的“普通型”紫云英品种的基因组序列（He等人，2024年），但大多数现代品种都是从“普通型”和“巨型”类型衍生而来的（Mora-Ortiz和Smith，2018年），因此该品种的遗传背景可能与目前广泛种植的品种有很大差异。为了提供一个更具代表性的参考序列，我们致力于为四倍体（2n = 4x = 28）紫云英品种AAC Mountainview生成基因组序列。该品种在加拿大通过多种品种/种质的杂交培育而成，具有更强的耐刈割性、更好的再生能力和更高的牧草产量（Acharya，2015年；Poudel等人，2023年）。由于紫云英是一种遗传复杂、必须异交繁殖且高度杂合的物种（Kempf等人，2015年），这一过程较为复杂。通过结合使用PacBio HiFi和Oxford Nanopore长读长序列、Illumina短读长序列以及来自不同组织的短读长和长读长转录组数据，我们成功生成了一个高度连续、考虑单倍型信息的参考基因组。同时还组装了叶绿体和线粒体基因组，为未来的系统发育和功能研究提供了完整的细胞器参考。我们还分析了全基因组DNA甲基化情况，并研究了组织特异性的基因表达、等位基因的优先表达以及非生物胁迫下的转录水平响应。这种基于单倍型解析的参考基因组，结合广泛的基因注释和表达数据，为紫云英的分子育种、基因组编辑、泛基因组评估和功能遗传学研究提供了新的机会，从而促进其农艺改良。
材料与方法
植物基因型和生长条件
所有基因组和RNA测序均使用来自单一基因型（通过无性繁殖获得）的O. viciifolia品种AAC Mountainview的植物进行（Acharya，2015年）。为了在整个实验过程中保持和扩增该基因型，采用了根插条的无性繁殖方法。植物在控制温室条件下生长，光照周期为16小时/8小时，昼夜温度约为22°C/18°C。在干旱处理中，土壤含水量先被调整至约50%，之后停止供水。土壤含水量使用ML3 ThetaKit土壤湿度计（Hoskin Scientific Ltd.，加拿大不列颠哥伦比亚省本纳比）进行测量。在水涝处理中，土壤含水量同样被调整至约50%，之后将花盆浸入水中至花盆高度的3/4。正常浇水的植物作为对照组。
Illumina、Oxford Nanopore和PacBio HiFi全基因组测序
对于Oxford Nanopore和Illumina测序，使用基于盐析的方法从选定的AAC Mountainview基因型的幼叶组织中提取高分子量基因组DNA（Doyle和Doyle，1987年）。根据制造商的建议（Oxford Nanopore Technologies，英国牛津），使用Ligation Sequencing Kit V14（SQK-LSK114）制备Nanopore文库，并使用R10.4.1流式细胞器（Oxford Nanopore Technologies）进行测序。Nanopore碱基调用使用Dorado版本0.9.6（https://github.com/nanoporetech/dorado）以双链模式完成，包括5-甲基胞嘧啶（5mC）、5-羟甲基胞嘧啶（5hmC）和N6-甲基腺嘌呤（6mA）的修饰。短读长Illumina数据由服务提供商在Illumina NovaSeq 6000平台上生成（150 bp配对末端读长；Illumina Inc.，美国加利福尼亚州圣地亚哥），使用Novogene的标准DNA文库构建流程（Novogene Corporation Inc.，美国加利福尼亚州萨克拉门托）。读长经过fastp v0.23（Chen等人，2018年）进行适配器修剪和质量过滤。
对于PacBio HiFi测序，从经过黑暗处理（53小时）的幼叶组织中提取样本，并迅速冷冻在液氮中。使用PacBio提供的LN2协议（https://www.pacb.com/wp-content/uploads/Procedure-checklist-Isolating-nuclei-from-plant-tissue-using-LN2-disruption.pdf）分离植物核，然后使用Nanobind Plant Nuclei试剂盒（Pacific Biosciences，美国加利福尼亚州门洛帕克）提取高分子量DNA。DNA的长度、质量和浓度通过脉冲场凝胶电泳和Qubit Flex荧光仪（Thermo Fisher Scientific Inc.，美国马萨诸塞州沃尔瑟姆）进行测定。文库使用SMRTbell Express Template Prep Kit 2.0（Pacific Biosciences）制备，并由服务提供商（Novogene Corporation）在PacBio Revio平台上进行测序，以生成HiFi读长。
Hi-C测序
从选定的AAC Mountainview基因型中采集幼叶组织，迅速冷冻在液氮中，并置于-80°C。分别使用Arima-HiC试剂盒（Arima Genomics，美国加利福尼亚州卡尔斯巴德）和Phase Genomics试剂盒（Phase Genomics，美国华盛顿州西雅图）与DpnII和HinfI以及DdeI、MseI和HinfI一起制备了两个Hi-C染色质相互作用文库。测序分别在Illumina NovaSeq 6000平台和Phase Genomics公司内部完成。
基因组大小和倍性估计
使用KMC v3（Kokot等人，2017年）从PacBio HiFi读长中计算21-mer频率，移除零计数条目，并使用GenomeScope v2.0（Ranallo-Benavidez等人，2020年）确定基因组大小、重复序列含量、k-mer和测序错误率，倍性设置为4。使用Smudgeplot v0.2.5（Ranallo-Benavidez等人，2020年）分析k-mer计数和等位基因多态性模式。覆盖度截止值通过cutoff函数自动确定（L = 30，U = 4700）。
基因组组装和支架构建
使用PacBio HiFi、Oxford Nanopore双链和Nanopore单链读长，分别设置最小长度为1,000 bp、1,000 bp和20,000 bp，以及最小质量分数为20、20和10。使用Verkko（Rautiainen等人，2023年）结合PacBio HiFi、Oxford Nanopore和Hi-C数据进行组装。双链Nanopore读长被指定为高精度（-hifi），同时使用PacBio HiFi读长。单链读长仅用于图形遍历（-nano）。组装图使用Bandage（Wick等人，2015年）进行检查。Hi-C读长使用Cutadapt（Martin，2011年）进行修剪，并使用Arima Genomics映射管线（https://github.com/ArimaGenomics/mapping_pipeline）进行对齐、过滤和合并。支架构建使用HapHiC和标志–RE GATC、GANTC（Zeng等人，2024年）完成。支架使用Juicebox（Robinson等人，2018年）进行手动整理。单倍体和小的未支架化contigs在Hi-C接触图和组装图连续性的支持下进行重新分配。剩余的重复序列（约200 Mb）根据组装图分配到染色体组（chr2un至chr7un）。使用Pilon（Walker等人，2014年）对组装结果进行Illumina polishing处理。完整性和连续性评估使用了BUSCO v5.3.2（embryophyta_odb10, fabales_odb10；Simao等人，2015年）。RNA-Seq和Iso-Seq分析分别针对叶片、茎、根和花组织进行，这些组织是在切割后大约3个月时采集的，此时在我们的温室条件下植物通常已经开始开花。组织样本来自4株具有相同AAC Mountainview基因型的生物重复植株，这些植株都是通过茎插条繁殖得到的独立克隆体。共收集了9种类型的组织，并立即用液氮快速冷冻保存，包括刚长出的叶子（仍然紧密折叠）、正在发育的叶子（已经展开但较小）、第一片完全展开的叶子（从茎尖开始数的第三片）、第三片完全展开的叶子、茎（从茎尖开始数的第三个节段）、花蕾（未开放）、开放的花朵、靠近根冠的根以及靠近根尖的根。同样，我们也从分别经历了7天干旱和积水处理的植株中采集了第三片完全展开的叶子，在这个时间点，植株开始显示出压力迹象（干旱时叶子变干，积水时出现黄化）。对于干旱处理，组织是在土壤湿度在2.6%到4.9%之间时采集的（而正常灌溉的植株土壤湿度约为50%）。压力和非压力组织的采集是在同一时间进行的。所有样本的总RNA都是使用Sigma Spectrum Total RNA试剂盒（协议A）提取的，并进行了柱上DNase处理（Sigma-Aldrich公司，美国密苏里州圣路易斯）。RNA随后根据制造商的说明使用Norgen RNA Clean-up and Concentration试剂盒（Norgen Biotek公司，加拿大安大略省索罗尔德）进行了纯化。RNA的质量、纯度和数量通过琼脂糖凝胶电泳以及NanoDrop 8000分光光度计（Thermo Fisher Scientific公司）、Qubit Flex荧光计和Agilent 2100 Fragment Analyzer（Agilent Technologies公司，美国加利福尼亚州圣克拉拉）进行了确认。
对于RNA-Seq和Iso-Seq分析，文库制备和测序工作由一家服务提供商（Novogene Corporation）完成。在RNA-Seq的情况下，使用Novogene的RNA-Seq文库构建流程制备了双链文库，并在Illumina NovaSeq 6000平台上进行了测序（150 bp配对末端读长）。读长经过fastp v0.23（Chen等人，2018年）进行质量过滤和修剪，然后使用HISAT2 v2.2.1（Kim等人，2019年）与我们的AAC Mountainview基因组序列比对。基因水平的计数使用StringTie v2.2.1（Pertea等人，2015年）进行量化。对于Iso-Seq分析，将每种组织类型的总RNA等摩尔比混合后进行测序。文库使用Kinnex Full-Length RNA试剂盒（Pacific Biosciences公司）制备，并在PacBio Revio系统上进行测序以生成HiFi读长。原始HiFi BAM文件使用PacBio Iso-Seq v3工具链（https://github.com/PacificBiosciences/IsoSeq）进行处理。
转座元件注释
转座元件（TEs）的注释使用了RepeatMasker v4.1.2和EDTA v2.2.2流程（Ou等人，2019年），该流程结合了从基因组组装中生成的定制重复文库，通过结构发现长末端重复（LTR）逆转录转座子、末端反向重复/微型反向重复（TIR/MITE）转座元件和Helitrons来进行新注释，随后进行分类和全基因组掩蔽。转座元件的分布和分类使用LTR组装指数（LAI）和LTR_retriever（Ou和Jiang，2018年）进行了总结，并评估了重复区域的基因组连续性。
基因和非编码RNA预测及功能注释
基因模型使用BRAKER3（Gabriel等人，2024年）在蛋白质-RNA模式下进行预测，结合了OrthoDB中的Fabales蛋白质数据（Zdobnov等人，2021年），一次运行中整合了9种组织的PacBio Iso-Seq数据，另一次运行中整合了9种组织和2种压力条件下的Illumina RNA-Seq数据。结果使用TSEBRA（Gabriel等人，2021年）进行了合并。基因模型质量使用BUSCO（Benchmarking Universal Single-Copy Orthologs）v5.4.3（fabales_odb10）进行了评估。功能注释使用了DIAMOND BLASTP（Buchfink等人，2021年）在NCBI nr（https://www.ncbi.nlm.nih.gov/）、Swiss-Prot（Boutet等人，2007年）和拟南芥TAIR（Rhee等人，2003年）数据库中进行搜索。结构域预测和同源性分配使用InterProScan v5.52（Jones等人，2014年）和eggNOG-mapper v2.1.9（Cantalapiedra等人，2021年）进行，分类范围限制在植物上。eggNOG-mapper中的基因本体（GO）术语和Pfam结构域被整合到综合注释集中。基因结构统计信息，包括每个转录本的外显子数量、每个基因的异构体数量以及平均外显子和内含子长度，直接从最终的GFF3注释中计算得出，通过总结每个染色体上的外显子、内含子和转录本特征得出。编码核糖体RNA（rRNAs）、转运RNA（tRNAs）和微小RNA（miRNAs）的序列使用barrnap（https://github.com/tseemann/barrnap）、tRNAscan-SE v2.0（Chan等人，2021年）和Infernal（使用Rfam协方差模型，Nawrocki和Eddy，2013年）进行识别。
DNA甲基化检测与分析
全基因组DNA甲基化谱是通过Oxford Nanopore读长生成的。读长使用Dorado aligner（Oxford Nanopore Technologies）与参考基因组比对，该工具保留了修饰标签。修饰频率数据使用modkit pileup工具汇总，得到BEDMethyl格式的每个位点的修饰和未修饰碱基计数。这些文件按基序分割并转换为bigWig格式，以便后续分析和可视化。使用deepTools v3.5.5（Ramírez等人，2014年）生成了基因和TEs的平均甲基化模式meta图。
同源性、系统基因组学和共线性分析
同源性分析使用OrthoFinder v2.5.4（Emms和Kelly，2019年）进行，数据来源于sainfoin和其他18种参考被子植物（其中14种是豆科植物，4种是外群）的蛋白质、编码序列和GFF3文件。233个高占用度同源群的物种树使用MAFFT v7.490（Katoh和Standley，2013年）进行比对，使用FastTree2（Price等人，2010年）进行系统发育重建，并使用STRIDE（Emms和Kelly，2017年）以拟南芥为根。共线性块使用DIAMOND BLASTP（Buchfink等人，2021年）和MCScanX（Wang等人，2012年）识别，并使用TBtools v1.123（Chen等人，2023年）分析同义替换率（Ks, Ka, Ka/Ks）。
细胞器基因组组装
叶绿体和线粒体基因组使用GetOrganelle（Jin等人，2020年）、Tiara（Karlicki等人，2022年）和TIPPo（Xian等人，2025年）以及Illumina和PacBio HiFi读长进行组装和环化。尽管我们的PacBio HiFi文库是从纯化的细胞核中制备的，但在植物细胞核制备过程中通常会存在低水平的质体/线粒体残留物，这些残留物常被用来从全基因组长读长数据中组装细胞器基因组（Uliano-Silva等人，2023年）。组装结果使用GeSeq（Tillich等人，2017年）和Viridiplantae细胞器参考序列进行注释，随后手动检查基因边界。编码tRNAs的基因使用tRNAscan-SE v2.0（Chan等人，2021年）进行预测。最终的特征表以GenBank格式导出，环状图使用OGDRAW（Greiner等人，2019年）生成。
差异表达和组织富集基因分析
原始基因计数矩阵使用StringTie（Pertea等人，2015年）和prepDE.py脚本根据用于RNA-Seq的9种组织类型和2种压力条件的转录本组装生成。这些矩阵随后被用作差异表达和组织富集分析的输入。差异表达分析在DESeq2 v1.34.0（Love等人，2014年）中进行。对于干旱条件下与对照组以及积水条件下与对照组相比的叶片组织的比较RNA-Seq分析，使用调整后的P值≤0.05和|log2倍变化|≥1作为差异表达基因（DEGs）的判断标准。对于组织富集基因表达分析，进行了包含组织类型作为因素的完整模型和排除组织效应的简化模型的似然比测试（LRT）。假发现率（FDR）<0.01的DEGs被保留。使用方差稳定转换获得标准化计数，并计算每个组织的平均表达量。为了识别组织富集基因，我们将每个基因分配到其表达水平最高的组织，并进一步要求最高表达组织与次高表达组织之间的log2倍变化差异≥1。将密切相关的组织合并为更广泛的类别（正在发育的叶子+刚长出的叶子=年轻叶子；第一片完全展开的叶子+第三片完全展开的叶子=成熟叶子；根冠+根尖=根），以提高统计功效和可解释性。如果至少在2个文库中每个基因的计数≥1，则认为该基因表达，如果所有样本中的变异系数（CV）≤0.20，则认为基因表达稳定。对于所有RNA-Seq分析，除非另有说明，否则基因计数指的是每个基因的各个等位基因。
差异表达和组织富集基因都进行了功能富集分析。GO术语富集使用GOATOOLS（Klopfenstein等人，2018年）进行，显著性在FDR调整后的p<0.05水平确定。使用REVIGO（Supek等人，2011年）进行富集GO术语的总结和冗余减少，推荐的不可替代性截止值为0.4。过表达通过Fisher的精确检验和Benjamini–Hochberg校正进行评估，可视化结果以R中的点图形式呈现，其中点大小反映了每个类别的DEGs数量，颜色表示调整后的显著性。对于应激响应实验，使用MapMan v3.6（https://mapman.gabipd.org/）进行通路分析，以AAC Mountainview基因组为参考，并使用Ora-Fisher测试进行PageMan分析以确定显著性，截止值为1.0。为了避免MapMan分析中的计数膨胀，我们使用OrthoFinder（Emms和Kelly，2019年）将基因模型按4个单倍型分组，然后将4个等位基因合并为一个代表。然后为每个同源群选择一个代表基因ID作为下游可视化的标识符。
黄酮类化合物通路分析
为了分析黄酮类化合物通路（包括CT通路基因），从拟南芥和V. vinifera中编译了一个包含20个基因的参考面板（表S1）。使用DIAMOND BLASTP（Buchfink等人，2021年）在sainfoin中搜索预测的完整蛋白质集合，并保留了最高得分的候选基因进行映射。然后从全局分析中提取这些候选基因的组织富集状态，并通过保留每个参考基因的前4个最高得分的同源匹配来进一步总结通路部署。为了比较AAC Mountainview基因组与先前发表的sainfoin基因组（He等人，2024年）中这些黄酮类生物合成基因的同源基因数量，使用DIAMOND BLASTP（Buchfink等人，2021年）对两个组装的预测蛋白质集合进行搜索，E值截止值为1e−10，并根据序列同一性（≥50%）和对齐覆盖率（≥70%）进行筛选。对于每个基因家族，合并了冗余的转录本异构体，仅保留每个位点的最高得分匹配。
等位基因特异性表达分析
通过识别高置信度的等位基因组来研究干旱和积水压力下的等位基因特异性表达模式。预测的蛋白质序列根据染色体标识符（例如chr1a-chr7a）分配到4个单倍型组（A、B、C、D），并使用OrthoFinder v2.5.4（Emms和Kelly，2019年）作为伪物种进行分析。选择每个单倍型中恰好包含一个基因的同源群（1:1:1:1组）进行进一步分析。过滤同源群以确保有足够的读长支持（每个单倍型总读长≥50，每种条件≥20个读长，每个基因≥10个读长，并且在每种条件下至少有2个可检测表达的重复样本）。对于每个等位基因组，等位基因特异性表达表示为相对于该组总表达的读长比例（即4个单倍型比例之和为1）。条件内的等位基因偏倚通过比较最高表达和最低表达单倍型之间的差异来衡量，接近0的值表示表达平衡，较大的值表示显性（轻度显性的阈值为0.25，强显性的阈值为0.5）。使用Python（SciPy，statsmodels）中的chi-square检验测试条件间（对照组与干旱或积水）的单倍型比例差异，P值使用Benjamini–Hochberg方法进行调整（FDR<0.05）。在MapMan分析中，使用OrthoFinder（Emms和Kelly，2019年）将基因模型按4个单倍型分组，然后将4个等位基因合并为一个代表。然后为每个同源群选择一个代表基因ID作为下游可视化的标识符。
对于黄酮类化合物通路分析（包括CT通路基因），从拟南芥和V. vinifera中编译了一个包含20个基因的参考面板（表S1）。使用DIAMOND BLASTP（Buchfink等人，2021年）在sainfoin中搜索预测的完整蛋白质集合，并保留了最高得分的候选基因进行映射。然后从全局分析中提取这些候选基因的组织富集状态，并通过保留每个参考基因的前4个最高得分的同源匹配来进一步总结通路部署。为了比较AAC Mountainview基因组与先前发表的sainfoin基因组（He等人，2024年）中这些黄酮类生物合成基因的同源基因数量，使用DIAMOND BLASTP（Buchfink等人，2021年）对两个组装的预测蛋白质集合进行搜索，E值截止值为1e−10，并根据序列同一性（≥50%）和对齐覆盖率（≥70%）进行筛选。对于每个基因家族，合并了冗余的转录本异构体，仅保留每个位点的最高得分匹配。选定基因型的测序使用了Illumina（全细胞）的短读长测序技术、Oxford Nanopore（全细胞）测序技术以及PacBio HiFi（核基因组）测序技术。Illumina测序产生了超过7亿对读长（总计14亿个读长，210Gb），相当于单倍体基因组的365倍覆盖度。在Oxford Nanopore测序中，我们获得了245.6Gb的长读长数据，平均读长在7-17kb之间，N50值为13-36kb，相当于426倍的覆盖度。此外，我们还获得了123.4Gb的PacBio HiFi长读长数据，平均读长为17.3kb，N50值为17.2kb，相当于214倍的覆盖度（见表S2，图S1）。最后，通过Illumina 150 bp配对末端读长技术共生成了4.395亿个Hi-C读长对（132.1Gb）（见表S2）。使用PacBio HiFi读长进行的K-mer分析表明，估计的单倍体基因组大小为575.8Mb，相应的四倍体基因组大小约为2.3Gb（见图1a）。杂合度估计为4.7-6.3%，主要由AAAB（2.9-3.0%）和AABB（1.5-2.0%）K-mer类别组成。4个不同的K-mer覆盖峰的存在以及Smudgeplot分析支持AAC Mountainview sainfoin具有自体四倍体特性（见图1b）。
图1：该图像的替代文本可能是通过AI生成的。
**AAC Mountainview四倍体sainfoin基因组的K-mer分析和基因组组装概览：**
a. GenomeScope分析得到的K-mer谱。
b. 使用21-mers进行的Smudge plot K-mer对频率分析。
c. 染色体级基因组组装的Hi-C接触矩阵。该矩阵显示了28个组装伪染色体之间的染色质相互作用频率，对角线上有明显的方形块，表示高质量的染色体级支架结构。对角线上的每个方块代表一个染色体，染色体间的相互作用有限。
d. 最终sainfoin伪染色体组装的Circos图。28个伪染色体以圆形布局显示，按颜色分为7组染色体（chr1a–chr7d）。显示了6个同心轨道：I. 染色体示意图，长度以Mb为单位，按染色体组进行颜色编码；II. Gypsy（红色）和Copia（绿色）转座元的密度，以堆叠直方图形式显示；III. 每1Mb窗口内的CG甲基化水平（% mC）；IV. 来自注释的基因密度；V. 作为线图的GC含量分布；VI. 来自全基因组比对的染色体间同线性链接，突出显示了染色体组之间的大规模结构一致性。为了清晰起见，仅显示了代表性A单倍型的同线性链接。
核基因组组装整合了PacBio HiFi、Nanopore和Hi-C数据，随后使用Hi-C数据进行支架构建和精细化处理，并使用短读长Illumina数据进行抛光。这一流程产生了2.15Gb的染色体级组装，涵盖28个伪染色体（占总组装量的91.1%），每个伪染色体代表7条同源染色体中的一个拷贝，其支架N50超过75Mb（见图1c，d；表1；表S3）。最终组装的GC含量为36.2%，各个单倍体组装的大小非常相似，范围从531.3Mb（单倍体C）到547.7Mb（单倍体D），平均约为541Mb，表明单倍体之间的覆盖度和连续性一致。检查Verkko组装图显示了7个不同的缠结，对应于7个不同的染色体组，并允许在支架构建后将未定位的contigs分配到染色体区间内。Hi-C接触图确认了支架结构的准确性，显示出强烈的染色体内信号和最小的染色体间相互作用（见图1c）。BUSCO分析评估了跨物种保守的单拷贝直系同源基因的存在，使用embryophyta_odb10和fabales_odb10数据集，分别显示出99.9%和97.7%的完整性，以及97.2%的重复基因（表S4）。各个单倍体的BUSCO得分相当，平均完整性约为95-97%，证实了单倍体之间的一致组装质量（表S4）。
**表1：AAC Mountainview sainfoin基因组组装指标摘要**
为了评估结构变异，我们对每个染色体的4个单倍体拷贝（例如chr1a-chr1d）之间以及我们的AAC Mountainview组装与He等人（2024年）发表的栽培品种之间进行了全染色体比对。AAC Mountainview中的单倍体间比较显示了广泛的同线性，大约85%的对齐序列是同源的，支持了这一四倍体基因组的成功分相（见图S2）。在品种间的比较中，每个AAC Mountainview单倍体与其在He等人（2024年）研究中的匹配单倍体进行了比对（例如chr1a-chr1a），同线性仍然很高，每个染色体为86.9%，全基因组为82.4%（按长度加权计算）。此外，品种间的结构变异总数仅略高于单倍体内的变异（平均每个染色体分别为4.22×10^4和4.07×10^4个结构变异）。最显著的重排包括倒位、易位和小的插入或删除，主要发生在chr5和chr6上，chr1和chr2上也观察到了额外的差异（见图S3）。此外，对于黄酮类生物合成基因家族（包括直接参与CT生物合成的基因；表S1）的同源基因数量比较，AAC Mountainview与之前发表的未知栽培品种/系谱的sainfoin基因组在大部分情况下显示出相似性（表S5）。一个显著的例外是CHALCONE ISOMERASE（CHI）家族，在AAC Mountainview中识别的同源基因数量明显多于之前报道的基因组（20个对比12个；表S5）。
**转座元件景观和LTR组装质量：**
我们共注释了2,559,894个重复元件，占组装的2.15Gb染色体基因组的64%以上（表2）。LTR逆转录转座子是最丰富的TE类别，占总掩蔽内容的38.6%，其中Gypsy占11.6%，Copia占5.4%，未分类的LTR占21.7%。TIR DNA转座子占另外12.0%的掩蔽序列，其中Mutator占4.6%，CACTA占3.1%，hAT占2.5%。Helitrons和长散布核元件（LINEs）分别占8.5%和4.5%，而短散布核元件（SINEs）和非LTR逆转录元件（例如pararetrovirus）的代表性较低。
这些完整的LTR结构使我们能够通过LAI评估基因组的连续性和组装质量。全基因组的LAI值为14.7（表S6），符合参考级组装标准（Ou等人，2018年）。各个染色体的LAI值在10.2-24.3之间，所有4个单倍体（A-D）的LAI值均高于14.7，超过了参考质量的阈值10（表S6；图S4）。这些结果共同提供了TE结构的详细视图，并证明了重复富集区域的组装连续性。
**AAC Mountainview sainfoin基因组包含大量的rDNA重复序列：**
由于rRNA编码序列的高度同质性和重复性，它们在基因组组装中仍然是一个重大挑战。我们在AAC Mountainview sainfoin基因组中共识别出74,779个rDNA特征。其中20.8%（15,524个）位于28个染色体支架中（chr1a–chr7d），但大多数位于染色体分组的contigs中，这与这4个单倍体之间的序列缺乏差异一致（表S7，图S5）。5S rDNA序列主要位于chr3和chr4上，而45S阵列位于chr5上（表S7，图S5）。在主要染色体中，rRNA编码基因主要分布在chr3（5,784个）、chr4（5,598个）和chr5（3,993个）上，5S位点的数量多于45S位点。相比之下，未锚定的染色体相关支架包含了大部分45S拷贝，特别是chr5un（39,487个rDNA注释，每个45S亚基大约有13,000个）。在chr3un和chr4un上也检测到了类似的串联阵列（表S7）。这些阵列显示出紧密连接的5.8S:18S:28S拷贝比例，支持rDNA重复单元的完整性。由于几乎所有未锚定的支架都对应于这些rDNA阵列，它们的组织进一步支持了组装的完整性。
**具有四个等位基因的全面、近乎完整的基因模型：**
我们整合了大量的转录组数据集以及同源蛋白质序列，为AAC Mountainview sainfoin基因组构建了全面的基因注释。为了最大化转录多样性并提高条件和组织特异性异构体的检测能力，使用了来自11种不同组织和条件的Illumina RNA-Seq文库，包括营养、生殖和受压样本，共44个文库（每个文库4个生物学重复），产生了大约24.2亿个链特异性150 bp配对末端读长（721.8Gb）。这些数据通过PacBio Iso-Seq全长转录本读长进行了补充，以进一步细化外显子-内含子边界并恢复可能被短读长数据遗漏的长链、可变剪接的转录本。这些数据包括大约2250万个全长非嵌合读长和116万个高质量异构体，平均转录本长度为1.7 kb，最长转录本长度为9.6 kb（1.99Gb）。
最终注释包括117,890个高置信度的蛋白质编码基因，平均编码序列长度为1,297 bp，每个基因平均有5.3个外显子（表S8）。每个单倍体总共包含大约28,700-29,000个基因，相当于每个染色体平均约4,100-4,140个基因，支持所有4个单倍体中编码序列的均匀分布（表S8）。基因集的BUSCO值非常高，完整性达到97.7%，少于0.1%的基因片段化，2.3%的基因缺失，表明保守基因内容的近乎完整表示，并证实了组装的高覆盖度和连续性（图S6）。根据严格的标准（单倍体内无重复），74.4%的直系同源群为四等位基因，11.4%为三等位基因，12.6%为二等位基因，1.7%为单等位基因（n=21,326个直系同源群）。更宽松的分析允许单倍体内的重复，得到了类似的比例（四等位基因：70.7%，三等位基因：13.4%，二等位基因：12.9%，单等位基因：3.0%；n=27,776）。这些分布表明，在四倍体基因组中大多数位点保留了所有4个等位基因。
**功能注释：**
我们结合了序列同源性、结构域预测和基于直系同源性的推断。总共72.1%到97.4%的基因与NCBI NR、TrEMBL、TAIR10和Swiss-Prot中的条目匹配，81.3%的基因获得了Clusters of Orthologous Genes（COG）分类（表S9）。总体而言，超过97%的预测基因至少被一个主要资源进行了注释，支持了基因模型的完整性和生物学相关性，并为sainfoin的下游基因表达、进化和功能研究提供了完整的基因组框架。除了蛋白质编码基因外，结构注释还识别出了1,031个假定的miRNA前体位点，代表了多个保守的植物miRNA家族（表S10）。
**sainfoin基因组中的DNA甲基化景观：**
为了表征sainfoin的表观基因组景观，我们使用Nanopore数据进行了全基因组DNA甲基化分析。我们在基因区域观察到了典型的CG甲基化（mCG）模式，转录起始和终止位点附近的甲基化水平较低，而基因体内部的甲基化水平较高（见图2a）。启动子附近区域（上游2 kb）的mCG信号较低至中等（平均每个位点的甲基化比例为45.9%），表明大多数启动子保持转录允许状态，而不是持续沉默（见图2a，图S7）。
**图2：**该图像的替代文本可能是通过AI生成的。
**sainfoin中基因和转座元件的DNA甲基化模式：**
a. 蛋白质编码基因的每个位点的胞嘧啶甲基化（mCG、mCHG和mCHH）信号，包括上游和下游2 kb的侧翼区域。
b. 三个主要转座元件（LTR Copia、LTR Gypsy和Helitron）家族的每个位点的mCG、mCHG和mCHH信号，包括1 kb的侧翼区域。
c. 蛋白质编码基因的每个位点的5-羟甲基胞嘧啶（5hmC）信号，包括2 kb的侧翼区域。
d. LTR Copia、LTR Gypsy和Helitron转座元件的每个位点的5hmC信号，包括1 kb的侧翼区域。每条线代表所有相关区域（基因或TE家族）的平滑平均值，按从开始到结束的标准化坐标对齐。5hmC和CH背景下的明显信号变化可能反映了相对于更密集的CG数据集的较低覆盖率和平滑效果。
相比之下，TE显示出均匀的高mCG水平，反映了强烈的转录沉默和表观遗传抑制。家族级分析揭示了不同的模式：LTR逆转录转座子（特别是LTR/Gypsy和LTR/Copia）在其全长上高度甲基化（平均每个位点的甲基化比例为89.5%和91.1%），而DNA转座子如Helitrons的平均mCG水平略低（平均每个位点的甲基化比例为81.2%），可能反映了不同的沉默效率或基因组背景（见图2b）。值得注意的是，着丝粒区域显示出最高的mCG水平，这与它们的异染色质结构和转录不活跃性一致（见图1d）。总体而言，这些结果突出了sainfoin中功能分区的表观基因组，其中基因丰富的区域保持中等的甲基化模式，而TE丰富的异染色质区域则高度甲基化以加强沉默效果。
**为了提供更广泛的胞嘧啶修饰视图，我们还检查了非CG甲基化背景（mCHG和mCHH）和5hmC。**mCHG和mCHH基因体的甲基化水平在平均每个位点的甲基化分数上分别达到了28.9%和24.2%的峰值，但在转录起始位点处明显减少，甲基化水平向外延伸至启动子两侧和基因体内部（图2a；图S7）。与转座元件（TEs）相关的甲基化在CHG背景下达到了55.0–58.8%的平均水平，在CHH背景下为49.0–59.7%（尽管在Helitrons中这一比例较低，为49.0%），这明显低于mCG的水平（图2b）。两个基因和TEs的平均5hmC信号都普遍较低（通常<平均每个位点的甲基化分数的10%），表明在营养组织中的生物学存在性有限或检测置信度较低（图2c，d；图s7）。 为了研究sainfoin在豆科植物中的进化关系并解释其基因组特征，我们使用了18种被子植物物种进行了比较基因组分析，其中包括14种豆科植物和4种非豆科植物作为外群。进行了正群分析，从233个广泛共享的正群构建的系统发育树得到了一个得到强支持的拓扑结构。o. viciifolia被置于hologalegina分支中，作为medicago谱系的姐妹群，与其他温带豆科植物如lotus japonicus、pisum sativum和cicer arietinum处于同一分支（图s8a）。 为了更好地理解sainfoin的进化动态，我们还分析了单倍型a基因组内的同义替换率（ks）以及选定的豆科和非豆科物种之间的同源基因的ks，以检测全基因组复制（wgd）事件。sainfoin同源基因对的ks分布在一个明显的峰值0.6附近（图s8b），这与之前在蝶形花科豆科植物中报道的全基因组复制事件一致（cannon等人，2015年）。o. viciifolia与其他豆科植物的同源基因ks比较显示出的模式与预期的分化时间相符。例如，sainfoin的近亲如medicago truncatula、m. sativa、c. arietinum和l. japonicus的ks峰值接近0.4–0.5，表明它们是在豆科植物全基因组复制之后分化的。而关系较远的豆科植物（如arachis hypogaea）则显示出稍高的ks峰值（大约0.55），这表明它们的分化时间稍早。相比之下，非豆科植物外群的ks峰值较高（>1.25），表明它们经历了更古老的分化事件和多倍体化事件（图S8b）。
除了核基因组组装外，我们还组装并环化了叶绿体和线粒体基因组。叶绿体基因组组装成一个122,397 bp的环状分子（GC=34.6%），包含78个蛋白质编码基因、24个tRNA编码基因和4个rRNA编码基因（表S11；图S9）。这些基因包括与光系统相关的PSA和PSB基因，以及推测编码核酮糖-1,5-二磷酸羧化酶/加氧酶大亚基（rbcL）、ATP合成酶（atpA-H）、RNA聚合酶（rpoA-C）、核糖体蛋白大/小亚基（RPL/RPS）、NADH脱氢酶亚基和细胞色素b6/f亚基（petB和petD；表S11）的基因。
线粒体基因组组装成一个321,160 bp的连续片段（GC=45.2%）。我们总共注释了50个基因，包括31个保守的蛋白质编码基因、16个tRNA编码基因和3个线粒体rRNA编码基因（rrn5、rns、rnl），涵盖了核心被子植物线粒体基因组（表S12；图S10）。蛋白质编码基因包括那些推测编码ATP合成酶亚基（atp1/4/6/8/9）、细胞色素b（cob）、细胞色素-c氧化酶亚基（cox1/2/3）、NADH脱氢酶（复合体I）亚基（nad1/3/4/4L/5/6/7/9）、细胞色素-c成熟因子（ccmB/C/FC/FN）、双精氨酸转运蛋白组分（tatC）和II类内含子成熟酶（matR）以及核糖体蛋白（rpl5、rpl16、rps1/3/4/12/14）的基因。此外，还存在许多额外的假想开放阅读框（ORFs；≥100个氨基酸），通常以多个拷贝的形式存在，这与植物线粒体基因组的重复富集特性一致。
组织富集的基因表达模式揭示了功能特化。为了研究与组织身份和功能相关的基因表达模式，我们使用了来自9种sainfoin组织的RNA-Seq数据进行了组织富集基因分析（图3a）。我们识别出86,084个表达基因（在至少2个文库中的计数≥1），其中33,989个基因在所有评估的发育阶段和组织中都有表达（数据S1）。在这些普遍表达的基因中，有219个基因在所有组织中的表达稳定性很高，变异系数（CV）≤0.20（数据S1）。应用更严格的截止标准（CV<0.15）后，这一数量减少到只有5个基因（被注释为ankyrin重复蛋白和DEAD-box RNA解旋酶），且没有基因满足CV<0.10的阈值（数据S1），这突显了在多组织数据集中识别普遍不变转录本的难度。
为了评估组织偏见，我们保留了FDR<0.01且在最高表达组织和第二高表达组织之间log2倍数变化至少为1的基因，从而得到了高置信度的组织特异性基因集（图S11，数据S2）。总的来说，我们的组织富集流程识别出了12,586个具有明显器官偏向性表达的基因。对这些组织富集基因集进行GO富集分析证实了它们在生物学上的连贯表达模式（数据S2）。例如，花蕾中表现出与花粉发育、细胞骨架动态和能量代谢相关的基因富集。在根组织中，与离子运输、应激信号传导和异生物质解毒相关的基因富集。茎组织中表现出与次生细胞壁生物合成相关的基因高表达，特别是木质素和木聚糖途径相关的基因。在幼叶中，与光合作用、质体功能和色素生物合成相关的基因富集，而在成熟叶中，与营养运输、衰老和水杨酸介导的反应相关的基因上调（数据S2）。
在这个框架下，我们特别研究了黄酮类生物合成途径中的基因，包括一组参与CT生物合成的核心基因。在从拟南芥和V. vinifera参考文献中筛选出的20个候选基因中，有9个在我们的图谱中表现出组织富集（图3b）。其中6个基因在幼叶中富集，包括那些推测编码花青素合成酶（ANS；也称为leucoanthocyanidin dioxygenase [LDOX]）、查尔酮合成酶（CHS）、花青素还原酶（ANR）、黄酮酮3'-羟化酶（F3H）、黄酮3', 5'-羟化酶（F3′5’H）和二氢黄酮醇还原酶（DFR）的基因。其余3个基因在根（UDP-糖基转移酶 [UGT]）和茎（O-甲基转移酶 [OMT] 以及4-香豆酸辅酶A连接酶 [4CL]）中富集（图3b）。
为了研究sainfoin对干旱和积水胁迫的转录组响应，我们还使用来自分别经历7天干旱和积水的植物叶片组织的RNA-Seq数据进行了分析，并与在正常灌溉条件下生长的相同阶段的叶片进行了比较（图4a）。干旱胁迫后，与正常灌溉对照相比，有14,188个基因模型（包括每个基因的所有等位基因）表现出差异表达，其中5,536个基因上调，8,652个基因下调。同样，在积水胁迫后，有19,040个差异表达基因（DEGs），其中8,489个上调，10,551个下调（图4b）。许多DEGs表现出较大的log2倍数变化（>10）和高度显著的调整后P值，表明它们对胁迫有强烈的转录响应（图4c；数据S3）。有趣的是，尽管干旱和积水代表不同的胁迫条件，但它们之间有相当数量的DEGs是共有的，有2,562个基因在两种处理下都上调，4,318个基因下调（图4c），表明对水相关胁迫有保守的核心响应。
为了补充GO术语分析，我们还进行了PageMan和MapMan通路分析，以进一步揭示干旱和积水条件下共享和独特的全球转录重编程。干旱和积水都诱导了与ABA、非生物胁迫、主要碳水化合物降解、细胞壁前体合成、染色质结构以及AP2/EREBP转录因子（上调基因占多数）相关的通路/家族的相似响应，以及生长素、乙烯和受体激酶信号传导（下调基因占多数）（图5a，b；图S13；图S14；数据S5）。然而，正如GO术语分析所示，某些通路/类别中也观察到了明显的胁迫特异性。例如，干旱实验中上调的DEGs在活性氧（ROS）代谢过程、过氧化氢生物合成过程的调节、氧化还原酶活性等GO术语中表现出富集。而在下调的DEGs中，干旱特异性富集的GO术语包括响应缺水反应的术语，如光系统II（PSII）修复、角质素生物合成过程的调节、苯丙素生物合成过程、类胡萝卜素生物合成过程和抗氧化活性等。此外，与植物细胞壁相关的GO术语也表现出富集，包括细胞壁组织或生物发生、植物类型次生细胞壁生物发生、次生细胞壁生物发生的调节等。另一方面，积水特异性富集的GO术语包括光系统II中的光捕获、黄酮类生物合成过程、气孔开闭调节等。
为了研究sainfoin对干旱和积水胁迫的转录组响应，我们还使用了来自分别经历7天干旱和积水的植物叶片组织的RNA-Seq数据（从茎尖数第三片叶子开始），并与在正常灌溉条件下生长的相同阶段的叶片进行了比较（图4a）。干旱胁迫后，与正常灌溉对照相比，有14,188个基因模型（包括每个基因的所有等位基因）表现出差异表达，其中5,536个上调，8,652个下调。同样，在积水胁迫后，有19,040个DEGs，其中8,489个上调，10,551个下调（图4b）。许多DEGs表现出较大的log2倍数变化（>10）和高度显著的调整后P值，表明它们对胁迫有强烈的转录响应（图4c；数据S3）。有趣的是，尽管干旱和积水代表不同的胁迫条件，但它们之间有相当数量的DEGs是共有的，有2,562个基因在两种处理下都上调，4,318个基因下调（图4c），表明对水相关胁迫有保守的核心响应。
为了研究sainfoin叶片在干旱和积水胁迫下的转录组响应，我们进行了主成分分析（PCA）（图3a）。样本按组织类型聚类，根（RT，RC）、花（FB，OF）、茎、叶（DL，EL，第一叶和第三叶）之间有明显的分离。b 黄酮类到缩合单宁的生物合成途径中，富集的酶用不同颜色突出显示。在sainfoin中鉴定为组织富集的基因编码的酶分别用绿色（幼叶）、蓝色（茎）和灰色（根）标出。PAL苯丙氨酸氨裂解酶、C4H肉桂酸-4-羟化酶、4CL 4-香豆酸：CoA连接酶、CHS查尔酮合成酶、CHI查尔酮异构酶、F3H黄酮酮3-羟化酶、F3′5’H黄酮类3’, 5’-羟化酶、DFR二氢黄酮醇还原酶、ANS花青素合成酶、ANR花青素还原酶、LAR白花青素还原酶、OMT O-甲基转移酶、UGT UDP-糖基转移酶。
为了评估组织偏见，我们保留了FDR<0.01且在最高表达组织和第二高表达组织之间log2倍数变化至少为1的基因，从而得到了高置信度的组织特异性基因集（图S11，数据S2）。总的来说，我们的组织富集流程识别出了12,586个具有明显器官偏向性表达的基因。对这些组织富集基因集进行GO富集分析证实了它们在生物学上的连贯表达模式（数据S2）。例如，花蕾中表现出与花粉发育、细胞骨架动态和能量代谢相关的基因富集。在根组织中，与离子运输、应激信号传导和异生物质解毒相关的基因富集。茎组织中表现出与次生细胞壁生物合成相关的基因高表达，特别是木质素和木聚糖途径相关的基因。与光合作用、质体功能和色素生物合成相关的基因在幼叶中富集，而在成熟叶中，与营养运输、衰老和水杨酸介导的反应相关的基因上调（数据S2）。
在这个框架下，我们特别研究了黄酮类生物合成途径中的基因，包括一组参与CT生物合成的核心基因。在从拟南芥和V. vinifera参考文献中筛选出的20个候选基因中，有9个在我们的图谱中表现出组织富集（图3b）。其中6个基因在幼叶中富集，包括那些推测编码花青素合成酶（ANS；也称为leucoanthocyanidin dioxygenase [LDOX]）、查尔酮合成酶（CHS）、花青素还原酶（ANR）、黄酮酮3'-羟化酶（F3H）、黄酮类3', 5'-羟化酶（F3′5’H）和二氢黄酮醇还原酶（DFR）的基因。其余3个基因在根（UDP-糖基转移酶 [UGT]）和茎（O-甲基转移酶 [OMT] 以及4-香豆酸辅酶A连接酶 [4CL]）中富集（图3b）。
为了研究sainfoin对干旱和积水胁迫的转录组响应，我们还使用了来自分别经历7天干旱和积水的植物叶片组织的RNA-Seq数据（从茎尖数第三片叶子开始），并与在正常灌溉条件下生长的相同阶段的叶片进行了比较（图4a）。干旱胁迫后，与正常灌溉对照相比，有14,188个基因模型（包括每个基因的所有等位基因）表现出差异表达，其中5,536个上调，8,652个下调。同样，在积水胁迫后，与正常灌溉对照相比，有19,040个DEGs，其中8,489个上调，10,551个下调（图4b）。许多DEGs表现出较大的log2倍数变化（>10）和高度显著的调整后P值，表明它们对胁迫有强烈的转录响应（图4c；数据S3）。有趣的是，尽管干旱和积水代表不同的胁迫条件，但它们之间有相当数量的DEGs是共有的，有2,562个基因在两种处理下都上调，4,318个基因下调（图4c），表明对水相关胁迫有保守的核心响应。
为了研究sainfoin对干旱和积水胁迫的转录组响应，我们还使用了来自干旱（D）和积水（W）条件下植物的叶片组织的RNA-Seq数据进行了主成分分析（PCA）（图4a）。b 火山图显示了干旱与对照（左）和积水与对照（右）之间差异表达的基因（DEGs）。log₂倍数变化>1且调整后P值<0.05的基因用红色突出显示。c 文氏图显示了干旱和积水处理之间上调和下调DEGs的重叠。大量的重叠表明这两种类型的水胁迫有共同的转录响应。
为了研究干旱和积水条件下DEGs的GO术语富集情况，我们比较了干旱和对照条件以及积水条件和对照条件之间的DEGs（图S12，数据S4）。对于上调的DEGs，干旱和积水植物的叶片在活性氧（ROS）代谢过程、过氧化氢生物合成过程的调节、氧化还原酶活性等GO术语中表现出富集。对于下调的DEGs，常见的GO术语包括光合作用和信号受体活性等。然而，许多显著富集的GO术语是每种胁迫特有的。例如，干旱实验中上调的DEGs在响应缺水反应的GO术语中表现出富集，如光系统II（PSII）修复、角质素生物合成过程的调节、苯丙素生物合成过程、类胡萝卜素生物合成过程和抗氧化活性等。此外，与植物细胞壁相关的GO术语也表现出富集，包括细胞壁组织或生物发生、植物类型次生细胞壁生物发生、次生细胞壁生物发生的调节等。另一方面，积水特异性富集的GO术语包括光系统II中的光捕获、黄酮类生物合成过程、气孔开闭调节等。CHO（碳水化合物）、OPP（氧化戊糖磷酸途径）、TCA（三羧酸循环）
等位基因特异性表达和单倍型偏差
为了研究四倍体紫花苜蓿中的等位基因表达偏差，我们使用了基于单倍型解析的基因模型以及来自在对照和胁迫（干旱和积水）条件下生长的植物的叶片组织的RNA-Seq数据进行了等位基因特异性表达分析。在代表所有4种单倍型（A–D）的12,424个高置信度同源基因组中，大多数表现出轻微的表达偏差，其分布偏向于较低的偏差值（偏差<0.25，定义为表达最高和最低的同源基因之间的差异），这与单倍型水平上的广泛平衡表达一致（图6a）。仅有2.6–3.0%的基因组对（320–370个基因模型）显示出强烈的表达优势（偏差>0.5，表明单个等位基因的功能优势），这表明在紫花苜蓿中固有的顺式调控优势相对较少。
图6
该图像的替代文本可能是使用AI生成的。
全尺寸图像
在对照和胁迫条件下紫花苜蓿中的等位基因特异性表达偏差分布。a. 对照（蓝色）、积水（橙色）和干旱（绿色）条件下，1:1:1:1直系群中等位基因表达偏差（最大值-最小值单倍型比率）的累积分布函数。红色虚线表示轻微优势阈值（0.25）。紫色虚线表示强烈优势阈值（0.5）。b. 散点图显示了等位基因偏差（x轴）与胁迫和对照之间最高等位基因比率变化（Δ，y轴）的关系。虚线表示阈值（偏差>0.5，Δ≥0.2）。红色点表示同时具有强烈偏差并发生等级顺序变化的基因组对（Δ≥0.2），而灰色点表示其他所有情况。
为了评估胁迫下的调控动态，我们研究了在胁迫条件下最高表达同源基因发生等级顺序变化的事件。值得注意的是，超过三分之一的基因组对（干旱条件下4,435个，积水条件下4,774个）表现出这种变化，表明存在广泛但微妙的调控变化。实际上，93.2%的干旱相关和91.2%的积水相关基因组对的表达比率变化Δ<0.1（即胁迫和对照之间最高表达等位基因比率的差异），这表明是细调而非主要的等位基因变化。在应用明确的效应大小过滤器后，只剩下一个小部分基因组对（积水条件下420个，Δ≥0.1；干旱条件下104个，Δ≥0.2）。值得注意的是，具有强烈表达偏差的基因与发生优势变化的基因之间的重叠非常有限（每种条件下分别为60个和26个），这意味着高度偏倚的等位基因通常不受环境扰动的影响（图6b）。GO术语分析也支持这一观点：等级顺序变化组（Δ≥0.1/0.2）在两种胁迫条件下均未显示出富集的术语（FDR<0.05），而强偏差组则产生了一些与胁迫相关的类别，包括单加氧酶活性、次级代谢物生物合成过程以及氧化还原酶活性（干旱与对照比较），以及细胞壁修饰（积水与对照比较）（表S13）。这些发现共同表明，紫花苜蓿采用了一种灵活而平衡的调控系统，其中大多数基因在轻微的单倍型偏差下运作，胁迫引发的是细调而非广泛的优势行为，从而增强了表达景观的总体稳定性。
讨论
尽管紫花苜蓿作为饲料作物具有优势，但其改良速度相对较慢，部分原因在于相关基因组资源的匮乏。虽然最近对一种“普通”紫花苜蓿的栽培品种进行了测序（He等人，2024年），但大多数现代紫花苜蓿栽培品种同时包含“普通”和“巨型”类型背景（Mora-Ortiz和Smith，2018年），因此我们旨在为更具代表性的紫花苜蓿栽培品种生成一个全面的基因组和转录组，以进一步支持基础生物学研究和实际育种应用。为此，我们使用PacBio和Oxford Nanopore长读长测序技术、Illumina短读长测序技术以及Hi-C数据进行基因组组装，生成了一个高度连续且完整的2.36 Gb基因组，该基因组被划分为28个假染色体，代表了7条基因组的4个同源基因（图1，表1）。
4种单倍型之间的全基因组比对显示了广泛的共线性，结构重排较少（图S2），表明其遗传性质较为平衡。有趣的是，AAC Mountainview基因组比之前发表的未知中国紫花苜蓿栽培品种的1.95 Gb基因组（He等人，2024年）要大得多，我们的基因组组装还在2号、5号和6号染色体上显示了一些其他显著差异，包括插入缺失、易位和倒位（图S3）。这些差异可能是由于栽培品种特异性、重复序列处理策略和/或组装方法的差异所致。系统基因组分析将紫花苜蓿归入Hologalegina分支，并表明它经历了大约5800万年前发生的蝶形花科范围内的全基因组复制（WGD），这与先前的研究结果一致（He等人，2024年；图S8）。我们的发现与其他豆科植物的基因组进化一致，在蝶形花科范围内的WGD之后，染色体共线性通常得以保留（Cannon等人，2015年）。
重复序列分析表明，基因组的64%以上由重复序列组成，其中LTR逆转录转座子所占比例最大（表2），这与之前在紫花苜蓿（He等人，2024年）以及其他温带多倍体豆科植物（如紫花苜蓿Shen等人，2020年）中的发现一致。在这些重复序列中，Gypsy（11.6%）和Copia（5.4%）元素在紫花苜蓿中占主导地位，21.7%的基因组由未分类的LTRs组成（表2），这可能反映了高度分化或特定谱系的逆转录转座子，这些逆转录转座子未被现有的重复序列库捕获。此外，当前的紫花苜蓿基因组组装具有14.7的LAI（基因组连续性指数），单个染色体的LAI得分在10.7到24.3之间，所有4种单倍型的LAI得分均超过了参考质量阈值10（Ou等人，2018年；图S4；表S6）。这一基因组连续性指标明显高于最近发表的“普通”紫花苜蓿基因组，后者的LAI得分在每个单倍型之间为4.8到7.1（He等人，2024年）。这可能反映了包含长读长PacBio HiFi测序的优势，有助于恢复完整的逆转录转座子，说明了测序技术对重复序列表示的影响（Pucker等人，2022年）。
我们紫花苜蓿基因组组装的一个意外但生物学上有趣的结果是在未锚定的染色体相关支架上识别出了大量高拷贝数的串联rRNA编码基因阵列（图S5；表S7）。即使有大量的长读长数据和Hi-C支架，大多数45S rDNA序列也无法被锚定，这可能是因为45S rDNA位点通常位于染色体末端区域（Roa和Guerra，2012年）。尽管这些阵列未锚定，但它们的结构和拷贝数——特别是chr5un上的大约13,000个拷贝单元，总计近40,000个rDNA单元，且5.8S:18S:28S的比例平衡——表明它们很可能是功能性核仁组织区（NORs），而不是碎片化的组装产物。虽然不同植物物种/品系之间的rDNA拷贝数可能存在显著差异（Hasterok等人，2006年），但即使是高质量的长读长组装也经常在rDNA阵列处断裂（Igolkina等人，2025年），这突显了它们的生物学规模以及解析它们的技术难度。虽然大多数其他豆科植物主要在细胞遗传学水平上评估了rDNA的组织结构（例如Abirached-Darmency等人，2005年），但已知异常大的阵列会影响其他植物的调控，包括Brachypodium hybridum中的核仁优势（Borowska-Zuchowska等人，2021年）和小麦中的NOR活性（Tulpová等人，2022年）。因此，紫花苜蓿中观察到的这些大规模阵列为研究饲料豆科植物中的rDNA组织和潜在调控效应提供了一个有价值的参考点。
基因注释产生了117,890个高置信度的蛋白质编码基因（表1），每种单倍型总共包含大约28,700–29,000个基因，BUSCO得分为97.7%（Fabales），超过97%的基因在多个数据库中得到了成功注释（表S4；图S6）。作为对比，最近的组装报告了M. truncatula中的44,623个单倍体等效基因（Pecrix等人，2018年）、P. sativum中的47,526个基因（Yang等人，2022年）、L. japonicus中的28,251个基因（Li等人，2020年）以及M. sativa中的大约41,000个基因（Chen等人，2020年）。相比之下，在之前的“普通”四倍体紫花苜蓿组装中鉴定出了109,998个高置信度基因，每种单倍型平均包含27,284个基因（He等人，2024年）。这种紫花苜蓿基因组之间的差异可能反映了生物学因素和/或注释技术方面的差异，我们研究中使用的注释流程可能更广泛地捕获了基因空间。此外，结合长读长Iso-Seq数据对我们的注释过程至关重要，因为它使得能够恢复近10,000个转录本和基因异构体，而这些是仅靠短读长RNA-Seq证据无法发现的。除了蛋白质编码基因外，AAC Mountainview紫花苜蓿基因组还包含1,031个假定的miRNA前体位点，代表了多个保守的植物miRNA家族（表S10；Guo等人，2020年），为未来的比较和功能分析提供了有用的资源。
在植物中，胞嘧啶甲基化发生在三种序列背景下，包括CG、CHG和CHH（H=A、T或C），并支持基因组的稳定性、转座元素的沉默和基因调控。其中，mCG通常是体细胞组织中最丰富且最稳定的标记（Law和Jacobsen，2010年）。在紫花苜蓿中，我们发现着丝粒区域显示出最高的mCG水平（图1d），这与它们的异染色质性质相符（Meyer，2010年）。此外，基因体中的mCG含量较高（非CG甲基化水平较低）（图2a；图S7），这可能在表达稳定、剪接准确性、防止转座元素插入和/或减少异常转录方面发挥作用，尽管这种甲基化的确切作用仍有争议（Muyle等人，2022年）。启动子附近区域（转录起始位点上游2 kb）的mCG水平较低（图2a），这与通常促进转录起始的宽松染色质状态一致（Qiao等人，2025年）。相比之下，转座元素，尤其是LTR逆转录转座子，显示出密集的mCG、mCHG和mCHH（图2b），表明强烈的表观遗传沉默（Liu和Zhao，2023年）。相比之下，DNA转座元素（如Helitrons）的mCG水平较低，这可能表明它们的转录活性较低或较新（例如Vonholdt等人，2012年）。在大豆中也观察到了类似的模式，其中mCG在基因体中占主导地位，而非CG甲基化（特别是mCHH）在转座元素和异染色质区域中富集（例如Song等人，2013年）。除了mCG、mCHG和mCHH之外，我们还研究了5hmC的存在，这是一种氧化形式的胞嘧啶甲基化；其在植物中的存在和生物学功能尚不完全清楚（例如Erdmann等人，2014年；Wang等人，2015b）。正如预期的那样，紫花苜蓿基因组中的5hmC水平较低（<平均每个位点的10%；图2c，d；图s7），这表明在幼叶中的普遍性较低，与其他植物物种的先前报告一致（例如wang等人，2015b）。
除了核基因组外，我们还组装了aac mountainview紫花苜蓿的完整叶绿体和线粒体基因组。122.4 kb的叶绿体基因组缺乏大型倒位重复序列（ir），因此没有表现出许多植物物种中典型的四部分结构：大型单拷贝（lsc）-ir-小型单拷贝（ssc）-ir（wicke等人，2011年；图s9）。这种结构与缺乏ir的豆科分支（irlc）以及onobrychis的质体一致（moghaddam等人，2022年）。该叶绿体基因组包含78个蛋白质编码基因、24个trna编码基因和4个rrna编码基因，gc含量为34.6%，这与之前的紫花苜蓿研究结果相符（jin等人，2021年）。这包括psi和psii的全部基因、rubisco、atp合成酶、核糖体蛋白和nadh脱氢酶复合体的完整基因组（表s11），所有这些都具有典型的被子植物质体内含子结构（daniell等人，2016年）。
另一方面，线粒体基因组被解析为一个321.2 kb的连续体（图s10），其大小介于其他已测序的蝶形花科豆科植物的线粒体基因组之间（choi等人，2019年）。该基因组的gc含量为45.2%，共包含50个基因，其中包括31个蛋白质编码基因、16个trna编码基因和多个假设的开放阅读框（≥100个氨基酸），这与植物线粒体基因组的重复序列丰富特性一致（wynn和christensen，2019年）。蛋白质编码基因包括那些编码呼吸链复合体和核糖体蛋白的基因（表s12），这与其他豆科植物的观察结果相似（choi等人，2019年）。 除了核基因组外，我们还组装了aac mountainview紫花苜蓿的完整叶绿体和线粒体基因组。122.4 kb的叶绿体基因组缺乏大型倒位重复序列（ir），因此没有表现出许多植物物种中典型的四部分结构：大型单拷贝（lsc）-ir-小型单拷贝（ssc）-ir（wicke等人，2011年；图s9）。这种结构与缺乏ir的豆科分支（irlc）以及onobrychis的质体一致（moghaddam等人，2022年）。该叶绿体基因组包含78个蛋白质编码基因、24个trna编码基因和4个rrna编码基因，gc含量为34.6%，这与之前的紫花苜蓿研究结果相符（jin等人，2021年）。这包括psi和psii的全部基因、rubisco、atp合成酶、核糖体蛋白和nadh脱氢酶复合体的完整基因组（表s11），所有这些都具有典型的被子植物质体内含子结构（daniell等人，2016年）。 另一方面，线粒体基因组被解析为一个321.2>
除了核基因组外，我们还组装了aac mountainview紫花苜蓿的完整叶绿体和线粒体基因组。122.4 kb的叶绿体基因组缺乏大型倒位重复序列（ir），因此没有表现出许多植物物种中典型的四部分结构：大型单拷贝（lsc）-ir-小型单拷贝（ssc）-ir（wicke等人，2011年；图s9）。这种结构与缺乏ir的豆科分支（irlc）以及onobrychis的质体一致（moghaddam等人，2022年）。该叶绿体基因组包含78个蛋白质编码基因、24个trna编码基因和4个rrna编码基因，gc含量为34.6%，这与之前的紫花苜蓿研究结果相符（jin等人，2021年）。这包括psi和psii的全部基因、rubisco、atp合成酶、核糖体蛋白和nadh脱氢酶复合体的完整基因组（表s11），所有这些都具有典型的被子植物质体内含子结构（daniell等人，2016年）。
另一方面，线粒体基因组被解析为一个321.2 kb的连续体（图s10），其大小介于其他已测序的蝶形花科豆科植物的线粒体基因组之间（choi等人，2019年）。该基因组的gc含量为45.2%，共包含50个基因，其中包括31个蛋白质编码基因、16个trna编码基因和多个假设的开放阅读框（≥100个氨基酸），这与植物线粒体基因组的重复序列丰富特性一致（wynn和christensen，2019年）。蛋白质编码基因包括那些编码呼吸链复合体和核糖体蛋白的基因（表s12），这与其他豆科植物的观察结果相似（choi等人，2019年）。>据我们所知，这是首次对山萝卜（sainfoin）的线粒体基因组进行特征分析，为该物种的后续基因组研究提供了额外信息。除了基因组资源外，我们还从9个不同的发育阶段/组织类型中生成了大量转录组数据（图3a），这些数据为后续实验提供了生物学背景。在多个组织中相对稳定表达的一小部分基因（CV≤0.20的基因有219个；CV<0.15的基因有5个；CV<0.10的基因没有）可以作为qRT-PCR的参考基因（数据S1），但严格阈值导致的基因数量急剧减少表明，在不同植物组织中普遍存在的管家基因控制因子较为罕见。这表明在特定的实验设计中可能需要多基因标准化。同时，我们的组织富集流程识别出12,586个具有明显器官偏向性表达的基因，这些基因的功能富集与最高表达组织高度一致（图S11；数据S2），例如花蕾中的花粉发育和细胞骨架重塑、根部的离子转运和应激信号传导、茎部的木质素/木聚糖生物合成以及幼叶中的光合作用和色素代谢（数据S2）。在这个框架下，我们还发现大多数核心类黄酮（CT）生物合成基因在叶片中表达偏重，包括编码CHS、F3H、F3′5’H、DFR、ANS和ANR的基因（图3b）。这种分布与山萝卜幼叶中类黄酮的优先积累现象相符（Li等人，2014年），并且这些类黄酮被认为具有抵御昆虫害虫、氧化应激和植物病原体的功能（Dixon和Sarnala，2020年）。有趣的是，两个参与将花青素前体转化为花青素的基因（OMT和UGT）分别在茎部和根部表达偏重（图3b），这也可能有助于类黄酮前体在山萝卜叶片中的优先转化。为了将这些发现置于基因组背景下，我们还比较了AAC Mountainview品系与之前发表的山萝卜基因组（He等人，2024年）中黄酮生物合成基因家族的拷贝数。这些家族中的绝大多数在两个基因组中的假定同源基因数量大致相似（表S5），表明黄酮生物合成基因库在这两个基因组之间具有高度保守性。基因拷贝数的微小差异可能反映了基因注释、组装分辨率或常见于植物基因家族的小规模复制和丢失事件（Panchy等人，2016年）。另一方面，CHI家族在AAC Mountainview基因组中有20个同源基因，而在之前测序的山萝卜基因组中只有12个同源基因（He等人，2024年；表S5），这可能反映了两种基因型之间的拷贝数差异较大。然而，需要进一步研究来确定这一发现的生物学意义。总体而言，这些结果证实了之前在山萝卜中检测到的类黄酮生物合成途径的扩展（He等人，2024年），并增加了跨组织的解析能力。
由于山萝卜的农艺表现受到干旱和水涝等非生物胁迫的制约（例如，Heinrichs，1970年；Biligetu等人，2014年），提高其对这些环境挑战的耐受性仍然是关键的育种目标，我们还对来自干旱、水涝和正常灌溉条件下的山萝卜叶片进行了RNA-Seq分析（图4a，b）。山萝卜叶片在干旱和水涝条件下的转录响应有显著重叠，有超过6,800个基因模型（等位基因）在这两种条件下都存在（图4c；数据S3）。总体而言，山萝卜在ABA相关信号通路中发生了变化，与ROS生物合成、碳水化合物降解、细胞壁前体以及AP2/EREBP转录因子表达增加相关基因的表达增加，而与光合作用、受体激酶信号传导和生长素/乙烯通路相关的基因表达减少（图5a，b；图S12；图S13；图S14；数据S4；数据S5）。这些变化之上还有一组更小的差异表达基因（DEGs），这些基因仅在干旱或水涝胁迫条件下特异性表达（图4c；数据S3）。在干旱情况下，这包括参与细胞壁强化的上调基因、角质层生成以及抗氧化剂及其活性的增加，以及参与气孔开放、黄酮生物合成、营养运输、钙和MAP激酶信号传导的基因下调（图5a；图S12；图S13；图S14；数据S4；数据S5）。另一方面，水涝导致参与水杨酸通路、离子和氨基酸运输的基因以及许多DOF、HB、MYB和bZIP转录因子的表达增加，而参与能量产生、细胞壁代谢和核糖体生物合成的基因表达减少（图5b；图S12；图S13；图S14；数据S4；数据S5）。虽然许多这些变化与干旱和水分胁迫下植物的常见转录响应一致（例如，Mukarram等人，2021年；Pan等人，2021年），但其他一些变化，如信号传导的减少（例如，Schulz等人，2021年）、黄酮生物合成的减少（例如，Baozhu等人，2022年）、能量相关通路的改变（例如，Zahra等人，2021年）以及可能的ROS过度积累（例如，Fujita和Hasanuzzaman，2022年），为未来提高山萝卜的非生物胁迫耐受性提供了潜在的目标。尽管这些发现为未来旨在提高山萝卜耐水胁迫能力的育种工作提供了分子框架，但将茎部和根部的转录组在各个时间点的整合、小RNA分析、代谢组数据以及植物生长/形态变化结合起来，对于将分子状态与决定山萝卜在不同水分条件下表现的生理适应联系起来至关重要。
山萝卜通常被认为是一种自四倍体（He等人，2024年）；然而，这一点存在争议，也有观点认为它可能是一种由非常接近的祖先物种衍生的异源四倍体（Özturk和Tek，2024年）。虽然异源多倍体基因组可能表现出强烈的亚基因组优势和同源基因表达偏倚，但这在自多倍体中通常不观察到的（Alger和Edger，2020年）。与山萝卜的自四倍体解释一致（同源单倍型而非不同的亚基因组），我们的等位基因特异性表达分析显示，山萝卜保持了相对平衡的表达谱，大多数同源基因四联体在正常灌溉、干旱和水涝胁迫条件下仅表现出轻微到中等的表达偏倚（图6）。然而，一小部分等位基因在对照和胁迫条件下表现出强烈的表达偏倚（大约3-4%的等位基因）（图6），在干旱条件下生物上合理的富集包括单加氧酶/氧化还原酶活性和次级代谢物生物合成，在水涝条件下包括细胞壁修饰（表S13），这表明只有少数位点具有显著的顺式不对称性。在其他自多倍体中也报告了类似的结果（例如，马铃薯、甘蔗、苜蓿；Pham等人，2017年；Chen等人，2020年；Xue等人，2022年），这与自多倍体表达的剂量平衡/化学计量限制一致。在这样的系统中，优先等位基因表达与高水平的基因组异质性有关（Pham等人，2017年），可能反映了微妙的、位点特异性的顺式差异或局部染色质特征，而不是广泛的亚基因组效应。实际上，有证据表明这种现象可能与等位基因特异性的DNA甲基化、TE插入的变异、局部染色质结构和/或顺式调控元件的变化有关（Springer和Stupar，2007年；Xue等人，2022年）。总的来说，我们的结果表明山萝卜具有自四倍体起源，并伴有适度的单倍型差异，但没有新异源多倍体典型的系统范围表达优势。然而，目前还不能完全排除它是由密切相关的祖先物种衍生的异源四倍体的可能性。
总之，我们为AAC Mountainview品系生成了一个分阶段的、染色体级别的参考核基因组，以及全面的基因模型和甲基化图谱。与之前测序的未公开山萝卜品系的基因组（He等人，2024年）相比，该基因组提供了更高的单倍型分辨率和重复序列连续性，以及扩展的基因注释。此外，我们还为AAC Mountainview品系生成了全长叶绿体基因组序列，以及该物种的第一个线粒体基因组序列。我们还提供了一个组织富集图谱，突出了幼叶中类黄酮生物合成等通路，并阐明了干旱和水涝胁迫下的保守性和特异性转录组变化。此外，我们发现山萝卜在等位基因之间的转录调控是平衡的，这表明它具有自四倍体的特性。这些数据为后续的山萝卜育种工作提供了坚实的基础，将有助于改善这种有前景的饲料作物的农艺性能。
打赏
