恩施碎米荠基因组--Cell Discovery

The Cardamine enshiensis genome reveals whole genome duplication and insight into selenium hyperaccumulation and tolerance

恩施碎米荠基因组揭示全基因组复制事件及硒超富集与耐硒机制

摘要

恩施碎米荠（Cardamine enshiensis）是知名的硒超富集植物。硒是一种对健康有多方面益处的必需微量元素。尽管该物种研究价值突出，但相关基因组信息仍然匮乏。本研究报道了恩施碎米荠的染色体水平基因组组装结果，基因组大小为 443.4 Mb，共 16 条染色体，scaffold N50 达 24 Mb。为阐明恩施碎米荠耐硒与硒超富集的分子机制，本研究整合基因组、转录组与代谢组数据集开展分析。结果表明，类黄酮、谷胱甘肽及木质素生物合成途径在保护恩施碎米荠抵御硒胁迫损伤中发挥重要作用。染色质互作的 Hi-C 分析显示，恩施碎米荠染色质可划分为 A、B 两个区室；每条染色体两端端粒间存在强互作，该特征与组蛋白修饰、表观标记、DNA 甲基化及 RNA 丰度相关。外源施硒可在染色质区室层面改变恩施碎米荠的三维染色质结构。硒处理后发生区室转换的基因参与硒化合物代谢；拓扑关联结构域（TAD）边界区域内的基因参与细胞硒响应、硒结合以及类黄酮生物合成过程。该多组学研究为解析恩施碎米荠耐硒与硒超富集的分子机制提供了理论依据。

引言

硒（Se）是人体必需微量元素，硒掺入硒蛋白后具备抗氧化、抗炎以及调控甲状腺代谢的功能\(^{[1]}\)。人体硒摄入不足会提升死亡风险，造成免疫功能低下、认知衰退以及不可逆脑损伤\(^{[2]}\)。补硒能够上调谷胱甘肽过氧化物酶 4（GPX4）的活性与转录水平，有效抑制 GPX4 依赖的铁死亡，改善出血性与缺血性中风的预后\(^{[3]}\)。近期研究发现，补硒可显著提升硫氧还蛋白还原酶 1（TXNRD1）活性，大幅增强顺铂及 TXNRD1 抑制剂的抗肿瘤效果（Huang 等人，未发表数据）。

恩施碎米荠最早在中国恩施硒矿区被发现\(^{[4]}\)，现已作为新型食用植物开展规模化种植。该植物在 400 μM 硒处理条件下培养 3 个月，生长未出现明显受抑，干重中硒积累量可达 3.7% 以上\(^{[5]}\)。因此，恩施碎米荠在硒污染水土的植物修复等方面具备潜在应用价值。事实上，硒相关产业经济贡献占恩施市年度 GDP 近 50%。由此可见，解析恩施碎米荠耐硒与硒超富集的分子机制，兼具环境价值与经济意义。

结果与讨论

基因组测序、组装及注释

本研究对恩施碎米荠开展基因组测序，将重叠群（contig）锚定组装得到 16 条假染色体，contig 可覆盖基因组 86.6%（图 1）。最终组装基因组总长度 443.46 Mb（染色体倍性 2n=32），contig N50 为 1.23 Mb，scaffold N50 达 24.41 Mb（表 1，补充表 1）。该组装基因组大小与 Kmer 分析预估基因组大小 481.37 Mb 较为接近（补充图 2）。通用单拷贝直系同源基因评估（BUSCO）结果显示，完整基因检出率达 97.7%\(^{[6,7]}\)（补充表 1）。

本研究采用 Maker 注释流程\(^{[8]}\)，整合同源蛋白比对、从头基因预测、转录组证据完成基因组注释，最终获得 52 725 个基因模型；基因平均长度 2.1 kb，编码序列平均长度 1.1 kb，单个基因平均包含 5.1 个外显子（补充表 1）。绝大多数预测基因（96.0%）获得同源序列支持并完成功能注释，76.3% 的基因可比对到 InterPro 数据库\(^{[9]}\)（补充表 1）。借助 InterProScan，基于基因本体（GO）分类完成 27 391 个基因的功能注释；同时预测得到 3324 个非编码 RNA（ncRNA）（补充表 1）。重复序列占基因组的 61.4%，其中长末端重复（LTR）反转录转座子占比 48.5%（补充表 1）。

a 恩施碎米荠 16 条染色体特征，窗口大小 200 kb。图中依次展示基因密度、重复序列密度、GC 含量、RNA 丰度以及共线性区块间关联。 b 硬毛碎米荠旁系同源与直系同源基因共线性区块的同义替换率（Ks）分布。 c 硬毛碎米荠与恩施碎米荠基因组之间共线性关系的 Circos 圈图。 d 基于 9 个物种的直系同源基因构建的推断系统发育树，标注物种分化时间与全基因组复制事件（WGD）。所有分支的后验概率均大于 0.99。α 与 β：古老 α、β 全基因组复制事件，分别发生于约 4700 万年前、1.24 亿年前。

