INTRODUCTION

Theoretical ecology has been frequently criticized for its lack of predictive power and practical application (^{Peters 1991}). This perception is in close association with the belief that ecology has no laws (^{Murray 2000}; ^{O’Hara 2005}). There seems to be several reasons to argue that ecology has no laws. Some authors blame the lack of predictive power and generality of most ecological studies (^{Lockwood 2008}), while others point to the complex nature of ecological phenomena (^{Benton et al. 2006}).

Many ecologists are prone to think that natural systems are so complex and that the number of fac- tors involved is always so large that understanding and predicting in ecology is only possible using complex and detailed models (^{Benton et al. 2006}). However, the argumentation about the complicated nature of ecology is something that cannot be determined a priori. For example, one of the main problems in population dynamics is to deduce the factors that are essential to understand and predict population growth. This often makes a model simpler. The critical question then, is not whether variation among individuals or any other factor is important, but rather, how to determine which factor is the most important in any particular situation. What we require, therefore, are unbiased methods to determine the forces that dominate or drive the dynamics of specific populations (^{Berryman 1999}).

Different models exist to understand and predict population changes. All of them have been used with more or less success in different contexts and for different systems. Some are as simple as the exponential growth equation (^{Malthus 1798}) that assumes a population will grow indefinitely without any (environmental) limitation. Others are very complex such as the agent based models (^{De Angelis & Gross 1992}) or those considering age and sex structure or cohorts of individuals, each with its own set of parameters (^{Leslie 1945}; ^{Caswell 2001}; ^{Lindström & Kokko 2002}). There are also many empirical regression-based statistical models that have been used to explain population change as a function of different variables. Somewhere in the middle are the theoretical models that include interactions with other populations (such as Lotka-Volterra’s predator-prey system), intra-specific competition and/or the physical environment such as Ricker or Royama’s models (^{Ricker 1954}; ^{Royama 1992}; ^{Berryman 1999}). With the growing availability of population counts time series, many authors have used these later models (or some form of them) to understand fluctuations, trends and main drivers of population change (e.g. ^{Lima & Berryman 2006}).

In Argentina, some of the longest time series available are those of rodents species in agroecosystems, which date back to the early 80’s (^{Andreo et al. 2009b}). Particularly, the exact same place (a railway bank) has been sampled since 1990 until 2014, first monthly for several years and lastly once per year at the moment of known high population abundance.

Former studies have used both statistical and theoretical models to understand the environmental drivers of rodent species in this area (^{Castellarini et al. 2002}; ^{Andreo et al. 2009a}; ^b; ^{Polop et al. 2012}; ^{Gomez et al. 2016}). One of these studies analyzed the effects of endogenous feedback structure (i.e., density dependence) and climatic/environmental factors on the dynamics of two of the most common rodent species in the assemblage, namely Akodon azarae and Calomys venustus (^{Andreo et al. 2009a}). The authors fitted simple population dynamics models (^{Royama 1992}) to 18 years of data (1990-2007) and found that the inter-annual numerical fluctuations of both rodent species showed first-order feedback structures; i.e., populations are regulated by intra-specific competition (^{Johst et al. 2007}). Furthermore, for A. azarae, plant cover and human induced changes in land-use seemed to represent the main exogenous perturbations, while C. venustus seemed to be more affected by rainfall and temperature (^{Andreo et al. 2009a}).

There are many studies that present candidate models and hypotheses to explain the population dynamics of different species. Not so many, however, provide further demonstration of the forecasting or predictive power of such models. One important question, therefore, is to determine whether previous models obtained by us (^{Andreo et al. 2009a}) are able to reproduce the dynamics of rodent species when new data is considered. We think that using independent data to test model predictions is the best approach to validate working hypotheses and determine which factor is the most important driving force. Moreover, using our previous models to predict independent data illustrates how general population dynamic principles can be employed to explain the population changes observed. We think that the question of model predictive power can only be solved through prediction. By this, we mean forecasting or predicting the future. The objective of this study is therefore to put previously fitted models for the two small mammal species, Calomys venustus and Akodon azarae (^{Andreo et al. 2009a}) under the test of independent data. For this purpose, we used new abundance data obtained by live trapping conducted once a year for 7 years (2008-2014) in border habitats of central Argentina agro-ecosystems.

