新暗色柱節(jié)孢H1全基因組測序和注釋

趙曉珍,齊 勇,王紅林,馬玉華,鄭乾明

(貴州省農(nóng)業(yè)科學(xué)院果樹科學(xué)研究所,貴陽 550006)

【研究意義】火龍果(Hylocereusspp.Britt.& Rose)屬仙人掌科(Cactaceae)量天尺屬(Hylocereus)多年生果樹,又名紅龍果、青龍果或仙蜜果等?；瘕埞a(chǎn)于中南美洲[1],目前在馬來西亞、越南和中國等的熱帶和亞熱帶地區(qū)種植,在我國主要分布在海南、廣西、廣東、臺灣、貴州和云南等地。目前,國內(nèi)外報道較多的火龍果真菌病害主要有潰瘍病[2-8]、炭疽病[9-10]等。潰瘍病是最嚴重的病害之一,主要危害植株的莖和果實,嚴重降低果實商品價值。探討火龍果潰瘍病病原的基因組信息,有助于了解其遺傳變異和感染機制,從而制定有效的防控策略和培育抗性品種。【前人研究進展】火龍果受潰瘍病侵染后,發(fā)病初期在莖和果實表面形成褪綠小斑點,嚴重時病斑連成片,后期病斑突起;濕度大時病斑邊緣形成水漬狀,不斷擴大和腐爛,果實外觀斑駁。近年來,火龍果潰瘍病在馬來西亞[2]、以色列[3]、美國佛羅里達州[4]等和中國臺灣[5]、廣東[6]、海南[7]及廣西[11]等地均有發(fā)生。該病在廣西和海南的發(fā)病率分別為50%以上和100%,在美國佛羅里達州發(fā)病率達70%以上[4,7,11]。自2012年在臺灣首次報道引起火龍果潰瘍病的病原為新暗色柱節(jié)孢(Neoscytalidiumdimidiatum)[5]以來,我國廣東、廣西和海南等地及馬來西亞、以色列和美國等也陸續(xù)報道由新暗色柱節(jié)孢引起的火龍果潰瘍病[2-4,6-7]。新暗色柱節(jié)孢屬葡萄座腔菌科(Botryosphaeriaceae),其異名包括Fusicoccumdimidiatum、Hendersonulatoruloides及Scytalidiumdimidiatum[12],可引起蘭花葉斑病[13]和李、杏枯萎病[14-15]等植物病害。除侵染植物外,新暗色柱節(jié)孢還引起皮膚病和鼻竇炎等人類疾病[16-17]。近年來,以Pacbio和Nanopore平臺為代表的單分子實時測序技術(shù)被廣泛用于基因組測序和組裝[18-20]。Illumina具有測序短讀長和高準確率等優(yōu)點[19],有學(xué)者分別采用Pacbio+Illumina混合策略對金針菇(Flammulinafiliformis)[21]、木賊鐮孢菌(Fusariumequiseti)D25-1[22]、菜豆殼球孢(M.phaseolina)11-12和Al-1[23]、葡萄座腔菌(Botryosphaeriadothidea)sdau11-99[24]和Nanopore+Illumina混合策略對Laburnicolarhizohalophila[25]和尖孢鐮刀菌(Fusariumoxysporum)[26]進行研究,均獲得高質(zhì)量的參考基因組序列。【本研究切入點】關(guān)于新暗色柱節(jié)孢的致病機制目前鮮有研究,也缺乏高質(zhì)量的參考基因組;且尚無分離自火龍果的新暗色柱節(jié)孢基因組序列報道?！緮M解決的關(guān)鍵問題】利用前期分離獲得的火龍果潰瘍病病原菌H1菌株,經(jīng)純培養(yǎng)、形態(tài)學(xué)和分子生物學(xué)鑒定為新暗色柱節(jié)孢。采用Nanopore PromethION+ Illumina NovaSeq混合策略對火龍果潰瘍病病原菌新暗色柱節(jié)孢H1菌株開展基因組測序、組裝和注釋分析,為了解其基因組構(gòu)成、遺傳變異、致病機理和防控策略奠定基礎(chǔ)。

1 材料與方法

1.1 材料

