基于S基因序列分析新型猪流行性腹泻病毒的遗传演化

猪流行性腹泻病毒(porcine epidemic diarrhea virus, PEDV)是导致猪群腹泻、仔猪死亡的主要病原之一，其引发疾病的临床特征以水样腹泻、呕吐、消瘦为主，严重降低新生仔猪存活率[1]。猪流行性腹泻病(porcine epidemic diarrhea, PED)具有接触传播的特性，可通过污染的饲料、动物制品传播，极大地限制了猪及其副产品的进出口贸易，带来巨大的经济损失[2-6]。因此对该病的防控和监测显得极为重要。

PEDV是甲型冠状病毒属成员[7]，其全基因约为28 kb。基因组3′端约三分之一的序列负责编码毒株的4个结构蛋白，分别为纤突(S)蛋白、包膜(E)蛋白、膜(M)蛋白和核衣壳(N)蛋白。其中，S蛋白负责结合病毒受体、介导宿主细胞的融合以及诱导机体产生中和抗体，是重要的毒力蛋白[8]。S基因负责编码S蛋白，具有较高的突变性，是毒株毒力的指征性基因，其变异和重组可能造成PEDV毒力的改变，产生新型毒株。因此监测S基因的遗传变异成为PEDV流行病学研究中的重要任务。目前，根据对S基因的遗传进化研究，可将PEDV分为G1和G2群，在此基础上又可进一步分为G1a、G1b、G2a和G2b亚群[9]，当下在中国流行的毒株以G2a和G2b亚群为主[10-11]。最近，本实验室在青藏高原藏猪腹泻粪便中发现了两种新G1和G2亚群的PEDV毒株，且已在高原地区流行[12]。为了进一步调查新型藏猪源PEDV是否在四川仔猪中存在或流行，本研究对本实验室2018—2019年保存的28个猪场的116份腹泻样本进行了PEDV的检测，并对其S基因进行遗传进化研究。

2018年2月—2019年6月期间西南民族大学预防兽医实验室保存的来自四川省西昌、内江、资阳、乐山、雅安、绵阳、崇州和自贡的28个规模化猪场116份仔猪腹泻的粪便或肠道组织。样本于-80 ℃冻存。

Trizol总RNA提取试剂、反转录试剂盒、DL2000TM DNA Marker购自TaKaRa生物技术有限公司；Quick Taq HS DyeMix购自TOYOBO生物科技有限公司；凝胶成像仪购自上海天能科技有限公司；TGL-16冷冻离心机购自四川蜀科仪器有限公司；PCR仪购自杭州博日科技有限公司。

剪碎腹泻病猪的部分小肠或粪便样品，加入去离子水制成匀浆，反复冻融。4 000 r·min-1离心10 min。取300 μL上清液，加入Trizol裂解液静置数分钟。再以氯仿分相，取上层水相加入异丙醇，低温高速离心沉淀核酸。以75% DEPC乙醇洗去杂蛋白，得到较纯的RNA。利用TaKaRa反转录试剂盒进行反转录。

利用PEDV检测引物[12]，以RT-PCR方法对116个样本进行检测，PCR反应体系：Quick Taq HS DyeMix 5 μL, 上、下游引物各0.5 μL, cDNA模板1 μL, ddH2O 3 μL。PCR反应条件：94 ℃ 2 min；94 ℃ 30 s, 53 ℃ 30 s, 68 ℃ 50 s，循环36次；72 ℃ 8 min。以琼脂糖凝胶电泳法对PCR产物进行检测，凝胶浓度为1%。

在每个PEDV阳性猪场中选取1~2个阳性样本的cDNA用于S基因的扩增。利用Qin等[12]设计的5对S基因扩增引物，对其中17份阳性病料的cDNA模板进行RT-PCR扩增。PCR反应体系：Quick Taq HS DyeMix 12.5 μL, 上、下游引物各1 μL, cDNA模板1.5 μL, ddH2O 9 μL。PCR反应条件：94 ℃ 2 min；94 ℃ 30 s，54 ℃/47 ℃/51 ℃ /57 ℃/56 ℃ 30 s，68 ℃ 1.5 min，循环36次；72 ℃ 8 min。将扩增产物送上海生工生物工程有限公司进行测序。

