Introduction

Behçet’s disease (BD) is a complex and systemic inflammatory syndrome; from a clinical point of view, it is a vasculitis that presents heterogeneous phenotypes. Genital and oral ulcers are the first and most common manifestations, although bilateral uveitis and cutaneous involvement are also present with a relatively high frequency, and there are occasional manifestations in other organs or systems, such as the nervous system1,2,3. Gender distribution is variable among countries, but unlike other immune-mediated diseases there is not a predominance of females among affected subjects4. BD is, in general, a rare disease although its prevalence is higher along the old “Silk Road”, especially in Turkey (up to 370 cases per 100000 population)4. The pathogenic mechanisms of BD remain unknown, although there are accumulated evidences of genetic predisposition and data suggesting that an aberrant immune response to certain infectious agents and some environmental factors may trigger the disease5. The first described genetic factor associated with the disease was HLA-B*516, which is present in up to 60% of patients. More recently, genome wide association studies (GWAS) reinforce HLA-B*51 as the most important genetic factor in susceptibility to this disease and, in addition, they support the contribution of other genes related with the immune system, such as IL23R, IL10, STAT4, CCR1-CCR3, KLRC4, ERAP1, TNFAIP3 and FUT27.

There is a considerable clinical overlap of BD with autoinflammatory syndromes (AID) which are characterised by a misregulation of the inflammasome proteins that lead to an increased release of IL-1b. AID are genetically conditioned diseases, caused by aberrant proteins encoded by pathogenic variants of genes related with both the inflammasome itself (for example, NLRP3) as well as with other sensing proteins involved in pathogen- and damage-associated molecular patterns (PAMPs and DAMPs) or in the signalling pathways8. Heat shock proteins (HSP), are a group of intracellular proteins with scavenger roles found from bacteria to mammals. Human NF-M is homolog of bacterial HSP-65, and BD patients’ sera show crossreativity with these two proteins, suggesting that bacterial HSP could be recognised by TLRs and act as a trigger antigen for BD development9,10. Based on the above, it is plausible that variants in genes involved in AID could also be related to BD. In this sense, in a previous study, which included CECR1, MEFV, MVK, NLRP3, NOD2, PSTPIP1 and TNFRSF1A, we reported that many rare variants in these AID related genes are found in BD patients11.

Gene-gene interactions are assumed to be a common factor in the genetic components of human diseases12, therefore, it is interesting to establish possible relationships among the different factors involved in the susceptibility to diseases. For instance, interaction between HLA and other genes associated with diseases has been demonstrated in BD and other immuno-mediated syndromes related to HLA class I, such as ankylosing spondylitis and psoriasis7,13,14. Our previous study with AID related genes was performed using next generation sequencing (NGS) with a pooled strategy which did not allow us to assign each variant to an individual. The aim of the present study was to determine the number of patients carrying variants, the proportion of individuals that accumulate variants and the relationship with other genetic factors, such as HLA-B*51.

Materials and Methods

Subjects

The study population consisted of 355 patients diagnosed with BD according the 1990 International Study Group classification criteria for this disease15,16. All the patients had Spanish European ancestry16, 43.7% were male and 36% were B51 positive. The clinical features of this cohort have been previously published11. The study was approved by the ethical committees of all centers and hospitals involved (Hospital Universitario Virgen del Rocío, Hospital Clínico San Cecilio, Hospital Universitari Clínic, Complejo Hospitalario Universitario A Coruña, Hospital Universitario de Valme, Hospital Universitari Son Espases, Hospital Vall d’Hebron, Hospital Universitario Marqués de Valdecilla, Complejo Hospitalario Torrecárdenas, Hospital Universitario Central de Asturias, Hospital Virgen del Camino, Hospital Universitari Mútua Terrassa, Hospital Clínico San Carlos, Hospital Regional Universitario de Málaga, Hospital de la Princesa, Hospital Universitario Doctor Peset, Instituto de Parasitología y Biomedicina López-Neyra). The healthy control population in the rare variants’ analysis, consisted on the whole genome sequences of 107 individuals from the Iberian population of Spain (IBS population), included in the phase 3 of the 1000 Genomes Project; whereas in polymorphisms analysis 363 in-house healthy controls were used.

AID Related Genes Genotyping

