Proteome and microbiota analysis reveals alterations of liver-gut axis under different stocking density of Peking ducks.

This study aimed to determine the impact of stocking density on the liver proteome and cecal microbiota of Peking ducks. A total of 1,200 21-day-old ducks were randomly assigned to 5 stocking density groups of 5, 6, 7, 8 and 9 ducks/m2, with 6 replicates for each group. At 40 days of age, duck serum and pectorals were collected for biochemical tests; liver and cecal contents of ducks were gathered for proteome and microbiota analysis, respectively. Serum MDA increased while pectorals T-AOC reduced linearly with enhancing stocking density. Duck lipid metabolism was altered under different stocking density as well. Serum LDL-C increased linearly with increasing stocking density. Proteome analysis revealed fatty acid biosynthesis proteins such as acyl-CoA synthetase family member 2 and fatty acid oxidation related proteins including acyl-CoA dehydrogenase long chain and acyl-coenzyme A oxidase were enriched in high stocking density group. Additionally, high stocking density increased oxidative response associated proteins such as DDRGK domain containing 1. Furthermore, increasing stocking density diminished proteins of anti-oxidant capacity including regucalcin and catalase. 16S rDNA analysis revealed that higher stocking density was accompanied with decreased microbial diversity, as well as depletion of anti-inflammatory bacterial taxa, including Bacteroidales, Butyricimonas and Alistipe. Besides, reduced bile acid metabolism-associated bacteria such as Ruminococcaceae, Clostridiales and Desulfovibrionaceae were found in the high-density group. Both proteome and 16S rDNA results showed inflammation and chronic liver disease trend in the high-density group, which suggests the involvement of the liver-gut axis in oxidative stress.

Introduction

In China, Peking Roast Duck is one of the most popular dishes, and the demand of Peking duck is increasing in the recent years. The pursuit of high economic interests by producers has accelerated the evolution of the Beijing duck feeding system from free-range to intensive [1]. Increasing stocking density can yield higher profits per kilogram of chicken. However, reduced space limits the ability to rest [2] and has a negative influence on performance [1,3,4], meat yield [3,4], immune status [5] and gut morphology [6]. In addition, high stocking density is associated with chronic oxidative stress [3]. However, these studies focused on broilers rather than ducks. Studies found the effects of high stocking density on ducks were similar with broilers. Increasing stocking density significantly decreased final body weight (BW) and weight gain of starter and growing ducks [1,7]. High stocking density also impaired immune function and anti-oxidative capacity of Peking ducks, high stocking density showed a lower spleen weight and lower total anti-oxidative capacity [8]. Although the transition from conventional free range outdoors to confinement in birdhouses is the development direction of duck production [7], the study on the stocking density of duck production is limited. The liver plays an important role in energy metabolism and it is the major site of triglyceride (TG) metabolism which is involved in TG digestion, absorption, synthesis, decomposition and transport. It has been demonstrated that high stocking density has been shown to have a deleterious effect on liver function [9] as indicated by increased activities with aspartate aminotransferase and alanine aminotransferase [10]. The liver total anti-oxidative capacity was decreased under high stocking density [8]. The gut microbiome is also highly connected to animal energy metabolism and health [11], and it has been frequently suggested that gut microbiota plays a critical role in chronic liver disorders through liver-gut axis [12]. High stocking density is associated with adverse effects on the chicken intestinal commensal bacteria [4]. However, it is unclear that whether stocking density could influence liver metabolism or gut microflora.

The emergence of novel proteomic techniques in recent years has greatly aided in the understanding of biological mechanisms. Tandem mass tag (TMT) [13] and isobaric tags for relative and absolute quantitation (iTRAQ) [14] methods have been widely used for analyzing the hepatic proteome. The TMT method can also be used to characterize liver proteome-wide changes in response to oxidative stress [15]. Further, due to progress in high-throughput next-generation sequencing, 16S rDNA analysis can be used to infer the structure and function of gut microbiota.

To investigate whether stocking density could affect liver function and gut microbiome, TMT-labeled quantitative proteomics combined with 16S rDNA analysis was used to identify changes in the protein profiles and microbiota of Peking ducks under low and high raising density.

Materials and methods

Ethics approval

All of the experiments were approved by the Institutional Animal Care and Use Committee of the China Agricultural University (Beijing, China).