MATERIALS AND METHODS

Study area

The study was carried out in the rural area of Chucul, southwestern Córdoba, Argentina (33°01’34”S; 64°11’21”W).The area is a typical undulating pampean plain (600-900 m a.s.l.) with temperate climate. The average temperature is 23 °C in January and 6°C in July. Annual rainfall averages 800 mm and it is mostly concentrated in summer months. The landscape mainly consists of individual crop fields surrounded by wire fences with borders dominated by weed species. In general, these border habitats are less disturbed than agricultural fields, maintaining relatively high plant cover throughout the year, thus providing good habitat conditions for small rodent species (^{Ellis et al. 1997}; ^{Simone et al. 2010}).

The system

Akodon azarae (Fisher 1829) and Calomys venustus (Thomas 1894) are two of the most abundant small rodents inhabiting the agro-ecosystems of the Pampean region in Argentina. They are usually found in relatively stable habitats with high vegetation cover, including crop field edges, roadsides, railway banks (borders), and remnant areas of native vegetation (^{Mills et al. 1991}). Both species are characterized as omnivorous, but A. azarae (25-30 g) consumes higher proportions of arthropods whereas C. venustus (55 g) eats higher proportions of leaves and seeds (^{Castellarini & Polop 2003}). Populations of both species show a strong seasonal variation in abundance, with a minimum in spring and peaks in autumn (C. venustus) or late autumn-early winter (A. azarae) (^{Mills et al. 1991}; ^{Castellarini et al. 2002}; Andreo et al. 2009b). Reproduction is also seasonal; the breeding season may last from September-October to April-June (^{Mills et al. 1992}; ^{Escudero et al. 2012}).

Importantly, the two species studied are hosts of different viruses: A. azarae is a reservoir of the virus Pergamino, one orthohantavirus genotype (^{Levis et al. 1998}), and C. venustus is a reservoir of the arenavirus Latino like (^{Calderón et al. 2011}). Though they have not yet been related to human diseases, it is relevant to understand their dynamics.

Population dynamics models

We used the non-linear logistic population model of discrete time proposed by ^{Royama (1992)}, derived from the logistic equation of ^{Ricker (1954)}. This model represents a randomly distributed population competing for a common resource:

Rt b exp^a·Xt−1+C ( Eq. 1)

where, Rt is the realized logarithmic per-capita population growth rate, b is a positive constant representing the maximum finite reproductive rate, X_t 1 is the logarithmic population density in t 1, C is a constant representing competition and resource depletion, and a indicates the effect of interference on each individual as density increases (^{Royama 1992}). Population density in each time period is

then obtained as:

X_t = X_t−1 + b − exp(a·Xt−1 +C) ( Eq. 2)

We used Eq. 2 to predict population density in the period 2008-2014. Exogenous effects of climate or vegetation cover over the endogenous feedback structure were included as extra terms:

Xt = Xt-1 + b − exp^{(a·Xt−1 +C+d·climate)}+ e·climate (Eq. 3)

where d and e represent lateral and vertical effects, respectively. Lateral effects operate over the carrying capacity of the system, while vertical effects act upon b, the maximum reproductive rate (^{Berryman 1999}). For further details about environmental variables tested and model ranking for each species, see Table 1 up to the fifth column, ∆AICc, in the Results section.

Small mammal data

To test our formerly fitted population models (see previous section), we used data from small mammal live-trapping conducted from 2008 to 2014 in a 6 x 10 grid (0.30 ha) placed in a railway bank (the exact same location of our previous study ^{Andreo et al. 2009a}). Stations in the grid were separated by 10 m. One Sherman-type live trap baited with a mixture of peanut butter and cow fat was placed in each station. Sampling was conducted for 4 consecutive nights during May-June of each year, when the maximum yearly abundance is usually registered (i.e., autumn in the Southern hemisphere). Population abundance was estimated as the minimum number of animals known to be alive (MNKA).

