皺邊喉毛花及其近緣種葉綠體基因組結構與進化分析

韓 霜,余靜雅,李孝平,牛 玉,張發起

(1 中國科學院西北高原生物研究所高原生物適應與進化重點實驗室,西寧 810001;2 中國科學院大學 生命科學學院,北京 100039)

基因組學與系統發育生物學相結合,建立新的學科——系統發育基因組學(Phylogenomics),也為正確解決多基因序列系統發育樹沖突等問題提供新的方法和見解[1,2-4]。基因組包含豐富的基因序列信息及基因復制、轉移、丟失等事件的重要信息,為系統發育研究提供大量證據。系統發育重建進入基因組時代,越來越多的學者將核基因組與質體基因組應用到系統發育學研究中[5-6],通過基因組數據構建更高分辨率、更成熟的系統發育樹,進而更準確地推斷植物類群的進化關系。對于高等植物和一些藻類植物而言,其體內固有的細胞器葉綠體在植物的光合作用過程中具有重要作用[7]。此外,由于葉綠體為母系遺傳,高度保守、變異位點多并具有獨立的基因組,是研究植物系統發育的理想目標[5]。測序技術的迅速發展與測序成本的降低,實現了許多植物的葉綠體基因組測序,為解決植物類群系統發育問題提供新的思路[8-9]。近年來,許多學者利用大量有價值的基因組數據解決了長期存在爭議的植物系統發育問題。在這些研究案例中,對薔薇科的123個葉綠體基因組進行測序,應用多種系統發育重建方法明晰薔薇科的深層系統發育關系,為其多樣化歷史提供了新的重要見解[10]。對麻黃屬3個物種的葉綠體基因組進行測序,為其系統發育研究提供有價值的信息[11]。

喉毛花屬(Comastoma)隸屬于龍膽科(Gentianaceae)獐牙菜亞族(subtribe Swertiinae),多為一年生或多年生草本植物,全球共15種,中國有11種[12]。該屬植物主要分布在青藏高原及其毗鄰地區,也是藏茵陳入藥基原植物之一,在保肝、抗氧化、抗病毒等方面具有顯著效果[13-15]。綜合以往研究分析,學者基于ITS序列對龍膽亞族進行系統發育分析,支持喉毛花屬為單系起源,胚胎學研究將喉毛花屬處理為1個獨立的屬,質體基因組與核基因組研究結果進一步支持該屬為1個單系類群[16-20]。目前,喉毛花屬屬級分類問題已明晰,但其屬下物種關系問題尚未得到解決。此外,研究者還發現不同數據集構建的系統發育樹表現出不一致的拓撲結構[16,20],主要表現在皺邊喉毛花(Comastomapolycladum)、鐮萼喉毛花(Comastomafalcatum)等物種上。目前,喉毛花屬下物種關系尚不清晰,有關喉毛花屬下物種的系統發育研究鮮見報道,對其屬下物種的認識有限。因此,后續研究應結合更多豐富、有價值的基因組數據闡明喉毛花屬下物種分類問題。鑒此,文章以皺邊喉毛花及其近緣種作為研究對象,通過生物信息學方法進行基因組比較、進化與密碼子偏好性等分析,構建不同數據集,重建喉毛花屬下物種的系統發育關系,為進一步展開喉毛花屬下物種分化研究提供可靠的分子依據。

1 材料和方法

1.1 材 料

選取12個喉毛花屬植物樣本(5個皺邊喉毛花、2個長梗喉毛花、2個鐮萼喉毛花和2個高杯喉毛花),這些樣本分別采自不同地區與不同群體(表1)。野外采集的新鮮幼葉放入硅膠中干燥,回到實驗室后,置于-20 ℃冰箱保存。憑證標本存放于中國科學院西北高原生物研究所青藏高原生物標本館(HNWP)內。

表1 物種信息與測序結果

1.2 DNA提取與測序

采用改良后的CTAB法從干燥葉片中提取基因組總DNA[21],瓊脂糖凝膠電泳檢測DNA純度和完整性,送至公司構建全基因組小片段文庫。構建完成后使用Agilent 2100對文庫插入片段進行檢測。文庫檢測合格后,利用Illumina高通量測序平臺NovaSeq 6000平臺進行150 bp長度的雙末端測序。

1.3 葉綠體基因組組裝與注釋

