Bovine respiratory syncytial virus and bovine coronavirus antibodies in bulk tank milk – risk factors and spatial analysis

Bovine respiratory syncytial virus (BRSV) and bovine coronavirus (BCoV) are considered widespread among cattle in Norway and worldwide. This cross-sectional study was conducted based on antibody-ELISA of bulk tank milk (BTM) from 1347 herds in two neighboring counties in western Norway. The study aims were to determine the seroprevalence at herd level, to evaluate risk factors for BRSV and BCoV seropositivity, and to assess how these factors were associated with the spatial distribution of positive herds. The overall prevalence of BRSV and BCoV positive herds in the region was 46.2% and 72.2%, respectively. Isopleth maps of the prevalence risk distribution showed large differences in prevalence risk across the study area, with the highest prevalence in the northern region. Common risk factors of importance for both viruses were herd size, geographic location, and proximity to neighbors. Seropositivity for one virus was associated with increased odds of seropositivity for the other virus. Purchase of livestock was an additional risk factor for BCoV seropositivity, included in the model as in-degree, which was defined as the number of incoming movements from individual herds, through animal purchase, over a period of five years. Local dependence and the contribution of risk factors to this effect were assessed using the residuals from two logistic regression models for each virus. One model contained only the x- and y- coordinates as predictors, the other had all significant predictors included. Spatial clusters of high values of residuals were detected using the normal model of the spatial scan statistic and visualized on maps. Adjusting for the risk factors in the final models had different impact on the spatial clusters for the two viruses: For BRSV the number of clusters was reduced from six to four, for BCoV the number of clusters remained the same, however the log-likelihood ratios changed notably. This indicates that geographical differences in proximity to neighbors, herd size and animal movements explain some of the spatial clusters of BRSV- and BCoV seropositivity, but far from all. The remaining local dependence in the residuals show that the antibody status of one herd is influenced by the antibody status of its neighbors, indicating the importance of indirect transmission and that increased biosecurity routines might be an important preventive strategy.


