基于PacBio測序數(shù)據(jù)的蜜蜂球囊菌SNP與InDel位點發(fā)掘及分析

蔡宗兵, 張文德, 隆 琦, 余岢駿, 孫明會, 吳 鷹, 許雅靜,劉佳美, 郭意龍, 徐細(xì)建, 陳大福,2,*, 郭 睿,2,*

(1. 福建農(nóng)林大學(xué)動物科學(xué)學(xué)院(蜂學(xué)學(xué)院), 福州 350002; 2. 福建農(nóng)林大學(xué)蜂療研究所, 福州 350002;3. 江西省養(yǎng)蜂研究所, 南昌 330000)

蜜蜂是一種具有代表性的真社會性昆蟲和不可替代的授粉昆蟲，但因群居性而易遭受病原微生物的侵襲(梁勤和陳大福, 2009)。蜜蜂球囊菌Ascosphaeraapis是一種專性侵染蜜蜂幼蟲的真菌病原，導(dǎo)致的白堊病可引起蜂群中成年工蜂數(shù)量和群勢的驟降，據(jù)統(tǒng)計，白堊病可致使蜂蜜產(chǎn)量下降5%～37%，每年給養(yǎng)蜂業(yè)造成較大損失(Zaghlouletal., 2005; Aronstein and Murray, 2010)。Shang等(2016)利用二代測序技術(shù)組裝和注釋了蜜蜂球囊菌ARSEF 7405菌株的參考基因組，為深入開展組學(xué)和分子生物學(xué)相關(guān)研究奠定了基礎(chǔ)。

單核苷酸多態(tài)性(single nucleotide polymorphism, SNP)和插入缺失(insertion-deletion, InDel)突變作為新型分子標(biāo)記被廣泛應(yīng)用于動植物和微生物的基因圖位克隆、基因型分型、遺傳圖譜繪制及遺傳多樣性等研究(Rafalski, 2002; Haraksingh and Snyder, 2013)。在球孢白僵菌Beauveriabassiana的研究中，Coates等(2002)鑒定了核小亞基中的InDel突變位點并將其應(yīng)用于昆蟲宿主中分離菌株的鑒定。在HapMap計劃中，研究者通過對人類HomosapiensSNP位點進(jìn)行基因分型繪制單倍型圖，為精確定位復(fù)雜疾病的致病基因提供了重要信息(Gonzaga-Jaureguietal., 2012)。此外，通過SNP和InDel位點與疾病或個體敏感性、耐藥性的相關(guān)性研究可指導(dǎo)遺傳育種和藥物設(shè)計，例如在結(jié)核分枝桿菌Mycobacteriumtuberculosis的研究中，Li等(2018)通過研究與利福平和異煙肼耐藥性相關(guān)的SNP和InDel位點信息檢測到大量耐藥性相關(guān)基因。

