LysG-driven transcriptional network rewiring underlies lineage-specific phenotypes in Mycobacterium tuberculosis

Domain: biology
Published: 2026-03-23
Editorial outcome: accept

Abstract
The Mycobacterium tuberculosis complex (MTBC) includes ten human-adapted lineages with varying geography and pathogenicity. Lineage 1 (L1) shows low virulence while Lineage 2 (L2) is hyper-virulent, more transmissible, and associated with drug-resistance. We performed comparative analyses integrating whole-genome sequencing with transcriptomic and proteomic profiling of L1 and L2 clinical strains under two in vitro growth conditions. Transcript-protein correlations varied by strain and gene category, suggesting lineage-specific post-translational regulation. Expression differences scaled with phylogenetic distance, one in three SNPs affected gene expression. A new transcriptional regulatory model identified master transcription factors, linked to the sigma factor network, whose targets were differentially expressed between L1 and L2. For instance, DosR regulon proteins had higher basal levels and exhibited a stronger nitric oxide response in L2. Time-course experiments involving LysG (Rv1985c) induction and wild-type H37Rv under hypoxia and subsequent reaeration confirmed that LysG contributes to reduced metabolic activity, thereby promoting increased tolerance to the novel tuberculosis drug bedaquiline in L2 strains relative to L1. Overall, our findings show how limited genetic variation in the MTBC can yield major phenotypic differences through differential regulation of key transcriptional networks.

Introduction
Adaptive evolution is mainly driven by mutations. While the phenotypic consequences of mutations can be easily mapped in the context of strong selection such as antibiotic pressure, the phenotypes emerging from more subtle evolutionary changes are more elusive.

Mycobacterium tuberculosis (Mtb) is the etiological agent of tuberculosis (TB), the major cause of human mortality due to a single infectious agent1. Mtb is an obligate human pathogen and comprises ten phylogenetic lineages, each with a distinct geographic distribution. Together with several animal-adapted lineages2,3, these make up the Mycobacterium tuberculosis complex (MTBC). Even though the MTBC harbors little genetic diversity compared to other bacteria4, MTBC clinical strains differ significantly in virulence and transmissibility5,6. However, linking this phenotypic variation to the limited genetic diversity in the MTBC has been challenging7. Among the ten human-adapted lineages of the MTBC, some have received particular attention because of their clinical significance. Of note, Lineage 2 (L2) strains are geographically widely distributed and have been associated with high virulence in macrophages and various animal models and rapid disease progression8, increased transmission9, and a propensity for antibiotic resistance10. Conversely, Lineage 1 strains (L1) are geographically localized mainly around the rim of the Indian Ocean, and have been linked to asymptomatic TB11, reduced virulence in infection models12, a reduced transmission potential9 and less likely to be antibiotic resistant5,13. Importantly, in geographical areas where both lineages are present, such as in Vietnam, L2 seems to be supplanting L1 as the predominant form of the MTBC9,14. These features point to important physiological differences between L1 and L2. Although, many studies have described phenotypic differences in MTBC clinical strains5,15,16, the molecular mechanisms underlying these differences remain poorly understood. To date, the majority of insights into the molecular processes have been gained from work on a few laboratory-adapted strains of the MTBC17–20 that do not capture the entire phylogenetic diversity of clinical strains21.

Microbiologists are increasingly applying systems biology approaches and multi-layer omics technologies to interrogate bacterial systems and their interactions with the human host22–25. These approaches profile various types of biomolecules, typically genome, transcriptome, proteome and, more recently, protein complexes that collectively describe the molecular makeup of specific cells and cellular states and predict interdependencies of different levels of gene expression along the axis of the central dogma26–31. In the MTBC, genomic profiling and measurements of the transcriptional response to a particular stress have been frequently used20,32,33. Yet, computational approaches that can extract mechanistic insights from the statistical associations in a multi-layer omics dataset need further development. Moreover, a growing body of evidence shows that direct functional information can be retrieved from the state of the proteome, making quantitative proteomic data an asset and a relatively recent addition to integrated multi-layered omics analyses of bacteria20,34. The mass spectrometric technique DIA-MS, combines data independent acquisition (DIA) with massively parallel targeted analysis of the acquired data35 and thus offers the high degree of reproducibility, consistency, and throughput needed for large cohort studies. It has been optimized in model experiments with Mtb36–38.

Here, we describe a multi-layer omics dataset that includes transcriptome and proteome data generated from a set of genome-sequenced clinical MTBC strains. Specifically, three L1 and L2 strains, respectively, were studied during two in vitro growth conditions. We use these data to systematically examine the strength of post-transcriptional and translational regulation in the MTBC across various functional categories. Moreover, we tailor a genome scale transcriptional framework – GenSTrans – to identify transcription factors (TFs) that regulate their target genes differently between L1 and L2 strains. We further show that ~40% of the observed significant changes can be explained by a small set of TFs and the number of SNPs between each pair of strains determines the extent of the significant changes at both mRNA and protein level. Our model highlights four master TFs (DosR (Rv3133c), LysG (Rv1985c), Lsr2 (Rv3597c), and Rv0691c) that are known to interact with several sigma factors, including SigB. We discuss how the observed differential expression of DosR proteins, most likely triggered by SigB in response to nitric oxide (NO) stress, correlates with shorter growth arrest in L2 strains compared to L1. To validate the biological significance of LysG in differentiating L1 and L2 clinical strains, we design and implement a hypoxia–reaeration time-course experiment in which LysG is induced, thereby supporting the predictive power of our transcriptional model. In summary, we demonstrate the potential of network modeling to describe molecular mechanisms that differentiate MTBC L1 strains from L2. We propose that these differences may explain some of the phenotypic changes observed in clinical settings where these strains circulate.

Results
Multi-layer omics profiling of MTBC clinical strains
To establish the basis for the integrated, multi-layer analysis of MTBC L1 and L2 clinical strains, we selected three drug-susceptible strains from each lineage39,40 to represent the greatest phylogenetic diversity within each lineage based on whole genome sequence data (Fig. 1A). We then grew the six strains in vitro to mid-log phase and profiled their transcriptomes and proteomes in technical triplicates using RNA sequencing and DIA-MS, respectively (Fig. 1B, Fig. S1). We quantified 3,933 transcripts across the samples (transcript per million, TPM > 1). DIA-MS quantified 21,184 peptides consistently across samples from which 2,479 proteins, corresponding to 62% of Mtb’s open reading frames, were inferred (Fig. 2A, Fig. S2A). Through the targeted data analysis pipeline OpenSWATH41 and using a Mtb proteome library37 as prior information, three or more proteotypic peptides per protein were quantified for ~76% of the detected proteins (Fig. 2B). Next, we examined the quality and reproducibility of each of the omics layer individually. The median Coefficient of Variation (CVs) across the technical triplicates was ~12% at both the transcript and protein level, indicating comparable data quality between the two platforms (Fig. S2B, C). Overall, we were able to quantify 94% of detected proteins across all samples by direct measurement. Where this was not possible, we inferred upper bounds of protein expression from local background values (Fig. 2C). The correlation between biomolecular features of the replicates was high for both the transcriptomic and proteomic data (Fig. S2D, E). We assessed the technical and biological variability of the respective datasets. The results indicated an ascending molecular variation with increasing genomic distance within and between the two MTBC lineages (Fig. 2D). Our analyses also resulted in a perfect unsupervised clustering of MTBC strains and lineages for both transcriptomic and proteomic data (Fig. S2F, G). Principal component analysis (PCA) suggested that the protein data were more informative with respect to separating strains by their phylogenetic relationship compared to transcriptomic data. Indeed, the first two principal components of protein data exhibited a good separation, resolving individual strains (Fig. 2E). By contrast, the resolution of PCA based on transcriptomic data was more limited and reflected the phylogenetic classification (Fig. 2F). In short, we have generated a dataset, which spans from the genome of fully sequenced clinical isolates of MTBC to transcriptome and proteome. The multifaceted nature and deep coverage of the high data quality provided a solid base on which to build our downstream analyses.

Post-translational regulation is prevalent in MTBC clinical strains
We used our combined multi-layer omics dataset to systematically explore the extent of correlation between mRNA and protein levels in MTBC clinical strains and to assess the relative impact across gene functional categories annotated by Mycobrowser42 between the two lineages. We computed the Spearman rank correlation between averaged intensities of mRNA and protein products of the same gene across technical triplicates in the six clinical strains studied. We found that the magnitude of the correlation, and hence the assumed extent of post-transcriptional and translational regulation, varied between functional categories. The overall correlation was 0.46, slightly higher than in previous reports20. The strongest and weakest signals of post-transcriptional and translational regulation were observed in the “virulence, detoxification, and adaptation” (Spearman’s rho: 0.33) and “lipid metabolism” (Spearman’s rho: 0.69) group, respectively. Remarkably, the L2 strains presented a significantly higher mRNA to protein correlation compared to L1 for both regulatory (T-Test’s p-value = 0.015) and virulence genes (T-Test’s p-value = 0.0036) (Fig. 3A). A growing body of evidence has revealed that the extent of post-transcriptional regulation in the MTBC is largely limited, as indicated by the high correlation (r = 0.93) observed between mRNAs and ribosome-protected fragments (RPFs)19. This suggests that the modest genome-wide correlation between mRNA and proteins may instead originate from post-translational regulation in the MTBC. Comparing the global correlation observed in the MTBC to that in other species revealed that the MTBC’s potential for post-translational regulation could significantly decouple its proteome from its transcriptome, to a degree comparable to that observed in humans, despite the generally lower complexity of bacterial systems23,43. An earlier study showed that mRNAs and proteins of 1,921 Escherichia coli genes were largely correlated (r = 0.8), supporting the notion that post-transcriptional and translational regulation are relatively limited in bacteria44. Other studies in yeast and metazoan species showed mRNA-protein correlations at the gene-to-gene basis of above 0.545–48. A possible explanation for the apparent increased post-translational regulation in the MTBC may be the presence of the Pup-proteasome system responsible for protein degradation and homeostasis that is rarely found in other bacteria including E. coli49. Due to the substantial post-translational regulation in the MTBC, RNA data in isolation can lead to incomplete or misleading mechanistic understanding. This underscores the need for high quality quantitative proteome data20. For instance, our data (Supplementary Data 1) revealed that T cell antigens (both MHC class one and two) are significantly regulated at the RNA level in L2 strains compared to L1 (MCHI’s p-value = 4.57E-5 and MHCII’s p-value = 0.0071). However, this pattern of the differential regulation could not be observed at the protein level (MHC I: p-value = 0.07 and MHC II: p-value = 0.55) (Fig. S3).

