Global epidemiology of 
Tuberculosis (TB) caused by Mycobacterium tuberculosis ( Mtb ) is a significant global health concern, with an estimated 10.6 million infections and 1.6 million casualties by 2021 . The COVID-19 pandemic has disrupted access to medical healthcare, including TB diagnosis and treatment programs, aggravating the TB burden, and compromising the progress in TB control achieved in recent decades , . Furthermore, the incidence of TB has been increasing in several countries, including Ecuador, where the rate reached 48 cases per 100,000 inhabitants in 2021 due, among other factors, to COVID-19 containment measures , . The increase in global migration has contributed to the spread of TB, mainly in high-income countries where migrants seek better economic, educational, or living opportunities. In the Ecuadorian context, migration is primarily directed toward the USA, Spain, Italy, Canada, and Chile – . Identifying the transmission routes of TB cases is essential for reducing potential transmission networks. Different public health organizations have implemented screening programs among household contacts to control the spread of the disease, including the TB-Directly Observed Treatment Short Course (TB‐DOTS) . Genomic approaches, such as whole-genome sequencing (WGS) of Mtb strains, have been a landmark in the traceability of transmitted TB cases and have provided invaluable information on drug resistance and sub-lineage patterns – . WGS and Bayesian phylogenetic approaches have reconstructed the historical patterns of TB spread in Central and South America, dating back to the introduction of Mtb strains in these regions. The apparent ancestral emergence and marked diversification of the L2 (ancestral and modern Beijing sub-lineages) and L4 lineages have also been emphasized, reflecting their close correlation with restricted geographic distribution, leading to the independent emergence of multiple sub-lineages and local adaptation to distinct human populations – . Studies conducted in Brazil, Paraguay, Mexico, and the United States have applied genomic approaches to identify frequent transmission clusters between prison inmates, drug users, migrants, and mixed groups – , highlighting the need to prioritize contact tracing to groups with a higher likelihood of retrospective clustering to improve TB control. The Mycobacterium tuberculosis complex (MTBC) comprises distinct phylogenetic lineages that have evolved over centuries – . The Euro-American lineage (L4) exhibits extensive sub-lineage diversity both within and across countries. In Central and South America, the Caribbean, Europe, and Middle Africa, the most prevalent L4 sub-lineages include Latin American (LAM), Haarlem, X-type, and T families – . Previous reports have revealed that the LAM and Haarlem families are the most prevalent in Ecuador, with few cases in the Beijing family – . Transmission clustering in these cases has often been linked to isolates from neighboring countries , . A pivotal WGS study on a limited dataset of Ecuadorian Mtb isolates identified the 4.3.2/3 (LAM) and 4.1.2 (Haarlem) sub-lineages, highlighting the significant genetic diversity present . This study underscores the importance of understanding the local transmission dynamics of TB within Ecuador, particularly considering the potential influence of events occurring outside national boundaries, associated with migration patterns, trade, and regional cooperation. Although studies have investigated transnational TB transmission between Ecuador and countries such as Colombia , evidence remains sparse. We analyzed the transmission networks of 88 Mtb isolates of L4 lineage collected from various cities across Ecuador in comparison with 415 publicly available genomes in 19 Latin American countries. To accomplish this, we used a combination of phylogenetics and composition analysis. Additionally, we explored the role of genetic variation in virulence genes in the recent transmission of TB in Ecuador. Our findings provide crucial insights into the genetic diversity of Mtb in Ecuador, highlighting how migration influences the TB burden and transmission dynamics locally. This study contributes to the development of more effective targeted TB control strategies tailored to the unique characteristics of Mtb populations. 

