Local adaptation in populations of Mycobacterium tuberculosis endemic to the Indian Ocean Rim

Background: Lineage 1 (L1) and 3 (L3) are two lineages of the Mycobacterium tuberculosis complex (MTBC) causing tuberculosis (TB) in humans. L1 and L3 are prevalent around the rim of the Indian Ocean, the region that accounts for most of the world’s new TB cases. Despite their relevance for this region, L1 and L3 remain understudied. Methods: We analyzed 2,938 L1 and 2,030 L3 whole genome sequences originating from 69 countries. We reconstructed the evolutionary history of these two lineages and identified genes under positive selection. Results: We found a strongly asymmetric pattern of migration from South Asia toward neighboring regions, highlighting the historical role of South Asia in the dispersion of L1 and L3. Moreover, we found that several genes were under positive selection, including genes involved in virulence and resistance to antibiotics. For L1 we identified signatures of local adaptation at the esxH locus, a gene coding for a secreted effector that targets the human endosomal sorting complex, and is included in several vaccine candidates. Conclusions: Our study highlights the importance of genetic diversity in the MTBC, and sheds new light on two of the most important MTBC lineages affecting humans.


Introduction
The global tuberculosis (TB) epidemic is a major public health emergency, disproportionately affecting vulnerable populations and worsening existing inequalities. Although the estimated incidence and the number of fatalities have slowly decreased over time, every year ten million people develop the disease, and 1.4 million TB patients die (WHO, 2020). Moreover, the Covid-19 pandemic will likely set back the progress toward TB eradication for several years (Glaziou, 2020;Hogan et al., 2020;McQuaid et al., 2020;Saunders & Evans, 2020). TB is spread worldwide, but not all regions are equally affected: 95% of new TB cases occur in Africa and Asia, and the eight countries with the largest burden account for two thirds of the total number of cases (WHO, 2020).
Human TB is caused by members of the Mycobacterium tuberculosis complex (MTBC), which includes nine phylogenetic lineages with different geographic distributions. Five of these lineages are restricted to Africa (L5, L6, L7, L8, and L9), while the remaining four (L1, L2, L3, L4) are more broadly distributed (Chihota et al., 2018;Coscollá et al., 2021;Gagneux, 2018;Ngabonziza et al., 2020). While L2 and L4 strains occur across the world, L1 strains are predominantly found around the rim of the Indian Ocean (East Africa, South Asia, and Southeast Asia). L3 has a distinct geographic range (East Africa, Central Asia, Western Asia, and South Asia) that overlaps partially with L1 (Couvin et al., 2019;Wiens et al., 2018). While many MTBC lineages also occur in Northern Europe, North America, Australia and New Zealand, in these low-burden regions, the majority of TB cases are imported through recent migrations, and local transmission is rare (White et al., 2017).
Beyond their geographic distribution, MTBC lineages also differ in virulence, transmissibility, association with drug resistance, and the host immune responses they elicit (Peters et al., 2020). There is an increasing notion that MTBC genetic variation should be considered in the development of new antibiotics and vaccines, and when studying the epidemiology and pathogenesis of the disease (Gagneux, 2017;Peters et al., 2020;Wiens et al., 2018). Yet, much of TB research to date has focused on the laboratory strain H37Rv (belonging to L4) and a few other strains belonging to L2 and L4, because of their broad geographic range, and their association with drug resistance. By contrast, L1 and L3 have largely been neglected, and only few studies have investigated the global populations of these two lineages (Couvin et al., 2019;O'Neill et al., 2019). This knowledge gap is particularly severe as L1 and L3 are endemic to some of the world regions with the heaviest TB burden. For example, L1 and L3 cause the majority of TB cases in India, the country with the highest number of new TB cases in the world, and L1 is by far the most important cause of TB in the Philippines, the country with the 4 th highest global TB burden (Netikul et al., 2021;Wiens et al., 2018). The aim of this study was to characterize the global population structure of L1 and L3 using large-scale population genomics, and to investigate the evolutionary history and selective forces acting on these two lineages.

Results
We screened a large collection of publicly available MTBC genome sequences and selected those belonging to L1 and L3. Additionally, we newly sequenced 767 strains to cover regions that were not well represented in the public data. We applied a series of bioinformatic filters to exclude low quality sequences and mixed infections (Methods), and obtained a curated genome dataset of 4,968 strains (2,938 L1, 2,030 L3). While the phylogenetic tree of L3 showed a ladder-like topology without an evident division in sublineages, L1 comprised five clearly distinct sublineages coalescing close to the root of the tree (Extended data: Figures 1 and 2).
The genome sequences included in our final dataset represented 69 countries and covered the known geographic range of L1 and L3 (Methods; Extended data: Table 1 and Figures 3  and 4). Because we were interested in studying L1 and L3 in their endemic range, we excluded strains originating from North America, Europe, and Australia from the biogeographic analysis, as in these low-burden regions, most TB patients are recent migrants who were infected in their country of origin (White et al., 2017). We assigned strains to different geographic regions following a modified version of the United Nation geographic scheme (Methods, Figure 1 and Figure 2). Mapping the regions of origin onto the phylogenetic trees revealed that two sublineages of L1 are found almost exclusively in Southeast Asia (L1.1.1 and L1.2.1), while the others are spread over the complete geographic range of L1 ( Figure 1; Extended Data: Figures 5-6).
We performed a phylogeographic analysis using two alternative approaches (Mascot and PASTML), which are based on different models and inference methods (Ishikawa et al., 2019;Müller et al., 2018). We found that South Asia was predicted to be the geographic range of the Most Recent Common Ancestor (MRCA) of both L1 and L3 (Extended Data: Information). Most interestingly, we found a strongly asymmetric pattern of migration. For L1, PASTML identified most migration events from South Asia toward other regions, followed by further dispersion, but almost no back migration to South Asia ( Figure 1; Extended data: Figure 7 and Files 1 and 2). When we estimated the migration rates with Mascot, we found that the forward migration rate from South Asia toward the rest of the world were 3 to 17 times larger than the migration rates in the opposite direction, confirming the results of PASTML (Figure 1; Extended data: Table 2). We found a similar scenario for L3: PASTML inferred the largest number of migrations from South Asia toward East Africa, further spread from East Africa

Amendments from Version 1
In this revised version we amended some sections to enhance clarity, and we provided more details about the methods.
We also added a short comparison of the biogeography of L1 and L3 with that of other MTBC lineages, and included some reference to relevant studies that were missing.
Any further responses from the reviewers can be found at the end of the article Figure 1. Results of the biogeography analysis: L1. a) Heatmap indicating the origin of the 2,061 L1 strains included in the dataset used for the biogeography analysis. b) The geographic regions used in the biogeography analysis of L1. c) Results of PASTML (compressed tree), the color code follow the legend of panel (b), the size of the circles is proportional to the number of tips, and the size of the arrows is proportional to the number of times the pattern of migration was observed on the tree. This plot excluded nodes that represented less than three strains. The same plot including all nodes can be found as Extended Data File 1. d) Tree obtained with the Mascot analysis; the branches are colored following the legend in panel (b) and represent the inferred ancestral region with the largest posterior probability. e) Representation of the relative effective population sizes (circles) and migration rates (lines connecting circles) estimated by Mascot. Lines representing migration rates are colored based on the region of origin (interpreted forward in time). For example, the dark blue line connecting the dark blue circle with the red circle represents the forward migration rate between South Asia and East Africa, which is the largest migration rate estimated by Mascot. The parameter estimates are reported in Extended Data Table 2. toward neighboring regions, but essentially no migration back to South Asia ( Figure 2; Extended data: Files 3-4). The Mascot analysis showed that the forward migration rates from South Asia toward East Africa were 26 times larger than the migration rates in the opposite direction ( Figure 2; Extended data: Table 2).
We performed tip dating to estimate the age of the trees but got inconsistent results for L1, due to a lack of a reliable temporal signal (Extended data: Information, Figures 8 and 9, Tables 3 and 4). However, the results of a previous study suggested a relatively fast evolutionary rate for L1 (~ 1.4×10 -7 nucleotide changes per site per year; Menardo et al., 2019), and that the MRCA of L1 lived around the 12 th century AD (Extended data: Information). By contrast, L3 had a good temporal signal, and different methods estimated that its MRCA lived between the 2 nd and the 13 th century AD. However, the uncertainty around all of these estimates was very large (Extended data: Information, Figure 8, Tables 5 and 6). Together, these results corroborate the findings of previous studies (Bos et al., 2014;Duchêne et al., 2016;Menardo et al., 2019). Calibrating MTBC trees that are hundreds or thousands of years old, with sequences sampled in the last few decades is notoriously challenging and subject to limitations (Menardo et al., 2019). Therefore, we refrain from any strong interpretation of the results of the molecular clock analyses of L1 and L3. We emphasize that the ages reported here are the most likely estimates supported by the available data. Additional data, or alternative methods, might result in different temporal scenarios.
With respect to the geographical aspects, we identified several interesting instances of potential long-range dispersal. First, we found a clade of L1, composed of 11 strains sampled in five different West African countries. This was surprising because West Africa is usually not considered part of the geographic range of L1. This clade is nested within sublineage L1.1.1, which is essentially only found in mainland Southeast Asia, and the PASTML analysis inferred a direct introduction from Southeast Asia to West Africa. This introduction is unlikely to have happened before the 16 th century, when the Portuguese established the maritime route between Europe and Asia by circumnavigating Africa. We benchmarked the molecular clock of L1 and found that this scenario is indeed compatible with a clock rate equal to, or larger than 1.4×10 -7 nucleotide changes per site per year, but not with lower ones (Extended data: Information, Figure 10).
Second, we found that L1 was introduced to South America on at least 11 independent occasions (Figure 1; Extended data: Files 1 and 2). Assuming a clock rate of 1.4×10 -7 , the earliest introduction was between 1620 and 1830 AD from East Africa, while subsequent introductions occurred from East Africa, South Africa, and South Asia. These results support the hypothesis that L1 was first introduced to Brazil through the slave trade from East Africa (Allen, 2013;Conceição et al., 2019), although, due to the lack of samples from other South American countries, we cannot exclude a more complex evolutionary history of L1 in this continent. Interestingly, this is in contrast with L5 and L6, which are endemic to West Africa, but did not establish themselves firmly in South America (De Jong et al., 2010;Rabahi et al., 2020).
Third, similar to the West African clade of L1.1.1, we found an East African clade embedded within sublineage L1.2.1, which otherwise is found almost exclusively in Southeast Asia. This East African clade is composed of 11 strains from five countries, and its sister clade is found in East Timor and Papua New Guinea. We inferred a direct migration from the islands of Southeast Asia to East Africa that occurred between the 13 th and the 16 th century AD (assuming a clock rate of 1.4×10 -7 ). This would be compatible with early Portuguese expeditions, which reached East Timor and Papua New Guinea in the early 16 th century. An alternative explanation could be the Austronesian expansion. Austronesians are thought to have reached the Comoros islands and Madagascar between the 9 th and the 13 th century AD, possibly through direct navigation from Southeast Asian islands. Malagasy speak an Austronesian language, and Austronesian genetic signatures are found in human populations in the Comoros, Madagascar, and to a small extent also in the Horn of Africa (Blench, 2010;Boivin et al., 2013;Brucato et al., 2016;Brucato et al., 2018;Brucato et al., 2019;Pierron et al., 2014).
L1 and L3 coexist in many regions around the Indian Ocean. Yet, in their evolutionary history, these lineages colonized areas occupied by different human populations. Human genetic The geographic regions used in the biogeography analysis of L3. c) Results of PASTML (compressed tree), the color code follows the legend of panel (b), the size of the circles is proportional to the number of tips, and the size of the arrows is proportional to the number of times the pattern of migration was observed on the tree. This plot excluded nodes that represent less than three strains, the same plot including all nodes can be found as Extended Data File 3. d) Tree obtained with the Mascot analysis; the branches are colored following the legend in panel (b) and represent the inferred ancestral region with the largest posterior probability. e) Representation of the relative effective population sizes (circles) and migration rates (lines connecting circles) estimated by Mascot. Lines representing migration rates are colored based on the region of origin (interpreted forward in time). The effective population size of South Asia was estimated to be much larger than the one of East Africa. Additionally, the forward migration rate between South Asia and East Africa was much larger than the one in the opposite direction. The parameter estimates are reported in Extended Data Table 2. variation has been shown to influence the susceptibility to TB (Qu et al., 2011). Most notably, human leukocyte antigen (HLA) genes play a crucial role in the activation of the immune responses to the MTBC by exposing bacterial peptides (epitopes) to the surface of an infected cell, where they can be recognized by T cells. HLA genes are extremely polymorphic in human populations, and several alleles of different HLA genes are associated with TB susceptibility in different populations (Brahmajothi et al., 1991;Salie et al., 2014;Sveinbjornsson et al., 2016;Vejbaesya et al., 2002;Yuliwulandari et al., 2010).
Previous studies have shown that T cell epitopes are hyperconserved in the MTBC, suggesting that immune escape does not provide an advantage and that, contrary to other pathogens, the MTBC needs to be recognized by the immune system and to cause disease in order to transmit (Comas et al., 2010;Coscolla et al., 2015). Our large dataset of L1 and L3 genome sequences from different geographic regions provided an opportunity to scan for lineage-and/or region-specific signatures of selection at T cell epitopes in L1 and L3.
We reconstructed the mutational history of T cell epitopes in L1 and L3, and found that 51% of all epitopes were variable (at the amino acid level) in at least one L1 strain, while only 20% were variable in at least one L3 strain (Extended data: Table 7). However, this difference can be explained by the different size and diversity of the two datasets (2,061 genome sequences with 136,023 variable sites for L1; 1,021 genome sequences with 36,316 variable sites for L3). The epitope that accumulated most mutations was located at the N terminus of esxH (Rv0288). This peptide is a T cell epitope for both classes, MHCI and MHCII; it is also a B cell epitope and was previously identified as one of the few T cell epitopes that were not hyper-conserved (Comas et al., 2010). We found 15 derived haplotypes (at the amino acid level), generated by 28 independent replacements at five positions in a peptide of seven amino acids ( Figure 3). By contrast, the second most variable epitope accumulated only seven amino acid changes (Extended data: Table 7). Interestingly, we did not find this signature in L3, as all strains carried the ancestral haplotype at the N-terminal esxH epitope. Moreover, when we extended the analysis to two large genomic datasets of L2 and L4 strains, we found a much Figure 3. The hypervariable epitope at the N terminus of esxH (aa 1-18). The ancestral epitope is reported in top position, the derived haplotypes are reported below: mutations that were present in more than one lineage are highlighted in yellow, haplotypes that were found exclusively in L1, L2 or L4 are highlighted in pink, blue and red, respectively. Asterisks indicate the position inferred to be under positive selection by PAML: * posterior probability > 0.95: ** posterior probability > 0.99. The table on the right reports for each lineage the number of strains harboring the corresponding haplotype, and between parentheses, the number of independent parallel occurrences of the mutations as inferred by PAUP.
weaker signal: while 21% of L1 strains carried a derived haplotype for this epitope, only 1% of L2 and L4 strains, and no L3 strain had a derived haplotype. Despite analyzing datasets with more strains (6,752 and 10,466 for L2 and L4, respectively), and more polymorphic positions (140,309 and 277,648) compared to L1, we found only three amino acid replacements at the N-terminal epitope of esxH in L2, and seven in L4 ( Figure 3). The most frequently mutated position was the tenth amino acid, where we found 12 independent replacements in L1, and two in L2 and L4 ( Figure 3). The amino acid replacement A10V occurred eight times in parallel in L1, once in L2 and twice in L4. The most abundant derived haplotype was caused by a different amino acid replacement at position ten (A10T), which occurred once in L1 and once in L2 ( Figure 3). Overall, the replacements in all lineages occurred in eight residues in a peptide of 13 amino acids ( Figure 3).
We evaluated the robustness of these results by formally testing for positive selection on the complete sequence of esxH using PAML (Yang, 2007). We found that esxH was indeed under positive selection in L1 (p-value = 0.00004) but not in the other lineages (p-values = 0.39, 1.00 and 0.65 for L2, L3 and L4, respectively). PAML identified four codons that have been under positive selection in L1 (posterior probability > 0.95), all of them within the N-terminal epitope (codons 7, 9, 10 and 13; Figure 3). Codon 76, which is part of a different T cell epitope, had a posterior probability > 0.9, mutated three times in parallel and was possibly also under positive selection.
Our results further revealed that the derived haplotypes of the N-terminal esxH epitope were not distributed randomly across the geographic range of L1. Twenty-two of the 28 (79%) amino acid replacements occurred in sublineages L1.1.1 and L1.2.1, which are almost exclusively present in Southeast Asia (Extended data: Figure 11). We constructed a statistical test of association similar to phyC (Farhat et al., 2013) to determine whether replacements in the hypervariable esxH epitope were significantly associated with a particular geographic region (Methods). We found that South African strains were less likely to harbor a derived haplotype in the N-terminal esxH epitope than expected by chance (empirical p-value = 0.013; Table 1). While East African L1 strains were not associated with the derived haplotypes (empirical p-value = 0.276), we noticed important differences between countries: of the 29 East African strains harboring a derived haplotype, 28 were sampled in Madagascar. When we excluded Madagascar, we found that East Africa had a strong negative association with the derived haplotypes (i.e. East African strains harbored less derived haplotypes than expected by chance; empirical p-value = 0.0004). We then tested the most frequently replaced position (position ten) alone, and again found that East African strains were negatively associated with the derived haplotypes, with and without excluding Malagasy strains (empirical p-values = 0.0176 and 0.034, respectively). Finally, we tested the derived haplotype caused by the most frequent amino acid replacement (A10V; 8 parallel occurrences). Again, we found a negative association with East African strains (empirical p-values = 0.046 and 3 All p-values are one-tailed empirical p-values. This column indicates the direction of the association: + we tested for positive association between the derived haplotype(s) and the geographic region; -we tested for negative association between derived haplotype(s) and the geographic region. 4 Number of strains with the derived haplotype caused by the replacement A10V, and results of the test of association. 5 Number of strains with the derived haplotypes caused by any replacement at position 10 (A10X), and results of the test of association. 6 Number of strains with the derived haplotypes caused by any non-synonymous mutations occurred in the first 18 codons of esxH, and results of the test of association.
We hypothesized that the accumulation of missense mutations at the N-terminal epitope of esxH was due to immune escape. Therefore, we performed in silico prediction of the binding affinity of the ancestral haplotype and of the two most frequent derived haplotypes (caused by the amino acid replacements A10V and A10T) with different HLA-A, HLA-B and HLA-DRB1 alleles (Methods). We performed this analysis for: 1) Alleles found at high frequency (> 10%) in South-and Southeast Asia, but not in East Africa.
2) Alleles found at high frequency (> 10%) in East Africa, but not in South-and Southeast Asia.
However, we found no differences in the predicted binding affinities between alleles with different geographic distributions (Extended data: Table 8).
While esxH was the most striking example of a selective pressure specific for one lineage, our analysis suggests that it was not an isolated case. We performed a genome-wide scan for selection with PAML, and identified 17 genes under positive selection, of which five in common between L1 and L3 (Extended data: Table 9). We found two genes coding for transmembrane proteins, members of the Esx-1 secretion system, which were under positive selection in L1 (eccB1 and eccCa1, Bonferroni corrected p-values = 0.02 and 0.03), and several genes involved in antibiotic resistance that were under positive selection in both lineages (Extended data: Information). We further characterized the profile of drug resistance mutations of L1 and L3, and found that L1 strains harbored a greater proportion of inhA promoter mutations (conferring resistance to isoniazid) compared to L3 strains, confirming previous findings (Fenner et al., 2012; Extended data: Information, Figure 12, Table 10).

