Altered endothelial dysfunction-related miRs in plasma from ME/CFS patients

Myalgic encephalomyelitis/chronic fatigue syndrome (ME/CFS) is a complex disease characterized by unexplained debilitating fatigue. Although the etiology is unknown, evidence supports immunological abnormalities, such as persistent inflammation and immune-cell activation, in a subset of patients. Since the interplay between inflammation and vascular alterations is well-established in other diseases, endothelial dysfunction has emerged as another player in ME/CFS pathogenesis. Endothelial nitric oxide synthase (eNOS) generates nitric oxide (NO) that maintains endothelial homeostasis. eNOS is activated by silent information regulator 1 (Sirt1), an anti-inflammatory protein. Despite its relevance, no study has addressed the Sirt1/eNOS axis in ME/CFS. The interest in circulating microRNAs (miRs) as potential biomarkers in ME/CFS has increased in recent years. Accordingly, we analyze a set of miRs reported to modulate the Sirt1/eNOS axis using plasma from ME/CFS patients. Our results show that miR-21, miR-34a, miR-92a, miR-126, and miR-200c are jointly increased in ME/CFS patients compared to healthy controls. A similar finding was obtained when analyzing public miR data on peripheral blood mononuclear cells. Bioinformatics analysis shows that endothelial function-related signaling pathways are associated with these miRs, including oxidative stress and oxygen regulation. Interestingly, histone deacetylase 1, a protein responsible for epigenetic regulations, represented the most relevant node within the network. In conclusion, our study provides a basis to find endothelial dysfunction-related biomarkers and explore novel targets in ME/CFS.

