基于轉錄組測序方斑東風螺LPS刺激的免疫機制研究

陳宗發，洪雨潔，趙明明，王 菁，王忠良

(廣東海洋大學水產學院，廣東湛江 524088)

方斑東風螺(Babyloniaareolata)俗稱花螺，屬軟體動物門腹足綱前鰓亞綱新腹足目蛾螺科東風螺屬[1]，主要分布于環西太平洋。其生產周期短，具有肉質勁道、味道鮮美、營養價值豐富等優點[2]。但是隨著方斑東風螺養殖的快速發展，育苗與養殖過程中病害暴發日趨嚴重[3]，其病害主要以革蘭陰性菌感染為主[4-6]。脂多糖(lipopolysaccharide，LPS)是革蘭陰性菌細胞壁外膜結構的重要組分，可作用于不同生物細胞而表現出多種生物活性[7]。朱璐璐[8]研究發現，飼料中添加LPS能提高凡納濱對蝦對WSSV和副溶血弧菌的抵抗力，周妍英等[9]研究發現，脂多糖可以誘導河南華溪蟹提高血淋巴細胞酸性磷酸酶、堿性磷酸酶、溶菌酶活性和酚氧化酶活性，從而增強其免疫能力，被廣泛應用于水生動物的免疫學研究，揭示LPS刺激對方斑東風螺轉錄組的影響，有助于培育抗病品種。

轉錄組測序(RNA-seq)是一種利用深度測序技術進行轉錄組分析的方法[10]。轉錄組是隨著生物生長發育階段、生理狀態和外界環境的改變而變化的[11]，因此轉錄組測序在水生生物的生長[12]、免疫[13]、適應[14]等方面的研究中得到了廣泛應用。吳勇等[15]通過RNA-seq測序研究馬氏珠母貝(Pinctadafucata)血細胞中存在的補體系統組分，結果表明其補體系統可能通過凝集素途徑或替代途徑發揮生物學功能；張中日[16]通過轉錄組測序分析得出，氰氟草酯使斑馬魚(Brachydanioreriovar)的差異表達基因在藥物等代謝途徑、代謝-細胞色素P450、抗生素的生物合成和外來物質的代謝等途徑明顯富集；張恒澤等[17]通過轉錄組測序，揭示了大黃魚(Larimichthyscrocea)對致病性假單胞菌(LAV)的免疫應答機制，LAV免疫有效激發魚體的細胞免疫應答,有助于快速清除胞內病原。筆者以方斑東風螺為試驗對象，經LPS注射后，利用 RNA-seq 技術獲得轉錄本后，通過DEGs的表達分析和功能富集分析，揭示LPS調控方斑東風螺DEGs的功能和通路。該研究有助于全面深入了解細菌對方斑東風螺發生發展的分子機制，為方斑東風螺的遺傳改良提供理論基礎，對其養殖產業健康持續發展有現實意義。

1 材料與方法

1.1 材料方斑東風螺購自廣東省湛江市某東風螺養殖場，平均體質量為25 g，于實驗室海水桶中暫養7 d(80 L,25 ℃)。將方斑東風螺隨機分為試驗組和對照組，試驗組對方斑東風螺注射100 μL 0.5 mg/mL LPS懸浮液，對照組注射同體積的pH 7.0的磷酸鹽緩沖液(phosphate buffer solution，PBS)，分別于刺激后4、8 h(分別記為LPS-4 h、LPS-8 h組)采集各組足組織樣品，于液氮中速凍，置-80 ℃下保存，用于轉錄組分析。

1.2 轉錄組測序轉錄組文庫的制備及測序由廣州基迪奧生物科技有限公司完成。使用TRIzol法提取總RNA，利用NanoDrop 2000檢測所提RNA的純度，利用Agilent 2100檢測所提RNA的濃度，同處理組10個個體樣本合并后進行建庫測序，用帶有Oligo(dT)的磁珠富集mRNA。加入fragmentation buffer將mRNA打斷成短片段，以mRNA為模板，用六堿基隨機引物(random hexamers)合成第1條cDNA鏈，然后加入緩沖液、dNTPs、RNase H和DNA polymerase I合成第2條cDNA鏈，經過QiaQuick PCR試劑盒純化并加EB緩沖液洗脫后進行末端修復、加poly(A)并連接測序接頭，然后用瓊脂糖凝膠電泳進行片段大小選擇，最后進行PCR擴增，建好的測序文庫用Illumina HiSeqTM2000進行測序，測序讀長為雙端150 bp。

