蜜蜂球囊菌菌絲和孢子中全長轉錄本的差異表達分析

杜 宇, 蔣海賓, 王 杰, 范小雪, 王秀娜, 馮睿蓉, 張文德,隆 琦, 熊翠玲, 鄭燕珍,2, 陳大福,2, 郭 睿,2,*

(1. 福建農林大學動物科學學院(蜂學學院), 福州 350002; 2. 福建農林大學, 蜂產品加工與應用教育部工程研究中心, 福州 350002;3. 福建農林大學生命科學學院, 福州 350002; 4. 福建農林大學, 福建省病原真菌與真菌毒素重點實驗室, 福州 350002)

蜜蜂球囊菌Ascosphaeraapis(簡稱“球囊菌”)是一種致死性真菌病原，專性侵染蜜蜂幼蟲而導致白堊病。該病于1992年在我國全面暴發，并往往發生在涼爽濕潤的季節，與蜂群大量繁殖的春季和秋季相吻合，可導致成年蜜蜂數量和蜂群群勢的急劇下降，嚴重影響蜂產品的可持續發展(Aronstein and Murray, 2010; 趙紅霞等, 2014)。低日齡蜜蜂幼蟲易受球囊菌感染，而成年蜜蜂可作為媒介在蜂群內傳播疾病。球囊菌孢子隨食物進入蜜蜂幼蟲中腸，因此時中腸與后腸隔絕，腸道內低氧環境使孢子處于低水平萌發狀態，直至幼蟲到預蛹的過渡期，中腸與后腸連通，球囊菌進入后腸接觸氧氣后劇烈萌發，菌絲開始大量生長，相繼穿透腸壁和體壁，最終蔓延包裹幼蟲全身，形成白色或棕黑色的堅硬蟲尸(Jensenetal., 2013)。諸多生態因子如溫度、濕度、酸堿度(pH)、二氧化碳(CO2)濃度、輻射及營養素條件等均對球囊菌的生長和白堊病暴發產生一定影響，外部溫度介于(31～35)±0.5℃，相對濕度在80%以上，pH值介于5～7.8，且CO2濃度高于10%的環境條件有利于球囊菌孢子的萌發(梁勤等, 2000, 2001; Vojvodicetal., 2011; Jensenetal., 2013; 董文濱, 2014)。球囊菌侵染可引起宿主的脅迫應答。Aronstein等(2010)利用cDNA-AFLP和qRT-PCR技術研究發現1日齡西方蜜蜂Apismellifera每頭幼蟲攝入1×105個球囊菌孢子后，其攝食率迅速下降，并影響宿主體內免疫信號途徑的激活以及抗菌肽相關基因的轉錄過程。王立立(2014)研究發現意大利蜜蜂Apismelliferaligustica5日齡幼蟲被5×107孢子/mL濃度的球囊菌侵染后，血淋巴內酚氧化酶原級聯反應為主的體液免疫被顯著激活，并催化黑化反應。

此前，由于基因組基因功能注釋信息的缺失，球囊菌的分子生物學研究進展緩慢。Shang等(2016)通過二代測序技術對球囊菌ARSEF 7405菌株進行了基因組測序、組裝和注釋，公布了scafford水平的球囊菌參考基因組(assembly AAP 1.0)，為其分子生物學及組學研究奠定了基礎。近年來，以RNA-seq為代表的第二代高通量測序技術發展迅猛，為探究蜜蜂與球囊菌的互作機制提供了強大技術手段。前人利用二代測序技術對球囊菌進行的轉錄組研究僅有一例報道，Cornman等(2012)對來自培養基的球囊菌菌絲和來自西方蜜蜂幼蟲感染組織的球囊菌菌絲進行了轉錄組測序，功能分析表明球囊菌的差異表達基因(differentially expressed genes, DEGs)參與了交配類型、細胞內信號轉導和應激反應。近年來，筆者所在團隊基于二代測序技術對球囊菌進行了一系列轉錄組研究，例如基于Illumina測序數據組裝和注釋了球囊菌的首個參考轉錄組，并大規模開發了球囊菌的SSR(張曌楠等, 2017)；利用二代轉錄組數據對球囊菌的參考基因組注釋進行了補充和完善(郭睿等, 2019)；通過比較轉錄組分析初步揭示了球囊菌對意蜂幼蟲和中蜂幼蟲的侵染機制(陳大福等, 2017; 郭睿等, 2017)；并在全基因組水平鑒定和分析了球囊菌的微小RNA(microRNA, miRNA)、長鏈非編碼RNA(long non-coding RNA, lncRNA)和環狀RNA(circular RNA, circRNA)(Chenetal., 2018; Guoetal., 2018; 郭睿等, 2018)。二代測序技術雖具有通量和準確性較高的優勢，但同時也存在測序讀段較短(不超過300 bp)的劣勢，需要通過生物信息學算法對測序得到的短讀段進行拼接和組裝，以得到較長的轉錄本(Boldogk?ietal., 2019)。因此，基于二代短讀段數據只能進行基因的定量和差異表達分析，無法開展轉錄本的結構、定量和差異表達分析(Tombáczetal., 2018)。此外，基于二代測序數據的基因定量實際上是對來源于同一個基因的不同剪接異構體的表達量取平均值，因而混淆了不同剪接異構體的真實表達量。因此，在轉錄本水平進行定量和差異表達分析對于深入理解物種的轉錄組的復雜性及動態變化尤為必要。

