Epigenetic impairment and blunted transcriptional response to Mycobacterium tuberculosis of alveolar macrophages from persons living with HIV

Study subjects

To investigate the influence of HIV on the ability of AM to mount a response to Mtb infection, we studied three groups of subjects: PLWH receiving ART (n=20), persons without HIV but receiving ART as pre-exposure prophylaxis (PrEP; n=14) and persons without HIV not receiving PrEP (healthy controls, HC; n=16, Table 1). We obtained AM by BAL from all subjects and challenged the cells in vitro with Mtb. RNA and DNA were isolated from uninfected and Mtb challenged cells and used to assess gene expression levels by RNA-seq, chromatin structure by ATAC-seq, and presence of the histone activation mark H3K27ac by ChIP-seq (Fig. 1A). However, due to different quality-control procedures and variable numbers of AM obtained per subject, different counts of subjects were used for the transcriptomics and epigenetics experiments (Fig. 1A).

• Download figure
• Open in new tab

Fig. 1 Study design and results for differential gene expression by alveolar macrophages (AM) after challenge with Mycobacterium tuberculosis (Mtb).

(A) Experimental design and sample number by group. Each subject underwent bronchoalveolar lavage (BAL), AM were obtained from BAL and challenged in-vitro with Mtb. RNA was obtained for RNA-seq, while DNA was used for ATAC-seq and ChIP-seq experiments. (B) Bar graph summarizing differentially expressed genes (DEG) by AM following Mtb challenge across the three phenotypic groups (healthy controls – HC; persons living with HIV on antiretroviral therapy – PLWH; subjects on pre-exposure prophylaxis with antiretroviral therapy - PrEP). The y-axis indicates the cumulative DEG count and group ID is indicated below each bar. Dark colour shades represent up-regulated genes (positive log fold-change; logFC) while light shades depict down-regulated genes (negative logFC). (C) Density plot presenting logFC as absolute values. All DEGs from panel B had their log2FC converted to absolute values and plotted using density function. Vertical coloured lines indicate the median of absolute values per group, with exact values given on top. Blue shade represents HC, red indicates PLWH and PrEP were coloured in green. (D) Dot plot presenting significance and overall effect for selected pathways and GO-terms. Pathways/GO-terms are listed on the left, dot sizes reflect the negative log10 FDR values for the enrichment test (bigger dots have smaller p-values) and colours from red (higher value) to blue (lower value) indicate the median log2FC across all detected genes in the corresponding pathway or GO term. To derive the average log2FC for a pathway/GO term, for each group we identified the corresponding DEGs, pooled all gene IDs from the three groups, retrieved the log2FC for each of these genes and established the median of log2FC per group. (E) Bar graph summarizing genes that displayed significantly different log2FC between groups (differential Mtb response). The y-axis indicates the cumulative DEG count while the x-axis indicates the pair-wise group contrast. Dark colour shading represents genes with higher response (e.g. Mtb effect in PLWH vs HC subjects resulted in a positive log2FC difference) while light shading depicts genes with lower response (negative log2FC difference). The union of all gene IDs identified across the three contrasts (n=401) results in 333 unique gene IDs. (F) Boxplot for three selected pathway/GO terms for genes with significant differential Mtb response between PLWH or PrEP against controls (n=302). For each term, the y-axis displays the mean log2FC and dots represent per-subject mean value for pooled DEG. Phenotypic groups are indicated below each box plot across the three terms. White dots represent mean log2FC for subjects who do not smoke cigarettes or used cannabis, while black dots depict smokers and/or cannabis users. To derive the subject mean log2FC, we used pooled DEG from the first two columns in panel E which were significant for interaction contrasts PLWH vs HC and PrEP vs HC, obtained subject-wise log2FC for each DEG and averaged the log2FC from all genes per subject. Results showed that group differences are not driven by outliers and that smoking, or cannabis consumption do not explain group differences. (G) Manhattan plot for enrichment test of pathways and GO-terms from merged DEG detected in the differential Mtb responses of PLWH vs HC and PrEP vs HC. The y-axis indicates the negative log10 false discovery rate (FDR) values while the tested terms are arranged along the x-axis. The horizontal dashed line represents the 10% FDR cut-off for significant pathways/GO-terms. Dots are sized as a function of the DEG number in a term and colors represent the database from which the terms were obtained. Identified terms represent innate immune processes and intracellular defense mechanisms.

