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, .
【非常牛的案例文章！隆重推薦！】

群體基因組學(xué)，轉(zhuǎn)錄組學(xué)和毒力的綜合分析為化膿性鏈球菌發(fā)病機制提供了新的見解

摘要

化膿性鏈球菌每年在全世界造成7億人感染，盡管經(jīng)過一個世紀的努力，但沒有針對這種細菌的許可疫苗。盡管已經(jīng)發(fā)表了許多關(guān)于細菌病原體的大規(guī)?；蚪M研究，但大型細菌群體中基因組，轉(zhuǎn)錄組和毒力之間的關(guān)系仍然知之甚少。我們對2,101個emm28化膿性鏈球菌侵染菌株的基因組進行了測序，我們從中選擇了492個系統(tǒng)發(fā)育不同的菌株用于轉(zhuǎn)錄組分析，50個菌株用于毒力評估。數(shù)據(jù)整合提供了對該模式生物的毒力機制的新穎理解。全基因組關(guān)聯(lián)研究、表達數(shù)量性狀基因座分析、機器學(xué)習(xí)和同基因突變菌株鑒定并證實了基因間區(qū)域中的單核苷酸插入缺失，其顯著改變?nèi)洲D(zhuǎn)錄物譜和最終毒力。我們使用的綜合策略通常適用于任何微生物，并可能導(dǎo)致許多人類病原體的新療法。

無論生態(tài)位還是宿主范圍，所有細菌種類都包含遺傳多樣性菌株。人們對微生物分子遺傳學(xué)知之甚少的領(lǐng)域是大量傳染性菌株的基因組，轉(zhuǎn)錄組和毒力之間復(fù)雜的相互作用。遺傳變異可能會影響基因轉(zhuǎn)錄水平，但這種情況的真實程度及其對發(fā)病機制的影響仍不明確【good！】。盡管大型基因組學(xué)研究已經(jīng)發(fā)表1，但在比較轉(zhuǎn)錄組領(lǐng)域的研究卻少得多7和涉及自然種群的毒力分析10,11。此外，除了一項涉及相對較小的大腸桿菌菌株樣本12的研究外，尚未研究基因組，轉(zhuǎn)錄組和毒力之間的關(guān)系。對基于種群的菌株樣本的各種數(shù)據(jù)進行綜合分析可能會對疫苗制劑，新療法和診斷以及公共衛(wèi)生等領(lǐng)域的轉(zhuǎn)化研究工作產(chǎn)生影響。
化膿性鏈球菌（A群鏈球菌（GAS））是一種嚴格的人類病原體，每年在全球兒童和成人中引起超過7億的感染13。人類感染的嚴重程度從相對輕微的咽炎（“鏈球菌性咽喉炎”）到極端嚴重和危及生命的感染，如敗血癥和壞死性筋膜炎/肌炎，通常被稱為“肉食”疾病。該生物體還會引起皮膚和軟組織感染，并導(dǎo)致感染后遺癥，如風濕熱和風濕性心臟病，這是全球發(fā)病率的重要原因[13,14]。
GAS已被用作模型生物，用于研究菌株類型和疾病表型之間的關(guān)系，以及流行病學(xué)[1,6,15-7]。 emm28菌株是美國和許多歐洲國家中與侵襲性GAS感染相關(guān)的前五種emm類型之一18？3。由于仍然無法解釋的原因，一些emm類型或M蛋白血清型的菌株與特定類型的人類感染非隨機相關(guān)17,24？0。例如，emm28 GAS菌株在產(chǎn)褥期敗血癥（兒童發(fā)燒）和新生兒感染17,31？4的病例中反復(fù)出現(xiàn)。
盡管在所選生物的基因組學(xué)方面取得了重大進展，但對于共享最近共同祖先的細菌菌株的克隆相關(guān)后代的轉(zhuǎn)錄組多樣性的性質(zhì)和程度知之甚少。關(guān)于這個問題的數(shù)據(jù)對于增強對自然群體中細菌進化，表型多樣化和微生物流行的理解至關(guān)重要。為了解決這些知識gap，我們對基于綜合種群的研究中回收的2,101種emm28 GAS菌株的基因組進行了測序，并使用所得的系統(tǒng)發(fā)育信息選擇代表性菌株來分析轉(zhuǎn)錄組（n = 492株）和毒力（n = 50株）。數(shù)據(jù)整合提供了對該模型生物的生物學(xué)的新的理解，包括在相對密切近緣的生物進化枝中具有驚人的轉(zhuǎn)錄組變異幅度。統(tǒng)計學(xué)方法和機器學(xué)習(xí)的應(yīng)用促進了一種新的分子遺傳過程的發(fā)現(xiàn)，該過程支持一些GAS菌株中增強的毒力。

結(jié)果

群體結(jié)構(gòu)和時間分布。

我們對在北美和歐洲的六個國家中分離的2,101個emm28 GAS菌株的基因組進行了測序，這些菌株在26年期間：1991年至2016年（表1，補充表1和補充圖1）。所有菌株都作為綜合群體研究的一部分被回收。將基因組測序至202倍平均覆蓋度（補充圖2a，補充表1和方法）35。遺傳關(guān)系的推斷是通過使用核心基因組中存在的SNP（即沒有移動遺傳元件（MGE）的基因組，例如前噬菌體和整合 - 結(jié)合元件（ICE））來進行的（補充表1）。主要的emm28 GAS群體通過貝葉斯聚類分布到兩個主要進化枝（進化枝1和2）和四個進化枝（指定為SC1A，SC1B，SC2A和SC2B）中（圖1a，b）。進化枝2生物與進化枝1的部分區(qū)別在于28.0 kb水平基因轉(zhuǎn)移（HGT）區(qū)塊，其貢獻520個核心SNP，以及位于該HGT區(qū)塊外的19個核心SNP（補充表2）。這個28.0-kb的HGT區(qū)塊包括編碼分泌的毒素NAD + - 糖水解酶（SPN）和鏈球菌溶血素O（SLO）——已知GAS毒力的主要貢獻者6,16,36？9的nga–ifs–slo操縱子 。 ifs編碼SPN40的內(nèi)源性抑制劑。重要的是，重組獲得nga–ifs–slo操縱子的高表達變體可以增加在靈長類動物上呼吸道的存活率，增強毒力，并引發(fā)洲際流行病[1,6,16,41]。進化枝2生物具有與糞鏈球菌亞種equisimilis42中存在的類似三基因操縱子具有99％序列同一性的nga–ifs–slo區(qū)域，并且可能通過重組事件被亞類2A 2A GAS獲得。
SC1B包含的菌株最多，占菌株的49.7％，其次是SC2A（26.8％），SC1A（22.3％）和SC2B（0.53％）。屬于子SC1A，SC1B和SC2A的菌株根據(jù)年份，地理位置和MGE含量而變化（圖1c和補充圖3）。 SC1A菌株的顯著時間位移同時發(fā)生，美國SC2A菌株激增，其中約55％的菌株是SC2A。 SC1B菌株主要分離自芬蘭，挪威，法羅群島和冰島的患者，而SC1A是加拿大流行的亞類。評估群組中的MGE多樣性（補充說明和補充表3？），并且20種最豐富的MGE-50基因型占菌株的90％（圖1c）。
接下來，我們使用綜合策略從種群角度研究動物感染模型中基因組變異，轉(zhuǎn)錄組變化和毒力差異之間復(fù)雜的相互作用。

