蔣海賓,杜 宇,范小雪,王 杰,祝智威,范元嬋,熊翠玲,付中民,2,徐國鈞,2,陳大福,2,郭 睿,2

(1.福建農(nóng)林大學(xué)動物科學(xué)學(xué)院(蜂學(xué)學(xué)院)， 福州 350002；2. 福建農(nóng)林大學(xué)蜂產(chǎn)品加工與應(yīng)用教育部工程研究中心， 福州 350002)

1 引 言

蜜蜂作為最重要的授粉昆蟲，在農(nóng)業(yè)生產(chǎn)、生態(tài)保護(hù)和生物多樣性維持等方面具有至關(guān)重要的作用.但蜜蜂的群居性特點(diǎn)使其健康易受到多種傳播性病蟲害的威脅.其中，白堊病是一種由蜜蜂球囊菌(Ascosphaeraapis，簡稱球囊菌)侵染蜜蜂幼蟲而導(dǎo)致的蜜蜂真菌病，現(xiàn)已成為影響?zhàn)B蜂生產(chǎn)的一大頑疾.近期，有研究表明球囊菌對中華蜜蜂(Apisceranacerana)雄蜂和工蜂幼蟲[1]、成年熊蜂[2]也具有侵染性.不同性別的球囊菌菌絲通過異宗配合進(jìn)行有性生殖產(chǎn)生孢子，孢子外部具有一層由幾丁質(zhì)和蛋白構(gòu)成的致密孢子壁，能夠抵御外界不良環(huán)境的影響.球囊菌作為寄生蜜蜂幼蟲的專性真菌病原，在與宿主的長期協(xié)同進(jìn)化過程形成了特殊的侵染模式.與白僵菌(Beauveria)[3]和綠僵菌(Metarhizium)[4]等子囊菌門的昆蟲病原真菌不同，球囊菌不能通過菌絲從宿主體表進(jìn)行侵染，其孢子經(jīng)成年蜜蜂的哺育行為被蜜蜂幼蟲經(jīng)口攝入是其侵染宿主的唯一方式，孢子在宿主腸道內(nèi)環(huán)境的刺激下低水平萌發(fā)，并生長為具有致病性的菌絲體，期間分泌幾丁質(zhì)酶、蛋白酶、脂酶及次級代謝物等毒力因子共同抑制宿主的免疫防御系統(tǒng)，進(jìn)而降解圍食膜、腸壁和體壁，最終在預(yù)蛹期穿透尾部，逐漸蔓延包裹整個蟲尸[1].到目前為止，尚缺乏有效治療白堊病的手段.有關(guān)球囊菌菌絲和孢子生長、發(fā)育和生殖的分子機(jī)理仍不明了.因此，在實(shí)驗(yàn)室條件下對球囊菌進(jìn)行純培養(yǎng)，通過比較轉(zhuǎn)錄組的方法分析孢子和菌絲的共有基因、特有基因和差異表達(dá)基因(differentially expressed gene，DEG)，并深入探討上述基因在菌絲和孢子中的潛在作用以及與球囊菌致病性的關(guān)系，可為闡明球囊菌菌絲和孢子的生長、發(fā)育和生殖的分子機(jī)理提供有價值的參考信息和必要的數(shù)據(jù)基礎(chǔ).

