Systems genetics identifies a macrophage cholesterol network associated with physiological wound healing

Wound healing is the repair of damaged and injured skin. It is a complex and highly dynamic process consisting of overlapping phases that include inflammation, tissue formation and remodeling, which can lead to the production of a nonfunctioning mass of fibrotic tissue known as the scar (1). Impaired wound healing is a critical socioeconomic problem in the Western world (2), as chronic nonhealing wounds represent a clinical burden associated with significant healthcare costs (3).

Among other cells, macrophages are pivotal effector cells during wound healing. They dynamically change their activity throughout the healing process as their immune/phagocytic properties gradually convert into a more reparative and immunomodulatory phenotype (4, 5). Macrophages are present during all stages of the repair process and their depletion with anti-macrophage serum results in impaired wound healing (6). In addition, conditional depletion of macrophages during the different phases of wound healing has further confirmed their context-dependent plasticity (7). Different macrophage polarization stages correspond to different phases of tissue repair: proinflammatory macrophages infiltrating during the initial inflammatory phase are later antagonized by immunomodulatory macrophages that promote tissue repair and fibrosis (5, 8).

Animal model systems can provide a useful tool to map pathways and genes underlying (patho)physiological traits (9, 10). For instance, several inbred mouse strains differ significantly in their rate of wound healing (11) and some of these strains were used in genetic linkage studies aiming to map quantitative trait loci (12, 13) controlling physiological wound healing. A significant advance in gene mapping strategies came from genome-wide association studies (GWAS) carried out in genetically outbred mice, which dramatically improved mapping resolution over traditional quantitative trait locus (QTL) studies (14). Outbred mice, which are bred to maintain maximum heterozygosity (15), offer the advantage of gene-level mapping resolution due to relatively restrictive haplotype blocks and a large catalog of variants that can be tested for association (16).

Previously, a GWAS of multiple complex traits using a commercially available outbred mouse population achieved single-gene resolution mapping for multiple traits, including behavioral, physiological, and tissue phenotypes (14). The outbred mouse population offers a lower significance threshold and a smaller sample size in comparison with human GWAS (14). Quantitative variation in wound healing was one of the 92 phenotypes measured in almost 2,000 outbred mice, and importantly, this trait was associated with a single gene (Bmp2) following genetic mapping using low-coverage sequencing (14). Bmp2, and more generally Bmp signaling, was previously found to be involved in wound healing (17, 18) but the significance of this result in terms of the effector cell function in wound healing remained to be determined. Genotype-phenotype relationships, often first identified by GWAS, could be integrated with systems genetics approaches aiming to reveal the cellular and molecular pathways underlying complex traits (19). Gene- and/or network-level strategies in the relevant effector cells could complement GWAS in functionally annotating genes (20, 21). We therefore hypothesized that naturally occurring germline sequence variation could affect macrophage function associated with physiological wound healing.

To test this hypothesis, we used the ear punch model to study the rate of wound healing in 1,378 outbred mice. We integrated phenotypic and genotypic analyses in this population (i.e., ear area closure after healing and genome-wide single-nucleotide polymorphism [SNP] genotyping) with whole-genome mRNA sequencing (RNA-seq) in primary macrophages from mice with extreme wound healing phenotypes (i.e., fast and slow healers from normal distribution of the rate of wound healing). This led to the identification of a gene network enriched for fatty acid and cholesterol biosynthesis that (a) associates with the rate of wound healing and (b) is under strong genetic control in trans by the transcription factor Runx2 located on mouse chromosome 17. We showed that inhibiting the DNA binding activity of Runx2 affects the expression of the genes in this network. We then targeted Fasn, which catalyzes the synthesis of long-chain fatty acids and also shows a regulatory role in cholesterol synthesis (22). In vivo pharmacological blockage of Fasn results in delayed wound healing and increased macrophage recruitment in the skin wound, suggesting the persistence of an unresolved inflammation. These results suggest the importance of fatty acid/cholesterol biosynthesis in the healing process.

