Association analysis of non-synonymous polymorphisms of interleukin-4 receptor-α and interleukin-13 genes in canine atopic dermatitis.

Interleukin-4 (IL4) and interleukin-13 (IL13) are involved in the initial response of T helper 2 lymphocytes through the activation of the IL4 receptor alpha (IL4RA), which is a common receptor chain for these cytokines. In humans, several single-nucleotide polymorphisms (SNPs) identified in the IL4R and in interleukin coding genes were associated with atopic disorders. However, the association between canine IL4R polymorphisms and atopic disorders has not been investigated yet. This study aimed to determine the associations between four non-synonymous SNPs and canine atopic dermatitis (CAD) in shiba inu and miniature dachshund populations. Polymerase chain reaction (PCR) and restriction fragment length polymorphism (RFLP) analysis were used to genotype four polymorphisms of canine IL4R and IL13 in 34 shiba inu and 19 miniature dachshund patients with CAD, as well as 29 shiba inu and 39 miniature dachshund patients without the condition. Results from miniature dachshunds revealed a potential association between the presence of minor A allele rs24378020 and CAD (odds ratio, 0.10; 95% confidence interval, 0.01-0.85; Poriginal=0.0062). This CAD resistance allele led to an amino acid substitution (Arg688Cys) that could impair IL4 and IL13 signaling. In shiba inu patients, rs24378020 was fixed by homozygosity of the major G allele. No association was found between the remaining three evaluated SNPs and CAD. Nevertheless, the study suggests that the IL4R Cys688 variant reduces the risk of CAD in miniature dachshunds.

Canine atopic dermatitis (CAD) is a common allergic inflammatory skin disease [23] that is very similar to human atopic dermatitis (AD). AD is characterized by skin barrier dysfunction that is triggered by environmental factors and altered immune system responses, resulting in eczematous and itchy lesions [21, 25]. Humans with allergies exhibit higher total serum immunoglobulin E (IgE) levels as consequence of T helper 2 (Th2) responses, which are promoted by cytokines, such as interleukin-4 (IL4) and interleukin-13 (IL13) [7, 39]. Dogs with CAD demonstrate, on average, significantly higher levels of allergen specific IgE compared with dogs without CAD [11]. The IL4 receptor (IL4R) gene encodes the alpha chain (IL4RA) of the IL4 and IL13 receptors, which is a key component for promoting the Th2 lymphocyte phenotype and IgE production [19]. Therefore, polymorphisms in the IL4R, as well as in the IL13 and IL4 coding genes, may contribute to the complex regulation of atopy phenotypes. In human IL4R, the coding single-nucleotide polymorphisms (SNPs) rs1801275 and rs1805010, which cause p. Gln576Arg and p. Ile75Val amino acid substitutions, are strongly associated with atopy [8, 14, 22, 27]. Additionally, the human IL13 coding SNP rs20541, which causes p. Arg130Gln substitution, is associated with susceptibility to allergies [4, 12, 34]. In the human IL4 gene, a SNP related to atopy was reported in the promoter region −590C/T, but coding SNPs associated with atopy have not been described [17, 18].

According to Ensembl genome browser records, there are 664, 29, and 49 genetic variants in the canine IL4R, IL4, and IL13, respectively [6, 16]. However, their associations with CAD have not been investigated. The majority of these variants are located on the upstream and downstream regions surrounding the gene locus and introns where it is difficult to estimate the impact on gene function [16]. Conversely, bioinformatics tools can estimate the impact of non-synonymous coding SNPs on gene function [1, 5, 29]. Therefore, we focused on non-synonymous coding SNPs that result in amino acid variations in the canine IL4R and IL13. We evaluated potential associations between four non-synonymous SNPs and CAD in shiba inu and miniature dachshund populations in Japan, using CAD cases and non-CAD controls that were diagnosed at veterinary clinics. This study provides new information on the contribution of polymorphisms for achieving resistance to CAD.

MATERIALS AND METHODS

Sample collection

Dog blood samples were donated by attending veterinarians to the canine bio-resource banking project (Azabu University) with the consent of the dog owners. Blood collection for research purposes was approved by the animal research committee of Azabu University (Permission number: 130607).

