A Next-Generation Sequencing of Plasma Exosome-Derived microRNAs and Target Gene Analysis with a Microarray Database of Thermally Injured Skins: Identification of Blood-to-Tissue Interactions at Early Burn Stage.

BACKGROUND: Plasma exosome-derived microRNA (miRNA) profiles following thermal injury and their relationship with gene expression derangements in burned skin remain unexplored. This study focused on the identification of key miRNA-mRNA axes in potential blood-to-tissue interactions at early burn stage.
METHODS: Plasma exosomes were obtained from 6 severe burn patients 4-7 days post injury and 6 healthy volunteers. Next-generation sequencing (NGS) of exosomal small RNAs presented the differentially expressed miRNAs (DEMs). Target genes of the DEMs were predicted in the mirDIP database. Dataset GSE8056 was enrolled to acquire differentially expressed genes (DEGs) in burned skin compared to normal skin. Overlap between the DEGs and target genes of the DEMs were focus genes. The protein-protein interaction (PPI) network and enrichment analyses of the focus genes demonstrated hub genes and suggested underlying mechanisms and pathways. The hub genes and upstream DEMs were selected to construct key miRNA-mRNA axes.
RESULTS: The NGS of plasma exosome-derived small RNAs identified 85 DEMs (14 downregulated miRNAs and 71 upregulated miRNAs) with 12,901 predicted target genes. Dataset GSE8056 exhibited 1861 DEGs in partial-thickness burned skins 4-7 days postburn. The overlap between DEGs and target genes of DEMs displayed 1058 focus genes. The top 9 hub genes (CDK1, CCNB1, CCNA2, BUB1B, PLK1, KIF11, AURKA, NUSAP1 and CDCA8) in the PPI network of the focus genes pointed to 16 upstream miRNAs in DEMs, including 4 downregulated miRNAs (hsa-miR-6848-3p, has-miR-4684-3p, has-miR-4786-5p and has-miR-365a-5p) and 12 upregulated miRNAs (hsa-miR-6751-3p, hsa-miR-718, hsa-miR-4754, hsa-miR-6754-3p, hsa-miR-4739, hsa-miR-6739-5p, hsa-miR-6884-3p, hsa-miR-1224-3p, hsa-miR-6878-3p, hsa-miR-6795-3p, hsa-miR-550a-3p, and hsa-miR-550b-3p). A key miRNA-mRNA network in potential blood-to-tissue interactions at early burn stage was therefore constructed.
CONCLUSION: An NGS and bioinformatic analysis in the study identified key miRNA-mRNA axes in potential blood-to-tissue interactions at early burn stage, suggesting plasma exosome-derived miRNAs may impact on the alteration patterns of gene expressions in a burn wound.

Introduction

Burn injury remains a major global health problem exerting heavy burdens on victims.1 Therapeutic intervention to primary tissue damage and following wound conversion is a constant topic of total burn care.2 Rapid transvascular fluid filtration features prominently in thermal injuries, causing generalized edema in remote tissues as well as directly burned areas.3 Kremer et al first reported endothelial activation and albumin leakage in unburned recipients induced by transfusion of plasma from donor rats with 30% total body surface area (TBSA) full-thickness burns.4 Then researchers illustrated that systemic edema after infusion of burn plasma donated by 30% TBSA full-thickness scald rats was attenuated with the application of cerium nitrate5 and methysergide.6 Accumulated studies suggested undefined circulating factors postburn, which can be inhibited by several experimental drugs, contributed to undesirable changes in tissues.5–7 These circulating factors, hampering wound recovery, haunted researchers for decades and were supposed to arise from general conditions of burn patients such as hypermetabolism, inflammation and immunosuppression.8 Complicated pathophysiological courses following burns stimulated endeavors to portray underlying molecular derangements. However, convincing identifications or intervention therapies of the circulating factors postburn were still absent.

We hypothesized that the exosome was a veiled circulating factor strongly associated with cutaneous changes postburn. Within 30–150 nm in diameter, the exosome is a protagonist in extracellular vesicles research.9 Exosomes exist stably in almost all body fluid, serving intercellular communication mainly by lipid bilayers-protected nucleotides and proteins.10 Dramatic responses to severe burn may turn the peripheral blood into a “crowded playground” of exosomes secreted by various activated cells. The identification of circulatory exosomes and their contents exhibited potential value in diagnosis and prognosis for many diseases.9–11 MicroRNAs (miRNAs) are small non-coding RNA molecules the length of about 22 nucleotides, and loaded as regular shipments of exosomes participating in post-transcriptional regulations of gene expression by silencing mRNAs.12 Soared vascular permeability postburn and consequent extravasation of plasma would advance the opportunities for exosomes to flood into local tissues, exerting long-lasting adverse effects and/or promoting reparative processes in wounds through miRNA-mRNA axes. To explore the potential blood-to-tissue interactions at early burn stage, we recruited clinical samples to conduct the enrichment and characterization of plasma exosomes. More importantly, portraits of plasma exosome-derived miRNAs at an early stage following severe burn were first described in our study. Next-generation sequencing (NGS) analysis presented the differentially expressed miRNAs (DEMs). Focus genes were selected by overlapping the predicted target genes of DEMs and the differentially expressed genes (DEGs) in thermally injured skins detected from Gene Expression Omnibus (GEO) database. Top 10 Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways with the greatest significance in enrichment analyses of genes helped to elaborate their functions and underlying pathways. Protein–protein interaction (PPI) network of the focus genes delivered candidates from hub genes and upstream DEMs to construct a miRNA-mRNA crosstalk network.