相比于基因組測序，轉(zhuǎn)錄組測序的成本較低、周期較短，為高通量鑒定SNP和InDel提供了方便快捷的工具(張倩倩等, 2019)。目前，以Illumina為代表的二代測序技術(shù)已應(yīng)用于人類(姜玥, 2013)、尼羅羅非魚Oreochromisniloticus(Yáezetal., 2020)和小麥Triticumaestivum(陳廣鳳等, 2015)等物種的SNP和InDel研究。但蜜蜂球囊菌的分子標(biāo)記相關(guān)研究較為滯后。筆者所在團(tuán)隊前期基于高質(zhì)量的Illumina測序數(shù)據(jù)組裝和注釋了蜜蜂球囊菌的參考轉(zhuǎn)錄組，并對其SSR位點進(jìn)行了大規(guī)模鑒定和分析(張曌楠等, 2017)。隨著高通量測序技術(shù)的不斷革新與發(fā)展，以PacBio單分子實時(single molecule real-time, SMRT)測序技術(shù)和納米孔(nanopore)長讀段測序技術(shù)為代表的第三代測序技術(shù)日趨成熟，已成功應(yīng)用于中華章魚Octopussinensis(Lietal., 2020)、擬南芥Arabidopsisthaliana(Cuiet

al., 2020)和東方蜜蜂微孢子蟲Nosemaceranae(陳華枝等, 2020, 2021)等動植物和微生物的全長轉(zhuǎn)錄組研究。相比于二代測序，三代測序具有超長讀長和直接讀取堿基修飾等顯著優(yōu)勢。目前，蜜蜂病原的三代組學(xué)研究匱乏。前期研究中，筆者所在團(tuán)隊利用PacBio SMRT測序技術(shù)對純化的蜜蜂球囊菌菌絲和孢子分別進(jìn)行測序，基于高質(zhì)量的長讀段測序數(shù)據(jù)構(gòu)建和注釋了蜜蜂球囊菌的首個全長轉(zhuǎn)錄組(杜宇等, 2021)。目前，對于包括蜜蜂球囊菌在內(nèi)的蜜蜂病原，尚沒有利用三代測序數(shù)據(jù)開發(fā)分子標(biāo)記的研究報道。

本研究擬利用已獲得的蜜蜂球囊菌菌絲PacBio SMRT測序數(shù)據(jù)發(fā)掘蜜蜂球囊菌的SNP和InDel位點，并分析其突變類型、基因組功能元件分布和密碼子突變類型，進(jìn)而通過功能和通路注釋探討SNP和InDel位點基因的潛在功能，以期豐富蜜蜂球囊菌的SNP和InDel位點信息，并為新型分子標(biāo)記的開發(fā)和應(yīng)用提供基礎(chǔ)。

1 材料與方法

1.1 蜜蜂球囊菌的PacBio SMRT測序數(shù)據(jù)來源

利用PacBio SMRT測序技術(shù)對蜜蜂球囊菌的純化菌絲樣品進(jìn)行測序，獲得了高質(zhì)量的全長轉(zhuǎn)錄組數(shù)據(jù)(Chenetal., 2020)，共測得13 302 489條subreads(約23.97 Gb)，平均讀長和居中長度(N50)分別為1 802和3 077 bp；經(jīng)多輪校正共得到174 095條高質(zhì)量全長轉(zhuǎn)錄本，總堿基數(shù)為474 928 820，平均讀長和N50分別為2 728和3 543 bp(Chenetal., 2020)。經(jīng)嚴(yán)格質(zhì)控的數(shù)據(jù)用于后續(xù)的SNP和InDel位點鑒定與分析。

1.2 SNP位點和InDel位點的鑒定與分析

基于蜜蜂球囊菌菌絲的PacBio SMRT測序數(shù)據(jù)，參照Li等(2009)的方法，采用SAMtools軟件進(jìn)行蜜蜂球囊菌菌絲中全長轉(zhuǎn)錄本識別。參照Wang等(2010)的方法，利用ANNOVAR軟件將識別的全長轉(zhuǎn)錄本比對到蜜蜂球囊菌參考基因組(assembly AAP 1.0)以檢測SNP位點和InDel位點，檢測內(nèi)容包括：(1)SNP位點的突變的類型和數(shù)量，發(fā)生轉(zhuǎn)換和顛換的SNP位點的數(shù)量與比率；(2)SNP位點和InDel位點在參考基因組7種功能元件(外顯子區(qū)、內(nèi)含子區(qū)、基因上游區(qū)、基因下游區(qū)、基因上下游重疊區(qū)、基因間隔區(qū)、剪接區(qū))上的分布數(shù)量和占比；(3)SNP位點和InDel位點中不同類型密碼子突變的數(shù)量和占比。

1.3 SNP位點和InDel位點基因的數(shù)據(jù)庫注釋