Alveolar macrophages from PLWH and PrEP subjects display blunted mRNA transcriptomic changes in response to Mtb

Across the three phenotypic groups, we identified 1,626 differentially expressed genes (DEG) in response to Mtb (Data File S1). However, there were pronounced DEG count differences in their transcriptional responses (Fig. 1B). AM from HC subjects showed 1,434 DEG in response to Mtb as compared to only 679 and 255 DEG in PLWH and PrEP subjects, respectively (Fig. 1B and Fig. S1A). The magnitude of AM transcriptional response to Mtb also differed among the three groups. The smallest mean absolute log2 fold change (log2FC) in response to Mtb infection was observed among PrEP subjects while HC subjects displayed the strongest transcriptional response. (Fig. 1C and Fig. S1A). To test if Mtb effects in PLWH and PrEP subjects were weaker for all DEGs or only for specific gene sets, we assessed the correlations of Mtb-triggered logFC across the three groups. We observed that log2FC values were strongly correlated among all three groups, but consistently higher for corresponding genes from HC subjects (Fig. S1B) which suggested a general transcriptional impairment in cells from PLWH and PrEP subjects.

Pathway and GO-term enrichment analysis revealed that Mtb-induced transcriptional changes were enriched in biological functions representative of macrophage anti-mycobacterial host responses, such as interferon (IFN) and Toll-like receptor signalling pathways (Fig. 1D). However, compared to HC subjects, the number of significant GO-terms/pathways was substantially lower for the PLWH and PrEP groups (Table S1). Moreover, even among shared significant terms, the cumulative transcriptional response of genes in these terms was stronger in the HC group, which further emphasized the blunted anti-mycobacterial response by PLWH and PrEP subjects (Fig. 1D). We also tested which genes had a significantly different magnitude of response to Mtb among phenotype groups (Fig. 1E and Data File S1). We identified 333 genes significantly different between the groups (Data File S2) with the majority displaying lower up-regulation in response to Mtb for PLWH and PrEP groups relative to the HC group (Fig. S1, C and D and Table S1). These results provided a formal confirmation that the transcriptional response to Mtb challenge by AM of HC was more vigorous compared to the one mounted by PrEP and PLWH subjects.

To explore if the group differences were mainly driven by strong outliers, we derived per subject average log2FC across DEG in significant pathways/GO-terms. We found that AM from HC subjects displayed a more pronounced group transcriptional response to Mtb. For example, the mean log2FC for the term “I-kappaB kinase/NF-kappaB signaling” was 3.22-fold (p = 2.4×10⁻³) and 2.39-fold (p = 3.8×10⁻²) higher in HC as compared to PLWH and PrEP, respectively (Fig.1F). GO-terms and pathways tagged by the DEG from PLWH and PrEP vs HC, which were significantly more induced in AM from HC subjects, revealed a predominance of innate-immune and intracellular defense response terms (Fig. 1G). This observation alongside the exclusion of possible outlier-driven group differences supported the reduced ability of AM from PLWH and PrEP subjects to mount an effective anti-mycobacterial transcriptional response to challenge with Mtb.

Mtb does not induce chromatin remodelling in AM from PLWH and PrEP subjects

Next, we evaluated Mtb-triggered epigenetic changes of AM chromatin accessibility and H3K27 acetylation. Analysis of Differentially Open Chromatin (DOC) in response to Mtb indicated that AM from HC had significant increased accessibility (opening) in 8,389 regions and repression (closing) in 3,971 regions (Fig. 2A; Data File S3). In contrast, AM isolated from PLWH or PrEP subjects displayed near complete absence of significant chromatin remodeling in response to Mtb (Fig. 2, A and B). A broad-based epigenetic impairment was confirmed for AM from PLWH subjects by tracking regions with differently H3K27 acetylation (DAc), a mark of active chromatin (Fig. 3A; Data File S4). Pathway and GO-term enrichment analyses of genes annotated to DOC and DAc in HC identified typical viral and bacterial response terms (Fig. 2C; Fig. 3B). As in the transcriptomics analyses, the epigenetic response to Mtb by AM from HC strongly implicated host immune response to pathogens and immune-related signalling pathways, which was mainly driven by opening DOCs and regions with increased acetylated H3K27 (Fig. 2D; Table S2). Examples of genes with DOCs and DAc at the gene transcription start site (TSS) were CXCL10, IFI44L, APOBEC3A and MX1 (Fig. S3 and S4).

