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.

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.

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).

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 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-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.

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.

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.

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.

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).

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.

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.

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.

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).

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.

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.

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).

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).

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.

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).

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.