Dietary butyrate glycerides modulate intestinal microbiota composition and serum metabolites in broilers

Butyrate can modulate the immune response and energy expenditure of animals and enhance intestinal health. The present study investigated changes in the intestinal microbiota composition and serum metabolites of young broilers in response to 3,000 ppm butyrate in the form of butyrate glycerides (BG) via pyrosequencing of bacterial 16S rRNA genes and nuclear magnetic resonance (NMR). The dietary treatment did not affect the alpha diversity of intestinal microbiota, but altered its composition. Thirty-nine key operational taxonomic units (OTUs) in differentiating cecal microbiota community structures between BG treated and untreated chickens were also identified. Bifidobacterium was, in particular, affected by the dietary treatment significantly, showing an increase in not only the abundance (approximately 3 fold, P ≤ 0.05) but also the species diversity. The (NMR)-based analysis revealed an increase in serum concentrations of alanine, low-density and very low-density lipoproteins, and lipids (P ≤ 0.05) by BG. More interestingly, the dietary treatment also boosted (P ≤ 0.05) serum concentrations of bacterial metabolites, including choline, glycerophosphorylcholine, dimethylamine, trimethylamine, trimethylamine-N-oxide, lactate, and succinate. In conclusion, the data suggest the modulation of intestinal microbiota and serum metabolites by BG dietary treatment and potential contribution of intestinal bacteria to lipid metabolism/energy homeostasis in broilers.


Effects of dietary butyrate glycerides on the ileal microbiota.
To determine the effects of dietary treatment with butyrate glycerides on the ileal microbiota, the within-community (α) diversity was firstly assessed. The rarefaction curves including the plots of Shannon index, Chao1 index, phylogenetic distance (PD), and number of unique OTUs approached a plateau as sequence numbers increased (Fig. S1), indicating that the sequence depth was sufficient for capturing the majority of OTUs in the ileal samples. No significant treatment effects on the diversity of ileal microbiota were observed (P > 0.05). The similarities between pairs of microbial communities (β-diversity) were also examined using a principal coordinate analysis (PCoA) based on the unweighted UniFrac distance (Fig. 1A). No distinguishable clustering of the samples appeared to be evident between the control and treatment group (Fig. 1B). This was further confirmed by the analysis of similarities (ANOSIM) (R = 0.02, P = 0.42).
All OTUs with abundance ≥0.005% were taxonomically assigned with the Ribosomal Database Project (RDP) classifier at 80% confidence threshold. More than 99% of the ileal sequences were assigned to bacterial phyla. Firmicutes was the most predominant phylum (>90%) followed by Proteobacteria (4%) and Cyanobacteria (2%) (Fig. S2). At the family level, Clostridiaceae was predominant (>49%), followed by Lactobacillaceae (>21%), Enterobacteriaceae (>3%), Ruminococcaceae (>2%), and Lachnospiraceae (>1%). Over 99% of the sequences from family Clostridiaceae could not be classified further to the genus level, while Lactobacillus was the only genus found in family Lactobacillaceae. Analyses using the Metastats method 17 revealed no significant difference (P>0.05) in the relative abundance (normalized as percentage of total number of sequences) of microbiota composition between the control and treatment group except for more than 10-fold reduction of family Bacillaceae/genus Bacillus (P ≤ 0.05) in BG-treated chickens.
Effects of dietary butyrate glycerides on the cecal microbiota. The rarefaction curves showing the α-diversity of cecal microbiota are displayed in Fig. S1. The majority of OTUs in the cecal samples appeared to be captured. Similarly, there was no significant difference in the diversity of cecal microbiota between the control and treatment group (P > 0.05). The microbial community structures between the control and treatment group (β-diversity) were compared using the PCoA of unweighted UniFrac distance. Birds from the control group were distinctly separated from those in the BG group ( Fig. 2A). The first and second axis of the PCoA explained 12.50% and 25.95% variation in microbial diversity, respectively. Further analysis with ANOSIM revealed a significant difference (R = 0.57, P = 0.006) in microbial communities of cecal microbiota between the two groups ( Fig. 2B).
To identify key OTUs in the discrimination of microbial community structures of cecal microbiota between the control and treatment groups, sparse partial least squares discriminant analysis (SPLS-DA) was conducted. The SPLS-DA yielded a correct classification rate of 85% for the cecal samples. Thirty-nine OTUs were identified by SPLS-DA as the key variables in the differentiation of cecal microbial profiles of the two groups (Table 1). Among the 39 OTUs, the abundance of 14 OTUs was enhanced in the cecal samples from BG-treated birds, among which 8 OTUs belonged to genera Anaerotruncus (1), Bacteroides (5), Blautia (1), and Lactobacillus (1), and 6 OTUs belonged to families Lachnospiraceae (4), Ruminococcaceae (1), and Clostridiales Family XIII Incertae Sedis (1). On the contrary, the abundance of 25 OTUs was increased in the cecal samples from birds in the control group, among which 9 belonged to genera Bacteroides (2), Clostridium (1), Oscillospira (4), and Subdoligranulum (2); 12 OTUs belonged to families Ruminococcaceae (10) and Lachnospiraceae (2); 3 OTUs belonged to order Clostridiales; and 1 OTU was unclassified. The heat map of the 39 OTUs is shown in Fig. 3. Effect of dietary butyrate glycerides on Bifidobacterium and butyrate-producing bacteria. In the present study, genus Bifidobacterium was not detected in the digesta samples by the PCR primers universal to eubacterial 16S rRNA genes during pyrosequencing, which possibly resulted from a base-pair mismatch or low amplification efficiency of PCR primers as reported previously 18,19 . To investigate the effects of BG supplementation on the population of Bifidobacterium, a pair of PCR primers 20 specific to the xylulose-5-phosphate⁄fructose -6-phosphate phosphoketolase (xfp) gene, which is characteristic of the Bifidobacterium genus 21 , were used for subsequent sequence analysis and qPCR assays. Approximately 5,000 bifidobacterial sequences per sample with an average length of 309 nucleotides were obtained from pyrosequencing. Among the sequences identified at the species level by comparison of the gene sequence similarities 21 , 46.9%, 49.9%, and 3.2% were assigned to B. gallinarum, B. saeculare, and B. pullorum, respectively, in the cecal samples from BG-treated birds. B. saeculare was the only Bifidobacterium species identified in the control group.
The effects of dietary treatment with BG on the relative abundance of Bifidobacterium and some groups of butyrate-producing bacteria were also assessed by qPCR assays. Supplementation of BG increased the Bifidobacterium abundance in the cecal samples by approximately 3 fold (Table 2). However, the supplementation had no effects (P > 0.05) on the abundance of Clostridial cluster IV and cluster XIVa that harbour  butyrate-producing bacteria 22,23 . In addition, no significant difference (P>0.05) in the abundance of the gene encoding butyryl-CoA: acetate CoA transferase in the acetyl-CoA pathway, which is the most prevalent in bacterial production of butyrate 24 , was detected in the cecal samples between the control and treatment group.