1.1.1 菌株 火龍果潰瘍病病原菌H1菌株分離自貴州省鎮(zhèn)寧縣種植的火龍果品種黔蜜龍植株莖部,經(jīng)鑒定為新暗色柱節(jié)孢(N.dimidiatum)[10]。

1.1.2 試劑 文庫構(gòu)建試劑盒(Nextera DNA Flex Library Prep,Illumina,USA),連接試劑盒(SQK-LSK109,Oxford Nanopore Technologies,UK)。

1.1.3 儀器設(shè)備 超聲波破碎儀(M220,Covaris,USA),熒光計(Qubit 2.0,Life Technologies,USA),生物分析儀(Agilent 2100,Agilent Technologies,USA),全自動核酸剪切儀(Megaruptor 2,Diagenode,Belgium),全自動核酸回收儀(BluePippin,Sage Science,USA),Oxford Nanopore測序儀(PromethION,Oxford Nanopore Technologies UK),Illumina NovaSeq平臺(Illumina,San Diego,CA,USA)。

1.1.4 軟件 GenomeScope分析平臺(https://github.com/schatzlab/genomescope),NECAT軟件(https://github.com/xiaochuanle/NECAT),Racon軟件(v.1.4.11,https://github.com/isovic/racon),Pilon軟件(v.1.23,https:// github.com/broadinstitute/pilon),purge_haplogs軟件(https://github.com/skingan/purge_haplotigs_multiBAM),BWA軟件(v.0.7.17,https://github.com/lh3/bwa),BUSCO軟件(v.3.0.1,https://github.com/RoyNexus/busco),RepeatModeler軟件(v.1.0.4,https://github.com/rmhubley/RepeatModeler),RepeatMasker軟件(v.4.0.5,http://www.repeatmasker.org/),BRAKER軟件(v.2.1.4,https://github.com/GaiusLatestVersionAugustus/ BRAKER)。

1.2 方法

1.2.1 菌株培養(yǎng)和基因組DNA提取 將H1菌株接種于PDA培養(yǎng)基,28 ℃培養(yǎng)5 d。收集菌絲后液氮速凍,采用CTAB法[27]提取基因組DNA。

1.2.2 Illumina測序和K-mer分析 基因組DNA樣品使用超聲波破碎儀處理,然后使用文庫構(gòu)建試劑盒構(gòu)建長度為150 bp的小片段插入文庫。測序文庫使用熒光計定量并稀釋,再使用生物分析儀檢測。合格文庫使用Illumina NovaSeq平臺測序。獲得Raw reads后,過濾低質(zhì)量的序列。過濾標準如下:去除不確定堿基含量超過5%的Raw reads;去除質(zhì)量值不大于5的堿基含量達50%的Raw reads;去除含有接頭污染的Raw reads;去除PCR擴增過程中產(chǎn)生的重復(fù)序列。過濾后的高質(zhì)量序列使用GenomeScope進行K-mer分析,以估算基因組的大小、雜合率和重復(fù)片段含量等信息。

1.2.3 Nanopore測序 基因組DNA樣品使用全自動核酸剪切儀隨機打斷,磁珠富集和純化大片段DNA。使用全自動核酸回收儀回收大片段,進行損傷修復(fù)、純化、末端修復(fù)并加A。利用連接試劑盒連接,使用熒光計定量。測序文庫加入R9.4流動槽,轉(zhuǎn)移到PromethION測序儀測序。獲得原始序列后,過濾去除平均質(zhì)量值≤7的序列。以上實驗操作委托武漢貝納科技服務(wù)有限公司開展。

1.2.4 基因組組裝 首先利用Nanopore PromethION平臺測序獲得的高質(zhì)量序列,使用NECAT軟件進行序列糾錯和基因組組裝,獲得初步組裝結(jié)果。然后使用Racon軟件結(jié)合Nanopore數(shù)據(jù)對拼接結(jié)果進行2輪自我糾錯;再使用Pilon軟件對Illumina數(shù)據(jù)進行2輪糾錯。最后使用purge_haplogs去除雜合獲得組裝結(jié)果。使用BWA程序?qū)llumina數(shù)據(jù)比對到基因組。使用BUSCO程序選取真核和真菌直系同源基因評估完整性。