恩施碎米荠基因组演化与比较基因组分析

共线性分析检测到一次全基因组复制（WGD）事件，以及全基因组复制之后发生的片段重复事件（图 1；补充图 2、补充表 1）。KEGG 富集分析显示，铁死亡通路为富集程度最高的条目（补充表 1）。系统发育分析结果表明，碎米荠属与拟南芥的分化时间约为 1600 万年前；恩施碎米荠与硬毛碎米荠（Cardamine hirsuta，2n=16）的分化时间约为 650 万年前\(^{[10]}\)（图 1）。结合 Ks 值与物种分化时间，计算得到每个位点每年的同义替换速率为\(8×10^{-9}\)，据此推算该全基因组复制事件发生于约 500 万年前（图 1；补充图 2）。Circos 圈图显示恩施碎米荠与硬毛碎米荠之间存在清晰的 2:1 共线性关系，硬毛碎米荠各条染色体与恩施碎米荠对应染色体高度共线性（图 1），该结果提示恩施碎米荠为四倍体。

硒对恩施碎米荠三维基因组结构的影响

为探究硒对恩施碎米荠三维基因组结构的影响，本研究分别对硒处理前后的材料开展 HiC 实验。HiC 技术已广泛应用于细菌、酵母、拟南芥、棉花、水稻、玉米、小鼠及人类的基因组研究\(^{[1120]}\)。本研究利用 HiSeq 测序平台，采用 150 bp 双端测序（150PE）模式，获得共计 2300 万条有效双端读段，用于三维基因组比较分析。同时，利用相同叶片组织样本开展基因组重测序、DNA 甲基化测序、转录组测序，以及针对两种组蛋白修饰的染色质免疫共沉淀测序（ChIPseq）。使用 ICE（Iterative Correction and Eigenvector）软件完成数据归一化处理\(^{[21]}\)，分别构建 400 kb 分辨率的全基因组 HiC 互作图谱与 100 kb 分辨率的染色体水平 HiC 互作图谱（图 2）。

恩施碎米荠基因组共线性分析鉴定得到 8 对推测的部分同源染色体（图 1）。部分同源染色体之间的互作占全部染色体互作的 10.63%27.29%（补充表 1）。结果显示，恩施碎米荠染色体端粒主要发生各染色体端粒之间的相互作用，该特征与拟南芥类似\(^{[12]}\)。本研究同时检测到着丝粒周边区域较强的染色体内独立互作以及染色体间互作。基于三维基因组图谱共识别出 25000 对显著互作位点：对照组包含 1186 个顺式互作、7050 个反式互作；硒酸盐处理组包含 774 个顺式互作、7345 个反式互作（图 2）。

为更加直观展示两组样本染色质互作模式差异，参照 Crane 等人的方法\(^{[22]}\)，将各组互作矩阵转换为 Zscore 矩阵，两组 Zscore 矩阵相减得到差值互作矩阵（补充图 2）。互作衰减指数（IDE）用于描述染色质互作频率随距离变化的衰减趋势\(^{[2324]}\)。为明确硒酸盐是否改变恩施碎米荠染色体空间结构，整合两组互作衰减曲线，比较全基因组互作衰减的整体变化趋势\(^{[15]}\)，两组样本间存在明显差异（图 2；补充图 2、补充表 1）。对照组与硒处理组全基因组 IDE 值分别为0.908、0.8518，该数值与拟南芥、水稻以及后生动物的 IDE 特征保持一致\(^{[2426]}\)。

通过全基因组特征向量分析鉴定染色质 A、B 区室\(^{[14]}\)，该分区模式同样存在于拟南芥\(^{[24]}\)、玉米、番茄\(^{[13]}\)和水稻\(^{[26]}\)中。利用 Cscore 指标解析恩施碎米荠基因组区室特征，可清晰区分 A、B 两类区室（图 2）。A 区室基因密度更高，富集活化表观标记（H3K4me2），转录活性高；B 区室抑制型表观标记（H3K27me3）占比更低，胞嘧啶甲基化水平更高，转座子密度更大（补充图 2）。值得注意，硒酸盐处理可重塑 2、6、8、9、10 号染色体的 A/B 区室结构（图 2；补充图 2）。

以 100 kb 分辨率比较硒酸盐处理组与对照组的染色质区室，共鉴定 755 个区室发生转换。保守 A 区室基因密度高、富集活化表观标记；保守 B 区室与发生转换的差异区室转座子密度升高，抑制型表观标记水平下降（补充图 2）。对发生区室转换区域内的基因集开展 GO 与 KEGG 功能富集分析（补充图 2），显著富集的通路主要包括乙醛酸与二羧酸代谢、脂肪消化吸收、植物激素生物合成以及硒化合物代谢。

