Summary
By gathering an initial collection of 680 public Sequence Read Archives from the Mycobacterium tuberculosis complex (MTC) including 190 belonging to the lineage 2 Beijing, and using a new bioinformatical pipeline, TB-Annotator, that analyzes more than 50,000 characters, we characterized a phylogenetically significant L2 sublineage from isolates found in Tochigi province, that we designate as Asia Ancestral 5. These isolates harbor a number of specific criteria (42 specific SNPs) and their intra-cluster pairwise distance suggested historical and not epidemiological transmission. These isolates harbor a mutation in rpoC, and do not fulfill, either the Modern Beijing lineage criteria (L2.2.1.2.2) or the other currently known Ancestral Beijing lineages described so far. Asia Ancestral 5 isolates do not possess mutT2 58 and ogt 12 characteristics of Modern Beijing, but possess ancient evolutionary characteristics. By looking into the literature, we found in a reference isolate ID381, found in Kobe and Osaka, and defined as “G3”, 36 out of 42 shared and specific SNPs in the same sublineage. Using this study, we also confirmed the intermediate position of the Asia Ancestral 4 recently described in Thailand and suggest an improved classification of the L2 that now includes Asia Ancestral 4 and 5. Increasing the recruitment to around 3000 genomes (including 642 belonging to L2) confirmed our results and suggest new historical ancestral L2 phylogenetically relevant branches that remain to be investigated in detail. We discuss some anthropological and historical data from Japan history and its link to Korea and China. This study shows that the reconstruction of the early history of tuberculosis pandemia in Asia is likely to reveal complex patterns since its emergence.
Introduction
With 10 millions new cases in 2020, and around 500,000 Multi-drug resistant cases, Tuberculosis (TB) is far from being eradicated [1]. There are nine recognized genetic lineages in the bacterial pathogen agent, Mycobacterium tuberculosis complex (MTC), and one of these, lineage 2 (L2), is of particular interest due to the emergence of multi-drug-resistance tuberculosis (MDR-TB) and extremely drug-resistant tuberculosis (XDR-TB), that were described in numerous studies run in Russia, Europe, Africa, Asia and America [2-6]. L2 finds its origin in China and is predominant in Asia [7, 8]. The adaptation ability of these bacilli to develop specific pathogenic features including compensatory mutations in genes that restore their fitness have largely contributed to their epidemic success and to the global spreading of multidrug resistant strains [6, 9-11].
The epidemic bursts of L2 in the 90s can be linked to contemporary historical and geopolitical events, such as: (1) the fall of the former USSR, (2) the poor individualized TB treatment and follow-up of prisoners in Former Soviet Union (FSU) republics, that likely contributed when freed, to the seeding of the general population in Russia with L2 isolates that had become multidrug resistant in prisoners populations, (3) the increasing global trade share of China in the world, and in particular in Middle-East, South and East-Africa, and more generally around the African coasts, (4) the current Covid-19 pandemia with its dramatic consequences on global health [12-15].
The discovery of the Beijing Lineage can be dated back to 1995, even if the existence of the “W” strain, a member of L2, was known since its discovery in New York in the early 90s [16, 17]. It was first characterized by its high IS6110 Insertion Sequence copy number (16 to 20 copies), with a diagnostic copy in the dnaA-dnaN also designated as the NTF region [18, 19]. No lineage such as L2 has been fostering so many genetic, genomic, molecular epidemiological and public health studies, whether for detection of multidrug resistance (a feature rapidly shown to be associated to L2), for improvement of its taxonomy through deletion region analysis, for a better understanding of its genetic population structure, and of course for its surveillance [4, 20, 21]. This was also made possible thanks to Variable Number of Tandem Repeats (VNTR) analysis, by Regions of Difference (RD) analysis using DNA hybridization, and lastly by the Whole-Genome Sequencing approach that led to single nucleotide polymorphisms (SNP) discoveries [22, 23]. Very large outbreaks inside this lineage were shown to have independently emerged in Russia, in Europe, in Central Asia, in South Africa, and of course in China [24].
L2 were previously classified as “typical” (modern) or “atypical” (ancestral/ancient) based on IS6110-RFLP patterns and later on Deletion Region analysis [25, 26]. Today, the L2 lineage is split into two main sublineages, L2.1 and L2.2. L2.1 is known as “Proto-Beijing” and is found mainly in Southern China, particularly in Guangxi province [27]. Among the major characteristics is the fact that it possesses the L2 characterizing RD105 deletion but not RD207 [28]. Recently a rare L2.1 clone was shown to have become ultra-drug resistant [28]. Such findings raise the issue of genetic background on adaptation capacity and acquisition of multi-drug resistance. L2.2, or the “Beijing” strain, as defined by spoligotyping SIT1 most of the time, is composed of several sublineages broadly categorized into the “Ancestral” and “Modern” Beijing strains. L2.2.2 defines the Asia Ancestral 1 group, whereas L2.2.1 defines all others, and that do not fit anymore with the “modern Beijing” designation [29]. The switch from “Ancestral” to “Modern” Beijing would be associated to the presence of at least one IS6110 copy in the NTF and the presence of the MutT2 G->C mutation in position 1286766 and ogt codon 12 (C1477596T, ogt12) [9, 30, 31]. The “Ancestral Beijing” strains form a cluster made up, until now, of three “Asia Ancestral” lineages (Asia Ancestral 1 to Asia Ancestral 3, AAnc1 to AAnc3) [31], to which a fourth “Asia Ancestral 4” (AAnc4), found in Northern Thailand, was recently added [32]. RD181 is not deleted in L2.1 and L2.2.2 whereas is deleted in all other sublineages [31].
The “Modern” Beijing represents the clade that is at the end of the phylogenetic tree and it encompasses at least 7 large and epidemiologically acknowledged sublineages (Asian African 1, Asian African 3, Europe/Russia B0/W148 outbreak, Central Asia, Asian African 2, Asian African2/RD142, Pacific RD150) [31]. Interestingly, apart from the SNPs discovered these recent years, other major genomic rearrangements were shown to have happened in some L2 genomes [33, 34]. Modern Beijing strains are responsible of most but not of all of the recent MDR-TB outbreaks and are more transmissible [4, 35, 36].
However, apart from these well MDR-TB characterized strains and lineages, there remain many L2 isolates that do not fit to any known sublineages [37]. As an example, 26 Clonal Complexes (Bmyc1 to Bmyc26) were described based on 3R gene polymorphisms, however some of these types (e.g. Bmyc6 or Bmyc26) remained unclusterized based on the most recent classification scheme [31]. This suggests that some improvements in the Shitikov L2 classification scheme are possible; meanwhile other L2 sublineages were described such as in Thailand or Japan, and their history should also be linked to the global L2 history by performing comparative genomics [6, 23, 28, 32, 38-42].
The difficulties to define L2 epidemiological clusters using IS6110-RFLP, the impossibility to study this lineage by spoligotyping (almost all strains harbor the SIT1 spoligotype with only spacers 34-43 left) promoted the use of deletion regions and MIRU-VNTR to decipher its complex population structure [43]. One example of these difficulties was the need to add hypervariable VNTR Loci for L2 molecular epidemiological studies [44, 45]. A 15-9 VNTR signature designation was also introduced to facilitate designation of MTC epidemiological clones [45]. Hence, the introduction of Whole Genome Sequencing (WGS) and the building of very large databases to create representative worldwide collections is the best way to develop a more comprehensive and precise knowledge of this complex lineage [4, 23, 46]. The Merker et al. study represents the greatest endeavour to date to get an in-depth characterization of the L2 lineage [23]. It showed the existence of 6 clonal complex (CC1 to CC6) and a basal sublineage (BL7) [23].
RD and SNP-based phylogeny provides the best tools to classify L2 isolates. However, an even more refined evolutionary scheme, can be based on IS6110 insertion sites that allows a fine study on microevolution of all lineages [47, 48]. Indeed, the adaptation to an IS6110 high copy number and diversity of insertion sites in L2 are likely to be linked to increased mutation rates and are likely to provide in some cases a higher fitness of some strains [36, 49]. The use of MIRU-VNTR results in association to WGS is also interesting but seems to reach some limits [39]. Recently, thanks to a very useful standardization of all previously identified sublineages naming, a very clear population structure of L2 emerged, with 11 main L2 sublineages [31].
The history of the evolution of L2 still contains many local and global shadow zones, in particular concerning the dating of the split between L2 and L3, and the precise dating of emergence and geographical origin of the various sublineages and on their specific phylodynamics [23, 30]. Hirsch et al. suggested that East Asian and Philippines human populations carrying M. tuberculosis belonging to distinct lineages may have split 240–1,000 years ago [50]. Merker et al. generated estimates for the time to the most recent common ancestor with a TMRCAs of ∼6,000 and 5,000 years for CC6 and BL7 (the most ancient estimated by coalescence), and a TMRCA of ∼1,500 years for CC5 (the most recently appeared) [23]. More recently however, different estimated emergence dates were given [8]. In this study, an estimate of 2200 Before Present (BP) was given for the coalescence between Proto-Beijing and ancestral Beijing, and an estimate of 1300 BP for the split between all “ancestral-Beijing” lineages [8]. 900 BP would correspond to “Proto-Beijing” expansion whereas expansion of “Modern Beijing” would have appeared only around 500 years ago. The place of Japan, and its relation to China was also raised by some authors, however no precise emergence date was suggested [51]. Hence, as seen above, there remain many uncertainties on the exact dating of origin, and various split and bursts of L2.
The geographic origin of the emergence of L2 is almost as blurry as its dating. According to some authors, north-central and northeast China could have been the initial spreading center [51]. According to others, based on large differences in the prevalence of the ancient L2 lineages in China and on the higher genetic diversity observed in the south-west province of Guizhou, South-China would be the geographical origin of L2 [27, 30]. Indeed, Guizhou counts 17 ethnical minorities, and most of the non-officially acknowledged ethnic groups of China are located in this province. South-East-Asia shows a higher human genetic diversity than the north-east-Asia [52]. Arguing in favour of a South-China origin, a recent paper described the existence of an “Asia-Ancestral 4” L2 sublineage, that was found in the north of Thailand in a city named Chiangrai [32]. Chiangrai has been inhabited since the 7th century and peopled by tribes now residing in Thailand, Myanmar and Lao [32]. In this place, Ancestral Beijing strains are associated with ethnic minorities, Akha and Lahu tribes, who mostly migrated from Southern China [32].
In Japan, one of the main characteristics of the tuberculosis global history is the presence of a high and still poorly characterized genetic diversity of Ancestral L2 Beijing strains [23,42,44]. Mokrousov had detected that MIRU-VNTR diversity was quite important and specific in Japan [51]. Indeed, contrarily to what is observed in many other countries, L2 outbreaks in Japan are often due to ancestral Beijing strains, especially in aged people, that were not vaccinated by BCG [39, 42, 44, 53]. Hence, tuberculosis in Japan could be as ancient as in China or almost. Indeed, peopling of Japan is known to have originated through the Korean peninsula very early in the history of human migrations, followed by what is known as the Jomon period [54]. During the following Yayoi period, technology and agriculture was imported from Korea and from China [54].
We are developing the “TB-Annotator” project, a new bioinformatical pipeline whose aim is to perform data-mining of MTBC genomic diversity for both evolutionary biology and epidemiology using most of the information contained in Sequence Read Archives (SRAs). During the advancement of this long-term research project, we studied a set of L2 isolates from a prefecture located in Central Honshu in Japan [42]. In this study, 169 clinical isolates of M. tuberculosis were obtained from as many TB patients residing in Tochigi Prefecture [42]. Patient strains were collected in 2007 or in 2013 and were more likely to belong to the Ancestral (Atypical) L2 Lineage (n=79) than to the Modern (typical) (n=32) (see also Table 1 in [42]). Three patients only, out of the 79 bearing an Atypical Beijing strains were not born in Japan [42]. Japanese patients suffering from an Atypical Beijing were more likely to be found in the Central part of the Tochigi prefecture [42]. All the isolates had been characterized by WGS and analyzed using a bioinformatical pipeline, coined as “CAST” [55], however an in-depth comparative study of these isolates with other ancestral L2 had not been performed and we thus decided to study these genomes in more detail using our pipeline.
As more and more MTC genomes from isolated population are described, i.e. mountainous areas or isolated ethnic tribes or islanders, it becomes more and more likely to find rare though interesting L2 isolates that have their own specific history [32]. We describe in this article such a finding, i.e. the characterization of a new L2 sublineage from Japan and that is exclusive from Japan for the time-being. We also confirm the “intermediate” (i.e. between ancestral and modern) character of the AAnc4 sublineage recently described at the light of this comparative genomic analyses, and link these two lineages to the latest public and unified L2 classification scheme [31].
Material and Methods
Selected genomes
We downloaded a set of 680 SRAs; these genomes were selected to include the maximum of diversity of the currently described MTC sublineages characterized so far, including recent landmark papers [56, 57]. The list of these SRAs is shown in Supplementary Table1 and includes 190 L2-assigned SRAs of all sublineages except for the terminal Pacific RD150 sublineage. From [32], we selected 28 SRAs labeled as “Asia ancestral 4”. From [42], we initially included a total of 158 representative genomes, however we discarded the SRAs for which the coverage was either too weak or for which it was impossible to rebuild the spoligotype using CRISPR-builder-TB [58, 59].
Brief description of the TB-Annotator pipeline
The full version of TB-Annotator is going to be released in another article. In brief, regarding the set-up of the platform, the processes are the following: SRAs of interest are selected and kept only provided a certain number of conditions that together reinforce the reliability of the data: they have read length>75bp, if clean reads file is at least 100 Mo, and if CRISPR could be reconstructed using CRISPRbuilder-TB [59]. For each SRA, in addition to reconstructing the CRISPR-Cas region using CRISPRbuilder-TB and apart from collecting NCBI information on the genomes, two scripts successively: (1) search for SNPs according to reference catalogues (Supplementary Table2) totaling more than 50,000 SNPs, including drug-resistance related SNPs, phylogenetic SNPs as per 26 studies contributing to SNP-based classification ([23, 29, 31, 56, 57, 60-80]; (2) search for additional SNPs in each new isolate based on H37Rv reference sequence; (3) look for the presence/absence of H37Rv genes as annotated in mycobrowser (https://mycobrowser.epfl.ch/), (4) look for the presence/absence of deletion regions; (5) identify insertion sites of all known Insertion sequences in MTBC. The CRISPR locus is rebuilt using a dedicated and previously published script (43/68/360 Spacers format) with assignment of a Spoligo-International-Type (SIT) tag; the application produces an ordered list of spacers/repeat with variants and IS6110 insertion sequences if present [58, 59, 81]. Final scripts allow to produce a phylogeny based on the list of studied characters and using RAxML and SplitsTree [82, 83].
All computations have been performed on the “Mésocentre de Franche-Comté” supercomputer facilities (141 nodes, 2292 cores, 9,27 To memory, 74,2 CPU Power TFlops, 66,4 GPU TFlops), using adequate commands. Apart from the Excel results file, the final phylogenetic tree on 680 selected genomes was built using SplitsTree to provide a graphical display of the results and using proprietary python scripts allowing queries to be made [83].
SNP-based classification of L2 sublineages
In order to first assign the 680 selected SRAs into known lineages and sublineages, we used the reference list defined in Supplementary Table2. Full origin and classification results on the 190 SRA of the L2 isolates, as produced by “TB-Annotator” are found in the Supplementary Table3.
Definition of an unclassified Asia Ancestral 5 sublineage (AAnc5)
After assessment of existing known L2 sublineages, the unknown branch (Figure 1B) was investigated in detail and was shown to originate only from Japan SRAs. The bioinformatical pipeline, based on its graphical user interface, allows to select and display new exclusive or shared SNPs and specific characteristics, that were further investigated.
Left part (1A); TB-annotator unrooted phylogenetic tree on 680 SRA-derived data. L2 samples are shown in blue. Right part (1B); close-up on the Lineage 2 with all known branches named except in red the new unknown ancestral Japan sublineage we design as Asia Ancestral 5.
Use of TB-annotator platform to explore diversity
The TB-Annotator interface is meant to explore diversity among TB samples. Presently, 680 samples representative of TB diversity were implemented. Other samples can be added upon request. The default layout presents the samples phylogeny. One specific function allows to identify the genetic markers that are specific to selected samples, either based on a graphic mode or on predefined filters. Useful metrics such as intra-cluster pairwise SNP distance can be computed. A second layout allows to identify the countries from which the selected samples were collected on a map. At the time of writing a new version with 3380 genomes confirmed our results.
Results
1. In depth analysis of known SNPs phylogenetical markers within the 190 L2 genomes, discovery of a new L2 Sublineage
We implemented a representative set of samples from L2 and the 169 samples from the Tochigi province study, Honshu, onto the web platform TB-annotator [42]. Regarding the presence of phylogenetic SNPs and clustering in the phylogeny, as expected, all samples carrying the L2 SNP were found in the same phylogenetic branch (Figure 1), and 186 of the 190 isolates inside this branch carried this SNP (4 L2.1 isolates were missed). This branch can thus be considered as L2. Among these L2 samples, nine were Proto-Beijing (L2.1), ten were AAnc1 1 (L2.2.2) and the 171 isolates left were confirmed to fit the former L2.2.1 definition using Coll’s classification scheme [29]. This L2.2.1 sublineage, as previously mentioned, should not be referred to as “modern” Beijing since two of its sub-branches, AAnc 2 and 3 have been defined since [31]. These 171 L2.2.1 isolates shared a list of 31 exclusive SNPs (results not shown).
Going further in the analysis, 93 isolates from L2.2.1 were harboring the MutT2 G1286766C “modern Beijing” characteristic SNP. Out of these 93, 65 samples were also harboring the second “modern Beijing” signature, i.e. C1477596T SNP in ogt [32] while the 28 remaining harbored Asia Ancestral SNPs.
The “modern” Beijing (n=65), could be split into Asian African 1 (n=2), Asian African 3 (n=5), Asian African 2 (AAfr2, n=16), Asian African 2-RD142 (AAfr2-RD142, n=1), Central Asia (n=13), B0/W148 (n=2) and there were a remaining 26 “unclassified modern Beijing” isolates, that did not fit any described modern sublineage definition, and that were not investigated further in this study. No Pacific RD150 isolates was included in this study. All suggested binary names fit with Shitikov et al. scheme and we present a new improved classification that includes the recent discovery of “Asia Ancestral 4 and 5” (Figure 2).
Mtb lineage 2 dendrogram representing the major phylogenetic groups and informative genetic markers (inspired and redrawn from Shitikov et al. 2017; note that on this new figure, SNP position of MutT2 and ogt 12 were corrected compared to the original figure where position were inverted). Bmyc6 and Bmyc26 were not considered in the binary naming suggested here.
As shown on Figure 1, the “ancestral” Beijing (n=124) could be split into Proto-Beijing (L2.1; n=9), Asia Ancestral 1 (L2.2.2-AAnc1, n=10), Asia-Ancestral 2 (L2.2.1.1-AAnc2, n=22), Asia-ancestral 3 (L2.2.1.2.1-AAnc3, n=9), “Asia-Ancestral 4” (L2.2.1.2.3*-AAnc4, n=27) and an unknown branch (suggested L2.2.1.3*-Asia Ancestral 5 or AAnc5, n=20); there still remains 27 unclustered and uncategorized ancestral Beijing.
The reanalysis of the phylogenetical SNP set described by Shitikov et al. confirmed all the phylogenetical SNPs specific of the Ancestral Beijing lineages up to the definition of the “modern” Beijing. The MutT2 SNP (G1286766C) is a good phylogenetic marker as it is present only in L2 modern sublineages whereas T1892017C, C4137829T and C4393590G (MutT4) are present in some ancestral sublineages. It is clear from the SNPs results and from the phylogenetical trees shown in Figure 1 that the AAnc4 branch is intermediate, as suggested by their inventors, and is neither a bona fide modern L2 isolate nor a truly ancestral lineage [32]. Conversely, the 20 isolates that we studied from Japan, further designated as Asia Ancestral 5 (AAnc5), all fulfill the ancestral Beijing criteria since they do not possess the expected MutT2 (58) SNP nor any other characteristic of modern Beijing. They also branch before AAnc4. In all cases, the clear-cut split of L2 isolates that either possess none, one, or the two MutT2 (58) and ogt (12) SNPs confirms that these two markers are excellent phylogenetical ones. Based on the results obtained, we decided to adapt the classification scheme of Shitikov et al. to take into account all recently described sublineages (Figure 2).
2. Description of the new “Asia-Ancestral 5” AAnc 5 Japan sublineage
2.1 SNP-based results
We further characterized the Japan “Asia-Ancestral 5” (AAnc5) (n=20) that branches between AAnc 1 and 3 (Figure 1). These isolates can all be defined by many exclusively shared SNPs. The complete results matrix produced by the TB-Annotator pipeline on the 190 L2 genomes can be found in Supplementary Table3. In the future versions, reports will be downloadable from a dedicated website. The list of 42 exclusive SNPs, found in 16 of these 20 isolates is shown in Supplementary Table4. An interesting observation is that all the AAnc5 strains harbor a non-synonymous mutation in 765140 (G->C) in rpoC. This results will be discussed in detail further. Three other genomes (2 in other L2 sublineages and one in a L4 sublineage) also carry this rpoC mutation suggesting independent acquisition.
This new phylogenetical branch of L2 is believed to have emerged after the birth of AAnc1 but before the ones of AAnc2 and 3. A pairwise distance matrix between these isolates was constructed (Supplementary Table5); the future pipeline will systematically compute the intra-branch SNP distance for clusters or selections of interest below 50 SRAs. The pairwise distance between these samples shows between a minimum of 166 SNPs (between DRR157280 and DRR157281) and a maximum of 439 SNPs distance (between DRR130203 and DRR034366) between AAnc5 isolates, thus excluding recent transmission (Supplementary Table 5). Assuming a 0.3 SNP mutation rate per year per genome, these strains might have diverged approximately between 250 to 600 years from their MRCA. Our results also confirmed the intermediate character of the “Asia Ancestral 4” (AAnc4) which was described in Chiangrai in Thailand [32]. If we accept the suggestion of timing of AAnc4 emergence or expansion around the 7th century in Thailand, (start of Chiangrai), then AAnc5 could have been introduced into Japan earlier than this century, in line with archeological information [32] [84, 85]. A second list of 46 exclusively shared SNPs, found between two very distant isolates on a specific subbranch of AAnc5 (DRR034381 and DRR130203, pairwise distance: 417 SNP) AAnc5 is also shown in (Supplementary Table5).
2.2 in Silico spoligotyping and reconstruction of the CRISPR locus structure of AAnc 5 using TB-Annotator, IS6110 copy number and insertion sites
One of the strengths of the TB-Annotator pipeline is that is also reconstructs the CRISPR by integration of the CRISPR-builder TB results [58, 59]. The 20 AAnc5 isolates showed 6 different spoligotype patterns, which was quite unexpected for an L2 sublineage (Supplementary Table6); most of these patterns have been previously described in the SITVIT database (SIT1, SIT190, SIT269, SIT1364, SIT1674), however one remained undefined as SIT”X”. No SNP variants were found in spacers and repeats, but three isolates exhibited duplications: a duplication of sp65 for DRR034478 and SRR130160, and of sp50 for DRR034476 (Supplementary Table6). The phylogeny that can be derived from reconstruction of the CRISPR-Cas structure confirmed the SNPs results: it reveals sporadic deletions of cas genes, Rv2807c, Rv2808c and Rv2813c in some isolates (Supplementary Table7).
IS6110 remains a very important micro-evolutionary marker particularly in L2 that accumulated so many copies [48]. AAnc 5 strains were harboring between 14 and 22 IS6110 copies, and two specific copies were found in almost all these L2 isolates and not in other L2 sublineages: one copy was found around position 1724419 in Rv1527c (found in 15 of these isolates) and the second one was found, around position 2041756 (found in 19 of these isolates) (Supplementary Table5). DRR034455, DRR034471 and DRR034476 were showing the same CRISPR structure, however harbored different missing genes (see next paragraph). Using TB-Annotator, 14 of the AAnc5 isolates were predicted to be drug-sensitive and four were harboring mono-resistant mutations, two were MDRs (Supplementary Table 8).
2.3 Missing genes
Six isolates among the 20 Aanc5, apart from showing classical deletions (RD105, RD207, RD181 and PhiRv1), were harboring specific missing genes: as an example, DRR034363 had Rv1081 to 1084c deleted, DRR034416 was missing Rv1523 to Rv1526c. (Supplementary Table7). These deletions confirm that phylogenetically linked MTB genomes sometimes harbor strain-specific dependent deletions regions, due to recombination events and in relation to the high number of IS6110 copies.
2.4 In Silico VNTR copy number computation using CAST and comparison with other isolates from previous studies
Since TB-Annotator does not yet produce in Silico VNTR, we used the CAST pipeline to compute the VNTR on the 20 AAnc5 clinical isolates [55]. Apart from getting in Silico VNTR results, our first aim was also to compare VNTR with a very large set of isolates representative of Japan L2 biodiversity that had been published earlier [44]. Results were always partial, and no specific 15-9 VNTR signature could be obtained for the AAnc5 isolates. ETRC, QUB26 and QUB4156 could never be predicted. Depending on SRA quality, between 6 to 20 VNTR could be predicted (Supplementary Table 9). The VNTR results showed slight variation between isolates; eleven VNTR Loci were invariant in this collection (MIRU04, MIRU10, MIRU16, MIRU20, Mtub29, Mtub30, ETRB, MIRU24, MIRU27, Mtub34, MIRU39) whereas nine loci showed variation (MIRU02, MIRU40, Mtub21, QuB11b, ETRA, MIRU23, MIRU26, MIRU31, Mtub39). When comparing with an in-depth VNTR study performed earlier, AAnc5 was shown to belong to M10 and M37 respectively found in Russia and Singapore [51]. When comparing with a set of 5 reference japanese isolates (A05N056, ID381, 4558, 4994, 4991/M) that were described to represent the main L2 sublineages found in Japan, ID381 was sharing the same copy number with AAnc5 on 9 loci (Supplementary Table 9) [39, 44]. However, when comparing in Silico VNTR results with previous VNTR results from published studies in Korea, on the “K” strain that had been found in Korea we retrieved relatively poor similarity [38] (Supplementary Table 9).
2.5 Comparison between TB-Annotator and CAST Server results on prospective Drug Susceptibility Testing and on spoligotyping results
When comparing the results obtained using TB-Annotator and CAST, the predictive DST results were identical (Supplementary Table8 and Supplementary Table9). Identical results were also obtained on the classical 43 spacers-format spoligotype reconstruction with a minor and yet unexplained discrepancy on a single spacer of a single isolate, DRR034366, for which CAST predicted SIT250 whereas TB-Annotator predicted SIT290 (Supplementary Table9).
2.6 AAnc 5 is part of the G3 endemic L2 Ancestral strains cluster in Japan
We compared the 42 SNPs table found in the Asian ancestral 5 sublineage (Supplementary Table4) with the ones found in sequences of 5 japanese reference isolates, compared to the K1-K2 epidemic strain [39]. According to the definition made by these authors 5 L2 subgroups (G1/2, G3, G4, G5/6 and M) could be defined in Japan based on 10 phylogenetically selected SNPs. From the ID381 strain, a member of the “G3” genetic group, described for the first time in Kobe and Osaka in 2006 [39], and by comparing the SNPs specific of each sublineage (Supplementary file of [39]), we came to the conclusion that the G3 group was new and did not fit with the Shitikov et al. classification scheme. AAnc 1, AAnc2 and AAnc3 were found to respectively match with G1/2, G4 and G5/6 in [39], however no equivalent was found for “G3”.
By comparing SNPs list, we also found that the Tochigi province AAnc5 cluster of strains was sharing 36 out of 42 SNPs with the G3:ID381 reference strain found in Kobe and Osaka (Supplementary Table4). 6 SNPs only (C587945T, G765140C, G1202113A, AGGGAG1476812A, G3148446C and G3820365A) were not found in the ID381 strain. We concluded that the Tochigi strains were highly likely to be historically related to the Kobe and Osaka 2006 G3 group and to the sequenced reference isolate ID381. Accordingly, we propose to retain the common SNPs described by Wada et al. and this study as characteristic of this sublineage, that we suggest to name Asia Ancestral 5 (AAnc5) to fit with Shitikov et al. nomenclature [31]. By digging more in-depth into the comparative SNPs list between our study and the former Kobe-Osaka study, we found that DRR034489 was the closest isolate to the G3 ID381 reference sharing 40 more SNP exclusive to the G3 group, whereas another cluster of 3 genomes were more distant but were sharing 15 more SNPs with ID381 (results not shown). As mentioned above, two very distant genomes, DRR034381 and DRR130203 (417 SNPs pairwise distance) where also sharing 46 supplementary SNPs that were not found elsewhere (Supplementary Table 5).
Discussion
We describe in this paper an historical sublineage of L2 based on samples collected in central Japan, Tochigi Prefecture, former Shimotsuke province, that we coined as AAnc5. The identified sublineage is highly related to the Japanese G3 group defined in 2012 [39]. This strengthens the relevance of this sublineage and shows that it was transmitted historically in several Japanese cities. The chronology of emergence of this sublineage relatively to other L2 sublineages was positioned in Shitikov’s L2 diversification scheme. Its position relative to previously described lineages clearly showed that it should be qualified as ancestral according to the current use of this terminology and appeared just after AAnc1.
Retracing the early history of L2 in China, and South-East and in East Asia, in relation to China’s TB history, might be quite difficult [8]. Early TB outbreak history could be related to migrations of people from the 5th century BC to the third century AD [84, 85]. Tuberculosis was known to be present in premodern time in Japan under the name of rôga, that was used in chinese medicine [86]. There are many traces in Chinese medical history texts of a disease that can be identified as tuberculosis [87]. The China-south Japan relation is known since very early in human history (1 century AD) through the Kan-no-Wa-no-Na-no-kokuōin, the seal of the King of Na, close to today’s Fukuoka; The King of Na was a vassal of the Han Emperor [88]. The hundred of small states that were present in Japan during these early years left the place to a progressive domination of the Yamato province during the Kofun period. In parallel, a south to north progression of rulers was seen in relation to the adoption of an organizational system based on the one of the Chinese empire [54]. However simultaneously, the Yamato kings were also in contact with the former Korean Baekje kingdom [89]. The Shimotsuke province, where the isolates we are describing are coming from, was in permanent contact with the Yamato court since the Kofun period [90]. This organization contributed to the creation of permanent capital cities in Japan instead of moving capitals, and a progressive development of the plains north of Edo (Tokyo), with a direct link to the former Shimotsuke known today as Tochigi province [54].
Our results on the historical connection between two locally endemic clones could only be obtained thanks to the availability of public Short Read Archives and to the development of a new bioinformatical pipeline that analyses more than 50,000 characters that builds highly resolutive phylogenetic trees, and includes repeated sequences information (PE-PPE genes, CRISPR and Insertions sequences). With the increasing power of computational analysis of genomes and using supercomputing centers, it now becomes possible to analyze not only Single-Nucleotide Polymorphisms, but also repetitive sequences, which paradoxically were at the start of molecular epidemiology [91-93], however whose analysis was almost totally supplanted by WGS [56, 57, 94]. Going more in-depth into the complex historical phylodynamics history of all MTC lineages (made up of various waves of expansion and extinction history) will be facilitated by the use of all available markers and the increase into study size, towards 100,000 genomes analysis, to better disentangle all the threads between ancient and recent historical event that shaped today’s pandemia, and to obtain an improved understanding of the historical pandemic in relation to ancient/modern population migrations [8, 95, 96].
In this study, we also wanted to compare the in Silico VNTR results of the AAnc 5 with a representative set of similar results from Kobe and Osaka [44]; however, since in Silico VNTR had produced sub-optimal results and since initial comparison between in Silico and in Vitro VNTR results provided limited results and is conceptually difficult using Illumina Short-reads sequencing, we did not deepen these results here since too many missing data would have jeopardize the results. The few results we obtained do not provide any information on a possible route of introduction of AAnc5 into Japan since M10 is found in Russia, whereas M37 is found in Singapore. These limited results also shows the current technical limits of in Silico VNTR results production; nanopore sequencing could solve this issue in a near future [97]. However when looking at SNPs dataset and the correlation between SNPs and VNTR results, it became evident that G3 and AAnc5 were sharing common ancestors.
According to its branching point, the AAnc5 lineage has emerged after AAnc1 but before AAnc 2 and 3. The very existence of a more important number of L2 ancestral clones had been guessed by the discovery of the AAnc4 [32]. The AAnc5 isolates harbored a set of unique characteristics and were shown to be distant one from the other by between 200 to 400 SNPs thus definitively excluding a contemporary outbreak history of these isolates.
Among this set of unique characteristics was the presence of a non-synonymous mutation in rpoC, of variations in spoligotypes as well as in VNTR structures, reinforcing progressively the evidence of historical links within the G3 and AAnc5 clusters. It is well known that rpoC mutations are compensatory mutations mainly found in epidemiologically-successful isolates, most of the time modern Beijing strains, that contain specific rpoB mutations [10, 98,Gygli, 2021 #12772]. The rpoC compensatory mutation could be an adaptative trait, that was suggested to explain the epidemiological success of MDR-TB L2 isolates in outbreaks such as the Central Asian or the B0/W148 outbreak in Russia [6, 99] or the epidemiological success in Georgian prisoners [36]. Merker et al. have shown recently in the central Asian clade that rpoC mutations may arise within the entire gene and authors have shown that these SNPs are found in epistasia with rpoB mutations, hence these mutations are designated as adaptative or compensatory [6, 10, 98]. However this is not the case here, since AAnc5 were all, except for the most recent last six isolates, drug-susceptible and hence, this rpoC mutation should not be considered here as a compensatory mutation, but rather as a phylogenetical marker of AAnc5 [36]. It is noteworthy that this SNP was not present in the G3 group from Kobe and Osaka. Hence, given this slight rpoC difference between G3 and AAnc5, it could be interesting in the future to try to search if there is a significant statistical difference between drug-resistance emergence in one or the other cluster group. We do not know if such mutations could also provide another selective advantage, independently of rifampin drug-resistance, for example on fitness, however the recent paper in prisoners in Georgia showed that compensatory mutations and patient incarceration were two independent factors associated with increased transmission, thus creating a “perfect storm” for MDR-TB transmission [36]. The precise physiological or fitness consequences of potential rpoC mutations in drug-sensitive isolates should be the focus of other in-depth studies; the rpoC SNP found in AAnc5 is non-synonymous but we did not investigate in this study if it produces a change in protein structure and has functional consequences [10]. Variations in genetic backgrounds inside L2 (intra-lineage diversity) could be linked to the early evolutionary history of AAnc5, in relation or not, to yet unknown mutator effects specific to this L2 sublineage [9].
Dating of strains diversification is important to identify the conditions that fostered past epidemics. Here, according to consensus molecular clock of TB on the middle term, the time frame of coalescence of AAnc5 sublineage landmarks could be between 280-310 years ago before present for the closest isolates, and 760-800 years ago for the farthest ones [73]. However we should be cautious in our dating estimation. Indeed, we based this estimate on L4 mutations rates [96]; we know that calibration of the molecular clock varies a lot according to samples, time-frames and lineages (between 0.04 to 2.2 SNPs per-genome-per-year) [100]. Dating can be comforted if adequate correlation exists between historical, genomic, epidemiological and demographic facts. In a former similar study, three scenarios have been tested for the dating of a L4.2 sublineage occurring in Japan and Turkey using historical, anthropological, human genetics, paleopathological events and genomics [73]. Here, for the time-being, we have few representatives of AAnc 5 within all Japan, with no independent data such as genetic backgrounds of the host, or paleontological data. We may only try to correlate our dating with historical epidemiology.
TB Transmission could have been low in early Japan during the Kofun period or later, until Japan reached a certain demography. During the 19th century, the transmissible nature of TB was not admitted, it was rather believed that the disease was hereditary [101]. Within the post-Meiji society however, interest in workers’ health emerged, more for strategical than for humanitarian reasons [101]. Apart from the lack of bacteriological knowledge, fights between western and Chinese medicine lasted during centuries in Japan and bacteriology was only slowly introduced [101]. Statistics on death causes had started to be collected by the Imperial state in 1880, however tuberculosis deaths were counted in the respiratory diseases category, and it is only in 1883 that “pulmonary consumption” was specifically recorded as a cause of death [101, 102]. A few medical doctors at the japanese central hygiene office in 1903 were pinpointing the fact that infectious diseases first start within workers and then spread towards the general population [86]. For tuberculosis, it was linked to women working in silk and cotton spinning factories, the major driver of Japanese economy at that time [86, 103]. Hence, governments started to take TB seriously when bad health of women working in textile industry was statistically significant [86]. Mortality rate were as high as 2300 per 100,000 inhabitants in these working women, to be compared with 700 per 1000,000 in young women of the same age in the general population [101]. For a province such as Tochigi, (former Shimotsuke), with a population varying between 400,000 and 500,000 between 1720-1850 and such an incidence, around 4600-5750 cases/year could be expected in women. Then, three TB epidemic peaks are described in Japan: a first one from 1880 to 1900, a second from 1900 to 1919 with an historical peak in 1918, a third one corresponding to stabilization at a high level between 1919 till 1950 and lastly a fast decline since 1950 [101]. Women working in textile industry were paying the biggest death toll as mentioned, but men working in mining were also hit, however death causes were not distinguished from silicosis [86]. Indeed, after silk industry developed in Japan (60% of the raw silk in the world was made in Japan at the beginning of the 20th century [104], mining industry became the cornerstone of the imperial state expansion, and hence, health of male workers also became of strategic national interest for the growing military empire [86]. In the Tochigi province, the copper Ashio mine, whose exploitation started beginning of the XVIIth century was responsible of one-fourth of Japanese copper production in 1884; it closed down only in 1973. One of our estimation of the expansion date of the AAnc 5 is compatible wih the opening of the mine in Ashio [105]. Ashio mine could thus have been the very location of AAnc5 expansion and diversification. Looking at the historical spreading dynamics of the AAnc5 group in Japan could be done in the future by looking at finer geographical gradients of all L2 ancestral sublineages found in Japan.
Another origin of Japanese AAnc5 could be the importation of L2-ancestral clones by forced workers from Korea or China who were made to work in the mining industry from 1939 till 1945 [86]. Indeed 300,000 Koreans and 38,935 Chinese workers, mostly men were forced to work for Japan during these years [86]. We could not investigate in this study the potential distribution of G3 isolates within Japan, however our formal proof that the G3 and AAnc5 are linked and deeply rooted in Japan and that, are likely to have created independent local outbreaks in the past, is a first insight into this complex history. Future study could try to dissecate further, using combination on recent statistics on the total of epidemic cases, in combination to genomic characterization and geographical distribution, the global and the local history of ancestral L2 lineages in Japan [106].
The existence of a relatively high diversity specific to Tochigi and likely Kobe and Osaka, reminds of historical transmission restricted to a circumscribed area. Of course, such a picture of an endemic past outbreak is more easily observed in islands settings as was shown recently in New Zealand in Maori people who are hosting a specific “CS” (Colonial S-type) L4.4 Sublineage [96]. Islands are indeed excellent settings to distinguish endemic from imported species and tuberculosis history, as all infectious diseases, is definitively linked to human migration and local and global demographic history [107]. Since Japan is now a country with a very low tuberculosis prevalence rate (13 per 100,000 in 2017), it is an excellent setting to try to identify the historical events linked to past tuberculosis outbreaks [73].
Conclusions
Thanks to TB-Annotator, a new bioinformatical pipeline that analyzes large amounts of the information contained in Sequence Read Archives, we mapped a local ancestral L2 sublineage onto the global MTC phylogeny, and we designated it as Asia Ancestral 5 (AAnc5), close to the formerly described G3 group in Kobe and Osaka. This sublineage now appears together with another one in Shitikov’s evolutionary scheme. This sublineage possesses many specific characters allowing to distinguish it from all other historical clones described so far for the “Beijing” (Lineage 2). This finding opens new ways of research, to look for the history of spreading of this cluster within all Japan, and further in relation to Korea, to South China or elsewhere in South-East Asia and East Asia. This pipeline will permit new investigations on the spatial and dynamics history of tuberculosis to be made.
Data Availability
All data are publicly available or contained in the article as supplementary files.
Funding
This study did not receive specific funding. The authors declare that they have no conflict of interest.
Legend to Supplementary Material
Supplementary Table1: List of all 190 SRAs of the L2 lineage investigated in this study, with a majority of L2 from central Tochigi, and other L2 representing the diversity of L2 sublineages described so far, with the exception of Pacific RD150, main SNPs assessed, and assignment to Shitikov classification scheme or newly assigned for Asia Ancestral 4 and 5.
Supplementary Table2: List of SNPs catalogs and/or other markers, used to classify SRAs in TB-Annotator
Supplementary Table3: Full Excel Matrix results produced on the L2 Lineage for the 190 L2 isolates and list of 680 SRAs used in the study (Sheet 1: General)
Supplementary Table4: list of 42 SNPs specific of Asia Ancestral 5
Supplementary Table5: SNPs Pairwise distance between the 20 Asia Ancestral 5 isolates, extra list of 46 SNPs found between two AAnc5 isolates.
Supplementary Table6: IS6110 insertion sites on the 20 Asia Ancestral 5 isolates
Supplementary Table7: summary of missing genes for the 20 Asia Ancestral 5 isolates
Supplementary Table8: SNPs detected in drug resistance genes in the 20 Asia Ancestral 5 isolates
Supplementary Table9: in Silico VNTR results on the 20 Asia Ancestral 5 and comparison of predictive DST and spoligotyping using CAST or TB-Annotator
References
- 1.↵
- 2.↵
- 3.
- 4.↵
- 5.
- 6.↵
- 7.↵
- 8.↵
- 9.↵
- 10.↵
- 11.↵
- 12.↵
- 13.
- 14.
- 15.↵
- 16.↵
- 17.↵
- 18.↵
- 19.↵
- 20.↵
- 21.↵
- 22.↵
- 23.↵
- 24.↵
- 25.↵
- 26.↵
- 27.↵
- 28.↵
- 29.↵
- 30.↵
- 31.↵
- 32.↵
- 33.↵
- 34.↵
- 35.↵
- 36.↵
- 37.↵
- 38.↵
- 39.↵
- 40.
- 41.
- 42.↵
- 43.↵
- 44.↵
- 45.↵
- 46.↵
- 47.↵
- 48.↵
- 49.↵
- 50.↵
- 51.↵
- 52.↵
- 53.↵
- 54.↵
- 55.↵
- 56.↵
- 57.↵
- 58.↵
- 59.↵
- 60.
- 61.
- 62.
- 63.
- 64.
- 65.
- 66.
- 67.
- 68.
- 69.
- 70.
- 71.
- 72.
- 73.↵
- 74.
- 75.
- 76.
- 77.
- 78.
- 79.
- 80.
- 81.↵
- 82.↵
- 83.↵
- 84.↵
- 85.↵
- 86.↵
- 87.↵
- 88.↵
- 89.↵
- 90.↵
- 91.↵
- 92.
- 93.↵
- 94.↵
- 95.↵
- 96.↵
- 97.↵
- 98.↵
- 99.↵
- 100.↵
- 101.↵
- 102.↵
- 103.↵
- 104.↵
- 105.↵
- 106.↵
- 107.↵