1.2.5 基因組注釋 (1)重復(fù)序列預(yù)測。通過TRF軟件預(yù)測串聯(lián)重復(fù)序列;使用RepeatModeler軟件重復(fù)序列文庫;結(jié)合repbase文庫,使用RepeatMasker軟件預(yù)測基因組中的轉(zhuǎn)座子。

(2)非編碼RNA預(yù)測。使用tRNAscan-SE(v.1.23)和RNAmmer(v.1.2)分別預(yù)測tRNA和rRNA,再使用Infernal程序(v.1.1.1)預(yù)測miRNA和snRNA。

(3)編碼基因預(yù)測和注釋。基因組的編碼基因預(yù)測使用BRAKER軟件,選擇“Fungus”選項。利用diamond blastp程序(v.2.9.0,e值=1e-5)在Uniprot、Refseq、NR、COG、GO和KEGG數(shù)據(jù)庫進行序列相似性檢索。將KEGG注釋獲得的結(jié)果,使用KOBAS程序關(guān)聯(lián)到KEGG orthology以及PATHWAY?；赨niprot數(shù)據(jù)庫中的Swiss-Port結(jié)果關(guān)聯(lián)性,獲得eggNOG數(shù)據(jù)庫的注釋信息,然后挑選COG注釋結(jié)果分類。使用Hmmscan(v.3.1;e值=1e-2)在Pfam數(shù)據(jù)庫和InterProScan數(shù)據(jù)庫進行Motif查找和結(jié)構(gòu)域預(yù)測。

2 結(jié)果與分析

2.1 Illumina測序及K-mer分析

從表1可知,通過Illumina NovaSeq平臺對N.dimidiatumH1測序共獲得6.04 Gb原始數(shù)據(jù)。去除低質(zhì)量的序列,共獲得6.03 Gb純凈數(shù)據(jù)。對純凈數(shù)據(jù)開展K-mer分析,N.dimidiatumH1基因組規(guī)模約45.06 Mb,雜合率為0.079%。

2.2 Nanopore測序和基因組組裝

Nanopore PromethION平臺測序獲得7.93 Gb原始數(shù)據(jù),序列數(shù)量為639 418條。進一步過濾低質(zhì)量序列,獲得7.67 Gb純凈數(shù)據(jù),序列數(shù)量為614 580條。平均長度為12.49 Kb,最大長度為156.40 Kb,N50長度為23.19 Kb,平均質(zhì)量值為10.25。

基于Nanopore平臺獲得的純凈數(shù)據(jù)進行基因組組裝,并利用Illumina平臺的純凈數(shù)據(jù)糾錯再去雜合,最終獲得的N.dimidiatumH1基因組大小為43.78 Mb,測序深度為175×。contig數(shù)量為14條,Gap數(shù)量為0。N50長度為3.92 Mb,平均長度為3.13 Mb,最小長度為175.69 Kb,最大長度為5.63 Mb。

2.3 組裝完整性評估

使用Pilon軟件將Illumina測序數(shù)據(jù)比對到組裝的N.dimidiatumH1基因組序列上,其比對重復(fù)率為98.86%,覆蓋率為99.94%。從表2看出,在真核基因集290個直系同源基因中,完整的單拷貝基因為286個(98.6%),完整的多拷貝基因為1個(0.3%),缺失基因僅3個(1%)。在真菌基因集758個直系同源基因中,完整的單拷貝基因為747個(98.5%),完整的單拷貝基因為2個(0.3%),基因片段僅2個(0.3%),缺失基因僅7個(0.9%)。上述結(jié)果表明,組裝獲得的N.dimidiatumH1基因組完整性較好。

2.4 基因組的結(jié)構(gòu)預(yù)測

2.4.1 編碼基因 對N.dimidiatumH1基因組中編碼基因預(yù)測表明,編碼蛋白的基因共有12 962個,其mRNA平均長度為1800 bp,絕大部分長度分布于300～9000 bp,以600～2100和6000～9000 bp基因較多(圖1)。編碼基因的序列(CDS)平均長度為1541.34 bp,絕大部分長度分布于300～6000 bp,以長度300～1800 bp的基因數(shù)量較多。外顯子數(shù)量共計41 149個,平均長度為485.52 bp;內(nèi)含子數(shù)量共計28 187個,平均長度為118.94 bp;內(nèi)含子合計長度為3.35 Mb,占基因組的7.65%。

