甘新軍,賓 粵,陳煥錦,朱韋光,熊露橋,余恩萍,4,王崢峰**,徐鳳霞,曹洪麟

(1.廣東從化陳禾洞省級(jí)自然保護(hù)區(qū)管理處,廣東廣州 510950;2.中國科學(xué)院華南植物園,廣東省應(yīng)用植物學(xué)重點(diǎn)實(shí)驗(yàn)室,中國科學(xué)院退化生態(tài)系統(tǒng)植被恢復(fù)與管理重點(diǎn)實(shí)驗(yàn)室,廣東廣州 510650;3.華南國家植物園,廣東廣州 510650;4.中國科學(xué)院大學(xué),北京 100049;5.中國科學(xué)院華南植物園,中國科學(xué)院植物資源保護(hù)與可持續(xù)利用重點(diǎn)實(shí)驗(yàn)室,廣東廣州 510650)

木蘭科(Magnoliaceae)是最原始的被子植物,其屬種豐富,類型多樣,是研究被子植物起源和進(jìn)化的重要類群[1-4]。我國是木蘭科起源地和避難所,為木蘭科植物分布中心,擁有很多古老、孑遺和特有物種[2,5-7]。我國的木蘭科物種現(xiàn)主要分布于我國熱帶亞熱帶地區(qū),如云南省、廣西壯族自治區(qū)、貴州省和湖南省。木蘭科中的很多物種樹形優(yōu)美,葉型秀麗、色澤鮮艷,花形態(tài)各異、花色豐富明艷且高貴典雅,觀賞性強(qiáng),是優(yōu)良的綠化樹種[8-11]。木蘭科物種的枝、葉、花含有揮發(fā)性有機(jī)物,作為綠化樹種可以凈化空氣,促進(jìn)人們身心健康[12-15],而且這些揮發(fā)性有機(jī)物與木蘭科產(chǎn)生的其他次生代謝物如生物堿,常被用作中藥[16-18]。另外,木蘭科的一些種類具有較強(qiáng)的光合能力、固碳能力以及土壤改良作用,可用于人工林改造,實(shí)現(xiàn)林業(yè)提質(zhì)增效[19,20];而且一些木蘭科植物樹干通直,木材材質(zhì)均勻、細(xì)密,是很好的用材樹種[ 8,21,22]。同時(shí),木蘭科植物最早記載于我國秦漢時(shí)代,歷經(jīng)千年,寄托了人們對(duì)美好生活的追求和向往,富有深厚的文化沉淀[23]。

厚葉木蓮(Manglietiapachyphylla)為木蘭科木蓮屬(Manglietia)的木本植物,為中國特有種,目前零星分布于我國廣東省和廣西壯族自治區(qū)海拔500 m以上的常綠闊葉林中[24,25]。厚葉木蓮的顯著特征是有光澤的革質(zhì)厚葉,可作為園林綠化和用材物種[24]。厚葉木蓮種群少,個(gè)體數(shù)量小,結(jié)實(shí)率低,種子易被動(dòng)物啃食,自然更新差[25,26],現(xiàn)被列為國家二級(jí)重點(diǎn)保護(hù)野生植物[24]。目前國內(nèi)外已開展厚葉木蓮群落學(xué)[25,26]、花粉形態(tài)[27]和光響應(yīng)生理[28]等方面的研究,但還未有關(guān)于厚葉木蓮的遺傳多樣性及其基因組的研究報(bào)道。基因組及基于基因組開展的物種遺傳多樣性研究,可以揭示物種進(jìn)化過程,從而了解物種的適應(yīng)性并進(jìn)一步應(yīng)用于品種改良[29]。為此,本研究采用二代和三代高通量測(cè)序手段,對(duì)厚葉木蓮基因組進(jìn)行測(cè)序,組裝其基因組草圖,為今后更好地開展其進(jìn)化、遺傳多樣性研究提供參考。

1 材料與方法

1.1 材料

在廣東省廣州市從化區(qū)陳禾洞省級(jí)自然保護(hù)區(qū)內(nèi)選擇1株厚葉木蓮成年大樹(生長(zhǎng)點(diǎn)地理位置為113°55′31.84″E,23°45′1.25″N),其胸高直徑(Diameter at Breast Height,DBH)為22.4 cm,采集其無蟲咬和病斑的樹葉3片。采集的樹葉用剪刀剪碎后,用錫箔紙包裹后立刻投入液氮罐中,之后將樣品送至武漢未來組生物技術(shù)有限公司進(jìn)行高通量測(cè)序。

1.2 方法

