Quantitative Microbial Risk Assessment for Spray Irrigation of Dairy Manure Based on an Empirical Fate and Transport Model.

BACKGROUND: Spray irrigation for land-applying livestock manure is increasing in the United States as farms become larger and economies of scale make manure irrigation affordable. Human health risks from exposure to zoonotic pathogens aerosolized during manure irrigation are not well understood.
OBJECTIVES: We aimed to a) estimate human health risks due to aerosolized zoonotic pathogens downwind of spray-irrigated dairy manure; and b) determine which factors (e.g., distance, weather conditions) have the greatest influence on risk estimates.
METHODS: We sampled downwind air concentrations of manure-borne fecal indicators and zoonotic pathogens during 21 full-scale dairy manure irrigation events at three farms. We fit these data to hierarchical empirical models and used model outputs in a quantitative microbial risk assessment (QMRA) to estimate risk [probability of acute gastrointestinal illness (AGI)] for individuals exposed to spray-irrigated dairy manure containing Campylobacter jejuni, enterohemorrhagic Escherichia coli (EHEC), or Salmonella spp.
RESULTS: Median risk estimates from Monte Carlo simulations ranged from 10-5 to 10-2 and decreased with distance from the source. Risk estimates for Salmonella or EHEC-related AGI were most sensitive to the assumed level of pathogen prevalence in dairy manure, while risk estimates for C. jejuni were not sensitive to any single variable. Airborne microbe concentrations were negatively associated with distance and positively associated with wind speed, both of which were retained in models as a significant predictor more often than relative humidity, solar irradiation, or temperature.
CONCLUSIONS: Our model-based estimates suggest that reducing pathogen prevalence and concentration in source manure would reduce the risk of AGI from exposure to manure irrigation, and that increasing the distance from irrigated manure (i.e., setbacks) and limiting irrigation to times of low wind speed may also reduce risk. https://doi.org/10.1289/EHP283.

Introduction

Land application of livestock manure by spray irrigation is becoming more common in the United States (Genskow and Larson 2016). In fact, state-specific guidance documents for manure irrigation are available for Colorado, Idaho, Iowa, Missouri, Kansas, Nebraska, New Hampshire, New Mexico, North Carolina, Pennsylvania, Utah, and Wisconsin (Genskow and Larson 2016).

For manure irrigation, liquid manure in storage lagoons is pumped through pipes to the field site, where it is applied to the land surface using conventional water irrigation equipment, such as traveling gun and center pivot (Figure 1). In contrast with conventional methods of manure application before planting and after crops are harvested, irrigated manure can be applied throughout the year, including applications onto growing crops.

Figure 1.

Photos of spray irrigation conducted using traveling gun (top) and center pivot (bottom) equipment. Photos taken by authors Mark Borchardt (traveling gun) and Rebecca Larson (center pivot).

