Metagenomic investigation of the equine faecal microbiome reveals extensive taxonomic diversity

Sample collection and storage

Faecal samples were from five, 12-month-old Thoroughbred racehorses from the same farm and field in Ireland. All samples were collected in April 2019 from horses raised on permanent pasture of mixed ryegrass. Horses were not being exercised at the time of sample collection. Feed supplementation whilst at pasture was proprietary post weaning cereal and trace element pellets plus an additional trace mineral and amino acid supplement. All horses had received ivermectin and praziquantel paste four weeks prior to sampling. Samples were collected as part of the Alborada Well Foal study, under the University of Surrey’s ethical review framework, project code: NERA-2017-007-SVM. 100 g of freshly evacuated faeces was collected from each horse in sterile tubes before immediate storage at 4 °C on site at the stud. All samples were shipped the same day at ambient temperature and received within 24 h. Upon receipt, samples were refrigerated before being aliquoted and stored at −80 °C until DNA extraction. Samples were thawed and homogenized before DNA extraction using the DNeasy PowerSoil kit (Qiagen), following manufacturer’s instructions. Extracted DNA was stored at −20 °C before further analysis.

Metagenomic sequencing and processing

Illumina sequencing libraries were constructed as previously described by Ravi et al. (2019). Paired-end metagenomic sequencing was performed on the Illumina NextSeq, before bioinformatic processing on the Cloud Infrastructure for Microbial Bioinformatics (CLIMB) (Connor et al., 2016). Output reads (2 × 150 bp) were assessed for quality using FastQC v0.11.8 and then trimmed using Trimmomatic v0.36 configured to a minimum read length of 40 (Andrews, 2019; Bolger, Lohse & Usadel, 2014). All metagenomic samples described here can be accessed on the Sequence Read Archive under BioProject ID PRJNA590977. Reads were aligned to the horse genome (GCF_002863925.1) using Bowtie2 v2.3.4.1 (Langmead & Salzberg, 2012), allowing removal of host reads with SAMtools v1.3.1 (Li et al., 2009).

Taxonomic profiling of sequencing reads was performed using Kraken 2 (Wood, Lu & Langmead, 2019) to search a microbial database built from archaeal, bacterial, fungal, protozoan, viral and univec_core sequences in Refseq in January 2021. Bracken was used to estimate taxon abundance from Kraken 2 profiles, accepting only those taxa with >1,000 assigned reads (Lu et al., 2017). Bracken-database files were generated using “bracken-build” on our microbial database and visualised using Pavian (Breitwieser & Salzberg, 2016).

Metagenomic assembly and binning

Host-depleted reads were assembled individually from each metagenomic sample with MegaHIT (Li et al., 2016), using kmer sizes 25, 43, 67, 87 & 101, before assessing the quality of resulting contiguous sequences (contigs) with anvi’o v7 (Eren et al., 2015). Filtered reads from each sample were mapped against the associated assembly to provide an estimate of contig abundance using Bowtie 2 (Langmead & Salzberg, 2012). Resulting Sequence alignment/map (SAM) files were converted to binary alignment/map (BAM) files before being sorted and indexed using SAMtools (Li et al., 2009). Contig coverage depth was translated from each BAM file, before separately binning contigs >1,000 bp with MaxBin v2.2.6 (Wu, Simmons & Singer, 2016) and CONCOCT v1.1.0 (Alneberg et al., 2014) and binning contigs >1,500 bp with MetaBAT 2 v2.12.1 (Kang et al., 2019).

DAS Tool was applied to the output from all three bin predictors, generating a catalogue of 196 bins from five samples (Sieber et al., 2018). All bins were profiled against the BAM file for their source metagenomic sample using the anvi’o ‘anvi-profile’ workflow (Eren et al., 2015). Using the ‘anvi-interactive’ tool, each bin was refined manually according to GC content, single copy core gene (SCG) taxonomy and coverage as well as detection statistics. CheckM v1.0.11 (Parks, Imelfort & Skennerton, 2015) was used for quality assessment of all bins using the lineage_wf function. Bins showing >50% completion and <10% contamination were assessed for quality score (defined as estimated genome completeness score minus five times estimated contamination score), a commonly used standard for defining acceptable bin quality (Parks et al., 2017). Bins with <70% completion and/or a quality score of <50 were categorised as low-quality metagenome-assembled genomes (MAGs) (n = 29); those with >70% completion, <10% contamination and quality score >50 were categorised as medium-quality MAGs (n = 68) and those with >90% completion, <5% contamination and quality score >50 were classified as high-quality MAGs (n = 55).

