Extensive Computational Screening Identifies Marco (SR-A6) and Abca1/2 as Key Regulators of Intracellular Lipid Metabolism to Compensate Low Cholesterol Biosynthesis in M1 Macrophages Hysteresis

1. Introduction

In the human body, macrophages perform fundamental functions such as maintaining hemostasis and resisting pathogen invasion [1]. Upon exposure to various environmental stimuli, macrophages polarize into distinct phenotypes within different tissues, including M1 macrophages and M2 macrophages [2]. For example, lipopolysaccharide (LPS) can drive polarization toward the M1 phenotype, while interleukin-4 (IL-4) promotes the polarization of the M2 subtype [3]. M1 macrophages are primarily responsible for triggering inflammatory responses, whereas M2 macrophages serve to balance inflammation and promote wound healing and tissue repair [4]. In inflammatory tissues, macrophages often undergo M1 polarization initially and are subsequently repolarized into the M2 phenotype for tissue repair. This transition from one phenotype to another is referred to as repolarization or reprogramming.

Numerous studies have detailed memory T and B cells[5,6,7], yet the concept of innate memory and specifically memory macrophages was only established about a decade ago. Recent research challenges the traditional view that relegates innate immune cells merely to the first line of defense, suggesting instead that these cells can also develop and maintain immunologic memory or hysteresis[8,9]. Multiple independent clinical and biological studies have confirmed such macrophage hysteresis[10,11,12,13,14]. This de-facto innate immune memory, referred to as trained immunity, is mediated through extensive metabolic rewiring and epigenetic modifications[15]. Energy and lipid metabolism, along with dietary compounds, influence enzymes that regulate chromatin compaction and structure, whereby lipid metabolites play a critical role in mediating epigenetic modifications that impact macrophage immune memory, contributing significantly to trained immunity[16]. Multiple experiential evidence underscore the significant influence of lipid metabolism on innate immune memory and macrophage hysteresis[17,18].

In our previous research, we reported significant hysteresis in bone marrow-derived macrophages following M1 stimulation and subsequent de-/repolarization, attributed to the regulation of various chromatin sites by factors such as AP-1 and CTCF. This sustained M1 hysteresis induces functional heterogeneity in macrophages, affecting inflammatory responses, cell cycle regulation, and cell migration. Remarkably, we observed a hysteresis phenomenon in a set of interferon-stimulated genes (ISGs)[9]. Research by Lisa Willemsen et al. indicates that lipid accumulation in mouse and human macrophages leads to differential expression of type-I interferons[19]. Furthermore, studies suggest that lipid metabolism is crucial for mounting effective inflammatory responses[20]. We hypothesize that the observed M1 hysteresis may be directly linked to macrophage lipid metabolism. In this study, through RNA-seq analysis of a total of 162 polarized and lipid-loaded macrophages (Figure 1), we demonstrate that the scavenger receptor Marco (SR-A6) sustains high expression during de-/repolarization, leading to persistent lipid intake. Concurrently, the low expression of ATP Binding Cassette Subfamily A members (Abca1 and Abca2) results in reduced lipid efflux, collectively creating an environment of high intracellular lipid levels. This environment compensates for the suppressed expression of lipid and cholesterol biosynthesis-related genes. The elevated lipid levels provide sustained energy to macrophages, prolonging the inflammatory state associated with M1 hysteresis. These findings significantly enhance our understanding of the regulation of macrophage hysteresis.

2. Results

2.1. Subsection

2.1.1. Trajectory Analysis and Hierarchical Clustering Confirmed the High Sensitivity of Macrophages to the Time of Stimulation

We initially categorized the processed RNA-seq data from macrophage polarization and lipid-loaded macrophages according to different stimulation conditions and obtained eight distinct groups: lipid (all lipid-loaded macrophages described in method section), M0 (unpolarized macrophages), M1 (M1 macrophages polarized from M0), M2 (M2 macrophages polarized from M0), deM0_M1 (depolarized M0 from M1), deM0_M2 (depolarized M0 from M2), reM2_M1 (repolarized M2 from M1), and reM1_M2 (repolarized M1 from M2). Then, we performed batch effect correction on the TPM-transformed RNA expression matrix containing 162 samples. We compared the deviation of expression levels before and after batch effect correction (Figure 2A). We observed that after batch effect correction using mouse housekeeping genes, both intra-group and inter-group expression variations were significantly reduced across all samples, indicating the successful elimination of batch effects. To further describe the distribution of all samples and validate sample quality, we compared the PCA plot generated from the original TPM matrix (Figure S1) with the PCA plot after batch removal (Figure 2B). In the PCA plot after batch removal, the trajectory of M1 polarization appeared longer than that of M2 polarization, consistent with previous findings by Liu et al[21].