Introduction
The overall health among Norwegian dairy cattle is good with few endemic infectious diseases present. Several infections, such as bovine tuberculosis, bovine brucellosis and bovine viral diarrhea (BVD), have been eliminated through successful control programs (Sviland et al., 2015a(Sviland et al., , 2015bÅkerstedt et al., 2015). However, bovine respiratory syncytial virus (BRSV) and bovine coronavirus (BCoV) are endemic and prevalent in the national herd (Gulliksen et al., 2009). The prevalence of these infections is considered high in most parts of the world, and they cause disease problems leading to reduced animal welfare, increased use of antibiotics and financial loss for the farmer (Valarcher and Taylor, 2007;Boileau and Kapil, 2010;Sacco et al., 2014). BRSV causes respiratory disease, most often in young animals, and bronchopneumonia due to secondary bacterial infection is common (Larsen, 2000). BRSV was the most commonly isolated agent in respiratory outbreaks in cattle herds in a recent Norwegian study (Klem et al., 2014a). BCoV is the cause of calf diarrhea, respiratory disease and winter dysentery (contagious diarrhea in adult cattle) (Boileau and Kapil, 2010)  have shown significant effects of BCoV infection on production in terms of decreased milk yield and poor growth rate (Tråvén et al., 2001;Beaudeau et al., 2010b) which both result in economic loss.
Bulk tank milk (BTM) serology is a cheap and effective method used to screen herds for infectious diseases. However, due to long lasting seropositivity after infection, a herd will stay test-positive for many years after circulation of virus in the herd (Alenius et al., 1991;Tråvén et al., 2001;Klem et al., 2014b). Likewise, testnegative herds might have been virus free for years and serology on bulk tank milk is therefore an indicator of herd status with an inherent time-lag.
Herd level risk factors previously found to be of importance for BCoV status in Swedish dairy herds are herd size, not providing boots for visitors and geographic location (Tråvén et al., 1999;Ohlson et al., 2010). For BRSV seropositivity, herd level risk factors found to be of importance both in Scandinavia and beyond are herd size, age profile of the herd, type of production and existence of bordering cattle herds (Norström et al., 2000;Solís-Calderón et al., 2007;Ohlson et al., 2010;Saa et al., 2012).
Previous studies in Scandinavia have indicated large variations in prevalence of BRSV and BCoV between regions (Elvander, 1996;Tråvén et al., 1999;Beaudeau et al., 2010a;Klem et al., 2013), but spatial analyses involving BRSV and BCoV infections are infrequently reported. For control-and eradication purposes, locating high and low risk areas is important in order to know which control strategies should be applied to different regions. Risk factors like herd size, animal movement between herds, and proximity to neighbors are likely to vary geographically. However, it is currently not known how geographical differences in risk factors are associated with the spatial variation in prevalence of positive herds for these two viruses. Because the spatial pattern of antibody-positive herds may be largely driven by the spatial patterns of herd characteristics, such as herd size and distance to neighbors, spatial clusters of positive herds might only be reflecting the geographical distribution of known characteristics. Hence, it is of major interest to determine if adjusting for these factors changes the appearance of the spatial clusters.
BRSV and BCoV can be spread between herds by direct animal contact and indirect transmission. Direct contact includes physical contact between animals from different herds, for instance through shared pasture, or by live animal trade. Indirect transmission happens through passive transfer of animal secretions and excretions between herds by fomites like clothing or equipment.
The topography in western Norway, where the area under investigation is located, is characterized by mountains and fjords separating the herds and limiting direct contact. However, animal movements between holdings might provide an important route of transmission. In-degree is a measurement of animal contact which is defined as the number of incoming animal movements from individual herds, through animal purchase, over a defined time period (Nöremark and Widgren, 2014). Livestock movements are often registered in central databases which allows for calculation of in-degree, but factors affecting indirect transmission can be more difficult to assess because information on movement of people and biosecurity routines are not readily available in central registries. Nevertheless movement of people is associated with herd size, because larger herds have more visitors (Nöremark et al., 2013).
The aim of this study was to determine the spatial variation in herd-level prevalence of BRSV and BCoV, as measured by BTM-antibodies, across the study region in western Norway. Furthermore, the effect of the risk factors herd size, location, animal movement, and proximity to neighbors were evaluated and the effect of these risk factors on the spatial distribution of positive herds was assessed.

Study population
This cross-sectional study was performed in "Sogn og Fjordane" and "Møre og Romsdal" counties on the west coast of Norway (Fig. 1). The region was thought to be a suitable study area because of an expected mix of BTM-positive and negative herds. One BTM sample from each of 1347 herds was collected by the dairy company (Tine, Norwegian Dairies SA), between December 2012 and June 2013. In 2013, 1854 herds delivered milk in the two counties, which means samples were collected from 73% of all eligible dairy herds. Milk samples were treated and analyzed as described in Section 2.2, and each herd was categorized as either positive or negative based on the BRSV and BCoV antibody test results, respectively. If a herd contributed more than one sample during the study period, only the result from the first sample was included. Prevalence estimates were calculated for the region as a whole and for each county separately. True prevalence was calculated using the Rogan-Gladen-estimator based on the sensitivity and specificity of the tests as specified by the manufacturer (Greiner and Gardner, 2000).
During the study period, 98% of all dairy herds were members of the Norwegian Dairy Herd Recording System (NDHRS) which provides reliable records on herd characteristics, production parameters and disease occurrence (Espetvedt et al., 2013). The medical company distributing the only registered BCoV vaccine in Norway was contacted to get information regarding the number of units sold. The use of the only registered BRSV vaccine was recorded by contacting the veterinary practitioners by phone. Veterinarians in all municipalities of the study area with more than 10 herds were contacted, covering 1295 of 1347 herds.

Laboratory analysis/outcome variable
BTM samples were collected by the milk truck drivers and transported at a temperature of 4 • C to the dairy plant were they were frozen at between −18 and −20 • C, and were kept at this temperature until thawing at the time of laboratory analysis (Tine mastittlaboratoriet in Molde). All BTM samples were analyzed using indirect ELISA (SVANOVIR ® BRSV-ab and BCoV-ab, Svanova Biotech, Uppsala, Sweden). The optical density (OD) reading of 450 nm was corrected by the subtraction of OD for the negative control antigen, and per cent positivity (PP-value) (Takiuchi et al., 2009) was calculated as (corrected OD/positive control corrected OD) × 100. The cut-off for a positive result was set at a PP-value of 10 for both tests (Anon., 2016a,b) and the dichotomized results (BRSV +−/and/or BCoV +−) were used as the outcome in all analyses. The sensitivity and specificity of the ELISAs provided by the manufacturer were 94% and 100% for BRSV, and 84.6% and 100% for BCoV respectively (Alenius et al., 1991;Elvander et al., 1995).

Explanatory variables
Test results were combined with production data and health recordings from the NDHRS. All statistical analyses were done using Stata/SE (StataCorp. 2011. Stata Statistical Software: Release 12. College Station, TX: StataCorp LP) unless otherwise specified. Permission to use the database was given by the owner, Tine, Norwegian Dairies SA. All recordings were from the year of 2012. To describe the general characteristics of the herds in the area the mean, standard deviation and range were calculated for the following herd parameters retrieved from tables of annual summary data in the NDHRS: herd size, milk production, somatic cell count (SCC), overall herd disease incidence, and replacement rate. Herd size was defined as the herds' mean number of cow-years in 2012 (one cowyear = 365 days for a cow in a herd, calculated for each cow from date of first calving). SCC was measured as mean somatic cell count in BMT and milk production was measured as kg milk produced per cow-year. Herd disease incidence was the combined incidence rate for all recorded diseases per 100 cow-years in 2012, where recorded diseases include all cases treated by a veterinarian as this is reported in on-herd health recordings. Replacement rate was the number of cows in first lactation divided by the herds' number of cow-years, multiplied by one hundred. Reports of respiratory disease was available at the individual level, and this information was dichotomized to whether or not the herd had one or more animals with reported respiratory disease during the year of 2012. Herds that were not NDHRS members, or had incomplete registrations during this time, had to be excluded from the risk factor analysis, but were still included in prevalence estimates, point maps and isopleth maps. For an overview of eligible, sampled and analyzed herds see Fig. 2.
Data on animal movements between holdings were provided by the Norwegian Food Safety Authority, and in-degree was calculated as a measure of animal purchase as described in Section 2.4. Access to recordings on the location of each herd, given by geographic coordinates (latitude, and longitude, projection: EPSG: 4326-WGS 84), was provided by Tine, Norwegian Dairies SA. No information on the location of non-dairy cattle holdings was available. As a measure of proximity to neighbors, the mean Euclidian distance to the five closest dairy herds was calculated. This calculation also included herds outside the study area to avoid biased values for herds close to the county borders. The date at which the sample was collected by the tank milk driver, was divided into two categories: winter; December 1st-March 31th vs. summer; April 1st-June 30th.

Animal movements
In this study, the term 'animal movement' refers to change in ownership of an animal. Registration of cattle purchases is mandatory in Norway. In-degree was used as a measure of livestock movement, and is the number of direct ingoing contacts, from individual herds, through animal purchase (Nöremark and Widgren, 2014). In-degree was calculated as a sum of purchases from individual herds for a period of almost five years; January 1st 2008-December 5th 2012. I.e. an in-degree of five for a given herd in this study indicates that the herd has purchased live animals from five different holdings during the five year period described. Purchases reported after December 5th were excluded from the in-degree calculations because this was the date of collection of the first BTM sample. A total of 347 holdings had not registered purchases during this time period. This was assumed to be true and no herds were excluded due to missing values, but herds with inconsistent duplicate registrations were omitted (Fig. 2).

Univariable analyses
A total of 1194 out of the 1347 sampled herds, or 89%, had complete records, and were included in the risk factor analysis. To assess the probability of selection bias the proportion of positive herds was calculated for the 153 herds (11%) that did not have complete NDHRS records as well as all sampled herds. Herds lacking geographic coordinates were excluded from the maps and spatial analyses (Fig. 2).
Univariable analyses for a set of 11 predictors were performed in order to select which variables to include in the multivariable models. These predictors were chosen from the available data based on a causal diagram and biological plausibility of an association with the dichotomized test result of BRSV and BCoV antibodies in BTM. The same variables were evaluated for both viruses. Continuous variables assessed for an effect on the outcome were: herd size, herd disease incidence, average milk production per cow-year, replacement rate, mean SCC in BTM, geographic coordinates and average distance to the five nearest herds. The continuous variables were included as such in the analyses unless otherwise mentioned. Dichotomous variables were: time of sampling (winter; December 1st-March 31th vs. summer; April 1st-June 30th) and whether or not the herd had reported respiratory disease the year before sampling. The association between in-degree and the outcome was assessed treating in-degree both as a continuous-and as a categorical variable. For analytical and interpretational reasons in-degree was eventually included as a categorical variable with three categories: category 1 for 0-1 direct ingoing contacts, category 2 for 2-9 direct ingoing contacts and category 3 for more than 9 direct ingoing contacts.
For all variables the association with the outcome was evaluated by simple logistic regression (Wald-test), and the predictor was included in the subsequent model-building process if the pvalue < 0.2.
Linearity of continuous predictors was assessed by grouping observations in groups of equal size, and making plots of the group means against the log odds of the outcome. In case of non-linearity, different transformations were evaluated. To avoid multicollinearity in the model, correlation coefficients between all pairs of two predictors were calculated before the multivariable analysis was performed (Dohoo et al., 2003).

Multivariable analyses
Based upon the significant associations from the univariable analyses, two logistic regression models with different outcomes were built: one with the BRSV antibody status of the herd as the outcome and the second with BCoV antibody status as the outcome. Large scale trends, also called first-order spatial effects, relate to variation in the mean value of a spatial process (Dohoo et al., 2003), and to control for possible first-order effects the x-coordinate (longitude) and y-coordinate (latitude) were added in the model as continuous variables. Biologically plausible pairwise interactions between significant variables from the final models were assessed by adding their cross-product in the model and then determining if the coefficient for the term was statistically significant. For interactions, a more stringent criterion was used for model inclusion (p < 0.02) in order to choose the most parsimonious model. Possible confounding factors were identified through a causal diagram and monitored by calculating the changes in other covariates when one factor was added and withdrawn from the model. The final models were fitted using a manual backward stepwise procedure, with a selection threshold of p < 0.05. The area under the curve (AUC) of the receiver operating characteristic (ROC) was used to evaluate overall model performance, and the Hosmer-Lemeshow test was used as a test for the modelís goodness of fit, with data grouped in ten groups on the basis of percentiles of estimated probability.
Pearson and standardized deviance residuals were calculated for both models. To detect possible influential observations, Q-Q plots of Pearson residuals were made, and the delta beta statistics were calculated. Observations with high residual values or delta beta value above 0.2 were omitted, and the analyses were rerun to evaluate their impact on the estimates.
2.6. Spatial patterns 2.6.1. Point maps and maps of the prevalence risk distribution All maps were created using QGIS 2.4.0 (QGIS Development Team, 2014). Point maps were created to show the point location of all study herds with respect to their antibody status for the two viruses. Kernel density estimation was used for both BRSV and BCoV positive herds in addition to all herds, using the isotropic Gaussian kernel function implemented in the "spatstat" library in R (Baddeley and Turner, 2005). Kernel density estimation is a weighted moving average method that can be used to estimate the intensity, or mean function, for point processes (Berke, 2005). The resulting values can be presented as a raster map with one density value for each grid cell. A common fixed bandwidth determined from the coordinate ranges from the study herds ((1/8) × min(x range ,y range )) was used. Dividing the range distance by eight was done to avoid over smoothing the intensity function, as reported elsewhere (Vanderstichel et al., 2015). Generating a risk map with spatial point data (locations of cases and non-cases) is based on the ratio of two intensities as described by Berke (2005). Thus, the isopleth map showing prevalence risk on a smoothed color scale was made by dividing the Kernel density raster layer for the cases by the Kernel density raster layer for the background population.

Local clusters
The spatial scan statistic test was applied to explore spatial clusters of positive herds by using the software SaTScan version 8.1.1 (Kulldorff, 2009). The spatial scan statistic can analyze spatial point data (Kulldorff and Nagarwalla, 1995;Kulldorff, 1997), and cluster detection is done by gradually scanning a window across space, noting the number of observed and expected observations inside the window at each location (Kulldorff, 2015). Clusters of positive herds were detected using the Bernoulli model with analysis settings as purely spatial, scanning for areas with high rates and maximum spatial cluster size 20% of population at risk. No overlap of clusters was allowed. Results from the analyses includes location of clusters, the value of observed/expected cases, the relative prevalence (not shown) and a p-value for each cluster obtained by the Monte Carlo method (999 iterations).
To evaluate spatial clusters of positive herds first after correcting for first order effects and then after adjusting for other herd level risk factors, two sets of logistic regression models were built using BRSV-and BCoV-test results (0/1) as the outcome. One set included only the x-(longitude) and y-(latitude) coordinates (called the xy-models). The other set also included the predictors of interest that remained in the model as described in section 2.5.2. After model diagnostics and evaluation of model fit, the deviance residuals from all four models were obtained and analyzed using the spatial scan test under the normal probability model with analysis settings as purely spatial, scanning for areas with high values of residuals, and maximum spatial cluster size 20% of population at risk. No overlap of clusters was allowed. The output reports key statistics, including the location, the number of herds, the loglikelihood ratio and a p-value for each cluster obtained by the Table 1 Mean, standard deviation (STD) and range for herd parameters obtained from NDHRS from the year 2012 in 1194 dairy herds, included in a study of BRSV and BCoV as measured by bulk tank milk antibodies in the study area of Sogn og Fjordane and Møre og Romsdal county on the northwest coast of Norway. Monte Carlo method (999 iterations). Because the same modelling approach was used for the spatial assessment of the residuals from both the xy-and the final model it was possible to compare the spatial clusters of residuals before and after correcting for the risk factors.

Study population
Mean values, standard deviations and ranges of descriptive parameters for the study population are presented in Table 1. The overall apparent prevalence of seropositive herds in the study area was 46.2% for BRSV, and 72.2% for BCoV. 40.7% of all herds were positive for both viruses and 22.3% were negative for both viruses on BTM. This means that a herd which is antibody positive for one virus had a 5.3 times increased odds of positivity for the other virus. The prevalence of positive herds was higher in the northern county (Møre og Romsdal), 54.4% for BRSV and 79.8% for BCoV, compared to the southern county (Sogn og Fjordane), where the prevalence of BTM positive herds was 36.7% for BRSV and 63.4% for BCoV.
Based on the sensitivity and specificity of the ELISA tests given by the manufacturer, the calculated true prevalence was 49.1% for BRSV and 85.3% for BCoV. For the 153 herds (11%) that provided milk samples but were not part of the multivariable analyses (see Fig. 2), the prevalence was 50.3% and 77.1% for BRSV and BCoV, respectively.
Vaccination against BRSV was known to have been used in a total of six herds before the time of sampling. Five of these where in the northern county (Møre og Romsdal) and one in the southern county (Sogn og Fjordane). It was decided not to exclude any herds due to vaccination because vaccination was so rarely used and because the herds that had used it reported a prolonged history of respiratory disease and were likely to be antibody positive on BTM sampling regardless of the use of vaccine. Regarding the BCoV vaccine, no units of the vaccine were sold to pharmacies in the study area during 2012. This is not a guarantee that it is not used, but strongly implies limited use, and thus the risk that use of vaccine would influence the results was considered negligible.

Animal movements
Incoming animal movements were registered from most parts of the country, but the majority of purchases were across short distances within the study region. For the 1194 herds that had complete records, the median in-degree over the period of almost five years (January 1st 2008-December 5th2012) was 2 -with a range of 0-181.

Multivariable model
Time of year for sampling was excluded from the model due to collinearity with the geographic coordinates (r = −0.79 for x and r = −0.81 for y). The x-and y-coordinates were also correlated (r = 0.71), which was expected due to the north-east slope of the coastline. Because no herds are located off-shore, an increase in y will tend to entail an increase in x. The stability of the models was tested by removing the coordinates one at a time, fitting the model with only the x-coordinate, only the y-coordinate and both. No substantial changes were observed in the estimates for the other covariates in the model, and it was decided to keep both coordinates despite the correlation in order to correct for large geographic trends (first order effects) so that any remaining geographic variation in the residuals could be attributed to local dependence. The distance to the five closest dairy herds showed lack of linearity with the log odds of the outcome, and was therefore log transformed (natural logarithm).

BRSV-model
Variables included in the final BRSV logistic regression model were: herd size, x-and y-coordinates and log of mean distance to the five closest dairy herds (in km). Results from the BRSV logistic regression model are shown in Table 2. The area under the ROC curve was 0.73, and the p-value of the Hosmer-Lemeshow goodness of fit test with ten groups was 0.91 indicating acceptable overall fit of the model. Calculation of the delta beta statistics revealed no obvious outliers and no observations had delta beta >0.2.

BCoV-model
Variables included in the final BCoV logistic regression model were: herd size, herd disease incidence, x-and y-coordinates, log Table 2 Estimated odds ratios with 95% CI and coefficients with standard errors, along with p-values based on a logistic regression model on factors associated with herd level BRSV-status as measured by antibodies in bulk tank milk in 1194 dairy herds in two counties on the west-coast of Norway.   (year 2012). ** The number of a herd's direct ingoing contacts through animal purchase from unique herds over a period for almost five years. Category 1 includes herds with in-degree 0-1, category 2 for in-degree 2-9 and category 3 for in-degree more than 9. of the distance to the five closest dairy herds and in-degree. After the introduction of in-degree the variables "replacement rate" and "reported respiratory disease" were no longer positively associated with BTM positivity. Results from the logistic regression model are shown in Table 3. The BCoV model had an area under the ROC curve of 0.81, and a Hosmer-Lemeshow goodness of fit test with ten groups gave a p-value of 0.63, indicating good overall fit of the model. Calculation of the delta beta statistic detected five possible influential observations (delta beta >0.2). However, omitting these did not substantially influence the model estimates. The model had lowest predictive ability for large BCoV negative herds, with a relatively short distance to the five nearest dairy herds, located in the northern county. These herds were BCoV-negative despite the high probability of a positive outcome predicted by the model.

Point maps and maps of prevalence risk distribution
The point location of all study herds are shown in Fig. 3. Kernel density estimation was used to make smoothed maps of the prevalence risk distribution for evaluation of large trends regard-ing spatial variation of positive herds for the two viruses. These maps show the density of positive herds over and above the density of the background population, and are shown in Figs. 4 and 5. The spatial distribution of risk is similar for the two viruses with the highest prevalence risk in the northwestern region, and the lowest prevalence risk in the south.

Local clusters
Application of the spatial scan test under the Bernoulli model identified five spatial clusters of BRSV-positive, and four of BCoVpositive herds (p < 0.05). The BRSV-positive clusters included from 15 to 182 herds and the ratio of observed/expected cases ranged from 1.91 to 2.17. For clusters of BCoV-positive herds the number of herds in a cluster ranged from 30 to 160 and the ratio of observed/expected cases ranged from 1.23-1.39.
The location of spatial clusters of high values of deviance residuals from the xy-models and the final models are shown in Fig. 6. Key statistics from the analyses are shown in Table 4. A spatial cluster of high values of residuals is an area with an excess of cases based on what is expected under the current model. For the xy-models cluster analysis using the spatial scan test identified several areas with Table 4 Key statistics from the cluster analyses of residuals from the logistic regression model with x-and y-coordinates as the only predictors, and the final logistic regression model with all risk factors included, for BRSV and BCoV antibodies in bulk tank milk in 1194 dairy herds in two counties on the northwest coast of Norway. "Mean inside" and "Mean outside" refers to the mean value of deviance residuals inside and outside the cluster, respectively.  high values of deviance residuals. These clusters consist of positive herds with a low probability of positivity predicted by the model, i.e. the herds were expected to be negative when correcting for large (first order) geographic trends. Clusters with a p-value > 0.05 were excluded. The cluster analyses of model residuals detected six spatial clusters of BRSV-positive, and four of BCoV-positive herds. The BRSV-positive clusters included from 10 to 180 herds, two clusters were located in Møre og Romsdal, one on the border between the two counties, and three were located in Sogn og Fjordane. For clusters of BCoV-positive herds the number of herds in a cluster ranged from 52 to 233, one cluster was located in Møre og Roms- dal and the other three in Sogn og Fjordane. The most northern cluster had approximately the same geographic location for both viruses, a peninsula in the northwest of Møre og Romsdal (Romsdalshalvøya). For the final models the deviance residuals were spatially clustered in four locations for both viruses (p < 0.05, Fig. 6). For the BRSV-model the number of clusters was reduced, but the changes in log-likelihood ratio of the remaining clusters were small. For BCoV the number of clusters remained the same, but there were substantial changes in the log-likelihood ratio, see Table 4.

Number of cases
The spatial scan statistic will search for clusters with high values of residuals, and what is considered high values is relative to the Fig. 6. Geographic map of the study area indicating the location of clusters of high values of deviance residuals from the BRSV (A and B) and BCoV (C and D) logistic regression models. Clusters from xy-models are shown in A and C, and spatial clusters of high values of deviance residuals after correcting for all the risk factors in the final logistic regression models are depicted in B and D. Analyses were performed for n = 1194 herds in the study area situated in the northwest part of Norway. Clusters were detected using the normal probability model of the spatial scan statistic, and all clusters have a p-value < 0.05. Clusters are sorted according to likelihood ratio, with the most likely cluster as number one (number displayed in the center).
rest. The reference values was different for the BRSV and the BCoV models because the BRSV models had higher values of residuals on average both for the xy-model and for the final model. This means that the evaluation of spatial clusters must be interpreted as clusters of unexplained variation in the outcome for that model, and comparison of the spatial clusters of BRSV positive herds and BCoV positive herds must be done with caution.

Discussion
The overall apparent prevalence of seropositive herds in the study area was 46.2% for BRSV and 72.2% for BCoV, which is low compared to reports worldwide (Paton et al., 1998;Uttenthal et al., 2000;Ohlson et al., 2010). This is also lower than estimates from previous studies in Norway using serologic methods (Gulliksen et al., 2009;Klem et al., 2013). The present study classified herds according to detection of antibodies measured in a milk sample taken from the BTM. This methodology generally increases the prevalence of a disease when compared to individual sampling of a group of young animals -the method used by the previous Norwegian studies. Hence, it makes the discrepancy between the present study and the previous ones even larger, and is most likely due to differences between study regions. The study region was selected as it was believed it would contain a mix of BTM positive and negative herds. The large variation in prevalence across regions is in agreement with a study performed by Klem et al. (2013).
For both models the odds of being BTM positive increased from south to north (latitude) and for BRSV from east to west (longitude). These large trends can be interpreted as first order effects, but because the time of sample collection was correlated with the geographic coordinates, and had to be omitted from the model, the observed geographic trends cannot with complete certainty be separated from a possible temporal effect. About 40% of the herds were positive against both BRSV and BCoV, and the odds of being positive for one virus were approximately five times larger if a herd was positive for the other virus. The large proportion herds with antibodies against both viruses was not surprising given known common risk factors.
Herd size was positively associated with seropositivity for both BRSV and BCoV. Increasing the herd size by one cow-year increased the estimated odds of being antibody positive by 5% for both viruses. This corresponds to a 72% and 84% increase in the odds of BTM positivity when increasing the herd size across the interquartile range for BRSV and BCoV, respectively. The association between herd size and BRSV and BCoV positivity is well documented (Tråvén et al., 1999;Norström et al., 2000;Solís-Calderón et al., 2007;Ohlson et al., 2010). The dairy production in Norway is typically organized in small units with a mean of 24.2 cow-years per herd in 2013 (Anon., 2015). Even though this is smaller than in most developed countries this association holds true. The reason for the association remains unclear; however, it may be linked to larger herds having more indirect contact for instance via visits from veterinarians, AI technicians, advisory personnel or others (Norström, 2001). Furthermore, herd size might be associated with differences in management, and larger herds might provide better conditions for intra-herd virus circulation.
It is interesting that in-degree was only a significant predictor in the BCoV model, and not for BRSV. For the study population, the majority of purchased livestock came from within the region, and the fairly low herd level prevalence of BRSV in this region could explain why the number of ingoing contacts (in-degree) was not associated with BRSV positivity. The most commonly purchased animals are calves and young-stock, and because the prevalence on calf level is lower than on herd level, the risk of buying a young animal with either current viral infection or antibodies might not be high enough to show an association with the outcome. The prevalence of BCoV is higher, and thus the risk of buying antibody positive or infected animals is also higher, and more likely to affect the BTM result. A biological explanation behind differences in the likelihood of direct transmission between the two viruses should also be considered. An important difference in the pathogenesis of the two viruses is that BCoV replicates both in cells in the respiratory tract and intestinal epithelial cells, leading to shedding of virus in nasal secretions as well as in feces (Boileau and Kapil, 2010). On the other hand, BRSV only replicates in cells of the respiratory tract (Valarcher and Taylor, 2007). However, several important aspects of the pathogenesis are common for the two viruses: Shedding of virus is highest in the acute stage of the infection and disease can vary from subclinical to severe (Larsen, 2000;Cho et al., 2001;Boileau and Kapil, 2010).
Increasing mean distance to the 5 nearest dairy herds was associated with a significant decrease in odds of BTM positivity for both viruses. Association between existence of bordering herds and BRSV was also found by Saa et al. (2012). Another study found that the odds of BCoV positivity decreased as the distance to the nearest cattle herd increased, but no association was found for BRSV (Ohlson et al., 2010). Norström (2001) found an increased risk of outbreak of BRSV if at least one positive herd was within a radius of 500 m of a herd. BRSV and BCoV are enveloped viruses, with relatively short survival time outside the host, depending on environmental factors such as temperature, humidity and light (Hall et al., 1980;Larsen, 2000;Wolff et al., 2005;Casanova et al., 2010). As the number of infective virions on equipment decreases over time, the likelihood of indirect transmission by fomites decreases with increasing travelling time and therefore distance. Distance to neighbors will also influence on the number of possible indirect contacts for a herd, and thus the likelihood of exposure. This effect might be more evident during periods of high infectious pressure (epidemics). It is also possible that the distance between herds is associated with the risk of direct transmission, if animals in herd dense areas have more contact during pasture time in the summer. However, we did not have any information on the location of pastures.
The results of this study show that the geographic distribution of BRSV and BCoV in the study area are far from uniform, and that there are both local high risk clusters (Fig. 6), and large geographic trends (Figs. 4 and 5). The cluster analyses on the residuals showed that some of the local dependence changed when correcting for other risk factors. In other words, local dependence seems to be partially explained by spatial variation in the distribution of risk factors included in the logistic regression model, such as proximity to neighbors, herd size and large geographic trends (xand y-coordinates). However, spatial clusters of high residual values from the final models indicates that there are still spatially dependent unmeasured risk factors. No information on biosecurity was available and good hygiene and husbandry practices could be an unmeasured preventive factor. A study by Ohlson et al. (2010) showed a preventive effect of using boot covers on BCoV positivity. Other non-measured potential risk factors that could be spatially dependent include the use of common grazing, and historical data on previous disease outbreaks. Both winter dysentery and respiratory disease typically occur as epidemics with years between (Boileau and Kapil, 2010). An epidemic spread of these viral infections might cause all, or the majority of, herds in an area to be antibody positive, and thus affect the spatial distribution of positive herds for years.
For BCoV the number of clusters remained the same after correcting for the risk factors in the final model, however with large changes in the likelihood ratio. These changes in the log-likelihood ratio mean that adjusting for geographic differences in herd size, proximity to neighbors and in-degree results in a more random distribution of the residuals. However, the effect was not uniform for all clusters, indicating that the effect of the risk factors might not be the same in all areas. Compared to the BRSV-model, the moderate values of observed/expected for the BCoV clusters from the Bernoulli model also support that local dependence might be more important for BRSV than for BCoV. The spatial clusters of BRSV antibody positive herds had high values of observed/expected from the Bernoulli model, and there were relatively small changes in log-likelihood ratio of the clusters between the xy-and the final model, which might indicate strong local dependence. This also agrees with the lower predictive ability of this model compared to the BCoV model (AUC values 0.73 and 0.81, respectively). For BRSV the results imply the existence of spatially dependent unmeasured risk factors and that each herd relies strongly on the status of its neighbors, thus indicating the importance of indirect transmission routes. The implementation of a high level of biosecurity could, therefore, be important to prevent virus introduction. The higher overall predictive ability of the BCoV model compared to the BRSV model means that despite a higher overall prevalence of BCoV it is easier to predict the serologic status of a herd, or to locate "high risk herds", for BCoV than for BRSV, based on the number of animals purchased and relatively constant factors like herd size, proximity to neighbors and location. The difficulty in finding strong associations between the investigated risk factors and BRSV positivity, and the strong local dependence, could mean that the spread of BRSV in this region has been of a more epidemic character, involving more stochasticity than what has been the case for BCoV.
Classification of herds in this study was based on a single BTM sample. The use of BTM serology cannot be relied on to give an updated picture of the infection status of a herd because animals shed antibodies for years after infection (Alenius et al., 1991;Tråvén et al., 2001;Klem et al., 2014b). The proportion of herds with ongoing or recent infection is therefore likely to be much lower than the prevalence of BTM-positive herds. Several other diagnostic options for classification of herds with respect to BRSV and BCoV status exists. Serology of individual animals, either using milk or blood samples, can give a more recent picture of the herds' infection history than BTM samples, depending on the age and number of animals sampled. (Ohlson et al., 2009;Klem et al., 2013). The ideal method for classifying herds with respect to detecting virus circulation would be to detect the virus; however, this is demanding on a larger scale (Klem et al., 2013).
The antibody ELISA tests used to classify herds as either positive or negative for either virus are imperfect. This means that there will be some misclassification of outcome, which could lead to an underestimation of prevalence. (Estimates of true prevalence are shown in Section 3.1). However, as previously mentioned, prevalence estimates based on bulk tank milk serology will be much higher than the proportion of herds that have circulating virus, so compared to this the inaccuracy introduced by an imperfect test is negligible. When evaluating the effect of the risk factors, misclassification of the outcome was considered non-differential because the performance of the tests was not believed to be associated with any of the risk factors. Hence, this is not likely to have influenced the results.
The internal validity of this study was deemed high as all dairy herds in the study area were equally likely to be sampled and included in the study. The high proportion of sampled herds (73% of all eligible herds) also minimizes the risk of severe selection bias that could have affected the validity of the study. However, the prevalence estimates in the 153 herds that were excluded from the multivariable analysis due to incomplete NDHRS registrations, were slightly higher than for the entire population of sampled herds, which could indicate differences in, for example, management. But because the excluded herds represented only 11% of the sampled herds, and the difference in prevalence was modest, the introduced bias is likely to be small. The unknown location and BRSV-/BCoV status of beef herds in the area could potentially bias the results if the proximity to neighbors variable is incorrectly specified. However, the authors believe it is unlikely that their geographic distribution differs substantially from the distribution of dairy herds. The x-and y-coordinates were included in the models to reduce spatial heterogeneity. However, spatial correlation structures in the data may be more complex than a simple latitudinal/longitudinal gradient. In case of overdispersion due to spatial autocorrelation, this could alter the effective sample size, leading to increased chance of Type I error. However, given the low p-values and the inclusion of the x-and y-coordinates in the models, it is unlikely that the significance or direction of the effect estimates would change. The results of this study are believed to be representative for the Norwegian dairy herd as a whole, because the management systems for dairy production are comparable across the country. The external validity is therefore deemed good, and the results might also be valid for other temperate areas of smallerscale dairy production.
The study demonstrates that the herd level prevalence of BRSV and BCoV as measured by antibodies in bulk tank milk varied considerably in the region investigated. Of all the herds, about 40% were positive for both viruses. Several herd level risk factors were of importance for both BCoV and BRSV, such as herd size, geographic location and distance to neighboring herds, and for BCoV also in-degree. Adjusting for these risk factors explains some of the spatial clusters of positive herds, but spatial clusters of unexplained variation in the outcome was also detected. The remaining local dependence indicates that the antibody status of one herd is influenced by the antibody status of its neighbors and that indirect transmission is likely to be important. This means that a joint effort in terms of implementing preventive measures in an area could be an effective way to lower the prevalence of these infections. Measures should involve caution when purchasing livestock, implementing a high level of biosecurity and increased awareness among farmers and other people travelling between herds in order to prevent between-herd transmission of virus.