Covariation of the gut microbiome with diet in Aves CURRENT STATUS: POSTED

Research over the past few decades has revealed a vital role for the gut microbiome in the health of various animals including birds. Multiple factors can influence the gut microbiome. Opportunistic feeding and multiple other environment factors can influence the results, and bias the conclusions, when wild animals are used to study the influence of phylogeny and diet on their gut microbiomes. Therefore, to study this question in this study, we collected fecal samples from 43 species of Aves at one time to avoid influences such as geography, weather, and season. Results: Approaches based on both 16S rRNA gene sequencing (135 samples) and whole metagenome shotgun sequencing (17 samples) were used. Our data show that diets containing native starch will increase the abundance of Lactobacillus in gut microbiome, while those containing plant-derived fiber will mainly enrich the levels of Clostridium . Greater numbers of Fusobacteria and Proteobacteria are detected in carnivorous birds, while in birds fed a commercial corn-soybean basal diet, a stronger inner-connected microbial community containing Clostridia and Bacteroidia was enriched. Furthermore, a microbial functional analysis based on the metagenomic sequences showed that the function of microbes was adapted to different food types to achieve the most beneficial state for the hosts. suggests that the fitness of the structure and function of microbiota assembly to food were beneficial to the growth and health of the host. Adaptive evolution of microbes to food types formed a dietary-microbiome-host interaction reciprocal state. In general, this study, by determining the modulation of

Approaches based on both 16S rRNA gene sequencing (135 samples) and whole metagenome shotgun sequencing (17 samples) were used. Our data show that diets containing native starch will increase the abundance of Lactobacillus in gut microbiome, while those containing plant-derived fiber will mainly enrich the levels of Clostridium. Greater numbers of Fusobacteria and Proteobacteria are detected in carnivorous birds, while in birds fed a commercial corn-soybean basal diet, a stronger inner-connected microbial community containing Clostridia and Bacteroidia was enriched.
Furthermore, a microbial functional analysis based on the metagenomic sequences showed that the function of microbes was adapted to different food types to achieve the most beneficial state for the hosts.

Conclusions:
The covariation of diet and gut microbiome identified in our study demonstrates modulation of the gut microbiome by dietary diversity and expands our knowledge of diet-microbiome-host interactions in birds.

Background
Digestive tracts of animals contain microbial communities that are composed of different bacterial groups with varying abundances and functional characteristics [1]. Our understanding of the role of the animal intestinal microbiome has now extended far beyond its importance for digestion and energy acquisition, with many recent studies showing that the microbiome contributes to detoxification, immune system development, and behavior [2,3]. Obviously, this means that the adaptive capacity and health status of an animal is not solely determined by the its genome, but must also include the vast genetic repertoire of the microbiome [4].
Birds exhibit the most diverse range of ecological functions among vertebrates [5] and represent a highly evolved lineage that plays important roles in natural ecosystems, and also in agriculture by disseminating crop seeds and capturing harmful insects [6]. Coevolution between host and microbial lineages has played key roles in the adaptation of mammals to their diverse lifestyles, but has been far less studied in birds [7]. Diet is likely an important correlate for microbial diversity during the longterm evolution of Aves. The diets of birds range from robust and opportunistic to strictly carrion feeding or nectar feeding [8]. Exploitation of a new dietary niche is a powerful driver for changes in gut microorganisms and their co-evolution with their animal host [9]. To date, however, studies of gut microbiota in birds have focused on ornamental and economically important birds [10,11], with fewer studies characterizing the gut microbiome in birds associated with highly diverse dietary habits.
The dynamic gut micro-ecological system that formed in the adaption of birds to their environment is influenced by many factors, such as sex, reproductive status, age, geography, and social structure [11][12][13][14]. The dominant drivers of gut microbiome diversity in birds appear to be host evolutionary history and diet [15]. Diet drives not only taxonomic diversity but also the functional content of the gut microbiome in mammals [16]. In Aves, however, most studies to date have focused on just the microbial taxa and has ignored the significance of microbial metabolic functional adaptation to dietary diversity. Moreover, it is difficult to address questions on the identity of host-specific microbes, as the microbial species were often collected in different habitat niches and environmental conditions (e.g., weather and season) that have roles in influencing the composition of the gut microbiota [17]. Therefore, to better reveal the complex relationship between diet and the gut microbiome in birds, it is necessary to minimize the influence of any interfering factors.
In this study, fecal material was collected from birds housed at Guangzhou Zoo representing 37 species with various dietary habits at one point in time to minimize the influence of external factors such as geography, weather, and season. In addition, we also compared the gut microbiota of 6 domestic poultry species that were fed with different types of food to identify differences in gut microbiota that can occur in a species as a response to different food types. To understand the dietary and phylogenetic effects on the taxonomic composition and metabolic function of gut microbiota in birds, a systematic analysis combining 16S rRNA amplification and metagenomic sequencing was used.