利用fastp v 0.23.1[22]對原始數據(raw reads)進行質量控制和過濾(參數設置:-poly_x_min_len:10;cut_front cut_tail cut_window_size:4;qualified_quality_phred:15;low_complexity_filter complexity_threshold:30;length_required:30 thread:4),過濾后的高質量數據(clean reads)用于后續生物信息學分析。使用SPAdes v 3.13.1 與Getorganelle v 1.7.7.0組裝12個個體的葉綠體基因組[23-24],提交至在線工具Geseq進行基因組注釋,得到Genbank文件[25],隨后提交至Sequin軟件,手動去除內含子及外顯子,調整起始/終止密碼子的位置,最終導出用于后續分析的Genbank文件并提交至NCBI數據庫中。此外,從本數據庫中下載已經發表的龍膽科4個屬獐牙菜屬(Swertia)、花錨屬(Halenia)、扁蕾屬(Gentianopsis)與肋柱花屬(Lomatogonium)的9個物種:蒙自獐牙菜(Swertialeducii,登錄號:NC_045301,153 015 bp)、顯脈獐牙菜(Swertianervosa,登錄號:NC_057596,153 690 bp)、祁連獐牙菜(Swertiaprzewalskii,登錄號:MZ261904,153 269 bp)、宿根肋柱花(Lomatogoniumperenne,登錄號:NC_050659,151 678 bp)、橢圓葉花錨(Haleniaelliptica,登錄號:NC_050657,153 305 bp)、Haleniacoreana(登錄號:NC_050657,NC_042674,153 198 bp)、大花扁蕾(Gentianopsisgrandis,登錄號:NC_049879,151 271 bp)、濕生扁蕾(Gentianopsispaludosa,登錄號:NC_050656,151 308 bp)和扁蕾(Gentianopsisbarbata,登錄號:MZ579704,151 123 bp),獐牙菜屬物種作為外類群。

1.4 密碼子偏好性分析

