Dissecting the contribution of single nucleotide polymorphisms in CCR9 and CCL25 genomic regions to the celiac disease phenotype.

PURPOSE AND
OBJECTIVES: Given their role in homing immune cells to the intestine, CC motif chemokine receptor 9 (CCR9) and its specific ligand CC motif chemokine ligand 25 (CCL25) are interesting candidate genes for celiac disease. These genes are located in regions previously shown to be associated with or linked to celiac disease, but no investigations on their association with various celiac disease phenotypes have so far been conducted. Here we studied such associations of both genotyped and imputed single nucleotide polymorphisms (SNPs) with either regulatory function or exonic location of the CCR9 and CCL25 loci.
RESULTS: Exploiting a carefully phenotyped cohort of 625 celiac disease patients and 1817 non-celiac controls, we identified that multiple SNPs with predicted regulatory function (RegulomeDB score ≤3a and/or eQTL effect) located between 100 kB upstream and downstream of CCR9 and CCL25 are associated with celiac disease and/or selected phenotypes. Of the genotyped SNPs in the CCR9 loci, rs213360 with an eQTL effect on CCR9 expression in blood was associated with celiac disease and all investigated phenotypes except high HLA risk. Rs1545985 with an eQTL on CCR9 expression and rs7652331 and rs12493471, both with RegulomeDB score ≤3a, were all associated with gastrointestinal symptoms and malabsorption and the latter additionally with anemia. The genotyped CCL25 SNPs rs952444 and rs882951, with RegulomeDB scores 1d and 1f respectively and eQTL effect on CCL25 expression in small intestine, were associated with gastrointestinal symptoms and malabsorption. The CCL25 SNP rs2303165 identified in sequencing followed by imputation was associated with partial villous atrophy. However, the association did not pass the permutation based multiple testing correction (PEMP2 > 0.05).
CONCLUSIONS: We conclude that SNPs in the region of CCR9 and CCL25 with predicted functional effect or exonic localization likely contribute only modestly to various celiac disease phenotypes.

Introduction

Celiac disease is an immune-mediated chronic condition where oral tolerance to dietary gluten has been lost. The multifaceted disease can present with varying gastrointestinal and/or extraintestinal signs and symptoms. In addition, patients may be completely asymptomatic [1]. A prerequisite for the development of the disease is the presence of human leukocyte antigen (HLA) molecules HLA-DQ2 or -DQ8 and individuals homozygous for HLA-DQ2 encoding alleles are at a particularly high risk [2,3]. Celiac disease, with its autoimmune nature, is characterized by an IgA class autoantibody response against endogenous enzyme transglutaminase 2 (TG2) measurable in serum. In patients with celiac disease, the ingestion of gluten leads to small-bowel mucosal villous atrophy and crypt hyperplasia of varying severity. Moreover, a profound mucosal inflammation, characterized by increased density of intraepithelial lymphocytes and infiltration of both T cells and plasma cells in the lamina propria is usually present [1]. The inflammatory cells arrive at the small bowel from secondary lymphoid tissue in a process involving the gut-homing chemokine CC motif chemokine ligand 25 (CCL25) and its specific receptor CC motif chemokine receptor 9 (CCR9) [4].

In addition to being involved in homing immune cells to the intestine, the main organ affected by celiac disease, the CCR9-CCL25 axis is implicated in celiac disease by several other studies. Firstly, CCR9-positive dendritic cells play an essential role in maintaining gut homeostasis and tolerance by regulating the phenotype and function of both innate and adaptive immune cells [5]. Secondly, there is data to suggest that the number of CCR9-positive type 2 conventional dendritic cells is increased in circulation in untreated celiac disease [6] while the number of CCR9-expressing T cells is diminished in the small intestines of patients [7]. Thirdly, a clinical trial to test the efficacy of an oral CCR9 inhibitor as an alternative form of treatment for celiac disease has been performed but not yet published [8]. Lastly, the genes coding for CCL25 and CCR9 are located in chromosomal regions previously linked to or shown to be associated with celiac disease [[9], [10], [11]]. However, in earlier genetic association studies the coverage of single nucleotide polymorphism (SNPs) at CCR9 and CCL25 loci, particularly those in exonic regions, has been sparse. Moreover, to the best of our knowledge no studies focusing on genomic variants at these loci on different celiac disease phenotypes exist.