Cases of CAD included dogs presenting with non-seasonal chronic pruritus, according to Favrot’s criteria [9]. Diagnosis was made by a veterinarian with clinical experience in dermatology. Causes of non-atopic dermatitis, such as bacterial and fungal skin infections and ectoparasites, were ruled out by routine examinations. Because CAD usually affects young adult dogs (onset before three years-old) [13], unaffected dogs over the age seven (double of susceptible age) were selected as control dogs. Control dogs had variable medical histories and were randomly selected. The cases of CAD consisted of 34 shiba inu and 19 miniature dachshund, and the control group comprised samples from 29 shiba inu and 39 miniature dachshund. Shiba inu case group included 10 males and 21 females, and three individuals without sex records, and their average age at the time of blood sampling was 7.7 ± 4.0 years. Shiba inu control group comprised 11 males and 10 females, and eight individuals without sex records, and their average age was 11.6 ± 3.1 years. Miniature dachshund case group included seven males and 11 females, and an individual without sex record, and their average age was 6.8 ± 3.4 years old. Miniature dachshund control group was comprised 26 males and 22 females, and an individual without sex record, and their average age was 11.4 ± 2.7 years.

Bioinformatics analysis

In the canine IL4R, seven non-synonymous substitutions were found in the Ensembl genome browser [6, 16]. Of these, rs24378020 and rs851400460 were present in both the Broad Institute’s Dog SNP library [20] and the European Variation Archive [32], and genotype was confirmed in 85% or more among 238 individuals analyzed on the high quality variant calls from multiple dog genome project−Run1 [32]. Therefore, we chose these SNPs as top priority for analysis. Additionally, rs9193906 on the C-terminus of IL4RA was included for analysis. In the canine IL13, we targeted the only non-synonymous substitution ever reported, rs22147008. The canine IL4 was excluded from the study as no non-synonymous substitutions have been reported.

Lastly, we selected four non-synonymous SNPs known to cause amino acid substitutions in mature IL4RA and IL13 molecules: rs9193906 G >C, rs24378020 G >A, and rs851400460 G >A for IL4R, and rs22147008 A >G for IL13 (Table 1). The impact of the analyzed SNPs on IL4RA and IL13 were estimated through PolyPhen-2 (http://genetics.bwh.harvard.edu/pph2/) [1], SIFT (https://sift.bii.a-star.edu.sg) [29], and PROVEAN (http://provean.jcvi.org/index.php) [5] bioinformatics tools. Combining prediction results from multiple tools can increase the chance of identifying functional variants that had been missed by other tools [5]. PolyPhen-2 predicts the effect of an amino acid substitution on the structure and function of a protein based on a number of features that characterize the substitution, such as sequence, phylogeny, and structural information [1]. The PolyPhen score ranges from 0.0 to 1.0 and represents the probability that a substitution is damaging, with values closer to 1.0 being more confidently predictive of a deleterious effect. The PolyPhen scores are categorized according to “probably damaging”, “possibly damaging”, and “benign”. SIFT predicts whether an amino acid substitution affects the protein function based on sequence homology and the physicochemical properties of the alternate amino acid. The SIFT score shows normalized probabilities for all possible substitutions from the alignment, with scores lower than 0.05 being “deleterious” and scores greater than or equal to 0.05 being “tolerable” [29]. PROVEAN predicts the functional effect of protein sequence variations based on the similarity of a semi-global pairwise sequence alignment score between the query sequence and each of the related sequences. PROVEAN can predict a “deleterious” or “neutral” effect based on scores equal or lower than, or greater than a pre-defined threshold, respectively. In this study, the threshold value used was −2.5 (default setting).

Table 1.

Details of the non-synonymous single nucleotide polymorphisms (SNPs) analyzed of IL4R and IL13, and in silico prediction of their possible effects

Gene	Dog chromosome	Reference SNP ID	SNP location	DNA strand	Amino acid substitution	In silico prediction a)
						PolyPhen	SIFT	PROVEAN
IL4R	6	rs9193906	NC_006588.3	Template	p.Thr220Arg	0.45	0.49	−0.673
			g.19263381 G>C			Benign	Tolerated	Neutral
		rs24378020	NC_006588.3	Template	p.Arg688Cys	0.76	0.06	−0.709
			g.19255617 G>A			Possibly damaging	Tolerated	Neutral
		rs851400460	NC_006588.3	Template	p.Ser781Leu	0.00	1.00	0.335
			g.19255337 G>A			Benign	Tolerated	Neutral
IL13	11	rs22147008	NC_006593.3	Coding	p.Thr81Ala	0.99	0.04	−4.86
			g.20960082 A>G			Probably damaging	Deleterious	Deleterious

