Kachroo P, Eraso J M, Beres S B, et al. 2019. Integrated analysis of population genomics, transcriptomics and virulence provides novel insights into Streptococcus pyogenes pathogenesis. Nature Genetics, .
【非常牛的案例文章！隆重推荐！】

群体基因组学，转录组学和毒力的综合分析为化脓性链球菌发病机制提供了新的见解

摘要

化脓性链球菌每年在全世界造成7亿人感染，尽管经过一个世纪的努力，但没有针对这种细菌的许可疫苗。尽管已经发表了许多关于细菌病原体的大规模基因组研究，但大型细菌群体中基因组，转录组和毒力之间的关系仍然知之甚少。我们对2,101个emm28化脓性链球菌侵染菌株的基因组进行了测序，我们从中选择了492个系统发育不同的菌株用于转录组分析，50个菌株用于毒力评估。数据整合提供了对该模式生物的毒力机制的新颖理解。全基因组关联研究、表达数量性状基因座分析、机器学习和同基因突变菌株鉴定并证实了基因间区域中的单核苷酸插入缺失，其显着改变全局转录物谱和最终毒力。我们使用的综合策略通常适用于任何微生物，并可能导致许多人类病原体的新疗法。

无论生态位还是宿主范围，所有细菌种类都包含遗传多样性菌株。人们对微生物分子遗传学知之甚少的领域是大量传染性菌株的基因组，转录组和毒力之间复杂的相互作用。遗传变异可能会影响基因转录水平，但这种情况的真实程度及其对发病机制的影响仍不明确【good！】。尽管大型基因组学研究已经发表1，但在比较转录组领域的研究却少得多7和涉及自然种群的毒力分析10,11。此外，除了一项涉及相对较小的大肠杆菌菌株样本12的研究外，尚未研究基因组，转录组和毒力之间的关系。对基于种群的菌株样本的各种数据进行综合分析可能会对疫苗制剂，新疗法和诊断以及公共卫生等领域的转化研究工作产生影响。
化脓性链球菌（A群链球菌（GAS））是一种严格的人类病原体，每年在全球儿童和成人中引起超过7亿的感染13。人类感染的严重程度从相对轻微的咽炎（“链球菌性咽喉炎”）到极端严重和危及生命的感染，如败血症和坏死性筋膜炎/肌炎，通常被称为“肉食”疾病。该生物体还会引起皮肤和软组织感染，并导致感染后遗症，如风湿热和风湿性心脏病，这是全球发病率的重要原因[13,14]。
GAS已被用作模型生物，用于研究菌株类型和疾病表型之间的关系，以及流行病学[1,6,15-7]。 emm28菌株是美国和许多欧洲国家中与侵袭性GAS感染相关的前五种emm类型之一18？3。由于仍然无法解释的原因，一些emm类型或M蛋白血清型的菌株与特定类型的人类感染非随机相关17,24？0。例如，emm28 GAS菌株在产褥期败血症（儿童发烧）和新生儿感染17,31？4的病例中反复出现。
尽管在所选生物的基因组学方面取得了重大进展，但对于共享最近共同祖先的细菌菌株的克隆相关后代的转录组多样性的性质和程度知之甚少。关于这个问题的数据对于增强对自然群体中细菌进化，表型多样化和微生物流行的理解至关重要。为了解决这些知识gap，我们对基于综合种群的研究中回收的2,101种emm28 GAS菌株的基因组进行了测序，并使用所得的系统发育信息选择代表性菌株来分析转录组（n = 492株）和毒力（n = 50株）。数据整合提供了对该模型生物的生物学的新的理解，包括在相对密切近缘的生物进化枝中具有惊人的转录组变异幅度。统计学方法和机器学习的应用促进了一种新的分子遗传过程的发现，该过程支持一些GAS菌株中增强的毒力。

结果

群体结构和时间分布。

我们对在北美和欧洲的六个国家中分离的2,101个emm28 GAS菌株的基因组进行了测序，这些菌株在26年期间：1991年至2016年（表1，补充表1和补充图1）。所有菌株都作为综合群体研究的一部分被回收。将基因组测序至202倍平均覆盖度（补充图2a，补充表1和方法）35。遗传关系的推断是通过使用核心基因组中存在的SNP（即没有移动遗传元件（MGE）的基因组，例如前噬菌体和整合 - 结合元件（ICE））来进行的（补充表1）。主要的emm28 GAS群体通过贝叶斯聚类分布到两个主要进化枝（进化枝1和2）和四个进化枝（指定为SC1A，SC1B，SC2A和SC2B）中（图1a，b）。进化枝2生物与进化枝1的部分区别在于28.0 kb水平基因转移（HGT）区块，其贡献520个核心SNP，以及位于该HGT区块外的19个核心SNP（补充表2）。这个28.0-kb的HGT区块包括编码分泌的毒素NAD + - 糖水解酶（SPN）和链球菌溶血素O（SLO）——已知GAS毒力的主要贡献者6,16,36？9的nga–ifs–slo操纵子 。 ifs编码SPN40的内源性抑制剂。重要的是，重组获得nga–ifs–slo操纵子的高表达变体可以增加在灵长类动物上呼吸道的存活率，增强毒力，并引发洲际流行病[1,6,16,41]。进化枝2生物具有与粪链球菌亚种equisimilis42中存在的类似三基因操纵子具有99％序列同一性的nga–ifs–slo区域，并且可能通过重组事件被亚类2A 2A GAS获得。
SC1B包含的菌株最多，占菌株的49.7％，其次是SC2A（26.8％），SC1A（22.3％）和SC2B（0.53％）。属于子SC1A，SC1B和SC2A的菌株根据年份，地理位置和MGE含量而变化（图1c和补充图3）。 SC1A菌株的显着时间位移同时发生，美国SC2A菌株激增，其中约55％的菌株是SC2A。 SC1B菌株主要分离自芬兰，挪威，法罗群岛和冰岛的患者，而SC1A是加拿大流行的亚类。评估群组中的MGE多样性（补充说明和补充表3？），并且20种最丰富的MGE-50基因型占菌株的90％（图1c）。
接下来，我们使用综合策略从种群角度研究动物感染模型中基因组变异，转录组变化和毒力差异之间复杂的相互作用。

