光荚含羞草基因组-文献精读26

Haplotype-resolved genome of Mimosa bimucronata revealed insights into leaf movement and nitrogen fixation

光荚含羞草单倍型解析基因组揭示了叶片运动和固氮的相关机制

摘要

背景

光荚含羞草起源于热带美洲，具有独特的叶片运动特征，其运动速度相对较慢。此外，这种植物还具有固氮能力。尽管这些特性引人入胜，但由于缺乏光荚含羞草的基因组资源，相关综合研究一直受阻。

结果

为了揭示叶片运动和固氮的复杂机制，我们成功组装了一个高质量的、单倍型解析的参考基因组，基因组总长648 Mb，并锚定在13条拟染色体上。共注释了32,146个蛋白质编码基因。其中，单倍型A注释了31,035个蛋白质编码基因，单倍型B注释了31,440个蛋白质编码基因。结构变异（SVs）和等位基因特异性表达（ASE）分析揭示了结构变异在光荚含羞草叶片运动和固氮中的潜在作用。我们检测到了两个全基因组复制（WGD）事件，分别发生在约290万年前和约7350万年前。转录组和共表达网络分析显示，水通道蛋白（AQPs）和钙离子（Ca2+）相关离子通道基因参与了叶片运动。此外，我们还鉴定了与根瘤形成相关的基因，并分析了在共生固氮（SNF）过程中关键基因NIN的结构和演化。

结论

详细的比较基因组和转录组分析提供了对光荚含羞草叶片运动和固氮机制的深入理解。本研究提供了基因组资源，并为光荚含羞草及其他豆科植物的功能基因组学研究提供了重要参考。

背景

含羞草属是豆科家族中最多样化的属之一，拥有超过530个物种，主要分布在新热带区，部分物种分布在马达加斯加，少数分布在东非和南亚[1, 2]。该属起源于中美洲和南美洲，大多数物种分布在低中海拔的热带地区，有些物种也分布在亚热带或温暖的温带地区[3, 4]。该属有两个多样性中心，一个位于巴西中部地区，特别是物种多样性高的塞拉多地区，另一个位于墨西哥中部和南部地区[5, 6]。该属的许多物种在经济上具有重要意义，可用于建筑和木工、生产纤维素和燃料[7]，以及由于其固氮能力，用于土壤改良和退化土地恢复[8]。含羞草属的物种对不同环境具有广泛的适应性，导致不同的生长形态和光需求，从而带来了极大的物种多样性[2]。因此，含羞草属长期以来被用作生物地理和宏观进化研究的研究系统。该属特别以在易受火灾影响的环境（如巴西塞拉多）中进化出频繁适应火灾的能力而著称，并经常在湿热带和干热带生物群落之间交替[10]。含羞草属最为人所知的是其地震感应性叶片运动的特性，并且地震感应性在含羞草属的八个谱系中独立发生[2]。例如，含羞草（Mimosa pudica）属于P谱系，对触摸有快速叶片运动反应，在几秒钟内即可完成。然而，这一特性在许多其他含羞草属物种中也存在，但反应程度较低。例如，光荚含羞草（Mimosa bimucronata）属于D谱系，其叶片在几分钟内对触摸产生缓慢反应[2]。地震感应性叶片运动的生物学意义仍未完全理解。最近的一项研究表明，这确实可以作为一种防御机制来对抗草食动物[11]，但是否作为一种保护机制来保护叶片在湿润气候地区不受损害还需要进一步研究。一项研究发现，含羞草属物种可以被视为地震前兆[12]。

光荚含羞草在引入中国后，已成为一种侵略性强的外来物种[13]。另外两种物种，包括含羞草（Mimosa pudica）和刺轴含羞草（Mimosa pigra），是分布最广的含羞草属物种，其中刺轴含羞草由于其众多有利于入侵的特性，是世界上最严重的入侵物种之一[14]。含羞草是分布最广的物种，以其快速叶片运动为特征，而光荚含羞草则以其作为护坡和堤防植物的有用性、快速生长和耐水淹特性而著称。光荚含羞草属于豆科的苏木亚科[15]。光荚含羞草被称为“篱笆王”，以其对环境的强适应性以及强抗逆性和快速生长的特性而闻名。光荚含羞草还具有一定的观赏价值，具有密集的球形花序和白色、芳香的花朵。相比含羞草，光荚含羞草的叶片对外部刺激的敏感性较低，这对研究叶片运动的分子机制有利，因为在触摸前后有足够的时间采集样本。此外，光荚含羞草和含羞草位于相对较远的谱系（分别为D和F谱系）[2]，对外部刺激的叶片运动反应速率不同，这表明这两种物种在叶片运动特性上具有不同的基因组基础。因此，需要获得更多关于光荚含羞草基因组的信息，以从进化角度分析这一特性。染色体数是物种分类和基因组组成的重要特征。先前的研究表明，x = 13是含羞草属的基本染色体数，为倍性确定提供了基础，并为组装基因组的锚定提供了指导[16]。

植物运动是一种有趣的现象，通常比较缓慢。然而，一些植物物种进化出了展示快速运动的能力，可与动物的运动相媲美，例如捕蝇草（Dionaea muscipula）和含羞草[17]。膨压是驱动许多植物运动的主要动力，弹性力是次要动力。通常，离子通道、钙离子（Ca2+）和细胞内外的水分是影响膨压的主要因素，参与植物运动[18, 19]。位于质膜和液泡膜上的水通道蛋白（AQPs）在含羞草的地震感应性叶片运动中发挥作用[20]。钙调蛋白样蛋白（CML）和木葡聚糖内转糖基酶/水解酶（XTH）基因对触摸和黑暗刺激作出反应，这些刺激具有部分重叠的信号转导途径[21]。与捕食者相关的机械和化学刺激可以诱导食虫植物毛毡苔叶片中钙信号的动态变化[22]。钙离子传感器将钙信号传递到下游，触发级联反应，调节植物的生长发育及其对环境的反应[23]。由于膨压导致的快速变形只是导致植物叶片运动的步骤之一。还需要快速检测和传递信号的系统。在某种程度上，这些研究提供了叶片运动的化学和分子基础。然而，地震感应性叶片运动在含羞草属中演化的频率仍然未知，而基因组学是研究植物进化的重要而有效的方法。

