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

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

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

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

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

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

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

1 材料與方法

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

利用PacBio SMRT測序技術對蜜蜂球囊菌的純化菌絲樣品進行測序，獲得了高質(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)嚴格質(zhì)控的數(shù)據(jù)用于后續(xù)的SNP和InDel位點鑒定與分析。

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

基于蜜蜂球囊菌菌絲的PacBio SMRT測序數(shù)據(jù)，參照Li等(2009)的方法，采用SAMtools軟件進行蜜蜂球囊菌菌絲中全長轉(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)的相關生物信息學工具對蜜蜂球囊菌菌絲中SNP位點和InDel位點基因進行GO數(shù)據(jù)庫(https:∥www.omicshare.com/tools/Home/Soft/gogsea)和KEGG數(shù)據(jù)庫(https:∥www.omicshare.com/tools/Home/Soft/pathwaygsea)比對，從而獲得相應的功能條目(term)和通路(pathway)的注釋。

2 結果

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

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

表1 基于PacBio SMRT測序數(shù)據(jù)的蜜蜂球囊菌SNP位點的詳細信息(僅展示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)。進一步分析發(fā)現(xiàn)，發(fā)生轉(zhuǎn)換和顛換的SNP位點數(shù)量分別為4 887和1 856個；雖然發(fā)生轉(zhuǎn)換的SNP數(shù)量較顛換更多，但是后者的SNP類型更為豐富，達到前者的2倍(圖1: A)?；蚪M功能元件分布的統(tǒng)計和分析結果顯示，分布在外顯子區(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)。進一步對SNP位點涉及的密碼子突變類型進行統(tǒng)計和分析，結果顯示發(fā)生同義單核苷酸突變的SNP位點數(shù)量最多，達到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ù)庫注釋結果顯示，蜜蜂球囊菌SNP位點基因可注釋到生物學進程、細胞組分和分子功能大類相關的34個功能條目，涉及細胞進程和代謝進程等13個生物學進程相關條目，細胞和細胞部分等8個細胞組分相關條目，催化活性和結合等4個分子功能相關條目(圖2: A)。KEGG數(shù)據(jù)庫注釋結果顯示，SNP位點基因可注釋到細胞進程、環(huán)境信息處理、遺傳信息處理和代謝相關的76條通路，其中注釋基因數(shù)量最多的通路是代謝通路、次生代謝物的生物合成、線粒體自噬-酵母、泛素介導的蛋白水解及內(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位點的詳細信息如表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)。進一步分析發(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位點的詳細信息(僅展示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ù)庫注釋結果顯示，上述InDel位點基因可注釋到細胞進程和代謝進程等17個生物學進程相關條目，細胞和細胞部分等12個細胞組分相關條目，催化活性和結合等10個分子功能相關條目(圖4: A)。KEGG數(shù)據(jù)庫注釋結果顯示，上述InDel位點基因還能注釋到細胞進程、環(huán)境信息處理、遺傳信息處理和代謝相關的87條通路，包括碳代謝、氨基酸的生物合成、細胞周期-酵母、真核生物中的核糖體生物發(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ā)生甲基化突變，進而自發(fā)脫去氨基而形成胸腺嘧啶。上述結果與蠶豆Viciafaba、水稻和椰心葉甲嚙小蜂Tetrastichusbrontispae等物種的研究結論(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，遠低于上述物種。Saranathan等(2017)在研究鮑氏不動桿菌Acinetobacterbaumannii時發(fā)現(xiàn)更高的SNP發(fā)生頻率有利于菌株的快速進化，并能導致菌株產(chǎn)生更強的毒力和耐藥性。這說明SNP發(fā)生頻率的高低可能會影響病原菌對宿主的侵染力及環(huán)境適應性。需要強調(diào)的是，測序深度和檢測方法等因素也可能對SNP發(fā)生頻率的檢出產(chǎn)生一定影響。

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

Chen等(2021)大規(guī)模開發(fā)了玉米Zeamays中的SNP標記，GO數(shù)據(jù)庫注釋結果顯示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ū)域，進一步的GO數(shù)據(jù)庫注釋結果顯示SNP和InDel基因可注釋到與蛋白質(zhì)結合、水解酶活性、轉(zhuǎn)錄和翻譯等相關的GO條目中。本研究中，SNP/InDel位點基因可注釋到生殖進程和發(fā)育進程等生長發(fā)育相關的功能條目；此外還可注釋到精氨酸生物合成和賴氨酸降解等15條氨基酸代謝相關通路，磷酸戊糖通路和半乳糖代謝等8條碳水化合物代謝相關通路，甘油磷脂代謝和甘油脂代謝等7條脂類代謝相關通路，以及三羧酸循環(huán)和氧化磷酸化等2條能量代謝相關通路(圖2和圖4)。以上結果表明SNP/InDel位點潛在與蜜蜂球囊菌的生長、發(fā)育、生殖及物質(zhì)和能量代謝密切相關。

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