Host-genotype-dependent cecal microbes are linked to breast muscle metabolites in Chinese chickens.

In chickens, the effect of host genetics on the gut microbiota is not fully understood, and the extent to which the heritable gut microbes affect chicken metabolism and physiology is still an open question. Here, we explored the interactions among chicken genetics, the cecal microbiota and metabolites in breast muscle from ten chicken breeds in China. We found that different chicken breeds displayed distinct cecal microbial community structures and functions, and 15 amplicon sequence variants (ASVs) were significantly associated with host genetics through different genetic loci, such as those related to the intestinal barrier function. We identified five heritable ASVs significantly associated with 53 chicken muscle metabolites, among which the Megamonas probably affected lipid metabolism through the production of propionate. Our study revealed that the chicken genetically associated cecal microbes may have the potential to affect the bird's physiology and metabolism.

Introduction

The gastrointestinal tracts of chickens are densely populated with complex microbial communities (bacteria, fungi, archaea, and viruses), which are dominated by bacteria (Wei et al., 2013). Similar to humans, the most abundant bacteria in the chicken gut belonged to four phyla: Firmicutes, Bacteroidetes, Actinobacteria, and Proteobacteria (Wang et al., 2017b). Over 600 different bacterial species are estimated to occur in the chicken cecum, which hosts the largest bacterial population (Torok et al., 2011).

The gut microbiota plays vital roles in nutrition absorption, immunity development, and disease resistance based on the production of vitamins, short-chain fatty acids (SCFAs), organic acids, and antimicrobial compounds (Shang et al., 2018). Gut bacteria together with the metabolites have also been demonstrated to contribute to skeletal muscle metabolism and function through the “gut-muscle axis” (Frampton et al., 2020; Grosicki et al., 2018). In chickens, the quantity and diversity of the gut microbiota are correlated with production performance (Bae et al., 2017; Siegerstetter et al., 2017). Manipulation of the gut microbiota by supplementing probiotics from Lactobacillus, Bacillus, Saccharomyces, and Citrobacter has positive effects on the live body weight and feed conversion ratio in both chickens and other animals (Adhikari and Kim, 2017). In addition, increasing evidence indicates the beneficial effects of probiotic bacteria in affecting chicken meat quality, a term that is used to describe meat traits including meat color, pH, texture, and so on (Maltin et al., 2003). For examples, supplemented with Enterococcus faecium, Streptococcus faecium, and Bacillus subtilis have all been evidenced to improve the chicken meat quality (Al-Shawi et al., 2020; Cramer et al., 2018; Ivanovic et al., 2012; Trabelsi et al., 2019; Wang et al., 2017b; Zheng et al., 2014). These results provide a basic knowledge that the gut bacteria have the potential to influence the chicken muscle physiology and metabolism.

The composition of the gut microbiota is highly variable among individuals and influenced by factors from both the host and the environment (Wang et al., 2016). Although diet is generally accepted as the leading factor in shaping the gut microbiota (Spor et al., 2011), increasing evidence from genome-wide association studies (GWAS) has revealed the important effect of host genetics on the gut microbiota. Wang et al. (2016) found that the vitamin D receptor could influence the overall microbial variation and individual taxa. In a cohort of 1,098 people, a total of 58 single nucleotide polymorphisms (SNPs) were found to be associated with the abundance of 33 fecal microbial taxa (Turpin et al., 2016). The effects of host genetics on the composition of gut microbiota have also been revealed in animals. The relative abundances of 39 fecal microbial taxa were reported to be associated with 42 quantitative trait loci in 27 different mouse genomic regions (Leamy et al., 2014). In fecal and cecum luminal samples from pig, heritability estimates (h2) greater than 0.15 were found for a total of 81 and 67 microbial taxa, respectively (Chen et al., 2018). In the rumen microbiota of beef cattle, 59 out of 174 microbial taxa were considered heritable and 19 bovine genomic SNPs were associated with 14 of these heritable taxa (Li et al., 2019a). The host genotype-determined heritable microbial features contribute to the variation of various host phenotypes, including the body mass index (Goodrich et al., 2014), visceral fat (Beaumont et al., 2016) and metabolic syndrome (Lim et al., 2017) in humans, feed efficiency in cattle (Li et al., 2019a), methane emissions, rumen and blood metabolites, and milk production efficiency in cows (Wallace et al., 2019).

In chickens, the diversity and abundance of gut microbiota in individuals with different genetic backgrounds have been investigated. Zhao et al. (2013) compared the composition of fecal microbiota in two lines (high body weight vs low body weight) from a common founder population and showed that 68 microbial species were affected by genotype (line), gender, and genotype by gender interactions. Ding et al. (2017) reported that 109 genera were significantly different among three native Chinese breeds reared under similar feeding conditions. Similar to studies in other animals (Chen et al., 2018; Fan et al., 2020; Leamy et al., 2014; Li et al., 2019a), these results confirmed the role of chicken host genetics in shaping their gut microbiota. However, research efforts exploring the microbiota-associated detailed genetic loci in chickens and the relationship between heritable microbes and chicken physiology or phenotype are at the nascent stage. In analogy to heritability, microbiability was introduced and defined as the extent to which the variations in phenotypes can be inferred from the changes in gut microbiota (Difford et al., 2018). Wen et al. (2019) suggested that the genera Methanobrevibacter and Mucispirillum are closely associated with fat deposition in chickens; however, they are largely independent of host genetics. Another study found a significant association between two chicken genetic loci and specific Methanobacterium species, and a locus within the doublesex and mab-3-related transcription factor (DMRT) gene cluster accounted for ∼21% of the variation in the chicken gut microbial community structure (Ji et al., 2019). A recent study demonstrated that the chicken genetic loci associated with the genes MTHFD1L and LARGE1 were related to the abundances of cecal Megasphaera and Parabacteroides, respectively (Wen et al., 2021). These findings provide clues for further investigations using integrated methods to explore the interaction among chicken genetics, gut bacteria and host physiology, and metabolism.

In this study, we performed an analysis of a cohort of 100 birds from ten breeds reared under the same husbandry and dietary regimes by integrating multiomics, including SNP-based whole genome genotyping, microbial amplicon sequencing, and metabolomics, to reveal the effect of genetics on the cecal microbiota and the contribution of heritable taxa to chicken muscle metabolites, which largely determine the muscle characteristics and meat quality (Muroya et al., 2020). Our study provides new insights into the interaction between chicken genetics and their cecal microbiota, which would be beneficial for future manipulation of both the host genome and the gut microbiota to improve chicken performance.

Results

Whole genome SNP-based phylogenetic relationships of 10 chicken breeds