However, the lipid-loaded macrophages, highlighted with a red dashed-line box, exhibited differential expression pattern compared to all samples across the macrophage polarization trajectory, although we observed that lipid-loaded macrophages were closer to the M1 or M1-like phenotypes (M1 and reM1_M2). This difference is not attributable to batch effects and is likely to indirectly affect the selection of gene clusters in subsequent WGCNA analysis, as the overall gene expression levels of lipid-loaded samples may exhibit a relatively low correlation with those of other polarized phenotypes. We then performed a preliminary WGCNA experiment using lipid samples as a phenotype (Figure S2). The results indicated that the lipid phenotype exhibited very low and non-significant correlations with all ten gene clusters compared to other phenotypes in Figure S2. Therefore, we considered not using the lipid samples directly in the WGCNA analysis but rather employing them in subsequent analyses to validate the significantly expressed lipid metabolism-related genes discovered in polarized samples.

To further ascertain the distribution among samples, we employed the hierarchical clustering method (Figure 2C). The results confirm that lipid samples marked with a red dashed-line box, clustered together and were distinctly different from other polarized samples. Interestingly, the results indicate that the distribution of all samples does not center from the M0-like phenotype and extends towards M1 and M2 directions, underscoring that in addition to the type of stimulation, the duration of stimulation also significantly influences the characteristics of macrophages (Figure 2C). Based on the analysis above, we have discovered that: 1. Lipid-loaded samples are not suitable to be directly analyzed as a distinct phenotype alongside other polarized samples, as this could significantly impact the analysis results. 2. Due to the high plasticity of macrophages and their extreme sensitivity to the duration of stimulation, polarized macrophages under different stimulation durations require further classification.

2.1.2. Considering the Stimulation Time Factor Significantly Improved the Accuracy of Sample Classification for Downstream Analysis

To accurately classify all macrophage samples and enhance the significance of subsequent analyses, we incorporated annotations of stimulation time to the hierarchical clustering tree in Figure 2C (Figure 3). We used the positions of M0, M1, and M2 samples in the hierarchical clustering tree as critical points, which indicate when macrophages are about to change their original phenotypes. The results show that macrophages polarized to M1 from 0 to 4 hours exhibit significant differential expression compared to those polarized for extended periods (4 to 24 hours). This explains why, in the previous classification, these M1-polarized macrophages were categorized under the M2-like group. Similarly, M2 macrophages undergoing M1 repolarization for less than 2 hours (reM1_M2, 24–26h) also show substantial differential expression compared to cells stimulated for longer periods (reM1_M2, 28–96h). Evidently, this is because these M0 and M2 cells have not completed polarization or repolarization within the short time frame and categorizing these cells as either M1-like or M2-like based solely on polarization type without considering the extent of polarization would severely impact subsequent analyses. Additionally, we observed that the depolarized M1 cells at 96 hours exhibit a state more similar to the M0 phenotype rather than the M1 phenotype, which led us to differentiate these samples from the depolarized M1 cells stimulated for 24 to 48 hours.

Based on our analyses and observations, we subsequently re-categorized all samples into the following 9 groups for downstream analysis to minimize the impact caused by incomplete cellular polarization or short period of stimulation: M0, 0~4h M1 (0_4M1), 4~96h M1 (4_96M1), M2, 24~26h repolarized M1 from M2 (24_26reM1_M2), 28 ~ 96h repolarized M1 from M2 (28_96reM1_M2), 24~48h depolarized M0 from M1 (24_48deM0_M1), 96h depolarized M0 from M1 (96deM0_M1), depolarized M0 from M2 (deM0_M2), and repolarized M2 from M1 (reM2_M1) (Table S1).

2.1.3. WGCNA Revealed Specific Gene Modules Associated with Macrophage Phenotypes