Taxonomic and phylogenetic profiling of MAGs

Medium- and high-quality MAGs from all five samples were de-replicated at 95% average nucleotide identity (ANI) with a default aligned fraction of >10% using dRep v2.0 (Olm et al., 2017), to create a non-redundant species catalogue. Clustering at 99% ANI was used to identify a non-redundant strain catalogue and select a representative MAG per strain. CompareM v0.1.1 (Oksanen et al., 2019) was used to assign Average Amino-acid Identity (AAI) values followed by AAI clustering at 60% to allow delineation at the genus level.

The Genome Taxonomy Database Toolkit (GTDB-Tk) v1.5.0 (Chaumeil et al., 2019), the Contig Annotation Tool (CAT/BAT) v5.2.3 (von Meijenfeldt et al., 2019) and ReferenceSeeker v1.4 (Schwengers et al., 2020) were used to perform taxonomic assignment of representative MAGs at strain-level compared to the ‘GTDB release 202’, ‘NCBI nr (2021-01-07)’ and ‘NCBI RefSeq release 201’ databases, respectively. Where taxonomic assignments differed between GTDB-Tk, CAT/BAT or ReferenceSeeker, GTDB-Tk assignments took precedence. Only when no species-level GTDB taxonomy was available did we adopt assignments according to CAT/BAT or ReferenceSeeker (6% of assignments). Phylogeny for our final de-replicated catalogue of MAGs was performed by aligning and concatenating a set of sixteen ribosomal protein sequences (ribosomal proteins L1, L2, L3, L4, L5, L6, L14, L16, L18, L22, L24, S3, S8, S10, S17 and S19), an approach previously used to reconstruct the tree of life (Hug et al., 2016). Ribosomal sequences were extracted using anvi’o before alignment using MUSCLE v3.8.155 (Edgar, 2004) and refinement using trimAl v1.4 (Capella-Gutiérrez, Silla-Martínez & Gabaldón, 2009). A maximum-likelihood tree was constructed using FastTree v2.1 (Price, Dehal & Arkin, 2010). All novel MAG species clusters were confirmed as monophyletic, drawing on all publicly available genomes from the genus to which they had been assigned by GTDB (with genomes retrieved from NCBI). Proteomes were predicted using Prodigal v2.6.1 (Hyatt et al., 2010) before comparison against 400 universal marker proteins using PhyloPhlAn v3.0.58 (Asnicar et al., 2020) in accordance with diamond v0.9.34 (Buchfink, Xie & Huson, 2015). Multiple sequence alignment and subsequent refinement was performed using MAFFT v7.271 (Katoh et al., 2002) and trimAl v1.4 (Capella-Gutiérrez, Silla-Martínez & Gabaldón, 2009; Stamatakis, 2014). All trees were subsequently visualised and manually annotated using iTol v5.7.

Abundance and metabolic profiling of MAGs

To estimate the proportion of reads within each BioSample represented by our final, de-replicated MAG catalogue, contigs from the non-redundant MAG catalogue were concatenated and filtered reads aligned back to this MAG database using Bowtie 2 (Langmead & Salzberg, 2012). Ordered BAM files were assessed using anvi’o (Parks, Imelfort & Skennerton, 2015) to calculate coverage statistics per-contig, allowing the calculation of mean coverage across each assembled genome according to methods available at: https://merenlab.org/data/2017_Delmont_et_al_HBDs/ and described by Delmont et al. (2018). Species accumulation and distribution analyses were conducted using the Vegan package in R (Oksanen et al., 2019) before visualisation using ggplot2 (Wickham, 2016).

Functional profiling of high- and medium-quality MAGs (n = 123) was performed using DRAM (Distilled and Refined Annotation of Metabolism) at a minimum contig length of 1,000 bp (Shaffer et al., 2020). Predicted amino-acid sequences identified by Prodigal in metagenome mode (Hyatt et al., 2010) were searched against KOfam, Pfam, and CAZy databases. tRNA and rRNA sequences were identified in MAGs using tRNAscan-SE (Chan & Lowe, 2019) and Barrnap v0.9, (Seemann, 2018) respectively.

Bacteriophage identification and characterisation