For the above mentioned reasons, we considered CCR9 and CCL25 to be potential positional and functional candidate genes in celiac disease. Consequently, we identified SNPs in CCR9 and CCL25 genomic regions that likely have functional effects or are located in the exons and tested their association with celiac disease and selected disease phenotypes.

Materials and methods

Patients and controls

The study was conducted at Tampere University and Tampere University Hospital. Altogether 1048 celiac disease patients with biopsy-proven small bowel mucosal damage were recruited with the assistance of national and local celiac societies and by nationwide media announcements. The patient information at diagnosis was collected from medical records and from supplementary interviews by a physician or a study nurse with expertise in celiac disease. In the case of children, the guardian was interviewed. The structured interviews included questions on celiac disease diagnosis, symptoms at the time of diagnosis and in childhood, and associated medical conditions. Whole blood samples were drawn for genetic analysis. In order to avoid bias (inflation of type 1 error) caused by the inclusion of several subjects from the same family, the present study considered only one randomly selected celiac case with full genotype available per family, resulting in 625 cases. The median age of the patients was 41 (range 0.5–79) years. Patients’ demographic data as well as their clinical characteristics in terms of selected celiac disease-associated phenotypes (gastrointestinal symptoms, malabsorption, anemia, severity of small bowel mucosal damage, HLA risk categories, and celiac disease antibodies) are presented in Table 1. DNA samples from 144 HLA-DQ2-positive biopsy-proven celiac disease patients and 144 non-celiac controls outside the genotyped celiac disease patient cohort were subjected to sequencing. The study design, patient recruitment, and collection of patient record data were approved by the Regional Ethics Committee of Tampere University Hospital. All participants, or in the case of children their legal guardians, gave written informed consent.

Table 1

Demographic data and selected celiac disease phenotypes of 625 celiac disease patients at diagnosis.

	N	%
Females	489	78
Gastrointestinal symptomsa	526	84
Malabsorptionb	267	43
Anaemia	157	25
Small bowel mucosal damagec
Total or subtotal villous atrophy	361	66
Partial villous atropthy	185	34
HLA riskd
High	98	16
Intermediate/low	527	84
Celiac disease autoantibodiese
Positive	279	94
Negative	19	6

Diarrhea, abdominal pain, flatulence, heartburn, nausea, vomiting.

Anemia, vitamin and micronutrient deficiencies.

Small bowel mucosal morphology data was available from 546 patients.

High risk (DQ2.5/DQ2.5; DQ2.5/DQ2.2), intermediate risk (DQ2.5/X, DQ2.2/D2.2, DQ2.2/X, DQ8/DQ8, DQ8/X), low risk (DQ7/X, DQ7/DQ7).

Autoantibody data (endomysial antibodies and/or tissue transglutaminase antibodies) was available from 298 patients.

As non-celiac controls, altogether 1817 subjects with information on gender (1032 males, 785 females) and HLA-genotype (Low HLA risk, N = 939; Intermediate HLA risk, N = 855; High HLA risk, N = 23) from the population representative cohorts FINRISK and Health 2000 [12] were included in the study. Ethics committee approvals were available from the National Public Health Institute's Ethics committee and the Ethics committee in Epidemiology and Public Health in the Hospital District of Helsinki and Uusimaa.

Association analysis of the genotyped CCR9 and CCL25 SNPs

