基于RNA-seq 數據分析長牡蠣(Crassostrea gigas)天然免疫的可變剪接事件*

李若暉 張琳琳 ①

(1. 中國科學院海洋研究所實驗海洋生物學重點實驗室 山東青島 266071; 2. 青島海洋科學與技術試點國家實驗室海洋生物學與生物技術功能實驗室 山東青島 266237; 3. 中國科學院大學 北京 100049)

免疫系統是生物體為了抵御病原體入侵, 在長期進化過程中形成的一套復雜的防御系統, 它可以維持機體的內穩態(蒲新明, 2008)。免疫系統通過一系列免疫應答反應參與免疫調節(王金星等, 2004)。免疫系統可以分為天然免疫系統(Innate immune system)和獲得性免疫系統(Adaptive immune system)兩種類型, 兩種類型對應著不同的作用機制(李穎, 2011)。普遍認為無脊椎動物只存在非特異性的天然免疫, 缺乏獲得性免疫(Litmanetal, 2005)。無脊椎動物應對病原微生物的入侵主要依賴天然免疫系統。為應對各種各樣的病原微生物入侵, 基因組需要通過不同的調控產生特定轉錄組和蛋白質組, 從而完成對病原體入侵做出有效和快速的天然免疫反應。

在漫長的自然進化過程中無脊椎動物進化出一套功效強大的先天防御體系, 通過細胞代謝的重編程和免疫系統的激活使它們能夠有效地應對病原體入侵。許多細胞過程的正常進行在很大程度上取決于RNA 代謝, 特別是RNA 的成熟和降解。在真核生物中, 利用不同的剪接位點, 將相同的前體 mRNA 轉化成多個信使 RNA (mRNA)的過程, 被稱作可變剪接(Alternative splicing, AS) (Kahlesetal, 2016)?？勺兗艚油ㄟ^使基因產生結構和功能不同的mRNA 和蛋白質, 從而增強轉錄組可塑性和蛋白質功能多樣性(Maniatisetal, 2002; Nilsenetal, 2010)?？勺兗艚涌梢詮臄盗肯鄬^少的基因中產生大量的mRNA 和蛋白質異構體?？勺兗艚舆€可以改變轉錄本及其編碼蛋白的結構, 從而調節蛋白質的功能, 而可變剪接產生的異構體甚至可能具有完全相反的功能(Stammetal,2005)。與哺乳動物中大量關于可變剪接的研究相比,這一機制在無脊椎動物天然免疫中的作用和調控卻知之甚少。

近年來, 高通量技術和生物信息學工具的快速發展, 特別是二代代測序技術, 如Illumina、PacBio和Oxford Nanopore 直接RNA 測序, 對研究基因表達變化和分析全基因組可變剪接事件做出巨大貢獻。目前, 利用RNA-Seq 技術分析了包括人(Homosapiens)、線蟲 (Caenorhabditiselegans)、果蠅(Drosophila melanogaster)、擬南芥(Arabidopsisthaliana)等多個物種基因組的AS 事件(Panetal, 2008; Graveleyetal,2011; Ramanietal, 2011; Marquezetal, 2012)。利用生物信息學工具進一步深入挖掘轉錄組學數據來識別可變剪接事件, 將有助于闡明無脊椎動物在生物應激條件下的基因調控機制。

長牡蠣(Crassostreagigas)既是世界上重要的經濟貝類, 也是一種常見的模式無脊椎動物, 已被廣泛用于研究(Zhangetal, 2012, 2015)。已有研究表明, AS在長牡蠣響應高溫、高鹽和空氣暴露等非生物應激中發揮重要作用(Huangetal, 2016), 而AS 在長牡蠣生物應激中的功能研究仍處于空白狀態。因此, 本研究基于Illumina 高通量測序的轉錄組數據(Zhangetal,2015), 分析了長牡蠣在不同病原體和混合弧菌不同感染時間誘導下的可變剪接差異, 發現混合弧菌誘導可變剪接事件的增加, 分析了在不同病原體和混合弧菌誘導中差異表達的可變剪接基因的功能, 從而探究可變剪接對長牡蠣天然免疫多樣性和特異性的調控作用。本研究為長牡蠣等無脊椎動物天然免疫多樣性和特異性的深入研究提供了理論依據, 為后續長牡蠣病害防御提供了數據支持。