a 恩施碎米荠染色体（1–16 号）内部及染色体间 Hi-C 互作图谱（400 kb 窗口）。所有染色体的常染色质臂之间均存在染色体内互作。颜色深浅代表两个 400 kb 位点之间的互作频率。对角线上方方框：同一条染色体内的顺式互作（cis）。白色行列代表无有效互作数据的窗口。 b 400 kb 分辨率下两组样本全基因组距离 - 互作频率关系图；有无硒处理条件下各染色体的互作衰减指数（IDE）。颜色代表不同样本的互作衰减曲线；横坐标：染色体不同位点间的相对距离；纵坐标：互作频率。 c Circos 圈图展示全基因组显著顺式互作位点（上）与反式互作位点（下）。 d 6 号染色体 Hi-C 互作图谱，以及基因组组分与多种表观标记的共定位。上方热图为 100 kb 窗口两两之间的染色体内 Hi-C 互作频率；下方依次展示 A/B 区室 PCA 特征向量、基因组与表观组信号轨道（均为 100 kb 窗口），包含基因与转座子丰度、不同组蛋白修饰、mRNA 表达量（FPKM 标准化）以及 CG、CHG、CHH 三种序列环境下的 DNA 甲基化。 e、f 一段基因组区域的 Hi-C 互作矩阵，展示拓扑关联结构域（TAD）（e：chr6:116240000–124160000；f：chr6:17280000–24840000）。上方：Hi-C 互作矩阵；下方：TAD 边界（竖线）与绝缘分值。图中纵轴与蓝线代表绝缘分值，灰线标记 TAD 边界。

拓扑关联结构域（TAD）是基因组空间结构的基本单元。这类结构广泛存在于哺乳动物基因组\(^{[28]}\)，植物中也鉴定出大量类 TAD 结构。TAD 边界富集启动子相关转录因子、转录起始位点、管家基因、tRNA 基因与短散在核元件（SINE），对维持 TAD 结构与稳定性至关重要\(^{[28]}\)。哺乳动物基因组 TAD 边界通常富集 CTCF 结合位点\(^{[29]}\)；而植物基因组 TAD 边界富集多重组蛋白修饰信号，染色质更为开放\(^{[17]}\)，基因表达水平更高\(^{[20]}\)。本研究基于 40 kb 分辨率 Hi-C 互作图谱在恩施碎米荠中鉴定 TAD（图 2e、f）：对照组与 400 μM 硒酸钠处理组分别得到 537、543 个 TAD，以及 521、527 个 TAD 边界（补充表 1）。有意思的是，TAD 边界区域的基因密度与基因表达水平高于 TAD 内部区域（补充图 2）。

为在 TAD 层面考察基因组甲基化位点变化，本研究采用滑动窗口法分析两组样本间的绝缘性差异\(^{[22]}\)。相比于 TAD 内部，TAD 边界富集活化表观标记 H3K4me2，且 CHG 环境下胞嘧啶甲基化水平更高（补充图 2）。基因表达数据分析显示，边界区域基因富集于细胞硒响应、硒结合以及类黄酮生物合成过程。对发生变化的边界区域基因开展注释、GO 与 KEGG 富集分析：硒酸钠处理后新增边界内的基因显著富集在细胞硒响应与硒结合通路（补充图 S9）。同时采用滑动窗口（窗口大小 11，步长 1）筛选两组绝缘分值存在差异的区域（所有重叠窗口 Pearson 相关系数＜0.6）\(^{[22]}\)，该区域基因显著富集于类黄酮生物合成通路（补充图 2）。综上结果表明，硒耐受与硒代谢在染色质层面存在关联。此外，硒酸钠可重塑metE、GS、PAL基因所在区域的 TAD 结构（图 4、图 5；补充图 2）。

高频互作区域（FIREs）是局部染色质互作热点，与染色质区室、TAD 和染色质环结构不同\(^{[22]}\)。FIRE 代表局部高频染色质互作区域，已被证实是增强子富集区，且超增强子出现概率更高\(^{[30]}\)。本研究在 10 kb 分辨率下共鉴定 1184 个显著 FIRE；这类区域显著富集在 2、9、13 号染色体、A 染色质区室以及 TAD 边界区域（补充图 2c 及补充图 2）。对 FIRE 集中区域的基因密度与 GC 含量分析显示，其与其他区域无显著差异，说明该区域并非被浓缩异染色质隔开的小型基因岛（补充图 2）。

恩施碎米荠硒耐受与硒超富集机制

