基于生物信息学探讨弓形虫ROP16Ⅰ对A549细胞表达谱的影响

苏雅静，乔 霞，刘学雷，杨宁爱，贾 伟,2，赵志军,2

刚地弓形虫(Toxoplasmagondii，T.gondii)于1908年由Nicolle等[1]发现，是一种寄生于真核细胞内的顶复门原虫，可在人和动物体内建立持久的感染[2]。最近研究显示，AIDS合并弓形虫病的情况愈发严峻，以撒哈拉以南的非洲地区最为显著，AIDS合并弓形虫感染可占AIDS感染的87.1%[3-4]。弓形虫除了会造成中枢神经系统感染外，亦可造成其他器官的感染，尤其是肺，但目前关于肺弓形虫病的报道较少[5]。

弓形虫在侵入宿主细胞、增殖及致病过程中，棒状体蛋白 (Rhoptry, ROP)等效应分子起到重要作用，其可在多点操纵宿主的抵抗机制、调节促炎反应，并对宿主基因表达和信号通路进行广泛的修饰[4]。ROP16属于ROP2家族中的丝氨酸/苏氨酸激酶，可直接定位于宿主细胞核，并参与宿主细胞内信号转录活化因子(STAT)相关的磷酸化。ROP16Ⅰ可调节宿主基因转录和蛋白质表达。研究发现，ROP16Ⅰ/Ⅲ可激活并维持STAT3/6的激活反应，而ROP16Ⅱ则无此功能。此外，在感染早期ROP16Ⅰ/Ⅲ可诱导宿主巨噬细胞向M2型极化，Th2 大量表达能够促使ariginase I 酶活性升高，L-arginine 代谢途径向着生成大量多胺类物质进行，促使弓形虫急性感染或复发感染[6-7]。但目前关于ROP16Ⅰ所调控的宿主基因和分子仍不完全清楚，因此从基因水平研究ROP16Ⅰ对弓形虫病的防治具有重要意义。本研究拟构建ROP16Ⅰ真核重组表达载体，通过对ROP16Ⅰ基因转染A549细胞株后的差异表达基因进行筛选，并联合生物信息学工具和数据库进行分析，从基因和分子水平探讨ROP16Ⅰ的致病机制。

1 材料与方法

1.1材料与试剂 大肠杆菌DH5α、真核表达质粒pEGFP-N1及弓形虫RH株为本实验室保存；Trizol(Invitrogen公司，美国)；PrimeScriptTMRT reagent Kit with gDNA Eraser(Perfect Real Time)、SYBR Premix Ex TaqTM荧光定量PCR试剂盒均购自大连Takala公司；TransLipid Transfection Reagent、Taq DNA聚合酶、T4连接酶均购自北京全式金公司；EcoRⅠ和SalⅠ限制性内切酶购自美国NEB公司；A549人肺腺癌细胞购自上海细胞库；细胞培养相关试剂购自Gibco公司；ICR小鼠购自宁夏医科大学动物中心；基因芯片检测由上海康成生物有限公司完成。

1.2 方 法

1.2.1弓形虫速殖子总RNA的提取 从液氮罐取出RH虫株，复苏后腹腔接种ICR小鼠。观察小鼠状态，出现背毛逆立、闭目及腹部膨大等典型症状时，折颈处死，收集腹水。传代3次后，抽取腹水，800 g离心去掉上清，PBS重悬后2000 g离心，收集沉淀即为RH虫株。采用 Trizol法提取总RNA，Nano Drop 2000测定总RNA的纯度和浓度。

1.2.2构建ROP16Ⅰ基因真核表达载体 根据GeneBank中ROP16Ⅰ基因全长编码序列设计引物，引物两端同时包含EcoR-Ⅰ，Sal-Ⅰ酶切位点，上游序列：5′-CCGGAATTCTGATGAAAGTGACCACGAAAGGG-3′，下游序列：5′- ACGGTCGACCATCCGATGTGAAGAAAGTTC-3′。建立40 μL PCR反应体系，退火温度为55 ℃。将PCR产物双酶切，回收目的片段。将其连接至真核表达载体pEGFP-N1上，连接产物转化至DH5α感受态细胞后提取质粒，双酶切鉴定并测序[8]。

1.2.3数据处理与差异基因筛选 应用Agilent Feature Extraction软件对表达谱结果进行处理，提取原始数据。而后，采用GeneSpring软件对原始数据进行Quantile标准化和随后的处理。差异表达基因(differential expression gene, DEG)通过Fold Change筛选，标准为差异倍数Fold change(FC)≥ 2.0，负值表示下调。

