凤丹牡丹巨型染色体与超大基因组--文献精读236

Genomic basis of the giga-chromosomes and giga-genome of tree peony Paeonia ostii

凤丹牡丹巨型染色体与超大基因组的基因组学基础

摘要

凤丹（Paeonia ostii）是我国原产、具备重要经济价值的观赏牡丹。其种子富含 α- 亚麻酸（ALA）等不饱和脂肪酸，油脂利用价值突出。本研究完成了凤丹染色体水平基因组组装，基因组大小 12.28 Gb。与大基因组单子叶植物不同，凤丹物种演化历程中未发生物种特异性全基因组复制事件；基因间区长末端重复序列（LTR）在约 200 万年的短期内爆发式扩增，是其超大基因组形成的关键诱因。此外，五类组蛋白编码基因发生扩增，在巨型染色体结构维持过程中发挥重要作用。

本研究基于 448 份种质开展全基因组关联分析（GWAS），证实脂肪酸合成通路关键节点基因（SAD、FAD2、FAD3）发生基因扩张且表达量显著上调，是凤丹种子高效富集 α- 亚麻酸的分子基础。通过与栽培牡丹（Paeonia suffruticosa）基因组对比发现，A 类花器官基因AP1异位表达、C 类基因AG表达下调，共同介导雄蕊瓣化性状形成。本研究发布的基因组数据可为牡丹染色体与基因组演化解析、新品种选育提供宝贵资源。

引言

凤丹（Paeonia ostii T. Hong et J.X. Zhang）为芍药科多年生中型灌木，原产我国华中地区。其花大色艳，在中国历史与传统文化中占据重要地位，素有 “花中之王” 的美誉，人工栽培史约 2000 年，至今栽培广泛，也是现代栽培牡丹（Paeonia suffruticosa）各类品种的原始亲本。相较于多数被子植物，凤丹染色体数目偏少（2n=10）、染色体体积偏大（长度 10–15 μm）、基因组巨型（＞12 Gb）。受巨型染色体自身繁育劣势，加上长期采挖根皮入药等因素影响，凤丹野生种群濒临灭绝。目前仅有栽培牡丹的碎片化基因组草图被发表，凤丹可用基因组资源十分匮乏。

形态层面，凤丹雄蕊数量繁多，雄蕊群呈离心式发育；从叶片到苞片、花萼、花冠器官呈连续螺旋式排列；胚胎发育阶段存在游离核分裂特征，上述性状与部分裸子植物相似。栽培牡丹品种普遍存在雄蕊瓣化现象，但野生原种凤丹无该性状，目前栽培牡丹雄蕊瓣化的分子调控机制尚不明确。

凤丹籽油不饱和脂肪酸含量超 90%，其中富含人体自身无法合成、必需从外界摄取的 α- 亚麻酸（ALA）。已有研究证实陆生植物脂肪酸合成通路进化保守，但凤丹种子为何能大量富集亚油酸与 α- 亚麻酸，相关分子机理仍不清晰。

本研究完成凤丹染色体水平基因组组装，证实该物种基因组大小突破 10 Gb，进化过程未发生被子植物类群特异性全基因组复制（WGD）；解析了巨型染色体形成与结构维持的潜在调控因素，筛选出不饱和脂肪酸合成关键候选基因，同时阐释了栽培牡丹雄蕊瓣化形成的潜在成因。

结果

凤丹基因组组装与注释

为构建高质量凤丹基因组，试验取材于河南洛阳原生凤丹植株。流式细胞术预估基因组大小为 12.76 Gb（附图 1a、1b）。测序共获得 2.97 Tb Illumina 干净短读长数据（测序深度 247×）、643.67 Gb PacBio 长读长数据（53.6×）以及 2.50 Tb Hi-C 测序数据（附表 1–3）。采用 SOAPdenovo 联合 wtdbg 开展从头组装，组装结果经 Pilon 纠错优化（附图 1c）；最终组装基因组总长 12.28 Gb，contig N50 为 228 kb，scaffold N50 达 2.43 Mb（表 1、附表 4–7）。

约 11.49 Gb（占全基因组 93.5%）的序列被锚定至 5 条染色体上，与细胞学鉴定结果吻合，等位序列间重叠区域极少（附表 8、附图 1d）。凤丹基因组平均 GC 含量为 32.8%（附图 1e）。采用 BUSCO 通用单拷贝同源基因评估组装完整度，组装完整基因占比 94.4%；结合转录组数据高比对覆盖率（附表 10–12），共同证实本次基因组组装质量优异（附表 9）。

多手段联合注释共预测得到 73177 个蛋白编码基因，其中锚定在染色体上的高可信度基因 54451 个，基因平均长度 8.22 kb（附表 13、附图 2a）。依托公共数据库开展功能注释，59768 个（84.53%）基因可匹配已知功能蛋白（附表 14）；BUSCO 评估显示 85.5% 的预测基因结构完整（附表 15）。受超大基因组特征影响，凤丹基因组共注释得到 330511 个假基因。本研究同时完成凤丹叶绿体与线粒体基因组的组装和注释（附图 2b、2c、图 3、图 4）。