VirSorter v1.0.5 (Roux et al., 2015) was applied to all contigs >5 kb within each BioSample. Contig sequences classified by VirSorter as Category 1 (“most confident”) or Category 2 (“likely”) were considered for further analysis. Candidate bacteriophage sequences were assessed for completeness and contamination, using CheckV v0.7.0 (Roux, Páez-Espino & Chen, 2021), retaining only the sequences classified as “High-quality” (>90% completeness) or “complete”. These sequences were collated and de-replicated using rapid genome pairwise clustering at 95% ANI with an aligned fraction of ≥70% to generate a catalogue of bacteriophage genome sequences. For dereplication clustering, all-vs-all genome comparisons were performed using BLASTn before ANI based clustering using the ‘anicalc’ and ‘aniclust’ CheckV scripts sequentially.

Bacteriophage contigs from the catalogue were used as queries in a BLASTn search against the NCBI non-redundant nucleotide database (conducted on 21/12/2020) using an e-value of ≤1e−5. Only matches with a query cover >50% and percentage ID >70% were selected as being significant. Initial taxonomic classification of phage genomes at order and family level was performed using https://github.com/feargalr/Demovir against a viral subset of non-redundant TrEMBL database with an e-value of ≤1e−5. For each viral contig, individual coding sequences were predicted using Prodigal (Hyatt et al., 2010), before concatenation for input into vCONTACT2 v0.9.19 (Bin Jang et al., 2019) for construction of a gene-sharing network incorporating a de-replicated RefSeq database of reference prokaryotic virus genomes. The resulting network was visualised using Cystoscape v3.8.0 (Shannon et al., 2003).

Reference-based profiling documents microbial diversity

Whole genome sequencing of five faecal samples derived from 12-month-old Thoroughbred horses, each yielded >6 ng/µl DNA and collectively generated >280 million paired reads or >84 Gbp of sequence data. Reads derived from the horse genome accounted for <1% of reads from each sample (Table S1). We initially analysed reads using the k-mer-based program Kraken 2, followed by refined phylogenetic analysis via the allied program Bracken. Such analyses revealed extensive novelty and diversity in the equine faecal microbiome, with >59% of sequence reads in each sample classified by Kraken as “unassigned”, i.e. from unknown organisms. Assignable reads represented all three domains of life, as well as viruses, although bacteria predominated, accounting for >89% of assigned reads in any sample (Table S2).

Bacterial reads were predominantly assigned to the four phyla in the NCBI taxonomy most associated with animal gut microbiomes—Proteobacteria, Firmicutes, Bacteroidetes and Actinobacteria. However, the Kraken 2 profiles also provided evidence of over thirty additional bacterial phyla in this ecosystem. Many of these appear to be novel in the context of the horse gut, including Deinococcus-Thermus, Thermotogae and the Candidatus phylum Cloacimonetes (also called WWE1), which has been reported almost exclusively from anaerobic fermenters and the aqueous environment (Calusinska et al., 2018; Limam et al., 2014). However, as this phylum has recently been detected in soil fertilised with manure from dairy cattle, chickens and swine and has been implicated in anaerobic digestion of cellulose, it may play important similar roles in the vertebrate gut (Limam et al., 2014; Laconi et al., 2021). Reads assigned to eukaryotes provided evidence of budding yeasts and apicoplexan parasites in these samples.

Remarkably, two samples showed a very high relative abundance of reads assigned to the genus Acinetobacter (44% and 66% of classified reads), mirroring similar findings on two healthy horses in a previous study using 16S rRNA gene sequences (Costa et al., 2012). Bracken assigns these reads to an implausible sixty-two species of Acinetobacter, which is more likely to represent misassignment of reads rather than genuine diversity within this genus in this context.

Over a hundred newly named bacterial species

We generated almost 200 non-redundant bins from single-sample assemblies using three different approaches to binning. 123 bins represent medium- or high-quality metagenome-assembled genomes (MAGs), 96 with ≥15 amino acid tRNAs (Tables S3 and S4). Genome sizes ranged from ~0.5 to 3.8 Mbp, while GC content ranged from 31% to 60%. De-replication at 95% ANI clustered MAGs into 110 species clusters, spanning ten phyla (Fig. 1A). An average of 18% of the initial, host-depleted metagenomic reads per sample were represented within the final, dereplicated MAG catalogue. According to GTDB, around half (48%) of the MAG species clusters belonged to the Bacteroidota, while just over a third (35%) belonged to the Firmicutes (split by GTDB into Firmicutes, Firmicutes_A and Firmicutes_C). Only fourteen of the bacterial species from the horse gut had been previously defined and delineated: nine with validly published Latin binomials and five simply with alphanumerical designations assigned by GTDB (these are placeholder names assigned when no well-formed Latin name exists for the species) (Table S5).