Discussion
Our results highlight the central role of South Asia in the dispersion of L1 and L3. First, we confirmed that the two lineages probably expanded from South Asia (O'Neill et al., 2019). Second, contrary to previous studies that assumed symmetric migration rates between regions (O'Neill et al., 2019), we found that these migrations occurred mostly in one direction. South Asia was a source of migrant strains that fueled the epidemics in other regions, especially in East Africa. Historically, a network of maritime trade, which followed the seasonality of the Monsoons, connected East Africa and South Asia. It is unclear how this could have promoted the spread of TB in one direction, but not in the other. A possible explanation is that strains originating in South Asia were more efficient in transmitting in East Africa, compared to East African strains that migrated to South Asia.
Another difference compared to previous studies is the temporal framework for the evolutionary history of L1. O'Neill et al.
(2019) estimated that the MRCA of L1 lived in the 4 th century BC and that the migration rates decreased after the 7 th century AD. Our results indicate that the MRCA of L1 did not exist before the 12 th century AD. This discrepancy is due to different assumptions about the clock rates of L1, but none of the two hypotheses can be excluded with the available data. Nonetheless, our discovery of a West African clade that likely originated through a direct introduction from Southeast Asia supports a faster clock rate of L1, and therefore a more recent MRCA (Extended data: Information We found evidence for a strong selective pressure acting on the N terminus of esxH in L1. In contrast, this selective pressure seems to be much reduced, if not absent, in the "modern" lineages (L2, L3 and L4). It is known that L1 strains interact differently with the infected hosts compared to other lineages. For example, L1 strains show a lower virulence in animal models (Bottai et al., 2020), transmit less efficiently in some clinical settings, and infect older patients (Holt et al., 2018). It has been shown recently that the increased virulence of the so-called "modern" MTBC strains is due to the loss of the genomic region TbD1, which remains present in L1 (Bottai et al., 2020). However, it was also reported that in some populations, L1 was associated with higher patient mortality (Smittipat et al., 2019). Given these differences with the "modern" lineages, it is likely that L1 is subject to different selective pressures, and it is possible that the greater selective pressure on esxH was caused by some epistatic effect specific to L1.
The signatures of selection at the N-terminal epitope of esxH were associated with strains sampled in South-and Southeast Asia, and were almost completely absent in East Africa (excluding Madagascar). Region-specific signatures of positive selection are a hallmark of local adaptation; in this case, adaptation of L1 strains to human hosts with South-and Southeast Asian genotypes. This corroborates previous studies reporting that in Taiwan, L1 is associated with indigenous populations with Austronesian ancestry (reviewed in Dou et al., 2015). This hypothesis is also supported by the L1 population in Madagascar. Madagascar is geographically linked to East Africa, however, Malagasies are genetically distinct from Africans, as they have mixed African and Southeast Asian ancestry due to the Austronesian colonization of Madagascar (Brucato et al., 2016). Madagascar was the region with the second highest prevalence of derived haplotypes in the N-terminal epitope of esxH after Southeast Asia (islands), 37.3% and 72.8% respectively, opposed to 0.3% (one single strain) in the rest of East Africa.
Although all the codons under positive selection in esxH were contained in one single T cell epitope, the selective pressure acting on esxH could be due to some other factor, and not to the binding of the epitope by the MHC, or to its recognition by T cells. A possible mechanism was suggested by experiments in a mouse model, showing that the N-terminal epitope of EsxH generates an immunodominant CD8 T cells response. The amino acid replacement A10T (which we found to be the most prevalent among the derived haplotypes of the N-terminal EsxH epitope) did not change MHC binding or T cell recognition, but it accelerated the degradation of the epitope, disrupting the immunodominant response (Sutiwisesak et al., 2020;Woodworth et al., 2008).
An additional driver of selection could be the effector function of EsxH. EsxH is a small effector secreted by the Esx-3 secretion system as dimer with EsxG (Ilghari et al., 2011). Within the host macrophages, the dimer EsxH-EsxG targets Hrs (a component of the human endosomal sorting complex), impairing trafficking, and hindering phagosomal maturation and antigen presentation, thus contributing to the survival of the pathogen (Mehra et al., 2013;Mittal et al., 2018;Portal-Celhay et al., 2016). The observed signatures of selection could be due to the adaptation of L1 strains to human genotypes of hrs prevalent in South-and Southeast Asia. A similar signature of selection was observed in another Esx effector (EsxW), in MTBC strains belonging to L2 (Holt et al., 2018). Holt and colleagues found evidence for parallel evolution at one residue in the N-terminal loop of EsxW, outside the region covered by known epitopes, suggesting that the selective pressure on EsxW was not due to antigen recognition.
The sampling effort in this study was considerable, and it provided a more complete picture compared to previous studies. Nevertheless, the sampling was not population-based, and for some regions the coverage was scarce (e.g. the Arabian Peninsula). Because of these limitations, our biogeographical analyses were limited to the subcontinental level. This approach revealed the global population structure and the main macroscopic patterns of diversity and migration of L1 and L3. However, MTBC populations are diverse also within subcontinental regions. For example, the MTBC population in Southern India is dominated by L1, while the most prevalent lineage in the North is L3 (Couvin et al., 2019). To investigate fine-scale processes in greater detail, including local adaptation, large population-based studies will be necessary.
In conclusion, the results presented here improve our knowledge about the TB epidemic around the Indian Ocean. A better understanding of the evolutionary dynamics of different MTBC populations might inform the development of control strategies for different regions. For example, esxH is part of several new vaccine candidates (Abel et al., 2010;Hoang et al., 2013;Radošević et al., 2007), and at least one of these, H4:IC31, is under clinical development in South Africa (Bekker et al., 2020;Nemes et al., 2018 were used to perform local realignment of reads around Indels. Reads with alignment score lower than (0.93*read_length)-(read_length*4*0.07)) were excluded: this corresponds to more than 7 miss-matches per 100 bp.
SNPs were called with Samtools v1.2 mpileup (Li et al., 2009) and VarScan v2.4.1 (Koboldt et al., 2012) using the following thresholds: minimum mapping quality of 20, minimum base quality at a position of 20, minimum read depth at a position of 7X, minimum percentage of reads supporting the call 90%. SNPs in previously defined repetitive regions were excluded (i.e. PPE and PE-PGRS genes, phages, insertion sequences and repeats longer than 50 bp) as described before (Brites et al., 2018).
We applied the following filters: genomes were excluded if they had 1) an average coverage <15x, 2) more than 50% of their SNPs excluded due to the strand bias filter, 3) more than 50% (or more than 1,000 in absolute number) of their SNPs having a percentage of reads supporting the call between 10% and 90%, 4) contained single nucleotide polymorphisms diagnostic for different MTBC lineages (Steiner et al., 2014), as this indicated that a mix of genomes was sequenced, 5) had more than 5,000 SNPs of difference compared to the reconstructed ancestral genome of the MTBC (Comas et al., 2013). Additionally, when multiple strains were sampled from the same patient, we retained only one. We further excluded all strains that had less SNPs than (average -(3 * standard deviation)) of the respective lineage (calculated after all previous filtering steps). We built SNP alignments for L1 and L3 separately, including only variable positions with less than 10% of missing data, and finally, we excluded all genomes with more than 10% of missing data in the alignment of the respective lineage. After all filtering steps, we were able to retrieve 4,968 strains with high quality genome sequences for further analyses (2,938 L1, 2,030 L3; Extended data: Table 1).

