Analysis of miRNA expression profiling in human umbilical vein endothelial cells affected by heat stress.

To investigate the regulation of endothelial cell (EC) microRNAs (miRNAs) altered by heat stress, miRNA microarrays and bioinformatics methods were used to determine changes in miRNA profiles and the pathophysiological characteristics of differentially expressed miRNAs. A total of 31 differentially expressed miRNAs were identified, including 20 downregulated and 11 upregulated miRNAs. Gene Ontology (GO) enrichment analysis revealed that the validated targets of the differentially expressed miRNAs were significantly enriched in gene transcription regulation. The pathways were also significantly enriched in the Kyoto Encyclopedia of Genes and Genomes analysis, and most were cancer-related, including the mitogen-activated protein kinase signaling pathway, pathways involved in cancer, the Wnt signaling pathway, the Hippo signaling pathway, proteoglycans involved in cancer and axon guidance. The miRNA-gene and miRNA‑GO network analyses revealed several hub miRNAs, genes and functions. Notably, miR‑3613-3p played a dominant role in both networks. MAP3K2, MGAT4A, TGFBR1, UBE2R2 and SMAD4 were most likely to be controlled by the altered miRNAs in the miRNA-gene network. The miRNA‑GO network analysis revealed significantly complicated associations between miRNAs and different functions, and that the significantly enriched functions targeted by the differentially expressed miRNAs were mostly involved in regulating gene transcription. The present study demonstrated that miRNAs are involved in the pathophysiology of heat-treated ECs. Understanding the functions of miRNAs may provide novel insights into the molecular mechanisms underlying the heat‑induced pathophysiology of ECs.

Introduction

Heat stroke (HS) is a hyperthermia-associated disease that may induce multi-organ dysfunction and death (1). Abnormalities in endothelial cell (EC) function have been implicated in the pathophysiology of HS (2). EC dysfunction initiates a number of events that trigger EC activation and predispose the vessel wall to increased endothelial permeability, leukocyte adherence, endothelial apoptosis, pro-oxidation and thrombosis (3–9). MicroRNAs (miRNAs) participate in the control of EC-mediated homeostasis (10,11).

miRNAs are small non-coding RNAs of 18–22 nucleotides in length that play an important role in the regulation of gene expression at the post-transcriptional level. Specific classes of miRNAs have been demonstrated to affect the function of ECs. For example, miR-126, miR-128, the miR-17-92 cluster, and the miR-23-24-27 cluster, have been reported to control blood vessel development, growth and differentiation; therefore, these miRNAs and clusters contribute to the regulation of angiogenesis (12–16). Previous studies demonstrated that certain inducible miRNAs also contribute to the regulation of inflammation. Inducible miRNAs, including miR-31, miR-17-3p, miR-126, miR-181b, miR-10a, miR-155, miR-221/222, miR-92, miR-125a/b, miR-365 and miR-663, may participate in the regulation of inflammatory cytokine-induced leukocyte adhesion to ECs and activation of leukocytes (17–27). Recent studies also indicated that certain miRNAs, including miR-155, miR-320, miR-222, miR-125b, miR-410 and miR-218, play an important role in controlling the barrier function of ECs (28–30). These miRNAs reportedly regulate adherens junction disassembly responses, cell migration and cell morphology, which contribute to changes in vascular permeability and integrity. Another set of endothelial miRNAs participate in the control of vascular tone. Several reports have indeed provided evidence that miR-155, miR-221, miR-222, miR-125a/b and miR-21 are associated with the regulation of blood pressure and blood flow by releasing vasodilators or vasoconstrictors (24,31–33). Upon exposure to certain stimuli, ECs may secrete miRNAs that mediate cell-to-cell communication. For example, ECs have been shown to produce microvesicles under shear stress. These microvesicles are enriched in miRNAs, particularly miR-143/145 or miR-150, which exert an atheroprotective effect on smooth muscle cells (34–36). This evidence suggests that miRNAs are important molecules that modulate the biological processes of ECs.

