五种肉苁蓉属植物叶绿体基因组-文献精读25

Structural mutations of small single copy (SSC) region in the plastid genomes of five Cistanche species and inter-species identification

五种肉苁蓉属植物叶绿体基因组中小单拷贝 (SSC) 区域的结构突变及物种间鉴定

摘要

背景

肉苁蓉属是列当科的重要属类，具有重要的药用、经济和防治荒漠化的价值。然而，肉苁蓉属的系统发育关系尚不清晰。迄今为止，还没有有效的分子标记能够有效区分报道的肉苁蓉属近缘种类。在本研究中，我们获得并表征了来自中国的四种肉苁蓉属植物的叶绿体基因组，以澄清该属内的系统发育关系，并开发用于物种鉴定的分子标记。

结果

四种肉苁蓉属植物（肉苁蓉、盐生肉苁蓉、管花肉苁蓉和中华肉苁蓉）通过Illumina深度测序。使用SPAdes组装它们的叶绿体基因组，并使用CPGAVAS2进行注释。详细分析了叶绿体基因组，发现它们都显示了保守的四分结构（LSC-IR-SSC-IR），完整大小范围从75到111 Kbp。我们观察到小单拷贝区域（SSC，范围为0.4–29 Kbp）的显著收缩和倒位重复区域（IR，范围为6–30 Kbp）的扩展，肉苁蓉和盐生肉苁蓉显示最小的SSC，仅有一个基因（rpl32）。与其他列当科植物相比，肉苁蓉属植物表现出极高的基因丢失和假基因化率，正如其他寄生列当科植物所报告的那样。此外，对蛋白质编码基因的序列变异性分析表明，在研究的肉苁蓉属植物中，三个基因（rpl22, clpP和ycf2）经历了正向选择。此外，通过比较所有可用的肉苁蓉叶绿体基因组，我们发现了25个高度变异的间隔区（IGS）区域，并使用这些区域预测了两个DNA条形码标记（基于IGS区域trnR-ACG-trnN-GUU的Cis-mk01和Cis-mk02）和十一种特异性DNA条形码标记。实验验证表明，两种类型的标记在物种鉴定中的成功率为100%。

结论

我们的研究结果表明，肉苁蓉属植物是研究寄生植物叶绿体基因组进化过程中结构变化、基因丢失和假基因化的理想模型，为肉苁蓉属植物的进化关系提供了新的见解。此外，开发的DNA条形码标记可以实现正确的物种鉴定，确保肉苁蓉属植物作为药用产品的有效和安全使用。

背景

肉苁蓉属植物是一类不进行光合作用的寄生植物，属于列当科。该科包含非寄生、半寄生和全寄生植物。列当科是研究植物寄生进化和表型可塑性的优秀模型系统。作为不进行光合作用的全寄生植物，肉苁蓉属完全没有叶绿素和功能性叶片。它们外观奇特，包含从裸露地面冒出的巨大、颜色鲜艳的花穗，因而被称为“沙漠风信子”。约有20至30种公认的肉苁蓉属植物，分布从西方（马卡罗尼西亚）到东方（中国西北）。肉苁蓉属植物生长在沙漠，有时也生长在沿海沙丘或盐沼中，寄生在各种耐盐灌木的根部。肉苁蓉属植物是用于固定沙地和防治荒漠化的重要工业栽培植物资源。

大多数肉苁蓉属植物传统上被用作滋补和药用草药，特别是被誉为“沙漠人参”的肉苁蓉。多项研究报道了肉苁蓉属植物的化学成分和药理作用。从该属植物中分离出100多种化合物，包括苯乙烯类、苷类、碳水化合物、木脂素、环烯醚萜类、紫草素、毛蕊花苷、绿原酸、毛蕊花苷、黄芩苷等。这些分离出的化合物表现出有趣的药理作用，如神经保护、免疫调节、抗衰老、抗炎、抗骨质疏松、护肝、抗氧化、抗菌、抗肿瘤和改善葡萄糖耐量的作用。由于这些药用特性和经济利益，人们对肉苁蓉属植物进行了过度开发和消费，这使得这些植物在本研究中具有重要性。

迄今为止，中国已报道了四种肉苁蓉属植物，即肉苁蓉、盐生肉苁蓉、中华肉苁蓉和管花肉苁蓉。然而，由于以下原因，肉苁蓉的分类仍有争议：首先，现有的凭证标本太少；其次，相关文献数量稀少；第三，新鲜和干燥标本中花的颜色差异显著，导致基于花色的物种分类混乱；第四，很难确定寄主植物的类型。Najibeh Ataeia等人报告了大规模取样的肉苁蓉属植物的分子系统发育分析，表明根据地理分布存在四个分支。然而，这项研究没有解决种间鉴定的问题。此外，最具经济重要性的物种肉苁蓉未被包括在内。因此，我们旨在收集中国的肉苁蓉属植物，并通过系统发育分析确定它们的进化地位。

