Tuberculosis (TB) is a chronic infectious disease with complex epidemiological characteristics (Abel et al., 2014; Barberis et al., 2017; Adhikari, 2022). It is an ancient disease causing illness and death in humans since the primordial times (Palomino et al., 2007). Over 10.9 and 1.4 million TB cases and deaths were reported in 2021 alone (WHO, 2022). Mycobacterium tuberculosis complex (MTBC), the etiological agents of TB, includes Mycobacterium tuberculosis (M. tuberculosis, Mtb) (Sakula, 1982), M. africanum (de Jong et al., 2010), M. bovis, M. microti, M. pinnipedii, M. orygis, M. mungi, and M. suricattae (Gagneux, 2018). The human-adapted members of MTBC (Mtb and M. africanum) are further classified into nine lineages with distinct geographic structures (Coscolla et al., 2021).

Ethiopia is a country in East Africa with high TB incidence and a complex Mtb population genetic structure (Mekonnen et al., 2019c). This complex population genetic structure is likely shaped by the consecutive entry of strains through long trade across the red sea, human migration over the millennia (Comas et al., 2013, 2015), and the long-established co-evolutionary trajectory (May and Anderson, 1983). The incidence of tuberculous lymphadenitis (TBLN), which is clinically characterized by the inflammation of the lymph node (cervical, axillary, or inguinal) with or without TB constitutional symptoms (cough, weight loss, fever, and night sweats), is exceptionally high in Ethiopia (Berg et al., 2015). However, causal factors leading to a high incidence rate of TBLN in Ethiopia remain speculative. Factors including ethnicity, Mtb strain type, HIV co-infection (Firdessa et al., 2013; Berg et al., 2015; Mekonnen et al., 2019b), over diagnosis (Iwnetu et al., 2009), and spill-over transmission of bovine TB (Ameni et al., 2011; Tadesse et al., 2017) were excluded. On the other hand, female sex, younger age (Mekonnen et al., 2019b), delayed diagnosis (Asres et al., 2019), and rural residency (Mekonnen et al., 2022) demonstrated a significant association.

The controversy about scientific evidence for “the incidence of TBLN with no lung involvement” is the subject of ongoing debate (Ganchua et al., 2020). For instance, Yew and Lee (1995) and Deveci et al. (2016) concluded that TBLN without pulmonary TB (PTB) is rare. To accept the theory that TBLN can exist as a distinct localized form of TB without pulmonary involvement, pathogenesis pathways other than lymphohematogenous spread, such as via tonsils and adenoids, must be proposed (Yew and Lee, 1995; Deveci et al., 2016). Herath and Lewis (2014) support the above claims and found a higher rate of concomitant PTB among extrapulmonary TB (EPTB) patients but with normal chest X-ray findings. In contrast to the above findings, a higher number of TBLN cases without concurrent PTB has been reported by multiple studies. Briefly, in Polesky et al.'s (2005) study, out of 106 TBLN cases, only 16 (15%) had evidence of pulmonary involvement. Of the total TBLN cases, concomitant PTB was reported among 20.8% of TBLN cases (Qian et al., 2018). In another study by Bomanji et al. (2020), only 100 (28%) out of 358 EPTB patients showed pulmonary involvement. Nevertheless, to disprove either of the two scientific narratives, a clinical, radiological, and pathological investigation should be coupled with bacteriological methods.

Whole-genome sequencing (WGS) technology demonstrated an incredible potential for investigating the population's genomic structure and transmission clusters (Gardy et al., 2011; Walker et al., 2013; and Meehan et al., 2018). Core genome multilocus sequence typing (cgMLST) and single-nucleotide polymorphism (SNP) with 5–12 SNP or allelic differences were suggested for identifying epidemiological links with temporal scales (Gardy et al., 2011; Walker et al., 2013; and Meehan et al., 2018). Recent transmission index (RTI) is the measure of the severity of the Mtb epidemic and indicates the rapid progression of the infection into clinical illness (Tanaka and Francis, 2005; Tessema et al., 2013).

