Integrative small and long RNA omics analysis of human healing and nonhealing wounds discovers cooperating microRNAs as therapeutic targets
Figures

Paired profiling of miRNA and mRNA expression in human wounds.
(a) Schematic of analysis in this study (n = samples used for RNA-seq and qRT-PCR validation in full-thickness tissues). (b) Principal component analysis (PCA) plots based on miRNA (upper panel) and mRNA (lower panel) expression profiles. Each dot indicates an individual sample. The numbers of differentially expressed (DE) miRNAs (c) and mRNAs (e) in venous ulcer (VU) (n = 5) compared to the Skin, Wound1, and Wound7 from five healthy donors are shown in Venn diagrams. False discovery rate (FDR) <0.05, fold change ≥2 for miRNAs, and ≥1.5 for mRNAs. (d) The heatmap depicts the 32 miRNAs specifically dysregulated in the VU with scaled expression values (Z-scores). Weighted gene coexpression network analysis (WGCNA) modules of each miRNA belongs to are marked with color bars. The miRNAs with experimentally validated expression changes are highlighted in red.
-
Figure 1—source data 1
miRNAs and mRNA with expression change specifically in venous ulcers.
- https://cdn.elifesciences.org/articles/80322/elife-80322-fig1-data1-v2.xlsx

Weighted gene coexpression network analysis (WGCNA) of miRNAs and mRNAs in wound healing.
(a) Cluster dendrogram shows miRNA coexpression modules: each branch corresponds to a module, and each leaf indicates a single miRNA. Color bars below show the module assignment (the first row) and Pearson’s correlation coefficients between miRNA expression and the sample groups (the second to the fifth row: red and blue lines represent positive and negative correlations, respectively). (b) Heatmap shows Pearson’s correlations between miRNA module eigengenes (MEs) and the sample groups. The correlation coefficients and the adjusted p values (false discovery rate, FDR) are shown where the FDRs are less than 0.05. For the venous ulcer (VU)-associated modules m8 (c) and m9 (d), bar plots (left) depict the ME values across the 20 samples analyzed by RNA-seq, and network plots (right) show the top 20 miRNAs with the highest kME values in each module. Node size and edge thickness are proportional to the kME values and the weighted correlations between two connected miRNAs, respectively. The miRs with their targetome enriched with VU-mRNA signature (see Figure 4b) are highlighted in red. Transcription factors (TFs) with their targets enriched in the m9 module (Fisher’s exact test: false discovery rate [FDR] <0.05) are listed below the network, and TFs differentially expressed in VU are underlined. (e) Heatmap shows Pearson’s correlations between mRNA MEs and the sample groups. (f) The gene expression pattern of each module across all the sample groups is depicted with line charts. Gene ontology analysis of biological processes enriched in each module is shown at the right.
-
Figure 2—source data 1
Weighted gene coexpression network analysis of miRNAs in wound healing.
- https://cdn.elifesciences.org/articles/80322/elife-80322-fig2-data1-v2.xlsx
-
Figure 2—source data 2
Top 20 driver miRNAs of each significant module in weighted gene coexpression network analysis (WGCNA).
- https://cdn.elifesciences.org/articles/80322/elife-80322-fig2-data2-v2.xlsx
-
Figure 2—source data 3
Transcription factors (TFs) regulating miRNA expression in each module.
- https://cdn.elifesciences.org/articles/80322/elife-80322-fig2-data3-v2.xlsx
-
Figure 2—source data 4
Weighted gene coexpression network analysis of mRNAs in wound healing.
- https://cdn.elifesciences.org/articles/80322/elife-80322-fig2-data4-v2.xlsx
-
Figure 2—source data 5
Transcription factors (TFs) with targets enriched in significant mRNA modules.
- https://cdn.elifesciences.org/articles/80322/elife-80322-fig2-data5-v2.xlsx

