It is time to mate: population-level plasticity of wild boar reproductive timing and synchrony in a changing environment.

On a population level, individual plasticity in reproductive phenology can provoke either anticipations or delays in the average reproductive timing in response to environmental changes. However, a rigid reliance on photoperiodism can constraint such plastic responses in populations inhabiting temperate latitudes. The regulation of breeding season length may represent a further tool for populations facing changing environments. Nonetheless, this skill was reported only for equatorial, nonphotoperiodic populations. Our goal was to evaluate whether species living in temperate regions and relying on photoperiodism to trigger their reproduction may also be able to regulate breeding season length. During 10 years, we collected 2,500 female reproductive traits of a mammal model species (wild boar Sus scrofa) and applied a novel analytical approach to reproductive patterns in order to observe population-level variations of reproductive timing and synchrony under different weather and resources availability conditions. Under favorable conditions, breeding seasons were anticipated and population synchrony increased (i.e., shorter breeding seasons). Conversely, poor conditions induced delayed and less synchronous (i.e., longer) breeding seasons. The potential to regulate breeding season length depending on environmental conditions may entail a high resilience of the population reproductive patterns against environmental changes, as highlighted by the fact that almost all mature females were reproductive every year.

Animals face changing environments throughout their whole life cycles. Individuals are adapted to the changes that are regular and predictable. The most common example is seasonality in temperate zones, for which photoperiod variation over the year represents a reliable and easily accessible predictor (Bradshaw and Holzapfel 2007). Other phenomena arise with irregular and usually unpredictable patterns, such as interannual weather variability and food or prey availability (e.g., fruit mast years) related to it (Nussbaumer et al. 2018). Whereas it is known that individuals and populations may react with plastic responses (e.g., Ruf et al. 2006; Ogutu et al. 2015), inter-individual phenotypic diversity may represent a further tool to deal with such irregular and unpredictable changes on a population level (Hertel et al. 2020).

A plastic reproductive phenology is a key ecological determinant of animal population sensitivity to changing environments as it represents the time dimension-linkage between reproduction and environment (Post et al. 2008; Ogutu et al. 2015). Such plasticity takes effect on several levels (ovulation, conception, and birth) on both individuals (Canu et al. 2015) and populations (Fernández et al. 2020). However, it is generally constrained by the reliance on rigid reproductive cues (i.e., photoperiod variations throughout the year, Bradshaw and Holzapfel 2007) that do not depend on the environmental conditions. Most ungulate populations, or at least those living at latitudes with clear day length variations throughout the year, typically show a tight reliance on photoperiod to trigger their reproduction (Zerbe et al. 2012). Nevertheless, evidence that favorable environmental and nutritional conditions facilitate a slightly earlier reproduction was frequently reported even in species whose predominant cue is photoperiodism (McGinnes and Downing 1977; Hamilton and Blaxter 1980; Flydal and Reimers 2002; Wolcott et al. 2015). Thus, a certain degree of plasticity in the reproductive timing (hereafter RT, always referred to the population level) seems to be quite spread among ungulate species and this can be expected to produce temporal displacements of breeding seasons among different years. In this context, the most plastic species have a reproductive output which is less constrained by environment, as they can respond to negative conditions by delaying the breeding season (Servanty et al. 2009).

The phenotypic diversity of reproductive phenology within a population (namely, “reproductive synchrony,” Findlay and Cooke 1982, hereafter RS) directly affects breeding (and, consequently, birth) season length. Higher RS (i.e., shorter breeding seasons) was observed in ungulate species and populations living in more seasonal and constant environments (English et al. 2012; Zerbe et al. 2012), relying on more specialist foraging strategies (English et al. 2012), showing gregarious habits associated with precocial young (Sinclair et al. 2000) and an even, rather than female-biased, sex ratio of adults (Milner et al. 2007). In a number of equatorial savanna ungulates, a substantial interannual RS variability in response to environmental conditions was reported, with longer breeding seasons observed during drought years (Ogutu et al. 2010, 2014). This phenomenon comes as no surprise in species mainly relying on environmental cues (i.e., rainfall patterns) to time their reproduction through a nutritional status mediation (Ogutu et al. 2015). Conversely, environment-driven interannual RS variability in ungulates of temperate regions (i.e., relying on photoperiod variations, Zerbe et al. 2012) is not obvious and so far has never been reported. On the one hand, as photoperiodism follows genetic heritability (Bradshaw and Holzapfel 2007; Zerbe et al. 2012), we may expect RS degree to remain substantially constant under different environmental conditions, at least assuming that they homogeneously affect all individuals. In this respect, Zerbe et al. (2012) reported unaltered RS between wild ungulates and those kept in captive conditions with high resources availability. On the other hand, resource-poor years may provoke a higher inter-individual variability in the time needed to achieve the nutritional condition required to reproduce and ultimately reduce RS.

