M. bovis AF2122/97 and M. bovis BCG-Russia (BCG-1) strains were grown in at 37 °C in Middlebrook 7H9 broth (Difco) containing 10% albumin- dextrose-catalase (ADC), 0.5% Pyruvate, and 0.05% Tween 80. To define the distinct in vitro growth phases, we generated standard growth curves for both M. bovis and BCG Russia by measuring CFU/mL over time at 24-hour intervals ( Supplementary Figures S1A, B ). This allowed consistent and biologically validated sampling at early-log, mid-log, and stationary phase.

For RNA preparation, bacterial cultures at OD600 0.5, 1 and 2 were snap-frozen on ice and centrifuged at 3,000 × g for 10 min at 4°C. Bacterial pellets were then stored at -80°C. Bacterial stocks for mouse infections were prepared as previously described (AL et al., 1997).

C3HeB/FeJ mice groups (N = 20/each) at 5–6 weeks age was infected with approximately 100 CFU per mouse were administered by aerosol using the Glas-Col inhalation system (Glas-Col, LLC, Terre Haute, IN) as outlined before (Abdelaal et al., 2019). The infectious dose for each group was confirmed by plating lungs of an infected mouse at 1-day post- challenge. Mice were sacrificed at 4- and 16-weeks’ post infection for both histopathology bacterial CFU enumeration as detailed before (Marcus et al., 2016; Abdelaal et al., 2019).

RNA was isolated form frozen cultures using a TRIzol based method described previously (Abomoelak et al., 2009; Ward et al., 2010a; Ward et al., 2010b). RNA from murine lung tissues was isolated by homogenizing in TRIzol Reagent (Molecular Research Center, Cincinnati) and centrifuging at 3,000 × g for 5 min at 4 °C. Bacterial pellets were resuspended in 0.5 ml TRIzol Reagent containing 1% polyacryl carrier (Molecular Research Center), transferred to screw-caped tubes with 0.25 ml of zirconia/silica beads, and broken in a bead beater (Biospec Products, Bartlesville, OK). RNA was isolated by using TRIzol per the manufacturer’s instructions and treated with DNase I (DNA-free kit, Ambion). RNA was treated with TURBO DNase until PCR negative to remove contaminating DNA.

RNA-seq data analysis was performed on the CLC Genomics Workstation 8.0. Sequence reads were aligned to the Mycobacterium bovis AF2122/97 parental reference genome (GenBank accession number NC_002755). Raw sequencing reads underwent quality control using FastQC, and adapter trimming was performed with Trimmomatic to remove low-quality bases (Bolger et al., 2014). Reads were then aligned to the reference genome using HISAT2 (Kim et al., 2015). The reads per kilobase per million (RPKM) value for each gene was generated. Normalization was conducted using the median of ratios method within DESeq2 to account for sequencing depth and RNA composition biases (Love et al., 2014). Differential gene expression analysis using the R DEseq2 package was performed for the following groups of data sets for OD600 of 0.5, 1.0 and 2.0. The exact-test function was applied to determine the association of the differences in expression read counts within each group, and corresponding P-values were adjusted using the default Benjamini & Hochberg procedure (Benjamini et al., 2001; Reiner et al., 2003). Their adjusted P-values, in -log10 scale on the y axis and fold changes in log2 scale on the x axis, were plotted as a volcano plot. Differential gene expression was determined by a false discovery rate (FDR) threshold of p < 0.05. Genes with read counts of less than 5 were eliminated.