使用基迪奧生物云平臺(https:∥www.omicshare.com/tools/home/index/index)的相關(guān)生物信息學(xué)工具對蜜蜂球囊菌菌絲中SNP位點和InDel位點基因進(jìn)行GO數(shù)據(jù)庫(https:∥www.omicshare.com/tools/Home/Soft/gogsea)和KEGG數(shù)據(jù)庫(https:∥www.omicshare.com/tools/Home/Soft/pathwaygsea)比對，從而獲得相應(yīng)的功能條目(term)和通路(pathway)的注釋。

2 結(jié)果

2.1 蜜蜂球囊菌SNP位點的鑒定與分析

共鑒定到蜜蜂球囊菌的6 743個SNP位點，包括6 091個純合位點和652個雜合位點。SNP位點的詳細(xì)信息如表1所示。

表1 基于PacBio SMRT測序數(shù)據(jù)的蜜蜂球囊菌SNP位點的詳細(xì)信息(僅展示10個)Table 1 Detailed information of SNP sites in Ascosphaera apis (only 10 presented) based on PacBio SMRT sequencing data

上述蜜蜂球囊菌的SNP位點的堿基突變類型共分為12種，包括C/T, G/A, A/G, T/C, G/T, A/T, C/A, T/A, T/G, A/C, C/G和G/C(圖1: A)；其中最豐富的突變類型為C/T型(1 296個)，最少的突變類型為G/C(173個)(圖1: A)。進(jìn)一步分析發(fā)現(xiàn)，發(fā)生轉(zhuǎn)換和顛換的SNP位點數(shù)量分別為4 887和1 856個；雖然發(fā)生轉(zhuǎn)換的SNP數(shù)量較顛換更多，但是后者的SNP類型更為豐富，達(dá)到前者的2倍(圖1: A)。基因組功能元件分布的統(tǒng)計和分析結(jié)果顯示，分布在外顯子區(qū)的SNP位點最多，共計3 860個(占57.24%)(圖1: B)；其次是分布在基因間隔區(qū)(1 117個，占16.57%)、基因下游區(qū)(781個，占11.58%)、基因上游區(qū)(614個，占9.11%)、基因上下游重疊區(qū)(207個，占3.07%)和內(nèi)含子區(qū)(160個，占2.37%)(圖1: B)；分布SNP位點數(shù)量最少的是剪接區(qū)，僅有4個(占0.06%)(圖1: B)。進(jìn)一步對SNP位點涉及的密碼子突變類型進(jìn)行統(tǒng)計和分析，結(jié)果顯示發(fā)生同義單核苷酸突變的SNP位點數(shù)量最多，達(dá)到2 892個(占74.92%)，其次是非同義單核苷酸突變(950個，占24.61%)和終止子增加(17個，占0.44%)；發(fā)生終止子減少的SNP數(shù)量最少，僅為1個(占0.03%)(圖1: C)。

圖1 基于PacBio SMRT測序數(shù)據(jù)的蜜蜂球囊菌SNP的突變類型(A)及基因組位置分布(B)和功能類型(C)Fig. 1 Mutation type (A), location distribution in genome (B) and functional type (C) of SNPsin Ascosphaera apis based on PacBio SMRT sequencing dataSNV: 單核苷酸突變Single nucleotide mutation.

2.2 蜜蜂球囊菌SNP位點所在基因的數(shù)據(jù)庫注釋

GO數(shù)據(jù)庫注釋結(jié)果顯示，蜜蜂球囊菌SNP位點基因可注釋到生物學(xué)進(jìn)程、細(xì)胞組分和分子功能大類相關(guān)的34個功能條目，涉及細(xì)胞進(jìn)程和代謝進(jìn)程等13個生物學(xué)進(jìn)程相關(guān)條目，細(xì)胞和細(xì)胞部分等8個細(xì)胞組分相關(guān)條目，催化活性和結(jié)合等4個分子功能相關(guān)條目(圖2: A)。KEGG數(shù)據(jù)庫注釋結(jié)果顯示，SNP位點基因可注釋到細(xì)胞進(jìn)程、環(huán)境信息處理、遺傳信息處理和代謝相關(guān)的76條通路，其中注釋基因數(shù)量最多的通路是代謝通路、次生代謝物的生物合成、線粒體自噬-酵母、泛素介導(dǎo)的蛋白水解及內(nèi)吞作用等(圖2: B)。

