Global changes explain the long-term demographic trend of the Eurasian common lizard (Squamata: Lacertidae).

The demographic trend of a species depends on the dynamics of its local populations, which can be compromised by local or by global phenomena. However, the relevance of local and global phenomena has rarely been investigated simultaneously. Here, we tested whether local phenomena compromised a species' demographic trend using the Eurasian common lizard Zootoca vivipara, the terrestrial reptile exhibiting the widest geographic distribution, as a model species. We analyzed the species' ancient demographic trend using genetic data from its 6 allopatric genetic clades and tested whether its demographic trend mainly depended on single clades or on global phenomena. Zootoca vivipara's effective population size increased since 2.3 million years ago and started to increase steeply and continuously from 0.531 million years ago. Population growth rate exhibited 2 maxima, both occurring during global climatic changes and important vegetation changes on the northern hemisphere. Effective population size and growth rate were negatively correlated with global surface temperatures, in line with global parameters driving long-term demographic trends. Zootoca vivipara's ancient demography was neither driven by a single clade, nor by the 2 clades that colonized huge geographic areas after the last glaciation. The low importance of local phenomena, suggests that the experimentally demonstrated high sensitivity of this species to short-term ecological changes is a response in order to cope with short-term and local changes. This suggests that what affected its long-term demographic trend the most, were not these local changes/responses, but rather the important and prolonged global climatic changes and important vegetation changes on the northern hemisphere, including the opening up of the forest by humans.

Ecological parameters, geographic barriers, and changes therein affect demographic processes, species distributions, genetic diversity, and speciation (Parmesan and Yohe 2003; Geffen et al. 2004; Hewitt 2004; Begon et al. 2006; Sharma et al. 2013; Horreo et al. 2016). Because some of these drivers are global and others are local, global and local processes may be important drivers of a species’ demographic trend. Global changes (e.g., Wallis et al. 2016), for example, the repeated cycles of climatic warming and cooling during the Late Quaternary (Hansen et al. 2013), affect large geographic areas and they may lead to similarities in the demographic trend of species with similar ecological requirements and to similarities among local populations belonging to the same species (Weiss and Ferrand 2007; Horreo and Fitze 2018). In contrast, changes at smaller geographic extents may lead to differences among allopatric genetic units and local populations belonging to the same species. But they may not necessarily affect a species’ demographic trend, because local phenomena (e.g., changes in regional or local climatic conditions, genetic drift and differences in sexual selection; Uyeda et al. 2009; Nosil 2012) may only affect a few local populations (e.g., Weiss and Ferrand 2007; Horreo et al. 2014). Moreover, longer-term changes and changes of big magnitude may affect a species’ genetic diversity, speciation (Wallis et al. 2016), and extinction rate (Sinervo et al. 2010), whereas short-term and local changes that rapidly feedback on the dynamics of local populations (Geffen et al. 2004; Romero-Diaz et al. 2017; Horreo et al. 2019; Masó et al. 2019) may not necessarily affect a species’ genetic diversity and demographic trend. Many studies investigated climatic effects on trends of local populations (e.g., Letnic et al. 2004; Chamaillé-Jammes et al. 2006), but only few studies simultaneously investigated the relevance of global and local events for a species’ demographic trend and they only analyzed short and recent time-intervals (Grimm-Seyfarth et al. 2018). Consequently, the relevance of considering a species’ entire geographic distribution to derive its long-term trends is unknown.