To identify gene co-expression modules associated with biological processes in M. bovis, we performed Weighted Gene Co-expression Network Analysis (WGCNA) using normalized gene expression data (Zhang and Horvath, 2005). Prior to network construction, lowly expressed genes were filtered, and variance-stabilized counts were used as input. A soft-thresholding power = 12 was selected using the scale-free topology criterion, ensuring optimal network fit and connectivity. A signed adjacency matrix was computed and transformed into a Topological Overlap Matrix (TOM), followed by hierarchical clustering to define co-expression modules. Gene module “eigengenes” were assigned distinct colors, and eigengenes were calculated to summarize module expression patterns. To investigate biological relevance, module-trait correlations were computed using sample metadata, identifying key modules associated with experimental conditions. Gene ontology (GO) enrichment and pathway analyses were performed using ClusterProfiler to assess functional relevance (Yu et al., 2012). Hub genes were determined based on intramodular connectivity, identifying core regulators within each module.

Inter-module relationships were visualized through eigengene expression analysis, highlighting shared and unique regulatory patterns. To enhance reproducibility, all analyses were performed in R using the DESeq2, WGCNA, and ClusterProfiler packages. This framework provided a robust systems-level approach for uncovering transcriptional signatures linked to growth-phase dynamics and virulence-associated differences in M. bovis strains.

Primers with similar melting temperatures (60–66°C) were designed by using PRIMER 3 software (Rozen and Skaletsky, 2000). The sequences of primers are available on Supplementary Table S1 . All primers were tested in PCRs with 100 AF2122/97 genome equivalents as template and the amplification products were evaluated by gel electrophoresis. A total of 1-2 μg RNA were used as template for cDNA using Superscript III. A SYBR green based qRT-PCR protocol utilizing GoTaq^® qPCR Master Mix (Promega, Madison, WI) and the Step One Plus TM Real-Time PCR System (Applied Biosystems^®, Foster City, CA) were used. For in vitro samples, gene expression was normalized to 16S rRNA (rrs). For in vivo samples, expression was normalized to sigA (BQ2027_MB0760), a stable internal control during mycobacterial infection in mice. Analysis was carried out using LinRegPCR (Ruijter et al., 2009). Two biological replicates with no less than three technical replicates each were completed. Gene expression levels were determined using the ΔΔCt method (Livak and Schmittgen, 2001) for murine lungs relative to in vitro culture at mid-log phase. Statistical significance between groups was assessed using a one-way ANOVA with Tukey’s post-hoc test for multiple comparisons or an unpaired Student’s t-test for two-group comparisons.

Statistical analyses for RT-qPCR, bacterial enumeration, and growth curve assays were conducted using GraphPad Prism 9.4.1 (GraphPad Software, San Diego, CA, USA). Differences between groups were evaluated using a one-way analysis of variance (ANOVA) with Dunnett’s post-test to control for multiple comparisons. Growth curves were fitted using a nonlinear regression model to determine key growth parameters, including lag phase duration, maximum growth rate, and stationary phase plateau.

All statistical analyses were performed under parametric assumptions, with normality and homogeneity of variance assessed using the Shapiro-Wilk test and Levene’s test, respectively. A threshold of p < 0.05 was considered statistically significant. **Significant differences are labeled accordingly in the figures as * p  < 0.05, ** p  < 0.01, *** p  < 0.001, and **p < 0.0001.

To profile M. bovis transcriptional machinery, we employed an RNA-Seq approach to identify genes associated with the growth of the virulent M. bovis AF2122/97 strain and the attenuated BCG-Russia strain. To capture the complete transcriptional landscape across different growth phases, cultures of M. bovis and M. bovis BCG (BCG) were harvested at early log (OD₆₀₀ = 0.5), mid-log (OD₆₀₀ = 1.0), and stationary phases (OD₆₀₀ = 2.0) Growth curves for both strains are shown in Supplementary Figure S1B . This design enabled us to assess the dynamic transcriptional changes associated with bacterial growth at early-log, mid-log and stationary phases. The RNA-Seq analysis generated an average of 138 ± 5.2 million paired-end reads per library, exceeding previously defined quality control criteria for sequencing depth and alignment (Tarazona et al., 2011; Castel et al., 2015). On average, 91 ± 1.3% of reads mapped uniquely to the M. bovis AF2122/97 reference genome (GenBank accession number NC_002755). The RNA-Seq datasets detected transcripts for 64–66% of the coding regions in each sample, providing robust transcriptome coverage for 98.2% of the predicted genes encoded in M. bovis under all examined growth phases. The sample distance matrix showed clear clustering by biological replicate ( Supplementary Figure 2A ), and the library size distributions across samples ( Supplementary Figure 2B ). Together, these metrics ensured robust differential gene expression analysis. The summary statistics of the RNA-Seq data for each replicate are summarized in Table 1 .

