Data collection
We measured 158 pine marten skulls collected between 1903 and 2020 originating from two sites: (1) The Białowieża Forest (NE Poland and W Belarus; N=93) and (2) Silesia (SW Poland; N=65; Fig. S1). The Białowieża Forest is the best-preserved forest area in the European lowland (Jaroszewicz et al., 2019), while Silesia includes Bory Dolnośląskie, which is one of the largest complexes of forests in Europe. The specimens were collected in the Museum für Naturkunde Berlin, the Zoological Institute of the Russian Academy of Sciences in Saint Petersburg, the Museum of Natural History, Wrocław University, and the Zoological Collection of the Mammal Research Institute Polish Academy of Sciences (MRI PAS). For each skull, measurements of 11 traits were made with an electronic calliper (with an accuracy of up to 0.01 mm): condylobasal length, maxillary tooth-row length, palatal length, zygomatic width, interorbital constriction width, postorbital constriction width, mastoid width, mandibular tooth-row length, mandible height, mandible length, and braincase height, as per Reig and Ruprecht (1989) (Fig. 1). All measurements were carried out by one observer in order to avoid measurement errors. We excluded young individuals with incomplete fusion of the cranium bones and juvenile teeth.
To analyse temporal sex ratio changes, in addition to the data set used for skull measurements, we also used data from sexed individuals, mostly roadkill, collected in the MRI PAS, whose skulls were damaged or not collected, which increased the dataset to 205 individuals (72 females and 133 males) between the years 1959 and 2020. The probability of being male was analysed for the last 61 years because of the collected data was more evenly distributed over the years analysed without large time gaps.
Climate data
We obtained monthly temperature (temp), precipitation (prec) and snow cover (snow) and averaged them for the winter period (from November of the previous year to February of the skull collection year; tempWINT, snowWINT) and/or summer (May–August of the skull collection year; tempSUMM, precSUMM) for each study site separately. We also used average carbon mass flux from the atmosphere due to Net Primary Production (NPP) in the months of peak vegetation, representative for Europe (May–August; nppSUMM; Hicke et al., 2002), to assess intensity of the influence of NPP on skull variation. Climate data were obtained from the Land Surface, Snow and Soil Moisture Model Intercomparison Project (LS3MIP; van den Hurk et al., 2016) developed under the sixth phase of the Coupled Model Intercomparison Project (CMIP6; Eyring et al., 2016). To obtain variables from across the entire study period, we used two LS3MIP models: Land-Hist and LFMIP-rmLC. The Land-Hist historical simulation covers the period 1850 to 2015 and provides a comprehensive assessment of land surface, snow and soil moisture–climate feedbacks. The LFMIP-rmLC model was designed to simulate the impact of a future climate on land surface and was extended to 2100 (van den Hurk et al., 2016). To cover the study period 1900–2020, we obtained climatic variables from Land-Hist from 1900 to 1970 and from LFMIP-rmLC from 1971 to 2020. Both data sets are highly correlated in the overlapping period of 1970–2015 (r = 0.93, p < 0.001); therefore, it was possible to combine them into one data set.
Statistical analysis
All analyses were conducted in R version 4.0.3 (R Core Team, 2020). We tested body size variation predicted by crania and mandible variables across a long-term trend (117 years), using a generalized additive model (GAM; package ‘mgcv’ version 1.8; Wood, 2017). To reduce the number of variables for skull measurements entering the model, which are often highly correlated, we performed a principal component analysis (PCA) of skull measurements using the PCA function (‘FactoMineR’ package; Lê 2008). Due to damage to some skulls and missing measurements (depending on the type of measurement, damage was from 4% to 17%), we conducted a PCA based on an estimated covariance matrix obtained by function imputePCA (‘missMDA’ package) which uses regularised iterative PCA to impute and estimate the covariance matrix of incomplete data (Josse and Husson, 2016; Nassiri et al., 2018). The GAMs were conducted separately for condylobasal length, first principal component (PC1), and postorbital constriction width, with year, sex, and site (two-level factor: Białowieża Forest and Silesia) as explanatory variables. We added sex and site as explanatory variables to control the effects of sexual dimorphism and differences in marten size between sites. We used condylobasal length separately, as this measurement is most often used to predict body size, to which it is highly correlated (Ochocinska and Taylor, 2003; Yom-Tov et al., 2003), and postorbital constriction width, which reflect chewing and biting force, as only this measurement was not correlated with PC1 (see Results). Isotropic thin plate regression splines with knot-based approximations were used for temporal smoothing. Models were specified using Gaussian distribution, an identity link function, and an REML approach.
In the next step, we tested crania and mandible variation in relation to climatic condition and primary productivity using four climatic variables and one primary productivity variable: tempWINT, tempSUMM, snowWINT, precSUMM, and nppSUMM. We were not able to obtain data about the years of birth of the martens, only their year of death. Pine martens can survive up to 11 years with an average lifespan of 3 years (Marchesi, 1989); therefore, climatic variables in the year of death do not represent the conditions during life, which most affect skull size (Geist, 1987; Searcy et al., 2004). In order to explain the effect of temporal changes in climate and NPP on body size, we used GAM-predicted climatic variable values by year as proxies of temporal changes of those variables (Fig. 2). Thus, we obtained smoothed values for each climatic variables and NPP corresponding to the annual temporal trend of each variable separately for each site (Fig. 2).
To test the influence of climatic conditions and NPP on skull size, we first used a GAM (as previously) but the relations were linear. Consequently, we used a linear model (LM) with climatic, NPP, sex, and site as explanatory variables. Variance inflation factors (VIFs) indicated collinearity among all climatic explanatory variables (Table S1); thus, models were conducted separately for each climatic and NPP explanatory variable. A set of five LM models, each with one relevant climatic or NPP variable, were performed for condylobasal length, PC1, and postorbital constriction width and their performance compared by Akaike information criterion (AIC) values. Models within ΔAIC < 2 were assumed to be equally parsimonious (Harrison et al., 2018). To account for climatic conditions during body growth, the influence of predicted climatic variable values on condylobasal length, PC1 and postorbital constriction width were considered with different time-lags (year t, year t-1, year t-2, year t-3) using LM and then compared by their AIC values.
Temporal trends in sex ratio were examined using logistic regression (GLM with binomial error distribution and logit link), with sex as the dependent variable and year and site (two-levels factor: Białowieża Forest or Silesia) as explanatory variables. We tested whether the probability that a collected individual would be a male (1) or a female (0) was predicted by time. Next, we tested probability of being male in relation to climatic and NPP variables, building five GLMs with sex as the dependent variable and site and one of the climatic or NPP variables as explanatory variables and then comparing by AIC values of those models.