Phenotypic and Transcriptomic Response of Auxotrophic Mycobacterium avium Subsp. paratuberculosis leuD Mutant under Environmental Stress

Mycobacterium avium subsp. paratuberculosis (MAP) is the causative agent of severe gastroenteritis in cattle. To gain a better understanding of MAP virulence, we investigated the role of leuD gene in MAP metabolism and stress response. For this, we have constructed an auxotrophic strain of MAP by deleting the leuD gene using allelic exchange. The wildtype and mutant strains were then compared for metabolic phenotypic changes using Biolog phenotype microarrays. The responses of both strains to physiologically relevant stress conditions were assessed using DNA microarrays. Transcriptomic data was then analyzed in the context of cellular metabolic pathways and gene networks. Our results showed that deletion of leuD gene has a global effect on both MAP phenotypic and transcriptome response. At the metabolic level, the mutant strain lost the ability to utilize most of the carbon, nitrogen, sulphur, phosphorus and nutrient supplements as energy source. At the transcriptome level, more than 100 genes were differentially expressed in each of the stress condition tested. Systems level network analysis revealed that the differentially expressed genes were distributed throughout the gene network, thus explaining the global impact of leuD deletion in metabolic phenotype. Further, we find that leuD deletion impacted metabolic pathways associated with fatty acids. We verified this by experimentally estimating the total fatty acid content of both mutant and wildtype. The mutant strain had 30% less fatty acid content when compared to wildtype, thus supporting the results from transcriptional and computational analyses. Our results therefore reveal the intricate connection between the metabolism and virulence in MAP.


Introduction
Mycobacterium avium subsp. paratuberculosis (MAP) is a multispecies pathogen and is the etiological agent of a severe gastroenteritis in ruminants, known as Johne's disease [1]. The clinical symptoms of MAP infection in cattle include diarrhea, weight loss, decreased milk production, and mortality. Clinically or subclinically infected animals shed MAP in feces and/or milk enabling easy dissemination to other animals [1,2]. MAP is also a concern in human health as it has been suspected as a possible cause of Crohn's disease in humans [3,4].
Considering its obvious importance in animal industry and possible importance in human health, different strategies including use of DNA and recombinant antigen vaccines have been used in attempts to reduce the economic loss against MAP infection [5,6,7,8]. Production of live attenuated strains of MAP with vaccine potential using transposon mutagenesis and allelic exchange has also been attempted [9,10,11]. Majority of these studies have utilized the advances in genetics research in related Mycobacteria, particularly Mycobacterium tuberculosis.
Several virulence and antigenic genes have been characterized in M. tuberculosis. The leuD gene which encodes isopropyl malate isomerase, an essential enzyme for leucine biosynthesis, is an example. In mice, the leuD mutant of M. tuberculosis has been demonstrated to be protective against challenge with virulent M. tuberculosis [12]. Similarly the leucine auxotroph of M. bovis, created by deletion of leuD using allelic exchange, was found to be protective in cattle against M. bovis challenge [13]. leuD knock-outs have been also used as auxotrophic selectable markers in M. avium and M. smegmatis [14,15]. In another study focused on M. bovis BCG strain, leuD was found to be involved in an oxidative stress response [16]. leuD has been identified as part of PhoPR system in M. tuberculosis [17] and as a component of codY regulon in Listeria monocytogenes [18]. Recently we created a leuD knockout (MAP-DleuD) strain using allelic exchange method [19]. We also demonstrated that MAPDleuD produced partial protection against MAP challenge [19]. Although these studies have clearly established the potential of knockout strain of leuD as a vaccine candidate, the role of leuD in virulence associated functions is largely unexplored. Since leuD is part of central intermediary metabolism, it will be particularly useful to understand how these pathways are involved in virulence associated functions. In this study, in order to elucidate the mechanism of leuD attenuation in MAP, the mutant and wildtype were characterized using phenotype array (PM) and the response of MAPDleuD was assessed by transcriptional profiling under different physiologically relevant stress conditions. We find that leuD metabolism is intricately connected with virulence associated genes and also fatty acid biosynthesis. Although we have detailed knowledge on molecular function of microbial virulence factors, our understanding of how nutrition of bacteria during infection controls these virulence genes is at its beginning. Therefore, the results of this study adds significantly to our understanding of how leuD biosynthetic pathway is connected to virulence functions.

