MicroRNAs Clustered within the 14q32 Locus Are Associated with Endothelial Damage and Microparticle Secretion in Bicuspid Aortic Valve Disease.

Background: We previously described that PECAM+ circulating endothelial microparticles (EMPs) are elevated in bicuspid aortic valve (BAV) disease as a manifestation of endothelial damage. In this study, we hypothesized that this endothelial damage, is functionally related to the secretion of a specific pattern of EMP-associated miRNAs.
Methods: We used a bioinformatics approach to correlate the PECAM+ EMP levels with the miRNA expression profile in plasma in healthy individuals and BAV patients (n = 36). In addition, using the miRNAs that were significantly associated with PECAM+ EMP levels, we inferred a miRNA co-expression network using a Gaussian graphical modeling approach to identify highly co-expressed miRNAs or miRNA clusters whose expression could functionally regulate endothelial damage.
Results: We identified a co-expression network composed of 131 miRNAs whose circulating expression was significantly associated with PECAM+ EMP levels. Using a topological analysis, we found that miR-494 was the most important hub within the co-expression network. Furthermore, through positional gene enrichment analysis, we identified a cluster of 19 highly co-expressed miRNAs, including miR-494, that was located in the 14q32 locus on chromosome 14 (p = 1.9 × 10-7). We evaluated the putative biological role of this miRNA cluster by determining the biological significance of the genes targeted by the cluster using functional enrichment analysis. We found that this cluster was involved in the regulation of genes with various functions, specifically the "cellular nitrogen compound metabolic process" (p = 2.34 × 10-145), "immune system process" (p = 2.57 × 10-6), and "extracellular matrix organization" (p = 8.14 × 10-5) gene ontology terms and the "TGF-β signaling pathway" KEGG term (p = 2.59 × 10-8). Conclusions: Using an integrative bioinformatics approach, we identified the circulating miRNA expression profile associated with secreted PECAM+ EMPs in BAV disease. Additionally, we identified a highly co-expressed miRNA cluster that could mediate crucial biological processes in BAV disease, including the nitrogen signaling pathway, cellular activation, and the transforming growth factor beta signaling pathway. In conclusion, EMP-associated and co-expressed miRNAs could act as molecular effectors of the intercellular communication carried out by EMPs in response to endothelial damage.

Introduction

Bicuspid aortic valve (BAV), the most common cardiac congenital malformation (occurs in 1–2% of the population), is associated with valve dysfunction and is a risk factor for aortopathy (Tzemos et al., 2008; Alegret et al., 2013). The progressive dilation of the aorta, if untreated, can lead to fatal consequences such as aortic dissection and/or rupture. The mechanisms that underlie aortic dilatation have been a matter of debate for years (Padang et al., 2013); the proposed causes include anomalous flow in the ascending aorta generated by the anomalous dynamics of BAV (Kim et al., 2012; Bissell et al., 2013) and genetic causes responsible for the anomalous structure of the aortic media (Biner et al., 2009; Pepe et al., 2014).

In a previous study, we reported that circulating endothelial microparticles (EMPs) are elevated in BAV patients and related to aortic dilation (Alegret et al., 2016). Circulating EMPs are a type of extracellular microvesicle (100–1,000 nm) that bud directly from the plasma membranes of endothelial cells upon activation, injury, or apoptosis and that are involved in intercellular communication. Specifically, we identified PECAM+ EMPs, which are a kind of EMP that express CD31 and are released in endothelial damage, as those related to BAV disease.

Micro-ribonucleic acids (miRNAs) are endogenously expressed, 19- to 23-nt-long noncoding RNAs that regulate gene expression at the post-transcriptional level, mostly via imperfect base-pairing interactions that occur preferentially within the 3′ untranslated regions (UTRs) of target mRNAs. MiRNA genes are distributed across diverse genomic locations, and although some miRNAs are isolated, ~50% are found in clusters transcribed as polycistronic miRNA transcripts (Mourelatos et al., 2002). MiRNAs are considered potent post-transcriptional regulators because each miRNA has multiple to several 100 target genes; therefore, inhibition of a single miRNA can lead to the activation of multifactorial physiological processes. In addition to functioning intracellularly, miRNAs can be exported or released by cells into the circulating blood in very stable forms. Microparticles are reportedly the major carriers of miRNAs in the blood (Diehl et al., 2012). MiRNA signatures have been proposed as potentials with the potential to improve disease diagnosis and prognosis in clinical practice and have been identified as useful biomarkers for a wide range of cardiovascular diseases. For example, in a previous study, we proposed miR-122, miR-130a, miR-486, and miR-718 as molecular features associated with BAV and aortic dilation (Martínez-Micaelo et al., 2017).

