Spread of porcine circovirus associated disease (PCVAD) in Ontario (Canada) swine herds: Part II. Matched case-control study

Background The emergence of porcine circovirus associated disease (PCVAD) was associated with high mortality in swine populations worldwide. Studies performed in different regions identified spatial, temporal, and spatio-temporal trends as factors contributing to patterns of the disease spread. Patterns consistent with spatial trend and spatio-temporal clustering were already identified in this dataset. On the basis of these results, we have further investigated the nature of local spread in this report. The primary objective of this study was to evaluate risk factors for incidence cases of reported PCVAD. Results A time-matched case-control study was used as a study design approach, and conditional logistic regression as the analytical method. The main exposure of interest was local spread, which was defined as an unidentified mechanism of PCVAD spread between premises located within 3 kilometers of the Euclidean distance. Various modifications of variables indicative of local spread were also evaluated. The dataset contained 278 swine herds from Ontario originally sampled either from diagnostic laboratory submissions or directly from the target population. A PCVAD case was defined on the basis of the producer's recall. Existence of apparent local spread over the entire study period was confirmed (OR = 2.26, 95% CI: 1.06, 4.83), and was further identified to be time-varying in nature - herds experiencing outbreaks in the later part of the epidemic were more likely than control herds to be exposed to neighboring herds experiencing recent PCVAD outbreaks. More importantly, the pattern of local spread was driven by concurrent occurrence of PCVAD on premises under the same ownership (OREXACTwithin ownership = 25.6, 95% CI: 3.4, +inf; OREXACToutside ownership = 1.3, 95% CI: 0.45, 3.3). Other significant factors included PRRSv status of a herd (OREXACT = 1.9, 95% CI: 1.0, 3.9), after adjusting for geographical location by including the binary effect of the easting coordinate (Easting > 600 km = 1; OREXACT = 1.8, 95% CI: 0.5, 5.6). Conclusions These results preclude any conclusion regarding the existence of a mechanism of local spread through airborne transmission or indirectly through contaminated fomites or vectors, as simultaneous emergence of PCVAD could also be a result of concurrent change in contributing factors due to other mechanisms within ownerships.


