Genome-wide association analysis of avian resistance to Campylobacter jejuni colonization identifies risk locus spanning the CDH13 gene.

The enteropathogen Campylobacter jejuni is a major worldwide health and economic burden, being one of the leading causes of bacterial gastroenteritis and commonly linked to postinfectious onset of autoimmune disease. Chickens are a major vector for human infection and even though variation in avian colonization level is heritable, no previous studies have identified regions of the genome associated with colonization resistance. We performed a genome-wide association study of resistance to C. jejuni colonization in the avian intestine by controlling for population structure, which revealed a risk locus with genome-wide significance spanning the T-cadherin (CDH13) gene. A second possible risk locus was also identified close to calmodulin (CALM1), a calcium-activated modulator of cadherin function. In addition, gene expression analysis of mRNA sequencing profiles revealed that the relative expression of the two genes is significantly associated with colonization resistance. Functional studies have previously demonstrated involvement of cadherins and calmodulin in C. jejuni intracellular invasion and colonization of human intestinal epithelial cells in vitro. Consistent with this finding, our analysis reveals that variation surrounding these genes is associated with avian colonization resistance in vivo and highlights their potential as possible targets for control of the bacterium in avian and human populations.

Campylobacter jejuni is one of the main causes of bacterial gastrointestinal infection worldwide (Scallan ; EFSA 2012) and represents a major public health burden. Campylobacteriosis is usually self-limiting but has been linked to postinfectious onset of the autoimmune diseases Guillain-Barré syndrome (McCarthy and Giesecke 2001; Tam ), Miller Fisher syndrome (Koga ), and reactive arthritis (Hill Gaston and Lillicrap 2003; Townes ). Consumption of contaminated poultry is a common source of human infection. Although pathogenic in humans, C. jejuni is an intestinal commensal in most mammalian and avian hosts and frequently colonizes the chicken intestine at extremely high densities of 108–1010 CFU/g (Van Deun ; Lamb-Rosteski ; Meade ).

Chickens are born free of C. jejuni and only begin to acquire the bacterium at an average age of 2−3 wk (Newell and Fearnley 2003; van Gerwe ). Once one individual in a flock becomes colonized, the bacterium spreads quite rapidly, with >95% of the flock colonized within several days (van Gerwe ). However, C. jejuni colonization levels in the gastrointestinal tract can vary substantially between individual chickens, both within and between populations. Several previous studies have established that this variation in C. jejuni colonization level can be heritable, with significant differences in susceptibility observed between chicken lines (Stern ; Boyd ; Li ). Although a quantitative trait loci mapping study has been referred to (Kaiser 2010), no attempts to identify the genes involved in this resistance by linkage mapping, genome-wide association or candidate-gene analysis have yet been published. The primary site of C. jejuni colonization in the chicken is within the cecum, where it populates the mucus layer overlying epithelial cells (Beery ; Meade ). Although epithelial cell attachment was not documented in vivo (Beery ), attachment and invasion do occur in vitro and are thought to be necessary for successful colonization (Byrne ; Hermans ) and several adhesins involved have been identified (Hermans ). The invasiveness of different strains of C. jejuni into chicken epithelial cells is correlated with systemic colonization (van Deun ), which commonly occurs in the spleen, liver, and bursa of Fabricius (Cox ; van Deun ; Lamb-Rosteski ; Meade ) even though commensal bacteria do not normally disseminate systemically. Both transcellular and paracellular translocation have been described in human epithelial cells in vitro (Konkel ; Grant ; Monteville and Konkel 2002; van Deun ; Hu ; Kalischuk ; Boehm ), and the translocation route in chicken epithelial cells is unknown. Strain invasiveness into human cells is correlated with chicken intestinal colonization potential (Hänel ).

