基于Ks分布的被子植物演化的時(shí)間尺度研究

焦貝貝， 王希胤

( 華北理工大學(xué) 生命科學(xué)學(xué)院， 河北 唐山 063210 )

被子植物的起源和早期快速演化及其發(fā)生時(shí)間一直是生物學(xué)的研究熱點(diǎn)。當(dāng)前估算物種演化時(shí)間的方法主要是基于分子鐘假設(shè)，即以某幾個(gè)特定類群的化石時(shí)間作為校正點(diǎn)，然后通過(guò)部分基因序列間的相似性，假設(shè)不同的物種擁有相同或相近的進(jìn)化速率，來(lái)估算系統(tǒng)發(fā)育樹上某一節(jié)點(diǎn)的時(shí)間，從而推斷出該類群的起源時(shí)間(唐先華等，2002；Donoghue & Yang，2016；Luo et al.,2020)。然而，近年的研究表明，不同物種的分子鐘通常具有顯著差異，即不同物種的進(jìn)化速率有顯著不同(Wang et al.，2017；2019)，不同年代具有不同的進(jìn)化速度(羅靜和張亞平，2000；Smith & Donoghue, 2008)，且在不同研究中，對(duì)分子進(jìn)化速率的估算有很大的差異(Lanfear et al., 2010)。此外，引入的化石時(shí)間對(duì)估算的時(shí)間影響很大，隨著更多化石且更準(zhǔn)確的年份測(cè)定，被子植物演化的時(shí)間尺度會(huì)隨之變動(dòng)(Hug & Roger，2007；Wang et al.，2015；Silvestro et al.，2021)。

基因組測(cè)序揭示了歷史上反復(fù)的多倍化事件(Ren et al., 2018)，多倍化事件使基因組內(nèi)所有基因發(fā)生重復(fù)，且基因組中的古老同源區(qū)域通常有相當(dāng)數(shù)目的重復(fù)基因保留下來(lái)，從而形成目前基因組內(nèi)或者基因組間的共線性同源基因(Jiao et al., 2011)。對(duì)共線性同源基因的分析，是揭示古代的多倍化或物種分歧事件并推定其發(fā)生時(shí)間和規(guī)模的重要途徑。多倍化發(fā)生后植物基因組通常會(huì)變得很不穩(wěn)定，進(jìn)化速率也變得顯著不同。由于減少了選擇性約束，因此這些重復(fù)基因通常以更快的速度進(jìn)化(Wang et al., 2016)。例如，在葫蘆科植物基因組的研究中發(fā)現(xiàn)，甜瓜的進(jìn)化速度最慢，西瓜和黃瓜的進(jìn)化速度分別快23.6%和27.4%(Wang et al., 2018)。

一般認(rèn)為，同義突變率(synonymous substitution rate，Ks)往往不會(huì)改變氨基酸的組成，不受自然選擇的影響。因此，Ks分布常常作為判定物種歷史上發(fā)生的多倍化或物種分歧事件的依據(jù)(Vanneste et al., 2013)。依據(jù)共享的演化事件應(yīng)該有相同的Ks峰值，Wang等(2015)首次提出了基于Ks峰值的矯正方法用以估算物種演化的時(shí)間尺度，得到了其他研究工作者的認(rèn)可，還被廣泛應(yīng)用于他們的研究中(Zhuang et al.，2019；Song et al.，2020；Song et al.，2021；Wang et al.，2021)。例如，兩個(gè)團(tuán)隊(duì)分別對(duì)睡蓮(Zhang et al., 2020a)和芡實(shí)(Yang et al., 2020)基因組分析，Yang等(2020)通過(guò)Ks峰值矯正的方式估算的芡實(shí)古老多倍化(被證實(shí)為睡蓮目共享)與另外的團(tuán)隊(duì)基于睡蓮目的轉(zhuǎn)錄組數(shù)據(jù)估算的時(shí)間尺度基本一致?；贙s峰值的矯正方法中，獲得準(zhǔn)確的Ks峰是準(zhǔn)確估算時(shí)間尺度的關(guān)鍵。然而，當(dāng)前獲取Ks分布的方式不統(tǒng)一且通常帶有長(zhǎng)尾現(xiàn)象(Tang et al., 2008)。為何Ks分布會(huì)有長(zhǎng)尾現(xiàn)象？長(zhǎng)尾現(xiàn)象對(duì)Ks峰是否有重要影響等問(wèn)題，也尚未有清晰的表述。