Effects of dietary butyrate glycerides on serum metabolites.
A total of 42 metabolites were unambiguously identified by the 1 H-NMR spectroscopic analysis in the samples from both treatments, and their chemical shifts, peak multiplicity, and the corresponding 1 H NMR signal multiplicities were determined ( Table 3). The assignment of metabolites was based on previous reports [25][26][27] . The spectra of the serum samples contained resonances from amino acids, organic acids, albumin, lipids, unsaturated lipids, choline, and creatine. Multivariate statistical analysis was used to detect subtle treatment-related metabolic differences 14,28 . Principal component analysis (PCA) was performed on the 1 H NMR spectra of serum samples between BG-treated and untreated birds. The PCA score plot of the 1 H NMR serum data is shown in Fig. 4A, with each point representing an individual spectrum of a sample showing separation of BG-treated birds from the untreated chicks. Partial least squares discrimination analysis (PLS-DA)-based profiling was also employed to explore the intrinsic differences between BG-and BD-fed groups. The samples from different groups were separated and classified into two distinct clusters as shown in the PLS-DA score plot (Fig. 4B). The model parameters (R 2 X = 0.563, R 2 Y = 0.932, Q 2 = 0.608) and the validated model (permutation number: 200) indicated no over fitting (Fig. 4C). All of the results indicated the existence of differences between the two groups. Furthermore, the spectral datasets were analyzed by orthogonal partial least squares discriminant analysis (OPLS-DA). The BG-fed group was clearly separated from the BD-fed chickens as shown in the OPLS-DA scores plot (Fig. 4D) and by permutation tests and ANOVA of the cross-validated residuals (CV-ANOVA) (P ≤ 0.05). The metabolites responsible for a significant contribution to the separation of two groups are indicated in the corresponding S-plot, including low-density lipoprotein (LDL), very low-density lipoprotein (VLDL), lipids, lactate, alanine, succinate, dimethylamine, trimethylamine, choline, glycerophosphorylcholine (GPC), and trimethylamine-N-oxide (TMAO) that are numbered from 1 to 10 in Fig. 4E. The VIP statistic of  More interestingly, the concentrations of intestinal bacteria-derived metabolites also increased (P ≤ 0.05), including choline, GPC, dimethylamine, trimethylamine, TMAO and succinate (Table 4).

