This study was conducted as a longitudinal observational experiment designed to determine the relationship between microbial community structure and beef spoilage. Ground beef production lots were used as the experimental subject for this experiment, as they facilitate repeated destructive sampling of a consistent product over the beef shelf life. Ground beef samples from 3 grinding lots (n = 3) were partitioned into subsamples for daily destructive sample collection. One subsample from each sample (lot) was collected daily for 14 d, including day 0 as the day of placement. These samples were used to analyze color, temperature, weight, nutrient composition, and aerobic plate counts (APC). After initial sampling, the remaining product was frozen at −18°C (FKFH21F7HW Frigidaire, Charlotte, NC, USA) for DNA extraction and thiobarbituric acid-reactive substances (TBARS) analysis to be conducted later in the experimental period.

Boxed beef shoulder clods (US Department of Agriculture institutional meat purchasing specification #114) were purchased from a commercial beef processing facility in the high plains region of the United States and commercially transported to the Auburn University Lambert-Powell Meat Laboratory under refrigerated conditions (1.5°C). The beef was ground using a 9.525-mm plate, followed by another grind with a 3.175-mm plate using a commercial meat grinder (Model AFMG-48, Biro Manufacturing Company, Marblehead, OH, USA). The grinder was disassembled, cleaned, and sanitized between grinding lots to simulate a different operation to produce 3 lots of ground beef. From each lot, subsamples were weighed out to 225 g, tray overwrapped, and randomly assigned sampling pull dates (N = 45). The samples were then placed into simulated retail storage in an LED-lit (2297 lux; ILT10C, International Light Technologies, Peabody, MA, USA), 3-tier display open-curtain retail case (Model TOM-60DX-BN, Turbo Air Inc., Long Beach, CA, USA). To best simulate retail storage, the curtain was closed at night and open during the day. To prevent bias by cooler position, samples were randomly shuffled daily throughout the experiment as recommended by the American Meat Science Association (AMSA) Guidelines for Meat Color Measurement (King et al., 2023). Case temperature was monitored via a built-in thermometer throughout the day to ensure proper temperature control (3.0°C).

Instrumental color analysis was measured with the HunterLab MiniScan EZ colorimeter, model 45/0 LAV (Hunter Associates Laboratory Inc., Reston, WV, USA), following the AMSA Meat Color Measurement Guidelines (King et al., 2023). The colorimeter was standardized using a black and white tile cover prior to each use and analyzed through the overwrap packaging. Instrumental color values were determined from the average of 3 readings per ground beef sample. Additionally, photographs were taken of each sample to monitor visible color changes during storage.

A meat compositional analysis was completed during sampling to determine protein moisture and fat content in each sample. The analysis was conducted using a near infrared with AOAC (2007.4) approved spectrophotometer (FoodScan^™, FOSS Analytical A/S, Hilleroed, Denmark), and data processing was determined using ISIscan^™ software to analyze the nutrient content in the meat (Martin et al., 2013).

APC were completed each day of storage to determine the rate of microbial growth. A 1.0-g portion of meat, representative of from the top, bottom, and inside of the sample, was placed in a sterile whirl-pak filter bag (Whirl-Pak, Pleasant Prairie, WI, USA). Sterile phosphate buffer (Criterion Dehydrated Culture Media, Hardy Diagnostics, Cat no. C9241, Santa Maria, CA, USA) (10 mL) was added to the meat and hand massaged for 60 s. Serial dilutions were done to measure aerobic microorganisms present in the samples. Dilutions were spread on tryptic soy agar media and incubated at 36°C for 48 h before counting (Criterion Dehydrated Culture Media, Hardy Diagnostics, Cat no. C7141, Santa Maria, CA, USA). Plate counts were log transformed and reported as log CFU/g. The plate counts demonstrated when the samples reached spoilage, indicated by 7-log CFU/g (Alexa et al., 2024).

The development of oxidative rancidity and lipid oxidation in the samples was evaluated throughout the sampling period using 2-TBARS. During each day of sampling, 2 g from each ground beef sample were frozen (−18°C) for later analysis. Samples were frozen and stored for less than 7 d postexperiment. This method provides the most accurate TBARS analysis to ensure methodological consistency (Aguilar Diaz De Leon and Borges, 2020). At the completion of all sampling days, the samples that were frozen after initial analyses were thawed and homogenized (Model PT 10-35, Brinkman Inst., Westbury, NY, USA) into a uniform sample in duplicates and mixed with 8 mL of chilled phosphate buffer (pH 7.0), containing 0.1% ethylenediamine-tetraacetic acid, 0.1% n-propyl gallate, and 2 mL of trichloroacetic acid. The homogenized samples were filtered through a Whatman No. 1 filter paper into glass tubes, and duplicate 2-mL aliquots of clear filtrate were transferred into 10-mL test tubes. Filtrate was mixed with 2 mL of 0.02 M 2-thiobarbituric acid reagent and placed into a hot water bath (100°C) for 20 min. Following the hot water bath, the tubes were transferred to an ice bath for 15 min. The absorbance of each sample was measured at 533 nm with a spectrophotometer (Model UV-1600PC, VWR International LLC, Randor, PA, USA). The values were measured by milligrams of malonaldehyde per kilogram of fresh meat and multiplied using the value of 12.21 according to the standard curve to obtain the TBARS value (Buege and Aust, 1978).

After all samples were collected, the 2.0-g frozen ground beef samples separated for microbiome analysis were thawed for DNA extraction. Prior to extraction, a 2-step centrifugation pre-extraction treatment was completed on the samples prior to DNA isolation to limit the host (bovine) DNA associated with the analysis, as recommended by previous literature (Hultman et al., 2015). Samples were diluted in 10 mL of peptone buffer and hand massaged for 30 s. Supernatant from the homogenization was then placed in a 15-mL conical tube and centrifuged at 200 rcf for 3 min (Eppendorf Centrifuge 5810 R). Following the first centrifugation, 10 mL of the supernatant were removed and centrifuged again at 10,000 rcf for 3 min (Eppendorf Centrifuge 5425) (Hultman et al., 2015). DNA was then extracted from the final pellet using Zymo Research Quick-DNA^™ Miniprep Kit (Lot No. 237579, Irvine, CA, USA), following manufacturer instructions with the optional steps for low biomass. DNA quantification was done with an Invitrogen Qubit^™ Flex Fluorometer.

Extracted DNA from all samples was frozen and shipped on dry ice to Novogene Corporation Inc. (Sacramento, CA, USA) for amplicon sequencing of the V3-V4 region of the ribosomal RNA (rRNA) gene at their China laboratory location (Guangzhou City, Guangdong Province, China). Sequencing was conducted on all samples in 1 lane on the Illumina NovaSeq platform with paired-end 250-bp chemistry (Illumina Inc., San Diego, CA, USA). Microbial 16S rRNA gene amplicons were analyzed using QIIME2 version 2024.2 (Bolyen et al., 2019). The sequences were denoised, and paired reads were joined using DADA2 (Callahan et al., 2016). Taxonomic classification was done using the SILVA 138-99 pretrained classifier with the QIIME feature-classifier plugin (Hung et al., 2022; Bokulich et al., 2018; Quast et al., 2013). Nonmicrobial DNA (chloroplast, mitochondria) was filtered from the dataset. Differential abundance analysis was conducted using the analysis of composition of microbiomes (ANCOM) with bias corrections, with Holm–Bonferroni multiple-test correction to suggest taxa with significant changes (Lin and Peddada, 2020). Alpha diversities (observed features, Faith’s phylogenetic diversity, and Shannon’s diversity) were statistically compared using a Kruskal–Wallis test, with Benjamini–Hochberg multiple-test correction (Faith, 1992). Beta diversities were analyzed using the Bray–Curtis dissimilarity metric, unweighted UniFrac, and principal coordinates analysis (PCoA) analysis was analyzed through QIIME2 Emperor (Lozupone et al., 2011). Furthermore, changes seen throughout the shelf life were calculated using the QIIME2 longitudinal plugin (Bokulich et al., 2018) with the volatility and linear mixed-effects tests (Xia and Sun, 2023; Xia et al., 2025), and the statistical assessments were made using the permutational multivariant analysis of variance test (Lozupone and Knight, 2005; Anderson, 2017).

The predictability of microbial succession during ground beef spoilage was assessed using a machine learning model constructed using the random forest algorithm (Zhang and Lu, 2012; Metcalf et al., 2015). The random forest algorithm is appropriate for this situation to manage the large number of predictive features while reducing, as much as possible, the likely overfitting issue due to a limited samples size. The microbiome relative abundance data were used as the predictor variable, with spoilage metric as the response variable. Models were trained to predict the following: (1) the day of retail display (day of spoilage), (2) APC levels, and (3) TBARS value. These were constructed using a K-fold cross-validation with hyperparameter tuning, as is appropriate for a situation with small sample sizes to reduce overfitting (Papoutsoglou et al., 2023). Grouping was applied such that all microbial data within a sample (production lot) of ground beef were placed in a validation fold together, and a random state parameter was applied to support feature stability across the folds. The accuracy of the models was assessed during cross-validation using the mean absolute error (MAE), defined as the deviation of the predicted from observed values represented in the same unit as the original data ([1] days; [2] log CFU/g; [3] mg of malonaldehyde/kg). Recursive feature elimination was used to identify the number and name of features (microbial amplicon sequence variants [ASV]) that contributed to the most accurate final model. For each model, the taxonomic names associated with the important ASV were assigned for comparison with other analyses. All modeling was conducted using the python machine learning package scikit-learn version 19.0 (Pedregosa et al., 2011). Results were visualized using ggplot2 in R (Wickham, 2016).

Statistics for nonmicrobiome analyses were analyzed using R Studio version 4.2.2 using the packages emmeans, tidyverse, and ggplot2 (R Core Team, 2022; Wickham, 2016; Wickham et al., 2019; Lenth et al., 2025). Means were compared using linear mixed-effects models with day as a fixed effect and sample as a random effect. Significance was identified at a P value less than .05. Statistics for microbiome analysis were analyzed using QIIME 2 and visualized using the ggplot2 package in R Studio version 4.2.2 (Bolyen et al., 2019; R Core Team, 2022; Wickham, 2016).

This experiment monitored the quality and microbial changes in ground beef products over the shelf life to understand the microbial succession patterns. Overall, quality data demonstrated progressive visual color changes and spoilage patterns throughout the experimental period (Supplemental Figure 1). As expected, there was little variation in the meat compositional analysis over the experimental period. There was significance among protein (P < .01; CI = 0.19–0.32) and fat (P < .01; CI = 0.17–0.27) by day, but no changes in moisture (P = .622) possibly due to the aerobic packaging and oxygen availability (Supplemental Table 2). The compositional analysis shows a statistical difference over time, but this shift is not biologically significant and does not impact product quality. Instrumental color, APC, and TBARS were monitored as a metric for product spoilage. There were no differences (P > .05) for full analysis of instrumental color, APC, or TBARS among samples of ground beef per each day of simulated storage. There was a decrease (P < .01; CI = 1.44–2.43) in the surface lightness (L*), redness (a*), and yellowness (b*) throughout storage until day 9, at which point the color became stabilized (Figure 1A, Supplemental Table 1). APC increased (P < .01; CI = 0.38–0.47) over time, crossing the threshold for the definition of spoilage (7 log CFU/g) on day 6 for all samples (Figure 1B). Similarly, TBARS, quantified by milligrams per kilograms of malonaldehyde, showed an increasing trend over time. The TBARS increased (P < .01; CI = 0.61–0.97) from day 0 to day 8, then decreased (P < .01) until day 11, before increasing until the end of the experimental period. TBARS crossed the threshold for spoilage defined by 2.0 mg/kg of malonaldehyde by oxidative rancidity on day 6 (Figure 1B) (Hernández-García et al., 2022).

The microbiome of ground beef throughout the shelf life was evaluated using 16S rRNA microbial gene sequencing. Sequencing a total of 45 samples resulted in a total of 8,951,852 demultiplexed sequences. The reads were denoised into 656 ASV, with 641 ASV remaining after filtering for nonmicrobial taxa (mitochondria and chloroplast). Samples were rarefied at 105,506 sequences for diversity analysis, using 35 out of 45 rinsate samples.

There were no differences among the samples for richness (mean = 39.9; P = .097), Shannon’s diversity (mean = 3.44; P = .115), and Faith’s diversity (mean = 13.1; P = .113) (Figure 2A). However, β diversity, as measured by unweighted UniFrac distances, increased over time (P = .009). Furthermore, there was a clear difference in β diversity between unspoiled and spoiled samples, as defined by APC (P < .05) (Figure 2B). When clustered using PCoA analysis, UniFrac dissimilarity shows clustering based on the spoilage state by each experimental day (Supplemental Figure 2).

These shifts in β diversity are also shown by the taxonomy. Ground beef has an initial microbiome that changes throughout the shelf life. As the shelf life progressed, the most dominant increase was in Pseudomonadaceae, which is a known leading spoilage organism (Taormina, 2021). There are also higher populations of Carnobacteriaceae, Listeriaceae, Streptococcaceae, Lactobacillaceae, and Yersiniaceae throughout the shelf life. Although many organisms are present on the day of grinding, there is a decrease in the number of bacterial families present throughout the experiment, as the organisms outcompete each other in this time, fighting for limited resources (Figure 2C). Additionally, when analyzing the differential abundance through ANCOM, there were significant changes among the taxa that are significantly enriched in the spoiled condition compared to the nonspoiled reference samples (P < .01). Results of ground beef samples demonstrate the top-enriched taxa displayed the greatest log fold change during storage from day 0 to 14 (Table 1).

Changes in the microbial community composition were predictable. Machine learning models were constructed and validated using K-fold cross-validation with hyperparameter tuning to determine the association of the microbiome with the changes in ground beef quality throughout shelf life. Random forest regression models using the microbial community with 642 microbial features as the predictor and the quality indicator and 3 samples per day measured for quality as the response variable were successfully constructed and evaluated. The MAE for the models were 1.09, 2.98, and 2.48 for predicting APC, TBARS, and days of spoilage, respectively. The most informative features (ASV) for each model were assigned taxonomic identifiers (Figure 3). The APC most optimal model contained 15 features, which included Methyloceanibacter, uncultured Micropepsaceae family, Flavobacteriaceae, Prevotella, Puniceicoccaceae, Enterococcaceae, Yersinia, Pseudomonadales clade SAR86, Acidobacteriaceae subgroup 1, Puniceicoccaceae, Latilactobacillus, Carnobacterium, Akkermansia, Varibaculum, and Longimicrobiaceae. The TBARS most optimal model contained 6 features, which included Methyloceanibacter, uncultured Micropepsaceae, Flavobacteriaceae, Prevotella, Puniceicoccaceae, and Enterococcaceae. The day model most optimal model contained 9 features, which included Methyloceanibacter, uncultured Micropepsaceae, Flavobacteriaceae, Prevotella, Puniceicoccaceae, Enterococcaceae, Yersinia, Pseudomonadales SAR86 clade, and Acidobacteriaceae subgroup 1.

This research study had many aspects that allowed for a new look into microbial meat spoilage, and we were able to see the different microbial communities changing throughout the shelf life of ground beef. The greatest limitation of this study was in the small sample number, which allowed preliminary findings in the area but not generalizable conclusions. It is difficult to define an optimal sample size for microbiome studies, especially when machine learning is employed; however, now that this study has demonstrated relatively large effects sizes associated with spoilage shelf life, these can be calculated with power analysis for future study. The sample size generally used to define the threshold for a “small study” in machine learning is 100 samples per treatment, which may provide practical challenges in conducting a study of this type. The use of cross-validation methods can help ameliorate this challenge, and the work could be done over several time periods to allow more samples even with limited physical space. Despite this limitation, the results presented here will serve as a basis for future work, investigating trends in microbial succession and their relationship with quality outcomes. Beyond increasing the sample size, more diverse sources for the meat samples that provide more variation in the initial microbiome may have allowed for more detailed results and models with less potential overfitting. Additional work may also employ other omics methods to investigate the fundamental trends underlying the observed patterns presented here.