Currently, there are no effective strategies to prevent the progression of BAV disease, including the aortic dilation, and the development of new strategies requires more detailed understanding of the molecular mechanisms associated with BAV and progressive dilation of the aorta. Therefore, in this study, we hypothesized that changes in blood flow in the ascending aorta caused by the bicuspid morphology of the aortic valve would induce endothelial damage, resulting in the secretion of a specific signature of EMP-associated miRNAs. We used a bioinformatics approach to integrate the PECAM+ EMP levels with the miRNA expression profile in plasma of healthy individuals and BAV patients. In addition, from the miRNAs that were significantly associated with PECAM+ EMP levels, we inferred a miRNA co-expression network using a Gaussian graphical modeling approach.

Methods

Study population

The patients included in this study belonged to a cohort of BAV patients who were prospectively included and followed-up in our facilities. Upon enrolment, the participants were informed and prospectively entered into a specific database, underwent a blood draw and provided informed written consent. The samples were stored in our biological samples bank (Biobanc IISPV—HUSJR) until they were needed. A BAV diagnosis was made when two aortic leaflets were clearly visualized, with or without a raphe, in the parasternal short-axis view of a transthoracic echocardiogram (Alegret et al., 2005), on a transesophageal echocardiogram (Alegret et al., 2005), or on a cardiac magnetic resonance image (Gleeson et al., 2008). Explorations were performed or supervised by the same observer (JMA). Our database and biobank also included a group of healthy tricuspid aortic valve (TAV) controls. This study was approved by the Institutional Review Board (the Clinical Ethics Committee) of our institution.

This study was designed in two evaluations in which EMP levels and the miRNA expression profile were determined in two independent cohorts of TAV individuals and BAV patients with or without aortic dilation (n = 60). In the first evaluation, plasma from the patients diagnosed with BAV, with or without aortic dilation, and healthy TAV control subjects (n = 24) was used to determine the levels of EMPs and the miRNA profile using microarrays. To improve the study's power, for this stage, we selected groups composed of patients with characteristics that were extremely homogeneous (Supplementary Table 1) to exclude possible confounding factors, such as age, sex, or BMI. We included subjects in whom high levels of circulating PECAM+ EMPs were expected, the BAV patients, and subjects in whom low levels of PECAM+ EMPs were expected, the healthy TAV controls, to study miRNAs whose expression levels could be related to PECAM+ EMPs. Patients diagnosed with cardiovascular diseases, Marfan syndrome, aortic stenosis, hypertension, or diabetes mellitus or who were receiving pharmacologic treatment (including statins, ACE/ARII, and/or β-blockers) were excluded. In the second evaluation, the EMP-associated miRNA candidates were validated by RT-qPCR in a new cohort (n = 36) of TAV healthy controls and BAV patients with or without aortic dilation (Supplementary Table 2).

Blood sampling

Blood samples were collected under overnight fasting conditions and were processed within 90 min after collection. The samples were centrifuged at 1,500 g for 15 min to obtain plasma, which was further centrifuged at 4,000 g for 10 min to obtain platelet-poor plasma. The plasmas were stored at −80°C in our biological samples bank (Biobanc IISPV—HUSJR) until they were needed.

Determination of circulating levels of EMPs

The circulating PECAM+ EMP levels were determined in the same cohorts used for the microarray and RT-qPCR validation analyses that were previously phenotyped (Alegret et al., 2016), but for this study, individuals were selected to maximize the power of the miRNA analysis. The levels of EMPs were characterized based on the presence of endothelial-specific surface antigens, the composition of which depends on the cellular origin of the microparticles and the generating process (Jimenez et al., 2003). In our previous study (Alegret et al., 2016), we determined that the presence of CD31 (PECAM) is a marker that can discriminate between microparticles released from endothelial cells subjected to endothelial damage produced by haemodynamic causes due to the anomalous aortic flow associated with BAV and microparticles released by other triggering stimuli, including cell activation or apoptosis, in TAV and BAV patients. The concentration of circulating PECAM+ EMPs was determined on an EPICS-XL (Beckman Coulter) flow cytometer at a low rate setting and a 30 s stop time. The Nano Fluorescent Particle Size Standard Kit (Spherotech) was used for instrument standardization, and Flow-Count fluorospheres (Beckman Coulter) were added as an internal calibrator to calculate microparticle amounts.

