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

Sample collection and storage

Faecal 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.

Sixty 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.

Sequencing and subsequent workflow

Workflow 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 (Perez-Sepulveda et al., 2020) 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 (Connor et al., 2016). 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; Bolger, Lohse & Usadel, 2014). Metagenomic sequences for all samples have been uploaded to the Sequence Read Archive under Bioproject ID PRJNA543206.

Reference-based metagenomic analysis

An initial analysis of our chicken faecal sequences using the Kraken 2 taxonomic classifier (Wood, Lu & Langmead, 2019) 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.

The 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 (O’Leary et al., 2016) in January 2020. Bracken (Lu et al., 2017) 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 (Breitwieser & Salzberg, 2016).

Metagenomic assembly

We 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 (Glendinning et al., 2020; Sergeant et al., 2014; Foster-Nyarko et al., 2020; Luiken et al., 2020)

All 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 (Puente-Sánchez, Aguirre & Parro, 2016) 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 (Hildebrand et al., 2014). 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.

Sequence datasets from our fifty samples—together with 582 samples from the selected BioProjects—were assembled using MegaHIT (Li et al., 2016) 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.

Bacteriophage identification and characterisation

Contig 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 (Roux et al., 2015). 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 (Fu et al., 2012). 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) (Menzel, Ng & Krogh, 2016). Only bacteriophage genomes with BLASTN hit E-Value < 0.05, percentage identity >70% and query covering >50% were selected as reliable hits.

A 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 (Sullivan, Petty & Beatson, 2011). Escherichia bacteriophage coverage per sample was determined using Anvi’o v6.1 (Eren et al., 2015) 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).

Gene catalogue

Complete genes identified by Prodigal v2.6.1 (Hyatt et al., 2010) were clustered at 95% nucleotide identity using CD-HIT-Est v4.6.1 (Fu et al., 2012). Incomplete genes were then mapped to this complete gene list using Bowtie2 v 2.3.4.1 (Langmead & Salzberg, 2012) and any mapping at 95% nucleotide identity were incorporated into the relevant gene clusters. Finally, genes representing the 40 conserved marker genes defined by Mende et al. (2013) 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)

Abundance estimates of contigs and genes

Prodigal (Hyatt et al., 2010) 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 (Langmead & Salzberg, 2012) 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 (Li et al., 2009). Only reads with mapping quality >20, >95% nucleotide identity and >75% overall alignment length were retained. BEDTools v2.21.0 (Quinlan, 2014) 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 (Kang et al., 2019).

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

Binning

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

Species-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 Hildebrand et al. (2019). 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 (Parks et al., 2015), so that we could be confident that each bin represents a single species.

A second approach de-replicated all MAGs at 95% ANI (species-level) and 99% ANI (strain-level) using dRep Version 2.0 (Olm et al., 2017) 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.

Taxonomy of metagenomic species

We used the Genome Taxonomy Database Toolkit (GTDB-Tk Release 95) to perform taxonomic assignments on strain-level dereplicated MAGs (Chaumeil et al., 2019). In addition, genes from each MGS were analysed through GTDB-Tk (Release 95), proGenomes resource (Mende et al., 2017) 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 (R-Core-Team, 2018), before visualisation using ggplot2 (Wickham, 2016) 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 (Page et al., 2015). Comparison of our derived metagenomes with those of Glendinning et al. (2020) was performed at 95% ANI using dRep and visualised using web-tool BioVenn (Hulsen, De Vlieg & Alkema, 2008).

Bacterial culture

To estimate species richness and diversity, the Phyloseq package of R (R-Core-Team, 2018) was applied to the output from Bracken (Lu et al., 2017) 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⁻⁵. 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.

For dilutions 10⁻³–10⁵, 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.

Genome sequencing and analysis

Genomic 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).

Paired-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 (Bankevich et al., 2012). 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.

Genome sequence taxonomic assignment