Continuing, using the reclassified RNA expression data, we utilized WGCNA to analyze the relationship between the obtained gene expression matrix and the nine macrophage polarization phenotypes. We chose a soft thresholding power of 6, and the scale-free fit index (signed R^2) was 0.9 (Figure 4A) and adjusted the minimum module size (minModuleSize) to 30 genes, then merged modules with a height cut-off (mergeCutHeight) less than 0.25. In this manner, we selected ten independent co-expression modules based on the generated hierarchical clustering dendrogram (Figure 4B), specifically named MEblack (51 genes), MEblue (1084 genes), MEbrown (466 genes), MEgreen (144 genes), MEpink (39 genes), MEred (69 genes), MEturquoise (1770 genes), MEyellow (203 genes), as significant gene modules for downstream analysis. MEgrey with 1174 genes represents unmapped genes in this context (Figure 4C). Through an investigation of these modules' eigengenes and eigengene associations, we found that despite their independence, MEgreen, MEturquoise, MEblack, and MEpink were relatively highly correlated. MEbrown, MEred, MEblue, and MEyellow were relatively highly correlated (Figure S3). Subsequently, we computed the correlation of each independent module with each macrophage phenotype (Figure 4C). To assess the correlation between gene clusters and M1/M2 hysteresis, we focused on phenotypes with depolarization or repolarization history. For instance, we observed that MEgreen shows a relatively high correlation with the deM0_M2 and M2 phenotypes, while it is not correlated with the M0 phenotype. This suggests that the genes in MEgreen retain the expression pattern of the M2 phenotype within depolarized M0 cells, indicating the presence of M2 hysteresis. Thus, MEgreen is associated with M2 hysteresis in M0 cells. Using a similar approach, we found that the clusters MEgreen, MEturquoise, MEblack, and MEpink exhibited strong correlations with the M0, M2, deM0_M2, and incompletely polarized phenotypes (0_4M1 and 24_26reM1_M2). These genes represent the features associated with the M0-like and M2-like phenotypes. We observed that MEgreen and MEturquoise are associated with early M1 repolarized (24_26reM1_M2) and M0 depolarized M2 (deM0_M2) hysteresis. However, as the stimulation duration increases, the correlation with M2 hysteresis diminishes following prolonged M1 repolarization, as evidenced by the low correlation observed between these two gene clusters and the late M1 repolarized macrophages (28_96reM1_M2). This aligns with our previous research findings that the intensity of M1 hysteresis is significantly larger than that of M2 hysteresis, and it is challenging to clearly define M2 hysteresis genes and chromatin regions [9].

On the other hand, the clusters MEbrown, MEred, MEblue, and MEyellow showed relatively high correlations with M1-like phenotypes (28_96reM1_M2, 4_96M1) and phenotypes displaying M1 hysteresis (24_48deM0_M1, reM2_M1). We specifically focused on the clusters MEbrown and MEblue, noting that MEbrown exhibits a stronger correlation with M1 hysteresis phenotypes (24_48deM0_M1 and reM2_M1) compared to MEblue. However, MEbrown does not show a correlation with the 4_96M1 and 28_96reM1_M2, which represent M1-like phenotypes. Upon comparing the gene expression profiles of these two gene clusters, we observed that, unlike the relatively consistent expression of genes in MEblue, most genes in MEbrown exhibit either high or low expression after more than 12 hours of M1 polarization and repolarization (Figure 4D). The initiation of polarization significantly impacts the RNA expression data from 0 to 12 hours, thereby affecting the results of the correlation analysis. This confirms that MEbrown also has a strong correlation with M1 cells and repolarized M1 cells stimulated for more than 12 hours. It again underscores the necessity of reclassifying macrophage phenotypes based on stimulation duration.

In summary, we identified two gene clusters, MEbrown and MEblue, that are significantly associated with M1 hysteresis. Notably, genes within the MEbrown cluster exhibit a late-onset phenomenon, distinguishing them from those in the MEblue cluster.

2.1.4. GO Term and KEGG Enrichment Analysis Unveiled the Functions of Identified Gene Modules

To further explore the functions of the identified gene modules by WGCNA, we performed GO term analysis (Figure 5A) and KEGG enrichment analysis (Figure 5B). The results revealed that all eight modules exhibited enrichments in specific GO terms and KEGG pathways. We observed that in the gene modules associated with M0 and M2, the top enriched biological processes and pathways for MEblack predominantly involve immune-related functions. These include defense response to the bacterium, regulation of inflammatory response, regulation of immune effector process, chemokine signaling pathway, and viral protein interaction with cytokine and cytokine receptor. MEpink is primarily associated with functions related to the extracellular matrix, including extracellular matrix and structure organization, connective tissue development, cell substrate adhesion, and focal adhesion. It is noteworthy that MEgreen and MEturquoise, which are associated with M2 hysteresis, are linked to leukocyte proliferation, chromosome segregation, nuclear division, DNA replication, and the cell cycle. This is consistent with our previous research findings, which concluded that M2 hysteresis potentially affects the cell cycle [9]. Subsequently, we examined gene modules highly associated with the M1-like phenotype. MEyellow represents hallmark M1 immunity and inflammatory responses, including, but not limited to, the NOD-like receptor signaling pathway, TNF signaling pathway, MAPK signaling pathway, NF-κB signaling pathway, and IL17 signaling pathway. MEred is associated with metabolic functions. Within the modules related to M1 hysteresis, MEblue is associated with M1 immune functions, including regulation of immune effector processes, response to viruses, and interferon-mediated signaling pathways.