1 材料與方法

1.1 數據收集與處理

收集來自NCBI 的13 個長牡蠣測序數據, 此數據由中國科學院海洋研究所2013 年上傳發布, 登錄號為PRJNA194079 和PRJNA194084。具體來說, 測序樣品來自: (1) 5 個混合弧菌感染(感染后0、6、12、24 和48 h 收集的感染組的鰓); (2) 8 個不同致病菌[4個革蘭氏陽性菌: 河口弧菌(Vibrioaestuarianus)、鰻弧菌(Vibrioanguillarum) 、 溶藻弧菌(Vibrio alginolyticus)和LPS, 1 個革蘭氏陰性菌: 藤黃微球菌(Micrococcusluteus), PBS 和對照組的鰓]。

為了保證信息分析質量, 對原始數據利用FastQC質控, 使用Trimmomatic 去除原始數據中帶接頭的、低質量的以及帶N 堿基的序列, 得到Clean reads 數據。

1.2 可變剪接事件的鑒定

將Clean reads 與長牡蠣基因組利用軟件STAR進行比對, 比對結果運用Cufflinks 軟件(Trapnelletal,2012)進行組裝, 并使用Cuffcompare 與長牡蠣參考注釋進行比較生成tmap 文件(Trapnelletal, 2010)。使用ASprofile 中的“extract-as”程序(Floreaetal, 2013), 利用gtf、tmap 文件在多個樣本中鑒定12 種可變剪接事件,包括單外顯子跳躍(skipped exon, SKIP)、單外顯子跳躍(模糊邊界) (approximate SKIP, XSKIP)、多外顯子跳躍(multi-exon SKIP, MSKIP)、多外顯子跳躍(模糊邊界) (approximate MSKIP, XMSKIP)、單內含子滯留(intron retention, IR)、單內含子滯留(模糊邊界)(approximate IR, XIR)、多內含子滯留(multi-IR,MIR)、多內含子滯留(模糊邊界)(approximate MIR,XMIR)、可變5′或3′端剪切(alternative exonends, AE)、可變 5′或 3′端剪切(模糊邊界) (approximate AE,XAE)、第1 個外顯子可變剪切(alternative 5′ first exon,TSS)、最后1 個外顯子可變剪切(alternative 3′ last exon, TTS)。

1.3 先天免疫相關基因分析及GO 功能富集分析

使用cuffdiff 程序對混合弧菌不同感染時間和不同病原體樣本進行不同表達基因分析(Trapnelletal,2013)。然后使用參數|log2(fold_change)|＞1 和P≤0.05和FDR≤0.05 對生成的差異表達基因數據進行過濾。

利用DAVID 數據庫, 對顯著差異表達的可變剪接基因ID 進行GO 分類, 獲取GO 編號, 并利用cluster Profiler 對其進行功能富集, 以P≤0.05 作為顯著富集標準。

1.4 預測的長牡蠣AS 基因亞型的表達分析

根據Cuffnorm 利用默認參數對每個樣本中鑒定出的異構體進行定量。通過Cuffdiff 計算混合弧菌感染條件下FPKM 的轉錄表達量, 然后對所有樣本進行歸一化處理。

2 結果與分析

2.1 長牡蠣感染混合弧菌后不同時間可變剪接事件的鑒定

為了研究混合弧菌感染下不同時間長牡蠣可變剪接的作用, 從 NCBI 數據庫下載了 5 個長牡蠣Illumina 高通量測序平臺的轉錄組測序數據。這些數據集來自5 個混合弧菌不同感染時間的RNA-seq 讀數, 包括0、6、12、24 和48 h, 數據集的總堿基讀數為15.27 Gb。

根據5 個NCBI RNA-seq 數據集, 計算了5 個不同感染時間的AS 事件的數量(表1), 其中, 感染48 h后有最豐富的AS 事件(80 566), 而感染6 h 后的數量最少(74 573) (圖1)。在經混合弧菌感染后, 長牡蠣發生的可變剪接事件類型仍以TSS、TTS、AE 和SKIP為主, 結果表明, 可變剪接在各時間點的數目存在差異, 隨著時間的增長, 可變剪接事件也隨著增多。