凤丹基因组共鉴定出 15238 个基因家族，在陆生植物中处于基因家族数量前列（附图 5a–d）。异黄酮等次生代谢物合成相关基因家族在凤丹基因组中显著富集（附表 16、附图 5e）；物种特异性扩张基因主要富集于 DNA 复制、脂肪酸代谢通路（附表 17、附图 5f）。研究还注释得到 42131 条非编码 RNA（附表 18）。最终将全部注释信息定位到凤丹 5 条假染色体上，各染色体长度介于 1.78～2.56 Gb 之间（图 1）。

a 假染色体圈图；b～d 依次为 CCGG 甲基化密度、CCWGG 甲基化密度、H3K27me3 修饰富集峰密度分布（统计窗口 10 Mb）；e～h 依次为 Copia、Gypsy、Del 三类转座子密度与基因密度分布（统计窗口 20 Mb）；i 共线性区块，连线代表筛选后的旁系同源基因共线性关联。

超大基因组演化

为解析凤丹演化历程，将其基因组与已发表的 12 个物种基因组开展比较分析（附表 19）。基于 13 个物种共 474 个单拷贝同源基因构建系统发育树，结果表明凤丹归于虎耳草目（附图 6a–d），与以往分子标记、质体基因组及转录组研究结论一致；佐证虎耳草目与葡萄目亲缘近，为蔷薇类的姊妹类群。分子钟推算芍药科与景天科分化时间约 1.09 亿年前（1.02~1.208 亿年前），和依托rbcL基因得到的 1.1 亿年分化时间相符。

对凤丹共线性区块内基因进行 ** 四简并位点颠换值（4DTv）** 分析：凤丹进化历程经历一次近期 LTR 转座子爆发峰值，同时发生真双子叶植物共同的 γ 全基因组复制事件（约 1.3 亿年前，图 2a、附图 6e）。在被子植物中凤丹基因组＞10 Gb 却未发生物种特有全基因组复制，该特征十分特殊（附表 20）。

a 4DTv 值分布，反映凤丹基因组发生一次近期 LTR 转座子爆发扩增与真双子叶共有 γ 全基因组复制事件；b 不同分类基因中 Ks＜0.2 的基因占比；c 选取的裸子植物、双子叶、单子叶物种全长 LTR 转座子亚类分类统计；d 凤丹基因组 LTR 整体及各亚类插入时间估算；e LTR 反转录转座子复制相关 5 种关键酶编码基因的富集丰度；f 凤丹、拟南芥、落地生根、葡萄基因间区相对长度，数据以平均值 ± 标准误表示；g 转座子总覆盖度、Gypsy 与 Del 亚型覆盖度在基因本体区及基因上下游 2 kb 侧翼区的分布；h 牡丹基因间区、基因上下游 2 kb 侧翼区内，表达基因、沉默基因、假基因对应的 CCGG、CCWGG 位点平均甲基化密度。

转座子（TE）是真核基因组的普遍组成元件。凤丹中侧翼序列与转座子发生重叠的基因里，83.51% 基因的 Ks＜0.2，说明转座子在巨型染色体与超大基因组形成过程中起到关键作用（图 2b、附表 21、22）。LTR 转座子占比最高（43.75%），其次为 DNA 转座子与长散在重复序列（LINE）；Gypsy 是丰度最高的 LTR 超家族（附表 23）。对比 16 个所选植物物种可见，转座子亚家族组成差异极大，其中 Gypsy 下属 Del 亚家族在凤丹中发生剧烈扩张，占全长 LTR 总量的 70.58%（图 2c、补充数据 1），证明 Gypsy、Copia 不同亚型在各类基因组中的扩增贡献差异显著。小麦、银杏的超大基因组同样经历过 LTR 爆发式扩增。凤丹 LTR 大规模插入发生在距今 100 万～200 万年前，插入峰值在 140 万年前（图 2b、d），该时间与 Gypsy、Copia 的插入窗口期吻合（附表 24）。进一步分析显示，参与 Del 元件合成的 5 种关键酶基因富集度，是其他类型 LTR 合成相关酶的 8.25 倍，这是凤丹 Del 元件爆发扩增的重要原因（图 2e）。

凤丹基因间区总长是葡萄、拟南芥等小基因组物种的 15 倍，但 mRNA、编码区（CDS）、外显子 / 内含子等基因本体区段长度在各物种间无明显差异（图 2f、附图 6f）。由此可知，大量 LTR 插入造成基因间区不断复制扩张，是凤丹染色体与基因组巨型化的核心驱动力。转座子及其 Gypsy、Del 亚型的基因组分布结果显示，转座子绝大多数插入富集在基因间区（图 2g），进一步佐证 LTR 爆发是基因组与染色体膨大的关键诱因。全基因组 MethylRAD 甲基化测序表明，凤丹表达基因、非表达基因与假基因在基因区、基因间区呈现截然不同的甲基化修饰模式（图 2h）；表达基因在 CCGG 位点甲基化富集、CCWGG 位点无明显富集，说明甲基化修饰调控基因表达类型。

