基于转录组测序分析细粒棘球绦虫原头蚴与巨噬细胞RAW264.7免疫互作相关基因

王正荣，马 勋，张艳艳，孟季蒙，薄新文*

(1.省部共建绵羊遗传改良与健康养殖国家重点实验室,石河子 832000;2.新疆农垦科学院畜牧兽医研究所，石河子 832000；3.石河子大学动物科技学院，石河子 832000)

细粒棘球绦虫(Echinococcusgranulosus)是囊型包虫病的主要病原，该病是一种危害严重的人兽共患寄生虫病，是一种被世界卫生组织严重忽视的疾病之一。这种疾病与贫穷和不良的卫生习惯有关，特别是在畜牧业养殖地区[1]。同时，包虫病也是导致我国西部农牧区牧民因病致贫、因病返贫的主要原因之一。目前，该病已在全国20个省份发生，主要分布在西北部分牧区，包括青海、西藏、新疆等地。据估计，我国包虫病患者人数约为100万，受威胁人口近6 600万，患者数量巨大。与现有的全球包虫病流行数据相比，中国包虫病受威胁人口数和患者人数仍居世界首位[1-2]。然而，包虫病病原与宿主的免疫互作的发病机制尚未明确，终末宿主疫苗的研究也相对滞后，这些都是影响该病彻底根除的主要因素和瓶颈问题。

巨噬细胞是一种重要的免疫细胞，是主要的抗原递呈细胞。在机体先天免疫和适应性免疫过程中具有重要的作用。现有的研究表明，巨噬细胞在宿主抗蠕虫感染中具有重要的调控作用[3]。现已有很多有关巨噬细胞在抗棘球绦虫感染中免疫功能研究的文献报道，研究的内容涉及巨噬细胞表达的一些功能分子[4-7]以及非编码的miRNA分子[8-10]和lncRNA分子[11-12]等。同时，也有一些研究棘球绦虫功能分子在调节宿主免疫反应方面的研究报道[13-15]，但截至目前有关棘球绦虫应对宿主免疫反应时，反映虫体整体基因差异表达的研究报道依然较少，致使人们对其缺乏系统的认识。基于此，本研究采用转录组学技术分析了细粒棘球绦虫原头蚴应对巨噬细胞免疫压力时，其基因的差异表达情况，以期系统解析其应对宿主细胞免疫反应时潜在的免疫调控机制，同时为疫苗和药物的研发挖掘新的靶标分子。

1 材料与方法

1.1 试验材料、试剂及主要仪器

在屠宰场收集被细粒棘球绦虫感染的带有包囊的羊肝，带回实验室，无菌条件下抽取含有原头蚴的包囊液，离心收集原头蚴，用加有双抗的PBS洗涤3～5次，最后将收集的原头蚴置于含有10%胎牛血清以及双抗的RPMI 1640培养液中，在37 ℃，5% CO2培养箱中培养备用。于细胞培养瓶中培养RAW 264.7细胞，待细胞生长稳定时将其传代培养至新的培养瓶，当培养瓶中细胞交汇率达到80%时，接种细粒棘球绦虫原头蚴，接种浓度为每毫升2 000个，设PBS为对照组，每组3个生物学重复，共培养24 h后，分别收集培养上清液中原头蚴以及RAW264.7细胞，用PBS洗涤3次，最后将收集的原头蚴保存于-80 ℃冰柜备用，建库测序时采用混池测序。

胎牛血清购自Gibco公司(美国)；RPMI 1640细胞培养液购自Sigma公司(美国)；小鼠单核巨噬细胞RAW 264.7购自ATCC细胞库(中国)；SMARTer Ultra Low RNA Kit for Illumina Sequencing试剂盒购自Clontech公司(美国)；AMPure XP beads试剂盒购自Beckman公司(美国)；Invitrogen Trizol Reagent购自Invitrogen公司(美国)；CW2582 cDNA合成试剂盒购自康为世纪生物技术有限公司(中国)；SYBR PremixExTaqMix购自TaKaRa公司(日本)；实时荧光定量PCR仪Lightcycle 2.0购自罗氏公司(德国)。

