Transcriptome analysis reveals normalization effect of nicotinamide and butyrate sodium on breast muscles of broilers under high stocking density

Background In recent years, increased attention has been focused on breast muscle yield and meat quality in poultry production. Supplementation with nicotinamide and butyrate sodium can improve the meat quality of broilers. However, the potential molecular mechanism is not clear yet. This study was designed to investigate the effects of supplementation with a combination of nicotinamide and butyrate sodium on breast muscle transcriptome of broilers under high stocking density. Methods A total of 300 21-d-old Cobb broilers were randomly allocated into 3 groups based on stocking density: low stocking density control group (L; 14 birds/m 2 ), high stocking density control group (H; 18 birds/m 2 ), and high stocking density group provided with a combination of 50 mg/kg nicotinamide and 500 mg/kg butyrate sodium (COMB; 18 birds/m 2 ), raised to 42 days of age. Results The H group signicantly increased cooking losses, pH decline and activity of lactate dehydrogenase in breast muscle while COMB showed a signicant decrease in these indices ( P < 0.05). The transcriptome results showed that key genes involved in glycolysis, proteolysis and immune stress were up-regulated whereas those relating to muscle development, cell adhesion, cell matrix and collagen were down-regulated in the H group. In contrast, genes related to muscle development, hyaluronic acid, mitochondrial function, and redox pathways were up-regulated while those associated with inammatory response, acid metabolism, lipid metabolism, and glycolysis pathway were down-regulated in the COMB group. Conclusions The combination of nicotinamide and butyrate sodium may improve muscle quality by enhancing mitochondrial function and antioxidant capacity, inhibiting inammatory response and glycolysis, and promoting muscle development and hyaluronic acid synthesis. Fluorescence quantitative PCR was performed according to the of TB GreenTM TM Plus) The reaction apparatus was a 7500 uorescence detection system (Applied Biosystems), and the PCR reaction conditions were as follows: after pre-denaturation at 95 °C for 30 s, 40 cycles of denaturation at 95 °C for 5 s and 60 °C for 34 s. After the completion of PCR amplication, the dissolution curve was observed, and agarose gel electrophoresis was performed to identify whether the amplied gene fragment conformed to the design length, and the specicity of the amplication product was veried. The target gene and the internal reference gene beta-actin primer sequence are shown in Table 2. The results of gene expression were analyzed and compared using 2 -ΔΔCT .

well (P < 0.05) (Figure 1). The 45-min pH value in the H group was higher than that in the other 2 groups (P < 0.05) while there was no signi cant difference in 24-h pH values among the groups. Thus, the pH decline during 45 min to 24 h in the H group was signi cantly higher than that in the other 2 groups, indicating that the H group had rapid pH drop rate, which was attenuated in the COMB group under high stocking density ( Figure 2).

Anti-oxidant capacity
The activity of LDH in the H group was higher (P < 0.05) than that in the L group. The COMB group had signi cantly decreased (P < 0.05) activity of LDH when compared with the H group. However, stocking density had no signi cant effect on the activities of CK, T-AOC, MDH, antisuperoxide anion and the content of hydroxyproline (Table 3). The gene sets were produced by DEGS. From Venn analysis of genes sets, we found that there were 1310 genes shared in common between the COMB group and the L group. Nevertheless, there were only 6 genes owed by both the COMB group and the H group. Similarly, from the iPath map of metabolic pathways, there were a total of 830 pathways annotated in common. In contrast, there was only 1 pathway owed by both the COMB group and the H group ( Figure 5).