To analyze the transcriptomic differences between the virulent and the attenuated bovine tubercle bacilli, we employed direct pairwise comparisons of the transcriptome of both organisms using a P-value threshold of < 0.05 and > ± 1.5-fold change. When early log (OD₆₀₀ = 0.5) versus mid-log (OD₆₀₀ = 1.0) cultures were compared in M. bovis, significant downregulation was observed for genes (n = 105) such as fadE23, fadE24, and fadE5 ( Figure 1A ), essential components of lipid metabolism pathways (Crick and Guan, 2016), suggesting metabolic reprogramming during growth progression. Additionally, the downregulation of cydB and ndh highlights adjustments in respiratory pathways (Mittal et al., 2018). In contrast, BCG at the same transition exhibited significant upregulation of BQ2027_MB1086 and scoA among the 126 upregulated and 95 downregulated genes ( Figure 1A ), pointing toward enhanced central metabolic activity in the attenuated strain (Houben et al., 2006).

When early log (OD₆₀₀ = 0.5) versus stationary phase (OD₆₀₀ = 2.0) cultures were compared in M. bovis, consistent downregulation of genes such as BQ2027_MB1650c, cydB, and cydA was observed, indicating shifts in respiratory and electron transport processes ( Figure 1B ) (Vaziri and Brosch, 2019). Upregulation of yrbE2A, fusA2b, and rnj suggests adaptive mechanisms involving nutrient acquisition and stress responses during exponential growth ( Figure 1B ) (Seshadri et al., 2009; Forrellad et al., 2013; Wang et al., 2022; Silva-Pereira et al., 2024). Similarly, BCG exhibited distinct metabolic changes with the upregulation of genes such as fadE22a, BQ2027_MB2662, and PE12, reflecting enhanced stress response mechanisms. Genes related to lipid metabolism, such as mce1D and espc, were significantly downregulated in BCG at this growth stage, consistent with its attenuated phenotype ( Figure 1E ) (Houben et al., 2006).

Finally, when comparing mid-log (OD₆₀₀ = 1.0) versus stationary (OD₆₀₀ = 2.0) phases, M. bovis exhibited upregulation of (n = 162) including dormancy- and survival-related genes, such as fadE23, fadE24, and BQ2027_MB0909, underscores its preparation for stationary phase and environmental adaptation ( Figure 1C ) (Houben et al., 2006). Concurrently, downregulation of genes (n = 189) such as BQ2027_MB2981 and BQ2027_MB1023 highlights the metabolic shifts favoring survival over replication (Vaziri and Brosch, 2019). In contrast, BCG demonstrated upregulation of stress response genes, such as desa3, desA1, and modA ( Figure 1F ) (Chang and Fox, 2006), while showing downregulation of genes like ino1, suggesting a diminished reliance on virulence-associated pathways (Movahedzadeh et al., 2004).

To further explore the transcriptional divergence of M. bovis AF2122/97 and M. bovis BCG-Russia, we conducted direct pairwise comparisons between both strains growing at different growth phases. This analysis reveals phase-specific expression patterns that define the M. bovis virulent and attenuated strains ( Figure 2 ).

