Combined Metallomics/Transcriptomics Profiling Reveals a Major Role for Metals in Wound Repair

Introduction

Skin is the first barrier to noxious environmental onslaught and as such, must be repaired quickly and efficiently when damaged. This wound healing process involves a cascade of cellular signalling events that actuate functional responses to confer tissue repair (Wilkinson and Hardman, 2020). These events are grouped into characteristic wound healing stages (haemostasis, inflammation, proliferation and extracellular matrix remodelling) that are tightly controlled but remain highly dynamic and often overlap (Greaves et al., 2013). Despite a high level of regulation and redundancy, aberrations to normal skin repair are surprisingly frequent. Over-exuberant healing leads to excessive scarring and keloids, while insufficient healing results in infection and non-healing “chronic” wounds. Chronic wounds remain a major socio-economic burden, with existing treatments often inadequate (Frykberg and Banks, 2015). Thus, there is an ongoing need to fully understand the molecular and cellular aspects of normal healing, in order to effectively address healing pathology.

Metallomics, the global characterisation of metals and metalloids in biological systems (Haraguchi, 2017), is a rapidly emerging field. When combined with existing proteomics, genomics or metabolomics approaches, it can provide an integrated, systematic understanding of biology (Singh and Verma, 2018). Metals are fundamental for almost all biological processes including metabolism, energy transduction, gene expression and cell signalling (Yannone et al., 2012). This is best demonstrated by the fact that almost half of all enzymes contain a transition metal core (e.g., copper, zinc or iron) that is necessary for appropriate enzymatic function (Andreini et al., 2008). Paradoxically, despite their clear biological importance, current understanding of the temporal and spatial distribution of metals in tissues, and subsequent integration with other omics data, remains fairly limited.

There is existing literature pertaining to the importance of specific endogenous metals in skin biology and wound healing. Calcium flashes occur within seconds to trigger a healing response (Razzell et al., 2013), while calcium is also widely reported to regulate epidermal stratification and barrier integrity (Elias et al., 2002), platelet activation (Elaïb et al., 2016), cellular migration (Mandeville and Maxfield, 1997; Wei et al., 2009), macrophage phagocytosis (Vashi et al., 2017) and cellular proliferation (Berridge, 1995; Schwarz et al., 2006). Each of these processes is crucial for effective skin wound healing. Wound-relevant roles for other metals are also emerging, such as iron, which is involved in haemostasis (Pretorius et al., 2013), inflammation (Haldar et al., 2014), and extracellular matrix (ECM) deposition and remodelling (Wilkinson et al., 2019b; 2019c). Given that these reductionist studies have revealed important mechanistic roles for individual metals in specific repair processes, it is somewhat surprising that little attempt has been made to systematically characterise the temporal distribution of metals across the wound repair continuum, and globally integrate these data with molecular and cellular level changes.

In the present study, we show for the first time that temporal alterations in the wound metallome correlate with both metal- and injury-linked transcriptional changes across normal murine skin repair. Specifically, wound-induced genes are highly enriched for metal functional annotation, clustering into particular biological processes that correlate with hitherto unreported temporal changes in wound metal distribution. Our findings are corroborated further by exploration of metal-linked transcriptional changes in human wounds.

Materials and Methods

In Vivo Wound Healing

Female C57BL/6J mice (8–10 weeks old) were purchased (Envigo Ltd., United Kingdom) and housed at the Biological Services Facility (The University of Manchester, United Kingdom). Animal experimentation was performed under full United Kingdom Home Office approval (project licence: 70/8136). Food and water was given ad libitum and cages were kept under a constant 12-h light-dark cycle. Mice were anaesthetised and 2 mm × 6 mm full thickness excisional dorsal wounds were created using trace metal free instruments. Wounds were positioned 2 cm caudal to the base of the skull, 0.5 cm either side of the midline. Mice were administered buprenorphine post operatively and wounds were left to heal via secondary intention. Skin (D0) and wounds were collected at days one (D1), three (D3), seven (D7) and fourteen (D14) post-wounding (n = 5 mice per time point). Wound samples were carefully dissected to include the full wound bed and 3–4 mm of peri-wound tissue. Skin and wounds were flash frozen in liquid nitrogen and stored at −80°C prior to ICP-MS and RNA isolation or fixed in neutral buffered formalin for histological assessment.

Metal Analysis

Frozen tissue was freeze dried and digested in a 50:50 mixture of 30% hydrogen peroxide (Sigma-Aldrich, Dorset, United Kingdom) and trace metal free nitric acid (Thermo Fisher Scientific, Loughborough, United Kingdom) in a MARS 6 microwave with a MARSXPRESS™ vessel (CEM Microwave Technology Ltd, Buckingham, United Kingdom) as previously described (Wilkinson et al., 2019c). DOLT-5 reference material (National Research Council, Canada) was prepared and analysed alongside digested tissue samples. Internal standards (Rh at ∼100 ppb) were used as calibration and a water reference standard (CRM 1643E) was included for internal standard verification. Samples were run on an Agilent 7800 inductively coupled plasma mass spectrometer (Agilent Technologies, Cheshire, United Kingdom). The operating conditions were: Carrier gas flow of 1.05 L/min and spray chamber temperature of 2°C; plasma mode was general purpose, RF forward power was 1,550 W with a sampling depth of 10 mm; He cell gas flow of 5 ml/min and Ni cones; the extraction lens was 1 V, with kinetic energy discrimination of 5 V. Major isotopes of the elements were monitored, with a mass to charge = 1. Working conditions were checked at the start of the day by injecting the standard solution of several isotopes at 1 ppb level, tuned and the ionisation signal compared to the previous day. The relative standard deviation among replicates was evaluated as standard. A precision of <5% with sensitivity of ∼ppb was routinely achieved.