Although the involvement of individual miRNAs in endothelial regulation has been characterized in several diseases, changes in the general miRNA expression profiles under heat stress have not been reported. The aim of the present study was to systematically analyze the profiles of miRNA expression in human umbilical vein ECs (HUVECs) exposed to acute heat stress using microarray and bioinformatics analyses, in order to evaluate the potential role of miRNAs in post-transcriptional regulation.

Materials and methods

HUVEC culture and treatment with heat

HUVECs (American Type Culture Collection, Manassas, VA, USA) were cultured in complete extracellular matrix supplemented with 5% fetal bovine serum (Invitrogen, Thermo Fisher Scientific, Inc., Waltham, MA, USA), 1% EC growth supplement (PromoCell, Heidelberg, Germany), 100 U/ml penicillin and 100 µg/ml streptomycin (Invitrogen, Thermo Fisher Scientific, Inc.), in a humidified atmosphere with 5% CO2 at 37°C. All the experiments were performed using 4–6 passages of HUVECs.

In the following experiments, the cells were randomly divided into two groups

The control cells were always placed in an incubator at 37°C. The cells in the heat treatment group were cultured in the incubator at 43°C for 1 h, and then at 37°C for a further 24 h.

RNA extraction

Total RNA from each sample was individually isolated using TRIzol (Thermo Fisher Scientific, Inc., Carlsbad, CA, USA) and miRNeasy mini kit (Qiagen, Inc., Valencia, CA, USA) according to the manufacturers' instructions. This procedure efficiently recovered all RNA species, including miRNAs. RNA quantity and quality were measured using a NanoDrop spectrophotometer (ND-1000; NanoDrop Technologies, Wilmington, DE, USA) and RNA integrity was detected by gel electrophoresis.

RNA labeling and array hybridization

The isolated RNA was labeled using the miRCURY™ Hy3™/Hy5™ Power Labeling kit (Exiqon A/S, Vedbaek, Denmark) according to the manufacturer's instructions. Each RNA sample (1 µg) was 3′-end-labeled with Hy3™ fluorescent label, using T4 RNA ligase. Following the array manual, the Hy3™-labeled samples were then hybridized to the miRCURY™ LNA array (v.18.0) (Exiqon A/S). The total mixture with hybridization buffer was first denatured for 2 min at 95°C, incubated on ice for 2 min, and then hybridized to the microarray for 16–20 h at 56°C in a 12-Bay Hybridization system (Nimblegen Systems, Inc., Madison, WI, USA). After the hybridization procedure was terminated, the slides were washed several times using the wash buffer kit (Exiqon A/S), and finally dried by centrifugation for 5 min at 340 × g. The slides were then scanned with the Axon GenePix 4000B microarray scanner (Axon Instruments, Foster City, CA, USA), which contains 3,100 capture probes annotated in miRBase 18.0, covering all human, mouse and rat miRNAs, as well as all viral miRNAs related to these species. In addition, the microarray scanner contains capture probes for 25 additional new miRPlus™ human miRNAs that are proprietary and not found in miRBase.

Array data analysis

Scanned images were imported into GenePix Pro 6.0 software (Axon Instruments) for grid alignment and data extraction. Replicated miRNAs were averaged, and miRNAs with intensities ≥30 were selected for calculation of the normalization factor. Expressed data were normalized using the median normalization. Following normalization, significantly differentially expressed miRNAs were identified through volcano plot change filtering. Finally, unsupervised hierarchical clustering was performed based on the normalized data of miRNAs in all samples using Multi Experiment Viewer software, version 4.6 (The Institute for Genomic Research, Rockville, MD, USA). A miRNA was defined as differentially expressed between the compared groups if the P-value was <0.05 and the fold change was >2.

Reverse transcription-quantitative polymerase chain reaction (RT-qPCR)

To validate the authenticity of miRNA expression screened by microarray assay, six differentially expressed miRNAs were analyzed by RT-qPCR. For this analysis, total RNA (10 ng) was reversely transcribed into cDNA using the TaqMan® MicroRNA Reverse Transcription kit (Applied Biosystems, Foster City, CA, USA) according to the manufacturer's instructions. qPCR was performed using SYBR-Green (Thermo Fisher Scientific, Inc.) according to the manufacturer's instructions in the ABI 7300 real-time qPCR system (Bio-Rad Laboratories, Inc., Hercules, CA, USA). The quantification of each miRNA was performed in relation to the U6 housekeeping miRNAs that was used as a standard. The relative abundance of each miRNA was calculated by the comparative Cq method (2−ΔΔCq), and the results were assessed by using a t-test.