为解析恩施碎米荠硒耐受机制，本研究用 400 μM 硒酸钠处理幼苗 24 h，清水作为对照。共鉴定 29671 个差异表达基因，占已注释基因总数的 66.6%（补充表 1）。为进一步研究硒处理后代谢组变化，取两组叶片开展代谢物定量，采用基于广靶液相色谱 - 串联质谱（LCMS/MS）的代谢组分析方法。在恩施碎米荠叶片中共鉴定 558 种代谢物（补充图 2、补充表 1），其中 127 种为差异代谢物。高浓度硒会使植物体内活性氧（ROS）与活性氮（RNS）大量累积，诱发氧化胁迫\(^{[31]}\)。类黄酮可作为金属螯合剂并清除 ROS，提升植物对重金属的耐受性\(^{[32]}\)。本研究发现两组样本间共有 10 个黄酮类代谢物含量发生变化，说明黄酮在硒耐受中起到关键作用（补充表 1）。对差异代谢物开展 KEGG 富集，显著富集通路主要集中于次生代谢物生物合成、黄酮 / 黄酮醇类化合物合成通路（补充图 2）。三蓟素 O - 丙二酰己糖苷与穗花杉双黄酮降至检测限以下，而木犀草素 O - 己糖基 - O - 己糖基 - O - 己糖苷上调 24444 倍（图 3；补充表 1）。全部代谢物主成分分析（PCA）结果显示两组样本转录组与代谢组存在明显差异（补充图 2）。

为解析转录组与代谢组的关联模式，基于两组数据开展相关性分析，设置 Pearson 相关系数阈值 r＞0.8，筛选与代谢物显著关联的基因。最终得到 105268 组表达相关关系，涉及 183 种代谢物与 3202 个基因（图 3）。进一步整合多组学数据构建关联网络，用于挖掘代谢通路与候选基因。选取 2000 个转录本与 2 种类黄酮代谢物开展 Pearson 相关分析，结果显示 175 个转录本与三蓟素 O - 丙二酰己糖苷、穗花杉双黄酮高度相关（\(R^2>0.96\)）（补充表 1）。

a 29 671 个差异调控基因与 127 个差异代谢物的关联网络，使用 Cytoscape 软件完成网络可视化。 b 采用 400 μM 硒酸钠处理恩施碎米荠幼苗，清水作为对照，处理时长 24 h；利用液相色谱串联质谱（LCMS/MS）测定三蓟素 O丙二酰己糖苷与穗花杉双黄酮的含量。

为探究硒超富集机制，本研究使用 400 μM 硒酸钠处理恩施碎米荠幼苗两周后开展转录组测序（补充表 1）。结合基因组与转录组分析，鉴定得到硒代谢通路的 8 个关键基因：硫酸盐转运蛋白基因（SULTR）、3′磷酸腺苷5′磷酰硫酸合成酶基因（PAPSS）、腺苷酰硫酸还原酶基因（APR）、亚硫酸盐还原酶基因（SiR）、半胱氨酸合酶基因（cysK）、胱硫醚γ合酶基因（metB）、甲硫氨酸合酶基因（MetE）以及甲硫氨酸 S甲基转移酶基因（MMT）（图 4a）。硒处理后，上述基因在恩施碎米荠体内表达水平发生改变：根组织中SULTR、SiR、cysK、metB、MetE、MMT基因上调；叶片组织中SiR、APR、MetE基因上调（图 4）。

硒的富集与挥发在硒污染环境的植物修复领域备受关注，该过程可将无机硒转化为气态的二甲基硒（DMSe）\(^{[33,34]}\)。二甲基硒是植物产生的主要挥发性硒化物，其毒性约为无机硒的 1/600\(^{[35]}\)。MMT是硒挥发酶促通路中的限速酶；拟南芥中该基因突变后，硒挥发能力几乎完全丧失\(^{[36]}\)。在鉴定到的 8 个硒代谢通路基因中，MMT表达量最高（补充表 1）；该基因由全基因组复制产生的拷贝与耐硒相关的扩张metE基因家族同源（图 4）。

谷胱甘肽（GSH）巯基（SH）对金属具有高亲和性，同时作为植物螯合肽前体，是清除金属的核心物质\(^{[37]}\)；谷胱甘肽还能够缓解重金属胁迫带来的氧化损伤\(^{[38]}\)。本研究通过转录组分析鉴定出谷胱甘肽代谢通路的 5 个酶编码基因（图 5）。葡萄糖6磷酸脱氢酶（G6PD）为胞质酶，参与生成烟酰胺腺嘌呤二核苷酸磷酸（NADPH），NADPH 负责将氧化型谷胱甘肽（GSSG）还原为还原型谷胱甘肽 GSH\(^{[39]}\)。硒处理可显著上调G6PD的表达，相较对照组上调 8.5 倍，提示硒胁迫促使 NADPH 大量生成；谷胱甘肽合成酶基因（GS）同样发生表达上调。

