Comparative Transcriptomics of Portuguese and Turkish Behçet Syndrome Patients: Are We Looking at the Same Behçet?

1. Introduction

Behçet syndrome (BS), also known as Behçet disease, is a chronic, relapsing, multisystemic inflammatory disorder of unknown etiology [1]. Initially described by Dr. Hulusi Behçet as a “triple symptom complex” of oral and genital ulcers and uveitis, BS is currently classified as a complex autoinflammatory disease characterized by a variable vessel vasculitis pathology [2–4]. In addition to the hallmark mucocutaneous manifestations of BS, the syndrome may also present with diverse ocular, musculoskeletal, gastrointestinal, pulmonary, cardiovascular, and central nervous system lesions [1,4,5].

Since its first definition in 1937, some important characteristics of BS have drawn considerable attention, albeit with no clear explanations [2]. The classical high prevalence geography of BS coinciding with the ancient trade route “Silk Road,” the discordant, sometimes paradoxical immunological findings present in BS studies, the divergent clinical phenotypes of the syndrome observed in different geographic regions, and the occurrence of predefined disease clusters with shared clinical manifestations among the BS patients are some of these characteristics waiting for explanations [4,6–8]. Still another interesting observation is that, even though BS occurs all around the globe, except for its characteristic “Silk Road” geography BS tends to follow a milder clinical course [9,10].

Previously, Oguz et al. were able to document significant gene expression and molecular disease mechanism differences among BS clinical subsets/phenotypes in two different populations (i.e., Portuguese and Turkish) [11,12]. The presence of such remarkable discrepancies, even among the BS patients of the same population, brings to mind the case of BS patients of different ethnicities.

The present study, which contrasted the gene expression profiling data sets of two populations of different ethnicities (i.e., Portuguese and Turkish), first documented the differentially expressed gene (DEG) sets and then analyzed these DEG sets further with multiple bioinformatics tools and thus aimed to display molecular level similarities and/or discrepancies between Portuguese and Turkish BS patients.

2. Materials and Methods

The design of the study enabled the researchers to proceed without receiving an ethics committee approval. The flowchart summarizing the study is presented in Figure 1.

2.1. The Gene Expression Profiling Studies by Xavier et al. and Oguz et al.

Although their designs showed certain differences, the profiling studies by Xavier et al. and Oguz et al. both included BS patients diagnosed according to the internationally recognized criteria, used GeneChip^® Human Genome U133 Plus 2.0 (Affymetrix, Santa Clara, CA, USA) microarrays for hybridization, and were conducted and reported in accordance with the Minimum Information About Microarray Experiment (MIAME) guidelines [12–14].

2.2. Retrieval of the Expression Data

The raw microarray data of the study by Xavier et al. was previously retrieved from Gene Expression Omnibus (GEO) data repository using the GEO Accession GSE17114 on August 25th, 2015 [13,15]. The metadata, raw, and final/normalized data of the study by Oguz et al. were deposited to GEO with the accession GSE209567 on July 22nd, 2022 [12,15].

2.3. The Demographic and Clinical Data of the Study Population and the Definition of the Clinical Subsets of Patients with Behçet Syndrome

The demographic and clinical data relevant to the study population is extracted from the papers of Xavier et al. (DOI: 10.1007/s00109-013-1022-4) and Oguz et al. (DOI: 10.1002/iid3.836) [12,13]. The BS patients included in these two studies were grouped into four clinical phenotypes as mucocutaneous BS (MB), ocular BS (OB), vascular BS (VB), and other BS (DB), according to their individual clinical characteristics and the criteria presented in Table 1. The group inclusive of all BS patients was named B, while the control group was given the name C. For Portuguese individuals, the initial P and for Turkish individuals, the initial T was also added.

2.4. Excluded Samples/Cases

Following Oguz et al.’s two studies, the sample GSM428038 from Xavier et al.’s dataset (as an outlier in the PVB group) and the samples GSM6376933, GSM6376934, GSM6376935, GSM6376936, GSM6376937, GSM6376938, and GSM6376939 from Oguz et al.’s dataset (as complex phenotype BS patients in the TDB group) were excluded from the study (Figure 1) [11,12].

2.5. Statistical Analysis of the Demographic Characteristics

As the demographic data of the study population did not follow a normal distribution, the median as a measure of central tendency, and the minimum and the maximum values were presented. For the comparisons between two independent groups, the Independent-Samples Median Test and the Chi-Square tests were used. The p values of the analyses were presented with their absolute numerical values and a p value of ≤0.05 was considered statistically significant. Statistical analyses of the demographic data were performed using “SPSS for Windows, Version 16.0” software (SPSS Inc. Released 2007. SPSS for Windows, Version 16.0. Chicago, SPSS Inc.).