Here, we used the Eurasian common lizard Zootoca vivipara (Lichtenstein 1823), the most widely and farthest north distributed terrestrial reptile (Takeuchi et al. 2013), as a model species. We investigate whether a species’ demographic trend is driven by climatic changes across the species’ entire range, or whether local genetic clades vary independently according to local conditions, which will reveal whether local or global phenomena drive the range-wide dynamics of a broadly distributed species. Zootoca vivipara exhibits 6 well-differentiated allopatric genetic clades (Surget-Groba et al. 2001; Horreo et al. 2018a), it inhabits several biogeographic regions (Horreo et al. 2018b), and its distribution ranges from Ireland and Western Spain in the West to Japan and Sakhalin in the East, and from Northern Spain and the Southern Balkans in the South to Northern Scandinavia and Siberia in the North. Zootoca vivipara originated approximately 4.4 million years ago and all 6 clades differentiated more than 2 million years ago (Mya; Horreo et al. 2018b). The clades exhibited important differences in range expansion/contraction and occupied different glacial refugia during the Pleistocene (location of the glacial refugia: Clade A, South of the Alps; Clade B, Northern Spain and Southwestern France; Clades C and F, Pannonian/Vienna Basin; Clade D, north of the black sea; Clade E, north of Pannonian/Vienna Basin; Horreo et al. 2018b). Zootoca vivipara is a strongly hydrophile species that lives in bogs and humid meadows (Grenot et al. 1987; Lorenzon et al. 1999), which is also reflected in its German (Mooreidechse), Spanish (lagartija de turbera), and French (lézard des tourbières) common species name, the literal translation of which means “Boag lizard.” Zootoca vivipara is highly sensitive to small differences in environmental parameters including climatic parameters (Grenot et al. 1987; Le Galliard et al. 2005; Bestion et al. 2015; Romero-Diaz et al. 2017; Masó et al. 2019). In this and many other ectothermic species, a bell-shaped relationship exists between ambient temperature and morphological, physiological, and behavioral parameters (Avery 1985). The performance of most ectotherms improves smoothly with increasing body temperature until reaching an optima, and thereafter it drastically decreases (Huey and Stevenson 1979). Zootoca vivipara’s preferred body temperature in May/June, that is, during reproduction, is between 30 °C and 34 °C in males and around 29.5 °C in gravid females, and the range over which its performance is at least 80% of maximal performance is quite narrow ±2–5 °C (Van Damme et al. 1991; Gvoždík and Castilla 2001; Rozen-Rechels et al. 2020). Its critical thermal maxima is between 43.8 °C and 44.1 °C. At this temperature its righting reflex is lost, leg spasms happen, and lethal or near lethal injury occurs (Gvoždík and Castilla 2001; Leon and Bouchama 2015). Operative temperatures, that is, the body temperature of a perfectly thermoconforming animal (Dzialowski 2005), reach more than 50 °C in the early afternoon at the beginning of May and on most study days, they were much above Z. vivipara’s preferred body temperature (measured in May 2018, Rozen-Rechels et al. 2020). However, in May, the daily maximal temperatures are much below the daily maxima occurring in summer (e.g., in May they were −7.2 °C lower than in summer 2018 in Nemours les Saint-Pierre, France, where Rozen-Rechels et al. [2020] recorded operative temperatures [click on this weblink to see the daily maximal temperatures at Nemours les Saint-Pierre]). As a consequence, Z. vivipara exhibits a unimodal activity pattern with an activity peak around midday when daily maximal temperatures are below 34 °C and evapotranspiration is not too high. Especially in southern populations, in summer, Z. vivipara exhibits a bimodal activity pattern with activity peaks in the morning and in the afternoon, and the hot period they spend below ground (or in the shade, Rozen-Rechels et al. 2020) where conditions are cooler, wetter, and water loss rates are lower (Rozen-Rechels et al. 2020).

Consequently, if ambient temperature is higher than optimal, lizards reduce their activity time (Rozen-Rechels et al. 2020) and gravid females will first respond given their lower thermal preferences (Gvoždík and Castilla 2001). Global warming as well alters surface water availability and reduced surface moisture enhances surface heating even more (Rozen-Rechels et al. 2020; Zhou et al. 2021). Both the increasing temperatures and the water restriction during summer reduce the length of Z. vivipara’s daily activity window (Rozen-Rechels et al. 2020), which may compromise its reproduction, and thereby lead to population decrease, and local extinction (Kearney 2013; Bestion et al. 2015; Sinervo et al. 2010). In contrast, in the presence of lower than optimal ambient temperatures ectotherms can achieve optimal body temperatures through thermoregulation and populations remain stable or increase. Moreover, Z. vivipara is a freeze tolerant species (Berman et al. 2016) that inhabits the northernmost latitudes of all terrestrial reptiles and thus habitats where competition with other reptiles is low.

