Functional Genetic Polymorphisms in the IL1RL1–IL18R1 Region Confer Risk for Ocular Behçet’s Disease in a Chinese Han Population

Single nucleotide polymorphisms (SNPs) in the IL1RL1–IL18R1 region are associated with various immune-mediated diseases. This study was carried out to investigate the causal variant for ocular Behçet’s disease (BD) and elucidate its target genes in the IL1RL1–IL18R1 region. Nine candidate functional SNPs were prioritized with bioinformatics analysis, followed by a two-stage association study in 694 ocular BD patients and 1,458 unaffected controls. Functional studies were performed in the peripheral blood mononuclear cells (PBMCs) of 45 healthy men and 16 active male BD patients. Genotyping was performed using the MassARRAY System. The mRNA expressions of IL1RL1, IL18R1, IL18RAP, and SLC9A4 were assayed by real-time PCR and secretion of cytokines was examined by ELISA. Significantly lower frequencies of the rs12987977 GG genotype/G allele (Pc = 8.93 × 10–7, OR = 0.39; Pc = 2.60 × 10–3, OR = 0.77, respectively), rs12999364 TT genotype/T allele (Pc = 3.15 × 10–4, OR = 0.51; Pc = 1.13 × 10–2, OR = 0.80, respectively), and rs4851569 AA genotype/A allele (Pc = 3.29 × 10–4, OR = 0.52; Pc = 9.72 × 10–3, OR = 0.80, respectively) were observed in BD patients compared with the controls. Functional experiments revealed a downregulation of IL1RL1, IL18R1, and SLC9A4 and a decreased secretion of IFN-γ in the anti-CD3/CD28 antibody-treated PBMCs as well as a decreased production of TNF-α in the lipopolysaccharide (LPS)-stimulated PBMCs in carriers of the protective homozygous rs12987977/GG genotype compared with the TT genotype. Our findings show that functional SNPs—rs12987977, rs12999364, and rs4851569—in the IL1RL1–IL18R1 region confer susceptibility to ocular BD in a Chinese Han population. And IL1RL1, IL18R1, and SLC9A4 may be the target genes of rs12987977.


INTRODUCTION
Behçet's disease (BD) is a systemic autoinflammatory/autoimmune condition with multi-tissue manifestations, such as oral aphthae, genital ulcers, arthritis, and uveitis. Being one of the most common uveitis entities in China, it mainly affects young adults and has a chronic relapsingremitting course (Yang et al., 2005). The exact etiology of BD remains largely unknown and is thought to be multifactorial, involving genetic, immunologic, and environmental factors (de Chambrun et al., 2012).
Genetic association studies are important since they may identify the pathological mechanisms causing BD and can be instrumental in the development of novel therapies. To date, HLA-B51 has been confirmed as the strongest genetic susceptibility factor for BD by several genome-wide association studies (GWAS) (Gul and Ohno, 2012;Kirino et al., 2013), although an increasing number of non-human leukocyte antigen (HLA) susceptibility variants have also been shown to be associated with this disease (Remmers et al., 2010;Lee et al., 2013). However, most susceptibility single nucleotide polymorphisms (SNPs) identified in GWAS were non-coding variants, and because of the incomplete coverage of SNPs in most GWAS and the phenomenon of co-inheritance of many variants, known as linkage disequilibrium (LD), functional studies addressing biological mechanisms which are aimed at pinning down the causal genetic variants and their gene targets in a given locus are still scarce (Gallagher and Chen-Plotkin, 2018). Fortunately, with the advent of public functional genomic databases such as GTEx (The GTEx Consortium, 2013), ENCODE (ENCODE Project Consortium, 2012), and Haploreg (Ward and Kellis, 2016), emerging studies are now reporting on this issue and promising results are obtained concerning the biological mechanisms that link susceptibility genetic variants with disease phenotypes (Allen et al., 2017;Cui et al., 2018).
An ongoing GWAS from our group suggested that the IL1RL1-IL18R1 region might harbor susceptibility loci for BD, although the p-values were below our GWAS threshold of 5 × 10 −8 (unpublished data). It has been reported that subthreshold SNPs could represent true disease risk loci that may have regulatory consequences and that an additional analysis of such SNPs may help to explain "missing heritability" when using a too stringent threshold (Wang et al., 2016). In addition, this region contains a cluster of immune-related genes that belong to the interleukin 1 receptor family (IL1RL2, IL1RL1, IL18R1, and IL18RAP). SNPs located in these genes were reported to be associated with various autoimmune or allergic diseases in earlier studies, such as inflammatory bowel disease and asthma (Liu et al., 2015;Zhu et al., 2018). Studies focusing on the biological function of the genes in this region also showed an important role in regulating the immune system. For instance, genetic ablation of IL1RL1 protected mice from developing experimental inflammatory bowel disease, and an antagonist of IL1RL1 could reduce the symptoms of colitis in these animals (Sedhom et al., 2013). The development of experimental autoimmune encephalomyelitis was attenuated in IL18R1-deficient mice (Gutcher et al., 2006). Taken together, this evidence highlights a crucial and pleiotropic role of this region in the development of immune-mediated diseases. The functional causal SNP and the exact susceptibility gene have, however, not yet been determined; therefore, despite the weak association in our GWAS, we decided to study this region in more detail in a large cohort of ocular BD patients and unaffected controls and also performed several functional experiments in peripheral blood mononuclear cells (PBMCs) from genotyped active BD patients and controls.