表1 N.dimidiatum H1基因組Illumina測序和K-mer分析

表2 組裝的N.dimidiatum H1基因組BUSCO評估

2.4.2 非編碼RNA 預(yù)測結(jié)果(表3)表明,N.dimidiatumH1基因組的非編碼RNA含有核糖體RNA(rRNA)、小的調(diào)控RNA(sRNA)、核仁小RNA(snRNA)及轉(zhuǎn)運RNA(tRNA)共計205個,總長度為98 500 bp,在基因組中占比0.2250%。其中,tRNA 108個,總長度10 166 bp,在基因組中占比0.0232%;rRNA 43個,總長度為79 409 bp,在基因組中占比0.1814%;snRNA 51個,總長度為8092 bp,在基因組中占比0.0185%;sRNA 3個,總長度為833 bp,在基因組中占比0.0019%。

2.4.3 轉(zhuǎn)座子和重復(fù)序列 從表4可知,N.dimidiatumH1基因組的轉(zhuǎn)座子總量為2158個,長度為697 214 bp,占基因組的1.59%。其中,逆轉(zhuǎn)座子中以長末端重復(fù)(LTR)和長散在元件(LINE)數(shù)量較多,其中,LTR類的Gypsy類型數(shù)量841個,長度為483 597 bp,占基因組的1.10%;Copia類型數(shù)量137個,長度為13 031 bp,占基因組的0.03%;LINE類數(shù)量585個,長度為150 268 bp,占基因組的0.34%。DNA轉(zhuǎn)座子數(shù)量492個,長度為42 843 bp,占基因組的0.10%。重復(fù)序列預(yù)測:N.dimidiatumH1基因組的重復(fù)序列總數(shù)量16 960個,長度為707 932 bp,占基因組的1.62%。其中,簡單重復(fù)序列數(shù)量最多,為15 318個,長度為624 367 bp,占基因組的1.43%;低復(fù)雜度重復(fù)序列數(shù)量1614個,長度為80 898 bp,占基因組的0.18%。

圖1 N.dimidiatum H1基因組中mRNAs和CDS的長度分布

2.5 基因功能注釋

從表5可知,預(yù)測獲得的12 962個基因經(jīng)Uniprot、Pfam、Refseq、Nr、Interproscan、GO、KEGG、Pathway和COG等9個數(shù)據(jù)庫比對,共計12 703個(98.00%)基因被注釋。其中,NR數(shù)據(jù)庫注釋的基因數(shù)最多,為12 643個,占總注釋量的97.54%;其次分別為Interproscan數(shù)據(jù)庫的10 443個(80.57%)和Pfam數(shù)據(jù)庫的10 408個(80.30%);然后是Uniprot數(shù)據(jù)庫的7956個(61.38%)和GO數(shù)據(jù)庫的7882個(60.81%)。

2.5.1 NR注釋 從圖2看出,NR注釋的物種分布中,以M.phaseolinaMS6的基因最多,占總注釋量的58.36%;其余依次為N.parvumUCR-NP2(23.12%)、Diplodiacorticola(6.23%)和D.seriata

表3 N.dimidiatum H1基因組的非編碼RNA預(yù)測

表4 N.dimidiatum H1基因組的轉(zhuǎn)座子和重復(fù)序列預(yù)測

表5 N.dimidiatum H1基因組編碼基因注釋

A:物種分布;B:E值分布;a,0～1e-150;b,1e-150～1e-100;c,1e-100～1e-50;d,1e-50～1e-20;e,1e-20～1e-5;C:序列一致性分布。

(5.82%)。NR注釋的序列一致率分布中,以80%～95%的基因數(shù)量最多,占51.88%;其余依次為60%～80%(27.26%)、40%～60%(9.77%)和95%～100%(8.76%)。NR注釋的E值分布中,0～1e-150、1e-150～1e-100和1e-100～1e-50的基因數(shù)量占比分別為59.85%、15.85%和15.38%。