The simpler method to investigate the variability of both RT and RS on a population level is to compare the temporal occurrence and duration of an adequate number of breeding seasons with one or more environmental variables (Ogutu et al. 2010, 2014; Fernández et al. 2020). Unfortunately, this approach requires the condensation of large datasets into 1 observation per year, with a substantial loss of statistical power. To overcome this limitation, analytical strategies aimed at evaluating the temporal variability of the individual reproductive status with respect to certain environmental conditions should be applied. A further constraint for specific investigations of RS variability in response to environmental changes is the typically short breeding season of mammal populations inhabiting temperate regions (Garel et al. 2009; Mason et al. 2011). We thus chose wild boar (Sus scrofa) as a model species because it presents the rare condition of living in temperate regions (i.e., in highly seasonal environments) and, at the same time, showing relatively long breeding seasons (Santos et al. 2006; Canu et al. 2015). The reproductive output of this species was widely investigated thanks to the large amount of data regarding culled individuals provided by hunting activities (e.g., Servanty et al. 2009; Fonseca et al. 2011; Canu et al. 2015; Bergqvist et al. 2018; Touzot et al. 2020). A high degree of individual plasticity was reported for several reproductive parameters of wild boar females, including their reproductive phenology, which tends to be anticipated in response to good environmental conditions (e.g., Servanty et al. 2009; Canu et al. 2015). Nevertheless, so far, the relationship between environmental drivers and population RT and RS has never been evaluated.

In this study, we aimed to evaluate age-specific wild boar population responses to such environmental factors as weather and resources availability in terms of both RT (anticipated or delayed breeding seasons) and RS (longer or shorter breeding seasons). In so doing, we aimed to determine whether:

wild boar shows an interannual variability of both population RT and, though inhabiting temperate regions, RS;

such interannual variability is the result of modifications of the overall individual likelihood of ovulating and getting pregnant, which in turn is affected by a number of environmental factors directly or indirectly related to resources availability; and

such environmental factors influence the population RT and RS.

Materials and Methods

Study area

We collected data in a mountainous area of 13,400 ha in Central Italy (Northern Apennines, Italy, 43° 48′ N, 11° 49′ E), which includes 2,700 ha of protected area (Oasi Alpe di Catenaia). Lowest and highest altitudes reach 330 and 1,414 m above the sea level, respectively. The climate is temperate continental with a marked seasonality. A mean temperature of 18.7 °C and a daily precipitation of 1.73 mm are recorded in summer, whereas winters are cold (mean temperature of 1.2 °C) and rainy (daily precipitation of 3.55 mm). Snowfalls are sporadic in winter and can also occasionally occur in spring. Mixed deciduous woods are the prevailing habitat category (67% of the total surface) and are mainly composed of Turkey oak Quercus cerris, beech Fagus sylvatica, and chestnut Castanea sativa. Agricultural crops (16%), mixed open-shrubs areas (10%), and conifer woods (7%) cover the rest of the surface. In the surroundings of the protected area, wild boar is unselectively hunted in drive hunts by teams of 25–50 people. During the study period, drive hunting was generally permitted 3 times a week from September to January, with an average of 58.3 hunting days per year. As a yearly average of 6.4 wild boar/km2 was harvested, the population underwent a high, but relatively constant, hunting pressure (Merli et al. 2017).