At early growth phase, M. bovis significantly upregulates (n = 326) including key virulence-associated genes, such as eccc2_2 (encoding ESX-II secretion-associated protein EccC2) and mmsA (encoding methylmalonate-semialdehyde dehydrogenase) ( Figure 2A ) (Gibson et al., 2022). The level of changes in gene expression ranged between log₂FC = 5.79 to log₂FC = 8.73 respectively. These Differentially Expressed Genes (DEGs) are integral to pathogenic mechanisms, including secretion systems and metabolic processes crucial for early-stage infection (Gibson et al., 2022). Although eccC2_2 expression was higher in M. bovis than in BCG, canonical ESX−1 genes (e.g., esxA/B, espA/C) were undetectable in BCG due to the stable deletion of the RD1 locus known to underlie its attenuation. The expression differences in eccC2_2 and mmsA do not themselves imply virulence roles but rather reflect strain-dependent variation in secretion-system activity and lipid metabolism.

At mid-log growth phase, the transcriptional landscape further differentiates (n = 374), with M. bovis continuing to show elevated expression of mmsA (MB0775c; methylmalonate−semialdehyde dehydrogenase), BQ2027_MB3287 (a putative metallopeptidase family protein), and BQ2027_MB0467 (exaC; NAD⁺−dependent acetaldehyde dehydrogenase), with log2 fold changes exceeding 6. Their regulation is consistent with known virulence strategies: metallopeptidases like Zmp1 are implicated in macrophage inflammasome suppression, and aldehyde dehydrogenation is required for detoxifying host-derived reactive aldehydes during persistent infection. Consequently, these genes are better framed as components of metabolic adaptation layers indirectly supportive of virulence (Guo et al., 2023). Meanwhile, M. bovis BCG demonstrates upregulation of stress-response and metabolic genes, indicative of its focus on environmental adaptability rather than virulence ( Figure 2B ).

At stationary phase, the virulent strain adapts to long-term survival, with pronounced upregulation of genes (n = 378) like tkt (encoding transketolase) and apa (encoding alanine-proline-rich antigen), highlighting a focus on dormancy and immune modulation (Sable et al., 2011; Fullam et al., 2012). In contrast, M. bovis BCG continues to prioritize stress response pathways, with limited expression of virulence-associated factors ( Figure 2C ). These findings illustrate the phase-specific transcriptional shifts underlying the divergent phenotypes of M. bovis and M. bovis BCG, providing insights into the molecular determinants of virulence and attenuation.

To further examine the overlap in differential gene expression across growth phases, we analyzed the intersection of differentially expressed genes between early, mid-log, and stationary phases for both M. bovis and M. bovis BCG ( Figure 2D ). The UpSet plot illustrates the extent of shared and unique transcriptional responses across growth transitions. Notably, M. bovis exhibits a greater number of overlapping DE genes between OD600 0.5 and OD600 1.0 (n = 31), suggesting a coordinated metabolic shift between early and mid-log phases. Conversely, M. bovis BCG shows fewer shared DE genes across phases, indicating a more gradual transcriptional adaptation with a focus on stress response mechanisms. The relatively high number of distinct genes (n = 42) differentially expressed only in the stationary phase in M. bovis further underscores its ability to enter a dormancy-like state, a feature less pronounced in M. bovis BCG. This comparative transcriptomic analysis highlights key regulatory mechanisms that distinguish virulence-driven metabolic shifts in M. bovis from the attenuation-associated transcriptional changes in BCG.

To better characterize gene underlying the M. bovis adaptation at different growth phases, we employed Weighted Gene Co-expression Network Analysis (WGCNA) to identify co-expressed gene networks and their relationships with experimental traits ( Figure 3 ). This analysis revealed 13 distinct gene modules ( Supplementary Table S2 , Supplementary Figure 3 ) that show correlation with bacterial growth phase progression, reinforcing the phase-specific regulatory mechanisms governing virulence and metabolic adaptation.