In a previously published study with the same cohort, 54 rare variants (minor allele frequency, MAF < 0.05 in IBS population) were found in the 7 genes included in the study by using an NGS approach of pooled samples11. In the present study, in order to identify those patients carrying each variant, all the samples included in the pools which contained each rare variant were individually studied. 51 of the 54 rare variants were sequenced by Sanger and the remaining three (NLRP3 Val198Met: rs121908147, C__64676050_10, NOD2 Asn289Ser: rs5743271, C__26935007_10 and TNFRSF1A Arg121Gln: rs4149584, C__12029224_20) were TaqMan genotyped. In addition, to determine whether close variants (in the same exon) are in the same (cis) or in different (trans) chromosome, a Hexaprimer Amplification Refractory Mutation System PCR (H-ARMS PCR) was used17. The pairs of variants studied were: MEFV p.Leu110Pro and p.Glu148Gln, using the primers MEFV-H_ARMS_110_148-F1, R1, F2, R2, F3 and R3 and NOD2 p.Arg311Trp and p.Arg703Cys with the primers NOD2-H_ARMS_311_703-F1, F2, R1 and R2 (Supplementary Table 1). In the case of MEFV H-ARMS, the following thermal-cycler conditions were used: 3 min at 95 °C followed by 35 cycles of: 20 s at 95 °C, 30 s at 65 °C and 45 s at 72 °C. In this case, an informative band of 145 bp is present only when both variants are in trans (Supplementary Fig. 1A). Regarding NOD2, a modified H-ARMS was used, since the 2 variants are more distant, a two-PCR reactions system was chosen (one with the pair of primers F1 and R1 and another with the pair F2 and R2). The thermal-cycler conditions were: 20 s at 95 °C, 30 s at 70 °C and 90 s at 72 °C. In this case, an informative band of 1228 bp in one of the PCR tubes is present only when both variants are in trans (Supplementary Fig. 1B).

Regarding the three functional polymorphisms found in the previous study, NLRP3 p.Gln703Lys, NOD2 p.Arg702Trp and p.Val955Ile were individually genotyped in all the samples included in those pools in which the variants had been detected. The following TaqMan assays were used: NLRP3 p.Gln703Lys (rs35829419, C__25648615_10), NOD2 p.Arg702Trp (rs2066844, C__11717468_20) and p.Val955Ile (rs5743291, C__25651076_10). Real-Time PCR assays were performed using TaqPath™ ProAmp™ Master Mix (Thermo Fisher Scientific, Waltham, MA), in a 7500 Fast Real-Time PCR System, according to manufacturer’s recommendations. No statistically significant differences in the allelic frequencies of the polymorphisms were observed between NGS and TaqMan assay results. The remaining 4 polymorphisms (MAF > 0.05 in IBS population) found in our previous study were excluded from the present study because they do not have a functional relevance.

Regarding the controls, genotypes of rare variants were retrieved from the IBS population of the 1000 Genomes Project database and they were imported to PLINK software for genotype counting and association testing. In addition, in the cases of the functional polymorphisms NLRP3 p.Gln703Lys, NOD2 p.Arg702Trp and p.Val955Ile, 363 in-house individuals were genotyped using the same TaqMan assays as in patients and used as local ethnically matched controls.

HLA-B Low Resolution Genotyping

The locus HLA-B was genotyped in the 343 BD patients and the 363 local controls by using a standard low resolution (first set of two digits) PCR-SSOP Luminex method (LABType SSO One Lambda Inc., Canoga Park, CA). The remaining 12 patients could not be genotyped in HLA because of lack of material or its poor quality. Regarding the IBS controls, HLA-B genotype of IBS population of the 1000 genomes project was called at the same level of resolution as the BD patients using a pipeline optimised for OptiType software18. Briefly, chr6:29910247-31325022 hg19 region was retrieved from the mapped exome BAM files and converted to FASTQ with bamtools; the readings were aligned to all the available HLA allele sequences of IMGT/HLA with razerS3 and these alignments were analysed by OptiType software for HLA class I type calling. Data were visually inspected for proper coverage of HLA class I exons. The allelic frequency of HLA-B*51 by this method in IBS population (N = 107) was 0.09, which is not significantly different from that obtained in the general Spanish population with a standard HLA typing method (0.07)19.

Statistical analysis

To analyse the relationship between HLA-B*51 and rare variants in the AID genes, BD patients were divided into two groups: B51 positive (assigned as affected phenotype) and B51 negative (assigned as non-affected phenotype) and the regions studied were those defined by the field refGene.name2. SKAT test implemented by Variant Tools software20,21 was used to analyse the distribution of rare variants in each of the genes included in the study in both patient groups. Additionally, to detect interactions between rare variants of each pair of AID genes as well as between each AID genes and HLA-B*51, a recently developed functional regression model was used with the R package FRGEpistasis22. This method has been used for quantitative traits23, but it is also suitable for the binary ones. Briefly, the cases/controls were encoded as 1/0, the AID genes variants were codified as 0 (wild-type) or 1 (variant) and also the B51 status (1 positive, 0 negative). P-values were adjusted (Padj) according to Bonferroni’s (considering 8 tests: 7 AID genes and B51), therefore P-values below 0.0063 were considered statistically significant. Type I error was estimated using a simulation with data retrieved from the 1000 genomes database (Supplementary Table 2). In addition, FRG analysis was also done in our BD patients and IBS controls with randomly assigned the phenotypes (affected and unaffected) (Supplementary Table 3). Q-Q plots of the p-values obtained in FRG with both, the real and the random data are displayed in Supplementary Fig. 2.

