基于全基因组测序技术筛选儿童金黄色葡萄球菌的致病相关分子标志物

陈建余 王旭麟 李文宇 陈敏琪 周俊立 姚振江 傅锦坚 叶小华

（1.广东药科大学公共卫生学院，广东广州 510310；2.柳州市妇幼保健院，广西柳州 545001）

金黄色葡萄球菌（Staphylococcus aureus，以下简称“金葡菌”）是儿童感染的重要致病菌之一［1］。金葡菌的多重耐药率较高，且社区和医院中的金葡菌交叉传播风险较高，对人类健康造成极大威胁，给全球带来了沉重的疾病经济负担，已成为国际上日益关注的临床问题和公共卫生问题［1］。既往研究表明，儿童人群中的金葡菌鼻腔定植现象较为严重且明显高于成人，在婴幼儿时期的定植率高达45%［2］；同时，儿童身体各方面尚未发育完善，对病原体的抵抗力较低，儿童金葡菌感染率明显高于成人，且一旦感染则更容易发展为菌血症、骨髓炎、坏死性肺炎等侵袭性疾病［3］，提示儿童是金葡菌定植和感染的重点高危人群。既往研究主要局限于研究医院感染人群金葡菌的耐药性、分子特征，而关于健康儿童鼻腔定植金葡菌的分子流行病学研究仍缺乏系统报告，且金葡菌致病相关标志物仍未阐明。金葡菌感染菌株与定植菌株的分子分型具有一定的相似性，仅通过克隆群、多位点序列分型等传统分子分型方法难以区分菌株间的微小遗传变异，有必要通过高通量全基因组测序技术阐明两类金葡菌间全基因组水平上的潜在遗传差异。金葡菌的致病力强弱与其自身含有的毒力因子高度相关，其中与金葡菌感染有关的毒力因子包括肠毒素、溶血素、中毒休克综合征毒素、杀白细胞毒素等。2020 年研究［4］表明，中毒休克综合征毒素基因tsst-1、表皮剥脱毒素基因eta和etb的表达可增强金葡菌的毒力及致病力，fnbA、cna、sdrE、sej、eta、hlg和ica基因在感染菌株中比定植菌株更加普遍，因此这些毒力因子为揭示金葡菌的致病相关标志物提供重要思路。随机森林是一种基于决策树的非参数模型，可以有效处理基因组数据的高维、高相关性等问题，其特点是可以对变量的重要性进行排序和识别致病相关因子，为揭示金葡菌致病相关标志物提供新策略［5］。本研究对512 株金葡菌（272株感染菌株和240株定植菌株）进行全基因组测序分析，全面检测金葡菌的耐药基因和毒素基因，运用随机森林筛选金葡菌致病相关分子标志物，为追溯高致病性金葡菌和开展精准的靶向干预提供遗传学证据。

1 资料与方法

1.1 研究对象

本研究采用横断面研究设计。（1）感染菌株来源：收集2014—2017 年在柳州市妇幼保健院就诊的感染儿童临床标本。感染儿童的纳入标准：年龄≤7 岁；家长签署知情同意书；有临床感染症状（如咳嗽、听诊异常、呼吸困难、发热超过38℃等），且从感染部位（包括痰液、脓液、血液、支气管肺泡灌洗液等）采集的临床样本中分离和鉴定出符合标准的阳性金葡菌，最终确定为金葡菌感染的患儿。排除标准：同一患儿的不同部位采集标本，取主要发病部位的标本；同一患儿对于同一疾病多次采样时，只取首次标本。感染儿童的信息采集包括：年龄、性别、病案号等基本信息；标本来源、采样时间、疾病诊断等临床信息。（2）定植菌株来源：采用分层整群抽样抽取2018 年4—6 月期间广西柳州市6 所幼儿园，开展健康幼儿园儿童的鼻咽拭子采样。健康儿童的纳入标准：年龄≤7 岁的健康儿童，且家长签署知情同意书。排除标准：调查前7 d 内患有肺炎、支气管炎等急性感染性疾病的儿童，以及患有血友病等易出血疾病导致不便采样的儿童。共调查1 702名健康儿童，检出277株金葡菌，其中240株成功进行全基因组测序分析。本研究获得柳州市妇幼保健院医学伦理委员会（批件号：2018-085）和广东药科大学伦理委员会（批件号：2015-20）批准。

1.2 菌株分离鉴定