Materials and Methods

Demographics and Samples Collection

Peripheral blood samples from 6 severe burn patients and 6 healthy volunteers were obtained with the written informed consent of participants. At least 6 mL whole blood of each individual was collected by using EDTA vacutainer tubes. The study was conducted according to the guidelines of the Declaration of Helsinki, and all procedures were approved by the Human Subjects Review Board of the First Affiliated Hospital of Anhui Medical University (number: PJ2019-06-07). Severity of burns was evaluated through the Abbreviated Burn Severity Index (ABSI).13 The inclusion criteria of severe burn patients included: thermal injury; % TBSA burn ≥ 60; ABSI score ≥ 10; 18 ≤ age ≤ 80; 4 ≤ days postburn ≤ 7. And the exclusion criteria included patients under uncorrected severe shock; sepsis; pregnancy; cancer; hematopathy. The demographics and clinical information of sampled patients including sex, age, % TBSA burn, % deep burn (deep second degree and higher-grade burns), ABSI score and sampling day postburn were described in Table 1. Within 10 min after collection, all samples were centrifuged at 3500 g for 15 min. Plasma aliquots were stored at −80 °C.

Table 1

Demographics for Participants Sampled for NGS and Patients from GSE8056

Research and Sample Type	Control Group	Sex	Age (Years)	Mean Age	Burn Group	Sex	Age (Years)	Mean Age	% TBSA Burn (% Deep Burn)	Mean%TBSA	ABSI Score	Day Postburn
NGS of plasma exosome-derived small RNA	Normal 01	F	24	28.3 ±5.28	Burn 01	F	54	39 ±16.69	60 (32)	79.83 ±19.99	11	5
	Normal 02	F	21		Burn 02	M	28		60 (46)		10	5
	Normal 03	F	29		Burn 03	F	25		65 (20)		11	5
	Normal 04	M	28		Burn 04	M	54		97 (80)		15	4
	Normal 05	M	33		Burn 05	F	54		98 (92)		16	4
	Normal 06	M	35		Burn 06	M	19		99 (95)		13	4
cDNA-microarray analysis of skin	GSM198875-01	F	35	38.4 ±14.37	GSM198869-01	M	8	30.8 ±14.83	6	14.67 ±12.94	No record	5
	GSM198875-02	F	26		GSM198869-02	M	19		8			4
	GSM198875-03	F	43		GSM198869-03	M	30		10			4
	GSM198875-04	F	20		GSM198869-04	F	11		5			5
	GSM198875-05	F	20		GSM198869-05	F	31		15			5
	GSM198876-01	F	32		GSM198870-01	M	26		30			4
	GSM198876-02	F	43		GSM198870-02	M	33		5			5
	GSM198876-03	F	52		GSM198870-03	F	37		4			6
	GSM198876-04	F	56		GSM198870-04	F	38		20			7
	GSM198876-05	F	44		GSM198870-05	M	18		20			7
	GSM198877-01	F	55		GSM198871-01	M	51		50			4
	GSM198877-02	M	43		GSM198871-02	M	25		10			5
	GSM198877-03	F	22		GSM198871-03	M	42		27			5
	GSM198877-04	M	22		GSM198871-04	M	28		5			7
	GSM198877-05	F	63		GSM198871-05	M	65		5			6

Abbreviations: TBSA, total body surface area; ABSI, abbreviate burn severity index; NGS, next-generation sequencing; F, female; M, male.

Exosome Isolation and Characterization

Enrichment and purification of plasma exosomes was conducted by using exoRNeasy Midi Kit (77144, Qiagen, Germany). Exosomes were isolated according to the protocols described in the user manual and then processed for verification. Tecnai G2 Spirit BioTwin Transmission electron microscopy (TEM) (FEI, USA) visualized the morphology of exosomes. Nanoparticle tracking analysis (NTA) was performed to track and measure the particle size and distribution using ZetaVIEW S/N 17–310 (Particle Metrix, Germany). Western blotting (WB) analyzed the exosome protein marker CD63 and endoplasmic reticulum protein marker Calnexin. Primary antibodies employed included rabbit anti-human CD63 antibody (EXOAB-CD63A-1, SBI, USA), and rabbit anti-human Calnexin antibody (YT0613, ImmunoWay, USA). Goat anti-rabbit HRP antibody (SBI, USA) was used as a secondary antibody.

Small RNA Library Construction and Next-Generation Sequencing

Extraction and purification of total RNA from exosomes was also conducted with the help of exoRNeasy Midi Kit (77144, Qiagen, Germany). Small RNA libraries were then prepared with the application of the QIAseq miRNA Library Kit (331502, Qiagen, Germany). Reverse transcription (RT) primers with unique molecular indices (UMIs)14 were introduced to analyze the quantification of miRNA expressions during cDNA synthesis and PCR amplification. The Next-generation sequencing (NGS) was performed on the NovaSeq 6000 System (Illumina, USA).