In regards to the disease association study of the 3 functional polymorphisms: NLRP3 p.Gln703Lys, NOD2 p.Arg702Trp and p.Val955Ile, each sample was encoded as 0, 1 and 2 according to the number of copies of the minor allele of each polymorphism carried by the individual. Allelic and genotypic distributions in dominant (for minor allele) and allelic models in BD patients and controls were compared with Χ2 in two by two contingency Tables with 1-degree freedom using PLINK software and P-values were corrected by Bonferroni’s adjustment considering 12 tests, P values < 0.004 were considered statistically significant and the Odds Ratios (ORs) were calculated with fmsb package for R. Finally, to analyse interaction between B51 and these functional polymorphisms, conditional logistic regression models were constructed by categorizing all the individuals according to the presence/absence of the B51 (dominant model) and the genotype of each polymorphism (dominant model for the minor allele) and they were analysed with Epi Info 2002. Additionally, the R package FRGEpistasis22 was also used to detect interaction of these polymorphisms with B51.

To quantify linkage disequilibrium (LD) when it was relevant, the R values were calculated in the patients and local controls by using PLINK, whereas in the EUR population, LD between variants was consulted in HaploReg v.4.124. A descriptive circular layout figure, performed with circlize R package, was used to graphically represent the rare variants that are found together in same patient25.

Ethical use of human participants statement

All participants signed a written informed consent to publish identified information prior to their enrolment in the study. All methods were performed after obtaining the signed informed consent from patients, and in accordance to our institution ethical committee: Hospital Universitario Virgen del Rocio’s ethical committee, and with the Helsinki Declaration of 1975, as revised in 1983.

Results

The individual genotyping of the samples revealed that 98 patients out of the total of 355 in our cohort (27.6%) carried at least one of the total of 54 rare variants found in the genes included in the study. From these 98 patients, 82 harboured only one rare variant, 13 had two and 3 had three, in the same or different genes. All the patients with one variant had a single copy, except two homozygous (one NOD2 p.Gly908Arg and the other TNFRSF1A p.Arg121Gln). In regard to the 13 patients carrying two variants, in 3 cases the variants were located in the same gene and in 10 cases they were in different genes. Regarding the 3 patients with 3 variants, 2 of them had the combined allele MEFV p.Arg408Gln/p.Pro369Ser and one variant in a different gene, whereas the remaining patient had each variant in a different gene. Figure 1 illustrates connections between variants and Table 1 displays the list of combinations found in our cohort. In order to know whether variants in the same gene were in cis or in trans, PCR-H-ARMS was performed in the patients with MEFV p.Glu148Gln/p.Leu110Pro and with NOD2 p.Arg703Cys/p.Arg311Trp. In the case of MEFV, the variants were in cis whereas in NOD2 they were in trans (Supplementary Fig. 1A,B). No additional study was made in the 3 individuals with the combination MEFV p.Arg408Gln/p.Pro369Ser since these variants have been reported as being in a strong LD in the European control population (r2 = 1), therefore cis configuration was assumed.

Figure 1
figure 1

Descriptive circular layout showing coexistence of rare variants in patients. Length of sectors corresponds to number of variants; width of lines corresponds to the number of connections, each tick in inner circle axis represents one variant in one patient.

Full size image
Table 1 List of rare variants found in combination in the seven AID genes studied in our cohort of BD patients.
Full size table

The results of the SKAT test suggested non-significant differences in the distribution of rare variants between B51 positive and B51 negative patients (p > 0.05, Table 2). Nevertheless, functional regression model analysis, suggests epistatic interaction between HLA-B*51 and MEFV (p < 0.05, Table 2) in BD. No evidences of interaction between pairs of AID related genes were detected.

Table 2 Analysis of the distribution of rare variants in seven genes related to AID in B51 positive and negative BD patients using two different approaches.
Full size table

Table 3 displays the frequencies of the three functional polymorphisms in BD patients and healthy controls. A suggestive association was found in the case of NOD2 p.Arg702Trp (p = 0.011 in the allelic model and p = 0.014 in the dominant model: TT + CT vs. CC) with a protective effect of the minor allele. For the two remaining polymorphisms, no-significant association with BD was found (p > 0.05). No LD was detected between the two NOD2 polymorphisms included (R < 0.1 in all cases). With respect to the relationship of these polymorphism with B51, conditional regression analysis suggests that the protective effect of NOD2 p.Arg702Trp is independent of the HLA-B51, because the OR is decreased in both, B51 positive and negative groups (Table 4). The functional regression model analysis does not suggest epistatic interaction between B51 and these polymorphisms (p > 0.05 in all the cases).

