Dysregulations of MicroRNA and Gene Expression in Chronic Venous Disease

Chronic venous disease (CVD) is a vascular disease of lower limbs with high prevalence worldwide. Pathologic features include varicose veins, venous valves dysfunction and skin ulceration resulting from dysfunction of cell proliferation, apoptosis and angiogenesis. These processes are partly regulated by microRNA (miRNA)-dependent modulation of gene expression, pointing to miRNA as a potentially important target in diagnosis and therapy of CVD progression. The aim of the study was to analyze alterations of miRNA and gene expression in CVD, as well as to identify miRNA-mediated changes in gene expression and their potential link to CVD development. Using next generation sequencing, miRNA and gene expression profiles in peripheral blood mononuclear cells of subjects with CVD in relation to healthy controls were studied. Thirty-one miRNAs and 62 genes were recognized as potential biomarkers of CVD using DESeq2, Uninformative Variable Elimination by Partial Least Squares (UVE-PLS) and ROC (Receiver Operating Characteristics) methods. Regulatory interactions between potential biomarker miRNAs and genes were projected. Functional analysis of microRNA-regulated genes revealed terms closely related to cardiovascular diseases and risk factors. The study shed new light on miRNA-dependent regulatory mechanisms involved in the pathology of CVD. MicroRNAs and genes proposed as CVD biomarkers may be used to develop new diagnostic and therapeutic methods.

could also be involved in the pathogenesis of CVD and may be relevant as potential diagnostic and therapeutic targets.
In our work, integrated miRNA and gene expression analysis was applied to find potential robust biomarkers of CVD and to show the impact of miRNA-regulated genes on pathological processes governing CVD development.

Study Participants Characteristics
The study was performed in accordance with the Declaration of Helsinki and approved by the Ethics Committee at Medical University of Lublin (approval No. KE-0254/341/2015). Participants were recruited between February 2016 and May 2017. All subjects gave their informed consent for inclusion before they participated in the study. The CVD group consisted of 34 patients diagnosed and hospitalized only due to CVD, without any other diagnosed vascular diseases or comorbidities, in the Independent Public Clinical Hospital No. 1 in Lublin, Poland. The control group comprised of 19 healthy volunteers without any visible CVD characteristics and lack of comorbidities during examination. Detailed characteristics of included participants are presented in Table 1. Inclusion and exclusion criteria were evaluated by a vascular surgeon. Included CVD patients were examined using tourniquet test, auscultation and duplex ultrasound scanning. Venous reflux lasting longer than 1 s was classified as pathological. Patients diagnosed with symptoms classified according to CEAP as varicose veins (C2) of superficial veins (As) with primary etiology (Ep) and reflux pathophysiology (Pr) were included. The exclusion criteria were previous vascular surgery of lower limbs, insufficiency of deep veins, acute and chronic inflammation of veins, lower extremities arterial disease, coronary artery disease, cerebrovascular disease, aneurismal disease, myocardial infarction, hypertension, stroke, diabetes mellitus type 2, and pregnancy.
Body Mass Index (BMI), pain symptoms, ankle-brachial index, smoking habits and applied medical treatment were also evaluated ( Table 1).
The control group consisted of 19 healthy, nonsmoking volunteers (Table 1). Only subjects without blood flow disturbances and with normal morphology of veins in lower limbs, confirmed by physical examination and duplex ultrasound scanning, were included to the study as controls. Any symptoms, comorbidities and treatment of vascular diseases were indicated in the analysis of medical history of control subjects.

Study Material Preparation
Isolation of Peripheral Blood Mononuclear Cells (PBMCs) was performed from whole blood samples by density gradient centrifugation using Gradisol L reagent (Aqua-Med, Łódź, Poland). Small RNA fractions were isolated from PBMCs samples of 34 CVD patients and 19 controls using MirVana microRNA Isolation Kit (Ambion, Austin, TX, USA), according to the manufacturer's protocol. Total RNA was isolated from PBMCs of seven randomly selected CVD patients and seven randomly selected control subjects, using TRI Reagent Solution (Applied Biosystems, Foster, CA, USA) according to the manufacturer's protocol. For a more detailed description of study material isolation and assessment refer to [30].