在进化过程中，根瘤菌和豆科植物建立了一种独特的互惠共生关系，其特征是共生固氮（SNF），最早的共同祖先可追溯到9000万年前[24]。SNF是豆科家族大多数物种的另一个重要现象，涉及复杂的机制和相互作用[25, 26]。通过各种正向和反向遗传学方法，发现了近200个与SNF相关的基因[27]。植物中与SNF相关的基因分为两类，包括共生基因（Sym基因）和根瘤素基因[28]。Sym基因主要参与响应细菌信号、共生信号转导和根瘤感染过程，而根瘤素基因主要参与根瘤器官形成、发育和固氮过程[24, 27]。在模式豆科植物如蒺藜苜蓿（Medicago truncatula）和光叶百脉根（Lotus japonicus）中，已经鉴定了超过50个Sym基因[27]。第一个共同的共生基因，即LjSYMRK/MtDMI2，是根瘤菌和菌根共生过程所必需的[29]。另外两个研究较多的共生基因是结瘤因子受体1（LjNFR1/MtLYK3）和结瘤因子受体5（LjNFR5/MtNFP），它们与LjSYMRK/MtDMI2合作启动根瘤器官形成和细菌感染[27, 30]。此外，LjNFR5可以通过磷酸化激活Nodule Inception（NIN）基因[31]。关键的根瘤基因，即LjNIN、MtRPG和LjNFR1/MtLYK3在非根瘤形成物种中已丧失，表明这三个基因在根瘤形成和固氮过程中发挥了关键作用，并为其他豆科植物通过根瘤固氮提供了参考[25, 32]。最近的一项研究表明，非根瘤形成物种中，LjNFR5/MtNFP和NIN基因已假基因化，表明这两个基因已丧失其原有功能，并在根瘤固氮物种中形成了新功能[33]。然而，光荚含羞草根瘤中的固氮尚未系统研究，特别是参与共生和固氮过程的关键基因NIN。本项目的目标是测序光荚含羞草的基因组，并揭示该物种中叶片运动和固氮的基因组基础。

结果

基因组测序、组装和注释

光荚含羞草是二倍体，x = 13（补充图1），其基因组大小约为660 Mb，通过流式细胞术（FCM）估计（补充表1，补充图2），并通过基于Illumina数据的K-mer分析估计约为654.60 Mb（补充表2, 3和4，补充图3）。随后，构建并测序了一个PacBio HiFi文库，生成了17.34 Gb的长读序列，覆盖率为26×（补充表4）。此外，构建并测序了一个Hi-C文库，产生了59.27 Gb的干净数据（补充表5和6）。组装后的基因组大小约为648 Mb，contig N50为19 Mb，scaffold N50为48 Mb。此组装作为后续分析的参考基因组（表1和补充表7）。生成了染色体级组装，13条拟染色体锚定了基因组的627 Mb（96.66%）（图1，补充图4，补充表8）。在染色体级别上解析的单倍型组装分别生成了contig N50为6.57 Mb和6.88 Mb的单倍型A和单倍型B（表1）。

基因组特征概况从外到内的轨道：（a）13条染色体（Mbchr01-Mbchr13），（b）基因密度，（c）基因表达，紫色表示高表达水平，黄色表示低表达水平。表达数据来自六个样本，包括根、顶端分生组织（SAM）、茎和三个状态的叶片，每个样本都有三个生物重复。（d）GC含量，（e）DNA转座子元素丰度，（f）LTR/Copia丰度，（g）LTR/Gypsy丰度。所有分布图的窗口大小为1 Mb。

Illumina序列（约41×）映射到组装的基因组，映射率为96.47%，覆盖率为99.99%（补充表9）。转录组序列映射到基因组以评估序列完整性，准确度达到99.99%（补充表10）。通过基准测试单拷贝直系同源基因（BUSCO，v3.0.2）评估基因组的完整性。在1614个植物特异性直系同源基因中，识别出98.70%（1593个）为完整基因，其中94.40%（1523个）为单拷贝基因（表1，补充表11）。基因组显示了长末端重复逆转录转座子（LTR）组装指数（LAI）值为15.27（10 ≤ LAI < 20）（表1），达到参考基因组的质量。对单倍型A和单倍型B的基因组进行BUSCO分析，分别获得97.4%和98.4%的完整性（表1，补充表12）。单倍型A和B的LAI值分别为13.54和13.40（表1）。两个单倍型基因组的共识质量值（QV）均超过40（表1）。

通过同源基因预测、de novo预测和基于转录组的方法，预测并注释了32,146个蛋白质编码基因。其中，31,137个（96.86%）基因在不同功能数据库中获得了功能注释（补充表13）。注释的BUSCO完整性为98.40%（补充表11）。我们进一步注释了非编码RNA基因，单倍体基因组中共注释了165个tRNA、99个miRNA、39个rRNA和336个snRNA（补充表14）。对于单倍型A和B的基因组，分别注释了31,035和31,440个基因，注释的BUSCO完整性分别为95.8%和97.4%（表1）。

对重复序列进行了全面特征分析，发现光荚含羞草组装的单倍体基因组中57.04%由重复元素组成。值得注意的是，长末端重复逆转录转座子（LTRs）是最丰富的，占基因组组装的22.99%（补充表15）。重复序列分别占单倍型A和B基因组的58.27%和58.57%（补充表16）。

基因组进化和系统发育分析