The expressed benefits of manure irrigation tend to be maximized at large economies of scale (Genskow and Larson 2016), and the trend of intensification in U.S. agricultural production practices includes consolidation to larger farms. For example, the percentage of dairy cows housed on U.S. farms with 1,000 or more head of cattle increased from 17% to 49% between 1997 and 2012 (http://www.nass.usda.gov/). Thus, increasing use of spray irrigation for land application of livestock manure seems likely to continue.

However, spray irrigation presents unique public health challenges relative to other technologies for manure land application. Spray irrigation equipment can apply manure for prolonged periods of time and throughout the growing season, potentially resulting in more continuous and more frequent exposure of nearby residents to aerosolized contaminants. Zoonotic pathogens are chief among the contaminants of concern in livestock manure (U.S. EPA 2013). During spray irrigation, wind may transport pathogens from application areas, leading to human exposure via inhalation or contact with contaminated vectors, food, and fomites.

While exposure to zoonotic pathogens from spray-irrigated livestock manure is plausible, the actual exposure levels and corresponding human health risks are not well understood. Concentrations of human pathogens and fecal microorganisms have been measured in air samples collected downwind from spray irrigation using untreated municipal wastewater (Katzenelson and Teltch 1976, U.S. EPA 1980). The incidence of enteric infections was higher among the residents of 77 communities that used partially treated domestic wastewater for spray irrigation than in residents of 130 communities that did not (Katzenelson et al. 1976). However, exposures and risks for spray-irrigated municipal wastewater may differ from those for spray-irrigated livestock manure because municipal wastewater contains human pathogens not found in livestock manure and undergoes treatment that differs from livestock fecal material.

We are aware of only two peer reviewed studies that estimated health risks posed by exposure to zoonotic pathogens in land-applied dairy manure, including a study that generated risk estimates for center pivot irrigation of dairy manure based on an air dispersion model (Dungan 2014) and one that estimated risks based on an empirical fate and transport model for coliphage spiked into groundwater that was sprayed using conventional tanker equipment (Brooks et al. 2012). Both studies used quantitative microbial risk assessment (QMRA), but their risk estimates varied over 10 orders of magnitude, and neither included data from empirical measurements during full-scale dairy manure spray irrigation.

Our objective was to estimate the human health risk from spray irrigation of dairy manure for residents with dwellings near spray-irrigated fields. We measured air concentrations of zoonotic pathogens and bovine commensal microorganisms downwind of traveling gun and center pivot spray irrigation equipment during full-scale irrigation events. We then modeled these measurements as a function of distance using empirical hierarchical models to account for variation in the concentration vs. distance relationship among trials. Finally, we used concentrations predicted from these models as inputs to a QMRA to estimate risk.

Methods

Field Sampling and Irrigation Equipment

We collected air samples for manure-related microbes during 21 manure irrigation trials. These trials were conducted during the growing season (May–October) over 3 y (2012–2014) on three Wisconsin dairy farms (herd sizes of approximately 400, 660, and 3,500 head) during routine manure application. These three farms were a convenience sample of local farms willing to participate in our study. Two farms were sampled throughout the 3-y study period; the third farm was sampled only during 2014. Manure used for irrigation was pumped directly from storage lagoons and applied to agricultural fields at 94,000 to . The storage lagoons contained untreated feces and urine (two smallest farms), the liquid fraction of manure that had passed through an anaerobic digester (the largest farm), and water flushed from animal pens and milking parlors (all farms). The two smallest farms used traveling guns (two trials and 11 trials, respectively), which moved at approximately , and released manure between 900 and . The third and largest farm used center pivot (eight trials) to irrigate the liquid fraction of anaerobically digested manure, which was diluted with groundwater (approximately 1:2) during irrigation. The pivots were approximately long, typically contained 80 to 120 nozzles two meters or more above ground surface, pivoted at , and released manure between 2,700 and .

Our ideal plan was to collect air samples at ten downwind stations at each farm during each trial. Two stations were placed apart at each of five target distances: 30, 61, 105, 152, and (100, 200, 350, 500, and ). Distance was measured from the irrigation-wetted perimeter, and stations were placed downwind from irrigation equipment along a line perpendicular to the center pivot’s arm (from the approximate midpoint) or to the direction traveled by the gun. Each station typically consisted of one button sampler (SKC Inc. SKC 225-360) and two single-stage Andersen impactor samplers (SKC Inc. SKC 225-9611) mounted 1.5 to (depending on crop height) above ground surface. Button samplers were used in all 21 trials and were loaded with gelatin filters (SKC Inc. SKC 225-9551) for quantitative polymerase chain reaction (qPCR) analyses. Impactor samplers were used in the 12 trials conducted during 2014 ( for farms A, B, and C). These were loaded with Petri plates of MacConkey, Cefoperazone Vancomycin-Amphotericin B Agar (CVA), 4-methylumbelliferyl-β-D-galactopyranoside (MUGal) and Indoxyl-β-D-glucuronide (IBDG) [MUGal and IBDG (MI)], and m-Enterococcus media (i.e., four types of media per distance). Calibrated air pumps for button samplers were set to ; impactor samplers were operated at . Two additional stations located approximately upwind from the irrigation equipment were used in all trials to collect background air samples before irrigation began and upwind samples during irrigation. Field conditions often forced us to adapt our ideal plan, so the exact number of air samples collected varied by trial (see Table S1).

A portable weather station logged meteorological data during every irrigation trial. The weather station consisted of a combined anemometer and wind vane (Model 03002, R.M. Young Co.), a pyranometer (Model CS300, Apogee Instruments) for measuring solar irradiance, a combined temperature and humidity probe (Model CS215, Campbell Scientific), and a CR800 data logger (Campbell Scientific).

Air sampling was timed to maximize airborne microbe detection. For center pivots, air sampling began when the system was fully pressurized and the nozzle arm began pivoting. For traveling guns, air was sampled for the 1-h period when the gun was most directly upwind from the samplers. Button samplers operated for 60 min (traveling guns) or 90 min (center pivots). Impactor samplers operated for 20 min (the recommended maximum period), then the plate was exchanged for a duplicate, and another 20-min sample was collected.

Liquid manure was sampled (100 to ) during irrigation from center pivots by holding a collection vessel over a nozzle and from traveling guns from a pump valve. Manure samples were analyzed by both qPCR and culture methods for the same genetic targets and microbes as the air samples. The number of individual manure samples collected per trial varied based on convenience ( samples per trial). Occasionally, we conducted two trials on the same farm on the same day. In these instances, we collected manure samples for just one trial and used those data for both trials. This occurred twice in 2014. Furthermore, we failed to collect manure samples for a single trial in 2012. Thus, we have 18 trials with unique qPCR analyses of manure samples instead of 21 and 10 instead of 12 for culture-based analyses.

Laboratory Analyses

Button sampler gelatin filters were placed into AVL buffer (Qiagen) immediately following sampling and the remainder of the DNA extraction procedure completed the next day [QIAamp DNA Blood Midi Kit (Qiagen)]. Liquid manure samples were diluted 1:10 in AE buffer (Qiagen) followed by DNA extraction [QIAamp DNA Blood Mini Kit (Qiagen)]. qPCR was performed with the LightCycler 480 instrument (Roche Diagnostics) and companion LightCycler 480 Probes Master kit. Microbial targets included Campylobacter jejuni, enterohemorrhagic Escherichia coli (EHEC), Salmonella spp., bovine polyomavirus, bovine Bacteroidales-like CowM3, bovine Bacteroides, Giardia lamblia, and Cryptosporidium parvum. Two to four qPCR technical replicates were performed per target. qPCR target genes, references for primers and hydrolysis probes, and quality assurance parameters for standard curves are reported in Table S2. Standard curves were created by serially diluting gBlocks or Ultramers (Integrated DNA Technology) of the target sequence. Cycle of quantification values (Cq) were determined by the second derivative maximum method. No-template controls were performed for every master mix and extraction batch. Data were acceptable only if all no-template controls in an analysis batch were negative (i.e., no Cq value). All 257 air samples (195 downwind, 33 upwind, and 29 background samples) and 55 manure samples were checked for qPCR inhibition following the approach of Gibson et al. (2012) using lambda phage DNA (New England BioLabs) as the inhibition control. No samples were inhibited based on this assessment (data not shown).

Bacterial culture procedures for manure and Andersen impactor samples were conducted to quantify gram-negative bacteria, Salmonella spp., C. jejuni, commensal E. coli, and Enterococcus spp. (see Table S3). Colony identity was confirmed by qPCR, except blue colonies on MI agar were counted as E. coli following U.S. Environmental Protection Agency (EPA) Method 1604 (U.S. EPA 2002). For culturable C. jejuni in manure samples, it was not possible to perform qPCR confirmations for all colonies on each plate, so the detection frequency for C. jejuni in manure based on culture procedures may be an underestimate. Concentrations of C. jejuni in manure based on culture procedures were estimated by multiplying the total concentration of colonies on each plate by the fraction of colonies testing positive as C. jejuni during qPCR confirmations. The positive-hole correction for Andersen samplers (Macher 1989) was not applied because the corrected plate counts were not substantially different than uncorrected plate counts.

C. parvum and G. lamblia were analyzed by qPCR for manure samples from irrigation trials in 2012 and 2013; manure samples from 2013 were additionally analyzed for both pathogens by immunofluorescence (Merifluor® Cryptosporidium/Giardia, Meridian Bioscience, Inc.). Positive controls for immunofluorescence microscopy were obtained from Waterborne Inc. Cryptosporidum and Giardia were never observed in any manure samples from 2012 and 2013, and as none of the three study farms had storage lagoons containing calf manure (the age most susceptible to Cryptosporidum and Giardia infections), samples collected in 2014 were not analyzed for these two pathogens.

Empirical Modeling

We modeled downwind air concentrations of four commensal microorganisms, bovine Bacteroides, Bacteroidales-like CowM3, Enterococcus spp., and gram-negative microorganisms, using two empirical models. The first set of models predicted concentration based on distance only (hereafter designated distance models) for use in acute gastrointestinal illness (AGI) risk calculations. The second set of models (hereafter referred to as distance plus trial-level variable models) were used to determine which of the following factors were significant predictors of downwind air concentrations: distance, microbe concentration in source manure, wind speed, solar irradiance, relative humidity, and temperature. We evaluated minimum, median, maximum, and mean values of the four meteorological parameters and used the version of each variable that best fit the data in final models, as described below. The microbes in source manure corresponded to the microbes for which air concentrations were modeled. We were missing measurements of manure microbe concentrations based on qPCR for one trial. This trial was excluded during model building for distance plus trial-level variable models. We were never missing meteorological data, distance, or manure microbe concentrations based on culture methods.

Hierarchical Modeling

Both modeling objectives entailed fitting separate two-level hierarchical models for each of the four commensal microorganisms. Concentration measurements of microorganisms measured using qPCR were grouped within trials (). Concentration measurements of microorganisms measured using culture methods were grouped within subtrials (). Subtrials are the 20-min time periods within trials during which Andersen impactors were operated. There are fewer subtrials than might be expected (12 trials using ) because field conditions prevented us from conducting two subtrials for three trials. Data points collected at the same distance for each trial (or subtrial) were handled as distinct data points during modeling.

Our hierarchical modeling approach implies that the trials (or subtrials) were independent and the concentration measurements within trials were correlated to the extent that they reflected the same underlying experimental conditions (farm, day, time, weather conditions). For distance models, the level 1 equation is: where, for the concentration of a given microbe in air, represents downwind measurement for trial , is the random intercept for trial , is the random slope for trial , is distance, and are the random residuals. The level 2 equations for distance models are: where, the terms are fixed regression coefficients, and the terms are random regression coefficients. In other words, the random intercepts and slopes of Equation 1 are predicted by two types of coefficients: a fixed overall coefficient and a random coefficient that varies by trial. The final distance model, with Equations 2 and 3 substituted into Equation 1, is:

For distance plus trial-level variable models, the level 1 equation is the same as Equation 1. However, the level 2 equations explicitly include terms for trial-level predictors: where, the x terms represent the k trial-level predictors (up to a maximum of five variables, including microbe concentration in source manure, wind speed, solar irradiance, relative humidity, and temperature). Their corresponding terms are fixed regression coefficients. Like the final distance model, the final distance plus trial-level variable model is derived by substituting Equations 5 and 6 into Equation 1.

The formulations presented above represent the most complex formulations that could be used for a given response variable in this study. However, a forward model-building process was used to evaluate simplifications to model structures and avoid potential problems with multicollinearity. For each response variable, models were developed by adding predictors one at a time and using nested model comparisons to determine the statistical significance of each additional term. These comparisons were made using a chi-squared test on the difference in deviance () between nested models. Additional terms were retained if they were statistically significant at the 0.05 level. Nested linear model comparisons were made using the appropriate estimation technique (restricted maximum likelihood or maximum likelihood).

For distance models, the sequence of terms added was a) fixed intercept and random intercepts; b) fixed slopes (i.e., fixed distance effects); and c) random slopes. Conceptually, this sequence corresponds to specifying that measured air concentrations a) vary only by trial; b) vary by both trial and distance; and c) vary by trial and distance, with the distance relationship itself also varying by trial.

