DNA Methylation Proles of Immune Cells From Tuberculosis-Exposed Individuals Overlap With BCG-Induced Epigenetic Changes

The mechanism of protection of the only approved tuberculosis (TB) vaccine, Bacillus Calmette Guérin (BCG) is poorly understood. In recent years, epigenetic modications induced by BCG have been demonstrated to reect a state of trained immunity. The concept of trained immunity is now explored as a potential prevention strategy for a variety of infections. Studies on human TB immunity are dominated by those using peripheral blood as surrogate markers for immunity. Here, we instead studied the lung compartment by obtaining induced sputum from subjects included in a TB contact tracing. CD3- and HLA-DR-positive cells were isolated from the collected sputum and DNA methylome analyses performed. Unsupervised cluster analysis revealed that DNA methylomes of cells from TB-exposed individuals and controls appeared as separate clusters, and the numerous genes that were differentially methylated were functionally connected. The enriched pathways were strongly correlated to previously reported epigenetic changes and trained immunity in immune cells exposed to the BCG vaccine in human and animal studies. We further demonstrated that similar pathways were epigenetically modied in human macrophages trained with BCG in vitro. Altogether, our study demonstrates that similar epigenetic changes are induced by M. tuberculosis and BCG.


Introduction
Tuberculosis (TB) is a pulmonary infection of pandemic rank and an expansion of the current toolkit for diagnosis, prevention and treatment is critical for reaching the United Nations' Sustainable Development Goals for 2030 of ending the TB epidemic 1 . TB is caused by Mycobacterium tuberculosis, which transmit via aerosols and target alveolar macrophages in exposed individuals 2 . In susceptible hosts, the bacteria replicate inside the macrophages and use them as trojan horses in order to disseminate in the tissues 3 .
Bacillus Calmette Guérin, a non-virulent derivative of Mycobacterium bovis, has been used for almost a century as a vaccine against TB, with variable e cacy. Numerous studies have failed to provide correlates of protection, leaving the vaccine mechanism elusive. In recent years, the concept of trained immunity has evolved as an epigenetically encoded immune memory that can be triggered by a variety of stimuli and is re ected in a reprogrammed immune state characterized by a higher magnitude of response to subsequent pathogen challenges. The discovery of epigenetically regulated antimicrobial defense mechanisms goes beyond the classical understanding of immune defense and opens up a new eld of research. Along this line, we have demonstrated that administration of the BCG vaccine to healthy subjects induced profound epigenetic alterations in immune cells, which correlated with enhanced antimycobacterial activity in macrophages isolated from the vaccinees 4 . The changes were re ected in the DNA methylome, with the strongest response being recorded within weeks after vaccination 4 . Our observation that BCG induces alterations of the DNA methylome of immune cells has later been con rmed by others 5,6 . Since BCG vaccination re ects an in vivo interaction between immune cells and viable mycobacteria, we here hypothesized that natural exposure to M. tuberculosis would induce similar changes not only in TB patients, but also in individuals who have been exposed to TB. Analyses of DNA methylomes of immune cells isolated from lungs and peripheral blood allowed us to identify distinct DNA methylation (DNAm) signatures in TB-exposed individuals. Pathway analyses revealed strong overlaps with previous studies on BCG-induced epigenetic signatures that could be correlated with protection against M. tuberculosis. We therefore set up an experiment to scrutinize the anti-mycobacterial mechanisms correlating with BCG-induced epigenetic reprogramming in human primary macrophages.
BCG-trained macrophages were subjected to DNA methylome analysis and in parallel monitored in a live cell imaging system during exposure to M. tuberculosis. Again, the DNAm change in vitro highly overlapped with the signatures of the TB-exposed individuals and the previously identi ed BCG-induced epigenetic changes.