Barrnap 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 (Altschul et al., 1990). 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 (Chun et al., 2018). We used ReferenceSeeker Version 1.6.2 (Schwengers et al., 2019) to determine ANI and conserved DNA values compared to RefSeq bacterial genomes (Completed March 2020) (O’Leary et al., 2016). 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 (Chaumeil et al., 2019). Genomes were clustered at 95% and 99% ANI before selection of a single representative isolate per species using dRep (Olm et al., 2017). 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).

Phylogenetic analysis

For phylogenetic analysis of all MGS and genome sequenced isolates we used PhyloPhlAn v3.0.58 (Asnicar et al., 2020) with the ‘diversity high’ and a proteome input predicted from all genome sequences using Prodigal v2.6.1 (Hyatt et al., 2010). Diamond v0.9.34 (Buchfink, Xie & Huson, 2015) was used to perform a search against 400 universal PhyloPhlAn markers. MAFFT v.7.271 (Katoh et al., 2002) was used to perform multiple sequence alignment before refinement with trimAl v.1.4 (Capella-Gutiérrez, Silla-Martínez & Gabaldón, 2009) and reconstruction into trees using FastTree v2.1 and RAxML v. 8.2.12 (Price, Dehal & Arkin, 2010; Stamatakis, 2014). All trees were visualised and annotated manually using the online iTOLv5.7 platform (Letunic & Bork, 2016). 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.

To 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 (Angiuoli & Salzberg, 2011), 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 (Stamatakis, 2014).

Reference-based profiling documents novel diversity

We 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).

We initially analysed the faecal samples using the k-mer-based programme Kraken 2, followed by refined phylogenetic analysis using the allied programme Bracken (Lu et al., 2017) (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.

When 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 (Hildebrand et al., 2019; Ravi et al., 2019), 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 (Narasinakuppe Krishnegowda et al., 2020) 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 (Almohaya et al., 2020)—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 (Kollarcikova et al., 2019).

Bracken 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 (Jamalludeen et al., 2008)—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 (Tsonos et al., 2012).

Although 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’.

The 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.

Metagenomic assembly uncovers a wealth of viral diversity

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

To 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 (Roux et al., 2015).

VirSorter 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).

Of 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).

Remarkable microbial genome diversity in the chicken gut

Next, 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).

Sequence 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) (Bowers et al., 2017). Genome sizes of the MAGs ranged from ~0.4 to 6.8 Mbp, while GC content ranged from 24% to 73%.

Then, 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 (Jain et al., 2018; Luo, Rodriguez-R & Konstantinidis, 2014). 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) (Hildebrand et al., 2019). Similarly, using canopy clustering (Nielsen et al., 2014), 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.

Analysis of bacterial metagenomic species, primarily using the Genome Taxonomy Database (GTDB) taxonomy (Parks et al., 2020), 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).

Of 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 Glendinning et al. (2020), with these sequences not currently available in public archives. However, only 158 of our metagenomic species possess validly published names based on Latin binomials.

We performed a search of PubMed with the species name and ‘chicken’, leaving aside the 33 species named by Glendinning et al. (2020). 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 (Yoon et al., 2003) and present in 197 chicken samples; Aliicoccus persicus, first isolated from a hypersaline lake (Amoozegar et al., 2014) and present in 241 chicken samples; and Bacteroides reticulotermitis, first isolated from the gut of a termite (Sakamoto & Ohkuma, 2013) and present in 39 chicken samples.

We 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.

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

Linnaean binomials for hundreds of new candidate species

Linnaeus first proposed the assignment of Latin binomials to provide a universal nomenclature for biological species (Linnaeus, 1759). The International Code of Nomenclature of Prokaryotes (ICNP) sets the rules for naming prokaryotic species (Parker, Tindall & Garrity, 2019), 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 (Parks et al., 2020) 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.

To 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 (Oren, 2017; Oren et al., 2020). 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 (Parks et al., 2020) 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).