For distance plus trial-level variable models, the sequence of terms was determined by screening all potential predictors individually against a null model, again using nested model comparisons. The initial null model consisted of fixed and random intercepts, and the list of potential predictors consisted of distance, microbe manure concentrations, wind speed, solar irradiance, relative humidity, and temperature. After screening each predictor individually, we selected the predictor with the lowest p-value, added it to the null model, then screened all remaining predictors individually against this new null model (original null model plus most significant predictor). This process was repeated until there were no more predictors with . For each of the four meteorological variables, we evaluated each possible form of the variable (minimum, median, maximum, or mean) and selected the measure that provided the best fit to the data. Note that it was possible to drop distance as a predictor in both types of models (distance models and distance plus trial-level variable models).

Sequential Probit and Linear Modeling

As part of the hierarchical modeling approach described above, each response variable was also decomposed into two components: a) a binary value indicating whether an organism was detected or not detected downwind; and b) for observations in which the organism was detected, the natural logarithm of the air concentration. We used probit models to estimate the probability of detection and linear models to estimate concentration conditional on a detectable value (Gelman and Hill 2007). This approach alleviated the need to impute values for nondetects prior to model fitting while simultaneously allowing the use of the full data set for each response variable. In the QMRA (see below), we used the probit model to predict the distribution of probabilities of detection at each distance, which was transformed to a distribution of 0s and 1s using a Bernoulli random number generator in R (version 3.2.1; R Foundation for Statistical Computing), and multiplied these values by the distribution of concentrations predicted by the linear model to derive a distribution of concentrations (including an appropriate proportion of nondetect samples with a concentration of 0) at each distance.

All probit and linear models used to estimate concentration distributions were fit using the lme4 package in R (Bates et al. 2015). Probit models were fit using LaPlace approximation, and the final linear models were fit using restricted maximum likelihood. Preliminary hierarchical modeling using an equipment type indicator variable (center pivot or traveling gun) indicated that equipment type was not a significant predictor of air concentration (, chi-squared test on difference in deviance for nested models), so equipment type was not included in models (data not shown). Probit models were checked for accuracy to quantify how well they reproduced the data they were fit to; accuracy here is defined as the proportion of data points for which the model correctly predicts the occurrence of detections and nondetections. Linear models were checked graphically to verify linearity, equal variance of errors, and normality of errors (data not shown). For distance plus trial-level variable models, we derived standardized coefficients by fitting models to standardized data (we divided each variable by one standard deviation). This was done to facilitate comparisons of effect size among variables with different ranges and units of measure.