Data collection

We collected and examined reproductive traits of 2,500 female wild boars culled from 1 September to 31 January during 10 consecutive hunting seasons (2006–2016). Culling date and live body mass were recorded for each individual. In so doing, we included the reproductive trait mass, though it accounted only for a negligible percentage of female live body mass (Brogi et al. 2021). All females were aged on the basis of their tooth eruption and abrasion (Briedermann 1990) and assigned to one of the following age classes: juvenile (< 1 year), subadult (between 1 and 2 years), and adult (> 2 years). In order to determine their reproductive status, we dissected ovaries and uterus of each female to check for the presence of corpora lutea and embryos/fetuses, respectively. Corpora lutea were used as a sign that ovulation occurred, whereas embryos and fetuses as a sign of ongoing pregnancy (e.g., Malmsten et al. 2017a). Over 823 culled juvenile females, only 30 ovulated and 3 pregnant individuals were identified. We thus decided to exclude the individuals belonging to this class from our analysis. The Regional Hydrological Service of Tuscany kindly provided weather data (average temperature and rain) daily recorded in a weather station located inside our study area (43° 42′ N, 11° 55′ E). We obtained local data on yearly seed productivity of beech, chestnut, and Turkey oak measured inside the Oasi Alpe di Catenaia from an online database (Chianucci et al. 2019) and used it as a measure of food availability.

Data analysis

Step 1: ovulation and pregnancy heterogeneity among years and classes

In order to assess interannual heterogeneity in ovulation and pregnancy patterns, we modeled both individual likelihood of ovulating and getting pregnant throughout the sampling period. We divided our dataset into 2 sub-datasets corresponding to subadult and adult females. By means of the glm() function of the stats package (R version 4.0.3, R Development Core Team, 2015) we ran a Generalized Linear Model (GLM) with a binomial distribution, with the individual reproductive states (ovulated or pregnant, alternatively) as binary-dependent variables, separately for the subadult and the adult female sub-datasets. The binary variable “ovulated” took the value 1 whenever at least 1 corpus luteum, embryo, or fetus was detected, and 0 otherwise; the binary variable “pregnant” took the value 1 whenever at least 1 embryo or fetus was detected, and 0 otherwise. In so doing, we built a total of 4 models, hereafter called “1S-ov” (model of ovulation in subadults), “1A-ov” (ovulation in adults), “1S-pr” (pregnancy in subadults), and “1A-pr” (pregnancy in adults). We included in all models the interaction term between the standardized culling date (expressed as days from 1 September) and the hunting season (categorical) as the only predictor to check for interannual variations in the effect of the standardized date. The hunting season was expressed as a sequential number from 1 (referring to the 2006–2007 hunting season) to 10 (2015–2016).

To check for statistical differences among age classes in ovulation RT and RS, we used the models 1S-ov and 1A-ov to predict yearly dates of onset, middle point, and end of ovulation on a population level for each monitored hunting season. The dates in which the proportion of ovulated females reached 0.025, 0.5, and 0.975 were used as onset, middle point (inflection point of the curve) and end date, respectively. In so doing, we included 95% of the predicted ovulation events between the onset and end dates. In order to test whether ovulation was significantly anticipated in a certain age class in respect to the other, we performed a paired t-test (t.test() function of the R package stats), which compared subadult and adult female middle point dates for each hunting season. Moreover, to check for inter-class differences in ovulation RS, we measured the duration of ovulation seasons (1 per year) as the number of days from the onset to the end dates in both subadult and adult females. We then calculated the average duration of the ovulation season and its associated variance, separately in subadult and adult females. Finally, we ran a 2 samples t-test for summary data implemented by the tsum.test() function (R package BSDA). The whole procedure was then exactly replicated on pregnancy RT and RS by using yearly predictions of the models 1S-pr and 1A-pr.

Step 2: factors influencing individual reproductive status

