To determine whether monocyte responses to Mtb include stimulation-specific eQTLs, we examined genome-wide genotyping (MEGA^EX genotyping array) and RNAseq transcriptional profiles of cells infected with H37Rv for 6h or uninfected controls in 101 participants, as part of a previously described Ugandan TB household contact study (14–16). We linked the genotypes and transcriptional profiles to define eQTLs among the 80 subjects from whom these datasets were available ( Figure 1 ). Although this primary study compared individuals who did (LTBI) or did not (“resister” or RSTR) convert their TST/IGRA after heavy exposure, we did not include these phenotypes in our primary analysis. Both the RSTR and LTBI groups were relatively young (mean age = 23 years), with no significant differences observed in sex, body mass index (BMI), or Bacille Calmette-Guérin (BCG) vaccination history, or epidemiologic exposure risk score between RSTR and LTBI groups ( Table 1 ) (17). Exposure risk scores were calculated based on index case exposure intensity (ranged 0-10), and both RSTR and LTBI groups had median risk scores of 6 (p= 0.316) (16).

To discover cis-eQTLs, we considered single nucleotide polymorphisms (SNPs) located 1 megabase (Mb) up- and down-stream of the transcription start site (TSS) and their association with gene expression in the Mtb infected and uninfected conditions. We observed 1,567 cis-eQTLs associated with the expression of 159 genes in Mtb-infected monocytes, whereas in uninfected monocytes, we identified 2,106 cis-eQTLs associated with the expression of 261 genes (false discovery rate [FDR] < 0.01, Supplementary Table S1 ). There was a significant overlap of eQTLs identified in infected and uninfected monocytes (n = 1,269; Wilcoxon signed rank test, p < 0.001). To discover “Mtb-dependent” cis-eQTLs, we next evaluated the interaction effect of genotypes between Mtb-infected and uninfected monocytes, while controlling for age and sex (FDR < 0.1). Among the 1,567 eQTLs identified in Mtb-infected monocytes, we found 32 eQTLs associated with the expression of 17 genes that were significant for an Mtb-infection:genotype interaction (FDR <0.1). Three of these eQTL were also identified under unstimulated conditions (FDR <0.01), and thus were excluded ( Figure 2A ). For further analysis, we focused on the remaining 29 eQTLs, referred to as “Mtb-dependent eQTLs”, which were associated with the expression of 16 genes ( Supplementary Table S2 , Supplementary Figure S1 ).

Among the 16 Mtb-dependent eQTL genes ( Table 2 , Supplementary Table S2 , Figure 2B , Supplementary Figure S1 ), 3 were associated with multiple eQTLs. Specifically, 12 eQTLs were associated with SIRPB1 expression ( Figure 3A ), 2 eQTLs with PRKAG2, and 2 eQTLs with NDUFAF4 ( Supplementary Figure S2 ). Notably, these eQTLs showed a high linkage disequilibrium (LD, r² > 0.8), indicating a potential single causative locus for each gene. Mtb infection resulted in clear genotype-dependent effects as compared to the uninfected condition ( Figures 3B–G ). In summary, our study identified 29 eQTLs that are associated with expression patterns of 16 genes in monocytes during Mtb infection, with genotype and condition-dependent patterns that vary in direction and magnitude.

To investigate the functional roles of 16 Mtb-dependent eQTL genes, we constructed a network of their protein-level interactions using the STRING database (18) and identified two high-confidence edges (STRING score > 700, Figure 4 ). To gain further insight into the functional roles of these genes, we performed hypergeometric mean enrichment pathway analysis of the 16 Mtb-dependent eQTL genes using the Human Molecular Signatures Database (MsigDB v2023.1, Supplementary Table S3 ) (19). In the C2 canonical pathway, we observed significant enrichment of gene sets related to the fatty acid and glutathione metabolism for GPX2, GSTO2, MLYCD, and PRKAG2 (FDR < 0.1). Additionally, within the C5 gene ontology biological pathways, we found associations of BMP6, CHMP1B, FANCA, MAS1, MLYCD, PRKAG2, and PIBF1 with 52 gene sets related to fatty acid, lipid, and organic hydroxy metabolism as well as mitotic spindle ( Supplementary Table S3 , FDR < 0.1). In summary, our network and enrichment analyses of the 16 Mtb-dependent eQTL genes revealed their involvement in multiple pathways associated with fatty acid metabolism, glutathione metabolism, and other cellular metabolic processes.