Quantitative Microbial Risk Assessment

Air concentrations of zoonotic pathogens that may be present in manure and cause AGI, specifically, C. jejuni, EHEC, Salmonella spp., G. lamblia, and C. parvum, were of primary interest for AGI risk assessment, but these pathogens were rarely or never detected in air samples (data not shown). Therefore, we also measured air concentrations of four potential surrogate microorganisms (bovine Bacteroides, gram-negative microorganisms, commensal E. coli, and Enterococcus spp.). We selected bovine Bacteroides and gram-negative bacteria as surrogates to represent pathogen fate during manure irrigation because the two represented the highest and lowest air concentrations, respectively.

Predicted air concentrations of the surrogate microbes were estimated at each distance of interest using distance models, with probit models used to determine the probability of detection and linear models used to determine concentrations conditional on detection, as described above. To make predictions based on the distance models, we used the random coefficients in a Monte Carlo simulation, as described below (see Table 1 for a summary of all simulation inputs for these models). To estimate air concentrations of the pathogens C. jejuni, EHEC, and Salmonella spp., which are common in U.S. dairy operations (USDA 2003; USDA 2011) from predicted air concentrations for the surrogate microbes, we used pathogen-to-surrogate ratios defined from the literature for EHEC and Salmonella spp. (Hutchison et al. 2004) and ratios based on culture data from our own study’s manure samples for C. jejuni (Table 1). Exposure was also considered at two levels of prevalence for each pathogen in dairy farm manure: a) 100%; and b) a pathogen-specific typical prevalence value in the United States, namely, 39%, 40%, and 90% for EHEC, Salmonella spp., and C. jejuni, respectively. The values for EHEC and Salmonella are based on published national data (USDA 2003; USDA 2011). The value for C. jejuni (Table 1) is based on the proportion of trials from our own study in which C. jejuni was detected in any manure sample using qPCR; it is consistent with published national data (USDA 2011). We define pathogen prevalence as the proportion of dairy farms on which a given pathogen can be found in any manure sample. The use of two prevalence values and two surrogates resulted in four AGI risk estimates for each pathogen, thus capturing some of the uncertainty in these estimates.

Table 1

Summary of two-dimensional Monte Carlo simulation inputs.

Description	Typea	Distribution	Value(s)b	Source
C. jejuni prevalence	V	Empirical discrete	Values={0, 1}	This studyc
			Probabilities={0.10, 0.90}
Enterohemorrhagic Escherichia coli (EHEC) prevalence	V	Empirical discrete	Values={0, 1}	USDA 2003
			Probabilities={0.61, 0.39}
Salmonella spp. prevalence	V	Empirical discrete	Values={0, 1}	USDA 2011
			Probabilities={0.60, 0.40}
Distance (feet)d	V	Uniform	Min=100	Specified
			Max=1,000
Bovine Bacteroides probit model intercepte	VU	Mixture of 20 normal distributions	Q1=1.54	This study
			Med=2.10
			Q3=2.65
Bovine Bacteroides probit model slope	U	Normal	Mean=−0.15	This study
			SD=0.07
Bovine Bacteroides linear model intercepte	VU	Mixture of 20 normal distributions	Q1=0.30	This study
			Med=1.32
			Q3=3.14
Bovine Bacteroides linear model slope	U	Normal	Mean=−0.37	This study
			SD=0.06
Ratio of C. jejuni to bovine Bacteroides	C	NAf	1×10^−4	This studyc
Ratio of EHEC to bovine Bacteroides	C	NAf	5×10^−5	Hutchison et al. 2004 and this study
Ratio of Salmonella spp. to bovine Bacteroides	C	NAf	5×10^−4	Hutchison et al. 2004 and this study
Gram negatives probit model interceptg	VU	Mixture of 21 normal distributions	Q1=−0.26	This study
			Med=0.52
			Q3=1.17
Gram negatives linear model interceptg	VU	Mixture of 20 normal distributions	Q1=−5.34	This study
			Med=−4.30
			Q3=−3.20
Gram negatives linear model slopeg	VU	Mixture of 20 normal distributions	Q1=−0.26	This study
			Med=−0.13
			Q3=−0.003
Ratio of C. jejuni to gram negatives	C	NAf	7×10^−3	This studyc
Ratio of EHEC to gram negatives	C	NAf	4×10^−3	Hutchison et al. 2004 and this study
Ratio of Salmonella spp. to gram negatives	C	NAf	4×10^−2	Hutchison et al. 2004 and this study
Age (years)	V	Mixture of 18 uniform distributions	Q1=20	U.S. Census Bureau 2015
			Med=39
			Q3=56
Exposure time (minutes)	V	Mixture of five age-dependent distributions, each age-dependent distribution is a mixture of nine uniform distributions	Q1=40	U.S. EPA 2011, Table 16–20
			Med=100
			Q3=189
Inhalation rate (cubic meters per minute)	V	Mixture of 14 age-dependent distributions, each age-dependent distribution is a mixture of 14 uniform distributions	Q1=1.1×10^−2	U.S. EPA 2011, Tables 6–17 and 6–19
			Med=1.2×10^−2
			Q3=1.4×10^−2
Ingestion-to-inhalation ratio	C	NAf	0.8	Hardy et al. 2006
C. jejuni dose–response model parameter (alpha)	U	Empirical continuous	Q1=1.3×10^−1	Schmidt et al. 2013
			Med=1.5×10^−1
			Q3=1.8×10^−1
C. jejuni dose–response model parameter (beta)	U	Empirical continuous	Q1=4×10^0	Schmidt et al. 2013
			Med=1×10^1
			Q3=4×10^1