2.6. Pre-Processing of the Microarray Data

Before continuing with the downstream analyses, pre-processing of the microarray data was performed using BRB-ArrayTools Version 4.6.1 Stable Release developed by Dr. Richard Simon and BRB-ArrayTools Development Team (Figure 1) [16]. The gene expression data present as individual raw .cel files was collated by the data import function of BRB-ArrayTools. The Robust Multiarray Average (RMA) algorithm was implemented for the normalization of the microarray data [17]. Following the normalization step, the replicate spots within each array were averaged and finally, gene filters were used which helped to exclude genes that were not likely to be informative. For this purpose, a “minimum fold-change filter” in addition to a “percent missing filter” were used, which were ready-to-use functions of the BRB-ArrayTools software [16].

2.7. Bioinformatics Analyses of the Gene Expression Data

Class comparison analyses were performed using BRB-ArrayTools (Version 4.6.1 Stable Release) [16]. Two-sample t-test with a random variance model was implemented for class comparison analysis of two classes. Differentially expressed genes were selected using a p value of ≤0.05 and fold changes (FC) of ≥1.5 and ≥2.

Comparisons of the DEG sets with Venn diagrams were performed using Venny 2.1 by Juan Carlos Oliveros [18]. For Venn diagram analyses, DEGs with p≤0.05 and FC≥2 were used.

BRB-ArrayTools’ (Version 4.6.1 Stable Release) built-in clustering tools, Cluster 3.0 (developed by Michael Eisen and the Stanford group and modified by Michiel de Hoon), and Java Treeview softwares were used for clustering analyses [16,19,20]. A hierarchical clustering algorithm using Euclidean distance metric and average linkage was implemented and both the samples/cases and the genes were clustered.

Functional enrichment analyses of the DEG sets were performed with the WEB-Based Gene Set Analysis Toolkit (WebGestalt) and the chosen method of Over-Representation Analysis (ORA) specifically focused on the functional database of gene ontology biological process [21]. In order to provide sufficient number of genes with the uploaded gene lists to WebGestalt, DEG lists with p≤0.05 and FC≥1.5 were used. For the reference set, “affy hg u133 plus 2” was chosen.

3. Results

3.1. Demographics of the Study Population

The essential demographic characteristics of the study population are presented in Table 2 (Supplementary Material Table S1). Regarding the basic demographic characteristics, the Portuguese and Turkish populations were similar.

3.2. Pre-Processing Results of the Microarray Data

The box-and-whisker plots of the normalized microarray expression data are presented in Figure 2.

Following the pre-processing steps performed on the collated 65 .cel files, 16060 probe sets out of the 54675 present on the GeneChip^® Human Genome U133 Plus 2.0 microarrays passed the adjusted filters and were used during the subsequent bioinformatics analyses.

3.3. Class Comparison Analyses

Details regarding the number of the DEGs obtained during the class comparison analyses are shown in Table 3 (country/population specific) (Supplementary Material Table S2 [FC≥1.5] and Supplementary Material Table S3 [FC≥2.0]) and Table 4 (between the two populations) (Supplementary Material Table S4 [FC≥1.5] and Supplementary Material Table S5 [FC≥2.0]).

The top 10 most differentially expressed genes documented during the class comparison analyses (p≤0.05 and FC≥1.5 for all) are listed according to their FC values and presented in Table 5 (country/population specific) and Table 6 (between the two populations). Remarkably, CA1 (Carbonic Anhydrase 1) gene was the one and only shared gene, when the country/population specific class comparisons’ top 10 lists were intersected (Table 5). Likewise, genes (i.e., CD69, IFI27, CLC, CLEC12A) which were previously shown to potentially contribute to BS pathogenesis in Portuguese and Turkish populations, were not shared among the country/population specific top 10 lists (Table 5).

3.4. Venn Diagram Analyses

The illustrations of the PB vs. PC ∩ TB vs. TC, PMB vs. PC ∩ TMB vs. TC, POB vs. PC ∩ TOB vs. TC, and PVB vs. PC ∩ TVB vs. TC Venn diagram analyses are presented in Figure 3 (Supplementary Material Table S3). It was striking to see the PB vs. PC ∩ TB vs. TC and the PMB vs. PC ∩ TMB vs. TC intersection sets empty (i.e., no shared genes) and the POB vs. PC ∩ TOB vs. TC intersection set with only 2 genes (Figure 3).