叶绿体是绿色植物特有的细胞器，在光合作用中起着重要作用，并具有自己的基因组。叶绿体基因组高度保守，包括基因组大小、结构、基因内容和组织。因此，叶绿体基因组是系统发育分析、遗传多样性评估和分子鉴定的优秀工具。近年来，完整的叶绿体基因组序列已成功用作植物超级条形码，以区分某些类群中的近缘种，例如芍药属。而在其他情况下，叶绿体衍生的DNA标记已被开发用于鉴定药用植物，例如在人参属植物叶绿体基因组中的SNP或插入-缺失突变。

在我们的研究过程中，四种肉苁蓉属植物的叶绿体基因组序列已被发布。讨论了肉苁蓉属植物叶绿体基因组中的严重基因丢失和假基因化，但仅作为补充材料展示。在这项工作中，我们首次使用五种完整的肉苁蓉属植物叶绿体基因组来：（1）探讨肉苁蓉属植物叶绿体基因组内的结构变异；（2）阐明肉苁蓉属植物中的基因丢失和假基因化现象；（3）开发特定分类群的分子标记和DNA条形码，以区分肉苁蓉属植物。这里获得的结果提高了对肉苁蓉属植物之间关系的理解，并将对确保肉苁蓉属植物药用产品的有效和安全使用具有重要价值。

结果

肉苁蓉叶绿体基因组的测序和组装

Illumina读取序列映射到组装的叶绿体基因组上，获得了100-300倍的平均覆盖深度（图S1）。四个叶绿体基因组显示出典型的圆形四分结构，并在基因组织和结构上表现出高度保守性（总结见图S2）。四分结构由一个大单拷贝（LSC）区域（32,470–52,005 bp）和一个小单拷贝（SSC）区域（398–29,719 bp）组成，两个倒位重复（IR）区域（6,593–30,352 bp）将其分开（表S1）。肉苁蓉（图S3）、盐生肉苁蓉（图S4）、管花肉苁蓉（图S5）和沙苁蓉（图S6）的叶绿体基因组示意图显示出来。肉苁蓉、盐生肉苁蓉、沙苁蓉和管花肉苁蓉的叶绿体基因组大小分别为109,454 bp、111,690 bp、111,500 bp和75,735 bp（表S1）。肉苁蓉、盐生肉苁蓉、沙苁蓉和管花肉苁蓉叶绿体基因组的总体GC含量分别为36.27%、36.11%、36.75%和34.95%（表S1）。总体而言，除管花肉苁蓉外，SSC区域的GC含量低于LSC和IR区域。五个叶绿体基因组的基因内容列在表S1中。所有基因的数量，包括蛋白质编码基因（PCGs）、假定假基因、tRNA基因和rRNA基因的预测以及可能丢失的基因数量（表S1）。总体而言，可能丢失了29-44个PCGs。丢失的基因数量与这些叶绿体基因组中识别的PCGs数量相似，表明肉苁蓉叶绿体基因组存在大规模的基因丢失和基因侵蚀。

肉苁蓉叶绿体基因组中SSC区域的扩展和收缩

与其他列当科植物的叶绿体基因组相比，肉苁蓉叶绿体基因组的SSC区域显著收缩，大小范围从仅27 bp到61,091 bp（图1A，表S2）。相反，IR区域显著扩展，大小范围从2,318 bp到45,796 bp，而其他15种列当科植物的IR区域长度超过25 kb（表S1，表S2）。值得注意的是，肉苁蓉（398 bp）和盐生肉苁蓉（435 bp）的SSC区域极短，只包含rpl32基因（图1B）。管花肉苁蓉的SSC区域最长，约为29,719 bp。管花肉苁蓉的SSC区域中有八个PCGs，包括ycf2、ycf15、rps7、rps12、rpl23、rpl32、ycf1和rps15（表S4）。这些观察结果表明，肉苁蓉的SSC区域是基因丢失、假基因化和重排的热点（图S10）。SSC区域的大小变化主要由ndh、rps和ycf基因的丢失和假基因化引起（表S3）。此外，由于IR的扩展/收缩，肉苁蓉叶绿体基因组中IR/LSC和IR/SSC的连接高度可变（图2）。肉苁蓉和盐生肉苁蓉的LSC/IRa（IRb）区域的连接高度保守，rpl32基因位于SSC区域（表S4）。但rpl32基因与边界的距离有所不同（图2）。与管花肉苁蓉相比，盐生肉苁蓉和肉苁蓉中rpl32和rps15基因的顺序相反。另一个有趣的观察是，ycf1基因在沙苁蓉和C. phelypea的IRb和SSC区域中都有存在（图2）。如图S7所示，这些肉苁蓉属植物的叶绿体基因组保守性较低。许多基因丢失，包括rpoA、rpoB、rpoC、rpoC2、psaB、psaD、psaI、petA、petB、petD、cemA等（表1和表S5）。此外，与肉苁蓉相比，在盐生肉苁蓉、沙苁蓉、C. phelypaea和管花肉苁蓉的叶绿体基因组中发现了一些大规模倒位（图S8）。相比非寄生的列当科植物地黄（R. glutinosa），五种肉苁蓉属植物显示出不同程度的基因组丢失（图S9）。