a 已报道的硒代谢调控通路及相关基因。红色标注基因为本研究在恩施碎米荠中鉴定得到。 b 硒代谢通路基因在不同组织中的组织特异性表达。L：对照组叶片；Lse：施硒组叶片；R：对照组根；Rse：施硒组根。 c 恩施碎米荠与其他植物的metE基因家族系统发育树。 d Hi-C 互作图谱（cen039979，4 号染色体，坐标 11482047–11488608，MMT基因上下游各 2 Mb 范围），展示互作信号与 TAD 结构（40 kb 分辨率）。上方：Hi-C 互作矩阵；下方：TAD 边界（竖线）与绝缘分值。图中纵轴与蓝线代表绝缘分值，灰线标记 TAD 边界。

a 已报道的谷胱甘肽（GSH）代谢通路相关基因。红色标注基因为本研究在恩施碎米荠中鉴定得到。 b GSH 代谢通路基因在不同组织中的组织特异性表达。L：对照组叶片；Lse：施硒组叶片；R：对照组根；Rse：施硒组根。 c HiC 互作图谱（cen025484，13 号染色体，坐标 1711372317116459，GS基因上下游各 2 Mb 区间），展示染色质互作与 TAD 信号（40 kb 分辨率）。上方：HiC 互作矩阵；下方：TAD 边界（竖线）与绝缘分值。图中纵轴和蓝线代表绝缘分值，灰线标记 TAD 边界。

细胞壁是金属离子重要的储存位点，在重金属超富集与超强耐受过程中发挥关键作用\(^{[40,41]}\)。研究发现，恩施碎米荠中木质素、果胶、纤维素、葡聚糖生物合成等细胞壁代谢相关基因发生显著的基因组扩张（补充表 1）。木质素在镉耐受与镉积累过程中具有重要作用\(^{[42]}\)。此外，比较基因组分析表明木质素生物合成通路基因数量显著富集，共鉴定 30 个相关基因，包括PAL、4CL、CCR、CAD、CCOAOMT、F5H、COMT与POX（补充表 1，补充图 2）。

苯丙氨酸解氨酶（PAL）是苯丙烷途径的关键酶，可催化苯丙氨酸脱氨生成反式肉桂酸，该产物是木质素与类黄酮生物合成的共同前体\(^{[43]}\)。本研究观察到PAL基因发生明显的基因组扩张；PAL基因家族系统发育分析表明，复制基因PAL1与PAL2来源于恩施碎米荠特有的全基因组复制事件（补充图 2）。有趣的是，KEGG 分析显示，13 号染色体上受全基因组复制产生的复制基因主要富集于铁死亡、MAPK 信号通路、苯丙氨酸生物合成、硒化合物代谢、硫代谢以及类黄酮生物合成通路（补充表 1）。

DNA 甲基化与染色质重塑参与调控植物响应非生物胁迫的基因表达\(^{[44]}\)。本研究采用 Illumina 测序平台、150 bp 双端测序（150PE）模式，对硒酸钠处理第 14 天的叶片开展全基因组 DNA 甲基化分析。结果表明，硒胁迫能够影响 DNA 甲基化水平；6 号染色体对照组平均甲基化水平低于硒胁迫组（0.6066 对比 0.6803）（补充图 2）。

综上所述，本研究鉴定出恩施碎米荠特有的一次全基因组复制事件，解析了该物种的演化特征；同时通过多组学研究揭示了恩施碎米荠耐硒与硒超富集的分子机制。

材料与方法

植株培养与处理

恩施碎米荠材料取自中国西南地区恩施市湖北硒产业技术研究院。植株于温室自然光条件下培养，日间温度 20 ℃30 ℃。取一月龄叶片用于 DNA 提取；两月龄植株采用 400 μM 硒酸钠处理 24 h，清水处理作为对照，收集叶片与根组织用于 RNA 提取与代谢物定量。另取两月龄植株，使用 400 μM 硒酸钠处理两周，清水作为对照，采集根、叶组织用于 DNA、RNA 提取，后续开展转录组测序与 ChIPseq 实验。

DNA 提取与测序

选取生长状态最优的恩施碎米荠单株，采集新鲜健康叶片，液氮速冻后置于80 ℃保存，用于 DNA 提取。采用改良 CTAB 法从叶片中提取 50 μg 高质量基因组 DNA\(^{[45]}\)，加入 RNase A 去除 RNA 污染。利用 NanoDrop 2000 分光光度计检测 DNA 浓度，0.8% 琼脂糖凝胶电泳检测 DNA 完整性，样品显示高分子量单一条带；进一步采用 Femto Pulse 仪器验证 DNA 片段长度大于 30 kb，证明 DNA 完整性良好，可用于构建 Illumina HiSeq X Ten 与 PacBio Sequel 平台测序文库。 依据试剂盒说明书构建插入片段 350 bp 的 Illumina HiSeq X Ten 文库，得到 29.9 Gb 短读段数据。使用 HTQC 软件过滤低质量碱基与读段，去除接头序列、N 碱基占比＞10% 或低质量碱基（≤5）占比＞50% 的读段，最终获得 24.8 Gb（约 56×）有效数据，用于基因组调研分析以及基因组序列的碱基水平校正。