Up-regulated genes in the H group
Compared with those in the L group, a total of 1894 genes were up-regulated in the H group (Figure 4), which were mainly involved in muscle contraction, cell localization, ion transport, lipid metabolism, glycolysis, proteolysis, and immune stress ( Figure 6).
Muscle contraction-related pathways were enriched in the H group. They involved vital genes including MYLK2, NOS1, TMOD4, and Six1 ( Table   4). The H group was enriched for cell-localization-related genes such as KEAP1, CDKN1A, ERBB4, and TMOD4 (Table 4). Additionally, highdensity up-regulated ion and amino acid transport-related genes included KCNJ12, KCNA7, SLC38A3 and SLC38A4, which are involved in ion transmembrane transport and transporter activity (Table 5). High-density enriched glycolysis-related pathways included fructose metabolism, fructose-2,6-diphosphate 2-phosphatase activity, and fructose 2,6-diphosphate metabolism ( Table 6). The lipid metabolism-related genes such as MID1IP1, ACACB and Lpin1 were up-regulated in H group, which are involved in lipid synthesis and lipid oxidation (Table 6).
Stress response pathways including non-biologically stimulated cellular responses, extracellular stimuli response and nutritional level response were also enriched in the H group. Furthermore, high-density up-regulated proteolysis-related genes include TINAG, USP24, OTUD1, KEAP1, KLHL34, and SMCR8. Also, high-density enriched immune pathways include the regulation of host defence responses to viruses and prostaglandin receptor-like binding (Table 7).
In Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis, genes involved in calcium signalling pathway (RYR), in ammatory mediator regulation of RTP channels (PLA2) and chemokine signalling pathway (SOS) (Fig. S1-S3) were enriched in the H group.

Down-regulated genes in the H group
Compared with those in the L group, a total of 1858 genes were down-regulated in the H group (Figure 4), which were involved in cell adhesion, cell matrix, and cell migration, etc ( Figure 7).
Up-regulated genes in the COMB group Compared with those in the H group, up-regulated genes in the COMB group were involved in muscle development, hyaluronic acid synthesis, mitochondrial function, and redox pathway (Figure 8).
The muscle development-related pathways enriched in the COMB group included positive regulation of muscle tissue development and muscle  cell decision processes, which involved key genes such as MYF6, LMCD1 and TRPC3. Besides, the COMB group was enriched for   mitochondria-associated pathways such as electron transport chains, mitochondrial respiratory chain complex I and mitochondrial protein   complex pathways, which involved genes including TOMM6, NDUFV1, NDUFS5, NDUFB2, NDUFA2, LMCD1, ZNF593 and COASY (Table 11).
The hyaluronic acid-related genes up-regulated in the COMB group included HYAL1 and HYAL3. Besides, the redox-related genes including LDHD, CPOX, SUOX, NDUFV1, GRHPR, DOHH and NDUFA2 were up-regulated in the COMB group, which were involved in the pathways such as redox process, NAD binding, NADPH binding and NADH dehydrogenase complex (Table 12). In KEGG enrichment analysis, up-regulated genes in the COMB group were involved in oxidative phosphorylation (NDUFS5, NDUFV1, NDUFA2, NDUFA13, NDUFB2, NDUFB7 and NDUFC2) (Fig.   S7).

Down-regulated genes in the COMB group
Compared with those in the H group, down-regulated genes in the COMB group were involved in the in ammatory response, acid metabolism, fatty acid metabolism, and glycolysis-related pathways ( Figure 9).
The in ammatory response-related genes down-regulated in the COMB group included CCR5 and ALOX5 while the immune response-related genes included C1S, BLK, CCR5 and MARCH1 (Table 13). The acid metabolism-related pathways include organic acid synthesis process, oxoacid metabolism process and carboxylic acid synthesis process, which involved genes such as PSAT1, SCD, MAT1A, ALOX5, ST3GAL1 and ALDOB. The genes involved in fatty acid metabolism pathways include SCD and ALOX5. In addition, down-regulated genes in the COMB group were involved in glycolytic and carbohydrate metabolism, which included GALNT16, ST3GAL1, ALDOB and MAT1A (Table 13).
Transcriptome differential gene veri cation The transcriptome differential genes were veri ed by real-time PCR, and the gene expression pattern was consistent with the transcriptome results ( Figure 10).