選取4個同科植物(顯脈獐牙菜、宿根肋柱花、扁蕾與橢圓葉花錨)與1個已發表的喉毛花MK331815,與皺邊喉毛花(Chensl-0749)、鐮萼喉毛花(Zhang2019002)、長梗喉毛花(Chensl-0876)和高杯喉毛花(Chen2013572-1)進行密碼子偏好性分析。使用codonW v 1.4.2(http://codonw.sourceforge.net/)計算相對同義密碼子使用度(relative synonymous codon usage,RSCU)、密碼子適應指數(codon adaptation index,CAI)、有效密碼子數(effective number of codons,ENc)、同義密碼子第三位GC含量(GC content of the synonymous third codons,GC3s)、同義第三密碼子胸腺嘧啶含量(synonymous third codon thymine content,T3s)、同義第三密碼子胞嘧啶含量(synonymous third codon cytosine content,C3s)、同義第三密碼子腺嘌呤含量(synonymous third codon adenosine content,A3s)和同義第三密碼子鳥嘌呤含量(homonymous third codon guanine content,G3s)等指數。在TBtools v 1.09中繪制密碼子偏好性分析圖[26]。

1.5 葉綠體基因組比較分析

為比較喉毛花屬及其近緣屬的葉綠體基因組相似性及差異性,利用mVISTA軟件[27]進行可視化分析。以鐮萼喉毛花作為參考,與其余8個物種(顯脈獐牙菜、宿根肋柱花、扁蕾、橢圓葉花錨、皺邊喉毛花、長梗喉毛花、喉毛花及高杯喉毛花)進行比較,選擇Shuffle-LAGAN模型,其他參數設為默認值。利用MISA v 1.0[28]進行短重復序列(short repeat sequence)分析,參數設置:1=10;2=5;3=4;4=3;5=3;6=3,利用Excel 2016繪制SSR數量及其分布情況。

在Phylosuite v 1.2.2[29]中提取9個GB文件的編碼基因(protein-coding gene,PCGs)以及基因間隔區(intergentic spacers,IGS),在MAFFT v 7.313[30]中進行基因及基因間隔區的序列比對,得到的比對結果為輸入文件進行核苷酸多態性分析。將提取的共有基因分為不同的數據集:accD(fatty acid synthesis)、atp(ATP synthase)、cemA(carbon metabolism)、clpP(proteolysis)、infA(translational initiation factor)、matK(RNA processing)、ndh(NADPH dehydrogenase)、pet(cytochrome b/f complex)、psa(photosystem I)、psb(photosystem II)、rbcL(rubisco)、rpl(large subunit of ribosome)、rpo(DNA dependent RNA polymera)、rps(small subunit of ribosome)。利用DNAsp v 6[31]的Nucleotide diversity模塊計算PCGs和IGS的Pi值(參數設置:滑動窗口大小為600 bp,為200 bp),利用Excel 2016繪制結果圖。從結果中篩選高度變異基因間隔區(Pi>0.06)用于系統發育分析。

以顯脈獐牙菜為參考,利用軟件KaKs_Calculator v 2.0[32]計算14個數據集的同義替代率(synonymous,Ks)、非同義替代率(nonsynonymous,Ka)及同義替代率/非同義替代率的比值(Ka/Ks)(method:NG;genetic code:11-Bacterial and plant plastid code)。在R v 4.2.4(https://www.R-project.org/)中利用ggplot2包繪制進化分析圖。

1.6 喉毛花屬下的系統發育分析

共構建了4個數據集:蛋白編碼序列(CDS),第一、二位密碼子位置(CDS1+2),第三位密碼子位置(CDS3)以及篩選出的高度變異基因間隔區(ndhC_trnV-UAC、psbC_trnS-UGA和psbK_psbI)序列用于喉毛花屬下系統發育關系重建。利用MAFFT v 7.313[30]對共有CDS基因以及高度變異IGS進行比對(密碼子比對策略),對24條葉綠體基因組序列共有的82個編碼基因進行多基因聯合分析。CDS基因的密碼子位點數據集采用Partition估算建樹最佳模型,其余數據集在Modelfinder中計算最佳模型[33-34](表2)。在IQ-tree v 1.6.8[35]中進行最大似然(maximum likelihood,ML)分析,選擇Ultrafast bootstrap,快速自然重復設置為1 000次。使用MrBayes v 3.2.6[36]進行貝葉斯(bayesian inference,BI)分析,每次分析進行2次獨立的馬爾可夫鏈蒙特卡羅(MCMC)鏈運行,迭代選擇1 000 000代,每1 000代選1次樣,舍棄前25%的樣本。利用ITOL(interactive tree of life, https://itol.embl.de/)在線軟件對結果圖進行可視化和修改。

表2 所有數據集的最佳建樹模型

2 結果與分析

2.1 葉綠體基因組結構特征

經測序后,對12個樣本的測序質量進行評估,其中,10個樣本的Q30均在90%以上(除高杯喉毛花的2個樣本外),GC含量范圍在35.34%～41.38%之間,基因組長度范圍在150 772～152 093 bp之間(表1)。經組裝及注釋后,獲得12個完整的葉綠體基因組,均由四分體結構組成,包括大單拷貝區(large single copy,LSC)、小單拷貝區(small single copy,SSC)及1對反向重復區(inverted repeat region,IRa、IRb)。選取9個葉綠體基因組進行比較分析,葉綠體基因組長度范圍在151 123～153 690 bp之間(表3)。

表3 9個葉綠體基因組結構信息

2.2 密碼子偏好性分析

基于82個共有基因對喉毛花屬物種的密碼子偏好性進行分析,結果顯示,喉毛花屬與其近緣屬之間略有差異,主要體現在編碼蛋白基因的密碼子數量不同,范圍為50 374～51 230。在20個氨基酸中,亮氨酸(Leu)具有最多的密碼子,而色氨酸(Trp)具有最少的密碼子。RSCU結果表明,有35個密碼子的RSCU值均大于1,其余密碼子均≤1(圖1)。密碼子指數計算結果顯示,所有物種的GC含量相差不大,均在0.39～0.394之間。密碼子適應指數CAI在0.15～0.162之間,有效密碼子數目ENc值在55.89～56.44之間(表4)。

圖1 9個葉綠體基因組的密碼子偏好性熱圖

表4 葉綠體基因組編碼基因的密碼子偏好性指數

2.3 葉綠體基因組比較分析

9個葉綠體基因組序列相似性可視化分析顯示編碼區比非編碼區更保守,相比于反向重復區單拷貝區差異性更大,與核苷酸多態性分析結果基本保持一致(圖2)。SSR分析結果顯示,喉毛花屬及其近緣屬中的重復類型總數在42～53之間,其中扁蕾和顯脈獐牙菜的數量最多。單堿基重復類型最豐富,范圍在28～37之間(圖3)。在10個基因數據集的核苷酸多態性結果中,其Pi值范圍為0.005～0.037。其中,cemA、matK基因具有較高的Pi值,說明它們在喉毛花屬植物中的變異率最高,而基因間隔區的結果中,筆者只展示Pi值大于0.05的基因間隔區,并選取Pi大于0.06的基因間隔區作為數據集構建系統發育樹(圖4)。

圖2 喉毛花屬及其近緣類群的葉綠體基因組比較

圖3 9個葉綠體基因組分SSR分析

圖4 蛋白編碼基因和基因間隔區的核苷酸多態性分析

2.4 葉綠體基因組進化分析

為進一步檢測葉綠體基因組編碼基因受到的選擇壓力,利用9個葉綠體基因組的共有編碼基因,將其分為不同的9個基因數據集,計算同義替代率、非同義替代率以及Ka/Ks值(圖5)。結果表明,非同義替換率均較低,范圍在0.002～0.035之間,而同義替代率的范圍值在0.003～0.11之間。非同義替代率與同義替代率的比值(Ka/Ks)結果中,psbA基因擁有最小值0.005,而ndhB基因的比值(1.24)最大。

圖5 蛋白編碼基因的同義替代率、非同義替代率以及Ka/Ks比率

2.5 喉毛花屬下的系統發育分析

基于CDS、密碼子第一、二位置及第三位置構建的ML和BI系統發育樹具有相似的拓撲結構,所有分支具有較高的支持率(圖6)。從系統發育樹可知,所有的喉毛花屬物種均聚為一支,其姊妹類群為肋柱花屬植物,二者共同與花錨屬的3個物種和獐牙菜屬的1個物種聚為1個大支。而基因間隔區的系統發育樹結果也具有相似的結果(圖7)。

圖6 基于葉綠體基因組蛋白編碼序列,第1、2位密碼子位置與第三位密碼子位置的喉毛花屬下系統發育關系

圖7 基于高度變異基因間隔區聯合喉毛花屬下系統發育關系

從喉毛花屬的分支結果可知,皺邊喉毛花、鐮萼喉毛花、長梗喉毛花均聚在了同一支上且相互嵌套在一起。在分支A中,2個鐮萼喉毛花(QXA0366和Zhang2019356-1)與1個皺邊喉毛花(Chensl-0749)聚在一支(分支a1),與分支b1構成姐妹類群。在分支b1中,皺邊喉毛花(Zhang2019427-1)與鐮萼喉毛花(MK331815)、長梗喉毛花(Zhang2019277)和皺邊喉毛花(Zhang2019004-1)構成姐妹類群。分支B由分支a2和b2構成,二者互為姐妹類群,在分支a2中,皺邊喉毛花(Zhang2018056)與皺邊喉毛花(Zhang2019179-1)、長梗喉毛花(Chensl-0876)聚為一支。在分支b2中,2個高杯喉毛花(Chen2013572-1/2)與鐮萼喉毛花(Zhang2019002)、2個喉毛花(NC_050654、MW324577)互為姐妹類群。基于以上結果可以得知,皺邊喉毛花及其近緣種相互嵌套,未按物種聚類。

此外,聯合變異程度較高的基因間隔區(ndhC_trnV-UAC、psbC_trnS-UGA和psbK_psbI)構建了ML和BI樹,拓撲結構進一步說明喉毛花屬的單系性(圖7)。在2種系統發育樹中,部分物種表現出不一致的拓撲結構。

在最大似然樹中,長梗喉毛花(Chensl-0876)與皺邊喉毛花(Zhang2018056)互為姐妹類群,而在貝葉斯樹中,其與皺邊喉毛花(Zhang2019004-1)聚在一起(圖7)。類似的結果進一步說明皺邊喉毛花、鐮萼喉毛花與長梗喉毛花無法按照物種聚類,與CDS系統發育樹結果一致。

3 討 論

測序技術的迅速發展與測序成本的降低,實現了許多植物的葉綠體基因組測序,其長期作為解決植物類群系統發育問題的有利手段[8-9,37]。對于系統發育位置存在爭議的喉毛花屬而言,葉綠體基因組數據構建出成熟、高分辨率的喉毛花屬系統發育樹,支持了該屬為單系起源。綜合歷史研究分析,喉毛花屬的系統發育位置問題已得到解決,但其屬下物種的系統發育關系尚不明晰[18]。一些系統發育研究案例中,不同數據構建的系統發育樹表現出不一致的拓撲結構,類似情況也在其他植物類群中發現[38-40]。為此,本研究構建了不同數據集重建喉毛花屬下的系統發育關系,以期為喉毛花屬下物種關系的研究提供有價值的分子信息。

對測序原始數據的質控結果顯示,其Q20均在96%以上,Q30值基本在90%以上,進一步說明其測序質量良好,可以進行后續的生物信息學分析。組裝及注釋后的葉綠體基因組為四分體結構,均由1個LSC區、SSC區和1對IR區組成。5個喉毛花屬物種在基因組長度差異較小,長度范圍為從151 526～152 093 bp不等,長梗喉毛花物種的基因組長度最短。同屬物種之間的基因組長度差異是葉綠體基因組研究中最常見的結果,類似的結果也在其他植物中發現[41-42]。這種差異主要與某些基因的擴張和收縮有關,從而影響基因組大小[43]。

葉綠體基因組比較分析表明,喉毛花屬下物種之間的差異較小。mVISTA以及核苷酸多態性結果表明,非編碼區比編碼區的變異程度高,反向重復區比單拷貝區的核苷酸多態性低,這與大多數植物的結果[44-45]保持一致。此外,篩選出變異程度高的基因間隔區或基因可作為進一步研究喉毛花屬物種鑒定及群體遺傳的分子標記候選。SSR結果表明,單堿基重復最為豐富。這與其他青藏高原地區的植物結果保持一致,也與多數植物的結果不一致,或與SSR分析所設定的參數、測序方法及物種本身有關[46]。這些結果可為進一步設計SSR引物及開發高度多態性的SSR標記提供基礎數據。

同義核苷酸替換率和非同義核苷酸替換率的計算,對于植物系統發育的重建具有重要意義,可以判斷植物葉綠體基因組編碼基因是否受到選擇壓力的作用[47-48]。同義核苷酸的替代并不改變蛋白質的結構,而非同義核苷酸的替代卻能改變蛋白質結構,因此非同義替換往往帶來有害性狀而被固定下來[49]。當Ks=Ka,說明編碼基因未受到自然選擇的影響[50];如果Ka/Ks>1,說明基因受正選擇的作用。喉毛花屬植物的進化分析結果顯示,幾乎所有編碼基因的非同義替代率與同義替代率比值(Ka/Ks)均小于1,表明Ks>Ka,進一步說明這些編碼基因受到純化選擇的作用。以上結果與其他植物的結果[11,50]保持一致。本研究篩選的正選擇基因,能夠為進一步研究喉毛花屬下物種的適應性進化提供基礎數據。

密碼子偏好性分析結果表明喉毛花屬的20個氨基酸均由61個密碼子編碼(3個終止密碼子除外)。計算多個衡量指標可以量化植物葉綠體基因組密碼子的偏好性,RSCU可以很直觀地展示密碼子使用的偏好性[51]。在喉毛花屬植物中,有1/2以上密碼子的RSCU值大于1,說明使用這些密碼子的頻率較高。ENc值可以反映出同義密碼子非均衡使用的偏好程度,其結果具有較高的參考價值[52]。當ENc值越小時,其密碼子的偏好程度高。因此,通過比較ENc值可以確定內源基因表達量的相對高低[52]。本研究中,喉毛花屬物種的ENc范圍在55～57之間(均大于35),說明該屬物種的葉綠體基因密碼子偏好性較弱。

喉毛花屬下物種系統發育關系重建結果顯示,基于CDS、密碼子位置與基因間隔區數據集的系統發育樹均表現出高度一致的拓撲結構。該結果進一步支持喉毛花屬植物為單系類群,肋柱花屬植物為該屬的近緣類群,與之前的研究結果[18]保持一致。在支持率良好的喉毛花屬下的系統發育關系結果中,皺邊喉毛花、長梗喉毛花和鐮萼喉毛花各自并未按物種進行聚類,而是相互嵌套,在聯合IGS區構建的系統發育樹中也出現類似結果,進一步支持ITS、matK基因系統發育樹的結果[20]。利用質體基因組數據重建喉毛花屬下物種系統發育關系的研究受到限制,無法明晰皺邊喉毛花及其近緣種之間的關系。在后續的研究中,需要對這些物種展開群體采樣并進行物種分化研究。

本研究利用二代測序技術對皺邊喉毛花及其近緣物種進行了葉綠體基因組的測序與組裝,獲得的測序數據用于分析葉綠體基因組比較、進化、密碼子偏好性等分析,并基于不同數據集和變異程度較高的基因間隔區聯合進行系統發育樹的構建,所呈現的系統發育拓撲結構一致。然而,皺邊喉毛花及其近緣物種未按物種分類,各自并未形成獨立的分支。對于喉毛花屬下的物種分化問題需要開展進一步的群體遺傳學研究,基于現有的成熟技術,揭示屬下物種關系,闡明物種形成機制。