Analysis of sublineages
We used the curated datasets and inferred phylogenetic trees based on all polymorphic positions (excluding the ones in repetitive regions, see above) with raxml 8.2.11 (Stamatakis, 2014; -m GTRCAT and -V options). We then identified sublineages following the classification (and using the diagnostic SNPs) of Coll et al. (2014).

Molecular clock analyses with LSD
We selected all strains for which the year of sampling was known (2,499 strains, 1,672 L1, 827 L3). For both lineages, we built SNP alignments including only variable positions with less than 10% of missing data. We inferred phylogenetic trees with RAxML 8.2.11 (Stamatakis, 2014), using the GTR model (-m GTRCAT -V options). Since the alignments contain only variable positions, we rescaled the branch lengths of the trees: rescaled_branch_length = ((branch_length * alignment_lengths) / (alignment_length + invariant_sites)). We rooted the trees using the genome sequence of a L2 strain as outgroup (accession number SAMEA4441446).
We used the least square method implemented in LSD v0.3beta (To et al., 2016) to estimate the molecular clock rate with the QPD algorithm and calculating the confidence interval (options -f 100 and -s). We also performed a date randomization test by randomly reassigning the sampling dates among taxa 100 times and estimating the clock rate from the randomized and observed datasets. All LSD analyses were performed on the two lineages individually (L1 and L3), and on the five sublineages of L1.