測(cè)序包括3個(gè)方面的內(nèi)容:一是采用Nanopore PromethION測(cè)序平臺(tái)對(duì)厚葉木蓮進(jìn)行三代基因組測(cè)序,二是采用MGI DNBSEQ-T7 測(cè)序平臺(tái)對(duì)厚葉木蓮進(jìn)行二代基因組測(cè)序,三是采用MGI DNBSEQ-T7 測(cè)序平臺(tái)對(duì)厚葉木蓮進(jìn)行轉(zhuǎn)錄組測(cè)序,具體實(shí)驗(yàn)流程參考Wang等[30]的研究。針對(duì)所得的測(cè)序數(shù)據(jù),利用不同程序開展數(shù)據(jù)處理和基因組組裝、分析。在數(shù)據(jù)處理過程中主要使用程序的默認(rèn)參數(shù),如在程序運(yùn)行過程中對(duì)默認(rèn)參數(shù)進(jìn)行改動(dòng),則會(huì)在文中具體說明。

1.2.1 測(cè)序數(shù)據(jù)前處理

本研究利用Sickle v1.33 (https://github.com/najoshi/sickle)對(duì)厚葉木蓮二代基因組測(cè)序數(shù)據(jù)進(jìn)行過濾,去除測(cè)序數(shù)據(jù)中堿基質(zhì)量小于30、片段長(zhǎng)度小于80 bp的測(cè)序數(shù)據(jù)。過濾后的二代測(cè)序數(shù)據(jù)用 RECKONER v1.1[31]進(jìn)行糾錯(cuò)。對(duì)于三代基因組測(cè)序數(shù)據(jù),由于低質(zhì)量數(shù)據(jù)在交付前測(cè)序公司已過濾,無需再過濾,本研究只利用Porchop 0.2.4 (https://github.com/rrwick/Porechop)對(duì)厚葉木蓮三代基因組測(cè)序數(shù)據(jù)進(jìn)行接頭過濾。

1.2.2 基因組大小預(yù)測(cè)與組裝

針對(duì)厚葉木蓮二代測(cè)序數(shù)據(jù),利用GenomeScope 2.0[32]進(jìn)行基因組大小預(yù)測(cè)(參數(shù)“-k 21”)。針對(duì)厚葉木蓮三代測(cè)序數(shù)據(jù),本研究選擇大于10 kb測(cè)序序列,利用NextDenovo 2.3.1(https://github.com/Nextomics/NextDenovo)進(jìn)行基因組組裝,組裝的基因組利用Pseudohaploid (https://github.com/schatzlab/pseudohaploid)和Purge_Dups v1.2.6[33]去除冗余序列(如雜合導(dǎo)致的拼接序列),之后利用racon v1.5.0[34],hapo-G v1.3.2[35]和polypolish v0.5.0[36]進(jìn)行組裝序列糾錯(cuò)。針對(duì)組裝完成的基因組,利用BUSCO (Benchmarking Universal Single-Copy Orthologs)v5.4.6[37]對(duì)照“eudicots_odb10.2020-09-10”和“embryophyta_odb10.2020-09-10”兩個(gè)單拷貝基因庫進(jìn)行組裝完整性的評(píng)估。

1.2.3 重復(fù)序列和基因的預(yù)測(cè)、注釋

本研究利用EDTA (Extensive De-novo TE Annotator)v2.1.0[38]和RED (REpeat Detector)v2.0[39]預(yù)測(cè)厚葉木蓮基因組中的重復(fù)序列。利用BEDTools v2.29.2[40]中的“merge”命令將兩個(gè)重復(fù)序列預(yù)測(cè)結(jié)果合并,再利用“maskfasta”命令將預(yù)測(cè)的重復(fù)序列屏蔽掉。針對(duì)屏蔽了重復(fù)序列的基因組,首先利用BRAKER2[41],同時(shí)結(jié)合轉(zhuǎn)錄組數(shù)據(jù)和9個(gè)物種的蛋白質(zhì)序列(表1)預(yù)測(cè)厚葉木蓮基因組組裝序列中的基因;然后將結(jié)果輸入funannotate pipeline v1.8.13(https://github.com/nextgenusfs/funannotate)中,同樣結(jié)合轉(zhuǎn)錄組數(shù)據(jù)和9個(gè)物種的蛋白質(zhì)序列,進(jìn)一步預(yù)測(cè)厚葉木蓮基因組組裝序列中的基因,這一過程包括3個(gè)步驟和命令:“funannotate train”“funannotate predict”和“funannotate update”,在“predict”步驟中使用的參數(shù)為“-max_intronlen 100,000 -busco_db embryophyta -organism other”。

表1 預(yù)測(cè)厚葉木蓮基因組基因的參考物種

基因預(yù)測(cè)結(jié)束后,利用8個(gè)不同的蛋白質(zhì)注釋平臺(tái)對(duì)預(yù)測(cè)的基因進(jìn)行功能分析,包括:dbCAN (DataBase for automated Carbohydrate-active enzyme ANnotation)v10.0[42],eggNOG-mapper (Evolutionary genealogy of genes:Non-supervised Orthologous Groups-mapper)v5.0.2[43],GO (Gene Ontology)[44,45],KEGG (Kyoto Encyclopedia of Genes and Genomes)[46],InterPro (The Integrated Resource of Protein Domains and Functional Sites)v5.60-92.0[47],MEROPS v12.0[48],Pfam (The protein families database)v35.0[49]和SignalP 5.0b[50]。

考慮到注釋的基因中有可變剪切的狀況,在進(jìn)行比較基因組的研究中去除了各物種基因中的可變剪切產(chǎn)生的基因,只保留其中最長(zhǎng)的基因序列進(jìn)行分析。

1.2.4 基因家族和系統(tǒng)發(fā)育分析

利用OrthoFinder 3.0.0[51,52]并結(jié)合其他9個(gè)物種(表1)的蛋白質(zhì)序列進(jìn)行基因家族分析。在OrthoFinder分析過程中,程序自動(dòng)分析獲得物種間的單拷貝基因,并用這些單拷貝基因進(jìn)行系統(tǒng)發(fā)育樹構(gòu)建。利用構(gòu)建的系統(tǒng)發(fā)育樹,使用MCMCTree[53]進(jìn)行物種間分化時(shí)間的估算。在這一分析過程中需要參考已有物種間的分化時(shí)間,本研究從http://timetree.org/獲得這些信息,并在表2中列出。在得到有分化時(shí)間的系統(tǒng)發(fā)育樹后,再利用cafe (Computational Analysis of gene Family Evolution)v5[54]進(jìn)行基因家族的擴(kuò)張和收縮分析;對(duì)其中顯著擴(kuò)張和收縮的基因家族,則利用 TBtools v1.115[55]進(jìn)行GO和KEGG的富集分析。

表2 從http://timetree.org/獲得的物種對(duì)間分化時(shí)間

1.2.5 基因重復(fù)(Gene duplications)

利用wgd v1.1.2[56]進(jìn)行基因組的全基因組重復(fù)(Whole genome duplication)事件分析,該過程是在基因組中進(jìn)行基因間的同源性分析,找到兩兩同源基因,并對(duì)兩兩同源基因進(jìn)行同義突變率(synonymous substitution rate,Ks)計(jì)算,之后查看Ks值的密度分布圖,其峰值出現(xiàn)的地方提示有全基因重復(fù)事件發(fā)生;進(jìn)一步利用 Doubletrouble v0.99.1 (https://github.com/almeidasilvaf/doubletrouble)開展全基因組重復(fù)基因、串聯(lián)重復(fù)(Tandem duplications)基因、近端重復(fù)(Proximal duplications)基因、轉(zhuǎn)座重復(fù)(Transposed duplications)基因、散在重復(fù)(Dispersed duplications)基因的分析[55]。其中,串聯(lián)重復(fù)基因是指彼此連續(xù)或中間間隔不超過5個(gè)其他基因的近緣基因;近端重復(fù)基因是連續(xù)分布、中間間隔不超過10個(gè)其他基因的近緣基因;轉(zhuǎn)置重復(fù)基因是由轉(zhuǎn)座子介導(dǎo)的重復(fù)基因;散在重復(fù)基因是隨機(jī)分布、彼此間不靠近的重復(fù)基因[57]。利用TBtools分別對(duì)全基因組重復(fù)基因、串聯(lián)重復(fù)基因、近端重復(fù)基因進(jìn)行GO和KEGG富集分析。

2 結(jié)果與分析

2.1 厚葉木蓮基因組測(cè)序與組裝

本研究中三代基因組測(cè)序共獲得大約118.1 Gb測(cè)序數(shù)據(jù),二代基因組測(cè)序獲得約264.1 Gb測(cè)序數(shù)據(jù),轉(zhuǎn)錄組測(cè)序獲得約31.4 Gb測(cè)序數(shù)據(jù)?；蚪M和轉(zhuǎn)錄組原始測(cè)序數(shù)據(jù)已上傳至GenBank,三代基因組測(cè)序數(shù)據(jù)序列號(hào)為SRR24423593、SRR24423594、SRR24423595;二代基因組測(cè)序數(shù)據(jù)序列號(hào)為SRR24390471、SRR24390472、SRR24390473;轉(zhuǎn)錄組測(cè)序數(shù)據(jù)序列號(hào)為SRR24415003。利用GenomeScope 2.0分析厚葉木蓮基因組大小為1 969 269 649 bp。

本研究利用NextDenovo 2.3.1進(jìn)行厚葉木蓮基因組組裝,得到基因組組裝大小為2 350 821 062 bp,包含1 436個(gè)拼接序列(contig),N50(將組裝的序列按照長(zhǎng)度由大到小進(jìn)行累加,當(dāng)累加到某個(gè)序列時(shí),累加的值為基因組50%的長(zhǎng)度時(shí),此序列的長(zhǎng)度即為N50)為6 839 193 bp,最長(zhǎng)拼接序列為26 073 671 bp,最短拼接序列為16 675 bp,平均為1 637 062.0 bp。去除冗余序列和糾錯(cuò)后,最終的基因組組裝大小為2 092 298 891 bp,包含676個(gè)拼接序列,N50為7 961 115 bp,最長(zhǎng)拼接序列為26 180 362 bp,最短拼接序列為27 281 bp,平均為3 095 117 bp,組裝的其他統(tǒng)計(jì)信息見表3。組裝的厚葉木蓮基因組數(shù)據(jù)上傳至GenBank,序列號(hào)為JASAUF000000000。

表3 厚葉木蓮基因組組裝結(jié)果

BUSCO對(duì)組裝的基因組完整性的評(píng)估顯示,比對(duì)“eudicots_odb10.2020-09-10”單拷貝基因庫,全部2 326個(gè)BUSCO單拷貝基因,96.6% BUSCO單拷貝基因在厚葉木蓮基因組中完整匹配。其中,能完整匹配到BUSCO單拷貝基因庫而且在厚葉木蓮基因組中也是單拷貝的基因占90.8%,能完整匹配到BUSCO單拷貝基因庫但在厚葉木蓮基因組中為多拷貝的基因占5.8%;不完整匹配到BUSCO單拷貝基因庫的基因有27個(gè),占1.2%;有54個(gè)BUSCO單拷貝基因未在厚葉木蓮基因組的基因中匹配到,占2.3%。比對(duì)“embryophyta_odb10.2020-09-10”單拷貝基因庫,全部1 614個(gè)BUSCO單拷貝基因中,98.8% BUSCO單拷貝基因在厚葉木蓮基因組中完整匹配,能完整匹配到BUSCO單拷貝基因庫而且在厚葉木蓮基因組中也是單拷貝的基因占95.3%,能完整匹配到BUSCO單拷貝基因庫但在厚葉木蓮基因組中為多拷貝的基因占3.5%;不完整匹配到BUSCO單拷貝基因庫的基因有12個(gè),占0.7%;有7個(gè)BUSCO單拷貝基因未在厚葉木蓮基因組的基因中匹配到,占0.4%。

2.2 厚葉木蓮基因組注釋

在厚葉木蓮基因組中,RED和EDTA分別檢測(cè)到1 342 604 397 bp (64.16%)和1 447 436 839 bp(69.17%)重復(fù)序列,兩者合并后得到1 601 108 919 bp的重復(fù)序列,占基因組的76.52%。EDTA檢測(cè)結(jié)果表明,厚葉木蓮基因組中重復(fù)序列最多的是Gypsy類的長(zhǎng)末端重復(fù)序列(Long terminal repeat),有475 508 695 bp,占基因組的22.73%(表4)。

表4 EDTA檢測(cè)的重復(fù)序列

本研究對(duì)厚葉木蓮基因組進(jìn)行基因預(yù)測(cè)共獲得37 900個(gè)基因,這些基因共編碼了41 675種蛋白質(zhì)。利用數(shù)據(jù)庫對(duì)這些蛋白質(zhì)序列進(jìn)行功能注釋,其中32 249個(gè) (77.4%)蛋白質(zhì)序列注釋到了8個(gè)蛋白質(zhì)注釋平臺(tái)數(shù)據(jù)庫中的一個(gè),具體如下:17 886個(gè)蛋白質(zhì)序列注釋到GO數(shù)據(jù)庫,21 430個(gè)蛋白質(zhì)序列注釋到InterPro數(shù)據(jù)庫,31 146個(gè)蛋白質(zhì)序列注釋到eggNOG-mapper數(shù)據(jù)庫,23 203個(gè)蛋白質(zhì)序列注釋到Pfam數(shù)據(jù)庫,1 160個(gè)蛋白質(zhì)序列注釋到dbCAN數(shù)據(jù)庫,1 016個(gè)蛋白質(zhì)序列注釋到MEROPS數(shù)據(jù)庫,3 078個(gè)蛋白質(zhì)序列注釋到SignalP數(shù)據(jù)庫,15 752個(gè)蛋白質(zhì)序列注釋到KEGG數(shù)據(jù)庫?；蜃⑨屛募焉蟼髦林袊茖W(xué)院華南植物園數(shù)據(jù)存儲(chǔ)庫(https://doi.org/10.57841/casdc.0001214)。

2.3 厚葉木蓮基因家族與基因組進(jìn)化分析

本研究對(duì)包括厚葉木蓮在內(nèi)的10個(gè)物種的蛋白質(zhì)序列進(jìn)行了基因家族分析,結(jié)果共得到24 616個(gè)基因家族。從表5統(tǒng)計(jì)結(jié)果可以看出,厚葉木蓮基因組中有34 319個(gè)基因歸為其中一個(gè)基因家族,占總基因的90.6%。對(duì)于木蘭科其他兩個(gè)物種鵝掌楸和望春玉蘭,它們分配到基因家族的基因比例分別為95.7%和93.4%,均高于厚葉木蓮;但厚葉木蓮基因分布于15 894個(gè)基因家族中,占全部基因家族的64.6%;而鵝掌楸和望春玉蘭的基因出現(xiàn)在13 858和14 935個(gè)基因家族中,占全部基因家族的56.3%和60.7%,低于厚葉木蓮。24 616個(gè)基因家族中有710個(gè)厚葉木蓮特有的基因家族,這些家族包含了3 373個(gè)厚葉木蓮基因。鵝掌楸和望春玉蘭特有的基因家族分別是437和999個(gè),雖然一個(gè)低于厚葉木蓮特有基因家族數(shù)目,一個(gè)高于厚葉木蓮特有基因家族數(shù)目,但是它們所包含的基因數(shù)(分別為6 845和4 543)均高于厚葉木蓮,而且這些特有基因數(shù)占兩個(gè)物種各自所有基因的比例也高于厚葉木蓮,特別是鵝掌楸,為19.4%(表5,圖1)。3種木蘭科植物共同出現(xiàn)的基因家族有345個(gè)(圖1)。

圖1 厚葉木蓮與其他物種基因家族交集圖

表5 基因家族及其基因信息統(tǒng)計(jì)

對(duì)厚葉木蓮特有的基因家族進(jìn)行GO和KEGG富集分析,結(jié)果表明這些基因家族所包含的基因的主要功能與細(xì)胞內(nèi)穩(wěn)態(tài)(GO:0055067、GO:0006885、GO:0030433、GO:0036503)、葡聚糖代謝(Glucan metabolic process)(GO:0044042、GO:0016762、GO:0046527)、原核生物抗性(Prokaryotic defense system)、苯丙氨酸代謝(Phenylalanine metabolism)和轉(zhuǎn)錄(Transcription)等有關(guān)(https://doi.org/10.6084/m9.figshare.23690001.v3,表S1和S2)。

系統(tǒng)發(fā)育分析表明厚葉木蓮與望春玉蘭聚在一起(圖2),兩個(gè)物種是在約10 500 000年前(95%CI:4 989 810-17 055 600年)從共同的祖先開始分化。對(duì)基因家族擴(kuò)張和收縮分析結(jié)果表明,厚葉木蓮基因組中有686個(gè)基因家族表現(xiàn)為收縮,1 295個(gè)基因家族表現(xiàn)為擴(kuò)張,其中有136個(gè)基因家族表現(xiàn)為顯著收縮,417個(gè)基因家族表現(xiàn)為顯著擴(kuò)張。GO和KEGG富集分析表明,厚葉木蓮基因組中顯著擴(kuò)張的基因家族與木質(zhì)部、韌皮部發(fā)育(GO:0010088,GO:0010087)、肌動(dòng)蛋白絲(Actin filament,GO:0030837、GO:0030833、GO:0061572、GO:0051017、GO:0030832、GO:0008064、GO:0030041)、細(xì)胞運(yùn)動(dòng)(Cell motility)、油菜素甾醇生物合成(Brassinosteroid biosynthesis)、類黃酮生物合成(Isoflavonoid biosynthesis)、維生素B6代謝(Vitamin B6 metabolism)、硫代葡萄糖苷(Glucosinolate biosynthesis)和萜類物質(zhì)生物合成有關(guān)(https://doi.org/10.6084/m9.figshare.23690001.v3,表S3和表S4);厚葉木蓮基因組中顯著收縮的基因家族主要與木質(zhì)素(Lignin)代謝(GO:0009808、GO:0046274)、類黃酮生物合成,以及二苯乙烯、二芳基庚酸類和姜酚生物合成(Stilbenoid,diarylheptanoid and gingerol biosynthesis)相關(guān)(https://doi.org/10.6084/m9.figshare.23690001.v3,表S5和表S6)。

Divergence times are shown below the tree;The "-/+" and the numbers beside the tree nodes and species represent the number of contracted and expanded gene families in Manglietia pachyphylla and other species.

2.4 厚葉木蓮全基因組和基因重復(fù)分析

全基因組重復(fù)事件分析表明,厚葉木蓮和鵝掌楸、望春玉蘭表現(xiàn)出同樣的Ks峰型(圖3),因此它們近期共同經(jīng)歷了一次全基因組重復(fù)事件。對(duì)厚葉木蓮基因組的基因重復(fù)研究表明,厚葉木蓮基因組中有4 769個(gè)基因與其全基因組重復(fù)有關(guān)(12.6%),3 317個(gè)基因?yàn)榇?lián)重復(fù)(8.8%),3 124個(gè)基因?yàn)榻酥貜?fù)(8.2%),136個(gè)基因?yàn)檗D(zhuǎn)置重復(fù)(0.4%),18 522個(gè)基因?yàn)樯⒃谥貜?fù)(48.9%)。GO和KEGG富集分析表明,厚葉木蓮基因組中全基因組重復(fù)基因主要功能與DNA內(nèi)復(fù)制調(diào)控(Regulation of DNA endoreduplication)、植物晝夜節(jié)律(Circadian rhythm-plant)、鳥嘌呤核苷酸結(jié)合蛋白(GTP-binding proteins)、檸檬酸循環(huán)[Citrate cycle (TCA cycle)],蛋白酶體(Proteasome)、轉(zhuǎn)錄因子(Transcription factors)(https://doi.org/10.6084/m9.figshare.23690001.v3,表S7和表S8)相關(guān)。GO和KEGG富集分析表明,厚葉木蓮基因組中串聯(lián)重復(fù)基因與NADH氧化(NADH oxidation,GO:0006116)、寡肽轉(zhuǎn)運(yùn)(Oligopeptide transport,GO:0006857)、RNA脫帽(RNA decapping,GO:0110154)、谷胱甘肽代謝過程(Glutathione metabolic process,GO:0006749)、過氧化氫分解代謝過程(Hydrogen peroxide catabolic process,GO:0042743、GO:0042744)、類黃酮生物合成、牛磺酸和低?；撬岽x(Taurine and hypotaurine metabolism)、苯并惡唑嗪酮類化合物生物合成(Benzoxazinoid biosynthesis)、玉米素生物合成(Zeatin biosynthesis)、萜類生物合成(Terpenoid biosynthesis)、苯丙烷類生物合成(Phenylpropanoid biosynthesis)、生物堿合成等相關(guān)(https://doi.org/10.6084/m9.figshare.23690001.v3,表S9和表S10)。厚葉木蓮基因組中近端重復(fù)基因主要和防御(GO:0043207、GO:0098542、GO:0051707、GO:0009607、GO:0050832、GO:0009605、GO:0006952)、次生代謝(GO:0044550、GO:0019748)、各類生物堿生物合成(Biosynthesis of various alkaloids)、花青素生物合成(Anthocyanin biosynthesis)、聚酮生物合成(Polyketide biosynthesis)、苯并惡唑嗪酮類化合物生物合成、萜類生物合成、硫代葡萄糖苷生物合成和油菜素甾醇生物合成等相關(guān)(https://doi.org/10.6084/m9.figshare.23690001.v3,表S11和表S12)。

圖3 全基因組重復(fù)分析中同義突變率(Ks)的密度分布

3 討論

3.1 厚葉木蓮的基因組組裝

高通量測(cè)序已成為研究物種基因組和遺傳多樣性的主要手段。目前高通量測(cè)序分為短片段測(cè)序和長(zhǎng)片段測(cè)序量?jī)煞N方式,前者測(cè)序平臺(tái)主要包括Illumia和MGI測(cè)序公司一系列儀器設(shè)備,后者測(cè)序平臺(tái)主要包括Nanopore和PacBio測(cè)序公司一系列儀器設(shè)備。短片段高通量測(cè)序的優(yōu)勢(shì)是測(cè)序的準(zhǔn)確性較高、數(shù)據(jù)量大和價(jià)格便宜;長(zhǎng)片段測(cè)序的優(yōu)勢(shì)是測(cè)得的序列長(zhǎng)度長(zhǎng)、完整性好,有利于后續(xù)基因組組裝的連續(xù)性,如Nanopore測(cè)序平臺(tái)可測(cè)100 kb或更長(zhǎng)的序列片段。但采用Nanopore和PacBio測(cè)序平臺(tái)測(cè)得的序列錯(cuò)誤率較高,為此PacBio公司推出了采用Hifi (High fidelity reads)模式的測(cè)序方式,使得測(cè)序錯(cuò)誤率大為降低,但是測(cè)序長(zhǎng)度有所限制(15-20 kb)。為得到高質(zhì)量基因組,可采用不同測(cè)序平臺(tái)測(cè)序,再進(jìn)行基因組組裝,以實(shí)現(xiàn)優(yōu)勢(shì)互補(bǔ)。本研究采用Nanopore長(zhǎng)片段測(cè)序平臺(tái)的測(cè)序結(jié)果進(jìn)行厚葉木蓮的基因組初步組裝,再用MGI短片段測(cè)序平臺(tái)的測(cè)序結(jié)果對(duì)組裝的基因組進(jìn)行糾錯(cuò),得到了較完整和準(zhǔn)確的組裝結(jié)果,但是還沒有組裝到染色體級(jí)別,后期本研究將繼續(xù)采用Hi-C(High-throughput chromosome conformation capture)測(cè)序方式進(jìn)一步對(duì)厚葉木蓮基因組進(jìn)行測(cè)序,測(cè)得的結(jié)果可用于進(jìn)一步組裝以完善厚葉木蓮基因組。

