基于PacBio測(cè)序數(shù)據(jù)的蜜蜂球囊菌轉(zhuǎn)錄因子、融合基因及RNA編輯事件的鑒定

許雅靜, 吳 鷹, 余岢駿, 孫明會(huì), 劉佳美, 郭意龍, 徐細(xì)建,鮑佳益, 康育欣, 陳大福,2, 郭 睿,2,*, 付中民,2,*

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

蜜蜂球囊菌Ascosphaeraapis專性侵染蜜蜂幼蟲而引發(fā)白堊病，該病可導(dǎo)致蜂群群勢(shì)和生產(chǎn)力的急劇下降，給養(yǎng)蜂業(yè)造成較大損失(郭睿等, 2017a)。此前，由于參考基因組的長期缺失，蜜蜂球囊菌的組學(xué)和分子生物學(xué)研究舉步維艱。Shang等(2016)利用二代測(cè)序技術(shù)測(cè)序和組裝了蜜蜂球囊菌的參考基因組，并公布了完整的基因序列和功能注釋，為其分子生物學(xué)和組學(xué)研究奠定了基礎(chǔ)。筆者所在團(tuán)隊(duì)前期利用Illumina測(cè)序技術(shù)對(duì)蜜蜂球囊菌開展了較為系統(tǒng)的轉(zhuǎn)錄組研究(陳大福等, 2017; 郭睿等, 2017b; 張曌楠等, 2017)。

轉(zhuǎn)錄因子(transcription factor, TF)是指能直接或間接與真核生物基因的啟動(dòng)子區(qū)域中的順式作用元件發(fā)生特異性相互作用，從而抑制或激活基因表達(dá)的一種DNA結(jié)合蛋白。作為轉(zhuǎn)錄水平調(diào)控基因表達(dá)的主要方式，TF可根據(jù)細(xì)胞類型、發(fā)育階段和疾病狀態(tài)來調(diào)節(jié)轉(zhuǎn)錄起始(Fengetal., 2020; Kimetal., 2020)。同時(shí)轉(zhuǎn)錄因子是由多個(gè)基因家族編碼，大大增加了它們的數(shù)量并使轉(zhuǎn)錄調(diào)控機(jī)制復(fù)雜化，Lamber等(2018)根據(jù)DNA結(jié)合區(qū)的特征將TF分為71個(gè)家族。所謂融合基因(fusion gene)，是指將兩個(gè)或多個(gè)基因的編碼區(qū)首尾相連并置于同一套調(diào)控序列控制之下，從而構(gòu)成的嵌合基因?；蛉诤峡赏ㄟ^轉(zhuǎn)錄介導(dǎo)來實(shí)現(xiàn)，例如相鄰基因的反式或順式剪接、讀取，其中順式剪接發(fā)生在相同的前體mRNA分子內(nèi)，而反式剪接可使兩個(gè)單獨(dú)的前體mRNA分子形成嵌合的非共線RNA，該RNA可以編碼新蛋白質(zhì)或充當(dāng)非編碼或調(diào)控RNA(Leietal., 2016)，也被稱為轉(zhuǎn)錄衍生基因融合；此外還能通過各種結(jié)構(gòu)重排來實(shí)現(xiàn)，例如基因復(fù)制或染色體易位、缺失、插入和倒位，也被稱為DNA介導(dǎo)的基因融合(McCartneyetal., 2019)。目前，融合基因的相關(guān)研究多集中在人類重大疾病方面(Dupainetal., 2017; Yoonetal., 2019)，也有少量研究涉及細(xì)菌(Baietal., 2006)和病毒(Valencia-Herreraetal., 2019)。但真菌的融合基因研究未見報(bào)道。RNA編輯(RNA editing)指在RNA分子中特定核苷酸序列發(fā)生變化的過程，包括核苷酸的插入、缺失和替換等，這種轉(zhuǎn)錄后修飾不僅可以增加基因產(chǎn)物的多樣性，還直接或間接地參與了基因表達(dá)調(diào)控，在諸多生命活動(dòng)中發(fā)揮重要作用(Leongetal., 2019)。人們?cè)阱F蟲的線粒體中首次鑒定到RNA編輯事件(Benneetal., 1986)，隨后在人類腸道的載脂蛋白基因轉(zhuǎn)錄本、小麥Triticumaestivum的線粒體中也相繼發(fā)現(xiàn)了RNA編輯現(xiàn)象(Powelletal., 1987; Covello and Gray, 1989)。隨著高通量測(cè)序技術(shù)的日趨成熟，為RNA編輯事件的大規(guī)模挖掘提供了技術(shù)手段，研究表明RNA編輯現(xiàn)象在動(dòng)物、植物和微生物中廣泛存在(Bar-Yaacovetal., 2018; Porathetal., 2019; Tangetal., 2019)。

