All applicable international, national, and/or institutional guidelines for the care and use of animals were strictly followed. All animal sample collection protocols complied with the current laws of China. According to the Law of the People’s Republic of China on the Protection of Wildlife (Order of the President of the People’s Republic of China, No. 16, 2018), all experimental procedures carried out in this study were under permitted. The non-invasive genetic sampling taken in this work enables collecting samples (fecal samples) without the need to directly disturb or contact individuals of oriental white storks.

Twelve fecal samples of oriental white storks (Ciconia boyciana) were collected at Tianjin Zoo (Z Group) and Tianjin Qilihai Wetland (W group), respectively in 2019 (Table 1). Tianjin Qilihai Wetland is one of the main stopover habitats of wild C. boyciana during their migration periods. All samples were carefully taken from the inside of the fecal matter using sterilized equipment, and then were placed on dry ice and transported to the laboratory as soon as possible.

The total genomic DNA of all fecal samples was isolated by CTAB extraction method; Sodium dodecyl sulfate (SDS) extraction method; Guanidine thiocyanate (GuSCN) extraction method and TIANamp Stool DNA Kit (TIANGEN, China). For all methods, DNA extraction was conducted according to the references and manufacturer’s instructions (Hosomi et al., 2017; Kachiprath et al., 2018). The concentration and quality of DNA samples were determined by a NanoDrop 2000 spectrophotometer (Thermo Fisher Scientific, Waltham, MA, United States).

The V3–V4 region of the 16S rRNA gene was amplified based on the literature using specific PCR primers (Wu et al., 2016). Forward primer 338F (5′-ACTCCTACGGGAGGCAGCAG-3′) and reverse primer 806R (5′- GGACTACHVGGGTWTCTAAT-3′) were used (Caporaso et al., 2012). Amplification was carried out in 50-μL reactions with 10 ng template DNA and 25 μL Pfu DNA polymerase mix. The PCR cycle comprised an initial denaturation at 94°C for 5 min, followed by 30 cycles of 94°C for 15 s, 60°C for 15 s, and 72°C for 30 s.

Purified PCR amplicons were sequenced on an Illumina MiSeq platform at Majorbio Bio-Pharm Technology Co., Ltd., Shanghai, China. PE amplicon libraries were constructed, and all raw reads were screened according to the primer sequences.

Raw files of the total genomic DNA from fecal samples were demultiplexed by QIIME (version 1.9.1), and high-quality reads were clustered as an operational taxonomic unit (OTU) by UPARSE pipeline (version 7.1) when the sets of sequences shared at least 97% identity (Edgar, 2013). In this study, gene sequence taxonomy was analyzed by RDP Classifier (version 2.2)¹ (Wang et al., 2007) using a confidence threshold of 80% against the Silva 16S rRNA database (Release 128)² (Quast et al., 2013).

Rarefaction curves were plotted for each sample to determine the abundance of communities and sequencing data of each fecal sample (Amato et al., 2013). α-diversity and species accumulation curve were calculated by mothur software (version v.1.30.1), which showed the complexity of gut flora species in two groups. β-diversities analysis was performed and principal coordinate analyses (PCoA) based on OTU compositions were determined, in order to investigate structural variation in microbial communities of two groups using unweighted and weighted UniFrac distance metrics (Lozupone and Knight, 2005). The one-way analysis of variance ANOVA (Permutational multivariate analysis of variance) was performed to compare α-diversity estimates and p < 0.05 was considered statistically significant. Differences in the UniFrac distances for pairwise comparisons were determined using the Student’s t-test.

Taxa abundances in two groups at the phylum, class, order, family, and genus levels were statistically compared by Metastats (an improved statistical method for analysis of metagenomic data) (White et al., 2009). Microbial functions were predicted by PICRUSt based on high-quality sequences and the linear discriminant analysis effect size (LEfSe) was used to present bacterial taxonomic distributions of sample communities (Segata et al., 2011).

DNA extraction is a very critical step for high-throughput sequencing and the choice of a DNA extraction method has an important impact on DNA yield and purity. As the results shown, the GuSCN extraction method and the TIANamp Stool DNA Kit produced much higher DNA yields than the CTAB and SDS methods (Table 2). The obtained DNA yields by all the four extraction methods were higher than 10 ng/μl. The 260/230 absorbance ratios GuSCN extraction method and TIANamp Stool DNA Kit were above 1.8, indicating a low contamination with proteins. However, the average 260/280 absorbance ratios detected for the CTAB and SDS were below 1.8 which might suggest protein contamination from the two methods. In the study of oriental white stork’s fecal sample, the DNA obtained by GuSCN extraction method and TIANamp Stool DNA Kit could both be used successfully in library preparation for high-throughput sequencing. Moreover, DNA extraction using the GuSCN protocol was more cost efficient compared to the commercially TIANamp Stool DNA Kit.