Genomic and functional analysis of Ecuadorian Mtb isolates This study focused on genomic analysis of raw reads of 88 M . tuberculosis isolates collected in Ecuador . To extend our investigation of genomic variability and epidemiology, we combined these data with a selection of 415 M . tuberculosis isolates from different countries (see Methods). Initially, we used the raw sequences (average coverage of 61X) of the 88 Ecuadorian isolates to perform quality control and genome assembly for each isolate using Unicycler and Pilon while discarding contigs smaller than 300 bp. The quality of the genomes was assessed and an average N50 of 65,747. Annotation was performed using Prokka software. These analyses yielded approximately 4298 genes per isolate, including coding sequences and RNA gene transfer. Descriptive statistics for this analysis are presented in Table and Supplementary Table . Genomic analysis of the Ecuadorian M. tuberculosis isolates revealed that they belonged to the Euro-American Lineage (lineage 4). Among the sublineages, 4.3.3 (27.3%, 24/88) and 4.1.1 (23.9%, 21/88) were the most prevalent, followed by 4.4.1.1 and 4.1.2.1 (11.4%, 10/88 each). Clades 4.1.2, and 4.3.4.1/2 were present, but in smaller proportions. Functional annotation and metabolic insights We performed functional annotation of each of the 88 Ecuadorian genomes. A substantial proportion (77%) of the annotated genes was associated with specific functional roles (e.g., metabolic processes, cellular processes, energy, protein processing, and stress response, Defense, and Virulence). We enriched our dataset with functional annotations, assigning Enzyme Commission (EC) numbers to 1,061 proteins, Gene Ontology (GO) classifications to 918 proteins, linking 814 proteins to specific KEGG pathways, and identifying approximately 816 virulence factors distributed across three databases. This enrichment analysis provided deep insights into the metabolic functions and potential pathogenic mechanisms of the isolates (Supplementary Table ). Our subsystem analysis, categorized according to the KEGG database, showed that approximately 1,961 genes per isolate are involved in various biological processes and structural complexes (Supplementary Table ). The analysis revealed that the majority of the annotated protein-coding genes were related to metabolism (40.3%), with 37.9% of these genes being associated with cofactors, vitamins, and prosthetic groups (300 genes). This was followed by genes related to the stress response, Defense, Virulence (10.8%), and energy (10.5%). In our detailed metabolic pathway analysis, we observed a predominance of genes involved in amino acid (369 genes), carbohydrate (312 genes), lipid (262 genes), and xenobiotic biodegradation and metabolism (240 genes). These pathways are crucial for the pathogenicity and survival of Mtb (Supplementary Table – ). Notably, variability was observed among different sub-lineages. Sub-lineage 4.1.2.1 shows a higher proportion of proteins related to metabolism and a lower proportion of proteins related to protein processing. Conversely, sub-lineage 4.1.2 presents a higher proportion of proteins related to protein processing and fewer proteins related to metabolism (Supplementary Fig. A). Lineage 4.1.2.1 exhibits a higher number of gene variations related to metabolism, suggesting significant metabolic diversity, indicating a complex and adaptable metabolic network, whereas lineage 4.3.2/3 showed greater variations in genes associated with protein processing, suggesting a diverse set of genes involved in protein synthesis, folding, and degradation, reflecting the complexity of protein processing mechanisms (Table ). Pangenome analysis of the Ecuadorian Mtb isolates We also performed pan-genome analysis to understand the genetic diversity and conservation of 88 Mtb isolates from Ecuador. Pangenome reconstruction of the 88 Mtb isolates revealed 4,397 gene families, including 3,104 classified as core (present in at least 99% of isolates), 666 as accessory, and 270 as cloud gene families (Supplementary Fig. B and Supplementary Table ). Interestingly, 70.5% of the isolates belonged to core gene families, suggesting minimal within-genome variability and emphasizing the high level of genetic conservation among the isolates. According to pangenome calculations, a b value of 0.086 in the power-law regression model indicated a close pangenome for Mtb . Three isolates (S1454, S1477, and S1453) from sub-lineages 4.1.1 and 4.3.2/3 showed the highest number of cloud genes (261, 200, and 177, respectively). Most gene families within the core and accessory partitions belong to metabolic subsystems. Conversely, the majority of the unique gene families were classified within the Environmental Information Processing subsystem (Supplementary Fig. C,D). Genotypic drug-resistance analysis Furthermore, when analyzing the phylogenetic relationships and genetic diversity associated with resistance in Ecuadorian isolates (Fig. ), our findings revealed that sub-lineages 4.3.2/3 and 4.1.1, exhibited a higher prevalence of resistant variants (46.6% and 26.1%, respectively). We identified 42 single-nucleotide variants (SNVs) in resistance-related genes, with rpoB showing the most mutations. Additionally, SNVs were more frequently identified in genes related to intermediary metabolism and respiration (21.71%), followed by cell wall processes (21.34%) and conserved hypotheticals (20.44%). Exhaustive details of the SNV are provided in Supplementary Table . To gain a deeper insight into the genomic epidemiology of TB in Ecuador, we incorporated an additional 415 Mtb samples from 19 countries categorized as both continental and transoceanic migratory nations (See “ ” section). All these 503 Mtb genomes belonged to the L4 lineage, including 4.3.2/3 (35.4%), 4.1.2 (22.7%), 4.4.1 (12.7%), and 4.1.1 (10.7%) sub-lineages (Supplementary Table ). Moreover, 63.8% (321/503) of Mtb samples were genotypically susceptible to all anti-TB drugs. The remaining 36.2% (180/503) were resistant to at least one antibiotic. The clinical classification of these drug-resistant samples showed that 16.7% were identified as MDR-TB, 8.5% as HR-TB, 4.2% as pre-XDR-TB, and 1.6% as RR-TB. Table summarizes the drug-resistant clinical classification for the 503 Mtb samples per region and Table provides an overview of their drug-resistant canonical mutations. An extended list of canonical mutations is provided in Supplementary Table . Since multidrug-resistant TB (MDR-TB) and extensively drug-resistant TB (XDR-TB) can increase domestic transmission , we determined the frequency of resistance-associated canonical mutations in 503 Mtb isolates. The most frequent mutations were katG Ser315Thr (n = 92), rpoB Ser450Leu (n = 68), rpsL Lys43Arg (n = 21), embB Met306Ile (n = 17), pncA Gln10Pro (n = 14), and gyrA Ala90Val (n = 11), which confer resistance to INH, RIF, STR, EMB, PZA, and FQs, respectively. Other common canonical mutations included fabG 1 − 15C > T (n = 16), ahpC -74G > A (n = 15), and gid 329_330delTG (n = 16), which are associated with INH and STR resistance, respectively. TB population structure and transmission clusters To better understand the genetic diversity and potential transmission dynamics of TB in Ecuador, we performed a population structure analysis of 503 Mtb samples (Fig. ). Clusterization based on the rhierBAPS approach revealed four primary clusters at the initial hierarchical level, which were further divided into 15 sub-clusters at the second hierarchical level. Phylogenetic reconstruction of 503 isolates using 1,437 SNPs revealed four main clusters with different branch colors. Additionally, the sub-clusters (level 2), lineages, geographical distribution, and resistance are represented by colored concentric lines. Cluster 1 (pink) encompassed 128 isolates divided into three sub-clusters (1–3), containing samples belonging to 4, 4.4.1.1, 4.2.2, 4.6.2, and 4.7/8 sub-lineages. Cluster 2 (in orange) comprises 125 isolates, including three sub-clusters (4–6), with the majority belonging to the 4.1.2.1 sub-lineage, and some isolates belonging to the 4.1.2 sub-lineage (sub-cluster 6). Cluster 3 (in green) was composed of 178 isolates, distributed across three sub-clusters (7–9), primarily associated with 4.3.2, 4.3.3, 4.3.4.1, and 4.3.4.2, and, unexpectedly, a 4.1.2 isolate within sub-cluster 9. Finally, Cluster 4 (in blue) encompasses six sub-clusters (10–15), including 71 isolates predominantly from the 4.1.1, 4.1.1.1, and 4.1.1.3 sub-lineages. Notably, the isolates formed subclusters 13 and 15, with a small portion of subcluster 10 corresponding to the 4.1 sub-lineage. In addition, we analyzed the isolates for transmission clusters by pairwise comparisons using a minimum distance of 12 SNPs with MTBseq. Thus, strains in the same transmission genomic cluster (TGC) had fewer than 12 SNPs. Using this strategy, 92.8% (n = 467) of the isolates were classified into 51 TGCs, ranging from two to 63 members, whereas 35 isolates were not classified. Based on the number of members per TGC (see “ ” section), we classified seven members as small , 20 as medium , and 24 as large (Supplementary Table ). Most of the 51 TGCs resolved using this methodology were within Cluster 3 (37.3%) and Cluster 1 (29.4%), as determined using rhierBAPS (Fig. A). In particular, Mtb samples identified as susceptible were more prevalent in larger TGC than in resistant samples. Most pre-XDR isolates exhibited clonal relationships either among themselves or with isolates from other countries. Of the 35 isolates not grouped by MTBseq, 28.6% belonged to sub-lineage 4.3.2/3/4 and 25.7% to sub-lineage 4.8 (Fig. A and Table ). Upon closer examination of the five largest TGCs, specifically, TGC_1, TGC_2, TGC_9, TGC_18, and TGC_19, distinct geographical profiles were identified. Notably, TGC_1 (n = 63 isolates) predominantly comprised Mtb isolates from Paraguay (46.3%), Peru (19.0%), and Spain (15.9%). TGC_18 (n = 56) included isolates from Ecuador (41.0%) and Paraguay (26.8%), whereas TGC_19 (n = 31) consisted solely of isolates from Spain (38.7%). TGC_2 (n = 30) was predominantly isolated from Argentina (53.3%), whereas TGC_9 (n = 24) was isolated mainly from Canada (54.1%). Among the 88 Ecuadorian Mtb isolates, 84 were grouped into 15 of 51 identified TGCs. These included five medium-sized and ten large TGCs. Notably, ten of these TGCs (TGC_2, TGC_3, TGC_5, TGC_14, TGC_15, TGC_18, TGC_20, TGC_21, TGC_37, and TGC_47) included isolates from Ecuador alongside isolates from various countries, unveiling previously unrecognized potential continental and intercontinental connections. In particular, TGC_37 and TGC_14 exhibited a close genetic relationship among their members, showing distances of ≤ 9 SNPs and ≤ 6 SNPs, respectively. In addition, TGC_37 comprised five isolates from Ecuador and one from an Ecuadorian migrant in Spain, highlighting the impact of migration on TB spread. Conversely, TGC_14 grouped two isolates from Ecuador with five from Brazil, indicating active transmission between these nations (Supplementary Table ). The Ecuadorian isolates within each TGC exhibited a high degree of clonality. For example, approximately 90% of the Ecuadorian isolates in TGC_11 and 70% in TGC_15 and TGC_18 had a genetic distance of zero SNPs (Supplementary Tables , and Supplementary Figs. , ). Similarly, in some TGCs, many isolates remained distinct and did not cluster with other isolates (singletons) in various proportions. (Supplementary Table ). Regarding the sub-lineage composition, the most representative sub-lineages in the Ecuadorian isolates among the TGCs were 4.3.3 (46.9%, 23/49; TGC_18), 4.1.1 (32.7%, 16/49; TGC_15), and 4.4.1.1 (20.4%, 10/49; TGC_11). Concerning the potential transmission of isolates, we identified four pre-XDR TB cases in TGC_15 and three pre-XDR isolates in TGC_15 and TGC_18, corresponding to clonal samples (genetic distance < 1 SNP). Among these three TGCs, 41% were from previously treated patients with TB, 34% were from untreated patients, and 25% were from patients currently undergoing treatment. TB transmission networks analysis TB transmission networks, that is, the suggested transmission routes of TB considering the genetic distances between samples, their geographical location, and year of isolation, were inferred for the 503 Mtb isolates using TransFlow. Using this strategy, 46 transmission networks were identified; however, only 18 of these networks had more than three members, as shown in Fig. A. As can be seen, most of these transmission networks (with more than three members) correspond to different TGCs, meaning they present more than 12 SNPs between them. However, six transmission networks (7, 11–14, and 16) corresponded to TGC_18, indicating the low genetic distance between them (Fig. A). Considering the three clustering strategies employed (rhierBAPS, MTBseq, and TransFlow), we identified 17 convergent transmission groups to define potential TB transmission networks between Ecuadorian and non-Ecuadorian isolates. The geographic distribution of networks involving Ecuadorian isolates is shown in Fig. B, where TGCs that were confirmed and congruent with TransFlow are highlighted in red. Notably, some of these Ecuadorian samples (collected between 2019 and 2021) showed clonality with isolates from Colombia (n = 2) and Latin American migrants in Spain (n = 4), which were detected in 2014 and 2015, respectively (Fig. A, Supplementary Figs. , and Supplementary Table ). We identified several transmission networks showing connections among countries, including Paraguay, Argentina, and Brazil, with specific samples acting as potential index cases in the corresponding transmission networks (Supplementary Fig. ). The distribution of local and migrant Ecuadorian isolates revealed diverse representations spread across 12 TGCs (Fig. C). Some of these isolates showed clonality with samples from Colombia, Latin America, and Ecuadorian migrants in Spain, or were identified as potential index cases within transmission networks (Fig. C and Supplementary Fig. ). Notably, TGC18 and 19 exhibit diverse nationalities, with 10 and 8 different origins represented within these clusters, respectively. The short genetic distance between different TGCs indicated a close relationship, suggesting a possible joint origin. This relationship was supported by the genetic distances observed between the Ecuadorian isolates, implying a shared ancestry among the TGCs (Fig. D). These findings highlight the complex transmission dynamics of TB, with evidence of cross-border and international transmission networks involving Ecuadorian populations and those from other Latin American countries. The diversity of isolates and potential index cases underscore the need for targeted public health interventions to address the multifaceted nature of TB epidemiology. Genetic diversity within the virulence-associated genes To characterize mutations in the genes involved in host adaptability, we analyzed SNPs in genes commonly associated with MTBC virulence in Ecuadorian isolates . A total of 303 SNPs were identified in these virulence genes, and all 88 isolates had at least one SNP each in the mce1F, mmpL4, phoR, ctpV, pepD, mce3F, fadD13, and nuoG genes. The pks12, fadD5, mce3C, pks12, nuoG, and katG genes were mutated with at least one SNP in more than 50% of the isolates. The top six genes with the most significant SNPs were plcA, plcB, pks7, pks12, phoR, and PPE46 (Supplementary Tables , ). Among the isolates, 46.6% (41/88) had more than 40 SNPs, with a maximum of 79 SNPs identified within virulence genes. Isolates corresponding to sub-lineages 4.1.1 (TGC_15 and TGC_37), 4.4.1.1 (TGC_11), and 4.3.2/3/4 (TGC_18, _20, _21, _26, _29, and _30) presented higher numbers of virulence-associated polymorphisms. Additionally, two members of TGC_14 (S0017 and S0039) and one ungrouped sample (S2193) showed a high number of mutations in the virulence genes (70, 71, and 74 SNPs, respectively). Interestingly, the five isolates associated with sub-lineage 4.7/8 showed the lowest number of SNPs in the virulence genes (15–25 SNPs). Phylodynamics of TB in Ecuador The divergence time to the most recent common ancestor (TMRCA) estimated for the Ecuadorian isolates fell within a density interval between 697 and 1,475 years, with an estimated time of 1,054 years. We also estimated TMRCA for different TGCs, identifying TGC_15 as the oldest (218 years before present, YBP) and TGC_03 as the youngest (130 YBP). Notably, isolates corresponding to TGC_02, TGC_21, and TGC_47, each of which included only one Ecuadorian isolate, had TMRCAs greater than 150 YBP (191, 208, and 229 YBP, respectively). TMRCAs were also estimated for the ungrouped strains, ranging from 157 to 404 YBP (Supplementary Fig. ). Because most Ecuadorian Mtb isolates were closely related to the 4.3.2/3 and 4.1.1 sub-lineages, we focused on all isolates belonging to these sub-lineages. The TMRCA estimates for the 4.3.2/3 and 4.1.1 sub-lineages were 470 YBP (95% HPD, 358 to 720) and 450 YBP (95% HPD, 349 to 647), respectively. Within the 4.1.1 sub-lineage, TGC_15 members have a TMRCA of 291 YBP, while those from TGC_37 have a TMRCA of 254 YBP. Isolates from TGC_14 had a TMRCA of 238 YBP, whereas ungrouped isolates S0516 and S2193 had TMRCA estimates of 482 and 354 YBP, respectively. In the 4.3.2/3 sub-lineage, TMRCA estimates ranged from 242 to 295 YBP, covering six TGCs: TGC_26 (289 YBP), TGC_18 (242–265 YBP), TGC_30 (249 YBP), TGC_20 (271–293 YBP), TGC_21 (295 YBP), and TGC_29 (292 YBP), and one ungrouped isolate (S2192, 289 YBP). Notably, isolate S2192 was included in the clade containing TGC_26, and TGC_18 was divided into three clades. In the 4.1.2.1 sub-lineage, isolates have a TMRCA estimated to be between 237 and 292 YBP, including 263 YBP for isolates corresponding to TGC_05, 237 to 273 YBP for TGC_03, and 292 YBP for TGC_02. The S-type sub-lineage had a TMRCA of 257 YBP (TGC_11), while isolates corresponding to the 4.8 sub-lineage had a TMRCA of 278 YBP (TGC_48) and 269 YBP (TGC_47). Isolate S0137, which formed a clade with TGC_48, had a TMRCA of 340 YBP (Supplementary Fig. ). 