We genotyped 94 chickens from seven local Chinese purebred lines and three crossbred lines by using the commercial chicken 55K SNP genotyping array [ten birds per breed, except Chahua (CH, n = 8), Green-eggshell (GE, n = 8), Tibetan (TC, n = 9), and Shouguang (SC, n = 9)] (Table S1). The average call rate for each breed ranged from 98.1% Nongda Yellow B (NYB) to 99.1% Subei Grass (SG) (Table S2). After quality control, a total of 39,063 SNPs were generated from these samples. The principal component analysis (PCA) indicated that individuals from the same breed were tightly clustered together, and different breeds were clearly separated from each other (Figure 1A, pairwise Wilcox test, p values <0.05 for PC1 values). Five chicken breeds, namely, TC, SG, Luxi Cockfighting (LC), SC, and CH, were more closely clustered, whereas the other five breeds were scattered. We then built a phylogenetic tree based on the identity by state (IBS) distance (Figure 1B), with CH considered the ancestral chicken breed (Zhang et al., 2019). The phylogenetic analysis confirmed the close genetic relationship between the CH breed and TC breed (Bao et al., 2008). The three crossbred lines displayed longer distances to the ancestor, indicating the role of artificial selection in shaping chicken genome. Collectively, these results suggested that the 10 chicken breeds used in this study had distant genetic backgrounds with less genetic variation within the same breed and more genetic differentiation among different breeds, which guaranteed the downstream analysis of the effects of host genetics on their gut microbiome diversity.

Figure 1

Phylogenetic relationship of ten chicken breeds based on whole-genome genotyping

(A) Principal component analysis of ten chicken breeds using SNP genotypes.

(B) Phylogenetic tree constructed based on the IBS distance generated from SNP genotypes. Individuals from the same breeds are collapsed. Individuals from different breeds are in different colors.

Different chicken breeds have distinct cecal microbial community structures and functions

We first investigated the community composition of the cecal microbiota in the 10 chicken breeds (Table S3). At the phylum level, Firmicutes and Bacteroidetes dominated the chicken cecal microbiota, accounting for ∼80% of the total microbial relative abundance (Figure S1A). The Firmicutes/Bacteroidetes (F/B) ratio, a potentially important biomarker for animal and human gut microbial composition (Mariat et al., 2009), differed in the 10 chicken breeds (Figure S1B), with the highest in NYB and the lowest in SG, which might suggest the different capacities of the cecal microbiota to extract energy from the diet in different chicken breeds. Increased F/B ratio could increase energy harvest by the fermentation of dietary fiber to SCFAs, which could be used as energy by the host (Turnbaugh et al., 2006). At the genus level, the top ten most abundant genera accounted for ∼60% of the cecal microbiota (Figure S1C). Bacteroides was the dominant genus in the cecal microbiota regardless of breed. Dirichlet multinomial mixture (DMM) models were used to group individuals based on bacterial community composition. Two distinct community types—type 1 and type 2—were generated, and they were driven by three common genera, namely, Bacteroides, the Rikenellaceae RC9 gut group, and a genus from Lachnospiraceae (Figures S1D and S1E), indicating that these genera were important for the structure of the cecal microbiota of chickens. All individuals in the Beijing Fatty (BF) breed belonged to type 1, whereas all chickens in SC and SG were community type 2, and the two types coexisted in the other seven breeds (Figure S1F). Breeds harboring a higher proportion of community type 1 were prone to having a higher F/B ratio, e.g., BF, NYB, CH, and GE (Figures S1B and S1F).

We then investigated the diversity and co-occurrence network of chicken cecal microbiota at the ASV level. The cecal microbiota displayed different Shannon indexes among breeds [false discovery rate (FDR) < 0.1, Figures 2A and S2A], and the structures of the microbiota from different breeds were separated from each other (PERMANOVA, FDR <0.1, Figures 2B, 2C, and S2B). The correlation between the Shannon index and the first Principal Component Analysis (PCoA) axis suggested that the structure of the bacterial community had a close connection with alpha diversity (Figure 2D). Of note, approximately three-quarters (15½00, 75.5%) of the ASVs with relative abundances over 0.1% were shared by all of the breeds (Figure 2E), suggesting a core role of these ASVs in different breeds under the same rearing conditions. However, the complexity of the bacterial co-occurrence network in different breeds is disparate. Bacteria in the CH breed had the most complex co-occurrence network compared with other breeds according to the index of degree (p = 5.53 × 10−6, Figures 2F and S2C). In addition, the bacterial networks from different breeds rarely shared the same keystone ASVs, although ASV144-Bacteroides was a keystone ASV for both Recessive White (RW) and NYB chickens. The most common keystone ASVs were Bacteroides, followed by Desulfovibrio and Rikenellaceae (Figure 2G).

Figure 2

Characterization of the cecal microbiota

(A) Shannon index of each breed.

(B) PCoA based on Bray-Curtis dissimilarity of all the individuals.

(C) Values of PCoA axis one for each breed.

(D) Correlation between the Shannon index and values of PCoA axis 1 (Spearman’s correlation, rho = 0.526, p = 1.88 × 10−8).

(E) UpSet graphs show the number of shared major ASVs among the breeds. The pie chart shows the ratio of “core” ASVs (present in all breeds) out of the total ASVs. Red: “core” ASVs; blue: other ASVs.

(F) Degree and betweenness indexes of the microbial co-occurrence networks in each breed. Dots are the median values of the two indexes. X axis: degree index (p = 5.53 × 10−8, Kruskal-Wallis test); Y axis: betweenness index (p = 0.304, Kruskal-Wallis test).

(G) Keystone ASVs inferred from the microbial co-occurrence networks in each breed.

We next explored the functions encoded by the chicken cecal microbiota. According to PICRUSt2, the microbial metabolic pathways in the ten breeds can be clustered into nine clusters (Figure 3A). Each cluster had a different category composition (Figure 3B). Cluster 1, cluster two, and cluster three had a higher ratio of nucleoside and nucleotide biosynthesis; cluster seven was mostly associated with SCFA production; cluster 1, cluster three and cluster nine contained more pathways of vitamin biosynthesis and approximately half of the pathways of cluster four and cluster nine were involved in amino acid biosynthesis. The ten breeds exhibited different abundance profiles of the pathway categories except for carbohydrate biosynthesis (p = 0.39, Kruskal-Wallis test, Figure 3C). Pathways of nucleoside and nucleotide degradation, vitamin biosynthesis, and nucleoside and nucleotide biosynthesis were more abundant in the SG breed; carboxylate degradation and SCFA production were enriched in the RW breed; phospholipid biosynthesis and amino acid biosynthesis were highly represented in the NYB breed, and cell wall biosynthesis and enzyme cofactor biosynthesis were abundant in the GE and SC breeds, respectively.