Bioinformatics analysis

Potential target genes of differentially expressed miRNAs were predicted using the three most popular databases, namely, TargetScan, miRanda and miRDB. Only the overlapping genes that were in all three databases were considered as the target genes of the differentially expressed miRNAs. The predicted genes subsequently underwent Gene Ontology (GO) analysis for functional annotation analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) database analysis to identify the involved enriched pathways. Fisher's exact and Chi-square tests were used to determine the significance of the GO terms and pathways, and the false discovery rate (FDR) was calculated to correct the P-value (37). Only the GO terms and pathways with an adjusted P-value of <0.05 and an FDR of <0.05 were selected.

miRNA-gene-network analyses

The target genes identified in the functional and pathway enrichment analyses were investigated. The 494 overlapping genes that were in both databases were selected. In addition, to improve the understanding of the associations between the key miRNAs and hub genes, networks for miRNA target genes were constructed according to the regulatory associations between miRNAs and target genes. The associations of the genes and miRNAs were evaluated by graph theory. The degree was used to measure the regulated effect of miRNAs on genes.

miRNA-GO-network analysis

A miRNA-GO-network was constructed according to the specifically expressed miRNAs and the results of the GO analysis. The control status of GO and miRNAs in the network was evaluated by graph theory, and the degree of miRNAs indicated the number of GO which was regulated by miRNAs in the network. In the same way, a higher degree of GO indicated that more miRNAs contributed to GO.

Statistical analysis

All quantitative data are presented as mean ± standard deviation. Comparisons between the two groups were performed using the Student's t-test. Comparisons of multiple group data were performed using one-way analysis of variance followed by Turkey's post hoc test. P-values <0.05 were considered statistically significant. Statistical values were calculated using SPSS software, version 20.0 (IBM Corp., Armonk, NY, USA).

Results

Screening differentially expressed miRNAs

To identify significantly differentially expressed miRNAs between heat-treated cells and controls, the screening threshold was set to a fold change of >2 and a P-value of <0.05. In total, 31 miRNAs (20 downregulated and 11 upregulated) were identified (Table I).

Table I

Differentially expressed microRNAs (miRNAs) between heat-treated cells and controls (P<0.05, fold change >2).

miRNA	Fold change	Type of regulation	P-value
hsa-miR-4448	0.240462612	Down	<0.05
hsa-miR-3180	0.243442001	Down	<0.05
hsa-miR-1231	0.309725485	Down	<0.05
hsa-miR-3620	0.330975034	Down	<0.05
hsa-miR-665	0.331296255	Down	<0.05
hsa-miR-4656	0.381063941	Down	<0.05
hsa-miR-4458	0.384188702	Down	<0.05
hsa-miR-3613-3p	0.385431663	Down	<0.05
hsa-miR-486-5p	0.391376978	Down	<0.05
hsa-miR-3180-3p	0.404438992	Down	<0.05
hsa-miR-378h	0.416819304	Down	<0.05
hsa-miR-1296	0.450624607	Down	<0.05
hsa-miR-616-3p	0.453329075	Down	<0.05
hsa-miR-4486	0.453329201	Down	<0.05
hsa-miR-4749-5p	0.46624915	Down	<0.05
hsa-miR-489	0.468301966	Down	<0.05
hsa-miR-3130-5p	0.469450744	Down	<0.05
hsa-miR-411-5p	0.470553479	Down	<0.05
hsa-miR-299-3p	0.476023136	Down	<0.05
hsa-miR-4665-5p	0.492758274	Down	<0.05
hsa-miR-1290	2.07284178	Up	<0.05
hsa-miR-4430	2.135374336	Up	<0.05
hsa-miR-4655-5p	2.191798238	Up	<0.05
hsa-miR-140-5p	2.318002303	Up	<0.05
hsa-miR-4500	2.324482718	Up	<0.05
hsa-miR-485-3p	2.399120465	Up	<0.05
hsa-miR-589-5p	2.399121462	Up	<0.05
hsa-miR-34c-5p	2.518895027	Up	<0.05
hsa-miR-3158-5p	3.033915491	Up	<0.05
hsa-miR-3613-5p	3.690735374	Up	<0.05
hsa-miR-1281	4.409294993	Up	<0.05