This study focused on genomic analysis of raw reads of 88 M . tuberculosis isolates collected in Ecuador . To extend our investigation of genomic variability and epidemiology, we combined these data with a selection of 415 M . tuberculosis isolates from different countries (see Methods). Initially, we used the raw sequences (average coverage of 61X) of the 88 Ecuadorian isolates to perform quality control and genome assembly for each isolate using Unicycler and Pilon while discarding contigs smaller than 300 bp. The quality of the genomes was assessed and an average N50 of 65,747. Annotation was performed using Prokka software. These analyses yielded approximately 4298 genes per isolate, including coding sequences and RNA gene transfer. Descriptive statistics for this analysis are presented in Table and Supplementary Table . Genomic analysis of the Ecuadorian M. tuberculosis isolates revealed that they belonged to the Euro-American Lineage (lineage 4). Among the sublineages, 4.3.3 (27.3%, 24/88) and 4.1.1 (23.9%, 21/88) were the most prevalent, followed by 4.4.1.1 and 4.1.2.1 (11.4%, 10/88 each). Clades 4.1.2, and 4.3.4.1/2 were present, but in smaller proportions. 

We performed functional annotation of each of the 88 Ecuadorian genomes. A substantial proportion (77%) of the annotated genes was associated with specific functional roles (e.g., metabolic processes, cellular processes, energy, protein processing, and stress response, Defense, and Virulence). We enriched our dataset with functional annotations, assigning Enzyme Commission (EC) numbers to 1,061 proteins, Gene Ontology (GO) classifications to 918 proteins, linking 814 proteins to specific KEGG pathways, and identifying approximately 816 virulence factors distributed across three databases. This enrichment analysis provided deep insights into the metabolic functions and potential pathogenic mechanisms of the isolates (Supplementary Table ). Our subsystem analysis, categorized according to the KEGG database, showed that approximately 1,961 genes per isolate are involved in various biological processes and structural complexes (Supplementary Table ). The analysis revealed that the majority of the annotated protein-coding genes were related to metabolism (40.3%), with 37.9% of these genes being associated with cofactors, vitamins, and prosthetic groups (300 genes). This was followed by genes related to the stress response, Defense, Virulence (10.8%), and energy (10.5%). In our detailed metabolic pathway analysis, we observed a predominance of genes involved in amino acid (369 genes), carbohydrate (312 genes), lipid (262 genes), and xenobiotic biodegradation and metabolism (240 genes). These pathways are crucial for the pathogenicity and survival of Mtb (Supplementary Table – ). Notably, variability was observed among different sub-lineages. Sub-lineage 4.1.2.1 shows a higher proportion of proteins related to metabolism and a lower proportion of proteins related to protein processing. Conversely, sub-lineage 4.1.2 presents a higher proportion of proteins related to protein processing and fewer proteins related to metabolism (Supplementary Fig. A). Lineage 4.1.2.1 exhibits a higher number of gene variations related to metabolism, suggesting significant metabolic diversity, indicating a complex and adaptable metabolic network, whereas lineage 4.3.2/3 showed greater variations in genes associated with protein processing, suggesting a diverse set of genes involved in protein synthesis, folding, and degradation, reflecting the complexity of protein processing mechanisms (Table ). 

