All animal handling procedures were approved by the Institutional Animal Care and Use Committee (IACUC) of Sichuan Agricultural University (Chengdu campus, Sichuan, China, Permit No. DKY20170913).

The same batch of male goslings, from the Sichuan Agricultural University Waterfowls Breeding Farm (Ya’an, Sichuan, China), were reared under the same rearing environment with free access to feed and water until 120 days. Afterwards, they were randomly divided into 2 groups: CR and FR. At 270 days old, 8 geese were randomly selected for slaughter from each group. The experimental geese were euthanised by carbon dioxide inhalation and cervical dislocation after approximately 12 h of fasting. The intestinal digesta and mid-ileum mucosa were quickly collected and rapidly frozen in liquid nitrogen. Digesta was used for 16S rRNA-seq and mucosa was used for RNA-seq, and both were stored at −80°C prior to sequencing.

In this study, microbial DNA extraction was conducted using the E.Z.N.A. Stool DNA Kit (Omega Bio-Tek, Norcross, GA). DNA concentration and purity were characterized by NanoDrop spectrophotometer (Thermo Fisher Scientific, Waltham, MA, United States) and electrophoresis on 1% agarose gels. DNA with (OD260/OD280) range from 1.8 to 2.0 and (OD260/OD230) range from 2.0 to 2.5 could be used for subsequent experiments. Based on the concentration, DNA was diluted to 1ug/μL with sterile water. And, the V3−V4 hypervariable region was targeted for amplification, utilizing primers 338-F (5′-ACTCCTACGGGAGGCAGCAG-3′) and 806-R (5′-GGACTACNNGGGTATCTAAT-3′). The PCR reaction was performed on a thermocycling PCR system (Bio-Rad T100, Germany) using high-fidelity polymerase according to the following procedure: 98°C for 60 s; 30 cycles of 98°C for 10 s, 50°C for 30 s, and 72°C for 30 s; 72°C for 5 min. The same volume of IX loading buffer (contained SYB green) was mixed with the PCR products and detected by electrophoresis on 2% agarose gels. The PCR products were mixed at an equidensity ratio. Then, the mixture PCR products were purified with Qiagen Gel Extraction Kit (Qiagen, Germany). Following the manufacturer’s recommendations, sequencing libraries were generated using TruSeq^® DNA PCR-Free Sample Preparation Kit (Illumina, United States) and index codes were added. Library quality was assessed by Qubit@2.0 Fluorometer (Thermo Scientific) and Agilent Bioanalyzer 2100 system (Agilent Technologies, CA, United States). Purified amplicons were sequenced on Illumina NovaSeq 6000 platform (2 × 250 paired ends) by Novogene Co., Ltd. (Beijing, China).

The assembly of paired-end reads, characterized by overlaps exceeding 10 bp, was accomplished using FLASH software (version 1.2.11) (18). Subsequently, reads of low quality, containing ambiguous characters and sequences shorter than 400 bp, were excluded. QIIME2 software (version 2023.5) (19) was employed for processing and assignment of these assembly readings. The denoise-paired method in DADA2 was applied to identify amplicon sequence variants (ASVs). Annotation of results utilized the SILVA 138 database (20), providing classifications at the kingdom, phylum, class, order, family and genus. QIIME2 (version 2023.5) facilitated the calculation of α and β diversity, and unweighted UniFrac distance metrics were used to generate principal coordinate analysis (PCoA). Differences in microbial composition between the 2 groups were assessed using Linear Discriminant Analysis Effect Size (LEfSe) (21), employing screening criteria of LDA score > 3 and p < 0.05. Additionally, we then predicted ileal microbial metabolic pathways using the PICRUSt2 (22) and used STAMP software (version 2.1.3) (23) to compare differences in function. Welch t-test (two-sided) was used for intergroup comparison, and the Welch’s inverted confidence interval (CI) method was used for CI calculation. p < 0.05 indicated a significant difference.