Validation by RT-qPCR

Based on their expression levels and fold change in expression, 3 downregulated miRNAs (miR-1281, miR-34c-5p and miR-140-5p) and 3 upregulated miRNAs (miR-4448, miR-1231 and miR-4656) were selected for validation by RT-qPCR. The results revealed a significant intra-group difference in the expression of miRNAs, in a manner consistent with the data obtained from the microarray analysis (P<0.05) (Fig. 1).

Figure 1

Validation of microRNA (miRNA) microarray results by reverse transcription-quantitative polymerase chain reaction. (A) Downregulated miRNAs in the heat-treated group; (B) upregulated miRNAs in the heat-treated group. All miRNAs exhibited statistically significant changes in expression in the experimental group compared with the control group. *P<0.05.

Target gene and bioinformatics analysis

A total of 2,024 target genes were identified for the validated miRNAs. To better understand the potential implications of the differentially expressed miRNAs, these target genes were subjected to GO analysis to evaluate their potential functions. In the present study, the top GO terms of the target genes of the downregulated miRNAs were transcription (DNA-dependent), regulation of transcription (DNA-dependent), positive regulation of transcription from RNA polymerase II promoter, positive regulation of transcription (DNA-dependent), and negative regulation of transcription from RNA polymerase II promoter. The top GO terms of the target genes of the upregulated miRNAs were transcription (DNA-dependent), positive regulation of transcription (DNA-dependent), positive regulation of transcription from RNA polymerase II promoter, nervous system development, and regulation of transcription (DNA-dependent) (Fig. 2).

Figure 2

(A) Top 20 functional Gene Ontology (GO) terms of the predicted downregulated microRNA (miRNA) target genes. (B) Top 20 functional GO terms of the predicted upregulated miRNA target genes.

The target genes were then subjected to pathway enrichment analysis using KEGG pathways to determine the canonical pathways controlled by the identified miRNAs. Pathways involved in cancer, proteoglycans involved in cancer, the Wnt signaling pathway, axon guidance and Hippo signaling pathways, were the most active pathways regulated by the downregulated miRNAs. The mitogen-activated protein kinase (MAPK) signaling pathway, pathways involved in cancer, the Wnt signaling pathway, the Hippo signaling pathway and proteoglycans involved in cancer, were the top pathways regulated by the upregulated miRNAs (Fig. 3).

Figure 3

(A) Top 20 signaling pathways predicted to be regulated by downregulated microRNAs (miRNAs). (B) Top 20 signaling pathways predicted to be regulated by upregulated miRNAs. PI3K, phosphoinositide-3 kinase; MAPK, mitogen-activated protein kinase; HTLV, human T-lymphotropic virus; TGF, transforming growth factor.

A miRNA-gene-network was constructed according to the results of the GO and pathway analyses. The core miRNAs of the interaction network included miR-3613-3p, miR-34c-5p and miR-485-3p (Table II). The network also revealed that MAPK kinase kinase 2 (MAP3K2), α-1,3-mannosyl-glycoprotein 4-β-N-acetylglu cosaminyltransferase A (MGAT4A), transforming growth factor β receptor 1 (TGFBR1), ubiquitin-conjugating enzyme E2 R2 (UBE2R2), and Sma- and Mad-related family 4 (SMAD4), were the most crucial target genes (Table III). The associations of miRNAs with genes are shown in Fig. 4.

Table II

Crucial microRNAs (miRNAs) in the miRNA-gene-network (degree >10).

miRNA	Degree
hsa-miR-3613-3p	237
hsa-miR-34c-5p	43
hsa-miR-485-3p	43
hsa-miR-4458	34
hsa-miR-4500	31
hsa-miR-140-5p	30
hsa-miR-1290	30
hsa-miR-665	27
hsa-miR-3158-5p	24
hsa-miR-3130-5p	14
hsa-miR-486-5p	13
hsa-miR-489	13
hsa-miR-616-3p	11
hsa-miR-3620	11
hsa-miR-4448	11
hsa-miR-378h	11