Histology

Sections (5 µM thick) were taken from paraffin embedded blocks. Slides were dewaxed and brought to dH₂O down an ethanol gradient. Masson’s trichrome staining, and immunoperoxidase staining for neutrophils and macrophages, was performed as described in Wilkinson et al. (2019a). For immunofluorescent staining, antigen retrieval was achieved with citrate buffer and sections blocked in goat serum and M.O.M block (Vector Laboratories, CA, United States). Rabbit anti-keratin 14 (clone: Poly19053; Biolegend, CA, United States), mouse anti-TNF-α (clone: 52B83; Abcam, Cambridge, United Kingdom) and rat anti-CD107b (clone: M3/84; BD Biosciences, NJ, United States) primary antibodies were detected with Alexa Fluor conjugated secondary antibodies (Thermo Fisher Scientific) and slides were mounted with MOWIOL 488 containing DAPI (Thermo Fisher Scientific). Slides were imaged on an LSM 710 confocal laser scanning microscope (Carl Zeiss, Oberkochen, Germany) using a 20× objective lens and 405-nm diode, 488-nm argon, and 561-nm diode-pumped solid-state lasers. Staining was quantified using ImageJ v.1.8.0 (National Institutes of Health, MD, United States).

Extraction of UniProt Metal Lists

To initially investigate the link between wound transcriptional changes and metals, we extracted lists of metal-linked genes from the UniProt Knowledge Base—a collated resource providing comprehensive protein functional annotation. UniProt entries are tagged with keywords to enable retrieval of specific protein subsets. To extract metal-associated entries, we retrieved subsets using the following keywords: calcium (KW-0106); magnesium (KW-0460); iron (KW-0408); zinc (KW-0862); copper (KW-0186) and; manganese (KW-0464). To be included in a metal keyword list, proteins must bind at least one metal atom, or be functionally dependent on the metal in question. Each metal keyword list was next filtered to include only reviewed annotations, and then filtered by species of interest (e.g., mouse for RNA-Seq data and human for microarray comparison). Finally, gene name entries for each protein were extracted for functional annotation and overrepresentation analysis (described below).

RNA Isolation

RNA was extracted from skin, wounds and keratinocytes. Tissue was homogenised (T10 ULTRA-TURRAX^®, IKA, Oxford, United Kingdom), and cells vortexed, in TRIzol^® reagent (Thermo Fisher Scientific) prior to chloroform phase separation. RNA was purified and eluted using the PureLink™ RNA Mini Kit (Thermo Fisher Scientific) following manufacturer’s instructions.

RNA Sequencing

RNA-Sequencing (RNA-Seq) was achieved using the paired-end method and Illumina platform. Error rate distribution among reads was <0.04%, while over 90% of reads from each sample passed quality control. The sequenced library was aligned to the reference genome using HISAT2 and assembled to map the transcriptome (service provided by Novogene Ltd., Cambridge, United Kingdom).

Differential expression analysis was performed on RNA-Seq count data using “DESeq2” (Love et al., 2014) in R.v.3.6.1 (R Core Team, 2020) with LFC shrinkage (Zhu et al., 2019). Genes were deemed differentially expressed when expression between experimental groups was above 1.5-fold different (log₂ = 0.58), with 5% alpha level significance following multiple correction to reduce false discovery rate (Benjamini-Hochberg procedure). PC analysis was performed following variance stabilising transformation of count data. For hierarchal clustering, the top 250 most significant genes were filtered from transformed counts and clustered using Euclidian distance and Ward D2’s method within R package “gplots” (Warnes et al., 2020).

Overlap in gene subsets was determined in InteractiVenn (Heberle et al., 2015), while volcano plots were made from gene lists using the R package “EnhancedVolcano” (Blighe, 2019). For functional annotation, Ensembl gene identifiers were scrutinised in the Database for Annotation, Visualisation and Integrated Discovery v.6.8 (Huang et al., 2009), with top relevant annotation (gene ontology, KEGG pathway, UniProt KB) plotted against Benjamini-Hochberg adjusted p value. Markov clustering was performed on count data filtered for the combined list of upregulated DEGs using the Graphia tool (Kajeka Ltd., United Kingdom). Here, Pearson’s correlation for similarity was set to 0.94 and nodes clustered using the Markov clustering algorithm based on expression profiles. Gene clusters were manually annotated using functional annotation (as in Nirmal et al., 2018). Metal overrepresentation was determined using software to calculate the significance in overlap between two sets of genes (see Plaisier et al., 2010). Ingenuity Pathway Analysis (Qiagen, Manchester, United Kingdom) was used to determine top upstream regulators in specifically identified clusters.

Quantitative Real-Time qPCR (qRT-PCR)