We examined the clinical significance of Mtb-dependent eQTL genes and variants by testing their association with resistance to TST/IGRA conversion in highly exposed HHCs of pulmonary TB subjects, a phenotype that could be regulated by macrophage responses. We conducted a candidate gene association study with a previously collected extended sample size (RSTR [case], n=74) versus (LTBI [control], n=189) for each Mtb-dependent eQTL gene, adjusting for age, sex, and kinship using a generalized linear mixed model (GENESIS) (11). One of the Mtb-dependent eQTLs, rs6599528 (associated with the cis-expression of CPSF1), was significantly associated with the RSTR phenotype where the adjusted odds of being RSTR increased by 1.63 times for each copy of the minor allele at this locus (p= 0.04, Table 3 ). We also examined whether SNPs within the up- and down-stream 5-kilobase (kb) flanking sequences of each of the 16 Mtb-dependent eQTL genes were associated with the RSTR outcome. We identified 26 SNPs in 7 genes associated with clinical RSTR status (p < 0.05; Supplementary Table S4 ), the majority of which were found in PRKAG2 (n= 20), with 6 having previously been reported (14, 20). Only these 6 previously reported SNPs in PRKAG2 remained significant after correction for multiple comparisons (FDR <0.2, 411 SNPs compared). Taken together, our results suggest that Mtb-dependent eQTLs and SNPs in their gene region may be associated with TB clinical outcomes.

Under inflammatory conditions, genetic polymorphisms of pro-inflammatory cytokines, such as tumor necrosis factor (TNF), interleukin (IL)-1β, IL-6, or IFN-β, have been associated with susceptibility to TB in humans (21–24). We examined whether Mtb-dependent eQTLs are associated with the expression levels of these cytokines. We identified 8 of the 29 Mtb-dependent eQTLs that were associated with Mtb-induced cytokine expression in monocytes (additive regression model adjusted for age and sex, p < 0.05, Figure 2A ). Among these, 4 eQTLs (rs55966428 [BMP6], rs7241199 [CHMP1B], rs1016193 [PRKAG2], rs12610448 [TIMM44]) were associated with altered IFNB1 gene expression in response to Mtb infection (p < 0.05, Figures 5A–C , 6C ). The remaining 4 eQTLs (rs247894 [MLYCD], rs7994794 [PIBF1], rs6599528 [CPSF1], rs6924573 [MAS1]) were significantly associated with altered IL1B, IL6 and/or TNF expression levels between Mtb-infected and uninfected monocytes (p < 0.05, Figures 5D–H , Supplementary Table S5 ). We also observed that 6 of the 8 eQTL genes additionally had variants that were associated with clinical RSTR outcomes ( Figure 2A ). These findings indicate that Mtb-dependent eQTL genes identify host genes that may regulate Mtb-induced cytokine induction in trans and potentially contribute to clinical resistance.

To further test our hypothesis that Mtb-dependent eQTL genes regulate cytokine pathways in monocytes, we used gene silencing methods to investigate potential mechanisms. Among 16 Mtb-dependent eQTL genes, we prioritized genes based on three criteria: (i) more pronounced slope changes between Mtb-infected and uninfected conditions, (ii) significant associations with cytokine gene expression, and (iii) significant associations with the clinical RSTR outcomes. From this analysis, BMP6, CHMP1B, KCNK5, PIBF1, PRKAG2, and TIMM44 were initially selected for gene silencing experiments. After assessing siRNA knockdown efficiency, we strategically narrowed our focus to BMP6, PRKAG2, and TIMM44 in subsequent cytokine signaling experiments, guided by their notable knockdown efficiency ( Supplementary Figure S3 ).

