单倍型解析的自交系大黄基因组组装为了解其基因组进化和蒽醌类化合物的大量积累提供了线索

The haplotype-resolved genome assembly of autotetraploid rhubarb Rheum officinale provides insights into its genome evolution and massive accumulation of anthraquinones

Plant Commun. 2024 Jan 8;5(1):100677. doi: 10.1016/j.xplc.2023.100677. Epub 2023 Aug 26.

1. 前言

药用大黄是大黄属、蓼科多年生草本植物，数千年来一直被列入许多经典的传统中药配方中，目前已在全球范围内使用。大黄根茎粗壮，茎中空且直立，基生叶大，叶片近圆形，种子宽卵形，白绿色或紫红色的小花簇生在枝条上（图1A）。药用大黄已被证明具有重要的治疗潜力，它具有抗菌、抗病毒、免疫抑制、止血和抗肿瘤活性，还能消除淤血。中国汉代（公元前 202 年至公元 220 年）的《神农本草经》中首次记载了 药用 大黄的药用价值。900 多种中成药中含有 药用 大黄，中国每年对 药用 大黄的需求量达 5500 吨。大黄还被用于预防严重急性呼吸系统综合症（SARS）、流感和肝炎病毒，甚至还被用于应对当前的 COVID-19 大流行，如 "莲花清瘟胶囊 "和 "清瘟排毒剂 "等产品。

Rheum officinale 是《中国药典》中记载的官方大黄中唯一的四倍体（2n = 4x = 44），因其有效成分丰富、品质优良、环境适应性强而被全球广泛使用。据推测，多倍体可在进化过程中掩盖有害的杂合突变，创造新的表型特征和变体，促进优势表型和物种多样性的形成，从而有助于扩大生物的生态位，提高其应对环境变化的灵活性。测序和组装技术的进步使得复杂基因组，如高度杂合、高度重复、多倍体甚至自多倍体基因组的解密成为可能。最近，已经组装了几个染色体级的自交四倍体基因组（如紫花苜蓿、甘蔗和马铃薯）和异源多倍体基因组（如小麦、花生和棉花）。以往的研究还表明，转座元件（ transposable elements， TE）在植物基因组进化、基因表达调控和环境适应方面发挥着关键作用。然而，目前仍不清楚 R. officinale 是自交系还是异交系，也不清楚多倍体化是否在其表型、代谢物生物合成和环境适应方面发挥了重要作用。

现代医学研究表明， R. officinale 富含蒽醌、蒽酮和链烯，表明它具有抗菌、抗病毒、抗炎和抗肿瘤活性，而且这些代谢物的积累在不同发育阶段和组织中存在明显差异。最近，多组学与生化分析相结合，在 Senna tora 中发现了一个与蒽醌生物合成有关的类查尔酮合成酶（CHS-L）基因。然而， R. officinale 中丰富的蒽醌类化合物的遗传基础仍然未知。最近，一些重要药材的基因组已经组装完成，为我们提供了有关三七内酯、马兜铃酸、小檗碱、皂苷和青蒿素合成的信息，从而加速了与药用植物功能基因组学相关的研究。

在本报告中，我们展示了自交四倍体草本植物药用大黄（ R. officinale ）的染色体水平、单体型解析基因组。比较基因组学分析揭示了 R. officinale 代谢物生物合成和抗逆性的基因组和基因家族进化以及同源染色体的基因组特征。我们结合转录组和代谢组数据，探讨了药用大黄不同组织中蒽醌含量差异的遗传基础。该参考基因组和多组学数据为了解大黄的进化提供了见解，并为该属植物未来的遗传研究提供了支持。

2. 材料与方法

2021年，河北大学生命科学学院将采集自云南丽江的一龄大黄置于25°C-30°C、每天14小时光照的植物温室中。用DNA Safe Plant Kit（TIANGEN）从嫩叶中分离出基因组DNA，制备CCS文库和Illumina短读长文库，分别在贝瑞基因公司（北京）的PacBio Sequel和Illumina HiSeq平台上进行测序。

用新鲜叶片制备与 DNA 交联并用甲醛固定的染色质，然后用 MboI 对交联的 DNA 进行消化。由此产生的粘性末端经生物素化和近似连接形成嵌合连接。处理后的 DNA 会进一步富集并物理剪切成 300 至 500 bp 的片段。随后，所有制备好的 DNA 片段都被处理成成对端测序文库。这些程序都是在贝瑞基因组公司（Berry Genomic Corporation）完成的。

