Extensive microbial diversity within the chicken gut microbiome revealed by metagenomics and culture

1The domestic chicken is the most abundant bird and most abundant food animal on Earth, accounting for a larger fraction of the planet’s biomass than all species of wild birds combined. Consumption of chicken meat is growing faster than any other type of meat and is seen as a cheaper, healthier, low-carbon alternative to meat from mammalian livestock. Chicken eggs remain a nutritious, affordable food across the globe.

2The chicken gastrointestinal tract is home to a complex community of microbes and their genes—the chicken gut microbiome—that underpins links between diet, health and productivity in poultry, as evidenced by the ability of antibiotics to promote growth in chicks. This microbial community also acts as a source of pathogens associated with disease in birds or in humans—including Campylobacter, Salmonella, and Escherichia coli—as well as providing a reservoir of antimicrobial resistance genes.

3Previous studies of this community have documented a rich variety of microorganisms (dominated by bacteria, but including viruses, archaea and microbial eukaryotes) and have shown that the taxonomic composition of this community varies with age, breed and disease status. However, these earlier efforts have largely relied on analyses of molecular barcodes (in particular short 16S rRNA gene sequences), which fail to provide species-level resolution, are unable to detect viruses and reveal nothing about the genome sequences, population structures or functional repertoires of microbial species.

4Two strategies have proven productive for exploring taxonomic and functional diversity in complex microbial communities. Culture-independent approaches rely on shotgun metagenomic sequencing of DNA extracted from relevant samples, followed by bioinformatics-based community profiling and analysis. Culture-dependent approaches combine large-scale isolation of microorganisms in pure culture with whole-genome sequencing and phylogenomic analysis. To explore taxonomic novelty in the chicken gut microbiome, we generated phylogenetic profiles to document known and unknown diversity and then exploited culture-dependent and culture-independent approaches to create an unprecedented high-quality reference collection of microbial genes and genomes from the chicken gut, revealing and naming hundreds of new candidate species from this commonplace but important ecological setting.

5Faecal samples were collected in South-East England from adult Lohmann Brown laying hens and adult Silkie hens in 2018. Birds were housed in a large outdoor run with a substrate of stone chippings and small turf enrichment beds during the day and kept in a coop overnight. They were fed a commercial layer feed, Farmgate Layer pellets and mash (ForFarmers UK Limited, Rougham, Bury St Edmunds), according to the manufacturer’s instructions and no antibiotics were used. Faecal sampling was approved by the University of Surrey’s NASPA ethics committee.

6Sixty faecal samples were collected from the Lohmann Brown laying hens and thirty samples from the Silkie hens (six and three samples per day, respectively, for 10 days). Freshly evacuated faeces from individual birds were collected in sterile containers and immediately stored at −20 °C. Samples were then transferred to the laboratory for culture and/or DNA extraction. DNA was extracted using DNeasy PowerSoil kit (Qiagen, Hilden, Germany), following manufacturer’s instructions and then stored at −20 °C.