圖2 基于PacBio SMRT測序數(shù)據(jù)的蜜蜂球囊菌SNP位點基因的GO(A)和KEGG(B)數(shù)據(jù)庫注釋Fig. 2 GO (A) and KEGG (B) database annotation of genes with SNP site in Ascosphaera apisbased on PacBio SMRT sequencing data

2.3 蜜蜂球囊菌的InDel位點的鑒定與分析

共鑒定到蜜蜂球囊菌的597個InDel位點，包括349個純合位點和248個雜合位點。InDel位點的詳細(xì)信息如表2所示。

此外，InDel位點分布最多的為基因下游區(qū)(182個，占30.49%)，其次是基因上游區(qū)(146個，占24.45%)、外顯子區(qū)(109個，占18.26%)、基因間隔區(qū)(63個，占10.55%)、基因上下游重疊區(qū)(60個，占10.05%)和內(nèi)含子區(qū)(37個，占6.20%)(圖3: A)。進(jìn)一步分析發(fā)現(xiàn)，最豐富密碼子突變類型為非移碼插入(37個，占33.94%)，其次是非移碼缺失(33個，占30.28%)、移碼缺失(22個，占20.18%)和移碼插入(16個，占14.68%)，最少的終止子增加僅有1個，占0.92%(圖3: B)。

表2 基于PacBio SMRT測序數(shù)據(jù)的蜜蜂球囊菌InDel位點的詳細(xì)信息(僅展示10個)Table 2 Detailed information of InDel sites in Ascosphaera apis (only 10 presented)based on PacBio SMRT sequencing data

圖3 基于PacBio SMRT測序數(shù)據(jù)的蜜蜂球囊菌中InDel位點所在基因組位置分布(A)和功能類型(B)Fig. 3 Location distribution in genome (A) and functional type (B) of genes with InDel sites in Ascosphaera apisbased on PacBio SMRT sequencing data

2.4 蜜蜂球囊菌的InDel位點基因的數(shù)據(jù)庫注釋

GO數(shù)據(jù)庫注釋結(jié)果顯示，上述InDel位點基因可注釋到細(xì)胞進(jìn)程和代謝進(jìn)程等17個生物學(xué)進(jìn)程相關(guān)條目，細(xì)胞和細(xì)胞部分等12個細(xì)胞組分相關(guān)條目，催化活性和結(jié)合等10個分子功能相關(guān)條目(圖4: A)。KEGG數(shù)據(jù)庫注釋結(jié)果顯示，上述InDel位點基因還能注釋到細(xì)胞進(jìn)程、環(huán)境信息處理、遺傳信息處理和代謝相關(guān)的87條通路，包括碳代謝、氨基酸的生物合成、細(xì)胞周期-酵母、真核生物中的核糖體生物發(fā)生和RNA轉(zhuǎn)運等(圖4: B)。

圖4 基于PacBio SMRT測序數(shù)據(jù)的蜜蜂球囊菌InDel位點所在基因的GO(A)和KEGG(B)數(shù)據(jù)庫注釋Fig. 4 GO (A) and KEGG (B) database annotation of genes with InDel site in Ascosphaera apisbased on PacBio SMRT sequencing data

3 討論