Across the 24 samples of oriental white storks, a total of 1,263,339 raw reads with an average length of 415 bp were obtained after denoising by Illumina MiSeq sequencing. The total number of raw reads, base pairs, and the mean length of the reads were compared in all the 24 samples (Supplementary Table 1). The number of sequences varied from 42,700 to 70,383, and a total of 2,390 OTUs with 97% sequence similarity threshold representing 43 bacterial phyla were identified in the feces of C. boyciana. The α-diversity indices for bacteria, including Sobs, Shannon, Chao1, ACE, and Good’s coverage are shown in Table 3. The total number of OTUs in the W group (2,134) was much higher than that in the Z group (1,555), these values include 1,299 OTUs shared by both groups. There was a significant difference in α-diversity indices between Z group and W group (p > 0.05) (Supplementary Figure 1). Average Good’s coverage was 99.28 ± 1.57% (mean ± SD) around all 24 fecal samples, and thus the majority of the bacterial diversity in feces was obtained thereby ensuring accuracy of the analyses (Table 3).

To assess whether the depth of sequencing in this study was large enough to estimate the richness of gut microflora, the observed species rarefaction curves were calculated. According to species accumulation curves and Good’s coverage estimator, the curve tended to gradually flatten, indicating that the sample size collected was large enough to identify the majority of OTUs in the gut of oriental white storks (Supplementary Figure 2), and the high Good’s coverage indicated that the majority of bacterial phylotypes in all the samples were identified.

The taxonomic composition of each fecal sample clustering at the 97% phylotype similarity level was outlined at phylum, class, order, family, and genus levels (Supplementary Table 2). Overall gut microbiomes of the oriental white storks comprised of 43 phyla, 96 classes, 188 orders, 338 families, 675 genera, and 1,222 species; microbial compositions at these different levels are shown in Figure 1 and Supplementary Figures 3–5, respectively.

At the phylum level, 43 prokaryotic phyla were identified from the 16S rRNA gene sequences (Figure 1A). Gut microbiota of C. boyciana from Z group was dominated by Firmicutes (65.59%), Actinobacteria (17.23%), and Proteobacteria (10.48%). These three dominant phyla accounted for 93.30% of all sequences across all the samples, while unclassified bacterial sequences only occupied 0.01% (Figure 1A and Supplementary Table 2). The gut microbiota from W group was dominated by Firmicutes (28.43%), Proteobacteria (23.58%), Actinobacteria (14.71%), Chloroflexi (11.57%), and unclassified bacterial sequences accounted for 0.04% of all sequences in the W group (Figure 1A and Supplementary Table 2).

At the genus level, 15 genera had abundances greater than 2% in all the sequences (Figure 1B). The most abundant genera in W group were Lactobacillus (8.33%), norank_f_Anaerolineaceae (5.92%), norank_f_Xanthomonadales_Incertae_Sedis (3.91%), norank_f_Bacteroidales_S24-7_group (3.38%), norank_c_ Cyanobacteria (3.35%), and Thiobacillus (2.943%). The most abundant genera in Z group were Paeniclostridium (23.82%), Lactobacillus (8.55%), Peptostreptococcus (7.78%), Clostridium_sensu_stricto_1 (7.17%), Actinomyces (3.99%), and Romboutsia (3.77%). Mean relative abundances of the 10 most abundant genera are listed in Table 4. Community heatmap analysis at the genus level of the 24 samples can be seen in Figure 2A.

α-diversity, which refers to the diversity within a particular region, was calculated for two groups to examine differences between oriental white storks from Tianjin Zoo and Tianjin Qilihai Wetland. Community richness measured by OTUs (Sobs and Chao indices) showed that the index of the W group was higher than that of the Z group, indicating that richness of gut microbiota was greater in W group (Supplementary Figure 1). Differences between two groups were statistically significant (p < 0.05) (Figure 2B). Bacterial community diversity (Shannon index) was measured in the two groups, and the index of W group was higher than that of Z group (Supplementary Figure 1). There was a significant difference in community diversity between the W and Z groups (p < 0.05) (Figure 2C).