A few number of studies examined whether the high incidence rate of TBLN in Ethiopia was associated with Mtb lineages and sub-lineages (Firdessa et al., 2013; Ejo et al., 2021). These findings indicated that lineages and sub-lineages were identical between PTB and TBLN forms. Furthermore, the rate of clustering was also similar among the two clinical forms (10% for PTB and 11% for TBLN) (Firdessa et al., 2013). Additionally, a comparative literature review in Africa showed no difference with regard to lineages and sub-lineages between the two forms of TB (Mekonnen et al., 2019a). However, as can be seen from the figures in our previous research studies (Mekonnen et al., 2019b, 2022), TBLN showed a geographic pattern similar to the MTBC geographic structuring. To explain further, TBLN is mainly reported in countries having a high prevalence of geographically restricted specialist MTBC (sub)-lineages (Mekonnen et al., 2019a). Hence, we speculated the presence of MTBC genomics influences on TB clinical phenotypes. Hence, this study was designed to investigate the genomic diversity and transmission dynamics of Mycobacterium tuberculosis (Mtb) isolates obtained from TBLN and PTB cases in Northwest Ethiopia using WGS data.

Whole-genome sequencing was performed for 161 Mtb isolates using Illumina NovaSeq 6000 technology. The overall sequence read quality was good. The sequence read length varied from 35 to 151 nucleotides (mean 148). The mean and median sequencing coverage were 267 and 201 times, respectively. The mean (Std.error) and median (interquartile range/IQR) of SNP events were 1,163.62 (32.795) and 1,027.0 (498), respectively. The maximum deletion event among isolates was 610, and the minimum number of deletion was 114. The highest SNP, deletion, insertion, and substitution events were observed in L3 and L7 isolates than in L4 sub-lineages (Supplementary Table 1). An independent sample t-test was carried out to compare the mean evolutionary events between PTB and TBLN isolates. The analysis did not reveal any significant correlation (p > 0.05, Supplementary Table 2).

A few studies used the WGS technique to determine the genomic diversity and transmission dynamics of Mtb in Ethiopia (Comas et al., 2015; Yimer et al., 2016; and Welekidan et al., 2021). However, the aims of these studies were different from each other and also different from the present study. This study assessed and compared the genomic diversity and transmission dynamic of Mtb isolates obtained from PTB and TBLN cases in Northwest Ethiopia.

The WGS data unveiled a high genetic diversity of the Mtb population, L4 being the predominant (100/146, 68.5%), followed by L3 (29.45%). While the LAM and Haarlem isolates showed a low rate of SNP, insertion, substitution, and deletion, the L3 and L7 isolates demonstrated a high rate of mutations (Supplementary Table 1). The difference might be due to the nature (sympatric vs. allopatric) and timing of the pathogens–host association. When the specific host and Mtb coevolve together for a long, the host presents constant metabolic supplies. Thus, the metabolic genes of the Mtb become rendered useless and deleted (Ochman and Moran, 2001; Moran, 2002) as a way to cellular economization (Martínez-Cano et al., 2015).

The comparative Mtb genetic diversity showed similar lineages and sub-lineages, but it showed distinct genomic clusters among the two clinical forms (PTB vs. TBLN). The clusters became mixed of TBLN and PTB isolates at a higher threshold of allelic difference (>5 allelic difference). Faksri et al.'s (2018) study found genetic variants of Mtb that are commonly found in TB meningitis patients compared to PTB cases. To further support our finding, we re-evaluated the TBLN geographic distributions mapped in the previous studies at Africa (Mekonnen et al., 2019b), Ethiopia (Mekonnen et al., 2019b), and Amhara Regional State level (Mekonnen et al., 2022). Figures 6A, B shows the geographic pattern of TBLN similar to Mtb geographical structuring. Furthermore, we assessed the Mtb lineages and sub-lineages distribution among countries that had TBLN reports (Mekonnen et al., 2019a). As such, TBLN was mainly reported from West and East African countries, having a complex MTBC population genetic structure and high prevalence of Mtb Genotypes that have narrow host ranges (specialist genotypes). For instance, in West Africa, L5/6/Maf and L4.6.2/LAM-CAM are dominant specialist lineages and sub-lineages, respectively. Similarly in East Africa, L4.6/SIT37, L4.2.2/SIT149, L4.10/SIT53, L4.6.1/T2_Uganda, and several CAS sub-lineages are prevalent (Mekonnen et al., 2019a). Hence, the high incidence rate of TBLN in Ethiopia and other Eastern and Western African countries is likely the result of Mtb genomic variation. Negrete-Paz et al. (2021) also demonstrated intra-(sub)-lineage differences among TB clinical phenotypes. They briefly stated that strains of a sub-lineage vary in terms of immune response and virulence. Hence, these variations perhaps affect disease severity and clinical presentation (Negrete-Paz et al., 2021). Tuberculous lymphadenitis is the less severe form of TB (Hasan et al., 2009). It is known that long-established host–parasite coevolution leads to benign disease with insidious onset (May and Anderson, 1983). Further computational co-evolution model intuition could prove this scenario if it works for TBLN.