借助 BluePippin 片段筛选系统构建插入片段 20 kb 的 SMRTbell 文库；将 SMRTbell 模板在 PacBio Sequel 平台的 8 个 SMRT cell 上完成测序，得到 803 万条 subreads，总数据量 61.38 Gb，用于基因组组装。

mRNA 的 Iso-Seq 分析

采用 IsoSeq 流程获取恩施碎米荠全长转录本\(^{[46]}\)。选取与基因组测序同一株植株，混合茎、根、叶组织，使用 TRIzol 试剂（Thermo Fisher Scientific，货号 15596018）提取总 RNA。利用 BluePippin 筛选 13 kb 以及＞3 kb 片段，构建 4 个 SMRT cell 文库，在 PacBio Sequel 平台测序，产出 42.654 Gb subreads 数据。同时按照 Illumina TruSeq RNA 文库试剂盒构建转录组文库，在 Illumina HiSeq 平台开展转录组测序。

转录组与代谢组检测

按照试剂盒说明书，使用 TRIzol 试剂提取恩施碎米荠根、叶组织 RNA，构建转录组文库并在 Illumina HiSeq 2500/X 平台测序。沿用基因组组装的质控流程过滤低质量读段；使用 Trinity 软件，以默认参数组装恩施碎米荠转录本\(^{[47]}\)；RSEM 软件计算基因表达量，表达水平以 FPKM（每百万比对读段中，每千碱基基因长度对应的读段数）表示\(^{[48]}\)。 代谢组检测由武汉迈维代谢生物技术有限公司完成，采用广泛靶向代谢组技术，基于液相色谱电喷雾电离串联质谱（LCESIMS/MS）开展代谢物相对定量；利用 Cytoscape 软件构建基因代谢物关联网络。

染色质免疫共沉淀测序（ChIP-seq）

取 3 g 恩施碎米荠样本，预冷 PBS 洗涤 2 次；1% 甲醛室温交联 10 min，加入甘氨酸至终浓度 125 mmol/L 终止交联。样品裂解后于冰上获取染色质，超声破碎得到 200500 bp 的可溶性染色质片段。取 20 μL 染色质20 ℃保存，作为 Input 对照 DNA；取 100 μL 染色质，分别使用 H3K27me3 抗体（CST9733）、H3K4me2 抗体（CST9725）进行免疫沉淀。 免疫沉淀得到的 DNA 采用 NEXTFLEX® ChIPSeq 文库制备试剂盒（NOVA514120）构建测序文库；由武汉爱基百客生物科技有限公司在 Illumina X Ten 平台以 150PE 模式完成测序。

数据分析

使用 BSseeker 软件将有效读段比对至参考基因组\(^{[49]}\)；CGmapTools 统计全基因组胞嘧啶位点测序深度\(^{[50]}\)。甲基化水平计算方式：覆盖该甲基化胞嘧啶（mC）位点的读段数 ÷ 覆盖该胞嘧啶位点的全部读段数，即每个胞嘧啶位点的 mC/C 比值。利用 CGmapTools 统计不同序列环境下胞嘧啶平均甲基化水平，重新计算不同样本 mC 的分布占比\(^{[51]}\)。MethGo 软件计算各样品基因拷贝数变异\(^{[50]}\)；circlize R 包绘制基因组上甲基化位点、差异甲基化区域（DMR）以及拷贝数变异（CNV）的分布图\(^{[52]}\)。

全基因组重亚硫酸盐测序（WGBS）

取 5 g 恩施碎米荠叶片提取 DNA；合格 DNA 使用 Bioruptor 超声破碎，片段平均长度 300500 bp。采用 EZ DNA MethylationGold™试剂盒完成 DNA 亚硫酸盐转化与 PCR 扩增；Pico MethylSeq™文库试剂盒构建亚硫酸盐测序文库，合格文库片段集中在 300 bp 左右。Illumina 150PE 平台测序，测序深度达到 30×。

基于 PacBio 长读长的基因组组装

使用长读长数据进行 contig 序列组装，采用 Falcon v0.30 软件，参数为默认\(^{[53]}\)。Falcon 基因组组装分为以下步骤：首先利用 daligner 进行读段比对并生成一致性读段；随后通过 daligner 识别纠错后读段之间的重叠区域；最后基于重叠信息构建有向字符串图，通过字符串图解析得到 contig 路径。组装得到的基因组序列先使用 arrow 基于 PacBio 长读长进行校正，再利用 Pilon 软件结合 Illumina 测序数据完成二次校正\(^{[54]}\)。

原位 Hi-C 文库构建与利用 Hi-C 数据进行染色体水平组装

