The experiments were authorized by the Animal Care and Use Committee of the Zhejiang Academy of Agricultural Sciences (2019ZAASLA37).

A population of newly-hatched 5000 Muscovy (C. moschata) ducks (Lanxi Hewang Breeding Duck Co., Ltd., Jinhua, China) were fed with a commercial diet and water ad libitum under standardized conditions for 70 days. The diet composition was as previously described (24). Two hundred ducks were randomly selected from a 5000 population. After being weighed, ducks were euthanized by cervical dislocation. The ileal segment and content of each duck were gathered immediately, frozen in liquid nitrogen, and stored at -80°C until the mRNA and DNA were separated.

Host genomic DNA was separated from each ileum sample using a QIAamp DNA Stool Mini Kit (Qiagen, Valencia, CA) following the manufacturers’ instructions. Use 1% agarose gel electrophoresis and NanoDrop ND-1000 (Thermo Fisher Scientific, Waltham, MA, USA) to evaluate the quality and concentration of the DNA extract. Use next-generation sequencing to sequence high-quality DNA (25). In detail, the V4-V5 region of the bacterial 16S rRNA gene was amplified by the barcode-fusion forward primer 515F (5’-GTGCCAGCMGCCGCGGTAA-3’) and the reverse primer 907R (5’- CCGTCAATTCMTTTRAGTTT -3’). The conditions of PCR were as previously described (26). After PCR, the 2% (w/v) agarose gel was used to separate and qualify the amplicons and then purified using a GeneJET Gel Extraction Kit (Thermo Scientific). An Illumina TruSeq DNA PCR-Free Library Preparation Kit (Illumina) was appropriated for sequencing library generation. An Agilent Bioanalyzer 2100 System and a Qubit 2.0 Fluorometer (Thermo Scientific) were used to evaluate the quality of the generated library. The qualified library was sequenced commercially using Mingke Biotechnology (Hangzhou) on an Illumina NovaSeq platform, producing a 250 bp paired-end read.

The Illumina paired-end reads were demultiplexed and filtered in Quantitative Insights into Microbial Ecology (QIIME) quality filters (27) to remove low-quality reads, merged into tags by FLASH (28), and assorted each sample according to a unique barcode. Analyze the label of each sample after removing redundancy, and the UPARSE and UCHIME were used to assign unique tags with ≥ 97% sequence similarity to the same operational taxonomic units (OTUs).

