Comprehensive Longitudinal Microbiome Analysis of the Chicken Cecum Reveals a Shift From Competitive to Environmental Drivers and a Window of Opportunity for Campylobacter.

Chickens are a key food source for humans yet their microbiome contains bacteria that can be pathogenic to humans, and indeed potentially to chickens themselves. Campylobacter is present within the chicken gut and is the leading cause of bacterial foodborne gastroenteritis within humans worldwide. Infection can lead to secondary sequelae such as Guillain-Barré syndrome and stunted growth in children from low-resource areas. Despite the global health impact and economic burden of Campylobacter, how and when Campylobacter appears within chickens remains unclear. The lack of day to day microbiome data with replicates, relevant metadata, and a lack of natural infection studies have delayed our understanding of the chicken gut microbiome and Campylobacter. Here, we performed a comprehensive day to day microbiome analysis of the chicken cecum from day 3 to 35 (12 replicates each day; final n = 379). We combined metadata such as chicken weight and feed conversion rates to investigate what the driving forces are for the microbial changes within the chicken gut over time, and how this relates to Campylobacter appearance within a natural habitat setting. We found a rapidly increasing microbial diversity up to day 12 with variation observed both in terms of genera and abundance, before a stabilization of the microbial diversity after day 20. In particular, we identified a shift from competitive to environmental drivers of microbial community from days 12 to 20 creating a window of opportunity whereby Campylobacter can appear. Campylobacter was identified at day 16 which was 1 day after the most substantial changes in metabolic profiles observed. In addition, microbial variation over time is most likely influenced by the diet of the chickens whereby significant shifts in OTU abundances and beta dispersion of samples often corresponded with changes in feed. This study is unique in comparison to the most recent studies as neither sampling was sporadic nor Campylobacter was artificially introduced, thus the experiments were performed in a natural setting. We believe that our findings can be useful for future intervention strategies and help reduce the burden of Campylobacter within the food chain.

Introduction

Chickens (Gallus gallus domesticus) are an important food source for humans with over 50 billion reared annually for meat and eggs (Part et al., 2016). Feed conversion and the health of chickens is heavily dependent on the largely unexplored complex gut microbial community which plays a role in nutrient assimilation, vitamin and amino acid production and prevention of pathogen colonization (Józefiak et al., 2004; Apajalahti, 2005; Mcnab, 2007; Sergeant et al., 2014). In chickens, the organ with the highest number and variety of bacteria is the cecum (1010-1011 cells/g) which plays an essential role in the digestion of non-starch polysaccharides (NSPs) found in chicken feed (Barnes et al., 1972; Józefiak et al., 2004; Bjerrum et al., 2006). The importance of this organ is demonstrated when up to 10% of energy needs can be recovered from a well-functioning cecum (Hegde et al., 1982; Józefiak et al., 2004). The cecum remains a source of bacterial human infection and a reservoir of antibiotic resistance determinants.

The chicken cecum contracts several times a day releasing contents toward the ileum and the cloaca (Pauwels et al., 2015). Notably the cecal drop contains Campylobacter, a Gram-negative spiral shaped bacterium which causes an estimated 400 million human infections each year (Friedman et al., 2000; Walker, 2005). Campylobacter causes bloody diarrhea, fever and abdominal pains in humans and can also cause post infectious sequelae such as Guillain-Barré syndrome which is a potentially fatal paralytic autoimmune illness. In low-resource areas, asymptomatic and occasionally persistent Campylobacter infections are common in children younger than 1 year and correlate with stunted growth and therefore life-long physical and cognitive deficits (Amour et al., 2016). Approximately 80–90% of these infections are attributed to Campylobacter jejuni, with poultry as the most important source of human campylobacteriosis within industrialized countries (Humphrey et al., 2007; Mullner et al., 2009; Sheppard et al., 2009). C. jejuni colonizes the chicken cecum with relatively high numbers (109 CFU per gram) and whereas traditionally was considered a commensal of the chicken gut, more recently has been demonstrated to be pathogenic to the chicken, with this dependent on the genetics of the host and the strain of infection (Van Deun et al., 2008; Hermans et al., 2012; Humphrey et al., 2014, 2015; Wigley, 2015). Natural colonization of chickens is reported to be at approximately day 14 of the chicken life cycle, although we do not know how and why this occurs, and what the impact of Campylobacter is on the microbiome (Neill et al., 1984; Kalupahana et al., 2013; Thibodeau et al., 2015).

The microbiome of chickens develop rapidly from days 1–3 where Enterobacteriaceae dominate, with Firmicutes increasing in abundance and taxonomic diversity from approximately day 7 onwards (Danzeisen et al., 2011; Ballou et al., 2016; Mancabelli et al., 2016). Bacterial populations within the chicken gut are subsequently driven by the rearing environment and from the bacteria present in food and water (Connerton et al., 2018). How and when Campylobacter appears and the impact on the chicken gut microbiome remains unanswered. The presence of Campylobacter has been noted to prompt an increase in Bifidobacterium and modify abundances of Clostridia and Mollicutes (Thibodeau et al., 2015). The identification of a number of hydrogenases within the ceca may lead to a potential hydrogen sink and provide an explanation as to the high abundance of genera such as Campylobacter (Sergeant et al., 2014). Comparison of broilers not exposed and exposed to C. jejuni at day 6 or day 20 revealed reductions in the relative abundance of operational taxonomic units (OTUs). These were within the taxonomic families Lactobacillaceae and the Clostridium cluster XIVa, with specific members of the Lachnospiraceae and Ruminococcaceae families exhibiting transient shifts in microbial community populations dependent upon the age at which the birds become colonized by C. jejuni (Connerton et al., 2018). These studies have enhanced our understanding of the chicken cecal microbiome, however the lack of day to day microbiome data, suitable replicate numbers, relevant metadata, and lack of natural infectivity studies have not allowed us to fully appreciate what is occurring in a natural habitat in relation to how and when Campylobacter appears within the chicken gut. To answer these questions, in this study we have performed a comprehensive analysis of the chicken cecal microbiome from days 3 to 35, with 12 replicates per day (final n = 379), correlating additional metadata such as chicken weight and feed conversion rates with Campylobacter detection in a natural environmental setting.

Materials and methods

Ethics statement

Approval to conduct the experiments were granted by Agri-Food and Biosciences Institute (AFBI) Establishment License 5002 for AFBI Veterinary Science Division. Euthanasia of birds were carried out by methods laid out in Animal Scientific Protection Act (ASPA) schedule 1.

