Comparative transcriptomic profiling of hydrogen peroxide signaling networks in zebrafish and human keratinocytes: Implications toward conservation, migration and wound healing

Abstract

Skin wounds need to be repaired rapidly after injury to restore proper skin barrier function. Hydrogen peroxide (H2O2) is a conserved signaling factor that has been shown to promote a variety of skin wound repair processes, including immune cell migration, angiogenesis and sensory axon repair. Despite growing research on H2O2 functions in wound repair, the downstream signaling pathways activated by this reactive oxygen species in the context of injury remain largely unknown. The goal of this study was to provide a comprehensive analysis of gene expression changes in the epidermis upon exposure to H2O2 concentrations known to promote wound repair. Comparative transcriptome analysis using RNA-seq data from larval zebrafish and previously reported microarray data from a human epidermal keratinocyte line shows that H2O2 activates conserved cell migration, adhesion, cytoprotective and anti-apoptotic programs in both zebrafish and human keratinocytes. Further assessment of expression characteristics and signaling pathways revealed the activation of three major H2O2-dependent pathways, EGF, FOXO1, and IKKα. This study expands on our current understanding of the clinical potential of low-level H2O2 for the promotion of epidermal wound repair and provides potential candidates in the treatment of wound healing deficits.

Figure 1. Whole transcriptome RNA-seq profile of larval zebrafish in response to low H₂O₂ treatment.

(a) Pools of ~500 larvae/set of 4 day-post-fertilized (dpf) zebrafish larvae were treated with 0.01% (3 mM) H₂O₂ for three hours and total RNA was subsequently collected followed by pair-end next generation RNA sequencing (n = 3 biological replicates). (b) H₂O₂ sensor (HPF) either alone or with H₂O₂ treatment shows that H₂O₂ is mostly retained in the skin epithelium (n = 5 fish). (c) Distribution of mapped reads in the zebrafish transcriptome. RNA-seq data sent to NCBI (GEO: GSE75728). (d) Transcript complexity between untreated and H₂O₂-treated larval zebrafish samples. Based on read CPM (counts per million), the left-most value on the X-axis represents the most highly expressed transcripts, which is incrementally summed with each successively lower expressed transcript (rightward). The y-axis (% contribution of the total transcripts) was calculated using: [CPM/sum of all CPM] x 100%. (e) RNA-seq data was normalized with the read CPM method of the number of mapped reads on gene exons. Transcript expression data transformed on M (log ratio of fold change) and A (mean average) scale. Boxed blue regions represent statistically significant transcripts (p < 0.05) returned by the test for differential expression. The MA-plot shows the log₂ fold changes from the treatment over the mean of normalized counts, i.e. the average of counts normalized by size factor. Cutoff set to >1 log₂ CPM averaged over all samples, and below the cutoff there is no real inferential power. Note: Statistical significance drops below the threshold. (f) Quantitative PCR validation of RNA-seq results of a sub-set of candidate targets. Full data set is presented in Table S7 (n = 3 biological replicates). (g) Heat map indicates unsupervised hierarchical clustering of the top (left) and bottom (right) most significantly enriched transcripts derived from the RNA-seq data after H₂O₂ treatment. Hierarchical clustering was performed between individual experiments and transcripts. The color key indicates the log₂ CPM expression values. Abbreviations: HPF (hydrogen peroxide fluorogenic probe or pentafluorobenzenesulfonyl-fluorescein), CPM: Counts per million, FC: Fold-change.

Figure 2. Transgenic NF-κB:EGFP zebrafish reveal peripheral NF-κB activation after H₂O₂ treatment Spatial patterns of increased NF-κB activation after H₂O₂ treatment of zebrafish larvae.

(a) 3-dpf Tg(NF-κB:EGFP) reporter zebrafish larvae at the start of the experiment. (b) Tg(NF-κB:EGFP) zebrafish were imaged after 2hr post 0.01% H₂O₂ treatment. Both whole larvae (upper) and tail fins (below) were imaged, whereby GFP⁺ cells were partially localized to the periphery (arrows). (c) Quantitative analysis using relative mean fluorescence of the z-stack projected images using ImageJ (n = 2, 10 fish). Observation of increased number of GFP-labeled cells and overall fluorescence intensity in the whole larvae. (d) Higher-resolution analyses of the tail fin revealed a peripheral tissue spatial pattern of increase NF-κB activated cells after H₂O₂ treatment. n = 3 independent experiments per condition.

Figure 3. Enriched disease and biological functions and classification of differentially regulated transcripts after H₂O₂ treatment.