Since the results of the functional analysis indicated that MEblue is significantly associated with interferon-mediated signaling pathways (with an adjusted p-value = 2.42E-17), we investigated the expression trends of all hysteresis genes within MEblue. We found that 13 ISGs (including Ccl5, Isg15, Ifi44, Ifi47, Ifi205, Ifit1, Ifit2, Ifit3, Ifitm1, Oasl1, Oasl2, Oas3, and Mx1) exhibiting hysteresis phenomena [9] in MEblue (Figure S4). Additionally, significant hysteresis genes including Ccr3 and Lcn2 were also identified within the MEblue module, indicating a strong correlation between MEblue and M1 hysteresis. MEbrown, which is another gene cluster that is significantly associated with M1 hysteresis, is enriched in pathways such as cholesterol and secondary alcohol biosynthesis, lipid metabolism, and atherosclerosis, as well as steroid biosynthesis. Additionally, four M1 hysteresis genes (Fn1, Timp1, Dcstamp, Serpinb2) were identified in MEbrown (Figure S5). These genes exhibit sustained non-expression or low expression following M1 stimulation. The above finding preliminarily confirms our hypothesis that M1 hysteresis in macrophages influences their lipid metabolism.

2.1.5. Lipid Biosynthesis-Related Genes Drive the Expression of MEbrowm to Influence Macrophage Function in Lipid Metabolism

The notable differences in gene expression confirm our findings in the functional analysis section, indicating a close association between MEblue and M1 immune functions and interferon-mediated signaling pathways. Willemsen, Lisa, et al. previously reported significant expression changes in these ISGs in lipid-loaded samples. The expression of ISGs in macrophages may be related to intracellular lipid levels. Based on our functional analysis, we speculate that regulatory factors present in the MEbrown cluster may influence macrophage lipid metabolism and, consequently, affect lipid levels within macrophages. To address this issue, we independently screened MEbrown genes for hub genes, which are highly connected to other genes in the gene regulatory network and are crucial for the overall regulation of the network, specifically in association with the 24_48h_deM0_M1 and reM2_M1 phenotypes. We utilized MM (module membership) > 0.7 and GS (gene significance) > 0.4 as filtering criteria for hub genes (Figure 6A). We found 117 genes and 42 genes as hub genes for the 24_48h_deM0_M1 and reM2_M1 phenotypes, respectively, within MEbrown. Through comparison, we identified 22 common hub genes between the described two phenotypes, including Hspa8, Tpm4, Crip1, Pros1, Fdft1, Idi1, Dhcr24, Gm3571, Gm5873, Cyp51, Fdps, Plaur, Gm43712, Slc11a1, Selenow, Tigd2, Ddit3, Ormdl3, Gadd45g, A530064D06Rik, Tnfrsf23, and Smg9 (Figure 6B). After calculating the Pearson correlation coefficients between the expression levels of these genes across all 144 samples, we found that these genes exhibit significantly high correlations with each other, distinctly forming two groups with either positive or negative correlations (Figure 6C). Then we investigated the interactions of these genes in the PPI network (Figure 4D). Of note, a sub-network including the proteins Dhcr24, Idi1, Fdps, Fdft1, Cyp51, Sc5d, and Por was identified, with four of these proteins being common M1 hysteresis hub genes that we screened for. Among them, Dhcr24 (24-dehydrocholesterol reductase), squalene synthase Fdft1 (farnesyl diphosphate farnesyltransferase 1), and Cyp51 (cholesterol biosynthetic enzyme lanosterol 14α-demethylase) all have been reported to have significant impacts on cholesterol biosynthesis [22]. Similarly, Idi1 has been confirmed to play a key role in the cholesterol biosynthesis pathway and lipid metabolism synthesis [23]. Fdps is another key cholesterol biosynthetic gene that has always been induced in colon cancer tissues [24]. Also, inhibition of POR activity would halt the production of cholesterol and its breakdown into bile acids [25]. Sc5D (sterol-C5-desaturase), which encodes an enzyme involved in cholesterol biosynthesis, converts lathosterol to 7-dehydrocholesterol in the cholesterol biosynthesis pathway [26]. The evidence suggest that this sub-network of hub genes plays a significant role in lipid metabolism and cholesterol biosynthesis.