Experimental design, broilers and sample collection

This study was performed using a total of 396 Ross-308 male broiler chickens provided by Moy Park (39 Seagoe Industrial Estate, Portadown, Craigavon, Co. Armagh, BT63 5QE, UK). The birds were divided into 12 pens; each pen contained 33 chickens (Supplementary Figure 1). Birds were raised on three phase diets from day 0 to day 35. Starter diets were offered to the birds from days 0 to 10, grower diets from days 11 to 25 and finisher diets from days 26 to 35. Every 24 h, a single chicken from each of the 12 pens was removed at random, and euthanized according to ASPA schedule 1 guidelines. Briefly, birds under 250 g were euthanized by dislocation of the neck whereas those over 250 g and up to 1 kg were euthanized by dislocation of the neck following anesthesia using isoflurane. Birds over 1 kg were euthanized by an overdose of anesthetic (isoflurane) followed by dislocation of the neck. Anesthesia was carried out using an anesthetic mask fitted over the bird's head to deliver the vapourised isoflurane with oxygen with death confirmed in all birds by the onset of rigor mortis. Following this, genomic DNA (gDNA) was extracted from the chicken cecum. Out of 396 samples, a total of 17 were removed from the final analysis due to poor gDNA quality giving a final n = 379.

Poultry growth and performance measurements

The performance parameters investigated were mean body weight (BW_mean), body weight gain (Gain), feed intake (FI) and feed conversion ratio (FCR). Measurements were taken at time points 3–7 days, 8–14 days, 15–24 days, and 25–35 days. These variables were then correlated with the microbial community's composition in various statistical analyses.

DNA extraction, 16S rRNA amplification and sequencing

Cecal gDNA was extracted using the QIAamp DNA Stool Mini Kit according to the manufacturer's instructions and stored at −20°C. 16S metagenomic sequencing library construction was performed using Illumina guidelines (Illumina, U.S.A). The 16S ribosomal primers used were V3 (tcgtcggcagcgtcagatgtgtataagagacagcctacgggnggcwgcag) and V4 (gtctcgtgggctcggagatgtgtataagagacaggactachvgggtatctaatcc) (Klindworth et al., 2013; D'Amore et al., 2016). A second PCR step was performed to attach dual indices and Illumina sequencing adapters using the Nextera XT Index kit. Sequencing was performed on the Illumina MiSeq at LSHTM using a v3 300 bp paired-end kit.

Bioinformatics