菌株分离鉴定试验包括：甘露醇氯化钠琼脂培养基显色反应、革兰染色、纯化培养、溶血试验及血浆凝固酶试验，依次进行上述试验，结果阳性则鉴定为金葡菌。

1.3 全基因组测序分析

经菌液制备、试剂盒提取DNA、质量控制后获得金葡菌DNA，经分光光度检测仪检测浓度和纯度及格的基因组DNA 进行全基因组测序，采用Illumina Hiseq 2000 测序仪对基因组DNA 进行双端测序，使用FastQC 进行原始测序数据的质控［6］，采用SPAdes 3.12.0 软件进行序列拼接。使用基因组流行病学中心（Center for Genomic Epidemiology,CGE） 网站（http://www.genomicepidemiology.org/）中的ResFinder 4.1 （https://cge.cbs.dtu.dk/services/ResFinder/）进行序列比对，对15 类2 208 个耐药基因进行检测［7］。使用CGE 中VirulenceFinder 2.0（https://cge.cbs.dtu.dk/services/VirulenceFinder/） 进行序列比对获得菌株的毒素基因，主要包括免疫逃逸基因、胞外酶基因、溶血毒素基因、肠毒素基因、表皮剥脱素基因、杀白细胞毒素基因、中毒休克综合征毒素基因。以S. aureussubsp.aureusMRSA252（accession number: NC_002952）为参考基因组，使用Snippy 软件（https://github.com/tseemann/snippy/）对菌株进行核心单核苷酸多态性鉴定，使用Gubbins软件去除基因重组或水平基因转移相关区域［8］。使用Fasttree 软件构建基于核心基因组的系统进化树，采用广义时间可逆替代（general time reversible, GTR）模型［9］，估计方法为极大似然估计法，Boostrap 值设置为100。使用Chiplot 网站（http://www.evolgenius.info/）对系统进化树进行可视化和美化。

1.4 统计学分析

采用EpiData 3.0 软件进行数据的双人双录入，并通过一致性检验保证数据的准确性。采用Stata 16.0 和R 4.1.2 软件对数据进行统计学分析。计数资料以频数和百分率（%）表示，采用Pearsonχ2检验或Fisher确切概率法比较金葡菌感染菌株与定植菌株的耐药基因、毒素基因携带率，从而初步筛选致病相关因子。P＜0.05 为差异有统计学意义。采用单因素logistic 回归或确切logistic 回归模型估计OR及95%CI。运用机器学习方法中的随机森林进一步筛选重要的致病相关特征变量，模型以菌株类型（1=感染菌株，0=定植菌株）作为结局变量，以分子特征（毒素基因和耐药基因）作为自变量。随机森林分析中，通过randomForestSRC 程序包中随机森林的变量捕获法（variable selection using random forests, VSURF）对重要特征变量进行筛选［10］；采用randomForest程序包对最终模型中的变量进行重要性评分，特征变量的重要性以平均基尼指数减少量（mean decrease in the Gini, MDG）排序，从而判断各个变量在模型中的重要性。采用组内回代和十折交叉验证方法评价最终模型的预测效果，评价指标包括正确率、灵敏度、特异度、阳性预测值、阴性预测值、受试者操作特征曲线（receiver operating characteristic curve, ROC曲线）。

2 结果

2.1 基本情况

有272株金葡菌（即感染菌株组）来源于医院感染儿童，年龄为0～7 岁，中位数为2 个月，172株（63.2%）来源于男性、100株（36.8%）来源于女性。有240 株鼻腔定植金葡菌（即定植菌株组）来源于健康儿童，年龄为3～7 岁，中位数为5 岁，124 株（51.7%）来源于男性、116 株（48.3%）来源于女性。

2.2 定植菌株组与感染菌株组的毒素基因携带情况

金葡菌感染菌株组的seb、sep、splA、splB、splE、edinC、lukD、lukE、lukF-PV、lukS-PV、eta和etb基因携带率高于定植菌株组（均P＜0.05）；但感染菌株组的sec、sec3、seg、seh、sei、sel、sem、sen、seo和seu基因携带率低于定植菌株组（均P＜0.05），其他毒素基因携带率比较差异无统计学意义（P＞0.05）。见表1。

表1 不同类型金葡菌的毒素基因携带情况比较 ［株（%）］

2.3 定植菌株组与感染菌株组的耐药基因携带情况

