Microbiome-Informed Food Safety and Quality: Longitudinal Consistency and Cross-Sectional Distinctiveness of Retail Chicken Breast Microbiomes

Chicken has recently overtaken beef as the most-consumed meat in the United States. The growing popularity of chicken is accompanied by frequent occurrences of foodborne pathogens and increasing concerns over antibiotic usage. Our study represents a proof-of-concept investigation into the possibility and practicality of leveraging microbiome-informed food safety and quality. Through a longitudinal and cross-sectional survey, we established the chicken microbiome as a robust and multifaceted food microbiology attribute that could provide a variety of safety and quality information and retain systematic signals characteristic of overall processing environments.

M icroorganisms and their communities on foods are important determinants and indicators of food safety and quality. There have been growing interests in studying food and food-related microbiomes. Deep sequencing of food and environment microbiomes allowed the detection of multiple foodborne pathogens and the monitoring of their relative abundances throughout the beef production chain (1). Microbiome analyses of meat and seafood products provided a microbiological dissection of food spoilage to study microbial origin and population dynamics of the complex process (2). Besides targeted investigations of pathogenic and spoilage microorganisms, overall food microbiome profiles have been investigated to suggest how agricultural and food processing practices affect resident food microbial communities in the context of food safety and spoilage (3)(4)(5). In addition to taxonomic compositions, particular microbiome constituents such as antimicrobial resistance genes (ARGs), or Nonpareil (62), the sequencing coverage of chicken breast microbiomes ranged from 22.5% to 93.8%. Two samples that had a coverage below 75% were considered insufficiently sequenced and excluded from subsequent analyses. The coverages of the remaining 35 samples ranged from 75.6% to 93.8%, with a mean of 83.7% ( Table 2). The average aerobic plate count (APC) of these samples ranged from 4.2 to 7.0 log 10 CFU/g, with a mean of 5.9 log 10 CFU/g ( Table 2). There was no significant difference in APC between vacuum-packing and air-packaging samples (P ϭ 0.90, t test). According to redundancy-based coverage projections by Nonpareil, ϳ5.6 million bacterial reads (ϳ1.1 billion bases) are required for 95% abundance-weighted average coverage of the chicken breast bacterial communities (see Fig. S1 in the supplemental material).
Longitudinal consistency and cross-sectional distinctiveness of chicken breast microbiomes. Principal coordinates and hierarchical clustering analyses based on overall microbiome similarities revealed consistent clustering of chicken breast microbiomes over the sampling period ranging from 6 weeks up to 3 months (Fig. 1). Four major clusters were identified (C1 through C4), of which C2 was most distinctive and exclusively composed of vacuum-packaged samples. Clustering of samples by processing establishments was observed in other clusters, of which C1 and C3 consisted of samples from only processing establishments I and IV, respectively. C4 was further divided into 2 subclusters, C4-1 and C4-2. C4-1 contained all the samples from processing establishment II, and C4-2 included all the air-permeable (nonvacuum) samples from establishment III and a single establishment II sample (Fig. 1b). C4-2 was the only cluster or subcluster that had samples under different brands (B, C, and D). It is worth noting that brand D samples in C2 (D-III-ABF-V) and C4 (D-III-ABF) were both antibioticfree (ABF) products processed at the same establishment, differing only by package types (vacuum versus air permeable). Two additional samples (D-III-ABF-V_A and D-III-ABF-V_V) were sampled 11 months later after the original D-III-ABV-V samples to confirm the packaging effect on the chicken breast microbiome as detailed below.
Vacuum-packaged samples had significantly higher alpha diversities than airpermeable samples (Shannon index, Fig. 2b; Simpson index, Fig. S2a). Group comparison of beta diversity (Bray-Curtis dissimilarity) also showed a significant discrimination between vacuum-packaged and air-permeable products (Table 3). To confirm the impact of vacuum packaging on the chicken breast microbiome, we compared the microbiomes of two additional vacuum-packaged samples by unwrapping one sample and exposing it to air (D-III-ABF-V_A) while keeping the other one's vacuum packaging intact (D-III-ABF-V_V) during storage. An apparent shift in microbiome composition was observed between the 2 samples, resulting in a Pseudomonas-dominated microbiome on the unwrapped sample ( Fig. 1 and 2a). The 2 additional samples were collected 11 months after the original D-III-ABF-V samples and had a different packaging barcode and design; their overall microbiome profiles diverged from the earlier D-III-ABF-V samples. Because of the significant impact of vacuum packaging on chicken breast microbiomes, vacuum-packaged samples were excluded from evaluating whether other production factors affected the microbiomes, including overall processing environments (different processing establishments), antibiotic usage (ABF versus CONV products), and seasonality (July 2017 to September 2017 versus November 2017 to January 2018). Consistent with overall microbiome distinctiveness associated with processing establishments according to Bray-Curtis dissimilarities (Fig. 1), significant compositional differences were observed among microbiomes from different establishments as measured by both alpha (Shannon index, Fig. 2b; Simpson index, Fig. S2a) and beta (Table 3) diversities. However, the numbers of detected genera were not significantly different among these microbiomes (Fig. S2b), suggesting that the observed microbiome differences were caused by genus diversity instead of genus richness.   showing compositional differences among chicken breast samples. "*" indicates two additional vacuum-packaged samples, of which D-III-ABF-V_A was unwrapped and exposed to air and D-III-ABF-V_V was kept intact during storage. (b) Comparisons of bacterial alpha diversities (Shannon index) based on genus profiles showing the differences between packaging types, antibiotic usage, processing establishments, and seasonality. Samples were colored by products. Circles represent CONV samples, and triangles represent ABF samples. January 2018 (https://www.wunderground.com/). According to FSIS, all poultry should be chilled immediately after processing and stored at 40°F (4.4°C) or less before packaging and shipping (19).
Chicken breast microbiomes were distinct from other poultry-related microbiomes. We compared chicken breast microbiomes with other poultry-related microbiomes including those of whole chicken carcasses, chicken litter, and chicken feces. Also included in this comparison were fecal microbiomes from nonpoultry animals and humans. Principal-coordinate analysis (PCoA) of the microbiomes based on Bray-Curtis dissimilarities identified 4 major clusters including chicken breasts, other poultryrelated sources, nonpoultry animal sources, and human feces (Fig. 3). Notably, the average pairwise Bray-Curtis dissimilarity between chicken breasts and other poultryrelated sources (0.12 Ϯ 0.03) was smaller than that between chicken breasts and nonpoultry animal sources (0.18 Ϯ 0.11) and that between chicken breasts and human feces (0.32 Ϯ 0.06).
Packaging type had greater impact on chicken breast resistomes than antibiotic usage. In total, 132 ARGs of 10 antibiotic classes were detected from chicken breast microbiomes (Table S2 and Fig. S3). The majority (53.8%) of the detected ARGs belonged to the classes of aminoglycosides (n ϭ 14) and beta-lactams (n ϭ 57) (Table S2).
Resistome comparison between vacuum-packaged and air-permeable samples showed significant differences in ARG abundance and composition (Fig. 5a). The  average normalized ARG abundance of the vacuum-packaged samples (n ϭ 5) was 4.5 times higher than that of air-permeable samples (n ϭ 28) ( Fig. 5a; Table S3). Significant differences in alpha diversity (Shannon index) (P ϭ 0.000, Mann-Whitney U test) of detected ARGs were found between vacuum-packaged and air-permeable samples. The average Shannon index was 3.33 Ϯ 0.27 for vacuum-packaged samples and 1.74 Ϯ 0.48 for air-permeably packaged samples (Table S4). The resistomes of vacuum-packaged samples were dominated by beta-lactam ARGs, with an average relative abundance of 69.6% Ϯ 7.4% (mean Ϯ SD). In comparison, aminoglycoside ARGs appeared to be more abundant in air-permeable samples (50.4% Ϯ 24.8%), accounting for the majority of the resistome in 21 out of 28 samples (61.3% Ϯ 16.5%) ( Fig. 5a and b).
Compared with packaging type, the impact of antibiotic usage on chicken breast resistomes, if any, was much milder. There was no significant difference in normalized ARG abundances (P ϭ 0.97, t test) between ABF samples (mean abundance ϭ 0.10) and CONV samples (mean abundance ϭ 0.10) with air-permeable packaging. Similarly, when comparing the diversity of ARGs, no significantly different Shannon indices (P ϭ 0.834, Mann-Whitney U test) were found between CONV and ABF samples (excluding vacuum-packaged samples). The average Shannon index was 1.76 Ϯ 0.58 for ABF samples and 1.72 Ϯ 0.42 for CONV samples (Table S4).
Sensitive and specific detection of Salmonella by shotgun metagenomic sequencing. Salmonella was detected by culture enrichment in 4 samples of 2 products (B-II-CONV_4, B-II-CONV_6, C-III-ABF_1, and C-III-ABF_3) ( Table 4). The Salmonellapositive rate among the 33 chicken breast samples was 12.1% (2 samples sequenced later for evaluating vacuum packaging effect were not included for Salmonella detection), which was similar to Salmonella prevalence on retail skinless chicken breasts (12.3%) that were surveyed spatio-temporally close to the samples used in the current study (Atlanta metropolitan area, Georgia, in 2017) (21). Salmonella detection in metagenomic DNA by real-time PCR yielded negative results for all samples (Table 4).
Salmonella detection was also carried out by identifying Salmonella sequences in shotgun sequencing data. According to MetaPhlAn2 (76) identification using reference marker genes, no Salmonella sequence was detected in any sample. Based on Salmonella-specific k-mers detected by KrakenUniq, 17 samples were positive for Salmonella sequences (Table 4). Notably, the metagenomes of vacuum-packaged sam-ples were particularly abundant in unique k-mers classified by KrakenUniq (60) as Salmonella (Table 4). After filtering out sequencing reads that were likely of plasmid or phage origin by PPR_Meta, 5 samples remained as Salmonella positive, 3 of which were Salmonella positive also by culture enrichment (Table 4).
Using Salmonella detection by culture enrichment as a benchmark, we evaluated the sensitivity and specificity of Salmonella detection from shotgun sequencing data ( Table 4). The best overall performance was delivered by reads classification followed by mobile genetic elements (MGE) filtering, with a sensitivity of 0.75 and specificity of 0.93. In comparison, reads classification by KrakenUniq alone delivered a sensitivity of 0.75 and specificity of 0.52 for Salmonella detection. Both real-time PCR and MetaPhlAn2 analyses of the microbiomes yielded a sensitivity of zero by detecting Salmonella in none of the samples.