1.3 轉錄組數據分析測序儀產生的原始圖像數據轉化為序列數據(raw reads)，經過初步的過濾后得到質控數據(clean date或clean reads)，再對下機的clean reads去除含有帶接頭、重復、測序質量很低的reads，最后得到高質量的質控數據(clean reads)。通過Trinity[18]軟件進行轉錄組從頭組裝。通過BLAST[19]將Unigene序列比對到蛋白數據庫nr、SwissProt、KEGG和COG/KOG(evalue<1×10-5)，得到與給定Unigene具有最高序列相似性的蛋白，從而得到該Unigene的蛋白功能注釋信息。

1.4 差異表達基因(DEGs)分析采用短reads比對軟件Bowtie2[20]將高質量質控數據與參考基因序列進行對比，根據對比結果，使用RPKM法(Reads Per kb per Million reads)計算Unigene的表達量。利用edgeR[21]進行RNA-seq數據成對樣本之間或組間的差異表達分析，獲得不同時間梯度的DEGs。在篩選過程中，利用FDR與log2FC篩選DEGs，篩選條件為FDR<0.05且|log2FC|>1。得到DEGs后，對DEGs進行GO功能注釋分析及KEGG Pathway富集分析。

2 結果與分析

2.1 測序結果統計及質量評估用Illumina HiSeqTM2000分別對方斑東風螺LPS空白對照組(BLANK)、LPS刺激4 h后取樣組(LPS-4 h)和LPS刺激8 h后取樣組(LPS-8 h)測序，共構建了3個轉錄組文庫，共獲得81 773條unigenes，序列平均長度681 bp，序列長度中位數(N50)為1 035 bp,每組GC含量為46.15%～47.81%，過濾后的clean reads在97.58%～97.74%， Q20均不小于98.49%， Q30均不小于95.40%(表1)，說明測序質量良好，能進行后續的數據分析[22]。從長度分布(圖1)與GC含量等方面對unigenes進行評估，數據顯示測序質量好，可信度高。

表1 轉錄組測序數據評估統計

圖1 方斑東風螺unigenes的長度分布Fig.1 Length distribution of unigene sequence in B.areolata

2.2 Unigene功能注釋使用4大數據庫(KOG、Nr、Swiss-Prot和KEGG)注釋方斑東風螺轉錄組unigenes，分別對應有13 256、21 702、15 335和11 417條unigenes獲得注釋。其中，9 280條unigenes在以上所有數據庫中同時注釋成功，21 938條unigenes至少被1個數據庫注釋。

通過轉錄組測序中的unigenes與GO數據庫的比對，得知在被unigenes注釋的47個功能條目中，包括生物功能(生長、定位等)、細胞組分(細胞器、高分子配合物等)、分子功能(裝訂、運輸等)3個亞類。分析顯示，生物功能含19個功能條目，細胞組分含17個功能條目，分子功能含11個條目(圖2)。

圖2 方斑東風螺GO功能注釋圖Fig.2 GO function annotation of B.areolata

2.3 DEGs分析KEGG數據庫結果顯示，DEGs富集在5個A級KEGG通路,33個B級KEGG通路,232個C級KEGG通路中。其中，5個A級KEGG通路分別為新陳代謝、遺傳信息處理、環境信息處理、細胞進程和有機系統。在232個C級KEGG通路(表2)中，共富集到5 141個基因，其中參與代謝通路的基因數量最多，共有1 455個(28.30%)，參與次生代謝物的生物合成的基因有447個(8.69%)，參與核糖體的基因有384個(7.47%)，參與內吞作用的基因有300個(5.84%)，參與抗生素生物合成的基因有276(5.37%)，參與泛素介導的蛋白水解的基因有261個(5.08%)。

表2 方斑東風螺轉錄組Unigene的KEGG功能分類統計(前10個)