运用Seqman软件拼接序列，得到13条完整的S基因序列。利用MegAlign软件将获得的序列与50条S基因参考序列进行比对分析。通过网站http://www.cbs.dtu.dk/services/NetNGlyc/预测S蛋白N-糖基化位点，通过网站http://crdd.osdd.net/raghava/abcpred/预测COE区域的B细胞抗原表位。利用MEGA 7.0软件的邻近法(Neighbor-Joining)构建系统进化树。以RDP 4.39软件的Chimaera、MaxChi、BootScan、RDP、GENECONV、SiScan和3Seq方法进行重组分析，利用SimPlot 3.5.1软件进行相似性分析和BootScan分析，对重组区和非重组区构建独立的进化树以验证重组。

以贝叶斯进化分析软件包(BEAST v1.8.0)的BEAUti程序设置序列的采样时间、核苷酸置换模型、分子钟模型、马尔科夫链长和先验参数，输出XML文件。以BEAST程序运算生成的XML文件，通过MCMC方法进行树的重复抽样，获得具有最大后验概率的进化树。通过TreeAnnotator程序设置Burnin值，舍弃老化样本，输出最终的进化树文件。以FigTree v1.4.4软件对进化树进行修饰和分歧时间的标定，分析毒株的演化过程。

利用PEDV检测引物[12]，以RT-PCR方法对28个规模化猪场的116份腹泻样本进行检测，结果表明：在粪便或肠组织样本中，PEDV阳性率为42.2% (49/116，95%CI=33.1%~51.8%)；28个规模化猪场中，PEDV检出率为57.1% (16/28，95%CI=37.2%~75.5%)。结果说明PEDV依然是导致猪群腹泻的主要原因。

利用5对S基因扩增引物[12]，对17份PEDV阳性的cDNA进行S基因的扩增，利用SeqMan软件拼接得到13条完整的S基因核苷酸序列，长度为4 149~4 170 bp。将序列命名并上传至GenBank，分别为SWUN-2-CH-SCXC-2018(GenBank accession No. MK 592416)、SWUN-1-CH-SCNJ-2019(No.MK685665)、SWUN-4-CH-SCXC-2018(No.MK598821)、SWUN-5-CH-SCXC-2018(No.MK598820)、SWUN-17-CH-SCZY-2018(No.MK598819)、SWUN-18-CH-SCLS-2019(No.MK820037)、SWUN-19-CH-SCZY-2018(No.MK592415)、SWUN-C6-CH-SCYA-2019(No.MN161584)、SWUN-MY-CH-SCMY-2019(No.MK820039)、SWUN-H1-CH-SCYA-2019(No.MK820040)、SWUN-Y1-CH-SCCQ-2019(No.MK820041)、SWUN-3CH-CH-SCZG-2019(No.MK820038)和SWUN-H3-CH-SCYA-2019 (No.MK820042)。

运用MegAlign进行序列同源性分析，结果表明：扩增得到的13条PEDV S基因核苷酸序列相似性为94.2%~99.9%，氨基酸序列相似性为94.0%~99.8%。13条PEDV S基因核苷酸序列与2014年的强致病株YC2014相似性最高，达到95.5%~99.1%，与经典毒株CV777相比相似性较低，为92.4%~94.9%；与本实验室的7条藏猪源PEDV S基因核苷酸序列相比，相似性为93.5%~98.6%。值得注意的是，SWUN-H3-CH-SCYA-2019与5株藏猪源新G1亚群PEDV的S基因序列相似性高达97.0%~98.6%，氨基酸序列相似性为95.3%~98.2%；SWUN-1-CH-SCNJ-2019、SWUN-3CH-CH-SCZG-2019、SWUN-4-CH-SCXC-2018和SWUN-19-CH-SCZY-2018与2株藏猪源新G2亚群PEDV的S基因序列相似性达到96.5%~98.5%，氨基酸序列相似性为95.2%~98.5%。

