2.1 SNPs selection

The DisGeNET (https://www.disgenet.org/) and GWAS Catalog (https://www.ebi.ac.uk/gwas/home) databases were used to identify the genetic biomarkers associated with IBD. The GWAS catalog has been used to find all SNPs associated with IBD. Non-synonymous missense variants showing association (p-value ≤5 × 10⁻⁸) with IBD risk in Caucasian populations were further filtered (Wijmenga and Zhernakova, 2018). The selected SNPs population frequencies were validated in Ensembl Human GRCh38.p13 genome (https://www.asia.ensembl.org/index.html), the Saudi Genome Project (SGP) database (https://www.sso.saudigenomeproject.com/), and the Greater Middle East (GME) Variome Project (GME Variome) (http://www.igm.ucsd.edu/gme/data-browser.php).

2.2 Subjects and Ethical approval

This research work was executed at the Department of Genetic Medicine, Faculty of Medicine, and PACER-HD labs of King Abdulaziz University. Ethical approval for this proposal was obtained through the Institutional Review Board (IRB) at King Abdulaziz University Hospital. The inclusion criteria for recruitment of IBD patients include confirmed clinical and pathological diagnosis of IBD, free from the family history of IBD, and Saudi nationality. All non-Saudi patients or Saudis with a positive family history of IBD were excluded. Healthy Saudis with no reported family history of IBD or any other autoimmune, gastrointestinal, or metabolic diseases were recruited as controls.

2.3 Sampling

IBD patients were recruited from pediatric and adult gastroenterology clinics based on their confirmed diagnosis, clinical presentations, endoscopic and histopathological characteristics, and gastrointestinal imaging. All study participants (97 IBD patients and 100 healthy controls) or their parents/guardians (for children less than 18 years of age) were thoroughly informed about the study’s aims and significance before recruitment and they signed the informed consent forms voluntarily. The clinical data, such as symptoms and age at the onset, gastrointestinal manifestations, and extra-intestinal manifestations, were also obtained from the hospital records of each participant.

2.4 Genomic DNA isolation

Two milliliters of peripheral blood were obtained from each participant in EDTA vacutainers and stored at −20°C until further processing was done. The DNA quality and quantity analysis was done by NanoDrop™ 2000c Spectrophotometer.

2.5 Genetic analysis

After extracting and quantifying DNA samples, Realtime PCR based TaqMan Allelic Discrimination Assay was carried out for genotyping 4 SNPs [rs76418789, rs11209026, rs2066842, rs2066844] using ABI7500HT Sequence Detection System (Applied Biosystems).

2.6 Validation

Sanger sequencing was done to verify the genotypes obtained by real-time PCR. Three samples from different genotype categories (homozygous major allele, homozygous minor allele and heterozygous alleles) of all four SNPs were sequenced using a 3500 Genetic Analyzer machine (Applied Biosystems). The PCR product sizes of NOD2 gene SNPs, i.e., rs2066842 and rs2066844, were 574 bp and 558 bp respectively, and IL23R gene SNPs, i.e., rs76418789 and rs11209026, were 390 bp, and 480 bp, respectively . The primers for the Polymerase Chain Reaction (PCR) for selected SNPs were designed using Primer3 (http://primer3.ut.ee/), and their specificity was validated using UCSC (https://genome.ucsc.edu/) and Blast (http://blast.ncbi.nlm.nih.gov/Blast.cgi) (.

2.7 Statistical analysis

GraphPad QuickCalcs, version 6.0 (GraphPad Software Inc., USA) online software was used to perform statistical analysis of the categorical (in percentage form) and continuous variable (in mean and standard deviation form) data. The statistical significance of allele and genotype frequency distributions was tested with the help of Pearson's standard chi-square test, odds ratio (OR), and 95% confidence interval (CI). Additionally, Hardy-Weinberg equilibrium was checked by the contingency table of observed versus predicted genotype distributions by using a chi-squared test with one degree of freedom. Epistasis between four examined polymorphisms and IBD occurrence was tested using the open-source multifactor dimensional reduction (MDR) program (MDR version 2.0 beta 5 http://www.epistasis.org/). The MDR (0.4.9 alpha) permutation test module performs more than 1000 permutation experiments. In this study, we compared the observed test accuracy with the predicted results to analyze the statistical significance. A p-value of <0.05 was considered as statistically significant.

2.8 Pathogenicity prediction

We used the variant effect predictor (VEP) webserver (http://asia.ensembl.org/Tools/VEP) hosted in Ensembldatabase to determine the variant pathogenicity with many computational prediction methods (McLaren et al., 2016). VEP accepts different input options including chromosome and nucleotide numbers, allelic variant, strand of the variant in human genome variation society (HGVS) identification format. From the VEP’s output, pathogenicity prediction scores of Sorting Intolerant from Tolerant (SIFT), Polyphen2, CADD, and Condel, were selected. The deleterious prediction (x) scoring criteria adopted by these different prediction tools is as follows; x < 0.05 for SIFT, x > 0.5 scores for PolyPhen2, x > 25 score for CADD and x > 0.5 score for Condel.

2.9 The 3-dimensional protein modelling and stability analysis

We have also explored the variant induced effects on structural features of proteins at molecular level. In this regard, to build the molecular 3D model of proteins we used I-Tasser webserver (Yang and Zhang 2015), which generates five models based on multiple-threading alignments and iterative template fragment assembly simulations. Out of these five models, the best model (selected based on the confidence (C) score, TM-score and RMSD values) was considered for downstream structural analysis. The selected models were further subjected to energy minimization to remove bad contacts in the protein structures by the SPDBV software. The energy minimized models of wild NOD2 and IL23R proteins have been used as a template to build the mutated NOD2 and IL23R structural models, followed by mutant protein stability analysis by the DUET web server (http://biosig.unimelb.edu.au/duet/) (Pires et al., 2014). The Pymol tool was used to evaluate the structural divergence between wildtype and mutant models of NOD2 and IL23R (Janson et al., 2017).

3.1 Clinical details

In this study, 97 patients (53 CD patients, 25 UC patients, and 19 indeterminate IBD) and 100 healthy control participants were recruited. The gender distribution of our IBD patients showed that there were more males (57.7%) than females (43.2%). The mean age of disease onset was 10.98 ± 7.39 years, and the average age was 16.84 ± 8.13 years. Most were diagnosed during their childhood (88.6%), and the rest (11.4%) in adulthood. Most of the CD patients had ileocolonic location (57%), and 16% had inflammation limited to the terminal ileum; but none of the participants had upper GI involvement. More than half (58%) of CD patients had a non-structuring and non-penetrating disease behavior, 28.6% structuring disease, and 13.4% penetrating disease behavior. In UC patients, (according to Paris classification) 33.8% had pancolitis, 33.8% Left-sided colitis, and about 32.4% with extensive UC. According to the Montreal classification, approximately 44.4% of UC patients had extensive UC, 38.3% ulcerative proctitis, and 17.3 % left-sided UC.

3.2 SNPs selection

Searching for the keywords “Inflammatory Bowel Disease” on the GWAS catalog database yielded 1150 SNPs. Using the DisGeNET database, we identified five genes or SNP with P-value (p-value ≤5 × 10⁻⁸). Concordance results between the DisGeNET and GWAS catalog resulted in 37 significantly associated missense SNPs with IBD. Out of 37 missense SNPs, four including two in NOD2 gene (rs2066842, rs2066844) and two in IL23R gene (rs76418789, rs11209026), which were previously reported in the Caucasian IBD patients were selected (Supplementary Table S1). The variant details like its chromosomal location, variant nature and minor allele frequency (MAF) were also obtained from Ensemble and SGP databases (Table 1).

3.3 Genetic analysis

The genotyping results of all four tested SNPs [rs2066842, rs2066844, rs11209026, rs76418789] were consistent with Hardy-Weinberg equilibrium in Saudi cohort. For rs2066842 (c.802C > T: p.P268S, NOD2) polymorphism, the distribution of the T allele in both monoallelic and biallelic forms is almost similar in IBD patients and control individuals (p values are > 0.05 for all tests). Additional testing of T-allele under recessive and codominant models have also confirmed the above findings (p values are >0.05 for all tests). Hence, it is concluded that T allele is not associated with IBD development in Saudi population (p = 0.71, X² = 0.13, OR (95% CI) = 1.09 [0.67–1.73]) (Table 2).

For rs2066844 (c.2023C > T: p.R702W, NOD2) polymorphism, the homozygous C allele is the predominant genotype (99%) in the Saudi Arabia population. Though the homozygous T allele was completely absent, one individual in both case and control groups has possessed this allele in heterozygous form. The extreme rare prevalence of the minor allele (T) in the Saudi population indicates that this polymorphism is not associated with IBD development (p = 1, X² = 0.03, OR (95% CI) = 1.07 [0.06–17.32]) (Table 2).

The major allele is predominant in both case and control group subjects for rs76418789 (c.445G > A: p.G149R, IL23R) polymorphism. The homozygosity of minor allele was found in one IBD patient, but absent in controls. Moreover, the heterozygous AG allele is also missing in both cohorts (p values are >0.05 for all tests). Both recessive and co-dominant models have further confirmed that there is no statistical association of this ‘A’ allele with disease risk in Saudi IBD patients (p = 0.07, X² = 3.23, OR (95% CI) = 0.98 [0.97–1.00]) (Table 2).

For rs11209026 (c.1142G > A: p.R381Q, IL23R) polymorphism, the minor allele is completely absent in the IBD patient group but found in one control subject. The distribution of heterozygous and homozygous (major allele) genotypes was comparable in study cohorts. Both dominant and recessive models have also not shown any significant difference in their distribution (p values are >0.05 for all tests). Hence, it is concluded that the A allele is not associated with the IBD risk in Saudi population (p = 0.58, X² = 0.30, OR (95% CI) 0.72 [0.23–2.26]) (Table 2).

3.4 Epistatic interaction between NOD2 and IL23R SNPs and IBD risk

The MDR test was used to identify SNP-SNP interactions in IBD cases and controls. MDR results suggested that none of the SNP-SNP combinations of NOD2 and IL23R genes is associated with a statistically significant risk of IBD (Table 3).

3.5 Sanger sequencing validation of TaqMan allelic discrimination assay results

The success rate of TaqMan assay genotyping for NOD2 SNPs i.e., rs2066842 (P268S) and rs2066844 (R702W) in IBD cases is 95.87% (93/97) and 100% (100/100) in controls. For IL23R, the success ratio of genotyping the rs76418789 (G149R) SNP in the case group is 91.75% (89/97), and in the control group, it is 98% (98/100). The rs2066844 (R702W) SNP was successfully genotyped in 87.62% (85/97) of cases and 100% (100/100) of control group subjects. The TaqMan genotype calls were further validated by directly sequencing target DNA regions harboring NOD2 and IL23R SNPs. The Sanger sequencing chromatograms of these 4 SNPs are shown in Figs. 1, 2. Our real-time PCR and Sanger sequencing results have shown a 100% match in genotype calling of IBD SNPs.

Close

Fig. 1

NOD2 gene and its exon localization, the Sanger sequencing chromatograms of NOD2 SNPs (rs2066842 and rs2066844).

Export to PPT

Close

Fig. 2

Sanger sequencing chromatograms of IL23R SNPs (rs76418789 and rs11209026) GG (wild-type homozygote); GA (heterozygote), and AA (variant homozygote).

Export to PPT

3.6 Computational pathogenicity predictions of SNPs

Table 4 shows the pathogenicity prediction scores of IL23R and NOD2 SNPs. The variant rs2066842 in the NOD2 gene, was predicted to be benign by SIFT, PolyPhen2, CADD and Condel prediction tools. Whereas the other variant (rs2066844), was tested by all the four tools and predict a deleterious functional effect on the NOD2 protein. All four prediction tools have attributed the deleteriousness to rs76418789 and rs11209026 to the structure and functional features of IL23R protein (Table 4).

3.7 3-D modeling of NOD2 and IL23R

NOD2 and IL23R protein crystal structures are not fully available in the Protein Data Bank (PDB). Therefore, the full-length NOD2 and IL23R protein was constructed depending on the hierarchical protein structure modeling approach by the I-TASSER webserver. The best NOD2 protein model predicted by this server was chosen based on the maximum C-score (0.08), template modeling (TM) score (0.72 ± 0.11), and RMSD (8.8 ± 4.6 Å), which indicates good structural similarities between template and the query proteins. The human NOD2 gene encodes a 1040 residue long protein of a 110 kDa that contains two domains in its N- terminal: caspase recruitment domain (CARD) and NACHT domain. A helical linker known as NOD2 winged-helix domain links the N-terminal domain with a central NOD domain that known as an NLRC4 helical domain HD2. The C-terminal consists of three repeated LRR domains. Fig. 3 shows the 3D molecular structure of NOD2 protein.

Close

Fig. 3

The 3-dimensional Protein Structure of NOD2 protein. In the middle Superpose of Pro268Ser and Arg702Trp in NOD2 protein. The wildtype Proline 268 and mutant Ser 268 amino acid interactions with surrounding amino acids in the NOD2 protein. The wildtype Arg 702 and mutant Trp 702 amino acid H- Bond interaction with surrounding amino acids in the NOD2 protein.

Export to PPT

The best model of IL23R protein has been chosen based on the maximum C-score (-1.39) that demonstrates the predicted model's reliability, Template Modeling (TM) score (0.54 ± 0.15), and RMSD (11.1 ± 4.6 Å), indicate the structural similarities between template and the query proteins. The C-score value ranged between −5 to 2, which reflect the good quality of predicted protein models. The human IL23R gene encodes a 629-residue protein that weighs 71.7 kDa and consists of Interleukin-6 receptor alpha chain binding domain, fibronectin type-III domain, cytokine interleukin-12p40 C-terminus, and Ig-like C2-type domain. Fig. 15-b shows the 3D molecular structure of IL23R protein. NOD2 and IL23R generated protein models were further energy minimized using steepest descent algorithm and models were minimized around 1500 cycles.

3.8 Stability analysis of NOD2 and IL23R missense mutations

Table 5 shows free energy (ΔΔG) values for NOD2 (P268S and R702W) and IL23R (G149R and R381Q) variants by DUET webserver, which calculates the combined and consensus predictions of both mCSM and SDM methods (Pires et al., 2014). As per DUET results, all the 4 variants were predicted to be destabilizing to their corresponding proteins owing to their free energy values. The rs2066842 variant results in the substitution of Proline to Serine amino acid at the 268th residue of NOD2 protein. Pro268 Relevant Solvent Accessibility (RSA) property is 66%, and the Absolute Surface Accessibility (ASA) degree is 93 Å. The Pro268 binds to Ala266 by a total of 2.78 hydrogen bonds. The RSA property percentage of Serine is 56%, the ASA degree is 65 Å, and it binds to Asp262 by 2.62 hydrogen bonds. The total RMSD values of the c-alpha atoms in each amino acid of wild and mutant NOD2 protein models have been calculated to estimate the structural drifts at both amino acid residues and polypeptide chain levels. The RMSD value of the normal NOD2 protein model is 0.0 Å, and the RMSD value of the mutant model is 0.02 Å. The alteration in hydrogen bonds and RMSD values may alter the structure and function of the protein (Fig. 3).

In case of rs2066844 variant, the native Arg is replaced by Trp at 702th residue of NOD2 protein. The RSA property percentage of the Arg702 is 35%, the ASA degree is 81Å, and this amino acid forms a strong H-bond (3.02Å and 3.01Å) with Lys698 and Cys706 , respectively. The RSA property percentage of the Trp702 is 18%. The ASA degree of this amino acid is 42 Å, and it binds to Lys698 and Cys706 by 3.05Å hydrogen bonds. The RMSD value of the NOD2 normal model is 0.0 Å, and the mutant model is 0.00 Å. No significant difference between the RMSD values in both models at the protein level (Fig. 3). The rs76418789 variant substitutes Gly to Arg at position 149 of the IL23R protein. The Gly149 RSA property percentage is 28%, the ASA degree is 22 Å, and this amino acid binds with Asp126 by 2.83Å hydrogen bonds. The RSA property percentage of the Arg149 is 43%, and the ASA degree is 98 Å. The Arg149 binds with Lys126 by 2.88Å hydrogen bonds. The normal model’s RMSD value is 0.0 Å, and the value of the mutant model is 0.01 Å. The alteration in hydrogen bonds and RMSD values may change the structure and function of the protein (Fig. 4).

Close

Fig. 4

The 3D view of IL23R protein structure. In the middle Superpose of Gly149Arg and Arg381Gln in IL23R protein. The Wildtype Gly 149 and mutant Arg 149 forming H-bonds with surrounding amino acids. c) The wildtype Arg 381 and mutant Gln 381 forming H-bonds with surrounding amino acids of IL23R.

Export to PPT

The normal amino acid in the rs11209026 variant is an Arginine at position 381of the IL23R protein, and the mutant amino acid is Glutamine. The RSA property percentage of Arg381 is 60%, and the ASA degree is 138 Å. The Arg381 links with Gln459 in two regions by 2.99 and 2.98 hydrogen bonds. On the other hand, the RSA property percentage of the Gln381 is 63%, the ASA degree is 112 Å, and this amino acid does not link with any other amino acids in the same protein. The RSMD value of the normal IL23R protein model is 0.0 Å, and the value of the mutant model is 0.01 Å. The changes in hydrogen bonds and RMSD values may alter the functional and structural properties of the protein (Fig. 4).