中国百日咳疫苗生产用菌株完成图基因组遗传分析

李喆，吴燕，王丽婵，卫辰，廖平安，朱德武，瞿明霞，马霄

1.中国食品药品检定研究院卫生部生物技术产品检定方法及其标准化重点实验室，北京 102629；2.武汉生物制品研究所有限责任公司，湖北 武汉 430027

《中国药典》三部（2020版）在生物制品生产检定用菌种管理及质量控制中规定，应建立生产用菌毒种种子批全基因序列的背景资料，生产用菌毒种种子批应进行全基因测序。应对生产用菌毒种已知的主要抗原表位的遗传稳定性进行检测，并证明在规定的使用代次内其遗传性状是稳定的［1］。因此对于菌毒种而言，获得其完整、准确的全基因组数据，并保证其安全性及遗传稳定性，对于生产过程中的源头质量控制具有重要意义。

中国百日咳疫苗生产菌株CMCC58003（CS 株）于1951年分离自北京 1名婴儿患者。2011年，实验室研究人员在继往研究中通过二代测序（454 测序平台）与一代测序（鸟枪法）相结合的方式获得了该菌株的基因组（GenBank：CP002695）［2］。目前细菌基因组的二代测序存在单个片段读长较短，片段末端测序准确度不高的问题，导致拼接过程易造成片段缺失，较难获得样本的完成图基因组，因此需要辅以一代测序进行查缺补漏。CP002695即是由287个重叠群（contigs）拼接形成，但即使如此多的重叠群仍未完全将该基因组完全测通，中间的空白区域需通过单独设计引物进行一代测序补充［2］。而当前兴起的三代测序技术一次读长可达几十Kb 级别，实现了基因组的原生测序，直接获得样品的完成图基因组，非常适用于细菌基因组测序［3-6］。

预防百日咳的最佳方法是接种百日咳疫苗［7-9］。尽管世界卫生组织的报告表明，百日咳疫苗在全球范围内具有较高的接种覆盖率（https：／／www.who.int／immunization／monitoring_surveillance／data／en／），但自1990 年以来，百日咳在全球范围内多个国家发病率明显上升，俗称百日咳重现［10-13］。中国疾病预防控中心历年发布的法定传染病疫情概况表明，百日咳在中国的发病数近年来呈显著上升趋势。引起百日咳重现这一现象的潜在原因较多，其中流行菌株的变异（基因突变、基因组结构变化等）造成的免疫逃逸是其中重要的原因［14-16］。流行株的变异研究一般需要疫苗株作为比对［17-18］，因此背景完整清晰、结构准确的疫苗菌株基因组尤为重要。

基于上述背景，本研究利用中国医学细菌菌种保藏管理中心（National Center for Medical Culture Collections，CMCC）保藏的百日咳疫苗生产用菌株CS株，通过PacBio 公司的三代测序技术辅以二代测序技术，获得该菌株的完整基因组序列和图谱。同时对基因组序列进行编码基因、碱基修饰、基因组结构等多维度的数据分析。随后通过与原二代测序拼接结果进行基因组水平上的比对。同时与武汉生物制品研究所有限责任公司合作，对百日咳疫苗生产发酵过程中的不同阶段进行菌体基因组的提取与三代测序，从抗原基因水平和基因组水平对百日咳生产菌株的遗传稳定性进行评价。

1 材料与方法

1.1 菌株 中国百日咳疫苗生产菌株CMCC58003（CS株）由中国食品药品检定研究院中国医学细菌保藏管理中心保存。

1.2 主要试剂无菌脱纤维羊血购自北京路桥技术股份有限公司；碳琼脂基础培养基购自英国Oxoid公司；Wizard Genomic DNA Purification Kit购自美国Promega公司；琼脂糖购自国药集团化学试剂有限公司。

1.3 基因组提取 百日咳疫苗原始菌株：在负压条件生物安全柜中开启冻存菌种，生理盐水复溶后，一次性接种环蘸取适量菌液划线接种于20%碳琼脂羊血培养基平板上，37 ℃培养72 h，作为第1代复苏菌种；培养结束后，取适量第1代复苏菌种接种于20%碳琼脂羊血培养基平板上，37 ℃培养72 h，作为第2 代复苏菌种；用一次性接种环刮取适量第2代复苏菌种，按照Wizard Genomic DNA Purification Kit 使用说明书提取菌种基因组。百日咳疫苗生产发酵过程菌体（工作种子、种子罐、扩增罐、发酵罐）：按照Wizard Genomic DNA Purification Kit使用说明书进行操作。通过1%琼脂糖凝胶电泳进行基因组DNA的纯度和完整性检测。

