Novel Insights Into Gene Signatures and Their Correlation With Immune Infiltration of Peripheral Blood Mononuclear Cells in Behcet’s Disease

Background Behcet’s disease (BD) is a chronic inflammatory disease that involves systemic vasculitis and mainly manifests as oral and genital ulcers, uveitis, and skin damage as the first clinical symptoms, leading to gastrointestinal, aortic, or even neural deterioration. There is an urgent need for effective gene signatures for BD’s early diagnosis and elucidation of its underlying etiology. Methods We identified 82 differentially expressed genes (DEGs) in BD cases compared with healthy controls (HC) after combining two Gene Expression Omnibus datasets. We performed pathway analyses on these DEGs and constructed a gene co-expression network and its correlation with clinical traits. Hub genes were identified using a protein–protein interaction network. We manually selected CCL4 as a central hub gene, and gene-set enrichment and immune cell subset analyses were applied on patients in high- and low-CCL4 expression groups. Meanwhile, we validated the diagnostic value of hub genes in differentiating BD patients from HC in peripheral blood mononuclear cells using real-time PCR. Results Twelve hub genes were identified, and we validated the upregulation of CCL4 and the downregulation of NPY2R mRNA expression. Higher expression of CCL4 was accompanied by larger fractions of CD8 + T cells, natural killer cells, M1 macrophages, and activated mast cells. Receiver operator characteristic curves showed good discrimination between cases and controls based on the expression of these genes. Conclusion CCL4 and NPY2R could be diagnostic biomarkers for BD that reveal inflammatory status and predict vascular involvement in BD, respectively.


INTRODUCTION
Behcet's disease (BD), considered a variable vessel vasculitis, is characterized by a wide spectrum of clinical features including mucocutaneous lesions, non-granulomatous uveitis, lower extremity vein/cerebral venous sinus thrombosis, aortic aneurysms, and neurological involvement (1). The diagnostic criteria of the international study group for BD requires oral aphthous as the entry symptom, with the concurrence of genital ulcers, ocular or vascular involvement, as well as pathergy test positivity, which renders an extreme specificity (97.5%) (2). Nonetheless, its dynamic manifestation with disease relapse and remission periods, distinct population and sexual clinical expression, and disease type subsets pose grand challenges in identifying novel biomarkers for the diagnosis of BD.
Genetic risk factors impact BD pathogenesis, as the occurrence of BD is 5.78 times higher in HLA B51/B5 carriers than in noncarriers (3). BD susceptibility loci, for instance IL-10, STAT, and IL-23R polymorphisms (4,5), indicate the predominant pathogenetic role of inflammation in the disease etiology. MicroRNA (miRNA), influenced by inflammatory signaling pathways, regulates mRNA expression through pre-and posttranslational processes, and miRNA variants have shed light on disease status (6). A recent study in Turkey identified the FAS rs1800682 polymorphism and mir146a as protective factors against BD (7). Ibrahim and Jadideslam confirmed the association between mir146a expression and disease activity, with involvement of the eyes and vasculature (8,9). Downregulation of miR-155 and that of miR146a likely exert their effects by disturbing Th1/Th17 cell and cytokine homeostasis (10). Profiling analysis of miRNA expression has also identified epigenetic pathways that are enriched in vascular biology, the immune response, and apoptosis (11). Long noncoding RNAs (lncRNAs) are also involved in BD etiology, with the homozygous T allele of rs9517723 in lncRNA LOC107984558 correlated with ocular and central nervous system (CNS) lesions, the induction of UBAC2 expression, and overactivation of the ubiquitination-related pathway (12). However, we could not only concentrate on individual gene and ignore the potential interactions between them which might be the backbone of intracellular communications and etiology pathways.
Hub genes are highly connected nodes between co-expression network, exploring of which is vital for seeking diagnostic biomarkers or therapeutic drug targets and clarify the biological regulation process in disease conditions. In this study, we utilized weighted gene co-expression (WGCNA) network and proteinprotein interaction (PPI) network to identify hub genes. Since WGCNA analysis could provide a systematic insight into coexpression gene sets and their relationships with clinical phenotypes, we invoked this algorithm to establish a scale-free network using adjacency matrix from which topological overlap matrix and corresponding dissimilarity were converted to differentiate gene co-expression modules (13). For subsequent analysis, we rigorously chose the turquoise module as a candidate module which gains highest correlation coefficient and significant corresponding p value with clinical features among patients suffering Behcet's disease and healthy controls. Recently, a few of studies have comprehensively utilize WGCNA to filtrate synergistically altered gene sets in autoimmune diseases (14)(15)(16)(17)(18). Up till now, amassed structure and function of single protein fragmented our comprehension of their intracellular interactions and pathogenesis pathways. By conducting PPI program, we could easily attain a tighter protein-based co-expressed DEG channel and a more integrated interaction network, the connectivity of which is ranked by combined score calculated from online sources including text-mining, experiments, databases and genomic context information, etc. (19). To predict hub genes bridged through networks, we employed MCODE plugin to visualize intensively related modules and CytoHubba package to screen out hub nodes avail of functional enrichment.
Despite the body of research elucidating the genetic variation contributing to BD, a deep exploration of its molecular mechanisms, interaction networks, and key pathways is lacking. In the present study, we aimed to uncover hub genes acting as biomarkers for the diagnosis, classification, and biological functions and networks of BD using bioinformatics methods and a validation cohort.