2.5.2 GO注釋 統(tǒng)計GO注釋(圖3)表明,在生物學(xué)過程的70個二級分類中,生物合成過程(1851個)、細胞含氮化合物代謝過程(1330個)、生物過程(1172個)、催化過程(1085個)和小分子代謝過程(962個)是數(shù)量最多的5個二級分類,占比分別為23.48%、16.87%、14.87%、13.77%和12.21%。在細胞組分的35個二級分類里,細胞核(2050個)、細胞質(zhì)(1671個)、細胞組分(1630個)、胞漿(1362個)和蛋白復(fù)合體(1357個)是數(shù)量最多的5個二級分類,占比分別為26.00%、21.20%、20.68%、17.28%和17.22%。在分子功能的41個二級分類里,離子結(jié)合(3030個)、分子功能(2054個)、 氧化還原活性(1472個)、跨膜轉(zhuǎn)運活性(744個)、DNA結(jié)合(620個)是數(shù)量最多的5個二級分類,占比分別為38.44%、26.05%、18.68%、9.44%和7.87%。

2.5.3 KEGG注釋 統(tǒng)計KEGG注釋(圖4)表明,31個二級分類中,碳水化合物代謝相關(guān)基因數(shù)量最多,為445個,占比為15.78%;其次分別為氨基酸代謝(417個)、翻譯(365個)、全局及概要圖(338個)、轉(zhuǎn)運和分解(311個),占比分別為14.79%、12.94%、11.99%和11.03%;折疊、分類和降解(296個)、信號轉(zhuǎn)導(dǎo)(253個)、脂類代謝(250個)、能量代謝(205個)及有害異物生物降解和代謝(202個),占比分別為10.50%、8.97%、8.87%、7.27%和7.16%。

3 討 論

葡萄座腔菌科屬子囊菌門(Ascomycota)座囊菌綱(Dothideomycetes)葡萄座腔菌目(Botryosphaeriales),是林木和園藝作物潰瘍病的重要真菌性病原,寄主種類多樣,地理分布廣泛[28-30]。葡萄座腔菌科中一些危害園藝作物的屬,如葡萄座腔菌屬(Botryosphaeria)[24,31-32]、殼球孢屬(Macrophomina)[23,33]、色二孢屬(Diplodia)[34]、新殼梭孢屬(Neofusicoccum)[35]和毛色二孢屬(Lasiodiplodia)[36-37]均有若干種或株系已報道基因組序列,為了解病原菌的系統(tǒng)分類、遺傳變異和致病機理奠定了基礎(chǔ)。新節(jié)格孢屬(Neoscytalidium)于2006年確立,目前僅包含N.dimidiatum、N.novaehollandiae和N.orchidacea-rum[38-39]。N.dimidiatum是新節(jié)格孢屬的典型種[38],近年來在園藝作物普遍發(fā)生[13-15]。研究表明,N.dimidiatum引起人類皮膚病或鼻竇炎等[16-17],Santana等[39]從患者皮膚分離N.dimidiatumUM880株系開展基因組測序。目前尚未有分離自園藝作物的N.dimidiatum基因組報道,本研究從火龍果潰瘍病中分離到病原新暗色柱節(jié)孢(N.dimidiatum) H1菌株并對其開展基因組測序和組裝分析。結(jié)果表明,火龍果潰瘍病病原N.dimidiatumH1基因組包含14個contig,N50長度為3.92 Mb,contig最長長度為5.63 Mb。與葡萄座腔菌科其他真菌測序策略和組裝結(jié)果比較后發(fā)現(xiàn),與采用Illumina測序策略(僅構(gòu)建0.5 Kb文庫)獲得的B.dothideaCBS115476[28]、LW030101[31]和PG45[32]、DiplodiaseriataDS831[34]、N.parvumUCR-NP2[35]和LasiodiplodiatheobromaeLA-SOL3[37]基因組相比,H1基因組的contig長度大幅增加,且數(shù)量明顯減少。與采用Illumina測序策略(構(gòu)建0.5 Kb文庫+5～6 Kb文庫)獲得的LasiodiplodiatheobromaeCSS-01s[36]和N.dimidiatumUM880[40]基因組相比,N50長度基本一致,而H1基因組的contig數(shù)量略少。與采用PacBio+Illumina或454+Illumina混合策略獲得的B.dothideasdau11-99[24]和LW-Hubei(GenBank登錄號:GCA_011064635.1)、M.phaseolina11-12[23]、Al-1[23]和MS6[33]基因組相比,N50長度也基本一致。進一步分析Illumina序列的比對率和覆蓋率,結(jié)合真核和真菌直系同源基因的BUSCO評估,均表明組裝的H1基因組完整性較高。因此,本研究獲得的N.dimidiatumH1基因組與采用PacBio+Illumina或454+Illumina策略獲得的基因組組裝水平基本一致,contig長度較長,完整性較高,為開展后續(xù)研究奠定了良好基礎(chǔ)。

