Multi-omics computational analysis unveils the involvement of AP-1 and CTCF in hysteresis of chromatin states during macrophage polarization

Macrophages display extreme plasticity, and the mechanisms and applications of polarization and de-/repolarization of macrophages have been extensively investigated. However, the regulation of macrophage hysteresis after de-/repolarization remains unclear. In this study, by using a large-scale computational analysis of macrophage multi-omics data, we report a list of hysteresis genes that maintain their expression patterns after polarization and de-/repolarization. While the polarization in M1 macrophages leads to a higher level of hysteresis in genes associated with cell cycle progression, cell migration, and enhancement of the immune response, we found weak levels of hysteresis after M2 polarization. During the polarization process from M0 to M1 and back to M0, the factors IRFs/STAT, AP-1, and CTCF regulate hysteresis by altering their binding sites to the chromatin. Overall, our results show that a history of polarization can lead to hysteresis in gene expression and chromatin accessibility over a given period. This study contributes to the understanding of de-/repolarization memory in macrophages.

1 Introduction

Macrophages, an essential cell type of the innate immune system, migrate into different tissues from the peripheral blood and undergo differentiation upon receiving various immune cues (1). The major function of macrophages is to bridge the innate immune system with the adaptive immune system through several biological processes, such as phagocytosis of pathogens/cell debris, antigen presentation on major histocompatibility complexes, and production of cytokines (1–3). Non-activated M0 macrophages polarize into either the classical pro-inflammatory M1 phenotype or the alternative M2 phenotype. M1 is typically induced by IFN-γ or Toll-like receptor (TLR) ligands and produces redox molecules (4). M2 macrophages are induced by exposure to a set of cytokines and are observed in conditions, such as helminth infection (5), wound healing, cancer, and obesity (6). The polarized macrophages express pro-inflammatory cytokine genes, including IL-1β, TNF, IL-12, and IL-18 in M1; and Chi3l1, Ym1, Rentla, and Arg1 in M2 (7).

Interestingly, macrophages can retain their plasticity. This means that M1 cells exposed to M2-inducing signals can differentiate into M2 states and vice versa, which often influences the outcomes in therapeutic treatment of diseases, such as cancers (8). The transcriptional phenotypes in the cellular microenvironment and subsequent immune response are influenced by not only the ongoing de-/repolarization (9) but also the memory of the past cellular state, which is referred to as hysteresis (10, 11). Specifically, such hysteresis manifests as a system reluctance of macrophages to immediately polarize into new phenotypes when changes occur in the stimulus cues, preserving its previous cell state for a period of time. This results in the macrophage system’s response curve exhibiting a hysteresis loop or closed-loop shape, rather than a simple linear relationship. For instance, trained immunity is a typical phenomenon that may be attributed to hysteresis (12), and recent experimental evidence supports the significant influence of hysteresis on immune responses (11).

However, the molecular mechanisms underlying macrophage plasticity and cell polarization-mediated hysteresis remain unclear. Therefore, in this study, we performed a large-scale computational analysis to investigate the regulatory mechanisms of macrophage hysteresis during polarization. Consequently, we identified hysteresis-related genes that are enriched among interferon-stimulated genes (ISGs) and involved in fundamental immune functions. Remarkably, the hysteresis induced by M1 stimulation was more pronounced than that induced by M2 stimulation. We demonstrated that specific DNA-binding proteins are tightly associated with different hysteresis strengths. AP-1 and CTCF were consistently enriched in accessible and inaccessible distal regions, respectively, during M1 hysteresis. These findings suggest an impact of chromatin remodeling on the regulation of hysteresis genes.

2 Materials and methods

2.1 Data preparation

The time-course bulk RNA-seq, single-cell RNA-seq (scRNA-seq), bulk ATAC-seq data, and single cell ATAC-seq (scATAC-seq) used in this study were sourced from the Liu Project (9) (GSE158094) and Daniel project (13) (GSE178526). The scRNA-seq data used to train logistic regression classifier were obtained from the Andres project (GSE161125) (14) and the Chuan project (GSE117176) (15). The ChIP-seq data employed for enhancer prediction were acquired from the ENCODE project (GSM1000074). The ChIP-seq data used to validate the role of AP-1 were obtained from the Gene Expression Omnibus (GEO) database (GSE38377, GSE36099, GSE36104, GSE84519, and GSE84520). All the data in this study were obtained from mouse bone marrow-derived macrophages (BMDM) in the in vitro condition.

2.2 NGS data processing

2.2.1 Bulk RNA-seq data analysis

The log₁₀ normalized fragments per kilobase million (FPKM) mRNA expression matrices previously processed using Cuffdiff (16) by Liu et al. from the primary (M0), polarized (M1 and M2), de-polarized (deM0_M1, deM0_M2), and re-polarized (reM1_M2 and reM2_M1) macrophages were directly downloaded from GEO and used for downstream analysis. Differentially expressed genes (DEGs) were selected based on average difference of log₁₀(FPKM) > 0.5, P values < 0.05, false discovery rate (FDR) < 0.05 between M0/M1/M2 (control) and 96 h de-polarized M0/re-polarized M1/M2, respectively. P-values were calculated using Student’s t-test, and FDRs were calculated using Benjamini-Hochberg correction.

2.2.2 Single-cell RNA-seq analysis

Seurat (version 4.2.0 (17)) was used for preprocessing and downstream analysis of the macrophage scRNA-seq data from the Liu project (9). First, the macrophage phenotype annotation provided in the original article was used to differentiate single-cell macrophages. Thereafter, the selected cells with 200–2500 detected genes, > 200 nUMIs, and < 10% mitochondria genes were used for quality control and performed downstream analyses. Subsequently, 22752 macrophage cells with totally 4085 gene expressions were extracted and normalized using the standard Seurat protocol. The top 3000 variable features were selected, and the cells were clustered and visualized in the uniform manifold approximation and projection (UMAP) plot. Canek (18) was used to remove batch effects from the single-cell data. Raw counting matrices from the Andres and Chuan projects were used directly for logistic regression without any additional processing.

2.2.3 Bulk ChIP-seq and ATAC-seq analyses

The raw ChIP-seq/ATAC-seq SRR files were downloaded from Sequence Read Archive (SRA). After converting to FASTQ format, reads were checked and trimmed by FastQC (19) and Trim Galore (20). Read alignments were performed against the GRCm38/mm10 reference mouse genome using Bowtie2 (version 2.4.5) (21), with the default settings. BAM (*. bam) files were generated and merged using samtools (1.3.1) (22). Low-quality alignments and mitochondrial reads were detected and removed using the Sambamba software (0.7.1) (23); peaks were called by MACS2 (2.1.1) (24). Plotting of ChIP-seq/ATAC-seq signal abundance in the transcription start site (TSS) region, conversion of files to BigWig format, and RPKM normalization were performed using deepTools (25) for subsequent analysis using default settings. Specifically, we utilized the bamCompare function to convert BAM files into BigWig format while normalizing to RPKM values. Then, we addressed the BigWig files using the function computeMatrix in deepTools, and generated the overall accessibility profiles of macrophages in different cellular states within the promoter regions (1000 bp region around TSS), which were visualized by the plotHeatmap function. The peaks in the BED files of the genome were visualized using IGV (26).

2.2.4 Single-cell ATAC-seq analysis

We downloaded five fragment files from the GSE178526 including CTR (0 h M0), Pre.IL4 (48 h deM0_M2), IL4 (24 h M2), Pre.IL4+IFGN (51 h reM1_M2), and IFGN (3 h M1). After merging fragment files, functions called CallPeaks and FeatureMatrix in Signac v1.7.0 R package (27) are applied to obtain significant accessible peaks from fragments and construct the peak matrix, respectively. To include the cells used in the original research, we extracted the cell barcodes of IL4 in fig5_cellColData_trim.csv file and combined them with the cell barcodes of CTR, Pre.IL4, Pre.IL4+IFGN, and IFGN in fig6_cellColData_trim.csv file. A subset of the peak matrix, which only contains 186,413 peaks and 30,474 cells, is then leveraged as the input of the downstream analysis. FindMarkers function in Signac is used to find differentially accessible peaks between two cell groups with all suggested parameters. GeneActivity function suggested by Signac is used to generate gene activity matrix by computing counts per cell in gene body and promoter region based on mm10 reference. As a result, gene activity matrix with 21,808 protein-coding genes and 30,474 cells is obtained. Singular Value Decomposition (SVD) was applied to the term frequency-inverse document frequency (TF-IDF)-coded peak matrix, with the second to thirtieth SVD components are visualized by UMAP.

2.3 Logistic regression on scRNA-seq data