Figure 3

Microbial metabolic pathways encoded by the cecal microbiota

(A) Heatmap of the differential metabolic pathways in different breeds. Gaps in rows indicate the divisions of each cluster. Red: higher abundance; blue: lower abundance.

(B) Counts of metabolic categories in nine clusters.

(C) Relative abundance distributions of ten categories in ten chicken breeds. Vertical dashed lines are the boundaries of 90% confidence intervals. p values for the ten categories were 0.003, 8.91 × 10−5, 1.54 × 10−5, 5.32 × 10−4, 4.93 × 10−4, 4.17 × 10−4, 0.389, 1.34 × 10−4, 0.002 and 1.18 × 10−5, respectively (Kruskal-Wallis test).

Taken together, these results suggested that the cecal microbiota from different breeds exhibited different compositions and functions, even under the same rearing conditions and the same diet, indicating the potential role of host genetics in shaping the cecal microbiota.

Identification of 15 heritable ASVs and associated host genetic loci

To explore the effects of host genetics on the cecal microbiota composition, we analyzed the correlation between the genomic relationship matrices and the Bray-Curtis dissimilarity of the cecal microbiota. We found a significant correlation between the genotype and the structure of the cecal microbiota, with 8.1% of the total variation in the cecal microbiota explained by genetics (Mantel test, p < 0.0001, Figure 4A). The procrustes analysis showed that the correlation between genetics and microbiota was weaker for NYB than for the other breeds, which is probably because of the intrabreed variation of the core taxa composition in NYB chickens, a phenomenon has also been observed in horses (Massacci et al., 2020). We then fit a linear model to find heritable ASVs by using the first five host genome principal components. We found that 15 out of 200 ASVs were significantly associated with genetics (FDR <0.05), with R2 values ranging from ∼20% to ∼45% (Figure 4B). Notably, the ASV150-Rikenellaceae RC9 gut group and ASV58-Methanobrevibacter reached R2 values over 40%, suggesting that these two ASVs were highly associated with chicken genetics. The 15 heritable ASVs were from seven microbial genera, most of which (n = 8) were from Bacteroides and two were from Rikenellaceae.

Figure 4

Association between chicken genetics and the cecal microbiota

(A) Procrustes analysis of the correlation between the genomic relationship matrices and the microbial community based on Bray-Curtis dissimilarity. Circle nodes: cecal microbiota of each individual; triangle nodes: the genetics of each individual.

(B) Distribution of genetic principal component R2 values for different ASVs in the cecal microbiota. The Y axis shows the variance explained, and the X axis shows p values for each ASV. Only the names of ASVs with FDR <0.05 are shown.

(C) Manhattan plots of genome-wide association p values for the relative abundance of ASV38-Megamonas.

(D) Boxplot of the relative abundance of ASV38-Megamonas.

(E) Manhattan plots of genome-wide association p values for the relative abundance of ASV142-Bacteroides salanitronis DSM 18170.

(F) Boxplot of the relative abundance of ASV142-Bacteroides salanitronis DSM 18170.

(G) Manhattan plots of genome-wide association p values for the relative abundance of the ASV45-Rikenellaceae RC9 gut group.

(H) Boxplot of the relative abundance of the ASV45-Rikenellaceae RC9 gut group.

To identify the host genetic loci that were associated with the cecal microbiota, we employed GWAS to identify associated genetic loci for the 15 heritable ASVs. We found a total of 170 loci that were associated with the relative abundance of at least one heritable ASV (p < 8.53 × 10−8, Figures 4C, 4E, 4G, and S3) except for ASV36-Lachnoclostridium and ASV176-Bacteroides sp. SB5 (Figures S3C and S3K). Among the 170 loci, 35 were located within 34 known genes, and 135 were in intergenic regions (Tables S4 and S5). Notably, many of the ASV-associated genetic loci were located within or near genes related to tight junctions, cell proliferation and differentiation, and histological development. For example, we found that the ASV38-Megamonas-, ASV142-Bacteroides salanitronis DSM 18170- and ASV45-Rikennellaceae RC9 gut group-associated loci were located within CDH11 and CDH6 (Figures 4D, 4F, and 4H); ASV26-Bacteroides sp. Marseille-P3166 was associated with SNPs in ADAM17 and near JAM2 and PTPRO, and ASV9-Bacteroides sp. Marseille-P3166 with loci near GATA3. In addition, there were also ASV-associated genetic loci located within or close to genes having the potential to regulate gut microbiota, such as SLIT3, MMP16, VAV2, and NOX1.

Metabolic profiling of the chicken breast muscle

Bacteria have the capacity to change the physicochemical characteristics of poultry meat, such as improving the pH, color, and chemical composition, etc (Popova, 2017). Therefore, we first examined the relationship between the physical characteristics of chicken muscles and their cecal microbiota. Interestingly, we showed that the cecal microbiota exhibited a closer correlation with the physical characteristics of breast muscle, which are represented by the pH, color, and shear force (p = 0.023, Mantel test), than with the physical characteristics of thigh muscle (p = 0.059, Mantel test). In detail, both the alpha diversity (Shannon index) and the beta diversity (Bray-Curtis distance) of the cecal microbiota displayed a significant correlation with the color of breast muscle, whereas only the alpha diversity showed a correlation with the thigh pH but not with other physical characteristics (Figure 5A and Table S6). This finding prompted us to explore the influence of chicken cecal microbiota on the metabolites in their breast muscle.

Figure 5

Profiling the metabolites in chicken breast muscle

(A) Spearman’s correlations between the cecal microbiota and physical characteristics of breast muscle and thigh muscle. Shannon: Shannon index; BC distance: Bray-Curtis distance; Brown: positive correlation; green: negative correlation. The stars indicate the associations that are significant after adjusting for multiple testing (FDR <0.05).

(B) Number of metabolites annotated by the HMDB database and Metlin database.

(C) Abundance and adjusted p values of the metabolites among the breeds. Top panel: density of the abundance of the metabolites; right panel: density of the adjusted p values.

(D) Top panel: Counts of metabolites that are significantly different among breeds. Bottom panel: Adjusted p values of the metabolites annotated by the HMDB database.

(E) Median abundance of the seventeen metabolites in different breeds that may come from the gut microbiota (FDR <0.05).