The large geographic distribution of Z. vivipara, the allopatry of its genetic clades, the different evolutionary histories, and its sensitivity to climatic parameters, suggest that changes in a single clade (i.e., local changes) may be determinants of its population size and population growth rate. To test whether the trajectory of a single clade affected the demographic trend of Z. vivipara, we first derived Z. vivipara’s demographic history (i.e., ancient effective population sizes and ancient population growth rates) from DNA sequences using the Bayesian skyline plot (BSP) method (Drummond et al. 2005). Thereafter, we tested whether the demographic patterns observed in the entire species (i.e., including all genetic clades, hereafter referred to as “global model”) were importantly affected by a specific clade. To this end, we compared the results of the global model with those of models excluding 1 clade at a time. We also tested whether clades that colonized huge areas of Northern Europe, Eastern Europe, and Asia after the last glaciation, affected the species’ effective population size and population growth. Additionally, we tested whether changes in effective population size and population growth are explained by global ecological changes, by correlating effective population size and population growth rate with the prevailing global surface temperatures (Hansen et al. 2013).

Materials and Methods

A recently published dataset including DNA sequences of fragments of 3 nuclear (ZV1–ZV3; Horreo et al. 2018a, 2018b) and 3 mitochondrial (16 s, cytB, and ND2) genes, and consisting of 2,600 base pairs from 231 Z. vivipara (Lichtenstein, 1823) specimens (Horreo et al. 2018b) was used for this study. The 231 sequences belong to the 6 described Z. vivipara clades (A–F) and stem from 185 localities covering the majority of Z. vivipara’s natural distribution (from Spain and Ireland in the West to Sakhalin [Russia] and Japan in the East and from Northern Greece in the South to Sweden and Finland in the North; Figure 1 and Online Appendix S1). Clades A and B correspond to the oviparous, and Clades C–F to viviparous clades (Horreo et al. 2018b). The number of sampled individuals per clade (Table 1) ranged from 24 (Clades C and F) to 63 individuals (Clade B).

Figure 1.

Maps showing the sampled populations and the geographic distribution of Z. vivipara (Lichtenstein, 1823) across Eurasia (black area of the upper left map). Letters and colors correspond to the 6 clades (A–F) and areas highlighted in color to the geographic distribution of the clades. In the species distribution, the square highlighted in blue corresponds to the location of the main map and the 3 orange dots to the sampled populations of Clade D, which are outside the range of the main map.

Table 1.

Genetic variability of the global dataset (“all” specimens) and of each clade (A–F)

Analyzed specimens	N	S	h	Hd	Nd
All	231	226	112	0.977	0.012
Clade A	38	82	25	0.937	0.004
Clade B	63	73	31	0.902	0.007
Clade C	24	17	16	0.949	0.002
Clade D	30	69	19	0.903	0.005
Clade E	52	63	37	0.960	0.005
Clade F	24	28	17	0.956	0.005

N, number of specimens per analyzed dataset; S, number of polymorphic sites; h, number of haplotypes; Hd, haplotype diversity; Nd, nucleotide diversity.

The genetic variability of the total dataset as well as per Z. vivipara clade was estimated and the number of polymorphic sites and haplotypes determined using DNAsp (Librado and Rozas 2009). Haplotype diversity (Hd) and nucleotide diversity (Nd) were calculated with Arlequin version 3.5 (Excoffier and Lischer 2010).

Ancient effective population sizes (Ne) were inferred using BEAST version 2.3.1 (Bouckaert et al. 2014), the coalescent BSP approach (Drummond et al. 2005), and the extended BSP (eBSP; Heled and Drummond, 2008). These analyses allow to infer ancient effective population sizes independent of a pre-specified parametric model of demographic history (Horreo et al. 2016). BSPs and eBSP were calibrated using molecular dating from a previous study (Horreo et al. 2018b) that dated the ancestor of Z. vivipara (i.e., the split between Clade A and the other clades) to 4.4 (4.2–4.6, 95% confidence interval) Mya using the 6 gene fragments used for this study. A strict clock and a generalized time reversible model of sequence evolution (the most general neutral and independent model) were used because JModelTest version 0.1.1 (Posada 2008) unraveled it as being the best-fit model of nucleotide substitution. Markov chain Monte Carlo methods implemented in BEAST were used and the number of chains determined based on the effective sample size (ESS). A 10% burn-in and 15 million chains were used in order to obtain ESSs > 200, which ensured reasonable confidence of the analyses. First, a global model was fit that included all 231 specimens and all 6 Z. vivipara clades (Horreo et al. 2018b) using BSP and eBSP. Both types of analyses rendered the same dynamics and estimates of effective population size were highly correlated (Spearman rank correlation: rs = 0.9999). Consequently, only BSPs were used for the subsequent analyses. To investigate whether global dynamics may represent the dynamics of a single clade, 1 clade was excluded from the global model and the BSP was re-fit. A difference in the chronology of Ne between this model and the global model would show that the species’ dynamics was importantly affected by the excluded clade. Given that 6 Z. vivipara clades exist (A, B, C, D, E, or F) 6 different BSPs were fit, each of them excluding another clade from the global model. In addition, previous phylogeographic analyses showed that Clades D and E colonized huge areas of Northern Europe, Eastern Europe, and Asia after the last glaciation (Surget-Groba et al. 2001; Horreo et al. 2018b) and their joint distribution spans the majority of Z. vivipara’s geographic distribution. Therefore, an additional BSP including Clades A, B, C, and F (i.e., all clades that did not colonize huge geographic areas after the last glaciation) was fit, to test whether the species’ dynamics may represent the dynamics of Clades D and E, or rather global dynamics.