根据制造商的说明，使用 TRIzol 试剂（Invitrogen 公司）从三个重复的大黄根、茎和叶样本中提取总 RNA。测序在 Berry Genomic 公司的 Illumina 平台上进行。

大黄的基因组大小是通过武汉杰瑞生物技术有限公司（中国四川成都）的 流式细胞仪 （BD FACSCalibur）和 k-mer 分布 分析估算的。在流式细胞仪中，根据相同比例的荧光强度和基因组大小，将番茄和玉米作为内部对照。番茄和玉米的平均荧光强度分别为 20.52 和 52.96，而大黄的平均荧光强度分别为 80.48 和 80.97。估计的 R. officinale 基因组大小分别为 6.90 和 7.04 Gb。使用 Jellyfish (v2.3.0) 和 Genome Characteristics Estimation (GCE) (v1.0.2) ，利用 Illumina 短读长创建 k-mer 集。基因组大小采用基于频率的 k-mer（k = 17）方法估算，k-mer 数量/k-mer 深度。根据 17 个 k-mer分布，估计 R. officinale 的基因组大小为 7.60 Gb (2c)，杂合率为 1.96% 。

使用 Canu (v2.1.1) ，以基因组大小 = 4g 为参数，用干净的 CCS 读数将 R. officinale 基因组组装成contigs。这样就得到了一个 7.68 Gb 的集合，N50 长度为 3.6 Mb。使用Burrows-Wheeler Aligner (BWA, 0.7.17-r1198)以默认参数将Hi-C读数映射到contigs级组装，然后用Juicer (v1.5) 进行处理。使用 ALLHiC 对这些contig进行重新排序生成 scaffold，并在 3D-DNA 的辅助下，根据染色质互作的可视化结果进行组装。根据Hi-C图谱和映射信息，使用Juicebox（https://github.com/aidenlab/Juicebox，v1.11.08）对潜在的错误组装进行人工检查和校正。使用 BWA (0.7.17-r1198)将所有 Illumina 短读数映射到已组装的染色体级基因组上，以估计组装质量。LAI 使用 LTR_retriever (v2.9.0) 计算。使用 BUSCO（v5.2.2）和 Embryophyta odb10 数据库评估了基因组组装和预测基因的准确性和完整性。

Switch error evaluation

采用 CTAB 法制备高分子量基因组 DNA，并使用基因组试剂盒进行纯化。使用武汉贝纳基因技术有限公司的牛津纳米孔高通量（ONT）测序平台生成超长读数。我们使用 N50 长度为 175 kb 的最长 10 000 个超长 ONT 读数来评估错误率。我们将每个读数分成 10 kb 的窗口，并用 BWA（0.7.17-r1198）搜索每个窗口在已组装基因组中的最佳比对结果。我们使用最近一项研究中描述的内部脚本对结果进行了处理：只有当拼接读长的长度大于其自身长度的 50%时，才将其定义为映射读 长 。如果一个读 长 中至少有三个连续窗口与其他单倍体基因组对齐，则该读 长 可能表明存在潜在的切换错误。