C. jejuni morbidity ratio	C	NAf	0.28	Soller et al. 2010
EHEC dose-response model parameter (alpha)	U	Empirical continuous	Q1=1×10^−2	Teunis et al. 2008
			Med=3×10^−1
			Q3=8×10^0
EHEC dose–response model parameter (beta)	U	Empirical continuous	Q1=2×10^0	Teunis et al. 2008
			Med=4×10^1
			Q3=8×10^2
Salmonella spp. dose–response model parameter (alpha)	U	Empirical continuous	Q1=3×10^−2	Teunis et al. 2010
			Med=3×10^−1
			Q3=2×10^1
Salmonella spp. dose-response model parameter (beta)	U	Empirical continuous	Q1=2×10^−3	Teunis et al. 2010
			Med=1×10^−1
			Q3=6×10^1
Salmonella spp. dose–response model parameter (eta)	U	Empirical continuous	Q1=7×10^−2	Teunis et al. 2010
			Med=8×10^−1
			Q3=1×10^2
Salmonella spp. dose–response model parameter (rho)	U	Empirical continuous	Q1=1×10^0	Teunis et al. 2010
			Med=7×10^0
			Q3=3×10^1

Simulation inputs are either constant (C), variable (V), uncertain (U), or variable and uncertain (VU).

Constants and parameters of parametric distributions are presented directly. Summary statistics of simulated values are provided for continuous empirical distributions and mixture distributions. Distributions that are both variable and uncertain are summarized at their median in the uncertainty dimension.

Parameters for Campylobacter jejuni are calculated based on manure data from the trials used to construct hierarchical models plus manure data from an additional two trials conducted during the same study period. The air data for these two trials are not presented here because manure application was by tanker spraying, not spray irrigation. However, the manure data from these two trials are representative of C. jejuni prevalence in our study and allow for a larger sample size.

Distance is treated as a random input when performing sensitivity analyses and assessing output stability. It is not treated as random when determining the relationship between risk and distance.

The empirically observed correlation () between random coefficients of the bovine Bacteroides probit model and random coefficients of the bovine Bacteroides linear model was reproduced in the variability dimension of the simulation using the cornode function in mc2d (Pouillot and Delignette-Muller 2010).

Not applicable.

The empirically observed correlations among random coefficients of the gram-negative models were reproduced in the variability dimension of the simulation using the cornode function in mc2d (Pouillot and Delignette-Muller 2010). Spearman’s correlation coefficients among these random model coefficients were (between the linear intercepts and linear slopes), 0.38 (between the linear intercepts and probit intercepts), and (between the probit intercepts and linear slopes).

Inhalation volumes were calculated from age-dependent distributions of inhalation rates and time spent outdoors at home (US EPA 2011) (Table 1). The age distribution (25th, 50th, and 75th percentiles of 20, 39, and 56 y, respectively) was defined using 2010 U.S. Census data for Wisconsin (U.S. Census Bureau 2015). Inhalation rates were for light intensity activity (e.g., walking) (U.S. EPA 2011). The duration of exposure depends on the time spent outdoors, the duration of irrigation, and the degree to which these two time periods overlap. However, data are lacking for the latter two items, so we assumed that the duration of irrigation was equal to the amount of time spent outdoors and that the two time periods overlapped perfectly.

Pathogen dose was the product of pathogen air concentration, inhalation volume, and an ingestion-to-inhalation ratio of 0.8 (Hardy et al. 2006), which is high compared with ratios of 0.1 to 0.5 used in similar QMRAs (Brooks et al. 2005b, Tanner et al. 2008, Brooks et al. 2012, Dungan 2014).

Risk of AGI was estimated using beta-Poisson dose–response models for EHEC (Teunis et al. 2008), Salmonella spp. (Teunis et al. 2010), and C. jejuni (Schmidt et al. 2013) (Table 1). The C. jejuni dose–response model estimates probability of infection, which was multiplied by a morbidity ratio of 0.28 (Soller et al. 2010). Random samples of dose–response model parameters were obtained from authors of the three models. To avoid overestimating risk at low doses (Teunis and Havelaar 2000), exact beta-Poisson models were evaluated using the Kummer confluent hypergeometric function in the gsl package for R (Hankin 2006).

Computations

The QMRA was evaluated by two-dimensional Monte Carlo simulation using the mc2d package for R (Pouillot and Delignette-Muller 2010). Two-dimensional Monte Carlo simulations consist of multiple simulations addressing variability (i.e., the variability dimension) nested within a single larger simulation that addresses uncertainty (i.e., the uncertainty dimension). In this context, the term uncertainty is akin to measurement error; it refers to heterogeneity in QMRA inputs that could theoretically be reduced by collecting more measurements. In contrast, our use of the term variability refers to heterogeneity in QMRA inputs that is not reducible because it is due to natural processes. Heterogeneity in risk factors (e.g., age, exposure time) was assigned to the variability dimension, while heterogeneity in uncertain model parameters (e.g., dose–response parameters) was assigned to the uncertainty dimension (see Table 1 for the dimension associated with each individual input variable). Pathogen prevalence was not treated as a distribution within Monte Carlo simulations. For each pathogen, the pathogen was considered to either be always present (using 100% prevalence) or present during a fixed percentage of simulated exposure events (using typical prevalence). The distributions of random coefficients (e.g., probit and linear intercepts) from distance models were simulated in the variability dimension. Thus, these coefficients represent the aggregate effects of trial-level variables (i.e., risk factors like wind speed and microbe manure concentration) on risk estimates. To determine the relationship between risk and distance, simulations were performed in a loop at each of 19 distances (100 to in increments). To perform sensitivity analyses and assess output stability, simulations were performed with distance specified as a random variable (see Table 1 for distribution). Each simulation consisted of 15,000 iterations in the variability dimension and 3,000 iterations in the uncertainty dimension, which resulted in 45 million risk estimates for each of the four scenarios () for each pathogen. The overall size of each simulation was constrained by available computing power; the size of each dimension within simulations was optimized within this constraint based on the stability of simulation outputs (see Table S4).

Results

Meteorology

