基于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峰值矯正的方式估算的芡實古老多倍化(被證實為睡蓮目共享)與另外的團隊基于睡蓮目的轉錄組數據估算的時間尺度基本一致?；贙s峰值的矯正方法中，獲得準確的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)。因此，在矯正過程中，以無油樟目為作為參考，不討論它和睡蓮目的關系，認為五大分支的分化時間尺度在同一個時間范圍內?；诤诵恼骐p子葉植物共享的γ事件，時間范圍為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)。本研究結果還表明，被子植物基因組的進化速率雖然差異顯著，但不同分支間的進化速率仍具有一致性。并且，不同譜系的被子植物具有同步的輻射進化和適應性進化現象。隨著更多高質量的被子植物基因組的公布和有效化石年份的準確測定，被子植物演化的時間尺度會越來越清晰，更有利于植物系統發育的構建和更深層次的理解物種的演化歷程。