苔藓泛基因组--文献精读174

Bryophytes hold a larger gene family space than vascular plants

苔藓植物拥有比维管植物更大的基因家族空间

摘要

经过 5 亿年的进化，现存陆生植物构成两大姐妹类群：苔藓植物和维管植物。尽管苔藓植物体型微小、结构简单，却能在包括极端环境在内的多种生境中繁茂生长。然而，其生态适应性与长期存活的遗传基础尚未被充分阐明。

本研究整合 123 个新测序的苔藓植物基因组进行全面的超级泛基因组分析，结果显示：苔藓植物的基因家族多样性显著高于维管植物。这一特征体现为更多独特的、类群特异性的基因家族 —— 这些基因家族源于其漫长进化历程中广泛的新基因形成，以及持续的微生物基因水平转移。

苔藓植物丰富多样的遗传工具包（含独特免疫受体等生理新特征）的进化，可能为其扩散至不同生物群落提供了助力。这些新测序的苔藓植物基因组，为探究陆生成功的替代进化策略提供了宝贵资源。

正文

植物登陆是地球生命进化史上的关键节点 ¹²，最终孕育了如今复杂的陆生生态系统 ³。目前，全球约 40 万种陆生植物（有胚植物）分为两大主要类群：苔藓植物（即角苔、地钱和 moss）和维管植物（ tracheophytes）⁴。后者以具分枝的孢子体和木质化输导细胞为显著特征。

苔藓植物与维管植物之间、以及苔藓植物类群内部的系统发育关系曾长期存在争议。直到近年通过转录组和基因组数据解析，才明确苔藓植物是一个单系类群，且与所有现存维管植物互为姐妹群⁵⁻⁹。这些发现表明，现存陆生植物起源于约 5 亿年前分化的两条独立进化支系。

苔藓植物与维管植物的核心差异在于优势营养世代：苔藓植物以配子体为优势世代，维管植物则以孢子体为优势世代。因此，这两类植物的多样化可能主要（但非完全）由对各自优势世代相关性状的选择所驱动⁷。然而，支撑这些多样化过程的基因组进化模式是否存在差异，目前仍知之甚少 —— 这主要是由于不同类群苔藓植物的基因组数据极为有限。

全球苔藓植物约有 2.2 万种，分布广泛，常能在极端生境中繁茂生长甚至占据优势 ¹⁰¹¹。它们的结构相对简单：缺乏木质化输导细胞、无根 ¹²，孢子体无分枝且依赖母体配子体获取营养 ¹³。苔藓植物之所以能在生态上广泛成功，可能得益于其独特且分化的生理适应性（如脱水耐旱性）—— 这些适应性与维管植物存在显著差异 ¹⁴，且可能与它们出人意料的分化基因组相关 ¹⁵⁻¹⁹。但这种基因组多样性是否在苔藓植物中普遍存在，仍有待验证。

事实上，大多数苔藓植物目尚未完成基因组组装，而在约 70 个仍缺乏基因组信息的陆生植物目中，苔藓植物目占了相当比例 ²⁰。截至目前，55 个苔藓植物目中仅 8 个有测序基因组，这一知识缺口严重阻碍了对陆生植物基因组进化的深入理解 ²¹。本研究通过新测序 123 个苔藓植物基因组，填补了这一空白，实现了对苔藓植物的全面 “泛基因组” 分析。研究结果还揭示了塑造苔藓植物和维管植物基因组的、差异显著的基因家族空间及独特进化过程。

结果

苔藓植物超级泛基因组：核心基因、附属基因与特有基因

本研究新测序、组装并注释了 123 个高质量、高完整性的苔藓植物基因组，并整合至集中化平台（www.bryogenomes.org）。样本涵盖 37 种地钱、82 种 moss 和 4 种角苔，代表了 55 个已知苔藓植物目中的 47 个（图 1、扩展数据图 1 及补充数据 1）。

与其他陆生植物相比，苔藓植物基因组相对较小，平均基因数量也更少（总计约 27,959 个，其中角苔 19,267 个、地钱 27,768 个、moss 28,701 个；补充数据 2、3 及补充说明 1），而维管植物平均含 34,794 个基因。

系统发育分析显示，苔藓植物与现存维管植物互为姐妹群，且苔藓植物类群内部的亲缘关系与近年系统基因组学研究结果基本一致⁴⁸²²⁻²⁴（图 1、扩展数据图 2 及补充数据 1、4），为追溯苔藓植物基因组的进化提供了可靠框架（补充说明 1、2）。

该系统发育树基于 190 个类群的绿色植物系统发育分析（详见补充数据 1-3 及扩展数据图 2）构建，分歧时间源于苔藓植物的分子定年分析（补充图 19）。分支及物种名称按苔藓植物三大类群进行颜色标注。

从内到外，图中五条通道分别展示：基因组大小、基因数量、重复序列比例（基因组中重复序列的占比，蓝色、紫色和绿色柱形分别代表长末端重复序列（LTR）、长散在核元件（LINE）和未知重复序列）、基因长度及外显子长度。

左侧至右侧展示的苔藓植物图像（比例不一致）包括：角苔类的平滑褐角苔 ¹；地钱类的小型无萼苔 ²、合子地钱 ²、拉氏瓣叶苔 ²、紫叶羽苔 ² 和博兰德裂叶苔 ²；以及藓类的沼生泥炭藓 ²、金发藓属某种 ¹、无叶牛毛藓 ³、尖叶双扇藓 ²、缘毛大帽藓 ²、叶衣藓属某种⁴、白八齿藓 ²、紫红叶苔 ²、小红叶藓 ²、紫萼藓 ¹、南非卷叶藓 ²、黄边垂枝藓 ²、加州反叶藓 ² 和挺叶苔 ²。

插图基于以下作者提供的照片绘制：董世伟（平滑褐角苔、金发藓属某种、紫萼藓）、丹・卡拉汉（小型无萼苔、合子地钱、拉氏瓣叶苔、紫叶羽苔、博兰德裂叶苔、沼生泥炭藓、尖叶双扇藓、缘毛大帽藓、白八齿藓、紫红叶苔、小红叶藓、南非卷叶藓、黄边垂枝藓、加州反叶藓、挺叶苔）、刘阳（无叶牛毛藓）和张磊（叶衣藓属某种）。