3.5. Cluster Analyses

The dendrogram and heatmap representations of the results of the cluster analyses of the Portuguese and Turkish BS patients and patient subsets are presented in Figure 4 (Supplementary Material Table S6). As shown in Figure 4, the gene sets used and the algorithm implemented successfully clustered both Portuguese and Turkish BS patients and Portuguese and Turkish BS clinical subsets. Important to note, the gene sets used for clustering purposes were constructed from the top 10 most differentially expressed genes observed during the population specific class comparison analyses.

3.6. Functional Enrichment Analyses

The results of the population specific functional enrichment analyses performed individually with DEG sets with increased and decresed expression are schematically presented in Figure 5 (DEGs with increased expression) and Figure 6 (DEGs with decreased expression) and given in detail in Supplementary Material Table S7. Strikingly and unexpectedly, not even a single gene ontology biological process term was common in the Portuguese and the Turkish populations’ functional enrichment analyses top 15 terms (Figure 5 and Figure 6).

4. Discussion

This study sought to investigate the molecular immunopathogenetic similarities and differences underlying BS across two different ethnic populations, Portuguese and Turkish. As Oguz et al.’s two previous studies pointed to the potential pathogenetic role of certain innate immunity genes (i.e., CD69, CLEC12A, CLC, and IFI27) in BS, specific focus was on BS’s fundamental innate immunity pathology [11,12]. Our main findings highlighted significant gene expression and molecular disease mechanism differences between Portuguese and Turkish BS patients, particularly regarding the anti-inflammatory and the pro-inflammatory genes. As the authors we comment that, the anti-inflammatory and the pro-inflammatory genes are the key to understand the immunopathogenesis of BS [12] and the documented expression differences in this study among the two different ethnicities may account for the observed discordant, sometimes paradoxical immunological findings present in BS studies and also the divergent clinical phenotypes of the syndrome observed in different geographic regions. Additionally, the distinct DEGs and the different molecular mechanisms observed across ethnicities, reinforces the concept of BS as a syndrome rather than a single disease entity.

Analysis of demographic characteristics revealed that BS and control groups, as well as the three BS clinical phenotypes (i.e., mucocutaneous, ocular, and vascular), were comparable in terms of age and gender distribution across the two different ethnicities (Table 2). While Turkish BS patients showed a slight male predominance, Portuguese BS patients had an equal Male/Female ratio, which aligns well with the existing literature [4,8,10]. Furthermore, a gender-based difference in severity was observed among Turkish BS patients, males were more prominent in vascular BS, while the mucocutaneous subset showed a Male/Female ratio of 1/7 (Table 2). This pattern was absent in the Portuguese population, consistent with the literature indicating a male predominance in more severe BS phenotypes [4,8].

Class comparison analyses revealed notable differences between populations. When comparing each population’s BS patients to its respective control, the DEG counts ranged from single to triple digits (min.: 4, max.: 549, p≤0.05, FC≥2.0) (Table 3). In contrast, comparisons between Portuguese and Turkish BS groups revealed DEG counts in the four-digit range (min.: 7466, max.: 8642, p≤0.05, FC≥2.0) (Table 4). Subtracting the Portuguese Control vs. Turkish Control class comparison’s DEGs (6549, p≤0.05, FC≥2.0) did not substantially lower these DEG counts and this finding suggests that, the gene expression differences observed in BS groups are not solely attributable to ethnic background. These substantial molecular discrepancies between the BS groups reinforce the idea that distinct pathways underlie BS across different populations [1]. Comparison of the top 10 DEGs demonstrated minimal overlap between Portuguese and Turkish BS patients, with only CA1 shared among the two populations (Table 5). Notably, previously identified innate immunity genes with potential pathogenetic contribution to BS (i.e., CD69, IFI27, CLEC12A, and CLC) were not shared in the top 10 DEG lists of the two populations (Table 5). This finding indicates that while gene functions may be similar, the exact molecular drivers differ across populations in BS. The same finding also aligns with our hypothesis that the gene expression and molecular mechanisms involved in BS likely vary by ethnicity, further supporting the notion of BS as a heterogeneous syndrome.