Discussion
The effects of butyrate on the intestinal microbiota of piglets and broiler chickens have been investigated previously by culture-based 29 and non-culture-based methods, such as restriction fragment length polymorphism (RFLP) and fluorescence in situ hybridization (FISH) 30,31 . Due to the limitation of the methods, limited knowledge has been generated from these studies. The use of pyrosequencing in the present study uncovered a more    Subdoligranulum, Holdemania, Anaerotruncus, and Bifidobacterium. Many Mollicutes members cause diseases (including colitis) in humans and animals 32 . Holdemania has also been reported to be associated with unhealthy ceca and the use of antibiotics in animals 33,34 . The reduction in the abundance of Mollicutes and Holdemania in the present study suggests the beneficial effects of dietary BG treatment on chicken intestinal health. According to Bergey's Manual of Systematic Bacteriology (the second edition), the genera Subdoligranulum 35 , Eubacterium 36 , and Anaerotruncus 37 contain 1 (Subdoligranulum variabile), 44 (Eubacterium), and 1 (Anaerotruncus colihominis) species, respectively, among which many can produce butyric and lactic acids 22,38 . The decrease in the Subdoligranulum population by BG supplementation in the present study implies a functional suppression by butyrate released from BG, possibly as the result of feedback from the high level of butyrate in the chicken cecum. One striking observation in the present study is the effect of BG on the significant increase of Bifidobacterium in both diversity and abundance, as revealed by amplification of the xfp gene with both pyrosequencing analysis and qPCR assays. Bifidobacterium spp. are well-recognized probiotic bacteria with a wide spectrum of benefits 39 . Thus, the supplementation of BG likely promotes chicken health and well-being. Recently, we reported the modulation of energy and lipid metabolism in young broiler chickens at the age of three weeks by dietary BG treatment 5 . There was up-regulation of genes involved in the biological processes for the reduction of synthesis, storage, transportation and secretion of lipids in the jejunum, and for the enhancement of the oxidation of ingested lipids and fatty acids in the liver. In particular, the gene expression of transcriptional regulators of thyroid hormone responsive (THRSP) and early growth response gene-1 (EGR-1) as well as some in the peroxisome proliferator-activated receptors (PPAR) signaling pathway was significantly affected by BG. Moreover, serum triglycerides and total cholesterol were lowered in BG-fed birds. While lipoprotein lipase was decreased in the jejunum, liver and adipose of BG-fed birds, fatty acid synthase levels were lowered in the serum, liver and adipose tissue of the same birds. In animals, energy and lipid metabolism is a complex process, which is regulated by a wide array of interdependent factors, including nutrients, hormones, nuclear transcription factors, and corresponding enzymes 40,41 . A change of the profiles of circulating metabolites may thus partially reflect the effects of dietary treatments on energy and lipid metabolism. There are two pathways, through which metabolites are involved in energy catabolism, including the pyruvate (alanine) production pathway and the "energy shift" pathway 42 . Alanine is most commonly produced by reductive amination of pyruvate with the consumption of ATP 43 ; the higher serum concentration of alanine in BG-fed birds in the current study indirectly indicates an increase of pyruvate consumption as well as energy expenditure 44 . Another interesting observation is the relative higher concentrations of LDL/VLDL and lipids in the BG-fed birds than those in the control group, suggesting that BG promoted the transportation of lipids; thus reducing the storage of fatty acids in peripheral tissue 14,45 .
The intestinal microbiota can also affect host metabolism 46 . In the current study, dietary supplementation of BG increased the serum concentrations of choline, GPC, dimethylamine, trimethylamine, TMAO, and succinate when compared to those of control birds. These metabolites have been suggested to be derived from intestinal bacteria, particularly from Bifidobacterium, and are associated with several metabolic pathways 46 . For example, choline, dimethylamine, trimethylamine, and TMAO are mainly involved in the choline metabolic pathway to modulate lipid catabolism and glucose homeostasis [47][48][49][50][51] . It has been reported that major adverse cardiovascular events are normally accompanied with an increased level of TMAO 52 . It is unclear if the serum TMAO level observed in the present study was sufficient to trigger corresponding diseases in broilers. The GPC is a putative acetylcholine precursor that potentially increases growth hormone secretion and enhances fat oxidation 53 . Succinate is a key intermediate in microbial propionate synthesis 54 , and is also identified as a substrate for intestinal gluconeogenesis, a biological process that improves glucose homeostasis 55 . Collectively, enhanced serum concentrations of these metabolites by BG supplementation indicate that chicken intestinal bacteria can also contribute to lipid and energy metabolism through their metabolites to benefit the host.
Discovery of the relationships among the intestinal microbiota composition, microbiota-dependant metabolism, and animal performance can provide opportunities to improve food animal production. There have been studies on a link between bacterial species in the chicken intestine to the energy metabolism of broilers 56,57 . Some bacterial groups in the ileum and cecum 58,59 or feces 60 were related to feed conversion efficiency in chickens. For example, OTUs representing 26 bacterial species or phylotypes related to Lactobacillus spp., Ruminococcaceae, Clostridiales, Gammaproteobacteria, Bacteroidales, Clostridiales/Lachnospiraceae, and unclassified bacteria/clostridia in the ileum and cecum were found to be associated with feed efficiency of broiler chickens 58 . Recently, we reported that the supplementation of BG increased feeding efficiency by 10% and reduced abdominal fat deposition in 3-week-old broilers 5 . In the present study, 39 OTUs in the cecum of young birds were identified, which demonstrated changes in the relative abundance responding to BG treatment. Among these OTUs, Ruminococcaceae (18 OTUs) was the most predominant family followed by Bacteroidaceae (7 OTUs), Lachnospiraceae (6 OTUs), and order Clostridiales (5 OTUs). Only one OTU was assigned to Lactobacillus. These data support the previous observation by Torok et al. 58 . Bifidobacterium and Faecalibacterium prausnitzii are two groups of intestinal bacteria producing choline metabolites that can modulate lipid metabolism and glucose homeostasis 46 , and thus possibly alter animal health status and performance. The observations on the reduction of abdominal fat deposition in broilers in our previous study, and on the increase in serum concentrations of choline metabolites and in the diversity and abundance of Bifidobacterium in the chicken intestine by BG supplementation in the present study suggest that Bifidobacterium may have contributed to the decrease of abdominal fat deposition through the influence over the production of choline metabolites.