结合使用了 ab initio 和同源方法来识别 TE 和其他重复区域。串联重复搜索器（TRF）（v.4.07b）用于识别串联重复序列。对于从头搜索，使用了 RepeatModeler (v.1.0.11) (http://www.repeatmasker.org/RepeatModeler/)、LTR_FINDER (v1.05) 、LTRharvest (v1.5.11)和 LTR_retriever (v1.5.11) 、LTR_retriever（v1.9）用于构建从头重复文库，RepeatMasker（v.4.1.1）（http://www.repeatmasker.org/RepeatMasker/）用于基于从头文库鉴定全基因组重复序列。对于基于同源性的搜索，RepeatMasker 使用已知的重复序列库（Repbase 15.02）（https://www.girinst.org/server/RepBase/index.php）识别重复序列。

重复序列标记的基因组用于基因注释。 R. officinale 基因组中的蛋白质编码基因是通过结合转录证据、相关物种的同源性支持和非初始方法进行注释的。对于非初始预测，使用 PASA（v.2.3）利用 Trinity（v.2.12）组装的转录本预测基因结构，然后在 AUGUSTUS（v.3.2.3）中用于训练基因模型。对于基于同源性的预测，GenomeThreader（v.1.7.3）被用来搜索近缘物种的蛋白质序列。在基于转录组的预测中，使用 Trinity 将 RNA 测序（RNA-seq）数据组装成转录本，并根据转录组组装结果使用 PASA 软件进行基因结构预测。此外，使用 HISAT2（v.2.2.1）将干净的 RNA-seq 数据与基因组进行比对，并使用 StringTie（v.2.1.6）组装转录本。随后使用 TransDecoder（v.5.1.0）（https://github.com/TransDecoder/TransDecoder）对组装的转录本进行开放阅读框（ORF）预测。使用 EVidenceModeler（EVM）将上述方法预测的所有基因结构整合到非冗余基因集中。BUSCO 与 embryophyta_odb10 数据集一起用于评估最终基因集的完整性。通过整合 InterProScan (v5.48-83.0)、Nr (ftp://ftp.ncbi.nih.gov/blast/db) 和 SWISS-PROT (https://www.uniprot.org/downloads) 数据库中的注释信息，对生成的蛋白质模型进行了功能注释。使用 KofamScan（v1.3.0）和默认参数。

我们使用Infernal（v1.1.4）（http://eddylab.org/infernal/）与Rfam（https://rfam.org/）数据库和预测的rRNA、snRNA、microRNA等进行比较，检测非编码RNA。为了对 lncRNA 进行注释，根据外显子数（>1）和预测 ORF 长度（>300 bp）对组装的转录本进行筛选。然后使用 CPC (http://cpc.cbi.pku.edu.cn) 和 CNCI (https://github.com/www-bioinfo-org/CNCI) 去除具有编码潜力的转录本。我们运行了 PfamScan (https://github.com/aziele/pfam_scan)，并与 Swiss-Prot 数据库进行 BLAST 搜索，以去除已知蛋白质。最后的数据集是根据每百万映射片段中外显子的每千基片段数（FPKM）（<0.1）过滤后得到的。

从 NCBI 数据库中获得了 11 个物种（ S. oleracea、B. vulgaris、F. tataricum、A. thaliana、Theobroma cacao、Glycine max、V. vinifera、Solanum lycopersicum、Solanum tuberosum、Oryza sativa、Zea mays ）的所有蛋白质序列，并从 He 等人（2022）的研究中获得了 F. dibotrys 的所有蛋白质序列构建系统发生树。使用 OrthoFinder（v2.5.4）和默认参数确定了这些植物中的单拷贝同源基因。使用 MUSCLE（v3.8.1551）创建了氨基酸序列比对，并使用 RAxML 软件构建了系统发生树，该软件采用 PROTGAMMAIJTTF 模型和 100 次引导重复；水稻和玉米被用作外群。12 个物种之间的分歧时间是用 r8s 估算的，并从 TimeTree 网站（http://timetree.org/）获得了 6 个校正分歧时间点。基因家族的扩展和收缩由 CAFÉ 确定。使用 ClusterProfiler R 软件包进行 GO 和 KEGG 功能富集分析。

使用 MCscan（Python 版本）和默认参数分析了 R. officinale 和 F. tataricum 基因组之间的共线性。用 BLASTP（E 值 = 1e-5）对蛋白质序列进行了全序列比对。使用 MCScanX以默认参数识别具有旁系对共线性的连接区。使用 ParaAT（v2.0）以参数（-m clustalw2 -f axt）构建多个蛋白质编码 DNA 的比对，并将比对结果作为 KaKs_Calculator （v2.0）的输入数据，以默认参数计算每个旁系配对的 Ks 值。

我们使用 DupGen_finder的默认参数来识别不同的基因重复模式：WGD、串联重复（TD）、近端重复（PD；同一染色体上距离少于 10 个基因）、转位重复（转位基因重复 [TRD]）和分散重复（DSD）。

使用 iTAK（v1.6）识别蛋白质序列中的转录因子（TFs）和转录调控因子（TRs），并将单个 TFs 和 TRs 归入不同的基因家族。

根据 processDrimm  (https://github.com/xjtu-omics/processDrimm) 中的说明，OrthoFinder 被用来识别 R. officinale、F. tataricum、B. vulgaris、S. oleracea 和 V. vinifera 之间的同源基因，经过过滤后生成 DRIMM-Synteny 的输入文件。随后，使用 Python 脚本 processDrimm.py (https://github.com/xjtu-omics/processDrimm/blob/master/processDrimm.py) 处理 DRIMM-Synteny 生成的合子和区块信息，并保留符合预期拷贝数的合子区块作为 IAGS 的输入文件。最后，根据 IAGS  (https://github.com/xjtu-omics/IAGS) 的说明推断 祖先基因组结构 。

我们使用 TEsorter 对 LTR 进行分类，包括 REXdb 数据库所涵盖的 I 类和 II 类元素。然后，我们利用RT域分析关系，并用FastTree（v2.1.10）构建逆转录转座子的系统发生树。系统发生树用 iTOL (https://itol.embl.de/) 显示。TBtools (https://github.com/CJ-Chen/TBtools) 用于显示注释的基因结构。

我们选择单倍体基因组 A 作为参照。使用 Mummer (v4.0) 将每个染色体的其他三个单倍体基因组与单倍体基因组 A 对齐，参数为 -c 1000。SNP 由 show-snps 获得，参数为 -C -H -I -T -r -l。使用参数 -m -i 90 -l 100 的 delta-filter 进一步过滤原始比对结果。过滤后的 delta 文件使用 SyRI 管道以默认参数检测结构变异。根据 SyRI 输出中的序列变异定义，我们将这些变异转换成三种类型的 SV：存在/不存在变异（PAV）、倒位和易位，参考 https://schneebergerlab.github.io/syri/fileformat.html。

不耐受与耐受排序（SIFT）算法被用来预测氨基酸替换是否可能是有害的。单倍体基因组 A 被用来创建一个带有 SIFT 预测的基因组数据库（https://github.com/pauline-ng/SIFT4G_Create_Genomic_DB），单倍体基因组间的 SNP 则使用 SIFT4G_Annotator.jar 进行注释。

我们参考其他自交系基因组，利用大黄的四个单倍体基因组进行转录组分析。用 HISAT2（v2.1.1）将来自 F. tataricum 和 R. officinale 的干净读长映射到它们的基因组。使用 StringTie（v2.1.6）的 -G 选项估算基因表达量。蛋白编码基因的表达以每百万映射读数（TPM）每千碱基外显子模型的转录本进行归一化。 F. tataricum 转录组数据（根、茎和叶）是从 NCBI 下载的。

每组基因的表达水平在四个单倍体基因组中以 1:1:1:1:1 的对应关系进行比较。我们对一组中每个基因的相对表达量进行了标准化处理（1:1:1:1:1）：每个基因的总表达量除以所有基因的总表达量之和，得出每个基因介于 0 和 1 之间的数值。这一归一化表达量是参照最近的一项研究计算的三个样本的平均值。

使用 ABScript III RT Master Mix（Abclonal，RK20429）将每个组织的 RNA（1 μg）反转录为 cDNA，反应体积为 20μl。使用 2X Universal SYBR Green Fast qPCR Mix（Abclonal，RK21203）进行 RT-qPCR。反应在 StepOne Plus Real-Time PCR 系统（赛默飞世尔）上进行：95°C 预变性一个循环 30 秒，然后 95°C 5 秒和 60°C 30 秒各 40 个循环。采用 2-ΔΔCt 法计算各组织中各基因的相对表达量。

从 1 年的大黄中收集根、茎和叶组织的三个重复样品，按照 Metware 生物技术公司（武汉，中国）提供的方法提取代谢物，如前所述。简而言之，冻干样品用混合研磨机（MM 400，Retsch）以 30 Hz 的频率研磨 1.5 分钟。然后，将 50 毫克冻干粉末溶于 1.2 毫升 70% 的甲醇提取液中，每隔 30 分钟搅拌 30 秒。这一混合步骤重复六次。在 12 000 rpm 和 4°C 下离心 3 分钟后，用微孔滤膜（SCAA-104，孔径 0.22 毫米；ANPEL，中国上海）过滤提取物。样品提取物在超高效液相色谱串联质谱（UPLC-MS/MS）系统上进行分析。最后，将所有样品等量混合，制备质量控制样品（混合物），以监测分析条件的稳定性。UPLC-ESI-MS/MS 系统用于检测代谢物概况。所有代谢物的定性分析均采用中国武汉 Metware 生物技术公司建立的数据库以及其他公共代谢物数据库，包括 MassBank、KNAPSAcK、HMDB 和 METLIN。质谱数据使用 Analyst 软件 v1.6.3 和 MultiQuant 软件 v3.0.2 进行分析和定量。

我们使用了 Mfuzz R 软件包，该软件包实现了用于分析表达数据的软聚类工具。提供的分区矩阵包含完整的成员值。这些信息被用来定义由高度相关基因组成的 10 个聚类核心。

从以前的研究和 NCBI 中获得了 TPS 基因、OSC 基因和其他编码假定参与蒽醌生物合成途径的关键酶的基因。使用 BLASTP（E 值 = 1e-5）对基因组进行搜索，并选择含有 Pfam 结构域的基因作为候选基因。使用 MAFFT（v7.310）（Katoh 等人，2002 年）对多序列进行了比对，并使用 FastTree（v2.1.0）构建了最大似然树。系统发生树用 iTOL（https://itol.embl.de/）显示和注释。

如果少于两个不参与蒽醌生物合成的基因位于任何两个参与蒽醌生物合成的基因之间，则后两个基因被定义为 "接近"。三个或更多参与蒽醌生物合成的接近基因形成一个基因簇，而两个参与蒽醌生物合成的接近基因被认为是一对。这一方法参考了之前的一份报告。

通过流式细胞仪估计， R. officinale 的基因组大小为 6.90-7.04 Gb，通过 k-mer 分布分析估计为 7.60 Gb。总共获得了 131 Gb 的 PacBio CCS 读长（contig N50 为 16 kb）和 335 Gb 的 Illumina 短读长。使用 Canu 进行的 De navo 装配产生了一个 7.68 GB 的初始 contig-level 装配，N50 为 3.46 Mb（表 1），与流式细胞仪和基于 k-mer 方法估计的基因组大小一致（图 1C）（2n = 4x）。我们使用 ALLHiC 算法为自多倍体基因组构建了 allele-aware 组装，通过整合 204 Gb 的 Hi-C 数据，将组装的contigs 搭桥。在 3D-DNA 的帮助下，进一步校正了组装结果。最终的装配包含 44 条染色体中的 7.68 Gb，其中包括 11 个同源染色体组，四个等位染色体组分别命名为 A、B、C 和 D，以及 6.07 Mb 的未定位contigs（图1D）。我们通过研究 Hi-C 接触矩阵评估了组装质量，结果显示染色体组的划分非常清晰（图 1B）。四个单倍体基因组分别含有 96.20%、96.40%、96.50% 和 94.30% 的完整基准通用单拷贝同源 ( benchmarking universal single-copy ortholog， BUSCO) 基因 。

图 1. R. officinale 基因组组装和特征概述。(A) R. officinale 的种子、叶、根和植株形态。(B) 已组装染色体的 Hi-C 热图概述；每个等位基因组包含四条染色体。(C) R. officinale 基因组组装的工作流程。(D) R. officinale 的基因组特征。轨迹表示（从外到内）a、染色体（Mb）；b、GC 含量分布；c、基因密度；d、以 ChrXA 为参照的 SNP 密度；e、Gypsy 密度；f、Copia 密度；g、LTR 密度；h、同源区块。红线表示 ChrXA 和 ChrXB 之间的同源关系，蓝线表示 ChrXA 和 ChrXC 之间的同源关系，绿线表示 ChrXA 和 ChrXD 之间的同源关系，X 代表 1-11。

共鉴定出1.536-1.675 Gb的重复序列，占每个单倍体基因组的83.17%-83.64%。长末端重复（LTR）反转座子占每个单倍体基因组的64.91%-65.70%，其中Gypsy和Copia超家族分别占38.43%-39.17%和7.10%-7.13%。我们采用综合策略，包括循证方法和 ab initio 基因预测，对大黄基因组的蛋白质编码基因进行了注释。最终获得了 190 916 个蛋白质编码基因，其中 190 643 个（∼99.86%）被分配到 44 条染色体上。预测基因的平均基因长度、蛋白质序列长度和外显子数量分别为 2520 bp、310 个氨基酸和 4.4 个。在注释的基因中，179 587 个（∼94.07%）在非冗余蛋白质序列数据库（NR）、 Gene Ontology （GO）、京都基因和基因组百科全书（KEGG）、Swiss-Prot 和 InterPro 数据库中进行了功能注释。BUSCO 分析表明，这些基因的完整率达到 98.1%，表明 R. officinale 基因注释的质量很高。我们还注释了非编码 RNA 基因，在四个单倍体基因组中获得了 561 个 microRNA 基因、13 775 个转移 RNA（tRNA）基因、38 328 个核糖体 RNA（rRNA）基因、807 个核小 RNA（snRNA）基因、3426 个核小 RNA（snoRNA）基因和 6423 个长非编码 RNA（lncRNA）基因。

单倍型解析基因组组装的质量评估

为了验证单倍型解析基因组组装的准确性，我们获得了 12 Gb Oxford Nanopore ultralong读长（N50 为 101 kb）；其中 99.3% 的超长读数正确地映射到了我们组装的基因组，平均同一性为 90.2%。 Oxford Nanopore ultralong读长 映射到四个单倍体基因组的数量没有明显差异，映射比例分别为 26.40%、25.20%、24.53% 和 23.88%。我们计算出单倍体基因组的错误连接率为 1.29%，估计每 Mb 有 0.07 个切换错误。通过这种高度连续的组装，我们确定了 44 条染色体上的 58 个端粒区域（AAACCCT 或 AGGGTTT）。为了进一步评估我们的组装质量，我们将 Illumina 短读长与组装后的基因组进行了比对，结果映射率为 99.88%，组装覆盖率为 99.84%。基因组区域的读数深度分布也呈现出类似的模式；大多数区域的平均深度为 41，只有 0.28% 的 100 kb 窗口的深度大于 82。LTR 装配指数（LAI）为 18.98。这些结果证实了我们的 R. officinale 基因组组装的高质量。

比较基因组学和基因家族进化分析

我们对 R. officinale 单倍体基因组 A 和 12 个代表性植物物种的基因组进行了比较基因组学分析，以确定 R. officinale 的系统发育位置。共鉴定了 209 个单拷贝基因，并利用这些基因构建了系统发育关系。结果表明， R. officinale 与 Fagopyrum tataricum 和 Fagopyrum dibotrys 的祖先具有最密切的进化关系，它们的估计分化发生在33.4亿年前。基因家族聚类分析显示，13个物种共有8705个基因家族，2617个基因家族为 R. officinale 所独有（图2A）。功能富集分析表明，大黄中的特定基因家族主要与甜菜素生物合成、糖基转移酶、叶酸生物合成、类固醇生物合成、植物激素信号转导和植物与病原体的相互作用有关。大黄中有 2473 个扩展基因家族和 3885 个收缩基因家族（图 2A）。扩增的基因家族主要集中在根系发育；代谢的特殊方面，如苯丙类生物合成、倍半萜类和三萜类生物合成、细胞色素 P450 活性和苹果酸酶（ME）活性；以及应激反应，如对氧化应激的反应、防御反应和植物与病原体的相互作用（图 2B）。总之，大黄中的特异性和扩展基因家族主要参与特异性代谢产物的合成以及对生物和非生物胁迫的响应，这可能部分解释了为什么大黄含有丰富的药用成分，并且作为野生物种表现出很强的环境适应能力。

图 2. R. officinale 的比较基因组学和基因家族进化。(A) 基于 13 个植物物种的 209 个单拷贝直向同源物的系统发生树。右图为各物种中基因家族的分布情况。黑色数字表示每个节点的分化时间（mya，百万年前），橙色和蓝色分别表示扩展和收缩的基因家族数量。(B) R. officinale 中扩展基因家族的基因本体分析。(C) 同源直向基因的同义分歧（Ks）值柱状图。(D) R. officinale 和 F. tataricum 的同源点图。橙色线条突出显示了主要 WGD 事件的例子，表明同源关系为 2:2。示意图显示了 R. officinale 和 F. tataricum 在分化后的序列分化。(E) R. officinale 和其他三个物种中不同复制模式产生的基因数量。(F) R. officinale 中不同复制基因的功能富集分析。(G) R. officinale 中显著扩展的 ME 基因家族在四个植物物种中的系统发育推断。分支和标签根据底部的物种配色方案着色。

ME 由一个扩展基因家族编码，是一种调节植物生长和发育的关键蛋白，在植物抗逆性方面发挥着重要作用。在 R. officinale 单倍体基因组 A 中发现了 38 个 ME 基因同源物，比 F. tataricum （10 个）、 F. dibotrys （8 个）和 Beta vulgaris （5 个）多得多（图 2G）。我们发现 39.47% 的 ME 基因（38 个基因中的 15 个）在所有三个测试组织（根、茎和叶）中都有表达，其中 9 个基因在根中高表达，4 个基因在茎中高表达，2 个基因在叶中高表达（图 2G）。

通过计算同义替换率（Ks）的密度分布，研究了大黄（ R. officinale ）进化过程中的全基因组重复（WGD）事件，发现了两个显著的峰值。第一个Ks峰为∼2.21，代表了古老的γ三重复制事件，与最近在荞麦和黄花蒿中的发现一致；第二个Ks峰为∼0.89，反映了蓼科植物在66.5 mya时的另一个WGD事件（图2C）。我们在 R. officinale 、 F. tataricum 和 Vitis vinifera 之间发现了许多 2:2:1 的同源关系，并发现大多数重复的片段在进化过程中经历了大规模的染色体重排、基因缺失或假基因化，导致 R. officinale 中的同源区块支离破碎。幸运的是，我们在 R. officinale 和 F. tataricum 之间发现了几个具有 2:2 同源关系的大型保留区段，它们可能在进化过程中得到了高度保守（图 2D）。对 R. officinale 单倍体基因组 A 的比对分析表明，4 号染色体（Chr4）上的一个区域（26.8 Mb）与 Chr8 上的另一个区域（29.4 Mb）存在明显的同源区，Chr2 上的一个区域（42.2 Mb）与 Chr5 上的另一个区域（29.8 Mb）存在明显的同源区；同源区分别包含 693 和 833 个基因对（图 2D）。所有单倍体基因组都保留了这些片段。GO 和 KEGG 富集分析表明，配对基因主要与蛋白质二聚化活性、聚半乳糖醛酸酶活性、肌醇磷酸代谢以及戊糖和葡萄糖醛酸的相互转化有关。WGD 基因的保留、缺失和分化可能导致了 R. officinale 的新表型和多样化。

基因复制为进化提供了原始遗传材料，并为物种进化提供了潜在的新基因。我们在 R. officinale A 基因组中发现了 6619 个分散重复基因（DDGs）、2073 个近端重复基因（PDGs）、2068 个串联重复基因（TDGs）、15842 个转座重复基因（TRDGs）和 17260 个 WGD 基因（WGDGs）。 R. officinale 的四个单倍体基因组的基因重复率相似。值得注意的是， R. officinale 的 DDGs 和 TRDGs 数量高于 F. dibotrys 、 F. tataricum 和 B. vulgaris （图 2E）。在 TRDGs 和 DDGs 中发现了更高的 Ka/Ks 比值和更小的 Ks 值，这表明与其他复制模式产生的基因相比，这些基因的序列分化更快，正向选择更强，是一个持续不断的过程。功能富集表明，WGDGs在植物与病原体相互作用、丝裂原活化蛋白激酶（MAPK）信号通路和转录因子中明显富集，这可能提高了大黄对生物和非生物胁迫的抗性。TDGs主要富集于苯丙类化合物的生物合成、甜菜碱类化合物的生物合成和黄酮类化合物的生物合成中（补图 10D），它们可能在大黄药用物质的生物合成中发挥关键作用。TRDGs主要富集于细胞对DNA损伤刺激的反应、根系发育和囊泡介导的运输中，表明它们在维持基因组稳定、运输特定分泌蛋白和根茎发育中发挥作用（图2F）。值得注意的是，包括 AP2/ERF 和 bZIP 家族在内的许多转录因子在欧芹中都得到了扩增，主要是通过 WGD。据报道，AP2/ERF 和 bZIP 基因家族参与了植物对非生物胁迫的反应以及药用植物主要活性成分的生物合成。因此， R. officinale 中扩大的重复基因可能对特定物质的合成和强大的环境适应能力非常重要。

通过构建 R. officinale 、 F. tataricum 、 B. vulgaris 、 Spinacia oleracea 和 V. vinifera 的祖先核型，研究了蓼科和藜科植物染色体的进化历史。结果表明，蓼科植物和藜科植物的共同祖先祖先1（Ancestor1）含有8条染色体。随后，它经历了 45 次染色体裂解和 47 次染色体融合，形成了含有 6 条染色体的祖先 3。祖先 3 经历了一次 WGD 事件，形成了 12 条染色体。79 次 染色体裂解和80 次 染色体融合形成了祖先3的11条 R. officinale 染色体，125 次 染色体裂解和129 次 染色体融合形成了 F. tataricum 的8条染色体（图3A）。这些结果有助于我们了解蓼科植物染色体的进化历史。此外，我们还发现，在 R. officinale 和 F. tataricum 分化之后，发生了广泛的染色体重排。有趣的是，这些重排经常发生在富含 TE 的区域，根据这些区域的 GC 含量、基因密度和 TE 覆盖率，我们推断这些区域为候选的同源染色体区（图 3A）。因此，我们推测染色体断裂和重排可能是 TEs 的转座或重组潜能的结果。

图 5. R. officinale 的代谢物和基因表达图谱。(A) 圆图显示了 R. officinale 的代谢物类别组成。每种颜色代表一种代谢物类别，色块的面积表示其所占比例。(B) 不同类别次生代谢物中每种代谢物的平均含量。(C) R. officinale 根、茎和叶中不同醌类化合物含量的层次聚类。(D) R. officinale 中第 4 组和第 9 组的功能富集分析。(E) R. officinale 根、茎和叶中不同萜类化合物含量的层次聚类。(F) R. officinale 和 A. thaliana 的 TPS 基因的系统发育推断。标签根据右下方的物种颜色方案着色。最外侧的轨迹代表五个类别：TPS-a、TPS-b、TPS-c、TPS-e/f 和 TPS-g。(G) 根、茎和叶中 TPS 和 OSC 基因表达的层次聚类。

在已鉴定的代谢物中，醌类化合物的平均含量最高（图 5B）。有趣的是，我们发现 62 种醌类化合物中有 51 种是蒽醌类化合物，它们在根中的含量远高于在茎和叶中的含量（图 5C）。萜类化合物的含量在不同组织中也有显著差异；叶中的单萜类化合物含量更高，而根和茎中的三萜类化合物含量更高（图 5E）。

通过对所有表达基因的共表达分析， 确定了 10 个共表达簇 （补图 27A）。簇 3、5 和 7 中的基因在根中的表达量高于在茎和叶中的表达量。功能富集分析表明，这些基因富集于磺基转移酶活性、微管马达活性、微管结合、O-甲基转移酶活性、UDP-糖基转移酶活性和防御反应中（补充图 27B）。第 4 和第 9 组基因在叶片中的表达量高于根和茎。功能富集分析表明，簇 4 中的基因主要富集于萜烯合成酶活性和单萜、倍半萜和三萜类化合物的生物合成，而簇 9 中的基因主要与光合作用有关（图 5D）。

4. 讨论

药用植物可合成多种特殊代谢物，在治疗疾病和维护人类健康方面发挥着重要作用。由于欧当归是一种具有治疗潜力的珍贵药材，因此探索其基因组进化及其活性物质的生物合成机制非常重要。

在这项研究中，我们组装了一个染色体组规模的、单体型解析的大黄自交系基因组。比较基因组学分析表明， R. officinale 的特异基因家族和扩增基因家族主要参与特殊代谢产物的生物合成以及对生物和非生物胁迫的响应。TEs 对 R. officinale 的基因组进化、基因拷贝数变异以及代谢物生物合成、抗逆性和根系发育相关基因的调控和表达做出了重要贡献。我们发现了许多基因具有等位基因差异和等位基因特异性表达，这可能是由于在其调控区域插入了不同的 TE。转录组和代谢组分析表明，参与蒽醌途径的基因显著扩增以及四倍体导致的基因数量增加，可能是大黄蒽醌含量高的原因之一。此外，与蒽醌生物合成有关的重要基因的高表达导致根中的蒽醌含量远高于其他组织。

分子技术的最新进展推动了对植物多倍体和多倍体基因组进化的研究。多倍体可产生新的表型性状和大量染色体级变体，其中一些变体因有助于多倍体物种的环境适应而被选择固定下来。多倍体还能引起基因表达和基因产物丰度的全基因组变化。大黄（ R. officinale ）是一种最近被证实的加倍自多倍体，其基因组将对未来研究其与二倍体大黄在基因组进化、表型和药用物质组成方面的差异具有重要价值。

我们的研究为了解药用植物，尤其是多倍体植物代谢途径的基因组基础提供了深入的见解。我们在大黄（ R. officinale ）中发现了几个可能与代谢物生物合成和根系发育有关的基因，这些基因通过转座复制经历了显著的扩增。我们假设 Scaffold_103.6 可能是 R. officinale 根茎发育的关键基因，RT-qPCR 验证了该基因在根中的最高表达量。对这些基因的进一步研究可能对提高大黄的质量和产量以及该物种的育种计划具有重要意义。我们的研究结果表明，转座复制产生的重要基因数量的增加以及这些基因的表达可为大黄的表型创新提供遗传基础。

众所周知， 代谢物可能由相邻基因编码 。我们鉴定了参与蒽醌生物合成途径的基因对和基因簇，发现基因簇往往位于 R. officinale 基因组的 Chr7 和 Chr10 上，而基因对往往位于 Chr1 和 Chr7 上。这些基因簇或基因对可能是蒽醌生物合成的重要候选基因。特别是 Chr7 上的基因可能是代谢工程和有益分子异源重组的明显候选基因。

目前已公布的自交系基因组包括重要的饲料作物紫花苜蓿、重要的糖料作物甘蔗、世界四大粮食作物之一的马铃薯以及中国著名的香料辣椒。在这项研究之前，还没有关于 药用植物自交系基因组的报道 。通过组装自交系药用植物基因组，我们的研究为基础研究和鉴定与农艺性状相关的基因提供了前所未有的基因组资源，并可指导大黄的利用、改良和多倍体育种。考虑到大黄中丰富的 TE 含量及其对基因拷贝数和表达的重要影响，本研究还为详细研究 TE 多态性在长期进化史中的功能效应和动态活动提供了重要参考。

高颜值免费 SCI 在线绘图 ( 点击图片直达 )

最全植物基因组数据库IMP ( 点击图片直达 )

往期精品 ( 点击图片直达文字对应教程 )

机器学习