Gut microbial diversity of birds
To assess microbial diversity, we sequenced the V3-V4 regions of the 16S rRNA gene to identify 2,459 OTUs (operational taxonomic units) in 135 fecal samples from 43 species (37 wild bird species and 6 domestic poultry, Supplementary Fig. S1 and Supplementary Table 1). We first assessed the impact of general feeding habits on gut microbiome diversity, noting however, that the food types of omnivores vary widely (Supplementary Fig. S2 and Fig. 1a). When we grouped according to six classes of food types (fruits, corn-soy, grains, foliage, flesh, and omnivore) we found that 215 OTUs were shared by all groups. The food type with the highest number of unique OTUs was the fruit food group (271 OTUs), followed by the omnivore group (221 OTUs) (Fig. 1b). In contrast, the fewest unique OTUs were detected in the corn-soy (33 OTUs), grain (22 OTUs) and foliage (13 OTUs) food groups. Moreover, most OTUs (90%) were only detected in less than 20% of samples ( Supplementary Fig. S3).
Next, we used Chao1, phylogenetic diversity (PD whole tree index) and Shannon index to illustrate bacterial richness and diversity within the communities based on the OTUs level. The alpha diversity index of the grain food group was the lowest, and significantly different from the fruit, corn-soy, flesh and omnivore food groups. For example, the grain group was significantly lower than the fruit group with the Chao1, PD and Shannon index. (FDR p < 0.05). No significant differences were observed among the other groups ( Fig. 1c and Supplementary Table S2). A PCoA (principal coordinates analysis) analysis based on the Bray-Curtis distances was used to assess the differences in bacterial community structure between the samples. The results of this analysis revealed a significant clustering of gut microbiota by all diet groups, except the omnivore food group (p < 0.001), with food types separating the microbial communities along the first principal coordinate (PC1, 17.6% of variance) (Fig. 1d).
To further assess the effects of host phylogeny and diet on gut microbiome diversity, we compared the host phylogenetic tree with a UPGMA tree of the gut microbiota ( Supplementary Fig. S4). Gut microbiomes of the different species mainly clustered based on food types. Furthermore, the gut microbiota of domestic poultry species fed with different food types (e.g., Gallus gallus, Meleagris gallopavo, Anas platyrhynchos and Cairina moschata) was diverse, and mainly clustered according to their food types ( Supplementary Fig. S4). Based on the mantel test, both host diet (r = 0.1938, pvalue = 0.0001) and phylogeny (r = 0.137, p-value = 0.0072) affect the gut microbiota of birds.
MaAsLin2 analysis (adjusted p-value < 0.05) further revealed that diet plays a major role, and accounts for 92% of the microbiota features ( Supplementary Fig. S5).