转录组签名和群体结构。

为了确定基因表达的不同模式是否与emm28群体结构非随机相关，我们首先对50个菌株的子集进行转录组RNA测序（RNA-seq）分析，这些菌株在遗传上代表三个数值显着的亚类（即SC1A，SC1B和SC2A）（补充表7）。根据方法中描述的标准，从2,095个M28菌株的主要样品中选择菌株进行RNA-seq分析。 50株来自不同国家，国家和地区，具有不同的MGE含量。 RNAseq分析在中指数和早期 - 静止生长期一式三份（三个生物学重复）进行（方法，补充图2b和补充表7）。
尽管50种菌株的基因组背景不同，但分析的菌株数量加上一式三份样品中转录水平之间极高的相关系数（补充图4），允许鉴定与种群结构有关的不同转录组改变（图1）。我们发现了两个出乎意料地具有“异常”转录组的菌株（图2a，b）。手动检查这两种异常菌株的基因组序列数据，在covS全局调控基因中发现了两个独立的大缺失事件（图2a，b）。编码全局转录调节因子如CovR / CovS，RopB和Mga的基因的突变可以改变转录组的显着比例（5~5％）43。
对于在所有已知主要调节基因中具有野生型等位基因的46个菌株（图2c，d），差异表达基因的最大数目发生在分配给不同亚类的菌株之间（补充图5）。在midexponential阶段，当我们比较SC2A菌株与SC1A和SC1B菌株的转录组（分别为32（2％）和15（0.9％））时，观察到最大数量的差异表达基因（补充图5a）和补充表8）。当比较三个遗传子片段的早期静止转录组时，类似的模式是明显的，但差异表达基因的数量相当大（5？倍）（补充图5a和补充表8）。与SC1A和SC1B菌株相比，SC2A菌株具有最大数量的差异表达基因（分别为318（19.9％）和83（5.2％））（补充图5a和补充表8）。
大部分差异表达的基因位于水平转移的28.0-kb区域（HGT区域）并包括nga–ifs–slo操纵子（补充表8）。在midexponential阶段，该HGT区域内基因（28个基因）的35.7％（SC2A与SC1A比较：P <0.0001）和25％（SC2A与SC1B比较：P <0.001）差异表达（P值通过Fisher精确检验评估） ，而早期 - 静止期差异表达21.4％（SC2A与SC1A比较：P = 0.81）和39.3％（SC2A与SC1B比较：P <0.001）差异表达。与SC1A和SC1B菌株相比，SC2A菌株中三个最强烈上调的基因是nga，ifs和slo，在指数期间转录水平增加约4倍，在早期稳定期增加约8倍（补充表8和补充图5b）。
近年来，SC1B菌株引起的感染在包括美国，芬兰，冰岛和挪威在内的一些国家有所增加，而SC1A菌株已大幅减少（补充图3a），从而提高了SC1B菌株可能进化为比SC1A菌株更合适。因此，我们检查了核心染色体上的遗传差异（补充表9），区分了所有SC1A和SC1B菌株，发现所有SC1B菌株在RivR中都有两个连续的非同义突变——grab基因的一个负调节因子（蛋白-G相关的α2-巨球蛋白结合）。 在中期指数期，与SC1A和SC2A菌株相比，grab是SC1B菌株中唯一上调的基因（补充表8和补充图5c）。类似地，与SC1A相比，在SC1B和SC2A菌株的早期 - 稳定期观察到更高的grab转录物丰度（补充图5c）。 GRAB是一种结合α2-巨球蛋白的细胞表面锚定蛋白，允许它在GAS表面保留蛋白水解活性形式的半胱氨酸蛋白酶SpeB，用于保护GAS免受抗菌肽LL-37的杀伤（参考文献48,49），并有助于小鼠模型中的侵袭性感染50。我们不能排除导致subclade位移的随机过程。

验证单体（RNAtag-seq）分析【无重复】。

传统上，通过使用生长至两个不同生长期的菌株的一式三份生物复制品进行细菌的转录组分析。然而，这种方法目前在研究数百种菌株方面经济上不可行。鉴于RNA-seq的准确性，灵敏度和重现性得到改善，我们假设使用单一菌株进行大规模转录组分析（即缺乏重复），与病原体的最佳测序深度（5？0百万测序读数）一致基因组大小约为2Mb（参考文献51），将大大增强对一组相对密切相关的菌株的转录组景观的理解。为了检验我们的假设，我们使用RNAtag-seq52来增加转录组分析的菌株通量。发现来自没有（单体）和具有生物学重复的菌株的表达数据是高度相关的（补充图6和补充说明）。因此，我们进行了通过k均值聚类统计策略（方法）选择的442个遗传代表性和多样化单一菌株（补充表10和补充图7）的群体转录组分析。

不同菌株的群体转录组分析。