轉(zhuǎn)錄組簽名和群體結(jié)構(gòu)。

為了確定基因表達的不同模式是否與emm28群體結(jié)構(gòu)非隨機相關(guān)，我們首先對50個菌株的子集進行轉(zhuǎn)錄組RNA測序（RNA-seq）分析，這些菌株在遺傳上代表三個數(shù)值顯著的亞類（即SC1A，SC1B和SC2A）（補充表7）。根據(jù)方法中描述的標準，從2,095個M28菌株的主要樣品中選擇菌株進行RNA-seq分析。 50株來自不同國家，國家和地區(qū)，具有不同的MGE含量。 RNAseq分析在中指數(shù)和早期 - 靜止生長期一式三份（三個生物學(xué)重復(fù)）進行（方法，補充圖2b和補充表7）。
盡管50種菌株的基因組背景不同，但分析的菌株數(shù)量加上一式三份樣品中轉(zhuǎn)錄水平之間極高的相關(guān)系數(shù)（補充圖4），允許鑒定與種群結(jié)構(gòu)有關(guān)的不同轉(zhuǎn)錄組改變（圖1）。我們發(fā)現(xiàn)了兩個出乎意料地具有“異常”轉(zhuǎn)錄組的菌株（圖2a，b）。手動檢查這兩種異常菌株的基因組序列數(shù)據(jù)，在covS全局調(diào)控基因中發(fā)現(xiàn)了兩個獨立的大缺失事件（圖2a，b）。編碼全局轉(zhuǎn)錄調(diào)節(jié)因子如CovR / CovS，RopB和Mga的基因的突變可以改變轉(zhuǎn)錄組的顯著比例（5~5％）43。
對于在所有已知主要調(diào)節(jié)基因中具有野生型等位基因的46個菌株（圖2c，d），差異表達基因的最大數(shù)目發(fā)生在分配給不同亞類的菌株之間（補充圖5）。在midexponential階段，當我們比較SC2A菌株與SC1A和SC1B菌株的轉(zhuǎn)錄組（分別為32（2％）和15（0.9％））時，觀察到最大數(shù)量的差異表達基因（補充圖5a）和補充表8）。當比較三個遺傳子片段的早期靜止轉(zhuǎn)錄組時，類似的模式是明顯的，但差異表達基因的數(shù)量相當大（5？倍）（補充圖5a和補充表8）。與SC1A和SC1B菌株相比，SC2A菌株具有最大數(shù)量的差異表達基因（分別為318（19.9％）和83（5.2％））（補充圖5a和補充表8）。
大部分差異表達的基因位于水平轉(zhuǎn)移的28.0-kb區(qū)域（HGT區(qū)域）并包括nga–ifs–slo操縱子（補充表8）。在midexponential階段，該HGT區(qū)域內(nèi)基因（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菌株中三個最強烈上調(diào)的基因是nga，ifs和slo，在指數(shù)期間轉(zhuǎn)錄水平增加約4倍，在早期穩(wěn)定期增加約8倍（補充表8和補充圖5b）。
近年來，SC1B菌株引起的感染在包括美國，芬蘭，冰島和挪威在內(nèi)的一些國家有所增加，而SC1A菌株已大幅減少（補充圖3a），從而提高了SC1B菌株可能進化為比SC1A菌株更合適。因此，我們檢查了核心染色體上的遺傳差異（補充表9），區(qū)分了所有SC1A和SC1B菌株，發(fā)現(xiàn)所有SC1B菌株在RivR中都有兩個連續(xù)的非同義突變——grab基因的一個負調(diào)節(jié)因子（蛋白-G相關(guān)的α2-巨球蛋白結(jié)合）。 在中期指數(shù)期，與SC1A和SC2A菌株相比，grab是SC1B菌株中唯一上調(diào)的基因（補充表8和補充圖5c）。類似地，與SC1A相比，在SC1B和SC2A菌株的早期 - 穩(wěn)定期觀察到更高的grab轉(zhuǎn)錄物豐度（補充圖5c）。 GRAB是一種結(jié)合α2-巨球蛋白的細胞表面錨定蛋白，允許它在GAS表面保留蛋白水解活性形式的半胱氨酸蛋白酶SpeB，用于保護GAS免受抗菌肽LL-37的殺傷（參考文獻48,49），并有助于小鼠模型中的侵襲性感染50。我們不能排除導(dǎo)致subclade位移的隨機過程。

驗證單體（RNAtag-seq）分析【無重復(fù)】。

傳統(tǒng)上，通過使用生長至兩個不同生長期的菌株的一式三份生物復(fù)制品進行細菌的轉(zhuǎn)錄組分析。然而，這種方法目前在研究數(shù)百種菌株方面經(jīng)濟上不可行。鑒于RNA-seq的準確性，靈敏度和重現(xiàn)性得到改善，我們假設(shè)使用單一菌株進行大規(guī)模轉(zhuǎn)錄組分析（即缺乏重復(fù)），與病原體的最佳測序深度（5？0百萬測序讀數(shù)）一致基因組大小約為2Mb（參考文獻51），將大大增強對一組相對密切相關(guān)的菌株的轉(zhuǎn)錄組景觀的理解。為了檢驗我們的假設(shè)，我們使用RNAtag-seq52來增加轉(zhuǎn)錄組分析的菌株通量。發(fā)現(xiàn)來自沒有（單體）和具有生物學(xué)重復(fù)的菌株的表達數(shù)據(jù)是高度相關(guān)的（補充圖6和補充說明）。因此，我們進行了通過k均值聚類統(tǒng)計策略（方法）選擇的442個遺傳代表性和多樣化單一菌株（補充表10和補充圖7）的群體轉(zhuǎn)錄組分析。

