金黄色葡萄球菌在万古霉素中的表型可塑性

郑旭阳,戎梦迪,白 珺,何晓青,谢响明,金 一

(北京林业大学生物科学与技术学院，北京 100083)

表型可塑性指同一基因型在不同环境下产生不同表型的能力[1].表型可塑性不仅是物种适应性进化的一个重要方面，也是选择进化的产物.表型可塑性在自然界中普遍存在，有关植物和动物的表型可塑性已有大量文献报道，但对微生物表型可塑性的研究较少[2-3].微生物在自然环境中分布广泛，对维持生态系统的稳定起着重要作用.环境变化对微生物造成的影响十分复杂，表型可塑性在一定程度上抵消了环境变迁对微生物生长所造成的影响[4].对微生物表型可塑性的研究有助于进一步解析微生物群落特征及功能稳定性[5].金黄色葡萄球菌(Staphylococcusaureus)作为一种常见病原菌能够引起各种感染，通常使用抗生素对金黄色葡萄球菌引起的感染进行治疗[6]；但由于抗生素的过度使用，金黄色葡萄球菌耐药性愈加严重[7].菌株在不同抗生素胁迫下表现出不同程度的可塑性，推动其适应性进化，进一步发展出耐药菌株.从野生动物和动物园动物分离的金黄色葡萄球菌菌株具有不同的耐药性，且在基因分型中表现出较强的多样性[8].研究表明，根据奶牛中分离的金黄色葡萄球菌的基因型差异可以预防疾病的发生[9].但目前有关金黄色葡萄球菌在抗生素环境中的表型可塑性尚未见报道.

针对表型可塑性的研究主要是将基因型与表型进行关联分析，从而探究基因在个体表型中起的作用[10].全基因组关联分析(genome-wide association analysis,GWAS)能够探测同一物种中不同个体全基因组范围内的基因变化，通过基因型与表型性状关联筛选相关的单核苷酸多态性(single nucleotide polymorphism,SNP)位点[11].GWAS为探索基因或染色体区域功能提供了新方法，使得大量功能未知的基因及区域被发现，为人类攻克疾病提供了新的方向[12-13].GWAS应用在微生物方面不仅能帮助人们了解微生物的结构和功能，还能进一步确定其与环境之间的作用关系[14]，同时能够快速识别调节表型变异的遗传因素，更容易鉴别出导致某种表型变化的遗传基础[15].目前的研究多为单一表型变量与基因型进行关联分析，使用多个相关变量对应一种基因型的分析较少.双变量GWAS可将2种表型在同一时间进行分析，与单变量相比，其统计效力加大，参数估计精确性提高[16].前人应用GWAS分析了金黄色葡萄球菌在单独培养和混合培养条件下的生长情况，定位到与金黄色葡萄球菌生长显著相关的SNP位点，表明GWAS 在研究细菌互作进化机制方面的潜力[17].本研究主要应用双变量GWAS的方法对金黄色葡萄球菌在2种不同浓度万古霉素中的表型可塑性基因进行定位，为研究金黄色葡萄球菌在抗生素压力下的耐药性生长提供依据，同时为应用双变量GWAS研究细菌表型可塑性提供参考.

1 材料与方法

1.1 菌株来源及其培养方法

本研究所用菌株见表1.41株原始金黄色葡萄球菌(编号为S1-S9、S11-S42)均对万古霉素敏感，经脑心浸液琼脂(brain heart infusion agar,BHIA)培养基复苏，挑取单克隆细菌，用0.5～4.0 μg·mL-1脑心浸液肉汤(brain heart infusion,BHI)—万古霉素选择性固体培养基测定其最低抑菌浓度(minimum inhibitory concentration,MIC)值.将菌液用生理盐水稀释至1.5×108个·mL-1左右，接种于含1/2 MIC万古霉素的BHIA平板上，37 ℃恒温培养24 h.将生长的单菌落接种于同浓度BHI-万古霉素平板上，每个药物浓度培养4 d，每2 d转接至同浓度的平板，每4 d测定一次细菌的万古霉素MIC，根据MIC决定下次接种的万古霉素平板浓度，连续培养60 d.采用E-test试纸条法再次验证所得菌株的MIC，并将诱导获得的菌株编号为S1′-S9′、S11′-S42′[18].所有菌株保存于-80 ℃.

1.2 试剂

BHI培养基、BHIA培养基均为英国Oxiod公司产品.细菌基因组提取试剂盒购自天根生化科技有限公司.万古霉素、溶葡萄球菌素购自美国Sigma公司.

1.3 生长曲线及最大生长速率的测定

将41株金黄色葡萄球菌平板划线，37 ℃恒温培养12 h进行菌株复苏.挑取单菌落在37 ℃、130 r· min-1过夜培养，测定菌液的D600 nm值.将菌液接种于20 mL不含抗生素或含4 μg·mL-1万古霉素的BHI培养基中，初始菌液浓度为5×103个·mL-1，培养48 h.每株菌设置3个重复.共设置14个时间点，分别为1、2、4、6、8、10、12、16、20、24、30、36、42、48 h，使用酶标仪(Infinite M200 PRO,TECAN，瑞士)检测其D600 nm值，应用GraphPad Prism 6.0软件绘制生长曲线，并计算其最大生长速率.