The RDP classifier was used to annotate the classification information of the selected OUT (29). Alpha-diversity (Observed Species, Chao 1 estimator, ACE, Shannon, and Simpson indices) was calculated and visualized using GraphPad Prism 8 (GraphPad Software, San Diego, CA, USA). The rooted phylogenetic tree was visualized in the Interactive Tree of Life (iTOL, https://itol.embl.de/; version 5.5).

Total RNA in the ileal tissue was extracted by Trizol^® Plus RNA Purification Kit (Thermo Fisher) in accordance with the manufacturer’s instructions. Briefly, 50-100 mg tissue samples were ground into a powder with a tissue homogenizer, and then transferred to a tube with 1 ml of Trizol and incubated at room temperature for 5 mins. Added 0.2 ml chloroform and shook the tube vigorously by hand for 15 s, then incubated at room temperature for 2-3 mins. After being centrifuged at 12000 × g for 15 mins at 4°C, the supernatant was transferred to a new tube, then RNA was precipitated with an equal volume of 70% ethanol. The quantity and quality of RNA were evaluated by a spectrophotometer (NanoDrop-2000, Thermo Fisher Scientific, MA, USA). cDNA was synthesized by SuperScript™ III First-Strand Synthesis SuperMix (Thermo Fisher). RT-qPCR was performed on an ABI prism 7700 Sequence Detector System (Applied Biosystems, Foster City, CA, USA). The reaction scheme was as follows: pre-incubation at 95°C for 1 min, then performing 40 cycles of denaturation at 95°C for 15 s and annealing at 60°C for 25 s. Then, a melting curve analysis was conducted to affirm the specificity and reliability of the PCR products. Using glyceraldehyde 3-phosphate dehydrogenase (GAPDH) as the internal reference gene, the relative expression of mRNA was calculated by the2^−ΔΔCt method (30). The primers designed in this study are listed in Table 1 .

To determine the relationship between growth-related microbiota and immunity. First, we precluded the taxa present in <30% of Muscovy duck ileum samples to detect which microorganisms were significantly related to BW. All the samples were then successively ranked by host BW and the relative abundance of each microorganism. The highest 20% and lowest 20% of the sorted ducks were regarded as two different groups, and all traits between the two groups were statistically analyzed. Subsequently, in the characterization of body weight-related bacteria (detected in at least 30% of ileum samples), static analysis of the BW between the highest 20% and the lowest 20% relative abundance of common bacteria was performed. Finally, the growth-related microbiota and immune-related values were used for correlation heatmap analysis to explore the relationship between them.

The correlation network analysis was performed based on Spearman’s correlation matrixes among the relative abundance of each genus. Networks were visualized by using Gephi v0.9.2 software (France) (31). The correlation heatmap analysis was performed based on the Pearson algorithm among the relative abundance of each genus and immunity value.

Data are expressed as the mean ± SD. All data analyses were visualized using GraphPad Prism 8 software. Unpaired two-tailed Student’s t-test analyzed the difference between the two groups.t. When the P-value did not exceed 0.05, it was considered significant.

The original sequencing reads of this study have been deposited in NCBI with the deposit number of BioProject PRJNA762153.

A total of 8,360,552 high-quality reads were generated with an average of 41,803 reads in each sample, and were classified into 7035 bacterial OTUs based on 97% sequence similarity. Then these OTUs were assigned into 49 phyla, 136 classes, 332 orders, 572 families, 1323 genus, and 1989 species. We first characterized the bacterial composition in the ileum of Muscovy ducks that may be related to growth.

To reveal the bacterial composition in the ileum of Muscovy ducks, the relative abundance was calculated at both the phylum and genus levels. The top six dominant phyla were Firmicutes, Bacteroidetes, Proteobacteria, Fusobacteria, Campilobacterota, and Planctomycetes, representing 69.84% to 99.95% of the total bacterial population, respectively. Among these phyla, Firmicutes accounted for 1.40% to 98.27% of the bacteria while Bacteroidetes amounted to 0.03% to 78.42% ( Figure 1A ). Planctomycetes was the minor phylum detected at an average of 1.53% in the ileum. Additionally, the minor phyla and unidentified bacteria were denoted as “Others”.

At genus level, ten most abundant genera were Candidatus Arthromitus, Bacteroides, Streptococcus, Vibrio, Romboutsia, Cetobacterium, Clostridium sensu stricto 1, Terrisporobacter, Escherichia-Shigella, and Lactobacillus ( Figure 1B ). Obviously, the ileal microbiota appeared to be very diverse, with the top 10 genera accounting for an average of 59.29%, ranging from 1.94% to 98.11%, in relative bacterial abundance in the ileum. Additionally, these genera had a similar relative abundance. The relative abundance of the top three genera was around 10%, while the abundance of the following seven genera was around 5%.

Next, to reveal the interaction among the common genera in the ileum of Muscovy ducks, we calculated the number of common bacteria genera presenting in 50% of the ileum samples, which accounted for 11.48% of all the genera in the ileum, and generated the phylogenetic tree and the co-occurrence network based on the Spearman’s correlation among the representative bacteria in the ileum. The result showed that among the phyla categorized from the 149 common genera, Firmicutes, Bacteroidetes and Proteobacteria were the top three phyla ( Figures 1C, D ). Meanwhile, among these positively related bacterial genera, most Firmicutes, Bacteroidetes and Proteobacteria formed an extensive co-occurrence network with other bacterial genera, and a few formed several independent, small and stable clusters. Some phyla with a small number of genera, such as Acidobacteria and Actinobacteria, were only associated with a few genera to form a small co-occurrence network. In contrast, Cyanobacteria and Planctomycetata were merely part of an extensive co-occurrence network and were related to multiple genera ( Figure 1D ). These results revealed that the relationship between ileal bacteria was diverse and complex.

To determine the growth performance of Muscovy ducks, we weighed the ducks individually and collected ileal segment and content after dissection. We first tested whether the BW data fit a normal distribution before further data analysis. As expected, the BW gave a classical inverted-bell curve (Gauss Fit, P < 0.0001, Figure 2A ). Next, we set the highest 20% (n = 40) and lowest 20% (n = 40) of the BW-ranked ducks as the high and low groups, respectively, and compared BW between these two groups. Not surprisingly, the result showed a significant difference in BW (3.05 vs. 2.49 kg) between the high and the low groups ( Figure 2B , P < 0.0001).

To detect the relative expression of tight junction proteins and cytokines in the high and low groups, we used RT-PCR to measure the gene expression level of tight junction (Claudin1, Claudin2, ZO-1, ZO-2) and cytokines (IL-1β, IL-2, IL- 10, TNF-α) in the ileum tissue. The results indicated that, compared with the low group, the tight junction proteins expression in the high group were higher, but the difference between them was not significant ( Figure 3A , P > 0.05). Similarly, as shown in Figure 3B , the relative expression levels of the pro-inflammatory factors IL-1β, IL-2 and TNF-α were higher in the low group. In contrast, the expression of anti-inflammatory factor IL-10 was lower in the high group. Of note, the expression of cytokines other than IL-2 is not significantly different in the high and low groups.

To investigate whether gut microbiota contributed to the BW difference between the high and low groups, we collected ileum content from ducks of high and low groups. After isolating genomic DNA, 16S rRNA sequencing was performed to study the ileal microbiota diversity and structure. As shown in Figure 4A , the OTUs obtained from ileal microbiota in ducks of the high group with an average of 728.23 was similar to the reads of ileal microbiota in the low group ducks with an average of 735.01 (P > 0.05). The α-diversity was determined by the Chao, Shannon, and Simpson indexes ( Figure 4 ). The Chao index, which accounts for the richness of a community, showed that ducks in the high group harbored a more prosperous microbial species than the low group ( Figure 4B , P > 0.05). Nevertheless, considering the evenness, there was no significant difference in Shannon and Simpson indexes of ileal microbiota between the two groups ( Figures 4C, D , P > 0.05).

For ileal microbial composition, the top 6 phyla were Firmicutes, Bacteroidetes, Proteobacteria, Fusobacteria, Campilobacterota, and Planctomycetes in both of the high and low groups, which was coincident with the bacterial composition in the whole population of 200 ducks ( Figure S1 ). The relative abundance of the top three phyla, Firmicutes, Bacteroidetes, and Proteobacteria, were similar in the high and low groups, which accounted for 89.05% and 84.77%, respectively. However, Bacteroidetes phyla showed a relatively higher abundance in the high group (28.77%) than in the low group (24.60%). Additionally, ducks in the high group had a relatively lower abundance of Firmicutes, Proteobacteria and Fusobacteria than the low group. At genus level, 7 out of the top 10 abundant genera were shared in the two groups with a pretty different abundance, namely Candidatus Arthromitus, Fusobacterium, Vibrio, Streptococcus, Clostridium sensu stricto 1, Romboutsia, and Bergeyella ( Figure S1 ). The other three genera in the high group were Bacteroides, Enterococcus, and Lactobacillus, while Escherichia-Shigella, Terrisporobacter, and Turicibacter were the other 3 of the top 10 genera in the low group.

To investigate the microbiota composition in the high and low groups, we calculated the relative abundance of each phylum and genus shown in a Sankey diagram ( Figure 5 ). In the figure, the lines of different colors represent different bacterial genera, and the height of the rectangles indicates the relative abundance in the groups (left), phyla (middle), and genera (right). At the phylum level, three predominant phyla were Firmicutes, Bacteroidetes, and Proteobacteria, same with the results of 200 ducks. At the genus level, the top eight genera were Candidatus Arthromitus, Bacteroides, Clostridium sensu stricto 1, Lactobacillus, Streptococcus, Romboutsia, Vibrio, and Escherichia-Shigella. In the top 10 genera, Bacteroides and norank were from the phylum Bacteroides, Candidatus Arthromitus, Romboutsia, Streptococcus, Lactobacillus, Clostridium sensu stricto 1 and uncultured from the phylum Firmicutes, Escherichia-Shigella and Vibrio were from the phylum Proteobacteria. Interestingly, the relative abundance of Bacteroides and Candidatus Arthromitus (16.90% and 19.10%) in the high group of ducks were much higher than the low group (2.76% and 5.57%). In contrast, the relative abundance of Escherichia-Shigella, Romboutsia, Streptococcus, Clostridium sensu stricto 1, and Vibrio in the high group were lower than the low group (1.15% vs. 4.83%, 2.04% vs. 3.06%, 2.69% vs. 12.18%, 2.13% vs. 2.97% and 3.37% vs. 10.33%; Figure S1 ).

To determine the relationship between growth-related microbiota and immunity. Above all, we explored the specific bacteria correlated with BW and identified the taxa of which the abundance was significantly different between the high and low groups. The six BW-related genera were Candidatus Arthromitus, Bacteroides, Faecalibacterium, Streptococcus, Escherichia-Shigella, and Oscillospira ( Figure 6A ). Among these, we found that the abundance of streptococcus (2.69%) and Escherichia-Shigella (1.15%) in the high group was significantly lower than that in the low group (12.18%, 4.83%), indicating that these two genera were negatively correlated with BW. However, the situation with the bacterial genera Candidatus Arthromitus, Bacteroides, Faecalibacterium, and Oscillospira in high and low groups was opposite, which indicated a positive correlation with BW.

Next, we investigated whether the tight junctions and cytokines were associated with BW-related genera. The correlation heatmap was used to show the correlation between them, respectively. As shown in Figure 6B , Escherichia-Shigella was positively correlated with IL-1β (P<0.05). Other BW-related bacteria were related to cytokines, but there was no significant difference. Similarly, we found that BW-related bacteria are related to tight junctions and tend to differ even though they are not significant ( Figure 6C , P<0.05). This result indicated that the bacteria in the ileum might be related to the intestinal barrier and immunity.