Introduction

The Mycobacterium tuberculosis complex (MTBC) is a group of closely related bacteria that cause tuberculosis (TB) in humans and other animals. Worldwide, TB is the 13th leading cause of death and is anticipated to once again become the top infectious killer following global efforts to combat the COVID-19 pandemic (WHO, 2023). The MTBC comprises both human-associated lineages (L1 to L9) and animal-associated lineages (La1, La2, La3, M. microti and M. pinnipedii) (Coscolla et al., 2021, Zwyer et al., 2021). The human-adapted MTBC exhibits a strong phylogeographical population structure, with some lineages occurring globally and others showing a strong geographical restriction (Gagneux, 2018). For example, L2 and L4 are widespread globally, with L2 dominating in Southeast Asia and L4 in Europe (Napier et al., 2020). L1 and L3 occur mainly in regions around the Indian Ocean (Menardo et al., 2021). L5 and L6 are highly restricted to West Africa (de Jong et al., 2010) and L7, L8 and L9 are exclusively found in East Africa (Blouin et al., 2012, Ngabonziza et al., 2020, Coscolla et al., 2021). Animal-associated lineages refer to lineages of the MTBC that primarily infect and cause disease in animals but can also infect humans in close contact with infected animals (Brites et al., 2018). Little is known about the geographic distribution and host range of animal-associated lineages (Zwyer et al., 2021).

Various factors such as human migration, local adaptation to host population, and host immune response are thought to have contributed to the current phylogeographical distribution of human-adapted MTBC members (Darwin and Stanley, 2022, Hershberg et al., 2008). Additionally, different lineages of MTBC are known to differ in virulence, metabolism and transmission potential, which can affect the spread of certain lineages and sub-lineages (Coscolla, 2017, Gagneux et al., 2006). Understanding the genetic diversity within MTBC is essential for gaining insights into how these bacteria adapt to different host populations and environments, and identify traits that may be important for survival, virulence, and antibiotic resistance.

MTBC is recognised for its clonal population structure and little genetic diversity between different strains, caused by the absence of horizontal gene transfer, plasmids, and recombination (Godfroid et al., 2018, Chiner-Oms et al., 2022). Single nucleotide polymorphisms (SNPs) and large sequence polymorphisms in the so-called regions of difference (RDs) represent the source of MTBC genetic diversity (Gagneux, 2018). RDs are known to differentiate different MTBC lineages from each other (Gagneux, 2018), but within lineage diversity is thought to be mostly driven by SNPs (Napier et al., 2020). While this has been shown to not be the case in some lineages (Sanoussi et al., 2021), these larger genomic structural differences have not been fully explored. Deletions in some RDs can surpass 10 kb in length and encompass groups of genes with diverse functions, which can be linked to factors like virulence and metabolic diversity (Bespiatykh et al., 2021, Bottai et al., 2020). Additionally, the IS6110 element has a notable impact on the genetic diversity of the MTBC. Transposition of IS6110 into various coding or regulatory regions, and genomic deletions through homologous recombination between proximal copies, can also translate into significant genotypic variation at strain level (Alonso et al., 2013, Antoine et al., 2021).

The pangenome represents the set of all genes present in a species. It comprises both the core genome, which consists of genes shared by all members of the species, and the accessory genome, which includes genes present in some but not all members (Brockhurst et al., 2019). Previous research has focused on pangenomic investigations of particular MTBC lineages (Baena et al., 2023, Ceres et al., 2022, Sanoussi et al., 2021). A recent study investigated the MTBC pangenome focusing on protein-coding sequences (Silva-Pereira et al., 2024). However, there is a gap in comprehensive pangenomic analyses encompassing both coding and non-coding regions. Despite the name, pangenome analyses tend to not focus on full genome comparisons, but rather on protein-coding sequences. This focus neglects non-coding RNAs and pseudogenes as well as essential elements in intergenic regions, such as promoters, terminators, and regulatory binding sites, despite their demonstrated role in selection and phenotypic impact (Tonkin-Hill et al., 2023). Moreover, the inference drawn from pangenomic analyses is often influenced by the use of short-read datasets. It has been shown that using short-read datasets can artificially inflate the size of the pangenome estimates, leading to incorrect assumptions about evolutionary scenarios (Marin, 2024). Here, we used a diverse collection of closed MTBC genomes including human and animal-associated lineages to develop our understanding of MTBC pangenome and evolution incorporating both coding and non-coding regions. The objectives of this study are: (i) to establish the true size of MTBC pangenome using a highly curated dataset, (ii) to assess specific biological functions that are prevalent in the core and accessory genome and (iii) to assess genetic variation between lineages and sub-lineages that may affect virulence, drug resistance, metabolism and evolutionary characteristics. We find that the MTBC has a small, closed pangenome and genome reduction is the driving evolutionary force that shapes the genomic evolution of the species. We identified distinct genomic features among lineages and sub-lineages, including sub-lineage specific RDs, that can account for variations in virulence, metabolism, and antibiotic resistance.

