基于Ks分布的被子植物演化的时间尺度研究

焦贝贝， 王希胤

( 华北理工大学 生命科学学院， 河北 唐山 063210 )

被子植物的起源和早期快速演化及其发生时间一直是生物学的研究热点。当前估算物种演化时间的方法主要是基于分子钟假设，即以某几个特定类群的化石时间作为校正点，然后通过部分基因序列间的相似性，假设不同的物种拥有相同或相近的进化速率，来估算系统发育树上某一节点的时间，从而推断出该类群的起源时间(唐先华等，2002；Donoghue & Yang，2016；Luo et al.,2020)。然而，近年的研究表明，不同物种的分子钟通常具有显著差异，即不同物种的进化速率有显著不同(Wang et al.，2017；2019)，不同年代具有不同的进化速度(罗静和张亚平，2000；Smith & Donoghue, 2008)，且在不同研究中，对分子进化速率的估算有很大的差异(Lanfear et al., 2010)。此外，引入的化石时间对估算的时间影响很大，随着更多化石且更准确的年份测定，被子植物演化的时间尺度会随之变动(Hug & Roger，2007；Wang et al.，2015；Silvestro et al.，2021)。

基因组测序揭示了历史上反复的多倍化事件(Ren et al., 2018)，多倍化事件使基因组内所有基因发生重复，且基因组中的古老同源区域通常有相当数目的重复基因保留下来，从而形成目前基因组内或者基因组间的共线性同源基因(Jiao et al., 2011)。对共线性同源基因的分析，是揭示古代的多倍化或物种分歧事件并推定其发生时间和规模的重要途径。多倍化发生后植物基因组通常会变得很不稳定，进化速率也变得显著不同。由于减少了选择性约束，因此这些重复基因通常以更快的速度进化(Wang et al., 2016)。例如，在葫芦科植物基因组的研究中发现，甜瓜的进化速度最慢，西瓜和黄瓜的进化速度分别快23.6%和27.4%(Wang et al., 2018)。

一般认为，同义突变率(synonymous substitution rate，Ks)往往不会改变氨基酸的组成，不受自然选择的影响。因此，Ks分布常常作为判定物种历史上发生的多倍化或物种分歧事件的依据(Vanneste et al., 2013)。依据共享的演化事件应该有相同的Ks峰值，Wang等(2015)首次提出了基于Ks峰值的矫正方法用以估算物种演化的时间尺度，得到了其他研究工作者的认可，还被广泛应用于他们的研究中(Zhuang et al.，2019；Song et al.，2020；Song et al.，2021；Wang et al.，2021)。例如，两个团队分别对睡莲(Zhang et al., 2020a)和芡实(Yang et al., 2020)基因组分析，Yang等(2020)通过Ks峰值矫正的方式估算的芡实古老多倍化(被证实为睡莲目共享)与另外的团队基于睡莲目的转录组数据估算的时间尺度基本一致。基于Ks峰值的矫正方法中，获得准确的Ks峰是准确估算时间尺度的关键。然而，当前获取Ks分布的方式不统一且通常带有长尾现象(Tang et al., 2008)。为何Ks分布会有长尾现象？长尾现象对Ks峰是否有重要影响等问题，也尚未有清晰的表述。

目前，已有400余种被子植物的基因组得到不同水平的测定，便于在全基因组的尺度上理解这些被子植物的演化历程(Kress et al., 2022)。全基因组数据能有效消除横向基因转移和类群间基因进化速率差异等因素对系统发育树的影响。因此，急需在全基因组数据层面上，利用新方法对被子植物的演化时间进行重新估计。本文拟对三种获取Ks分布的方式进行比较，明确哪种方式获得的Ks峰值更接近真实情况；对于Ks分布中常见的长尾现象，采用模拟仿真的方式，探究出现长尾现象的原因；区分共享多倍化和共享早期分化两种情况，创建基于全基因组数据的Ks分布矫正模型，对44个代表性被子植物基因组演化事件的时间尺度进行重新估计，得到相对可靠的被子植物演化时间轴。这有助于更深层地了解被子植物多样性和系统发育以及被子植物基因组的进化模式。

1 材料与方法

1.1 基因组数据材料

收集44个高质量染色体水平的被子植物基因组(主要来自NCBI和PHYTOZOME)，共包含43科39目(表1)。

表 1 研究所用的44个被子植物及基因组数据来源Table 1 List of the 44 angiosperms involved and the genome data sources

续表1

1.2 方法

1.2.1 共线性分析 使用WGDI v0.5.3(Sun et al., 2021)软件进行共线性分析。首先，使用BLASTP来识别基因组内或基因组间的基因相似性。随后，用WGDI软件的‘-d’子程序绘制同源点阵图，并运行‘-icl’子程序获得共线性基因。