Abundance tables were obtained by constructing OTUs (a proxy for species) as follows. Paired-end reads were trimmed and filtered using Sickle v1.200 (Joshi and Fass, 2011) by applying a sliding window approach and trimming regions where the average base quality drops below 20. Following this we applied a 10 bp length threshold to discard reads that fall below this length. We then used BayesHammer (Nikolenko et al., 2013) from the Spades v2.5.0 assembler to error correct the paired-end reads followed by pandaseq (v2.4) with a minimum overlap of 20 bp to assemble the forward and reverse reads into a single sequence spanning the entire V3-V4 region. The above choice of software was as a result of author's recent work (Schirmer et al., 2015; D'Amore et al., 2016) where it was shown that the above strategy reduces the substitution rates (main form of error) significantly. After having obtained the consensus sequences from each sample, we used the VSEARCH (v2.3.4) pipeline (all these steps are documented in https://github.com/torognes/vsearch/wiki/VSEARCH-pipeline) for OTU construction. The approach is as follows: we pool the reads from different samples together and add barcodes to keep an account of the samples these reads originate from. We then dereplicate the reads and sort them by decreasing abundance and discard singletons. In the next step, the reads are clustered based on 97% similarity, followed by removing clusters that have chimeric models built from more abundant reads (–uchime_denovo option in vsearch). A few chimeras may be missed, especially if they have parents that are absent from the reads or are present with very low abundance. Therefore, in the next step, we use a reference-based chimera filtering step (–uchime_ref option in vsearch) using a gold database (https://www.mothur.org/w/images/f/f1/Silva.gold.bacteria.zip). The original barcoded reads were matched against clean OTUs with 97% similarity (a proxy for species level separation) to generate OTU table (a total of 18,588 unique sequences) for n = 379 samples.

The representative OTUs were then taxonomically classified against the SILVA SSU Ref NR database release v123 database with assign_taxonomy.py script from the Qiime (Caporaso et al., 2010) workflow. To find the phylogenetic distances between OTUs, we first multisequence aligned the OTUs against each other using Kalign v2.0.4 (Lassmann and Sonnhammer, 2005) (using the options -gpo 11 -gpe 0.85) and then used FastTree v2.1.7 (Price et al., 2010) to generate the phylogenetic tree in NEWICK format. Finally make_otu_table.py from Qiime workflow was employed to combine abundance table with taxonomy information to generate biome file for OTUs. Tax4Fun (Aßhauer et al., 2015) was used to predict the functional capabilities of microbial communities based on 16S rRNA datasets (all prokaryotic KEGG organisms are available in Tax4Fun for SILVA v123 and KEGG database release 64.0) and then utilizing ultrafast protein classification (UProC) tool (Meinicke, 2015) to generate metabolic functional profiles after normalizing the data for 16S rRNA gene copy numbers. In Tax4Fun, we used MoP-Pro approach (Asshauer and Meinicke, 2013) to give pre-computed 274 KEGG Pathway reference profiles. Although Tax4Fun based metabolic prediction is constrained by the taxa available in the reference database, it gives a statistic called fraction-of-taxonomic-units-unexplained (FTU) which reflects the amount of sequences assigned to a taxonomic unit and not transferable to KEGG reference organisms. This can be used as a measure of confidence in trusting the predictions. Summary statistics of FTUs returned in this study are as follows: 1st Quantile:0.09129; Median:0.13995; Mean:0.14902; and 3rd Quantile:0.19800 (Figure 1H). Thus, on average metabolic profiles of ~86% of the taxa were present and therefore with this high representation, we used the pathways in the statistical analysis.

Figure 1

Day-wise statistical measures calculated on the microbiome data. (A) Shannon entropy with first appearance of Campylobacter (≥5 sequences) highlighted as triangles. (B–D) Local contribution to beta diversity (LCBD) calculated by using Hellinger transform on the microbial counts, Unweighted Unifrac dissimilarity (phylogenetic distances only), and Weighted Unifrac dissimilarity (phylogenetic distances weighted with abundance counts) respectively (E,F) Nearest-Taxon-Index (NTI) and nearest-relative-index (NRI) considering presence/absence of OTUs in samples (G) Richness calculated as exponentiation of Shannon entropy on the proportional representation of KEGG pathways on samples, and (H) fraction-of-taxonomic-units-unexplained (FTU) calculated on each sample. In all subfigures, the mean value is represented by solid blue line with 95% confidence interval of standard deviation given as dark shaded region around the mean. The samples are colored with respect to the pens they originate from. Based on the analysis given in this study, we have identified days 12–20 of importance and are thus highlighted as lighter shaded regions.

Statistical analysis

Statistical analyses were performed in R using the tables and data generated as above as well as the meta data associated with the study. For community analysis (including alpha and beta diversity analyses) we used the vegan package (Oksanen et al., 2015). For alpha diversity measures, we calculated: Richness, estimated number of species/features per sample; and Shannon entropy: a commonly used index to measure the balance of a community within a sample. Exponentiating Shannon entropy gives the richness profile. These alpha diversity measures are calculated after rarefying the abundance table to minimum library size, as is the norm. To calculate Unifrac distances (that account for phylogenetic closeness), we used the phyloseq (McMurdie and Holmes, 2013) package. Nonmetric Distance Scaling (NMDS) plot of community data (OTUs) used different distance measures (Vegan's metamds() function): Bray-Curtis, considers the species abundance count; Unweighted Unifrac, considers the phylogenetic distance between the branch lengths of OTUs observed in different samples without taking into account the abundances; and Weighted Unifrac, unweighted unifrac distance weighted by the abundances of OTUs. The samples are grouped for different treatments as well as the mean ordination value and spread of points (ellipses were drawn using Vegan's ordiellipse() function that represent the 95% confidence interval of the standard errors of the groups).

To understand multivariate homogeneity of groups dispersion (variances) between multiple conditions, we used Vegan's betadisper() function in which the distances between objects and group centroids are handled by reducing the original distances (BrayCurtis, Unweighted Unifrac, or Weighted Unifrac) to principal coordinates and then performing ANOVA on them. We used Vegan's adonis() for analysis of variance using distance matrices (BrayCurtis/Unweighted Unifrac/Weighted Unifrac) i.e., partitioning distance matrices among sources of variation (Grouping type i.e., weeks, body weight, feed intake, feed conversion ratio etc.). This function, henceforth referred to as PERMANOVA, fits linear models to distance matrices and uses a permutation test with pseudo-F ratios.

To find OTUs that are significantly different between multiple conditions (days/weeks), we used DESeqDataSetFromMatrix() function from DESeq2 (Love et al., 2014) package with the adjusted p-value significance cut-off of 0.05 and log2 fold change cut-off of 2. This function uses negative binomial GLM to obtain maximum likelihood estimates for OTUs log fold change between two conditions. Then Bayesian shrinkage is applied to obtain shrunken log fold changes subsequently employing the Wald test for obtaining significances. To find KEGG pathways significantly up/down-regulated between multiple conditions (days/weeks), the Kruskal-Wallis test was used with p-values adjusted for multiple comparisons using the fdrtool package (Klaus and Strimmer, 2013, 2015).

We performed Local Contribution to Beta Diversity (LCBD) analysis (Legendre and De Cáceres, 2013) by using LCBD.comp() from adespatial package (Dray et al., 2018). We used the Hellinger distance (abundances), unweighted (phylogenetic distance) and weighted Unifrac (phylogenetic distance weighted by abundance) dissimilarities. LCBD gives the sample-wise local contributions to beta diversity that could be derived as a proportion of the total beta diversity. In the context of this longitudinal study, it provides a mean to show how markedly different the microbial community structure of a single sample is from the average (with higher LCBD values representing outliers), and also provides a mean to show when the community structure has stabilized in a temporal setting.

To characterize the phylogenetic community composition within each sample whether the microbial community structure is stochastic (driven by competition among taxa) or deterministic (driven by strong environmental pressure i.e. host environment), we quantified: mean-nearest-taxon-distance (MNTD) and the nearest-taxon-index (NTI) using mntd(), and ses.mntd(); and mean-phylogenetic-diversity (MPD) and nearest-relative-index (NRI) using mpd() and ses.mpd() function from the picante (Kembel et al., 2010) package. NTI and NRI represent the negative of the output from ses.mntd() and ses.mpd(), respectively. They also quantify the number of standard deviations that the observed MNTD/MPD is from the mean of the null distribution (999 randomization by using null.model = “richness” in the ses.mntd() and ses.mpd() functions and only considering the taxa as present/absent without taking their abundances). We used the top 1,000 most abundant OTUs for calculation of these measures based on the recommendations given in (Stegen et al., 2012).

We used the “BVSTEP” routine (Clarke and Ainsworth, 1993), an algorithm that searches for highest correlation (Mantel test) between dissimilarities of a fixed and variable multivariate datasets using bvStep() from sinkr package (Taylor, 2014) by permuting through 2n-1 possible combinations of features in the variable dataset. Testing all feature combinations is unrealistic and computationally intractable when the feature space is high (18,588 OTUs in our case). Thus, we used the abundance table with 1000 most abundant OTUs (with the premise that the most abundant species that may have a significant role to play) to best correlate with the overall similarities given all the OTUs (18,588 in our case). This analysis is complimentary to the differential analysis and identified the OTUs that were causing the major shifts in beta diversity.

The phylogenetic tree and annotations summarizing the findings of this study were drawn using Evolview (http://www.evolgenius.info/evolview/).

We considered analyses on two different groupings of the sample data, comparison of microbial profiles on a daily basis to reveal temporal patterns, and on a weekly basis (4 weeks), primarily because the poultry growth and performance parameters were recorded on a weekly basis. The statistical scripts and workflows for all above can be found at http://userweb.eng.gla.ac.uk/umer.ijaz#bioinformatics.

Results

Daily diversity patterns converge to a stable community as we go forward in time

Although alpha diversity (Shannon) on microbial counts (Figure 1A) shows a rapid increase over the first ten days, it follows a plateauing effect where the microbiome normalizes at approximately day 12. This is in line with previous reports whereby the gastrointestinal (GI) tract of poultry comes into contact with exogenous microorganisms immediately after hatch and as the host grows, this microbiome becomes highly diverse until it reaches a relatively stable yet dynamic state (Pan and Yu, 2014). The same temporal phenomenon can be observed when considering local contributions to beta diversity based on abundance count (Hellinger distance; Figure 1B). When considering phylogenetic distances only (Unweighted Unifrac; Figure 1C), although the decrease in beta diversity contributions is marginally slower than the abundance counts counterpart, there is a sudden increase around day 20. Using both abundances and phylogenetic distances this seems to disappear (Weighted Unifrac; Figure 1D). It should be noted that a higher LCBD value suggests the diversity patterns of a sample is markedly different from the rest of the samples in an average sense. In contrast, the level of microbial diversity between the different pens was relatively stable (results not significant and thus not shown) suggesting less or no variability amongst pens. Campylobacter was detected in three chickens from the 12 pens at day 16 (Figure 1A). This is in line with previous reports where natural colonization of chickens has been reported at approximately day 14 of the chicken life cycle (Neill et al., 1984; Hermans et al., 2011; Kalupahana et al., 2013; Thibodeau et al., 2015). Campylobacter was also identified in one of the chickens at day 3 and previously it has also been reported that chickens between 0 and 3 days of age can become infected with Campylobacter (Cawthraw et al., 1996).

Window of opportunity for Campylobacter between day 12 and day 20

Next, we explored ecological drivers of microbial community to determine whether there is any environmental pressure (host environment) responsible for assemblage of microbial community or if it is driven purely by competition. Using NTI and NRI (Figures 1E,F), one can observe a step function response around day 12. For a single community, NTI/NRI greater than +2 indicates strong phylogenetic clustering (driven by environmental filtering) and less than −2 indicates phylogenetic overdispersion (environment has little or no role to play). Since chicken ceca are already a constrained environment to begin with (as opposed to real environmental datasets), the lower bound of −2 may not be feasible and hence the values should be taken relatively with an increasing value implying increasing host environmental pressure. It should be noted that NRI reflects the phylogenetic clustering in a broad sense (whole phylogenetic tree) with the negative values representing evenly spread community. On the other hand, NTI focuses more on the tips of the tree with positive values of NTI indicating that species co-occur with more closely related species than expected, and negative values indicating that closely related species do not co-occur. We have chosen presence/absence of species while calculating these measures without taking into account the abundances as they mask the phenomenon similar to LCBD profiles (Figures 1C,D). When we consider differential analysis of OTUs (Supplementary Table 1), we can notice that between days 9 and 11 there is a high proportion of OTUs that were log2 fold different. After day 20, we also observe the same between days 26 and 28 with the changes in phylogenetic structure responsible for peaks in NTI/NRI. Interestingly, chickens were raised on three phase diets; starter diets (days 0–10), grower diets (days 11–25) and finisher diets (days 26–35). The high proportion of OTUs that were log2 fold different between days 26 and 28 may be attributed to the change in feed from grower to finisher feed. Since the NTI/NRI are already significantly higher than 2, we do not consider this as an upper bound and revert back to day 20 as an upper bound for the window. Based on beta dispersion analysis (Table 1), we observe days 11–13 and days 19–21 when the dispersions of the microbial communities are changing significantly. The alteration in the chicken feed from starter diet (days 0–10) to grower diets (days 11–25) may also play a role in the significant beta dispersion between days 11 and 13, although the feed change does not seem a likely explanation for days 19–21. For completeness we also generated differential analysis of genus level where Campylobacter was identified as being significantly down-regulated between day 16 and day 17 (Supplementary Table 2).

Table 1

Statistics for beta dispersion comparison on daily microbiome data.

Day comparisons	Bray-curtis	Unweighted unifrac	Weighted unifrac
3–4	p = 0.018071 (*)	p = 0.18436	p = 0.085112
4–5	p = 0.85255	p = 0.18547	p = 0.25546
5–6	p = 0.60961	p = 0.1225	p = 0.73468
6–7	p = 0.82972	p = 0.94104	p = 0.21369
7–8	p = 0.71257	p = 0.88392	p = 0.47401
8–9	p = 0.060007	p = 0.94453	p = 0.36231
9–10	p = 0.9966	p = 0.11357	p = 0.53314
10–11	p = 0.20247	p = 0.20845	p = 0.13289
11–12	p = 0.38794	p = 0.014818 (*)	p = 0.62198
12–13	p = 0.88847	p = 0.064143	p = 0.013623 (*)
13–14	p = 0.63766	p = 0.16696	p = 0.41304
14–15	p = 0.9467	p = 0.64383	p = 0.46855
15–16	p = 0.89972	p = 0.055618	p = 0.79989
16–17	p = 0.59807	p = 0.37379	p = 0.41167
17–18	p = 0.70773	p = 0.66013	p = 0.30413
18–19	p = 0.40112	p = 0.92525	p = 0.5994
19–20	p = 0.020548 (*)	p = 0.087076	p = 0.56858
20–21	p = 0.033097 (*)	p = 0.12251	p = 0.52086
21–22	p = 0.29506	p = 0.055585	p = 0.90226
22–23	p = 0.24688	p = 0.90221	p = 0.99695
23–24	p = 0.79886	p = 0.71275	p = 0.34913
24–25	p = 0.21019	p = 0.67687	p = 0.11096
25–26	p = 0.14334	p = 0.20716	p = 0.97116
26–27	p = 0.96286	p = 0.044866 (*)	p = 0.80425
27–28	p = 0.50377	p = 0.096107	p = 0.1382
28–29	p = 0.91052	p = 0.87339	p = 0.69398
29–30	p = 0.34265	p = 0.60245	p = 0.11773
30–31	p = 0.61843	p = 0.55324	p = 0.20403
31–32	p = 0.24674	p = 0.082761	p = 0.50328
32–33	p = 0.73392	p = 0.53114	p = 0.62586
33–34	p = 0.7431	p = 0.36694	p = 0.57642
34–35	p = 0.16111	p = 0.20181	p = 0.77382

Asterisks denote a statistically significant difference (.

If we consider the richness of metabolic pathways (Figure 1G), we notice that they achieve stability before the microbial community at around day 6 with no obvious patterns to suggest anything apparent between day 12 and day 20 other than a marginal decrease to day 16 and increasing again onwards. However, if we consider the differential expression analysis of pathways (Supplementary Table 3), we can notice a large proportion of these pathways changing between day 14 and 15, a day before Campylobacter was first observed. We identified a reduction in lysine degradation (ko00310) from day 14 to day 15, and an increase in D-Alanine metabolism (ko00473) from day 14 to day 15. C. jejuni typically cannot utilize sugars as a carbon source as it lacks the glycolytic enzyme phosphofructokinase and so depends on the availability of free amino and keto acids scavenged from the host or from the intestinal microbiome (Parkhill et al., 2000; Velayudhan and Kelly, 2002; Lee and Newell, 2006). C. jejuni utilizes serine, aspartate, glutamate and proline preferentially as nutritional substrates in vitro with serine catabolism required for colonization of the intestinal tract (Elharrif and Mégraud, 1986; Leach et al., 1997; Hendrixson and DiRita, 2004; Velayudhan et al., 2004). Amino acids can also potentially be deaminated to a small number of intermediates that can directly feed into the central metabolism, including pyruvate (from serine and alanine), oxaloacetate (from aspartate), and 2-oxoglutarate (from glutamate) (Velayudhan et al., 2004). The variation of such metabolic pathways may give an indication as to the appearance of Campylobacter at this time point. We also identified a reduction from days 14 to day 15 of a number of pathways relating to specific bacteria; Vibrio cholerae pathogenic cycle (ko05111; Biofilm formation - Vibrio cholerae), Escherichia coli (ko05130; Pathogenic Escherichia coli infection), Salmonella species (ko05132; Salmonella infection). In addition, we identified a reduction from day 14 to 15 of Bacterial secretion systems (ko03070). Future studies are needed to elucidate and confirm the predicted pathways. In view of these findings, Camplyobacter appears at day 16 within this window of opportunity (Figure 1) where there exists a shift from competitive to environmental drivers of microbial community, with day 16 lying immediately after the most substantial changes in metabolic profiles observed over the whole period.

Analysis of dominant bacterial taxa over time

Analysis of the 50 most abundant genera (Supplementary Figure 2) have identified trends that were reported previously in the literature i.e., chicken microbiome contains Enterobacteriaceae at early days of development, and that Firmicutes increase in abundance and taxonomic diversity over time (Danzeisen et al., 2011; Ballou et al., 2016; Mancabelli et al., 2016). Escherichia.Shigella (Phylum Proteobacteria; Family Enterobacteriaceae) was identified as being highly abundant at day 3 and showed a general reduction up to approximately day 7. Escherichia.Shigella was also noted to be present after day 28. This pattern was observed for Eisenbergiella (Phylum Firmicutes; Family Lachnospiraceae) which displayed a decrease from early time points, but remained present throughout. This pattern was also observed for Ruminiclostridium (Phylum Firmicutes; Family Ruminococcaceae) which however was not in the abundant genera after day 23. Flavonifractor (Phylum Firmicutes; Family -) was identified consistently at early time points, but was rarely abundant after day 19. Enterobacter (Phylum Proteobacteria; Family Enterobacteriaceae) was only observed at days 3 and 4 and was not abundant at any other time points. Here we identified that Ruminiclostridium.5 and Ruminiclostridium.9 (Phylum Firmicutes; Family Ruminococcaceae) which were consistently present throughout at a relatively significant level of abundance. This was also the case for Anaerotruncus (Phylum Firmicutes; Family Clostridiaceae), but at a lower level of abundance, especially before day 7. Faecalibacterium (Phylum Firmicutes; Family Clostridiaceae) was rarely abundant at early time points, however was observed consistently at a relative high abundance after day 14. Lachnoclostridium (Phylum Firmicutes; Family Lachnospiraceae) was found to be present throughout with a relatively high level of fluctuation. Certain genera such as Ruminococcaceae.UCG.005 and Ruminococcaceae.UCG.014 (Phylum Firmicutes; Family Ruminococcaceae) were not abundant at high levels at early time points however increased significantly at approximately days 16-19. Finally, Megamonas (Phylum Firmicutes; Family Veillonellaceae) and Intestinimonas (Phylum Firmicutes; Family -) were not abundant throughout most time points, before appearing post day 22-25 onwards.

Weekly microbial profiles and analysis of poultry performance metadata

The metadata collected here included Bird Weight (BW_Mean; grams), Body Weight Gain (Gain; g/bird), Feed Intake (FI), Feed Conversion Ratio (FCR), and was recorded on a weekly basis where we have considered grouping the microbiome samples accordingly; days 03–07 (week 1), days 08–14 (week 2), days 15–24 (week 3), and days 25–35 (week 4). As is the case with the daily microbiome profile, alpha diversity (rarefied richness and Shannon; Figure 2A) increases over time, however, due to the nature of this grouping, we lose the plateauing effect over time. In accordance with daily analysis, we can see a major shift in the parameters as we transition from days 08–14 to days 15–24 (Figure 2B). FCR in particular increases substantially in this period remaining stable for week 4 (days 25–35). Gain is also significantly elevated in this transition period (days 08–14 to days 15–24) when compared to other periods. In terms of beta diversity (Figure 2C), we observe the samples more sparsely spread in the first week (days 03–07) as compared to other weeks on abundance (Bray-Curtis) alone. The phylogenetic dispersion (Unweighted Unifrac) on the other hand is more preserved. We can also notice a gradient forming with later weeks more or less close to suggest convergence as we established in the case of daily profiles. Based on beta dispersion analysis (Table 2), we can notice that the dispersion in week 1 is significantly different to other weeks with 16, 6, and 17% variability in microbial community explained by PERMANOVA using counts alone (Bray-Curtis), phylogenetic distance alone (Unweighted Unifrac), and combination of the two (Weighted Unifrac), respectively. With this grouping, main sources of variation are then the distribution of species rather than their phylogenetic relatedness. The metadata explains 10–12% variability (all significant) in terms of counts alone (Bray-Curtis) with 3–6% in terms of phylogeny (Unweighted Unifrac). For the sake of completeness, we also performed differential analysis of OTUs and pathways on a consecutive weekly basis (lower halves of Supplementary Tables 1–3); however, these should be interpreted with great care as main source of variability are the daily changes and grouping samples on weekly basis will always return more significant OTUs and pathways.

Figure 2

Week-wise measures calculated on the microbiome data (A) Alpha diversity measures: richness (after rarefying the samples to minimum library size) and Shannon entropy (B) Extrinsic parameters calculated on weekly basis were mean body weight (BW_mean), body weight gain (Gain), feed intake (FI), feed conversion ratio (FCR), and (C) Beta diversity measures using Bray-Curtis (counts), Unweighted Unifrac (phylogenetic distance), and Weighted Unifrac (phylogenetic distance weighted by abundance counts). In (A,B) we have performed pair-wise ANOVA and where significant the pairs were connected with p-values drawn on top. In (C) the ellipses represent the 95% confidence interval of the standard error of the ordination points of a given grouping with labels drawn at the center (mean) of the ordination points.

Table 2

Statistics for pairwise beta dispersion and PERMANOVA when using different dissimilarity measures on weekly microbiome data.

Beta dispersion		Bray-curtis	Unweighted unifrac	Weighted unifrac
Day 03–07	Day08–14	p = 0.0061142 (**)	p = 0.00014712 (***)	p = 9.6914e−05 (***)
	Day15–24	n.s.	p = 0.010418 (*)	p = 2.5203e−09 (***)
	Day25–35	p = 0.042066 (*)	p = 0.00015112 (***)	p = 3.5789e−12 (***)
Day08–14	Day15–24	p = 0.00077017 (***)	n.s.	p = 0.019953 (*)
	Day25–35	n.s.	n.s.	p = 0.0011717 (**)
Day15–24	Day25–35	p = 0.0075651 (**)	p = 0.020128 *	n.s.
PERMANOVA
Groups		R2 = 0.16763 (p = 0.001) (***)	R2 = 0.06048 (p = 0.001) (***)	R2 = 0.17577 (p = 0.001) (***)
BW_Mean		R2 = 0.11721 (p = 0.001) (***)	R2 = 0.03964 (p = 0.001) (***)	R2 = 0.08723 (p = 0.001) (***)
FI		R2 = 0.11856 (p = 0.001) (***)	R2 = 0.04069 (p = 0.001) (***)	R2 = 0.09301 (p = 0.001) (***)
FCR		R2 = 0.1086 (p = 0.001) (***)	R2 = 0.03842 (p = 0.001) (***)	R2 = 0.11787 (p = 0.001) (***)
Gain		R2 = 0.11886 (p = 0.001) (***)	R2 = 0.04146 (p = 0.001) (***)	R2 = 0.0998 (p = 0.001) (***)

Asterisks denote a statistically significant difference (

p < 0.05,

p < 0.01,

p < 0.001).

In beta dispersion analysis, the pair-wise differences in distances from group center/mean were subjected to ANOVA after performing Principle Coordinate Analysis, and if significant (p ≤ 0.05) the values are shown. In PERMANOVA analysis, R.

Key species representing majority of the shift in community dynamics

In addition to differential analysis on OTUs (Supplementary Table 1) which returned OTUs that were log2 fold different between consecutive days, we also considered the subset analysis where we imploded the abundance table to the minimum set of OTUs, the resulting reduced-order abundance table correlated highly with the full table by preserving the beta diversity between the samples (Table 3). To see how much variability is lost, the PERMANOVA with full OTU table (18,588 OTUs) is provided as a reference. The 17 OTUs listed represent only ~2% (Subset S1 in Table 3) loss in variability and thus represent the main OTUs that are driving the community dynamics. In terms of metadata, the loss in variability is ~1% (Subset S1 in Table 3). The subset of the phylogenetic tree of these OTUs, in addition to those selected in the differential analysis (daily comparisons), a total of 110 OTUs were then extracted and annotated in Figure 3 along with taxonomy information. It can be seen that majority of these (>50%) belong to Firmicutes (Bacillaceae, Ruminococcaceae, Lachnospiracaeae, Lactobacillaceae, Peptostreptococcaceae, and Clostridiales vadin BB60 group), with a small proportion belonging to Actinobacteria (Coriobacteriacaea), Tenericutes (Mollicutes RF9), and Proteobacteria (Enterobacteriaceae including Escherichia.Shigella as mentioned before).

Table 3

Subset analysis from BVSTEP routine listing top 18 subsets with highest correlation with the full OTU table considering Bray-Curtis distance done on weekly basis.

Subsets of top 1000 most abundant OTUs		Correlation with full OTU table (R)	PERMANOVA (full OTU table)
			Groups	BW_Mean	FI	FCR	Gain
			R2 = 0.16763 (p = 0.001) (***)	R2 = 0.11721 (p = 0.001) (***)	R2 = 0.11856 (p = 0.001) (***)	R2 = 0.1086 (p = 0.001) (***)	R2 = 0.11886 (p = 0.001) (***)
			PERMANOVA (subsets)
			Groups	BW_Mean	FI	FCR	Gain
S1	OTU_2165 + OTU_2448 + OTU_33 + OTU_1121 + OTU_23 + OTU_2474 + OTU_6 + OTU_28 + OTU_157 + OTU_15 + OTU_24 + OTU_3028 + OTU_2496 + OTU_1024 + OTU_10 + OTU_3 + OTU_2555	0.833	R2 = 0.14768 (p = 0.001) (***)	R2 = 0.10732 (p = 0.001) (***)	R2 = 0.10784 (p = 0.001) (***)	R2 = 0.10117 (p = 0.001) (***)	R2 = 0.11143 (p = 0.001) (***)
S2	OTU_2165 + OTU_2448 + OTU_33 + OTU_1121 + OTU_23 + OTU_2474 + OTU_6 + OTU_28 + OTU_157 + OTU_15 + OTU_24 + OTU_3028 + OTU_2496 + OTU_1024 + OTU_3 + OTU_2555	0.83	R2 = 0.13977 (p = 0.001) (***)	R2 = 0.10028 (p = 0.001) (***)	R2 = 0.09987 (p = 0.001) (***)	R2 = 0.09203 (p = 0.001) (***)	R2 = 0.1025 (p = 0.001) (***)
S3	OTU_2165 + OTU_2448 + OTU_33 + OTU_1121 + OTU_23 + OTU_2474 + OTU_6 + OTU_28 + OTU_157 + OTU_15 + OTU_24 + OTU_3028 + OTU_2496 + OTU_1024 + OTU_3	0.827	R2 = 0.14186 (p = 0.001) (***)	R2 = 0.10205 (p = 0.001) (***)	R2 = 0.10165 (p = 0.001) (***)	R2 = 0.09344 (p = 0.001) (***)	R2 = 0.10435 (p = 0.001) (***)
S4	OTU_2165 + OTU_2448 + OTU_33 + OTU_1121 + OTU_23 + OTU_2474 + OTU_6 + OTU_28 + OTU_157 + OTU_15 + OTU_24 + OTU_2496 + OTU_1024 + OTU_3	0.823	R2 = 0.14241 (p = 0.001) (***)	R2 = 0.10262 (p = 0.001) (***)	R2 = 0.10228 (p = 0.001) (***)	R2 = 0.0939 (p = 0.001) (***)	R2 = 0.10509 (p = 0.001) (***)
S5	OTU_2165 + OTU_2448 + OTU_33 + OTU_1121 + OTU_23 + OTU_2474 + OTU_6 + OTU_28 + OTU_15 + OTU_24 + OTU_2496 + OTU_1024 + OTU_3	0.816	R2 = 0.14289 (p = 0.001) (***)	R2 = 0.10313 (p = 0.001) (***)	R2 = 0.10279 (p = 0.001) (***)	R2 = 0.09436 (p = 0.001) (***)	R2 = 0.1056 (p = 0.001) (***)
S6	OTU_2165 + OTU_2448 + OTU_33 + OTU_1121 + OTU_23 + OTU_2474 + OTU_6 + OTU_28 + OTU_15 + OTU_24 + OTU_2496 + OTU_1024	0.809	R2 = 0.14742 (p = 0.001) (***)	R2 = 0.10587 (p = 0.001) (***)	R2 = 0.10556 (p = 0.001) (***)	R2 = 0.098 (p = 0.001) (***)	R2 = 0.1084 (p = 0.001) (***)
S7	OTU_2165 + OTU_2448 + OTU_33 + OTU_1121 + OTU_2474 + OTU_6 + OTU_28 + OTU_15 + OTU_24 + OTU_2496 + OTU_1024	0.799	R2 = 0.17779 (p = 0.001) (***)	R2 = 0.12583 (p = 0.001) (***)	R2 = 0.12598 (p = 0.001) (***)	R2 = 0.11968 (p = 0.001) (***)	R2 = 0.12936 (p = 0.001) (***)
S8	OTU_2165 + OTU_2448 + OTU_33 + OTU_1121 + OTU_2474 + OTU_6 + OTU_28 + OTU_15 + OTU_24 + OTU_2496	0.789	R2 = 0.14605 (p = 0.001) (***)	R2 = 0.0875 (p = 0.001) (***)	R2 = 0.08646 (p = 0.001) (***)	R2 = 0.08927 (p = 0.001) (***)	R2 = 0.08866 (p = 0.001) (***)
S9	OTU_2165 + OTU_2448 + OTU_33 + OTU_1121 + OTU_2474 + OTU_28 + OTU_15 + OTU_24 + OTU_2496	0.777	R2 = 0.14132 (p = 0.001) (***)	R2 = 0.10216 (p = 0.001) (***)	R2 = 0.10313 (p = 0.001) (***)	R2 = 0.0989 (p = 0.001) (***)	R2 = 0.10736 (p = 0.001) (***)
S10	OTU_2165 + OTU_2448 + OTU_33 + OTU_1121 + OTU_2474 + OTU_28 + OTU_15 + OTU_24	0.763	R2 = 0.13922 (p = 0.001) (***)	R2 = 0.10051 (p = 0.001) (***)	R2 = 0.10121 (p = 0.001) (***)	R2 = 0.09676 (p = 0.001) (***)	R2 = 0.10532 (p = 0.001) (***)
S11	OTU_2165 + OTU_2448 + OTU_33 + OTU_1121 + OTU_2474 + OTU_28 + OTU_15	0.746	R2 = 0.13173 (p = 0.001) (***)	R2 = 0.09081 (p = 0.001) (***)	R2 = 0.09291 (p = 0.001) (***)	R2 = 0.09023 (p = 0.001) (***)	R2 = 0.09828 (p = 0.001) (***)
S12	OTU_2165 + OTU_2448 + OTU_33 + OTU_1121 + OTU_2474 + OTU_28	0.723	R2 = 0.09574 (p = 0.001) (***)	R2 = 0.04819 (p = 0.001) (***)	R2 = 0.04936 (p = 0.001) (***)	R2 = 0.05663 (p = 0.001) (***)	R2 = 0.05163 (p = 0.001) (***)
S13	OTU_2165 + OTU_2448 + OTU_33 + OTU_1121 + OTU_2474	0.696	R2 = 0.0952 (p = 0.001) (***)	R2 = 0.04875 (p = 0.001) (***)	R2 = 0.05019 (p = 0.001) (***)	R2 = 0.05606 (p = 0.001) (***)	R2 = 0.05246 (p = 0.001) (***)
S14	OTU_2165 + OTU_2448 + OTU_1121 + OTU_2474	0.661	R2 = 0.10232 (p = 0.001) (***)	R2 = 0.05483 (p = 0.001) (***)	R2 = 0.05623 (p = 0.001) (***)	R2 = 0.06108 (p = 0.001) (***)	R2 = 0.05869 (p = 0.001) (***)
S15	OTU_2448 + OTU_33 + OTU_1121 + OTU_2474	0.655	R2 = 0.06994 (p = 0.001) (***)	R2 = 0.02489 (p = 0.001) (***)	R2 = 0.02416 (p = 0.001) (***)	R2 = 0.03155 (p = 0.001) (***)	R2 = 0.02375 (p = 0.001) (***)
S16	OTU_2448 + OTU_33 + OTU_2474	0.604	R2 = 0.03489 (p = 0.001) (***)	R2 = 0.01348 (p = 0.003) (**)	R2 = 0.01238 (p = 0.005) (**)	R2 = 0.01154 (p = 0.006) (**)	R2 = 0.01053 (p = 0.012) (*)
S17	OTU_33 + OTU_1121 + OTU_2474	0.599	R2 = 0.06662 (p = 0.001) (***)	R2 = 0.01995 (p = 0.001) (***)	R2 = 0.02047 (p = 0.001) (***)	R2 = 0.03201 (p = 0.001) (***)	R2 = 0.02183 (p = 0.001) (***)
S18	OTU_1121 + OTU_2474	0.538	R2 = 0.07 (p = 0.001) (***)	R2 = 0.02571 (p = 0.001) (***)	R2 = 0.02628 (p = 0.001) (***)	R2 = 0.03796 (p = 0.001) (***)	R2 = 0.02793 (p = 0.001) (***)

Asterisks denote a statistically significant difference (

p < 0.001).

For each subset, PERMANOVA was performed against different sources of variations.

OTU_2165,Bacteria;Firmicutes;Clostridia;Clostridiales;Lachnospiraceae.

OTU_2448:Bacteria;Firmicutes;Clostridia;Clostridiales;Ruminococcaceae;Ruminiclostridium.

OTU_33:Bacteria;Firmicutes;Clostridia;Clostridiales;Ruminococcaceae;Ruminiclostridium 5.

OTU_1121:Bacteria;Firmicutes;Clostridia;Clostridiales;Lachnospiraceae;Eisenbergiella.

OTU_23:Bacteria;Firmicutes;Clostridia;Clostridiales;Ruminococcaceae;Ruminiclostridium 9.

OTU_2474:Bacteria;Firmicutes;Clostridia;Clostridiales;Ruminococcaceae;Ruminiclostridium 5.

OTU_6:Bacteria;Firmicutes;Clostridia;Clostridiales;Ruminococcaceae.

OTU_28:Bacteria;Firmicutes;Clostridia;Clostridiales;Ruminococcaceae.

OTU_157:Bacteria;Firmicutes;Bacilli;Lactobacillales;Lactobacillaceae;Lactobacillus.

OTU_15:Bacteria;Firmicutes;Clostridia;Clostridiales;Ruminococcaceae.

OTU_24:Bacteria;Firmicutes;Clostridia;Clostridiales;Ruminococcaceae;Ruminiclostridium.

OTU_3028:Bacteria;Firmicutes;Clostridia;Clostridiales;Lachnospiraceae.

OTU_2496:Bacteria;Firmicutes;Clostridia;Clostridiales;Lachnospiraceae;Tyzzerella.

OTU_1024:Bacteria;Firmicutes;Clostridia;Clostridiales;Ruminococcaceae;Faecalibacterium.

OTU_10:Bacteria;Firmicutes;Clostridia;Clostridiales;Lachnospiraceae.

OTU_3:Bacteria;Firmicutes;Clostridia;Clostridiales;Ruminococcaceae;Ruminiclostridium 5.

OTU_2555:Bacteria;Firmicutes;Clostridia;Clostridiales;Lachnospiraceae.

Figure 3

Phylogenetic tree of the subset of OTUs selected as significant on differential analysis (based on Table 3 and Supplementary Table 1). Next to the OTU labels are descriptive text representing where the OTUs were found to be significant, for example, the first entry for OTU 231, “u 26-27 d 27-28 u 30-31,” can be read as upregulated going from day 26 to 27 and then from day 30 to 31 and downregulated going from day 27 to 28. “b” represents the OTUs selected in the subset analysis. The next two columns are a pictorial representation of the above-mentioned descriptive text with pink color representing OTUs selected in subset analysis, red color for upregulated OTUs, blue for downregulated OTUs, and purple for OTUs which show the both trends (up/down regulation). The next column shows the taxonomy of the OTUs according to SILVA v123 with coloring at unique family level. The heatmap was drawn by collating the mean values of OTUs for samples from the same day after performing proportional standardization on the full OTU table using wisconsin() function.

Discussion

Comprehensive investigation of the chicken cecal microbiome at a day to day level revealed a rapid increase in diversity up to day 12, with microbial variation observed both in terms of genera and abundance. We suspect this early variation is due to competitive factors determined by space and available food resources. Post day 20 there exists a considerable stabilization of the chicken cecal microbiome where the relative microbial diversity and abundances are standardized, with environmental factors (in this case the host chicken) exerting a greater influence on any change in the microbial diversity. Between days 12 and 20 we observe a shift from competitive to environmental drivers of microbial community creating a window of opportunity whereby Campylobacter appears. We identified Campylobacter at day 16 with this day lying immediately after the most substantial changes in metabolic profiles observed over the whole period. Whilst we identified Campylobacter within 25% of the pens on day 16, we would naturally expect Campylobacter to spread to other chickens and pens and also be identified on subsequent days. We suspect that the experimental set-up here was such that following random selection of birds from each pen on each day, sacrificing the bird (to perform gDNA extraction from the ceca) did not allow for an opportunity for Campylobacter to spread to other chickens or pens. Clearly in a typical farm set-up this would not be the case and Campylobacter would spread naturally.

Microbial variation over time is most likely influenced by diet of the chickens whereby significant shifts in OTU abundances and beta dispersion of the samples often corresponded with changes in feed. Notably, the relatively high proportion of OTUs that were log2 fold different between days 9 and 11, and days 26 and 28, and beta dispersion for days 11–13 corresponded with changes in feed from grower to finisher. Further studies investigating different feed content is required to ascertain the complete impact on chicken cecal microbiome.

Previous microbiome studies of chicken ceca have often lacked the day to day sampling points, replicate numbers, relevant metadata and have often provided external Campylobacter infection that may potentially perturb the natural habitat. These have not allowed us to fully appreciate what is occurring in a natural environment in relation to how and when Campylobacter appears within the chicken gut. Thus, we believe the major strength of this study is that we have filled these gaps by performed the most comprehensive analysis of the chicken cecal microbiome to date. This was made possible by sampling from days 3 to 35, with 12 replicates per day (final n = 379), correlating additional metadata such as chicken weight and feed conversion rates and with Campylobacter detection in a natural environmental setting giving the most comparable experimental design to a farm set-up. As we were not able to sample the same chicken for all time points, future studies should investigate this further with added dietary information than what we have considered here, with experimental designs also to investigate and confirm the predicted pathways.

Conclusions

Industry has endeavored to reduce the burden of Campylobacter within chicken production lines with supplements often administered with the aim of performance enhancing and/or reducing bacteria such as Campylobacter, typically post day 25. The relative stability of the chicken cecal microbiota at this time point may explain the efficacy of such products, however the identification of a window of opportunity for bacteria such as Campylobacter may call for intervention strategies between days 12 and 20, or even earlier. This study can act as a baseline for future intervention strategies and help reduce the burden of Campylobacter within chickens.

Availability of supporting data

The raw sequence files supporting the results of this article are available in the European Nucleotide Archive under the project accession number PRJEB25776.

Author contributions

AM, AR, UL, BW, ND, NC, and OG contributed to the study design. ND, NC, and OG managed the study. CK, AM, AS, and ML performed the sample collection and DNA extraction. LS, AM, AE, and OG performed the library preparation and Illumina MiSeq sequencing at the LSHTM. UZI wrote the analysis scripts to generate the figures and tables in this paper. UZI and OG performed the bioinformatics and statistical analysis. UZI, ND, NC, and OG drafted the initial version of the manuscript with all authors contributed to redrafting.

Conflict of interest statement

AM, AR and UL are employed by company Moy Park. AM is also a PhD student between Moy Park and AFBI. The remaining authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.