We also performed pan-genome analysis to understand the genetic diversity and conservation of 88 Mtb isolates from Ecuador. Pangenome reconstruction of the 88 Mtb isolates revealed 4,397 gene families, including 3,104 classified as core (present in at least 99% of isolates), 666 as accessory, and 270 as cloud gene families (Supplementary Fig. B and Supplementary Table ). Interestingly, 70.5% of the isolates belonged to core gene families, suggesting minimal within-genome variability and emphasizing the high level of genetic conservation among the isolates. According to pangenome calculations, a b value of 0.086 in the power-law regression model indicated a close pangenome for Mtb . Three isolates (S1454, S1477, and S1453) from sub-lineages 4.1.1 and 4.3.2/3 showed the highest number of cloud genes (261, 200, and 177, respectively). Most gene families within the core and accessory partitions belong to metabolic subsystems. Conversely, the majority of the unique gene families were classified within the Environmental Information Processing subsystem (Supplementary Fig. C,D). 

Furthermore, when analyzing the phylogenetic relationships and genetic diversity associated with resistance in Ecuadorian isolates (Fig. ), our findings revealed that sub-lineages 4.3.2/3 and 4.1.1, exhibited a higher prevalence of resistant variants (46.6% and 26.1%, respectively). We identified 42 single-nucleotide variants (SNVs) in resistance-related genes, with rpoB showing the most mutations. Additionally, SNVs were more frequently identified in genes related to intermediary metabolism and respiration (21.71%), followed by cell wall processes (21.34%) and conserved hypotheticals (20.44%). Exhaustive details of the SNV are provided in Supplementary Table . To gain a deeper insight into the genomic epidemiology of TB in Ecuador, we incorporated an additional 415 Mtb samples from 19 countries categorized as both continental and transoceanic migratory nations (See “ ” section). All these 503 Mtb genomes belonged to the L4 lineage, including 4.3.2/3 (35.4%), 4.1.2 (22.7%), 4.4.1 (12.7%), and 4.1.1 (10.7%) sub-lineages (Supplementary Table ). Moreover, 63.8% (321/503) of Mtb samples were genotypically susceptible to all anti-TB drugs. The remaining 36.2% (180/503) were resistant to at least one antibiotic. The clinical classification of these drug-resistant samples showed that 16.7% were identified as MDR-TB, 8.5% as HR-TB, 4.2% as pre-XDR-TB, and 1.6% as RR-TB. Table summarizes the drug-resistant clinical classification for the 503 Mtb samples per region and Table provides an overview of their drug-resistant canonical mutations. An extended list of canonical mutations is provided in Supplementary Table . Since multidrug-resistant TB (MDR-TB) and extensively drug-resistant TB (XDR-TB) can increase domestic transmission , we determined the frequency of resistance-associated canonical mutations in 503 Mtb isolates. The most frequent mutations were katG Ser315Thr (n = 92), rpoB Ser450Leu (n = 68), rpsL Lys43Arg (n = 21), embB Met306Ile (n = 17), pncA Gln10Pro (n = 14), and gyrA Ala90Val (n = 11), which confer resistance to INH, RIF, STR, EMB, PZA, and FQs, respectively. Other common canonical mutations included fabG 1 − 15C > T (n = 16), ahpC -74G > A (n = 15), and gid 329_330delTG (n = 16), which are associated with INH and STR resistance, respectively. 

To better understand the genetic diversity and potential transmission dynamics of TB in Ecuador, we performed a population structure analysis of 503 Mtb samples (Fig. ). Clusterization based on the rhierBAPS approach revealed four primary clusters at the initial hierarchical level, which were further divided into 15 sub-clusters at the second hierarchical level. Phylogenetic reconstruction of 503 isolates using 1,437 SNPs revealed four main clusters with different branch colors. Additionally, the sub-clusters (level 2), lineages, geographical distribution, and resistance are represented by colored concentric lines. Cluster 1 (pink) encompassed 128 isolates divided into three sub-clusters (1–3), containing samples belonging to 4, 4.4.1.1, 4.2.2, 4.6.2, and 4.7/8 sub-lineages. Cluster 2 (in orange) comprises 125 isolates, including three sub-clusters (4–6), with the majority belonging to the 4.1.2.1 sub-lineage, and some isolates belonging to the 4.1.2 sub-lineage (sub-cluster 6). Cluster 3 (in green) was composed of 178 isolates, distributed across three sub-clusters (7–9), primarily associated with 4.3.2, 4.3.3, 4.3.4.1, and 4.3.4.2, and, unexpectedly, a 4.1.2 isolate within sub-cluster 9. Finally, Cluster 4 (in blue) encompasses six sub-clusters (10–15), including 71 isolates predominantly from the 4.1.1, 4.1.1.1, and 4.1.1.3 sub-lineages. Notably, the isolates formed subclusters 13 and 15, with a small portion of subcluster 10 corresponding to the 4.1 sub-lineage. In addition, we analyzed the isolates for transmission clusters by pairwise comparisons using a minimum distance of 12 SNPs with MTBseq. Thus, strains in the same transmission genomic cluster (TGC) had fewer than 12 SNPs. Using this strategy, 92.8% (n = 467) of the isolates were classified into 51 TGCs, ranging from two to 63 members, whereas 35 isolates were not classified. Based on the number of members per TGC (see “ ” section), we classified seven members as small , 20 as medium , and 24 as large (Supplementary Table ). Most of the 51 TGCs resolved using this methodology were within Cluster 3 (37.3%) and Cluster 1 (29.4%), as determined using rhierBAPS (Fig. A). In particular, Mtb samples identified as susceptible were more prevalent in larger TGC than in resistant samples. Most pre-XDR isolates exhibited clonal relationships either among themselves or with isolates from other countries. Of the 35 isolates not grouped by MTBseq, 28.6% belonged to sub-lineage 4.3.2/3/4 and 25.7% to sub-lineage 4.8 (Fig. A and Table ). Upon closer examination of the five largest TGCs, specifically, TGC_1, TGC_2, TGC_9, TGC_18, and TGC_19, distinct geographical profiles were identified. Notably, TGC_1 (n = 63 isolates) predominantly comprised Mtb isolates from Paraguay (46.3%), Peru (19.0%), and Spain (15.9%). TGC_18 (n = 56) included isolates from Ecuador (41.0%) and Paraguay (26.8%), whereas TGC_19 (n = 31) consisted solely of isolates from Spain (38.7%). TGC_2 (n = 30) was predominantly isolated from Argentina (53.3%), whereas TGC_9 (n = 24) was isolated mainly from Canada (54.1%). Among the 88 Ecuadorian Mtb isolates, 84 were grouped into 15 of 51 identified TGCs. These included five medium-sized and ten large TGCs. Notably, ten of these TGCs (TGC_2, TGC_3, TGC_5, TGC_14, TGC_15, TGC_18, TGC_20, TGC_21, TGC_37, and TGC_47) included isolates from Ecuador alongside isolates from various countries, unveiling previously unrecognized potential continental and intercontinental connections. In particular, TGC_37 and TGC_14 exhibited a close genetic relationship among their members, showing distances of ≤ 9 SNPs and ≤ 6 SNPs, respectively. In addition, TGC_37 comprised five isolates from Ecuador and one from an Ecuadorian migrant in Spain, highlighting the impact of migration on TB spread. Conversely, TGC_14 grouped two isolates from Ecuador with five from Brazil, indicating active transmission between these nations (Supplementary Table ). The Ecuadorian isolates within each TGC exhibited a high degree of clonality. For example, approximately 90% of the Ecuadorian isolates in TGC_11 and 70% in TGC_15 and TGC_18 had a genetic distance of zero SNPs (Supplementary Tables , and Supplementary Figs. , ). Similarly, in some TGCs, many isolates remained distinct and did not cluster with other isolates (singletons) in various proportions. (Supplementary Table ). Regarding the sub-lineage composition, the most representative sub-lineages in the Ecuadorian isolates among the TGCs were 4.3.3 (46.9%, 23/49; TGC_18), 4.1.1 (32.7%, 16/49; TGC_15), and 4.4.1.1 (20.4%, 10/49; TGC_11). Concerning the potential transmission of isolates, we identified four pre-XDR TB cases in TGC_15 and three pre-XDR isolates in TGC_15 and TGC_18, corresponding to clonal samples (genetic distance < 1 SNP). Among these three TGCs, 41% were from previously treated patients with TB, 34% were from untreated patients, and 25% were from patients currently undergoing treatment. 