To test whether the effective population sizes (Ne ; calculated from present to 4.4 Mya) and population growth rates (Δ effective population size/Δ time; with Δ effective population size = Ne−Ne+1 and Δ time = t−t+1; i corresponds to the ith estimate of Ne from present and t corresponds to the time of this estimate) were correlated to the global climate, global surface temperatures (Ts; in degree Celsius) of the last 4.4 million years (i.e., because the evolution of Z. vivipara; Horreo et al. 2018b) were obtained from Hansen et al. (2013, i.e., from its Figure 4B). We used global surface temperatures because the temporal resolution of the reconstructed values and the time horizon (present to 60 Mya) were large enough to allow for correlations with Ne of Z. vivipara. Global surface temperatures are of high interest to humanity (Hansen et al. 2013) and highly relevant for ectothermic species and encountered reconstructed global precipitation datasets did not meet the resolution or time horizon required to correlate them with Ne of Z. vivipara. For each date (t) rendered by BSP, the corresponding Ts was attributed (Figure 2). If the dates (measured in million years and rounded to 4 decimal places) rendered by BSP (tBSP,) did not coincide with a precision of 4 decimal places with the dates of Hansen’s climatic data (tHansen,), a rule of proportions (i.e., rule of 3; Cajori 1993, p. 164) was applied to reconstruct Ts at tBSP,. For the correlation between population growth rate and surface temperature, average Ts () were derived for each time interval (Δtime) over which population growth rates were calculated. Previous to the statistical analyses, for Figure 4A and B, effective population size and both variables were log transformed, respectively. Generalized additive models (GAMs) were used to test for significant correlations using the “gam” function of the “mgcv” package in R (Wood and Augustin 2002; R Core Team 2018).

Figure 4.

Correlations between (A) global surface temperatures (Ts) and effective population size of Z. vivipara and (B) average global surface temperatures () and population growth rate measured between 4.336 Mya and present. Shown are data and GAM model predictions (blue lines).

Figure 2.

BSP showing the temporal variation of Z. vivipara’s effective population size (Ne, pink), 95% confidence intervals (light pink), and growth rate (ΔNe/Δtime, blue) for the entire species (global dataset) during the last 4.4 Mya. The y-axis is given on a logarithmic scale.

Results

The genetic variability in the global dataset and in each clade is shown in Table 1. The number of polymorphic sites per clade ranged between 17 (Clade C) and 82 (Clade A), and the number of haplotypes between 16 (Clade C) and 37 (Clade E). Hd was high and it ranged between 0.902 (Clade B) and 0.960 (Clade E), and Nd was between 0.002 (Clade C) and 0.007 (Clade B).

The BSP shows that Z. vivipara’s effective population size started to increase from 2.3 Mya and it started to steeply increase until now from approximately 0.7 Mya onward (pink line, Figure 2). Population growth rate exhibited 2 pronounced maxima, between 0.531 and 0.487 Mya (this flatter peak also included 0.487–0.443 Mya, when growth rate was only slightly smaller) and between 0.133 and 0.088 Mya (blue line; Figure 2). The chronology of the species’ past demography was very similar between the analysis based on all clades and the analyses excluding 1 clade (Figure 2 versus Figure 3) or the 2 clades with current huge geographic distribution (Clades D and E; Supplementary Figure S1).

Figure 3.

BSP showing effective population size (Ne, solid black line) and 95% confidence intervals (blue line) during the last 4.4 Mya based on 5 of the 6 Z. vivipara clades (in each analysis another clade was excluded).