After the analysis to test potential differences among seasons within age classes, we aimed to identify internal and external factors which influenced ovulation and pregnancy ratios. We modeled the individual likelihood of ovulating and, alternatively, of getting pregnant by means of 4 GLMs with a binomial distribution (ovulation in subadults, ovulation in adults, pregnancy in subadults, and pregnancy in adults). The standardized culling date (days from 1 September) was used as predictor to consider photoperiod-mediated seasonal variations of the individual reproductive status. We also included such internal factors as individual age (months) and live body mass (kg) as predictors. Among external factors, 4 season average temperature and rain precipitation calculated on a yearly basis were used as predictors to account for the potential effect of weather. Because all individuals were culled between September of year x and January of year x + 1, winter weather variables were averaged from December of year x−1 to February of year x, spring ones from March to May of year x, summer ones from June to August of year x, and autumn ones from September to November of year x. Moreover, we used current year seed productivity of Turkey oak, beech, and chestnut (t/ha) measured on a yearly basis to check for potential effects of food availability on ovulation and pregnancy patterns. To summarize the influence of the 3 deciduous species in a single variable, we included a further global forest productivity index in the models, which we calculated following the protocol described by Bisi et al. (2018). Finally, we calculated the yearly average number of adult males per female as the number of culled adult males (>3 years; Brogi et al. 2021) divided by the total number of adult and subadult females. We added this yearly variable as a predictor within our models to take into account the potential effects of reproductive male relative abundance on female reproductive status (Milner et al. 2007). We recognize that, by measuring adult male availability on the basis of culling data, we may obtain an unreliable approximation of the real population structure. However, in this study, we were only interested in the variation of male availability throughout different years.

Separately for each sub-dataset, we screened all available predictors for collinearity and multicollinearity by means of a Pearson correlation matrix (rp) and the variance inflation factor (VIF), setting thresholds to rp = ± 0.7 and VIF = 3, respectively (Zuur et al. 2009). Weather variables of the same season (particularly spring and autumn) were the most recurring pairs of collinear variables. We performed a random forest calculation (random.Forest package) to rank all predictors on the basis of their potential to explain the dependent variable (Breiman 2001). The worst predictor variable of each collinearity and multicollinearity condition was dropped until all rp and VIF were below the corresponding thresholds. Finally, we included the remaining predictor variables in a full GLM and used the dredge() function (MuMln package) to run a set of models with all possible combinations of predictor variables. We followed the minimum Akaike's information criterion (AIC) and selected the most parsimonious in terms of number of predictors among groups of models with ΔAIC < 2 (Symonds and Moussalli 2011), identifying the 4 best models: “2S-ov” (ovulation in subadults), “2A-ov” (ovulation in adults), “2S-pr” (pregnancy in subadults), and “2A-pr” (pregnancy in adults).

Step 3: effects of internal and external factors on RT and RS

In the last step of our analysis, we aimed to assess whether the factors affecting ovulating and pregnant female ratios (Step 2 of our analysis) may also provoke modifications in ovulation and pregnancy temporal patterns. We thus built 4 further GLMs, 1 for each combination of dependent variables and age classes (e.g., ovulation in subadult females). We included the set of predictor variables of the corresponding best model selected in Step 2 (e.g., 2S-ov) and added them all their interactions with the standardized date. Following the same protocol described in Step 2, we then screened this enlarged sets of predictor variables for collinearity and multicollinearity, ran full models, and processed them with dredge() function to finally select 4 new best GLMs including single and interaction terms: “3S-ov,” “3A-ov,” “3S-pr,” and “3A-pr.”

Results

Step 1: ovulation and pregnancy heterogeneity among years and classes

Interannual ovulation and pregnancy patterns predicted by 1S-ov, 1A-ov, 1S-pr, and 1A-pr are summarized in Figure 1. A marked interannual heterogeneity affected temporal patterns of both reproductive statuses considered, although the predicted portion of females achieving ovulation or pregnancy within the sampling period was always equal or close to 1 in both age classes. A number of reproductive seasons were relatively early and short (hunting seasons 2, 5, and 7), whereas others showed either a late onset (3, 6, and 10) or a longer duration (1 and 9). Likewise, the temporal distance between ovulation and pregnancy curves varied among the years, with the minimum value observed in hunting season 2 and the maximum in 5 and 8. Finally, subadult and adult females showed completely overlapped reproductive patterns in a number of hunting seasons (2, 6, and 7) and markedly divergent in other ones (3 and 4).