不同菌株的群體轉(zhuǎn)錄組分析。

為了檢驗442個單一菌株的轉(zhuǎn)錄組之間的關(guān)系，我們首先對歸一化的表達數(shù)據(jù)使用主成分分析（PCA）并鑒定了兩個主要的簇，稱為簇A（n = 339）和簇B（n = 83）。 （圖3a）?；诿芏鹊木哂性肼暤膽?yīng)用的空間聚類驗證了這兩個聚類的存在（補充圖8a）。野生型樣菌株（生物信息學(xué)評估的菌株具有所有主要調(diào)節(jié)基因的野生型等位基因）主要與簇A相關(guān)，除了10個異常菌株（圖3a）。通過使用Pilon（方法）重新檢查這些異常菌株的基因組數(shù)據(jù)，在這10個菌株中的8個中鑒定了covS全局轉(zhuǎn)錄調(diào)控基因中未檢測到的插入缺失。因此，轉(zhuǎn)錄組引導(dǎo)的多態(tài)性發(fā)現(xiàn)確定了這些異常轉(zhuǎn)錄組的潛在遺傳原因。
考慮到8種covS突變體菌株的轉(zhuǎn)錄組與野生型菌株顯著不同，并且與之前的結(jié)果53？8一致，我們假設(shè)在特定主要全局調(diào)節(jié)因子中具有突變的菌株可能具有可被利用的基因表達的不同潛在模式以區(qū)分特定類別的調(diào)節(jié)基因突變體。為了驗證這一假設(shè)，我們使用隨機森林機器學(xué)習(xí)59來確定四種類別標簽中的一種（即野生型，covR，covS或ropB突變體）是否可以高概率地分配給異常菌株。簡而言之，在分析283個單一菌株的轉(zhuǎn)錄組譜（參見方法）的基礎(chǔ)上，使用隨機森林分類來預(yù)測8個異常菌株的類別標記。基于轉(zhuǎn)錄組的分類正確地將所有八種生物鑒定為covS突變菌株（補充表11）。在81個covRS和21個ropB菌株中，分別有85.2和61.9％的菌株被準確分類（補充表11）。 covRS菌株被錯誤分類為野生型表型（轉(zhuǎn)錄物譜）類似于野生型菌株，與群A菌株分組（圖3c）。因此，轉(zhuǎn)錄物譜的機器學(xué)習(xí)分類準確地預(yù)測了基因型（調(diào)節(jié)基因突變狀態(tài)）并預(yù)測了具有主要調(diào)節(jié)基因突變的菌株的轉(zhuǎn)錄物表型（突變樣或野生型）。

調(diào)節(jié)基因突變和轉(zhuǎn)錄組變化。

將異常菌株重新分配為covS突變體導(dǎo)致簇A具有野生型和突變菌株，而簇B僅由突變菌株組成（圖3b）。檢查群集A和B的轉(zhuǎn)錄組學(xué)和基因組數(shù)據(jù)產(chǎn)生了五個發(fā)現(xiàn)。

• 首先，所有簇B菌株在covS或covR中具有突變，
• 其次，具有covS（68.5％）或covR（37.5％）突變的大多數(shù)菌株被分配至簇B（圖3c）。
• 第三，分配給簇A的大多數(shù)菌株是類似野生型或在除了covRS之外的主要調(diào)節(jié)基因中具有突變（如下所述）。
• 第四，群集B菌株具有比群集A菌株顯著更多的差異表達基因（圖4）。
• 第五，關(guān)于兩個主要轉(zhuǎn)錄組簇，沒有明顯的簡單基因組亞片特異性關(guān)聯(lián)（補充圖8b）。

CovRS雙組分系統(tǒng)負調(diào)節(jié)15％轉(zhuǎn)錄組的表達，包括關(guān)鍵的毒力因子44。與此一致，CovRS的失活增強了毒力53,56。我們比較了由188個預(yù)測的野生型菌株組成的442個菌株的轉(zhuǎn)錄組，在covRS中具有不同類型突變的132個菌株，以及在其他主要調(diào)節(jié)基因中具有不同突變的122個菌株。盡管簇B僅含有covRS突變株（圖3c），但相當大比例的covS（20.4％）和covR（45.8％）突變株與簇A中的幾乎所有野生型樣株組合（圖3c）。 covRS突變株主要由兩個不同的轉(zhuǎn)錄組聚類組成的發(fā)現(xiàn)表明covRS中的多態(tài)性不相同，如前所述[57,58]，并且群集A covRS突變體與野生型菌株的分組表明某些多態(tài)性可能比其他多態(tài)性具有更少的功能性后果。僅A群（n = 33）和群集B（n = 83）中的covRS突變株的PCA將該分組重現(xiàn)為兩個不同的簇（圖5a）。接下來，我們使用基于距離的聚類來檢驗這樣的假設(shè)：PCA（圖5a）不明顯的額外亞結(jié)構(gòu)可能存在于轉(zhuǎn)錄組數(shù)據(jù)中（圖5b）。調(diào)查結(jié)果載于補充說明。
為了檢驗B群covRS菌株與A群covRS菌株相比具有獨特轉(zhuǎn)錄組的假設(shè)，我們檢測了兩組基因的轉(zhuǎn)錄，發(fā)現(xiàn)了142個差異表達的基因（補充表12）。這兩個簇在差異表達基因的互補性以及差異表達基因的改變的轉(zhuǎn)錄物變化（向上或向下）的大小上不同。許多（32％）差異表達的基因具有五倍或更高的改變的轉(zhuǎn)錄物水平，包括編碼毒力因子的關(guān)鍵CovRS調(diào)節(jié)的基因，例如SPN和SLO，Mga調(diào)節(jié)子，HasABC和SpeB（補充表12）。此外，與群集A covRS菌株相比，群集B covRS菌株具有顯著增加的移碼誘導(dǎo)插入缺失和無義突變的頻率（雙尾Fisher精確檢驗，P <0.0001），可能使該調(diào)節(jié)系統(tǒng)失活（喪失 - 功能突變）。
鑒于編碼B群中多個關(guān)鍵毒力因子的基因轉(zhuǎn)錄本增加，我們假設(shè)B群菌株比A群covRS突變菌株毒性更強。與該假設(shè)一致，通過使用小鼠感染模型分析來自每個簇的四種菌株的毒力，顯示簇B菌株導(dǎo)致顯著更高的死亡率（圖5c）和更大的病變，具有更多的組織破壞（圖5d）。比較ropB突變可變地影響SpeB半胱氨酸蛋白酶毒力因子的表達和活性的類似研究在補充說明，補充圖8d環(huán)和補充表13和14中給出。