The relationship between the feces samples from two breeding conditions were represented by a dendrogram using Bray-Curtis distances (Figure 3A). This dendrogram indicated that the gut microbiota of the oriental white stork in Tianjin Zoo clustered together, while those in Tianjin Qilihai Wetland located on similar branches. To evaluate the difference in β-diversity of the 24 samples, PCoA was used to visualize similarities or dissimilarities between the Z and W groups (Figures 3B,C). Each symbol represents one gut microbiota on the PCoA plot, and most of the fecal samples in two groups were distinguishable by PCoA when considering the relative abundances of gut microbiomes. Congruent with the cluster analysis, bacterial communities of group Z clustered tightly and were separated from those of group W under both weighted and unweighted UniFrac metrics. The two groups tended to cluster separately from each other, indicating that each group had distinctive microbial populations. Using weighted UniFrac distance, bacterial communities of group Z were separated from those of group W along principal coordinate axis 1 (PC1) with the largest amount of variation (42.10%) (Figure 3B). Using unweighted UniFrac distance, bacterial communities of the 24 fecal samples explained the largest amount of variation (47.89%) along PC1 (Figure 3C). This clustering pattern was confirmed by an analysis of similarities (ANOSIM) and unweighted UniFrac showed a clearer separation of communities (ANOSIM, R = 0.71). Microbial composition between two groups at the phylum level (ANOSIM, R = 0.5621, p = 0.001) and genus level (ANOSIM, R = 0.6196, p = 0.001) were significantly different (Supplementary Figure 6).

Unique and shared bacterial taxa of the gut microbiota from two groups were analyzed. A Venn diagram demonstrated that 1,299 OTUs from the 24 samples were shared as core bacterial OTUs and the unique OTUs in each group were 256 (Z group) and 839 (W group), respectively (Supplementary Figure 7).

There were some similarities in gut microbiota of oriental white storks from two groups. The top five abundant core phyla were Firmicutes, Proteobacteria, Actinobacteria, Chloroflexi, and Bacteroidetes, while the top five abundant core genera were Paeniclostridium, Lactobacillus, Peptostreptococcus, Clostridium_sensu_stricto_1, and Actinomyces. Next, LEfSe was used to identify OTUs differentially represented between oriental white storks from two groups (Figure 4A). LEfSe identified 10 and 28 taxa (LDA > 4.0) with discrepancies in relative abundance in the Z group and W group, respectively (Figure 4A).

The cladogram revealed that the core bacterial species in all 24 fecal samples were significantly different at all levels (Supplementary Figure 8). At the phylum level, the abundance of Firmicutes was significantly higher in the Z group than tin the W group, while the abundance of Bacteroidetes was significantly higher in the W group than in the Z group (p < 0.001) (Figure 4B). At the genus level, Paeniclostridium had a significantly higher abundance in the Z group than in the W group (p < 0.01) (Figure 4C).

A total of 296 KOs (KEGG Orthologs) and 42 KEGG pathways were mapped based on the high throughput sequencing, which were then classified into secondary KEGG pathways. The secondary KEGG pathways associated with the gut microbiotas included cell processes, metabolism, and genetic information processing, indicating that the differences in gut microbiota had important influences on the metabolism of C. boyciana (Figure 5).

Characterization of bacterial communities in the feces of oriental white storks provided evidence that the microbiota compositions were similar between the captivity and wild breeding groups, but abundances of the bacteria were significantly different between two groups. Gut microbiota of C. boyciana were predominantly composed of Firmicutes, Proteobacteria, and Actinobacteria at the phylum level, and Firmicutes was the most abundant phylum of Gram-positive bacteria in both Z group and W group. This microbiota composition was consistent with that of other wild bird species such as Canada goose (Lu et al., 2009), red-crowned cranes (Xie et al., 2016), Darwin’s finches (Michel et al., 2018), chinstrap penguins (Barbosa et al., 2016), northern bald ibis (Spergser et al., 2018), and kakapo (Waite et al., 2014), where members of the phylum Firmicutes dominated the microbiota of avian guts.

