Systems Biology in Chronic Heart Failure—Identification of Potential miRNA Regulators

Heart failure (HF) is a complex disease entity with high clinical impact, poorly understood pathophysiology and scantly known miRNA-mediated epigenetic regulation. We have analysed miRNA patterns in patients with chronic HF (cHF) and a sex- and age-matched reference group and pursued an in silico system biology analysis to discern pathways involved in cHF pathophysiology. Twenty-eight miRNAs were identified in cHF that were up-regulated in the reference group, and eight of them were validated by RT-qPCR. In silico analysis of predicted targets by STRING protein-protein interaction networks revealed eight cluster networks (involving seven of the identified miRNAs) enriched in pathways related to cell cycle, Ras, chemokine, PI3K-AKT and TGF-β signaling. By ROC curve analysis, combined probabilities of these seven miRNAs (let-7a-5p, miR-107, miR-125a-5p, miR-139-5p, miR-150-5p, miR-30b-5p and miR-342-3p; clusters 1–4 [C:1–4]), discriminated between HF with preserved ejection fraction (HFpEF) and HF with reduced ejection fraction (HFrEF), and ischaemic and non-ischaemic aetiology. A combination of miR-107, miR-139-5p and miR-150-5p, involved in clusters 5 and 7 (C:5+7), discriminated HFpEF from HFrEF. Pathway enrichment analysis of miRNAs present in C:1–4 (let-7a-5p, miR-125a-5p, miR-30b-5p and miR-342-3p) revealed pathways related to HF pathogenesis. In conclusion, we have identified a differential signature of down-regulated miRNAs in the plasma of HF patients and propose novel cellular mechanisms involved in cHF pathogenesis.


Introduction
Heart failure (HF), affecting 1-2% of the world's population over the age of 65, is one of the leading causes of mortality and morbidity [1,2]. It is a complex disease where the heart is unable to meet the metabolic demands of the body [3,4]. Chronic HF (cHF) has several underlying aetiologies (ischaemic disease, dilated cardiomyopathy, hypertrophy cardiomyopathy, hypertension, heart valve disease, etc.) that produce alterations in heart structure (left ventricle hypertrophy or dilation, cardiac fibrosis, heart remodelling) and function (abnormal cardiac output, diastolic dysfunction), leading to common symptomatology [3]. The underlying pathological mechanisms that conduce to cHF are not completely understood, but several biological pathways are involved in the progression of the disease [5].
Plasma sample analysis allows the easy acquisition of biologic material, such as mi-croRNAs (miRNAs) that may assist in disease characterisation and patient risk stratification.

Differential miRNA Signature in cHF
We comparatively analysed the plasma levels of 377 miRNAs in cHF patients and reference subjects (RS). Of these 377 miRNAs, 167 were consistently detected in both groups. Differential expression analyses between cHF and RS showed that 28 of these miRNAs presented statistically significant differences and a fold-change (log 2 ) higher than ±0.38 between groups. Of these, the top 15 miRNAs presenting an area under the curve (AUC) ≥0.7 after C-statistics analysis and a fold-change >1.7 (0.75 in log 2 ) were selected for further analysis. Three of the top-fifteen miRNAs were of the same family (let-7a-5p, let-7f-5p and let-7g-5p) and have been reported to be involved in the regulation of common pathways [17]. Therefore, only miRNA let-7a-5p, with the best AUC among them, was further analysed. In addition, and considering the probability of a low detection level, we excluded those miRNAs that according to miRBase presented a number of reads in the lowest 25% range ( Figure 1 and Supplemental Table S1).