Plasma EMPs were labeled by incubating 50 μl of platelet-poor plasma with the corresponding antibody, anti-CD31-PE (Beckman Coulter), anti-CD42b-FITC (Beckman Coulter), or anti-CD45-PE (Beckman Coulter), at room temperature in the dark for 20 min as previously described (Ci et al., 2013). Then, 500 μl of PBS was subsequently added, and the EMP levels were determined as previously described (Jimenez et al., 2003; Sutherland et al., 2010; Ci et al., 2013). EMPs were defined as particles >0.1 and < 1 μm in size, and their endothelial origin was identified based on their affinity to specific cell surface antigens, namely, CD31 and CD42b. To evaluate the extent of possible contamination with leukocyte-derived microparticles, the circulating levels of CD31+CD45+ microparticles were determined in all samples. We found that < 4.5% of CD31+ microparticles co-expressed CD45+; this result is consistent with previous reports by other authors (Amabile et al., 2005; Pirro et al., 2006). EMP levels were measured by trained technicians who were blind to the clinical status of the patients as well as to the results.

RNA isolation and preparation of miRNA microarrays

Total RNA was extracted from 250 μl of plasma using TRIzol reagent according to the manufacturer's instructions (Invitrogen) and purified using an RNeasy minikit (Qiagen). To increase RNA recovery, 1 μg of MS2 carrier RNA was added to each plasma sample. The quality of total isolated RNA was determined using the Agilent 2100 Bioanalyzer.

The plasma miRNA expression levels were assessed using Sure Print G3 human 8 × 60 k miRNA microarrays (Agilent Technologies) covering 1,205 human miRNAs (Sanger miRBase release 16). The miRNAs were dephosphorylated and labeled with cyanine 3-cytidine biphosphate including a labeling spike-in solution (Agilent Technologies) to assess labeling efficiency. The samples were hybridized on the arrays with the inclusion of a hybridization spike-in solution (Agilent Technologies) to monitor hybridization efficiency. The arrays were scanned with a G2565CA Microarray Scanner System with SureScan High-Resolution Technology (Agilent Technologies) using Scan Control software. The Feature Extraction 11.5.11 (Agilent Technologies) and GeneSpring 12.6.1. software packages were used for data processing.

Analysis of miRNA microarray expression data and integrative analysis and gene co-expression network

The data from microarrays (deposited at Gene Expression Omnibus under GEO Series accession number GSE101616) were normalized using the robust multi-array average (RMA) method (Irizarry et al., 2003) implemented in the AgiMicroRna Bioconductor package, and the fold changes in circulating miRNAs were determined using the linear model implemented in the limma Bioconductor package (Ritchie et al., 2015). The Benjamini and Hochberg methods were used to adjust p-values for multiple testing and to control the false discovery rate (Benjamini and Hochberg, 1995). The miRNA co-expression network was constructed from the expression profiles of those miRNAs whose expression was significantly associated with PECAM+ EMPs (Spearman p < 0.05). The co-expression network was inferred using graphical Gaussian models (GGMs) implemented in the R package GeneNet (Schäfer and Strimmer, 2005). Briefly, a partial-correlation matrix was estimated by computing the partial correlation between the expression profiles of each miRNA pair. Bayesian posterior edge probability >0.95 (corresponding to a local false discovery rate of < 5%) was used to determine the significance of the resulting pairwise partial correlations. In the resulting co-expression network, which was visualized using Cytoscape software (Cline et al., 2007), the nodes represent the set of miRNAs that were significantly correlated with PECAM+ EMP levels, and the edges link the pairs of miRNAs whose expression was not conditionally independent, defined as the pairwise partial correlation once the common effects of the other miRNAs in the subset were removed (Opgen-Rhein and Strimmer, 2007).

MiRNA quantification by real-time qRT-PCR