Transcript and protein abundances differ substantially between MTBC L1 and L2 strains
Next, we performed differential protein and transcript expression analyses of the set of L1 and L2 strains. We identified 578 differentially abundant transcripts and 390 differentially abundant proteins (fold change > 1.5 and adjusted p-value < 0.01) (Supplementary Data 1). The vast majority of the changes occurred in non-essential genes involved in survival, infection and persistence of the bacilli, whereas the expression of essential genes50 remained largely unchanged (p-values = 6.9E-7 and 2.17E-6 on mRNA and protein level, respectively). Gene set enrichment analyses of the differentially abundant transcripts and proteins highlighted various regulons and functional categories that distinguish the two lineages. These included the Dormancy survival Regulator (DosR) and Iron-dependent Repressor (IdeR) regulons. DosR regulon proteins contribute to maintenance of cell viability during transition to a non-replicating state induced by various stresses such as macrophage engulfment, hypoxia and nitric oxide (NO)51. The abundance of these genes was significantly differently regulated in L2 strains compared to L1, both on the transcript and protein level (Fig. 3B, S4A). For transcripts this observation has been reported previously26,52,53 and it is confirmed here using both transcriptome and proteome data from the same samples. Chao JD et. al. showed that the abundance of DosR regulon proteins is increased upon knockout of the serine/threonine-protein kinase PknH54. Since we found the abundance of PknH decreased in L2 strains, we decided to examine whether there was a link between the PknH downregulation and DosR upregulation. We hypothesized that PknH could be involved in the differential regulation of DosR regulon proteins, in addition to the causative variant C601T55, through either DosR or an alternative transcription factor (TF). To test this hypothesis, we profiled the proteome of a PknH knockout and overexpressing strain generated from the laboratory strain H37Rv, along with the wild type strain in technical triplicates, using the same method as described above. However, we could not replicate the results shown by Yossef Av-Gay and colleagues54 and DosR regulon proteins remained unchanged upon PknH knockout and overexpression (Supplementary Data 1).

Our analyses also pointed to the regulon IdeR, (repressed in iron-replete conditions) as being significantly differentially expressed between the two lineages (Fig. 3C, S4B). Iron acquisition systems in pathogenic bacteria are critical virulence factors and the MTBC requires iron to be able to grow in culture and during infection56,57. Hence, the increased expression of IdeR regulon proteins in L2 strains could boost the potential of the iron acquisition system and therefore strengthen L2’s success during infection58. Lipid metabolism related proteins that showed the highest mRNA-protein correlation also appeared significantly differentially abundant between the two lineages (Fig. 3D, S4C). This functional category contains 272 genes in Mtb whereas, in comparison, the E. coli genome encodes only ~50 such proteins. The amount of energy the MTBC spends to transcribe and translate these genes indicates their critical roles in the life cycle of the MTBC. Proline-glutamic acid (PE)/proline-proline-glutamic acid (PPE) genes, many of which are considered as virulence factors, also displayed differential regulation between the lineages at both RNA and protein level (Fig. 3E, S4D). These proteins are challenging to quantify since they are, to a large extent, either secreted or localized on the cell surface. We next examined TFs of which 138 out of 214 (64%) known factors could be detected at the protein level. We found that at least 24 TFs showed changes in their abundance at the protein level between the two lineages (Fig. 3F, S4E). Taken together, our multi-omics data revealed important differences in the functional and regulatory organization of MTBC clinical isolates from L1 and L2.

The number of single nucleotide polymorphisms (SNPs) is linked to the magnitude of transcript and protein expression differences between strains.
The large extent of the observed transcriptional and translational changes between the two lineages led us to examine the role of SNPs, the most common form of genetic variation in the MTBC. To this end, we first calculated the genomic distance based on number of SNPs between each pair of the tested strains. The average pairwise genomic distances within L1 and L2, and between the two lineages, were 873, 446 and 1,761 SNPs, respectively. These results are consistent with results reported previously5 that showed that the genetic diversity within L1 was about twice that of L2. We then examined the extent to which the observed changes at the transcript and protein level correlated with the number of SNPs in the corresponding genomes. Our analysis revealed a Spearman correlation between the genomic distance, measured in number of SNPs, and the differentially regulated mRNAs and proteins of 0.73 and 0.83, respectively (Fig. 4A, B). The lower proteomic coverage diminished the slope of the corresponding regression line but it almost met the slope calculated for transcripts following a coverage-based normalization. The elucidation of the one-to-one relationship and exact mechanism of each causative SNP goes beyond the scope of this analysis due to the complexity of potential mechanisms. For instance, genetic variation might affect regulatory proteins such as TFs by either changing the affinity of a given protein to DNA (both synonymous and non-synonymous SNPs), or modifying a regulatory protein (non-synonymous SNPs), and therefore influence gene expression. Hence, we included both synonymous and non-synonymous SNPs in our analysis. Overall, our data indicate that, on average, every three SNPs (1/slope of regression line) caused a significant expression change that can be detected at either the transcript or protein level or both.

A genome-scale transcriptional model identifies divergent regulators between L1 and L2 strains.
The differential expression of a high number of TFs between the two lineages observed in our dataset led us to ask to what extent changes in transcript or protein abundance can be explained by direct TF regulation. To address this question, we leveraged an existing transcriptional framework, the Mtb’s Environment and Gene Regulatory Influence Network (EGRIN) model, to identify TFs causally linked to the observed changes in transcript and protein abundance59. The EGRIN model defines clusters of co-regulated transcripts (modules) and their corresponding TFs. It requires dual enrichment analysis to connect a given regulated gene set to the modules and the modules to TFs (Fig. S5A)60. Applying the model to our transcriptomic dataset highlighted 15 modules that were i) associated with at least one TF and ii) significantly enriched in transcripts that were differentially abundant between the two lineages (adjusted p-value < 0.05) (Fig. S6A). The TF DosR was identified with high significance through four different modules as a putatively causative TF for divergent gene expression in L1 and L2 strains (Fig. S6B). In addition, we also found that KstR2 (Rv3557c) could regulate its target genes differently between the two lineages. This TF was identified through two different modules and controls the expression of a small set of genes involved in the utilization of cholesterol, the main carbon source used by the MTBC during infection61. Whereas EGRIN provides a high degree of sensitivity and is well suited to pathway-level analysis provided that the model includes pathways of interest, the model could not satisfactorily address our quest for the identification of lineage-specific regulators for two reasons. First, the EGRIN model contains only 66 TFs and is therefore not a genome-scale model. Second, each module encompasses only 13 genes on average, which mitigates the specificity of the model. Consequently, the Mtb EGRIN based analysis can suggest several TFs regulating a given module, making it impractical to identify the causative TF with confidence (Fig. S6A). To improve on these limitations, we devised a new genome-scale transcriptional framework.

We repeated the analyses with a new genome-scale transcriptional network model, GenSTrans, that was developed based on ChIP-seq experiments62. The network consists of 143 TFs and their corresponding sets of target genes represented as specific sub-networks. Each sub-network on average consists of 64 genes. The high density of information in the network is expected to boost the specificity of the analysis compared to the EGRIN model. We overlaid the transcript data of this study on the GenSTrans network and examined the degree of enrichment of differentially abundant transcripts in specific sub-networks using hypergeometric test followed by Benjamini-Hochberg based multiple testing correction (Fig. S5B). This analysis suggested 28 TFs that differentially regulated gene expression between L1 and L2 strains (Fig. 5A). Of the 28 TFs identified by GenSTrans, 17 were not identified by the EGRIN model due to insufficient genome coverage (Fig. 5A). Our new model could recapitulate at least one TF corresponding to 11 modules also identified by the EGRIN model (Fig. S6A). Of note, GeneSTrans more precisely pointed to one or two TFs per EGRIN’s module in cases where the EGRIN model has suggested several equally likely candidates. For instance, the EGRIN model identified module 446 that could potentially be regulated by either Rv0081 or Rv1460, or a combination of both (Fig. S6A). Rv0081 and Rv1460 bind to 647 and 11 genes, respectively, according to ChIP-seq experiments62. Deconvoluting the signals is further complicated because the module overlaps with both sub-networks significantly. Nevertheless, the non-overlapping genes pointed to Rv1460 as the most likely causative TF explaining the observed pattern, as 9 out of 11 genes in the corresponding sub-network appeared regulated in the comparison of L1 and L2 strains.

Next, we returned to the original question about causation of differential gene expression between L1 and L2 strains. Specifically, we examined the extent to which the identified TFs could describe the observed expression differences. Our results revealed that at least 37% of the modulated transcripts could be statistically significantly linked to a relatively small set of TFs. A set of only four TFs – DosR, LysG, Lsr2, and Rv0691c – were sufficient to explain a quarter of the transcriptomic changes between L1 and L2 strains (Fig. 5B). DosR on both transcript and protein level and LysG63 on transcript level – not identified through our proteome profiling – were differentially regulated in their abundance very significantly between the two lineages. Notably, LysG was deleted in one of the L1 strains (N0153)64, which contributed significantly to the observed upregulation of LysG in L2 strains compared to L1. However, this does not compromise the results, as GenSTrans identifies TFs based on the differential expression patterns of their target genes between the two lineages. To further assess this effect, we performed the same network analysis across all possible L2 versus L1 strain comparisons (nine in total). LysG was not significantly identified in only two comparisons (N0155/N0072 and N0155/N0153), suggesting that its effects are consistent at the lineage level. Interestingly, and in contrast to DosR and LysG, Lsr2 and Rv0691c did not show significantly different transcript/protein abundance between the lineages. We speculated that Lsr2 and Rv0691c might undergo post-translational modifications (PTMs) regulation that affect their respective level of activity. In the MTBC, various types of PTMs have been identified on at least 2,617 proteins30,65. For instance, PknB, a serine/threonine-protein kinase, can phosphorylate Lsr266. The sub-networks of these four TFs shared several genes with each other, indicating their close functional relationship (Fig. 5C, D). Of note, DosR and LysG regulate 31 and 28 DosR regulon genes, respectively (Fig. 5E), according to ChIP-seq experiments, with a significant gap between them and other TFs in terms of their involvement in regulating DosR regulon genes62 (Fig. S7). These highly overlapping target genes belonging to the DosR regulon might indicate LysG as an alternative TF contributing to the DosR target regulation under certain conditions.