In several studies, researchers have analyzed gene expression in the chicken cecum after C. jejuni infection and have demonstrated a proinflammatory response (Borrmann ; Smith ; Larson ), lymphocyte involvement (Li ; Shaughnessy ; Li ), and activation of Toll-like receptors (de Zoete ), the mitogen-activated protein kinase pathway, and small GTPase-mediated signal transduction (Li ). We have previously described the separation of a single population of 255 4-wk-old Barred Rock chickens into those resistant and susceptible to C. jejuni colonization 48 hr postinfection and have investigated the caecal gene expression profiles of resistant and susceptible chickens by using high-throughput short-read sequencing of mRNA (RNAseq) (Connell ). Resistance was associated with significantly increased expression of genes involved in the innate immune response, cytokine signaling, B-cell and T-cell activation, immunoglobulin production, and the renin-angiotensin system. In this complementary study, a genome-wide analysis of association with colonization status was performed with 194 chickens, selected from both ends of the colonization spectrum and including the 28 transcriptionally profiled birds. These were genotyped on a chicken 60K single-nucleotide polymorphism (SNP) chip (Groenen ). Using a case-control model incorporating population stratification, we identified a significantly associated resistance locus comprising four SNPs in a region of chromosome 11, which suggests involvement of the CDH13 (T-cadherin) gene in resistance to colonization. A second putative locus was also identified on chromosome 5 close to the CALM1 (calmodulin) gene, a mediator of calcium signaling that modulates cadherin function. Investigation of gene expression in cases and controls revealed that the relative expression of these two genes is associated with colonization status. Together, this information strongly suggests a functional role for these genes in establishment of C. jejuni colonization within the avian host.

Materials and Methods

Population history

The population of chickens used in the study was obtained from a line maintained at the University of Saskatchewan, Saskatoon, Canada. The Barred Rock is a dual-purpose breed which was developed in the middle of the 19th century. The line was acquired in the 1920s, and was initially selected but was subsequently maintained unselected from 1965 to 2003. Mild selection for egg production was performed after 2003. The procedure for selection involved maintaining 80 hens per generation. A total of 50 hens and roosters from the same dams were then selected for reproduction based on egg production. Some genetic substructure within the population is therefore expected due to violations of several assumptions of Hardy-Weinberg equilibrium. This population substructure was characterized before association analysis, and tests used to detect significance of association were adjusted to incorporate genetic stratification identified within the population.

Sample selection, SNP genotyping, and quality control

Caecal colonization levels were estimated for all 255 chickens in the population 48 hr postinoculation and caecal samples collected as reported previously (C. jejuni strain NCTC11168v1) (Connell ). Caecal samples from 196 individuals were selected for phenol-chloroform genomic DNA extraction (Sambrook and Russell 2001) and genotyping with the 60K chicken SNP chip (Groenen ). Samples were genotyped commercially using the standard protocol from Illumina Infinium iSelect Beadchips implementing BeadStudio Genotyping v3.0.19.0. Forty birds with no or very minimal colonization (C. jejuni−resistant cases) and 128 birds with very high levels of colonization >107 CFU/g (C.jejuni-susceptible controls) were selected for case−control association analysis (Figure 1). Of the total number of 57,636 SNP positions queried on the chip, 25,385 SNPs were polymorphic in this population. Average identity by state (IBS) was calculated from all SNP genotypes in the R package GenABEL (Aulchenko ). Strict quality control filtering of the data was carried out before association analysis. In association studies of complex diseases, where single alleles are expected to only have a small effect on the trait under study, genotyping errors are more likely to result in false-positive or false-negative findings of association; therefore, SNPs with low minor allele frequency (≤1.8%) were removed. All SNPs with call rate <95% were excluded. Two individuals with call rate <95% were excluded because they may represent poor-quality samples with unreliable genotypes. All SNPs exhibiting substantially higher or lower heterozygosity than expected (P value ≤ 10−7) under Hardy-Weinberg equilibrium were also excluded. Autosomal heterozygosity was also checked to ensure all individuals had moderate autosomal heterozygosity but no individuals had excessively high heterozygosity [false-discovery rate (FDR) <1%]. Genetic relatedness of all possible pairs of individuals was calculated to identify individuals too closely related to be used in the association study. A matrix of IBS values was created from all filtered SNPs. No pair of individuals were >95% IBS. Subsequently, SNPs were filtered so that SNPs in complete linkage disequilibrium (LD) on the same chromosome were represented by only one tag SNP. Each set of SNPs in complete LD was tested as a unit using these tag SNPs.

Figure 1