为探究蛋白DNA 互作关系，本研究开展 ChIPseq 测序。结果显示：相较于基因本体区，转座子区段旁系同源基因的启动子区显著富集 H3K27me3 修饰；且非表达基因启动子的 H3K27me3 富集程度高于表达基因（附图 6g）。上述结果说明，巨型染色体与超大基因组并未破坏绝大多数功能基因的结构与正常表达。

巨型染色体的形成演化

为解析凤丹古演化历程，将其基因组与桃、拟南芥、葡萄开展比较基因组分析。重构得到包含 7343 个原始基因的 7 条真双子叶祖先染色体，筛选出凤丹与祖先基因组间 171 个共线性区块（图 3a）。依托重建的两套祖先核型（ARK，n=7、n=21）推算：凤丹从祖先核型演化至现存 5 条染色体，至少经历4 次染色体断裂、20 次染色体融合。凤丹现存核型演化过程未发生额外多倍化事件；而拟南芥（现存 5 条染色体 = 21 条祖先染色体 + 85 次断裂−101 次融合）在从十字花科祖先分化时，经历了更为复杂的多轮基因组加倍事件（图 3a）。

a 从真双子叶祖先到凤丹及其他植物的染色体重排与结构演化历程。染色体以不同颜色区分，用以追溯片段源自含 7 条原始染色体的重构祖先核型（ARK，n=7、n=21）；染色体长度做 100 倍缩放处理。图中标注全基因组复制（WGD）事件，物种分化预估时间标注于进化分支上。b 五类组蛋白（H1、H2A、H2B、H3、H4）基因家族扩张情况（**p＜0.01，***p＜0.001）。c H2、H3 组蛋白基因扩张与表达模式热图，标尺颜色代表基因相对表达量。

染色体主要由 DNA 与组蛋白构成。本研究发现，与其他植物相比，凤丹中编码五类核心组蛋白（H1、H2A、H2B、H3、H4，维系染色质结构稳定的关键组分）的基因家族发生显著扩张，据此推测该现象与凤丹巨型染色体的结构维持密切相关；同样拥有大染色体的玉米也出现 H2A 组蛋白基因扩张的特征（图 3b、附图 7a–e、补充数据 2）。凤丹基因组共注释得到 208 个组蛋白基因，均匀分布在 5 条巨型染色体上（附图 7f、补充数据 2）。其中，促进染色质浓缩的H2A.W、参与染色体组装与 DNA 复制的关键基因H3.1，相较于同亚家族其他亚型基因表现出基因扩张且表达量显著上调（图 3c）。已被证实参与染色体组装、基因沉默调控的 H3.1 蛋白第 90 位突变位点存在特异性 G 变异，该变异在凤丹及另外 4 种植物中被检出（附图 7g）。有意思的是，38 个H3.1旁系同源基因里，有 36 个基因侧翼序列与转座子区域重叠（附表 25）。

脂肪酸性状全基因组关联分析（GWAS）

为从群体水平解析凤丹种子高富集多不饱和脂肪酸（PUFA）的分子机制，本研究对源自国内不同产区的 448 份凤丹种质开展全基因组重测序（附图 8a）。采用 SLAF-seq 简化基因组测序，共获得约 1.34 Tb 原始测序数据，单样本标签测序深度约 12×（补充数据 3）。

围绕脂肪酸合成相关 14 个性状（各类脂肪酸含量表型等）开展全基因组关联分析（附表 26、补充数据 4）。C18:0/C18:1△9、C18:1△9/C18:2△9,12 是脂肪酸合成通路两个关键中间产物比值，与上述两个性状显著关联的 SNP 位点分别定位在 2 号、3 号染色体区段（图 4a）。3 号染色体一段约 20.02 Mb 区间（655617 k～675640 k，先导 SNP 位于 673208760）密集分布大量显著性 SNP，与 C18:0/C18:1△9 性状强关联（图 4b、附图 8b）；该区间注释到 9 个SAD同源基因，SAD 酶负责将 18:0-ACP 转化生成 18:1-ACP（图 4c、附图 8c、补充数据 5）。其中Pos.gene65901定位于 900 kb 连锁不平衡区块内唯一关键拷贝，转录组与荧光定量 PCR（RT-qPCR）共同验证该基因高表达；该区段核苷酸多态性 π 值偏低，暗示人工选育过程中该位点受到定向选择（附图 8d）。2 号染色体上候选基因 * FAD2（Pos.gene24209）* 催化 18:1 – 磷脂生成 18:2 – 磷脂，与 C18:1△9/C18:2△9,12 性状高度连锁（图 4d、e，附图 8e、f）；其余候选基因详见补充数据 6。

a、d 分别为 C18:0/C18:1△9、C18:1△9/C18:2△9,12 性状的全基因组关联分析曼哈顿图；标注出与本研究已鉴定SAD、FAD2基因重合的显著关联峰位点。b、e 从上至下依次为 3 号、2 号染色体显著峰值区间的局部曼哈顿图、连锁不平衡（LD）热图。c 3 号染色体小段区间内成簇排布的 9 个SAD基因。

