Variations in gut bacterial communities of hooded crane (Grus monacha) over spatial-temporal scales

Ethics statement

Fecal samples were collected after foraging to ensure that the hooded cranes were devoid of human disturbance. It does not involve direct hunting or capture. Permission for the collection of fecal samples was obtained from the local government for wildlife management.

Study areas

The study was carried out in Shengjin (30.25–30.50°N, 116.92–117.25°E) and Caizi (30.75–30.97°N, 117.00–117.15°E) Lakes, which were located in the middle and lower Yangtze River floodplain (Fig. 1). Both lakes are globally important wintering and stopover habitats for migratory wading birds on the East Asian-Australasian Flyway (Cao & Fox, 2009; Fox et al., 2011a). Shengjin and Caizi Lakes are river-connected shallow lakes with a northern subtropical monsoon climate. The average annual temperature is 14–18 °C and the average annual rainfall is approximately 1,000–1,400 mm.

The seasonally inundated wetlands provide plenty of food (e.g., Oryza sativa, Polygonum criopolitanu, and Potentilla supina) for the wading birds (Zhao et al., 2013; Zheng et al., 2015). There are approximately 600 individuals of hooded cranes wintering in the two lakes each year. The habitat utilization rate of the cranes has shifting over spatial and time scales (Zhao et al., 2013; Zheng et al., 2015). During the wintering period, shifts in food density and resources triggers cranes to adjust their foraging habitats (Wan, Zhou & Song, 2016; Zhao et al., 2013; Zheng et al., 2015). The wintering period can be divided into three stages according to the migratory rhythm and hydrological processes of the cranes: the early stage from November to December; the middle stage from January to February; and the late stage from March to April (Zhao et al., 2013; Zheng et al., 2015). The foraging habitats of hooded cranes are relatively stable at certain wintering stage in Shengjin and Caizi Lakes (Cao & Fox, 2009).

Sample collection and preservation

Fecal samples were collected during the three-wintering periods. The foraging sites of the crane flocks were observed with a telescope before sampling to ensure that there were no other species. To avoid human disturbance, the fecal samples were collected immediately upon the completion of foraging. Fecal samples were collected at a minimum distance of approximately 1–2 m to avoid individual sampling repetition (Zhang & Zhou, 2012). All samples were obtained from inside the feces to avoid soil contaminants. The samples were kept in a cooler, transported under refrigeration to the lab within several hours and stored at −80 °C. A total of 87 fecal samples were used in this study. Forty-three samples were collected from Shengjin Lake (SJL) and 44 samples from Caizi Lake (CZL). In Shengjin Lake, the samples were 13 in the early wintering period (SJL-E), 14 in the middle wintering period (SJL-M) and 15 in the late wintering period (SJL-L). In Caizi Lake, the samples were 14 in the early wintering period (CZL-E), 15 in the middle wintering period (CZL-M), and 15 in the late wintering period (CZL-L).

DNA extraction and species determination

DNA was extracted from the fecal samples using the Qiagen QIAamp^® DNA Stool Mini Kit following the DNA isolation protocol for pathogens. DNA was eluted in 150 µL of sterilized deionized water and then stored at −80 °C. The COI gene was amplified using the BIRDF1–BIRDR1 primer pair to determine host species. PCR amplification was performed in a 50 µL reaction volume containing 100 ng of fecal DNA, 25 µL of 2 × EasyTaq mix (TransGen), 10 µM BIRDF1 (TTCTCCAACCACAAAGACATTGGCAC) and 10 µM BIRDR1 (ACGTGGGAGATAATTCCAAATCCTG). The cycling conditions were as follows: 95 °C for 5 min, followed by 30 cycles of 95 °C for 30 s, 55 °C for 45 s and 72 °C for 90 s, with a final extension period at 72 °C for 10 min. The PCR products were sequenced by Sangon Biotech Co. Ltd. (Shanghai). The sequences were aligned by NCBI. All of the samples were confirmed to contain hooded cranes DNA based on the sequence analysis.

16S rRNA sequencing

Primer sets 338F/806R equipped with sequencing adapters and unique identifier tags were used to amplify the V3-V4 hypervariable regions of the bacterial 16S rRNA genes fragments for the Illumina Hiseq platform (PE250). PCR was carried out in 50 µL volume containing 60 ng of fecal DNA, 25 µL of 2 × Premix Taq and 1 µL each of the forward and reverse primers (10 µM). The cycling conditions were as follows: 94 °C for 3 min, followed by 30 cycles of 95 °C for 30 s, 55 °C for 45 s and 72 °C for 45 s, with a final extension at 72x °C for 10 min. Triplicate reaction mixtures per sample were pooled together and purified using the EZNA Gel Extraction Kit (Omega, Norcross, GA, USA). Quantification was performed with a NanoDrop. Sequences were processed using Quantitative Insights Into Microbial Ecology (QIIME). After denoising, poor quality (below an average quality score of 30) and short sequences (shorter than 200 bp) were removed and those sequences which exactly matched with their barcode and primer were kept. The barcode and primer sequences were then removed. Sequences were clustered into OTUs using a 97% identity threshold by UPARSE. Singleton OTUs were deleted and chimers were filtered to remove erroneous OTUs due to sequence errors using the USEARCH in QIIME. The most abundant sequence within each cluster was selected as the representative sequence for that OTU. The taxonomic identify of each OTU was checked against the ribosomal database project (RDP) Classifier. The raw data has been submitted to the SRA database of NCBI with accession number SRP095274.

