Quantitative microbiome profiling reveals the developmental trajectory of the chicken gut microbiota and its connection to host metabolism

Abstract Revealing the assembly and succession of the chicken gut microbiota is critical for a better understanding of its role in chicken physiology and metabolism. However, few studies have examined dynamic changes of absolute chicken gut microbes using the quantitative microbiome profiling (QMP) method. Here, we revealed the developmental trajectory of the broiler chicken gut bacteriome and mycobiome by combining high‐throughput sequencing with a microbial load quantification assay. We showed that chicken gut microbiota abundance and diversity reached a plateau at 7 days posthatch (DPH), forming segment‐specific community types after 1 DPH. The bacteriome was more impacted by deterministic processes, and the mycobiome was more affected by stochastic processes. We also observed stage‐specific microbes in different gut segments, and three microbial occurrence patterns including “colonization,” “disappearance,” and “core” were defined. The microbial co‐occurrence networks were very different among gut segments, with more positive associations than negative associations. Furthermore, we provided links between the absolute changes in chicken gut microbiota and their serum metabolite variations. Time‐course untargeted metabolomics revealed six metabolite clusters with different changing patterns of abundance. The foregut microbiota had more connections with chicken serum metabolites, and the gut microbes were closely related to chicken lipid and amino acid metabolism. The present study provided a full landscape of chicken gut microbiota development in a quantitative manner, and the associations between gut microbes and chicken serum metabolites further highlight the impact of gut microbiota in chicken growth and development.


Animals and sample collection
Two hundred male chicks (Arbor Acres) that hatched on the same day were reared on the same poultry farm (Zhuozhou, China) with free access to water and corn-soybean-based diets.The diet formula of broilers was formulated according to the feeding standards of Chinese chickens (NY/T33-2004, Supplementary Table S1), and contained 3.0 kcal/kg metabolizable energy and 21.5% crude protein.
The broilers were divided into 10 flocks (20 birds per flock).The birds were kept at ambient temperature of 20°C, relative humidity of 50% and a lighting program of 16L/8D.The flocks were sampled eight times (1, 4, 7, 14, 21, 28 35 and 42 DPH) according to the study design, and one bird per flock was randomly selected each time.Whole blood of the selected birds was collected from the wing vein, and serum samples were obtained by centrifuging at 4000 × g for 15 min, which was then promptly frozen at -20°C for further analyses.The luminal content of four gut segments (the duodenum, jejunum, ileum, and cecum) was collected and frozen in liquid nitrogen.All of the samples were transported to the laboratory and stored at -80°C for further analyses.