a) Bioinformatics tools used to predict the impact of an amino acid substitution on protein function.

DNA extraction and genotyping

We performed polymerase chain reaction (PCR) and restriction fragment length polymorphism (RFLP) analysis to genotype the polymorphisms of the canine IL4R and IL13 (Table 1). Genomic DNA was extracted from peripheral blood using the QuickGene DNA whole blood kit L (KURABO, Osaka, Japan). The extraction of genomic DNA was confirmed using a NanoDrop ND-1000 Spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). PCR amplifications were performed in 20 µl reaction volumes according to the instructions of the KAPA Taq DNA Polymerase kit (Kapa Biosystems, Wilmington, MA, USA). The primers and protocol conditions used for the amplification of IL4R and IL13 polymorphisms are listed in Table 2. The effectiveness of the primers used was confirmed by Sanger sequencing of the amplicons using BigDye Terminator v3.1 (Thermo Fisher Scientific) and ABI 3130 genetic analyzers (Applied Biosystems, Foster City, CA, USA). For RFLP analysis, each PCR product was incubated with two units of specific restriction endonucleases (Table 2) under the conditions of the manufacturer’s recommendations (New England Biolabs, Ipswich, MA, USA). Digested amplicons were separated through electrophoresis on 2.5% agarose gel. The genotypes were defined according to generated fragment patterns, as summarized in Table 2.

Table 2.

Summary of the primers, polymerase chain reaction (PCR) conditions, and genotyping methods used for the four single-nucleotide polymorphisms

Gene	Polymorphism	Primer pairs	PCR conditions	Restriction enzyme	RFLP a) Genotype
IL4R	rs9193906; G >C	5′-GGCTCAGACCTACAACAGCA-3′	94°C, 30 sec; 57°C, 15 sec; 72°C, 25 sec for 30 cycles	PmlI	C: 259 bp (uncut); G: 210 & 49 bp
		5′-GCTCATTTGTTACCGCAGCC-3′
	rs24378020; G >A	5′-CCCCTGTTCACCTTTGGACT-3′	94°C, 30 sec, 60°C, 5 sec; 72°C, 10 sec for 30 cycles	PvuII	G: 309 bp (uncut); A: 191 & 118 bp
		5′-TGAGGACCTGTCTCCACAGC-3′
	rs851400460; G >A	5′-GTGGCAAGGCCCACATAGT-3′	94°C, 30 sec; 60°C, 10 sec; 72°C, 20 sec for 30 cycles	AvaI	A: 311 bp (uncut); G: 160 & 151 bp
		5′-CCTGCAGACATCAGCAACCA-3′
IL13	rs22147008; A >G	5′-TCTCAAACCCCACCTCCTGT-3′	94°C, 30 sec; 65°C, 15 sec; 72°C, 25 sec for 30 cycles	PspOMI	A: 230 bp (uncut); G:168 & 62 bp
		5′-AGGACAGAGGGCCTTACCC-3′

a) RFLP: restriction fragment length polymorphism.

Statistical analysis