We then profiled the metabolites in chicken breast muscle via UPLC-Q-TOF/MS. A total of 1755 compounds were identified according to the annotation of the Metlin database, and 544 of them could be further annotated by the Human Metabolome Database (HMDB, Figure 5B). After filtration by removing metabolites present in less than 50% of the samples, 848 were considered in further analysis, among which the abundance of 534 was significantly different among the breeds (FDR <0.05, Kruskal-Wallis test, Figure 5C). According to the HMDB, 272 metabolites were assigned to 13 metabolic superclasses, among which lipids and lipid-like molecules were the most prevalent (89 metabolites), followed by organic acids and derivatives (55 metabolites) (Figure 5D). A total of 178 metabolites showed different abundances among the ten breeds, the majority of which (64 metabolites) were lipids and lipid-like molecules (FDR <0.05, Kruskal-Wallis test). Interestingly, 17 metabolites with different abundances (FDR <0.05, Kruskal-Wallis test, Figure 5E) in the breast muscle could also be found in fecal samples according to the annotation in the HMDB. Four of the 17 metabolites belonged to dipeptides and three of them, i.e., tyrosyl-alanine, phenylalanyl-alanine, and isoleucyl-alanine, existed in all breeds and were highly represented in the CH, RW, and Nongda Yellow A (NYA) breeds. Aminocaproic acid, a derivative of the amino acid lysine (Markowska et al., 2021), exhibited the highest abundance among the 17 compounds in all chicken breeds and was also enriched in the CH, RW, and NYA breeds. This compound can only be found in the diet or the gut microbiota (Atkinson et al., 2013; Xu et al., 2006) but not formed naturally in chickens. These results suggested that the metabolic profiles of the breast muscles differed in different chicken breeds and that the gut microbiota might exert effects on the metabolites in breast muscles.

Microbiota together with metabolites in breast muscle increased the ability to discriminate chicken breeds

As both the cecal microbiota and the metabolites in breast muscle were different among the chicken breeds, we tried to combine the features of microbiota and metabolome to discriminate different breeds. We compared three random forest models using the abundance profiles of the cecal microbiota, metabolome, or both. In the model constructed based on the profiles of the cecal microbiota, using the top 20 important features (Figure 6A) led to a classification accuracy of 57% (Figure 6B). Among these 20 features, eight were the heritable ASVs mentioned before. When using metabolites to construct the model, 10 of the top 20 features were dipeptides or tripeptides, and the accuracy of the model increased to 63% (Figures 6C and 6D). When combining the features of both the microbiota and metabolome of breast muscle together, the accuracy of the classification reached 72% (Figures 6E and 6F). In this feature-combined model, the top 20 features consisted of five ASVs and 15 metabolites, among which three were heritable ASVs and eight were peptides. These results suggested that heritable ASVs and peptides contributed greatly to or even acted as biomarkers for the discrimination of chicken breeds.

Figure 6

Classification models based on random forest

(A) Mean decrease in the Gini coefficient for the ASV-based model.

(B) Confusion matrix of the classification results using the ASVs.

(C) Mean decrease in the Gini coefficient for the metabolite-based model.

(D) Confusion matrix of the classification results using the metabolites.

(E) Mean decrease in the Gini coefficient for the ASV-based and metabolite-based models.

(F) Confusion matrix of the classification results using the ASVs and the metabolites together.

Heritable cecal ASVs contribute to variations in chicken breast muscle metabolites

To further explore the relationship between breast muscle metabolites and the cecal microbiota, we constructed a co-occurrence network based on the abundance of the microbiota, muscle metabolites, and microbial metabolic pathways (Figure S4). The network contained 374 nodes, including 194 ASVs, 148 metabolites, and 32 microbial metabolic pathway categories. The top three ASVs with the greatest number of degrees were ASV163-Rikenellaceae RC9 gut group, ASV121-Bacteroidales, and ASV116-Bacteroides gallinaceum, whereas the top three metabolites were stearyl citrate, 5-L-glutamyl-taurine, and 3b,6a-dihydroxy-alpha-ionol 9-[apiosyl-(1->6)-glucoside]. Microbial metabolic pathways belonging to fatty acid biosynthesis, secondary metabolite degradation, and phospholipid biosynthesis were the three categories most frequently connected to ASVs or metabolites. These results highlighted that the cecal microbiota might exert an influence on muscle metabolites through different microbial metabolic pathways.

We then focused on heritable ASVs and applied a linear mixed model to quantify the contribution of heritable ASVs, i.e., microbiability (m) (Difford et al., 2018), and host genetics, i.e., heritability (h) (Visscher et al., 2008), to the abundance variations of the metabolites. Although the estimated microbiability for most of the metabolites was not high compared with the estimated heritability (median 0.15 vs 0.42), the m of 66 metabolites (ranging from 0.08 to 0.21) was even higher than the h (Figures 7A and 7B, Table S7). The top three metabolites that were more influenced by the heritable ASVs than the host genetics were (2S)-2-butanol O-[b-D-apiofuranosyl-(1->6)-b-D-glucopyranoside], dioscoretine, and valyl-proline (Figure 7B). These three compounds were considered to be derived from the diet according to the HMDB and the previous study (Iwu et al., 1990). These may suggest that the gut microbiota play an important role in the bioavailable of these compounds, considering the same diet of the birds.

Figure 7

Contributions of host genetics and heritable cecal ASVs to metabolites in breast muscle

(A) Estimated heritability (h) and estimated microbiability (m) of metabolites calculated by LMM. The color bar indicates the degree of variance in the metabolites explained.

(B) Estimated microbiability (m) was higher than the estimated heritability (h) for 66 metabolites.

(C) Co-occurrence network of the five heritable ASVs and metabolites. Brown: heritable ASVs; green: metabolites; red edges: positive correlations; blue edges: negative correlations. Background colors indicate the metabolites from different categories.

(D) Heatmap of the correlations between the abundance of Megamonas and the predicted microbial SCFA metabolic pathways in cecal microbiota.

Next, we investigated the relationship between heritable ASVs and specific metabolites in breast muscle (Figure 7C). Five heritable ASVs — ASV38-Megamonas, ASV171-Bacteroides, ASV45-Rikenellaceae RC9 gut group, ASV118-Bacteroides coprocola DSM 17136, and ASV26-Bacteroides sp. Marseille-P3166 — were found to have connections with 53 muscle metabolites (52 negative correlations vs one positive correlation). These metabolites mainly belong to lipids and lipid-like molecules, organic acids and derivatives (dipeptides), benzenoids, and organoheterocyclic compounds, which are frequently associated with meat quality and flavor. Besides, metabolites, such as beta-guanidinopropionic acid and diadenosine polyphosphates, having other physiological functions were also found correlated with the heritable ASVs. Strikingly, among the five ASVs, ASV38-Megamonas negatively correlated with 37 metabolites, 13 of which belonged to lipids and lipid-like molecules, suggesting its potential role in modulating muscle lipid metabolism. As members of Megamonas are known as SCFAs-producers (Polansky et al., 2016), we correlated the abundance of Megamonas with the predicted microbial SCFA metabolic pathways in cecal microbiota (Figure 7D). We showed that Megamonas was positively associated with propionate, acetate, and acetate/lactate fermentation pathways. The strongest correlation was found between Megamonas and propanoic acid synthesis pathway, i.e., L-glutamate degradation VIII (to propionate) (rho = 0.757, p = 8.44 × 10−20, Spearman’s correlation).