DNA extraction and quantification
Microbial DNA extraction was performed according to the protocol described previously [1].Because mechanical lysis and repeated head beating could increase the DNA extraction efficiency [2], we performed three successive rounds of bead beating with glass beads (1 min, 30 s, and 30 s).Microbial genomic DNA from the luminal contents of each segment was extracted using the QIAamp DNA Stool Mini Kit (catalog no.51504, Qiagen, Hilden, Germany).The weight of the samples was recorded during the process of DNA extraction (Supplementary Table S2).
For microbial enumeration, we first constructed a 10-log-fold standard curve that ranged from 10 3 to 10 8 copies using appropriate reference organisms with known copy numbers of 16S or internal transcribed spacer (ITS) rRNA genes (Escherichia coli for total bacteria and archaea, Komagataella pastoris for fungi) as described previously [3].The genomes of E. coli and K. pastoris were extracted using a TIANamp Bacteria DNA Kit (DP302, Tiangen, Beijing, China) and a TIANamp Yeast DNA Kit (DP307, Tiangen, Beijing, China), respectively.The DNA concentration of the extracted microbial genomes was quantified by a Qubit â fluorometer (Invitrogen, Carlsbad, CA, USA).The number of bacteria and fungi was calculated per ng of DNA using an online calculator (https://cels.uri.edu/gsc/cndna.html),according to the following Euqation (1): Quantification of total bacteria, archaea, and fungi was carried out by quantitative PCR using an ABI 7500 Real-time PCR system (Applied Biosystems, Waltham, MA, USA) with a 20-µL reaction volume containing 10 µL TB Green â Premix Ex Taq TM (RR420A, Takara, Dalian, China).The primers used in this study were synthesized by Sango Biotech Co., Ltd.(Shanghai, China).The primers used for the quantification of total bacteria were 5'-TCCTACGGGAGGCAGCAGT-3' (forward) and 5'-GGACTACCAGGGTATCTAATCCTGTT-3' (reverse) [4].The primers used for the quantification of total fungi were 5'-CTTGGTCATTTAGAGGAAGTAA-3' (forward) and 5'-TCCTCCGCTTATTGATATGC-3' (reverse) [5].The absolute abundance of taxon A was calculated by Equation (2): Estimated absolute abundance of taxon A = Relative abundance of taxon A in bacteria/archaea/fungi × Total number of gene copies of bacteria/archaea/fungi The relative abundance of taxon A in bacteria or archaea was calculated from the following amplicon sequencing data of the V3-V4 regions of the 16S rRNA gene and from the ITS2 region in fungi.
ASVs were filtered out when they were present in fewer than two samples or had a feature count of less than 10.Taxonomy assignment was performed using a pretrained naïve Bayes classifier on the basis of the SILVA database (v132) [10] and UNITE database (v8) [11].For taxonomic annotation, all unassigned sequences and sequences annotated as mitochondria and chloroplasts were removed.
Samples were rarefied to the same number of reads for the downstream analyses.The Shannon index and principal coordinate analysis (PCoA) based on Bray-Curtis dissimilarity were calculated using the R package vegan [12].Statistical analyses for alpha diversity and beta diversity were performed using the Kruskal-Wallis test/Wilcoxon rank-sum test and permutational multivariate analysis of variance (PERMANOVA) [13], respectively.The symmetric Procrustes correlation coefficients between the microbiome based on relative abundance and absolute abundance and P values were obtained using the "protest" function in vegan.

Enterotype clustering and the stochasticity of the gut microbiota assembly
Dirichlet multinomial mixture (DMM) models were used to assign the samples to the community types [14], and bin samples based on the log-transformed absolute abundance.The appropriate number of clusters was determined based on the lowest Laplace approximation score.
The fast expectation-maximization microbial source tracking (FEAST, v1.0.0) algorithm was used to track the sources of microbial populations [15].FEAST is a highly efficient expectationmaximization-based method that estimates the fraction of a microbial community contributed by a potential source environment.Each sampling site was identified as a sink, starting with the jejunum, and all other segments were treated as potential sources.Only the proportion of microbes sourced from the adjacent anterior gut segment was considered in the present study.

Microbial cooccurrence network analysis
Microbial taxon-taxon cooccurrence networks were constructed using partial Spearman correlation using the R package ppcor (v1.1) in different segments.To reduce noise and false-positive predictions, network inclusion was restricted to genera that were present in at least 30% of samples.The false discovery rate (FDR) was used in each microbial network, and correlations with adjusted P values above 0.05 were filtered out [18].Topological features were estimated using the R package igraph (v1.2.11) [19].The edge number, vertices number, clustering coefficient, mean betweenness centrality, and average separation were calculated with functions embedded in the R package igraph.Edges present in only one subnetwork were specialist edges, and edges in more than one subnetwork were considered generalist edges.Keystone taxa are highly connected taxa that exert a considerable influence on microbiome structure and functioning [20].We identified the genera at the top degree from the four subnetworks as keystone taxa.

Metabolic profiling of chicken serum and data processing
We performed untargeted metabolomics to analyze the metabolites in chicken serum.The frozen serum samples were melted and mixed with 400 µL of a methanol/acetonitrile (1:1, v/v) mixture.The mixture was vortexed for 30 s, followed by ultrasonic extraction at 4°C for 10 min and incubation at -20°C for 60 min.After centrifugation at 12,000 × g for 15 min, the residue was removed, and the supernatant was collected.The supernatant was dried under nitrogen.The dried extractions were redissolved in 100 µL of a methanol/acetonitrile (1:1, v/v) mixture.The mixture was ultrasonically extracted at 4°C for 10 min.The supernatant was collected after centrifugation at 12,000 × g for 15 min.
To analyze the metabolites in chicken serum, an Agilent 1200 series high performance liquid chromatography (HPLC) system (Agilent Technologies, Santa Clara, CA, USA) equipped with an ACQUITY ultraperformance LC (UPLC) BEH C18 column (100 mm × 2.1 mm i.d., 1.7 µm, Waters, Milford, MA, USA) was used to separate the metabolites.For the positive mode, mobile phases A and B were water (containing 0.1% formic acid) and acetonitrile (containing 0.1% formic acid), respectively.
For the negative mode, mobile phases A and B were water (containing 5 mM ammonium formate) and acetonitrile (containing 5 mM ammonium formate), respectively.The following gradient was used: 0 min, 15% B; 2 min, 15% B; 10 min, 30% B; 30 min, 90% B; 35 min, 90% B; 36 min, 15% B; and 45 min, 15% B. The flow rate was 400 µL/min, the injection volume was 3 µL, and the column temperature was 45°C.An Agilent 6510 ESI-Q-TOF system from Agilent Technologies was used for mass spectral (MS) acquisition.The positive and negative ionization modes were used.Agilent Mass Hunter workstation software was used for data analyses and compound identification.The downstream analyses were based on the annotations generated from the METLIN [21] and HMDB [22] databases.
We filtered the features in two steps: 1) we retained the metabolic features present in > 30% of the samples; and 2) we filtered out the unknown metabolic features according to the annotations from the databases.A total of 1,109 different kinds of metabolites were used for downstream analyses.
The metabolites of serum at different time points were assessed using orthogonal partial least-squares discriminant analysis (OPLS-DA) in the R package ropls [23].The variable importance in projection (VIP) was calculated, which reflects the loading weight and the variability of the response explained by this component.We clustered the metabolites into six clusters using the R package pheatmap (v1.0.12).
Pathway analysis of metabolites in each cluster was performed using MetaboAnalyst 5.0 (http://www.metaboanalyst.ca)[24].The category and functional enrichment analysis of the metabolites were analyzed using the MetOrigin [25].

Microbe-metabolite interactions
To examine the association between the microbiota and the metabolome, we analyzed the symmetric Procrustes correlation coefficients between the microbiome based on relative/absolute abundance and metabolic profiling.P values were obtained using the "protest" function in vegan.The residuals of the correlations were generated from the results of the "procrustes" function from the R package vegan.
Microbiota from different segments and metabolome feature tables were analyzed using microbiomemetabolite vectors (mmvec) [26] to identify microbe-metabolite interactions based on their cooccurrence probabilities as predicted by neural networking.The top 1% microbe-metabolite interactions with the highest conditional probabilities were regarded as positive co-occurrences.

Statistical analysis and visualization
The correlations not mentioned above were calculated using Spearman's correlations using the R function cor.test.When multiple hypotheses were considered simultaneously, P values were adjusted to control the FDR [18].The intersecting sets were analyzed using the R package UpSetR (v1.4.0) [27].
The significantly changed metabolites during chicken growth were examined using general linear models (GLMs) on time series metabolic data.

Figure
Figure S1 The archaeal community in the chicken gut.(A) Prevalence of archaea in the gut four

Figure S2
Figure S2 Microbial community composition of the chicken gut microbiota at the phylum level.

Figure
Figure S3 Beta-diversity of the chicken gut microbiota based on quantitative microbiome

Figure S6 20 Figure
Figure S6 Longitudinal occurrence patterns of the bacterial genera in the chicken gut.A total of

Figure S8
Figure S8 Absolute abundance of the genera belonging to the "core" pattern.The changes in the

Figure S9
Figure S9 The four submicrobial networks generated from the genera in different gut segments.

Figure S10 24 Figure
Figure S10 Cooccurrence profiles of microbes and metabolites.Heatmap of the cooccurrence