Effective population size and global surface temperatures were significantly and negatively correlated (Ts, ; Figure 4A) and population growth rate was significantly and negatively correlated with average global surface temperatures measured over the same time interval as population growth (, ; Figure 4B).

Discussion

Zootoca vivipara’s effective population size (Ne) started to increase from approximately 2.5 Mya onward (Figure 2), after the trend of average global surface temperatures dropped below 14 °C (see Figure 4B in Hansen et al. 2013). The increase started to be steep from 0.531 Mya onward until present (Figure 2), that is, shortly after the trend of average global surface temperatures dropped below 12 °C (see Figure 4B in Hansen et al. 2013). There was a significant negative association between global surface temperatures and effective population size (Figure 4A) and population growth rate of Z. vivipara (Figure 4B), which is in line with the important cold tolerance of Z. vivipara (Costanzo et al. 1995) and the fact that it is the reptile with the northernmost distribution (Berman et al. 2016).

Zootoca vivipara’s population growth rate exhibited 2 pronounced maxima (Figure 2). The first growth rate peak (between 0.531 and 0.443 Mya) appeared just before the mid-Brunhes event (MBE; 0.43 Mya; Tzedakis et al. 2009; Meckler et al. 2012), which was a climatic transition between 2 climatic modes: from low sea levels, large ice sheets, and moderate warmth from early to mid-Pleistocene interglacials (0.78–0.45 Mya; Helmke et al. 2003; Meckler et al. 2012) to higher sea levels, greater warmth, and smaller ice-sheets from mid- to late-Pleistocene interglacials (occurring after 0.450 Mya; Tzedakis et al. 2009; Blain et al. 2012). During the MBE, the intertropical convergence zone displaced northward and intensified and expanded precipitation in widespread northern hemisphere regions (compared with previous interglacials; Ao et al. 2020). The MBE was the most extreme climatic shift in the last 0.80 Mya, with average interglacial temperatures reaching levels comparable with, or higher than present temperatures (Jansen et al. 1986; Candy et al. 2010). The MBE was a global event with a significant impact on many components of the climate system (Barth et al. 2018), including an increase of mean global temperatures during Northern Hemisphere winters (Yin and Berger 2010). However, higher temperatures are generally associated with lower population growth of Z. vivipara (Figure 4), suggesting that the increased humidity may have buffered their negative effects (Rozen-Rechels et al. 2020). Additionally, warming may have allowed for Z. vivipara’s northward and altitudinal expansion. Moreover, habitats with closed tree canopies are unsuitable for Z. vivipara, and forest clearings, forest verges (Dent and Spellerberg 1987; Glandt 2006) are preferred habitats (Grenot et al. 1987; Pilorge 1987), suggesting that Z. vivipara may also have benefited from early human activities, because novel hominin behavioral developments including fire control, hand axe industries, and Levallois techniques (appearing from ∼500 ka onward) may have opened up the forest (Ao et al. 2020; Moncel et al. 2020) and increased the amount of suitable habitat.

Population growth rate of Z. vivipara exhibited a second peak between 0.088 and 0.133 Mya, a period that mainly includes the time following the maximal global surface temperature Ts of the Eemian interglacial (i.e., after 0.119 Mya; Hansen et al. 2013), and thus a period when temperatures started to decrease toward the last ice age (Hansen et al. 2013). During the Eemian, dramatic habitat changes occurred (Rioual et al. 2001; Guiter et al. 2003). At this time, European forests opened up and a physiognomically steppe-like vegetation (Bell 1969) predominated by herbaceous arctic tundra or steppe vegetation (Woillard 1978; de Beaulieu and Reille 1992) appeared in large parts of Europe, the so-called steppe-tundra (Tzedakis and Bennett 1995; Van Andel and Tzedakis 1996). The steppe-tundra was rich in open vegetation species, for example. Poaceae, Artemisia and accompanying Pinus and Betula (Klotz et al. 2004; Helmens 2014) and the appearance of this new habitat has been related to the expansion of several species, especially of large herbivorous mammals such as mammoth, giant deer, horse, and bison (Zimov et al. 1995; Řičánková et al. 2014). At the end of the Pleistocene, there was a shift from steppe-tundras toward more humid moss tundras (Zimov et al. 1995), from which Z. vivipara may have additionally benefited, because of its preferences of hydric preferences (Grenot et al. 1987; Pilorge 1987) and its high densities in bogs, humid meadows, and tree less areas.