Hi-C 技术可用于 scaffold 构建，本研究采用该技术将 contig 锚定到染色体上\(^{[55]}\)。文库构建与 Hi-C 分析所用植株与前文为同一株恩施碎米荠，实验流程参照已发表方法\(^{[56]}\)。文库在 Illumina HiSeq X Ten 平台采用 150PE 模式测序，共得到 46.54 Gb 双端读段；经过滤后保留 45.43 Gb 数据用于后续 Hi-C 分析。使用 Bowtie 软件将双端读段分别比对至恩施碎米荠组装基因组\(^{[57]}\)。为提升有效 Hi-C 互作读段占比，采用已报道的迭代比对策略；仅保留两端均唯一比对的读段用于后续分析。参照已发表 Hi-C 流程，从比对结果中过滤自连、非连接以及其他无效读段，包括酶切位点附近起始片段、PCR 扩增产物、随机断裂片段、过大 / 过小片段以及极端片段\(^{[56]}\)。通过追踪酶切位点，计算并标准化 contig 之间的互作计数。基于 contig 互作频率矩阵对 contig 聚类，同时校正 Falcon 组装中的部分错误；对存在错误的 contig 进行拆分，得到较短的 contig。

采用 Lachesis 软件，利用凝聚层次聚类方法，基于唯一比对的双端读段对 contig 进行聚类、排序与定向，构建染色体序列\(^{[57]}\)。

基因组注释

利用 Repbase 数据库与从头构建重复序列文库对重复序列进行注释。Repbase 数据库下载自Repbase - GIRI；从头重复文库由 RepeatModeler（open−1.0.8 版）、Piler、RepeatScout、串联重复序列查找工具 TRF 以及 LTRFINDER 构建\(^{[58,59,60,61,62]}\)。使用 RepeatMasker 识别恩施碎米荠从头重复文库与 Repbase 文库中的重复元件。

蛋白编码基因注释整合同源预测、从头预测以及基于三代测序的全长转录本预测。同源预测部分：从 Ensembl 数据库下载拟南芥、甘蓝型油菜、甘蓝、碎米荠、白菜的蛋白序列，通过 TBLASTN 比对至恩施碎米荠基因组\(^{[63]}\)。使用 Augustus 与 GlimmerHMM 在屏蔽重复序列的基因组上开展基因从头预测\(^{[64,65]}\)。将 IsoSeq 获得的全长转录本通过 GMAP 比对到基因组，再用 TransDecoder 预测转录本的开放阅读框（ORF），确定候选编码序列（CDS）\(^{[66,67]}\)。最后，使用 MAKER 软件整合全部基因模型，辅以大量人工检查，共预测得到 52725 个恩施碎米荠蛋白编码基因\(^{[8]}\)。基因数量、基因长度分布、CDS 长度分布、外显子与内含子长度分布和其他物种相近（补充图 2）。通过序列同源检索完成基因功能注释，利用 BLAST2GO 在 SwissProt、TrEMBL、KEGG、InterPro 与 GO 蛋白数据库中完成蛋白功能注释\(^{[68]}\)。采用 BUSCO 评估预测得到的蛋白编码基因完整性\(^{[7]}\)。

采用多款软件与数据库注释非编码 RNA（ncRNA）：使用 tRNAscanSE，选择真核生物参数识别 tRNA 序列\(^{[69]}\)；使用 Infernal 的 cmscan 程序检索 Rfam 数据库，鉴定 miRNA、rRNA、核小 RNA 与核仁小 RNA\(^{[70,71]}\)。

基因家族演化分析

收集 9 个已测序植物物种的蛋白序列：拟南芥、琴叶拟南芥、白菜、甘蓝、萝卜、盐芥、碎米荠、荠与番木瓜。采用全对全 BLASTP 鉴定不同物种间的同源基因\(^{[72]}\)。

使用 OrthoMCL 基于序列相似性对各物种蛋白序列聚类，筛选直系同源基因\(^{[73]}\)（补充图 2）。基于基因家族聚类得到的单拷贝直系同源基因构建进化树：采用 Muscle 对每个单拷贝同源基因家族进行多序列比对\(^{[74]}\)；合并多序列比对结果，生成 phylip 格式的超基因比对矩阵。长度过滤后共保留 1124 个单拷贝直系同源基因，使用 RAxML 基于最大似然法构建系统发育树\(^{[75]}\)。

基于构建的进化树，结合 TimeTree 数据库与文献获取时间校正点；使用 PAML 软件包中的 MCMCTREE，采用贝叶斯分子钟与惩罚似然法估算物种分化时间\(^{[76]}\)。

利用 CAFÉ 软件模拟进化树上各支系的基因家族扩张与收缩事件\(^{[77]}\)。基于各单拷贝直系同源基因家族，使用 PAML 软件分支位点模型检测基因是否受到正选择；采用 Fisher 精确检验，Bonferroni 多重检验校正，FDR 阈值设为 0.05 进行显著性检验（补充图 2）。

恩施碎米荠全基因组复制事件（WGD）