Total RNA extraction from the ileum was accomplished using Trizol (Invitrogen, Carlsbad, CA, United States), following the manufacturer’s instructions. RNA integrity was assessed using the Fragment Analyzer 5400 (Agilent Technologies, CA, United States). RNA integrity values range from 6.8 to 8.9. The library construction utilizing the obtained RNA was undertaken by Novogene Co., Ltd. (Beijing, China). All Illumina PE libraries were constructed and 2 × 150 bp RNA-seq was completed using the Illumina sequencing platform (NovaSeq 6000). The datasets presented in this study can be found in the National Center for Biotechnology Information (NCBI) under BioProject ID PRJNA1054312.

Standard quality control measures were implemented using Fastp software (version 0.23.1) (24) to filter out low-quality reads, ensuring the retention of clean reads for subsequent analyses. The obtained clean reads were aligned to the goose reference genome (BioProject ID PRJNA801885, data not released) using HISAT2 software (version 2.2.1) (25). The resulting SAM (sequencing alignment/mapping) files were subsequently converted to BAM (binary alignment/mapping) files and sorted using SAMtools (version 1.15.1) (26). The expression levels of each transcript were computed using featureCounts (version 2.0.3) (27), while gene expression was quantified using the transcripts per million (TPM) method.

DEseq2 (28) was employed to identify differentially expressed genes (DEGs) between groups. Genes exhibiting |Log₂(FC)| ≥ 1 and p < 0.01 were designated as DEGs. And the online tool KOBAS (version 3.0) (29) was utilized for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) functional analysis, accessible at http://kobas.cbi.pku.edu.cn/kobas3/?t=1. The functional gene analysis was conducted using G. gallus as a reference.

Statistical analysis was performed using SPSS 27.0 software. Spearman’s correlation coefficients were calculated to analyze the correlation. Differences were considered statistically significant at p < 0.05.

Following quality control and filtering procedures, a total of 957,019 effective reads were generated from 16 samples, averaging 59,814 reads per sample (Supplementary Table S1). A total of 745 ASVs were identified through QIIME2 analysis. Of these, 321 ASVs were found to be common to both rearing systems, while 221 ASVs were exclusive to the CR and 203 ASVs were unique to the FR (Figure 1). Subsequent taxonomic classification categorized these ASVs into 17 phyla, 34 classes, 88 orders, 167 families, 215 genera.

The complexity of the ileal microbiota was estimated on the basis of α-diversity indices (Observed_features, Shannon, Chao1, Faith’s PD, Evenness). The analysis revealed that the Faith’s PD α-diversity of CR was significantly lower than that of FR (Figure 2A). PCoA analyses demonstrated that the first principal coordinate (PCo1) explained 30.02% of the variations among samples and the second principal coordinate (PCo2) explained 18.34% of the variations (Figure 2B). The ANOSIM test further confirmed significant differences in ileal microbial communities between CR and FR (R = 0.22, p = 0.04). Comparison analysis of ileal microbiota at the phylum level and genus levels revealed broad similarities in dominant microbiota between the 2 rearing systems. At the phylum level, the dominant phyla in CR were Firmicutes (72.11%), Proteobacteria (10.39%) and Fusobacteriota (6.27%); in FR the dominant phyla were Firmicutes (54.64%), Fusobacteriota (19.76%) and Bacteroidetes (10.62%) (Figure 2C). At the genus level, in CR the top 3 genera were Romboutsia (45.03%), Clostridium_sensu_stricto_1 (12.26%) and Fusobacterium (6.78%); in FR, the top 3 genera were also Romboutsia (23.38%), Fusobacterium (21.13%) and Clostridium_sensu_stricto_1 (8.85%) (Figure 2D).