圖1 長牡蠣在混合弧菌脅迫下不同時間點的可變剪接事件Fig.1 Alternative splicing events of C. gigas under Vibrio challenge at different time points

表1 長牡蠣混合弧菌入侵不同時間點的可變剪接事件Tab.1 Alternative splicing events of C. gigas under Vibrio challenge at different time points

2.2 長牡蠣感染不同病原體可變剪接事件的鑒定

為了研究不同病原體長牡蠣可變剪接的作用,從NCBI 數據庫下載了7 個長牡蠣Illumina 高通量測序平臺的轉錄組測序數據。這些數據集來自不同病原體入侵的長牡蠣的RNA-seq 讀數, 包括LPS、河口弧菌(Vibrioaestuarianus)、鰻弧菌(Vibrioanguillarum)、溶藻弧菌(Vibrioalginolyticus)、塔氏弧菌(Vibrio tubiashii)、藤黃微球菌(Micrococcusluteus)。

根據7 個NCBI RNA-seq 數據集, 計算了7 個不同病原體誘導的AS 事件的數量(表2)。結果表明, 不同病原體誘導后長牡蠣的可變剪接的總數目存在差異, 但長牡蠣發生的可變剪接事件類型仍為TSS、TTS、AE 和SKIP 四種主要類型(圖2)。

圖2 長牡蠣不同病原體入侵下的可變剪接事件Fig.2 Alternative splicing events of C. gigas challenged with different types of bacteria

表2 長牡蠣不同病原體入侵下的可變剪接事件Tab.2 Alternative splicing events of C. gigas challenged with different types of bacteria

2.3 可變剪接事件差異表達分析

對感染混合弧菌的長牡蠣不同時間的可變剪接基因進行差異表達分析, 分別得到528、330、302 和313 個顯著差異表達基因。對不同病原體誘導時發生可變剪接的基因進行差異表達分析, 分別得到528、330、302 和313 個顯著差異表達基因。在這些顯著差異基因中, TSS 和TTS 是主要的可變剪接方式, 這表明它們在響應弧菌的基因表達中可能起著重要的作用。

發生可變剪接事件的差異表達基因有多個基因為天然免疫相關基因家族成員, 混合弧菌不同感染時間發生可變剪接的顯著差異表達基因中分別有38、28、30 和31 個天然免疫基因(表3), 不同病原體入侵發生可變剪接的顯著差異表達基因中分別有25、26、25、23、29 和29 個天然免疫基因(表4), 這些基因家族包括RLR、TLR、SRCR、C1q等。但是混合弧菌不同感染時間發生可變剪接的天然免疫相關基因分別有450、466、482 和501 個, 占長牡蠣天然免疫相關基因總數的三分之一左右, 說明在病原體感染時, 機體中三分之一的天然免疫基因通過可變剪接發揮作用。

表3 弧菌誘導發生可變剪接的差異基因中的天然免疫相關基因Tab.3 Natural immunity-related genes of differentially expressed alternative splicing genes under Vibrio challenge at different time points

表4 不同病原體入侵發生可變剪接的顯著表達差異基因中的天然免疫相關基因Tab.4 Natural immunity-related genes of differentially expressed alternative splicing genes challenged with different types of bacteria

2.4 差異可變剪接基因的GO 富集分析

為了解可變剪接對長牡蠣先天免疫的調控, 對弧菌誘導的5 個不同時間基因組中發生AS 事件的顯著差異表達基因進行了GO 富集, 以P≤0.05 為顯著富集標準, 發現得到143、179、153 和164 個條目(圖3)。

圖3 弧菌誘導下不同時間差異表達可變剪接基因GO 富集Fig.3 Statistics of GO annotations for alternative splicing genes with variable expression in response to Vibrio challenge at different points

在0 h 和6 h 的差異基因中, 生物學過程、細胞構成和分子功能中顯著富集的GO 條目分別有36、35、72 個。在生物學過程中, 發生可變剪接的差異基因主要集中于有機物的分解代謝過程(GO: 1901575)等代謝相關過程; 在分子功能中, 較多可變剪接基因被富集于催化活性(GO: 0003824)相關GO 項中; 在細胞構成中, 發生可變剪接的基因主要集中于膜的固有成分(GO: 0031224)與膜相關條目。