在采用Nanopore測(cè)序數(shù)據(jù)進(jìn)行厚葉木蓮基因組組裝過程中,本研究采用NextDenovo 2.3.1軟件進(jìn)行序列拼接,這一軟件采用先糾錯(cuò)再組裝(Correction then assembly)的模式進(jìn)行基因組組裝,保證了組裝序列的連續(xù)性和正確性,優(yōu)于其他組裝軟件。同時(shí)本研究利用多種方法去除了原始組裝序列中的冗余序列并進(jìn)行后期糾錯(cuò):其中racon v1.5.0和hapo-G v1.3.2分別使用Nanopore測(cè)序獲得的長(zhǎng)片段和MGI測(cè)序獲得的短片段進(jìn)行組裝基因組內(nèi)堿基和插入缺失錯(cuò)誤的糾錯(cuò),polypolish v0.5.0主要針對(duì)組裝序列中的同聚體長(zhǎng)度(Polypolish-length,如“AAAAAA”這種重復(fù)序列的長(zhǎng)度)進(jìn)行糾錯(cuò),3個(gè)糾錯(cuò)軟件的使用可以很好地保證組裝基因組的準(zhǔn)確性。

在基因組注釋環(huán)節(jié),本研究利用EDTA和RED軟件相互組合的方法進(jìn)行基因組中重復(fù)序列的查找。EDTA主要用于基因組中轉(zhuǎn)座子的查找,而RED查找包括轉(zhuǎn)座子在內(nèi)的所有重復(fù)序列類型,如微衛(wèi)星體(Microsatellite)等,兩者結(jié)合保證了重復(fù)序列查找的全面性,但RED分析的結(jié)果并沒有對(duì)不同重復(fù)序列進(jìn)行歸類。本研究首先使用BRAKER2,結(jié)合轉(zhuǎn)錄組數(shù)據(jù)和其他物種的蛋白質(zhì)序列預(yù)測(cè)厚葉木蓮基因組中的基因,得到初步基因序列注釋結(jié)果。由于這一結(jié)果還存在較多錯(cuò)誤預(yù)測(cè),本研究進(jìn)一步利用funannotate軟件對(duì)初步預(yù)測(cè)的結(jié)果進(jìn)行整合。Funannotate是個(gè)軟件整合工具,被編譯為一個(gè)管段(Pipeline)流程。該軟件可以進(jìn)行不同方式的基因預(yù)測(cè)(包括de novo注釋、轉(zhuǎn)錄組注釋和同源蛋白質(zhì)注釋),再結(jié)合其他分析軟件的結(jié)果(本研究中為BRAKER2基因預(yù)測(cè)結(jié)果),可以獲得統(tǒng)一、可靠的高質(zhì)量基因預(yù)測(cè)結(jié)果。

