Spatial and temporal dynamics of microbiomes and resistomes in broiler litter stockpiles.

Farmers apply broiler chicken litter to soils to enrich organic matter and provide crops with nutrients, following varying periods of stockpiling. However, litter frequently harbors fecal-derived microbial pathogens and associated antibiotic resistance genes (ARGs), and may be a source of microbial contamination of produce. We coupled a cutting-edge Loop Genomics long-read 16S rRNA amplicon-sequencing platform with high-throughput qPCR that targeted a suite of ARGs, to assess temporal (five time points over a 60-day period) and spatial (top, middle and bottom layers) microbiome and resistome dynamics in a broiler litter stockpile. We focused on potentially pathogenic species from the Enterobacteriaceae, Enterococcaceae and Staphylococcaceae families associated with food-borne disease. Bacterial diversity was significantly lower in the middle of the stockpile, where targeted pathogens were lowest and Bacillaceae were abundant. E. coli was the most abundant Enterobacteriaceae species, and high levels of the opportunistic pathogen Enterococcus faecium were detected. Correlation analyses revealed that the latter was significantly associated with aminoglycoside (aac(6')-Ib(aka aacA4), aadA5), tetracycline (tetG), vancomycin (vanC), phenicol (floR) and MLSB (mphB) resistance genes. Staphylococcaceae were primarily non-pathogenic, but extremely low levels of the opportunistic pathogen S. aureus were detected, as was the opportunistic pathogen S. saprophyticus, which was linked to vancomycin (vanSA, vanC1), MLSB (vatE, ermB) and tetracycline (tetK) resistance genes. Collectively, we found that stockpile microbiomes and resistomes are strongly dictated by temporal fluctuations and spatial heterogeneity. Insights from this study can be exploited to improve stockpile management practice to support sustainable antimicrobial resistance mitigation policies in the future.

Introduction

Broiler litter (BL), a by-product of meat production, is a mixture of fecal droppings, bedding (usually sawdust, wood shavings, rice hulls or straw), feathers and waste feed [1], [2]. It enriches soil organic matter and provides crops with nitrogen, phosphorus and other essential nutrients, and is therefore globally applied as an organic fertilizer [3]. For certain crops, application of BL has been reported to be more effective than synthetic fertilizers [4], [5]. However, it may harbor fecal-derived microbial pathogens and associated antibiotic resistance genes (ARGs), especially when considering the high volume of antibiotics used in intensive poultry farming in many parts of the world [6], [7], [8], [9]. Previous reports have documented the presence of foodborne pathogens such as E. coli, Salmonella spp. and Campylobacter spp. in BL [9], [10], [11]. It can therefore be a source of transmission to animals, humans and the environment, considering the capacity of at least some of these pathogens to persist for months in water, soil and crops [12], [13], [14], [15]. Furthermore, several studies observed slightly increased levels of both antibiotic resistance bacteria (ARB) and ARGs (even in the absence of antibiotics) in manure-amended soils, however, the level of persistence of these bacteria/genes between different studies is contradictory [16], [17], [18]. Thus, while, animal manure is traditionally used as a soil conditioner and fertilizer for crop production, the application of BL to soil is a potential route for environmental exposure of antimicrobial resistance (AMR) [19].

The presence of potentially pathogenic bacteria and associated ARGs in BL, necessitates manure stabilization practices such as anaerobic digestion and composting, to mitigate dissemination of AMR to the environment [20], [21]. Composting is a practical method for eliminating human and zoonotic pathogens from BL [22], but proper composting protocols that generate high temperatures in compost stockpiles are required for inactivating pathogens in BL or BL-based organic fertilizers prior to land application [23], [24], [25]. Based on international regulations, to ensure effective pathogen elimination, all compost particles need to be exposed to a minimum temperature of 55 °C for a period of at least three consecutive days [26], [27]. However, pathogens surviving the composting process or any phase of stabilization may regrow during storage or following land application. Indeed, several studies have demonstrated the persistence of zoonotic pathogens in the finished compost at different levels of maturity and in compost-amended soils [28], [29], [30].

Most farmers apply BL without processing or after partial stabilization in a stockpile [24], [31], [32]. Although static stockpiling facilitates partial thermal inactivation, such minimal management regimes are not as effective in pathogen removal as controlled thermophilic composting due to high levels of heterogeneity [10]. Holistic understanding of temporal and spatial dynamics of microbiomes and resistomes in commercial BL stockpiles is imperative in order to minimize risks and enable development of best management practices.

This study comprehensively analyzed the diversity and composition of bacterial communities and ARGs in a commercial BL stockpile at three depths over a 60-day period. We focused on bacterial pathogen-encompassing species from the Ennterobacteriaceae, Enterococcaceae and Staphylococcaceae families using a unique Loop Genomics long-read sequencing platform that couples standard Illumina short-read sequencing with a unique molecule barcoding technology (Loop Genomics, USA), generating full length (∼1500 bp) 16S rRNA gene sequences. In tandem, we applied High Throughput qPCR (HT-qPCR) that targeted a suite of ARGs and mobile genetic elements (MGEs).

Materials and methods

Stockpile materials and experimental design

The field study was conducted in Newe Ya’ar Research Center, Israel, between August-September 2018 (ambient temperatures 26.5–28.4 °C). An open static pile ('stockpile') of ca. 35 m3 was set with BL from the poultry farm of Moshav Balfouria, Northern Israel. It was obtained at the end of an extensive indoor rearing period of six weeks, and was characterized for water content (21.4%; SD = 0.81), pH (6.95; SD = 0), electrical conductivity, EC (9.9 dS/m; SD = 0), and ash content (14.87%; SD = 0.66); see analytical methods in Section 2.2.

The stockpile was prepared by unloading BL from a large truck onto a concrete floor outside. This pile essentially composed of three sub-piles, which were monitored in five sampling campaigns over week 1 (day 7), week 2 (day 14), week 3 (day 21), week 4 (day 28) and week 8 (day 60) after the initial BL unloading. During each campaign samples were collected from the top, middle and bottom layers of the stockpile for physiochemical analyses and DNA extraction (see below), as shown in Fig. S1A. The average stockpile height decreased from 113 cm to 77 cm within this period (Fig. S1B).

Physicochemical analyses

Temperature, water-content, pH, EC, gas phase oxygen and gas-phase ammonium were measured in three replicates at three depths for each of the sampling campaigns. Temperature was measured by a mobile type K thermocouple that was constructed on 80-cm long stainless-steel rod and connected to a 305 thermometer (Elcon LTD, Israel). BL samples were collected from the designated locations using a soil auger. Water content (w/w) was determined on ca. 10 g samples by oven drying at 70 °C for 48 h. The pH and EC were measured in aqueous extracts at 1:9 dry weight: deionized water using a reciprocal shaker at 200 RPM for 1 h. The pH was analyzed directly in the suspension (LL-Ecotrode Plus WOC; Metrohm, Herisau, Switzerland), while the EC was determined in the supernatant after centrifugation at 6000 RPM for 20 min at 25 °C (CyberScan CON 11, Eutech Instruments, Thermo Fisher Scientific Inc., Waltham, Massachusetts, USA). Gas-phase oxygen and ammonium were monitored in the same locations by pumping the headspace using two pocket pumps (SKC 210 and SKC 220-1000TC, 84, PA, USA). The air was pumped through a Teflon tube (1/4″ inner diameter) and was distributed to an acid trap (for ammonia analysis) and to an oxygen sensor (SO-210, Apogee Instruments, Inc., Logan, UT, USA). Ammonia (trapped in 0.1 M sulfuric acid) was analyzed using the spectrophotometric method with slight modifications after Willis et al. (1996) [33] and similar to Avidov et al. (2021) [10].

DNA extraction

Samples from the BL stockpile were collected at three depths (top, middle and bottom layers) at the five sampling campaigns described above. Triplicate samples for each depth were combined into composite samples, resulting in a total of 15 composite samples (5-time points, 3 sampling depths) for DNA extraction and subsequent analyses of bacterial community and ARG composition (see below). BL samples from each of the above points were collected in 50 mL sterile tubes, transferred on ice to the laboratory and stored at −80 °C until DNA extraction. After thawing on ice, DNA was extracted from 0.5 g of BL sample with the GenALL DNA extraction kit (GeneAll Biotechnology Co. Ltd., Seoul, South Korea). Briefly, samples were lysed twice by bead beating (FastPrep FP120, BIO101 Savant, La Jolla, USA), for 45 s, at 5.5 intensity level in two cycles, and then centrifuged at 20 °C, 14,000 rpm, for 10 min. Supernatant was combined and DNA was extracted according to the kit’s protocol. DNA concentration and quality were determined using a NanoDrop 2000c Spectrophotometer (Thermo Fisher Scientific, Wilmington DE., US) and a Qubit Fluorimeter (Invitrogen, Carlsbad, CA, US), and validated by gel electrophoresis. Finally, the DNA extracts of all samples were stored at −20 °C until further analysis.

High-throughput sequencing and data analysis

Compared to commonly applied amplicon sequencing that targets one or two variable regions of the 16S rRNA gene, Loop Genomics synthetic long-read sequencing provides high quality base-resolution and accurately classifies species and sometimes potential strain levels by reducing false positives [34]. Full-length 16S rRNA gene libraries were constructed using 10 ng of genomic DNA using the LoopSeq 16S long-read kit from Loop Genomics (protocols available at www.loopgenomics.com). Two DNA tags (i) a Unique Molecular Identifier (UMI) to each unique “parent” molecule and (ii) a sample-specific tag (i.e. a Sample Index) equally to all molecules in the same sample. Barcoded molecules were amplified, multiplexed and each UMI was distributed to a random position within each parent molecule. Molecules were then fragmented into smaller units at the position of UMI insertion, creating a library of UMI-tagged fragments compatible with an Illumina sequencing platform run in 150 bp paired-end (PE) mode. 100–150 M PE reads (50–75 M clusters passing filter) were used for each sequencing run, yielding ∼20 gigabases (Gb) of data. The short-read raw data were collected in real-time on Illumina’s BaseSpace, which generates FASTQ file and then were uploaded to the Loop Genomics cloud-based platform for processing raw short-reads into assembled synthetic long contigs. Reads were trimmed employing Trimmomatic [35] to remove adapter sequences before they were de-multiplexed based on their Loop Sample Index, which groups reads by the sample from which they originated. Sample-specific reads were then binned by UMI and processed collectively through SPAdes [36]. Processed full-length 16S rRNA genes were then clustered into ASV (Amplicon Sequence Variant) bins of 100% sequence homology by applying the DADA2 package [37]. The generated representative ASVs were then used to assign bacterial taxonomy using the SILVA rRNA database [38], and the database entry with the highest mapping score was assigned to that ASV. The full-length 16S amplicon sequencing data generated for this study were uploaded to the NCBI data repository, under the primary accession number BioProject ID: PRJNA777722.

High-throughput qPCR analysis

To characterize the temporal and spatial BL stockpile resistome, HT-qPCR was performed using a WaferGen SmartChip Real-Time PCR System (Fremont, CA, USA) on each of the composite samples. A total of 296 verified primer sets were selected for HT-qPCR reactions as reported previously [39], [40] to specifically target 283 ARGs conferring resistance to nearly all major types of antibiotics, 12 MGE marker genes including 8 transposase encoding genes, 4 integron-integrase genes and a 16S rRNA gene used for normalization. The DNA samples were transferred to the Key Lab of Urban Environment and Health, Institute of Urban Environment, Chinese Academy of Sciences, where HT-qPCR was performed. A dye-based PCR assay was employed. Briefly, PCR mixtures (100 nL per well) consisted of 1 × LightCycler 480 SYBR Green I Master (Roche Applied Sciences, Indianapolis, IN), nucleasefree PCR-grade water, 1 mg/mL bovine serum albumin (New England Biolaboratories, Beverly, MA), 500 nM of each primer and a DNA template of 4 ng/μL. The PCR cycle consisted of 10 min at 95 °C, followed by 40 cycles of denaturation at 95 °C for 30 s and annealing at 60 °C for 30 s. The Wafergen software automatically generated melting curve analyses. The qPCR results were analyzed using SmartChip qPCR software (2.7.0.1 version), and wells with multiple melting peak and amplification efficiencies beyond the range (1.8–2.2) were discarded. A threshold cycle (Ct) of 31 was used as the detection limit, and only samples with all three technical replicates being amplified were regarded as positive and used in further data analysis. Relative copy numbers of the targeted genes were calculated as described previously [41], [42].

Statistical and network analysis

All statistical analyses were performed in the R software package (4.0.3 version). Kruskal-Wallis tests followed by multiple pairwise comparisons using Nemenyi post-hoc were applied to the alpha‐ diversity metrics to assess statistically significant differences in the bacterial diversity within samples between groups. One-way analysis of variance (ANOVA) and non-parametric Mann-Whitney U tests were conducted in GraphPad Prism 8 for significance analysis (p < 0.05 was considered as significant). Nonmetric multi-dimensional scaling (NMDS) based on Bray-Curtis dissimilarity matrices, PERMANOVA test, Mantel tests and Procrustes analysis was performed using “vegan” package in R [43]. Graphs were generated using “ggplot2” or “pheatmap” packages in R.

For the co-occurrence network modeling: a correlation matrix was first calculated based on normalized abundance data using Spearman’s rank correlation (rho) between resistance genes obtained from HT-qPCR array analysis. Only statistically strong (rho > 0.80) and significant (padj < 0.01) correlations were selected for this purpose [44]. The co-occurrence network was visualized employing Gephi software package (v0.9.2) [45].

Results

Physiochemical parameters

The temporal and spatial dynamics of temperature, water-content, pH, EC, gas phase oxygen and ammonia in the stockpile are described in Figure S2. The average temperature was highest in the middle layer and decreased towards the top and bottom of the stockpile. Over time, the temperature significantly dropped in the top layer from 50.6 to 39.7 °C (p < 0.01). Temperature also slightly decreased in the middle layer (from 54.9 to 51.0 °C), while it increased in the bottom layer over time from 45.3 to 55.1 °C and then decreased to 51.3 °C. However, these results were not statistically significant (p > 0.05). Thus, although the pile was not at optimal water-content and was only passively aerated, it is apparent that microbial activity still generated increased temperatures. Water-content decreased over time, from 22.2 to 10.5%, 34.2 to 29.2%, and from 35.1 to 30.3% in the top, middle and bottom layers, respectively. The decrease in water-content with time was only statistically significant for top layer samples (p < 0.05). The vertical gradient (from dry to moist) indicates both evaporation from the top layer, and water production via organic matter oxidation and accumulation in the bottom layer. The EC ranged from 8.7 to 12.3 dS/m and it was linearly correlated with water content (r = 0.89, p < 0.01) presumably due the extraction of manure salts during water movement in the stockpile. The pH did not significantly change in any of the analyzed depths In contrast, gas phase oxygen measurements showed a strong vertical gradient, ranging from 19.9% at the top of the pile to 4.1% at the bottom on day 7. On days 21–28, oxygen values were below detection limits in the bottom layer. Gas-phase ammonium concentrations did not show a clear spatial trend but were generally higher on day 60.

Microbial community composition

An average of 82,995 synthetic long-reads (SLR) per sample were generated using the Loop Genomics SLR NGS platform, with an average of ∼45,000 high quality SLR counts that were phylogenetically characterized using the SILVA reference database (Table S1), which were grouped into 648 unique bacterial taxa. The effect of time and depth on BL stockpile bacterial (α-) diversity and community composition (β-diversity) was assessed using Shannon diversity and Bray-Curtis distance indices, respectively. The age of the pile did not significantly influence bacterial diversity, as only a slight non-significant reduction was observed on day 60. In contrast, bacterial diversity was significantly lower in the middle of the stockpile relative to the top (p < 0.05; Fig. 1A). Correlation analysis revealed that this spatial reduction in α-diversity was primarily associated with higher temperature (r = −0.51, p < 0.05). Higher water-content (r = −0.26) and EC (r = −0.38), and lower gas-phase oxygen levels (r = 0.32) were also correlated to this reduction in diversity, but they were not statistically significant (p > 0.05; Fig. 1B). These associations are in line with co-correlations observed between the physicochemical properties (EC vs. water-content, r = 0.89, p < 0.01; water-content vs. temperature, r = 0.64, p < 0.05; EC vs. temperature, r = 0.67, p < 0.05; and oxygen vs. temperature, r = −0.45, p < 0.05).

Fig. 1

Microbial community diversity, structure and dynamics in broiler litter stockpiles. (A) Temporal dynamics and spatial differences in bacterial (alpha) diversity as measured by the Shannon diversity index in the stockpile. Values are presented as the median (black horizontal line), lower and upper hinges correspond to the 25th and 75th percentiles. Outliers are displayed by small black dots. Different letters represent statistically significant differences between groups determined by the Kruskal-Wallis test followed by pairwise comparisons using the Nemenyi post-hoc method, (B) Correlations between bacterial diversity (Shannon index) and selected physiochemical properties, (C) Non-metric multidimensional scaling (NMDS) plot of litter stockpile samples based on Bray-Curtis distances showing the temporal and spatial differences in bacterial community composition (Stress = 0.17, PERMANOVA; R2 = 0.419, F = 1.805 and P < 0.01), and (D) Temporal and spatial impact on the relative abundance of bacterial families in the litter stockpile over time. Unclassified taxa and bacterial families with relative abundances of<0.5% were grouped into “Others”. T = top, M = middle, and B = bottom samples.

In contrast to the α-diversity described above, the bacterial community composition (β-diversity) significantly varied with stockpile age (Bray-Curtis, PERMANOVA test by adonis2; R2 = 0.419 and p < 0.01), and this temporal effect was more distinct than the stockpile depth effect (Fig. 1C). Among the 648 unique bacterial taxa identified, 214 (33.02%) and 293 (45.22%) were ubiquitous in all time points and depths, respectively (Figure S3 and Table S2). In contrast, a substantial proportion of unique bacteria were detected in the litter stockpile on Day 7 (77 unique taxa, 17.2%) and in the Bottom (135 unique taxa, 25.1%) samples, suggesting the time and depth could be key drivers in shaping the stockpile microbiota.

Taxonomic profiling at the family level revealed that Bacillaceae (32.35%), Dermabacteraceae (14.73%) and Staphylococcaceae (12.56%) represented over 50% of the classified bacterial families in the BL stockpile (Table S3; Fig. 1D). The relative abundance of Bacillaceae significantly increased over time (5.83 ± 0.68% on Day 7 vs. 46.17 ± 8.82% on Day 60; p < 0.05), whereas Dermabacteraceae decreased (30.7 ± 3.11% on Day 7 vs. 9.97 ± 4.37% on Day 60; p < 0.05). The relative abundance of Bacillaceae was higher in the middle of the stockpile, whereas that of Dermabacteraceae was lower. However, this depth-based effect was not significant (p > 0.05). We specifically focused on the relative abundance and composition of Staphylococcaceae, Enterococcaceae and Enterobacteriaceae, because these families contain opportunistic or obligatory pathogens. The relative abundance of Enterococcaceae (1.87% vs. 0.79%), Enterobacteriaceae (1.11% vs. 0.54%), and Staphylococcaceae (17.55% vs. 13.54%) decreased between Day 7 and Day 60 (p > 0.05), but the temporal dynamics were not linear (Fig. S4). Spatial effects on the abundance of these families were also apparent. While Enterococcaceae (2.24 ± 0.71% in top, 0.92 ± 0.24% in middle, and 1.34 ± 0.72% in bottom; p < 0.05) and Enterobacteriaceae (0.82 ± 0.42% in top, 0.62 ± 0.93% in middle, and 0.7 ± 0.56% in bottom; p > 0.05) abundance was lower in the middle of BL stockpile, depth did not substantially influence the relative abundance of Staphylococcaceae (12.72 ± 1.99% in top, 13.19 ± 5.75% in middle, and 11.77 ± 4.44% in bottom, p > 0.05).

Several opportunistic or obligatory pathogen-associated species were identified within the three targeted families described above. Escherichia coli (78.25% of Enterobacteriaceae and 0.57% of total characterized bacterial reads) was the most dominant Enterobacteriaceae species in litter stockpiles, but significantly lower abundances of Proteus mirabilis (8.92% of Enterobacteriaceae and 0.08% of total bacterial reads), Salmonella enterica (1.81% of Enterobacteriaceae and 0.004% of total bacterial reads) and Klebsiella pneumoniae (1.71% of Enterobacteriaceae and 0.005% of total bacterial reads) were also detected. (Fig. 2A). Most Enterococcus species were non-pathogenic fermenters (Fig. 2B), however, significant levels of the opportunistic enteric human pathogen Enterococcus faecium (29.18% of Enterococcaceae and 0.45% of total bacterial reads) were detected, as well as lower abundances of Enterococcus faecalis (1.95% of Enterococcaceae and 0.024% of total bacterial reads) and the zoonotic pathogen Enterococcus hirae (8.74% of Enterococcaceae and 0.11% of total bacterial reads), involved in growth depression and endocarditis in chickens, and blood infections in humans [46], [47], [48]. Non-pathogenic species, including Jeotgalicoccus sp. (25.05% of Staphylococcaceae and 3.10% of total bacterial reads), Staphylococcus cohnii (19.27% of Staphylococcaceae and 2.68% of total bacterial reads), Staphylococcus sp. (16.36% of Staphylococcaceae and 2.00% of total bacterial reads) and Staphylococcus lentus (13.82% of Staphylococcaceae and 1.61% of total bacterial reads), represented ∼75% of all Staphylococcaceae (Fig. 2C). However, extremely low levels of the opportunistic pathogen Staphylococcus aureus (0.07% of Staphylococcaceae and 0.009% of total bacterial reads) were also detected in the BL stockpile.

Fig. 2

Relative distribution of bacterial species from three targeted pathogen-associated families in broiler litter stockpile samples. (A) Enterobacteriaceae, (B) Enterococcaceae, and (C) Staphylococcaceae. Unclassified, and bacterial species with a relative abundance per family of<0.5% were grouped into “Others”. T = top, M = middle, and B = bottom samples.

Using Spearman’s rank correlation analysis, we investigated potential correlations between selected physicochemical factors (i.e., temperature, water-content, pH, EC, ammonium and oxygen) and the abundance of species from the three-priority pathogen-harboring families described above (Fig. 3). Collectively, the abundance of most species correlated negatively with EC, water-content and temperature, and positively with oxygen and pH. Among the priority pathogens E. faecium showed a significant positive correlation with pH (r = 0.63, p = 0.011) and oxygen (r = 0.65, p = 0.008), and a negative correlation with EC (r = −0.53, p = 0.043). Other enterococci such as E. durans (r = 0.52, p = 0.044), E. casseliflavus (r = 0.65, p = 0.009), and E. eurekensis (r = 0.68, p = 0.005) also showed strong positive correlation with oxygen. A similar, albeit not statistically significant (p > 0.05) pattern was observed for E. coli and Shigella spp. (also facultative anaerobes). In contrast, S. enterica, K. pneumoniae and P. mirabilis showed a slight negative correlation with oxygen and a positive correlation with EC and water-content. A similar but not significant trend was also observed for Staphylococcus aureus and S. lentus.

Fig. 3

Correlation between selected physiochemical parameters and detected bacterial species from three pathogen-harboring families in the broiler litter stockpiles. Only bacteria detected in at-least 50% of the samples (8 out of 15 samples) were used for correlation analysis applying Spearman’s rank correlation method. Statistically significant correlations (p < 0.05) are marked by asterisk (*) in the heatmap.

Occurrence and diversity of antibiotic resistance genes

A total of 180 unique genes including 173 ARGs (conferring resistance to nine major classes of antibiotics including aminoglycosides, β-lactams, multidrug, macrolide-lincosamide-streptogramin B (MLSB), phenicols, sulfonamides, tetracycline, trimethoprim, and vancomycin) and 7 MGEs (4 integrase and 3 transposase genes) were detected in the different depth and time points of BL stockpile (Table S4). Collectively, the relative abundance (copy number/16S rRNA gene copy number) of ARGs conferring resistance to phenicols (2.01 ± 0.51), aminoglycosides (1.15 ± 0.28), MLSB (1.19 ± 0.31) and tetracycline (0.55 ± 0.15), and multidrug efflux pumps (0.51 ± 0.14), were most abundant. However, a variety of β-lactam (28 ARGs with a relative abundance of 0.005 ± 0.002) and vancomycin resistance genes (22 ARGs with a relative abundance of 0.001 ± 0.0009) were also detected, albeit at very low concentrations (Table S4). The diversity of the litter stockpile resistome and mobilome increased during the first three weeks, but remained relatively stable between Day 21 and Day 60. In contrast, there was no significant effect of depth on ARG/MGE diversity (Fig. 4A). Similar to the bacterial community profile, NMDS analysis showed a clear and robust temporal effect (Bray-Curtis, PERMANOVA test by adonis2; R) on ARGs that was most prominent between day 7 and 14, contrasted by a very subtle spatial effect (Fig. 4B). Despite the general increase in gene diversity over time, for most stockpile depths, the relative abundance of ARGs and MGEs decreased with the time (7.52 ± 0.41 on Day 7 vs. 6.06 ± 1.19 on Day 60; p > 0.05; Fig. 4C).

Fig. 4

Temporal dynamics and spatial differences in antibiotic resistance gene diversity, structure and composition in the broiler litter stockpile. (A) ARG richness analysis showing unique ARGs and MGEs detected in the stockpiles over time and at different depths. The ARG richness values are presented as the median (black horizontal line), lower and upper hinges correspond to the 25th and 75th percentiles. The outliers are displayed by small black dots, (B) Non-metric multidimensional scaling (NMDS) plot of litter stockpiles samples based on Bray-Curtis distances showing the ARG composition over time (Stress = 0.068, PERMANOVA; R2 = 0.454, F = 2.079 and p < 0.05), (C) relative abundance (normalized to per copy of 16S-rRNA gene) of ARGs and MGEs observed in litter stockpiles over time at different depths, and (D) The number of persistent and fluctuating resistance genes detected for each resistance classes. The ARG and MGE detected in 13 samples and more (13 to 15 samples) are classified as Persistent resistance genes, however ARGs and MGEs detected in at-least in 1 sample and<13 samples (1 to 12 samples) are categorized as fluctuating resistance genes. T = top, M = middle, and B = bottom samples.

We characterized 112 (62.2%) ARGs as persistent, because they were detected in at least 13 of the 15 samples. The 68 (37.8%) remaining genes were classified as fluctuating ARGs (Table S4; Fig. 4D). The “persistent resistome” included all phenicol (7) and trimethoprim (3) resistance genes and a large fraction of tetracycline (19, 86.4%), aminoglycoside (17, 76.9%), and MLSB (21, 72.4%) ARGs, as well as integrases (3, 75%). In contrast, more than half of the genes conferring resistance to β-lactams (16, 57.1%), sulfonamide (2, 66.7%) and vancomycin (13, 59.1%) were part of the fluctuating resistome. Persistent ARGs were more abundant in early time points (Figure S5), and cmx(A) (phenicol resistance), tnpA (transposase), qacEdelta1 (multidrug resistance), ermC (MLSB resistance) and tetL (tetracycline resistance) were the most abundant persistent ARGs. Interestingly, multidimensional scaling analysis based on Bray-Curtis distances revealed that the temporal effects observed in BL stockpiles (as described above) was strongly dictated by persistent resistance genes (PERMANOVA; p = 0.036, Figure S6A), and not by fluctuating genes (PERMANOVA; p = 0.227, Figure S6B), suggesting that these genes are more deterministic vs. the fluctuating genes that are more stochastic.

Co-occurrence network analysis of persistent genes

Co-occurrence patterns among persistent ARGs and MGEs were explored using network inference modeling based on statistically strong (rho > 0.08) and significant (adjusted p-value, padj < 0.01) correlations as described previously [49], [50]. The generated network consisted of 80 nodes (persistent resistance genes) and 173 edges (connections/associations) (Fig. 5). The network modularity index was 0.776, and based on the modularity class, the entire network could be parsed into 12 modules (i.e. clusters of nodes), with half of the total nodes (40 out of 80) occupied by three large modules: Module I, II and III. The co-occurring genes of the module hubs (i.e. most densely connected node) are listed in Table S5. Aminoglycoside resistance-conferring genes aadA and aadD were hubs of modules I and II, respectively, and a multidrug efflux pump-encoding gene, yceL_mdtH was the hub of module III (Fig. 5). Network analysis shed light on potential linkage between MGEs and associated ARGs. For example, module I contained the transposase gene tnpA, which clustered with tetD, aadA2, qacEDelta1, cmx(A) and the hub gene aadA. The class 2 integrase gene, intI2 clustered with the blaPSE β-lactamase [51], previously detected in Salmonella enterica [52], and the tetracycline resistance gene tetH, encoding for an efflux pump frequently associated with MGEs. Interestingly, the sulfonamide-resistance gene, sul2 clustered with three tetracycline resistance genes (tetX, tetG, tetR), the β-lactamase blaOXA10, a multidrug efflux pump (floR), and an aminoglycoside-resistance gene (aac(6′)-Ib(aka aacA4)). While the co-occurrence network analysis relies on strong statistical associations between gene abundance data, it cannot validate gene linkage, and therefore should be taken carefully when attempting to link explicit ARGs and MGEs.

Fig. 5

Co-occurrence network linking persistent ARGs and MGEs detected in broiler litter stockpiles. Nodes represent ARG/MGE subtypes and edges (i.e., connections in between) indicate strong and significant (padj < 0.01) pairwise correlations (rho > 0.80). The size of each node is proportional to the number of connections (i.e. degree). The obtained network consisted of 80 nodes and 173 edges with a modularity index of 0.776, indicating that obtained network had a modular structure.

Correlation between resistance genes and microbial communities in broiler litter stockpiles

Mantel tests and Procrustes analysis were performed to understand the relationship between ARGs (from HT-qPCR) and taxon-based bacterial community structure (at species or where not possible genus level) based on amplicon sequencing data (Table S6). Mantel tests failed to reveal correlations between specific ARGs and bacterial community composition (r = 0.294, p > 0.05, 9999 permutations). However, correlations between bacterial communities and persistent ARGs were significant (p < 0.05, 9999 permutations), albeit at low r values (r = 0.2985) suggesting that although particular resistance genes may be linked to specific microbial communities correlations are most likely weak. A similar, but stronger trend was observed for the persistent resistome using Procrustes analysis when applying a goodness-of-fit test (sum of squares M2 = 0.5484, r = 0.672, p < 0.05, and 9999 permutations) using Bray–Curtis dissimilarity metrics. To reveal associations between bacterial species from the selected pathogen-associated families and specific ARGs, we performed correlation network analysis, in which the entire network consisted of 32 nodes and 28 edges (Fig. S7). The analysis identified a significant positive association between the commensal human-pathogen Enterococcus faecium (see above) and 6 ARGs including two aminoglycoside (aac(6′)-Ib(aka aacA4), aadA5), and tetracycline (tetG), vancomycin (vanC), phenicol (floR) and MLSB (mphB) resistance-conferring genes. The opportunistic pathogen Staphylocuccus saprophyticus (a notable portion of Staphylococcaceae) was linked to 5 ARGs including two vancomycin (vanSA & vanC1), two MLSB (vatE & ermB) and one tetracycline (tetK) resistance genes.

Discussion

The present study provided insights into potential microbial hazards of broiler litter applied to agricultural fields as fertilizer, and comprehensive understanding of spatial and temporal dynamics of microbiomes and resistomes in broiler litter stockpiles. BL stockpiling is extremely common throughout the world, and the study provided fundamental data that can be translated into best management practices that can be applied to meet local and international regulations and reduce epidemiological risks. BL stockpile age had a critical impact on microbiomes and resistomes, and bacterial diversity decreased over time, whereas ARG diversity increased. Although as a whole, stockpile depth did not affect ARG/MGE composition or diversity, the middle of the stockpile had considerably lower bacterial diversity relative to the top and bottom. This was most likely facilitated by high temperature [53], salinity [54] and a decrease in oxygen saturation [55]. Avidov et al. (2021) [10] demonstrated that the time required to achieve a 7 log10 reduction in Salmonella enterica serovar Infantis under controlled temperature, water-content, and initial pH, ranged between 13.7–27.2, 6.5–15.6, 1.2–4.7, and 1.3–1.5 days at 30, 40, 50, and 60 °C, respectively. This clearly demonstrates the capacity of temperature to reduce Gram-negative enteric bacteria. They proposed that antagonism by native communities is the primary cause of S. Infantis reduction observed under mesophilic temperatures (30 and 40 °C). This highlights the effect of both temperature and compitition in mitigating microbial pathogens in compost and static BL stockpiles.

Despite the decrease in the relative abundances of ARGs in the BL stockpile over time, more than half of the detected genes persisted (although sometimes at low abundances), and therefore may proliferate under favourable “gut-like” conditions [49], following ingestion of stockpile dust or fertilized produce. The pathogen-associated families Enterobacteriaceae, Enterococcaceae and Staphylococcaceae were profuse in the stockpile. Although the relative abundance of the first two decreased with time, they were detected in all time points, suggesting that broiler litter may be a vector of potentially antibiotic resistant pathogens under certain conditions [56], [57], [58], [59]. Although ARG diversity increased, the relative abundance of ARGs decreased with time. This may be linked to the temporal reduction in Enterococcaceae and especially Enterobacteriaceae abundance, given the fact that these two families harbor a large fraction of annotated ARG [44], [49]. In contrast to Enterobacteriaceae and Enterococcaceae, the relative abundance of Bacillaceae significantly increased with time, and was significantly higher in the middle of the BL stockpile. This may be attributed to the fact that they are spore formers and are more adept at surviving in natural environments than enteric bacteria [60], [61].

In microbiome studies, long-read sequencing technology has increased our capacity to achieve species-level characterization of 16S rRNA gene amplicons [62], [63]. Hence, we exploited the Loop Genomics full-length 16S rRNA gene sequencing platform [34] to target pathogen-associated families (described above), and identify priority pathogens in the broiler litter stockpile. A large fraction of Enterobacteriaceae were E. coli, which were detected in all depths of the stockpile samples. This is consistent with previous studies suggesting the long-time survival of E. coli in poultry manure and/or amended soils [64], [65]. We also detected presence of S. enterica and K. pneumoniae (at very low relative abundances), which are associated with food poisoning [66], [67]. The occurrence of different Salmonella species including S. enterica in poultry litter was previously documented [68], and these strains often harbor class 1 and 2 integrons that can capture and integrate antimicrobial resistance genes [69], [70], [71].

Taxonomic extrapolation also revealed a diverse array of enterococci in the BL stockpile including commensal human pathogens (E. faecium and E. faecalis) and avian strains (E. durans and E. galinarium), which generally decreased with time. Clinically relevant enterococci are often associated with vancomycin resistance [72], but HT-qPCR analysis revealed low abundance of vancomycin resistance genes in general, and vanA, which is frequently associated with multidrug resistance in opportunistic enterococci pathogens was not detected, suggesting that many of these strains are not resistant to vancomycin. This finding was recently supported by our comprehensive metagenomic assessment of chicken litter following prophylactic and therapeutic antibiotic treatment, which revealed that while certain vancomycin resistance genes (i.e. vanC) were detected in non-pathogenic Enterococcaceae strains, vanA was not detected and MAGs (metagenome-assembled genomes) of opportunistic human pathogens (i.e. E. faecium and E. faecalis lacked vancomycin resistance genes [9]. Interestingly, the relative abundance of the commensal human pathogen E. faecalis in the stockpile slightly increased with time. This is of concern, considering a recent study that detected E. faecium and E. faecalis in soil after litter amendment. Considering this data, we recommend targeting E. faecalis (especially vancomycin-resistant strains) in routine monitoring and source tracking experiments that evaluate transmission between BL, the environment and humans, through exposure to contaminated soil or consumption of poorly washed crops from manure fertilized farms [73]. While most of the identified staphylococci are not pathogenic, low abundances of the opportunistic pathogen S. aureus (an etiological agent of arthritis and gangrenous dermatitis in chickens and foodborne gastroenteritis in humans) were detected in the BL stockpile, especially at the later sampling times. This is supported by previous studies that also detected S. aureus in poultry litter, [74], [75]. Interestingly, HT-qPCR analysis detected persistent prevalence albeit low abundances of the methicillin resistance gene mecA, which is associated with methicillin resistance S. aureus (MRSA). Previous reports have demonstrated the prevalence and frequency of these difficult to treat pathogens in poultry litter, whose presence is linked to antibiotic usage [76], [77].

Correlation analysis indicated that most of the enterococci including the priority pathogens E. faecium and E. faecalis detected in the BL stockpile showed positive correlation with oxygen levels. This suggests that despite the fact that enterococci are generally considered to be facultative anaerobes, they persisted better in the stockpile under aerobic conditions. Similar phenomena have been observed in anaerobic digestion [20] and/or composting [19] of chicken litter, and hence reducing oxygen levels in chicken litter stockpiles may assist in mitigating antibiotic resistant pathogens in litter stockpiles in the future. This is supported by controlled simulations, that found higher reduction of S. Infantis strains under anaerobic conditions at mesophilic temperatures (<45 °C) [10].

Although many of the ARGs observed in this study are ubiquitous in the environment, several genes conferring resistance to clinically important antibiotics, e.g. β-lactams (i.e. ampC, blaCMY2, blaOXA1/OXA30, blaOXA10, blaPSE, blaZ, blaTEM, mecA, blaCTX-M, blaSHV) were observed at low relative abundances. Although it is difficult to evaluate the risk of these genes, they can potentially be horizontally transferred to soil and crop-associated microbiomes under favorable conditions. Co-occurrence analysis identified statistically strong associations between class D β-lactamases such as blaPSE (known to be found in S. enterica) and the class 2 integron-integrase gene, intI2 [78], and of the extended spectrum beta-lactamase blaOXA10 with the sulfonamide resistance gene sul2 [79], implying potential dissemination of these genes between environments via horizontal gene transfer. Several class 1 and 2 integrons have been identified in different Salmonella serotypes including S. enterica isolated from varied sampling sources including food animals [71]. Although abundance of S. enterica in the stockpile was relatively low, its presence suggests it may constitute an epidemiological risk under certain conditions in edible crops, considering the capacity of S. enterica strains to persist in the field [80]. As shown by Avidov et al [14], S. Infantis persisted over 90 days at 30 °C in a Vertisol soil amended with inoculated BL, presumably due to reduced antagonistic activity compared to the BL alone.

The transposases gene tnpA was strongly associated with certain ARG types including aadA, a aminoglycoside resistance gene that is often associated with various pathogens including E. faecalis and E. coli [81], [82]. A recent study showed increased abundance of aadA gene during 10 weeks of poultry manure composting [83]. Interestingly aadA was a hub gene in our network analysis, and it co-occurred with several different ARGs. This can facilitate co-selection of resistance genes, which can maintain high prevalence of AMR in the absence of antibiotics [84].

It is important to stress that in contrast to cultivation, the molecular methods (16S rRNA amplicon sequencing and HT-qPCR) applied in this study cannot determine if the detected bacteria and ARGs are viable due to the fact that DNA of lysed cells can bind to organic material and clay [85] in the BL stockpile. In contrast, Salmonella and other Gram-negative enteric bacteria can enter a viable but not culturable (VBNC) state in certain environments, and thus may not be detected using culture-based assays [86]. Given the limitations of both methods, we strongly recommend that future evaluation of BL microbiomes complement molecular analyses with cultivation-based methods that target potentially hazardous antibiotic resistant bacteria in BL, as recently demonstrated in our study that evaluated the effect of antibiotics on BL microbiomes and resistomes [9].

Conclusions

Collectively, we found that broiler litter may be a reservoir of opportunistic bacterial pathogens, including E. faecalis, S. enterica, K. pneumonia, S. aureus and S. saprophyticus, which are potentially associated with various ARGs. While stockpiling partially reduces bacterial pathogens and associated ARGs, there is a significant amount of spatial heterogeneity, with lower removal efficiency at the top and bottom of the stockpile, and therefore some of these antibiotic resistant opportunistic pathogens seemingly persist in the pile. Under particular scenarios these bacteria may be transmitted and survive in fertilized produce, posing a threat to public health. These results necessitate rethinking of stockpile regulations to reduce epidemiological risks. This should include requiring periodical screening for selected bacterial pathogens (and specifically for antibiotic resistant strains) described above using standard cultivation methods, mixing to increase litter exposed to the more optimal “inner-pile” conditions, and minimal BL stockpile maturation time of at least two months prior to BL amendment in the field.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.