After silencing expression of BMP6, PRKAG2, and TIMM44, we stimulated PMA-differentiated THP-1 cells with DNA/RNA, toll-like receptor (TLR), or inflammasome ligands. We observed a decrease in IFNB1 expression after 4h of stimulation with 4 μg/mL sheared calf thymus DNA in BMP6-knockdown cells compared to the control group ( Figure 6A ). This expression pattern aligns with the HHC primary monocyte data, where the minor allele (G) of rs55966428 was associated with a significant reduction in the expression of both BMP6 ( Figure 6B ) and IFNB1 ( Figure 6C ) after 6h of Mtb infection. Although the other DNA/RNA ligands showed a similar trend of lower IFNB1 expression, the results were not statistically significant. In contrast, we did not find significant differences in IFNB1 expression for PRKAG2 or TIMM44 ( Supplementary Figure S4A ). There were no significant differences in IL-1β, IL-6, or TNF secretion observed between the BMP6, PRKAG2, and TIMM44-silenced cells and the siRNA controls ( Figures 6D–G , Supplementary Figures S4B-E ). Collectively, our data suggest that BMP6 regulates IFNB1 expression in response to Mtb infection in primary monocytes and DNA ligand stimulation in THP-1 cells, while no such regulation is observed in TLR- or inflammasome-induced cytokine expression in either cell type.

Between 2002 and 2012 (phase 1), we enrolled 872 culture-confirmed pulmonary TB index cases and their 2,585 HHCs in the Kawempe Community Health Study in Kampala, Uganda (92). HHCs were defined as individuals who had lived in the same household as the TB index case for at least 7 consecutive days in the preceding 3 months. At baseline and every 3-6 months thereafter for up to 24 months, all HHCs were evaluated for evidence of LTBI using TST. Between 2014 and 2017 (phase 2), a subgroup of these HHCs were re-contacted after 8-10 years of follow-up, and serial TST and IGRA were completed for TB screening when PBMCs were also collected. Detailed recruitment and screening procedures for this Ugandan HHC cohort have been previously described (11, 16, 92).

Cryopreserved PBMCs (n=115) were thawed and resuspended in RPMI 1640 medium (Gibco) supplemented with 10% heat-inactivated fetal bovine serum (FBS) and 50 ng/mL recombinant human monocyte colony-stimulated factor (M-CSF) for 24h. CD14⁺ monocytes were enriched from PBMCs by magnetic separation (Monocyte Isolation Kit II, Miltenyi Biotec) and incubated in RPMI-10 with M-CSF for 24h. Monocytes were then infected with H37Rv Mtb at a multiplicity of infection (MOI) of 1 or media alone (uninfected control) and incubated at 37°C with 5% CO² for 6h. After incubation, both the media-only and Mtb-infected monocytes were lysed in TRIzol (Invitrogen), and RNA was isolated using the RNeasy Mini Kit according to the manufacturers’ instructions (Qiagen). The RNA samples were quantified using Nanodrop 8000 instrument (Thermo Scientific), and their quality was measured by TapeStation (Agilent) (RNA Integrity Number ≥ 8.0). Next, cDNA libraries were prepared with random hexanucleotide primers and rRNA depletion using SMARTer RNAseq Kit (Takara), followed by sequencing on Illumina HiSeq 2500 and Novaseq 6000 platforms as previously described (15).