Unlike other zones of Amhara Regional State shown in Figure 6C, the low incidence rate of TBLN in the North Shewa Zone (Figure 6C light green color) could be due to environmental and host genetic differences besides Mtb genomics. Unlike other zones of Amhara, North Shewa is high land with cold weather.

Perhaps the other most compelling finding of this study is the higher CR and RTI in TBLN cases which are even slightly higher compared with PTB cases (p > 0.05). This is direct evidence showing that the high incidence rate of TBLN in the Amhara region is not just a result of the reactivation of latent TB. Rather, the finding indicated the rapid progression of clinical illness from primary infection or short latency. The high incidence rate of TBLN in younger age people, such as students (Table 2), was another evidence of the short latency nature of TBLN in Amhara Regional State. This rapid clinical progression of TBLN from primary infection was not due to immunodeficiency. This is because the immunohematological values of TBLN are relatively higher among TBLN compared with PTB (Mekonnen et al., 2023). Rather, it might be due to immunological tolerance resulting from sympatric-host pathogen association (Seal et al., 2021).

Due to the poor health system and socio-demographic factors, unrecognized extensive TB transmission might have taken place in Ethiopia for a long period. Such hidden extensive micro epidemics could probably push the co-evolution toward immune tolerance and then TBLN clinical phenotype. Host controls Mtb infections by activating innate and adaptive cellular immunities. Immune response against Mtb causes immunopathology in both the host and the pathogen. Hence, the host's ability to invest in self-toxic immune responses might have a limit. Beyond this threshold, the host switches from active resistance to tolerance strategy. Furthermore, tolerance to invading pathogens also has a threshold. Hence, the host's immune response strategy against Mtb is perhaps fine-tuned by the fitness effects of both the host and Mtb. This ability of the host for maintaining a balanced immunity and homeostasis during infection results in less severe forms of TB, such as TBLN, smear-negative PTB, lower cavitation, and negative chest X-ray finding (Seal et al., 2021). Taken together, the insights gained from this study may be of assistance to re-evaluate the current prevailing dogma which says- “no strain difference between PTB and TBLN isolates.”

L4.2.2.ETH forms two major branches containing a mix of clustered and unique strains (Figure 2). The L4.2.2.ETH topology in the cgMLST and SNP phylogeny is similar to Comas et al.'s (2015) phylogenetic tree. The phylogenetic position of L4.2.2.ETH is close to the Ural family, indicating their introduction into Ethiopia from Central Asian states (Mokrousov, 2012). Compared with L4.6, the arrival of L4.2 into Ethiopia is likely a recent phenomenon and possibly linked with the major north–south human migration during the reign of the Queen of Sheba (Comas et al., 2015). The LAM (L4.3.4.1 & L4.3.4.2) sub-lineages showed low deletion and substitution rates in this study. LAM is considered a generalist sub-lineage that is found in over 47 countries and highly prevalent in Europe with high genetic diversity. Its global success was likely driven by European migration and colonization (Stucki et al., 2016). The other highly clustered sub-lineage was Haarlem/L4.1.2.1 which is found in over 49 countries (Stucki et al., 2016). This sub-lineage is estimated to be introduced into Ethiopia at the time of early Portuguese contact with Ethiopia in the 16th century or later (Comas et al., 2015).