Study design
To determine epigenetic changes in the immune cells in TB-exposed individuals, we recruited subjects enrolled in a TB contact tracing at Linköping University Hospital, Sweden, according to standard contact tracing routines (household contacts and persons with >8 hours of interaction with the index case). First, we included one index case and three contact tracing subjects who were ethnicity-matched and not relatives ( Table 1). The patient with active tuberculosis (referred to as 'index case') was diagnosed with drug-sensitive pulmonary TB and had completed two out of six months of standard treatment at the time of sample collection. Age-and ethnicity-matched individuals were included as controls (Table 1). In the second part of the study, we included three subjects (one TB patient and two subjects in other contact tracings) who were not ethnicity-matched (Table 1). All included subjects except one (a TB contact) had been BCG-vaccinated at least two years before participation (Table 1). Interferon-Gamma Release Assay (IGRA) status was determined and among the exposed individuals, two were positive (including the index case) and among the controls, one individual (C2) was classi ed as 'borderline'-positive 7 (Table 1). From induced sputum, we used an established protocol for the separate isolation of HLA-DR-positive and CD3 positive cells 8 , whereas the PBMC fraction extracted from blood was kept as a mixed population. For the in vitro experiment with BCG and M. tuberculosis exposure, human primary monocytes were collected from blood samples from healthy blood donors and differentiated into macrophages (Fig. 1).
DNA methylome data from TB-exposed individuals form a separate cluster DNA isolation from the studied cell populations was followed by global DNAm analysis using the Illumina 450K protocol. After curation of the data 9 , the datasets were subjected to unsupervised hierarchical cluster analysis based on DNA CpG methylation β-values. This unsupervised approach accurately clustered the participants into TB-exposed and controls based on the DNA methylome data derived from both HLA-DR-and CD3-positive cell populations ( Fig. 2a-b). On the other hand, in the PBMCderived dataset, the TB index case appeared outside the clusters and two of the controls clustered with the other exposed individuals, one of them ("Con_2") being the individual identi ed as IGRA borderlinepositive (Fig. 2c).
Next, we identi ed the differentially methylated CpG sites (DMCs) and differentially methylated genes (DMGs) by comparing the TB-exposed and control groups for each cell population. To lter out the most signi cantly altered DMGs in the dataset, the stringency criteria of log 2 >|0.3| fold increased or decreased β-values and Benjamini-Hochberg (BH)-corrected p-value <0.05 (HLA-DR), <0.1 (CD3) and <0.2 (PBMC) were applied. The results are depicted as volcano plots, which show that DNA methylomes of TB-exposed most strongly differ in the HLA-DR cells as compared to control subjects, followed by the CD3 population, whereas PBMC datasets reveled fewer DMGs ( Fig. 2d-f, Suppl. Table 1). To highlight the locus position of the DMGs, chromosome maps were constructed (Suppl. Fig. S1a-c). Using the same stringency criteria as for the HLA-DR analysis, we tested whether DMGs would emerge when the datasets were arranged in other possible groups as derived from the demographics (>/< median age, sex, IGRA status). Neither age nor IGRA status generated any signi cant DMGs with these settings, and gender rendered three DMGs (X and Y chromosomes were removed in the initial ltering prior the analysis) in the HLA-DR positive and in the CD3 positive cells (Suppl. Table 1).