Subsequently, we screened the co-expression data from the STRING database and found that, among all common hub genes, only the co-expression relationships involving Fdft1, Idi1, Fdps, Cyp51, and Dhcr24 have been consistently observed across numerous experiments in the database. Specifically, the co-expression of Cyp51 and Dhcr24 has been observed not only in Mus musculus but also in Homo sapiens, Bos taurus, and Gallus gallus base on the record from string database, highlighting the conservation co-expression relationship between Cyp51 and Dhcr24 across different species (Figure 6E). The roles of these hub genes in lipid metabolism were further confirmed in lipid-loaded macrophages which we collected from GEO database. By comparing to the RNA expression level of DMSO-treated BMDM, significant expression differences were observed during lipid loading processes (AcLDL, GW3965, and KLA) for Selenow, Slc11a1, Hspa8, Dhcr24, Idi1, Fdps, Fdft1, Plaur, Pros1, and Crip1 (Figure 6F). The expression of genes related to cholesterol biosynthesis, including Dhcr24, Idi1, Fdps, and Fdft1, was suppressed in lipid-loaded samples. This suggests that the accumulation of externally ingested lipids in macrophages significantly suppresses the expression of cholesterol biosynthesis-related genes, reducing endogenous lipid synthesis to maintain cellular lipid levels [27].

In summary, our results further confirm that the hub genes of MEbrown are significantly associated with macrophage lipid metabolism and cholesterol biosynthesis. These genes drive the overall gene expression of the MEbrown module to influence macrophage function in lipid metabolism.

2.1.6. Hysteresis of Marco (SR-A6) and Abca1/Abca2 Results in Elevated Intracellular Lipid Levels and Suppression of Lipid Biosynthesis-Related Genes

Due to the abnormal expression of ISGs observed in lipid-loaded macrophages in Lisa Willemsen et al.'s study [19], and our previous finding of ISG hysteresis in macrophages with a history of M1 polarization [9], We hypothesize that ISGs may be related to the expression of lipid biosynthesis-related genes, which are also the hub genes regulating the MEbrown module that we identified in M1 hysteresis. To investigate their relationship, we first calculated the Pearson correlation coefficients between the RNA expression levels of ISGs and hub genes across all macrophage polarization datasets. We found that, regardless of stimulus condition, subsets of hub genes were significantly positively or negatively correlated with the expression of ISGs (Figure 7A). Interestingly, we found that the hub genes associated with cholesterol biosynthesis, including Cyp51, Dhcr24, Fdps, Idi1, Sc5d, and Fdft1, exhibited a significantly negative correlation with the expression levels of all ISGs (Figure 7B). To investigate the negatively correlated gene expression during the repolarization process, we examined the expression profile of ISGs (Figure 7C, F), negatively correlated hub genes (Figure 7D, J), and positively correlated hub genes (Figure 7E, H) in the transitions of M0 -> M1 -> M0 and M0 -> M1 -> M2 in 96 hours. We found that, unlike the positively correlated hub genes which maintain a relatively consistent expression trend with ISGs, the negatively correlated hub genes exhibited a significant decrease in expression 6 to 12 hours after M1 stimulation. During the de-/repolarization process, these genes displayed hysteresis, remaining at relatively low expression levels and not fully recovering even after 96 hours. This trend is opposite to the ISGs, which exhibit sustained high expression that gradually decreases during de-/repolarization.

From the above results, we found that during the occurrence of M1 hysteresis in macrophages, the expression of cholesterol biosynthesis-related genes remains suppressed. In contrast, pro-inflammatory genes, represented by ISGs, continue to be expressed. This suggests that the main energy source driving the pro-inflammatory response associated with macrophage M1 hysteresis is not endogenous cholesterol biosynthesis. Studies indicate that another major pathway for increasing lipid levels in macrophages is through direct uptake of lipids from the extracellular environment [27]. To further investigate how cells replenish lipids levels during the M1 hysteresis process to secure energy sources, we examined the expression of scavenger receptors which facilitate macrophage lipid uptake to a large extent [20]. It has been reported that besides serving as pattern recognition receptors and being known to act in coordination with other co-receptors such as Toll-like receptors to generate immune responses, a critical function of scavenger receptors is to uptake LDL, VLDL, and oxidized lipoproteins [28]. After being digested in the lysosome, this process increases levels of free cholesterol and free fatty acids in the cell, subsequently promoting the activation of multiple downstream transcription factors (TFs) including LXRs, PPARs, SREBPs, and NF-κB to enhance inflammatory responses [27]. We examined the expression profiles of several major scavenger receptors including, Class A receptors SR-A1 and Marco (SR-A6), Class B receptor Cd36, Class D receptor Cd68, and Class E receptor Lox-1. We observed that SR-A1, Cd36, and Cd68 maintained high expression patterns similar to the unpolarized state following M1 stimulation, while Lox-1 quickly returned to its original low expression pattern after M1 stimulation under M0->M1->M0 (Figure 8A) and M0->M1->M2 (Figure S6). Interestingly, Marco exhibited a hysteresis effect during the depolarization process, only returning to its original expression level after 96 hours under M0->M1->M0 (Figure 8B) and M0->M1->M2 (Figure S7). Previous research has shown that in both human and conditional mouse models of prostate cancer with macrophage infiltration, the overexpression of Marco is linked to lipid intake and lipid droplet accumulation, while inhibition of Marco promotes tumor suppression in these models [29]. So the hysteresis in Marco gene expression may lead to an increase in macrophage lipid intake and accumulation of intracellular lipid levels.

