The yellow-feathered broilers were randomly collected from the packaging site of a commercial slaughtering-line in a large-scale plant in China. The carcasses that had not been packaged and stored (0 day) were regarded as control samples. The carcasses were packaged with penetrated-air packaging (PAP, air-filling) or modified-atmospheres packaging (80% CO₂/20% N₂) and then were stored at 4°C for 4 and 8 days, respectively. MAP was applied with a polyethylene/polyamide/low-density polyethylene film, with O₂ permeability of 24 cm³/(m²⋅day⋅atm), CO₂ permeability of 78 cm³/(m²⋅day⋅atm), and water vapor permeability of 44 g/(m²⋅day). PAP was applied with a high-density polyethylene film, with O₂ permeability of 14483 cm³/(m²⋅day⋅atm), CO₂ permeability of 63683 cm³/(m²⋅day⋅atm), and water vapor permeability of 54 g/(m² day⋅atm). The whole research proposal of this study and more detailed information is shown in Figure 1.

The total microorganisms from broiler carcasses were collected by carcass-washing methods: a prepared half carcass (about 700 g weight) was transferred into a stomacher blender bag containing 300 mL 0.85% NaCl and 0.1% peptone water, and then the sample was shaken for 10 min with 100 rpm. Total viable counts of each sample were determined with plate count agar (PCA, Beijing Land Bridge Technology, Co., Ltd., Beijing, China) according to the China National Food Safety Standard methods (GB 4789.2-2010). Lactic acid bacteria were determined with MRS agar (Oxoid) according to the method proposed by Chouliara et al. (2007). Four independent carcasses were replicated for the numeration of total viable counts and lactic acid bacteria.

The homogenized liquid of each carcass was centrifuged with 1500 ×g (10 min, 4°C) to remove various meat tissues, and the supernatant was further filtered with a 0.22 μm pore diameter filter as described by Takahashi et al. (2017). Subsequently, the filter was inserted into a tube containing several beads, and was vortexed and washed three times with 1 mL of 0.85% NaCl, and then microorganisms were gathered by centrifugation with 10,000 ×g for 5 min. Total DNA of microorganisms was directly extracted using DNA Microbiome Kit (Qiagen) following the manufacturer’s recommendations. Three replications were performed for DNA extraction and metagenome sequencing.

The integrity, purity and concentration of total DNA was analyzed using 1% agarose gel electrophoresis, NanoPhotometer spectrophotometer (IMPLEN, Westlake Village, CA, United States) and Qubit dsDNA HS Assay Kit in Qubit 2.0 Fluorometer (Life Technologies, Carlsbad, CA, United States). The total DNA of each sample that passed quality control was used for library construction. A total amount of 1 μg DNA per sample was used as input material for the DNA sample preparations. Sequencing libraries were generated using NEBNext Ultra^TM DNA Library Prep Kit for Illumina (NEB, United States) following the manufacturer’s recommendations and index codes were added to attribute sequences to each sample. Briefly, the DNA sample (n = 3) was fragmented by sonication to a size of 300 bp, and then DNA fragments were end-polished, A-tailed, and ligated with the full-length adaptor for Illumina sequencing with further PCR amplification. At last, the PCR products were purified (AMPure XP system), and the concentrations of amplification were primarily quantified by Qubit2.0, and then the concentrations were diluted to 2 ng/uL. All libraries were analyzed for insert size distribution by Agilent 2100 Bioanalyzer, and finally, the concentration of library (>3 nM) was accurately quantified using real-time PCR to ensure the quality of libraries. The clustering of the index-coded samples was performed on a cBot cluster generation system according to the manufacturer’s instructions. After cluster generation, the library preparations were sequenced on an Illumina HiSeq 4000 platform in Novogene Bioinformatics Technology, Co., Ltd. (Beijing), and the 150 bp paired-end reads were generated.

The raw reads were filtered by removing low-quality sequences to obtain clean reads according to the following criteria: (1) removing reads containing N base more than 10 bp; (2) removing reads overlapping with adapter more than 15 bp; (3) removing reads with low-quality (quality values ≤ 38) more than 40 bp; and (4) removing reads probably from host by blasting the genome using SOAPAligner software. A total of 56,524.61 Mbp clean data were obtained from Illumina HiSeq system, with an average of 6,280.51 Mbp (Supplementary Table S1), indicating that a deep coverage was achieved, and the detail information of sequencing depth of metagenomics in each sample was shown in Supplementary Figure S1. All the downstream analyses were performed using high-quality clean data. Metagenomic sequence data is available at SRA site from NCBI with SAMN08123132 of biosample identifier (BioProject accession: PRJNA420874; SRA: SRS2729591).