本研究利用前期已獲得的蜜蜂球囊菌菌絲的PacBio SMRT測序數(shù)據(jù)，共鑒定到6 743個SNP位點，包含12種堿基突變類型，其中最豐富的類型為C/T型；發(fā)生轉(zhuǎn)換和顛換的SNP分別為4 887和1 856個；SNP位點在蜜蜂球囊菌參考基因組的7種功能元件上均有分布，分布數(shù)量最多的元件為外顯子；SNP位點涉及8種密碼子突變類型，以同義單核苷酸突變最為豐富(圖1)。此外還鑒定到597個InDel位點，在基因下游區(qū)的分布數(shù)量最多；InDel位點涉及5種密碼子突變類型，其中最豐富的類型為非移碼插入突變(圖3)。

本研究檢測到的SNP顛換類型共有8種，為轉(zhuǎn)換類型的2倍，推測這是由于基因組中不同堿基的突變頻率不同所致。另外，在SNP的12種堿基突變類型中，發(fā)生頻率最高的為C/T突變，這可能是由于CpG二核苷酸上的胞嘧啶殘基極易發(fā)生甲基化突變，進(jìn)而自發(fā)脫去氨基而形成胸腺嘧啶。上述結(jié)果與蠶豆Viciafaba、水稻和椰心葉甲嚙小蜂Tetrastichusbrontispae等物種的研究結(jié)論(Temnykhetal., 2001; Ocaaetal., 2015; 張潔慧等, 2021)相似。不同物種的SNP發(fā)生頻率存在差異，例如在小麥條銹菌Pucciniastriiformis基因組中為1個/0.67 kb(Kiranetal., 2017)，在向日葵銹菌Pucciniahelianthi基因組中為1個/2.8 kb(王妍等, 2018)。本研究中，蜜蜂球囊菌的SNP發(fā)生頻率約為1個/70 kb，遠(yuǎn)低于上述物種。Saranathan等(2017)在研究鮑氏不動桿菌Acinetobacterbaumannii時發(fā)現(xiàn)更高的SNP發(fā)生頻率有利于菌株的快速進(jìn)化，并能導(dǎo)致菌株產(chǎn)生更強(qiáng)的毒力和耐藥性。這說明SNP發(fā)生頻率的高低可能會影響病原菌對宿主的侵染力及環(huán)境適應(yīng)性。需要強(qiáng)調(diào)的是，測序深度和檢測方法等因素也可能對SNP發(fā)生頻率的檢出產(chǎn)生一定影響。

外顯子與氨基酸翻譯過程直接相關(guān)，其堿基的改變極易引起編碼氨基酸的改變，最終影響蛋白質(zhì)的功能。本研究中，蜜蜂球囊菌的SNP位點主要分布在基因組的外顯子區(qū)(3 860個，占57.2%)，而在其他區(qū)域的分布相對較少；進(jìn)一步分析發(fā)現(xiàn)對于蜜蜂球囊菌的SNP位點，同義單核苷酸突變的數(shù)量最多(2 892個，占74.9%)(圖1: C)，鑒于同義單核苷酸突變并不會造成氨基酸序列的變化，上述結(jié)果表明蜜蜂球囊菌在長期的進(jìn)化過程中需要保持基因組編碼蛋白質(zhì)的穩(wěn)定，這對蜜蜂球囊菌的存活與進(jìn)化至關(guān)重要。InDel可能發(fā)生移碼突變，致使mRNA在翻譯時遇到錯誤的終止密碼子。Zhang ZK等(2020)將17株來自不同宿主的球孢白僵菌Beauveriabassiana與參考基因組比較，共鑒定出10 098個非同義突變基因，功能注釋結(jié)果表明其中的大部分都涉及毒力蛋白的生物學(xué)功能。Tambong等(2014)以采自加拿大安大略省南部查爾斯頓湖附近林地的黏質(zhì)沙雷氏菌Serratiamarcescens為研究材料，通過與黏質(zhì)沙雷氏菌E-15參考菌株的金屬蛋白酶基因比較后鑒定到位于該基因上的72個SNP位點和3個InDel位點，進(jìn)而鑒定得到8個非同義核苷酸突變，對蛋白質(zhì)變異影響的分析表明，由蛋白質(zhì)結(jié)構(gòu)中的非同義核苷酸突變產(chǎn)生的新天冬氨酸殘基可能對其生物學(xué)功能產(chǎn)生最顯著的影響。鑒于此，SNP和InDel位點對蜜蜂球囊菌的生長發(fā)育具有潛在影響。