TaqMan microRNA assays (Applied Biosystems) were used to quantify the expression of selected miRNAs. Briefly, reverse transcription was performed using the TaqMan MicroRNA Reverse Transcription Kit (Applied Biosystems) and the miRNA-specific oligonucleotides provided with the TaqMan MicroRNA Assay (Applied Biosystems). The final concentration of total RNA used was 2.5 ng/μL. The reaction was performed at 16°C for 30 min, 42°C for 30 min, and 85°C for 5 min. We used 1.33 μL of the obtained cDNAs in a subsequent quantitative qRT-PCR amplification using the TaqMan Universal PCR master mix (Applied Biosystems, Madrid, Spain) and the associated specific probe provided in the TaqMan® MicroRNA Assay Kit (Applied Biosystems). Specific TaqMan probes were used for each gene: let-7d (hsa-let-7d-5p), let-7g (hsa-let-7g-5p), miR-122 (hsa-mir-122-5p), miR-130a (hsa-mir-130a-3p), miR-337 (hsa-mir-337-5p), miR-409 (hsa-mir-409-3p), miR-486 (hsa-mir-486-5p), miR-494 (hsa-mir-494), and miR-718 (hsa-mir-718). The results were normalized to the expression of the U6 small nuclear RNA (U6 snRNA), which was used as an endogenous control. Amplification was performed in a 7900HT thermocycler (Applied Biosystems) at 95°C for 10 min, followed by 40 cycles of 95°C for 15 s and 60°C for 1 min. The fold change in the miRNA level was calculated by the log 2 scale according to the equation 2−ΔΔCt, where ΔCt = Ct miRNA-Ct U6 and ΔΔCt = ΔCt treated samples-ΔCt untreated controls. Negative control reactions, with no RNA, had undetectable quantification cycle values (Cq).

Genomic region overrepresentation analysis

The genetic loci visualization of co-expressed miRNAs across chromosomes was performed using PhenoGram (Wolfe et al., 2013). To determine the significance of the overrepresentation of chromosome regions in the generated co-expressed miRNA set, we applied a hypergeometric test.

Target prediction and functional pathway analysis

Putative miRNA target sites for the miRNA candidates were identified by bioinformatics analysis using miRNA target prediction databases, including DIANA-microT-CDS (Paraskevopoulou et al., 2013), DIANA-TarBase v7.0 (Vlachos et al., 2015a), and TargetScan v6.2 (Agarwal et al., 2015). Functional enrichment analysis of miRNA target genes was conducted using miRPath v3.0 (Vlachos et al., 2015b). Visualization and topological analysis of the network was done using R and Cytoscape software.

Statistical analysis

Because of the right-skewed distribution of the values, the PECAM+ EMP plasma levels underwent a natural logarithmic transformation and were expressed as log-transformed counts per μl (log PECAM+ EMPs/μl). EMP levels are presented as the mean ± SEM. Chi-squared tests, or Fisher exact tests when appropriate, were used to compare the frequencies of the categorical variables. The effects of valve morphology and aortic root dilation were assessed using ANOVA. Tukey's test was utilized for pairwise comparisons. p < 0.05 were considered significant. The statistical analysis was performed using SPSS software, version 21.0 (IBM, Chicago, IL, USA).

Results

Identification of circulating miRNA sets associated with PECAM+ EMP levels

To explore the circulating PECAM+ EMP-associated miRNAs, we first corroborated our previously reported finding that BAV patients have increased levels of PECAM+ EMPs in plasma compared with TAV healthy controls (Supplementary Figure 1A). Furthermore, the results of this analysis also supported the role of PECAM+ EMP circulating levels as a biomarker of aortic dilation in BAV disease, as an increased diameter of the aortic root or ascending aorta was associated with higher levels of PECAM+ EMPs (Supplementary Figure 1B).

We used a microarray-based screening approach to determine the miRNA expression pattern that could be related to endothelial damage. The expression of 1,205 miRNAs were evaluated in plasma samples from BAV patients and TAV healthy controls, and after the miRNA microarray expression data were processed, 277 miRNAs were identified as expressed in at least 5% of the samples. We prioritized the potential miRNA candidates based on the linear relationship between PECAM+ EMP circulating levels and the miRNA expression (Supplementary Figure 1C), and based on this analysis, we found that the expression of 175 miRNAs was significantly associated with the circulating levels of PECAM+ EMPs (p < 0.05).

Construction of a miRNA co-expression network based on the EMP-associated miRNAs