Preservation analysis of miRNA and mRNA coexpression modules.
The preservation analysis for miRNAs (a) and mRNAs (b) modules. The standardized Z-scores were calculated for each module by permutating 200 times using the 20 samples analyzed by RNA-seq as reference and test datasets. Modular preservation is strong if Z-summary >10, weak to moderate if 2 < Z-summary < 10, no evidence of preservation if Z-summary ≤2. miRNA modules.

Weighted gene coexpression network analysis (WGCNA) of miRNAs in human skin wound healing (related to Figure 2a–d).
(a) The scale-free topology fit index (left) and mean connectivity (right) for various soft-threshold powers. The red line indicates signed R2 = 0.8. (b–i) Bar plots (left) depict the module eigengene (ME) values across the 20 samples analyzed by RNA-seq and network plots (right) show the top 20 miRNAs with the highest kME values in each module. Node size and edge thickness are proportional to the kME values and the weighted correlations between two connected miRNAs, respectively. Transcription factors (TFs) with their targets enriched in the modules (Fisher’s exact test: false discovery rate [FDR] <0.05) are listed below the networks.

Transcription factors (TFs) with targets enriched in the miRNA m9 module (related to Figure 2d).
(a) Pearson’s correlations between the first principal component (PC1) of miRNAs module m9 and the expression of GATA3 or KLF4. (b) GATA3 and KLF4 expression in the skin, day 1 and 7 acute wounds from five healthy donors and in five venous ulcers analyzed by RNA sequencing. Mann–Whitney t-test, **p < 0.01.

Weighted gene coexpression network analysis (WGCNA) of mRNAs in human skin wound healing (related to Figure 2e).
(a) The scale-free topology fit index (left) and mean connectivity (right) for various soft-threshold powers. The red line indicates signed R2 = 0.8. (b) Cluster dendrogram shows mRNA coexpression modules: each branch corresponds to a module, and each leaf indicates a single miRNA. Color bars below show the module assignment (the first row) and Pearson’s correlation coefficients between mRNA expression and the sample groups (the second to the fifth row: red and blue lines represent positive and negative correlations, respectively). (c) Network plots of mRNA modules: the top 20 mRNAs with the highest kME values in each module were plotted in the networks. Node size and edge thickness are proportional to the kME values and the weighted correlations between two connected mRNAs, respectively. Transcription factors (TFs) with their targets enriched in the modules (Fisher’s exact test: false discovery rate [FDR] <0.05) are listed below the networks.

Functional enrichment analysis for mRNA modules (related to Figure 2f).
(a) The top 10 gene ontology (GO) terms with false discovery rate (FDR) less than 0.05 are shown for each module. (b) MetaCore analysis based on a curated database identified the process networks (left) and pathway maps (right) enriched in the M8 or the combined M11.M12 modules, respectively.

Correlation analysis between miRNA and mRNA expression changes in venous ulcer (VU).
(a–c) Correlations between the first principal component (PC1) of VU-associated differentially expressed (DE) miRNAs, and the PC1 of VU DE mRNAs predicted as miRNA targets. (d) PC1 correlations between the hub miRNAs and their predicted targets in the VU-associated miRNA and mRNA modules. Pearson’s correlation coefficients (r) and p values are shown. The mRNA networks are plotted with the top 20 most connected module genes. Transcription factors (TFs) with targets enriched in the VU-associated modules are listed below the networks.

Integrative analysis of miRNA and mRNA expression changes in venous ulcer (VU).
Heatmaps show the enrichment for VU-affected mRNAs and mRNA modules in the top targets of (a) VU-associated differentially expressed (DE) miRNAs and miR modules and (b) the individual candidate miRNAs. Odds ratio (OR) and p values are indicated when OR >1 and p value <0.05 (Fisher’s exact test). The miRNAs with experimentally validated expression changes are highlighted in bold. (c) miR-mediated gene regulatory networks in VU are constructed with the miRNAs in (b), the mRNAs predicted as the strongest targets and with anticorrelated expression patterns (Pearson’s correlation, p value <0.05 and r < 0) with these miRNAs in human wounds, and the transcription factors (TFs) regulating these miRNAs' expression from the TransmiR v2 database. An enlarged version of these networks can be found in Figure 4—figure supplements 1 and 2.
-
Figure 4—source data 1
Gene set enrichment analysis for VU-affected DE mRNAs and mRNA modules in the strongest targets of VU-associated DE miRNAs and miRNA modules.
- https://cdn.elifesciences.org/articles/80322/elife-80322-fig4-data1-v2.xlsx
-
Figure 4—source data 2
Individual candidate miRNAs with their targets enriched for the VU mRNA signature.
- https://cdn.elifesciences.org/articles/80322/elife-80322-fig4-data2-v2.xlsx