Table 3 Distribution of frequencies of NLRP3 and NOD2 functional polymorphisms in BD cases and healthy controls.
Full size table
Table 4 Conditional logistic regression model of HLA-B51 and NOD2 p.Arg702Trp.
Full size table

Discussion

In the present study, a substantial percentage (about 30%) of the BD patients was found to have some rare variant in at least one out the seven AID related genes that were included. The proteins encoded by these AID genes are related with the IL-1b release, mainly through sensing bacterial or viral components and activating the NLRP3 inflammasome and, therefore, it would be reasonable to expect some degree of interaction among them. In this way, a relatively high ratio of our BD patients carried variants in different AID genes, suggesting that interplay between variants might exist. Hypothetically, these variants in AID genes could modify the association of B51 with BD or vice versa. It is necessary to note that, genetic interaction is difficult to detect because multilocus genotype combinations for gene-gene interaction increase exponentially and require a larger sample size as well as more computation burden. In addition, most of the methods developed to detect gene-gene interactions have been designed for common variants and it is difficult to apply them to data with rare variants because the low frequency of these variants leads to a lack of statistical power26,27. In spite of the use of a method designed for rare variants, no-interactions between AID genes were detected, although this result could be attributable to the low sample size, since only 13 patients accumulate variants in different genes. In addition, it is necessary to take into account that, from a biological point of view, gene-gene interaction is the phenomenon in which the effects of one gene on a trait (for example a disease) are modified by another gene or several other genes; whereas statistical epistasis is described as a deviation from additivity in a linear statistical model28. In order to assess interaction between rare variants of AID genes and HLA-B*51, two different approaches were used; in the first, patients B51 positive and negative were compared with SKAT whereas in the second, patients and controls were compared using FRG. The only discrepancy observed between both tests was in the case of MEFV in which, both tests have apparently contradictory results. According to SKAT, no differences in the distribution of rare variants were detected between B51 positive and negative individuals whereas our results with FRG suggest interaction between HLA-B*51 and MEFV. The estimated type I error in our simulation is even lower for FRG than for SKAT therefore, discrepancy could be based in the characteristics of each test, SKAT was designed to study association of rare variants in complex diseases by grouping variants and checking the association of each entire gene, whereas FRG was specifically designed to study gene-gene interactions with rare variants. Contradictory results in the association of MEFV with BD have been reported29,30,31,32, in fact, in our previous study in AID genes, MEFV was found associated only with SKAT but not with the rest of the statistical tests used11. To our acknowledgement, no previous work on interaction between HLA and MEFV has been published for BD. Our results must be interpreted with caution, nevertheless, there is some evidence suggesting that involvement of MEFV in autoimmune diseases could be associated with HLA molecules33,34,35 and conversely, differential effects of HLA class I molecules in Japanese patients with familial Mediterranean fever (FMF), even in those with high-penetrance MEFV mutations, have been identified36. MEFV protein interacts with cellular components of inflammasome whose activity is augmented by the pathogenic mutations37. Therefore, it is reasonable to consider that FMF-associated MEFV mutations enhance the cellular responses in the inflammatory diseases which are also associated with certain HLA polymorphisms. In other words, mutations and polymorphisms in the MEFV gene are associated with subclinical inflammation therefore, it is possible that the presence of MEFV mutations increase the baseline level of inflammation, and the presence of a certain HLA molecule (such as HLA-B*51 in BD) modulates the disease. Although the same type of interaction could be hypothesized for any of the other AID genes studied, no more interactions were detected.

Regarding the polymorphisms in NLRP3 and NOD2, only a protective association of NOD2 p.Arg702Trp was found. The loss of function variants of NOD2 have been reported as possibly protective in BD30,38. For this reason, in our previous study, the frequencies of NOD2 p.Asn289Ser, p.Arg702Trp and p.Gly908Arg in our patient cohort were compared with those of the 1000 genomes IBS control population11. Although no significant differences were found, a trend was detected for p.Arg702Trp (p = 0.09) that in the present study reaches statistical significance by increasing the number of controls by including local ethnically matched individuals. Therefore, the present study suggests a protective effect of this NOD2 variant at least in Caucasian patients. This effect seems to be independent of HLA-B*51.

In conclusion, a high percentage of patients with BD have one or several variants in AID genes. Additionally, although our data require replication, they suggest that the association of MEFV with BD could be modulated by the HLA molecules; whereas the protective effect of the loss of function of NOD2 p.Arg702Trp polymorphism would be independent of HLA.