The module-trait correlation analysis ( Figure 3A ) revealed that five modules (turquoise, blue, brown, yellow, and greenyellow) exhibited significant correlations with bacterial growth phase progression. While not all modules exceeded an R-value of 0.6, their eigengenes expression (co-expressed genes based on WGCNA) trends ( Figure 3B ) aligned with observed DEGs and strain-specific transcriptomic behavior, prompting their inclusion in subsequent functional analyses.

The turquoise module I (1,626 genes) exhibited the strongest correlation with bacterial growth phase (r = 1.0, p < 0.01), suggesting its role in adaptive responses to metabolic shifts during transition from exponential to stationary phase. Genes within this module included several lipid metabolism regulators, such as fadD21, fadD22, fadD23, fadD25, fadD26, fadD28, fadD29, and fadD30, which are involved in fatty acid degradation and lipid utilization ( Figures 3A, B, J ), processes essential for energy metabolism and adaptation to stationary phase (Simeone et al., 2010). Additionally, the presence of yrbE2A, fusA2b, rnj, and mmpL3 suggests that this module also contributes to membrane transport and stress adaptation (Su et al., 2021), enabling bacterial survival under nutrient-limited condition which is considered key factors for mycobacterial virulence during infection.

The blue module II, (695 genes), exhibited a strong negative correlation with the bacterial strain factor (r = -0.40, p < 0.05), indicating an inverse relationship—higher expression in M. bovis relative to BCG—which is a valid and interpretable outcome within the WGCNA framework (Langfelder and Horvath, 2008). This module was enriched with genes associated with virulence and host-pathogen interactions. Notably, eccC2_2 and mmsA, key components of the ESX-1 secretion system and methylmalonate metabolism, respectively, were strongly upregulated in M. bovis compared to M. bovis BCG ( Figures 3A, D, G ) (Cole et al., 1998). The presence of PPE and PE_PGRS family genes further suggests that this module plays a crucial role in immune evasion and intracellular persistence, consistent with the virulent phenotype of M. bovis (Domenech et al., 2005).

Conversely, the brown module III, (466 genes), which exhibited a moderate correlation with growth phase (r = 0.35, p < 0.05), contained genes linked to central metabolism and oxidative stress response, including scoA, BQ2027_MB1086, and fadE22a ( Figures 3A, C, H ) (Bitter et al., 2009). The upregulation of these genes in M. bovis BCG suggests an enhanced reliance on metabolic pathways that compensate for the attenuation of virulence-associated factors (Gey van Pittius et al., 2006).

The yellow module IV, (193 genes), displayed a strong association with stationary phase adaptation (r = 0.73, p < 0.01) and contained genes such as bfrB and desa2, which are known to be involved in oxidative stress protection and lipid homeostasis ( Figures 3A, F ) (Geiman et al., 2006). The presence of these genes in this module suggests an essential role in long-term survival strategies, particularly for M. bovis BCG, which exhibits increased dormancy-associated gene expression (Brosch et al., 2007).

The green-yellow module V, (49 genes), was enriched with genes involved in cell wall remodeling and transport processes, such as mce1D, espc, and yrbE2A, which were significantly upregulated in M. bovis compared to BCG ( Figures 3A, F, I ) (Reddy et al., 2012). This module showed a moderate negative correlation with bacterial strain (r =-0.74, p = 0.03), indicating strain-specific expression differences rather than direct regulation by bacterial growth phase. The higher expression of these genes in M. bovis suggests a potential role in cellular integrity and nutrient uptake during early exponential growth. However, their expression declines over time, indicating that they are primarily utilized in the early stages of growth rather than throughout bacterial adaptation (Boon and Dick, 2002).

Overall, these findings provide a comprehensive view of the transcriptional programs governing M. bovis and M. bovis BCG physiology. Among the 13 identified co-expression modules, five were strongly associated with virulence, metabolic adaptation, and growth phase progression; turquoise, blue, brown, yellow, and green-yellow. These modules encapsulate key gene networks related to lipid metabolism, immune evasion, oxidative stress response, dormancy, and cell wall remodeling. Together, they highlight distinct strategies employed by the virulent and attenuated strains to adapt to environmental shifts and host-related pressures.