PEDV S基因的S1区负责编码S蛋白的1—789 aa，包含了多个重要的抗原表位，是PEDV S基因中变异程度最高，遗传多样性最大的区域[13-14]。为进一步观察序列的变异情况，本研究利用MEGA 7.0软件的ClustalW方法进行S蛋白序列的比对分析。结果发现推导出的13条氨基酸序列均在CV777和YC2014的基础上出现了氨基酸的替换，在S1区的氨基酸变异情况最为明显(图 1)。其中SWUN-1-CH-SCNJ-2019出现了1处氨基酸的插入：376RQSE。值得注意的是，SWUN-H3-CH-SCYA-2019在S1区的变异特征与藏猪源新G1亚群的PEDV非常相似，存在39处共同的氨基酸替换，推断SWUN-H3-CH-SCYA-2019的S基因序列可能在S1区与藏猪源新G1亚群PEDV重组。

PEDV的S蛋白包含4个主要的中和表位：S1区的COE(499—638 aa)、SS2(748—755 aa)、SS6(764—771 aa)和S2区的2C10(1 368—1 374 aa)[15]。本研究的毒株在SS2和2C10区相对保守，主要在COE和SS6区出现氨基酸的替换：在COE区，新G2亚群的SWUN-19-CH-SCZY-2018、SWUN-4-CH-SCXC-2018、SWUN-1-CH-SCNJ-2019、SWUN-3CH-CH-SCZG-2019和藏猪源新G1亚群的SWUN-H3-CH-SCYA-2019出现一致突变524P/H和529I/L；SWUN-19-CH-SCZY-2018和SWUN-4-CH-SCXC-2018出现一致突变611S/G；SWUN-1-CH-SCNJ-2019出现独特的突变570K/T、612E/D和637K/T；SWUN-19-CH-SCZY-2018出现独特的突变636E/V；SWUN-3CH-CH-SCZG-2019出现独特的突变520S/A。在SS6区，SWUN-19-CH-SCZY-2018、SWUN-4-CH-SCXC-2018和SWUN-H3-CH-SCYA-2019出现一致突变766L/P；SWUN-1-CH-SCNJ-2019和SWUN-3CH-CH-SCZG-2019出现一致突变766L/S。

糖基化位点预测结果表明：与藏猪源PEDV相比，本研究的新亚群PEDV毒株SWUN-4-CH-SCXC-2018、SWUN-19-CH-SCZY-2018、SWUN-1-CH-SCNJ-2019、SWUN-3CH-CH-SCZG-2019和SWUN-H3-CH-SCYA-2019均出现了糖基化位点的引入或缺失。其中新G2亚群的SWUN-4-CH-SCXC-2018、SWUN-19-CH-SCZY-2018、SWUN-1-CH-SCNJ-2019和SWUN-3CH-CH-SCZG-2019引入了预测的糖基化位点118NATA121，SWUN-4-CH-SCXC-2018、SWUN-19-CH-SCZY-2018、SWUN-1-CH-SCNJ-2019和藏猪源新G1亚群的SWUN-H3-CH-SCYA-2019缺失了预测的糖基化位点131NKTL134，SWUN-1-CH-SCNJ-2019引入了预测的糖基化位点867NISS870，SWUN-3CH-CH-SCZG-2019引入了预测的糖基化位点1 196NHTA1 199，SWUN-H3-CH-SCYA-2019引入了预测的糖基化位点1 192NHTA1 195。本研究G2b亚群PEDV毒株与YC2014相比，仅SWUN-5-CH-SCXC-2018缺失了预测的糖基化位点743NCTE746，其余毒株均未出现糖基化位点的引入或缺失。B细胞抗原表位预测结果表明：13株PEDV均在COE区域的23—28、36—39、63—74和110—114 aa(S蛋白的521—526、534—537、561—572、608—612 aa)处预测到存在较大的成为B细胞抗原表位的可能。有所不同的是，G2b亚群的SWUN-2-CH-SCXC-2018、SWUN-5-CH-SCXC-2018和新G2亚群的SWUN-19-CH-SCZY-2018因一处突变(636E/V)，导致COE区域的129—138 aa(S蛋白的627—636 aa)增加一处预测的抗原表位：LITGTPKPLV。这些结果提示，毒株的S蛋白已发生了不同程度的突变，可能会导致毒株抗原性的改变。