内质网定位型 α- 亚麻酸（ALA）生物合成

借助质谱检测对凤丹种子脂肪酸含量进行动态测定，必需脂肪酸（尤以 α- 亚麻酸为主）在胚乳成熟期积累量达到峰值（图 5a）。不同发育时期转录组分析显示，差异表达基因显著富集于不饱和脂肪酸生物合成通路（附表 27）。该通路多个关键节点基因发生家族扩张且表达量显著上调，包括SAD（硬脂酰 – ACP 去饱和酶，13 个拷贝）、FAD2（脂肪酸去饱和酶 2，4 个拷贝）、KAS1（酮脂酰 – ACP 合酶 1，3 个拷贝）与 FAD3（脂肪酸去饱和酶 3，4 个拷贝）（图 5b、附图 9a~f、补充数据 7），是凤丹种子高富集 α- 亚麻酸的重要分子基础。

a 牡丹种子发育过程连续五年 α- 亚麻酸（ALA）平均含量变化；DAF 为受精后天数。红色标记样品已开展转录组测序，数据以平均值 ± 标准误表示。b 脂肪酸合成通路简图并标注候选基因；热图依据凤丹转录组数据展示候选基因在 6 个发育阶段（受精后 35、49、63、77、91、119 天）的表达变化，关键主基因标红突出。c 凤丹及代表性植物 ω-3 去饱和酶系统发育树，区分内质网定位FAD3、质体定位FAD7/FAD8两大亚家族。

被子植物普遍合成亚油酸、α- 亚麻酸（ALA）等必需脂肪酸；进化上，这类脂肪酸先在质体内合成前体，再转运至光面内质网（ER）完成组装。亚油酸向 α- 亚麻酸转化的关键反应由 ω-3 脂肪酸 Δ15 去饱和酶（ω-3 FAD）催化，该酶分为两类：一类定位于质体的FAD7/8，主要在叶片等绿色组织表达，不参与种子油脂中 α- 亚麻酸合成；另一类为种子内质网定位的FAD3，是种子油脂 α- 亚麻酸合成的关键功能基因。目前 ω-3 脂肪酸 Δ15 去饱和酶基因家族的演化规律尚不明确。被子植物祖先经历一次全基因组复制（WGD）事件，推动该基因家族后续发生功能分化与亚细胞定位分化（图 5c）。系统发育结果显示，FAD7/8与FAD3在单子叶、双子叶各类陆生植物中广泛分布（图 5c、补充数据 8）。已有研究表明，相较于双子叶植物，玉米等大宗粮食作物在内的绝大多数单子叶种子 α- 亚麻酸含量偏低；而牡丹、亚麻等双子叶作物依靠种子发育特定阶段内质网型FAD3上调表达，实现多不饱和脂肪酸（PUFA）大量积累。本研究克隆 3 个FAD3基因开展功能验证，亚细胞定位证实其蛋白定位于内质网（附图 10a）；酵母异源表达实验证明FAD3_4（Pos.gene25040）为优势高表达功能拷贝（附图 10b、附表 28）。综上，脂肪酸合成通路关键节点基因家族扩张、内质网型FAD3高表达，共同驱动凤丹种子大量积累 α- 亚麻酸。

雄蕊瓣化花器官的形成机制

野生原种凤丹由叶原基分化形成花原基后，花原基进一步分化出心皮、雄蕊、花瓣、花萼、苞片与叶片，天然不存在雄蕊瓣化性状（附图 11a）；但栽培牡丹（P. suffruticosa）品种花芽分化过程中，部分雄蕊可转变为花瓣，约 60% 品种在现蕾期、透色期出现该性状（图 6a）。凤丹与栽培牡丹比较进化分析表明，AP1、AG等 MADS-box 家族多拷贝同源基因是花器官建成、雄蕊瓣化形成的关键调控基因（图 6b、附图 11b）。两类牡丹花发育相关基因整体表达模式相近（图 6b、附图 11c），符合经典 ABCE 花发育模型：A 类基因AP1在花瓣高表达，C 类基因AG在雄蕊高表达；在瓣化雄蕊组织中，AP1出现异位低量表达，AG表达水平相较正常雄蕊小幅下调（图 6b）。由此说明：雄蕊中AP1异位表达叠加AG表达下调，诱发花器官同源异型转换，雄蕊发育为瓣状结构。该性状增大花冠体积，提升栽培牡丹观赏价值。

上述表达特征完善了经典 ABCE 模型在核心双子叶中的演化规律：拟南芥等模式植物严谨 ABCE 模型由花菱草等基部双子叶 “边界渐变（fading borders）” 演化而来；凤丹作为栽培牡丹的野生祖先种遵循严谨 ABCE 调控模式，而栽培牡丹保留 “边界渐变” 调控特征。两类牡丹 ABCE 家族基因进化速率存在差异，人工驯化进程中栽培牡丹 A 类基因受正向选择、C 类基因受纯化负选择（附图 11d）。A、C 类基因表达的相互拮抗发生改变，是凤丹与栽培牡丹花器官分化差异、栽培牡丹形成雄蕊瓣化的核心诱因（图 6c）。利用荧光定量 PCR 检测AP1、AG在两类牡丹不同组织表达量，验证该调控模型（图 6d、附图 12a–e）：A 类基因AP1（Pos.gene28418）在花瓣高表达、瓣化雄蕊低表达、正常雄蕊几乎不表达；C 类基因AG（Pos.gene74744）表达趋势完全相反（图 6d）。该结果证实AP1异位表达与AG下调直接促成雄蕊瓣化，也为牡丹瓣花新品种定向选育提供理论依据。

