利用加权基因共表达网络分析牛口蹄疫病毒感染通路变化

杨生海，刘西兰，张 勇

(1.甘肃农业大学 动物医学院，甘肃 兰州 730070； 2.中农威特生物科技股份有限公司，甘肃 兰州 730046)

口蹄疫(foot-and-mouth disease，FMD)是一种在偶蹄动物间常见的具有急性、热性、高度接触性传染的病毒性疾病，易感动物达70多种，包括牛、羊、猪、鹿、骆驼等其他偶蹄动物，一旦暴发，就会迅速在易感动物之间传播，引发大范围疫情，造成暴发疫情地区的重大经济损失，被世界动物卫生组织(OIE)列为A类动物疫病之首[1]。口蹄疫病毒(FMDV)有 O、A、C、SAT1、SAT2、SAT3、AsiaⅠ型7个血清型，每个血清型又包含若干血清亚型。该病毒不仅在各血清型之间无交叉免疫性，相同血清型的各亚型之间也仅有部分交叉免疫性[2]。

FMDV由衣壳蛋白包裹着包含约8 500个核苷酸的单链RNA组成[3]。FMDV单股正链RNA包括5′非编码区(5′-UTR)，一个开放阅读框(ORF)以及 3′-UTR。ORF编码一种多聚蛋白，被病毒和细胞的蛋白酶切割，产生构成病毒衣壳的成熟结构蛋白及复制依赖的非结构蛋白[4]。ORF根据其不同的功能被区分为4个功能区域：L、P1、P2 和P3，其中L区域编码Lpro蛋白，该蛋白可水解自身多聚蛋白[5]；P1区域编码衣壳蛋白，经蛋白酶分解后得到组成衣壳的成熟蛋白，包含VP1、VP2、VP3和VP4，VP2与VP4是VP0水解得到的；P2和P3区域分别编码2A、2B、2C蛋白以及3A、3B、3Cpro、3Dpol蛋白，整个病毒衣壳由60个原聚体构成，这60个原聚体在细胞内组装成一个二十面体的衣壳，衣壳直径大约30 nm[6-7]。但目前涉及FMDV持续感染的机制仍不清楚[8-9]，因此，了解FMDV感染的分子机制，对寻找分子靶标和研发新型高效疫苗具有重要意义。

本研究采用生物信息学方法，通过加权基因共表达网络分析(WGCNA)，筛选可能与FMDV 感染相关的通路及生物学功能。通过以上生物信息学方法从整体基因角度，了解FMDV 感染发生进展过程中的关键基因、信号通路，为治疗FMDV 寻找新的理论基础。

1 材料与方法

1.1 数据集的标准化