The rate of wound healing is normally distributed in a large outbred mouse population. Outbred mice are large, genetically heterogeneous colonies that commercial mouse breeders maintain. Typically, these colonies present the necessary genetic structure and variation to allow high resolution mapping of complex traits, which can be exploited for gene identification by GWAS. A 2-mm ear punch was performed in 16-week-old CFW outbred mice and wound healing monitored after 5 weeks (ear area after healing) in 1,378 outbred animals as part of a larger study combining phenotyping of multiple traits and genotyping by low-coverage sequencing (14). We reasoned that the major molecular determinants underlying this process could be more easily captured by selective RNA-seq of the macrophages of fast and slow healers (23, 24). One hundred forty-six mice were selected for this purpose at the extreme of the ear area distribution (62 fast healers and 84 slow healers) and their bone marrow–derived macrophages (BMDMs) were analyzed by RNA-seq (Supplemental Table 1; supplemental material available online with this article; https://doi.org/10.1172/jci.insight.125736DS1). The available resources for these 146 mice were (a) quantitative differences in wound healing indicated as ear area, (b) BMDM mRNA expression by RNA-seq, and (c) genome-wide SNP genotyping imputed from low-coverage sequencing (Figure 1).

Figure 1

A systems genetics approach for deciphering macrophage-dependent wound healing. Study design showing 1,378 outbred mice subjected to a 2-mm ear punch and wound healing measured after 5 weeks. One hundred forty-six mice were selected at the extreme of the ear area distribution (62 fast healers and 84 slow healers) and their macrophages analyzed by RNA sequencing (RNA-seq). The available resources for these 146 mice are (i) quantitative differences in wound healing indicated as ear area following punch, (ii) RNA-seq in bone marrow–derived macrophages (BMDMs), and (iii) genome-wide SNP genotyping by low-coverage sequencing.

A macrophage-mediated healing network maps strongly to the Runx2/Supt3h locus. Instead of mapping the genetic control of single gene expression, we used a gene coexpression module/network identification and mapping strategy, which we have previously applied to identify trans regulators of networks in complex diseases (25–27). First, we inferred macrophage gene coexpression modules from the transcriptome of fast and slow healers. This analysis identified 40 coexpression modules, each containing 30 to 1,151 genes, listed in Supplemental Table 2.

To identify key regulators of the transcriptional programs captured by the coexpression modules inferred in the fast and slow healers, we carried out network QTL mapping (28). This analysis resulted in a unique module that was strongly associated (Bayes factor [BF] = 211,096.6) with a SNP located on mouse chromosome 17 (Figure 2A). The strong association of this module (designated as “macrophage-mediated healing network” [MMHN]) is also illustrated by its 195-fold increase in BF when compared with the second strongest associated module (BF = 1,084.5). Rather than adopting a nearest gene approach to investigate the putative target gene for the regulatory SNP, we interrogated publicly available high-resolution chromatin interaction data (29). This revealed that the regulatory SNP (chr17_45131552) is located within a highly conserved and well-characterized topologically associated domain (TAD), containing both Runx2 and Supt3h genes (Figure 2B). TADs have been proposed to represent the regulatory domains of the genes located within them (30), primarily by delimiting the scope of enhancer-promoter interactions (31, 32). The TAD that contains the SNP associated with the MMHN is highly conserved across tissues and species (33–35) and found in human macrophage Hi-C data (36) (Supplemental Figure 1). The analysis of additional candidate trans-acting regulators upstream the macrophage gene coexpression networks associated with wound healing is reported in Supplemental Figure 2.

Figure 2

Strong association of a macrophage-mediated healing network with the Runx2/Supt3h locus. (A) Genome-wide mapping results from multivariate Bayesian mapping of 40 macrophage coexpression modules in the outbred mouse population showing fast and slow healing (n = 146). The coexpression module with the highest association (Bayes factor [BF] of 211,096) is M30 (also designated as macrophage-mediated healing network, MMHN) through SNP (chr17_45131552). The results for the top 5 modules with strongest genetic association are presented in a table (left). (B) M30 is associated with chr17_45131552 on mouse chromosome 17, which is located within a well-conserved and tissue invariant (Supplemental Figure 1) topologically associated domain (TAD) containing Runx2 and Supt3h in mouse embryonic stem cells. As TADs can represent the regulatory domain of the genes they contain, the location of chr17_45131552 within this TAD implicates it in the transcriptional regulation of either Runx2 or Supt3h. See supplemental methods for more details.

The MMHN correlates with wound healing and is under Runx2 trans regulation. Functional enrichment analysis of the MMHN showed strong and specific enrichment for cholesterol and fatty acid biosynthetic processes (Figure 3A). Closer inspection of the network genes showed the presence of Hmgcr, the rate-limiting enzyme in cholesterol synthesis, together with Lss and Dhcr7, which are positioned downstream of Hmgcr. Furthermore, Fasn, which catalyzes the synthesis of palmitate, has been recently linked to cholesterol biosynthesis in macrophages (22).

Figure 3

The macrophage-mediated healing network is enriched for cholesterol/lipid biosynthesis and is partly controlled by Runx2 activity. (A) Overview of coexpression module M30 (MMHN). Top left: Gene Ontology (GO) functional enrichment of genes in MMHN (177 genes). Network graph with the genes (nodes) in the MMHN that have any known protein-protein, coexpression, and database interactions (edges) in the STRING protein database (67). Node size represents the degree of probability of association for each gene to the regulatory SNP chr17_45131552 (computed by Bayesian multivariate mapping). Genes with a predicted RUNX2 transcription factor binding site (TFBS) in their promoter are shown with a yellow border. Genes annotated with the GO functional term “lipid biosynthetic process” are presented in green. Runx2 (trans regulator) is colored in orange. See supplemental methods for additional details. (B) Genes whose expression levels correlate with Runx2 also correlate with the rate of healing. Spearman’s correlation of the expression level of each gene in the MMHN with Runx2 expression levels (normalized variance-stabilized counts), against the correlation between the expression levels of each gene in the MMHN and the rate of healing (ρ = 0.74, P = 1.04^–31, n = 146 mice). (C and D) qRT-PCR results of a subset of MMHN genes showing positive and negative correlation with Runx2, respectively, following Runx2 blockage with CADD522 (20 μM) in mouse BMDMs. Data are expressed as mean ± SEM (n = 4). *P < 0.05, **P < 0.01, ***P < 0.005, ****P < 0.001 by 2-tailed Student’s t test.

We investigated Runx2 transcription factor binding sites (TFBS) in the promoter of MMHN genes and identified 70 genes with a predicted Runx2 TFBS (Figure 3A), further suggesting a trans-acting genetic regulation of this network by Runx2.

We set out to inspect the relationship between the transcriptional regulation of the MMHN by Runx2 and wound healing using a 2-step approach. We first correlated the variance-stabilized transformed (VST) gene counts of all transcripts in the network with Runx2 expression. We then correlated the VST gene counts of all the transcripts in the network with the rate of wound healing (see supplemental methods). These analyses showed that the transcripts in the network that positively correlate with Runx2 expression also positively correlate with healing, and the same is true for those correlating negatively (ρ = 0.74, P = 1.04^–31) (Figure 3B). Taken together, these results showed that the genetically controlled MMHN is enriched for genes implicated in cholesterol and fatty acid biosynthetic processes, whose expression correlates with Runx2 mRNA expression and the rate of wound healing.

In order to validate Runx2 regulation of MMHN genes in vitro, we cultured BMDMs and incubated them for 48 hours with 20 μM CADD522, a drug previously shown to selectively block Runx2 DNA binding activity (37). Expression of genes positively and negatively correlated with Runx2 was then assessed by quantitative reverse-transcription PCR (qRT-PCR). In line with the correlation identified in silico, we found significant downregulation of transcripts showing a positive correlation with Runx2 (Med23, Med17, Rsrp1, P4ha1, Rbm5, Nfkbiz, Kdm1a, and Jmjd4) (Figure 3C), and significant upregulation of the genes having negative correlation with Runx2 (Fasn, Acly, Hmgcr, Prkch, Fdps, Dhcr7, Mvk, Lss, and Poc1b) (Figure 3D). These results confirm the role of Runx2 as a trans-genetic regulator of MMHN genes in murine BMDMs.

Pharmacological targeting of Fasn in vivo delays wound healing in rats. In order to assess the effect of the modulation of the MMHN in vivo, we pharmacologically targeted one of its genes, Fasn, using cerulenin, a fatty acid synthase blocker (38). We reasoned that the modulation of a regulator of both fatty acid and cholesterol pathways (22) might accurately reflect MMHN in the healing process.

Full-thickness, 10-mm splinted excisional wounds were inflicted on the dorsal skin of 13-week-old Lewis (LEW) rats. By splinting the wounds, contraction is minimized to favor granulation and reepithelialization, which characterizes wound healing in humans. Skin wounds were treated with 300 μg of cerulenin in propylene glycol at 0, 3, and 6 days after wounding. Histological analysis was performed at day 8 after wounding (Figure 4A). Wounds treated with cerulenin showed a significant delay in healing at both day 6 and day 8 after wounding (Figure 4, B–D) and were associated with increased fatty acid synthesis/cholesterol genes (Supplemental Figure 3). Immunohistochemical analysis of the macrophage marker CD68 showed a significantly increased number of CD68⁺ cells in the wounds treated with cerulenin (Figure 4, E–G), suggesting the delayed wound healing observed with cerulenin treatment can be accompanied by an unresolved inflammation characterized by persistent macrophage infiltration. Overall, these results suggest that the modulation of the MMHN in vivo can affect the rate of wound healing in rats.

Figure 4

FASN inhibition with cerulenin affects wound healing in the skin. (A) Experimental design. Full-thickness, 10-mm splinted excisional wounds were performed on the dorsal skin of 13-week-old Lewis (LEW) rats. Wounds were treated with 300 μg of cerulenin in propylene glycol at 0, 3, and 6 days after the excision or with propylene glycol alone (vehicle). Tissues were collected 8 days following wounding. (B) Representative images of wounds immediately following surgery (day 0) and after 3, 6, and 8 days for control and cerulenin-treated rats. (C) Quantification of wound area at days 0, 3, 6, and 8 after wounding. Values are expressed as the healing ratio in comparison to day 0 and are represented as mean ± SEM (n = 6 controls and n = 8 cerulenin treated). (D) Representative image of H&E-stained wounds 8 days after wounding in control and cerulenin-treated rats. Black arrows indicate the edges of granulation tissue. (E) Representative CD68 (rat ED-1) immunohistochemistry within the wound edges in control and cerulenin-treated rats. (F) CD68 staining showing the center of the granulation tissue. (G) Quantification of CD68⁺ cells per high-power field (HPF, ×20) in control and cerulenin-treated wounds. Data are expressed as mean ± SEM (n = 5 controls and n = 8 cerulenin treated. Scale bars: 1,000 μm (D and E) and 100 μm (F). *P < 0.05 by 2-tailed Student’s t test.

Physiological wound healing involves highly interactive and dynamic processes including inflammation, angiogenesis, matrix deposition, and cell recruitment, which can only be under polygenic and complex genetic control in genetically outbred populations. Evidence from several studies conducted in murine and porcine models support the idea that physiological wound healing is indeed heritable (39–41). While there is genetic association evidence for the pathological forms of wound healing such as diabetic foot ulcers (42), keloids (43), and chronic venous disease (44), the only GWAS for physiological wound healing was conducted in a commercially available outbred population, which allowed single-gene-resolution mapping and identification of causal variants (14). Here, we performed RNA-seq of macrophages derived from fast and slow healers identified from the same outbred population. We show that a MMHN controlled by Runx2 associates with wound healing through cholesterol/fatty acid biosynthesis. These results complement the GWAS carried out in the same population where Bmp2 was the only gene identified as being associated with wound healing. Runx2 is downstream of Bmp2 in osteoblasts (45, 46), which suggests that our systems genetics approach broadened the GWAS results by identifying variants that could affect wound healing through macrophage function.

In our gene network mapping strategy, we used BMDMs as a proxy for wound macrophages in the outbred mice. Besides the methodological advantages (i.e., cell number, control of environmental factors), there is evidence based on the green fluorescent protein–bone marrow chimeric mouse models of wound healing showing that macrophages at the wound healing originate predominantly from the bone marrow (47). In addition, it is possible that the RUNX2 trans-acting genetic regulation of the network is conserved across the myeloid lineage, given previous studies showing the conservation of cis-expression QTLs in monocytes/macrophages in general (48).

The evidence for Runx2 involvement in skin wound healing is scarce (49), despite a large body of work implicating its role in skeletal development (50) and more specifically in osteoblast and terminal chondrocyte differentiation (51). In keeping with this, SNPs located at the Runx2-Supt3h locus have been associated with bone and cartilage phenotypes by GWAS (52), which suggests shared genetic control of physiological wound healing rates and bone and cartilage phenotypes. The role of Runx2 in immune cells and more particularly in macrophages is not well documented. It has been reported that RUNXs are involved in NOD2 stimulation in human macrophages (53) and that RUNXs mediates egress of plasmacytoid dendritic cells into the circulation (54), controlling their development and migration (55). The exact role of Runx2 in macrophage function remains to be established in the context of macrophage plasticity during the healing process. In addition, the previously established physical proximity between the promoters of Runx2 and Supt3h (34) suggests a coordinated trans regulation of the MMHN. Thus, we cannot rule out a subgroup of genes belonging to the network that are under Supt3h control, especially given its ability to regulate ADAMTS13 activity (56), which has been linked to wound healing (57).

Recent advances in the metabolic control of innate immune function will have a considerable impact on manipulating the rate of wound healing, given the known plasticity of macrophage function depending on context-dependent metabolic pathways (5). The importance of targeting HIF-1α, a hallmark of proinflammatory, glycolytic macrophage activation in wound healing (58) and, more recently, the lactic acid–driven delivery of repair agents such as CXCL12 (59) show how the healing process is tightly dependent on the metabolic status of the wound microenvironment.

The regulation of the macrophage inflammatory response by cholesterol metabolism previously revealed a role of lanosterol (the first sterol intermediate in the cholesterol biosynthetic pathway) in modulating the TLR4 response (60). Similarly, desmosterol, the last intermediate in the Bloch pathway of cholesterol biosynthesis, has been linked to the suppression of inflammatory-response genes in macrophages (61). Furthermore, the oxysterol 25-hydroxycholesterol (25-HC) acts by antagonizing sterol response element–binding protein (SREBP) processing to reduce IL-1β transcription and to broadly repress IL-1–activating inflammasomes (62). Selective deletion of FASN in macrophages protects again inflammation (63), and interestingly, macrophages lacking FASN had reduced cholesterol levels (63), suggesting the regulatory role of FASN in both cholesterol and fatty acid synthesis, a finding recently confirmed in proinflammatory macrophage activation (22). Since FASN is required for proinflammatory macrophage activation (22), our in vivo results obtained with cerulenin suggest that inhibiting proinflammatory macrophage activation delays wound healing. This is in line with the role of proinflammatory macrophages in promoting angiogenesis, an essential process during early repair (64). However, further studies will be required to establish the role of macrophage-specific fatty acid synthesis and the respective contribution of each pathway (cholesterol vs. fatty acid synthesis) in wound healing. This is particularly important in view of a study showing a beneficial effect of statins (simvastatin) in wound healing (65). The negative correlation between Runx2 and Fasn at the mRNA level within the module could be the result of a feedback mechanism (66), which may reflect the differences between mRNA and protein levels in terms of their correlation with each other.

In summary, we report a MMHN under Runx2 trans regulation, revealing the fatty acid/cholesterol pathways that could explain differences in physiological wound healing. Our systems genetics approach complements previous GWAS findings and highlights the BMP/Runx2 axis as a major genetic determinant of physiological wound healing. Understanding the genetic basis of wound healing will be pivotal in targeting the pathologies characterized by defective repair mechanisms.

The mouse macrophage RNA-seq data have been deposited in the NCBI’s Gene Expression Omnibus database (GEO GSE112171). Phenotype data are available in Supplemental Table 1. See full methodological description of data, analyses, and experimental procedures, including Supplemental Tables 3 and 4, in the supplemental information.