• Download figure
• Open in new tab

Fig. 2 Alveolar macrophage chromatin remodeling in response to Mtb.

(A) Significance of chromatin changes in response to Mtb is plotted against the magnitude of the change. Significantly opened or closed chromatin regions at a false discovery rate (FDR) of 5% after Mtb challenge are represented by dots marked in red or blue, respectively. The counts of differential open/closed chromatin (DOC) regions are given at the top left and right corners. Dashed lines indicate the 5% FDR on the y-axis and the minimum log2FC of 0.2 on the x-axis. In healthy controls (HC), 23.3% of the tested regions displayed significant chromatin remodelling in response to Mtb, while subjects on pre-exposure prophylaxis (PrEP), and persons living with HIV receiving antiretroviral treatments (PLWH) lacked significant chromatin changes. (B) Density of the absolute log2FC effect for 12,360 DOC regions observed in the HC response to Mtb is shown for HC, PrEP, and PLWH subjects in blue, green, and red, respectively. The mean log2FC for each group is highlighted by dashed horizontal lines. Both PrEP and PLWH groups show strongly impaired chromatin remodelling in response to Mtb compared to HC. (C) Pathway/GO-term enrichment analysis of genes assigned to DOC regions for the HC response to Mtb. GO biological process (BP), KEGG, and Reactome pathways to which at least five genes had been assigned are plotted against the negative log10 FDR. In total, 724 out of the 8066 tested terms are significant. Enrichment analysis indicated that in AM of HC subjects, Mtb challenge promoted chromatin remodelling in regions assigned to genes belonging to IFN and TNF signalling pathways as well as response to viral and bacterial pathogens. (D) Rank plot of selected GO BP, KEGG, and Reactome pathways. DOC regions in AM of HC subjects tagged genes assigned to “I−kappaB kinase/NF−kappaB signaling”, “TNF signalling”, and “Interferon signalling” pathways. DOC are plotted according to their log2FC on the x-axis and ranked according to the magnitude of change on the y-axis. Genes corresponding to opening or closing DOC regions are shown in red and blue, respectively. Horizontal doted lines indicate the log2FC cut off of -0.2 and 0.2. In the selected pathways, the majority of the chromatin changes observed indicated increased accessibility in response to Mtb.

• Download figure
• Open in new tab

Fig. 3 Histone H3K27 acetylation changes in response to Mtb.

a, Volcano plots for H3K27 acetylation response to Mtb. Chromatin regions significantly acetylated at a false discovery rate (FDR) < 5% after Mtb challenge are marked in red and blue, respectively. The numbers of acetylated regions are given at the top. Significant changes in H3K27 acetylation, indicating active enhancers, were observed mostly for the healthy control (HC) group, which encompassed five subjects compared to 11 persons living with HIV on antiretroviral therapy (PLWH). While many peaks for H3K27ac did not pass FDR 5%, the magnitude of log2FC was higher in HC compared to PLWH. b, Enrichment analysis of genes assigned to regions with increased H3K27 acetylation for the HC response to Mtb. GO biological processes (BP), KEGG, and Reactome pathways with at least five assigned genes are plotted against the negative log10 FDR. A total of 167 out of 906 tested terms were significant. As observed for the chromatin changes, enrichment analysis indicated that in AM of HC subjects, Mtb challenge promoted higher acetylation in regions assigned to genes belonging to IFN and TNF signalling pathways as well as response to viral and bacterial pathogens.

Covariates do not explain Mtb-triggered transcriptional and epigenetic differences among groups

