- Research article
- Open Access
Rearing pattern alters porcine myofiber type, fat deposition, associated microbial communities and functional capacity
BMC Microbiologyvolume 19, Article number: 181 (2019)
The Chinese believe that the meat of pigs reared in the past with free range tastes better than that of the pigs reared indoor on a large scale today. Gastrointestinal microflora is closely associated with the main factor of meat flavour, including fibre characteristics and lipid metabolism. Our method in this study involved different raising patterns within the semi free-grazing farm (FF) or indoor feeding farm (DF), the measurement of fat deposition and myofiber type by paraffin section and reverse transcription polymerase chain reaction and the identification of microbiome and functional capacities associated with meat quality through metagenomic sequencing.
Results showed that the fat area in muscle and adipose tissue and the myofiber density significantly increased in the pigs of the FF group. The relative abundance of bacteria associated with lipid metabolism, such as g_Oscillibacter, in the feces of the FF group was higher than that in DF group, and the relative abundance of some bacteria with probiotic function, including g_Lactobacillus and g_Clostridium, was lower than that in DF group. The abundance of g_Clostridium was significantly positively correlated with the intramuscular fat area, whereas health-related bacteria, such as g_Butyricicoccus, g_Eubacterium, g_Phascolarctobacterium and g_Oribacterium, was significantly negatively correlated with abdominal fat area, myofiber density and adipose triglyceride lipase (ATGL) mRNA expression. KEGG analysis showed that pigs raised in semi free-grazing farm can activate the pathway of inosine monophosphate (IMP) biosynthesis, glycolysis/gluconeogenesis and alanine, aspartate and glutamate metabolism.
Free range feeding improves meat quality by changing the fibre type, myofiber density and metabolic pathways related to flavour amino acids, IMP or glycolysis/gluconeogenesis in muscle. However, prolonged feeding cycle increases fat deposition and associated microbial communities.
Except for the influences of nutrition, feeding management and slaughter mode, muscle fibre and intramuscular fat are the main factors affecting muscle quality, such as tenderness, juiciness and flavour. The number and volume of muscle fibres mainly determine the meat quantity, and the type of muscle fibre is closely related to the pH change, colour, and water-holding capacity of post-mortem muscles [1,2,3]. The content of intramuscular fat is closely related to tenderness, juiciness and flavour of the meat .
Intestinal flora and its metabolites are closely associated with the immune system, bone, cardiovascular system and nervous system and thus affect the physical and mental health [5, 6]. Most research on intestinal microorganism and lipid metabolism focus on metabolic disorders and obesity. Adult germ-free mice with microbiota harvested from the distal intestine of conventionally raised animals produces a 60% increase in body fat content and insulin resistance despite reduced food intake . Gut microbiota is an important factor affecting energy harvest from the diet and energy storage of the host [8, 9]. Yan et al. (2016) demonstrated significant difference of the gut microbiota, which contributed in regulating fibre characteristics and lipid metabolism in skeletal muscle between obese and lean pigs . Specific gut microbiome in obese individuals can enhance ectopic fat deposition in skeletal muscle and inhibit muscle growth . These findings provide new approaches to intervene with host metabolism and animal phenotypes. Intestinal flora is an independent environmental factor for regulating host metabolism, and its changes are closely related to the host phenotype. However, the biological markers of intestinal microorganism related to phenotypic traits, such as muscle fibre and fat, have not been determined. A reference gene catalogue of pig gut microbiome has been established by Xiao et al. (2016) . Their study would accelerate research about the complex interactions between intestinal microbiota and host physiological function.
Almost all the Chinese believe that pork raised in the past with free range tastes much better than pork raised indoor on a large scale today. On June 18, 2016, the famous international media New York Times published large-scale reports on Linan Sun Commune and China’s rural construction. They followed past pig raising model, where pigs were running in the yard, swimming in the pool, rolling in the mud, and using their mouths to arch the fruit tree fences. Hence, in this study, pigs in the same litter were raised separately. One part of the pigs was raised in traditional indoor farm, while the others were raised above Linan Sun Commune with the semi free-grazing model. The main factors affecting meat flavour, including fibre characteristics and lipid metabolism, and the correlation between porcine gut microbiome and muscle fibre or fat deposition would be tested, which will provide a theoretical basis for regulating pork quality through intestinal microflora and its metabolites.
Measurement of fat deposition and myofiber type traits
The intramuscular fat area and fibre number in the longissimus dorsi muscles of animals in the semi free-grazing farm (FF) group were significantly (p < 0.01) bigger than that in the indoor feeding farm (DF) group (Fig. 1a). Back, abdominal and lard fat area of pigs in the FF group were significantly (p < 0.01) bigger than that in the DF group (Fig. 1b). The mRNA expression of MyHC1, MyHC2b and MyHC2x in the longissimus dorsi muscles of animals in the FF group was significantly (p < 0.05) higher than that in the DF group (Fig. 1d). The percentage of MyHC2x and MyHC1 in FF group was higher than that in DF group (Fig. 1c). The mRNA expression of adipocyte development genes, including FAS, ATGL and HSL, in the longissimus dorsi muscles of animals in the FF group was significantly (p < 0.05) higher than that in the DF group (Fig. 1e). The mRNA expression of MyOD1 and SOSC3 in the longissimus dorsi muscles of animals in the FF group was significantly (p < 0.05) lower than that in the DF group (Fig. 1e).
Taxonomic profiles in pig feces
Metagenome sequencing by using HiSeq Illumina platform produced 64,035.32 Mbp of raw data from 10 faecal samples (Additional file 4: Table S1). A total of 63,409.46 Mbp clean data were obtained after quality control. A total of 1,282,750,898 bp scaftigs were obtained, which were assembled from the single or mixed sample assembly. MetaGeneMark software predicted 1,992,770 open reading frames (ORFs). A total of 1,069,209 ORFs with 644.58 Mbp lengths, which include 348,172 complete genes, were obtained after redundanting. MicroNR database was used to blast with no-redundant gene catalogue, and 78.76, 75.37, 69.83, 69.35, 58.72, 55.00 and 41.31% at the kingdom, phylum, class, order, family, genus, or species level, respectively, were annotated using the lowest common ancestor (LCA) method. The observed numbers of non-redundant genes curves for all samples of FF and DF groups are shown in Additional file 1 :Figure S1.
A total of 1646 annotated genera were found in all samples. Differences were observed at the genus level between the DF and FF groups, indicating different microbial community profiles and abundances. The most dominant genera in the DF group were Prevotella (8.59%), Clostridium (6.05%), Oscillibacter (3.74%), Lactobacillus (3.28%), Ruminococcus (3.03%), Bacteroides (2.15%), Treponema (1.84%), Streptococcus (1.44%), Alistipes (1.24%) and Methanobrevibacter (0.68%), while those in the FF group included Prevotella (8.64%), Oscillibacter (6.81%), Clostridium (4.72%), Ruminococcus (2.90%), Bacteroides (2.65%), Streptococcus (1.73%), Treponema (1.25%), Alistipes (0.60%), Lactobacillus (0.48%) and Methanobrevibacter (0.10%) (Fig. 2a).
A total of 7016 species were annotated from the MicroNR database, accounting for an average of 9.71% or 11.40% of the total NR from the DF or FF group, respectively (Additional file 5: Table S2). Firmicutes bacterium CAG:110 exhibited high abundances and was dominant in sample from DF (3.30%)) and FF (2.84%). The species were followed in abundances by Lactobacillus johnsonii (1.19%), Prevotella sp. P5–92 (0.89%) and Prevotella sp. P2–180 (0.74%) in the DF group and by Ruminococcus sp. CAG:177 (1.65%), Oscillibacter sp. PC13 (1.35%) and Clostridium sp. CAG:138 (1.34%) in the FF group (Additional file 2 :Figure S2).
Principal component analysis (PCA) and ANOSIM distances were utilised to visualise the differences in taxa composition between groups. The PCA cluster showed an obvious separation between the DF and FF group (Fig. 2b), and the ANOSIM enhanced this dissimilarity (R = 0.373, p < 0.05) in the species level (Fig. 2c), wherein identical communities are given R statistic near 0, and completely distinct communities are given a value of + 1. LDA indicated that the bacteria were significantly reduced in the DF group compared with those in the FF group, which was composed of s_Succinatimonas_sp_CAG_777 in species level and g_Megasphaera and g_Anaerovibrio in genus level (Fig. 4c). The bacteria were significantly increased in the DF group compared with those in the FF group, consisting of s_Mycoplasma_sp_CAG_472, s_Eubacterium_sp_CAG_841, s_Firmicutes_bacterium_CAG_176, s_Eubacterium_coprostanoligenes, s_Clostridium_sp_CAG_533, s_Oscillibacter_sp_CAG_241_62_21, s_Firmicutes_bacterium_CAG_83, s_Ruminococcus_flavefaciens, s_Lactobacillus_reuteri and s_Lactobacillus_johnsonii in the species level (Fig. 2d).
Correlation between the phenotypic data and microbial abundance
We used the canonical correlation analysis (CCA) to detect the correlation between the eight selected environmental variables and the overall microbial community (Fig. 3a). CCA analysis showed that the abdominal fat area (p < 0.05), back fat area (p = 0.1) and MyHC2b (p = 0.1) had the greatest effect on microbial community composition, followed by myofiber density, leaf fat area, intramuscular fat area, ATGL mRNA, SOSC3 mRNA and MyOD1 mRNA. After applying CCA, all significant correlation between the selected phenotypic data and gut microbial by using the spearman rank correlation at the genus level are shown as a heat map in Fig. 3b. The effect size and direction of the correlation was presented by the fold change value and colour. In total, 27 significant correlations between individual variables and gut microbial were observed in multivariate analysis, including 8 positive correlations (red square) and 19 negative (blue square) correlations.
At the genus level, a significant positive correlation was observed in the intramuscular fat area of g_Clostridium (R = 0.65, p = 0.04), and a significant negative correlation was observed in the intramuscular fat area of g_Streptococcus (R = − 0.93, p = 0.0001) (Fig. 3b). A significant positive correlation was observed in myofiber density and the abdominal fat area of g_Sphaerochaeta (R = 0.72, p = 0.01; R = 0.64, p = 0.04) and g_Succinatimonas (R = 0.70, p = 0.02; R = 0.66, p = 0.04), whereas a significant negative correlation was observed in myofiber density and abdominal fat area of g_Oribacterium (R = − 0.64, p = 0.05; R = − 0.72, p = 0.02), g_Phascolarctobacterium (R = − 0.72, p = 0.02; R = − 0.87, p = 0.001), g_Eubacterium (R = − 0.70, p = 0.02; R = − 0.68, p = 0.03) and g_Butyricicoccus (R = − 0.76, p = 0.01; R = − 0.85, p = 0.001). A significant negative correlation was observed in ATGL mRNA expression of g_Oribacterium (R = − 0.77, p = 0.009), g_Butyricicoccus (R = − 0.68, p = 0.03), g_Faecalibacterium (R = − 0.67, p = 0.02) and g_Butyrivibrio (R = − 0.68, p = 0.03) (Fig. 3b). A significant positive correlation was observed in the myofiber density of g_Parabacteroides (R = 0.77, p = 0.009) (Fig. 3b). A significant positive correlation in MyOD1 mRNA expression (R = 0.73, p = 0.01) and a significant negative correlation in leaf (R = − 0.77, p = 0.009) and back fat area (R = − 0.83, p = 0.003) were observed on g_Angelakisella (Fig. 3b). A significant positive correlation in MyOD1 mRNA (R = 0.72, p = 0.02) expression and a significant negative correlation in leaf (R = − 0.72, p = 0.02) and abdominal fat area (R = − 0.68, p = 0.03) were observed on g_Eubacterium (Fig. 3b).
Functional composition of the metagenome analysis
Functional classification of annotated NR genes by using Kyoto Encyclopedia of Genes and Genomes (KEGG) revealed a predominance of pathways related to metabolism (carbohydrates, amino acid and nucleotide metabolism), genetic information processing (translation) and environmental information processing (membrane transport) as shown in Additional file 3 :Figure S3. PCoA clustered the two feeding models (DF vs. FF) distinctly with two PCoA axes describing 89% of the total variation, wherein 81.33% of variation was in the first principal axis (Fig. 4a). Chief differences between groups were reads that were functionally annotated as being involved in carbohydrate and amino acid metabolism. Interesting, samples from the FF group were significantly enriched for pathway to alanine, aspartate and glutamate metabolism by the action of glutamate dehydrogenase (K00262, EC:188.8.131.52), glutamate synthase and NADPH large chain (K00265, EC:184.108.40.206) and 2-oxoglutarate/2-oxoacid ferredoxin oxidoreductase subunit alpha (K00174, EC:220.127.116.11, 18.104.22.168). Moreover, the samples remarkably enriched pathway to glycolysis/gluconeogenesis by the action of aldose 1-epimerase (K01785, EC:22.214.171.124) and 2,3-bisphosphoglycerate-independent phosphoglycerate mutase (K15633, EC:126.96.36.199) and remarkably enriched the module to inosine monophosphate (IMP) biosynthesis by the action of phosphoribosylformylglycinamidine synthase (K01952, EC:188.8.131.52) and amidophosphoribosyltransferase (K00764, EC:184.108.40.206) (Fig. 4b).
Raising pattern alter fibre density, MyHC I proportion and intramuscular fat area
The histological characteristics of muscle fibre are closely related to meat quality (flavour, juiciness and tenderness). Rearing conditions lead to substantial changes in the MyHC transformation as evidenced by the different proportions of myofiber types and differences in their myofiber numbers. Muscle tissue is composed of myofibers that can be subdivided into four interchangeable types, marked by the expression of four myosin heavy chain (MyHC) isotypes: one slow (MyHC I) and three fast MyHC isoforms (MyHC IIa, IIx and IIb) . The shortening velocities increase in the following order: I < IIa < IIx < IIb  Muscle MyHC composition can be influenced by additional factors, such as animal nutrition, physical activity, age and environmental temperature . Transition between MyHCs expressions follows reversible obligatory rules: I ⇄ IIa⇄IIx⇄IIb, and the type of muscle fibres changes to left or right with the external conditions .
In comparison with the DF group, percentage changes of MyHCs in the FF group were only found between oxidative myofibers: I ← IIa → IIx, whereas no changes were observed in the percentage of MyHC IIb. The main reason for the increased percentage of slow-twitch type I myofibers in the FF group in our study may be explained by the different raising condition. As a comprehensive factor, raising pattern had a substantial influence on muscle fibre composition. In our study, pigs in the DF group were kept under indoor breeding conditions and fed with a conventional commercial diet ad libitum, while the pigs in FF group were raised in a large yard and were highly physically active daily. This result was consistent with the report of Fazarinc et al. (2017) , wherein the percentage of MyHC-I was higher in wild pigs breeding outside, which were physically active than domestic pigs indoors. Higher percentage of slow-twitch type I myofibers was also found in wild pigs kept in group housing with diminished activity and was fed ad libitum  or was obtained from a zoological garden . MyHC differentiation was accompanied by the myofiber density. Fazarinc et al., (2017)  found myofiber hypertrophy in domestic pig of all myofiber types, which was especially more intense than in the wild pig feeding outside. The same results were found in our study, wherein the fibre number in the longissimus dorsi muscles of animals in the FF group was significantly higher than that in the DF group.
The composition of skeletal muscle fibre type was related to meat quality, lipid content and distribution in muscle, which directly affects muscle quality . The small muscle fibre area and high density are consistent with the prominent meat quality. The diameter of type I muscle fibres is small, and the number of muscle fibres per unit area is large. Hence, the intramuscular fat (IMF) content is high, and IMF affects the tenderness and flavor of muscle. The expression of MyHC I and IIa was positively correlated with IMF, whereas the expression of MyHC IIb was negatively correlated with IMF . Hence, an increase in IMF can enhance the eating quality of pig meat , and MyHC I myofibers and fibre density are beneficial to pork quality .
Intramuscular accumulation of lipids and myofiber type of skeletal muscles is the result of interactions between genetical background and environmental factors. Other explanation for increasing IMF and levels of type I fiber in longissimus dorsi muscles of animals in FF group maybe influence by some bioactive compounds in diet, the longer feeding cycle. For example, the bioactive compounds in mulberry leaves may play a regulatory role in fat deposition and muscle fiber types, like resveratrol, anthocyanin or 1-deoxynojirimycin . A decreased expression of fast and an increased expression of slow MyHC occurred in many physiological situations including ageing was demonstrated in studies in variety of species . Prolonged feeding cycle of pigs in the FF group increased fat deposition, as indicated by increasing back, abdominal, lard and intramuscular fat area in the longissimus dorsi muscles. In brief, pigs reared in semi free-grazing farm increased fibre density, the MyHC I proportion and intramuscular fat area in this study, which could partly explain the better taste of free range pork.
Raising pattern alter intestinal microbial
The effect of feeding pattern on intestinal microorganisms was studied. Prolonged of feeding cycle of pigs in the FF group increased fat deposition. The ratio of Firmicutes/Bacteroidetes (F/B) was associated with the energy harvest and fat deposition of host , and the F/B ratio seemed as a biomarker indicative to obesity. However, with the rapid accumulation of data by meta-analyses, a clear trend between the F/B ratio and obesity status could not be found , which suggested that the complexity of how gut microbiome modulates obesity is way more than a simple imbalance in the status of these phyla. The abundance of Oscillibacter was the major determinant of obese or normal status , when Bacteroides and Faecalibacterium were equally abundant. As a conditional pathogenic bacterium, Oscillibacter was positively associated with obesity-related metabolic pathways . Consistent with our study, g_Oscillibacter in feces was relative higher in higher fat deposition group. In the top 10 genera, the relative abundance of g_Lactobacillus was relative higher in the DF group. Linear discriminant analysis effect size (LEfSe) analysis found that s_Lactobacillus_johnsonii and s_Lactobacillus_reuteri were significantly higher bacterial species in the DF group. Lactobacillus is recognised probiotics in animal husbandry, which compete for nutrients with existing gut microbiota, which reduce body weight and fat mass . The relationship of g_Lactobacillus and obesity are controversial. Among the bacterial genera, Lactobacillus increased after a weight loss program in adolescents , but it was also abounded in obese and overweight children . Million et al. (2012)  demonstrated that s_Lactobacillus_reuteri was associated with obesity in adults. These findings suggest a possible role of Lactobacillus at the species level in body weight and obesity. From the LEfSe analysis at the species level, significantly higher bacterial species was found in the FF group s_Firmicutes_bacterium_CAG_176, s_Oscillibacter_sp_CAG_241_62_21 and s_Firmicutes_bacterium_CAG_83, which all belong to Clostridium cluster IV. Clostridium cluster IV is a dominant bacterial group known as butyrate-producing bacteria , which have beneficial functions, including nutrient absorption, epithelial cell maturation, production of short chain fatty acids and maintenance in human . The abundance and diversity of Clostridium cluster IV decreases in vegetarians . s_Ruminococcus_flavefaciens was another higher bacterial species in the DF group. Ransom-Jones et al. (2012)  found that R. flavefaciens and Clostridium cluster IV degraded dietary fibre. Above all, bacteria related probiotic function, including g_Lactobacillus and Clostridium cluster IV abound in DF group, and the relative abundance of g_Oscillibacter, which was associated with obesity, was higher in the FF group.
Relationship between faecal microbial and myofiber types or fat deposition
Reports about the relationship between faecal microbial and muscle fibre types or fat deposition were rare. Yan et al. (2016) demonstrated abundance levels of major bacteria that produce butyrate and acetic acid, such as species from Roseburia and Blautia, were higher in obese Rongchang pigs and its mouse recipients . In our study, bacterial genera g_Butyricicoccus, g__Eubacterium, g_Phascolarctobacterium and g_Oribacterium had significant negative correlation with abdominal fat area, myofiber density and ATGL mRNA expression. Butyricicoccus is a butyrate-producing clostridial cluster IV genus. The abundance of Butyricicoccus increased significantly with age . Genera Eubacterium of the gut microbiota belongs to Clostridium clusters XIVa, and are known for polysaccharide fermentation and bile acids dehydroxylation . Patients with poorer nutritional status revealed a much lower abundance of Clostridium cluster XIVa compared to healthy siblings with poorer nutritional status . Many Clostridium clusters IV and XIVa bacteria produce butyrate, a metabolite that acts as a primary source of energy for colon cells, helps eliciting an anti-inflammatory response, establishment and maintenance of the GI barrier, and reduction of intestinal permeability [39, 40]. Clostridial cluster XIVa bacteria have potential beneficial effects with respect to the development of obesity and associated metabolic disorders . Clostridium cluster IV, which has been associated with both obesity and weight loss [29, 42]. The higher levels of Clostridium cluster IV were detected in the Swedish with a high consumption of fish and meat , and lower abundance of Clostridium cluster IV in faecal microbiota of vegetarians . g_Oribacterium is a health-associated genus , which decreased in chronic obstructive pulmonary disease patients. Phascolarctobacterium may exert a beneficial role in the human gastrointestinal tract and can produce short-chain fatty acids, including acetate and propionate [45, 46]. Besides, Phascolarctobacterium spp. specialise in the utilisation of succinate produced by other bacteria. The abundance of Parabacteroides, which is a major producer of succinate, was increased by high fat diet and was positively correlated with body weight . Phascolarctobacterium faecium-like bacteria in younger individuals were maintained at a high level with a gradual increase with increasing ages (between 1 and 60 years old) but with a decrease in elderly individuals (> 60 years old) . In brief, health related bacteria had negative correlation with abdominal fat area, myofiber density and ATGL mRNA expression. The report about relationships between relevant microbiota and phenotypes of meat quality is rare, further research is needed.
Semi free-grazing feeding activated metabolic pathway related with meat flavour
Sixteen functional KO were detected to be differentially abundant in at least one of the two raising patterns. Ten of these KO were significantly enriched in FF group compared with the DF group. Interestingly, two metabolic KO (K01952 and K00764) were involved in the modules of IMP biosynthesis (M00048), which is the pathway module for phosphoribosyl pyrophosphate + glutamine = > IMP . Additionally, five KO detected in FF group were responsible for alanine (Ala), aspartate (Asp) and glutamate (Glu) metabolism (K00262, K00265 and K00174) and glycolysis/gluconeogenesis (K01785 and K15633). ABC transport pathway (K16785 and K19350), deoxyribonucleic acid (DNA)-damage-inducible protein (K07473) and transposase (K07483) were significantly overrepresented in the DF group. They are all associated with DNA repair processes.
Glu, Asp and Ala are important flavour amino acids, which can act on the precursor amino acid as a fundamental substance forming delicate flavour of meat; especially, Glu is the most pivotal amino acid affecting meat flavour and acid-base buffering capacity [49, 50]. By the KEGG pathway analysis, the Clostridium butyricum -treated group enriched alanine, aspartate and glutamate metabolism provided supporting evidence for the increased meat quality in breast muscle of ducks [51, 52]. Significant alterations were observed in the alanine, aspartate and glutamate metabolism pathway with the treatment ofβ-Methylamino-l-alanine . IMP is the key indicator for meat flavor [54,55,56]. IMP and inosine degrade into active ribose and ribose phosphate, which participate in the Maillard reaction in meat during cooking and produce various volatile flavour compositions . IMP and sodium glutamate can co-participate in the regulation of meat flavour . The genes encoding the modules of IMP biosynthesis (M00048) in the cells of Se-enriched C. utilis were up-regulated according to KEGG pathway modules analysis . More sequences related to inosine monophosphate (required for the de novo synthesis of purines) were identified from O. flexuosa transcripts than from the nearly-complete genome of B.malayi . Glycolysis/gluconeogenesis is also a muscle-development-related pathway, which is strongly connected and belongs to the most important energetic processes that influence the muscle-to-meat conversion . Drip loss strongly depends on post mortem energetic processes in the muscle. Julia et al. (2016)  found that glycolysis/gluconeogenesis significantly influences drip loss. Glycolytic potential can also predict meat quality, which is related with high drip loss and low pH24h . A report aimed to examine the specific genetic contribution of each breed to meat quality found that most of the differentially expressed genes were grouped in the Glycolysis/Gluconeogenesis pathway, which was over-expressed in Maremmana and down-expressed in Chianina . They hypothesized that the high expression levels of genes involved in Gluconeogenesis due to that the Maremmana is a breed, which adapted to the harsh environment o and to poor forage. The glycolytic pathway was enhanced in oysters with high-glycogen content, and these oysters have a high-energy metabolism, as well as an increased antioxidant capacity and stress resistance . Notably, the enrichment of KO involved in flavour amino acid metabolism, IMP biosynthesis and glycolysis/gluconeogenesis pathway in the longissimus dorsi muscle provides supporting evidence for the increased porcine meat quality in the FF group. However, more metagenome samples are needed for accurately statistical validation of these characterised KO as promising biomarkers for meat quality.
In conclusion, pigs reared in semi free-grazing farm increased myofiber density, MyHC I proportion and intramuscular fat area. It also activated the metabolic pathways related to meat flavour, including flavour amino acid metabolism, IMP biosynthesis and glycolysis/gluconeogenesis. These findings provide evidence for the good taste of free range pork. Prolonged feeding cycle of pigs in the FF group increased the conditional pathogenic bacteria associated with obesity.
Animal sample collection
Animal studies were conducted in accordance with the guidelines of the Zhejiang Farm Animal Welfare Council of China and approved by the ethics committee of Zhejiang Academy of Agricultural Sciences. Twenty Bajiazhe black piglets from Jicheng Animal Husbandry Co., Ltd. were divided into two groups of 10 at 3 months of age with equal representation of littermates and sex. Piglets in the group of the indoor feeding farm (DF) was fed with general compound feed (Additional file 6 :Table S3) as recommended by National Research Council (NRC, 1998) and were housed in Jicheng Animal Husbandry Co., Ltd. Piglets in the semi free-grazing farm (FF) group were fed with compound feed with seasonal pastures and vegetables (Additional file 6 :Table S4) and were housed in Hangzhou Sun Commune. Faecal samples in the DF group and FF group were collected from each animal at 7 and 10 months of age respectively. Feces was collected rectally in each animal by using a sterile cotton swab wet, and then the tube with feces were rapidly frozen in liquid nitrogen and stored at − 80 °C until DNA extraction. Animals in the DF group or the FF group were slaughtered in the fixed-point slaughter house according the operating procedures of pig slaughtering (GB/T 17236–2008) with approximately 100 kg body weight at 7 or 10 months of age, respectively. Approximately 1 m3 back fat or abdominal fat was collected from the junction of the sixth and seventh thoracic vertebra. Approximately 1 m3 leaf fat was collected from perirenal fat. Exactly 5-cm-thick longissimus dorsi muscle from the last third or fourth rib backwards was collected from the carcass within 30 min after exsanguination. One piece of longissimus dorsi muscle was placed in the RNA locker immediately for RNA isolation. The other piece of longissimus dorsi muscles, back fat, abdominal fat and leaf fat were fixed in 10% formalin for histology.
The muscle and fat samples were placed into 10% formaldehyde (pH 7.4) for 12–24 h and then dehydrated and subjected to routine haematoxylin eosin staining. The muscle tissue was randomly selected at 5–10 visual fields at 100× (intramuscular fat area measurement) or 200× magnification (muscle fibre number and single muscle fibre area measurement) visual field. Then, the intramuscular fat area and single muscle fibre area were measured using Image-ProPlus 5 image analysis software. The adipose tissue was randomly selected at five visual fields at a 50× magnification visual field. Then, the Image-ProPlus 5 image analysis software was used to measure the adipocyte diameter of each sample (15 cells per specimen), as well as the area of the single fat.
Reverse transcription polymerase chain reaction (RT-PCR)
Total RNA was extracted according to the RNAiso Plus instructions. According to the PrimeScript® RT-reagent Kit manual, the total RNA was used for first strand cDNA synthesis. The primers used are shown in Table S5. Real-time PCR was carried out with ABI plus one (Life Technologies, Carlsbad, CA, USA) by using SYBR Premix Ex TaqTM. Two-step PCR amplification was performed under the following conditions: 95 °C for 30 s, 40 cycles of 95 °C for 5 s and 62 °C for 34 s. Melting curve analysis was performed to verify the single-product generation at the end of the assay. Standard curves were generated based on the data obtained from the standards of the 2~2− 6 dilution series template. The amplification efficiency ranged from 90%~ 110%. Therefore, the relative level of gene expression between different groups can be calculated using the 2-△△Ct method.
Total genomic DNA from was extracted using the hexadecyltrimethylammonium bromide method. The concentration of the extracted DNA was measured using 1% agarose gel electrophoresis. DNA purity was checked using Nano Photometer® spectrophotometer (IMPLEN, CA, USA). OD260/OD280 ≈ 1.8 was a qualified DNA sample. Then, the sample was diluted to 1 ng/μL by using sterile water.
DNA library construction and sequencing
Exactly 1 μg DNA per sample was used as input material for the DNA sample preparations. Sequencing libraries were generated using NEBNext® UltraTM DNA Library Prep Kit for Illumina (NEB, USA) following the manufacturer’s recommendations, and index codes were added to attribute sequences to each sample. Briefly, the DNA sample was fragmented by sonication to a size of 350 bp. Then, DNA fragments were end-polished, A-tailed and ligated with the full-length adaptor for Illumina sequencing with further PCR amplification. At last, PCR products were purified (AMPure XP system), and libraries were analysed for size distribution in an Agilent2100 Bioanalyzer and quantified using real-time PCR. Index-coded samples were clustered on a cBot Cluster Generation System according to the manufacturer’s instructions. After cluster generation, the library preparations were sequenced on an Illumina HiSeq platform, and paired-end reads were generated.
Quality control of reads
Clean data were obtained from raw data through quality control by using Readfq (V8, https://github.com/cjfields/readfq). To prevent host DNA contamination, it needs to be compared with the host database to filter out the reads that may originate from the host by using SoapAligner software (soap2.21, https://omictools.com/soapaligner-tool) with parameters  ‘identity ≥90%, -l 30, -v 7, -M 4,-m 200,-x 400’.
De novo assembly of the short reads. Clean data were assembled using the SOAPdenovo software  (V2.04, https://omictools.com/soapdenovo-tool) with parameters  ‘-d 1, -M 3, -R, -u, -F, -K 55’. Scaftigs without N were obtained by interrupting assembled scaffolds from the joint of N. Unused reads were obtained by mapping clean data against scaftigs by using SoapAligner software (soap2.21) with parameters ‘identity ≥90%, −m 200 -× 400’ . Then, unused data were pooled using SOAPdenovo software (V2.04) with the same parameters of single sample . Mixed assembled scaffolds were split into scaftigs, and scaftigs of lengths less than 500 bp were discarded.
Gene prediction and gene abundance profiling
MetaGeneMark(V2.10, http://topaz.gatech.edu/GeneMark/) was used to predict ORFs in scaftigs from the single or mixed sample (> = 500 bp) . Information was screened using a 100 nt cut-off with default parameters. Non-redundant initial gene catalogue was obtained using the CD-HIT software (V4.5.8, http://www.bioinformatics.org/cd-hit/) from predicted ORFs with parameters ‘-c 0.95, -G 0, -aS 0.9, -g 1, -d 0’ . SoapAligner (soap2.21) was used to align paired-end clean data against the initial gene catalogue with parameters ‘-m 200, −× 400, identity ≥95’ . Unigenes were obtained from paired-end reads by removing genes with reads ≤2. We counted the gene abundance according to both numbers of paired-end reads and gene length, which was calculated as follows:
Core-pan gene analysis, correlation analysis and Venn diagram of gene numbers were employed based on the gene abundance in each sample.
Bacteria, fungi, archaea and viruses in the non-redundant (NR) library (Version: 20161115, https://www.ncbi.nlm.nih.gov/) of NCBI were used as reference database, and DIAMOND software (V0.7.9, https://github.com/bbuchfink/diamond/)  was used in sequence blast with parameters ‘blastp,-e 1e-5’. Results of evalue > minimum evalue*10 were removed . To avoid multiple alignment results, we used the LCA (https://en.wikipedia.org/wiki/Lowest_common_ancestor) was used to identify species annotation of the sequence. Then, the species abundance and gene numbers of the six-level taxonomic classification from phylum to species of each sample were obtained. Then, Krnoa analysis , relative abundance profile, heat map of abundance clustering, PCA (R ade4 package, Version 2.15.3)  and NMDS (R vegan package, Version 2.15.3) . Dimension reduction analysis was conducted based on the abundance of each classification. Difference between groups were found using ANOSIM analysis (R vegan package, Version 2.15.3). Differential species were found using Metastat and LEfSe statistical analysis. The p value of Metastat statistical analysis was obtained using permutation test of each classification between groups; then, the q value was obtained by the correcting p value by using the Benjamini and Hochberg false discovery rate method . LEfSe software was used for LEfSe statistical analysis with the LDA score 2 .
Gene functional classification
DIAMOND software (V0.7.9) was used to blast unigenes with function databases with parameters ‘blastp, -e 1e-5’. The function databases included the KEGG database (Version 201,609, http://www.kegg.jp/kegg/). Best blast hit (one HSP > 60 bits) from above blast was selected for subsequent analysis. The relative abundance of different function levels (five levels in KEGG database) was counted. Based on the abundance of different function levels, the numbers of annotated gene, relative abundance profile and heat map of abundance clustering were obtained, and PCA and NMDS dimension reduction analysis, ANOSIM analysis, metabolic pathway analysis and Metastat and LEfSe statistical analysis were carried out.
Correlation between the phenotypic data and microbial abundance. The correlation and p value between the phenotypic data and microbial alpha diversity or species were obtained using the Spearman correlation model. Windrose figures were constructed using the correlation data between phenotypic data and the top 35 genus or species abundance. Count value in figure was correlation coefficient ‘r’. r < 0 indicates significant negative correlation, whereas r > 0 indicates significant positive correlation. The absolute value of ‘r’ was found in the windrose figure. p < 0.05 was considered statistically significant.
Data of the paraffin section and real-time PCR were expressed as mean ± standard deviation. Statistical differences between groups were determined by Analysis of Variance with Duncan’s multiple comparison. p < 0.05 was considered statistically significant.
Availability of data and materials
The datasets used and/or analysed during the current study are available from the corresponding author on reasonable request. All data generated or analysed during this study are included in this published article and its supplementary information files.
canonical correlation analysis
indoor feeding farm
semi free-grazing farm
Kyoto Encyclopedia of Genes and Genomes
lowest common ancestor
linear discriminant analysis effect size
myosin heavy chain
open reading frames
principal component analysis
reverse transcription polymerase chain reaction
Ryu YC, Kim BC. The relationship between muscle fiber characteristics, postmortem metabolic rate, and meat quality of pig longissimus dorsi, muscle. Meat Sci. 2005;71:35.
Choi YM, Ryu YC, Kim BC. Influence of myosin heavy- and light chain isoforms on early postmortem glycolytic rate and pork quality. Meat Sci. 2007;76:281–8.
Choe JH, Choi YM, Lee SH, Shin HG, Ryu YC, Hong KC, et al. The relation between glycogen, lactate content and muscle fiber type composition, and their influence on postmortem glycolytic rate and pork quality. Meat Sci. 2008;80:355–62.
Davoli R, Braglia S. Molecular approaches in pig breeding to improve meat quality. Brief Funct Genomic Proteomic. 2007;6:313–21.
Dorrestein PC, Mazmanian SK, Knight R. Finding the missing links among metabolites, microbes, and the host. Immunity. 2014;40:824.
Greenhalgh K, Meyer KM, Aagaard KM, Wilmes P. The human gut microbiome in health: establishment and resilience of microbiota over a lifetime. Environ Microbiol. 2016;18:2103–16.
Backhed F, Ding H, Wang T, Hooper LV, Koh GY, Nagy A, et al. The gut microbiota as an environmental factor that regulates fat storage. P Natl Acad Sci USA. 2004;101:15718.
Turnbaugh PJ, Ley RE, Mahowald MA, Magrini V, Mardis ER, Gordon JI. An obesity-associated gut microbiome with increased capacity for energy harvest. Nature. 2006;444:1027–31.
Thaiss CA, Itav S, Rothschild D, Meijer M, Levy M, Moresi C, et al. Persistent microbiome alterations modulate the rate of post-dieting weight regain. Nature. 2016;540:544–51.
Yan H, Diao H, Xiao Y, Li W, Yu B, He J, et al. Gut microbiota can transfer fiber characteristics and lipid metabolic profiles of skeletal muscle from pigs to germ-free mice. Sci Rep-UK. 2016;6:31786.
Xiao L, Estellé Jordi Kiilerich P, Ramayo-Caldas Y, Xia Z & Feng Q, et al. A reference gene catalogue of the pig gut microbiome Nat Microbiol 2016;1:16161.
Schiaffino S, Reggiani C. Fiber types in mammalian skeletal muscles. Physiol Rev. 2011;91:1447–531.
Bottinelli R, Betto R, Schiaffino S, Reggiani C. Maximum shortening velocity and coexistence of myosin heavy chain isoforms in single skinned fast fibres of rat skeletal muscle. J Muscle Res Cell Motil. 1994;15:413–9.
Herpin P, Lossec G, Schmidt I, Cohen-Adad F, Duchamp C, Lefaucheur L, et al. Effect of age and cold exposure on morphofunctional characteristics of skeletal muscle in neonatal pigs. Pflueg Arch Eur J Phy. 2002;444:610–8.
Lefaucheur L, Milan D, Ecolan P, Callennec CL. Myosin heavy chain composition of different skeletal muscles in large white and Meishan pigs. J Anim Sci. 2004;82:1931–41.
Fazarinc G, Vrecl M, Skorjanc D, Cehovin T, Candek-Potokar M. Dynamics of myosin heavy chain isoform transition in the longissimus muscle of domestic and wild pigs during growth: a comparative study. Animal. 2017;11:164–74.
Muller E, Rutten M, Moser G, Reiner G, Bartenschlager H, Geldermann H. Fibre structure and metabolites in m. longissimus dorsi of wild boar, pietrain and Meishan pigs as well as their crossbred generations. J Anim Breed Genet. 2002;119:125–37.
Losel D, Franke A, Kalbe C. Comparison of different skeletal muscles from growing domestic pigs and wild boars. Archiv Tierzucht. 2013;56:766–77.
Karlsson AH, Klont RE, Fernandez X. Skeletal muscle fibres as factors for pork quality. Livest Prod Sci. 1999;60:255–69.
Guo J. Comparison of muscle fiber types in Jinhua and landrace pigs and the effect of ERK gene on muscle fiber transition. Doctoral dissertation: Zhejiang University; 2011.
Fernandez X, Monin G, Talmant A, Mourot J, Lebret B. Influence of intramuscular fat content on the quality of pig meat-2. Consumer acceptability of m. longissimus lumborum. Meat Sci. 1999;53:67–72.
Choi YM, Ryu YC, Kim BC. Effect of myosin heavy chain isoforms on muscle fiber characteristics and meat quality in porcine longissimus muscle. J Muscle Foods. 2006;17:15.
Zhu Z, Junjie J, Jie Y, Xiangbing M, Bing Y & Daiwen C. Effect of dietary supplementation with mulberry(morus alba l.) leaves on the growth performance,meat quality and antioxidative capacity of finishing pigs. J Integr Agr 2019;18(1):147–155.
Lina C, Silvia Q, Luana T, Carlo R, Francesco M, Luca M, et al. Age-dependent variations in the expression of myosin isoforms and myogenic factors during the involution of the proximal sesamoidean ligament of sheep. Res Vet Sci. 2019;124:270–9.
Walters WA, Xu Z, Knight R. Meta-analyses of human gut microbes associated with obesity and IBD. FEBS Lett. 2014;588:4223–33.
Hae-Jin H, Sin-Gi P, Byul JH, Min-Gyu C, Kyung-Hee P & Heon KJ, et al. Obesity alters the microbial community profile in korean adolescents. PLos One. 201;10:e0134333.
Kong C, Gao R, Yan X, Huang L, Qin H. Probiotics improve gut microbiota dysbiosis in obese mice fed a high-fat or high-sucrose diet. Nutrition. 2019;60:175–84.
Okeke F, Roland BC, Mullin GE. The role of the gut microbiome in the pathogenesis and treatment of obesity. Global Advances in Health and Medicine. 2014;3:44–57.
Santacruz A, Marcos A, Wärnberg J, Martí A, Martinmatillas M, Campoy C, et al. Interplay between weight loss and gut microbiota composition in overweight adolescents. Obesity. 2012;17:1906–15.
Bervoets L, Hoorenbeeck KV, Kortleven I, Caroline VN. Differences in gut microbiota composition between obese and lean children: a cross-sectional study. Gut Pathog. 2013;5:10.
Million M, Maraninchi M, Henry M, Armougom F, Richet H, Carrieri P, et al. Obesity-associated gut microbiota is enriched in lactobacillus reuteri and depleted in bifidobacterium animalis and methanobrevibacter smithii. Int J Obesity. 2012;36:817–25.
Louis P, Flint HJ. Diversity, metabolism and microbial ecology of butyrate-producing bacteria from the human large intestine. FEMS Microbiol Lett. 2009;294:1–8.
Woodmansey EJ. Intestinal bacteria and ageing. J Appl Microbiol. 2007;102:1178–86.
Liszt K, Zwielehner J, Handschur M, Hippe B, Thaler R, Haslberger AG. Characterization of bacteria, clostridia and bacteroides in faeces of vegetarians using qpcr and pcr-dgge fingerprinting. Ann Nutr Metab. 2009;54:253–7.
Ransom-Jones E, Jones DL, Mccarthy AJ, Mcdonald JE. The fibrobacteres: an important phylum of cellulose-degrading bacteria. Microb Ecol. 2012;63:267–81.
Xi Y, Shuling N, Kunyuan T, Qiuyang Z, Hewen D, ChenCheng G, et al. Characteristics of the intestinal flora of specific pathogen free chickens with age. Microb Pathogenesis. 2019;132:325–34.
Philippe G. Metabolism of cholesterol and bile acids by the gut microbiota. Pathogens. 2014;3:14–24.
Duytschaever G, Huys G, Bekaert M, Boulanger L, De Boeck K, Vandamme P. Dysbiosis of bifidobacteria and clostridium cluster Xiva in the cystic fibrosis fecal microbiota. J Cyst Fibros. 2013;12(3):206–15.
Devriese S, Eeckhaut V, Geirnaert A, Bossche LVD, Hindryckx P, Wiele TVD, et al. Reduced mucosa-associated butyricicoccus activity in patients with ulcerative colitis correlates with aberrant claudin-1 expression. J Crohns Colitis. 2017;11:229–36.
Canani RB, Costanzo MD, Leone L, Pedata M, Meli R, Calignano A. Potential beneficial effects of butyrate in intestinal and extraintestinal diseases. World J Gastroenterol. 2011;7:1519–28.
Neyrinck AM, Possemiers S, Verstraete W, Backer FD, Cani PD, Delzenne NM. Dietary modulation of clostridial cluster Xiva gut bacteria (roseburia spp.) by chitin–glucan fiber improves host metabolic alterations induced by high-fat diet in mice. J Nutr Biochem. 2012;23(1):51–9.
Geurts L, Lazarevic V, Derrien M, Everard A, Van Roye M, Knauf C, et al. Altered gut microbiota and endocannabinoid system tone in obese and diabetic leptin-resistant mice: impact on apelin regulation in adipose tissue. Front Microbiol. 2011;2:149.
Mueller S, Saunier K, Hanisch C, Norin E, Alm L, Midtvedt T, et al. Differences in fecal microbiota in different European study populations in relation to age, gender, and country: a cross-sectional study. Appl Environ Microbiol. 2006;72:1027–33.
Wu X, Chen J, Xu M, Zhu D, Wang X, Chen Y, et al. 16s rDNA analysis of periodontal plaque in chronic obstructive pulmonary disease and periodontitis patients. J Oral Microbiol. 2017;9:1324725.
Li L, Su Q, Xie B, Duan L, Zhao W, Hu D, et al. Gut microbes in correlation with mood: case study in a closed experimental human life support system. Neurogastroent Motil. 2016;28:1233–40.
Wu F, Guo X, Zhang J, Zhang M, Ou Z, Peng Y. Phascolarctobacterium faecium abundant colonization in human gastrointestinal tract. Exp Ther Med. 2017;14:3122–6.
Lecomte V, Kaakoush NO, Maloney CA, Raipuria M, Huinao KD, Mitchell HM, et al. Changes in gut microbiota in rats fed a high fat diet correlate with obesity-associated metabolic parameters. PLoS One. 2015;10:e126931.
Patterson D, Bleskan J, Gardiner K & Bowersox J. Human phosphoribosylformylglycineamide amidotransferase (FGARAT): regional mapping, complete coding sequence, isolation of a functional genomic clone, and DNA sequence analysis. Gene. 1999;239:0–391.
Cameron ND, Enser MB. Fatty acid composition of lipid in longissimus dorsi muscle of duroc and british landrace pigs and its relationship with eating quality. Meat Sci. 1991;29:295–307.
Heyer A, Lebret B. Compensatory growth response in pigs: effects on growth performance, composition of weight gain at carcass and muscle levels, and meat quality. J Anim Sci. 2007;85:769–78.
Liu YH, Jia YX, Liu C, Limin Ding LM, Xia ZF. RNA-Seq transcriptome analysis of breast muscle in Pekin ducks supplemented with the dietary probiotic Clostridium butyricum. BMC Genomics. 2018;19:844.
Liu YH, Li YY, Feng XC, Wang Z, Xia ZF. Dietary supplementation with Clostridium butyricum modulates serum lipid metabolism, meat quality, and the amino acid and fatty acid composition of Peking ducks. Poult Sci. 2018;97(9):3218–29.
Engskog MKR, Ersson L, Haglöf J, Arvidsson T, Pettersson C. & Brittebo e.β-n-methylamino-l-alanine (BMAA) perturbs alanine, aspartate and glutamate metabolism pathways in human neuroblastoma cells as determined by metabolic profiling. Amino Acids. 2017;49(5):905–19.
Tikk M, Tikk K, Torngren MA, Meinert L, Aaslyng MD, Karlsson AH, et al. Development of inosine monophosphate and its degradation products during aging of pork of different qualities in relation to basic taste and retronasal flavor perception of the meat. J Agr Food Chem. 2006;54:7769–77.
Vani ND, Modi VK, Kavitha S, Sachindra NM & Mahendrakar NS. Degradation of inosine-5-monophosphate (IMP) in aqueous and in layering chicken muscle fibre systems: effect of pH and temperature. LWT-Food Sci Technol 2006;39:0–632.
Yamaguchi S. Roles and efficacy of sensory evaluation in studies of taste. J Japan Soc Food Sci. 1991;38:972–8.
Arya SS, Parihar DB. Changes in free nucleotides, nucleosides and bases during thermal processing of goat and sheep meats. Part i. effect of temperature. Die Nahrung. 2010;23:1–7.
Mónica F, Grimm CC, Toldrá F, Spanier AM. Correlations of sensory and volatile compounds of spanish “serrano” dry-cured ham as a function of two processing times. J Agr Food Chem. 1997;45:2178–86.
Zhang GC, Wang DH, Wang DH, Wei GY. The mechanism of improved intracellular organic selenium and glutathione contents in selenium-enrichedcandida utilisby acid stress. Appl Microbiol Biot. 2017;101(5):2131–41.
Mcnulty SN, Sahar A, Simon GM, Makedonka M, Mcnulty NP, Kerstin F, et al. Transcriptomic and proteomic analyses of a wolbachia-free filarial parasite provide evidence of trans-kingdom horizontal gene transfer. PLoS One. 2012;7(9):e45777.
Scheffler TL, Gerrard DE. Mechanisms controlling pork quality development: the biochemistry controlling postmortem energy metabolism. Meat Sci. 2007;77:7–16.
Julia W, Christiane N, Hanna H, Mehmet C, Christian L, Karl S, et al. Integrative analysis of metabolomic, proteomic and genomic data to reveal functional pathways and candidate genes for drip loss in pigs. Int J Mol Sci. 2016;1:1426.
Zybert A, Sieczkowska H, Antosik K. Krzęcio-Nieczyporuk Elżbieta Adamczyk G & Koćwin-Podsiadła M. relationship between glycolytic potential and meat quality of duroc pigs with consideration of carcass chilling system. Ann Anim Sci. 2013;13:645–54.
Bongiorni S, Gruber CEM, Chillemi G, Bueno S, Failla S & Moioli B, et al. Skeletal muscle transcriptional profiles in two italian beef breeds,chianinaandmaremmana, reveal breed specific variation. Mol Biol Rep 2016;43(4):253–268.
Busu L, Kai S, Jie M, Li L, Zhang GF. Integrated application of transcriptomics and metabolomics provides insights into glycogen content regulation in the pacific oyster crassostrea gigas. BMC Genomics. 2017;18(1):713.
Law J, Jovel J, Patterson J, Ford G, O’Keefe S, Wang W, et al. Identification of hepatotropic viruses from plasma using deep sequencing: a next generation diagnostic tool. PLoS One. 2013;8:e60595.
Luo R, Liu B, Xie Y, Li Z, Huang W, Yuan J, et al. Soapdenovo2: an empirically improved memory-efficient short-read de novo assembler. Giga Science. 2012;1:18.
Qin N, Yang F, Li A, Prifti E, Chen Y, Shao L, et al. Alterations of the human gut microbiome in liver cirrhosis. Nature. 2014;513:59–64.
Karlsson FH, Tremaroli V, Nookaew I, Bergström G, Behre CJ, Fagerberg B, et al. Gut metagenome in european women with normal, impaired and diabetic glucose control. Nature. 2013;498:99–103.
Sunagawa S, Coelho LP, Chaffron S, Kultima JR, Labadie K, Salazar G, et al. Structure and function of the global ocean microbiome. Science. 2015;348:1261359.
Buchfink B, Xie C, Huson DH. Fast and sensitive protein alignment using diamond. Nat Methods. 2014;12:59–60.
Oh J, Byrd AL, Deming C, Conlan S. NISC comparative sequencing program, Kong HH. Segre JA. Biogeography and individuality shape function in the human skin metagenome. Nature. 2014;514:59–64.
Phillippy AM, Bergman NH, Ondov BD. Interactive metagenomic visualization in a web browser. BMC Bioinformatics. 2011;12:385.
Avershina E, Frisli T, Rudi K. De novo semi-alignment of 16s rRNA gene sequences for deep phylogenetic characterization of next generation sequencing data. Microbes Environ. 2013;28:211–6.
Noval Rivas M, Burton OT, Wise P, Zhang YQ, Hobson SA, Garcia Lloret M, et al. A microbiota signature associated with experimental food allergy promotes allergic sensitization and anaphylaxis. J Allergy Clin Immun. 2013;131:201–12.
White JR, Nagarajan N, Pop M, Ouzounis CA. Statistical methods for detecting differentially abundant features in clinical metagenomic samples. PLoS Comput Biol. 2009;5:e1000352.
Segata N, Izard J, Waldron L, Gevers D, Miropolsky L, Garrett WS, Huttenhower C. Metagenomic biomarker discovery and explanation. Genome Biol. 2011;12:1–18.
This work was supported by grants from the Earmarked Fund for Modern Agro-industry Technology Research System (Grant No. CARS-36), Zhejiang Major Scientific and Technological Projects for Breeding of New Breed of Livestock and Poultry (Grant No. 2016C02054).
Animal studies were conducted in accordance with the guidelines of the Zhejiang Farm Animal Welfare Council of China and approved by the ethics committee of Zhejiang Academy of Agricultural Sciences.
Consent for publication
The authors declare that they have no competing interests.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Figure S1. Observed number of non-redundant genes of the 10 total samples. (TIFF 168 kb)
Figure S2. Relative abundances of the bacteria in the top 10 species level taxa of the DF and FF groups. (TIF 275 kb)
Figure S3. Functional classification of annotated NR genes by using KEGG analysis. (TIF 304 kb)
Table S1. Total statistic information of the metagenome sequencing by using HiSeq Illumina platform from 10 faecal samples. (XLS 57 kb)
Table S2. Absolute abundance of unigenes in species level. (XLS 1105 kb)
Table S3. Composition of the basal diets in the traditional feeding farm. Table S4. Composition of the basal diets the semi free-grazing farm. Table S5. Primer sequences (5′to 3′) used for the quantitative polymerase chain reaction. (DOC 59 kb)