利用edgeR[21]進行RNA-seq數據成對樣本之間或組間的差異表達分析，利用FDR與log2FC篩選DEGs，篩選條件為FDR<0.05且|log2FC|>1。在BLANK vs.LPS-4 h時，顯著DEGs的表達情況為：9 096個顯著DEGs中4 754個上調，4 342個下調；在BLANK vs.LPS-8 h時，顯著DEGs的表達情況為：10 335個顯著DEGs中4 631個上調，5 704個下調；在LPS-4 h vs.LPS-8 h時，表達情況為：8 756個顯著DEGs中3 475個上調，5 281個下調(圖3)。

圖3 樣品間DEGs統計Fig.3 Statistics of DEGs between samples

2.4 DEGs的GO和KEGG注釋在DEGs的GO分類圖(圖4)中，橫坐標為二級GO term，縱坐標為該term的基因數量，紅色代表基因上調，綠色代表基因下調。結果顯示，在BLANK vs.LPS-4 h組中，在生物功能中，細胞過程(cellular process)、代謝過程(metabolic process)和單組織過程(Single organism process)基因的上調與下調較多；在細胞組分中，細胞(Cell)、細胞部分(cell part)和細胞器(organelle)基因的上調與下調較多；而在分子功能中，基因的上調與下調主要集中在綁定(binding)、催化活性(catalytic activity)方面(圖4a)。在BLANK vs.LPS-8 h組中，在生物功能中，細胞過程(cellular process)、代謝過程(metabolic process)和單組織過程(single organism process)基因的上調與下調較多；在細胞組分中，細胞(cell)、細胞部分(cell part)、大分子配合物(macromolecular complex)和細胞器(organelle)的上調與下調較多；而在分子功能中，基因的上調與下調主要集中在綁定(binding)、催化活性(catalytic activity)方面(圖4b)。在LPS-4 h vs.LPS-8 h組中，在生物功能中，細胞過程(cellular process)、代謝過程(metabolic process)和單組織過程(single organism process)基因的上調與下調較多；在細胞組分中，細胞(cell)、細胞部分(cell part)、大分子配合物(macromolecular complex)和細胞器(organelle)的上調與下調較多；而在分子功能中，基因的上調與下調主要集中在綁定(binding)、催化活性(catalytic activity)方面(圖4c)。

圖4 BLANK vs.LPS-4 h(a)、BLANK vs.LPS-8 h(b)和LPS-4 h vs.LPS-8 h(c)的差異基因GO分類圖Fig.4 GO classification maps of differential genes BLANK vs.LPS-4 h (a),BLANK vs.LPS-8 h (b) and LPS-4 h vs.LPS-8 h(c)

對BLANK、LPS刺激4 h和8 h的DEGs進行KEGG通路富集分析，結果顯示，在BLANK vs.LPS-4 h組中，在心肌收縮、神經活性配體受體相互作用和代謝途徑等201個通路中，共富集到506個DEGs(圖5a)；在BLANK vs.LPS-8 h組中，在谷胱甘肽代謝、核糖體和神經活性配體受體相互作用等191個通路中，共富集到551個DEGs(圖5b)；在LPS-4 h vs.LPS-8 h組中，在DNA復制、神經活性配體受體相互作用和鞘糖脂生物合成等173個代謝通路中，共富集到414個DEGs(圖5c)。

圖5 BLANK vs.LPS-4 h(a)、BLANK vs.LPS-8 h(b)和LPS-4 h vs.LPS-8 h(c)的KEGG富集氣泡圖Fig.5 KEGG classification maps of differential genes BLANK vs.LPS-4 h (a),BLANK vs.LPS-8 h (b) and LPS-4 h vs.LPS-8 h(c)

根據DEGs的GO功能分類和KEGG富集分析，篩選到NOD樣受體信號通路、趨化因子信號通路、B細胞受體信號通路、補體和凝血級聯反應、Fc蛋白RI信號通路、自然殺傷細胞介導的細胞毒性、T細胞受體信號通路、toll樣受體信號通路、rig-i樣受體信號通路、血小板激活、造血細胞系通路和抗原的處理和呈遞等免疫信號通路。

3 討論