Methods

Summary of sequencing data and quality control

Previous studies have shown that closed genomes should be used for best estimation of pangenome size (Marin et al., 2024). From a combination of previously published research and publicly available data, 322 non-duplicated complete and closed (single contig) genome assemblies were selected for having (i) long-read sequencing (Nanopore or PacBio) or (ii) hybrid sequencing with Illumina across the MTBC lineages. Information related to the geographical location, source of isolation, and assembly method can be found in Supplementary File 1. The quality of all assemblies was assessed using BUSCO v5.4.7 (Manni et al., 2021) and only assemblies with ≥ 95% completeness were used for further analysis. The assemblies were set to start on dnaA gene using the fixstart command of Circlator v1.5.5 (Hunt et al., 2015). Complete circular assemblies were annotated using PGAP v6.4 with default options (Tatusova et al., 2016). This approach is recommended by previous analyses (Marin et al., 2024) as PGAP tends to annotate pseudogenes in cases of pseudogenisation and frameshift, rather than annotating them as new coding sequences. Visualisations of the PGAP annotation were made using the DNA Features Viewer Python library (Zulkower and Rosser, 2020). Lineage classification was performed using TB-profiler v4.4.2 (Phelan et al., 2019).

Nanopore sequencing and assembly

We sequenced lineages and sub-lineages that either had no sequencing data available in public databases or were under-represented in the collection. These isolates belonged to L5 (n=2, L5.1 and L5.2), L4.6 (n=3, L4.6.1, L4.6.1.1 and L4.6.1.2), L4.7 (n=3), L9 (n=1), La2 (M. caprae, n =1) and M. pinnipedii (n=1). A list of genomes sequenced in this study and their assembly information can be found in Table S1. Assembled genomes were deposited at the NCBI under the Bioproject PRJNA1085078.