前人對球囊菌的研究涉及培養(yǎng)方法、流行病學(xué)、生態(tài)學(xué)、增殖方式、病理學(xué)、疾病防治及免疫防御等方面[5-8].溫度、酸堿度(pH)、維生素等營養(yǎng)生態(tài)條件顯著影響球囊菌的萌發(fā)、生長及產(chǎn)孢過程[9].球囊菌是第一個完成基因組測序的的昆蟲真菌病原體，Qin等[10]于2006年組裝并公布了球囊菌基因組的序列信息，但作者當(dāng)時未公布基因注釋信息.2012年，Cornman等[11]利用Roche 454焦磷酸測序技術(shù)對純培養(yǎng)的球囊菌菌絲和被球囊菌感染的西方蜜蜂(Apismellifera)幼蟲體表長出的菌絲進(jìn)行轉(zhuǎn)錄組測序和比較分析，結(jié)果顯示DEG涉及病原的交配繁殖、應(yīng)激反應(yīng)、毒力因子和信號傳導(dǎo)等通路.2016年，Shang等[12]對包括球囊菌在內(nèi)的7種昆蟲病原真菌的基因組進(jìn)行denovo組裝和注釋，通過系統(tǒng)分析揭示了真菌致病性的演化模式.前期研究中，筆者所在課題組圍繞球囊菌開展了一系列的組學(xué)和分子生物學(xué)研究，基于高質(zhì)量轉(zhuǎn)錄組數(shù)據(jù)組裝和注釋了球囊菌的參考轉(zhuǎn)錄組，并大規(guī)模挖掘了病原的SSR位點(diǎn)[13]；利用鏈特異性建庫的lncRNA-seq和small RNA-seq技術(shù)對純培養(yǎng)的球囊菌菌絲和孢子混合樣品進(jìn)行測序，基于高質(zhì)量的全轉(zhuǎn)錄組數(shù)據(jù)對球囊菌的微小RNA(microRNA,miRNA)[14]、長鏈非編碼RNA(long-non coding RNA,lncRNA)[15]和環(huán)狀RNA(circular RNA,circRNA)[16]及調(diào)控網(wǎng)絡(luò)進(jìn)行了全基因組鑒定和分析.近期，Chen等[17]利用PacBio單分子實(shí)時測序技術(shù)對球囊菌菌絲進(jìn)行了深度測序，基于全長轉(zhuǎn)錄組數(shù)據(jù)重構(gòu)了球囊菌的參考轉(zhuǎn)錄組，可為球囊菌組學(xué)和分子生物學(xué)研究的深入開展提供關(guān)鍵數(shù)據(jù)支撐.

球囊菌的菌絲和孢子是病原的兩種不同形態(tài)，二者的共有基因、特有基因和DEG與病原的生長、發(fā)育、生殖和致病性密切相關(guān).本研究利用鏈特異性建庫的RNA-seq技術(shù)對純培養(yǎng)的球囊菌菌絲和孢子分別進(jìn)行測序，通過比較轉(zhuǎn)錄組分析篩選菌絲和孢子的共有基因、特有基因和DEG，并進(jìn)行GO分類和KEGG通路富集分析，以期明確球囊菌菌絲和孢子的基因表達(dá)差異，挖掘球囊菌生長、發(fā)育、生殖及病原致病性相關(guān)的重要信息.

2 材料和方法

2.1 供試生物材料

球囊菌菌株[13]由福建農(nóng)林大學(xué)動物科學(xué)學(xué)院(蜂學(xué)學(xué)院)蜜蜂保護(hù)學(xué)實(shí)驗(yàn)室分離和保存.

2.2 球囊菌培養(yǎng)及測序樣品制備

將實(shí)驗(yàn)室保存的球囊菌孢子轉(zhuǎn)接至10塊PDA固體培養(yǎng)基上，置于生化培養(yǎng)箱，在33±0.5 ℃條件下培養(yǎng).為避免樣本菌絲和孢子之間的交叉污染，首先在超凈臺中將培養(yǎng)7 d的培養(yǎng)基表面最上層的白色菌絲小心刮至RNA-free的EP管(共計100 mg)，避免接觸與孢子囊接觸的菌絲，然后進(jìn)行顯微制片和鏡檢(不含黑色孢子)，即為球囊菌的純凈菌絲.將菌絲樣品(AaM)放入液氮速凍后迅速移至-80 ℃超低溫冰箱保存?zhèn)溆?

利用接種環(huán)將覆蓋在孢子囊上的薄層菌絲小心刮去，更換接種環(huán)并將剩余孢子囊刮至RNA Free的EP管中，加入適量無菌水震蕩混勻；7 000 g離心3 min，因球囊菌孢子密度較大而沉在底部，棄上清，無菌水重懸孢子沉淀，重復(fù)三次；取少量孢子懸液進(jìn)行顯微制片和鏡檢，400倍視野下可見較多游離的成熟孢子且未見菌絲，即為球囊菌的純凈孢子；將各管孢子集中于一個RNA Free的EP管(共計100 mg)；再將孢子樣品(AaS)放入液氮速凍后迅速移至-80 ℃超低溫冰箱保存?zhèn)溆?

2.3 RNA提取、cDNA文庫構(gòu)建及深度測序