3.2 厚葉木蓮的比較基因組學(xué)研究

本研究中厚葉木蓮的基因組組裝大小為2 092 298 891 bp(約2.09 Gb),大于木蘭科的日本厚樸(Magnoliahypoleuca,修正名為Houpoeaobovata)基因組(約1.64 Gb)[58]、厚樸(Magnoliaofficinalis,修正名為H.officinalis)的基因組(約1.68 Gb)[59]和鵝掌楸的基因組(約1.74 Gb)[60],小于望春玉蘭的基因組(約2.22 Gb)[29]。基于“embryophyta_odb10.2020-09-10”單拷貝基因庫,厚樸基因組的組裝完整性評(píng)估為86.20%的BUSCO單拷貝基因能被完整匹配到,日本厚樸基因組的組裝完整性為98.6%,望春玉蘭基因組的組裝完整性為95.7%(上述數(shù)值為本研究利用望春玉蘭基因組分析的結(jié)果),鵝掌楸基因組的組裝完整性為98.8%(此值為本研究利用鵝掌楸基因組重新分析的結(jié)果),而厚葉木蓮的BUSCO組裝完整性為98.8%,與鵝掌楸的結(jié)果相同,但高于厚樸、日本厚樸和望春玉蘭3個(gè)物種。從重復(fù)序列比例看,厚葉木蓮基因組中重復(fù)序列占基因組的76.5%,高于望春玉蘭(66.48%)[29]、日本厚樸(64.54%)[58]以及鵝掌楸(63.81%)[60],但小于厚樸(81.44%)[57],厚葉木蓮基因組重復(fù)序列比例在已有木蘭科物種基因組重復(fù)序列比例的范圍內(nèi)。上述結(jié)果說明,雖然本研究報(bào)道的厚葉木蓮基因組還是草圖狀態(tài),沒有組裝到染色體級(jí)別,但是組裝已經(jīng)很完整,結(jié)果可靠。

