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.

1. 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-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 miRNAs 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.

2. Materials and Methods

2.1. 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.

Table 1

Patients utilized for gene array study	6 (100%)
Sex
Male	4
Female	2
Clinical features
Aphthous stomatitis	6 (100%)
Genital ulcers	2 (33%)
Erythema nodosum-like lesions	1 (17%)
Papulopustular lesion	5 (83%)
Uveitis	2 (33%)
Epididymitis	0
Neurological symptoms	0
Vasculitis	4 (67%)
Joint manifestation	3 (50%)
Gastrointestinal involvement	0
Association with HLA-B51	4 (67%)

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.

2.2. Microarray Analysis

Blood samples were collected in BD Vacutainer K2EDTA 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). MicroRNAs 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/) [26].

2.3. Protein-Protein Interaction (PPI) Network Construction and Network Clustering

Differentially expressed genes (DEGs) in BD samples from our previous analyses (Puccetti et al. unpublished observations) were mapped to the STRING database (version 1.0; http://string-db.org/) [27] to detect protein-protein interactions (PPI) pairs validated by experimental studies. A score of ≥0.7 for each PPI pair was selected to construct the PPI network. For the topological analysis of the built network, Cytoscape software was used and network clustering analysis was performed with the MCODE plugin of Cytoscape, based on the thresholds of module score > 1.5 [28].

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 nuclease-free 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.

3. Results

3.1. 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.

Table 2

miRNAs modulated in BD patients versus healthy subjects.

ID	Probe set name	Transcript ID (array design)	Fold change	FDR p value	Accession
20518901	MIMAT0019041_st	hsa-miR-4505	9.46	0.0001	MIMAT0019041
20500781	MIMAT0004609_st	hsa-miR-149-3p	6.05	0.0003	MIMAT0004609
20500119	MIMAT0000065_st	hsa-let-7d-5p	−1.55	0.0066	MIMAT0000065
20500444	MIMAT0000256_st	hsa-miR-181a-5p	−1.76	0.0095	MIMAT0000256
20500778	MIMAT0000449_st	hsa-miR-146a-5p	−1.87	0.0032	MIMAT0000449
20501197	MIMAT0000703_st	hsa-miR-361-5p	−1.88	0.0013	MIMAT0000703
20503908	MIMAT0004780_st	hsa-miR-532-3p	−1.92	0.0032	MIMAT0004780
20502123	MIMAT0004748_st	hsa-miR-423-5p	−2	0.0052	MIMAT0004748
20501036	MIMAT0000617_st	hsa-miR-200c-3p	−2.12	0.0032	MIMAT0000617
20501182	MIMAT0000692_st	hsa-miR-30e-5p	−2.17	0.0081	MIMAT0000692
20500158	MIMAT0000085_st	hsa-miR-28-5p	−2.2	0.0016	MIMAT0000085
20500422	MIMAT0000244_st	hsa-miR-30c-5p	−2.49	0.0019	MIMAT0000244
20501276	MIMAT0000751_st	hsa-miR-330-3p	−2.53	0.0089	MIMAT0000751
20500797	MIMAT0000460_st	hsa-miR-194-5p	−2.6	0.0028	MIMAT0000460
20502124	MIMAT0001340_st	hsa-miR-423-3p	−2.61	0.0025	MIMAT0001340
20500159	MIMAT0004502_st	hsa-miR-28-3p	−2.66	0.0012	MIMAT0004502
20500718	MIMAT0000417_st	hsa-miR-15b-5p	−2.94	0.0005	MIMAT0000417
20500424	MIMAT0000245_st	hsa-miR-30d-5p	−2.99	0.0003	MIMAT0000245
20500795	MIMAT0004614_st	hsa-miR-193a-5p	−3.01	0.0015	MIMAT0004614
20500385	MIMAT0000222_st	hsa-miR-192-5p	−3.06	0.0019	MIMAT0000222
20500758	MIMAT0000438_st	hsa-miR-152-3p	−3.06	0.0046	MIMAT0000438
20500151	MIMAT0000081_st	hsa-miR-25-3p	−3.2	0.0003	MIMAT0000081
20503811	MIMAT0002821_st	hsa-miR-181d-5p	−3.3	0.0005	MIMAT0002821
20500123	MIMAT0000067_st	hsa-let-7f-5p	−3.34	0.0041	MIMAT0000067
20504274	MIMAT0003218_st	hsa-miR-92b-3p	−3.43	0.0018	MIMAT0003218
20500162	MIMAT0000087_st	hsa-miR-30a-5p	−3.74	0.0006	MIMAT0000087
20500488	MIMAT0000280_st	hsa-miR-223-3p	−4.32	0.0043	MIMAT0000280
20503887	MIMAT0002876_st	hsa-miR-505-3p	−4.58	0.0023	MIMAT0002876
20500733	MIMAT0000424_st	hsa-miR-128-3p	−5.07	0.0029	MIMAT0000424
20501291	MIMAT0000759_st	hsa-miR-148b-3p	−5.43	0.0084	MIMAT0000759
20501278	MIMAT0000752_st	hsa-miR-328-3p	−5.63	0.0029	MIMAT0000752
20500798	MIMAT0000461_st	hsa-miR-195-5p	−5.8	0.0012	MIMAT0000461
20500121	MIMAT0000066_st	hsa-let-7e-5p	−5.98	0.0052	MIMAT0000066
20500187	MIMAT0004514_st	hsa-miR-29b-1-5p	−6.81	0.0093	MIMAT0004514
20504378	MIMAT0003297_st	hsa-miR-628-3p	−6.83	0.0032	MIMAT0003297
20500170	MIMAT0004507_st	hsa-miR-92a-1-5p	−6.84	0.0005	MIMAT0004507
20500723	MIMAT0000419_st	hsa-miR-27b-3p	−7.22	0.0041	MIMAT0000419
20504553	MIMAT0004819_st	hsa-miR-671-3p	−8.16	0.0006	MIMAT0004819
20501287	MIMAT0000757_st	hsa-miR-151a-3p	−8.47	0.0084	MIMAT0000757
20503105	MIMAT0002177_st	hsa-miR-486-5p	−17.98	0.0066	MIMAT0002177
20500400	MIMAT0000232_st	hsa-miR-199a-3p	−21.75	0.0049	MIMAT0000232
20500458	MIMAT0004563_st	hsa-miR-199b-3p	−21.75	0.0049	MIMAT0004563
20500769	MIMAT0000445_st	hsa-miR-126-3p	−26.32	0.0023	MIMAT0000445
20504312	MIMAT0003249_st	hsa-miR-584-5p	−51.18	0.002	MIMAT0003249
20500399	MIMAT0000231_st	hsa-miR-199a-5p	−63.55	0.0023	MIMAT0000231
20500432	MIMAT0000250_st	hsa-miR-139-5p	−66.67	0.0027	MIMAT0000250
20500752	MIMAT0000435_st	hsa-miR-143-3p	−83.47	0.0046	MIMAT0000435

Selected miRNAs significantly deregulated in the microarray analysis were validated by real-time PCR in the entire series of patients analyzed (see Figure 1).

Figure 1

Validation of differentially expressed miRNAs by real-time PCR. Real-time PCR for the indicated miRNAs were performed in healthy controls (Healthy) and in BD patients (Behçet). Values are calculated as fold change with respect to healthy samples with the ΔΔCt method. miR-16-5p and miR-26-5p were used as endogenous controls for miRNA expression (see Material and Methods). Histograms indicate mean values; bars indicate standard deviation (SD). p = 0.05 (Mann–Whitney test).

3.2. 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, platelet-derived 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.

Figure 2

Histogram representing pathways enriched in BD-modulated miRNA targets genes and in BD DEGs. y axis: −log10(p value) (hypergeometric test).

3.3. 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, up-regulated 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.

Figure 3

Selection of miRNAs that targeted genes differentially expressed in BD. (a) Pie chart showing the percentage of BD DEGs that were targeted by miRNAs modulated in BD. (b) BD-modulated miRNAs and their respective targeted BD DEGs. Genes associated with Th17 cells are written in bold characters.

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 microRNAs that may control the gene modulation involved in the disease pathogenesis. Table 3 shows all the targeted genes and their corresponding targeting miRNAs.

Figure 4

Functional classification of BD DEGs targeted by selected miRNAs. Pie charts showing the different GO biological processes in which BD DEGs targeted by selected deregulated miRNAs in BD can be classified.

Table 3

List of all BD DEGs that were targeted by selected miRNAs modulated in BD and their respective targeting miRNAs.

DEGs	MicroRNA
ACTR2	miR-181a-5p	miR-151-3p	miR-199a-3p	miR-199b-3p	miR-15b-5p	miR-195-5p
AKIRIN2	miR-181a-5p	miR-27b-3p	miR-139-5p
APLP2	miR-181a-5p	miR-139-5p	miR-199a-3p	miR-199b-3p
ARF3	miR-181a-5p	miR-27b-3p
BCL2L11	miR-181a-5p	miR-199a-3p	miR-199b-3p	miR-30e-5p	miR-30c-5p	miR-30d-5p	miR-181d-5p	miR-92b-3p	miR-30a-5p	miR-148b-3p
BTG2	miR-181a-5p	miR-27b-3p	miR-199a-5p	let-7d-5p	miR-15b-5p	miR-25-3p	let-7f-5p	miR-92b-3p	miR-195-5p	let-7e-5p	miR-146a-5p
CALM1	miR-181a-5p	miR-139-5p	miR-143-3p	miR-151-3p	miR-199a-3p	miR-199b-3p	let-7d-5p	miR-181d-5p	let-7f-5p	miR-195-5p	let-7e-5p
CD4	miR-181a-5p
CXCL2	miR-181a-5p
DAZAP2	miR-181a-5p	miR-27b-3p	miR-139-5p	miR-151-3p	miR-29b-1-5p	miR-146a-5p
DDX3X	miR-181a-5p	miR-27b-3p	miR-139-5p	miR-199a-5p	miR-671-3p	miR-195-5p	miR-146a-5p
DNAJB1	miR-181a-5p	miR-27b-3p	miR-143-3p	miR-486-5p
EGR1	miR-181a-5p	miR-143-3p	miR-199a-3p	miR-199b-3p	miR-146a-5p
ERAP1	miR-181a-5p
EREG	miR-181a-5p	miR-199a-3p	miR-199b-3p	miR-192-5p
F5	miR-181a-5p
HIPK1	miR-181a-5p	miR-139-5p	miR-143-3p	miR-199a-3p	miR-199a-5p	miR-199b-3p	miR-486-5p	miR-146a-5p
HSP90B1	miR-181a-5p	miR-152-3p	miR-181d-5p	miR-223-3p	miR-148b-3p
IL1A	miR-181a-5p	miR-181d-5p	miR-505-3p
IL6ST	miR-181a-5p	miR-27b-3p	miR-143-3p	miR-152-3p	miR-223-3p	miR-505-3p	miR-148b-3p
INPP4A	miR-181a-5p	miR-223-3p
KMT2A	miR-181a-5p	miR-27b-3p	miR-199a-3p	miR-199b-3p	miR-148b-3p	miR-195-5p	miR-146a-5p
LMAN1	miR-181a-5p	miR-152-3p	miR-148b-3p
MARCKS	miR-181a-5p	miR-27b-3p	miR-143-3p	miR-199a-5p	miR-584-5p	miR-200c-3p	miR-128-3p
MCL1	miR-181a-5p	miR-139-5p	miR-151-3p	miR-199a-3p	miR-199b-3p	miR-486-5p	miR-584-5p	miR-628-3p
NKTR	miR-181a-5p	miR-27b-3p	miR-151-3p
NR4A2	miR-181a-5p	miR-361-5p	miR-30e-5p	miR-30c-5p	miR-30d-5p	miR-30a-5p
NR4A3	miR-181a-5p	miR-25-3p	miR-92b-3p
PRRC2C	miR-181a-5p	miR-27b-3p	miR-486-5p	miR-195-5p
SOCS3	miR-181a-5p	miR-30e-5p	miR-30c-5p	miR-30d-5p	miR-152-3p	miR-30a-5p	miR-148b-3p
TNF	miR-181a-5p	miR-181d-5p
YY1	miR-181a-5p	let-7f-5p	let-7e-5p
ZFP36L2	miR-181a-5p	miR-27b-3p	miR-128-3p
TNFAIP3	miR-27b-3p	let-7e-5p
PLK2	miR-27b-3p	miR-126-3p	miR-128-3p	miR-146a-5p
MAP4K5	miR-27b-3p	miR-128-3p
HSPD1	miR-27b-3p	miR-143-3p	miR-628-3p	miR-505-3p
IKZF1	miR-27b-3p	miR-199a-5p
IL10	miR-27b-3p	let-7d-5p	let-7f-5p	let-7e-5p	miR-146a-5p
CD28	miR-27b-3p	miR-143-3p	miR-15b-5p	miR-195-5p
MLL	miR-27b-3p	miR-143-3p	miR-199a-3p	miR-199b-3p
PIK3R1	miR-27b-3p	miR-126-3p	miR-199a-3p	miR-199b-3p	miR-486-5p	miR-15b-5p	miR-128-3p	miR-195-5p
PTGER4	miR-27b-3p	miR-25-3p	miR-92b-3p	miR-146a-5p
STAT1	miR-27b-3p	miR-584-5p	miR-128-3p	miR-146a-5p
ICOS	miR-27b-3p	let-7d-5p	let-7f-5p	let-7e-5p
KAT6A	miR-27b-3p	miR-139-5p	miR-143-3p
VEGFA	miR-126-3p	miR-139-5p	miR-199a-3p	miR-199a-5p	miR-199b-3p	miR-486-5p	miR-361-5p	miR-15b-5p	miR-195-5p
CXCR4	miR-139-5p	miR-143-3p
JUN	miR-139-5p	miR-200c-3p
AP2M1	miR-139-5p	miR-143-3p	miR-423-5p
TCF3	miR-139-5p	miR-195-5p
CREM	miR-139-5p	miR-199a-5p
RBL2	miR-139-5p	miR-199a-3p	miR-628-3p
TGIF1	miR-139-5p	miR-27b-3p	miR-199a-3p	miR-199b-3p	miR-25-3p	miR-92b-3p	miR-146a-5p
THBD	miR-139-5p
THBS1	miR-139-5p	miR-143-3p	miR-151-3p	miR-199a-3p	miR-199b-3p	miR-628-3p	let-7d-5p	miR-194-5p	let-7f-5p	let-7e-5p
MMP8	miR-139-5p
DDX6	miR-139-5p	miR-143-3p	miR-199a-5p	miR-152-3p	miR-148b-3p
PGK1	miR-27b-3p	miR-143-3p	miR-146a-5p
TSC22D3	miR-143-3p	miR-195-5p
IL6R	miR-143-3p	let-7e-5p
DUSP1	miR-143-3p	miR-199a-3p	let-7d-5p	miR-200c-3p	miR-152-3p	miR-25-3p	let-7f-5p	miR-148b-3p	let-7e-5p
CD84	miR-143-3p
CTLA4	miR-143-3p
STAT2	miR-143-3p
TLR2	miR-143-3p	miR-146a-5p
PPP2R1A	miR-143-3p	miR-199a-5p	miR-15b-5p	miR-195-5p
SIKE1	miR-143-3p	miR-199b-3p	miR-199a-3p	miR-146a-5p
DUSP2	miR-143-3p
FLOT1	miR-143-3p
FOSL2	miR-143-3p	miR-30e-5p	miR-30c-5p	miR-30d-5p	miR-25-3p	miR-92b-3p	miR-30a-5p	miR-195-5p
IKZF3	miR-151-3p	miR-146a-5p
PTPN12	miR-151-3p	miR-486-5p	miR-194-5p
NAMPT	miR-199a-3p	miR-199b-3p
PDE4B	miR-199a-3p	miR-199b-3p	miR-152-3p	miR-148b-3p
STK4	miR-199a-5p	miR-486-5p
IKBKB	miR-199a-5p	miR-200c-3p	miR-15b-5p	miR-195-5p
HSPA4	miR-199a-5p
GNAI2	miR-29b-1-5p	miR-30e-5p	miR-30c-5p	miR-30d-5p	miR-30a-5p
ZNF331	miR-486-5p	miR-146a-5p
PTX3	miR-628-3p
CCL3	let-7d-5p	let-7f-5p	let-7e-5p
OLR1	let-7d-5p	let-7e-5p
S100A8	let-7d-5p	let-7f-5p	let-7e-5p
SOCS1	let-7d-5p	miR-30e-5p	miR-30c-5p	miR-30d-5p	let-7f-5p	miR-30a-5p	let-7e-5p
DGKA	miR-200c-3p
NAPA	miR-330-3p
ARHGDIA	miR-15b-5p	miR-195-5p
SNX27	miR-152-3p	miR-148b-3p
NLRP3	miR-223-3p
RHOB	miR-223-3p
FLOT2	miR-148b-3p	miR-146a-5p
RAB35	miR-148b-3p	miR-195-5p
REL	let-7e-5p	miR-146a-5p
KIR2DL4	miR-146a-5p

Selected miRNAs targeted many DEGs that sustained the inflammatory response typically associated with BD including TNF, IL1A, IL10, IL6R, CXCL2, CXCR4, TNFAIP3, OLR1, S100A8, HSP90B1, and CCL3 (see Figure 3(b) and Table 3). Interestingly, TNF (targeted by hsa-miR-181d-5p and -181a-5p), CCL3 (targeted by hsa-let-7e-5p, -7d-5p, and -7f-5p), IL10 (targeted by hsa-let-7f-5p, -7e-5p, -7d-5p, -miR-146a-5p, and -27b-3p), and IL1A (targeted by hsa-miR-505-3p, -181d-5p, and -181a-5p) have been detected at increased concentrations in sera or plasma of BD patients when compared to normal subjects [29-31]. Moreover, down-modulated miRNAs targeted DEGs involved in the adaptive immune response including genes that played a role in T and in B cell immune response. Among these genes, several miRNAs targeted CD28 (hsa-miR-143-3p, -27b-3p, -15b-5p, and -195-5p), CD4 (hsa-miR-181a-5p), ICOS (hsa-miR 27b-3p and -let-7f-5p, -7e-5p, and -7d-5p), CTLA4 (hsa-miR-143-3p), EGR1 (hsa-miR-143-3p, -146a-5p, -181a-5p, -199b-3p, and -199a-3p), and AKIRIN2 (hsa-miR-139-5p, -27b-3p, and -181a-5p). Interestingly, 18 miRNAs targeted genes of the TH-17 gene signature including CD28, CD4, ICOS, SOCS3, IL6ST, YY1, IL6R, TNF, CXCL2, and SOCS1 (Figure 3(b)).

In addition, selected miRNAs could also control DEGs involved in the innate immune response including, for example, NKTR (hsa-miR-181a-5p, -486-5p, -151-3p, and -27b-3p), KIR2DL4 (hsa-miR-146a-5p), STAT1 (hsa-miR-584-5p, -27b-3p, -128-3p, and -146a-5p), STAT2 (hsa-miR-143-3p), and DEGs belonging to the Toll-like receptor (TLR) pathway such as TLR2, IKBKB, JUN, REL, NAMPT, HSP90B1, HSPA4, and ARF3. In particular, we have to mention that TLR2 (targeted by hsa-miR-143-3p and -146a-5p) is thought to be up-regulated in BD patients [32, 33]. As many as 22 miRNAs targeted genes of the type I interferon signature like VEGFA, DDX3X, and the above-mentioned STAT1, STAT2, SOCS1, and IL10. The coexpression of miRNAs targeting DEGs involved in the type I IFN signaling and Th-17 related genes may reflect the presence of a synergy between IFN and Th17 pathways that is typical of autoimmune diseases [34-43], thus suggesting that an autoimmune mechanism can be involved in the BD pathogenesis.

Noteworthily, several miRNAs targeted all the eight genes involved in the JAK/STAT signaling pathway that we found up-regulated in our previous work (Puccetti et al., unpublished observations) including PIK3R1 and the above-mentioned STAT1, STAT2, IL6ST, IL6R, IL10, SOCS1, and SOCS3. The downmodulation of miRNAs targeting this pathway further supports the hypothesis of an autoimmune origin of BD since its activation is very frequently associated with autoimmune diseases [44]. Moreover, the JAK/STAT signaling pathway is active in CD4+ T cells of patients with BD [45].

Interestingly, modulated miRNAs targeted DEGs involved in the vascular damage associated with BD that is characterized by myointimal proliferation, fibrosis, and thrombus formation [46]. Indeed, several miRNAs targeted DEGs associated with blood coagulation (i.e., THBS1, F5, and LMAN1; see Table 3), a process whose alteration is typically associated with BD vasculitis. Noteworthily, 10 down-modulated miRNAs (hsa-let-7d/7e/7f-5p, hsa-miR-151-3p, -628-3p, -139-5p, 143-3p, -194-5p, -199a-3p, and -199b-3p) targeted THBS1, suggesting a crucial role of THBS1 in the pathogenesis of vascular damage. Other modulated miRNAs targeted DEGs that played a role in angiogenesis including MMP8, VEGFA, NR4A3, and NAPA. Moreover, NR4A3 is a transcription factor involved in vascular development [47] and NAPA promotes vascular endothelial- (VE-) cadherin localization at endothelial junctions [48]. Interestingly, increased serum levels of soluble VE-cadherins have been detected in BD patients [49]. Among the above-mentioned genes, we have to notice that VEGFA was predictively targeted by a high number (9) of down-modulated miRNAs, including hsa-miR-199a-3p, -15b-5p, -195-5p, -361-5p, -126-3p, -486-5p, -199a-5p, 199b-3p, and hsa-miR-139-5p. Other DEGs involved in vascular damage and targeted by modulated miRNAs were FOSL2, THBD, and PTX3. THBD, targeted by hsa-miR-139-5p, is increased in sera of BD patients compared with healthy subjects, and PTX3, targeted by hsa-miR-628-3p, is considered as a marker of small vessel vasculitis [50].

Moreover, DEGs involved in the apoptotic process were targeted, and among these, we can cite MCL1, DNAJB1, BCL2L11, ARHGDIA, ZNF331, and IER3 (Table 3). Finally, several modulated miRNAs targeted DEGs that played a role in cell proliferation, a process that can be induced in response to both apoptosis and to skin ulcer formation. Among these miRNAs, we can mention those that targeted BTG2, EREG, PRRC2, and APLP2. In particular, EREG (targeted by hsa-miR-199a-3p, -199b-3p, and -181a-5p and by hsa-miR-192-5p) and APLP2 (targeted by hsa-miR-199a-3p, -199b-3p, -139-5p, and -181a-5p) were involved in the proliferation of corneal epithelial cells and in corneal epithelial wound healing, respectively [51, 52]. Thus, this gene regulation may have a role in ocular manifestations of BD like keratitis.

3.4. 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 protein-protein 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).