金葡菌感染菌株组的lnuA、aadD、tetK和dfrG基因携带率明显高于定植菌株组（P＜0.05）。但感染菌株组的blaTEm-1A基因携带率低于定植菌株（P＜0.05）。其余耐药基因携带率比较差异均无统计学意义（P＞0.05）。见表2。

表2 不同类型金葡菌的耐药基因携带情况比较 ［株（%）］

2.4 随机森林筛选金葡菌致病相关分子特征

基于核心单核苷酸多态性的系统进化树（图1）表明：感染菌株与定植菌株在系统进化树上的分布没有明显差异，提示两者可能具有相似的遗传背景，需进一步揭示致病相关分子特征。以菌株类型（1=感染菌株，0=定植菌株）为因变量，建立随机森林模型筛选儿童金葡菌致病相关的重要特征变量（69 个），采用VSURF 法最终筛选出16 个特征变量。对变量筛选前后的模型进行预测效果评价，结果表明：变量筛选前后预测模型的交叉验证正确率分别为69%、68%（表3）；变量筛选前模型的曲线下面积（area under the curve,AUC）稍高于筛选后（0.75 vs 0.70）（图2A～B）。最终筛选出的16 个特征变量为肠毒素基因（sem、sep、ser、sea）、表皮剥脱毒素基因etb、胞外酶编码基因splE、β-内酰胺类耐药基因（mecA、blaZ、blaTEm-1A）、大环内酯类耐药基因（ermA、ermB、lnuA）、氯霉素类耐药基因cat(pC233)、氨基糖苷类耐药基因［aph(3')-Ⅲ、ant(9)-Ⅰa、ant(6)-Ⅰa］；变量重要性排序结果显示，前5个最重要的特征变量分别为sem、etb、splE、sep、ser（图2C）。使用随机森林模型进一步对16 个重要致病相关标志物进行风险预测，结果显示：有10 个致病相关标志物位于对角线上方，提示携带这些致病相关标志物会增加金葡菌的致病风险；而有6个致病相关标志物位于对角线下方，提示携带这些致病标志物会降低金葡菌的致病风险（图2D）。

图1 512 株金葡菌的系统进化树 从内环到外环的变量分别为菌株来源（红色为感染菌株，蓝色为定植菌株）、毒力基因数、耐药基因数、16个致病相关特征变量（深红表示有，浅红表示无）。

图2 随机森林模型的拟合效果和风险得分图 A：变量筛选前模型的ROC曲线。B：变量筛选后模型的ROC曲线。C：最终模型中16个特征变量的重要性排序。D：最终模型中16个特征变量的风险得分图。

表3 基于分子特征的随机森林模型预测效果评价

3 讨论

关于金葡菌的耐药机制研究较多，常见于以下几类耐药机制［11-14］：（1）酶降解抗生素，如β-内酰胺酶；（2）产生拮抗物质，减少抗生素的有效浓度；（3）改变细胞靶标，使抗生素不能正常发挥作用；（4）主动外排泵作用。blaZ基因受blaR1-BlaI 系统的调节控制，blaR1 蛋白受到β-内酰胺类抗生素的刺激后，导致抑制蛋白BlaI水解脱离结合位点，编码产生β-内酰胺酶，继而破坏β-内酰胺从而使金葡菌对青霉素类抗生素耐药［11］。mecA编码PBP2a 蛋白，产生过量的PBP 蛋白消耗β-内酰胺类药物，从而导致金葡菌对β-内酰胺类药物耐药。既往研究［12］显示，虽然blaTEm基因编码的β-内酰胺酶主要通过丝氨酸残基灭活β-内酰胺环产生耐药性，但是blaTEm在多重耐药菌株和敏感菌株中均被检出且差异无统计学意义，提示blaTEm可能不是导致其广泛耐药的主要原因；值得注意的是，本研究的随机森林预测模型揭示，blaTEm-1A基因是金葡菌致病相关标志物，提示其可能影响菌株的致病能力。erm编码核糖体甲基化酶催化金葡菌23S rRNA 发生甲基化反应使核糖体靶位点发生变化，进而阻止大环内酯类药物与核糖体靶位结合，从而导致金葡菌对大环内酯类药物耐药［13］。lnuA基因编码核苷酸转移酶，从而介导对林可酰胺类抗生素的耐药性［14］。

细菌的毒力因子在感染性疾病中起着关键作用，它是由多种毒素基因决定的。金葡菌具有数十种毒素因子，包括肠毒素、表皮剥脱毒素、胞外酶等。本研究随机森林结果显示，金葡菌致病相关的毒素基因有肠毒素基因（sem、sep、ser和sea）、表皮剥脱毒素基因（etb）、胞外酶编码基因（splE）。葡萄球菌肠毒素在食源性中毒中最常见，它不仅可引起食物中毒，还可引起中毒性休克综合征等严重的侵袭性疾病［15］。本研究中，感染菌株组sep基因的携带率明显高于定植菌株组，但是sem基因的携带率明显低于定植菌株组，提示感染菌株与定植菌株存在毒力差异，与中国台湾的一项研究结果一致［16］。表皮剥脱毒素A 和表皮剥脱毒素B是导致大多数人类葡萄球菌烫伤样皮肤综合征的重要原因［17］；本研究中，感染菌株组较定植菌株组的eta（7.4% vs 2.9%）、etb（9.2% vs 2.5%）基因携带率高，以上结果提示表皮剥脱毒素与菌株的致病能力密切相关。splE编码丝氨酸蛋白酶样蛋白，可以在健康人，尤其是哮喘患者中诱导Ⅱ型超敏反应［18］；西班牙一项研究提示，splE基因可作为金葡菌菌血症的致病相关标志物［19］。以上结果提示不同类型疾病可能存在致病相关标志物的差异，因此后续研究可针对不同疾病类型进一步探讨潜在的特异性标志物。

基因组数据具有“高维度、小样本”的数据特点，且存在“非线性、高相关”的复杂关联，使得传统分析方法不再适用，这对基因组关联分析提出了新挑战。随机森林是利用多棵决策树对样本进行训练和预测的机器学习算法，能有效处理线性、非线性、具有交互作用的复杂数据，具有抗噪声、防止过拟合、不受共线性影响的重要特征。随机森林不仅具有良好的预测能力，且具有精度高、稳定性好、易操作等优点，因此已成为处理医学领域基因组数据的精确机器学习方法之一。为了获得更准确的风险预测模型，本研究以金葡菌的常见致病相关基因（32 个耐药基因和37 个毒素基因）为预测变量、菌株类型（感染和定植菌株）为因变量建立随机森林模型，最终模型的16个预测因子分别为肠毒素基因（sem、sep、ser和sea）、表皮剥脱毒素基因etb、胞外酶编码基因splE、β-内酰胺类耐药基因（mecA、blaZ和blaTEm-1A）、大环内酯类耐药基因（ermA、ermB和lnuA）、氯霉素类耐药基因cat(pC233)、氨基糖苷类耐药基因［aph(3')-III、ant(9)-Ⅰa 和ant(6)-Ⅰa］，最终模型的交叉验证正确率为68%，AUC 为0.70，提示模型的拟合效果较好。随机森林模型揭示：前5 个最重要的预测变量为sem、etb、splE、sep、ser；其中sep（OR=3.97）、etb（OR=3.95）、ser（OR=1.68）和splE（OR=1.68）明显提高金葡菌的致病风险，提示肠毒素、表皮剥脱毒素和胞外酶是金葡菌的重要致病相关标志物。以上结果揭示这些致病相关分子标志物有较大的潜力预测金葡菌致病株。

综上，本研究使用高通量全基因组测序技术全面检测金葡菌基因组的耐药基因和毒素基因，采用随机森林筛选出金葡菌的16 个致病相关标志物（6 个毒素基因和10 个耐药基因），且模型预测效果较优，为追溯高致病性金葡菌和开展精准的靶向干预提供遗传学证据。本研究也存在一定局限性：（1）本研究提示致病相关特征变量与结局之间存在统计学关联，但筛选的致病相关标志物与疾病状态之间不一定存在因果关系，且仍未能深入阐明潜在的生物学功能与作用机制，因此今后可以结合生物学功能富集分析、通路富集分析对其潜在作用机制进行深入研究。（2）本研究于2014—2018 年间对金葡菌感染患儿和健康儿童进行金葡菌采样，菌株的采集时间较早；但是，本研究主要从基因组水平上探究金葡菌的致病相关标志物，且既往研究表明菌株分子特征与菌株的采样时间无关［20］。今后可进一步开展多中心调查来增大样本量和补充新菌株。

利益冲突声明：所有作者均声明不存在利益冲突。