Behçet's: A Disease or a Syndrome? Answer from an Expression Profiling Study

Behçet’s disease (BD) is a chronic, relapsing, multisystemic inflammatory disorder with unanswered questions regarding its etiology/pathogenesis and classification. Distinct manifestation based subsets, pronounced geographical variations in expression, and discrepant immunological abnormalities raised the question whether Behçet’s is “a disease or a syndrome”. To answer the preceding question we aimed to display and compare the molecular mechanisms underlying distinct subsets of BD. For this purpose, the expression data of the gene expression profiling and association study on BD by Xavier et al (2013) was retrieved from GEO database and reanalysed by gene expression data analysis/visualization and bioinformatics enrichment tools. There were 15 BD patients (B) and 14 controls (C). Three subsets of BD patients were generated: MB (isolated mucocutaneous manifestations, n = 7), OB (ocular involvement, n = 4), and VB (large vein thrombosis, n = 4). Class comparison analyses yielded the following numbers of differentially expressed genes (DEGs); B vs C: 4, MB vs C: 5, OB vs C: 151, VB vs C: 274, MB vs OB: 215, MB vs VB: 760, OB vs VB: 984. Venn diagram analysis showed that there were no common DEGs in the intersection “MB vs C” ∩ “OB vs C” ∩ “VB vs C”. Cluster analyses successfully clustered distinct expressions of BD. During gene ontology term enrichment analyses, categories with relevance to IL-8 production (MB vs C) and immune response to microorganisms (OB vs C) were differentially enriched. Distinct subsets of BD display distinct expression profiles and different disease associated pathways. Based on these clear discrepancies, the designation as “Behçet’s syndrome” (BS) should be encouraged and future research should take into consideration the immunogenetic heterogeneity of BS subsets. Four gene groups, namely, negative regulators of inflammation (CD69, CLEC12A, CLEC12B, TNFAIP3), neutrophil granule proteins (LTF, OLFM4, AZU1, MMP8, DEFA4, CAMP), antigen processing and presentation proteins (CTSS, ERAP1), and regulators of immune response (LGALS2, BCL10, ITCH, CEACAM8, CD36, IL8, CCL4, EREG, NFKBIZ, CCR2, CD180, KLRC4, NFAT5) appear to be instrumental in BS immunopathogenesis.


Introduction
Behçet's disease (BD) is a multisystemic inflammatory disorder with a strong and complex genetic background [1]. Now, nearly 80 years after its initial description in 1937, many important questions regarding BD still remain unanswered, not only in relation to its etiology/pathogenesis but also to its classification [2]. Besides its significant morbidity profile, BD is reported to be a cause of increased mortality among the young male patients [3].
The hallmark of BD is its recurrent mucocutaneous lesions. Nevertheless, patients with BD also display a diverse spectrum of clinical manifestations including ocular, vascular, gastrointestinal, musculoskeletal, and central nervous systems [4]. The presence of well-defined clusters/subsets of BD patients with distinctly associated manifestations of the disease, marked regional variations in the expression of BD around the globe, and the proposal stating that distinct immunological abnormalities are underlying distinct classification groups of BD raised the question whether Behçet's is "a disease or a syndrome" [4][5][6]. At present, despite the massive amount of available data, this question is still not answered conclusively.
The introduction of microarray technology and its implementation for whole-genome expression analysis allowed scientists to study the differentially expressed genes (DEGs) in health and disease states at a genome-wide level [7]. With such a huge amount of highthroughput expression data in hand, bioinformatic and pathway analysis tools help researchers to delineate the pathways responsible for the development of diseases. Furthermore, the accumulation of research data in public repositories (i.e., databases open to public access) creates an opportunity for meta-analysis and data mining and thus analyzing data from different perspectives and condensing it into useful information.
With the purpose of answering the question of whether BD is a disease or a syndrome, we aimed to clarify and compare the molecular mechanisms underlying different expressions of BD. In this context, we used the expression data of the key gene expression profiling and association study on BD by Xavier et al [8]. The gene expression data provided by Xavier et al was retrieved from Gene Expression Omnibus (GEO) and reanalysed by implementation of gene expression data analysis/visualization and bioinformatics enrichment tools [9]. We obtained evidence of apparent expression profile discrepancies among BD patients with distinct expressions of the disorder. Furthermore, our findings supported the potential role of four gene groups (i.e., negative regulators of inflammation, neutrophil granule proteins, antigen processing and presentation proteins, and regulators of immune response) in BD immunopathogenesis.