對(duì)基因家族的研究表明,厚葉木蓮中與次生代謝相關(guān)的基因家族顯著擴(kuò)張,如萜類物質(zhì),這與日本厚樸、望春玉蘭基因家族分析結(jié)果相似[29,58]?；虮磉_(dá)分析表明,望春玉蘭花中萜類相關(guān)基因的表達(dá)高于其葉子,說明萜類物質(zhì)是望春玉蘭花香的主要原因[29],厚葉木蓮花是否也有類似的結(jié)果還需結(jié)合基因組進(jìn)一步研究確定。除了次生代謝物質(zhì)相關(guān)的基因家族在厚葉木蓮基因組中有顯著擴(kuò)張外,本研究還發(fā)現(xiàn)與木質(zhì)部、韌皮部發(fā)育相關(guān)的基因家族,以及與光合作用[光合作用光反應(yīng)(Photosynthesis light reaction,GO:0019684)、光合作用光收獲(Photosynthesis light harvesting,GO:0009765)以及KEGG中的光蛋白(Photosynthesis proteins)]、溫度[熱反應(yīng)(Response to heat,GO:0009408)、溫度刺激反應(yīng)(Response to temperature stimulus,GO:0009266)]等相關(guān)的基因家族也有擴(kuò)張。厚葉木蓮生長(zhǎng)在較高海拔的地區(qū),木材致密通直[26],這些擴(kuò)張的基因家族是否與其適應(yīng)高海拔的溫度、光照、強(qiáng)風(fēng)等生境有關(guān)還需進(jìn)一步研究。厚葉木蓮基因組中的擴(kuò)張基因家族中的一部分與肌動(dòng)蛋白絲(也稱“微絲”)相關(guān),而肌動(dòng)蛋白是植物細(xì)胞骨架的主要元素,相關(guān)研究表明它與植物抗真菌功能密切相關(guān)[61,62]。