Figure 5

Network analysis of genes modulated in Behçet's disease. (a) PPI network of differentially expressed genes in BD. (b–g) Graphical representation of the six clusters that were extracted from the PPI network of modulated genes in BD. Red dots indicate DEGs that are targeted by modulated miRNAs in BD.

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.

Table 4

Pathways enriched in DEGs targeted in each cluster.

Panther pathways	p value
CL1
Toll-like receptor signaling pathway	8.44E − 07
Apoptosis signaling pathway	2.34E − 05
Inflammation-mediated by chemokine and cytokine signaling pathway	3.09E − 05
Interleukin signaling pathway	1.68E − 04
PDGF signaling pathway	1.14E − 03
Oxidative stress response	1.46E − 03
B cell activation	1.80E − 03
T cell activation	3.49E − 03
JAK/STAT signaling pathway	3.81E − 03
Angiogenesis	1.79E − 02
CL2
Wnt signaling pathway	1.20E − 02
Hypoxia response via HIF activation	1.87E − 02
Blood coagulation	2.03E − 02
Insulin/IGF pathway-protein kinase B signaling cascade	2.31E − 02
p53 pathway feedback loops 2	2.87E − 02
PI3 kinase pathway	3.15E − 02
VEGF signaling pathway	3.64E − 02
Endothelin signaling pathway	4.41E − 02
T cell activation	4.74E − 02
p53 pathway	4.90E − 02
CL3
5HT1-type receptor-mediated signaling pathway	2.07E − 02
PI3 kinase pathway	2.63E − 02
Oxidative stress response	2.72E − 02
CL4
Glycolysis	6.63E − 03
T cell activation	2.79E − 02
CL5
PDGF signaling pathway	2.71E − 02
Angiogenesis	3.03E − 02
Integrin signaling pathway	3.04E − 02
CL6
Cadherin signaling pathway	2.10E − 02