1.2 总RNA的提取、转录组文库的建立以及测序

使用SMARTer Ultra Low RNA Kit for Illumina Sequencing试剂盒进行细粒棘球绦虫原头蚴和RAW264.7细胞的裂解和第一条cDNA链的合成，使用AMPure XP beads试剂盒纯化后，用Advantage 2 PCR kit试剂盒扩增第一条cDNA链，并再次纯化，最终得到双链cDNA。Qubit 2.0定量检测，合格后,用Covaris系统超声打断双链cDNA短片段进行末端修复，加A尾并连接测序接头后，纯化，选择片段大小在200 bp左右的文库进行PCR富集得到最终的cDNA文库。先使用Qubit 2.0进行文库的初步定量，稀释文库至2 ng·μL-1，随后用Agilent 2100对文库的insert size进行检测，符合预期后，使用qPCR方法对文库的有效浓度进行准确定量(文库有效浓度>2 nmol·μL-1)，以保证文库质量。文库检验合格后，把不同文库按照有效浓度及目标数据量的需求，进行Hiseq测序。

1.3 差异基因的表达分析

使用Stringtie软件对基因在转录水平上进行表达量分析，皮尔逊相关系数表示样品间基因的表达水平相关性，使用R语言的procmp函数，根据表达量对各样品进行PCA主成分分析。采用DESeq对基因表达进行差异分析，筛选条件为表达差异倍数|log2(Fold Change)|>1，显著性Pvalue<0.05。

使用R语言ggplots 2软件包绘制差异表达基因的火山图，Pheatmap软件包根据同一基因在不同样品中的表达水平和同一样品中不同基因的表达模式进行聚类，采用Euclidan方法计算距离，使用层次聚类最长距离法进行聚类。

1.4 差异表达基因的GO分析

利用GO(gene ontology)对差异表达基因进行功能聚类分析。GO分为分子功能(molecular function)、生物过程(biological process)和细胞组成(cellular component)3个部分。差异表达的基因可以通过ID对应或者序列注释的方法找到与之对应的GO编号，而GO编号可用于对应到TERM，即功能类别或者细胞定位。如P-value≤0.05，则表示显著性富集。

1.5 差异表达基因的KEGG分析

采用KEGG(kyoto encyclopedie of genes and genomes)数据库对差异表达的基因进行了信号通路分析，通过Pathway显著性富集能确定特异基因参与的最主要生化代谢途径和信号转导途径。具体方法以KEGG Pathway为单位，应用超几何检验，找出与整个基因组背景相比，在特异基因中显著性富集的Pathway，最终将FDR≤0.05的Pathway定义为在差异表达基因中显著富集的Pathway。

1.6 部分差异表达基因的转录水平分析

将无菌采集的细粒棘球绦虫原头蚴保存于液氮中，备用。总RNA提取采用Invitrogen Trizol Reagent试剂，提取方法参照试剂说明书步骤进行。反转录采用康为世纪公司10 μL反转录体系的CW2582 cDNA合成试剂盒。

以上述cDNA为模板，以GAPDH为内参基因，对随机筛选的差异表达基因(引物见表1)进行qRT-PCR扩增，反应体系：SYBR PremixExTaqMix(2×) 10 μL，上下游引物(10 mmol·L-1)各0.4 μL，cDNA 1 μL，加ddH2O 至20 μL。反应条件：95 ℃预变性30 s；95 ℃ 5 s，60 ℃ 20 s，共40个循环。每个基因做3个生物学重复，采用相对比较ΔCt(Qr=2-ΔΔCt)法计算目的基因的相对表达量。

表1 实时定量PCR试验相关引物

2 结 果

2.1 细粒棘球绦虫原头蚴差异基因表达分析

当细粒棘球绦虫原头蚴与巨噬细胞RAW264.7共培养24 h后，与0 h相比，共培养24 h后的原头蚴共有435个基因的表达出现显著变化，其中，227个基因呈上调表达，208个基因呈下调表达。对差异表达基因进行火山图以及前20位(top20)聚类热图可视化分析，本研究将-lg(FDR)>1的差异表达定义为显著差异表达基因(图1)。