Discussion

The chicken gut harbors numerous microbes that influence host nutrient utilization, gut morphology, immunity, and physiology. Factors affecting the chicken gut microbiota include diet, sex, stress, medications, husbandry, genetics, etc (Falony et al., 2016). The gut microbiota can be remarkably different between animals kept under the same conditions, thus highlighting the importance of genetics in influencing the microbiota (Dabrowska and Witkiewicz, 2016; Kurilshikov et al., 2017). However, because of the complexity of the associated factors as well as the gut microbiota, consistent findings have not been obtained on the effects of genetics on the gut microbiota, especially in chickens (Ding et al., 2017; Wen et al., 2019; Zhao et al., 2013). In this study, to control for confounding factors, we maintained uniform diet, sex, and husbandry conditions for ten chicken breeds with different genetic backgrounds, which allowed us to better understand the genetic effects.

Our results showed that the cecal microbiota composition and the microbial interaction network were very different among the ten chicken breeds despite the animals being raised under the same environmental conditions with the same diet (Figure 2). The difference in cecal microbiota structure and interaction led to the difference in cecal microbiota function, which was reflected by the differentially enriched microbial metabolic pathways (Figure 3). The cecum of chickens is the most important site for fermentation, and the microbiota within this intestinal segment has the capacity to digest foods that cannot be digested in the foregut. The resulting products, such as SCFAs, vitamin B groups, and essential amino acids, have various physiological functions related to the health and performance of the host (Lan et al., 2005; Torrallardona et al., 2003; Yan et al., 2017). By using chickens with distinct genetic backgrounds, our study confirmed the role of chicken genetics in shaping their cecal microbiota structures and functions (Ding et al., 2017; Pandit et al., 2018; Zhao et al., 2013).

In chickens, GWAS are frequently used to identify SNPs and candidate genes associated with their quantitative traits. A large number of SNPs have been reported to be associated with chicken disease resistance (Raeesi et al., 2017), egg production (Liu et al., 2019a), growth (Pértille et al., 2017), fat deposition (Moreira et al., 2018), and meat quantitative traits (Felício et al., 2013). However, genetic loci associated with the microbial community structure in chickens have rarely been reported. Only loci within three genes, DMRT1, pleomorphic adenoma, and lck/yes-related novel tyrosine kinase, were identified to be associated with the microbial community (Ji et al., 2019, 2020). In the present study, we identified 35 loci within 34 known genes, and 135 in intergenic regions were closely associated with 15 ASVs (Tables S4 and S5). The related genes are involved in different functions from tight junctions, cell proliferation, and differentiation to histological development. CDH6, CDH8, and CDH11 belong to the cadherin superfamily, which plays pivotal roles in important morphogenetic and differentiation processes during development and in maintaining tissue integrity and homeostasis (Stemmler, 2008). JAM2, CLDND1, and FLOT2 are all junction protein encoding genes, which are closely related to the integrity of the mucosal barrier (Gadde et al., 2017; Ohnishi et al., 2017; Xu et al., 2020). ADAM17 also contributes to the intestinal barrier integrity and regulates intestinal inflammation by cleaving proinflammatory cytokines TNF-α and IL-6 (Hsiao et al., 2013; Wang et al., 2005). PTPRO inhibits the proliferation of intestinal epithelial cells (Zhao et al., 2020), and GATA3 plays a key role in the homeostasis, development, and function of ILC3 subsets (Zhong et al., 2016). As the mucosal barrier is the most important interface for the interaction between gut bacteria and host cells, it is reasonable that the chicken genes affect the gut microbiota by tuning the gut mucosal barrier function.

We also found genetic loci within or near genes SLIT3, MMP16, VAV2, and NOX1 that have been reported to be associated with gut microbes. Variants of gene SLIT3 are associated with the abundance of unclassified Clostridiaceae and Dermacoccus in human gut and nasal microbiomes, respectively (Goodrich et al., 2016; Igartua et al., 2017), and MMP16 correlates with Succinivibrionaceae family in cattle (Wang et al., 2017a). In addition to correlation studies, direct evidence regarding gene function has also been provided. VAV2 was found to be required for the invasion of Campylobacter jejuni and involved in the uptake of Yersinia and Chlamydia (Krause-Gruszczynska et al., 2011; Lane et al., 2008; McGee et al., 2003). NOX1-derived H2O2 plays a critical role in maintaining homeostasis of the gut microbiota by governing the growth of bacteria like segmented filamentous bacteria (SFB) (Kumar et al., 2016; Miller et al., 2020). These findings together with our results highlight the important roles of these genes in shaping gut microbiota and warrant further investigation of the mechanisms involved in the interaction between these host genes and gut microbes.

In this study, we observed that the ASV150-Rikenellaceae RC9 gut group and ASV58-Methanobrevibacter were the top two ASVs most associated with chicken genetics (Figure 4). Coincidentally, both Rikenellaceae and Methanobrevibacter have been found to be associated with human genetics (Goodrich et al., 2016; Turpin et al., 2016). Moreover, the abundance of Methanobacterium from the same family as Methanobrevibacter was found to be associated with chicken genetic loci (Ji et al., 2019). Other heritable bacteria we identified included ASVs from Bacteroides and Prevotella, which is consistent with a study in swine in which certain taxa from these two genera were heritable across different sampling points (Bergamaschi et al., 2020). Collectively, these findings indicate that the same heritable bacterial taxon can be found in different hosts, suggesting the existence of common interactions between these gut bacteria and host genetics, which deserves to be better known. However, inconsistent results regarding heritable bacteria have also been described. For example, in the chicken cecum, Wen et al. (2019) identified 21 operational taxonomic units mainly from Firmicutes but not Bacteroidetes, and showed that the fat deposition-correlated Methanobrevibacter is nonheritable. Thus, large-scale sampling and sequencing efforts and strict trial designs that exclude confounding factors are still needed to confirm these results.

