Animal experiments
A total of 613 animals were sampled in the three experimental replicates (trials A, B, and C) performed in 2019 within the H2020 project HoloFood [24]; 205, 182 and 226 birds from A, B and C trials respectively. Broiler chickens from two genetic lines (Cobb 500 and Ross 308) and both sexes were reared in intensive farm conditions for 35-37 days. Each trial comprised 24 pens (12 groups replicated twice distributed in 3 treatments x 2 genetic lines x 2 sexes), each pen containing 40 animals. More details about the experimental design, diet and performance results are available in Tous et al. [25]. Chickens were euthanized, weighed, and sampled at days 7–8, 21–22 and 35–37 (multiple days were necessary due to workload, and these differences have been accounted for in the statistical analyses), hereafter simplified to three time points (days 7, 21 and 35). Caecal pathogens detection (Salmonella spp., Campylobacter spp., and Clostridium spp.) procedures using PCR are explained in detail in Tous et al. [25]. Molecular data was obtained from three different sections of the caecum. In short, the end of one of the caecum bags was isolated and longitudinally opened to gently collect ca. 100 mg of digesta for each metagenomic and metatranscriptomic analysis. After carefully washing the intestinal surface with saline solution, the mucosal layer was scraped and ca. 100 mg of mucosa were collected for host transcriptomic analyses. All types of samples were preserved in DNA/RNA Shield buffer (Zymo) and stored at -20 ºC until nucleic acid extractions.
Data generation
DNA and RNA extraction
A total of 613 metagenomic, 125 metatranscriptomic and 169 host transcriptomic datasets were generated. Both nucleic acids were extracted using a custom purification method optimised for samples preserved on DNA/RNA Shield buffer [45]. The protocol consisted in a bead-beating for tissue disruption, followed by digestion, nucleic acid separation (DNA and RNA) and purification steps. Samples were processed in batches of 90 samples, along with 6 extraction, library preparation, and library indexing blanks (2x2x2). Samples within each batch were randomised using a custom script, but different sample types were not mixed to minimise the risk of cross-contamination due to DNA concentration differences.
Library preparation of metagenomic DNA
The nucleic acids extracted were fragmented to obtain an average length of 400 bp using a Covaris LE220 ultrasonication device. A standard amount of 200 ng of DNA was used for the library preparation. We used the BEST [46] ligation-based library preparation protocol to prepare sequencing libraries. In order to evaluate the success of the libraries, we conducted quality controls using qPCR assays. The optimal number of cycles was estimated to achieve the desired DNA molarity while reducing clonality. Any libraries that exceeded 12 cycles were repeated for library preparation due to potential technical biases. Subsequently, libraries were indexed using unique dual tags, along with the necessary PCR cycles. Before final quality-checks were performed with a DNA Fragment Analyser (Agilent), bead purification was carried out. Libraries with expected fragment-size distributions and molarities were equimolarly pooled for sequencing. Libraries that showed too low molarities were re-indexed to achieve the desired molarity, and the ones exhibiting unusual fragment distributions and large adaptor dimers were re-built. Sequencing was performed in multiple BGIseq runs with 150bp paired-end chemistry. Sequencing effort per sample typically varied between 8GB and 16GB, equivalent to 26 and 52 million reads.
Library preparation of metatranscriptomic RNA
rRNA depletion was performed using TIANSeq rRNA Depletion Kit (Animal) (Cat.No. NR101-T1), and the remaining RNA was fragmented into 250-300bp, to finally reverse-transcribe into double stranded cDNA with random hexamers. Total cDNA was sent to Novogene for library preparation and sequencing. In short, cDNA libraries were constructed using Novogene NGS RNA Library Prep Set (PT042), which comprises end repair, A-tailing and adapter ligation steps. The libraries were checked with Qubit and real-time PCR for quantification and Bioanalyzer (Agilent) for size distribution detection. Quantified libraries were pooled and sequenced on an Illumina NovaSeq 6000 platform with 150bp paired-end chemistry, aiming for 5GB of protein-coding gene data.
Library preparation of chicken transcriptomic RNA
Total RNA was quantified using Nanodrop (Thermo Scientific) and Bioanalyzer 2100 (Agilent), as well as analysed for integrity (Agilent 2100) and purity (agarose electrophoresis and Nanodrop). Samples were subjected to rRNA removal step by poly-A enrichment, using poly-T oligo-attached magnetic beads. After fragmentation, the first strand cDNA was synthesised using random hexamer primers. During the second strand cDNA synthesis, dUTPs were replaced with dTTPs in the reaction buffer. Total cDNA was sent to Novogene for library preparation and sequencing. The directional libraries were ready after end repair, A-tailing, adapter ligation, size selection, USER enzyme digestion, amplification, and purification. Libraries were checked with Qubit and real-time PCR for quantification, as well as with bioanalyzer for size distribution detection. Quantified libraries were pooled and sequenced on NovaSeq 6000 (Illumina), according to effective library concentration and data amount required.
Bioinformatic data processing
Generation of the MAG catalogue
Details on the procedures employed to create the MAG catalogue used in this study will be published, and code can also be accessed at Workflowhub (https://workflowhub.eu/programmes/28). In short, data from 261 caecal metagenomic samples collected from chickens from the three experiments were used to generate the caecal MAG catalogue. We performed de novo metagenomic assemblies using the MGnify assembly pipeline [47]. The assembly tool MetaSPAdes [48] was used preferentially for single-run assemblies, while MEGAHIT [49] was used for co-assemblies when memory requirements for MetaSPAdes were too high. Samples prioritised for co-assembly were selected by hierarchical clustering based on Jaccard distance between low-quality bins generated by single assembly. Contigs shorter than 500 base pairs were excluded, and further host, human and PhiX decontamination was performed post-assembly with blastn [50]. Contig binning was performed using ‘binning’ and ‘bin_refinement’ modules of metaWRAP’s. Genome quality assessment was conducted using checkM [51], with the criteria of retaining genomes with completeness >50%, contamination <5%, and a quality score (QS) >50 (where QS = completeness - 5*contamination). Genomes were de-replicated using an Average Nucleotide Identity (ANI) of 95%, and 30% alignment fraction to generate species-level clusters using dRep [52]. Lastly, GUNC [53] was employed to identify potentially chimeric genomes for subsequent removal, utilising specific parameters that included a clade separation score >0.45, contamination >0.05, and reference representation score >0.5.
Functional annotation and distillation of MAG catalogue
Taxonomy annotation and phylogenetic tree construction was carried out using GTDB-Tk [54]. Functional annotation of the MAGs was performed through an ensemble approach implemented in DRAM [55]. This approach incorporates data from various databases, including Pfam [56], KEGG [57], UniProt [58], CAZY [59] and MEROPS [60]. To distil these annotations into quantitative genome-inferred functional traits (GIFTs) representing metabolic capacities provided by the microbiota to its host, we used the R package DistillR, which can be found at the following link: (https://github.com/anttonalberdi/distillR). DistillR contains a set of >300 metabolic curated metabolic pathways and modules derived from KEGG and MetaCyc [61] databases, which are used to obtain quantitative estimates of the metabolic capacities of microorganisms through quantifying the relative representation of genes required for accomplishing a metabolic task. GIFTs range between 0-1, the zero indicating none of the genes defined in the pathway are present in the genome and one indicating that all genes are present. In cases where a step within a pathway requires the presence of two Identifiers, the step is considered full if both Identifiers are present, half full if one is present, and empty if none is present. We quantified 170 GIFTs per genome (complete detailed list can be found in Supplementary Table S1), whose values were first corrected by MAG genome completeness to reduce functional biases [62], and then averaged to obtain a genome-level overall metabolic capacity metric, hereafter referred to as Metabolic Capacity Index (MCI). We also distilled microbial gene expression data into 170 GIFTs per genome, by weighing gene expression values according to the weight of each gene in each metabolic pathway.
Genome-Scale Metabolic Networks
A genome-scale metabolic network (GSMN) is a comprehensive representation of all the metabolic reactions that occur within an organism. It is constructed based on the genomic information of the organism and integrates biochemical and genetic knowledge to capture the complexity of the organism's metabolism. We employed the software metage2metabo [17], which in turns relies on Pathway Tools [63] to reconstruct the GSMNs of every bacterial genome in our study using a custom snakemake [64] pipeline. Shortly, due to software dependencies, bacterial genomes were re-annotated using eggnog-mapper2 [65] against the eggNOG 5.0 database [33]. The annotation files were transformed into Genbank annotation files (gbk) using ‘emapper2gbk’ and SBML files generated using ‘m2m recon’, as implemented in metage2metabo. SBML files were analysed using the package COBRApy [66] and custom python scripts to quantify source, transit and sink metabolites as explained below.
Metagenomic data processing and read mapping
Sequencing adapters and identical duplicates were filtered out using AdapterRemoval 2.2.4 [67] and seqkit 0.7.1 [68]. Sequences were mapped to the latest chicken reference genome (galGal6, NCBI Assembly accession GCF_000002315.6) using bwa [69] increasing the minimum seed length to 25 to minimise the likelihood of incorrect read pair alignments from the metagenomic fraction. To evaluate the quality of the alignment, mapping statistics including depth and breadth of coverage, and percentage of mapped reads were calculated using SAMtools 1.11 [70]. Aligned reads were sorted and the metagenomic fraction was isolated using SAMtools. Metagenomic reads were mapped to the MAG catalogue using bwa at 90% identity and 60% coverage threshold and further summarised with samtools. Read-mapping counts resulting in < 30% genome coverage per sample were removed from further analysis. Retained read mapping counts were divided by the total number of paired-reads per sample, and multiplied by 100 to give the percentage of reads mapped to the MAG catalogue for each sample. Relative abundance was estimated by adapting the RPKG (Reads Per Kilobase per Genome equivalent) formula provided by Nayfach and Pollard. It is referred to as RPMM (Reads Per Million bases of genome, per Million mapped reads), as reads mapped to MAGs were normalised both by genome length (divided by 1M) and by read length (divided by 1M).
Metatranscriptomic data processing and read mapping
We employed a custom snakemake pipeline for preprocessing metatranscriptomic data (https://github.com/anttonalberdi/holoflow/tree/EisenRa/workflows/metatranscriptomics). In short, reads were trimmed and quality controlled using fastp [71], keeping reads >60 bp and with Phred scores >20. Processed reads were then mapped against the host genome (galGal6) using STAR [72]. The unmapped reads were subsequently mapped to a combined database containing SILVA 16S rRNA SSU and LSU NR 99 [73], as well as 5SRNAdb [74] using Bowtie2 [75] with default parameters. Unmapped reads were then mapped to the MAG catalogue genes (outputted from DRAM; genes.fna.gz) using Bowtie2 with default parameters. Finally, gene read counts were calculated using CoverM (https://github.com/wwood/CoverM), requiring both pairs of reads to hit the gene (--proper-pairs-only flag).
Chicken transcriptomic data processing and read mapping
Raw transcriptomic reads were quality-filtered using fastp, mapped against the host reference genome using STAR, and gene count data extracted using the gene count option. Each sample yielded on average 12.5±3.2 million reads against the 24,131 genes annotated in the chicken reference genome (galGal6).
Data analysis
Metagenomic data analysis
Metagenomic counts were standardised by MAG length and sequencing depth. Dirichlet Multinomial Mixtures (DMM) [76] were utilised to profile and identify enterotypes in chicken microbial communities. Models were run by setting the maximum allowed number of community types to 5. A total of three runs were performed, one for each sampling day. At day 7, two community profiles were defined, where half of individuals from trial C formed an enterotype, and the rest of individuals from trials C grouped together with animals from trials A and B, forming another enterotype. At day 21, the model detected two clearly defined enterotypes, where in one of them all animals from trials A and B clustered together, and in the other enterotype all animals from trial C. At day 35, three enterotypes that were consistent with trials were detected. Thus, we defined microbial enterotypes as the distinct and standard enterotypes. The distinct enterotype comprised chickens from trial C that grouped separately from the rest at day 7, and the rest of the chickens from trials C at days 21 and 35. The rest of enterotypes were grouped together under the standard enterotype. Top community driver bacteria (i.e. MAGs with highest contribution to discriminate between enterotypes) were identified by selecting the 3% of MAGs with the highest posterior probabilities in DMM analysis.
To assess the temporal development of the composition of microbial communities across time, the MAG sequence count table was transformed using centred log-ratio (CLR) [77] and submitted to the constrained ordination from R package vegan [78]. MAG sequence count table was constrained by the factor trial (categorical variable with three levels: trials A. B and C), sampling time (categorical variable with three levels: days 7, 21 and 35) and their interaction. The significance of the constraining variables was tested using 999 permutations. To test the null hypothesis of no differences between enterotypes in microbial composition and function at day 7, PERMANOVAs were fitted through the ‘adonis2’ function of R package vegan. Euclidean distance matrices of CLR-transformed microbial abundances and functional profiles were included as responses in the PERMANOVAs and trial, enterotype, chicken age, sex, genetic line and treatment were included as explanatory variables. P-values were generated with 999 permutations.
A cladogram derived from the GTDB [79] tree constructed by GTDB-tk for taxonomic annotation was built with R package ggtree [80]. Tips of reference genomes were pruned using the ‘keep.tips’ function included in the R package ape [81]. Relative abundances of each bacterial genome of the catalogue for both enterotypes at each sampling time were illustrated with barplots after counts were standardised by MAG length and sequencing depth.The order of the bacteria was based on the tree obtained with GTDB.
To assess the drivers of the relative abundance of B. fragilis_A (the main indicator of the distinct enterotype at day 7) on trial C and day 7 we used linear mixed effect models as implemented in the R package lme4 [82]. P-values for the fixed effects were computed with the R package lmerTest [83]. CLR-transformed abundance of B. fragilis_A was used as response variable and trial, treatment, chicken age, sex and genetic line were included as fixed explanatory variables. To account for the fact that chickens were nested within pens we included a pen-level random intercept (1|pen). Then, we calculated the marginal and conditional R2 using the R package MuMIn [84]: marginal R2 captures the variance explained by fixed effects whereas the conditional R2 quantifies the variance explained by fixed and random effects together. Therefore, the variance associated with random effects (between-pen variance in relative abundance of B. fragilis_A) was calculated by subtracting the marginal from the conditional R2.
Genome-Scale Metabolic Network analysis
General statistics of metabolic properties were calculated for each bacterial genome using custom python functions. These included listing and quantifying source, transit and sink metabolites. Source metabolites were defined as metabolites that a given bacteria is able to use as reactant in at least one metabolic reaction inferred from the genomic information, but that the bacteria is unable to produce itself. Therefore, source metabolites have to be acquired from elsewhere, and can potentially be provided by other bacteria. Transit metabolites were defined as any metabolite that a bacteria can produce and utilise, while sink metabolites are a special case of transit metabolites, defined as metabolites that a bacteria can produce but is unable to use itself. Sinkmetabolites and, potentially transit metabolites, are therefore metabolic by-products that may become available for other bacteria.
Using the GSMNs of all 882 characterised bacteria, we calculated the community-weighted average number of source metabolites for Campylobacter that the microbial community associated with each chicken at day 7 could potentially produce. Additionally, we calculated the capacity of each enterotype to produce specific source metabolites for Campylobacter. We did so by calculating the cumulative relative abundance of bacteria in each sample collected at day 7, capable of producing each metabolite. First, to test the null hypothesis of no difference between enterotypes to produce source metabolites for Campylobacter,we used linear mixed models as implemented in the R package lme4. The capacity of each enterotype to produce source metabolites was used as response variable and enterotype, trial, treatment, chicken age, sex and genetic line were included as fixed explanatory variables. A pen-level random intercept (1|pen) was used to account for the nested design. The assumptions of homoscedasticity and normal distribution of errors were verified with visual inspection of residual plots. Then, the null hypothesis of no difference in capability of producing specific source metabolites for Campylobacter between enterotypes was tested using generalised linear mixed models using the function glmmPQL() of R package MASS [85]. In this case, the response variables were fractional (i.e. they take values between 0 and 1), thus we used a quasibinomial distribution with logit link function, which allowed us modelling the fractional response variables while accounting for under- or overdispersion and obtaining robust standard errors [86]. The set of fixed explanatory variables and random effects were the same as in above linear mixed models. Since multiple metabolites were tested consecutively, P-values were corrected for multiple testing using the Benjamini-Hochberg false discovery rate procedure [87].
Metatranscriptomic data analysis
Quantitative GIFTs were calculated with ‘distillq’ function using the R package distillR. To explore the temporal development of the functional expression profile of the standard and distinct chicken caecum enterotypes, the CLR-transformed community level quantitative GIFT profiles were ordinated with the constrained ordination RDA using ‘rda’ command from R package vegan. The ordination was constrained with the continuous variable chicken age, the categorical variable enterotype and their interaction. The significance of the factors was assessed using 999 permutations.
To identify specific quantitative GIFTs differentially expressed between standard and distinct chickens at different time points, linear mixed effect models were used with the R package lme4. CLR-transformed community level quantitative GIFTs were used as response variables in linear mixed models. As fixed explanatory variables in the models we used the enterotype, chicken age, sex, genetic line and treatment, and a pen-level random intercept (1|pen) was included to account for the nested design. P-values were corrected for multiple testing using the Benjamini-Hochberg false discovery rate procedure.
Lastly, to assess which bacterial strains were contributing the most to the expression of specific functions at different time points and enterotypes, we calculated the average relative expressions (given as gene counts per million) of specific functions for each MAG, at different combinations of enterotype and sampling time.
Chicken transcriptomic data analysis
Gene counts were processed with the tidybulk R package [88]. Briefly, gene counts were imported with the tidyverse metapackage [89]. Then, counts were normalised using the TMM method from edgeR [90]. Samples were clearly differentiated by sex and age. Thus, sex was a controlled variable for subsequent analyses. No statistical differentiation between breeds could be observed. We set as confounding variables the sex, breed, laboratory (two sample extraction batches) and sequencing batch (three sequencing batches). Then, samples were compared for differential expression between trials A and B versus C at the three sampling times. The method chosen for differential expression was the one implemented in edgeR. p-values were corrected for multiple testing using the Benjamini-Hochberg method [87]. Analysis of overrepresented Gene Ontologies [91] and KEGG pathways was done with clusterProfiler [92].