單核苷酸插入顯著改變毒力。

分泌的R28蛋白是一種GAS毒力因子，已被研究作為潛在的候選疫苗60,61。該蛋白質(zhì)由位于ICE樣元件上的Spy1336 / R28基因編碼，注釋為差異區(qū)域2（參考文獻）62,63（圖6a）。該37.4kb的DNA區(qū)段與B組鏈球菌63的染色體中存在的區(qū)域> 99％相同。在我們對最初50個菌株的轉(zhuǎn)錄組研究中，我們觀察到大約三分之一的菌株表達低水平的Spy1336 / R28轉(zhuǎn)錄物，而三分之二的菌株表達高水平的Spy1336 / R28轉(zhuǎn)錄物。相鄰基因（Spy1337）具有相同的表達模式（圖6b）。 Spy1336 / R28和Spy1337轉(zhuǎn)錄物水平與遺傳結(jié)構(gòu)，地理位置，分離年份或MGE-50基因型之間沒有相關(guān)性。這一令人費解的發(fā)現(xiàn)促使我們使用SEER64,65對所有442株菌株的從頭組裝進行了全基因組關(guān)聯(lián)研究（GWAS），因為這442株我們有相關(guān)的轉(zhuǎn)錄組數(shù)據(jù)。根據(jù)Spy1336 / R28和Spy1337基因的轉(zhuǎn)錄水平，我們檢查了具有低或高轉(zhuǎn)錄物表型的菌株是否與任何遺傳事件（例如，SNP，插入缺失或重組）顯著相關(guān)。對于兩種表型（高和低轉(zhuǎn)錄水平），100％的重要 k-mers映射到Spy1336 / R28和Spy1337之間的基因間區(qū)域（補充圖10a），這導(dǎo)致鑒定位于Spy1336 / R28之間的基因間區(qū)域的poly(T)同聚物區(qū)域中的變體和Spy1337基因（圖6c）。與高轉(zhuǎn)錄物表型正相關(guān)的顯著k-mers具有10個T殘基，與高轉(zhuǎn)錄本表型負相關(guān)的那些具有9個T殘基。還通過50-和442-菌株數(shù)據(jù)集的表達數(shù)量性狀基因座（eQTL）分析66,67（補充圖10b）鑒定了10T變體與Spy1336 / R28和Spy1337的增加的轉(zhuǎn)錄物水平的關(guān)聯(lián)。
與親本菌株（9個T殘基）相比，同基因突變菌株（10個T殘基）具有顯著增加的Spy1336 / R28和Spy1337的轉(zhuǎn)錄水平（圖6d）;如在小鼠壞死性肌炎感染模型中評估的那樣，引起顯著更大的肉眼和微觀病變，以及更接近的死亡率（圖6e，f）;對人體多形核白細胞的離體殺傷具有顯著的抗性（P <0.05;圖6g）;并產(chǎn)生更多分泌的和細胞相關(guān)的Spy1336 / R28蛋白（圖6h）。通過對生長至中指數(shù)或早期穩(wěn)定期的同基因菌株的RNA-seq分析進一步證實由該同聚物核苷酸區(qū)中T殘基數(shù)目的變化引起的Spy1336 / R28和Spy1337的改變水平（補充表15）。因此，在該同聚物區(qū)中插入或缺失單個T殘基顯著改變了Spy1336 / R28和Spy1337的轉(zhuǎn)錄水平，轉(zhuǎn)錄組和菌株毒力。

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

全基因組序列和轉(zhuǎn)錄組數(shù)據(jù)顯示emm28菌株之間存在顯著差異，我們推斷這些基因組和轉(zhuǎn)錄組學(xué)變化可能導(dǎo)致毒力的顯著變化。為了驗證這一假設(shè)，我們在壞死性肌炎小鼠模型中評估了相對于SC1A參考菌株MGAS28426的50個emm28菌株（上述初始轉(zhuǎn)錄組研究中使用的相同系統(tǒng)發(fā)育不同菌株組;補充表7）的毒力16,39， 68。事實上，對于所有已知的主要調(diào)節(jié)基因，所有這些菌株（96％）都是野生型。作為群體，SC1A和SC1B菌株的毒力沒有顯著差異（圖7a，b）。與之形成鮮明對比的是，SC2A菌株的毒力顯著高于SC1A和SC1B菌株（圖7a，b）。我們假設(shè)SC2A菌株的毒力增加可能至少部分是由于nga，ifs和slo基因的表達顯著增加，從而導(dǎo)致SC2A生物分泌的SPN和SLO毒素的產(chǎn)生增加，如圖所示。其他GAS血清型1,6,15,16。與該假設(shè)一致，SC2A菌株具有比SC1A或SC1B菌株顯著更高的nga轉(zhuǎn)錄物水平和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啟動子的同基因突變體菌株在體外產(chǎn)生顯著較低的SPN活性（圖7d），并且在小鼠壞死性肌炎感染模型中引起顯著較低的死亡率和組織破壞（圖7e，f）。

討論