Ethical Considerations
This study was approved by the Clinical Ethical Committee of the First Affiliated Hospital of Chongqing Medical University (permit no. 2009-201008). The tenets of the Declaration of Helsinki were strictly followed in all procedures. Informed consent was collected from all subjects.

Subjects
A total of 694 ocular BD patients and 1,458 ethnically matched unaffected Chinese Han volunteers independent of our ongoing GWAS cohorts (unpublished data) were recruited from the uveitis clinic in the Ophthalmology Department of the First Affiliated Hospital of Chongqing Medical University (Chongqing, China) (Yang et al., 2018). Patients were diagnosed following the criteria for the diagnosis of Behçet's disease from the International Study Group for Behçet's Disease (International Study Group for Behçet's Disease, 1990). For the two-stage association study, the subjects were enrolled from 2008 to 2019. Another 45 healthy male volunteers and 16 active male ocular BD patients of Han ethnicity were also included for subsequent functional studies.

Prioritization of Candidate Functional SNPs
The regional plot of GWAS susceptibility loci was generated by the Locuszoom website (LOCUSZOOM, RRID:SCR_009257). The plot of pairwise LD matrix analysis was created using LDlink 1 . Candidate functional SNPs that were in strong LD with the chosen index SNPs and their epigenomic annotations and transcription factor (TF) binding information were obtained via the Haploreg database (HaploReg, RRID:SCR_006796). Candidate target genes (protein coding) of functional SNPs with FDR ≤ 0.05 for SNP-gene pairs were queried from the GTEx portal (Genotype-Tissue Expression, RRID:SCR_013042) and 3DSNP databases 2 (Lu et al., 2017).

DNA Extraction and Genotyping
DNA was extracted from venous blood with the QIAamp DNA Blood Mini Kit (QIAGEN, Valencia, CA, United States) according to the manufacturer's instructions. SNPs were genotyped with the MassARRAY system (Sequenom Inc., San Diego, CA, United States). The call rates of the SNPs tested in our study in the cases and controls were all above 95%.

Cell Isolation and Culture
The PBMCs of 45 healthy male volunteers were isolated from fresh peripheral blood by Ficoll-Hypaque density gradient centrifugation, then cultured in 24-well plates with complete RPMI 1640 medium (consisting of 10% fetal bovine serum, 100 U/ml penicillin, and 100 µg/ml streptomycin) at a density of 2 × 10 6 cells per well. The PBMCs of each individual were treated with 100 ng/ml lipopolysaccharide (LPS) (Sigma, MO, United States) for 1 day or a combination of anti-CD3 and anti-CD28 antibodies (5:1) (Miltenyi Biotec, Palo Alto, CA, United States) for 3 days, respectively, in an incubator with 5% CO 2 at 37 • C.

Enzyme-Linked Immunosorbent Assay
The concentrations of IL-1β, TNF-α, and IL-6 in the culture supernatants of the LPS-stimulated PBMCs as well as IFN-γ, IL-10, and IL-17 in the anti-CD3/CD28 antibody-stimulated PBMCs were quantified with the human Duoset enzyme-linked immunosorbent assay (ELISA) development kit (R&D Systems, Minneapolis, MN, United States) according to the instructions of the manufacturer. The representative standard curve of ELISA is shown in Supplementary Figure S5.

Statistical Analysis
For genetic association analysis, the genotype and allele frequency data were analyzed using Typer4.0 software from the MassARRAY system. The Hardy-Weinberg equilibrium (HWE) of all tested SNPs in the controls was performed using the SHEsis online tool (SHEsis: Analysis Tools For Random Samples, RRID:SCR_002958) (Shi and He, 2005). Statistical power of sample size was calculated with the online tool of power and sample size calculator 3 . χ 2 test, P-value, odds ratio (OR), 95% confidence intervals (95% CIs), and age-and sex-adjusted multivariate logistic regression analysis were calculated by SPSS (SPSS, RRID:SCR_002865). The P c values were calculated from the original P-values in the χ 2 test and logistic regression analysis with the Bonferroni correction in order to adjust for multiple comparisons. A correction factor of 9 was chosen, representing the number of SNPs analyzed in our study. The threshold of statistical significance of the P c values was set as 0.05. For in vitro functional analysis, unpaired t-test or non-parametric Mann-Whitney test was used for comparisons between two independent groups, while one-way ANOVA or non-parametric Kruskal-Wallis test was chosen for comparisons among three independent groups according to different applicable conditions using GraphPad Prism 8 (GraphPad Prism, RRID:SCR_002798). 3 http://powerandsamplesize.com/ Statistical significance was recognized as P-values lower than 0.05 (two-tailed).

Data Availability Statement
The raw data supporting the conclusions of this manuscript will be made available by the authors, without undue reservation, to any qualified researcher upon request.

Characteristics of the Study Subjects
The gender, age, and clinical characteristics of the 694 ocular BD patients and 1,458 unaffected controls are shown in Tables 1, 2. All BD patients had uveitis. The BD group consisted of 585 males (84.3%) and 109 females (15.7%), whereas 746 males (51.2%) and 712 (48.8%) females compose the control group. The proportion of females in the BD group is smaller than that of the control group (P < 0.05). The mean age of the BD group and the control group is 34.3 and 39.7 years, respectively (P < 0.05). Given the difference between these groups in terms of gender and age, we used multivariate logistic regression analysis to adjust for possible confounding effects.

Bioinformatics Analysis of Susceptibility Loci and Prioritization of Candidate Functional SNPs
An ongoing Behçet's disease GWAS in our laboratory identified a potential susceptibility region on chromosome 2q12, containing 25 non-coding SNPs with P-values below the GWAS threshold, but with a suggestive level of significance (4.63 × 10 −5 < P < 9.88 × 10 −4 ). These SNPs cover about 0.5 Mb on the chromosome and show a varying degree of LD with each other (Figure 1). They are located close to or in the genes encoding IL1RL2, ILIRL1, IL18R1, and IL18RAP (Figures 1, 2). Pairwise LD analysis revealed that rs2160202 and rs1420106 could capture almost all the other SNPs in this locus, and they were therefore chosen as the index SNPs for further bioinformatics analysis (Figure 2). The r 2 between these two SNPs is 0.481. It is now becoming apparent that causal variants that were thought to confer disease risk are not necessarily the tag variants detected by the microarrays used in GWAS, but variants that are in LD with the tag SNPs (Oldridge et al., 2015;Painter et al., 2016). Causal SNPs often reside in regulatory elements (promoter and enhancer) and could bind to certain TFs (Painter et al., 2016). They influence the disease susceptibility partly by regulating the expressions of its target genes in disease-relevant tissues and cell types (Farh et al., 2015). In addition, by long-range chromosome interactions, causal variants could even regulate distal gene targets. Thus, all SNPs that were in strong LD (r 2 > 0.8) with these two index SNPs were analyzed using genomic structure data from the Han Chinese in Beijing (CHB) and Southern Han Chinese (CHS) databank 4 . Candidate functional SNPs were prioritized according to the following criteria: (1) whether they were expression quantitative trait loci (eQTLs) in the GTEx portal; (2) whether they had three-dimensional interacting genes in the 3DSNP database; (3) whether they were located in regulatory 4 https://www.ncbi.nlm.nih.gov/variation/tools/1000genomes/ elements (promoter/enhancer) in immune cells with typical histone acetylation markers such as H3K27ac, H3K4Me1, and H3K4Me3; and (4) whether they were occupied by TFs or could change the binding motif of certain TFs in the Haploreg database. With the criteria mentioned above, we prioritized nine candidate functional SNPs for further association studies (Table 3).

Two-Stage Association Study
To confirm the association with the ocular BD phenotype and fine-map the true genetic signal in the IL1RL1-IL18R1 region, we next performed a two-stage association study for the nine candidate functional SNPs with 694 BD patients and 1,458 unaffected controls that were independent of our GWAS groups. None of the tested SNPs deviated from the HWE. A χ 2 test in the first-stage analysis, where 434 BD patients and 558 controls were included, showed significantly decreased frequencies of the rs12987977 (NC_000002.11:g.102975336T > G) GG genotype/G allele (P c = 2.30 × 10 −5 , OR = 0.32; P c = 2.88 × 10 −4 , OR = 0.66, respectively), rs12999364 (NC_000002.11:g.102974129C > T) TT genotype/T allele (P c = 4.76 × 10 −3 , OR = 0.46; P c = 2.30 × 10 −3 , OR = 0.70, respectively), and rs4851569 (NC_000002.11:g.102983247C > A) AA genotype/A allele (P c = 3.14 × 10 −3 , OR = 0.46; P c = 2.12 × 10 −3 , OR = 0.70, respectively) in the BD group compared to the controls ( Table 4). We were not able to detect an association with ocular BD for the other six SNPs (Supplementary Table S2), and these SNPs were not investigated in our further studies.
FIGURE 1 | Regional plot of nominal susceptibility SNPs identified in our ongoing genome-wide association study (GWAS) in the IL1RL1-IL18R1 region.
Frontiers in Genetics | www.frontiersin.org To expand the sample size in the first-stage study, another independent cohort consisting of 260 ocular BD patients and 900 unaffected controls was added to the association study. The results in the second-stage study also showed suggestive lower frequencies of the rs12987977/GG genotype (P = 9.20 × 10 −3 , OR = 0.51), rs12999364/TT genotype (P = 2.82 × 10 −2 , OR = 0.59), and the rs4851569/AA genotype (P = 3.22 × 10 −2 , OR = 0.60) before multiple corrections in BD compared to the controls ( Table 4). After combining the two stages together, it was confirmed that BD patients had lower frequencies of the rs12987977 GG genotype/G allele (P c = 8.93 × 10 −7 , OR = 0.39; P c = 2.60 × 10 −3 , OR = 0.77, respectively), rs12999364 TT 3 | Prioritization of nine candidate functional variants in susceptibility to Behçet's disease (BD) in the IL1RL1-IL18R1 region.
To adjust for possible confounding effects of gender and age, a multivariate logistic regression analysis of various genetic models was performed for the three susceptibility SNPs. The results of the co-dominant model were in agreement with the previous χ 2 tests. In addition, it also suggested that the associations of the rs12999364-T allele, rs12987977-G allele, and the rs4851569-A allele with BD behaved as recessive models, respectively ( Table 5).
Effect of rs12987977 on IL1RL1, IL18R1, IL18RAP, and SLC9A4 Expression rs12987977 showed the strongest association with BD and was therefore chosen for further investigations on possible causal biological mechanisms. Previous bioinformatics analysis with the 3DSNP database implied that rs12987977 could interact with IL1RL1, IL18R1, and IL18RAP in the 3D chromosome. Analysis with the GTEx portal tool also indicated that it correlated with the expressions of IL1RL1, IL18R1, IL18RAP, and SLC9A4 in different cell types and tissues. The GTEx portal lacks data for PBMCs or any other specific immune cells that are considered relevant to autoimmune disease. Donors in the GTEx project were of various ancestries. To investigate the target gene(s) of rs12987977 that might play an important role in the development of BD in Chinese Han, we performed a real-time quantitative PCR experiment to determine the effect of the different genotypes of rs12987977 on the expressions of those candidate target genes in the PBMCs obtained from a newly recruited group of genotyped healthy individuals. To rule out a possible confounding effect of gender, all the 45 volunteers enrolled for these functional studies are men. The role of various stimuli was tested and included unstimulated control, LPS-stimulated (non-specific inflammation signal), and anti-CD3/CD28 antibody-stimulated (antigen-mimicking signal). The results showed that, in unstimulated cells, there was a slightly significant downregulation of the mRNA level of SLC9A4 in rs12987977/GG carriers compared with TT carriers (P = 0.044; Figure 3A). There was no significant difference between the different genotypes for the expressions of the other three genes in the unstimulated controls (Supplementary Figures S1A-C). When the PBMCs were stimulated with the anti-CD3/CD28 antibody, the expressions of three out of four genes tested showed a significant difference between the different genotypes ( Figures 3B-D), while no significant difference was seen between the different genotypes and the expression of IL18RAP (Supplementary Figure S1D). A significantly lower expression of IL1RL1 in rs12987977/GG carriers was observed as compared to TT carriers (P = 0.021; Figure 3B), and a significantly lower expression of IL18R1 in rs12987977/GG carriers in contrast to the GT/TT carriers ( * P = 0.022 and * * P = 0.002, respectively; Figure 3C). Similar to the findings in non-stimulated PBMCs, the expression of SLC9A4 in rs12987977/GG carriers was significantly decreased compared with the TT carriers following the anti-CD3/CD28 antibody stimulation (P = 0.018; Figure 3D). When stimulated with LPS, none of the genes tested showed a significant effect of genotype on the mRNA expression (Supplementary Figures S1E-H).

Differential Expression Analysis Between Active BD Patients and Healthy Controls
To investigate whether the expressions of the target genes of rs12987977 were dependent on the disease state, we performed a real-time quantitative PCR assay in the PBMCs of 16 active BD patients and 45 healthy controls. All the patients and controls tested are men and individuals of the same genotype were compared. No differences were observed between the patients and the controls (Supplementary Figures S2A-H).

Effect of rs12987977 on Cytokine Production
To investigate whether the functional SNP could indirectly affect disease risk by exerting an influence on the release of proinflammatory or anti-inflammatory cytokines that had been wellstudied in autoimmune diseases, we tested cytokine secretion by PBMCs from healthy male volunteers with various genotypes of rs12987977. We failed to collect the cell culture supernatant of one subject by accident. Individuals carrying the TT genotype had an elevated secretion of IFN-γ (P = 0.002; Figure 4A) when stimulated with anti-CD3/CD28 and an increased TNF-α (P = 0.021; Figure 4B) production when stimulated with LPS as compared to those with the GG genotype. No effect of genotype was observed on the secretion of IL-17, IL-10, IL-1β, or IL-6 ( Supplementary Figures S3A-D).

DISCUSSION
In the current study, we identified nine candidate functional variants in the IL1RL1-IL18R1 region on chromosome 2 using public bioinformatics tools. Significant associations with ocular BD were found in three of these SNPs in a combined two-stage association study, including rs12987977, rs12999364, and rs4851569. The second-stage study by itself lacked statistically significant results, which may be due to the fact that the sample size of BD patients in the second stage was too small to provide sufficient statistical power. These three SNPs are novel susceptibility variants for BD. Since rs12987977 showed the strongest association with BD, we focused on this genetic variant to study the effect of genotype on the causal biological mechanisms. In vitro allelespecific expression analysis showed that various genotypes of the susceptibility SNP rs12987977 were correlated with the mRNA expressions of IL1RL1, IL18R1, and SLC9A4. Significant downregulation of the SLC9A4 mRNA level was observed in unstimulated PBMCs along with IL1RL1 and SLC9A4 in the anti-CD3/CD28 antibody-stimulated PBMCs of rs12987977/GG carriers compared with the TT carriers. The results also showed a significant downregulation of IL18R1 in rs12987977/GG carriers compared with the GT and TT carriers. Furthermore, cytokine production analysis between various genotypes of rs12987977 indicated significantly lower secretions of IFN-γ in the anti-CD3/CD28 antibody-stimulated PBMCs and TNF-α in the LPS-stimulated PBMCs of rs12987977/GG carriers compared with the TT carriers. These findings are in agreement with the fact that the TT genotype is the risk genotype for developing BD. The SNP rs4851569 discovered in our study is a novel genetic variant, and to our knowledge, no disease association has been reported for this SNP. SNP rs12987977 has been reported to be associated with both Crohn's disease and ulcerative colitis in a European population (Ellinghaus et al., 2016). The BD susceptibility variant rs12999364 found in our study was also shown to confer risk of asthma in a Dutch cohort (Reijmerink, 2008). Two SNPs (rs1420098 and rs1420101) that are in strong LD (r 2 > 0.8) with rs12987977 were reported to be susceptibility factors for IBD and asthma in other studies (Liu et al., 2015;Demenais et al., 2018). The evidence indicates a pleiotropic role for this region in the development of immune-mediated diseases. Chromosome 2q12, which contains the IL1RL1-IL18R1 region, harbors several genes of the interleukin 1 receptor family (IL1RL2, IL1RL1, IL18R1, and IL18RAP). IL1RL1 (also known as ST2) was identified as a target gene of rs12987977 in our study and is a receptor of interleukin-33 (IL-33) (De la Fuente et al., 2015). It is expressed in various cell types, including dendritic cells, T helper 1 (Th1) cells, T helper 2 (Th2) cells, CD8 + T cells, NK cells, B cells, and neutrophils (De la Fuente et al., 2015;Cayrol and Girard, 2018). The function of the IL-33/ST2 axis has been investigated in asthma, suggesting an important role in FIGURE 3 | Significant differential expressions of candidate target genes between various genotypes of rs12987977 in non-stimulated peripheral blood mononuclear cells (PBMCs) (A) and anti-CD3/CD28-stimulated PBMCs (B-D) from healthy male controls (GG = 10, GT = 17, TT = 18). Data are shown as the mean ± SEM. * P < 0.05; * * P < 0.01.
allergic disease (Grotenboer et al., 2013;Snelgrove et al., 2014). Nevertheless, emerging evidence suggests that this signaling pathway also plays an important role in chronic inflammatory diseases (Verri et al., 2010;Li et al., 2015). Increased expressions of IL-33 and ST2 have been observed in skin lesions of vitiligo patients . Administration of a blocking anti-ST2 antibody at the onset of collagen-induced arthritis attenuated the severity of the disease in mice and reduced joint destruction and was associated with a marked decrease in IFN-γ production by draining lymph node cells (Palmer et al., 2009). In a mouse model of methylated bovine serum albumin (mBSA)-induced arthritis, suppression of IL-33R (IL1RL1) expression in neutrophils could prevent IL-33-mediated neutrophil migration into the knee joint (Verri et al., 2010). Besides IL1RL1, we showed that rs12987977 also targets IL18R1. The latter forms a receptor complex with IL18RAP, responsible for the binding of the pro-inflammatory cytokine IL-18 (Yasuda et al., 2019). It can mediate IFN-γ synthesis from Th1 cells and contributes to Th1 differentiation (Tominaga et al., 2000). An increased IL18R1 mRNA expression has been shown in both cerebrospinal fluid cells and PBMCs in multiple sclerosis (Gillett et al., 2010). Blocking IL18R1 in MRL-Fas lpr mice, which develop a systemic lupus erythematosuslike disease, results in significant reductions in renal pathology and the IFN-γ and TNF-α mRNA levels (Kinoshita et al., 2004). Further studies showed that deletion of IL-18R1 in intestinal epithelial cells inhibited mucosal damage and colitis in a mouse model of ulcerative colitis (Nowarski et al., 2015). The third gene that we identified as a target for rs12987977 is SLC9A4 (Solute Carrier Family 9 Member A4). SLC9A4 is a sodium/hydrogen exchanger involved in intracellular pH regulation and cell volume (Beltrán et al., 2008;Arena et al., 2012). It also functions as a transporter of glucose, bile salts, organic acids, and amine compounds (Beltrán et al., 2008;Arena et al., 2012). Its role in the pathogenesis of immune-mediated diseases has, to our knowledge, not yet been reported, although recent studies have pointed out a role of basic molecules in T cell differentiation (Mcgaha et al., 2012). Taken together, the evidence suggests a pro-inflammatory role for IL1RL1 and IL18R1 in autoimmune disease. The genetic control of these genes by rs12987977 is supported by the observation that individuals with the protective homozygous GG genotype of rs12987977 had significantly lower IL1RL1 and IL18R1 mRNA expressions and decreased secretions of IFN-γ and TNF-α from their PBMCs treated with certain extracellular stimuli. This is in line with the hypothesis that certain environmental factors are needed to trigger the development of BD in a genetically predisposed individual, whereby a multiplicity of genes, all leading to a proinflammatory state, may eventually lead to a full-blown disease. Until now, we were not able to detect differential expressions for IL1RL1, IL18R1, and SLC9A4 when comparing active BD patients and controls carrying the same genotype. There are a number of explanations for this observation. We investigated whether the target genes of rs12987977 (IL1RL1, IL18R1, and SLC9A4) were differentially expressed between BD and controls using 16 active BD patients and 45 controls. For the controls, we planned to recruit at least 10 cases for each of the three genotypes. After we collected 45 controls and genotyped them, we found that we had sufficient samples from each genotype. However, it was much more difficult to collect the same number of active BD patients during the same time period because most BD patients had already consulted other hospitals and had the inflammation under control before they came to us and could therefore not be included in the "active" BD group. The unequal sample size may therefore be one reason, and another possible reason is that most active BD patients were under immunosuppressive medication, which might affect the expressions of these genes. It has been reported that patients with BD have a lower incidence of allergic diseases, such as atopic dermatitis, allergic rhinitis, bronchial asthma, and food/drug allergies, than the entire population in Japan (Horie et al., 2016). In our study, none of the controls had allergic disease and we only recorded a self-reported atopic dermatitis in 89 of our BD patients. A stratification analysis in BD cases with atopic dermatitis did not reveal a significant genetic correlation, which may be due to the fact that the numbers were too small to perform a meaningful genetic association analysis (data not shown). Future research with more complete information concerning the presence of allergic diseases in BD is required to adequately address this issue.
In our current study, we only used GTEx and 3DSNP to find the candidate target genes of rs12987977 and focused on protein-coding genes. Other databases may also provide useful information. For instance, rs12987977 has been shown to be associated with the expressions of AC007278.3, AC007278.2, IL1RL1, IL18R1, and MFSD9 in the whole blood and/or PBMCs of Europeans (see website of eQTLGen 5 ) (Võsa et al., 2018). In future studies, we would like to investigate whether rs12987977 is also associated with the expressions of other genes (including non-coding genes) using more databases.
How an intron variant of IL18R1 such as rs12987977 affects the expressions of the target genes mentioned above is not yet clear. Analysis with bioinformatics tools suggests that it resides in the enhancer region of various cells of the immune system, including T cells, B cells, monocytes, natural killer cells, and neutrophils, and could bind to a well-studied transcription factor named CTCF that has been reported to be instrumental in chromosome looping (Tang et al., 2015;Nora et al., 2017). We suggest that the susceptibility variant rs12987977 may regulate the expressions of these gene targets by possibly changing the activity of the enhancer region, affecting chromosome interactions, or by altering the affinity of binding to CTCF when shifted from the T allele to the G allele. However, more functional studies should be carried out to clarify these possible mechanisms.
An earlier study of our group has identified rs2058660 as a susceptibility genetic variant for ocular BD in a Chinese Han cohort with G as the risk allele (P c = 5.94 × 10 −6 , OR = 1.33) in a total of 1,063 patients and 1,872 controls that did not exclude the GWAS cohorts (Tan et al., 2016). In the current study, this SNP showed a weak association with BD whereby the G allele was also identified as the risk factor (P = 0.047, OR = 1.198) in the firststage association study, but the significance disappeared after Bonferroni correction for multiple comparisons. The samples used in the first stage of these two studies are different, and the sample size was also not the same. The difference in sampling could partly account for this discrepancy. Therefore, we cannot entirely exclude the possibility that the SNP rs2058660 may also be associated with ocular BD when this is tested in a larger cohort.
There are several limitations to our study. Firstly, our BD patients were recruited from an ophthalmology department and all had uveitis. However, not all patients with BD necessarily have uveitis, and they may consult different medical departments depending on their main complaints. It would therefore be very interesting to see whether our findings concerning the associations found in the IL1RL1-IL18R1 region confer susceptibility to Behçet's disease itself or to the development of ocular lesions in this disease. Future studies including BD patients from other medical departments (e.g., rheumatology) and other subpopulations are required to address this issue. Secondly, our results showed a correlation between various genotypes of rs12987977 and the expressions of candidate susceptibility genes and predicted the biological mechanisms underlying this correlation with bioinformatic methods, but further confirmation studies using enhancer reporter assays, CHIP-seq, and chromosome conformation capture should be carried out to elucidate the exact mechanism. Third, due to economic reasons, we did not expand the sample size for those SNPs that showed no association with BD after stringent corrections in the first-stage study. Larger cohorts are needed to study the association with BD for these SNPs in future studies. Lastly, we used primers for total isoforms for each gene instead of separate splice variants. More detailed studies regarding certain gene isoforms should be performed in the future to help us understand the function of the various gene isoforms.

CONCLUSION
Our findings show that several functional SNPs-rs12987977, rs12999364, and rs4851569-residing in the IL1RL1-IL18R1 region confer susceptibility to ocular BD in a Chinese Han population. We focused on rs12987977 and showed that it targets IL1RL1, IL18R1, and SLC9A4. Further studies are needed to reveal the exact mechanism whereby this variant affects the expressions of these target genes and how this eventually may lead to the development of ocular BD.

DATA AVAILABILITY STATEMENT
The associations of the three susceptibility SNPs with Behcet's disease from our manuscript are available at the ClinVar database of NCBI (accession numbers: SCV001334435 -SCV001334437).

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by the Clinical Ethical Committee of the First Affiliated Hospital of Chongqing Medical University (permit no. 2009-201008). Written informed consent to participate in this study was provided by the participants' legal guardian/next of kin.

AUTHOR CONTRIBUTIONS
XT conceptualized and designed the study, carried out the analyses, drafted the initial manuscript, and revised the manuscript. QZ and ML designed the study, collected data, and reviewed the manuscript. HT, QW, and LZ reviewed the manuscript for important intellectual content. QC and GY collected samples and reviewed the manuscript. GS, AK, and PY reviewed the manuscript for important intellectual content and revised it. All authors approved the final manuscript as submitted and agreed to be accountable for all aspects of the work.