Discussion
In the current study, the H group showed signi cantly increased cooking loss of breast muscle when compared with the L group. The muscle disease such as PSE (Pale, Soft and Exudative) meat [36] and wooden breast [37] have higher cooking loss than normal meat.
Stress is an essential cause of the decline in meat quality. In this study, the activity of LDH in the H group was higher than that in the L group. In transcriptome analysis, the enriched genes in the H group were involved in stimuli response pathway. In the H group, genes encoding nitric oxide synthase 1 (NOS1), Kelch-Like ECH-associated protein 1 (KEAP1) and cyclin-dependent kinase inhibitor 1A (p21, Cip1) (CDKN1A) were up-regulated. High levels of NO reduce the antioxidant capacity of post-mortem muscles, increasing the accumulation of ROS and reactive nitrogen, resulting in high levels of protein oxidation. Studies have shown that inhibition of nitric oxide synthase can signi cantly reduce protein carbonyl content and protein oxidation [38]. Inhibition of CDKN1A expression by miRNAs promotes myoblast proliferation [39]. Upregulation of KEAP1 expression increases the degradation of Nrf2 in cells, making cells more susceptible to free radical damage [40]. Heat stress can reduce the oxidative stability of broiler muscle protein and reduce the strength of the myo brillar gel, resulting in increased drip loss and cooking loss in broilers [41]. A study has shown that genes involved in the stimulation response pathway are signi cantly enriched in muscles with high drip loss [42]. Therefore, increased expression of stress pathway-related genes such as KEAP1 and CDKN1A may be one of the causes of muscle quality deterioration.
This study found that the H group had the fastest pH decline rate. The rapid decline in pH is usually accompanied by an increase in the rate of glycolysis and the accumulation of lactic acid, resulting in a decrease of muscle function [43]. In this study, high stocking density led to upregulation of genes involved in glycolysis and fat metabolism pathways. Anaerobic glycolysis is a vital energy metabolism pathway for postmortem broilers. Under anaerobic conditions, muscle glycogen degradation occurs through glycolysis, which causes pyruvate to synthesize lactic acid, thus leading to a decrease in muscle pH due to the accumulation of lactic acid [44,45]. High stocking density in this study also caused up-regulation of striated muscle contraction pathway-related genes such as SIX homeobox 1 (Six1). It has been found that white streak muscles have up-regulated expression of striated muscle contraction-related genes compared with normal meat [46]. Six1 converts slow muscle bres into fast muscle bres [47,48]. The proportion of fast muscle bres was negatively correlated with post-mortem pH [49]. Besides, the enriched genes in the H group were involved in calcium transport, sodium transport, and cation transport. Importantly, ion balance is the basis for maintaining normal physiological functions. Abnormal metabolism caused by high concentrations of calcium ions may be associated with the incidence of turkey PSE [50]. Furthermore, changes in muscle cation homeostasis may mark the beginning of muscle degeneration [51] and cause a reduction in meat quality [52].
Dietary supplementation with niacin (nicotinamide precursor) at 60 mg/kg was reported to reduce the drip loss of breast muscles in broilers [14]. In our study, the COMB group showed signi cantly reduced drip loss and cooking loss compared with the H group. Further, the COMB group showed signi cantly decreased activity of LDH compared to the H group. Besides, the COMB group showed inhibited expression of glycolytic and in ammation genes [43].
In KEGG enrichment analysis, the enriched genes in the H group were involved in in ammatory mediator regulation of RTP channels and chemokine signalling pathway. In contrast, the up-regulated genes in the COMB group were involved in the in ammatory response. Macrophage in ltration in the pectoral muscle might cause muscle damage [53]. The muscle disease such as white striped muscle is usually accompanied by elevated expression of immune-related genes [46]. During tissue degeneration, immune cells immediately enter the site of injury, triggering an in ammatory response, and attracting more immune cells to the damaged area. It can cause phagocytosis of cell debris and release of cytokines, prostaglandins and other signalling proteins, resulting in interstitial spaces [54].
We found that key genes down-regulated in the H group, such as MYOZ2, were involved in muscle development, cell adhesion, cell matrix, collagen, and cytoskeleton. MYOZ2 belongs to sarcomeric family and links calcineurin to alpha-actinin at the Z-line of skeletal muscle sarcomere and can play a role in skeletal muscle differentiation and growth [55]. It was suggested that MYOZ2 knockout mice had neuromuscular disease [56]. Also, genes down-regulated in the H group were involved in cell matrix and collagen pathways. Extracellular matrix (ECM) is a major macromolecule in skeletal muscle and has a substantial effect on meat quality. The remodelling of ECM is mainly regulated by matrix metalloproteinases. The expression of matrix metalloproteinase-1 is negatively correlated with cooking loss and positively correlated with hydraulic performance [57]. Collagen is an abundant connective tissue protein that is an important factor in the tenderness and texture of the meat and is well resistant to physical damage during cooking [58]. The addition of collagen increases the ability of pork [59] and poultry [60] to combine with water and reduces cooking losses. Furthermore, high stocking density downregulates cell adhesion, cytoskeletal and integrin binding-related genes such as integrin subunit alpha 8 (ITGA8), integrin subunit beta 8 (ITGB8) and integrin subunit beta like 1