Data Analysis and DEM Identification

Sequence reads were trimmed, filtered and then screened with the aid of UMIs. Subsequent RNA-seq reads were mapped to miRbase, piRNAbank and Rfam databases. Expression matrix of quantified UMI counts of miRNAs was normalized to counts per million (CPM) and calculated to relative log expression via the EdgeR package.15 Differentially expressed miRNAs (DEMs) were identified according to the threshold set at FDR < 0.05 and |log2FC (fold change)|≥1.

Target Genes Prediction and Enrichment Analyses

Target genes of DEMs were predicted in mirDIP (), which integrated information of miRNAs and genes in Homo sapiens from 30 major databases including miRbase, miRDB and TargetScan.16 The filter option of minimum integrated score, which combined confidence scores from all available predictions, was set “High (Top5%)” to retrieve reliable target genes of DEMs. Enrichment analyses of GO and KEGG were conducted using the ClusterProfiler17 package in R software (version 4.1.2). All significantly enriched GO terms and KEGG pathways in this study were identified based on the criterion set at adj.P.value < 0.05. Top 10 GO terms and KEGG pathways with the lowest adj.P.value helped to elaborate enriched functions and underlying pathways.

Gene Expression Microarray Data Matching

In order to narrow the spectrum of DEMs’ target genes to study the possible blood-to-tissue interactions at the early burn stage, we searched gene expression profiles of burned skin in the GEO repository (). Only one dataset under the accession number GSE8056 was well matched. Microarray data documented in the GSE8056 dataset consisted of gene expression profiles of normal skins from 15 patients as the control group and thermally injured skins from 45 burn patients.18 The burn patients were clustered into 3 groups according to the post-burn days (PBD), including PBD 0–3 group (15 patients), PBD 4–7 group (15 patients) and PBD 7–17 group (15 patients), to reveal the temporal gene expression profiles in thermally injured skin. All the burn patients admitted into the microarray study suffered deep partial-thickness to full-thickness burns warranting operative excision. Burn skin specimens were sampled from a partial-thickness burned area since a full-thickness burned skin undergoing necrosis obviated the detection of gene expressions. Additionally, seriously infected specimens were excluded. More detailed clinical information such as ABSI score of the enrolled patients, however, was undescribed.

Gene expression matrices of burn wounds in the PBD 4–7 group and the normal human skin performed on GPL570 (Affymetrix Human Genome U133 Plus 2.0 Array) were downloaded and analyzed. Specifically, each sub-dataset in accession number of GSM198869, GSM198870 and GSM198871 pooled RNA samples from 5 burn patients, while GSM198875, GSM198876 and GSM198877 pooled RNA samples from 5 normal skin specimens respectively. Details of sampled patients are also described in Table 1. Differentially expressed genes (DEGs) were identified using the LIMMA19 package with enrichment analyses conducted using the ClusterProfiler package in R software.

Venn and Enrichment Analysis of the Focus Genes

Venn analysis presented shared genes as an overlap between the target genes of DEMs and the DEGs. These shared genes were named focus genes suggesting focused targets of DEMs involved in blood-tissue interactions at the early burn stage. Enrichment analyses of GO and KEGG were repeated on the focus genes.

Construction of PPI and miRNA-mRNA Networks

A protein–protein interaction (PPI) network of the focus genes was constructed in Cytoscape software (version 3.7.2) based on the STRING database ().20 Hub genes in the PPI network were selected by degree using the cytoHubba21 plugin of Cytoscape. The hub genes were finally selected as key genes to construct a network of key miRNA-mRNA axes to determine the effect of plasma exosomal miRNAs postburn on partial-thickness burned skins.

Results

Characterization of Isolated Exosomes

Extracellular vesicles were observed in the transmission electron microscopic image, and the size of the vesicles met the specified diameter of exosomes (Figure 1A). NTA confirmed the size distributions measured up to exosomes (Figure 1B). Positive expression of exosomal marker protein CD63, along with the negative expression of endoplasmic reticulum marker protein calnexin, assured the credible presence of isolated exosomes (Figure 1C).

Figure 1

Characterization of plasma exosomes from burn patients and control group. (A) Visualized exosomal morphology in transmission electron microscopy (TEM). Scale bar = 200 nm. (B) Size distribution of plasma exosomes by nanoparticle tracking analysis (NTA). (C) Western blotting (WB) analyses of the exosome protein marker CD63 and the endoplasmic reticulum protein marker Calnexin.

NGS of Exosomal Small RNA and Identification of DEMs

The NGS presented the plasma exosome-derived small RNA profiles. Over 16 million clean reads were yielded from all samples. Mapped miRbase mature reads counted as 16.05% of all reads, and the ratios of other classifications were presented in Table 2.

Table 2

Genome Alignment Ratios of RNA Read Counts Concerning Plasma Exosome-Derived Small RNAs