(a) Downstream Effects Analysis in the Ingenuity’s Pathway Analysis was used to visualize, via color-coded heatmaps, putative biological and disease trends in H₂O₂-treated zebrafish larvae. Within the cell movement category (boxed in red) are 24 differentially expressed genes (q < 0.1, false discovery rate). The color intensity of the squares in the heatmaps reflects the strength of the absolute z-score for predictions (orange = positive, blue = negative). The categories are assembled with the most significant p-values displayed on the left of the heatmap. The size of the squares reflects the z-score values. (b) Protein classification of transcripts affected by H₂O₂ treatment of larval zebrafish. Bars represent the % of mapped transcripts to appropriate annotations. Absolute gene numbers are shown in parentheses. (c) Shown is the % of mapped transcripts up (red) or down (blue) regulated for three different categories: oxidoreductase, cell adhesion and ECM. (d) Classification of the molecular functions of affected transcripts after H₂O₂ treatment. Graphs represent the % of mapped transcripts to appropriate annotations. Analyses were performed using PANTHER.

Figure 4. Upstream pathway analysis of larval zebrafish treated with H₂O₂.

(a) Causal upstream networks determined using Ingenuity’s Upstream Regulatory Analysis after H₂O₂ treatment based on the literature compiled in the Ingenuity^® Knowledge Base. Fisher’s exact test p-values were calculated to assess the significance of enrichment of the RNA-seq data for the genes downstream of an upstream regulator (Worksheet 4). (b) Most significant downstream genes within the H₂O₂ upstream network in zebrafish samples (p < 0.05, n = 3 biological replicates). Top numerical value for each transcript represents log₂(fold change) (H₂O₂ vs. untreated), while the lower value represents the p-value. Shades of red indicate the degree of upregulation, while shades of green represent the degree of downregulation. The edges connecting the nodes are colored orange when leading to activation of the downstream node, and yellow if the findings underlying the relationship are inconsistent with the state of the downstream node. Pointed arrowheads indicate that the downstream node is expected to be activated, while blunt arrowheads indicate that the downstream node is expected to be inhibited. (c) Overlap among the differentially expressed genes in H₂O₂-treated zebrafish and curated chemical-gene interactions for H₂O₂ derived from the Comparative Toxicogenomics Database (http://ctdbase.org). (d) Functional annotation analysis using DAVID of the H₂O₂-downstream genes. Enrichment scores ≥1 are considered significant.

Figure 5. Comparison of ARE/EPRE:GFP activation in zebrafish.

(a) The caudal fin of an uninjured EPRE:GFP larval zebrafish was imaged over the course of 12 hours. First and last images of the time-lapse sequence are shown. (b) Matching surface plots and quantification, comparing the fluorescence means of 4 individual fish (n = 4). Statistical significance was tested between first and last time points, showing lack of EPRE:GFP activation by 12 hours. (c) First and last image of a time-lapse sequence showing the amputated caudal fin (arrows) of an EPRE:GFP larval zebrafish. (d) Matching surface plots and quantification, comparing the fluorescence means of 6 individual fish, show that injury fails to activate EPRE:GFP.

Figure 6. Comparative whole transcriptomic analysis between human epidermal keratinocytes and zebrafish larvae treated with H₂O₂.

(a) Comparison of H₂O₂-mediated gene expression between HaCaT cells and zebrafish larvae reveals 23 congruent and 18 non-congruent overlapping genes. (b) Distribution of the overlapping genes based on logFC values. Congruent transcripts are boxed in gray. (c) Functional annotation clustering of the congruent genes using DAVID reveals biological processes, which are conserved between human epidermal cells and zebrafish treated with H₂O₂. (d) Detailed analysis of the HaCaT response to H₂O₂ [100 μM] over the entire time course within the GSE46343 study. Data represented as the % gene expression of treated vs. untreated samples. (e) qPCR analyses of H₂O₂-treated HEK01 (6hr) consistent with activation of epithelial cell migration and adhesion. One-way ANOVA at an alpha = 0.05 (95% confidence interval) and Tukey’s multiple comparison post-tests were utilized. Significance is denoted with asterisks: *p < 0.05, n = 3–5 experiments. Abbreviations: Fold change (FC), beta actin (ATCB), involucrin (IVL), integrin beta 4 (ITGB4), insulin-like growth factor-binding protein 1 (IGFBP1), heat shock 70 kDa protein 1 L (HSPA1L), not significant (ns).

Figure 7. Inhibition of IKK/Ikk delays scratch closure and intracellularly accumulates within injured epidermal keratinocytes.

(a) IKK is necessary for H₂O₂-induced HEK01 keratinocyte migration after wounding. Scratch assays were performed with H₂O₂ (0.1 μM) and Wedelolactone (50 μM) using HEK01 cells. Two-way ANOVA at an alpha = 0.05 (95% confidence interval) and Bonferroni’s multiple comparison post-tests were utilized. Significance is denoted with asterisks: *p < 0.05, **p < 0.01 (n ≥ 3–5 cell culture experiments). (b) Rapid H₂O₂ production using 1 μM HPF (hydrogen peroxide fluorogenic probe) at the scratch (sc) margin of HEK01 keratinocytes within 30 minutes. Bar = 100 μm (c) Rapid subcellular accumulation of IKKα within injured HEK01 cells at the scratch (sc) margin (white arrows) after 30 minutes compared to unscratched cells. Orthogonal views (red arrows) of an injured cell show peri-nuclear and cytoplasmic distribution and accumulation of IKKα. Bar = 20 μm. (d) Schema of overall findings.