Once we identified the 175 miRNAs whose expression patterns in plasma significantly correlated with the levels of PECAM+ EMPs, we took into consideration that about half of all described miRNAs are co-expressed in clusters of miRNAs. For this reason, we used a Gaussian graphical model approach to map the simultaneous expression of miRNA pairs into a co-expression network.

The resulting miRNA co-expression network is composed of a single connected component formed by 131 miRNAs (FDR < 5%) linked by 391 edges (Figure 1A). We first evaluated the biological implications of these co-expressed miRNAs in terms of intercellular communication in BAV-induced endothelial damage (e.g., whether they are expressed in endothelial cells and whether their expression levels have been previously associated with cardiovascular diseases; Figure 1B). More interestingly, we also focused on miRNAs that were described in previous studies as related to BAV or whose expression might be sensitive to blood flow. We found in the literature that 98 of the 131 co-expressed miRNAs were previously reported as expressed in endothelial cells. Moreover, 88 miRNAs included in the co-expression network had been related to cardiovascular disease. More specifically, modulation of the expression of 4 of them had been associated with BAV, and 9 miRNAs had been described as blood flow sensitive. We also carried out a topological analysis of the network and identified miR-494 as the most important hub within the co-expression network, that is, the miRNA with the largest degree; the highest betweenness, stress, and closeness scores; and the highest radiality coefficient (Figure 1C).

Figure 1

Co-expression network of miRNAs associated with the circulating levels of EMPs. (A) Representation of the miRNA co-expression network. The node color represents the Spearman correlation between the expression of the corresponding miRNA and PECAM+ EMP levels; red nodes correspond to miRNAs that are positively correlated with PECAM+ EMPs, and green nodes correspond to miRNAs that are negatively correlated with PECAM+ EMPs. The edge shape is related to the direction of the partial correlation; the continuous lines represent partial positive correlations, and dotted lines refer to partial negative correlation. (B) Biological implications of the miRNAs included in the co-expression network based on previous knowledge. The border color of the node represent miRNAs involved in cardiovascular diseases or in bicuspid aortic valve. The hexagon node shape corresponds to flow-sensitive miRNAs. The node color represents miRNA expression in endothelial or hematological cells. (C) Topological analysis of the miRNA co-expression network. miR-494 was identified as the master switch based on its role as the most important hub in the co-expression network. The color and the size of the node refer the topological importance of the miRNA within the network. Node positions are conserved between networks.

To further validate the results obtained in the microarray analysis and the EMP-miRNA correlations, we determined the expression of 7 of miRNAs included in the miRNA co-expression network (let-7d, let-7g, miR-130a, miR-337, miR-409, miR-494, and miR-718; Figure 2) by RT-qPCR in a new cohort of TAV individuals and BAV patients (n = 36). In this way, by correlating the RT-qPCR-determined expression data for each of the selected miRNAs with the circulating levels of PECAM+ EMPs, we validated and reaffirmed the sensitivity and power of our bioinformatics approach for the identification and selection of PECAM+ EMP-associated miRNA candidates.

Figure 2

(A-G) The microarray results and the integrative analysis were further validated by determining the expression of 7 miRNAs included in the co-expression network. The expression of these miRNAs was determined by RT-qPCR and correlated with circulating PECAM+ EMP levels.

Genomic region enrichment analysis for the miRNA co-expression network

Because co-expressed miRNA pairs tend to reside in close genomic proximity, we performed a positional gene enrichment analysis to identify chromosome regions that were overrepresented among the PECAM+ EMPs-associated and co-expressed miRNAs.

We constructed chromosome ideograms to map the genomic locus of each of the miRNAs included in the co-expression network (Figure 3A). We found that 21 of the 131 PECAM+ EMP-associated co-expressed miRNAs were located on chromosome 14 and that 19 of them were located in the same chromosome region, the 14q32 locus. We confirmed the statistical significance of this 14q32 genomic region overrepresentation using a hypergeometric test (p = 1.90 × 10−7). This genomic location consists of the miRNA clusters A and B and is also known as the Dlk1-Dio3 miRNA cluster. Specifically, we identified 4 of the 8 miRNAs located within miRNA cluster A (miR-127-3p, miR-337-3p, miR-431-5p, and miR-493-5p) and 15 of the 42 miRNAs located within miRNA cluster B (miR-154-5p, miR-376a-3p, miR-376c-3p, miR-377-3p, miR-379-5p, miR-381-3p, miR-382-5p, miR-409-3p, miR-410, miR-487b, miR-494, miR-495-3p, miR-543, miR-654-3p, and miR-654-5p; Figures 3B–C).

