Weather influences M. arvalis reproduction but not population dynamics in a 17-year time series

Rodent outbreaks have plagued European agriculture for centuries, but continue to elude comprehensive explanation. Modelling and empirical work in some cyclic rodent systems suggests that changes in reproductive parameters are partly responsible for observed population dynamics. Using a 17-year time series of Microtus arvalis population abundance and demographic data, we explored the relationship between meteorological conditions (temperature and rainfall), female reproductive activity, and population growth rates in a non-cyclic population of this grassland vole species. We found strong but complex relationships between female reproduction and climate variables, with spring female reproduction depressed after cold winters. Population growth rates were, however, uncorrelated with either weather conditions (current and up to three months prior) or with female reproduction (number of foetuses per female and/or proportion of females reproductively active in the population). These results, coupled with age-structure data, suggest that mortality, via predation, disease, or a combination of the two, are responsible for the large multi-annual but non-cyclic population dynamics observed in this population of the common vole.


Material and Methods
Study area. Data were collected from August 1979 to October 1996 at Septfontaines -Le Souillot (6.18°E, 46.97°N), in an area of 1400 ha (800 ha of farmland, 600 ha of forest), at an average altitude of 750-800 m above sea level 38 . There, 100% of the farmland was permanent grassland used for pasture and (high grass) meadow for cattle feeding in winter (minimum 6 months, November-March), with a productivity ranging from 2.3 tonnes of dry matter.ha −1 .an −1 for areas of lower productivity and up to 9.0 tonnes of dry matter.ha −1 .an −1 for rich meadows 37 . Legumes (alfalfa, clover) and silage were (and still are) not allowed with respect to the European Protected Geographical Indication specifications of the locally produced "Comté" cheese. A KML file with the bounding box of the study area is provided in the Supplementary Material. Trapping and dissection. Rodents were captured using INRA traplines. INRA live traps (15 × 5 × 5 cm) are suitable for species of body mass less than 50 g. Each line consisted of 34 live traps spaced 3 meters apart. Distances between traps corresponds to a standard design [39][40][41] for which a minimum of 2 traps should (theoretically) be available in each M. arvalis home range. Trap lines were set up for three nights and checked every day. Traplines were all set in grasslands but their position was changed from one session to the other in order to avoid over-trapping. A minimum distance of 100 m was kept between two traplines. Animals were euthanized by cervical dislocation, weighed and dissected for sex, reproductive status and age determination. Relative age was estimated by weighing the crystalline eye lenses after drying 42,43 . Pregnant or lactating females, as well as those with visible corpus luteum, were considered reproductive.
Trapping and animal handling was carried out in full accordance with the relevant European guidelines (Directive 86/609/EEC) and national regulations. The rodent species investigated in this study does not have protected status (see IUCN and CITES lists), and is listed as a pest, subject to control, under Article L201-1 of the Code Rural et de la Pêche Maritime of French law. INRA (Institut National de la Recherche Agronomique), the umbrella organization under which the field work was carried out, created its first ethical committee in 1998, 2 years after the end of our study. It was therefore impossible to get formal ethical approval prior to the study, and INRA does not issue retroactive ethical approvals. A similar research protocol used from 2014 to 2017 44  Meteorological data. Meteorogical data were kindly provided by Meteo-France. Monthly temperatures from 1983 to 1996 were from the Levier station (altitude 713 m), 6 km from the study area. From 1979 to 1982, temperature data from this station were not available, thus we estimated monthly values in Levier based on temperatures recorded at the Pontarlier station (altitude 831 m, 20 km from Levier and 15 km from the study area), after calibrating the two stations over the period 1983-1986 with linear regression (adjusted R 2 = 0.99, p < 0.00001). Assuming that low temperature was a limiting factor with possibly delayed effects, and following Krebs's opinion 1 that "the time scale of rodent dynamics is monthly or even weekly" (due to population turnover), we chose to investigate the effect of the deviation from the monthly average minimum temperature of the current month and of the previous three months before a focal event on reproduction parameters and population dynamics.
Monthly rainfall data were from the Levier station from 1979 to 1996. The variable under study was the deviation from the average rainfall of the month during the study period. As for temperature, hypothesizing delayed effects, we considered rainfall of the current month and of the previous three months of a focal event.
Statistical analyses. We investigated density dependence after log-transforming population estimates 1,45,46 obtained from the mean of the whole sampling region. In brief, we used relative density estimates from autumn and from spring separately, calculated autocorrelation coefficients and fitted autoregressive models of order 1, 2 and 3 to the data.
N t = abundance estimate for year t, and t = time in year ε t = the regression residuals.
Other data were either counts (e.g. a number of individuals) or presence/absence (e.g. reproductive/not reproductive) and were modelled against independent variables (e.g. relative age, temperature at given months, rainfall, etc.) using Generalized Linear Models (GLM) with Poisson log or Binomial logit link functions, respectively 47 since 1982 (all estimates before 1982 but one -May 1981 -are just indicative, since they are based on very small sample sizes).
Delayed effects of temperature and rainfall were investigated by modelling the response variable (either the reproductive status or the population growth rate) against meteorological variables of the current and 3 previous months using GLM (Binomial or Gaussian link function, according to case), as where x t = the response variable (reproductive status or population growth rate) y t = the deviation from the average at time t (minimum temperature of the month or rainfall) t = time in month ε t = the regression residuals All statistical analyses were performed in R (version 3.4.1) 48 .