Hierarchical organization for interactions between TFs and sigma factors
In addition to 214 TFs, the tubercle bacillus has one housekeeping sigma factor, SigA, and 12 accessory sigma factors leading to the reprogramming of RNA polymerase (RNAP) and consequently the initiation of transcription of particular gene sets. The sigma subunit ensures specificity of the RNAP holoenzyme for special promoter sequences67,68. Gennaro and colleagues reconstructed a sigma factor regulatory network of Mtb. It consists of 13 sigma factors, seven anti-sigma factors, two anti-anti-sigma factors and 50 TFs69, altogether 72 genes. Of these, we were able to quantify the abundance of all 72 genes at the transcript level and 48 genes at the protein level, allowing us to ask how significantly the four master TFs could interact with the sigma factor network and generate differential responses to a given stress. We assessed whether the four master TFs identified in our analysis were present in the sigma factor network and how their interactions with sigma factors/anti-factors might affect gene expression (Fig. S8). Three of the four master TFs, Lsr2, DosR and Rv0691c, were among the 50 TFs represented in the sigma factor network. Therefore, our four TFs of interest were significantly enriched using hypergeometric test (p-value = 0.0027). According to the model, Lsr2 and DosR both interact with three sigma factors directly, whereas Rv0691c interacts with four sigma factors indirectly through their corresponding anti-sigma factors. To examine the potential of these three TFs of regulating the whole sigma factor network, we inferred the hierarchical organization of the sigma factor network using the hierarchy score maximization algorithm70. This algorithm shows how a given directed network is organized at multiple levels and elucidates the position of each component at the inferred levels. The results indicated that the sigma factor network has either three or four levels, and that DosR, Lsr2 and Rv0691c fall into the top-level of the network propagating regulations towards the lower levels (Supplementary Data 2). Moreover, SigB, controls the TF DosR expression69. Of note, only the protein RsfB (anti-anti-sigma factor) showed significant differential expression between the two lineages among all the sigma factors and their corresponding anti-sigma factors and anti-anti sigma factors quantified on transcript and protein level (Supplementary Data 1). Overall, these three TFs located at the top-level of the sigma factor network interact with eight sigma factors that could, in part, facilitate stress responses of L2 strains compared to L1. As the DosR regulon showed a higher basal expression in L2 strains, we next considered the extent to which DosR regulon proteins would respond differently respond to a relevant stress.

Differential proteomic and phenotypic response of L1 and L2 strains to nitric oxide stress
To determine whether the molecular differences observed between L1 and L2 strains, such as in the DosR regulon, are functionally relevant, we exposed the strains to NO stress for 24 hours and analyzed the samples with respect to their growth status and protein expression levels. NO stress simulates the physiological environment that the bacilli experience following engulfment by macrophages71. Specifically, we (i) quantified DosR target proteins and (ii) verified the growth arrest duration after NO stress (see Materials and Methods). The expression of DosR regulon proteins increased more strongly in response to NO in L2 strains compared to L1 (Wilcoxon signed rank test’s p-value = 0.00022) (Fig. 6A). Given the tight interaction with DosR and sigma factors, we also examined differential expression of sigma factors in L1 and L2 on treatment with NO. Our differential protein expression analysis showed that the strength of SigB induction in response to NO differed between L1 (fold change = 1.24, adjusted p-value = 0.006) and L2 (fold change = 1.94, adjusted p-value = 3.37E-6) strains, while other quantified sigma factors remained unchanged in both lineages (Supplementary Data 3). Notably, RsdA and UsfX also exhibited differential NO responses between L2 and L1 strains; however, they did not interact with the four major TFs identified through our transcriptional model (Fig. S8, Supplementary Data 3). This observation highlights a potential molecular mechanism underlying the stronger induction of DosR proteins in L2 strains following NO stress compared to L1 strains. Altogether, DosR regulon proteins showed not only a higher basal expression but also a stronger induction following NO challenge in L2 strains when compared to L1.

Based on these molecular data, we hypothesized that following NO exposure, L2 strains, would be better capable to restart their growth compared to L1 strains. We tested this hypothesis by measuring the duration of growth arrest in two strains from each lineage in technical triplicates. The results demonstrated a significantly shorter growth arrest in L2 strains compared to L1 (Wilcoxon signed rank test’s p-value = 0.0181), supporting our hypothesis (Fig. 6B). In summary, our findings show that TF network-based analysis can generate mechanistic hypotheses that are not directly apparent from the multi-omics data. We then validated the network-driven hypothesis experimentally and found that differential regulation of bio-molecular networks might drive distinct phenotypes in MTBC clinical isolates.

LysG induction accounts for one-third of the expression differences between L1 and L2 strains
Among the four master TFs identified by our transcriptional model, LysG has been largely understudied in the scientific literature, and we could not leverage publicly available data to highlight its biological significance. However, we consider it a master regulator differentiating the two lineages for three reasons: first, it was identified by GenSTrans; second, its affected subnetwork, comprising 58 genes, was the largest (Fig. 5A); and third, it is heavily involved in regulating DosR regulon genes (Fig. 5E). Hence, we conducted a transcriptional regulator-induced phenotype (TRIP)72 screen in the H37Rv background, in which ~200 Mtb TFs were individually induced, to examine the extent to which each of Mtb TFs reveals phenotypes under hypoxia and normoxia. The results indicated that only a handful (~5 TFs) of TFs, including LysG, showed a significant phenotype (altered growth rate) when the two conditions were compared. These screening results go beyond the scope of this study. Nevertheless, we designed and implemented a time-course experiment in which the LysG induction and control (empty plasmid) strains were induced one day prior to the experiment, followed by 7 days of hypoxia and 5 days of reaeration (Fig. 7A). We measured CFU on days 0, 7, 8, 9, 10, 11, and 12, and performed RNA seq analysis on days 0, 2, 4, 7, 8, 9, 10, 11, and 12 (Fig. 7A, S9). The CFU results revealed a significant difference in the growth rate of the strains during the reaeration phase. Of note, LysG is upregulated under hypoxia, and therefore this difference in growth rate was not expected during the hypoxic phase (Fig. 7B). Next, we conducted differential expression analysis between the induced and control strains at each time point to examine how many genes were up- or downregulated as a result of LysG overexpression under different conditions. The results showed that approximately 500 genes were upregulated and 500 downregulated (fold change > 1.5 and adjusted p-value < 0.01) during hypoxia, and around 1,000 genes were upregulated and 1,000 downregulated during the reaeration phase (Fig. 7C; Supplementary Data 4). Given that LysG displayed significantly higher basal expression in the control strain under hypoxia compared to reaeration, the more modest transcriptional response to LysG overexpression during hypoxia was expected.

We next asked how many of the genes regulated by LysG overlap with those differentially expressed between L2 and L1. For example, 51, 37, 24, 54, 36, 103, 111, 99, and 99 L2-upregulated genes (out of 275) were also upregulated on D0, D2, D4, D7, D8, D9, D10, D11, and D12, respectively (Fig. 7C). Similarly, 50, 60, 41, 83, 60, 91, 104, 109, and 101 L2-downregulated genes (out of 303) were also downregulated at the same time points (Fig. 7C & Fig. S10). We further performed hypergeometric testing to evaluate whether the L2/L1 up- and downregulated gene sets were significantly enriched within the LysG-regulated sets at each time point. The results showed that the downregulated L2/L1 genes were significantly enriched at all time points, while the upregulated genes were enriched in nearly all time points during the normoxic/reaeration phase, with the exception of D8 (Supplementary Data 4). Overall, these results validated that LysG is indeed a master regulator differentiating L2 from L1 strains, explaining approximately one-third of the significant transcriptomic changes (~200 out of 578 genes) - a substantially higher proportion than what was inferred from the transcriptional model alone. A possible explanation for this higher proportion could be that DosR and LysG appear to share some targets, making it difficult to disentangle the individual contribution of each TF, especially since LysG also affects DosR.

We further explored the function of the LysG up- and downregulated genes at each time point to identify which Mtb pathways were enriched. We retained only significantly enriched pathways (adjusted p-value < 0.05) that appeared in at least two time points. Major upregulated pathways enriched in the LysG-upregulated gene sets during hypoxia included “Information pathways” and “ABC transporters,” while “DosR regulon” and “Transcription factor” categories were highly enriched during the reaeration phase. This could explain why LysG overexpression significantly altered ~2,000 genes during reaeration. It appears that LysG is an upstream TF that upregulates ~80 other TFs (Fig. 7D). These findings further support the notation that LysG and DosR work closely together, in line with previously published ChIP-seq data62. (Fig. 5E). Interestingly, the “DosR Regulon” was also significantly enriched among LysG-downregulated genes during hypoxia. In addition, several metabolic pathways, such as the “TCA cycle” and “Microbial metabolism in diverse environments,” were significantly downregulated during reaeration (Fig. 7E). This observation is further supported by a recent study that demonstrated reduced intracellular ATP levels in L2 strains compared to L1 strains73. This reduced metabolic activity may contribute to increased drug tolerance, which could help explain why L2 strains are more associated with drug tolerance - such as to bedaquiline - compared to L1 strains (Fig. S11)73. Overall, this time-course experiment validated LysG as a master regulator differentiating the two lineages, demonstrated its close functional association with DosR, and provided a potential explanation for the increased drug tolerance observed in L2 strains.

Methods
Mtb strains and bacterial culture
We used six clinical isolates of Mtb (CT2018-00082, CT2018-00083, CT2018-00084, CT2018-00085, CT2018-00088, and CT2018-03241), three from L1 and three from L2, that are part of a recently described set of reference strains40. These strains were chosen in an attempt to capture the global genetic diversity of Mtb, both within and between lineages. We grew bacteria in a modified 7H9 medium supplemented with 0.5% w/v pyruvate, 0.05% v/v tyloxapol, 0.2% w/v glucose, 0.5% bovine serum albumin (Fraction V, Roche) and 14.5 mM NaCl. With respect to conventional 7H9 culture medium, we omitted glycerol, tween 80, oleic acid and catalase. For global expression profiling experiments we cultured the bacteria in 1l bottles containing large glass beads to avoid clumping and 100 ml of media. We incubated the cultures at 37°C and rotated them continuously on a roller. For NO stress experiments we grew the bacteria in 50ml conical screwcap tubes containing either 17ml (proteomic profiling) or 10ml of culture medium. The cultures were incubated at 37°C on an orbital shaker. All strains were inoculated and cultivated under identical conditions, and cultures were harvested at a matched physiological state (mid-log phase; OD₆₀₀ = 0.5 ± 0.1) to ensure comparability across RNA-seq and proteomics datasets (Fig. S1). We have cultivated each MTBC strain under a given condition in technical triplicates, unless stated otherwise. Technical replicates refer to independent cultures and measurements of the same strain under identical conditions, whereas biological replicates correspond to different MTBC strains.