目前，已有400余種被子植物的基因組得到不同水平的測(cè)定，便于在全基因組的尺度上理解這些被子植物的演化歷程(Kress et al., 2022)。全基因組數(shù)據(jù)能有效消除橫向基因轉(zhuǎn)移和類群間基因進(jìn)化速率差異等因素對(duì)系統(tǒng)發(fā)育樹的影響。因此，急需在全基因組數(shù)據(jù)層面上，利用新方法對(duì)被子植物的演化時(shí)間進(jìn)行重新估計(jì)。本文擬對(duì)三種獲取Ks分布的方式進(jìn)行比較，明確哪種方式獲得的Ks峰值更接近真實(shí)情況；對(duì)于Ks分布中常見(jiàn)的長(zhǎng)尾現(xiàn)象，采用模擬仿真的方式，探究出現(xiàn)長(zhǎng)尾現(xiàn)象的原因；區(qū)分共享多倍化和共享早期分化兩種情況，創(chuàng)建基于全基因組數(shù)據(jù)的Ks分布矯正模型，對(duì)44個(gè)代表性被子植物基因組演化事件的時(shí)間尺度進(jìn)行重新估計(jì)，得到相對(duì)可靠的被子植物演化時(shí)間軸。這有助于更深層地了解被子植物多樣性和系統(tǒng)發(fā)育以及被子植物基因組的進(jìn)化模式。

1 材料與方法

1.1 基因組數(shù)據(jù)材料

收集44個(gè)高質(zhì)量染色體水平的被子植物基因組(主要來(lái)自NCBI和PHYTOZOME)，共包含43科39目(表1)。

表 1 研究所用的44個(gè)被子植物及基因組數(shù)據(jù)來(lái)源Table 1 List of the 44 angiosperms involved and the genome data sources

續(xù)表1

1.2 方法

1.2.1 共線性分析 使用WGDI v0.5.3(Sun et al., 2021)軟件進(jìn)行共線性分析。首先，使用BLASTP來(lái)識(shí)別基因組內(nèi)或基因組間的基因相似性。隨后，用WGDI軟件的‘-d’子程序繪制同源點(diǎn)陣圖，并運(yùn)行‘-icl’子程序獲得共線性基因。

1.2.2 Ks分布 Ks分布主要是通過(guò)WGDI軟件完成的。首先，使用WGDI軟件的‘-ks’子程序調(diào)用PAML(Yang, 2007)軟件計(jì)算共線性基因?qū)Φ腒s值。通過(guò)‘-bi’子程序整合共線性和Ks值的結(jié)果，并使用WGDI軟件的‘-bk’子程序查看共線基因的Ks值的分布，結(jié)果以點(diǎn)圖的形式展示(圖1：A)。根據(jù)物種內(nèi)或種間已知的多倍化或分歧事件，通過(guò)WGDI的‘-c’子程序?qū)簿€性片段進(jìn)行過(guò)濾，只保留多倍化事件或分歧事件產(chǎn)生的共線性片段。然后，通過(guò)WGDI的‘-kp’子程序獲取Ks分布(圖1：B)。最后，使用WGDI中的“-pf”子程序?qū)Σ煌录謩e進(jìn)行擬合并獲取Ks分布(圖1：C)。

2 結(jié)果與分析

2.1 Ks分布和長(zhǎng)尾現(xiàn)象解析