注：LINE 为长散在核元件（long interspersed nuclear element）的缩写。

本研究将 343 个原始色素体生物物种（138 种苔藓植物、146 种维管植物和 59 种藻类；补充数据 5）的蛋白质组划分为 1,113,359 个直系同源群（基因家族）。

尽管维管植物的样本类群数量略多，但苔藓植物的非冗余基因家族累积数量显著高于维管植物（637,597 个 vs 373,581 个；图 2a）。苔藓植物中仅存在于单个类群的特有基因家族平均数量也更高（3,862 个 vs 2,223 个；图 2a 及补充数据 6），附属基因家族（存在于某一支系≥2 个样本但 < 80% 样本中的基因家族）数量同样更多（4,021 个 vs 1,583 个），而核心基因家族（存在于某一支系≥80% 样本中的基因家族）数量则较少（6,233 个 vs 6,647 个）。

此外，苔藓植物中特有或 “孤儿” 基因家族的占比虽较低（84% vs 87%；图 2a），但其绝对数量远超维管植物（532,840 个 vs 324,552 个）。总体而言，单个苔藓植物基因组中特有基因家族与附属基因家族的总数平均为 7,883 个（占基因家族总数的 56%），而维管植物中仅为 3,806 个（占 36%）。

a，138 种苔藓植物（左）和 146 种维管植物（右）单个基因组中各类基因家族（核心基因家族、附属基因家族、特有基因家族）的数量。x 轴上的类群按系统发育关系排序（扩展数据图 3）。虚线标注三类基因家族的累计平均值。饼图展示所有取样苔藓植物和维管植物中三类基因家族的总数。b，苔藓植物（蓝点）和维管植物（紫点）核心基因组大小减少（上图）及泛基因组大小增加（下图）的模拟结果。每个物种数量均采用 100 次随机组合分析。c，散点图显示，基于扩展数据图 3 的定年结果（x 轴），苔藓植物（上图）和维管植物（下图）系统发育树各祖先节点重建的基因家族总数（y 轴）。MRCA：最近共同祖先。

苔藓植物基因组具有极高的基因家族多样性 —— 随着基因组取样数量增加，直系同源群数量急剧上升（图 2b，下图），而核心基因家族在取样 30 个基因组后趋于稳定。维管植物基因组呈现类似模式，但基因家族累积曲线斜率更平缓（图 2b），这可能反映了苔藓植物更高的基因创新速率，或对苔藓植物的系统发育取样更全面（扩展数据图 3）。这些发现表明，尽管取样广泛，陆生植物中仍有大量基因和基因家族尚未被发现。

陆生植物基因家族获得与丢失的估计数量（扩展数据图 3 及补充数据 7）远高于此前报道⁷²⁵，这说明取样规模对重建基因家族进化具有深远影响（补充说明 3）。苔藓植物有着漫长的基因家族创新历史，尤其在白垩纪早期（约 1 亿年前）的藓类（如真藓亚纲；图 2c 及扩展数据图 3）中更为显著 —— 这一现象可能与藓类中连续全基因组复制（WGD）的累积效应相关（补充图 1 及补充说明 5）。与之形成对比的是，维管植物在过去 6500 万年分化出的支系中，祖先基因家族多样性较低，且总基因家族数量保持稳定且较少（图 2c 及扩展数据图 3），这可能与白垩纪 - 古近纪界线的大规模灭绝事件有关 ²⁶。长期来看，苔藓植物祖先节点的基因家族多样性估计高于维管植物（图 2c），凸显了苔藓植物中根深蒂固的基因家族创新过程。

苔藓植物的附属基因家族和特有基因家族是其基因组的重要组成部分：它们的鸟嘌呤 - 胞嘧啶（GC）含量与核心基因相似（约 50%；扩展数据图 4d 及补充图 2），且有相当比例（50%–80%）会表达（扩展数据图 4b 及补充图 3）—— 包括小立碗藓（Physcomitrium patens）、地钱（Marchantia polymorpha）等模式物种（扩展数据图 4a）。此外，小立碗藓中这些基因在多种胁迫条件下存在差异表达（扩展数据图 4c）。然而，苔藓植物中大多数附属基因家族和特有基因家族的功能仍不明确：基于蛋白质结构域的功能注释比例仅为 27% 和 16%，而核心基因家族的注释比例高达 91%（扩展数据图 4e）。这凸显了在理解苔藓植物基因功能方面存在巨大的知识缺口。

苔藓植物中孤儿基因的从头起源

在所有取样苔藓植物中（补充图 2），包括小立碗藓、地钱等模式物种（扩展数据图 4f、g），附属直系同源群和特有直系同源群中的基因比核心直系同源群中的基因含更少内含子，且编码区更短。这种特征在水稻 ²⁷、拟南芥 ²⁸²⁹等被子植物的新进化或 “年轻” 基因中也有发现，表明苔藓植物的大多数附属基因家族和特有基因家族可能起源较晚。此外，小立碗藓和地钱中特有基因家族的基因表达水平通常低于附属基因家族和核心基因家族（扩展数据图 4h 及补充图 3），这与水稻中新基因的表达模式一致 ²⁷。苔藓植物中的大多数特有直系同源群可能通过多种分子机制以孤儿基因的形式起源，包括快速序列进化和从头基因起源 ²⁷³⁰。在小立碗藓和地钱中，仅有不到 15% 的特有基因与其他直系同源群的基因存在序列相似性（扩展数据图 4i），提示它们源于基因复制；而 85% 的特有基因与现有直系同源群无序列相似性，可能通过不同机制起源，其中主要是从非编码区从头产生 ³¹。