Table III

Crucial target genes in the microRNA-gene-network (degree >3).

Gene symbol	Gene description	Degree
MAP3K2	Homo sapiens mitogen-activated protein kinase kinase kinase 2 (MAP3K2), mRNA	5
MGAT4A	Homo sapiens mannosyl (α-1,3-)-glycoprotein β-1,4-N-acetylglucosaminyltransferase, isozyme A (MGAT4A), transcript variant 1, mRNA	4
TGFBR1	Homo sapiens transforming growth factor, β receptor 1 (TGFBR1), transcript variant 1, mRNA	4
UBE2R2	Homo sapiens ubiquitin-conjugating enzyme E2 R2 (UBE2R2), mRNA	4
SMAD4	Homo sapiens SMAD family member 4 (SMAD4), mRNA	4
PRPF40A	Homo sapiens PRP40 pre-mRNA processing factor 40 homolog A (S. cerevisiae) (PRPF40A), mRNA	3
PARVA	Homo sapiens parvin, α (PARVA), mRNA	3
ACSL4	Homo sapiens acyl-CoA synthetase long-chain family member 4 (ACSL4), transcript variant 2, mRNA	3
ENAH	Homo sapiens enabled homolog (Drosophila) (ENAH), transcript variant 2, mRNA	3
JAG1	Homo sapiens jagged 1 (JAG1), mRNA	3
PDE4D	Homo sapiens phosphodiesterase 4D, cAMP-specific (PDE4D), transcript variant 2, mRNA	3
E2F5	Homo sapiens E2F transcription factor 5, p130-binding (E2F5), transcript variant 1, mRNA	3
DCC	Homo sapiens deleted in colorectal carcinoma (DCC), mRNA	3
PPM1A	Homo sapiens protein phosphatase, Mg2+/Mn2+-dependent, 1A (PPM1A), transcript variant 3, mRNA	3
IGF1R	Homo sapiens insulin-like growth factor 1 receptor (IGF1R), mRNA	3
CDC25A	Homo sapiens cell division cycle 25 homolog A (S. pombe) (CDC25A), transcript variant 2, mRNA	3
CPEB1	Homo sapiens cytoplasmic polyadenylation element binding protein 1 (CPEB1), transcript variant 1, mRNA	3
NLK	Homo sapiens nemo-like kinase (NLK), mRNA	3
FZD3	Homo sapiens frizzled family receptor 3 (FZD3), transcript variant 2, mRNA	3
PPP1CB	Homo sapiens protein phosphatase 1, catalytic subunit, β isozyme (PPP1CB), transcript variant 3, mRNA	3
UBOX5	Homo sapiens U-box domain containing 5 (UBOX5), transcript variant 2, mRNA	3
BBC3	Homo sapiens BCL2-binding component 3 (BBC3), nuclear gene encoding mitochondrial protein, transcript variant 4, mRNA	3
ACVR1C	Homo sapiens activin A receptor, type IC (ACVR1C), transcript variant 1, mRNA	3

Figure 4

miRNA-gene-network showing the associations between key microRNAs (miRNAs) and the relevant hub genes. The square nodes represent miRNAs (red nodes, upregulated miRNAs; blue nodes, downregulated miRNAs), and the circular nodes represent target genes. The lines indicate the regulatory effect of miRNAs on target genes. The size of the square nodes indicates the degree to which miRNAs contribute to the network.

The network analysis was useful for determining regulatory associations between the key miRNAs and hub GO. In this network, miR-3613-3p, miR-4458 and miR-4500, which contributed more than the other specifically expressed miRNAs, exhibited 479, 217 and 216 GO functions, respectively (Table IV). The most significantly regulated function-cluster of the total 20 categories was gene expression, whereas others included positive regulation of transcription, positive regulation of transcription by RNA polymerase II promoter and signal transduction (Table V). The significantly complicated associations of the top 20 functions with miRNAs are shown in Fig. 5.

Table IV

Crucial microRNAs (miRNAs) in the miRNA-GO-network (degree >100).