Figure 1.

Ovulation (continuous lines) and pregnancy (dashed lines) patterns of subadult (red) and adult (blue) females throughout 10 hunting seasons in Northern Apennines, Italy. Values were predicted by 4 GLMs with the interaction between date and hunting season as the only predictor variable (see the text for more details). Color-shaded areas represent 95% confidence intervals.

On average, the date when the proportion of ovulated females reached 0.5 corresponded to 82.46 (21 November) ± 14.67 (mean ± SD) and 83.77 (23 November) ± 13.60 days from 1 September in subadults and adults, respectively, without a statistically significant difference between the 2 age classes (t = −0.55, P-value = 0.60). A similar result was detected for pregnancy, as subadult females reached the middle point at 109.60 (19 December) ± 14.82 days from 1 September and adult females at 115.61 (25 December) ± 17.88 days from 1 September, with the paired t-test returning a non-significant difference (t = −1.70, P = 0.12). Conversely, the duration of the ovulation season (a measure of RS) was shorter in subadult (96.54 ± 9.46 days) than in adult females (114.00 ± 10.85 days) and this difference was statistically significant (t = −3.84, P = 0.0012). As 95% of subadult females got pregnant in 94.20 ± 10.65 days, whereas adult females in 121.13 ± 16.01 days from the onset, pregnancy season duration was significantly shorter in subadult females (t = −4.43, P = 0.0004).

Step 2: factors influencing individual reproductive status

Predictor variable sets included in the best model for the 4 GLMs explaining the individual likelihood of ovulating and getting pregnant are summarized in Table 1, whereas those selected for random forest analysis and dredge are summarized in Supplementary Table S1. Standardized date and average spring temperature were included in all 4 best GLMs and positively affected both ovulation and pregnancy rates in both age classes. Individual body mass only increased the likelihood of subadult females ovulating, whereas its positive effect on pregnancy ratio concerned both age classes. As for food availability, at least 1 predictor reflecting seed productivity was included in each best GLM. The relative abundance of adult males was not selected for any best GLM.

Table 1.

Sets of explanatory variables included in the best GLM on the individual likelihood of: subadult females ovulating (2S-ov); adult females ovulating (2A-ov); subadult females getting pregnant (2S-pr); and adult females getting pregnant (2A-pr).

Model	Sub-dataset	Reproductive state	Best model formula
2S-ov	Subadult females	Ovulation	Ovulated ∼ standardized date + body mass + spring temperature + autumn rain + global productivity index
2A-ov	Adult females	Ovulation	Ovulated ∼ standardized date + spring temperature + summer rain + autumn rain + chestnut productivity + beech productivity
2S-pr	Subadult females	Pregnancy	Pregnant ∼ standardized date + body mass + spring temperature + summer rain + chestnut productivity + global productivity index
2A-pr	Adult females	Pregnancy	Pregnant ∼ standardized date + body mass + spring temperature + chestnut productivity + beech productivity

Standardized culling date, culling date expressed as days from 1 September; body mass, individual body mass (kg); season x temperature, average environmental temperature recorded during the season x; season x rain, average daily rain precipitation recorded during the season x; productivity of species y, mast productivity of the tree species y during the current year expressed as t/ha; global productivity index, index summarizing all tree species productivity during the current year (see the text for more details).

Step 3: effects of internal and external factors on RT and RS

The model subadult female ovulation (3S-ov) included individual body mass, spring average temperature, and autumn rain as single variables in addition to the 2 interaction terms composed of [global productivity index: date] and [spring temperature: date], all showing a positive effect on the dependent variable (Supplementary Table S2a). The increase of global productivity index did not cause a substantial displacement of the ovulation onset. However, it was related to a marked shortening of the ovulation season (higher RS) from ∼110 days predicted for low productive years to ∼70 days predicted for highly productive years (Figure 2A). Likewise, in years with higher average spring temperature, subadult female ovulation season was shorter, though with a markedly anticipated RT (Figure 2B).