a 凤丹（P. ostii）与 4 个栽培牡丹品种（从上至下：二乔、乌龙捧盛、玉楼春、白雪塔）分别在现蕾期、透色期、盛花期的花器官形态（从左至右）；比例尺：1 cm。b 花发育相关基因亚家族扩张与表达量热图，取样时期为现蕾期、透色期。Po：凤丹；EQ：二乔；WL：乌龙捧盛；YLC：玉楼春；BXT：白雪塔；st_pe：瓣化雄蕊；色标代表基因相对表达水平。c 凤丹与栽培牡丹 ABCE 花发育调控模式示意图。d 荧光定量 PCR 检测 A 类AP1、C 类AG基因在花瓣、瓣化雄蕊、正常雄蕊中的相对表达量；每个基因设置 4 次生物学重复，数据以平均值 ± 标准误表示。

讨论

自然界中拥有超大基因组的植物种类较多，但兼具超大基因组与巨型染色体的物种十分罕见。祖先核型重构结果表明，凤丹从祖先核型演化至现存 5 条巨型染色体，演化历程未发生物种特异性多倍化，仅依靠至少 4 次染色体断裂、20 次染色体融合完成核型演变。本研究证实，凤丹巨型染色体形成主要由基因间区长末端反转录转座子（LTR）大规模扩张驱动。以往巨型染色体研究多局限于细胞学观察，其基因组层面的形成与维持机制尚不清晰。本研究结果表明，凤丹单条巨型染色体大小介于 1.78~2.56 Gb，五类组蛋白基因家族扩张与巨型染色体结构维持密切相关；后续仍需借助组蛋白突变体功能试验，进一步解析组蛋白在巨型染色体基因组稳定性维系中的具体作用。

植物超大基因组的形成主要存在两大驱动途径：转座子扩张与全基因组复制（WGD）。小麦、大蒜等单子叶超大基因组物种依靠全基因组复制实现基因组膨大，而凤丹未发生物种特异性全基因组复制，在超大基因组被子植物中演化模式独特；其基因组与染色体巨型化源于基因间区转座子大量插入扩张，且驱动Del型转座子增殖的 5 类关键酶编码基因优势富集，进一步助推Del亚型爆发扩增。

经典 ABCE 模型阐释了双子叶植物花器官分化发育规律：A 类基因单独调控花萼形成，A+B 类基因共同调控花瓣发育，B+C 类基因协同决定雄蕊分化，C 类基因单独控制心皮生成。本研究证实栽培牡丹雄蕊瓣化由AP1与AG的表达模式改变主导，为牡丹及其他观赏花卉瓣花新品种选育提供候选靶标基因。依据 ABCE 理论，单独 A/E 类基因决定花萼（绿色保护性花器官），B+C 基因共同调控雄蕊发育。本研究发现，AP1异位表达、AG表达下调，加之驯化过程中 A、C 类基因分别受到正向选择与纯化选择，使得野生凤丹遵循严谨的 ABCE 调控模型，而栽培牡丹呈现 “边界渐变型（fading borders）” 花器官调控特征。

基因组与转录组数据明确了 ω3 脂肪酸去饱和酶基因家族（FAD3/FAD7/8）的复制事件。本研究系统解析凤丹脂肪酸保守合成通路全部基因，证实通路各关键节点至少 1 个基因拷贝呈高表达；SAD、FAD2、FAD3等关键基因发生家族扩张且优势高表达，是凤丹种子高富集 α亚麻酸（ALA）的分子基础。后续可通过基因编辑或转基因技术，将内质网定位型 FAD 基因导入谷物作物，提升粮食作物多不饱和脂肪酸合成能力。

材料与方法

植物材料与基因组测序

试验用河南洛阳原生凤丹种植于上海辰山植物园，取材用于全基因组从头测序。采集幼叶、幼芽等组织，采用天根 DNAsecure 植物基因组提取试剂盒提取基因组 DNA。构建 61 个插入片段梯度文库用于 Illumina 双端测序，原始数据经 SOAPnuke（v1.5.5）过滤质控；构建 20 kb 大片段文库，依托 PacBio RSII 与 Sequel 平台开展单分子实时测序（SMRT-seq），保留长度＞500 bp 的测序读长用于后续组装。HiC 文库采用四碱基内切酶MboI酶切基因组 DNA，酶切片段连接后进行双端高通量测序。

基因组组装