Sequences were processed as described previously (15). Briefly, sequences were aligned to the GRCh38 reference genome using STAR 2.6.0 (93) and quantified read counts using RSEM 1.3.0 (94). We narrowed our analysis down to 18,130 genes (68% of the initial 26,485 genes expressed in our monocyte samples), specifically focusing on those expressed in monocytes. We excluded 4,691 genes ( Figure 1 ) that were non-protein coding, not reliably annotated, or located on a sex chromosome. The latter would have complicated assessment of deviation from Hardy-Weinberg equilibrium (HWE) due to males being hemizygous. This left 13,439 genes, which we then trimmed-mean of M-values (TMM) normalized and converted to log2 counts per million (CPM) using limma (95). In parallel, genome-wide genotyping of 263 HHCs was performed using the Illumina MEGA^EX array. We successfully genotyped 1,002,430 SNPs, of which we filtered out 731,662 SNPs that deviated from HWE (P < 10^-6) and had a minor allele frequency (MAF) of less than 0.05.

To identify genetic variants that regulate mRNA expression in a cis-acting manner, we tested for associations between transcript expression levels and genotypes at SNPs located within a 1Mb window centered on the genes’ transcription start sites (TSS). As we lacked substantial power for trans analysis due to the small sample size, we focused solely on cis-eQTLs. Initially, we performed an additive linear regression of genotypes on log-transformed gene expression level, adjusted for age and sex in either Mtb-infected or uninfected monocytes. The R package MatrixEQTL was used to perform all regressions (96). We estimated FDR by using a Fisher’s exact test to correct for multiple testing by the Benjamini–Hochberg method to control for false positives. Significant eQTLs were defined as those with an FDR < 0.01 in Mtb-infected or uninfected monocytes.

For the 1,567 eQTLs identified in Mtb-infected monocytes (comprising 1,269 eQTLs identified in both Mtb-infected and uninfected monocytes, and 298 eQTLs identified in Mtb-infected monocytes), we further introduced an interaction term (Mtb-infection:genotype) in the main additive regression model to uncover whether the genetic effect of eQTLs on the level of gene expression differed by Mtb infection status. We applied a loose cut-off of FDR < 0.1 in the interaction model, given the small number of eQTLs included in this analysis. From the 32 eQTLs that showed significance in Mtb-infected condition and the interaction model, we removed 3 eQTLs that overlapped with those identified in uninfected condition to capture only those that regulate target genes in response to Mtb infection, which we named “Mtb-dependent eQTLs”. We defined the lead eQTLs per gene as the SNP that was most significantly associated (with the lowest FDR value) with the expression of that gene ( Figure 2 ). Finally, we assessed the genetic model fit across additive, dominant, and recessive models using the Akaike Information Criterion (AIC). The change in AIC was calculated for each eQTL by comparing the additive model to the dominant or recessive model (i.e., AIC of the additive model minus AIC of the dominant/recessive model; a negative value indicates that the additive model is a better fit, and vice versa). The change in AIC ranged from -59.31 to 7.41, with a median of -24 ( Supplementary Figure S6 ). The dominant or recessive model did not improve the model fit for the majority of eQTLs. Only 2 eQTLs exhibited enhanced fit with dominant and recessive models, with minimal changes in AIC (7.42 at rs55966428 and 3.89 at rs550516418). Considering the limited improvement of model fit from dominant or recessive models, we opted to utilize eQTLs identified through the additive model for our findings.

The associations between genotype and cytokine expression were assessed for each of the 16 lead Mtb-dependent eQTLs. Specifically, we examined the expression levels of IFNB1, IL1B, IL6, and TNF in Mtb-infected and uninfected monocytes, as genetic variation in these cytokine genes has previously been linked to TB susceptibility (21–24). We used a linear regression model to regress the number of minor alleles on the differences in cytokine expression between Mtb-infected and uninfected monocytes (normalized log2 [Mtb-infected – uninfected relative expression]), while adjusting for age and sex, at a significance level of 0.05.