In Silico System Biology Analysis and Target Prediction
Gene-targets of the differential miRNA signature (577 genes) were obtained with the combination of four different miRNA-target prediction algorithms (TargetScan, miRmap, miRWalk and miRDB [18][19][20][21][22]). All validated miRNAs were included in target data mining (let-7a-5p, miR-107, miR-125a-5p, miR-139-5p, miR-150-5p, miR-30b-5p and miR-342-3p), except miR-335-5p, given that, for this miRNA, only one target was obtained after combining predictions of all four algorithms used. The resulting targets were analysed using the ShinyGO web-based software [23] for pathway enrichment analysis (Table 1). Gene ontology (GO) annotations for biological processes of predicted targets, sorted by average and rank (FDR and fold), revealed significant enrichment in pathways related to coronary vasculature morphogenesis and development (GO:0060977 and GO:0060976, respectively; FDR < 0.001, both), TGF-β signalling regulation (GO:0030511 and GO:1903846; FDR < 0.001, both) and with ventricular and cardiac septum development (GO:0003281 and GO:0003279, respectively; FDR < 0.001, both). Additionally, a protein-protein interaction (PPI) network, including all 577 predicted targets, was built according to data from the STRING database [24] using Cytoscape [25]. The resulting network revealed a highly connected nodule involving 273 targets from all included miRNAs (Figure 3a). We then applied a community cluster strategy to identify PPIenriched nodules [26] and ran a functional overrepresentation analysis upon pathways from the KEGG database [27,28]. Clusters with fewer than five nodes were excluded from the analysis. As a result, eight clusters were identified from the main connected component of the network (Figure 3b). Noteworthy, four clusters included targets from all query miRNAs (clusters 1-4) and were the ones selected for more specific pathway analysis. As depicted in Supplemental Table S3, the identified clusters aggregate functionally related genes. Highly significant enrichment for focal adhesions (hsa04510; FDR < 0.001), chemokine, Ras and PI3K-AKT signalling pathways (hsa04062, hsa04014 and hsa04151, respectively; FDR < 0.001, all) was found for cluster 1. Also, strong enrichments in cell cycle (hsa04110; FDR < 0.001), cellular senescence (hsa04218; FDR < 0.001) and TGF-β signalling (hsa04350; FDR < 0.001) were detected for cluster 2. For clusters 3 and 4, only marginal enrichments were found in genes involved in RNA transport (hsa03013; FDR < 0.001), Wnt and Hippo signalling pathways (hsa04310 and hsa04390, respectively; FDR < 0.001). Additionally, a protein-protein interaction (PPI) network, including all 577 predicted targets, was built according to data from the STRING database [24] using Cytoscape [25]. The resulting network revealed a highly connected nodule involving 273 targets from all included miRNAs (Figure 3a). We then applied a community cluster strategy to identify PPI-enriched nodules [26] and ran a functional overrepresentation analysis upon pathways from the KEGG database [27,28]. Clusters with fewer than five nodes were excluded from the analysis. As a result, eight clusters were identified from the main connected component of the network (Figure 3b). Noteworthy, four clusters included targets from all query miRNAs (clusters 1-4) and were the ones selected for more specific pathway analysis. As depicted in Supplemental Table S3, the identified clusters aggregate functionally related genes. Highly significant enrichment for focal adhesions (hsa04510; FDR < 0.001), chemokine, Ras and PI3K-AKT signalling pathways (hsa04062, hsa04014 and hsa04151, respectively; FDR < 0.001, all) was found for cluster 1. Also, strong enrichments in cell cycle (hsa04110; FDR < 0.001), cellular senescence (hsa04218; FDR < 0.001) and TGF-β signalling (hsa04350; FDR < 0.001) were detected for cluster 2. For clusters 3 and 4, only marginal enrichments were found in genes involved in RNA transport (hsa03013; FDR < 0.001), Wnt and Hippo signalling pathways (hsa04310 and hsa04390, respectively; FDR < 0.001). Only interactions with a confidence score higher than 0.8 were considered. Node colours represent targeting miRNA. Grey nodes correspond to STRING-added nodes. (b) Community clusters identified within the main connected component of the PPI network. Clusters with less than 5 genes were excluded. (c) Table with number of target genes regulated by each miRNA, per cluster. cHF: chronic heart failure; PPI: protein-protein interaction.