The three groups differed with respect to several covariates (Table 1). Among the HC group, the proportion of females was higher, and the use of recreational drugs was lower compared to the PLWH and PrEP groups. Subjects in the PrEP group on average were younger compared to the two other groups. The duration of ART was longer for PLWH (mean = 12.12 years, standard deviation (SD) = 7.64) than PrEP subjects (mean = 3 years, SD = 1,47). All PrEP subjects received Emtricitabine and Tenofovir disoproxil fumarate (TDF) while among PLWH there was a wider spread of ART drug regimens. However, only 4 out of 20 subjects had not received TDF at the time of enrolment (Table 1). Most of the inter-subject variability was accounted for by the analytical approach employed. Still, given the strong heterogeneity of study subjects we investigated the impact of smoking, recreational drug use and age on Mtb-triggered chromatin and transcriptional responses. We found that recreational drug use did not explain the reduced transcriptional activity or absence of chromatin changes in response to Mtb among PLWH and PrEP subjects (Fig. 1F and Fig. 4A). Age and duration of ART/PrEP were also unlikely confounders since PLWH and PrEP subjects differed significantly in age and ART duration, yet both groups showed similar reduced transcriptional activity and lack of chromatin changes relative to HC subjects following challenge with Mtb (Fig. 4, B and C). Moreover, when stratifying HC and PLWH by age, HC subjects of 40 years and older showed pronounced Mtb-triggered chromatin changes that were not observed for PLWH subjects (Fig. 4D). In addition to the covariates listed in Table 1, we also considered possible group-dependent differences in phagocytic activity. Since the yields of AM were generally not sufficient to conduct independent assessment of phagocytic activity, we used presence of Mtb reads in ATAC-seq libraries as proxy and failed to detect any significant group-dependent differences (Fig. 4E).

• Download figure
• Open in new tab

Fig. 4 Impact of covariates on Mtb-triggered transcriptional and epigenetic changes in AM.

(A) Volcano plots for chromatin accessibility changes by AM in response to Mtb after removing subjects who smoked cigarettes, cannabis, or both. Included in this analysis were 10 HC, 9 PLWH, and 6 PrEP. Significantly opened or closed chromatin regions at a false discovery rate (FDR) of 5% after Mtb challenge are represented by dots marked in red or blue, respectively. The counts of differential open/closed chromatin (DOC) regions are given at the top left and right corners. Dashed lines indicate the 5% FDR on the y-axis and the minimum log2FC of 0.2 on the x-axis. (B) Boxplot for three selected pathway/GO terms for genes with significant differential Mtb response between PLWH or PrEP against HC (n=302). For each term, the y-axis displays the mean log2FC and dots represent per-subject mean value for pooled DEG from the three groups. White dots represent mean log2FC for subjects younger than 40 years at enrolment, while black dots depict subjects ≥ 40 years of age. To derive the subject mean log2FC, we used DEG from Fig. 1e which were part of significant pathway/GO terms, pooled the genes for each group, obtained the log2FC per subject for each DEG and averaged the log2FC from all genes per subject. (C) Boxplot showing the distribution of age at enrolment across the three studied groups. Subjects of the PLWH group were significantly older than the subjects of the healthy controls and PrEP groups. (D) Volcano plots for chromatin accessibility changes by AM in response to Mtb stratified by age for subjects of the HC and PLWH groups. Subjects tested were ≥ 40 years at enrolment. In this age class the HC, PLWH and PrEP groups encompassed 9, 15 and 3 subjects, respectively. Due to low subject count PrEP subjects were excluded from this stratified analysis. Significantly opened or closed chromatin regions at a false discovery rate (FDR) of 5% after Mtb challenge are represented by dots marked in red or blue, respectively. The counts of differential open/closed chromatin (DOC) regions are given at the top left and right corners. Dashed lines indicate the 5% FDR on the y-axis and the minimum log2FC of 0.2 on the x-axis. (E) Boxplot of distribution for Mtb-aligned reads in ATAC-seq experiments. The y-axis indicates the proportion of total unique paired-end reads aligned to the Mtb H37rv genome relative to the Human genome hg38. Unique alignments to Mtb genome quantify DNA traces of phagocyted/lysed mycobacteria. In the absence of Mtb challenge, little to no alignment was observed to Mtb genome. In the challenge condition, the three groups (HC, PLWH, and PrEP) presented a range of proportions of Mtb unique reads. No statistical significance was observed between the proportion of Mtb reads between the phenotype groups.

Chromatin accessibility and binding of transcription factors