为了检验442个单一菌株的转录组之间的关系，我们首先对归一化的表达数据使用主成分分析（PCA）并鉴定了两个主要的簇，称为簇A（n = 339）和簇B（n = 83）。 （图3a）。基于密度的具有噪声的应用的空间聚类验证了这两个聚类的存在（补充图8a）。野生型样菌株（生物信息学评估的菌株具有所有主要调节基因的野生型等位基因）主要与簇A相关，除了10个异常菌株（图3a）。通过使用Pilon（方法）重新检查这些异常菌株的基因组数据，在这10个菌株中的8个中鉴定了covS全局转录调控基因中未检测到的插入缺失。因此，转录组引导的多态性发现确定了这些异常转录组的潜在遗传原因。
考虑到8种covS突变体菌株的转录组与野生型菌株显着不同，并且与之前的结果53？8一致，我们假设在特定主要全局调节因子中具有突变的菌株可能具有可被利用的基因表达的不同潜在模式以区分特定类别的调节基因突变体。为了验证这一假设，我们使用随机森林机器学习59来确定四种类别标签中的一种（即野生型，covR，covS或ropB突变体）是否可以高概率地分配给异常菌株。简而言之，在分析283个单一菌株的转录组谱（参见方法）的基础上，使用随机森林分类来预测8个异常菌株的类别标记。基于转录组的分类正确地将所有八种生物鉴定为covS突变菌株（补充表11）。在81个covRS和21个ropB菌株中，分别有85.2和61.9％的菌株被准确分类（补充表11）。 covRS菌株被错误分类为野生型表型（转录物谱）类似于野生型菌株，与群A菌株分组（图3c）。因此，转录物谱的机器学习分类准确地预测了基因型（调节基因突变状态）并预测了具有主要调节基因突变的菌株的转录物表型（突变样或野生型）。

调节基因突变和转录组变化。

将异常菌株重新分配为covS突变体导致簇A具有野生型和突变菌株，而簇B仅由突变菌株组成（图3b）。检查群集A和B的转录组学和基因组数据产生了五个发现。

• 首先，所有簇B菌株在covS或covR中具有突变，
• 其次，具有covS（68.5％）或covR（37.5％）突变的大多数菌株被分配至簇B（图3c）。
• 第三，分配给簇A的大多数菌株是类似野生型或在除了covRS之外的主要调节基因中具有突变（如下所述）。
• 第四，群集B菌株具有比群集A菌株显着更多的差异表达基因（图4）。
• 第五，关于两个主要转录组簇，没有明显的简单基因组亚片特异性关联（补充图8b）。

CovRS双组分系统负调节15％转录组的表达，包括关键的毒力因子44。与此一致，CovRS的失活增强了毒力53,56。我们比较了由188个预测的野生型菌株组成的442个菌株的转录组，在covRS中具有不同类型突变的132个菌株，以及在其他主要调节基因中具有不同突变的122个菌株。尽管簇B仅含有covRS突变株（图3c），但相当大比例的covS（20.4％）和covR（45.8％）突变株与簇A中的几乎所有野生型样株组合（图3c）。 covRS突变株主要由两个不同的转录组聚类组成的发现表明covRS中的多态性不相同，如前所述[57,58]，并且群集A covRS突变体与野生型菌株的分组表明某些多态性可能比其他多态性具有更少的功能性后果。仅A群（n = 33）和群集B（n = 83）中的covRS突变株的PCA将该分组重现为两个不同的簇（图5a）。接下来，我们使用基于距离的聚类来检验这样的假设：PCA（图5a）不明显的额外亚结构可能存在于转录组数据中（图5b）。调查结果载于补充说明。
为了检验B群covRS菌株与A群covRS菌株相比具有独特转录组的假设，我们检测了两组基因的转录，发现了142个差异表达的基因（补充表12）。这两个簇在差异表达基因的互补性以及差异表达基因的改变的转录物变化（向上或向下）的大小上不同。许多（32％）差异表达的基因具有五倍或更高的改变的转录物水平，包括编码毒力因子的关键CovRS调节的基因，例如SPN和SLO，Mga调节子，HasABC和SpeB（补充表12）。此外，与群集A covRS菌株相比，群集B covRS菌株具有显着增加的移码诱导插入缺失和无义突变的频率（双尾Fisher精确检验，P <0.0001），可能使该调节系统失活（丧失 - 功能突变）。
鉴于编码B群中多个关键毒力因子的基因转录本增加，我们假设B群菌株比A群covRS突变菌株毒性更强。与该假设一致，通过使用小鼠感染模型分析来自每个簇的四种菌株的毒力，显示簇B菌株导致显着更高的死亡率（图5c）和更大的病变，具有更多的组织破坏（图5d）。比较ropB突变可变地影响SpeB半胱氨酸蛋白酶毒力因子的表达和活性的类似研究在补充说明，补充图8d环和补充表13和14中给出。

单核苷酸插入显着改变毒力。