Previous studies have shown that lipid homeostasis in macrophages is regulated at multiple levels, including receptor-mediated uptake, metabolism, and efflux [30]. To determine how cholesterol efflux-related genes are expressed during M1 hysteresis and how they influence intracellular lipid levels, we examined the expression of Lxra and its target genes Abca1/Abca2 which are critical factors that regulate the efflux of free cholesterol[27]. Clinically, an increase in the expression of the cholesterol transporter ABCA1 in macrophages, and in the liver, which enhances cholesterol efflux, is always indicative of an anti-atherosclerotic effect [30,31]. Lack of functional Abca2 generates abnormalities in intracellular lipid distribution/trafficking in macrophages [33]. We found that Lxra remained continuously activated after M1 stimulation, even maintaining a period of hysteresis during depolarization under M0->M1->M0 (Figure 8C) and M0->M1->M2 (Figure S8). And Lxra's target genes, Abca1 and Abca2, exhibited sustained repression post-stimulation and did not fully return to their initial expression states even after 96 hours under M0->M1->M0 (Figure 8D) and M0->M1->M2 (Figure S9). The decrease in the expression of Abca1 and Abca2 would reduce lipid efflux in macrophages, maintaining high intracellular lipid levels. To further validate the decrease of lipid efflux in macrophages, we examined two other important Lxr targets, Cd5l and Apoc1, which are highly associated with macrophage lipid transportation [33,34]. We observed that these two genes consistently exhibited reduced expression following M1 stimulation and maintained a state of hysteresis (Figure 8D).

Thus, we determined that in the M1 hysteresis, the expression of cholesterol biosynthesis-related genes is suppressed, while the upregulation of Marco leads to increased lipid uptake. Additionally, the low expression of genes such as Abca1/2, Cd5l, and Apoc1 results in inhibited lipid efflux, providing sufficient energy supply for the activation of IFN signaling driven by the high expression of ISGs and other pro-inflammatory genes we found in M1 hysteresis.

3. Discussion

Extensive prior research, as well as our studies, have confirmed that macrophages, like adaptive immune cells, can develop an immunological memory, which is reflected not only in gene expression but also in epigenetic modifications. In previous studies, we utilized in vitro data from BMDM polarization and depolarization to ascertain significant M1 hysteresis in vitro, a phenomenon closely associated with the regulation of different chromatin regions by the AP-1 family and CTCF. We are aware that sustained M1 inflammatory responses require external lipid intake or endogenous biosynthesis to continue as a source of energy. Concurrently, we observed a hysteresis phenomenon in the gene expression of numerous ISGs, which may be linked to lipid metabolism. We propose that macrophage hysteresis may be associated to some extent with macrophage lipid metabolism.

In this study, we reanalyzed RNA-seq data from macrophages during polarization and depolarization over a period of 0 to 96 hours, incorporating RNA-seq from lipid-loaded macrophages for comparative validation. Our findings include 1. Due to the substantial plasticity of macrophages and their extreme sensitivity to stimulation duration. 2. Through WGCNA analysis, we identified two gene clusters significantly related to M1 hysteresis, MEblue and MEbrown. Compared to MEblue, most MEbrown genes exhibited a late onset of gene upregulation or downregulation after 12 hours of stimulation. Functionally, MEblue genes are enriched in M1 immune functions and interferon-mediated signaling pathways, whereas MEbrown genes are highly related to lipid metabolism, cholesterol biosynthesis, and secondary alcohol biosynthesis. 3. Further exploration within MEbrown identified 22 core driver hub genes within the gene co-expression network generated through WGNCA, including those regulating cholesterol biosynthesis such as Dhcr24, Idi1, Fdps, Fdft1, and Cyp51, which exhibit either significant positive or significant negative correlations and robust PPI relationships. These genes were also found to be significantly downregulated in lipid-loaded samples. 4. Further analysis of the correlation between hub genes which regulate cholesterol biosynthesis and ISGs revealed significant negative relationships. 5. We observed that during the M1 hysteresis process, the scavenger receptor Marco, which regulates macrophage lipid intake, exhibited a hysteretic RNA expression pattern. Regulators of lipid efflux, such as Abca1 and Abca2, along with their upstream regulators, showed persistently low expression levels. These expression patterns compensate for the reduction in intracellular cholesterol biosynthesis, maintaining lipid homeostasis within macrophages and providing the necessary energy supply for the sustained activation of pathways like IFN signaling.