图1 原头蚴与RAW264.7共培养24 h后差异表达基因的火山图(A)、top20热图(B)分析

其中,上调表达前10的基因包括Novel00223、EgrG_000549200、EgrG_001026200(NHRF3)、EgrG_000207400、EgrG_000381700、Novel01656、Novel00517、Novel00271(EGR_00497)、EgrG_000464200和Novel01519(EG95 protein precursor)。下调表达前10的基因包括EgrG_000155900(ALDR)、EgrG_001072700(DPYD)、EgrG_002026800(SYNC)、EgrG_001085050(HSP70)、EgrG_000654600(CATL)、EgrG_000355900、Novel00327(EF-hand domain-containing protein)、EgrG_002320300、pathogen_EgG_scaffold_0325_ncRNA_4499-4601(LORF2)和EgrG_002059100(表2)。

表2 感染后24 h Top10差异表达基因

此外，在上调表达的基因中还有EgrG_000320800(HSP10)、EgrG_000515900(AgB)、EgrG_000550000(FABP1)、EgrG_000253300(囊泡运输蛋白SC22B)、EgrG_000460300(钙结合蛋白SM16)和EgrG_000715700(PDCD6)等。下调表达的基因中还包括EgrG_001018500(跨膜蛋白TM144)、EgrG_000802400(内固醇类受体)、Novel01542(MAPK7)、EgrG_000317500(丙氨酸氨基转移酶2，ALAT2)、EgrG_000589100(谷氨酸脱氢酶2，GLUD2)和EgrG_000644000(谷氨酸合酶，NADH)等。

2.2 细粒棘球绦虫原头蚴差异表达基因的GO分析

GO分析的结果显示，当细粒棘球绦虫原头蚴与巨噬细胞RAW264.7共培养24 h后，其差异表达基因的功能主要富集在生物学过程(BP)的血红素转运、铁配位实体运输、辅因子输送、α-氨基酸生物合成过程以及细胞形态发生等。细胞组分(CC)的细胞外间质、核糖核酸酶MRP复合物、核糖核酸酶P复合物以及胞内核糖核酸复合体等。分子功能(MF)的丝氨酸型内肽酶抑制剂活性、α-半乳糖苷酶活性、辅因子结合以及磷脂酰肌醇乙酰氨基葡萄糖转移酶活性等过程(图2，表3)。

图2 与RAW264.7共培养24 h后原头蚴中差异表达基因的GO分析

2.3 细粒棘球绦虫原头蚴差异表达基因的KEGG分析

KEGG分析的结果显示，当细粒棘球绦虫原头蚴与巨噬细胞RAW264.7共培养24 h后，原头蚴中差异表达的基因主要参与剪接体、内吞作用、内质网的蛋白质处理、吞噬体、MAPK信号通路以及钙信号通路等。其中，上调表达的基因主要参与剪接体、内质网的蛋白质处理、核糖体、吞噬体以及MAPK等信号通路。下调表达的基因主要参与丙氨酸、天冬氨酸和谷氨酸代谢、钙信号通路、氨基酸生物合成以及内吞作用等信号通路(图3，表4)。

图3 与RAW264.7共培养24 h后原头蚴中差异表达基因的KEGG分析

2.4 原头蚴差异表达基因中各种酶相关基因的表达谱