Differentially Expressed Gene and Enrichment Analysis
Differentially expressed gene (DEG) analysis was conducted using the "limma" package in R with the thresholds of an adjusted P < 0.05 and |log Fold Change (FC)| > 0.5 (21). To identify the biological functions and corresponding pathways involving these DEGs, we performed Gene Ontology (GO) enrichment and visualized the biological process, cellular component (CC), and molecular function (MF) categories and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis of signaling pathways using the R packages "clusterProfiler" and "enrichplot" (22). The significance cutoff criteria were P < 0.05 and FDR < 0.05. into co-expression modules and ultimately to identify the specific gene module which is most relevant to clinical phenotype of Behcet's disease. After shearing samples below the abline (h=20000), we constructed the co-expression network with soft thresholding power to obtain a higher level of scale free R 2 and mean connectivity (13). In dynamic tree cut and module identification section, we altered 10 as the minimum number of gene modules. Utilizing the clinical traits data containing BD patients and healthy controls from GSE70403 and GSE17114, the gene significance (GS) and module membership (MM) were calculated. The relationship between the gene modules and clinical traits was represented in the form of a heatmap, from which we selected a module that was most clinically relevant.
With a threshold score of 0.150, we constructed a proteinprotein interaction (PPI) network consisting of the DEGs using the STRING database (19). We have realized from previous reference (23) that the combined score from STRING database is meant to express an approximate confidence of the association between proteins being true based on every channel of evidence (depended on genomic context information, co-expression, text-mining, experiments and databases). Thus, we have ranked the first 10 interactions by the combined score in the PPI network. To attain a higher degree of connectivity, we applied CytoHubba and then extracted the key Hubba nodes ranked by maximum clique centrality (MCC) in Cytoscape software (24). Using the MCODE plugin (25), we sought dense regions divided by the critical value of a degree cut-off = 2, node score cut-off = 0.2, k-core = 2, and max. Depth = 100.
Furthermore, we united the first 10 genes ranked by the combined score in the PPI network and node genes filtered with an MCC > 5 in CytoHubba analysis, together with genes exported using the MCODE program, and 12 hub genes were ultimately identified for principal component analysis