36种列当科植物叶绿体基因组中SSC区域的变异。A）GS、LSC、SSC和IR区域的长度绘制在最大似然树（ML树）的右侧。B）SSC区域中的蛋白质编码基因。SSC区域的示意图绘制在最大似然树（ML树）的右侧。

五种肉苁蓉属植物倒位重复区扩展和收缩的分析。轨道顶部的基因逆时针转录，而轨道下方的基因顺时针转录。肉苁蓉属植物LSC（浅蓝色）、IRa（黄色）、SSC（绿色）和IRb（黄色）区域边界的示意图显示在左侧。LSC和IRb（JLB）、IRb和SSC（JSB）、SSC和IRa（JSA）、IRa和LSC（JLA）之间的连接位点用垂直虚线标记。基因rpl22（深蓝色）、rps19（浅红色）、ycf1（蓝色）、ycf2（深蓝色）、rps15（浅橙色）、rpl32（浅绿色）、rpl2（绿色）显示在叶绿体基因组上方。基因特征上方的数字表示基因边界与连接位点之间的距离（起始或结束）。

Photosynthesis and energy production genes						Ribosomal RNA genes							Transfer RNA genes
	Cd	Cs	Cc	Cp	Ct		Cd		Cs	Cc	Cp	Ct		Cd	Cs	Cc	Cp	Ct
atpA	Ψ	Ψ	Ψ	○	○	rrn16	•		•	•	•	•	trnA-UGC	•	•	•	•	○
atpB	Ψ	Ψ	Ψ	○	○	rrn23	•		•	•	•	•	trnC-GCA	•	•	•	•	•
atpE	Ψ	Ψ	Ψ	○	○	rrn4.5	•		•	•	•	•	trnD-GUC	•	•	•	•	•
atpF	Ψ	Ψ	○	○	○	rrn5	•		•	•	•	•	trnE-UUC	•	•	•	•	•
atpH	○	Ψ	○	○	○								trnF-GAA	•	•	•	•	•
atpI	○	Ψ	Ψ	○	○								trnG-GCC	•	•	•	•	•
ndhA	○	○	○	○	○								trnG-UCC	○	○	○	•	○
ndhB	Ψ	Ψ	Ψ	Ψ	Ψ								trnH-GUG	•	•	•	•	•
ndhC	○	○	○	○	○								trnI-CAU	○	○	○	•	○
ndhD	○	○	○	○	○	RNA polymerase and intron maturase genes							trnI-GAU	○	○	•	•	•
ndhE	○	○	○	○	○		Cd	Cs		Cc	Cp	Ct	trnK-UUU	•	•	•	•	•
ndhF	○	○	○	○	○	matK	•	•		•	•	•	trnL-CAA	•	•	•	•	•
ndhG	○	○	○	○	○	rpoA	Ψ	Ψ		Ψ	○	○	trnL-UAA	•	•	•	•	•
ndhH	Ψ	Ψ	○	Ψ	○	rpoB	Ψ	Ψ		Ψ	Ψ	Ψ	trnL-UAG	•	•	•	•	•
ndhI	○	○	○	○	○	rpoC1	○	○		○	○	○	trnM-CAU	•	•	•	•	•
ndhJ	○	○	Ψ	○	○	rpoC2	Ψ	Ψ		Ψ	Ψ	○	trnN-GUU	•	•	•	•	•
ndhK	○	○	Ψ	○	○								trnP-UGG	•	•	•	•	•
petA	○	○	○	○	○								trnQ-UUG	•	•	•	•	•
petB	○	○	Ψ	○	○								trnR-ACG	•	•	•	•	•
petD	○	○	Ψ	○	○								trnR-UCU	•	•	•	•	○
petG	Ψ	Ψ	Ψ	○	○	Ribosomal protein and initiation factor genes							trnS-CGA	•	•	○	○	○
petL	○	○	○	○	○		Cd	Cs		Cc	Cp	Ct	trnS-GCU	•	•	•	•	•
petN	○	○	○	○	○	infA	•	•		•	•	•	trnS-GGA	•	•	•	•	•
psaA	Ψ	Ψ	Ψ	Ψ	Ψ	rpl14	•	•		•	•	•	trnS-UGA	•	•	•	•	•
psaB	Ψ	Ψ	Ψ	Ψ	Ψ	rpl16	•	•		•	•	•	trnT-GGU	•	•	•	•	•
psaC	○	○	○	○	○	rpl2	•	•		•	•	•	trnT-UGU	•	•	•	•	•
psaI	○	○	○	○	○	rpl20	•	•		•	•	•	trnV-GAC	•	•	•	•	•
psaJ	○	○	○	•	○	rpl22	•	•		•	•	•	trnW-CCA	•	•	•	•	•
psbA	Ψ	Ψ	Ψ	Ψ	Ψ	rpl23	•	•		•	Ψ	•	trnY-GUA	•	•	•	•	•
psbB	Ψ	Ψ	○	○	○	rpl32	•	•		•	○	•
psbC	Ψ	Ψ	○	○	Ψ	rpl33	•	•		•	•	•
psbD	Ψ	Ψ	○	○	○	rpl36	•	•		•	•	•
psbE	Ψ	○	○	Ψ	Ψ	rps11	•	•		•	•	•
psbF	Ψ	Ψ	○	○	○	rps12	•	•		•	•	•	Other essential genes
psbH	○	○	○	○	○	rps14	•	•		•	•	•		Cd	Cs	Cc	Cp	Ct
psbI	○	○	Ψ	Ψ	○	rps15	•	•		•	•	•	accD	•	•	•	•	•
psbJ	○	○	○	Ψ	○	rps16	•	•		•	•	○	ccsA	○	○	○	Ψ	○
psbK	○	○	○	Ψ	○	rps18	•	•		•	•	•	cemA	○	○	○	○	○
psbL	Ψ	Ψ	○	○	○	rps19	•	•		•	•	•	clpP	•	•	•	•	•
psbM	Ψ	Ψ	Ψ	○	○	rps2	•	•		•	•	•	ycf1	•	•	•	Ψ	•
psbN	○	○	○	○	○	rps3	•	•		•	•	•	ycf15	•	•	Ψ	○	•
psbT	○	○	○	○	○	rps4	•	•		•	•	•	ycf2	•	•	•	•	•
psbZ	○	○	○	Ψ	○	rps7	•	•		•	•	•	ycf3	Ψ	Ψ	○	○	○
rbcL	Ψ	Ψ	Ψ	Ψ	Ψ	rps8	•	•		•	○	•	ycf4	○	○	○	Ψ	Ψ