Colonization status of chickens selected for the genome-wide association study. (A) The initial population of 255 chickens had highly variable colonization levels; 196 of these were selected for genotyping using the 60K chicken SNP chip. Their colonization levels are illustrated in (B). A total of 168 chickens from opposite ends of the colonization spectrum (i.e., all individuals with colonization level ≤400 CFU/g and all individuals with colonization level ≥107 CFU/g) were selected to conduct the case-control association study. Two individuals were filtered due to genotyping call-rate <95%. Colonization levels of all 166 individuals which passed quality control filtering for the case-control study are plotted in (C). Forty of these had no or very minimal colonization (resistant cases), and 126 had high colonization levels (susceptible controls). Birds for which RNAseq expression profiles are available are highlighted in red at the extremes of the distribution. The lower limit of detection = 400 CFU/g.

Detection of population stratification and case−control association analysis

To investigate any population stratification that may be present in the population, a matrix of genomic kinship was generated from all pairs of 194 individuals in the population that passed quality control filtering using genotypes of all filtered tag SNPs. This matrix of genomic kinship was converted to a distance matrix that was used to carry out classical multidimensional scaling analysis. All analysis of population stratification was conducted in GenABEL (Aulchenko ). The qtscore function in GenABEL was implemented to carry out logistic regression of the trait onto SNP genotypes using the first three principal components of the genomic kinship matrix as covariates in the regression model, accounting for the genetic substructure within the population (Price ). Overinflation of the test statistic was investigated by comparison of the distribution of the test statistic to that expected under the null hypothesis [inflation factor λ (Amin )]. Consistent deviation is indicative of inflation due to population substructure.

Gene expression analysis

Read counts of gene expression were first corrected using upper quartile normalization and relative log expression in edgeR (Robinson ). No correlation between SNP genotypes and relative expression levels was found. Similarly, no correlation between SNP genotypes in both genes was found. Significant differences in the ratios of CDH13/CALM1 expression (Figure 6) were observed between cases and controls (Mann-Whitney U-test, P = 0.004). All calculations were performed in R using the “lm” and “wilcox.test” functions.

Figure 6

Relative expression levels of CDH13 and CALM1. Boxplots of resistant vs. susceptible birds reveal separation of the population based on CDH13/CALM1 expression (Mann-Whitney U-test, P = 0.004) and no bias in relative ratios due to sequencing depth. The ratio of CDH13/CALM1 expression for all 28 transcriptionally profiled birds is plotted, with shading reflecting depth of sequencing.

Results

Genotyping and quality control

Querying of the study population of 196 chickens using a 60K Illumina iSelect chip [37] revealed that of 57,636 SNP positions, just 25,385 were polymorphic in this population. All SNP genotypes and colonization status of individuals are detailed in Supporting Information, File S1 and File S2, respectively. Average IBS calculated using all SNP genotypes was 91.5% for all pairs of individuals, attesting to the high inbreeding present in the population due to its population history and maintenance as an inbred population since the 1920s. A total of 16,871 SNPs and 194 chickens passed quality control filtering (call-rate >95%, minor allele frequency < 1.8%, Hardy-Weinberg P-value < 10−7). After filtering out SNPs on the same chromosome in complete LD, 5609 nonredundant tag SNPs remained and were used in the analysis of population substructure and genome-wide association testing.

Multidimensional scaling reveals population stratification

Population substructure can induce spurious associations in genetic association studies. It is therefore essential that population stratification is identified and adjusted for, especially in genome-wide association studies of complex traits in which SNPs may have only a small effect on the trait being analyzed. To characterize genetic substructure within the population, genomic kinship was investigated using the R package GenABEL (Aulchenko ). First, a matrix of IBS values was calculated between all pairs of 194 individuals which passed quality control filtering using genotypes of the 5609 tag SNPs. Classical multidimensional scaling analysis was performed on the distance matrix derived from this matrix of kinship values. The major stratification within the population was resolved by the first three dimensions and these are plotted in Figure 2. This plot reveals some substructure within the population which is most likely due to breeding procedures, particularly those used to select for increased egg production (see Materials and Methods). Resistant cases and susceptible controls are distributed randomly throughout the plots indicating that they do not cluster together genetically.

Figure 2