In Portuguese BS patients, CD69 showed increased expression across patient subgroups relative to Turkish counterparts (Table 6). The CD69 gene, previously associated with milder disease expression in Portuguese BS patients by Oguz et al., suggests a protective, anti-inflammatory role that might contribute to milder BS manifestations in Portuguese patients [11,22–24]. Conversely, the involvement of genes such as IFI27, CLEC12A, and CLC in Turkish BS patients suggests a more severe clinical presentation in this population (Table 5). Together, these findings support a model where distinct genetic profiles influence BS severity across ethnicities, potentially providing a molecular basis for population-specific therapeutic strategies in future.

Venn diagram analyses revealed no shared DEGs in the intersections of PB vs. PC ∩ TB vs. TC and PMB vs. PC ∩ TMB vs. TC (Figure 3). It was striking, especially not to see even one shared gene in the PB vs. PC ∩ TB vs. TC set. When the relatively more severe forms of BS (i.e., ocular and vascular BS) were analyzed, the number of shared DEGs increased (2 for POB vs. PC ∩ TOB vs. TC [out of 166 genes] and 32 for PVB vs. PC ∩ TVB vs. TC [out of 733 genes]) (Figure 3). This finding may be reflecting gene expression modifications secondary to prolonged and severe tissue damage and therapeutic effects. In cases of severe BS, especially in those patients with organ and/or life threatening disease, common immunosuppressive or immunomodulatory treatments are often used to limit the disease activity with potential consequences in gene expression profiles [25].

Clustering analyses successfully distinguished Portuguese and Turkish BS patients and their respective clinical phenotypes using only around 20 genes per analysis (Figure 4). Notably, IFI27 emerged as the only shared gene across clustering DEG lists, underlining its potential importance as a key pro-inflammatory gene in Turkish BS patients [12]. This finding may (1) support the view that IFI27 may play a pivotal role in distinguishing BS patients across ethnicities, (2) emphasize its relevance in BS pathogenesis in Turkish population, and (3) further underline innate immunity and inflammatory mechanisms as central in BS pathogenesis. Also, as the authors we hypothesize that, given the role of IFI27 in driving pro-inflammatory pathways, it may represent a key therapeutic target for BS management [12,26,27].

The functional enrichment results further confirmed significant molecular differences between Portuguese and Turkish BS patients. No shared enriched gene ontology terms were identified among the top 15 in either population whether using DEGs with increased or decreased expression (Figure 5 and Figure 6). Functional enrichments differed significantly between Portuguese and Turkish BS patients, with Portuguese patients showing pathways related to neutrophil degranulation and immune response, while Turkish patients displayed enrichment in Type I interferon signaling and T-cell receptor pathways (Figure 5 and Figure 6). This lack of common functional enrichments between the two populations again supports the hypothesis that BS’s molecular mechanisms vary by ethnicity.

We must acknowledge several limitations regarding the study. First, this in-silico study was relying on microarray data, which, while valuable, lacks the advantages of next-generation sequencing methods. As is well known, RNA sequencing could capture rare transcripts and identify genetic variants which are not detectable in microarrays, potentially expanding our understanding of BS pathogenesis [28]. Furthermore, this study was limited to Portuguese and Turkish populations, which clearly is not able to represent BS heterogeneity globally. Additionally, treatment effects on gene expression cannot be ruled out as BS patients enrolled in the study were known to continue their pharmacological therapies. Finally, the present study did not include BS patients with complex clinical phenotypes involving multiple organs/organ systems. Future research with larger treatment-naive BS patients, across more diverse geographic regions could provide further insights into BS’s immunopathogenesis and immunopathogenetic diversity.

5. Conclusions

In conclusion, this study demonstrated significant gene expression and molecular disease mechanism differences between Portuguese and Turkish BS patients. This finding strongly supported the complex autoinflammatory disorder classification of BS and added that, in addition to its complex genetic nature, BS is also a heterogeneous syndrome. The observed interethnic differences of gene expression and molecular mechanisms may account for the inconsistent, sometimes contrasting immunological findings observed in BS research and the distinct clinical phenotypes of BS cumulating in certain geographic regions of the globe. Decreased expression of anti-inflammatory genes was a common theme documented in both Portuguese and Turkish BS patients, albeit the identities of these genes differed across populations. This finding occurred to such an extent that, important anti-inflammatory genes with potential pathogenetic roles in one population were not detected in the other. Additionally, if present, pro-inflammatory gene expression differences seemed to enhance disease severity. Future research should explore the anti-inflammatory pathways, especially involving innate immunity proteins like CD69, CLEC12A, and CLC, as potential therapeutic targets to reduce inflammation in BS. These insights signal the forthcoming era of personalized therapeutics based on molecular profiles for BS patients.