We conducted association analyses to determine whether a particular eQTL genotype co-occurs with a RSTR clinical phenotype more often than would be expected by chance. We calculated odds ratios (ORs) using a quasi-likelihood approximation to the generalized linear mixed model (GLMM) (GENESIS, R package) (97). The model, incorporating adjustments for age, sex, and kinship, demonstrated minimal variations observed in close proximity to p=0.05. This included the loss of 2 SNPs in PRKAG2 and 1 in KCNK5, and the gain of 2 different SNPs in PRKAG2 and 1 in PIBF1. While no significant differences were identified, findings from the adjusted model were reported due to the biological relevance of age and sex in the context of TB risk. We calculated FDR values to adjust for multiple comparisons. Statistical comparison of demographic and epidemiologic data between RSTR and LTBI groups were conducted at a 2-sided significance level of 0.05 using Stata/IC 15.1 (StataCorp LLC) (98).

To examine potential protein-level interactions, we used the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING) database (18) against the 16 Mtb-dependent eQTL genes with strong interaction defined at a combined score > 700. We also conducted hypergeometric mean enrichment pathway analysis using the Molecular Signatures Database (MSigDB) C2 canonical pathways (BIOCARTA, KEGG, PID, REACTOME, WikiPathways) and C5 gene ontology biological pathways (BP) (19) with the R package clusterProfiler (99). Pathways were significant at FDR < 0.1 and overlap > 1. GO pathways were grouped by semantic similarity > 0.85 using rrvgo.