To link Mtb-triggered epigenetic changes in AM from HC subjects with increased gene expression, we evaluated the enrichment of transcription factors (TFs) motifs in Mtb-triggered DOC and DAc regions and found an enrichment of IFN-regulatory factor (IRF) motifs in those regions (Table 2). Next, we compared the average difference in TF footprint depth between Mtb challenged and non-challenge AM to estimate TF activity. Of the 682 TFs obtained from the JASPAR catalog, 21 had significantly higher activity (increased footprint depth) in DOC regions after Mtb challenge (Fig. 5A). Among these significant TFs were IRF9 (part of the anti-mycobacterial host response) and ZNF684 (a TB biomarker). Of note, Mtb challenge had triggered a higher fold change for IRF9 and ZNF684 RNA expression in AM from HC subjects compared to AM from PrEP or PLWH subjects (Fig. 5B). Moreover, IRF9 and ZNF684 footprints responsive to Mtb challenge were enriched in the transcription starting site (TSS) of genes that are part of the IFN signalling pathway (p = 8.9 x 10⁻⁶ and p = 0.03, respectively) (Fig. 5C). Unexpectedly, we also observed an enrichment of IRF9 active footprints in the TSS of genes of the TNF signalling pathway (p = 0.01) (Fig. 5C). ZNF684 active footprints were observed in the core promoter of TNF after Mtb challenge providing a direct link to a major pro-inflammatory cytokine (Fig. 5D). Moreover, in support of a coordinated effect of these two TFs, an active footprint for ZNF684 was observed in the promoter region of IRF9 (Fig. 5D).

• Download figure
• Open in new tab

Fig. 5. Transcription factor (TF) footprint activity in differential open/closed chromatin (DOC) regions.

(A) Activity scores denoting differences in depth of TF footprints between healthy controls (HC) and HC + Mtb AM are plotted on the y-axis for 682 TFs from the JASPAR catalog showing a minimum of 50 footprints in DOC regions. A positive activity score indicates that footprints for a TF are more pronounced in HC after Mtb infection. Twenty-one TF with significant activity score at FDR < 5% are displayed as red dots. (B) mRNA expression levels for the transcription factors IRF9 and ZNF684 are shown as violin plots. The residual log2 copy per million (CPM) reads, after removing inter-individual variability, is plotted for each group i.e., healthy controls (HC), subjects on pre-exposure prophylaxis (PrEP), and persons living with HIV receiving antiretroviral treatments (PLWH) and condition (Mtb-challenged or not). (C) Percentage of genes in the TNF or IFN pathways that have IRF9 or ZNF684 footprints in their transcription starting site (TSS; ± 5kb from the promoter) was compared to all genes with TSS in DOC regions. As expected, we observed that IRF9 footprints are more frequent in genes assigned to Reactome’s IFN signalling pathway compared to random TSS in DOCs. Surprisingly, IRF9 footprints were also more frequent than by chance in KEGG’s TNF pathway. (D) ZNF684 and IRF9 footprints in selected IFN and TNF genes. The mean normalized fragment pileup is plotted on the y-axis as red and blue lines for non-infected and Mtb-challenged AM from HC, respectively. The standard deviation of the means is shown by shades for each condition. The location of the motif assigned to the footprint is shown by a black bar at the bottom. We detected a pronounced footprint for ZNF684 in the core promoter of the TNF gene and also in the promoter of IRF9. Motifs for IRF9 were detected in footprints located close to the promoter of the NOD2 and ISG15 genes.

To further assess the link between Mtb-dependent TF binding and RNA expression, we evaluated if the 135 expressed genes with active IRF9 footprint at the TSS were differentially expressed in HC following Mtb challenge. We observed that 65 out of the 135 genes were DEGs (57 up-regulated and 8 down-regulated) following Mtb challenge (p = 7.9 x 10⁻²³, Table S3 and Data File 5). AM from PLWH and PrEP had lower logFC in response to Mtb for all 65 DEGs (Data File 5). Of the 57 up-regulated DEGs with IRF9 TSS footprint, 17 were part of the IFN or TNF pathways, including NOD2 and ISG15 (Fig. 5D). Collectively, these data showed that the epigenetic response observed in AM from HC subjects promoted the binding of key transcriptions factors required for IFN and TNF signalling. The locked state of chromatin in AM from PLWH and PrEP subjects was therefore likely to account, at least in part, for the reduced transcriptional activity of those cells in response to Mtb.

Subjects