We evaluated the contribution of the 15 identified heritable cecal ASVs to the abundance variations of metabolites in chicken breast muscle. Interestingly, we showed that the contributions of the 15 heritable ASVs were higher than those of the host genetics to 66 metabolites, with m values ranging from 0.08 to 0.21 (Figure 7). This contribution level is comparable to that in previous studies on the contribution of the gut microbiota community to animal phenotypes. For example, the contribution of microbes to chicken abdominal fat deposition ranges from 0 to 0.24 in different digestive tract segments (Wen et al., 2019); the microbiability of horse behavior traits is approximately 0.15 (Mach et al., 2020); the microbiability of feed-related traits in Japanese quail ranges from 0.09 to 0.27 (Vollmar et al., 2020), and the swine belly weight has a microbiability estimate of 0.29 (Khanal et al., 2021). Although these findings reinforce the important roles of gut microbiota, especially heritable microbes, we should stress that the contribution of gut bacteria to host phenotypes cannot be overinterpreted. As we showed, host genetics still play a decisive role in the variation of most of the metabolites (more than 75% of the total, with h values ranging from 0.08 to 0.68). This finding was also observed in previous studies, wherein the microbiability estimates were usually lower than the heritability estimates for specific host traits (Sanglard et al., 2020; Wen et al., 2019).

We found that five heritable cecal ASVs were significantly associated with 53 specific metabolites in breast muscle, including lipids, dipeptides, benzenoids, etc (Figure 7). The ASV38-Megamonas was most associated with lipids and lipid-like molecules. Megamonas has been reported to produce different SCFAs, including propionic acid, acetic acid, lactic and succinic acids, but mainly propionic acid as fermentation ends (Polansky et al., 2016; Sakon et al., 2008). The SCFA-producing ability in Megamonas was confirmed by our finding that the bacterial SCFA fermentation metabolic pathways, especially propionic acid-generating pathways, were positively correlated with Megamonas (Figure 7D). We inferred that the negative correlation between the ASV38-Megamonas and lipid compounds was probably because of the propionic acid it produced. This finding can be supported by previous studies concerning the role of propionate in lipid metabolism. In HFD-fed mice, oral sodium propionate supplementation reduced inguinal white adipose tissues by regulating lipogenesis gene expression (Song et al., 2019). The infusion of propionate sodium through a cecal fistula in a pig model decreased serum and liver triglyceride levels by reducing lipid metabolism-related metabolites (Yu et al., 2019). By using isolated rat and human hepatocytes, propionate was verified to inhibit fatty acid de novo synthesis and thus affect lipid metabolism (Lin et al., 1995; Nishina and Freedland, 1990). These results together with our findings suggest that the genus Megamonas may have the ability to reduce host fat deposition. However, as lipid metabolism is important for normal muscle function (Cortright et al., 1997), whether there is a possibility that Megamonas expands in birds with poor metabolism of muscle needs to be further addressed.

The heritable ASV-linked breast muscle metabolites we identified here may play multiple roles in meat quality or flavor. Intramuscular lipids are known as important flavor-forming precursors; they not only act as solvents for volatile compounds generated during processing but also produce distinct flavors, such as volatile carbonyls, when degraded or reacted with other meat components (Khan et al., 2015). Peptides, such as dipeptides and tripeptides, are involved in modulating the flavor and physiological functions of meat by the taste of the amino acids they contain or by their antioxidative and pH-buffering capacity (Flores et al., 2000; Liu et al., 2003; Ramalingam et al., 2019; Wu and Shiau, 2002). We also showed that peptides could be essential markers in the random forest model to discriminate chicken breeds (Figure 6). Benzenoids are important organic molecules related to aromatic hydrocarbons that are involved in fragrance (Kirk and Othmer, 1983). In addition to flavor components, compounds modulating muscle physiology were also found related to heritable ASVs. Beta-guanidinopropionic acid was evidenced to regulate energy metabolism in skeletal muscle (Oudman et al., 2013). Diadenosine polyphosphates are found to be extracellular signaling molecules and cellular function regulators (Pliyev et al., 2014). Taken together, we propose that the chicken genetically associated cecal microbiota has the potential to influence the host; moreover, our findings support the practice of manipulating chicken gut microbiota to improve meat physical and sensory characteristics (Popova, 2017).

Conclusions

We revealed the potential effects of genetics on chicken cecal microbiota, showed that chicken genetics may be involved in the assembly of gut microbiota via different genes, and revealed that heritable cecal microbes could exert an influence on the metabolites of breast muscle and thus affect meat quality. We reinforce the notion that host genetics are the primary determinants for phenotypes in animals, although the gut microbiota can act as an additional factor because of its interaction with the host (Lu et al., 2018).

Limitations of the study

The main findings in this study are summarized in Figure S5. Although the multiomics data in this study permitted the investigation of the interaction between chicken genetics and the cecal microbiota, the sample size is still limited. Moreover, we omitted the effect of sex and only targeted the cecal microbial community, thus excluding the communities in the whole gastrointestinal tract. These limitations should be considered in future work. In addition, we should stress that the interaction among the chicken genetics, the cecal microbiota, and the metabolites of breast muscle we revealed here are mainly based on multi-omic correlation analysis. Future in vitro and in vivo experiments are highly needed to confirm these results and to reveal the underlying mechanisms, such as using pure culture of microbial species to explore the metabolic changes induced by Megamonas in chicken models.

STAR★Methods

Key resources table

Resource availability

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Yongfei Hu (huyongfei@cau.edu.cn).

Materials availability

This study did not generate new unique reagents.

Method details

Animal experiments and sample collection

A total of 100 60-week-old chickens (female) were randomly selected from 10 breeds (ten birds per breed, Table S1), including seven purebred lines, i.e., BF, CH, GE, LC, SC, SG and TC; and three crossbred lines, i.e., NYA, NYB and RW. All birds were raised with the same corn-soybean-based diet in the Poultry Genetic Resource and Breeding Experimental Unit of China Agricultural University. Feed was purchased from Charoen Pokphand Group Co., Ltd. (Ningbo, China), and its ingredients were shown in Table S8. The birds were kept at ambient temperature of 20°C, relative humidity of 50% and a lighting program of 16L/8D. To avoid the bias caused by the cage effect, all birds were randomly sampled from different cages. Blood samples were collected from the wing vein and stored at −20°C until DNA extraction. The birds were then euthanized by cervical dislocation. After slaughter, cecal contents were collected and immediately frozen in liquid nitrogen and then stored at −80°C. The whole breast and two legs were removed from the chicken carcasses, with the left breast and leg subjected to physical characteristic determination and the right breast and leg used for untargeted metabolomics analysis.