Facilities and experimental animals

A total of 1,200 mixed-sex 21-d-old white Peking ducks were randomly assigned to 5 stocking density treatments, each having 6 replicates. Each replicate corresponds one pen with 40 ducks (20 males, and 20 females). The raising area in different treatments was 8.20 m2 (2.50 × 3.28 m), 6.88 m2 (2.50 × 2.75 m), 5.93 m2 (2.50 × 2.37 m), 5.20 m2 (2.50 × 2.08 m) and 4.65 m2 (2.5 × 1.86 m). The corresponding stocking density was 5 (low stocking density represented as L group), 6, 7, 8, and 9 (high stocking density represented as H group) ducks/m2, respectively. All ducks were raised in a plastic wire-floor pen and provided water and feed ad libitum from 21 to 40 d of age. In the house, lighting, temperature and ventilation programs followed commercial practices. During the experimental period, all ducks were raised with diets based on the NRC (1994) feeding standard.

Growth performance and carcass traits

Initial body weight, final BW and feed intake (FI) were recorded. The weight gain and feed/gain ratio were calculated.

Sample collection

On day 40, a total of 6 male ducks per treatment were collected and sacrificed by stunning after blood collection. The liver samples of ducks from L and H density treatments were collected, flushed with cold PBS, frozen using liquid nitrogen, and stored at -80°C for proteomic analysis. Cecum contents of ducks from L and H density treatments were also collected, transferred into Eppendorf tubes, and immediately frozen in liquid nitrogen and stored at -80°C for microbiota analyses.

Abdominal fat and left pectorals of each duck were removed manually from the carcass and weighed. For each duck, a piece of pectorals fixed position was separated and put on the ice bag immediately, then stored at -80°C for biochemical analysis.

Serum and pectorals biochemical parameters