將上述球囊菌菌絲樣品和孢子樣品分別置于兩個干凈的研缽(DEPC水預(yù)處理)，加入液氮充分研磨.利用RNA提取試劑盒(AxyPrepTM Multisource Total RNA Miniprep Kit)(TaKaRa公司，日本)分別提取AaM和AaS的總RNA，去除其中的核糖體RNA(rRNA)，加入Fragmentation buffer將其打斷成短片段，再以短片段的RNA為模板，用六堿基隨機(jī)引物(Random hexamers)合成cDNA第1鏈，加入緩沖液、dNTPs(dUTP代替dTTP)、RNase H2O和DNA polymerase Ⅰ合成cDNA第2鏈，經(jīng)過QiaQuick PCR試劑盒純化并加EB緩沖液洗脫經(jīng)末端修復(fù)、加堿基A，加測序接頭，然后通過UN(Uracil-N-Glycosylase)酶降解第二條鏈.用瓊脂糖凝膠電泳進(jìn)行片段大小選擇，進(jìn)行PCR擴(kuò)增.將建好的cDNA文庫委托廣州基迪奧生物科技有限公司進(jìn)行雙端(paired-end)測序，測序平臺為Illumina HiSeq 4000.測序數(shù)據(jù)已上傳到NCBI SRA數(shù)據(jù)庫，BioProject號：PRJNA560452.

2.4 測序數(shù)據(jù)質(zhì)控及參考基因組比對

為保證測序數(shù)據(jù)的質(zhì)量，對于測序下機(jī)的原始讀段(raw reads)進(jìn)行過濾，去除含有接頭序列的reads，全部都是A的reads，含未知堿基(N)比例大于10%的reads及低質(zhì)量的reads(質(zhì)量值Q≤20的堿基數(shù)占整條reads的50%以上)，得到有效讀段(clean reads).使用短reads比對工具bowtie 2將clean reads比對到核糖體數(shù)據(jù)庫(不允許錯配)，去除比對上的clean reads；將未比對上的clean reads繼續(xù)比對球囊菌基因組(assembly AAP 1.0)，將比對上的clean reads用于轉(zhuǎn)錄本的組裝及分析.

2.5 共有和特有基因的篩選及分析

根據(jù)FPKM(Fragment Per Kilobase of per Million Mapped Reads)法計算基因表達(dá)量，公式為FPKM=Total exon fragment/[Mapped fragment (millions) × exon length (KB)].利用OmicShare在線工具集合(http://www.omichsare.com)對AaM和AaS中的基因進(jìn)行Venn分析，篩選出共有基因和特有基因.利用WEGO軟件對共有基因和特有基因進(jìn)行GO(Gene Ontology)分類.利用Blast工具將共有基因和特有基因分別映射到KEGG(Kyoto Encyclopaedia of Genes and Genomes)數(shù)據(jù)庫(http://www.kegg.jp)，進(jìn)而統(tǒng)計基因顯著性富集的通路(pathway).

2.6 差異表達(dá)基因(DEG)分析

利用edgeR軟件對AaM vs AaS比較組的DEG進(jìn)行篩選，篩選標(biāo)準(zhǔn)為P≤0.05且|log2fold change|≥1.按上述方法對分別對DEG進(jìn)行GO分類和KEGG通路富集分析.

2.7 DEG的RT-qPCR的驗(yàn)證

隨機(jī)選取6個DEG(5個上調(diào)基因和1個下調(diào)基因)進(jìn)行RT-qPCR驗(yàn)證.其中，上調(diào)基因分別為：鐵載體轉(zhuǎn)運(yùn)蛋白酶編碼基因(KZZ91599.1)、核糖體生物發(fā)生蛋白酶編碼基因(KZZ88689.1)、核糖體生物合成ATP酶 RIX7編碼基因(KZZ97232.1)、ABC轉(zhuǎn)運(yùn)體編碼基因(KZZ88904.1)和Fe3+ABC 轉(zhuǎn)運(yùn)蛋白底物結(jié)合蛋白編碼基因(KZZ89047.1)；下調(diào)基因?yàn)椋鸿F載體轉(zhuǎn)錄因子SreA酶編碼基因(KZZ89485.1).根據(jù)所選DEG的序列信息，利用Primer Primer 6軟件設(shè)計特異性的上游(F)引物和下游(R)引物，委托上海生物工程有限公司進(jìn)行引物合成.選擇球囊菌5.8s rRNA基因(GenBank：U68313.1)作為內(nèi)參.利用RNA抽提試劑盒(天漠生物，中國)分別提取孢子和菌絲的總RNA，利用oligo(dT)引物進(jìn)行反轉(zhuǎn)錄得到相應(yīng)的cDNA，作為模板進(jìn)行RT-qPCR.RT-qPCR反應(yīng)體系(20 μL)包含：SYBR Green Dye 10 μL，F(xiàn)引物1 μL，R引物1 μL，DNA模板1 μL，DEPC水補(bǔ)至20 μL.進(jìn)行PCR擴(kuò)增，反應(yīng)條件：95 ℃ 1 min，95 ℃ 15 s，55 ℃ 60 s，共40個循環(huán)，溶解曲線默認(rèn)系統(tǒng)程序.采用2-△△Ct法計算DEG在AaM vs AaS比較組的相對表達(dá)量.每個反應(yīng)進(jìn)行3次生物學(xué)重復(fù)和3次技術(shù)重復(fù).最后利用GraphPad Prism 7軟件進(jìn)行qPCR結(jié)果的檢驗(yàn)及繪圖.本研究使用的引物序列信息詳見表1.