PacBio單分子實(shí)時(shí)(single molecule real-time, SMRT)測(cè)序技術(shù)是近年來興起的新一代測(cè)序技術(shù)，與一代和二代測(cè)序技術(shù)相比，PacBio SMRT測(cè)序技術(shù)具有超長讀長及可直接檢測(cè)堿基上的化學(xué)修飾等顯著優(yōu)勢(shì)，已成功應(yīng)用于跳鐮猛蟻Harpegnathossaltator(Shieldsetal., 2018)、衛(wèi)氏并殖吸蟲Paragonimuswestermani(Oeyetal., 2019)和小麥(Moolhuijzenetal., 2020)等動(dòng)植物的相關(guān)研究。近期，筆者所在團(tuán)隊(duì)利用PacBio SMRT測(cè)序技術(shù)對(duì)蜜蜂球囊菌的純化菌絲樣品進(jìn)行測(cè)序，構(gòu)建和注釋了蜜蜂球囊菌的首個(gè)高質(zhì)量全長轉(zhuǎn)錄組，解析了基因的可變剪接與可變腺苷酸化，并鑒定和分析了長鏈非編碼RNA(long non-coding RNA, lncRNA)和TF(Chenetal., 2020)。目前，蜜蜂球囊菌的TF相關(guān)信息十分有限，融合基因和RNA編輯的相關(guān)研究仍然缺失。本研究擬利用已獲得的蜜蜂球囊菌菌絲和孢子的PacBio SMRT測(cè)序數(shù)據(jù)對(duì)TF、融合基因和RNA編輯事件進(jìn)行系統(tǒng)鑒定和分析，以期豐富蜜蜂球囊菌的相關(guān)信息，并為深入探究它們的功能提供依據(jù)和基礎(chǔ)。

1 材料與方法

1.1 蜜蜂球囊菌菌絲和孢子的PacBio SMRT測(cè)序數(shù)據(jù)來源

前期研究中，筆者所在團(tuán)隊(duì)已在實(shí)驗(yàn)室條件下獲得蜜蜂球囊菌的純培養(yǎng)，并制備了純化的菌絲樣品(AaM)和孢子樣品(AaS)；此外，已利用PacBio SMRT測(cè)序技術(shù)對(duì)AaM和AaS進(jìn)行測(cè)序，AaM共測(cè)得13 302 489條subreads，平均長度和居中長度(N50)分別為1 802 bp和3 077 bp；檢測(cè)到464 043條環(huán)形一致性序列(circular consensus sequence, CCS)，平均長度為2 970 bp；經(jīng)嚴(yán)格校正最終得到174 095條全長轉(zhuǎn)錄本，平均長度和N50分別為2 728 bp和3 543 bp(Chenetal., 2020)。AaS共測(cè)得9 911 345條subreads，平均長度和N50分別為1 742 bp和2 731 bp；檢測(cè)到315 135條CCS，平均長度為2 733 bp；經(jīng)校正最終得到103 845條全長轉(zhuǎn)錄本，平均長度和N50分別為2 502 bp和3 262 bp(未發(fā)表數(shù)據(jù))。高質(zhì)量的PacBio SMRT測(cè)序可為本研究中TF、融合基因和RNA編輯事件的鑒定與分析提供可靠的數(shù)據(jù)基礎(chǔ)。

1.2 轉(zhuǎn)錄因子的鑒定及分析

