Modelling the dispersal of the two main hosts of the raccoon rabies variant in heterogeneous environments with landscape genetics

Predicting the geographic spread of wildlife epidemics requires knowledge about the movement patterns of disease hosts or vectors. The field of landscape genetics provides valuable approaches to study dispersal indirectly, which in turn may be used to understand patterns of disease spread. Here, we applied landscape genetic analyses and spatially explicit models to identify the potential path of raccoon rabies spread in a mesocarnivore community. We used relatedness estimates derived from microsatellite genotypes of raccoons and striped skunks to investigate their dispersal patterns in a heterogeneous landscape composed predominantly of agricultural, forested and residential areas. Samples were collected in an area covering 22 000 km2 in southern Québec, where the raccoon rabies variant (RRV) was first detected in 2006. Multiple regressions on distance matrices revealed that genetic distance among male raccoons was strictly a function of geographic distance, while dispersal in female raccoons was significantly reduced by the presence of agricultural fields. In skunks, our results suggested that dispersal is increased in edge habitats between fields and forest fragments in both males and females. Resistance modelling allowed us to identify likely dispersal corridors used by these two rabies hosts, which may prove especially helpful for surveillance and control (e.g. oral vaccination) activities.


Introduction
Understanding spatiotemporal patterns of pathogen spread is crucial to implement effective actions to contain epidemics (Ostfeld et al. 2005;Vander Wal et al. 2014). Wildlife pathogens, including some that can be very harmful to humans and livestock, are transmitted when infected hosts come in direct or indirect contact with uninfected individuals. In both directly and indirectly transmitted diseases, the extent and speed of propagation is expected to be linked to the dispersal ability of the hosts (Biek and Real 2010). Thus, information about movement and dispersal of hosts is required to identify potential spread pathways. As an example, rivers and highways appear to slow the spread of chronic wasting disease in white-tailed deer (Odocoileus virginianus), most likely because they act as barriers to dispersal and gene flow for this species (Blanchong et al. 2008). Similarly, large rivers hamper gene flow in raccoons (Procyon lotor) and may reduce the propagation of the raccoon rabies variant (RRV; Cullingham et al. 2009). Control operations that aim at containing and eventually eradicating a given disease are thus likely to be more efficient if positioned alongside these barriers to strengthen their effect (Russell et al. 2006). This strategy was adopted and prevented the northward spread of RRV in Ontario (Canada), in 1999 (Rosatte et al. 2001). Distribution of oral vaccine baits along major rivers to control rabies was performed as early as the 1980s for red foxes (Vulpes vulpes), eventually contributing to the elimination of rabies from Switzerland (Wandeler et al. 1988).
Not all environmental barriers to host dispersal and pathogen dissemination, whether natural or of anthropogenic origin, are spatially discrete or easily identifiable, such as rivers and roads, but may instead be continuous or follow a gradient of biotic or abiotic conditions (Storfer et al. 2007). Climate (Geffen et al. 2004), elevation (Shirk et al. 2010) and presence of unsuitable habitats (Goldberg and Waits 2010) are all examples of such limiting conditions. Based on the ecology, behaviour and dispersal ability of host species, these features may restrict pathogen dispersal or promote it through dispersal corridors. Integrating environmental features in models of disease spread can help in predicting the spread and geographic expansion of a disease (Ostfeld et al. 2005).
Different approaches are available to understand the effects of habitat composition on animal movement and dispersal. The first relies on trapping of animals to determine their resource selection and density (Manly et al. 2002), which requires important time and resources to gather large sample sizes. Studies have also been conducted using very high frequency (VHF) transmitters and, more recently, global positioning system (GPS) radio-telemetry to track animal movement and analyse habitat use (Cagnacci et al. 2010). Despite constant technological improvements, collecting large GPS data sets remains very costly and logistically challenging for several species. Spatial simulations can also be used to characterize factors affecting movement and connectivity among individuals in a population (Russell et al. 2006;Rees et al. 2013). While these models can bring insights on the links between habitat and dispersal, the quality of model outputs will depend on an appropriate characterization of ecological processes, which can only be obtained through empirical evaluation. Finally, another approach relies on tools provided by landscape genetics, a discipline integrating aspects of population genetics, landscape ecology and spatial analysis. This field has tremendously progressed in the past 10 years (Manel and Holderegger 2013). Typically, landscape geneticists are interested in describing how gene flow among populations or subpopulations is influenced in often heterogeneous, fragmented landscapes, leading to estimates of functional connectivity (Manel and Holderegger 2013). However, measuring gene flow among groups of individuals imposes limitations on the interpretability of results in terms of functional connectivity because (i) there may be important discrepancies between gene flow and ecological dispersal, that is, movement among habitat patches may not necessarily be associated with opportunities for mating (Garant et al. 2007), and (ii) gene flow measured among populations reflects migration that has occurred for several generations in the past and may not accurately reflect current ecological processes (Epps et al. 2007), including sex-specific differences. Ideally, the operational unit in landscape genetics should be the individual (Manel et al. 2003), in which case estimates of pairwise genetic relatedness can be used as the response variable to model landscape connectivity according to habitat features (Segelbacher et al. 2010;Etherington 2011;Shafer et al. 2012).
Rabies is enzootic to many species of bats and carnivores throughout the world and has a relatively long average incubation period (between 30 and 90 days) in comparison with a short morbidity period (2-10 days) that almost always leads to death (Leung et al. 2007). In eastern North America, the predominant terrestrial rabies strain is the RRV, which has spread in wild populations of both raccoons and striped skunks (Mephitis mephitis, hereafter skunks, Guerra et al. 2003). This rabies variant was historically restricted to Florida, but infected raccoons were inadvertently moved to Virginia in the late 1970s and the virus has since expanded northward at a rate of 30-50 km/ year (Rupprecht et al. 2002). In Canada, it was first detected in southern Ontario in 1999 (Rosatte et al. 2001), then in New Brunswick in 2000 and finally in Qu ebec in 2006 (Rees et al. 2011). Here, we used estimates of genetic relatedness derived from multilocus microsatellite genotypes to determine which landscape features promoted or limited dispersal of the two main hosts of RRV in an intensively studied area of southern Qu ebec where this viral disease is still under surveillance, control and research activities (Boyer et al. 2011;Houle et al. 2011;Rees et al. 2011;Côt e et al. 2012;Mainguy et al. 2012;Talbot et al. 2012).
Previous work in the study area (south-eastern Qu ebec) has shown very little genetic structuring in resident raccoons and skunks, with highways and rivers generally generating a rather weak effect on patterns of genetic differentiation Talbot et al. 2012). These equivocal results may conceal the effect of unmeasured spatial variables and do not allow modelling mesocarnivore dispersal at the landscape scale. Our main objective here was to build on the population genetic results obtained in the previous work, using an approach that applies landscape genetic analyses and spatially explicit models, to predict the most likely pathways of skunk and raccoon dispersal and, by extension, terrestrial rabies spread in this area. Based on the ecological knowledge of habitat use by these two hosts, we expected dispersal in both species to be reduced in agricultural fields, but did expect movement to be increased in habitats characterized by a high density of edges (e.g. Glueck et al. 1988;Dijak and Thompson 2000;Larivi ere and Messier 2000). We expected skunks to be more sensitive to the presence of fields and residential areas than raccoons, as raccoons typically show a greater affinity for dispersal and use cornfields and other human-related food sources (Riley et al. 1998;Prange et al. 2004). We also expected females to be more sensitive to landscape structure than males in both species, as dispersal is usually male-biased in mammals in general (Greenwood 1980), including raccoons and skunks (Cullingham et al. 2008;Côt e et al. 2012;Talbot et al. 2012). To our knowledge, this is the first empirical work that addresses movement of these two important rabies hosts in a spatially explicit landscape genetics framework and also the first attempt to quantify the effect of habitat composition on their dispersal. Such work is important to refine predictive models of rabies propagation that use rivers (or other discrete barriers such as mountain chains) and human density indices to predict the rate of propagation of rabies (e.g. Smith et al. 2002;Russell et al. 2005Russell et al. , 2006.

Study area
Our study area was located in southern Qu ebec, Canada (45°23 0 N, 72°43 0 W), where all known positive cases of rabies in the province of Qu ebec have been recorded between 2006 and 2009 ( Fig. 1). We used biological samples collected on raccoons and skunks over 3 years (2008)(2009)(2010) in this RRV epizootic region, over an area of approximately 22 000 km 2 (Fig. 1). This area corresponds to the so-called RRV monitoring area where rabies-related surveillance has been increased since 2006. From east to west, the study area follows a gradient of increasing agricultural intensification and urbanization (Ghilain and B elisle 2008) where hayfields and pastures in the east are gradually replaced by large-scale, continuous row cropping for corn, cereals and soyabean. Forest cover also follows this gradient, as it becomes smaller and more fragmented along the gradient of agricultural intensification.

Sampling
Tissue samples used in this study were collected by the Qu ebec Minist ere des Ressources Naturelles et de la Faune (MRNF) and its partners between 2008 and 2010 during surveillance and control activities (see Rees et al. 2011;Côt e et al. 2012;Talbot et al. 2012). They included samples either taken from recovered road-killed animals (Rees et al. 2011) or collected during post-oral rabies vaccination campaigns . All samples were georeferenced using a hand-held GPS device. A skin biopsy was collected from the ear with a 2-mm punch for genetic analyses. Samples were stored in 95% ethanol until DNA extraction. Because many animals were sampled over small areas due to field activities conducted in specific zones (% 20 individuals/100 km 2 zones), a random subset of individuals (up to 3 individuals/25 km 2 ) was selected for each year to generate a sample as evenly distributed as possible over the study area (Fig. 1). A total of 330 raccoon samples (192 males and 138 females) and 345 striped skunk samples (195 males and 150 females) were retained for this study. There was no geographic bias in the locations of samples for both sexes.

Genetic analyses
DNA extraction was conducted using a protocol described in Chambers and Garant (2010). Microsatellite polymorphism was assessed at ten loci developed for raccoons (see Côt e et al. 2012 for the list of loci used and PCR protocols) and at nine loci developed for skunks (see Talbot et al. 2012). DNA amplification was performed using GENEAMP SYSTEM 9700 thermocyclers (Applied Biosystems, Foster City, CA, USA). Genotyping was performed using an AB 3130 DNA sequencer (Applied Biosystems), and allele size was scored using GENEMAPPER 4.0 (Applied Biosystems). All microsatellite loci were tested for departures from Hardy-Weinberg equilibrium and linkage disequilibrium with a sequential Bonferroni correction, and indices of genetic diversity (number of alleles and observed and expected heterozygosity) were computed using GENEPOP 4.0 (Raymond and Rousset 1995) and POPGENKIT (Rioux Paquette 2011). Finally, molecular sexing was conducted using a modification of Shaw et al. (2003) as described in Côt e et al. (2012).

Landscape genetics
We applied two types of analyses that both aim at finding the model that best explains patterns of genetic distances among individuals. The first approach (multiple regression on distance matrices or MRM; Legendre et al. 1994;Lichstein 2007) provides a way to statistically test the influence of a set of variables (e.g. geographic distance and land cover types) on pairwise genetic distances in a regression model and retain only those that have a significant effect. The second approach (isolation-by-resistance models or IBR, McRae 2006) is based on circuit theory and treats the landscape as a surface with various electrical resistances, while dispersing individuals (or gene flow) are analogous to electrical current. Any number of models with different landscape resistances can be tested, but we restricted their number to a small number of competing hypotheses that integrated results from MRM. The model that best fits genetic distances can further be used to illustrate dispersal at the landscape scale. The following paragraphs describe these statistical analyses in details.

Calculation of genetic distance
An estimator of pairwise genetic relatedness (r xy ; Wang 2002) was calculated using SPAGEDI 1.3 (Hardy and Vekemans 2002). We converted relatedness coefficients to genetic distances using 1r xy for easier interpretation. In all analyses, these were used as an index of genetic dissimilarity among all pairs of sampled individuals.
Multiple regression on distance matrices MRM is derived from partial Mantel tests of matrix correlations, in which predictor matrices (in this case, landscape variables, see paragraph below and Appendix A) are used to explain variation in a response matrix (genetic distance). Because of nonindependence of observations in pairwise distance matrices, MRM comprises a permutation procedure that takes into account the structure of distance matrices (i.e. keeping observations belonging to a given individual together) to assess the statistical significance of explanatory variables (Legendre et al. 1994).
To generate landscape variables, southern Qu ebec was mapped using aerial photographs provided by the MRNF (Fig. 1) in ARCGIS 9.2 (Environmental Systems Research Institute, Redlands, CA, USA). Landscape was characterized by assigning each pixel to one of five categories of land type (see Table 1). We selected these 5 habitat categories based on a previous study on raccoon habitat conducted by Houle et al. (2011) in this region. A buffer zone was traced around each pair of individuals, within each species. The length of the buffer was equal to the Euclidian distance between the two individuals (ranging from 0 to 305 km;  mean AE SD: 84.0 AE 47.4 km). The width of the buffer was equal to 4 km. This arbitrary value was selected as a reasonable width since the home range of raccoons in rural Ontario is typically < 4 km 2 (Rosatte 2000). An illustration of the buffer tracing method is provided in Appendix A for clarity. Landscape composition (percentage of a given buffer covered by each habitat type) was calculated using HAW-TH'S TOOLS (Beyer 2004) in ARCGIS 9.2. The length of edges between forested and agricultural areas (fields) was calculated within each of the buffer zones, also using HAWTH'S TOOLS, and divided by the area of the buffer zone, to obtain the edge density separating each pair of individuals of the study. In addition to the continuous landscape variables, major rivers and highways (i.e. landscape barriers) were also included in the analyses. These barriers were also mapped with aerial photographs provided by the MRNF using ARCGIS 9.2. We selected rivers that are known to maintain an important discharge throughout the year (water flow > 0.6 m/s, as in Talbot et al. (2012): Richelieu, St. Franc ßois, Châteauguay, Yamaska and Magog rivers). We selected highways based on high speed limit (100 km/ h) and the absence of crossroads. We calculated how many times the Euclidian distance between each pair of samples intersected with a section of a major river or a highway using HAWTH'S TOOLS. MRM models were computed for the two species when including all individuals and then separately for both sexes using the R package ECODIST v. 1.2.2 (Goslee and Urban 2007), and statistical significance was assessed with 10 000 permutations in each case. Standard errors on model coefficients were estimated by jackknifing individual data. Initial MRM models included geographic distance (km), number of major river crossings, number of highway crossings, proportion of fields, proportion of residential areas, proportion of wetlands and edge density. A backward selection procedure was then applied (P = 0.05) to progressively eliminate nonsignificant variables from the models. While the suitability of stepwise methods (including backward elimination) to select regression variables has been debated (see Whittingham et al. 2006;but Murtaugh 2009 for a counterargument in support of their use), we mainly chose this approach because common alternatives (e.g. comparing models on the basis of Akaike information criterion [AIC] values) cannot be applied in a dissimilarity framework in which independence of observations is not respected. This is a common issue in landscape genetics because pairwise genetic distances may not be independent (Goldberg and Waits 2010). The proportion of forested areas was not included in the models, as it was very strongly correlated (negatively) with the proportion of agriculture areas (r < À0.9) and explained a smaller proportion of the variance. The proportion of open water was also excluded because the number of rivers was already included as a discrete landscape variable.

Isolation-by-resistance models
In the past, a large number of studies have relied solely on Euclidian distance (isolation-by-distance model, or IBD) or on the number of discrete barriers (isolation-by-barrier model, or IBB) between samples to explain genetic differentiation. However, the IBR framework (McRae 2006) has recently emerged as a valuable approach in landscape genetics to assess 'effective distance', that is, the actual distance an individual would need to travel between two points, assuming he chooses a path of least resistance (Amos et al. 2012). After assigning a resistance value to every element of the landscape, an algorithm computes resistance distances between pairs of points on the surface. In several cases, resistance distances explain patterns of genetic differentiation much more accurately than regular Euclidian distances (e.g. McRae and Beier 2007).
We used CIRCUITSCAPE (McRae and Shah 2009) to compute pairwise resistance distances. This program uses a raster file (map) and a list of coordinates (sampling points) as input. Seven competing models of landscape resistance, each characterized by specific resistance values assigned to landscape elements, were evaluated (see Table 2 for After running CIRCUITSCAPE and obtaining pairwise resistance distances for each model, they were compared by fitting linear models regressing genetic distance against resistance distance (i.e. IBR models), with R version 3.0.1 (R Development Core Team 2013). Because of nonindependence of pairwise distances, we could not rely on the use of an information theoretic method (e.g. AIC values) to compare models. Instead, we applied the pseudo-bootstrap approach of Worthington Wilmer et al. (2008) to select the best IBR model. This approach consists in retaining, for each pseudo-bootstrap replicate, only a randomly chosen subset of distance values that are completely independent from one another. For instance, for our data set of 330 raccoons, a maximum of 165 independent distance values can be obtained. For each replicate, linear models for the different resistance models are fitted and the one with the lowest AIC is selected as the best model. After a large number of replicates, the model most often selected is considered the best IBR model. In our case, 50 000 pseudobootstrap replicates were performed for each species and each sex within species with a custom-made R script (available on request from the authors). The best models were used to produce maps of landscape resistance to visualize zones of high dispersal (i.e. corridors) and zones of high resistance across the study area.

Results
For the two species, no pair of loci exhibited significant linkage disequilibrium. A single locus showed a significant departure from Hardy-Weinberg equilibrium (locus PLM20 in raccoons). All loci were highly variable: the number of alleles per locus in raccoons ranged from 5 to 26, H O between 0.685 and 0.888, and H E between 0.754 and 0.931. In skunks, the number of alleles per locus ranged from 7 to 15, while H O varied between 0.588 and 0.881 and H E between 0.620 and 0.896. Polymorphism information for all loci is provided in Appendix B. These summary results indicated that the microsatellite data sets were suitable for further analyses.

MRM analyses
Globally, MRM models showed that raccoon dispersal was generally less sensitive to landscape composition than in skunks. These models also showed that females of both spe-cies were more sensitive to landscape structure than males. All results from MRM models, including the final models retained after backward selection of variables, are reported in Table 3. In raccoons, the analysis including all individuals and the one restricted to males showed that the only significant predictor of genetic distance among individuals was geographic distance. In females, the only significant variable in the final model was the proportion of agricultural fields, as genetic distance among females increased with the proportion of fields in the landscape. In skunks, for models computed either with all individuals or with females only, genetic distance increased with geographic distance and proportion of agricultural fields, but decreased as the amount of forest edges increased in the landscape (Table 3). In contrast, male skunks did not exhibit detectable sensitivity to agriculture.

Isolation-by-resistance models
In raccoons, the pseudo-bootstrap procedure revealed the greatest support for a simple model of isolation by distance (model 6 from Table 2), both when considering all individuals or males only (Fig. 2). In females, model 4 was best supported, indicating that agricultural fields show a greater resistance to dispersal than forested patches or residential areas. However, support for models 6 (IBD) and 7 (edge effect) was only slightly lower in female raccoons, and support for model 7 was found to in fact surpass that of model 4 when edge resistance was fitted as being only 1/50th of forest resistance (see Appendix C). In skunks, all three set of analyses [(i) all individuals, (ii) males and (iii) females] indicated the greatest support for model 7 (Fig. 2), the most complex model of landscape resistance considered (edges with the lowest resistance, followed by forest fragments, residential areas and then fields; Table 2). For both species, none of the selected models included highways or rivers as greater barriers to dispersal than agricultural fields.
Using resistance values from the selected models for each species and sex combination, we generated current maps, defined as such because of the analogy with electrical resistance (McRae and Shah 2009), illustrating the hypothetical paths of a raccoon or skunk dispersing from the last known record of rabies in Qu ebec to the northernmost boundary of our study site. Figure 3A shows the case of a female raccoon, while Fig. 3B represents the same hypothesis for a skunk (either male or female). The male raccoon scenario is not shown because a simple model of IBD best explained male raccoon dispersal. On these maps, dispersal corridors are visible from south to north in the central part of the study area, which is associated with the presence of forest patches in an otherwise agriculturally dominated landscape. Table 3. Values of MRM regression coefficients (b) in models explaining species-specific genetic distance among individual raccoons or skunks. For both species, analyses were conducted with all individuals at first and then separately for males and females. For each analysis, the complete regression model is shown (with all variables), as well as the final model resulting from backward elimination of nonsignificant variables.

Discussion
In this study, we used a combination of MRM and resistance analyses to identify landscape features affecting dispersal in raccoons and striped skunks and to predict the most likely pathways of terrestrial rabies spread by these hosts. As expected, we found that raccoon dispersal was generally less sensitive to landscape composition than Figure 2 Distribution of pseudo-bootstrap replicates according to which model of isolation by resistance (IBR) had the lowest Akaike information criterion (AIC) values; among seven different models considered (see Table 2 for landscape resistance values of all models), IBR models were fitted to explain genetic distances among raccoons (A-C) and striped skunks (D-F), first by combining data from both sexes and then by analysing each sex separately. A total of 50 000 pseudo-bootstrap replicates were performed for each analysis.
(A) (B) in skunks and that females of both species were more sensitive to landscape structure than males. Our results also suggest that the most likely pathways for northward host dispersal in the region are in the centre or our study area in corridors associated with forested fragments where the landscape sharply shifts from predominantly extensive farming to mostly intensive agriculture. Previous studies have shown no difference in dispersal distance or general behaviour between rabies-infected and other individuals in raccoons (Rosatte et al. 2006) and skunks (Greenwood et al. 1997) and no effect of the presence of rabid individuals on the resulting genetic population structure (Talbot et al. 2014). Thus, we argue that modelling the dispersal of hosts in heterogeneous landscapes can, by extension, allow a better understanding and forecasting of rabies spread, especially when combined with ecological epidemiological models (e.g. Russell et al. 2005;Rees et al. 2013).

Dispersal in raccoons and striped skunks
Our initial prediction that dispersal in raccoons would be less sensitive to landscape heterogeneity than in skunks was supported globally, and other studies support frequent long-distance dispersal in raccoons. For instance, 10% of natal dispersal distances were >10 km in southern Ontario (Cullingham et al. 2008), whereas in an agriculturally fragmented landscape in Indiana (USA), long-distance dispersers accounted for 10% of the population and there was no pattern of IBD among habitat patches (Dharmarajan et al. 2009). In the similarly heterogeneous landscape of our study site, geographic distance was the main predictor of pairwise relatedness in raccoons. This indicates that intrinsic raccoon behaviour and propensity to disperse are more important than landscape composition in shaping patterns of interindividual genetic distance at the scale of the study. The only significant landscape effect was observed in female raccoons, for which an increase in the proportion of agricultural fields for a given area leads to an increase in pairwise genetic distances. This result is suggestive of reduced dispersal in that habitat type. In raccoons, females are the most philopatric sex (Cullingham et al. 2008;Côt e et al. 2012), and thus, we expected dispersal to be more affected by landscape resistance in this sex. Our results also help explaining the results of Dharmarajan et al. (2009), who found that in an agricultural landscape, 50% of males dispersed over distances larger than 5 km, whereas this proportion was only 19% in females. Further support for this comes from a radio-telemetry study conducted by Beasley and Rhodes (2010) that indicated that female raccoons, as opposed to males, tend to concentrate their activities in remnant forested patches within agriculturally dominated areas.
Interpreting landscape resistance in terms of habitat quality can be misleading without field-based ecological evidence, because dispersal is influenced by individual and population conditions (Ronce et al. 2001). For instance, a high-quality habitat may reduce dispersal if all resources needed by an organism are found within a small area, but may stimulate dispersal if a habitat patch becomes overcrowded. Indeed, dispersal in raccoons is often considered to be partly driven through a density-dependent process (Russell et al. 2006), and local patch dynamics may explain patterns of dispersal in fragmented landscapes (Dharmarajan et al. 2009). In the case of the effect of agricultural fields on female raccoon dispersal, field studies of raccoon habitat use have suggested that raccoons generally avoid field interiors (Fritzell 1978;Glueck et al. 1988;Beasley and Rhodes 2010). This would support the idea that agricultural lands reduce raccoon dispersal because they represent low-quality habitat. Conversely, edges between forest patches and fields may still represent particularly suitable habitat for raccoons (Dijak and Thompson 2000;Barding and Nelson 2008), where they tend to be found in higher abundance when compared to agriculturally dominated areas, as previously documented in our study area (Houle et al. 2011).
In skunks, the model of greater dispersal in habitat edges received the strongest support in all three analyses (all individuals combined, males and females). We expected skunks to be more sensitive to landscape heterogeneity based on previous studies of their habitat use that emphasized the preference of skunks for ecotones (Larivi ere and Messier 2000;Frey and Conover 2006). While it is interesting that this preference was captured in our IBR analyses derived from microsatellite-based relatedness, its impact on the resulting current maps is rather weak, as modelling of dispersal across the landscape is very similar for skunks and female raccoons (Fig. 3A,B). These maps emphasize that, globally, raccoons and skunks respond similarly to landscape heterogeneity, and the IBR model of edge effects (no. 7, i.e. the best model in skunks) was the second-best model in all three raccoon analyses ( Fig. 2A-C), especially in females (see Appendix C also). This reinforces the proposed idea that control operations against RRV spread should particularly target areas where agricultural fields and forest patches are greatly interspersed (Boyer et al. 2011;Houle et al. 2011).
Perhaps one of the most surprising results of our study is the lack of support for rivers and highways as barriers to dispersal within our study area. Previous population genetic work covering the same geographic area suggested that in both species, rivers significantly constrained gene flow, albeit weakly Talbot et al. 2012). These analyses did not include landscape composition data. We found no evidence that rivers and highways act as greater barriers to dispersal than agricultural fields when accounting for landscape heterogeneity. This discrepancy emphasizes the need for caution when interpreting tests of IBD or IBB, as collinearity between either geographic distance or the number of rivers and other missing or excluded landscape variables may lead to spurious conclusions about the real processes driving spatial genetic structure (Cushman and Landguth 2010). Additionally, rivers that were included as possible barriers were selected based on their mean discharge, but it is still possible that some sections are much easier to cross than others, which could weaken the signal for their effect as natural barriers.

Value of combining MRM and resistance analyses
The approach we used in this work allowed us to describe the dispersal of two terrestrial vectors of rabies living in a heterogeneous landscape from pairwise relatedness estimates (converted to distance) among sampled individuals. Analyses of processes driving spatial genetic structure are only reliable if these processes are sufficiently homogeneous at the scale of the study area (Rioux Paquette et al. 2010;Segelbacher et al. 2010). For this reason, it is necessary to ensure that the studied populations represent single, genetically homogeneous clusters (Born et al. 2008). Through the use of Bayesian clustering algorithms, this verification has already been carried out, for both the raccoon and the skunk populations of the RRV monitoring area in southern Qu ebec Talbot et al. 2012Talbot et al. , 2014. Although microsatellite data can lead to somewhat imprecise estimates of pairwise relatedness (Van de Casteele et al. 2001;Csill eri et al. 2006), they should still reflect 'true' genetic relatedness (e.g. Mainguy et al. 2009) and thus should allow detecting species-and sex-specific dispersal patterns without relying on a very large data set. We believe that this makes our approach a valuable tool for future investigations of possible pathways of dispersion in the context of disease management.
Concerns about the statistical power and biases related to the application of the partial Mantel test and its derived forms (e.g. MRM) have been raised (Legendre and Fortin 2010;Guillot and Rousset 2013). Nonetheless, in landscape genetics, the problem often lies in choosing among several statistically supported models, in which case MRM may be especially helpful (Wagner and Fortin 2013). Simulation studies indeed suggest that partial Mantel tests and causal modelling with MRM are reliable to distinguish hypotheses of IBD, IBB or IBR (e.g. Balkenhol et al. 2009; Rioux Paquette and Lapointe 2009; Cushman and Landguth 2010; but see Cushman et al. 2013). A possible limitation of the method is the potential sensitivity of the results to different sizes of geographic buffers around pairs of sampled individuals. However, we evaluated this possibility and found that habitat composition values were highly correlated for different buffer sizes (e.g. correlations of values obtained for 2-km vs 4-km buffers were 0.940 and 0.936 for forest and fields, respectively, while they were 0.853 and 0.820 for 4-km vs 8-km buffers and remained above 0.750 until buffers reached 20 km). Thus, results should not be very sensitive to variations within this range of buffers. Here, we mostly used MRM to provide an empirical evaluation of the effect of various landscape elements on dispersal to restrict the number of considered IBR models to a small number. Expert opinion models of landscape resistance are often subjective and can be misleading (Shirk et al. 2010;Spear et al. 2010), so relying on empirical evidence to establish a list of hypothesis-based models is advised (Spear et al. 2010). It is possible to perform optimization procedures to obtain the resistance values that maximize the fit between genetic and resistance distances (Shirk et al. 2010). However, the calculation of IBR models might be computationally intensive, especially when using individuals as the sampling unit rather than populations, which increases the data set size. Furthermore, optimizing the fit of resistance and individual genetic distances would probably lead to overfitting considering the uncertainty associated with pairwise relatedness estimates derived from a relatively small number of markers (Van de Casteele et al. 2001) and knowing that the relative importance of landscape features may vary in different areas (Short Bull et al. 2011). Our objective was to compare the support of a restricted set of possible hypotheses regarding landscape resistance, but different resistance values may affect results (Spear et al. 2010). Most notably, in the case of female raccoons, for which IBR results were equivocal, lowering the value of edge resistance to 2% of the forest resistance (instead of 20% as in Table 1) leads to that model receiving slightly greater bootstrap support than the one that was previously best supported (Appendix C). In all other cases where one model clearly had greater support, results were not sensitive to edge resistance variation (Appendix C). Final results from IBR analyses and those from MRM were globally consistent.
Finally, the pseudo-bootstrap method of Worthington Wilmer et al. (2008) that we used for selection of IBR models provided a simple solution to the issue of inflated sample sizes when working with distance matrices (Wagner and Fortin 2013). It also provides an information theoretic method (i.e. using AIC) to select among competing IBR models without relying on partial Mantel tests between resistance and genetic distances, which can be problematic because of very high autocorrelation among IBR models (Cushman et al. 2013).

Implications on rabies surveillance and control
Most researchers who have investigated raccoon dispersal patterns have found that major rivers represent key semipermeable barriers (e.g. Smith et al. 2002;Cullingham et al. 2009). These provide ideal opportunities for efficient oral vaccination campaigns through reinforcement of natural barriers (Wandeler et al. 1988;Rosatte et al. 2001). Our landscape genetic results show that in this case, models accounting for landscape composition perform better than those that only include discrete barriers and provide further insights about previous population genetics work that had revealed weak effects of rivers on the gene flow of raccoons and skunks in the RRV monitoring zone at the same geographic scale Talbot et al. 2012). Depending on the study area, additional landscape features may influence dispersal, such as 'landscape shape', for example, a terrestrial constriction between water bodies ). At the landscape scale, being able to accurately model host dispersal is required to develop appropriate surveillance and/or control strategies, including the distribution of oral vaccine baits (e.g. Boyer et al. 2011). Our approach has allowed us to identify potential dispersal corridors between the last known rabies-positive cases close to the US-Canada border and the north of the RRV monitoring area. These corridors could be targeted as areas of high importance should RRV reach southern Qu ebec again. It is important to note that the current maps illustrate resistance between the source (the location where rabies cases have been detected) and ground node, which we chose, in this case, to place at the north end of the study area, considering the northward trend in RRV expansion. If the ground node was modified (e.g. if we wished to model dispersal to the north-eastern boundary of the map on Fig. 3), the resulting current map would change, illustrating the paths of least resistance in that direction. Recent developments of the CIRCUITSCAPE model may allow visual representation of connectivity in all directions across the landscape , so this limitation could be avoided depending on the research questions. In addition, landscape resistance values for the different habitat types could be incorporated in models of rabies propagation that already account for discrete barriers and human density (e.g. Russell et al. 2005Russell et al. , 2006. It would be especially interesting to combine our findings with the recent simulation work of Rees et al. (2013), which is the first to integrate habitat quality and heterogeneity in spatial simulations of rabies vaccination efficacy. Models combining landscape genetics with ecological epidemiology could be used to determine the most likely path of disease spread, at the landscape scale, between the location of a documented positive case and one or many nearby cities to assess, among other things, health risks for human populations.
Following the first RRV-positive case recorded in the province in 2006, the control campaigns performed by the Qu ebec MRNF in the past years have been effective: since 2009, there have been no reported cases despite the mainte-nance of intensive surveillance operations (Rees et al. 2011). Nevertheless, considering that 218 rabid raccoons and 88 rabid skunks we recorded in the four US states sharing a border with Qu ebec in 2011 alone (Blanton et al. 2012), the province likely will continue to be vulnerable to the possibility of receiving dispersing rabid animals. There are no rivers or other discrete barriers to dispersal along the US-Qu ebec border; consequently, applying resistance models that integrate habitat composition to model dispersal pathways, like the ones we presented here, may enhance current operations aimed at preventing RRV from re-entering the province. We believe this is especially relevant in the light of analyses showing that the economic efficiency of the provincial rabies management programme in the future will mainly depend on minimizing programme costs (Shwiff et al. 2013).