我們使用GAS作為模型病原體來研究emm28化膿性鏈球菌菌株中群體基因組學(xué)，轉(zhuǎn)錄組學(xué)和毒力之間的復(fù)雜相互作用。我們發(fā)現(xiàn)轉(zhuǎn)錄組變化的幅度（差異表達基因的數(shù)量）與總體基因組 - 遺傳距離之間沒有簡單的相關(guān)性（補充圖9和補充說明）?？偟膩碚f，我們的分析表明，涉及基于種群的菌株樣本的多種類型的高維數(shù)據(jù)的整體方法可以提供對不太容易通過較少綜合和綜合方法發(fā)現(xiàn)的病原體環(huán)境相互作用的新理解。通過利用轉(zhuǎn)錄組特征分析，我們發(fā)現(xiàn)常用生物信息學(xué)方法會遺漏~5％的在的全局調(diào)節(jié)因子中具有突變菌株。這一發(fā)現(xiàn)可作為調(diào)查基因組轉(zhuǎn)錄關(guān)系研究的警示。
我們開發(fā)了基于群體的多維基因組和轉(zhuǎn)錄組數(shù)據(jù)集，以及包括GWAS和eQTL在內(nèi)的統(tǒng)計方法，以確定負責改變Spy1336 / R28轉(zhuǎn)錄水平的分子事件 - 一種編碼R28蛋白毒力因子和候選疫苗的基因60,63 69。這些發(fā)現(xiàn)增加了細菌中一個新興的主題，即基因間區(qū)域看似適度的變化可以改變基因表達并具有適應(yīng)性15,70？4。關(guān)于如何增加Spy1336 / R28和Spy1337的轉(zhuǎn)錄水平可能增強毒力存在幾種可能性。一種可能性是改變的毒力表型單獨或主要由Spy1336 / R28的產(chǎn)生增加引起，Spy1336 / R28是已知的毒力因子60,75。確切地說，這種蛋白質(zhì)如何有助于毒力尚不清楚，盡管據(jù)報道它可以促進人類上皮細胞的粘附[60,76]。第二種可能性是Spy1337調(diào)節(jié)蛋白直接或間接地改變其自身和可能影響毒力的其他基因的轉(zhuǎn)錄。為了驗證這一假設(shè)，我們對含有9或10個Ts的同基因菌株進行了RNA-seq分析，發(fā)現(xiàn)在中指數(shù)期，只有兩個基因（Spy1336 / R28和Spy1337）顯著上調(diào)，而在早期 - 穩(wěn)定期。Spy1336 / R28和Spy1337以及33個其他基因被上調(diào)，165個被下調(diào)（補充表15）。與9T親本生物相比，10T同基因突變株中的三基因fruRBA操縱子（Spy0641，Spy0642和Spy0643編碼蛋白參與果糖利用）高度上調(diào)。 fruR或fruB基因的失活顯著降低了人全血中或多形核白細胞存在下突變菌株的存活率77。與該發(fā)現(xiàn)一致，10T同基因突變株在壞死性肌炎小鼠模型中具有顯著更高的毒力（圖6e，f），對人多形核中性粒細胞的殺傷具有顯著增強的抗性（圖6g），并產(chǎn)生更多的分泌物。和細胞相關(guān)的Spy1336 / R28蛋白質(zhì)（圖6h;圖6i中的模型）。
最簡單的假設(shè)解釋了這個基因間區(qū)域中一個T殘基的插入或缺失如何改變Spy1336 / R28和Spy1337的轉(zhuǎn)錄水平是由Spy1337編碼的轉(zhuǎn)錄調(diào)節(jié)因子直接與該基因間區(qū)域結(jié)合并同時增加兩個基因的轉(zhuǎn)錄。在該假設(shè)下，同聚物核苷酸區(qū)可以是：（i）是整個Spy1337共有結(jié)合位點的一部分或構(gòu)成，或（ii）位于兩個共有結(jié)合位點側(cè)翼的間隔區(qū)中。在第一種情況下，具有10 Ts（與9 Ts相比）的同聚物區(qū)域?qū)?gòu)成更好的共有結(jié)合位點，而在第二種情況下，間隔區(qū)域中存在10 Ts（與9 Ts相比）將位于側(cè)翼推定位置在DNA螺旋中更有利的空間取向的共有結(jié)合位點，用于結(jié)合Spy1337轉(zhuǎn)錄調(diào)節(jié)因子（圖6c）。正在進行其他研究來解決這個問題。
本研究以新穎的方式使用機器學(xué)習(xí)和eQTL分析。我們的轉(zhuǎn)錄組數(shù)據(jù)集促進了機器學(xué)習(xí)的使用，以分析和正確分類基于單獨的基因組序列數(shù)據(jù)分析錯誤識別的調(diào)節(jié)突變株。直到最近，大多數(shù)使用機器學(xué)習(xí)的病原微生物的研究都使用了基于DNA序列的數(shù)據(jù)，通常側(cè)重于預(yù)測對抗菌藥物的抗性4,78？3。 eQTL分析已用于人和其他真核生物中的表達數(shù)據(jù)84？0;該研究將eQTL分析應(yīng)用于細菌數(shù)據(jù)集，通過生成的廣泛轉(zhuǎn)錄組數(shù)據(jù)實現(xiàn)。我們發(fā)現(xiàn)，基因間區(qū)域的插入缺失與五個基因的表達改變顯著相關(guān) - 兩個在順式（Spy1336 / R28和Spy1337）和三個在反式（Spy1338，Spy1339和Spy1340） - 通過使用50的轉(zhuǎn)錄數(shù)據(jù)在指數(shù)階段的應(yīng)變。類似地，通過使用來自早期穩(wěn)定期的442個菌株的轉(zhuǎn)錄物數(shù)據(jù)，發(fā)現(xiàn)順式（Spy1336 / R28和Spy1337）和與47個額外基因的反式關(guān)聯(lián)（錯誤發(fā)現(xiàn)率<0.0005）（補充圖10b和補充表16）。重要的是，通過eQTL分析鑒定的49個基因中的60％也通過等基因（10T）突變體和（9T）親本菌株的RNA-seq分析差異表達。
轉(zhuǎn)錄組數(shù)據(jù)還使我們得出結(jié)論，盡管HGT事件可以并確實改變了轉(zhuǎn)錄組;大多數(shù)轉(zhuǎn)錄組的變化是由SNP（錯義或無義突變）和影響主要調(diào)節(jié)基因如covR / covS，ropB和mga的短插入引起的，并導(dǎo)致同源編碼蛋白的截短53,5,57。這些發(fā)現(xiàn)與觀察結(jié)果一致，即這些調(diào)節(jié)基因是GAS菌株群體基因組分析中具有最高密度多態(tài)性的基因之一[6,91,92]。
由于未知原因，emm28 GAS菌株在產(chǎn)褥期敗血癥（兒童發(fā)燒），女性生殖道感染和新生兒感染17,31？4中的比例過高。雖然我們的研究并非旨在解決細菌種群結(jié)構(gòu)與感染菌株的詳細臨床表型之間非常復(fù)雜的關(guān)系，但一項觀察值得評論。對于美國患者的951株分離株，可獲得合理詳細的感染類型信息。我們發(fā)現(xiàn)，與來自美國的非SC2A菌株相比，來自美國的SC2A菌株的顯著更高比例與產(chǎn)褥敗血癥，新生兒感染和女性生殖道感染相關(guān)（2（1）= 5.854; P = 0.015;補充表1）。在這方面，我們注意到作為一組，SC2A菌株在小鼠壞死性肌炎實驗中也顯著更具毒性（圖7a，b）。
總之，我們的研究可以作為一個范例，說明如何有效整合基于種群的樣本生成的多維數(shù)據(jù)集，以產(chǎn)生關(guān)于微生物遺傳和病原體環(huán)境相互作用的新知識。與三種類型數(shù)據(jù)中的任何一種或兩種數(shù)據(jù)的研究相比，三種不同類型數(shù)據(jù)的整合導(dǎo)致對病原體的分子遺傳學(xué)的更多理解。該策略通常適用于任何微生物 - 致病性或其他 - 并可能導(dǎo)致新的治療方法。