Reads	Counts	Ratio
All clean reads	16,534,063	100%
Noadapter reads	400,031	2.42%
Short UMIs reads	8252	0.05%
Mapped miRBase mature reads	2,653,251	16.05%
Mapped miRBase hairpin reads	295,149	1.79%
Mapped piRNA reads	107,862	0.65%
Mapped Rfam reads	998,098	6.04%

According to the threshold set at FDR<0.05 and | log2FC (fold change) | ≥ 1, 85 DEMs including 14 downregulated miRNAs and 71 upregulated miRNAs were rendered (Table 3). All DEMs were displayed in a volcano plot and heatmap respectively (Figure 2A and B).

Table 3

Differentially Expressed microRNAs (DEMs) in Plasma Exosomes at the Early Burn Stage

DEMs	Log2FC	P-value	FDR	Regulation	Targets Number in Focus Genes
hsa-miR-6873-3p	3.6939	2.79E-24	7.41E-21	Up	79
hsa-miR-6751-3p	2.8164	2.98E-14	3.95E-11	Up	50
hsa-miR-1281	2.6124	2.33E-11	1.65E-08	Up	51
hsa-miR-6741-3p	2.6198	2.49E-11	1.65E-08	Up	51
hsa-miR-7111-3p	2.4319	4.00E-10	2.08E-07	Up	68
hsa-miR-6776-3p	3.0672	4.69E-10	2.08E-07	Up	61
hsa-miR-6809-3p	2.4835	6.62E-09	2.51E-06	Up	79
hsa-miR-6830-3p	2.1379	8.48E-09	2.82E-06	Up	87
hsa-miR-6754-3p	2.5296	1.03E-08	3.03E-06	Up	62
hsa-miR-6880-3p	2.2448	1.55E-08	4.12E-06	Up	51
hsa-miR-718	2.9613	2.05E-08	4.94E-06	Up	50
hsa-miR-6739-5p	2.5980	2.27E-08	5.01E-06	Up	69
hsa-miR-4646-3p	2.2586	3.31E-08	6.76E-06	Up	65
hsa-miR-6894-3p	2.1834	8.96E-08	1.70E-05	Up	68
hsa-miR-6739-3p	2.6488	1.53E-07	2.71E-05	Up	74
hsa-miR-10398-5p	2.1290	1.67E-07	2.77E-05	Up	27
hsa-miR-4667-5p	2.1068	2.27E-07	3.55E-05	Up	67
hsa-miR-1224-3p	1.9389	5.33E-07	7.87E-05	Up	45
hsa-miR-4434	2.7177	6.80E-07	9.18E-05	Up	72
hsa-miR-150-5p	−1.7822	7.03E-07	9.18E-05	Down	87
hsa-miR-6823-3p	1.8625	7.26E-07	9.18E-05	Up	48
hsa-miR-6759-3p	2.1785	9.26E-07	1.12E-04	Up	57
hsa-miR-12116	1.6649	1.96E-06	2.26E-04	Up	42
hsa-miR-6884-3p	1.7718	2.56E-06	2.83E-04	Up	47
hsa-miR-4716-5p	1.8601	5.63E-06	5.98E-04	Up	61
hsa-miR-6883-3p	1.7636	6.57E-06	6.71E-04	Up	59
hsa-miR-7847-3p	1.4644	7.91E-06	7.78E-04	Up	61
hsa-miR-4516	1.7274	9.01E-06	8.54E-04	Up	72
hsa-miR-6797-3p	1.9338	1.18E-05	1.08E-03	Up	51
hsa-miR-335-3p	2.1197	1.38E-05	1.23E-03	Up	67
hsa-miR-6510-5p	1.7775	1.61E-05	1.38E-03	Up	57
hsa-miR-7107-3p	1.8947	1.76E-05	1.46E-03	Up	66
hsa-miR-550b-3p	2.8765	2.03E-05	1.63E-03	Up	72
hsa-miR-4667-3p	1.5593	3.42E-05	2.67E-03	Up	56
hsa-miR-657	1.5340	4.03E-05	3.05E-03	Up	46
hsa-miR-6886-3p	1.6442	4.21E-05	3.10E-03	Up	51
hsa-miR-6833-3p	1.5770	5.57E-05	4.00E-03	Up	71
hsa-miR-4646-5p	1.7058	5.73E-05	4.01E-03	Up	56
hsa-miR-10401-3p	2.2197	7.39E-05	5.03E-03	Up	26
hsa-miR-6756-3p	1.5415	7.73E-05	5.13E-03	Up	56
hsa-miR-365b-5p	−5.0048	9.60E-05	6.20E-03	Down	60
hsa-miR-6760-3p	1.4502	9.80E-05	6.20E-03	Up	54
hsa-miR-6867-3p	1.5707	1.04E-04	6.43E-03	Up	68
hsa-miR-550a-3p	2.1836	1.21E-04	7.32E-03	Up	68
hsa-miR-10400-5p	3.7818	1.45E-04	8.36E-03	Up	40
hsa-miR-4522	2.0049	1.45E-04	8.36E-03	Up	66
hsa-miR-4786-5p	−2.4633	2.01E-04	1.11E-02	Down	56
hsa-miR-5006-5p	1.7414	2.01E-04	1.11E-02	Up	61
hsa-miR-4684-3p	−3.7701	2.31E-04	1.25E-02	Down	75
hsa-miR-1539	−3.4463	2.37E-04	1.26E-02	Down	60
hsa-miR-6891-3p	1.5370	2.46E-04	1.27E-02	Up	57
hsa-miR-375-3p	−2.9319	2.49E-04	1.27E-02	Down	63
hsa-miR-4649-3p	2.5969	2.95E-04	1.45E-02	Up	62
hsa-miR-125a-5p	−1.4474	2.96E-04	1.45E-02	Down	64
hsa-miR-8052	−4.4882	3.01E-04	1.45E-02	Down	58
hsa-miR-7114-3p	1.7273	3.25E-04	1.54E-02	Up	44
hsa-miR-6795-3p	1.5891	3.47E-04	1.62E-02	Up	44
hsa-miR-7110-3p	1.3498	3.90E-04	1.78E-02	Up	76
hsa-miR-365a-5p	−3.0986	3.95E-04	1.78E-02	Down	55
hsa-miR-4685-3p	1.6637	4.35E-04	1.92E-02	Up	58
hsa-miR-4728-5p	1.7391	4.41E-04	1.92E-02	Up	63
hsa-miR-4739	1.2684	5.12E-04	2.18E-02	Up	65
hsa-miR-6747-3p	1.3237	5.16E-04	2.18E-02	Up	74
hsa-miR-3059-5p	1.4961	5.28E-04	2.19E-02	Up	62
hsa-miR-5193	1.1997	5.68E-04	2.32E-02	Up	65
hsa-miR-6882-3p	1.4466	6.14E-04	2.47E-02	Up	59
hsa-miR-10398-3p	3.1236	6.31E-04	2.50E-02	Up	37
hsa-miR-122b-3p	2.0180	7.49E-04	2.92E-02	Up	63
hsa-miR-1184	2.5324	7.84E-04	3.02E-02	Up	49
hsa-miR-6848-3p	−1.9548	8.17E-04	3.06E-02	Down	66
hsa-miR-6513-5p	1.4819	8.19E-04	3.06E-02	Up	81
hsa-miR-6826-3p	1.2170	8.49E-04	3.13E-02	Up	56
hsa-miR-6726-5p	−3.1241	9.51E-04	3.41E-02	Down	50
hsa-miR-4488	1.5293	9.58E-04	3.41E-02	Up	50
hsa-miR-4754	2.9626	9.63E-04	3.41E-02	Up	54
hsa-miR-6817-3p	3.3331	9.81E-04	3.43E-02	Up	76
hsa-miR-4644	1.1645	1.03E-03	3.56E-02	Up	59
hsa-miR-520g-5p	−1.9903	1.12E-03	3.82E-02	Down	60
hsa-miR-663b	1.6953	1.17E-03	3.94E-02	Up	48
hsa-miR-6878-3p	1.3082	1.25E-03	4.17E-02	Up	74
hsa-miR-6727-3p	1.3258	1.35E-03	4.43E-02	Up	60
hsa-miR-6816-5p	−3.4400	1.45E-03	4.66E-02	Down	52
hsa-miR-4783-5p	1.4944	1.45E-03	4.66E-02	Up	37
hsa-miR-4478	1.3517	1.56E-03	4.93E-02	Up	73
hsa-miR-631	−2.6145	1.59E-03	4.97E-02	Down	45

