Different sets of traits explain abundance and distribution patterns of European plants at different spatial scales

Aim: Plant functional traits summarize the main variability in plant form and function across taxa and biomes. We assess whether geographic range size, climatic niche size, and local abundance of plants can be predicted by sets of traits (trait syndromes) or are driven by single traits. Location: Eurasia. Methods: Species distribution maps were extracted from the Chorological Database Halle to derive information on the geographic range size and climatic niche size for 456 herbaceous, dwarf shrub and shrub species. We

local abundance depends on factors operating at the local scale of species assemblages, such as habitat suitability, the local combination of environmental conditions, and biotic interactions (Peterson et al., 2011;Staniczenko et al., 2017). Under the assumption that species' functional traits reflect the mechanisms through which species respond to abiotic and biotic conditions to maximize their fitness, these traits are expected to predict both broad-scale distribution and local abundances (Suding et al., 2008;Heino & Tolonen, 2018).
Species can be rare or common (i.e. less or more abundant) within a local plant community. Similarly, some species have restricted distribution ranges while others are geographically widely distributed (Rabinowitz, 1981;Gurevitch et al., 2002;Enquist et al., 2019). It has been observed that species with larger geographic range sizes tend to have broader environmental tolerances (i.e. broader climatic niches), while geographically narrowly distributed species are also more likely to be narrowly distributed in climatic space (Slatyer et al., 2013;Sporbert et al., 2020). A positive relationship between climatic niche size and geographic range size across species thus seems to be a general macroecological pattern (Gaston, 2000;Slatyer et al., 2013;Cardillo et al., 2019). A species' local abundance results from population growth and demographical performance (Peterson et al., 2011). Within the geographic distribution range of a species, its local abundance at the community level is often highly variable. At the local scale, species abundance values are frequently used as descriptors of species performance and are an important characteristic of the composition of herbaceous plant communities (Kent and Coker, 1992;Chiarucci et al., 1999). In general, locally rare species tend to have a sparse cover in plant communities (Murray & Lepschi, 2004). Thus, potentially, local cover could also be considered a proxy for local rarity or commonness. However, local cover is in general low at most sites and high at only a few sites across a species' distribution range (Murphy et al., 2006). In contrast to "everywhere sparse" species, these "somewhere abundant" species are reflected in right-skewed species abundance distributions, a common pattern in plant community ecology (McNellie et al., 2019). This skewness in local abundance might be caused by the distribution of optimal ecological conditions, and thus, might be causally linked to functional traits. As mean abundance across the species range itself does not capture the full variability of skewed frequency distributions, it should be considered together with the skewness of a species' cover value across its distribution range as proxies for rarity or commonness.
Functional traits have been used as proxies for species' dispersal abilities (Greene and Johnson 1993;Thompson et al., 2011), environmental tolerances (Loehle, 1998;Bohner & Diez, 2020) or competitiveness (Kunstler et al., 2016). Specific functional traits have been linked to commonness and rarity on both local and large scales (see Table 1). For example, studies have found plant height, used as a proxy for competitive ability, to be positively correlated with range size, with taller species more widespread than shorter ones (Lavergne et al., 2004;Kolb et al., 2006). Similarly, on the local scale, common (i.e. more abundant) species have been associated with taller stature and with other traits that are proxies for species' physiological activity and productivity, including larger specific leaf area (SLA) and higher leaf nitrogen (N) content (Grime et al., 1997;Hegde & Ellstrand, 1999;Lavergne et al., 2004;Mariotte, 2014;Lachaise et al., 2020). Nitrogen (N) and phosphorus (P) availabilities limit plant growth in most terrestrial ecosystems (Güsewell, 2004).
Low nutrient availability (e.g. phosphorus limitation) may weaken the relationship between productivity-related traits and macroclimate (Bruelheide et al., 2018). As a consequence, there might be a negative correlation between species' N:P ratio and both their local abundance and broad-scale distribution. Regarding species' persistence, locally more abundant species have been associated with perennial life cycle and clonal growth (Eriksson & Jakobsson, 1998;Kolb et al., 2006). In contrast, at large spatial scales, rare TA B L E 1 Traits used in this study, their function in the community, and their reported correlation with local abundance and broad-scale distribution being unimodal (─), positive (↑) or negative (↓) SPORBERT ET al. species have been associated with prevailing clonal growth (Kelly & Woodward, 1996) and woodiness (Oakwood et al., 1993).
Several studies have investigated the relationships linking dispersal or regeneration-related traits with species' local abundance and broad-scale distribution patterns. On the local scale, more abundant species were found to produce fewer and lighter seeds than rare species (Hedge & Ellstrand, 1999;Guo et al., 2000;Kolb et al., 2006). In contrast, at large spatial scales, geographically widespread species have been found to produce significantly more and heavier seeds than small-ranged plant species (Lavergne et al., 2004;Kolb et al., 2006;Van der Veken et al., 2007).
While some studies have found relationships between functional traits and local abundance and/or broad-scale distribution patterns, others have failed to detect a clear correlation (see Table 1). So far, the majority of studies have focused on single traits rather than on trait combinations or trait syndromes (but see Díaz et al., 2016;Guo et al., 2018) as predictors of large and local distribution patterns. However, no single trait can completely describe a species' ecological strategy (Winemiller et al., 2015;Marino et al., 2020). Rather, species' local abundance and broadscale distribution patterns might be affected by different sets of traits (Marino et al., 2020). It has been suggested that locally rare and geographically restricted plant species differ systematically from more common species in functional traits that are related to species' productivity, competitive ability, dispersal, regeneration and persistence (Murray et al., 2002). However, the different states and values of traits cannot be unconditionally combined. Díaz et al. (2004) highlighted that the functional space occupied by vascular plant species is strongly constrained by trade-offs between traits. On the one hand, the leaf economics spectrum describes a productivity-persistence trade-off and contrasts species with a set of successful trait combinations for quick returns on investments of nutrients and dry mass in leaves to species with a slower potential rate of return of more persistent leaves (Wright et al., 2004). On the other hand, the size spectrum reflects the species' life cycle, with small stature species, smaller seeds and short lifespans vs long-lived woody plants (Díaz et al., 2016; Table 1).
In this study, we aimed at unravelling the relationships between traits (single traits or trait syndromes) and species distributions at broad spatial scale and abundances at local scale. Specifically, we focused on 20 traits that are expected to respond to bioclimatic drivers and capture the essence of plant life forms and functions (Wright et al., 2004;Petchey & Gaston, 2006;Díaz et al., 2016;Bruelheide et al., 2018). We tested for these relationships across 456 European herbaceous, dwarf shrub and shrub species by interrelating existing data on functional traits with the species' (a) geographic range size; (b) climatic niche size; and (c) local abundance, which was measured as (i) mean cover from all the vegetation plots in which a species was present and (ii) skewness of cover values. We expected climatic niche size and geographic range size to be driven by the same underlying environmental factors and ecological processes (Colwell & Rangel, 2009), and therefore to be positively correlated, and to be predicted by many of the same single traits or trait syndromes (Table 1). We aimed to answer the following research questions: (a) can single plant functional traits or sets of traits (trait syndromes) best explain the local abundance including their neophytic occurrences. A list of these species can be found in Appendix S1. In total, CDH stores information on species distribution ranges for more than 17,000 vascular plant species but expert-drawn range maps were compiled for 5,583 taxa based on national and floristic databases and maps from the floristic literature (Tralau, 1969(Tralau, -1981Lundquist & Nordenstam, 1988;Lundquist, 1992;Lundquist & Jäger, 1995. These data are published as distribution range maps (Meusel et al., 1965(Meusel et al., , 1978Meusel & Jäger, 1992). We used the subset of these species that met the criteria described below. Data stored in CDH can be requested for research objectives via choro logie.biolo gie.uni-halle.de/choro/.
We aggregated species' point and polygon distribution data using a raster grid layer of 2.5 arc-min resolution, which corresponds to grid cells covering approximately 15 km 2 each across Central Europe. As a measure of range size for each species, we counted the number of grid cells occupied (approximating the area of occupancy in the geographical space).
We determined the multi-dimensional climatic space (or climatic niche) of each geographic range based on principal component analysis (PCA) of 19 bioclimatic variables from the WorldClim 2.0 database (Fick & Hijmans, 2017), also at 2.5 arc-min resolution. The resulting global background climatic space was well represented by the first two principal components, which accounted for 70.75% of the total climatic variance. The two-dimensional PCA space was rasterized into 100 × 100 PCA grid cells, considered as the background climatic niche, as explained in Appendix S2. The species' niche size was then calculated as the number of PCA grid cells occupied in the climatic space (i.e., the area of occupancy in the bioclimatic niche space; for detailed information see Appendix S2).