Ks分布常常用來(lái)判定物種歷史上發(fā)生的多倍化或物種分歧事件的依據(jù)。目前獲取Ks分布主要有三種方式。方式一：先通過(guò)OrthoMCL(Li et al., 2003)等聚類軟件獲取旁系同源基因?qū)?，再?jì)算這些同源基因?qū)Φ腒s值并繪制Ks分布圖。方式二：先進(jìn)行基因組共線性分析，再計(jì)算共線性基因?qū)s值并繪制Ks分布。方式三：在方式二的基礎(chǔ)上，提取共線性區(qū)塊上Ks值的中位數(shù)并繪制Ks分布。三種方式中，方式一由于沒(méi)有共線性分析，因此所獲取的旁系同源基因?qū)νǔ?huì)有大量串聯(lián)重復(fù)基因從而影響Ks分布。方式二和三都經(jīng)過(guò)了共線性分析，當(dāng)把共線性區(qū)塊(長(zhǎng)度大于5)上同源基因?qū)Φ腒s值以點(diǎn)圖的形式展示出來(lái)時(shí)(圖1：A)，這里以水稻為例，可以看到大部分由綠色的點(diǎn)組成的片段，如8號(hào)與9號(hào)染色體，這與水稻近期的一次多倍化事件相符。Ks點(diǎn)圖中大部分點(diǎn)的顏色相近，說(shuō)明Ks值波動(dòng)很小。對(duì)共線性區(qū)域的Ks值的中位數(shù)(方式三)、平均值和所有的基因?qū)?方式二)進(jìn)行正態(tài)分布擬合(帶寬為0.01，homo范圍0.3～1)(圖2：B)，可以看到方式二并沒(méi)有產(chǎn)生明顯的峰，而且Ks分布整體帶有長(zhǎng)長(zhǎng)的尾巴。方式三和區(qū)塊的平均值的Ks分布有明顯峰值，數(shù)據(jù)更為集中。由于中位數(shù)是對(duì)總體中心很好的估計(jì)，且穩(wěn)健性更強(qiáng)，中位數(shù)的峰值顏色和Ks點(diǎn)圖的顏色更為接近，因此區(qū)塊的Ks值的中位數(shù)更接近Ks真正的峰值，對(duì)方式三的Ks分布按照正態(tài)分布擬合來(lái)提取Ks峰值(圖1：C)。

A. 水稻基因組的共線性區(qū)塊; B. 共線性區(qū)塊上Ks值的擬合分布; C. 共線性區(qū)塊Ks值的核密度估計(jì)。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

當(dāng)v為常數(shù)值時(shí)，X′=X×vn；

當(dāng)v服從正態(tài)分布時(shí)，X′=X×Xvn。

當(dāng)假設(shè)Ks值的時(shí)間累積系數(shù)v為一個(gè)常數(shù)值時(shí)，設(shè)置假設(shè)的Ks分布為X～N(μ,σ2)，依據(jù)Ks分布數(shù)據(jù)特征，設(shè)定μ=0.2，σ=0.01，v=1.02，n=100。每迭代10次，繪制Ks分布結(jié)果(圖2：A)。隨著進(jìn)化事件的推移，Ks峰值也逐漸變大， Ks分布依舊完美符合正態(tài)分布且沒(méi)有明顯的長(zhǎng)尾分布現(xiàn)象。

A. Ks分布在恒定進(jìn)化速率下的模擬; B. Ks分布在進(jìn)化速率服從正態(tài)分布的模擬。A. Simulation of Ks distribution at a constant evolution rate; B. Simulation of Ks distribution under a normal distribution of evolution rates.圖 2 Ks分布在不同進(jìn)化速率下的模擬結(jié)果Fig. 2 Simulation results of Ks distribution at different evolution rates

當(dāng)假設(shè)Ks值的時(shí)間累積系數(shù)(v)服從正態(tài)分布時(shí)，最初設(shè)置假設(shè)的Ks分布為Xv～N(μv,σv2)，其中μ=0.2，σ=0.01，μv=1.02，σv=0.01，n=100。每迭代10次，繪制Ks分布結(jié)果(圖2：B)。隨著進(jìn)化事件的推移，Ks峰值逐漸變大，Ks分布不再是正態(tài)分布，并帶有明顯的長(zhǎng)尾現(xiàn)象。由于這種假設(shè)所得到的結(jié)果更接近于真實(shí)情況，因此基因的進(jìn)化速率不是相對(duì)恒定的，它可能并非等速進(jìn)行，而是在不同年代具有不同的進(jìn)化速度，這可能符合正態(tài)分布。對(duì)模擬的Ks分布通過(guò)高斯擬合獲取峰值時(shí)，發(fā)現(xiàn)Ks峰值與進(jìn)化速率勻速時(shí)的沒(méi)有明顯差異(表2)。因此，Ks分布中長(zhǎng)尾現(xiàn)象對(duì)提取到的Ks峰值的影響較小。