表1 RT-qPCR引物信息

3結(jié)果與分析

3.1 深度測序數(shù)據(jù)的質(zhì)控與評估

球囊菌菌絲(AaM)和孢子(AaS)的RNA-seq分別得到108 614 646和105 675 408條raw reads，經(jīng)質(zhì)控得到的clean reads分別為107 780 032和104 621 402條.AaM和AaS兩端Q30(99.9%堿基正確率)均在95.78%及以上.此外，AaM和AaS的序列比對參考基因組的mapping率分別為85.81%和84.94%.上述結(jié)果共同表明本研究獲得的測序數(shù)據(jù)質(zhì)量良好，可滿足后續(xù)分析.

3.2 共有基因和特有基因的GO和KEGG分析

在AaM和AaS中分別檢測到7 557和7 640個表達(dá)基因.Venn分析結(jié)果顯示共得到7 835個基因，其中AaM和AaS的共有基因數(shù)為7 362個，二者的特有基因數(shù)分別為195和278個.

AaM和AaS的共有基因涉及44個功能條目，其中共有基因富集數(shù)前10位的條目分別為代謝進(jìn)程(1 284)、細(xì)胞進(jìn)程(1 284)、催化活性(1 161)、細(xì)胞(1 055)、細(xì)胞組件(1 054)、單組織進(jìn)程(977)、結(jié)合(836)、細(xì)胞器(700)、大分子復(fù)合物(395)、定位(320).AaM的特有基因涉及13個條目，包括催化活性(4)、細(xì)胞(3)、細(xì)胞組件(3)、結(jié)合(3)、單組織進(jìn)程(2)、細(xì)胞器組件(2)、細(xì)胞器(2)、代謝進(jìn)程(2)、細(xì)胞進(jìn)程(2)、轉(zhuǎn)運(yùn)器活性(1)、細(xì)胞膜(1)、大分子復(fù)合物(1)和定位(1).AaS的特有基因涉及生物學(xué)進(jìn)程和分子功能相關(guān)的6個條目，分別為單組織進(jìn)程(2)、催化活性(2)、多細(xì)胞生物進(jìn)程(1)、發(fā)育進(jìn)程(1)、代謝進(jìn)程(1)和細(xì)胞進(jìn)程(1).括號內(nèi)的數(shù)字代表富集在該條目的基因數(shù).

此外，AaM和AaS的共有基因涉及21類通路，富集共有基因數(shù)最多的是代謝總覽(35.9%)，其次為代謝的輔助因子和維生素(9.1%)、翻譯(8%)、氨基酸代謝(7.5%)、折疊、分類和降解(6.1%)以及運(yùn)輸和分解代謝(5.1%)等(圖1).進(jìn)一步分析發(fā)現(xiàn)，上述21類KEGG通路涉及122條通路，其中在次級代謝物的生物合成(362)、內(nèi)吞作用(79)、甘氨酸、絲氨酸和蘇氨酸的代謝(41)和MAPK信號通路-酵母(33)等通路上富集大量基因.AaM的特有基因富集在51條通路，富集特有基因數(shù)最多的是代謝途徑(13)，其次為泛素介導(dǎo)的蛋白水解(4)和內(nèi)吞作用(4)等.AaS特有的基因富集在55條通路，富集特有基因數(shù)最多是代謝途徑(16)，其次為減數(shù)分裂-酵母(5)和嘧啶代謝(5)等.