分析发现，在435个差异表达的基因中，有79个各种酶相关的基因，其中，上调表达的26个，主要包括乙酰氨基葡萄糖转移酶、神经鞘脂脱饱和酶、特异性二酰甘油脂肪酶、线粒体输入性内膜转位酶、谷氨酸脱氢酶、谷氨酸半胱氨酸连接酶、胞苷脱氨酶、L-天冬酰胺酶、组蛋白赖氨酸甲基转移酶、蛋白质糖基转移酶、信号肽酶、甘露糖-1-磷酸鸟苷转移酶、甘油激酶、脱氧胞嘧啶合酶、葡萄糖胺6-磷酸N-乙酰基转移酶、精氨酸转氨酶、蛋白转移酶、甲戊二酸单酰辅酶A合酶、磷脂促翻转酶以及亚精胺合成酶。同时还包括3个抑制剂相关的基因，分别为Kunitz型丝氨酸蛋白酶抑制剂6、Kunitz型丝氨酸蛋白酶抑制剂Bt-KTI以及蛋白酶抑制剂I2(表3)。下调表达有53个，主要包括天冬酰胺-tRNA连接酶、二氢嘧啶脱氢酶、醛糖还原酶、腺苷酸琥珀酸裂解酶、谷氨酸脱氢酶2、谷氨酸合酶(NADH)、鸟氨酸转氨酶、L-苏氨酸3-脱氢酶、L-乳酸脱氢酶A、尿苷磷酸化酶1等，图4显示了各种酶相关的基因差异表达top20的分析结果。

图4 与RAW264.7共培养24 h后原头蚴中差异表达酶相关基因top20的表达谱分析

表3 与巨噬细胞RAW264.7共培养24 h后原头蚴中上调表达的各种酶相关基因

2.5 巨噬细胞RAW264.7基因差异表达分析

当巨噬细胞RAW264.7与细粒棘球绦虫原头蚴共培养24 h时，和PBS对照组相比，原头蚴处理组的RAW264.7细胞共有3 745个基因的表达出现显著变化，其中，1 159个基因出现上调表达，2 586个基因出现下调表达(图5)。

图5 RAW264.7与原头蚴共培养24 h后差异表达基因的火山图(A)、top20热图(B)分析

其中，表达上调前10的基因包括ENSMUSG0011463(Cpb1)、ENSMUSG0058579(Cela2a)、ENSMUSG0029882(2210010C04Rik)、ENSMUSG0036938(Try5)、ENSMUSG0035896 (Rnase1)、ENSMUSG0062478(Ctrc)、ENSMUSG0084174(Sycn)、ENSMUSG0031896(Ctrl)、ENSMUSG0031957(Ctrb1)和ENSMUSG0074268(Amy2a5)；下调表达前10的基因包括ENSMUSG0022678(Nde1)、ENSMUSG0067367(Lyar)、ENSMUSG00110686、ENSMUSG0057657、ENSMUSG0041506(Rrp9)、ENSMUSG0031527(Eri1)、ENSMUSG0018040(Rrp7a)、ENSMUSG0069682、ENSMUSG0046330(Rpl37a)和ENSMUSG0078965(表4)。

表4 共培养24 h后 RAW264.7 Top10差异表达基因

2.6 巨噬细胞RAW264.7差异基因表达的GO分析

当巨噬细胞RAW264.7与细粒棘球绦虫原头蚴共培养24 h时，与PBS对照组相比，GO分析结果显示，差异表达基因主要富集在GO功能聚类的生物过程(BP)中的代谢过程、初级代谢过程、有机物代谢过程以及细胞代谢过程等。细胞组分(CC)中的细胞内组分、细胞内、细胞器以及膜结构细胞器等(图6)。

图6 与原头蚴共培养24 h后RAW264.7中差异表达基因的GO分析

2.7 巨噬细胞RAW264.7差异基因表达的KEGG分析

当巨噬细胞RAW264.7与细粒棘球绦虫原头蚴共培养24 h时，与PBS对照组相比，KEGG信号通路的分析结果表明，差异表达的基因主要参与代谢通路、核糖体通路、剪接体通路、RNA转运以及泛素介导的蛋白水解等通路(图7)。

图7 与原头蚴共培养24 h后RAW264.7中差异表达基因的KEGG分析

2.8 差异表达基因的qRT-PCR验证

为了验证从转录组测序数据中获得差异基因结果的可靠性，随机选取10个差异表达的基因进行qRT-PCR验证，其中，包括6个上调表达基因，以及4个下调表达基因，结果显示，qRT-PCR验证结果与RNA-seq测序结果一致，进一步说明转录组测序结果可信(图8)。

图8 差异表达基因的qRT-PCR验证分析

3 讨 论