1.2.2 Ks分布 Ks分布主要是通过WGDI软件完成的。首先，使用WGDI软件的‘-ks’子程序调用PAML(Yang, 2007)软件计算共线性基因对的Ks值。通过‘-bi’子程序整合共线性和Ks值的结果，并使用WGDI软件的‘-bk’子程序查看共线基因的Ks值的分布，结果以点图的形式展示(图1：A)。根据物种内或种间已知的多倍化或分歧事件，通过WGDI的‘-c’子程序对共线性片段进行过滤，只保留多倍化事件或分歧事件产生的共线性片段。然后，通过WGDI的‘-kp’子程序获取Ks分布(图1：B)。最后，使用WGDI中的“-pf”子程序对不同事件分别进行拟合并获取Ks分布(图1：C)。

2 结果与分析

2.1 Ks分布和长尾现象解析

Ks分布常常用来判定物种历史上发生的多倍化或物种分歧事件的依据。目前获取Ks分布主要有三种方式。方式一：先通过OrthoMCL(Li et al., 2003)等聚类软件获取旁系同源基因对，再计算这些同源基因对的Ks值并绘制Ks分布图。方式二：先进行基因组共线性分析，再计算共线性基因对Ks值并绘制Ks分布。方式三：在方式二的基础上，提取共线性区块上Ks值的中位数并绘制Ks分布。三种方式中，方式一由于没有共线性分析，因此所获取的旁系同源基因对通常会有大量串联重复基因从而影响Ks分布。方式二和三都经过了共线性分析，当把共线性区块(长度大于5)上同源基因对的Ks值以点图的形式展示出来时(图1：A)，这里以水稻为例，可以看到大部分由绿色的点组成的片段，如8号与9号染色体，这与水稻近期的一次多倍化事件相符。Ks点图中大部分点的颜色相近，说明Ks值波动很小。对共线性区域的Ks值的中位数(方式三)、平均值和所有的基因对(方式二)进行正态分布拟合(带宽为0.01，homo范围0.3～1)(图2：B)，可以看到方式二并没有产生明显的峰，而且Ks分布整体带有长长的尾巴。方式三和区块的平均值的Ks分布有明显峰值，数据更为集中。由于中位数是对总体中心很好的估计，且稳健性更强，中位数的峰值颜色和Ks点图的颜色更为接近，因此区块的Ks值的中位数更接近Ks真正的峰值，对方式三的Ks分布按照正态分布拟合来提取Ks峰值(图1：C)。

A. 水稻基因组的共线性区块; B. 共线性区块上Ks值的拟合分布; C. 共线性区块Ks值的核密度估计。A. Synteny blocks of the Oryza sativa genome; B. Fitted distribution of Ks values for synteny blocks; C. Kernel density of Ks values for synteny blocks.图 1 Ks分布Fig. 1 Ks distribution

当v为常数值时，X′=X×vn；

当v服从正态分布时，X′=X×Xvn。

当假设Ks值的时间累积系数v为一个常数值时，设置假设的Ks分布为X～N(μ,σ2)，依据Ks分布数据特征，设定μ=0.2，σ=0.01，v=1.02，n=100。每迭代10次，绘制Ks分布结果(图2：A)。随着进化事件的推移，Ks峰值也逐渐变大， Ks分布依旧完美符合正态分布且没有明显的长尾分布现象。

A. Ks分布在恒定进化速率下的模拟; B. Ks分布在进化速率服从正态分布的模拟。A. Simulation of Ks distribution at a constant evolution rate; B. Simulation of Ks distribution under a normal distribution of evolution rates.图 2 Ks分布在不同进化速率下的模拟结果Fig. 2 Simulation results of Ks distribution at different evolution rates

当假设Ks值的时间累积系数(v)服从正态分布时，最初设置假设的Ks分布为Xv～N(μv,σv2)，其中μ=0.2，σ=0.01，μv=1.02，σv=0.01，n=100。每迭代10次，绘制Ks分布结果(图2：B)。随着进化事件的推移，Ks峰值逐渐变大，Ks分布不再是正态分布，并带有明显的长尾现象。由于这种假设所得到的结果更接近于真实情况，因此基因的进化速率不是相对恒定的，它可能并非等速进行，而是在不同年代具有不同的进化速度，这可能符合正态分布。对模拟的Ks分布通过高斯拟合获取峰值时，发现Ks峰值与进化速率匀速时的没有明显差异(表2)。因此，Ks分布中长尾现象对提取到的Ks峰值的影响较小。

表 2 不同进化速率模拟下的Ks峰值Table 2 Ks peaks under simulations at different evolution rates

2.2 Ks分布矫正方法

被子植物基因组常常经历不止一次多倍化事件，不同物种的进化速率显著不同，从而导致共享的多倍化事件的Ks峰值也大不相同。而Ks分布矫正方法的核心理念就是将这些共享事件的Ks峰矫正到一起。根据共享事件的不同，Ks分布矫正方法可分为共享多倍化和共享分化两种情况。

A. 共享多倍化事件; B. 共享早期分化。 A. Shared polyploidy events; B. Shared early divergence.图 3 Ks分布矫正方法的原理Fig. 3 Principle of the Ks distribution correction method