Figure 3

Genomic region enrichment analysis of the miRNA co-expression network. (A) The locus of each of the miRNAs included in the co-expression network was mapped in the chromosome ideograms to visualize the genomic loci of co-expressed miRNAs. (B) Detailed genomic localization of the co-expressed miRNAs located on chromosome 14. The genomic region overrepresentation was corroborated by the results of a hypergeometric test (p = 1.9 × 10−7). (C) The miRNAs located in the 14q32 locus are highlighted within the miRNA co-expression network.

Functional analysis of the 14q32 miRNA cluster

To gain deeper insight into the biological significance of these findings, we determined the biological significance of the regulation of the 14q32 miRNA cluster based on the functional implications of the genes targeted by the selected miRNAs. We carried out a functional enrichment analysis of the set of genes putatively targeted by the 19 EMP-related miRNAs located in the 14q32 cluster. Interestingly, this functional analysis revealed that the genes targeted by the 14q32 co-expressed miRNAs were categorized with the “cellular nitrogen compound metabolic process” (p = 2.34 × 10−145), “immune system process” (p = 2.57 × 10−6), and “extracellular matrix organization” (p = 8.14 × 10−5) GO terms as well as the “TGF-β signaling pathway” (p = 2.59 × 10−8) KEGG term, among others (Figure 4).

Figure 4

Functional enrichment analysis of the genes targeted by the 14q32 miRNAs represented in the miRNA co-expression analysis. (A) A network representing the genes targeted by each of the miRNAs included in the miRNA co-expression network and located within the 14q32 locus and their association with the highlighted biological process gene ontologies. (B) The significance of the functional enrichment analysis of the genes targeted by the 14q32-located miRNAs.

Discussion

Using a bioinformatics approach, our study reported that the endothelial damage observed in BAV disease results in the differential regulation of a post-transcriptional regulatory miRNA network associated with the increased release of endothelial-derived microparticles. By determining the circulating miRNAs associated with PECAM+ EMP levels, we inferred a miRNA co-expression network. Furthermore, based on the results of genomic enrichment analysis, we focused on a cluster of highly co-expressed miRNAs located at the 14q32 locus on chromosome 14 that could be acting as molecular effectors in BAV-related pathophysiological processes, including endothelial damage. We showed that the miRNAs contained within the 14q32 miRNA cluster could mediate crucial biological processes in BAV disease, such as nitric oxide biosynthesis, immune activation, reorganization of the extracellular matrix and TGF-β signaling.

Bicuspid morphology of the aortic valve is associated with haemodynamic abnormalities that result in increased wall shear stress on the endothelial layer, which contributes to dilation of the ascending aorta (Braverman et al., 2005). At the cellular level, this disturbed blood flow is associated with rearrangement of the extracellular matrix and with endothelial damage and dysfunction (Davignon and Ganz, 2004). Furthermore, in our previous studies, we demonstrated that BAV and dilation of the aorta are associated with endothelial-mediated release of microparticles (Alegret et al., 2016). Thereby, when designing this study, we hypothesized that the characteristics of blood flow through the ascending aorta in patients with BAV disease are related to aortic endothelial damage and play key a role in EMP generation.

It is important to consider that on the one hand, circulating miRNAs take very stable forms, mainly packaged in transport particles; circulating microparticles have been reported as the major carriers of miRNAs in the blood (Diehl et al., 2012). On the other hand, although most published studies have focused mainly on differences in the expression of single miRNAs instead of studying miRNA co-expression networks, miRNAs are significantly enriched in clusters in discrete genomic regions (Lau, 2001; Kim and Nam, 2006), and miRNAs in the same cluster might be transcribed in a polycistronic manner (Baskerville, 2005; Wang et al., 2016) and likely regulate functionally related genes (Kim et al., 2009; Wang et al., 2011). We considered the determination and analysis of the miRNA co-expression network resulting from the integration of PECAM+ EMP levels with the expression of circulating miRNAs as an advantageous strategy to unravel the molecular mechanisms underlying the pathophysiological processes in BAV disease for the following reasons: changes in haemodynamic forces in the vascular system might alter the expression of miRNAs in endothelial cells (Marin et al., 2013); BAV promotes endothelial dysfunction, which results in the release of PECAM+ EMPs in plasma; microparticles are the main carriers of miRNAs in plasma; and miRNA-coding genes are enriched in discrete genomic regions.