Figure 2.

Predicted effect of the interaction between environmental variables and the standardized date on the proportion of: ovulating subadult females (A and B), ovulating adult females (C and D), pregnant subadult females (E and F), and pregnant adult females (G and H), expressed by the chromatic scale (white = low; black = high). Blue lines represent 0.025 (ovulation and pregnancy season onset), 0.25, 0.5, 0.75, and 0.975 (ovulation and pregnancy season end) isolines. Spring temperature: average air temperature of previous spring (°C); Global productivity index: mast tree global productivity index (see the text for more details); Beech productivity: beechnut productivity (t/ha).

For adult female ovulation patterns, model 3A-ov included spring average temperature, autumn rain, and chestnut productivity as single variables and [beech productivity: date] and [spring temperature: date] as interaction terms (Supplementary Table S2b). Beech productivity only accounted for a slight shortening of the ovulation season (higher RS), with no effect on the timing of its onset (Figure 2C). Conversely, warmer spring temperatures were associated to both anticipated RT and higher RS of ovulation seasons (Figure 2D).

The model 3S-pr, which explained subadult female pregnancy patterns, included individual body mass and chestnut productivity as single variables in addition to the same interaction terms selected for ovulation patterns of the same age class, that is, [global productivity index: date] and [spring temperature: date]. When seed productivity was higher, subadult female pregnancy showed an anticipated RT and a higher RS (Figure 2E). A similar pattern was observed for average spring temperature, though with a stronger effect in anticipating pregnancy RT (Figure 2F).

The model 3A-pr, which accounted for adult female pregnancy patterns, included individual body mass and chestnut productivity as single predictor variables in addition to the same interaction terms selected for ovulation patterns of the same age class, that is, [beech productivity: date] and [spring temperature: date]. Their effects on RT and RS were similar to those shown on adult female ovulation, though isolines showed an overall delay (Figures 2G,H).

Discussion

We showed that, in an ungulate species inhabiting temperate latitudes, breeding seasons can change in timing and duration, depending on environmental conditions. Both population RT and RS widely varied among different years and our analytical approach enabled to properly evaluate their dependence on the environment. These phenomena were essentially due to the individual tendency to reproduce even when a harsh environment made the investment risky in terms of offspring survival. Such population-level features likely entail a high resilience of the population reproductive patterns against ecological perturbations and environmental changes as confirmed by the extremely high average likelihood of females ovulating or getting pregnant by the end of the reproductive season in every sampling year.

We observed a high temporal heterogeneity among yearly reproductive patterns (Figure 1). However, in accordance with Servanty et al. (2009), the model described in Step 1 predicted an average individual likelihood of ovulating which reached values close to 1 before 31 January every year and in both age classes considered. Pregnancy followed similar patterns, thus proving that ovulation rates represent a good wild boar pregnancy proxy. Interannual pregnancy delay variability in respect to ovulation was likely the effect of a variable proportion of ovulated females failing to get pregnant. However, thanks to their ability to repeat the estrus (Henry 1968; Barrett 1978; Macchi et al. 2010), all female wild boar (subadult and adult) were predicted to achieve pregnancy even in the years with the highest delays (e.g., hunting seasons 5 and 8). Although minor reproductive events may occur all year round in other wild boar populations (relying on artificial food, Macchi et al. 2010; Malmsten et al. 2017b; Bergqvist et al. 2018), our results showed that, for adult and subadult females, an actual breeding season existed and was included within our sampling period. The minimal number of reproductive juvenile females detected in our study (823 culled juvenile females, 30 ovulated, and 3 pregnant) may be a sign of their contribution to reproduction being negligible or the consequence of the 5 months sampling period duration being insufficient to detect juvenile reproduction, which has been shown to occasionally occur in other wild boar populations (Šprem et al. 2016; Gamelon et al. 2017). Collecting data all year round (possible in cases of wild boar hunting being performed during the whole year) would be necessary to properly evaluate the reproductive contribution of different classes of females outside the core reproductive period, but it is worth noting that this was not the objective of this study.

