MicroRNA Expression Profiling in Behçet's Disease

Background Behçet's disease (BD) is a chronic inflammatory multisystem disease characterized by oral and genital ulcers, uveitis, and skin lesions. MicroRNAs (miRNAs) are key regulators of immune responses. Differential expression of miRNAs has been reported in several inflammatory autoimmune diseases; however, their role in BD is not fully elucidated. We aimed to identify miRNA expression signatures associated with BD and to investigate their potential implication in the disease pathogenesis. Methods miRNA microarray analysis was performed in blood cells of BD patients and healthy controls. miRNA expression profiles were analyzed using Affymetrix arrays with a comprehensive coverage of miRNA sequences. Pathway analyses were performed, and the global miRNA profiling was combined with transcriptoma data in BD. Deregulation of selected miRNAs was validated by real-time PCR. Results We identified specific miRNA signatures associated with BD patients with active disease. These miRNAs target pathways relevant in BD, such as TNF, IFN gamma, and VEGF-VEGFR signaling cascades. Network analysis revealed several miRNAs regulating highly connected genes within the BD transcriptoma. Conclusions The combined analysis of deregulated miRNAs and BD transcriptome sheds light on some epigenetic aspects of BD identifying specific miRNAs, which may represent promising candidates as biomarkers and/or for the design of novel therapeutic strategies in BD.


Introduction
Behçet's disease (BD) is a rare and chronic multisystem disease characterized by a triple-symptom complex of recurrent oral aphthous ulcers, genital ulcers, and uveitis. Moreover, manifestations of vascular, articular, neurologic, urogenital, gastrointestinal, pulmonary, and cardiac involvement may occur. Hippocrates described BD in the fifth century BCE. In 1930, the Greek ophthalmologist Benediktos Adamantiades reported a patient with inflammatory arthritis, oral and genital ulcers, phlebitis, and iritis. The disease is named after the Turkish dermatologist Hulusi Behçet, who identified it in a patient in 1924 and published a description of the disease in 1937.
As virtually no unique histological or laboratory features have been identified to help in the diagnosis of the disease, clinical features are used to define and diagnose Behçet syndrome. An international study group on Behçet's disease has recently revised the criteria for the classification/ diagnosis of BD [1].
There are sporadic cases of BD all around the world, but it is most frequently seen along the ancient Silk Route because of its frequency in the Middle East and far-east Asia (prevalence of 14-20/100,000 inhabitants), and these regions have traditionally been considered the endemic areas for the condition [2].
BD is a sporadic disease, but a familial aggregation is well known. Carriers of HLA-B51/HLA-B5 have an increased risk of developing Behçet's disease compared with noncarriers. HLA-B51 is the strongest associated genetic factor, and it has been shown to be more prevalent in Turkish, Middle Eastern, and Japanese populations, with a higher prevalence of Behçet's disease in these populations [3]. Non-HLA genes also contribute to the development of BD [3]. Genome-wide association studies have shown that polymorphisms in genes encoding for cytokines, activator factors, and chemokines are associated with increased BD susceptibility. Among cytokines, IL-10 polymorphisms cause a reduction in the serum level of IL-10, an inhibitory cytokine that regulates innate and adaptive immune responses; on the other hand, IL-23 receptor polymorphism, which reduces the response to IL-23 stimulation, is associated with protection from BD [3][4][5]. Recent data reported also associations with CCR1, STAT4, and KLRC4 encoding for a chemokine receptor, a transcription factor implicated in IL-12 and IL-23 signaling and a natural killer receptor [6,7]. Finally, susceptibility genes involved in the innate immune response to microbial exposure have recently been identified by Immunochip analysis [8].
Increased Th1, CD4 + and CD8 + T cell, γδ + T cell, and neutrophil activities have been found both in the serum and in inflamed tissues of BD patients, suggesting the involvement of innate and adaptive immunity in the pathogenesis of BD [2,9]. Studies on T lymphocytes have suggested a Th1-predominant response. Both CD4 + and CD8 + lymphocytes are higher in the peripheral blood, with increased levels of IL-2 and interferon-(INF-) γ cytokines [10]. The cytokine Th17 may also play an important role in the pathogenesis of the disease [2,11]. This hypothesis is supported by the observation of high IL-21 and IL-17 levels in sera of patients affected by BD with neurologic involvement [12,13]. Another study has reported a higher Th17/Th1 ratio in peripheral blood of patients with BD compared to healthy controls, and this ratio was higher in patients with uveitis or folliculitis compared with patients without these manifestations [14,15].
MicroRNAs (miRNAs) are small noncoding RNAs that play an important role in the regulation of several biological processes through their interaction with cellular messenger RNAs [16]. Inflammatory responses have an impact on miRNA expression, regulating their biogenesis by altering the transcription and processing of precursor transcripts or influencing the stabilization of mature miRNAs [17,18]. In recent years, the number of miRNAs implicated in immune system development and function has dramatically enhanced, and there has been widespread discussion of their potential use as therapeutics for immunological diseases [16]. Indeed, the aberrant expression of miRNAs frequently occurs in human diseases, including hematological disorders and autoimmunity [19,20]. The concept that miRNAs participate in the pathogenesis of diseases, especially refractory diseases with unidentified mechanisms, might lead to a novel effective treatment. A number of studies have reported a differential expression of miRNAs in several inflammatory autoimmune diseases, such as in rheumatoid arthritis (RA), multiple sclerosis, systemic lupus erythematosus, psoriasis, and systemic sclerosis [21]. These studies highlighted a deep implication of miRNAs as regulatory molecules in autoimmunity and the intriguing possibility to use miRNAs as disease biomarkers in these immunological disorders.
As far as BD concerns, little is known about miRNA expression; moreover, no high-throughput miRNA expression studies have been conducted to identify miRNAs specifically associated with the disease and no study has been so far performed which combines the analysis of blood microRNAs with transcriptional profiles in patients with BD.
In the present study, we performed a miRNA microarray analysis on peripheral blood mononuclear cells (PBMCs) of BD patients. We identified specific miRNA signatures associated with patients with active BD. These deregulated miR-NAs target signaling pathways typically implicated in BD pathogenesis, such as TNF, interferon gamma, and VEGF and VEGFR signaling cascades.
The modular analysis of differentially expressed genes in BD revealed pathogenetically relevant networks that are possibly targeted by the identified miRNAs. This study sheds light on some aspects of BD pathogenesis identifying deregulated miRNAs as promising candidates for the discovery of disease biomarkers and/or as molecular tools for designing novel therapeutic strategies in BD.