Integrating these findings, we can conclude our hypothesis that following depolarization or repolarization, the sustained high expression, i.e. hysteresis, of the scavenger receptor Marco on the macrophage surface membrane leads to higher intracellular lipid uptake. Concurrently, the continuous expression of the upstream mediator Lxra leads to the suppression of genes such as Abca1, Abca2, Cd5l, and Apoc1, which are in charge of lower lipid efflux in macrophages. Studies have shown that lipids uptaken by macrophages are broken down into free cholesterol (FC) and fatty acids (FA), which are then transported to the endoplasmic reticulum (ER) membrane for further processing. When deregulated lipid uptake occurs, the excessive accumulation of FC and other lipids in the ER activates sensors of the unfolded protein response, triggering ER stress, which in turn regulates the gene expression on the lipid uptake, efflux, and biosynthesis [30]. So it is highly likely that the increase in lipid uptake, combined with the reduction in lipid efflux, leads to lipid accumulation within the cells. Elevated lipid levels subsequently influence gene expression, resulting in the suppression of genes involved in lipid biosynthesis, such as Dhcr24, Idi1, Fdps, and Fdft1, to balance overall intracellular lipid levels. This phenomenon was also confirmed in our validation analysis using data from lipid-loaded samples. The accumulated intracellular lipids then serve as an energy source for the inflammatory responses generated by M1 hysteresis. Ultimately, this results in the activation of regulators such as IRFs, STATs, AP-1, and CTCF in M1 hysteresis, contributing to cellular functional heterogeneity, which has been extensively explored in our previous research [9] (Figure 8E).

Over recent years, the role of lipid metabolism has been extensively studied within the classic M1/M2 macrophage polarization paradigm. Numerous studies underscore the critical role of lipid metabolism in inflammatory macrophage polarization, supporting proper inflammatory responses through membrane remodeling and as precursors for various inflammatory mediators [35,36,37,38]. Lipids, sourced from diverse origins such as native LDL, oxidized LDL particles, or free fatty acids, and existing in different forms, are internalized by macrophages through various scavenger receptors, including CD36, MARCO, and SR-A1. Once inside the cell, most lipids—including those bound to lipoproteins and fatty acids—are transported to the lysosomal-endosomal compartment. Here, lysosomal acid lipase (LAL) catalyzes the conversion of LDL-derived cholesteryl esters into fatty acids and free cholesterol, which are then transferred to the endoplasmic reticulum where excess lipids are stored in lipid droplets [39,40,41]. These lipids continuously provide energy to the macrophage. Additionally, macrophages can actively export lipids through ATP-binding cassette (ABC) transporters. Both lipid scavenging receptors and ABC transporters play pivotal roles in dictating the metabolic fates of scavenged lipids [42,43,44].

Multiple lines of evidence substantiate the concept that lipid remodeling following macrophage activation is crucial for maintaining appropriate inflammatory responses and host defense [45,46,47]. Certain lipids also function as precursors for the production of inflammatory lipid mediators [48]. MARCO, a class A scavenger receptor, is widely expressed across various myeloid populations, including peritoneal, alveolar, splenic macrophages, Kupffer cells, and dendritic cells [49,50,51]. It is highly induced by stimuli such as LPS or granulocyte/macrophage colony-stimulating factor (GM-CSF) [52], which are associated with lipid accumulation during M1 hysteresis. Rapid induction of MARCO expression in response to infectious stimuli plays a pivotal role in mediating appropriate Toll-like receptor (TLR)-dependent inflammatory responses[51,53]. Further, evidence indicates that lipid-loaded macrophages exhibit heightened MARCO surface expression. Notably, murine and human tumor-associated macrophages (TAMs) show elevated MARCO expression, potentially linked causally to their increased lipid accumulation[54,55,56]. Experiments demonstrate that the accumulation of modified and oxidized lipids is a primary driver of inflammatory responses and is implicated in numerous diseases including atherosclerosis, steatohepatitis, obesity-induced insulin resistance, and neurodegeneration[57,58,59]. While these studies collectively validate our findings, several aspects still require attention, such as potential batch effects due to variations in experimental conditions across different datasets, and the limitation of RNA-seq data capturing only the endpoints of de-/repolarization at 96 hours. Therefore, further data and experimental validation are needed in the future to deepen our understanding of the hysteresis phenomenon in macrophages.