分泌的R28蛋白是一种GAS毒力因子，已被研究作为潜在的候选疫苗60,61。该蛋白质由位于ICE样元件上的Spy1336 / R28基因编码，注释为差异区域2（参考文献）62,63（图6a）。该37.4kb的DNA区段与B组链球菌63的染色体中存在的区域> 99％相同。在我们对最初50个菌株的转录组研究中，我们观察到大约三分之一的菌株表达低水平的Spy1336 / R28转录物，而三分之二的菌株表达高水平的Spy1336 / R28转录物。相邻基因（Spy1337）具有相同的表达模式（图6b）。 Spy1336 / R28和Spy1337转录物水平与遗传结构，地理位置，分离年份或MGE-50基因型之间没有相关性。这一令人费解的发现促使我们使用SEER64,65对所有442株菌株的从头组装进行了全基因组关联研究（GWAS），因为这442株我们有相关的转录组数据。根据Spy1336 / R28和Spy1337基因的转录水平，我们检查了具有低或高转录物表型的菌株是否与任何遗传事件（例如，SNP，插入缺失或重组）显着相关。对于两种表型（高和低转录水平），100％的重要 k-mers映射到Spy1336 / R28和Spy1337之间的基因间区域（补充图10a），这导致鉴定位于Spy1336 / R28之间的基因间区域的poly(T)同聚物区域中的变体和Spy1337基因（图6c）。与高转录物表型正相关的显着k-mers具有10个T残基，与高转录本表型负相关的那些具有9个T残基。还通过50-和442-菌株数据集的表达数量性状基因座（eQTL）分析66,67（补充图10b）鉴定了10T变体与Spy1336 / R28和Spy1337的增加的转录物水平的关联。
与亲本菌株（9个T残基）相比，同基因突变菌株（10个T残基）具有显着增加的Spy1336 / R28和Spy1337的转录水平（图6d）;如在小鼠坏死性肌炎感染模型中评估的那样，引起显着更大的肉眼和微观病变，以及更接近的死亡率（图6e，f）;对人体多形核白细胞的离体杀伤具有显着的抗性（P <0.05;图6g）;并产生更多分泌的和细胞相关的Spy1336 / R28蛋白（图6h）。通过对生长至中指数或早期稳定期的同基因菌株的RNA-seq分析进一步证实由该同聚物核苷酸区中T残基数目的变化引起的Spy1336 / R28和Spy1337的改变水平（补充表15）。因此，在该同聚物区中插入或缺失单个T残基显着改变了Spy1336 / R28和Spy1337的转录水平，转录组和菌株毒力。

SC2Asubclade菌株在小鼠中毒性更强。

全基因组序列和转录组数据显示emm28菌株之间存在显着差异，我们推断这些基因组和转录组学变化可能导致毒力的显着变化。为了验证这一假设，我们在坏死性肌炎小鼠模型中评估了相对于SC1A参考菌株MGAS28426的50个emm28菌株（上述初始转录组研究中使用的相同系统发育不同菌株组;补充表7）的毒力16,39， 68。事实上，对于所有已知的主要调节基因，所有这些菌株（96％）都是野生型。作为群体，SC1A和SC1B菌株的毒力没有显着差异（图7a，b）。与之形成鲜明对比的是，SC2A菌株的毒力显着高于SC1A和SC1B菌株（图7a，b）。我们假设SC2A菌株的毒力增加可能至少部分是由于nga，ifs和slo基因的表达显着增加，从而导致SC2A生物分泌的SPN和SLO毒素的产生增加，如图所示。其他GAS血清型1,6,15,16。与该假设一致，SC2A菌株具有比SC1A或SC1B菌株显着更高的nga转录物水平和SPN活性（图7c）。对于ifs（P <0.001，Mann璚hitney U-test）和slo（P <0.001，Mann璚hitney U-test） - 同一操纵子中的另外两个基因也观察到相同的情况。
为了明确证明SC2A菌株的毒力显着高于SC1A染色，部分原因是由于更大的nga璱fs璼lo启动子活性，我们用SC1A nga启动子替换了SC2A参考菌株MGAS27961的nga启动子。具有SC1A启动子的同基因突变体菌株在体外产生显着较低的SPN活性（图7d），并且在小鼠坏死性肌炎感染模型中引起显着较低的死亡率和组织破坏（图7e，f）。

讨论