由于序列快速分化 ³² 以及近缘类群基因组资源匮乏，鉴定祖先非编码序列存在一定挑战。本研究选取粗裂地钱（Marchantia polymorpha ssp. ruderalis），与两个同物种亚种及其姐妹物种进行比较 —— 这些类群在过去 500 万–700 万年内分化 ³³³⁴。粗裂地钱孤儿基因家族中约 70%–80%（3120–3583 个；扩展数据图 4j）的基因可与近缘物种的非编码区比对上，表明它们可能起源于这些区域（扩展数据图 4j）。

从头基因起源此前仅在水稻 ²⁷、竹子 ³⁰等被子植物中有报道，其过程涉及突变的逐步累积。例如，在粗裂地钱的多个特有基因中，可重建从非编码区衍生的原基因 ³⁵到成熟基因的进化轨迹，体现了原基因中终止密码子的逐步消除或转换（补充图 4）。以粗裂地钱的孤儿基因 Mp7g02380 为例，它与地钱属近缘亚种非编码区的直系同源序列相似度逐渐提高，同时终止密码子数量相应减少；另一孤儿基因 Mp4g17190 可能起源于地钱属亚种间高度相似的非编码区，通过单碱基突变将终止密码子转换为有义密码子。在藓类（尤其是真藓亚纲）中，约 26% 的特有基因与至少两个其他物种的非编码区高度相似（补充图 5，含示例性比对结果）。这是目前为主要植物支系报道的最大规模从头基因候选集，共计 36481 个基因。这些发现不仅证实了被子植物之外也存在从头基因创新，还凸显了其在苔藓植物中的重要性。该机制 ³⁶在推动苔藓植物适应性进化和生态多样化方面发挥关键作用，揭示了其基因组的动态特征及塑造基因组的进化过程。

频繁的水平基因转移（HGT）塑造苔藓植物进化

49 种链形植物的基因组中，包含 1809 个起源于原核生物、病毒、真菌或动物的直系同源基因家族，其中 1306 个为新鉴定家族（补充数据 8–11、图 3 及补充图 6–13）。水平基因转移的最大爆发发生在链形植物祖先中，涉及 115 次转移事件，其次是有胚植物祖先中的 90 次转移事件（图 3a）。平均而言，单个苔藓植物物种获得的水平转移基因多于维管植物（229 个 vs 163 个；图 3b）。对苔藓植物基因组的全面系统发育分析显示，其物种特异性水平基因转移事件数量更高（图 3b）—— 这一特征在小立碗藓近缘物种中也有近期报道 ³⁷。这表明水平基因转移可能是苔藓植物整个进化历程中的一个持续特征，而维管植物的多样化过程中未呈现该模式（图 3a）。

绿色植物中的水平转移基因通常比核心基因含更少内含子、长度更短，这与它们主要来源于细菌或真菌供体生物的基因特征一致 ³⁶（补充图 14）。在苔藓植物中，水平转移获得的基因常位于长末端重复序列（LTR）密度高于核心基因的区域，这暗示 LTR 元件可能促进了水平转移基因在其基因组中的整合或增殖（补充图 14）。

a，链形植物时间树（苔藓植物为蓝色，维管植物为紫色）上水平基因转移（HGT）事件数量汇总（补充数据 10）。灰色圆圈中的数字代表在对应节点重建的总 HGT 事件数。饼图展示选定祖先节点及特定支系（含苔藓植物或维管植物所有内部节点）中，HGT 推定供体的大类分布。末端分支旁的数字代表本研究数据集中该末端类群特有的 HGT 事件数。b，各取样物种中水平转移直系同源群（OGs）的累计数量柱状图，每个柱子拆分显示源于祖先 HGT 事件的基因家族占比。维管植物节点总和 = 所有维管植物内部节点的汇总数据；苔藓植物节点总和 = 所有苔藓植物内部节点的汇总数据。

苔藓植物通过水平基因转移获得的基因中，很大一部分功能未知，其余则参与复制、碳水化合物代谢、氨基酸代谢及次生代谢等过程（补充数据 12）。在小立碗藓中，大量水平转移基因在高温、干旱和紫外线胁迫下存在差异表达：分别有 151 个（60%）、148 个（59%）和 190 个（76%）水平转移直系同源群表现出表达变化（补充数据 13–22）。总计 233 个水平转移直系同源群的表达与非生物胁迫响应相关，占小立碗藓所有水平转移直系同源群的 93%，这暗示它们可能增强了苔藓植物在多样环境中的生态适应性。

此前已有研究报道其他水平转移基因具有胁迫响应功能。例如，藓类获得的细菌来源 HPA3 基因编码 N - 乙酰转移酶，可通过 N - 乙酰化作用解毒有毒 D - 氨基酸，这是苔藓植物在营养匮乏或胁迫生态系统中维持细胞稳态和存活的关键适应性特征 ³⁸³⁹。此外，地钱获得的真菌来源 RDS1 基因可能与胁迫响应相关⁴⁰，凸显了水平转移基因对植物环境韧性的贡献。综上，这些基因表明水平基因转移在塑造苔藓植物适应能力方面发挥关键作用，使其能在恶劣生态条件下繁茂生长。

本研究新鉴定的水平转移基因可能还在与其他生物的相互作用中具有多种适应性功能。例如，藓类中细菌来源的 LEC 基因编码甘露糖特异性凝集素，对链格孢菌（Alternaria alternata）和炭疽菌属（Colletotrichum）物种具有抗真菌活性⁴¹，凸显其在病原菌防御中的潜在作用⁴²。类似地，地钱中真菌来源的 YANB 基因编码参与 yanuthone D 生物合成的脱羧酶，对真菌和细菌具有抗菌活性，可能保护苔藓植物免受微生物竞争者的侵害⁴³。地钱中另一个真菌来源基因 EUPF 可生物合成 eupenifeldin，这种次生代谢产物具有广谱抗真菌、抗疟和驱虫特性⁴⁴，可能增强苔藓植物在多样生态压力下的存活能力。