Principal components revealing population stratification. All 194 individuals that passed quality control filtering are represented. Colonization-susceptible individuals from the case-control study are shaded blue and colonization-resistant birds shaded red. Gray circles represent individuals from the center of the colonization spectrum that were not used in the association study. All autosomal markers were used to construct a kinship matrix of genetic relatedness or identity by descent (IBD). This matrix was converted to a distance matrix and classical multidimensional scaling performed to detect any population substructure. There is evidence for some stratification of the population which was resolved by the first three principal components (PC1, PC2, and PC3). The percentage variation which is represented by each of these principal components is 18.6% (PC1), 16.8% (PC2), and 11.1% (PC3). The first two principal components reveal separation of the population into two subpopulations—one large group with 157 individuals and one smaller group of 37 individuals.

Case−control study identifies loci associated with colonization resistance

The qtscore function in GenABEL was applied to perform logistic regression of colonization status onto SNP genotypes. Population stratification was corrected for by using the first three principal components of the kinship matrix as covariates in the regression model. This yields a genomic inflation factor λ = 1.06 (compared with λ = 1.42 without adjustment for population stratification) indicating that the major substructure has been adjusted for. The quantile-quantile plot of observed P-values vs. the expected null distribution (Figure 3) did not reveal deviation from the null, apart from 12 SNP outliers at the tail of the distribution exhibiting the most significant P-values, which all lie within three regions of the genome—on chromosomes 11, 5, and 1. When these SNPs are not included the genomic inflation factor λ = 0.995. Combined, this finding suggests that these SNPs may represent truly associated regions of the genome and deviations are not due to over-inflation of the test statistic.

Figure 3

Quantile-Quantile plot of observed vs. expected P-values. The q-q plot of observed P-values shows they follow the expected null distribution apart from 12 SNPs with the most significant P-values. These are located within three regions: chromosomes 11 (red squares), 5 (blue triangles), and 1 (yellow circles).

The genome-wide plot of SNP association significance is shown in Figure 4, where these three regions are again visible outliers in this plot, displaying high association P-values. The most significant associations identified involved SNPs on chromosome 11, where four SNPs reached genome-wide significance, after we controlled for multiple testing (FDR = 0.04), and odds ratios of 2.33−2.50 (Table 1). These SNPs are all in high LD, with two adjacent SNPs being in complete LD. They are located within a region spanning 80 kbp, are intergenic, and surround the CDH13 gene on chromosome 11 (Figure 5).

Figure 4

Manhattan plot of genome-wide association. All 16,871 SNPs that passed quality control are shown. Genomic location is plotted against −log10(P). Four SNPs on chromosome 11 reached genome-wide significance (FDR = 0.04).

Table 1

Case-control study: top SNPs

SNP ID	Chr	Position	N	Allele Frequency	Odds Ratio (95% Confidence Interval)	P	FDR
Gga_rs15622247	11	17594419	166	0.40/0.21	2.50 (1.45−4.30)	1.61E-05	0.0426
GGaluGA079086	11	18376977	166	0.40/0.22	2.33 (1.35−3.99)	1.78E-05	0.0426
Gga_rs15622362	11	17655435	166	0.40/0.21	2.50 (1.41−4.19)	2.28E-05	0.0426
Gga_rs14027234	11	17657890	166	0.40/0.21	2.50 (1.41−4.19)	2.28E-05	0.0426

SNPs exhibiting the most significant association with colonization status are listed, along with chromosome position. Four SNPs exhibited genome-wide significance after FDR control for multiple testing (FDR = 0.04). Minor allele frequencies in cases/controls and odds ratios with 95% confidence intervals are shown. Two SNPs (Gga_rs15622362 and Gga_rs14027234) are in complete LD. SNP, single-nucleotide polymorphism; FDR, false-discovery rate; LD, linkage disequilibrium.

Figure 5

Regional plot of significant association on chromosome 11. Four SNPs with genome-wide significance FDR < 0.05 are highlighted red. These four SNPs are intergenic but surround the T-cadherin (CDH13) gene. Two of these SNPs are in complete LD, situated very close to each other, and are almost indistinguishable on the plot. The physical positions of all Refseq genes within the region are depicted below the plot.

Association with gene expression profiles