Association of Clusters with cHF Subtypes
Considering that the networks obtained in the PPI analysis associated with specific pathways, we analysed if clustered miRNAs provided new information to differentiate
To study if the miRNAs implicated in the clusters with discriminatory potential associated with specific pathways, we performed two other PPI networks. Specifically, we compared if the miRNAs included in C:5+7 (miR-107, miR-139-5p and miR-150-5p) and the miRNAs not included in these clusters but present in C:1-4 (let-7a-5p, miR-125a-5p, miR-30b-5p and miR-342-3p; from now on "Differential miRNAs [DiffmiR]") were enriched in different biological pathways. Further, we analysed if those divergent pathways could be associated with the discriminatory differences observed between C:1-4 and C:5+7.
Thus, two PPI networks were generated in Cytoscape with data from STRING [24]. Networks for C:5+7 and DiffmiR were created with 408 and 189 predicted target genes, respectively. The new resulting networks included 204 predicted targets for C:5+7 and 63 for DiffmiR (Figures 4a and 5a). To identify PPI-enriched nodules, we applied a community cluster strategy [26], as described above, and pursued a functional overrepresentation analysis upon pathways from the KEGG database [27,28]. Only clusters with more than five nodes were considered for analysis. To study if the miRNAs implicated in the clusters with discriminatory potential associated with specific pathways, we performed two other PPI networks. Specifically, we compared if the miRNAs included in C:5+7 (miR-107, miR-139-5p and miR-150-5p) and the miRNAs not included in these clusters but present in C:1-4 (let-7a-5p, miR-125a-5p, miR-30b-5p and miR-342-3p; from now on "Differential miRNAs [DiffmiR]") were enriched in different biological pathways. Further, we analysed if those divergent pathways could be associated with the discriminatory differences observed between C:1-4 and C:5+7.

miR-107 and miR-139 Were Overrepresented in All Clusters
When the weight of each miRNA in each cluster of the main PPI network was considered (Figure 3a), we observed that miR-107 and miR-139-5p were overrepresented, with 81 (26.4%) and 95 (31.25%) of all clustered protein-coding genes targeted by these miRNAs, respectively (Figure 3c). Therefore, miR-107 and miR-139-5p were analysed individually to better understand their specific contribution to cHF pathophysiology.
In silico analysis querying "heart failure" to STRING Pubmed in Cytoscape [25] resulted in a repository-based list of 837 genes related to HF (repository-based gene database [RBGD]). The 81 and 95 targets of miR-107 and miR-139-5p, respectively, were matched individually to the RBGD list (Supplemental Figure S1). Using this approach, we found three protein-coding genes regulated by miR-107 and nine regulated by miR-139-5p that were also identified as genes with a role in HF, according to the repository-derived data included in the RBGD list (Table 3). Table 3. Clustered miRNA targets detected in the STRING-generated heart failure database.

Association of Individual Plasma miRNAs with Aetiology and with Left Ventricular Ejection Fraction in cHF Patients
A study of individual miRNAs showed that miR-139-5p presented differential expression patterns between cHF patients with and without underlying ischaemic aetiology (p = 0.022), but not for left ventricular ejection fraction (LVEF) levels (Figure 8 and Supplemental Tables S6 and S7). miRNA plasma levels did not show statistically significant changes between patients with different degrees of severity according to the New York Heart Association (NYHA) classification. These 12 common protein-coding genes and their RBGD-based neighbouring interacting genes (101 for miR-107 and 176 for miR-139-5p) were selected and analysed with ShinyGO [23] to determine their putative pathways. Using the KEGG database [27,28], we observed that both miRNAs participated in the pathways PI3K-AKT signalling (hsa04151; FDR < 0.001; both miRNAs) and MAPK signalling (hsa04010; FDR < 0.001; both miRNAs) (Figure 7a,b).  In blue: miR-139-5p target genes and their neighbouring genes; in yellow: miR-107 target genes and their neighbouring genes; in green: miR-139-5p and miR-107 common target genes and their neighbouring genes; in white: additional nodes for pathway completion. Unframed gene names correspond to targeted genes included under a broader category such as "growth factor", "RTK" or "cytokine". Adapted from pathways "MAPK signalling pathway (hsa04010)" and "PI3K-AKT signalling pathway (hsa04151)" generated in the KEGG database.