在0 h 和12 h 的差異基因中, 在生物學過程、細胞構成和分子功能中顯著富集的 GO 條目分別有101、20、58 個。在生物學過程中, 發生可變剪接的差異基因主要集中于α-氨基酸代謝過程(GO: 1901605)等代謝相關過程; 在分子功能中, 較多可變剪接基因被富集于輔助因子結合(GO: 0048037)相關GO 項中; 在細胞構成中, 發生可變剪接的基因主要集中于膜的固有成分(GO: 0031224)與膜相關條目。

在0 h 和24 h 的差異基因中, 在生物學過程、細胞構成和分子功能中顯著富集的 GO 條目分別有101、10、72 個。在生物學過程中, 發生可變剪接的差異基因主要集中于有機物的分解代謝過程(GO:1901575)等代謝相關過程; 在分子功能中, 較多可變剪接基因被富集于催化活性(GO: 0003824)相關GO項中; 在細胞構成中, 發生可變剪接的基因主要集中于膜的固有成分(GO: 0031224)與膜相關條目。

在0 h 和48 h 的差異基因中, 在生物學過程、細胞構成和分子功能中顯著富集的 GO 條目分別有112、13、39 個。在生物學過程中, 發生可變剪接的差異基因主要集中于單一有機體代謝過程(GO:0044710)等代謝相關過程; 在分子功能中, 較多可變剪接基因被富集于催化活性(GO: 0003824)相關GO項中; 在細胞構成中, 發生可變剪接的基因主要集中于膜的固有成分(GO: 0031224)與膜相關條目。

由此可見, 這些差異基因主要包括以下主要生物過程: (1) 細胞和代謝過程; (2) 生物過程的調控;(3) 細胞分泌物的產生; (4) 生物黏附過程; (5) 對刺激反應響應先天性免疫系統的啟動; (6) 先天性免疫信號轉導; (7) 細胞的生長等。其中與免疫相關的比例占據較小, 細胞和代謝過程以及生物調節過程的比例占據較高。這說明在長牡蠣生物應激過程中, 細胞、代謝過程、催化活性和生物調節等途徑可能參與機體自身的免疫防御, 與機體的免疫途徑共同響應病原體感染。

2.5 差異可變剪接基因的KEGG 富集分析

對混合弧菌入侵不同時間發生可變剪接的差異表達基因進行KEEG 富集分析(圖4), 結果表明富集到與免疫相關通路, 在12 h 富集到細胞凋亡(ko04210)、CD 分子(ko04090)、信號分子與相互作用(ko09133)、細胞生長和死亡(ko09143)等; 在24 h 富集到NOD-like受體信號通路(ko04621)和免疫系統(ko09151)等。

圖4 弧菌誘導下不同時間差異表達可變剪接基因KEGG 富集Fig.4 Statistics of KEGG annotations for alternative splicing genes with variable expression in response to Vibrio challenge at different points

對不同病原體發生可變剪接的差異表達基因進行KEEG 富集分析(圖5), 結果表明富集到與免疫相關通路, 在LPS 富集到PI3K-AKT 通路(ko04151)等;在V.aes富集到 NF-kappaB 細胞信號傳導通路(ko04064); 在V.ang富集到Toll 與Imd 信號通路(ko04624)、NOD 樣受體信號通路(ko04621)、免疫系統(ko09151)和 NF-kappaB 細胞信號傳導通路(ko04064); 在V.alg富集到NOD 樣受體信號通路(ko04621); 在V.tub富集到免疫系統(ko09151); 在M.lut富集到NOD 樣受體信號通路(ko04621)、Toll與Imd 信號通路(ko04624)。

圖5 不同病原體誘導下差異表達可變剪接基因KEGG 富集Fig.5 Statistics of KEGG annotations for alternative splicing genes with variable expression in response to various Vibrio challenges