Subadult females were significantly more synchronous than adults, likely on account of an overall higher homogeneity of their individual conditions. Differently from the older class, all subadult females belonged to the same cohort and most of them were at their first reproductive attempt (as confirmed by the almost null reproductive rate observed in juvenile females) and had not to cope with previous parental reproductive costs. Conversely, adult females had different ages and might have coped with different costs related to their previous reproduction (Hamel et al. 2010).

The fact that the average likelihood of ovulating and getting pregnant reached values close to 1 within our sampling period enabled an unambiguous interpretation of the Steps 2 and 3 analyses: the effects of the environmental factors identified only either anticipated or delayed changes of the reproductive status, without truly affecting the individual likelihood of ovulating and getting pregnant by the end of the reproductive season. This evidence helps to understand environmental influence on female wild boar reproductive status, which so far was widely investigated by focusing on the overall proportion of reproductive females (Fonseca et al. 2011; Bergqvist et al. 2018; Touzot et al. 2020) and seldom considering the temporal dimension (Servanty et al. 2009). In this context, a yearly proportion of reproductive females estimated without taking into account culling dates is prone to be substantially underestimated. In fact, females culled early in the hunting season with no sign of ongoing ovulation or pregnancy and considered “not reproductive” (Fonseca et al. 2011; Bergqvist et al. 2018; Touzot et al. 2020) should rather be considered “not reproductive yet.”

The influence of the standardized date was included in all the best models selected in Steps 2 and 3 (as single predictor and in interaction with environmental variables, respectively). Thus, it is suggested that photoperiodism still constrained wild boar RT, though its influence was not so strong if compared with that exerted over most ungulates inhabiting temperate regions. This evidence places wild boar at an intermediate position along an ideal continuum between temperate ungulates (which rigidly rely on photoperiodism to time their reproduction, with minor environmental influence, Zerbe et al. 2012) and equatorial, seasonal breeding ungulates (whose reproductive phenology mainly relies on environmental cues, Ogutu et al. 2015).

The approach adopted to build Step 3 models enabled to evaluate ovulation and pregnancy temporal patterns of the population in respect to the environment, that is, to monitor the breeding season temporal onset, progress, and duration at varying environmental conditions. Ovulation and pregnancy RTs were substantially anticipated under good environmental conditions (i.e., higher resources availability and warmer spring temperatures) in both age classes (Figure 2), thus showing the high degree of ecological plasticity of wild boar reproductive phenology. The physiological phenomenon was likely mediated by individual nutritional conditions (McGinnes and Downing 1977; Hamilton and Blaxter 1980; Flydal and Reimers 2002; Wolcott et al. 2015), which were directly improved either by resource abundance or by the advanced vegetation growth due to high spring temperatures.

The possibility to either plastically anticipate or delay breeding seasons maximizes population reproductive outcomes under optimal conditions, whereas increasing its resilience against ecological perturbations. During favorable years, anticipated breeding seasons produce earlier births, which are known to increase offspring survival in ungulates (Côté and Festa-Bianchet 2001). In the case of wild boar, earlier births may directly reduce the young mortality caused by red fox (Vulpes vulpes) predation (Bassi et al. 2012) by producing a beneficial mismatch between the time when piglets are of vulnerable size and the time when fox food requirement is most intense (young raising, from May onwards in Southern Europe, Cavallini and Santini 1995). The potential to plastically anticipate breeding seasons may result extremely beneficial also when facing global change by softening or even preventing mismatches between births and the most favorable nutritional conditions for offspring. In this respect, wild boar may represent an exceptional case of a species “pre-adapted” to global change, as already suggested (Vetter et al. 2015; Touzot et al. 2020). Conversely, when less resources are available, a delayed breeding season gives individuals more time to get the nutritional condition needed to reproduce. In so doing, a higher proportion of mature individuals can achieve reproduction at the cost of an increased offspring mortality. The high hunting pressure may have increased the advantage of such a risky investment, as individuals counting on a short life expectancy have to exploit every reproductive opportunity to maximize their fitness (Festa-Bianchet 2003). We observed no relationship between the number of culled adult males per female and ovulation and pregnancy temporal patterns. This result is surprising in a heavily hunted population (i.e., subject to adult male scarcity, Fernandez-Llario and Mateos-Quesada 2003; Toïgo et al. 2008) and appears in contrast with the results obtained for other ungulate species (Milner et al. 2007). Nonetheless, it is consistent with the findings proposed by Diefenbach et al. (2019) on white tailed deer (Odocoileus virginianus) as well as with Brogi et al.’s (2021) hypothesis regarding the flexible reproductive involvement of subadult male wild boar. As we did not consider other population traits, such as density and structure, further investigations are needed to evaluate their potential effect on wild boar temporal reproductive patterns.