Background
Porcine circovirus associated disease (PCVAD), also known as porcine circovirus disease and previously as post weaning multisystemic wasting syndrome (PMWS), emerged in the early 1990's and soon became a major animal health problem in many swine-producing regions worldwide [1,2]. The disease has several clinical presentations and can result in high mortality in growing pigs. Most frequently, PCVAD clinically manifests as excessive and rapid weight loss in growing pigs, respiratory illness, and increased mortality. Post-mortem findings frequently include generalized lymphadenopathy on macroscopic examination, and lymphoid depletion and/ or histiocytic replacement of follicles in lymphoid tissues on microscopic examination [3,4]. Porcine circovirus type 2 (PCV2) has been recognized as a necessary cause of PCVAD, but several other infectious and non-infectious factors have been identified as component causes acting through various mechanisms [4,5]. Studies of epidemics at the herd-level within distinct geographical regions [6][7][8][9][10] are disproportionately rare, relative to the impact that PCVAD had on production, possibly due to diagnostic uncertainty and availability of geographical and temporal information required to make inferences about the spread. Reports describing large-scale extensive investigations of PCVAD involving trace-back and trace-out are limited [6,10], likely reflecting the epidemiology of PCV2 infection, and the non-reportable nature of the disease. Management changes [11,12] and vaccination have been used to control the disease. In North America, commercial vaccines first became available during 2006, and since then became one of the most commonly used vaccines in growing pigs.
Studies conducted in different regions before vaccine introduction identified spatial, temporal, and spatio-temporal trends as factors contributing to patterns of the disease spread [8,9], although not under all conditions. In our accompanying article [13], we have identified spatial trends in PCVAD risk and described the nature of spatio-temporal dependence. This dependence was interpreted as the existence of a pattern of local spread. On the basis of this exploratory analysis, we have postulated further hypotheses about local spread and investigated them here. The primary objective of this study was to evaluate risk factors for incidence cases of reported PCVAD. The main exposure of interest was local spread, which in this manuscript was defined as the unidentified mechanism of PCVAD spread from infectious to susceptible premises located within 3 km of the Euclidean distance. Table 1 contains descriptive statistics for a subset of variables that were included in an inferential analysis based on the time of PCVAD outbreaks (in case herds) and on time of selection as a time-matched control (in the applicable risk set). We will refer to the sampling times as herds, although they reflect the status of a herd at a particular time (month) since December of 2003. The high proportion of PRRSv-positive herds in each PCVAD status is a reflection of the sampling strategy that favored selection of PRRSv-positive herds. Of note is that the proportion of herds that had "high risk local exposure" (HRLE) at the time of sampling was relatively low for case and control herds. Aside from being a characteristic of the study population, this might be a characteristic of the target population as well, because pork producers in Ontario generally try to maintain distance between swine premises for biosecurity purposes. In addition, it is noteworthy that control herds did not have exposure to PCVAD-infectious herds within ownership, unlike the case herds. In the case group, a higher proportion of local exposure was from PCVAD-infectious premises that were in different ownership than from PCVAD-infectious herds within the same ownership (Table 1, 4.7% versus 2.9%, respectively). Table 2 details measures of univariable associations obtained from conditional logistic regression based on assumption of Gaussian distribution (GD) for the test statistic. Additionally, measures of univariable associations obtained from the exact logistic regression were obtained for nominal variables ( Table 2). The exact method provided more reasonable point estimates and test of significance once ownership similarity between case status and HRLE was taken into account. This risk factor appeared to follow a "dose-response" pattern because the odds increased as the number of infectious herds in the 3-km zone increased ( Table 2, LNHRH). Local spread seemed to be driven primarily by high-risk local exposure to infectious herds within the ownership ( Table 2; OWHRLE). In addition to variables that met the criteria for inclusion into multivariable analysis, Table 2 also contains some non-significant associations of interest, such as cumulative local exposure (CLE), number of direct contacts/month, number of indirect contacts, and number of breeding contacts/month for farms where the breeding part of the swine operation was present. Although not shown in Table 2, non-significant associations were also detected for frequency of direct contacts/month through semen, gilts and boars, growing pig shipments, veterinary visits, veterinary supplies visits, and rodent control visits. A variable depicting frequency of feed delivery (Table 2, number of feed contacts/month) was not statistically significant (P > 0.10) once herds with frequency of feed contacts of >25/ month were excluded from the analysis (n = 4, n strata = 1). Similarly, a candidate final model that contained frequency of visits by the feed suppliers also indicated that the association was not statistically significant once this influential point was removed from analysis. On this basis, a model containing frequency of feed contacts was not considered a reasonable candidate for the final model. Interactions between local high-risk exposure  Table 3 lists the final multivariable model fitted using the likelihood method and the exact method. The exact logistic regression contained only the variables linked with local spread as exposure of primary interest, for which coefficients were provided, and herd PRRSv status and binary variable of easting (>600 km East) as confounders. A final exact logistic regression model including easting as a continuous explanatory variable could not be fit.