轉錄組學是研究特定細胞、組織或器官在特定生長發育階段或某種生理狀況下所有轉錄本的學科[10]，基于RNA-seq的轉錄組分析在許多水生生物研究中得到應用，如李喜蓮等[23]通過RNA-seq研了羅氏沼蝦(Macrobrachiumrosenbergii)不同組織的轉錄組差異；王菁等[24]通過轉錄組測序挖掘了方斑東風螺(Babyloniaareolata)的單核苷酸多態性位點。該研究利用RNA-seq技術，對方斑東風螺進行了無參轉錄組測序，得到81 773個unigenes，其中得到注釋的基因數總共為21 938個，注釋率為26.83%，注釋率較低的原因可能是與方斑東風螺相近的物種在公共數據庫中收錄較少，其他無脊椎動物也有注釋率較低的情況，如克氏原螯蝦(Procambarusclarkii)的卵巢、肝胰腺和肌肉組織只有30.1%的基因得到注釋[13]。該研究得到的unigenes與Nr數據庫比對結果顯示，與在公共數據庫中已注冊的方斑東風螺序列匹配的unigenes僅有0.014%，這說明目前公共數據庫中仍缺乏方斑東風螺的基因信息，該研究是對方斑東風螺基因信息庫的有力補充。免疫調控機制一直是貝類研究的熱點，已有研究報道，貝類的免疫受到免疫細胞、體液和化學遞質等的綜合調控[22]。

該試驗中，BLANK vs.LPS-4 h組的DEGs顯著富集至NOD樣受體信號通路，在NOD樣受體信號通路中篩選到半胱天冬酶8(caspase 8; CASP 8)、半胱天冬酶募集域家族9(caspase recruitment domain-containing protein 9; CARD 9)和含有桿狀病毒IAP重復序列的蛋白質2/3(baculoviral IAP repeat-containing protein 2/3; BIRC2_3)。caspase是一類具有天冬氨酸特異性的半胱氨酸蛋白酶，參與所有凋亡通路的級聯反應，是細胞凋亡的核心[25]。Caspase 8是外源性細胞凋亡途徑的關鍵啟動子，caspase 8酶原通過其上的DED與FADD(fasassoclated death domain protein)的DED相結合形成多聚體，將caspase 8酶原招募到細胞質膜上，形成死亡誘導信號復合物DISC(death inducing signaling complex)，在復合物中進而形成有活性的caspase 8，發生caspase級聯反應，啟動細胞凋亡[26]。CARD 9可與開放的Apaf-1結合，進而與大小亞基結合，使caspase 9得到活化，啟動細胞凋亡[27]。BIRC2_3是凋亡相關基因,其編碼的蛋白可抑制細胞凋亡，其中凋亡抑制因子(IAP)是通過抑制半胱天冬酶(caspases)的活性而起到抑制凋亡的作用[28]。有研究者從盤紋羅鮑(Haliotisdiversicolor)中鑒定出了caspase-3-like(saCASP-3lp)基因，并且進一步研究表明saCASP-3lp可能是sacaspase-8上游的調節分子，并參與細胞凋亡過程[29]。在BLANK vs.LPS-4 h組的BIRC2_3呈上調，其通過受體相互作用的絲氨酸/蘇氨酸蛋白激酶2(receptor-interacting serine/threonine-protein kinase 2；RIPK2)激活CASP 8和CARD 9，CASP 8和CARD 9同樣參與了細胞凋亡過程，其中CARD 9激活了p38、JNK，最終激活MAPK通路，促進促炎因子的釋放，參與了細胞的凋亡過程[27]。在BLANK vs.LPS-8 h組，免疫相關DEGs同樣顯著富集至NOD樣受體信號通路，但在NOD樣受體信號通路中除了篩選到CARD 9和BIRC2_3，還篩選到了熱休克蛋白90α(molecular chaperone 90 A；HSP90α)。HSPs是一類進化上高度保守的蛋白，可以增強細胞對應激的反應能力[30]。有研究發現，縊蟶(Sinonovaculaconstricta)[30]、菲律賓蛤仔(Ruditapesphilippinarum)[31]在受到脅迫后HSP90發生顯著變化。在該試驗中，HSP90α在方斑東風螺受到LPS刺激后亦發生顯著變化，說明HSP90α在LPS刺激后的免疫反應中發揮了作用。