DISCUSSION
Microbial populations on chicken experience drastic shifts in composition and abundance over multiple stages during production including live production, poultry processing, retail display, and storage. While microbiome characterization throughout the "farm-to-fork" continuum is possible and has been attempted (22), a routine practice of microbiome monitoring, if justified, would require a more targeted ap- Links were drawn when Spearman's correlation coefficient () between two nodes (ARG or microbial genus) was larger than 0.8 (P Ͻ 0.01). The width of each link is proportional to the value of the correlation coefficient. The size of each node is proportional to the relative abundance of ARG or microbial genus.
proach that is scientifically effective and logistically viable. Our study captured postprocessing and preconsumption snapshots of chicken breast microbiomes. Our results suggest that chicken microbiomes at this junction allow robust metagenomics characterization through shotgun sequencing, retain sufficient information for microbiological evaluation of poultry processing, and provide multiple food safety and quality characteristics such as the presence of Salmonella in the final products that directly matter to consumers.
Overall, we observed distinct and temporally consistent chicken breast microbiomes associated with individual processing establishments. The observation indicates that variation in processing practices may lead to quantifiable and stable characteristics in finished product microbiomes. We envision that microbiome monitoring of poultry products can help establish a baseline of the regular chicken microbiome, which may be characteristic of individual processing establishments and potentially useful in signaling irregularities or perturbations to normal processing routines, such as change or failure of antimicrobial agents or protocols. As we have shown that poultry micro-  biomes inform the presence of Salmonella and the composition of spoilage organisms among other potential metagenomics markers of food safety and quality, the accumulation of baseline and irregular microbiome data may help classify or rate poultry microbiomes according to these markers, which may eventually lead to the evaluation of processing practices and environments. Shotgun metagenomics sequencing promises unbiased and comprehensive characterization of food microbiomes unmatched by 16S rRNA amplicon sequencing. However, its efficiency is challenged by food matrices that have high ratios of food to microbial cells. In such cases, large proportions of sequencing output are food DNA sequences uninformative to food safety and quality. For instance, in a recent study, over 99% of shotgun sequencing data of beef microbiomes came from the bovine genome, leaving less than 1% of sequencing data attributed to bacteria (6). Chicken breast rinsates sampled in this study were much less contaminated by animal cells, resulting in an average bacterial yield being about half of the sequencing output. The relatively high efficiency of chicken microbiome sequencing allowed the coverage of the vast majority of microbial taxa in each microbiome sample by sequencing the equivalent of 5 or 6 samples in a MiSeq sequencing run (after combing reads from 2 replicate samples). These results suggest that sufficient characterization of chicken breast microbiome is possible through multiplexed shotgun sequencing, which can be relatively cost-effective and potentially practical for industrial applications.
In this study, we synchronized chicken breast microbiomes by holding samples under refrigeration until their designated "sell by" dates. During the holding storage, the microbiomes shifted from their earlier states off processing lines as indicated by increasing APC over time (chicken breast APC before "sell by" dates not shown). Despite the shift, overall microbiome profiles clustered by their association with particular processing establishments, and the clustering largely persisted throughout the sampling period with products up to 3 months apart. Such longitudinal consistency and cross-sectional distinctiveness suggest that chicken breast microbiomes on "sell by" dates carry over and retain sufficient metagenomic information that was reflective and characteristic of the processing environment at individual establishments.
Poultry processing converts live birds to finished products in multiple steps, including prescald brushing, scalding, defeathering, evisceration, rehang, on-and off-line reprocessing, carcass washing, chilling, and postchill treatment. Such steps can have substantial impact on the microbial ecology of processing environments and chicken meats. Some examples include microbial contamination due to gut rupture and gut content spillage during evisceration (23) and antimicrobial treatment during equipment spraying (24) and carcass washing (25). Complete microbial monitoring of all these steps can therefore be complex and logistically challenging for routine practices. Alternatively, after chicken meats go through the processing steps and pick up environmental influences along the way, the chicken meat microbiome constitutes an endpoint sample that represents the eventual consequences of processing environments on the finished products. It is possible that microbiome characterization of chicken meats right off processing provides closer indicators of microbial ecology of poultry processing. However, obtaining sufficient metagenomic DNA for shotgun sequencing from samples immediately after processing remains a challenge as we experienced in the current study. Such samples were much less abundant in microorganisms than samples held until "sell by" dates. According to our preliminary analysis of 3 brands of chicken breast, their APC ranged from 0.9 to 2.2 log 10 CFU/g when assayed on the same day of purchase. Similarly, APCs of beef trimmings in the United States are routinely between 2 and 3 log 10 CFU/g (26). The low microbial load on beef likely also contributed to the low output of microbial DNA in beef metagenomes (less than 1% compared with a mean of 49.9% in this study) as previously mentioned (6).
Interestingly, microbiomes of ABF samples under the same brand, from the same processing establishment, but collected 11 months apart diverged from each other (D-III-ABF-V_V versus D-III-ABF-V, Fig. 1). Different packaging barcodes and designs indicate different products. It is unknown whether any changes in the processing of the later product occurred. In contrast, samples of the same product collected in warm and cold seasons from a North Carolina processing establishment did not show any significant difference in microbial composition. This finding was inconsistent with the seasonal differences observed in 16S amplicon sequencing-characterized chicken microbiomes from South Korea (27). In the South Korean study, seasonal samples were collected from months of maximum and lowest rates of processing (July and January), while year-round broiler production in the United States does not feature noticeable seasonality (28). It is unclear whether the disparity in production volumes, any other structural difference in the production system, or actual climatic seasonality caused the seasonal differences found in South Korean chicken microbiomes.
Our study also demonstrated the significant influence of vacuum packaging on chicken breast microbiomes. Microbiomes of vacuum-packaged samples were more diverse and characteristic of a variety of facultative anaerobic bacteria, compared with those of air-permeable samples that were dominated by Pseudomonas. Similarly, Pseudomonas was found to become dominant on chicken cuts under modified atmosphere packaging (MAP) when there was a higher percentage of oxygen in the MAP (29) and on beef chops during refrigerated storage in air (30). While inhibited, Pseudomonas, an obligate aerobe, was still present on vacuum-packaged samples, which is consistent with a previous observation that Pseudomonas was still able to grow on foal meat under refrigerated vacuum storage (30) and suggests the presence of residual oxygen in vacuum packages. Pseudomonas and Shewanella, which were dominant on chicken breast samples with air-permeable packaging, are major psychrotrophic bacteria that cause meat spoilage under aerobic storage conditions (31). On the other hand, Aeromonas, Buttiauxella, Carnobacterium, Hafnia, and Serratia, which were found to be more abundant on vacuum-packaged samples, commonly occur during chilled storage of raw meat in MAP or vacuum packaging (32). Our data suggest the potential of food microbiome characterization in evaluating and developing packaging technologies for extending product shelf life against spoilage microorganisms.
Excessive use of antibiotics in food animal production is perceived to be an important contributor to the rise of antibiotic resistance throughout food production environments (33). Shotgun metagenomics has been suggested to help quantitative risk assessment of antibiotic resistance in the food supply chain (33)(34)(35). Our resistome comparison between CONV and ABF samples identified no significant difference in ARG abundance and composition between the two production practices. Furthermore, ARG abundance of the chicken breast samples was considerably lower than that of any other livestock samples and all but one (U.S. soil) environmental samples analyzed in the current study. These results indicate a low risk of ARG accumulation on chicken breasts regardless of antibiotic usage in live production. Our finding is consistent with antimicrobial resistance surveys in beef production. It was reported that no ARG was detected on beef products and that interventions during slaughtering and beef processing might reduce the risk of ARG transmission to consumers (6). Similar levels of antimicrobial resistance were measured between conventional and "Raised Without Antibiotics" ground beef products (36) and beef cattle (37), which led to the doubt on the assumed impact of antimicrobial use in U.S. beef production and whether significant reductions in antimicrobial resistance could be yielded by reducing antimicrobial use. Although a recent study identified differences in resistomes between farms with different practices regarding antibiotic use (conventional versus "Raised Without Antibiotics"), the authors did not conclude a cause-effect relationship between antimicrobial use and the resistome differences because other factors (e.g., location, cattle source, management practices, and diet) could have affected the resistomes (38). It should be noted that our findings regarding the resistomes of ABF versus CONV samples were limited to finished chicken products. It remains to be investigated whether resistome differences related to antibiotic usage exist in any production environments, and whether such differences, if any, can be largely neutralized during poultry processing, similar to cattle processing (6).
Between chicken breasts, the resistomes of vacuum-packaged products had higher normalized ARG abundances than the resistomes of air-permeable samples, largely due to the higher levels of beta-lactam-class ARGs in the vacuum-adapted resistomes (Fig. 5). The higher ARG abundances were unrelated to antibiotic use because all the vacuum-packaged products were also ABF. The perceived resistome differences between vacuum-packaged and air-permeable samples were caused by different microbiome compositions as the consequence of oxygen availability. For example, Aeromonas was abundant on vacuum-packaged samples and outcompeted on air-permeable samples, whereas Pseudomonas was dominant on air-permeable samples and inhibited on vacuum-packaged samples (Fig. 2a). Notably, Pseudomonas species have been reported to be antagonistic against Aeromonas (39). Network analysis revealed the cooccurrence of Aeromonas and 3 beta-lactam ARGs (cphA1, cphA5, and bla OXA-427 ) in chicken breast microbiomes, which was substantiated by the fact that cphA genes were commonly found on Aeromonas chromosomes (40) and confirmed by ARG analysis of metagenome contigs (see Table S5 in the supplemental material). Therefore, metagenomics profiling of food resistome needs to be interpreted with caution when it comes to evaluating product safety, quality, and the impact of antibiotic usage during production. ARGs naturally occur; their presence and even elevated abundance in food microbiomes do not necessarily indicate selection by antibiotics or occurrence of multidrug-resistant organisms. As we showed, normalized ARG abundances (ARG copies per prokaryotic cell) vary among different bacterial populations, with densely populated populations such as livestock soil and gut (fecal) microbiomes displaying higher abundances. Abundant ARGs on food samples may be associated with certain production environments naturally rich in ARGs, which does not directly translate into higher food safety or quality risks. Therefore, food ARG abundances should not be interpreted in isolation but considered along with their bacterial hosts and environmental origins. While source attribution of ARGs using short-read metagenomics data is challenging, cooccurrence identification by network analysis as shown in the current and previous studies (41) can help reveal prominent associations. Low abundances of mcr genes (variants of mcr-3, mcr-4, mcr-7, and mcr-9) were identified in 15 out of 33 retail chicken microbiomes. Notably, mcr-3-positive microbiomes (n ϭ 7) had significantly higher relative abundances of Aeromonas than mcr-3negative microbiomes (n ϭ 26), indicating that Aeromonas might be a source of mcr genes in chicken microbiomes. Indeed, chromosome-mediated mcr-3 variants have been found in Aeromonas isolates from chicken in China (42,43). Aeromonas has also been speculated as a reservoir of mcr-7 (44). Other bacterial hosts of mcr genes have also been identified in our chicken microbiomes, including Klebsiella (45) and Shewanella (46). Together, our results provide metagenomics evidence for the prevalence of mcr genes in poultry production (47), although we were not able to determine whether the mcr genes were plasmid-borne and/or carried by common foodborne pathogens such as Escherichia coli and Salmonella.
Finally, we showed that shotgun metagenomics data supported sensitive and specific detection of S. enterica from chicken breasts as benchmarked by the gold standard of culture enrichment. Shotgun metagenomics considerably outperformed real-time PCR in detecting S. enterica from the same DNA samples; the latter method failed to detect the pathogen in any of the 4 culture-confirmed chicken breast samples. These results suggest that taxonomic classification of sequencing reads increased detection sensitivity by allowing comprehensive screening for S. enterica signals that may come from multiple loci of the pathogen's genome. In comparison, real-time PCR targeted only a single locus and was less sensitive in detecting S. enterica when only fractional S. enterica genomes were present in DNA samples extracted from chicken breast rinsates.
On the other hand, the specificity of taxonomic classification was compromised by mobile genetic elements such as plasmids and phages. Using k-mer-based taxonomic classification alone, 17 out of the 33 samples were called S. enterica positive. A similar issue was observed in metagenomics detection of S. enterica from cattle feces and attributed to the presence of plasmid sequences in the reference database (48). In our study, we found that exclusion of reads of likely plasmid and phage origins by PPR_Meta substantially increased detection specificity (from 0.52 to 0.93, Table 4). This method was particularly effective in correcting false-positive calls on vacuum-packaged samples by sole taxonomic classification, which yielded large quantities of S. entericalike reads and unique k-mers for these samples (Table 4). It is possible that vacuumpackaged samples, whose microbiome compositions were distinct from those of air-permeable samples, contained non-Salmonella bacteria that shared plasmids or phages with S. enterica.
Sequencing coverage of metagenome can affect the sensitivity of Salmonella detection. In our study, except for 2 outliers that were excluded from analysis due to exceedingly low metagenomics coverage (22.5% and 55.5%), the estimated coverage of samples ranged from 75.6% to 93.8% (mean ϭ 83.7%) ( Table 2). These samples yielded a sensitivity of 75% in detecting S. enterica using the metagenomics approach. The results suggest that an estimated metagenomics coverage above 80% (ϳ5.3 million reads) has the potential to detect S. enterica from the vast majority of contaminated chicken breast microbiomes. Higher coverages may improve detection sensitivity, but it cannot be statistically determined in our study due to the limited sample size (n ϭ 35), and especially the small number of Salmonella-positive samples (n ϭ 4). In fact, the only false-negative sample by metagenomics detection had a higher coverage (87.6% for CIII-ABF_3) than 2 of the 3 true-positive samples (77.6% for B-II-CONV_4 and 79.8% for B-II-CONV_4), suggesting stochastic coverage of a low-abundance organism in the microbiome. Future studies of larger scale may help determine optimal levels of metagenomics coverage that balance the sensitivity of pathogen detection and the demand for sequencing output.
The performance of metagenomics detection of Salmonella was benchmarked by culture enrichment that is commonly regarded as the gold standard for foodborne pathogen detection. However, the equivalency between metagenomics and culture detection may be skewed by the presence of Salmonella cells that are difficult or unable to culture. The so-called viable but nonculturable (VBNC) state of Salmonella cells can be induced by adverse environmental conditions such as chlorine treatment (49) and starvation (50). Dead pathogen cells whose DNA molecules are still present in food samples are considered to be a source of false-positive results by DNA-based detection methods incapable of distinguishing between signals from live and dead cells (51). In our study, Salmonella was detected by metagenomics sequencing in 2 samples that were Salmonella negative by culture enrichment, which could be attributed to difficultto-culture or dead Salmonella cells. Further validation of metagenomics detection of foodborne pathogens should evaluate the likelihood that DNA from dead cells could be sequenced from food samples and whether removal of extracellular DNA during sample preparation could mitigate or avoid sequencing DNA from dead cells. While a Salmonella-positive chicken microbiome does not definitively suggest the presence of live and infectious Salmonella cells on the product, it indicates the occurrence of Salmonella in the processing environment. Recent advances in quantitative metagenomics analysis, such as the use of synthetic internal DNA standards (52), may allow robust qualification of pathogen or indicator DNA in food microbiomes, which could help develop probabilistic or quantitative risk assessment of microbial contamination in processing establishments.
Metagenomics detection of Salmonella has also been performed by 16S rRNA amplicon sequencing using universal primers for 16S rRNA genes on retail chicken products from France (29) and South Korea (8). Neither study detected Salmonella; it is not known whether the negative results reflected low Salmonella prevalence on retail chicken in these countries or a difference in Salmonella detection sensitivity between 16S rRNA amplicon and shotgun sequencing. A recent comparison between the two metagenomics approaches in detecting foodborne bacteria suggested superior performance by shotgun sequencing (53). Using a mock bacterial community, 16S rRNA amplicon sequencing based on universal primers showed a lack of consistency and missed several pathogenic species (53). Nonetheless, it is possible to use specific primers for targeted amplicon sequencing to improve the detection of particular pathogens. 16S rRNA amplicon sequencing is relatively cost-effective and able to characterize communities with low DNA levels, which can have advantageous applications.
Conclusion. Our study represents a proof-of-concept investigation into the possibility and practicality of leveraging microbiome-informed food safety and quality. Through a longitudinal and cross-sectional survey of a household item, we established the chicken microbiome as a robust and multifaceted food microbiology attribute that could provide a variety of safety and quality information.
The methodology described here can apply to studying other food commodities. The diversity of foods and their microbiomes warrants expanding the investigation to more food matrices to better explore and evaluate metagenomics applications in food safety and quality.

