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

  1. Zhuang Liu
  2. Letian Zhang
  3. Maria A Toma
  4. Dongqing Li
  5. Xiaowei Bian
  6. Irena Pastar
  7. Marjana Tomic-Canic
  8. Pehr Sommar  Is a corresponding author
  9. Ning Xu Landén  Is a corresponding author
  1. Dermatology and Venereology Division, Department of Medicine Solna, Center for Molecular Medicine, Karolinska Institutet, Sweden
  2. Key Laboratory of Basic and Translational Research on Immune-Mediated Skin Diseases, Chinese Academy of Medical Sciences; Jiangsu Key Laboratory of Molecular Biology for Skin Diseases and STIs, Institute of Dermatology, Chinese Academy of Medical Sciences and Peking Union Medical College, China
  3. Wound Healing and Regenerative Medicine Research Program, Dr Phillip Frost Department of Dermatology and Cutaneous Surgery, University of Miami Miller School of Medicine, United States
  4. Department of Plastic and Reconstructive Surgery, Karolinska University Hospital, Sweden
  5. Ming Wai Lau Centre for Reparative Medicine, Stockholm Node, Karolinska Institute, Sweden
11 figures, 4 tables and 4 additional files

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
Figure 2 with 5 supplements
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
Figure 2—figure supplement 1
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.

Figure 2—figure supplement 2
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.

Figure 2—figure supplement 3
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.

Figure 2—figure supplement 4
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.

Figure 2—figure supplement 5
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.

Figure 4 with 2 supplements
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
Figure 4—figure supplement 1
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.

Figure 4—figure supplement 2
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.

Figure 5 with 1 supplement
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
Figure 5—figure supplement 1
RNA-sequencing results for the miRNAs selected for experimental validation.

Mann–Whitney t-test, *p < 0.05, **p < 0.01, ns: not significant.

Figure 6 with 1 supplement
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).

Figure 6—figure supplement 1
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
Figure 9 with 2 supplements
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
Figure 9—figure supplement 1
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.).

Figure 9—video 1
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.

Appendix 1—figure 1
Summary of microRNAs reported to regulate skin wound healing.
Author response image 1
Immunostaining of LGR4 in the skin, normal wound (NW), and VUs (n = 3).

Adapted from our recent publication PMID: 31376385.

Tables

Table 1
Characteristics of the healthy donors and the venous ulcer patients.
CharacteristicsHealthy donorsPatients with VU
Study population (n)1012
Age, years (mean ± SD)65.3 ± 3.275.8 ± 12.0
EthnicityCaucasianCaucasian
Gender (male:female)2:85:7
Biopsy locationLower legLower leg
Wound durationAcute (1 or 7 days after injury)3.7 ± 5.3 years
  1. SD, standard deviation; VU, venous ulcer (n = samples used for RNA-seq and qRT-PCR validation).

Table 2
Characteristics of the patients with venous ulcer.
PatientSexAge(years)EthnicityWound size (cm)Wound duration (years)Wound locationExperiment
1M86Caucasian6 × 54Lower legRNA-seq and qRT-PCR
2F68Caucasian15 × 152Lower legRNA-seq and qRT-PCR
3M70Caucasian3 × 0.50.3Lower legRNA-seq and qRT-PCR
4F78Caucasian15 × 121.5Lower legRNA-seq and qRT-PCR
5F87Caucasian3 × 2 + 8 × 42.5Lower legRNA-seq and qRT-PCR
6F73Caucasian3 × 43.5Lower legqRT-PCR
7M99Caucasian20 × 100.5Lower legqRT-PCR
8F71Caucasian20 × 2020Lower legqRT-PCR
9F77Caucasian2.5 × 3 + 15 × 154.5Lower legqRT-PCR
10M51Caucasian2 × 1.53Lower legqRT-PCR
11M69Caucasian12 × 151Lower legqRT-PCR
12F81Caucasian7 × 2.51Lower legqRT-PCR
13M65Caucasian3 × 32.2Lower legLCM
14M78Caucasian10 × 50.9Lower legLCM
15M63Caucasian4 × 40.3Lower legLCM
16M50Caucasian4.5 × 10.4Lower legLCM
17M71Caucasian13 × 101Lower legLCM
18M92Caucasian12 × 122Lower legLCM
19F87Caucasian3 × 32.5Lower legLCM
  1. M, male; F, female.

Table 3
Characteristics of the healthy donors.
DonorSexAge(years)EthnicityWound locationExperiment
1F66CaucasianLower legRNA-seq
2M69CaucasianLower legRNA-seq
3F67CaucasianLower legRNA-seq
4M69CaucasianLower legRNA-seq and qRT-PCR
5F64CaucasianLower legRNA-seq and qRT-PCR
6F60CaucasianLower legqRT-PCR
7F66CaucasianLower legqRT-PCR
8F60CaucasianLower legqRT-PCR
9F67CaucasianLower legqRT-PCR
10F65CaucasianLower legqRT-PCR
11F66CaucasianLower backCell isolation; RNA-seq
12M69CaucasianLower backCell isolation; RNA-seq
13F67CaucasianLower backCell isolation; RNA-seq
14M69CaucasianLower backCell isolation; RNA-seq
15F64CaucasianLower backCell isolation; RNA-seq
16F42CaucasianLower backLCM
17M58CaucasianLower backLCM
18M42CaucasianLower backLCM
19M28CaucasianLower backLCM
20M30CaucasianLower backLCM
21F46CaucasianLower backLCM
22F47CaucasianLower backLCM
  1. M, male; F, female.

Author response table 1
Summary of published microRNA profiling data of human skin wounds.
ReferencesStudy 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-1439Skin 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

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Zhuang Liu
  2. Letian Zhang
  3. Maria A Toma
  4. Dongqing Li
  5. Xiaowei Bian
  6. Irena Pastar
  7. Marjana Tomic-Canic
  8. Pehr Sommar
  9. Ning Xu Landén
(2022)
Integrative small and long RNA omics analysis of human healing and nonhealing wounds discovers cooperating microRNAs as therapeutic targets
eLife 11:e80322.
https://doi.org/10.7554/eLife.80322