The single-piece and mixture group assembly analysis was performed using SOAPdenovo software, and total of 820,974,610 bp scaffolds and 775,303,201 bp scaftigs were, respectively, obtained (Supplementary Table S1). Based on these assembly scaftigs (> =500 bp), total of 1,171,588 open reading frames (ORFs), with an average of 117,159 ORFs in each sample (Supplementary Table S2), were predicted using the MetaGeneMark software, and following by removing the redundancy of ORFs using CD-HIT software, total of 403,555 gene catalogs (unigenes) were obtained using SOAPAligner software (Supplementary Table S2). The abundances of each unigenes in samples were calculated based on the numbers of reads blasted to unigenes and length of unigenes (Cotillard et al., 2013; Oh et al., 2014). Then each unigene was compared with the NR database in NCBI, and total of 88.48% (357,061) of unigenes could be identified with the criteria of E value less than 10⁵, and after filtering the identified unigenes with a standard of minimum E value multiply ten (Qin et al., 2010), the annotation of unigenes was carried out using LCA algorithm analysis with MEGAN version, and finally, the taxonomic profiles were obtained based on the abundances and annotations, with 92.57% on family level and 84.00% on genus level, the detail information was shown in Supplementary Table S2.

A raw-broiler meat juice agar (RJA) plate prepared by mixing raw chicken juice (the fresh broiler breasts were homogenized with deionized water and then filtered through two layers of gauze) and 1% agar, was used to identify the dominant isolates in PAP and MAP samples. Both ingredients were maintained at 45°C before mixing. RJA plates were sterilized by irradiation at a dose of 6 KGy via the ⁶⁰Co source at Hangyu (Hangyu Irradiation Technology, Co., Ltd., Nanjing, China).

Each 25 g sample of surface meat from spoilage carcasses of PAP and MAP was aseptically weighed and homogenized in 75 mL sterilized 0.85% NaCl solution, and then the solution was streaked on an RJA plate and incubated aerobically at 25°C for 48 h. Each colony was obtained from the plates and then streaked on an RJA plate twice to obtain a single isolate. Total of 23 independent colonies were obtained from PAP samples, and 21 colonies were obtained from MAP samples. The obtained isolates were identified by 16S rRNA gene sequencing. Amplification of the 16S rRNA gene was performed using region-specific primers 27F (5′-aga gtt tga tcc tgg ctc ag-3′) and 1492R (5′-ggt tac ctt gtt acg act t-3′). The amplification products were purified and sequenced by Invitrogen (Invitrogen Biotechnology, Co., Ltd., Shanghai, China). The identified isolates were further confirmed using a VITEK 2 automated system (BioMerieux, France).

Four grams meat showing typical spoilage features and a 4 mL supernatant of each isolate identified by a combination of 16S DNA sequencing and VITEK 2 (Pseudomonas fragi and Pseudomonas fluorescens from PAP meat, Aeromonas salmonicida and Serratia liquefaciens from MAP meat) incubated for 48 h at 25°C was tested for volatile organic compounds by GC–MS. The approach of supernatant collection was performed as previous reported (Beristain-Bauza et al., 2016). Each sample was performed on an Agilent 7890A gas chromatograph coupled to an Agilent 5973C mass spectrometer. A slight modification of the method described by Jaffres et al. (2011) was used. Briefly, a 4 g or 4 mL portion was transferred into a 20-mL vial with a polypropylene screw-on cap and a PTFE/silicone septum to make it airtight. The vial was heated at 40°C for 40 min to equilibrate the system. The SPME fiber, an 85-μm Carboxen/polydimethylsiloxane StableFlex^TM (Supelco), was inserted through the septum and exposed in the headspace of the vial for 30 min to allow absorption of the volatile compounds into the fiber. The SPME was then introduced into the injector port of the gas chromatograph for 5 min in splitless mode, set at 250°C, to desorb the volatile compounds. The desorbed components were analyzed on an Agilent DB-WAX (30 m × 0.25 mm × 0.25 μm) capillary column. Helium was used as a carrier gas at a constant flow rate of 0.8 mL/min, and the oven temperature was programmed as follows: 40°C for 5 min, ramped at 3°C/min until reaching 140°C, ramped at 10°C/min until reaching 250°C, and then held for 5 min. The mass spectrometer was operated in the electron impact mode with the electron energy set at 70 eV and a scan range of 50–500 m/z. The NIST library and comparison with spectra and retention times of the standards were used to identify the volatile components.

All identified isolates were cultured for 24 h at 25°C. Cell density was adjusted to an optical density of 600 nm of 0.4 value. Two microliter aliquots of isolates were spotted onto an RJA plate. Proteolytic decomposition was measured after 3 days of incubation at 25°C. Four independent replications were performed for each isolate.

The results of numeration for total viable counts and lactic acid bacteria and the decomposition zone diameters were expressed as the mean ± standard deviation. Statistical significance was determined by a one-way Duncan’s ANOVA procedure using SPSS 13.0. The level of statistical significance was p < 0.05.

Principal component analysis (PCA) was used to dissociate the patterns of volatile organic compounds between each packaging pattern. Additionally, PCA was used to determine the relationship between the volatile organic compounds of spoilage meat in situ and the volatile organic compounds of isolates incubated in standard TSB. The PCA analysis was conducted by SIMCA-P software (Umetri AB, Sweden).