将光荚含羞草的基因组与其他七个已测序的植物基因组进行了比较，包括四个豆科植物基因组（费氏金合欢、尖尾叶、降香黄檀和蒺藜苜蓿），以及拟南芥、番木瓜和水稻的基因组（补充表17）。选择这七个物种的原因是四个物种与光荚含羞草有亲缘关系，另外三个物种用于校准。共识别出920个单拷贝直系同源基因，并用于构建物种系统发育树（图2a，补充表18）。光荚含羞草与费氏金合欢和尖尾叶较为接近，分化时间分别为30.6百万年前（mya）和52.1百万年前（图2a），这主要是因为这三个物种都属于苏木亚科，而其他两个豆科物种属于蝶形花亚科。在五个豆科植物物种中，共享了11,546个基因家族，光荚含羞草基因组中特有252个基因，分属于70个基因家族（图2b，补充表19和20）。这些特有基因分布在光荚含羞草的所有13条染色体上（补充图5）。对这些特有基因进行了基因本体（GO）和京都基因与基因组百科全书（KEGG）富集分析。特有基因的功能主要富集在生物过程的GO术语中，如质体醌生物合成过程、内肽酶活性的负调控和肽酶活性的负调控。在KEGG分析中，这些特有基因富集在核质转运通路（补充图6）。

系统发育和基因家族分析光荚含羞草及其他代表性植物基因组的系统发育和基因家族分析 a 八个物种的分化时间估计。蓝色条的长度表示分化时间范围。内部节点上的数字显示中位分化时间（百万年前，Mya），并附有95%置信区间（CI）。一个最近的全基因组重复（WGD）事件发生在6.8百万年前的光荚含羞草和费氏金合欢中，用红色标示。

b Upset图显示了光荚含羞草与其他四个物种之间共享基因家族的分布。比较了光荚含羞草、费氏金合欢、尖尾叶、降香黄檀和蒺藜苜蓿之间的基因家族补体。共有11,546个共同基因家族和70个光荚含羞草特有的基因家族。

c 费氏金合欢、尖尾叶和光荚含羞草基因组的全基因组重复（WGD）分析。分析了光荚含羞草、尖尾叶和费氏金合欢内部平行基因的同义置换率（Ks）分布以及光荚含羞草与其他两个基因组之间的直系同源基因。直系同源基因的Ks分布中的峰值暗示了物种分化的潜在时间，而平行基因的峰值则指示了全基因组重复（WGD）事件。

d 尖尾叶和光荚含羞草的共线性显示了它们共享区域的范围。灰线连接匹配的基因对。

基因家族扩展/收缩分析

我们检查了系统发育分析中包括的八个物种的基因家族扩展和收缩。在光荚含羞草基因组中，1,817个基因家族表现出扩展，而1,078个基因家族表现出收缩。此外，222个基因家族经历了快速进化（补充图7和补充表21）。GO和KEGG富集分析进一步揭示了快速进化基因家族的功能意义。在GO分析中，快速进化基因家族主要与萜类代谢过程和蛋白质磷酸化相关。在KEGG富集分析中，快速进化基因家族主要在倍半萜和三萜生物合成以及谷胱甘肽生物合成途径中进化（补充图8）。

全基因组重复（WGD）事件分析和共线性

全基因组重复（WGD）事件对植物基因组进化具有重大影响[25, 34]。我们在光荚含羞草、费氏金合欢和尖尾叶中调查了WGD事件（图2c）。光荚含羞草和费氏金合欢属于含羞草类群，并共享一个古老的WGD事件。光荚含羞草中的Ks值分布显示出两个峰值，分别在Ks值约为0.03和0.76。第一个峰值0.03表明最近的WGD事件发生在约290万年前（Mya），在光荚含羞草和费氏金合欢分化之后，晚于费氏金合欢的最近WGD。第二个峰值发生在7350万年前，为光荚含羞草和费氏金合欢所共享，代表了一个更古老的WGD事件。

光荚含羞草和尖尾叶染色体之间的共线性显得零散，揭示了2:2的共线深度比（补充图9）。具体来说，光荚含羞草的染色体1与尖尾叶的染色体7和10的部分区域对齐，而光荚含羞草的染色体2主要与尖尾叶的染色体2的部分区域对齐。总体而言，光荚含羞草的大多数染色体与尖尾叶的染色体部分呈现一对一的关系（图2d）。

两种单倍型之间的结构变异和等位基因表达分析

在单倍型A和B之间共检测到7,606个结构变异（SVs），包括3,202个重复（DUPs）、20个倒位（INVs）、2,698个倒置重复（INVDPs）、815个倒置易位（INVTRs）和841个易位（TRANSs）（图3a和补充图10）。在这些SVs中，选择并手动验证了INVs，使用PacBio长读序列和Illumina序列，确认了20个倒位（补充图11a）。这些SVs影响了光荚含羞草基因组中的714个基因。对这些SV相关基因（位于SV区域的基因）进行了GO分析，这些基因主要富集在蛋白水解调节、水解酶活性调节、草食反应、伤口反应和外部生物刺激反应过程中（补充图12a）。KEGG分析结果显示，SV相关基因主要富集在核质转运、脂肪酸生物合成和苯丙素类生物合成过程中（补充图12b）。在上述GO和KEGG分析之后，我们继续注释显示显著富集的SV相关基因，以增强我们对这些SVs生物学意义的理解。我们观察到，不同类型的SVs中有些基因参与了叶片运动和根瘤形成过程。例如，来自INV的MbHA12g09480属于多抗菌挤出（MATE）家族，可能参与植物-微生物相互作用，并在根瘤形成过程中发挥重要作用，而来自INVDP的MbHA04g15050是Tre-2/Bub2/Cdc16（TBC）结构域家族的成员，可能通过调节参与囊泡运输和膜动态的小GTP酶来调节叶片运动。

两种单倍型之间的结构变异和等位基因比较 a 在两个单倍型之间检测到的结构变异和等位基因的共线性分析。从外到内的轨道：（a）染色体核型，（b）重复（DUPs），（c）倒位（INVs），（d）倒置重复（INVDPs），（e）倒置易位（INVTRs），（f）易位（TRANSs），（g）等位基因之间的共线性分析。所有分布图的窗口大小为1 Mb。b 等位基因对之间Ka/Ks的分布。c 两个单倍型之间等位基因FPKM（每百万映射片段的每千碱基转录片段数）值的箱线图。浅蓝色和深蓝色分别代表单倍型A和单倍型B。d 在两个单倍型中表达模式不同的四个基因。红色表示高表达，蓝色表示低表达。使用三个生物重复的平均FPKM作为表达水平，并通过Log2（FPKM + 1）进行归一化。e 与SV相关的基因和等位特异性表达基因的Venn图。