M. bovis encodes nearly 200 transcriptional factors, similar to M. tuberculosis (Smith, 2003). The differential gene expression observed between the virulent and the attenuated bovine tubercle bacilli may be a consequence of differences in global transcriptional regulators between the two species. To address this hypothesis, a curated transcription factor enrichment analysis was performed and revealed the significant association of 10 transcription factors (Turkarslan et al., 2015) ( Figure. 4A ). The transcription factors included in our analysis were selected based on known regulatory networks within the MTBC, and their differential expression (list fold of change used here) observed in our dataset and their known associations with virulence. This targeted approach allowed us to focus on the most relevant transcription factors impacting pathogenicity.

The association of transcription factors such as alternate sigma factors SigK and SigF along with cytoplasmic redox sensor WhiB3 with the differentially expressed gene (DEG) lists indicates that disparate expression of virulence- relate pathways regulated by these transcription factors between the two pathogens could have significant consequence for infection (Gebhard et al., 2008; Veyrier et al., 2008). Importantly, PhoP and EspR were significantly differentially expressed, these transcription factors are important for adaptation of M. bovis to the intracellular environment and are functionally linked by such processes (Sherman et al., 2001; Singh et al., 2009). PhoP and EspR regulate the expression of ESX-1 secretion system-related genes. Furthermore, PhoP was expressed to a higher level in BCG; this may represent an attempt at a compensatory mechanism for aberrant PhoP signaling and supports previous reports of suboptimal PhoP signaling in M. bovis (Gonzalo Asensio et al., 2006; Walters et al., 2006; Gonzalo-Asensio et al., 2014). Additionally, 8 of 48 canonical DosR regulon, namely BQ2027_MB1632c, BQ2027_MB1633c, BQ2027_MB1634c, BQ2027_MB1636c, BQ2027_MB1650, BQ2027_MB1653, BQ2027_MB2030, and BQ2027_MB2033 (Sivaramakrishnan and de Montellano, 2013; Belardinelli et al., 2025) were significantly upregulated in M. bovis virulent strain, consistent with the role of DosR in early adaptation to hypoxic or stress conditions. This differential expression pattern may reflect the integration of DosR responses with other transcriptional regulators such as PhoP, EspR, and WhiB3 (Gonzalo Asensio et al., 2006; Gonzalo-Asensio et al., 2014). DosR regulon is important during conditions that do not allow aerobic respiration (Leistikow et al., 2010), like in the lungs of M. tb-infected mice (Voskuil et al., 2003) and in interferon-gamma-activated murine macrophages (Schnappinger et al., 2003).

Interestingly, several genes that were in the deleted regions from BCG were highly expressed in M. bovis ( Supplementary Figure 4 ) including the major antigens ESAT-6 and CFP10, secreted by the ESX-1 secretion system of the MTBC, a system which has been implicated in mycobacterial escape from the phagosome to the cytosol that results in a Type-I interferon response within the infected macrophage (Simeone et al., 2009; Simeone et al., 2012; Simeone et al., 2015a). As EspR is induced in M. bovis AF2122/97, there is a significant induction of the ESX-1 secretion system in M. bovis, including esx-1-related proteins such as esxA, espA, espC, espD (Raghavan et al., 2008; Simeone et al., 2009; Chen et al., 2013; Solans et al., 2014; Cao et al., 2015) ( Figure 4B ). Additionally, all genes related to RD1 region ( Figures 4B, C ) was observed to be diminished expression in M. bovis BCG-Russia, this is an emphasis on the identity of each strain as RD1 is absent from the vaccinal strain. Alternate transcriptional regulation between the M. bovis BCG-Russia and the M. bovis AF2122/97 may represent differential priming events in preparation for the initial interactions of both species with their respective host immune systems. Increased expression of the ESX-1 secretion system may facilitate faster escape of M. bovis AF2122/97 from the phagosome into the cytosol in contrast to BCG, hence triggering DNA-sensing pathways and increased IFN response seen in our data (Simeone et al., 2012; Simeone et al., 2015b).