miRNA	Degree
hsa-miR-3613-3p	479
hsa-miR-4458	217
hsa-miR-4500	216
hsa-miR-665	182
hsa-miR-1290	180
hsa-miR-140-5p	173
hsa-miR-34c-5p	173
hsa-miR-485-3p	165
hsa-miR-486-5p	141
hsa-miR-3158-5p	135

GO, Gene Ontology.

Table V

Crucial GO functions in the microRNA-GO-network (degree >15).

GO (name)	Degree
Gene expression	21
Positive regulation of transcription, DNA-dependent	21
Positive regulation of transcription from RNA polymerase II promoter	20
Signal transduction	20
Negative regulation of transcription from RNA polymerase II promoter	19
Small molecule metabolic process	19
Transforming growth factor-β receptor signaling pathway	19
Axon guidance	18
Innate immune response	18
Transcription, DNA-dependent	18
Virus-host interaction	17
Negative regulation of apoptotic process	16
Negative regulation of transcription, DNA-dependent	16
Regulation of transcription, DNA-dependent	16
Synaptic transmission	16
Transcription initiation from RNA polymerase II promoter	16
Heart development	15
Mitotic cell cycle	15
Positive regulation of cell proliferation	15
Protein phosphorylation	15

GO, Gene Ontology.

Figure 5

Interaction between the top 20 targeted functions and the differentially expressed microRNAs (miRNAs). The red square nodes represent upregulated miRNAs, the blue square nodes represent downregulated miRNAs, and the violet circular nodes represent target functions. The lines represent the inhibitory effect of miRNAs on target functions. The size of square nodes indicates the degree to which the miRNAs contribute to the network.

Discussion

In the present study, the miRNA expression profiles of heat-treated HUVECs were evaluated by microarray analysis. A number of miRNAs were obtained from both heat-treated cells and normal controls, and 31 differentially expressed miRNAs were selected, including 20 downregulated and 11 upregulated miRNAs (Table I). The target genes of 31 differentially expressed miRNAs were predicted, and 2,024 target genes were validated. GO and pathway enrichment analysis revealed the significantly different functions and pathways targeted by the differentially expressed miRNAs. A total of 755 notable GO terms and 136 pathways were identified. The most active functions and pathways are shown in Figs. 2 and 3. The wide target range indicates that these miRNAs may play important roles in the biological regulation of HUVECs under heat stress. These results also demonstrated that the differentially expressed miRNAs were predicted to control various functions that were most relevant to the transcription-related process. Moreover, the miRNAs regulated gene expression in ECs under heat treatment via multiple signaling pathways, and biological regulatory mechanisms were integrated in several main pathways. The majority of the pathways were cancer-related. Following exposure to environmental stressors or cytokines, these pathways may be activated to regulate cell proliferation, migration, differentiation, evasion and angiogenesis.

The network analysis linked the differentially regulated miRNAs to target genes and the relevant biological functions, and facilitated a better understanding of the role of dysregulated miRNAs. In the networks, certain miRNAs functioned as network hubs, such as miR-3613-3p, miR-34c-5p, miR-485-3p, miR-4458 and miR-4500. Of these miRNAs, miR-3613-3p was the most active miRNA in the networks, targeting 237 genes and 479 biological functions (Tables II and IV and Figs. 4 and 5). Previous studies have not focused on miR-3613-3p, and it has not been previously described as participating in EC physiology and pathology. However, plasma miR-3613-3p has been reported as a specific biomarker of lung adenocarcinoma (38). Another study indicated that circulating miR-3613-3p may predict persistent severe axial pain following exposure to common stressors, and may be involved in the pathogenesis of post-traumatic musculoskeletal pain (39). Liu et al found that miR-3613-3p was upregulated in left atrial appendages in atrial fibrillation (40). Moreover, Wang et al reported that urinary expression of miR-3613-3p was downregulated in patients with immunoglobulin A nephropathy, and the levels of miR-3613-3p were correlated with disease severity (41). Therefore, although the specific mechanism of function of miR-3613-3 remains unknown, decreased miR-3613-3p expression may exert a modulatory effect on ECs under heat stress. Additional studies are required to elucidate the detailed role of miR-3613-3p in ECs.