从NCBI的GEO数据库(https://www.ncbinlm.nih.gov/geo/)下载编号GSE83514的芯片数据集。GSE83514芯片数据集是基于GLP22061平台，包含20头荷斯坦奶牛样本，其中FMDV 阳性(FMDV+)成年病牛11头，FMDV 阴性(FMDV-)成年牛7头， FMDV 阴性(FMDV-)幼崽2头。样品后续均进行了cy3和cy5的荧光标记。从GEO数据库下载原始数据，利用R软件中的Affy包和limma包读取原始文件，使用RMA算法进行标准化。

1.2 WGCNA分析

利用变异系数为筛选条件，挑选最高的5 000个基因进行后续分析。使用R软件的WGCNA包构建基因共表达网络分析，为了确保符合无尺度网络分布，WGCNA需要选择合适的加权系数β。利用WGCNA包中的“pickSoftThreshold”函数计算β值的相关系数和基因连接度均值，选择适当的软阈值β使得构建的网络更符合无标度网络的标准。采用一步法构建相关的基因网络，将邻接矩阵转化为拓扑重叠矩阵TOM，利用层次聚类产生一个基因的层次聚类树。计算基因显著性(gene significance， GS)和模块显著性(module membership， MS)，用以衡量基因与临床信息的显著性，并分析模块及模型的显著关联。

1.3 模块内基因网状图构建及枢纽基因筛选

利用WGCNA包中的“exportNetworkToCytoscape”函数导出基因间的网络关系，并使用Cytoscape软件作图，利用WGCNA包中的“chooseTopHubInEachModule”函数挑选每个模块唯一的枢纽基因。

1.4 GO和KEGG富集分析

使用R软件中的“Clusterprofiler”包对特定的模块基因进行GO注释和KEGG通路富集分析，以P值≤0.05为差异具有统计学意义。

2 结果与分析

2.1 软阈值选择

以变异系数最高5 000个基因，为构建符合无尺度的网络模块，同时降低计算的复杂性和时间，使用pick Soft Threshold函数来权重参数β值，在正常样本和感染FMDV 的样本中，选取了1～20个软阈值，对每个软阈值计算相关系数和平均连通性。结果显示，当β=8时，相关系数r>0.8，此时，该软阈值对应的网络平均连接度接近零，说明此网络近似于无尺度网络。因此，本文选取软阈值β=8来构建基因共表达网络(图1)。

图1 软阈值的筛选

2.2 基因模块识别

按照在选择加权系数β之后，获得基因之间的相异度矩阵dissTOM，对该相异度矩阵dissTOM进行层次聚类得到基因聚类树(图2-A)，使用混合动态树剪切法从该基因聚类树中识别出30个模块(图2-B)，对这30个模块以0.80的相关性合并相似度高的模块后得到20个模块(图2-C)。即最后将5 000个基因分成了20个共表达网络模块。结果表明，成功构建了无尺度网络并完成了基因的模块划分。

A， 基因聚类树；B，动态树；C，动态树合并。

2.3 FMDV感染相关的模块及枢纽基因

WGCNA中模块与性状相关性热图的提示，有8个基因模块与FMDV 显著相关，我们选择正相关系数最高(Green模块，相关系数r=0.57，P<0.05)和负相关系数最高(Darkred模块相关系数r=-0.48，P<0.05)的两个模块进行后续研究(图3-A)。图3-B表示Darkred模块内连通性和基因重要性的相关性(P<0.05)，表明这个模块内每个基因与整个模块呈显著负相关(P<0.05)。图3-C提示Green模块内连通性和基因重要性的相关性，也同样表明该模块内每个基因与整个模块呈显著(P<0.05)正相关。接着，挑选WGCNA中导出的weight值最大的150条edge构建模块内基因的网络图，并将枢纽基因用黄色标注。其中，Darkred模块内的枢纽基因为MFSD4(图3-D)，Green模块内的枢纽基因为RHOH(图3-E)。上述结果表明，Darkred模块和Green模块内的基因与FMDFV感染显著相关，且MFSD4和RHOH可能在其中发挥着重要作用。

2.4 模块的富集分析

对上述Darkred模块和Green模块中的基因进行富集分析，寻找与FMDV感染相关的生物学过程和KEGG通路。Darkred模块GO富集结果提示，FMDV与免疫系统过程、细胞周期、免疫反应、DNA代谢过程、细胞分裂、有丝分裂细胞周期、DNA复制、B细胞活化和增殖等生物学过程相关(图4-A)。Green模块GO富集结果提示，FMDV与中枢神经系统发育、中枢神经系统神经元分化、磷脂酰肌醇3激酶信号转导、中枢神经系统神经元发育、中枢神经系统投射神经元轴突形成、中枢神经系统神经元轴突形成、神经系统过程的负性调节、次生代谢过程的调控、白脂肪细胞分化和多巴胺能神经元分化等GO注释相关(图4-B)。结果表明，以上生物学过程与FMDV相关。

Darkred模块KEGG通路富集结果提示，FMDV感染与Rap1信号通路、黏附连接、DNA复制、核苷酸切除修复等信号通路相关(图5-A)。Green模块KEGG通路富集结果提示，FMDV 感染细胞因子-细胞因子-受体相互作用、细胞周期、B细胞受体信号通路、肠道免疫网络IgA、DNA复制和原发性免疫抑制等信号通路相关(图5-B)，表明以上KEGG通路与FMDV相关。

3 讨论

本研究利用WGCNA对GSE83514芯片进行生物信息学分析，将基因划分成了20个模块并分析了各个模块与FMDV 感染的相关性。结果显示，Darkred模块与FMDV 感染呈显著负相关，Green模块与FMDV 感染呈显著正相关；且Darkred模块的枢纽基因为MFSD4，Green模块的枢纽基因为RHOH。对以上两个模块内的基因进行GO功能注释和KEGG通路富集，共得到20个GO注释及14条KEGG通路，提示以上生物学过程和通路与FMDV 感染密切相关。

上述的多条GO注释、KEGG通路已被报道参与FMDV 感染过程中。在免疫过程方面，FMDV感染早期，原发性免疫会被抑制，使得病毒能在动物系统中快速增殖并传播至多个感染部位[10-11]。FMDV感染可以导致宿主内的淋巴细胞减少[12]、抑制树突状细胞的功能[13]、NK细胞的功能紊乱[14-15]。与原发免疫系统密切相关的B细胞受体通路在FMDV感染或FMDV疫苗研发中也发挥着重要作用。血浆内的浆细胞及记忆B细胞被激活后可以产生抗体，后者可抵抗病毒感染，产生长期的血清学保护。FMDV的衣壳蛋白以TI-2抗原形式直接与B细胞受体结合，促使细胞产生抗体[16-17]。

在分子机制方面，最为经典的FMDV 感染机制是：FMDV感染细胞时，TLR受体被激活，通过TANK/TRAF3-NEMOTBK1IKKε-IRF3/7通路，最终促进IFNα、IFNβ、IFNγ及IL21等细胞因子的表达。产生的Ⅰ型、Ⅱ型IFN还可以与细胞膜上的IFNR1、IFNR2相互作用，激活下游JAK/STAT通路，由此证实细胞因子-细胞因子受体作用通路在其中发挥着重要作用[18-19]。

当然也有学者进行了FMDV的体外感染实验，寄希望于寻找新的感染机制：将感染了FMDV的BHK-21细胞与未感染细胞进行了转录组测序，通过对4 686个差异基因富集分析发现，细胞周期通路可以被显著富集，其中细胞周期相关蛋白CCND1存在明显差异[20]

然而目前，模块枢纽基因MFSD4和RHOH基因在FMDV 感染中的直接作用尚未报道。RHOH基因所编码的蛋白质是鸟苷三磷酸(GTP)代谢酶的Ras超家族的成员，其GO注释包括GTPase活性、蛋白质结合相关，其主要参与白细胞跨内皮迁移(llleukocyte transendothelial migration)信号通路和沙门氏菌感染(Salmonellainfection)通路[21-22]。既往有文献报道，白细胞中的外周血单核细胞(PBMC)跨内皮迁移参与口蹄疫疫苗的免疫应答，并在其中起着重要的免疫调节作用[23]。MFSD4是主要易化子超家族成员之一，目前已知该基因编码的蛋白主要定位于质膜上，其GO注释主要包括参与膜成分的构成，如质膜的组成部分(integral component of plasma membrane)、膜(membrane)和膜的组成部分(integral component of membrane)[24]，但具体功能尚无完善细致的研究。这两个基因仍有待学者进一步研究。

虽然以上两个枢纽基因并没有被报道发挥重要功能，但是Green模块和Darkred模块内富集到多个与FMDV相关的经典基因，如Green模块内的TLR9、TRAF3IP3、IRAK1BP1、MAP3K6、OAS1和IL1A，Darkred模块内的MAP3K13和NFKB1等。以上基因早已被证实在FMDV感染过程中发挥着重要作用[14，25-26]。当细胞感染FMDV后，病毒RNA由胞浆中的RLRs和/或膜结合TLRs(主要是TLR3、7、8、9)检测，触发信号转导级联，最终激活转录因子，促使下游如IFN、促炎细胞因子(如IL1等)和ISGs(如OAS1、PKR和Mx1等)的分泌和产生。首先，TLR9可以与TRIF和MyD88相互作用，组成复合物检测胞内游离的RNA。MyD88继续与下游TRAF3/6与 IRAK1/4蛋白相互作用，促使NEMO感知并招募TBK1或IKK家族蛋白，或者直接激活MAPK3。其中NEMO/TBK1复合物主要导致IRF3/7的磷酸化并入核发挥转录因子的作用。IKK复合物可磷酸化IκB，导致其降解和NF-κB脱离，然后转移到胞核内。核NF-κB可与IFNβ、TNFα、IL-1、OAS1等的启动子结合促进转录。同时，TRAF3/6可激活MAPK3，ERK1/2和JNK，这些激酶可促进IFNβ的转录。总的来说，以上经典基因的出现证实了本研究方法的可靠性。

综上所述，本研究通过构建基因共表达网络，筛选出与FMDV 感染相关的2个枢纽基因(MFSD4和RHOH)，20个GO注释及14条KEGG通路。以上结果均有助于揭示FMDV 感染的分子机制，为FMDV感染的诊断及治疗提供新的靶点。