Molecular clock analyses with Beast
Bayesian molecular clock analyses are computationally demanding and they would be impossible to apply onto the complete datasets. Therefore, we sub-sampled the L1 and L3 datasets used for the LSD analysis to 400 genomes with two different strategies: 1) random subsampling 2) random subsampling keeping at least 25 genomes for each year of sampling (where possible; "weighted subsampling"). For this second subsampling strategy, we used Treemmer v0.3 (Menardo et al., 2018) with options -X 400, -pr, -lm, and -mc 25. For both subsampling schemes, we generated three subsets of the original datasets, resulting in six subsets for each lineage. We assembled SNP alignments including only variable positions with less than 10% of missing data, and used jModel-Test 2.1.10 v20160303 (Darriba et al., 2012) to identify the best fitting nucleotide substitution model (according to the Akaike information criterion) among 11 possible schemes including unequal nucleotide frequencies (total models = 22, options -s 11 and -f).
We performed Bayesian inference with Beast 2.5 (Bouckaert et al., 2019). We corrected the xml files to specify the number of invariant sites as indicated here: https://groups.google.com/ forum/#!topic/beast-users/QfBHMOqImFE, and used the tip sampling year as calibration. We assumed the best fitting nucleotide substitution model as identified by jModelTest, a relaxed lognormal clock model (Drummond et al., 2006) and an exponential population size coalescent prior. We chose a 1/x prior for the population size [0-109], a 1/x prior for the mean of the lognormal distribution of the clock rate [10−10-10−5], a normal(0,1) prior for the standard deviation of the lognormal distribution of the clock rate [0 -infinity]. For the exponential growth rate prior, we used the standard Laplace distribution [-infinity-infinity]. For all datasets, we ran at least two runs, we used Tracer 1.7.1 (Rambaut et al., 2018) to identify and exclude the burn-in, to evaluate convergence among runs and to calculate the estimated sample size. We stopped the runs when at least two chains reached convergence, and the ESS of the posterior and of all parameters were larger than 200.
Since we detected a strong temporal signal only for L3, we performed a set of additional analyses of the subsets of L3. We repeated the Beast analysis with an extended Bayesian Skyline Plot (BSP) prior instead of the exponential growth prior, and performed a nested sampling analysis (Russel et al., 2019) to identify which of these two models (exponential growth and extended BSP) fitted the data best. The nested sampling was run with chainLength = 20000, particleCount= 4, and subChainLength = 10000.
All xml input files are available as Supplementary Files in Extended data.