Bioscript™ (Bioline, London, United Kingdom) and random primers (Promega, Southampton, United Kingdom) were used for reverse transcription of purified RNA. qRT-PCR was performed on a CFX Connect™ thermocycler (Bio-Rad, Hertfordshire, United Kingdom) with Takyon™ SYBR mastermix (Eurogentec, Hampshire, United Kingdom). Primer sequences for RNA-Seq validation and keratinocyte experiments are provided in Supplementary Table S1.

Keratinocyte Experiments

Keratinocytes were isolated from mouse and human skin for TNFα inhibitor experiments. Human skin was collected from Castle Hill Hospital (Cottingham, Hull, United Kingdom) under full, informed patient consent and institutional approval (LREC: 17/SC/0220). Human and mouse skin was incubated in 0.2% Dispase II (Merck Life Science United Kingdom Ltd., Dorset, United Kingdom) at 4C overnight to separate epidermis and dermis. The epidermis was then dissociated in 0.25% Trypsin (Thermo Fisher Scientific). Trypsin was neutralized with heat-inactivated Foetal Bovine Serum (Thermo Fisher Scientific) and the cell suspensions pelleted. Human epidermal keratinocytes (HEKs) were cultured in Epilife media containing Epilife defined growth supplement on coating matrix (Thermo Fisher Scientific). Mouse epidermal keratinocytes (MEKs) were cultured in CnT-07 medium (CELLnTEC, Bern, Switzerland) on collagen IV (Corning, NY, United States). Confluent keratinocytes were scratched with a 1 ml filter tip to induce a wounding response. Treatment conditions were low calcium (Epilife or CNT07 alone), high calcium (media plus 2mM CaCl₂) and high calcium with TNFα inhibitor (10 μM C87, Biotechne, Oxford, United Kingdom). HEKs (n = 3 donors across three independent experiments) and MEKs (n = 5 pooled, three independent experiments) were collected for qRT-PCR.

Statistical Analysis

Mean ± standard deviations of the mean (SEM) were used for non-transcriptional data. One-way ANOVAs were performed on this data with Tukey’s post-hoc using R. v.3.6.1. Significance was determined at the alpha 0.05 level. Pearson’s correlation was performed on RNA-Seq qRT-PCR validation data using R. v.3.6.1.

Results

Global Profiling Reveals Temporal Changes in the Wound Metallome During Normal Repair

Inductively coupled plasma mass spectrometry (ICP-MS), a highly sensitive elemental analysis technique (Lee et al., 2014), was used to measure metal abundance across an acute wound healing time-course. Tissue concentrations of eight skin-relevant metal elements were compared between unwounded skin (D0) and wounds at days 1, 3, 7 and 14 (D1, D3, D7 and D14) post-wounding (PW; Figure 1; absolute values provided in Supplementary Table S2). Healing stages were confirmed via histological analysis (Supplementary Figure S1).

FIGURE 1. Temporal profiling reveals global changes in the metallome throughout normal wound repair. ICP-MS was used to measure metal abundance in skin and acute healing wounds. Absolute abundance of individual metals in tissue (A–H). Heat map relative to normal skin (D0) where the saturated red value is 7 (I). D1 = day 1 post-wounding. Schematic showing temporal metal changes alongside the main wound healing phases (J). Mean ± SEM. n = 5 mice per group. One-way ANOVA with Tukey’s post-hoc analyses were performed. * = p < 0.05, ** = p < 0.01 and *** = p < 0.001. Asterisk alone versus D0.

Calcium levels were significantly elevated at D3 (p < 0.05) and D7 (p < 0.01) PW compared to D0 but dropped substantially at D14 PW to levels lower than observed at D0 (p < 0.01 versus D7; Figure 1A). Although magnesium concentration was not significantly altered between wounds and skin (Figure 1B), levels increased at D3 (versus D1; p < 0.05) and remained higher than D1 wounds throughout the 14 day period. Iron was substantially elevated at D7 and D14 versus D0 (Figure 1C, as previously reported in Wilkinson et al., 2019c), while zinc levels were not altered throughout the wound healing time course (Figure 1D). Aluminium (Figure 1E) was not significantly changed between skin and wounds, however levels at D14 were increased versus D1 (p < 0.01) and D3 (p < 0.01). Copper (Figure 1F) was highest at D7 PW (compared to D0, p < 0.05, and D1, p < 0.01). By contrast, manganese (p < 0.05; Figure 1G) and cobalt (p < 0.001; Figure 1H) showed the greatest elevation in levels at D14. The relative fold change in each metal versus D0 is shown in a heat map (Figure 1I), while the temporal profile for each metal is illustrated against specific stages of wound repair in Figure 1J. Collectively, these data demonstrate clear temporal changes in endogenous metals throughout normal murine wound healing.

Independent of Context, Metal-Linked Genes Are Highly Annotated for Wound-Associated Biological Processes