通过同义替换率 Ks 估算，检测恩施碎米荠基因组中的全基因组复制事件。首先采用全对全 BLASTP（E 值＜1e−5）鉴定不同物种的同源基因\(^{[72]}\)；使用 MCSCAN 识别共线性旁系同源区块\(^{[78]}\)（补充图 2）。提取区块内所有旁系同源与直系同源基因对，利用 PAML yn00 NG 模型计算 Ks 值；根据公式 Ks/2r 计算分化时间。最终基于 Ks 分布特征评估基因组中潜在的全基因组复制事件（补充图 2）。

原位 Hi-C 文库构建

取约 2 g 植物叶片样品，2% 甲醛固定，裂解液提取，0.1% SDS 裂解；使用 200 U MboI（NEB）酶切，通过 biotin-14-dCTP 引入生物素标记的胞嘧啶核苷酸。T4 DNA 连接酶进行平末端连接；加入 200 μg/mL 蛋白酶 K（Thermo Fisher）解除交联。纯化后的 DNA 超声打断至约 400 bp。使用 Dynabeads® MyOne™ Streptavidin C1 磁珠捕获带生物素标记的连接交界处\(^{[79]}\)。采用 NEBNext® Ultra™ II DNA 文库制备试剂盒构建 Illumina Hi-C 测序文库。选取 400–600 bp 片段，在 Illumina HiSeq X Ten 平台采用 150PE 双端测序；每组样本设置两个生物学重复。

Hi-C 数据分析

基于 Hi-C 数据的染色体组装

Hi-C 文库共产出 46.54 Gb 原始双端读段；使用 BWA（bwa-0.7.17）默认参数将读段比对至校正后的基因组\(^{[79]}\)。双端读段两端分别比对至不同 contig（scaffold）的读对用于 Hi-C 辅助 scaffold 构建。利用 Lachesis 软件的凝聚层次聚类算法，将 2267 条 contig（总长 443.45 Mb）聚类为 16 个群组\(^{[80]}\)；再通过 Lachesis 对聚类后的 contig 进行排序与定向；使用 Juicebox（v1.8.8）人工校正组装错误\(^{[81]}\)。最终成功将 1089 条 contig 组装，总长 383.27 Mb，获得恩施碎米荠首个染色体水平高质量基因组；染色体长度介于 15.82 Mb 至 29.81 Mb 之间，组装序列覆盖基因组总长度 86.65%。

互作图谱构建

使用 Trimmomatic（0.38 版）质控过滤\(^{[82]}\)；叶片两组、根两组生物学重复的 Hi-C 有效读段，通过 ICE 软件迭代比对至基因组。过滤悬垂末端等无效数据，保留有效读对；使用 QuASAR-Rep（3DChromatin-ReplicateQC v0.0.1）评估生物学重复间的相关系数\(^{[83]}\)。合并每组重复数据用于后续分析。

Hi-C 图谱记录 Hi-C 实验得到的 DNA-DNA 互作信息。合并后的有效读对，按 500 kb、200 kb、100 kb、40 kb、20 kb、10 kb、5 kb 互不重叠窗口分箱，生成互作矩阵。原始 Hi-C 互作图谱存在多种系统偏差，如可比对性、GC 含量、酶切位点分布不均；本研究采用迭代归一化方法校正互作图谱，消除系统偏差。

图谱分辨率分析

Hi-C 矩阵分辨率定义：构建互作矩阵所使用的窗口大小；图谱分辨率指最小窗口，满足 80% 窗口至少存在 1000 个互作事件\(^{[29]}\)。分辨率代表能够可靠识别局部染色质结构的最小尺度。

染色质区室（Compartment）分析

染色质区室指位于同一条或不同染色体上、内部互作强度更高的结构域集合。在 200 kb 窗口热图中呈现典型棋盘格样式，高低互作交替区块对应 A、B 区室。主成分分析（PCA）第一主成分可区分两类区室；对每条染色体臂，根据第一特征向量（PC1）正负，将基因组窗口划分为 A 区室或 B 区室。A 区室为基因富集的常染色质活跃区域；B 区室为基因稀疏的异染色质沉默区域。

TAD 分析

TAD 是连续基因组区域，内部互作强度高，与相邻区域被清晰边界分隔。使用 40 kb 窗口的互作数据识别 TAD；采用绝缘分值算法，识别各样品 TAD 边界，确定 TAD 位置与数量\(^{[22]}\)。

染色体内部与染色体间互作计算

使用 Ay 团队开发的 Fit-Hi-C 软件（v1.0.1，参数 L 20000 –U 2000000 –p 2 –b200）分析各样品 10 kb 窗口的染色体内、染色体间互作，计算对应的累积概率 P 值与 FDR（q 值）\(^{[84]}\)。筛选 P 值、q 值均＜0.01 且互作计数＞2 的互作为显著互作。