Methods
Animals, diets, and sample collection. The protocol for the animal trial was approved by the Animal Care and Use Committee of University of Guelph (AUP No. 3176) in accordance with the Canadian Council on Animal Care's Guidelines. The feed ingredients and levels, including both BG additives and monensin, were the same as those used in the previous study [5]. The BG additives were mono-butyrin (mono-C4) and a mixture of 30% mono-, 50% di-, and 20% triglycerides of n-butyric acid (Baby C4), commercially available from SILO (Industria Zootecnica, Florence, Italy). No antibiotics were used throughout the trial except for monensin for the purpose of coccidiosis prevention.
Forty eight 1-d-old male birds (Ross 308; Stratford Chick Hatchery, Stratford, ON, Canada) were allocated equally into two dietary treatments: 1) basal diet-fed group (control), and 2) BG diet (basal diet + BG)-fed group. There were 24 birds per treatment with four birds each unit on the floor. The birds in the control group were fed a commercial diet for starter phase (0 to 20 d), while those in the BG diet-feed group consumed their assigned diet in a two phases program, namely a starter diet containing 3,000 ppm each of Mono C4 and Baby C4 from 0 to 7 d, and 3,000 ppm of Mono C4 only from 8 to 20 d. Feeding program, environmental temperature, and lighting schedules were the same as previously reported 5 . Feed and water were freely available.
All the birds appeared healthy and grew well throughout the entire experimental period. On d 20, six birds were randomly selected from each treatment (one bird per cage), and approximately 4 mL blood samples were collected from the jugular vein into 5 mL heparin-free vacutainer tubes (Becton Dickinson Vacutainer Systems, Franklin Lakes, NJ, USA). Samples were centrifuged at 750 g for 10 min at 4 °C, the supernatant (serum) was immediately collected and placed into test-tubes and stored at −20 °C until a NMR-based analysis. After blood collection, the same birds (six chickens per treatment) were killed by cervical dislocation, and digesta samples from the cecum and a 10-cm-long section of the ileum (5 cm away from the ileocecal junction) of each bird were aseptically collected, and stored at −80 °C until further analysis. There were 24 digesta samples collected in total with 12 samples from each treatment (6 birds per treatment with 2 segmental digesta samples of each bird).
1 H-NMR spectroscopic analysis of serum. 1 H-NMR spectroscopic analysis of chicken serum samples was performed by the NMR Laboratory (nmr@chemistry.mcmaster.ca) in the Department of Chemistry and Chemical Biology, McMaster University (Hamilton, Ontario, Canada). One hundred microlitres of 0.9% saline in D 2 O with 1 mg/ml sodium [2,2,3,3-d4]3-trimethylsilyIpropanoate (TSP-d4) was mixed with 500 µl serum in high quality 5-mm NMR tubes. 1 H NMR spectra of serum were recorded on a Bruker Avance III 700 MHz NMR spectrometer (Bruker Biospin, Rheinstetten, Germany) equipped with a 5 mm QNP cryo-probe and SampleJet autosampler, and operating at a proton frequency of 700.17 MHz. A carr-Purcel-Meiboom-Gill (CPMG) spin-echo pulse sequence [recycle delay−90°− (τ-180°-τ) n -acquisition] was used to emphasize resonances from low molecular-weight metabolites 14,28,61 . 1 H NMR data for each sample were acquired using 128 scans (64k data points) with a 2.5 second relaxation delay. Chemical shifts were referenced to the internal reference (TSP-d 4 : 0.00 ppm). 1 H NMR data with water suppression using excitation sculpting with gradients were subsequently acquired using the same number of scans and time of relaxation delay as that for 1 H NMR data.
Analysis of NMR data. The NMR spectra of serum samples were Fourier-transformed, phase adjusted, and baseline corrected using Mnova-6.1.1 (Mestrelab, Santiago de compostela, Spain) as previously described 14 .
Chemical shifts were referenced to the internal reference (TSP-d 4 : 0.00 ppm). Each spectrum (δ 0.8-8.5) was segmented into contiguous segments having an equal width (0.01 ppm) and integrated over the region from equal width. The region δ 4.69-5.20 was removed to avoid the influence of the water signal. The integral of each region was determined. Resultant data sets were then imported into SIMCA-P 13.0 (Umetrics, Sweden) for multivariate statistical analysis. The PCA was applied to discern the presence of inherent similarities of spectral profiles and identify possible outliers within the dataset 62 . The PLS-DA was conducted to distinguish BG and BD diet-fed groups in a supervised manner. Parameters for model fitness (R 2 ) and predictive ability (Q 2 ) with leave-one-out cross validation and the response of the permutation test (200 times) needed to be used to evaluate whether the model could be established because of the small number of samples. Furthermore, a supervised pattern recognition approach known as OPLS-DA was used to improve the classification of the BG and BD diet-fed group while screening biomarkers. With an aim to discover the potential variables contributing to the differentiation, we generated an S-plot for the OPLS-DA model to define metabolites significantly contributing to the separation of the two groups. On the basis of the variable importance in the project (VIP) threshold of 1 (VIP ≥ 1.00), a number of metabolites responsible for the difference in metabolic profiles of the two groups could be obtained. In parallel, the metabolites identified by the OPLS-DA were validated at a univariate level using t-test (SPSS 17.0) with the critical P-value of 0.05 in order to detect the main metabolites that were significantly different in leading to the class discrimination 63 .
Bacterial DNA extraction. Each digesta sample was processed individually for DNA extraction and each individual DNA sample was used for subsequent analyses of pyrosequencing and qPCR assays. Bacterial DNA was extracted using QIAamp DNA Stool Mini Kit (Qiagen, Valencia, CA, USA) in accordance with the manufacturer's instructions. Briefly, approximate 0.2 g of digesta sample was lysed by incubation at 95 °C in Buffer ASL. PCR inhibitors and DNA-degrading substances were adsorbed to InhibitEX. Proteins were digested by incubation with Proteinase K at 70 °C. DNA was bound to the QIAamp silica-gel membrane and impurities were washed away. Purified DNA was eluted from the QIAamp spin column. The extracted bacterial DNA was stored at −80 °C until further analysis.
Quantitative PCR assays. The qPCR assays were used to determine the relative abundance of genus Bifidobacterium and some bacterial groups relating to butyrate production, including Clostridial clusters IV and XIVa and bacteria harboring the gene of butyryl-CoA:acetate CoA transferase 22 . The assays were performed with individual digesta DNA samples separately using a S1000 Thermocycler (1852148, Bio-Rad Laboratories, Hercules, CA, USA). The following was the thermal cycle profile: 95 °C for 3 min; 40 cycles of 95 °C for 15 sec, a specific annealing temperature for 30 sec, and 72 °C for 30 sec. A thermal melt curve was generated by heating at 95 °C for 1 min, 55 °C for 30 sec, and ramping back to 95 °C in 0.5 °C increments. The 25-μl reaction mixture contained 1.0 μl of DNA template, 12.5 μl of 2× iTaq SYBR Green Supermix (Bio-Rad, Hercules, CA, USA), and forward and reverse primers (Table S1). The abundance of the different groups of bacteria was normalized using the amplicon from the universal PCR primers towards eubacteria 64 , which served as an internal housekeeping control. Fold changes in the abundance of different group bacteria were calculated by the 2 −∆∆Ct method 65 . The PCR primers were those used for Bifidobacterium 20 , Butyryl-CoA:acetate CoA transferase 66 , Colostridial cluster IV 67 , and Colostridial cluster XIVa 68,69 , respectively (Table S1).
Bioinformatics and statistical analysis of microbiota. The data generated from each sample (either pyrosequencing or qPCR assays) were analyzed separately before combination for the analysis of treatment effects. Raw sequence data were denoised using the software USEARCH 70 . Detection and removal of chimeras were performed using the software UCHIIME 71 . Sequences with low quality and less than 250 bp in length were removed. Sequences passing the quality control screening were processed by the Quantitative Insights Into Microbial Ecology (QIIME) pipeline (http://qiime.sourceforge.net). Barcode and primers were removed during the demultiplexing step. The sequences were grouped into OTUs at least 97% similarity using UCLUST. For the 16S rRNA gene, a representative sequence was picked from each OTU and assigned taxonomy using the Ribosomal Database Project (RDP) naïve Bayesian classifier 72 with the confidence threshold as 80%. The sequences were aligned against the Greengenes-imputed core reference alignment using PyNAST, and the concatenated alignment of OTU was filtered to remove gaps and hypervariable regions using the Greengenes Lane mask. The filtered sequences alignment was then used to build a phylogenetic tree for calculation of UniFrac diversity. Data were rarefied to the lowest counts of sequences per sample for calculation of alpha and beta diversities. The PCoA based on the unweighted UniFrac distance was implemented in the software QIIME. For the sequences of the xfp gene, taxonomic assignments of the OTUs were made by comparison to a reference database of the xfp gene 21,73 . The ANOSIM provides a way to test whether similarities within groups are higher than those between groups, thus allowing for testing whether bacterial communities were different between two or more groups 74 . This nonparametric test was implemented in R package vegan (version 2.0-4), and UniFrac distance was used as a measure of dissimilarity of bacterial communities in the ileal or cecal digesta 75 . When the ANOSIM test indicated a difference of bacterial communities between the control and BG-treated group, SPLS-DA was utilized to select OTUs that could be used for separation of the two groups using the R package spls 76,77 . Sparse partial least squares regression can achieve variable selection and dimension reduction simultaneously. For our data, the binary response was the dietary treatment, i.e. 6 birds were selected from each of the control and BG-treated group. The relative abundance of OTUs was normalized as percentage of total number of sequences for an individual bird and related to the corresponding dietary treatment. A heat map of the relative abundance of selected OTUs in the cecal digesta was created by the R package gplots. Relative abundances of taxonomy and OTU were compared between the control and BG-treated group using a permutation test via online software Metastats 17 . The qPCR data were subjected to simple t-test 78 . Significance level was set at 0.05.