Bacterial strains
As the genome of M. paratuberculosis strain K10 has been sequenced, we used this strain to develop a mutant of leuD gene by allelic exchange, using methods described previously [19]. MAP Wild type (MAP-WT) and MAPDleuD generated was grown in 7H9 medium supplemented with 10% oleic acid-albumindextrose-catalase(OADC) (Becton Dickinson Co., Sparks, MD) and 1.25 mg/L mycobactin J (Allied Monitor, Inc., Fayette, MO).

Phenotype microarrays
Phenotype microarray experiments were performed following standard Biolog Inc. (http://www.biolog.com/) protocols [20,21]. Briefly, MAP strains were grown on Middlebrook 7H10 agar at 30uC with supplementation of casitone and were resuspended in Biolog inoculating fluid so as to have 81% transmittance. Aliquots (100 ml) of this suspension were used to inoculate Phenotype microarray panels. Full metabolic profiles covering 760 tests were performed. Metabolite utilization at 37uC was measured every six hours up to 138 hours as color changes of redox dye G using an ELISA reader at 630 nm. Tests were performed in duplicate, and the mean signal intensity was calculated for each replicate. Each test well was then compared against negative control well of respective PM plate using ANOVA. Signal intensity increase of 25% or more from negative control at 5% significance level was considered positive phenotype.

Stress experiments
Bacteria were grown to a cell density of 0.7 at OD 600 in 7H9 broth with supplementation of casitone (7H9-casitone) and were subjected to the following three stressors for 3 hours (i) shift from 30uC to 42uC, (ii) shift from pH 7.0 to pH 5.5 and 9.0 and (iii) anaerobic incubation. Growth in medium without casitone was used as the fourth stress condition (minimal medium). Since leuD mutant could not grow without the supplementation of casitone in 7H9 broth, it was grown in 7H9-casitone broth as described above and was shifted to 7H9 without addition of casitone and was incubated for 2 days. Three biological replicates were conducted for each stress condition. At the end of stress experiments, twice the volume of RNA protect reagent was added to bacterial cultures; these were incubated at room temperature for 5 minutes. Ten milliliters of cells were then pelleted by centrifugation at 13,000 g for 5 minutes and were stored in 280uC.
RNA extraction, cDNA synthesis, template labeling, array hybridization and data analysis Cell pellets were transferred to a 2 ml sterile screw-cap microcentrifuge tube. 1 ml of TRIzol reagent (Invitrogen, Life Technologies) and 0.6 ml of sterile RNase-free 0.1 mm zirconia beads were added to each tube. Cells were homogenized and lysed by bead beating four times in a Mini Bead-Beater (BioSpec Products, Inc., Bartlesville, OK) for 30 seconds with a gap of 30 seconds. Cell lysate was transferred to a 1.5 ml sterile RNasefree microcentrifuge tube. RNA from the lysate was cleaned up using RNeasy MinElute Cleanup Kit (Qiagen, Valencia, CA). Before elution of RNA from the MinElute columns, following manufacturers protocols, possible contaminating DNA was removed using Qiagen RNase-Free DNase Set. Integrity of isolated RNA was estimated using Agilent Bioanalyzer 2100. Only those samples with an RNA Integrity Number (RIN).8 was used for cDNA synthesis. 3.0 mg of RNA from each sample were reverse-transcribed using SuperScript III reverse transcriptase with random hexamers (Invitrogen) according to the manufacturer's instructions. The corresponding cDNA was labeled with 25 U of Klenow fragment (New England Biolabs, Ipswich, MA), 7 mg of exo-resistant random primers (Fermentas, Glen Burnie, MD), dNTP mix (0.12 mM each dATP, dTTP, and dGTP and 0.03 mM dCTP), and 0.1 mM of Cy3/Cy5-dCTP (Amersham Biosciences, Piscataway, NJ) at 37uC for 12 hours. Labeling reaction was cleaned up using a QIAquick PCR Purification Kit (Qiagen). A custom expression array for MAP K10 genome was designed using Agilent E-array algorithm (https://earray.chem. agilent.com/earray/). This array consisted of 60mers representing each CDS in MAP K10 genome [Gene Expression Omnibus (GEO) platform ID# GPL15269]. Following standard Agilent protocols, labeled cDNA was hybridized against this custom array. Cy3 labeled cDNA from wildtype strain served as control channel and Cy5 labeled leuD mutant cDNA from the corresponding stress condition served as the test channel. All three biological replicates were hybridized in duplicates. Arrays were then scanned using Agilent G2565BA microarray scanner and data was extracted using Agilent's feature extraction image analysis software v 10.2. The complete set of microarray files has been deposited in GEO (accession # GSE36053). Resulting data files were then loaded to GeneSpring GX v.11 (Agilent Technologies, CA). To allow interpretation of the data, values below 0.01 were set to 0.01. The log2 ratio of the intensity of signal channel to that of control channel was calculated and was then normalized to each slide's median intensity ratio. The default interpretation mode was set to 'log of ratio'. Each time point was then tested separately for statistical differences using ANOVA at 5% confidence level, in conjunction with Benjamini-Hochberg multiple testing correction. From this list, those genes changing more than 1.5 fold and 2.0 fold were selected as differentially expressed genes (DEGs).
Quantitative real-time PCR (qRT-PCR) 1.0 mg of the total RNA was converted to cDNA by using SuperScript III reverse transcriptase (Invitrogen) with random hexamers, following the manufacturer's instructions. The real-time reaction mixture included 5 ml of cDNA template, 200 nM of each of both forward and reverse primers, and 16 Power SYBR Green master mix (Applied Biosystems, Foster, CA). Primers used in this study are listed in the Table S1. RT-PCR was performed in 96well optical plates using the ABI 7500 Sequence Detection System instrument and software (Applied Biosystems). The thermal cycling consisted of an initial denaturing step at 95uC for 10 min followed by 40 cycles consisting of 95uC for 15 s and 60uC for 1 min. For C t determination, three parallels were assayed for each gene and MAP_4078 gene was used as the internal reference gene. qRT-PCR was performed in triplicate to evaluate reproducibility of data.

Computational analysis
The interaction data for MAP was extracted from STRING database [22]. Since online version of STRING database supports only a maximum of 1,000 nodes, complete interaction data from STRING v 9.0 was downloaded as a full database dump. MAP specific interaction data was then extracted using custom Linux shell scripts. The interaction data in STRING database includes both experimental as well as predicted interaction information [22]. Interactions in STRING are provided with a confidence score. The interactions with a confidence score less than 700 were removed and this interaction network was then imported into Cytoscape [23] for visualization and overlaying of microarray expression data. Effect of leuD deletion in MAP pathways was analyzed using Pathway tools software [24,25]. An Omics viewer [25] was used to visualize the differentially expressed genes on MAP pathways. Enrichment analysis implemented in Pathway tools was used to indentify statistically significant pathways associated with each stress condition. The mycolic acid biosynthesis pathway in MAP was deduced using M. tuberculosis pathway as reference. Takayama et al. has defined the M. tuberculosis genes associated with mycolate biosynthesis pathway [26]. This gene list was used as input for ortholog mapping tool in Integrated Microbial Genomes [27] to find the MAP orthologs. Result of this mapping is given in Table S2. For analyzing the effect of leuD deletion on the virulence, MAP virulence associated gene definitions were downloaded from the Virulence factor database (VFDB) [28] and were combined with microarray expression data.

Analysis of fatty acids and fatty acid methyl esters (FAME)
Three 100-ml biological replicates of MAP-WT and MAPDleuD were grown to mid-exponential-phase in 7H9 medium. MAP-DleuD was then shifted to 7H9 medium without casitone and incubated for 12 hours. MAP-WT was incubated in fresh 7H9 medium for the same duration. The respective replicates of both MAP-WT and MAPDleuD were combined, and the frozen pellets were submitted for FAME analysis at Microbial ID, Newark, DE (http://www.microbialid.com/). The FAME data was analyzed using SherlockH V 6.1.

Growth characteristics of MAPDleuD
Previously, we reported the development of MAPDleuD strain and its potential as a vaccine candidate [19]. MAPDleuD is an auxotroph of leuD which requires supplementation of 7H9 with Casitone in order to culture in that medium. Here, casitone was used as a leucine source since our strain had an impaired leucine biosynthetic pathway and cannot grow in un-supplemented 7H9 (Fig. 1A). With the addition of casitone, growth was restored albeit to a slightly less degree than the wild type strain (Fig. 1A). To confirm that the impairment of the growth was only due to the disruption of leuD but not due to the so called ''polar effect'' on the genes upstream or downstream, we compared the mutant and wildtype strains using qRT-PCR on leuD locus and adjacent genes. The results show that while the expression of leuD is completely abolished in MAPDleuD, expression of all other genes upstream and downstream are unaffected (Fig. 1B).

Deletion of leuD leads to an impaired metabolic phenotype
Biolog phenotype microarray (PM) for microbial cells is a highthroughput assay for testing utilization of various nutrients. The assay chemistry uses a redox dye to measure respiration of cells [20]. When the bacterial cell is able to utilize a specific test nutrient in a given PM well, redox dye is reduced producing purple color. Since the dye reduction is mostly irreversible under physiological conditions, the color accumulates in the well over a period of hours, amplifying the signal and integrating the amount of respiration over time [20]. There are eight PM plates which test for the utilization of the following nutrient sources; PM1, 2carbon sources, PM3-nitrogen sources, PM4-phosphorus and sulphur sources, PM5-nutrient supplements, PM6, 7 and 8peptide nitrogen sources. In this study, we have utilized PMs 1-8 to establish the complete metabolic profile of MAP WT and MAPDleuD. Of the total of 760 phenotypes tested, 182 tests were positive for MAP WT (Table S3).PM results show that the most preferred carbon source for MAP was Tween 80. Most preferred simple sugars were D-ribose, D-tagatose, D-fucose, D-trehalose and a-D-glucose respectively. Tween 40 and Tween 20 were utilized at the rate similar to D-ribose. Among amino acids, Lleucine was utilized most efficiently as carbon source. The most used nitrogen sources were alloxan and D-mannosamine. The overall results for positive phenotypes in PM plates 1-5 are given in Fig. 2. As shown in Fig. 2, when compared to MAP WT, the ability to utilize most of the nutrients was impaired in MAPDleuD. A similar trend was observed for peptide nitrogen sources (Fig. 3). Therefore, the deletion of leuD gene resulted in the impairment of overall nutrient utilization capacity of MAP. However, one 1A. Both MAP-WT and leuD mutant were grown in 7H9 medium and 7H9 medium supplemented with casitone as source of leucine. Growth was measured as increase in turbidity of culture medium at OD600. 1B. qRT-PCR measurement of transcriptional levels of genes upstream and downstream of leuD. qRT-PCR products for each gene was loaded into a 1.2% agarose gel and was visualized using a gel documentation system. doi:10.1371/journal.pone.0037884.g001 exception was the ability of MAPDleuD to utilize phosphonoacetic acid while MAP WT did not use this chemical effectively. leuD deletion has a large impact on MAP transcriptional response In order to understand the role of leuD in MAP stress response, using microarrays we compared the transcriptional profile of mutant and wildtype after subjecting these to different stressors such as shift in pH, nutrient, temperature and anaerobic growth. Similar conditions have been used previously to define the ''stressome'' of MAP [29]. Anaerobic growth and pH shift were used to mimic the exposure of MAP inside macrophages and in the host abomasum. Growth in the absence of casitone (minimal medium) was used to define the impact on metabolic pathways and temperature shift was used to define general stress response. Two days of growth in the absence of casitone was selected based on the growth pattern of mutant as shown in Fig. 1A. We used three hours of exposure for all other stress conditions as it has been demonstrated to be adequate in a previous study [29]. After transcriptional profiling, two levels of cutoff for defining differentially expressed genes (DEGs) at 5% confidence level were considered. The first set was DEGs with at least 1.5 fold change when compared to control (Fig. 4A) and second set was DEGs with at least 2.0 fold change (Fig. 4B). At both cutoff levels, a large number of DEGs were present in all stress conditions. There were 26 DEGs present in all conditions tested at 1.5 fold cutoff. However, this number diminished to one when the cutoff was raised to two fold levels (Fig. 4B). A comparison of the total number of DEGs at these two cutoff levels is given in Fig. 4C. An overview of expression pattern of 26 DEGs that were found changing in all stress conditions at 1.5 fold cutoff is given in Fig. 4D. The details about the DEGs at both of these cutoff levels are given in Table S4 and S5.

Deletion of leuD affects transcriptional response of known virulence genes in MAP
Comparative phylogenomics is a very useful method for deducing the functions of less studied bacterial genomes by extrapolating data from related species. Although there are not many functional genomics studies in MAP, many related species such as M. tuberculosis, M. smegamatis, M. bovis, M. avium, and M. leprae have been studied in more detail. We used the curated phylogenomics data on other mycobacteria from Virulence Factor Database (VFDB) [28] to identify virulence associated genes in MAP (Table S6). We then combined this virulence gene list with microarray expression data from different stress conditions (Fig. 5). As Fig. 5 shows, leuD deletion influenced the expression profile of most of the virulence associated genes. It is interesting to note that leuD gene itself is annotated as a virulence associated gene in VFDB. Some of the virulence genes in VFDB have also been studied extensively in MAP. For example, the functions of MAP_0216 (fbpA), MAP_0187c (sodA) and MAP_3921(sodC) which were found to be modulated in this study are well defined.

DEGs are dispersed throughout the MAP gene network
Since the Biolog phenotype microarrays and functional microarray data revealed a global effect on MAPDleuD, we used genome wide systems approach to analyze the DEGs. The interaction data for MAP was derived from STRING [22]. This database contains known and predicted protein-protein interactions which include physical or functional associations [22]. Evidence for these is derived from genomic context, highthroughput experiments, coexpression and literature mining [22]. After combining different types of evidences for each interaction, only interactions with score of 700 or more was considered. This score cut off represented medium and high confidence interactions for MAP. The final interaction network contained 3,350 nodes (genes) and 28,932 edges (interactions). When 2.0 fold DEGs were visualized on this network, we find that they are scattered throughout the interaction network (Fig. 6). If 1.5 fold DEGs are considered, this effect is more pervasive (data not shown). This reveals that leucine metabolism is tightly integrated with the whole cellular networks and would explain the severe impairment of MAPDleuD's ability to utilize most of the nutritional sources in Biolog phenotype arrays.
We further attempted to elucidate the role of 26 DEGs that were found changing at least 1.5 fold in all stress conditions. We used manual search of STRING database to find the interaction partners for these genes (Fig. 7A). This analysis predicted the following interaction partners with very high confidence evidence scores; groL, MAP_3239, MAP_3058c, dnaK and MAP_2204c (Fig. 7B). Cumulative results of this analysis suggest that 26 DEGs found to be changing in all stress conditions are involved in heat shock response and maintenance of cell viability.
leuD deletion affects several MAP cellular pathways, particularly the fatty acid biosynthesis pathways Since leuD gene is part of super pathway of leucine, valine and isoleucine biosynthesis, it is natural to reason that deletion of this gene will have major impact on other cellular pathways of MAP. For analyzing the DEGs in the context of cellular pathways, we used the Biocyc [30] pathway reconstruction of MAP. Pathway tools software [25] was used to visualize 2.0 fold changing DEGs on the cellular pathways (Fig. S1). As expected, fatty acid biosynthesis II, palmitate biosynthesis II, mycolate biosynthesis, fatty acid beta-oxidation, acetate degradation and oleate betaoxidation pathways were found to be differentially expressed. Other differentially expressed pathways included pathway of assimilatory sulfate reduction, malonate degradation, acetyl-CoA fermentation to butyrate I, glutaryl-CoA degration, arginine biosynthesis II, and pantothenate and coenzyme A biosynthesis I. We used the pathway enrichment tool in Pathway tools to identify statistically significant pathways changes associated with each stress condition. The result of this analysis is given Table 1. As evident in Table 1, there were a total of 53 pathways that were associated with different stress conditions. The following main compounds are either produced or consumed by leucine biosynthesis pathway; (2R, 3S)-3-isopropylmaleate, (2S)-2-isopropyl-3-oxosuccinate, (2S)-2-isopropylmaleate, 2-oxoisovalearate, 4methyl-2-oxopentanoate, and L-leucine. In addition, the following side compounds in leucine biosynthesis pathway forms links with other pathways; 2-oxoglutarate, acetyl-CoA, coenzyme A, NAD+, NADH and L-glutamate. As shown in Fig.S1, the main effect of leuD deletion is on the pathways dependent on or linked to altered levels of needed primary or side compounds produced or consumed by the leucine biosynthetic pathway. For example, we found that key mycolate biosynthetic pathway genes are modu- lated in different stress conditions, particularly during growth in minimal medium. Our results measuring mRNA levels mirrored changes in gene activity observed in microarray tests by performing qRT-PCR on the same set of genes in this particular pathway (Fig. 8) showing that both microarray and qRT-PCR test results are in agreement with each other.

MAPDleuD has reduced cellular lipid profile
Since the transcriptional analysis suggested that the fatty acid synthesis pathways in MAPDleuD are affected, we verified this possibility using FAME analysis of MAP-WT and MAPDleuD. SherlockH FAME analysis system uses gas chromatography and sophisticated pattern recognition software to identify over 300 fatty acids in the range of 9 to 20 carbons in length. The total amount of fatty acids detected by FAME analysis of MAP-WT and MAPDleuD were 396831 and 275561 units, respectively. This corresponds to 30% reduction of fatty acid content in mutant cells in equivalent amount of biomass. The most abundant fatty acids were 18:1 w9c, 16:00 and 16:1 w9c (Table 2). Some peaks that were present in MAP-WT were missing in MAPDleuD. The FAME spectrum analysis result for MAP-WT and MAPDleuD is given in Fig. S2. Overall, the FAME analysis results are in agreement with computational analysis of transcriptional data suggesting that deletion of leuD leads to alteration of cellular lipid profile of MAP.

leuD is intricately linked to global MAP metabolism
In this study, for the first time, we report the complete nutrient utilization profile of MAP using phenotype arrays. One of the surprising results we obtained is that Tween 80 was used by MAP as carbon source and was even more preferred than most common simple sugars such as D-ribose and a-D-glucose (Table S3). Other Tween types such as Tween 40 and 20 also were utilized well as carbon sources. Generally mycobacterial culture media are  supplemented with small amounts of Tween as a surfactant to avoid cell clumping and more quantities are believed to be toxic. Addition of Tween 80 in the medium has been reported to change the susceptibility of M. avium-M. intracellulare complex [31]. On the other hand, Tween 80 has been reported to be hydrolyzed by M. bovis [32] and M. avium subsp. paratuberculosis [33]. Further, M. smegmatis grown in the presence of Tween 80 synthesizes series-2 type glycolipids [34]. It has also been reported that bacteria such as Microthrix parvicella uses the oleic acid moiety of Tween 80 as carbon and energy source [35]. Therefore, it is reasonable to assume that MAP also uses oleic acid component of Tween 80 in the same manner. PM results showed that while L-leucine is used as carbon source by MAP-WT, the MAPDleuD lost this ability. This indicates that external supplementation of leucine is not as efficient as native biosynthesis of this amino acid and could explain why supplementation of casitone did not restore growth rates of MAPDleuD to those of wildtype MAP. Among pentose sugars such as D-trehalose that was shown to be used as carbon source have well recognized role in mycobacteria. For example, trehalose is an essential component of mycolic acid [26]. However, we could not find annotated enzymes for the metabolism of sugars such as Llyxose, D-tagatose, and D-fucose. Rather than the absence of transferases and isomerases required for the utilization of these sugars, it most likely indicates the need for annotating these genes correctly in MAP genome. Reason for this suggestion is that these enzymes have been identified in related mycobacteria such as M. smegmatis [36] and M. avium [37]. Since there is very limited literature regarding nutrient sources such as alloxan, xanthine and peptide nitrogen sources, we could not link the metabolism of these compounds and the enzymes and pathways in MAP. The overall observation that MAPDleuD lost the ability to use most of the nutritional sources reveals the intricate connection between MAP metabolism and virulence. The only instance of reversal of nutrient utilization was the ability of MAPDleuD to utilize phosphonoacetic acid. Although it has been reported that some mycobacteria utilizes phosphonoacetate as sole carbon, phosphorus and energy source [38], further studies are required to understand the importance of this ability reversal. In addition to revealing the intricate link between leuD and global metabolism of MAP, the PM data we report here can have several other uses, including (i) improvement of bacterial metabolic and systems models [39,40,41,42], (ii) assessing the accuracy of genome annotation [43], and (iii) developing new improved culture media for bacteria [44].

Effect of leuD deletion on MAP transcriptional response to stress
In addition to studying the effect of leuD deletion at the phenotype level, we also compared the transcriptional response of MAP-WT and MAPDleuD. In the past there have been a very limited number of studies examining the MAP transcriptional responses to deletion mutants [29,45,46,47,48]. However, this is the first report where MAP transcriptional response is analyzed after deletion of a metabolically important MAP gene. As expected, we found that a large number of genes were modulated under different stress conditions. The 26 genes that were found to be modulated in all stress conditions in our study can be candidates of further characterization. Although the fold change Figure 6. Visualization of Differentially Expressed Genes (DEGs) on MAP gene interaction network. Interaction information of MAP genes was retrieved from STRING database. Interactions with confidence score more than 700 (medium and high confidence interactions) were considered. DEGs with 2.0 fold or more change were then visualized on this network using Cytoscape. DEGs that were changed at different stress conditions are indicated using different color code given at top right side of the figure. Nodes (genes) are indicated by locus tags. Nodes in grey color were not differentially expressed. doi:10.1371/journal.pone.0037884.g006 cutoff for these genes is 1.5 fold, the fact that these genes are differentially expressed in all conditions indicates their importance in MAP environmental responses. The majority of these was hypothetical proteins and could be candidates for further characterization. While groES (MAP_4264) and groEL (MAP_4265) were up-regulated during temperature shift, they were found to be down-regulated several fold in all other conditions. The heat shock proteins were also found to be among the DEGs in the previous MAP ''stressome'' study [29]. Two other annotated DE genes common to all conditions were pfpl (MAP_0372) and thil (MAP_1699). It has been reported that pfpl plays an anti-mutator role and provides general stress protection Pseudomonas aeruginosa [49]. In Salmonella enterica serovar Typhimurium it has been shown that ThiL has dual role in de novo biosynthesis and in salvage of exogenous thiamin [50]. We used an integrated approach of using systems level gene interaction networks and phylogenomics to understand the link between leuD deletion and metabolic and transcriptional responses. As Fig. 6 shows, the DEGs were distributed throughout the MAP gene network. This would explain why leuD deletion had such huge effect on the global phenotype. The dispersed nature of DEGs in Figure 7. Gene interaction network of 26 genes that were found to be changing at all stress conditions. 7A. STRING database was searched using list of 26 genes that were found to be differentially expressed 1.5 fold or more at all stress conditions tested. Nodes (genes) are indicated by locus tag or gene name. Edges (interactions) are denoted by the color code for source of evidence and are indicated at the top side of Fig. 7B in the inset. 7B. List of predicted interaction partners for these 26 DEGs. List of predicted interaction partners are given as locus tag or gene name in the left most column. The evidence source for each interaction is represented as columns. Cumulative evidence score is given in the right most columns. doi:10.1371/journal.pone.0037884.g007 Figure 8. Comparison of gene expression fold change in key genes associated with mycolate biosynthesis pathway. Genes associated with mycolate biosynthesis were identified using pathway annotation in Biocyc and by ortholog mapping of genes as defined in mycolate biosynthesis pathway of M. tuberculosis. Gene expression fold change was measured using both microarrays, and qRT-PCR. doi:10.1371/journal.pone.0037884.g008  the MAP gene network also explains why a large number of virulence associated genes were also modulated in the stress conditions tested.

Effect of leuD deletion on MAP pathways
Since leuD deletion leads to the abolition of leucine biosynthesis, it is not surprising to see that a total of 53 pathways significantly modulated in the stress conditions studied. Previously using a mouse model, we have demonstrated that vaccination of MAPDleuD strain provided significant protective efficacy against a MAP virulent strain challenge [19]. The pathway analysis results suggested that this protective efficacy could be due the effect of leuD deletion on metabolic pathways associated with fatty acids, particularly mycolate biosynthesis (Fig. S1). Our results suggest that in MAPDleuD, reduced activity of leucine biosynthesis pathway leads to diminished levels of 2-oxoglutarate, acetyl-CoA, coenzyme A, NAD+, and NADH. These compounds are required to complete several steps in synthesis of mycolic acid and other fatty acids. We verified this possibility experimentally by estimating the levels of fatty acids in MAP WT and MAPDleuD. Although the overall spectrum (Fig. S2) and fatty acid peaks appear similar in both strains, FAME analysis results shows that in MAPDleuD there is about 30% reduction in fatty acid levels ( Table 2). The reason for similarity in peaks in the FAME spectrum is that initially MAP-WT and MAPDleuD were grown to mid-exponential-phase in 7H9 medium and MAPDleuD was switched to minimal medium. This initial culturing of both strains in the same medium seems to be resulting in similar spectrum appearance. However, the 30% reduction in overall fatty acid levels does support our pathway analysis results that leuD deletion has a major effect on pathways associated with fatty acids. It is well known that the membrane lipid profile is important in mycobacterial survival and virulence [26,51,52,53]. In addition to this, leuD deletion also appears to be influencing mycolate biosynthesis pathway be modulating genes such as fbpA. Mycolic acids are synthesized as two alkyl chains by enzymes similar to those that assemble fatty acids and these chains are fused and modified numerous times to generate a complex family of lipids that are incorporated into the cell envelope [54,55,56]. Secreted proteins FbpA, FbpB, and FbpC attach mycolic acids to arabinogalactan, generating mycolic acid methyl esters (MAME) or trehalose, generating a,a9-trehalose dimycolate (TDM) and is essential for mycobacterial virulence [56]. Therefore, the mechanism behind the attenuation of MAPDleuD strain could be due to reduced levels of cellular fatty acids and mycolic acids. Since mycolic acid biosynthesis pathway is generally conserved in Mycobacteria, it is natural to assume that the link between leuD metabolism and mycolic acid biosynthesis we indentified in MAP could be also be present in other Mycobacteria such as M. tuberculosis. Further, several of the MAP phenotypes we defined here, could also be conserved in other Mycobacteria.

Supporting Information
Table S1 Details about primers used for qRT-PCR. (XLSX)   Figure S1 Visualization of 2.0 fold changing DEGs on MAP pathways. The DEGs are highlighted on pathways according to the following color codes; (i) Blue-pathways not associated with DEGs, (ii) red-pathways mapped to pH 5.5 associated DEGs, (iii) orange-pathways mapped to pH 9.0 associated DEGs, (iv) greenpathways mapped to temperature associated DEGs, (v) pinkpathways mapped to anaerobic associated DEGs, (vi) brownpathways mapped to minimal medium associated DEGs. Pink lines indicate pathways consuming metabolites from leucine pathway. (TIF) Figure S2 FAME spectrum of MAP-WT and MAPDleuD. Raw gas chromatography spectra for both mutant and wildtype strains are shown. Peaks in the GC spectrum were identified using SherlockH V 6.1.