As metals were temporally altered throughout normal wound repair, we next asked whether metal-linked genes were independently linked to skin- and wound-associated pathways and processes (i.e., would there be global association between the metallome and transcriptome). Specifically, “metal” gene lists were extracted from the UniProt Knowledge Base (KB) and scrutinised by functional annotation analysis to identify overrepresented skin- and wound-related biological processes (Figures 2A–F). Here, genes from the (mouse-specific) calcium (KW-0106), magnesium (KW-0460), iron (KW-0408), zinc (KW-0862), copper (KW-0186) and manganese (KW-0464) UniProt KB lists were highly annotated for many wound-linked functional annotation terms, in a metal-specific manner (Supplementary Table S3). For example, cell adhesion (UniProt KB: KW-0130) and ECM (UniProt KB: KW-0272) annotations were overrepresented only in the calcium (KW-0106) gene list, while platelet activation (KEGG pathway: mmu04611) was highly enriched in magnesium (KW-0460) and manganese (KW-0464) lists.

FIGURE 2. Metal-specific genes are associated with biological processes independent of context. Functional annotation was performed on metal gene lists from the UniProt knowledge base [KB (A–F)], demonstrating high annotation of processes linked to skin and wound healing. KW, keyword identifier; Numbers of genes, left of bars; ECM, extracellular matrix; P., pathway; Act., activity. The top 1,000 expressed skin and organ genes show differential overrepresentation of UniProt KB metal genes [observed vs expected fold change (G)]. Metal genes (combined UniProt KB list) are significantly overrepresented in the wound healing gene ontology group [GO: 0042060; (H)]. NS, non-significant; Padjust, Benjamini-Hochberg adjusted p value.

Given the strong context independent correlation between metal-linked genes and wound-associated processes, we then asked whether metal-linked genes were overrepresented in specific murine tissues; skin (our data) and other organs (Söllner et al., 2017). To do this, we evaluated the number of metal-linked genes found in the top 1,000 expressed genes (based on RNA-Seq count) for each organ. In most cases, we observed more UniProt KB metal genes than would be expected by chance, with clear organ-specific differences in relative metal overrepresentation (Figure 2G; Supplementary Table S4). For example, in the skin, genes restricted to the manganese UniProt KB list were >7 fold overrepresented, while genes in the cobalt UniProt KB list displayed <1.5 fold overrepresentation. In line with known biological association (Graham et al., 2007), the most overrepresented (>8 fold) metal UniProt KB list across all organs was iron in the liver.

At an even broader level, genes linked to any one of the seven UniProt metals (combined UniProt KB lists: KW-0106, KW-0460, KW-0408, KW-0862, KW-0186, KW-0464 and KW-0170) showed significant overlap with known wound healing genes (GO group: 0042060; Figure 2H). Here, 92 out of 390 wound healing genes (23.6%) were also metal-linked, a 2.85-fold overrepresentation (−log₁₀ p = 19.56; Supplementary Table S5). Taken together, these data clearly show a hitherto unappreciated convergence of genes associated with both metals and wounds. Indeed, metal-linked genes are overrepresented in the most highly expressed genes across a range of biological tissues, including the skin.

RNA-Seq Profiling Reveals Temporal Alterations in Tissue Transcriptional Profiles Throughout Wound Repair

We performed detailed RNA-Seq profiling of gene expression change across normal murine healing, followed by in-depth analysis to explore the contribution of metals across the wound transcriptome (flow chart of approach provided in Supplementary Figure S2). RNA-Seq identified 10,294 differentially expressed genes (DEGs) between wounds and skin. Principal component analysis demonstrated clear clustering of individual samples by healing time-point, with some overlap between D3 and D7 (Figure 3A). Hierarchal clustering analysis of the top 250 most significant DEGs, enriched for keratins and chemokines (Supplementary Table S6), also showed high intra-group similarity (Supplementary Figure S3).

FIGURE 3. Metal-linked genes are considerably overrepresented in mouse and human wound transcriptomes. Murine wounds were collected at days 1, 3, 7 and 14 post-wounding for RNA-Sequencing. D0 = normal skin. n = 3 mice per group. Principal component (PC) analysis (A). Percentage (%) contribution. Venn diagrams show differentially expressed genes (DEGs) in wounds versus D0. Log₂ fold change ≥0.58 and Benjamini-Hochberg adjusted p value (Padjust) < 0.05 (B). DEGs linked to metals by functional annotation (C). Overrepresentation analysis demonstrates metal-linked (UniProt KB) genes are also injury-associated DEGs (D). Wound-upregulated DEGs further separated by metal and time-point in (E). Metal enrichment was confirmed in DEGs upregulated in human wounds via microarray (F–G). AW, acute wound; D, day post-injury; NS, non-significant.

Interrogation of DEG-level regulation revealed significant overlap in DEG lists between wound time points (compared to skin; Figure 3B). A total of 6,175 genes were upregulated in wounds (any time-point) compared to skin. Of these, 503 were specific to D1, 415 specific to D3, 238 specific to D7 and 1,670 specific to D14. Fewer genes were downregulated in wounds (4,288 at any time-point), with the highest number of genes again specific to the D14 comparison (980). To independently validate the RNA-Seq data, qRT-PCR was performed on a subset of upregulated and downregulated DEGs at each wound stage, which strongly correlated with RNA-Seq data (R > 0.94 and p < 0.001; Supplementary Figure S3).