MATERIALS AND METHODS
Sample collection and preparation. Eight different retail boneless and skinless chicken breast products under 5 brands (A, B, C, D, and E) were purchased from grocery stores in Griffin (Atlanta metropolitan area), Georgia. For 7 products, a total of 35 samples were collected between July 2017 and January 2018; these chicken samples were collected every other week for 3 to 6 months depending on their retail availability. Two more samples of an additional product were later collected in June 2018 to confirm the effect of vacuum packaging on the product microbiome. All 37 samples of the 8 products were traced back to 4 processing establishments (I, II, III, and IV) in Georgia, Arkansas, and North Carolina using the FSIS establishment numbers (54) on the packages. Five of the 8 products were considered antibiotic free (ABF) according to the "No Antibiotics Ever" package label. The other 3 without any antibiotic usage claim were classified as conventional products (CONV). Two products of 8 ABF samples under brand D and from processing establishment III were vacuum packaged (D-III-ABF-V and D-III-ABV-V_V/A); these products were sealed in thick oxygen-impermeable film. Other samples had air-permeable packaging. Detailed sample information is summarized in Table 1. For each sample, two packages of the same product were purchased and prepared as biological replicates. With the exception of sequencing coverage estimation, for any quantitative metagenomics measurement detailed below, each replicate was individually analyzed and the average of two measurements was taken. For the coverage estimation and the k-mer-based clustering analyses of chicken microbiomes, sequencing reads from the two replicates of each sample were pooled and analyzed.
Breast samples were stored at 4°C until their designated "sell by" dates. On the "sell by" date, each sample was inspected; no visual sign of spoilage and off-odor was detected. Each aliquot of 500 g meat was rinsed in 150 ml buffered peptone water (BPW; Acumedia, USA) by hand massage for 2 min. Thirty milliliters of the rinsate was collected for aerobic plate count and culture-based detection of Salmonella.
One hundred milliliters of the rinsate was pelleted by a two-step centrifugation before DNA extraction. The rinsate was first centrifuged at 100 ϫ g for 10 min to remove heavier food particles. The supernatant was transferred to another tube and centrifuged at 5,000 ϫ g for 15 min to collect bacterial cells. Each cell pellet was washed by resuspension in 10 ml of 1ϫ phosphate-buffered saline (PBS; Fisher Scientific, USA) followed by another round of centrifugation at 5,000 ϫ g for 15 min. Finally, the pellet was resuspended in 5 ml of 1ϫ PBS and used for DNA extraction.
Aerobic plate count. Aerobic plate count was determined for each rinsate sample. Serial dilutions (10 Ϫ1 to 10 Ϫ4 ) of the rinsates were prepared in BPW, and enumerations were performed on plate count agar (PCA; Becton, Dickinson, USA) according to the Bacteriological Analytical Manual (BAM) (55). Plate count agars were incubated at 35°C for 48 h, and single colonies were enumerated from the appropriate dilutions.
Culture-based detection of Salmonella. Culture-based detection of Salmonella was conducted according to the Microbiology Laboratory Guidebook (USDA-FSIS) (56) with modifications. Briefly, each 30-ml aliquot of BPW chicken rinsate was incubated at 37°C for 24 h. An 0.1-ml amount of each preenriched rinsate was transferred to 10 ml of Rappaport-Vassiliadis soya (RV) broth (Acumedia, USA) and incubated at 42°C for 24 h. After the selective enrichment in RV, a loopful of the enriched sample was streaked on xylose-lysine-Tergitol 4 agar (XLT4 agar; Becton, Dickinson, USA) and incubated at 37°C for 24 h. The presumptive Salmonella colonies from the XLT4 agar were plated on triple sugar iron agar (TSI; Becton, Dickinson, USA) and lysine iron agar (LIA; Acumedia, USA) and incubated at 37°C for 24 h, followed by confirmation using real-time PCR. Real-time PCR for the detection of Salmonella was conducted as previously described (57). Samples with a threshold cycle (C T ) value of Ն40 or higher than that of negative control were considered negative results.
DNA extraction and real-time PCR detection of Salmonella in DNA samples. Pellets of chicken rinsate samples were subjected to DNA extraction using the PowerSoil extraction kit (Mo Bio, USA) according to the manufacturer's instructions. DNA concentrations were determined using a Qubit BR double-stranded DNA (dsDNA) assay kit (Invitrogen, USA) and were diluted to 0.2 ng/l for library preparation. Real-time PCR for Salmonella detection in DNA samples was performed as described above.
Shotgun metagenomic sequencing. DNA libraries were prepared according to the Illumina Nextera XT DNA library prep kit reference guide. Shotgun metagenomic sequencing was performed on an Illumina MiSeq platform with the paired-end sequencing strategy (2 ϫ 250 bp). Initial quality inspection of the sequencing data was conducted using FastQC (58). Low-quality reads were removed or trimmed by Trimmomatic (59). The leading three and the trailing three nucleotides were removed from the reads, and a 4-nucleotide sliding window was used to remove nucleotides from the 3= ends when the average Phred score dropped below 20. After trimming, reads that were shorter than 75 bp were discarded.
Taxonomic classification of sequencing reads. Trimmed and filtered reads were taxonomically classified using KrakenUniq (60) with the NCBI nucleotide collection (nt) database (as of August 2018), which yielded the highest average recall according to KrakenUniq's performance test (60). Reads classified as plasmid or phage sequences were discarded after taxonomic classification. The relative abundances of microbial genus and species in sequencing samples were then determined using Bracken (61).
Estimation of microbial diversity coverage by shotgun metagenomics sequencing. Reads classified as "Bacteria" by KrakenUniq were used for further microbiome analyses. Nonpareil (v3.0) (62) was used to estimate average sequencing coverage of bacterial diversity in each microbiome sample.
Similarity and clustering analyses of chicken breast microbiomes. Pairwise Bray-Curtis dissimilarity was calculated between shotgun metagenomics samples using Simka with default settings (63). Simka computed the dissimilarities based on k-mer counts from sequencing reads classified as bacteria. Hierarchical clustering and principal-coordinate analysis (PCoA) were performed using Bray-Curtis dissimilarities with the base R package (64).
Alpha and beta diversities of chicken breast microbiomes. Genus-level read counts obtained from analyses by KrakenUniq and Bracken were used for alpha and beta diversity calculations. Alpha diversities were calculated using the Shannon index in the base R package and the Simpson index in the Vegan R package (65). Beta diversities were calculated using the Bray-Curtis dissimilarities in the Vegan R package. Genus richness was also calculated using the Vegan R package.
ARG identification from chicken breast microbiomes. Trimmed and filtered sequencing reads were screened for ARGs using the ARGs-OAP v2.0 (66) pipeline with default settings except using a customized reference database. ARGs were identified by aligning reads against the ResFinder database (67) (as of June 2019) using BLASTϩ blastn (68) with E value Ͻ1 ϫ 10 Ϫ7 . A read was annotated as an ARG sequence if its best hit had Ͼ80% identity and the alignment length was Ͼ75 bp (66). ARG sequence abundances were normalized against a curated list of single-copy genes known as the universal essential single-copy marker genes available in the ARGs-OAP pipeline (66). The abundance of ARG was expressed as copies of ARG per prokaryote cell.
ARGs were also identified from assembled metagenome sequences. Metagenomes were de novo assembled using classified bacterial reads (including plasmid and phage reads) with MEGAHIT (69). Only contigs over 1 kb were kept for analysis and merged at 99% identity with CD-HIT-EST (70). The resulting contigs were assembled again using Minimus2 from the AMOS package (71). A median of 6,719 contigs was obtained for each sample, and the median N 50 for these contigs was 2,751 bp (data not shown). These contigs were screened for ARGs against ResFinder (67) with ABRicate (72). Only hits with sequence identity Ͼ80% and query coverage Ͼ80% were kept. To determine the occurrence of ARGs in chromosomes or plasmids, we used BLASTϩ blastn to align the metagenome contigs against the NCBI nt database. Hits with sequence identity Ͼ70% and query coverage Ͼ70% were kept.
Statistical and network analyses. Statistical comparisons of microbiome diversities were performed using Mann-Whitney U test (2 levels) or Kruskal-Wallis test (Ͼ2 levels) in the base R package. Permutational analysis of variance (PERMANOVA) was performed to compare community difference among different categories (e.g., ABF versus CONF) using the "adonis" function in the Vegan R package.
For network analysis of microbial and ARG abundances, we first screened 48 ARGs that occurred in at least 5 retail chicken samples, and 65 microbial genera that occurred in at least 10 retail chicken samples. This filtering step was performed to remove the poorly represented ARGs and microbial genera with limited occurrence across samples and thus reduce artificial association bias (41). We then constructed a correlation matrix by calculating all possible pairwise Spearman's rank correlations among ARGs and microbial genera. Only correlations with the Spearman correlation coefficient () Ͼ0.8 and the P value Ͻ0.01 were considered statistically robust (41). The P values were adjusted using the Benjamini-Hochberg method (73). The robust pairwise correlations were then used to build a cooccurrence network between ARGs and microbial genera. The correlation analysis was performed using the psych R package (74), and the cooccurrence network was visualized by Gephi (75).
Salmonella detection from chicken breast microbiomes. We compared different approaches for Salmonella detection from the metagenomic shotgun sequencing data. (i) Taxonomic classification of the sequencing reads was performed by MetaPhlAn2 with its default marker database (76). Only samples with an assignment of "s__Salmonella_enterica" were considered positive results. (ii) Taxonomic classification was performed by KrakenUniq with the NCBI nucleotide collection (nt) database. Only samples with unique k-mer numbers Ն1,000 for Salmonella enterica were considered positive results (60). (iii) To reduce the false-positive Salmonella detection caused by Salmonella-like mobile genetic elements, such as plasmids and phages, that were present in non-Salmonella organisms, reads classified as "Salmonella enterica" by KrakenUniq were analyzed by PPR_Meta with default settings (77) to determine whether these sequences came from chromosomes, plasmids, or phages. Only samples with at least one chromosomal read classified as S. enterica were considered Salmonella positive.
Availability of data and materials. The sequencing data from this study are available in the NCBI Sequence Read Archive (SRA) under BioProject accession number PRJNA612140. Accession numbers for the publicly available shotgun metagenomics data used in the study are listed in Table S6.

SUPPLEMENTAL MATERIAL
Supplemental material is available online only.