DNA extraction suitable for PacBio SMRT sequencing was performed using the Genomic-tip 100/G (Qiagen Inc) as described previously (Ngabonziza et al., 2020). DNA was enriched in samples with low DNA concentration (1-60 ng/μl) using AMPureXP Beads (SPRIselect, Beckman Coulter) in a 1:1 ratio. Samples were incubated for 10 minutes at room temperature followed by pelleting the beads on a magnetic stand. The supernatant was discarded, and DNA was eluted from the bead using 10-20 μl of nuclease-free water (ThermoFisher Scientific). DNA concentration was quantified using a Qubit fluorometer and 1X dsDNA High Sensitivity (HS) assay kit (ThermoFisher Scientific). Nanopore sequencing was performed using the native barcoding kit 24 V14 (SQK-NBD114.24) according to the manufacturer’s instruction on the MinION Mk1C platform with R10.4.1 flow cells. Sequence adaptors were trimmed using Porechop v0.2.4 with “--end_threshold 95” and “-- middle_threshold 85” options (https://github.com/rrwick/Porechop.git). Sequences were assembled using flye v2.9.2 with default options (Kolmogorov et al., 2019) and iteratively polished using Racon v1.5.0 (Vaser et al., 2017) and the “consensus” subcommand of Medaka v1.8.0 (https://github.com/nanoporetech/medaka.git). Complete assemblies underwent Circlator, BUSCO and PGAP steps as above.

Pangenome analysis

Panaroo was used to characterise the MTBC protein-coding pangenome (Tonkin-Hill et al., 2020) using the “--merge_paralogs”, “--clean-mode strict” and “--core_threshold 0.95” options. Ceres et al. (2022) showed that genes impacted by pseudogenisation could emerge as accessory genes. To find split genes due to pseudogenisation, the reference pangenome created by Panaroo was compared to the H37Rv genome using BLASTn with a similarity e-value cutoff of 1×10⁻¹⁰ (Zhang et al., 2000). The BLASTn result was then filtered to only include queries with ≥90% identity and ≥75% of the H37Rv gene length covered using in-house Python scripts (https://github.com/conmeehan/pathophy). A core genome alignment was constructed using the built-in MAFFT aligner included in Panaroo (Katoh et al., 2002). The core phylogeny was created using RAxML-NG v1.2.0 with a GTR substitution model (Kozlov et al., 2019). Tree visualisations were created using the Treeio package v1.20.0 (Wang et al., 2020).

Pangenome graph construction

To look at total pangenome content (both coding and non-coding), multiple whole genomes were aligned into a graph using Pangraph v0.7.2 (Noll et al., 2023). The DNA sequences from gene blocks present in at least one sub-lineage but completely absent in others were extracted. These sequences were subjected to BLAST against the H37Rv genome, along with RvD1 and TbD1 genes to identify their corresponding Rv loci within the H37Rv genome and determine sub-lineage-specific regions of difference.

Functional genome analysis

The reference pangenome created by Panaroo was translated using transeq v5.0.0 (Rice et al., 2000) to obtain the protein sequences. COG functional categories were assigned using eggNOG Mapper v5.0.2 using default options (Huerta-Cepas et al., 2017). Genes that did not return a match within the eggNOG database were included in the not found category.

Functional domain annotation was performed with InterProScan v5.63-95 (Jones et al., 2014). For the gene association analysis, the accessory genome (genes present in less than 95% of genomes) was used as input for Coinfinder v1.2.1 (Whelan et al., 2020). STRING database (https://string-db.org/) was used to infer interaction networks of genes. Predicted functional partners with a high confidence interaction score (≥0.7) by STRING v12.0 were considered in this study.

Statistical analysis

Statistical analysis was conducted using R studio v2023.12.0+369 (R Core Team, 2023). Principal component analysis (PCA) based on gene presence and absence data from Panaroo and Pangraph derived accessory genome contents was performed using the prcomp package in R. The fviz_contrib function in Factoextra (v1.0.7) was used to identify the contribution of important variables from the PCA results (Kassambara and Mundt, 2020).

Results

Overview of the genomic dataset

We combined a total of 335 whole-genome sequences covering all known MTBC lineages. This included 324 genomes published previously, as well as 11 genomes from sub-lineages with few or no representatives in our collection to better cover the known phylogeny of MTBC (Supplementary Table S1). Except Antarctica, our collection includes genomes from every continent, and different host species including humans, cattle, voles and BCG vaccine strains (Figure 1-A). Africa represents the highest diversity of MTBC genomes among other continents, which is consistent with the fact that the continent is known to have the most diverse MTBC strains (Figure 1-B) (Gagneux, 2018). Most L1 genomes have been collected from East Africa, South and Southeast Asia, where this lineage is geographically common. L3 genomes were mostly isolated from India, while L2 and L4 genomes have been collected from different continents, representing the global distribution of these lineages. L5 and L6 genomes were exclusively found in West Africa, where they contribute significantly to the overall burden of TB in this region (Coscolla et al., 2021).

Figure 1:

Pangenome characteristics and functional diversity

Panaroo estimated a protein-coding pangenome consisting of 4,252 genes, including 3,639 core genes (present in more than 99% of genomes) and 613 accessory genes including 258 soft core, 218 shell and 137 cloud genes. After mapping all the genes in the pangenome to H37Rv, we identified 127 genes that were misclassified as accessory genes, including 25 out of 258 soft core, 24 out of 218 shell and 78 out of 137 cloud genes (Supplementary Table S2). This approach increased the core genome size from 3,639 to 3,746 and reduced the accessory genome from 613 to 506. This shift occurred because 107 genes out of 127 genes detected as false accessories were re-classified as part of the core genome. To examine genomic diversity among lineages, we investigated variations in the size of the accessory genome among different lineages. This analysis revealed significant differences in accessory genome size among different lineages (Kruskal-Wallis test, p= 2.2e-16). Post-hoc Dunn tests with Bonferroni correction identified significant pairwise differences between lineages (Supplementary Table S3). L6, La1 and M. microti have significantly smaller accessory genomes compared to others (Figure 2, Supplementary Table S3). L2 has a significantly smaller accessory genome in comparison with L1, L3, and L4.

Figure 2:

Functional annotation of core and accessory genes showed that genes of unknown function were the largest functional group in both (26.8 % and 12.0% respectively). The core genome of MTBC encodes a high proportion of genes involved in lipid transport and metabolism (category I) and cellular transport and secretion systems (category N), indicating the importance of these biological functions in defining the species. Functional breakdown of the accessory genome also revealed variations in genes involved in DNA replication, repair and mobile elements (category L) and lipid transport and metabolism (category I) across the species. A breakdown of the functional composition of the core and accessory genome is shown in Figure 3.

Figure 3:

Population structure analysis

The population structure was assessed using SNPs in the core genome and PCA plots were constructed based on presence-absence patterns in the accessory genome determined with Panaroo and Pangraph (Figure 4). In the PCAs, genomes group together based on their lineage, indicating variations in accessory genome are driven mainly by lineage structure (Figure 4-B and 4-C). The first principal component (PC1) highlights the contrast between L2 from other members of MTBC. In both PCA plots, L2 formed a distinct group, defined primarily by RD105 and RD207 deletions. RD105 serves as a classic marker for all L2 genomes and RD207 deletion is found in all L2.2 sub-lineages; the single L2.1 genome is separated from the L2.2 cluster due to the absence of RD207. The second principal component (PC2) demonstrates the difference between modern human-associated lineages (L2, L3 and L4) from other lineages. Important variables driving this separation are the TbD1 gene deletion in modern lineages, and gene deletions associated with RD11 (prophage phiRv2) in L6, L7, La1 and La3 genomes and some L1 sub-lineages. L1 genomes exhibit considerable genetic diversity primarily due to RD3 (related to prophage phiRv1) and RD11. The latter was identified in all L1.1 sub-lineages and certain L1.2 sub-lineages (L1.2.1, L1.2.2, and L1.2.2.2). The positioning of the single L8 genome, marked with a positive sign at the centre of PCAs, implies that L8, serving as the ancestor of the MTBC, exhibits similarities in the accessory genome with other lineages.

Figure 4:

The accessory genome of MTBC

Among the 506 filtered accessory genes from Panaroo output, 107 were identified as being specific to certain lineages and sub-lineages (Figure 5-A, Supplementary Table S4). Similarly, we identified 82 lineage and sub-lineage-specific genomic regions ranging from 270 bp to 9.8kb, encompassing a total of 116 genes and 4 non-coding RNA using the Pangraph approach (Figure 5-B, Supplementary Table S4 and Supplementary File S1). Our analysis revealed both previously known RDs and new deletions in the accessory genome, suggesting that it originates from genome deletions. Lineage and sub-lineage-specific genomic regions were consistently identified using both methods, including TbD1, RD9, RD10, RD105, RD303, RD312, RD316, RD702, RD711, RD713, RD750, ND1, RDoryx, and N-RD25bov/cap (Figure 6).

Figure 5:

Figure 6:

Panaroo failed to identify partial gene deletions caused by RD9, RD10, RD13, RD105, RD207, and RD702. For example, RD702 is present in all L6 and L9 genomes and covers a 641 bp region (Figure 7). RD702 affects less than a quarter of the bglS gene (470 bp from 2,076 bp) and mymT gene (126 bp from 162 bp gene) in the 5’ to 3’ direction as well as a regulatory RNA (ncRv0186c) in the 3’ to 5’ direction. Panaroo did not identify the 470 bp deletion downstream of the bglS and this gene was classified as being in the core genome. These results could be explained by the fact that Panaroo uses a relaxed alignment threshold along with neighbourhood gene information to cluster gene families (Tonkin-Hill et al., 2020). This can lead to the grouping of partially deleted genes with those that are complete. The Pangraph alignment approach identified partial gene deletion and non-coding regions of the DNA that were impacted by genomic deletion. However, the graph approach failed to identify four structural variants, including a gene deletion corresponding to RD8 (Rv3618) in L6 and animal-associated genomes. Despite these limitations, both methods yielded consistent results indicating that the accessory genome of MTBC is small and is acquired vertically from a common ancestor within the lineage. Gene association analysis using coinfinder confirmed that most accessory genes (as determined by Panaroo) appear in pairs or groups in larger multi-gene RDs (as determined by Pangraph) (Supplementary Figure S1). Thus, the small accessory genome is contained within an even smaller number of RDs. Coincident gene blocks that were associated with lineages were RD8, RD9, RD10, and RD178, whereas those that were distributed across the phylogeny belonged to RD3 and RD11.

Figure 7:

Discussion

MTBC accessory genome is small and is acquired vertically

The MTBC pangenome size has been reported to vary across studies. In a recent investigation conducted by Marin (2024), the accessory genome ranged from 314 to 2951 genes, as a function of the assembly quality, software and specific parameters used. Genome sequencing using short-read technology has some limitations that can result in incomplete assemblies, particularly in regions with repetitive sequences that are difficult to sequence (Meehan et al., 2019). About 10% of the coding sequences of the MTBC genome are composed of repetitive polymorphic sequences known as the PE (Proline-Glutamic acid) and PPE (Proline-Proline-Glutamic acid) families, implicated in host-pathogen interactions and virulence (Cole et al., 1998). Incomplete assemblies can impact the accurate detection of genes and may lead to artificial inflation of the estimated pangenome size (Ceres et al., 2022, Marin et al., 2024). We used a curated dataset to address these issues and estimate the true size of the MTBC pangenome. We showed that MTBC has a compact accessory genome that is acquired vertically. This finding is consistent with recent research, which found a small accessory genome with a clonal genetic structure in MTBC (Baena et al., 2023, Ceres et al., 2022, Silva-Pereira et al., 2024). Using a BLASTn analysis approach, we identified core genes that were designated as accessory and further reduced the accessory genome from 613 to 506. For example, indel mutations observed in the lppA gene across 128 genomes resulted in the misclassification of this gene as an accessory (labelled as group_1700). Notably, over 50% of the cloud genes (genes present in less than 15% of strains) were found to be core genes. These genes mapped to the reference genes using the BLAST approach and appeared to be fragmented genes caused by misassemblies.

Attributes of the accessory genome

We employed a pangenomic approach to analyse the genetic diversity within the MTBC, aiming to understand its evolution and functional potential. In addition to inferring a gene-based pangenome with Panaroo, we also analysed variation in non-coding regions by inferring a pangenome graph with Pangraph. Examining non-coding regions is particularly important for species with limited genetic diversity, as they play a vital role in controlling gene expression, stress response and host adaptation (Arnvig et al., 2014). Our study demonstrated that, although both methods produced reliable outcomes, the whole-genome alignment approach could capture higher genetic variations, especially those involving large structural variants (>200 bp). Moreover, this method identified non-coding regions of the genome that were affected by genomic deletions. Results from both methodologies revealed that the accessory genome of MTBC is a product of gene deletions, that can be classified into lineage-specific deletions and lineage-independent deletions, with the former category being the more common. Coinfinder results also showed that coincident genes in the accessory genome either belonged to particular lineages or were dispersed throughout the phylogeny. Region of difference, which are large deletions often encompassing multiple genes, can act as distinctive markers not only for lineages, as demonstrated previously, but also for sub-lineages. These warrant further investigation to understand the biological consequences of such deletions on traits like virulence, drug resistance, and metabolism, along with SNP only analyses. On the other hand, we observed several genomic deletions that appear to be lineage-independent and may indicate gene instability (Brosch et al., 2002). This includes RD3 and RD11 corresponding to the two prophages phiRv1 and phiRv2. One drawback of our study was that many of the lineages were poorly represented, especially in the L5, L6, L9 and all animal-associated lineages. The production of more closed genomes from the geographic regions where these are abundant (e.g. East and West Africa) will provide further information about the RD distributions within the MTBC and better understanding of the driving evolutionary forces.

Differential traits between lineages and sub-lineages

The size of the accessory genome content varied significantly among different lineages. L2 genomes exhibited a high degree of clonality and a smaller accessory genome. L2 strains are characterised by a large number of IS6110 copies, possibly due to the RD207, which affects the CRISPR-Cas type III-A system, responsible for adaptive immunity against mobile genetic elements (Alonso et al., 2013). While RD207 encompasses seven genes (Rv2814c-Rv2820c) (Bespiatykh et al., 2021), our analysis revealed that only two genes, Rv2819c and Rv2820c, are specifically absent in L2 and can serve as a genetic marker for this lineage. Other RDs affecting all L2 genomes are RD3 (9.2 kb), and RD105 (3.4 kb). The latter is exclusively present in all L2 genomes and is caused by the insertion of IS6110 (Alonso et al., 2013). Previous research proposed a potential role of RD105 in conferring drug resistance in L2.2 strains (Qin et al., 2019). RD207 (7.4 kb), linked to the insertion of the IS6110 element, is also present in all L2 genomes except for L2.1. L3 genomes exhibited significantly less genomic deletion in comparison with L2, and RDs affecting all L3 genomes were RD316 (1.2 kb) and RD750 (765 bp). The genes affected by these RDs are involved in lipid metabolism (category I) and secondary metabolite biosynthesis and transport respectively (category Q), which may contribute to metabolic variations in L3 strains. L4 exhibits substantial genetic diversity primarily due to variations in the presence and absence of specific RDs at the sub-lineage level. For example, RD145 (4.7 kb) and RD182 (6.5 kb) were detected in L4.1.2.1 genomes isolated from Africa, Asia, Europe and South America. However, none of the affected genes appears advantageous or essential for MTBC growth and survival and this sub-lineage is successful worldwide (https://tbdr.lshtm.ac.uk/sra/browse). Other L4 sub-lineage-specific RDs detected in this study with no known advantage are RD115 in L4.3.4, RD727 in L4.6.1, RD219 in L4.8, RvD1 in L4.8 and L4.9. RD4 in La1, L5.3 Del in L5.3, and RD750 in L3, RD147c in L1.

Nonetheless, some genomic deletions can affect the growth advantage gene group in various lineages. Disruption of these genes has been shown to provide a fitness advantage for the in vitro growth of H37Rv by analysis of saturated Himar 1 transposon libraries (DeJesus et al., 2017). For example, RD174 specific to L4.3.4 intersects with RD743 specific to L5, and both deletions affect two genes (Rv1995 and Rv1996) including cmtR transcriptional regulator. The latter regulates bacterial stress response to metal deficiency and is important during persistence in phagocyte (Flores-Valdez et al., 2020). While the absence of these genes may provide a fitness advantage for in vitro growth, the impact of this deletion is likely to be different in in vivo conditions. The low prevalence of L4.3.4 strains in the TB Profiler database may suggest either the lower fitness of this sub-lineage, or its presence in under-sampled geographic areas. RDs impacting the growth advantage gene group were prevalent in L5 genomes (L5.3 Del, RD743 and RD711). For instance, RD711 affects the Rv1335 gene, which is annotated as a sulphur carrier protein involved in cysteine biosynthesis by PGAP. This may imply that the genome reduction observed in L5 has evolved to depend on host factors for essential nutrients, consequently making these strains difficult to grow in a culture medium (Sanoussi et al., 2017).

L6 and animal-associated lineages accumulated a higher number of deletions including RD7, RD8, RD9 and RD10. Comparative analysis of the accessory genome size in L6, La1, and M. microti revealed significant differences compared to other lineages. The identified genomic regions in L6, L9, and animal-associated lineages impact genes associated with virulence, potentially explaining their reduced virulence in human hosts. RD7 contains the mammalian cell entry (mce3) operon implicated in the invasion of host cells, while RD9 includes genes related to biliverdin reduction, potentially aiding in protection against oxidative stress inside macrophages (Szklarczyk et al., 2023). RD10 impacts genes related to lipid storage crucial for MTBC transition to a dormant state (Deb et al., 2009). RD1BCG and RD1mic affect genes encoding the ESX-1 secretion system, crucial for virulence and phagosomal escape, with RD1BCG being deleted during BCG vaccine development (Szklarczyk et al., 2023). Overall, the findings of this study suggest that, despite having a compact accessory genome, it can account for certain biological distinctions among MTBC. Previous research has highlighted the significance of random genetic drift in shaping the evolutionary dynamics of the species (Hershberg et al., 2008). Gene loss, driven by genetic drift, is likely to be a key contributor to the observed genetic diversity within the MTBC.

A need for more high-quality genomes

While this dataset represents high-quality MTBC genome sequences that cover most of its known diversity, it suffers from sampling bias as some lineages are over-represented whereas others have a few representatives. The over-representation of L2 and L4 in the dataset is primarily attributed to their association with Bioprojects analysing TB outbreaks or originating from countries with extensive sampling, highlighting a bias in sequencing efforts towards more prevalent or historically studied lineages. Many sub-Saharan African countries do not have a representative sample in our dataset, and this may explain the under-representation of L5 and L6 in the dataset. These lineages contribute significantly to the overall burden of TB in this region (Coscolla et al., 2021). Additionally, we found that certain RDs that were previously assumed to be lineage-specific were identified in other lineages (Supplementary table S4). The observed variation in RD distribution emphasizes the need for an expanded genomic dataset to provide a more comprehensive understanding of their prevalence, specificity, and implications in the evolution of the MTBC. Future sequencing efforts should target high-burden TB countries, ensuring a comprehensive representation of the population structure to enhance our understanding of MTBC diversity and evolution.

Conclusion

In conclusion, this study reveals several significant findings. Firstly, the MTBC has a small, closed pangenome, predominantly driven by genome reduction. This observation underlines the importance of genomic deletions in the evolution of the MTBC. The study also highlights distinct genomic features across MTBC, including sublineage-specific RDs, which may play a crucial role in variations in virulence, metabolism, and antibiotic resistance. The analysis further reveals that the accessory genome of MTBC is a product of gene deletions, which can be classified into lineage-specific and independent deletions. Finally, the study underscores the impact of non-coding regions and large structural variations on MTBC’s genetic diversity, offering a new perspective on MTBC pangenomics.

Acknowledgements

We would like to thank all the staff at the BCCM/ITM collection for the growth of strains for sequencing in this study. We thank Alasdair Hubbard for assistance with Nanopore sequencing.

C.J.M. and M.B. are supported by the Academy of Medical Sciences (AMS), the Wellcome Trust, the Government Department of Business, Energy and Industrial Strategy (BEIS), the British Heart Foundation and Diabetes UK and the Global Challenges Research Fund (GCRF) via a Springboard grant [SBF006\1090].

Author contributions

M.B. contributed to conceptualization, methodology, formal analysis, investigation, data curation, writing-original draft, writing – reviewing & editing, and visualization.

M.M. contributed to methodology, formal analysis, resources, data curation and writing – reviewing & editing.

M.F. contributed to methodology, data curation, writing – reviewing & editing and supervision.

J.C.T contributed to methodology, resources, and writing – reviewing & editing.

M.R.D.S. contributed to methodology, validation, formal analysis, investigation, writing – reviewing & editing and visualisation.

C.J.M. contributed to conceptualization, methodology, formal analysis, investigation, writing-original draft, writing – reviewing & editing, supervision, project administration and funding acquisition.