Climatic and NDVI data

Data on temperature and rainfall for the period 2008-2014 were provided by the laboratory of Agro-Meteorology from the University of Río Cuarto (Argentina). Normalized difference vegetation index (NDVI) time series for the mice trapping area was drawn from the MODIS MOD13A2 Collection 6 NDVI product available at the Land Processes Distributed Active Archive Center (LP-DAAC) site (https://lpdaac.usgs.gov/products/mod13a2v006/).

The MODIS vegetation index series (MOD13 series) was designed to continue the 20+ years of NOAA-AVHRR time series (^{Huete et al. 2002}) that we used in our previous work. In order to match the spatial resolution of the NOAA- AVHRR derived NDVI data series which we by then, we re-sampled MODIS NDVI images to a spatial resolution of 8 km. Remote sensing data was downloaded by means of the pymodis library (http://www.pymodis.org) and further processed with GRASS GIS 7.6 (^{GRASS GIS 2019}). Climatic and NDVI data was temporally aggregated as in the previous work, i.e., spring rainfall: accumulated rainfall from October, November and December (t 1); summer rainfall: accumulated rainfall from January, February and March (t); spring average temperature: average mean temperature from October, November and December (t 1); and annual minimum NDVI as the minimum NDVI value from winter t 1 to autumn t, with t being the time steps in years. Since our aim was to test our former models with new data, we estimated only the same environmental variables that appeared significant in our previous study (^{Andreo et al. 2009a}).

Prediction with new data

We used the first five models (according to the ranking given by AICc values) fitted for each species in our previous study (^{Andreo et al. 2009a}) to run total trajectory (TT) and one-step-ahead (OSA) predictions and compare with the 7 years of new data (2008-2014). Additionally, we also generated predictions with the endogenous model for A. azarae, since our previous work and others have indicated that this species’ inter-annual abundance changes might be mostly determined by intra-specific competition. Because TT and OSA predictions operate differently and might give different results, we assessed them both. Total trajectory predictions use only the first abundance value, in this case 2007 abundance, along with environmental variables in t or t- accordingly, to generate predictions for the whole 7 years period. OSA predictions, on the other hand, use the 2007 observed abundance value as starting point, and then it uses its own predicted X_t-1 in each time step t. In order to account for variability in abundances, we included the 95% confidence intervals to our deterministic predictions. These were estimated by adding a normal noise of µ = 0 and σ = 0.1 in 10 000 iterations and determining the 0.025 and 0.975 quantiles. While different, this is in line with ^{Sæther et al. (2009)}’s Population Prediction Interval (PPI) to quantify uncertainties in population fluctuation projections. To assess models’ predictions and compare them, we estimated the root mean squared deviation (RMSD) that represents the mean deviation of predicted vs. observed values in the units of the variable under study (^{Kobayashi & Salam 2000}; ^{Gauch et al. 2003}). Lower values of RMSD imply better agreement between observed and predicted values.

RESULTS

The predictions obtained for C. venustus and A. azarae (Fig. 1 and 2, respectively), provided a different model ranking than the one obtained formerly. This new ranking is also evident from RMSD values (Table 1).

The model that showed the smallest RMSD in the case of C. venustus was model C that included summer rainfall as a vertical effect. However, it fails to represent the observed decrease in abundances of the last 3 years (2015-2017). The only model that showed a slight tendency to display a decrease in C. venustus’ abundance was the endogenous model without any covariate (model D, Table 1) in OSA predictions. All the rest, failed to represent the decrease in numbers. Indeed, most models predicted the opposite.

Population fluctuations of A. azarae from 2008 to2014 were only captured correctly by the predictions of the endogenous model without covariates (model F.1, Fig. 2, Table 1). OSA predictions from this model showed a lower RMSD than TT, but still failed to predict 2008 and 2010 within 95% confidence. All the other models tested, that were deemed better than the endogenous according to ∆AICc, predicted much lower abundance values than observed and were difficult to rank since RMSD for TT and OSA provided different results.

DISCUSSION

We have used 7 years of new data to test population dynamics models that were fitted with 18 years of mice records. This exercise has allowed us to challenge our previous results and inferences regarding the population dynamics of C. venustus and A. azarae in agro-ecosystems from central Argentina (^{Andreo et al. 2009a}). The different model rankings that we obtained based on forecasting (Table 1) reinforce the need to keep testing our models with independent datasets, as prediction is key to verify/refute our hypotheses about population dynamics. Moreover, these results highlight the need to continue studying the forces that drive population changes and how these may change (strengthen or weaken) in different periods (^{Deitloff et al. 2010}).

In our previous work, we suggested that C. venustus’ dynamics was influenced by spring-summer rainfall and spring mean temperature acting as lateral effects over the system carrying capacity (^{Andreo et al. 2009a}). This hypothesis was supported by the fact that simulations of models C and D (Table 1) were outperformed by those of model E, even if the latter was the last in the ranking according to ∆AICc (Table 1, ∆AICc was slightly over the 2 units commonly used as rule of thumb). Using new data, however, we found that model C yielded better predictions than model E, which is now ranked in third place according to the RMSD values (Table 1). Model C includes a vertical effect of summer rainfall. This implies that climate is operating on the population growth rate by affecting b, the maximum finite reproductive rate (^Berryman1999). Thus, new results imply a partial agreement regarding one of the variables involved, but a disagreement on how the variable operates (over the maximum finite reproductive rate instead of overthe carrying capacity). The hypothesis supportedby model C seems ecologically plausible, i.e., higher summer rainfall might provide good conditions for higher reproduction levels in the cohort born late in this period (^{Polop et al. 2005}) and thus increasing theoverall b. It is however difficult to suggest that the dynamics of C. venustus is indeed affected by summerrainfall since the predictions of model C did not capture the decreasing trend starting in 2012 (Fig. 1); actually none of the models tested did. In a former work with monthly data, (^{Castellarini et al. 2002}) did not find rainfall effects either, only temperatures below 4 degrees were associated with C. venustus abundance fluctuations after a 4 to 6-months delay.

Fig. 1. Observed (points) vs. predicted (continuous lines) population abundance for Calomys venustus in the rural area of Chucul (Córdoba, Argentina). Each plot shows total trajectory and one-step-ahead predictions with the 95% confidence interval (dashed line). Runs began with an initial population abundance (MNKA) of 33 individuals corresponding to the observed value in 2007. 

Fig. 2. Observed (points) vs. predicted (continuous lines) population abundance for Akodon azarae in the rural area of Chucul (Córdoba, Argentina). Each plot shows total trajectory and one-step-ahead predictions with the 95%confidence interval (dashed line). Runs began with an initial population abundance (MNKA) of 32 individuals corresponding to the observed value in 2007. 

Table 1 Population dynamics models for Akodon azarae and Calomys venustus used for prediction and Root Mean Squared Deviation (RMSD) for total trajectory (TT) and one step ahead (OSA) predictions. References: Sum, summer (January, February and March); Spr, spring (October, November and December); Tmean, mean temperature; NDVImin, annual minimum NDVI. ∆AICc differences in Akaike information criterion corrected for small sample bias. These values come from ^{Andreo et al. (2009a)} and are reproduced here as to inform about previous models’ ranking. Models without lateral or vertical effects correspond to endogenous models. In italics, we mark the best models from ^{Andreo et al. (2009a)}. Models with the lowest RMSD in this work are in bold. 

For the case of A. azarae, there is a clear best model in terms of prediction of observed abundances and RMSD values; i.e., the endogenous model that repre- sents the effect of density dependence only (model F.1 in Fig. 2 and Table 1). In our previous work, we concluded that the minimum annual NDVI might act as a vertical forcing over A. azarae’s dynamics, i.e., determining the maximum finite reproductive rate b, since model A.1 (Table 1) showed the lowest ∆AICc and was able to reflect the dramatic fall that the species showed by mid 90’s (^{Andreo et al. 2009a}). In that opportunity however, we also fitted models without those low abundance “odd” points since there was evidence of strong exogenous perturbations in those years, and by doing so, we found that the dynamics was mainly endogenous, i.e., population changes were regulated by a first order feedback structure (^{Royama 1992}; ^{Berryman 1999}). Results obtained now are consistent with our hypothesis back then. Apparently, the minimum annual NDVI in the former study was useful to explain dynamics in those years of very low abundance, but for the rest, the inter-annual changes are regulated by intra-specific competition. This finding is supported by several studies that highlight the importance of spacing behavior and intra-specific competition in this species. In fact, A. azarae has a polygynous mating system that operates through female defense (^{Bonatto et al. 2012}), i.e., a minority of males controls multiple females leaving other males without access to them. Moreover, females of A. azarae are known to be strongly territorial and aggressive (^{Hodara et al. 2000}; ^{Suárez & Kravetz 2001}) during the breeding season, when they seem to mainly compete for highly covered areas (^{Hodara et al. 2000}; ^{Busch et al. 2001}) with higher densities of in sects (^{Bilenca & Kravetz 1998}). On a seasonal (intra-annual) basis, however, statistical models suggested that both rainfall and vegetation cover determine the abundances of A. azarae with lags of 3 to 6 months (^{Andreo et al. 2009b}).

Testing models against independent data helps uncovering periods governed by different forces. Clearly, long time series are essential for this kind of exercises. A related change in the same rodent assemblage was observed in the late 80’s - early 90’s when A. azarae had a significant increase in abundance and A. dolores almost disappeared (^{Polop et al. 2012}). It was suggested that environmental modifications mainly represented by changes in agricultural practices and crops sown surface may have affected the competition outcome of these two Akodon species (^{Polop et al. 2012}), favouring A. azarae. Are we facing a similar phenomenon between A. azarae and C. venustus? Further research should be carried out to answer such a question and discard migration for example, since it has been shown that unfavourable habitats pose little or no resistance to dispersal in C. venustus (^{Chiappero et al. 2016}). Furthermore, other recent studies in the region have not detected such a decrease in C. venustus’ abundances (^{Serafini et al. 2019}). In any case, it seems likely that human induced exogenous changes might affect the factors that are key for individuals’ persistence and survival (^{Deitloff et al. 2010}).

Even though in ecology and population dynamics, a 7-years period seems a long time, in statistical terms, 7 points it’s a (very) small sample. Therefore, it might as well be that our observations and predictions have a certain degree of uncertainty that we cannot account for with the data available. For example, we did not measure migration nor sampled predators as to assess the effects of these factors over C. venustus or A. azarae’s dynamics. We attempted to mitigate this uncertainty with the inclusion of confidence intervals around our deterministic predictions in line with ^{Sæther et al. (2009)} suggestion. Moreover, to keep our analysis as consistent as possible with the previous study, performed when the full Landsat archives were not yet open, we used satellite data of rather coarse resolution (i.e., MODIS data resampled to 8 km) to relate to a relatively small area where the mice trapping was carried out. The opening of the Landsat archive would allow a more detailed view of the area over the whole study period (^{Wulder 2016}) and the possibility to fit and test our models again from scratch. This is something we plan to address in the future.

CONCLUSION

In this work we have put our previously fitted population models under the test of independent data. This has allowed us to give further support to the hypothesis of intra-specific competition as the only regulatory mechanism of A. azarae’s abundances and challenge our previous inferences regarding C. venustus’ dynamics. None of the models obtained was able to reproduce population changes of C. venustus during 2008-2014. Either previous models were not good enough or the dynamics of the species was subject to some unrecorded local effect. Indeed, other researchers in the same study area have not observed such a decrease. Despite model ranking for A. azarae has also changed, the results of the predictions are consistent with our previous study and with existing literature, too.

This forecasting exercise highlights the relevance of challenging our results with new data to increase or decrease support for our hypothesis and improve our understanding of population dynamics, especially in the light of human induced environmental changes. In this sense, one of the advantages of theory based models such as those tested here is that they allow for a more mechanistic understanding of population changes.