近年來，以牛津納米孔(Oxford Nanopore)長讀段測序技術為代表的三代測序技術逐漸興起和成熟，因具有超長讀段(平均約15 kb)、較高準確性和直接讀取核酸修飾等優勢而成功應用于人、兔和橄欖實蠅Dacusoleae等物種的全長轉錄組研究(Chenetal., 2017; Bayegaetal., 2018; Workmanetal., 2019)。此外，由于三代測序得到的超長讀段無需拼接便能獲取轉錄本序列全長，使轉錄本的精確定量和差異表達分析成為現實(Weiratheretal., 2017)。目前，真菌的全長轉錄組研究十分有限(Piroonetal., 2018)；對于蜜蜂病原，迄今僅在東方蜜蜂微孢子蟲Nosemaceranae中有過相關研究報道(陳華枝等, 2020a, 2021a)。近期，筆者利用Nanopore長讀段測序技術對球囊菌的純化菌絲(Aam)和純化孢子(Aas)分別進行了轉錄組測序，基于混合數據構建和注釋了球囊菌的高質量全長轉錄組(杜宇等, 2021b)；系統鑒定和分析了Aam和Aas中基因的可變剪接(alternative splicing, AS)和可變腺苷酸化(alternative polyadenylation, APA)(杜宇等, 2021a)。然而到目前為止，球囊菌全長轉錄本定量及差異表達轉錄本(differentially expressed transcripts, DETs)分析未見報道。

本研究在前期研究基礎上對Aam和Aas中全長轉錄本進行定量分析，通過比較分析篩選出DETs并進行功能和通路注釋，進一步對Aam和Aas中DET的結構進行比較，并構建和分析DETs對應蛋白的互作網絡。研究結果可為闡明DETs在球囊菌的孢子萌發、菌絲生長和有性生殖等方面的分子功能提供理論依據和數據基礎。

1 材料與方法

1.1 納米孔長讀段測序數據來源

前期已利用Oxford Nanopore技術對球囊菌Aam和Aas進行轉錄組測序，分別測得6 321 704和6 259 727條原始讀段，經過濾分別得到5 669 436和6 233 159條有效讀段，并分別鑒定9 859和16 795條非冗余全長轉錄本(杜宇等, 2021b)。

1.2 全長轉錄本的表達量計算

利用minimap2軟件(Li, 2018)將1.1節Aam和Aas有效轉錄本讀段與球囊菌參考基因組(assembly AAP 1.0)注釋的已知轉錄本進行序列比對，獲取全長轉錄本與已知轉錄本的對應信息。Nanopore全長轉錄組測序可以模擬成一個隨機抽樣的過程，為了讓片段數目能真實地反映轉錄本表達水平，首先對各樣品中的比對上球囊菌參考基因組的有效讀段數目進行歸一化處理，然后采用CPM(counts per million mapped reads)算法(Zhouetal., 2014)計算每一條全長轉錄本的表達量，公式如下：

CPM=比對到轉錄本上的讀段數(reads mapped to transcript)/樣本中對齊的總讀段數(total reads aligned in sample)×1 000 000。