1.4 全基因组测序

取2 mL过夜培养的菌液，使用溶葡萄球菌素和细菌基因组提取试剂盒提取其基因组，检测合格后用Illumina HiSeq 4000进行全基因组重测序(北京奥维森基因科技有限公司).基因组检测、建库、测序等环节都可能对数据质量和数量产生影响，进而影响后续分析结果.因此对原始数据进行质量控制，包括去除低质量读数和接头，计算测序错误率，统计Q20、Q30、GC含量等，得到有效数据.选择金黄色葡萄球菌亚种NCTC 8325为标准菌株进行比对，得到BAM格式的比对结果.运用SAMtools对比对结果进行排序，标记重复序列，对得到的SNPs进行过滤.

1.5 GWAS

以41株金黄色葡萄球菌在无抗生素培养基和含4 μg·mL-1万古霉素培养基中的生长数据为表型，结合全基因组测序数据进行关联分析.为检验SNP与2个表型的相关性，在GWAS中群体分层被视为随机效应.本研究将基于单变量GWAS线性混合模型开发的多变量线性混合模型用于双变量GWAS，借助GEMMA软件的多变量线性混合模型计算P值，将多个表型同时与遗传标记相关联[10].分析结果经1 000次重组排列，在任意一次排列中，基因型和多变量表型之间的关系被随机重组，将这1 000个P值由小到大排列，以0.05作为显著性的评价标准，将第50个值作为GWAS的阈值，确定每个时间点显著的SNPs.显著的SNPs显示在由软件R3.3.0绘制的曼哈顿图中.

表1 细菌来源1)Table 1 Sources of bacteria

1)S1-S9、S11-S42为原始菌株；S1′-S9′、S11′-S42′为本研究中的菌株；CGMCC为中国普通微生物菌种保藏中心，CFCC为中国林业微生物菌种保藏管理中心，CICC为中国工业微生物菌种保藏中心，CMCC为中国医学微生物菌种保藏管理中心，ACCC为中国农业菌种保藏中心，CCTCC为中国典型培养物保藏中心，CPCC为中国药学微生物菌种保藏管理中心，CAU为中国农业大学动物医学院；MIC为菌株经60 d诱导后在万古霉素压力下的最低抑菌浓度；序列号来自http://bigd.big.ac.cn.

2 结果与分析

2.1 生长曲线

如图1所示，41株菌在无抗生素压力下生长时，除S6′和S8′，其他菌株均呈现明显的S型曲线.26株菌培养6 h左右(时间点4)进入对数生长期，其余15株菌在培养10 h左右(时间点6)进入对数生长期，前者对数生长期比后者略陡.在含4 μg·mL-1万古霉素的条件下生长时，菌株进入对数生长期的时间分散.10株菌(S7′、S8′、S9′、S15′、S16′、S28′、S32′、S35′、S39′、S41′)在2种培养条件下的生长曲线基本重合，说明菌株生长未受到万古霉素的影响；10株菌(S3′、S11′、S13′、S20′、S23′、S33′、S34′、S37′、S38′、S40′)在万古霉素影响下对数期生长略缓慢，但迟滞期持续时间未变；8株菌(S5′、S18′、S19′、S21′、S26′、S27′、S29′、S31′)在万古霉素影响下进入对数生长期耗时较长且生长十分缓慢，说明万古霉素影响了菌株的生长；其余13株菌在含4 μg·mL-1万古霉素的培养基中没有生长.

2.2 最大生长速率

将无抗生素压力下菌株的最大生长速率按照由大到小排列，如图2所示.7株菌(S7′、S8′、S15′、S26′、S27′、S28′、S32′)在2种培养条件下的最大生长速率基本一致；9株菌(S16′、S19′、S23′、S33′、S34′、S35′、S37′、S40′、S41′)在含万古霉素培养基中的最大生长速率略低于无抗生素培养基；23株菌(S1′、S2′、S3′、S4′、S5′、S6′、S11′、S12′、S13′、S14′、S17′、S18′、S20′、S21′、S22′、S24′、S25′、S29′、S30′、S31′、S36′、S38′、S42′)的最大生长速率明显降低；而S9′和S39′在含万古霉素培养基中的最大生长速率略高于无抗生素培养基.由此可知，S7′、S8′、S15′、S26′、S27′、S28′、S32′、S9′和S39′这9株菌较其他32株菌更具表型可塑性.

图1 41株金黄色葡萄球菌生长曲线Fig.1 Growth curves for 41 strains of S.aureus

续图1Continued Fig.1

图2 41株金黄色葡萄球菌最大生长速率Fig.2 Maximum growth rates for 41 strains of S.aureus

2.3 全基因组测序结果

金黄色葡萄球菌重测序结果见表2.每株菌的有效数据量在941～2 581 Mb之间.本研究所用数据保存在生命与健康大数据中心的组学原始数据归档库中，序列号为CRA001125.41株金黄色葡萄球菌共得到61 869个SNPs.

2.4 GWAS结果