Genotypes for all subjects had been produced by Illumina 610-Quad BeadChip array (Illumina Inc., San Diego, CA, USA) [9]. First, variants spanning between 100 kB upstream and downstream of CCR9 and CCL25 genes were identified resulting in 105 SNPs. Thereafter, FUMA (Functional Mapping and Annotation of GWAS) platform [13] was used to annotate SNPs fulfilling our functional annotation criteria. Publicly available genome-wide association study (GWAS) summary statistic results in which our cohort has been included [9] and a pre-defined list with the 105 SNPs of interest were uploaded. The RegulomeDB 2.0 was used to identify all the SNPs with known and predicted regulatory elements and to assign them a score ranging from 1a to 7 in RegulomeDB ranking [14]. SNPs with RegulomeDB score between 1a and 3a likely to affect the gene expression, were selected for our study. In addition, the expression quantitative trait loci (eQTL) mapping using public eQTL data was used to select SNPs having significant eQTL effects (FDR <0.05) on the expression of CCR9 gene in blood (GTEx whole blood [15], Blood eQTL [16], BIOS QTL [17], and eQTLGen) and on the expression of CCL25 gene in the intestine (GTEx data for small intestine terminal ileum, colon sigmoid and colon transverse [15]; and of the “CEDAR” study [18], terminal ileum). Altogether 41 SNPs with RegulomeDB ranking from 1a to 3a, and/or tissue-specific eQTL effects were identified. After applying quality control (QC) filtering for missing genotype rate <5% and missing genotype rate differences between the cases and controls (<3%), and minor allele frequency (MAF >5%), 348 cases and all the 1817 controls remained in the analysis, and 39 SNPs passed QC. All markers were in Hardy-Weinberg equilibrium (HWE) (P > 1 × 10−6) in the controls. Allelic associations of the 39 genotyped SNPs with celiac disease and the selected disease phenotypes were tested as described in Section 2.5.

Sequencing exonic regions of CCR9 and CCL25

DNA was extracted from whole blood samples or from leukocyte enriched buffy coats. Extractions were performed using FlexiGene DNA kit (Qiagen, Hilden, Germany). All coding exons and exon-intron boundaries of both CCR9 and CCL25 were amplified with PCR using primer pairs represented in Supplementary Table 1. The PCR reactions (20 μL) contained 10 mM Tris-HCl (pH 8.8 at 25 °C), 1.5 mM MgCl2, 50 mM KCl, 0.1% Triton X-100 in 10X Optimized DyNAzyme Buffer (Thermo Fisher Scientific, Waltham, MA, USA), 0.5 mM dNTPs, 1.0 μM of each primer, 80 ng genomic DNA and 1.0 U DyNAzyme II DNA Polymerase (Thermo Fisher Scientific). PCR conditions included initial denaturation at 95 °C for 5 min, 40 cycles of 95 °C for 30 s, 55–60 °C (depending on the PCR product) for 40 s and 72 °C for 1 min. The final extension was performed at 72 °C for 5 min. Successful PCR amplification was confirmed by analyzing PCR products with agarose gel electrophoresis and by UV–Vis spectrophotometry (NanoDrop, Thermo Fisher Scientific). Sequencing was performed by Macrogen Europe B.V. (Amsterdam, the Netherlands) exploiting ABI 3730 (Applied Biosystems, Thermo Fisher Scientific). The sequences were analyzed for SNPs using Sequencher 4.10.1 software (Gene Codes Corporation, Ann Arbor, MI, USA). Sequence data was analyzed further by calculating allele frequencies for found SNPs in the study material.

Imputation of the identified exonic SNPs and their association analysis with celiac disease and its phenotypes