1.3 DETs的篩選

將同一全長轉錄本在Aam和Aas中表達量的比值定義為相對變化倍數(fold change)。通過對差異顯著性P值進行校正得到錯誤發現率(false discovery rate, FDR)。由于Nanopore全長轉錄組測序的差異表達分析是對大量的全長轉錄本表達量進行獨立的統計假設檢驗，會存在假陽性問題，因此在進行差異表達分析過程中，采用了公認的Benjamini-Hochberg校正方法對原有假設檢驗得到的顯著性P值進行校正，并最終采用FDR作為DETs篩選的關鍵指標。 根據fold change≥2且FDR<0.01的標準，采用EBSeq軟件(Lengetal., 2013)從AasvsAam比較組中篩選DETs。

1.4 DETs的生物信息學分析

利用百邁克云平臺(https:∥international.biocloud.net)的相關工具繪制DETs的曼哈頓圖及表達量聚類熱圖。 利用BLAST工具將DETs比對GO數據庫(https:∥www.geneontology.org)和KEGG數據庫(https:∥www.genome.jp/kegg/)，獲得相應的功能和通路注釋。再根據DETs與KEGG通路的注釋關系構建調控網絡，并利用Cytoscape軟件(Smootetal., 2011)進行網絡可視化。參照Sylvain等(2007)的方法，通過IGV(Integrative Genomics Viewer)瀏覽器對DETs的結構進行可視化。

2 結果

2.1 球囊菌菌絲和孢子中全長轉錄本的表達量分析

Aam和Aas共有的表達轉錄本數量為19 966條，特有的表達轉錄本數量分別為1 273和2 856條。Aam中全長轉錄本的表達量介于0.0010～76 721.3655，其中表達量最高的為ONT.6598.1, ONT.4715.1和ONT.3612.1等(表1)。Aas中全長轉錄本的表達量介于0.0010～20 116.2770，其中表達量最高的為ONT.7346.1, ONT.2508.21和ONT.2150.1(表2)。此外，部分全長轉錄本同時在Aam和Aas中高量表達，例如ONT.705.2(CPM值分別為2 290.0785和9 827.5333), ONT.2558.1(CPM值分別為2 917.3948和6 693.9654)和ONT.4966.2(CPM值分別為6 640.7930和5 544.5046)。

表1 蜜蜂球囊菌菌絲(Aam)中CPM值前10位的全長轉錄本信息概要Table 1 Summary of full-length transcripts in Ascosphaeraapis mycelium (Aam) with the top 10 CPM values

2.2 球囊菌菌絲和孢子中DETs的篩選及表達量聚類

在AasvsAam比較組中共篩選到3 230條DETs，包含3 072條上調轉錄本和158條下調轉錄本。 上述DETs的log2(fold change)介于-6.5506～14.4956。上調倍數最大的全長轉錄本是ONT.4715.1，其次為ONT.7176.3和ONT.3781.2(表3)；下調倍數最大的全長轉錄本是rna5277，其次為ONT. 2275.2和ONT.7416.7(表4)。

表3 蜜蜂球囊菌Aas vs Aam比較組中前10位上調轉錄本的信息概要Table 3 Summary of the top 10 up-regulated transcripts in the Aas vs Aam comparison group of Ascosphaera apis

表4 蜜蜂球囊菌Aas vs Aam比較組中前10位下調轉錄本的信息概要Table 4 Summary of the top 10 down-regulated transcripts in the Aas vs Aam comparison group of Ascosphaera apis

2.3 球囊菌菌絲和孢子中DETs的功能和通路注釋

GO數據庫結果顯示，上述DETs涉及生物學進程(biological process)、細胞組分(cellular component)和分子功能(molecular function) 3個大類，其中有1 251條DETs注釋到生物學進程大類的16個功能條目，包括細胞進程(cellular process)、代謝進程(metabolic process)及單組織進程(single-organism process)等(圖1)；有1 170條DETs注釋到細胞組分大類的14個功能條目，包括細胞(cell)、細胞部分(cell part)及細胞器(organelle)等(圖1)；有1 157條DETs注釋到分子功能大類的12個功能條目，包括催化活性(catalytic activity)、結合(binding)和結構分子活性(structural molecule activity)等(圖1)。