我们使用GAS作为模型病原体来研究emm28化脓性链球菌菌株中群体基因组学，转录组学和毒力之间的复杂相互作用。我们发现转录组变化的幅度（差异表达基因的数量）与总体基因组 - 遗传距离之间没有简单的相关性（补充图9和补充说明）。总的来说，我们的分析表明，涉及基于种群的菌株样本的多种类型的高维数据的整体方法可以提供对不太容易通过较少综合和综合方法发现的病原体环境相互作用的新理解。通过利用转录组特征分析，我们发现常用生物信息学方法会遗漏~5％的在的全局调节因子中具有突变菌株。这一发现可作为调查基因组转录关系研究的警示。
我们开发了基于群体的多维基因组和转录组数据集，以及包括GWAS和eQTL在内的统计方法，以确定负责改变Spy1336 / R28转录水平的分子事件 - 一种编码R28蛋白毒力因子和候选疫苗的基因60,63 69。这些发现增加了细菌中一个新兴的主题，即基因间区域看似适度的变化可以改变基因表达并具有适应性15,70？4。关于如何增加Spy1336 / R28和Spy1337的转录水平可能增强毒力存在几种可能性。一种可能性是改变的毒力表型单独或主要由Spy1336 / R28的产生增加引起，Spy1336 / R28是已知的毒力因子60,75。确切地说，这种蛋白质如何有助于毒力尚不清楚，尽管据报道它可以促进人类上皮细胞的粘附[60,76]。第二种可能性是Spy1337调节蛋白直接或间接地改变其自身和可能影响毒力的其他基因的转录。为了验证这一假设，我们对含有9或10个Ts的同基因菌株进行了RNA-seq分析，发现在中指数期，只有两个基因（Spy1336 / R28和Spy1337）显着上调，而在早期 - 稳定期。Spy1336 / R28和Spy1337以及33个其他基因被上调，165个被下调（补充表15）。与9T亲本生物相比，10T同基因突变株中的三基因fruRBA操纵子（Spy0641，Spy0642和Spy0643编码蛋白参与果糖利用）高度上调。 fruR或fruB基因的失活显着降低了人全血中或多形核白细胞存在下突变菌株的存活率77。与该发现一致，10T同基因突变株在坏死性肌炎小鼠模型中具有显着更高的毒力（图6e，f），对人多形核中性粒细胞的杀伤具有显着增强的抗性（图6g），并产生更多的分泌物。和细胞相关的Spy1336 / R28蛋白质（图6h;图6i中的模型）。
最简单的假设解释了这个基因间区域中一个T残基的插入或缺失如何改变Spy1336 / R28和Spy1337的转录水平是由Spy1337编码的转录调节因子直接与该基因间区域结合并同时增加两个基因的转录。在该假设下，同聚物核苷酸区可以是：（i）是整个Spy1337共有结合位点的一部分或构成，或（ii）位于两个共有结合位点侧翼的间隔区中。在第一种情况下，具有10 Ts（与9 Ts相比）的同聚物区域将构成更好的共有结合位点，而在第二种情况下，间隔区域中存在10 Ts（与9 Ts相比）将位于侧翼推定位置在DNA螺旋中更有利的空间取向的共有结合位点，用于结合Spy1337转录调节因子（图6c）。正在进行其他研究来解决这个问题。
本研究以新颖的方式使用机器学习和eQTL分析。我们的转录组数据集促进了机器学习的使用，以分析和正确分类基于单独的基因组序列数据分析错误识别的调节突变株。直到最近，大多数使用机器学习的病原微生物的研究都使用了基于DNA序列的数据，通常侧重于预测对抗菌药物的抗性4,78？3。 eQTL分析已用于人和其他真核生物中的表达数据84？0;该研究将eQTL分析应用于细菌数据集，通过生成的广泛转录组数据实现。我们发现，基因间区域的插入缺失与五个基因的表达改变显着相关 - 两个在顺式（Spy1336 / R28和Spy1337）和三个在反式（Spy1338，Spy1339和Spy1340） - 通过使用50的转录数据在指数阶段的应变。类似地，通过使用来自早期稳定期的442个菌株的转录物数据，发现顺式（Spy1336 / R28和Spy1337）和与47个额外基因的反式关联（错误发现率<0.0005）（补充图10b和补充表16）。重要的是，通过eQTL分析鉴定的49个基因中的60％也通过等基因（10T）突变体和（9T）亲本菌株的RNA-seq分析差异表达。
转录组数据还使我们得出结论，尽管HGT事件可以并确实改变了转录组;大多数转录组的变化是由SNP（错义或无义突变）和影响主要调节基因如covR / covS，ropB和mga的短插入引起的，并导致同源编码蛋白的截短53,5,57。这些发现与观察结果一致，即这些调节基因是GAS菌株群体基因组分析中具有最高密度多态性的基因之一[6,91,92]。
由于未知原因，emm28 GAS菌株在产褥期败血症（儿童发烧），女性生殖道感染和新生儿感染17,31？4中的比例过高。虽然我们的研究并非旨在解决细菌种群结构与感染菌株的详细临床表型之间非常复杂的关系，但一项观察值得评论。对于美国患者的951株分离株，可获得合理详细的感染类型信息。我们发现，与来自美国的非SC2A菌株相比，来自美国的SC2A菌株的显着更高比例与产褥败血症，新生儿感染和女性生殖道感染相关（2（1）= 5.854; P = 0.015;补充表1）。在这方面，我们注意到作为一组，SC2A菌株在小鼠坏死性肌炎实验中也显着更具毒性（图7a，b）。
总之，我们的研究可以作为一个范例，说明如何有效整合基于种群的样本生成的多维数据集，以产生关于微生物遗传和病原体环境相互作用的新知识。与三种类型数据中的任何一种或两种数据的研究相比，三种不同类型数据的整合导致对病原体的分子遗传学的更多理解。该策略通常适用于任何微生物 - 致病性或其他 - 并可能导致新的治疗方法。

方法

全基因组测序和多态性分析。

如前所述，使用Illumina NextSeq 550仪器进行菌株生长，染色体DNA的分离，成对末端文库的产生和多重测序，如前所述[1,6,93]。用于细菌基因组分析的管道显示在补充图2a中。通过使用Musket94校正修剪的序列读数，并通过使用SMALT将其定位至参考血清型M28菌株MGAS6180（GenBank登录号CP000056）63的基因组（参见URL）。通过使用描述的FreeBayes（参见URL）和用于检测插入缺失的Pilon95鉴定SNP和插入和缺失（插入缺失）。 SNPfx.pl（内部开发的PERL脚本;参见URL）用于确定SNP的性质（编码与非编码，同义与非同义，等等）。或者，结合使用，SPAdes算法用于从头基因组装配96。短读序列分型2（SRST2）用于鉴定基因，等位基因和多位点序列类型97。算法Gubbins用于检测HGT事件98。 hierBAPS99（群体结构的分层贝叶斯分析;参见URL）用于确定种群结构，SplitsTree用于估计系统发育树和网络100。通过使用范围在50和200之间的簇的数量的先前上限值来运行hierBAPS，其具有估计算法的五个重复，每个运行收敛到群体结构的相同后验模式估计。
通过使用具有R9.5流动池和快速条形码试剂盒（SQK-RBK004）的Oxford Nanopore MinION仪器进行24个菌株的长读序列（补充表19）。这些菌株选自50个菌株，其中RNA-seq表达分析一式三份进行，以数字代表主要的遗传进化枝，包括多种原噬菌体和ICE含量（MGE基因型），并包括在R28中不同的菌株。启动子区T-核苷酸同聚物束。如先前所述122，使用Unicycler101生成纳米孔长读数和Illumina短读数据的混合组件。使用Trimmomatic进行读取质量过滤和修剪的质量指标，使用Musket完成读取错误纠正，使用SPAdes / Unicycler生成从头组装，使用SRST2进行多位点序列类型和emm测定，使用SMALT进行读取映射，并且用FreeBayes完成多态性发现（补充表20）。
基于连锁序列核心染色体SNP通过邻近连接推断菌株之间的系统发育，并且通过使用hierBAPS基于整个核心基因组序列的贝叶斯分析来定义相关菌株的进化枝。通过使用k-means（在下面和在补充说明中描述）选择用于转录组分析的参考菌株。