綜合差異表達篩選、GO 富集分析和KEGG 富集分析發現, 在混合弧菌感染6 h 時富集到代謝相關通路(圖4a), 推測在感染初期, 機體對能量的需求增加,迫使體內養分重新分配, 積極提供應激所需的能量物質和氨基酸; 在混合弧菌感染12 h 時富集到細胞凋亡相關通路(圖4b), 推測在感染中期, 由于混合弧菌對長牡蠣鰓組織的急性感染, 導致鰓細胞出現凋亡, 以至于在感染6 h、12 h 發生可變剪接事件總個數出現減少; 在混合弧菌感染24 h 時富集到免疫相關通路、48 h 時富集到異源物質代謝降解相關通路(圖4c, 4d), 推測隨著感染時間推移, 鰓組織的天然免疫系統發揮作用, 消除了混合弧菌的急性感染, 并且及時清除機體內由于病原體入侵產生的異物, 逐步恢復鰓組織的免疫穩態。

2.6 混合弧菌脅迫下先天免疫基因的可變剪接事件異構體變化

天然免疫在進化上是一個相對保守的系統, 主要通過PRRs 識別PAMPs 進行信號傳遞。RLRs 是重要的RNA 傳感器分子。對兩個RLR 基因的分析顯示,Cg02656 和Cg26493 基因都擁有六個可變剪接異構體。圖6a 和圖7a 顯示了這兩個基因的異構體經歷了多種可變剪接事件, 表明可變剪接的調節可能非常復雜和多樣, 例如, 異構體Cg02656E 出現了TSS、IR和TTS。另外, 為了探索特定可變剪接異構體與致病菌脅迫之間的關系, 計算了不同感染時間可變剪接轉錄本的表達量(圖6b 和圖7b)。RNA-seq 數據的FPKM 值表明, Cg02656 和Cg26493 基因的可變剪接異構體表達量在不同的感染時間存在差異表達。Cg02656 和Cg26493 基因的表達量均在48 h 后達到最高值, 其中Cg02656 基因在48 h 時的表達量約是感染0 h 時的2 倍。

圖6 Cg02656 基因的AS 異構體(a)及其FPKM 表達量(b)Fig.6 The AS isoforms of the Cg02656 gene (a) and their FPKM expressions (b)

圖7 Cg26493 基因的AS 異構體(a)及其FPKM 表達量(b)Fig.7 The AS isoforms (a) of the Cg26493 gene and their FPKM expressions (b)

基于這些發現, 我們認為長牡蠣可能具有復雜的可變剪接調節系統, 并且可變剪接可能是長牡蠣產生天然免疫多樣性和特異性的重要途徑。例如, 這些異構體可能作為Cg02656 和Cg26493 編碼的完整蛋白的調節因子, 直接參與免疫反應或其他重要的生物學功能。不同感染時間中異構體表達的差異也表明可變剪接調控在不同免疫過程之間存在差異。

3 討論

機體的免疫是在長期進化過程中形成的, 免疫機制復雜性與各種生物大分子的種類與數量有關,而在有限長度的遺傳信息基礎上產生更多的轉錄本,無疑能在轉錄水平上大幅度提高免疫反應的效率,進而實現對機體最大程度的保護。目前已有大量研究指出許多動物的天然免疫與可變剪接高度相關(Tanetal, 2019; Guoetal, 2022), 這種轉錄后水平上的調控與免疫應答中多種關鍵因子的生物合成過程密不可分。此外, 這種調控方式也會造成物種之間代謝通路的特異性。

本研究對感染混合弧菌后長牡蠣的鰓組織中的可變剪接事件進行檢測。隨著時間的增長, 可變剪接事件的數量也均呈現增加的趨勢, 表明混合弧菌的感染會誘導長牡蠣可變剪接事件的產生。關于脅迫反應引起可變剪接事件數量增加的現象在長牡蠣中已有相關報道。Huang 等(2016)的研究結果表明在長牡蠣在高溫、高鹽、空氣暴露脅迫下, 可變剪接事件的數量都有所增加。由此可見, 可變剪接事件的產生可能是長牡蠣響應各種脅迫反應的一種重要調控機制。在本研究中, 12 種類型的可變剪接事件均被鑒定出,其中可變剪接事件數量最多的是TSS 和TTS, 其次是AE 和SKIP。