Patients.
A group of 6 subjects with BD was used for the gene array study. All the patients attended the Unit of Autoimmune Diseases at the University Hospital in Verona, Italy.
All patients fulfilled the International Criteria for Behçet Disease (ICBD): oral aphthosis, genital ulcers, and ocular lesions were each given 2 points, whereas 1 point was assigned to each of skin lesions, vascular manifestations, and neurological manifestations. A patient scoring 4 points or above was classified as having BD [22,23]. At enrollment, none of the patients had active infections or was affected by malignancies.
The clinical features of the patients are reported in Table 1 that also includes a description of the BD patients selected for the gene array study. A written informed consent was obtained from all the participants to the study. The study was approved by the local Ethical Committee of the Azienda Ospedaliera Universitaria of Verona, Verona, Italy. All investigations have been conducted according to the principles expressed in the Helsinki declaration.

Microarray
Analysis. Blood samples were collected in BD Vacutainer K 2 EDTA tubes using a 21-gauge needle. Peripheral blood mononuclear cells (PBMC) were obtained upon stratification on Lympholyte® cell separation density gradient (Cedarlane, Burlington, Canada). Total RNA extraction from PBMC was performed with miRNeasy Mini Kit following the manufacturer's protocol (Qiagen GmbH, Hilden, Germany). RNA concentration and purity were evaluated by spectrophotometric analysis (NanoDrop 2000; Thermo Fisher Scientific, Wilmington, DE, USA), and a further check for RNA integrity was performed with 2100 Bioanalyzer (Agilent Genomics, Santa Clara, CA, USA) before microarray hybridization. Sample hybridization and scanning were performed as recommended by the Affymetrix (Affymetrix; Thermo Fisher Scientific Inc.) supplied protocols, by the Cogentech Affymetrix Microarray Unit (Campus IFOM-IEO, Milan, Italy), and Affymetrix GeneChip® miRNA 4.0 (Affymetrix; Thermo Fisher Scientific Inc., Waltham, MA, USA) was used. The GeneChip miRNA 4.0 arrays contain more than 30,000 probes including those encoding for 2578 mature human miRNAs, according to Sanger miRBase v.20.
The arrays were analyzed employing the Transcriptome Analysis Console (TAC) 4.0 software (Applied Biosystems, Foster City, CA USA, by Thermo Fisher Scientific, Waltham, MA, USA). The Signal Space Transformation-(SST-) Robust Multiarray Average (RMA) algorithm was applied to background-adjust, normalize, and log-transform signal intensity.
Relative expression levels of each microRNA were validated applying a one-way analysis of variance (ANOVA) and false discovery rate (FDR) correction (p ≤ 0 01). Micro-RNAs with an expression level of at least 1.5-fold different in the test sample versus control sample were analyzed.
Targeted genes of significantly modulated miRNAs were identified using the integrative database for human microRNA target predictions mirDIP [24]. All the source filters and a very high confidence class were applied for our analyses.
Pathway enrichment analysis of miRNA gene targets and differentially expressed genes (DEGs) in Behçet's disease was carried out using FunRich (http://www.funrich.org/) [25], and only Bonferroni-corrected enriched p values ≤0.01 calculated by the hypergeometric test were considered.
Pathways enrichment analysis of DEGs in BD was also performed employing the Panther expression analysis tools (http://pantherdb.org/) [ 2.4. Real-Time PCR. Mature miRNA expression was assayed by TaqMan® Advanced miRNA assay chemistry (Applied Biosystems, Foster City, CA, USA). Briefly, 10 ng of total RNA was reverse transcribed and preamplified with TaqMan Advanced miRNA cDNA synthesis kit following the manufacturer's instructions (Applied Biosystems, Foster City, CA, USA). Preamplified cDNA was diluted 1/10 in nucleasefree water, and 5 μL of diluted cDNA for each replicate was loaded in PCR. 20 μL PCR reactions were composed by 2x Fast Advanced Master Mix and TaqMan Advanced miRNA assays for hsa-miR-143-3p (477912_mir), hsa-miR-199a-5p (478231_mir), and hsa-miR-4505 (477842_mir). The mean of Ct for hsa-miR-16-5p (477860_mir) and hsa-miR-26a-5p (477995_mir) expression was used to normalize miRNA expression. Real-time PCR was carried out in triplicate on a QuantStudio 6 Flex instrument (Applied Biosystems, Foster City, CA, USA). Expression values were reported as fold change with respect to healthy controls by the ΔΔCt method using QuantStudio Real-Time PCR system software versus 1.3.

High-Throughput MicroRNA Expression Profiling in
Peripheral Blood Mononuclear Cells of Behçet's Disease. Since a global miRNA expression analysis in BD with an updated coverage of miRNA sequences has not been performed yet, we wanted to provide a careful description of miRNAs associated with BD by interrogating the transcription of a large amount of different microRNA sequences in BD PBMCs by microarray strategy.
Therefore, PBMCs derived from 6 patients with BD and 6 healthy age-and sex-matched donors were analyzed using a dedicated and high-density array with a coverage for more than 2500 human microRNA transcripts and all mature miRNA sequences in miRBase Release 20. The clinical features of patients included in the microarray study are reported in Table 1.
Microarray analysis revealed a high number (269) of modulated miRNAs that satisfied the FDR-corrected p value criterion (p ≤ 0 01) and the fold change criterion (FC ≥ |1 5|), showing a robust and statistically significant variation between BD and healthy control samples (Supplementary Table 1). Such a large number of modulated transcripts clearly reflect the high performance of the array in the detection of a wide range of microRNA sequences. We thereafter sharpened our analysis by selecting only modulated miRNAs annotated as "high confidence" in miRBase 21 (http://www.mirbase.org), and to make our results more informative, we further narrowed down the analysis to miRNAs for which gene targets were annotated in FunRich.
By these criteria, we selected 47 modulated miRNAs that are shown in Table 2. Interestingly, almost all (45/ 47) miRNAs were down-modulated with only two up-regulated microRNAs.
Selected miRNAs significantly deregulated in the microarray analysis were validated by real-time PCR in the entire series of patients analyzed (see Figure 1).

Pathway Enrichment Analysis of miRNAs Deregulated in
Behçet's Disease. In a second part of our analysis, we wanted to identify all the molecular pathways that were targeted by the selected miRNA performing a pathway enrichment analysis based on annotated gene targets in FunRich. The software allows to evaluate the miRNA regulatory effect and to identify controlled pathways based on predicted and/or validated miRNA-target interactions. We also applied a pathway enrichment analysis on the dataset of differentially expressed genes (DEGs) obtained in our previous study of gene expression profiling in BD (Puccetti et al., unpublished observations). In this study, we were able to select modulated genes that may play an important role in BD pathogenesis since they are involved in biological processes strongly connected to the typical features of the disease.
Despite the strong statistical stringency applied to the two datasets (p ≤ 0 01), we obtained a high number of significantly enriched pathways both from the miRNA target datasets (277) and from the BD DEGs (164) (Supplementary Tables 2 and 3). Notably, we found that a large proportion of pathways from the BD DEGs was also enriched in the list of miRNA-validated targets (64%, 106/ 164; Supplementary Table 4), indicating that the selected miRNAs exert a strong impact on the molecular pathways altered in the disease. Moreover, the large number of enriched pathways clearly reflected the multisystemic involvement typical of Behçet's disease. Figure 2 shows a graphical representation of selected commonly enriched pathways. Interestingly, the enriched categories were involved in vascular biology (i.e., glypican pathway, vascular endothelial growth factor, VEGF and VEGFR network, endothelin pathway, PAR1-mediated thrombin signaling events, thrombin/protease-activated receptor (PAR) pathway, EGFR-dependent endothelin signaling events, plateletderived growth factor (PDGF) receptor signaling network, and urokinase-type plasminogen activator (uPA), and uPAR-mediated signaling) and in apoptosis (TRAIL, p53, and FAS signaling pathways). In addition, other relevant pathways enriched in the two datasets were related to the immune response (i.e., IL6-mediated signaling events, TCR and BCR signaling, calcineurin-regulated NFAT-dependent transcription in lymphocytes, and toll-like receptor cascades) and to the inflammatory response (i.e., tumour necrosis factor, TNF receptor, IFN-gamma, p38 MAP kinase, pathway, and IL1-and CXCR4-mediated signaling events.

Comparative Analysis of Selected miRNA Gene Targets and Differentially Expressed Genes in BD.
To better define the role played by miRNAs in BD pathogenesis, we wanted to select miRNAs that were able to target genes modulated in BD. Therefore, we used a more sophisticated integrative database for human microRNA target predictions (mirDIP) [24] to obtain the lists of genes that were targeted, with a very high score class, by each of the selected miRNAs. Then, we compared the resulting target lists to genes differentially expressed in BD from our previous study (Puccetti et al., unpublished observations) and we observed that 65% of DEGs were targeted by the selected miRNAs. Interestingly, the vast majority of these DEGs showed an opposite modulation with respect to the relative targeting miRNA (Figure 3(a)), consistently with the typical role of miRNAs as negative regulator of gene expression (i.e., typically, upregulated genes are targeted by down-modulated miRNAs and underexpressed genes are targeted by overexpressed miRNAs). All the above-mentioned miRNAs are presented in Figure 3(b) along with the compiled lists of their targets DEGs.
All these selected miRNAs targeted genes involved in biological processes implicated in the disease pathogenesis including apoptosis, immune response, inflammation, and vascular damage ( Figure 4). Thus, we could identify micro-RNAs that may control the gene modulation involved in the disease pathogenesis. Table 3 shows all the targeted genes and their corresponding targeting miRNAs.

Network Analysis of Genes Differentially Expressed in BD.
We performed a network analysis in which the functional interactions between the protein products of modulated genes in BD were evaluated. By this approach, a proteinprotein interaction (PPI) network comprising 171 genes (nodes) and 3272 pairs of interactions (edges) was constructed ( Figure 5(a)). We then performed a clustering analysis to identify areas of densely interconnected nodes (clusters/modules; CL) that are predicted to be involved in common biological processes and to have a crucial role in the disease pathogenesis. We could detect six clusters that are graphically represented in Figures 5(b)-5(g). Comparing the list of miRNA targets to DEGs included in the six clusters (Supplementary Table 5), we found that, in each cluster, a large number of DEGs were targeted (Supplementary Table 6 and Figures 5(b)-5(g)). In particular, we observed that several of these genes were involved in immune and inflammatory responses including TNF (CL2), IL10 (CL2), and IL1A (CL1). Moreover, many of such genes played a role in B cell response, like EGR1 (CL1), or in T cell response like, for example, CTLA4, CD4, and MLL (CL1); DUSP2 (CL3); and CD28 (CL4) and NR4A2 (CL5). Other DEGs targeted in clusters were involved in TLR signaling including HSPA4, IKBKB, JUN, REL, S100A8, TLR2 (CL1), and ARF3 and NAMPT (CL2). Interestingly, in the six clusters, we also observed that various genes involved in angiogenesis and/or in vascular damage were targeted by selected miRNAs including THBS1, MMP8, VEGFA (CL1), FOSL2, THBD (CL2), HIPK1 and DDX6 (CL3), PGK1 (CL4), and ACTR2 (CL6). Given the well-known biological significance of highly connected gene clusters, to gain insight on the most relevant signaling networks that were controlled by the deregulated miRNAs, we performed a pathway enrichment analysis on targeted genes that were present in the six clusters. All the above-mentioned enriched pathways are listed in the Supplementary Table 7, and Table 4 shows a selection of the most relevant enriched signaling networks.

Discussion
A systematic analysis of miRNA expression profiles in BD has not been performed yet. The aim of our work is therefore to provide a compiled description of miRNAs associated with BD using an array able to query a very large number of transcripts in the attempt to dissect the possible regulatory effects exerted by these molecules on the molecular pathways relevant for the disease pathogenesis.
In a previous work, we investigated BD-associated transcriptional profiles by a gene expression analysis of PBMC derived from BD patients and identified a gene modulation strictly connected to the disease pathogenesis (Puccetti et al., unpublished observations). Moreover, we showed the presence of a type I interferon and Th-17 gene signature, which suggests an autoimmune origin of BD. In this study, we aimed to complement this gene expression analysis detecting modulated miRNAs that may target differentially expressed genes (DEGs) identified in our previous analysis.
To this purpose, we used a sophisticated target prediction system (mirDIP) [24] to obtain the list of miRNA targets, which was then compared to the list of DEGs identified in blood samples of BD patients. The comparison of the two analyses (gene expression and miRNAs) allows selecting within deregulated miRNAs only those related to the modulation of genes that would be effectively altered in the PBMCs of BD patients and might have therefore a pivotal role in the pathogenesis of the disease.
We observed that there was a good overlap (60%) between the selected miRNA targets and the BD-modulated genes, indicating that the majority of the identified miRNAs regulate genes differentially expressed in BD. Such target genes belonged to functional classes strictly connected to typical features of BD. Indeed, several miRNAs targeted proinflammatory genes such as TNF and IL1A. Moreover, transcripts involved in both adaptive and innate immune response were also targeted. In this regard, we have to mention that selected down-modulated miRNAs controlled several Th17 cell-associated genes and transcripts involved in type I interferon response that we found up-regulated in our previous analysis. This may suggest a loss of control in BD on two synergistic mechanisms typically associated with an autoimmune response. We moreover observed a downregulation of miRNAs that control members of TLRs and JAK/STAT pathways, two molecular signalings involved in autoimmune diseases [44,59] that are also active in BD [32,60,61]. In addition, miRNAs target genes involved in angiogenesis, and in blood coagulation, two processes commonly associated with BD vasculitis were also down-modulated in BD samples.
Consistently with the tight correspondence between deregulated miRNA and DEGS in BD, we observed a good overlap (64%) between the pathways enriched in the deregulated miRNA targets and in the genes differently expressed in BD. This indicated that, globally, the pathways identified by the selected miRNAs reflected fairly well the gene regulation described in our gene expression study, suggesting that the criteria applied for the selection of miRNAs had effectively identified miRNAs able to explain the gene modulation which we had previously described.
Interestingly, we found that both in BD-miRNA and in BD-DEGs dataset, meaningful signaling networks including,  for example, VEGF and VEGFR, endothelins, PDGF receptor, TNF receptor, and IL1-and IL6-mediated TCR and BCR pathways were enriched.
A pivotal task in molecular biology is to understand gene modulation in the context of biological networks. Indeed, proteins achieve their functions in the protein-proteininteracting network, and interestingly, the dynamics of such networks can be influenced by miRNAs [62]. We therefore wanted to identify significant relationships among modulated miRNAs and the PPI network in which the protein product of modulated genes in BD can be involved. Moreover, since it is known that the repressive effect of a miRNA may lead to more severe biological effects when it is exerted on proteins with more interacting partners [62], to highlight more efficient interrelations, we focused our attention on the most connected proteins of the network, inspecting the presence of miRNA targets inside the six clusters extracted from the PPI network. We found that the six clusters identified were extensively targeted by several modulated miRNAs, thus indicating that the deregulation of selected miRNAs may have a meaningful effect on the dynamics of a protein network that control the disease pathogenesis. Indeed, genes targeted in the six clusters played an important role in vascular biology (i.e., VEGFA, THBS, THBS1, etc.), in inflammation (i.e., TNF, IL1A, IL6R, CXCR4, etc.), and in immune response (i.e., CD28, CTLA4, EGR1, TLR2, etc.). Moreover, the analysis of pathways that were enriched in the target genes present in the clusters confirmed the essential role of these transcripts and their relative targeting miRNA in BD pathogenesis.
In conclusion, this work represents the first analysis performed on such a large number of miRNAs and integrated with the study of the profiles of gene expression in BD. The study allowed correlating the expression of miRNAs and the modulation of genes important for the pathogenesis of the disease. Using this approach, we have been able to identify the specific molecular pathways on which the regulation of these miRNAs may occur.
This study sheds light on some epigenetic aspects of BD identifying specific miRNAs, which may represent promising candidates for the identification of disease biomarkers and/or the design of novel therapeutic strategies in BD.

Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
The authors declare that they have no conflicts of interest.