Functional enrichment analysis reveals common and unique interactomes in the datasets
Using the PANTHER Database, we investigated whether the identi ed DMGs were enriched in known pathways ( Fig. 3a-c). The analysis revealed enrichment in pathways with relevance for TB infection, including hypoxia-inducible factor (HIF)1-α activation, Vitamin D metabolism and p38, Wnt, Notch, interleukin, chemokine, and cytokine signaling pathways 10−17 . Common pathways shared between at least two of the cell populations included B cell activation, glycolysis, angiotensin II signaling, and cholecystokinin signaling. Notably, several pathways named after their known functions in the nervous system were enriched in the studied cell populations, including pathways involved in axon guidance and adrenaline, acetylcholine and glutamate signaling. In the PBMC population but not in the lung cell populations, the interferon-γ signaling pathway was identi ed among the enriched pathways.
Comparisons across cell populations and species reveals the existence of a common DNA methylomebased biosignature in mycobacteria-exposed immune cells Given the fact that the interaction between mycobacteria and eukaryotes is evolutionary ancient, we predicted that highly conserved pathways exist that are common among the studied cell populations. Comparing the identi ed DMGs from the HLA-DR, CD3 and PBMCs in a Venn analysis, we discovered 185 common DMGs (Fig. 4a). We expanded the Venn analyses to include data from our previous work on BCG vaccine-induced DMGs that correlated with enhanced mycobacterial control 4 , as natural exposure to TB and BCG vaccination both represent in vivo encounters between mycobacteria and host immune cells. Even though the routes of mycobacterial exposure differ profoundly in these settings, a set of 151 DMGs could be identi ed as overlapping between our previous BCG study and all cell populations studied here (Fig. 4b), suggesting that a highly conserved epigenetic response to mycobacterial challenge exists.
In 2018, Hasso-Agopsowicz et al. described alterations in DNAm patterns in PBMCs from BCG-vaccinated individuals, with concomitant enrichment in many immune-related pathways 5 . In order to compare that study with ours, we performed PANTHER analysis with the 185 common DMGs and matched the identi ed enriched pathways, revealing that 75% of those pathways were the same as in the present study ( Fig. 5a), further corroborating the relationship of the altered DNAm patterns induced through TB exposure and BCG-induced changes. In a study by Kaufmann et al., BCG-induced alterations of the epigenome in mice was correlated with protection against M. tuberculosis infection 18 . In order to translate our human DNA methylome signature to the signature identi ed in the mouse study, we searched for pathway overlaps between the two using Gene Ontology (GO) enrichment analysis (Suppl.  Figure 5b demonstrates that for our PBMC data, the GO terms "biological processes" overlapped 100% with the mouse study (same cell population) and to 31% and 65% for HLA-DR-and CD3positive cells, respectively. In 2014, Saeed et al 19 demonstrated the induction of trained immunity pathways by another immune-training agent, β-glucan. We assessed possible pathway overlap with that study and although there were fewer overlaps as compared to the BCG-induced pathways described above. Again, the strongest correlation was found in the PBMC fraction, in this case in the GO terms "cellular components" (Fig. 5c).
To validate how well the 284 CpG sites corresponding to the 185 overlapping DMGs performed in an unsupervised cluster analysis, we included one additional TB patient and two contacts, and collected HLA-DR-positive cells from induced sputum, as the DNA methylome data from this cell type was clearly outperforming the others with respect to accurate separation of the groups. Figure 6 shows a k meansbased dendrogram with a heatmap of the β values of the 284 CpG sites from both previous and new subjects' DNAm data, revealing a distinct separation of the subjects in accordance with TB exposure.
In vitro BCG training of macrophages induced DNAm changes corresponding to exposure to TB To investigate whether the BCG-induced epigenetic changes can be mimicked in vitro, we set up a BCG training experiment with macrophages isolated from donor blood that were trained with BCG and subsequently infected with the virulent M. tuberculosis strain H37Rv (expressing GFP) and monitored during the course of infection. In a subset of donors, BCG-trained cells displayed an increased capacity to kill M. tuberculosis ( Fig. 7a-b) while cell viability was retained (Fig. 7c-d). DNAm analyses was performed on DNA isolated from these donors' macrophages 24 hours after BCG training. We identi ed 7471 DMGs with the stringency criteria of DNAm difference > 30% and BH corrected p-value < 0.01. A PANTHER pathway analysis based on the identi ed DMGs showed signi cant enrichment in the Wnt signaling pathway, cadherin signaling pathways and angiogenesis, overlapping with the over-represented pathways in the HLA-DR+ and CD3+ cells from the TB-exposed individuals (Fig. 7e).

Discussion
In this study, we present data suggesting that exposure to TB generates a distinct DNAm signature in pulmonary immune cells. The signature was found not only in those with active or latent TB infection, but also in individuals who are exposed but IGRA-negative. The nding that healthy, TB-exposed individuals also carry the signature opens the possibility that the epigenetic alterations re ect a host-bene cial reprogramming of the immune mechanisms rather than being induced by M. tuberculosis as a step to evade the immune defense. This notion is supported by our observation that the DMGs identi ed in the present study strongly overlapped with the epigenetic alterations identi ed in the in vitro BCG-trained macrophages, and the previously reported DNAm changes induced during BCG vaccination, which correlated with increased anti-mycobacterial capacity of macrophages 4 . In addition, we demonstrate that the GO data derived from our dataset display a strong overlap with data from a study on protective BCG vaccination in mice 6 . BCG vaccination has convincingly been shown to induce heterologous immunity protecting against childhood mortality from other causes than TB 20,21 . Based on our nding that natural TB exposure and BCG vaccination trigger similar epigenetic changes we propose the hypothesis that a "bene cial exposure" to TB exists, which protects against other infections through heterologous immunity. Along the same line, it has been shown that a substantial fraction of individuals exposed to TB can be de ned as 'early clearers', since they remain tuberculin skin test or IGRA negative 22 , suggesting effective eradication of the infection 22 . Identifying these early clearers and understanding the biology behind their resistance to TB infection could move the eld forward towards novel strategies of TB prevention.
In concordance with the fact that macrophages constitute the main niche for mycobacterial replication, the strongest enrichment of DNAm changes was observed in the HLA-DR-positive cell population, which is dominated by alveolar macrophages. The pathways identi ed to be enriched in the HLA-DR-positive population have been described in the context of trained immunity, BCG exposure and TB. For example, activation of Hypoxia-Inducible Factor 1 α and glycolysis pathways (P00030 and P00024, respectively) are hallmarks of macrophages that have undergone the epigenetic changes re ective of trained immunity (reviewed in 23,24 ), which is induced in myeloid cells upon BCG-stimulation 22,23 . VEGF-release (P00056) by macrophages has been shown to recruit immune cells during granuloma formation 25 . Further, vitamin D has been shown to strengthen the anti-mycobacterial activity of macrophages 11,26 , and upregulation of components of the vitamin D pathway is linked to the production of anti-microbial peptides 12 , providing a possible effector mechanism for mycobacterial control. Recent literature on immune regulation through T cell-derived acetylcholine 27,28 attributes relevance to the acetylcholine receptor pathway identi ed among the HLA-DR pathways.
Although macrophages and lymphocytes are not generally viewed as having many similarities, we found 34 of the identi ed pathways to overlap between HLA-DR-and CD3-positive cells. In data derived from the CD3 and PBMC populations, both of which represent lymphocytes, overlaps were identi ed for glycolysis, glutamate receptor and angiotensin II pathways. Interestingly, a metabolic shift towards increased glycolysis, representative of the Warburg effect, has been strongly associated with trained immunity 23 .
However, the literature is dominated by the view that this event takes place in trained myeloid cells, while we identi ed this circuit in CD3 cells (T lymphocytes) and not in the HLA-DR cells (mainly macrophages). The glutamate receptor is widely expressed on immune cells and has been described as having an important regulatory role in T cells, which can also produce and release glutamate 29 . The role for angiotensin II pathway in TB remains elusive, while Angiotensin II Converting Enzyme 2 has been in the spotlight due to fact that the SARS-CoV2 virus utilizes it as a receptor for entry into host cells 30 . In the PBMC population, which over all showed a weaker epigenetic response, we found the interferon-γ signaling pathway, which has a well-established role in anti-mycobacterial defense (reviewed in 31  Taken together, we present data supportive of mycobacteria exposure-induced DNAm changes that correlate with previous ndings from studies on BCG vaccination including TB protection and trained immunity.

Study design and ethics
Patients with pulmonary TB, healthy participants with a history of TB-exposure and healthy controls, with an age ranging from 18 to 53 years, were enrolled at Linköping University Hospital and Linköping University, respectively. Included subjects (Table 1) donated peripheral blood and induced sputum samples 8 following oral and written informed consent (ethical approval obtained from the Regional Ethics Review Board (EPN) in Linköping, no. 2016/237-31). The study protocol included questionnaires on respiratory and overall health, the evaluation of IGRA-status and sputum samples for DNA extraction. The subjects' samples and questionnaires were not linked to any personal information at any stage of the study. The sputum sample collections were performed in accordance with guidelines from the Department of Respiratory Medicine at Linköping University Hospital. The IGRA samples (QuantiFERON-TB Gold) were collected and analyzed by medical personnel according to the guidelines at the Department of Clinical Microbiology at Linköping University Hospital. For in vitro experiments, deidenti ed buffy coats were purchased from the blood facility at Linköping University Hospital. The buffy coats were obtained from healthy blood donors, who gave written consent of research use for the blood products.

Induced sputum and pulmonary immune cell isolation
Induced sputum is a well-tolerated, non-invasive method to collect cells from the surface of the bronchial airways after inhalation of a hypertonic saline solution. The procedure of sputum induction takes approximately 30 minutes and is both cost effective and safe with minimal clinical risks 32 . Sputum specimens were collected as described by Alexis et al 33 , with the following modi cations: premedication with an adrenergic β2-agonist, salbutamol (Ventoline, 1ml 1mg/ml) was administrated before the inhalation of hypertonic saline, using a nebulizer (eFlow, PARI). The subsequent steps of sputum processing were adopted from Alexis et al. (2005) 34 and Sikkeland et al 35 . The HLA-DR and CD3-positive cells were isolated using superparamagnetic beads coupled with anti-human CD3 and Pan Mouse IgG antibodies and HLA-DR/human MHC class II antibodies (Invitrogen Dynabeads, ThermoFisher, cat no. 11041 and 14-9956-82, respectively). An initial positive selection was done with CD3 beads followed by a positive HLA-DR selection. Bead-coating and cell isolation were performed according to manufacturer's protocol.

PBMC isolation
PBMCs were isolated from whole blood (TB-exposed and controls) or from leukocyte rich fractions of blood obtained from healthy volunteers (Linköping University Hospital blood bank, Linköping). Isolation was performed by the method of density gradient centrifugation using Lymphoprep (Axis-Shield) and Sepmate-50 tubes (Stemcell Technologies) according to manufacturer's protocol. IGRA status was determined on the whole-blood samples using QuantiFERON-TB Gold (Cellestis) following the manufacturer's instructions.

Cell culture and in vitro BCG training
Following PBMC isolation cells were seeded in cell culture treated asks in Dulbecco's Modi ed Eagle's Medium (DMEM, Invitrogen) with 25 mM hepes, 100 U/ml penicillin and 100 µg/ml streptomycin (PEST, Gibco). Cells were incubated in 37 °C for 2 h before the non-adherent lymphocytes were washed away using warm Krebs-Ringer Glucose buffer (made in house). Complete DMEM supplemented with 10% pooled human serum, 2mM L-glutamine, 100 U/ml penicillin and 100 µg/ml streptomycin (all from Gibco) was then added to the cells along with immune training agents, or medium only as negative control for 24h. For the DNA methylation analysis and data presented in gure 7a-d we used 10 µg/ml of BCG (freeze dried M. bovis bacillus Calmette-Guérin Danish Strain). Training agents were washed off and the cells were incubated for 6-7 days in complete medium with media change every 2-3 days. The cells were washed with PBS followed by trypsinization and reseeding of 5000 cells/ well in a 384-well plate (Corning Falcon) for infection experiments. DNAm analysis DNA was isolated at day 6 after training with BCG. with subset-quantile within array (SWAN) normalization 36, 37 and ChAMP (v2. 19.3) with beta-mixture quantile normalization (BMIQ) packages 38, 39 . The type I and type II probes were normalized using the quantile normalization method. Using the default setup of the ChAMP package, following probes were ltered out: i) probes below the detection p-value (>0.01), ii) non-CpG probes, iii) multi-hit probes, and iv) all probes of X and Y chromosomes. Cell type heterogeneity was corrected for the PBMC cell types using the Houseman algorithm 40 and batch effects were xed using ComBat from the SVA package (v3.38.0). Differential methylation analysis was performed with the linear modeling (lmFit) using the limma package 41 (v3.46.0) in a contrast matrix of the TB-exposed and TB-non-exposed (Control) individuals. All Differentially methylated CpGs (DMCs) were considered signi cant at the Bonferroni-Hochberg (BH) corrected p-value < 0.05 (for HLA-DR cell types), <0.1 (for CD3 cell types) and <0.2 (for PBMC cell types). The DNA from the human primary macrophages was sequenced using the Diagenode ś RRBS due to a lower DNA yield. The sequenced reads were quality checked using the FastQC 42 (v0.11.9). The sequences were trimmed to remove arti cially lled-in cytosines at the 3′ end using the TrimGalore (v.0.6.5)( https://github.com/FelixKrueger/TrimGalore) with a phred score cutoff of 20 and quality checked again after trimming. The trimmed sequences were aligned with the human reference genome (hg38.13) using Bowtie2 43 and removed the duplicates using the Bismark (v.0.22.3) 44 . The methylation extractor from Bismark was used to extract the CpG methylation data from the sequences. The SAMtools (v1.7) package 45 was used to sort the binary alignment les (BAM) les on CpG-site chromosomal location and converted to sequence alignment map (SAM) les. The methylated and unmethylated CpG counts were extracted and combined using the DMR nder 46 (v0.3) package in R (v4.0.2). The CpG-sites located in the X and Y chromosome as well as CpG-sites from mitochondrial DNA were ltered out. For differential methylation analysis the methylKit 47 (v1.18.0) package was used. CpG-sites with a read coverage > 10 with both methylated and unmethylated reads were removed from the analysis. The function "calculateDiffMeth" from methylKit package using logistic regression was used as statistical test with "before" and "after" BCG-training as treatment used to identify differentially methylated CpG-sites DMCs. The DMCs were annotated to the o cial gene symbol using org.Hs.eg.db (v3.12) and AnnotationDbi (v1.52) packages using the human genome version hg38. DMGs were de ned as sited with a methylation difference of 30% between the conditions and an BH corrected p-value of <0.01.

Unsupervised cluster analysis
Hierarchical clustering of the all TB-exposed and control individuals was performed with the normalized β-values obtained after the data ltration in each cell type individually. The distance was calculated using the Euclidean distance matrix. The dendextend 48 (v1.14.0) and ape 49 (v5.4-1) packages in R were used to construct the horizontal hierarchical plots from the three different cell populations using the hclust and dendrogram functions.

Structural annotations
The EnhancedVolcano package 50 (v1.8.0) was used to generate the individual volcano plots from all cell populations. The ChromoMap package 51 (v0.3) was used to annotate and visualize the genome-wide chromosomal distribution of the DMGs. The interactive plots were generated using the plotly (v4.9.3) package.
The heatmaps were generated from the ltered DMGs with their respective CpGs for each cell type using the ComplexHeatmap (v2.6.2) package 52 . The clustering dendrogram in heatmaps were plotted using the Euclidean distance matrix.

Pathway and functional enrichment analyses
We used the PANTHER database (PantherDB v15, 16) 53 to identify the enriched pathways related to our identi ed DMGs. In addition, to assess functional enrichment, we used the ReactomePA (v1. 34

Venn and overlap analyses
Venn analyses were performed to detect the DMGs overlapping between cell populations and between studies. We constructed the Venn diagrams by using matplotlib-venn package (https://github.com/konstantint/matplotlib-venn) using in-house python script. The overlap analyses were calculated and plotted using the go.Sunburst function from plotly using an in-house python script.

Statistical analyses
All differences with a p-value < 0.05 were considered signi cant if not otherwise stated. We calculated family-wise error rate (FWER) using the BH correction method. All analyses were performed in R (v4.0.2) with the mentioned packages.

Declarations Data Availability
The datasets generated and analysed during the current study are not publicly available due to Intellectual property rights but are available from the corresponding author on reasonable request.  Figure 1 Schematic overview of the project work ow. Sputum and blood were isolated from TB-exposed individuals and controls and the DNA methylomes were analyzed for the different cell types.
DNA methylome analyses of immune cells from TB-exposed individuals. Volcano plots of DMGs from d. HLA-DR, e. CD3 and f. PBMCs. Red dots represent DMGs above cut-offs (±0.3 Log 2 fold change and BHcorrected p-value < 0.05, <0.1 or 0.2 as indicated). DMGs, differentially methylated genes.  Venn analyses comparing DMGs, pathways and GO terms between different datasets. a. Overlapping DMGs derived from the HLA-DR (orange), CD3 (green) and PBMC (red) DNA methylomes. b. Overlapping DMGs from this study and from our previous work on BCG induced DMGs (light green).