表 2 不同進(jìn)化速率模擬下的Ks峰值Table 2 Ks peaks under simulations at different evolution rates

2.2 Ks分布矯正方法

被子植物基因組常常經(jīng)歷不止一次多倍化事件，不同物種的進(jìn)化速率顯著不同，從而導(dǎo)致共享的多倍化事件的Ks峰值也大不相同。而Ks分布矯正方法的核心理念就是將這些共享事件的Ks峰矯正到一起。根據(jù)共享事件的不同，Ks分布矯正方法可分為共享多倍化和共享分化兩種情況。

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

如果兩個(gè)物種A、B雖不存在共享的多倍化事件但存在共享的早期分化事件，就通過(guò)尋找外類群來(lái)輔助矯正(圖3：B)。物種C、D、E是外類群，物種C和D的祖先在P點(diǎn)與物種A、B的祖先分化，所以物種C與A、B間的Ks峰值應(yīng)該相等，物種D與A、B間的Ks峰值也應(yīng)該相等，即KsCA=KsCB，KsDA=KsDB。同樣，由于物種間不同的進(jìn)化速率，因此現(xiàn)實(shí)情況下它們大多不相等。按照前面的假設(shè)，

同理，

2.3 被子植物系統(tǒng)發(fā)育樹時(shí)間矯正

目前，很多用系統(tǒng)發(fā)育樹的方法推測(cè)被子植物的演化時(shí)間，認(rèn)為被子植物的起源為三疊紀(jì)225百萬(wàn)年至240百萬(wàn)年前(Magallón, 2010)，這與起傳粉作用的核心植食性鱗翅目昆蟲的起源時(shí)間(約230百萬(wàn)年前)一致(Li et al., 2019)。由于無(wú)油樟目和睡蓮目、核心被子植物五大分支之間的關(guān)系仍然沒(méi)有完全解析，且已有多個(gè)證據(jù)暗示核心被子植物祖先可能發(fā)生了快速輻射分化(Yang et al., 2020)。因此，在矯正過(guò)程中，以無(wú)油樟目為作為參考，不討論它和睡蓮目的關(guān)系，認(rèn)為五大分支的分化時(shí)間尺度在同一個(gè)時(shí)間范圍內(nèi)?；诤诵恼骐p子葉植物共享的γ事件，時(shí)間范圍為115～130百萬(wàn)年(Million years ago，Mya)，對(duì)44個(gè)被子植物基因組(表1)進(jìn)行了時(shí)間尺度矯正(圖4)。從矯正后的時(shí)間尺度來(lái)看，被子植物在130百萬(wàn)年前附近，單子葉植物、真雙子葉植物、木蘭類植物祖先都發(fā)生了快速輻射進(jìn)化，與Zhang等(2020b)的結(jié)論一致。此外，在早白堊世(130百萬(wàn)年)時(shí)期，白堊紀(jì)-古新世(K-Pg)邊界時(shí)期(66百萬(wàn)年)和中新世(20百萬(wàn)年，靠近冰川期)很多被子植物發(fā)生的多倍化事件，研究發(fā)現(xiàn)WGD的時(shí)間在被子植物的系統(tǒng)發(fā)育中并不是隨機(jī)分布與Wu等(2020)的結(jié)論一致。

圖 4 時(shí)間矯正后的被子植物系統(tǒng)發(fā)育樹Fig. 4 Angiosperm phylogenetic tree after time correction