圖3 N.dimidiatum H1基因組編碼基因的GO分類

組裝獲得N.dimidiatumH1的基因組大小為43.78 Mb,與B.dothideaCBS115476[28]、L.theobromaeCSS-01s[36]和LA-SOL3[37]基本一致,但略高于分離自人體的N.dimidiatumUM880[40]。H1基因組含有的編碼基因數(shù)量和基因密度均與L.theobromaeCSS-01s[36]和LA-SOL3[37]較一致,略高于N.dimidiatumUM880[40]。H1基因組含有108個tRNA,低于B.dothideaPG45[32]和N.dimidiatumUM880[40];HI基因組含有43個rRNA,介于N.dimidiatumUM880[40]和B.dothideasdau 11-99[24]之間。N.dimidiatumH1基因組中的重復(fù)序列和轉(zhuǎn)座子占比均明顯低于M.phaseolinaMS6[33],但轉(zhuǎn)座子占比顯著高于N.dimidiatumUM880[40]。在逆轉(zhuǎn)座子類型中,N.dimidiatumH1基因組以LTR類Gypsy亞族和LINE類為主要類型,與N.dimidiatumUM880、M.phaseolinaMS6及其他真菌中主要逆轉(zhuǎn)座子的類型均一致[33,40-41]。但N.dimidiatumH1基因組中的LTR類Gypsy亞族和LINE類的占比均明顯高于N.dimidiatumUM880[40]?；谏鲜龌蚪M結(jié)構(gòu)上的初步比較表明,分離自火龍果的N.dimidiatumH1菌株與分離自人類皮膚的N.dimidiatumUM880的基因組存在一定差異。

N.dimidiatumH1的基因組中含有12 962個編碼基因,共計98%的基因獲得注釋。NR注釋的物種主要為葡萄座腔菌科的M.phaseolinaMS6和N.parvumUCR-NP2;NR注釋的序列一致率主要分布在60%～95%;E值主要分布在0～1e-50,說明同源性較高,注釋結(jié)果可靠。N.dimidiatumH1的基因組GO注釋基因數(shù)量與N.dimidiatumUM880的基本一致,氧化還原酶活性、跨膜轉(zhuǎn)運蛋白活性等基因數(shù)量基本一致,但離子結(jié)合相關(guān)基因數(shù)量遠大于N.dimidiatumUM880[40]。KEGG注釋表明,N.dimidiatumH1的基因組含有碳水化合物代謝、氨基酸代謝、脂類代謝、有害異物生物降解和代謝等基因,與N.dimidiatumUM880基本一致[40]。下一步需深入開展基因功能注釋,如預(yù)測真菌致病相關(guān)的分泌蛋白及其包括的效應(yīng)蛋白、碳水化合物酶類[42-44],為篩選致病基因和開展致病機理研究奠定基礎(chǔ)。

4 結(jié) 論

利用Nanopore+Illumina混合測序策略對分離自火龍果潰瘍病的N.dimidiatumH1開展基因組測序、組裝和注釋,獲得N.dimidiatumH1的基因組大小為43.78 Mb,組裝獲得14個contig,N50長度為3.92 Mb,最長長度為5.63 Mb。N.dimidiatumH1的基因組含有12 962個編碼基因,12 703個編碼基因被注釋,與同屬于葡萄座腔菌科(Botryosphaeriaceae)的MacrophominaphaseolinaMS6和NeofusicoccumparvumUCR-NP2具有較高的同源性。該研究獲得高質(zhì)量N.dimidiatumH1的基因組序列,為了解N.dimidiatum的遺傳變異、致病機理和制定防控策略奠定了理論基礎(chǔ)。