Predominant Bacteria Are Influenced By Dietary Differences
OTUs and genera were sparsely distributed in all samples ( Supplementary Fig. S3). However, the predominant bacterial phyla present in the feces of all birds were Firmicutes (mean abundance ranged from 32.97-72.46%), Proteobacteria (12.21%~37.78%), Acitinobacteria (1.62% ~ 9.8%) and Bacteroidetes (0.07% ~ 14%) ( Fig. 2a and Supplementary Table S3). At the genus level, the five most abundant genera were Lactobacillus, Clostridium, Enterococcus, Escherichia, and Turicibacter ( Fig. 2b and Supplementary Table S3). Consistent with the α-diversity characteristic, the number of genera with mean relative abundance > 1% in the omnivore food group was higher than in the other groups.
Differences in abundance of the bacterial taxa was determined through a LEfSe analysis and a total of 28 taxa at different classification levels were found to have significant differnces (p < 0.05, LDA > 4) ( Fig. 2c and Supplementary Fig S6). At the genus level, Lactobacillus, Leuconostoc, Clostridium and Cetobacterium were the dominant genera, and were significantly abundant in the grain, fruit, foliage and flesh food groups, respectively (p < 0.05). At the order level, a significant enrichment of Bacteroidales, was detected in the corn-soy food group. No significantly enriched taxa were observed in the omnivore food group.

Microbial Co-occurrence Association Patterns Are Influenced By Dietary Differences
We next examined how bacterial species co-occur among the birds, which might be due to dietary differences or microbe-microbe interactions. A network contained 344 nodes and 2559 edges was constructed ( Fig. 3a), with the set of topological metrics including average degree, average weighted degree, average path degree, density and average clustering coefficient listed in Supplementary   Table S4. Based on the layout structure, this integrated network could be divided into 6 sub-networks, each with differing taxonomic compositions at the class level (Fig. 3a). The taxonomic information represented by each node is listed in Supplementary Table S5. At the class level, the co-occurrence network was mainly composed by interactions of Clostridia, Bacteroidia, Gammaproteobacteria and Bacilli. Furthermore, this network consisted of six subnetworks, each with differing taxonomic compositions. Subcommunity a (SC-a), subcommunity b (SCb) and subcommunity c (SC-c) were dominated by Clostridia and Bacteroidia; subcommunity d (SC-d) and subcommunity e (SC-e) possessed more members of Actinobacteria, Gammaproteobacteria and Alphaproteobacteria; SC-f was mostly composed of taxa from Bacilli. In addition, most of the interactions were positive, while negative interactions only appeared between Pseudomonas (classified in Gammaproteobacteria) in SC-d with Ruminococcaceae UCG-014, and between Intestinimonas and Christensenellaceae R-7 group (classified in Clostridia) in SC-a (Fig. 3a), suggesting that competitive inhibition between these pairs of communities. Interestingly, we found Parasutterella was the only gram-negative genera that co-occurred with other gram-positive bacteria in SC-a.
We then calculated the co-occurrence percentage and total abundance of each sub-microbial community to estimate the microbial coexistence in the different groups. The presence and abundance of OTUs from each sub-network differed substantially among the groups (Fig. 3b). SC-d was generally the most prevalent in all birds, while wide differences in the prevalence and abundance of SC-a and SC-f occurred among the groups, suggesting that certain host dietary specificity of these microbial consortiums. The abundance of SC-a in the corn-soy group was significantly higher than in the other groups, and the abundance of SC-d in the corn-soy group was significantly lower than in the other groups (p < 0.05) (Supplementary Table S6). This phenomenon is consistent with antagonism between SC-a and SC-b as described above.