Transcriptomic profiling of L1 and L2 strains
We transferred a 40 ml aliquot of bacterial culture in mid-log phase (OD600 = 0.5 ± 0.1) into a 50ml Falcon conical tube containing 10 ml ice. We harvested the cells by centrifugation (3,000×g, 7 min, 4°C), re-suspended the pellet in 1 ml of RNApro solution (MP Biomedicals) and transferred the suspension to a Lysing matrix B tube (MP Biomedicals). We disrupted the bacterial cells using a FastPrep24 homogeniser (40s, intensity setting 6.0, MP Biomedicals). We clarified the lysate by centrifugation (12,000×g, 5 min, 4°C), transferred the supernatant to a clean tube and added chloroform. We separated the phases by centrifugation (12,000×g, 5 min, 4°C) and precipitated the nucleic acids from the aqueous phase by adding ethanol and incubating at -20C overnight. We performed a second acid phenol extraction to enrich for RNA. We treated our samples with DNAse I Turbo (Ambion), and removed stable RNAs by using the RiboZero Gram Positive ribosomal RNA depletion kit (Epicentre). We prepared the sequencing libraries using the TruSeq stranded Total RNA kit (Illumina) and sequenced them on a HiSeq2500 run in high output mode (50 cycles, single end).

The resulting reads were mapped to the Mtb H37Rv reference (NC_000962.3) genome using BWA79 (ver. 0.7.13); the resulting mapping files were processed with samtools80 (ver. 1.3.1). Per-feature read counts were performed using the Python module htseq-count (ver. 0.6.1p1) and Python (ver. 2.7.11). We performed differential expression analysis using the R package edgeR81 and R (ver. 3.4.0) to identify Lineage-specific transcriptional changes.

Environment and gene regulatory influence network (EGRIN) based transcriptional model
We used the EGRIN model to analyze our transcriptional data in the exact same way that described before 60. The model describes a set of modules that each one contains a dozen of co-regulated genes. The ChIP-seq and TFOE data paved the way to correspond each module to its potential regulators namely TFs59. Here, we projected our significant differentially regulated mRNAs (fold change > 1.5 and adjusted p-value < 0.01) between L1 and L2 strains onto the modules to figure out how significantly each module was enriched. The enrichment analysis was performed using hypergeometric test followed by Benjamini-Hochberg (BH) multiple testing correction (adjusted p-value < 0.05). To visualize a given identified TFs, we included all its corresponding modules even if they were not significantly enriched.

Genome-scale transcriptional model and networks integration
To reconstruct the Mtb transcriptional network, we retrieved a ChIP-seq data publicly available at the Mtb portal V262. We included both operon and direct interactions in the final network, which offered a great framework for genome-scale analysis. In contrast to the EGRIN model, the TFOE data was totally excluded as they largely feature indirect effects of a TF overexpression to the Mtb transcriptome. The reconstructed network comprised 143 TFs and 2,943 genes and can easily be extended upon the availability of ChIP-seq data for the remaining TFs. Next, a new set of regulons/modules, each contained the target genes of a given TF, were defined. Such a definition provided a one-to-one relationship between the modules and the TFs. Moreover, the average size of each module became ~65 genes, five times that of the EGRIN model, offering a larger specificity. The significantly regulated transcripts in L2 strains compared to L1 were tested against each module using a hypergeometric test followed by Benjamini-Hochberg (BH) multiple testing correction (adjusted p-value < 0.05). The identified TFs were compared to the EGRIN model results (see the color codes in the corresponding figures). Four TFs with the largest regulated gene size explaining 25% of the significant changes were introduced as the master regulators. This approach is easily portable to other bacterial systems upon the availability of respective ChIP-seq data. The TFs of the identified sub-network was integrated into the Mtb sigma factors network constructed by Gennaro and colleagues69.