Two of the species with validly published names, Ligilactobacillus hayakitensis (synonym Lactobacillus hayakitensis) (Morita et al., 2007) and Limosilactobacillus equigenerosi (synonym Lactobacillus equigenerosi) (Endo et al., 2008), have been previously cultured from the faeces of thoroughbred racehorses and are thought to be positively associated with equine intestinal health (Morita et al., 2009). Similarly, the species Streptococcus equinus was named in the early twentieth century after its association with horse dung and has been repeatedly isolated from this source (Andrewes & Horder, 1906; Smith & Shattock, 1962). Another of the named species found among our MAGs, Treponema succinifaciens, has been reported from the equine gut by 16S studies (Daly et al., 2001), but ours represents the first report of a genome from this species in this setting.

The recently named species Acinetobacter lanii (Zhu et al., 2021) has been isolated from the Tibetan wild ass Equus kiang, but our MAG represents the first report of an association between this species and the domesticated horse. Although the genus Phascolarctobacterium is known to inhabit the horse gut (Metcalf et al., 2017; Edwards et al., 2020), here we provide the first evidence of a specific link between the horse and the species P. succinatutens, previously found in human and pig faeces (Watanabe, Nagai & Morotomi, 2012). Our MAG catalogue provides the first report in the horse of the species Pseudomonas lundensis, first isolated from meat, but now recognised as an emerging pathogen of humans (Molin, Ternström & Ursing, 1986; Scales et al., 2018).

Among our MAG species clusters, ninety-six represent new candidate species within sixty-one bacterial genera previously delineated by GTDB. The majority of these novel species had <85% ANI to their closest known representative within GTDB databases (Fig. 1B). Sixty of these genera occur in the gut microbiota of at least one additional mammalian host species. Eleven of our species that could be assigned only to the level of family fell into ten clusters (delineated at 60% AAI) representing novel candidate genera from seven different families (Table S6). The archaeal genus Methanocorpusculum is thought to play a role in methane production in the equine gut (Murru et al., 2018). Here, we have delineated a novel species from this ecosystem: Candidatus Methanocorpusculum equi.

Building on our recent efforts with the chicken gut microbiome (Gilroy et al., 2021) and with the automated creation of well-formed Latin names, we have created Candidatus names (abbreviated as Ca.) for all the unnamed taxa revealed by our metagenomic analyses (Table 1). We also created Latin names for species and genera recognised by GTDB, but previously assigned only alphanumeric designations. For taxa found only in the horse, we created names that incorporated Greek or Latin roots for this host (e.g., Ca. Equimonas). However, if searches of the GTDB and NCBI databases suggested that genera had representatives in other gut microbiomes, we opted for names that specified gut or faeces as habitat (e.g., Ca. Limimonas).

A novel class within the Armatimonadetes

One of our MAGs—and the associated species cluster, which we have called Ca. Hippobium faecium—was assigned to the family of alphanumeric designation UBA5829. Assigned by GTDB to its own class, order and family, all members of this family belong to the recently named phylum Armatimonadetes (also called Armatimonadota; previously known as OP10) (Tamaki et al., 2011). Scrutiny of the NCBI database in August 2021 reveals that no genome assemblies linked to this phylum originate from the vertebrate gut, instead being metagenome-assembled genomes largely derived from bioreactors. Ca. Hippobium faecium was found at >1× coverage in two samples (SAMN13344080 & SAMN13344082), with relative abundance of this species across both samples being 94% and 4% respectively.

Distribution and metabolism

Our de-replicated high- and medium-quality MAGs account for 18% (±5%) of our host-depleted metagenomic reads. Distribution analysis identified 17 species present at ≥1× coverage in all samples (core MAGs represent 15% of our dereplicated MAG catalogue), spanning four bacterial phyla and the archaea (Fig. 2A and Table S7). No species were present at ≥10× coverage in all samples. While the majority of identified MAG species clusters had predominant relative abundance in only one sample, species including Ca. Methanocorpusculum equi, Acinetobacter lanii, Ca. Colimonas fimequi and Ca. Colisoma equi had more uniform distribution across all sample indicating a more central function in equine health. Species quantification shows a steady incline in the cumulative number of species identified when successively adding each of the five separate horse faecal samples (Fig. 2B). While species in the genus Acinetobacter represented seven of the ten most abundant species in our samples according to the Kraken2 analysis, only one Acinetobacter species was identified within out MAG catalogue. It seems likely that the presence of multiple closely related strains and species from this genus increases the likelihood of chimeric, unresolved bins. More generally, high abundance of a genus in the Kraken2/bracken analysis is not strongly linked to the recovery of medium or high quality MAG from that genus.