用于转录组分析的菌株选择。

为了从2,101个emm28菌株的测序群体中选择一个分离子集，这些菌株对于给定大小的子集大致跨越尽可能多的遗传变异，我们使用基于投影的方法，类似于细菌中使用的群体结构校正GWAS方法SEER103。首先，通过使用核心SNP，分别计算属于给定谱系（SC1A，SC1B，SC2A或SC2B）的所有emm28分离物的基于SNP的成对距离矩阵。我们排除了由Gubbins98鉴定的MGE和潜在的重组区域。随后，我们按照每个谱系代表的总种群大小的比例对分离物进行比例取样。然后使用每个谱系的距离矩阵将所有谱系分离物投射到具有多维缩放的三维欧几里德空间中，如SEER103中所示。对于从谱系中选择的子集的给定总大小k，使用具有200个随机重启104的k均值算法来识别k个质心的最佳集合以跨越存在于三维多维缩放投影中的变化。谱系中存在遗传变异。选择具有到每个质心的最小欧几里德距离的分离物以确定k个代表性分离物的最终子集。

RNA-seq文库制备和测序。

emm28菌株一式三份生长并在两个时间点收获。代表四种主要遗传子片段的50个菌株（图1a，b和补​​充表7）在生长的中指数和早期 - 静止期一式三份进行测定。根据制造商的说明用RNeasy试剂盒（Qiagen）提取RNA。核糖体RNA（rRNA）用革兰氏阳性细菌（Illumina）的Ribo-Zero rRNA去除试剂盒耗尽，如前所述[1105]。分别用RNA Nano和Pico芯片（Agilent Technologies）和Agilent 2100生物分析仪评估总RNA和rRNA耗尽的RNA的质量。用来自ScriptSeq Index PCR Primers试剂盒（Illumina）的索引反向引物制备互补DNA（cDNA）文库，并用AMPure XP珠（Beckman Coulter）纯化。用高灵敏度DNA芯片（Agilent Technologies）评估cDNA文库的质量。对于每个样品，用Qubit dsDNA HS测定试剂盒（Invitrogen）荧光测定法测量cDNA文库浓度。将cDNA文库稀释，合并，并用Illumina NextSeq仪器测序。该相同方案用于9T和10T同基因菌株的比较RNA-seq分析（补充表15）。
emm28菌株作为单一培养物生长并在一个时间点收获。 461单体的转录组分析（补充表1和10）如描述52中的RNAtag-seq进行，具有本文所述的修饰。通过使用Qubit RNA BR测定试剂盒（Life Technologies）荧光测定从每个菌株分离的总RNA。将来自每个样品的约400ng在94℃下以16μl在1μlFastAP缓冲液中片段化3分钟，用FastAP碱性磷酸酶（Thermo Fisher Scientific）去磷酸化12分钟，37℃，最终体积为20μl，并在5℃下磷酸化。以T4多核苷酸激酶（New England Biolabs）在37℃终止30分钟，最终体积为82μl。根据制造商的说明书，在1.5ml微量离心管和DynaMag-2磁体（Invitrogen）中，使用体积（164μl）的Agencourt RNAClean XP顺磁珠纯化片段化的去磷酸化的总RNA。最终的洗脱体积为12μl。通过条形码的连接实现RNAtag-seq过程期间总RNA的汇集，使得来自相同菌株的所有RNA片段用单独的条形码明显标记。我们使用了早期研究中描述的16种独特的条形码寡核糖核苷酸52（补充表17）。对于每个菌株，通过使用体积为20.11的T4 RNA连接酶1（ssRNA连接酶; New England Biolabs）将5μl片段化的磷酸化总RNA连接到1μl各自的寡核糖核苷酸上，终浓度为5μM。反应在22℃下进行90分钟。结扎后，通过加入59.9l RLT缓冲液（RNeasy Mini试剂盒; Qiagen）将每个样品的体积增加至80μl，并与含有80μlRNA结合缓冲液（RNA Clean＆Concentrator-5; Zymo）的1：1混合物混合。研究）和在150ml微量离心管中的80μl100％乙醇。因此，通过将对应于构成一个特定池的8个菌株的总RNA依次通过一个Zymo柱并将它们浓缩在一起，为每组48个样品制备6个样品，每个样品包含8个样品，如补充图7a所示。每个池的最终洗脱体积为32μl。
用RNA Pico芯片（Agilent Technologies）评估总RNA库的质量。我们制备了57个池，其中含有来自8个菌株的总RNA和来自5个菌株的总RNA的一个额外的池，总共461个菌株。 Ribo-Zero rRNA去除试剂盒（革兰氏阳性细菌）用于从池中消除不需要的rRNA。通过使用RNA Pico芯片（Agilent Technologies）分析ribulplepleted RNA的质量。如前所述进行第一链cDNA合成52。对于每个库，将12μl核糖蛋白RNA与2μlAR2寡核苷酸（其与本研究中使用的所有16个条形码寡核糖核苷酸中存在的区域互补（补充表17））混合，并在70℃下变性2分钟。然后，使用AffinityScript逆转录酶（Agilent）以20μl的体积在55℃下进行第一次cDNA合成55分钟。随后将RNA在0.09N NaOH中在70℃下降解12分钟，并用终浓度为76.9mM的乙酸中和，最终体积为26μl。加入14μl水后，用2.5μl体积（100μl）的Agencourt RNAClean XP顺磁珠纯化单链cDNA（sscDNA），将sscDNA与珠子一起重悬于5μl水中。在珠子中，将sscDNA与2μl3Tr3衔接子52混合，并使用T4 RNA连接酶1以20μl的体积连接。将反应物在22℃温育过夜，然后使用2.5倍体积的Agencourt RNAClean XP珠子进行两次连续的净化反应，并用25μl水洗脱。
用通用引物univP5（参考文献52）和四种不同的P7条形码衔接子之一（补充表17）进行文库扩增。在PCR富集试验后，sscDNA库以四个组的形式组织并在50μl的最终体积中单独扩增，以确定正确的扩增条件;然后使用Agencourt RNAClean XP珠子进行两次净化步骤，并在20μl低TE缓冲液（10mM Tris和0.1mM EDTA）中洗脱。对于每个样品，通过使用HighSensitivity DNA芯片（Agilent Technologies）测定cDNA文库的平均大小，并用Qubit dsDNA HS测定试剂盒（Invitrogen）以荧光测定法测量cDNA文库浓度。
在方案的这一点上将样品再合并一次。将对应于四个库的cDNA文库以等摩尔量混合在一起，每个库对应于八个菌株。该过程再重复14次;因此，我们最终得到了15个超级池，共计58个池，代表461个菌株（补充图7a）。对应于每个超级池的文库分别掺入10％PhiX文库以改善簇多样性并用Illumina NextSeq仪器测序。

