Human alveolar macrophage response to Mycobacterium tuberculosis: immune characteristics underlying large inter-individual variability

Background: Mycobacterium tuberculosis (M.tb), the causative bacterium of tuberculosis (TB), establishes residence and grows in human alveolar macrophages (AMs). Inter-individual variation in M.tb-human AM interactions can indicate TB risk and the efficacy of therapies and vaccines; however, we currently lack an understanding of the gene and protein expression programs that dictate this variation in the lungs. Results: Herein, we systematically analyze interactions of a virulent M.tb strain H37Rv with freshly isolated human AMs from 28 healthy adult donors, measuring host RNA expression and secreted candidate proteins associated with TB pathogenesis over 72h. A large set of genes possessing highly variable inter-individual expression levels are differentially expressed in response to M.tb infection. Eigengene modules link M.tb growth rate with host transcriptional and protein profiles at 24 and 72h. Systems analysis of differential RNA and protein expression identifies a robust network with IL1B, STAT1, and IDO1 as hub genes associated with M.tb growth. RNA time profiles document stimulation towards an M1-type macrophage gene expression followed by emergence of an M2-type profile. Finally, we replicate these results in a cohort from a TB-endemic region, finding a substantial portion of significant differentially expressed genes overlapping between studies. Conclusions: We observe large inter-individual differences in bacterial uptake and growth, with tenfold variation in M.tb load by 72h.The fine-scale resolution of this work enables the identification of genes and gene networks associated with early M.tb growth dynamics in defined donor clusters, an important step in developing potential biological indicators of individual susceptibility to M.tb infection and response to therapies.


Background
Mycobacterium tuberculosis (M.tb) infects approximately one-quarter of the world's population, and nearly 1.3 million people die annually from tuberculosis (TB), a continuing worldwide public health problem (1). Since only a fraction of infected individuals progress to active TB, we must understand the immunopathogenesis of M.tb infection, including host susceptibility or resistance factors, to develop effective diagnostic, therapeutic and vaccine strategies tailored to different individuals and populations (2).
Aerosolized M.tb enters the alveolar lung space where it is phagocytosed by alveolar macrophages (AMs), unique resident cells with a complex immunologic pro le, to become an intracellular pathogen. M.tb infection activates several macrophage immunobiological pathways involved in phagocytosis, vesicle tra cking, and triggering of in ammatory cytokines, oxidants, and cell death pathways -all processes of innate immunity. Yet, critical factors in the human host that dictate M.tb infection and progression to TB remain uncertain -a roadblock to understanding an individual's susceptibility to TB (3)(4)(5).