MiR-mediated gene regulatory network in VU, which is constructed with the upregulated miRNAs in Figure 4b, the mRNAs predicted as the strongest targets and with anticorrelated expression patterns (Pearson’s correlation, p value <0.05 and r < 0) with these miRNAs in human wounds, and the transcription factors (TFs) regulating these miRNAs’ expression from the TransmiR v2 database (related to Figure 4b, c and Figure 4—source data 2).
The mRNAs cotargeted by miRNAs with unrelated seed sequences are highlighted in yellow.

MiR-mediated gene regulatory network in VU, which is constructed with the downregulated miRNAs in Figure 4b, the mRNAs predicted as the strongest targets and with anticorrelated expression patterns (Pearson’s correlation, p value <0.05 and r < 0) with these miRNAs in human wounds, and the transcription factors (TFs) regulating these miRNAs’ expression from the TransmiR v2 database (related to Figure 4b, c and Figure 4—source data 2).
The mRNAs cotargeted by miRNAs with unrelated seed sequences are highlighted in yellow.

Experimental validation of miRNAs’ expression in human skin wounds.
(a) A schematic diagram of experimental validation in this study. (b–j) qRT-PCR analysis of venous ulcer (VU)-associated differentially expressed (DE) miRNAs in the skin, day 1 and 7 acute wounds from seven healthy donors and VU from 12 patients. (k) VU-associated DE miRNA expression in RNA-seq of epidermal keratinocytes from human skin and wound day 7. (l) Representative photographs of laser capture microdissection (LCM) -isolated epidermis from skin, wound7 and VU. qRT-PCR analysis of (m) miR-34a-5p and (n) miR-424-5p in LCM-isolated epidermis (n = 7). The data are presented as mean ± standard deviation (SD). Wilcoxon signed-rank test was used for the comparison between Skin, Wound1, and Wound7; Mann–Whitney U-test was used for comparing VU with the skin and acute wounds. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001.
-
Figure 5—source data 1
Experimental validation of miRNAs’ expression in human skin wounds (related to Figure 5b–j), or in LCM-isolated epidermis (related to Figure 5m, n), enrichment analysis of the experimentally validated miRNA targets for the venous ulcer (VU) gene signature (related to Figure 6k), miRNAs’ overexpression efficiency (related to Figure 6—figure supplement 1), and miRNA targets validated by the microarray and in VU gene signature (related to Figure 7a and b, lower panels).
- https://cdn.elifesciences.org/articles/80322/elife-80322-fig5-data1-v2.xlsx

RNA-sequencing results for the miRNAs selected for experimental validation.
Mann–Whitney t-test, *p < 0.05, **p < 0.01, ns: not significant.

Targetome analysis of the dysregulated miRNAs in human primary keratinocytes or fibroblasts.
(a) A schematic diagram of targetome analysis. Microarray analysis was performed in keratinocytes (KC) or fibroblasts (FB) with miR-218-5p (b), miR-34a-5p (c), miR-34c-5p (d), miR-7704 (e), miR-96-5p (f), miR-424-5p (g, h), miR-450-5p (i), or miR-516-5p (j) overexpression (OE). Density plots of mRNA log2(fold change) are shown. Wilcoxon t-tests were performed to compare the TargetScan predicted strongest targets (purple) with the nontargets (blue) for each of these miRNAs. The conserved and experimentally validated targets are marked with red and green colors, respectively. Dotted lines depict the average log2(fold change) values for each mRNA group. (k) A heatmap shows the enrichment for venous ulcer (VU)-affected mRNAs and mRNA modules in the targets of miR-218-5p, miR-34a-5p, miR-34c-5p, or miR-7704 validated by the microarray analysis. Odds ratio (OR) and p values are shown when OR >1 and p value <0.05 (Fisher’s exact test).

