Different Transcriptional Profiles of RAW264.7 Infected with Mycobacterium tuberculosis H37Rv and BCG Identified via Deep Sequencing

Background The Mycobacterium tuberculosis H37Rv and BCG effects on the host cell transcriptional profile consider a main research point. In the present study the transcriptome profiling analysis of RAW264.7 either infected with Mycobacterium tuberculosis H37Rv or BCG have been reported using Solexa/Illumina digital gene expression (DGE). Results The DGE analysis showed 1,917 different expressed genes between the BCG and H37Rv group. In addition, approximately 5% of the transcripts appeared to be predicted genes that have never been described before. KEGG Orthology (KO) annotations showed more than 71% of these transcripts are possibly involved in approximately 210 known metabolic or signaling pathways. The gene of the 28 pathways about pathogen recognition receptors and Mycobacterium tuberculosis interaction with macrophages were analyzed using the CLUSTER 3.0 available, the Tree View tool and Gene Orthology (GO). Some genes were randomly selected to confirm their altered expression levels by quantitative real-time PCR (qRT-PCR). Conclusion The present study used DGE from pathogen recognition receptors and Mycobacterium tuberculosis interaction with macrophages to understand the interplay between Mycobacterium tuberculosis and RAW264.7. Meanwhile find some important host protein which was affected by Mycobacterium tuberculosis to provide evidence for the further improvement of the present efficacy of existing Mycobacterium tuberculosis therapy and vaccine.


Introduction
Tuberculosis (TB) is one of the most common and widely distributed bacterial diseases where one third of the world population infected and causing 1.7 million deaths annually [1,2]. The causative agent of TB is Mycobacterium tuberculosis which can survive in host under a dormant state to delay its multiplication under adverse metabolic conditions for years [3]. TB continue to increase annually in spite of efficient chemotherapies which have been developed [4]. More seriously, strains of Mycobacterium tuberculosis which is resistant to a range of antimicrobials have emerged [5].
To face this problem, much research effort has been made in recent years towards understanding Mycobacterium tuberculosis and especially the interplay between Mycobacterium tuberculosis and host cells, expecting it will lead to the development of novel therapies and vaccine. DNA microarrays has been used to study the transcriptional profiling of Mycobacterium tuberculosis to investigate how Mycobacterium tuberculosis can destroy the host immune system response by Ward SK, 2010 [6]. However, the opportune response of host to Mycobacterium tuberculosis and the antigens they confront are still unclear. Simultaneously, the innate immune response to pathogenic bacteria and the subsequent adaptive immune response to eliminate the pathogens which can protect the host from pathogenic bacterial invasion have not been expound systematically [7].
As we known, H37Rv is a virulent Mycobacterium tuberculosis laboratory strain, where BCG is a widely used vaccine against Mycobacterium tuberculosis [8]. Thoroughly understanding the bacterial factors that affect some important pathways such as phagosomal antigen processing may help in improving strategies to eliminate Mycobacterium tuberculosis. H37Rv and BCG may have different effects on genes under adverse metabolic conditions positively selected in Mycobacterium tuberculosis. The functional discussion of these findings may imply that these selected genes may be involved in antigen variations and immune evasions of Mycobacterium tuberculosis. However, limited studies have reported on the large-scale the identification of immune-related genes at the genome or transcriptome levels in RAW264.7, which become a general cell model for TB studies nowadays, because current highthroughput deep sequencing technologies are insufficient [9,10]. The DGE system is a tag-based transcriptome sequencing approach in the production of technical innovation of RNA deep sequencing technologies; this system removes the limitation of study about immune-related genes and functional transcription complexes [11,12]. During the DGE approach, clean tags strained through raw tags with certain end nuclease should be mapped to the reference database. The degrees of expression of the actual total genes in the sample are the criterion defined by the quantity of each mRNA molecule yield from each gene in this approach [13]. Due to reasonable agreement on cell biology, DGE is mainly used to study the comparative gene expression [10,14]. The DGE approach can also be applied in the field of cellular, disease, and immune defense.
The present study is the first to conduct a transcriptome profiling to characterize the relationship between RAW264.7 gene expression profiles after Mycobacterium tuberculosis H37Rv and BCG infection using the DGE system. Furthermore, a round survey of the anti-bacterial immune defense gene activities in RAW264. 7 can contribute to the in-depth investigation of candidate genes in the host immunity defense against Mycobacterium tuberculosis. The results of the present study may help in improving the current understanding on host-pathogen interactions.