3.3 DEG的GO和KEGG分析

AaM vs AaS比較組共含有1 664個DEG，包括977個上調(diào)基因和687個下調(diào)基因.GO分類結(jié)果顯示，上調(diào)基因涉及15條生物學(xué)進(jìn)程相關(guān)條目；11條細(xì)胞組分相關(guān)條目以及6條分子功能相關(guān)條目；下調(diào)基因涉及12條生物學(xué)進(jìn)程相關(guān)條目；12條細(xì)胞組分相關(guān)條目以及8條分子功能相關(guān)條目(圖2).

進(jìn)一步對上述DEG進(jìn)行KEGG通路富集分析，顯著富集在113條通路上，按顯著性排序依次為：糖酵解/糖異生(31)、次級代謝物的生物合成(99)、不同環(huán)境下微生物代謝(65)、氨基酸的生物合成(45)、丙酮酸代謝(21)、代謝途徑(182)、抗生素的生物合成(72)、碳代謝(40)、脂肪酸代謝(14)、苯丙氨酸代謝(9)等(表2).

圖1 球囊菌菌絲和孢子共有基因的KEGG數(shù)據(jù)庫注釋Fig.1 KEGG database annotation of common genes in A.apis mycelium and spore

圖2 AaM vs AaS比較組中DEG的GO數(shù)據(jù)庫注釋Fig.2 GO database annotation of DEGs in AaM vs AaS comparison group

表2 AaM vs AaS比較組中DEG富集數(shù)前20位通路

3.4 DEG的通路富集網(wǎng)絡(luò)構(gòu)建及分析

根據(jù)Nr數(shù)據(jù)庫的蛋白功能注釋信息，篩選內(nèi)吞作用，次級代謝物的生物合成，賴氨酸、絲氨酸和蘇氨酸代謝，MAPK信號通路相關(guān)的DEG，并利用Cytoscape軟件構(gòu)建及可視化上述DEG與通路之間的富集網(wǎng)絡(luò)，分析發(fā)現(xiàn)KZZ98057.1、KZZ87090.1、KZZ90651.1、KZZ90058.1、KZZ90024.1、KZZ92664.1、KZZ98218.1和KZZ92771.1等8個DEG同時富集在甘氨酸、絲氨酸和蘇氨酸代謝及次級代謝物的生物合成，KZZ89869.1同時富集在次級代謝的生物合成與內(nèi)吞作用，而KZZ96974.1同時富集在內(nèi)吞作用與MAPK信號通路(圖3).

圖3 AaM vs AaS比較組中DEG及富集通路的關(guān)系網(wǎng)絡(luò)Fig.3 Relationship network between DEGs and their enriched pathways in AaM vs AaS comparison group

3.5 DEG的RT-qPCR的驗(yàn)證

RT-qPCR檢測結(jié)果顯示6個DEG的表達(dá)水平變化趨勢與測序數(shù)據(jù)中相應(yīng)的表達(dá)水平變化趨勢一致(圖4)，證明了本研究中的轉(zhuǎn)錄組數(shù)據(jù)和基因差異表達(dá)的可靠性.

圖4 DEG的RT-qPCR驗(yàn)證Fig.4 RT-qPCR validation of DEGs

4 討 論

為在全基因組水平最大限度鑒定球囊菌的ncRNA，GUO等[14-16]利用鏈特異性建庫的RNA-seq和sRNA-seq技術(shù)，對球囊菌菌絲和孢子的混合樣品進(jìn)行了全轉(zhuǎn)錄組測序，系統(tǒng)預(yù)測、分析和驗(yàn)證了球囊菌的miRNA、lncRNA和circRNA.由于測序數(shù)據(jù)同時包含菌絲數(shù)據(jù)和孢子數(shù)據(jù)，導(dǎo)致無法準(zhǔn)確區(qū)分僅來源于菌絲的數(shù)據(jù)或僅來源于孢子的數(shù)據(jù)，限制了菌絲生長、孢子萌發(fā)、有性生殖等方面的數(shù)據(jù)挖掘和進(jìn)一步深入分析.本研究通過比較轉(zhuǎn)錄組分析和挖掘球囊菌生長、發(fā)育、生殖及病原致病性相關(guān)的共有基因、特有基因、DEG以及功能條目和通路富集信息，并探討相關(guān)分子機(jī)理.