The study was approved by the Research Ethics Board of the McGill University Health Center (MUHC) (MP-CUSM-15-406). Participants were recruited at the McGill University Health Centre (MUHC, Montreal, Canada). HIV negative participants represented a control group (n=23) and a high-risk group on preventive therapy (PrEP, n=14). HIV positive participants (PLWH, n=27) were receiving antiretroviral therapy (ART) to suppress peripheral viral load for at least 3 years and without any active respiratory symptoms or infections. Exclusion criteria included past history of TB, pregnancy, chronic cardiovascular or active pulmonary infection and immunosuppressive medication. All participants were tested with QuantiFERON Gold, Cat#0594-2010/Cat#T0590-0301 (Qiagen, Germany) and one subject each of the HC and PrEP groups tested positive. Since exclusion of these two subjects had no significant impact on reported results, we left them included in the present presentation of the data. Information on participant characteristics and antiretroviral drugs was extracted from the clinical database and participant electronic records. Of these enrolled subjects, we excluded 7 HC and 7 PLWH from the analyses due to low yield of BAL cells, or inability to pass quality control for sample or library preparation.

Bronchoalveolar lavage and sample processing

Bronchoscopies were performed at the Centre for Innovative Medicine of the Research Institute of the MUHC following the guidelines of the American Thoracic Society (30). Briefly, standard flexible bronchoscopy of the middle lobe was conducted under local anesthesia with lidocaine with additional conscious sedation with intravenous midazolam. A total of 200 ml sterile saline instilled in 50-ml volumes resulted in an average return of 115 ml (range of 50-150 ml) of bronchoalveolar lavage (BAL) fluid. BAL fluid in non-adherent tubes was stored on ice for a maximum of 30 min prior to the isolation of AM. All BAL fluids were strained through a 100μm filter to remove any mucus and cell clumps. Cytology assessment of BAL samples indicated that AM represented >90% of all recovered mononuclear cells.

Alveolar macrophage preparation

AM were isolated by adhesion. Collected BAL cells were spun for 10 minutes at 300g, 10°C, and washed twice in RPMI-1640 with L-glutamine (Wisent, Canada), containing 2% human serum (heat-inactivated AB+ off the clot, Wisent, Canada), 1% penicillin/streptomycin (Gibco, USA), 10mM HEPES (Gibco, USA), 1% Non-essential amino acids (Gibco, USA) and 1.75ng/ml Amphotericin B (Wisent, Canada). After the cell count, 1.5×10⁶ AM were seeded per well of a 6-well plate (Falcon, USA) in RPMI-1640 supplemented with 10% human serum and incubated for 1.5h-2h at 37°C, 5% CO2, and 95% relative humidity. Non-adherent cells were removed by thorough washings with pre-warmed (37°C) RPMI-1640 supplemented with 2% human serum. Washed adherent AM were then infected with freshly prepared Mtb. An estimate for the number of adherent AM was derived by subtracting the number of detached from the number of seeded cells.

Mycobacterial cultures and macrophage challenge in vitro