Allele frequency, Hardy-Weinberg equilibrium, linkage disequilibrium, and genotypic associations for all the SNPs were statistically assessed using the SNPStats (https://www.snpstats.net) [30]. Association analysis for atopic dermatitis was only performed when the minor allele frequency was more than 5%. Original P-values of less than 0.0084 were considered to be statistically significant. We selected this particular P-value because we tested six hypotheses, and the adjusted α=0.05 became 0.05/6 based on Bonferroni correction.

RESULTS

In silico prediction of the effect of SNP-related amino acid substitution

All three bioinformatics tools used predicted the rs22147008, which results in the amino acid substitution p. Thr81Ala in the 63rd residue in mature IL13, to be harmful to the structure of the protein (Table 1). This result strongly suggests that Ala81 causes loss- or gain-of-function on canine IL13. Moreover, the rs24378020 (p. Arg688Cys in IL4R) was predicted by PolyPhen to be possibly damaging; however, the SIFT and PROVEAN scores did not reach statistical significance. Thus, this substitution potentially affects IL4R function. The other two SNPs on IL4R−rs9193906 (p. Thr220Arg) and rs851400460 (p. Ser781Leu)−were predicted by all three bioinformatics tools to be benign to the protein structure and function.

Genotyping and association analysis

Figure 1 shows examples of the result of PCR-RFLP genotyping for the four SNPs evaluated on canine IL4R and IL13 genes. The genotypic and allelic frequencies in each population are shown in Table 3. No polymorphisms were detected in shiba inu and in miniature dachshund samples related to rs24378020 and rs851400460, respectively, as all individuals were homozygous for the major allele (Table 3). Moreover, for rs851400460 the minor allele frequency was evidently below 5% in shiba inu samples (Table 3). Thus, we excluded rs24378020 from further analysis in shiba inu samples, and rs851400460 in both breeds. The remaining three SNPs exhibited polymorphisms with a minor allele frequency above 5%, and their genotypic distributions did not deviate from Hardy-Weinberg equilibrium, except for rs9193906 in shiba inu samples. No significant linkage disequilibrium (P<0.05) was observed for all combinations of the three SNPs in both breeds.

Fig. 1.

PCR-RFLP results of non-synonymous coding SNPs for canine IL4R and IL13 genes. Digested PCR products were separated by electrophoresis on a 2.5% agarose gel and stained with 0.5 µg/ml ethidium bromide. M indicates the position of the DNA Ladder molecular weight marker. The genotype is indicated above each lane. Panels A to D show the results for rs9193906 G>C, rs24378020 G>A, rs851400460 G>A, and rs22147008 A>G, respectively. The bands marked with * in panels A and C are primer dimers that do not affect genotyping.

Table 3.

Association of genotypes and allelic frequencies in the canine IL4R and IL13 single nucleotide polymorphisms (SNPs) with atopic dermatitis in shiba inu and miniature dachshund

Gene	SNP ID	Breed	Group	Genotype (n)			MAF a)	HWE b)	Allele 2 vs. 1	Dominant model genotype12+22 vs. 11
				11	12	22		P value c)	P value ORd) (95% CIe))	P value OR (95% CI)
IL4R	rs919390			GG	GC	CC
	(G >C)	Shiba inu	Control	9	19	1	0.362	0.05	0.45	0.40
			Case	14	20	0	0.294	0.03
		Miniature dachshund	Control	15	26	8	0.429	0.77	0.171	0.37
			Case	8	11	0	0.289	0.25
	rs24378020			GG	GA	AA
	(G >A)	Shiba inu	Control	29	0	0	0.000	N/A f)	N/A	N/A
			Case	34	0	0	0.000	N/A
		Miniature dachshund	Control	32	16	1	0.184	1.00	0.024	0.0062
			Case	18	1	0	0.026	1.00	OR: 0.12 (0.02–0.74)	OR: 0.10 (0.01–0.85)
	rs851400460			GG	GA	AA
	(G >A)	Shiba inu	Control	29	0	0	0.000	N/A	N/A	N/A
			Case	33	1	0	0.015	1.00
		Miniature dachshund	Control	49	0	0	0.000	N/A	N/A	N/A
			Case	19	0	0	0.000	N/A
IL13	rs22147008			AA	GA	GG
	(A >G)	Shiba inu	Control	26	3	0	0.052	1.00	0.505	0.60
			Case	29	4	1	0.088	0.21
		Miniature dachshund	Control	12	24	13	0.510	1.00	0.445	0.22
			Case	5	6	8	0.579	0.16

a) MAF: minor allele frequency, b) HWE: Hardy-Weinberg Equilibrium, c) P=probability level, d) OR: Odds ratio, e) CI: confidence interval, f) N/A: not applicable.

Association analysis with CAD was performed for each breed, as the allele frequencies differed substantially between the shiba inu and miniature dachshund samples (Table 3). The rs24378020 G>A (IL4R p. Arg688Cys) was found to be significantly associated with the risk of CAD in miniature dachshunds in the dominant model for minor (non-reference) allele, with an odds ratio of 0.10 (95% confidence interval=[0.01–0.85]; Poriginal=0.0062). This result indicates that rs24378020 in IL4R is associated with CAD susceptibility, and the minor allele A (Cys688) can promote resistance to CAD. No significant association was found between the other two SNPs, rs9193906 and rs22147008, and CAD in allelic frequencies and dominant models. Due to the shortage of dogs homozygous for the minor allele, the association analysis in co-dominant and recessive models could not be performed for both shiba inu and miniature dachshund samples. Table 4 shows the detected haplotypes constructed by the four SNPs and their estimated frequencies. The haplotype distributions were significantly different between shiba inu and miniature dachshund samples (P<0.01, Fisher’s exact test). In total, eight haplotypes were detected; however, no significant combination effect was observed between canine IL4R and IL13 SNPs.

Table 4.