KEGG數據庫注釋結果顯示分別有36, 639和832條DETs注釋到細胞進程(celluar process)、遺傳信息處理(genetic information processing)和新陳代謝(metabolism)相關通路(圖2)。注釋DETs數量最多前5位通路分別為核糖體(ribosome)、抗生素的生物合成(biosynthesis of antibiotics)、碳代謝(carbon metabolism)、氨基酸的生物合成(biosynthesis of amino acids)及糖酵解/糖異生(glycolysis/gluconeogenesis)(圖2)。

進一步分析發現64條DETs注釋到糖酵解/糖異生，包括ONT.1193.1, ONT.1218.4和ONT.29.1等62條上調轉錄本，以及2條下調轉錄本(ONT.7961.1和ONT.7961.5)；ONT.5566.7, ONT.5566.3和ONT.4500.10等26條上調轉錄本注釋到內吞作用。5條上調轉錄本(ONT.1358.7, ONT.1474.3, ONT.750.14, ONT.750.2和rna3824)以及2條下調轉錄本(ONT.5822.1和rna735)注釋到MAPK信號通路；2條上調轉錄本(ONT.750.14和ONT.750.2)同時注釋到內吞作用和MAPK信號通路(圖3)。

圖3 蜜蜂球囊菌Aas vs Aam比較組中差異表達轉錄本(DETs)與糖酵解/糖異生、內吞作用及MAPK信號通路的注釋關系網絡Fig. 3 Annotation relationship networks of glycolysis/gluconeogenesis, endocytosis and MAPK signaling pathway associateddifferentially expressed transcripts (DETs) in the Aas vs Aam comparison group of Ascosphaera apis

2.4 球囊菌菌絲和孢子中DETs的結構比較

對Aam和Aas中DETs的結構進行比較，發現685個DETs(ONT.3612.1, ONT.1081.2, ONT.1218.6和ONT.1534.4等)在菌絲和孢子中除具有不同的表達量外，同時也具有不同的結構，例如全長轉錄本ONT.3612.1由分別在Aam和Aas中發現的剪接異構體ONT.2373.1和ONT.3407.1合并得到，ONT.3612.1在Aam和Aas中的CPM值分別為33 850.23和47.22；此外，ONT.2373.1和ONT.3407.1的結構存在差異，二者均含有1個外顯子，但前者的外顯子明顯短于后者(圖4)。

圖4 球囊菌菌絲(Aam)和孢子(Aas)中差異表達轉錄本 ONT.3612.1結構的IGV瀏覽器視圖Fig. 4 IGV browser view of the structure of differentially expressed transcript ONT.3612.1in Ascosphaera apis mecylium (Aam) and spore (Aas)白色矩形表示scafford，上面的數字表示scafford上的堿基位置；紅色矩形、藍色矩形和綠色矩形表示外顯子，箭頭表示轉錄方向；Final代表在Aam和Aas中鑒定到的轉錄本的合并結果。White rectangle indicates scafford, and the numbers above scafford indicate base position. Red rectangle, blue rectangle and green rectangle indicate exons, and arrows indicate the direction of transcription. Final represents the integration of identified transcripts in Aam and Aas.

3 討論與結論

前期研究中，筆者所在團隊基于二代短讀段數據對球囊菌菌絲和孢子轉錄組中的表達基因進行了定量和比較分析，分別檢測到7 557和7 640個表達基因，篩選出7 362個共有表達基因以及195和278個特有表達基因，并鑒定到977個上調基因和687個下調基因(蔣海賓等, 2020)。由于二代測序產生的短讀段限制，無法對轉錄本進行結構、精確定量和差異表達分析。同一個基因轉錄而來的前體mRNA(pre-mRNA)通過AS可形成不同的剪接異構體，最終形成不同的蛋白質而表現出不同的表型(Modrek and Lee, 2002; Wangetal., 2008)。