Discussion
In the accompanying article [13], spatio-temporal clustering at the small spatial and temporal scale was detected, suggesting that local spread was a possible pattern of PCVAD spread. Specifically, using the space-time K-function, we detected considerable spatio-temporal clustering of cases up to approximately 4.5 km of spatial and 5 mo of temporal distance to a PCVAD case (proportional increase of >100% of what was expected), with the highest increase at 2 km of spatial and 1 mo of temporal distance. In contrast to the previous method, the approach taken here required that we define local exposure using fixed spatial and temporal distance, recognizing that this could lead to some misclassification of exposure. We aimed to use cutpoints that would provide a balance between specificity of defining exposures (higher if shorter distance used) and number of herds that would be classified as exposed (higher if longer distance used). Thus, three km of spatial distance was selected subjectively and used as the only analytical approach for geographical distance, while a temporal distance of 3 mo was used as the most important, but not the only, temporal cutpoint. This spatial cutpoint was in concordance with how local spread was defined in a field during the 2001 outbreak of foot-and-mouth disease in the UK [14]. Our assumption was that the causative agent of PCVAD from the initially infected premises would not be more contagious to neighboring farms than the foot-and-mouth disease virus. The limitation of our original finding of the existence of apparent local spread was that the results were not adjusted for factors that could facilitate such a conclusion, even in the complete absence of radial spread of a causative agent from an infected farm to a neighboring susceptible herd. One such factor is ownership of the premises. In at least some Ontario regions, premises under the same ownership cluster geographically, and finding a pattern of spatial or spatio-temporal clustering of infectious disease might be facilitated by introducing infected animals from the same source, concurrent change in management, or frequency of other component causes, or simply by recall bias of the interviewee. It was therefore of interest to explore the nature of local spread after adjusting for ownership, and a time-matched case-control study was considered a suitable design for this objective.
Of interest for this investigation was finding that CLE was not associated with the emergence of incident PCVAD cases. Thus, even if local spread existed, the farms that reported a "historical" PCVAD problem did not contribute to the emergence of incident PCVAD cases in the neighborhood, at least in this study. On the other hand, the existence of spatio-temporal clustering at the small geographical and temporal scale -described in the accompanying article [13] -was confirmed in this study by having a significant variable that was indicative of HRLE. The results suggested that PCVAD-positive herds were more likely to have neighboring herds (up to 3 km) involved in a recent PCVAD outbreak (0 to 2 mo prior, total = 3 mo), than the PCVAD-negative herds. It is also noteworthy that this local exposure showed some limited dose-response relationship, with odds of being a case higher for herds that had two infectious herds within a 3-km distance than for herds with no infectious herds or only one such herd. Interaction of HRLE with time was biologically plausible and in agreement with previously reported results of PCVAD spread. Woodbine et al reported that radial spread was important in the later period of PCVAD spread in Great Britain [8,9]. When considered as the only term in the model, local spread has a time-varying nature (statistical interaction with time). The variable indicative of local spread (HRLE) became statistically significant in October of 2005, corresponding approximately to the 50 th percentile of reported PCVAD outbreak times, and remained significant until the end of the study.
However, when HRLE was partitioned into withinownership and outside-ownership, it was obvious that the likelihood of local spread was determined primarily by the co-existence of neighboring herds under the same ownership that experienced a PCVAD outbreak at concurrent times. The point estimate for the odds ratio of "High-risk local exposure within ownership" in the final model was very high and with very wide confidence intervals. This certainly is a limitation of this analysis and is a consequence of absence of this exposure in the control herds, as visible from descriptive statistics. However, collectively these results suggest that it is possible that such a mechanism of "local spread" could be driven by mechanisms that do not involve any type of spatial interaction between the neighboring farms (e.g., introduction from a common source, simultaneous change in management or prevalence of other component causes, or recall bias). We have no information that would allow us to further investigate these hypotheses. Our results, however, suggest that exploratory spatial and spatio-temporal analysis should be substantiated by additional analyses that will take ownership or pig movement into consideration, particularly in zones where spatial clustering of premises under the same ownership exists. We defined ownership at the highest possible level, and this might have some impact on the results. Large systems might have had several independent pig flows, and, theoretically, neighboring premises might belong to different flows (have different pig sources). In such situations, the nature of high-risk exposure would likely be incorrectly classified as withinownership. We believe that the latter scenario is relatively infrequent and would likely have a low impact on overall results and general interpretation. This issue of ownership, however, does warrant access to suitable ownership and movement databases that would allow epidemiologists to adjust their analysis for the proper ownership and the contact structure. In addition it is important to note that the interaction of HRLE with time was not significant once the ownership was accounted for. Although the interaction of HRLE with time was a biologically plausible term to include in the final model, the OWHRLE was selected instead because we considered these results were more insightful in explaining the underlying biological processes. The final model also contained two additional confounders: herd PRRS status and the linear effect of easting. PRRS infection in individual pigs has been identified as a factor contributing to PCVAD due to its effect on the immune system, and as a risk factor for PCVAD expression in observational studies [15][16][17].
Also of interest was the finding that the frequency of direct contacts, including through breeding animals for a subset of suitable herds, was not identified as a risk factor. This was a surprising finding because introduction of infectious agents through live animals is considered an important component of PCVAD emergence [8,9,12]. Two factors might have contributed to this. First, the frequency of visits used as a covariate in this analysis does not necessarily correspond to the time prior to a PCVAD outbreak in a herd. It is possible that the frequency of direct contacts changed in response to the PCVAD outbreak or an outbreak of another disease (e.g., PRRS). While short-term change in frequency of movements can be expected as a part of outbreak control, a long-term change is much less likely because it would require a change in the type of farming system, typically associated with high cost. The second possible explanation is that emergence of PCVAD is linked with the specific PCV2 source, and not with the general measure of frequency of incoming visits or movements. This is of note because either PCV2 as a necessary cause, or PRRSv or other agents as component causes, could spread through movement of animals and have an impact on disease emergence.
Although frequency of feed deliveries was identified as the only significant factor related to visits, this was most likely a chance finding, on the basis of the following rationale. The significance of this association was influenced by one herd that reported a very high frequency of feed deliveries, but the feed was delivered from the farm's own feed mill. This would contrast with the logical explanation that a causative agent was distributed either through feed from an external point source, or indirectly through delivery trucks between farms. Feed has been discussed as a possible source previously [12]. Evaluation of PCV2 transmission through feed containing contaminated porcine plasma yielded opposing results under experimental conditions [18,19].
Case definition was the largest limitation of this study, and the rationale for including it in the analysis is given in the accompanying article [13]. Briefly, preliminary analysis with more specific case definitions gave us similar results. Additionally, the evaluation of epidemic spread under current conditions using complete set of diagnostic criteria is not possible because of high vaccination usage and endemic presence of PCV-2 infection in many herds. Introduction of commercial PCV2 vaccines likely influenced the development of this outbreak. A major supplier gradually introduced its product to Canada in April of 2006 [20,21]. We fitted the final model only with cases that reported PCVAD prior to April of 2006 and found results (not reported) similar to our final model. We thus believe that results reported here were not influenced in a substantial manner by vaccination.