原头蚴是细粒棘球绦虫的重要发育阶段，具有双向发育的特征，即在中间宿主体内可以经无性繁殖发育为新的包囊，在终末宿主体内可以经有性繁殖发育为成虫[16]。原头蚴也是中间宿主形成二次感染的重要原因，同时还是细粒棘球绦虫感染终末宿主，在其肠道定植发育为成虫的最初形态。可以看出，原头蚴不仅是细粒棘球绦虫生长发育的关键环节，同时也是感染宿主的重要阶段。因此研究细粒棘球绦虫原头蚴生长发育关键分子以及原头蚴与终末宿主免疫互作关键分子均具有重要意义，研究结果将不仅揭示其独特的发育现象，同时还可能挖掘出潜在疫苗候选靶标，进而为疫苗的研发奠定基础。基于以上原因，本研究通过转录组学技术解析了细粒棘球绦虫原头蚴与重要免疫细胞-巨噬细胞的免疫互作关键分子，以期初步阐明其免疫互作机制，同时发掘疫苗候选抗原。

本研究的结果表明，当细粒棘球绦虫原头蚴与巨噬细胞RAW264.7共培养24 h后，和0 h相比，共培养24 h后的原头蚴共有435个基因的表达出现显著变化，其中227个基因出现上调表达，208个基因出现下调表达。由此可见，当原头蚴在面对巨噬细胞RAW264.7免疫压力时，其自身基因表达会发生变化，包括与抗逆反应相关的热激蛋白HSP70、HSP10，与免疫调节相关的AgB、FABP1和Kunitz型丝氨酸蛋白酶抑制剂，还有与胞外囊泡运输相关的蛋白SC22B，这些基因的表达都出现显著上调，由此推测其可能在原头蚴抗巨噬细胞RAW264.7免疫压力过程中起到重要的调控作用，但其具体的调控机制还有待进一步研究。反观巨噬细胞RAW264.7的基因表达情况，当巨噬细胞RAW264.7与细粒棘球绦虫原头蚴共培养24 h后，和PBS对照组相比，原头蚴处理组的RAW264.7细胞共有3 745个基因的表达出现显著变化，其中，1 159个基因出现上调表达，2 586个基因出现下调表达。在这些差异表达的基因中，作者发现大量与天然免疫相关的基因的表达也发生了变化，诸如C型凝集素受体Reg1、Reg2、Reg3b、Reg3d、Reg3a和Reg3g的表达出现上调。病原模式识别受体下游相关效应分子Fam132a、Cbln2、和Tnfrsf12a的表达也出现上调，其中，Fam132a和Cbln2属于肿瘤坏死因子样的结构域，而Tnfrsf12a属于肿瘤坏死因子受体12，都属于炎型细胞因子调节因子。并且在表达上调前10的分子中有6个属于丝氨酸蛋白酶家族，包括Cela2a、2210010C04Rik、Try5、Ctrc、Ctrl和Ctrb1[7]。提示这些基因可能在巨噬细胞RAW264.7抗原头蚴感染过程中发挥重要功能。

在这些差异表达的基因中，AgB基因的表达出现上调，AgB是棘球绦虫中研究最为广泛的抗原之一[17]，其在细粒棘球绦虫包囊液中普遍存在，是目前人类包虫病诊断的主要抗原之一。已有的研究证实，细粒棘球绦虫分泌的AgB可以通过中性粒细胞分泌的弹性蛋白酶干扰中性粒细胞的活性，进而有利于寄生虫逃避宿主的免疫反应。另外，AgB可以通过调节巨噬细胞的活性和抑制有效细胞因子的产生来调节宿主的免疫系统[18]。此外，本研究还发现FABP1基因的表达同样出现上调。现有的研究结果表明，FABP对于寄生的扁形动物具有重要的意义，因为这些寄生虫不能合成它们自身所需的大部分脂质，特别是长链脂肪酸和胆固醇[19-20]。因此，这种脂质必须从其寄生宿主体内获得，然后通过载体蛋白质运送到寄生虫体内的特定部位。FABP蛋白可能参与了其细胞外脂质的获取，以及传递到宿主细胞的过程。同时FABP蛋白家族也是极具潜力的药物传递以及疫苗候选靶标[21-24]。研究者在棘球蚴囊液和原头节分泌物中发现了FABP1分子[24-25]。细粒棘球绦虫的FABP1被认为是其亚家族的重要成员之一[26]，其成员参与脂质氧化过程。研究表明，细粒棘球绦虫的FABP1通过一种机制与膜相互作用，这种机制涉及到与膜直接接触以它们的物质交换[27]。由此推断，以上分子可能在细粒棘球绦虫应对巨噬细胞免疫反应时起到重要免疫调控作用，以创造有利于虫体寄生的免疫微环境。