List of haplotypes constructed using four single nucleotide polymorphisms in IL4R and IL13 and their frequencies in two dog breeds

Haplotype	IL4R			IL13	Haplotype frequencies estimation
					Shiba inu		Miniature dachshund
	rs9193906	rs24378020	rs851400460	rs22147008	Control	Case	Control	Case
1 (Reference)	G	G	G	A	0.6379	0.6502	0.1888	0.3037
2	G	G	G	G	0	0.0557	0.2397	0.3805
3	G	A	G	A	0	0	0.1019	0
4	G	A	G	G	0	0	0.0411	0.0263
5	C	G	G	A	0.3103	0.2469	0.1991	0.1174
6	C	G	G	G	0.0517	0.0325	0.1888	0.1721
7	C	G	A	A	0	0.0147	0	0
8	C	A	G	G	0	0	0.0407	0

DISCUSSION

According to a previously reported genome-wide association analysis of CAD, 40 SNPs were identified to be significantly associated with this disease [38]. Among these candidate SNPs, rs24318716 (NC_006588.3: g.1223683T>C), rs24327271 (NC_006588.3: g.5498220C>A), and rs24332727 (NC_006588.3: g.13248759G >A) were located on canine chromosome 6. However, these reported SNPs are located 6 to 18 mega base pairs from the IL4R locus. Therefore, the effect of rs24378020 on IL4R is unlikely to be the result of linkage disequilibrium with previously reported loci.

Currently, more than 200 dog breeds exist worldwide, and each of them is bred separately as pure breeds [10]. As the breeding population for each breed is not very large, they are easily affected by genetic drift. Therefore, allele frequencies are likely to differ between breeds [24, 31]. In this study, the CAD resistant allele in the miniature dachshund (Cys688) was absent in the shiba inu (Table 3), probably due to random fixation. As random fixation can occur independently at any locus in any breed, the analysis of pooled data from multiple dog breeds can weaken the ability to detect associations between genetic polymorphisms and a specific disease. Furthermore, the incidence rate of CAD is significantly different among individual breeds. For example, boxers and French bulldogs are reported to have very high incidence of CAD, while dachshunds and poodles exhibit a significantly lower incidence of CAD [26]. Additionally, it is known that there are substantial breed-associated differences in the clinical phenotype of CAD [36], which may be due to genetic differences among breeds. Thus, for multifactorial diseases such as CAD, it is preferable to perform genetic association analysis within a single dog breed rather than mixing multiple breeds.

This is the first report on the effect of IL4R polymorphism (rs24378020) as a susceptibility factor for CAD. How the IL4RA amino acid substitution p. Arg688Cys affects the function of this receptor remains unclear. However, it may contribute for resistance to CAD if we assume that the minor allele (Cys688) can impair the IL4 and IL13 signals. In humans, blocking IL4RA with dupilumab, a monoclonal antibody against this protein, modulates IL4 and IL13 signals and reduces the levels of Th2 biomarkers [28, 35]. Dupilumab was shown to be very effective in treating refractory human AD [2, 33]. If Cys688 directly affects the function of IL4R, similar effects are expected in dog breeds other than miniature dachshunds. Therefore, additional research is required on other breeds.

In silico analysis predicted that Ala81 of canine IL13 could cause significant changes in protein conformation or function (Table 1). It was previously reported that 78.5% of prediction results by PROVEAN, SIFT, and PolyPhen-2 for the 19,898 disease-associated human variants were in agreement and shared by all three tools [5]. Thus, canine IL13 variant Ala81 was presumed to affect the function of the protein. However, we did not find a significant association between canine IL13 polymorphisms and CAD. Genetic variants of IL13 and IL4R were associated with atopy and asthma risk in humans [4, 15]. In mice, IL13 is necessary for the expression of allergen-induced airway inflammation through a mechanism independent of IgE and eosinophils [37]. Thus, canine IL13 variants p. Thr81Ala and IL4RA p. Arg688Cys may represent candidate genes for allergic disease risk in dogs. Moreover, in humans, the synergistic effects of the IL13 and IL4R variants on IL13-dependent gene induction were previously reported [3, 12]. However, due to the small number of dogs included in this study, the statistical power was insufficient to measure the interaction between IL13 and IL4R polymorphisms. Therefore, large-scale studies are warranted to elucidate the interaction between IL13 and IL4R polymorphisms.

In conclusion, our study suggests that the IL4R Cys688 variant reduces the risk of CAD in miniature dachshunds.