Total RNA isolation. Plasma aliquots were thawed in a 37 °C bath for 1 min. Then, 200 µL were used for total RNA isolation using the miRNeasy Serum/Plasma Advanced Kit (#217,204, QIAGEN) in accordance with the operating manual. C. elegans miR-39 miRNA mimic (#219,610, QIAGEN) was added during the RNA extraction to monitor RNA recovery and reverse transcription efficiency. RNA aliquots of 20 µL were stored at -80 °C.
First-strand cDNA synthesis and real-time PCR. miRCURY LNA RT Kit (#339,340, QIAGEN) was used to synthesize cDNA. 2 μL of RNA were used as a template for the reverse transcription reaction (60 min at 42 °C). Then, the samples were incubated for 5 min at 95ºC to heat inactivate the reverse transcriptase (peqSTAR 2X Gradient Thermocycler, Peqlab Biotechnologie GmbH). miRCURY SYBR Green PCR Kit (#339,346, QIA-GEN) and miRCURY LNA miRNA PCR Assays (#339,306, QIAGEN) were used to evaluate the expression of hsa-miR-21-5p (YP00204230) (Gene ID: 406,991), hsa-miR-34a-5p (YP00204486) (Gene ID: 407,040), hsa-miR-92a-3p (YP00204258) (Gene ID: 407,048), hsa-miR-126-3p (YP00204227) (Gene ID: 406,913), and hsa-miR-200c-3p (YP00204482) (Gene ID: 406,985) relative to cel-miR-39-3p (YP00203952) by means of real-time PCR. 3 μL of cDNA (diluted 1:20) were used as a template for each reaction (final volume 10 μL). PCR cycling conditions included initial heat activation (2 min at 95 °C), 2 step cycling, denaturation (10 seg at 95 °C), combined annealing/extension (1 min at 56 °C), 40 cycles, and melting curve analysis (60-95 °C). The real-time PCR was performed in a 96-well thermocycler (StepOnePlus Real-Time PCR System, Applied Biosystems-Life Statistical analysis of the UKMEB data derived from plasma. To compare basic characteristics between HC and patients with ME/CFS, t-tests for two independent groups were used for normally distributed data. The Mann-Whitney rank test for two independent groups was used otherwise. Given that CRP data had many values below the detection limit of 1, a linear regression based on truncated log-normal distribution was used to compare the different groups under analysis. In this analysis, it was fitted to two models, one including an intercept and a group variable and another one including an intercept only. These models were then compared with Wilks' likelihood ratio test. P value < 0.05 indicated no differences between the three groups.
A pairwise correlation analysis between the five miR under analysis were conducted using Pearson's correlation coefficients given that the respective log 10 -transformed data followed Normal distributions and the statistical relationship between any pair of miR was linear approximately. A similar analysis was also performed to assess the pairwise correlation between each miR and each clinical, blood and SF-36 variable but using Spearman's correlation coefficient.
A principal component analysis on the miR data was used to assess the similarity of the study participants in terms of these variables. In this analysis we only used the first principal components because these components could account for more than 90% of the miR data variability. In addition, two linear discriminant analyses based on the miR data were performed to distinguish patients with ME/CFS from healthy controls and to distinguish severely affected patients from patients with mild or moderate symptoms. In general, linear discriminant analysis when applied to two-group case aims to find the best linear combination of a given set of variables that could distinguish two groups. The key output of this analysis is the probability of classifying each individual in the respective group given the estimated linear combination of the variables.
To detect differential expression of each miR between the three study groups, we estimated linear regression models defining the log10-transformed data of each miR as the outcome variable, and age, gender, body mass index, and two possible disease status variables (healthy vs ME/CFS and healthy vs ME/CFSmm vs ME/CFSsa) as the respective covariates. Three and four-way ANOVA were used to test the statistical significance of the covariates in the models including three and four covariates, respectively. Raw p values were then adjusted in order to ensure a global false discovery rate of 5% using the Benjamini-Hochberg procedure. Finally, the assumptions associated with linear regression modelling were confirmed by a residual analysis including normality assessment by Lilliefors and Shapiro-Wilk tests and different diagnostic plots (e.g., scatterplots of fitted values against standardized residuals).
The last step of the analysis was to estimate an empirical statistical power for detecting a group effect using the different estimated models including such effect as the true models for the respective data. With this purpose, we used the so-called parametric Bootstrap 36 based on the following algorithm: (1) for each miR, simulate new values for each study participant according to the Normal distributions associated with the respective linear regression models estimated from the real data; (2) estimate the same models in the simulated data; (3) for each miR, calculate the p value to test the significance of the group effect as performed for the real data; (4) adjust the p values associated with each miR using the Benjamini-Hochberg procedure; (5) repeat the previous steps 1000 times in order to obtain 1000 data sets and to perform the respective analysis; and (6) estimate the power associated with each miR as the proportion of data sets which the group effect was statistically significant using adjusted p values as described for the analysis of the real data.
Analysis of public microarray data derived from PBMCs. We re-analyzed a microarray data set from a previously published study 17 . In a nutshell, this data set consisted of 15 patients with ME/CFS and 29 HC matched for age and gender. Mean age of patients and HC was 43.8 and 44.76 years old, respectively. The proportion of males was 0.20 and 0.17 in patients and HCs, respectively. Patients with ME/CFS were diagnosed by the Fukuda 1994-CDC, but also fulfilled the 2003 Canadian Consensus Criteria. Additional information of the study participants can be found in the original reference 17 . The microarray data referred to log10-transformed intensities for around 385 miR as available in the Ambion Bioarray V1.0. We focused our analysis on data concerning the five miR of interest. Again, to detect differential expression of each miR between HC and patients with ME/ CFS, we estimated linear regression models where the intensity data associated with each miR was the outcome variable while age, gender, and disease status (HC vs ME/CFS) were considered as the covariates. The diagnostics of the estimated regression model were conducted as described above for the analysis of plasma-derived data sample. A three-way ANOVA was used to test the statistical significance of the covariates in the estimated models. Again, raw p values were adjusted for multiple testing using the Benjamini-Hochberg procedure. The adjusted p values ensured a global false discovery rate of 5%. The whole data set is publicly available from the National Center for Biotechnology Information (NCBI) Geo data sets under the reference code GSE70371. The annotation file of the respective microarray can be found under the accession code GPL3444 in the same data repository. This file was used to identify each miR in the respective data set.
Statistical software. All statistical analysis was performed in the R software version 4.0.2 37 (https:// www.rproje ct. org/). The following R packages were used: "tcensReg" to analyze the CRP data, "nortest" to test the normal distribution of quantitative variables and residuals of the linear regression analyses, "PerformanceAnalytics" www.nature.com/scientificreports/ to produce the correlation matrix shown in Fig. 2, and "MASS" to adjust raw p values for multiple testing. The R scripts used in this paper are freely available from NS upon request.
Integrative network analysis. Potential miR-mRNA interactions were determined using miRTarBase 7.0 (http:// mirta rbase. mbc. nctu. edu. tw/ php/ index. php), a database that estimates interactions based on experimentally validated assays (e.g. reporter gene assay, qPCR, and western blot). Protein-protein interactions (PPI) were also estimated using the GeneMANIA prediction server, which uses experimentally confirmed physical associations 38 . In order to analyze the PPIs topology network, a graphical representation was generated with Cytoscape software (version 3.7.2) using the Organic layout. The degree and betweenness centrality of proteins in the PPIs network were also calculated. The degree score (DS) corresponds to the count of how many interactions a node has in the whole network. On the other hand, betweenness centrality (BC) measures the frequency with which a node lies between the shortest communication path of all other possible pairs of nodes within a network, determining which are the key nodes for communication within the network. In this context a protein with a high DS and BC is likely an essential protein.
Over-representation analysis. Over representation analysis (ORA), an approach to determine whether known biological functions or processes are enriched in an experimentally-derived gene list, was performed calculating a p value using hypergeometric distribution 39 . The list of miRs targets previously obtained using miRTarBase 7.0 were used as a target group. In parallel, a full list of genes obtained from the NCBI database was established as a background set of genes that function together in a known biological pathway. The analysis was performed using five random lists of genes obtained from miRTarBase 7.0 to ensure the biological processes are represented in the gene list more often than expected by chance. Random lists were generated using the random python module to determine if a significative gene ontology term could be obtained by chance.

Results
Study participant characteristics. The demographic and clinical characteristics of the participants are summarized in Table 1. The final age-matched cohort comprised 58 patients with ME/CFS divided into mild/ moderate (mm) (ME/CFSmm: 28; females: 71%) and severely affected (sa) ones (ME/CFSsa: 30; females: 80%). Healthy controls included 29 participants (Females: 41%). Fatigue severity scores were higher both in ME/ CFSmm and ME/CFSsa in relation to HC (p < 0.0001). As expected, clinical assessments related to pain and fatigue were higher in patients with ME/CFS compared to HC, with higher scores in ME/CFSsa comparing to ME/CFSmm (p < 0.0001). In contrast, waist circumference (p = 0.0124), body metabolic rate (p = 0.0044), and body muscle impedance (p = 0.0011) were lower, while body fat bioimpedance was higher (p = 0.0235) both in ME/CFSsa and ME/CFSmm patients compared to HC participants. Diastolic and systolic pressure, along with BMI were similar across the groups. Blood tests showed increased platelets (p = 0.0264), basophils (p = 0.0359), erythrocyte sedimentation rate (ESR) (p = 0.0110) and reduced levels of creatinine (p = 0.0045) and creatine phosphokinase (CPK) (p < 0.0001) both in ME/CFSmm and ME/CFSsa in relation to HC. The median obtained from the physical functioning scale variables (SF-36 questionnaire) such as physical functioning, role physical, bodily pain, general health, vitality, social functioning, role emotional, mental health, and both physical and mental component summary were significantly reduced in both ME/CFSsa and ME/CFSmm patients when compared to HC.
Analysis of miR expression data from plasma samples. We then used plasma samples from our cohort to analyze the relative abundance of the set of miRs reported to modulate the Sirt1/eNOS axis. We detected that the levels/quantities of miR-21, miR-34a, miR-92a, miR-126, and miR-200 were all increased in both ME/CFS groups compared to HC (Fig. 1), and that these miRs are highly and positively correlated with each other irrespective of the study group (Fig. 2). In line with this result, the first two principal components explained more than 90% of the variation in the data (Fig. 3). Interestingly, the principal component analysis showed that some patients with ME/CFS were clearly different from HC in the abundance of these miRs (Fig. 3A), but it seems not to be related to disease severity (Fig. 3B). To complement this visual interpretation of the principal components, linear discriminant analysis suggested that miR data could correctly classify 60.8% (n = 31) of patients with at least 80% probability (Fig. 3C). In contrast, only 30.4% (n = 7) of HC could be classified as such with the same minimum probability. When linear discriminant analysis was used to distinguish between ME/CFSsa and ME/CFSmm patients, the corresponding classification probabilities were in the vicinity of 50% for individuals from each disease group (Fig. 3D). In conclusion, miR expression data may be a supportive biomarker to discriminate between ME/CFS patients and HCs, but cannot predict disease severity. Spearman correlation coefficients were used to assess for an association between the increased Sirt1/eNOSrelated miRs and each clinical parameter, in order to establish novel relationships that could help in better discriminating between patients and controls. In HC, the correlation coefficient involving ESR tends to be positive for miR-34, miR-92a and miR-200c. Diastolic pressure tends to be positive for miR-34, and miR-92a. The amount of white blood cells and neutrophils tends to be positively correlated with the level of miR-34a (Fig. 4A). The correlation estimates associated with ME/CFS patients (Fig. 4B) and all the participants from our cohort (HC and ME/CFS patients) (Fig. 4C) tends to be around zero for all miRs and clinical parameters, suggesting that the heterogeneity of ME/CFS patients might mask stronger associations related to the analyzed covariates. The association analysis shows the comparison between groups after adjusting by either age/gender, BMI or age/ gender/BMI (Fig. 5A). All associations are statistically significant irrespective of their adjustment for possible confounders (Fig. 5A, Supplementary Table 1), highlighting the need to consider these analyses to obtain more reliable findings. The empirical power of detecting a group effect assuming the estimated models as the true ones www.nature.com/scientificreports/   www.nature.com/scientificreports/ for the data was in most cases higher than 80% in unadjusted analysis and in the analysis controlling for age and gender (Fig. 5B). This power tended to be lower in the analysis controlling for age, gender, and BMI. As expected, the lowest power was obtained for miR-34a which was already on the borderline of statistical significance after adjusting for multiple testing. In the worst-case scenario of this miR, the power of detecting a group effect was estimated at 35% for the analysis controlling age, gender, and BMI when using ME/CFSsa, ME/CFSmm and healthy controls as the group covariate.
Analysis of microarray data from PBMCs. To provide additional evidence for the above findings on the UKMEB cohort of participants, we re-analyzed previously published microarray data from 15 patients with ME/CFS and 29 HC matched for age and gender 17 . We found that the five miRs under analysis were increased in patients with ME/CFS in comparison with HC (Fig. 6A, Supplementary Table 1). This increase was statistically significant for miR-21, miR-34a, miR-92a, and miR-126 after adjusting for multiple testing and controlling for possible confounding effects of age and gender (Fig. 6B, Supplementary Table 1). With respect to miR200c, the increase in expression was in the borderline of statistical significance (FDR-adjusted p value = 0.06). In conclusion, the five miRs are increased in both plasma (present data) and PBMCs 17 from different cohorts, supporting their altered expression in ME/CFS.

Bioinformatics analysis.
Bioinformatics analyses were carried out to visualize targets and biological processes associated with the set of miRs. Based on miR-mRNA interaction analyses we observed 182 genes associated with miR-21 (n = 49), miR-34a (n = 61), miR-92a (n = 17), miR-126 (n = 33), and miR-200c (n = 38) (Fig. 7, Supplementary Table 2), indicating they are eventually regulated by more than one miR. After integration with protein-protein interaction data, we observed a cluster of 135 proteins with experimental evidence of physical interaction. Thus, our set of increased miRs is predicted to target highly interconnected proteins, suggesting a functional module (Fig. 7). Interestingly, histone deacetylase 1 (HDAC1) represented the most relevant node based on centrality measures (Fig. 7, Supplementary Table 2). Later, we identified several biological processes associated with miR-21-5p, miR-34a-5p, miR-92a-3p, miR-126-3p, and miR-200c-3p from a gene list validated experimentally. Functional pathway enrichment analyses showed that 5 out of 20 most over-represented biological processes are linked to regulation of vasculature development, followed by reproductive structure development, response to oxidative stress, ERK1/ERK2 cascade, and response to oxygen levels, respectively (Fig. 8).
Additionally, 232 biological processes derived from the same analysis are also provided (Supplementary Table 3).
Overall, the set of miRs selected for our study is mainly associated with endothelial function through regulation of highly interconnected proteins, where HDAC1 and Sirt1 are predicted as central players.

Discussion
Despite ME/CFS has been studied for decades with significant progress achieved in the field, its etiology remains elusive. Accordingly, ME/CFS is diagnosed based on symptoms-related criteria which also leads to either undiagnosed or misclassified cases due to the heterogeneity of symptoms. Consequently, there is no single objective diagnostic or therapeutic biomarker, and current treatment is mainly focused on alleviating the complex symptomatology 1 . Therefore, there is an urgent need to identify reliable biomarkers for ME/CFS that allow a more precise diagnose and follow-up of its progression 40 . Thus, strategies to improve the reproducibility among studies such as the utilization of uniform clinical and research criteria, standardization of sample collection along with proper statistical analysis are certainly crucial to overcome this challenge. We used plasma samples from healthy controls and age-matched mild/moderate and severely affected ME/ CFS patients, obtained from the UKMEB 41 . In agreement with the stratification of groups in this study, all clinical parameters related with pain and fatigue were higher in ME/CFSsa than ME/CFSmm than HCs. In direct contrast, parameters associated with muscle mass (body muscle impedance) and muscle activity (creatinine and CPK) were lower in ME/CFSsa than ME/CFSmm than HCs, whereas body fat mass showed the opposite pattern, as expected due to the debilitating condition of the syndrome and the possibility of a gradual filling of the space of atrophic muscle tissue by fat tissue, as it occurs in other myopathies 42 , which, however, does not alter BMI despite a tendency towards a reduced body weight. As in other reports 43 , blood pressure and CRP were not different within our study population, in agreement with the immunological heterogeneity across ME/CFS cohorts. Furthermore, the non-specific marker of inflammation ESR, was found increased in ME/CFS patients compared to HCs, as well as reported in a different cohort from UKMEB 44 .
Notwithstanding the current knowledge that inflammation promotes ED in different diseases 6 , few studies have addressed endothelial function in vivo in patients with ME/CFS [7][8][9]43 . From a clinical perspective, Newton and colleagues 9 reported ED in large and small arteries from ME/CFS patients using flow-mediated dilatation (FMD), linked to increased serum levels of CRP. Conversely, Moneghetti et al. neither report significant differences in CRP levels nor endothelial function assayed by peripheral arterial tonometry (EndoPAT) between ME/ CFS patients still able to exercise and sedentary controls 43 . Scherbakov et al. also evaluated endothelial function by EndoPAT, showing ED associated with disease severity and immune symptoms in 51% of ME/CFS patients 7 .
In line with these findings, Sørland et al. have recently reported ED in 40% of ME/CFS patients by using FMD 8 . Interestingly, our results show that ~ 61% of ME/CFS patients from our cohort might be classified as patients with at least 80% probability of suffering from ED, based on the increased set of miRs reported to regulate endothelial function via the Sirt1/eNOS axis. Although these findings also suggest that ED might be a trait observed in a subset (~ 40-60%) of ME/CFS patients, further studies focused on evaluating both endothelial function in vivo and this set of Sirt1/eNOS-related miRs in plasma from the same cohort might be of great value to support this assumption, offering a more sensitive approach to detect ED earlier. www.nature.com/scientificreports/ The enzyme eNOS has a crucial role in endothelial function as a key mediator of vasodilatation by converting L-arginine into NO and L-citrulline. On the contrary, uncoupled eNOS and/or reduced eNOS activity induced by oxidative stress and inflammation triggers ED 6 . Sirt1 is another key modulator of endothelial homeostasis by activating eNOS to promote endothelium-dependent vasodilatation 11 , whose expression/activity is reduced during inflammation 10 . Over the last years, miRs have emerged as interesting candidates to identify biomarkers in ME/CFS [13][14][15][16][17] , since these noncoding RNAs play widespread roles in regulating gene expression 12 . Microarray analyses have identified several miRs differentially expressed in ME/CFS patients compared to HCs, also supporting by in silico analyses that inflammation-related processes and immune abnormalities play a role in its     www.nature.com/scientificreports/ Figure 5. Association analysis between miR expression and study groups. (A) Association analysis adjusting for a global false discovery rate of 5% and controlling for age, gender, and BMI. Each dot represents the −log 10 (p value) of a specific statistical association test while the dashed line represents −log 10 (0.05) above which it was consider a statistically significant association (i.e., −log 10 (p value) > −log 10 (0.05)). P values were adjusted for a false discovery rate (FDR) of 5%. (B) Power analysis using a parametric Bootstrap approach under the assumption that the estimated models in A were the true ones for the data. In this analysis, the probability of detecting a group effect was estimated by the proportion of simulated data sets in which the group effect was statistically significant after adjusting for an FDR of 5%. The statistical analysis was performed in the R software version 4.0.2 (https:// www.r-proje ct. org/). www.nature.com/scientificreports/ pathophysiology 14 . Then, we compared the expression of miR-21, miR-34a, miR-92a, miR-126, and miR-200c in plasma samples from ME/CFSmm and ME/CFSsa versus HCs, obtained from the UKMEB 41 . We detected increased plasma levels of each miR in ME/CFS patients compared to HCs, independent of disease severity. In order to provide additional evidence, we also re-analyzed available microarray data in PBMCs from a different cohort of ME/CFS patients 17 . miR-21, miR-34a, miR-92a, and miR-126 were also found increased in patients with ME/CFS in comparison with HC, while miR-200c showed the same trend but at the borderline of statistical significance. Altogether, these findings not only support the fact that evaluating this set of miRs, reported to regulate the endothelial function via the Sirt1/eNOS axis, might be used to differentiate a subset of ME/CFS patients from HCs but also highlight the importance of available data sets to compare findings across studies.
Taking into account the immunological heterogeneity across patients during the course of the illness 2,3 , several studies have investigated immune patterns in ME/CFS, often with conflicting results 45 . Nevertheless, pro-inflammatory cytokines, such as TGF-β, TNF-α, IL-1β, and IL-6 seem to play a role related to illness duration and severity 2,3 . In agreement with a recent study carried out in PBMCs and extracellular vesicles from ME/ CFSsa patients 15 , we found increased levels of miR-21 in plasma samples from both ME/CFSmm and ME/CFSsa patients compared to HCs. miR-21 is described to down-regulate the Sirt1/eNOS axis via TGF-β 20 and TNF-α 28 pathways in ECs and endothelial progenitor cells (EPCs), respectively. Increasing experimental and clinical data have also suggested that asymmetric dimethylarginine (ADMA), an endogenous eNOS inhibitor, promotes ED by reducing NO bioavailability 46 . Indeed, TNF-α-induced miR-21 up-regulation is linked to increased ADMA concentration and reduced eNOS activity in ECs 27 . L-citrulline, the by-product of NO synthesis from L-arginine by eNOS, ameliorated ADMA-induced ED by enhancing eNOS function and reducing oxidative stress 47 . Interestingly, L-citrulline was found significantly lower in plasma from ME/CFS patients compared to HCs 48 . Since the main source of L-citrulline is the gut, its circulating levels are used as an indicator of gut function 49 . Gut microbial dysbiosis has been proposed to be involved in ME/CFS 50 , and the degree of fatigue associated with the decrease of plasma L-citrulline in other pathological settings of dysbiosis 49 . Thus, increased plasma levels of ADMA 47,51 and decreased plasma levels of L-citrulline 48 may negatively affect NO production in ME/CFS patients. Additionally, miR-21 is thought to act as a positive but indirect regulator of Foxp3 expression in human regulatory www.nature.com/scientificreports/ CD4 + T cells (Tregs) 52 , a key cellular population of the adaptive immune system that controls autoimmunity and autoreactive T cells 53 . Accordingly, miR-21 ameliorated the clinical severity of experimental autoimmune encephalomyelitis in mice via a mechanism involving Tregs 53 . Therefore, it is possible to hypothesize that an overexpression of miR-21 could promote the expression of Foxp3 in naive CD4 + T cells and their conversion to the Treg pool. This putative Treg-cell conversion could explain the increase in the proportion of these cells in patients with ME/CFS observed in some studies 54 .
In line with this observation, miR-34 also regulates Tregs via NF-κB signaling 55 , a well-described pathway that induces inflammation by inhibiting Sirt1 10 . miR-34a also seems to play a key role in the regulation of endothelial cell inflammation by differentially modulating the antagonistic crosstalk between NF-κB and Sirt1 26 . Moreover, the negative regulation of miR-34a over the Sirt1/eNOS axis is associated with the premature senescence of ECs 25 . Interestingly, a randomized clinical study reported a negative correlation between miR-34 and Sirt1 expression Figure 7. Visualization of miR-target interaction network. miR-21-5p, miR-34a-5p, miR-92a-3p, miR-126-3p, and miR-200c-3p are displayed in yellow squares, while white ovals represent their targets. Dotted yellow lines and solid blue light ones indicate miR-mRNA and protein-protein interaction, respectively. The graphical representation was generated with Cytoscape software version 3.7.2 using the Organic layout (https:// cytos cape. org/). www.nature.com/scientificreports/ in patients with coronary artery disease 29 . Thus, considering we found miR-34a increased in ME/CFS patients compared to HCs, further studies are certainly needed to elucidate its potential contribution to the immune and endothelial abnormalities observed in ME/CFS patients. Our results also showed increased circulating miR-92a both in ME/CFSmm and ME/CFSsa patients, in agreement with a recent study using PBMCs 16 . In vitro, inhibition of miR-92a attenuates ED by decreasing the release of TNF-α and IL-6 by ECs 24 . In addition, miR-92a overexpression decreased eNOS activity and release of NO 23 , whereas its inhibition restored eNOS-mediated endothelial function by increasing Sirt1 expression in ECs 22 . In observational clinical studies, the circulating miR-92a positively correlated with microvascular coronary ED 30 . Moreover, miR-92a was inversely correlated with FMD and positively associated with IL-1β and CRP in serum from patients with coronary artery disease 22 .
In line with our findings, miR-126 was also found increased in extracellular vesicles from ME/CFSsa patients 15 . miR-126 is described to exert an endothelial protective role against hypoxia/reoxygenation-induced injury, oxidative stress, and TNF-α by activating the Sirt1/eNOS axis in ECs 21 . The circulating miR-126 was correlated to the improvement of endothelial function in male obese adolescents after exercise and diet control 31 . In addition, miR-126 reduces oxidative stress, IL-6, and TNF-α and activates both eNOS and vascular endothelial growth factor (VEFG) in ECs 21 . Recently, VEGF was found decreased in serum from ME/CFS patients 56 , an interesting finding considering that VEGF promotes survival and stability of ECs 57 . Altogether, higher plasma levels of miR-126 in ME/CFS patients might involve a compensatory mechanism against ED.
A large body of evidence suggests an autoimmune etiology in a subset of ME/CFS patients 7,58 . For instance, the M 3 muscarinic acetylcholine receptor increases NO synthesis mediated by acetylcholine in ECs 59 , but autoantibodies against this receptor were found significantly higher in ME/CFS patients compared to HC 60 . Several autoantibodies, including those against the M 3 muscarinic acetylcholine receptor, have been identified in people suffering from postural orthostatic tachycardia syndrome (POTS), characterized by tachycardia after moving from the supine to an upright position 61 , a condition also linked to ME/CFS 62 . Interestingly, miR-200c overexpression attenuates acetylcholine-induced endothelium-dependent relaxation in human renal arteries 19 . Thus, miR-200c might be related to acetylcholine-related ED. In this regard, circulating miR-200c was reported increased in atherosclerotic patients associated with reduced Sirt1/eNOS expression 32 . Additionally, miR-200c up-regulation induced by oxidative stress reduced Sirt1/eNOS axis activity 18 , causing apoptosis and senescence in ECs 63 .
Several bioinformatics analyses were also performed to visualize additional targets and biological processes associated with the set of miRs. Based on in silico analyses we validated that endothelial function-related Figure 8. Biological processes related to miR-21-5p, miR-34a-5p, miR-92a-3p, miR-126-3p, and miR-200c-3p. Dot plot of functional enrichment analysis for the top 20 over-represented biological processes (BPs) related to our selected miRs. Dot sizes represent the number of genes (count) related to a particular BP. Gene Ratio is the number of genes found enriched in each category over the number of total genes associated to that BP. Dot colors represent the adjusted p values. Dot plots were generated in R software version 3.6.1 using the clusterProfiler package version 3.14.3 (https:// bioco nduct or. org/ packa ges/ relea se/ bioc/ html/ clust erPro filer. html). www.nature.com/scientificreports/ signaling pathways are closely related to these miRs, including oxidative stress and oxygen regulation. Interestingly, HDAC1 was the most relevant node within our network built with miR-21, miR-34a, miR-92a, miR-126, and miR-200c. Histone deacetylases remove acetyl-lysine marks on proteins, especially in histones, leading to chromatin condensation and reduced levels of gene expression 64 . Changes in chromatin compaction through altered HDAC1 64 along with other epigenetic alterations previously detected in ME/CFS patients 65 , might play also a role in this syndrome. Sirt1 deacetylates and activates HDAC1, but cell stress induces both Sirt1 degradation and HDAC1 acetylation, rendering HDAC1 less active 66 . Other HDAC isoforms have been reported increased in PBMCs from ME/CFS patients 67 . Notably, up-regulation of HDAC1 has been shown to reduce NO production by deacetylating eNOS in ECs in vitro 68 . Despite these findings might represent potential strategies for preventing ED via HDAC1, they should be interpreted with caution because both up-regulation of HDAC1 68 and pharmacological inhibition of HDACs 69 negatively affect endothelial function. Thus, we need further studies to confirm any alteration in HDAC1 amount and/or activity directly related to ED in our cohort. However, based on the mutual regulation between Sirt1 and HDAC1 66 , a better analysis of endothelial function should address HDAC1/Sirt1/eNOS axis as a novel platform to dissect the altered regulatory pathways leading to ED. Noteworthy, we found similarities in the expression of our target miRs in plasma compared to PBMCs from another cohort 17 , but their origin is unknown. From a biological aspect, miRs are produced, released, and incorporated by different cell types, playing widespread roles in regulating gene expression throughout the body 70 . NF-κB is also involved in the modulation of a wide range of pleiotropic responses 10 . Thus, it is not surprising that this master transcriptional factor, well-known to modulate both immune and inflammatory responses 71 , might also participate in puzzling signaling pathways ranging from metabolic syndrome 72 to depression 73 in crosstalk with several miRs 74 , including miR-21 75 , miR-34a 55 , miR-92a 22 , and miR-126 76 . Based on the heterogeneous nature of ME/CFS 1 , we are aware that these increased miRs are not exclusively related to this syndrome, but they reflect an alteration, not related to severity, that might be further exploited to identify novel pathways with diagnostic and therapeutic value in ME/CFS. Since the ECs-related function varies from acute to chronic inflammation and even during the switch from innate to adaptive immunity 4 , further basic and clinical studies are certainly required to elucidate whether our studied miRs are involved in immune abnormalities in ME/CFS.
To our knowledge, this is the first report that collectively evaluates a set of circulating miRs associated with the Sirt1/eNOS axis, providing a basis for further studies to find ED-related biomarkers in ME/CFS. Overall, lowinvasive detection of this set of miRs in plasma might help to identify ME/CFS patients, but their abundance per se cannot discriminate between ME/CFSsa and ME/CFSmm patients in this cohort. Our findings also support the hypothesis that endothelial homeostasis is an underestimated and partially addressed process which might play an important role in the complex pathophysiology of ME/CFS. Considering the pivotal function exerted by ECs and the fact that they are ubiquitous in tissues throughout the body, our study not only supports other reports associating ME/CFS with ED 7-9 but also advances our understanding about potential new players not reported so far. In conclusion, we propose that a combination of clinical evaluation of endothelial function by FMD and/or EndoPAT, along with the detection of circulating Sirt1/eNOS-related miRs, might allow a more sensitive characterization of ED in a subset of ME/CFS patients for better stratification, which will certainly translate into improved treatment possibilities.
We are aware that these findings are not yet transferable into the clinical setting. Accordingly, observational studies with follow-up measurements in a bigger cohort are certainly needed to improve the power of the analysis and confirm our findings. In this line, we agree with Jason et al. 77 that it is crucial to include minimum data in ME/CFS reports regarding patient characteristics, sampling, statistical methods, clinical and research assessments, as we intended to do in this report. Moreover, the use of more specific and restrictive case definitions represents a pivotal strategy to reduce misclassification bias 78 and, as a consequence, improves the chances to identify reliable biomarkers in ME/CFS 40 .
For instance, despite the Fukuda case definition is the most used in ME/CFS research, it may misclassify individuals with major depression as ME/CFS patients 79 . In summary, including standardized diagnosis and research criteria into the study design will increase reproducibility/comparability among studies and positively impact ME/CFS patients' health care. www.nature.com/scientificreports/