虫体表达的各种酶相关的蛋白以及抑制剂，不仅在虫体自身的代谢中起重要的作用，同时可以降解对其有害的宿主蛋白以及抑制宿主来源的各种酶，保护虫体自身不受危害[28]。为此，本研究系统分析了细粒棘球绦虫原头蚴在应对巨噬细胞免疫压力时，其各种酶相关基因以及抑制剂相关分子的表达情况。结果发现，有79个酶相关的基因其表达出现显著差异，其中，包括3个上调表达的蛋白酶抑制剂，分别为Kunitz型丝氨酸蛋白酶抑制剂6、Kunitz型丝氨酸蛋白酶抑制剂Bt-KTI以及蛋白酶抑制剂I2。有学者在细粒棘球绦虫中分离鉴定了2个Kunitz型蛋白酶抑制剂，其中，一种为Kunitz型丝氨酸蛋白酶抑制剂，是一种强有力的糜蛋白酶和中性粒细胞弹性蛋白酶抑制剂，在局部炎症中结合钙以及减少中性粒细胞浸润。同时也可能通过抑制中性粒细胞弹性蛋白酶和组织蛋白酶G，参与宿主免疫逃避。而另一个是一种强有力的胰蛋白酶抑制剂。作为哺乳动物肠道蛋白酶的强有力抑制剂，细粒棘球绦虫Kunitz型蛋白可能在防止蛋白水解酶攻击方面发挥重要的保护作用，从而确保细粒棘球绦虫在其哺乳动物宿主体内的存活[29]。近期，有研究发现，细粒棘球绦虫的Kunitz型蛋白酶抑制剂可以抑制巨噬细胞炎型细胞因子的产生，并且抑制巨噬细胞的增殖[15, 30]。肝片吸虫的排泄-分泌产物中存在两种Kunitz型丝氨酸蛋白酶抑制剂，这些抑制剂可以保护寄生虫在入侵时免受宿主蛋白酶，如糜蛋白酶等有害蛋白的水解作用。同时研究还发现此类丝氨酸蛋白酶抑制剂可以通过干扰激肽释放酶-激肽系统，进而调节宿主的炎症和通透性反应[31]。由此作者推测，这3种蛋白酶抑制剂可能在细粒棘球绦虫应对巨噬细胞免疫反应时起重要的作用，一方面，抑制宿主细胞释放的蛋白酶活性，保护虫体在寄生过程中免受宿主蛋白酶水解作用；另一方面，抑制宿主免疫细胞分泌炎性细胞因子，调节宿主免疫反应，达到寄生的目的。

4 结 论

本研究系统解析了细粒棘球绦虫原头蚴应对巨噬细胞RAW264.7免疫压力时其基因的差异表达谱规律及其巨噬细胞RAW264.7面对虫体抗原刺激时其免疫相关基因的表达情况，筛选了与细粒棘球绦虫抗宿主免疫反应相关分子。研究结果显示，当细粒棘球绦虫原头蚴面对宿主巨噬细胞免疫压力时，原头蚴中的AgB、FABP1、Kunitz型丝氨酸蛋白酶抑制剂等具有免疫调节作用的分子表达明显上调，进而有利于虫体在宿主的寄生以及免疫逃避。这些数据将为进一步揭示细粒棘球绦虫抗宿主免疫反应以适应寄生的机理提供理论依据，同时为包虫病防治提供新的疫苗和药物靶点奠定基础。