在該試驗中，BLANK vs.LPS-8 h組的DEGs顯著富集到趨化因子信號通路。趨化因子(chemokine)是機體重要的免疫調節因子，可協調免疫細胞遷移[32]。范婷婷[33]在短蛸(Octopusocellatus)幼體對鰻弧菌(Vibrioanguillarum)感染應答中發現有3個差異表達基因(Gai、PLCβ、β-arrestin)在趨化因子信號通路富集。該試驗中，有PKA、SOS(son of sevenless)、CDC42(cell division control protein 42)差異表達基因在趨化因子信號通路富集。其中CDC42和Rac1是PAK的主要激活劑，通過調控PKA的表達來抑制肌動蛋白細胞骨架的調節[34]。從而達到清除病原體和促進創傷組織愈合，維持組織細胞的平衡的目的[33]，表明該信號通路可能在方斑東風螺免疫應答過程中發揮重要作用。

在BLANK vs.LPS-4 h組的DEGs除了顯著富集至趨化因子信號通路，還富集至B細胞受體信號通路、補體和凝血級聯反應、Fc蛋白RI信號通路、自然殺傷細胞介導的細胞毒性、T細胞受體信號通路、toll樣受體信號通路、rig-i樣受體信號通路、血小板激活和抗原的處理和呈遞等免疫相關KEGG通路，這些通路均與化學遞質的傳遞有關，在這些通路中篩選出了一些免疫相關的基因，如Ras 相關C3肉毒桿菌毒素底物1(Ras-related C3 botulinum toxin substrate 1；Rac1)、熱休克蛋白70(heat shock 70kDa protein 1/2/6/8；HSPA1s)、α-1-抗胰蛋白酶(alpha-1-antitrypsin；Serpina1b)和蛋白激酶A(protein kinase A[EC:2.7.11.11]；PKA)。其中，Rac1不僅與PAK參與了MARK通路的激活，還通過JNK通路來激活炎癥反應，從而在由免疫球蛋白受體所介導的巨噬細胞吞噬作用中起著關鍵作用[35]。

在LPS-4 h vs.LPS-8 h組的DEGs還富集至抗原的處理和呈遞通路、NOD樣受體信號通路、造血細胞系和Toll樣受體信號通路等免疫信號通路。其中，抗原的處理和呈遞通路篩選到HSP70、HSP90α差異表達基因。其中，HSP90α基因還被富集至NOD樣受體信號通路、PI3K-Akt信號通路等。與試驗組LPS-4 h組相比，HSP90α表達量下降。在NOD樣受體信號通路中，除了篩選到HSP90α基因，還篩選到BIRC2_3基因。與試驗組LPS-4 h組相比，BIRC2_3的表達量差異不顯著，在4 h時表達量最大。在Toll樣受體信號通路中篩選到Rac1。與試驗組LPS-4 h組相比，Rac1表達量下降，在4 h 時表達量最大。在造血細胞系通路中篩選到低親和力免疫球蛋白εFc受體(low affinity immunoglobulin ε Fc receptor；FCER2)。Ren等[36]在鰻弧菌感染下菲律賓蛤仔肝胰腺轉錄組分析中，同樣發現了8個DEGs在造血細胞系通路中。

該試驗注釋到了僅存在于脊椎動物中的特定通路和同源基因，如造血細胞系通路、抗原處理和呈遞通路。這些通路和基因的進化關系和免疫相關功能將是未來幾年的主要研究領域。

4 結論

根據GO功能分類和KEGG信號通路分析篩選到NOD樣受體信號通路、趨化因子信號通路、B細胞受體信號通路、補體和凝血級聯反應、Fc蛋白RI信號通路、自然殺傷細胞介導的細胞毒性、T細胞受體信號通路、toll樣受體信號通路、rig-i樣受體信號通路、血小板激活、造血細胞系通路和抗原的處理和呈遞等免疫信號通路以及大量與方斑東風螺免疫相關的候選基因，如Birc3、Rac1、hsps-1、Serpina1b、Casp8、BIRC2、TPK3、TP02_0244、ACT1、Fcer2、Diap2、Cdc42-、CARD11、actc1和HSP70等。該研究結果豐富了方斑東風螺的基因資源，可為方斑東風螺的免疫研究提供基礎數據。