Materials and Methods
Cell and Bacterial Cultures RAW264.7 (ATCC TIB-71) was purchased from the National Institute for the Control of Pharmaceutical and Biological Products of china. RAW264.7 cells were maintained in supplemented Dulbecco's modified eagle medium (DMEM) composed of 10% fetal bovine serum, 0.2 mM L-glutamine, 100 mg streptomycin/mL, and 100 U penicillin/mL (Sigma, St. Louis, Mo.) at 37uC with 5% CO 2 [15].

Bacterial Cultures and Bacterial Growth Assay in vitro
H37Rv (ATCC27294) and BCG (ATCC19274-50) were purchased from the National Institute for the Control of Pharmaceutical and Biological Products of china. The H37Rv and BCG were grown as previously described using standard methods [16]. Bacterial growth assay in RAW264.7 was done as previously described [17]. showed that the distribution of tags matched the distribution of genes for each group. Furthermore, the frequencies of tags or genes decreased with increasing number of tags or gene expression. Results suggested that these DGE dates are significant. doi:10.1371/journal.pone.0051988.g001

Mycobacterium tuberculosis Infection
Bacteria were grown to mid-log phase before being pelleted and resuspended in RPMI 1640 medium plus 10% fetal bovine serum and HEPES. RAW264.7 cells were plated in 75 cm flasks in supplemented DMEM without fetal bovine serum and antibiotics 1 day prior to infection at a concentration of 1610 6 per flask. Macrophages (6610 6 cells per flask) were infected with different bacterial strains at a multiplicity of infection (MOI) of 5 for 4 h [18]. Following this, colony counts were done as previously described [18].

RNA Extraction
The total RNA was extracted from RAW 264.7 cells infected with Mycobacterium tuberculosis H37Rv and BCG using standard protocols (Trizol) [19]. To eliminate the potential genomic DNA contamination, RNase-free DNase I (Qiagen) was used according to the manufacturer's protocol. The RNA concentrations were evaluated using Agilent 2100 Bioanalyzer (Agilent Technologies). A 1.5% (w/v) agarose gel was used to analyze the integrity of RNA.

Library Construction and Sequencing
Illumina Gene Expression Sample Prep Kit and Solexa Sequencing Chip (flowcell) were used for the library construction and sequencing [20]. The main instruments used were the Illumina Cluster Station and Illumina HiSeq TM 2000 System. Library construction and sequencing were done as previously described [20].

Analysis and Mapping of DGE Tags
Firstly, dirty tags were filtered as previously described [21], to acquire clean tags containing CATG and 21 bp tag sequences for mapping the DGE tags. Subsequently, the clean tags were mapped to the Mus musculus reference database, allowing one nucleotide mismatch at most, for the study of tag annotation. The clean tags were labeled as clear-cut clean tags. The saturation analysis of each sequence was performed as described previously [20] to verify whether the number of detected genes continue to increase  with increasing sequencing amount (total tag number). For the gene expression analysis, the number of clear-cut clean tags for each gene was calculated and normalized to tags per million (TPM).

Evaluation of the DGE Libraries
To compare the gene expression differences, the tag frequency in each DGE library was statistically analyzed according to the method described by Audic and Claverie [22]. To determine which genes are significantly differentially expressed in the libraries, a P-value of ,0.005, a false discovery rate (FDR) of ,0.01, and an estimated absolute log2-fold change of .0.5 were selected as threshold.