相對于菌絲，真菌孢子是一種休眠態，僅保持低水平的新陳代謝和生命活動(Liuetal., 2016)。本研究利用已獲得的Nanopore長讀段測序數據對Aam和Aas中表達的轉錄本進行定量，結果顯示Aam中全長轉錄本的CPM值介于0.0010～76 721.3655(表1)，而Aas中全長轉錄本的CPM值介于0.0010～20 116.2770(表2)，說明球囊菌菌絲中全長轉錄本的整體表達水平更高，與實際情況相符。此前，我們利用二代測序技術在球囊菌孢子中鑒定到392個miRNA(陳華枝等, 2020b)、850個lncRNA(陳華枝等, 2021c)和2 995個circRNA(陳華枝等, 2021b)，并對部分ncRNA的表達進行了驗證，結果說明球囊菌孢子中存在一定水平的基因轉錄活動。 本研究的結果進一步提供了球囊菌孢子中存在轉錄活動的證據。本研究發現，部分轉錄本在菌絲中特異性高表達，例如ONT.1954.6(CPM=1 421.3261), ONT.5520.1(CPM=6 245.3525)和ONT.7176.1(CPM=2 651.4502)等；部分轉錄本在孢子中特異性高表達，例如ONT.7346.1(CPM=20 116.2770), ONT.2508.21(CPM=17 482.8947)和ONT.2150.1(CPM=16 328.9977)等；且ONT.705.2, ONT.2558.1和ONT.4966.2等轉錄本同時在菌絲和孢子中高量表達。推測菌絲中特異性高表達的轉錄本與菌絲的生長和發育息息相關，孢子中特異性高表達的轉錄本與孢子的萌發和生長存在密切關系，而菌絲和孢子均維持高量表達的轉錄本在球囊菌生活史的不同階段皆發揮必要功能。相較于孢子，球囊菌菌絲中分別有3 072和158條全長轉錄本上調和下調表達，上調轉錄本的數量遠多于下調轉錄本，說明多數轉錄本在菌絲中的表達更為活躍，與實際情況相符。本研究還發現一些轉錄本上調和下調的幅度很大，例如ONT.4715.1[Aam中CPM=75 415.9773, Aas中CPM=11.0906, log2(fold change)=14.4956]和ONT.2430.18[Aam中CPM=1.5644, Aas中CPM=364.8805，log2(fold change)=-5.9176]等，鑒于這些轉錄本在球囊菌菌絲和孢子中均維持不低的表達量，值得進一步深入研究。

球囊菌在增殖過程中需要從外界攝取多種碳源、氮源、維生素和礦質元素，其中果糖和葡萄糖分別是球囊菌菌絲體生長和雜交產孢的最佳碳源(梁勤等, 2001)。真菌能通過糖酵解途徑將攝取到的葡萄糖轉化為丙酮酸，從而滿足自身生長的能量需求(劉天明等, 2008)。因此，球囊菌在菌絲生長過程中必然進行活躍的營養攝取和能量轉化活動。糖酵解與糖異生途徑是方向相反的兩個能量代謝途徑，二者的多數酶促反應可逆。糖異生是一種從氨基酸或三羧酸(TCA)循環的中間體等非碳水化合物底物中生成葡萄糖或糖原的代謝過程，在生物體發育和適應營養饑餓方面發揮著重要作用(Marreroetal., 2010)。真菌在面對葡萄糖或其他碳源不足時能夠通過糖異生途徑維持自身發育所需的能量，并在侵染期間通過該途徑應對攝取外源營養不足等情況(Brakhage, 2013)。Liu等(2018)研究發現糖異生途徑相關的PCK1基因可啟動灰霉菌Botrytiscinerea分生孢子萌發和侵染結構形成等過程，并協助灰霉菌侵染寄主細胞以及增強病原毒力。本研究中，64條DETs同時注釋到糖酵解/糖異生相關通路，包括62條上調轉錄本(ONT.1193.1, ONT.1218.4和ONT.29.1等)以及2條下調轉錄本(ONT.7961.1和ONT.7961.5)(圖3)，上調轉錄本數量遠高于下調轉錄本的數量，暗示球囊菌菌絲可通過上調表達糖酵解/糖異生途徑相關轉錄本，影響自身能量轉化過程，以滿足其旺盛的代謝活動需要。但由于糖酵解和糖異生在KEGG數據庫中的pathway信息為糖酵解/糖異生，所以尚無法區分分別注釋到糖酵解和糖異生的DETs。