| Local abundance metrics in vegetation plots: mean cover and skewness of cover values
As a measure of local abundance, percentage cover values were obtained for each of the study species in 740,113 vegetation plots from | 5 of 14  Jansen & Dengler, 2008) and to four taxonomic reference lists available via the R package taxize (Chamberlain & Szöcs, 2013;R Core Team, 2018), i.e.

Journal of Vegetation Science
Encyclopedia of Life (EOL), International Plant Names Index (IPNI), Integrated Taxonomic Information Service (ITIS) and Tropicos. In cases where no exact match was found, taxon names were resolved using the Taxonomic Name Resolution Service (TNRS) and all names matched or converted from a synonym were considered accepted taxon names when probabilities were ≥95%. We merged the data for subspecies at the species level following the taxonomic hierarchy in TNRS. The selected study species occurred within at least 100 vegetation plots in the EVA dataset. Vegetation plots with a geographic location uncertainty of more than 10 km were removed prior to this selection. The median occurrence (i.e. number of vegetation plots a species occurred in) per species was 2,162 (interquartile range 846 to 5,137). Information on source databases that provided vegetation-plot data can be found in Appendix S3. Cover or cover abundance values that were based on ordinal scales (e.g. Domin, 1928;Braun-Blanquet, 1951) were converted to percentage cover (van der Maarel, 1979).
For each species, we measured two aspects of the species' abundance across the vegetation plots. First, we calculated its "mean cover": the arithmetic mean of the percentage cover values from all the vegetation plots at 2.5 arc-min raster cells in which the species was present in EVA. Second, we evaluated the frequency distribution of these percentage cover values (see Figure 1 for details on the procedure for three example species). For this, we computed the shape of the distribution function of the percentage cover values. In general, those values are not normally, exponentially or log-normally distributed (Figure 1a-c); thus, we developed a non-parametric approach for measuring the shape of the distribution function. This was achieved by calculating the distribution quantiles in 5% steps, resulting in 20 quantile values. We then fitted a non-linear model on those 20 quantile values and obtained the estimate and the credible F I G U R E 1 Examples of distribution of species' cover values from vegetation plots and calculated mean cover value for (a) Achillea nobilis, the species with the lowest mean cover value, (b) Atriplex portulacoides, a species with intermediate mean cover value and (c) Carex elongata, the species with the highest mean cover value. Note the log scale for frequency. Distribution quantiles from species' cover values were calculated and used to compute the shape of the frequency distribution function for each species, respectively (d-f). Non-linear models on the extracted quantile values were applied to calculate the area under the histograms of cover values (AUH), ranging from 0 to 1, with values close to 0 indicating a strongly right-skewed distribution whereas values close to 1 point to a strongly left-skewed distribution of cover values SPORBERT ET al. interval of the area under the histogram (AUH) (Figure 1d-f). We applied a Bayesian Markov chain Monte Carlo (MCMC) method following Feng et al. (2017), using an exponential distribution, 0.95 confidence level and 10,000 iterations. The resulting AUH value for a given species ranged from 0 to 1, with values lower or higher than 0.5 meaning that the distribution of cover values for a focal species is right-or left-skewed, respectively. The lower the AUH value, the higher was the rarity (i.e. the proportion of relatively low cover values). Thus, the AUH values are suitable as proxies for abundance structure across the vegetation plots. Hereafter, we refer to the AUH values as "skewness of cover values" and use it as an alternative metric, additional to mean cover, to assess across-plot species abundance.