對(duì)日本厚樸、望春玉蘭的研究均表明木蘭科都共同經(jīng)歷了一次近期的全基因組重復(fù)事件[29,58]。厚葉木蓮有著與日本厚樸相似的全基因組重復(fù)基因、串聯(lián)重復(fù)基因和近端重復(fù)基因比例。與全基因組重復(fù)相關(guān)的基因占厚葉木蓮基因組基因的比例為12.6%,日本厚樸為13.4%[58];串聯(lián)重復(fù)基因在厚葉木蓮基因組基因的比例為8.8%,日本厚樸為7.6%[58];近端重復(fù)基因在厚葉木蓮基因組基因的比例為8.2%,日本厚樸為9.4%[58]。基因串聯(lián)重復(fù)和近端重復(fù)是基因形成的重要方式,與物種適應(yīng)性密切相關(guān)[57]。日本厚樸中串聯(lián)重復(fù)和近端重復(fù)基因的功能主要與苯丙烷類、萜類、類黃酮的生物合成等有關(guān)[58],這與厚葉木蓮這兩類基因的研究結(jié)果相同,但厚葉木蓮可產(chǎn)生更多的次生代謝合成物,如苯并惡唑嗪酮類化合物、生物堿、聚酮等。厚樸基因組中,萜類合成相關(guān)的基因也呈串聯(lián)重復(fù)狀況[59]。值得注意的是日本厚樸原產(chǎn)日本,具有很好的抗寒性,研究發(fā)現(xiàn)苯丙烷類生物合成與日本厚樸抗寒性相關(guān)。