Cd 肉苁蓉, Cs 盐生肉苁蓉, Cc 沙苁蓉, Cp C. phelypaea, Ct: 管花肉苁蓉.•: genes present; Ψ: pseudogene; ○: genes missing

基因丢失和假基因化在肉苁蓉叶绿体基因组中

SSC区域的收缩和IR区域的扩展伴随着基因丢失和假基因化。我们基于34种列当科植物的基因丢失和假基因比例构建了一个层次聚类树，并使用拟南芥和烟草作为外类群（图3）。列当科包括九个属，如列当、肉苁蓉或铁线草等。与拟南芥相比，除地黄外，大多数列当科植物的叶绿体基因组中都存在基因丢失（表S5）。值得注意的是，列当科植物的叶绿体基因组特征是高比例的基因丢失，范围从41.25%到75%。丢失的基因主要集中在光合作用和能量生产基因组，如ndh或atp基因。在肉苁蓉叶绿体基因组中，C. phelypaea和管花肉苁蓉的基因丢失比例最大，而肉苁蓉和盐生肉苁蓉的基因丢失比例最小（图3A，表1和表S5）。

推测的36种列当科植物的基因丢失和假基因化。树枝长度与丢失基因或假基因的最大数量成正比。热图显示的是丢失基因与假基因的比例。层次聚类树使用dendro进行。A) 基于基因丢失比例的层次聚类树。B) 基于假基因比例的层次聚类树。

Pedicularis ishidoyana、Orobanche austrohispanica、O. densiflora、O. pancicii、O. rapum-genistae、Lindenbergia philippensis、Rehmannia glutinosa和R. solanifolia的所有叶绿体基因都是功能性的。假基因最常见于Lathraea squamaria和Aphyllon californicum，其比例分别为41.25%和37.5%。在肉苁蓉属植物的叶绿体基因组中，假基因的比例从0.16到0.29不等，依次为管花肉苁蓉（0.16）、沙苁蓉（0.19）、肉苁蓉（0.23）、C. phelypaea（0.28）、盐生肉苁蓉（0.29）（图3B和表S6）。肉苁蓉属植物中的大多数光合作用和能量生产基因都是假基因，包括ndh、psa、psb、rps、rpl基因（表S6）。

肉苁蓉叶绿体基因组中PCGs的核苷酸置换率

dN/dS比率（ω = dN/dS）和aBSREL模型分析通过检查相关物种间的多样化选择过程，提供了对基因进化的全面理解。总体而言，肉苁蓉属叶绿体基因组的置换率较低，表明叶绿体基因高度保守。然而，我们在ycf2基因中检测到正选择信号（肉苁蓉 ω = 1.24；盐生肉苁蓉 ω = 1.22；管花肉苁蓉 ω = 1.17；沙苁蓉 ω = 1.10和C. phelypaea ω = 1.12），以及在C. phelypaea和盐生肉苁蓉中的chpP和rpl22基因（图4，表2）。其他PCGs主要有同义置换，rps4基因显示最高的同义置换率（ω = 0.11）（图4）。

五种研究中的肉苁蓉属植物各个保留的叶绿体基因间成对dN/dS值的箱线图