1.4 全基因组测序及分析 全基因组DNA电泳检测合格后，构建三代测序文库，利用PacBio Sequel 三代测序平台进行测序，测序服务委托北京诺禾致源公司完成。测序完成后，利用Prodigal和Augustus软件进行基因组编码基因检测［19-20］；SMRT Link 软件进行碱基修饰检测［21］；RepeatMasker 软件进行重复序列检测［22］；IslandPath-DIOMB 软件进行基因岛检测［23］；tRNAscan-SE 软件进行非编码RNA 检测［24］；phiSpy软件进行前噬菌体预测［25］；COG（Cluster of Orthologous Groups）和KEGG（Kyoto Encyclopedia of Genes and Genomes）数据库进行功能与信号通路注释；Circos软件绘制全基因组图谱［26］；Mauve软件（version 2015-0226）进行基因组比较分析［27］；以二代测序的CS 株序列（CP002695）与本研究完成三代测序的百日咳疫苗原始菌株CMCC58003（CS 株）、武汉生物制品研究所有限责公司（厂家A）和成都欧林生物技术股份有限公司（厂家B）生产用百日咳疫苗生产菌株进行全基因组水平的共线性比较；对厂家A 百日咳疫苗发酵过程不同代次的菌种全基因组进行共线性比较。

1.5 三代测序准确性验证 采用PCR 法。以提取的基因组为模板，进行插入序列IS481、百日咳毒素启动子（pt promoter）、百日咳毒素（pt）、百日咳菌毛蛋白2（fim2）和百日咳菌毛蛋白3（fim3）基因的扩增和Sanger 测序，扩增引物及扩增片段信息见表1，引物由苏州金唯智生物科技有限公司合成。将Sanger 测序与三代测序结果进行比对，比较其一致性。

表1 引物序列及扩增产物理论大小Tab.1 Theoretical size of primer sequence and amplified products

1.6 数据采集及分析 使用Pacbio软件采集数据，使用Hifiasm软件处理数据。

2 结果

2.1 基因组DNA 的提取 基因组提取物经1%琼脂糖凝胶电泳分析，可见清晰的电泳条带，见图1。为满足三代测序要求，3份平行样品送测基因组浓度均>100 ng／μL，2.0>A260／280>1.8，2.2>A260／230>2.0。

图1 CMCC58003（CS株）基因组DNA琼脂糖凝胶电泳图Fig.1 Agarose gel electrophoretic profile of genome DNA of CMCC58003（CS）strain

2.2 CMCC58003（CS株）全基因组测序

全基因组为环形，由1 个重叠簇拼接组成，长度为4 141 631 bp，GC 含量67.72%。对获得的基因组序列进行了编码基因检测、碱基甲基化修饰、重复序列检测、基因岛检测、非编码RNA 检测、前噬菌体检测，结果如下。将序列提交至GenBank，获得编号为CP086368.1。

2.2.1 编码基因 编码基因4 031个，平均长度921 bp，基因区总长度占全基因组长度的89.66%。

2.2.2 碱基修饰 在细菌基因组中，碱基甲基化修饰的主要类型为m4C、m5C 和m6A。CMCC58003 中m4C 修饰在基因组中总位点数为6 736 个，其中编码区6 233个（92.53%），非编码区503个（7.47%）。m5C修饰在基因组中总位点数为8 672 个，其中编码区8 335 个（96.11%），非编码区337 个（3.89%）。m6A修饰在基因组中总位点数为1 688 个，其中编码区1 543个（91.40%），非编码区145个（8.60%）。

2.2.3 重复序列 重复序列在细菌基因组上的分布分为散在和串联重复序列。散在重复序列是散在方式分布于基因组内的重复序列，串联重复序列即相邻的、重复2 次或多次特定核酸序列模式的重复序列。散在重复序列：长末端重复序列（long terminal repeats，LTR）数量为17，总长度1 912 bp，平均长度112 bp；短散在重复序列（short interspersed nuclear elements，SINE）数量为10，总长度639 bp，平均长度64 bp；长散在重复序列（long interspersed nuclear elements，LINE）数量为6，总长度562 bp，平均长度94 bp。串联重复序列（tandem repeats，TR）数量为269，总长度56 829 bp，平均长度211 bp。

2.2.4 基因岛 某些区段基因密度比全基因组平均密度高，形成基因岛。CMCC58003（CS 株）包含基因岛数量26，总长度298 323 bp，平均长度11 474 bp。

2.2.5 非编码RNA 对于细菌，非编码RNA 的类型主要指tRNA、rRNA及sRNA 3种。CMCC58003（CS株）tRNA 数量为51，总长度3 918 bp，平均长度77 bp；5SrRNA 数量为3，总长度336 bp，平均长度112 bp；16SrRNA数量为3，总长度4 557 bp，平均长度1 519 bp；23SrRNA数量为3，总长度8 640 bp，平均长度2 880 bp；未检测到sRNA。

2.2.6 前噬菌体 整合在宿主基因组上温和噬菌体的核酸称为前噬菌体。检测发现，CMCC58003（CS株）基因组包含前噬菌体数量为7，总长度244 270 bp，平均长度34 895.7 bp。

2.3 三代测序准确性验证 成功对IS481、pt promoter、pt、fim2、fim3基因完成扩增，见图2。利用Bioedit 软件对不同基因片段的Sanger 测序结果与三代测序相应片段的序列进行比对，未发现差异。