检测到390万个SNP和222,701个小的插入缺失变异（< 50 bp），包括112,097个插入和110,604个删除（补充表22），这些变异广泛且随机分布在13条拟染色体上（补充图13）。还检测到总共8,677个大插入缺失变异（> 50 bp），并随机选择了染色体1上的100个大插入缺失变异进行验证，使用PacBio长读序列验证了其中95个，使用PCR验证了随机选择的20个大插入缺失变异中的16个（补充图11b-d，补充表23）。还检测了每对同源染色体之间的序列同一性，没有间隙比对，这两个单倍型之间的序列同一性为99.4%（表1）。

通过分析两个单倍型之间的共线性，分析了等位基因和单倍型特异性基因（图3a）。其中，识别出14,685对具有两个等位基因（至少有一个氨基酸差异）的基因对，9,929对具有两个相同基因型的拷贝。这些等位基因对的平均Ka/Ks比率为0.46（图3b），表明大多数等位基因经历了净化选择。单倍型特异性基因包括单倍型A中的1,653个基因和单倍型B中的2,061个基因（表1）。单倍型A特异性基因主要富集在响应其他生物、响应外部生物刺激和响应外部刺激过程中，而单倍型B特异性基因主要富集在核酸磷酸二酯键水解、RNA修饰和RNA磷酸二酯键水解过程中（补充图14a和b）。这些结果表明，单倍型A在响应外部环境刺激方面贡献更大。此外，可以推断出单倍型A在光荚含羞草的叶片运动和根瘤固氮两个过程中可能发挥更大作用。

基于相位基因组，我们检查了等位基因之间的表达差异，发现大多数等位基因在两个单倍型之间的表达水平没有显著差异（图3c）。有趣的是，在六种组织中，有930对基因对在两个等位基因之间显示出差异表达（补充表24）。除了染色体6、10和13上的差异表达基因（DEGs）较少外，高表达水平的等位基因在两个单倍型和染色体之间随机分布（补充表24）。对DEGs的GO分析揭示了在触摸关闭的叶片（TCL）和根（R）中对生物刺激和外部生物刺激的响应途径富集（补充图15）。随后，这些基因进一步功能注释，揭示了一些基因可能参与叶片运动或根瘤固氮过程（补充表25）。例如，四个基因，包括MbCML、MbTIP、MbNRT1和MbGST，在两个单倍型中具有不同的表达模式，预测它们可能参与不同的过程。MbCML和MbTIP可能参与叶片运动过程，而MbNRT1和MbGST可能与根瘤固氮过程相关，需要进一步研究（图3d）。通过对SV相关基因和显示等位基因差异表达的基因进行关联分析，研究了SVs对等位基因表达差异的影响。我们的发现显示，89个SV相关基因在两个单倍型之间显示出差异表达（图3e，补充图16），表明SVs通过影响等位基因的差异表达参与了一些重要的生物过程。

分析光荚含羞草中参与叶片运动的关键基因

光荚含羞草的叶片运动包括地震感应运动和昼夜运动。对光荚含羞草的转录组测序以识别参与叶片运动的候选基因。对三种状态的叶片进行了测序，包括触摸前的开放叶片（OL）、触摸关闭的叶片（TCL）和由生物钟调节的夜间关闭的叶片（NCL）。每个样本进行了三次生物重复测序（补充表26）。在TCL vs. OL和NCL vs. OL中分别鉴定了574个和3,336个差异表达基因（DEGs）（补充表27）。

对三种叶片状态的DEGs进行表达聚类分析，揭示了六种不同的表达谱（图4a和4b）。值得注意的是，簇4和簇6中的DEGs在TCL中表现出高表达水平，而簇1和簇5中的DEGs在NCL中表现出高表达水平（图4b）。有趣的是，簇2中的一组基因在OL中高度表达，而在TCL和NCL中低表达（图4b），对应于叶片的开放和关闭状态。同时，簇3中的一组基因在OL和TCL中高度表达，而在NCL中低表达（图4b），对应于白天和夜间的叶片。

不同基因调控图谱在叶片运动中的差异表达基因 a 三种不同叶片阶段（OL、NCL、TCL）差异表达基因的热图。红色表示高表达，蓝色表示低表达。使用三个生物重复的平均FPKM作为表达水平，并通过Log2（FPKM + 1）进行归一化。 b 通过TCseq分析和模糊C-均值聚类确定的差异表达基因的六种表达模式类型。对簇质心具有高亲和力的显示为紫色，低亲和力的显示为绿色。

对这些在不同表达模式中的差异表达基因（DEGs）进行了GO和KEGG富集分析。在GO富集分析中，簇4和簇6的DEGs主要参与响应刺激、响应压力和响应含氧化合物的生物过程（补充图17a）。簇1和簇5的DEGs主要参与RNA生物合成过程的调控、节律过程和昼夜节律的生物过程（补充图17b）。簇2的DEGs富集在离子运输、跨膜运输和无机离子跨膜运输的生物过程中，而簇3的DEGs富集在响应非生物刺激、响应光刺激和响应温度刺激的生物过程中（补充图17c和d）。在KEGG富集分析中，簇4和簇6的DEGs主要参与次生代谢物生物合成、MAPK信号通路-植物和植物激素信号转导的生物通路，而簇1和簇5的DEGs主要参与植物激素信号转导、淀粉和蔗糖代谢以及缬氨酸、亮氨酸和异亮氨酸降解的生物通路（补充图17b和18a）。簇2的DEGs富集在次生代谢物生物合成、植物激素信号转导和昼夜节律-植物的生物通路中，而簇3的DEGs富集在代谢途径、光合作用和光合作用-天线蛋白的生物通路中（补充图18c和d）。总体而言，这些结果揭示了叶片运动是一个涉及许多生物过程和通路的复杂过程。地震感应运动主要与应激反应的生物过程相关，而昼夜运动主要与节律过程的生物过程相关。此外，离子跨膜运输和植物激素信号转导过程可能在地震感应和昼夜运动中都起作用，而光和温度也可能是影响昼夜运动的潜在因素。