Data analysis

Non-metric multidimensional scaling (NMDS) based on the Bray–Curtis dissimilarity (calculated from the relative abundance matrix) and Analyses of Similarities (ANOSIM; permutations = 999) were performed to compare the community composition in different treatments in R 3.4.1 (vegan 2.4-3) (Dixon, 2010). In addition, permutational analysis of multivariate dispersions (PERMDISP; permutations = 999) was used to test the gut bacterial community structure heterogeneity varied over spatial–temporal scale. This test was using a Bray–Custis similarity matrix. The significance of PERMDISP was determined via the ANOVA F-statistic to compare among-group differences in the distance from observations to their corresponding group centroids. One-way and two-way ANOVA was used to test observed species, PD whole tree, Shannon and Chao 1 indices among different treatments. The nearest taxon index (NTI) and betaNTI are used to test the assembly processes of the gut bacterial community. NTI was calculated to assess the spatial and temporal changes in bacterial phylogenetic structure using Picante software. The NTI can be used to test for phylogenetic clustering or overdispersion. Positive NTI values and low quantiles (P < 0.05) indicate that co-occurring species are more closely related than expected by chance (clustering), whereas positive values and high quantiles (P > 0.05) indicate that the co-occurring species are less closely related than expected by chance (overdispersion). Positive betaNTI values indicate greater than expected phylogenetic turnover, and betaNTI was calculated in phylocom. LEfSe (linear discriminant analysis effect size) was used to identifies genomic features characterizing the differences between two or more biological conditions (Segata et al., 2011). To detect KEGG pathways with significantly different abundances between the two lakes, LEfSe analysis was used according to the online protocol (http://huttenhower.sph.harvard.edu/galaxy/). Nearest sequenced taxon index (NSTI) was the sum of phylogenetic distances for each organism in the OTU table and the closest genetic relationship of the sequencing reference genome, as measured by the substitutions per site in the 16S rRNA gene and weighted by the frequency of that organism in the OTU table (Langille et al., 2013). Functional predictions were made based on the 16S rRNA OTU membership using PICRUSt (Phylogenetic Investigation of Communities by Reconstruction of Unobserved States) according to the online protocol (http://picrust.github.io/picrust/). One-way ANOVA was used to detect the influence of the sampling site and season on bacterial taxonomy (phylum) variation. All univariate statistical analyses were conducted using SPSS 20.0.

Bacterial alpha-diversity

Gut bacterial alpha-diversity was calculated at a depth of 10,000 randomly selected sequences per sample. Alpha-diversity was estimated by the observed species index, phylogenetic diversity, Shannon and Chao 1 index. Two-way ANOVA showed that all alpha diversity was significantly different with time, and only phylogenetic diversity and Chao 1 index were significantly different in spatial variation. Alpha diversity was significantly different under the interaction of wintering period and spatial distance (Table S1). At the early wintering periods, the alpha-diversity of the SJL samples were significantly higher than that of the CZL samples (P < 0.001) (Fig. 2), but there was no significant difference in middle and late wintering periods between the two lakes (P > 0.05 in both causes). Alpha-diversity was the highest in the middle period across the temporal changes in CZL, and the Chao 1 index of the SJL samples showed significant different across the temporal change (Fg. S1).

Bacterial community composition

In total, we obtained 4,717,488 reads of high-quality sequences with an average of 46,193 reads was found. Potentially chimeric sequences and singleton OTU (698,694) were then discarded. The dominant gut bacterial phyla were Firmicutes (59.82% ± 32.35%), Proteobacteria (26.82% ± 24.87%), Actinobacteria (5.43% ± 8.19%), Fusobacteria (3.78% ± 12.22%), and Bacteroidetes (2.24% ± 6.21%) (Figs. 3A, 3C). Within Firmicutes, the dominant bacterial classes were Clostridia (41.79 ± 31.94%) and Bacilli (17.96 ± 21.17%). The classes of Epsilonproteobacteria (12.61 ± 21.80%), Gamaproteobacteri a (8.69 ± 13.87%) and Alphaproteobacteria (4.28 ± 8.04%) were dominant in Proteobacteria (Figs. 3B, 3D). However, the distribution of each taxon among the four groups was uneven, as indicated by Fig. 3. At the family level, the top 30 families were shown for all samples (Fig. S2). At the lower level, only 84.7% of sequences could be assigned to 995 genera. The dominate genera were Clostridium (15.5%), SMB53 (13.62%), Helicobacter (11.86%), Lactobacillus (7.22%), Epulopiscium (4.82%), Enterococcus (4.76%), Fusobacterium (3.67%), Pseudomonas (2.54%), Turicibacter (2.08%), Serratia (2.03%), Agrobacterium (1.77%) and Lysinibacillus (1.59%).

One-way ANOVA showed that relative abundance of Firmicutes was significantly higher in SJL samples than in CZL samples, whereas the relative abundance of Proteobacteria from the SJL samples was significantly lower than the CZL samples (Fig. 4, Fg. S3). At Shengjin Lake, the relative abundance of the Firmicutes increased significantly along the winter period. In order to explore the differences in spatial and temporal scale, we conducted LEfSe tests to detect the difference in relative abundance of microbial taxa. At Shengjin Lake, four indicator bacterial taxa (i.e., Pseudomonadaceae, Pseudomonas, Sphingobacteriales, Sphingobacterii) were found in the gut of the hooded crane in SJL-E samples and five indicator bacterial taxa (i.e., Clostridiales, Clostridia, Lachnospiraceae, Epulopiscium, Peptostreptococcaceae) were found in the gut of hooded crane in SJL-L samples (Figs. 5A, 5B). At Caizi Lake, three indicator bacterial taxa (Fusobacteriaceae, Fusobacterium, Peptostreptococcaceae) were found in the CZL-E samples, four indicator bacterial taxa (i.e., Sphingobacteriaceae, Sphingobacterii, Sphingobacteriales, Microccaceae) were found in the CZL-M samples and 11 indicator bacterial taxa (Gammaproteobacteria, Enterobacteriaceae, Enterobacteriales, Bacilli, Lactobacillales, Enterococcaceae, Enterococcus, Lactobacillales, Enterococcus, Lactobacillus, Lactobacillaceae, Enterobacter) were found in the CZL-L samples (Figs. 5C, 5D).

The bacterial community compositions were significantly different between the guts of hooded crane from SJL samples and CZL samples. NMDS revealed that SJL samples tended to be less different compared to the CZL samples (Fig. 6E). ANOSIM analysis confirmed that seasonal changes had a significant impact on the microbial composition of the hooded cranes in two lakes (Table 1). The bacterial community composition of SJL samples and CZL samples were significantly different (P = 0.001) in the three wintering periods (Table 1), the bacterial community compositions in different wintering period were significantly different in both Shengjin and Caizi Lake (R = 0.238, P = 0.001) (Table 1).

According to PERMDISP, average distance to the corresponding group centroid have no significant differences among different wintering periods at Shengjin Lake (P = 0.167). However, mean distance to the corresponding group centroid have significant differences among different wintering periods at Caizi Lake (P < 0.001). Furthermore, no significant differences could be detected among average distance to the corresponding group centroid of gut bacterial communities in different sample collected sites.

Assemblage processes of the gut bacterial community

NTI was used to evaluate the gut bacterial phylogenetic structure. All of the NTI values were positive, which showed that the bacterial communities were phylogenetically clustered (Fig. 7, Fig. S4). At the early wintering stage, the CZL-E samples NTI values were less positive compared to the SJL-E samples, which indicated that phylogenetic clustering was weakest in the CZL-E samples (Fig. 7). In addition, the NTI values of SJL-L samples were less positive compared to the CZL-L samples at late wintering stage, which indicated that phylogenetic clustering was weakest in the SJL-L samples. Phylogenetic clustering was similar in SJL samples during three wintering periods, whereas CZL-L samples phylogenetic clustering were more similar in CZL samples.

Variation in predicted metagenomes between the two lakes

In addition, NSTI also influenced PICRUSt accuracy. The NSTI for our samples was 0.18 ± 0.07 similar to the previously reported analyses in soils (mean NSTI = 0.17 ± 0.02). In this study, a total of 41 functional genes were predicted in hooded crane population. The majority of functions were membrane transport (12.63%), carbohydrate metabolism (9.36%), amino acid metabolism (9.19%), replication and repair (8.08%), energy metabolism (7.10%) and translation (5.18%) during the wintering period. ANOSIM analysis revealed that the potential functions of the gut bacterial communities of the two lakes were significantly different during early and late wintering periods (P = 0.001) (Table 2). However, during the middle wintering period, there was little difference in the potential functions of the gut bacterial in the two lakes (P = 0.116). In early wintering period, hooded crane from Caizi Lake were predicted to have a microbiota with enhanced capability for metabolic disease, carbohydrate metabolism, amino acid metabolism, circulatory transport, cellular processes and singling and a lower rate of microbial genes associated with genetic information processing (Table S2). In middle wintering period, hooded crane from Caizi Lake were predicted to have a microbiota with enhanced capability for amino acid metabolism, sensory system, transport and catabolism, membrane transport and genetic information processing. When comparing predicted metagenomes in the late period, the microbiota of hooded crane from Caizi Lake was predicted to have a capability for metabolic disease, carbohydrate metabolism, amino acid metabolism, replication and repair and glycan biosynthesis and metabolism.