Materials and Methods
The gene expression profiling study by Xavier et al Fifteen patients with BD all diagnosed according to the revised International Criteria and 14 healthy control subjects were enrolled in the study by Xavier et al [8,10]. Total RNA was isolated from peripheral blood mononuclear cells and GeneChip 1 Human Genome U133 Plus 2.0 (Affymetrix, Santa Clara, CA, USA) microarrays were used for hybridization [8]. According to the specifications in its product description, the GeneChip 1 Human Genome U133 Plus 2.0 array is a comprehensive whole human genome expression array which covers >47,000 transcripts for expression profiling. The study by Xavier et al was conducted and reported in accordance with the Minimum Information About Microarray Experiment (MIAME) guidelines and both the raw and the processed microarray data were deposited on GEO database with the Series ID GSE17114 [8,9,11,12].

Retrieval of the microarray data
The raw microarray data of the study by Xavier et al was retrieved from GEO by using the GEO accession GSE17114 on August 25 th 2015 [8,9,12]. The relevant file was a TAR file with the name "GSE17114_RAW" including 29 individually compressed CEL files (GSM428037. CEL.gz-GSM428065.CEL.gz).

Definition of subsets of patients with Behçet's disease
By using the principal clinical characteristics (major clinical symptoms) of BD patients briefly summarized in the article by Xavier et al, we subgrouped the BD patients according to their disease manifestations [8]. Three subsets were generated, namely, BD patients with mucocutaneous involvement (MB), BD patients with ocular involvement (OB), and BD patients with vascular involvement (VB). BD patients with isolated mucocutaneous manifestations (i.e., oral aphtosis, genital aphtosis, skin aphtosis, pseudofolliculitis, erythema nodosum, positive Pathergy test) were grouped as MB; BD patients with any kind of ocular involvement (i.e., anterior uveitis, posterior uveitis, retinal vasculitis) were grouped as OB, and BD patients demonstrating large vein thrombosis were grouped as VB. The group inclusive of all of the 15 BD patients was named as B, while the control group was given the name C.
Pre-processing of the microarray data Before obtaining the transcriptomic-level measurements and continuing with the downstream analysis, pre-processing of the microarray data was performed using BRB-ArrayTools v4.4.1 Stable Release developed by Dr. Richard Simon and BRB-ArrayTools Development Team [13]. The gene expression data present as raw CEL files was collated by the data import function of BRB-ArrayTools. The Robust Multiarray Average (RMA) algorithm, including background correction, log base 2 transformation, and quantile normalization was used for normalization of the microarray data [14]. Following normalization, the replicate spots within each individual array were averaged. Finally, gene filters were implemented which excluded genes if less than 20% of the genes' expression values had at least a 1.5 fold change in either direction from the genes' median expression values or if the genes' missing expression values exceeded 50%.

Verification of the manifestation based subgrouping of Behçet's disease patients
For the verification of our manifestation based subgrouping of BD patients, an initial cluster analysis was performed. BRB-ArrayTools' built-in clustering tools, Cluster 3.0 and TreeView softwares developed by Michael Eisen and the Stanford group were used for clustering [13,15]. A hierarchical clustering algorithm using Euclidean distance metric and average linkage was implemented and both the patients and the genes were clustered. The gene sets used for clustering were constituted from the DEGs identified during the class comparisons MB vs OB, MB vs VB, OB vs VB (two-sample t-test, p0.001, fold change 3 for MB vs VB and 4 for MB vs OB and OB vs VB), and MB vs OB vs VB (ANOVA, p0.001).