Results
Over the 17 years of our study, we set up 1188 traplines and captured 8504 common voles (Table 1). Sample sizes ranged from as low as 1 animal per trapping session (July 1980 and August 1981) to as high as 72 animals (April 1988).

Population dynamics patterns.
Population dynamics showed a classical pattern with characteristic seasonal population declines almost every winter due to concomitant reproduction stoppages in combination with winter mortality, with the exception of winter 1982-1983 and perhaps winter 1981-1982 in which populations didn't decline over winter (Fig. 1, upper panel). Inter-annual variations of large amplitude were also observed. Figure 1, lower panel, shows log transformed population growth rate; we observed 7 spring and 2 summer population declines during the reproductive period. In detail: • 14 winters show clear decreases in population abundance, 2 show increases or stability; the first of these 2 latter paradoxical exceptions corresponds to an imprecise estimate of population change in 1981-1982, from close to zero in the autumn to 0.3 individuals.trapline −1 in the next spring; this is geometrically exaggerated on the chart by the log scale. The second in 1982-1983 corresponds to more reliable estimates and provides clear evidence of a lack of population decline in winter (stability or even a slight increase). • At least 6 winter declines are followed by a spring decline (7 if 1981 is included).
• At least two summer declines could be observed: • One (1987) in late summer after a winter and spring decline and early summer increase.
• One (1989) following a winter and spring decline.
• One (1994) spring-to-autumn decline (no summer estimate available) could also be observed.
Except for winter declines, attributable to seasonal reproductive patterns and subsequent uncompensated mortality, we did not observed temporal regularity. ( www.nature.com/scientificreports www.nature.com/scientificreports/ Density dependence. We did not detect delayed density dependence, and found evidence of direct density dependence in autumn only (Table 2): the larger the population growth at year t, the more likely a decline between year t and t + 1.
We also detected positive correlations between autumn and early spring densities (adjusted R 2 = 0.70, p < 0.0001) and between spring and autumn densities (adjusted R 2 = 0.34, p = 0.01); this might however be a trivial result, quantifying something that can be observed visually in Fig. 1: the time series considered is a sequence of periods of several successive years where densities are on average high, and of periods where densities are on average low, whatever the season.
Age structure and reproduction parameters. The age structure of the population varied in accordance with seasonal reproduction; lack of breeding over winter led to population aging until early spring when reproduction started again, and breeding throughout the summer until at least October. This led to rejuvenated populations from spring to summer, and sometimes into autumn (Fig. 2).
The proportion of reproductive females increased with age ( Fig. 3 and Table 3), from a baseline of 37% in the younger female age class (Fig. 3). Notably, the minimum crystalline lens weight for a pregnant female in our sample was 1.1 mg, which corresponds to an animal of about 13 days old 42,49 . Its body mass was 15 g (3/4 of the average body mass of an animal fully grown) and it carried 5 foetuses. Moreover, we found that the proportion of reproductive females increased with age but that they produced fewer foetuses (Table 3); however, the effect size could be considered nominal since models accounted for an extremely small proportion of the total deviance (R 2 = 0.01 for the two models) with an extremely low proportion of old females in the total population (Fig. 3).
Reproduction parameters and declines. The proportion of reproductive females was highly variable between years at the beginning of the reproductive season (April). Declines did not correlate with variations in the ratio of reproductive females; declines could (and most often did) occur with proportions of reproductive females similar to and sometimes higher than those observed during population growth (Fig. 4A). We observed similar patterns for the other months (July and August, Fig. 4A).  www.nature.com/scientificreports www.nature.com/scientificreports/ In October, close to the end of the reproductive season, large variations in the proportion of reproductive females could be observed between years (Fig. 4A). Moreover, the average number of foetuses per female and the proportion of reproductive females in the population were positively correlated (Fig. 4B).
Furthermore, we did not detect correlations between the population growth rate and the ratio of reproductive females and/or the number of foetuses, alone or combined, including interactions, at the beginning of the time span on which the growth rate was computed (e.g. between the reproduction parameters of April and the growth rate computed from April to July, etc.). We also took the sum of the proportion of reproductive females in April and October as a proxy of the reproduction season duration, with the assumption that the larger the ratio in April, the earlier the reproduction began that year, and the larger the ratio in October, the longer reproduction continued into the autumn. Here again we did not detect any correlation between the population growth rate from April to October and this proxy.
Meteorological conditions and reproduction. The ratio of reproductive females in April was significantly lower when previous months were colder than the average (Table 4). For the other periods and meteorological variables, correlations between the ratio of reproductive females, temperature and rainfall, although statistically highly significant, appeared complex, often time-lagged, and hard to interpret unequivocally (Tables 4  and 5). Clearly, meteorological conditions have an impact on reproduction, but differ among circumstances (e.g. rainfall as well as temperature may impact reproduction positively or negatively from one month to the other).

Meteorological conditions and population dynamics.
We did not detect any effect, delayed or direct, of temperature or rainfall, alone or combined (data not shown), on population growth (Table 6). For example, Fig. 1 shows that winters 1988-1989, 1989-1990, and 1993-1994 were among the warmest, but voles still underwent the "classical" winter declines, which in 1989 and 1994 were even followed by spring declines. Winters from 1984 to 1986 were among the coldest but this was the period when vole densities were higher, though there were some spring declines that were largely compensated by further population growth. 1989-1991 and 1994-early 1996 were periods when vole densities were lower, but these periods were far from the coldest; and while spring and summer 1984 were among the coldest, we observed one of the most dramatic increases in abundance in the entire time series. In short, population growth rates, whether they be positive or negative, appear wholly uncorrelated with meteorological conditions, or at least temperature and rainfall.

Discussion
Using a 17-year time series of Microtus arvalis abundance and reproductive activity and concurrent weather data, we have explored the relationship between meteorological conditions and reproduction in this species, as well as reproduction and population growth rates. We have detected significant but complex relationships between meteorological conditions (temperature and rainfall) and female reproduction, but we did not detect any clear relationship between female reproduction and population growth rates.  Table 3. Parameters of GLM: reproductive status (1/0) (Binomial logit link function) and number of foetuses of reproductive females (Poisson log link function) against relative age (expressed as crystalline lens weight, cryst).  31 have demonstrated consistent dampening of cycle amplitude associated with a reduction in winter population growth on a continental scale, suggesting a common climatic driver; Imholt et al. 34 , using linear methods, reported that weather parameters in winter and early spring are related to regional outbreak risk in eastern Germany; and Esther et al. 35 , using regression trees, found that between 12 and 20 weather parameters, depending on population and location, could influence relative vole density. In this last study, the number of days with snow cover in December and March, rainfall in spring, and the maximum temperature in October were all identified as key predictors of spring population densities the following year, and monthly maximum temperatures between February and June and the amount of precipitation in April and July were found to be correlated with population densities in fall. In short, meteorological conditions are generally considered to be important in M. arvalis populations dynamics.
Somewhat contrary to the studies listed above, we did not detect any impact of weather variation, direct or delayed, on population dynamics. However, we found strong direct and delayed correlations between weather parameters and reproduction: low temperature in winter delayed reproduction, decreasing the proportion of reproductive females the following spring. While we found other parameters and their combinations to have an impact on reproduction, these results appear complex and are quite hard to interpret in terms of the underlying   Table 5. Parameters of GLM (Binomial logit link function): reproductive status (1/0) against the deviation from the average rainfall of the month and each of the previous months. "R04" reads deviation from the average rainfall of April of the study period, etc. Last terms have been dropped when statistically not significant. We did to detect regression coefficients significantly different from zero in August (data not shown).
Scientific RepoRtS | (2019) 9:13942 | https://doi.org/10.1038/s41598-019-50438-z www.nature.com/scientificreports www.nature.com/scientificreports/ ecological processes. We also acknowledge that the list of weather variables we considered was not exhaustive, and that linear models might not always capture the observed complexity. Pucek et al. 50 suggested that ambient temperature may play a role in Polish bank vole (Myodes glareolus) population dynamics via it's effect on plant growth and subsequent food availability, and Tkadlec et al. 36 found that between North Atlantic Oscillation (NAO) index data and crop yield indices, the latter were better predictors of population change in voles in the Czech Republic. Blank et al. 51 found that topography and soil properties affect the local risk of common vole outbreaks in Eastern Germany. In the present study, we were unable to include plant productivity as a predictor because reliable data (NDVI, etc.) were not available at the relevant temporal and spatial scales and resolution over the study period; however, this may be a mechanism by which temperature and rainfall can influence reproduction in this population of M. arvalis, and deserves further study.
We found that population dynamics were completely disconnected from variation in reproduction (and thus also from variation in weather conditions). We did not find correlations, direct or delayed, between female reproduction parameters and population growth rate. Strikingly, apart from regular annual winter population decreases, we even observed population declines in spring and summer in periods when reproduction, as measured by the proportion of reproductive females and the number of foetuses, was among the most intense. Here it does not appear that reproductive change is an important component of population density change in our system. This is in contrast to some of the recent work on M. arvalis populations elsewhere; Pinot et al. 30 reported on a M. arvalis population in which steeper over-winter declines were attributed to more pronounced reductions in winter reproduction and recruitment following higher October densities; density appeared to have a slightly positive effect on survival, and it was concluded that mortality did not drive the steeper declines observed at high density. This is consistent with one of the competing paradigms regarding cyclic rodent populations: population cycling is a result of changes in reproductive parameters and not mortality. Modeling work has suggested that phase-specific changes in female age or size at sexual maturity, aided by changes in juvenile survival, are sufficient to generate multi-annual population cycles in rodents 26,27 , and empirical work in a variety of systems and rodent species has suggested that the length of the breeding season also plays a role ( 52,53 , reviewed in 1 ). If changes in reproduction are not responsible for the changes in population growth rates, then mortality must be considered. The analysis of the relative age structure of our population of M. arvalis each year conforms to seasonal reproductive patterns and does not indicate, for example, a lack of recruitment of specific age categories (e.g. young individuals during the reproductive season). This would have been the case if differential mortality occurred according to age categories (e.g. infanticide, differential predation or diseases, etc.). This rules out, for example, specific nest mortality as a unique parameter of the decline. Weather being excluded, social interactions, predation and diseases, alone or combined remain as alternative sources of mortality. Since the 1990s, most population ecologists have rejected intrinsic or self-regulation hypotheses as a possible explanation for population fluctuations 1 , thus predation and/or diseases become the most likely factors driving population dynamics here. High population densities may create conditions favourable to the spread of parasites and diseases, both through increased transmission but also increased host susceptibility [54][55][56] , and also to predator concentration in the area.
We did not detect delayed density dependence whatever the season considered; however, large but non-cyclic variations in density, with high density peaks >10 individuals.trapline −1 of 2-6 years, were observed over the 17 year of our study. Based on Spitz's 41 calibrations, a peak of 10 individuals.trapline −1 should correspond to approximatively 100 individuals.ha −1 on a landscape scale (5-25 km 2 ). Peaks of up to 20-30 individuals.trapline −1 on average were observed in our study, which corresponds to 200-300 individuals.ha −1 with large variability between plots (1-5 ha) and peaks of a magnitude of 1000 to 2000 individuals.ha −1 locally. This is corroborated by 2-3 km of transects walked across grassland which noted the presence of indices (holes, corridors, feces) in every 10 m interval (see e.g. 18,19 ). This pattern, of large but non-cyclic variations in density, is consistent with the range of dynamics observed in M. arvalis. Mackinrogalska et al. 57 , in an extensive review (36 data sets), found 6 populations with cycles shorter than 3 years, 26 ranging between 3 and 4.9 years, 6 larger than 5 years and 4 with no evidence of cyclicity in  Table 6. Parameters of GLM (Gaussian link function): population growth (ln(N t /N t−1 )) against the deviation from the average minimum temperature of the month and each of the previous months, and against the deviation from the average rainfall of the month and each of the previous months. "T0" and "R0" read deviation from the average minimum temperature of April and from the average rainfall at t 0 , "T1" and "R1" at t −1 (in month), etc.  24  Though Mackinrogalska's review 57 used the s-index to determine cyclicity 59,60 , a measure of cyclicity which has been criticized elsewhere 1 , it and the other studies listed above clearly show the considerable spatial and temporal variation in demographic behaviour within the same species 23 . They also show that simple unified explanations for all of those patterns have yet to be found. As discussed in Lambin et al. 24 , comparisons between the amplitude of population density variations in M. arvalis are quite difficult because most studies are grounded on various index-based trapping methods or qualitative historical records, and where density estimates exist they are scale-dependent. Furthermore, landscape descriptions generally do not provide quantified information on habitat composition and structure, and sampling design descriptions are rarely complete enough to ensure full comparability between studies. To our knowledge, single-factor hypotheses have considerably increased our understanding of possible causes for vole population regulation, but did not provide comprehensive explanations regarding the diversity of vole population dynamics (see 1 for a review); this is demonstrated by the fact that none of those single-factor frameworks have provided practical and efficient ways to control vole pests. On the contrary, vole control methods based on multifactorial hypotheses, that take specificities of the local environment into account, have proven to be efficient when implemented at relevant temporal and spatial scales e.g. for A. terrestris in Franche-Comté, France [61][62][63] . This approach consists of synergistically manipulating the several factors that are known to directly or indirectly effect the spatial and temporal dynamics of the vole populations (e.g. mole control, habitat disturbance such as cattle trampling, short grass, plough rotation, targeted chemical control and trapping on early vole colonies, installation of buzzard roosts, vole-predator protection, etc.).
Managed grassland ecosystems all over Europe are suitable for M. arvalis and provide an extremely large range of landscape configurations and habitat productivities. However, comparisons are dependent on a clear description of context, hence on landscape-and habitat-description quality. This can be quite complex but essential in heterogeneous landscapes (Table 7), and scale dependencies must be taken into account 64 . This should lead to considering landscape characteristics on a large scale relevant to vole dispersal and predator populations (e.g. ≃100 km −2 ). A relatively large number of studies deal with vole cyclicity patterns, but few of them have monitored populations parameters (reproduction, mortality, etc.) associated with density variation in the long term. We believe that grassland vole species with a similar ecology can benefit from such studies and also provide relevant comparisons where their range offers a variety of landscapes. Based on our field experience, this is for instance the case for A. terrestris, and maybe M. agrestis in Europe (in places where M. arvalis is absent for the latter), M. obscurus, M. gregalis, Ellobius tancrei, Lasiopodomys brandtii in central Asia, M. limnophilus, Phaiomys leucurus/ Lasiopodomys fuscus on the Tibetan plateau and its margins, all species extremely common with evidence of population outbreaks on large areas among a diversity of land uses and landscapes [65][66][67] .