Overexpression efficiency in primary keratinocytes or fibroblasts transfected with respective miRNAs. (a–g) qRT-PCR analysis of miRNAs in keratinocytes or fibroblasts transfected with miRNA mimics for 24 hr (n = 3–4).
The data are presented as mean ± standard error of the mean (SEM). **p < 0.01, ***p < 0.001, and ****p < 0.0001 by two-tailed Student’s t-test.

Analysis of the dysregulated miRNAs and their targets in vivo and in vitro.
Pearson’s correlation analysis of expression of the dysregulated miRNAs: (a) upper panels: miR-218-5p, miR-34a-5p, miR-34c-5p, and miR-7704; (b) upper panels: miR-96, miR-424, miR-450b, and miR-516b and their targets in RNA-seq data of human skin and wound samples. Gray circles indicate correlation coefficients >0. The gene expression of miRNA-regulated targets in epidermal CD45− cells from human skin and wound day 7 (a, b: middle panels). The target mRNAs are significantly changed (fold change ≤−1.2 and p value <0.05) by the indicated miRNAs in the microarray analysis (a, b: lower panels).

Cooperativity of venous ulcer (VU) pathology-relevant miRNAs.
For the microarray data in keratinocytes (KC) or fibroblasts (FB) with miR-218-5p, miR-96-5p, miR-34a-5p, miR-34c-5p, miR-7704, miR-424-5p, and miR-516-5p overexpression, we analyzed Kyoto Encyclopedia of Genes and Genome (KEGG) pathways, biological processes (a), and Molecular Signatures Database (MSigDB) hallmarks (b) enriched by the genes up- or downregulated by these miRNAs. (c) Venn diagram shows the number of genes in the E2F target pathway regulated by miR-34a/c-5p and/or miR-424-5p. The direct targets of each miRNA are depicted as a STRING Protein Network in (d). The number of inflammatory response-related genes regulated by miR-34a/c-5p and miR-516b-5p and ranked log2fold changes of overlapped genes from microarray analysis are shown in (e, f).
-
Figure 8—source data 1
Gene ontology analysis of the miRNA-regulated genes in the microarray data.
- https://cdn.elifesciences.org/articles/80322/elife-80322-fig8-data1-v2.xlsx

Cooperation of miR-34a, miR-424, and miR-516 in regulating keratinocyte proliferation, migration, and inflammatory response.
Ki-67 expression was detected in keratinocytes transfected with miR-34a-5p or miR-424-5p mimics alone or both mimics for 24 hr (n = 3) by qRT-PCR (a) and immunofluorescence staining (n = 4–6) (b). (c) The growth of the transfected keratinocytes (n = 3) was analyzed with a live cell imaging system. The migration of keratinocytes transfected with miR-34a-5p or miR-424-5p mimics alone or both mimics (d, e) (n = 8–10) or miRNA inhibitors (f, g) (n = 10) was analyzed with a live cell imaging system (scale bar, 300μm). (h) qRT-PCR analysis of CCL20 in keratinocytes transfected with miR-34a-5p or miR-516b-5p mimics alone or both mimics for 24 hr (n = 3). (i) Proposed mechanism by which venous ulcer (VU)-dysregulated miRNAs cooperatively contribute to the stalled wound healing characterized with failed inflammation–proliferation transition (schematic generated by BioRender). *p < 0.05; **p < 0.01; ***p < 0.001; and ****p < 0.0001 by one-way analysis of variance (ANOVA) with uncorrected Fisher’s least significant difference (LSD) (a, b, h) and two-way ANOVA (c, d, f). Data are presented as mean ± standard deviation (SD) (a–c, h) or mean ± standard error of the mean (SEM) (d, f).
-
Figure 9—source data 1
Cooperation of miR-34a, miR-424, and miR-516 in regulating keratinocyte proliferation and inflammatory response.
- https://cdn.elifesciences.org/articles/80322/elife-80322-fig9-data1-v2.xlsx