As SNPs identified by sequencing were not available in the Illumina 610-Quad BeadChip array, and thus not genotyped, those having MAF ≥5% were selected to be phased and imputed using a Finnish population-specific panel of 3775 high-coverage (25–30 ×) whole-genome sequences (SISu v3). Sample-wise, variant-wise, and post-imputation QC was applied as previously described [19]. Phasing of genotyped data was performed with Eagle 2.3.5 (https://data.broadinstitute.org/alkesgroup/Eagle/) and imputation was carried out with Beagle 4.1 (version 08Jun17.d8b, https://faculty.washington.edu/browning/beagle/b4_1.html) as described in the following protocol: dx.doi.org/10.17504/protocols.io.nmndc5e. In the post-imputation QC, SNPs with good imputation quality metrics (INFO score ≥0.8) were included. The same QC filtering criteria applied in the association analysis of genotyped SNPs (section 2.2) were applied to the imputed genetic data. Out of the sequenced exonic variants, nine SNPs were found in the Finnish reference panel, resisted QC filtering, and were thus included in the association analyses.

Statistical analyses

Association analyses of genotyped SNPs were performed using PLINK v1.07 (https://zzz.bwh.harvard.edu/plink/). Results are presented as odds ratios (OR) with 95% confidence intervals (95% CI). In these analyses, associations were adjusted for multiple testing and small sample size groups by using 104 permutation analysis. The generated empirical PEMP2 value ≤ 0.05 (uncorrected P value ≤ 0.001) was assumed to be statistically significant. In the imputation analysis, the genotypes’ probability dosages were handled by BCFtools (https://samtools.github.io/bcftools/) and the association analysis performed using PLINK 2.0 (https://www.cog-genomics.org/plink/2.0/). The associations of the post-imputation genotype probabilities of the nine SNPs with celiac disease and selected celiac disease phenotypes were tested using the frequentist likelihood score method implemented in SNPTEST v2.5.2. Associations reaching our permutation threshold described above were considered to be statistically significant [19,20].

Results

Association of genotyped functionally annotated CCR9 and CCL25 SNPs with celiac disease and its distinct phenotypes

Of the genotyped CCR9 SNPs, rs2133660 with an eQTL effect on the expression of CCR9 in whole blood and peripheral blood mononuclear cells (PBMCs) was associated with the presence of malabsorption (OR = 1.45, 95% CI = 1.14–1.83, P = 0.002) although this did not resist correction for multiple testing (PEMP2 = 0.064) (Fig. 1, Supplementary Table 2). The same SNP had nominal associations (PEMP2 > 0.05) with celiac disease (OR = 1.20, 95% CI = 1.02–1.42, P = 0.031), with the presence of gastrointestinal symptoms (OR = 1.29, 95% CI = 1.08–1.54, P = 0.005), anemia (OR = 1.35, 95% CI = 1.01–1.81, P = 0.041), partial villous atrophy (OR = 1.35, 95% CI = 1.01–1.79, P = 0.039), and with negative celiac disease autoantibodies (OR = 3.48, 95% CI = 1.38–8.74, P = 0.005) (Fig. 1). Further, rs12493471 and rs7652331 with RegulomeDB scores ≤ 3a and rs1545985 with an eQTL effect on CCR9 expression were all associated with the presence of gastrointestinal symptoms and malabsorption and rs12493471 also with anemia. Four SNPs (rs12983784, rs952444, rs882951 and rs11667975) in the CCL25 gene region were nominally associated (PEMP2 > 0.05) with more than one phenotype, and all had a RegulomeDB score ≤ 3a (Fig. 2). Of these, rs882951 and rs952444, with the lowest RegulomeDB scores (1d and 1f respectively), had significant eQTL effects on CCL25 expression in small intestine and were both associated with the presence of gastrointestinal symptoms and malabsorption (for both SNPs OR = 1.20, 95% CI = 1.01–1.43, P = 0.034 and OR = 1.32, 95% CI = 1.04–1.66, P = 0.021 respectively). Rs12983784 was associated with celiac disease (OR = 1.20, 95% CI = 1.01–1.44, P = 0.041) and positive serum autoantibodies (OR = 1.43, 95% CI = 1.12–1.82, P = 0.004). Rs1129763 and rs11667975 were further associated with the presence of total/subtotal villous atrophy (OR = 0.65, 95% CI = 0.48–0.89, P = 0.006 and OR = 1.31, 95% CI = 1.03–1.68, P = 0.027, respectively), rs11667975 with high HLA risk (OR = 4.20, 95% CI = 1.21–14.53, P = 0.015), and positive serum autoantibodies (OR = 1.36, 95% CI = 1.04–1.79, P = 0.025). The detailed results are presented in Fig. 1, Fig. 2 and Supplementary Table 2.

Fig. 1

Forest plot representing association of SNPs in the genomic region of the candidate gene CCR9 with celiac disease and its different phenotypes. eQTL; expression quantitative trait loci, OR; odds ratio, CI; confidence interval, PEMP2; empirical P value at 10,000 permutation threshold, PBMC; peripheral blood mononuclear cell, PVA; partial villous atrophy, CeD; celiac disease.

Fig. 2

Forest plot of representing association of SNPs in the genomic region of the candidate gene CCL25 with celiac disease and its different phenotypes. SNPs identified in sequencing analysis and imputed are indicated by an asterisk. eQTL; expression quantitative trait loci, OR; odds ratio, CI; confidence interval, PEMP2; empirical P value at 10,000 permutation threshold, PBMC; peripheral blood mononuclear cell, TVA/SVA; total or subtotal villous atrophy, PVA; partial villous atrophy, CeD; celiac disease.

Association analysis using sequenced and imputed data

By sequencing the exons of CCR9 and CCL25, we identified two variants in CCR9 and eight in CCL25 with MAF ≥5% (Supplementary Table 3). One of the variants found in CCR9 was previously unreported (referred to as CCR9ex3snp7 in Supplementary Table 3) and changed threonine to alanine at position 100. However, this variant could not be imputed and therefore association analyses were carried out with the remaining nine SNPs. None of the SNPs were associated with celiac disease, but of the CCL25 SNPs, the synonymous variant rs2303165 was nominally associated (PEMP2 > 0.05) with partial villous atrophy (OR = 0.39, 95% CI = 0.15–1.00, P = 0.038) (Fig. 2), but the SNP had no eQTL effect on CCL25 expression in the intestine.

Discussion

Our study identified several nominal associations for SNPs in the genomic regions of CCR9 and CCL25 genes with celiac disease or distinct disease phenotypes. As regards CCR9, associations were detected with several genotyped SNPs with potential or proven functional effect but not with any of the exonic imputed ones. None of associated CCR9 SNPs had both RegulomeDB score less than 3a and an eQTL effect on CCR9 expression. Of the SNPs with RegulomeDB score less than 3a, rs12493471, associated with gastrointestinal symptoms, malabsorption, and anemia in the present study has previously been associated with celiac disease [21]. SNPs with RegulomeDB scores ranging from 1a to 1f are by definition predicted to affect different DNA regulatory elements and have an eQTL effect based on the Encyclopedia of DNA Elements (ENCODE) project, Gene Expression Omnibus data, and published literature [14,22] and are thus functionally particularly interesting. Rs7652331 at the CCR9 locus associated with gastrointestinal symptoms and malabsorption falls into this category with RegulomeDB score 1f. However, as the SNP did not have an eQTL on CCR9 expression, it is likely that the effect of the SNP on the phenotype is due to other genes than CCR9. Of the SNPs with an eQTL on CCR9, rs2133660 mapping to intronic region of FYCO1 gene downstream of CCR9 was associated with malabsorption along with an association with celiac disease and several other phenotypes. However, as the RegulomeDB score of rs2133660 was over 3a, it is possible that the phenotypes associated with this variant rather reflect the effect of a proxy SNP causing increased CCR9 expression in lymphocytes. Moreover, rs2133660 also had eQTL effects on other nearby genes, many of them belonging to different chemokine receptor families. Thus, the effect of the SNP on different phenotypes is hardly attributable solely to increased CCR9 expression.

As regards the SNPs in the genomic region of CCL25, genotyped SNPs rs952444 and rs882951 both showed nominal associations with the presence of gastrointestinal symptoms and malabsorption. These variants had the highest likelihood (RegulomeDB score ≤ 1f) of being connected to functional transcriptional effects in gene expression and they both had eQTL effects on the expression of CCL25 in the small intestine [15,18]. Although as far as we know no studies have addressed CCL25 expression in the small intestine in different celiac disease phenotypes, the lack of evidence for altered CCL25 expression in celiac disease [[23], [24], [25]] would suggest that the SNPs exert their effects through other mechanisms than CCL25 expression. Sequencing CCL25 exons followed by imputation revealed that rs2303165 located in CCL25 gene region was associated with partial villous atrophy, without eQTL CCL25 effects in intestine. This exonic CCL25 SNPs could exert its effects on the phenotypes by yet to be determined mechanisms that do not directly involve CCL25.

The associations identified in our study did not pass the 10,000 permutation based multiple testing correction. Although there is a debate regarding the number of permutations required, 10,000 permutations provide empirically adjusted P-values with strong evidence of association [26,27]. As the permutation threshold used in the present study is quite stringent and because of our small sample size, we had limited statistical power to detect strong associations. Thus, our findings need to be confirmed in further studies with larger cohorts of carefully phenotyped celiac disease patients. In any case, CCR9 and CCL25 variants likely make only a minor contribution to the generation of celiac disease phenotypes addressed in the present study. We have previously reported associations of SNPs outside the CCR9 and CCL25 genomic regions with the same phenotypes that were resistant to correction for multiple comparisons. One of these was the association of rs13010713 in integrin subunit alpha 4 gene (ITGA4) with the presence of gastrointestinal symptoms, total or subtotal villous atrophy, and intermediate HLA risk [19]. ITGA4 codes for the α4 subunit of heterodimeric integrin molecules involved in adhesion and the pairing of α4 subunit with β7 subunit promotes homing of T cells to intestinal sites [28]. Interestingly, expression of CCR9 and α4β7 on T cells and their subsequent localization to the gut is required for induction of oral tolerance, at least in mice [29]. Thus, due to this solid firm connection between CCR9 and α4β7, variants found in the integrin locus may possess changes of functionality in the gut-homing pathway which are reflected in celiac disease phenotypes.

Undoubtedly, a given celiac disease phenotype likely also has nongenetic determinants. These may include environmental factors such as microbes, infections or the amount of gluten consumed by a patient, all of which have been associated with the development of celiac disease [30,31]. Moreover, delay in the diagnosis of celiac disease may allow the disease to progress to a more severe form, thus affecting some of the phenotypes at the time of diagnosis. Data showing that at least diarrhea, abdominal pain, and malabsorption are associated with long diagnostic delay lends credence to this hypothesis [32].

The main strength of the present study is the carefully phenotyped cohort of celiac disease patients. In addition, the exploitation of the imputed data allowed us to study the exonic SNPs that were not present in the Illumina 610-Quad BeadChip array previously used to study celiac disease associations. As a limitation, the sequencing approach was performed using Sanger sequencing, which has its shortcomings in accuracy compared to Next-Generation Sequencing (NGS) techniques with more comprehensive coverage and higher sensitivity to detect low-frequency variants. To overcome this, we used MAF ≥5% as our cut-off. However, this may have resulted in missing interesting variants, including rs12721497 in CCR9 with MAF = 0.01, according to the 1000 Genomes Project, associated with acute and chronic stage graft versus host disease [33]. In addition, unfortunately eQTL data, blood or small bowel mucosal samples from celiac disease patients were not available which precluded us from addressing eQTL effects particularly in patients or studying the effects of the SNPs on immunological changes in celiac disease more generally. Moreover, the number of individuals in our study cohort was rather small, particularly for genotype-phenotype association analyses.

Conclusions

We conclude that SNPS in the region of CCR9 and CCL25 having probable functional effect or being located in exons are weakly associated with various celiac disease phenotypes. Our results thus suggest that, regardless of the importance of CCR9 and CCL25 in maintaining gut homeostasis and tolerance, variation within these genomic regions likely makes a minor contribution to the phenotype of celiac disease.

Funding

This study was supported by the , the , the , the , the , and the Competitive State Research Financing of the Expert Area of . Funding sources were not involved in the design of the study, in the collection, analysis or interpretation of data, or in the decision to submit the article for publication.

Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.