Heterogeneity in M.tb uptake, adaptation and growth in human AMs among donors
We isolated human AMs from 28 adult healthy donors (Table 1) and soon after (within 6h) infected them with a virulent M.tb strain (H 37 R v ) engineered to emit light (expressed in relative luminescence unit or RLU values) over 2h, and followed RLUs over 72h (Fig. 1A,B). Given the focus of this study on host variation in response to infection, we selected a single, well-studied strain (H 37 R v ). The assays presented throughout were conducted at a multiplicity of infection (MOI) of 2:1 (M.tb/AM cells), and macrophage monolayers were maintained throughout the assays. RLUs correlate with M.tb colony forming units (CFUs), jointly re ecting the number of intracellular bacilli and microbial metabolic activity (29). RLU levels at 2-24h post-infection largely re ect M.tb cellular uptake and early intracellular adaptation. After correcting for batch effects, large variations in RLU values were detectable at each time point (Additional le 1: Table  S1). Since large RLU differences at 2h post infection could have arisen from variations in early adaptation of M.tb to the intracellular environment, they were not further considered. At 24h, early events largely stabilized, and RLU increases from 48 to 72h accurately re ected M.tb growth rates. AMs with the highest and lowest RLUs displayed a nearly tenfold difference at 72h. RLU ratios (48/24h and 72/48h) were similar using an MOI of 10:1 (Additional le 1: Table S1). During exponential growth, the generation time is the time interval over which a bacterial population doubles in number, and represents the growth rate. We estimated the generation time for M.tb in each donor from RLU values at 48-72h post infection, nding these to be broadly consistent with prior estimates of ~ 15h while varying two-fold between donors. To identify proteins and genes that in uence M.tb-AM interactions, we focused on generation times (48-72h) in this study, since RLU levels prior to 48h post infection re ect M.tb uptake, adaptation, and early growth combined and, in addition, since some values at the early time points fell below the detection limit of the instrument.  Table S2). At 2h, 6 of the 27 proteins (22.2%) were signi cantly differentially produced after M.tb infection (Fig. 1C), increasing to 23 proteins (85.2%) at 72h (Fig. 1E), all with increased secretion in infected cells. Differential secretion varied over time, with IL6 and MMP2 displaying the highest levels at 2h, and IL6, CCL13, and CCL3 at 72h. Pronounced M.tb-stimulated secretion of MMP2 only at 2h (without stimulated mRNA expression; not shown) suggests a mechanism of rapid release from the cell into the extracellular matrix rather than de novo transcription. In ammatory mediators including IL1β, TNF, and IL6 displayed signi cant differential secretion across the time points, indicating their contribution to a pro-in ammatory reaction to infection.
We then tested whether protein levels secreted from infected cells correlated with M.tb generation times (Fig. 1F-H; Additional le 3: Table S3). Three proteins were signi cantly correlated with M.tb generation times (FDR adj. p values < 0.05): IL7 (2h, 72h), IL10 (2h, 24h), and MMP10 (24h), with Pearson correlation reaching R 2 = 0.5. However, IL7 levels (typically produced by lymphoid cells) were at the lower range of the assay sensitivity and may not be reliable. High IL10 levels signi cantly correlated with longer generation times (slower growth) (r = 0.46 and 0.49 at 2h and 24h, respectively). MMP10 levels were signi cantly correlated with slower growth at 24h, though switch to become positively (though nonsigni cantly) correlated with rapid growth at 72h. These results support the notion that protein expression pro les over the rst 3 days of infection in uence M.tb-AM interaction dynamics.
Transcriptomes of uninfected control and M.tb -infected human AMs reveal differentially expressed gene pro les Transcript pro les were assessed for each AM sample at 2, 24, and 72h post infection using AmpliSeq (30). This method has been shown to be scalable across a large concentration range (30). At each time point 10,000-14,000 mRNAs were detectable, and reads per million (RPM) from replicate assays were highly correlated (r 2 ≥ 0.99), enabling sensitive detection of differentially expressed (DE) genes. A principal component analysis (PCA) of all datasets from control and infected AMs revealed that 52% of the variance in gene expression between time points resulted from exposure to ex vivo culture conditions alone, while M.tb-induced RNA expression changes were smaller but increased over time (Additional le 4: Supplementary Fig. 1). We previously showed that ex vivo cultured HAMs rapidly (starting at approximately 6h) begin to display the transcriptional pro le of monocyte-derived macrophages (MDM) suggesting a de-differentiation process is driving the changes we see (30). The dataset we explore here supports this (Additional le 5: Supplementary Fig. 2). Canonical human AM markers are highly expressed at 2h in both control and infected samples, but downregulated over 72h. In contrast, MDM markers are upregulated during ex vivo culture.
Although there is a strong impact of ex vivo culture on gene expression pro les at 24 and 72h, infected samples are clearly distinguished from control samples by PCA of gene expression and by differential expression of secreted proteins. Macrophage polarization is often simpli ed to the trajectory towards M1like or M2-like states, the former representing a pro-in ammatory state. To explore the relevance of these polarized states to differences between control and infected samples we rst tested whether the phenotype of the cultured macrophages changed over time by deconvolution of the AmpliSeq pro les with CIBERSORT-X. Across both control and infected samples we saw a shift to an M0 macrophage population -likely re ecting the proposed shift of cells towards an MDM-like state. We identi ed a higher proportion of M1-like macrophages in infected samples than control samples. Conversely, in control AMs a shift toward M2 cell type was observed by 72h ( Fig. 2; Additional le 5: Supplementary Fig. 2).
To identify signi cant DE genes (FDR adjusted p ≤ 0.05) speci c to infection, we analyzed RNA expression by comparing uninfected control cells to infected cells at each time point (30). This yielded 62 DE genes at 2h, 2,177 at 24h, and 3,662 at 72h post infection (Additional le 6: Table S4). In a previous study also using AmpliSeq (30) on macrophages from three donors, we had detected 16 DE genes at 2h, 899 at 24h, and 174 at 72h, showing substantial overlap with DE genes in this larger study. Thirty-six DE genes were signi cant at all-time points ( Fig. 3A and Fig. 4), including those encoding in ammatory cytokines (IL1A, IL1B, TNF) and chemokine receptors (e.g., CCR7) characteristic of a spectrum representing "M1 type" cell states.
Our cell type deconvolution identi ed changes in the abundance of M1-like and M2-like macrophages as a major difference between control and infected cells. As deconvolution is contingent on the reference data sets we sought to further investigate the relationship to M1/M2-like macrophages and differential expression. Using an independent assignment of M1-and M2-associated genes ( and M2-like gene signatures (32). To assess the states of the macrophages responding to infection we performed gene set enrichment analysis of the 72h DE gene pro les, using gene co-expression modules previously identi ed from the full spectrum of macrophage differentiation (33). This showed a signi cant enrichment of M1-like modules (Fig. 3F), supporting our initial observation, but also enrichment of gene modules from other macrophage states.
As one example, induction of the highly signi cant DE transcript STAT1 (a transcription factor associated with an M1-like state) varied between AM donors from a 2-fold reduction to a 25-fold increase. These results indicate that classically de ned M1 and M2 cell types play a major role in the early macrophage responses. Enrichment in gene co-expression modules correlated with a spectrum of macrophage phenotypes [(e.g., modules correlated with TNF, PGE2 and Pam3CysSK4 or P3C (TPP)] stimulation which correspond to chronic in ammation (Fig. 3F). These results are consistent with AMs representing a unique macrophage phenotype that does not t neatly into a classically de ned M1 or M2 cell type but exists along a spectrum of macrophage types (Additional le 5: Supplementary Fig. 2) (32, 34-36).
Differentially expressed genes (DE genes) displaying the most variable expression in AMs between donors (VE genes) The DE genes we identi ed capture differences in the transcriptional response between infected and uninfected control cells, although this did not identify genes that show high variability between individuals -one important goal of this study. To capture inter-individual variation in gene expression we performed variance analysis of gene expression in both control and infected AMs at all time-points, yielding a set of 324 genes with highly variable expression (VE genes; Levene's test, ratios of variances < 0.15; Additional le 7: Table S5A). All 324 VE genes are also signi cant DE genes which were found to be upregulated in and speci c to the infected AMs . Examples include IFI6, IL1B, CCL4, IDO1, GBP5, IRF1,  JAK3, UBD, CXCL5, CCL20, VDR, CD80, IFI44L, NLRP3, and IL7R, several of which were previously implicated in TB pathogenesis (37). Correlations between expression of these genes with M.tb generation times are listed in Additional le 7: Table S5A.
Among the VE genes, expression of IDO1 (a well characterized marker of M.tb infection) was substantially stimulated by M.tb in only 10 of the 28 donor AMs (Fig. 5A). IL1B was more broadly expressed, but with high expression mostly coinciding with high IDO1 expression (Fig. 5A), similar to several other key genes (e.g., IFNG, UBD, GBP5, IFI44L, CXCL9-11; Fig. 5B-H). These most variably expressed VE genes, including STAT1, are likely relevant to inter-individual variability following infection of AMs. A degree of co-expression of genes with proposed opposing effects on M.tb, e.g., IL1B and IFNG (better control) versus IDO1 and Type 1 IFNs (reduced control) -suggests counteracting in uence on growth during the early phase of infection and requires further experimentation.

Reactome analysis yielded several signi cant pathways, including Immunological Diseases/Cellular
Function and Maintenance/In ammatory Response, and Function/Immune Cell Tra cking (Additional le 7: Table S5B). Ingenuity Pathway Analysis (IPA) of the VE genes yielded the highest scoring hierarchical network with IL1B at the top, and with STAT1, IRF1, and IDO1 as connected nodes (Fig. 6, full gene annotation in Additional le 8: Supplementary Fig. 3). Both STAT1 and IRF1 interact with the IDO1 promoter (38) and play a role in macrophage polarization.
Identi cation of human AM gene networks associated with M.tb growth Single RNA pro les across the donor AM dataset displayed high correlations with M.tb growth dynamics but failed to yield signi cance upon FDR corrections. For future studies, it is likely that increased sample size will be valuable in the identi cation of correlations between growth and single gene expression. Several candidate genes strongly correlated between reads per millions (RPMs) at 24h and generation time, suggesting biologically relevant relationships. As naïve FDR approaches do not account for the correlation between the expression of groups of genes a common approach is to collapse the expression of genes which show highly similar patterns into gene expression modules, with each module randomly assigned a color for identi cation. Eigengenes representing the expression level of genes within a module for an individual can then be correlated with a phenotype. To this end we built three consensus networks using batch-corrected RNA pro les from infected and control cells at 2, 24 and 72h, with RNAs having RPMs > 35 using WGCNA. This partitioned 9,193, 9,098 and 9,115 genes into 26, 27, and 37 modules for the 2, 24 and 72h time points, respectively (Additional le 9: Table S6A). Correlation between network eigengene modules for infected cells and M.tb generation times (during 48-72h) at 2, 24 and 72h identi ed four signi cant associations between bacterial growth rate and eigengene modules, all in the 24h network (FDR-corrected p-values < 0.05; blue r=-0.52, yellow r=-0.61, red r=-0.55 and tan r = 0.52; Additional le 10: Table S7, columns D,E, modules highlighted in red). Reactome pathway over-enrichment analysis (Additional le 11: Table S8) showed Processing and Metabolism of RNA, and DNA Damage/Repair as signi cantly associated with these gene modules. Regulation of genes involved in DNA repair and recombination had been linked to M.tb regrowth from its non-replicating, persistent state (39). This result supports a nexus between gene networks and M.tb generation times, supporting the notion that M.tb-stimulated AM gene expression patterns modulate bacterial growth that vary between individuals.
To prioritize individual genes within signi cantly associated modules, we measured correlations between generation time and gene expression for each gene within the four signi cant modules (Additional le 9: Table S6B). We lowered the stringency of our approach by performing multiple test correction based upon the module size, rather than for all expressed genes. This approach yielded 143 signi cant genes, with 139 genes belonging to the yellow module in the 24h network (97.2%). The yellow module genes had mostly negative correlations with generation times (117 of 139 genes, 84.2%); i.e., elevated expression of these genes was associated with shorter generation times (higher growth rates), perhaps signaling a key M.tb-induced AM response favoring bacterial growth.
We performed hierarchical clustering of gene expression pro les and AM donors for genes within the yellow module from the 24h network, limiting to genes where we found a signi cant correlation between generation time and expression level (Fig. 7A). This revealed four groups of genes, two of which contained genes with a positive correlation to generation time (clusters I and II), and two with a negative correlation to generation time (clusters III and IV). Hierarchical clustering also identi ed three groups of AM donors, with groups 1 and 2 displaying signi cantly longer generation times compared to group 3 ( Fig. 7B). The correlation coe cients (r 2 values) were also determined from 24h gene expression data to assess their distribution across all modules (Fig. 7C).
Gene expression sub-groups displayed opposite pro les between groups 1 and 3, with group 2 in between. We performed KEGG pathway enrichment of gene expression clusters I and II (positively correlated with generation time) and clusters III and IV (negatively correlated with generation time). This identi ed the TNF signaling pathway as signi cantly enriched in negatively correlated genes (Fig. 7D), and oxidative phosphorylation and thermogenesis signi cantly enriched in positively correlated genes (Fig. 7E), both were unexpected results. The co-expressed gene group representing a TNF signaling pathway with rapid growth does not include TNF itself, and is likely a consequence of rapid M.tb growth. The slower M.tb growth rate in group 1 is likely re ective of less M.tb-mediated inhibition of oxidative phosphorylation. Only some genes involved in oxidative phosphorylation genes fell into group 3, suggesting overlapping regulatory pathways with differential effects on M.tb growth. Nevertheless, our results further support a biological relationship between gene expression pro les and M.tb generation times.
Eigengene modules and the secreted protein response to M.tb infection in human AMs Among the measured proteins implicated in TB pathogenesis, IL10 and MMP10 were signi cantly correlated with bacterial growth (Fig. 1). Although in our in vitro system the other measured proteins did not show an association with M.tb growth, they are all known to impact TB disease progression through interactions with the immune system and the regulation of their expression is of major interest (40). To identify the transcriptional modules that may regulate or respond to these proteins, we assessed correlations between eigengene modules and secreted protein levels (Additional le 12: Table S9), identifying 29 signi cant correlations after Bonferroni correction. We opted for Bonferroni correction at this point for greater stringency in capturing correlations. Six of twenty-nine signi cant correlations were in control AMs, 3 of which were at the 2h time point where there was little divergence in gene expression between control and infected conditions. Of these, CCL3 expression was strongly correlated with the royal blue module [Additional le 11: Table S8 (column b) and Additional le 12: Table S9] and may be driven by ex vivo adaptation of macrophages.
In the 24h network, we found three modules with a signi cant link to M.tb growth rate and to both host transcriptional and protein pro les (Additional le 9: Table S6B and Additional le 10: Table S7). These are CCL3 (blue, red and yellow modules) and IL10 (red module). Signi cant correlations between these gene modules and individual proteins also suggest a causative relationship between gene expression and M.tb growth. At 72h, a single gene expression module (tan) correlated with the expression of multiple proteins (IL15, IL18, IL1A and TNF). This module overlapped substantially with our highest scoring network from IPA analysis of variably expressed VE genes (Additional le 8: Supplementary Fig. 3 Table S10 and Additional le 14: Supplementary Fig. 4. At 2h post infection, fewer RNAs were assigned DE status (45), and only 3 upregulated RNAs overlapped with the DE genes displayed above, likely owing to the small number of samples and the less precise and less targeted RNA-Seq methodology. However, at 24h, 360 of 773 DE genes (46.6%), and at 72h, 89 of 148 DE genes (60.1%) overlapped with the DE genes reported in the main study, although the donors of two different demographic cohorts were genetically very diverse and of different ages. Pathway analyses revealed a similar spectrum of enriched functional terms for upregulated genes, including cytokine signaling in immune system, interferon signaling, and innate immune system. Inspection of key DE genes, irrespective of AM donors and infection time points, replicates the nding of substantial inter-individual differences, for example IDO1 and its co-expressed genes (Additional le 15: Supplementary Fig. 5).

Discussion
The protein and RNA expression pro les parallel large inter-individual variation in interactions between healthy human AMs and a single virulent strain of M.tb (Figs. 1 and 2, measured with RLUs). For the rst time this study dissects biological mechanisms underlying individual differences. Sousa et al. (9) had demonstrated that secretion of IL1B is a surrogate marker distinguishing between TB cases with mild disease and severe TB, attributing differences in IL1B induction to virulence of the strains tested, independent of the host. However, here we show that a single virulent M.tb strain elicits substantial differences in IL1B expression between individual AMs, highlighting host factors, with IL1B one possible predictive marker of TB susceptibility. A smaller replication study identi es an overlapping set of DE genes and con rms substantial inter-individual variability of critical DE genes such as IDO1 (Additional le 15: Supplementary Fig. 5).
The nexus between M.tb growth rates and secreted proteins and eigengene modules further supports the hypothesis that gene expression networks, both in uninfected and M.tb-infected cells, affect M.tb-AM interactions. While we have applied stringent FDR adjustments to reveal signi cant correlations, strong candidate genes and networks can serve to extract a deeper understanding of the various phases of M.tb-AM interactions and their possible signi cance in the pathogenesis of TB.
We rst measured 27 candidate proteins that are secreted by M.tb-infected AMs and uninfected controls, nding robust stimulation by M.tb and inter-subject differences (Fig. 1C-H). Importantly, two proteins (IL10 and MMP10) are signi cantly correlated with bacterial growth over 48-72h (Additional le 3: Table   S3), suggesting the hypothesis that proteins relevant to TB pathogenesis affect early M.tb-AM interactions, possibly presaging individual susceptibility to TB. Production of both IL10 and MMP10 have been reported to promote the disease caused by another bacterial pathogen, e.g., Helicobacter pylori (41,42).
RNA expression pro les also appear to re ect or in uence M.tb-AM interactions. The precision of the AmpliSeq RNA assay, together with use of control uninfected AM expression at each time point (30), enabled sensitive detection of numerous DE RNAs across the 28 donors for AMs ( Fig. 3; Additional le 5: Table S4). DE genes signi cant at each time point (Fig. 4) are consistent with previously reported results [e.g., (22,28,30)], including in ammatory cytokines CSF2, IL6, IL1B, and IFNG, but also anti-in ammatory factors such as CCL22. Gene Ontogeny (GO) Pathway analysis reveals the expected prevalence for Innate Immune System, Cytokine Signaling in Immune System, and more.
Macrophages were initially characterized as polar phenotypes, i.e., M1 or M2 (43). However, AMs represent a unique spectrum of phenotypes and express M1 and M2 markers (Additional le 5: Supplementary Fig. 2) (35). In this study, we used published modules (31) that have de ned classical M1 or M2 type genes as the starting point. We found that the DE genes include both a classically activated M1 type induced by in ammatory agents (e.g., IFNG) and alternatively activated M2 type induced by IL4 and IL13 (20). The M1 type is strongly enriched in infected samples, while M2 markers increase at 72h particularly in the uninfected controls ( Fig. 2; Additional le 5: Supplementary Fig. 2). Substantial differences between M1 and M2 markers occur between individuals, reinforcing variability of an individual's response to infection. The balance of this mixed mode of activation may account in part for inter-subject differences in M.tb-AM interactions, but no signi cant trends were observed with respect to growth rates. In previous reports, M.tb ESAT-6-induced macrophage polarization to M1 phenotype occurred early, then switched to M2 phenotype at a later stage of infection (44). In our study, a switch to the M2 phenotype at a later stage (72h) was also observed in uninfected control macrophages compared to infected AMs, which could be due to increased expression of MMPs in control cultures. Upregulation of certain MMPs was found to be associated with macrophage polarization to both M1 and M2 phenotypes in M.tb-infected, cigarette smoke-exposed macrophages (45). The majority of individuals in our primary cohort were non-smokers. Future work on the impact of smoking and MMPs can be best assessed in the African replication cohort which were predominantly smokers.
Among the gene co-expression modules, one (yellow) signi cantly correlates with M.tb generation times and, in addition, contains numerous individual genes that also signi cantly correlate with M.tb growth (Additional le 9: Table S6). As most correlations are negative, these genes appear to favor M.tb growth. The module GO terms indicate in uence over translation processes, a potential mechanism by which M.tb usurps the cellular machinery to its advantage (46). VE genes can reveal markers of individual susceptibility or resistance to disease phenotypes (47) Supplementary Fig. 3). This IL1B-dominated network connects IL1B, STAT1, and IRF1 to IDO1 (Fig. 6), consistent with previous reports that the IRF1/STAT1 transcription complex binds to the IDO1 promoter (38, 48).
We observed a signi cant increase in GBP2 expression in M.tb-infected macrophages versus controls at 24 and 72h (highly signi cant at 72h). GBP2 is one of the genes in the blood based RISK6 transcriptome proposed as a biomarker for TB disease and treatment response (49). However, GBP2 expression did not correlate with M.tb growth generation times (Additional le 7: Table S5A). Expression of the other 5 genes in the RISK6 signature was not signi cantly different in M.tb-infected AMs, suggesting some speci city for blood vs tissue responses. IDO1 has been proposed recently as a sensitive and selective biomarker for active TB (26). Its metabolic products, immunosuppressive kynurenins, act by preventing access of cytotoxic T cells to infected macrophages in TB lung granulomas (27). IDO1 is robustly expressed in only 10 of 28 AMs after M.tb stimulation, with a correlated expression pro le for IL1B, CXCL9-11, IFNG, UBD, IFI44L and GBP5 (an IFN/IL1B activated GTPase mediating antibacterial defense) (Fig. 5). We propose that M.tb-induced expression of both IL1B and IDO1 could be early predictive indicators of susceptibility to M.tb for a given individual if similar results are obtained in more readily accessible cells (e.g., MDMs). Indeed, increased macrophage IL1B expression was observed in non-human primate granulomas that are more permissive to M.tb growth (50). However, increased IL1B could simply be increased secondary to bacterial burden. Increased IL1B activity is typically associated with M.tb control (51).
The nexus we observe between secreted proteins reported relevant to TB and gene expression modules, with overlap to variably expressed VE/DE genes and, on the other hand, to M.tb generation times, strengthens our hypothesis that early M.tb-AM interactions could signal an individual's susceptibility to infection and potentially risk of overt TB.

Conclusions
In summary, our results identify known and novel key genes and their encoded transcripts and proteins in the early phase of M.tb infection of freshly isolated human AMs. Previous work shows that infection of macrophages with M.tb stimulates multiple pathways, including both interferon type I and II pathwaysthought to exert opposing effects on infections (20,21). We show here that variably and differentially expressed genes and their networks affecting M.tb-AM interactions are associated with inter-individual variability during the early phase of infection. Results of this study provide a foundation for identifying predictive biological indicators of an individual's susceptibility or resistance to M.tb infection, and response to therapies and vaccines (52).

Limitations
The current study employs freshly isolated AMs (used within 6h of harvest) and hence is distinct from MDMs isolated from peripheral blood after 5 days of in vitro incubation, accounting for substantial differences of RNA pro les between them (30). However, during AM incubation for 72h, we observed substantial transcriptome changes even in the absence of infection, with a tendency towards the pro le of MDM transcriptomes as we have found previously (30). Nevertheless, the cultured AMs maintain distinct characteristics from MDMs over the incubation period. We have mitigated any in vitro artifacts with use of AM controls at each time point for each donor. In addition, very large inter-individual differences occur at 2h post infection, likely caused by variable cellular adaptation of the M.tb-lux strain and thus we opted not to focus on this time period but can pursue further in the future. While the study group of 28 healthy AM donors was su cient to reveal signi cant genes associated with M.tb-AM interactions (this is a relatively large study given the di culty of obtaining fresh AMs from healthy individuals), supported by a smaller replication cohort, functional studies will be required to establish causal relationships. Also, larger cohorts will be needed to address contributions of single genes and networks to uptake, adaptation, and growth, reveal genetic factors, and assess ethnicity, sex and age, and environmental factors. Nevertheless, the study approach provides rich datasets to aid in facilitating development of predictive biomarker panels of an individual's risk for TB. The design of our study precluded repeat measures from each donor. However, our longstanding studies with healthy donor human MDMs has demonstrated that individual donors provide consistent results for comparable measurements over time.

Measurements of uptake, adaptation, and growth rate of M.tb in infected AMs
Fresh AMs were obtained, prepared and cultured within 6h from 28 tuberculin skin test (TST)-negative healthy donors from Caucasians, Asians, and Africans, according to the demographics of the Columbus, Ohio area (Table 1), under an approved IRB protocol at the Ohio State University Wexner Medical Center. Isolation and culture of human AMs from bronchoalveolar lavage (BAL) was done as described (53,54).
Brie y, BAL uid was centrifuged and washed once in cold RPMI at 4 0 C, and the cell pellet was resuspended in RPMI medium. A portion of the cell suspension was subjected to cytospin followed by staining and microscopy to determine macrophage content (94 ± 5%; mean ± SD, N = 10). AMs were adhered for 2h in either a 24-well plate (1.5x10 5 cells/well) or 96-well plate (5x10 4 cells/well) in RPMI containing 10% human AB serum and Penicillin G (10,000 U/ml) (~ 99% pure). The resultant uninfected HAM monolayers were then washed, and cultured in RHH media for another 2h (this is the time period where M.tb was added to separate wells for the infection). Supernatants were collected for protein and cell monolayers lysed for RNA analysis. For longer time points, monolayers were washed, medium replaced, and incubations continued for 24 or 72h for repeat protein and RNA analysis. AmpliSeq transcriptome analysis incorporates a targeted, amplicon-based (~ 110 bps, spanning exons) work ow, and is quantitative over orders of magnitude. The precision of AmpliSeq analysis detects M.tbinduced expression changes with high sensitivity in human MDMs and AMs infected with M.tb (30). Genomic DNA and total RNA (TRIzol® Reagent (Ambion™, Austin, TX)) were prepared from AMs using published procedures (11). RNA was puri ed, DNase-treated, RNA concentration measured, and RNA integrity assessed as described (30). After reverse transcription of 10 ng total RNA, using the AmpliSeq primers with the SuperScript® VILO™ cDNA Synthesis kit the cDNAs were ampli ed for 12 cycles with Ion AmpliSeq™ primers and barcoded adapters, resulting libraries puri ed and pooled in equal amounts for emulsion PCR on an Ion OneTouch™ 2 instrument, followed by sequencing with the Ion Proton™ sequencer (30). Reads were aligned to BED (Browser Extensible Data) le speci c for AmpliSeq amplicons. Typically, we obtain 5-9 million mappable reads per sample, with ~ 50-60% of RNA targets detected (30). Repeat experiments in the same sample yield correlation values of r 2 > 0.99 (for both independent replicates and sequencing chip replicates). The AM sample from D17 at 2h and 72h yielded < 1 million reads and were excluded from the analyses. The AmpliSeq reads were normalized to mapped fragments per million reads for quantifying transcript expression levels (57), yielding relative abundance for predicted transcripts in each AM.

Differentially-expressed (DE) genes between control and M.tb -infected human AMs
To identify DE genes, we employed DESeq2 (58). FDR adjusted p values of 0.05 were used as a cutoff for identifying DE genes at each time point. We estimated size factors using the "poscount" approach to correct for different sequencing depth and performed independent analysis at each time point, specifying the individual/donor and condition (infected or control) in the model.

Weighted gene co-expression network analysis) (WGCNA)
We derived three gene co-expression works using WGCNA (59) with the following parameters; power = 8, TOMtype = signed, minModuleSize = 30, mergeCutHeight = 0.25, minKMEtoStay = 0. These were blockwise consensus modules with the gene expression data from the control and infected cells grouped independently. We correlated the module eigengenes to either bacterial growth rate, or to proteins expressed/secreted from the same time point and performed FDR correction on p values.

Detection of genes with variable expression (VE) in control and M.tb -infected AMs, characterized by variance measures
To identify the most variably expressed (VE) RNAs in control AMs and in those separately after M.tb infection, all AM transcriptome data at each time were subjected to Levene's test (60), with ratios of variances as test statistics, reported adjusted p-value (FDR) for selected RNAs. A second variability test assesses whether the entropy of a given RNA's expression is higher than expected given the total entropy in the M.tb-treated AMs (61). A permutation test yields p-values for signi cance of the entropy computations.

Gene pathway and ontogeny analysis
We performed over representation analysis of gene ontology (GO) terms and reactome pathways using WebGestalt (62) taking the relevant gene lists for each comparison as an input and using default parameters. We retained GO terms or pathways which reached an FDR corrected p value of 0.05.
Comparison to genes associated with M1 and M2 states was performed using a curated list of genes from Viola et al.

Replication study of AMs in South Africa
AMs were collected from 11 close contacts of TB patients (n = 11) in Cape Town, South Africa (Table 2), under the approval of the Health Research Ethics Committee of Stellenbosch University. Close contacts were de ned as individuals who shared a closed space with a newly diagnosed TB patient for at least 5 hours per week (all contacts were QuantiFERON positive). These TB household close contacts were used as healthy controls. The processing of BAL uid to isolate AMs and their downstream uses were performed as described above for the main cohort from Columbus, Ohio with minor variations. AMs were adhered and cultured in 96-well plate (1.3x10 5 cells/well, in triplicate) in RPMI medium containing 20% AB serum. Cells were infected with wild-type M.tb H 37 R v at an MOI of 1:1. CFU, RNA and protein were measured at 2, 24 and 72h.    previously determined as operational in macrophage differentiation. Modules annotated as 'M1-like' or 'M2-like' are shown in blue and red, another major grouping of modules was related to ex vivo TPP (TNF, PGE2 and P3C) stimulation in (33) and are highlighted by a line.

Figure 4
Page 29/32 The transcriptional response of 28 human AMs to M.tb infection. Shown are the 36 signi cant DE genes common to all time points (see Figure 3A), and their average log 2 fold change.  Network of DE genes with highly variable RNA expression (VE/DE genes). These VE/DE genes (n=324) were identi ed in the 28 AMs after M.tb infection across all time points. Standard pathway enrichment program, Ingenuity Pathway Analysis (IPA) (https://www.qiagenbioinformatics.com/products/ingenuitypathway-analysis/), generates a top scoring gene network with IL1B, UBD, STAT1, and IRF1 as key hub genes. STAT1 and IRF1 cooperatively bind to the promoter of IDO1, which is also highlighted (fully annotated network is depicted in Additional le 8: Supplementary Fig. 3). These key gene names are highlighted in red.