Our measurements of airborne microorganisms during manure irrigation were conducted across a wide range of meteorological conditions (Figure 2). For instance, mean wind speed varied from 0.9 to among irrigation trials, maximum wind speed varied from 2 to , and solar irradiance varied between 0 (complete darkness) and 1,000 (bright sunshine) .

Figure 2.

Mean temperature (T), mean relative humidity (RH), mean wind speed (WS), maximum wind speed (), and mean solar irradiance (SI) during all 21 field trials. Each symbol represents a field trial. Reprinted from Borchardt and Burch (2016), with permission.

Microbes in Source Manure

C. jejuni was the only zoonotic pathogen detected in manure samples (in 16 of 18 field trials). Commensal microorganisms (i.e., pathogen surrogates) were detected in all manure samples (see Table S5). The detection frequencies of microorganisms were consistently higher in manure samples than in air samples. For example, commensal E. coli were detected in 100% of manure samples, but only 11% of air samples. Likewise, the concentrations of microbes in manure were consistently higher than in air. For instance, median concentrations of bovine Bacteroides were over 100 million times higher in manure () than in air from the irrigated wetted perimeter ().

Microbes in Air

Detection frequencies and concentrations of airborne microbes tended to decrease with distance (Figure 3). Background (i.e., before irrigation) and upwind (i.e., during irrigation) concentrations of airborne microorganisms are reported in Tables S6 and S7. Random slopes did not significantly improve model fit in distance models for Bacteroides and Bacteroidales-like CowM3, so these models included fixed slopes only. Consequently, the estimated rate of decrease with distance for these organisms was constrained by the models to be equal for all trials. In contrast, random slopes were retained in the final models for Enterococcus spp. and gram-negative bacteria, and the rates at which predicted concentrations decreased with distance varied among trials. Random intercepts were retained in distance models for all four organisms; thus, air concentrations at the minimum distance of (i.e., the intercept) varied among trials. The distance models explained the relationship between distance and microbe air concentrations well. The accuracy of probit models varied between 0.85 and 0.89. For linear models, the values were 0.80, 0.76, 0.89, and 0.58 for gram-negative bacteria, bovine Bacteroides, Enterococcus spp., and Bacteroidales-like CowM3, respectively. We attribute the poor fit of the Bacteroidales-like CowM3 linear model to apparent outliers (see Figure S1).

Figure 3.

Air concentrations of microorganisms downwind of full-scale dairy manure irrigation. Panels on the left represent quantitative polymerase chain reaction (qPCR) data, while panels on the right represent culture data. Points represent the median concentration measured at each distance; error bars are the first and third quartiles. Percentages represent the detection frequency at each distance. Each point represents 11 to 42 measurements. Based on our ideal sampling plan, each point should represent 24 or 42 measurements for culture measurements and qPCR measurements, respectively ( for each measurement type). However, field conditions often forced us to deviate from this plan.

In the distance plus trial-level variable models, distance remained significantly associated with downwind air concentrations for each microorganism (). It was retained in linear models for Enterococcus spp. and gram-negative microorganisms; it was retained in both probit and linear models for bovine Bacteroides and Bacteroidales-like CowM3. In addition, four trial-level variables were also significantly associated with air concentrations (Table 2). Wind speed and microbe concentration in manure were positively associated with downwind concentration () for three and two of four microorganisms, respectively. Relative humidity and solar irradiance were negatively associated with downwind concentration () for two and one of four microorganisms, respectively.

Table 2

Summary of distance plus trial-level variable models.

Microorganism	Modeld	Standardized coefficienta,b,c
		Distance	Wind speede	Microbe manure concentration	Relative humidityf	Solar irradianceg
Bovine Bacteroides	Probit	−0.27±0.13	–	–	−0.63±0.19	−0.72±0.18
	Linear	−0.31±0.05	0.43±0.13	–	–	–
Bacteroidales-like CowM3	Probit	−0.60±0.14	–	–	–	–
	Linear	−0.18±0.08	0.54±0.12	–	−0.27±0.13	–
Enterococcus spp.	Probit	–	0.83±0.42	1.08±0.44	–	–
	Linear	−0.35±0.09	0.41±0.19	–	–	–
Gram negatives	Probit	–	–	0.67±0.23	–	–
	Linear	−0.34±0.11	–	0.37±0.14	–	–

Estimate of fixed error. Dashes indicate that the relevant predictor was not included in the final model due to a lack of statistical significance (at the 0.05 level) during the stepwise model building process.

We standardized coefficients by fitting models to standardized data. Data were standardized by dividing by one standard deviation for each variable.

In addition to the five predictors indicated in this table, temperature was also evaluated, but was not significant (at the 0.05 level) in any model.

The probit model predicted microorganism detection (Yes/No); the linear model predicted microorganism concentration conditional on detection.

Wind speed is represented as median wind speed for the bovine Bacteroides and Bacteroidales-like CowM3 linear models. It is represented as minimum wind speed for the Enterococcus spp. probit and linear models.

Relative humidity is represented as median relative humidity for the bovine Bacteroides probit model. It is represented as minimum relative humidity for the Bacteroidales-like CowM3 linear model.

Solar irradiance is represented as maximum solar irradiance for the bovine Bacteroides probit model.

Risk Estimates

Risk estimates for AGI decreased with distance and varied substantially based on pathogen prevalence and surrogate (Figure 4). Median risk estimates were derived by finding the median of risk distributions produced by the Monte Carlo simulations, and varied between roughly and per irrigation event. Risk estimates were highest when pathogen prevalence in manure was assumed to be 100% and when bovine Bacteroides was used as the pathogen surrogate. They were lowest when pathogen prevalence was assumed to equal typical values (C. jejuni: 90%, Salmonella: 40%, EHEC: 39%), and gram-negative bacteria were used as the pathogen surrogate. In terms of pathogens, Salmonella spp. presented the highest risk when pathogen prevalence was modeled as 100%, while C. jejuni presented the highest risk when pathogen prevalence was modeled using typical values.

Figure 4.