對(duì)厚葉木蓮中與全基因組重復(fù)相關(guān)基因的富集分析表明,這些基因與植物最基本的生長(zhǎng)發(fā)育密切相關(guān),如植物晝夜節(jié)律相關(guān)基因[63]、與碳代謝相關(guān)的檸檬酸循環(huán)基因[64]、參與細(xì)胞的多種生命活動(dòng)(如細(xì)胞通訊、核糖體與內(nèi)質(zhì)網(wǎng)的結(jié)合、小泡運(yùn)輸、蛋白質(zhì)合成等)的GTP結(jié)合蛋白基因[65,66]等,說明全基因組重復(fù)在木蘭科植物適應(yīng)性方面具有重要作用。

4 結(jié)論

本研究首次報(bào)道了木蘭科木蓮屬物種的基因組,這一結(jié)果對(duì)全面深入了解木蘭科物種以及厚葉木蓮的進(jìn)化、適應(yīng)具有重要作用。本研究組裝的厚葉木蓮基因組大小為2 092 298 891 bp,基因組序列中76.5%為重復(fù)序列。通過基因預(yù)測(cè),在組裝的厚葉木蓮基因組中共注釋到37 900個(gè)基因,它們編碼了41 675種蛋白質(zhì)。對(duì)厚葉木蓮基因組研究發(fā)現(xiàn),全基因組重復(fù)與木蘭科植物進(jìn)化適應(yīng)性密切相關(guān),很多與植物基本生長(zhǎng)發(fā)育相關(guān)的基因通過全基因組重復(fù)在木蘭科得到加強(qiáng);厚葉木蓮富含有與次生代謝相關(guān)的多種基因,由此產(chǎn)生的次生代謝物質(zhì)有利于厚葉木蓮在高海拔生長(zhǎng)以及抵御病蟲害,但其(種子)香味也使得其易受啃食傷害。因此,厚葉木蓮基因組的組裝為從遺傳角度深入了解其瀕危機(jī)制提供了可能,也為科學(xué)合理利用厚葉木蓮以及提取其次生物質(zhì)作為生物醫(yī)藥提供了參考。