內吞作用涉及細胞的生長和分化以及分子信號傳遞。真菌菌絲可通過內吞作用將酶類與質膜等關鍵分子運送至菌絲頂端和其他部位以促進菌絲生長和極性延伸(Wangetal., 2019)。羅伯茨綠僵菌Metarhiziumrobertsii體內的MrArk1作為內吞作用的關鍵調控因子參與分生孢子和毒力水平的調節過程(Wangetal., 2019)。稻瘟病菌Magnaportheoryzae內吞作用基因MoARK1和MoPAN1可顯著影響其細胞壁完整性、產孢量以及致病性(王佳妹, 2012)。但內吞作用在球囊菌增殖過程中的生物學功能目前仍未明確。在球囊菌侵染蜜蜂幼蟲腸道過程中，其菌絲分泌各類水解酶抑制或降解宿主的防御因子，并吸收利用宿主的營養物質以促進生長過程(Jensenetal., 2013)。本研究發現ONT.5566.7, ONT.5566.3, ONT.4500.10和ONT.1179.1等26條全長轉錄本注釋到內吞作用(圖3)，且在Aam中均為上調表達。推測上述內吞作用相關全長轉錄本在促進菌絲生長方面發揮一定作用，并與病原毒力水平存在潛在關聯。

絲裂原活化蛋白激酶(mitogen-activated protein kinase, MAPK)級聯反應可引起下游靶基因的表達變化，參與真菌細胞壁完整性、交配繁殖、孢子形成以及致病性等方面的調控過程(Zhaoetal., 2007; Igbariaetal., 2008)。新月彎孢菌Curvularialunata的MAPK信號通路相關基因Clk1和Clm1分別與病原孢子形成以及細胞壁形成密切相關(Nietal., 2018)。筆者所在團隊前期利用二代測序技術發現球囊菌菌絲和孢子中的lncRNA和circRNA差異表達，并分別通過順式(cis)作用和調控來源基因轉錄，影響MAPK信號通路相關基因的表達(陳華枝等, 2021b, 2021c)。本研究發現5條上調轉錄本(ONT.1358.7, ONT.1474.3, ONT.750.14, ONT.750.2和rna3824)以及2條下調轉錄本(ONT.5822.1和rna735)皆可注釋到MAPK信號通路(圖3)。其中，ONT.750.14和ONT.750.2可同時注釋到內吞作用和MAPK信號通路(圖3)。上述結果表明MAPK信號通路相關全長轉錄本可能與球囊菌的孢子形成、細胞壁完整性維持以及雜交產孢密切相關。

三代測序技術為深入研究轉錄本結構提供了強大工具(Boldogk?ietal., 2018, 2019)。基于二代短讀段測序數據只能進行基因表達量的計算和差異表達分析，但基于三代長讀段測序數據不僅能夠同時進行基因和轉錄本表達量的計算和差異表達分析，還能對基因(轉錄本)的結構進行精確的AS和APA分析(Abdel-Ghanyetal., 2016; Moldovánetal., 2018)。相對于較早出現的PacBio SMRT測序技術，Nanopore長讀段測序技術由于誕生較晚，在轉錄本結構方面相關研究還較為有限。本研究發現，部分來源于同一基因的剪接異構體在球囊菌菌絲和孢子中具有不同的結構，例如全長轉錄本ONT.3612.1在菌絲和孢子中經過AS形成兩條表達量和結構均不同的剪接異構體ONT.2373.1和ONT.3407.1(圖4)，表明球囊菌在生長和發育階段不僅伴隨著轉錄本表達量的變化，同時也伴隨著轉錄本結構的改變；同一基因來源的全長轉錄本在菌絲和孢子中經歷不同的AS，形成具有不同結構的剪接異構體。這為深入理解球囊菌的菌絲生長、孢子萌發和有性生殖提供了一個嶄新的維度，也為其他物種的轉錄本結構研究提供了可供參考的思路和方法。

球囊菌的轉基因操作技術平臺缺失導致其基因(剪接異構體)功能研究嚴重滯后。近期，Tauber等(2020)通過體外轉錄合成Ras家族編碼基因和β-葡聚糖合成蛋白編碼基因的雙鏈RNA(dsRNA)并處理球囊菌，發現球囊菌孢子可能在萌發早期吸收上述dsRNA，孢子萌發率明顯降低。該研究為深入開展球囊菌的基因(剪接異構體)的功能研究提供了方法借鑒。