Genes	Branch	B	LRT	Test p-value	Uncorrected p-value	ω distribution over sites
rpl22	C. phelypaea	0.0035	21.7266	0.0000	0.0000	ω1 = 1.00 (99%) ω2 = 10,000 (1.4%)
clpP	盐生肉苁蓉	0.0217	9.1315	0.0254	0.0036	ω1 = 10,000,000,000 (100%)

Genes

Branch

LRT

Test p-value

Uncorrected p-value

ω distribution over sites

rpl22

C. phelypaea

0.0035

21.7266

0.0000

ω1 = 1.00 (99%)

ω2 = 10,000 (1.4%)

clpP

盐生肉苁蓉

0.0217

9.1315

0.0254

0.0036

ω1 = 10,000,000,000 (100%)

B 优化的分支长度
LRT 用于选择的似然比检验统计量
Test p-value 校正后的多重检验 P 值
Uncorrected p-value 未校正的原始 p 值
ω 在位点上的分布 推断的 ω 估计值及其各自的位点比例

肉苁蓉属植物的系统发育分析和分歧时间估算

我们使用最大似然（ML）方法，通过36种列当科植物中共享的蛋白质编码基因（PCGs）推断物种系统发育树（表S2，图S11），然后使用BEAST通过拟南芥和烟草的化石记录数据进行校准（图S12）。结果表明，这些研究中的肉苁蓉属植物显示出单系起源，与之前确定的结果一致，且分支支持值较低（BS: 66），分布在三个主要的分支中，以沙苁蓉为基底物种/分支（图S11）。第二个分支包含肉苁蓉和盐生肉苁蓉（BS: 100）（图S11），两者具有相似的叶绿体基因组结构（表S1）。第三个分支包含管花肉苁蓉和C. phelypaea（BS: 100）（图S11），它们的SSC区域具有相似的长度和组成（表S1，表S4）。在初步的ML系统发育重建之后，确定了分歧时间，以Pedicularis cheilanthifolia作为列当科的根物种，估计其分歧时间为1.55亿年前（Mya）（图S12）。肉苁蓉属的单系群体在约9570万年前分化，沙苁蓉在约7730万年前分化，其余的肉苁蓉属物种分支在约4520万年前分化。管花肉苁蓉和C. phelypaea之间的分裂发生在约2290万年前，而肉苁蓉和盐生肉苁蓉相对较年轻，其分裂估计在约1230万年前（图S12）。

肉苁蓉属植物中高度可变的间隔区的鉴定

为了识别最适合开发肉苁蓉属物种特异性分子标记的区域，使用Kimura双参数距离在同源间隔区中进行了确定（使用EMBOSS软件包中的dismat脚本；图5A）。五个叶绿体基因组中最具差异的间隔区是clpP-rps11（40.62）、rpl33-rps18（14.97）、infA-rps8（14.64）、rpl36-infA（14.38）、rpl16-rps3（12.81）和rps18-rpl20（10.18）（图5A）。

肉苁蓉属植物中超可变间隔区的鉴定，以及开发的Cis-mk01和Cis-mk02 DNA条形码标记的验证。A）比较肉苁蓉、盐生肉苁蓉、管花肉苁蓉、沙苁蓉和C. phelypaea的叶绿体基因组中IGS区域的变异性。X轴表示IGS区域；Y轴显示不同物种对间的K2p距离范围，平均K2p距离用菱形标出。B）Cis-mk01标记区域的测序色谱图，包括沙苁蓉（Cisi）、肉苁蓉（Cide）、盐生肉苁蓉（Cisa）和管花肉苁蓉（Citu），以及共识序列和比对。C）Cis-mk02标记区域的测序色谱图，包括沙苁蓉（Cisi）、肉苁蓉（Cide）、盐生肉苁蓉（Cisa）和管花肉苁蓉（Citu），以及共识序列和比对。

开发和验证用于区分肉苁蓉属植物的分子标记

我们选择了先前识别出的一个超可变间隔区（IGS），trnR-ACG-trnN-GUU，开发了两个DNA条形码标记，分别命名为Cis-mk01和Cis-mk02，设计的引物列于表S7。通过PCR扩增测试了所有四个肉苁蓉样本的总DNA（图S13）。对每个PCR产物的序列分析显示，标记Cis-mk01有九个特异性SNP位点和三个Indel位点（图5B），而标记Cis-mk02有十四个特异性SNP位点和三个Indel位点（图5C）。这些SNP和Indel位点成功区分了所有四个肉苁蓉属植物。

肉苁蓉属植物特异性DNA条形码标记设计

如在介绍中所述，叶绿体基因组中的高度可变区域可用于开发DNA条形码标记。在这项工作中，我们使用EcoPrimer软件识别这些超可变区域作为物种特异性条形码，并为每个区域设计一对引物，以准确识别每个肉苁蓉属植物（表S7）。验证了十一个位点（图6和图S14），其中Cis-mk03显示了15个特异性SNP（图6A），Cis-mk04显示了14个特异性SNP（图6B），能够成功区分本文使用的所有四个肉苁蓉属植物。而另外九个标记在图S15中显示。