Functional annotation was next performed on the full list of genes induced or repressed at each wound time-point, focusing on metallome-associated keywords and pathways (Figure 3C, Supplementary Table S7). At all wound time points, upregulated genes were highly annotated for calcium (UniProt KB: KW-0106), while enrichment for other metal-linked terms was wound stage specific. For example, in the D14 wound upregulated DEG group there was significant annotation for the terms metal-binding (UniProt KB: KW-0479; −Log₁₀ Padjust = 5.08), zinc (UniProt KB: KW-0862; −Log₁₀ Padjust = 2.11) and magnesium (UniProt KB: KW-0460; −Log₁₀ Padjust = 1.38). Downregulated genes were instead highly enriched for oxidoreductase (UniProt KB: KW-0560), metal-binding, iron, and magnesium from D3 onwards. In summary, our newly generated wound RNA-Seq data show, for the first time, a close association between metals and wound repair at the global transcription level.

Metal-Linked Genes Are Significantly Overrepresented in Murine and Human Wounds

In our RNA-Seq data set, 6,175 genes were upregulated upon injury (at any time-point), with 967 (15.7%) also present in the combined UniProt KB metal list (455 more genes than expected by chance, 1.89-fold overrepresentation, −log₁₀ p = 94.3; Figure 3D). Hierarchal clustering of the 250 most significant of these overlapping genes revealed four distinct sub-clusters with unique transcriptional signatures, annotated to specific functional processes (Supplementary Figure S4; Supplementary Tables S8, S9). Ingenuity Pathway Analysis (IPA) identified key upstream regulators for each sub-cluster (Supplementary Figure S4; Supplementary Tables S10, S11). As a key identified process was ECM, we extracted subsets of ECM-linked genes from UniProt and compared their temporal expression (RNA-Seq data) to metal profile (Supplementary Figure S5).

Of the 4,288 downregulated genes, 684 (16%) were also in the combined UniProt KB metal list (329 more genes than expected by chance, 1.92-fold over-representation, −log₁₀ p = 66.8; Figure 3D; Supplementary Table S12). Overrepresentation analysis was subsequently performed for wound upregulated DEGs for each individual metal (UniProt KB lists), displayed as a heat map (Figure 3E). Interestingly, manganese and calcium (particularly at D7) were the most significantly overrepresented metal-gene lists across healing, while zinc was overrepresented in wound DEGs only at D14 (Supplementary Table S13).

Due to the high overrepresentation of metal-linked genes in mouse wounds, we next interrogated a published human wound microarray data set comparing skin, acute wound (immediately after biopsy), and wounds at day 3 and 7 post-injury (GSE28914; Nuutila et al., 2012). Similar to the mouse data set, samples clustered by time point, with the highest similarity in expression profiles observed between D3 and D7 (Supplementary Figure S6). We then assessed metal gene enrichment (using the combined human UniProt KB metal lists), where of the genes that were upregulated in human wounds, 827 were also found in the combined metal gene list (398 more genes than expected by chance, 1.93-fold overrepresentation, −log₁₀ p = 85.3; Figure 3F; Supplementary Table S14). DEGs downregulated in human wounds were also significantly enriched for metals (74 more genes than expected by chance, 2.12-fold overrepresentation, −log₁₀ p = 17.4; Figure 3F). Additionally, we performed overrepresentation analysis on human wound-upregulated DEGs separated by wound time point and metal (Figure 3G; Supplementary Table S15). Taken together, these data clearly highlight metal-specific wound gene enrichment profiles that are temporally altered in mouse and human wounds.

DEGs Cluster Into Wound-Relevant Subsets That Are Also Highly Enriched for Metals

Wound-upregulated DEGs from RNA-Seq analysis (combined list from all wound time points) were next independently clustered (Markov Cluster Algorithm) based on expression signature, and each cluster assigned a dominant wound-linked functional annotation (Figure 4A; Supplementary Table S16). The identified upregulated gene clusters were linked to inflammation, proliferation, keratinisation and ECM. Overrepresentation analysis, displayed as a heat map (Figure 4A), demonstrated specific and significant enrichment for different metal-linked gene groups (UniProt KB) across 9 out of 10 identified wound-upregulated DEG clusters (Supplementary Table S17). Normalised count data were then used to generate mean cluster expression profiles (Figure 4B). For example, genes in cluster four are overrepresented for both calcium and zinc, have a profile of expression that peaks at D14, and are associated with keratinisation. Hence, for each cluster we show linked function, metal association (Figure 4A) and expression profile (Figure 4B). Next, we extracted subsets of genes from the largest of the assessed UniProt metal lists (calcium, magnesium, iron and zinc) that were associated with protein binding (GO:0005515), transporter activity (GO:0005215) and signalling receptor activity (GO: 0038023). The protein binding GO group was the only group to feature a high proportion of genes from all four metal lists (Supplementary Table S18). Interestingly, cluster-specific overrepresentation was altered for metal protein binding genes versus all metal genes (Figure 4C; Supplementary Tables S18, S19). For example, calcium- and zinc-linked genes were no longer significantly overrepresented in cluster 4 (keratinisation) while zinc became overrepresented in cluster 6 (ECM) when focusing on protein binding subsets. Stepping down a level of GO annotation revealed that a high proportion of protein binding metal-linked genes were enzyme and signalling receptor binding (Figure 4D).

