Development of the Caecal Microbiota in Three Broiler Breeds

The development of the caecal microbiota plays a role in the metabolism and immune competence of chickens. A detailed understanding of normal succession in the caecal microbiota can inform the use of probiotics and other interventions to optimize the caecal microbiota. The development of the microbiota in caecal mucus and lumen samples from three breeds of broiler chicken (Cobb 500, n = 36; Hubbard JA87, n = 38; and Ross 308, n = 36) was observed between 0 and 42 days post hatch. Chicks were housed in the same room of a climate-controlled, biosecure chicken housing unit. Between 0 and 14 days post hatch, chicks were kept in brooder pens ensuring a mixture of breeds in each brooder. From 22 days post hatch, chicks were removed from the brooders and kept in the same room. DNA was extracted from a pooled sample of caecal mucus and luminal contents from five birds of each breed at 0, 3, 7, 14, 21, 28, and 42 days post hatch. High-throughput Illumina sequencing was performed for the V4 hypervariable region of the 16S rRNA gene. The early caecal microbiota was characterized by poor diversity and dominance by one or two bacterial species. Early colonizers of the caecum included Bifidobacteriaceae, Lachnospiraceae, Bacteroidaceae and Burkholderiaceae with some amplicon sequence variants (ASVs) assigned to Ruminococcaceae. Later colonizers of the caecal microbiota were most apparent from 14 d.p.h and included Ruminococcaceae, Clostridiales vadin BB60 group, Christensenellaceae and Bacillaceae. The caecal microbiota continued to change until 42 d.p.h when the microbiota was characterized by a high abundance of Bacteroidaceae, Lachnospiraceae and Ruminococcaceae. The lumen microbiota was significantly different to the mucus with some ASVs assigned to Lachnospiraceae, Ruminococcaceae, Christensenellaceae and Bacillaceae showing increased abundance in the mucus. ASVs assigned to Bacteroidaceae, Lactobacillaceae and Burkholderiaceae showed a preference for the lumen. Analysis of five caecal mucus samples from each breed at 42 days post hatch showed differences in microbiota composition between Ross and Cobb as well as between Ross and Hubbard. Since performance data was not collected no functional inferences as to the significance of this finding can be made.

The development of the caecal microbiota plays a role in the metabolism and immune competence of chickens. A detailed understanding of normal succession in the caecal microbiota can inform the use of probiotics and other interventions to optimize the caecal microbiota. The development of the microbiota in caecal mucus and lumen samples from three breeds of broiler chicken (Cobb 500, n = 36; Hubbard JA87, n = 38; and Ross 308, n = 36) was observed between 0 and 42 days post hatch. Chicks were housed in the same room of a climate-controlled, biosecure chicken housing unit. Between 0 and 14 days post hatch, chicks were kept in brooder pens ensuring a mixture of breeds in each brooder. From 22 days post hatch, chicks were removed from the brooders and kept in the same room. DNA was extracted from a pooled sample of caecal mucus and luminal contents from five birds of each breed at 0, 3, 7, 14, 21, 28, and 42 days post hatch. High-throughput Illumina sequencing was performed for the V4 hypervariable region of the 16S rRNA gene. The early caecal microbiota was characterized by poor diversity and dominance by one or two bacterial species. Early colonizers of the caecum included Bifidobacteriaceae, Lachnospiraceae, Bacteroidaceae and Burkholderiaceae with some amplicon sequence variants (ASVs) assigned to Ruminococcaceae. Later colonizers of the caecal microbiota were most apparent from 14 d.p.h and included Ruminococcaceae, Clostridiales vadin BB60 group, Christensenellaceae and Bacillaceae. The caecal microbiota continued to change until 42 d.p.h when the microbiota was characterized by a high abundance of Bacteroidaceae, Lachnospiraceae and Ruminococcaceae. The lumen microbiota was significantly different to the mucus with some ASVs assigned to Lachnospiraceae, Ruminococcaceae, Christensenellaceae and Bacillaceae showing increased abundance in the mucus. ASVs assigned to Bacteroidaceae, Lactobacillaceae and Burkholderiaceae showed a preference for the lumen. Analysis of five caecal mucus samples from each breed at 42 days post hatch showed differences in microbiota composition between Ross and Cobb as well as between Ross and Hubbard. Since performance data was not collected no functional inferences as to the significance of this finding can be made.