Correlation of Hub Genes and Vasculitis
Genecards (https://www.genecards.org/) was screened using the keyword "vasculitis." We ranked the relevance score given by website and selected top 15 genes, then revisited the normalized gene expression data merged from Method 2.1 to see which genes from Genecards have the intersections with GEO datasets. Thus, we deleted 3 genes which were not documented in normalized gene expression data merged from GEO datasets (ADA2, HLA-DRB, C4A). Finally, the top 12 genes (PRTN3, PTPN22, CTLA4, DNASE1L3, MPO, MEFV, HLA-B, HLA-DPA1, HLA-DPB1, IL-10, TNF and CRP) were selected. According to the relevance score, the top 12 genes were selected. Moreover, we visualized the correlation heatmap using the R package "corrplot" (28).

Gene-Set Enrichment Analysis
To explore the biological functions of hub genes, we downloaded Gene-Set Enrichment Analysis (GSEA) software (http://software. broadinstitute.org/gsea/index.jsp). The patients were divided into high/low-expression groups according to the median expression of the central hub genes. Afterward, we selected "c2.cp.kegg.v6.1.symbols.gmt" as the reference gene-set in the Molecular Signatures Database (29) and diagrammed the KEGG enrichment plot of distinct expression groups by normalized P < 0.05 and FDR < 0.2.

Correlation Analysis of Gene Signatures and Immune Cells
CYBERSORT is a deconvolution algorithm used to determine the relative proportions of 22 human immune cells from bulk expression data (30). Using the R packages "e1071," "preprocessCore," and "limma," we performed immune cell profiling among the BD samples. Samples with a P < 0.05 were included. To estimate the diverse immune cell patterns in different central hub gene expression groups, we used "vioplot" in R studio.

Diagnostic Prediction and Validation of Hub Genes
To validate the expression pattern and diagnostic value of the selected hub genes, we collected blood samples from 16 BD patients and 16 age-and sex-matched HCs under a protocol approved by the Medical Ethics Committee of Peking Union Medical College Hospital with informed consent acquired from the enrolled subjects (demographic and clinical data of the validation population and controls in online Supplementary Materials, Table S3). We then isolated PBMCs by Ficoll density gradient centrifugation. Further, total RNA was extracted from PBMCs using the Trizol Reagent, and cDNA was reversetranscribed following the manufacturer's instructions (Mei5 Biotechnology, Co., Ltd). For amplification, SYBER green (Mei5 Biotechnology, Co., Ltd) was utilized for real-time PCR on a Roche LightCycler 480 system. Finally, the relative expression levels were determined using the 2 −ΔΔCt formula after normalization with endogenous GAPDH mRNA expression. The primer sequences were as follows: CCL4 forward CTTTTCTTACACCGCGAGGA and reverse GCTT GCTTCTTTTGGTTTGG, AGTR2 forward TTCCCTTCCA TGTTCTGACC and reverse AAACACACTGCGGAGCTTCT, NPY2R forward CAAAACTTCTCCTCCAGTCCCC and reverse AAGGTGGGAGGATCAGAGATGG, ZIC1 forward GCATC CCAGTTCGCTGCGCAAA and reverse GGAGACACGAT GGTGGGAGGCG, and GADPH forward CCTCAAGATC ATCAGCAAT and reverse CCATCCACAGTCTTCTGGGT. Receiver operator characteristic (ROC) curves were rendered using the R package "pROC," and DEGs between BD and HC were visualized using "ggpubr."

Screening for DEGs Between BD and HC
The workflow illustrated in Figure 1 demonstrates the procedure of exploring hub genes and diagnostic markers in BD. Two eligible datasets GSE70403 and GSE17114 with 56 BD and 14 HC combined were included. In total, 82 genes were differentially expressed in BD compared with HC (details provided in Table S1). Meanwhile, 11 DEGs were upregulated, and 71 were downregulated (Figures 2A, B).

GO and KEGG Annotation of DEGs
We investigated the pathogenesis of BD by applying functional and pathway analysis. Several DEGs were prominently enriched in "clathrin-coated endocytic vesicle membrane" and "clathrincoated endocytic vesicle" terms in the context of CC. For MF, DEGs were preferentially included in "IgG binding," "immunoglobulin binding," and "carbohydrate binding" demonstrating the immunopathogenic background of BD ( Figure 3A). Various disease-associated KEGG pathways were hyperactivated, including "Fc gamma R-mediated phagocytosis," "Staphylococcus aureus infection," and "Leishmaniasis"

Identification of Hub Genes
Five clusters of modules were constructed with WGCNA coexpression network analysis ( Figures 4A, B and Table S2). A sample dendrogram and trait heatmap calculated using Pearson's correlation coefficients indicated that there were no outlier samples ( Figure 4C). When incorporating the expression profile with clinical traits, we found that the turquoise module was significantly associated with clinical phenotypes and contained 72 genes (r = 0.60, P < 0.05) ( Figures 4D, E). A Venn plot intersecting the DEGs and the turquoise module narrowed the hub gene candidates down to 71 ( Figure 4F). Next, we input these 71 genes into PPI network analysis to investigate the interrelationship between the proteins encoded by these genes and discovered 70 interactions ( Figure 5A). We have extracted the first 10 interactions by the combined score in the PPI network. CCL4, NPY2R, AGTR2, TAS2R1, ASB14, ASB17, C1orf110, SOX14, MAGEA1, NPAS4, EYA1 and HOXA11 was elected (Set 1). We filtered the top 10 genes ordered by a combined score and an MCC > 5 filtered by CytoHubba, including ZIC1, CTXN3, NPY2R, AGTR2, LRRC3B, EYA1, HOXA11, SOX6, CCL4, CAMKV, ANGPTL3, TAS2R1, SLC6A3 and SOX14 (Set 2). ( Figure 5B and Table 1), along with subnetwork gene nodes exported using the MCODE program where NPY2R, LRRC3B, AGTR2, CTXN3, CAMKV, CCL4 and TAS2R1 was clustered as subset 1 (shown in Figure 5C-ii) and HOXA11, EYA1, SOX6 as subset 2 (shown in Figure 5C-i), they are all named as Set 3. ( Figure 5C and Table 2). Eventually, 12 hub genes (AGTR2, CAMKV, CTXN3, EYA1, HOXA11, LRRC3B, NPY2R, SOX14, SOX6, TAS2R1, ZIC1, and CCL4) were identified (shown in Figure S1: a union of the intersection of set 2 and 3 plus the intersection of set 1 and 2). ZIC1 is the top 1 significant node gene ranked by MCC, therefore we could not neglect it for further analysis and also add it into hub genes. These hub genes were identified for PCA analysis in GSE17114 to divide BD patients based on clinical manifestations (26). These hub genes could well discriminate isolated mucocutaneous manifestations (MB), ocular involvement (OB), and large vein thrombosis (VB) patients diagnosed with BD ( Figure 5D).

GSEA Analysis of Central Hub Genes
Through an intensive literature search in PubMed, we selected CCL4 as the central hub gene because of its known pathogenic role in autoimmune diseases including rheumatoid arthritis (RA) (31), multiple sclerosis (MS) (32) and Behcet's disease (33). Accordingly, our patients were categorized into high-and low-expression groups of CCL4 for GSEA function and pathway analysis. Twelve KEGG pathways, "natural killer cell-mediated cytotoxicity," "glycosaminoglycan biosynthesis chondroitin surface," "antigen processing and presentation," "regulation of autophagy," "graft versus host disease," "type I diabetes," "T cell receptor signaling pathway," "cytosolic DNA sensing pathway," "TGF-b signaling pathway," "chemokine signaling pathway," "cytokine-to-cytokine receptor interaction," and "toll-like receptor signaling pathway," were enriched in the CCL4 high-expression BD patients (Figure 7).

Landscape of Immune Cell Proportions
We calculated the relative percentage of immune cells among BD patients and HC ( Figure 8A). In comparison with the low-

Diagnostic Performance of Hub Genes and Validation
We then performed ROC analysis to detect the performance of hub genes in differentiating patients with BD from HC. The area under the curve (AUC) ranged from 0.81 to 0.88 ( Figure 9A).
Additionally, CCL4 was highly expressed in our validation cohort (P = 0.0234), whereas NPY2R was downregulated (P = 0.024) ( Figures 9B, C and Table 3); however, no significant difference was found in AGTR2 and ZIC1 expression among BD and HC individuals (data not presented).

DISCUSSION
In this current study, we knit in-silico results and the in vitro data closely then disclosed the biological functions as well as pathogenic pathways involved in Behcȩt's disease (BD). To begin with, we downloaded and merged 2 GEO datasets  (26,27). Wondering which hub gene is previously considered as the potential pathogenic manipulator for autoimmune disease in references, we retrieved PubMed and only to find CCL4, as a central hub gene, stands out in RA, MS and BD etiology. Therefore, we classified the BD patients from microarray datasets following the expression level of CCL4 for GSEA function, twelve KEGG pathways intensively enriched in intracellular and intercellular biological communications among immune cells.
Hence, CYBERSORT was utilized to visualize the immune cell profiling based on the expression level of CCL4 in BD patients. The upregulated CCL4 could recruit a higher proportion of CD8 + T cells (P = 0.019), NK cells (both resting and activated, P = 0.041 and 0.010, respectively), M1 macrophages (P = 0.044), and activated mast cells (P = 0.041) into inflammatory sites (skin, eyes and even vascular), which likely exacerbates proinflammatory situation, immune cell function and intensify disease activity in BD patients. Last but not least, the performance of hub genes exerts a great ability in discerning BD and HC (AUC ranging from 0.81-0.88) where CCL4 and NPY2R are proved as a moderate diagnostic biomarker in disease recognition from healthy subjects in our validation cohort (clinical characteristics are provided in the online Supplementary Material, Table S1).
Taken together, we have identified the immune-activated, inflammation-induced, and carbohydrate-mediated pathogenesis of BD and identified DEGs and pathways involved in its etiology.
Our results indicate that the identified hub genes could be of diagnostic value and reveal new insights into the etiological analysis of BD.
Subsequent validation disclosed that CCL4 expression is significantly upregulated in BD patients compared to HC. Additionally, we found that there was a higher proportion of

Node1
Node2 Combined_score   (37). From snap shots of CCL4/MIP-1b participation in clotting formation, there is supposed to exist a hidden mechanism in causing thrombosis of BD patients. Previous evidence supports that upregulated CCL4 could be considered to be a circulating inflammatory cytokine marker in primary and systemic lupus erythematosus (SLE)-related autoimmune hemolytic anemia, indicating erythroid compensation and disease severity due to its positive relationship with reticulocyte count (38). Being a major macrophage attractant, universal involvement of CCL4 has been found in diabetes mellitus, since its circulating concentration increased regardless of disease stage, which may contribute to destruction of islet cells and the progression of diabetes. By providing an inflammatory microenvironment, CCL4 induces cell adhesion to endothelial cells through oxidative stress and tends to augment atherosclerosis and vascular injury (39). Additionally, elevated CCL4/MIP-1b in the cerebrospinal fluid is a characteristic of IgG4 anti-neurofascin 155-positive chronic inflammatory demyelinating polyneuropathy, which implies exacerbation of protein leakage, damage of encephalic and CCL4/MIP-1b decreases with combined immune-treatment response (40). Della-Torre et al. demonstrated that B cells and fibroblasts might orchestrate inflammatory infiltration through CCL4, which in turn amplified the recruitment of monocytes, CD4 + cytotoxic T cells through CCR5, and IgG4-related disease lesions (41). In discriminating between relapsing-remitting and progressive clinical phenotypes of multiple sclerosis (MS), the plasma CCL4/MIP-1b level showed a good diagnostic value (AUC = 0.702, P < 0.0002) and was a protective factor preventing disease progression (odds ratio = 0.9873, P = 0.0171) (42). However, the pathogenic impact of CCL4 in BD has not been reported before. As the similar pathogenesis existing among autoimmune diseases, we could deduce that the expression of CCL4 might be a circulating inflammatory chemokine which aggravates inflammation and disease progression via inducing immune cell infiltration and damaging surrounding tissue in BD patients.
Neuropeptide Y receptor subtype 2 (NPY2R) is a G-protein binding receptor for neurotransmitter NPY, which triggers vascular smooth cell proliferation and vasoconstriction, and results in ischemic angiogenesis (43). In our study, NPY2R mRNA expression in PBMCs of BD patients, especially in those with severe symptoms including thrombosis and gastrointestinal involvement, was distinctly lower than that in HC. A genome-wide association study demonstrated that the NPY2R rs1902491 polymorphism was associated with a severe diabetic retinopathy subgroup without end-stage renal disease (44); however, no significant association was replicated in a small cohort with rapidly proliferative diabetic retinopathy in Lithuania (45). This may imply NPY2R is potentially correlated with microvascular leakage and obstruction which could stands for vascular destruction within BD patients. The gene encoding the type 2 angiotensin II (Ang II) receptor (AGT2R) plays a pivotal role in atherosclerosis, vascular inflammation, and remodeling and has anti-inflammatory and antiproliferative influences, and it exhibited a female advantage by counteracting the AT1 receptor in mice (46). In SLE, AGT2R polymorphisms containing the A allele (CA + AA) in exon 3 or G allele (AG + GG) in intron 1 are of apparent significance compared with HC (47). AGT2R expression is also strongly enhanced in RA synovial tissue, an exogenous agonist that could ameliorate inflammation (48). On the contrary, AGT2R blockers are postulated to dampen the clinical situation and laboratory metrics in RA patients (49). In idiopathic pulmonary fibrosis, AGT2R is highly expressed in interstitial fibroblasts and serves as a determinant of fibroblast proliferation and migration (50). In 883 Austrian MS patients and 972 controls. The ZIC1 rs1841770 polymorphism was associated with an earlier disease onset of MS (51), but a Spanish replication study did not support it as an MS predisposition factor (52). Nevertheless, we did not find any statistical difference in AGT2R and ZIC1 mRNA expression between BD and HC; this may be attributed to the small population of our validation cohort. NK cells are innate immune lymphocytes that have a cytolytic immune-regulatory effect in autoinflammatory disease, particularly in BD. We confirmed that in patients with upregulated CCL4 expression, there is a higher fraction of NK cells infiltrating the inflammatory zone. However, conflicting studies have reported decreased, normal, and increased NK cells in BD compared with those in HC (53). It has been hypothesized that impaired NK cytotoxicity and degranulation facilitate and active disease progress (54,55). The NK1/NK2 ratio is significantly increased, resulting in CD16 + IFN-g + NK1-induced IFN-g secretion and suppressed NK regulation in mucocutaneous BD. Moreover, increasing IFN-g produced by NK cells reflects active disease and relapse of disease, which demonstrated the dominant interaction between NK cells and interferons (56). Accordingly, NKMHC-I modules encoding HLA-B*51 could also contribute to NK cell inhibition through binding killer immunoglobin-like receptors, further supporting hereditary factors as key drivers of BD etiology. Hyperactivation of monocytes and neutrophils is a feature of BD (57). Our results identified a high level of monocytes in CCL4 low expression group, instead of neutrophils in BD. Monocytes in BD patients appear to have intensified oxidative burst activity compared with those in RA patients and HC (58). The enriched pathways JAK/STAT, IL-6, and interferon signaling are upregulated in CD14 + monocytes isolated from BD patients (59), whereas toll-like receptors 2, 4, and 5 are highly expressed in monocytes, causing a cascade of activation of nuclear factor-kB and production of TNF-a, unveiling a proinflammatory mechanism of BD (60). Bacterial infection might be a trigger of BD pathogenesis. Lipopolysaccharide (LPS) stimulation increases TNF-a secretion in vivo in monocytes of active BD rather than in quiescent ones (61). A recent study demonstrated a perturbed FcgRs balance in monocytes from BD patients that represented downregulated inhibitory FcgRIIb and upregulated activating FcgRIII, the level of which fluctuated along with BD treatment and disease severity. Moreover, LPS and IgG-stimulated monocytes produced more IL-6 than those from HC (57). IgG binding activates FcgRs and FcgRIIb, mainly leading to endocytosis and internalization of immune complexes; therefore, these studies corroborate our finding of "IgG binding," "clathrin-coated endocytic vesicle membrane," and "clathrin-coated endocytic vesicle" in our GO analysis and the enriched pathways "Fc gamma R-mediated phagocytosis" and "Staphylococcus aureus infection" in the KEGG analysis.
M1 macrophages stimulate an inflammatory reaction by triggering cytokines such as TNF-a and IL-6, resulting in Th1 reactivity, whereas BD patients with skin lesions conferred M1 macrophage dominance in systemic sclerosis (SSc). It confers that M1 macrophages is more like a predominant pathogenesis factor in BD than in SSc. A higher expression of CCR1 ligands for CCL3/ MIP-1a than that in HC again demonstrated that a skewed M1 and M2 macrophage balance led to MIP-1a-induced inflammation (62). Macrophage incubation with active BD serum facilitated CD64 (FcgRI) positivity and IL-8 mRNA expression compared to inactive BD, suggesting an overwhelming proinflammatory response and aberrant endothelial alteration (63). Mast cells are prominently elevated in the skin lesions of active Behcet's disease compared with those uninvolved skin (64). It seems that mast cells are attracted by chemotaxis influenced by complement cascades in BD patients with high-expressed CCL4, which is analogical to basophil hyperactivation process.
The profound role of IL-17A and CD8 + T cell-mediated immune response is pronounced in BD with papulopustular skin lesions. Investigation into Tc17 cell activation confirmed that CD8 + T cells are infiltrated in peripheral vascular walls and the subcutis in BD, and, concurrently, IL-17A + CD8 + T cells were elevated in BD compared to psoriasis vulgaris lesions and were a major source of IL-17A rather than CD4 + T cells (65). IFN-gproducing CD8 + T cells are enhanced in the peripheral blood of active BD patients and are even higher in patients with pulmonary manifestations and neurological involvement. Furthermore, the Tc1/Tc2 ratio is dramatically increased in the active and remission phases of BD (66). Increased IFN-g and CD8 + Leu7 + T cells together with an increased titer of IgG antibodies to HSV-1 demonstrated an antigen-driven T celldominant etiology characteristic of BD (67), which is in accordance with our GSEA pathway analysis results.
Our GO analysis highlighted the "carbohydrate binding" pathway, and accumulating evidence suggests that galectins regulate immune cell homeostasis, in which galectin-1 inhibits T cell proliferation, enhances apoptotic activity, and thus suppresses inflammation, whereas galectin-3 amplifies the inflammatory response and survival of T lymphocytes, which subsequentially activates NADPH oxidase, the superoxidative reaction of neutrophils, and recruits monocytes (68). This coincides with our finding of "alcohol dehydrogenase [NAD(P)+] activity" among the MF GO terms. The combination of lectins and epithelial cells and lectins binding to specific tissues triggers antigen humoral responses leading to autoimmune rheumatic responses, for instance, in RA or type I diabetes (69). S-nitrosated mannose binding lectin (SNO-MBL) is upregulated in RA patients and was dysfunctional in depositing C4, phagocytosing bacteria, and binding to apoptotic cells, which consequently caused enhancement of anti-SNO-MBL antibodies (70). There are several limitations in the present study. First, although we have verified significant difference between BD and HC group in mRNA expressions of CCL4 and NPY2R, we could not associate our results with disease activity index (Behcȩt's Disease Current Activity Form) attributed to lack of integrated clinical information of each patient. Second, we aimed to identify a diagnostic marker in discerning BD patients and healthy subjects, thus we did not add disease controls such as takayasu arteritis, inflammatory bowel disease, ANCA-associated systemic vasculitis and even patients with either coronary or venous thrombosis which shares semblable clinical feature with BD. Third, via MCODE program we have discovered the protein interaction network between CCL4 and NPY2R, future analysis should more concentrate on underlying up/down regulation between them which might be the mechanism leading to vasculitis moreover vascular deteriorations in BD. Further studies using more comprehensive data and larger validation cohort should be performed to elucidate the association of CCL4 and NPY2R with BD activity as well as clinical phenotypes.
In conclusion, as a proinflammatory factor synthesized in the inflammation zone, higher mRNA expression of CCL4 might be an upstream switch to attract more immune cells into inflammatory sites, these cells in turn positively amplify cytokine secretion and exacerbates BD disease stage. In vasculitis related to BD, NPY2R expression is potentially impaired and may be a marker of vascular complications. Future studies are required to elucidate the potential mechanisms underlying the dysregulated expression of target genes and their correlation with the immune microenvironment and clinical phenotypes in BD.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material. Further inquiries can be directed to the corresponding author.

ETHICS STATEMENT
The studies including human participants were approved by Medical Ethics Committee of Peking Union Medical College Hospital with informed consensus acquired from enrolled subjects. The patients/participants provided their written informed consent to participate in this study.

AUTHOR CONTRIBUTIONS
YL conceived and designed the research. HZ extracted data, performed software analysis and visualized graphs and tables. LC, SY, HL, and HZ categorized and supervised graphs and tables. WZ provided the clinical data of each patient. HZ and HL went on validation experiments. HZ wrote the paper. All authors are accountable for all aspects of the study, and attest to the accuracy and integrity of the results. All authors contributed to the article and approved the submitted version.