Conclusion
High stocking density may cause oxidative stress, abnormal muscle contraction, and abnormal metabolism of glycolipids; destroy ion channels and cell matrix; reduce muscle strength by inhibiting muscle development, and cell adhesion and collagen synthesis, all of which result in reduced muscle function. Supplementation with NAM and BA in combination can improve mitochondrial function and antioxidant capacity, and inhibit in ammatory response and glycolysis by promoting muscle development and hyaluronic acid synthesis, thereby reducing drip loss of the breast muscle and improving muscle quality ( Figure 11).

Methods
Experimental birds, diets, and management A total of 300, 21-day-old Cobb broilers, were randomly divided into 3 groups: low stocking density (L, 14 birds/m 2 ), high stocking density (H, 18 birds/m 2 ) and combination of NAM and BA (COMB, 18 birds/m 2 ), with 6 replicates for each group. The L and H groups were fed a basal diet. The COMB group was fed basal diet supplemented with 50 mg/kg NAM and 500 mg/kg BA. NAM (99% purity, Jiangxi Brothers Medicine Co. Ltd., China) and BA (encapsulated, 30% effective content, Hangzhou King Technology Feed Co. Ltd., China) were purchased from the market. The composition and nutrient levels of basal diet are shown in Table 1. Experimental diets were formulated to meet or exceed the minimum nutrient requirements recommended by the National Research Council (1994) [29].
This study was conducted in an experimental chicken farm of the College of Animal Science and Technology, China Agricultural University.
Broilers were raised from 21 to 42 days of age, and feed and water were provided ad libitum. The temperature was maintained at 20-21 °C throughout the experiment while the illumination period was 18 h per day.

Sample collection
At 42-day, after 5 h of starvation, 1 broiler per replicate was randomly selected and euthanized by intravenous injection of pentobarbital sodium (390 mg/ml) at a dose of 300 mg/kg. The breast muscle of each broiler was collected and put into liquid nitrogen immediately, then stored at -80°C until further analysis. Each group had six replicates for the determination of meat quality, enzyme activities and mRNA relative expression; there were three biology replicates in each group for RNA-sequencing.

Meat quality analysis
After slaughtering, the right side of the major pectoral muscle was quickly removed for meat quality evaluation, including drip loss, cooking loss and pH. For the determination of drip loss, approximately 10 g muscle was weighed (W1) and placed in a sealed polyethylene bag at 4°C.
The muscle was reweighed (W2) after 24 h, and drip loss was expressed as (W1-W2) / W2 * 100% [30]. Cooking loss was determined according to the method described by Cai et al. [31]. Cooking loss of samples was calculated as: (initial weight-nal weight)/initial weight × 100%. The pH values of the pectoral muscle at 45 minutes and 24 hours were measured by a pH meter (testo 205; Germany). Each sample was tested at 3 different locations (top, middle and bottom) and the average of 3 measurements was calculated.

Library preparation and Illumina Hiseq xten Sequencing
RNA-seq transcriptome library was prepared following TruSeqTM RNA sample preparation Kit from Illumina (San Diego, CA, USA) using 5μg of total RNA. Shortly, messenger RNA was isolated according to the polyA selection method by oligo (dT) beads and then fragmented by fragmentation buffer rstly. Secondly, double-stranded cDNA was synthesized using a SuperScript double-stranded cDNA synthesis kit (Invitrogen, CA, USA) with random hexamer primers (Illumina). Then the synthesized cDNA was subjected to end-repair, phosphorylation and 'A' base addition according to Illumina's library construction protocol. Libraries were size selected for cDNA target fragments of 200-300 bp on 2% Low Range Ultra Agarose followed by PCR ampli ed using Phusion DNA polymerase (NEB) for 15 PCR cycles. After quanti ed by TBS380, paired-end RNA-seq sequencing library was sequenced with the Illumina HiSeq xten (2 × 150bp read length).