采用德国 PARTEC CyFlow 流式细胞仪，以烟草基因组（约 4.5 G）为内标测定凤丹基因组大小；利用 GenomeScope（v1.0）开展 Kmer 分析。Illumina 短读序列采用 SOAPdenovo（v2.04）设定参数63mer -K 57 -z 20000000000 -R -M 1 -k 31 -F完成初步组装；PacBio 长读序列使用 wtdbg（v1.2.8）-t 50 -i PEO.fa.gz –tidy-reads 5000 -fodbg -k 0 -p 21 -S 4 –rescue-low-cov-edges组装得到原始 contig。通过 BLASR、SPARC、BWA 将 Illumina 数据比对至 PacBio 校正后的 contig，Pilon（v1.23）纠错，SSPACE（v2.1.1，参数-x 0 -m 32 -o 20 -z 0 -k 5 -a 0.7 -n 15 -v 0）利用大片段文库挂载得到 scaffold。组装质量从三方面验证：随机 Illumina reads 用 BWA 比对、RNAseq 数据用 HISAT2 比对、最终 scaffold 利用 BUSCO（v3.0.2）评估完整度。HiC 数据经 JUICER、3DDNA 挂载 contig 至染色体，JuiceBox 可视化手动微调得到伪染色体。

基因组注释

重复序列注释

结合从头预测 + 同源比对两套策略注释重复序列：LTR_FINDER、LTRharvest、LTRdigest、RepeatModeler 构建物种特异性从头重复序列库；RepeatMasker、RepeatProteinMask 比对 Repbase 数据库注释已知转座子，两套结果整合，TRF 注释串联重复序列。

基因结构注释

整合同源蛋白比对、从头预测、转录本辅助预测三种证据集：

高可信度 / 低可信度基因划分标准

基因功能注释

基因蛋白序列 BLASTP（E＜1e-10）比对 SwissProt/TrEMBL、KOG、NR 数据库，优选最优匹配注释基因功能；InterProScan 检索 Pfam、PANTHER、SMART 等结构域数据库注释保守基序与结构域；依托 NR 注释映射 GO 条目，基因序列比对 KEGG 完成通路注释。

假基因注释

参考蛋白经 Exonerate 比对屏蔽重复后的基因组，查询序列比对覆盖度＞70%、存在移码突变或提前终止密码子判定为假基因；最终共注释 330511 个假基因，其中 77372 个源自高可信度基因，253139 个源自低可信度基因。

非编码 RNA 注释

tRNA：tRNAscan-SE（v1.3.1）真核生物参数预测；rRNA：拟南芥 5S/5.8S/18S、水稻 28S rRNA 为模板 BLASTN（E＜1e-5）检索；miRNA、snRNA：INFERNAL（v1.1.2）比对 Rfam 数据库预测。

超大基因组进化分析

下载 12 个物种基因组注释基因集用于比较进化分析：拟南芥（TAIR10）、无油樟、耧斗菜、番木瓜、葡萄、番茄、莲、麻风树、猕猴桃、落地生根、向日葵、大花红景天。

基因家族聚类：蛋白序列两两 BLASTP 比对（E＜1e-5），OrthoMCL (v2.0.9) 采用默认参数（膨胀系数 1.5）完成基因家族分簇；单拷贝同源家族蛋白经 MUSCLE (v3.8.31) 多序列比对。提取各物种单拷贝基因四简并位点（4DTv）与第一密码子位点序列，首尾拼接构成超级基因用于系统发育构建；trimAl (v1.4，-gappyout) 修剪拼接后的第一位点序列。基于修剪后的序列，RAxML 结合 GTRGAMMA 模型构建最大似然进化树；同时采用 PhyML 重构进化树进行结果互证。依托全部单拷贝基因第一位点，使用 PAML 软件包 MCMCTree 的马尔可夫链蒙特卡罗贝叶斯算法估算物种分化时间与进化速率；CAFÉ(v2.1) 解析基因家族扩张与收缩。共线性分析：牡丹、落地生根、葡萄基因集 BLASTP 比对，MCscanX（-k 50 -s 5 -e 1e−05 -m 25）检索基因组共线性区块。

为解析牡丹全基因组复制事件与基因重复，经 BLASTP (E＜1e-5) 筛选双向最优匹配（RBH）旁系同源基因；PAML (v4.8) 中 yn00 程序采用 Nei–Gojobori 算法计算配对基因 Ks 值；筛选 Ks＜0.2 的基因对，开展转座子插入位置与侧翼序列关联分析。

计算 13 个物种单拷贝同源基因的 Ka/Ks：MUSCLE 比对编码序列，PAML codeml 自由比率模型估算各分支 Ka、Ks 与 Ka/Ks；配对 Wilcoxon 秩和检验比较牡丹与其余物种间同源基因平均 Ka/Ks 差异。针对牡丹支上 Ka/Ks＞1 的候选基因，采用分支位点模型开展正选择验证，先后进行 M1a 模型 vs 分支位点模型、分支位点模型（model=2, NSsites=2）vs 空模型（固定 ω=1）两次似然比检验（LRT）；贝叶斯经验贝叶斯（BEB）筛选正选择位点，卡方检验（df=1，p＜0.05）判定显著性。