We created a catalogue of 228,125 genes from our medium- and high-quality MAGs. All 123 MAGs encoded known carbohydrate-active enzymes (CAZymes), with an average of 69 CAZymes per genome (Table S8). Most (>70%) MAG species clusters with a higher-than average repertoire of CAZymes belonged to the Bacteroidota. Of the ~8,500 CAZyme genes reported, most were associated with classes devoted to assembly (glycosyltransferases [GT] 29%) and breakdown (glycoside hydrolases [GH] 51%) of carbohydrate complexes, with far fewer from other groups of CAZymes; being the polysaccharide lyases (PL) and carbohydrate esterases (CE) alongside two further non-enzymatic groups being the carbohydrate-binding modules (CBM) and the auxiliary activities (AA). (Fig. 2C). Recovery of 93 classes of glycoside hydrolases from the equine gut mirrors similar enzymes in the sheep rumen linked to fibre degradation (He et al., 2019). Over half of our equine MAGs encode CAZymes with presumed involvement in degradation of hemi-cellulose (58%), cellulose (51%) or pectin or soluble fibre (>60%).

Many novel bacteriophage genomes

The program VirSorter classified 2,500 contigs as “highly likely” or “likely” to originate from bacteriophages (Table S9). Of these, 190 bacteriophage genomes were identified as “high-quality” (n = 181) or “complete” (n = 9) after de-replication (Fig. 3A). However, as none showed close identity to known viral sequences, they all represent novel bacteriophage species. Genome sizes ranged from 5 to 145 kb, including 42 genomes ranging from 5 to 15 kb in length. Using the viral taxonomy tool Demovir (https://github.com/feargalr/Demovir), we could assign 150 of these new phages to known viral families. An additional 29 could be assigned to taxonomically informative viral clusters, based on similarities between predicted proteins from our contigs and proteins from the viral component of the RefSeq94 database (Table S10). Just under half (n = 14) of these viral clusters contained at least one reference genome, thus expanding the known diversity of four viral families (Fig. 3B).

Almost all of our viral genomes represented tailed dsDNA phages from the order Caudovirales (Babenko et al., 2020) and could be sub-classified into the families Siphoviridae (73%), Podoviridae or the newly delineated Schitoviridae (Wittmann et al., 2020) (21%) and Myoviridae (6%). Seven genomes were assigned to ssDNA viruses from the family Microviridae, four of which cluster as part of the subfamily Gokushovirinae. Weak connections of three viral genomes to a viral cluster of Obolenskvirus, whose known members all infect Acinetobacter sp., likely indicates the presence of novel bacteriophage genomes predating on the prominent population of bacterial Acinetobacter within the equine hind-gut. Present within the viral cluster network but notably absent within our bacteriophage catalogue included the model Escherichia coli phages T4 (Tevenvirinae) and T7 (Studiervirinae) or the Mycobacterium infecting actinophages. We observed several novel viral clusters comprising only genomes assembled in this study, which could be classified as the first representatives of new horse hindgut-associated phage families. Based on the proteome comparisons (Fig. 3B), we predict at least three new families.

Over three quarters of the recovered phage genomes were found at >1 × coverage in just a single sample (Fig. 3C), with observed phage genomes ranging in coverage from 27 to 65. Similar inter-individual variation in phage abundance and diversity has been described within the human gut microbiome despite evidence of strong temporal stability within individuals (Ogilvie & Jones, 2015). Variation in phage composition between samples is probably driven by environmental and host-derived factors, although the balance of these influences is yet to be defined (Duerkop, 2018). Only one phage was found in all five samples, with coverage ranging from 1.9×–29× and forming a cluster with Lactococcus phage P087 of the family Siphoviridae. The small sample size makes it impractical at this stage to define a core virome for the horse using criteria applied to the human gut microbiome (Manrique et al., 2016).