Read mapping
The raw paired-end reads were trimmed and quality controlled by SeqPrep (https://github.com/jstjohn/SeqPrep) and Sickle (https://github.com/najoshi/sickle) with default parameters. Then clean reads were separately aligned to the reference genome with orientation mode using TopHat version2.0.0 (http://tophat.cbcb.umd.edu/) [32] software. The mapping criteria of bowtie was as follows: sequencing reads should be uniquely matched to the genome allowing up to 2 mismatches, without insertions or deletions. Then the region of the gene was expanded following depths of sites and the operon was obtained. Also, the whole genome was split into multiple 15kbp windows that share 5kbp. New transcribed regions were de ned as more than 2 consecutive windows without the overlapped region of genes, where at least 2 reads mapped per window in the same orientation.

Differential expression analysis and Functional enrichment
To identify DEGs (differentially expressed genes) between two different samples, the expression level of each transcript was calculated according to the fragments per kilobase of exon per million mapped reads (FRKM) method. RSEM (http://deweylab.biostat.wisc.edu/rsem/) [33] was used to quantify gene abundances. R statistical package software EdgeR (Empirical analysis of Digital Gene Expression in R, http://www.bioconductor.org/packages/2.12/bioc/html/edgeR.html) [34] was utilized for differential expression analysis. Besides, functionalenrichment analysis including GO and KEGG were performed to identify which DEGs were signi cantly enriched in GO terms and metabolic pathways at P-value ≤0.05 compared with the whole-transcriptome background. GO functional enrichment and KEGG pathway analysis were carried out by Goatools (https://github.com/tanghaibao/Goatools) and KOBAS (http://kobas.cbi.pku.edu.cn/home.do) [34].

Muscle developmental gene
The cDNA was synthesized by using a reverse transcription kit PrimeScript TM RT reagent Kit with gDNA Eraser (Perfect Real Time) (RR047A, TAKARA, Japan). Then it was stored in a -80 °C refrigerator. Fluorescence quantitative PCR was performed according to the instructions of TB GreenTM Premix Ex Taq TM (Tli RNaseH Plus) (RR420A, TAKARA, Japan). The reaction apparatus was a 7500 uorescence detection system (Applied Biosystems), and the PCR reaction conditions were as follows: after pre-denaturation at 95 °C for 30 s, 40 cycles of denaturation at 95°C for 5 s and 60 °C for 34 s. After the completion of PCR ampli cation, the dissolution curve was observed, and agarose gel electrophoresis was performed to identify whether the ampli ed gene fragment conformed to the design length, and the speci city of the ampli cation product was veri ed. The target gene and the internal reference gene beta-actin primer sequence are shown in Table 2. The results of gene expression were analyzed and compared using 2 -ΔΔCT .

Statistical analysis
The results are expressed as means with their standard error mean (SEM). One-way ANOVA was used for single factor analysis by SPSS 20.

Consent for publication
Not applicable.

Competing interests
The authors declare that they have no competing interests.     The pH values of breast muscle. Data are shown as the means ± SEM. Different letters a, b indicate that there are signi cant differences (P < 0.05) among these groups. L, low stocking density (14 birds/m2); H, high stocking density (18 birds/m2); COMB, combination of NAM and BA (18 birds/m2).

Figure 3
Principal Component Analysis (PCA) and Wayne (VEEN) analysis of gene sets. For the PCA graph, the distance between each sample point represents the distance of the sample. The closer the distance means higher the similarity between samples; for the VEEN graph, the numbers inside the circle represents the sum of the number of expressed genes in the group. The crossover region represents the number of consensus expressed genes for each groups.

Figure 4
Volcanic map of differential expression genes. The abscissa is the fold change of the gene expression difference between the two samples and the ordinate is the statistical test value of the gene expression. Each dot in the gure represents a speci c gene, the red dot indicates a signi cantly up-regulated gene, the green dot indicates a signi cantly down-regulated gene, and the grey dot is a non-signi cant differential gene.

Figure 5
The Veen diagram and the map of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways analysis of gene sets. For VEEN diagram: the sum of all the numbers inside the circle represents the total gene of the set. The number, circle intersection area represents the number of shared genes among the gene sets. For the map of KEGG metabolic pathway, the red represents the pathway of the common annotation of the genes in the gene sets of two groups. We thank Kanehisa Laboratories for providing the copyright permission of KEGG pathway maps [75].

Figure 10
The mRNA relative expression of DEGs quanti ed by quantitative reverse transcription-PCR. Data presented as means ± SEM.

Figure 11
The graphic description of the normalization effect of nicotinamide and sodium butyrate on breast muscle.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.