Conclusions
In conclusion, the spread of PCVAD due to high-risk local exposure was primarily driven by herd ownership. It was therefore impossible to distinguish local spread from common direct or indirect sources that contributed to the emergence of disease concurrently in different premises under the same ownership. Surprisingly, frequency of direct or indirect contacts did not differ between case and control herds. This also applies to frequency of feed delivery, which, although statistically significant, was driven by a single highly influential point, which in itself did not align with a biological rationale for spread through feed deliveries. Two risk factors that remained stable were herd PRRSv status and directional spread in a western direction. Multivariable analysis following exploratory spatial analyses proved to be a useful complement for this study. After completion of this analysis, our conclusions about disease spread have changed. Accounting for ownership and contact structure will become increasingly important in swine populations because investigation of spatial and temporal information alone would have yielded incomplete results and different recommendations.

Dataset
The full description of the data sources and data management steps is provided in the accompanying article [13]. Briefly, 278 herds located in southern Ontario were included in this study initially using two different sampling mechanisms: (i) sampling of herds positive for porcine reproductive and respiratory syndrome virus (PRRSv) using data from a diagnostic laboratory (Animal Health Laboratory, University of Guelph, Guelph, ON, Canada), based on rtPCR-positive diagnostic findings, and (ii) sampling of PRRSv-negative herds, based on assessment of a herd veterinarian. Herds were classified into PCVAD-positive and negative herds on the basis of producers' assessments, and results obtained during exploratory spatial analysis were used to guide the design of this case-control study.

Exposures of interest Local exposure
Estimation of the spatio-temporal K-function suggested significant spatio-temporal clustering over small spatial and temporal distance [13]. The risk was proportionally higher (relative increase >1) for up to~4.5 km of spatial and up to~5 mo of temporal distance from an arbitrary herd experiencing a PCVAD outbreak, and was increasing as both distances (in Euclidean distance and in time) to an arbitrary selected PCVAD-positive herd decreased. This result led us to generate a binary variable depicting "high-risk local exposure" (HRLE) to a PCVAD herd, whereby a herd would be considered positive for HRLE if it was within 3 km of a spatial distance from a herd that started experiencing a PCVAD outbreak that month or up to 2 mo earlier (total of 3 mo). Maximum relevant spatial and temporal distance was partly based on the spatio-temporal K-function, while considering that a shorter distance to the nearest herd could increase data scarcity and complicate analysis. A binary variable depicting HRLE to PCVAD was generated from available spatial and temporal data in the following way. In the first step, the herd-level dataset was expanded to include 52 mo duration of the period for each herd. This time period corresponds to the duration of time between December of 2003 (Month = 0), after which date producers started reporting outbreaks, and April of 2008 (Month = 52), when the last producer was interviewed. Each PCVAD-positive herd was considered infectious in the month when PCVAD was first reported and in the two subsequent months for a total duration of 3 mo, as non-infectious otherwise, and as censored after the month of the interview. In the second step, data for each month were exported to GIS software (ArcGIS 9.1) and buffering for each infectious herd was performed. Each buffer had a radius of 3 km and all herds that were contained within each buffer were selected and considered exposed to a high-risk PCVAD herd. In the third step, only herds exposed to another high-risk PCVAD herd were classified as exposed. This was done by removing herds that were used to generate buffers (infectious herds) using basic geoprocessing (PBSMapping; R 2.8.0) and data processing operations (R 2.8.0). Thus, a herd was considered positive for HRLE in a particular month of the study period only if it had another infectious PCVAD herd (herd was infectious over a period of 3 mo) within a 3km distance. This operation of generating binary variable HRLE was repeated for each month, and the resulting dataset served to generate a new set of variables that investigated different aspects of local exposure.
The second variable of local exposure was generated by counting the number of herds that experienced an outbreak of PCVAD in the month of interest or the preceding 2 mo and were located within a 3-km radius around a herd of interest. We refer to this variable as "local number of high-risk herds" (LNHRH). The third variable of local exposure was generated by taking ownership of the exposing (infectious) herd and the exposed herd into account. A list of ownerships and contact names was available for each herd, where ownership was considered a membership in a particular production system at the highest possible level. If ownership was unknown, we used a contact name as a proxy for ownership. This list was used in order to depict whether an exposing (infectious) herd was under the same ownership as a herd that was exposed. This resulted in generation of a nominal variable, "ownership-adjusted highrisk local exposure" (OWHRLE) with three levels: (i) no high-risk local exposure (OWHRLEno), (ii) high-risk local exposure within ownership (OWHRLEin), and (iii) high-risk local exposure outside of ownership (OWHR-LEout). In addition, a binary variable that depicted exposure of a herd to another herd that had experienced PCVAD at any point in the past (not only during the high-risk period) and that was ≤3 km distant was also created during data processing. We refer to this variable as "cumulative local exposure" (CLE).