水平转移获得的基因可能在苔藓植物防御植食性动物方面发挥重要作用。例如，在地钱中，水平转移获得的萜烯合酶（TPSs）被认为主要作为抗植食性动物物质⁴⁵⁴⁶。此外，许多地钱和藓类中存在从真菌获得的 FBT 基因，该基因编码含 FB 凝集素结构域（PF07367）的蛋白质，具有杀虫效果（补充图 15）。地钱合成的该蛋白质对注射 5.41 μg FB 凝集素的小菜蛾（Plutella xylostella）个体造成 80% 的死亡率⁴⁷。本研究进一步证实了该基因的抗植食性功能：即使使用较低剂量（0.45 μg・mg⁻¹，相当于每个个体 2.86 μg）的砂藓（Acrolejeunea sandvicensis）来源 FB 凝集素，也能导致棉铃虫（Helicoverpa armigera）幼虫 97.62% 的死亡率（扩展数据图 5b 及补充数据 23），其效力接近苏云金芽孢杆菌（Bacillus thuringiensis）分泌的 Bt 毒素⁴⁸。含 FB 凝集素的饲料还显著降低了草地贪夜蛾（Spodoptera frugiperda）的生长速率（扩展数据图 5c 及补充数据 24）。

因此，苔藓植物拥有大量源于微生物水平基因转移的多样化基因家族，其中许多为特有基因家族（图 3），这些基因构成了筛选潜在抗细菌或真菌病害基因（如 KP4⁴⁹）、探索其在植物发育⁵⁰及植物 - 真菌相互作用⁵¹ 中功能的巨大基因库。这凸显了水平基因转移在苔藓植物进化和生态多样化中的重要性。

讨论

与维管植物相比，苔藓植物具有显著更多的特有基因家族和附属基因家族（图 2），这些家族在其整个进化历程中不断出现（图 2）。尽管苔藓植物基因组的基因家族多样性更高（图 2），但其平均基因数量少于维管植物（27,937 个 vs 34,813 个）—— 维管植物的基因家族主要通过基因复制扩张（补充图 16）⁵²。苔藓植物的特有基因家族可能源于从头基因⁵³，这一基因创新过程此前仅在水稻、竹子等被子植物中有报道 ²⁷³⁰，本研究证实苔藓植物中也存在该过程。从头起源似乎是苔藓植物基因组多样化的关键机制。

本研究对植物界的全面基因组取样，大幅拓展了从头基因的系统发育分布范围，凸显了从头基因起源在植物基因进化和功能创新中的重要作用。新基因多样性的模式表明，在苔藓植物进化过程中，新基因持续出现并固定。苔藓植物中 “新” 基因的固定效率可能高于维管植物⁵⁴，因为这些基因在苔藓植物单倍体营养体中会立即受到选择压力。此外，苔藓植物配子体组织的全能性可能进一步促进遗传创新的传播，从而推动新基因的固定。

植物的许多新功能通过水平基因转移获得，本研究的密集取样表明，水平基因转移并非植物中零星发生的事件，而是持续的新基因来源，在苔藓植物中尤为显著。值得注意的是，苔藓植物中近 70% 的水平转移基因家族属于附属基因家族或特有基因家族（补充数据 25），且大多数取样类群均携带特有水平转移基因（图 3）。与种子植物相比，苔藓植物的水平转移基因积累存在差异，这可能与生殖细胞直接接触微生物的概率相关⁵⁵。

苔藓植物的结构本质上较为简单，但其遗传库的多样化（形成大量附属基因家族和特有基因家族，图 2a）可能是由对新生态生理适应性的选择及其优化所驱动 ¹⁴。例如，在小立碗藓中，63% 的附属基因家族和 27% 的特有基因家族对紫外线、高温、干旱等非生物胁迫有响应（扩展数据图 4c）。参与耐旱性的 PpARDT 基因⁵⁶最初被鉴定为小立碗藓的 “孤儿” 基因，而本研究发现它在更多藓类中普遍存在。其他生理创新可能包括对超冷环境的适应，使苔藓植物能在极地生态系统中繁茂生长⁵⁷。

苔藓植物还以其多样化的次生代谢产物为特征，尤其是萜类化合物⁵⁸，这些化合物由丰富的萜烯合酶（TPS）基因编码，特别是 TPS-c 型基因（补充图 17）。苔藓植物中超过 70% 的 TPS 基因属于附属基因家族和特有基因家族（补充数据 25），且在多个陆生植物支系中通过水平基因转移多次独立获得⁵⁹。在苔藓植物中，萜类化合物参与众多关键生理过程，包括化感作用⁶⁰、植食性动物 deterrence⁴⁶⁶¹、病原菌防御⁶² 及激素活性⁵⁸。此外，苔藓植物与被子植物在免疫受体方面的多样性和差异（扩展数据图 6、补充说明 4 及补充图 18），也表明苔藓植物具有独特的进化防御策略。

综上，本研究对陆生植物系统发育树中各类群基因组的全面分析表明，与维管植物不同，苔藓植物基因组的进化和多样化在很大程度上由持续的从头基因创新和微生物基因获取所驱动。尽管苔藓植物基因组的平均基因数量较少，但这些基因属于更多样化的功能基因家族。这种遗传多样性可能是这些结构简单、单倍体占优势的生物具有强韧性的基础，也可能解释了它们在过去 5 亿年中持续存在于陆地、广泛分布于各类生态系统，并成功占据全球极端环境的原因。

研究方法

类群取样、DNA 及 RNA 提取

苔藓植物样本采集自南极、比利时、加拿大、智利、中国、德国、新西兰、俄罗斯、南非及美国的野外生境，后续置于人工气候箱中培养。凭证标本存放于中国深圳仙湖植物园标本馆（SZG）和康涅狄格大学乔治・萨福德・托里标本馆（CONN）。

除 3 种角苔采用无菌培养材料（DNA 编号 639、899、902；补充数据 1）外，其余所有苔藓植物均以野外采集材料用于 DNA 提取和基因组测序。苔藓植物样本经蒸馏水至少清洗 3 次，在解剖显微镜下选取样本以避免潜在污染。