如果两个物种A、B虽不存在共享的多倍化事件但存在共享的早期分化事件，就通过寻找外类群来辅助矫正(图3：B)。物种C、D、E是外类群，物种C和D的祖先在P点与物种A、B的祖先分化，所以物种C与A、B间的Ks峰值应该相等，物种D与A、B间的Ks峰值也应该相等，即KsCA=KsCB，KsDA=KsDB。同样，由于物种间不同的进化速率，因此现实情况下它们大多不相等。按照前面的假设，

同理，

2.3 被子植物系统发育树时间矫正

目前，很多用系统发育树的方法推测被子植物的演化时间，认为被子植物的起源为三叠纪225百万年至240百万年前(Magallón, 2010)，这与起传粉作用的核心植食性鳞翅目昆虫的起源时间(约230百万年前)一致(Li et al., 2019)。由于无油樟目和睡莲目、核心被子植物五大分支之间的关系仍然没有完全解析，且已有多个证据暗示核心被子植物祖先可能发生了快速辐射分化(Yang et al., 2020)。因此，在矫正过程中，以无油樟目为作为参考，不讨论它和睡莲目的关系，认为五大分支的分化时间尺度在同一个时间范围内。基于核心真双子叶植物共享的γ事件，时间范围为115～130百万年(Million years ago，Mya)，对44个被子植物基因组(表1)进行了时间尺度矫正(图4)。从矫正后的时间尺度来看，被子植物在130百万年前附近，单子叶植物、真双子叶植物、木兰类植物祖先都发生了快速辐射进化，与Zhang等(2020b)的结论一致。此外，在早白垩世(130百万年)时期，白垩纪-古新世(K-Pg)边界时期(66百万年)和中新世(20百万年，靠近冰川期)很多被子植物发生的多倍化事件，研究发现WGD的时间在被子植物的系统发育中并不是随机分布与Wu等(2020)的结论一致。

图 4 时间矫正后的被子植物系统发育树Fig. 4 Angiosperm phylogenetic tree after time correction

尽管不同物种的进化速率数值显著不同，但是同一类群中的进化速率往往具有部分一致性。由矫正方法可知，矫正后的Ks峰值应该相等。因此，Ks峰值越大，表明进化速率越快。对木兰类植物、真双子叶植物和单子叶植物与无油樟的Ks峰值的比较发现，木兰类植物(大多数为木本)进化速率最慢，真双子叶植物(大多数为灌木)次之，单子叶植物(大多数为草本)进化速率最快(表3)，这与多年生木本植物比草本植物的分子进化速率慢的结论相符(Lanfear et al., 2013)。此外，对多倍化事件发生的时间与矫正前后的Ks峰值比较(图5)发现，矫正前的Ks峰值与时间并不是线性关系。随着Ks峰值的增大，多倍化事件发生的时间并没有更古老； 由于矫正后的Ks峰值与时间成正比，因此对Ks峰值进行矫正之后估算物种演化事件的时间是十分必要的。

表 3 部分核心被子植物与无油樟之间的Ks峰值Table 3 Ks peaks between some species of mesangiospermae and Amborella trichopoda

图 5 矫正前后Ks峰值与时间的关系Fig. 5 Relationship between Ks peaks and time before and after correction

3 讨论与结论

长期以来，估算被子植物演化的时间尺度主要是基于分子钟假设，然而分子进化异速现象的广泛存在严重影响其准确性，Wang等(2015)提出的基于Ks分布的矫正方法，获得了令人信服的时间尺度。本文对获取Ks分布三种常见的方式进行了比较分析， 明确了通过提取共线性区块上Ks值的中位数更能代表真实的Ks峰值。此外，还进一步解析了Ks分布中常见的长尾现象，本研究模拟结果表明基因的进化速率并非相对恒定和等速进行。当假设进化速率并非相对恒定，而是符合正态分布的时候，Ks分布出现了有明显的长尾现象，但这并不影响提取到的Ks峰值的准确性。Vanneste等(2013)研究表明，当Ks值大于1时，容易受到饱和效应的影响，并且随着Ks值增大，这种效应越明显。模拟的Ks峰值范围接近于1，随着Ks峰值增大，估计的Ks峰值可能会受到饱和效应的影响。

本研究还详细描述了基于Ks峰值的矫正方法的矫正过程。先前的研究只对共享多倍化和共享早期分化两种情况分开进行了描述，这是首次全面的描述，有助于深入理解和传播。基于该方法，还对44个高质量的被子植物基因组演化事件的时间尺度进行了重新估计，估计结果与近期发表的时间尺度基本一致(Li et al., 2019；Wu et al., 2020)。本研究结果还表明，被子植物基因组的进化速率虽然差异显著，但不同分支间的进化速率仍具有一致性。并且，不同谱系的被子植物具有同步的辐射进化和适应性进化现象。随着更多高质量的被子植物基因组的公布和有效化石年份的准确测定，被子植物演化的时间尺度会越来越清晰，更有利于植物系统发育的构建和更深层次的理解物种的演化历程。