Analysis of the gene expression data
Class comparison analysis among BD and C groups were performed using BRB-ArrayTools [13]. For class comparison analysis of two classes, two-sample t-test with random variance model was implemented. DEGs were selected using a p-value 0.05 and a fold change (FC) 2. In only two cases, namely, for the class comparisons B vs C and MB vs C, an FC of 1.5 was also used. The Venn diagram representation of class comparisons was drawn with Venny 2.0.2 by Juan Carlos Oliveros [16].
The tools used (i.e., BRB-ArrayTools' built-in clustering tools, Cluster 3.0 and TreeView softwares) and the methodology implemented (i.e., hierarchical clustering using Euclidean distance metric and average linkage) for cluster analysis of the gene expression data were exactly the same as described above [13,15]. Similarly, patients and genes were clustered together and again, the gene sets used for clustering were constituted from the DEGs identified during the class comparisons MB vs OB, MB vs VB, OB vs VB (two-sample t-test, p0.001 and FC4 for all), and MB vs OB vs VB (ANOVA, p0.001).
Gene Ontology (GO) term enrichment analysis of the DEGs were performed with Web-Based Gene Set Analysis Toolkit (WebGestalt) and the enrichment analysis specifically focused on the sub-root of biological process (BP) [17,18]. For the enrichment analysis of GO terms, the DEGs retrieved by the class comparisons MB vs C (p0.05 and FC1.5), OB vs C and VB vs C (p0.05 and FC2 for both) were implemented and the setup included the hypergeometric test, the Benjamini-Hochberg procedure for multiple test adjustment, and 2 as the minimum number of genes for a category.
The flow diagram of the study is shown in Fig 1. Matching of the loci of the differentially expressed genes with the loci identified in the genome-wide association and the genome-wide linkage studies of Behçet's disease In order to document the matches between the loci identified in the genome-wide association (GWA) and the genome-wide linkage (GWL) studies of BD and the loci of the DEGs defined in the present study, the linkage study by Karasneh et al, and the two association studies by Meguro et al and Kirino et al were employed [19][20][21]. A total of 25 non-HLA loci were included and the DEGs identified during the class comparisons MB vs C (p0.05 and FC1.5), OB vs C and VB vs C (p0.05 and FC2 for both) were utilized.

Statistical analysis
Demographic data analysis of the study population was carried out using the SPSS 17 software (SPSS Statistics for Windows, Version 17.0, Chicago, SPSS Inc.). Ages were expressed as mean ±SD and gender as ratios (M/F). For comparing the means of two independent groups the Mann-Whitney U test was implemented, whereas comparison of the ratios of two independent groups was performed by the chi-square (χ2) test. p0.05 was considered to be statistically significant.

Demographic and clinical characteristics
Demographic and basic clinical characteristics of the study population is presented in Table 1.

Verification cluster analysis
During initial clustering, the sample with the experiment name VB1 (study ID 2, GEO accession number GSM428038) appeared to belong to a different BD subset (i.e., MB) instead of its originally assigned BD subset (i.e., VB).   this sample was excluded prior to further analysis and thus the number of samples in group B and VB dropped to 14 and 3 respectively. Except the case with VB1, the expression profiling based clustering results exactly matched the manifestation (clinical) based clustering of BD patients. The experiment descriptor file (EDF) of the study is given in the S1 File, while the gene lists used for initial clustering are presented in the S2 File. . For both cases, hierarchical clustering using Euclidean metric and average linkage was employed and both patients and genes were clustered. For sake of simplicity, only the dendrograms for clustering of the patients are shown in the heatmaps. Take note of the position of VB1 (study ID 2, GSM428038) in the MB branch of the dendrograms. Based on these clustering results, VB1 was excluded prior to further analysis. Also, as can be seen in the figure, the cluster analysis successfully clustered distinct expressions of Behçet's disease and with the exception of VB1, the expression profiling based clustering results were in accordance with the manifestation based clustering of Behçet's disease patients. The gene sets used for clustering were constituted from the DEGs identified during the corresponding class comparisons (i.e., MB vs VB and MB vs OB vs VB).

Class comparison analysis
The number of the DEGs found during class comparison analysis between BD and C groups are summarized in Table 2. The class comparison B vs C yielded a DEG number of 4 and when an FC of 1.5 was implemented, the number of the DEGs increased to 20. Similarly, the class comparison MB vs C also documented a very low number of DEGs (i.e., 5) which increased to 71 again with an FC of 1.5. Interestingly, class comparison analysis between the two other BD subsets and C (i.e., OB vs C and VB vs C) and also class comparison analysis of BD subsets among themselves yielded significantly higher numbers of DEGs ( Table 2). The number of the DEGs for class comparisons OB vs C, VB vs C, MB vs OB, MB vs VB, and OB vs VB were 151, 274, 215, 760, and 984 respectively. The gene lists of the DEGs are presented in the S3 File.
The Venn diagram representation of the class comparisons MB vs C, OB vs C, and VB vs C is shown in Fig 3. As can be seen in Fig 3, the number of the common DEGs in the intersection of MB vs C and OB vs C and VB vs C is zero.
The top 20 DEGs with respect to their FC values are listed in Table 3 with their gene symbols, probe set IDs, FC and p values. Worthy of note, while the absolute maximum and minimum FC values of the top 20 DEGs of the class comparison MB vs C were between 2.49 and 1.52, they were between 7.37 and 2.83 for OB vs C, and 6.31 and 2.66 for VB vs C. Another particularly important finding was the appearance of Epiregulin (EREG) in the top 20 DEGs list since it was a key finding of Xavier et al and was also among the "top genes differentially expressed between BD cases and controls" in their study (Table 3) [8].

Cluster analysis
The results of the cluster analysis are displayed in

Gene Ontology term enrichment analysis
Summary of the key findings of the GO term enrichment analysis is presented in Table 4. For each class comparison (i.e., MB vs C, OB vs C, and VB vs C), the top 10 GO categories with  Remarkably, several GO categories with relevance to interleukin-8 (IL-8) production (for MB vs C) and immune response to microorganisms (for OB vs C) were enriched. Additionally, biological processes of cytokine production regulation, leukocyte activation, and immune response gained substantial prominence.

Discussion
For BD, the appropriateness of the use of "Behçet's syndrome" instead of "Behçet's disease" has been previously suggested [22,23]. In support of this recommendation, it was proposed by   Lehner et al that distinct immunological abnormalities were underlying distinct classification groups of BD [6]. To date, the epidemiological basis of and the genetic linkage in BD has been pretty well studied. Nevertheless, in the era of "multi-omics", omics data, particularly genomewide transcription data is scarce in BD. In the present study, by borrowing the expression profiling data of Xavier et al and implementing a "data mining" approach, it was demonstrated that: (1) BD patients demonstrate distinct expression profiles in distinct disease subsets; (2) Different disease associated pathways seem to be functional in different disease expressions of BD; and (3) Four functionally related gene groups, namely, negative regulators of inflammation, neutrophil granule proteins, antigen processing and presentation proteins, and regulators of immune response are differentially expressed in BD patients with respect to healthy controls [8].
The immunological aberrations underlying the clinical manifestations of BD is comprehensively studied and reviewed elsewhere [24,25]. As previously stated, BD is a chronic relapsing multisystemic inflammatory disorder with a strong genetic background. HLA-B51 allele, a MHC class I gene, is shown to be a causal risk determinant for BD. Infectious agents including Table 5. Matches between the loci identified in the genome-wide association and the genome-wide linkage studies of Behçet's disease and the loci of the differentially expressed genes documented in the present study. a

GWAS/GWLS Loci
Differentially Expressed Genes with Overlapping Loci Remarks some common bacteria (e.g., streptococci) and viruses (e.g., HSV) seem to have a role in triggering the immune responses in BD patients [24,25]. While neutrophil leukocyte hyperactivity is a well-documented and central theme in BD, γδ T lymphocytes, which have functions in both innate and adaptive immune responses show distinct expansion patterns during periods of increased BD activity [24,25]. With regard to adaptive immune responses, increased levels of IL-12 and a consequent Th1 response is characteristic of BD and human heat shock proteins seem to be targets for adaptive immune responses as a result of their homologies with certain streptococcal and/or mycobacterial peptides [24,25]. Although not implicated in disease pathogenesis, various autoantibodies (e.g., antibodies against α-enolase, α-tropomyosin, kinectin) are also detected in certain subsets of BD patients [24]. Finally, endothelial cell injury is another important heading in the immunopathogenesis of BD which probably is responsible for the well-known prothrombotic state of this disorder [24]. Before going deep into the discussion of the findings, it is necessary to touch on the aberrant behaviour of the sample VB1. This 40 year old female patient with a BD diagnosis of 5 years duration was reported to demonstrate "oral and genital aphtosis, pseudofolliculitis, erythema nodosum, and large vein thrombosis" as her "major clinical symptoms" and therefore was subgrouped as VB initially [8]. Unexpectedly, during initial verification cluster analysis, VB1 consistently clustered with MB (Fig 2). As a potential explanation for this finding, we propose the possible association of a hereditary and/or acquired hypercoagulable state as the primal cause of vascular thrombosis in this mucocutaneous BD patient. It is a well-documented fact that various hypercoagulable states associated with increased risk of thrombosis contribute to the intrinsic prothrombotic state of BD [26]. Therefore, in the case of VB1 a search for thrombophilia seems relevant and may prove worthwhile.
The results of the class comparison analysis revealed strong evidences of an immunogenetic heterogeneity in different disease expressions of BD. First of all, pooling and collectively comparing the BD patients with controls (i.e., B vs C) seemed to have a pronounced attenuating effect on the number of DEGs (Table 2). Conversely, the class comparisons of BD subsets both with C and among themselves yielded substantially increased number of DEGs (Table 2). When taken together, these two findings point to a reciprocal gene expression pattern in different subsets of BD patients which was exactly the case for some of the DEGs (e.g., CD69, LTF, CEACAM8, OLFM4) as documented in Tables 3 and 6. This pattern of opposite immunological findings is a well-known concept in BD (e.g., conflicting reports of increased, normal or decreased neutrophil functions) [27].
When taken together with the relatively limited number of DEGs found in the class comparison MB vs C, the modest FC values observed may implicate that, among BD subsets, MB has the least difference in gene expression patterns compared to controls (Tables 2 and 3). Consistently, BD patients with only the mucocutaneous manifestations of the disease are widely recognized as having the mildest presentation of the disease.
Another evidence of immunogenetic heterogeneity came from the Venn diagram analysis of the class comparisons. As shown in Fig 3, the number of the common DEGs in the binary intersections of the class comparisons were markedly limited (i.e., 3, 17, and 11), while the same number in the intersection "MB vs C" \ "OB vs C" \ "VB vs C" was zero. This was a striking finding demonstrating that not even a single DEG was shared among the class comparisons of BD subsets with C; again indicating an important degree of pathogenetic heterogeneity among BD subsets.
An additional evidence was provided by the results of the cluster analysis. Using a gene set of 373 DEGs (ANOVA, p0.001), the clustering experiment effectively clustered BD patients into three clusters which exactly matched the manifestation based clusters of BD patients (Fig  4). The chosen level of significance, the number of the DEGs employed, and the success of clustering offered supporting evidence for an immunogenetic heterogeneity in distinct disease expressions of BD.
The results of the GO term enrichment analysis were also supportive. It appeared that GO categories with relevance to IL-8 production (MB vs C) and immune response to microorganisms (OB vs C) were prominently and differentially enriched (Table 4). IL-8, which is also known as "neutrophil chemotactic factor", plays a central role in neutrophil functions by inducing both chemotaxis and phagocytosis [28]. The mucocutaneous lesions (e.g., Pathergy reaction, pustular folliculitis) which are the hallmarks of BD, characteristically demonstrate significant neutrophilic infiltrates [29][30][31]. Additionally, IL-8 has previously shown to be increased in BD patients [32,33]. Thus, enrichment of the GO terms relevant to IL-8 production was consistent with the literature.
Infectious agents, especially streptococci have long been pointed to as etiologic/triggering factors in BD [2,[34][35][36][37]. With regard to an infectious etiology, an indirect mechanism involving heat shock proteins (HSP) and a cross-reactivity/molecular mimicry etiology have been postulated among many others [38]. An important potential link between ocular involvement in BD and the streptococci may be the streptococcus-related bes-1 gene derived peptides, which are shown to demonstrate a high level of homology with human retinal protein Brn-3b and HSP60 [39,40]. As such, it was noteworthy to find the immune response to microorganisms related GO categories enriched for the class comparison OB vs C.
Although many different definitions exist, "disease" can be defined as "a definite pathological process having a characteristic set of symptoms and signs" while "syndrome" as "the aggregate of symptoms and signs associated with any morbid process" [41,42]. While currently no molecular level differentation of the terms "disease" and "syndrome" is possible, as authors we strongly recommend the preference of "Behçet's syndrome" (BS) instead of "Behçet's disease", based on the findings of the present study.
As is known, the distinguishing features of BS are its recurrent mucocutaneous lesions. Additionally, whether the International Study Group (ISG) or the International Team for the Revision of the International Criteria for Behçet's Disease (ITR-ICBD) criteria are used, the diagnosis / classification of BS requires the presence of certain characteristic mucocutaneous lesions (4 out of 5 and 4 out of 6 criteria are mucocutaneous in origin in the ISG and the ITR-ICBD criteria sets respectively) [10,43]. Furthermore, while oral aphthous ulcer is a "must" in the ISG criteria, genital aphthous ulcers have more diagnostic value than other criteria in the ITR-ICBD set [10,43]. As summarized in Table 1, all BS patients in the current study were also sharing mucocutaneous manifestations irrespective of their BS subsets (i.e., MB, OB or VB). We believe that this resemblance of BS patients is of importance, not only from a diagnostic / classification perspective but also from an etiopathogenetic point of view. When the Venn diagram representation of the class comparisons MB vs C, OB vs C, and VB vs C is taken into consideration (Fig 3), the close resemblance of BS patients with respect to their mucocutaneous manifestations is inexplicable with no common DEGs in the intersection "MB vs C" \ "OB vs C" \ "VB vs C". Therefore, it is crucial to elucidate the molecular pathogenetic mechanisms responsible for the common mucocutaneous manifestations observed in different BS subsets displaying disparate sets of DEGs and distinct disease pathways.
When the FC values and the genomic loci of the DEGs were specifically taken into consideration, negative regulators of inflammation (CD69, CLEC12A, CLEC12B, TNFAIP3), neutrophil granule proteins (LTF, OLFM4, AZU1, MMP8, DEFA4, CAMP), antigen processing and presentation proteins (CTSS, ERAP1), and regulators of immune response (LGALS2, BCL10, ITCH, CEACAM8, CD36, IL8, CCL4, EREG, NFKBIZ, CCR2, CD180, KLRC4, NFAT5) were found to be differentially expressed in BS patients with respect to controls (Tables 3, 5 and 6). If the fundamental pathogenetic mechanism of BS is defined as a pro-inflammatory, innateimmune system derived response sustained by acquired immune system responses against environmental and/or self-antigens, it is motivating to note the congruences between the definition and the above-listed gene groups [34].
Albeit deserving a comprehensive and rigorous discussion which is beyond the limits and the main theme of this manuscript, the negative regulators of inflammation merit special consideration. CD69, CLEC12A, CLEC12B, and TNFAIP3 are among the well-documented inhibitors of inflammation/immune response [44][45][46][47][48]. Recently, Zhou et al reported a novel autoinflammatory disorder (haploinsufficiency of A20 [HA20]) occurring as a result of loss-offunction mutations in the TNFAIP3 gene [48]. It is remarkable to note that, phenotypically, HA20 closely resembles BS with the occurrence of recurrent oral ulcers, dermal abscesses, positive pathergy response, and retinal vasculitis [48]. In the present study, we documented that, at the transcriptomic level TNFAIP3 was downregulated in the OB and VB subsets of BS patients (Table 7 and Fig 5). CD69, yet another gene responsible for the regulation of inflammation was also shown to be significantly underexpressed in the VB subset (Table 7 and Fig 5). CLEC12A followed a similar trend with underexpression in the VB subset when compared to both the MB and OB subsets. Recently, this finding was communicated as a preliminary result to support a hypothesis about the role CLEC12A plays in BS [49]. Taken together, these findings indicate that, in patients with severe forms of BS, negative regulators of inflammation are underexpressed compared to controls and/or patients with milder presentations of the syndrome. Such a downregulation of inhibitors of inflammation may well be responsible for the pro-inflammatory milieu which is characteristic of BS [34]. Conversely stating, in patients with BS, increased expression of negative regulators of inflammation may serve a protective role against severe forms of the syndrome.
This study may be a good example for data mining. By borrowing the gene expression profiling data of Xavier et al, keeping a different perspective, and implementing a novel strategy, our group was able to document significant molecular level discrepancies among BS patient subsets [8,12]. In their original paper Xavier et al combined gene expression profiling with association studies to elucidate BS's genetic background [8]. Finally, they were able to document that EREG-AREG and NRG1 (members of the epidermal growth factor family), seemed to modulate BS susceptibility through both by direct effects and by gene-gene interactions [8]. Their study strongly emphasized the value of combining "omics" strategies (integration of "omics" data) to reveal the genetic background of complex diseases. Nevertheless, in their study Xavier et al collectively analysed BS patients and did not implement any manifestation based clinical grouping [8]. We believe that their approach has at least two important justifications; the first one being to keep a time honored approach implemented in BS research. With the exception of a limited number of studies which mainly investigate ocular BS in isolation, the current literature harbors research which analyze BS patients in a collective manner regardless of their clinical picture. The second one is due to the "omics" integration design of their study. The genome-wide association studies and the validation patient dataset used by Xavier et al belonged to collective / inclusive sets of BS patients [8]. As such, for their integration study, Xavier et al used the gene expression profiling data of a collective set of BS patients [8].
As is the case with any scientific research, the present study also is not without limitations. The well-known fact about the marked regional variations in the expression of BS necessitates the interpretation of our findings in the context that they belong to a Portuguese population [5,8]. Because of ethical considerations, continuing therapeutic regimens of BS patients had not been interrupted with potential implications in their expression profiles [8]. Also, in addition to a limited number of patients in each BS subset, BS patients with gastrointestinal, musculoskeletal, or central nervous system involvement were not represented [8]. Therefore, new expression profiling studies enrolling large numbers of treatment-naive BS patients with a wide spectrum of manifestations are clearly needed. The authors also think that in BS, eQTL (expression quantitative trait loci) analysis which simultaneously explore genome-wide expression and genetic variation data will prove worthwhile [50,51]. The need for validation of the findings in an independent BS cohort may be another issue regarding limitations. Regrettably, this validation was not possible due to lack of an independent BS cohort's peripheral blood mononuclear cells expression profiling data [9,52]. Furthermore, the marked regional variations observed in the expression of BS complicate the matter and necessitate that, such a validation data should belong to a Portuguese population.

Conclusions
BS patients display distinct expression profiles and different disease associated pathways in distinct subsets of the disorder. IL-8 production and immune response to microorganisms categories are differentially enriched among BS patient subsets. Future research, especially the studies focusing on a molecular level should take into account the immunogenetic heterogeneity of BS subsets. Based on these discrepancies, the designation as "Behçet's syndrome" should be encouraged.
Negative regulators of inflammation, neutrophil granule proteins, antigen processing and presentation proteins, and regulators of immune response appear to be instrumental in BS immunopathogenesis. Some of these genes/gene products may prove to be specific, effective, and low toxicity therapeutic targets in BS.