基因组 DNA 和 RNA 分别采用 FastPure 植物 DNA 提取迷你试剂盒（Vazyme）和 RNA-easy 提取试剂（Vazyme）提取。通过 1% 琼脂糖凝胶电泳、Qubit 荧光计（赛默飞世尔科技）及 NanoDrop 2000 分光光度计（赛默飞世尔科技）对 DNA 和 RNA 进行定量与质量检测。

本研究取样覆盖苔藓植物系统发育树的多个分支，涵盖了目前已确认的 55 个现存苔藓植物目中的 47 个，仅 8 个目未包含（即 Catoscopiales、Dendrocerotales、Disceliales、Leiosporocerotales、Neohodgsoniales、Oedipodiales、Phymatocerotales 和 Scouleriales）²⁴⁶³。

基因组测序与 k-mer 分析

短读长全基因组测序文库采用华大智造 MGIEasy FS DNA 文库制备试剂盒（货号 1000006988）或 Illumina TruSeq DNA 无 PCR 文库制备试剂盒构建，片段大小为 300-500 bp，在 MGI-SEQ 或 Illumina NovaSeq 6000 平台进行测序，生成 150 bp 双末端读长数据。

转录组文库经 polyA 筛选后，采用 Illumina TruSeq RNA 文库制备试剂盒 v2 构建，插入片段大小为 200-400 bp，在 MGI-SEQ 或 Illumina NovaSeq 6000 平台进行 150 bp 双末端读长测序。

长读长测序方面，选取长度大于 30 kb 的 DNA 片段构建文库，先用 NEBNext FFPE 修复混合液（新英格兰生物实验室）和 NEBNext Ultra II 末端修复 /dA 加尾模块（NEB）进行末端修复，再使用 SQK_LSK109 连接测序试剂盒制备测序文库。采用纳米孔 PromethION 测序仪为 38 个物种生成长读长数据；为 85 个物种生成单管长片段读长（stLFR）和 10x Genomics 读长数据。

stLFR 文库采用华大智造 MGIEasy stLFR 文库制备试剂盒（货号 1000005622）构建，在 MGI-SEQ 平台测序；10x Genomics Chromium 基因组文库按照制造商 protocol，采用 Chromium 基因组试剂试剂盒（v2 化学试剂）构建，在 MGI-SEQ 平台测序。

Hi-C 文库采用 DpnII 限制性内切酶，依据原位连接 protocol 构建，在 Illumina NovaSeq 6000 平台测序，生成 150 bp 读长数据。参照文献 64 的 protocol 制备 MethylC-seq 文库，在 Illumina NovaSeq 6000 平台测序，生成 150 bp 读长数据。

利用 TRIMMOMATIC v0.39⁶⁵对全基因组测序短双末端读长进行过滤，去除重复序列、低质量读长及接头序列。经质量控制后的读长进一步通过 GC 含量筛选（补充图 20）去除异常读长，过滤后的读长用于 k-mer 分析以估算基因组大小。使用 JELLYFISH v2.3.0⁶⁶和 KMERFREQ⁶⁷，设置参数 “-m -C” 进行 k-mer 频率计数；采用 GENOMESCOPE v2.0⁶⁸，设置参数 “-p 1”，其余参数默认，估算基因组大小。

基因组组装、Hi-C scaffolding 及去污染

采用 NEXTDENOVO v2.4.0（https://github.com/Nextomics/NextDenovo）组装纳米孔长读长数据，参数设置为 “seed_cutoff=45”。随后分别使用 RACON v1.4.7⁶⁹和 PILON v1.23⁷⁰，结合纳米孔长读长和短读长数据对组装得到的重叠群（contig）进行抛光优化。

对于 stLFR 数据，采用 STLFR2SUPERNOVA 流程，默认参数将 stLFR 读长组装为支架（scaffold）；使用 SUPERNOVA v2.1.1⁷¹，默认参数组装 10x Genomics 读长数据（补充图 21）。

原始 Hi-C 读长首先用 TRIMMOMATIC v0.39⁶⁵默认参数过滤（补充图 22），再通过 JUICER v1.6⁷² 提取有效数据。利用 3D-DNA 流程 v180922⁷³，基于交联 Hi-C 数据对组装错误的 contig 进行校正、锚定、排序和定向（补充图 22）；使用 JUICEBOX 组装工具 v1.11.08⁷⁴手动校正连接错误、易位和倒位（结果见补充图 23）。

对基因组草图进行潜在污染序列过滤：首先采用 BLASTN⁷⁵搜索，设置期望值阈值为 1×10⁻⁵（补充图 22），参考数据库为美国国家生物技术信息中心（NCBI）核苷酸数据库和非冗余蛋白质序列数据库（nt/nr 数据库，2020 年 5 月发布版），并为每个 BLAST 比对结果分配分类学归属。

具体流程为：将基因组草图序列与 nt 数据库进行 BLAST 比对，去除非有胚植物比对位点占比超过 50% 的支架（计算方式：非有胚植物比对位点数量 /nt 数据库总比对位点数量）；随后将保留的支架通过 DIAMOND v0.9.25⁷⁶采用 BLASTX 方式与 nr 数据库比对，去除非有胚植物比对位点占比超过 75% 的支架（计算方式：非有胚植物比对位点数量 /nr 数据库总比对位点数量）。

通过与 3 种模式苔藓植物（地钱M. polymorpha、小立碗藓P. patens、狭叶角苔Anthoceros angustus）的细胞器序列进行 BLASTN 比对，去除线粒体和质体序列。为评估去污染流程的可靠性，采用 10 kb 滑动窗口比较基因组草图和去污染后洁净基因组的 GC 含量 - 测序深度分布，结果表明该去污染流程有效（补充图 20）。

基因组注释

重复序列注释

结合从头预测和同源比对方法构建定制化重复序列文库：使用 PILER v0.4.1⁷⁷、LTR_FINDER v1.0.5⁷⁸、REPEATSCOUT v1.0.5⁷⁹和 REPEATMODELER v2.0⁸⁰构建定制重复序列文库，用于从头重复序列鉴定；采用 REPEATMASKER v4.1.2，基于 Repbase 数据库⁸¹ 进行同源比对重复序列鉴定。