| Explanatory variables: plant functional traits
We compiled a complete species trait matrix with 20 plant functional traits (see Table 2 and Appendix S1). The trait matrix included nine binary variables: five for life form following Raunkiaer (1934); three for life cycle (derived from BiolFlor database; Kühn et al., 2004); and one for clonality (derived from the CLO-PLA database; Klimešová et al., 2017). We included information on 11 continuous trait variables from the global plant-trait database TRY (Kattge et al., 2020 et al., 2015) to fill gaps in the observed species-by-trait matrix data received from TRY. Continuous trait variables were ln-transformed prior to analysis.
A PCA of the 20 traits included in this study was generated using the package factoextra (Kassambara & Mundt, 2017)

| Statistical modelling: linking plant functional traits to mean cover values, skewness of cover values, geographic range size, and climatic niche size
We used the function phylo.maker from the package V. PhyloMaker (Jin & Qian, 2019) to create a phylogenetic tree of the studied species. The function phylo4d from the package phylobase (Hackathon et al., 2013) was applied to link trait data to the species' phylogeny. We applied Pagel's Lambda statistic (Pagel, 1999) (Barton, 2019). We allowed for a maximum of three predictor terms to be included in a given candidate model (m.max = 3); with this, univariate models were applied for single traits and multivariate models for trait combinations (see Appendix S4 for the lists of candidate models for the four response variables).
Finally, the Akaike Information Criterion (AIC), with ∆AIC < 2 was used to identify the most parsimonious candidate model with a maximum of three predictor terms for each of the four studied response variables.
We computed the variance inflation factor (VIF) for each predictor term in the most parsimonious models to check for potential multicollinearity issues among the continuous predictor variables, using the function vif PCs as predictor terms in the models described above and tested for F I G U R E 2 Principal component analysis of the 20 traits included in this study. Colour represents the trait contributions (%) to the PCA (first and second components). The first and second components accounted for 18.8% and 14.6% of the total variation in trait values, respectively. For abbreviation of the trait names see Table 2 SPORBERT ET al. combinations and interactions between PCs in the same way as described for traits. By applying the dredge function we tested all possible combinations of the predictor variables that contributed the most, including two-way interaction between PCs, for each of the four response variables. AIC with ∆AIC < 2 was used to identify the most parsimonious candidate model with a maximum of three predictor terms for each of the four studied response variables (see Appendix S4 for the trait contributions [loadings] to all 20 PCs).