The significant and negative associations between global surface temperatures and effective population size (Figure 4A) or population growth rate of Z. vivipara (Figure 4B) leveled off at surface temperatures >14 °C, when the population growth was close to 1 (i.e., when populations remain stable and do neither grow nor decline; Figure 4B). Likewise, steep population increase peaked when average global surface temperature ( below 14 °C, namely at 12.281 °C ± 0.503 sd during peak growth previous to the MBE, and at 12.990 °C ± 0.973 sd during peak growth around the Eemia. Higher population growth at cooler temperatures (below global air temperatures of 14 °C) is in line with the important cold tolerance of Z. vivipara (Costanzo et al. 1995). The average global surface temperatures prevailing during the 2 peaks of population growth were quite similar and they can hardly explain the huge differences in the peaks’ magnitude (Figure 2). Moreover, during the other interglacials, no growth rate peaks existed (Figure 2), showing that global temperature change alone cannot explain the steep growth rate increase. Both peaks are best explained by the opening of the forests and increased humidity, and the difference in the height of the growth peaks most likely depends on the magnitude of their effect. Moreover, the population growth continued between the MBE and the Eemian, and between the Eemian and the present (Figure 2), when temperatures and most likely as well the precipitation varied. This suggests that the opening of the northern hemisphere forests induced by humans and large animals from the start of the Eemian until present (Ao et al. 2020, p. 2) may have been the most important factors for its demographic trend.

Both the MBE and the Eemian were global rather than regional phenomena (Van Kolfschoten et al. 2003; Barth et al. 2018) during which important climatic and vegetation changes happened. Their effect on Z. vivipara’s ancient demographic trends is in line with the finding that Z. vivipara’s population growth rate did not result from an increase of a single clade, but rather from multiple clades (Figures 2 and 3, and Supplementary Figure S1). Similarly, the colonization of huge areas by Clades D and E (i.e., colonization of the areas which were covered by ice during the last glaciation; Last Glacial Maximum: 21,000 years ago; Horreo et al. 2018b) also did not affect Z. vivipara’s effective population size importantly.

Zootoca vivipara exhibits a very large geographic distribution, its genetic clades are allopatric, and the abundant literature shows that it is highly sensitive to small differences in environmental parameters (Grenot et al. 1987; Le Galliard et al. 2005; Romero-Diaz et al. 2017; Masó et al. 2019). Consequently, the finding that multiple clades have been involved in the steep effective population size increases in the past suggests that only important and long-lasting global changes rather than local changes affected Z. vivipara’s long-term demographic trend (Figures 2–4). This suggests that effective population size and population growth rates of Z. vivipara were not importantly affected by warm/cold periods between the MBE and the Eemian or the last ice age, probably due to the species’ high plasticity and the fact that it lives in very humid habitats. What affected population size most were the long-term and prolonged global climatic and the vegetation changes on the northern hemisphere. This suggests that for species which are not close to extinction (i.e., vulnerable to go extinct on the short run), predictive models dealing with short-term effects or a small subsample of a species’ geographic range may not be able to unravel longer-term dynamics.

In summary, our analyses revealed that global phenomena were most important for the species’ demographic trend and it suggests that Z. vivipara’s reactions to short-term ecological changes and local phenomena enabled Z. vivipara to deal with them. Consequently, our results indicate that the current climatic change and the associated vegetation changes will most likely affect the species’ demography much more than local phenomena. These findings are important for safeguarding many species and they suggest that the most effective means for species conservation might be those slowing the currently observed climatic change and the associated vegetation change, which is in line with proposed conservation actions for a wide range of organisms (insects: Lancaster 2016; corals: Osborne et al., 2017; plants: Di Marco et al., 2019).

Funding

J.L.H. was supported by a Spanish MINECO grant IJCI-2015-23618. Project funds were provided by the Spanish Ministry of Education and Science [CGL2012-32459, CGL2016-76918 AEI/FEDER, UE to P.S.F.].

Author Contributions

Both authors designed the study, generated and analyzed the data, and wrote the manuscript. Four anonymous reviewers provided helpful comments.

Conflicts of interest

No conflicts of interest exist.

Supplementary Material

Supplementary material can be found at https://academic.oup.com/cz.

Click here for additional data file.