Mycobacterium tuberculosis-dependent monocyte expression quantitative trait loci, cytokine production, and TB pathogenesis

Introduction The heterogeneity of outcomes after Mycobacterium tuberculosis (Mtb) exposure is a conundrum associated with millennia of host-pathogen co-evolution. We hypothesized that human myeloid cells contain genetically encoded, Mtb-specific responses that regulate critical steps in tuberculosis (TB) pathogenesis. Methods We mapped genome-wide expression quantitative trait loci (eQTLs) in Mtb-infected monocytes with RNAseq from 80 Ugandan household contacts of pulmonary TB cases to identify monocyte-specific, Mtb-dependent eQTLs and their association with cytokine expression and clinical resistance to tuberculin skin test (TST) and interferon-γ release assay (IGRA) conversion. Results cis-eQTLs (n=1,567) were identified in Mtb-infected monocytes (FDR<0.01), including 29 eQTLs in 16 genes which were Mtb-dependent (significant for Mtb:genotype interaction [FDR<0.1], but not classified as eQTL in uninfected condition [FDR≥0.01]). A subset of eQTLs were associated with Mtb-induced cytokine expression (n=8) and/or clinical resistance to TST/IGRA conversion (n=1). Expression of BMP6, an Mtb-dependent eQTL gene, was associated with IFNB1 induction in Mtb-infected and DNA ligand-induced cells. Network and enrichment analyses identified fatty acid metabolism as a pathway associated with eQTL genes. Discussion These findings suggest that monocyte genes contain Mtb-dependent eQTLs, including a subset associated with cytokine expression and/or clinical resistance to TST/IGRA conversion, providing insight into immunogenetic pathways regulating susceptibility to Mtb infection and TB pathogenesis.


Introduction
Despite ongoing efforts to eradicate tuberculosis (TB), it remains one of the leading infectious diseases worldwide (1).The wide spectrum of clinical states, ranging from asymptomatic latent TB infection (LTBI) to active TB disease, presents a challenge in identifying immune events providing protection early after exposure.While various inter-individual and environmental factors can influence TB susceptibility, such as bacillary load, proximity, duration of contact, age, malnutrition, and hygienic conditions, host genetics and immunologic factors are believed to play a significant role in increasing risk for TB disease (2,3).Therefore, it is critical to identify immunogenetic determinants of TB susceptibility to develop effective prevention strategies and provide insights for new TB treatments.
Mycobacterium tuberculosis (Mtb), the primary causative pathogen of TB, is an ancient infectious agent and has co-evolved with humans for over 10,000 years (4).Host and pathogen genetic pressure over this long time period selects genes and variants that may regulate protective responses in innate and adaptive immune cells.The primary niche of this evolutionary battle occurs within myeloid cells where Mtb resides in a phagosome.Such pressures can contribute to inter-individual variation in susceptibility to Mtb infection and TB disease (5); however, the specific major susceptibility genes and variants remain largely unknown.The search for the genetic determinants of clinical TB susceptibility has included case-control studies of candidate genes and genomewide association studies (GWAS) (6)(7)(8)(9)(10)(11).Despite extensive efforts, the major genes and specific functional effects of the polymorphisms associated with TB remain poorly understood.One approach to address this knowledge gap is expression quantitative trait loci (eQTL) mapping (12).eQTL mapping allows for the identification of genetically programmed genes responding to Mtb, regulating essential antimicrobial and inflammatory responses that contribute to Mtb clearance and susceptibility to Mtb infection.Given the crucial role of myeloid cells in phagocytosis, antigen presentation, and cytokine production during early Mtb infection (13), associations between genetic variants and Mtb-dependent changes in intermediate traits within myeloid cells, particularly cytokine gene expression, may provide insights into the underlying immunogenetic mechanisms of Mtb infection and TB disease.
To investigate the functional consequences of genetic variants within macrophages responding to Mtb infection, we examine our hypothesis that a specific subset of monocyte genes has eQTLs that regulate monocyte function including pro-inflammatory cytokine responses to Mtb or that are associated with clinical TB phenotypes.Through genome-wide genotyping and transcriptional profiling of Mtb-infected and uninfected CD14 + monocytes from 80 household contacts (HHCs) of pulmonary TB in Uganda, we identified 29 eQTLs associated with the expression of 16 genes under Mtb-stimulated and not unstimulated conditions.Some of these "Mtb-dependent" eQTLs were also associated with Mtb DNA sensor-dependent interferon-b (IFN-b) expression, resistance to tuberculin skin test (TST) and interferon-g (IFN-g) release assay (IGRA) conversion after Mtb exposure, and modulation of immunometabolic pathways.These findings shed light on the immunogenetic factors underlying TB susceptibility and hostpathogen interactions in the context of Mtb infection.