4. Materials and Methods

4.1. Data Preparation

We collected a total of 162 sets of bulk RNA-seq data from polarized/depolarized/repolarized BMDM (Bone Marrow Derived Macrophage) (GSE158094) and lipid-loaded macrophage (GSE118656, GSE193118) from the NCBI Gene Expression Omnibus (GEO). For the polarized/depolarized/repolarized data, BMDM were treated and polarized/depolarized/repolarized with either M1 stimulus (LPS + IFNγ) or M2 stimulus (IL-4). For the lipid-loaded data, BMDM were treated by DMSO (Dimethyl sulfoxide, a chemical that dissolves organic and inorganic substances), 1.0 μM GW3965 (a selective liver X receptor agonist), 100 ng/ml KLA (a natural endotoxin and a combination of oligosaccharide and lipid), and 50 μg/ml AcLDL (Acetylated-low density lipoprotein).

4.2. Sequencing Data Preprocessing

We assessed the quality of all sequence data using FastQC [21] and MultiQC [22]. After confirming the quality of the sequence data, we employed Trim Galore [23] to remove low-quality reads and process adaptors. Following a thorough quality check, we used STAR [24] to map the processed FASTQ files to the mm10 reference genome. Subsequently, we utilized the cuffnorm functionality within Cufflinks [25] to obtain normalized FPKM values for genes.

4.3. Batch Effect Correction, Quality Control, and Data Processing

We initially transformed the gene's FPKM (Fragments Per Kilobase of transcript per Million mapped reads) expression matrix into a TPM (transcripts per million) matrix. Using the 3167 mouse housekeeping genes (https://housekeeping.unicamp.br/Housekeeping_Genes_Mouse.RData) as a background, we applied RUVSeq (Risso et al., 2014) to remove batch effects from bulk RNA-seq data originating from different projects. Through further analysis, including the measurement of relative gene expression levels across all samples and examination of the distribution of PCA, we confirmed the successful elimination of batch effects. Subsequently, we calculated and kept the top 5000 genes with the highest Median Absolute Deviation (MAD) for downstream analysis.

4.4. Weighted Gene Co-Expression Network Analysis (WGCNA)

Based on the processed gene expression matrix, we constructed co-expression networks using the WGCNA [26] and performed gene co-expression analysis to determine the correlation between gene modules and macrophage phenotypes. After conducting hierarchical clustering of samples to remove outliers, in order to ascertain the link strengths between two genes, a co-expression network was first built using a suitable soft thresholding power (6, scale-free R² = 0.9), chosen in accordance with scale-free fitting. Then, the link between gene pair expression, or Topological Overlap Matrix (TOM) matrix, is defined to reflect the similarity between genes at both the expression and network topology levels. The dissimilarity matrix (1-TOM) derived from TOM was constructed as an additional adjacency matrix that takes topological similarity into account. By using a sensitive module detection parameter of 3 and a minimal module size of 30, hierarchical clustering of genes was used to identify important modules. High similarity gene modules were found and combined (height threshold of 0.25). Potential hub genes correlated with M1 hysteresis clusters (GS > 0.4 and MM > 0.7) were identified for downstream analysis using the gene significance (GS), which is defined as the association with each gene expression level and the clusters, and the module membership (MM), which is defined as the association with the gene expression profile and the module eigengene.

4.5. Gene Ontology Analysis and Pathway Analysis

To elucidate the functions of different gene modules, we conducted enrichment analysis for Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways using clusterProfiler [27].

4.6. Other Bioinformatics Analysis and Statistical Analysis

We used the protein-protein interaction (PPI) and gene co-expression data from the STRING database[28]. We calculated the Pearson correlation coefficient with hub gene expression and generated the PPI network graph using the igraph package [29]. Benjamini-Hochberg correction was used to compute P-values. Statistical significance was defined as P < 0.05. The remaining p-value results were as follows: “-”: P ≥ 0.05, “*”: P< 0.05, “**” :P< 0.01, and “***”: P< 0.001.

5. Conclusions

Taken together, our research elucidates a potential regulatory relationship between M1 macrophage hysteresis and macrophage lipid metabolism, which contributes to advancing our understanding of the nature of innate immune memory. Hysteresis has also been demonstrated in other immune-related cells and tissue-resident cells. Further efforts are required to detail the mechanisms of hysteresis and assist in the design of novel vaccines.