Association of Individual Plasma miRNAs with Aetiology and with Left Ventricular Ejection Fraction in cHF Patients
A study of individual miRNAs showed that miR-139-5p presented differential expression patterns between cHF patients with and without underlying ischaemic aetiology (p = 0.022), but not for left ventricular ejection fraction (LVEF) levels (Figure 8 and In blue: miR-139-5p target genes and their neighbouring genes; in yellow: miR-107 target genes and their neighbouring genes; in green: miR-139-5p and miR-107 common target genes and their neighbouring genes; in white: additional nodes for pathway completion. Unframed gene names correspond to targeted genes included under a broader category such as "growth factor", "RTK" or "cytokine". Adapted from pathways "MAPK signalling pathway (hsa04010)" and "PI3K-AKT signalling pathway (hsa04151)" generated in the KEGG database.

Discussion
Heart failure is a highly frequent disease with a complex pathophysiology [3]. Circulating miRNAs are emerging as tools to facilitate disease diagnosis, prognosis, risk stratification and management [6][7][8][9][10]. In addition, miRNAs have a role in gene expression regulation and are, therefore, involved in disease pathogenesis [12,13]. In this study, we have used a non-targeted approach to miRNA identification and investigated plasma samples of cHF patients. We have studied the putative role of miRNAs in cHF pathogenesis through gene regulation by pursuing an in silico enrichment pathway analysis. Additionally, we analysed if these miRNAs were individually related to pathways associated with cHF pathogenesis.

Discussion
Heart failure is a highly frequent disease with a complex pathophysiology [3]. Circulating miRNAs are emerging as tools to facilitate disease diagnosis, prognosis, risk stratification and management [6][7][8][9][10]. In addition, miRNAs have a role in gene expression regulation and are, therefore, involved in disease pathogenesis [12,13]. In this study, we have used a non-targeted approach to miRNA identification and investigated plasma samples of cHF patients. We have studied the putative role of miRNAs in cHF pathogenesis through gene regulation by pursuing an in silico enrichment pathway analysis. Additionally, we analysed if these miRNAs were individually related to pathways associated with cHF pathogenesis.
In contrast to our findings at the systemic level, let-7a-5p was found up-regulated in plasma from the coronary sinus of cHF patients when compared to controls [29]; although differences between both studies could relate to the plasmatic sample origin. Our study could reflect systemic changes in cHF, whereas those miRNA patterns in the coronary sinus could be more related to pathological alterations in the heart [29]. In agreement with our results, Ellis and colleagues [30] detected down-regulated levels of miR-150-5p in the plasma of cHF in comparison to controls [30]. Similarly, Scrutinio et al. [31] reported miR-150-5p to be down-regulated in cHF patients with severe symptomatology in comparison to cHF with mild symptomatology and to controls, being miR-150-5p associated with an increased risk of death, heart transplantation or left ventricular assist devices implantation [31]. Furthermore, the serum levels of miR-150-5p have been reported to be significantly decreased in acute myocardial infarction (AMI) patients compared to controls, with AMI patients who developed cHF presenting the lowest levels [32]. The same study suggested that miR-150-5p was able to predict post-AMI cHF when combined with brain natriuretic peptide (BNP) [32] and that miR-150-5p was independently associated with post-AMI cHF through multivariate analysis including age, BNP and LVEF [32]. Although we did find a negative correlation of this miRNA with LVEF, we did not observe an association of miR-150-5p with ischaemic disease. This could be due to the sample type (plasma instead of serum) [33], or to the fact that we compared those patients with cHF due to an AMI with those patients that developed cHF due to other reasons, such as heart valve disease or cardiac hypertrophy, which could also affect miR-150-5p levels.
In contrast, we observed that miRNA miR-139-5p was up-regulated in the plasma of those patients with underlying ischaemic disease in comparison with those with underlying non-ischaemic disease. Divergent results were found by Voellenkle et al., who observed similar levels of this miRNA (from peripheral blood mononuclear cells [PBMCs]) between cHF with ischaemic and non-ischaemic backgrounds [34]. Nevertheless, in agreement with our findings, the authors described down-regulated levels of miR-139-5p in cHF [34]. The difference in ischaemia-related miR-139-5p response between blood cells and plasma samples suggests that a tissue origin other than blood accounts for circulating miR-139-5p in patients who suffered an AMI. In this respect, Ming et al. have reported decreased levels of miR-139-5p in left ventricle tissue of patients with hypertrophic cardiomyopathy and demonstrated an anti-hypertrophic effect mediated by this miRNA through the regulation of IGF-R1 and c-Jun [35]. It is known that IFG-R1 is a component of both MAPK and PI3K-Akt signalling pathways [36]. Interestingly, the assessment of miRNAs in-cluster individual weight in our study evidenced the presence of miR-139-5p and miR-107 in all clusters. Further, the in silico analysis revealed that both miRNAs were implicated in the MAPK and PI3K-AKT signalling pathways. Voellenkle et al. also detected miR-107 down-regulation in cHF compared to their controls [34]. Some reports have associated this miRNA with metabolism and cardiac function [37], as well as angiogenesis [38] or hypoxia [39]. Further, some reports indicate that miR-107 regulates gene expression through the PI3K-AKT and MAPK pathways [40,41], in correlation with our in silico results. Moreover, miR-30b-5p was found to be deregulated in HF compared to both, healthy subjects and patients with shortness of breath [30].
To our knowledge, this is the first time miR-335-5p has been identified in patients recruited in a heart failure clinic, that is, patients with clinical evidence and diagnosis of chronic heart failure. Some studies including ischaemic patients detected a down-regulation of this miRNA in acute ischaemia [6] or during a three-month post-myocardial infarction period [42]. The latter study suggested miR-335-5p together with five other miRNAs as potential markers of deleterious left ventricular remodelling, but it did not provide evidence of changes in patients with established chronic heart failure. Studies in animal models have also related this miRNA with atherosclerosis [43] or ventricular remodelling [42,44,45], with apparently controversial findings. Thus, studies in mice have associated miR-335-5p overexpression with a reduction of atherosclerotic plaque vulnerability and angiogenesis promotion [46]. In contrast, this miRNA was found up-regulated in rat and mouse models of pulmonary hypertension, where the induced right ventricular remodelling was linked to increased miR-335-5p expression and its blockade alleviated remodelling [45]. Interestingly, contrasting results were observed in a model of sepsis-induced myocardial injury [44]. Long and colleagues observed that up-regulated levels of this miRNA ameliorated myocardial injury by reducing apoptosis and Ca 2+ overload [44]. Nevertheless, we need to take into consideration that miRNA patterns in animal models may not reflect a human profile [47] and that all those studies reflect an acute injury rather than a chronic disease. Similarly, little is known about miR-342-3p. To this respect, Ellis et al. observed down-regulated levels of miR-342-3p in cHF in comparison to healthy subjects or patients with shortness of breath [30], but no association of this miRNA with background aetiology or ejection fraction in cHF has been reported.
The in silico analysis we performed with the validated miRNAs (main PPI network) revealed an enrichment of several pathways, including cell cycle, cellular senescence, and Ras, chemokine, PI3K-AKT and TGF-β signalling pathways. Multiple pathophysiological mechanisms are involved in HF. Among them, cHF patients are characterised by presenting an unresolving inflammatory state, endothelial dysfunction and cardiac remodelling (including cardiac fibrosis, hypertrophy and apoptosis). In agreement with our results, Ras/MAPK and PI3K-AKT signalling have been previously related to cardiac hypertrophy and HF [36,48]. Similarly, the TGF-β and PI3K-AKT pathways have been related to cardiac fibrosis [49][50][51]. Further, HF patients are characterised for presenting an unresolved inflammatory state, where chemokine levels are increased [52] and innate immunity cell increase their release of extracellular vesicles [53].
It is widely accepted that in HF there is a change in cardiac metabolism. In the healthy heart, cardiomyocytes produce energy mainly through fatty acid oxidation. However, in the failing heart, energy production swifts to a glucose-based metabolism (glycolysis and anaplerosis) and other forms of metabolism, such as the use of lactate, branched-chain amino acids and ketone bodies [54]. Although by in silico analysis we evidenced that miRNA-target genes in the DiffmiR cluster were associated with metabolism (glucagon signalling pathway, lipolysis regulation and carbohydrate digestion and absorption), a specific implication of these metabolic pathways in metabolic remodelling in cHF is, as far as we know, still undetermined. However, alterations in the glucagon signalling pathway (compromising cAMP signalling in cardiomyocytes) and increased levels of adipocyte lipolysis, have been detected in HF [55,56]. Further, suppression of HF development was observed in a mice model fed with a strict carbohydrate-restricted diet [57].

Study Design
Patients with cHF treated by guideline-directed medical therapy (GDMT) were prospectively recruited in the outpatient HF unit of the Hospital de la Santa Creu i Sant Pau (Barcelona, Spain) from September 2016 to July 2018 (n = 58). Patients with LVEF between 40% and 50% (mildly reduced LVEF; HFmrEF), past history of cancer, inflammatory disorders, sepsis or infection were excluded from the study. Pregnant women were also excluded. A group of reference subjects of age-and sex-matched subjects without cHF (n = 26) was also analysed.
A quantitative Low-density TaqMan™ Advanced miRNA Human A Card was run on 28 patients, 14 patients with preserved LVEF (HFpEF; LVFE ≥ 50%), and 14 with reduced LVEF (HFrEF; LVFE < 40%) ( Table 4). The reference-subject group comprised 16 sex-and age-matched subjects without cHF. Baseline demographic and biochemical data, classical CV risk factors, background medication and clinical outcomes of cHF patients and reference subjects are included in Tables 4 and 5. A schematic diagram of the study design is shown in Figure 9.  Figure 9. Schematic diagram of the study design. cHF: chronic heart failure; Reference Sub n = 26 Figure 9. Schematic diagram of the study design. cHF: chronic heart failure; RS: reference subjects.
The identified differential miRNAs were validated using the quantitative Low-density Custom TaqMan ® Array Advanced MicroRNA Cards in cHF patients (n = 58). Patients consistently expressing miRNA species in plasma were further investigated (n = 46). Samples of the reference subjects were similarly processed (n = 26). HFrEF affected 52.2% of patients, and 30.4% of cHF had ischaemic disease as the underlying cause of cHF. Clinical and biochemical characteristics, as well as clinical outcomes of these groups, are shown in Tables 6 and 7.  AMI: acute myocardial infarction; cHF: chronic heart failure; HTx: heart transplantation; NT-proBNP: N-terminal prohormone of brain natriuretic peptide; SD: standard deviation. 1. Includes patients that suffered a stroke, an AMI, a cardiovascular death (mainly due to cHF), or were admitted to the emergency department. It does not include patients that underwent HTx or patients that were rehospitalised. 2. Includes patients that died due to a septic shock, a haemorrhage, or a non-successful HTx.

Blood Sampling
After 10-14 h of fasting, venous blood was withdrawn without tourniquet into 3.8% sodium citrate tubes. All samples were processed identically and within the first 2 h. Firstly, blood was centrifuged at 1560× g for 20 min at room temperature to eliminate blood cells and generate platelet-poor plasma (PPP). Secondly, to obtain platelet-free plasma (PFP), samples were centrifuged at 1500× g for 10 min at room temperature. Finally, aliquots were stored at −80 • C until miRNA analyses were performed.

miRNA Extraction, cDNA Synthesis and Analysis
The isolation of miRNA from PFP was performed using the miRNeasy Serum/Plasma Advanced kit (217201; Qiagen; Hilden, Germany) following the provider's instructions. During miRNA extraction, and according to the manufacturer's directions, all samples were spiked-in with 25 fmol of the exogenous miRNA cel-miR-39-3p.
miRNAs were transcribed and pre-amplified using the TaqMan™ Advanced miRNA cDNA Synthesis Kit (A28007; Applied Biosystems; Foster City, CA, USA), and then profiled. For the initial miRNA identification and selection studies, miRNAs were profiled with the Low-density TaqMan™ Advanced miRNA Human A Card (A34714; Life Technologies; Carlsbad, CA, USA). This array focuses on highly characterised miRNAs [58]. Specifically, it contains 377 different human miRNAs, including 1 endogenous control (hsa-miR-16-5p). Additionally, it also comprises 1 exogenous control (cel-miR-39-3p) and 1 negative control (ath-miR-159a) that does not amplify in human samples. Validated miRNAs were profiled using the Low-density Custom TaqMan ® Array Advanced MicroRNA Cards (Life Technologies; Carlsbad, CA, USA) (Supplemental Table S8). Taqman RT-qPCR amplification was performed in an ABI PRISM 7900HT Sequence Detection System (Applied Biosystems; Foster City, CA, USA) and the amplification curves were analysed using the Thermo Fisher Connect Dashboard (Life Technologies; Carlsbad, CA, USA). Only miRNAs with an expression threshold below 35 Ct cycles were considered. Ct values were obtained and processed by the standard curve method and normalised with the exogenous miRNA cel-miR-39-3p, according to the equation 2 −(Ct [target]−Ct [exogenous]) .

In Silico Systems Biology Analysis
Putative target genes of the selected miRNAs were determined using the databases TargetScan, miRmap, miRWalk and miRDB [18][19][20][21][22]. Considered target genes were those with a high likelihood to be regulated by miRNA according to their search algorithms. For each miRNA, only those targets predicted in all four databases were selected for gene-network analysis.
Initial enrichment pathway analysis of curated target lists of the identified miRNAs was pursued with the web-based software ShinyGO [23]. The top-10 Gene Ontology terms, sorted by average and rank (FDR and fold), were selected.
Additionally, the curated target lists were imported into Cytoscape 3.0 (online free available platform) [25] to build protein-protein interaction (PPI) networks according to the STRING database interaction data [24]. A confidence cut-off was set to 0.8, and the maximum number of additional interactions was limited to 30. In order to identify PPI-enriched clusters, a community cluster strategy (GLay) [26] was applied to the main connected component of each network. The resulting clusters were then analysed for functional enrichments according to the pathway annotations from the KEGG database [27,28].
For individual miRNA pathway analysis, a repository-based gene database (RBGD) was created by querying "heart failure" to the STRING Pubmed database interaction data (bioinformatics tool of the Cytoscape platform). Predicted targets of selected miRNAs were then matched to the RBGD list. Common genes and their neighbour interacting genes were imported to ShinyGO for pathway enrichment analysis. A schematic diagram of the protocol used can be found in Supplemental Figure S1.

Statistical Analysis
Sample normality was assessed with the Shapiro-Wilks test. Qualitative variables were described using the number of cases and percentages and analysed with the Chi-squared test, while quantitative variables were described using the median [interquartile range] and analysed with the U Mann-Whitney test. Receiver Operating Characteristic (ROC) curve analyses and the corresponding area under the curve (AUC), along with its 95% confidence interval (CI) were calculated to determine the discrimination potential of miRNAs (patients vs. reference subjects), and AUCs values (with cut-offs ≥0.7 and ≤0.3) were used for miRNA selection. Additionally, ROC curves for predicted probabilities were pursued to identify expression thresholds of miRNAs between patients and reference subjects, between cHF types (HFpEF and HFrEF), background aetiology (ischaemic and non-ischaemic), or need for rehospitalisation during the 4-6 years of follow-up, and the corresponding area under the AUC with its 95% CI was calculated. A p < 0.05 was considered statistically significant. SPSS Statistical Analysis System (version 26.0, IBM Corp; Armonk, NY, USA) was utilised.

Conclusions
In conclusion, in this hypothesis-generating study, we have identified and validated a differential pattern of eight miRNAs that are down-regulated in the systemic circulation of cHF patients. Additionally, we have shown by in silico analysis and systems biology approaches that these miRNAs are, when considered as a cluster, involved in molecular mechanisms and signalling pathways related to cHF. Nevertheless, further research on this topic, including pursuing mechanistic studies to validate our results, is warranted to corroborate these findings as well as to improve the understanding of cHF pathophysiology.

Study Limitations
This study is not free of limitations. Comparisons were pursued between cHF and a reference group of healthy subjects. This latter group was used to determine baseline miRNA levels for comparison purposes, as there are no established "normal" miRNA levels. Further studies to identify the effect of miRNAs in other populations are warranted.  Institutional Review Board Statement: The local ethics committee for clinical research approved the study (Ref. 16/44). Written informed consent was obtained from all participants before the initiation of the study and all research was conducted following the Declaration of Helsinki.
Informed Consent Statement: Informed consent was obtained from all subjects involved in the study before inclusion.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author upon reasonable request.