Adaptive Evolution Of Microbial Functions To Fit Food Types
To further investigate the functional capacities of the gut microbial communities in birds, a metagenomic analysis was conducted. In total, 2,425,998 assembled genes (92.02%, 2,636,348) were identified from the prokaryotic microbes and fungi by searches against the NCBI NR database. Of this total, 1,733,474 (65.75%) and 52,110 (1.98%) were annotated in the KEGG database and CAZy database, respectively.
Detailed annotation KOs information is listed in Supplementary Figure S7. Notably, 2,508 (28.47%) of the KOs were annotated in global and overview metabolism maps and 741 KOs were annotated in carbohydrate metabolism. The number of KOs with an average relative abundance higher than 0.01% in the different groups ranged from 1,329 to 2,672, and the total abundance of those KOs in each group was higher than 85% ( Supplementary Fig. S8). This indicates that the high abundance KOs cover most of the microbial functions.
To compare the microbial functions between each group, we first tested the KEGG pathway enrichment analysis based on the top abundance KOs (mean relative abundance > 0.01%) in each group. The top 20 significantly enriched KEGG metabolism pathways in each group are shown in Moreover, we explored the distribution of CAZymes in the different groups. The top 20 abundant CAZymes in each group are shown in a z-score normalized heatmap (Fig. 4b). This data showed that most of the high abundance CAZymes were detected in the corn-soy, flesh and grain groups. Groups that have diets containing plant-derived fiber (fruit, omni and foliage) had similar enzyme profiles and clustered together. Discussion community composition in the vertebrate gut [15,18]. When wild animals are used to study the influence of diet and phylogeny in the composition of the gut microbiome, other factors, such as habitat, weather, and season can bias the analysis and conclusions. In this study, we studied multiple species of Aves at one time to minimize these confounding factors.
In our study, different species of birds eating the same food composition tended to have similar gut microbiomes. For example, the gut microbiomes of birds with fish and meat as their diet, such as Pelecanus onocrotalus, Ciconia boyciana, and Aquila chrysaetos, clustered together based on the UPGMA analysis ( Supplementary Fig. S4). Although not all samples clustered according to food characteristics, our results support the conclusion that food source is a major factor determining differences in intestinal microbial composition [18]. The MaAsLin2 analysis further supported the conclusion that diet has a greater effect than phylogeny on the gut microbiome ( Supplementary Fig.   S5).
The predominant bacteria found in the intestinal tracts of birds fed different types of food vary greatly. The relative abundance of the Lactobacillus, classified in Bacilli, in the grain fed group reached more than 60% abundance and was significantly higher than for any other groups (p < 0.05).
Indeed, the ranking of food types, from high to low in Lactobacillus abundance was grain, corn-soy, fruit, foliage and then flesh ( Fig. 2 and Supplementary Table S3). Starch is the major storage polysaccharide in cereal grains, legumes, and in many roots and tubers [19], and the above result indicates that the starch content of the food was positively correlated with Lactobacilli abundance. It has been reported that starch is the only polysaccharide hydrolysised by the extracellular enzymes (amylopullulanase) of Lactobacilli [20]. In addition, analysis of the metabolic pathways from the metagenome analysis showed that birds with greater plant-derived polysaccharide intake (grain, corn-soy, fruit, foliage, and omnivore food groups) have microbiota that were more highly enriched in genes that function in starch and sucrose metabolism compared to carnivores (group flesh) (Fig. 4a).
Although birds can secrete pancreatic amylase, they have limited ability to digest native starch as it is highly organized [21]. Taken together, we conclude that a high abundance of Lactobacillus in the gut microbiota of birds is essential for the metabolism of native starch.
Plant-derived fiber is the main energy source for folivores, but dietary fiber utilization by birds is inefficient and variable due to the absence of enzymes that can digest fiber [22,23]. We found that Clostridium was significantly enriched in the foliage food group (p < 0.05), followed by the fruit and omnivore food groups (Fig. 2c, Supplementary Fig S6). A similar CAZymes profile appeared in these 3 groups (foliage, fruit and omnivore) whose diets contain high levels of plant cellulose (Fig. 4b).
Enzymes in Clostridium digest fibers and produce various metabolites such as short-chain fatty acid (SCFA) that can be used by the host [24,25], therefore, the enriched levels of Clostridium in the guts of birds that intake plant-derived fiber could make up for the lack of fiber digestive enzymes in the host.
However, some of predominant microbes might simply be due to food intake. For example, the relative abundance of Cetobacterium, which was identified as a gut microbe of various freshwater fish [26], was significantly higher in flesh eating group than in other groups. Since that food eaten by carnivorous birds contains freshwater fish, the enriched levels of Cetobacterium in these birds might directly be from fish intake.
Consistent with previous studies, the phyla Firmicutes and Proteobacteria were predominant among all food groups [27][28][29]. However, significant difference in the microbiota diversity and richness of these phyla was observed between the groups ( Fig. 1c and Fig. 1d). In the corn-soy group, the prevalence of Clostridia and Bacteroidia resulted in a relatively high α-diversity ( Fig. 1 and Fig. 3).
These strong inner-connected microbes could display a complex interaction web and reveal specific microbiome diversity and function within certain environments [30,31]. Various genera, such as Bacteroides, Butyricicoccus, Faecalibacterium and Ruminiclostridium, in the classes Clostridia and Bacteroidia are beneficial symbiots and act as SCFAs producer [32]. Similarly, microbes with functions enriched for lipid metabolism, including glycerolipid metabolism and fatty acid biosynthesis, are higher in the corn-soy group than in other groups (Fig. 4a). Previous work had demonstrated that diets high in soybean proteins significantly elevate the production of SCFAs by the gut [33], thus, high soybean protein content might have led to an elevation of SCFAs in the corn-soy group.
In addition, the amino acid composition of natural foods is lower than in commercial feeds, which have a well-balanced amino acid profile when soybean meal and other crude protein are used as the basic ingredients [34]. Moreover, essential amino acids (EAAs), which cannot be synthesized by vertebrates but can be synthesized by gut microbes [35], are usually added as additives to commercial feeds. Our data showed that gut microbiota have a strong ability to synthesize amino acids in birds that eat natural foods, with metagenomes significantly enriched for genes involved in the biosynthesis of branched-chain amino acids including valine, leucine, and isoleucine (Fig. 4a). The above result suggests that when the direct source of amino acids is limited in natural foods, then adapting the microbiome could allow enhanced synthesis of these amino acids.
In addition to digestion, gut microbes also play important role in the health of the host. Our data showed that Proteobacteria was prevalent in all birds (Fig. 3), a result similar to previous studies in Passerines [36]. Most bacteria in this phylum, such as Pseudomonas, Ralstonia and Acinetobacter, are gram-negative and act as pathogens or opportunistic pathogens in their hosts [37]. The levels of these gram-negative bacteria (in SC-d) were negatively related to the levels of Clostridium and Bacteroidia (in SC-a) (Fig. 3a). Moreover, pathogens can be inhibited by SCFAs produced by Clostridium and Bacteroides [38]. These characteristics might be the reason why the richness of these pathogens was significantly lower in the corn-soy group compared to other groups ( Fig. 3b and Supplementary Table S6), further substantiating the conclusion that dietary composition influences the interactions between microbes to modulate intestinal microbiome to reduce enteric pathogens [39].