Monocytes express 16 genes with Mtb-dependent eQTLs
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)(15)(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-Gueŕin (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 downstream 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 2 > 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.

Network and enrichment analyses implicate Mtb-dependent eQTL genes in fatty acid and glutathione metabolism, along with other cellular metabolic processes
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 highconfidence 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 Mtbdependent eQTL genes revealed their involvement in multiple pathways associated with fatty acid metabolism, glutathione metabolism, and other cellular metabolic processes.Study Design and eQTL Mapping Data Processing in the Uganda Household Contact Study.Household contacts (HHCs) of pulmonary TB index cases were studied, enrolled, and followed longitudinally for assessment using serial TST and IGRA testing.Laboratory studies included RNAseq measurements in Mtb-infected monocytes and genome-wide genotyping with the MEGA EX genotyping array.eQTL, expression quantitative trait loci; HHC, household contact; HWE, Hardy-Weinberg equilibrium; IGRA, interferon-g release assay; LTBI, latent tuberculosis infection, MAF, minor allele frequency, MOI, multiplicity of infection; Mtb, Mycobacterium tuberculosis; PBMC, peripheral blood mononuclear cell; RNAseq, RNA sequencing; RSTR, resister or resistant to TST/IGRA conversion after high TB exposure; SNP, single nucleotide polymorphism; TST, tuberculin skin test.Dotted arrows indicate the exclusion of samples or data for the listed reasons.

Polymorphisms in Mtb-dependent eQTL genes and association with resistance to TST/IGRA conversion in TB household contacts
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 5kilobase (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.

Mtb-dependent eQTLs associate with monocyte cytokine expression
Under inflammatory conditions, genetic polymorphisms of pro-inflammatory cytokines, such as tumor necrosis factor (TNF), interleukin (IL)-1b, IL-6, or IFN-b, have been associated with susceptibility to TB in humans (21)(22)(23)(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 Mtb-dependent eQTL mapping: (A) Flowchart of Mtb-dependent eQTL identification.Using an additive regression model (FDR< 0.01), we identified 2,106 eQTLs in uninfected monocytes and 1,567 eQTLs in Mtb-infected monocytes after adjusting for age and sex.Mtb-dependent eQTLs were specifically defined as cis-eQTLs that showed strong statistical evidence only in Mtb-infected monocytes.We performed additive linear regression with an interaction term for the Mtb stimulation status to evaluate the interaction effect at an FDR of 10%.Out of the 2,106 eQTLs in Mtb-infected monocytes, 32 eQTLs showed significant linear associations with gene expression levels after adjusting for age and sex (FDR <0.1).Three eQTLs were excluded as they were identified in both Mtb-infected and uninfected monocytes.We identified 29 Mtbdependent eQTLs regulating the expression levels of 16 genes.Additionally, 8 Mtb-dependent eQTLs were associated with altered levels of IFNB1, IL1B, IL6 or TNF expression in trans (orange circle).regulate Mtb-induced cytokine induction in trans and potentially contribute to clinical resistance.

BMP6 expression is associated with Mtb and DNA-ligand induced IFNB1 expression in monocytes
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 Mtbdependent 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, tolllike receptor (TLR), or inflammasome ligands.We observed a decrease in IFNB1 expression after 4h of stimulation with 4 mg/ 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-1b, 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.

Discussion
Myeloid cells play a critical role to restrict Mtb replication through intrinsic microbicidal pathways and by stimulating downstream cellular responses, but the genetic factors that The expression profiles of these Mtb-dependent eQTL genes and their associated pathways indicate the involvement of diverse inflammatory and cellular metabolic processes, potentially influencing the mechanisms underlying resistance or susceptibility to Mtb infection.Our primary findings suggest that BMP6 (bone morphogenetic protein 6), a member of the transforming growth factor (TGF)-b superfamily, is involved in the genetic regulation of IFN-b and may influence antimicrobial responses.BMP6 plays a critical role in tissue remodeling and organogenesis (25,26).The TGF-b and BMP pathways have an antagonistic relationship due to shared receptor structures and signaling mechanisms (27,28), which have been observed in conditions like pulmonary fibrosis (29, 30) and malignancies (31)(32)(33).In myeloid cells, BMP6 distinguishes itself with its unique pro-inflammatory function, enhancing antimicrobial activities, while other BMPs (BMP2, BMP4, and BMP7) contribute to anti-inflammatory functions and the promotion of the M2 macrophage phenotype (34).In murine models, BMP6 treatment activates NF-kB signaling, induces IL1B and TNF expression, and increases the expression of inducible nitric oxide synthase (iNOS), resulting in elevated nitric oxide production and enhanced microbicidal effects (35,36).In human monocytes, Ma et al. found that BMP6 expression is stimulated by extracellular HSP70 (an endogenous TLR4 ligand), which is released from LAMP3overexpressing salivary gland epithelial cells in Sjogren's syndrome (37).Although the direct association between BMP6 expression and    polyadenylation specific factor 1) is crucial for mRNA maturation and alternative 3' untranslated region isoform generation (45).The involvement of CPSF1 in various human diseases, including cancers (46)(47)(48) and ocular disorders (49, 50) has been studied, but its role in TB susceptibility remains unexplored.Importantly, our genetic analysis revealed that the minor allele (C) at rs6599528 was associated with a 1.63-fold increase in the odds of being RSTR, potentially indicating a protective property of the C allele against TST/IGRA conversion.Upon Mtb infection, the majority of individuals who carry the major allele (A) has decreased expression of CPSF1; while no changes in its expression among individuals with recessive genotype (C/C).These findings suggest that the genotypic regulation of CPSF1 via rs6599528 may potentially mediate alternative splicing of host transcripts that influence susceptibility to Mtb infection.Among potential host pathways that mediate this protection, recent attention has focused on the alternative role of CPSF1 as an E3 ubiquitin ligase that degrades the transcription factor hypoxia-inducible factor (HIF)-1a (51).HIF-1a is essential for the IFN-g-driven macrophage immunometabolic response (52,53) and its activation results in increased IL6 expression (54)  result in HIF-1a degradation, is consistent with IFN-g-independent protective mechanisms (55,56).Interestingly, in our data, IL6 expression is elevated in resting monocytes carrying the minor allele (Supplementary Figure S5G).However, the genetic effects of rs6599528 on IL6 expression become non-significant following Mtb infection, indicating that IL-6 regulation is influenced by more complex and multifaceted factors.Therefore, further investigations are needed to elucidate the complex influence of rs6599528 on proinflammatory signaling, antimicrobial mechanisms against Mtb, and clinical resistance outcomes.Mtb has evolved strategies to manipulate macrophage metabolism by utilizing host fatty acids and cholesterol as its primary energy sources (57)(58)(59).We identified 5 Mtb-dependent eQTL genes, including PRKAG2, MLYCD, BMP6, PIBF1, and GPX2, that are involved in fatty acid and lipid metabolism including GO pathways clustered in purine-containing processes (Supplementary Table S3).The heterotrimeric AMP-activated protein kinase (AMPK) contains one kinase (AMPK-a) and two regulatory domains (AMPK-b and AMPK-g), which upon activation promotes catabolic processes such as fatty acid oxidation (FAO) (60,61) by phosphorylation of central metabolic enzymes (62).Genetic variants in PRKAG2, which encodes one AMPK-g isoform, are associated with the RSTR phenotype (14).Moreover, we identified a conserved enrichment of gene sets associated with carbon metabolism and the transcriptional response to free fatty acids (FFAs) across independent RSTR cohorts in Uganda and South Africa (14).Although the PRKAG2 Mtb-dependent eQTLs identified here are distinct from PRKAG2 polymorphisms that are associated with the RSTR phenotype, our findings further support a role for PRKAG2 and AMPK in measurable in vitro and clinical TB outcomes that is under selective pressure.To further support host lipid metabolism during monocyte Mtb responses, the MLYCD Mtbdependent eQTL indicates host-genotype dependent responses to Mtb that may impact lipid synthesis.Malonyl-CoA decarboxylase (MCD), encoded by MLYCD, catalyzes the decarboxylation of malonyl-CoA (63)(64)(65).Since malonyl-CoA is both an intermediate during lipid synthesis and an allosteric inhibitor of carnitine palmitoyltransferase 1 (CPT1)-mediated FAO, MCD activity results in the inhibition of lipid synthesis and the promotion of FAO.Although the direct impact of MLYCD downregulation in response to Mtb infection has not been studied, it is likely to disrupt FAO.Consistent with the predicted effect that MLYCD downregulation (and thus malonyl-CoA accumulation) may have towards Mtb virulence, Mehrotra et al. ( 66) measured increased malonyl-CoA levels following Mtb infection that was specific to virulent strains and resulted in net fatty acid synthesis.Further inhibition of FAO induces HIF-1a and excessive production of reactive oxygen species (ROS) from the mitochondria (67).As a result, the excessive HIF-1a and ROS promote the transcription of glycolytic enzymes and inflammatory mediators, including IL-1b, IL-6 and TNF, facilitating additional antimicrobial responses (68,69).Our data demonstrate that monocytes with homozygous dominance at rs247894 exhibit reduced MLYCD expression in cis but increased IL1B and TNF expression in trans in response to Mtb infection.These findings suggest that the downregulation of MLYCD contributes to the host protective pro-inflammatory properties.However, reduced FAO also leads to the accumulation of lipids that can serve as a nutrient source for Mtb growth (70,71).The conflicting associations observed in MLYCD suggest an ongoing evolutionary battle within host cells, necessitating further investigation into its role in immunometabolic function and host-protective response.Paradoxically, macrophages also employ redox regulation mechanisms, primarily reliant on glutathione metabolism, to prevent oxidative damage caused primarily by NADPH-oxidase-driven hydrogen peroxide (72).Our findings support the involvement of glutathione-mediated antioxidant enzymes, GPX2 (glutathione peroxidase 2) (73-75) and GSTO2 (glutathione S-transferase omega 2) (76,77), in redox regulation during Mtb infection.These enzymes, through genetically regulated expression that is induced by Mtb, may influence antioxidant defense or bacterial clearance (78,79).Together, our study highlights the complex interplay between immunometabolic pathways, genetic factors, and Mtb infection outcomes, emphasizing the necessity for further research to elucidate the roles of Mtb-dependent eQTL genes and their impact on immunometabolic function and host protection in the context of Mtb infection.
Prior studies of Mtb-dependent eQTLs illuminate potential areas of convergence of findings.Previously, Barreiro et al. (80) identified 102 eQTLs in Mtb-uninfected and 96 in Mtb-infected monocyte-derived dendritic cells (DCs) using a similar study design with some significant differences, including cell type (DCs versus monocytes), age (adults versus range of ages), sex (male versus both male and female), population background (White versus Ugandan), time point (18h versus 6h), and transcriptomic method (expression array versus RNAseq).In addition, different statistical criteria were used to define "response eQTLs" (FDR significance in one condition but not the other versus Mtb-infection:genotype interaction model) and distance of cis-regulatory effects (200kb proximity from the TSS, versus 1Mb).Despite no identical eQTLs or eQTL genes being found, both studies found several genes within the same family.For example, BMP1, a response eQTL gene in uninfected DCs, shared the same gene family with our Mtb-dependent eQTL gene, BMP6.Unlike the rest of the BMP family, BMP1 does not belong to the TGF-b superfamily; instead, it acts as a procollagen C-proteinase in collagen maturation and indirectly contributes to TGF-b upregulation by activating matrix metalloproteinase 2 (MMP2) and MMP9, influencing tissue remodeling (81,82).In addition, Mtb-infected DCs exhibited CPSF3 and NDUFAF2 as response eQTL genes, belonging to the same gene family as our Mtbdependent eQTL genes, CPSF1, and NDUFAF4, respectively.Little research has focused on CPSF3, but like CPSF1, it is thought to be involved in mRNA processing and may contribute to the regulation of gene expression in response to diverse pathogenic microorganisms, such as Toxoplasma gondii (83), Plasmodium falciparum (84), and Cryptosporidium (85) infection.Finally, NDUFAF4 plays an essential role in the assembly and maturation of the NADH:ubiquinone oxidoreductase complex (complex I) in the mitochondrial respiratory chain (86,87), linked with the latestage assembly factor NDUFAF2 (88, 89), crucial for energy production through oxidative phosphorylation (90).The convergence of these findings adds weight to the results and encourages further exploration of shared genes or pathways.
This study has several limitations, primarily stemming from the constraints imposed by a small sample size and limited availability of primary cells from the HHC cohort.The "Mtbdependent eQTL" identified may lack specificity to Mtb infection since our stimulation condition did not include other pathogens or ligands, an area for future investigation.The limited availability of peripheral blood mononuclear cells (PBMCs) precluded using multiple time points to understand dynamic transcriptional changes during various stages of Mtb infection.Furthermore, the candidate gene association analysis with the RSTR phenotype is underpowered.Additionally, our findings using peripheral blood monocytes may not extend to other tissues and cell types at the site of Mtb pathogenesis such as alveolar macrophages where variants or Mtb-dependent eQTL genes may be distinct.The supernatants for protein analysis were not collected from HHC samples, limiting the assessment of the impact of Mtb-dependent eQTLs and SNPs on cytokine release.Finally, although our siRNA silencing experiments and cytokine gene association analysis suggest a plausible causal effect that Mtb-dependent eQTL gene expression has on cytokine responses following Mtb-infection, future advanced functional assays will explore the regulatory roles of eQTL promoter/enhancer elements that mediate these responses.Interestingly, Souza de Lima et al. (91) identified that the imbalance of IL-1b and IFN-a is associated with TB resistance in TB-endemic areas, suggesting a need for a multifaceted approach to various cytokines and their interconnection.In future investigations, considering the inclusion of other cytokines, such as type I IFNs beyond IFN-b, as well as IL-18 and other members of the IL-1 family, along with their secretion levels, will be crucial to enhance the specificity of Mtb-dependent eQTL genes and variants in response to Mtb infection.
In summary, this study demonstrates the robust application of genome-wide eQTL mapping to investigate genetic regulation of the host response against Mtb.Our findings highlight significant associations between Mtb-dependent eQTL genes and variants and alterations in the immune response to Mtb infection, as well as potential regulatory interactions between specific eQTLs and Mtb-induced cytokines.These insights into immunogenetic mechanisms provide valuable targets for host-directed therapies and TB control.

Overview of study design and sampling
Between 2002 and 2012 (phase 1), we enrolled 872 cultureconfirmed 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).

CD14 + Monocyte isolation, Mtb infection, and RNA sequencing
Cryopreserved PBMCs (n=115) were thawed and resuspended in RPMI 1640 medium (Gibco) supplemented with 10% heatinactivated 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 2 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).

Sequence data processing
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 trimmedmean 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.

eQTL analysis
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 (Mtbinfection: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 cutoff 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 "Mtbdependent 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.

Mapping cytokine expression levels
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)(22)(23)(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-infecteduninfected relative expression]), while adjusting for age and sex, at a significance level of 0.05.

Mtb-dependent eQTL gene association study
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).