Datasets for biogeography analyses
For the biogeography analysis, we considered only genome sequences obtained from strains for which the locality of sampling was known. When the country of sampling did not correspond to the country of origin of the patient (or was unknown), we considered as sampling locality the country of origin of the patient (this affected 187 strains, 121 L1 and 66 L3). Furthermore, similarly to other studies (O'Neill et al., 2019), we excluded all genomes that were sampled from Europe, North America and Australia, as most contemporary infections in these regions affect recent migrants. The final dataset for the biogeography analyses comprised 3,082 strains (2,061 L1 and 1,021 L3). We assigned the different isolates to different subcontinental geographic regions according to sampling locality. To do this, we followed a modified version of the United Nations geographic scheme (https://unstats.un.org/unsd/methodology/m49/; Extended data: Table 1 and Figures 1-3).

Phylogeography analysis with PASTML
For both lineages, we built SNP alignments including only variable positions with less than 10% of missing data, and all strains with known sampling locality (excluding strains from North America, Europe and Australia; see above, 2,061 L1 genomes and 1,021 L3 genomes). We inferred phylogenetic trees with raxml 8.2.11 (Stamatakis, 2014), using the GTR model (-m GTRCAT -V options). Since the alignments contain only variable positions, we re-scaled the branch lengths of the trees: rescaled_branch_length = ((branch_length * alignment_ lengths) / (alignment_length + invariant_sites)).
Since PASTML needs a time tree as input, we calibrated the phylogenies with LSD, assuming a clock rate of 1.4×10 -7 for L1, and 9×10 -8 for L3. In this analysis, genomes for which the sampling date was not known were assumed to have been sampled between 1995 and 2018, which is the period in which all strains with known date of isolation were sampled. Importantly, using different clock rates for this analysis would only change the time scale of the trees, but not the reconstruction of the ancestral characters.
We assigned to each strain the subcontinental geographic region of origin as character, and used PASTML (Ishikawa et al., 2019) to reconstruct the ancestral geographical ranges and their changes along the trees of L1 and L3. We used the MPPA as prediction method (standard settings) and added the character predicted by the joint reconstruction even if it was not selected by the Brier score (option -forced_joint). Additionally, we repeated the PASTML analysis for the sublineages of L1 individually.

Phylogeography analysis with Mascot
As a complementary method to reconstruct the ancestral range and the migration pattern of different populations, we used the Beast package Mascot (Müller et al., 2018). We assumed that strains sampled in the different subcontinental regions represent distinct subpopulations, and we considered only populations for which we had at least 75 genome sequences: four populations for L1 (East Africa, South Asia, Southeast Asia (mainland) and Southeast Asia (islands)), and two populations for L3 (East Africa and South Asia). For computational reasons, we subsampled the two datasets (L1 and L3) to ~300 strains. We sampled an equal number of strains from each geographic region (where possible), and within regions, we randomly sampled an equal number of strains from each country (where possible). This sub-sampling scheme resulted in two subsets of 303 and 300 strains for L1 and L3, respectively.
We assembled SNP alignments including only variable positions with less than 10% of missing data, and used jModelTest 2.1.10 v20160303 (Darriba et al., 2012) to identify the best fitting nucleotide substitution model as described above.
We performed Bayesian inference with Beast2.5 (Bouckaert et al., 2019). We corrected the xml files to specify the number of invariant sites as indicated here: https://groups.google.com/ forum/#!topic/beast-users/QfBHMOqImFE. We assumed a lognormal uncorrelated clock and we fixed the mean of the lognormal distribution of the clock rate to 1.4×10 -7 (L1) and 9×10 -8 (L3). We assigned the tip sampling years to the different strains, and when the sampling time was unknown, we assumed a uniform prior from 1995 to 2018 (similarly to what done in the PASTML analysis). We further assumed the best fitting nucleotide substitution model as identified by jModelTest, a gamma prior for the standard deviation of the lognormal distribution of the clock rate [0 -infinity], and a lognormal prior for the population size with standard deviation = 0.2, and mean estimated in real space. Finally, we used an exponential distribution with mean = 10 -4 as prior for the migration rates. For each analysis, we ran at least two runs. We used Tracer 1.7.1 (Rambaut et al., 2018) to identify and exclude the burn-in, to evaluate convergence among runs, and to calculate the estimated sample size. We stopped the runs when at least two chains reached convergence, and the ESS of the posterior, prior and of the parameters of interest (population sizes and migration rates) were larger than 200.
The xml input files are available as Supplementary Files in Extended data.

L2 and L4 datasets
For the analysis of esxH, we wanted to compare the results obtained for L1 and L3 with the other two major human-adapted MTBC lineages L2 and L4. Therefore, we compiled two datasets of publicly available genome sequences for these two lineages. We applied the same bioinformatic pipeline described above: genomes were excluded if they had 1) an average coverage <15x, 2) more than 50% (or more than 1,000 in absolute number) of their SNPs having a percentage of reads supporting the call between 10% and 90%, or 3) contained single nucleotide polymorphisms diagnostic for different MTBC lineages (Steiner et al., 2014), as this could indicate mixed infection. The final dataset consisted of 6,752 L2 genome sequences (with 140,309 polymorphic positions) and 10,466 L4 genome sequences (with 277,648 polymorphic positions; Extended data: Table 1). We reconstructed the phylogenetic tree of L2 with raxml 8.2.11 (Stamatakis, 2014) as described above. Due to the large size of the dataset, we used FastTree (Price et al., 2010) with options -nt -nocat -nosupport to reconstruct the phylogenetic tree of L4.