Cooperation of miR-34a and miR-424 in regulating keratinocyte proliferation.
(a) In the Kyoto Encyclopedia of Genes and Genome (KEGG) cell cycle pathway map, genes regulated by miR-34a-5p or miR-424-5p or both miRNAs as shown in the microarray analysis of keratinocytes are marked with green, blue, and orange colors, respectively. Among the regulated genes, the miRNA targets are highlighted in red and labeled with * for miR-34a targets, # for miR-424 targets, *# for the gene cotargeted by both miRNAs. (b) Ki-67 expression was detected in keratinocytes transfected with miR-34a-5p or miR-424-5p mimics alone or both mimics for 24 hr by immunofluorescence staining. Cell nuclei were costained with 4′,6-diamidino-2-phenylindole (DAPI). Ki-67 + cell nuclei are highlighted with white arrows in the representative photographs (scale bar: 50 μm). (c) The growth of the transfected keratinocytes (n = 3) was analyzed with a live cell imaging system, and representative photographs of cells at 0 and 96 hr are shown (scale bar, 300μm.).
Keratinocyte growth was analyzed with a live cell imaging system.
(1) Keratinocytes cotransfected with miR-34a-5p and control mimics. (2) Keratinocytes cotransfected with miR-424-5p and control mimics. (3) Keratinocytes cotransfected with miR-34a-5p and miR-424-5p mimics.
Tables
Characteristics of the healthy donors and the venous ulcer patients.
Characteristics | Healthy donors | Patients with VU |
---|---|---|
Study population (n) | 10 | 12 |
Age, years (mean ± SD) | 65.3 ± 3.2 | 75.8 ± 12.0 |
Ethnicity | Caucasian | Caucasian |
Gender (male:female) | 2:8 | 5:7 |
Biopsy location | Lower leg | Lower leg |
Wound duration | Acute (1 or 7 days after injury) | 3.7 ± 5.3 years |
-
SD, standard deviation; VU, venous ulcer (n = samples used for RNA-seq and qRT-PCR validation).
Characteristics of the patients with venous ulcer.
Patient | Sex | Age(years) | Ethnicity | Wound size (cm) | Wound duration (years) | Wound location | Experiment |
---|---|---|---|---|---|---|---|
1 | M | 86 | Caucasian | 6 × 5 | 4 | Lower leg | RNA-seq and qRT-PCR |
2 | F | 68 | Caucasian | 15 × 15 | 2 | Lower leg | RNA-seq and qRT-PCR |
3 | M | 70 | Caucasian | 3 × 0.5 | 0.3 | Lower leg | RNA-seq and qRT-PCR |
4 | F | 78 | Caucasian | 15 × 12 | 1.5 | Lower leg | RNA-seq and qRT-PCR |
5 | F | 87 | Caucasian | 3 × 2 + 8 × 4 | 2.5 | Lower leg | RNA-seq and qRT-PCR |
6 | F | 73 | Caucasian | 3 × 4 | 3.5 | Lower leg | qRT-PCR |
7 | M | 99 | Caucasian | 20 × 10 | 0.5 | Lower leg | qRT-PCR |
8 | F | 71 | Caucasian | 20 × 20 | 20 | Lower leg | qRT-PCR |
9 | F | 77 | Caucasian | 2.5 × 3 + 15 × 15 | 4.5 | Lower leg | qRT-PCR |
10 | M | 51 | Caucasian | 2 × 1.5 | 3 | Lower leg | qRT-PCR |
11 | M | 69 | Caucasian | 12 × 15 | 1 | Lower leg | qRT-PCR |
12 | F | 81 | Caucasian | 7 × 2.5 | 1 | Lower leg | qRT-PCR |
13 | M | 65 | Caucasian | 3 × 3 | 2.2 | Lower leg | LCM |
14 | M | 78 | Caucasian | 10 × 5 | 0.9 | Lower leg | LCM |
15 | M | 63 | Caucasian | 4 × 4 | 0.3 | Lower leg | LCM |
16 | M | 50 | Caucasian | 4.5 × 1 | 0.4 | Lower leg | LCM |
17 | M | 71 | Caucasian | 13 × 10 | 1 | Lower leg | LCM |
18 | M | 92 | Caucasian | 12 × 12 | 2 | Lower leg | LCM |
19 | F | 87 | Caucasian | 3 × 3 | 2.5 | Lower leg | LCM |
-
M, male; F, female.
Characteristics of the healthy donors.
Donor | Sex | Age(years) | Ethnicity | Wound location | Experiment |
---|---|---|---|---|---|
1 | F | 66 | Caucasian | Lower leg | RNA-seq |
2 | M | 69 | Caucasian | Lower leg | RNA-seq |
3 | F | 67 | Caucasian | Lower leg | RNA-seq |
4 | M | 69 | Caucasian | Lower leg | RNA-seq and qRT-PCR |
5 | F | 64 | Caucasian | Lower leg | RNA-seq and qRT-PCR |
6 | F | 60 | Caucasian | Lower leg | qRT-PCR |
7 | F | 66 | Caucasian | Lower leg | qRT-PCR |
8 | F | 60 | Caucasian | Lower leg | qRT-PCR |
9 | F | 67 | Caucasian | Lower leg | qRT-PCR |
10 | F | 65 | Caucasian | Lower leg | qRT-PCR |
11 | F | 66 | Caucasian | Lower back | Cell isolation; RNA-seq |
12 | M | 69 | Caucasian | Lower back | Cell isolation; RNA-seq |
13 | F | 67 | Caucasian | Lower back | Cell isolation; RNA-seq |
14 | M | 69 | Caucasian | Lower back | Cell isolation; RNA-seq |
15 | F | 64 | Caucasian | Lower back | Cell isolation; RNA-seq |
16 | F | 42 | Caucasian | Lower back | LCM |
17 | M | 58 | Caucasian | Lower back | LCM |
18 | M | 42 | Caucasian | Lower back | LCM |
19 | M | 28 | Caucasian | Lower back | LCM |
20 | M | 30 | Caucasian | Lower back | LCM |
21 | F | 46 | Caucasian | Lower back | LCM |
22 | F | 47 | Caucasian | Lower back | LCM |
-
M, male; F, female.
Summary of published microRNA profiling data of human skin wounds.
References | Study design |
---|---|
J Clin Invest. 2015 Aug 3;125(8):3008-26. a study from our group | Skin from healthy donors (n = 5) at 0 hours and 24 hours after injury. |
Burns. 2012 Jun;38(4):534-40. | Denatured dermis and paired normal skin of burn patients (n=3). |
J Cell Physiol. 2022 Feb;237(2):1429-1439 | Skin and day 14 wounds from both nonlesional and lesional sites from patients with vitiligo (n = 5). |
Exp Dermatol. 2021 Aug;30(8):1065-1072. doi: 10.1111/exd.14405. | The parental control vs. iPSC-derived fibroblasts from non-healing diabetic foot ulcers, the intact foot skin of diabetic patients, or the healthy foot skin of non-diabetic patients (n = 2 or 3). |
Additional files
-
Supplementary file 1
Quality control of small RNA-sequencing data.
- https://cdn.elifesciences.org/articles/80322/elife-80322-supp1-v2.docx
-
Supplementary file 2
Quality control of rRNA-depleted total RNA-sequencing data.
- https://cdn.elifesciences.org/articles/80322/elife-80322-supp2-v2.docx
-
Supplementary file 3
Summary of microRNAs reported to regulate skin wound healing.
- https://cdn.elifesciences.org/articles/80322/elife-80322-supp3-v2.xlsx
-
Transparent reporting form
- https://cdn.elifesciences.org/articles/80322/elife-80322-transrepform1-v2.docx