长末端重复序列（LTR）插入时间计算流程：使用 LTR_STRUC⁸² 提取完整 LTR 元件的 5' 端和 3' 端序列；通过 EMBOSS v6.5.7.0⁸³ 的 DISTMAT 工具计算逆转录转座子 5'-LTR 和 3'-LTR 序列的 K 值；最终采用公式 T=K/2r 计算 LTR 插入时间（T），其中 r 为同义位点年均替换率，取值为 9.4×10⁻⁹⁸⁴。

基因结构注释

重复序列屏蔽后，采用 BRAKER v2.1.5⁸⁵流程，整合转录组证据、同源比对证据和从头预测证据进行基因结构注释。同源比对证据采用 6 种绿色植物的蛋白质序列，均来源于 Phytozome v13 数据库（https://phytozome-next.jgi.doe.gov/），包括拟南芥（Arabidopsis thaliana）、满江红（Azolla filiculoides）、地钱（M. polymorpha）、小立碗藓（P. patens）、槐叶苹（Salvinia cucullata）和卷柏（Selaginella moellendorffii）。采用 BUSCO v5.3.2⁸⁶，基于 Viridiplantae_odb10 数据库评估基因组组装完整性。

基因功能注释

通过 NCBI BLASTP⁷⁵，设置期望值为 1×10⁻⁵，将注释得到的蛋白质序列与 KEGG（https://www.genome.jp/kegg/）和 SWISSPROT（https://www.uniprot.org/）数据库进行比对；使用 INTERPROSCAN v5.51-85.0⁸⁷预测蛋白质结构域；利用 EGGNOG-MAPPER 在线服务器⁸⁸⁸⁹，基于直系同源基因簇功能分类对基因进行功能注释与分类（补充图 24）。

甲基化分析

为明确苔藓植物的甲基化模式，采用 BISMARK v0.23.1⁹⁰将全基因组亚硫酸氢盐测序（WGBS）生成的双末端读长比对至对应的代表性基因组。使用 bismark_methylation_extractor 工具，设置参数 “—CX_context –comprehensive”，提取 CG、CHG 和 CHH 类型的甲基化胞嘧啶（C），仅保留覆盖度≥5 条读长的甲基化胞嘧啶用于甲基化水平计算。

单个候选胞嘧啶的甲基化水平计算公式为：甲基化读长数 / 总比对读长数，即 mC/(mC+T)。基因甲基化分析中，将基因本体及基因侧翼区域（±2 kb）分别均等划分为 20 个区间（bin），每个区间的加权甲基化水平计算公式为：（mC1+mC2+…+mCn）/（mC1+mC2+…+mCn+T1+T2+…+Tn）。

全基因组复制（WGD）分析

使用 JCVI v1.1.8⁹¹ 生成散点图，设置共线性得分阈值（cscore）为 0.99，展示基因组内共线性区块的关系（结果见补充图 25）。采用 PAML 软件包⁹³ 的 yn00 程序，通过 Nei–Gojobori 法⁹² 计算共线性区块上旁系同源基因对的同义替换率（Ks）（结果见补充图 26）。

采用 KSRATES v1.1.1⁹⁴分析代表性苔藓植物物种旁系同源基因和直系同源基因的 Ks 年龄分布（补充图 27-31），并对同义替换率进行校正。藓类、地钱类和角苔类数据集的最大三元组数量分别设置为 14、8 和 7；对每个全旁系同源组的 Ks 年龄分布拟合指数 - 对数正态混合模型。

按照 KSRATES 的分析要求，使用 I-ADHORE v3.0.01⁹⁵推断锚定基因对用于后续 Ks 分析；通过锚定基因对 Ks 聚类，将其划分为不同组别，潜在对应不同的全基因组复制事件。选取 10 种代表性藓类物种，验证所有候选全基因组复制事件，避免重复取样（补充图 32-36；方法细节见补充说明 5）。

系统发育与分子定年分析

分别基于 54 个类群（补充数据 8）、190 个类群（补充数据 4 及扩展数据图 2）和 343 个类群（补充数据 5）的数据集进行系统发育分析。使用 ORTHOFINDER v2.4.0⁹⁶，默认参数对所有蛋白质序列进行基因家族聚类；采用 KINFIN v1.0.3⁹⁷，默认参数筛选核心单拷贝基因家族。

选取的基因家族蛋白质序列经 MAFFT v5.0⁹⁸比对后，使用 GBLOCKS v0.91b⁹⁹修剪模糊区域，参数设置为：-b3（最大连续非保守位点数量）=8，-b4（区块最小长度）=5，-b5（允许的缺口位置）= 半保留，-b6（使用相似性矩阵）= 是。分别采用 RAxML v8.2.12¹⁰⁰（串联法）和 ASTRAL v5.7.3¹⁰¹（合并法）进行最大似然树优化，RAxML 树推断采用 PROTOGAMMAAUTO 模型和 200 次快速自举检验（比对结果和系统发育树见 Figshare 数据库，https://doi.org/10.6084/m9.figshare.23528667）。使用 PHYPARTS v0.0.1¹⁰² 映射所有单基因树，观察每个节点的基因树不一致性。

分子定年分析中，通过 SORTADATE v1.0¹⁰³ 筛选进化历史一致的基因，将所得比对结果串联后，利用 BEAST v2.4.6¹⁰⁴的 BEAUti 工具生成 XML 文件。采用化石标定（补充数据 26），设置非相关对数正态松弛分子钟模型和 LG+G 模型（4 个位点速率类别），在 BEAST v2.4.6¹⁰⁴中进行定年分析。不同基因座单独分区，各分区模型参数独立。使用 TREEPL¹⁰⁵生成满足所有化石先验约束的带分支长度起始树。

进行两次独立的马尔可夫链蒙特卡洛（MCMC）模拟，每次运行 10,000 万代，燃烧期（burn-in）为 25%，每 5 万代取样一次，评估参数后验分布的可信度。使用 Tracer v1.5¹⁰⁶评估收敛性和混合度，确保所有参数的有效样本量（ESS）均大于 200。