结合功能注释分析和基因表达水平的考虑，帮助我们确定了最终的候选基因。对于地震感应运动，识别出57个在触摸关闭的叶片中高表达或特异表达的候选基因，主要是应激反应基因，包括锌指蛋白、乙烯响应转录因子和Nudix水解酶（补充图19c和20a，补充表27）。对于昼夜运动，识别出107个在夜间关闭的叶片中高表达或特异表达的候选基因，包括伤口诱导蛋白、水通道蛋白MIP、通用应激蛋白和UDP-糖基转移酶（补充图19d和20b，补充表27）。

研究了NCL vs. OL的DEGs表达谱，揭示了11个AQP基因、12个Ca2+相关基因和20个离子通道基因的显著表达差异（补充图19a, b和21b）。此外，与OL相比，在NCL中高表达且上调的两个AQP基因为Mb12g014040和Mb01g016320。随后，识别了光荚含羞草基因组中的33个AQP基因，并构建了包含拟南芥和蒺藜苜蓿等已广泛研究物种的AQP基因的系统发育树。MbAQPs被分类为5个亚家族（补充图21a，补充表28）。所有33个MbAQPs的表达谱显示PIP和TIP亚家族的成员在几乎所有组织中都有表达（补充图21b）。有趣的是，Mb12g014040和Mb01g016320都属于PIP亚家族，而另一个高表达且上调的基因Mb11g013010属于TIP亚家族。

分析光荚含羞草叶片运动的基因共表达网络

为了更好地理解光荚含羞草叶片运动相关基因的共表达动态，进行了差异基因表达分析，并使用基于表达基因的加权基因共表达网络分析（WGCNA）进行了加权基因共表达网络构建。该分析得到了9个簇（补充图22）。模块-组织关联分析表明，红色和黄色模块分别与TCL和NCL相关（补充图22）。为了识别可能参与叶片运动的重要基因，计算了红色和黄色模块中每个基因的模块特征基因连接性（kME）。分别识别出红色模块中的23个和黄色模块中的341个基因作为核心基因，并用于构建共表达网络（图5a和c，补充表29和30）。在TCL和NCL的共表达网络中，分别识别了8个和22个转录因子（TF）基因，属于不同的TF家族。热图显示了这些TF的表达谱（图5b和d）。有趣的是，在TCL的共表达网络中，识别了两个基因，MbCML27（Mb13g000970）和MbCaCA（Mb06g009710），而在NCL的共表达网络中，发现了一个名为MbCBL（Mb04g013530）的基因。这些基因可能在叶片运动中发挥重要作用。

叶片运动关键基因的共表达网络和表达谱 a 与TCL相关的WGCNA“红色”模块核心基因的共表达网络。包含kME-p值大于0.9的23个基因在Cytoscape生成的图中。橙色三角形代表转录因子，绿色圆形代表某些通路基因。

b 23个基因的热图，包括8个转录因子和15个通路基因。红色表示高表达，蓝色表示低表达。使用三个生物重复的平均FPKM作为表达水平，并通过Log2（FPKM + 1）进行归一化。

c 与NCL相关的WGCNA“黄色”模块核心基因的共表达网络。包含kME-p值大于0.9的341个基因在Cytoscape生成的图中。橙色三角形代表转录因子，绿色菱形代表参与植物激素信号转导通路的基因，粉色六边形代表参与植物昼夜节律通路的基因，紫色矩形代表参与缬氨酸、亮氨酸和异亮氨酸降解通路的基因。

d NCL共表达网络中22个转录因子基因的热图。红色表示高表达，蓝色表示低表达。使用三个生物重复的平均FPKM作为表达水平，并通过Log2（FPKM + 1）进行归一化。

光荚含羞草中关键根瘤基因的分析

豆科植物中，大多数物种的根瘤固氮是一个重要特性。近年来发表的30个关键根瘤基因被作为查询对象，共识别出光荚含羞草中的189个直系同源基因（补充表31）。189个候选根瘤基因的热图显示大多数基因在根中高度表达（图6a）。在这些高度表达的基因中，选择了对于根瘤形成至关重要的转录因子NIN（Nodule Inception）进行进一步分析。我们在拟南芥和其他六个物种中识别了NIN和NIN样蛋白（NLP）基因，这六个物种包括四个根瘤物种（光叶百脉根、蒺藜苜蓿、大豆和光荚含羞草）和两个非根瘤物种（C. canadensis和N. schottii）并构建了系统发育树（图6b）。结果显示，在七个物种中共识别出46个NIN和NLP基因，所有这些基因都包含RWK-RK和PB1结构域（补充图23）。然而，在多序列比对中，识别出光荚含羞草NIN的保守区域中一个与NLP对应区域不同源的区域（补充图24）。此外，NIN在非根瘤物种中缺失，而NLP在根瘤和非根瘤物种中均有发现。

结瘤相关基因的表达文件和NIN及NLP的进化分析 a 光荚含羞草中189个同源基因的热图。红色表示高表达，蓝色表示低表达。使用三个生物重复的平均FPKM作为表达水平，并通过Log2 (FPKM + 1) 进行归一化。 b 七个物种中NIN和NLP的系统发育树分析。

讨论

高质量的植物基因组和泛基因组组装促进了植物结构变异研究【35, 36】。单倍型解析基因组使得研究单倍型之间的SVs和等位基因特异性表达（ASE）成为可能【37】。最近的一项研究表明，SVs可以影响基因表达和基因型-表型关联，从而改良作物【38】。在光荚含羞草中，识别了714个SV相关基因。GO富集分析揭示这些基因主要与响应草食动物、伤口、外部生物刺激和其他生物过程相关。这表明SVs可能在塑造光荚含羞草的叶片运动和根瘤固氮特性中起作用。

单倍型A和单倍型B特异性基因的GO分析揭示了在不同生物过程中的富集。单倍型A基因组可能对光荚含羞草基因组中的叶片运动和根瘤固氮这两个典型特性贡献更大。此外，对ASE基因进行了GO分析，发现其在响应其他生物、生物刺激和外部生物刺激、草食动物和单萜生物合成及代谢过程的生物过程中的富集。这些发现为理解光荚含羞草叶片运动的机制提供了新的见解。