利用BLASTx工具(http:∥www.ncbi.nlm.nih.gov/BLAST/)將AaM和AaS中鑒定到的全長轉(zhuǎn)錄本序列比對(duì)到Nr(http:∥www.ncbi.nlm.nih.gov), Swiss-Prot(http:∥us.expasy.org/sprot/)和KEGG(http:∥www.genome.jp/kegg)數(shù)據(jù)庫，獲得一致性最高的蛋白序列。再利用hmmscan軟件將上述蛋白序列比對(duì)到Plant TFdb數(shù)據(jù)庫(http:∥planttfdb.cbi.pku.edu.cn/)，從而獲得TF的分類及注釋信息?；駻為TOFU軟件分析過程中對(duì)組成該融合基因中一個(gè)基因的編號(hào)，基因B為TOFU軟件分析過程中對(duì)組成該融合基因中另一個(gè)基因的編號(hào)。

1.3 融合基因的鑒定及分析

采用TOFU分析套件中的fusion_finder.py程序進(jìn)行融合基因的預(yù)測(cè)，參數(shù)設(shè)置為：dun-merge-5-shorter: turned off, min_locus_coverage: 0.05, min_locus_coverage_bp: 1, min_total_coverage: 0.99, min_dist_between_loci: 10 000。然后根據(jù)預(yù)測(cè)結(jié)果分析融合基因的序列和位置信息。

1.4 RNA編輯事件的預(yù)測(cè)及分析

首先使用SAMtools工具(Lietal., 2009)預(yù)測(cè)AaM和AaS中的RNA編輯事件。然后利用ANNOVAR軟件(Wangetal., 2010)對(duì)RNA編輯事件進(jìn)行注釋，注釋結(jié)果包括同義單核苷酸突變(synonymous single nucleotide mutation)、非同義單核苷酸突變(nonsynonymous single nucleotide mutation)和終止子獲得(stop-gain)3種功能類型。最后利用基迪奧在線云平臺(tái)(www.omicshare.com)的相關(guān)軟件對(duì)RNA編輯位點(diǎn)基因進(jìn)行功能和通路注釋。

2 結(jié)果

2.1 蜜蜂球囊菌孢子中的轉(zhuǎn)錄因子

AaM中的轉(zhuǎn)錄因子的相關(guān)分析結(jié)果前期已另文發(fā)表(Chenetal., 2020)。在AaS中共鑒定到來源于17個(gè)TF家族的213個(gè)TF；其中來源于C2H2家族的TF數(shù)量最多，達(dá)到72個(gè)；其次為來源于bZIP和bHLH家族，TF數(shù)量均為25個(gè)；來源于TALE和Trilhelix家族的TF數(shù)量最少，均僅為1個(gè)(圖1)。

圖1 基于PacBio測(cè)序數(shù)據(jù)蜜蜂球囊菌孢子中鑒定到的轉(zhuǎn)錄因子

2.2 蜜蜂球囊菌菌絲和孢子中的融合基因

在AaM中共鑒定到921個(gè)融合基因，基因A來源于正鏈和負(fù)鏈的數(shù)量和占比分別為315個(gè)(34.2%)和606個(gè)(65.8%)，基因B來源于正鏈和負(fù)鏈的數(shù)量和占比分別為369個(gè)(40.1%)和552個(gè)(59.9%)。此外，基因A分布于39條scaffold，其中分布基因A數(shù)量最多的scaffold為AZGZ01000008.1(320個(gè), 占34.7%)，其次為AZGZ01000001.1(108個(gè), 占11.7%)和AZGZ01000002.1(84個(gè), 占9.1%)；分布基因A數(shù)量最少的scaffold為AZGZ01000038.1(1個(gè), 占0.1%)；基因B分布于57條scaffold，其中分布基因B數(shù)量最多的scaffold為AZGZ01000008.1(202個(gè), 占21.9%)，其次為scaffold AZGZ01000014.1(37個(gè), 占4.0%)和scaffold AZGZ01000011.1(36個(gè), 占3.9%)，分布基因B數(shù)量最少的scaffold為AZGZ01000049.1(1個(gè), 占0.1%)。在AaS中共鑒定到510個(gè)融合基因，基因A來源于正鏈和負(fù)鏈的數(shù)量和占比分別為160個(gè)(31.4%)和350個(gè)(68.6%)，基因B來源于正鏈和負(fù)鏈的數(shù)量和占比分別為214個(gè)(42.0%)和296個(gè)(58.0%)。此外，基因A分布于33條scaffold，其中分布基因A數(shù)量最多的scaffold為AZGZ01000008.1(200個(gè), 占39.2%)，其次為scaffold AZGZ01000001.1(51個(gè), 占10.0%)和scaffold AZGZ01000002.1(43個(gè), 占8.4%)，分布基因A數(shù)量最少的scaffold為AZGZ01000034.1(1個(gè), 占0.2%)；基因B分布位于54條scaffold，其中分布基因B數(shù)量最多的scaffold為AZGZ01000008.1(118個(gè), 占23.1%)，其次為scaffold AZGZ01000029.1(27個(gè), 占5.3%)和AZGZ01000013.1(15個(gè), 占2.9%)，分布基因B數(shù)量最少的scaffold為AZGZ01000054.1(1個(gè), 占0.2%)。進(jìn)一步分析發(fā)現(xiàn)，AaM和AaS共有的融合基因?yàn)?10個(gè)，AaM特有的融合基因?yàn)?11個(gè)，而AaS無特有的融合基因。