方法

全基因組測序和多態(tài)性分析。

如前所述，使用Illumina NextSeq 550儀器進行菌株生長，染色體DNA的分離，成對末端文庫的產(chǎn)生和多重測序，如前所述[1,6,93]。用于細菌基因組分析的管道顯示在補充圖2a中。通過使用Musket94校正修剪的序列讀數(shù)，并通過使用SMALT將其定位至參考血清型M28菌株MGAS6180（GenBank登錄號CP000056）63的基因組（參見URL）。通過使用描述的FreeBayes（參見URL）和用于檢測插入缺失的Pilon95鑒定SNP和插入和缺失（插入缺失）。 SNPfx.pl（內(nèi)部開發(fā)的PERL腳本;參見URL）用于確定SNP的性質(zhì)（編碼與非編碼，同義與非同義，等等）?；蛘撸Y(jié)合使用，SPAdes算法用于從頭基因組裝配96。短讀序列分型2（SRST2）用于鑒定基因，等位基因和多位點序列類型97。算法Gubbins用于檢測HGT事件98。 hierBAPS99（群體結(jié)構(gòu)的分層貝葉斯分析;參見URL）用于確定種群結(jié)構(gòu)，SplitsTree用于估計系統(tǒng)發(fā)育樹和網(wǎng)絡(luò)100。通過使用范圍在50和200之間的簇的數(shù)量的先前上限值來運行hierBAPS，其具有估計算法的五個重復(fù)，每個運行收斂到群體結(jié)構(gòu)的相同后驗?zāi)Ｊ焦烙嫛?br> 通過使用具有R9.5流動池和快速條形碼試劑盒（SQK-RBK004）的Oxford Nanopore MinION儀器進行24個菌株的長讀序列（補充表19）。這些菌株選自50個菌株，其中RNA-seq表達分析一式三份進行，以數(shù)字代表主要的遺傳進化枝，包括多種原噬菌體和ICE含量（MGE基因型），并包括在R28中不同的菌株。啟動子區(qū)T-核苷酸同聚物束。如先前所述122，使用Unicycler101生成納米孔長讀數(shù)和Illumina短讀數(shù)據(jù)的混合組件。使用Trimmomatic進行讀取質(zhì)量過濾和修剪的質(zhì)量指標，使用Musket完成讀取錯誤糾正，使用SPAdes / Unicycler生成從頭組裝，使用SRST2進行多位點序列類型和emm測定，使用SMALT進行讀取映射，并且用FreeBayes完成多態(tài)性發(fā)現(xiàn)（補充表20）。
基于連鎖序列核心染色體SNP通過鄰近連接推斷菌株之間的系統(tǒng)發(fā)育，并且通過使用hierBAPS基于整個核心基因組序列的貝葉斯分析來定義相關(guān)菌株的進化枝。通過使用k-means（在下面和在補充說明中描述）選擇用于轉(zhuǎn)錄組分析的參考菌株。

用于轉(zhuǎn)錄組分析的菌株選擇。

為了從2,101個emm28菌株的測序群體中選擇一個分離子集，這些菌株對于給定大小的子集大致跨越盡可能多的遺傳變異，我們使用基于投影的方法，類似于細菌中使用的群體結(jié)構(gòu)校正GWAS方法SEER103。首先，通過使用核心SNP，分別計算屬于給定譜系（SC1A，SC1B，SC2A或SC2B）的所有emm28分離物的基于SNP的成對距離矩陣。我們排除了由Gubbins98鑒定的MGE和潛在的重組區(qū)域。隨后，我們按照每個譜系代表的總種群大小的比例對分離物進行比例取樣。然后使用每個譜系的距離矩陣將所有譜系分離物投射到具有多維縮放的三維歐幾里德空間中，如SEER103中所示。對于從譜系中選擇的子集的給定總大小k，使用具有200個隨機重啟104的k均值算法來識別k個質(zhì)心的最佳集合以跨越存在于三維多維縮放投影中的變化。譜系中存在遺傳變異。選擇具有到每個質(zhì)心的最小歐幾里德距離的分離物以確定k個代表性分離物的最終子集。

RNA-seq文庫制備和測序。