基因组和转录组数据的整合为基因发现提供了宝贵的资源，并促进了植物生长、次生代谢和生理过程各个方面的探索【39】。在识别的DEGs中，高表达和特异表达的基因被选为参与叶片运动的候选基因。对于地震感应运动，候选DEGs的功能主要涉及转录因子和细胞表面蛋白。这些基因可能介导信号传导，并在细胞发育、活化、增长和运动中发挥调控作用。对于昼夜运动，除了转录因子和细胞表面蛋白，我们还发现了一些与水分运输相关的基因和植物脱水蛋白家族的基因。这些基因可能通过细胞内外的水流调节叶片运动。此外，我们预测了在TCL vs. OL DEGs中，通过细胞壁果胶的去甲酯化修饰细胞壁的基因。同样地，在NCL vs. OL DEGs中预测了与微管和细胞壁构建相关的基因。这些发现与之前的研究一致，表明酪氨酸残基上的肌动蛋白磷酸化可以影响肌动蛋白丝的重组，从而导致含羞草中的地震感应运动【40, 41】。此外，次生细胞壁生物合成的变化可以影响叶片关节的物理强度，导致水稻旗叶叶片角度的变化【42】，这表明某些调控细胞壁形成的基因在光荚含羞草的叶片运动中起关键作用。

叶片运动是一个复杂的过程，涉及各种基因、生物过程和通路，作为植物对外部刺激的响应【43, 44】。Ca2+已被公认为是一个关键信使，在通过钙依赖的信号通路发挥植物应激反应中起重要作用【11, 45】。在植物细胞中，三种主要类型的传感器接收Ca2+信号，包括钙调蛋白（CaM/CAM）和CaM样（CML）蛋白、钙依赖性蛋白激酶（CDPKs/CPKs）、钙调磷酸酶B样蛋白（CBL）和CBL-互作蛋白激酶（CIPKs）【46】。其中，CaM和CMLs作为信号分子在信号传导和响应环境刺激中发挥作用。值得注意的是，GsCML27，一个Ca2+结合EF手蛋白，已被证明在植物对碳酸氢盐、盐和渗透胁迫的反应中发挥作用【47】。

转运蛋白在Ca2+和其他阳离子的运输中起关键作用，对矿物营养、离子胁迫耐受性和信号传导至关重要【48】。先前的研究揭示了CBL蛋白在植物非生物胁迫耐受性中的重要新作用【49, 50】。这些发现强调了Ca2+在信号传导和植物抗逆性中的潜在作用。我们在TCL的共表达网络中识别了一个钙结合蛋白MbCML27（Mb13g000970）和一个属于Ca（2+）阳离子反向转运蛋白（CaCA）（TC 2.A.19）家族的转运蛋白MbCaCA（Mb06g009710）。这一发现表明Ca2+和阳离子转运蛋白对叶片运动的潜在贡献。我们还在NCL的共表达网络中识别了一个钙调磷酸酶B样蛋白MbCBL（Mb04g013530），一个Ca2+传感器。这一发现表明Ca2+传感器可能感知外部刺激，传递信号，并可能参与叶片运动。此外，在TCL和NCL组织的基因共表达网络中还识别了AQP（水通道蛋白）基因。总的来说，这些发现使我们推测AQP基因和钙信号可能作为叶片运动和根瘤固氮的重要调控因子。

豆科家族中大多数物种共享根瘤固氮这一重要特性，92.27%的豆科物种表现出这一特性【25】。这种高频率的根瘤固氮是促进物种多样性增加的原因之一，使豆科植物成为第三大被子植物家族【25】。基于先前发表的30个关键根瘤基因的知识，我们在光荚含羞草中识别了这些关键根瘤基因的189个同源基因。通过进化分析，我们研究了根瘤和非根瘤物种中NIN和NLP基因的关系，确认NIN在非根瘤物种中丧失【25, 32】。此外，多序列比对揭示了NIN和NLP的N端区域的差异，表明它们在功能上的变化。之前在光叶百脉根中的研究表明，NIN的N端区域不响应硝酸盐，这导致了NIN丧失硝酸盐响应可能是Fabaceae共生固氮出现的必要条件【51】。这一适应使得在氮缺乏条件下诱导根瘤形成成为可能，这是共生固氮的特征【51】。保守结构域分析揭示了所有NIN和NLP基因中存在RWK-RK和PB1结构域，突显了它们的密切关系，并支持NIN作为NLP家族创始成员的观点【52】。

叶片昼夜运动和根瘤固氮是大多数豆科物种的显著特征，代表了两个耗能过程，而叶片地震感应运动主要存在于一些含羞草属的谱系中【2】。这种保留突显了这些特性的意义，并表明它们在使豆类适应环境方面的必要性。叶片地震感应运动被假设在含羞草属的八个谱系中独立进化，而含羞草和光荚含羞草被发现位于不同的谱系中【2】，表明这两个物种在基因组基础上具有不同的特性。最近的一项研究表明，含羞草是四倍体植物，并且其基因组已被发布【24】。我们的研究表明，光荚含羞草是二倍体物种，并且其基因组已组装完成。结合已发布的含羞草基因组，新基因组为从基因组和比较基因组学的角度提高对含羞草属叶片运动的理解提供了新的见解。此外，不同染色体倍性引起的剂量效应可能是叶片运动对外部刺激敏感性差异的原因。最近的研究集中在蒺藜苜蓿中的时钟基因LHY，揭示了其在根瘤形成中的作用。发现MtLHY功能的丧失会减少根瘤形成并损害氮同化，进而

影响根瘤中的内源性昼夜节律，最终影响昼夜叶片运动和生物量减少【53】，建立了昼夜叶片运动和根瘤固氮之间的关联。

结论