为进一步探究毒株的亲缘关系，运用MEGA 7.0软件的Neighbor-Joining法，设置自展值为1 000，采用p-distance模型，对63条国内外PEDV S基因核苷酸序列构建遗传进化树(图 2)。将建好的进化树分为G1和G2大群，再将G1大群分为1a和1b亚群，G2大群分为2a和2b亚群。由图 2可见，在本研究的13株PEDV中，有12株位于G2大群，唯独SWUN-H3-CH-SCYA-2019位于G1大群，与藏猪源新G1亚群PEDV亲缘关系较近。12株PEDV中，有8株分布在G2b亚群，与中国其他省区毒株亲缘关系较近。另外4株PEDV(SWUN-4-CH-SCXC-2018、SWUN-19-CH-SCZY-2018、SWUN-1-CH-SCNJ-2019和SWUN-3CH-CH-SCZG-2019)与CH/SCBC/2018处在一个新分支中，且与藏猪源PEDV亲缘关系较近，可能是一个新的G2亚群。以上结果表明四川地区已经出现了与藏猪源PEDV S基因遗传关系较近的PEDV毒株。

运用RDP 4.39软件的Chimaera、MaxChi、BootScan、RDP、GENECONV、SiScan和3Seq方法进行重组分析，发现仅SWUN-H3-CH-SCYA-2019具有明显的重组现象，其余12株PEDV序列由于无法在SimPlot软件中断点，重组现象不明显。有趣的是，所有方法均支持SWUN-H3-CH-SCYA-2019与主要亲本株CH/TP-D15/2018和次要亲本株KNU-1707发生重组，断点范围为679—869 bp，同藏猪源新G1亚群的CH/TP-A9/2018、CH/TP-4-3/2018、CH/TP-3-1/2018、CH/TP-2-2/2018和CH/TP-E5/2018有着一致的重组模式[12]。为进一步验证重组，选取重组序列SWUN-H3-CH-SCYA-2019、主要亲本序列CH/TP-D15/2018、次要亲本序列KNU-1707和一条外群序列P5-V，利用SimPlot 3.5.1软件进行相似性分析(图 3A)和BootScan分析(图 3B)，对重组片段进行断点。图中竖线(breakpoint)为SWUN-H3-CH-SCYA-2019的断点位置，结果为742 bp，处在RDP 4.39软件计算的679—869 bp内，验证了重组检测的结果。以断点位置为界，将S基因划分成重组(1—741 bp)和非重组(742 bp—3′end)两个区域，通过MEGA 7.0软件的Neighbour-Join法对两个区域的核苷酸序列分别构建进化树(图 3C)。在重组区进化树中，SWUN-H3-CH-SCYA-2019归属于藏猪源新G1亚群；而在非重组区进化树中，SWUN-H3-CH-SCYA-2019远离了藏猪源亚群，归属于G2b亚群。这些结果表明SWUN-H3-CH-SCYA-2019与藏猪源新G1亚群的PEDV重组模式一致。