Figure 2

DEMs in plasma exosomes and DEGs in skin. (A) 85 DEMs identified from NGS of plasma exosomal small RNAs, including 14 downregulated miRNAs in green and 71 upregulated miRNAs in red, are exhibited in the volcano plot. (B) Heatmap of all DEMs with the intensity plot showing a relatively low expression in green and a relatively high expression in red. (C) 1861 DEGs identified from microarray data of skin gene expressions (GSE8056), including 1030 downregulated genes in green and 831 upregulated genes in red, are exhibited in the volcano plot. (D) Heatmap of top 100 DEGs with greatest significance, with the intensity plot showing a relatively low expression in green and a relatively high expression in red.

Target Genes Prediction and Enrichment Analysis

In the mirDIP database, results of all DEMs targeted genes were rated combined confidence scores from all available predictions in 30 reference databases. After the filter option of minimum score was set “High (Top5%)”, 12,901 target genes were produced.

The ClusterProfiler package in R software enriched 209 GO terms and 165 KEGG pathways of DEMs’ target genes. Top 10 GO terms and KEGG pathways are shown in Figure 3A and B. GO analysis indicated that the target genes enriched mainly in multiple GTPase activities, “protein serine/threonine kinase activity”, and activities regarding DNA-binding transcription. KEGG analysis demonstrated several classic function pathways such as “MAPK signaling pathway” and “Hippo signaling pathway”.

Figure 3