Lineage 3 is the second leading genotype found in this study, and this report is in line with previous studies (Tessema et al., 2013; Biadglegne et al., 2015; Mekonnen et al., 2019c; and Welekidan et al., 2021). There is a high genomic substitution, deletion, or SNP in L3 similar to L7. Our phylogenetic tree showed that the genetic diversity of L3 is very complex and likely shaped by several independent evolutionary events or multiple entries into Ethiopia (Figure 3). The L3 contains both specialist (isolates restricted in some countries including Ethiopia) and generalist genotypes (Comas et al., 2015; Tulu and Ameni, 2018). With a 90.2% probability, South Asia was predicted as the origin of all L3 strains with multiple (at least four) independent introductions into East and North Africa (Shuaib et al., 2022). L7 was restricted to Ethiopia and branched off around the time of initial human migrations out of Africa (Comas et al., 2015). Collectively, except L7, the other (sub) lineages followed the human migration; out of Africa (northern route via Egypt) (Pagani et al., 2015) and back migration into Africa in the different demographic histories (Pankhurst, 1964; Molinaro et al., 2019). Taken together, the population genomic structure of Mtb in Ethiopia is likely driven by long and stable host-Mtb reciprocal evolutionary changes (May and Anderson, 1983; Freschi et al., 2021) and consecutive entry of new variants via trade, missionaries, and human migration (Curtis, 2008; Comas et al., 2015). In addition to that, the poverty and high population density might have contributed to local adaptation and expansion of newly imported Mtb.

This WGS study gave high-resolution information about the role of Mtb genotypes in the clinical presentation of TB, genomic diversity, and transmission dynamics. However, it has some limitations. This study was carried out at the institutional level and included patients coming from different zones. Thus, clustered isolates did not mean they are epidemiologically linked, rather, it might be due to the widespread availability of the genotypes in the regional state. Similarly, the unique isolates might not be the result of reactivation, rather it might have been linked to cases in the community where the patient came. Hence, there might be a sampling bias. Furthermore, the sample size was small which precluded us from the in-depth assessment of the socio-demographic-related variables.

TBLN and PTB showed no difference at lineage and sub-lineage levels. However, the isolates from the two clinical forms formed a distinct cluster, indicating the presence of pathogen genomic factors. The higher rate of clustering and RTI among TBLN confirmed the rapid progression of TBLN from primary infection and/or short latency. Hence, the high incidence rate of TBLN in Amhara Regional State is likely shaped by Mtb genetic diversity and rapid progression from primary infection. Comparative community-based studies containing data from the environment and host genome could give a detailed view. Additionally, the characterization of lineage-independent genetic components might give further insight.

The original contributions presented in the study are publicly available. This data can be found here: NCBI, PRJNA975069.

The study was approved by research and Ethical Review Committee of Science College of Bahir Dar University, with reference number PGRCSV/111/2012. Informed written consent was obtained from each participant before data collection. This study was conducted in accordance with the Declaration of Helsinki. All the information obtained from the study subjects were coded to maintain confidentially. The patients/participants provided their written informed consent to participate in this study.

DM: substantial contributions to the conception of the work, data collection, laboratory works, analyses, interpretation of data for the work, and drafting of the work. AM and EN: substantially contributed to the conception, methodology, validation, monitoring of the work, secured partial fund for the work, and revising it critically for important intellectual content. BA: substantially contributed to the bioinformatics analysis and revised it critically for important intellectual content. AAA and AB: substantially contributed to supervision, monitoring of the laboratory work, and revising it critically for important intellectual content. EA and SH-L: participated in the DNA extraction, data analysis and interpretation, and revising the draft critically for important intellectual content. CJ: prepared the sequencing library, substantially contributed to the data analyses, interpretation of data for the work, and revised the draft critically for important intellectual content. AA: substantially contributed to the conception, methodology, and revising the draft critically for important intellectual content. LH-L: monitoring and coordinating the DNA extraction and WGS, substantially contributed to the bioinformatics, and revising the draft critically for important intellectual content. All authors critically revised and approved the version to be published.