为深入了解PEDV的遗传演化过程，对63条确定采样时间的S基因序列以BEAST v1.8.0软件包进行分歧时间的估算。以BEAUti程序设置各序列的采样时间、核苷酸置换模型、马尔科夫链长和先验参数，选取松散的对数正态分布分子钟模型，输出XML文件。以BEAST程序的MCMC法对BEAUti输出模型的拓扑结构、分支长度以及核苷酸置换模型等参数进行重复抽样，获得稳定的参数分布，输出具有最大后验概率值的贝叶斯进化树。通过TreeAnnotator程序舍弃老化样本，输出最终的进化树文件。以FigTree v1.4.4进行进化树的修饰和分歧时间的标定后(图 4)，估算出藏猪源新G1亚群毒株SWUN-H3-CH-SCYA-2019的分歧时间约为2012.3年，早于其余5株藏猪源新G1亚群PEDV CH/TP-2-2/2018、CH/TP-3-1/2018、CH/TP-4-3/2018、CH/TP-E5/2018和CH/TP-A9/2018的最早分歧时间(2015.7年)；新G2亚群毒株SWUN-4-CH-SCXC-2018、SWUN-19-CH-SCZY-2018和SWUN-3CH-CH-SCZG-2019的分歧时间约为2014.2年，早于藏猪源新G2亚群PEDV CH/TP-4-4/2018和CH/TP-D15/2018的分歧时间(2014.7年)。结果表明，两个亚群的藏猪源PEDV的分歧时间都晚于亲缘关系较近的四川毒株。

PEDV是引起猪腹泻病的主要病原之一，严重影响猪的生产性能，且在新生仔猪上表现出高致死性，给我国规模化养猪业带来巨大的经济损失。为了防控PED，人们通过预防接种的方法，取得过较好的免疫效果。然而，当下PEDV正发生着变异和重组[12, 16-17]，导致新型毒株不断出现，抗原性不断改变，从而使该病的防控变得越来越艰难[18]。随着对PEDV毒株毒力分子机制研究的深入[8, 19]，监测PEDV毒力基因的变化成为其流行病学研究中的重要任务。

最近，本实验室在青藏高原发现两种新亚型PEDV[12]，为进一步调查新型PEDV是否在四川地区腹泻猪群中存在或流行，本研究对2018—2019年间保存的四川地区28个规模化猪场的116份腹泻或肠道样本进行PEDV的检测，选取部分检测为阳性的样本进行S基因的扩增、测序和遗传进化研究。结果发现四川地区的样本阳性率为42.2% (49/116，95%CI=33.1%~51.8%)，猪场检出率为57.1% (16/28，95%CI=37.2%~75.5%)，与全国范围的样本阳性率49.58%(686/1 383，95%CI=46.9%~52.3%)和猪场检出率57.83% (4 060/7 021，95%CI=56.7%~59.0%)接近[19]，说明PEDV依然是导致猪群腹泻的主要原因。将本研究13条推导的氨基酸序列与藏猪源序列进行同源性分析，结果显示SWUN-4-CH-SCXC-2018、SWUN-19-CH-SCZY-2018与藏猪源PEDV CH/TP-4-4/2018的氨基酸序列相似性高达98.3%~98.5%，SWUN-H3-CH-SCYA-2019与藏猪源PEDV CH/TP-3-1/2018、CH/TP-E5/2018和CH/TP-2-2/2018的氨基酸序列相似性高达98.1%~98.2%，由此猜测藏猪源PEDV已经在四川地区出现。对S基因推导的氨基酸序列进行对位排列，发现其S1区(1—789 aa)差异最大，并且可以看到SWUN-H3-CH-SCYA-2019的氨基酸序列与藏猪源毒株CH/TP-A9/2018、CH/TP-4-3/2018、CH/TP-3-1/2018、CH/TP-2-2/2018和CH/TP-E5/2018高度相似，明显有别于本研究其余12条序列。