Logistic regression classification was performed using the scikit-learn (28) Python package. Raw scRNA-seq counts from the Andres and Chuan projects were used as inputs. Eighty percent of the samples were selected as training data and the remaining twenty percent was selected as testing data. An OVR (one vs. rest) strategy was used to identify the subtypes of macrophages. The hyperparameters were selected using the GridSearchCV strategy to minimize the loss function and achieve the best performance. The cell-type label provided by the original research was considered the ground truth to calculate the accuracy and F1 score. Finally, the feature genes and their feature importance were selected from the trained classifier.

2.4 Network construction, graph embedding, and clustering analysis

To construct the ground truth of the TF-Gene interaction network for mice, gene co-expression data were collected from COXPRESdb (29), and the co-expression data were filtered based on the co-expression score provided by the database. Specifically, gene pairs with co-expression scores > 5 were retained. Confidence in the TF-Gene regulation data was acquired from the TRRUST database (30). Subsequently, the information from the two data sources was merged and a sub-network of genes, which only existed in our time-course bulk RNA-seq dataset from the Liu (9) project, was selected to construct a TF–gene interaction network.

In our network, 6776 gene nodes within 382 transcription factors (TFs) were identified. If two nodes had either a co-expression or a regulatory relationship, an undirected edge connecting between the nodes was added to the graph. Given the undirected graph G = [V, E], where V represents all selected genes and E represents all co-expression and TF regulation relationships; the weight of E was not considered in our study. Graph embedding analysis was performed using DeepWalk (31). For an unweighted graph, by performing a random walk method, DeepWalk could resemble the node sequences in the graph, which are then input into Word2Vec to perform the graph embedding. OpenNE (https://github.com/thunlp/OpenNE), which is an open-source Python package, was used to apply the DeepWalk embedding method. The K-means clustering method was directly applied to the embedding matrix, where the K value was chosen based on the number of genes in each cluster, with 16 selected as the final K value. Subsequently, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis were performed using ClusterProfiler (32) for cluster function annotation. GO terms and KEGG pathways with corrected P-values < 0.05 were selected as significant annotations. The enrichment of the hysteresis genes in each cluster was calculated using Fisher’s Exact Test.

2.5 Transcription factor binding site, enhancer prediction, and hysteresis peaks identification

UniBind (33) was used to predict the potential TFBS among all selected ATAC peaks. The selected DEGs were uploaded to the ImmunoNavigator (34) online tool to predict potential TFs that regulate gene expression on promoter regions. The enhancer-gene interactions were predicted using the activity-by-contact (ABC) method (35), following documentation from the ABC repository in GitHub and setting up the conda (Anaconda, 2020) environments and suggested packages. Subsequently, both the processed ATAC-seq peaks and H3K27ac peak files were input. Candidate Region function with the recommended parameters (nStrongestPeaks 150000 and peakExtendFromSummit 250) was used. Enhancer activity was quantified using the ABC run.neighborhoods function with the mm10 gene as a reference. Finally, the ABC scores of enhancer-gene interactions were computed with the ABC prediction function setting –scale hic_using_powerlaw and a threshold of 0.2.

To define hysteresis peaks based on ATAC-seq, we selected all significant peaks from M0(0 h), M1(24 h) and deM0_M1(96 h) data. Among all significant peaks, the following two types of hysteresis regions were defined: the open-close-close (OCC) region, representing the regions that are highly accessible at M0(0 h), but significantly less accessible at M1(24 h) and deM0_M1(96 h); and the close-open-open (COO) region, representing the regions that are significantly less accessible at M0(0h), but significantly accessible at M1(24 h) and deM0_M1(96 h). For the OCC region, significant peaks exist in M0(0 h), but do not exist in M1(24 h) and deM0_M1(96 h) were kept and were taken the union. For the COO region, we kept and took the union of the peaks that exist in M1 (24 h) and deM0_M1 (96 h) but do not exist in M0 (0 h). Then, we selected ATAC peaks that overlapped these regions as OCC/COO peaks for downstream analysis using BedTools (v2.30.0) (36). The significance of the results was determined using the chi-square test.

2.6 Other bioinformatics analysis and statistical analysis

The CTCF-binding annotation was downloaded from ENCODE (https://www.encodeproject.org/data/annotations/). Protein–protein interaction (PPI) analysis was performed based on the PPI data from the STRING database (37). P<0.05 was considered statistically significant, and the other p-value conclusions were as follows: “-”: P> 0.05, “*”: P< 0.05, “**”: P< 0.01, “***”: P< 0.005, and “****”: P< 0.001.

3 Results

3.1 Differential gene expression analysis reveals the existence of hysteresis genes after de-/repolarization

To ensure that the macrophage data we used had undergone thorough de-/repolarization well, as an example, we presented the RNA expression of key genes and selected M1/M2 marker genes during the M0->M1->M0 polarization process (Figures S1, S2) and the trajectory plot base on overall transcriptomic profile (Figure S3). These results perfectly replicate the findings of Liu et al.’s research. Then, we used the time-course RNA-seq data (Figure 1A) obtained from the Liu project (9) to investigate the influence of de-/repolarization on the macrophage transcriptome and chromatin accessibility. This dataset employed M0 cells as the starting point (0 h), followed by stimulating the cells with either M1 or M2 cues for 24 h to yield M1 and M2 cells (24 h). After this 24-hour stimulation, the cells were either washed with wash buffer for 72 h to facilitate de-polarization (de-polarization group) or subjected to the opposite stimulus for 72 h to facilitate re-polarization (re-polarization group) (Figure 1B). Thus, we established the following four cellular states: de-polarized M0 from M1 (deM0_M1), de-polarized M0 from M2 (deM0_M2), re-polarized M1 from M2 (reM1_M2), and re-polarized M2 from M1 (reM2_M1). The expression alterations observed in the transcriptome due to de-/repolarization within the same macrophage phenotype were collectively defined as hysteresis in gene expression. Genes exhibiting significant expression discrepancies were identified as hysteresis genes, which likely mirror the past state of the cells.

Figure 1

Figure 1 Multi-omics data used for macrophage analysis and explanation of macrophage hysteresis. (A) Schematic illustration depicting the macrophage response to cytokine treatment and acquisition of multi-omics data. The horizontal time axis signifies distinct time points (in hours), and the arrows symbolize various types of stimuli. The labels 0 h, 24 h, and 96 h correspond to the pivotal macrophage states. Positioned beneath the time axis, the dots and the accompanying text indicate instances of data collection and respective acquired data types, respectively. (B) Schematic representation of macrophage polarization in the context of the original research and delineation of macrophage hysteresis. Each distinct pattern displayed in the colored cells symbolizes the cellular state of macrophages. The colored trajectories and arrows denote the various stimuli (M1 stimulus, M2 stimulus, and wash buffer) and polarization trajectories encountered by the macrophages. The double-headed arrow signifies the variation in expression between the macrophage after de-polarization/re-polarization and baseline macrophage states, encapsulating the concept of macrophage hysteresis. (deM0_M1: de-polarize M0 from M1, deM0_M2: de-polarize M0 from M2, reM1_M2: re-polarize M1 from M2, and reM2_M1: re-polarize M2 from M1.).

We compared the expression profiles of primary macrophages (M0, M1, and M2) with those of the corresponding de-/repolarized macrophages. To elucidate the effects on the transcriptional profile of distinct polarization histories, we compared the primary M0 cells with deM0_M1 and deM0_M2 conditions. By taking into account an absolute log₁₀(FPKM) change > 0.5 and P-value threshold of 0.05, 75 (Supplementary Excel Sheet 1) and 31(Supplementary Excel Sheet 2) genes were defined as significant hysteresis genes for the M0 vs. deM0_M1 and M0 VS. deM0_M2 group, respectively (Figure 2A). We showed the RNA expression profiles of 10 hysteresis genes in the M0 vs. deM0_M1 group, and compared to marker genes, these hysteresis genes did not return to their original expression levels within 96 hours (Figure S4). Within these two sets of significant hysteresis genes, 24 common genes, including S100a8, S100a9, Retnlg, Plac8, Ccr2, and Bgn, were identified. In addition, 51 genes, represented by Ccl5, Irf7, and Ccr3, displayed significant differential expression in the deM0_M1 group compared to that in the M0 group, but not in the deM0_M2 group, thereby highlighting the existence of a stronger M1 hysteresis. In contrast, seven hysteresis genes, including Retnla, were found to be present only in the deM0_M2 group.

Figure 2

Figure 2 Identification of Hysteresis through differential expression analysis. (A–C). Expression pattern of significant hysteresis genes within (A) M0/M1/deM0_M1/M2/deM0_M2, (B) M1/reM1_M2, and (C) M2/reM2_M1. The FPKM was normalized to the z-score scale. (D) Venn diagram illustrating the intersection of hysteresis genes across 4 gene sets. (E) Expression patterns of Ccr2 over time in macrophage polarization/de-polarization. The dashed lines indicate the key cell states at 24-h and 96-h time points. (M0->M2->M1: cell transition of M0 (at 0 h) to M2 (at 24 h) then to M1(96 h), M0->M1->M2: cell transition of M0 (at 0 h) to M1 (at 24 h) then to M2 (96 h), M0->M1->M0: cell transition of M0 (at 0 h) to M1 (at 24 h) then to M0 (96 h), M0->M2->M0: cell transition of M0 (at 0 h) to M2 (at 24 h) then to M0 (96 h)).

Simultaneously, we applied identical criteria to screen for hysteresis genes between M1 and reM1_M2, and M2 and reM2_M1, illustrating two distinct types of hysteresis upon repolarization from M1 and M2. Our analysis identified 36 (Figure 2B, Supplementary Excel Sheet 3) and 93 (Figure 2C, Supplementary Excel Sheet 4) significantly differentially expressed hysteresis genes, respectively. These results underscore that macrophages in groups with a history of M1 polarization (deM0_M1 and reM2_M1) generated a greater number of hysteresis genes following de-/repolarization.

Subsequently, in addition to the previously mentioned 24 common genes between deM0_M1 and deM0_M2, an additional 26 genes, including Ccl5, Fn1, and Irf7, were found to be common between deM0_M1 and reM2_M1 (Figure 2D). In contrast, the overlap of shared genes between the groups with an M2 polarization history (deM0_M2 and reM1_M2) was limited. This observation may stem from the relatively low number of hysteresis genes generated during M2 polarization. Remarkably, Ccr2 was the only gene that showed hysteresis across all four groups. In our study, Ccr2 expression peaked in the M0 state, significantly decreased after receiving M1 or M2 stimulation, and continued to decrease during the de-/repolarization process (Figure 2E). Notably, the blockade of the Ccl2–Ccr2 axis has been associated with an elevated expression of M1 polarization-associated genes and a reduced expression of pro-M2 markers (38). Ccr2 is also a known marker for immature monocytes, implying that M0 macrophages contain immature cells expressing Ccr2, and its expression decreases after maturation and is downregulated even after de-polarization. However, few studies have investigated the role of Ccr2 in macrophage hysteresis.

Among the common hysteresis genes with M1 history, S100a8 and S100a9, which belong to the S100 family and function as Ca2+-binding heterodimers that modulate the inflammatory response of macrophages (39), were significantly downregulated. Existing evidence indicates that the impairment of S100a9 expression during macrophage differentiation can cause an imbalance in inflammation and the wound healing processes (40). A significant reduction in the expression of Plac8, which can bind to Cebpb and control macrophage development and polarization (41), was observed in both de-polarized M0 and reM1_M2 groups. Similarly, as markers for M2 macrophages, the expression of Retnla and Fn1 were downregulated in reM2_M1 cells compared to that in the M2 cells. In contrast, Ccl5 and Irf7 were the highly expressed hysteresis genes in deM0_M1 and reM2_M1 cells. These results suggest that the hysteresis gene produced by de-/repolarization may have an impact on the inflammatory response and developmental function of macrophages. Furthermore, the hysteresis induced by the M1 polarization history may be greater at the transcriptional level than that by M2 polarization history.

3.2 The expression pattern of the hysteresis genes is similar to that of traditional macrophage marker genes

In addition to the established marker genes, we examined the expression patterns of other genes at the single-cell level. To achieve this, we leveraged the scRNA-seq data from the Liu project (9), encompassing unstimulated M0 cells at 0, 6, and 24 h, along with M1 cells with LPS & IFN-γ and M2 cells with IL-4 & IL-13 at 6 h and 24 h (Figure 1A) to validate their expression levels. After implementing the data pre-processing and quality control procedures, we amassed a total of 22,881 cells. We identified seven distinct cell subpopulations by appending their phenotype annotations (M0, M1, and M2). These subpopulations comprised M0 (0 h, 6 h, 24 h), M1 (6 h, 24 h), and M2 (6 h, 24 h) (Figure 3A).

Figure 3

Figure 3 Hysteresis genes exhibit properties similar to traditional marker genes of macrophage phenotypes. (A) Umap visualization depicting gene expression patterns across different macrophage phenotypes at different stimulation time points. Hysteresis genes (Ccl5, Irf7, Ifi205, Saa3) exhibit significantly high expression in the M1 cell population, whereas Retnla and Mgl2 show high expression in the M2 population. (B) Expression levels of hysteresis genes linked to M1 history. (C) Changes in chromatin accessibility at promoter regions and gene expression for selected M1 and M2 markers, along with hysteresis genes, during the M0->M1->deM0_M1 transition. (D) Scaled feature importance in the logistic regression model for each of macrophage subtypes, was evaluated using the Andres dataset (14) and Chuan dataset (15).

Subsequently, we evaluated the expression levels of hysteresis genes, particularly those with a history of M1 polarization. Some hysteresis genes (Ccl5, Oas3, Oasl1, Oasl2, Ifi47, Ifit2, Isg15, Ifit1, and Mx1) shared by both the deM0_M1 and reM2_M1 groups with M1 polarization history exhibited notable upregulation in the M1 phenotype (Figure 3B). Particularly, Ccl5 and Irf7 were consistently highly expressed in the M1 phenotype. Similarly, Ifi205 (a conventional M1 marker) and Saa3, which are associated with M2 polarization, showed elevated expression within the M1 cluster (Figure 3A). The elevated expression of these genes within the M1 cell subset underscores the influence of the M1 polarization history on the comprehensive transcriptome of de-/repolarized macrophages. Interestingly, among the common hysteresis genes in cells with an M2 polarization history (deM0_M2 and reM1_M2), fewer hysteresis genes (Retnla and Mgl2) exhibited a notably strong inclination toward high expression in the M2 cell population (Figure 3A). This aligns with our previous observation that the hysteresis effect is more pronounced in cells with an M1 polarization history than in those with an M2 polarization history.

We further analyzed whether the changes in the accessible chromatin regions of hysteresis genes were consistent with changes in gene expression levels by analyzing the bulk ATAC-seq dataset (Figure 1A). After processing the data, we performed a comparative analysis of the changes in the chromatin accessibility of the promoter regions and RNA expression levels among the three states M0, M1, and deM0_M1. Although hysteresis genes exhibit sustained expression profiles, changes in chromatin accessibility in their promoter regions do not perfectly mirror the observed gene expression trends, akin to conventional marker genes. This suggests that the regulation of hysteresis gene expression may involve mechanisms beyond promoter region activation (Figure 3C).

Given the absence of apparent hysteresis in gene expression and chromatin accessibility among the marker genes, we investigated the significance of hysteresis genes in macrophage subtype classification. Employing a logistic regression classifier, we aimed to determine if hysteresis genes can be used to accurately classify macrophage phenotypes. We initially selected sixteen hysteresis genes (Ccl24, Retnla, Mgl2, Fn1, Cst7, Irf7, Plac8, Oas3, Lcn2, Ctla2a, Sell, Clec4n, Ccl5, Saa3, Rtp4, and Cxcl9) as M0, M1, and M2 subtype features. Subsequently, a logistic regression model was trained using single-cell data from the Andres (14) and Chuan (15) research projects. Sixteen M1 and M2 markers (Cxcl2, Ccl3, Stat1, Nfkb2, Tap1, Cxcl9, Cxcl10, Ccl8, Tlr1, Il1R, Tnf, Arg1, Mrc1, Nos2, Tgm2, and Cd163), along with the PCA results, were used as positive controls. To establish a baseline, we randomly selected 16 genes from a single-cell dataset 100 times and averaged their performances. A grid search strategy was employed to optimize the hyperparameters for improved performance and mitigate overfitting concerns (refer to Methods). To achieve the optimal performance, we compared the prediction accuracy and F1 score (Table 1) derived from the confusion matrix (see Figure S5). The results indicate that hysteresis (genes) yielded high accuracy (0.87 in Andres data and 0.81 in Chuan data) and F1 scores (0.88 in Andres data and 0.79 in Chuan data) in comparison with randomly selected genes. However, this performance was surpassed by that of traditional M1/M2 macrophage markers (ACC: 0.89, 0.84; F1: 0.90, 0.83, respectively). Consistent feature importance trends (Figure 3D) across these genes in the model highlighted the reproducibility of classification outcomes; this verified that, while not achieving parity with marker genes, high accuracy underscored the unique expression profiles of hysteresis genes in M0, M1, and M2 macrophages.

Table 1

Table 1 Result of the logistic regression model.

These findings demonstrated that similar to conventional macrophage markers, most hysteresis genes exhibit a propensity for elevated expression in either M1 or M2 subtypes. Notably, what distinguishes them is their ability to sustain expression patterns from the preceding state following de-/repolarization. This sustained expression/non-expression phenomenon could potentially exert an additional influence on functional modifications in de-/repolarized macrophages.

3.3 Embedding of gene regulatory network revealed the functional impact of hysteresis genes on essential immunity, cell cycle, and cell migration in macrophages

To elucidate the roles of hysteresis genes, we performed functional enrichment analyses of the four distinct gene sets introduced in Section 3.1. GO and KEGG analyses (Figure 4A) revealed that hysteresis genes within the deM0_M1 and reM2_M1 groups exhibited pronounced enrichment in immune-related functions and pathways associated with M1 polarization, such as Type I interferon, NOD-like receptor, and TLR signaling pathways. Within the reM2_M1 group, we observed a notable enrichment of genes implicated in nuclear division and cell cycle functions. Conversely, cells with a history of M2 polarization showed limited enrichment, potentially due to the comparatively small gene pool. However, the study by Daniel, Belk et al. also demonstrates that the cell cycle affects macrophage plasticity. Specifically, macrophages lose their plasticity to IFNG during the S-G2/M phases, while IL-4 can induce a specific gene signature during S-G2/M (13). Notably, a shared enrichment of genes linked to cell chemotaxis was discernible across all four groups, implying a plausible connection between cellular migratory capacity and polarization history.

Figure 4

Figure 4 Hysteresis genes are enriched in gene clusters that influence multiple functions of macrophages. (A) GO and KEGG enrichment results measured in adjusted P-value in negative log-10 scale (-log₁₀ Padj) for the four groups of hysteresis genes. (B) Schematic representation of graph embedding method and downstream analysis. A gene co-expression network was constructed by integrating gene–gene interaction (GGI) data from the COXPRES database and TF–gene interaction (TGI) data from the TRRUST database. Subsequently, the gene network was subjected to DeepWalk algorithm and K-means clustering to identify gene clusters for downstream analysis. (C) Enrichment levels of the four groups of hysteresis genes were measured by P-value in negative log-10 scale (-log₁₀ P) in each of the 16 gene clusters. The gene clusters are represented on the x-axis, and the number of genes and transcription factors in each cluster on the y-axis. (D) t-SNE (TSNE) plot reveals the distribution of all gene clusters within a reduced-dimensional space following graph embedding. The accompanying bar plot demonstrates significant GO terms and KEGG pathway enrichment outcomes corresponding to each gene cluster.

The genes were then embedded into a gene regulatory network (Figure 4B). This process harnessed the existing knowledge of macrophage TFs regulation and gene–gene co-expression from publicly available databases. Our approach initially involved constructing a network graph using both the gene co-expression matrix from the COXPRES database and the TF gene regulation matrix from the TRRUST database. Subsequently, using DeepWalk, we projected this gene network graph into a lower-dimensional vector space to enable k-means clustering for predicting sub-gene regulatory clusters. Determining the number of sub-clusters entailed a joint consideration of the gene/TF numbers and cluster functional annotations. Oversized (N > 800) or undersized (N < 30) gene clusters were meticulously avoided to ensure mutually exclusive cluster annotations. Finally, we chose 16 clusters (denoted as C1, C2,…, C16) for subsequent analysis. Among these 16 clusters (Figure 4C), all four hysteresis gene groups displayed significant enrichment in C15, which also incorporated critical TFs involved in type I interferon signaling (Irf1, Irf2, and Irf9) (42) and IFN response (Stat1, Stat2) (43). Although less pronounced than that in C15, all four hysteresis gene groups were enriched in C7, within Pparg which regulates macrophage intracellular lipid metabolism (44). Hysteresis genes linked to a common M1 polarization history (deM0_M1 and reM2_M1) were notably enriched in C4, encompassing three additional members of the IRF family (Irf4, Irf5, and Irf8). Intriguingly, within all groups, significant enrichment at C8 was exclusively exhibited by the hysteresis genes of the reM2_M1 group.

Upon visualizing the network in a low-dimensional space (Figure 4D), we observed the spatial distribution of the 16 gene clusters using UMAP. Subsequently, we examined the GO annotations and KEGG pathways associated with each sub-cluster (Figure 4D). Our observations highlighted that C4, C7, and C15 are affiliated with antigen presentation, IFNγ response, Nod-like signaling/detoxification, PPAR signaling/macrophage activation, and NF-κB signaling, respectively. Notably, hysteresis genes from the reM2_M1 group were significantly enriched in C8, which was a cluster that exerted a substantial influence on DNA replication and the cell cycle, as indicated by the functional annotations.

Subsequently, we quantified transcriptomic disparities across the four groups (M0 vs. deM0_M1, M0 vs. deM0_M2, M1 vs. reM1_M2, and M2 vs. reM2_M1). We used the Manhattan distance as a reference to gauge the overall gene expression differences. Treating each of the 16 sub-gene clusters as a distinct subset of overarching expression, we autonomously computed the expression disparities between the sub-clusters in de-/repolarized cells and the sub-cluster in the initial macrophage state. In our analysis, a sub-cluster was deemed to contribute to de-/repolarization-associated expression differences (i.e., hysteresis) if the Manhattan distance in a given comparison group exceeded the overall gene expression. Conversely, if a cluster’s Manhattan distance was low, it indicated a lack of or negative contribution to the expression differences. The results indicated that all four groups, C4, C7, C8, and C15, exhibited positive contributions to hysteresis (Figure S6). Comparatively, C8 showed substantial contributions to the overall differences compared to C15 in the deM0_M2 and reM2_M1 groups. However, the cumulative impact of minor expression differences among the genes, resulting from individual gene variations, significantly influenced overall expression disparities, despite the absence of conspicuous clustering of hysteresis genes. Collectively, these findings highlight the profound influence of differential hysteresis gene expression on essential immunity, cell cycle dynamics, and cell migration within macrophages.

3.4 AP-1 motif is enriched in hysteresis chromatin region after M1 polarization and negatively correlated with the binding of CTCF

To investigate the underlying factors responsible for the hysteresis phenomenon between primary cells and de-/repolarized cells, we employed ImmunoNavigator online tools to predict enriched TFBS within the promoter regions of hysteresis genes across the four groups. Notable enrichment of TFBS was observed in the deM0_M1 and reM2_M1 groups (Table 2), whereas no significant output was observed in the deM0_M2 and reM1_M2 groups, both of which had previously undergone M2 polarization. This aligned with our gene expression analysis, which indicated a lower expression of hysteresis genes compared with that of macrophages with an M1 polarization history. Furthermore, common TFBS were identified between significant DEGs in the deM0_M1 and reM2_M1 groups. These TFBS include Irf1, a pivotal factor in M1 polarization (45); Irf2, known for inducing pro-inflammatory responses (46); and Stat1, a key mediator of M1 polarization induced by IFN-γ (47). To evaluate the specificity of this enrichment of hysteresis genes only, we established a set of M1-specific DEGs using identical screening criteria. This yielded 365 M1-specific DEGs that were used as the control group (Supplementary Excel Sheet 5). In addition to Irf-related TFBS and Stat1, we identified enriched TFBS for NF-κB and Cebpa, which enable macrophages to modulate inflammatory responses via TLR4 stimulation. Our findings imply that IRFs-associated transcription factors and Stat1, rather than Nfkb1 or Cebpa, are potential upstream mediators of M1 hysteresis.

Table 2

Table 2 TFBS motifs enriched in promoter regions.

Next, we aimed to identify the cis-regulatory elements responsible for regulating the transcriptional identity of M1 hysteresis (Specifically, M0->M1->M0), focusing on chromatin accessibility. Liu et al. (9) involved three distinct ATAC-seq datasets corresponding to the M0, M1, and deM0_M1 states (Figure 1A). We aimed to detect changes in chromatin accessibility between primary macrophages (M0) and de-polarized macrophages (deM0_M1) by analyzing these datasets. Our investigation of the accessible regions using the three ATAC-seq datasets consistently revealed strong signals for active cis-regulatory elements near the TSS regions (Figure 5A). However, we did not observe significant changes in chromatin accessibility following de-polarization, which is consistent with the findings of the previous Liu project (9). Interestingly, both the overall ATAC-seq signal on TSS regions across the genome and TSS regions of the hysteresis-related genes with significantly differential expression showed minimal differences in accessibility. Consequently, as discussed in Section 3.2, we propose that the expression of these hysteresis-related genes may not be primarily governed by direct binding events at the TSS regions, but rather by indirect regulation through distal regions.

Figure 5

Figure 5 AP-1 and CTCF exert opposing regulatory effects on hysteresis. (A) ATAC signals in the upstream and downstream regions (3 kb) of TSSs (transcription start sites) of overall genes (heatmap and line plot) and hysteresis genes (line plot) were compared across the three cell states 0 h M0, 24 h M1, and 96 h deM0_M1. (B) The genome browser displays the potential gene–enhancer (the light blue bars) interactions (the red curves) within a 40-kb genomic region, including the genes S1008 and S100a9. It shows the ATAC signal at 0 h, 24 h, and 96 h, as well as the H3K27ac signal in untreated (UT) macrophages, macrophages stimulated with LPS for 4 h and 24 h, untreated Junb knockdown macrophages, and Junb knockdown macrophages stimulated with IFN-γ for 4 h. (C) Schematic representation of hysteresis region on chromatin. The regions that exhibit low chromatin accessibility at 0 h and high chromatin accessibility at 24 h and 96 h are defined as the close-open-open regions. Vice versa, the regions that exhibit high chromatin accessibility at 0 h and low chromatin accessibility at 24 h and 96 h are defined as the open-close-close regions. (D) Log₂(fold change) (validated by chi-square test) of the binding of AP-1 members, Irf1/2, Stat1, and Ctcf/Ctcfl in the open-close-close regions and close-open-open regions. (E) Gene expression levels of the AP-1 members and Ctcf from 0 h to 96 h. (F) Alterations in Ctcf and AP-1 members binding sites within all the accessible regions of 0 h M0 macrophages and 96 h deM0_M1 macrophages. (PLS (p): promoter-like signatures, ELS (e): enhancer-like signatures, DNase regions (d), and CTCF: CTCF-binding regions).

To validate our hypothesis, we employed the ABC method to predict enhancer regions within macrophages. Among the hysteresis genes, S100a8 and S100a9 were the most significant candidates (Figure S7). Notably, the promoter region of S100a9 showed a significant decrease in the peak signal from M0 to M1, which was maintained even after de-polarization. This pattern was consistent with the ATAC and H3K27ac signal from 0 h to 24 h. In contrast, S100a8 did not exhibit any notable changes in it’s promoter region (Figure 5B). Intriguingly, we observed sustained closure of the predicted enhancer regions near the S100a8 locus during the de-polarization phase. This observation suggested a potential link between hysteresis and the cumulative impact of sustained opening and closure dynamics within gene-associated distal regions. However, upon further investigation of the peak variations of other hysteresis genes in this dataset, we discovered that this pattern could not be generalized to all genes.

Subsequently, we performed motif enrichment analysis on the open chromatin regions, discerning two distinct types of hysteresis peak regions based on selected significant peaks (Figure 5C). The first type, termed open-close-close (OCC), demonstrated heightened chromatin accessibility in the M0 state, but diminished accessibility in both M1 and deM0_M1 states, suggesting sustained closure following de-polarization. In contrast, the close-open-open (COO) region indicated the regions of persistent openness after de-polarization. Using the UniBind TFBS extraction tool, we derived the predicted motifs associated with each peak region. We then assessed the enrichment levels of these motifs in the OCC and COO regions relative to the overall peaks of high accessibility in the M0 and deM0_M1 states. Significantly enriched or depleted motifs within the OCC and COO regions were identified, and their log₂(fold changes) compared to the initial state were computed (Figure 5D). Members of the IRFs (Irf1, Irf2), Stat1, and AP-1 (Maf, Mafb, Maff, Jun, Junb, Jund, Fos, Fosl1, Fosl2, Atf2, and Atf3) family were markedly depleted in the OCC regions, but significantly enriched in the COO regions. Interestingly, Ctcf and Ctcfl exhibited a reverse enrichment pattern, being prominently enriched in the OCC regions and depleted in the COO regions. The enrichment of Irf1/2 and Stat1 aligns with findings in this study on hysteresis gene promoters as well as the network analysis results implicating Irf1, Irf2, and Stat1 from C15 in hysteresis regulation. Collectively, these findings underscore the pivotal roles of the OCC and COO regions in the regulation of macrophage function, with AP-1, Stat1, IRFs, and CTCF occupying central regulatory positions. However, our investigation of the hysteresis gene promoter regions did not reveal the presence of AP-1 or CTCF, indicating that these transcription factors predominantly bind to distal regions for regulation, rather than to promoters. Notably, the expression levels of the AP-1 family and Ctcf showed no significant changes at 0 h and 96 h (Figure 5E), suggesting alterations in their binding activity, rather than their relative abundance during de-polarization.

To gain deeper insights into the functional interplay between the AP-1 family and CTCF, we annotated mouse CTCF-bound cis-regulatory elements (cCREs) using SCREEN (48). This annotation categorizes cis-regulatory regions within the mouse genome into distinct signatures, namely promoter-like signatures (PLS), enhancer-like signatures (ELS), DNase regions, and CTCF-binding regions (CBRs), based on distinct histone marker signals. We focused on elucidating the changes in the binding patterns of AP-1 family motifs within all peak regions, comparing the 0 h M0 and 96 h de-polarized M0 states. We observed a notable shift in the motif-binding sites of members of the AP-1 family, including Junb, Fos, Fosl1, Fosl2, Jun, Jund, and Maf (Figure 5F). Specifically, in the M0 state, these motifs were predominantly enriched in regions devoid of Ctcf binding (non-Ctcf binding regions or non-CBR regions), with Junb accounting for 73.4% (Figure S8). However, in the deM0_M1 state, these motifs were concentrated within the CBRs, with Junb occupancy of 88.1%. In comparison, the binding regions associated with Ctcf and Ctcfl exhibited no significant differences between the two states (Figure S9). These observations suggest that under de-polarization conditions, AP-1 family members progressively bind to Ctcf-associated binding sites.

To further investigate the mechanism underlying hysteresis in M2 polarization, we utilized scATAC-seq data to observe chromatin accessibility changes between M0 and deM0_M2. This dataset encompassed a total of 30,474 cells in five different cell states: 0h M0, 3h M1, 24h M2, 48h deM0_M2, and 51h reM1_M2 in a total of 186,413 significant peaks (Figure 1A). The UMAP visualized the distribution of chromatin accessibility across these five cell states (Figure S10). Density plots (Figures S11, S12) illustrated the completion of polarization to M1, M2, and depolarization. Subsequently, we compared the differentially accessible peaks between M0 and deM0_M2. We identified 57 differential peaks, of which 34 were located in gene promoter regions with distance less than 1000bp to TSS (Supplementary Excel Sheet 6). Between deM0_M2 and reM1_M2, there were only 5 differential accessible peaks, with three peaks occurring in gene promoter regions (Supplementary Excel Sheet 7). This number is significantly lower than the 151 observed between M0 and M2 (Supplementary Excel Sheet 8) and the 588 observed between M0 and M1 (Supplementary Excel Sheet 9). These findings indicate that M2 hysteresis is not only weak at the gene expression level but also at the chromatin accessibility level.

3.5 ChIP-seq analysis revealed the role of AP-1 in hysteresis

To further dissect the regulatory interplay between the AP-1 family, CTCF, and a cohort of 75 M1 hysteresis genes (M0 -> M1 -> M0 group), we extracted PPI data from the STRING database. The network graph highlights the pivotal regulatory functions of some of AP-1 family members (Jun, Junb, Jund, Fos, Fosl1, Fosl2, Maf, Atf3, Atf4, and Batf3) and Ctcf within this interaction network (Figure 6A). Notably, we also identified a cluster of ISGs (Ifi47, Ifit3, Osl1, Osl2, Oas1, Oas2, Oas3) in this network. These genes exhibited hysteresis in terms of promoter region accessibility (Figure 6B) and RNA expression levels (Figure 6C).

Figure 6

Figure 6 ChIP-seq analysis revealed the role of AP-1 in hysteresis. (A) Protein–protein interaction of hysteresis genes and AP-1/CTCF. (B) ATAC signal near promoter regions of ISGs across M0, M1, and deM0_M1. (C) Alterations in the expression level of ISGs from 0 h to 96 h. (D) The binding of Junb to the promoter regions of overall genes and hysteresis genes in BMDM changes upon LPS stimulation. (E) The binding of Junb to the promoter regions of overall genes and hysteresis genes in BMDM changes upon 4 h LPS and 4 h LPS+IL4 treatment. (F) H3k27ac signal to the promoter regions of overall genes and hysteresis genes in BMDM changes upon 4 h/24 h LPS treatment and 4 h IFN-γ+Junb knockdown.

We compiled ChIP-seq datasets from BMDMs, encompassing Junb binding profiles from 0 to 2 h post LPS stimulation, Junb binding at 4 h post LPS stimulation, Junb binding at 4 h following combined LPS and IL-4 stimulation, H3K27ac signals in BMDMs at 4 h and 24 h post LPS stimulation, and H3K27ac signals in Junb knockdown BMDMs at 4 h post IFN-γ. Through rigorous ChIP-seq data analysis, we established that while Junb did not exhibit an evident enrichment tendency within the promoter regions of hysteresis genes, it did show an amplified binding signal in the promoter regions of hysteresis genes relative to the binding signals across all gene promoters after stimulation with LPS (Figure 6D). Interestingly, further validation using another set of ChIP-seq data revealed that the binding trend of Junb within the promoter regions of hysteresis genes after LPS stimulation remained unaltered by IL-4 (an M2 polarization cue) stimulation (Figure 6E). Moreover, during LPS stimulation, a notable enhancement in the H3K27ac signal was observed within the promoter region of hysteresis genes. Notably, in macrophages with Junb knockdown, the augmented H3K27ac intensity in the promoter region of hysteresis genes was no longer evident upon another type of M1 stimulation (IFN-γ) (Figure 6F). These findings collectively underscore the robust regulatory influence of Junb on hysteresis genes during M1 stimulation, which cannot be fully mediated by M2 stimulation. This dynamic interplay contributes to the emergence of the hysteresis phenomenon.

4 Discussion

We sought to understand the molecular characteristics underlying macrophage de-/repolarization by integrating and re-analyzing large-scale publicly available datasets employing computational approaches. Liu et al. (9) demonstrated no significant global alteration in the transcriptome and chromatin accessibility of M0, M1, or M2 macrophages following de-/repolarization. Our findings reaffirmed this observation. However, our investigation specifically concentrated on a subset of genes with limited, yet notable, differential expressed after de-/repolarization (designated hysteresis genes), as well as regions exhibiting substantial differences in chromatin accessibility. Our observations encompass the following facets: 1. A number of genes maintained their expression patterns during polarized state (persistently high or persistently low expression). This phenomenon is termed “macrophage hysteresis”, and the genes manifesting this behavior are categorized as hysteresis genes. Notably, the persistence of M1 polarization gave rise to a more pronounced form of hysteresis than that induced by M2 polarization, reflected in the differences in RNA expression and chromatin accessibility. 2. These hysteresis genes exhibited expression trends analogous to those of conventional macrophage marker genes and played a pivotal role in macrophage phenotype classification. 3. The functional attributes of these hysteresis genes were interconnected with fundamental immunological and basic cell functions, such as inflammatory responses, cell cycle, and cell migration within macrophages. 4. In macrophages undergoing M0 -> M1 -> M0 polarization, TFBS of IRFs and Stat1 were enriched in the promoter regions of hysteresis genes. Conversely, AP-1 and CTCF bounded to mutually distinct distal regions (the COO and OCC regions, respectively). Collectively, these factors orchestrate regulating macrophage hysteresis (Figure 7) and trigger the cellular functional heterogeneity of macrophages.

Figure 7

Figure 7 Potential regulatory mechanisms underlying M1 hysteresis. Upon withdrawal of the M1 stimulus (LPS & IFN-g) with wash buffer in polarized M1 macrophages, IRFs and STAT exhibit binding affinity towards the promoter regions of hysteresis genes, whereas AP-1 and CTCF engage with the COO (Close-Open-Open) and OCC (Open-Close-Close) regions of chromatin, respectively. The conjoint orchestration of regulatory elements spanning both promoter and distal regions underpins the sustained expression profile of hysteresis genes. This transcriptional divergence subsequently reverberates across functionalities encompassing cell migration, antiviral response, antigen presentation, lipid metabolism, etc.

Given the opposing enrichment patterns, it is conceivable that AP-1 and CTCF may exert antagonistic effects on macrophage hysteresis by competing for binding sites, thereby inducing extensive alterations in the regulatory landscape. Although the relationship between AP-1 and CTCF remains unclear, in the human A549 cell line, CTCF depletion occurred at the hotspots enriched with AP-1 motifs (49). Within the context of inflammatory responses, the TLR4-NOS1-AP1 signaling axis governs macrophage polarization by facilitating AP-1 interactions and inducing IL-12/IL-23 (50). The AP-1 family including Fos and Jun (51) can co-modulate with NF-κB to regulate macrophage inflammatory responses (52). Furthermore, CTCF depletion attenuates acute inflammatory responses in induced macrophages (53). CTCF also controls proinflammatory gene expression through TLR signaling in macrophages (54). In addition to macrophages, through the regulation of the 3D enhancer network, CTCF participates in the inflammatory response of bone marrow-derived dendritic cells. Thus, both AP-1 and CTCF regulate inflammatory responses in cells. We identified a potential antagonistic relationship between AP-1 and CTCF; however, the precise mechanistic underpinnings require further investigation.

In our study, inducible ISGs exhibited hysteresis phenomena in response to M1 stimulation and displayed a robust correlation with AP-1 and CTCF within the PPI network of hysteresis genes. The downregulation of these type-I IFN genes serves as a protective mechanism against excessive lipoprotein uptake in human and murine foamy macrophages (55). Notably, another significant hysteresis gene, Ccl5, demonstrated tight PPIs with these ISGs. Elevated Ccl5 expression has been linked to epigenetic modifications induced by AP-1 in adipose tissue macrophages (ATMs) of mice fed a high-fat diet (56). Furthermore, in the context of obesity, S100a8/S100a9, which are significant hysteresis genes recognized as damage-associated molecular patterns (DAMPs), are presumed to trigger ATM-mediated immune responses (57–59). Interestingly, our findings suggested that the hysteresis of S100a8/S100a9 is controlled by the accessibility of the promoter and enhancer regions. This evidence suggests that, apart from the aforementioned cellular and immune functions, the differential expression of hysteresis genes is potentially linked to macrophage lipid metabolism. This intriguing avenue is a potential direction for future in-depth studies.

Polarization is thought to have little impact on macrophage memory, as traditional macrophage marker genes showed minimal differences in transcriptome and chromatin accessibility (9). However, we found that even subtle differences in expression at the cellular level can lead to functional changes in macrophages. Functional changes are likely to be mediated through epigenetic reprogramming, which is critical for adaptive immunity. For instance, a history of diet-induced obesity permanently alters chromatin accessibility (56).

However, several important points need to be considered. Firstly, our study has been built upon the analysis of data from multiple distinct projects, and as such, we were unable to completely control the batch effect due to various experimental conditions (dose of M1/M2 stimulus, etc.) and ensure that all macrophages can be fully de-/re-polarized within 96 hours. Secondly, given the limited number of differentially accessible peaks identified, we were unable to analyze the potential mechanisms underlying M2 hysteresis in the same manner as we did for M1 hysteresis. Furthermore, we could only consider the 96-hour as the end point of macrophage de-/repolarization, the hysteresis beyond the 96-hour stimulation period has not been discussed. Lastly, all the data were derived from mouse BMDM in vitro conditions. Therefore, the results may not directly correspond to tissue-resident macrophages in vivo conditions.

Taken together, our approach allows the establishment of a multi-omics profile of macrophage differentiation and promotes our understanding of cellular hysteresis, which contributes to the future development of novel macrophage repolarization therapies. Hysteresis has been hypothesized in other immune-related cells, and further efforts are needed to understand the general and cell-type-specific hysteresis mechanisms.

Data availability statement

Publicly available datasets were analyzed in this study. This data can be found here: GSE158094(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE158094) GSE117176(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE117176) GSE161125(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE161125) GSM1000074(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM1000074) GSE38377(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE38377) GSE36099(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE36099) GSE36104(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE36104) GSE84519(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE84519) GSE84520(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE84520) GSE178526(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE178526.

Author contributions

YZ: Conceptualization, Data curation, Formal Analysis, Methodology, Validation, Visualization, Writing – original draft, Writing – review & editing. WY: Formal Analysis, Visualization, Writing – review & editing. YK: Investigation, Writing – review & editing. ML: Formal Analysis, Writing – review & editing. WZ: Formal Analysis, Visualization. S-JP: Investigation, Writing – review & editing. KN: Conceptualization, Funding acquisition, Investigation, Project administration, Supervision, Writing – review & editing.

Funding

The author(s) declare financial support was received for the research, authorship, and/or publication of this article. The expenses for publication were covered by the annual running budget of the KN laboratory.

Acknowledgments

Computational resources were provided by the supercomputer system SHIROKANE at the Human Genome Center, Institute of Medical Science, the University of Tokyo.

Conflict of interest

The authors declare that this study was conducted in the absence of any commercial or financial relationships that could be construed as potential conflicts of interest.

The author(s) declared that they were an editorial board member of Frontiers, at the time of submission. This had no impact on the peer review process and the final decision.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fimmu.2023.1304778/full#supplementary-material

References

1. Epelman S, Lavine KJ, Randolph GJ. Origin and functions of tissue macrophages. Immunity (2014) 41(1):21–35. doi: 10.1016/j.immuni.2014.06.013

PubMed Abstract | CrossRef Full Text | Google Scholar

2. Murray RZ, Stow JL. Cytokine secretion in macrophages: SNAREs, rabs, and membrane trafficking. Front Immunol (2014) 5:538. doi: 10.3389/fimmu.2014.00538

PubMed Abstract | CrossRef Full Text | Google Scholar

3. Wynn TA, Chawla A, Pollard JW. Macrophage biology in development, homeostasis, and disease. Nature (2013) 496(7446):445–55. doi: 10.1038/nature12034

PubMed Abstract | CrossRef Full Text | Google Scholar

4. Murray PJ. Macrophage polarization. Annu Rev Physiol (2017) 79:541–66. doi: 10.1146/annurev-physiol-022516-034339

PubMed Abstract | CrossRef Full Text | Google Scholar

5. Rolot M, Dewals BG. Macrophage activation and functions during helminth infection: Recent advances from the laboratory mouse. J Immunol Res (2018) 2018:2790627. doi: 10.1155/2018/2790627

PubMed Abstract | CrossRef Full Text | Google Scholar

6. Kadomoto S, Izumi K, Mizokami A. Macrophage polarity and disease control. Int J Mol Sci (2021) 23(1):144. doi: 10.3390/ijms23010144

PubMed Abstract | CrossRef Full Text | Google Scholar

7. Orecchioni M, Ghosheh Y, Pramod AB, Ley K. Macrophage Polarization: Different Gene Signatures in M1(LPS+) vs. classically and M2(LPS–) vs. Alternatively Activated Macrophages. Front Immunol (2019) 10:1084. doi: 10.3389/fimmu.2019.01084

PubMed Abstract | CrossRef Full Text | Google Scholar

8. Degboé Y, Poupot R, Poupot M. Repolarization of unbalanced macrophages: Unmet medical need in chronic inflammation and cancer. Int J Mol Sci (2022) 23(3):1496. doi: 10.3390/ijms23031496

PubMed Abstract | CrossRef Full Text | Google Scholar

9. Liu SX, Gustafson HH, Jackson DL, Pun SH, Trapnell C. Trajectory analysis quantifies transcriptional plasticity during macrophage polarization. Sci Rep (2020) 10(1):12273. doi: 10.1038/s41598-020-68766-w

PubMed Abstract | CrossRef Full Text | Google Scholar

10. Smith TD, Tse MJ, Read EL, Liu WF. Regulation of macrophage polarization and plasticity by complex activation signals. Integr Biol (Camb) (2016) 8(9):946–55. doi: 10.1039/c6ib00105j

PubMed Abstract | CrossRef Full Text | Google Scholar

11. Veremeyko T, Yung AWY, Anthony DC, Strekalova T, Ponomarev ED. Early growth response Gene-2 is essential for M1 and M2 macrophage activation and plasticity by modulation of the transcription factor CEBPβ. Front Immunol (2018) 9:2515. doi: 10.3389/fimmu.2018.02515

PubMed Abstract | CrossRef Full Text | Google Scholar

12. Netea MG, Joosten LAB, Latz E, Mills KHG, Natoli G, Stunnenberg HG, et al. Trained immunity: A program of innate immune memory in health and disease. Science (2016) 352(6284):aaf1098. doi: 10.1126/science.aaf1098

PubMed Abstract | CrossRef Full Text | Google Scholar

13. Daniel B, Belk JA, Meier SL, Chen AY, Czimmerer Z, et al. Macrophage inflammatory and regenerative response periodicity is programmed by cell cycle and chromatin state. Mol Cell (2023) 83(1):121–138.e7. doi: 10.1016/j.molcel.2022.11.017

PubMed Abstract | CrossRef Full Text | Google Scholar

14. Muñoz-Rojas AR, Kelsey I, Pappalardo JL, Chen M, Miller-Jensen K. Co-stimulation with opposing macrophage polarization cues leads to orthogonal secretion programs in individual cells. Nat Commun (2021) 12(1):301. doi: 10.1038/s41467-020-20540-2

PubMed Abstract | CrossRef Full Text | Google Scholar

15. Li C, Menoret A, Farragher C, Ouyang Z, Bonin C, Holvoet P, et al. Single cell transcriptomics based-MacSpectrum reveals novel macrophage activation signatures in diseases. JCI Insight (2019) 5(10):e126453. doi: 10.1172/jci.insight.126453

PubMed Abstract | CrossRef Full Text | Google Scholar

16. Trapnell C, Hendrickson DG, Sauvageau M, Goff L, Rinn JL, Pachter L, et al. Differential analysis of gene regulation at transcript resolution with RNA-seq. Nat Biotechnol (2013) 31(1):46–53. doi: 10.1038/nbt.2450

PubMed Abstract | CrossRef Full Text | Google Scholar

17. Stuart T, Butler A, Hoffman P, Hafemeister C, Papalexi E, Mauck WM, et al. Comprehensive integration of single-cell data. Cell (2019) 177(7):1888–1902.e21. doi: 10.1016/j.cell.2019.05.031

PubMed Abstract | CrossRef Full Text | Google Scholar

18. Loza M, Teraguchi S, Standley DM, Diez D. Unbiased integration of single cell transcriptome replicates. NAR Genom Bioinform (2022) 4(1):lqac022. doi: 10.1093/nargab/lqac022

PubMed Abstract | CrossRef Full Text | Google Scholar

19. FASTQC. A quality control tool for high throughput sequence data . BibSonomy. Available at: https://www.bibsonomy.org/bibtex/f230a919c34360709aa298734d63dca3 (Accessed 5.28.23).

Google Scholar

20. Martin M. Cutadapt removes adapter sequences from high-throughput sequencing reads. EMBnet.journal (2011) 17(1):10–2. doi: 10.14806/ej.17.1.200

CrossRef Full Text | Google Scholar

21. Langmead B, Salzberg SL. Fast gapped-read alignment with Bowtie 2. Nat Methods (2012) 9(4):357–9. doi: 10.1038/nmeth.1923

PubMed Abstract | CrossRef Full Text | Google Scholar

22. Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. 1000 Genome project data processing subgroup. The sequence alignment/Map format and SAMtools. Bioinformatics (2009) 25(16):2078–9. doi: 10.1093/bioinformatics/btp352

PubMed Abstract | CrossRef Full Text | Google Scholar

23. Tarasov A, Vilella AJ, Cuppen E, Nijman IJ, Prins P. Sambamba: Fast processing of NGS alignment formats. Bioinformatics (2015) 31(12):2032–4. doi: 10.1093/bioinformatics/btv098

PubMed Abstract | CrossRef Full Text | Google Scholar

24. Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, et al. Model-based analysis of chip-seq (MACS). Genome Biol (2008) 9(9):R137. doi: 10.1186/gb-2008-9-9-r137

PubMed Abstract | CrossRef Full Text | Google Scholar

25. Ramírez F, Ryan DP, Grüning B, Bhardwaj V, Kilpert F, Richter AS, et al. deepTools2: A next generation web server for deep-sequencing data analysis. Nucleic Acids Res (2016) 44(W1):W160–5. doi: 10.1093/nar/gkw257

PubMed Abstract | CrossRef Full Text | Google Scholar

26. Thorvaldsdóttir H, Robinson JT, Mesirov JP. Integrative Genomics Viewer (IGV): High-performance genomics data visualization and exploration. Brief Bioinform (2013) 14(2):178–92. doi: 10.1093/bib/bbs017

PubMed Abstract | CrossRef Full Text | Google Scholar

27. Stuart T, Srivastava A, Madad S, Lareau CA, Satija R, et al. Single-cell chromatin state analysis with Signac. Nat Methods (2021) 18(11):1333–41. doi: 10.1038/s41592-021-01282-5

PubMed Abstract | CrossRef Full Text | Google Scholar

28. Pedregosa F, Varoquaux G, Gramfort A, Michel V, Thirion B, Grisel O, et al. Scikit-learn: Machine learning in python. J Mach Learn Res (2011) 12:2825–30.

Google Scholar

29. Obayashi T, Kodate S, Hibara H, Kagaya Y, Kinoshita K. COXPRESdb v8: An animal gene coexpression database navigating from a global view to detailed investigations. Nucleic Acids Res (2023) 51(D1):D80–7. doi: 10.1093/nar/gkac983

PubMed Abstract | CrossRef Full Text | Google Scholar

30. Han H, Cho JW, Lee S, Yun A, Kim H, Bae D, et al. TRRUST v2: An expanded reference database of human and mouse transcriptional regulatory interactions. Nucleic Acids Res (2018) 46(D1):D380–6. doi: 10.1093/nar/gkx1013

PubMed Abstract | CrossRef Full Text | Google Scholar

31. Perozzi B, Al-Rfou R, Skiena S. DeepWalk: Online learning of social representations in Proceedings of the 20th ACM SIGKDD International Conference on Knowledge Discovery and Data Mining, KDD’14. New York NY: Assoc Computing Machinery (2014). doi: 10.1145/2623330.2623732

CrossRef Full Text | Google Scholar

32. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data. Innovation (Camb) (2021) 2(3):100141. doi: 10.1016/j.xinn.2021.100141

PubMed Abstract | CrossRef Full Text | Google Scholar

33. Puig RR, Boddie P, Khan A, Castro-Mondragon JA, Mathelier A. UniBind: Maps of high-confidence direct TF-DNA interactions across nine species. BMC Genomics (2021) 22(1):482. doi: 10.1186/s12864-021-07760-6

PubMed Abstract | CrossRef Full Text | Google Scholar

34. Vandenbon A, Dinh VH, Mikami N, Kitagawa Y, Teraguchi S, Ohkura N, et al. Immuno-Navigator, a batch-corrected coexpression database, reveals cell type-specific gene networks in the immune system. Proc Natl Acad Sci U.S.A. (2016) 113(17):E2393–402. doi: 10.1073/pnas.1604351113

PubMed Abstract | CrossRef Full Text | Google Scholar

35. Fulco CP, Nasser J, Jones TR, Munson G, Bergman DT, Subramanian V, et al. Activity-by-contact model of enhancer-promoter regulation from thousands of CRISPR perturbations. Nat Genet (2019) 51(12):1664–9. doi: 10.1038/s41588-019-0538-0

PubMed Abstract | CrossRef Full Text | Google Scholar

36. Quinlan AR. BEDTools: The Swiss-army tool for genome feature analysis. Curr Protoc Bioinf (2014) 47:11.12.1–11.1234. doi: 10.1002/0471250953.bi1112s47

CrossRef Full Text | Google Scholar

37. Szklarczyk D, Gable AL, Nastou KC, Lyon D, Kirsch R, Pyysalo S, et al. The STRING database in 2021: Customizable protein–protein networks, and functional characterization of user-uploaded gene/measurement sets. Nucleic Acids Res (2021) 49(D1):D605–12. doi: 10.1093/nar/gkaa1074

PubMed Abstract | CrossRef Full Text | Google Scholar

38. Sierra-Filardi E, Nieto C, Domínguez-Soto A, Barroso R, Sánchez-Mateos P, Puig-Kroger A, et al. CCL2 shapes macrophage polarization by GM-CSF and M-CSF: Identification of CCL2/CCR2-dependent gene expression profile. J Immunol (2014) 192(8):3858–67. doi: 10.4049/jimmunol.1302821

PubMed Abstract | CrossRef Full Text | Google Scholar

39. Wang S, Song R, Wang Z, Jing Z, Wang S, Ma J. S100A8/A9 in inflammation. Front Immunol (2018) 9:1298. doi: 10.3389/fimmu.2018.01298

PubMed Abstract | CrossRef Full Text | Google Scholar

40. Franz S, Ertel A, Engel KM, Simon JC, Saalbach A. Overexpression of S100A9 in obesity impairs macrophage differentiation via TLR4-NFkB-signaling worsening inflammation and wound healing. Theranostics (2022) 12(4):1659–82. doi: 10.7150/thno.67174

PubMed Abstract | CrossRef Full Text | Google Scholar

41. Jimenez-Preitner M, Berney X, Uldry M, Vitali A, Cinti S, Ledford JG, et al. Plac8 is an inducer of C/EBPβ required for Brown fat differentiation, thermoregulation, and control of body weight. Cell Metab (2011) 14(5):658–70. doi: 10.1016/j.cmet.2011.08.008

PubMed Abstract | CrossRef Full Text | Google Scholar

42. Günthner R, Anders HJ. Interferon-regulatory factors determine macrophage phenotype polarization. Mediators Inflammation (2013) 2013:731023. doi: 10.1155/2013/731023

CrossRef Full Text | Google Scholar

43. Ivashkiv LB, Hu X. Signaling by STATs. Arthritis Res Ther (2004) 6(4):159–68. doi: 10.1186/ar1197

PubMed Abstract | CrossRef Full Text | Google Scholar

44. Rigamonti E, Chinetti-Gbaguidi G, Staels B. Regulation of macrophage functions by PPAR-α, PPAR-γ, and LXRs in mice and men. Arterioscler Thromb Vasc Biol (2008) 28(6):1050–9. doi: 10.1161/ATVBAHA.107.158998

PubMed Abstract | CrossRef Full Text | Google Scholar

45. Chu YB, Li J, Jia P, Cui J, Zhang R, Kang X, et al. Irf1- and Egr1-activated transcription plays a key role in macrophage polarization: A multiomics sequencing study with partial validation. Int Immunopharmacol (2021) 99:108072. doi: 10.1016/j.intimp.2021.108072

PubMed Abstract | CrossRef Full Text | Google Scholar

46. Cui H, Banerjee S, Guo S, Xie N, Liu G. Interferon regulatory factor 2 inhibits expression of glycolytic genes and lipopolysaccharide induced pro-inflammatory responses in macrophages. J Immunol (2018) 200(9):3218–30. doi: 10.4049/jimmunol.1701571

PubMed Abstract | CrossRef Full Text | Google Scholar

47. Lawrence T, Natoli G. Transcriptional regulation of macrophage polarization: Enabling diversity with identity. Nat Rev Immunol (2011) 11(11):750–61. doi: 10.1038/nri3088

PubMed Abstract | CrossRef Full Text | Google Scholar

48. Project Consortium ENCODE, Moore JE, Purcaro MJ, Pratt HE, Epstein CB, Shoresh N, et al. Expanded encyclopaedias of DNA elements in the human and mouse genomes. Nature (2020) 583(7818):699–710. doi: 10.1038/s41586-020-2493-4

PubMed Abstract | CrossRef Full Text | Google Scholar

49. Seo J, Koçak DD, Bartelt LC, Williams CA, Barrera A, Gersbach CA, et al. AP-1 subunits converge promiscuously at enhancers to potentiate transcription. Genome Res (2021) 31(4):538–50. doi: 10.1101/gr.267898.120

PubMed Abstract | CrossRef Full Text | Google Scholar

50. Srivastava M, Saqib U, Naim A, Roy A, Liu D, Bhatnagar D, et al. The TLR4-NOS1-AP1 signaling axis regulates macrophage polarization. Inflammation Res (2017) 66(4):323–34. doi: 10.1007/s00011-016-1017-z

CrossRef Full Text | Google Scholar

51. Fujioka S, Niu J, Schmidt C, Sclabas GM, Peng B, Uwagawa T, et al. NF-κB and AP-1 connection: Mechanism of NF-κB-dependent regulation of AP-1 activity. Mol Cell Biol (2004) 24(17):7806–19. doi: 10.1128/MCB.24.17.7806-7819.2004

PubMed Abstract | CrossRef Full Text | Google Scholar

52. Dorrington MG, Fraser IDC. NF-κB signaling in macrophages: Dynamics, crosstalk, and signal integration. Front Immunol (2019) 10:705. doi: 10.3389/fimmu.2019.00705

PubMed Abstract | CrossRef Full Text | Google Scholar

53. Stik G, Vidal E, Barrero M, Cuartero S, Vila-Casadesús M, Mendieta-Esteban J, et al. CTCF is dispensable for immune cell transdifferentiation but facilitates an acute inflammatory response. Nat Genet (2020) 52(7):655–61. doi: 10.1038/s41588-020-0643-0

PubMed Abstract | CrossRef Full Text | Google Scholar

54. Nikolic T, Movita D, Lambers ME, Ribeiro de Almeida CR, Biesta P, Kreefft K, et al. The DNA-binding factor CTCF critically controls gene expression in macrophages. Cell Mol Immunol (2014) 11(1):58–70. doi: 10.1038/cmi.2013.41

PubMed Abstract | CrossRef Full Text | Google Scholar

55. Willemsen L, Chen HJ, van Roomen CPAA, Griffith GR, Siebeler R, Neele AE, et al. Monocyte and macrophage lipid accumulation results in down-regulated Type-I interferon responses. Front Cardiovasc Med (2022) 9:829877. doi: 10.3389/fcvm.2022.829877

PubMed Abstract | CrossRef Full Text | Google Scholar

56. Hata M, Andriessen EMMA, Hata M, Diaz-Marin R, Fournier F, Crespo-Garcia S, et al. Past history of obesity triggers persistent epigenetic changes in innate immunity and exacerbates neuroinflammation. Science (2023) 379(6627):45–62. doi: 10.1126/science.abj8894

PubMed Abstract | CrossRef Full Text | Google Scholar

57. Erridge C. Endogenous ligands of TLR2 and TLR4: Agonists or assistants? J Leukoc Biol (2010) 87(6):989–99. doi: 10.1189/jlb.1209775

PubMed Abstract | CrossRef Full Text | Google Scholar

58. Vandanmagsar B, Youm YH, Ravussin A, Galgani JE, Stadler K, Mynatt RL, et al. The NLRP3 inflammasome instigates obesity-induced inflammation and insulin resistance. Nat Med (2011) 17(2):179–88. doi: 10.1038/nm.2279

PubMed Abstract | CrossRef Full Text | Google Scholar

59. Wen H, Gris D, Lei Y, Jha S, Zhang L, Huang MT-H, et al. Fatty acid–induced NLRP3-ASC inflammasome activation interferes with insulin signaling. Nat Immunol (2011) 12(5):408–15. doi: 10.1038/ni.2022

PubMed Abstract | CrossRef Full Text | Google Scholar