2.3 蜜蜂球囊菌菌絲和孢子中的RNA編輯事件的鑒定與比較

在蜜蜂球囊菌中共鑒定到738次RNA編輯事件，涉及3種功能類型(圖2: A)。在AaM中鑒定到的RNA編輯事件為547次，其中同義單核苷酸突變的數(shù)量最多，達(dá)到360次(占65.8%)，其次為非同義單核苷酸突變(171次, 占31.3%)和終止子獲得(16次, 占2.9%)；在AaS中鑒定到的RNA編輯事件為191次，其中非同義單核苷酸突變的數(shù)量最多，為119次(占62.3%)，同義單核苷酸突變(63次, 占33%)和終止子獲得(9次, 占4.7%)次之(圖2: B)。

圖2 基于PacBio測(cè)序數(shù)據(jù)蜜蜂球囊菌菌絲和孢子中RNA編輯事件的功能類型統(tǒng)計(jì)

在蜜蜂球囊菌中共鑒定到738次RNA編輯事件，涉及12種堿基替換類型(圖3: A)。在AaM中鑒定到12種堿基替換類型，其中發(fā)生C->T的RNA編輯事件數(shù)量最多(158次，占28.8%)，其次為T->C(89次, 占16.2%)和C->A(58次, 占10.6%)。在AaS中鑒定到9種堿基替換類型，其中發(fā)生C->T和G->T的RNA編輯事件數(shù)量最多，均為42次(占21.9%)，其次為C->A(39次, 占20.4%)，A->T(27次, 占14.1%)和T->C(16次, 占8.3%)(圖3: B)。

圖3 基于PacBio測(cè)序數(shù)據(jù)蜜蜂球囊菌菌絲和孢子中RNA編輯事件堿基替換類型統(tǒng)計(jì)

GO數(shù)據(jù)庫注釋結(jié)果顯示，AaM中RNA編輯位點(diǎn)基因可注釋到生物學(xué)進(jìn)程大類相關(guān)的8個(gè)功能條目，包括代謝進(jìn)程(10)、細(xì)胞進(jìn)程(8)和單細(xì)胞進(jìn)程(6)等；分子功能大類相關(guān)的4個(gè)功能條目，包括結(jié)合(4)、催化活性(4)和結(jié)構(gòu)分子活性(1)等；細(xì)胞組分大類相關(guān)的7個(gè)功能條目，包括細(xì)胞(5)、細(xì)胞部分(5)和細(xì)胞器(4)等(圖5)。AaS中RNA編輯位點(diǎn)基因可注釋到細(xì)胞進(jìn)程(12)和代謝進(jìn)程(11)等11個(gè)生物學(xué)進(jìn)程大類相關(guān)功能條目，催化活性(13)和結(jié)合(8)等5個(gè)分子功能大類相關(guān)功能條目，細(xì)胞(10)和細(xì)胞部分(10)等8個(gè)細(xì)胞組分大類相關(guān)功能條目(圖4)。括號(hào)內(nèi)的數(shù)字代表注釋在該條目的基因數(shù)量。