Gene Ontology and Kyoto Encyclopedia of Genes and Genomes (KEGG) Functional Enrichment Analysis for Differentially Expressed Genes (DEG)
During gene ontology analysis, the gene sets with P,0.1 should be deemed significantly enriched. An FDR of ,0.1 must be considered for the KEGG enrichment analysis [21] [28].

Hierarchical Clustering and Tree View
The genes and their library expressions were analyzed using the CLUSTER 3.0 (http://bonsai.hgc.jp/˜mdehoon/software/ cluster/software.htm) available on the Internet and the Tree View tool (http://www.yiiframework.com/forum/index.php?/ topic/9180-the-tree-view/).

Analysis of DGE Libraries
The global gene expression profiles of RAW264.7 to Mycobacterium tuberculosis were measured using the Solexa/Illumina's DGE system. Two separate RAW264.7 libraries were sequenced using massively parallel sequencing on the Illumina platform after infection with H37Rv or BCG. The summary of these two libraries is shown in Table 1 and Figure S1, where the total clean sequence tags numbers of the two libraries was 5.6 million ( Figure 1) and 120074 of which were distinct clean tags sequences. The proportion of the amount of distinct tags to the total tags in H37Rv challenged was 5.26% and in BCG challenged libraries was 5.33%. In saturation analysis of the capacity of libraries, with more and more sequence tags, the amount of new distinct tags was decreasing, in the condition that the amount of total sequence tags was big enough. As the amount of the total sequence tags was 1 million, the new distinct tags were not detected any more ( Figure  S2 and Figure S6). Both results of the distribution of the clean tags in the two libraries showed that highly expressed genes had a broad distribution in the clean tags, while narrow in distinct clean tags (Table 1). By contrast, the fraction of genes with low expression in the clean tags is minimal, but was found high in distinct clean tags.

Analysis of Tag Mapping
The mapping fraction of these distinct tag sequences to Mus musculus UniGene reference sequences could be seen in Figure S3. As polymorphism may exist in specimen, tolerances were defined to permit one mismatch in each alignment. Based on the criteria, 56.63% and 64.74% of distinct clean tags were mapped to the UniGene virtual tag database; 51.51% and 58.34% of the distinct clean tags were mapped unambiguously to the UniGene; 4.66% and 3.80% of the distinct clean tags did not match to the UniGene virtual tag database in the H37Rv and BCG-treated library, respectively ( Table 1). The distribution of sample tag positions is shown in Figure S4. As essential qualities of sequencing tags, the sense and antisense strands of the transcripts can be distinguished using the Solexa sequencing. By comparison, the ratio of sense to antisense strand of the transcripts was approximately 1.99:1 in the BCG-treated libraries, and in H37mRv-treated libraries was 1.75:1 ( Figure S5). In spite of the high number of antisense mapping events detected, the sense strand is of first importance on the part of the transcriptional regulation in the Mycobacterium tuberculosis-induced immune response.

Gene Expression Variations between BCG-challenged and H37Rv-challenged Libraries
In order to search for gene expression variation between BCGchallenged and H37Rv-challenged libraries, the amount of clean tags for each gene was calculated, and the differently expressed genes between two samples were identified using a rigorous formula developed by Audic et al [22]. According to Figure 2, the number of significantly DEGs between two samples was 1,917, of which 1,135 genes and 782 genes were up-regulated and downregulated, respectively. Specific information about DEGs can be seen in Table S2.

Validation of DGE Data via qPCR
We randomly selected six DEGs: MyD88, Cd68, ULK1, TLR2, Arpc1b and TNF-a, to validate the correctness and confidence of DGE data. Data were presented as fold changes in gene expression normalized to GAPDH in Figure 3. Pearson's correlation coefficient (r) showed that both the DGE and qPCR data were highly correlated, where the DEGs had a high consistency and the r values ranging from 0.681 (TLR2) to 0.995 (MyD88) between the two methods ( Figure 2). The results show that the DGE results are dependable.

KEGG and DEGs with Pathway Annotation
Among the 1,917 DEGs, 1,371 can be mapped into 210 pathways (Table S3). Based on the KEGG function classification, these mapped DEGs can be grouped into seven kinds of pathways: global pathway, metabolism, genetic information processing, environmental information processing, cellular processes, organismal systems, and human diseases. The pathways distribution showed in Figure 4. As can be seen, the human diseases pathways were showed as the highest amount of genes measured as 603 genes, while organismal systems pathways were 562. Meanwhile, genetic information processing pathways were 385, metabolism pathways were 383, cellular processes were 361. For the environmental information processing pathways and global pathways, the number was 224 and 212, respectively.
As the pathways are too much to study it, so we from our point of view selected 28 pathways based on the potential role,  Figure 4: among the 1,917 DEGs, 1,371 can be mapped into 210 pathways which be classified functionally into seven categories: global pathways, metabolism pathways, genetic information processing pathways, environmental information processing pathways, cellular processes pathways, organismal systems pathways, human diseases pathways based on the KEGG function classification, where the human diseases pathways showed highest amount of genes measured as 603 genes, while Organismal systems pathways were 562, Genetic information processing pathways were 385, metabolism pathways were 383, cellular processes were 361, environmental information processing pathways were 224 and Global pathways were 212. Results suggested that after H37Rv or BCG infecting macrophages or possibly DCs, the global pathway, metabolism, genetic information processing, environmental information processing, cellular processes, and organismal systems all adopt some changes as shown above. doi:10.1371/journal.pone.0051988.g004 interaction with Mycobacterium tuberculosis and trying to make the selection cover all the stages of infection and immune response, then made further study on the selected pathways, where we grouped it into six categories based on KEGG, and illustrate its distribution as shown in Figure 5.

Hierarchical Clustering and Tree View
The genes involved in the above 28 pathways and their expressions in three libraries (BCG treated, H37Rv treated and non-treated group as a control group) were analyzed using CLUSTER 3.0 available on the Internet and Tree View tool. The up regulated genes showed value over zero, while the down regulated genes showed value range from 0 to 22. The most up regulated genes belongs to the following pathways: cell cycle and apoptosis; MAPK signaling pathway, TGF-beta signaling pathway; ribosome; microbial metabolism in diverse environment; T cell receptor signaling pathway, antigen processing and presentation molecules. Meanwhile, the most down regulated genes belongs to the following pathways: phagosome, endocytosis and autophagy; calcium signaling pathway; metabolic pathways; Fc epsilon RI signaling pathway (as shown in Figure 6, Table S3, and  Table S5).

Gene Ontology (GO) Classification
We used the WEGO database to analysis our DEGs output, we made analysis for the whole 1917 DEGs, the analysis showed that 1767, 1730, 1722 genes maybe annotated in GO component, GO function, and GO process based on sequence homologies, respectively.
The DEGs from the selected 28 pathways were analyzed by the WEGO analysis. About 468, 466, and 468 genes could be annotated in GO component, GO function, and GO process based on sequence homologies, respectively. In the three main categories (cellular component, molecular function, and biological process) of the GO classification, ''cell,'' ''cell part,'' ''binding,'' ''catalytic,'' ''metabolic process,'' and ''cellular process'' were dominant (Figure 7, Figure S7 and Table S6).

Pathogen Recognition Receptors
Toll-like receptors (TLR) and NOD-like receptors are main members in the pathogen recognition receptors family (PRRs). Compared with BCG, H37Rv can upregulate TLR1, TLR2, and TLR7 expression ( Table 2). TLR2 is important for host to recognize lipoproteins. The expression of some nuclear factor just as Nfkbid, Nfkbiz, and Nfrkb were also upregulated in the H37Rvtreated library (Table S2).

Mycobacterium tuberculosis Interaction with Macrophages
In order to survival in macrophage, Mycobacterium tuberculosis interacts with macrophage through several different receptors, including complement receptor (CR), FccR, and TLR [29]. For the expression of complement receptor (CR), Cr1l was downregulated, and for the expression of FccR, Fcgr1 can be upregulated in the H37Rv-treated library compared with the BCG-  (Table 2). MyD88 and Hmha1 involved in the TLR signaling pathway were downregulated in the H37Rv-treated library in comparison with the BCG-treated library ( Table 2).
Pro-inflammatory cytokine signaling is crucial for an effective immune response against Mycobacterium tuberculosis [30]. The expression of IL-6, TNF-a, IL-15, tumor necrosis factor (Tnfaip3, Figure 6. Hierarchical clustering and Tree View analysis. Figure 6A. The metabolism pathways analysis. Figure 6B. The genetic information processing pathways analysis data. Figure 6C. The Environmental information processing analysis. Figure 6D. The cellular processes analysis Figure 6E. The Organismal systems analysis. Number~t heexpressionofgeneinH37Rvgroup-theexpressionofgeneinBCGgroup theexpressionofgeneincontrolgroup(uninfected) The up regulated genes showed value over zero, while the down regulated genes showed value range from 0 to 22. Results showed that the most up regulated genes belongs to the following pathways: cellular process pathways, such as cell cycle and apoptosis; environmental information processing pathways, such as MAPK signaling pathway, TGF-beta signaling pathway; genetic information processing pathways such as ribosom; global pathways such as microbial metabolism in diverse environment; organismal systems pathways, such as T cell receptor signaling pathway, antigen processing and presentation molecules, while the most down regulated genes belongs to the following pathways: cellular process pathways, such as phagosome, endocytosis and autophagy; environmental information processing pathways, such as calcium signaling pathway; global pathway such as metabolic pathways; organismal systems pathways, such as Fc epsilon RI signaling pathway (as shown in Figure 6, Table S3, and Tnfaip6, and Tnfaip8), IFN-b1 as well as the cytokine receptor (TGF-bR1, Tnfrsf12a, Ripk1, Traf3, Traf6, and Traf7) were different between the H37Rv-treated and BCG-treated libraries ( Table 2).
The comparative analysis between H37Rv-treated and BCGtreated library revealed significant expression for nine genes that involved in the delivery inhibition of vacuolar H+-ATPases to the phagosomes (Table S2). All these nine genes were downregulated in H37Rv-treated library compared with BCG-treated library. Among these nine down-regulated genes, seven involved in the ATPase H+ transporting (Atp5a1, Atp5d, Atp5h, Atp6ap, Atp6ap2, Atp6v0e2, and Atp5j2) ( Table S2). The expression of actin-related protein expressing genes (Arpc1b and Arpc4) is also lower in the H37Rv library compared with that in the BCG library (Table S2). The disruption of actin can inhibit the phagosomal maturation, which is the best-characterized mechanism for Mycobacterium tuberculosis to manipulate the macrophage for survival and create a favorable environment for replication [31].
Compared with those in BCG-treated library, Cep110 and Cep55, which may contain early endosomal antigen 1 (EEA1), were down regulated in H37Rv-treated library (Table S2). In addition the expression of Rab5b is up regulated and Rab5c is down regulated in H37Rv challenged library (Table S2). Both Rab5 and EEA1 are necessary for phagosomal maturation to proceed [32,33].
The expression of Med8 and Hmhal were down-regulated in H37Rv-treated library in comparison with BCG-treated library (Table S2), while the expression of Rab7 and Cd68 were up regulated and the expression of Scarb2 and Cd63 was down regulated in H37Rv-treated library in comparison with the BCG one (Table S2).
The expression of serine/threonine kinase, such as Ripk1, Pisd, Pxk, and Psat1, are different between the two libraries (Table S2). Among them, only Psat1 was highly downregulated, and all others were up regulated in the H37Rv-treatd library. Lysosomal hydrolases and serine/threonine kinase play an important role in regulating phagosomal maturation. Lysosomal hydrolases, such as cathepsins, glycosidases, sulfatase, lipase, and ceramidase were differently expressed between the two libraries. In 2007 Rohde K et al. demonstrated that PknG can regulate phagosomal maturation [34].
Autophagy regulation can eliminate Mycobacterium tuberculosis that resides within the phagosomes [16,35]. Results showed that UlK1, Med8 and Hmhal were down regulated, while Atg7 and Atg12 were up regulated in H37Rv-treated group in compare with BCG one.

Putative Novel Immune/stress Response Genes
The homology analysis showed that only 1828 transcripts had high sequence homology with the known sequences in the public database, while 89 transcripts didn't show homology and couldn't annotated clearly (Table S4), which maybe recommended it to be a putative novel immune related genes in RAW264.7 and involved in the Mycobacterium tuberculosis responses, which need further investigations and researches to make it more clear. Discussion A timeless Solexa/Illumina's DGE system was used to fully understand the whole transcriptional changes in the BCG-and H37Rv-challenged RAW264.7 with qPCR done to ensure the accuracy of the data. The transcriptome can explain all expressed RNA transcripts in a cell [9,10]]. Furthermore, some cellular activities in organisms, such as disease and immune defense, can be well studied [9,10].
There were 1,917 DEGs between H37Rv-challenged and BCGchallenged libraries (Figure 2, Table S2), around 1,135 genes up regulated while 782 genes down regulated, which suggested that these DEGs might be putative novel immune-relevant genes in RAW 264.7 involved in the response to Mycobacterium tuberculosis H37Rv and BCG challenges. Among the 1,917 DEGs, 1,371 mapped into 210 pathways (Table S3), grouped into seven kinds of pathways ( Figure 4). Results suggested that after H37Rv or BCG infecting macrophages or possibly DCs, the global pathway, metabolism, genetic information processing, environmental information processing, cellular processes, and organismal systems all adopt some changes as shown above.
TLRs (TLR1-TLR13) are one of the signaling families of PRRs which can identify pathogen-associated molecular patterns (PAMPs) [36]. DGE analysis showed that the expression of TLR1, TLR2, and TLR7 were up-regulated at different levels in H37Rv-challenged library compared with that in BCG-challenged library ( Figure 6E). TLR2 can interact with TLR1 or TLR6 to form heterodimers which can engage in the recognition of mycobacterial cell wall glycolipids, such as triacylated lipoproteins (TLR2/TLR1), or diacylated lipoproteins (TLR2/TLR6) [37][38][39]. This suggests that TLR-based immunity may participates in RAW264.7 defense against Mycobacterium tuberculosis.
In 2004 Gutierrez MG et al. demonstrated that autophagy was a defense mechanism inhibiting BCG and Mycobacterium tuberculosis survival in infected macrophages [16]. In 2007 Xu Y et al. demonstrated that the MyD88-independent gene was shown to be involved in lipopolysaccharide-induced autophagy [40]. Our DGE analysis showed the MyD88 transcript was downregulated in the H37Rv-challenged library compared with that in BCG-challenged library. We make a forecast that as a consequence of their co-evolution, pathogens mimic host activities and  TPM represented tags per million. H37Rv (QD) represented the group in which RAW264.7 were treated with H37Rv; BCG (RD) represented the group in which RAW264.7 were treated with BCG. TPM represented tags per million. Symbol, Gene symbol (Brief gene description can be found in Table S2). Results showed that, the expression of some DEGs which involved in recognition and ingestion, complement receptor, cytokine and cytokine receptor, Endosome, Lysosome, lysosomal membrane protein, Lysosomal acid hydrolase and Fc gamma Rmediated phagocytosis, suggesting that these genes may be the involved in the main host protein which was affected by H37Rv or BCG. doi:10.1371/journal.pone.0051988.t002 regulation mechanisms to influence the autophagy by downregulation of the MyD88 transcript (unpublished data). The pro-inflammatory cytokine signaling is important for an accurate and effective immune response against Mycobacterium tuberculosis. DGE analysis demonstrated that the expression of IL-6, TNF-a and IL-15 were up-regulated in the H37Rv-treated library compared with those in BCG-treated library. If we want to improve present tuberculosis therapy, we may pay more attention to the function of these cytokine to effectively control Mycobacterium tuberculosis infection.
The well-known mechanism for Mycobacterium tuberculosis to evade the immunologic surveillance and the elimination by macrophage is to inhibit phagosome maturation. The delivery inhibition of vacuolar H+-ATPases to the phagosome may be the main reason for the delay of phagosome maturation [34]. The results showed that the expression of seven H+ transporting genes (Atp5a1, Atp5d, Atp5h, Atp6ap, Atp6ap2, Atp6v0e2, and Atp5j2) were all downregulated in H37Rv-challenged library. Thoroughly understanding the function of these important host genes may have important role in increasing the efficient of existing Mycobacterium tuberculosis vaccine (unpublished data).
During the interaction of phagosome and endosome, actin is essential for the scaffolding of endosomes [41][42][43]. DGE analysis demonstrated that the expression of actin-related protein expressing genes (Arpc1b and Arpc4) was also lower in H37Rvchallenged library. This suggests that Arpc1b and Arpc4 may be part of the reason for Mycobacterium tuberculosis-induced disruption of actin and the delay of phagosomal maturation [44,45].
EEA1 and Rho-GTPase Rab5 are also essential for phagosome maturation [46,47]. The present results in H37Rv group showed that the expression of Cep110 and Cep55, which may be involved with EEA1, and Rab5c were down regulated, while Rab5b was up regulated in comparison with BCG one. We suspect that Cep110 and Cep55 may be part of the tethering factors, in concert with Rab proteins, specifically bridging membranes to be fused with Rab5 supporting early and Rab7 supporting late fusion events. Meanwhile we implicated that Rab5b and Rab5c have antagonistic roles on the trafficking of endosomal and reconciling the fusion of phagosome with other organelles (unpublished data).
Our results showed that there was down regulation for the UlK1 expression in H37Rv treated group in comparison with BCG ones (Table S2). Regulation of autophagy pathway analysis showed Ulk1 involved in the ATG 1 expression which has an important role in autophagy induction [48] (Table S3). This suggests that U1K1 may be part of the reason for Mycobacterium tuberculosis-induced the inhibition of autophagy induction.
The expression of Med8 and Hmhal were down-regulated in H37Rv-treated library in comparison with BCG-treated library (Table S2). We predicted according to regulation of autophagy pathway analysis (Table S3) that this may result PI3K down regulation which will inhibit the phagosomal maturation by influencing many fusion and fission events [49].
The Atg7 and Atg12 have important roles in docking and fusion, showed up regulation events in H37Rv-challenged library compared with that in BCG-challenged library (Table S2) [50]. This suggests that H37Rv do not make any further inhibion of the Vesicle expansion and completion in compared with BCG.
This study mainly used DGE to understand the interplay between Mycobacterium tuberculosis and RAW264.7. In a word, after H37Rv or BCG infecting macrophages or possibly DCs, the global pathway, metabolism, genetic information processing, environmental information processing, cellular processes, and organismal systems all adopt some changes as shown above. As Mycobacterium tuberculosis has evolved an array of mechanisms to allow survival within the macrophage, we just want to find some important host protein which was affected by Mycobacterium tuberculosis to provide evidence for the further improvement of the efficacy of existing Mycobacterium tuberculosis therapy and vaccine. Figure S1 Distribution of total clean tags and distinct clean tags in experimental and control library. The numbers in square brackets indicate the range of copy numbers of each tag category. The data in parentheses indicate the percentage of corresponding tags among the total clean tags and distinct clean tags. The results confirmed that the quality of the data were up to specification suggesting that DGE data are significant. (TIF) Figure S2 Saturation of DGE libraries. The saturation analysis of capacity of libraries showed that new emerging distinct tags were gradually decreased with increasing total sequence tags when the number of sequencing tags was high enough. (A) H37Rvtreated library; (B) BCG-treated library; (C) control library. In saturation analysis of the capacity of libraries, with more and more sequence tags, the amount of new distinct tags was decreasing, in the condition that the amount of total sequence tags was big enough. As the amount of the total sequence tags was 1 million, the new distinct tags were not detected any more. Results suggested that the saturation of DGE libraries fulfill the need of DGE's requirement. (TIF) Figure S3 Mapping of total and distinct clean tags. The mapping fraction of these distinct tag sequences to Mus musculus UniGene reference sequences could be seen in Figure S3. (TIF) Figure S4 The positions of tags in the gene. A) the H37Rvinfected library; (B) the BCG-infected library; (C) the control library. Ideally, the tag is the 3 most one. But for alternative splicing or incomplete enzyme digestion, the tag may be the 2 nd or 3 rd from the 3 most. (TIF) Figure S5 Sense vs antisense transcripts expression in the experimental and control library. (A) H37Rv-infected library; (B) BCG-infected library; (C) control library. As essential qualities of sequencing tags, the sense and antisense strands of the transcripts can be distinguished using the Solexa sequencing. By comparison, the ratio of sense to antisense strand of the transcripts was approximately 1.99:1 in the BCG-treated libraries, and in H37Rvtreated libraries was 1.75:1 as shown in Figure S5. Results suggested that in spite of the high number of antisense mapping events detected, the sense strand is of first importance on the part of the transcriptional regulation in the Mycobacterium tuberculosisinduced immune response. (TIF) Figure S6 Effect of library size on the number of genes identified. (A) H37Rv-infected library; (B) BCG-infected library; (C) contral library. The increasing tendency in the rate of increase in all genes identified and the genes identified by unambigous tags declined with increasing library size. The rate of increase of all genes identified and genes identified by unambigous clean tags declined drastically as the size of the library increased. When the library size reached one million, we could identify 35% and 30% all genes and genes identified by unambigous clean tags, in (A) H37Rv-infected library and (B) BCG-infected library, respectively. Simultaneously, we could identify 48% and 30% all genes and genes identified by unambigous clean tags, in (C) control library. At this time, library capacity approached saturation. Results suggested that the saturation of DGE libraries fulfill the need of DGE's requirement. (TIF) Figure S7 Signaling pathways of the DEGs. The pathway analysis was mainly based on the KEGG database. During gene ontology analysis, the gene sets with P,0.1 should be deemed significantly enriched. An FDR of ,0.1 must be considered for the KEGG enrichment analysis. The vertical axis is the pathway category, and the horizontal axis is the log10 (p value) of these pathways. We from our point of view selected 28 pathways, which may have some potential function for in-depth studies on Mycobacterium tuberculosis. Results suggested that in these 28 selected pathways the DGEs mainly distributed in global pathways, cellular processes pathways and organismal Systems pathways. (TIF)

Table S2
The whole information about DEGs in the H37Rvinfected library and BCG-infected library. Each column stands for: the number in ''gene'' represented Gene ID; Raw Intensity, Raw clean tag number; TPM (Transcripts Per Million clean tags) is a standardized indicator, pointing out number of transcript copies in every 1 million clean tags. This is Sum of the normalized clean tag number in Table S2 Table S4. These 28 pathways were selected from our point of view which may have some potential function for in-depth studies on Mycobacterium tuberculosis. Results suggested that in these 28 selected pathways the DGEs mainly distributed in global pathways, cellular processes pathways and organismal Systems pathways. In addition, more attention may be given to these Mycobacterium tuberculosis-affected immune-related genes in RAW264.7. (XLS )   Table S5 KEGG categorization of potential pathway based on KEGG. Results showed the number of DEGs in every pathway which belongs to the KEGG categorization among these 28 pathways selected from our point of view which may have some potential function for the in-depth studies on Mycobacterium tuberculosis.

(XLS)
Table S6 Complete information on the GO annotation on DEGs belonging to the potential pathway. What each column stands for just as shown in Table S2. The DEGs from the selected 28 pathways were analyzed by the WEGO analysis. About 468, 466, and 468 genes could be annotated in GO component, GO function, and GO process based on sequence homologies, respectively. In the three main categories (cellular component, molecular function, and biological process) of the GO classification, ''cell,'' ''cell part,'' ''binding,'' ''catalytic,'' ''metabolic process,'' and ''cellular process'' were dominant. (XLS)