TB transmission networks, that is, the suggested transmission routes of TB considering the genetic distances between samples, their geographical location, and year of isolation, were inferred for the 503 Mtb isolates using TransFlow. Using this strategy, 46 transmission networks were identified; however, only 18 of these networks had more than three members, as shown in Fig. A. As can be seen, most of these transmission networks (with more than three members) correspond to different TGCs, meaning they present more than 12 SNPs between them. However, six transmission networks (7, 11–14, and 16) corresponded to TGC_18, indicating the low genetic distance between them (Fig. A). Considering the three clustering strategies employed (rhierBAPS, MTBseq, and TransFlow), we identified 17 convergent transmission groups to define potential TB transmission networks between Ecuadorian and non-Ecuadorian isolates. The geographic distribution of networks involving Ecuadorian isolates is shown in Fig. B, where TGCs that were confirmed and congruent with TransFlow are highlighted in red. Notably, some of these Ecuadorian samples (collected between 2019 and 2021) showed clonality with isolates from Colombia (n = 2) and Latin American migrants in Spain (n = 4), which were detected in 2014 and 2015, respectively (Fig. A, Supplementary Figs. , and Supplementary Table ). We identified several transmission networks showing connections among countries, including Paraguay, Argentina, and Brazil, with specific samples acting as potential index cases in the corresponding transmission networks (Supplementary Fig. ). The distribution of local and migrant Ecuadorian isolates revealed diverse representations spread across 12 TGCs (Fig. C). Some of these isolates showed clonality with samples from Colombia, Latin America, and Ecuadorian migrants in Spain, or were identified as potential index cases within transmission networks (Fig. C and Supplementary Fig. ). Notably, TGC18 and 19 exhibit diverse nationalities, with 10 and 8 different origins represented within these clusters, respectively. The short genetic distance between different TGCs indicated a close relationship, suggesting a possible joint origin. This relationship was supported by the genetic distances observed between the Ecuadorian isolates, implying a shared ancestry among the TGCs (Fig. D). These findings highlight the complex transmission dynamics of TB, with evidence of cross-border and international transmission networks involving Ecuadorian populations and those from other Latin American countries. The diversity of isolates and potential index cases underscore the need for targeted public health interventions to address the multifaceted nature of TB epidemiology. 

To characterize mutations in the genes involved in host adaptability, we analyzed SNPs in genes commonly associated with MTBC virulence in Ecuadorian isolates . A total of 303 SNPs were identified in these virulence genes, and all 88 isolates had at least one SNP each in the mce1F, mmpL4, phoR, ctpV, pepD, mce3F, fadD13, and nuoG genes. The pks12, fadD5, mce3C, pks12, nuoG, and katG genes were mutated with at least one SNP in more than 50% of the isolates. The top six genes with the most significant SNPs were plcA, plcB, pks7, pks12, phoR, and PPE46 (Supplementary Tables , ). Among the isolates, 46.6% (41/88) had more than 40 SNPs, with a maximum of 79 SNPs identified within virulence genes. Isolates corresponding to sub-lineages 4.1.1 (TGC_15 and TGC_37), 4.4.1.1 (TGC_11), and 4.3.2/3/4 (TGC_18, _20, _21, _26, _29, and _30) presented higher numbers of virulence-associated polymorphisms. Additionally, two members of TGC_14 (S0017 and S0039) and one ungrouped sample (S2193) showed a high number of mutations in the virulence genes (70, 71, and 74 SNPs, respectively). Interestingly, the five isolates associated with sub-lineage 4.7/8 showed the lowest number of SNPs in the virulence genes (15–25 SNPs). 

The divergence time to the most recent common ancestor (TMRCA) estimated for the Ecuadorian isolates fell within a density interval between 697 and 1,475 years, with an estimated time of 1,054 years. We also estimated TMRCA for different TGCs, identifying TGC_15 as the oldest (218 years before present, YBP) and TGC_03 as the youngest (130 YBP). Notably, isolates corresponding to TGC_02, TGC_21, and TGC_47, each of which included only one Ecuadorian isolate, had TMRCAs greater than 150 YBP (191, 208, and 229 YBP, respectively). TMRCAs were also estimated for the ungrouped strains, ranging from 157 to 404 YBP (Supplementary Fig. ). Because most Ecuadorian Mtb isolates were closely related to the 4.3.2/3 and 4.1.1 sub-lineages, we focused on all isolates belonging to these sub-lineages. The TMRCA estimates for the 4.3.2/3 and 4.1.1 sub-lineages were 470 YBP (95% HPD, 358 to 720) and 450 YBP (95% HPD, 349 to 647), respectively. Within the 4.1.1 sub-lineage, TGC_15 members have a TMRCA of 291 YBP, while those from TGC_37 have a TMRCA of 254 YBP. Isolates from TGC_14 had a TMRCA of 238 YBP, whereas ungrouped isolates S0516 and S2193 had TMRCA estimates of 482 and 354 YBP, respectively. In the 4.3.2/3 sub-lineage, TMRCA estimates ranged from 242 to 295 YBP, covering six TGCs: TGC_26 (289 YBP), TGC_18 (242–265 YBP), TGC_30 (249 YBP), TGC_20 (271–293 YBP), TGC_21 (295 YBP), and TGC_29 (292 YBP), and one ungrouped isolate (S2192, 289 YBP). Notably, isolate S2192 was included in the clade containing TGC_26, and TGC_18 was divided into three clades. In the 4.1.2.1 sub-lineage, isolates have a TMRCA estimated to be between 237 and 292 YBP, including 263 YBP for isolates corresponding to TGC_05, 237 to 273 YBP for TGC_03, and 292 YBP for TGC_02. The S-type sub-lineage had a TMRCA of 257 YBP (TGC_11), while isolates corresponding to the 4.8 sub-lineage had a TMRCA of 278 YBP (TGC_48) and 269 YBP (TGC_47). Isolate S0137, which formed a clade with TGC_48, had a TMRCA of 340 YBP (Supplementary Fig. ). 