Using a bioinformatics approach, we identified 175 miRNAs that were significantly associated with PECAM+ EMP levels. Of these miRNAs, 131 exhibited significant pairwise partial correlations and were inferred into a co-expression network to analyse miRNA co-expression patterns. Co-expression network analysis is considered a powerful method to extract strong regulatory associations that could be responsible for modulating transcriptional networks underlying biological processes. Accordingly, we first analyzed previous knowledge regarding the biological implications of the inferred network and found that 75% of the co-expressed miRNAs had been previously described as miRNAs expressed in endothelial cells and that 67% of them had been associated with cardiovascular diseases. We also analyzed the topology of the co-expression network and found that miR-494 was the most important hub of the network because it was the most highly connected miRNA. miR-494 was also the most influential miRNA within the network based on centrality parameters, which showed that miR-494 was most often found on the shortest path between two other miRNAs. These results indicate that miR-494 might have a large regulatory impact on the biological functions that could be modulated by the co-expression network.

We further analyzed the miRNA co-expression network in terms of enrichment for specific genomic locations, and we found that the 14q32 locus on chromosome 14 was significantly overrepresented as a highly co-expressed miRNA cluster for the PECAM+ EMP-associated and co-expressed miRNAs in BAV. This genomic region contains the largest miRNA cluster found in the human genome (Benetatos et al., 2013). Specifically, we found that 19 of the 131 miRNAs inferred in the co-expression network, including miR-494, are located within this chromosome region. This 14q32 miRNA cluster is also known as the imprinted DLK1-MEG3 genomic region. This region contains the protein-coding genes DLK1, RTL1, and DIO3, which are expressed from the paternally inherited chromosome. The region also contains multiple long and short non-protein coding RNAs, including MEG3, MEG8, the miRNA cluster and small nucleolar RNA (snoRNA) genes, which are expressed from the maternally inherited chromosome. The imprinting of a genomic region refers to biased expression of the genes contained in either the paternally or maternally inherited chromosome instead of the more common biallelic expression. Although, miRNAs located within the 14q32 region have been proposed as candidates for various diseases, including cancer, psychiatric illness, alcoholism, and non-alcoholic fatty liver disease, to date, there is no published data relating the 14q32 miRNA cluster with cardiovascular diseases (Benetatos et al., 2013; Okamoto et al., 2016).

Based on the functional enrichment of the genes targeted by the PECAM+ EMP-associated and co-expressed miRNAs located within the 14q32 locus, we could determine the biological significance of the regulation of this miRNA cluster. We found that the genes targeted by the 14q32 co-expressed miRNAs are strongly involved in the regulation of genes in categories such as “cellular nitrogen compound metabolic process,” “immune system process,” “extracellular matrix organization,” and “TGF-β signaling pathway.” The “cellular nitrogen compound metabolic process” gene ontology term includes the “nitric oxide biosynthetic pathway” subcategory. Nitric oxide (NO) is a pivotal endothelium-derived substance that plays a crucial role in the homeostasis of the cardiovascular system by modulating endothelium-dependent vasodilation; in fact, impairment of NO production or activity has been proposed as a major mechanism of endothelial dysfunction (Davignon and Ganz, 2004). Moreover, modulation of the expression and activity of endothelial nitric oxide synthase (eNOS), the specific enzyme that produces NO in endothelial cells, by fluid shear stress and the implication of this modulation in the development of BAV are well established (Ranjan et al., 1995; Aicher et al., 2007; Vion et al., 2013). Thus, in the context of BAV, the valvulopathy and the associated abnormal haemodynamics are related to the altered distribution of aortic wall shear stress, which promotes regional disruption of the eNOS pathway in a manner that is primarily mediated by the differential expression and activity of eNOS. Endothelial cells release EMPs in response to cellular stress and cell activation (Dignat-George and Boulanger, 2011). Consistent with this fact and the role of PECAM+ EMP-related and co-expressed miRNAs in the regulation of endothelial dysfunction, we found that among the genes targeted by the 14q32-located miRNAs, those with functional implications related to modulation of immune system processes and to extracellular matrix organization were significantly enriched. In fact, interactions between inflammatory activation and endothelial dysfunction have been previously described in patients with BAV (Ali et al., 2014). In addition, we also found that the 14q32-located miRNAs present in our co-expression network targeted genes that were related to the TGF-β signaling pathway. The TGF-β family is composed of several cytokines with diverse functions, including the regulation of tissue repair and fibrosis, extracellular matrix remodeling, and inflammation as well as cell proliferation and migration (Bobik, 2006). Previous studies have implicated the impairment of TGF-β signaling in BAV pathology and the development of BAV-associated aortopathy (Forte et al., 2016). Additionally, functional interactions between TGF-β and NO have been demonstrated in endothelial cells (Saura et al., 2005). All these data suggest that the 14q32 miRNA cluster may play a pivotal role in the induction of BAV-associated endothelial damage.