| Local abundance metrics: mean cover and skewness of cover values
Species' mean cover from all the vegetation plots in which a species was present ranged from 2.4% for Achillea nobilis to 24.8% for Carex elongata (Figure 1a and c). The interquartile range of was 4.6% to 8.1% and the median was 5.9%. Species' skewness of cover values ranged from 0.081 (strongly right-skewed distribution of low cover values) in Achillea nobilis to 0.385 in Atriplex portulacoides (Figure 1d and e). The interquartile range was 0.158-0.226 and the median was 0.180. Species' mean cover was positively related to species' skewness of cover values in a phylogenetically corrected model (R 2 = 0.763, p-value < 0.001; Appendix S4).

| The contribution of functional traits to explaining values and skewness of cover values
All response variables were better explained by trait combinations than by single traits, e.g. among the list of candidate models for mean cover as response variable, the best univariate model with plant height as predictor variable was ranked 17 (see Appendix S4).
SLA was the strongest predictor for all response variables; with species with high SLA having larger range sizes, broader climatic niche sizes, and higher local abundances. For each specific response variable, SLA interacted with other traits to give scale-specific and different trait responses. Specifically, geographic range size was larger in species with taller stature and lower leaf N:P ratio. In contrast, climatic niche size was larger in species that had lower leaf C content. The mentioned functional traits were significantly, though not strongly related to species' geographic range size (R 2 = 0.090, p-value < 0.001) and climatic niche size (R 2 = 0.069, p-value < 0.001) in the phylogenetic generalized least-squares models (Table 3, Figure 3a and b). Species' mean cover and the skewness of cover values was higher in species with higher SLA value and with higher leaf area values. The interaction of the variables SLA and leaf area was Journal of Vegetation Science SPORBERT ET al. positive and these functional traits were significantly related to both species mean cover (R 2 = 0.211, p-value < 0.001) and the AUH measure of the skewness of cover values (R 2 = 0.169, p-value < 0.001; Table 3, Figure 4a and b).
The specific traits identified in the final multivariate models for the four response variables had also high loadings on the PCs that were identified as important predictors in the multivariate PCAbased models (see Appendix S4 for the trait contributions [loadings] for all 20 PCs). In addition, the axes captured some more traits with maximum absolute loadings that were not selected in the final trait-based models, such as leaf area for geographic range size, leaf P content for climatic niche size, clonal growth for mean cover and leaf area for the skewness of cover values. However, for each of the four response variables, the three PCs in the final models explained less variation than the trait-based models: geographic range size (R 2 = 0.063, p-value < 0.001; in the sequence of importance, the model included PCs 12, 2 and 1), climatic niche size (R 2 = 0.069, pvalue < 0.001; based on PCs 1, 12 and 19), mean cover (R 2 = 0.123, p-value < 0.001; based on PCs 13, 2 and 4) and skewness of cover values (R 2 = 0.094, p-value < 0.001; based on PCs 6, 2 and 4).