Human-adapted MTB strains show a high degree of genomic conservation but vary in geographic distribution, virulence, transmissibility, and drug resistance patterns . To better understand the transmission pathways within our study population, we analyzed the sequences of 503 Mtb isolates of the L4 lineage, mostly from neighboring regions in Ecuador and other areas worldwide. It should be noted that the number of isolates varies significantly between countries. For instance, El Salvador and Chile have only one sample each, while countries such as Brazil and Paraguay have more than one hundred. This heterogeneity may affect the analysis of geographic distribution and limit the interpretation of the results. Our analysis of high-quality genomes revealed that 63.8% of the isolates were sensitive to all drugs used for TB treatment, whereas the remaining 36.2% were resistant to at least one drug. The genomes of all the isolates studied belonged to the Euro-American lineage, with the most common sub-lineages being 4.3.2/3 (35.4%), 4.1.2.1 (22.7%), 4.4.1 (12.7%), and 4.1.1. (10.7%). These results were congruent with those of previous studies in Ecuador that used the MIRU-VNTR strategy for genotyping circulating MTBC strains and also identified other Mtb lineages, including L2.1, and L4.2 , , . The distribution of these sub-lineages suggests that both historical European roots , , and migratory processes , , contributed to their presence in South and Central America and the Caribbean. Furthermore, we identified 19 genes that harbor mutations associated with resistance, predominantly linked to resistance to first-line drugs, thus enhancing the local genetic data for TB research , , . Functional characterization of genes is crucial for understanding how microorganisms such as M. tuberculosis adapt and survive in their hosts. Annotation of protein-encoding genes revealed that genes mainly associated with Cofactors, Vitamins, and Prosthetic groups, fatty acids, Lipids, and Isoprenoids, and Amino acids, and derivatives were most representative. Similar findings have been reported in other M. tuberculosis populations, where genes related to energy production and conversion, amino acid transport and metabolism, and lipid transport and metabolism are highly represented – . The wide conservation of these genes indicates their importance in interactions between bacteria and their human hosts. Specifically, during mycobacterial persistence, when the host–pathogen struggles for nutrient and immune recognition, these genes play a crucial role in ensuring the survival and adaptability of bacteria. On the other hand, we found that certain genes associated with the virulence of Mycobacterium tuberculosis like ESAT-6-like protein EsxS (63.6%), Acid and phagosome-regulated protein Apr AB (69.3%), and Chorismate mutase I were absent in the Ecuadorian isolates which would suggest possible affectation in adaptation and persistence of Mtb , due ESAT-6-like protein EsxS is related with the modulation of host immune responses , ; the Apr ABC locus modulates pH-driven adaptation to the macrophage phagosome , , and Chorismate mutase I is involved in inhibiting intrinsic apoptotic cell death of macrophages, playing a key role in the pathogenesis of TB , . Utilizing WGS in spanning network analysis enables researchers to connect TB transmission events , , providing insights into the spatial and temporal dynamics of TB transmission and identifying individuals and locations that play a critical role , . A recent population-based sequencing approach realizes a critical analysis of the utility of a pairwise distance threshold of < 12 SNPs and suggests that in scenarios with higher transmission rates, it is necessary to comprehend long-term transmission dynamics, adhering to a strict transmission 12 SNPs threshold is not advisable ; however, contrary to several studies supporting the utility of pairwise distance threshold of < 12 SNPs to identify recent transmission events – . We examined the possible transmission networks between Mtb isolates from Ecuador and from other 19 countries. Evidence suggests the potential spread of Ecuadorian Mtb isolates to individuals in different countries, based on the most significant clustering proportion confined to specific geographical locations. Combining genetic and epidemiological data could facilitate TB transmission management, particularly in migrant communities where socio-epidemiological changes due to migration may increase transmission complexity . Migration fluxes to and from high-burden countries significantly influence TB incidence, potentially leading to disease reactivation , . TB case clusters may involve autochthonous, mixed multinational, or cases among foreign-born individuals concentrated in a specific country , . Furthermore, Mtb can exhibit clonal transmission between hosts and establish clonal infections within a single host, with limited genetic diversification during infection or reactivation – . Our study observed clonality among Ecuadorian isolates with Colombian or Latin American migrants in Spain, indicating a potential transmission route involving direct contact between migrants, particularly in shared workplaces, such as plantations, factories, or restaurants. Surprisingly, three isolates corresponding to female patients were found in both scenarios, indicating that there might have been a relationship among the individuals, which might have been the cause of disease transmission. These findings suggest possible transnational transmission events involving Ecuador and its border countries and frequent migration destinations, highlighting the need to strengthen disease surveillance to reduce the possibility of more dangerous strains of tuberculosis entering the country , . Genome comparisons offer valuable insights into the molecular mechanisms that bacteria employ to survive and multiply within intracellular or extracellular host environments and to induce lesions and diseases . However, our understanding of the virulence factors expressed by Mtb is limited, and genetic variations resulting from selection pressure may affect their expression. Despite these challenges, performing this type of analysis may contribute to our understanding of how these factors function under local conditions . Our study identified 303 SNPs located within 103 genomic regions associated with virulence, with particular emphasis on genes such as plcA, plcB, pks7, pks12, phoR , and PPE46 which displayed the highest number of SNPs. These genes play pivotal roles in lipid metabolism, which is crucial for Mtb virulence of Mtb and is an integral component of the complex mycobacterial cell envelope , . Variations in these genes could affect their functionality, affecting the ability of the bacterium to evade host immune defense, induce necrosis in macrophages, and modulate virulence, particularly the PhoP-PhoR two-component system, potentially contributing to the pathogenic capabilities of the bacterium , . Phylodynamic analysis is useful for understanding TB strain variation and dynamics in various countries . The TMCRA and many TGC datasets suggest that TB isolates sampled in Ecuador can trace their ancestry back hundreds of years, indicating a complex evolutionary situation and implying the possibility that different TB isolates have caused infections at various time intervals from varied sources. The Ecuadorian Mtb data reported frequencies of 4.3.2/3 and 4.1.1, with TMRCA estimated to range from 450 to 470 YBP. These data indicate the possible ancient roots of the isolates. Nonetheless, it should be mentioned that these calculations only depict TB ancestral isolates within each sub-lineage, and the origin of TB in Ecuador could have been even earlier, before the arrival of Europeans, indicating a long history of human-pathogen co-evolution in the region, consistent with that reported in the ancient Andean population . While some isolates appear to have signs of long-living existence and evolution, and thus historical transmission, others seem to have appeared more recently, suggesting current transmission. Genetic diversity among Ecuadorian Mtb isolates is likely to arise from transcontinental interactions, human migration, and various factors that promote the circulation of local and global TB isolates. Thus, additional studies integrating genomic data from various geographic sites as well as comprehensive epidemiological information could help trace the origin and dissemination routes of TB isolates in Ecuador more accurately. Although our study offers valuable insights, it had several limitations. The primary limitation is the modest sample size of 88 Mtb isolates from cultured samples, a small fraction compared to the total number of TB cases reported in Ecuador in 2021 (5595). Furthermore, the lack of comprehensive epidemiological information, including contact tracing details, restricted our analysis of transmission dynamics. To enhance our understanding, future research should aim to combine genomic data with additional epidemiological details to uncover potential Mtb transmission pathways for Mtb . Moreover, our sequencing focused on cultured sputum isolates, a standard approach in Mtb genomic epidemiology that may only partially capture the complete spectrum of mycobacteria in the lungs, thus limiting our ability to capture the intricacies of within-host variations in individual infections Longitudinal studies have the potential to enhance sequencing analyses and to uncover variations in resistance genes and virulence factors, thereby assisting in refining treatment strategies. These insights hold promise for alleviating the global TB burden. To the best of our knowledge, this is the first study to establish local Mtb transmission networks in Ecuador using whole-genome analysis. Our findings reinforce and contribute to the knowledge of transmission networks previously characterized in Ecuador based on the MIRU-VNTR approach , , , . Our study provides valuable insights into the genomic and epidemiological characteristics of Mtb isolates from Ecuador. Additionally, our analysis identified drug-resistant isolates and transmission events between individuals and across borders, underscoring the need for more extensive whole-genome sequencing and network analyses to guide public health interventions. 