S蛋白是PEDV重要的毒力蛋白，其氨基酸序列决定了蛋白的亲水性、可塑性等性质[20]，影响病毒与抗体和免疫细胞受体结合的过程。蛋白质糖基化是糖链与蛋白质连接的重要的修饰过程，很大程度上影响着机体的免疫反应[21-22]。糖基化位点是蛋白质上结合糖链的氨基酸序列，其突变可能引起糖基化作用的改变，导致毒株抗原性改变。本研究对PEDV S蛋白氨基酸全序列进行了N-糖基化位点预测，对S蛋白的核心抗原表位区(COE)进行了B细胞抗原表位的预测，结果表明：本研究5株新亚群的PEDV(其中4株新G2亚群PEDV SWUN-4-CH-SCXC-2018、SWUN-19-CH-SCZY-2018、SWUN-1-CH-SCNJ-2019、SWUN-3CH-CH-SCZG-2019和1株藏猪源新G1亚群PEDV SWUN-H3-CH-SCYA-2019)与其余藏猪源毒株相比，均出现了糖基化位点的引入或缺失；其余8株G2b亚群的PEDV与强致病株YC2014相比，仅SWUN-5-CH-SCXC-2018出现了糖基化位点的缺失，部分毒株虽出现相关位点的氨基酸突变，但对糖基化位点的预测结果影响不大。B细胞抗原表位的预测结果表明：13株PEDV在COE区预测到了4处共同的B细胞的抗原表位。值得注意的是，SWUN-2-CH-SCXC-2018、SWUN-5-CH-SCXC-2018和SWUN-19-CH-SCZY-2018因一处氨基酸替换(636E/V)，导致COE区域的129—138 aa处增加一处预测的抗原表位——LITGTPKPLV。这些结果提示：四川地区的PEDV中，部分毒株的抗原性可能已经发生改变，其具体的生物学功能有待进一步研究。

对国内外63株PEDV的S基因进行遗传进化研究，结果表明本研究13株PEDV中有12株位于G2大群，其中又有8株分布在G2b亚群，与其他省区毒株亲缘关系较近，另外4株PEDV(SWUN-4-CH-SCXC-2018、SWUN-19-CH-SCZY-2018、SWUN-1-CH-SCNJ-2019和SWUN-3CH-CH-SCZG-2019)临近藏猪源新G2亚群，聚成一簇，可能形成了一个新的G2亚群。值得注意的是，SWUN-H3-CH-SCYA-2019远离其余12株PEDV，单独处在藏猪源新G1亚群中，与5株藏猪源PEDV表现出较近的亲缘关系。再对本研究13株PEDV的S基因进行重组分析，发现仅SWUN-H3-CH-SCYA-2019的重组检测结果可以通过SimPlot断点和BootScan的验证，并且与新G1亚群藏猪源PEDV有一致的重组模式[12]。这些结果表明，本研究在四川地区发现了与藏猪源PEDV有着相似S基因和相同重组模式的毒株。

“分子钟”是遗传学和分子生物学建立的一种模型，为流行病的溯源、生物界的考古等研究提供理论依据，且在病毒的演化研究中已经得到证实[23]，使得学者能够通过分析遗传数据解决大量流行病学问题[22, 24-25]。贝叶斯遗传进化分析软件包(BEAST)运用贝叶斯统计推断的MCMC方法(一种类似bootstrap但更为严格的抽样法)，从随机建立的初始树开始，按照一定规则，对不断产生的新树予以接受或否认，经过大量的循环抽样，以具有最大后验概率的进化树作为最终的结果[26]，具有较准确的计算分歧时间的能力。本研究利用BEAST v1.8.0对63株PEDV的S基因进行了分歧时间的估算，结果表明藏猪源新G1亚群毒株SWUN-H3-CH-SCYA-2019的分歧时间约为2012.3年，早于其余5株藏猪源新G1亚群PEDV(CH/TP-2-2/2018、CH/TP-3-1/2018、CH/TP-4-3/2018、CH/TP-E5/2018和CH/TP-A9/2018)的最早分歧时间(2015.7年)；新G2亚群毒株SWUN-4-CH-SCXC-2018、SWUN-19-CH-SCZY-2018和SWUN-3CH-CH-SCZG-2019的分歧时间约为2014.2年，早于藏猪源新G2亚群毒株CH/TP-4-4/2018和CH/TP-D15/2018的分歧时间(2014.7年)。这些结果说明两个亚群的藏猪源PEDV分歧时间都晚于亲缘关系较近的四川毒株，是PEDV由四川传入藏区的有力证据。

对2018—2019年保存的四川地区116份腹泻样本进行PEDV的分子学检测，对病毒的S基因进行扩增和测序。从序列的相似性分析、遗传进化研究和重组分析3个方面都发现存在与藏猪源PEDV相似的四川毒株。最后通过比较估算的分歧时间，证明了青藏高原的PEDV源自四川地区。