| D ISCUSS I ON
Species' local abundances (i.e. a measure of commonness) were more strongly related to traits than were species' broad-scale distribution patterns in the geographic and climatic space. This indicates that plant traits better capture local processes acting at the community level (such as biotic processes) than broad-scale macroecological processes. Both local abundances and broad-scale distribution patterns were better predicted by combinations of traits than by single traits.
Specific leaf area had a significant positive effect and explained most of the observed variation in all four models predicting species' local abundance and broad-scale distribution patterns. SLA is a productivity-and competitive ability-related trait, that reflects species strategies for rapid acquisition of resources, with higher SLA values allowing a species to capture more light for a given biomass investment in leaves, for example (Díaz et al., 2004;Wright et al., 2004;Mariotte, 2014). In line with our findings, several studies state common species to be associated with higher SLA (Grime et al., 1997;Díaz et al., 2004;Mariotte, 2014;Lachaise et al., 2020). While species' local abundances were best predicted by the interaction between leaf area and SLA, reflecting the leaf economics spectrum trait syndrome (Díaz et al., 2004), broad-scale distribution metrics were best predicted by different combinations of traits. While geographic range size increased with increasing plant height, climatic niche size decreased with increasing leaf carbon content, and both increased with increasing SLA.
At the local scale, leaf area showed a significantly positive effect on species abundance. This result offers a functional explanation that species with larger leaves, allowing better light capture, are able to attain higher local abundances than species with smaller leaves (Mariotte, 2014). Moreover, leaf area was particularly important in F I G U R E 3 Scatter plot of observed values and regression lines from phylogenetic generalized least-squares models, showing the effects of the three most predictive terms on species' (a) geographic range size and (b) climatic niche size. For geographic range size, coloured and dashed lines represent the 5th, 50th and 95th percentile in values for plant height and leaf N:P ratio, respectively; for climatic niche size, coloured lines, solid and dotted, represent the 5th, 50th and 95th percentile in leaf C content values for therophytic and non-therophytic species, respectively SPORBERT ET al.
interaction with SLA values, as species' local abundance was higher in species with large leaves and high SLA values. Like SLA, leaf area is interpreted as a trait that is positively related to species' productivity and competitive ability (Diaz et al., 2004;Wright et al., 2004). This indicates that being nearer the "fast" end of the leaf economics spectrum (i.e. trait syndrome) tends to increase local abundance. In contrast, the ability to grow clonally was not selected in any of the final models based on original trait variables. The ability of species for clonal growth plays an important role both in short-distance spread and in persistence within habitats (Benot et al., 2013) and previous studies found clonality to be positively associated with local abundance (Eriksson & Jakobsson, 1998;Kolb et al., 2006). Accordingly, this trait had high axis loadings in the PC-based model with mean cover value as response variable. The absence of clonal growth from our final models based on original trait values is probably due to its negative correlation with SLA in our set of species.
At large spatial extent, geographic range size of species was positively related to plant height (i.e. with taller species being more widespread). High stature is known to have a competitive advantage and to be associated with common species. Greater plant height of widely distributed species suggests that these species may have higher competitive ability for space and light than narrowly distributed species. A similar positive correlation between plant height and geographic range size was found for herbaceous species in the French Mediterranean region (Lavergne et al., 2004) and in temperate forests in Germany (Kolb et al., 2006). We found leaf N:P ratio to be negatively correlated with geographic range size. Nitrogen (N) and phosphorus (P) availability can limit plant growth in terrestrial ecosystems, and N:P ratios are on average higher in stress-tolerant species compared to ruderals (Güsewell, 2004). Ruderal species are characterized by rapid growth and they establish much quicker and thrive better in disturbed habitats than stress-tolerant and competitor species (Grime, 1979;Wright et al., 2004;Guo et al., 2018) and generally undergo long-distance dispersal (Baker, 1965). Thus, a plausible explanation is that ruderality has a positive effect on species' geographic range size. As shown in a global study by Bruelheide et al. (2018), species' leaf N:P ratio declines at higher latitudes. Many species primarily found in boreal regions obtain broader geographic range sizes in comparison to species mainly found further south (e.g. in Mediterranean regions), presumably because of post-glacial re-expansion. This might be another plausible explanation for the negative relationship between geographic range size and leaf N:P ratio in our study.
In our study, the distribution range in climatic space was larger in species with lower leaf carbon content, even when accounting for SLA.
In general, carbon content is expected to be negatively related to SLA (Reich, 2014), but both traits seem to explain independent variation in climatic niche size. This was brought about by species with broad climatic niche sizes, for which SLA alone was a poor predictor, such as species with a tendency to succulence (e.g. Plantago major), which have leaves with low SLA but yet low leaf carbon content. This indicates that species following a "fast" strategy, according to the leaf economics spectrum, are better adapted to obtain broad climatic niche sizes. species that generally undergo long-distance dispersal (Baker, 1965).
Finally, traits related to dispersal, regeneration and persistence were not significantly correlated with local abundance or broad-scale distribution in our models. These results confirm previous studies that found no significant relationship between species' abundance or broad-scale distribution and seed mass (Thompson et al., 1999;Leishman & Murray, 2001;Lavergne et al., 2004) or life cycle length (Kolb et al., 2006).
Our results largely confirm trends previously reported about the existing association between species' geographic range and climatic niche size, with widely distributed species also having broad climatic tolerances and geographically narrowly distributed species also narrowly distributed in climatic space (Gaston, 2000;Slatyer et al., 2013;Cardillo et al., 2019). We found an overall right-skewed distribution in cover values for most of the studied species, with species exhibiting low cover at most sites and high cover in only a few sites across their distribution range. The species mean cover values were positively related to the skewness of cover values. Therefore, for our species set, we consider the measure of skewness, calculated as the AUH, a robust tool to capture both the mean and the variability of cover values across a species' whole distribution range.
Nevertheless, single plant traits and trait syndromes only weakly explained the total observed variation in species' broad-scale distribution metrics. We see two plausible explanations for the weak predictive power of functional traits on species' broad-scale distribution metrics.
First, our study was carried out on the species' whole Eurasian distribution range, which includes a wide range of habitat types and bioclimatic zones. Species functional traits are expected to be related to those environmental conditions under which the species occurs (Lavorel & Garnier, 2002;McGill et al., 2006). With this, both widely distributed and geographically restricted species might be characterized by the same traits in different habitats, vegetation types or geographic regions, depending on the local conditions (Aerts & Chapin, 2000). Therefore, future studies should incorporate habitat variability, by means of comparisons among single habitat types (e.g. by applying EUNIS habitat classification; Chytrý et al., 2020), and test for consistency of the role of traits for patterns of species commonness at different spatial scales.
Second, in this study, we used mean trait values derived from trait databases. Several studies have provided evidence that functional traits express not only species-specific characteristics, but also intraspecific variability in leaf traits (Reich & Oleksyn, 2004;Albert et al., 2011;Moles et al., 2014;Niinemets, 2015;Wright et al., 2017). This intraspecific trait variation may influence the interactions among and between species and their environment and, therefore, might influence species performance (Bolnick et al., 2003;Siefert et al., 2015). Therefore, we encourage future studies to include intraspecific trait variation in addition to mean values for species traits when investigating studies over large geographic scales.

ACK N OWLED G EM ENTS
We thank all the scientists who collected vegetation-plot data and traits in the field and/or converted them to electronic databases, the custodians of the databases represented in EVA and TRY, the EVA database managers Stephan Hennekens, Borja Jiménez-Alfaro and Ilona Knollová, the TRY database managers Jens Kattge and Gerhard Boenisch, and sPlot database manager Francesco Maria Sabatini whose contributions were essential for this broad-scale study.

AUTH O R CO NTR I B UTI O N S
MS, EW and HB conceived the study. GS and MS harmonized data retrieved from EVA and CDH. MS harmonized data retrieved from TRY with data retrieved from EVA and CDH. EW, GS and MS developed the measures for niche properties and abundance skewness.
MS and HB carried out statistical analyses. MS produced the graphs.
MS and HB wrote the paper, MS led the writing. All other authors contributed data, discussed the results and commented on, and/or substantially edited, the manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data used in this paper are from large multi-contributor databases (EVA, TRY). They cannot be made publicly available because of the third-party ownership issues. The data selections released for this study are stored in internal repositories of the source databases and can be made available for re-analyses upon request. The EVA dataset is stored in the EVA repository with reference to project no. 24.