Cis-mk03和Cis-mk04的验证用于区分肉苁蓉属植物。A）沙苁蓉（Cisi）、肉苁蓉（Cide）、盐生肉苁蓉（Cisa）和管花肉苁蓉（Citu）的Cis-mk03区域的测序色谱图，以及共识序列和比对。B）沙苁蓉（Cisi）、肉苁蓉（Cide）、盐生肉苁蓉（Cisa）和管花肉苁蓉（Citu）的Cis-mk04区域的测序色谱图，以及共识序列和比对。

讨论

肉苁蓉叶绿体基因组的巨大结构变异

大多数陆地植物的叶绿体基因组在基因组结构、基因顺序和实际基因含量等方面高度保守。植物叶绿体基因组的长度主要在115-165 kb之间，具有典型的四分结构，并编码105–135个基因，包括70–90个蛋白质编码基因（PCGs）、8个rRNA基因和30–40个tRNA基因，这些基因通常位于LSC或SSC区域（大约15–20个基因位于IR区域）。在我们的研究中，我们确定了四种研究中的肉苁蓉叶绿体基因组的大小较小（小于115 kb），并编码的基因数量较少。这些观察结果与肉苁蓉属植物的全寄生生活方式一致。由于它们从寄主的根或茎中获取碳、水以及其他营养物质，这些植物不需要光合作用能力，它们的叶绿体基因组容易出现基因减少而不影响植物的生存。先前的研究表明，五种肉苁蓉属植物的结构相似性相对较高。然而，我们的结果显示，研究的肉苁蓉属植物的IR区域在不同程度上收缩和扩展。其中，肉苁蓉和盐生肉苁蓉显著扩展（达到30,352 bp），而管花肉苁蓉显著收缩（降至6,593 bp）。值得注意的是，除管花肉苁蓉外，研究中的肉苁蓉属植物的IR区域中有四个rRNA基因（rrn16、rrn23、rrn4.5和rrn5）是重复的。在管花肉苁蓉中，这四个rRNA基因位于SSC区域，解释了其相对较长的SSC和较短的IR（见图S5）。此外，作为叶绿体基因组中最保守的区域，IR的急剧收缩表明管花肉苁蓉经历了更剧烈的基因丢失，反映在其叶绿体基因组大小上（管花肉苁蓉，75,375 bp），是研究中五个物种中最小的。我们推测这可能是由于肉苁蓉属植物的生长环境差异巨大，导致肉苁蓉属植物之间的叶绿体基因组结构显著不同。

肉苁蓉叶绿体基因组的广泛基因丢失和假基因化

与超过7372个完全测序的陆地植物叶绿体基因组相比，非光合植物的完全测序叶绿体基因组数量非常少（到本研究结束时约为100个）。迄今为止，仅报道了一些全寄生植物的叶绿体基因组，如一些列当科物种、菟丝子亚科和水晶兰属植物。在这项研究中，我们测序了四种肉苁蓉属植物的叶绿体基因组，这是一种列当科的全寄生植物。与其他非寄生植物的系统基因组比较将进一步揭示寄生叶绿体基因组的进化。

在四个新的肉苁蓉叶绿体基因组中，所有与光合作用和能量生产相关的基因都是假基因或丢失的，这与肉苁蓉属植物的全寄生和异养生活方式一致。所有ndh基因要么缺失要么假基因化，缺失基因的比例高于假基因。叶绿体ndh基因编码烟酰胺腺嘌呤二核苷酸-质体醌氧化还原酶复合物的亚基，认为在高光强超过光合作用能力时，减少光氧化应激，导致有害的活性氧物质产生。作为旋花科的全寄生植物，菟丝子似乎有一个完整的叶绿体基因组，但也丢失了所有ndh基因，而atp、pet、psa和psb基因要么丢失要么假基因化。其他寄生物种也有atp、pet、psa和psb基因作为丢失的基因或假基因。因此，这种在肉苁蓉属植物中的基因丢失和假基因化现象应该是其简化的叶绿体基因组的原因，与其他寄生植物一样，以适应其生活方式和环境。

合并两个物种：肉苁蓉和盐生肉苁蓉

关于肉苁蓉和盐生肉苁蓉是否应被视为两个独立的物种存在争议。之前的研究表明，盐生肉苁蓉表现出与肉苁蓉相似的化学成分和药理活性，并已在中国宁夏和甘肃省作为当地药用草药使用。我们的研究结果显示，肉苁蓉和盐生肉苁蓉具有非常相似的叶绿体基因组结构、长度和组成（肉苁蓉 109,454 bp，GC含量为36.27%；盐生肉苁蓉 111,690 bp，GC含量为36.11%；见表S1）。具有相似长度的典型四个亚区域（LSC，IRs，SSC），包括似乎正在消失的极短的SSC区域。ML系统发育树显示，肉苁蓉和盐生肉苁蓉作为一个独立的群体，与其他肉苁蓉属植物分开。此外，这两个姐妹物种的平均分歧时间估计为1230万年前。通过实地调查，我们发现盐生肉苁蓉与肉苁蓉有相同的生长环境，但在自然界中的储量更大。因此，我们建议将这两个物种视为一个物种。盐生肉苁蓉可能被用作中国药典中的肉苁蓉草药，以解决资源短缺的紧迫问题。

