Comparative Analysis of the Duodenum Microbial Communities of Zhedong White Goose (Anser cygnoides) at Different Reproductive Stages by High-Throughput Sequencing

Wanqiu Zhao1,2, Taoyan Yuan1, Li Chen2, Xue Du2, Weihu Chen4, Yan Fu1, Dong Niu3,* and Lizhi Lu2,* 1College of Animal Sciences, Zhejiang University, Hangzhou 310058, China 2Institute of Animal Husbandry and Veterinary Science, Zhejiang Academy of Agricultural Sciences, Hangzhou 310021, China 3College of Animal Science and Technology, Zhejiang A&F University, Key Laboratory of Applied Technology on Green-Eco-Healthy Animal Husbandry of Zhejiang Province, Hangzhou 311300, China 4Xiangshan County Agricultural and Rural Bureau, Ningbo 315700, China Article Information Received 01 November 2019 Revised 01 March 2020 Accepted 14 March 2020 Available online 10 December 2020


INTRODUCTION
G oose is an economically important herbivorous waterfowl, supplying humans with eggs, nutritious meat, as well as down and feathers. The Zhedong White goose, which is characterized by strict seasonality and a high tendency to broodiness and incubation, as well as a low rate of laying, is one of the typical seasonal reproductive poultry in China (Zhao et al., 2013). Their regular reproductive season usually starts in September and lasts until the May of the following year, including 3 or 4 laying cycles, and each cycle contains three stages of laying, broody, and recovery. The food intake, metabolism, and neuroendocrine vary at each stage. During the broodiness of female geese, the behavioral tendency to sit on the nest and incubate eggs is considered to be largely responsible for their low productivity. This behavior always accompanies by a decreased appetite. We thus hypothesize that the gut microbiota differs at different reproductive stages and is associated with the reproductive behaviors of geese.
Recently, there have been increasing studies on the intestinal microorganisms of poultry due to the availability of enhanced technologies. For instance, the bacterial community in Turkey feces was identified to be dominated by low G+C Gram-positive bacteria using a combination of 16S rRNA gene and metagenome approach (Lu and Domingo, 2008). In a similar way, Lim et al. (2015) conducted sequencing of 16S rRNA libraries of chicken fecal samples at two developmental time points to investigate the microbial population dynamics. The differences of the microbial community structures and compositions of functional genes between wild geese and artificially-bred geese were thoroughly described,