At the genus level, LEfSe analyses identified 2 and 14 differential microbiota from CR and FR, respectively (Figures 3A,B). In CR, the abundance of hgcI_clade and Faecalibacterium was significantly higher than that of FR. Meanwhile, Bacteroides, Proteiniphilum, Proteiniclasticum, Syner_01, Phascolarctobacterium, Sutterella, Thermovirga, Colidextribacter, Allorhizobium_Neorhizobium_Pararhizobium_Rhizobium, Paraclostridium, Petrimonas, Succinivibrio, Desulfovibrio, and Campylobacter had higher abundance in FR. In addition, PICRUSt2 predictive function analyses indicated that 4 metabolic pathways differed between the 2 groups, including Alzheimer disease, biosynthesis of amino acids, nicotinate and nicotinamide metabolism and pantothenate and CoA biosynthesis (Figure 3C).

A comprehensive set of 292,451,214 raw reads were generated across the 7 samples, and subsequent stringent filtering yielded an average of 40,961,820 clean reads for each sample. The quality metrics, including Q20, Q30, GC content, and mapping rates, exhibited favorable ranges of 97.15–97.81%, 92.82–94.08%, 44.72–52.03%, and 86.86–93.52%, respectively (Supplementary Table S2). These results attested to the robust sequencing quality essential for subsequent analyses. We identified a total of 1,198 DEGs (Figure 4A), of which 241 were up-regulated and 957 were down-regulated (Figure 4B). The clustering heatmap of TPM illustrated the expression profiles of DEGs in the ileum of CR and FR geese in a visually informative manner (Figure 4C).

Compared to FR, 957 genes were up-regulated in CR. These genes were enriched to 220 GO terms, including 129 biological processes (BP), 44 cellular components (CC), and 47 molecular functions (MF) (p < 0.05) (Figure 5A). The top 3 GO terms were cytosol, protein homodimerization activity, and lysosome. Based on KEGG pathway enrichment analysis, 28 pathways were identified significantly (p < 0.05) (Figure 5B). Metabolic pathways, lysosome, and glycosaminoglycan degradation were included. Further analysis revealed that 15 of the 28 pathways were related to metabolism, 5 to genetic information processing, 5 to cellular processes, 2 to environmental information processing, and 1 to organismal systems. Metabolism was the most up-regulated function, with 46.67% of the pathways associated with amino acid metabolism. These pathways included glycosaminoglycan degradation, biosynthesis of amino acids, glycosphingolipid biosynthesis—ganglio series, glutathione metabolism, taurine and hypotaurine metabolism, selenocompound metabolism, and glycine, serine and threonine metabolism.

Compared to CR, 241 genes were up-regulated in the FR and they included 119 BP, 39 CC, and 41 MF (p < 0.05) (Figure 5C). Integral component of membrane, cytosol, and GTPase activator activity were the GO terms with the highest significance. Through KEGG enrichment analysis, we identified 19 signaling pathways (p < 0.05), with the highest percentage of pathways associated with cellular processes (Figure 5D). Pathways associated with cellular processes included apoptosis, regulation of actin cytoskeleton, necroptosis, cellular senescence, and gap junction were up-regulated in FR.

Combining these findings, we generated heat maps based on the results of Spearman correlation analysis. Illustrated in Figure 6A, we performed Spearman correlation analysis between up-regulated DEGs involved in amino acid metabolism and differential microbiota in CR. The abundance of hgcI_clade was significant positively correlated with the expression of MARS, NAGLU, ACY1, IDH2, ST6GALNAC6, HYAL1, SGSH, SCLY, GGT7, and GRHPR, while demonstrating a significant negative correlation with the expression of GLB1 (p < 0.05). Simultaneously, we analyzed the correlation between up-regulated DEGs involved in cellular processes and differential microbiota in FR, and the results showed that FGF10 was significantly and positively correlated with Phascolarctobacterium (p < 0.05), and, PIK3R1 was significantly and positively correlated with Sutterella (p < 0.01) (Figure 6B).