Conclusions
In conclusion, our study identified food source as the main factor modulating gut microbiome diversity, rather than bird phylogeny, after minimizing the other complex interfering factors such as weather, season and geography. The covariation of diet and gut microbiome identified in our study suggests that the fitness of the structure and function of microbiota assembly to food were beneficial to the growth and health of the host. Adaptive evolution of microbes to food types formed a dietarymicrobiome-host interaction reciprocal state. In general, this study, by determining the modulation of the bird gut microbiome by dietary change, expands our knowledge of dietary-microbiome-host interactions in birds.

Materials And Methods Sample collection
Sampling was conducted in May 2019. In total, 108 fresh fecal samples were collected from 37 species of birds housed at Guangzhou Zoo, Guangzhou city, China. No antibiotic drug use was recorded for any of these birds within the 6 months prior to sampling. In addition, 27 fresh fecal samples from 6 species of domestic poultry (such as chicken, duck and goose) were collected at the same time from a farm in Guangzhou. The host phylogeny and sample metadata associated with this study are shown in Supplementary Figure S1 and Supplementary Table S1, respectively. Each fecal sample was collected in a 5 mL sterile sampling tubes and transported with dry ice. Samples were stored at -80 °C until DNA extraction.
The bird samples were classified into 6 groups (omnivorous, frugivorous, granivorous, folivore, carnivorous and piscivorous) according to their general feeding habits. In addition, based on their food composition, the species were also divided into 6 food types: fruits (mainly include apple, banana), corn-soy (commercial corn-soybean basal diet supplement with feed additive), grains (unprocessed rice and wheat), foliage (plant stem and leaves), flesh (mainly include fish, mouse, rabbit) and omnivore (multiple types of food). Detailed information on these groups is listed in Supplementary Table S1.