We have previously analyzed the gene expression profiles of 28 of these chickens (14 resistant cases with no C. jejuni colonization and 14 controls with the highest C. jejuni colonization levels) through mRNA sequencing (i.e., RNAseq) (Connell ). All data from RNA sequencing gene expression analysis, including raw sequence files, can be downloaded from NCBI GEO accession no. GSE44341. As the second most significant region identified in this study approaches significance and surrounds another Ca2+-regulated gene, CALM1, associations between gene expression and SNP genotypes in both regions were investigated. No correlation with SNP genotypes was detected. Neither of these genes was differentially expressed between cases and controls. However relative CDH13/CALM1 gene expression is significantly associated with resistance to colonization (Mann-Whitney U-test, P = 0.004), with resistant birds displaying lower CDH13 and higher CALM1 relative expression levels (Figure 6).

Discussion

Several studies have demonstrated that resistance to C. jejuni colonization of the chicken gastrointestinal tract can be heritable, but this is the first to identify specific regions of the genome associated with colonization resistance. Association analysis treating colonization status as a binomial trait was performed and identified a resistance locus for C. jejuni colonization on chromosome 11 containing four significantly associated SNPs. This locus centers on the CDH13 (T-cadherin) gene, suggesting a functional role for this gene in the colonization process. In previous studies, authors have demonstrated a direct interaction between C. jejuni HtrA protease and another member of the cadherin superfamily, E-cadherin (cadherin 1, CDH1), in which cleavage of the protein facilitates transmigration across polarized epithelial cells in vivo (Hoy ; Boehm ). Because T-cadherin is expressed on the apical surface of chicken intestinal epithelial cells in vivo (Koller and Ranscht 1996) and has a highly similar ectodomain to E-cadherin, the interaction of T-cadherin with C. jejuni in the intestinal lumen is highly likely and variations altering the function or expression of this gene may plausibly affect caecal colonization status.

Cadherins are a superfamily of calcium-dependent proteins with prominent roles in homophilic cell−cell adhesion and maintenance of structural and functional tissue integrity (Halbleib and Nelson 2006). Classical cadherins are composed of an ectodomain containing 5 extracellular cadherin repeats, a transmembrane anchor, and an intracellular domain. The most extensively studied of the cadherins, E-cadherin (cadherin 1, CDH1), is a cell adhesion molecule and tumor suppressor that localizes to adherens junctions on the basolateral membrane of epithelial cells. Despite its typical location below tight junctions, it interacts with several bacteria, facilitating internalization and colonization. The protein can act as a receptor for Listeria monocytogenes and Candida albicans, allowing these pathogens to gain access into the cell (Mengaud ; Phan ). Cleavage of the protein is also a common feature of epithelial cell invasion. Candida albicans (Frank and Hostetter 2007; Villar ), Helicobacter pylori (Hoy ; Hoy ), Shigella flexneri (Sansonetti ; Hoy ), Porphyromonas gingivalis (Katz ; Katz ), Bacteroides fragilis (Wu ), and enteropathogenic Escherichia coli (Hoy ) all cleave E-cadherin on the surface of human epithelial cells, leading to loss of the adherens junction and cell−cell adhesion, thus disrupting the epithelial barrier. Streptococcus pneumoniae (Anderton ) and Clostridium botulinum (Sugawara ) also produce proteins that bind and disrupt E-cadherin−mediated cell−cell adhesion without cleavage. A direct interaction with C. jejuni HtrA protease has recently been shown in which E-cadherin is cleaved by the bacterial protease, resulting in inhibition of intercellular epithelial adhesion (Hoy ). Deletion of C. jejuni HtrA leads to loss of E-cadherin cleavage and highly defective invasion and paracellular transmigration (Hoy ; Boehm ).

Inactivation of E-cadherin in the mouse intestine has indicated that the protein plays a crucial role in maintenance of homeostasis not only through development of epithelial cell adherens junctions and desmosomes but also in maturation and localization of Paneth and goblet cells (Schneider ). These cells are involved in the first line of host defense against pathogenic invasion and secrete a variety of protective antimicrobial molecules into the intestinal lumen in response to pathogenic challenge (Bevins and Salzman 2011; Gill ; Forman ). Inactivation of E-cadherin subsequently results in increased susceptibility to intestinal bacterial infection (Schneider ).

