DNA methylation changes associated to interferon-gamma during classical to non-classical monocyte differentiation
To investigate the regulatory mechanisms leading monocyte subsets differentiation, we initially purified by flow cytometry the three monocyte subsets (Figure 1A) from healthy donors following a negative gating strategy similar to those published previously [4,19,20] (Figure 1B). Bisulphite modified DNA obtained from the resulting sorted cells was used to analyze the DNA methylation profiles of the three subsets by Infinium MethylationEPIC V1.0 bead arrays. The analysis revealed a total of 2625 CpG positions differentially methylated between at least two of the subsets (FDR < 0.05 and differential of beta > 0.2) (Figure 1C and Additional File 1. Table S1). Hierarchical clustering of the methylation patterns at these positions enabled the discrimination of two distinct clusters. Cluster M1 comprises 1905 differentially methylated positions (DMPs) that experience DNA demethylation during the differentiation from cMO to ncMO. Cluster M2 consists of 720 DMPs that undergo gains in DNA methylation from cMO to ncMO. In both clusters, the iMO subset displays intermediate methylation levels between the other two subsets (Figure 1C).
Principal component analysis of the DMPs revealed a higher similarity between iMO and ncMO than with cMO (Figure 1D). CpGs from both clusters annotated preferentially to regions outside CpG islands while annotation to CpG islands was clearly underrepresented in comparison to the background (Figure 1E, left). Concordantly, annotation to genomic regions showed a clear reduction in promoter regions representation in both clusters in favour of an increase in distal intergenic and intron regions (Figure 1E, right). Enrichment analysis of chromatin functional states showed that DMPs from both clusters were significantly enriched in enhancers as well as significantly underrepresented in transcription start sites (TSS) and repressed regions (Additional File 2. Figure S1A)
Gene ontology analysis revealed that CpGs from both clusters are associated to immune response functions such as leukocyte activation or regulation of immune response (Additional File 3. Table S2). CpGs from cluster M1 were repeatedly associated with myeloid activation functions. Notably, the ontology terms “response to interferon-gamma" and “positive regulation of interferon-gamma production” were significantly enriched with an adjusted P-value < 0.05 exclusively in the DMPs from cluster M1 (Figure 1F). Among the genes annotated by proximity to these DMPs we found some with a functional relevance in immunological context such as FLT1, HCK, PRDM1 or PTGS2 for cluster M1 (Additional File 2. Figure S1B). Also, CpGs from cluster M2 were enriched for “positive regulation of fractalkine biosynthetic process” function, which has been described to be relevant for reducing apoptosis levels in the different subsets of monocytes. In cMO particularly, it achieves this by reducing intracellular levels of reactive oxygen species (ROS)[46]. Some example genes annotated to these DMPs were SPI1, TLE1, ITGB7 and SBNO2 (Additional File 2. Figure S1B).
TF binding motif enrichment analysis of the regions surrounding the identified DMPs showed a clear representation of known myeloid lineage TFs (Additional File 4. Table S3). For example, cluster M1 presented an overrepresentation of factors such as PU.1 and KLF4 [47,48] as well as several interferon regulatory factor family members. In cluster M2 we observed a redundant representation of CEBPs [47], as well as NFIL3 [49]. It is interesting to note that many of the TFs enriched in both clusters have previously been described to play a crucial role in myeloid lineage and differentiation. However, different TFs appear to have major roles in the different subsets. Importantly, interferon regulatory factors (IRF) such as IRF1, IRF2, IRF3 and IRF8 were enriched only in cluster M1 indicating, again, a potential demethylation-associated activation of the pathways regulated by interferon in the positions associated with demethylation in ncMOs (Figure 1G). Together, these results reinforce the idea that the three monocyte subsets are both part of a continuous differentiation process and, at the same time, different enough to deserve individual characterization.
Transcriptomic reprogramming from classical to non-classical monocytes unveils antagonistic roles for HIF-1a and T-bet
We further characterized the three monocyte subsets by profiling their transcriptomes. RNA-sequencing analysis of the paired subsets from three healthy donors revealed a behaviour similar to that already observed in DNA methylation. That is, cMO and ncMO exhibited the most distinct profiles, while iMO displayed an intermediate phenotype. In particular, 3035 genes presented differential expression between at least two of the subsets (differentially expressed genes, DEG, adjusted P-value < 0.05 and fold change < 0.5 for downregulated or > 2 for upregulated) (Figure 2A, Additional File 5. Figure S2A and Additional File 6. Table S4). Specifically, 1617 genes had decreased expression levels in the differentiation from cMO to ncMO (cluster E1), while 1412 DEGs showed increased expression levels (cluster E2).
Gene ontology (GO) analysis of the DEGs showed that cluster E1 genes (overexpressed in cMO) were enriched for pathways related to functions previously described for cMO subset. Response to bacterium, response to lipopolysaccharide or angiogenesis were among the GO terms enriched in this cluster. Cluster E2 genes, on the other hand, were enriched in functions repeatedly related to T cell activation and cell-cell adhesion which is also concordant with previously described functions for the ncMO subset (reviewed in [6] and [7]). Interestingly, in agreement with the analysis of ontology terms of the DNA methylation data, interferon gamma-production related terms appeared significantly enriched in cluster E2 (Figure 2B and Additional File 7. Table S5).
Using DoRothEA (Discriminant Regulon Expression Analysis [41]), we predicted TFs activity through the expression of their target genes (Figure 2C and Additional File 8. Table S6) in the differentiation from cMO to ncMO. This analysis identified Hypoxia Inducible factor 1 (formed by a heterodimer of HIF-1α and ARNT) and PU.1 (SPI-1) as the TFs with the highest activity in cMO. PU.1 is a key regulator of monocyte differentiation, also predicted by motif enrichment analysis of CpGs hypomethylated in cMO (Figure 1G, cluster M2). HIF1 is a crucial factor that regulates cellular response to low oxygen concentrations or hypoxia. Also in monocytes, HIF-1α regulates many aspects of the cell’s phenotype and metabolism under these stress conditions [50]. The observed decrease in HIF-1 regulons, as well as the expression of the factors themselves (Figure 2D and Additional File 5. Figure S2B), could potentially be explained by the increase in oxygen concentrations encountered by monocytes during their differentiation process. In the bone marrow, where cMO originate, there is less than 5% of oxygen concentration. In the bloodstream cMO only spend 24h, during which the oxygen concentration varies between 12-20%. Conversely, ncMOs have resided in the bloodstream for at least 5 days, during which they never encounter hypoxic conditions [51,52]. Therefore, it is expected that cMOs have a phenotype adapted to hypoxic conditions, regulated by HIF1, which is lost during the differentiation towards ncMOs. Other TFs that appear to have an enriched activity in cMO based on the differential expression between cMO and ncMO include NFκB, AhR and MAFB. These TFs are very important for the activation of monocytes under different conditions.
Interestingly, regulon analysis of the ncMO transcriptome revealed the participation of SNW1-associated protein 4 (SNAPC4), a factor relevant for mRNA splicing, T-box TF TBX21, a well-known lineage defining TF from naïve Th precursors to Th1 [53], and a Max-binding protein from the network Myc/Max/Mad called MNT. Considering these results, reports associating ncMOs with T cell regulation and communication and given that TBX21 is the only one among the top predicted TFs to exhibit differential expression between cMOs and ncMOs (see Figure 2D and Additional File 5. Figure S2B), we decided to further investigate the role of TBX21 in this differentiation process.
TBX21 encodes the TF known as T-bet. This TF is preferentially expressed in T cells and natural killer (NK) cells. In fact, previous studies showed that it is not expressed in monocytes at basal state, but its expression can be induced upon treatment with IFN-γ [54]. Other studies indicate that it is overexpressed in human M1 macrophages [55]. TBX21 was shown to be significantly overexpressed in ncMO in comparison to cMO in our transcriptomic data (Figure 2D). These results were validated at the protein level, indicating a potential contribution of T-bet in ncMO, and maybe in iMO, that would be differential to cMO (Figure 2E). The antagonistic activity of HIF1 and T-bet activity in the subsets was also validated by the determining the differential expression of most of their known target genes (Additional File 5. Figure S2C), suggesting a contribution of these factors to the differential transcriptome of the monocyte subsets.
We then performed an integration of our DNA methylation and expression data to test whether DNA methylation can regulate gene expression in our datasets. A significant inverse correlation was seen between DMPs-associated genes by proximity and DEGs (P-value = 3.92e-13, correlation coefficient [r] = - 0.262) (Figure 2F). Genomic context of the DMPs has been proved to be an important factor to consider when studying its impact in gene expression regulation [56]. For this reason, we studied correlation in each of the tiles we use to divide genomic annotations. Only promoter and intronic associated DMPs showed significant correlation with their annotated gene’s expression (promoters: P-value = 7.576e-9, r = - 0.3961, see Additional File 5. Figure S2D for example genes; introns: P-value = 2.771e-8, r = - 0.3164). Considering this, we studied the methylation status of the CpGs proximal to HIF1-α and TBX21 gene. Many of them showed significant differences in methylation (p-value < 0.05) indicating that the differential expression of these factors might be regulated at the epigenetic level (Figure 2G and Additional File 5. Figure S2E).
T-bet as an important TF for the epigenomic reprogramming from cMO to ncMO
Based on our results, which highlight the significance of the TF T-bet in the transcriptomic profile of ncMO (Figure 2C, right panel) and the ontology terms relating this transcriptomic profile to that of T cell activation (Figure 2B, cluster E2), we hypothesized that T-bet might have a role in the differentiation from cMO to ncMO. T-bet can be induced by a variety of extracellular signals. In NK cells, it has been proven to be induced in response to cytokines such as interferon-γ (IFN-γ), interleukin-12 (IL-12), IL-15 or IL-21 (reviewed in [57]). To test whether monocytes were equally responsive to these stimuli in terms of induction of T-bet, we performed stimulations of primary CD14+ monocytes in vitro for 6 and 24 hours (Figure 3A). Our results showed that T-bet in monocytes is mainly induced in response to IFN-γ, both at the mRNA (Figure 3B) and protein (Figure 3C) levels, which is in agreement with previous results [54]. The mRNA increase is higher at 6 than 24 hours while the protein increase is more pronounced at 24 hours. Moreover, this induction probably takes place through the IFN-γ archetypical pathway, which involves STAT1 phosphorylation (Figure 3C).
Taking advantage of available high throughput public data of the response of monocytes to IFN-γ, we studied whether part of the differentiation process from cMO to ncMO could be explained by this stimulus. In terms of DNA methylation, we used the data from GSE134425, where they stimulated PBMCs with several stimuli for the posterior analysis of the methylome of the sorted stimulated monocytes. In this study, they tested three stimuli: IFN-γ, IFN-α and TNF-α. By comparing their results to our data, we observed that ncMO are significantly more similar to the stimulated samples in all instances. In this regard, the positions that become demethylated in either of the three stimuli also display demethylation in ncMO. This is particularly relevant for IFN-γ (Figure 3D, E). Unexpectedly, this was not the case for expression data. For this analysis, we made use of the data from GSE130567 where they compared the effect of IFN-γ to primary monocyte-derived macrophages. In this case, we found that the genes significantly upregulated during this stimulation were significantly more expressed in cMO than in ncMO (Figure 3F, G). From this analysis, we conclude that ncMO have an epigenetic profile similar to stimulated monocytes. However, transcriptomically their phenotype is not as straightforwardly explained as an individual stimulus response.
T-bet deficiency affects the whole monocyte compartment
To investigate the role of T-bet in the monocyte compartment, we investigated the impact of a homozygous mutation in T-bet in an individual in comparison to a heterozygote male progenitor, used as a control, by inspecting the single-cell RNA-seq datasets from their peripheral blood (Figure 4A) [42]. The affected individual suffers from Mendelian susceptibility to mycobacterial disease (MSMD), which was directly attributed to the indel variant in homozygosity in the TBX21 sequence, replacing two amino acids that are highly conserved across species. The control progenitor was not reported to present any affectation attributed to its heterozygotic profile. This causes a loss-of-function allele that, although it allows expression of T-bet itself, it is unable to induce the expression of its target genes. The first effect we observed at the single-cell RNAseq level was a variance in the proportions of different cell types between the affected individual and his control progenitor (Additional File 9. Figure S3A). As expected, there is a marked reduction in the contribution of the mutant individual on CD8 Naive T cells and NK cell populations. These are two cell types where T-bet factor is known to be crucial for their development [57–59]. More surprisingly, we discovered an equal reduction in the myeloid compartment, suggesting a role for this factor to the proper development of this lineage as well.
By focusing on the deficient myeloid compartment, we identified clusters of cells with a transcriptomic signature similar to that observed in our bulk RNA-seq analysis. The cMO and ncMO signatures were defined as the top 200 genes with the lowest p-value in the comparison between the two subsets, in each direction of fold change. These signatures identified a region of the object clearly expressing ncMO signature, while the vast majority of the myeloid compartment was more enriched for the cMO signature (Figure 4B). Interestingly, this region forms a distinct cluster that separates from the rest of the myeloid compartment even at very low resolutions, and it is notably less abundant in the affected individual compared to the control progenitor (see Figure 4C). Additionally, there is another cluster of cells in the cMO region that becomes apparent at higher resolution and is also deficient in the mutant individual.
We then performed a discriminant regulon expression analysis between the control and the affected individual in the whole myeloid compartment with the DoRothEA dataset [41] (Figure 4D). This analysis revealed that one of the top TF responsible for the differences was HIF-1α. This suggest that the deficiency of T-bet activity also affects the proper activation of the HIF-1 pathway so that in the absence of T-bet there is also a reduction of the expression of HIF-1α target genes (Figure 4E and Additional File 9. Figure S3B). Other TFs that appeared to be negatively altered in the absence of T-bet included the NF-kβ pathway with components such as NFKB1, RELB and RELA. Additionally, NFKB2 was presented as the main TF affected in the mutant individual. This pathway is very tightly associated with HIF signalling in macrophages [60]. Also, ATF6 appeared as a top TF differentially implicated across the samples. This TF is activated upon endoplasmic reticulum stress which associates with NF- kβ activation in many circumstances [61].
Finally, we studied what impact could T-bet dysfunction have on the phenotype of these monocytes. To this end, we studied the predicted interactions they would have with the rest of the immune cells. We focused our study on the effect T-bet deficiency has on the ncMO-like cluster of cells, due to our demonstration in bulk of T-bet expression in cells with this transcriptomic profile. For this analysis, the myeloid object was divided into several clusters that were annotated based on their top differential markers and the functions enriched within them (Additional File 9. Figure S3C-E). By using CellChat software [45], which explores cell-cell communication based on the expression of ligand-receptor pairs, we identified a few cell communication pathways that were altered in the individual carrying the homozygous mutations. In particular, the affected individual’s ncMO cells expressed higher levels of CCL signalling pathways receptors, making them more susceptible to receiving this type of signalling from a broad plethora of immune cell types. This includes a potential positive feedback loop among themselves. On the other hand, these affected ncMO appeared to have fewer levels of TNF receptors than their control counterparts, making them less able to respond to this inflammatory signalling (Figure 4F and Additional File 9. Figure S3F). Interestingly, one of the pathways responding to TNF signaling is the canonical NF-kβ pathway [62,63] which would correlate with the previous analysis where multiple components of this signaling pathway were altered in the comparison between the control and the affected individual.