1.2.4文献挖掘已知候选基因 Genecards是一个收集并展示人类基因及其产物和相关疾病等综合信息的知识平台[9]，本研究以“Toxoplasma”作为搜索词，进入Genecards(http://www.genecards.org/)搜索已知弓形虫疾病基因，也采用OMIM(http://omim.org/)在线工具搜索弓形虫疾病相关基因。

1.2.5ToppGene筛选候选基因 ToppGene(https://toppgene.cchmc.org/)是一个基于GO数据库、KEGG数据库、BIOCYC数据库及REACTOME 信号通路等数据库的基因列表功能富集分析的一站式门户网站[10]。将Genecards及OMIM两个数据库搜索到的共有基因作为“Training Gene Set”，而以所获得的差异表达基因作为“Test gene set”，利用ToppGene进行候选基因的筛选，并对差异表达基因进行GO 分析及KEGG分析。

1.2.6蛋白互作网络 STRING(https://string-db.org/)蛋白质相互作用数据库，主要通过高通量实验数据、共表达数据、基因组上的关联及之前文献的累积挖掘整合收录了已知蛋白间和预测蛋白质间的相互作用关系[2]。将得到的差异表达基因上传至STRING数据库，通过在线分析差异表达基因所编码的蛋白间的相互作用图，找出处在关键节点的蛋白质。而后，将ToppGene筛选得到的弓形虫疾病候选基因也上传至STRING在线工具，绘制出相应的蛋白互作网络图。

2 结 果

2.1差异表达基因的筛选 经分析后发现，ROP16Ⅰ基因转染A549细胞后差异性表达的基因共有483个，其中表达上调的基因有252个，表达下调的基因有231个，上调和下调幅度最大的前10个基因见表1。在上调表达的差异基因分布中，BC015119表达量最大(32.71)，其次是 LOC646999(16.83)；而在下调表达的差异基因分布中，OR51G2(-117.86)下调幅度最大。

2.2Genecards和OMIM挖掘已知弓形虫疾病相关基因 在Genecards中搜索已知“Toxoplasma” 疾病基因，共获得191条gene card记录；通过OMIM查询获得39条疾病相关基因。对这两种方法获得的基因进行交集分析，共获得已知弓形虫疾病相关基因9个，见表2。

2.3筛选Ι型弓形虫疾病候选基因 采用ToppGene候选基因筛选策略，对候选基因进行优化排序，表3显示的是筛选出的前30个Ι型弓形虫疾病候选基因。通过对ToppGene和ToppNet获得的候选基因进行交集分析，发现共有14个基因同时被筛出。

表1 FC值(上调/下调)表达量最大的部分基因
Tab.1 Representative list of(up-regulated/down-regulated)genes

RankUp regulated gene SymbolFCRankDown regulated gene SymbolFC1BC01511932.711OR51G2-117.862LOC64699916.832RGS3-6.913LOC1001298908.223APCDD1-5.524CMPK27.754DIO3-5.335FGF56.805LIPG-4.906OAS26.096C18orf26-4.787BATF25.927C1orf180-4.698OASL5.908OR51F2-4.489IFIT25.899GIMAP1-4.4110IFIT35.7710SCARNA4-4.26

表2 通过Genecards和OMIM获得的已知弓形虫疾病相关基因
Tab.2 Toxoplasma related genes using Genecards and OMIM

RankSymbolGene ID1STAT667782MYD8846153CCR112304TNF71245IL435656IL12B35937CXCL1036278ALOX52409IDO13620

表3 前30个ToppGene筛选的Ι型弓形虫疾病候选基因
Tab.3 Top 30 candidate genes of ToppGene forT.gondiitype Ι

ToppGeneToppNetRankGene symbolAverage score(×10-1)Overall P-value(×10-3)RankGene symbolInteractant countScore(×10-4)1STAT110.0001.6171CXCR3520.3502TLR39.9972.4362CCL54113.3703IFIH19.9573.3073STAT11949.7084DDX589.9713.8084CCL899.5195STAT29.6203.8565ESR17668.8466HLA-B9.8544.0456CXCL1146.8927IFNB110.0004.4507HDAC53564.2718CCL59.9985.2828AMOT1014.0799CXCR39.9955.6469NMI693.804表3(续)ToppGeneToppNetRankGene symbolAverage score(×10-1)Overall P-value(×10-3)RankGene symbolInteractant countScore(×10-4)10FLT39.3086.19710PML2493.70711TREX19.5246.62511IFI163573.65012FGFR18.8637.51512IRF7483.43213CXCL119.9617.54513STAT2493.36614EIF2AK29.9687.54514TLR3323.18515CASP19.9527.54515LMO22083.15916IL229.9457.66916HTT2553.09817OAS19.6018.31617PARP14132.88618FGFR28.3178.35118CD93352.65319ESR18.3089.94919BST2122.40320IRF77.4989.95520ISG152142.27221CCL89.71111.28021PIAS11312.19422ANGPT28.24411.60022PCDH7401.87923TNFSF109.46612.05023FGFR11041.55524TAP17.82912.36024FGFR2921.43425TAP27.82812.36025EIF2AK21061.41526IFNA169.98914.45026RGS3581.27127OSMR9.16715.10027PLSCR11371.24328OAS29.51916.48028LAMP3641.03529BMP78.69116.57029IFIT3690.99930PIAS19.64917.33030HLA-B780.990

2.4GO功能聚类分析 利用ToppGene对差异表达基因进行GO分析。发现，ToppGene共获得95条有统计学意义的生物学过程，其中MF(Molecular Function，MF)有37条，BP(Biological Process，BP)有50条，CC(Cellular Component，CC)有8条。GO 富集结果主要表现在BP和MF两大类。表4显示的是前10条GO功能注释，主要富集在免疫应答、炎症反应、防御反应、免疫系统进程、白细胞游走控制、趋化因子合成调节和细胞因子产生的调控等方面(表4)。

表4 前10条差异表达基因的GO功能注释
Tab. 4 Top 10 enriched GO terms for DEGs

RankCategoryGO IDGO nameGene numberP(×10-9)1BPGO:0006954inflammatory response80.3782BPGO:0002684positive regulation of immune system process80.4753BPGO:0032103positive regulation of response to external stimulus62.2924BPGO:0050777negative regulation of immune response52.7495BPGO:0002685regulation of leukocyte migration58.7036BPGO:0045080positive regulation of chemokine biosynthetic process39.3507BPGO:0001819positive regulation of cytokine production69.7908BPGO:0002683negative regulation of immune system process610.0709BPGO:0034097response to cytokine710.88010BPGO:0002682regulation of immune system process815.020

2.5KEGG通路分析 利用ToppGene对差异表达基因进行KEGG分析，共获得50条有统计学意义的信号通路，表5显示的是前10条信号通路。主要涉及白介素IL-10, IL-4, IL-13信号通路，细胞因子-细胞因子受体相互作用通路、细胞因子网络及Toll样受体信号通路等(表5)。

表5 前10条富集的信号通路
Tab. 5 Top 10 enriched KEGG pathways for DEGs

RankIDPathway namenumberGene nameP(×10-11)11269318Signaling by Interleukins8ALOX5,TNF,MYD88,IL12B,STAT6,CXCL10,IL4,CCR19.01721470924Interleukin-10 signaling5TNF,IL12B,CXCL10,IL4,CCR19.53431269310Cytokine Signaling in Immune system8ALOX5,TNF,MYD88,IL12B,STAT6,CXCL10,IL4,CCR1163.7004194384African trypanosomiasis4TNF,IDO1,MYD88,IL12B652.80051470923Interleukin-4 and 13 signaling5ALOX5,TNF,IL12B,STAT6,IL4721.2006842771Inflammatory bowel disease 4TNF,IL12B,STAT6,IL48360.0007144181Leishmaniasis4TNF,MYD88,IL12B,IL413400.0008M17406Cytokine Network3TNF,IL12B,IL434520.000983051Cytokine-cytokine receptor in-teraction5TNF,IL12B,CXCL10,IL4,CCR154250.0001083076Toll-like receptor signaling pathway4TNF,MYD88,IL12B,CXCL1056070.000

2.6蛋白互作网络 通过STRING 在线分析工具对差异表达基因进行蛋白互作分析，结果如图1所示。发现， STAT1、STAT2、HLA-B、CXCR3、CCL5、TLR3 等存在于蛋白相互作用网络的中心节点，删除这些关键节点后，蛋白互作网络结构涣散。而有趣的是，ToppGene 所筛选得到的14个候选基因所调控的蛋白恰好处于网络的中心节点处。随后我们将这14个候选基因也上传至STRING 在线工具，绘制出相应的蛋白互作图，结果如图2所示。

图2 Toppgene 筛选的候选基因蛋白互作图Fig. 2 PPI network of the candidate genes screened from the ToppGene

3 讨 论

弓形虫入侵宿主细胞过程复杂，虫体蛋白和宿主蛋白等多种因子相互作用并参与该过程。ROP16是弓形虫入侵过程中的重要毒力因子，可参与宿主细胞内 STAT3/6 相关的磷酸化，并抑制宿主细胞IL-12的生成，干扰宿主细胞的增殖、分化和凋亡等功能[11-12]，这是弓形虫发挥其致病性的一个重要途径。Saeij等[13]研究发现，I 型弓形虫入侵宿主细胞后可诱导 STAT3/6 发生强烈的磷酸化，而 II 型虫株则无此能力。Butcher等[14]研究也发现弓形虫基因型是决定 STAT 信号通路能否被激活的重要因素。ROP16I感染小鼠巨噬细胞后能够持续诱导 STAT3 的活化，并驱动巨噬细胞向 M2 极化，同时抑制IL-12生成。而 ROP16II虫株则不能持续激活STAT3的磷酸化，但可导致下游IL-12的大量合成。目前关于弓形虫 ROP16I入核后的具体分子机制仍不清楚，涉及多个基因、多条通路的相互作用和相互影响。

图1 差异表达基因蛋白互作图Fig.1 PPI network of the screened DEGs

生物信息学在表达谱的基础上能够多层次的了解基因表达情况，对基因功能富集分析、通路富集分析及网络调控分析的探讨具有重要作用[15]。本研究运用 ToppGene 在线分析工具，筛选了Ι 型弓形虫疾病候选基因；同时利用 STRING 分析工具构建了ROP16I分子调控网络。发现弓形虫 ROP16Ⅰ在人肺腺癌A549细胞过表达后，共有483个基因表达差异显著，其中上调基因252个，下调基因231个。在上调表达差异的基因分布中，BC015119表达量最大，其次是LOC646999；而在下调表达差异的基因分布中，OR51G2变化量最大，下调倍数达-117.86。LOC646999属Akirin1 家族，而 Akirins 是天然免疫反应信号通路所必需的核内因子。研究发现，Akirin 在NF-κB 转录调控中起着重要作用，而 NF-κB 与发育、免疫应答及细胞增生和细胞凋亡等多种生物过程密切相关[16]。OR51G2 是嗅觉受体家族中的一员，参与G蛋白偶联受体蛋白介导的气味信号转导途径，从而引发对气味的感知。而BC015119属于未知基因，在转染ROP16I基因后BC015119上调显著，其生物学功能有待进一步研究。

弓形虫感染宿主后通过复杂的细胞因子网络进行免疫调节[17]。本研究通过 ToppGene 生物信息学工具最终筛选得到CXCR3、CCL5、STAT1、CCL8、ESR1、CXCL11、IRF7、STAT2、TLR3、PIAS1、FGFR1、FGFR2、EIF2AK2及HLA-B 等14个 Ι 型弓形虫相关的疾病候选基因。GO 分析和 KEGG 分析发现，这些候选基因与白介素IL-10、IL-4、IL-13信号通路，细胞因子-细胞因子受体相互作用通路、细胞因子网络及Toll 样受体信号通路相关。在蛋白-蛋白互作网络中，我们发现这14个候选基因恰好处于网络的最核心位置。目前，关于CCL8、CXCL11、PIAS1及EIF2AK2 在弓形虫疾病方面的研究报道很少。CCL8 作为趋化因子的一员，在机体炎症反应和免疫应答等方面发挥重要作用。有研究显示，寄生虫感染可引起CCL8表达的增加，CCL8 与 CXCL12通过协调发挥免疫抑制作用，导致持续感染。通过降低 CCL8水平，能够显著抑制疾病的发生发展[18-19]。CXCL11正常情况下表达量较少。CXCL11-CXCR3是免疫炎症反应的一个重要环节，参与炎症、感染等多种生理病理的发生发展。CXCL11与受体CXCR3 结合后可促进下游一系列促炎因子IL-1β、IL-6及TNF-α等的表达，从而抑制病原体的增殖和扩散。但目前关于CXCL11相关信号通路的研究则较少[20]。PIAS1是STAT1特异性抑制子，具有小泛素样修饰剂SUMO-E3链激酶的活性，参与多种转录因子的活性调控和多种蛋白质的SUMO化修饰。另外，PIAS1 也可以作为抗炎因子负调控NF-κB/STAT1信号通路，从而抑制炎症损伤[21]。EIF2AK2 是由IFN诱导产生，在IFN抗病毒通路中扮演着重要的角色。EIF2AK2 又被称为蛋白激酶R(PKR)。PKR 最初是以非活性状态存在的，但在ds RNA或5′端磷酸化的 ss RNA等病原体RNA出现后，作为一种应答反应，PKR通过自磷的方式被激活，接下来PKR磷酸化 EIF2AK2 的51位丝氨酸，抑制其翻译的起始，进而阻止病原体复制[22]。这些基因可能是Ι型弓形虫感染过程中的重要靶点，可作为今后研究的关键分子。但由于本研究采用的是肿瘤细胞而非原代细胞，有些通路不一定反应体内真实的细胞事件，而且仅限于对弓形虫肺炎的机制解释。

综上，转染ROP16Ⅰ基因后影响宿主细胞基因表达谱，CCL8、CXCL11、PIAS1及EIF2AK2 等候选基因的发现为进一步研究Ι型弓形虫疾病的致病机理提供了新的思路，但目前其具体的作用机制仍不清楚，文献也未见报道。关于这些候选基因在弓形虫疾病发生发展中的作用仍需进一步证实。

利益冲突：无