对41株金黄色葡萄球菌在2种培养条件下的生长情况进行监测，通过GWAS确定表型与基因型的关联.金黄色葡萄球菌的61 869个位点中，387个位点位于阈值之上.时间点7得到的显著位点最多，-logP值最大的位点位于时间点4.将14个时间点的显著SNPs按照-logP值由大到小排列，筛选出现频率高且-logP值大的SNP，最终在时间点3、4、13分别得到3个显著SNPs(表3、图3).其中，-logP值最大的位点382 501位于SAOUHSC_00375上，在所有时间点中共出现5次，SAOUHSC_00375编码一磷酸鸟苷合酶；位点2 535 773共出现3次并在时间点4出现最大值，其所在基因编码谷氨酸合酶α亚基，参与谷氨酸合成过程；位点264 897所在基因与药物转运蛋白有关；位点615 798和1 842 166分别位于编码H+反向转运蛋白亚基A和Ⅰ型限制修饰系统亚基的基因上；位点965 494所在基因与双功能自溶素有关.

表2 金黄色葡萄球菌重测序数据1)Table 2 Resequencing data for S.aureus

1)插入片段300 bp指建库时构建的文库片段大小;读长150∶150指采用双端测序,两端各测150 bp;原始数据产量为所测序列的个数乘以所测序列的长度；有效数据量为原始数据产量过滤后序列的个数乘以序列的长度；Q20含量为碱基正确识别率在99%以上的碱基数目百分比(测序仪在碱基识别时，会给出每个碱基正确识别的概率)；Q30含量为碱基正确识别率在99.9%以上的碱基数目百分比.

3 讨论

由41株金黄色葡萄球菌的MIC结合生长表型数据分析可得，在含万古霉素的培养基中生长时，大部分菌株随抗生素浓度的增大，迟滞期时间延长，指数期坡度变缓，最大生长速率降低，且达到稳定期的最大D600 nm值略低于无抗生素培养基.部分菌株(S1′、S17′、S25′、S42′)的MIC均在4 μg·mL-1以上，但在4 μg·mL-1万古霉素培养基中未生长；部分菌株(S18′、S23′、S29′、S33′)的MIC值均为3 μg·mL-1，但在4 μg·mL-1万古霉素培养基中均有生长，可知这些菌株在万古霉素环境中适应性更好.另外，同一MIC值的金黄色葡萄球菌生长情况也有差异，可知MIC的大小与金黄色葡萄球菌的表型可塑性无直接关联.

表3 9个候选基因和蛋白Table 3 Nine candidate genes and proteins

A、B、C分别代表时间点4、3、13.每一个点代表一个SNP.图3 9个显著SNPs曼哈顿图Fig.3 Manhattan maps for 9 significant SNPs

本研究应用双变量GWAS研究金黄色葡萄球菌在万古霉素中的表型可塑性，最终得到9个可能与其表型可塑性有关的SNPs，通过基因注释发现这些SNPs主要与氨基酸形成及离子转运等有关.在培养6 h出现-logP值最大的SNP位点382 501，该位点位于基因SAOUHSC_00375上，此基因编码一磷酸鸟苷合酶；在使用抗生素治疗奶牛乳腺炎时，检测到SAOUHSC_00375的表达对细菌存活具有重要作用，并在感染期间可能帮助识别药物靶标[19].SNP 2 535 773所在基因SAOUHSC_02760编码谷氨酸合成酶亚基，参与铵离子与有机化合物的结合——生产氨基酸的重要步骤[20].SNP 738 836位于非编码区，在所有位于非编码区的SNP中-logP值最大，出现在培养后4 h，非编码区不编码蛋白质，但可以调节生物合成过程中基因的表达.SNP 264 897所在基因编码转运蛋白，离子转运能够维持金黄色葡萄球菌菌体中的抗生素浓度，以促进其在抗生素中的生长，可能与金黄色葡萄球菌的耐药机制有关[21].SNP 1 394 043在培养4 h出现，-logP值为30.10，此位点所在基因功能仍有待发掘.SNP 615 798所在基因SAOUHSC_00625编码H+反向转运亚基，在金黄色葡萄球菌面对环境压力时，能够通过离子的转运调节菌体内的渗透压和pH值等[22].SNP 1 842 166编码Ⅰ型限制修饰系统M亚基，限制性系统修饰的甲基转移酶能够维持序列的甲基化方式并切割不同的外源DNA以保持较低水平基因的转移效率[23].SNP 965 494所在基因SAOUHSC_00994(atl)编码双功能自溶素，金黄色葡萄球菌可产生N-乙酰-胞壁酰-L-丙氨酸酰胺酶和内切-Δ-N-乙酰氨基葡萄糖苷酶2种不同的肽聚糖水解酶，从而参与金黄色葡萄球菌的自溶、细胞壁更新、细胞分离等过程[24].SNP 1 871 317所在基因为SAOUHSC_01966，编码假定蛋白.

综上可知，氨基酸合成、离子转运等相关基因很可能参与菌株的可塑性调控.另外，MIC与金黄色葡萄球菌在万古霉素中的表型可塑性没有直接关联.