RNA-seq数据分析。

用于处理RNA-seq数据的生物信息学管道在补充图2b中给出。
分析50个emm28菌株，一式三份生长并在两个时间点收获。对于每次测序运行，Illumina bcl2fastq软件（参见URL）用于将Illumina生成的BCL碱基调用文件转换为FASTQ文件。通过使用FASTQC软件评估测序数据的读取质量（参见URL）。使用Trimmomatic106进行适配器污染和读取质量过滤。通过使用EDGE-pro107将读数定位于参考菌株MGAS6180的基因组，并且从随后的分析中排除映射至rRNA和转移RNA基因的读数。此外，根据补充说明中描述的策略，将具有低表达的基因排除在下游分析之外。通过使用DESeq2进行差异表达分析（参考文献108）。该相同的管道用于在Spy1336 / R28和Spy1337之间的同聚物区中含有9或10个Ts的同基因菌株的比较RNA-seq分析。
在一个时间点分析442个emm28单一菌株。使用bcl2fastq将从超级池的读取解复用到单独的池中。分别使用FASTQC和Trimmomatic完成了读取质量评估和读取质量修剪。通过使用FASTX-Toolkit（参见URL），根据内联条形码将每个池中的读取解复用为对应于各个样本的单独fastq文件。每个池的读数中位数为92.6百万（补充图7b），每个样品每个样品的读数中位数在800万到2400万之间（补充图7c）。补充说明中描述了低表达基因的排除。没有发现显着的批次效应与表达数据相关联。通过使用R包variancePartition109估计，发现批次效应解释的方差百分比<2％。
通过使用在NOISeq包110中实施的NOISeq-sim进行最终442个单一菌株组的差异表达分析。与参考菌株MGAS28737相比，在441个菌株的每一个中鉴定差异表达的基因，参考菌株MGAS28737如补充说明中所述选择。 DESeq2（参考文献108）用于差异表达分析，用于进行两组比较的情况，其中每组包含一种以上的菌株。

机器学习与随机森林分析。

使用函数TreeBagger（MATLAB R2016b，Statistics and Machine Learning Toolbox，MathWorks），使用MATLAB进行随机森林分析。目的是根据菌株的转录组谱训练随机森林59，将异常菌株分类为四类（covR，covS，ropB和野生型）。对随机森林进行了训练，其中转录组数据针对仅在单一主要全球调控基因（covR，covS或ropB）中突变的菌株生成，并且已知所有已知主要调控基因的野生型菌株（n = 283）。因此，训练数据由283个菌株组成，其中已知1,614个基因和类别标签上的转录组谱。测试数据由8个异常菌株组成，其中已知1,614个基因的转录组谱，但类别标记未知。在学习最终模型之前，应用以下特征选择过程。建立了所有基因的初始随机森林（1,000棵树），并通过使用内置的特征重要性测量来估计基因的预测重要性。该过程重复十次，最终的特征重要性值作为各次运行的平均值。从最重要的功能开始，我们根据给定的重要顺序依次包含更多功能。对于每个特征子集，我们进行了双重交叉验证，通过使用三分之二的训练数据构建随机森林（100棵树），并对三分之一的训练数据进行评估。通过100次交叉验证迭代的平均样本外分类准确度来测量子模型的样本外性能。随着更多特征的增加，分类准确度的提高迅速趋于稳定，并且在此基础上，为最终模型选择了十个最具信息性的基因。最终模型用于预测8个异常菌株的类概率（补充表11）。