GO terms and KEGG pathways enriched in DEMs’ target genes, DEGs and focus genes. (A) Top 10 GO terms enriched in target genes of plasma exosome-derived DEMs. (B) Top 10 KEGG pathways of DEMs’ target genes. (C) Top 10 GO terms enriched in DEGs of skins (from GSE8056). (D) Top 10 KEGG pathways of DEGs. (E) Top 10 GO terms enriched in the focus genes. (F) All 7 KEGG pathways of the focus genes.

Identification of DEGs in GSE8056 and Enrichment Analysis

Compared to normal skins, microarray gene expressions of partial-thickness burned skins exhibited 1861 DEGs meeting the criteria of FDR < 0.05 and |log2FC (fold change)|≥1. The volcano plot presented total DEGs, while a heatmap showed top 100 DEGs with the greatest significance among 831 upregulated genes and 1030 downregulated genes (Figure 2C and D). The ClusterProfiler package in R software found 48 GO terms and 19 KEGG pathways of DEGs. The top 10 GO terms and KEGG pathways are shown in Figure 3C and D.

Venn and Enrichment Analyses of Focus Genes

By overlapping target genes of DEMs predicted in the mirDIP and the DEGs in GSE8056, 1058 shared genes were identified as focus genes (Figure 4). Intriguingly, the focus genes made up about 56.9% of the DEGs and all DEMs were predicted to have targets in the focus genes.

Figure 4

Venn analysis of DEMs’ target genes and DEGs. Up/down DEGs, differentially expressed genes significantly upregulated/downregulated in thermally injured skins from GSE8056. Up/down DEM targets, target genes of differentially expressed miRNAs significantly upregulated/downregulated in plasma exosomes from burn patients.

GO and KEGG analyses of focus genes suggested possible functions and mechanism of mRNA alterations related to blood-to-tissue interactions postburn mediated by exosomal miRNAs. For GO enrichment, top 10 terms were “immune receptor activity”, “extracellular matrix structural constituent”, “carbohydrate binding”, “cytokine activity”, “glycosaminoglycan binding”, “sulfur compound binding”, “cytokine receptor activity”, “cytokine binding”, “heparin binding” and “fibronectin binding”. KEGG pathway analysis indicated 7 pathways: “cell cycle”, “cell adhesion molecules”, “Glutathione metabolism”, “ECM-receptor interaction”, “hematopoietic cell lineage” and “JAK-STAT signaling pathway” (Figure 3E and F).

Based on STRING interactions with the highest confidence (0.900) and FDR < 0.05, the main PPI network of focus genes was plotted in Cytoscape (Figure 5A), containing 274 nodes and 755 edges. Nodes’ scores calculated and ranked by degree recommended top 9 hub genes using the cytoHubba. Summaries based on the GeneCards database ()22 and detailed information of the hub genes are listed in Table 4. These 9 hub genes were predicted to be target genes of 16 miRNAs in DEMs. A miRNA-mRNA network was accordingly constructed in Figure 5B.

Figure 5

Networks of potential blood-to-tissue interactions. (A) The protein–protein interaction network of focus genes. Size of the node is scaled by | log2FC (fold change) | of gene expressions; color of the node and its border represents the expression regulation of mRNA and predicted upstream miRNA, respectively. (B) An integrated miRNA-mRNA and PPI network of the key miRNAs and their target genes (including 9 hub genes) in the focus genes. Edge line between genes (proteins) is solid; edge line between miRNA (round rectangle) and mRNA (round) is in dots. Red: upregulated. Green: downregulated.

Table 4

Top 9 Hub Genes in the PPI Network of Focus Genes

Gene Symbol	Gene Name	GeneCards Summary of Protein Encoded	Degree Score	Log2FC in GSE8056	Regulation	P-value	FDR	Downregulated Upstream DEMs	Upregulated Upstream DEMs
CDK1	Cyclin dependent kinase 1	Ser/Thr protein kinase essential for G1/S and G2/M phase transitions	51	1.481599	Up	6.31E-04	1.99E-03		hsa-miR-4754
CCNB1	Cyclin B1	Cyclin member necessary to the G2/M transition	41	2.535214	Up	6.42E-06	6.90E-05		hsa-miR-6884-3p, hsa-miR-718
CCNA2	Cyclin A2	Cyclin member promoting G1/S and G2/M transition	39	2.514916	Up	1.45E-05	1.22E-04		hsa-miR-550b-3p, hsa-miR-550a-3p, hsa-miR-6795-3p
BUB1B	BUB1 Mitotic Checkpoint Serine/Threonine Kinase B	Checkpoint kinase delaying anaphase and ensuring proper chromosome segregation	35	2.196643	Up	3.54E-06	4.65E-05		hsa-miR-6754-3p
PLK1	Polo like kinase 1	Serine/threonine-protein kinase prominent throughout M phase	31	2.341512	Up	7.04E-05	3.84E-04		hsa-miR-1224-3p, hsa-miR-4754
KIF11	Kinesin family member 11	Motor protein required for establishing a bipolar spindle	29	2.446886	Up	9.00E-04	2.61E-03	hsa-miR-6848-3p	hsa-miR-4739
AURKA	Urora kinase A	Mitotic serine/threonine kinase regulating microtubules at the spindle pole during chromosome segregation	27	1.387036	Up	9.31E-06	8.90E-05	hsa-miR-4684-3p	hsa-miR-4739, hsa-miR-6739-5p, hsa-miR-6751-3p
NUSAP1	Nucleolar and spindle-associated protein 1	Regulator to bundle and stabilize microtubules	27	2.091017	Up	1.41E-04	6.41E-04	hsa-miR-4786-5p, hsa-miR-365a-5p	hsa-miR-4754
CDCA8	Cell division cycle associated 8	Component of the chromosomal passenger complex essential for mitosis and cell division	27	1.059758	Up	7.16E-04	2.19E-03		hsa-miR-4739, hsa-miR-6878-3p

