焦貝貝, 王希胤
( 華北理工大學(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)化模式。
收集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 共線性分析 使用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)。
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
被子植物基因組常常經(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è),
同理,
目前,很多用系統(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
長(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)建和更深層次的理解物種的演化歷程。