我们报告了一个高质量的、单倍型解析的、染色体级别的光荚含羞草基因组，提供了一个宝贵的遗传资源，用于进一步的研究，并作为豆类比较基因组学研究的关键参考。通过基于单倍型解析组装的结构变异分析和等位基因特异性表达研究，我们揭示了结构变异在光荚含羞草叶片运动和固氮中的潜在作用。此外，详细的转录组分析识别了许多有前途的候选基因，使得在基因水平上揭示植物叶片运动和固氮的分子机制成为可能。

材料与方法

植物材料和生长条件

光荚含羞草的种子来自中国西南地区的野生种质库（中国西南野生生物种质资源库-首页）。为了促进发芽，种子进行了95°C热水处理一分钟，然后放置在带滤纸的培养皿中。随后，将发芽的种子小心种植在温室环境中，温度保持在25°C，光周期为16/8小时。

染色体数目和倍性评估的统计

光荚含羞草植物在温室中培养约一个月后，取根尖进行染色体数目分析。我们遵循Xin等人的方法【54】，并对植物样品进行了轻微调整。酶解过程包括用含有4%纤维素酶、2%果胶酶和1%果胶分解酶的0.01 M柠檬酸缓冲液（pH 4.5）在37°C下处理根尖，酶解时间从1小时减少到30分钟。选择分散良好的有丝分裂染色体制备物，并用1.5 μg/mL 4,6-二脒基-2-苯基吲哚（DAPI）染色，然后进行染色体数目观察。用Lecia显微镜拍摄染色体扩展的图像。

基于K-mer分析，使用Smudgeplot（v0.2.2）进行倍性估算【55】。

DNA提取、文库构建和基因组测序

新鲜健康的叶片被迅速采集并在液氮中快速冷冻。随后，这些冷冻样品用干冰在运输过程中送至生物技术公司BioMarker。使用CTAB法从光荚含羞草叶片中提取高质量基因组DNA【56】。用Qubit2.0荧光计评估分离DNA的质量和数量。为了便于全面的基因组分析，按照每种技术的说明构建了PacBio HiFi（CCS）、Illumina和Hi-C文库。对于PacBio HiFi文库，整个基因组在PacBio Sequel II系统上进行单分子实时（SMRT）测序。Illumina文库按照Illumina提供的标准协议在Illumina HiSeq X Ten平台上进行测序，插入片段大小为350 bp，双端模式为150（PE150）。Hi-C文库在与Illumina数据相同的平台上进行测序，使用PE150模式。

通过FCM和K-mer方法估算基因组大小

使用流式细胞术（FCM）估算光荚含羞草的基因组大小，使用已知基因组大小为2300 Mb的玉米作为内部参考。

使用GCE（v1.0.2）计算基因组的K-mer深度和频率，精确估算基因组大小和异质性【57】。随后，使用R包进行K-mer直方图分析和K-mer分布图绘制，增强对光荚含羞草基因组特征的全面理解。

基因组组装和质量评估

使用HiFiasm（v0.15.5）组装基因组，通过Khaper（https://github.com/lardo/khaper）去除初步基因组组装中的冗余。使用Illumina数据通过NextPolish（v1.3.1）对共识基因组进行抛光【59】。为了将contig锚定到染色体并精炼以实现染色体级基因组，使用ALLHi-C将contig锚定到染色体，使用JUICEBOX工具进行组装错误校正【60, 61】。使用Ragtag（v2.1.0）将单倍型解析基因组与单倍体基因组作为参考进行拼接【62】。所有上述工具均以默认参数运行。

使用BUSCO（v3.0.2）和OrthoDB（embryophyta_odb10）评估最终基因组的质量【63, 64】。为了确认单倍体基因组的质量，用BWA-MEM（0.7.13-r1126）将Illumina数据和RNA测序数据映射到最终基因组【65】。LTR组装指数（LAI）也用于进一步评估包括单倍体和单倍型解析基因组在内的基因组组装质量【66】。使用Merqury（v1.3）基于K-mer分析进一步评估基因组组装质量【67】。

RNA提取、文库构建、RNA测序和分析

从不同组织中提取总RNA，包括开放叶片、触摸关闭叶片、夜间关闭叶片、茎、根和茎顶分生组织（SAM），使用多糖多酚植物总RNA提取试剂盒（TIANGEN, DP441），按照制造商的说明进行操作。每个样本有三个生物重复。用NanoDrop和Qubit2.0荧光计对分离的RNA进行质量评估。随后，将高质量RNA送至贝瑞基因公司进行文库构建和测序。

使用HISAT2将测序读段比对到组装的基因组上【68】。Stringtie2【69】用于通过每百万映射片段的每千碱基转录片段数（FPKM）计算所有基因的表达水平。DeSeq2【70】用于识别显著的差异表达基因（DEGs），阈值设为假发现率（FDR）≤0.05和绝对log2（Fold Change）≥1。

基因预测、结构和功能注释

GETA是一种自动化的全基因组注释工具，集成了基于同源性的方法、de novo方法和基于转录组的方法，用于预测光荚含羞草基因组的蛋白质编码基因（GitHub - chenlianfu/geta）。使用BUSCO（v3.0.2）评估注释的完整性【63, 64】。

为了识别非编码RNA元素，使用INFERNAL（v1.1.4）基于Rfam数据库预测rRNA、miRNA和snRNA，使用默认参数【71】。使用tRNAscan-SE（v2.0.9）识别tRNA【72】。

使用RepeatMasker（v4.0.9）和RepeatModeler（v1.0.8）对比RepBase数据库注释光荚含羞草基因组的重复序列【73, 74, 75】。使用TEclass对未知的重复序列进行分类【76】。使用LTR_FINDER识别完整的LTR逆转座子，使用Tandem Repeats Finder（TRF）检测串联重复序列【77, 78】。

使用EggNOGmapper在线软件对蛋白质编码基因进行功能注释，使用默认参数【79, 80】。此外，进行BLASTP（E-value=1e-5）搜索比对UniProt/SwissProt【81】、非冗余蛋白质序列数据库（NR）、真核生物完整基因组直系同源基因簇（KOG）【82】和京都基因与基因组百科全书（KEGG）数据库【83】。

结构变异检测和等位基因表达分析