Abbreviations: PPI, protein–protein interaction; DEMs, differentially expressed miRNAs.

Discussion

Investigations into exosomes and exosome-derived miRNAs from the peripheral blood of trauma patients are scarce compared to cancers. The enrichment and identification of circulating exosomes are hitherto absent to our knowledge in the field of burns. Although it is difficult to trace the exact source of circulating exosomes among multiple burn-activated cells, exosomal contents well-protected in stable vesicles may reflect the general patterns of fierce systematic turbulence postburn, especially in patients with greater TBSA. Enrichment of plasma exosomes and exosomal small RNA-seq at early burn stage were conducted for the first time in our investigation. The sampling time was set at 4–7 days postburn for clinical considerations. Treatment of burn shock necessitates massive fluid resuscitation, and fresh frozen plasma (FFP) is a preferred option of colloid administration which ameliorates endothelial dysfunction as well as restores circulation volume.23,24 Plasma samples of burn patients in the shock stage would inevitably bring in excess FFP transfusion-induced exogenous exosomes. Moreover, within this time frame, the hemorrhagic shock of burn patients had been generally addressed with low risk of severe infection caused by various formidable pathogens,25 therefore the results may accurately refer to the burn injury instead of shock or sepsis. A total of 6 severe burn patients with ABSI score ≥10 and an average burned area of about 80% TBSA hospitalized in our burns department were recruited into the study. Firstly, TEM, NTA and WB proved the reliable existence of exosomes in the peripheral blood of burn patients sampled 4–7 days postburn as well as healthy volunteers. Then next-generation sequencing drew the profile of plasma exosomal small RNAs, showing that small RNAs were abundant in plasma exosomes. Compared to the control group, 85 miRNAs were significantly differentially expressed in the enrolled 6 severe burn patients. Enrichment analyses of predicted DEMs’ target genes could reveal the courses of postburn systemic response in sampled patients.

It is extensively illustrated that circulating exosomes secreted by tumor cells promoted the formation of the local microenvironment facilitating metastatic lesions.26,27 However, the effect of circulating factors, especially exosomal cargoes, on the healing processes of burn wounds was unexplored. Locally, massive and rapid activation of fibroblast proliferation is detected to repair tissue defects shortly after thermal impairments.28 But meanwhile, adverse general conditions of burn individuals, such as immunosuppression state which is proven prevalent at the early burn stage, would greatly hinder the healing process via deranged mediators.29–31 Since plasma transfusion of burn animals invoked the undesired conversion of remote tissues in repeated reports,4,5,7 we hypothesized that circulating exosomes restrained the repairing functions of local cells at the early burn stage, resulting in delayed healing of burn wounds and following infectibility. Furthermore, whether the tendency of the hyperplastic scar in the recovery stage of burns is a “resilience” to the previous suppressive circulating factors shall be an interesting topic of great potential.

To discover the possible blood-tissue interactions at the early stage postburn, we compared the target genes of exosome-derived DEMs with the DEGs in burn wounds. Gene expression matrices of normal skins and partial-thickness burned skins 4–7 days post thermal injuries in dataset GSE8056 were enrolled to obtain the DEGs. The overlap between 12,901 target genes of DEMs and 1861 DEGs consisted of 1058 focus genes. To exploit the leading candidates to key miRNA-mRNA axes, hub genes were rendered in Cytoscape after the PPI network construction of the focus genes. Calculated by degree scores, CDK1, CCNB1, CCNA2, BUB1B, PLK1, KIF11, AURKA, NUSAP1 and CDCA8 ranked top 9 hub genes. These hub target genes led to the identification of 16 predicted upstream DEMs, including 4 downregulated miRNAs (hsa-miR-6848-3p, has-miR-4684-3p, has-miR-4786-5p and has-miR-365a-5p) and 12 upregulated miRNAs (hsa-miR-6751-3p, hsa-miR-718, hsa-miR-4754, hsa-miR-6754-3p, hsa-miR-4739, hsa-miR-6739-5p, hsa-miR-6884-3p, hsa-miR-1224-3p, hsa-miR-6878-3p, hsa-miR-6795-3p, hsa-miR-550a-3p, and hsa-miR-550b-3p). A potential key miRNA-mRNA network in the blood-to-tissue interactions network at the early burn stage was therefore constructed.