FIGURE 4. Wound-induced genes cluster by biological function and show process-specific metal overrepresentation. Wound upregulated differentially expressed genes (DEGs) were clustered (A) by similarity in expression signatures (VST data) using a Markov clustering algorithm. Clusters named based on functional annotation. Overrepresentation of UniProt Knowledge Base metal genes in defined clusters (heat map, right). Expression profile for each cluster shown in (B). Overrepresentation performed on metal-linked genes from the protein binding gene ontology (GO) group (GO:0005515) in (C), with specific subset contribution in (D). NS, non-significant; None, no genes in metal group; VST, variance stabilised transformation; R, Pearson’s correlation.

Calcium Shows Temporal Wound Level Regulation That Converges on Key Time-Point Specific Upstream Regulators

Given their high representation in wound DEGs (Figures 3, 4), we explored calcium-specific genes (i.e., only genes found in the calcium UniProt KB group) within our murine RNA-Seq data set and the published human microarray data set in more detail. See Supplementary Figure S6 for volcano plots of calcium-specific changes in wound gene expression at each healing time-point (versus D0; Supplementary Table S20), with corresponding time-point dependent functional annotation (Supplementary Figure S7, Supplementary Tables S21, S22).

Focusing on the 250 most significant calcium-specific genes (filtered by fold change <1.5 and significance) in both mouse and human data sets, hierarchical clustering revealed greatest expression similarity between D3 and D7 wounds (Figure 5A, Supplementary Tables S23, S24). Note, D3 and D7 were also the time points where murine wound tissue calcium levels were strongly increased (ICP-MS; Figure 1A). Hierarchical clustering at the gene level revealed four distinct clusters for both heat maps (labelled 1–4; Figure 5A). Functional annotation analysis showed that calcium-linked genes in cluster 3 (high expression at D3 and D7 in mouse and human wounds) were strongly associated with a range of wound processes, including ECM (−log₁₀ p = 18.5 in mouse; −log₁₀ p = 6.8 in human) and cell adhesion (−log₁₀ p = 4.9 in mouse; −log₁₀ p = 5 in human) (Figure 5B; Supplementary Tables S25, S26). Cluster 4 (high expression at D1, D3 and D7 in mouse and D3 and D7 in human) contained calcium-linked genes involved in collagen degradation (−log₁₀ p = 11.7 in mouse; −log₁₀ p = 5 in human) and immunity (−log₁₀ p = 2.5 in mouse; −log₁₀ p = 8 in human).

FIGURE 5. Transcriptional regulation of calcium is wound stage specific and converges on key upstream regulators. Hierarchal clustering of the top 250 most significant genes in mouse RNA-Seq and human microarray filtered from the calcium UniProt knowledge base (KB) list for each species [KW-0106; (A)]. D0 = normal skin. D1 = day 1 post-wounding. Cut-off = Benjamini-Hochberg adjusted p value (Padjust) < 0.05 and log₂ fold change = 0.58. Functional annotation of identified clusters and their expression profiles (B). VST = variance stabilised transformation. Network analysis identifying TNF as a top transcriptional regulator for cluster 4 (C). m, mouse; h, human.

Next, IPA analysis was performed on RNA-Seq and microarray data to determine putative upstream regulators for each calcium-specific cluster. Interestingly, multiple potential upstream regulators, common to both mouse and human, were identified for clusters 1, 3 and 4 (Supplementary Tables S27, S28). These included SMARCA4 (cluster 1), TGFB1, alpha catenin, VEGFA (cluster 3) and TNF, EGF and IL6 (cluster 4). TGFB1, the top identified upstream regulator for cluster 3 in both mouse (activation score: 3.562, p = 2.96E-08) and human (activation score: 2.531, p = 3.19E-05), was linked to gene networks that contained multiple ECM components and ECM remodeling enzymes. TNF, one of the top predicted activators for cluster 4 for both mouse (activation score: 4.2, p = 6.83E-06) and human (activation score: 3.7, p = 5.31E-08), was used to construct networks dominated by epidermal and immune-linked genes (Figure 5C; Supplementary Tables S29, S30). Collectively, these data reveal high transcriptional regulation of calcium in wounds, particularly at D3 and D7 post-injury, suggesting that calcium may have new unappreciated roles in later stages of healing.

Metal-Led Transcriptomic Analysis Identifies TNF as an Important Regulator of Calcium-Induced Epidermal Differentiation

Finally, we investigated the link between calcium and one of the top predicted upstream regulators (TNF) in more detail. We first performed immunofluorescence on mouse wound tissue to determine the temporal expression profile of TNFα in the epidermis and dermis throughout healing (Figure 6A). TNFα expression peaked at D1 in wound edge epidermis (p < 0.001; Figure 6B), while the thickness of the wound-edge neo-epidermis peaked at D7 (p < 0.001; Figure 6C). In the dermis, TNFα levels were significantly elevated at D1 and D3 versus D0 (both p < 0.001; Figure 6D), with expression correlating closely with macrophage numbers at these time points (Figure 6E).