结论

我们对五种肉苁蓉属植物的叶绿体基因组进行了比较分析，显示由于基因丢失和假基因化导致的结构变异，主要影响光合作用基因。这些肉苁蓉属植物的基因组大小小于其他被子植物。特别是，肉苁蓉和盐生肉苁蓉经历了显著的SSC收缩和IR扩展，而管花肉苁蓉经历了SSC扩展和IR收缩。系统发育关系表明肉苁蓉属是一个约9570万年前建立的单系群体，肉苁蓉和盐生肉苁蓉被认为是两个公认的物种。物种间的序列变异使得我们开发了两个DNA条形码标记，Cis-mk01和Cis-mk02，以及在高度变异的间隔区中开发了11个DNA条形码标记（Cis-mk03至Cis-mk13），这些标记可以区分中国的肉苁蓉属植物。这里获得的结果可以提高我们对肉苁蓉属分类系统、叶绿体基因组进化及其药用产品鉴定的理解。

方法

植物材料和DNA提取

从中国的阿拉善盟（内蒙古自治区）、塔城市（新疆维吾尔自治区）、青铜峡市（宁夏回族自治区）、和田地区（新疆维吾尔自治区）分别采集了肉苁蓉、盐生肉苁蓉、沙苁蓉和管花肉苁蓉的新鲜样品。样品取自肉苁蓉种植基地，没有损害野生物种，符合地方和国家的伦理要求。样品由林余霖教授鉴定，并存放在中国医学科学院&北京协和医学院标本馆（标本登记号CMPB13484、CMPB13485、CMPB13486和CMPB13487）。新鲜样品用锡箔覆盖，液氮冷冻，并保持在-80°C直到使用。使用植物基因组DNA提取试剂盒（Tiangen Biotech，北京，中国）进行总DNA提取。通过1%（w/v）琼脂糖凝胶电泳评估总DNA质量，并按照制造商说明在Qubit 3.0（Life Technologies，Carlsbad，CA，USA）上定量。

四种肉苁蓉属植物的地理分布和生长环境。红点：内蒙古自治区阿拉善盟的肉苁蓉，也称为“肉苁蓉”。黄点：新疆维吾尔自治区塔城市的盐生肉苁蓉，也称为“盐生肉苁蓉”。蓝点：宁夏回族自治区青铜峡市的中华肉苁蓉，也称为“沙苁蓉”。绿点：新疆维吾尔自治区和田地区的管花肉苁蓉，也称为“管花肉苁蓉”。

基因组测序、组装和验证

约500 ng DNA用于构建插入片段大小为500 bp的末端配对文库，按照制造商的建议进行操作。文库在Illumina HiSeq 4000平台（Illumina Inc., San Diego, CA, USA）上进行测序，按照制造商的说明进行操作。因此，每个物种获得了约5G数据。下载所有记录在GenBank中的植物叶绿体基因组，建立本地数据库。使用BLASTn以1e-5为截止值过滤清洁的末端配对读取序列。使用SPAdes（v. 3.10.1）对过滤后的读取序列进行下游基因组组装。为了验证完整草案叶绿体基因组的正确性，我们使用Bowtie 2（v. 2.0.1）将所有原始读取序列映射到草案叶绿体基因组，并使用ggplot2构建覆盖图。构建并可视化每个叶绿体基因组的点图以评估保守的四分结构。

基因组注释、丢失基因和假基因的识别

使用CpGAVAS2网络服务对四个叶绿体基因组进行注释。BLASTn和BLASTx的E值截止设置为1E-10，以分配每个基因，并使用Apollo基因组编辑器手动编辑。使用图图云平台科研好伙伴绘制圆形基因图。使用CGView Server确定每个基因和叶绿体基因组的GC含量。组装的叶绿体基因组已存储在GenBank中，登录号为：MN614127（肉苁蓉）、MN614128（盐生肉苁蓉）、MN614129（沙苁蓉）和MN614130（管花肉苁蓉）。从NCBI下载Cistanche phelypaea的叶绿体基因组数据（NC_025642.1），并将其纳入我们的下游分析。

与已知蛋白质编码基因（PCGs）相似但相比非寄生植物列当科的地黄（Rehmannia glutinosa，NC_034308）被截断或包含一个或多个移码突变的基因定义为假基因。与列当科的非寄生植物地黄（R. glutinosa，NC_034308）相比缺失的基因定义为丢失基因。使用以下公式计算丢失基因和假基因的比例：丢失基因（假基因）比例=丢失基因（假基因）数量/地黄中的PCG数量。然后，我们使用R包“PD”基于丢失基因和假基因的比例构建拓扑树。