圖4 基于PacBio測(cè)序數(shù)據(jù)蜜蜂球囊菌菌絲(AaM)和孢子(AaS)中RNA編輯位點(diǎn)基因的GO數(shù)據(jù)庫注釋

KEGG數(shù)據(jù)庫注釋結(jié)果顯示，AaM中RNA編輯位點(diǎn)基因可注釋到11條通路，包括內(nèi)吞作用(2)、代謝通路(2)、磷脂酰肌醇信號(hào)系統(tǒng)(1)、核糖體(1)和次級(jí)代謝產(chǎn)物的生物合成(1)等(圖5: A)；AaS中RNA編輯位點(diǎn)基因可注釋到20條通路，包括代謝通路(3)、次級(jí)代謝產(chǎn)物的生物合成(2)、內(nèi)吞作用(2)、丙酮酸代謝(1)和RNA轉(zhuǎn)運(yùn)(1)等(圖5: B)。括號(hào)內(nèi)的數(shù)字代表注釋在該通路的基因數(shù)量。

圖5 基于PacBio測(cè)序數(shù)據(jù)蜜蜂球囊菌菌絲(A)和孢子(B)中RNA編輯位點(diǎn)基因KEGG數(shù)據(jù)庫注釋

3 討論

本研究中，在蜜蜂球囊菌孢子中鑒定到來源于17個(gè)TF家族的213個(gè)成員，總數(shù)少于菌絲中的TF；其中包含TF最多的家族也為C2H2(72)，與菌絲中的情況一致；但包含TF最少的家族為TALE(1)和Trihelix(1)(圖1)，與菌絲中的情況存在差異；另外，菌絲中C2H2家族包含的TF數(shù)量多于孢子(Chenetal., 2020)。上述結(jié)果說明C2H2家族作為蜜蜂球囊菌最大的轉(zhuǎn)錄因子家族，其成員具有潛在的重要性。在動(dòng)植物中，C2H2家族也是最重要的TF家族之一，Wu等(2019)研究闡明BmBlimp-1是家蠶Bombyxmori翅發(fā)育的重要調(diào)節(jié)因子。有研究表明該家族參與植物細(xì)胞生理狀態(tài)變化期間細(xì)胞活動(dòng)的調(diào)節(jié)(Maternaetal., 2006; Kametal., 2008)。棉花黃萎病菌VerticilliumdahliaeVdMsn2與酵母C2H2家族的轉(zhuǎn)錄因子Msn2具有同源性，有研究顯示VdMsn2的缺失能夠?qū)е旅藁S萎病菌菌絲生長緩慢，隔膜和菌絲分枝增加(Tianetal., 2017)。因此，推測(cè)轉(zhuǎn)錄因子C2H2家族可能與蜜蜂球囊菌菌絲和孢子的生長發(fā)育和細(xì)胞活動(dòng)等生物學(xué)過程密切相關(guān)。

基因的融合現(xiàn)象在人類疾病中廣泛存在。此外，融合基因也是許多癌癥的重要驅(qū)動(dòng)因素，因而可作為抗癌治療中潛在的診斷標(biāo)記和治療靶點(diǎn)(Yakushinaetal., 2018)。研究表明BRAF和KIAA1549基因的融合存在于80%的細(xì)胞性星形細(xì)胞瘤(Jeuken and Wesseling, 2010)；融合基因TRIM52-RACK1可促進(jìn)口腔鱗狀細(xì)胞癌(OSCC)細(xì)胞的增殖、遷移和侵襲，因此該融合基因有望成為治療OSCC的有效靶點(diǎn)(Panetal., 2020)。相比于人類等少數(shù)哺乳動(dòng)物，真菌的融合基因研究嚴(yán)重滯后。在蜜蜂病原中還沒有相關(guān)研究報(bào)道。本研究綜合蜜蜂球囊菌菌絲和孢子的PacBio SMRT測(cè)序數(shù)據(jù)共鑒定到921個(gè)融合基因，其中菌絲和孢子共有的融合基因?yàn)?10個(gè)，菌絲和孢子特有的融合基因分別為411和0個(gè)。推測(cè)上述共有融合基因在蜜蜂球囊菌的不同形態(tài)中發(fā)揮重要作用，而菌絲的特有融合基因在菌絲的生長發(fā)育中扮演特殊角色，但仍需要進(jìn)一步探究。進(jìn)一步分析發(fā)現(xiàn)蜜蜂球囊菌菌絲和孢子中的融合基因均主要分布在負(fù)鏈，且主要由染色體AZGZ01000008.1產(chǎn)生，體現(xiàn)出菌絲和孢子中融合基因的共性。

