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

1 Vol.:(0123456789) Scientific Reports | (2021) 11:18789 | https://doi.org/10.1038/s41598-021-98226-y www.nature.com/scientificreports Genome‑wide estimation of recombination, mutation and positive selection enlightens diversification drivers of Mycobacterium bovis Ana C. Reis1,2 & Mónica V. Cunha1,2* Genome sequencing has reinvigorated the infectious disease research field, shedding light on disease epidemiology, pathogenesis, host–pathogen interactions and also evolutionary processes exerted upon pathogens. Mycobacterium tuberculosis complex (MTBC), enclosing M. bovis as one of its animal‑adapted members causing tuberculosis (TB) in terrestrial mammals, is a paradigmatic model of bacterial evolution. As other MTBC members, M. bovis is postulated as a strictly clonal, slowly evolving pathogen, with apparently no signs of recombination or horizontal gene transfer. In this work, we applied comparative genomics to a whole genome sequence (WGS) dataset composed by 70 M. bovis from different lineages (European and African) to gain insights into the evolutionary forces that shape genetic diversification in M. bovis. Three distinct approaches were used to estimate signs of recombination. Globally, a small number of recombinant events was identified and confirmed by two independent methods with solid support. Still, recombination reveals a weaker effect on M. bovis diversity compared with mutation (overall r/m = 0.037). The differential r/m average values obtained across the clonal complexes of M. bovis in our dataset are consistent with the general notion that the extent of recombination may vary widely among lineages assigned to the same taxonomical species. Based on this work, recombination in M. bovis cannot be excluded and should thus be a topic of further effort in future comparative genomics studies for which WGS of large datasets from different epidemiological scenarios across the world is crucial. A smaller M. bovis dataset (n = 42) from a multi‑host TB endemic scenario was then subjected to additional analyses, with the identification of more than 1,800 sites wherein at least one strain showed a single nucleotide polymorphism (SNP). The majority (87.1%) was located in coding regions, with the global ratio of non‑synonymous upon synonymous alterations (dN/dS) exceeding 1.5, suggesting that positive selection is an important evolutionary force exerted upon M. bovis. A higher percentage of SNPs was detected in genes enriched into “lipid metabolism”, “cell wall and cell processes” and “intermediary metabolism and respiration” functional categories, revealing their underlying importance in M. bovis biology and evolution. A closer look on genes prone to horizontal gene transfer in the MTBC ancestor and included in the 3R (DNA repair, replication and recombination) system revealed a global average negative value for Taijima’s D neutrality test, suggesting that past selective sweeps and population expansion after a recent bottleneck remain as major evolutionary drivers of the obligatory pathogen M. bovis in its struggle with the host. The Mycobacterium tuberculosis complex (MTBC) is one of the most successful taxon of bacterial pathogens and a paradigmatic case in bacterial evolution, revealing a strikingly high nucleotide identity at the genome level (> 99%) among its members1,2. The different MTBC ecotypes cause tuberculosis (TB), an infectious granulomatous disease, in a broad group of host species, ranging from micro-mammals to humans3–5. Currently, the OPEN 1Centre for Ecology, Evolution and Environmental Changes (cE3c), Faculdade de Ciências, Universidade de Lisboa, Campo Grande, C2, Room 2.4.11, 1749-016 Lisbon, Portugal. 2Biosystems and Integrative Sciences Institute (BioISI), Faculdade de Ciências da Universidade de Lisboa, Lisbon, Portugal. *email: [email protected] 2 Vol:.(1234567890) Scientific Reports | (2021) 11:18789 | https://doi.org/10.1038/s41598-021-98226-y www.nature.com/scientificreports/ complex encompasses human [M. tuberculosis (Mtb), M. africanum] and animal-adapted pathogens (M. bovis, M. caprae, M. pinnipedii, M. microti, M. mungi, M. orygis, M. suricattae, “chimpanzee bacillus” and “dassie bacillus”)5,6. M. canettii (also known as “smooth tubercle bacilli”) has an average nucleotide identity of 98% with the aforementioned mycobacteria and comparative genomic works suggest that M. canettii and the rest of MTBC have diverged very recently from a common ancestor7. Considering this notion, several authors refer to M. canettii as an MTBC member8. The MTBC has been systematically described as a strictly clonal complex, with population structure being apparently dominated by reductions in diversity, bottlenecks, selective sweeps and genetic drifts9,10. Assuming the strictly clonal evolution of the complex, polymorphisms such as deletions cannot be restored by recombination9. Based on this premise, the successive events of genomic deletions of the regions of difference (RD) and TbD1 (Mtb specific deletion 1 region) have been proposed as molecular markers of MTBC evolution2,5,11. Comparative genomics and whole genome sequencing (WGS) works support the division of human-adapted members into nine lineages (M. tuberculosis L1 to L4, L7 and L8; and M. africanum L5, L6 and L9), with lineages L2 to L4 sharing the deletion of TbD1 region2,11–13. Moreover, animal-adapted members have been proposed to share a common ancestor and are defined by clade-specific deletions in the RD7, RD8, RD9 and RD102,5,14. Events of horizontal gene transfer (HGT) and recombination are assumed to be rare and to have occurred in the ancestors of MTBC, rather than throughout the diverging history of MTBC members15–17. Two early reports by Hughes and collaborators (2002) and Gutacker and collaborators (2006) suggested that recombination events might have helped to shape the polymorphisms marking specific loci of M. tuberculosis strains18,19. The apparent absence of recombination in MTBC has been attributed to: (1) loss of mechanistic processes and ability for HGT; (2) rareness of HGT events; and (3) no opportunity for recombination events within MTBC ecological niches14,17. More recently, a few Whole Genome Sequencing (WGS) studies applied to MTBC strains20 and M. bovis21 provided evidences of recombination, with the first suggesting that MTBC strains frequently exchange small DNA fragments, but because of the limited nucleotide sequence variation, these events remain unnoticed. Mycobacterium bovis is the MTBC member most frequently recovered from livestock, mainly cattle, although it can also be isolated from free-ranging and fenced wildlife4,22–24. M. bovis evolved to five main clonal complexes [European 1 (Eu1), European 2 (Eu2), European 3 (Eu3), African 1 (Af1) and African 2 (Af2)], defined based on spoligotyping profile, specific deletions and single nucleotide polymorphisms (SNPs) in specific genes25–29. These clonal complexes evidence the diversity structure of M. bovis population and association with geographic regions. Furthermore, a recent WGS work by Zimpel and collaborators (2020) devised an M. bovis SNP-based phylogeny with over 1900 genomes, which suggested the existence of at least four distinct lineages in the world (named Lb1 to Lb4), that are not entirely concordant with the previous defined clonal complexes, although geographic specificities may also be confirmed30. These authors performed phylogenetic and molecular dating divergence analyses but did not investigate recombination30. Previous works employing different molecular techniques such as spoligotyping, MIRU-VNTR (Mycobacterial Interspersed Repetitive Unit-Variable Number of Tandem Repeat) and, more recently, SNP typing, revealed a certain level of genetic diversity among M. bovis strains31–35. The differentiation of genetic variants has become a crucial tool to study disease epidemiology, contributing to gain insights into pathogenesis, virulence and disease transmission. The arrival of WGS methodologies opened the possibility to shed light into the evolutionary drivers exerted upon M. bovis genomes during adaptation and persistence to different hosts and epidemiological scenarios. In this work, we take advantage of a comparative genomic analysis of a diverse M. bovis dataset (n = 70), including isolates from different clonal complexes to gain insights into the evolutionary processes of M. bovis, specifically addressing phylogenetic relationships and recombination events. Complementary to this analysis, the sub-dataset of M. bovis isolates (n = 42) obtained from a well characterized multi-host TB endemic region in Portugal31,36 was further explored to infer the balance between the relative rates of nonsynonymous (dN) to synonymous (dS) nucleotide substitution, and the evolutionary contribution of specific groups of genes referred to in the literature as having been acquired though HGT by the MTBC ancestor37,38, as well as genes encoding 3R (DNA repair, replication and recombination) system components39. The genes proposed to be acquired through HGT were selected since they may represent ancient polymorphisms, and so it is expected that they might contain a higher fraction of synonymous alterations. The genes included in the 3R system were selected since previous work performed with M. tuberculosis strains suggest a general negative/purifying selection acting upon these genes and that they might play an important role in evolution39. Another objective of the work was to infer the presence of recombination events. For this purpose, and considering that our dataset from Portugal only had genomes included in European clonal complex 2 and strains without a clonal complex assigned, we decided to include publicly available genomic data to end up with representatives from all clonal complexes and to increase robustness and breadth of results. Methodology Mycobacterium bovis isolates dataset. Forty-two newly sequenced M. bovis genomes from an endemic multi-host TB scenario in Portugal (details below), previously characterized from an epidemiological point of view36, were at the centre of this work. Considering that the dataset from Portugal only has representatives of European 2 clonal complex and strains without complex assigned, publicly available whole genome sequencing data was added in order to enlarge the dataset with representatives from all M. bovis clonal complexes. Therefore, three sources of whole genome sequencing data were used in this work: complete/draft genome assemblies up to a maximum of 10 scaffolds deposited at NCBI (National Center for Biotechnology Information) (n = 15 isolates); Illumina fastq files deposited at SRA (Sequence Read Archive) representative of M. bovis clonal complex diversity (n = 12 isolates)30; and 42 newly sequenced genomes from Portugal. Mycobacterium bovis BCG (bacil- 3 Vol.:(0123456789) Scientific Reports | (2021) 11:18789 | https://doi.org/10.1038/s41598-021-98226-y www.nature.com/scientificreports/ lus Calmette-Guérin) was excluded from the NCBI search. M. bovis AF2122/97 commonly used as reference genome was included in the dataset. Due to the public unavailability of whole genome sequences from representatives of African 1 clonal complex, and the low numbers of genomes from representative strains of Af2 and Eu1, raw sequencing data available at SRA was used in those cases. The work of Zimpel and collaborators (2020) helped in the identification of genomes from the aforementioned clonal complexes and in the selection process of M. bovis to include in the dataset. For Eu3, only one type genome is described (Branger etal., 2020), thus the genome that we included is the solo representative of the Eu3 complex. Globally, the dataset included 70M. bovis isolated from eight host species, distributed by 12 countries between 1985 and 2016. Thirty-six were assigned as Eu2, seven as Eu1, one as Eu3, three as Af1, four as Af2 and 19 were not attributed to any clonal complex (details below). Detailed information about the M. bovis used in this study (including accession numbers) can be found in Table1 and Supplementary Table1. Newly sequenced genomes (dataset from Portugal). Forty-two newly sequenced M. bovis whole genomes originating from animal TB hotspots in Portugal and scattering a period of over 12years were at the centre of this study, as the underlying wildlife-livestock disease system has been monitored regularly31,36 (Supplementary Fig.1). These strains were isolated from cattle (n = 14), red deer (n = 16) and wild boar (n = 12) from 2003 to 2015, according to the ensuing procedure: animal tissue samples were pooled and processed following the protocol guidelines recommended in the OIE Manual for Terrestrial Animals and inoculated onto Stonebrink and Löwenstein-Jensen pyruvate solid media and liquid medium. Cultures were incubated at 37°C and inspected weekly for growth for a minimum period of 12weeks. Colonies were directly stored at glycerol solution at -80ºC. The DNA for the WGS procedure was obtained after a single invitro passage of original archived samples in mycobacteria selective medium (Middlebrook 7H9, BD Diagnostics). For that purpose, frozen culture stocks were re-cultured on Middlebrook 7H9 supplemented with 5% sodium pyruvate and 10% ADS enrichment (50g albumin, 20g glucose, 8.5g sodium chloride in 1 L water) at 37°C. After four weeks’ growth, the culture medium was renewed, and the cultures were monitored regularly until growth was observed. Cells were harvested by centrifugation, the pellet was resuspended in 500 µL phosphate buffer saline (PBS), heat-killed at 99°C during 30min, centrifuged, and the supernatant stored at -20°C until WGS. All procedures were performed on a level 3 biosecurity facility. WGS paired-end genomic libraries were prepared with unique indexing of each DNA sample and sequenced using Illumina MiSeq (2 × 250 pb) (40 samples) and HiSeq (2 × 150 pb) (two isolates) technology (Eurofins Genomics, Germany). The genomic DNA was sequenced using the Illumina Genome Analyser with the pairedend module attachment and libraries were constructed with Nextera XT DNA Library Prep Kit from Illumina, according to the manufacturer’s specifications. Clonal complex assignment. Considering the data recovered from SRA (n = 12), the clonal complex identification was available as metadata of the corresponding publications30,41,43. When considering complete genomes, with the exception of M. bovis AF2122/97 and M. bovis 3601 that are recognized members of Eu1 and Eu3 clonal complexes, respectively25,29, whole genome alignment with M. tuberculosis H37Rv (NCBI accession NC_000962.3) was performed using MAFFT (Multiple alignment program for amino acid or nucleotide sequences, version 7.458) with parameter–addfragments48. Then, the presence of the deletions and/or SNP characteristic of the different clonal complexes was searched. The newly sequenced M. bovis (n = 42) and raw reads from draft assembly genomes (n = 3) were aligned with reference genome M. tuberculosis H37Rv via vSNP pipeline and the presence of the deletions and/or SNP characteristic of the different clonal complexes was searched. Information from the presence/absence of characteristic deletions and/or SNP and spoligotyping profile were gathered to assign the genomic data to the corresponding clonal complex. For four draft assemblies it was not possible to infer the spoligotyping profile, and so they were included in the “without complex” group. Bioinformatics analysis. The bioinformatics workflow followed in this work started from de novo assembly and map to reference strategies, with the purpose to explore recombination events and the polymorphisms of specific gene groups. Figure1 provides a flowchart of the steps followed. For the recombination analysis, all the genomes were used to increment the robustness of inferencesand the associated metrics. De novo genome assembly. In order to mitigate errors in the generation of genome consensus sequences, we first obtained de novo assemblies and, then, the core multi-alignment. The Unicycler pipeline, currently available at https:// github. com/ rrwick/ Unicy cler49, was implemented to perform de novo assembly for 54 sequenced genomes (42 newly sequenced and 12 fastq files recovered from SRA). Briefly, before de novo assembly, reads quality analysis was performed in FastQC version 0.11.7 (https:// github. com/sandre ws/ FastQC), and whenever necessary cleaned with Trimmomatic version 0.36 (options “cut adapter and other illumina-specific sequences from the read” and “cut bases off the end of a read, if bellow a threshold quality of 20” were applied) (http:// www. usade llab. org/ cms/? page= trimm omatic)50. Then, SPAdes optimiser49 was used for genome assembly and Pilon version 1.1851 for post-assembly optimization. A conservative bridging mode was selected to avoid misassemble and the k-mer size was searched and selected between 20 and 95% of read length. Following SPAdes guidelines and considering reads’ size, contigs with less than 300bp were removed and a 20 read depth coverage cut-off was established52. In the de novo assembly strategy, no genome regions, such as the highly repetitive ProlineGlutamate (PE) and Proline-Proline Glutamate (PPE) paralogous genes, were removed. 4 Vol:.(1234567890) Scientific Reports | (2021) 11:18789 | https://doi.org/10.1038/s41598-021-98226-y www.nature.com/scientificreports/ M. bovis ID Clonal complex(a) Country Year Host species References Type of sequence Mb0220 w/o CC Portugal 2003 Cattle 40 Newly sequenced Mb0261 Eu2 Portugal 2006 Red deer 40 Newly sequenced Mb0601 Eu2 Portugal 2007 Cattle 40 Newly sequenced Mb0769 Eu2 Portugal 2008 Cattle 40 Newly sequenced Mb0783 Eu2 Portugal 2008 Wild boar 40 Newly sequenced Mb0865 Eu2 Portugal 2008 Cattle 40 Newly sequenced Mb0891 Eu2 Portugal 2009 Red deer 40 Newly sequenced Mb0893 Eu2 Portugal 2008 Wild boar 40 Newly sequenced Mb1317 Eu2 Portugal 2010 Cattle 40 Newly sequenced Mb1339 Eu2 Portugal 2010 Cattle 40 Newly sequenced Mb1458 w/o CC Portugal 2010 Wild boar 40 Newly sequenced Mb1480 w/o CC Portugal 2010 Cattle 40 Newly sequenced Mb1654 Eu2 Portugal 2011 Cattle 40 Newly sequenced Mb1670 w/o CC Portugal 2011 Red deer 40 Newly sequenced Mb1711 Eu2 Portugal 2011 Red deer 40 Newly sequenced Mb1712 Eu2 Portugal 2011 Red deer 40 Newly sequenced Mb1714 Eu2 Portugal 2011 Cattle 40 Newly sequenced Mb1744 w/o CC Portugal 2012 Wild boar 40 Newly sequenced Mb1746 Eu2 Portugal 2012 Red deer 40 Newly sequenced Mb1758 Eu2 Portugal 2012 Cattle 40 Newly sequenced Mb1769 Eu2 Portugal 2012 Wild boar 40 Newly sequenced Mb1785 Eu2 Portugal 2012 Red deer 40 Newly sequenced Mb1789 Eu2 Portugal 2012 Cattle 40 Newly sequenced Mb1841 Eu2 Portugal 2012 Cattle 40 Newly sequenced Mb1870 Eu2 Portugal 2012 Wild boar 40 Newly sequenced Mb1915 Eu2 Portugal 2013 Red deer 40 Newly sequenced Mb1948 w/o CC Portugal 2013 Red deer 40 Newly sequenced Mb1960 Eu2 Portugal 2013 Red deer 40 Newly sequenced Mb2026 Eu2 Portugal 2013 Cattle 40 Newly sequenced Mb2043 Eu2 Portugal 2013 Red deer 40 Newly sequenced Mb2067 Eu2 Portugal 2013 Wild boar 40 Newly sequenced Mb2206 Eu2 Portugal 2014 Cattle 40 Newly sequenced Mb2235 w/o CC Portugal 2014 Red deer 40 Newly sequenced Mb2277 w/o CC Portugal 2014 Red deer 40 Newly sequenced Mb2300 Eu2 Portugal 2014 Wild boar 40 Newly sequenced Mb2310 Eu2 Portugal 2015 Red deer 40 Newly sequenced Mb2313 Eu2 Portugal 2015 Wild boar 40 Newly sequenced Mb2325 Eu2 Portugal 2015 Red deer 40 Newly sequenced Mb2328 Eu2 Portugal 2015 Red deer 40 Newly sequenced Mb2347 w/o CC Portugal 2015 Wild boar 40 Newly sequenced Mb2395 Eu2 Portugal 2015 Wild boar 40 Newly sequenced Mb2397 Eu2 Portugal 2015 Wild boar 40 Newly sequenced Mb502499 Af1 Ghana NA Human 30,41 SRA deposited Mb502526 Af1 Ghana NA Human 30,41 SRA deposited Mb1203064 Af1 Ghana NA Human 30,41 SRA deposited Mb4117155 Af2 France NA Wild boar 30,42 SRA deposited Mb1791710 Af2 Tanzania NA Chimpanzee 30,43 SRA deposited Mb1791712 Af2 Tanzania NA Chimpanzee 30,43 SRA deposited Mb1792006 Eu1 USA 2006 Cattle 43 SRA deposited Mb1792127 Eu1 USA 2008 Cattle 43 SRA deposited Mb1792361 Eu1 USA 2013 Cattle 43 SRA deposited Mb7240242 Eu1 USA 2016 Cattle 43 SRA deposited Mb7240415 Eu1 USA 2014 Cattle 43 SRA deposited Mb1791984 Eu1 USA 2005 Cattle 43 SRA deposited MBE1 w/o CC Egypt 2014 Cattle NA assemble/draft genomes NCBI MBE3 w/o CC Egypt 2014 Cattle NA assemble/draft genomes NCBI Continued 5 Vol.:(0123456789) Scientific Reports | (2021) 11:18789 | https://doi.org/10.1038/s41598-021-98226-y www.nature.com/scientificreports/ Table 1. Characteristics of Mycobacterium bovis genomes used in this work. Eu1: European 1, Eu2: European 2, Eu3: European 3, Af1: African 1, Af2: African 2, and w/o CC: without clonal complex. NA: non-available information. M. bovis ID Clonal complex(a) Country Year Host species References Type of sequence MBE4 w/o CC Egypt 2014 Cattle NA assemble/draft genomes NCBI MBE10 w/o CC Egypt 2015 Cattle NA assemble/draft genomes NCBI Mb0077 w/o CC Canada 2006 Elk NA assemble/draft genomes NCBI Mb0565 w/o CC Canada 2011 Cattle NA assemble/draft genomes NCBI BMR25 w/o CC Canada 1985 Bison NA assemble/draft genomes NCBI Mb3601 Eu3 France 2014 Cattle 29 assemble/draft genomes NCBI Mb0476 Eu2 Canada 2002 Cattle NA assemble/draft genomes NCBI MbSP38 Eu2 Brazil 2010 Cattle 44 assemble/draft genomes NCBI Mb1595 w/o CC Korea 2012 Cattle 45 assemble/draft genomes NCBI Mb0030 w/o CC China NA NA 46 assemble/draft genomes NCBI Mb0001 Eu2 Brazil 2015 Tapirus terrestris NA assemble/draft genomes NCBI Mb0003 w/o CC India 1986 Cattle NA assemble/draft genomes NCBI Mb31150 Af2 Uganda NA Chimpanzee 30,47 assemble/draft genomes NCBI Figure1. Bioinformatics workflow followed in this study. 6 Vol:.(1234567890) Scientific Reports | (2021) 11:18789 | https://doi.org/10.1038/s41598-021-98226-y www.nature.com/scientificreports/ The quality of de novo assemblies was assessed by QUAST pipeline (http:// quast. sourc e f orge. net/ q uas t. html), which promotes the remapping of contigs with M. bovis AF2122/97 reference genome (NCBI accession number LT708304.1) (quality parameters presented in Supplementary Table1). Genome map to reference. The FASTQ files from the newly sequenced M. bovis obtained from Illumina sequencing were aligned with M. bovis AF2122/97 reference genome (LT708304.1) with the help of vSNP pipeline (https:// github. com/ USDAVS/ vSNP). The standard filtering parameters or variant quality score recalibration were applied according to Genome Analysis Toolkit (GATK)’s Best Practices recommendations53–55. Results were filtered using a minimum SAMtools quality score of 150 and AC = 2. Reads were also examined using Kraken (http:// ccb. jhu. edu/ softw are/ kraken/) to exclude contamination. The vSNP pipeline used for the map to sequence strategy in our work examines a series of defining SNPs and targets also to exclude mixed infection scenarios. Genome coverage by reads was superior to 99% (Supplementary Table1). To avoid mapping errors and false SNPs, a variant was filtered out if: (1) it was supported by less than 20 reads, (2) it was found in a frequency of less than 0.9, (3) it was registered in at least one strain but also with a gap in at least another strain. SNPs and positions with mapping issues or alignment problems were visually validated with Integrated Genomics Viewer (IGV) version 2.4.19 (ht t p:// so ftw a r e . b road insti t ute. or g/ so ftw a r e/ igv/)56. Since Proline-Glutamate (PE) and Proline-Proline Glutamate (PPE) genes are highly repetitive and part of multigene families, they are prone to misreading by Illumina sequencing and mis-mapping and so are preferentially removed from the bioinformatics workflow of Mycobacterium tuberculosis complex members when a strategy of map to sequence is used to confirm SNPs. We thus filtered PE/PPE genes out from the analysis, as well as indels. All SNPs were grouped into functional categories according with Bovilist (h t t p:// g eno l i st. pas te ur. fr/ BoviL ist/). The SnpEff pipeline (htt ps:// pcing ola. github . io/ SnpEff/) was employed to infer SNP consequences (synonymous or non-synonymous alterations). A new database for M. bovis AF2122/97 genome (LT708304.1) was created. Global core genome multi-alignment. The core genome multi-alignment was performed with Parsnp v1.2, currently available at https:// github. com/ marbl/ parsnp57, using the 69 complete genomes/draft assemblies (with option -c) and M. bovis AF2122/97 (LT708304.1) as reference. Four core multi-alignment were performed: including only members of Eu2 clonal complex (n = 37), including all members of European clonal complexes (n = 44), including a junction of European and African clonal complexes (n = 51), and including all M. bovis from this study (n = 70). The core alignments generated by Parsnp were used to infer maximum-likelihood (ML) phylogenetic trees using RAxML, via CIPRES Science Gateway v3.3 (http:// www. phylo. org/)58, with 1000 bootstrap replications. Estimation of recombination events. The presence of recombination events was examined using three different algorithms and bioinformatics tools in parallel: SplitsTree4 software, Gubbins (Genealogies Unbiased By recomBinations In Nucleotide Sequences) pipeline and RDP4 (Recombination Detection Program, version beta 4.101) software. The split decomposition method implemented in SplitsTree4 v4.15.1 (http:// www. split s tr ee. o r g/)59 was implemented to compute unrooted phylogenetic networks, which were validated statistically using the Phi test, with a significance threshold of p = 0.05. The core multi-alignments from Parsnp analysis were used as input and the split decomposition as network criteria was implemented. Gubbins pipeline v2.3.1 (https:// github. com/ sangerpatho gens/ gubbi ns60 was run using default parameters, as another way to assess the impact of recombination on M. bovis. The algorithm implemented in the pipeline reconstructs the clonal genealogy relating the complete genomes/draft assemblies of our dataset and the reference genome (M. bovis AF2122/97, LT708304.1) to each other; and scans the positions of SNPs across each branch of the tree in order to detect clusters of SNPs that would indicate recombination events. The null hypothesis for branch assumes the absence of any recombination events, therefore implying that the SNPs occurring on the branch should be evenly distributed. The core multi-alignments from Parsnp and the best scoring ML tree from RAxML were used as input files. Finally, to confirm the recombination events suggested by the Gubbins pipeline, six algorithms (RDP61, GENECONV62, Bootscan63, Maxchi64, Chimaera65, and SiScan66) implemented in RDP467 were applied to the core multi-alignments from Parsnp under default settings. We established that at least three of the algorithms implemented in RDP4 had to concordantly evidence a significant signal to validate each recombination event. Considering that both Gubbins and RDP software seek recombination signals by inspecting the core multialignment in windows of 500bp maximum, and to confirm that the inclusion of PE/PPE genes in the de novo assembly process did not interfere with the recombination signals found, the neighbourhood of genes in which recombination events were identified were further inspected through a synteny analysis. Synteny maps, using complete genomes, were constructed with MAUVE—multi-genome alignment (http:// darli ngl ab. org/ mauve/ mauve. html) to exclude local genome translocations or inversions. Furthermore, a synteny analysis with aminoacidic sequences was performed via SyntTax webserver (https:// archa ea. i2bc. parissaclay. fr/ SyntT ax/) using complete genomes. Gene diversity analyses. The genome dataset obtained from a multi-host TB system in Portugal was subjected to deeper analyses with the objective to examine the polymorphisms in the genes referred in the literature as having been acquired through HGT by the MTBC ancestor37,38 and in the genes encoding 3R (DNA repair, replication and recombination) system components39. Gene sequences of the 42M. bovis, together with gene sequence from the reference genome (M. bovis AF2122/97, NC_002945.4), were aligned using ClustalX v2.1 7 Vol.:(0123456789) Scientific Reports | (2021) 11:18789 | https://doi.org/10.1038/s41598-021-98226-y www.nature.com/scientificreports/ (http:// www. clust al. org/ clust al2/) and used as an input for the calculation of gene diversity, nucleotide diversity (π) and Tajima’s D neutrality test parameters via DnaSP v6.12.03 (http:// www. ub. edu/ dnasp/). Results and discussion Global phylogenetic analysis. A Maximum Likelihood (ML) phylogenetic tree based on the 69M. bovis isolates and reference genome was obtained (Fig.2A). This strategy allows the generation of a more robust tree, when comparing with single gene based trees or multi-locus based trees, that do not capture the variability across the entire genome and consequently present low inter-specific discriminatory power68,69. The resulting topology of the ML tree generally agrees with clonal complex classification, with genomes of Eu2 clustering in one tree branch and genomes of Af1 also clustering together (Fig.2A). Results are also in agreement with the known M. bovis evolutionary relationships that present a large division between Eu1 members and a group composed by all the other clonal complexes and genomes without assigned clonal complex30. Small inconsistencies between clonal complex and the relationships observed at the phylogenetic tree can be explained by the fact that clonal complexes are described based on specific genomic regions, while the phylogenetic tree is based on core genome multi-alignment representing the whole genomes. Evidences of recombination in Mycobacterium bovis. Mycobacterium tuberculosis complex is described to have clonally evolved, and most evidences accumulated over the years support the idea that ongoing HGT and recombination events do not occur at detectable levels in the MTBC15,17,18. Previous works have suggested that there might be limited recombination among MTBC strains20,21, while others were not successful to identify measurable recombination events70,71. To revisit this issue with focus on M. bovis, and unlike previous works that only accounted for M. tuberculosis70,71; or that accountedMTBC as a whole, with few M. bovis representatives20; or that only considered a restrict M. bovis dataset21, in this work a total of 70 strains, with representatives from all clonal complexes, was used to screen for recombination. The dataset was scaled in four cumulative levels: (1) Eu2 members, (2) all European clonal complexes members (i.e. European), (3) bothEuropean and African clonal complexes (Eu + Af) and (4) the entire dataset (encompassing the genomes that are not included in any of the clonal complexes already described). To investigate this postulate further, a split decomposition network was performed to assess for the absence of recombination events between genomes, since this method enables the visualization of ancestral relationships between individuals and displays conflicting phylogenetic signals. The presence of cycles in the network (i.e. regions that do not converge into a single tree), was confirmed in all four datasets under analysis, however none was supported statistically by the Phi test (Eu2, p = 0.0956; European, p = 0.1637; Eu + Af p = 0.2774; entire dataset p = 0.2451), providing poor evidence for the presence of recombination events (Fig.3A-D). Following this analysis, and considering the observation of cycles in all networks, the reconstruction algorithm implemented in Gubbins pipeline was applied in order to reconstruct the clonal genealogy and to perform a complementary estimation of the impact of recombination in M. bovis genomes. A cumulative number of recombination events was inferred with the majority occurring in terminal branches (i.e. occurring in a single genome) (Table2). The metrics showed consistency across the datasets and revealed that recombination events occurred two hundred to three hundred times less frequently than mutations, once the rho/theta parameter Figure2. Maximum likelihood phylogenetic tree (GTR) built based on core-genome alignment of M. bovis genomes before (A) and after (B) the removal of recombination sites. Branch colors represent M. bovis clonal complexes: purple for European 1, red for European 2, blue for European 3, orange for African 1 and green for African 2. The tree is rooted and drawn to scale with branch lengths measured as the number of substitutions per site. 8 Vol:.(1234567890) Scientific Reports | (2021) 11:18789 | https://doi.org/10.1038/s41598-021-98226-y www.nature.com/scientificreports/ that represents the relative rates of recombination and point mutation on a branch presented an average value between 0.0037 and 0.0056 (Table3). Recently, a published work with 38 M. bovis strains evidenced a higher rho/ theta value (rho/theta = 0.1) than the one obtained for this dataset21, however the work by Patané and co-workers used reference-based assemblies to infer recombination parameters, a procedure detail that was already associated with enrichment of putative recombination events at terminal branches due to the assembly procedure70. Following, the r/m parameter, which represents the ratio of diversity introduced by recombination and mutation, revealed an average value between 0.025 and 0.037, pointing that recombination has a lower overall Figure3. Visualization of conflicting phylogenetic signals at unrooted phylogenetic trees by the split decomposition method in European 2 genomes (n = 37) (A), in European genomes (n = 44) (B), in a combination of European and African genomes (n = 51) (C) and in the entire dataset (n = 70) (D). Table 2. Number of recombination events inferred by the Gubbins pipeline and RDP4. Dataset No. Gubbins events (% in terminal branches) No. RDP4 events (% in terminal branches) European 2 (n = 37) 4 (50%) 1 (0%) European (n = 44) 5 (60%) 2 (0%) European and African (n = 51) 6 (66.7%) 2 (0%) Entire dataset (n = 70) 8 (75%) 3 (33.3%) Table 3. Recombination metrics obtained through the Gubbins pipeline analysis. Dataset r/m Rho/theta European 2 (n = 37) 0.025 0.0037 European (n = 44) 0.034 0.0046 European and African (n = 51) 0.037 0.0056 Entire dataset (n = 70) 0.037 0.0044 9 Vol.:(0123456789) Scientific Reports | (2021) 11:18789 | https://doi.org/10.1038/s41598-021-98226-y www.nature.com/scientificreports/ effect in M. bovis genetic diversity when comparing with mutation (Table3). To make a broad comparison, the r/m parameter was estimated using a similar methodology for an MTBC dataset composed by 23 genomes, revealing a mean value of 0.48620, while for the 38M. bovis dataset of Patané and co-workers21 it evidenced a mean value of 0.98. In the first study there were only two M. bovis (M. bovis BCG and reference strain) within the 23 genomes included in the work, so the obtained value might be biased by the overrepresentation of M. tuberculosis genomes. In the second report, the M. bovis population under analysis was mainly recovered from American countries and livestock hosts. In contrast, in our dataset, a higher number of geographic locations and host species is represented, and genomes grouped into different clonal complexes with distinct population genetic signatures were also used, enabling a deeper and wider population knowledge. The differential r/m average values obtained with our dataset are consistent with the notion that the extent of recombination vary widely among lineages assigned to the same taxonomical species, so these results suggest that M. bovis clonal complexes might exhibit a differential impact of recombination, as also suggested by Didelot & Maiden72. Nevertheless, enlarging significantly this dataset with the inclusion of a higher number of M. bovis genomes would allow further clarification of this point. Both r/m and rho/theta parameters present variability among the tree branches, a result that is in agreement with reports concerning other bacterial species72,73. Finally, to confirm the recombination events identified by Gubbins pipeline, the different core multi-alignments were also independently tested in RDP4 software with six different algorithms. Globally, less than half of the events identified by Gubbins were confirmed by RDP4 (Tables4, 5). Considering the entire dataset, three recombination events were confirmed, two involving internal nodes and another one involving a single genome in a terminal branch and for which a clonal complex could not be assigned (Tables4, 5). The identification of events in terminal branches might be a sign that recombination is still ongoing in contemporary M. bovis strains or the result of misalignment70. In this putative recombination region, circa 20% of positions have an undefined nucleotide (N), which can therefore influence the recombination signal (Supplementary Fig.2). Moreover, this region affects the rrs gene, encoding the 16S ribosomal RNA that is expected to be highly conserved, so this putative recombination signal could be the result of a sequencing error or wrong alignment. Whole genome alignment between Mb0003 and M. bovis AF2122/97 was thus then performed and the presence of undefined nucleotides and of SNPs was confirmed, so the likely issues related to wrong alignment did not arrive as a consequence of the bioinformatics procedure implemented in this work. No gaps or undefined nucleotides were identified in the recombination regions of internal nodes (Figs.4, 5). With respect to these events, one encompasses exclusively Eu2 genomes, affecting the pks12 gene that encodes a probable polyketide synthase; while the other one is registered across Eu1 genomes and affects narX gene Table 4. Detailed information concerning the recombination events identified by Gubbins and RDP4 in the entire dataset. Genome positions according with M. bovis AF2122/97. Recombination event Identification Core-alignment positions Genome positions(a) Gene name Mb gene name Classification of gene function M. bovis isolate ID #1 Gubbins 945,923–945,950 1,220,297–1,220,324 PE PGRS22 Mb1121 PE-PGRS family protein Mb2026 #2 Gubbins; RDP4 1,176,674–1,177,221 1,475,305–1,475,975 rrs Mb5019 Ribosomal RNA 16S Mb0003 #3 Gubbins; RDP4 1,532,736–1,532,787 1,953,495–1,953,548 narX Mb1765c Probable nitrate reductase NarX Mb1792361 Mb7240415 #4 Gubbins 1,532,751–1,532,781 1,953,840–1,953,870 narX Mb1765c Probable nitrate reductase NarX Mb1792361 #5 Gubbins; RDP4 1,794,609–1,794,714 2,283,200–2,283,315 pks12 Mb2074c Probable polyketide synthase pks12 Mb0891 Mb1711 Mb1789 Mb1870 Mb1758 Mb2043 Mb1960 #6 Gubbins 1,794,627–1,794,780 2,283,713–2,285,136 pks12 Mb2074c Probable polyketide synthase pks12 Mb0003 #7 Gubbins 2,242,002–2,242,098 2,839,474–2,839,570 tatA Mb2121 Probable Sec-independent protein translocase membrane-bound protein tatA Mb0565 #8 Gubbins 3,244,551–3,244,556 4,003,420–4,003,425 espa Mb3646c Conserved hypothetical alanine and glycine rich protein Mb2043 Table 5. Statistical values associated with different algorithms implemented in RDP4 for the confirmed recombination events. Recombination event Alignment positions RDP (p-value) GENECONV (p-value) Bootscan (p-value) MaxChi (p-value) Chimaera (p-value) #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) Scientific Reports | (2021) 11:18789 | https://doi.org/10.1038/s41598-021-98226-y www.nature.com/scientificreports/ 80. Boritsch, E. C. et al. pks5-recombination-mediated surface remodelling in Mycobacterium tuberculosis emergence. Nat. Microbiol. 1, 15019 (2016). https:// doi. org/ 10. 1038/ nmicr obiol. 2015. 19. Acknowledgements This work was funded by Fundação para a Ciência e a Tecnologia, IP (FCT) / MCTES through national funds (PIDDAC) and co-funded by the European Regional Development Fund (FEDER) of the European Union, through the Lisbon Regional Operational Program and the Competitiveness and Internationalization Operational Program for Portugal 2020 or other programs that may succeed (project ‘Colossus: Control Of tubercuLOsiS at the wildlife/livestock interface uSing innovative natUre-based Solutions’, ref. PTDC/CVT-CVT/29783/2017, LISBOA-01-0145-FEDER-029783, POCI-01-0145-FEDER-029783). Strategic funding to cE3c and BioISI Research Units (UIDB/00329/2020 and UIDB/04046/2020) from FCT is acknowledged. ACR was supported by FCT through a doctoral grant (PD/BD/128031/2016). Author contributions M.V.C. conceived this work and secured resources and funding. A.C.R. performed the bioinformatic analyses under the guidance of M.V.C. and explored the data under MVC supervision. A.C.R. wrote the first draft of the manuscript and M.V.C. critically revised all versions. Both authors approved the final version. Competing interests The authors declare no competing interests. Additional information Supplementary Information The online version contains supplementary material available at https:// doi. org/ 10. 1038/ s4159802198226-y. Correspondence and requests for materials should be addressed to M.V.C. Reprints and permissions information is available at www.nature.com/reprints. Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/. © The Author(s) 2021