Serum levels of malondialdehyde (MDA, cat#A003-1), TG (cat#A110-2), total cholesterol (TC, cat#A111-2), high-density lipoprotein cholesterol (HDL-C, cat#A112-2), low-density lipoprotein cholesterol (LDL-C, cat#A113-2), total antioxidant capacity (T-AOC, cat#A015-1) and the activities of lactate dehydrogenase (LDH, cat#A020-1), creatine kinase (CK, cat#A032), lipoprotein lipase (LPL, cat#A067) along with MDA, T-AOC, LPL and protein concentration (cat#A045-2) of pectorals were determined by using commercial analytical kits according to the manufacturer’s recommendations (Jian Cheng Bioengineering Institute, Nanjing, China).

Statistical analysis

One-way analysis of variance (ANOVA) models of SPSS (v.20.0, SPSS Institute, Chicago, IL) was fitted to determine the relationships between the stocking density groups and serum parameters, pectoral redox and lipid metabolism indices. Means were compared using Duncan’s multiple comparison procedure of SPSS software when density treatment was significant (P < 0.05), and curve estimation was used to assess the linear and quadratic effects of increasing stocking density on final body weight, body weight gain, feed gain ratio, and pectorals percentage. P < 0.05 indicated statistical significance, and results were considered as a significant trend at P < 0.1.

Proteome analysis

Protein extraction

Liver samples (~100mg each) were ground to a fine powder in liquid nitrogen and then transferred into a centrifuge tube. Four times volumes of lysis buffer containing 8 M urea (Sigma) and 1% Protease Inhibitor Cocktail (Calbiochem) was added, followed by sonication on ice for three times using a high-intensity ultrasonic processor (Scientz), and centrifuged for 10 min at 4°C and 12,000g. Finally, supernatants containing soluble proteins were collected, and the protein concentration was quantified with BCA kit (Beyotime Biotechnology) according to the manufacturer’s instructions.

Trypsin digestion

For trypsin digestion, the protein solution was reduced with 5 mM dithiothreitol (Sigma) for 30 min at 56°C and alkylated with 11 mM iodoacetamide (Sigma) for 15 min at room temperature under dark conditions. The protein sample was diluted with 100 mM triethylammonium bicarbonate (TEAB, Sigma) to make urea concentration less than 2 M. Finally, the sample was digested with trypsin (Promega) overnight with a 1:50 trypsin-to-protein ratio and then further digested for 4 h with 1:100 w: w trypsin-to-protein ratio (37°C).

TMT labeling

After trypsin digestion, the peptides were desalted in a Strata X C18 SPE column (Phenomenex) and vacuum dried. The peptide was dissolved in 1 M TEAB and processed according to the manufacturer’s protocol for TMT kit. Briefly, one unit of TMT reagent (defined as the amount of reagent requirement for labeling 100 μg of protein) was thawed and dissolved in acetonitrile (Fisher Chemical). The peptide mixtures were then incubated for 2 h at room temperature, pooled, desalted and dried by vacuum centrifugation.

HPLC fractionation

Peptides were fractionated by high pH reverse-phase HPLC with an Agilent 300 Extend C18 column (5 μm particles, 4.6 mm ID, 250 mm length). Briefly, peptides were firstly separated by 8% to 32% acetonitrile (pH 9.0) over 60 minutes into 60 fractions. Subsequently, the peptides were combined into 18 fractions and dried by vacuum centrifugation.

Mass spectrometry and TMT data analysis

Peptides were dissolved in 0.1% formic acid (Fluka), and separated by EASY-nLC 1000 Ultra Performance Liquid Chromatography (UPLC) system. Solvent A is an aqueous solution containing 0.1% formic acid and 2% acetonitrile; solvent B is an aqueous solution containing 0.1% formic acid and 90% acetonitrile. Running conditions included a liquid-phase gradient using a linearly increasing gradient of 5% to 24% solvent B for 38 min, 24% to 35% of solvent B for 14 min, then climbing to 80% in 4 min, and then maintaining at 80% for the last 4 min, all at a constant flow rate of 800 nl/min.

The peptides results were subjected to nano electrospray ionization (NSI) source followed by tandem mass spectrometry (MS/MS) in Q Exactive PlusTM (Thermo Fisher Scientific) coupled to the UPLC. The electrospray voltage applied was 2.0 kV. For MS scans, the m/z scan range was 350 to 1800. Intact peptides were detected in an orbitrap at a resolution of 70,000. Peptides were then selected for MS/MS in the orbitrap at a resolution of 17,500. Fixed first mass was set to 100 m/z. A data-dependent procedure that alternated between one MS scan was used, followed by 20 MS/MS scans. To improve the effective utilization of the mass spectrometer, the automatic gain control (AGC) was set to 5E4, the signal threshold was set to 10000 ions/s, the maximum injection time was set to 200 ms, and the dynamic exclusion time of the tandem mass scan was set to 30 seconds to avoid repeating a scan of the parent ion.

Database search

The resulting MS/MS data were processed using Maxquant with an integrated Andromeda search engine (v.1.5.2.8). Tandem mass spectra were searched against the uniprot Anas platyrhynchos database concatenated with reverse decoy database. Trypsin/P was specified as cleavage enzyme allowing up to 2 missing cleavages. The mass tolerance for precursor ions was set as to 10 ppm in First search and 5 ppm in Main search, and the mass tolerance was set as 0.02 Da for fragment ions. TMT 10-plex was selected for protein quantification. False discovery rate (FDR) was adjusted to < 1% for protein identification.

TMT quantification

For TMT quantification, the ratios of the TMT reporter ion intensities in MS/MS spectra (m/z 126–131) from raw data sets were used to calculate fold changes between samples. For each sample, the quantification was mean-normalized at the peptide level to center the distribution of quantitative values. Protein quantitation was then calculated as the median ratio of corresponding unique or razor peptides for a given protein. The mass spectrometry proteomics data have been deposited to the ProteomeXchange Consortium via the PRIDE [16] partner repository with the dataset identifier PXD010952. Two-sample, two-sided T-tests were used to compare the expression of proteins. A significance was set at P < 0.05 and the significant trend was set at P < 0.1 for statistical analysis.

Gene Ontology (GO) annotation

GO annotation proteome was derived from the UniProt-GOA database (www.http://www.ebi.ac.uk/GOA/). Firstly, Converting identified protein ID to UniProt ID and then mapping to GO IDs by protein ID. If some identified proteins were not annotated by UniProt-GOA database, the InterProScan soft would be used to annotated protein’s GO functional based on protein sequence alignment method. Then proteins were classified by Gene Ontology annotation based on three categories: biological process, cellular component and molecular function.

16S rDNA sequencing

DNA was extracted from 180–220 mg of the cecal samples using a QIAampTM Fast DNA Stool Mini Kit (Qiagene, No. 51604) according to the manufacturer’s instructions. Total DNA was quantified using a Thermo NanoDrop 2000 UV microscope spectrophotometer and 1% agarose gel electrophoresis. 16S rDNA high-throughput sequencing was performed by Realbio Genomics Institute (Shanghai, China) using the Illumina Hiseq PE250 platform. The V3-V4 region of the 16S rDNA gene was amplified using the universal primers, 341F (CCTACGGGRSGCAGCAG) and 806R (GGACTACVVGGGTATCTAATC). The raw pair-end reads were merged and quality-filtered to remove tags with lengths < 220 nt, an average quality score of <20, and tags containing >3 ambiguous bases using PANDAseq (v2.9) [17]. Singletons and chimeras were removed, and the resulting quality-filtered sequences were clustered into 97% operational taxonomic units (OTUs) using USEARCH (v7.0.1090) in QIIME software. The Ribosomal Database Project (RDP) algorithm trained on the Greengenes database was used to classify each OTU (http://greengenes.lbl.gov). The open source software package QIIME (http://qiime.org) was used to measure alpha diversity (including the chao1, observed species and PD whole tree indices). The metagenomic reads have been submitted to the NCBI-SRA database under accession number SRP159441.

16S rDNA analysis

For 16S rDNA, the Wilcoxon rank sum test to evaluate changes in alpha diversity between H and L groups. Venn diagrams were constructed in R v3.1.0 using the Venn Diagram package. Principal component analysis (PCA) was used to assess the relationships between samples based on the composition of the microbiota [18]. The bacterial relative abundance profiles were further compared between two groups by Linear Discriminant Analysis (LDA) Effective Size (LEfSe) in the literature [19]. Taxa with the logarithmic LDA score >2 were identified as differentially abundant.

Results

Duck performance and carcass traits

In our study, both the final body weight and the body weight gain over the study period decreased linearly with the increase of stocking density (Fig 1A and 1B). The feed/gain ratio increased linearly with increasing density (Fig 1C). Pectorals percentage was affected by stocking density as well, decreasing linearly with increasing stocking density (Fig 1D).

Fig 1

Duck performance and carcass traits.

Data presented as means ± SE. n = 6.

Biochemical indices of serum and pectorals

Increasing stocking density suppressed anti-oxidant capacity. In redox indices, raising density was significantly associated with serum MDA and LDH, as well as pectoral T-AOC (Table 1). Serum MDA increased linearly, while pectoral T-AOC decreased linearly with increasing stocking density (Table 1).

Table 1

Redox related biochemical indices in serum and pectorals.

	Pectorals		Serum
density/m2	T-AOC	MDA	T-AOC	MDA	LDH
5	0.71±0.15AB	2.12±0.60	15.79±2.30	7.39±0.54B	9186.19±888.64A
6	0.50±0.14AB	1.42±0.59	13.03±5.35	7.47±0.51B	9224.32±195.84A
7	0.48±0.28AB	1.74±0.77	10.77±5.51	7.10±0.58B	10192.19±732.17A
8	0.46±0.16AB	1.37±0.44	14.55±1.23	8.68±0.05AB	7864.87±814.26B
9	0.22±0.08B	1.87±0.46	11.51±3.95	9.57±1.56A	9225.23±273.26A
P-value
ANOVA	0.015	0.176	0.278	0.008	0.001
Linear	0.023	0.490	0.253	0.038	0.238
Quadratic	0.069	0.148	0.269	0.068	0.470

Note: Total antioxidant capacity (T-AOC), malondialdehyde (MDA), lactate dehydrogenase (LDH).

Stocking density was also significantly altered lipid metabolism indices including serum LPL, TG, TC, HDL-C and LDL-C, as well as pectoral LPL (Table 2). Additionally, serum LDL-C increased linearly, and HDL-C had a linear increasing trend with increasing density. Serum TG and pectoral LPL had quadratic trends associated with stocking density (Table 2).

Table 2

Lipid metabolism related biochemical indices.

	Pectorals	Serum
Group	LPL	LPL	TG	TC	HDL-C	LDL-C
5	0.67±0.10BC	2.76±0.62A	1.01±0.22C	4.14±0.32B	4.01±0.30	3.24±0.60B
6	0.99±0.23AB	1.15±0.17B	1.54±0.07B	4.17±0.30B	4.76±0.68	3.59±0.59B
7	1.13±0.27A	2.67±0.39A	1.31±0.12BC	4.13±0.15B	4.41±0.83	4.83±0.73A
8	0.84±0.30ABC	2.52±0.79AB	2.09±0.20A	5.03±0.36A	5.10±0.84	4.88±0.85A
9	0.55±0.09C	2.80±0.36A	1.58±0.16B	4.74±0.10A	5.30±0.11	5.02±0.69A
P-value
ANOVA	0.003	<0.001	<0.001	<0.001	0.111	0.001
Linear	0.990	0.311	0.089	0.118	0.095	<0.001
Quadratic	0.098	0.506	0.050	0.302	0.250	<0.001

Note: lipoprotein lipase (LPL) triglyceride (TG), total cholesterol (TC), high-density lipoprotein cholesterol (HDL-C) and low-density lipoprotein cholesterol (LDL-C).

Proteins identification and comparison

In GO analysis, the most striking different pathways between H and L group were small molecule metabolic process, organic acid metabolic process, oxoacid metabolic process and carboxylic acid metabolic process, high stocking density elevated all of these pathways. Pathways including oxidation-reduction process, cellular amino acid metabolic process and fatty acid metabolic process were significantly increased in H group as well (Fig 2). Many proteins involved in redox, lipid metabolism, protein turnover, DNA repair, and immunity were associated with stocking density (Table 3).

Fig 2

GO pathways enriched by H group.

Log-transformed P-value indicates the degree of difference between the H group and L group.

Table 3

Differentially expressed proteins identified under different stocking density.

Accession IDa	Annotation	Gene	H:Lb	P-value
	Redox related proteins
U3I9D6	DDRGK domain containing 1	DDRGK	1.218	0.010
MT	Metallothionein	N/A	1.332	0.014
U3ID81	Regucalcin	RGN	0.706	0.025
R0M3U9	Catalase	Anapl_09675	0.882	0.051
	Lipid metabolism related proteins
U3IAY7	Acyl-CoA dehydrogenase, long chain	ACADL	1.154	0.040
U3I7T9	Acyl-CoA synthetase family member 2	ACSF2	1.103	0.032
U3J928	Acyl-coenzyme A oxidase	ACOX1	1.098	0.007
U3I7I1	2-hydroxyacyl-CoA lyase 1	HACL1	1.143	0.095
	Protein turnover related proteins
U3I9U9	Pyridoxal kinase	PDXK	1.324	0.011
U3IHG9	Ribosomal protein S20	RPS20	1.103	0.045
U3I350	Aspartyl-tRNA synthetase	DARS	1.119	0.096
U3J7Z7	Asparaginyl-tRNA synthetase	NARS	0.904	0.068
	DNA repair related proteins
U3IUN7	RNA binding motif protein, X-linked	RBMX	1.136	0.068
U3IUB7	Lamin B2	LMNB2	1.119	0.076
	Immunity related proteins
R0K747	Macrophage mannose receptor 1	Anapl_08603	1.170	0.004
Q6JWQ2	MHC class I antigen alpha chain	Anpl-U	0.785	0.021

Note

a UniProt accession ID (http://www.uniprot.org/).

b ratio of protein abundance in the H group compared to the protein abundance in the L group

Antioxidant related proteins such as regucalcin and catalase were downregulated in H group. Oxidative response associated proteins like DDRGK domain containing 1 and metallothionein were enhanced in the H group.

Fatty acid synthesis protein acyl-CoA synthetase family member 2 (ACSF2) was enriched in high stocking density. Fatty acid oxidation connected proteins including acyl-CoA dehydrogenase long chain (ACADL), acyl-coenzyme A oxidase (ACOX1) and 2-hydroxyacyl-CoA lyase 1 (HACL1) were enriched in the H group.

Protein turnover related proteins including Ribosomal protein S20, Pyridoxal kinase and aspartyl-tRNA synthetase (DARS2) were elevated while Asparaginyl-tRNA synthetase was reduced in the H group.

In DNA repair associated proteins, both RNA-binding motif protein X-linked (RBMX) and Lamin B2 showed enhancing trends in H group. The immunity associated macrophage mannose receptor (MR) 1 was elevated, while major histocompatibility complex (MHC) class I antigen α chain was decreased in the H group.

16S rDNA analysis of cecal microbiota

Microbiota PCA analysis showed a clear separation of samples from the H group and samples from the L group (Fig 3A). The H group and L group had 290 bacterial taxa in common, while the H group and L group had 35 and 13 unique taxa, respectively (Fig 3B).

Fig 3

PCA and Venn chart of microbiota.

a, PCA diagram; b, Venn chart.

LEfSe analysis, revealed a higher relative abundance of Firmicutes and Phascolarctobacterium in the H group, while Bacteroidia, Bacteroidales, Thermoplasmata, Methanomassiliicoccales, Methanomassiliicoccus and Methanomassiliicoccaceae were more abundant in the L group (Fig 4A and 4B). Additionally, the H group had elevated Lachnospiraceae and Ruminococcaceae lower abundance of Butyricimonas, Desulfovibrionaceae, Alistipes and Clostridiales (Fig 5). Overall, the ratio of Firmicutes to Bacteroidetes in the H group was higher than in the L group (P = 0.058).

Fig 4

Bacterial taxa in H or L groups by LEfSe analysis.

Phylogenetic relationships among significant bacterial biomarkers are indicated in the cladogram (top). Log-transformed linear discriminant analysis (LDA) scores of the significant biomarkers are stated in the bar chart (bottom).

Fig 5

Relative abundance of anti-inflammatory bacterial taxa in H or L groups.

Discussion

Previous studies have indicated that high raising density has negative influences on physiology and overall meat production [3,4,7]. This study showed a linear decrease of BW with increasing stocking density. Similarly, breast muscle was depressed by increasing stocking density, which was consistent with the former research that increasing stocking density decreased breast fillet weight and its relative yield [20].

Increased density is associated with impaired antioxidant capacity, including decreasing total glutathione concentration and the glutathione (GSH): oxidized glutathione (GSSG) ratio [3]. Liver antioxidants including superoxide dismutase (SOD), catalase (CAT), GSH have previously been found to be reduced under high raising density [21]. In this study, enhancing stocking density increased serum MDA and decreased pectoral T-AOC. In the proteomic analysis of this study, regucalcin and catalase expression were reduced in high density group. Catalases are involved in protecting cells from the damaging effects of ROS [22], while regucalcin is a calcium-binding protein with multiple roles that include calcium homeostasis, anti-oxidative, anti-apoptosis, and anti-proliferation [23]. Decreases in these proteins indicate that high stocking density may reduce duck anti-oxidant capacity.

High stocking density increases blood oxidative stress, including acute phase proteins, heat shock protein 70, and circulating corticosterone [24]. In the current study, high stocking density had higher expressions of DDRGK domain containing 1, metallothionein, DARS2 and RBMX. DDRGK domain containing 1 and metallothionein were enriched in high density group. DDRGK1 is an endoplasmic reticulum membrane protein and plays an important role in maintaining the homeostasis of endoplasmic reticulum [25]. Metallothionein is an essential protein for the protection of cells against reactive oxygen species (ROS) [26]. Therefore, metallothionein can neutralize ROS and protect host from oxidative stress. A previous study found that loss of mitochondrial DARS2 leads to the activation of various stress responses in a tissue-specific manner independently of respiratory chain deficiency [27]. The RBMX is a nuclear protein that is involved in alternative splicing of RNA [28]. It can confer resistance to DNA damage [29]. Enhancement of these proteins suggests oxidative stress of Peking ducks under high stocking density.

Lipid storage is essential for protection against ROS toxicity [30]. The lipogenesis occurs in liver most exclusively, especially to waterfowls [31]. High stocking density was previously associated with higher hepatic TG storage [32]. This study showed serum LDL-C increased with density incrementing. Moreover, high density group elevated ACSF2 expression. Acetyl-CoA synthetase catalyzes the formation of acetyl-CoA from acetate, CoA and ATP and participates in various metabolic pathways, including fatty acid and cholesterol synthesis and the tricarboxylic acid cycle [33]. ACSF2 belongs to the acyl-CoA synthetase family, activating fatty acids by forming a thioester bond with CoA [34]. Therefore, the increasing lipid biosynthesis may be the self-protection mechanism of Peking ducks under oxidative stress.

ROS are considered to be involved in the progression of non-alcoholic fatty liver disease (NAFLD) [35]. Interestingly, fatty liver production in waterfowls is relatively similar to human NAFLD [36]. Therefore, waterfowl can be a good model for liver steatosis research. In current study, ACADL, ACOX1 and HACL1 were enriched in high stocking density group. ACADL is a crucial enzyme participating in fatty acid oxidation [37]. Similarly, ACOX1 are involved in β-oxidation in the liver [38]. HACL1 has two critical roles in α-oxidation, the degradation of phytanic acid and shortening of 2-hydroxy long-chain fatty acids so that they can be degraded further by β-oxidation [39]. It has been reported that inhibition of β-oxidation decreases NADPH levels and increases ROS levels [40]. Therefore, the elevation of these proteins may protect ducks from oxidative stress. It is increasingly recognized that the composition of the gut microbiota plays a critical role in influencing predisposition to chronic liver disorders such as NAFLD [12]. Samples from high stocking density had high levels of Lachnospiraceae. Interestingly, high Lachnospiraceae abundance was observed in patients with nonalcoholic steatohepatitis [41]. Ruminococcaceae and Alistipes were depleted in the high density group. Cirrhotic patients were previously shown to have lower Ruminococcaceae (7α-dehydroxylating bacteria) abundance compared to healthy patients [42]. Also, Alistipes was significantly more abundant in the gut microbiota of healthy subjects compared to NAFLD patients [43].

Oxidative stress can regulate immunity. A decreased broiler bursa weight was previously reported to be associated with higher stocking densities [3]. DDRGK domain-containing protein 1 and Macrophage mannose receptor (MR) expression was higher in high density compared to the low density group. DDRGK1 is also an essential regulatory protein of NF-κB [44]. A recent study found MR can protect against ROS burst [45]. The increase of these proteins reflects the immune response status of Peking ducks under high stocking density. MHC class I antigen α chain was decreased under high stocking density. MHC class I plays a crucial role in immunity by capturing peptides for presentation to T cells and natural killer (NK) cells [46]. Dysregulation of MHC class I was correlated with unfolded protein response (UPR) and endoplasmic reticulum (ER) stress [47]. Upregulation of these proteins in the high density group suggests the immune adaption to high stocking. Moreover, high stocking density has been previously associated with adverse effects on the chicken intestinal commensal bacteria [4]. In the current study, Phascolarctobacterium was enriched, while Bacteroidales, Butyricimonas and Alistipe are depleted in high density group. A previous study found that Phascolarctobacterium was significantly correlated with systemic inflammatory cytokines [48]. Besides, depletion of Bacteroidales has previously been associated with disease status [49]. Additionally, reduction in Butyricimonas is associated with increased proinflammatory gene expression [50]. Studies confirm that patients with IBD and Clostridium diffcile infection have a lower abundance of Alistipe than their healthy counterparts [51]. The current study showed higher Firmicutes to Bacteroidetes ratio in high density group. The microbiota of irritable bowel syndrome (IBS) patients, compared with controls, had a 2-fold increased ratio of the Firmicutes to Bacteroidetes [52]. Bile acids are ensential metabolites of the microbiome and can modulate the composition of the gut microbiota directly or indirectly through the activation of the innate immune system [53]. Bile acid plays a crucial role in the control of inflammation and NAFLD [54]. Samples from the high stocking density group had higher levels of Ruminococcaceae and lower abundance of Desulfovibrionaceae, and Clostridiales compared to low density group. Cholesterol 7α-hydroxylase (CYP7A1) is the enzyme responsible for catalyzing the first and rate-limiting step in the classical bile acid synthetic pathway [55]. Furthermore, an inverse relationship between Cyp7a1 expression and Ruminococcaceae abundance has been previously demonstrated [21]. A study found enriched Desulfovibrionaceae was accompanied by increased hepatic taurine-conjugated cholic acid and β-muricholic acid, which were the main constituent of bile acid pool [56]. Clostridiale was positively correlated with muricholic acid as well [21]. A loss of bacteria belonging to the Clostridiales order was connected with a disturbance in the bile-microbial axis [57]. Therefore, alternations of bacteria showed a decrease trend in bile acid synthesis.

In conclusion, high stocking density caused alternations of liver proteome and gut microbiome, which may cause oxidative stress of Peking ducks (Fig 6).

Fig 6

Graphic summary of liver proteome and gut microbiota alternations under high stocking density.