FIGURE 6. TNFα inhibitor treatment impairs calcium-dependent differentiation in keratinocytes. Mouse epidermal staining for keratin 14 (K14, red) and TNFα (green) and dermal staining for CD107b (green) and TNFα [red; (A)]. DAPI = blue nuclei. Bar = 50 µm. Quantification of epidermal TNFα (B), epidermal thickness (C), dermal TNFα (D) and macrophages [CD107b+ve (E)]. n = 5 mice. qRT-PCR for differentiation markers (F–H) and TNF-linked markers (I–K) in primary murine and human keratinocytes. C87 = TNFα inhibitor. n = 3 donors. Mean +SEM. One-way ANOVA with Tukey’s post-hoc analyses were performed. */# = p < 0.05, ** = p < 0.01 and ***/### = p < 0.001. Asterisk alone versus D0 or no calcium. Hashtag alone versus D1 in (B–E).

We next focused on keratinocyte differentiation as a process of interest, based on our transcriptomic analysis, tissue staining, and the widely documented association between calcium and keratinocyte differentiation (Elias et al., 2002). Here, the influence of TNFα on calcium-induced keratinocyte differentiation was assessed in primary MEKs and HEKs. As expected, high calcium (2 mM CaCl₂) significantly increased the expression of epidermal differentiation markers in mouse and human keratinocytes. However, this increase was strongly attenuated by co-treatment with the TNFα inhibitor, C87 (Figures 6F–H). Interestingly, treatment with C87 led to significant upregulation of Mmp1/MMP1 (Figure 6I) and downregulation of Plscr1/PLSCR1 (Figure 6J), genes predicted to be regulated by TNF in calcium cluster 4. The expression of Tnf/TNF itself was significantly increased following calcium treatment (p < 0.01—p < 0.001), an effect which was reversed by C87 (Figure 6K). Thus, these data suggest an important wound-relevant functional linkage between TNF and calcium in keratinocytes, validating our global transcriptional interrogation approach.

Discussion

This study is the first to comprehensively evaluate changes in multiple endogenous metals across a murine healing time course and link these temporal metal fluctuations to the global wound transcriptome. Metallomics is an emerging, but still relatively understudied, area of biology. Recent reports have linked specific metals to pathologies, including ageing (Möller et al., 2019) and diabetes (Fernández-Cao et al., 2019). However, understanding of how these endogenous bio-metals contribute to normal function and disease states remains limited. In the skin, almost all metallome research to date has focussed on measuring metal bioaccumulation following treatment or occupational exposure (Lidén et al., 2008; Al-Dayel et al., 2011; Bianco et al., 2014), rather than matching endogenous metal levels to function in normal skin biology.

Data presented in the current study reveal metal-specific, time-dependent variations in calcium, magnesium, iron, zinc, aluminium, copper and manganese across murine skin healing. Crucially, each wound metal fluctuated in a manner that provides an indication into their role(s) in cellular aspects of healing. For example, ICP-MS-measured calcium levels peaked at D3 and D7 (the inflammatory and proliferative stages of healing). Note, the majority of existing wound literature focuses on measuring calcium flashes or waves induced within minutes of injury (Xu and Chisholm, 2011; Yoo et al., 2012; Razzell et al., 2013; Tu et al., 2019). A tissue calcium peak at D3/D7 is much more in line with the tissue-relevant roles for calcium in regulating immune cell function and cellular proliferation (Schwarz et al., 2006; Vashi et al., 2017). Intriguingly, we report a substantial reduction in tissue calcium levels during remodelling (D14), a stage characterised by a reduction in cell infiltration and the cessation of many calcium-dependent processes (e.g., MMP production, Singh et al., 2015). Indeed, calcium-specific genes upregulated in D14 wounds were not associated with any upstream regulators, proposing they have no obvious wound-relevant cellular functions.

Copper significantly accumulates at day 7 PW, in line with its reported role as an essential cofactor in angiogenesis (Urso and Maffia, 2015; Das et al., 2016), where high numbers of newly formed blood vessels are observed (Hattori et al., 2009). Copper is also an essential cofactor for lysyl oxidase (Lucero and Kagan, 2006), a collagen stabilising enzyme (Huang et al., 2019), suggesting a role in early granulation tissue maturation. The transition metals iron, aluminium and manganese were all elevated in the later stages of healing (D7 and D14 PW). This pattern of temporal induction fits with the ECM deposition and remodelling stages of wound repair. In fact, we recently demonstrated that wound tissue iron is required for ECM deposition and remodelling, and significantly perturbed in diabetic healing (Wilkinson et al., 2019b; 2019c). Moreover, Craven et al. (2001) showed that overexpression of manganese superoxide dismutase led to collagen accumulation, while Treiber et al. (2009) demonstrated that manganese promoted the contraction of collagen gel lattices, a proxy for in vivo wound closure.

Interestingly, we did not observe significant alterations in the overall abundance of zinc throughout skin healing, despite its crucial roles in many wound-associated processes (e.g., metalloproteinase activity, Caley et al., 2015). This could indicate that the high basal levels of zinc observed in skin are sufficient to support wound healing requirements. Our results remain in line with previous studies in rat wounds (Coger et al., 2019; Lansdown et al., 1999), but crucially highlight a caveat of measuring total metals in digested tissue samples without assessing the state of existence or intracellular localisation of metals. Indeed, our tissue extraction process and measurement method was not compatible with preserving free versus bound metals, but a number of approaches could be utilised to gain this valuable information. To spatially differentiate metals in situ, laser-ablation ICP-MS (Cruz-Alonso et al., 2019) or synchrotron radiation x-ray fluorescence (Zhang et al., 2018) can be used. Moreover, elemental analysis alone does not allow assessment of free versus bound metal but combining ICP-MS with sophisticated proteomic approaches would permit identification of metal-containing proteins (Wang et al., 2019). While technologically challenging, future studies should thus focus on characterising the spatial localisation and biological availability of metals in tissues.