Demographic variables
Other variables considered for analysis were already present in the dataset and were considered as time-invariant and explored in the accompanying article [13]: (i) PRRSv status, (ii) herd type, (iii) capacity for growing pigs on premises (nursery and finisher capacity). Additionally, monthly incoming contact rates from outside sources were calculated for: (i) frequency of monthly contacts through animal movement -direct contacts -and number of direct sources (animal movements of nursery and finisher pigs to premises), (ii) frequency of indirect monthly contacts (through feed shipments, veterinary visits, veterinary supplier visits, and service visits), (iii) number of indirect monthly contacts through each of the indirect source types, (iv) number of monthly contacts through breeding stock supply (gilts and boars) for sites that had breeding (sows) present on-site. An assumption was made while making this calculation that contact rates and other time-variant variables were unchanged since the time of the PCVAD outbreak. In cases where a variable had missing values, we imputed the values using routines available in commercial software (Stata 10), and using all other variables with available information. In the final dataset, values of four observations from two different herds were imputed (out of 680 observations in the matched dataset).

Matching
Incidence density sampling was used for matching. Strata were constructed so that each PCVAD case herd was individually time-matched to three randomly selected control herds from the available risk set on the month of PCVAD occurrence in a case herd. The available risk set in a specific month were all herds that: (i) were not PCVAD-positive (infectious) in this month or earlier (current or previous case herds), (ii) were not interviewed prior to this month, and (iii) were not already matched as controls to another case in that specific month (simple random sampling without replacement). Herds could have been, however, selected as controls to another PCVAD-case in subsequent months. Exposure status with respect to HRLE, LNHRH, OWHRLE, and CLE, and other time-variant variables, was assigned to each case and control herd using the approach described.

Statistical analysis
In the univariable analysis, conditional logistic regression based on Gaussian distribution was used to evaluate each factor of interest. Variables indicative of local exposure and other categorical variables were additionally tested using exact conditional logistic regression. Univariable models were first examined and variables that were either significant using a liberal P value of <.15, or were considered confounders in the target population, were included in the multivariable modeling. Although the main effect of time (matching variable) cannot be included as an explanatory variable in this analysis, the significance of interaction between time and each explanatory variable can be investigated [22]. Such interaction terms were therefore tested with each variable that, on its own, was significant at P < 0.15. A likelihood ratio test was used to test interactions with linear and quadratic effects of time. Testing was performed at this stage, instead of at the end of the multivariable analysis, because it was anticipated that the ratio of number of observations and strata per variable evaluated would be more desirable during this univariable stage. The choice of the final multivariable model was based partly on biological rationale and partly on statistical significance. A final model could not be fitted using exact conditional logistic regression. Therefore, a full model containing all variables was fitted using GD assumptions, and additionally, a "full model" with binary variable for easting (>600 km East) was fitted using exact conditional logistic regression. Since residual diagnostics were not available for conditional logistic regression, a model containing identical covariates as the final model in the conditional logistic regression was refitted using ordinary logistic regression and influence statistics and residuals were evaluated. Identification of any influential points prompted refitting the model using conditional logistic regression, omitting influential points and assessing the change in coefficients and their statistical significance. The exact P-value and the associated 95% confidence intervals were based on the score method. The best median unbiased estimates were used when point estimates could not be obtained for coefficients. Models were run using commercial software (Stata 10 SE, College Station, TX).