Evaluation of the physical characteristics of the breast and thigh muscles

The physical characteristics of the breast muscles and thigh muscles were determined according to previously described methods (Wan et al., 2017). Briefly, the pH of breast and thigh muscle (average value of three measurements) was measured at 24 h postmortem using a pH meter (pH meter FE28, METTLER TOLEDO, Shanghai, China). Meat color attributes, including lightness (L∗), redness (a∗), and yellowness (b∗), were measured using a CHN Spec Colorimeter CS-10 (CHNSpec, Zhejiang, China). The average value of 3 measurements, with one from the middle and two from the corners of the muscle samples, were recorded. To determine the shear force of the muscle, the meat sample was divided into three strips with dimensions of 3 cm (length) × 1 cm (width) × 1 cm (thickness) and then heated in a water bath at 72°C until the central temperature reached 70°C. The samples were cooled to room temperature and stored at 4°C for 12 h. The texture was treated using a texture analyzer (TMS-PRO, VA, USA) with a Warner-Bratzler cutting blade head to measure the shear force (average value of the three strips).

Genotyping of chickens

Host genomic DNA was isolated from the plasma of all chickens using the TIANamp Blood DNA Kit (Tiangen, China). Six DNA samples that did not pass the quality control before genotyping were not included in this analysis. The chickens were genotyped with the commercial chicken 55K SNP genotyping array (Liu et al., 2019a2019b). This 55K SNP chip contains 52,184 SNP probes across 28 autosomes and a sex chromosome (chrZ). After converting the genome coordinates to the chicken reference genome (galGal5), 28 autosomes were extracted for further analyses. During the process of genotype quality control, SNPs with minor allele frequencies (MAFs) smaller than 5% and genotyping call rates smaller than 97% were removed. Finally, a total of 39,063 SNPs were generated from the 55K SNP chip. The population structure was calculated by PLINK 1.90 (Chang et al., 2015) and visualized by PCA in R software. A neighbor-joining tree was constructed from the IBS distance to infer the possible evolution of the breeds using MEGA X (Kumar et al., 2018). The CH breed was used as an ancient root according to previously reported findings (Zhang et al., 2019). The phylogenetic tree was edited and exported using the iTOL tool (Letunic and Bork, 2019). To make the conclusion more robust, we performed pairwise Wilcox test (function pairwise.wilcox.test package stats in R) for the PC1 values and the breeds.

Sequencing and processing of bacterial 16S rRNA sequences from the cecum

Bacterial genomic DNA from cecal contents was extracted using the QIAamp DNA Stool Mini Kit (Qiagen, Germany). For all samples, the V3-V4 regions of the 16S rRNA gene were amplified using the primers Pro341F (5′-CCTACGGGNBGCASCAG-3′) and Pro805R (5′-GACTACNVGGGTATCTAATCC-3′) as suggested (Takahashi et al., 2014). The pooled library was sequenced on the Illumina HiSeq platform (2 × 250 bp). Raw fastq files were quality-filtered and taxonomically analyzed using QIIME2 (V2019.7) (Mueller et al., 2016). In brief, primers of imported sequences were removed via Cutadapt (V2.4) (Martin, 2011), and then DADA2 was used to filter and denoise sequences, remove chimeras, infer ASVs and generate the abundance table (Callahan et al., 2016). The ASVs were taxonomically annotated using a pretrained naïve Bayes classifier on the basis of the SILVA V132 database (Quast et al., 2012). 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 analysis. Inferring functional profiles of microbiota profiling were performed using PICRUSt2 implemented in QIIME2 (Douglas et al., 2020). Microbial metabolic pathways were annotated using the MetaCyc database (Caspi et al., 2020).