基于基因的泛基因组分析

基于 343 个原始色素体生物基因组数据集（详见补充数据 5）进行泛基因组分析。使用 ORTHOFINDER v2.4.0⁹⁶，设置默认膨胀值（I=1.5），进行基因家族（直系同源群）聚类。

定义 “核心基因家族” 为某一支系中存在于 > 80% 样本的基因家族；“附属基因家族” 为存在于某一支系≥2 个但 < 80% 样本的基因家族；“特有基因家族” 为仅存在于单个样本的基因家族。采用 COUNT¹⁰⁷，基于多洛简约法（DOLLO PARSIMONY）模型推断系统发育树中各分支的基因家族获得与丢失事件。

利用维恩图工具（http://bioinformatics.psb.ugent.be/webtools/Venn/）分析 5 个原始色素体生物类群（藻类、维管植物、角苔类、地钱类、藓类）的特有和共有基因家族。采用 BLASTN，设置期望值阈值为 1×10⁻¹⁰，将粗裂地钱（M. polymorpha subsp. ruderalis）的特有基因与自身及其他 3 种地钱（M. polymorpha ssp. polymorpha BR5、M. polymorpha ssp. montivagans SA2、M. paleacea ssp. diptera MPA）的非编码基因组序列进行比对。

对于粗裂地钱的每个特有基因，提取其他 3 种地钱基因组中对应最佳比对结果的重叠群，使用 MAFFT v5.0 进行比对；在 Geneious v10.0.2（www.geneious.com）中分析比对结果，追溯粗裂地钱特有基因的进化历史。3 种地钱基因组（粗裂地钱、BR5、SA2）下载自https://marchantia.info/，MPA 基因组下载自 NCBI 数据库。

关键功能基因鉴定

利用 iTAK v1.6 软件（参考文献 108）鉴定新测序苔藓基因组的转录因子，并结合系统发育信息和结构域信息进行人工校正。通过 BLAST 比对（期望值 e≤1×10⁻⁵）拟南芥的植物激素相关基因，筛选目标基因组中的同源基因，再经系统发育分析和功能注释进行人工验证。从文献和公共数据库中挖掘与植物适应陆地环境相关的发育及生理生化网络关键基因：根 / 假根、根毛、侧根、维管组织和气孔的信号传导与发育相关基因参考参考文献 25；角质层生物合成与调控相关基因参考参考文献 109；三维生长相关基因参考参考文献 110；共生相关基因参考参考文献 111；叶片发育相关基因参考参考文献 112；耐旱相关基因参考参考文献 113；细胞壁生物合成与调控相关基因参考参考文献 114、115、116；孢粉素合成与调控相关基因取自参考文献 117；紫外防护相关基因取自 GO 数据库（GO:0009650；紫外防护）和参考文献 118；苯丙烷类、香豆素、类黄酮、木质素的生物合成相关基因，以及植物 - 病原菌互作、节律相关基因取自 KEGG 通路和 KO 数据库（https://www.genome.jp/kegg/）。以拟南芥相应蛋白质序列为查询序列，通过 NCBI BLASTP75 软件在 190 个物种的数据集中进行搜索，设定期望值阈值为 1×10⁻⁵以获取同源序列。进一步根据目标基因的特征 PFAM 结构域筛选结果，要求结构域比对期望值≤1×10⁻⁵且覆盖度超过查询序列的 60%。将各通路的基因家族导入 COUNT107 软件，推断这些基因在系统发育树中的获得与丢失事件。

采用参考文献 119 描述的方法鉴定 TPS 基因。利用 HMMER 3.0 软件（参考文献 120）的 hmmsearch 工具，以期望值 1×10⁻⁵为阈值，在本地 Pfam-A 数据库中搜索蛋白质序列。仅保留在四个 HMM 模型中获得最佳比对结果的序列作为候选 TPS 基因，包括萜烯合成酶 C 端结构域（PF03936）、TPS N 端结构域（PF01397）、TRI5 结构域（PF06330）和 SmMTPSLs 结构域（基于蛇足石杉中鉴定的 48 个微生物型 TPS 构建）。所有候选 TPS 蛋白质序列通过默认参数在 NCBI 非冗余数据库中进行 BLASTP 搜索。若前 10 个最佳比对结果均来自细菌和 / 或真菌，或与 SmMTPSLs 高度相似，则该 TPS 基因注释为 “微生物型 TPS 样”（MTPSL）蛋白。使用 MAFFT v5.0 软件（参考文献 98）对所有注释的 TPS 蛋白质序列进行比对，再通过 TRIMAL v1.4 软件（参考文献 121）修剪比对结果。将处理后的比对序列导入 FASTTREE122 软件，采用 JTT 模型构建最大似然树，并在系统发育树上进一步确认和鉴定植物型 TPS 和 MTPSL 蛋白。

NOD 样受体基因的鉴定详见补充说明 6。小立碗藓的胁迫实验及差异基因表达分析详见补充说明 7。

水平基因转移（HGT）分析

选取 54 个代表性绿色植物物种（补充数据 8）进行水平基因转移鉴定。尽管已通过严格流程对苔藓基因组组装结果进行去污染处理（“基因组组装、Hi-C scaffold 构建及去污染”）以排除潜在污染，但仍采用严格策略和严谨分析流程鉴定水平基因转移事件：剔除位于小 scaffold 上的候选水平转移基因，并确认其侧翼基因为真实植物基因（补充图 37）。

根据 NCBI 分类谱系，将 NCBI 非冗余数据库中的类群大致分为七类：细菌、真菌、古菌、病毒、后生动物、绿色植物及其他（除真菌、后生动物和绿色植物外的真核生物物种）。设定供体类群为细菌（蓝细菌除外）、真菌、古菌、病毒或后生动物的物种；受体类群为每个采样物种所代表的不同分类等级类群（如链形植物门、胚植物门、维管植物门、种子植物门等），具体列表见补充数据 9。