O n l i n e F i r s t A r t i c l e
including that the wild geese had a significantly higher abundance of Proteobacteria, while Bacteroidetes was more abundant in the artificially-bred geese (Wang et al., 2016(Wang et al., , 2017. A recent study has found that Firmicutes and Bacteroidetes were more abundant respectively in the geese fed with all-grass and high-grain diets (Xu et al. 2017). It was also revealed that the dominant microbiota communities mainly involved in the carbohydrate metabolism in all-grass-fed geese (Xu et al., 2017). Gut bacteria also have been shown to play a critical role in multiple physiological changes related to health and traits of the host, such as obesity (Cho et al., 2012), gut immune maturation (Chung et al., 2012), metabolic homeostasis, and reproductive performance (Gioacchini et al., 2010;Dai et al., 2015). Located at the beginning of the intestine, the duodenum is crucial for food digestion and absorption, has a lower pH than the hindgut, and is the region that absorbs most glucose and other nutrients within the small intestine (Muir and Hopfer, 1985;Heard and Annison, 1986). Accordingly, the present study aimed to investigate, for the first time, the effect of reproductive stages on the intestinal microbiota of Zhedong White goose and attempted to correlate the changes in microbial communities' profiles to the production performance of geese. We analyzed and compared the duodenum microbial communities of geese at different reproductive stages by high-throughput sequencing of the 16S rRNA gene with an Illumina Hi-seq platform.

Ethics statement
This study was performed according to the recommendations of the Guide for the Care and Use of Laboratory Animals of the National Institutes of Health, and all the protocols were approved by the Institutional Animal Care and Use Committee (IACUC) of Zhejiang Academy of Agricultural Sciences and Zhejiang University.

Experimental design and sample collection
The experiment was conducted in the Xiangshan Goose Breeding Farm, Zhejiang Province, China. Experimental geese were approximately 44 weeks old, from the same hatch, and were raised according to the farms standard practice: had free access to feed, given water and green forages ad libitum, at the same time, supplemented with a basal goose diet comprised 6%~10% soybean meal, 50% corn and 40% milled wheat, and exposed to natural light and temperature throughout our study. We collected the duodenum contents from a total of 15 Zhedong White geese (5 replicates for each group) at three different reproductive stages (laying, broody, and recovery) via exsanguination and necropsy. The body weight of laying geese (5.356±0.51 kg) was significantly higher (P<0.05) than that of broody geese (3.978±0.31 kg) and recovery geese (3.944±0.49 kg). The intestinal samples collected in sterile tubes were immediately transferred into liquid nitrogen for temporary storage and then transported to the laboratory where the samples were kept at -80°C until processed.

DNA extraction and 16S rRNA sequencing
DNA was extracted from the duodenal content samples using a QIAamp ® Fast DNA Stool Mini Kit (QIAGEN, Germany) according to the manufacturer's instructions. DNA quality and concentration were assessed on 1% agarose gel and then diluted to 1ng/μl using sterile water. The V4-V5 region of the bacterial 16S rRNA gene was amplified using the primer set of 515F (5'-GTGCCAGCMGCCGCGGTAA-3') and 907R (5'-CCGTCAATTCCTTTGAGTTT-3'), with the reverse primer containing an 8-bp barcode sequence unique to each sample. PCR reactions were performed in 20μl mixture containing 1μl DNA template, 0.5μl of each primer (10μM), and 10μl Phusion ® High-Fidelity PCR Master Mix (New England Biolabs) with the following conditions: 94°C for 3 min, then 30 cycles of 98°C for 10 s, 55°C for 30 s, 72°C for 30 s, and a last extension step of 72°C for 5 min. The PCR products were separated by electrophoresis in a 2% agarose gel and purified with a Qiagen Gel Extraction Kit (QIAGEN, Germany).
Subsequently, we conducted library sequencing using a TruSeq ® DNA PCR-Free Sample Preparation Kit (Illumina, San Diego, California, USA) and evaluated the quality on a Qubit@ 2.0 Fluorometer (Thermo Scientific) and an Agilent Bioanalyzer 2100 system. Finally, the samples were sequenced on an Illumina Hiseq 2500 platform according to the standard protocols that generated 250bp paired-end reads. The sequencing data were deposited into the Sequence Read Archive (SRA) of NCBI (https://www.ncbi.nlm.nih.gov/sra) and can be accessed via accession number SRP219168.

Assembly and quality control of paired-end reads
In order to obtain precise and reliable results, we preprocessed and demultiplexed the original sequencing reads with the following procedure: i) the paired-end reads were assigned to each sample according to their unique barcodes and primer sequences; ii) the barcodes and primer sequences were trimmed; iii) the paired-end reads were merged using a fast and accurate analysis tool-FLASH (V1.2.7; http://ccb.jhu.edu/software/FLASH/) (Magoc and Salzberg, 2011), which can assemble reads by O n l i n e F i r s t A r t i c l e their overlapping regions, and the merged sequences were called raw tags. After initial trimming, additional quality controls were performed by applying the QIIME (V1.7.0; http://qiime.org/scripts/split_libraries_fastq.html) (Caporaso et al., 2010) software package (Quantitative Insights into Microbial Ecology) to filter tags containing ≥5 Ns and low-quality tags. Next, the chimeric sequences were removed, and the effective tags were acquired for the downstream analysis.
Operational taxonomic units (OTUs) clustering and species annotation UPARSE (v7.0.1001; http://drive5.com/uparse/) was used to cluster the effective reads with more than 97% similarity into operational taxonomic units (OTUs) and to screen the representative sequences of each OTU, which were used to annotate taxonomic information utilizing the MOTHUR and SILVA (http://www.arb-silva.de/) SSU database (Wang et al., 2007;Quast et al., 2013). The taxon abundance of each sample at phylum, class, order, family, and genus levels was displayed in histograms and heatmaps.
In order to study phylogenetic relationships among different OTUs and the differences of the dominant species in the different samples (groups), a multiple sequence alignment was conducted using the MUSCLE software (version 3.8.31; http://www.drive5.com/muscle/) (Edgar, 2004). OTU abundance information was normalized using the sequence number corresponding to the sample with the least sequences. The subsequent analyses of alpha diversity and beta diversity were all performed based on the outputted normalized data.

Alpha diversity analysis
Alpha diversity analysis includes three metrics: Chao 1, Shannon index, and Observed Species, which estimate the community richness, the community diversity, and the amount of unique OTUs found in each sample, respectively. All these indices of our samples were calculated with QIIME (version 1.7.0), and Rarefaction curves were generated based on these three metrics using R software (version 2.15.3)

Beta diversity analysis
For the beta-diversity analysis, which was used to evaluate differences of samples in species complexity, both the weighted and unweighted unifrac were calculated by QIIME software (version 1.7.0). We used unweighted unifrac for principal component analysis (PCA), principal coordinate analysis (PCoA), non-metric multi-dimensional scaling (NMDS), and unweighted pair group method with arithmetic mean (UPGMA) clustering.  Tags) show the number of tags that could not be clustered to the OTUs; the blue bars (Taxon Tags) show the number of tags clustered to the OTUs that had the annotation information in each sample; the purple bars (OTUs) represent the number of OTUs obtained from each sample. L1~5, B1~5, and R1~5 represent 5 biological replicates at egg-laying, broody and recovery stages, respectively.  Raw PE, original paired-end reads obtained from the sequencing platform; Raw tags, merged paired-end reads; Clean tags: tags after filtering low-quality or short sequences from raw tags; Effective tags, tags after removing the chimera sequences from the clean tags; AvgLen, the average length of the effective tags; Q20/Q30, the percentage of bases with quality higher than 20/30 in the effective tags; GC%, the percentage of GC base in the effective tags; Effective%, The percentage of Effective tags in Raw PE; L1~L5, B1~B5, and R1~R5 represent the five biological replicates of the laying, broody, and recovery groups, respectively.
Analysis of similarity (ANOSIM) was performed to test significant differences in community structures among the groups. The two-sided Student's t-test was used to identify the different bacterial phylum. Biomarker discovery analysis was achieved through Linear Discriminant Analysis Effect Size (LEfSe). All the statistical analysis described above were conducted with R package vegan (version 2.15.3).

A general view of the 16S rRNA sequencing data
In the present study, we characterized the duodenum bacterial communities of the geese via 16S rRNA sequencing and obtained totally 1,167,120 raw paired-end reads from 15 samples. After data preprocess and quality control, including trimming the barcodes and primers, filtering low-quality reads and chimaeras by QIIME, finally a total of 764,344 effective tags with Q20>98% and Q30>96% were obtained from the 15 samples. The average number of the effective tags for each sample was 50,956 (the minimum one was 34,379, and the maximum one was 60,597), and the average effective rate was 65.3%. The detailed statistical information of each sample is listed in Table I. As shown in Figure 1, we totally identified 4,748 OTUs (average=317) from all the samples. The effective tags with ≥ 97% sequence similarity were assigned to the same OTUs. The Venn Graph ( Supplementary Fig. S1) showed that the laying group had the most specific OTUs, while the three groups had 263 shared OTUs.

Taxonomic compositions of goose gut microbial communities at the three reproductive stages
According to the OTUs annotation results produced by MOTHUR and SILVA SSUrRNA database, we calculated the relative abundance (phylum, class, order, family, and genus levels) of all the samples. The taxonomic distribution of each group at the phylum level is shown in Figure 2A.

F i r s t A r t i c l e
Comparative Analysis of the Duodenum Microbial Communities of Zhedong White Goose 5 Firmicutes in the laying group was significantly increased (P<0.05), and the Actinobacteria also was more abundant, but did not reach a significant level (P>0.05). Meanwhile, the abundance of Cyanobacteria in the broody group and Spirochaetes in the recovery group were specifically higher (P<0.05) than the corresponding one in the laying group. No significant differences were observed in the proportions of Proteobacteria in the goose duodenum at three reproductive stages. At the genus level, based on the filtering criteria of each genus (more than 0.005% of the total sequences) (Bokulich et al., 2013), the sequences of the 15 duodenal content samples represented 304, 216, and 233 genera, respectively (Supplementary Table I). The top 10 genera were listed in Table II. The average abundance of these 10 major genera in the laying group geese was 67.59%, ranging from 54.34% to 78.57%, and was 83.19% in the broody group geese, ranging from 75.05% to 93.36%. Simultaneously, in the recovery group samples, these genera contributed an average of 80.82% of the total microbial abundance, ranging from 42.92% to 92.04%. Turicibacter, unidentified-chloroplast, and Unidentifiedspirochaetes were the most abundant genera in the laying, broody, and recovery group, respectively. Moreover, we also found that the genera of Turicibacter, Lactobacillus, Trichococcus (phylum Firmicutes), and Brevibacterium (phylum Actinobacteria) were more abundant (P>0.05) in the laying group geese than other two groups, while unidentified-chloroplast (phyla Cyanobacteria) and unidentified-mitochondria (phylum Proteobacteria) in the broody group and unidentified-Spirochaetes (phylum Spirochaetes) in the recovery group were significantly more enriched (P<0.05) than those genus in the laying group. Furthermore, the clustering heatmaps of taxa were also constructed to detect the differences of bacterial community compositions as displayed in Figure 2B. It is evident that taxonomic distributions (phylum level) among the groups were significantly different.  . Differences in bacterial community richness and diversity among the three groups. Asterisk denote significant differences between groups. Significant differences are defined at the 95% confidence level. The Chao1 index of broody group was significantly lower than laying group; no significant difference was detected in Shannon index.

Differences in diversity of bacterial communities
Alpha diversity analysis was applied to describe the diversity of microbial communities within the samples.
The alpha diversity indices of four aspects were shown in Supplementary Table II, including community richness (Chao1, ACE, and observed species), community diversity (Shannon, Simpson), sequencing depth (Good's coverage), and phylogenetic diversity (PD whole tree). The Chao1 value of the laying geese was higher than the other two groups, and the difference between the laying and broody geese is statistically significant (P<0.05), but the difference between the laying and recovery geese is not significant (Fig. 3A). On the other hand, no significant differences were observed in the Shannon index among the three groups as shown in Figure 3B. Rarefaction curves were created with R software (Fig. 4). Using a rarefaction depth of 25,252 sequences per sample, the number of OTUs have almost approached a plateau except two samples from the laying group, meaning that although the depth was enough to draw our conclusion, it will be better to have a deeper sequencing to fully estimate the bacterial OTU diversity in these samples.

Differences in microbial community compositions
To compare the duodenal microbial community compositions of the geese at different reproductive stages, PCoA analysis was carried out. Figure 5A and B showed the PCoA plots based on the unweight and weighted Unifrac distance matrixes. Although exceptional samples (e.g. L4, L5, B3, R3) clustered separately from other members within their groups, the replicates of the same group tended to gather closer in general. Furthermore, the Anosim test supported the between-group distances were higher than the within-group distances (R>0), and significant differences O n l i n e

F i r s t A r t i c l e
Comparative Analysis of the Duodenum Microbial Communities of Zhedong White Goose 7 in the community structures were observed among the three groups (P<0.05, Supplementary Table IV). These data suggest that the microbial communities of the geese at the three reproductive stages were significantly different. Fig. 5. Principal Coordinate Analysis (PCoA) plots of the duodenal bacteria community structures of the three groups using unweight (A) and weight (B) UniFrac distance metrics. Red, green and blue colors represent the laying, broody and recovery samples, respectively. PC1 and PC2 represent two principal components, and the percentage represents the contribution of the principal component to the sample difference. Samples with high community structure similarity tend to cluster together, while samples with large community differences tend to be far apart.

Analysis of differences in bacterial taxa among the groups
We also performed LEfSe (LDA score≥4) analysis to detect specific species whose abundance was statistically different among the groups. Figure 6A-C showed that 10 bacterial taxa were significantly more abundant in the laying group geese (e.g., Firmicutes, Actinobacteria, Micrococales, P<0.05), while only 5 taxa were overrepresented in the broody group geese (e.g., Mitochodria, Rickettsiales, P<0.05). Meanwhile, 6 taxa were more enriched in the recovery group geese (e.g., Spirochaetes, Spironema_culicis, P<0.05).

DISCUSSION
In this study, we characterized the changes in duodenal bacterial communities during the laying, broody and recovery stages of Zhedong White geese by the high-throughput sequencing technology. Our results have shown that the representative taxonomic phyla in the laying, broody, and recovery geese were Firmicutes, Cyanobacteria, and Spirochaetes, respectively. Furthermore, the microbial community diversity and richness of the laying geese showed a tendency to increase compared with the other two groups, supported by the higher values of OTUs, Chao1, and Shannon indices in this group of geese. Several studies revealed that a rich bacterial community is associated with a healthy and productive status, while a deficient microbial community is associated with several disorders of metabolic and physiological functions (Kosiewicz et al., 2011).
Comparison of the gut microbiota indicated that Firmicutes were the most prominent phylum in the laying group with an abundance of 52.29%, almost twice of the other two groups. Numerous studies have demonstrated that Firmicutes primarily predominate the bacterial profile in various animals, including human (Turnbaugh et al., 2006), turkeys (Lu andDomingo, 2008), broiler chickens (Lim et al., 2015;Wang et al., 2018), and geese (Wang et al., 2016). Compared with low egg-laying hens, the hens of high egg-laying performance had a significantly higher relative abundance of Firmicutes, whose potentially useful effects on diets fermentation and digestion led to an increase in small water-soluble food molecules for absorption into the watery blood plasma (Delzenne and Cani, 2011;Elokil et al., 2019). Firmicutes is considered to play a major role in the cellulose decomposition and is supposed to be involved in the process of energy harvest and nutrient uptake from the feed (Daly et al,. 2001). Additionally, Firmicutes also has been associated with short-chain fatty acid metabolism, which contribute to the synthesis propitiate and butyrate (Polansky et al., 2015). Therefore, we speculate that a high concentration of Firmicutes present in the gut of the laying geese may significantly promote digestive efficiency and assimilation of feed energy, which may contribute to providing energy and nutrients for egg-laying activities.
Furthermore, the detailed compositions of phylum Firmicutes at the genus level in the three groups were not the same. For the laying geese, the gut microbiota was characterized by an unknown genus Turicibacter and genus Lactobacillus. As reported in recent 16S rRNA analysis studies, researchers detected the presence of Turicibacter in the gastrointestinal (GI) tracts of humans (Cuiv et al., 2011), pigs (Rettedal et al., 2009), rats (Licht et al., 2007, and goats (Liu et al., 2014), implying this genus might be an important member of the gut microbiota. The weaning piglets fed on a diet with chlortetracycline, which was used to promote growth through increasing feed intake, had a decrease in Turicibacter (Rettedal et al., 2009). In contrast, the high-grain feed increased the abundance of Turicibacter and concurrently resulted in caecal mucosal epithelial damages in goats (Liu et al., 2014). Consequently, they speculated that the Turicibacter bacteria present in animals may cause subclinical infections or have some other deleterious effects on the GI tract. Nevertheless, our results showed Turicibacter had a high proportion (16.2%) in the laying group, consistent with what Wang et al.
(2016) discovered in the Bar-Headed geese. The potential effects of Turicibacter on goose health or reproduction remain to be elucidated by further studies.
On the other hand, Lactobacillus bacteria are known to produce a range of bacteriocins, and many Lactobacillus isolates have been shown as probiotics (Parvez et al., 2006;Kim et al., 2007). Some Lactobacillus species are considered to be harmless and are possibly effective in preventing bacterial infections, as they can generate lactic acid, which may produce an acidic environment that could restrain the growth of many pathogenic genera (Kim and Isaacson, 2015). They tend to have a beneficial effect on broiler performance, including the modulation of intestinal microflora, pathogen inhibition, intestinal historical changes, immunomodulation, certain haemato-biochemical parameters (Kabir et al., 2005). Previous studies also suggested that the high egg-laying performance in hens was due to the increasing abundance of Lactobacillus, which were growth promoters and have antimicrobial against pathogenic microbes (Choe et al., 2012;Elokil et al., 2019). The Lactobacillus has been studied and used in medicine and the food industry for years. It is reported that Lactobacillus was highly related with the host feed efficiency, as more Lactobacillus enriched in better feed efficiency group than poor feed efficiency group, which could generally improve the gastrointestinal tract and thus protect the gut from pathogens and promote O n l i n e

F i r s t A r t i c l e
Comparative Analysis of the Duodenum Microbial Communities of Zhedong White Goose efficient nutrient and energy extraction in the host (Yan et al., 2017). There was more Lactobacillus in the duodenum of the egg-laying group than the other two groups, suggesting a higher digestion level in the laying geese. It is noteworthy that the fecundity significantly increased in zebrafish after probiotic Lactobacillus rhamnosus administration, which may act indirectly by activating a potent metabolic hormone, leptin (Gioacchini et al., 2010). Hence, we propose that the Lactobacillus bacteria may play a positive role in the egg-laying process, but more detailed investigations of the specific bacterium are still needed. An increase in the abundance of Actinobacteria was detected in the laying geese compared to the other two groups, which largely consisted of genus Brevibacterium. Brevibacterium, a Gram-staining-positive, aerobic, non-motile bacterium, has been isolated from diverse habitats, such as human skin (Roux and Raoult, 2009) and poultry manure (Tonouchi et al., 2013). It was found to be associated with the metabolism of aromatic chemicals (Rattray and Fox, 1999). In human, Actinobacteria had a higher proportion in the microbiome of obese persons compared with their lean twins. Among the obesityenriched genes (involved in carbohydrate, lipid, and amino acid metabolism), 75% of them were from Actinobacteria (Turnbaugh et al., 2008).
We particularly identified that Cyanobacteria, a phylum that is not generally considered to be gut microbes, were remarkably dominant in the broody geese. It was mainly represented by the genus of an unidentified-Chloroplast. Previous studies also have observed that Cyanobacteria was more abundant in the feces of low egg-laying hens than that of high egg-laying hens (Elokil et al., 2019). Ley et al. (2008) detected the deep-rooting Cyanobacteria in mammal gut microbiotas, and this phylum may represent the descendants of the non-photosynthetic ancestor of cyanobacteria, which have adapted to the life in animal GI tracts. Similarly, Liu et al. (2014) in his article also introduced the appearance of Cyanobacteria in the caecal luminal or mucosal samples of high-grain fed goats, indicating that environmentally resistant organisms reside within the caecum of goats. Likewise, Cyanobacteria were also detected as a dominant phylum in the porcine ileum and showed a significant increase as the dietary protein content reduced (Qiu et al., 2018), which is in accordance with the previous report that Cyanobacteria increased when fishmeal was progressively replaced by soybean meal in diets (Parma et al., 2016). Geese are herbivorous waterfowl. Therefore, it is reasonable to infer that the observed Cyanobacteria come from the plants (e.g., water hyacinth, maize) in the diet, similar to previous studies (Oakley and Kogut, 2016;Rimoldi et al., 2018). During the broody stage, the goose has reduced activities and energy demands, which leads to the decreased Firmicutes abundance and low digestion level. All of these might be responsible for the accumulation of large number of Cyanobacteria bacteria in the broody geese.
An overwhelming percentage of Spirochaetes was observed in the recovery geese, but this phenomenon has not been found in other studies on poultry, confirming again the finding that there can be a high variation in the microbiota composition among flocks (Stanley et al., 2013). The spirochetes are free-living or host-associated helical bacteria, some of which are pathogenic to men and animals (Paster and Dewhirst, 2000). In the current study, the major species belonging to this phylum was Spironema culicis, a new spirochete isolated from the mosquito, and its pathogenicity is unknown. The geese were healthy when they were sacrificed for the sample collection. The reason for the Spirochaetes increase in the recovery geese is not clear and needs to be further studied.
In summary, this is the first description of the Zhedong White goose duodenum microbial community using 16S rRNA sequencing profiling. Comparative analysis identified differences in the structures and functions of gut microbial communities among the laying, broody, and recovery geese, which may be related to the changes in the physiological state and dietary level of the geese. Our results suggest that the digestion process may differ in the laying and broody geese, resulting in changes in the microbiota composition. To obtain a more reliable result and fully understand the specific bacteria and their effects on the host, we need to conduct further largesample experiments or adopt other approaches to verify our findings here.

O n l i n e F i r s t A r t i c l e
W. Zhao et al.

Statement of conflict of interest
We certify that there is no conflict of interest with any financial organization regarding the material discussed in the manuscript.   O n l i n e O n l i n e