本研究識別了感染混合弧菌前后長牡蠣鰓組織中的發生可變剪接的差異表達基因, 并對這些基因進行了富集分析。其中GO 的富集分析結果主要包含天然免疫應答、代謝過程。值得注意的是, 長牡蠣感染混合弧菌后12 h 和24 h 發生可變剪接的差異表達基因富集到了免疫系統通路上, 但在6 h 的樣本中并未發現此類富集, 這表明長牡蠣感染后可變剪接事件不會迅速響應。天然免疫應答是非特異性的、機體固有的、反應迅速的、廣泛發生的免疫方式, 表明在感染混合弧菌后長牡蠣主要通過天然免疫應答清除入侵病原菌, 對機體進行保護。有機物代謝過程等GO 詞條被認為參與了機體能量代謝過程, 表明免疫應答過程會進行大量能量消耗與物質代謝。在KEGG富集分析結果中, 發生可變剪接的差異表達基因被顯著富集在細胞凋亡、NOD-like 受體信號通路和免疫系統等信號通路上。細胞凋亡是為了維持內環境的穩定, 由基因控制的非被動的細胞死亡, 并不只在機體受到脅迫時才發揮作用, 也是正常細胞代謝的生理現象之一(Elmore, 2007)。細胞凋亡也是免疫系統的一種作用方式。Nagata 等(2017)的研究表明, 許多凋亡調節因子均在免疫系統中發揮著重要作用。NODlike 受體(NLR)信號通路是生物體的一種重要免疫防御機制。作為天然免疫系統的關鍵組成部分之一,NLR 信號通路能夠識別并響應各種外來病原體、細胞損傷和自身核酸等危險信號, 從而啟動免疫防御反應,維護機體的健康穩定(Van Gorpetal, 2014)。NLR 信號通路主要包括模式識別受體(PRR)、炎癥小體、caspase-1 及相關炎性因子等多個組成部分, 這些部分在免疫防御反應中相互協作, 發揮著重要的作用。

本研究識別了不同病原體入侵長牡蠣鰓組織中的發生可變剪接的差異表達基因, 并對這些基因進行了富集分析。在KEGG 富集分析結果中, 發生可變剪接的差異表達基因被顯著富集在PI3K-AKT 通路、NOD-like 受體信號通路、NF-kappaB 細胞信號傳導通路和Toll 與Imd 信號通路等信號通路上。核轉錄因子κB (NF-κB)是轉錄因子蛋白家族中的一員, 在免疫系統中發揮非常重要的調節作用(Haydenetal,2006)。NF-κB 可以調控細胞因子、生長因子以及效應蛋白酶的表達, 從而參與到免疫應答中。在哺乳動物中, Toll 樣受體(TLRs)直接識別病原微生物的病原體相關分子模式, 通過NF-κB 觸發信號反應, 啟動先天免疫和適應性免疫反應(Takedaetal, 2003)。Toll信號途徑最先被發現于果蠅, 它主要與機體對真菌和革蘭氏陽性細菌的免疫反應有關(Kimetal, 2003;Kanekoetal, 2005)。果蠅Toll 蛋白能與多種PMAPs結合, 從而活化胞內Tube、Pele 等多種胞內信號分子,進而調控免疫相關基因的表達。Imd 信號轉導途徑是一條與機體響應革蘭氏陰性細菌的免疫反應相關的信號轉導途徑, 其在果蠅體內可調控220 種基因的表達, 其中就有革蘭氏抗陰性細菌的多肽(Romeoetal,2008)。最近的研究發現, Imd 信號通路還可以激活JNK (c-Jun N-terminal kinase), 實現細胞骨架蛋白的誘導表達。這表明Imd 信號轉導途徑不僅在天然免疫應答中, 而且也在組織修復中起著關鍵作用。

4 結論

本研究基于RNA-seq 測序數據, 鑒定了長牡蠣病原體誘導過程中發生的可變剪接事件類型及數目,篩選并分析在不同病原體和混合弧菌誘導中差異表達的可變剪接基因的功能, 從而探究可變剪接對長牡蠣天然免疫多樣性和特異性的調控作用。本研究為長牡蠣等無脊椎動物天然免疫多樣性和特異性的深入研究提供了理論依據, 為后續長牡蠣病害防御提供了數據支持。

致謝 本文獲得中國科學院海洋研究所海洋大數據中心的支持, 謹致謝忱。