盡管不同物種的進(jìn)化速率數(shù)值顯著不同，但是同一類群中的進(jìn)化速率往往具有部分一致性。由矯正方法可知，矯正后的Ks峰值應(yīng)該相等。因此，Ks峰值越大，表明進(jìn)化速率越快。對(duì)木蘭類植物、真雙子葉植物和單子葉植物與無(wú)油樟的Ks峰值的比較發(fā)現(xiàn)，木蘭類植物(大多數(shù)為木本)進(jìn)化速率最慢，真雙子葉植物(大多數(shù)為灌木)次之，單子葉植物(大多數(shù)為草本)進(jìn)化速率最快(表3)，這與多年生木本植物比草本植物的分子進(jìn)化速率慢的結(jié)論相符(Lanfear et al., 2013)。此外，對(duì)多倍化事件發(fā)生的時(shí)間與矯正前后的Ks峰值比較(圖5)發(fā)現(xiàn)，矯正前的Ks峰值與時(shí)間并不是線性關(guān)系。隨著Ks峰值的增大，多倍化事件發(fā)生的時(shí)間并沒(méi)有更古老； 由于矯正后的Ks峰值與時(shí)間成正比，因此對(duì)Ks峰值進(jìn)行矯正之后估算物種演化事件的時(shí)間是十分必要的。

表 3 部分核心被子植物與無(wú)油樟之間的Ks峰值Table 3 Ks peaks between some species of mesangiospermae and Amborella trichopoda

圖 5 矯正前后Ks峰值與時(shí)間的關(guān)系Fig. 5 Relationship between Ks peaks and time before and after correction

3 討論與結(jié)論

長(zhǎng)期以來(lái)，估算被子植物演化的時(shí)間尺度主要是基于分子鐘假設(shè)，然而分子進(jìn)化異速現(xiàn)象的廣泛存在嚴(yán)重影響其準(zhǔn)確性，Wang等(2015)提出的基于Ks分布的矯正方法，獲得了令人信服的時(shí)間尺度。本文對(duì)獲取Ks分布三種常見(jiàn)的方式進(jìn)行了比較分析， 明確了通過(guò)提取共線性區(qū)塊上Ks值的中位數(shù)更能代表真實(shí)的Ks峰值。此外，還進(jìn)一步解析了Ks分布中常見(jiàn)的長(zhǎng)尾現(xiàn)象，本研究模擬結(jié)果表明基因的進(jìn)化速率并非相對(duì)恒定和等速進(jìn)行。當(dāng)假設(shè)進(jìn)化速率并非相對(duì)恒定，而是符合正態(tài)分布的時(shí)候，Ks分布出現(xiàn)了有明顯的長(zhǎng)尾現(xiàn)象，但這并不影響提取到的Ks峰值的準(zhǔn)確性。Vanneste等(2013)研究表明，當(dāng)Ks值大于1時(shí)，容易受到飽和效應(yīng)的影響，并且隨著Ks值增大，這種效應(yīng)越明顯。模擬的Ks峰值范圍接近于1，隨著Ks峰值增大，估計(jì)的Ks峰值可能會(huì)受到飽和效應(yīng)的影響。

本研究還詳細(xì)描述了基于Ks峰值的矯正方法的矯正過(guò)程。先前的研究只對(duì)共享多倍化和共享早期分化兩種情況分開進(jìn)行了描述，這是首次全面的描述，有助于深入理解和傳播?；谠摲椒ǎ€對(duì)44個(gè)高質(zhì)量的被子植物基因組演化事件的時(shí)間尺度進(jìn)行了重新估計(jì)，估計(jì)結(jié)果與近期發(fā)表的時(shí)間尺度基本一致(Li et al., 2019；Wu et al., 2020)。本研究結(jié)果還表明，被子植物基因組的進(jìn)化速率雖然差異顯著，但不同分支間的進(jìn)化速率仍具有一致性。并且，不同譜系的被子植物具有同步的輻射進(jìn)化和適應(yīng)性進(jìn)化現(xiàn)象。隨著更多高質(zhì)量的被子植物基因組的公布和有效化石年份的準(zhǔn)確測(cè)定，被子植物演化的時(shí)間尺度會(huì)越來(lái)越清晰，更有利于植物系統(tǒng)發(fā)育的構(gòu)建和更深層次的理解物種的演化歷程。