We used siRNA to silence Mtb-dependent eQTL genes, and 300 nM of two oligos per gene (siRNA identification s2032 and s2033 for BMP6; s3275a and s32752 for CHMP1B; s16450 and s16451 for KCNK5; s20481 and s20482 for PIBF1; s28111 and s28112 for PRKAG2; s20496 and s20497 for TIMM44, Life technologies) were used. As a control, Silencer Select negative control siRNAs (Silencer Select Negative Control No. 1 siRNA catalog # 4390843, Life technologies) was used at 300 nM per well. Transfection of the siRNAs into THP-1 cells was performed using Amaxa™ Cell Line Nucleofector™ Kit V (Lonza) following the manufacturer’s protocol for “Amaxa™ 4D-Nucleofector™ Protocol for THP-1 [ATCC].” After 2h, the nucleofected cells were replated in RPMI-10 (i) at a concentration of 3.0 × 10⁵ cells per well in a 24-well plate and differentiated for 24h in phorbol 12-myristate 13-acetate (PMA, Invitrogen) for DNA/RNA sensor stimulation experiments; (ii) at a concentration of 1.0 x 10⁵ cells per well in a 96-well pate and differentiated for 24h in PMA for stimulation with TLR ligands, NLRP3/NLRC4-specific proteins, as described below. Cells were washed and rested overnight in RPMI-10 without PMA at 37°C in a humidified incubator.

siRNA knockdown experiments were repeated for BMP6 with a 2nd transfection technique. THP-1 cells were plated at 5 x 10⁵ cells per well in a 24-well plate with RPMI-10 and PMA and rested at 37°C in a humidified incubator. After 24h incubation, differentiated cells were washed and replated with RPMI-10 without PMA. We used 10 nM each of the two oligos (s2032 and s2033) for BMP6 per well. As a control, Silencer Select negative control siRNAs (Silencer Select Negative Control No. 1 siRNA 4390843 and Silencer Select Negative Control No. 2 siRNA 4390846) were used at 10 nM each per well. Transfection of the pooled siRNAs into THP-1 cells was performed using Lipofectaime™ RNAiMAX transfection reagent (ThermoFisher) following the manufacturer’s protocol for “Universal Lipofectamin^® RNAiMAX Reagent Protocol 2013.” After 24h of treatment with BMP6 siRNA, cells were washed and replated with fresh RPMI-10 to each well, and treatment with DNA/RNA ligands was performed as described below.

For DNA/RNA sensing-specific stimulation, nucleofected THP-1 cells were treated with a medium and high dose of each DNA ligand, complexed with an equal dose of Lipofectamine 2000 for 4h. DNA ligands include: 1μg/mL and 4μg/mL of Sheared Calf Thymus DNA (Sigma Aldrich, #D1501) complexed with 1μg/mL and 4 μg/mL of Lipofectamine 2000 (Thermo Fisher Scientific); 0.1 μg/mL and 1μg/mL of Super Coiled Plasmid DNA (Invivogen) complexed with 2.5 μg/mL and 4 μg/mL of Lipofectamine 2000; 5 μM and 20 μM of cGAMP complexed with 0.4μg/mL and 4 μg/mL of Lipofectamine 2000; and 80μM cGAMP without Lipofectamine 2000. RNA ligands include: 5 μg/mL, 10μg/mL, and 20 μg/mL of poly (I:C) (Invivogen) complexed with 1 μg/mL, 2 μg/mL, and 4 μg/mL of Lipofectamine 2000, respectively. As controls, 10 ng/mL lipopolysaccharide (LPS, Invivogen); RPMI-10 complexed with an equivalent dosage of Lipofectamine 2000, RPMI-10 without Lipofectamine 2000 were used in each well. Following a 4h stimulation, cells were directly lysed in Buffer RLT Plus (Qiagen) and homogenized. RNA was isolated from lysates using the RNeasy^® Plus Mini Kit (Qiagen) following the manufacturer’s recommendations. cDNA synthesis was performed using the High Capacity cDNA RT kit (Applied Biosystems). IFNB1 and selected Mtb-dependent eQTL gene expression were measured by qRT-PCR using pre-designed TaqMan gene expression assays (Applied Biosciences).

For repeated BMP6 knockdown experiments, RNA was isolated from lysates using the same protocol with RNeasy^® Plus Mini Kit (Qiagen). Synthesis of the first strand cDNA was performed using SuperScript II reverse transcriptase and oligo (dT) primer (Invitrogen). qRT-PCR was performed with the CFX96 real-time system (Bio-Rad) using the SsoFast EvaGreen Supermix with the LOW ROX kit (Bio-Rad). The following primers designed from PrimerBank were used. The PrimerBank identifications are BMP6 (133930782c1) and IFNB1 (50593016c1).

BMP6 forward: AGCGACACCACAAAGAGTTCA

BMP6 reverse: GCTGATGCTCCTGTAAGACTTGA

IFNB1 forward: ATGACCAACAAGTGTCTCCTCC

IFNB1 reverse: GGAATCCAAGCAAGTTGTAGCTC

The expression of GAPDH was used to standardize the samples. For knockdown efficiency analysis, mRNA levels of siRNA-treated cells were normalized to control siRNA-treated cells using the 2−ΔΔCT (cycle threshold) method to calculate fold induction. For DNA/RNA sensing-specific stimulation, the results were expressed as the normalized ratio of each expression relative to the RPMI-10 with an equivalent dosage of Lipofectamine 2000 control.

For TLR-specific stimulation, nucleofected THP-1 cells were treated with LPS (100ng/mL TLR4, Invivogen), Mtb whole cell lysate (25 μg/mL NR-14822, BEI Resources), Pam2CSK4 (0.25μg/mL TLR2/6, Invivogen), and Pam3CSK4 (0.25μg/mL TLR2/1, Invivogen) for 24h. To induce NLRP3-specific activation, cells were treated with nigericin (10 μm) for 4h following a 2h pre-stimulation with LPS at 0.1 ng/mL. NLRC4-specific stimulation was achieved using recombinant proteins consisting of Burkholderia thailandensis needle protein combined with the Bacillus anthracis lethal factor (generously provided by Russel Vance, University of California Berkeley) for 4h (100). The supernatant from these TLR and inflammasome pathway experiments was harvested and analyzed for TNF-α, IL-6, and IL-1β cytokines using ELISA (R&D Systems).