Genome database of the Mtb samples This study analyzed raw reads from the genomes of 88 clinical Mtb isolates of the L4 lineage, sequenced in a previous study , corresponding to BioProject PRJNA827129. These isolates were collected conveniently between 2019 and 2021 from private laboratories and the National Reference of Mycobacteria at the National Institute of Public Health Research "Leopoldo Izquieta Pérez" (INSPI-LIP) across the different provinces in Ecuador. A significant proportion of the samples came from Guayaquil, accounting for 81.8% of the isolates. This city represents epidemiologically more than half of the TB cases in the country . The remaining samples were collected from Babahoyo (5.6%), El Empalme (3.4%), and Quito (2.3%), highlighting the geographical spread and prevalence of TB in these regions. A small number of isolates came from Chone, Duran, Guaranda, Machala, and Nueva Loja, accounting for 1.1% of the total isolates. Additionally, the study included 415 publicly available sequences previously characterized as Mtb isolates of the L4 lineage from 19 countries identified as significant in the Ecuadorian context of continental and transoceanic migration. Among the continental countries, the study included isolates from Argentina (n = 18), Brazil (n = 84), Canada (n = 41), Colombia (n = 8), Guatemala (n = 16), Mexico (n = 35), Panama (n = 4), Paraguay (n = 67), Peru (n = 44), and the USA (n = 6). Isolates from Hungary (n = 5), the Netherlands (n = 12), Portugal (n = 8), Spain (n = 59), and the United Kingdom (n = 5) were included. Notably, among the 59 Spanish isolates, some were identified previously from Latin American migrants who had settled in Spain years earlier, including individuals from Bolivia (n = 20), Colombia (n = 8), Ecuador (n = 8), Chile (n = 1), and Honduras (n = 1), as detailed in reference . The accession numbers and distributions of these countries are listed in Supplementary Table , . All protocols used in this article were approved by The University Espiritu Santo Review Board under code 2022-001A. Genome assembly and annotation To ensure high-quality genomic data, raw reads from the 503 Mtb genomes were initially processed with rigorous quality control using Fastp v0.23.4 , and Kraken v2 for species confirmation and contamination screening, ensuring that only Mtb -specific reads were processed. Pangenome construction Pangenomic analysis was carried out using the Panaroo pipeline from the GFF archives annotated by Prokka v1.14.16 using the H37Rv (NC_000962.3) M. tuberculosis reference genome and default parameters for clustering to define gene families and identify core genes present in 99% of the isolates. The pangenome was divided into core, accessory, and unique genes, which were categorized based on their presence in isolates. Variant calling analysis The cleaned reads were processed using the MTBseq pipeline with standard input parameters to map the reads to the Mtb H37Rv reference genome (NC_000962.3). Briefly, this tool involves BWA-mem and SAMtools for mapping, GATK v3 for base call recalibration, and realigning reads around insertions and deletions (InDels), followed by SAMtools mpileup for Single Nucleotide Polymorphisms (SNPs) and InDel calling (using parameters B and d 1000). Only high-quality genomes were processed after accomplishing the following criteria: mean coverage greater than 20x, read depth less than 5x, and reference genome coverage of > 95%. MTBseq analysis facilitated sub-lineage classification and the construction of a genetic distance matrix to identify transmission groups. Virulence analysis Ecuadorian Mtb isolates were analyzed to identify mutations in the genes involved in host adaptability and virulence. Genomes from the isolates were assembled from high-quality reads using the Unicycler assembly pipeline with a minimum contig size of 300 bp and polished using Pilon . SNPs within high-confidence and repetitive genomic regions, such as the PE/PPE gene family, were confirmed via visualization using Integrated Genome Viewer , . Virulence-related proteins were analyzed using the Virulence Factor Database (VFDB) in the Pathosystems Resource Integration Center (PATRIC) online ( https://www.patricbrc.org/ ). Drug resistance prediction and lineage classification The TB-Profiler v4 pipeline uses the BCF tool for variant calling and predicts drug resistance-associated mutations. These include resistance to first-line drugs such as rifampicin (RIF), isoniazid (INH), pyrazinamide (PZA), and ethambutol (EMB), and second-line drugs such as fluoroquinolones (FQs), streptomycin (STR), ethionamide (ETH), and aminoglycosides (including second-line drugs, amikacin, kanamycin, and capreomycin). According to World Health Organization guidelines, isolates exhibiting resistance solely to INH are classified as isoniazid-resistant TB (HR-TB), those resistant to RIF as rifampicin-resistant TB (RR-TB), and those with resistance to both INH and RIF as multidrug-resistant TB (MDR-TB). Isolates that demonstrated resistance to any FQ in addition to MDR or RR status were categorized as pre-extensively drug-resistant TB (pre-XDR-TB). Phylogenetic analysis and population structure of Mtb isolates Phylogenetic reconstruction was conducted using concatenated SNPs aligned genomically using the MTBseq tool. Repetitive regions (PPE/PE-PGRS genes), consecutive indels, and genes involved in antibiotic resistance were excluded from further phylogenetic analysis. The optimal substitution model for this SNP alignment was determined using ModelTest-NG v0.1.7 . Following this, the phylogenetic tree was reconstructed using the maximum-likelihood (ML) method via RAxML-NG , utilizing the general time-reversible model with gamma-distributed rate heterogeneity (GTR + GAMMA) and supported by 1000 bootstrap replicates. To address the potential ascertainment bias resulting from the exclusive use of polymorphic sites and consequent rescaling of tree branches, ascertainment bias correction was applied by specifying the number of invariant sites , as outlined at https://github.com/conmeehan/pathophy . The phylogenetic tree was visualized using the Interactive Tree Of Life (iTOL) v6.6 . The genome of Mycobacterium microti genome (accession number: SRR3647357) was used as an outgroup to root the tree. Additionally, a population structure analysis was conducted using the rhierBAPS package , which was set to a maximum depth of two and n.pops of 20. This analysis employed a hierarchical nested clustering approach based on genetic data, specifically prioritizing SNP loci that displayed a minor allele in at least two sequences, to effectively identify subpopulations or clusters. Transmission cluster analysis and modeling of the genetic clustering network Transmission genomic clusters (TGCs) were identified using concatenated high-quality SNPs processed using MTBseq. We applied a pairwise distance threshold of 12 SNPs between isolates to achieve optimal population clustering, a standard established in previous genomic TB studies – , . TGCs were categorized based on their size into small (fewer than three isolates), medium (three to five isolates), and large (more than five isolates). Visualization of the resulting SNP alignment for each cluster was used to infer a genetic network. We used a parsimony-based algorithm for network reconstruction using PopART software because of the monomorphic and non-recombining behavior of Mtb, as well as the potential of the sample dataset, including the original genotype. We compared the distribution network of isolates using a median-spanning network (MSN) and median-joining network (MJN). For clusters of at least three samples, we utilized TransFlow to reconstruct the transmission network, enhancing our understanding of local transmission dynamics. TransFlow integrates genomic data with epidemiological factors, such as sampling dates and geographic coordinates, to map spatial connectivity among isolates. Additionally, to mitigate the bias arising from the use of lineage-specific reference genomes, TransFlow incorporates the PANPASCO pipeline . This pipeline employs a computational pan-genome consisting of 146 complete MTBC genomes from major lineages 1–4, facilitating accurate pairwise SNP distance calculations. This methodology ensures a thorough and representative analysis of genetic variation across MTBC populations. Phylodynamics of TB transmission clusters To study the temporal dynamics of the different sub-lineages identified in the Ecuadorian Mtb L4 isolates, with a particular focus on uncovering historical introduction events, a time-calibrated phylogeny was inferred using BEAST v1.10.4 , utilizing collection and tip dates from the isolates. The XML input file necessary for analysis was generated by concatenating SNPs derived from MTBseq and processed using BEAUTi. This file was adjusted to specify the number of invariant sites following guidance from the BEAST User Forum. To assess the temporal signal of the sequence alignments, we used TempEst v1.5.3 . The dating analysis employed the general time-reversible plus gamma distribution (GTR + GAMMA) substitution model coupled with a strict molecular clock and coalescent constant-size demographic model. Markov chain Monte Carlo (MCMC) simulations were conducted for 250 million iterations, with a 10% burn-in phase and samples taken every 10,000 generations. This approach facilitates independent evaluation of chain convergence. The analysis results were summarized and convergence was confirmed using Tracer v1.6 , ensuring that all essential parameters achieved an effective sample size (ESS) of over 200. The Maximum Clade Credibility (MCC) tree was then computed using TreeAnnotator v2.5.0, providing a statistically supported phylogenetic tree with time-calibrated estimates. 

This study analyzed raw reads from the genomes of 88 clinical Mtb isolates of the L4 lineage, sequenced in a previous study , corresponding to BioProject PRJNA827129. These isolates were collected conveniently between 2019 and 2021 from private laboratories and the National Reference of Mycobacteria at the National Institute of Public Health Research "Leopoldo Izquieta Pérez" (INSPI-LIP) across the different provinces in Ecuador. A significant proportion of the samples came from Guayaquil, accounting for 81.8% of the isolates. This city represents epidemiologically more than half of the TB cases in the country . The remaining samples were collected from Babahoyo (5.6%), El Empalme (3.4%), and Quito (2.3%), highlighting the geographical spread and prevalence of TB in these regions. A small number of isolates came from Chone, Duran, Guaranda, Machala, and Nueva Loja, accounting for 1.1% of the total isolates. Additionally, the study included 415 publicly available sequences previously characterized as Mtb isolates of the L4 lineage from 19 countries identified as significant in the Ecuadorian context of continental and transoceanic migration. Among the continental countries, the study included isolates from Argentina (n = 18), Brazil (n = 84), Canada (n = 41), Colombia (n = 8), Guatemala (n = 16), Mexico (n = 35), Panama (n = 4), Paraguay (n = 67), Peru (n = 44), and the USA (n = 6). Isolates from Hungary (n = 5), the Netherlands (n = 12), Portugal (n = 8), Spain (n = 59), and the United Kingdom (n = 5) were included. Notably, among the 59 Spanish isolates, some were identified previously from Latin American migrants who had settled in Spain years earlier, including individuals from Bolivia (n = 20), Colombia (n = 8), Ecuador (n = 8), Chile (n = 1), and Honduras (n = 1), as detailed in reference . The accession numbers and distributions of these countries are listed in Supplementary Table , . All protocols used in this article were approved by The University Espiritu Santo Review Board under code 2022-001A. 

To ensure high-quality genomic data, raw reads from the 503 Mtb genomes were initially processed with rigorous quality control using Fastp v0.23.4 , and Kraken v2 for species confirmation and contamination screening, ensuring that only Mtb -specific reads were processed. 

Pangenomic analysis was carried out using the Panaroo pipeline from the GFF archives annotated by Prokka v1.14.16 using the H37Rv (NC_000962.3) M. tuberculosis reference genome and default parameters for clustering to define gene families and identify core genes present in 99% of the isolates. The pangenome was divided into core, accessory, and unique genes, which were categorized based on their presence in isolates. 

The cleaned reads were processed using the MTBseq pipeline with standard input parameters to map the reads to the Mtb H37Rv reference genome (NC_000962.3). Briefly, this tool involves BWA-mem and SAMtools for mapping, GATK v3 for base call recalibration, and realigning reads around insertions and deletions (InDels), followed by SAMtools mpileup for Single Nucleotide Polymorphisms (SNPs) and InDel calling (using parameters B and d 1000). Only high-quality genomes were processed after accomplishing the following criteria: mean coverage greater than 20x, read depth less than 5x, and reference genome coverage of > 95%. MTBseq analysis facilitated sub-lineage classification and the construction of a genetic distance matrix to identify transmission groups. 

Ecuadorian Mtb isolates were analyzed to identify mutations in the genes involved in host adaptability and virulence. Genomes from the isolates were assembled from high-quality reads using the Unicycler assembly pipeline with a minimum contig size of 300 bp and polished using Pilon . SNPs within high-confidence and repetitive genomic regions, such as the PE/PPE gene family, were confirmed via visualization using Integrated Genome Viewer , . Virulence-related proteins were analyzed using the Virulence Factor Database (VFDB) in the Pathosystems Resource Integration Center (PATRIC) online ( https://www.patricbrc.org/ ). 

The TB-Profiler v4 pipeline uses the BCF tool for variant calling and predicts drug resistance-associated mutations. These include resistance to first-line drugs such as rifampicin (RIF), isoniazid (INH), pyrazinamide (PZA), and ethambutol (EMB), and second-line drugs such as fluoroquinolones (FQs), streptomycin (STR), ethionamide (ETH), and aminoglycosides (including second-line drugs, amikacin, kanamycin, and capreomycin). According to World Health Organization guidelines, isolates exhibiting resistance solely to INH are classified as isoniazid-resistant TB (HR-TB), those resistant to RIF as rifampicin-resistant TB (RR-TB), and those with resistance to both INH and RIF as multidrug-resistant TB (MDR-TB). Isolates that demonstrated resistance to any FQ in addition to MDR or RR status were categorized as pre-extensively drug-resistant TB (pre-XDR-TB). 

Phylogenetic reconstruction was conducted using concatenated SNPs aligned genomically using the MTBseq tool. Repetitive regions (PPE/PE-PGRS genes), consecutive indels, and genes involved in antibiotic resistance were excluded from further phylogenetic analysis. The optimal substitution model for this SNP alignment was determined using ModelTest-NG v0.1.7 . Following this, the phylogenetic tree was reconstructed using the maximum-likelihood (ML) method via RAxML-NG , utilizing the general time-reversible model with gamma-distributed rate heterogeneity (GTR + GAMMA) and supported by 1000 bootstrap replicates. To address the potential ascertainment bias resulting from the exclusive use of polymorphic sites and consequent rescaling of tree branches, ascertainment bias correction was applied by specifying the number of invariant sites , as outlined at https://github.com/conmeehan/pathophy . The phylogenetic tree was visualized using the Interactive Tree Of Life (iTOL) v6.6 . The genome of Mycobacterium microti genome (accession number: SRR3647357) was used as an outgroup to root the tree. Additionally, a population structure analysis was conducted using the rhierBAPS package , which was set to a maximum depth of two and n.pops of 20. This analysis employed a hierarchical nested clustering approach based on genetic data, specifically prioritizing SNP loci that displayed a minor allele in at least two sequences, to effectively identify subpopulations or clusters. 

Transmission genomic clusters (TGCs) were identified using concatenated high-quality SNPs processed using MTBseq. We applied a pairwise distance threshold of 12 SNPs between isolates to achieve optimal population clustering, a standard established in previous genomic TB studies – , . TGCs were categorized based on their size into small (fewer than three isolates), medium (three to five isolates), and large (more than five isolates). Visualization of the resulting SNP alignment for each cluster was used to infer a genetic network. We used a parsimony-based algorithm for network reconstruction using PopART software because of the monomorphic and non-recombining behavior of Mtb, as well as the potential of the sample dataset, including the original genotype. We compared the distribution network of isolates using a median-spanning network (MSN) and median-joining network (MJN). For clusters of at least three samples, we utilized TransFlow to reconstruct the transmission network, enhancing our understanding of local transmission dynamics. TransFlow integrates genomic data with epidemiological factors, such as sampling dates and geographic coordinates, to map spatial connectivity among isolates. Additionally, to mitigate the bias arising from the use of lineage-specific reference genomes, TransFlow incorporates the PANPASCO pipeline . This pipeline employs a computational pan-genome consisting of 146 complete MTBC genomes from major lineages 1–4, facilitating accurate pairwise SNP distance calculations. This methodology ensures a thorough and representative analysis of genetic variation across MTBC populations. 

To study the temporal dynamics of the different sub-lineages identified in the Ecuadorian Mtb L4 isolates, with a particular focus on uncovering historical introduction events, a time-calibrated phylogeny was inferred using BEAST v1.10.4 , utilizing collection and tip dates from the isolates. The XML input file necessary for analysis was generated by concatenating SNPs derived from MTBseq and processed using BEAUTi. This file was adjusted to specify the number of invariant sites following guidance from the BEAST User Forum. To assess the temporal signal of the sequence alignments, we used TempEst v1.5.3 . The dating analysis employed the general time-reversible plus gamma distribution (GTR + GAMMA) substitution model coupled with a strict molecular clock and coalescent constant-size demographic model. Markov chain Monte Carlo (MCMC) simulations were conducted for 250 million iterations, with a 10% burn-in phase and samples taken every 10,000 generations. This approach facilitates independent evaluation of chain convergence. The analysis results were summarized and convergence was confirmed using Tracer v1.6 , ensuring that all essential parameters achieved an effective sample size (ESS) of over 200. The Maximum Clade Credibility (MCC) tree was then computed using TreeAnnotator v2.5.0, providing a statistically supported phylogenetic tree with time-calibrated estimates. 

Supplementary Information 1. Supplementary Information 2. 