祖先染色体重构与组蛋白拷贝数统计

以葡萄基因组为参考，BLASTP (E＜1e−5) 筛选牡丹同源基因，每个基因保留最优 10 条比对结果；依托同源基因，MCScanX 检索葡萄基因组共线性区块，结合葡萄染色体信息重构真双子叶 7 条祖先原始染色体，筛选牡丹与祖先基因组间共线性区段；同步完成桃、葡萄、拟南芥与祖先核型的共线性比对。

选取 46 种植物（2 种裸子植物、2 种基部被子植物、10 种单子叶、32 种双子叶）统计各类组蛋白基因拷贝；拟南芥、玉米、银杏、牡丹组蛋白各亚型序列经 CLUSTALW (v2.1) 比对，MEGA7.0 构建最大似然进化树，自展值设 500。

脂肪酸含量测定与全基因组关联分析（GWAS）

2014 年将源自洛阳、铜陵、亳州、菏泽、邵阳 5 个产区的 448 份凤丹种质引种至上海辰山植物园芍药资源圃（北纬 31°4′52″，东经 121°10′14″）。2016–2019 连续 4 个生长季，每份材料自花人工授粉，逐年单株收种；种子 60℃烘干至恒重，液氮研磨成粉。称取 0.2 g 干粉装入螺口玻璃管，加入 3 mL 氯仿：甲醇（V/V=1:2），35℃水浴 120 rpm 振荡萃取 1 h；补加 1 mL 氯仿混匀，再加 1.8 mL 纯水，终溶剂比例氯仿：甲醇：水 = 1:1:0.9，4000×g 离心 15 min，吸取下层氯仿相；重复萃取 2 次，合并有机相，氮吹浓缩，−20℃避光保存。

油脂样品溶于 2 mL 4% 硫酸甲醇溶液，充氮密封，涡旋 1 min 后 90℃水浴甲酯化 1 h；冷却后加 1 mL 纯水 + 1 mL 正己烷，混匀后 4000×g 离心 15 min，吸取上清，氮吹定容，4℃待测。每份样品添加 20 μL 50 mg/mL 十九烷酸正己烷溶液作为内标。采用安捷伦 78905975 气质联用仪、HP88 毛细管柱（60 m×0.25 mm，膜厚 0.2 μm）检测；柱温程序：70℃保持 1 min，10℃/min 升至 210℃，继续 10℃/min 至 220℃，再 10℃/min 升至 235℃保温 8 min；进样口 250℃，分流比 5:1，进样量 1 μL；载气高纯氦气，流速 1 mL/min，EI 电离 70 eV。依托 NIST 质谱库、37 种脂肪酸甲酯混标（Sigma）定性；内标标准品外标法定量，最小二乘法建立标准曲线；脂肪酸含量以 mg/g 干重表示，每份样品 3 次生物学重复。单因素方差分析（ANOVA，p＜0.05）结合 Tukey 多重比较分析脂肪酸含量差异。

对 448 份嫩叶构建 SLAF 简化基因组文库测序，BWA (v0.7.17) 比对参考基因组；samtools+bcftools 调用 SNP，保留质控得分＞20 的变异位点，剔除次等位基因频率 MAF＜0.05、缺失率＞25% 的低质量 SNP；SHAPEIT (v4.0) 完成 SNP 填充与单倍型定相。基于定相 SNP 计算 p距离，PHYLIP 构建邻接（NJ）进化树；GCTA 开展主成分分析（PCA），ADMIXTURE 解析群体结构。

剔除性状偏离均值 ±2 倍标准差的极端离群样本；TASSEL (v5.0) 计算个体亲缘矩阵；利用 GAPIT (v3.0) 的 FarmCPU 模型开展 GWAS，PCA.total=4，分别采用多年均值、分年份表型数据关联分析；LDBlockshow 筛选连锁不平衡区块（-SeleVar 1 -BlockType 3 -MerMinSNPNum 3 -BlockCut 0.7:0.8）；自编脚本以 100 kb 无重叠滑窗统计全基因组核苷酸多态性 π 值。

种子发育期脂肪酸合成通路基因转录组分析

人工授粉后，2016–2019 连续 4 个生长季每周（7 d）从 7 株凤丹单株采收种子，共划分 13 个发育时期（图 5a），选取其中 6 个时期（红色标注样本）开展转录组测序。种子液氮速冻后−80℃低温保存。采用 Omega E.Z.N.A. HP 植物 RNA 提取试剂盒提取 6 个时期（受精后 35、49、63、77、91、119 d）共 12 份样品总 RNA，经 Qiagen RNeasy 植物微量试剂盒纯化；Agilent 2100 检测 RNA 浓度与完整性，所有样品 OD₂₆₀/OD₂₈₀介于 2.0~2.1、RNA 完整值 RIN＞7.0。总 RNA 经 DNase I 除 DNA 污染，Oligo (dT) 富集 mRNA，反转录合成 cDNA；纯化短片段、末端修复、加 A 尾、连接测序接头，筛选合适片段 PCR 富集建库；Agilent 2100 与 ABI StepOnePlus 实时荧光定量仪质控文库，文库在深圳华大依托 Illumina HiSeq 4000 平台上机测序。HISAT 将 clean reads 比对至凤丹参考基因组，用于差异基因筛选与不饱和脂肪酸合成通路富集分析。