The miRNA-gene network analysis revealed a range of hub genes, such as MAP3K2, MGAT4A, TGFBR1, UBE2R2 and SMAD4 (Table III). MAP3K2 encodes MAPK kinase kinase 2 (MEKK2), which is a component of a protein kinase signal transduction cascade that preferentially regulates the c-Jun N-terminal kinases (JNKs) and extracellular signal-regulated kinase 5 (ERK5) pathways by phosphorylating and activating MAP2K5 as well as MAP2K7 in a similar manner (42,43). The MAP2K5/ERK5 pathway is required for normal cardiovascular development and vascular integrity, improves EC viability and reduces apoptosis (44,45), whereas JNK MAPKs promote apoptosis of ECs under most conditions (46). MEKK2 coordinately activates signaling through MEK5/ERK5 and MEK7/JNK protein kinases to respond to different stimuli. MEKK2 has also been shown to directly phosphorylate and activate IκB kinases and, thus, plays a role in the nuclear factor-κB (NF-κB) signaling pathway (47,48). MEKK2 has also been found to bind and activate protein kinase C-related kinase 2, which suggests its involvement in a regulated signaling process (49). MGAT4A encodes the glycosyltransferase N-acetylglucosaminyltransferase IVa, and catalyzes the transfer of N-acetylglucosamine (GlcNAc) from UDP-GlcNAc to the GlcNAcβ1-2Mana1-3 arm of the N-glycan core via a β1-4 linkage (50,51). As an important post-translational modification of target proteins, N-glycosylation plays an important role in cell proliferation, differentiation, apoptosis and development, during embryogenesis as well as in mature tissues (52–54). The functions of N-glycans associated with cell surface receptors are of great importance, since N-glycans are directly involved in controlling cellular functions. The target proteins of glycosyltransferases include epidermal growth factor receptor, human epidermal growth factor receptor 3, transforming growth factor-β (TGF-β) receptor, T-cell receptors, β-site APP cleaving enzyme 1, glucose transporter 2, E-cadherin and α5β1 integrin (52). TGFBR1 encodes TGF-β type I receptor (TβRI), which is the sole known cell surface receptor serine/threonine kinase in humans and participates in TGF-β signaling. The binding of TGF-β to TβRII dimers stabilizes the interaction of TβRII dimers with two TβRI molecules to form a ligand receptor-complex involving a heterotetrameric receptor and a ligand dimer (55). Ligand binding to TβRII causes the phosphorylation of TβRI, which activates Smad2 and Smad3 (55). The activated Smad proteins form complexes with Smad4 and are translocated into the nucleus, where they regulate the transcription of target genes (56,57). TGF-β1 exerts its effects on a variety of cells through activation of the Smad proteins and other signaling pathways, including MAPKs, small GTPases, protein kinase C and phosphoinositide 3-kinase (58,59). Under most conditions, the activa tion of the TGF-β1 signaling pathway inhibits EC proliferation and migration and induces EC apoptosis (60). Smad4 is a central mediator of TGF-β1/Smad signaling. Smad4 belongs to the evolutionarily conserved family of Smad proteins that are crucial for transmitting TGF-β1 superfamily signals from the cell surface to the nucleus (61). UBE2R is involved in the ubiquitin cascade. The ubiquitin-proteasome system (UPS) comprises ubiquitin, ubiquitin-activating enzyme (E1), ubiquitin-binding enzyme (E2), ubiquitin-ligating enzymes (E3) and the 26S proteasome (62). Ubiquitin is attached to substrates by a complex three-step enzymatic cascade, which involves E1 ubiquitin activation (62,63), E2 ubiquitin conjugation and a variety of E3 ubiquitin-ligating enzymes (64,65). The UPS plays an important role in the regulation of cellular processes via degradation of damaged proteins and the control of protein quantity in cells (66,67). Although the common effect of ubiquitylation is degradation, the degradation-independent Ub-signaling system also affects various biological processes, which include cytokine production, proliferation, apoptosis, inflammation, canalization, pain sensitization to host defense by regulation of protein kinase activation, signal transduction, membrane transport, transcriptional regulation, gene silencing, DNA injury repair and DNA replication (62,68). All these genes receive physicochemical stimuli and inputs from a diverse array of receptors and coordinate the activation of multiple responses.

