Genetic control of Campylobacter colonisation in broiler chickens: genomic and transcriptomic characterisation

Campylobacter is the leading cause of bacterial foodborne gastroenteritis in many countries. Source attribution studies unequivocally identify the handling or consumption of contaminated poultry meat as the primary risk factor. One potential strategy to control Campylobacter is to select poultry with increased resistance to colonisation. We conducted genomic and transcriptomic analyses of commercial pedigree broilers exposed to Campylobacter to examine persistent colonisation of the caecum as a quantitative trait. 3,000 broilers were genotyped using a 50K single nucleotide polymorphism (SNP) array and imputed to 600K SNPs. Genotypes were analysed for associations with the number of viable Campylobacter in the caeca. Heritability of the trait was modest but significantly greater than zero (h2=0.11 ± 0.03). Genome-wide association analyses confirmed quantitative trait loci (QTL) on chromosomes 14 and 16 previously identified using the progeny of crosses of inbred lines differing in resistance, and detected two additional genome-wide significant QTLs on chromosomes 19 and 26. RNA-Seq analysis of the transcriptome of caecal tonsils from birds at the low and high extremes of C. jejuni colonisation phenotype identified differentially transcribed genes, mainly located within the QTL on chromosome 16 and proximal to the major histocompatibility complex (MHC) locus. We also identified strong cis-QTLs located within the MHC suggesting the presence of cis-acting variation in both MHC class I, class II and BG genes. Multiple other cis-acting variants were identified in association with key immune genes (COPS3, CCL4, CR1L, C4BP, PLGR) in the other QTLs. Pathway and network analysis implicated cooperative functional pathways and networks in colonisation, including those related to antigen presentation, innate and adaptive immune responses, calcium, and renin-angiotensin signalling. While co-selection for enhanced resistance and other breeding goal traits is feasible, the frequency of resistance-associated alleles was high in the population studied and non-genetic factors significantly influence Campylobacter colonisation in poultry. Author summary Campylobacter infection is estimated to cause 95 million illnesses in people worldwide each year. Human infections mostly involve gastroenteritis, but can have severe complications. The handling or consumption of contaminated poultry meat is a key risk factor for human campylobacteriosis. The bacteria reach high numbers in the intestines of chickens reared for meat (broilers) and are frequently found on carcasses after slaughter. Effective vaccines against Campylobacter are not yet available, and treatments to reduce carcass contamination (e.g. chlorination) are not acceptable in some markets. One alternative is to breed for chickens with improved resistance to Campylobacter colonisation. To test the feasibility of this option in commercial birds, we analysed the genetic make-up of 3,000 pedigree broilers and determined the number of Campylobacter in their gut. There were associations between specific regions of the chicken genome and resistance to Campylobacter. Within some of these regions, expression of certain genes differed between birds at the low and high extremes of Campylobacter colonisation, providing a potential explanation for genetic variation in resistance. Selection of poultry with increased resistance to Campylobacter colonisation may be a complementary strategy to improved biosecurity, management, handling and processing procedures to reduce the burden of Campylobacter on human health.


68
Human campylobacteriosis exerts profound societal and economic costs. The World Health Organisation 69 estimated that Campylobacter caused 95 million illnesses, 21,000 deaths and loss of 2.1 million disability-70 adjusted life years globally in 2010 [1]. In the United Kingdom alone, there were 63,946 laboratory-confirmed 71 human cases in 2017, most of which were due to C. jejuni The ratio of unreported cases of human 72 campylobacteriosis to those captured by national surveillance in the United Kingdom is estimated to be 9.3 73 to 1 [2] therefore the true burden may exceed half a million cases per annum at a projected median cost to 81 relatively modest 2 log 10 reduction in the number of Campylobacter on broiler carcasses could reduce the 82 incidence of human disease due to contaminated chicken by up to 30-fold [10]. A pressing need therefore 83 exists for strategies to reduce the entry of Campylobacter-contaminated poultry meat into the food chain.