脂肪酸去饱和酶编码基因进化分析

从 NCBI 下载 31 种代表性绿色植物基因组数据，BLAST 检索油脂合成关键基因：SAD（催化 C18:0→C18:1）、ω6 型FAD2/FAD6（C18:1→C18:2）、ω3 型FAD3/FAD7/FAD8（C18:2→C18:3）。31 个物种 FAD 家族氨基酸序列经 MUSCLE v3.8.31 多序列比对，PhyML v3.0 构建最大似然进化树，iTOL 在线工具可视化进化树。

α- 亚麻酸合成关键基因荧光定量 PCR（qRT-PCR）验证

筛选 8 个参与脂肪酸合成与三酰甘油组装、重点关联 ALA 合成的关键基因：FAD3_4、FAD7/8、SAD_3、FAD2_1、FAD6、CALO_5、STERO_4、OLE_1，检测其在种子 6 个发育时期（35/49/63/77/91/119 DAF，共 12 份样品）的表达量。以Actin为内参基因，ABI StepOnePlus 平台上机，每个基因 4 次生物学重复。Omega 试剂盒提取各组织总 RNA，天根 FastKing 反转录试剂盒（带 gDNA 去除酶）取 1 μg 总 RNA 合成第一链 cDNA。20 μL 反应体系：TB Green Premix Ex Taq II 10 μL、上下游引物（10 μM）各 0.8 μL、cDNA 模板 0.3 μL、无酶水 8.1 μL；扩增程序：95℃预变性 30 s；40 个循环（95 ℃ 5 s，64 ℃ 30 s）；95℃ 15 s，熔解曲线验证引物特异性。采用 2⁻ΔΔCt 法计算基因相对表达量。

FAD3 亚细胞定位与体外功能验证

亚油酸去饱和酶 FAD3、FAD7/8 是催化亚油酸（LA）生成 α- 亚麻酸（ALA）的限速酶，其中 FAD7/8 定位于质体，FAD3 定位于内质网。亚细胞定位试验：构建 FAD 与绿色荧光蛋白 GFP 融合载体，冻融法转入农杆菌，侵染本氏烟草叶片，荧光观察亚细胞位置。酵母异源表达：依托 Gateway 克隆体系，构建FAD3_1、FAD3_2、FAD3_4酵母表达载体 pDonr207、Pyes-DEST52（载体由上海辰山植物园赵卿研究员馈赠），转入尿嘧啶缺陷型酿酒酵母 INVSc1（菌株由上海海洋大学周志刚教授馈赠）；重组菌株 pY31（FAD3_1）、pY33（FAD3_3）、pY34（FAD3_4）在缺尿嘧啶 SC-U 筛选培养基（2% 葡萄糖为碳源）筛选阳性转化子；野生型与转基因酵母外源添加底物亚油酸，2% 半乳糖诱导基因表达；参照前述 GC-MS 方法检测酵母脂肪酸组分。

ChIP-seq 与全基因组甲基化测序

委托上海云序生物开展 3 组生物学重复 H3K27me3 ChIP-seq，使用 H3K27 三甲基化特异性抗体富集染色质片段；Qunat-IT 荧光定量试剂盒测定富集 DNA 总量，qPCR 验证富集效率；NEBNext 试剂盒构建 Illumina 文库，Agilent 2100 质控后 HiSeq 平台 150 bp 双端测序。选取 2 组凤丹种子生物学样本开展 MethylRAD 简化甲基化测序，解析全基因组 CpG、CHG 甲基化位点分布及其在基因功能元件上的定位；II 型限制性内切酶 FspEI 酶切建库，Illumina 2000 平台 100~150 bp 双端测序。

花器官性状转录解析

采集野生原种凤丹、4 个栽培牡丹品种（二乔、乌龙捧盛、玉楼春、白雪塔）现蕾期、透色期花芽；同一植株选取发育均一花芽置于冰上，徕卡体视显微镜下分离花瓣、瓣化雄蕊、正常雄蕊、露色期花瓣，每类组织 4 组生物学重复；液氮速冻 20~30 min 后−80℃保存，用于 RNA 提取、测序与荧光定量验证。HMMER 基于保守 HMM 模型全基因组鉴定 MADS-box 家族基因；PhyML 构建两类牡丹 ABCDE 花发育基因进化树；PAML 分支位点模型（runmode=0，fix_omega=0）计算各进化分支 ω 值（替换速率比）。检测 A/B/C/D/E 类花发育基因在花瓣、瓣化雄蕊、正常雄蕊、初绽花瓣 4 种组织的差异表达，解析现蕾、透色两时期基因亚家族扩张与表达特征；选取栽培牡丹与野生凤丹共 56 份组织样品开展 qRT-PCR 验证，每个基因设置 4 次重复。