In addition to the miRNA-gene network analysis, the miRNA-GO network analysis provides a more detailed view of the integrated association between the regulation of miRNAs and heat stress. Among the top 10 functions, gene expression was the most active and was regulated by 21 differentially expressed miRNAs. Furthermore, 4 transcription-related functions were regulated by 78 differentially expressed miRNAs (Table V). These findings indicated that heat stress may affect EC gene expression directly via differentially expressed miRNAs, or induce alterations in transcription-related functions indirectly via miRNA profile changes. This finding was not surprising, as most of the hub genes participate in gene expression and transcription via certain downstream molecules. For example, MAP3K2 was one of the top genes revealed by the miRNA-gene network. The most significant role of ERK5, which is a downstream target of MAP3K2, appears to be the regulation of a number of downstream transcription factors, including the myocyte-enhancer factor (MEF) family, MEF2A, C and D, which are the best characterized members of this family (69,70). ERK5 also directly enhances the transcription of c-Myc, cAMP-response element-binding protein, and Sap1a (71–73). Enhancement of Sap1a occurs via a serum response element that may also be involved in activating the c-Fos promoter (74). JNKs, another target of MAP3K2, translocate from the cytoplasm to the nucleus once activated, and regulate a series of transcription factors, such as c-Jun, Ets-like protein-1, activating transcription factor 2, p53, c-Myc, and nuclear factor of activated T-cells (75–79).

In the present study, the miRNA-GO network analysis revealed a significantly complicated association between miRNAs and different functions (Fig. 5 and Tables IV and V). First, more functions than genes that were targeted by miRNAs were found. Second, each differentially expressed miRNA modulated a large variety of functions, and each function was controlled by several miRNAs, thus resulting in a complicated network. For example, in the present study, miR-3613-3p targeted 479 functions, and the function of gene expression was regulated by 21 differentially expressed miRNAs (Tables IV and V). Third, one miRNA may simultaneously regulate opposing functions. For example, certain miRNAs, such as miR-3613-3p, miR-4458 or miR-4500, were not only involved in the positive regulation of transcription, but also in the negative regulation of transcription (data not shown). Finally, certain miRNAs contributed to the positive regulatory function, and other miRNAs contributed to the corresponding negative regulatory function. For example, miR-665 participated in the negative regulation of the apoptotic process, but miR-4500 not only participated in the anti-apoptotic process, but also promoted the cellular apoptosis process (data not shown). Moreover, each function was affected by both upregulated and downregulated miRNAs (data not shown). Our results suggested that the functions regulated by miRNAs were not an isolated event, and this regulation depends on the balance between the miRNAs that promote the function and the miRNAs that suppress the function. The miRNA-GO network provided a global view of the molecular pathophysiology of EC heat stress and helped identify the key 'players' in the whole network.

There were several limitations to this study. First, we were limited by the use of HUVECs. Although the baseline miRNA expression levels are similar between ECs from different vascular beds, some unique miRNA phenotypic diversity is present among the different types of ECs (80). As only one cell type was assessed, 'cell-specific' miRNAs could not be identified in our analysis, which could be determined using other cell types. Second, the number of samples utilized for the miRNA microarray was relatively small and, therefore, less predominantly differentiated miRNAs were difficult to detect with the microarrays due to low power. Third, the miRNA microarray technique had certain limitations, as some of the microarray results may have been affected by technical artifacts introduced by sample selection, the hybridization procedure, or the microarray platform.

In conclusion, this study demonstrated that 31 miRNAs were differentially expressed in ECs under heat stress. The bioinformatics analysis predicted the target genes, functions and signaling pathways of these miRNAs. The miRNA-gene and miRNA-GO network analyses revealed a range of hub miRNAs, genes and functions, which may consolidate our understanding of the involvement of miRNAs in the pathophysiology of heat-treated ECs. miR-3613-3p may be involved in heat stress by regulating relevant genes and functions. Although further mechanistic investigation is required, the findings presented herein provide novel insights into the molecular mechanisms of heat stress in ECs.