Spatiotemporal dynamics of mesocarnivore populations

Mammalian mesocarnivores play critical roles in ecosystems via trophic interactions. The fluctuation of mesocarnivore abundance may cause trophic cascading throughout the ecosystems. However, little was known about density dependence and spatiotemporal dynamics of mesocarnivore populations. Northern raccoon Procyon lotor is a common mammalian mesocarnivore in North America, and is the host of many human infectious diseases. Few studies have investigated density dependence and hierarchical spatiotemporal dynamics of raccoon populations. We used 23-year time series of raccoon relative abundance from 14 wildlife management areas in Mississippi, USA, to test for spatial synchrony of raccoon populations with nonparametric correlation functions. We developed non-Gaussian state space models to detect density dependence of raccoon populations, and also used dynamic factor analysis (DFA) to determine the structure of the spatiotemporal dynamics of raccoon populations. The 14 raccoon populations lacked common trends, and were not synchronized. Strength of density dependence varied among raccoon populations, but was not related to the amount of hardwood forests. Differences in the structure of density dependence probably prevented populations from being synchronize by climatic variability. The raccoon populations exhibited greater local or idiosyncratic variability than regional variability in Mississippi. Northern raccoons have plastic life history traits permitting their population dynamics to closely track local variations in resource availability.

Animal populations of the same species inhabiting heterogeneous landscapes may become distinct in their temporal patterns with increasing distance between populations (Kareiva et al. 1990). Differences in climate, habitat conditions, and interspecific interactions may differentiate demography and dynamic patterns among multiple populations over space (Kareiva et al. 1990, Ranta et al. 1998, Michel et al. 2016). On the other hand, dispersal, climatic changes, and predation that operate on large spatial scales may synchronize the dynamics of many populations over landscapes (Bjørnstad et al. 1999, Zuur et al. 2003, Liebhold et al. 2004). The regional dynamics of many geographically distinct populations may be scale-dependent, including one or several common trends or population growth trajectories and local or population-specific variability (Zuur et al. 2003, Fauchald et al. 2017. Therefore, studies of the regional dynamics of multiple populations may shed light on the spatially hierarchical patterns as well as regional and location management of animal populations. Population spatial synchrony is a phenomenon where abundances or growth rates of many populations are positively correlated within certain geographic distances or are more correlated if the populations are geographically closer (Liebhold et al. 2004). Long-distance dispersal may homogenize the dynamics of multiple populations. Likewise, mobile predators can synchronize the dynamics of multiple geographically separated populations of prey (Haynes et al. 2009, Koenig andLiebhold 2016). Furthermore, climatic variability may also synchronize the dynamics of many populations (Moran's effect) on large spatial scales (Moran 1953, Ranta et al. 1997. Nonetheless, the Moran effect requires that the synchronized populations share the similar structure of density dependence with the same lagged terms of density of similar strengths or values of coefficients (Moran 1953). We explicitly extended the statistical assumptions of the Moran effect to: 1) the shared or similar population growth trajectory; and 2) residual correlations among synchronized population time series and with climatic drivers.
Northern raccoons Procyon lotor (hereafter referred to as raccoons) are a common mammalian mesocarnivore in the United States (US) (Lotze and Anderson 1979). The rise and fall of raccoon populations and other mesocarnivores may cause trophic cascading throughout ecosystems (Gehrt and Clark 2003, Ritchie and Johnson 2009, Suraci et al. 2016. Additionally, raccoons are the host of infectious diseases, such as rabies and yellow distemper, causing public health concerns (Jones et al. 2003, Smith et al. 2002. Therefore, understanding the regional dynamics of raccoon populations are important for regional wildlife management. However, studies of the spatiotemporal dynamics of raccoon populations have been hindered by the lack of multiple long-term population time series (Gehrt 2002). Despite the critical roles, few studies have investigated the density dependence of mesocarnivores compared to large herbivores and small mammals (Troyer et al. 2014a). Little was known regarding the population ecology of raccoons (Troyer et al. 2014b). In this study, we used 14 23-year capture per unit effort (CPUE) time series of raccoons collected from 14 wildlife management areas (WMAs) to investigate the spatiotemporal dynamics of raccoon populations in Mississippi, USA. First, we tested if 14 raccoon CPUE time series were synchronized. Second, we determined the structure of density dependence of the 14 raccoon populations. Findings of this study enhance understanding the population ecology of raccoons, but also help plan regional management of raccoons.

Raccoon capture per unit effort
We used the raccoon CPUE time series of 14 WMAs across Mississippi from 1983 to 2005 (see Fig. 1 of Davis et al. 2017 for the geographic locations of the 14 WMAs). The 14 WMAs ranged from 11 140 ha to 86 910 ha in size (Table 1). Annual raccoon hunting seasons of Mississippi lasted from 1 July to 30 September. Each hunter was required to purchase a hunting permit from the Mississippi Dept of Wildlife, Fisheries and Parks (MDWFP) with a bag limit of one raccoon. It was mandatory that hunters completed and returned permit cards to self-service stations to report the number of raccoons harvested and the number of hunting days during a season. The MDWFP maintained the database of annual raccoon hunting statistics, including the total number of raccoons harvested by hunters and the total number of hunter days, for each WMA. We divided the annual total number of harvested raccoons by the annual total number of hunter-days to calculate annual CPUE for each WMA. We indexed the relative abundance of raccoon populations using CPUE.

Data on landscape forest composition
We computed proportions of forest by WMA using the National Land Cover Database (NLCD) 1991(NLCD) , 2001(NLCD) and 2006 classified by the Multi-Resolution Land Characteristics Consortium (< www.mrlc.gov/ >). The NLCD 1991, 2001and 2006 were developed based on the LandSat images taken in 1990, 2000and 2005, respectively (Fry et al. 2011, Homer et al. 2007). We reclassified NLCD deciduous forest and woody wetlands as hardwood forest (i.e. deciduous trees as the dominant form of vegetation).
Average dispersal distance of northern raccoon ranged from 10 km to 20 km (Gehrt andFritzell 1998, Rosatte et al. 2010). We calculated proportion of hardwood forest within a 10-km and 20-km circular buffer centered at the centroid of each WMA, respectively. We averaged proportion of hardwood forest over the NLCD 1992NLCD , 2001NLCD and 2006 for each buffer size (hereafter, 10-km and 20-km proportion of hardwood forest, respectively).

Statistical analysis of spatial synchrony
We used non-parametric spatial correlation to detect the regional synchrony of the 14 raccoon CPUE time series (Bjørnstad et al. 1999). Spatial correlation between CPUE time series was computed using the function Sncf() in the R package ncf (Bjørnstad 2016). The Sncf function estimates spline spatial correlograms (i.e. spatial correlation as a function of geographic distance between WMAs) to demonstrate the pattern of between-population correlation with increasing geographic distance (Bjørnstad 2016). We used the bootstrap option of 2000 iterations to compute the 95% confidence interval (CI) of regional synchrony. If the bootstrapped 95% CI of regional synchrony excluded zero, we concluded that raccoon CPUE time series were synchronized among the 14 WMAs.

Non-Gaussian state space models of raccoon population dynamics
We assumed the number (N t ) of raccoons harvested per season to have a Poisson distribution with the number of hunter days (i.e. effort) as offset (E t λ t is the Poisson parameter that quantifies mean CPUE. Furthermore, we used the framework of generalized linear mixed effect models, assuming log λ t t X = , where X t is the non-observable true state of population dynamics (Thorson and Minto 2015). We used the Gompertz model of order-1 autoregression (i.e. AR (1)) to represent the dynamics of raccoon population (Gompertz 1825, Royama 1992: where a is intercept; b¢ represents the effect of direct density dependence; and e t is an independent normal variate, e N t e ∼ 0 2 ,σ ( ) . The variance σ e 2 measures the strength of environmental stochasticity. When the sign of b¢ is negative, the coefficient b (= 1 + b¢) of term X t-1 is less than 1.0. When b¢ < 0 ([1+b¢]<1.0), increases in its magnitude (i.e. absolute value) indicate increasingly strong, negative feedback between population density and per capita population growth rate. We also considered a state space model of density independence (Eq. 2): X a X e where a is intercept; and the coefficient of term X t-1 is 1 with b¢ = 0. We fit Eq. 1 and 2, respectively, to the time series N t and E t of raccoons for each WMA, and used model selection to detect the density dependence for each raccoon population. Previous simulations and empirical studies have demonstrated that the Gompertz model is more powerful than the Ricker model and other population models of nonlinear density dependence (Herrando-Pérez et al. 2012).
We used the template model builder (TMB) to fit Eq. 1 and 2 to the time series N t and E t (Kristensen et al. 2016). The R function optim was used to maximize the likelihood function of state space models to estimate unknown parameters a, b (=1 + b¢), and σ e 2 for each raccoon population (Bolker et al. 2013, Kristensen et al. 2016. We checked if the optimization algorithm converged for each fitting. The maximized likelihood value was used to compute Akaike information criterion corrected for small sample size (AIC c ) and ΔAIC c for each of the two models (Burnham and Anderson 2002). The best approximating model has the lowest AIC c . A model with ΔAIC c < 2.0 is a competing model for a raccoon population.

Dynamic factor analysis
Dynamic factor analysis (DFA) represents N (N = 14) population time series with M (1 ≤ M ≤ N) common latent trends (Zuur et al. 2003). A common trend is modeled by random walk. Each population time series is a linear combination of the M common trends with a factor loading (or coefficient) for each time series on each common trend (Zuur et al. 2003). Dynamic factor analysis can be expressed in the form of: where X t is the latent states of the dimension M × 1. Therefore, DFA is a dimension-reduction analysis of multivariate time series similar to traditional multivariate factor analysis (Holmes et al. 2012, Zuur et al. 2003. The number of latent trend M indicates the hierarchical structure or the number of clusters of regional wildlife populations. The variances of process error indicates local population variability. We used the natural log transformation to normalize the CPUE time series for DFA. Then, the transformed CPUE time series were standardized or z-scored before DFA to facilitate model convergence, following Zuur et al. (2003) and Ohlberger et al. (2016). There are five physiographic regions in Mississippi; thus, we fit five DFA models having 1-5 common trends to each of the 14 CPUE time series, respectively (Davis et al. 2017, Pettry 1977, Strickland and Demarais 2008. Each of the five models was fit with two different types of measurement error covariance matrices R, equal measurement error and no covariance (i.e. the form 'diagonal and equal') and unequal measurement error and no covariance (i.e. 'diagonal and unequal'), respectively (Holmes et al. 2012). The R package MARSS was used in R ver. 3.2.2 to analyze the time series to estimate the unknown parameters, including factor loadings, variance of measurement error, and variance of process error (Holmes et al. 2012). Information-theoretic approaches were used to select the number of common trends with AIC c and ΔAIC c .
We regressed coefficient b against 10-km and 20-km proportions of hardwood forest, respectively, using linear models. Significance of linear regression slope was tested at the significance level of 0.05.

Results
Long-term means of raccoon CPUE averaged 0.69 and ranged from 0.31 to 0.81 over WMAs. Raccoon CPUE time series were not synchronized (regional synchrony = 0.08, 95% CI = -0.02-0.19; Fig. 1). Model selection indicated that all raccoon populations but the O'Keefe, Old River and Sandy Creek WMAs had density dependence, with the AIC c of density dependent models being less than that of density independent models by 2.0 or more ( Table 2). The state space model had good fit with the R 2 of the regression of observed CPUEs against predicted CPUEs being >0.95 except for the O'Keefe (0.58) and Upper Sardis (0.81) WMAs (Supplementary material Appendix 1 Fig. A1). The strength of density dependence was different among the raccoon populations (Fig. 2). The 95% CIs ( b SE ± 1 96 . ) of the coefficient b of several raccoon populations did not overlap (Fig. 2). Strength of density dependence was not related to 10-km proportion of hardwood forest (slope = 0.36, p = 0.73, df = 12), nor to 20-km proportion of hardwood forest (slope = 0.13, p = 0.93, df = 12). Additionally, three out of the 14 raccoon populations had non-negative intercept (Table 2). Therefore, raccoon populations had different structures of density dependence.
Dynamic factor analysis suggested one or two common trends for the 14 raccoon CPUE time series (Table 3, Fig. 3). The two-trend DFA had the lowest AIC c ; however, single-trend DFA had its ΔAIC c of 0.65 and was a competing model. The R 2 of the regression of the observed CPUEs combined over the 14 WMAs against the combined CPUE predictions of the single-trend DFA was 0.22, but was 0.36 for the two-trend DFA. The residuals of the single-trend DFA were not synchronized (regional synchrony = 0.01, 95%CI = -0.07 -0.06). Furthermore, five of the 14 raccoon CPUE time series were related to the single common trend with negative factor loadings, whereas nine raccoon CPUE time series were positively related to the single common trend (Fig. 4). Therefore, the raccoon populations appeared to have substantially different population growth trajectories.

Discussion
The 14 raccoon populations, indexed by capture per unit effort, of Mississippi exhibited more localized or site-specific population variability than regional trends. The 14 raccoon populations were not synchronized, differing in the deterministic trajectories of population growth. Although we found evidence for density dependence in 11 of the 14 raccoon populations, variability in the strength of density dependence suggested a localized population dynamics shaped by site-specific carrying capacities. Consequently, regional climatic variability did not synchronize raccoon populations in Mississippi. Our non-Gaussian state space models for density dependence allowed for missing observations of bag size and used the Poisson distribution for counts or numbers of harvested raccoons to avoid transformation for data normalization. Our non-Gaussian state space models represent a natural way to model capture per unit effort data.
Raccoons are habitat generalists and are distributed in nearly all types of terrestrial ecosystems (Fritzell 1978). Raccoons have plastic, flexible life history traits, which allow raccoons to adapt to local ecological conditions. Owing to high plasticity, raccoon populations are sensitive to locally (e.g. at the patch scale) spatial and temporal variations in resource availability and landscape structure (Beasley et al.   2011). In this study, dynamic factor analysis identified two competing models: single-and two-trend models. The twotrend model had slightly higher explanatory power than the single-trend DFA. However, more than 60% of the spatiotemporal dynamics of the 14 raccoon populations was not explained by the common trends of DFA, indicating more local population variability than regional variability. Raccoons use tree cavities for resting and maternal denning during the reproductive season (Owen et al. 2015). However, our study suggested that raccoon populations were not limited by the overall availability of hardwood forests in Mississippi. It is large, mature trees, as a critical habitat attribute, that provide raccoons with sufficiently large tree cavities for reproduction and denning (Chamberlain et al. 2002, Henner et al. 2004). Abundance of female raccoons was positively related to den tree density in the agricultural landscape of north central Indiana, USA (Beasley et al. 2012). Raccoons are the effective predators of avian nests and small mammals, and also scavenge for food from human residence in the rural areas. It is reasonable to expect raccoon populations increase in the landscape disturbed by human residential establishment and agriculture owing to enhanced food availability (Beasley et al. 2011, Troyer et al. 2014b. With increasing raccoon population size and limited number of mature trees for maternal denning, raccoon populations may become density dependently regulated and become stabilized at the carrying capacity (Beasley et al. 2012). Raccoon populations reached an equilibrium with the finite rate of increase being stabilized at 1.0 in a protected area in Florida, USA (Troyer et al. 2014b). Model simulations demonstrated that raccoon populations of density dependence were characterized with the stabilized dynamics (Broadfoot et al. 2001). Furthermore, population sizes from 1994 to 2007 suggest that raccoon populations may have undergone density dependence in the Niagara Falls, ON, Canada (Rosatte et al. 2010). Raccoon populations on the 14 WMAs of Mississippi may have different carrying capacities, exhibiting variable strength of density dependent regulation.
Several studies estimated the population density and demography of raccoons using capture recapture methods and rigorous statistics (Beasley et al. 2011, Troyer et al. 2014b. However, these studies collected raccoon population data either from 10 or more sites only for a few years or only from a couple sites for 12-13 years. Precise estimates of spatiotemporal dynamics often require longterm (e.g. 10 or more years) data from 10 or more studies sites located across a large spatial scale. Long-term relative abundance indices, such as spot light survey count, harvest indices (often without harvest effort), and capture per unit effort from hunting bags and hunter's efforts like those used in this study, provide a useful source of multiple long-term   time series of relative abundance (Rolley and Lehman 1992, Beasley et al. 2012, Hagen et al. 2014). However, we caution that relative abundance indices may underestimate population abundance (Beasley et al. 2012, Leclerc et al. 2016. Capture per unit effort is assumed to be positively related to true population abundance. Violation of this assumption leads to biased index of population abundance. Furthermore, non-random failure to report hunting efforts may result in a biased index of wildlife abundance. Hunting bags without correction for hunting effort may reflect variation in hunting effort. Nevertheless, we used CPUE time series to index the long-term population dynamics of raccoon populations. Additionally, the conservative bag limit of one raccoon per hunter per season did not appear to alter survival, movement, and spacing behavior of raccoons in central Mississippi (Chamberlain et al. 1999(Chamberlain et al. , 2007. Our findings also have important management implications. First, raccoon populations exhibited substantial population variability from site to site (this study; Rosatte et al. 2010, Beasley et al. 2011. Thus, population monitoring needs to be conducted at many survey locations in a region for better understanding the spatiotemporal dynamic patterns of raccoon populations (Rosatte et al. 2010). Second, population control needs to be prescribed based on local ecological conditions at individual sites because substantial local dynamic components and high plasticity of raccoon populations. Last, substantial spatial variability in raccoon population productivity may create a metapopulation dynamics under managemental control (Broadfoot et al. 2001). Highly productive populations may rescue locally eradicated raccoon populations within annual dispersal distance (Broadfoot et al. 2001, Rosatte et al. 2010). Therefore, a systematic metapopulation control for regional raccoon management is desired for the effective control of high-density raccoon populations (Broadfoot et al. 2001).
In summary, we demonstrated a violation of the assumption of the Moran effect in the unsynchronized raccoon populations of Mississippi. The raccoon populations had different structures of density dependence (Table 2, Fig. 2), preventing climatic variability from synchronizing raccoon populations. Dynamic factor analysis also demonstrated that the raccoon populations had different patterns of population growth trajectories. The residuals of the single-trend DFA were not synchronized among the 14 sites. Our study represented one of few studies that explicitly modeled the discordant trajectories of unsynchronized populations using a hierarchical modeling approach.