水平基因转移鉴定包括三个主要步骤：全基因组筛选候选获得性基因、排除污染、高可信度水平转移基因过滤，详细分析流程见补充图 37。针对每个物种，以其蛋白质序列为查询序列，通过 NCBI BLASTP75 软件搜索 2021 年 12 月下载的 NCBI 非冗余数据库，参数设置为 “-e-value 1e-5 –outfmt 6 –max_target_seqs 50,000”，并采用参考文献 47 描述的方法鉴定潜在水平转移基因。简要来说，若 BLASTP 最佳比对结果属于候选供体类群，且来自其他类群（受体类群除外）的比对结果数量超过 5 个，则该查询基因为候选水平转移基因。随后，通过检测邻近基因的垂直遗传模式、候选基因在大 scaffold 上的位置以及同源基因在近缘类群中的分布，降低候选获得性基因源于序列污染的可能性。最后，通过一系列参数筛选高可信度水平转移基因：供体类群与受体类群的序列一致性需 > 20%；供体类群中的物种数量需 > 50；除受体类群外，BLAST 结果前 50 个物种应主要属于候选供体类群（细菌、真菌、古菌、病毒或后生动物）；序列一致性超过 80% 的候选转移基因需仔细检查其近缘类群中的同源基因；源于后生动物的水平转移基因应仅分布于后生动物和受体类群；排除源于蓝细菌或在受体类群起源前获得的候选转移基因。将校正后的水平转移基因映射到基于 54 个物种系统发育树、通过 COUNT 软件推断的基因家族演化结果中。

对于每个潜在的水平转移基因，将其在定制数据库中通过 BLAST 搜索获得的同源序列，根据 NCBI 分类数据库分为七类：细菌、真菌、古菌、病毒、后生动物、绿色植物及其他真核生物物种。进一步将这七类划分为更低分类等级的亚类群（补充数据 9）。每个亚类群至少选择一个序列，以确保采样充分且均衡。对物种数量少于 3 个的亚类群进行人工检查，排除测序污染（如非冗余数据库中的序列标注错误）。使用 MAFFT v5.0 软件（参考文献 98）默认参数进行多序列比对，通过人工或 TRIMAL v1.4 软件（参考文献 121）（参数 “-automated1”）去除比对质量较差的区域和缺口。利用 IQTREE v2.2.0 软件（参考文献 123）构建最大似然树，自动选择最佳氨基酸替代模型，并通过 1000 次超快自举检验（参数 “-bb 1,000”）评估分支支持度。

FBT 基因功能的实验验证

为进一步确认 FBT 基因（含真菌凝集素结构域 PF07367 的基因）的功能，通过体外表达该基因获得 FBT 蛋白，将其添加到昆虫饲料中，统计昆虫的生长率和存活率。优化桑福德阿氏藓（A. sandvicensis）FBT 基因的全长编码序列密码子，以提高其在原核系统中的表达效率，合成优化后的序列并连接至 pET-30a 载体。将 pET-30a-FBT 质粒转化至大肠杆菌 BL21（DE3）pLysS 感受态细胞，利用 IPTG 诱导型 T7 启动子调控重组蛋白表达。将过夜培养的菌液按 20 倍比例稀释到 200 毫升新鲜 LB 培养基中，37℃、220 转 / 分钟培养至 OD600 值达到 0.5。加入 IPTG 至终浓度 1mM 诱导表达，28℃孵育 3-4 小时。收集细胞并悬浮于 20 毫升预冷的 50mM Tris-HCl 缓冲液（pH 8.0，含 200mM NaCl）中，4℃条件下超声破碎。细胞裂解液在 4℃、12,000×g 条件下离心 40 分钟，取上清液上样至经 Tris-NaCl 缓冲液预平衡的 Ni-NTA 琼脂糖柱（4℃）。用含 20mM 咪唑的 Tris-NaCl 缓冲液充分洗涤柱子，再用含 250mM 咪唑的 Tris-NaCl 缓冲液洗脱 6×His 标签蛋白。将含纯蛋白的洗脱产物用 Tris-NaCl 缓冲液洗涤三次，通过 Centricon 浓缩管（Millipore，PM10）浓缩。采用 HRP 标记的单克隆抗体（HRP-66005）通过蛋白质印迹法（Western blot）检测纯化的 His 标签蛋白，该蛋白分子量约为 20kDa。

实验所用的秋黏虫（S. frugiperda）和棉铃虫（H. armigera）均来自深圳农业基因组研究所，采用人工饲料饲养（参考文献 124）。为避免自相残杀，三龄后幼虫单独饲养在 5.0 厘米 ×1.5 厘米的玻璃管中。选取两天龄的三龄秋黏虫和棉铃虫幼虫进行实验，其平均体重分别为 5.15 毫克 / 头和 6.35 毫克 / 头。在幼虫饲料中添加毒素，设置毒素初始浓度梯度为 0、1、2、4、8 和 16 微克 / 毫升饲料，以添加 PBS 缓冲液的饲料作为对照。每个处理组设置 14 头幼虫，重复三次。记录幼虫取食含毒素饲料 48 小时后的死亡率。称量幼虫初始体重（M1）并转移至含毒素的饲料中，48 小时后再次称量体重（M2），相对生长率（%）计算公式为 [（M2-M1）×100/M1]。浓度为 4 微克 / 毫升饲料的毒素（对应每毫克幼虫体重 0.47 微克毒素）可导致 97.62% 的棉铃虫死亡；而该 FBT 毒素对秋黏虫无致死作用，8 微克 / 毫升饲料的毒素浓度（对应每毫克幼虫体重 1.5 微克毒素）可显著抑制秋黏虫幼虫的生长率。

统计学分析

采用双侧 Welch t 检验分析 PBS 对照组与 FBT 蛋白处理组之间，秋黏虫幼虫体重变化动态和棉铃虫幼虫死亡率的统计差异。该方法适用于处理组间方差不等的情况，可确保比较结果的可靠性。所有分析均使用 R 4.2.2 版本（R 统计计算基金会）完成。