T-cadherin is a unique member of the cadherin family which lacks the highly conserved cytoplasmic and transmembrane domains, and attaches to the cell membrane through a glycosylphosphatidylinositol anchor (Philippova ). In vivo, T-cadherin is expressed on the apical surface of chicken intestinal epithelial cells (Koller and Ranscht 1996) (unlike classical cadherins, which are expressed basolaterally). The ectodomain (extracellular domain) shares high sequence homology with those of classical cadherins such as E-cadherin. It also lacks the cytoplasmic region, which is usually required for classical cadherin-mediated cell−cell adhesion but induces homophilic cell−cell adhesion through an alternative binding mechanism (Vestal and Ranscht 1992; Ciatto ). The protein may play a major role in signal transduction due to its location within lipid rafts which regulate signal transduction (Simons and Toomre 2000). Due to its expression on the apical surface of polarized epithelial cells at the interface with the intestinal lumen and the high sequence homology of its ectodomain with that of E-cadherin, we hypothesize that it may interact with C. jejuni and consequentially affect resistance to colonization either by facilitating internalization or inciting a protective host response.

Although the locus exhibiting the second-most significant association with colonization resistance did not reach genome-wide significance after multiple testing adjustment (FDR = 0.24), it localized close to the calmodulin (CALM1) gene, which encodes a modulator of cadherin-mediated cell-cell adhesion (Li ). Calmodulin is a mediator of calcium signaling that binds calcium Ca2+ ions and subsequently regulates the activity of calcium-dependent enzymes, including cadherins, involved in a wide variety of cellular processes. An increase in cytosolic Ca2+ concentration is associated with intracellular invasion and pathogenicity of a wide range of bacteria including C. jejuni (Pace ; Hu ; Gekara ; Kim ; Bandyopadhaya ; Asmat ; Konar and Ghosh 2012). Invasion of human intestinal epithelial cells by C. jejuni in vitro requires release of Ca2+ from intracellular stores and inhibition of calmodulin using a calmodulin agonist inhibits C. jejuni intracellular epithelial cell invasion in vitro (Hu ). This finding suggests there may be a relationship with the association observed surrounding the T-cadherin gene, although no correlation between SNP genotypes from the two regions was found. This association study was carried out subsequent to a study on a subset of the birds analyzed here which characterized differences in gene expression between resistant and susceptible chickens. Correlation between the SNP genotypes observed here and gene expression levels was not found and neither of these genes displayed significantly differential expression. However, a significant correlation between colonization status and the ratio of expression levels of the two genes is seen. Interestingly, the gene displaying the most significant differential expression in this study, which demonstrated fourfold higher expression in resistant birds (p = 5.6e−20), encodes an epithelial calcium-activated chloride channel (LOC424523)(Connell ).

To our knowledge, this is the first study to identify a gene significantly associated with resistance to C. jejuni colonization of the chicken. The study was based on a population of Barred Rock chickens, a dual-purpose breed that originated in the United States in the 1800s and one of the founder breeds of the broiler industry. The breed has contributed significantly to the genetics of modern commercial broiler and brown egg layer breeds, most notably in extensive use of the White Plymouth Rock, which was developed from the Barred Rock, in maternal commercial broiler grandparent lines (Crawford 1990). Due to the high inbreeding present within the population, the likelihood of shared mechanisms of resistance to colonization and repetition of these findings with commercial populations is uncertain. The resistance locus identified centers on T-cadherin, and although of marginal significance and requiring further validation, suggests this gene plays a prominent role in the colonization process. This is supported by the importance of the related protein, E-cadherin in epithelial cell invasion of several bacteria, including C. jejuni and in the maintenance of intestinal homeostasis. A putative risk locus centering on the calmodulin gene, which can regulate cadherin function, and previous evidence for differential expression of a calcium activated chloride channel associated with colonization provide additional evidence for Ca2+-regulated cadherin involvement. C. jejuni is one of the most commonly reported causes of food-borne illness worldwide and chickens are a major vector for infection. We believe the results reported here highlight the significance of cadherins, and specifically T-cadherin, in control of the bacterium within the chicken cecum.