Chen等(2021)大規(guī)模開發(fā)了玉米Zeamays中的SNP標(biāo)記，GO數(shù)據(jù)庫注釋結(jié)果顯示SNP基因可注釋到大部分初級和次級代謝通路。Zhang Y等(2020)在全基因組范圍鑒定了鳳頭鴨的SNP和InDel位點，使用GO和KEGG數(shù)據(jù)庫注釋后發(fā)現(xiàn)SNP和InDel位點基因可富集到包括骨化、軟骨發(fā)育、大分子生物合成過程等條目和Hedgehog信號通路、甘油脂代謝、磷脂酰肌醇信號系統(tǒng)等通路。Calarco等(2018)曾檢測了犬新孢子Neosporacaninum中的SNP和InDel位點，鑒定并驗證了核基因組編碼區(qū)的SNP熱點區(qū)域，進(jìn)一步的GO數(shù)據(jù)庫注釋結(jié)果顯示SNP和InDel基因可注釋到與蛋白質(zhì)結(jié)合、水解酶活性、轉(zhuǎn)錄和翻譯等相關(guān)的GO條目中。本研究中，SNP/InDel位點基因可注釋到生殖進(jìn)程和發(fā)育進(jìn)程等生長發(fā)育相關(guān)的功能條目；此外還可注釋到精氨酸生物合成和賴氨酸降解等15條氨基酸代謝相關(guān)通路，磷酸戊糖通路和半乳糖代謝等8條碳水化合物代謝相關(guān)通路，甘油磷脂代謝和甘油脂代謝等7條脂類代謝相關(guān)通路，以及三羧酸循環(huán)和氧化磷酸化等2條能量代謝相關(guān)通路(圖2和圖4)。以上結(jié)果表明SNP/InDel位點潛在與蜜蜂球囊菌的生長、發(fā)育、生殖及物質(zhì)和能量代謝密切相關(guān)。

絲裂原活化蛋白激酶(mitogen-activated protein kinase, MAPK)級聯(lián)反應(yīng)在真菌的生長發(fā)育、致病性、繁殖等方面發(fā)揮關(guān)鍵的調(diào)控作用(Igbariaetal., 2008; 張楠等, 2017)。Di Pietro等(2003)靶向敲除了尖孢鐮刀菌FusariumoxysporumMAPK基因簇中的Fmk1基因后，發(fā)現(xiàn)產(chǎn)生的突變體可以在人工培養(yǎng)基上生長而不能侵入番茄Lycopersiconesculentum根部，導(dǎo)致其喪失了對番茄的致病能力。張楠等(2017)敲除了膠孢炭疽菌ColletotrichumgloeosporioidesMAPK信號通路上的CgSho1基因，通過與野生型相比發(fā)現(xiàn)其營養(yǎng)生長緩慢，菌絲稀疏且產(chǎn)孢量下降，致病力也明顯減弱。Jin等(2014)通過基因敲除技術(shù)對蝗綠僵菌Metarhiziumacridum中MAPK信號通路相關(guān)基因MaMk1進(jìn)行靶向敲除后，發(fā)現(xiàn)其不能穿透蝗蟲的角質(zhì)層，完全喪失了致病性。本研究中，共有4個SNP位點基因和6個InDel位點基因可注釋到MAPK信號通路，說明相關(guān)SNP/InDel位點與蜜蜂球囊菌的MAPK信號通路具有潛在關(guān)聯(lián)，未來可參考上述前人已報道的方法嘗試對注釋到MAPK信號通路的SNP和InDel位點基因進(jìn)行敲除以揭示其功能，進(jìn)而為白堊病的治療提供分子靶點。