The estimated probability of acute gastrointestinal illness (AGI) plotted against distance for pathogens modeled using different combinations of pathogen prevalence and surrogate. Prevalence values are either typical or 100%. Typical prevalence values are 90% for Campylobacter jejuni (determined in this study), 40% for Salmonella (USDA 2011), and 39% for enterohemorrhagic Escherichia coli (EHEC) (USDA 2003). Surrogate microorganisms are either gram-negative bacteria (GN) or bovine Bacteroides (BB). Plotted values represent the median in the variability dimension at the median in the uncertainty dimension. The U.S. Environmental Protection Agency (EPA)’s acceptable risk levels for drinking water (1 in 10,000 per y; blue line; U.S. EPA 1989) and recreational water (32 in 1,000 per swimming event; green line; U.S. EPA 2012) are depicted for context because acceptable risk has not been established for manure irrigation. The acceptable risk level for drinking water is visible in three panels: the top right, bottom left, and bottom right. The acceptable risk level for recreational water is only visible in one panel: the bottom right. EHEC plots are masked by the lowest risk values for Salmonella spp. (typical prevalence) or C. jejuni (100% prevalence). All the models used to generate risk estimates allowed risk to vary by distance. Risk estimates that do not appear to vary by distance are either much lower than risk estimates for other pathogens on the same plot or are the result of low pathogen prevalence (i.e., the pathogen is not present, so , regardless of distance). Please note that vertical axis scales vary among panels.

Sensitivity Analysis on Risk Estimates

The relationship between risk estimates and simulation inputs depended on the assumed level of pathogen prevalence in dairy manure (Figure 5). For EHEC and Salmonella spp., both of which have relatively low typical pathogen prevalence values (39–40%), risk estimates were most correlated with whether the pathogens were simulated to be present in dairy manure, which is determined by assumed prevalence values. For C. jejuni (), however, risk estimates were roughly equally correlated with simulated pathogen presence in manure, distance, exposure time, and the model intercepts, which represent the aggregate effects of trial-level conditions (i.e., microbe manure concentration and meteorological conditions). Individual risk factors like age and inhalation rate were not strongly correlated with estimated risks, although the full risk assessment cannot explicitly account for the effects of age on risk. Age is included explicitly in the exposure estimate as a predictor of time spent outdoors and inhalation rate, but age-specific dose–response models are not available. Sensitivity analysis using gram-negative bacteria as the pathogen surrogate produced similar results to those for bovine Bacteroides (data not shown).

Figure 5.

Sensitivity analysis of risk estimates to Monte Carlo simulation inputs defined in the variability dimension for Campylobacter jejuni, enterohemorrhagic Escherichia coli (EHEC), and Salmonella spp. modeled using typical pathogen prevalence on dairy farms and bovine Bacteroides as a pathogen surrogate. Typical prevalence values are 90% for C. jejuni (determined in this study), 40% for Salmonella (USDA 2011), and 39% for EHEC (USDA 2003). Bars are the median of Spearman’s correlation coefficients between risk estimates and the inputs listed on the vertical axis. Error bars are the 2.5 and 97.5 percentiles of these correlation coefficients. Coefficients are calculated from simulations that included all of the inputs shown with distributions as defined in Table 1. Probit and linear model intercepts represent the aggregate effects of trial-level conditions (i.e., meteorological conditions, microbe manure concentrations).

Discussion

Factors Affecting Risk of Illness

We estimated the risk of AGI for people downwind of dairy manure spray irrigation for three pathogens common in U.S. dairy operations. Overall, risk decreased with distance and tended to be most strongly controlled by the assumed level of pathogen prevalence in manure. Exposure time and the combined effects of trial-level conditions (i.e., microbe manure concentration and weather conditions) were also important when pathogen prevalence was assumed to be 100%, but no single factor dominated risk in this case. Acceptable risk for manure irrigation has not been established, so comparison to EPA’s acceptable risk levels for recreational water and drinking water are used for context: median risk estimates were generally low relative to acceptable risk levels in the United States for recreational water (32 illnesses per 1,000 primary contact recreators per event; U.S. EPA 2012). Compared to the acceptable risk level for U.S. drinking water (1 infection per 10,000 people per year; U.S. EPA 1989), risk estimates were higher when pathogens were modeled using 100% prevalence and similar, or lower when pathogens were modeled using typical prevalence values.

Based on how often they were significant predictors of air concentrations for individual microorganisms (Table 2), distance downwind and wind speed were the most important variables associated with air concentrations. Even though a brown cloud of suspended manure is visible during manure irrigation (Figure 1), microorganism concentrations and detection frequencies decreased from source manure to air, and concentrations of airborne microorganisms decreased with distance downwind. The decrease in microbial detections and concentrations likely resulted from aerosolization, deposition, dispersion, and microbial decay (i.e., inactivation). Our experimental design precludes isolating the individual contribution of each mechanism. However, the reduction in microbe concentration from the original manure by a factor of (or more) over the first of airborne transport demonstrates the substantial effects of these mechanisms.

Comparison to Epidemiological Studies

To our knowledge, there are no published reports of outbreaks following exposure to spray irrigation of dairy manure, and we are aware of only one early epidemiological study of the association between illness and spray irrigation of partially treated, nondisinfected municipal wastewater (Katzenelson et al. 1976). Katzenelson et al. (1976) estimated the excess risk of shigellosis and salmonellosis in communities practicing spray irrigation of partially treated, nondisinfected municipal wastewater to be and , respectively, which is higher than our median risk estimates for Salmonella-related AGI, assuming typical Salmonella prevalence in manure () and similar to or lower than our estimates at 100% prevalence ( to ). However, these comparisons are tenuous, as our risk estimates assume exposure is via airborne pathogens, while Katzenelson et al. (1976) did not consider exposure routes.

Comparison to Other Quantitative Microbial Risk Assessments