目前僅有禾谷鐮孢菌Fusariumgraminearum(Liuetal., 2016)、大孢指疫霉菌Sclerophthoramacrospora(Teichertetal., 2017)和粗糙鏈孢菌Neurosporacrassa(Liuetal., 2017)等少數(shù)幾種真菌有RNA編輯的研究報(bào)道。在芽殖酵母Saccharomycescastellii、靈芝Ganodermalucidum和松生擬層孔菌Fomitopsispinicola中，大部分RNA編輯事件的堿基替換類型為T->C, G->A, C->T, A->G(Zhuetal., 2014; Wangetal., 2016; Wuetal., 2018)。本研究在蜜蜂球囊菌中共鑒定到738次RNA編輯事件，涉及12種堿基替換類型(圖3: A)。其中最常見的堿基替換類型有4種(C->T, T->C, G->T, C->A)(圖3: B)，與上述其他真菌的編輯類型并不完全一致，表明不同物種中RNA編輯事件的編輯類型具有物種特異性(馬艷莉和俞嘉寧, 2009)，且RNA編輯酶復(fù)合物的識(shí)別具有特異性(Bahnetal., 2012)。此外，還發(fā)現(xiàn)蜜蜂球囊菌RNA編輯的功能類型中同義單核苷酸突變最為豐富(圖2: A)，但在其他已報(bào)道真菌中非同義單核苷酸突變最為常見(Bianetal., 2019)，說明蜜蜂球囊菌在進(jìn)化過程中需要保證基因編碼蛋白的穩(wěn)定性，從而利于生存和繁衍。進(jìn)一步的比較分析發(fā)現(xiàn)，蜜蜂球囊菌菌絲中的不同類型RNA編輯事件數(shù)量普遍多于孢子中的(圖2: B)。此前有研究表明RNA編輯具有階段特異性(Liuetal., 2017)。在禾谷鐮孢菌(Liuetal., 2016)和大孢指疫霉菌(Teichertetal., 2017)中，RNA編輯被證實(shí)參與了子實(shí)體的發(fā)育。在子囊菌中，RNA編輯可能是形成有性孢子所必需的(Liuetal., 2016)。推測(cè)一方面是由于孢子是一種休眠態(tài)，其新陳代謝較菌絲緩慢(Guoetal., 2018)，另一方面是由于在不同的真菌形態(tài)中RNA編輯發(fā)揮的功能不同(王淮和楊健康, 2020)。此外，AaM中C->T類型的RNA編輯事件數(shù)量最多，而AaS中G->T類型最為豐富(圖3: B)。這可能是AaM中催化C->T突變的PPR蛋白含量較高，但背后的分子機(jī)理尚未明確(Chu and Wei, 2020)。本研究發(fā)現(xiàn)，AaM中RNA編輯位點(diǎn)基因可注釋到代謝進(jìn)程、細(xì)胞進(jìn)程和催化活性等19個(gè)GO功能條目，以及內(nèi)吞作用、核糖體和次級(jí)代謝產(chǎn)物的生物合成等11條KEGG通路；AaS中RNA編輯位點(diǎn)基因也可注釋到細(xì)胞部分等24個(gè)功能條目以及RNA轉(zhuǎn)運(yùn)等20條KEGG通路(圖4和5)。上述結(jié)果說明RNA編輯與蜜蜂球囊菌菌絲和孢子的生長發(fā)育、物質(zhì)和能量代謝的潛在關(guān)聯(lián)，值得進(jìn)一步研究。