emm28菌株一式三份生長并在兩個時間點收獲。代表四種主要遺傳子片段的50個菌株（圖1a，b和補??充表7）在生長的中指數(shù)和早期 - 靜止期一式三份進行測定。根據(jù)制造商的說明用RNeasy試劑盒（Qiagen）提取RNA。核糖體RNA（rRNA）用革蘭氏陽性細菌（Illumina）的Ribo-Zero rRNA去除試劑盒耗盡，如前所述[1105]。分別用RNA Nano和Pico芯片（Agilent Technologies）和Agilent 2100生物分析儀評估總RNA和rRNA耗盡的RNA的質(zhì)量。用來自ScriptSeq Index PCR Primers試劑盒（Illumina）的索引反向引物制備互補DNA（cDNA）文庫，并用AMPure XP珠（Beckman Coulter）純化。用高靈敏度DNA芯片（Agilent Technologies）評估cDNA文庫的質(zhì)量。對于每個樣品，用Qubit dsDNA HS測定試劑盒（Invitrogen）熒光測定法測量cDNA文庫濃度。將cDNA文庫稀釋，合并，并用Illumina NextSeq儀器測序。該相同方案用于9T和10T同基因菌株的比較RNA-seq分析（補充表15）。
emm28菌株作為單一培養(yǎng)物生長并在一個時間點收獲。 461單體的轉(zhuǎn)錄組分析（補充表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。根據(jù)制造商的說明書，在1.5ml微量離心管和DynaMag-2磁體（Invitrogen）中，使用體積（164μl）的Agencourt RNAClean XP順磁珠純化片段化的去磷酸化的總RNA。最終的洗脫體積為12μl。通過條形碼的連接實現(xiàn)RNAtag-seq過程期間總RNA的匯集，使得來自相同菌株的所有RNA片段用單獨的條形碼明顯標記。我們使用了早期研究中描述的16種獨特的條形碼寡核糖核苷酸52（補充表17）。對于每個菌株，通過使用體積為20.11的T4 RNA連接酶1（ssRNA連接酶; New England Biolabs）將5μl片段化的磷酸化總RNA連接到1μl各自的寡核糖核苷酸上，終濃度為5μM。反應(yīng)在22℃下進行90分鐘。結(jié)扎后，通過加入59.9l RLT緩沖液（RNeasy Mini試劑盒; Qiagen）將每個樣品的體積增加至80μl，并與含有80μlRNA結(jié)合緩沖液（RNA Clean＆Concentrator-5; Zymo）的1：1混合物混合。研究）和在150ml微量離心管中的80μl100％乙醇。因此，通過將對應(yīng)于構(gòu)成一個特定池的8個菌株的總RNA依次通過一個Zymo柱并將它們濃縮在一起，為每組48個樣品制備6個樣品，每個樣品包含8個樣品，如補充圖7a所示。每個池的最終洗脫體積為32μl。
用RNA Pico芯片（Agilent Technologies）評估總RNA庫的質(zhì)量。我們制備了57個池，其中含有來自8個菌株的總RNA和來自5個菌株的總RNA的一個額外的池，總共461個菌株。 Ribo-Zero rRNA去除試劑盒（革蘭氏陽性細菌）用于從池中消除不需要的rRNA。通過使用RNA Pico芯片（Agilent Technologies）分析ribulplepleted RNA的質(zhì)量。如前所述進行第一鏈cDNA合成52。對于每個庫，將12μl核糖蛋白RNA與2μlAR2寡核苷酸（其與本研究中使用的所有16個條形碼寡核糖核苷酸中存在的區(qū)域互補（補充表17））混合，并在70℃下變性2分鐘。然后，使用AffinityScript逆轉(zhuǎn)錄酶（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的體積連接。將反應(yīng)物在22℃溫育過夜，然后使用2.5倍體積的Agencourt RNAClean XP珠子進行兩次連續(xù)的凈化反應(yīng)，并用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文庫濃度。
在方案的這一點上將樣品再合并一次。將對應(yīng)于四個庫的cDNA文庫以等摩爾量混合在一起，每個庫對應(yīng)于八個菌株。該過程再重復(fù)14次;因此，我們最終得到了15個超級池，共計58個池，代表461個菌株（補充圖7a）。對應(yīng)于每個超級池的文庫分別摻入10％PhiX文庫以改善簇多樣性并用Illumina NextSeq儀器測序。

RNA-seq數(shù)據(jù)分析。

用于處理RNA-seq數(shù)據(jù)的生物信息學(xué)管道在補充圖2b中給出。
分析50個emm28菌株，一式三份生長并在兩個時間點收獲。對于每次測序運行，Illumina bcl2fastq軟件（參見URL）用于將Illumina生成的BCL堿基調(diào)用文件轉(zhuǎn)換為FASTQ文件。通過使用FASTQC軟件評估測序數(shù)據(jù)的讀取質(zhì)量（參見URL）。使用Trimmomatic106進行適配器污染和讀取質(zhì)量過濾。通過使用EDGE-pro107將讀數(shù)定位于參考菌株MGAS6180的基因組，并且從隨后的分析中排除映射至rRNA和轉(zhuǎn)移RNA基因的讀數(shù)。此外，根據(jù)補充說明中描述的策略，將具有低表達的基因排除在下游分析之外。通過使用DESeq2進行差異表達分析（參考文獻108）。該相同的管道用于在Spy1336 / R28和Spy1337之間的同聚物區(qū)中含有9或10個Ts的同基因菌株的比較RNA-seq分析。
在一個時間點分析442個emm28單一菌株。使用bcl2fastq將從超級池的讀取解復(fù)用到單獨的池中。分別使用FASTQC和Trimmomatic完成了讀取質(zhì)量評估和讀取質(zhì)量修剪。通過使用FASTX-Toolkit（參見URL），根據(jù)內(nèi)聯(lián)條形碼將每個池中的讀取解復(fù)用為對應(yīng)于各個樣本的單獨fastq文件。每個池的讀數(shù)中位數(shù)為92.6百萬（補充圖7b），每個樣品每個樣品的讀數(shù)中位數(shù)在800萬到2400萬之間（補充圖7c）。補充說明中描述了低表達基因的排除。沒有發(fā)現(xiàn)顯著的批次效應(yīng)與表達數(shù)據(jù)相關(guān)聯(lián)。通過使用R包variancePartition109估計，發(fā)現(xiàn)批次效應(yīng)解釋的方差百分比<2％。
通過使用在NOISeq包110中實施的NOISeq-sim進行最終442個單一菌株組的差異表達分析。與參考菌株MGAS28737相比，在441個菌株的每一個中鑒定差異表達的基因，參考菌株MGAS28737如補充說明中所述選擇。 DESeq2（參考文獻108）用于差異表達分析，用于進行兩組比較的情況，其中每組包含一種以上的菌株。

機器學(xué)習(xí)與隨機森林分析。