As shown in Figure 2, The initial number of bacteria was approximately 4.0 log CFU/g, indicating good hygienic operation of the slaughter plant, which was in line with the production requirements of chilled broiler carcasses. The numbers of total plate count and lactic acid bacteria on carcasses exceeded 7.0 log CFU/g after storage for 4 and 8 days for PAP and MAP, respectively, which was widely considered as the threshold of microorganism counts during meat spoilage. No difference in the numbers of both kinds of bacteria between PAP and MAP was found at the end of shelf life. Compared to PAP, MAP prolonged the shelf life by 4 days.

The sequence statistics and data summary are shown in Table 1. The sufficient amount and high quality of clean data could guarantee the accuracy of metagenomics. The clean data were used to determine diversity analysis. The heat map shown in Figure 3 revealed the dynamics of the 35 most abundant bacteria along with the packaging patterns. The packaging patterns greatly influenced the bacterial diversity of broiler carcasses at the family level and genus level. The dominant bacteria genera in control samples mainly included Streptococcus, Enterobacter, Empedobacter, Macrococcus, Enhydrobacter, and Aeromonas. Compared with control treatment, PAP caused a shift in microorganism diversity, with an increase in abundance of Pseudomonas, Arthrobacter, and Janthinobacterium and a decrease in Aeromonas, Citrobacter, and Klebsiella. In contrast, MAP dramatically increased the relative abundance of Carnobacterium, Vagococcus, Enterococcus, Lactobacillus and Serratia, as well as Hafnia. Shewanella and Weissella, which are commonly associated with food spoilage, showed lower abundance in both control and PAP treatments. Unexpectedly, high abundances of several food pathogens (Figure 3), such as Salmonella, Shigella, Vibrio and Yersinia, were observed in all groups, especially for control samples, indicating that more effective interventions should be applied to control pathogens in the slaughter plant.

The clustering tree based on the Bray–Curtis distance obtained from the relative abundances in family and species levels is shown in Figure 4. The results indicated that a high similarity in relative abundance was found between each duplicate in individual treatments, and low similarity was found between each treatment, directly suggesting that the repeatability of microorganism diversity within the same treatment was fairly good. Additionally, the dynamic of bacterial diversity was dramatically influenced, in particular at the species levels, which showed that MAP significantly increased the diversity of bacteria and altered the dominant bacterial species. The results mentioned above indicated that the bacterial community shift was greatly influenced by packaging patterns, and this finding was consistent with the results reported in Figure 5, which come from the PCA based on the bacterial diversity at the family and species levels.

Additionally, the structure of the top 10 dominant bacteria at genus and species levels was dramatically affected by the packaging patterns (Figure 6). Compared to control treatment, the abundances of Escherichia coli and Aeromonas veronii were dramatically reduced in MAP and PAP treatments, whereas the abundances of Lactococcus raffinolactis and Pseudomonas sp. (including fragi, CF149, and fluorescens) were sharply increased in MAP and PAP, respectively, especially for Pseudomonas sp., which occupied a dominant abundance in PAP at the end of shelf life. In contrast, Serratia liquefaciens, L. raffinolactis, Aeromonas salmonicida, and Carnobacterium maltaromaticum together become the dominant bacteria in MAP at the end of shelf life.

The plate streaking methods (RJA plates) was used to confirm the species of dominant spoilage bacteria. The bacteria of P. fragi (eight isolates), P. fluorescens (nine isolates) and C. shigense (six isolates), were identified in PAP samples, and the isolates of A. salmonicida (eight isolates), A. hydrophila (five isolates), A. media (three isolates), and S. liquefaciens (five isolates) were identified in MAP samples.

To explore the metabolic profiles of microorganisms subjected to packaging patterns, the volatile organic compounds of broiler meat under corresponding packaging at the end of shelf life were investigated. All compounds were further analyzed by PCA. Packaging patterns significantly altered the volatile organic compounds of meat (Figure 7A), which were mainly associated with the metabolic profiles of spoilage bacteria. To explore the relationship between the metabolic profiles of dominant spoilage isolates and bacteria community subjected to packaging, the PCA based on the volatile organic compounds of spoilage meat in situ and the supernatant of identified isolates incubated by GC–MS was determined (Figure 7B). Great differences in volatile organic compounds between A. salmonicida isolates and the other three isolates were observed. No difference in volatile organic compounds was found among P. fragi, P. fluorescens, and S. liquefaciens. Unfortunately, no relationship was determined between the volatile organic compounds of meat in situ and the volatile organic compounds from the supernatants of A. salmonicida, P. fragi, P. fluorescens, and S. liquefaciens. This finding suggests that the metabolic profiles of spoilage bacteria grown in standard culture could not represent the spoilage potential of dominant bacteria grown in the meat matrix.

To evaluate the potential spoilage ability of isolates, the identified isolates were individually assessed with decomposition zone diameters. Figure 8 indicates that there were great differences in proteolytic activity among these tested isolates, with Pseudomonas fragi showing the largest proteolytic activity, reaching more than 18 mm in decomposition zone diameter. On the other hand, S. liquefaciens had hardly any proteolytic activity in broiler meat. Other isolates showed comparatively strong proteolytic activity, with decomposition zone diameters ranging from 13.5 to 17 mm.