The results of transcriptomic analysis indicated subtle difference in lipid metabolism related genes between the virulent and the attenuated bovine tubercle bacilli. M. bovis AF2122/97 showed higher expression of the Pks13/FadD32 pair ( Figure 4D ), which are involved in the Biosynthesis of Mycolic Acids (Gavalda et al., 2009). FadD32 gene is adjacent to pks13, this genetic loci, fadD32-pks13, is conserved also in M. tb and M. leprae (Marrakchi et al., 2014), as this gene cluster is restricted to mycolic-acid-producing bacterial species (Gande et al., 2004). FadD32, a fatty acyl- AMP ligase, is involved in catalyzing the formation of acyl-adenylates, the activated form of meromycolic acid substrate in the mycolic condensation reaction (Trivedi et al., 2004). It also assists the transfer of the meromycoloyl chain onto the N-terminal acyl carrier protein (ACP) domain of the condensing enzyme Pks13 (Le et al., 2016). While Pks13 is a unique polyketide synthase (PKS) forms the a-alkyl β-ketoesters which is the direct precursors of mycolic acids (Gavalda et al., 2009). Importantly, Polyketide synthase Pks13 and its acyl−AMP ligase partner FadD32, encoded within the same locus, are universally regarded as essential for mycolic acid biosynthesis and bacterial viability in the M. tb complex (Portevin et al., 2004). Also, Pks13/FadD32 pair have been shown to be required for virulence in M. tb (Sassetti and Rubin, 2003; Mukhopadhyay et al., 2012). The M. bovis virulent strain orthologues may also, therefore, play a role in virulence. Upregulation of Pks13/FadD32 pair (log₂ fold change ranging from 2.3 to 4.3 across growth phases) may indicate divergent expression of Mycolic acid between the virulent and attenuated strains of M. bovis, an observation that can be reflected on the composition of the cell wall of each bacillus.

MmaA4 was found also to be induced in the virulent M. bovis. MmaA4, a hydroxy- mycolate synthase, is also involved in mycolic acid modification by converting it to hydroxy mycolic acid, a precursor of methoxy‐ or keto‐mycolic acid (Alahari et al., 2009). It was found that MmaA4 modulates IL-12 production. MmaA4 knockout mutant induced more IL-12 from murine macrophages and were attenuated for virulence in mice (Dao et al., 2008). As a result, the ΔmmaA4 mutant strain induced significantly elevated levels of this critical Th1-type cytokine in macrophage cultures. Additionally, MmaA4 knockout mutant of BCG vaccine induced higher levels of mycobacterial-specific multifunctional T cells, is more protective than BCG vaccine, and, surprisingly, may be safer than BCG when used in immunocompromised animals (Derrick et al., 2012; Derrick et al., 2016). Although BCG-Russia exhibited overall lower expression of certain lipid metabolic enzymes, the functional consequences of these differences require further investigation.

Interestingly, BQ2027_Mb2982c, encoding a glycosyltransferase involved in the synthesis of the trisaccharide phenolic glycolipid (PGL) that is derived from phthiocerol dimycocerosates (PDIM), was found to be expressed at higher levels in M. bovis BCG-Russia (data not shown). However, it should also be noted that BQ2027_Mb2982c is non-functional in M. bovis, so the higher expression of its gene may simply be due to loss of negative feedback inhibition. Previously, it was found that loss of PIDM/PGL reduces the protective efficacy of BCG vaccine (Tran et al., 2016). Since the loss of PDIM and PGL occurs naturally in a subset of BCG strains (Chen et al., 2007), it also suggests that these strains may have been over-attenuated, which compromises their effectiveness. Overall, most of these genes are participated in energy metabolism, including the fatty acid, cholesterol and glycolipid metabolism. It is postulated that pathogens downregulated its metabolic activity to reduce energy consumption and to persist in a prolonged dormant state.