All the hub genes are closely centered on promoting cell proliferation as summarized in the GeneCards database which is founded on numerous papers (Table 4). In contrast, function annotation of miRNAs remains a promising project requiring more endeavors. Only part of the identified key miRNAs was reported in previous studies to act as novel biomarkers in diagnosis and prognosis of various cancers. Functions of these starring miRNAs in proliferation and migration of tumor cells shall be valuable references for explorations of burn wound conversion, given the fact that tumorigenesis and the wound healing process are similarly in a close relationship to the cell cycle distribution.32,33 For downregulated key miRNAs, have has-miR-4684-3p expressions reduced in radioresistant atypical meningiomas patients comparing to radiosensitive group.34 But whether the tumor’s ability to survive the adverse environmental impact of radiotherapy was induced by downregulation of has-miR-4684-3p or its underlying mechanism was not further explored. Expression slump of has-miR-365a-5p was found in the bone tissues of patients with infected tibial non-union, and was possibly related to increased extracellular matrix degradation and endplate cartilage degeneration.35 In addition, significantly depressed has-miR-365a-5p in non-small cell lung cancer promoted the viability and migration of tumor cells.36,37

Among upregulated key DEMs, 7 miRNAs were elaborated in previous literatures. Hsa-miR-718 can serve as a tumor promoter or a suppressor in different types of cancer. Elevated has-miR-718 and subsequent lowered phosphatase and tensin homolog (PTEN) expressions caused accelerated proliferation and invasion of gastric cancer cells.38 Pro-tumor roles of has-miR-718 were also disclosed in papillary thyroid cancer39 and hepatocellular carcinoma.40 But meanwhile, the tumor-suppressor effect of has-miR-718 was confirmed in ovarian cancer as it inhibited expressions of vascular endothelial growth factor (VEGF).41 Higher has-miR-4754 levels were found in the dysembryoplastic neuroepithelial tumor42 and positive lymphatic metastasis of the primary gastric tumor.43 And an NGS of plasma RNAs discovered MIR6754 family were prominent in newly diagnosed multiple myeloma patients.44

Hsa-miR-4739 was a well-studied miRNA and predicted to have most targets in our hub genes. LncRNA-sponged has-miR-4739 enhanced tumorigenesis and progression of prostate cancer.45 In type II diabetes patients, Denis et al46 first descried urinary exosomal has-miR-4739 highlighted in diabetic nephropathy, and a microarray profiling of plasma miRNAs revealed that has-miR-4739 levels were independently associated with critical limb ischemia.47 Moreover, has-miR-4739 increased collagen-I synthesis in pleural mesothelial cells and facilitated pleural fibrosis.48 The has-miR-6884-3p/CCNB1 axis lay in the tumorigenic behaviors of hepatocellular carcinoma cells boosted by lncRNA RP11-295G20.2.49 Significant overexpressions of urinary hsa-miR-1224-3p were observed in radiation-induced renal tubular injury,50 and it showed high sensitivity for bladder cancer diagnosis.51 MiR-1224-3p/PLK1 and miR-1224-3p/ETV1 axes were modulated in circular RNA ZNF609-induced promotion of glioma and lung adenocarcinoma, respectively.52,53 Hsa-miR-550a-3p downregulated TIMP2 and resulted in promoted non-small cell lung cancer cell proliferation and metastasis,54 but it exerted antitumor effect in breast cancer by alleviating levels of ERK1 and ERK2.55

The accordantly elevated expressions of hub genes in GSE8056 indicated the rapid and vigorous regeneration process of local tissues surrounding partial-thickness wounds at the early burn stage. However, only 4 miRNA-mRNA axes (hsa-miR-6848-3p/KIF11, has-miR-4684-3p/AURKA, has-miR-4786-5p/NUSAP1 and has-miR-365a-5p/NUSAP1) exhibited reversed expression regulations of DEMs in our NGS result and DEGs in GSE8056, which conformed to the general pattern that miRNA functions as a silencer of target mRNA.56 On the contrary, upregulations of 12 other DEMs were supposed to alleviate hub gene expressions in potential blood-to-tissue interactions, offsetting the stimulated healing process of burned skins. The key miRNA-mRNA network in the possible blood-to-tissue interactions we managed to construct, may suggest that suppression of cell proliferation could outweigh the promotion of wound healing induced by plasma exosome-derived miRNAs at the early-stage after severe burn injury.

We do acknowledge that many more investigations await accomplishments to perfect our study. A larger capacity of samples and verifications of identified key miRNA-mRNA axes were strongly desired. However, limited funds as well as a paucity of eligible samples for miRNAs screening and existing burn-related databases practically restricted this account. Future studies shall mainly focus on the source and target cells of burn-induced exosomes and functional mechanisms. Hopefully our study would offer a thoughtful interpretation of plasma exosome-derived miRNA profiles and inspire a novel insight into burn wound care.

Conclusion

A small RNA sequencing revealed significant derangements of plasma exosome-derived miRNAs following severe burns. Microarray dataset GSE8056 exhibited gene expression alterations of partial-thickness burned skins. Bioinformatic analysis identified 20 key miRNA-mRNA axes in potential blood-to-tissue interactions at the early burn stage. A constructed miRNA-mRNA network suggested circulating exosomal miRNAs possibly influence the healing process of burn wounds, offering candidate targets of therapeutic interventions.