4.1 球囊菌菌絲和孢子的共有基因與特有基因廣泛參與球囊菌細(xì)胞活動和新陳代謝等諸多生物學(xué)過程

真菌孢子應(yīng)對外界不良環(huán)境時處于低代謝水平的休眠態(tài)，而適宜條件可刺激孢子萌發(fā)生長為代謝旺盛的菌絲體.Aronstein等[5]曾對球囊菌孢子和菌絲的超薄切片進(jìn)行掃描電鏡觀察，發(fā)現(xiàn)二者的細(xì)胞質(zhì)中均存在大量的核糖體和線粒體.本研究發(fā)現(xiàn)，AaM和AaS的共有基因涉及代謝途徑(813)、次級代謝物的生物合成(362)、碳代謝(131)、核糖體(126)及氧化磷酸化(70)等122條通路，表明相關(guān)共有基因廣泛參與和維持了球囊菌菌絲和孢子狀態(tài)下生存所需的物質(zhì)和能量代謝.

對于真菌病原，其孢子具有一層由幾丁質(zhì)和蛋白質(zhì)構(gòu)成的致密孢子壁，能夠有效抵御宿主免疫系統(tǒng)的進(jìn)攻，而菌絲能夠高效利用宿主營養(yǎng)滿足自身繁殖所需的物質(zhì)和能量.幾丁質(zhì)合成酶參與真菌細(xì)胞壁的合成過程，并在維持真菌細(xì)胞形態(tài)、致病性和抗逆性方面發(fā)揮關(guān)鍵作用[18].Kong等[19]研究發(fā)現(xiàn)幾丁質(zhì)合成酶參與了稻瘟病病原(Magnaportheoryzae)附著胞的形成過程，并能維持病原的生長發(fā)育和致病性.球囊菌在侵染意蜂6日齡幼蟲腸道過程中，部分幾丁質(zhì)合成酶相關(guān)基因差異表達(dá)，暗示它們與球囊菌的致病性存在潛在聯(lián)系[20].本研究中，AaM和AaS的13個共有基因涉及幾丁質(zhì)的合成過程，包括幾丁質(zhì)合成酶-1編碼基因TCONS_00003150(log2fold change=0.37,P=8.40E-12)、真菌幾丁質(zhì)合成酶基因KZZ93915.1(log2fold change=-13.98,P=6.71E-99)、幾丁質(zhì)合成酶VI類編碼基因KZZ93066.1(log2fold change=2.10,P=1.04E-09)等，推測上述幾丁質(zhì)酶合成相關(guān)基因在球囊菌孢子和菌絲的生長、發(fā)育、致病性以及細(xì)胞形態(tài)維持等方面具有一定功能.

本研究中，AaM的特有基因有195個，包括細(xì)胞極性蛋白編碼基因KZZ92240.1(log2fold change=-14.22,P=3.57E-215)、P型ATP酶編碼基因KZZ94064.1(log2fold change=-14.72,P=3.60E-225)和細(xì)胞極性蛋白編碼基因KZZ92240.1(log2fold change=-14.22,P=3.57E-215)等，可富集到代謝途徑(13)、泛素介導(dǎo)的蛋白水解(4)、內(nèi)吞作用(4)、次級代謝物的生物合成(4)、核糖體(4)等51個通路；AaS的特有基因?yàn)?78個，包括葡萄糖苷酶編碼基因KZZ93178.1，(log2fold change=15.07,P=2.18E-134)、二烯內(nèi)酯水解酶編碼基因TCONS_00000462(log2fold change=16.08,P=2.90E-84)和含F(xiàn)HA結(jié)構(gòu)域蛋白編碼基因KZZ90006.1(log2fold change=14.70076481,P=2.33E-72)等，可富集到代謝途徑(16)、減數(shù)分裂-酵母(5)、嘧啶代謝(5)、次級代謝物的生物合成(5)、MAPK信號通路-酵母(3)等55個通路.上述結(jié)果表明菌絲和孢子的特有基因在球囊菌的不同生長階段發(fā)揮不同的作用，并通過調(diào)節(jié)物質(zhì)和能量合成及代謝、生殖相關(guān)通路滿足菌絲和孢子的生長、發(fā)育和生殖過程的物質(zhì)和能量需求.