图2 IS481、pt promoter、pt、fim2、fim3基因扩增产物电泳图Fig.2 Electrophoretic profile of PCR products of IS481，pt promoter，pt，fim2 and fim3

2.4 功能及信号通路注释 COG 将4 031 个蛋白编码基因按照蛋白功能注释分为24组（A~X），其中占比前5 位的依次为氨基酸转运和代谢、能量生产和转换、通用功能预测、前噬菌体与转座子、转录，见图3。KEGG 将4 031 个蛋白编码基因按照蛋白信号通路注释分为6个大组和35个小组，其中代谢、环境信息处理和细胞进程相关信号通路占比前3，见图4。

图3 CMCC58003（CS）株基因组蛋白编码序列COG注释Fig.3 COG functional annotation of protein coding sequences of CMCC58003（CS）strain

图4 CMCC58003（CS）株基因组蛋白编码序列KEGG注释Fig.4 KEGG signaling pathway annotation of protein coding sequences of CMCC58003（CS）strain

2.5 全基因组图谱分析 基因组环状图谱见图5，由外到里分别为编码基因（外侧代表正链编码基因，内侧代表负链编码基因）、COG基因功能注释结果、KEGG基因功能注释结果、GO 基因功能注释结果，ncRNA、基因组GC 含量［以窗口（染色体长度／1 000）bp，步长（染色体长度／1 000）bp来统计GC含量，向内的红色部分表示该区域GC 含量低于全基因组平均GC 含量，向外的绿色部分与之相反，且峰值越高表示与平均GC 含量差值越大］、基因组GC skew 值分布［窗口（染色体长度／1 000）bp，步长（染色体长度／1 000）bp，具体算法为（G-C）／（G+C），向内紫色部分表示该区域G 的含量低于C 的含量，向外的浅绿色部分与之相反］。

2.6 比较基因组分析二代测序的CS株序列（CP002695）与三代测序的CMCC58003（CS 株）（CP086368.1）在宏观上存在较多部位的倒置（浅绿色）、易位／倒置（蓝色，红褐色）区域，而CP086368.1与疫苗生产厂家A、B 在基因组水平上保持了较好的一致性，见图6。CP002695 和CP086368.1 的GC 百分比分别为67.3%和67.72%，预测编码基因数量分别为3 456和4 031，平均基因长度分别为981 和921 bp，tRNA 基因数量均为51，rRNA 数量均为3。单核苷酸多态性分析（single nucleotide polymorphism，SNP）表明，二者之间存在131处位点差异，但均在非编码区域。

图6 CMCC58003（CS 株）二代测序与三代测序共线性分析Fig.6 Whole genome collinearity analysis of second-generation sequencing and third-generation sequencing for CMCC58-003（CS）strain

2.7 菌种基因组遗传稳定性 武汉生物制品研究所有限责任公司百日咳疫苗生产菌株发酵过程不同代次菌种全基因组比较分析发现，不同代次菌种间未发现易位、倒置、插入和缺失等基因组水平的差异，见图7。同时，从基因组中单独提取出pt、pha、prn、fim2、fim3、act、tct、ompA、BrKA、rplD、BP1569等具有免疫原性或潜在免疫原性的基因进行不同代次序列的比对，未发现变异。

图7 CMCC58003（CS 株）发酵过程不同代次基因组共线性分析Fig.7 Whole genome collinearity analysis of different generations during fermentation of CMCC58003（CS）strain

3 讨论

疫苗生产菌株遗传背景清晰是疫苗生产的前提和保障。对于百日咳疫苗这种多组分疫苗，单个抗原基因的检测不能满足新型百日咳疫苗的发展趋势，新型百日咳疫苗将会纳入更多的潜在抗原组分。因此，完成图基因组可一劳永逸地解决该问题。此外，发酵过程不同代次基因组的获得，可同时比较任意所需抗原基因、启动子及基因组结构的遗传稳定性。

本研究结果显示，不同厂家和中国食品药品检定研究院中国医学细菌保藏管理中心保藏的百日咳生产菌株三代测序结果之间具有良好的基因组结构和序列的一致性。百日咳疫苗生产菌株二代测序与三代测序共线性的结果显示，二者之间存在基因组内部易位、倒置等结构性变化。表明产生差异的原因是不同测序方式的原理造成的。相比二代测序，三代测序的结构更能准确地表现百日咳疫苗生产原始菌株基因组的内部结构。

同时随着国内百日咳年度报告数持续维持高位，分子流行病学研究的重要性日益凸显。新的流行株相对于疫苗株的基因变异和结构变异程度以及优势菌株的进化方向成为新的关注点。结构完整，背景清晰的百日咳疫苗生产用菌株完成图基因组无疑会起到重要作用。

本研究为我国百日咳疫苗生产用菌种的质量控制提供了理论依据，也为百日咳杆菌的分子流行病学研究提供了参考。