Members of the phylum Firmicutes play important roles in metabolism, digestion, and absorption of protein and other nutrients (Bernini et al., 2016; Berry, 2016). For example, Firmicutes are associated with breakdown of fatty acids, and Firmicutes produce more butyrate are which is considered a health-promoting molecule, since it can increase insulin sensitivity, regulate energy metabolism, and increase leptin gene expression. Bacteroidetes mainly produce acetate and propionate. Propionate could reduce the expression of enzymes involved in the de novo synthesis of fatty acids and Acetate could promote the secretion of insulin and ghrelin (Kong et al., 2013; Million et al., 2013; Lee et al., 2015). Although some studies speculated that the Firmicutes-to-Bacteroidetes (F/B) ratio may be related to obesity both in many animal models (Turnbaugh et al., 2006) and humans (Ley et al., 2006), the evidence suggesting an association between obesity and alterations of the Firmicutes/Bacteroidetes ratio is not convincing (Magne et al., 2020). In our study, the F/B ratio in Z group was higher than the ratio in W group, indicating that the gut microbiota of C. boyciana living in Tianjin Zoo could help the host use the nutrients more efficiently and collect calories from food since Firmicutes are more effective. It is difficult to obtain the physical indicators of oriental white storks living in Tianjin Qilihai Wetland, such as weight and health status, thus the relationship between the Firmicutes/Bacteroidetes ratio and obesity of two groups were still unknown.

Paeniclostridium was the most common genus in Z group, while Lactobacillus was the most common genus in W group and was also a major component of the gut microbiota in Z group. The abundance of Lactobacillus was relatively high in both Z group and W group (8.55% and 8.33%, respectively). In previous reports, Lactobacillus were represented in gut microbiota communities that had relatively high β-xylosidase and β-glucosidase levels (Ruggiero et al., 2009; Cabre et al., 2012). The similarity in relative abundance of Lactobacillus in both groups in the current study might be due to the formulation of the diet. The diet of oriental white storks is relatively large, and although the dietary structure is different in two groups, the amounts of sugar consumed in the different environments might be similar. Thus, it would be interesting to study the relationship between dietary structure with habitation and explore the functions of these bacteria in the metabolism pathways of C. boyciana.

The UniFrac method that focused on the phylogenetic relationships was developed by Lozupone and Knight (2005). Both weighted and unweighted UniFrac metrics are very useful for analyzing the differences and associations of microbial communities (Lozupone and Knight, 2005; Campbell et al., 2015). In this study, both weighted and unweighted UniFrac metrics suggested that each fecal sample harbors different microbial communities, and the majority of gut microbiotas were conserved and clustered together among different breeding conditions of oriental white storks (Z group and W group).

PICRUSt was used to analyze the microbial functions of oriental white storks. Gut microbial taxa of C. boyciana were associated with functions such as cell processes and metabolism. These data facilitate understanding of the relationship between gut microbial taxa and metabolism, as well as the influence of gut microbes on host health (Sommer and Backhed, 2013). However, as in other avian studies on gut microbiota, one limitation of the current study is the acquirement of DNA from each individual bird, as this affects the study of gut microbiota in relation to host genetics and physiology at the individual level (Liu et al., 2020).

Analysis of the 24 fecal samples demonstrated that the gut microbiota of oriental white storks was similar to those of red-crowned cranes and black-necked cranes at the phylum level, while Firmicutes was the predominant phylum (Xie et al., 2016; Wang et al., 2020). Comparison of the relative abundance of gut microbiota in C. boyciana between Z group and W group revealed that community structure and abundance of gut microbiota were relatively stable at the phylum level, but were variable and complex at the genus level. There was a significant difference in the richness and diversity of microbial composition, and host intestinal microbiota may be related to diet (Turnbaugh et al., 2009), geography, and seasonal changes (Mueller et al., 2006). The differences between Z group and W group were living conditions and dietary conditions. Oriental white storks inhabiting Tianjin Qilihai Wetland required higher energy and had a more diverse diet, so W group had a greater abundance of gut microbiota than Z group. This study suggests that biogeographical and dietary factors may contribute to the community structure of the gut microbiota of C. boyciana since the metabolism and energy harvest needs were different between two groups. This is consistent with previous related avian findings that community structure and abundance of the gut microbiota are mainly affected by dietary conditions (Wang et al., 2020). The difference on diet patterns may contribute to the differences of fecal microbiota on two groups, and more work is needed to understand the functions of diet patterns and living condition on fecal microbiota. So next, whether living condition changes in the fecal microbiota will be investigated using a combination of metagenomics and metabolomics. Furthermore, since the different microbiota could have different effects on the digestion ability of carbohydrates, protein, and cellulose, the differences in feed composition will have important influences on the changes in host gut microbiota. This may affect nutrient utilization and could cause gastrointestinal diseases in oriental white storks. Hence, studies on the gut microbiota of C. boyciana will facilitate understanding of the specific functions of different bacteria related to digestion and absorption in the host, and could help improve conservation for threatened species both in artificial breeding groups and in wild groups.