4.2 球囊菌菌絲和孢子的DEG參與信號轉(zhuǎn)導(dǎo)及次級代謝產(chǎn)物的等關(guān)鍵途徑

球囊菌孢子和菌絲的DEG參與調(diào)控代謝進(jìn)程、細(xì)胞進(jìn)程、應(yīng)激反應(yīng)以及催化活性等過程，其中多糖分解進(jìn)程(GO:0000272)、長鏈脂肪酸代謝進(jìn)程(GO:0001676)、氨基酸跨膜轉(zhuǎn)運(yùn)(GO:0003333)等GO條目顯著富集.另一方面，KEGG通路富集分析發(fā)現(xiàn)DEG顯著富集在糖酵解/糖異生，次級代謝物的生物合成以及氨基酸的生物合成等關(guān)鍵通路，并參與甲烷代謝、氧化磷酸化、氮代謝等能量代謝途徑，以及脂代謝、磷酸戊糖途徑、氨基酸代謝等物質(zhì)代謝途徑.

昆蟲的病原真菌可通過分泌絲裂原活化蛋白激酶(MAPK)、絲氨酸/蘇氨酸蛋白激酶和金屬蛋白酶等蛋白酶，配合各類脂酶以及次級代謝產(chǎn)物共同抑制或降解宿主防御分子，并維持真菌的細(xì)胞壁完整性、雜交產(chǎn)孢以及毒力水平[21].當(dāng)病原真菌受到外界刺激后，其MAPK信號級聯(lián)反應(yīng)可引起下游靶基因的表達(dá)變化，從而調(diào)控真菌的交配繁殖和孢子的形成過程[22].本研究發(fā)現(xiàn)，絲裂原活化蛋白激酶MKC1編碼基因KZZ86589.1(log2fold change=2.25,P=6.36E-13)、Db1結(jié)構(gòu)域蛋白編碼基因KZZ92335.1(log2fold change=1.22,P=3.82E-40)、Phox/Bem1p編碼基因KZZ92336.1(log2fold change=1.35,P=4.24E-50)和小GTPase超家族編碼基因KZZ96974.1(log2fold change=13.89,P=5.15E-15)富集在MAPK信號通路且均顯著上調(diào)表達(dá)，說明上述4個基因在球囊菌的生長發(fā)育過程中被顯著激活，進(jìn)而通過調(diào)節(jié)MAPK信號通路影響球囊菌的生長、發(fā)育和生殖.絲氨酸/蘇氨酸蛋白激酶作為真菌細(xì)胞生長、代謝和自噬的關(guān)鍵調(diào)控因子，主要通過調(diào)節(jié)核糖體生物合成、營養(yǎng)攝入和代謝進(jìn)程以維持真菌細(xì)胞內(nèi)穩(wěn)態(tài)[23].李遠(yuǎn)東等[24]在M.robertsii中發(fā)現(xiàn)敲除絲氨酸/蘇氨酸蛋白激酶Ark1基因可抑制其毒力水平并降低其產(chǎn)孢量.本研究中，天冬氨酸激酶結(jié)構(gòu)域蛋白編碼基因KZZ98218.1(log2fold change=1.39,P=1.22E-112)、絲氨酸羥甲基轉(zhuǎn)移酶編碼基因KZZ92664.1(log2fold change=1.05,P=8.80E-57)和色氨酸合成酶編碼基因KZZ87090.1(log2fold change=1.02,P=2.87E-23)等12個DEG富集在甘氨酸、絲氨酸/蘇氨酸通路，暗示絲氨酸/蘇氨酸代謝途徑對球囊菌的孢子形成、菌絲生長具有調(diào)控作用，與球囊菌的毒力和致病性存在潛在關(guān)聯(lián).