DNA Extraction And Sequencing
For 16S rRNA gene sequencing, PCR amplicons for the V3-V4 region were generated with primers 338F-806R [40]. Sequencing libraries were generated using TruSeq Nano kit (Illumina) and assessed on the Qubit @ 2.0 Fluorometer (Thermo Fisher Scientific) following manufacturer's recommendations.
Sequencing was performed on an Illumina MiSeq using 2 × 250 bp paired-end V2 MiSeq reagent kits (Illumina) and an average of 50,994 tags per sample were produced.
Metagenomic sequencing was performed by Novogene Bioinformatics Technology Co., Ltd, China.
Briefly, bacterial cells were separated from undigested food particles and recovered through differential centrifugation before cell lysis. DNA was isolated with the Qiagen QIAamp DNA Stool Mini libraries were prepared with an insert size of 350 bp generated with NEBNext® Ultra™ DNA Library Prep Kit for Illumina® (New England Biolabs). Sequencing was performed on an Illumina Novaseq 6000 platform, with an average of 6.4 G raw data per sample produced.

16S rRNA Data Handling
Paired-end reads were assigned to samples based on their unique barcodes and truncated by cutting off the barcode and primer sequences using the Cutadapt (version 1.18) pipeline [41]. Usearch (v11) [42] was used for remove redundant sequences and then compared with the reference Gold database and SILVA databases (v132) [43] to remove chimeric and non-bacterial sequences. Sequences were clustered into operational taxonomic units (OTUs) at a similarity threshold of 97% using the UPARSE algorithm [44]. OTUs were subsequently mapped to the SILVA database to determine taxonomies, and samples with ≥ 1,000 mapped reads were used for the downstream analyses. General manipulation and basic analyses of the dataset were performed in QIIME and R with the phyloseq, vegan, and ggplot2 packages [45,46].
For with-in-sample microbial diversity analysis, Shannon, PD whole tree and Chao 1 indexes were used. A non-parametric t test was used to calculate significance of the α-diversity. Overall differences in the bacterial community structures were evaluated using PCoA based on the Bray-Curtis dissimilarity values, and the significance in the difference of the community compositions between groups was determined by permutational multivariate analysis of variances (PERMANOVA). UPGMA (Unweighted Pair Group Method with Arithmetic Mean) was performed with upgma_cluster.py in QIIME. The host phylogeny was obtained from http://timetree.org/. LefSE was used to show the comparison and identify significantly different bacterial species between each group by performing a linear discriminant analysis (LDA) effect size (p < 0.01, LDA > 4) [47].
The Mantel test was performed with the Vegan R package. Effect size and significance is derived by comparing the true data to randomly generated permutations (n = 9999 for all analyses). The host phylogeny was represented by the patristic distance (branch lengths). The OTUs table was converted to a Bray-Curtis distance matrix among all pairwise sample comparisons. Euclidean distances for the diet and phylogeny data were calculated. Multivariable associations between diet, phylogeny and OTUs were determined using MaAsLin2 (http://huttenhower.sph.harvard.edu/maaslin2) with default parameters. All tests with a Benjamini-Hochberg adjusted p value of < 0.05 were considered to be significant.

Co-occurrence Network Analysis
To reduce rare OTUs in the dataset, OTUs with a relative abundance < 0.01% were removed, leaving 506 OTUs for the following analysis. A co-occurrence network analysis was conducted as described in a previous study [48]. In detail, spearman's rank coefficients (ρ) between the OTUs were calculated pairwise using the R package co-occur [49]. Subsequently, significant and robust correlations (FDR pvalue < 0.01, |ρ| ≥ 0.6) were used to construct a network using the R package igraph [50]. Gephi with the layout algorithm of Fruchterman-Reingold was used for network visualization [51]. The relationships among the microbial taxa were estimated by establishing a correlation network, which considered both positive (Spearman's ρ > 0.6) and negative (Spearman's ρ < −0.6) edges. To reduce the network complexity, unclosed loops and closed loops having less than 4 nodes were not be presented.
In total, 6 subcommunities (SC) were observed in the co-occurrence network. To calculate the cooccurrence percentage of each microbial community in each sample, the presence of OTUs number of each subcommunity was divided by the total number of OTUs for the subcommunity. To calculate total abundance of a subcommunity in each sample, the relative abundances of present OTUs were added together.

Shotgun Metagenomics Data Analysis
De novo metagenome assembly was according to previous research [52]. Raw reads were cleaned to exclude adapter and low-quality sequences with fastp (v 0.19.7) [53]. Contamination reads were be discard after mapping the high-quality reads to a nucleotide dataset containing the chicken, maize, rice and zebrafish genomes with BWA-MEM (v 0.7.17) [54]. For each sample, clean reads were assembled by Megahit (v1.0.3) under the pair-end mode [55]. Contigs larger than 500 bp were used for ORF prediction with Prodigal (v2.6.3) [56], and the predicted CDS with lengths less than 102 bp were filtered out. The initial non-redundant gene set was produced by CD-HIT [57]. Taxonomic assignment of protein sequences was made on the basis of a DIAMOND (v0.8.28.90) [58] alignment against the NCBI-NR database and genes classified as eukaryote were excluded.
For functional profile generation, the protein sequences of all remaining genes were aligned to the Kyoto Encyclopedia of Genes and Genomes (KEGG) database and the carbohydrate-active enzyme (CAZy) database. Briefly, KEGG annotation was made on the basis of the DIAMOND alignment by taking the best hit with the criteria of an E value < 1e-5. The carbohydrate-active enzyme (CAZymes) annotation was made on the basis of an HMMER search against the dbCAN HMM (hidden Markov model) database [59].
Sequenced-based gene abundance profiles were calculated as previously described [60]. Relative gene abundance profiles were summarized into KEGG and CAZy functional profiles for further analysis. KEGG pathway enrichment analysis was done with the R package clusterProfiler [61]. All annotated microbial KEGG orthology (KOs) were set as background KOs. Enriched pathways in each group were calculated based on the KOs with mean relative abundance > 0.01%. Mean abundance differences between each group in functional terms were visualized as a heatmap and clustered by hierarchical clustering using the pheatmap package in R.

Statistical analysis
A t test was employed to determine the significant differences in the α-diversity between all groups.
The PERMANOVA test was carried out to test whether the β-diversity was significantly different between all groups. A Wilcoxon rank-sum test was used to statistically compare the co-occurrence percentage and abundance difference of sub-microbiota communities. P-values were adjusted for multiple testing using the Benjamini-Hochberg (BH) false discovery rate (FDR) controlling procedure.

Ethics approval and consent to participate
Feces sample collection was approved by from the Forestry Bureau of Guangdong province, China.          showing the taxonomic levels are represented by rings with phyla in the outermost ring and genera in the innermost ring. Circles with non-yellow color indicate that there is a significant differences in the relative abundance at the different taxa levels (P < 0.05; LDA score > 4) and yellow circles indicate non-significant differences.