As we were able to identify in vitro differentially expressed genes between attenuated and virulent M. bovis strains. To validate the biological relevance of key genes identified in vitro, we analyzed gene expression of a limited set of genes in a susceptible murine model infected with the virulent strain. Lung tissues were analyzed using qRT-PCR, a sensitive assay for gene expression, to confirm their expression during in vivo infection. Our analysis focused on key in vitro expressed differentially expressed genes change their expression during the infection of the host tissue using the C3HeB/FeJ mouse model. We chose this model as it recapitulates the hallmark of bovine tuberculosis lung lesions following M. bovis aerosol infection (Boute et al., 2017). Mouse groups were sampled at 4- and 16- weeks post-infection to represent early and progressive stages of infection following aerosol infection with M. bovis AF2122/97. colonization levels increased at 4 w.p.i but peaked by 16 w.p.i ( Figure 5a ). Histologically, type I lesions resulted from the occlusion of alveolar spaces by a cellular infiltrate were noticeable by 4 w.p.i ( Figure 5b ). However, previous studies of M. bovis infection in C3HeB/FeJ murine model showed that type I lesions that evolved into an organized granuloma with a central accumulation of foamy macrophages were only visible by 5 w.p.i (Boute et al., 2017). This might explain the less organized lesion observed in the examined mouse lungs. As infection progressed, the inflammatory responses were intensified in lungs by 16 w.p.i with the observation of central necrosis in type I lesions ( Figure 5c ). These histological findings agree with the increase of bacterial burden by progression of infection.

Several transcripts were selected to be quantified using quantitative RT-PCR ( Figure 6 ), from bacterial RNA purified from the murine lung tissue based on their differential expression magnitude, known virulence associations, relevance to secretion or lipid metabolism pathways and their potential involvement in adaptation of M. bovis into host microenvironment a as suggested previously (Brosch et al., 2007; Malone et al., 2018). Among the highly regulated genes Mb3614c (BQ2027_MB3614c), pstS3 (BQ2027_MB0951), PPE40 (BQ2027_MB2377c), fbpB (BQ2027_MB1918c), whiB6 (BQ2027_MB3892c) and espR (BQ2027_MB3910c). Interestingly, most of the genes whose expression dramatically changed between virulent and attenuated strains during in vitro culture showed similar differential expression during lung tissue infection ( Figure 5 ). For example, Mb3614c is a putative transcription factor reported to play a regulatory role under starvation conditions (Ramos et al., 2020), was found to be also induced during early stage of lung tissue infection. Similarly, PPE40, a hypothetical PPE-family protein predicted to be an outer membrane protein, which is a part of the ESX 5 secretion system only found in pathogenic slow growing mycobacteria (Baena et al., 2019). It was found that PPE40 was induced more at early phase of infection. While pstS3 (BQ2027_MB0951) a periplasmic phosphate-binding lipoprotein (Garnier et al., 2003) was found to be induced to more extent at earlier stages of infection, while downregulated in In-Vitro culture. The pstS3 is a known as a component of this primary phosphate uptake system and reported previously to be highly expressed in mouse lungs (Shi et al., 2004). Moreover, PstS3 is an excellent immunogen inducing CD8+ T-cell activation and both Th1 and Th17 immunity (Palma et al., 2011). Also, mice vaccinated with DNA coding for pstS3 demonstrated significant and sustained reduction in bacterial load in lungs after M. tb challenge (Tanghe et al., 1999). The observation of pstS3 induction during lung infection while downregulated during In-Vitro growth highlights the differential expression of virulence factors specific for the adaptation of M. bovis to host microenvironment.