GWAS分析。

使用SEER103的GWAS分析使用442个菌株的从头组装进行，其中我们还具有RNA-seq数据。 GWAS用于鉴定与二元表型分组显着相关的遗传变体（高转录物表达= 1;低转录物表达= 0），根据Spy1336 / R28基因的转录水平定义。绘制442个菌株的Spy1336 / R28基因的标准化转录物水平（计数），得到两个视觉上非常不同的组 - 低（三分之一的菌株）和高表达者（三分之二的菌株） - 类似于我们对50株菌的观察结果（图6b）。我们确定了阈值，发现归一化计数小于261.5的菌株被认为是低表达者（编码为0），等于或大于261.5计数的菌株被认为是高表达者（编码为1）。二元表型文件和de novo组装的fasta文件提供给SEER，k-mers通过使用fsm-lite从组装的读数计算（参见URL）。为了考虑种群结构，使用了由Mash111计算的距离矩阵。运行SEER分别产生17和13个显着的k聚体（调整后的P值<10-8），其与高或低转录物表达呈正相关或负相关。

eQTL分析。

使用R包MatrixEQTL112进行eQTL分析。对于eQTL分析，通过使用前十个主成分作为协变量，在模型中考虑了种群结构。如果多态性在所考虑的基因的1kb内，则认为关联是顺式的。
通过使用坏死性肌炎的小鼠模型对50种天然存在的和同基因型血清型M28 GAS突变株进行毒力研究。如前所述进行用血清型M28 GAS菌株进行的小鼠坏死性肌炎研究39,68。使用在初始RNA-seq实验中使用的50种天然存在的菌株，包括毒力参考菌株MGAS28426。菌株MGAS28426用作毒力参考菌株，因为它在基因组上代表SC1A菌株，并且它用作分泌的SPN测定的野生型参考。制备每种菌株的冷冻原种，并通过在连续稀释后制备从解冻的培养物中回收的菌落形成单位（c.f.u.）进行定量。将免疫活性的4周龄雌性远交CD1小鼠（Envigo）随机分配到菌株处理组，并在5'08c.f.u的右下肢接种。每种细菌菌株（n = 20mice perstrain;总共1,000只小鼠）。该剂量是在两个试验实验的基础上选择的，该实验表明毒力参考菌株MGAS28426在接种剂量为5？08c.f.u时引起约50％的近死亡率。该策略有助于鉴定具有显着增加或降低的毒力的比较菌株。通过使用具有以下变量的功率计算来选择小鼠样本大小：= 0.05;功率（1？）= 0.8;组间存活率差异= 0.4;组大小比= 1。
对于ccovRS突变体菌株比较，使用4个簇A菌株和4个簇B菌株（n = 45mice perstrain），剂量为5？08c.f.u。对于亲本野生型（27961）和同基因启动子突变体（27961-SC1A nga启动子）比较（n = 40miceper菌株），5？08c.f.u。被使用了。对于ropB突变体菌株比较，以1？09c.f.u的剂量使用3组I组和4组II组（n = 40mice perstrain）。对于9T和10T同基因菌株（27961-9T和27961-10T）和对照菌株28085-10T（n = 20miceper菌株），我们使用5？08c.f.u的剂量。获得了从分配给组织病理学分析的小鼠取得的肢体的代表性显微图像和显微图像。用于产生同基因突变体的寡核苷酸显示在补充表18（补充说明）中。
所有动物研究均按照由Methodist Hospital Research Institute的机构护理和使用委员会审查和批准的方案（AUP0615？041）进行。每天至少监测一次小鼠，并根据国家研究委员会（美国）2011年实验动物护理和使用指南更新委员会和护理指南提供的国际公认标准指南确定接近死亡率。和实验动物的使用113。存活数据表示为Kaplan璏曲线，并且用对数秩检验（Prism6; GraphPad）确定统计学显着性差异。

统计分析。

除非另有说明，误差条代表s.d.，Pvalues用Fisher精确，Mann璚hitney U或对数秩检验计算。 MatrixEQTL程序包报告使用了错误发现率。贝叶斯聚类用于定义emm28群体中的进化枝和亚进化枝。 k-means和基于距离的聚类用于识别由PCA产生的应变的二维空间聚类中的质心，并分别在聚类中找到另外的子结构。随机森林分析是在MATLAB中使用函数TreeBagger完成的。 R包variancePartition用于确认没有显着的批次效应。使用确定系数（R2）统计来研究遗传距离和转录组重构程度之间是否存在相关性。使用比例的onetailed测试确定II组RopB菌株含有显着比例的影响功能结构域的突变。 Pearson相关系数（r）用于比较来自一式三份分析的菌株的RNA-seq数据与从作为单体生长的菌株收集的RNAtag-seq数据。

入藏代码。

本文报道的序列已经保藏在国家生物技术信息中心（NCBI）序列阅读档案（见URL）中，其中BioProject登录号为PRJNA434389，NCBI基因表达综合登记号为GSE113058。

伦理。

所有小鼠研究均按照休斯顿卫理公会研究所的Institutional Animal Care and Use Committee批准的方案（AUP0318？016）进行。所有具有人血液和血液成分的研究均按照国家过敏和传染病研究所的人类受试者机构审查委员会批准的方案（01-I-N055）进行。所有研究志愿者都签署了书面知情同意书报告摘要。