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 (H37Rv) 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 (H37Rv). 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 reflecting the number of intracellular bacilli and microbial metabolic activity (29). RLU levels at 2-24h post-infection largely reflect 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 file 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 reflected 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 file 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, finding these to be broadly consistent with prior estimates of ~ 15h while varying two-fold between donors. To identify proteins and genes that influence M.tb-AM interactions, we focused on generation times (48-72h) in this study, since RLU levels prior to 48h post infection reflect 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 1
Donor demographics of AM samples used in this study (n = 28)
Donor | Sex | Age (Years) | Race | Donor | Sex | Age (Years) | Race |
D1 | Female | 24 | Asian (Chinese) | D15 | Male | 18 | Hispanic |
D2 | Male | 24 | African-American | D16 | Female | 25 | Caucasian |
D3 | Male | 27 | Caucasian | D17 | Female | 20 | Caucasian |
D4 | Female | 26 | Caucasian | D18 | Female | 26 | Caucasian |
D5 | Female | 26 | Caucasian | D19 | Male | 19 | African-American |
D6 | Male | 21 | Caucasian | D20 | Female | 22 | Hispanic |
D7 | Male | 24 | Caucasian | D21 | Male | 20 | Asian (Chinese) |
D8 | Female | 22 | Caucasian | D22 | Female | 19 | Not stated |
D9 | Male | 23 | Caucasian | D23 | Female | 19 | Caucasian |
D10 | Male | 21 | Caucasian | D24 | Male | 23 | Caucasian |
D11 | Female | 46 | Caucasian | D25 | Male | 24 | Hispanic |
D12 | Female | 22 | African-American | D26 | Male | 23 | Caucasian |
D13 | Male | 27 | Asian (Indian) | D27 | Female | 20 | Caucasian |
D14 | Male | 33 | Asian (Chinese) | D28 | Female | 19 | Caucasian |
Male/Female = 14/14; Caucasian = 60% (17/28), Asian = 14% (4/28), African-American = 11% (3/28), Hispanic = 11% (3/28) |
Human AM-secreted proteins capture the immune response to M.tb infection
To assess the inflammatory mediator response by macrophages to M.tb infection, we measured 27 secreted proteins, previously implicated in cellular responses to M.tb infection, in all AMs over 72h (Additional file 2: Table S2). At 2h, 6 of the 27 proteins (22.2%) were significantly 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. Inflammatory mediators including IL1β, TNF, and IL6 displayed significant differential secretion across the time points, indicating their contribution to a pro-inflammatory reaction to infection.
We then tested whether protein levels secreted from infected cells correlated with M.tb generation times (Fig. 1F-H; Additional file 3: Table S3). Three proteins were significantly 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 R2 = 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 significantly correlated with longer generation times (slower growth) (r = 0.46 and 0.49 at 2h and 24h, respectively). MMP10 levels were significantly correlated with slower growth at 24h, though switch to become positively (though non-significantly) correlated with rapid growth at 72h. These results support the notion that protein expression profiles over the first 3 days of infection influence M.tb-AM interaction dynamics.
Transcriptomes of uninfected control and M.tb -infected human AMs reveal differentially expressed gene profiles
Transcript profiles 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 (r2 ≥ 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 file 4: Supplementary Fig. 1). We previously showed that ex vivo cultured HAMs rapidly (starting at approximately 6h) begin to display the transcriptional profile 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 file 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 profiles 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 simplified to the trajectory towards M1-like or M2-like states, the former representing a pro-inflammatory state. To explore the relevance of these polarized states to differences between control and infected samples we first tested whether the phenotype of the cultured macrophages changed over time by deconvolution of the AmpliSeq profiles with CIBERSORT-X. Across both control and infected samples we saw a shift to an M0 macrophage population – likely reflecting the proposed shift of cells towards an MDM-like state. We identified 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 file 5: Supplementary Fig. 2).
To identify significant DE genes (FDR adjusted p ≤ 0.05) specific 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 file 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 significant at all-time points (Fig. 3A and Fig. 4), including those encoding inflammatory cytokines (IL1A, IL1B, TNF) and chemokine receptors (e.g., CCR7) characteristic of a spectrum representing “M1 type” cell states.
Our cell type deconvolution identified 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 (31) we confirmed a high overlap of DE genes at each time point with either M1 or M2 genes (Fig. 3B-E), accounting for 44.0% of M1 genes and 34.0% of M2 DE genes at 72h. While individual donors varied significantly in the expression of several characteristic M1 or M2 genes, significant correlations were not detectable between M1/M2 profiles and M.tb growth dynamics. While our data support the relevance of ‘M1’-like states to M.tb infection, polarization into either ‘M1’ or ‘M2’-like states is increasingly recognized as a false dichotomy with macrophages residing across a spectrum of states. For example, AMs express both M1 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 profiles, using gene co-expression modules previously identified from the full spectrum of macrophage differentiation (33). This showed a significant 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 significant 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 defined 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 inflammation (Fig. 3F). These results are consistent with AMs representing a unique macrophage phenotype that does not fit neatly into a classically defined M1 or M2 cell type but exists along a spectrum of macrophage types (Additional file 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 identified 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 file 7: Table S5A). All 324 VE genes are also significant DE genes which were found to be upregulated in and specific 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 file 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 influence on growth during the early phase of infection and requires further experimentation.
Reactome analysis yielded several significant pathways, including Immunological Diseases/Cellular Function and Maintenance/Inflammatory Response, and Function/Immune Cell Trafficking (Additional file 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 file 8: Supplementary Fig. 3). Both STAT1 and IRF1 interact with the IDO1 promoter (38) and play a role in macrophage polarization.
Identification of human AM gene networks associated with M.tb growth
Single RNA profiles across the donor AM dataset displayed high correlations with M.tb growth dynamics but failed to yield significance upon FDR corrections. For future studies, it is likely that increased sample size will be valuable in the identification 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 identification. 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 profiles 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 file 9: Table S6A). Correlation between network eigengene modules for infected cells and M.tb generation times (during 48-72h) at 2, 24 and 72h identified four significant 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 file 10: Table S7, columns D,E, modules highlighted in red). Reactome pathway over-enrichment analysis (Additional file 11: Table S8) showed Processing and Metabolism of RNA, and DNA Damage/Repair as significantly 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 significantly associated modules, we measured correlations between generation time and gene expression for each gene within the four significant modules (Additional file 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 significant 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 profiles and AM donors for genes within the yellow module from the 24h network, limiting to genes where we found a significant 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 identified three groups of AM donors, with groups 1 and 2 displaying significantly longer generation times compared to group 3 (Fig. 7B). The correlation coefficients (r2 values) were also determined from 24h gene expression data to assess their distribution across all modules (Fig. 7C).
Gene expression sub-groups displayed opposite profiles 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 identified the TNF signaling pathway as significantly enriched in negatively correlated genes (Fig. 7D), and oxidative phosphorylation and thermogenesis significantly 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 reflective 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 profiles 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 significantly 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 file 12: Table S9), identifying 29 significant correlations after Bonferroni correction. We opted for Bonferroni correction at this point for greater stringency in capturing correlations. Six of twenty-nine significant 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 file 11: Table S8 (column b) and Additional file 12: Table S9] and may be driven by ex vivo adaptation of macrophages.
In the 24h network, we found three modules with a significant link to M.tb growth rate and to both host transcriptional and protein profiles (Additional file 9: Table S6B and Additional file 10: Table S7). These are CCL3 (blue, red and yellow modules) and IL10 (red module). Significant 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 file 8: Supplementary Fig. 3) with 47 of 184 module eigengenes (25.4%) falling within this network, including IDO1, IL1A and IL1B.
Replication study of freshly isolated AMs in a different population using RNA-Seq
The goal of this part of the study was to test whether application of full RNA sequencing in AM samples from a different population in South Africa provides an overlapping set of DE genes in AM samples of our current study population. This population included 11 donors used as healthy controls with QuantiFERON positive status (Table 1 cohort A). M.tb H37Rv was the infecting strain. DE genes are presented in Additional file 13: Table S10 and Additional file 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 finding of substantial inter-individual differences, for example IDO1 and its co-expressed genes (Additional file 15: Supplementary Fig. 5).