次級代謝產(chǎn)物能顯著影響真菌病原的形態(tài)分化、孢子萌發(fā)和致病性過程[25].前人研究發(fā)現(xiàn)玉米赤霉烯酮作為禾谷鐮刀菌(Fusariumgraminearum)的次級代謝產(chǎn)物，可誘導(dǎo)F.graminearum分生孢子周圍菌落的產(chǎn)生過程.擬幼蟲芽孢桿菌(Paenibacilluslarvae)是一種嚴(yán)重威脅西方蜜蜂幼蟲健康的致死性革蘭氏陽性菌，可導(dǎo)致傳染性極強(qiáng)的美洲幼蟲腐臭病.P.larvae可合成與分泌非核糖體肽和肽-聚酮混合物等次級代謝產(chǎn)物對其他細(xì)菌和真菌產(chǎn)生抑制，從而在宿主營養(yǎng)物質(zhì)的爭奪中占據(jù)絕對優(yōu)勢[26].本研究中，葡萄糖苷酶編碼基因KZZ93178.1(log2fold change=15.07,P=2.18E-134)、組氨醇磷酸酶編碼基因KZZ90682.1(log2fold change=11.72,P=3.01E-06)和核糖/半乳糖異構(gòu)酶編碼基因KZZ87542.1(log2fold change=2.54,P=8.13E-143)等99個DEG富集在次級代謝產(chǎn)物的生物合成通路，表明上述DEG通過表達(dá)水平的調(diào)節(jié)對葡萄糖苷酶、組氨醇磷酸酶和丙酮酸激酶等各類次級代謝產(chǎn)物的合成進(jìn)行調(diào)節(jié)，從而調(diào)控球囊菌菌絲形態(tài)分化和孢子萌發(fā)；上述次級代謝產(chǎn)物相關(guān)DEG與球囊菌的致病性存在潛在聯(lián)系，并在其與其他病原微生物的種間競爭中可能扮演特定角色.需要強(qiáng)調(diào)的是，本研究的轉(zhuǎn)錄組數(shù)據(jù)來源是實(shí)驗(yàn)室條件下獲得的球囊菌純培養(yǎng)，從中挖掘出的相關(guān)信息僅能推測與病原致病性和毒力的潛在關(guān)聯(lián).下一步將基于球囊菌侵染蜜蜂幼蟲腸道的全轉(zhuǎn)錄組數(shù)據(jù)，通過連續(xù)比對核糖體數(shù)據(jù)庫、宿主基因組和病原基因組可獲得處于侵染狀態(tài)的球囊菌的組學(xué)數(shù)據(jù)，基于該組學(xué)數(shù)據(jù)分析得到的病原致病性和毒力相關(guān)信息更為真實(shí)可靠.

根癌農(nóng)桿菌介導(dǎo)(AtMT)的真菌遺傳轉(zhuǎn)化為真菌基因的分離克隆和功能研究提供了新的途徑.因其具有成本低、效率高、操作易和重復(fù)性好等優(yōu)點(diǎn)，而廣泛應(yīng)用于白僵菌、綠僵菌和蠟蚧輪枝菌(VerticilliumLecanii)等昆蟲病原真菌的基因功能研究[27].WANG等[28]利用根癌農(nóng)桿菌介導(dǎo)技術(shù)對稻瘟病菌的致病基因MoCPS1進(jìn)行研究，發(fā)現(xiàn)MoCPS1可顯著影響稻瘟病菌的產(chǎn)孢及分生孢子的形態(tài)發(fā)生過程.目前，球囊菌的分子生物學(xué)研究相對滯后，其遺傳轉(zhuǎn)化技術(shù)體系還未建立起來，嚴(yán)重限制了球囊菌的基因功能研究.本研究篩選出的共有基因、特有基因和DEG的數(shù)量依然較多，下一步將結(jié)合球囊菌菌絲、孢子的miRNA、lncRNA、circRNA組學(xué)數(shù)據(jù)構(gòu)建ncRNA與mRNA之間的調(diào)控網(wǎng)絡(luò)，并根據(jù)本研究關(guān)注的重點(diǎn)通路及富集基因進(jìn)一步構(gòu)建和分析局部調(diào)控網(wǎng)絡(luò)，篩選出若干居于核心位置的關(guān)鍵基因以及球囊菌菌絲和孢子的特異性表達(dá)基因，進(jìn)而利用根癌農(nóng)桿菌介導(dǎo)的真菌遺傳轉(zhuǎn)化技術(shù)深入開展相關(guān)功能研究.

綜上所述，本研究通過比較轉(zhuǎn)錄組分析揭示球囊菌菌絲和孢子的共有基因、特有基因和DEG可能通過調(diào)節(jié)細(xì)胞活動，甘氨酸、絲氨酸和蘇氨酸代謝，次級代謝產(chǎn)物合成，MAPK信號通路，以及內(nèi)吞作用等關(guān)鍵途徑參與球囊菌菌絲和孢子的生長、發(fā)育和生殖過程，研究結(jié)果為探明相關(guān)分子機(jī)制提供了基礎(chǔ).