使用函數(shù)TreeBagger（MATLAB R2016b，Statistics and Machine Learning Toolbox，MathWorks），使用MATLAB進行隨機森林分析。目的是根據(jù)菌株的轉(zhuǎn)錄組譜訓(xùn)練隨機森林59，將異常菌株分類為四類（covR，covS，ropB和野生型）。對隨機森林進行了訓(xùn)練，其中轉(zhuǎn)錄組數(shù)據(jù)針對僅在單一主要全球調(diào)控基因（covR，covS或ropB）中突變的菌株生成，并且已知所有已知主要調(diào)控基因的野生型菌株（n = 283）。因此，訓(xùn)練數(shù)據(jù)由283個菌株組成，其中已知1,614個基因和類別標簽上的轉(zhuǎn)錄組譜。測試數(shù)據(jù)由8個異常菌株組成，其中已知1,614個基因的轉(zhuǎn)錄組譜，但類別標記未知。在學(xué)習(xí)最終模型之前，應(yīng)用以下特征選擇過程。建立了所有基因的初始隨機森林（1,000棵樹），并通過使用內(nèi)置的特征重要性測量來估計基因的預(yù)測重要性。該過程重復(fù)十次，最終的特征重要性值作為各次運行的平均值。從最重要的功能開始，我們根據(jù)給定的重要順序依次包含更多功能。對于每個特征子集，我們進行了雙重交叉驗證，通過使用三分之二的訓(xùn)練數(shù)據(jù)構(gòu)建隨機森林（100棵樹），并對三分之一的訓(xùn)練數(shù)據(jù)進行評估。通過100次交叉驗證迭代的平均樣本外分類準確度來測量子模型的樣本外性能。隨著更多特征的增加，分類準確度的提高迅速趨于穩(wěn)定，并且在此基礎(chǔ)上，為最終模型選擇了十個最具信息性的基因。最終模型用于預(yù)測8個異常菌株的類概率（補充表11）。

GWAS分析。

使用SEER103的GWAS分析使用442個菌株的從頭組裝進行，其中我們還具有RNA-seq數(shù)據(jù)。 GWAS用于鑒定與二元表型分組顯著相關(guān)的遺傳變體（高轉(zhuǎn)錄物表達= 1;低轉(zhuǎn)錄物表達= 0），根據(jù)Spy1336 / R28基因的轉(zhuǎn)錄水平定義。繪制442個菌株的Spy1336 / R28基因的標準化轉(zhuǎn)錄物水平（計數(shù)），得到兩個視覺上非常不同的組 - 低（三分之一的菌株）和高表達者（三分之二的菌株） - 類似于我們對50株菌的觀察結(jié)果（圖6b）。我們確定了閾值，發(fā)現(xiàn)歸一化計數(shù)小于261.5的菌株被認為是低表達者（編碼為0），等于或大于261.5計數(shù)的菌株被認為是高表達者（編碼為1）。二元表型文件和de novo組裝的fasta文件提供給SEER，k-mers通過使用fsm-lite從組裝的讀數(shù)計算（參見URL）。為了考慮種群結(jié)構(gòu)，使用了由Mash111計算的距離矩陣。運行SEER分別產(chǎn)生17和13個顯著的k聚體（調(diào)整后的P值<10-8），其與高或低轉(zhuǎn)錄物表達呈正相關(guān)或負相關(guān)。

eQTL分析。

使用R包MatrixEQTL112進行eQTL分析。對于eQTL分析，通過使用前十個主成分作為協(xié)變量，在模型中考慮了種群結(jié)構(gòu)。如果多態(tài)性在所考慮的基因的1kb內(nèi)，則認為關(guān)聯(lián)是順式的。
通過使用壞死性肌炎的小鼠模型對50種天然存在的和同基因型血清型M28 GAS突變株進行毒力研究。如前所述進行用血清型M28 GAS菌株進行的小鼠壞死性肌炎研究39,68。使用在初始RNA-seq實驗中使用的50種天然存在的菌株，包括毒力參考菌株MGAS28426。菌株MGAS28426用作毒力參考菌株，因為它在基因組上代表SC1A菌株，并且它用作分泌的SPN測定的野生型參考。制備每種菌株的冷凍原種，并通過在連續(xù)稀釋后制備從解凍的培養(yǎng)物中回收的菌落形成單位（c.f.u.）進行定量。將免疫活性的4周齡雌性遠交CD1小鼠（Envigo）隨機分配到菌株處理組，并在5'08c.f.u的右下肢接種。每種細菌菌株（n = 20mice perstrain;總共1,000只小鼠）。該劑量是在兩個試驗實驗的基礎(chǔ)上選擇的，該實驗表明毒力參考菌株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的劑量。獲得了從分配給組織病理學(xué)分析的小鼠取得的肢體的代表性顯微圖像和顯微圖像。用于產(chǎn)生同基因突變體的寡核苷酸顯示在補充表18（補充說明）中。
所有動物研究均按照由Methodist Hospital Research Institute的機構(gòu)護理和使用委員會審查和批準的方案（AUP0615？041）進行。每天至少監(jiān)測一次小鼠，并根據(jù)國家研究委員會（美國）2011年實驗動物護理和使用指南更新委員會和護理指南提供的國際公認標準指南確定接近死亡率。和實驗動物的使用113。存活數(shù)據(jù)表示為Kaplan璏曲線，并且用對數(shù)秩檢驗（Prism6; GraphPad）確定統(tǒng)計學(xué)顯著性差異。

統(tǒng)計分析。

除非另有說明，誤差條代表s.d.，Pvalues用Fisher精確，Mann璚hitney U或?qū)?shù)秩檢驗計算。 MatrixEQTL程序包報告使用了錯誤發(fā)現(xiàn)率。貝葉斯聚類用于定義emm28群體中的進化枝和亞進化枝。 k-means和基于距離的聚類用于識別由PCA產(chǎn)生的應(yīng)變的二維空間聚類中的質(zhì)心，并分別在聚類中找到另外的子結(jié)構(gòu)。隨機森林分析是在MATLAB中使用函數(shù)TreeBagger完成的。 R包variancePartition用于確認沒有顯著的批次效應(yīng)。使用確定系數(shù)（R2）統(tǒng)計來研究遺傳距離和轉(zhuǎn)錄組重構(gòu)程度之間是否存在相關(guān)性。使用比例的onetailed測試確定II組RopB菌株含有顯著比例的影響功能結(jié)構(gòu)域的突變。 Pearson相關(guān)系數(shù)（r）用于比較來自一式三份分析的菌株的RNA-seq數(shù)據(jù)與從作為單體生長的菌株收集的RNAtag-seq數(shù)據(jù)。

入藏代碼。

本文報道的序列已經(jīng)保藏在國家生物技術(shù)信息中心（NCBI）序列閱讀檔案（見URL）中，其中BioProject登錄號為PRJNA434389，NCBI基因表達綜合登記號為GSE113058。

倫理。

所有小鼠研究均按照休斯頓衛(wèi)理公會研究所的Institutional Animal Care and Use Committee批準的方案（AUP0318？016）進行。所有具有人血液和血液成分的研究均按照國家過敏和傳染病研究所的人類受試者機構(gòu)審查委員會批準的方案（01-I-N055）進行。所有研究志愿者都簽署了書面知情同意書報告摘要。