Epitopes analysis
We downloaded the amino acid sequence of all MTBC epitopes described for Homo sapiens from the immune epitope database (https://www.iedb.org/; downloaded on the 03.08.2020). We considered separately MHCI epitopes and MHCII epitopes, and we analyzed overlapping epitopes individually. We mapped the epitope sequences onto the H37Rv genome (GCF_000195955.2) using tblastn and excluded sequences mapping equally well to multiple loci in the H37Rv genome. Additionally, we retained only epitopes that mapped with two mismatches or less over the whole epitope length. This resulted in a final list of 539 MHCI epitopes, and 1,144 MHCII epitopes. (Extended data: Table 7). We used the datasets obtained for the biogeography analysis (2,061 genome sequences for L1, and 1,021 genome sequences for L3).
For each lineage, we independently assembled a multiple sequence alignment for each epitope. We then translated the sequence to amino acids and used PAUP 4.a (Wilgenbusch & Swofford, 2003) to reconstruct the replacement history of all polymorphic positions on the rooted phylogenetic trees. We used two maximum parsimony algorithms (ACCTRAN and DELTRAN) and considered only the events reconstructed by both algorithms.

Analysis of esxH
We considered the first 18 amino acids of esxH (MSQIMY-NYPAMLGHAGDM), which was by far the epitope with most amino acid replacements in L1. We expanded the analysis of this epitope to the L2 and L4 datasets, so that we could compare the results with L1 and L3. For the PAML analysis, we randomly selected 500 strains from each MTBC lineage. We used the phylogenetic tree reconstructed by RAxML (same settings as above), and the gene alignment to estimate the branch lengths of the tree using the M0 codon model implemented in PAML 4.9e (Yang, 2007). This step was necessary to obtain a tree with the branch length in expected substitutions per codon. We then fitted two alternative codon models (M1a and M2a) to the trees and alignments. M1a allows ω to be variable across sites, and it assumes two different ω (0 < ω 0 < 1, and ω 1 = 1), modeling nearly neutral evolution; M2a assumes one additional ω (ω 2 > 1) compared to M1a, thus modeling positive selection. We performed a likelihood ratio test between the two models as described in Zhang et al. (2005). Templates for the control files of the codeml analyses (M0, M1a, M2a) are available as Supplementary Files in Extended data. The codon under positive selection were identified with the Bayes empirical Bayes method (Yang et al., 2005).
To test whether the derived haplotypes of this epitope were associated with a specific geographic region, we constructed a statistical test similar (but not identical) to PhyC (Farhat et al., 2013), which is normally used to test for association between a variant and phenotypic drug resistance. For each aminoacid replacement (or combination of aminoacid replacements): we considered the number of independent mutational events that occurred along the tree of L1 (as estimated by PAUP). We then randomly redistributed the same number of replacements on the phylogenetic tree of L1 (based on the branch length; i.e. on longer branches there is a higher probability of a mutational event), and counted the number of strains with a derived state for each geographic region. We built null distributions by repeating this procedure 10,000 times. Under the null hypothesis, mutations occur randomly on the tree, and therefore independently from the geographic region where the tips were sampled. For each region, we ranked the observed value (the number of tips with a derived state) on the corresponding null distribution to obtain empirical p-values. We used the same geographic region used for the biogeography analysis. Additionally, we considered East Africa excluding Madagascar.
R code and input files to perform this test are available as Supplementary Files in Extended data.

Prediction of binding affinity between T cell epitopes and HLA alleles
We considered three HLA loci: HLA-A, HLA-B and HLA-DRB1. We used the allele frequency database (Gonzalez-Galarza et al., 2020) to identify alleles that are prevalent in East Africa and not in South Asia and Southeast Asia, or the other way around. Because the coverage of the allele frequency database is patchy, we focused on the following countries: Kenya and Zimbabwe as representatives of East Africa; India, Thailand and Taiwan as representatives of South-and Southeast Asia. We identified: 1) Alleles that had a frequency of 10% or more in at least one population in South-and Southeast Asia, but had frequencies lower than 10% in all East African populations.
2) Alleles that had a frequency of 10% or more in at least one population in East Africa, but had frequencies lower than 10% in all populations in South-and Southeast Asia.
3) Alleles that had a frequency of 10% or more in at least one population both in South-and Southeast Asia and in East Africa.
For all these alleles, we performed in silico binding prediction with three epitopes: the ancestral epitope at the esxH N terminus (MSQIMYNYPAMLGHAGDM), and the two most frequently observed derived epitopes (MSQIMYNYPTMLGH-AGDM, and MSQIMYNYPVMLGHAGDM). For HLA-A and HLA-B alleles, we used the NetMHCPan4.1 server (Reynisson et al., 2020) with standard settings. For HLA-DRB1 alleles, we used the prediction tool of the immune epitope database (https://www.iedb.org/).

Genome wide scan for positive selection with PAML
For this analysis, we used the subsets generated for the Mascot analysis. These datasets are representative of the populations of L1 and L3 in their core geographic ranges, and are computationally treatable. We generated sequence alignments for all genes individually, excluding genes in repetitive regions of the genome (see above). Because some genes are deleted in L1 but not in L3, or the other way around, we obtained a slightly different number of gene alignments for the two lineages (L1: 3,623, L3: 3,622). For each gene, we performed a test for positive selection with PAML as described above for esxH. We considered as under positive selection all genes, for which the likelihood ratio test resulted in a Bonferroni-corrected p-value < 0.05.

Drug resistance mutations profiles
We considered 196 mutations conferring resistance to different antibiotics (Payne et al., 2019). We extracted the respective genomic positions from the vcf file of the 4,968 genomes of the complete curated data set (2,938 L1, 2,030 L3) and assembled them in phylip format. To determine the number of independent mutations, we reconstructed the nucleotide changes on the phylogenetic tree rooted with the L2 strain (SAMEA4441446).
To do this, we used the maximum parsimony ACCTRAN and DELTRAN algorithms implemented in PAUP 4.0a (Wilgenbusch & Swofford, 2003), and considered only the events reconstructed by both algorithms. • The folder Dating_Beast contains the xml files for the dating analyses.

Data availability
• The folder EsxH_haplotype_association_test contains the code and input files to perform the test of association between different haplotypes and geographic regions.
• The folder EsxH_PAML contains the control files and input files to perform the positive selection analysis on EsxH.
• The folder Mascot contains the xml files for the Mascot analyses. In the paper "Local adaptation in populations of Mycobacterium tuberculosis endemic to the Indian Ocean Rim", Menardo and colleagues explore diversity, phylogeographic relationships and patterns of adaptation in TB lineages 1 and 3, both of which are distributed mainly around the Indian Ocean. The paper demonstrates a strongly asymmetric trend of South Asian spread to other regions. Further, it describes a search for positive selection throughout these genomes, and identify one gene in particular, exsH, as being a target for T-cell dependent selection. The authors additionally calibrates molecular clocks for both L1 and L3 and use this to date certain events in the tree.
The methodology is solid and well suited to answer the questions posed. I found the paper to be well written and helpful in shedding light on important questions related to the evolution of L1 and L3. I have only a few minor discussion points. Now, the clock analyses are not heavily emphasized, but I thought it might be a potential flaw that the calibration rests on the assumption that the foundation of a single sub-lineage in West Africa could not have happened prior to the establishment of Portuguese sailing routes around the tip of Africa. While I'm certainly not an expert on historic migrations or anthropology, I would assume there must be many possible alternative scenarios. (For example, pilgrimages across the Sahara are known to have happened much prior to this). As far as I can tell, this calibration is the major driving force for the difference in L1 TMRCA from O'Neill. It would be interesting to know how the degree to which this calibration affects the TMRCA and clock rate. On a somewhat related note, this article fails to cite some key references on the phylogeography of TB.
In one section the authors evaluate association between esxH haplotypes and the geographic region of origin. They calculate this association using a home-made statistical test based on phyC. I could not gather from the description exactly what they measure. It is a simulation-based test that measures "how many strains from each region resulted to have a derived state according to the simulation". Is there a test statistic being ranked here? Is it single-region vs the rest or is the statistic a composite measure of multiple regions? How is the simulation for the null distribution set up? E.g. is it a forward-in-time simulation using a specified mutation rate or a permutation or something else? I see that the code for running this test is posted on github -However, I still think a few more details in the actual paper is in order for the biologists that don't read R. Speaking of this test, It is a bit unconvential for Table 1 to highlight all results with P<0.1, especially since a large number of associations are being tested here. I would not consider p<0.1 to be particularly strong evidence of association.
I'd like to commend the authors for making so much of their data publicly available. Datasets are in immediately-usable formats which makes it simple to download all the raw data, and as mentioned code is also readily available.
Very minor: HLA abbreviation never explained in manuscript.

○
Tip date randomisation performed, but the results of this test is never stated.
○ Why is PRJNA630228 posted? As far as I can tell it doesnt contain any SRA data, and the other accession contains data for all 767 strains.

If applicable, is the statistical analysis and its interpretation appropriate? Yes
Are all the source data underlying the results available to ensure full reproducibility? Yes

Are the conclusions drawn adequately supported by the results? Yes
Competing Interests: No competing interests were disclosed. Thank you very much for your comments, and for reviewing our manuscript. We address all your points below. We prepared a revised version of the manuscript based on the feedback of all reviewers.
Regarding the clock analysis of L1: We performed a tip dating analysis but found no temporal signal for L1. Therefore we used the clock rate estimate resulted from a previous study (Menardo et al. 2019) where we analyzed a different L1 dataset and found a good temporal signal and a clock rate of ~1.4x10 -7 . So far Menardo et al. 2019 is the only published study that estimated the clock rate for a L1 dataset.
In a second step we tested whether this clock rate and two alternative values, 5x10 -8 and 1.2x10 -7, , (the lowest one was used by O'Neill et al. 2019) were compatible with the hypothesis of a West African introduction of L1.1.1 posterior to the Portuguese circumnavigation of Africa. We found that only the largest clock rate was compatible with this hypothesis.
We are quite explicit in the manuscript about the limitation of this analysis, and how it should be interpreted in comparison to the results of O'Neill et al. 2019: "This discrepancy is due to different assumptions about the clock rates of L1, but none of the two hypotheses can be excluded with the available data".
Finally, it is possible that the West African L1.1.1 clade did not originate through a direct migration from South East Asia. However, a series of short range migrations seems less likely, as we did not find strains belonging to L1.1.1 in South and Central Asia, or in East Africa.
As we write in the manuscript, the dating analyses for MTB in general, and for this L1 dataset in particular, have limitations.
Here we discussed what we consider the most likely scenario given the data, we did not exclude alternative scenarios.
We added a section in the discussion comparing our findings with the biogeography of other MTBC lineages, and referencing the relevant studies.
The results of the DRT are reported in the supplementary Information (now amended to enhance clarity), in Sup. Figs 8 and 9 and in Sup. Table 3. In the main text we refer more generally to the lack of temporal signal (or to a good temporal signal), as this reflects the results of the DRT, but also of the Beast analyses.
We now refer only to PRJNA670836, which contains data for all the newly sequenced strains. Thanks.
complex (MTBC) strains belonging to Lineages 1 and 3, which are the most prevalent lineages in South Asian countries with a high TB burden, such as India. The authors compiled an excellent dataset of nearly 5,000 MTBC genomes from these two lineages, comprising publicly available as well as newly-sequenced MTBC strains. Their phylogenomic analyses showed that L1 comprises a further five distinct sublineages, whereas L3 does not contain distinct sublineages. The authors used two different methods (PASTML and Mascot) for conducting phylogeographic analyses, with results from both methods suggesting an interesting asymmetric pattern of migration: these lineages originated in and migrated from South Asia to other regions, with very little backmigration into South Asia. This underscores the important role South Asia has played in the dispersal of these lineages. The authors found several genes involved in virulence and antibiotic resistance to be under positive selection in L1 and L3 strains. Lastly, the authors also found evidence of local adaptation at the esxH epitope in L1 strains from South and Southeast Asia, which requires further study as this epitope is a part of several vaccine candidates.
Overall, I found this to be a well-designed and thorough study, and a very well-written manuscript. I have no major concerns and only a few minor suggestions for aesthetic improvement: The Methods section titled "Strain cultures, DNA extraction and genome sequencing" is missing citations for the methods used (eg. CTAB extraction method). expertise to confirm that it is of an acceptable scientific standard.

Author Response 19 Mar 2021
Fabrizio Menardo, Swiss Tropical and Public Health Institute, Basel, Switzerland Thank you very much for your comments, and for reviewing our manuscript. We prepared a revised version of the manuscript based on the feedback of all reviewers.
We added the missing reference for the DNA extraction method, and we revised the supplementary material for errors and typos, thanks.
Regarding Figures 1A and 2A: in both figures the other three panels have the same color code, where different colors correspond to different world regions. We decided to use a gray scale to avoid potential confusion of having two different color codes in the same figure. We did not modify the figures.