It is currently unclear how observed murine changes in wound metals will translate into human healing. However, we note that basal metal levels are comparable between human and murine skin (unpublished data), with metal transport proteins generally conserved between both species (Danks, 1986). Our transcriptional analysis also reveals key parallels in metal-linked genes and functions in mice and humans. Furthermore, our data focussed on acute wounds, yet a major future consideration will be to understand how metals may be altered in/or contribute to chronic wound healing. Indeed, the influence of excessive iron to chronic wound inflammation has been assessed previously (Sindrilaru et al., 2011). Metal composition may also modulate bacterial presence and pathogenicity (Ma et al., 2015), highlighting potential microbiome-metallome interactions as an attractive avenue for further exploration.

We took a global transcriptional approach to investigate the relationship between metals and the transcriptome across murine (our RNA-Seq data) and human acute wound healing (published microarray). The RNA-Seq data alone provides an important new resource for the research community, as there are no comprehensive publicly available murine wound time course RNA-Seq data sets. Our detailed interrogation of wound DEGs delivers vital new insight into the global ties between the metallome and transcriptome, revealing numerous metal-linked genes that are induced or repressed in a temporal manner. Markov clustering of wound-induced genes independently confirmed that metal enrichment was dependent on biological function. For example, calcium was strongly linked to proliferation, while both calcium and zinc, metals widely associated with epidermal differentiation (Elias et al., 2002; Inoue et al., 2014), were linked to keratinisation. Focusing only on the metal-linked genes that were associated with protein binding provided further information into their relevance in key wound-linked biological processes. Indeed, these data highlight the crucial importance of investigating associations between metals and biological function but are limited by virtue of assessing only transcriptional changes, which do not always correlate with protein levels and function (Takemon et al., 2021). Proteomic approaches should therefore be incorporated in future studies to fully elucidate the links between metals and the biology of wound repair.

IPA network analysis identified the wound-associated transcription factor, TNF, important in mediating inflammation in early repair (Ritsu et al., 2017), as a primary upstream regulator for calcium-linked gene clusters in mouse and human wounds. Finally, we performed cell-based studies to confirm a putative role for TNFα in calcium-induced keratinocyte differentiation. Here we demonstrated that inhibition of TNFα attenuated calcium-induced differentiation in mouse and human keratinocytes. It is important to note that previous studies show a role for TNFα in negatively augmenting the epidermal barrier in skin pathology (Basile et al., 2001; Marble et al., 2007; Danso et al., 2014). However, to our knowledge, this is the first time endogenous TNFα has been directly linked to calcium-induced epidermal differentiation in a normal reparative context in mice and humans.

Conclusion

In summary, our data provide the first comprehensive documentation of the metallomic/transcriptomic landscape of a healing wound, suggesting exciting new roles for metals in skin repair (overview provided in Figure 7). We note that previous studies have evaluated only a single metal (iron) in the context of healing pathology; which was shown to be perturbed in both diabetic murine and human wounds (Sindrilaru et al., 2011; Wilkinson et al., 2019c). Thus, our study not only reveals previously unappreciated and potentially wide-reaching roles for metals in wound repair, it also paves-the-way for essential follow on work to explore the cellular consequences of metallome dysbiosis during injury and disease in skin and other tissues.

FIGURE 7. An overview of the potential roles for metals in wound repair. Schematic outlining the key metal associations identified in this work, for specific wound healing stages. Associations are based on metal gene list overrepresentation and functional annotation analysis of genes differentially expressed during murine wound healing.

Data Availability Statement

The data that support the findings of this study are available as Supplementary Material or from the corresponding author upon reasonable request. RNA-Seq data are deposited in the Sequence Read Archive database (https://www.ncbi.nlm.nih.gov/sra) under BioProject accession number PRJNA686364.

Ethics Statement

The studies involving human participants were reviewed and approved by the South Central Oxford C Research Ethics Committee. The patients/participants provided their written informed consent to participate in this study. The animal study was reviewed and approved by the United Kingdom Home Office.

Author Contributions

HW and MH conceived and designed the study. HW collected and analysed the data and wrote the manuscript. B-AG and MH provided critical review of the manuscript and supervised the study. All authors reviewed and approved the final manuscript.

Funding

An MRC (United Kingdom) CASE PhD Studentship supported this work (grant number: MR/M016307/1).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

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.

Acknowledgments

We would like to thank Mr. Paul Lythgoe, Dr. Karen Theis and Dr. Kimberly Mace (University of Manchester), and Professor Tim Palmer (Hull York Medical School), for sharing of equipment and expertise. We would also like to thank Dr. Dave Lunt (University of Hull) for bioinformatics advice and Mr. Paolo Matteucci and the surgical team at Castle Hill hospital for access to human skin.

Supplementary Material

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

References