Genome-wide estimation of recombination, mutation and positive selection enlightens diversification drivers of Mycobacterium bovis

1 Vol.:(0123456789) Scien i ic Repo s | (2021) 11:18789 | h ps://doi.o g/10.1038/s41598-021-98226-y www.na u e.com/scien i ic epo s Genome‑wide es ima ion o ecombina ion, mu a ion and posi i e selec ion enligh ens di e si ica ion d i e s o Mycobac e ium bo is Ana C. Reis1,2 & Mónica V. Cunha1,2* Genome sequencing has ein igo a ed he in ec ious disease esea ch ield, shedding ligh on disease epidemiology, pa hogenesis, hos –pa hogen in e ac ions and also e olu iona y p ocesses exe ed upon pa hogens. Mycobac e ium ube culosis complex (MTBC), enclosing M. bo is as one o i s animal‑adap ed membe s causing ube culosis (TB) in e es ial mammals, is a pa adigma ic model o bac e ial e olu ion. As o he MTBC membe s, M. bo is is pos ula ed as a s ic ly clonal, slowly e ol ing pa hogen, wi h appa en ly no signs o ecombina ion o ho izon al gene ans e . In his wo k, we applied compa a i e genomics o a whole genome sequence (WGS) da ase composed by 70 M. bo is om di e en lineages (Eu opean and A ican) o gain insigh s in o he e olu iona y o ces ha shape gene ic di e si ica ion in M. bo is. Th ee dis inc app oaches we e used o es ima e signs o ecombina ion. Globally, a small numbe o ecombinan e en s was iden i ied and con i med by wo independen me hods wi h solid suppo . S ill, ecombina ion e eals a weake e ec on M. bo is di e si y compa ed wi h mu a ion (o e all /m = 0.037). The di e en ial /m a e age alues ob ained ac oss he clonal complexes o M. bo is in ou da ase a e consis en wi h he gene al no ion ha he ex en o ecombina ion may a y widely among lineages assigned o he same axonomical species. Based on his wo k, ecombina ion in M. bo is canno be excluded and should hus be a opic o u he e o in u u e compa a i e genomics s udies o which WGS o la ge da ase s om di e en epidemiological scena ios ac oss he wo ld is c ucial. A smalle M. bo is da ase (n = 42) om a mul i‑hos TB endemic scena io was hen subjec ed o addi ional analyses, wi h he iden i ica ion o mo e han 1,800 si es whe ein a leas one s ain showed a single nucleo ide polymo phism (SNP). The majo i y (87.1%) was loca ed in coding egions, wi h he global a io o non‑synonymous upon synonymous al e a ions (dN/dS) exceeding 1.5, sugges ing ha posi i e selec ion is an impo an e olu iona y o ce exe ed upon M. bo is. A highe pe cen age o SNPs was de ec ed in genes en iched in o “lipid me abolism”, “cell wall and cell p ocesses” and “in e media y me abolism and espi a ion” unc ional ca ego ies, e ealing hei unde lying impo ance in M. bo is biology and e olu ion. A close look on genes p one o ho izon al gene ans e in he MTBC ances o and included in he 3R (DNA epai , eplica ion and ecombina ion) sys em e ealed a global a e age nega i e alue o Taijima’s D neu ali y es , sugges ing ha pas selec i e sweeps and popula ion expansion a e a ecen bo leneck emain as majo e olu iona y d i e s o he obliga o y pa hogen M. bo is in i s s uggle wi h he hos . The Mycobac e ium ube culosis complex (MTBC) is one o he mos success ul axon o bac e ial pa hogens and a pa adigma ic case in bac e ial e olu ion, e ealing a s ikingly high nucleo ide iden i y a he genome le el (> 99%) among i s membe s1,2. The di e en MTBC eco ypes cause ube culosis (TB), an in ec ious g anu- loma ous disease, in a b oad g oup o hos species, anging om mic o-mammals o humans3–5. Cu en ly, he OPEN 1Cen e o Ecology, E olu ion and En i onmen al Changes (cE3c), Faculdade de Ciências, Uni e sidade de Lisboa, Campo G ande, C2, Room 2.4.11, 1749-016 Lisbon, Po ugal. 2Biosys ems and In eg a i e Sciences Ins i u e (BioISI), Faculdade de Ciências da Uni e sidade de Lisboa, Lisbon, Po ugal. *email: [email p o ec ed] 2 Vol:.(1234567890) Scien i ic Repo s | (2021) 11:18789 | h ps://doi.o g/10.1038/s41598-021-98226-y www.na u e.com/scien i ic epo s/ complex encompasses human [M. ube culosis (M b), M. a icanum] and animal-adap ed pa hogens (M. bo is, M. cap ae, M. pinnipedii, M. mic o i, M. mungi, M. o ygis, M. su ica ae, “chimpanzee bacillus” and “dassie bacillus”)5,6. M. cane ii (also known as “smoo h ube cle bacilli”) has an a e age nucleo ide iden i y o 98% wi h he a o emen ioned mycobac e ia and compa a i e genomic wo ks sugges ha M. cane ii and he es o MTBC ha e di e ged e y ecen ly om a common ances o 7. Conside ing his no ion, se e al au ho s e e o M. cane ii as an MTBC membe 8. The MTBC has been sys ema ically desc ibed as a s ic ly clonal complex, wi h popula ion s uc u e being appa en ly domina ed by educ ions in di e si y, bo lenecks, selec i e sweeps and gene ic d i s9,10. Assuming he s ic ly clonal e olu ion o he complex, polymo phisms such as dele ions canno be es o ed by ecombina ion9. Based on his p emise, he successi e e en s o genomic dele ions o he egions o di e ence (RD) and TbD1 (M b speci ic dele ion 1 egion) ha e been p oposed as molecula ma ke s o MTBC e olu ion2,5,11. Compa a- i e genomics and whole genome sequencing (WGS) wo ks suppo he di ision o human-adap ed membe s in o nine lineages (M. ube culosis L1 o L4, L7 and L8; and M. a icanum L5, L6 and L9), wi h lineages L2 o L4 sha ing he dele ion o TbD1 egion2,11–13. Mo eo e , animal-adap ed membe s ha e been p oposed o sha e a common ances o and a e de ined by clade-speci ic dele ions in he RD7, RD8, RD9 and RD102,5,14. E en s o ho izon al gene ans e (HGT) and ecombina ion a e assumed o be a e and o ha e occu ed in he ances o s o MTBC, a he han h oughou he di e ging his o y o MTBC membe s15–17. Two ea ly epo s by Hughes and collabo a o s (2002) and Gu acke and collabo a o s (2006) sugges ed ha ecombina ion e en s migh ha e helped o shape he polymo phisms ma king speci ic loci o M. ube culosis s ains18,19. The appa - en absence o ecombina ion in MTBC has been a ibu ed o: (1) loss o mechanis ic p ocesses and abili y o HGT; (2) a eness o HGT e en s; and (3) no oppo uni y o ecombina ion e en s wi hin MTBC ecological niches14,17. Mo e ecen ly, a ew Whole Genome Sequencing (WGS) s udies applied o MTBC s ains20 and M. bo is21 p o ided e idences o ecombina ion, wi h he i s sugges ing ha MTBC s ains equen ly exchange small DNA agmen s, bu because o he limi ed nucleo ide sequence a ia ion, hese e en s emain unno iced. Mycobac e ium bo is is he MTBC membe mos equen ly eco e ed om li es ock, mainly ca le, al hough i can also be isola ed om ee- anging and enced wildli e4,22–24. M. bo is e ol ed o i e main clonal complexes [Eu opean 1 (Eu1), Eu opean 2 (Eu2), Eu opean 3 (Eu3), A ican 1 (A 1) and A ican 2 (A 2)], de ined based on spoligo yping p o ile, speci ic dele ions and single nucleo ide polymo phisms (SNPs) in speci ic genes25–29. These clonal complexes e idence he di e si y s uc u e o M. bo is popula ion and associa ion wi h geog aphic egions. Fu he mo e, a ecen WGS wo k by Zimpel and collabo a o s (2020) de ised an M. bo is SNP-based phylogeny wi h o e 1900 genomes, which sugges ed he exis ence o a leas ou dis inc lineages in he wo ld (named Lb1 o Lb4), ha a e no en i ely conco dan wi h he p e ious de ined clonal complexes, al hough geog aphic speci ici ies may also be con i med30. These au ho s pe o med phylogene ic and molecula da ing di e gence analyses bu did no in es iga e ecombina ion30. P e ious wo ks employing di e en molecula echniques such as spoligo yping, MIRU-VNTR (Mycobac e- ial In e spe sed Repe i i e Uni -Va iable Numbe o Tandem Repea ) and, mo e ecen ly, SNP yping, e ealed a ce ain le el o gene ic di e si y among M. bo is s ains31–35. The di e en ia ion o gene ic a ian s has become a c ucial ool o s udy disease epidemiology, con ibu ing o gain insigh s in o pa hogenesis, i ulence and dis- ease ansmission. The a i al o WGS me hodologies opened he possibili y o shed ligh in o he e olu iona y d i e s exe ed upon M. bo is genomes du ing adap a ion and pe sis ence o di e en hos s and epidemiological scena ios. In his wo k, we ake ad an age o a compa a i e genomic analysis o a di e se M. bo is da ase (n = 70), including isola es om di e en clonal complexes o gain insigh s in o he e olu iona y p ocesses o M. bo is, speci ically add essing phylogene ic ela ionships and ecombina ion e en s. Complemen a y o his analysis, he sub-da ase o M. bo is isola es (n = 42) ob ained om a well cha ac e ized mul i-hos TB endemic egion in Po ugal31,36 was u he explo ed o in e he balance be ween he ela i e a es o nonsynonymous (dN) o synonymous (dS) nucleo ide subs i u ion, and he e olu iona y con ibu ion o speci ic g oups o genes e e ed o in he li e a u e as ha ing been acqui ed hough HGT by he MTBC ances o 37,38, as well as genes encoding 3R (DNA epai , eplica ion and ecombina ion) sys em componen s39. The genes p oposed o be acqui ed h ough HGT we e selec ed since hey may ep esen ancien polymo phisms, and so i is expec ed ha hey migh con ain a highe ac ion o synonymous al e a ions. The genes included in he 3R sys em we e selec ed since p e ious wo k pe o med wi h M. ube culosis s ains sugges a gene al nega i e/pu i ying selec ion ac ing upon hese genes and ha hey migh play an impo an ole in e olu ion39. Ano he objec i e o he wo k was o in e he p esence o ecombina ion e en s. Fo his pu pose, and conside ing ha ou da ase om Po ugal only had genomes included in Eu opean clonal complex 2 and s ains wi hou a clonal complex assigned, we decided o include publicly a ailable genomic da a o end up wi h ep esen a i es om all clonal complexes and o inc ease obus ness and b ead h o esul s. Me hodology Mycobac e ium bo is isola es da ase . Fo y- wo newly sequenced M. bo is genomes om an endemic mul i-hos TB scena io in Po ugal (de ails below), p e iously cha ac e ized om an epidemiological poin o iew36, we e a he cen e o his wo k. Conside ing ha he da ase om Po ugal only has ep esen a i es o Eu opean 2 clonal complex and s ains wi hou complex assigned, publicly a ailable whole genome sequencing da a was added in o de o enla ge he da ase wi h ep esen a i es om all M. bo is clonal complexes. The e o e, h ee sou ces o whole genome sequencing da a we e used in his wo k: comple e/d a genome assemblies up o a maximum o 10 sca olds deposi ed a NCBI (Na ional Cen e o Bio echnology In o ma ion) (n = 15 iso- la es); Illumina as q iles deposi ed a SRA (Sequence Read A chi e) ep esen a i e o M. bo is clonal complex di e si y (n = 12 isola es)30; and 42 newly sequenced genomes om Po ugal. Mycobac e ium bo is BCG (bacil- 3 Vol.:(0123456789) Scien i ic Repo s | (2021) 11:18789 | h ps://doi.o g/10.1038/s41598-021-98226-y www.na u e.com/scien i ic epo s/ lus Calme e-Gué in) was excluded om he NCBI sea ch. M. bo is AF2122/97 commonly used as e e ence genome was included in he da ase . Due o he public una ailabili y o whole genome sequences om ep e- sen a i es o A ican 1 clonal complex, and he low numbe s o genomes om ep esen a i e s ains o A 2 and Eu1, aw sequencing da a a ailable a SRA was used in hose cases. The wo k o Zimpel and collabo a o s (2020) helped in he iden i ica ion o genomes om he a o emen ioned clonal complexes and in he selec ion p ocess o M. bo is o include in he da ase . Fo Eu3, only one ype genome is desc ibed (B ange e al., 2020), hus he genome ha we included is he solo ep esen a i e o he Eu3 complex. Globally, he da ase included 70M. bo is isola ed om eigh hos species, dis ibu ed by 12 coun ies be ween 1985 and 2016. Thi y-six we e assigned as Eu2, se en as Eu1, one as Eu3, h ee as A 1, ou as A 2 and 19 we e no a ibu ed o any clonal complex (de ails below). De ailed in o ma ion abou he M. bo is used in his s udy (including accession numbe s) can be ound in Table1 and Supplemen a y Table1. Newly sequenced genomes (da ase om Po ugal). Fo y- wo newly sequenced M. bo is whole genomes o igi- na ing om animal TB ho spo s in Po ugal and sca e ing a pe iod o o e 12yea s we e a he cen e o his s udy, as he unde lying wildli e-li es ock disease sys em has been moni o ed egula ly31,36 (Supplemen a y Fig.1). These s ains we e isola ed om ca le (n = 14), ed dee (n = 16) and wild boa (n = 12) om 2003 o 2015, acco ding o he ensuing p ocedu e: animal issue samples we e pooled and p ocessed ollowing he p o- ocol guidelines ecommended in he OIE Manual o Te es ial Animals and inocula ed on o S oneb ink and Löwens ein-Jensen py u a e solid media and liquid medium. Cul u es we e incuba ed a 37°C and inspec ed weekly o g ow h o a minimum pe iod o 12weeks. Colonies we e di ec ly s o ed a glyce ol solu ion a -80ºC. The DNA o he WGS p ocedu e was ob ained a e a single in i o passage o o iginal a chi ed samples in mycobac e ia selec i e medium (Middleb ook 7H9, BD Diagnos ics). Fo ha pu pose, ozen cul u e s ocks we e e-cul u ed on Middleb ook 7H9 supplemen ed wi h 5% sodium py u a e and 10% ADS en ichmen (50g albumin, 20g glucose, 8.5g sodium chlo ide in 1 L wa e ) a 37°C. A e ou weeks’ g ow h, he cul u e medium was enewed, and he cul u es we e moni o ed egula ly un il g ow h was obse ed. Cells we e ha es ed by cen i uga ion, he pelle was esuspended in 500 µL phospha e bu e saline (PBS), hea -killed a 99°C du ing 30min, cen i uged, and he supe na an s o ed a -20°C un il WGS. All p ocedu es we e pe o med on a le el 3 biosecu i y acili y. WGS pai ed-end genomic lib a ies we e p epa ed wi h unique indexing o each DNA sample and sequenced using Illumina MiSeq (2 × 250 pb) (40 samples) and HiSeq (2 × 150 pb) ( wo isola es) echnology (Eu o ins Genomics, Ge many). The genomic DNA was sequenced using he Illumina Genome Analyse wi h he pai ed- end module a achmen and lib a ies we e cons uc ed wi h Nex e a XT DNA Lib a y P ep Ki om Illumina, acco ding o he manu ac u e ’s speci ica ions. Clonal complex assignmen . Conside ing he da a eco e ed om SRA (n = 12), he clonal complex iden i ica- ion was a ailable as me ada a o he co esponding publica ions30,41,43. When conside ing comple e genomes, wi h he excep ion o M. bo is AF2122/97 and M. bo is 3601 ha a e ecognized membe s o Eu1 and Eu3 clonal complexes, espec i ely25,29, whole genome alignmen wi h M. ube culosis H37R (NCBI accession NC_000962.3) was pe o med using MAFFT (Mul iple alignmen p og am o amino acid o nucleo ide sequences, e sion 7.458) wi h pa ame e –add agmen s48. Then, he p esence o he dele ions and/o SNP cha ac e is ic o he di e en clonal complexes was sea ched. The newly sequenced M. bo is (n = 42) and aw eads om d a assembly genomes (n = 3) we e aligned wi h e e ence genome M. ube culosis H37R ia SNP pipeline and he p esence o he dele ions and/o SNP cha ac e is ic o he di e en clonal complexes was sea ched. In o ma ion om he p esence/absence o cha ac e is ic dele ions and/o SNP and spoligo yping p o ile we e ga he ed o assign he genomic da a o he co esponding clonal complex. Fo ou d a assemblies i was no possible o in e he spoligo yping p o ile, and so hey we e included in he “wi hou complex” g oup. Bioin o ma ics analysis. The bioin o ma ics wo k low ollowed in his wo k s a ed om de no o assem- bly and map o e e ence s a egies, wi h he pu pose o explo e ecombina ion e en s and he polymo phisms o speci ic gene g oups. Figu e1 p o ides a lowcha o he s eps ollowed. Fo he ecombina ion analysis, all he genomes we e used o inc emen he obus ness o in e encesand he associa ed me ics. De no o genome assembly. In o de o mi iga e e o s in he gene a ion o genome consensus sequences, we i s ob ained de no o assemblies and, hen, he co e mul i-alignmen . The Unicycle pipeline, cu en ly a ail- able a h ps:// gi hub. com/ wick/ Unicy cle 49, was implemen ed o pe o m de no o assembly o 54 sequenced genomes (42 newly sequenced and 12 as q iles eco e ed om SRA). B ie ly, be o e de no o assembly, eads quali y analysis was pe o med in Fas QC e sion 0.11.7 (h ps:// gi hub. com/s- and e ws/ Fas QC), and whene e necessa y cleaned wi h T immoma ic e sion 0.36 (op ions “cu adap e and o he illumina-speci ic sequences om he ead” and “cu bases o he end o a ead, i bellow a h eshold quali y o 20” we e applied) (h p:// www. usade llab. o g/ cms/? page= imm oma ic)50. Then, SPAdes op imise 49 was used o genome assembly and Pilon e sion 1.1851 o pos -assembly op imiza ion. A conse a i e b idging mode was selec ed o a oid misassemble and he k-me size was sea ched and selec ed be ween 20 and 95% o ead leng h. Following SPAdes guidelines and conside ing eads’ size, con igs wi h less han 300bp we e emo ed and a 20 ead dep h co e age cu -o was es ablished52. In he de no o assembly s a egy, no genome egions, such as he highly epe i i e P oline- Glu ama e (PE) and P oline-P oline Glu ama e (PPE) pa alogous genes, we e emo ed. 4 Vol:.(1234567890) Scien i ic Repo s | (2021) 11:18789 | h ps://doi.o g/10.1038/s41598-021-98226-y www.na u e.com/scien i ic epo s/ M. bo is ID Clonal complex(a) Coun y Yea Hos species Re e ences Type o sequence Mb0220 w/o CC Po ugal 2003 Ca le 40 Newly sequenced Mb0261 Eu2 Po ugal 2006 Red dee 40 Newly sequenced Mb0601 Eu2 Po ugal 2007 Ca le 40 Newly sequenced Mb0769 Eu2 Po ugal 2008 Ca le 40 Newly sequenced Mb0783 Eu2 Po ugal 2008 Wild boa 40 Newly sequenced Mb0865 Eu2 Po ugal 2008 Ca le 40 Newly sequenced Mb0891 Eu2 Po ugal 2009 Red dee 40 Newly sequenced Mb0893 Eu2 Po ugal 2008 Wild boa 40 Newly sequenced Mb1317 Eu2 Po ugal 2010 Ca le 40 Newly sequenced Mb1339 Eu2 Po ugal 2010 Ca le 40 Newly sequenced Mb1458 w/o CC Po ugal 2010 Wild boa 40 Newly sequenced Mb1480 w/o CC Po ugal 2010 Ca le 40 Newly sequenced Mb1654 Eu2 Po ugal 2011 Ca le 40 Newly sequenced Mb1670 w/o CC Po ugal 2011 Red dee 40 Newly sequenced Mb1711 Eu2 Po ugal 2011 Red dee 40 Newly sequenced Mb1712 Eu2 Po ugal 2011 Red dee 40 Newly sequenced Mb1714 Eu2 Po ugal 2011 Ca le 40 Newly sequenced Mb1744 w/o CC Po ugal 2012 Wild boa 40 Newly sequenced Mb1746 Eu2 Po ugal 2012 Red dee 40 Newly sequenced Mb1758 Eu2 Po ugal 2012 Ca le 40 Newly sequenced Mb1769 Eu2 Po ugal 2012 Wild boa 40 Newly sequenced Mb1785 Eu2 Po ugal 2012 Red dee 40 Newly sequenced Mb1789 Eu2 Po ugal 2012 Ca le 40 Newly sequenced Mb1841 Eu2 Po ugal 2012 Ca le 40 Newly sequenced Mb1870 Eu2 Po ugal 2012 Wild boa 40 Newly sequenced Mb1915 Eu2 Po ugal 2013 Red dee 40 Newly sequenced Mb1948 w/o CC Po ugal 2013 Red dee 40 Newly sequenced Mb1960 Eu2 Po ugal 2013 Red dee 40 Newly sequenced Mb2026 Eu2 Po ugal 2013 Ca le 40 Newly sequenced Mb2043 Eu2 Po ugal 2013 Red dee 40 Newly sequenced Mb2067 Eu2 Po ugal 2013 Wild boa 40 Newly sequenced Mb2206 Eu2 Po ugal 2014 Ca le 40 Newly sequenced Mb2235 w/o CC Po ugal 2014 Red dee 40 Newly sequenced Mb2277 w/o CC Po ugal 2014 Red dee 40 Newly sequenced Mb2300 Eu2 Po ugal 2014 Wild boa 40 Newly sequenced Mb2310 Eu2 Po ugal 2015 Red dee 40 Newly sequenced Mb2313 Eu2 Po ugal 2015 Wild boa 40 Newly sequenced Mb2325 Eu2 Po ugal 2015 Red dee 40 Newly sequenced Mb2328 Eu2 Po ugal 2015 Red dee 40 Newly sequenced Mb2347 w/o CC Po ugal 2015 Wild boa 40 Newly sequenced Mb2395 Eu2 Po ugal 2015 Wild boa 40 Newly sequenced Mb2397 Eu2 Po ugal 2015 Wild boa 40 Newly sequenced Mb502499 A 1 Ghana NA Human 30,41 SRA deposi ed Mb502526 A 1 Ghana NA Human 30,41 SRA deposi ed Mb1203064 A 1 Ghana NA Human 30,41 SRA deposi ed Mb4117155 A 2 F ance NA Wild boa 30,42 SRA deposi ed Mb1791710 A 2 Tanzania NA Chimpanzee 30,43 SRA deposi ed Mb1791712 A 2 Tanzania NA Chimpanzee 30,43 SRA deposi ed Mb1792006 Eu1 USA 2006 Ca le 43 SRA deposi ed Mb1792127 Eu1 USA 2008 Ca le 43 SRA deposi ed Mb1792361 Eu1 USA 2013 Ca le 43 SRA deposi ed Mb7240242 Eu1 USA 2016 Ca le 43 SRA deposi ed Mb7240415 Eu1 USA 2014 Ca le 43 SRA deposi ed Mb1791984 Eu1 USA 2005 Ca le 43 SRA deposi ed MBE1 w/o CC Egyp 2014 Ca le NA assemble/d a genomes NCBI MBE3 w/o CC Egyp 2014 Ca le NA assemble/d a genomes NCBI Con inued 5 Vol.:(0123456789) Scien i ic Repo s | (2021) 11:18789 | h ps://doi.o g/10.1038/s41598-021-98226-y www.na u e.com/scien i ic epo s/ Table 1. Cha ac e is ics o Mycobac e ium bo is genomes used in his wo k. Eu1: Eu opean 1, Eu2: Eu opean 2, Eu3: Eu opean 3, A 1: A ican 1, A 2: A ican 2, and w/o CC: wi hou clonal complex. NA: non-a ailable in o ma ion. M. bo is ID Clonal complex(a) Coun y Yea Hos species Re e ences Type o sequence MBE4 w/o CC Egyp 2014 Ca le NA assemble/d a genomes NCBI MBE10 w/o CC Egyp 2015 Ca le NA assemble/d a genomes NCBI Mb0077 w/o CC Canada 2006 Elk NA assemble/d a genomes NCBI Mb0565 w/o CC Canada 2011 Ca le NA assemble/d a genomes NCBI BMR25 w/o CC Canada 1985 Bison NA assemble/d a genomes NCBI Mb3601 Eu3 F ance 2014 Ca le 29 assemble/d a genomes NCBI Mb0476 Eu2 Canada 2002 Ca le NA assemble/d a genomes NCBI MbSP38 Eu2 B azil 2010 Ca le 44 assemble/d a genomes NCBI Mb1595 w/o CC Ko ea 2012 Ca le 45 assemble/d a genomes NCBI Mb0030 w/o CC China NA NA 46 assemble/d a genomes NCBI Mb0001 Eu2 B azil 2015 Tapi us e es is NA assemble/d a genomes NCBI Mb0003 w/o CC India 1986 Ca le NA assemble/d a genomes NCBI Mb31150 A 2 Uganda NA Chimpanzee 30,47 assemble/d a genomes NCBI Figu e1. Bioin o ma ics wo k low ollowed in his s udy. 6 Vol:.(1234567890) Scien i ic Repo s | (2021) 11:18789 | h ps://doi.o g/10.1038/s41598-021-98226-y www.na u e.com/scien i ic epo s/ The quali y o de no o assemblies was assessed by QUAST pipeline (h p:// quas . sou c e o ge. ne / q uas . h ml), which p omo es he emapping o con igs wi h M. bo is AF2122/97 e e ence genome (NCBI accession numbe LT708304.1) (quali y pa ame e s p esen ed in Supplemen a y Table1). Genome map o e e ence. The FASTQ iles om he newly sequenced M. bo is ob ained om Illumina sequencing we e aligned wi h M. bo is AF2122/97 e e ence genome (LT708304.1) wi h he help o SNP pipe- line (h ps:// gi hub. com/ USDA- VS/ SNP). The s anda d il e ing pa ame e s o a ian quali y sco e ecalib a- ion we e applied acco ding o Genome Analysis Toolki (GATK)’s Bes P ac ices ecommenda ions53–55. Resul s we e il e ed using a minimum SAM ools quali y sco e o 150 and AC = 2. Reads we e also examined using K aken (h p:// ccb. jhu. edu/ so w a e/ k aken/) o exclude con amina ion. The SNP pipeline used o he map o sequence s a egy in ou wo k examines a se ies o de ining SNPs and a ge s also o exclude mixed in ec ion scena ios. Genome co e age by eads was supe io o 99% (Supplemen a y Table1). To a oid mapping e o s and alse SNPs, a a ian was il e ed ou i : (1) i was suppo ed by less han 20 eads, (2) i was ound in a equency o less han 0.9, (3) i was egis e ed in a leas one s ain bu also wi h a gap in a leas ano he s ain. SNPs and posi ions wi h mapping issues o alignmen p oblems we e isually alida ed wi h In eg a ed Genomics Viewe (IGV) e sion 2.4.19 (h p:// so w a e . b oad ins i u e. o g/ so w a e/ ig /)56. Since P oline-Glu ama e (PE) and P oline-P oline Glu ama e (PPE) genes a e highly epe i i e and pa o mul i- gene amilies, hey a e p one o mis eading by Illumina sequencing and mis-mapping and so a e p e e en ially emo ed om he bioin o ma ics wo k low o Mycobac e ium ube culosis complex membe s when a s a egy o map o sequence is used o con i m SNPs. We hus il e ed PE/PPE genes ou om he analysis, as well as indels. All SNPs we e g ouped in o unc ional ca ego ies acco ding wi h Bo ilis (h p:// g eno l i s . pas e u . / Bo iL is /). The SnpE pipeline (h ps:// pcing ola. gi hub . io/ SnpE /) was employed o in e SNP consequences (synonymous o non-synonymous al e a ions). A new da abase o M. bo is AF2122/97 genome (LT708304.1) was c ea ed. Global co e genome mul i-alignmen . The co e genome mul i-alignmen was pe o med wi h Pa snp 1.2, cu - en ly a ailable a h ps:// gi hub. com/ ma bl/ pa snp57, using he 69 comple e genomes/d a assemblies (wi h op ion -c) and M. bo is AF2122/97 (LT708304.1) as e e ence. Fou co e mul i-alignmen we e pe o med: including only membe s o Eu2 clonal complex (n = 37), including all membe s o Eu opean clonal complexes (n = 44), including a junc ion o Eu opean and A ican clonal complexes (n = 51), and including all M. bo is om his s udy (n = 70). The co e alignmen s gene a ed by Pa snp we e used o in e maximum-likelihood (ML) phylogene ic ees using RAxML, ia CIPRES Science Ga eway 3.3 (h p:// www. phylo. o g/)58, wi h 1000 boo s ap eplica ions. Es ima ion o ecombina ion e en s. The p esence o ecombina ion e en s was examined using h ee di e - en algo i hms and bioin o ma ics ools in pa allel: Spli sT ee4 so wa e, Gubbins (Genealogies Unbiased By ecomBina ions In Nucleo ide Sequences) pipeline and RDP4 (Recombina ion De ec ion P og am, e sion be a 4.101) so wa e. The spli decomposi ion me hod implemen ed in Spli sT ee4 4.15.1 (h p:// www. spli s ee. o g/)59 was imple- men ed o compu e un oo ed phylogene ic ne wo ks, which we e alida ed s a is ically using he Phi es , wi h a signi icance h eshold o p = 0.05. The co e mul i-alignmen s om Pa snp analysis we e used as inpu and he spli decomposi ion as ne wo k c i e ia was implemen ed. Gubbins pipeline 2.3.1 (h ps:// gi hub. com/ sange - pa ho gens/ gubbi ns60 was un using de aul pa ame e s, as ano he way o assess he impac o ecombina ion on M. bo is. The algo i hm implemen ed in he pipeline econs uc s he clonal genealogy ela ing he comple e genomes/d a assemblies o ou da ase and he e e ence genome (M. bo is AF2122/97, LT708304.1) o each o he ; and scans he posi ions o SNPs ac oss each b anch o he ee in o de o de ec clus e s o SNPs ha would indica e ecombina ion e en s. The null hypo hesis o b anch assumes he absence o any ecombina ion e en s, he e o e implying ha he SNPs occu ing on he b anch should be e enly dis ibu ed. The co e mul i-alignmen s om Pa snp and he bes sco ing ML ee om RAxML we e used as inpu iles. Finally, o con i m he ecombina ion e en s sugges ed by he Gubbins pipeline, six algo i hms (RDP61, GENECONV62, Boo scan63, Maxchi64, Chimae a65, and SiScan66) implemen ed in RDP467 we e applied o he co e mul i-alignmen s om Pa snp unde de aul se ings. We es ablished ha a leas h ee o he algo i hms implemen ed in RDP4 had o conco dan ly e idence a signi ican signal o alida e each ecombina ion e en . Conside ing ha bo h Gubbins and RDP so wa e seek ecombina ion signals by inspec ing he co e mul i- alignmen in windows o 500bp maximum, and o con i m ha he inclusion o PE/PPE genes in he de no o assembly p ocess did no in e e e wi h he ecombina ion signals ound, he neighbou hood o genes in which ecombina ion e en s we e iden i ied we e u he inspec ed h ough a syn eny analysis. Syn eny maps, using comple e genomes, we e cons uc ed wi h MAUVE—mul i-genome alignmen (h p:// da li ngl ab. o g/ mau e/ mau e. h ml) o exclude local genome ansloca ions o in e sions. Fu he mo e, a syn eny analysis wi h ami- noacidic sequences was pe o med ia Syn Tax webse e (h ps:// a cha ea. i2bc. pa is- saclay. / Syn T ax/) using comple e genomes. Gene di e si y analyses. The genome da ase ob ained om a mul i-hos TB sys em in Po ugal was subjec ed o deepe analyses wi h he objec i e o examine he polymo phisms in he genes e e ed in he li e a u e as ha ing been acqui ed h ough HGT by he MTBC ances o 37,38 and in he genes encoding 3R (DNA epai , eplica ion and ecombina ion) sys em componen s39. Gene sequences o he 42M. bo is, oge he wi h gene sequence om he e e ence genome (M. bo is AF2122/97, NC_002945.4), we e aligned using Clus alX 2.1 7 Vol.:(0123456789) Scien i ic Repo s | (2021) 11:18789 | h ps://doi.o g/10.1038/s41598-021-98226-y www.na u e.com/scien i ic epo s/ (h p:// www. clus al. o g/ clus al2/) and used as an inpu o he calcula ion o gene di e si y, nucleo ide di e si y (π) and Tajima’s D neu ali y es pa ame e s ia DnaSP 6.12.03 (h p:// www. ub. edu/ dnasp/). Resul s and discussion Global phylogene ic analysis. A Maximum Likelihood (ML) phylogene ic ee based on he 69M. bo is isola es and e e ence genome was ob ained (Fig.2A). This s a egy allows he gene a ion o a mo e obus ee, when compa ing wi h single gene based ees o mul i-locus based ees, ha do no cap u e he a iabili y ac oss he en i e genome and consequen ly p esen low in e -speci ic disc imina o y powe 68,69. The esul ing opology o he ML ee gene ally ag ees wi h clonal complex classi ica ion, wi h genomes o Eu2 clus e ing in one ee b anch and genomes o A 1 also clus e ing oge he (Fig.2A). Resul s a e also in ag eemen wi h he known M. bo is e olu iona y ela ionships ha p esen a la ge di ision be ween Eu1 membe s and a g oup com- posed by all he o he clonal complexes and genomes wi hou assigned clonal complex30. Small inconsis encies be ween clonal complex and he ela ionships obse ed a he phylogene ic ee can be explained by he ac ha clonal complexes a e desc ibed based on speci ic genomic egions, while he phylogene ic ee is based on co e genome mul i-alignmen ep esen ing he whole genomes. E idences o ecombina ion in Mycobac e ium bo is. Mycobac e ium ube culosis complex is desc ibed o ha e clonally e ol ed, and mos e idences accumula ed o e he yea s suppo he idea ha ongo- ing HGT and ecombina ion e en s do no occu a de ec able le els in he MTBC15,17,18. P e ious wo ks ha e sugges ed ha he e migh be limi ed ecombina ion among MTBC s ains20,21, while o he s we e no success ul o iden i y measu able ecombina ion e en s70,71. To e isi his issue wi h ocus on M. bo is, and unlike p e ious wo ks ha only accoun ed o M. ube culosis70,71; o ha accoun edMTBC as a whole, wi h ew M. bo is ep esen a i es20; o ha only conside ed a es ic M. bo is da ase 21, in his wo k a o al o 70 s ains, wi h ep esen a i es om all clonal complexes, was used o sc een o ecombina ion. The da ase was scaled in ou cumula i e le els: (1) Eu2 membe s, (2) all Eu opean clonal complexes membe s (i.e. Eu opean), (3) bo hEu opean and A ican clonal complexes (Eu + A ) and (4) he en i e da ase (encompassing he genomes ha a e no included in any o he clonal complexes al eady desc ibed). To in es iga e his pos ula e u he , a spli decomposi ion ne wo k was pe o med o assess o he absence o ecombina ion e en s be ween genomes, since his me hod enables he isualiza ion o ances al ela ionships be ween indi iduals and displays con lic ing phylogene ic signals. The p esence o cycles in he ne wo k (i.e. egions ha do no con e ge in o a single ee), was con i med in all ou da ase s unde analysis, howe e none was suppo ed s a is ically by he Phi es (Eu2, p = 0.0956; Eu opean, p = 0.1637; Eu + A p = 0.2774; en i e da ase p = 0.2451), p o iding poo e idence o he p esence o ecombina ion e en s (Fig.3A-D). Following his analysis, and conside ing he obse a ion o cycles in all ne wo ks, he econs uc ion algo- i hm implemen ed in Gubbins pipeline was applied in o de o econs uc he clonal genealogy and o pe o m a complemen a y es ima ion o he impac o ecombina ion in M. bo is genomes. A cumula i e numbe o ecombina ion e en s was in e ed wi h he majo i y occu ing in e minal b anches (i.e. occu ing in a single genome) (Table2). The me ics showed consis ency ac oss he da ase s and e ealed ha ecombina ion e en s occu ed wo hund ed o h ee hund ed imes less equen ly han mu a ions, once he ho/ he a pa ame e Figu e2. Maximum likelihood phylogene ic ee (GTR) buil based on co e-genome alignmen o M. bo is genomes be o e (A) and a e (B) he emo al o ecombina ion si es. B anch colo s ep esen M. bo is clonal complexes: pu ple o Eu opean 1, ed o Eu opean 2, blue o Eu opean 3, o ange o A ican 1 and g een o A ican 2. The ee is oo ed and d awn o scale wi h b anch leng hs measu ed as he numbe o subs i u ions pe si e. 8 Vol:.(1234567890) Scien i ic Repo s | (2021) 11:18789 | h ps://doi.o g/10.1038/s41598-021-98226-y www.na u e.com/scien i ic epo s/ ha ep esen s he ela i e a es o ecombina ion and poin mu a ion on a b anch p esen ed an a e age alue be ween 0.0037 and 0.0056 (Table3). Recen ly, a published wo k wi h 38 M. bo is s ains e idenced a highe ho/ he a alue ( ho/ he a = 0.1) han he one ob ained o his da ase 21, howe e he wo k by Pa ané and co-wo ke s used e e ence-based assemblies o in e ecombina ion pa ame e s, a p ocedu e de ail ha was al eady associ- a ed wi h en ichmen o pu a i e ecombina ion e en s a e minal b anches due o he assembly p ocedu e70. Following, he /m pa ame e , which ep esen s he a io o di e si y in oduced by ecombina ion and mu a ion, e ealed an a e age alue be ween 0.025 and 0.037, poin ing ha ecombina ion has a lowe o e all Figu e3. Visualiza ion o con lic ing phylogene ic signals a un oo ed phylogene ic ees by he spli decomposi ion me hod in Eu opean 2 genomes (n = 37) (A), in Eu opean genomes (n = 44) (B), in a combina ion o Eu opean and A ican genomes (n = 51) (C) and in he en i e da ase (n = 70) (D). Table 2. Numbe o ecombina ion e en s in e ed by he Gubbins pipeline and RDP4. Da ase No. Gubbins e en s (% in e minal b anches) No. RDP4 e en s (% in e minal b anches) Eu opean 2 (n = 37) 4 (50%) 1 (0%) Eu opean (n = 44) 5 (60%) 2 (0%) Eu opean and A ican (n = 51) 6 (66.7%) 2 (0%) En i e da ase (n = 70) 8 (75%) 3 (33.3%) Table 3. Recombina ion me ics ob ained h ough he Gubbins pipeline analysis. Da ase /m Rho/ he a Eu opean 2 (n = 37) 0.025 0.0037 Eu opean (n = 44) 0.034 0.0046 Eu opean and A ican (n = 51) 0.037 0.0056 En i e da ase (n = 70) 0.037 0.0044 9 Vol.:(0123456789) Scien i ic Repo s | (2021) 11:18789 | h ps://doi.o g/10.1038/s41598-021-98226-y www.na u e.com/scien i ic epo s/ e ec in M. bo is gene ic di e si y when compa ing wi h mu a ion (Table3). To make a b oad compa ison, he /m pa ame e was es ima ed using a simila me hodology o an MTBC da ase composed by 23 genomes, e ealing a mean alue o 0.48620, while o he 38M. bo is da ase o Pa ané and co-wo ke s21 i e idenced a mean alue o 0.98. In he i s s udy he e we e only wo M. bo is (M. bo is BCG and e e ence s ain) wi hin he 23 genomes included in he wo k, so he ob ained alue migh be biased by he o e ep esen a ion o M. ube culosis genomes. In he second epo , he M. bo is popula ion unde analysis was mainly eco e ed om Ame ican coun ies and li es ock hos s. In con as , in ou da ase , a highe numbe o geog aphic loca ions and hos species is ep esen ed, and genomes g ouped in o di e en clonal complexes wi h dis inc popula ion gene ic signa u es we e also used, enabling a deepe and wide popula ion knowledge. The di e en ial /m a e age alues ob ained wi h ou da ase a e consis en wi h he no ion ha he ex en o ecombina ion a y widely among lineages assigned o he same axonomical species, so hese esul s sugges ha M. bo is clonal complexes migh exhibi a di e en ial impac o ecombina ion, as also sugges ed by Didelo & Maiden72. Ne - e heless, enla ging signi ican ly his da ase wi h he inclusion o a highe numbe o M. bo is genomes would allow u he cla i ica ion o his poin . Bo h /m and ho/ he a pa ame e s p esen a iabili y among he ee b anches, a esul ha is in ag eemen wi h epo s conce ning o he bac e ial species72,73. Finally, o con i m he ecombina ion e en s iden i ied by Gubbins pipeline, he di e en co e mul i-align- men s we e also independen ly es ed in RDP4 so wa e wi h six di e en algo i hms. Globally, less han hal o he e en s iden i ied by Gubbins we e con i med by RDP4 (Tables4, 5). Conside ing he en i e da ase , h ee ecombina ion e en s we e con i med, wo in ol ing in e nal nodes and ano he one in ol ing a single genome in a e minal b anch and o which a clonal complex could no be assigned (Tables4, 5). The iden i ica ion o e en s in e minal b anches migh be a sign ha ecombina ion is s ill ongoing in con empo a y M. bo is s ains o he esul o misalignmen 70. In his pu a i e ecombina ion egion, ci ca 20% o posi ions ha e an unde ined nucleo ide (N), which can he e o e in luence he ecombina ion signal (Supplemen a y Fig.2). Mo eo e , his egion a ec s he s gene, encoding he 16S ibosomal RNA ha is expec ed o be highly conse ed, so his pu a- i e ecombina ion signal could be he esul o a sequencing e o o w ong alignmen . Whole genome alignmen be ween Mb0003 and M. bo is AF2122/97 was hus hen pe o med and he p esence o unde ined nucleo ides and o SNPs was con i med, so he likely issues ela ed o w ong alignmen did no a i e as a consequence o he bioin o ma ics p ocedu e implemen ed in his wo k. No gaps o unde ined nucleo ides we e iden i ied in he ecombina ion egions o in e nal nodes (Figs.4, 5). Wi h espec o hese e en s, one encompasses exclusi ely Eu2 genomes, a ec ing he pks12 gene ha encodes a p obable polyke ide syn hase; while he o he one is egis e ed ac oss Eu1 genomes and a ec s na X gene Table 4. De ailed in o ma ion conce ning he ecombina ion e en s iden i ied by Gubbins and RDP4 in he en i e da ase . Genome posi ions acco ding wi h M. bo is AF2122/97. Recombina ion e en Iden i ica ion Co e-alignmen posi ions Genome posi ions(a) Gene name Mb gene name Classi ica ion o gene unc ion M. bo is isola e ID #1 Gubbins 945,923–945,950 1,220,297–1,220,324 PE PGRS22 Mb1121 PE-PGRS amily p o ein Mb2026 #2 Gubbins; RDP4 1,176,674–1,177,221 1,475,305–1,475,975 s Mb5019 Ribosomal RNA 16S Mb0003 #3 Gubbins; RDP4 1,532,736–1,532,787 1,953,495–1,953,548 na X Mb1765c P obable ni a e educ- ase Na X Mb1792361 Mb7240415 #4 Gubbins 1,532,751–1,532,781 1,953,840–1,953,870 na X Mb1765c P obable ni a e educ- ase Na X Mb1792361 #5 Gubbins; RDP4 1,794,609–1,794,714 2,283,200–2,283,315 pks12 Mb2074c P obable polyke ide syn hase pks12 Mb0891 Mb1711 Mb1789 Mb1870 Mb1758 Mb2043 Mb1960 #6 Gubbins 1,794,627–1,794,780 2,283,713–2,285,136 pks12 Mb2074c P obable polyke ide syn hase pks12 Mb0003 #7 Gubbins 2,242,002–2,242,098 2,839,474–2,839,570 a A Mb2121 P obable Sec-independ- en p o ein anslocase memb ane-bound p o ein a A Mb0565 #8 Gubbins 3,244,551–3,244,556 4,003,420–4,003,425 espa Mb3646c Conse ed hypo he ical alanine and glycine ich p o ein Mb2043 Table 5. S a is ical alues associa ed wi h di e en algo i hms implemen ed in RDP4 o he con i med ecombina ion e en s. Recombina ion e en Alignmen posi ions RDP (p- alue) GENECONV (p- alue) Boo scan (p- alue) MaxChi (p- alue) Chimae a (p- alue) #2 1,176,674–1,177,221 7.524 × 10−22 1.871 × 10−20 1.004 × 10−15 9.926 × 10−05 9.753 × 10−05 #3 1,532,736–1,532,787 3.771 × 10−09 5.216 × 10−08 5.634 × 10−03 – – #5 1,794,609–1,794,714 1.338 × 10−11 2.324 × 10−10 6.200 × 10−12 – – 16 Vol:.(1234567890) Scien i ic Repo s | (2021) 11:18789 | h ps://doi.o g/10.1038/s41598-021-98226-y www.na u e.com/scien i ic epo s/ 80. Bo i sch, E. C. e al. pks5- ecombina ion-media ed su ace emodelling in Mycobac e ium ube culosis eme gence. Na . Mic obiol. 1, 15019 (2016). h ps:// doi. o g/ 10. 1038/ nmic obiol. 2015. 19. Acknowledgemen s This wo k was unded by Fundação pa a a Ciência e a Tecnologia, IP (FCT) / MCTES h ough na ional unds (PIDDAC) and co- unded by he Eu opean Regional De elopmen Fund (FEDER) o he Eu opean Union, h ough he Lisbon Regional Ope a ional P og am and he Compe i i eness and In e na ionaliza ion Ope a ional P og am o Po ugal 2020 o o he p og ams ha may succeed (p ojec ‘Colossus: Con ol O ube cuLOsiS a he wildli e/li es ock in e ace uSing inno a i e na U e-based Solu ions’, e . PTDC/CVT-CVT/29783/2017, LIS- BOA-01-0145-FEDER-029783, POCI-01-0145-FEDER-029783). S a egic unding o cE3c and BioISI Resea ch Uni s (UIDB/00329/2020 and UIDB/04046/2020) om FCT is acknowledged. ACR was suppo ed by FCT h ough a doc o al g an (PD/BD/128031/2016). Au ho con ibu ions M.V.C. concei ed his wo k and secu ed esou ces and unding. A.C.R. pe o med he bioin o ma ic analyses unde he guidance o M.V.C. and explo ed he da a unde MVC supe ision. A.C.R. w o e he i s d a o he manusc ip and M.V.C. c i ically e ised all e sions. Bo h au ho s app o ed he inal e sion. Compe ing in e es s The au ho s decla e no compe ing in e es s. Addi ional in o ma ion Supplemen a y In o ma ion The online e sion con ains supplemen a y ma e ial a ailable a h ps:// doi. o g/ 10. 1038/ s41598- 021- 98226-y. Co espondence and eques s o ma e ials should be add essed o M.V.C. Rep in s and pe missions in o ma ion is a ailable a www.na u e.com/ ep in s. Publishe ’s no e Sp inge Na u e emains neu al wi h ega d o ju isdic ional claims in published maps and ins i u ional a ilia ions. Open Access This a icle is licensed unde a C ea i e Commons A ibu ion 4.0 In e na ional License, which pe mi s use, sha ing, adap a ion, dis ibu ion and ep oduc ion in any medium o o ma , as long as you gi e app op ia e c edi o he o iginal au ho (s) and he sou ce, p o ide a link o he C ea i e Commons licence, and indica e i changes we e made. The images o o he hi d pa y ma e ial in his a icle a e included in he a icle’s C ea i e Commons licence, unless indica ed o he wise in a c edi line o he ma e ial. I ma e ial is no included in he a icle’s C ea i e Commons licence and you in ended use is no pe mi ed by s a u o y egula ion o exceeds he pe mi ed use, you will need o ob ain pe mission di ec ly om he copy igh holde . To iew a copy o his licence, isi h p:// c ea i eco mmons. o g/ licen ses/ by/4. 0/. © The Au ho (s) 2021