84
As effective vaccines and treatments for pre-slaughter control of Campylobacter in poultry are lacking, 85 much interest exists in the potential for breeding chickens with improved resistance to intestinal colonisation 86 by C. jejuni. The scope for genetic control of Campylobacter in poultry has been demonstrated in commercial 87 broiler lines that vary in resistance [11][12][13][14][15][16][17] and inbred layer lines exhibiting heritable differences in C. jejuni 88 colonisation following experimental inoculation [18,19]. Breeding for Campylobacter resistance may also 89 benefit avian intestinal health and productivity. In some commercial broiler chicken lines, experimental 90 inoculation of C. jejuni elicits damage to the intestinal mucosa, diarrhoea and impaired weight gain [20,21]

97
A previously published genome-wide association study (GWAS) on C. jejuni intestinal colonisation, where 98 phenotypes were analysed as a binary trait in a novel dual-purpose chicken breed, revealed a resistance-99 associated locus on chromosome 11 near the CDH13 gene encoding T-cadherin, and a second candidate locus 100 on chromosome 5 was identified close to calmodulin (CALM1), a calcium-activated modulator of cadherin 101 function [17]. Both genes differed in relative expression in a manner associated with resistance [17]. Studies 102 in inbred layer lines 6 1 and N found heritable differences in caecal colonisation by various C. jejuni strains 103 [18]. Initial low-resolution mapping studies using reciprocal (6 1 ♀×N♂) and (N♀×6 1 ♂) F1 crosses and the 104 progeny of an (N♀x6 1 ♂) F1♂ x N♀ backcross indicated that resistance was controlled by a single autosomal 105 dominant locus [18], but subsequent analysis of a backcross population using 1,243 fully-informative single 106 nucleotide polymorphism (SNP) markers revealed quantitative trait loci (QTL) on chromosomes 7, 11 and 14 107 [19]. Using a ninth-generation advanced intercross (6 1 xN) line and a 600K genome-wide SNP array, the 108 location of the QTL on chromosome 14 was confirmed and refined, and additional QTLs were identified on 109 chromosomes 4 and 16, indicating potential involvement of the Major Histocompatibility Complex (MHC) 110 region [19]. QTL for resistance of chickens to enteric carriage of Salmonella have been reported at the same 111 regions of chromosome 14 [27] and 16 [27][28][29] that were implicated in host resistance to C. Jejuni 112 colonisation. Analysis of caecal gene expression in chicken lines of varying susceptibility to Campylobacter 113 colonisation has identified transcriptional signatures associated with differential resistance, including genes 114 influencing the immune response [13][14][15][16].

115
In the present study, we conducted a comprehensive genome-wide association study on a commercial 116 pedigree broiler population exposed to Campylobacter naturally present in the litter. The genomic 117 architecture of resistance was analysed using imputed high-density SNP genotyping, and resistance was 118 further characterised by transcriptome analyses of intestinal tissue from birds of the predicted resistant or 119 susceptible genotypes at the corresponding extremes of caecal Campylobacter colonisation phenotype. Our 120 data provide valuable insights into the prospects for genetic control of Campylobacter in poultry.

123
Descriptive statistics and genetic parameters affecting Campylobacter levels 124 The mean number of caecal C. jejuni detected ± standard deviation was 7.057 ± 1.023 log 10 colony-forming 125 units (CFU) per gram (g), with a maximum of 10.64 log 10 CFU/g and minimum of 2 log 10 CFU/g, equal to the 126 limit of detection by direct plating. Sex had a marginal effect on Campylobacter levels (P <0.05), with males 127 having higher Campylobacter load (7.178 ± 0.034 log 10 CFU/g) compared to females (6.930 ± 0.032 log 10  Table).