Furthermore, in addition to studying the post-transcriptional implications of the regulatory functions of the 14q32 miRNA cluster, we thought it would be interesting to understand the regulatory mechanisms that could affect the activation or repression of the expression of this miRNA cluster. The imprinted status of the maternally expressed RNAs of the DLK1-MEG3 locus, including the 14q32 miRNA cluster, is regulated by epigenetic mechanisms, specifically by the methylation of two differentially methylated regions (DMRs) located upstream of the transcription activation site of MEG3 (Murphy et al., 2003; Kameswaran et al., 2014). The hypermethylation of either of these two DMRs has been associated with decreased expression of the maternal transcript (Kagami et al., 2010). In various models, such as in human type 2 diabetic islets, repression of the 14q32 miRNA cluster is strongly correlated with hypermethylation of the MEG3-DMR, and modifications at this region increase susceptibility to disease (Kameswaran et al., 2014). Interestingly, Boon et al. (2016) proposed that aging induces the expression of MEG3 and that MEG3-mediated changes in the epigenetic regulation of gene expression contributes to aging-related endothelial dysfunction. Moreover, Gordon et al. (2010) suggested that MEG3 may regulate the expression of VEGF; they showed that loss of MEG3 leads to up-regulation of the expression of genes in the VEGF and Notch signaling pathways in mouse brains. Furthermore, DNA methylation mechanisms are responsive to disturbed flow (Dunn et al., 2014); in fact, endothelial gene expression can be regulated by flow-dependent epigenetic mechanisms (Jiang et al., 2015). Therefore, the disturbed flow associated with BAV may promote not only the secretion of EMPs but also the MEG3-mediated epigenetic regulation of the 14q32 miRNA cluster and thus may play a key role in the regulation of endothelial damage in BAV disease. On the other hand, DLK1 is the best-studied non-canonical Notch ligand and acts as an inhibitor of Notch signaling in vitro (Baladrón et al., 2005). In this manner, regulation of the expression of genes located in the 14q32 locus could also modulate the crosstalk between TGF-β and the VEGF and Notch signaling pathways, which have been previously described to act co-ordinately in space and time in the regulation of vascular morphogenesis (Holderfield and Hughes, 2008).

Limitations

We studied the miRNAs associated with the secretion of PECAM+ EMPs in the context of BAV-induced endothelial damage, but we did not determine if this pattern of expression of PECAM+ EMPs-associated miRNAs is also deregulated when the endothelial damage is caused by other conditions.

In summary, we propose that the changes in hemodynamic forces in the vascular system that are caused by the bicuspid morphology of the aortic valve might result in the expression of a specific pattern of PECAM+ EMP-associated miRNAs that regulates a potent post-transcriptional network that might be involved in the regulation of endothelial damage. We reported that the highly co-expressed and PECAM+ EMP-associated miRNAs clustered within the 14q32 locus might modulate a functional regulatory miRNA co-expression network that targets functionally related genes, which could be the molecular effectors linking the impairment of various signaling pathways (VEGFA, TGF-β, NOTCH) involved in the pathophysiology of BAV disease. Thus, we propose that the 14q32 miRNAs may act as master switches for endothelial damage and that inhibition of either individual or a combination of 14q32 miRNAs might offer a new therapeutic approach for BAV disease.

Author contributions

NM: Performed the experiments, analyzed the data, and wrote the manuscript; RB and GA: Performed the experiments; MF: Contributed in sample collection and writing the manuscript; JA: Conceived and designed the study, analyzed the data and revised the manuscript. All authors read and approved the final manuscript.

Conflict of interest statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.