使用Nucmer比对工具对两个单倍型的染色体对进行全基因组比对，参数为-c 500 -b 500 -l 50。使用子程序delta-filter过滤比对，参数为-i 90和-l 100，并使用show-coords将比对转换为制表符分隔文件，使用默认参数【84】。随后，使用SYRI（v1.6）检测结构变异，并使用circos（v0.69-8）【86】和plotsr（v0.5.4）【87】可视化结构变异的分布。使用子程序show-snps识别SNPs，使用Assemblytics（v1.2）基于Nucmer的比对识别INDELs。使用RIdeogram包描绘染色体上的SNP和INDELs分布【88】。

为了生成高质量的SVs，参考在狨猴二倍体基因组中的SV检测方法手动检查倒位，并通过以下步骤进行改进： (1) 修剪两个单倍型之间每个断点的上游/下游2000 bp序列； (2) 使用minimap2（v2.26-r1175）【90】和BWA-MEM（v0.7.13-r1126）【65】将局部PacBio HiFi读段和Illumina读段比对到带有侧翼序列的断点上； (3) 使用整合基因组浏览器（IGV）【91】可

视化PacBio读段和Illumina读段的比对。此外，使用GitHub上发布的代码（https://github.com/comery/marmoset）随机选择染色体1上的100个大INDELs，并使用PacBio长读段检查其准确性。按照茶基因组中的单倍型变异分析方法，比较同源染色体对之间的序列同一性，无比对间隙【37】。

为了识别等位基因，参考茶基因组中的单倍型比较和多样性分析方法，并使用GitHub上发布的代码（https://github.com/sc-zhang/AlleleFinder/）。

基因共表达网络分析

为了探索叶片运动相关基因之间的关系，使用R中的WGCNA（v1.71）【92】包进行加权基因共表达网络分析（WGCNA）。预过滤表达数据，标准是至少在一个组织中FPKM超过3的基因被认为是表达的。计算每个模块的特征基因值，并用来测试与每种组织类型的关联，以进行模块-组织关联分析。计算总连接度和模块内连接度，kME（基因表达水平与模块成员关系的Pearson相关系数）和kME-p值。对于核心基因选择，kME-p值超过0.9作为阈值，最终识别出组织特异性模块。使用Cytoscape（v3.9.1）可视化生成的网络【93】。

基因家族和系统发育分析

使用Orthofinder（v2.3.3）【94】以默认参数识别和聚类光荚含羞草及其他七个物种的基因家族，包括费氏金合欢、尖尾叶、降香黄檀、蒺藜苜蓿、拟南芥、番木瓜和水稻。基因组序列和蛋白质序列从公共数据库下载（补充表17）。

从Orthofinder（v2.3.3）结果中提取单拷贝直系同源基因，并用MAFFT（v7.307）【95】比对相应的蛋白质序列。随后，基于蛋白质比对使用RaxML（v8.2.12）构建最大似然系统发育树【96】。使用PAML包中的MCMCTree程序估计八个物种之间的物种分化时间，主要参数为：burn-in=2,000，sample number=20,000，sample frequency=2。此外，一个校准点是拟南芥和番木瓜之间的分化时间（68-72百万年前），另一个是单子叶植物和双子叶植物的分化时间（120-140百万年前），显示为拟南芥和水稻【98】。

基因家族扩展/收缩分析和光荚含羞草染色体上的物种特异性基因分布

使用CAFE（v4.2.1）【99】识别八个测序物种中发生扩展或收缩的基因家族，显著性阈值设定为P<0.05。从光荚含羞草中提取快速进化的基因家族，并使用OmicShare平台进行GO富集和KEGG通路分析（OmicShare 基迪奥生物信息云平台）。同时，提取光荚含羞草中的物种特异性基因，并在OmicShare平台上进行GO富集分析和KEGG通路分析。这些物种特异性基因在染色体上的分布由TBtools显示【100】。所有后续分析中的GO和KEGG富集均在OmicShare平台上进行。

WGD事件和共线性

光荚含羞草和尖尾叶同属于豆科苏木亚科，具有相同的基本染色体数（x=13）。为了探索光荚含羞草基因组的进化动态，使用WGD管道计算每同义置换位点的同义置换数（Ks）分布，旨在识别潜在的WGD事件【101】。基于值5.17×10^-3和公式T=K/2r（其中K代表基因组之间每个碱基的置换数），计算分化时间【102】。对于光荚含羞草和尖尾叶的比较基因组分析，使用MCScanX识别和可视化共线性块【103】。

AQP基因家族的研究

采用三种方法识别AQP基因家族，包括从光荚含羞草基因组注释中提取，基于AQPs的HMM轮廓（PF00230）和基于BLASTP。此外，从Phytozome数据库下载拟南芥和蒺藜苜蓿的AQP基因家族的蛋白质序列（Phytozome v13）。随后，使用ClustalW进行多序列比对，并在MEGA-X【104】中使用1000次bootstrap重复执行邻接聚类法。系统发育树使用在线工具Evolview进行可视化（https://evolgenius.info//evolview-v2/#login）。

光荚含羞草中关键根瘤基因的识别

使用30个公开的关键结瘤相关基因作为查询对象（https://github.com/Genomic-docker/Evolution-of-key-nodulaiton-genes）【25】。为了识别光荚含羞草基因组中这些关键结瘤相关基因的同源基因，进行了BLASTP搜索，阈值设为1e-5以下，并且至少有20%氨基酸序列识别度。用R显示这些同源基因的热图。识别了根瘤和非根瘤物种中NIN（Nodule Inception）和NLP（NIN-like Protein）基因的蛋白质序列。使用MEGA-X构建系统发育树，使用邻接聚类法和1000次bootstrap重复【104】。

数据和材料的可用性

Pacbio HiFi读段、HiC读段和RNA-seq读段的原始测序数据，以及基因组组装和注释已存储在国家基因组数据中心（Home - National Genomics Data Center），项目编号为PRJCA018116。原始DNA和RNA序列数据的登录号分别为CRA011718和CRA011721，公开访问地址为Genome Sequence Archive - CNCB-NGDC。全基因组序列数据的登录号为GWHDODZ00000000，公开访问地址为https://ngdc.cncb.ac.cn/gwh。