INTRODUCTION
The intestinal microbiota of an individual chicken may be composed of between 200 and 350 different bacterial species (1) while around 640 bacterial species have so far been identified in the chicken gastrointestinal tract (2). In recent years it has become apparent that this diverse range of bacteria are not innocuous bystanders but play a range of roles in the host, from metabolism to immune maturation (3). With this realization has come the desire to identify beneficial bacteria and modulate their abundance to accentuate their effects. However, before successful interventions can be implemented a more detailed understanding of the normal development of the intestinal microbiota is required. The caeca are often chosen as a site for microbiota investigation due to their importance in both metabolism and immune maturation. Digestive caecal processes, such as bacterial fermentation to produce short chain fatty acids, provide up to 10% of a chicken's metabolizable energy (4). In terms of importance in immunity, the caecum is an important site of colonization by pathogens such as Campylobacter jejuni as well as receiving constant bacterial challenge by environmental bacteria introduced by reflux from the urodeum and cloaca (5,6).
Several studies have investigated the microbiota of day-old chicks, but when considering day-old chicks, there may be a problem with terminology. Day-old chicks sold by hatcheries may not be less than 24 h old and can be over 72 h old at point of sale. Microbes inhabiting the gut are derived from the environment during incubation, hatching and handling during delivery. This initial colonization is variable, and the composition of the initial intestinal microbiota will differ significantly between chicks from different hatcheries (7). This likely explains the wide variety of results obtained by different groups when examining the intestinal microbiota in day-old chicks. For example, Ballou et al. (8) identified Gammaproteobacteria as contributing roughly 85% of sequences recovered from caecal samples from day old chicks while Pedroso et al. (9) found that Pelotomaculum, a genus of Clostridiales, and Enterococcus were the most abundant genera at the same age.
Once chicks have arrived on the farm they are exposed to a more diverse microbial environment. Bacteria are ingested from litter, feed and water. Succession occurs rapidly and differentiation of populations in different gut compartments occurs at a young age. The caecal and ileal microbiota begin to diverge from about 3 days old (10, 11). Ballou et al. (8) describes the day-old chick microbiota as dominated by Gammaproteobacteria, mainly Enterobacteriaceae, with a smaller population of Enterococcus. By 3 days old, sequences belonging to Ruminococcaceae and other Firmicutes commonly found in the caecum, such as Clostridiales, are present. These were the dominant bacteria by 14 and 28 days old (8). Oakley et al. (12) found a caecal microbiota dominated by Clostridiales such as Flavonifractor, Pseudoflavonifractor and Lachnospiraceae by 7 days old. The microbiota remained dominated by Clostridiales, however the dominant genus switched to Faecalibacterium and Roseburia at 21 and 42 days old, respectively. Most importantly, the diversity of the caecal microbiota continued to develop, with a doubling of the number of genera between 7 and 42 days old (12). Wise and Siragusa (13) obtained similar results showing a shift from Enterobacteriaceae to Clostridiales by 14 days old. Zhu and Joerger (14), using FISH, describe a similar change in the caecal microbiota between 2 days and 6 weeks old with an increase in diversity. Initially, the microbial community was composed of Enterobacteriaceae, Lactobacillus and Bifidobacterium. These groups give way over time to produce a microbiota mainly composed of Clostridiales (14).
There is also some evidence that the intestinal microbiota can be influenced by host genotype in chickens (15)(16)(17). These studies have used divergent or inbred lines rather than commercial breeds making the results of limited direct use to the poultry industry. Another potential source of variation is that the experimental groups are often housed separately. Since environmental factors such as litter type and diet (18,19) are known to affect the composition of the microbiota it is reasonable to question whether housing experimental groups separately introduces a confounding variable.
This study aims to revisit the topic of normal caecal microbiota development using the increased resolution of next generation sequencing to shed light on microbial succession. The second objective of this study was to observe the development of the caecal microbiota in three common breeds of broiler chicken (Cobb 500, Hubbard JA87, and Ross 308) whilst they are housed together.

Animals and Housing
One hundred and ten (36 Cobb 500, 38 Hubbard JA87 and 36 Ross 308) day-old chicks were obtained from a single commercial hatchery. Chicks were distributed across three circular brooder pens (2 metres diameter) in the same room of a climatecontrolled, biosecure chicken housing unit. Each brooder used a wood shaving substrate and contained the same number of chicks from each breed. Chicks were tagged with colored wing tags to allow accurate identification of the different breeds but were not individually identifiable. Water and feed were provided ad libitum by a drinker and feeder in each brooder. Chicks were fed a pelleted vegetable protein-based starter diet (Special Diet Services, Witham, Essex, UK) until 14 days post hatch (d.p.h). From 14 d.p.h a pelleted vegetable proteinbased grower diet (Special Diet Services, Witham, Essex, UK) was provided until the end of the experiment. Nutritional composition of the starter and grower diets is displayed in Table 1 with a full list of ingredients and additives provided in the Supplementary Table 1. No coccidiostats or antimicrobials were added to either diet due to the high biosecurity levels maintained in the housing. At 22 d.p.h the birds no longer required brooder lamps and, as such, they were removed from the brooders and housed together in the same room on wood shavings. Temperature in the birds pens was maintained between 25 and 30 • C. No mortality was observed during the study. All experimental protocols were conducted in accordance with the Animals (Scientific Procedures) Act 1986 under project licence 40/3652 and was approved by the University of Liverpool