Integrative pathway enrichment analysis
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.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 5 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 5 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.
4.9 Cytokine expression and secretion following stimulation of cytoplasmic DNA/ RNA sensing, TLR, and inflammasome pathways 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: 1mg/mL and 4mg/mL of Sheared Calf Thymus DNA (Sigma Aldrich, #D1501) complexed with 1mg/mL and 4 mg/mL of Lipofectamine 2000 (Thermo Fisher Scientific); 0.1 mg/mL and 1mg/ mL of Super Coiled Plasmid DNA (Invivogen) complexed with 2.5 mg/mL and 4 mg/mL of Lipofectamine 2000; 5 mM and 20 mM of cGAMP complexed with 0.4mg/mL and 4 mg/mL of Lipofectamine 2000; and 80mM cGAMP without Lipofectamine 2000.RNA ligands include: 5 mg/mL, 10mg/mL, and 20 mg/mL of poly (I:C) (Invivogen) complexed with 1 mg/mL, 2 mg/mL, and 4 mg/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 −DDCT (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 mg/mL NR-14822, BEI Resources), Pam2CSK4 (0.25mg/ mL TLR2/6, Invivogen), and Pam3CSK4 (0.25mg/mL TLR2/1, Invivogen) for 24h.To induce NLRP3-specific activation, cells were treated with nigericin (10 mm) 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-a, IL-6, and IL-1b cytokines using ELISA (R&D Systems).informed consent for participation in this study was provided by the participants' legal guardians/next of kin.Ethical approval was not required for the studies on animals in accordance with the local legislation and institutional requirements because only commercially available established cell lines were used.

Finally, 7
Mtb-dependent eQTL genes contained SNPs associated with the clinical RSTR phenotype (blue circle).Importantly, 6 Mtb-dependent eQTL genes and variants were associated with both clinical RSTR phenotypes and the regulation of cytokine gene expression.(B) Manhattan plot of Mtb-dependent eQTLs.Significance of eQTLs in the Mtb-infected condition are plotted across the genome, and the horizontal dashed line represents the significance threshold (FDR = 0.01).Red dots mark lead eQTLs that are also significant in the interaction model (Mtb-infection:Genotype, FDR > 0.1), yet not in the uninfected condition (FDR > 0.01); these are labeled with their respective cis genes.Grey dots positioned above this line represent eQTLs in Mtb-infected conditions (FDR < 0.01), though they do not demonstrate significance in the primary interaction model.eQTL, expression quantitative trait loci; FDR, false discovery rate; Mtb, Mycobacterium tuberculosis; SNP, single nucleotide polymorphism.

3
FIGURE 3 Mtb-dependent eQTL Monocytes Expression Plots.eQTLs were defined by significance in the Mtb condition (FDR < 0.01) and interaction term (FDR < 0.1) models but not in uninfected condition.(A, top) eQTL significance for SNPs within 1 Mb of SIRPB1.X-axis indicates position in chromosome 20.Horizontal dashed line indicates FDR = 0.1 in the interaction model and the lead SNP is labeled.(A, bottom) Heatmap indicating R 2 linkage disequilibrium (LD) for SNPs in this region.(B-G) eQTL plots illustrate the relationship between genotype (x-axis) and log2 expression (y-axis) of the indicated genes.Each line indicates a linear fit derived from an additive regression model that includes an interaction term (Mtb_infection : Genotype), adjusting for age and sex.FDRs for each eQTL are reported for Mtb infected and uninfected conditions, not the interaction model, and are presented to illustrate the effects of Mtb stimulation and host genotype.The lead eQTLs (B) rs6136375 for SIRPB1 and (C) rs6599528 for CPSF1 may have minor genotype effects independent of Mtb stimulation (i.e., FDR <0.1 in the Mtb-uninfected condition).(D) rs247894 for MLYCD highlights a locus where no genetic effect was observed in the uninfected condition (FDR > 0.1), but becomes significant upon Mtb infection.(E) rs940746 for GSTO2, (F) rs7144688 for GPX2 and (G) rs1016193 for PRKAG2 highlight loci where the directionality of genotype effects differed between Mtb-infected and uninfected conditions.ALT, alternative allele; eQTL, expression quantitative trait loci; FDR, false discovery rate; Mb, megabase; Mtb, Mycobacterium tuberculosis; REF, reference allele; SNP, single nucleotide polymorphism.

FIGURE 4
FIGURE 4Network and Gene Enrichment Pathway Analysis.The left-side box in the figure illustrates the STRING network of protein-protein interactions among 16 Mtb-dependent eQTL genes.Each node denotes a Mtb-dependent eQTL gene and each line indicates a protein-level interaction.Nodes are colored by hypergeometric mean enrichment pathway analysis with MSigDB canonical pathways and gene ontology.Significant pathways are colored by groups determined by semantic similarity (threshold = 0.9) and labeled with the parent term determined by rrvgo.PRKAG2, MLYCD, and GPX2 are functionally related in regulation of fatty acid oxidation and metabolic process.GPX2 and GSTO2 are functionally related in glutathione metabolism and cellular oxidant detoxification.GOBP, Gene Ontology Biological Process; KEGG, Kyoto Encyclopedia of Genes and Genomes; STRING, Search Tool for the Retrieval of Interacting Genes/Proteins.

5
FIGURE 5 Mtb-dependent eQTLs associated with Mtb-induced cytokine expression in monocytes.Three Mtb-dependent eQTLs (A-C), rs7241199, rs1016193, and rs12610448, were associated with the expression of CHMP1B, PRKAG2, and TIMM44 in cis, respectively, and with the changed IFNB1 expression upon Mtb infection.(D) rs247894 and (E) rs7994794 were associated with the cis-regulation of MLYCD and PIBF1 expression, respectively, as well as with the regulation of IL1B.(F) rs6599528 and (G) rs6924573 were associated with the cis-regulation of CPSF1 and MAS1 expression, respectively, along with the regulation of IL6.(H) rs247894 was associated with the cis-regulation of MLYCD and the regulation of TNF expression.The slope of the lines indicates the ratio of the effect estimate to the standard error derived from a linear regression model, which estimates the relationship between the number of minor alleles and the differences in cytokine expression (normalized log2 [Mtb-infecteduninfected relative cytokine expression]), after accounting for age and sex (p < 0.05).The gene names in parentheses beside the rsID on the X-axis title indicate the target gene of each eQTL.ALT, alternative allele; REF, reference allele; Mtb, Mycobacterium tuberculosis.

6 BMP6
FIGURE 6    BMP6 expression is associated with host genotype in Mtb-infected cells and is required for maximal IFNB1 responses in monocytes following DNAligand stimulation.(A) The fold change expression of IFNB1 was significantly reduced in BMP6-silenced THP1 cells after 4h of stimulation with 4 mg/ mL calf thymus DNA ligand, compared to siRNA control cells.BMP6-silenced cells stimulated with 1 mg/mL of super coiled plasmid DNA, 20 mM of cGAMP, or 10 mg/mL of poly(I:C) exhibited statistically non-significant decreases in IFNB1 induction compared to siRNA control cells.The fold change expression of IFNB1 was measured by qRT-PCR and normalized against the IFNB1 expression from lipofectamine without DNA/RNA ligand to control background IFNB1 induction from lipofectamine.(B) In the Ugandan HCC cohort, the minor allele [G] of rs55966428 was associated with decreased expression of BMP6 following 6h of Mtb infection in primary monocyte data from HHCs.Host genotype (x-axis) and normalized log2 expression (y-axis) of BMP6 for the rs55966428 Mtb-dependent eQTL with false discovery rate (FDR) reported separately for Mtb-infected and uninfected conditions.(C) In the same cohort, the minor allele [G] of rs55966428 was associated with decreased IFNB1 expression following 6h of Mtb infection in primary monocyte data.This plot shows genotype on the x-axis and differences in IFNB1 expression on the y-axis (normalized log2 [Mtb-infecteduninfected relative IFNB1 expression]), adjusting for age and sex.To assess (D) inflammasome-mediated IL-1b response, nucleofected cells were primed with LPS for 2h and treated with nigericin for NLRP3-specific stimulation, Burkholderia thailandensis needle protein for 4h for NLRC4-specific stimulation.Needle protein was administered at 8 ng/ml in conjunction with 16 ng/ml Bacillus anthracis protective antigen.The IL-1b supernatant levels were measured by ELISA.For TLR-specific stimulation, BMP6-silenced and siRNA control THP-1 cells were stimulated with LPS, PAM2/PAM3, Mtb whole cell lysate, or media for 24h.(E) IL-1b, (F) IL-6, and (G) TNF supernatant levels were measured using ELISA.cGAMP, cyclic guanosine monophosphate-adenosine monophosphate; CT, calf thymus DNA; LPS, lipopolysaccharide; media, media control (RPMI + 10% FBS); PAM2, Pam2CSK4; PAM3, Pam3CSK; poly(I:C), polyinosinic acid-polycytidylic acid; SC, super coiled plasmid DNA; TBWCL, Mtb whole cell lysates; TLR, toll-like receptor; KD, knockdown; Mtb, Mycobacterium tuberculosis.**p < 0.01 tested by paired t-test, ns, not significant.

TABLE 1
Demographic characteristics of participants.
An additive linear regression model, including an interaction term (Mtb_infection : Genotype), is employed to calculate the t-statistics of eQTLs.The selection of Mtb-dependent eQTLs is based on FDR values to account for false positives resulting from multiple testing.eQTLswithan FDR value below 0.1 are considered statistically significant Mtb-dependent eQTLs, and only the lead eQTLs are listed in this table.FDR, false discovery rate; MAF, minor allele frequency; SNP, single nucleotide polymorphism.aGRCh37.bTheFDR values of the interaction model.

TABLE 3
Association of CPSF1 with clinical resistance to TST/IGRA conversion.
(42)(40)(41)ve model, the odds ratio (OR) measures the change in the odds of being RSTR with each additional copy of the minor allele.The R 2 value of linkage disequilibrium (LD) from Mtb-dependent eQTL to the designated SNP within the Mtb-dependent eQTL gene.IGRA, interferon-g release assay; kb, kilobase; LTBI, latent tuberculosis infection; OR, odds ratio; p, p-value; RSTR, resister or resistant to TST/IGRA conversion after high TB exposure; SNP, single nucleotide polymorphism; TST, tuberculin skin test.activation of cytosolic cGAS-or RLR-mediated sensing pathways triggers a robust type I IFN response that impairs host resistance to Mtb infection(39)(40)(41).Conversely, activation of other cytosolic pathways during Mtb infection, such as AIM2, NOD2, and NLRP3, promotes the production of protective inflammatory cytokines(42).investigation is needed to explore BMP6's pro-inflammatory function and its regulatory dynamics with IFN-b during Mtb infection in human innate immune cells.We also found evidence suggesting Mtb may regulate monocyte gene expression through mRNA processing.CPSF1 (cleavage and