Libraries Preparation and Sequencing
Small RNA libraries were prepared from 53 small RNA samples isolated from PBMCs of 34 CVD patients and 19 healthy controls. Technical limitations did not allow to perform transcriptome sequencing for all subjects included in the study, therefore transcriptome libraries were constructed from 14 total RNA samples isolated from randomly selected, representative subsets of PBMCs samples (seven from CVD patients and seven from healthy controls).
Small RNA and transcriptome libraries were prepared using Ion Total RNA-Seq Kit v2, Magnetic Bead Cleanup Module kit and barcoded with Ion Xpress RNA-Seq Barcode 01-16 Kit (all Life Technologies, Carlsbad, CA, USA), according to the manufacturer's protocol "Ion Total RNA-Seq Kit v2" revision B.0. Libraries were sequenced on Ion 540 Chips (Life Technologies) using Ion S5 XL System (ThermoFisher Scientific, Waltham, MA, USA). Small RNA and transcriptome raw sequencing data were aligned to 2792 human miRNAs from miRBase v21 (http://www.mirbase.org) and to 55,765 genes and splicing variants of hg19 human genome, respectively.
Detailed description of libraries preparation and sequencing procedures were included in our previous study [30].

Statistical Analysis
The differences between CVD and control groups were evaluated in terms of age and BMI using a two-sided Mann-Whitney U test (wilcox.test function in R), and in terms of sex and smoking using Fisher exact test (fisher.test function in R).
A Receiver Operating Characteristics (ROC) analysis implemented in pROC package 1.12.1 [41] (https://cran.r-project.org/web/packages/pROC/index.html) was used to evaluate the predictive value of selected miRNAs and gene transcripts for CVD classification.
Functional analysis of networked genes was performed using Database for Annotation, Visualization and Integrated Discovery (DAVID) 6.8 tool (https://david.ncifcrf.gov/) [47,48]. Default whole genome of Homo sapiens was applied as a background. All the terms of Kyoto Encyclopedia of Genes and Genomes (KEGG), Reactome and Genetic Association Database (GAD) databases associated with analyzed genes were harvested. Enrichment analysis of functional terms was proceeded for Gene Ontology (GO) terms separately for up-and downregulated genes.
All statistical procedures applied to this study were previously described in detail in [30]. Statistical analysis and visualizations were performed according to R code available in published reference manuals of used packages.
The summarized methodology applied in our study is presented on Figure 1.  [40] packages, respectively. Correlation analysis was performed using Spearman rank correlation test implemented in the Hmisc package 4.4-0 https://cran.r-project.org/web/packages/Hmisc/index.html.
A Receiver Operating Characteristics (ROC) analysis implemented in pROC package 1.12.1 [41] https://cran.r-project.org/web/packages/pROC/index.html was used to evaluate the predictive value of selected miRNAs and gene transcripts for CVD classification.
Functional analysis of networked genes was performed using Database for Annotation, Visualization and Integrated Discovery (DAVID) 6.8 tool https://david.ncifcrf.gov/ [47,48]. Default whole genome of Homo sapiens was applied as a background. All the terms of Kyoto Encyclopedia of Genes and Genomes (KEGG), Reactome and Genetic Association Database (GAD) databases associated with analyzed genes were harvested. Enrichment analysis of functional terms was proceeded for Gene Ontology (GO) terms separately for up-and downregulated genes.
All statistical procedures applied to this study were previously described in detail in [30]. Statistical analysis and visualizations were performed according to R code available in published reference manuals of used packages.
The summarized methodology applied in our study is presented on Figure 1.

Study Population Analysis
Characteristics of 34 patients with CVD and 19 non-CVD controls are presented in Table 1. A statistically significant difference between these groups was observed in age (p = 8.387 × 10 −3 ) and smoking history (p = 1.296 × 10 −4 ), which probably results from inclusion of healthy CVD-negative non-smoking individuals in control group. No statistically significant differences between CVD and control groups concerning gender and BMI were found (Table 1, Figure S1).

Primary Results
Detailed description of small RNA samples and small RNA libraries as well as the results of primary analysis of small RNA libraries sequencing data are presented in Table S1. Detailed description of transcriptome libraries and results of primary analysis of transcriptome libraries sequencing data are presented in Table S2. Plots depicting sequencing data quality, including boxplot of Cook's distances, MA plot and histogram of p values, regarding small RNA and transcriptome sequencing are presented in Figure S2 and Figure S3, respectively.

Differential Expression Analysis of miRNA
The comparison of miRNA expression levels between 34 CVD patients and 19 non-CVD controls was performed using DESeq2 and UVE-PLS methods and significantly dysregulated miRNAs selected by both methods were chosen.
DESeq2 analysis revealed 1034 differentially expressed miRNA transcripts in CVD subjects compared to controls. Ninety-six miRNA transcripts were differentially expressed with statistical significance p < 0.05 (Table S3). DESeq2 method is characterized by relatively high sensitivity, therefore a set of 49 differentially expressed miRNA transcripts (for 41 miRNAs) of higher significance (p < 0.01) was selected to limit the number of potentially false positive results.
In order to optimally filter uninformative miRNAs, the UVE-PLS method was applied to miRNA expression data of 1034 differentially expressed miRNA transcripts in CVD subjects. UVE-PLS analysis returned 48 informative miRNA transcripts (Table S4). Figure S4 shows the arrangement of prediction error and PLS components as well as cross-validated predictions versus measured values.
The set of 49 differentially expressed miRNA transcripts identified by DESeq2 method (with p < 0.01) and the set of 48 differentially expressed miRNA transcripts identified by UVE-PLS method were compared on the Venn diagram, revealing 34 miRNA transcripts common for both sets (Figure 2a). These 34 miRNA transcripts result in 31 miRNAs (22 upregulated and 9 downregulated), which constitute a proposed panel of potential miRNA biomarkers of CVD (Table 2). Differential expression of common 34 miRNA transcripts in CVD and control group is visualized on 3D PCA plot and in heatmap with Euclidean clustering (Figure 2b,c, respectively).     ROC analysis revealed that areas under ROC curves for 34 selected miRNA transcripts were covered in the range 0.930-0.757, indicating a high ability to distinguish patients with CVD from healthy subjects ( Table 2, Table S5 and Figure S5).
Correlation analysis between age and expression data of 34 selected miRNA transcripts was performed in CVD group in order to evaluate effect of age on these miRNAs (Table S6). Only hsa-miR-548ac was statistically significantly correlated with age (p = 0.0395) exhibiting weak and negative correlation (R = −0.35). The lack of statistically significant correlation of remaining 33 miRNA transcripts suggests their independency from age; however, further studies with larger populations should be performed to confirm this result.

Differential Expression Analysis of Genes
From CVD patient and non-CVD control groups, seven patients and seven controls were randomly selected for gene expression analysis. Similarly to miRNA, DESeq2 and UVE-PLS methods were applied to perform differential gene expression analysis. Significantly dysregulated genes revealed by both methods were selected.
DESeq2 analysis disclosed 23,204 differentially expressed genes in CVD subjects, comparing to controls. In total, 2719 genes presented statistical significance (p < 0.05). The risk of false positive results was decreased by selection of 183 differentially expressed genes with p < 0.00001 (Table S7).
Application of UVE-PLS analysis to gene expression data of 23,204 differentially expressed genes in CVD subjects compared to controls disclosed 74 informative genes (Table S8). Plot presenting the arrangement of prediction error and PLS components as well as plot of cross-validated predictions versus measured values were shown on Figure S6.
The set of 183 differentially expressed genes identified by DESeq2 method with p < 0.00001 and the set of 74 informative genes identified by UVE-PLS method were compared on a Venn diagram, revealing 62 genes common for both sets (Figure 3a). These 62 common genes constitute a proposed panel of potential biomarkers of CVD (Table 3). Clustering patterns of 62 genes in CVD subjects and controls was visualized on 3D PCA plot and in heatmap with Euclidean clustering (Figure 3b,c, respectively).
The ROC analysis showed that areas under ROC curves obtained for 62 selected genes were equal to 0.98 for RAC1P2 and RP11-318C24.1, and equal to one for the remaining 60 genes, indicating good precision of CVD classification (Table 3, Table S9, Figure S7).
Correlation analysis between age and expression data of 62 selected genes was performed in CVD group in order to evaluate effect of age on these genes (Table S10). HSPA8P1 (R = −0.82, p = 0.024), PTBP1P (R = −0.80, p = 0.031), TSC2 (R = −0.79, p = 0.033) and UBA52P5 (R = −0.77, p = 0.043) were statistically significantly and negatively correlated with age. Among these four genes, HSPA8P1, PTBP1P and UBA52P5 were disclosed as downregulated in CVD group, pointing to them as possible age-associated risk factors of CVD. A lack of statistically significant correlation of remaining 58 genes suggests their independency from age; however, further studies with larger populations should be performed to confirm this result.
To estimate the influence of cell subpopulations diversity in the PBMCs samples on results of gene expression analysis, the deconvolution procedure was carried out using "quanTIseq" and "MCPcounter" methods implemented in immunedeconv package. Both methods enable to perform comparisons between samples and "quanTIseq" allows also to make comparisons between cell types. Application of both methods to the gene expression data showed estimated proportions of 11 cell subpopulations in studied samples ( Figures S8 and S9). Although differences in proportions of particular cell subpopulations could be observed between samples, our data suggests that there is no significant impact of cell subpopulations composition in PBMCs samples on the study results.  P (FDR with Benjamini-Hochberg correction) and fold change values were obtained from DESeq2 analysis. PLS coefficients were obtained from UVE-PLS analysis. Areas under Receiver Operating Characteristics (ROC) curves (ROC-AUC) were received from ROC analysis. Genes were ordered according to increasing p values across groups of upregulated and downregulated genes. Genes without names assigned by HUGO Multi-symbol checker were termed as "Unmatched". Synonyms or previous gene symbols were put into brackets.

In Silico Identification of miRNA:Gene Interactions
Identification of miRNA:gene interactions between 31 selected miRNAs and 62 selected genes was performed in silico by multiMiR package. Twelve validated (Table S11) and 51 top 10%-predicted miRNA:gene pairs (Table S12) were returned from analysis. Interactions between miRNAs and their targets were visualized as a regulatory network containing 22 miRNAs and seven genes (Figure 4).

In Silico Identification of miRNA:Gene Interactions
Identification of miRNA:gene interactions between 31 selected miRNAs and 62 selected genes was performed in silico by multiMiR package. Twelve validated (Table S11) and 51 top 10%-predicted miRNA:gene pairs (Table S12) were returned from analysis. Interactions between miRNAs and their targets were visualized as a regulatory network containing 22 miRNAs and seven genes (Figure 4).

Discussion
High prevalence, multifactorial character and diverse symptomatology of CVD make this disease one of the major health problems worldwide [49]. There is a lack of sensitive and specific biomarkers for early detection of CVD and for monitoring disease progression. Altered expression of miRNA and its effects on regulation of gene expression make it a promising candidate for novel diagnostic and treatment approaches.
In the presented study, we conducted integrated, comparative analysis of miRNA and gene expression in PBMCs of patients with CVD and healthy controls. We applied Next Generation Sequencing and various statistical and bioinformatic tools to analyze expression profiles of CVD patients vs. controls and to search for genetic signatures of CVD (Figure 1). Most promisingly, 31 miRNAs (Table 2) and 62 genes (Table 3), which potentially may serve as biomarkers for CVD, were selected. The discriminative character of proposed biomarkers was reinforced by decreasing the p value of statistical significance in DESeq2 analysis (p < 0.01 for miRNAs and p < 0.00001 for genes, refer to Tables 2 and 3, respectively) and elimination of uninformative variables using the UVE-PLS method. In ROC analysis, we confirmed very solid diagnostic value of proposed biomarkers (Tables 2  and 3, Tables S5 and S9, Figures S5 and S7). This multi-step procedure with strict criteria was applied to obtain scientifically valid results and to eliminate RT-qPCR validation.
The presented miRNA and gene signatures of CVD expand the diagnostic perspective; however, the study population is relatively minor, especially in case of transcriptomic analysis and evaluation with larger cohorts is required to confirm diagnostic utility and discriminative ability of proposed CVD biomarkers.
Integration of miRNA and gene expression analysis in this study allowed us to determine the framework of miRNA regulatory network in CVD, introducing experimentally validated and predictive interactions between significantly differentially expressed miRNA and genes found in patients with CVD ( Figure 4). The findings of this study expand our understanding of miRNA functions and provide more insights into a complex network of post-transcriptional control in CVD.
To date, the topic of miRNA and gene expression analysis in venous disorders has not been sufficiently studied [33,34]. Cui et al. used microarray and RT-PCR to research upregulation of miR-34a, miR-202 and downregulation of miR-155 in vein biopsies from patients with CVI [33]. Our study confirmed that upregulated expression of miR-34a is indicative for chronic venous disorder.
Significantly higher expression of miR-16, miR-20a, miR-21, miR-106a, miR-203 and miR-130a was reported in skin biopsies from venous ulcers, comparing to normal skin specimens. Overexpression of these miRNAs was associated with inhibition of wound healing through targeting mRNA for EGR3, Vcl and LepR. Administration of miR-21 mimics in rat acute wound healing models leads to increase in infiltration of immune cells and reduction of epithelialization [34]. In our study, we demonstrated statistically significant upregulation of miR-21 (Table S3) in CVD subjects with varicose veins, which can be classified as an early stage of the disease. Therefore, enhanced signaling of miR-21 appears to be a factor involved in CVD progression both in early and advanced stages of disease development.
In other studies, miR-34a and miR-21 were significantly upregulated in plasma and in atherosclerotic plaques of patients with CAD (Coronary Artery Disease) [50,51], suggesting that overexpression of miR-34a and miR-21 also found in our study to is common for both arterial and venous pathology. Diagnostic and prognostic potential of increased level of miR-21 was also reported as oncomiR in many types of cancer, including Hodgkin lymphoma [52], head and neck squamous cell carcinoma [53] and lung cancer [54,55], while miR-34a acts as tumor suppressor and undergo epigenetic silencing in carcinogenesis [56].
Some studies indicated altered expression of genes in veins obtained from patients with CVD and venous ulcers [57][58][59][60] but none of the genes indicated in abovementioned research were found in our study. These differences could be related to the biological material type, inclusion criteria and methodological approach.
Presented miRNA regulatory network shows that upregulation of WNK1 is associated with downregulation of miR-181a-2-3p and miR-106b-3p (Figure 4). WNK1 is important for proper proliferation, chemotaxis and invasion of endothelial cells [61]. Deficiency of the WNK1 gene in mice induces embryonic lethality due to angiogenic and cardiovascular defects [61]. MiRNAs belonging to miR-181 family play a key role in vascular inflammation through regulation of NF-κB signaling, activation of endothelial cells and homeostasis of immune cells [62]. Upregulation of miR-181a in THP-1 cells exposed to oxidated LDL, lipopolysaccharides and phorbol myristate leads to reduction of foam cells formation and inflammatory cytokines levels, lower production of reactive oxygen species and inhibition of apoptosis [63,64]. A higher level of miR-181a in plasma was proposed as a biomarker of myocardial infarction [65]. Overexpression of miR-106b-3p was previously described in the late stage of endothelial replicative senescence [66]. Downregulation of miR-106b-3p was observed in human retinal pigment epithelium cells (ARPE-19) exposed to hydrogen peroxide (H 2 O 2 ) [67]. Therefore, we hypothesize that upregulation of WNK1 observed in our study, associated with downregulation of miR-181a-2-3p and miR-106b-3p is a possible mechanism promoting inflammation and aging in response to oxidative stress in CVD.
WNK1 was identified in our study as a putative target for downregulated miR-128-3p ( Figure 4). Zhou et al. showed that downregulation of miR-128-3p promotes endothelial cell proliferation through Ca 2+ and ERK1/2-Akt signaling [68]. Upregulation of WNK1 as a consequence of miR-128-3p downregulation may constitute the mechanism promoting endothelial cell proliferation in CVD. Downregulation of miR-128-3p could also be a factor alleviating inflammation, since it was reported to stimulate inflammation via targeting the TNFAIP3 gene and enhancing NF-κB signaling [69].
MiRNA regulatory network of CVD ( Figure 4) shows that miR-33a-5p targets PABPC3, which is a gene encoding a putative regulator of MEG3 stability [70]. MEG3 is a suppressor of tumor growth, inhibiting proliferation and promoting apoptosis in cancer cells [71] via activation of p53 [72]. Our study points to downregulation of PABPC3 with accompanying overexpression of miR-33a-5p as a factor which may affect MEG3 stability and in consequence promote cell proliferation in CVD.
Another upregulated miRNA found in CVD subjects, associated with MEG3 is miR-183-5p. This miRNA was reported to alleviate symptoms of hypoxia (such as a decrease in cell viability and migration) in H9c2 rat cardiomyocytes via targeting genes encoding MEG3 and p27 [73]. Thus, the upregulation of miR-183-5p observed in our study may be considered as another factor enhancing cell proliferation and hypertrophy in CVD. Additionally, in the presented miRNA regulatory network, upregulated miR-183-5p interacts with the upregulated WNK1 gene. This relationship may also lead to a hyper-proliferative effect, since WNK1 has been reported to be required for proliferation of endothelial cells and vascular smooth muscle cells [61,74].
In our study, downregulation of miR-769-5p in CVD patients was observed. Lower expression of this miRNA was reported in cancer cells and was associated with enhanced proliferation and reduced apoptosis [75,76]. Downregulation of miR-769-5p was also found in arterial samples of abdominal and popliteal arterial aneurysm [77]. In the miRNA regulatory network presented here, HDAC5 is proposed as a putative target of miR-769-5p. Xu et al. reported that HDAC5 mediates angiotensin II-induced MEF2 activation and vascular smooth muscle cells hypertrophy [78]. These findings suggest that HDAC5 upregulation, most likely resulted from downregulation of miR-769-5p, may be responsible for pathological vascular hypertrophy in CVD.
We indicated upregulation of miR-206 in CVD patients as an important target, which was already revealed to impair viability and migration of endothelial progenitor cells and to promote apoptosis through targeting VEGF [79]. MiR-206 mediate silencing of VEGF leading to inhibition of angiogenesis during Danio rerio development [80]. Upregulation of miR-206 in our study may thus point to inhibition of angiogenesis; however, two genes encoding downstream effectors of VEGF signaling, WNK1 [81] and CDS2 [82], were upregulated in CVD subjects. Upregulation of WNK1 and CDS2 may be a consequence of increase in VEGF expression mediated by chromatin remodeling caused by upregulation of HDAC5 [83,84], which was also observed in the current study. These findings may suggest an increase in angiogenesis in CVD independently from the VEGFR signaling. On the other hand, a different study indicated HDAC5 as a negative regulator of angiogenesis acting through genes for pro-angiogenic factors FGF2 and Slit2 [85], showing that reciprocal interaction network in CVD is more complex and requires further studies to answer all questions.
Downregulation of miR-30e-3p was previously observed in coronary microembolization and was associated with myocardial injury mediated by impairment of autophagy in cardiomyocytes [86]. Therefore, the downregulation of miR-30e-3p found in this study may be involved in vascular dysfunction related to the thrombosis in CVD.
Kim and collaborators reported at least 5-fold higher level of miR-33a in plasma from individuals with high risk of atherosclerosis, compared to non-high-risk subjects [87]. A significantly higher level of this miRNA was also described in PBMCs isolates from patients with CAD and was positively correlated with higher lipid levels and risk of atherosclerosis [88,89]. Therefore, the upregulation of miR-33a-5p found in our study may indicate higher risk of atherosclerosis in subjects with CVD, compared to non-CVD subjects.
Upregulation of circulating miR-92a-3p was observed in the plasma of both mild cognitive impairment and Alzheimer disease subjects [90] and was proposed as a biomarker of schizophrenia [91]. In our study, downregulated miR-92a-3p putatively targets upregulated PRRC2B (Figure 4), which is a gene encoding protein involved in brain development [92]. It suggests that some elements of neurodevelopmental process may also participate in CVD pathology, showing connections between vascular and neurodegenerative disorders.
The preliminary functional analysis of proposed transcriptomic biomarkers provides useful information on the pathogenesis of CVD. Functional analysis showed the association of three genes (CDS2, TBC1D22A and WNK1) with tobacco use disorder. The discriminative character of these genes in our study may be, to a certain extent, a result of smoking prevalence in the CVD group; however, smoking is also an established risk factor of CVD [2] and altered expression of CDS2, TBC1D22A and WNK1 could be involved in CVD development through smoking-related mechanisms. Many smoking-induced miRNAs were previously identified and associated with various pathological conditions, including carcinogenesis [93]. Presence of other smoking-related diseases was limited by inclusion patients diagnosed exclusively for CVD, with no other, especially cardiovascular, conditions detected during examination.
PBMCs constitute an important element of inflammation process in vascular diseases [11]; thus, transcriptional profiling of this cell pool should provide reliable information pertinent for vascular pathology. Another advantage of PBMCs is their accessibility through minimally invasive procedures, facilitating utilization in basic and clinical research. Despite of all advantages, co-prevalence of other diseases like chronic obstructive pulmonary disease, hypertension, diabetes mellitus and accompanying complications may bias the evaluation of expression profiles on a systemic scale. To cope with this bias, only subjects without co-existing conditions mentioned in the experimental section were included to the study. Such strict evaluation helped us to find systemic regulatory changes in miRNA and gene expression, which potentially are reflective of local changes in CVD. However, application of these exclusion criteria entailed construction of CVD and control groups with statistically significant differences in age and smoking history, which can introduce some biological bias to our results; therefore, further investigations with more balanced population groups should be performed.
We are aware that our research has several limitations. Studied PBMCs samples may differ in proportions of cell subpopulations (lymphocytes, monocytes), which may introduce a bias in miRNA and gene expression patterns. To assess the impact of this factor, we performed a deconvolution procedure using two methods: "quanTIseq" and "MCPcounter" implemented in an immunedeconv package, enabling comparison between cell types and between samples. The results of deconvolution do not indicate that differences in proportions of particular cell subpopulations across samples which included the CVD and control group had a significant impact on the study outcome ( Figures S8 and S9).
Due to technical constraints (server capacity), gene expression analysis was performed on the subset of participants subjected to miRNA expression analysis, which may affect the described transcriptomic effect of miRNA expression alterations in CVD patients. Despite this imbalance, we confirmed some previously validated interactions (Table S11) and determined with high probability other connections in signaling network (Table S12). Future studies should also include in vitro and in vivo validation of predictive interactions of presented miRNA regulatory network. Moreover, in the group selected for gene expression analysis, both miRNA and total RNA were isolated from PBMC specimens obtained from exactly the same subjects and represent the same physiological conditions probed at the same time and circumstances.
Despite applying strict criteria to select the most promising signatures of CVD, the biomarker role of particular miRNAs and gene sets should be confirmed in further studies with larger cohorts and with the application of other validation methods, such as RT-qPCR. Further experiments will involve investigations regarding clinical factors (e.g., stage of disease or medications) and elucidation of whether dysregulations of miRNAs and genes indicative for CVD were predictive of or responsive to the disease development.
The results obtained in our current study confirm the significance of miRNA-dependent epigenetic regulation in the pathogenesis of CVD. Although we proposed a novel biomarker panel of CVD, there is still a need for further research on the role of miRNA regulation in the CVD due to small sample size as well as a clear exploratory and hypothesis-generating character of the presented experiments.
The presented discoveries may be used in further fundamental research and may prove to be useful for clinicians and practitioners, providing new paths in diagnosis, differentiation and treatment procedures for CVD patients in future.

Conclusions
Owing to broad and detailed Next Generation Sequencing analysis one is able to draw some general conclusions about CVD. Analysis of microRNAs and genes dysregulated in CVD unveiled numerous terms related to general physiological processes and traits like inflammation, metabolism, aging as well as more specific ones like lipidomics, cardiovascular diseases and chemodependencies. Future research on much numerous groups of patients and controls would broaden our knowledge about cardiovascular diseases, enabling personalized approach to individual patients.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2077-0383/9/5/1251/s1, Figure S1: Distribution of age (A), BMI (B), sex (C) and smoking (D) in CVD and control group. Figure S2: Control plots for miRNA sequencing data. Figure S3: Control plots for transcriptome sequencing data. Figure S4: Plots for UVE-PLS analysis of miRNA expression in CVD group in comparison to control group. Figure S5: Receiver Operating Characteristics (ROC) curves for 34 miRNA transcripts selected as signatures of CVD. Figure S6: Plots for UVE-PLS analysis of genes expression in CVD group compared to control group. Figure S7: Receiver Operating Characteristics (ROC) curves for 62 genes selected as signatures of CVD. Figure S8: Results of gene expression deconvolution procedure for seven CVD patients and seven control subjects performed using "quanTIseq" method implemented to immunedeconv 2.0.0 package. Figure S9: Results of deconvolution procedure performed on gene expression datasets of 7 CVD patients (CVD) and 7 control subjects (Control) using "MCPcounter" method implemented to immunedeconv 2.0.0 package. Table S1: Measurements of small RNA samples and small RNA libraries as well as results of small RNA sequencing data analysis received from Ion Torrent small RNA Plugin v5.0.5r3. Table S2: Measurements of transcriptome libraries and results of transcriptome sequencing data analysis received from Ion Torrent RNASeqAnalysis plugin v.5.0.3.0. Table S3: The set of 96 differentially expressed microRNA transcripts resulted from DESeq2 analysis with p < 0.05 in 34 patients with CVD compared to 19 controls. Table S4: The set of 48 differentially expressed microRNA transcripts resulted from UVE-PLS analysis in in 34 patients with CVD group compared to 19 controls. Table S5: Results of ROC analysis for 34 miRNA transcripts selected as indicative for CVD. Table S6: Correlation analysis between age and expression of 34 selected miRNA transcripts in CVD group. Table S7: 183 differentially expressed genes resulted from DESeq2 analysis with p < 0.00001 in seven patients with CVD, compared to seven controls. Table S8: The set of 74 differentially expressed genes resulted from UVE-PLS analysis in seven patients with CVD compared to seven controls. Table S9: Results of ROC analysis for 62 genes selected as indicative for CVD. Table S10: Correlation analysis between age and expression of 62 selected genes in CVD group. Table S11: Twelve experimentally validated miRNA:gene pairs found among 31 miRNAs and 62 genes indicative for CVD using multiMiR package. Table S12: Top 10%-predicted miRNA:gene pairs found among 31 miRNAs and 62 genes indicative for CVD using multiMiR package.