The TLR and the JAK/STAT signaling networks, two pathways notoriously involved in immune response and frequently associated with autoimmunity, were enriched in CL1-targeted genes (CL1-TGs). Other pathways implicated in the immune response were enriched in several genes targeted in clusters including T cell activation (CL1-TGs, CL2-TGs, and CL4-TGs), B cell activation (CL1-TGs), and PI3K signaling (CL2-TGs and CL3-TGs). This pathway is a crucial element in the regulation of adaptive and innate immune response [53, 54] and is a key player in inflammatory response. Therefore, in recent years, highly selective inhibitors of PI3K have been developed for anti-inflammatory treatments [55]. An enrichment in signalings involved in the inflammatory response was also found in CL1-TGs (inflammation mediated by chemokine and cytokine and interleukin pathways) and in CL5-TGs (integrin pathway).

In CL1-TGs and in CL3-TGs, the oxidative stress response pathway was enriched. Interestingly, an increased oxidative stress has been described in BD and it has been correlated to the severe inflammatory and degenerative clinical manifestations of the disease [56]. Several networks clearly connected to vascular damage were enriched in cluster TGs including PDGF, angiogenesis (CL1-TGs and CL5-TGs), Wnt, blood coagulation, endothelin, VEGF (CL2-TGs), and cadherin (CL6-TGs) signaling pathway. VE-cadherins have a crucial role in endothelial barrier integrity, and noteworthily, in CL6 the ACTR2 gene was targeted, known to stabilize adherens junctions between endothelial cells of the vascular wall interacting with cadherins [57]. Moreover, the Wnt signaling pathway is critically involved in vascular biology [58]. Finally, in CL1-TGs and in CL2-TGs the apoptosis and the p53 pathways were enriched, respectively.

4. 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-protein-interacting 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.