DMM models were applied to assign the samples to the community types (Holmes et al., 2012). The Shannon index and PCoA based on Bray-Curtis dissimilarity were calculated using the R package Vegan. The statistical analysis for alpha diversity and beta diversity was performed by using the Kruskal-Wallis test (pairwise) and PERMANOVA (Zapala and Schork, 2006), respectively. The intersecting sets were analyzed using the R package UpSetR (Conway et al., 2017). We also applied SParse InversE Covariance Estimation for Ecological ASsociation Inference (SPIEC-EASI), a statistical method for the inference of microbial ecological networks from amplicon sequencing datasets for each breed (Kurtz et al., 2015). The keystone ASVs were inferred from the number of edges of each node according to the theory by Banerjee and colleagues (Banerjee et al., 2018). The degree and betweenness statistics were calculated using the R package igraph (https://github.com/igraph/rigraph), and the Kruskal-Wallis test was applied to compare the differences among breeds.

Liquid chromatography-mass spectrometry (LC-MS)-based untargeted metabolomics analysis and data processing

Untargeted metabolomics was performed to analyze the metabolites in the breast muscle. Frozen muscle samples (100 mg) were ground with liquid nitrogen, and 800 μL of a methanol/acetonitrile (1:1, v/v) mixture was added, vortex-mixed and ultrasonically agitated for 10 min. The mixture was incubated at −80°C overnight and then centrifuged at 15,000 rpm for 15 min at 4°C. The supernatant containing metabolites was then concentrated using a vacuum centrifuge (Jiaimu. CV200, China). Finally, 100 μL of methanol/acetonitrile (1:1, v/v) was added to the concentrated extracts and ultrasonically incubated for 10 min. After ultrasonic treatment, the supernatant was centrifuged at 15,000 rpm for 15 min. The supernatant was then injected into the LC-MS/MS system for metabolomic analysis. The quality control (QC) samples, including the muscle tissue samples, were inserted throughout the entire experiment to ensure the stability of the system.

For metabolite analysis, an Agilent 1290 Infinity II UPLC system (Agilent Technologies, USA) equipped with an Agilent Eclipse Plus C18 column (2.1 mm × 100 mm, 1.8 μm) was used to separate the metabolites based on a previously published protocol (Li et al., 2019b). Mobile phases A and B were water and acetonitrile, respectively. The gradient was as follows: 0-2 min, 5% B; 2-20 min, 5-100% B; 20-25 min, 100% B. The flow rate was 300 μL/min, the infection volume was 2 μL, and the column temperature was 40°C. An Agilent 6545 ESI-Q-TOF high-resolution accurate-mass spectrometer was applied for MS acquisition. Both positive and negative ionization modes were performed. Mass spectrometry parameters were set as follows: 3.5 kV, collision energy; 325°C, gas temperature; 11 L/min, drying gas flow rate; 35 psi, nebulizer pressure; and 100-1700 m/z, mass range. Nitrogen was utilized for drying, nebulization, and collision. Agilent Mass Hunter Professional software was applied for data analysis and compound identification. The downstream analysis was based on the annotation from the Metlin (Guijas et al., 2018) and the HMDB (Wishart et al., 2018).

Integrated multiomics analysis

Association between genetics and the microbiota. To explore the correlation between genetics and the microbiota, we first built genomic relationship matrices based on alleles at molecular genetic markers based on the IBS distance (Meuwissen et al., 2014). The correlation between the genomic relationship matrices and the Bray-Curtis dissimilarity of the cecal microbiota was then calculated by the Mantel test in R software. We performed a procrustes analysis on the Bray-Curtis dissimilarity of the cecal microbiota and the genomic relationship matrices using the “Procrustes” function from the R package Vegan. Procrustes analysis (Grower, 1975) is a commonly used tool for analyzing the congruence of two-dimensional shapes produced from superimposition of principal component analyses from two datasets. To associate genomic PCA with the microbiota, we fit a linear model to each of the microbial features using the first five principal components (Kolde et al., 2018). Based on the residuals, we calculated the amount of variance explained (R2). GWAS of the major ASVs (n = 200) were performed with GEMMA (Zhou and Stephens, 2012), a genome-wide efficient mixed model association algorithm. The Manhattan plots were drawn by using the R package CMplot (Yin et al., 2021). The cutoff was set according to the total number of multiple hypothesis tests. In the present study, the cutoff was set to 8.53 × 10−8 (0.05/(39,063 × 15)).

Association between the microbiota and physical characteristics. The relationship between the cecal microbiota and physical characteristics of muscles was determined by calculating Spearman’s correlations. The values of the first PCoA axis based on Bray-Curtis dissimilarity were used as representatives of the structure of the cecal microbiota. The Shannon index values for each individual were used as representatives of the α-diversity.

Classification based on the microbiota and metabolites. To distinguish different chicken breeds according to the abundance of the microbiota and the metabolites, we applied the R package randomForest. Feature selection was performed by cross-validation using the rfcv function in the R package randomForest. A confusion matrix was used to evaluate the classifier’s predictions.

Association among microbiota, metabolites and microbial metabolic pathways. We first calculated the pairwise Pearson’s correlation among multiomics measures and constructed a cooccurrence network. The FDR (Benjamini and Hochberg, 1995) was applied in the R package psych after normalization. The correlations between each pair of variables were considered only when FDR <0.2. We then calculated Spearman’s correlation between the heritable ASVs and the metabolites and between the microbial metabolic pathways and the heritable ASVs, and a strict FDR value of less than 0.05 was considered. Network visualization was performed using Cytoscape (V3.6.1) (Shannon et at., 2003).

To evaluate the contribution of the heritable ASVs to muscle metabolite variations, we used a linear mixed model (LMM) in the R package lme4 (Bates et al., 2015). The log-transformed abundance of associated metabolites was fit according to (Equation 1):where y is the log-transformed abundance of the metabolites in breast muscle, β is the vector of fixed effects with the corresponding design matrix X (abundance of the heritable ASVs), μ is the vector of random effects (the first five host genetic PCs), is the residual term, and X and Z are the corresponding design matrices. The contributions of genetics (h) and heritable ASVs (m) were estimated with the r.squaredGLMM function from the R package MuMIn. The estimated heritability (h) was calculated according to (Equation 2):

The estimated microbiability (m) of heritable ASVs was calculated (Difford et al., 2016) according to (Equation 3):where are the heritable gut microbial variance, host genetic variance and residual variance, respectively.

Quantification and statistical analysis

All details of quantification methods, statistical analysis and significance calculations are discussed in each figure legend.

REAGENT or RESOURCE	SOURCE	IDENTIFIER
Biological Samples
Cecal content	This study	NA
Blood	This study	NA
Muscle	This study	NA
Critical Commercial Assays
TIANamp Blood DNA Kit	Tiangen	Cat#DP348-02
Chicken 55K SNP genotyping array	Liu et al. (2019b)	NA
Deposited Data
Amplicon sequencing data	This paper	NCBI SRA: PRJNA693673
Microarray source dasta	This paper	figshare: 10.6084/m9.figshare.16882417.v1
Original code	This paper	figshare: 10.6084/m9.figshare.19424789.v1
Oligonucleotides
Pro341F: 5′-CCTACGGGNBGCASCAG-3′	Takahashi et al. (2014)	NA
Pro805R: 5′-GACTACNVGGGTATCTAATCC-3′	Takahashi et al. (2014)	NA
Software and Algorithms
PLINK Version 1.90	Chang et al. (2015)	https://www.cog-genomics.org/plink/
R 3.6.3	R Core Team	https://www.r-project.org/
MEGA X	Kumar et al. (2018)	https://www.megasoftware.net/
QIIME2-2019.9	Bolyen et al. (2019)	https://qiime2.org/
Cutadapt Version 2.4	Martin et al. (2011)	https://cutadapt.readthedocs.io/en/stable/
PICRUSt2	Douglas et al. (2020)	https://huttenhower.sph.harvard.edu/picrust/
SParse InversE Covariance Estimation for Ecological ASsociation Inference	Kurtz et al. (2015)	https://github.com/zdk123/SpiecEasi
Mass Hunter Professional software	Agilent Technologies Inc.	https://www.agilent.com/en/product/software-informatics/mass-spectrometry-software/data-analysis/mass-profiler-professional-software
GEMMA	Zhou and Stephens (2012)	https://bioinformaticshome.com/tools/gwas/descriptions/GEMMA.html
Cytoscape	Shannon et at., 2003	https://cytoscape.org/
SILVA 132	Quast et al. (2012)	https://www.arb-silva.de/
Metlin	Guijas et al. (2018)	https://metlin.scripps.edu/landing_page.php?pgcontent=mainPage
HMDB	Wishart et al. (2018)	https://hmdb.ca/
MetaCys	Caspi et al. (2020)	https://metacyc.org/
Other
pH meter FE28	METTLER TOLEDO	https://www.mt.com/cn/zh/home/products/Laboratory_Analytics_Browse/pH-meter/benchtop-pH-meter/fiveeasy/fiveeasy-plus-fe28-meter-cn.html
CHN Spec Colorimeter CS-10	CHNSpec	https://www.lonroy.com/product/cs-10-digital-handheld-colorimeter
Texture analyzer	TMS-PRO	https://www.foodtechcorp.com/products
Agilent 1290 Infinity II UPLC system	Agilent Technologies	https://www.agilent.com.cn/zh-cn/product/liquid-chromatography/hplc-systems/analytical-hplc-systems/1290-infinity-ii-lc-system