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 and Liebhold 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.

Methods

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.

Table 1.

Area, dominant vegetation type, and geographic coordinates for 14 wildlife management areas (WMAs), Mississippi, USA.

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), N_t ∼ Poisson(E_tλ_t), where λ_t is the Poisson parameter that quantifies mean CPUE. Furthermore, we used the framework of generalized linear mixed effect models, assuming logλ_t = X_t, 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_f ∼ N (0,σ). The variance 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):

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 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² 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 difierent among the raccoon populations (Fig. 2). The 95% CIs ( ± 1.96SE ) 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.

Figure 1.

Spatial correlogram of capture per unit effort time series of northern raccoons in 14 wildlife management areas, Mississippi, USA, from 1983 to 2005.

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² 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.

Table 2.

Model selection of non-Gaussian state space models for density independence and density dependence of raccoon populations in 14 wildlife management areas of Mississippi, USA. Letters ‘a’ and ‘b’ represent intercept and coefficient measuring density dependence. Symbol AIC_c is Akaike information criterion corrected for small sample size.

Figure 2.

Coefficient b measuring strength of density dependence in the Gomperz population models for 14 northern raccoon populations in Mississippi, USA, from 1983 to 2005.

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 two-trend 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.

Table 3.

Model selection of dynamic factor analysis with different combinations of measurement error covariance matrix structures and different numbers (1–5) of common trends.

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, 2003, 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.

Figure 3.

(a) Singe- and (b) two-latent common trends of dynamic factor analysis for capture per unit effort time series of northern raccoons in 14 wildlife management areas, Mississippi, USA, from 1983 to 2005.

Figure 4.

Factor loading of the single-trend dynamic factor analysis for the capture per unit effort time series of northern raccoons in 14 wildlife management areas, Mississippi, USA, from 1983 to 2005. Factor loadings are the coefficients relating the common trend to observed capture per unit effort at each wildlife management area.

Several studies estimated the population density and demography of raccoons using capture recapture methods and rigorous statistics (Beasley et al. 2011, 2012, 2013, 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 long-term (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 hunters 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, 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, Rosarte 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.