基因组比较

我们使用mVISTA软件在Shuffle-LAGAN模式下进行了肉苁蓉叶绿体基因组的比较基因组分析。在分析中使用注释过的地黄叶绿体基因组作为参考。使用IRscope分析五种肉苁蓉属植物在LSC/IRs和SSC/IRs边界的叶绿体基因组结构。使用E值截止为1e-10的BLASTN在五个肉苁蓉属叶绿体基因组之间识别保守序列。同源区域和基因注释使用基于网络的基因组共线性查看器GSV进行可视化。使用Mauve分析五种肉苁蓉属叶绿体基因组的基因重排。

超可变区域的识别

为了识别最具差异的区域，我们编写了一个自定义脚本，从五个叶绿体基因组的GenBank文件中提取IGS区域的起点和终点。共识别出25个由五个肉苁蓉叶绿体基因组共享的IGS。使用ClustalW2（v.2.0.12）程序对序列进行提取和比对，选项为“-type=DNA –gapopen=10 –gapext=2”。使用EMBOSS包中的distmat程序实现的K2p进化模型计算成对距离。

选择压力分析

使用Hyphy中实现的自适应分支位点随机效应可能性（aBSREL）模型进行选择压力分析。我们使用RevTrans v2.0通过其蛋白质序列引导的18个PCG的DNA序列比对，选项为CLUSTALW2。在实施程序的多重测试校正后，使用可能性比测试（p ≤ 0.05）进行显著性分析。我们使用PAML v 4.9中的yn00程序，使用F3X4密码子模型计算18个PCG的非同义置换率（dN）和同义置换率（dS）。将共享PCG推断的ML树用作输入树。

系统发育分析

对于系统发育分析，使用ClustalW程序比对包括本文肉苁蓉属植物在内的36种列当科植物的共享叶绿体PCG DNA序列。使用软件Randomized Axelerated Maximum Likelihood（RAxML）构建最大似然（ML）树，以拟南芥和烟草作为外类群。详细参数为“raxmlHPC-PTHREADS-SSE3 -fa -N 1000 –m PROTGAMMACPREV/GTRGAMMA—×551,314,260 -p 551,314,260 -o Arabidopsis_thaliana, Nicotiana_tabacum -T 20”。使用1,000次重复的自举测试评估系统发育树的显著性水平，并在每个节点之外显示自举值。

分子钟分析

我们使用软件BEAST对共享叶绿体PCG比对进行分子钟分析，使用拟南芥和烟草的化石信息，以及叶绿体DNA的cpREW氨基酸置换模型。在严格时钟方法下，使用默认设置执行MCMC分析（20,000,000代，Yule物种形成树先验的置换率，每1,000代采样一次）。使用TRACER软件检查树的接受度和收敛性，以达到平稳分布，使用TREEANNOTATOR软件在设置10%的燃烧后生成最大信度树。使用FigTree（v.1.4.3；FigTree，2019年10月）对树进行可视化。

使用EcoPrimer进行DNA条形码标记鉴定

使用EcoPrimer软件基于完整的叶绿体基因组序列鉴定DNA条形码标记。为此，我们下载了C. phelypaea的叶绿体基因组，并将其与本研究中获得的四个肉苁蓉叶绿体基因组序列结合。使用命令"ecoPCRFormat.py-g-n Cistanche.Fo-t Taxonomy Cistanche.gb"构建数据库。然后，在构建的数据库上运行命令"ecoPrimer-d Cistanche.Fo-l 100-L 1000-e 0-t species > Cistanche.Po"，为每个DNA条形码标记找到特定引物。

物种区分的分子标记和生态引物的鉴定与验证

我们使用从可变间隔区和生态引物鉴定出的标记来区分四种肉苁蓉属植物。在可变间隔区上设计的引物使用Snapgene 6.0（Insightful Science公司提供的Snapgene，网址SnapGene | Software for everyday molecular biology，最后使用于2022年），或使用EcoPrimer软件选择的条形码区域引物（表S7和S8）。PCR扩增在最终体积为20 μL的体系中进行，包括10 μL 2 × Taq PCR Master Mix，0.5 μM的每个引物，5 μL模板DNA和4 μL ddH2O，按照制造商的说明进行操作（Mei5 Biotechnology, Co., Ltd）。所有扩增在Pro-Flex PCR系统（Applied Biosystems, Waltham, MA, USA）中进行，条件为：95 ℃预变性3分钟，随后36个循环，94 ℃变性25秒，55℃退火10秒，72 ℃延伸2分钟，按照制造商的说明进行操作（Mei5 Biotechnology, Co., Ltd）。PCR扩增子在1%琼脂

糖凝胶上可视化，纯化后使用用于PCR扩增的相同引物进行双向Sanger测序，在ABI 3730 XL仪器（Applied Biosystems, USA）上使用BigDye v3.1化学试剂（Applied Biosystems）按照制造商的说明进行操作。