7Workflow from this point forward is summarised in Fig. 1.The fifty samples yielding >20 ng DNA were processed according to the Low Input, Transpose Enabled (LITE) library construction pipeline before being subjected to paired-end (2 × 150 bp) metagenomic sequencing on the Illumina Novaseq 6000 platform. Bioinformatics analyses were performed on the Earlham Institute’s High Performance Computing cluster and on the Cloud Infrastructure for Microbial Bioinformatics. Sequences were assessed for quality using FastQC Version 0.11.8 and trimmed using Trimmomatic Version 0.36, configured to a minimum read length of 40, ‘leading’ and ‘trailing’ settings of 3 (SLIDINGWINDOW:4:20) (http://www.bioinformatics.babraham.ac.uk/projects/fastqc; ). Metagenomic sequences for all samples have been uploaded to the Sequence Read Archive under Bioproject ID PRJNA543206.

8An initial analysis of our chicken faecal sequences using the Kraken 2 taxonomic classifier was performed on custom databases representing the domestic chicken genome (GenBank assembly accession GCF_000002315.6) and the food plants Triticum aestivum (wheat), Aegilops tauschii (diploid progenitor of the D genome of hexaploid wheat) and Glycine max (soy bean): GenBank assembly accessions GCF_001957025.1, GCA_900519105.1, GCF_000004515.5. Kraken 2 revealed that 8% of reads originated from the chicken and at least 16% originated from the diet. These sequences were filtered from our dataset and excluded from subsequent analyses by keeping only reads ‘Unclassified’ by Kraken 2 after comparison with each database in turn.

9The remaining dataset underwent taxonomic profiling using Kraken 2 against a microbial database built from all complete/representative archaeal, bacterial, fungal, protozoan, viral and UniVec_Core sequences in RefSeq in January 2020. Bracken was used to estimate taxon abundance from the Kraken 2 profiles, accepting only those taxa with ≥1,000 assigned reads. Bracken-database files were generated using ‘bracken-build’ on our microbial database and visualised using Pavian.

10We searched the NCBI BioProjects database (https://www.ncbi.nlm.nih.gov/bioproject/) in November 2019 with the term ‘chicken gut microbiome’ and then selected nine publicly available projects that contained at least one metagenomic sequence dataset >1 GByte in size: PRJEB33338, PRJNA193217, PRJNA291299, PRJNA375762, PRJNA415593, PRJNA417359, PRJEB22062, PRJNA543206, PRJNA417359, PRJNA385038, PRJNA616250. Only four of these studies were linked to research publications at the time of publication

11All shotgun metagenomic reads were quality-filtered by removing reads shorter than 70% of the maximum expected read length (100 bp, 250 bp for MiSeq data), an estimated accumulated error >2.5 with a probability of ≥0.01 or with an observed accumulated error >2, or >1 ambiguous position to assist assembly. If base quality dropped below 20 in a window of 15 bases at the 3′ end, or if the accumulated error exceeded 2, reads were trimmed. All these filter steps are integrated in sdm. Reads mapping to the chicken genome and diet were removed from the metagenomic data as described previously, classifying reads with Kraken 2 against custom databases built on the aforementioned genomes.

12Sequence datasets from our fifty samples—together with 582 samples from the selected BioProjects—were assembled using MegaHIT under the option ‘--k-list 25,43,67,87,101,127’. To avoid artefacts that sometimes result from co-assembly of sequences from different samples and different sources, we performed individual assemblies on each sample, with the exception of BioProject PRJNA417359. For that BioProject, as multiple metagenomic samples had been sourced from different tissues of the same individual bird, we co-assembled reads from the 120 BioSamples from that project.

13Contig sequences from the MegaHIT assemblies of our 50 samples that were ≥10 kb were analysed with VirSorter v1.0.5 with the ‘-db 2’ option to identify viral genomes. VirSorter Category 1 and 2 contig sequences were collapsed at 95% nucleotide identity over 70% of the sequence length using CD-Hit Est v4.6.1. Classification of bacteriophage sequences relied on nucleotide searches using BLASTN against the NCBI NT database (Completed April 2020) and protein searches using Kaiju Version 1.7.3 against the RefSeq database (Completed April 2020). Only bacteriophage genomes with BLASTN hit E-Value < 0.05, percentage identity >70% and query covering >50% were selected as reliable hits.

14A taxonomic assignment was drawn from the highest scoring BLASTN (or in rare cases BLASTP) hit ranked by query cover and percentage ID. Synteny between predicted coliphages and their respective reference genomes were visualised using EasyFig. Escherichia bacteriophage coverage per sample was determined using Anvi’o v6.1 using default parameters and visualised in R using the Pheatmap package (https://www.rdocumentation.org/packages/pheatmap). Remaining viral genomes were filtered for completeness, retaining those that were circular and encoded a complete terminase gene (as predicted by VirSorter). Taxonomic assignments to the level of family were performed on viral genomes using Demovir (https://github.com/feargalr/Demovir).

15Complete genes identified by Prodigal v2.6.1 were clustered at 95% nucleotide identity using CD-HIT-Est v4.6.1. Incomplete genes were then mapped to this complete gene list using Bowtie2 v 2.3.4.1 and any mapping at 95% nucleotide identity were incorporated into the relevant gene clusters. Finally, genes representing the 40 conserved marker genes defined by were clustered separately and then merged with the existing set of gene clusters. We thus obtained a gene catalogue of >20 million genes, defined as non-redundant at 95% average nucleotide identity (ANI). The final gene catalogue was uploaded to FigShare (https://doi.org/10.6084/m9.figshare.13116809.v4)

16Prodigal was applied in metagenome-mode to all contigs from the MegaHIT assemblies. Unfiltered reads from each sample were mapped against their respective assembly to provide an estimate of contig and gene abundance using Bowtie2 with the options ‘--no-unal--end-to-end -score-min L, −0.6,−0.6’. Samtools 1.3.1 was used to sort and index all resulting Bam files. Only reads with mapping quality >20, >95% nucleotide identity and >75% overall alignment length were retained. BEDTools v2.21.0 was used to create depth profiles from the Bam files. These depth profiles were then translated with rdCover (https://github.com/hildebra/rdCover) into average coverage (in a 50 bp window) per contig or per gene predicted from each contig. Bam files were translated to abundances using the ‘jgi_summarize_bam_contig_depths’ script from the MetaBAT 2 package.

17Gene abundances were linked to their respective gene clusters and originating samples. Redundant genes representing the same orthologue were removed.

18We identified metagenomic species (MGSs) using the combinatorial approach described by, incorporating single-assembly binning in the creation of metagenome-assembled genomes (MAGs), gene catalogue binning in the creation of canopy clusters and hierarchical clustering of candidate genes using the R function hclust, method = complete. To start with, we used MetaBAT 2 v2.15 to bin contigs ≥400 bp. These were quality filtered using CheckM v1.0.11 to obtain 5,595 bins at >80% completeness and <5% contamination.

19Species-level clusters were formed using a combination of two distinct approaches. One approach removed redundancy between samples by pre-clustering bins if ≥30% of their genes overlapped with a higher-quality bin to create a set of pre-MGS bins. Lower-quality bins (>60% completeness and <10% contamination) were also included in the analysis but were not used to form new species clusters. To recover prokaryotic species usually obscured using single-sample assemblies and conventional binning techniques, we refined all species bins into ‘hcl-clusters’ using gene correlations and hierarchical clustering, as described by. We chose genes occurring in ≥10% of all associated MAGs as representatives for each pre-MGS bin and used these to fish for additional co-occurring genes from the gene catalogue, using a threshold of >0.75 Pearson correlation and >0.85 spearman rho to identify gene co-occurrences within this core gene set. We then merged MetaBAT 2 bins, canopy bins and co-occurring genes into our species bins. We used the presence of 40 known single-copy marker genes, without duplicates, as a quality criterion in selection of sub-clusters, before extracting the final set of MGS gene representatives using MATAFILER (https://github.com/hildebra/MATAFILER). The final collection of MGS bins (canopy clusters + hcl-clusters) was re-assessed for contamination and completeness using CheckM, so that we could be confident that each bin represents a single species.

20A second approach de-replicated all MAGs at 95% ANI (species-level) and 99% ANI (strain-level) using dRep Version 2.0 and only species not identified in approach one were added to the resulting non-redundant species catalogue. The minimum aligned fraction used during ANI genome alignment was 60%. A single representative MAG for each novel species cluster was uploaded to NCBI SRA under BioProject PRJNA543206 and all MAGs generated were uploaded to FigShare (https://doi.org/10.6084/m9.figshare.13116809.v4.). CompareM Version0.1.1 (https://github.com/dparks1134/CompareM) was used to calculate average amino acid identity (AAI) when identifying novel genera, using a cut-off of 60% for the percentage identity and 70% for the minimum alignment length used to delineate genus boundaries.

21We used the Genome Taxonomy Database Toolkit (GTDB-Tk Release 95) to perform taxonomic assignments on strain-level dereplicated MAGs. In addition, genes from each MGS were analysed through GTDB-Tk (Release 95), proGenomes resource and underwent k-mer-based taxonomic profiling using Kraken 2. In assigning taxonomy, we allowed GTDB assignments to take precedence—only when no GTDB taxonomy was available would we adopt taxonomies assigned by ProGenomes and Kraken 2 and, then, only where genus and family assignments from these sources matched. When exploiting the taxonomy assigned according genes from metagenomic species, we applied a least-common-ancestor approach to unplaced taxa at higher taxonomic levels. Species distribution analyses were conducted using the Vegan package in R, before visualisation using ggplot2 and Pheatmap R packages (https://www.rdocumentation.org/packages/pheatmap). Pan-genome analysis was conducted using Roary v3.11.2 and visualised using the roary2svg.pl script. Comparison of our derived metagenomes with those of was performed at 95% ANI using dRep and visualised using web-tool BioVenn.

22To estimate species richness and diversity, the Phyloseq package of R was applied to the output from Bracken on all of our chicken faecal metagenomic datasets. The six faecal samples that showed highest species richness and taxonomic diversity were selected for culture-based studies. Frozen faecal samples were thawed, vortexed and two 0.5 g aliquots (once processed aerobically, the other anaerobically) from each sample were suspended in 5 ml PBS. Each aliquot was vortexed until homogenised, before performing serial dilutions in duplicate down to 1 × 10−5. Processing of samples for aerobic and anaerobic culture was identical, except that, for anaerobic culture, all culture media, diluent and consumables were pre-reduced to anaerobic conditions for at least 24 h before faecal samples were processed in a Whitley A95TG workstation.

23For dilutions 10−3–105, 200 µl was plated directly on to a set of three agar plates for each culture medium (Brain Heart Infusion, Colombia Blood Agar, Yeast extract, casitone and fatty acid) with or without vancomycin supplementation at a concentration of 6 µg/ml (Table S1). Cultures were incubated at 37 °C for 72 h in their respective conditions before assessment of colony growth. Well-isolated colonies were picked according to colonial morphotype distinctive in colour, shape and size, before being re-streaked on to the growth medium from which they were sourced to confirm purity. Individual colonies were subsequently used to inoculate 2 ml of broth based on the source culture medium, incubated at 37 °C for a further 24 h before bacterial DNA extraction. All isolates were archived at −80 °C in glycerol at 20% concentration.

24Genomic DNA was extracted using a DNeasy UltraClean DNA isolation kit according to the manufacturer’s instructions (Qiagen, Hilden, Germany). DNA was quantified using a Qubit® fluorometer (Invitrogen, Carlsbad, CA, USA) high-sensitivity assay, before dilution to the required concentration in RNase-free water and purification on AMPure XP beads (Beckman Coulter, Brea, CA, USA). Sequencing libraries were prepared from 0.5 ng/µl of RNA free genomic DNA. A total of 282 isolates were included for genomic sequencing using the Nextera-XT DNA sample preparation kit (Illumina, San Diego, CA, USA) and whole-genome sequencing performed using the Illumina NextSeq sequencing platform, generating paired-end reads (2 × 150 bp).

25Paired-end reads were quality-assessed and trimmed using FastQC and Trimmomatic as described above. Trimmed reads were assembled into contigs using SPAdes version 3.13.1. Contigs shorter than 500 bp were discarded from analysis. Genome contamination and completeness was assessed using CheckM version 1.0.13. To confirm assembly quality, only genomes conforming to all the following criteria were included in further analysis: (i) contig N50 of >20 kbp (ii) 90% of assembled bases at >5× read coverage (iii) completeness of >95% (iv) contamination of <5% (v) complete 16S rRNA gene sequence.

26Barrnap Version 0.9 (https://github.com/tseemann/barrnap) was applied to all genomes that passed the quality filters to extract full-length 16S rRNA gene sequences. These were then compared to NCBI 16S rRNA gene sequences from RefSeq genomes using the NCBI’s web-based BLASTN facility. 16S rRNA gene sequences that showed an identity of <98.7% to known sequences were assigned to novel species, using the conservative approach in proposed minimal standards. We used ReferenceSeeker Version 1.6.2 to determine ANI and conserved DNA values compared to RefSeq bacterial genomes (Completed March 2020). Genomes that showed ANI ≤95% and conserved DNA ≤69% to the closest relative were designated novel species. The Genome Taxonomy Database Toolkit (GTDB-Tk Release 95) was used to perform taxonomic assignments on isolate genomes. Genomes were clustered at 95% and 99% ANI before selection of a single representative isolate per species using dRep. Where a genome previously designated as novel clustered with a genome of assigned taxonomy, this taxonomy was then applied to the previously designated ‘novel’ genome. Final taxonomic assignments were based on genome-based ANI values derived from RefSeq and GTDB—with GTDB assignments taking precedence. A single representative genome for each novel or renamed species cluster was uploaded to NCBI SRA under BioProject PRJNA543206 and all genomes alongside respective 16S rRNA gene sequences generated were uploaded to FigShare (https://doi.org/10.6084/m9.figshare.13234556).

27For phylogenetic analysis of all MGS and genome sequenced isolates we used PhyloPhlAn v3.0.58 with the ‘diversity high’ and a proteome input predicted from all genome sequences using Prodigal v2.6.1. Diamond v0.9.34 was used to perform a search against 400 universal PhyloPhlAn markers. MAFFT v.7.271 was used to perform multiple sequence alignment before refinement with trimAl v.1.4 and reconstruction into trees using FastTree v2.1 and RAxML v. 8.2.12. All trees were visualised and annotated manually using the online iTOLv5.7 platform. Trees were scrutinised to confirm that species and genera were monophyletic. Phylogeny for all cultured genomes unassigned at species level was confirmed as previously described against all available reference proteomes of that respective genus downloaded from NCBI.

28To investigate the phylogenetic placement of cultured isolates designated as Escherichia marmotae and Escherichia sp001660175 by GTDB, we constructed a core genome phylogenetic tree. The genomes from cultured isolates were compared to genomes representing the full diversity of the genus Escherichia. Three Salmonella genomes were included as an outgroup. The genome sequences were aligned using Mugsy, and alignment blocks conserved across all genomes were concatenated to produce a core genome alignment. A phylogenetic tree was constructed by maximum likelihood with 100 rapid bootstrap replicates, using the general time reversible model of nucleotide substitution with gamma correction for rate heterogeneity, as implemented in RAxML version 8.2.12.

29We collected faecal samples from 90 chickens reared in the UK belonging to two breeds: Lohman Browns (n = 60) and Silkies (n = 30). Short-read sequencing of 50 of these faecal samples generated a metagenomic dataset in excess of a billion paired-end reads or three hundred billion base pairs (Table S2).

30We initially analysed the faecal samples using the k-mer-based programme Kraken 2, followed by refined phylogenetic analysis using the allied programme Bracken (Table S3). Unsurprisingly, these programmes assigned sequence reads from the faecal samples to all three domains of life, as well as to viruses (Table S4), although relative abundance assignments show that bacteria predominate in this environment. Sequences were assigned to a wide range of bacterial phyla, including the three expected as predominant in the vertebrate gut (Bacteroidetes, Firmicutes, Proteobacteria), but also including over twenty additional phyla. Searches of the PubMed database with each phylum name and the term ‘chicken’ reveal that round half of these have been previously documented in the chicken gut. However, at least a dozen appear to be novel in this setting, including the Aquificae, Balneolaeota, Calditrichaeota, Chlorobi, Dictyoglomi, Fibrobacteres, Gemmatimonadetes, Ignavibacteriae, Kiritimatiellaeota, Lentisphaerae, Nitrospirae, and the Thermodesulfobacteria.

31When we rank-ordered the species identified by Bracken according to maximum abundance in any one sample, we found, as expected, that species from the family Lactobacillaceae dominated among the top 20 most abundant organisms. However, we found that two species of Escherichia— Escherichia coli and Escherichia marmotae—accounted for ≥5% of reads in nearly half of the samples (22/50) and in two samples, accounted for more than 50%. Such monodominance of the gut microbiome by bacterial species has been described in diseased humans, but is surprising in the context of poultry reported as apparently healthy by their handlers. We also noted a high relative abundance of the recently described chicken pathogen Gallibacterium anatis in most birds (with four birds showing >5% reads assigned to this organism), despite their healthy status. Similarly, Fusobacterium mortiferum—an opportunistic pathogen of humans—accounted for >10% of sequences in 11 birds, corroborating a recent report of high abundance of 16S rRNA gene sequences from this organism obtained from the chicken caecum.

32Bracken assigned sequences to over a hundred bacteriophage genomes, predominately phages infecting members of the Enterobacteriaceae assigned to the families Myoviridae and Podoviridae. Particularly noteworthy was the high abundance of reads in some samples from two distinct bacteriophages that prey on E. coli: phiEcoM-GJ1—a lytic bacteriophage isolated in Canada from pig sewage—which accounted for 6.5% reads in a single sample and phAPEC8—a lytic bacteriophage with a large 147 kb genome, isolated from a Belgian poultry farm—which accounted for 10% of reads in a single sample and for >1% of reads in three others.

33Although these k-mer-based analyses can provide interesting insights into taxonomic diversity within the chicken gut, we quickly realised that they provide an incomplete and misleading picture of this important microbiome for several reasons: (1) they often report the presence of highly implausible organisms—for example, Kraken 2 reported the presence of human pathogens such as Shigella flexneri and Plasmodium falciparum that are simply not credible in this context on clinical grounds; (2) as with studies on 16S rRNA gene sequences, they fail to provide genomic data or insights into the functional diversity or population structure of the microbial species that they identify and; (3) they rely on a reference database and so can only report previously known organisms and can never uncover ‘unknown unknowns’.

34The scale of the problem of unknown diversity is clear from the observation that nearly three quarters (73%) of sequence reads from our chicken samples cannot be confidently classified by Kraken 2 to species level and more than half of the reads (52%) cannot be classified at all and are simply designated as ‘Unassigned’. We therefore sought to extend our understanding of this community through two powerful reference-free approaches: assembly-based metagenome analyses and high-throughput culture.

35Assembly of metagenomic sequences is a reference-free approach that involves aligning and merging short sequence reads into long contiguous sequences (contigs) contigs.

36To confirm the presence of bacteriophages inferred through the reference-based analysis and to identify novel viral genomes, we assembled sequence reads from our fifty chicken faecal samples into contigs. Contigs ≥10 kb were analysed with VirSorter—a programme designed to detect viral signals in microbial sequence data to find novel viruses.

37VirSorter identified 184 of our chicken faecal contigs as Category 1 (‘most confident’) bacteriophage sequences and identified an additional 1,840 contigs as Category 2 (‘likely’) bacteriophage sequences. This was de-replicated to 1,455 genomes using similarity thresholds of 95% ANI over 70% of the genome (Table S5). BLASTN analysis revealed only 10 of these bacteriophage genomes showed high similarity (percentage identity > 70%; query covering > 50%) to known phages at the nucleotide level (Table S6). These included close relatives of the two phages (phiEcoM-GJ1 and phAPEC8) found highly abundant in the Bracken analyses (Fig. 2). Interestingly, more than one genus of coliphage (e.g. Jilinvirus, Phapecoctavirus, or Gamaleyavirus) was often detected in the same sample, along with an abundance of reads from their predicted prey (Escherichia) suggesting interesting dynamics in phage-host and phage-phage interactions (Fig. 3; Table S7).

38Of the remaining 1,445 unclassified bacteriophage genomes, nearly 600 encoded either an obvious terminase region or were circular and as such were suggested as being near-complete. Classification of these genomes revealed all genomes were predicted to belong to the order Caudovirales of tailed phages, with the majority belonging to the family Siphoviridae (n = 429), but we also found representatives from the Myoviridae (n = 87) and Podoviridae (n = 27), plus some bacteriophages unclassified at family level (n = 28) (Table S8).

39Next, we subjected our samples to computational binning—a process of grouping contigs on the basis of sequence composition and depth of coverage into discrete population bins representing metagenome-assembled genomes (MAGs). However, to carry out a definitive survey of bacterial and archaeal diversity in the chicken gut microbiome—in addition to analysing the fifty faecal samples mentioned and before we started the binning—we retrieved all publicly available chicken gut metagenomic datasets, to create an expansive dataset representing >630 samples, drawn from ten studies and twelve countries (Belgium, China, France, Germany, Italy, Malaysia, Netherlands, Poland, Spain, The Gambia, UK, USA) (Figs. S1A and S1B; Table S9).

40Sequence assembly and binning on all these samples generated 5,595 MAGs that passed our quality threshold of ≥80% completion and ≤5% contamination (Fig. S1C). Of these 3,131 could be considered high-quality draft genomes, with >90% completion and <5% contamination, as judged by recently published criteria (Table S10). Genome sizes of the MAGs ranged from ~0.4 to 6.8 Mbp, while GC content ranged from 24% to 73%.

41Then, we grouped the MAGs into metagenomic species (MGSs). Initially, this involved de-replicating MAGs at the widely accepted 95% average nucleotide identity (ANI) for defining bacterial and archaeal species and 99% ANI for defining bacterial and archaeal strains. De-replication of MAGs at 95% ANI resulted in 846 clusters representing bacterial and archaeal species, while de-replication at 99% ANI resulted in 2182 clusters, representing strains. However, to improve recovery of MAGs, MGSs and associated gene sets, we used gene correlations to identify species-representative genes and then applied hierarchical clustering to co-occurring genes across the samples. This allowed us to identify additional genes from the core genome of a species, even when they show divergent nucleotide compositions (such as genes from genomic islands and plasmids). Similarly, using canopy clustering, we could identify commonly occurring species of low abundance. Using these approaches, we were able to identify an additional seven MGSs (Table S11). These MGS were prevalent at >1× coverage in 53% of all analysed samples spanning at least 4 different BioProjects.

42Analysis of bacterial metagenomic species, primarily using the Genome Taxonomy Database (GTDB) taxonomy, confirmed and extended the taxonomic novelty uncovered by reference-based community profiling (Fig. 4), recovering species spanning nineteen of the bacterial phyla defined by GTDB (Table S12). These include Cyanobacteria (12 species; 32 strains); Deferribacterota (1 species; 1 strain) Synergistota (2 species; 5 strains) and the Verrucomicrobiota (7 species; 8 strains).

43Of the 853 de-replicated bacterial metagenomic species, 321 represented previously delineated species catalogued in publicly available databases (Table S13). Following direct comparison, a further 165 metagenomic species had been previously identified by, with these sequences not currently available in public archives. However, only 158 of our metagenomic species possess validly published names based on Latin binomials.

44We performed a search of PubMed with the species name and ‘chicken’, leaving aside the 33 species named by. This suggested that our study provides the first-evidence-in-chickens for the majority (81/125) of these species (Table S14). Examples include: Jeotgalicoccus halophilus, first isolated from the traditional fermented seafood, Jeotgal and present in 197 chicken samples; Aliicoccus persicus, first isolated from a hypersaline lake and present in 241 chicken samples; and Bacteroides reticulotermitis, first isolated from the gut of a termite and present in 39 chicken samples.

45We found that 309 of our metagenomic species could be assigned a taxonomy only at the level of genus and so represent novel candidate species. A further 56 species could be assigned a taxonomy only at the level of family and, after AAI clustering at 60%, were assigned to 36 novel candidate genera. One candidate bacterial species could be assigned a taxonomy only at the level of order (Oscillospirales) and so represent a new family.

46Three MAGs were assigned to the domain Archaea. One represents the species Methanobrevibacter woesei—which is already known to inhabit the chicken gut—while the other two represent novel species within the genera Methanocorpusculum and UBA71, which we have renamed Candidatus Methanospyradousia.

47Linnaeus first proposed the assignment of Latin binomials to provide a universal nomenclature for biological species. The International Code of Nomenclature of Prokaryotes (ICNP) sets the rules for naming prokaryotic species, but currently precludes the valid publication of names of uncultivated organisms, represented by MAGs or other sequences. Furthermore, high-throughput generation of MAGs and of sequence-based taxonomies for bacteria, such as the GTDB is often assumed to preclude the detailed attention usually given to one-by-one construction of Linnaean binomials. As a result, most uncultured taxa, as well as many taxa defined on sequence-based criteria, have been assigned unstable, confusing and hard to-remember alphanumerical identifiers.

48To provide a stable, clear and memorable nomenclature for novel and/or previously unnamed bacterial and archaeal species from the chicken gut, we exploited the provision within the ICNP for naming uncultivated taxa via Candidatus assignments, which, although provisional, provide the scientific community with well-formed Latin binomials. However, this prompted us into an unprecedented effort to create hundreds of new names for the purpose of this single research study—an effort that required us to devise a scalable combinatorial system for the creation of binomials. Here, we made extensive combinatorial use of several dozen Latin and Greek roots pertaining to poultry (avi-, galli-, pulli-, alektryo, ptero, kotto-, ornitho-), intestines (intestini- entero-), faeces (faec-, kakke, merd-, kopro-, excrement-) or microbial life (-monas, -bacterium, -microbium, -coccus, -bacillus, -bium, -cola)—twinned with addition of these roots (singly or in tandem) and/or prefixes (allo, hetero, meta-, para-, crypto-) to existing genus names—to create over 150 Candidatus genus names. For genera with alphanumeric designations in GTDB Release 05-RS95 known to occur also in gut microbiomes of other animals, we adopted a similar combinatorial approach, but avoided roots pertaining to poultry and stuck instead with combinations that simply meant ‘gut or faecal microbe’, for example Fimicola, Caccocola. An additional source of diversity stemmed from repetitive use of around forty Candidatus species epithets built from similar roots, which when combined with genus names gave us a total of over 650 distinctive binomials for new Candidatus species (Table 1; Table S15).

49To extend our metagenomics analyses, we applied culture-based methods to six faecal samples that appeared species-rich in Kraken 2 analyses and in so doing obtained 282 isolates from aerobic culture (~80% of isolates) and anaerobic culture (~20% of isolates) (Table S16). All isolates underwent genome sequencing on the Illumina platform and phylogenetic analysis to enable taxonomic assignment. The resulting chicken gut culture collection was found to contain 56 genera, 93 species and 162 strains drawn from five phyla. These included thirty novel species, with all novel species confirmed to originate from a monophyletic group through phylogenetic analysis against all available reference genomes of their respective genus (Fig. S2). Curiously, there was no overlap between the species that we obtained and those reported by, suggesting that we are far from exhausting the set of species that can be cultured from this habitat. As with the metagenomic species, all novel or previously unnamed genera and species from cultured isolates were assigned Linnaean binomials (Table 2; Table S17). Species-level ANI clustering of all MAGs and all cultured isolates according to phylum is provided in Fig. S3.

50Interestingly, alongside ten cultured isolates of the well-characterised species Escherichia coli, we recovered three isolates from Escherichia marmotae (a species recently described in Himalayan marmots). As previously reported, the E. marmotae strains cluster closely with the Escherichia Clade V, so all members of this clade should be considered members of this species (Fig. 5; Table S18). Further analysis of the GTDB species designated Escherichia sp001660175 (https://gtdb.ecogenomic.org/searches?s=al&q=sp001660175) confirmed that this species forms a monophyletic lineage that corresponds to Clade II, among the cryptic environmental clades described by, which has subsequently been documented in birds. As Clade II is comparable in divergence to the other Escherichia spp. and cryptic clades, we have therefore assigned the Linnaean binomial Escherichia whittamii to designate a new species (Table 2), honouring the outstanding contribution of Thomas S. Whittam to the study of Escherichia spp..

51We found that only 16 species were common to our cultured isolates and our MGS. Subsequent sequence mapping allowed us to detect a further two cultured species at ≥1× coverage in at least one metagenomic sample (Fig. 6A; Table S19). The genomes from cultured isolates were on average 20% larger than the corresponding MAG sequences retrieved from the same source sample (Table S20), which is in line with the completeness threshold of 80% we adopted in quality assurance of the MAGs. However, when we performed detailed gene content analyses on three abundant species in both cultured and metagenomic datasets—Lactobacillus reuteri (with the synonym Limosilactobacillus reuteri), Escherichia coli (including the synonym Escherichia flexneri) and Enterococcus faecium—we found that >99% of the genes from the core genomes and nearly half of the genes in the accessory genomes of cultured species were represented in at least one MAG. These observations suggest that our high-quality MAGs are sufficiently complete to warrant Candidatus names.

52We analysed our chicken faecal metagenomes with a Kraken 2 database derived from genomes representing our candidate metagenomic and cultured species; this yielded a considerable improvement in the number of reads that can be classified through rapid phylogenetic profiling (Fig. 6B).

53An analysis of the distribution of 820 MGSs across the entire metagenomic dataset revealed marked variation between samples, with not a single species present at ≥1× coverage in all samples and only 39 species present in >90% of samples—although 441 species were present in >50% of samples at ≥1× coverage (Table S21). At ≥1× coverage, co-occurrence of nearly 300 species (n = 295) was identified across all 10 BioProjects (Fig. 7A), with no species identified in all BioProjects at ≥10× coverage (Fig. 7B). Focusing on samples from distinct anatomical or physiological sites (faeces, caecum etc.), we found no species present in all faecal samples at >1× coverage and only two species were found in all caecal samples at >1× coverage: both of them newly named in this study: Candidatus Paraprevotella stercoravium and Candidatus Blautia pullistercoris (the latter identified but not named by ). These findings rule out the concept of a core chicken gut microbiome. Studies on the human gut microbiome provide a useful comparison in that, in a recent study, only 14 genera were found to be shared across 95% of samples from the human gut.

54Among the species with high coverage, frequency is clearly linked to Bioproject. Although species quantification curves showed that the number of species identified increased rapidly with the number of samples, species discovery appeared to plateau at approximately 230 species after including only 50 metagenomes (Fig. S2A). Only two species appeared to be restricted (at ≥1× coverage) to just a single sample: Aliarcobacter thereius and Candidatus Avibacteroides faecavium. Correlation clustering confirmed structure in the data linked to BioProject (Fig. S2B)—for example, the BioProject from the study by clearly shows enhancement of clostridial species compared to other BioProjects, which reflects the fact that samples in that study were sourced from chicks with no post-hatching contact with an adult bird. However, the BioSamples do not appear to cluster by country (Fig. S2C) and show only limited clustering by anatomical/physiological sample site (Fig. S2D). Unfortunately, there is insufficient metadata for other potentially important factors, such as breed, age or diet to draw conclusions on how these might influence clustering.

55Given the dominance of chickens in the planetary biomass, the chicken gut microbiome ranks as one of the most abundant microbial communities on the planet. Here, we have exploited two complementary approaches—metagenomics and culture—to create an extensive catalogue of genes, genomes and isolates from this important ecosystem. Our work illustrates the value of combining culture-dependent and culture-independent approaches in analysing microbiomes.

56We have clearly demonstrated the advantages of shotgun metagenomic sequencing when applied to the chicken gut microbiome, providing catalogues of genes and genome sequences that takes us well beyond what can be achieved using 16S ribosomal RNA gene sequences. Similarly, the current study is much wider in scope than the previous study by, not just including analyses of viral genomes and cultured isolates, while also incorporating MAGs built from data not just from that study but from all publicly available metagenomic datasets. Furthermore, the limited overlap between bacterial species represented among our cultured isolates and in our MGS reinforces the utility of the combined approach. Nonetheless, the substantial co-linearity between genomes obtained by the two approaches—and with those from another similar metagenomic study—confirms the reliability of our binning approaches.

57We were surprised to find such a remarkable phylogenetic diversity within this commonplace livestock ecosystem—diversity that rivals that associated with the human gut. Our work has more than doubled the number of bacterial species known to reside in the chicken gut and has resulted in the creation of an unprecedented number of new Candidatus species. By including well-formed Latin binomials with the genomes we have uploaded into public repositories, we have ensured that the new proposed names and associated sequences will be integrated into commonly used online taxonomies and databases and will provide a stable taxonomic nomenclature for future studies. In addition, we have provided proof-of-principle for a scalable approach to Linnaean nomenclature that could be applied to species recovered from other metagenomic assembly projects.

58Given that we did not recover by culture some of the organisms that appear most abundant by metagenomics, there is clearly scope for additional culture-based investigations, using a wider range of cultural conditions—perhaps drawing on the precedent of the Human Microbiome Project to create and target a list of the ‘most-wanted-for-culture’ organisms documented by metagenomics. The fact that novel metagenomic species are still being recovered from human gut datasets that include tens of thousands of metagenomes—twinned with the promise of novel long-read and proximity-capture approaches to metagenome analyses—make it clear that our attempts here to analyse all currently available chicken gut metagenomes provide far from the last word on microbial diversity in this abundant and important ecosystem. Nonetheless, the availability of so many novel genes, genome and species represents a substantial step forward.

59The extensive catalogue of genes, genomes and isolates we have created here substantially improves the coverage of the chicken gut microbiome in the public databases and will make it possible to profile sequences from the chicken gut much more rapidly, easily and comprehensively, providing a valuable resource that lays the ground-work for future comparative and intervention studies. This study also sets a provocative precedent—relevant not just to animal microbiomes, but to studies on all microbiomes—assigning well-formed Latin binomials to hundreds of metagenomic species in a scalable alternative to the automated use of bland, unstable, user-unfriendly alphanumerical designations. Drawing on the precedent set by the current study, we have recently extended this approach to encompass creation of more than a million new names for Bacteria and Archaea. Thus, the time is now ripe to bring Linnaeus right into the heart of microbiome studies.