Proteomic profiling using DIA-MS
We aliquoted 20 OD equivalents from mid-log phase (OD = 0.5 ± 0.1) bacterial cultures (e.g. 40ml of OD 0.5) into 50ml conical tubes. In the case of NO stress we harvested 8ml of bacterial culture. We harvested the bacteria by centrifugation (3,000×g, 7 min, 4°C) and washed the pellet twice with cold PBS to remove tyloxapol from the samples. The washed pellets were re-suspended in lysis buffer, which contained 0.1M ammonium bicarbonate, 8M urea and 0.1% RapiGest (#186001861, Waters). We then transferred the suspension to a Lysing matrix B tube (MP Biomedicals) and disrupted the bacterial cells using a FastPrep24 homogeniser (40s, intensity setting 6.0, MP Biomedicals). We clarified the lysate by centrifugation (12,000×g, 5 min, 4°C) and sterilized it by filtering it twice through a 0.22 μm syringe filter (Milipore) prior to further processing.

Next, we measured the protein concentration of each sample using BCA assay (A53225, Thermo Fisher Scientific). Considering the measured concentrations, we started with 60 ug protein and appropriate volume of lysis buffer was added up to 100 ul to equalize the final concentration of each sample. Then, 5 mM tris(2-carboxyethyl)phosphine (TCEP) was added to reduce protein disulfide bonds while the sample were incubated at 37oC for 30 min. Afterwards, the free cysteine residues were alkylated by adding 40 mM iodoacetamide and incubating for 30 min in the dark. The samples were diluted 6 times with 0.05 M ammonium bicarbonate to reach a final Urea concentration below 2 M. To digest proteins, 1.2 ug sequencing grade modified trypsin (V5113, Promega) was added (w/w 1:50). The samples were incubated overnight (~16 hours) at 37oC with gentle shaking of 300 rpm. To stop protein digestion, we acidified the samples (PH < 2) using formic acid followed by and incubation for 30 min with shaking of 500 rpm. We desalted the clear peptide solution using C18 MicroSpin columns (The Nest Group, 30-300 ug loading capacity). Next, we spiked in iRT peptides (Ki-P3002, Biognosys) which allowed us to normalize retention time during the data processing step.

We measured the samples on a TripleTOF 5600 mass spectrometer (AB Sciex) coupled to a nanoLC system (Eksigent) as described before38. The data acquisition was performed in DIA-MS mode using the 64 variable windows scheme.

The OpenSWATH workflow, containing three software tools, paved the way to process the data. We used the pan Mtb spectral library generated on the same type of mass spectrometry, TripleTOF 5600 37. OpenSWATH extracted chromatograms and assigned peak groups according to the prior knowledge, the pan Mtb spectral library41. Next, we performed PyProphet-cli to score the assigned peak groups using a semi-supervised algorithm. The data was filtered to global 1% FDR on both peptide and protein level. Eventually TRIC was used to integrate various information of each run and align the extracted and scored peak groups between runs offering a high degree of consistency between measurements82. To perform absolute (based on top 3 peptides and top 5 transitions) and relative quantification, we used the R package SWATH2stats83 and eventually aLFQ84 and MSstats85. Of note, MSstats applies a Linear Mixed Model (LMM) framework to perform differential expression analysis while explicitly accounting for variance both across technical replicates and between strains.

Assessment of the response to nitric oxide stress
To measure the proteomic changes following NO exposure, we grew 17ml bacterial cultures to mid-log phase (OD approximately 0.5) and harvested 8ml of the culture for proteomic characterization. We then added Diethylenetriamine/NO adduct at a final concentration of 1mM to the remaining culture and incubated it for a further 24h at 37oC. After this time, we harvested 8ml of the culture for proteomic characterization.

To assess the duration of the NO induced growth arrest, we grew duplicate bacterial cultures to mid-log phase (OD approximately 0.5) and added Diethylenetriamine/NO adduct at a final concentration of 1mM to the treatment tube with nothing added to the control. We followed the subsequent growth dynamics by measuring OD600.

Genetic, transcriptomic and proteomic distance analysis between clinical isolates
To compute strain distances on various levels, we assigned a particular cutoff to digitize them. On genome layer, we ignored various genetic insertions and deletions and only counted the numbers of distinct SNPs between a pair of Mtb strains. For mRNA and protein level, we considered a given cutoff (fold change > 1.5 and adjusted p-value < 0.01) and then the number of the significant changes were calculated. For each layer, 15 comparisons consisting of six intra lineage (three pairs within each lineage) and nine inter lineages comparisons were assumed. Afterwards Spearman correlation was computed between the genomic and either transcriptomic or proteomic distance. The slope of the respective regression line showed how many significant changes on protein and transcript level are caused by a given distinct SNP.

Sigma factor network hierarchical properties
To infer the hierarchical organization of the sigma factor network recently described69, we used the hierarchy score maximization algorithm70. We performed the algorithm for different number of levels k (2-6). In each analysis, it computes probability scores to assign a given node to layers of the hierarchical organization. Two criteria were assessed to determine the optimal choice of k. First, the enrichment of the downward direction compared to expectation that is summarized in the corrected hierarchy score. Second, an ambiguity score which quantifies the uncertainty of a node assigning to a given layer of the organization. These two measurements revealed that k=4 is the optimal.

Time-course experiment under hypoxia and reaeration
To investigate the LysG-dependent response of MTBC to hypoxia and subsequent re-aeration, a time-course experiment was performed using the laboratory strain H37Rv (wild type) and a strain overexpressing the TF LysG. Cultures were induced with anhydrotetracycline (ATc) at a final concentration of 100 ng/mL, 18 h prior to hypoxia induction. Early log-phase cultures were diluted to an optical density (OD₆₀₀) of 0.1, and 400 mL aliquots were transferred into hypoxia flasks. Day 0 (D0) samples were collected immediately before the onset of low oxygen gas flow and the induction of hypoxia.

Samples were collected across nine time points: D0, D2, D4, D7, D8 (+1), D9 (+2), D10 (+3), D11 (+4), and D12 (+5). D0 represented normoxic conditions, D2–D7 corresponded to hypoxic culture, and D8-D12 reflected re-aeration following reintroduction to aerobic conditions. RNA-seq data were obtained for all nine time points, while CFU measurements were collected at D0, D7, D8, D9, D10, D11, and D12. All samples were processed in technical triplicates, except for the H37Rv strain harboring the empty plasmid (control) at Days 10 and 11, which was conducted in two technical replicates.

Hypoxic conditions were established using a 1 L three-armed spinner flask (Corning), equipped with a magnetic stir paddle and two side arms sealed with O-ring caps to prevent unwanted gas exchange. One arm was used for sampling, while the second was fitted with two metal tubes capped with 0.2 µm filters for biocontainment and controlled gas exchange. Hypoxia was achieved by continuously flowing a low-oxygen gas mixture (0.2% O₂, balance N₂) at 0.15 ft³/min. Cultures were maintained at 37 °C with constant stirring (60 rpm) inside an incubator. After seven days under hypoxia, re-aeration was initiated by transferring the cultures to standard aerobic growth conditions as described above and incubating for an additional five days.

At each time point, 20 ml bacterial cultures were collected and pelleted by centrifugation (3,000×g, 5 min, 4 °C). Pellets were resuspended in Trizol reagent, transferred to lysing matrix B tubes (MP Bio), and mechanically disrupted using an Omni Bead Ruptor 24 Elite (3 cycles of 30 s at 6.5 m/s, with cooling on ice between runs). RNA extraction followed phase separation with chloroform using Heavy Phase Lock Gel tubes (Eppendorf), precipitation with isopropanol/high-salt solution (0.8 M sodium citrate, 1.2 M NaCl), and purification using the Qiagen RNeasy kit. Ribosomal RNA was depleted using the Ribo-Zero Bacterial rRNA Removal Kit (Illumina). RNA quality and integrity were confirmed using an Agilent Bioanalyzer.

cDNA libraries were prepared with the NEBNext Ultra RNA Library Prep Kit for Illumina, and library concentrations were quantified using KapaQuant (Kapa Biosystems). Single-end sequencing was performed on an Illumina GAIIx platform at the Northwest Genomics Center (University of Washington).

Raw reads were assessed for quality using FastQC, trimmed using Trim Galore, and aligned to the Mtb H37Rv reference genome (NC_000962.3) using STAR86 with the parameter --alignIntronMax 1 to prevent spliced alignment. Read counts were assigned using FeatureCounts, and transcript abundances were quantified as TPM values. Low-count genes were filtered using filterByExpr (edgeR)81, and normalization and differential expression analysis were performed using the voom-limma87,88 pipeline, which models the mean–variance relationship and computes observation-level weights for linear modeling.

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "This study provides a comprehensive multi-omics analysis of Mycobacterium tuberculosis Lineage L1 and Lineage L2 clinical strains, revealing how genetic differences translate to significant phenotypic variations. The research integrates genomic, transcriptomic, and proteomic data to elucidate lineage-specific regulatory mechanisms, particularly in stress response pathways. The findings highlight how limited genetic diversity can lead to substantial functional differences through differential regulation of key transcriptional networks.\nThis work is interesting and contributes to the development of new bioinformatics tools to understand how certain genetic variations are reflected in molecular differences that underlie phenotypic changes such as increased or decreased virulence of pathogenic bacterial strains. However, there are some areas where greater clarity or additional details could improve the overall understanding of the results.\nComments:\n1) The L1 strain 0072 exhibits greater variability compared to the others based on Figures S2B and S2C. The authors should consider discussing this variability in more detail and evaluate whether it could affect the accuracy of the results obtained for Lineage L1.\n2) In the legend of Figure 2, the strains represented by the different symbols should be specified. As discussed above, it would be interesting to see if the Lineage 1 strain represented by the light blue cross symbol, which deviates from the other two strains, corresponds to the strain named 0072. Additionally, line 898 has (D) repeated twice - the second occurrence should be corrected to (E) and the following (E) to (F).\n3) At line 168, it should be specified that the modest correlation refers to Lineage L1.\n4) It is unclear where the information on T cell antigens from lines 181 to 185 is derived from. Table S1? However, this table is cited for the first time in the following paragraph. It would be helpful to include a table with the gene names and respective values used for the graphs in Figure S3.\n5) From lines 190 to 193, it is stated that most changes occur in non-essential genes. Adding a column to the table specifying which genes are essential and non-essential would clarify this statement.\n6) At line 210, the work by Yossef Av-Gay mentioned should be cited.\n7) Clarify whether the numerical values in Table S1C refer to proteins or mRNAs.\n8) There is an empty space between parentheses at line 310.\n9) The fold change values reported at lines 360-361 for SigB do not match the values listed for SigB in Table S3. Furthermore, Table S3 highlights RsdA and UsfX as proteins that are significantly differentially expressed between Lineages L1 and L2, in addition to SigB. However, the text does not provide any discussion about the potential roles or significance of the differential expression of RsdA and UsfX between the two lineages. Additionally, at line 341-342 they say that ”only the protein RsfB (anti-anti sigma factor) showed significant differential expression between the two lineage among all the sigma factors …”, however RsfB shows no significant p-value in table S3. They should contextualize the differential expression of SigB, RsdA, and UsfX, and their potential connection to the lineage-specific differential regulation mediated by the DosR regulon, which is a major focus of the study. The current section on sigma factors appears somewhat disconnected and lacks integration of these related findings\n10) The section from line 382 to 387 is very confusing. The source of this information is unclear, and the experimental conditions of hypoxia/normoxia for the Transcriptional Regulator Induced Phenotype (TRIP) assay are not described in the materials and methods section. The authors refer to reference 69, where the TRIP screen was conducted using isoniazid as a stress condition. Therefore, it is unclear how the authors reached the conclusion that based on a TRIP assay under hypoxic and normoxic conditions, only 5 transcription factors, including Rv1985c, showed a significant phenotype in terms of altered growth rate. Without this crucial information, it is unclear where the authors obtained the evidence that led them to select Rv1985c for further time-course experiments under hypoxia and reaeration conditions. Additionally, there is no information about the plasmid vector where Rv1985c was cloned in front of an indubible promoter. This entire paragraph needs to be rewritten in a clearer and more transparent manner.\n11) It is not clear from which data the authors conclude that Rv1985c regulates approximately 80 other transcription factors upstream (lines 425-426).\n12) The part from lines 581 to 587 should be moved to the results or discussion section.\n13) In Figure 6, it should be specified what the numerical values on the y-axis correspond to.\n14) In the legends of Figure 7 and Figure S9, the meaning of OE and EP should be explained.\n15) Some bibliographic references are incomplete, missing page numbers.", "reviewer_2": "This study from Banaei-Esfahani et al reported correlations between transcript and protein across three L1 and three L2 M. tuberculosis complex representative clinical isolates, indicating lineage-specific post-translational regulation. Rv1985c has been identified as main master transcription regulator responsible for differential target gene expression in L2 compared to L1 strains.\nMinor comments\n- Methods, Mtb strains and bacterial culture: please provide rationale and/or references related to the use of modified 7H9 medium.\n- Methods, BWA, samtools, Python and R: please consider including appropriate references. I would also suggest considering including any specific settings used within each tool/script/module (e.g. default vs custom).\n- Methods, H37Rv refence genome: please specify the accession number and version considered.\n- Methods: the Authors mentioned the use of three biological replicates. Please consider to specify the number of technical replicates for each method for improved clarity.\n- Methods/Results, EGRIN: given the insufficiency of the EGRIN model, please consider either to remove or reduce the Results section on its output. Also, given the relatively poor overlap between EGRIN and GenSTrans, it does not seem to be meant as a further supporting evidence.\nMajor comments\n- Methods: because bacterial gene expression can vary significantly with growth rate, growth phase, and inoculum history, the comparability of RNA-seq data depends on ensuring that cultures are physiologically matched. I therefore recommend clarifying in “Mtb strains and bacterial culture” that the inoculum source and growth phase at sampling were standardized across strains.\n- Results: there are a number of sentences that would probably fit better under the discussion section. E.g. lines 169-181, but also those referred to previous studies and literature. Alternatively, cross-check whether the Journal allows for a joint results and discussion section and expand discussion points.\n- Results, lines 201-211: According to the Chao’s study mentioned, under normal growth (aerobic, early log phase), the basal expression of the DosR regulon was unaffected by the pknH deletion. Under NO induction, ΔpknH showed reduced induction of DosR regulon relative to WT, suggesting PknH is required for full induction under stress. It is not fully clear whether the Authors were not able to replicate the results during basal or NO induction conditions.\n- Results, “The number of SNPs is linked to the magnitude of transcript and protein expression differences between strains”: based on the data provided, it seems that only SNPs have been considered; however, I would assume that indels are also part of the genetic diversity of the strains under study. I would recommend to explain somewhere how indels have been considered and, since excluded from this analysis, this should also somehow covered under a “limitations of the study” paragraph.\n- Results, “The number of SNPs is linked to the magnitude of transcript and protein expression differences between strains”: while I understand that elucidating the precise one-to-one relationship and underlying mechanism of each causative SNP is beyond the scope of this analysis, given the complexity of possible regulatory interactions, I believe that a more detailed description of the relationship between SNPs and transcription factors (TFs) would be valuable. Since the authors focus on differentially expressed TFs, it would be informative to assess whether SNP distributions across TF genes differ significantly between lineages. Regardless of whether such differences are observed, this analysis could strengthen the link between genetic variation and transcriptional regulation. Similarly, this should somehow be applied to sigma and anti-sigma factors later used to support the role of the TFs selected in “Hierarchical organization for interactions between TFs and sigma factors”. Additionally, as other functional categories also include a high number of differentially expressed genes, the rationale for prioritizing TFs in the downstream analysis should be clarified.\n- Results: looking at Table S1, it appears that strain 0153 is primarily driving the large expression differences observed for Rv1985c between L1 and L2. Given that some L1 strains from the same repository have been reported to harbor a deletion of Rv1985c (doi.org/10.1038/s41598-018-22237-5), it would be important to verify the presence or absence of this gene across all isolates included in the analysis. Notably, strain 0153 also displays higher diversity on the PCA for protein data, as well as in the intra-lineage mRNA expression profiles. Interestingly, the magnitude of its differential expression (|log2FC| approx. 10) is comparable both within its own genotype and in inter-genotype comparisons, suggesting that the expression pattern is strain-specific rather than genotype-driven. This reinforces the need to confirm that gene Rv1985c is present and intact in all strains before interpreting genotype-level effects.\n- The possibility that gene Rv1985c is deleted in some L1 isolates should be addressed in the discussion. While the Authors conclude that it functions as a validated master regulator distinguishing lineages L2 and L1, the presence of a deletion in some L1 isolates challenges this interpretation and suggests that the observed differential expression may not uniformly reflect a lineage-level effect.\n- Please consider reporting BCCM/ITM IDs for all strains to avoid misunderstanding.\n- Discussion: please consider including a “limitations of the study” paragraph"}, "author_response": "Point-by-point response to reviewers’ comments\nManuscript title: Rv1985c-Driven Transcriptional Network Rewiring Underlies Lineage-Specific Phenotypes in Mycobacterium tuberculosis.\nWe would like to sincerely thank both reviewers for their time and thoughtful evaluation of our manuscript. We deeply appreciate the reviewers’ constructive and insightful comments, which have greatly contributed to improving the clarity, depth, and overall impact of our study.\nWe have carefully considered and addressed all comments point by point, as detailed below.\nIn the revised manuscript, newly added text is highlighted in green and underlined, while deleted text is indicated using strikethrough formatting. We believe these revisions have substantially strengthened the manuscript.\nReviewer #1 (Remarks to the Author):\nThis study provides a comprehensive multi-omics analysis of Mycobacterium tuberculosis Lineage L1 and Lineage L2 clinical strains, revealing how genetic differences translate to significant phenotypic variations. The research integrates genomic, transcriptomic, and proteomic data to elucidate lineage-specific regulatory mechanisms, particularly in stress response pathways. The findings highlight how limited genetic diversity can lead to substantial functional differences through differential regulation of key transcriptional networks. This work is interesting and contributes to the development of new bioinformatics tools to understand how certain genetic variations are reflected in molecular differences that underlie phenotypic changes such as increased or decreased virulence of pathogenic bacterial strains. However, there are some areas where greater clarity or additional details could improve the overall understanding of the results.\nWe would like to thank the reviewer for the positive and encouraging evaluation of our work. We are pleased that the study is considered comprehensive and relevant, and we greatly appreciate the thoughtful and constructive comments provided. These insights have been invaluable in improving the clarity, completeness, and overall quality of the manuscript. We have carefully addressed all points raised in detail below and implemented the corresponding revisions in the manuscript to enhance its presentation and scientific precision.\nComments:\n1) The L1 strain 0072 exhibits greater variability compared to the others based on Figures S2B and S2C. The authors should consider discussing this variability in more detail and evaluate whether it could affect the accuracy of the results obtained for Lineage L1.\nWe thank the reviewer for this comment. Below, we provide the exact coefficients of variation (CVs) across technical triplicates for each strain at both the protein and mRNA levels.\nProtein-level CVs (median):\nmRNA-level CVs (median):\nN0052: 0.1486; N0072: 0.1224; N0145: 0.1183; N0153: 0.1278; N0155: 0.1247; N0157: 0.1150; Overall CV: 0.1252\nAs shown, at the mRNA level (Fig. S2C), the median CV for strain N0072 is slightly below the overall median CV, indicating consistent transcriptional measurements across replicates. Since the transcriptomic data form the foundation of our downstream transcriptional network analyses, this observation supports the robustness of our conclusions. At the protein level (Fig. S2B), the median CV of N0072 is somewhat higher than the overall median. To ensure that this variability did not bias the results, we processed the proteomic data using MSstats, which explicitly models variance both between strains and across technical replicates using a Linear Mixed Model (LMM). Therefore, this higher variability was already accounted for in the statistical analysis. We thank the reviewer for highlighting this point and have added a clarifying sentence to the Methods section to explicitly indicate that variance across technical replicates was incorporated during MSstats-based processing (Lines 655-657).\n2) In the legend of Figure 2, the strains represented by the different symbols should be specified. As discussed above, it would be interesting to see if the Lineage 1 strain represented by the light blue cross symbol, which deviates from the other two strains, corresponds to the strain named 0072. Additionally, line 898 has (D) repeated twice - the second occurrence should be corrected to (E) and the following (E) to (F).\nWe appreciate the reviewer’s attention to this detail. The symbols corresponding to each strain have been indicated between Fig. 2E and Fig. 2F; this information was included in the original version. We have also corrected the panel labeling in the figure legend (the second (D) changed to (E), and the subsequent (E) changed to (F)) in the revised manuscript (Lines 1001-1003).\n3) At line 168, it should be specified that the modest correlation refers to Lineage L1.\nThe modest correlation mentioned in the text refers to the genome-wide mRNA-protein correlation observed across both L1 and L2 strains (modest in comparison with E.coli for instance), rather than to L1 alone. For reference, the strain-specific Spearman correlation coefficients between mRNA and protein abundances are as follows:\nN0052: 0.453; N0072: 0.454; N0145: 0.473; N0153: 0.448; N0155: 0.503; N0157: 0.450;\nOverall average: 0.463.\nWe have revised the corresponding sentence in the manuscript to clarify this point and avoid possible misinterpretation (Line 172).\n2\n4) It is unclear where the information on T cell antigens from lines 181 to 185 is derived from. Table S1? However, this table is cited for the first time in the following paragraph. It would be helpful to include a table with the gene names and respective values used for the graphs in Figure S3.\nWe thank the reviewer for this comment. The data underlying all volcano plots (Figs. 3, S3, and S4) are derived from Table S1A and B, which contains the results of the differential expression analyses at both the mRNA and protein levels. Specifically, Table S1 includes the gene names, fold changes, P values, and adjusted P values for all relevant pairwise comparisons across transcriptomic and proteomic datasets. For greater clarity, we have now cited Table S1 earlier in the text, at the point where the data were first used, ensuring that readers can directly trace the source of the information (Line 186).\n5) From lines 190 to 193, it is stated that most changes occur in non-essential genes. Adding a column to the table specifying which genes are essential and non-essential would clarify this statement.\nWe thank the reviewer for this observation. Throughout the manuscript, we have integrated multiple sources of gene annotations. Including all these annotations in Table S1 could potentially create confusion for readers. The classification of essential genes was based on the work by Christopher Sasse and colleagues (PMID: 21980284), who systematically identified essential genes in M. tuberculosis. We have now cited this reference in the main text to clarify the source of this information so that readers can readily access the complete list of essential genes if desired (Line 196).\n6) At line 210, the work by Yossef Av-Gay mentioned should be cited.\nWe thank the reviewer for this comment. The corresponding paper by Yossef Av-Gay and colleagues (PMID: 20630871) was already cited earlier in the manuscript (Line 207). However, we have now cited it again in the concluding sentence of the relevant section to ensure proper clarity (Line 214).\n7) Clarify whether the numerical values in Table S1C refer to proteins or mRNAs.\nThe numerical values in Table S1C refer to protein-level measurements. As clarified in the manuscript, we profiled the proteome of the PknH knockout, PknH overexpression (induced), and wild type M. tuberculosis H37Rv strains in technical triplicates. The relevant sentence in the text reads:\n“To test this hypothesis, we profiled the proteome of a PknH knockout and overexpressing strain generated from the laboratory strain H37Rv, along with the wild type strain in biological triplicates, using the same method as described above.” (Lines 211-213)\n8) There is an empty space between parentheses at line 310.\nWe thank the reviewer for this comment. The unnecessary space between the parentheses has been removed in the revised manuscript (Line 322).\n3\n9) The fold change values reported at lines 360-361 for SigB do not match the values listed for SigB in Table S3. Furthermore, Table S3 highlights RsdA and UsfX as proteins that are significantly differentially expressed between Lineages L1 and L2, in addition to SigB. However, the text does not provide any discussion about the potential roles or significance of the differential expression of RsdA and UsfX between the two lineages. Additionally, at line 341-342 they say that ”only the protein RsfB (anti-anti sigma factor) showed significant differential expression between the two lineage among all the sigma factors …”, however RsfB shows no significant p-value in table S3. They should contextualize the differential expression of SigB, RsdA, and UsfX, and their potential connection to the lineage-specific differential regulation mediated by the DosR regulon, which is a major focus of the study. The current section on sigma factors appears somewhat disconnected and lacks integration of these related findings.\nWe are grateful to the reviewer for several points that were raised here, and we address them one by one below.\nFirst, the apparent discrepancy in fold change values arises from the difference in how values were reported. All fold changes in the supplementary tables (e.g. Table S3) were expressed as log₂ fold changes, whereas those in the main text were presented as actual fold changes (i.e., 2 to the power of the log₂ values reported in Table S3).\nSecond, the sigma factor analysis in the manuscript primarily focused on those sigma factors that directly interacted with the four transcription factors - Rv3133c (DosR), Rv1985c, Rv3597c (Lsr2), and Rv0691c - identified as master regulators distinguishing L1 and L2 strains (Fig. S8). RsdA and UsfX were not discussed in detail because they did not directly interact with these four transcription factors. They appear in Table S3 highlighted because they also exhibited differential responses to nitric oxide (NO) between L1 and L2 strains, which was not the main focus of that section. We have now clarified this point in the revised manuscript (Lines 375-377).\nThird, it is important to note that Table S3B presents differential expression analyses comparing NO-perturbed strains (both L1 and L2) to their respective unperturbed controls, whereas Table S1 reports differential expression between L1 and L2 under control (non-perturbed) conditions. RsfB is significantly differentially expressed between L1 and L2 strains under control conditions (as reported in Table S1), which is why Table S1 is cited in that context (Line 357).\nFinally, among the sigma factors, SigB is the only candidate that shows a lineage-specific differential response to NO exposure (Table S3) while remaining unchanged between L1 and L2 under control conditions (Table S1). This observation suggests that SigB most likely plays a key regulatory role in mediating the lineage-specific activation of the DosR regulon in response to NO, as summarized in Fig. 6A.\n10) The section from line 382 to 387 is very confusing. The source of this information is unclear, and the experimental conditions of hypoxia/normoxia for the Transcriptional Regulator Induced Phenotype (TRIP) assay are not described in the materials and methods section. The authors refer to reference 69, where the TRIP screen was conducted using isoniazid as a stress condition. Therefore, it is unclear how the authors reached the conclusion that based on a TRIP assay under hypoxic and normoxic conditions, only 5 transcription factors, including Rv1985c, showed a significant phenotype in terms of altered growth rate. Without this crucial information, it is unclear where the authors obtained the evidence that led them to select Rv1985c for further time-course experiments under hypoxia and reaeration conditions. Additionally, there is no information about the plasmid vector where Rv1985c was cloned in front of an indubible promoter. This entire paragraph needs to be rewritten in a clearer and more transparent manner.\nWe thank the reviewer for this detailed and constructive comment. Rv1985c was selected as a candidate master regulator for three independent reasons, as outlined in the manuscript:\n“However, we consider it a master regulator differentiating the two lineages for three reasons: first, it was identified by GenSTrans; second, its affected subnetwork, comprising 58 genes, was the largest (Fig 5a); and third, it is heavily involved in regulating DosR regulon genes (Fig 5e).” (Lines 395-398)\nWe additionally referred to results from a Transcriptional Regulator Induced Phenotype (TRIP) screen performed to evaluate the impact of transcription factor induction on growth under hypoxic and normoxic conditions. This assay was used solely to demonstrate that Rv1985c, among a small subset of transcription factors (~5 TFs), significantly affected growth under these contrasting conditions, underscoring its biological relevance. The complete dataset and extended analyses from the TRIP screen fall beyond the scope of the present manuscript. Nevertheless, all data related to Rv1985c (Fig. 7A) are included here, encompassing both its effect on growth rate (Fig. 7B) and the time-course RNA-seq experiment (Fig. 7C-E), which validates the transcriptional network analysis results for this regulator.\nWe thank the reviewer for noting that the description of the time-course experiment under hypoxia and reaeration was missing from the Methods section. We have now added a dedicated subsection titled “Time-course experiment under hypoxia and reaeration” (Lines 687-728), which details the experimental setup, conditions, and the inducible expression construct used for Rv1985c. Of note, the selection of Rv1985c was based on the three reasons particularly the result of our transcriptional network analysis (first reason) as discussed clearly in the manuscript.\nThese revisions clarify both the rationale for selecting Rv1985c and the experimental procedures used to validate its role.\n11) It is not clear from which data the authors conclude that Rv1985c regulates approximately 80 other transcription factors upstream (lines 425-426).\nTo enhance clarity, we have now explicitly cited Fig. 7D immediately following the sentence referring to the regulation of ~80 transcription factors in the revised manuscript. (Line 442)\n12) The part from lines 581 to 587 should be moved to the results or discussion section.\nWe have already discussed the interaction between sigma factors and transcription factor networks in detail in the Results section (Lines 341-347). However, we agree that the lines mentioned should not appear in the Methods section. Accordingly, we have revised that part, ensuring that the methodological description remains concise (Lines 611-617).\n13) In Figure 6, it should be specified what the numerical values on the y-axis correspond to.\nWe thank the reviewer for this helpful comment. We have now revised Figure 6B to clarify that the numerical values on the y-axis represent the number of hours each strain experienced growth arrest following nitric oxide (NO) exposure. The corresponding figure legend has also been updated for clarity (Line 1044-1045).\n14) In the legends of Figure 7 and Figure S9, the meaning of OE and EP should be explained.\nWe thank the reviewer for drawing our attention to this issue. OE and EP refer to the OverExpressor and Empty Plasmid strains, respectively, corresponding to the inducible and uninduced (control) Rv1985c strains. We have revised the legends of Figure 7 and Figure S9 to explicitly define these abbreviations.\n15) Some bibliographic references are incomplete, missing page numbers.\nWe thank the reviewer for this comment. We have carefully reviewed all bibliographic references and ensured that all entries now include complete information, including page numbers and other missing details.\nReviewer #2 (Remarks to the Author):\nThis study from Banaei-Esfahani et al reported correlations between transcript and protein across three L1 and three L2 M. tuberculosis complex representative clinical isolates, indicating lineage-specific post-translational regulation. Rv1985c has been identified as main master transcription regulator responsible for differential target gene expression in L2 compared to L1 strains.\nMinor comments\n1) Methods, Mtb strains and bacterial culture: please provide rationale and/or references related to the use of modified 7H9 medium.\nThe modified 7H9 medium was supplemented with 0.5% (w/v) pyruvate, 0.05% (v/v) tyloxapol, 0.2% (w/v) glucose, 0.5% bovine serum albumin (Fraction V, Roche), and 14.5 mM NaCl. We employed this modified medium to ensure optimal and consistent growth of genetically diverse Mycobacterium tuberculosis complex (MTBC) strains, including M. tuberculosis, M. bovis, and M. africanum.\nThis formulation also provides cross-platform consistency across our multi-omics studies. For instance, tyloxapol was used instead of Tween 80, as the latter is less compatible with proteomics workflows. Using the same culture conditions across studies facilitates direct comparison and integration of multi-omics datasets. Pyruvate was included as an alternative carbon source to support the growth of strains inhibited by glycerol. Tyloxapol prevents cell clumping while avoiding metabolically active or variable surfactants that could interfere with downstream proteomic and metabolomic analyses. Glucose supplied an additional readily metabolizable carbon source, while bovine serum albumin and NaCl maintained nutrient stability and osmotic balance.\n2) Methods, BWA, samtools, Python and R: please consider including appropriate references. I would also suggest considering including any specific settings used within each tool/script/module (e.g. default vs custom).\nWe thank the reviewer for this helpful comment. We have now included the appropriate references for all tools mentioned in the Methods section.\n3) Methods, H37Rv refence genome: please specify the accession number and version considered.\nWe thank the reviewer for this comment. The reference genome used for M. tuberculosis H37Rv was NC_000962.3, and we have now included the corresponding accession number and version in the Methods section (Lines 577 and 723).\n4) Methods: the Authors mentioned the use of three biological replicates. Please consider to specify the number of technical replicates for each method for improved clarity.\n5) Methods/Results, EGRIN: given the insufficiency of the EGRIN model, please consider either to remove or reduce the Results section on its output. Also, given the relatively poor overlap between EGRIN and GenSTrans, it does not seem to be meant as a further supporting evidence.\nThe EGRIN model represents the former transcriptional regulatory network of MTBC (PMID: 27573104), and our intention was to provide a comparative overview between EGRIN and our GenSTrans model. As shown in Fig. 5A, several transcription factors identified by both models, including DosR (Rv3133c), are highlighted (green), demonstrating decent concordance between the two approaches. Of note, TFs highlighted in blue were not identifiable through the EGRIN model (Lines 292-295), which underscores one of the major ways our modeling strategy provides additional insight and value to the TB research community. We believe that including this comparison provides readers with a broader perspective on transcriptional network analyses available for MTBC.\nMajor comments\n6) Methods: because bacterial gene expression can vary significantly with growth rate, growth phase, and inoculum history, the comparability of RNA-seq data depends on ensuring that cultures are physiologically matched. I therefore recommend clarifying in “Mtb strains and bacterial culture” that the inoculum source and growth phase at sampling were standardized across strains.\nWe appreciate the reviewer’s attention to this methodological detail. Our experimental design explicitly accounted for potential variability due to growth rate, growth phase, and inoculum history. As mentioned in the Methods section under “Transcriptomic profiling”, all strains were cultivated under identical conditions and harvested at the same physiological state (mid-log phase; OD₆₀₀ = 0.5 ± 0.1). The inoculum preparation and culture handling were standardized across all strains to ensure comparability of RNA-seq and proteomics measurements. This is described in the main text and Methods section:\n“We then grew the six strains in vitro to mid-log phase and profiled their transcriptomes and proteomes in biological triplicates using RNA sequencing and DIA-MS, respectively (Fig. 1b, Fig. S1).”\nMethods part:\n“We transferred a 40 ml aliquot of bacterial culture in mid-log phase (OD₆₀₀ = 0.5 ± 0.1) into a 50 ml Falcon conical tube containing 10 ml ice.”\nFurthermore, the OD₆₀₀ measurements recorded at harvest are summarized in Fig. S1, confirming that all samples were collected at comparable physiological states. Nevertheless, we have now reiterated this point under ‘Mtb strains and bacterial culture’ for full clarity (Lines 556-559).\n7) Results: there are a number of sentences that would probably fit better under the discussion section. E.g. lines 169-181, but also those referred to previous studies and literature. Alternatively, cross-check whether the Journal allows for a joint results and discussion section and expand discussion points.\nWe thank the reviewer for this suggestion. We agree that several of the sentences mentioned contain discussion-like elements. However, this particular section presents a key result, namely, the overall mRNA-protein correlation (r = 0.46) across MTBC strains, which is best interpreted when placed in context with similar analyses in other bacteria (e.g., E. coli, r ≈ 0.8). Including this contextual information directly after the result helps convey that the relatively modest correlation observed in MTBC likely reflects pronounced post-translational regulation rather than post-transcriptional effects, consistent with the high correlation (r = 0.93) observed between mRNA and ribosome-protected fragment (RPF) levels in MTBC.\nMoving this information entirely to the Discussion section might mislead readers, as they would only realize at the very end of the manuscript that the modest correlation is not driven by post-transcriptional events. Therefore, we prefer to retain this discussion-like explanation immediately following the corresponding result to ensure clarity and avoid potential confusion, of course, provided this is acceptable within the journal’s formatting guidelines.\n8) Results, lines 201-211: According to the Chao’s study mentioned, under normal growth (aerobic, early log phase), the basal expression of the DosR regulon was unaffected by the pknH deletion. Under NO induction, ΔpknH showed reduced induction of DosR regulon relative to WT, suggesting PknH is required for full induction under stress. It is not fully clear whether the Authors were not able to replicate the results during basal or NO induction conditions.\nWe thank the reviewer for this comment. The description in our manuscript is derived from the study by Chao et al. (PMID: 20630871). As shown in Table 1 (both untreated and NO-treated ratios are summarized) of that paper, the DosR untreated ratio (ΔpknH / WT) was 1.69.\n“Further examination of the iTRAQ data revealed that of the 48 genes commonly regulated by DosR in response to NO, hypoxia, and the Wayne model of dormancy, 13 gene products were identified by iTRAQ, all of which had higher protein levels in the ΔpknH mutant after 48-h standing conditions (Table 1). With the addition of NaNO₂, these 13 DosR-regulated proteins were induced to similar levels in WT and ΔpknH (Table 1).”\nIn contrast, when we performed the same analysis under normal (non-NO) conditions, we did not observe differential expression of DosR regulon members between the ΔpknH and wild-type strains, and thus could not replicate the elevated basal levels reported by Chao et al.. For this reason, we extended our analyses to nitric oxide perturbed conditions, where lineage-specific differences in DosR regulon responses became evident (please also see the sigma factor part particularly SigB and the relevant network in Fig. S8).\n9) Results, “The number of SNPs is linked to the magnitude of transcript and protein expression differences between strains”: based on the data provided, it seems that only SNPs have been considered; however, I would assume that indels are also part of the genetic diversity of the strains under study. I would recommend to explain somewhere how indels have been considered and, since excluded from this analysis, this should also somehow covered under a “limitations of the study” paragraph.\nWe thank the reviewer for highlighting this point. Indeed, insertions and deletions (indels) represent an important component of genetic diversity within MTBC strains. However, indels were not included in the present analysis. The rationale for focusing exclusively on SNPs is that assessing the relative contributions of SNPs and indels to transcriptional and proteomic variation would require analysis across a substantially larger cohort of MTBC strains, likely hundreds, to achieve sufficient statistical power. Such an investigation clearly extends beyond the scope of the current study, as also mentioned in the following comment by the reviewer.\nIn this work, we focused on SNPs because they represent the most common and well-characterized form of genetic variation and allow a conceptual framework for understanding how genetic diversity influences transcript and protein abundance in a one-to-one manner, analogous to expression quantitative trait loci (eQTL) and protein quantitative trait loci (pQTL) analyses in other organisms. This analysis was primarily intended to provide the TB research community with a conceptual foundation for designing future large population-scale MTBC studies.\nIn addition to indels, differences in gene content could also contribute to functional diversity, and emerging long-read sequencing technologies are increasingly enabling such analyses. We have now added a statement in the Discussion section acknowledging that indels and gene content variations were not considered here and that their potential impact remains to be systematically explored in future large-scale MTBC population studies (Lines 475-476).\n10\n10) Results, “The number of SNPs is linked to the magnitude of transcript and protein expression differences between strains”: while I understand that elucidating the precise one-to-one relationship and underlying mechanism of each causative SNP is beyond the scope of this analysis, given the complexity of possible regulatory interactions, I believe that a more detailed description of the relationship between SNPs and transcription factors (TFs) would be valuable. Since the authors focus on differentially expressed TFs, it would be informative to assess whether SNP distributions across TF genes differ significantly between lineages. Regardless of whether such differences are observed, this analysis could strengthen the link between genetic variation and transcriptional regulation. Similarly, this should somehow be applied to sigma and anti-sigma factors later used to support the role of the TFs selected in “Hierarchical organization for interactions between TFs and sigma factors”. Additionally, as other functional categories also include a high number of differentially expressed genes, the rationale for prioritizing TFs in the downstream analysis should be clarified.\nWe thank the reviewer for this multifaceted comment.\nFirstly, single-nucleotide polymorphisms (SNPs) can influence gene expression either through cis-regulatory effects (when the SNP resides within or near the gene it regulates) or trans-regulatory effects (when it impacts other genes via intermediate factors, such as transcription factors). However, a gene may harbor multiple SNPs, and while some variants may exert measurable effects, their cumulative impact on gene expression could be minimal, as the effects of individual SNPs may counteract one another. Therefore, the distribution of SNPs across functional categories, including transcription factors (TFs) and sigma factors (SFs), does not by itself provide clear biological insight without large-scale statistical power, typically achieved through very large population-level studies encompassing hundreds of strains. Such one-to-one association studies between genetic variants and molecular phenotypes clearly go beyond the scope of the present work, as the reviewer also noted.\nNevertheless, we examined our comparative genomic data between L1 and L2 strains, focusing on SNPs that are conserved within each lineage but differ between them (i.e., lineage-specific SNPs). We then analyzed the distribution of these SNPs across functional categories. Interestingly, both TFs and SFs were among the categories with relatively few nonsynonymous SNPs compared to other functional groups (Table 1). However, for reasons of thematic coherence, we feel that including this analysis may divert attention from the central objectives of the manuscript. We therefore prefer not to incorporate it in the main text but would be happy to provide the data upon request.\nCategory Gene Syn. Nonsyn. Syn. ratio Nonsyn. ratio\nVirulence, detoxification, adaptation 239 8 14 0.0335 0.0585\nSigma factor 13 1 0.0769\nTranscription Factor 214 16 20 0.0748 0.0934\nRegulatory proteins 198 19 24 0.0960 0.1212\nInformation pathways 242 27 31 0.1116 0.1280\nIntermediary metabolism and respiration 936 79 134 0.0844 0.1431\nCell wall and cell processes 772 73 118 0.0945 0.1528\nEssential Genes - Growth - In vitro 617 62 97 0.1005 0.1572\nmycolic acid related protein 62 4 12 0.0645 0.1935\nLipid metabolism 272 29 57 0.1066 0.2095\nTable 1. Distribution of lineage-specific SNPs across functional categories. The table summarizes the number of lineage-specific SNPs identified between MTBC lineages 1 and 2 across different functional gene categories. Ratios represent the number of synonymous (Syn.) and nonsynonymous (Nonsyn.) SNPs normalized by the number of genes within each category.\nSecondly, it is important to note that these two lineages represent different levels of genetic diversity, with L1 exhibiting greater genomic variation than L2. To provide broader evolutionary context, we extracted pN/pS values from the recently published paper by Sarah M. Fortune and colleagues (PMID: 41033311), which includes ~55,000 clinical isolates spanning lineages 1–4. We specifically evaluated TFs and SFs within this dataset.\nThe ratio of nonsynonymous to synonymous polymorphisms (pN/pS) reflects the selective pressure acting on protein-coding genes, where pN/pS < 1 typically indicates purifying (negative) selection. As summarized in the violin plots below (Fig. 1), both TFs and SFs exhibit pN/pS distributions comparable to the genome-wide averages across all lineages, suggesting that these gene classes are under similar purifying selection rather than lineage-specific adaptive diversification. In L2, slightly higher median and mean pN/pS values were observed for both TFs and SFs, but this trend parallels the overall distribution across all genes, indicating that some genes, including TFs and SFs, may experience relatively weaker purifying selection. While these results are of evolutionary interest, we believe they are not central to the main focus of this study and could potentially be misinterpreted if overemphasized.\nFigure 1. Distribution of pN/pS ratios across MTBC lineages and functional categories. Violin plots show the distribution of the ratio of nonsynonymous to synonymous polymorphisms (pN/pS) across\n12\n~55,000 MTBC strains. Lower pN/pS ratios indicate stronger purifying (negative) selection. Each panel represents a distinct gene subset: (A) all coding genes, (B) transcription factors, and (C) sigma factors. Within each panel, the first violin corresponds to all strains combined, followed by lineage-specific distributions (L1–L4).\nFinally, regarding the transcriptional model, we emphasize that it incorporated all differentially regulated genes, not only TFs. The TFs discussed in the manuscript were identified as outputs of the model. The main objective was to assess the extent to which transcriptional variation between L1 and L2 strains could be explained by TF-mediated regulatory control. Please see “Regulated target genes” in Fig. 5A, which illustrates how many genes differentially regulated between L1 and L2 strains are regulated by the corresponding TFs identified through our transcriptional network analysis.\n11) Results: looking at Table S1, it appears that strain 0153 is primarily driving the large expression differences observed for Rv1985c between L1 and L2. Given that some L1 strains from the same repository have been reported to harbor a deletion of Rv1985c (doi.org/10.1038/s41598-018-22237-5), it would be important to verify the presence or absence of this gene across all isolates included in the analysis. Notably, strain 0153 also displays higher diversity on the PCA for protein data, as well as in the intra-lineage mRNA expression profiles. Interestingly, the magnitude of its differential expression (|log2FC| approx. 10) is comparable both within its own genotype and in inter-genotype comparisons, suggesting that the expression pattern is strain-specific rather than genotype-driven. This reinforces the need to confirm that gene Rv1985c is present and intact in all strains before interpreting genotype-level effects.\n12) The possibility that gene Rv1985c is deleted in some L1 isolates should be addressed in the discussion. While the Authors conclude that it functions as a validated master regulator distinguishing lineages L2 and L1, the presence of a deletion in some L1 isolates challenges this interpretation and suggests that the observed differential expression may not uniformly reflect a lineage-level effect.\nWe thank the reviewer for these closely related comments (11 and 12), which provided a great opportunity to further clarify and strengthen the interpretation of our findings.\nWe confirm that the Rv1985c gene is absent in strain N0153 (L1), consistent with previously reported deletions in certain L1 isolates. However, we maintain that interpreting Rv1985c as a lineage-associated transcriptional regulator, rather than a purely strain-specific effect, remains valid for several reasons.\nFirstly, the identification of Rv1985c originated from our transcriptional network analysis, in which four transcription factors, including Rv1985c, collectively explained ~25% of the differentially expressed genes between L1 and L2 strains. This network analysis was based on experimentally validated TF-target relationships derived from ChIP seq experiments in H37Rv (PMID: 25581030). Importantly, the model does not select only TFs that are themselves\n13\ndifferentially expressed. Instead, it identifies TFs whose direct downstream targets exhibit significant differential regulation between lineages. Consequently, even if a TF’s own expression does not change, the model can detect its regulatory footprint if its target genes collectively show significant differential expression. In this context, the size of the target set, the subset of differentially expressed targets, and the overall number of significant changes are all relevant.\nFor example, Rv3597c emerged as a major TF in our analysis despite not being differentially expressed between lineages (Fig. 5B). This aligns with established evidence that many bacterial TFs undergo post-translational modifications or other mechanisms that modulate activity independently of transcript abundance. Therefore, the identification of Rv1985c as a key TF is not invalidated by its absence in a single L1 strain; rather, it reflects lineage-level transcriptional differences that are consistent with its regulatory network footprint.\nTFs Degree L2 / N0052 / N0052 / N0052 / N0145 / N0145 / N0145 / N0155 / N0155 / N0155 /\nL1 N0072 N0153 N0157 N0072 N0153 N0157 N0072 N0153 N0157\nRv0275c 10 1 1 1 1 1 1 1 1 1 1\nRv2017 10 1 1 1 1 1 1 1 1 1 1\nRv2506 10 1 1 1 1 1 1 1 1 1 1\nRv3183 10 1 1 1 1 1 1 1 1 1 1\nRv1985c 8 1 1 1 1 1 1 1 1\nRv2359 8 1 1 1 1 1 1 1 1\nRv3164c 8 1 1 1 1 1 1 1 1\nTable 2. Consistency of transcription factor (TF) identification across all pairwise lineage comparisons. The table summarizes the top TFs identified by network-based transcriptional analysis across all nine possible pairwise comparisons between L1 and L2 strains. “Degree” indicates the number of times a given TF is identified among all nine pairwise comparisons, in addition to the overall L2/L1 comparison.\nSecondly, to rigorously assess whether the signal for Rv1985c reflects a lineage-level or strain-specific effect, we repeated the network-based transcriptional analysis across all nine pairwise comparisons between L1 and L2 strains (3 × 3 combinations). In six of these comparisons, N0153 (the Rv1985c-absent strain) was excluded. The top TFs consistently identified across multiple comparisons are summarized in Table 2 and Fig. 2 (below). Notably, Rv1985c remained among the most consistently identified regulators even when N0153 was excluded, supporting that the observed signal reflects a lineage-level effect rather than being driven solely by a single strain. Interestingly, Rv1985c was not significantly identified in the N0155/N0153 comparison, despite being largely upregulated in N0155 relative to N0153 (Log FC = 10.33, adjusted P-value = 4.22 × 10⁻⁷²), which further illustrates the principles of the\n2\nnetwork-based analysis.\nThirdly, we assessed the expression of Rv1985c using TPM (transcripts per million) data, excluding N0153. Mean TPM values per strain (across technical replicates) were 37.17 (N0052), 43.41 (N0145), and 45.93 (N0155) in L2 strains, and 32.01 (N0072) and 34.32 (N0157) in L1 strains. Based on the median of these strain-level means, expression in L2 strains was ~1.4-fold higher than in L1 strains.\nOverall, these analyses further support the conclusion that Rv1985c represents a lineage-level transcriptional regulator, and we have discussed this point in the relevant Results and Discussion sections of the manuscript (Lines 312-319, 498-504).\nFigure 2. Overlaps of transcription factors (TFs) identified across L2 and L1 strains. UpSet plot illustrating the number of TFs shared between different strain-specific L2 vs L1 comparisons (3 L2 strains × 3 L1 strains), as well as TFs identified in the overall L2 vs L1 comparison.\n13) Please consider reporting BCCM/ITM IDs for all strains to avoid misunderstanding.\nWe thank the reviewer for drawing our attention to this issue. We have now added the BCCM/ITM strain identifiers to the corresponding Methods section to avoid any potential misunderstanding. Of note, the full description of the used clinical strains has been discussed in following paper as already cited (PMID: 30908506) (Lines 545-546)\nL1 strains:\n15\n N0072: ITM-500942; CT2018-00083\n N0153: ITM-500943; CT2018-00084\n N0157: ITM-500941; CT2018-00082\nL2 strains:\n N0052: ITM-601826; CT2018-03241\n N0145: ITM-500944; CT2018-00085\n N0155: ITM-500946; CT2018-00088\n14) Discussion: please consider including a “limitations of the study” paragraph\nWe thank the reviewer for this suggestion. We have now added several sentences to highlight the limitations of the study in the Discussion section.\n16"}]

Decision letter