Illumina MiSeq Sequencing
Extracted DNA was sent for paired-end sequencing of the 16S rRNA gene at the Centre for Genomic Research (University of Liverpool) using an Illumina MiSeq run. The V4 hypervariable region (515F/R806) was amplified to yield an amplicon of 254 base pairs (20). Library preparation was performed using a universal tailed tag design with subsequent amplification performed using a two step PCR with a HiFi Hot Start polymerase (Kapa) (21). The first round of PCR was performed using the primers 5 ′ -ACACTCTTTCCCTACACGACGCTC TTCCGATCTNNNNNGTGCCAGCMGCCGCGGTAA-3 ′ (forward) and 5 ′ -GTGACTGGAGTTCAGACGTGTGCTCTT CCGATCTGGACTACHVGGGTWTCTAAT-3 ′ (21

Data Analysis and Statistics
Alpha and beta diversity analyses were performed at a sampling depth of 23,000 using the alignment (27), phylogeny (28) and diversity (https://github.com/qiime2/q2-diversity) plugins. Alpha diversity, a metric used to assess species richness, was measured using an observed ASVs metric and compared between samples using a Kruskal Wallis test with a false discovery rate (FDR) correction. Taxa plots were produced using the q2-taxa plugin (https://github.com/qiime2/q2-taxa). Beta diversity, a metric used to compare species diversity and abundance between samples, was calculated with a weighted UniFrac metric. The beta diversity matrix was used to draw principal coordinate analysis (PCoA) plots and an ANOSIM test was used to determine the significance of differences in beta diversity between groups. Gneiss analysis was chosen to analyse differential abundance between groups since it overcomes challenges created by the compositional nature of microbiota data. To facilitate interpretation of results, less common taxa were filtered by excluding ASVs with less than the median frequency. Next, a dendrogram of ASVs is prepared using correlation clustering. Each node in the dendrogram is treated as a "balance" with taxa on one side of the balance termed numerators and on the other, denominators. Gneiss analysis examines the log ratio of abundances between numerator and denominator taxa at each balance. Each log ratio's final numerical value is dependent on the balance between the taxa composing the numerator and those composing the denominator of the ratio. Differences in the log ratio of a balance can be compared between sample groups to determine differences in microbiota composition. A significant difference between samples allows hypotheses to be formulated regarding changes in the absolute abundance of numerator and denominator taxa but gives no further information as to which hypothesis is correct. For example, if balance y0 is found to be significantly lower at Time A compared to Time B the following hypotheses could explain the result: (i) The numerator taxa have increased between times A and B; (ii) The denominator taxa have an decreased between times A and B; (iii) A combination of hypotheses (i) and (ii); (iv) Both numerator and denominator taxa have increased between times A and B, but numerator taxa have increased more; (v) Both numerator and denominator taxa have decreased between times A and B, but denominator taxa have decreased more. Further investigations, such as quantitative PCR, are required to discern which hypothesis is correct (29).
Gneiss analysis (29) was run using the gneiss plugin (https://biocore.github.io/gneiss/) to identify taxa which were differentially abundant between time points and area sampled. Principal balances for use in Gneiss were obtained via Ward's hierarchical clustering using the correlation-clustering command. Isometric log ratios for each balance were calculated using the ilr-transform command. A multivariate response linear regression model of log ratios balances was constructed with area, breed and days post hatch as covariates using the ols-regression command. Results were visualized through a regression summary and dendrogram heatmaps. Balances significantly affected by the covariates "days post hatch, " "breed, " and "area" were identified as those with a p-value less than 0.05. The results of this analysis were used to select taxa for further analysis using quantitative PCR.

Quantitative PCR
Taxa were selected for further exploration using quantitative PCR based on results from Gneiss analysis. A literature search was conducted to find suitable primers. Where suitable primers were not available, the sequences retrieved from Illumina sequencing were used to produce taxa specific primers. The sequence was input into Primer-BLAST and a suitable primer pair was chosen.
To test specificity of primers, each primer pair was input into TestPrime for comparison against the SILVA database SSU-r132. Further testing of primers was conducted using PCR. The primers were tested against known positive and negative samples to check for the correct amplicon size and non-specific amplification. A gradient PCR was conducted to establish the correct annealing temperature for quantitative PCR. Primers used are displayed in Table 2.
The real-time quantitative PCR assay was conducted on a 1:10 solution of extracted DNA with a Rotor-Gene Q (Qiagen) and Rotor-Gene SYBR Green PCR kits (Qiagen). The V4 region of the 16S rRNA gene was used as a reference gene. Rotor-Gene Q software (version 2.3.1.49) was used to produce melting curves and identify the cycle threshold (Ct), the point at which fluoresence above the background level is detectable. Each sample was run in triplicate with an averaged Ct used in further analysis. The Ct, defined as the difference between the Ct value for taxa specific primers and the Ct value for the reference gene, was calculated for each sample. Results were expressed as 40 − Ct. Results from qPCR were compared to the relative abundance of the corresponding taxonomic group using a Spearman rank-order correlation coefficient to assess correlation between abundance determined by qPCR and sequencing.

Alpha and Beta Diversity
At hatch, the caecum was populated by an average of 35 ASVs. By 3 d.p.h, this increased to an average of 60 ASVs although this change was not significant likely due to the small sample size used. There was no significant change in alpha diversity between 3 and 7 d.p.h although the average number of ASVs increased to 79 nor was there a significant change between 21 and 28 d. From the β-diversity plots a pattern of succession becomes apparent (Figure 2). When measured with a weighted UniFrac metric, beta diversity was significantly affected by days post hatch (R-statistic = 0.50, p = 0.001). The caecal community at 0 and 3 d.p.h was different when compared to other time points with more variation between samples taken at the same time. At 7 d.p.h the microbiota was more similar to later time points. Samples from 14 d.p.h continue to cluster separately, however, samples from 21 d.p.h cluster together with 28 and 42 d.p.h samples. There was also a noticeable separation between mucus and luminal samples from 14 d.p.h (Figure 2).

Differential ASVs Between Time Points
Gneiss analysis and taxa plots revealed differential ASV abundance between time points. The feature table was filtered to exclude ASVs with a frequency of less than 102, reducing the number of ASVs included from 847 to 424. For Gneiss analysis, the overall linear regression model fit was R2 = 0.654  Figure 3) these balances revealed waves of colonization in the caecal microbiota. The taxonomic assignments at the level of family of ASVs in these balances are displayed in Table 3 and full taxonomy plots are provided in the Supplementary Material. The log ratio of balance y0 had a high value at 0 d.p.h suggesting that some ASVs in y0 numerator were more abundant at this time ( Figure 4A). Between 3 and 42 d.p.h, the log ratio decreased as there was a shift in relative abundance from    y0 numerator to y0 denominator ASVs. The initial abundance and subsequent decline of y0 numerator ASVs can be observed in the dendrogram heatmap (Figure 3) and the taxa plot ( Figure 5). However, examination of further balances is required to fully distinguish ASVs prevalent at 0 d.p.h from those present at 3, 7, and 14 d.p.h.
Balance y2 is a subdivision of y0 numerator . The log ratio of balance y2 was lower at 0 d.p.h as there was a higher abundance of y2 denominator ASVs ( Figure 4B). There was an increase in log ratio from 7 to 14 d.p.h as the relative abundance of numerator ASVs increased (Figure 3). From 21 d.p.h, the log ratio began to decrease as the balance was shifted back toward 0 by a decrease in relative abundance of y2 numerator ASVs. It can be concluded that ASVs which colonized the caecum between 7 and 14 d.p.h are represented by y2 numerator . Further resolution is provided by balance y6 whose log ratio increased from 0 to 7 d.p.h due to an increase in y6 numerator ASVs at 3 and 7 d.p.h before declining at 14 d.p.h as y6 denominator ASVs began to colonize the caecum (Figure 4C).
Balance y5 is a subdivision of y2 denominator which allows the distinction between ASVs which colonized at 0 and 3 d.p.h. The log ratio of balance y5 increased from 0 to 3 d.p.h (Figure 4D) suggesting that y5 denominator ASVs were more abundant at 0 d.p.h with a higher prevalence of y5 numerator ASVs at 3 d.p.h. This pattern is confirmed by the dendrogram heatmap (Figure 3).
As previously discussed, ASVs associated with later time points are described by y0 denominator , but further balances are required to discern at which time point ASVs became prevalent. Balance y1 is a subdivision of y0 denominator . The log ratio of balance y1 was close to 0 at 0 and 3 d.p.h as both numerator and denominator ASVs were absent. There was a slight decrease at 7 d.p.h as some denominator ASVs colonized the caecum. This decrease continued at 14 and 21 d.p.h before the log ratio began to increase again at 42 d.p.h ( Figure 4E). This suggests that y1 denominator ASVs were associated with colonization between 7 and 21 d.p.h while y1 numerator ASVs represent later colonizers.
Balance y3 is a subdivision of y1 denominator . The log ratio of balance y3 was close to 0 at 0 and 3 d.p.h as discussed for balance y1. There was a decrease in log ratio between 7 and 14 d.p.h followed by an increase from 21 d.p.h ( Figure 4F). This shows that y3 denominator ASVs were associated with colonization between 7 and 14 d.p.h while y3 numerator ASVs represent later colonizers. Balance y7 provides further resolution of y3 denominator ASVs and shows that y7 numerator ASVs were associated with colonization at 7 d.p.h while y7 denominator ASVs began to colonize the caecum from 14 d.p.h ( Figure 4G). Equally, balance y8 provides further resolution of y3 numerator ASVs and shows that y8 numerator ASVs colonized the caecum at 14 d.p.h while y8 denominator colonized from 21 d.p.h ( Figure 4H).
Finally, balance y4 is a subdivision of y1 numerator . The log ratio of balance y4 was close to 0 at 0, 3 and 7 d.p.h. There was a slight decrease in log ratio at 14 d.p.h with a further decrease noticeable at 21 d.p.h (Figure 4I). The dendrogram heatmap shows that this decrease in log ratio can be attributed to an increase in the relative abundance of y4 denominator ASVs. There was an increase in log ratio at 28 and 42 d.p.h attributable to an increase in the relative abundance of y4 numerator ASVs.
In conjunction with the taxa plot (Figure 5), these results demonstrate a pattern of succession within the caecum. ASVs identified as more abundant at 0 d.p.h were assigned to Clostridiaceae 1 (n = 14), Enterobacteriaceae (n = 12) and Enterococcaceae (n = 3). This pattern is visible in the taxa plot where these families are prevalent at 0 d.p.h. One ASV assigned to Lachnospiraceae was differentially abundant at 0 d.p.h. The taxaplot shows a small percentage of Lachnospiraceae and Ruminococcaceae detected in the microbiota at this early time point. Additionally, one ASV assigned to Peptostreptococcaceae was differentially abundant at 0 d.p.h.
Fewer ASVs were differentially abundant at 7 d.p.h but, again, the majority were assigned to Lachnospiraceae (n = 9) and Ruminococcaceae (n = 4). Two ASVs assigned to Bacteroidaceae, a major constituent of the caecal microbiota in this experiment, were first seen colonizing the caecum at 7 d.p.h. One ASV assigned to Atopobiaceae was differentially abundant at 7 d.p.h.

Differences Between the Luminal and Mucosal-Associated Microbiota
Samples from 21, 28, and 42 d.p.h were used to compare lumen and mucus samples from a mature caecal microbiota. Samples from 14 d.p.h were excluded as there were significant differences in alpha diversity between 14 and 42 d.p.h suggesting that the microbiota was not fully developed at 14 d.p.h. This is confirmed by analysis of beta diversity as time post hatch had a small but

Alpha and Beta Diversity
There were no significant differences in alpha diversity between lumen and mucus samples. However, there were significant differences in beta diversity. Mucus and luminal samples formed distinct clusters on a PCoA plot (Figure 2), although there is one exception to this pattern. The lumen sample which clustered with the mucus samples is from Cobb at 42 d.p.h. With reference to the taxa plots (Figure 5), it's clear that the composition of this sample was very similar to mucus taken from Cobb at 42 d.p.h. An ANOSIM test (R-statistic = 0.79, p = 0.001) showed that area sampled had a significant effect on beta diversity.

Differential ASVs Between the Lumen and Mucus Microbiota
Gneiss analysis (Figure 6 and Table 4) and taxa plots (Figure 5) revealed differential ASV abundance between mucus and lumen communities. The feature table was filtered to exclude ASVs with a frequency of less than 81, reducing the number of ASVs included from 748 to 375. The overall linear regression model fit was R2 = 0.56 with covariate "area" accounting for 7.9% of variance. Log ratio balance y2 (β = −11.5, p < 0.001) was a significant predictor for the covariate of "Area." The log ratio of balance y2 was lower in mucus samples compared to lumen samples showing that the relative abundance of y2 denominator ASVs was higher in the mucus than the lumen (Figure 7) with this increased abundance visible in the dendrogram heatmap (Figure 6). ASVs identified as more abundant in mucus were assigned to Ruminococcaceae (n = 25), Lachnospiraceae (n = 20), Christensenellaceae (n = 5) and Bacillaceae (n = 4). Other ASVs were assigned to Eggerthellaceae (n = 2), Enterococcaceae (n = 1), Erysipelotrichaceae (n = 1), Peptococcaceae (n = 1), Coriobacteriaceae (n = 1), Peptostreptococcaceae (n = 1) and Mollicutes RF39 (n = 1).
The dendrogram heatmap does not show a pattern of increased relative abundance of y2 numerator ASVs in lumen samples compared to mucus samples. As a result, it is difficult to discern which y2 numerator ASVs are truly more abundant in the lumen as no other balances are significantly different allowing for greater resolution. However, the taxa plot can aid in this distinction. A paired Student's t-test of relative abundance of Burkholderiaceae (lumen average = 3.97%, mucus average = 2.8%, test statistic = 3.78, p = 0.005) and Bacteroidaceae (lumen average = 45.4%, mucus average = 21.4%, test statistic = 5.97, p < 0.001) show significant differences in relative abundance between lumen and mucus samples. On the other hand, some taxa in y2 numerator show a significantly higher relative abundance in mucus including Bifidobacteriaceae (lumen average = 2.8%, mucus average = 4.8%, test statistic = 4.11, p = 0.003), Atopobiaceae (lumen average = 0.68%, mucus average = 1.53%, test statistic = 4.04, p = 0.003) and Eubacteriaceae (lumen average = 0.02%, mucus average = 0.06%, test statistic = 3.11, p = 0.01). It's possible that since these taxa are represented by one or two ASVs that the limited sample size doesn't provide sufficient statistical power to demonstrate a significant difference in their abundance between mucus and lumen samples using Gneiss analysis. As a result, there is some doubt as to whether these taxa are differentially abundant between the mucus and the lumen.

Differences Between the Lumen and Mucus Microbiota Confirmed by Quantitative PCR
Based on the results from Gneiss analysis and the taxa plot, quantitative PCR assays for Bacteroidaceae, Lachnospiraceae-Ruminococcaceae, Bifidobacteriaceae and Bacillaceae were performed. A combined primer pair for Lachnospiraceae and Ruminococcaceae was used as well as specific primers for Clostridium cluster XIVa&b and Clostridium cluster IV with the former corresponding roughly to Lachnospiraceae and the latter to Ruminococcaceae. Previously published primers were used for Bacteroides and Bifidobacterium detection as well-trialed primer pairs were available in the literature. The primers for the Bacillaceae ASVs were generated as described in the Materials and Methods section. Results showed the same pattern as detected by Gneiss analysis (Figure 8). Bacteroidaceae were more abundant in luminal samples (test statistic = 6.3, p < 0.001) and Bacillaceae (test statistic = −2.4, p = 0.03) were more abundant in the mucus. Bifidobacteriaceae were more significantly abundant in the mucus, although this pattern breaks down at 21 d.p.h (test statistic = −2.5, p = 0.03). There was a significantly higher abundance of Lachnospiraceae-Ruminococcaceae, Clostridium cluster IV and Clostridium cluster XIVa&b in the mucus compared to the lumen (test statistic = −6.7, p < 0.001, test statistic = −9.1, p < 0.001 and test statistic = −5.8, p < 0.001, respectively).
The correlation between abundance determined by qPCR and relative abundance determined by sequencing varied between primers. There was a strong, positive correlation between results from qPCR and relative abundance determined by sequencing for Bacillus (rs (22) = 0.89, p < 0.001), Bacteroidaceae (rs (22) = 0.88, p < 0.001), Clostridium cluster XIVa&b (rs (22) = 0.80, p < 0.001) and Bifidobacteriaceae [rs (22) = 0.72, p < 0.001]. A weaker positive correlation was found between results from qPCR and relative abundance determined by sequencing for Clostridium cluster IV [rs (22) = 0.49, p = 0.02] and Lachnospiraceae-Ruminococcaceae [rs (22) = 0.40, p = 0.05]. FIGURE 7 | Log ratios of balance y2 which was significantly different between lumen and mucus samples. The lower log ratio in mucus samples suggests that the relative abundance of y2 denominator ASVs was higher in mucus samples compared to lumen samples.

Differences Between Breeds
The taxa plots of pooled samples showed some early differences in the caecal microbiota between breeds (Figure 5). At 0 d.p.h, there were large differences in composition between the breeds Bacteroides shows a high affinity for the lumen while Bacillus is found almost exclusively in the mucus. A non-specific Lachnospiraceae-Ruminococcaceae primer pair showed consistently higher abundance in the mucus as do the more specific primers for Clostridium Cluster IV and XIVa&b. Bifidobacterium appears more abundant in the mucus; however this is not consistent across time points.
with Hubbard and Ross mainly colonized by Enterobacteriaceae while Cobb were colonized by Enterococcaceae and Clostridiaceae (Figure 5). By 3 d.p.h, there was increased homogeneity between the breeds although Hubbard had proportionally more Bifidobacteriaceae and less Enterobacteriaceae than the other two breeds. From 7 d.p.h, there were no visible differences between breeds as seen in the taxa plots.

Alpha and Beta Diversity
Using data from 5 samples of caecal mucus from each breed taken at 42 d.p.h, some differences were noted. There was a significant difference in alpha diversity between Cobb and Ross samples (H = 4.84, p = 0.04) and Hubbard and Ross samples (H = 4.84, p = 0.04) but no significant difference between Cobb and Hubbard samples. Breed had a significant effect on beta diversity (R-statistic = 0.25, p = 0.03). Pair-wise tests between breeds revealed that there were no significant differences in beta diversity between Cobb and Hubbard samples and Ross and Hubbard samples. There was a significant difference between Cobb and Ross (ANOSIM: R-statistic = 0.57, p = 0.009).

Differential ASVs Between Breeds
The feature table was filtered to exclude ASVs with a frequency of less than 66, reducing the number of ASVs included from 790 to 398. Pairwise Gneiss analysis between breeds revealed differential ASV abundance between caecal mucus microbiota in Cobb and Ross ( Table 5) and Hubbard and Ross but found no appreciable difference between Hubbard and Cobb.
For the comparison between Cobb and Ross, the covariate Breed accounted for 15.3% of variance. Log ratio balances y2 (β = −13.7, p = 0.01) and y14 (β = 12.2, p = 0.001) were significant predictors for the covariate of breed. The log ratio of balance y2 was, on average, lower in Ross showing a higher relative abundance of denominator ASVs. This difference in abundance is visible on the dendrogram heatmap (Supplementary Figure 3). The log ratio of balance y14 was lower in Cobb. The dendrogram heatmap shows that this is likely due to an increased relative abundance of y14 denominator ASVs in Cobb samples. ASV taxonomy of these balances is displayed in Table 5.
For the comparison between Hubbard and Ross, the covariate breed accounted for 12.4% of variance. Log ratio balance y1 was approaching significance (β = 8.22, p = 0.07) and on review of the dendrogram heatmap (Supplementary Figure 4) was included in the analysis. The log ratio of balance y1 was lower in Hubbard samples. This difference was due to an increased relative abundance of y1 denominator ASVs which can be seen in the dendrogram heatmap. These ASVs were mainly assigned to Ruminococcaceae (n = 21), Lachnospiraceae (n = 7) and Clostridiales vadin BB60 group (n = 6). The remaining ASVs were assigned to Christensenellaceae (n = 2) and one each to Peptococcaceae, Enterococcaceae, Clostridiales Family XIII and Eggerthellaceae.
For the comparison between Hubbard and Cobb, the covariate breed accounted for 9.7% of variance. Log ratio balance y15 (β = 4.13, p = 0.01) was a significant predictor for the covariate of breed. The log ratio of balance y15 was lower in Cobb suggesting that y15 denominator ASVs are more abundant. This is confirmed by the dendrogram heatmap (Supplementary Figure 5). The ASVs composing y15 denominator were assigned to Ruminococcaceae (n = 3) and Coriobacteriaceae (n = 1).

Bacterial Colonization of the Caecum
This study aimed to characterize the succession of bacteria in the caecal microbiota in three broiler breeds between 0 and 42 days post hatch. Although the results of this study must be interpreted in light of small sample sizes which limit statistical power, and the pooling of samples, which obscures the inherent variability of individual microbiota composition, broad patterns are visible within the data and merit further discussion. The initial microbiota observed in Hubbard and Ross breeds was very similar to that described in other studies (8). Large differences in microbiota composition were observed between breeds immediately post hatch. This is most likely due to different bacterial exposure during handling of the chicks from hatch to collection (7). These early differences in microbiota composition did not affect the subsequent development of a mature microbiota. Differences between breeds were still visible at 3 and 7 d.p.h but were no longer visible at 14 d.p.h. This suggests that environmental exposure, diet and management practices are more important than genetics in shaping the caecal microbiota. The rise of Bifidobacteriaceae at 3.d.p.h may represent an important step in the maturation of the caecal microbiota as a stimulus for Bacteroidaceae growth. Recent studies have shown that Bacteroides fragilis and other intestinal bacteria can metabolize exopolysaccharides, a complex carbohydrate produced by some Bifidobacterium strains (36,37). It's possible that an initial rise in Bifidobacterium produces polysaccharides which are then used by other bacteria as a substrate and the basis for the expansion of their populations. Salazar et al. (37) also showed that exopolysaccharides from different strains of Bifidobacterium supported populations of different bacteria such as Faecalibacterium prausnitzii. This suggests that the initial colonizing strains of Bifidobacteriaceae could influence future microbiota development in the caecum.
As well as promoting the growth of other caecal microbes Bifidobacteriaceae play an important role in pathogen exclusion and intestinal barrier function (38,39). Most of the evidence for this comes from mammalian studies, however, a Bifidobacterium probiotic has been shown to improve epithelial integrity in chickens (40). The mechanism may be associated with Bifidobacteriums production of acetate which can protect mice epithelial cells in the face of E. coli O157 infection (41). Bifidobacterium also plays a role in dendritic cell maturation and the balancing of regulatory T cell and T helper 17 cell development (42,43), similar to the role of Candidatus Savagella in the ileum.
One of the caecums main functions is bacterial fermentation of indigestible polysaccharides to produce short chain fatty acids which can be utilized by the host's epithelial cells. The taxa which fulfil this role are certain classes of Firmicutes (such as Lachnospiraceae or Ruminococcaceae) and Bacteroidetes. In this study, they were equally abundant with some Firmicutes localizing to mucus and Bacteroides to the lumen. Previous studies have found that one or the other is more abundant in the caecum (12)(13)(14). This discrepancy between studies may be explained by differences in methodology or sampling. However, it has been reported that even among chickens from the same flock, the Firmicutes:Bacteroidetes (F/B) ratio can vary substantially (44). There is conflicting evidence as to the F/B ratio's impact on metabolism and feed conversion ratio (FCR). Stanley et al. (44) reported that the differences found between chickens in the same flock didn't appear to have an effect on apparent metabolizable energy or FCR. In contrast to this result, a comparison of the faecal microbiota between high and low FCR broiler chickens at 49 d.p.h found significant differences in microbiota composition between the two groups. Low FCR birds had a higher F/B ratio than high FCR birds suggesting that a higher relative abundance of Firmicutes can be associated with increased metabolic efficiency (45). Similar patterns have been noted in the human microbiota with a higher F/B ratio linked to obesity (46). However, Firmicutes is a diverse phylum and closer inspection at a lower taxonomic level reveals a more complex pattern. The two major families of caecal Firmicutes have been identified as Ruminococcaceae and Lachnospiraceae. Singh et al. (45) found that the relative abundance of Lachnospiraceae was nearly two times higher in high FCR birds while the relative abundance of Ruminococcaceae was 15 times higher in low FCR birds.
Lachnospiraceae was the first Clostridia to colonize the caecum in this study. In contrast, Ruminococcaceae was a relatively late colonizer which appeared to replace Lachnospiraceae, especially in the lumen, at later time points. These two families are poorly classified, and material related to their potential role in the microbiota is scarce. Differences in gene abundance between genomes of Lachnospiraceae and Ruminococcaceae have been noted (47). However, the functional significance of these differences in relation to the pattern of succession is not clear. One exception to this rule is Faecalibacterium prausnitzii, a butyrate producing member of Ruminococcaceae (48,49), which has been identified as a potentially beneficial microbe (50). Aside from its prominence as a late stage member of the caecal microbiota little is known about the impact of F. prausnitzii on the chicken. Some inferences can be made from observations reported in mouse and human studies. In these species, F. prausnitzii is thought to have an anti-inflammatory effect since lower counts are observed in various inflammatory diseases (51,52). Additionally, inflammation and intestinal barrier function are improved in a mouse IBD model by the addition of F. prausnitzii (51,53). While the production of butyrate may play a role in F. prausnitziis beneficial effects, the possibility of a more complex mechanism involving the production of other molecules should not be ruled out (54). The results also showed later colonization of the caecum by other members of Firmicutes such as Christensenellaceae, Clostridiales vadin BB60 group, Peptococcaceae and Bacillaceae, as well as Mollicutes RF39. Some, for example Clostridiales vadin BB60 group, are poorly classified and little is known about their metabolism or role in the microbiota. One exception is the family Christensenellaceae which has been linked to lower BMI and overall gut health in humans (55,56).
Another feature of the developing microbiota in this study was the decrease in Enterobacteriaceae over time. This bacterial family is considered of importance in poultry production, not just as pathogens, but also due to their carriage of genetic factors responsible for antimicrobial resistance such as extended-spectrum beta-lactamase genes (57). The relative abundance of Enterobacteriaceae was highest at 0 and 3 d.p.h before decreasing to its lowest relative abundance at 21 d.p.h. A similar pattern of replacement of Enterobacteriaceae by other taxa over the few weeks of life has been observed in previous studies (8,15,58). Reduced abundance and growth of Enterobacteriaceae have been linked to levels of SCFAs like acetate, butyrate and propionate both in vitro and in vivo (59). As such, the reduction in Enterobacteriaceae observed in this study and others could be attributed to caecal colonization by Lachnospiraceae, Ruminococcaceae and Bacteroidaceae. This highlights the importance of these taxa in areas other than metabolism and suggests that interventions which promote early colonization by these taxa should be prioritized.
The order of caecal bacterial succession described in this paper may not be representative of a commercial setting where environmental exposure and other factors will allow for a differing microbiota to develop. However, it does raise questions about the general mechanisms of succession. It's possible that the late arrival of taxa such as Clostridiales vadin BB60 group and Bacillaceae was a question of environmental exposure, that these taxonomic groups were not present in the environment until 14 or 21 d.p.h. However, it also raises the possibility that these taxa have some prerequisite conditions that must be fulfilled, either by the host or earlier bacterial colonizers, before they can establish a population in the caecum. Identifying factors or interventions that accelerate the development of a mature microbiota may provide further benefits to chicken production in terms of increased yields or reduced losses to infectious disease.

The Mucus and Lumen-Associated Microbiota Are Different
Similar differences between mucus and lumen microbiota have been described in other species. Higher levels of Lachnospiraceae and a lower abundance of Bacteroidaceae have been reported in mouse colon mucus when compared to the lumen (32). This may be due to the differing energy sources of these bacteria. As in other animals, bacterial degradation of indigestible polysaccharides into short chain fatty acids (SCFA) plays an important role in host nutrition. In the human gut, Bacteroides have the highest number and diversity of genes associated with polysaccharide metabolism (60). This is also true of chickens with metagenomics studies finding polysaccharide utilization systems associated with caecal Bacteroidetes (1). It would make sense for Bacteroides abundance to be higher in the lumen where their energy source is most abundant. The same reasoning can explain the higher abundance of Lachnospiraceae in mucus. A recent metagenomic analysis of Lachnospiraceae and Ruminococcaceae genomes showed differences in gene abundance related to carbohydrate metabolism. Ruminococcaceae have higher numbers of cellulase and xylanase genes related to fermentation of substrates which are more abundant in the lumen. Lachnospiraceae were found to have a higher number of genes for glycoside hydrolase 13, an enzyme associated with cleavage of α-amylase bonds present in starch and glycogen (47).
Members of the Lachnospiraceae family have also been shown to utilize mucin glycans as a sole carbon source, producing propanol and propionate (61). Due to the preference for chicken Lachnospiraceae to reside in the mucus, it is likely that these strains also possess enzymes which allow them to utilize host mucins as an energy source.
Equally, strains of Coriobacteriaceae and Bifidobacteriaceae have also been found with the ability to degrade mucin in other species (62,63). The ability of Bifidobacteriaceae to adhere to mucus is well documented and also explains higher abundance in the mucus than the lumen (64)(65)(66) The host likely benefits from abundant Firmicutes in mucus. Differences in SCFA production have been observed between Firmicutes and Bacteroides in the chicken caecum with the former producing butyrate while the latter produces mainly propionate (67). Since, epithelial cells utilize butyrate as their principal energy source, a ready source in the mucus is likely to support early epithelial maturation and growth (67).
While these results can be explained by bacterial adaptation to niches based on substrate availability, an alternative explanation may lie in host regulation of the microbiota and maturation of the immune system. In this study, there was a pattern of mucosal colonization by ASVs which appear in the caecum before 14 d.p.h. While the association between appearance in the microbiota and the ability to colonize caecal mucus may be incidental since Lachnospiraceae colonize the caecum before Ruminococcaceae, it may be that early presence in the caecum results in host tolerance of certain ASVs allowing them to occupy the mucus layer.
These results have implications for other studies of the chicken caecal microbiota. For example, studies observing the effects of probiotics containing Bifidobacteriaceae or Bacillus which observe only the luminal microbiota may miss increased abundance of these taxa. Additionally, the Firmicutes:Bacteroides ratio is often used to compare microbiota between groups (46,68,69). One notable study using this parameter found a large range in Firmicutes:Bacteroides ratio among individuals of the same breed housed under the same conditions (44). It is possible that some of this variation could be due to different ratios of lumen content to mucus in samples. This would be particularly noticeable in chickens which had empty caeca at the time of sampling since any content present is likely to have a high proportion of mucus. Many studies of the caecal microbiota do not detail the sampling methodology enough to discern whether samples were taken from the lumen or the mucus. These results further emphasize the importance of a detailed methodology which discerns between sampling from different compartments of the caecum.

Differences Between Breeds
Some differences between breeds were visible in the pooled samples. These differences occurred at 0 and 3 d.p.h with the greatest difference at 0 d.p.h. At subsequent time points there were no significant differences between pooled samples from different breeds. The difference in early microbiota can be explained by the previous finding that chicks from different hatcheries are colonized by different microbes (7). Although the chicks used in this experiment were sourced from a single hatchery, chicks from different breeds are hatched and reared in different buildings resulting in different environmental exposure. As such, the differences between breeds at 0 d.p.h cannot be attributed to genotype since it is likely that chicks are colonized by whichever environmental bacteria are present at the time of hatch and this will vary between hatcheries and within hatcheries.
Analysis of five samples of caecal mucus taken from each breed at 42 d.p.h showed differences in microbiota between breeds. The results of alpha and beta diversity, as well as Gneiss analysis, suggest that Cobb and Hubbard had the most similar microbiota composition with the greatest differences found between Cobb and Ross. Many of these ASVs that were differentially abundant between Cobb and Ross were assigned to Ruminococcaceae and Lachnospiraceae, perhaps as these were among the most prevalent taxa in the caecum. It is of interest that most Lactobacillaceae were more abundant in Ross samples while a high proportion of Christensenellaceae were more abundant in Cobb. Both of these taxa have been linked to metabolic effects on the host. Dietary supplementation with Lactobacillus strains has been previously demonstrated to improve body weight gain and FCR in chickens (70) while Christensenellaceae is more abundant in the microbiota of lean compared to obese humans (55,56). However, in the absence of detailed performance data such as FCR or body composition analysis, it is not possible to ascribe a functional or practical significance to these results. Greater taxonomic resolution or metagenomic analysis would also be required as the metabolic effects of bacteria is likely to be strain dependent.
Another factor which may have affected these results is the relative composition of lumen contents to mucus of the samples. It was noted earlier that the pooled lumen sample from Cobb at 42 d.p.h was the only lumen sample with a similar beta diversity to mucus samples. Since it has already been demonstrated that there are significant differences between lumen and mucus microbiota, the observed difference could be attributed to a random sampling error amplified by the small sample sizes used in this study. It is possible that samples taken from Ross had more content or less mucus compared to those from Cobb, accounting for the differences between them. Caecal emptying occurs several times a day in chickens. A recently emptied caecum would have a higher proportion of mucus to luminal contents when expressed during sampling which may yield different results had the caecum been full. With this in mind, while this experiment showed some differences in microbiota composition between genotypes, clear results are hampered by small sample size. However, the influence of genotype on caecal microbiota cannot be discounted.

CONCLUSION
This study observed the development of the caecal microbiota between hatch and 42 d.p.h, covering the lifespan of a modern broiler chicken. After hatch, the microbiota had poor diversity and was mainly composed of environmental bacteria. Between hatch and 21 d.p.h, the microbial community became more complex and matured to a stable, diverse microbiota. It is worth noting that for more than half of the production period the microbiota was developing and would be more susceptible to changes brought about by external factors. Earlier interventions in the microbiota are likely to be more successful, both in terms of altering the mature microbiota and optimizing the microbiota's impact on immune system maturation and metabolic function. Interventions should focus on promoting early maturation particularly with respect to Bacteroidaceae, Lachnospiraceae and Ruminococcaceae. Although some significant differences were found in the microbiota between breeds, it is possible that these differences were caused by other factors such as a higher relative composition of luminal content to caecal mucus.

AUTHOR CONTRIBUTIONS
PR conducted the sampling, sample processing, 16S rRNA gene analysis, interpretation of the results, and wrote the manuscript. MB advised on the design of the experiment and edited the manuscript. JF advised on experimental design and data analysis relating to 16S rRNA gene sequencing. PW advised on experimental design and helped conduct the experiment. All