A number of environmental factors in interaction with the standardized date were included as predictors in Step 3 best models, thus showing that good environmental conditions (higher spring temperatures, higher food availability) enhanced RS and ultimately led to shorter breeding seasons (Figure 2). We thus showed that, as previously reported only for equatorial ungulates (Ogutu et al. 2010, 2014), photoperiodic species inhabiting temperate regions also have the potential to adjust breeding season length depending on environmental conditions. In the monitored population, RS was enhanced by higher spring temperatures in both age classes and by global seed and beechnut productivity in subadult and adult females, respectively. Breeding seasons following hot springs were ∼40% shorter in respect to those following cold springs in both age classes. Global seed productivity had a similar impact (shortening of ∼36%) on subadult female ovulation seasons, whereas years with a high beechnut productivity reduced adult female ovulation season length of ∼20% in respect to less productive ones. These environmental factors likely induced a plastic anticipation of individual reproductive phenology but heterogeneously affected each individual. Conversely, only the average population RT would have been modified, with no effect on inter-individual differences and, therefore, on RS (as in the case of other environmental factors included as single predictors in Step 3 best models). We can suppose that, when the main food resources were more abundant, all females reached the threshold nutritional condition needed to reproduce early and achieved ovulation as soon as their internal photoperiodism enabled them to. This optimal nutritional condition induced a quite homogenous distribution of ovulation within the population. Conversely, in case of low resource availability, the pre-existing variability of individual conditions would be unaltered or even enhanced. For instance, foraging strategies would be more diversified, with a number of individuals either being able to outcompete the others for the scarce resources available or better exploiting secondary food items. The whole breeding season RT would be delayed (as observed, for example, when the global productivity index was low), though a number of individuals would be less affected than others by resource scarceness and still be able to pursue an early reproduction, thus inducing a substantial RS reduction. In this context, spring temperatures may have acted as a proxy of the vegetation growth season and regulated abundance and temporal occurrence of food resources other than mast seeds.

The possibility to regulate RS in respect to the environmental conditions may provide several advantages to the population reproductive outcomes. In particular, birthdates may be highly concentrated when, during the mating season, environmental conditions are good (and likely induced a high nutritional condition of females). When favored by resource availability, the advantageous (Côté and Festa-Bianchet 2001) phenotypic trait of early reproduction may thus be evenly distributed within the population. We can hypothesize that a higher birth synchrony may also reduce predation risk by saturating the number of newborns that predators (wolves, Canis lupus, and foxes in the monitored study area, Bassi et al. 2012) can catch per time unit (dilution effect, Darling 1938). Conversely, under suboptimal environmental conditions, the enhanced phenotypic diversity showed by the population reproductive phenology may produce more scattered birthdates. This may result in a more efficient resource partitioning among individuals that are raising young (Ims 1990). However, more scattered birthdates amount to a population trait and therefore may not be shaped directly by evolution and, as explained above, rather seems the consequence of the combination of individual adaptive features.

We provided the first evidence of breeding season length adjustment depending on environmental conditions in a species living in temperate regions and relying on photoperiodism to trigger its reproduction. This feature likely represents a key factor for wild boar renowned ecological plasticity and ultimately contributes to its high success and worldwide spread (Massei et al. 2015; Markov et al. 2019).

Click here for additional data file.