陈宗发，洪雨洁，赵明明，王 菁，王忠良

(广东海洋大学水产学院，广东湛江 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等。该研究结果丰富了方斑东风螺的基因资源，可为方斑东风螺的免疫研究提供基础数据。