Our risk estimates are consistent with or higher than those obtained in two other QMRAs that estimated health risks due to aerosolized dairy manure. Brooks et al. (2012) estimated that the risk of infection was higher for L. monocytogenes than C. jejuni, EHEC, or Salmonella spp., with estimated probabilities of infection from a single exposure downwind from the source between and . Dungan (2014) reported that the estimated risk of infection with C. jejuni was greater than the risk of infection with EHEC, Salmonella spp., Listeria monocytogenes, or non-O157 E. coli, with estimated probabilities of infection from a single exposure downwind from the source (the closest distance considered) between and . We estimated that the risk of AGI from exposure to irrigated dairy manure was higher for C. jejuni than for EHEC or Salmonella spp. when we assumed typical pathogen prevalence levels, with median probabilities of AGI between and for a single exposure from the irrigation source, and to at (the farthest distance we considered). Assumptions about the risk of infection following exposure to a given pathogen dose were similar for our study and the two previous QMRAs. Therefore, it is likely that differences in estimated risks are related to differences in the assumptions and methods used to estimate exposure.

Our exposure assessment was based on empirical observations of microbe air concentrations during full-scale dairy manure spray irrigation and differs from the two previous studies. Brooks et al. (2012) based their exposure assessment on an empirical model derived from data collected during tanker application of groundwater seeded with a viral pathogen surrogate (Brooks et al. 2005a). Dungan (2014) based his exposure assessment on assumed inputs to an air dispersion model. As a result, the pathogen doses estimated by these two previous studies (which can be back-calculated from their reported risk estimates and dose–response models) vary significantly from ours. Specifically, L. monocytogenes doses estimated by Brooks et al. (2012) at varied between and colony-forming units (CFU), while our median estimated C. jejuni doses, based on typical pathogen prevalence, varied between and CFU at roughly . C. jejuni doses estimated by Dungan (2014) at varied between and , while our median estimated C. jejuni doses based on typical pathogen prevalence varied between and CFU at . Thus, high doses compared to previous work likely account for our relatively high risk estimates, although the exact aspects of the exposure assessment driving differences in dose are unclear (e.g., sampling approach and conditions, choice and use of surrogate, models, assumed pathogen levels in source).

Our risk estimates are also high compared to QMRAs for aerosolized human fecal material. Brooks and others estimated risk for coxsackievirus A21 and Salmonella spp. (Brooks et al. 2005a and 2005b) aerosolized from land application of treated municipal biosolids. They estimated the single exposure probability of infection at to vary between and , depending on the pathogen and assumed pathogen content of the source biosolids, which is much lower than our estimates at a similar distance. These differences in risk may be due to methodological differences in exposure assessment (e.g., air sampling), or they may represent true differences in downwind pathogen exposure due to differences between treated municipal biosolids and the dairy manure in our study (e.g., pathogen levels, degree of aerosolization).

Considerations for Policy Decisions

Using these results to guide risk mitigation policies is a three-step process. First, acceptable risk must be clearly defined. Any definition of acceptable risk consists of both a risk threshold and the statistic used to summarize risk. We have summarized our risk estimates using medians, though other summary statistics are possible (e.g., 90th percentile) and would result in different levels of public health protection. Likewise, we compare our risk estimates to standards for water (standards do not exist for spray irrigation), but these standards may not be appropriate for manure irrigation. Second, pathogen prevalence values for manure must be selected. Using 100% pathogen prevalence implies that risk-management decisions are being made for the population of exposure events where zoonotic pathogens are always present in dairy manure. Alternatively, using typical pathogen prevalence implies that risk-management decisions are being made for the population of exposure events where pathogens may or may not actually be present. Third, the pathogen surrogate must be defined. In general, downwind concentrations and detection frequencies of bovine Bacteroides were higher than culturable Gram-negative bacteria in our study. Thus, AGI risk estimates were higher when measured air concentrations of bovine Bacteroides were used as the surrogate for pathogen concentrations. However, assumptions about the prevalence of each pathogen had a greater influence on risk estimates (variation by a factor of 100) than the commensal organism used as the surrogate (variation by a factor of 10).

Interpreting these results for policy decisions also depends on several important considerations. First, we estimated risks for single exposure events. However, spray irrigation allows multiple manure applications to a single field in one growing season, which increases potential exposure events. As a result, cumulative risk estimates for the entire season will be higher than those presented here. Second, we estimated risk for C. jejuni, EHEC, and Salmonella spp., but less prevalent pathogens like L. monocytogenes, C. parvum, and G. lamblia can be found in dairy manure (U.S. EPA 2013). While the low prevalence of such pathogens likely renders risk lower than for the high-prevalence pathogens that we considered, a risk of illness is possible. Third, estimated risks from exposure to spray-irrigated dairy manure may not reflect risks associated with exposure to spray-irrigated manure from other livestock, which may differ with regard to pathogen concentrations and specific organisms (e.g., hepatitis E virus in swine manure). Finally, our risk estimates only consider exposure via ingestion of aerosolized pathogens and do not consider other exposure routes (e.g., direct contact by workers, runoff, vectors) or other impacts of manure irrigation (e.g., odor).

Limitations

Three study limitations must be considered when interpreting our risk estimates. First, the 21 field trials represent conditions on only three farms. Selection bias could have occurred if the participating farms, or the trials included in the study, were not representative of farms or conditions during spray irrigation in general. Second, with regard to model assumptions, we assumed that there were no systematic differences among the farms, since we did not include an additional model level to nest trials within farms. Third, because the pathogens of interest were rarely detected in manure or air samples, we used concentrations of common commensal microorganisms as surrogate measures of pathogen concentrations, and applied pathogen-to-surrogate ratios to estimate pathogen risks based on the observed data for the surrogates.

Conclusions

Our median risk estimates, which varied between and , were generally lower than acceptable risk levels set by the EPA for AGI as a consequence of exposure to recreational water (U.S. EPA 2012), but were sometimes higher than acceptable risk levels for drinking water (U.S. EPA 1989). It is important to note that acceptable risk levels have not been established for AGI as a consequence of exposure to spray irrigation of manure. Risk estimates were most strongly controlled by pathogen prevalence in manure. When pathogen prevalence was assumed to be 100%, no single factor completely controlled risk. Distance, exposure time, and the aggregate effects of weather conditions and microbe manure concentration were equally important in this case. Our model-based estimates suggest that reducing the prevalence and concentrations of zoonotic pathogens in source manure, increasing distances from irrigation sources to points of potential human exposure, and restricting applications based on wind speed may reduce potential human health risks from spray irrigation of livestock manure.

Click here for additional data file.