202
Differential gene expression analysis. After false discovery rate (FDR, P<0.05) correction for multiple 203 testing, 3 protein-coding genes were found to exhibit significant differential expression ( In order to identify more subtle patterns of 207 differential expression, a relaxed significance threshold of unadjusted P value of 0.001 was implemented and 208 a total of 33 genes exhibited differential expression between high-, average-, and low-colonised birds at this 209 threshold (  Table). Of those, 90 were associations between SNPs in high LD, located in the same 223 QTL region on chromosome 16, and the expression of a single gene, BG-like antigen 1 (BG1) (Fig. 4A Table).
230 Trans-analysis. We detected a total of 13 significant trans-eQTLs within the QTL region on chromosome 231 19 and two on chromosome 26 (S5 Table). Most of these predicted trans-acting elements are for genes 232 related with metabolic processes. The trans-acting SNP on chromosome 26 is for a microRNA (gga-mir-1553) 233 located on chromosome 7, close to the peak of a previously identified QTL for C. jejuni resistance identified 234 using a back-cross population of inbred lines 6 1 and N [19].

235
Allele-specific expression analysis. If an individual is heterozygous for a cis-acting SNP it is expected that 236 the two alleles of the gene will be expressed unequally causing allelic expression imbalance. To verify the cis-

237
QTLs detected above, and identify additional ones, we identified genetic variation within the QTL regions 238 identified by the GWAS and RHM using the RNA-Seq data and performed allele-specific expression (ASE) 239 analysis for all the SNPs located within these regions. Several significant ASE events were identified in all QTL 240 regions (mean P value ≤ 0.05 with at least 4 heterozygous animals). Specifically, 14 significant ASE events 241 were identified for 3 genes located on chromosome 14, 30 for 6 genes on chromosome 16, 11 for 2 genes 242 and one microRNA (gga-mir-142) on chromosome 19, and 35 for 5 genes on chromosome 26 (S6 Table). A 243 highly significant ASE event was identified on chromosome 14 (for the QTL located at 12MB) for chloride   (Fig. 6). Moreover, three networks of molecular interactions related to 'immunological diseases', 'cell death and survival', and 'molecular transport and protein trafficking' were constructed using the list of 270 genes in the candidate regions (Fig. 7). We subsequently extracted the gene ontology terms for each of these 271 genes and performed functional annotation clustering analysis. The genes were organised into 41 clusters, 272 each given an enrichment score (ES). The first (ES = 4) and the second (ES = 3.5) clusters were both enriched 273 for genes functionally annotated as involved in 'antigen processing and presentation via MHC class I and class 274 II molecules' (including BF1, BF2, BLB1, BLB2, DMB1, DMB2, TAP1, and TAP2) (S7 Table). Bonferroni correction for multiple testing, significance thresholds for analysis with the 50K array were P ≤ 494 1.33 × 10 −6 and P ≤ 2.66 × 10 −5 for genome-wide significant levels (i.e., P ≤ 0.05) and suggestive significant 495 levels (namely one false positive per genome scan), respectively, corresponding to −log 10 (P) of 5.87 and 4.47.

496
The significance thresholds for the 600K array after Bonferroni correction were P ≤ 1.73 × 10 −7 and P ≤ 3.46 effects against a base model that excluded the latter effect. Only the 50K data was analysed with this method.

518
A total of 1,915 regions were tested across the genome. After adjustment for multiple testing, using the 519 Bonferroni correction, significance thresholds were P ≤ 2.63×10 −5 and P ≤ 5.26×10 -5 for genome-wide levels 520 (P ≤ 0.05) and suggestive levels (namely one false positive per genome scan), respectively, corresponding to 521 −log 10 (P) of 4.57 and 3.27.

526
Moreover, all the genes that were located within the 20 SNP windows found to be significant by RHM; and 527 the 250 kb 5' and 3' regions of the significant SNP markers identified by the GWAS were annotated using 528 GalGal5 data obtained by the BioMart data mining tool (http://www.ensembl.org/biomart/martview/). This 529 allowed us to catalogue all the genes that were located in the vicinity of the identified significant SNP markers 530 for Campylobacter colonisation resistance.

532
Transcriptomic analyses 533 Total RNA was prepared from the caecal tonsils of 23 broilers, selected on their genotype (allele 534 combination in the significant identified markers) and caecal Campylobacter load, after correction for other 535 sources of systematic variation (sex and date of sampling). Details of the birds selected are shown in S3 Table. 536 RNA was extracted using the RNeasy Mini Kit (Qiagen Hilden, Germany) according to manufacturer's 537 instructions. The resultant RNA was checked for quality using the Agilent Tapestation 2200, and all samples   538 were of high quality with RNA Integrity Numbers (RIN) greater than 9. Library preparation was performed by Fisher's exact test (IPA score or P score = -log 10 (P value)).

621
The gene list for Campylobacter colonisation resistance was also analysed using the Database for   with Campylobacter colonisation resistance in commercial chickens. All the significant SNP markers were 1000 in high LD, illustrated with red colour, and were located within the same LD block (230kb) marked with 1001 triangle.