Virulent Mtb strain H37Rv was grown in a liquid culture of Middlebrook 7H9 medium (Difco, Detroit, USA, #295939) containing 0.2% glycerol (Fisher), 0.05% Tween-80 (Sigma-Aldrich) and 10% albumin-dextrose-catalase (ACD, Fisher) at 37°C in rolling incubators. Bacteria were grown to log phase, as determined by an optical density of 0.3 to 0.8 at 600nm, prior to inoculum preparation. Further, bacterial cultures were spun for 15 minutes at 3700 rpm, resuspended in RPMI-1640 and dislodged with a 22G needle. Cell suspensions were filtered through 5μm filters (Millipore, USA) to ensure single mycobacteria suspensions for the challenge experiments. Bacterial counts of inocula were done using disposable Neubauer hemocytometer (C-Chip, INCYTO, South Korea). Bacterial loads were confirmed by colony-forming unit (CFU) counts by plating serial dilution of inoculum in 7H9 growth medium on Middlebrook 7H10 agar (Becton Dickinson, Franklin Lakes, NJ, USA, #254520) plates containing 0.5% glycerol (Fisher) and 10% oleic acid-albumin-dextrose-catalase (OADC, Becton Dickinson, Franklin Lakes, USA, #211886). Colonies were counted 3-4 weeks post-plating. AM were infected at a multiplicity of infection (MOI) of 5:1 for 18-20hrs 37°C, 5% CO2, and 95% relative humidity.

Preparation of RNASeq and ATACSeq libraries

For downstream RNA and ATAC experiments AM were rinsed with warm PBS to remove damaged cells and debris before further processing. For RNA, cells were lysed with Trizol and kept at -80°C until extraction. For RNA isolation the miRNeasy kit (Qiagen, Germany) was used and RNA integrity (RIN) was assessed with the Agilent 2100 Bioanalyzer (Agilent Technologies, Waldbronn, Germany). Only samples with RIN >8 were selected for library preparation using TruSeq RNA Library Preparation Kit v2, Set A (Illumina, USA).

For chromatin accessibility using the ATACseq method (31, 32); AM were lifted with pre-warmed CellStripper solution (Corning, USA) and gentle pipetting. Viability was assessed with trypan blue staining, and only samples with viability > 85% were used. We followed the protocol by Buenrostro et al. (32). Briefly, 50,000 cells were used for generating challenged and non-challenged libraries. AM were permeabilized with lysis buffer containing 0.05% IGEPAL. Nuclei were incubated with Tn5 transposase (Illumina, USA) and further ATAC-seq libraries were amplified for a total of 5-10 PCR cycles. Barcoded and multiplexed RNA or ATAC-seq libraries were sequenced on an Illumina HiSeq4000 sequencer at Genome Quebec, Montreal, Canada.

ChIP-seq library preparation and sequencing

Samples from challenged and non-challenged AM were crosslinked with 1% w/v formaldehyde for 10 min at RT and immediately quenched for 5 min with 125 mM Glycine at RT. The formaldehyde fixed samples were then sonicated to 150-500 bp using a S220 (Covaris) and then ChIP-DNA prepared using a manual chromatin immunoprecipitation method with an Antibody-Antigen incubation of 18 hrs., followed by bead incubation for 135 minutes, and 6 x 5-minutes washing steps.

Approximately 0.5 million cells were used for each ChIP and ∼50,000 cells for the input. We used H3K27ac antibody (Abcam, Ab4729, GR3211959-1). ChIP and Input libraries were prepared using the MicroPlex Library Preparation Kit (C05010012, Diagenode), with alterations including PCR enrichment (7 to 20 cycles) prior to size selection and use of Ampure beads for size selection (250-350 bp). Libraries for ChIPSeq were sequenced on an Illumina NovaSeq 6000 S2 SR100.

Data processing and statistical analysis

Detailed description for raw sequences processing to statistical analyses are provided in Supplementary Information. Briefly, Baseline and challenged transcriptome and epigenomic landscapes of AMs were determined following published protocols (31, 33). Gene expression readouts were generated using alignments and annotated set of transcripts as input to Salmon (34). Differences in chromatin accessibility and H3K27ac levels were assessed with a count-based peak quantification as shown in (35). Quantification files from each approach were merged as a single matrix, genes and peaks with low counts were filtered out and libraries were normalized, scaled and log2 transformed using edgeR and Limma-voom (36, 37).

For the statistical analysis, we defined linear models with blocking design on subject ID to account for inter-individual differences and modelled the challenge effect nested within phenotypic groups. For all experiments we tested the challenge effect for HC, PrEP and PLWH subjects, while for RNAseq we also contrasted the phenotype groups against each other (interaction test). ATACseq and ChIPseq p-values were adjusted with Storey’s false discovery rate (FDR) and peaks considered significant if FDR < 0.05 and absolute log fold-change (logFC) > 0.2. These were termed differentially open chromatin (DOC) for ATACseq and differentially acetylated (DAc) when referring to ChIPseq. For RNA-seq we employed a two-stage FDR procedure and we considered genes to be differentially expressed genes (DEG), if absolute logFC ≥ 0.2 and stageR p-val ≤ 0.05 (38).

We performed pathway and gene ontology enrichment analysis using three independent databases (KEGG, Reactome and AmiGO). Enrichment results were merged and the Benjamini-Hochberg FDR method was applied to raw p-values. We tested TF motif enrichment in DOCs and DAc regions using HOMER (39). Mtb-triggered differences in TF footprint activity over DOC regions were estimated with HINT-ATAC (40). TF motifs with FDR < 0.05 and activity score > 0.005 were considered significant.